Abstract

Multisubsystem co-design refers to the simultaneous optimization of physical plant and controller of a system decomposed into multiple interconnected subsystems. In this paper, two decentralized (multilevel and bilevel) approaches are formulated to solve multisubsystem co-design problems, which are based on the direct collocation and decomposition-based optimization methods. In the multilevel approach, the problem is decomposed into two bilevel optimization problems, one for the physical plant and the other for the control part. In the bilevel approach, the problem is decomposed into subsystem optimization subproblems, with each subproblem having the optimization model for physical plant and control parts together. In both cases, the entire time horizon is discretized to convert the continuous optimal control problem into a finite-dimensional nonlinear program. The optimality condition decomposition method is employed to solve the converted problem in a decentralized manner. Using the proposed approaches, it is possible to obtain an optimal solution for more generalized multisubsystem co-design problems than was previously possible, including those with nonlinear dynamic constraints. The proposed approaches are applied to a numerical and engineering example. For both examples, the solutions obtained by the decentralized approaches are compared with a centralized (all-at-once) approach. Finally, a scalable version of the engineering example is solved to demonstrate that using a simulated parallelization with and without communication delays, the computational time of the proposed decentralized approaches can outperform a centralized approach as the size of the problem increases.

1 Introduction

For many engineering systems, the solution to simultaneous optimization of the physical plant and controller design can outperform the solution to sequential optimization of those systems [1]. This simultaneous optimization problem is referred to as co-design. Existing literature in co-design is mainly focused on single-system problems, i.e., physical plant and control are parts of a single system, including analysis and solution approaches of co-design problems of single systems (e.g., Refs. [13]), and a wide range of mechanical (e.g., Refs. [4,5]) and control (e.g., Refs. [6,7]) system applications.

In a recent paper [8], the authors of this paper presented a class of multisystem co-design problems with a finite time horizon linear quadratic regulator. The class of problems in Ref. [8] consisted of multiple subsystems coupled through dynamic equations as well as plant design constraints. A decomposition approach based on an indirect method was proposed in Ref. [8], in which control variables were optimized using the equations derived from the necessary optimality conditions. The proposed techniques can be regarded as an extension of the decomposition-based approach of a single-system co-design problem, as in Ref. [3]. In this paper, two decentralized approaches (i.e., multilevel and bilevel) are developed for general multisubsystem co-design problems based on the direct collocation method and optimality condition decomposition, which differs from the analytical target cascading framework [9] for co-design problems. Using a direct method in the proposed approaches, a continuous-time optimal control problem can be transformed into a finite-dimensional optimization problem by the discretization of the time horizon.

Direct methods have been used to numerically solve optimal control problems (e.g., Refs. [10,11]). Some of the most popular direct methods include the direct shooting method [12,13], pseudospectral optimal control [1416], and direct collocation method [17]. In co-design, the direct transcription method has been applied to single-system problems [4,1820]. In this paper, the direct collocation method is selected for two main reasons. First, from the numerical solution of the state and control variables, a reliable estimate of co-state variables of the control problem can be obtained [17]. Second, when the complexity of the co-design problem increases, the direct collocation method is expected to perform better than other numerical methods [11]. It should be noted that given the proposed decentralized approaches in this paper, the direct collocation method can be replaced by any direct transcription scheme (e.g., direct transcription method in Ref. [4]) to solve the co-design problem. The direct collocation method has been employed in many optimal control applications. For example, Geiger et al. [21] investigated the application of the direct collocation method to optimal path planning of unmanned aerial vehicles (UAVs).

The main contribution of this paper is twofold. First, following the decentralized framework in Ref. [8], a multilevel decentralized approach is proposed to solve multisubsystem co-design problems. It is possible to use decentralized optimal control techniques to solve certain classes of co-design problems. For example, the control and formation optimization of multiple UAVs can be regarded as a co-design problem [22]. Some other methods (e.g., Refs. [23,24]) can be applied to co-design problems with modifications. However, to the best of our knowledge, the papers in the control field do not focus on the physical design of multiple subsystems or the entire system. In this paper, the proposed problem is more general wherein both the physical plant design and control are optimized using decomposition-based methods. In the proposed approach, the direct collocation method is extended and integrated with the optimality condition decomposition technique to solve the control part. Compared with Ref. [8], the new approach has relaxed the model restrictions on the control problem, which makes it possible to obtain an optimal solution for more general co-design problems, including those with nonconvex objective and state constraint functions, as well as nonlinear dynamics in the control part. Second, a bilevel framework is developed for solving multisubsystem co-design problems. For the two cases that the shared plant design variables exist or do not exist, the proposed approaches formulate and solve co-design subproblems for each subsystem. In particular, when there are shared plant design variables, both complicating variables and complicating constraints (as defined in Ref. [25]) exist in the converted optimization problem by the direct collocation method. To the best of our knowledge, no previous literature has addressed a decomposition-based technique for such problems. It is shown that the bilevel approach has a better performance than the multilevel one in terms of both solution quality and computational time. These two decentralized approaches are compared against a centralized (all-at-once) approach. For the tested examples, the bilevel approach obtains exactly the same optimal solution as the centralized one, while the multilevel one obtains a close but slightly worse solution. It is also shown that the decentralized approaches can outperform the centralized approach in terms of computational time as the size of the entire problem increases.

The rest of the paper is organized as follows. In Sec. 2, formulations for the multisubsystem co-design problem with relevant assumptions are detailed. In Sec. 3, the direct collocation method and the optimality condition decomposition are briefly introduced, and a decentralized implementation of the direct collocation method is described. Next, the multilevel and bilevel decentralized approaches are proposed in Sec. 4. The approaches are applied to a number of numerical and engineering examples, and the results are shown in Sec. 5. Finally, conclusion is presented in Sec. 6.

2 Problem Description

The co-design optimization problem considered is composed of N interconnected subsystems. As an example, Fig. 1 shows the structure of such a problem with three subsystems. In this paper, the ith subsystem is denoted as SSi. Each subsystem is formulated as a bi-objective optimization problem, with an objective function for its plant design (block “P”) as well as its control (block “C”) part.

Fig. 1
Problem structure—three coupled subsystems with plant (P) and control (C)
Fig. 1
Problem structure—three coupled subsystems with plant (P) and control (C)
Close modal

Four types of variables are considered in this problem, as shown in Fig. 1: local plant design variables yi, which exist only in the ith subsystem; shared plant design variables ysi,j, which exist in both the ith and jth subsystems; and state variables xi(t) and control variables ui(t) of the ith subsystem, which are time-variant variables in the control part. Subsystems are coupled in three ways. First, shared plant design variables are to have equal values in the related subsystems. Second, local and shared plant design variables can be coupled through the inequality and equality constraints on the physical plant. Third, state, control, and plant design variables are coupled through dynamic equations. A set V is defined such that if variables of the jth subsystem exist in the dynamic equations of the ith subsystem, then the integer pair (i, j) ∈ V; otherwise, (i,j)V.

A decentralized formulation of the ith subsystem for a multisubsystem co-design problem is as follows:
minui,yi,ysi,jZi=wPizPi(yi,ysi,j)+wCi0tfzCi(xi(t),ui(t))dtsubjecttotheconstraints,gi(yi,ysi,j)0,hi(yi,ysi,j)=0x˙i(t)=fi(xi,ui,xj,yi,ysi,j)=fii(xi,ui,yi,ysi,j)+fij(xj,yi,ysi,j)pi(xi(t),ui(t))0,qi(xi(0),xi(tf))=0(i=1,,N)
(1)
The objective function of the ith subsystem is a weighted sum of the plant objective, zPi, and control objective, zCi, for which the associated weights are, respectively, wPi and wCi, with wPi[0,1], wCi[0,1], and wPi+wCi=1 for all i ∈ {1, …, N}. The objective function of the entire co-design problem is to minimize the sum of all subsystems’ objective functions, i.e., Z=i=1NZi. gi and hi are physical plant constraints, x˙i(t)=fi represents the dynamic constraints, and pi and qi are path and control constraints, respectively. According to this formulation, the following assumptions are considered:
  • (A1) The coupling is unidirectional [26], i.e., the plant design objective functions and constraints are dependent only on plant design variables, while the control objective functions, dynamics, and constraints are dependent on both plant design and control variables.

  • (A2) In the plant design part of all subsystems, the objective functions, and inequality constraints are convex and equality constraints are linear functions of local and shared plant design variables. This assumption is made to satisfy the requirement of the dual decomposition method [25].

  • (A3) There exists an independent group of controllers (can be one or many) in each subsystem, with an overall time horizon [0, tf]. The couplings in the dynamics are additively separable (see Eq. (1)), while the control constraints are not coupled among subsystems. The control objective functions, dynamics, and state constraints can be nonlinear and nonconvex.

  • (A4) The optimal control subproblems of all subsystems are feasible, regardless of the values of plant design variables. It is also possible to explicitly impose the controllability constraints in the proposed optimization approach. However, given the format of the co-design problem as Eq. (1), the controllability constraints are not considered in this paper.

3 Decentralized Implementation of Direct Collocation Method

In this section, the key formulations of a decentralized implementation of the direct collocation method [17] using the optimality condition decomposition [27] is described. The detailed derivation can be found in the  Appendix. The implementation is then utilized in both the multilevel and the bilevel approaches proposed in Sec. 4.

Consider the control part of the co-design problem (Eq. (1)) without the plant design variables:
minuii=1NwCi0tfzCi(xi(t),ui(t))dtsubjecttotheconstraints,x˙i(t)=fi(xi,ui,xj)pi(xi(t),ui(t))0,qi(xi(0),xi(tf))=0(i=1,,N)
(2)

With the direct method, the entire time horizon [0, tf] is discretized into M subintervals by grid points, 0 = t0 < t1 < · · · < tM = tf, for all subsystems. The time difference s between two subsequent grid points is considered to be constant, i.e., a uniform grid with s = tktk−1 = tf /M, (k = 1, …, M). The control and state variables at all grid points are the optimizing variables in the converted problem. For convenience, state and control variables of SSi at tk are denoted as xi,k and ui,k.

Following the direct collocation method, the control variables between two grid points are approximated by linear interpolation, and the state variables are approximated by a cubic polynomial function. It can be shown [17] that the largest difference between the derivative of the approximation and that of the actual curve occurs at the midpoint of each interval, which is called a collocation point. Hence, a defect function ϕi,k for interval k in SSi is defined as follows:
ϕi,k=x˙i,k+1/2fi(xi,k+1/2,ui,k+1/2,xj,k+1/2)
(3)
and is constrained as ϕi,k=0 in the converted problem to enforce that the approximating curve of each state variable matches with the actual one. Furthermore, following Ref. [17], ϕi,k can be expressed as follows:
ϕi,k=3(xi,k+1xi,k)2s14(fi,k+fi,k+1)fi(xi,k+1/2,ui,k+1/2,xj,k+1/2)=0
(4)
The variables ui,k+1/2 and xi,k+1/2 of SSi at tk+1/2 can be evaluated from xi,k and ui,k. The evaluation of values of xj,k+1/2 starts from the equations of variables at the collocation points:
xj,k+1/2=12(xj,k+xj,k+1)+s8(fj,kfj,k+1)
(5)
Since x˙j=fj, and with a two-point approximation of this differential equation [28],
fj,kfj,k+1={2xj,k+1xj,kxj,k+22sk=1xj,k+1+xj,kxj,k+2xj2s1<k<M12xj,kxj,k+1xj,k12sk=M
(6)
Hence, xj,k+1/2 can be evaluated by plugging Eq. (6) into Eq. (5). If Xi is defined as the vector for the state and control variables at all grid points for SSi:
Xi=(ui,0T,,ui,MT,xi,0T,,,xi,MT)T
(7)
then ϕi,k can be represented as a function of Xi and Xj, i.e., ϕi,k(Xi,Xj)=0. Hence, the control problem (Eq. (2)) can be converted as a nonlinear program:
minXii=1NwCiz¯Ci(Xi)subjecttotheconstraints,ϕi,k(Xi,Xj)=0pi(Xi)0,qi(Xi)=0(i=1,,N)
(8)
where z¯Ci is the converted control objective function, in which the trapezoidal rule is used to approximate the integrals in the original objective function.
The converted problem (Eq. (8)) can be solved by the optimality condition decomposition (OCD) method [27], which is an iterative bilevel approach. The decomposed subproblem SSi (i = 1, …, N) is formulated as follows:
minXiwCiz¯Ci(Xi)+j=1,jiNϱ^jTϕj(Xi)subjecttotheconstraints,ϕi(Xi)=0,pi(Xi)0,qi(Xi)=0
(9)
where Xi=(X^1,,X^i1,Xi,X^i+1,,X^N), ϱi are the Lagrange multipliers associated with the complicating (or coupling) constraints ϕi, and the “^” (hat) notation stands for some fixed value of the variables.

In OCD, first the values X^i and ϱ^i are initialized. At the bottom level, Xi and ϱi can be optimized independently in SSi for i = 1, …, N, as shown in Eq. (9). At the top level, the difference of the optimized values Xi and the previous values X^i is calculated. If the difference is sufficiently small, the iteration stops. Otherwise, X^i and ϱ^i are updated with Xi and ϱi, respectively, and the procedure continues to the next iteration.

4 Decentralized Approaches for Multisubsystem Co-Design Problems

In this section, a multilevel and a bilevel decentralized framework to solve multisubsystem co-design problems are described. The direct collocation method [17] and the optimality condition decomposition method [27] are used in both frameworks.

4.1 Approach (A1): Multilevel Approach.

To solve the multisubsystem co-design problem (Eq. (1)) in a decentralized manner, both the plant and control parts are decomposed and solved. A similar multilevel framework as in Ref. [8] is adopted, as shown in Fig. 2. In the plant design part, as in Ref. [8], a dual decomposition approach [25] is used to solve for the plant design variables. In the control part, a decentralized implementation of the direct collocation method presented in Sec. 3 is employed.

Fig. 2
Decentralized approach (A1)
Fig. 2
Decentralized approach (A1)
Close modal

In the multilevel framework, the plant and control parts are integrated by gradients of the Hamiltonian of the control part with respect to plant design variables at y^i and y^si,j. By using a direct method, the state and control variables are discretized; hence, the gradients are calculated by summation functions of variable values at the grid points (see step (S3)).

4.1.1 Approach Steps.

The detailed steps, (S1)(S6), are as follows:

  • (S1) Initialize the plant design variable values yi0 and ysi,j0 of the ith subsystem. Set the current variable values, y^i and y^si,j, equal to yi0 and ysi,j0 for all i.

  • (S2) Apply the decentralized implementation of the direct collocation method in Sec. 3. Define Xi for the control and state variables of SSi at all grid points:
    Xi=(ui,0T,,ui,MT,xi,0T,,xi,MT)T
    (10)
    The variable and multiplier values X^i and ϱ^i, respectively, are initialized. At the bottom level of the control part, y^i and y^si,j are the current (fixed) values of the plant design variables. The SSi is then solved,
    minXiwCiz¯Ci(Xi)+(j,i)Vk=0Mϱ^j,kTϕj,k(Xi,X^j,y^j,y^sj,i)subjecttotheconstraints,ϕi,k(Xi,X^j,y^i,y^si,j)=0pi(Xi)0,qi(Xi)=0
    (11)
    to obtain Xi and ϱi for all SSi, where Xi and ϱi denote the optimal results of Xi and ρi, respectively, in this step. At the top level, the difference between Xi versus X^i is evaluated,
    e=i=1N(XiX^i2)
    (12)
    If e is greater than some tolerance value ϵD, then X^i and ϱ^i are updated with Xi and ϱi, and the iteration continues. Otherwise, the approach is considered to be converged for the control block. The approach continues to step (S3). The optimal variables and multipliers of the control part in the ith subsystem are denoted as Xi and ϱi, respectively.
  • (S3) The gradient of the Hamiltonian of the control part with respect to plant design variables at y^i and y^si,j can be obtained as follows:
    Hyi|y^i,y^si,j=i=1NwCi[k=0M1λi,k+1/2fi,k+1/2(Xi)yi]
    (13)
    Hysi,j|y^i,y^si,j=i=1NwCi[k=0M1λi,k+1/2fi,k+1/2(Xi)ysi,j]
    (14)
    λi,k+1/2 represents the approximating value of co-state variable at collocation points tk+1/2,
    λi,k+1/2=32ψi,kρi,k(k=0,1,,M1)
    (15)
    where ψi,k=z¯Ci,k+1/2/2s represents the scaling factor, which depends on the discretization [17].
  • (S4) The dual variable vector is initialized: μ = 0. For the ith subsystem, the partial Lagrangian of the plant part is as follows:
    Li=wPizPi(yi,ysi,j)+ηiTgi(yi,ysi,j)+γiThi(yi,ysi,j)+μTysi,j
    (16)
    The plant design variables yi, ysi,j are then solved by the necessary optimality conditions in the ith subsystem:
    Liyi+Hyi|y^i,y^si,j=0,Liysi,j+Hysi,j|y^i,y^si,j=0Liγi=0gi(yi,ysi,j)0,ηi0,ηigi(yi,ysi,j)=0
    (17)
    where “” denotes the Hadamard product.
  • (S5) Let y¯s be the vector of all ysi,j obtained in all subsystems. The consistency constraints can be represented as y¯s=Eσ [25], where E is a matrix with binary entries. Eij = 1 when a plant variable is shared between SSi and SSj, and 0 otherwise. For the lth iteration in the plant part, the approach computes the common values of shared plant design variables and evaluates the consistency constraint residual by
    σ=(ETE)1ETy¯s
    (18)
    e=y¯sEσ2
    (19)
    If e>ϵR, it updates the dual variable vector via a subgradient method with step size β,
    μ=μ+βy¯s
    (20)
    and repeat from Eq. (16). Otherwise, the results in the plant part are obtained as yinew and ysi,jnew. Calculate the overall objective function value Znew accordingly and continue to the next step.
  • (S6) If ZnewZcurrent2>ϵZ, update y^i, y^si,j with the new plant design values yinew, ysi,jnew, and repeat from step (S2). Otherwise, the approach has converged and the solution is obtained.

4.2 Approach (A2): Bilevel Approach.

In this framework, as shown in Fig. 3, the entire problem is decomposed in terms of subsystems and co-design of each subsystem. The solution steps are described on both cases when shared plant variables ysi,j exist or do not exist.

Fig. 3
Decentralized approach (A2)
Fig. 3
Decentralized approach (A2)
Close modal

4.2.1 Approach Steps (No Shared Plant Variables).

First, if there is no shared plant design variable, the decentralized implementation in Sec. 3 can be directly used. The bilevel framework is adopted from the optimality condition decomposition. At the bottom level, the co-design problem of each subsystem is solved by employing the direct collocation method. The optimizing variables are coordinated at the top level. The iterations stop when all variables converge.

The detailed steps (S1)(S3) are described as in the following:

  • (S1) For each subsystem at the bottom level, define vector Xi for the optimizing variables in SSi as follows:
    Xi=(ui,0T,,ui,MT,xi,0T,,xi,MT,yiT)T
    (21)
    which is the combining vector of the local control and state variables at all grid points, as well as the local plant design variables. Initialize the vectors of optimizing variables X^i and the vectors of Lagrange multipliers ϱ^i,k for all SSi.
  • (S2) By the decentralized implementation in Sec. 3, the formulation of SSi can be expressed as follows:
    minXiwPizPi(Xi)+wCiz¯Ci(Xi)+(j,i)Vk=0Mϱ^j,kTϕj,k(Xi,X^j)subjecttotheconstraints,gi(Xi)0,hi(Xi)=0ϕi,k(Xi,X^j)=0pi(Xi)0,qi(Xi)=0
    (22)
    where z¯Ci is the converted control objective function. Xi and ϱi,k are then solved, denoting the optimal optimizing variables and the Lagrange multipliers associated with the complicating constraints ϕi,k for all SSi.
  • (S3) At the top level, check the difference of the previous optimizing variable values X^i and the current optimized values Xi. If
    i=1N(X^iXi)2<ϵD
    (23)
    the converged co-design results are obtained and the iteration stops. Otherwise, update X^i and ϱ^i,k with Xi and ϱi,k, and repeat from step (S2).

4.2.2 Approach Steps (With Shared Plant Variables).

When there are shared plant variables ysi,j(i) between subsystems, both complicating variables and complicating constraints exist in the converted optimization problem by the direct method. Hence, a bilevel framework incorporating the optimality condition decomposition and the dual decomposition is devised. At the bottom level, the local state, control, and plant design variables as well as a local copy of the shared plant variables are optimized in each subsystem. At the top level, the dual variables associated with the shared plant variables are updated. The iterations continue until the shared plant variable values converge.

The detailed steps (S1)(S3) are described as in the following:

  • (S1) Similarly as in Sec. 4.2.1, in each subsystem at the bottom level, define vector Xi for the variables in the ith subsystem as follows:
    Xi=(ui,0T,,ui,MT,xi,0T,,xi,MT,yiT)T
    (24)
    In addition, the copy of the shared plant design variables is denoted as ysi,j(i) in the SSi. Initialize the vectors of optimizing variables X^i, the vectors of Lagrange multipliers ϱ^i,k, and the dual variable vectors νi,j for all subsystems.
  • (S2) By the decentralized implementation in Sec. 3, the formulation of SSi (i = 1, …, N) can be expressed as follows:
    minXi,ysi,j(i)wPizPi(Xi,ysi,j(i))+wCiz¯Ci(Xi)+(i,j)Vν¯i,jTysi,j(i)+(j,i)Vk=0Mϱ^j,kTϕj,k(Xi,X^j,y^si,j(j))subjecttotheconstraints,gi(Xi,ysi,j(i))0,hi(Xi,ysi,j(i))=0ϕi,k(Xi,X^j,ysi,j(i))=0pi(Xi)0,qi(Xi)=0
    (25)
    where ν~i,j=νi,j if i < j, and ν~i,j=νi,j if i > j. Hence, Xi and ysi,j(i) for the ith subsystem can be solved.
  • (S3) At the top level, the consistency constraints ysi,j(i)=ysi,j(j) are checked for the shared plant design variables. If
    (i,j)V(ysi,j(i)ysi,j(j))2<ϵys
    (26)
    the shared plant variable values have converged, and the co-design results are obtained. Otherwise, update the dual variable vectors νi,j by
    νi,j=νi,jβ(ysi,j(i)ysi,j(j))
    (27)
    where β is selected as a fixed step size value. Also, X^i and ϱ^i,k are updated with Xi and ϱi,k, and repeat from step (S2).

4.3 Notes on the Approaches.

For the proposed multilevel approach (A1), since the co-design problem is decomposed into both plant and control parts, it is assumed that there exists an objective function for both parts for each subsystem to ensure all the solution steps are valid. If the decomposition techniques are not applied in both parts, the approach has the same structure as the “nested” co-design strategy [1]. Hence, the proposed approach can be regarded as an extension of the nested strategy.

The bilevel approach (A2) aims to improve the performance of the multilevel approach. If no shared plant design variable exists, the convergence of the approach can be obtained by the convergence of the optimality condition decomposition method. In the presence of shared physical design variables, a bilevel approach is proposed in this paper to solve problems with both complicating variables and complicating constraints in a decomposed manner.

For both approaches, there is currently no convergence proof. The solution obtained can be a local optimum due to the nonconvexity property of the considered problem. On the implementation side, the proposed approaches can be implemented as black-box functions. While a guideline for selecting the tolerance values is not included, in general, if the tolerance values are made smaller, better optimized solutions are obtained, but then a longer computational time is required.

5 Examples

The decentralized approaches proposed in Sec. 4 are applied to a numerical and a spring-mass-damper system example in this section. In the numerical example, the dynamic system is linear and the objective function is quadratic; hence, the decentralized approach in Ref. [8] is applicable. In the engineering example, a three-subsystem co-design problem with nonlinear dynamics is considered, which the approach in Ref. [8] is not applied. Also, a scalable version of the engineering example is studied. Results by the proposed decentralized approaches are compared against the centralized solution.

5.1 Numerical Example.

The numerical example is slightly modified from Example 4.1 in Ref. [8]. The problem has two subsystems, and m1tom5 are the five physical plant design variables. The subsystems SS1 and SS2 are formulated as follows:

(SS1:)
minm1,m2,m3,m4,u1(t)Z1=wP1((m11)2+(m22)2+12(m31)2+12(m42)2)+12wC101(x1T(t)Q1x1(t)+R1u12(t))dt
subject to the constraints, g1(m1,m2,m3,m4):=m12+m22+m32+m4280
x˙1(t)=[2m1m12m2]x1(t)+[12]u1(t)+[0.01m400.010.02]x2(t)
(28)
x10=[1,0.1]T,Q1=diag(2,1),R1=1
(SS2:)
minm3,m4,m5,u2(t)Z2=wP2(12(m31)2+12(m42)2+(m53)2)+12wC201(x2T(t)Q2x2(t)+R2u22(t))dt
subject to the constraints, h2(m3,m4,m5):=m3+m4+2m58=0
x˙2(t)=[m50.52m5]x2(t)+[25]u2(t)+[0.020.0100.01m3]x1(t)
(29)
x20=[1,0.5]T,Q2=diag(1,2),R2=2

The weights in the objective function are set as wP1=wP2=wC1=wC2=0.5. Local design variables in the two subsystems are defined as y1 = [m1, m2]T and y2 = m5, respectively. The shared design variable between the subsystems ys1,2 = ys2,1 = [m3, m4]T. It can be verified that for all values of m1m5, the two subsystems are controllable.

The problem is then solved by the proposed decentralized approaches, as well as the centralized (or all-at-once) method by tomlab [29], which solves for the variables from all subsystems simultaneously. Since the dynamic equations is linear in the example, it can also be solved with the decentralized approach proposed in Ref. [8]. All the tolerance values are selected as 0.01.

5.1.1 Optimization Results and Discussion.

The optimal solutions by the centralized and decentralized approaches are shown in Table 1. It is shown that the optimal solution obtained by each of the decentralized approach is close to the centralized solution. In particular, the solution by the bilevel approach is exactly identical to the centralized one.

Table 1

Example 1: optimal solutions

y1*y2*ys1,2*Z*
Centralized[1.11, 1.80]T2.75[0.79, 1.70]T0.91
Decentralized [8][1.13, 1.78]T2.75[0.79, 1.71]T0.93
Decentralized (A1)[1.14, 1.78]T2.75[0.79, 1.71]T0.93
Decentralized (A2)[1.11, 1.80]T2.75[0.79, 1.70]T0.91
y1*y2*ys1,2*Z*
Centralized[1.11, 1.80]T2.75[0.79, 1.70]T0.91
Decentralized [8][1.13, 1.78]T2.75[0.79, 1.71]T0.93
Decentralized (A1)[1.14, 1.78]T2.75[0.79, 1.71]T0.93
Decentralized (A2)[1.11, 1.80]T2.75[0.79, 1.70]T0.91

In Fig. 4, the optimal control variables of SS1 and SS2 are plotted. As is observed, the decentralized optimal control solutions are again close to the centralized one, which demonstrates the validity of the proposed approaches. Finally, various feasible and infeasible initial points of the physical variables are used in the simulation, and the approximately identical converged results are obtained for all tests.

Fig. 4
Example 1: optimal controllers: (a) centralized, (b) decentralized: Ref. [8], (c) decentralized: multilevel, and (d) decentralized: bilevel
Fig. 4
Example 1: optimal controllers: (a) centralized, (b) decentralized: Ref. [8], (c) decentralized: multilevel, and (d) decentralized: bilevel
Close modal

5.2 Engineering Example: Spring-Mass-Damper System Model.

In this section, a similar engineering example in Ref. [8] is considered, as shown in Fig. 5; however, here, nonlinearity is introduced into dynamic equations. The plant design part of each subsystem consists of one mass (m), one spring (k), and one damper (c). The control part has a quadratic objective function. The state variable xi(t) = [xi1(t), xi2(t)]T denotes the displacement and velocity of the mass mi, and the control variable ui(t) is the force applied to the mass from 0 to 5 s.

Fig. 5
Example 2: spring-mass-damper system model
Fig. 5
Example 2: spring-mass-damper system model
Close modal
The formulation of the plant design part has been adopted from the spring design optimization formulation [30]. The wire diameters of springs are selected as the plant design variables, i.e., yi = di for all i. The objective function in the ith subsystem, zPi, is to maximize energy storage capacity. It is shown in Ref. [30] that if the spring index C is fixed at the maximum value, then the optimal spring design can be found by solving the following optimization problem:
minyizPi=(10.24G)(π12.8Fu)2C2yi4subjecttotheconstraints,C1κ1C1yi11(Insidediameter)κ2yi11(Lowerboundonwirediameter)0.8FuDsGC3yi21(Clashallowance)
(30)
where G is shear modulus (MPa), Fu is maximum allowable force (N), C is spring index (dimensionless), κ1 is minimum allowable inside diameter (m), κ2 is the lower bound on the plant design variables (m), and Ds is a clearance constant (dimensionless). The spring constants can be calculated as [4] follows:
ki=di4G8Di3Na(1+12C2)
(31)
where Di = Cdi is the coil diameter and Na is the number of active coils of the spring. With the proposed decentralized approaches using the direct collocation method, nonlinear dynamic equations can be handled in the co-design problems. In this example, a nonlinear force versus deflection for the springs is considered. In other words, it is assumed that each spring is increasingly less resilient as being elongated or depressed; hence, the response force of the ith spring can be expressed as follows:
Fki=kixi1+θxi13
where θ=0.5 is selected in this example. Then, the co-design problem is formulated as follows:
minyi,ui(t)Z=i=13wPi(10.24G)(π12.8Fu)2C2yi4+12i=13wCi05(xiT(t)Qixi(t)+Riui2(t))dtsubjecttotheconstraints,(C11)yiκ1C10,κ2yi00.8FuDsGC3yi21,xi0=[1,1]Tx˙11=x12,x˙21=x22,x˙31=x32u1=m1x˙12+c1x12+c2(x12x22)+k1x11+k2(x11x21)θx113θ(x11x21)3u2=m2x˙22+c2(x22x12)+c3(x22x32)+k2(x21x11)+k3(x21x31)θ(x21x11)3θ(x21x31)3u3=m3x˙32+c3(x32x22)+k3(x31x21)θ(x31x21)3Qi=diag(1,1),Ri=1(i=1,2,3)
(32)

The parameters are selected as c1 = 10, c2 = 15, c3 = 20, m1 = m2 = m3 = 5, G = 30 (for chrome silicon), C = 8, Fu = 6, κ1=0.05, κ2=0.1, Ds = 0.5, and Na = 200. All the tolerance values are selected as 0.01.

5.2.1 Optimization Results and Discussion.

The problem is solved by the decentralized approaches and by the centralized approach using tomlab [29]. The solutions for the optimal wire diameters are listed in Table 2. It is shown that the optimal solution of the bilevel approach is identical to that of the centralized method. The optimal solution of the multilevel approach is close but slightly worse than the other two. The profiles of optimal controllers ui(t) obtained by all approaches are shown in Fig. 6. As is observed, the optimal control variables of the decentralized approaches are comparable with the ones of the centralized approach.

Fig. 6
Example 2: optimal controllers: (a) centralized, (b) decentralized: multilevel, and (c) decentralized: bilevel
Fig. 6
Example 2: optimal controllers: (a) centralized, (b) decentralized: multilevel, and (c) decentralized: bilevel
Close modal
Table 2

Example 2: optimal solutions

y1*y2*y3*Z*
Centralized0.130.100.1010.14
Decentralized (A1)0.100.100.1010.58
Decentralized (A2)0.130.100.1010.14
y1*y2*y3*Z*
Centralized0.130.100.1010.14
Decentralized (A1)0.100.100.1010.58
Decentralized (A2)0.130.100.1010.14

5.2.2 Computational Cost for Increasing Number of Subsystems.

In this section, the previous spring-mass-damper model is reformulated as a scalable problem, as shown in Fig. 7. For the simplicity of the implementation, linear dynamic equations are used in this case, with all masses mi = 5 and damping coefficients ci = 10. The objective function of the plant design part of each subsystem is assumed to be (yi − 0.1)2. There exists a control input for each mass. The number of subsystems considered for each test problem is 5, 10, and 20–80. Each test problem is solved in both the centralized (by tomlab [29]) and the proposed decentralized approaches.

Fig. 7
Example: series of spring-mass-damper subsystems
Fig. 7
Example: series of spring-mass-damper subsystems
Close modal

Figure 8 shows a comparison of the computational (CPU) time for the test problems when no communication costs and delays are considered. The centralized results are given by the CPU time obtained from tomlab [29], which, as expected, grow nonlinearly as the number of variables is increased.

Fig. 8
Comparison of computational time between centralized and decentralized approaches
Fig. 8
Comparison of computational time between centralized and decentralized approaches
Close modal

By using the decentralized approach, the results show that the computational time increases about linearly as the number of subsystems is increased. However, if all steps are simulated on a single machine, i.e., subsystems are optimized sequentially one after another, the computational time values are greater than the centralized ones, as shown in Fig. 8. Hence, a “parallelization” procedure is simulated, which assumes some number of “machines” can be used simultaneously. The machines are used to solve the decomposed subproblems in a parallel manner. As an example, all the decomposed subproblems are divided into groups of ten machines, so each group can be solved “in parallel”. The largest value of the running time among the ten machines is measured. By this time, all other nine machines in the same group are assumed to have finished the tasks of optimizing the corresponding subproblems. Following this procedure, a second computational time plot is presented for this test problem, which shows that both the multilevel and the bilevel decentralized approaches are able to perform better when the number of subsystems increases. In addition, it is shown that the computational time of the bilevel approach is better than the multilevel one.

In a parallel computing environment, communication costs and delays may be important factors to consider; therefore, a simulation is implemented to include communication costs among ten machines in this example. A constant value, simulating the time to transfer information together with delays, is added in each iteration for the decentralized approaches. In Fig. 9, the constant value is selected as 0.05 s and 2 s. As is observed, when the number of subsystems is large, the decentralized approaches can outperform the centralized one.

Fig. 9
Comparison of computational time between centralized and decentralized approaches with various communication costs
Fig. 9
Comparison of computational time between centralized and decentralized approaches with various communication costs
Close modal

Finally, the number of machines can also influence the computational time of the decentralized approaches. If more subsystems can be optimized in parallel simultaneously, the total computational time to obtain optimal solution can be reduced. In the last simulation, the number of “machines” is set as 10 and 20, and the time to transfer information with delays in each iteration is assumed to be 1 s. It is shown in Fig. 10 that as more machines are used at the subsystem level, the total computational time can be further reduced. On the contrast, despite that there is no communication cost during the optimization process, a centralized approach lacks the capability of utilizing more computational power when the problem has large number of variables. In all the tests, the crossing points of the computational time curves can provide useful information for making decisions.

Fig. 10
Comparison of computational time between centralized and decentralized with different computational power and constant communication costs
Fig. 10
Comparison of computational time between centralized and decentralized with different computational power and constant communication costs
Close modal

6 Conclusion

In this paper, two new decomposition-based numerical approaches are developed to solve more general multisubsystem co-design problems. The direct collocation method is used to discretize the continuous-time dynamic constraints and is incorporated with the optimality condition decomposition method to solve the problems in a decentralized manner. The proposed approaches are applied to a numerical and an engineering example, and the solutions are compared with the centralized one.

Finally, a scalable test problem is considered to show a potential advantage of the proposed decentralized approaches. It is shown that by simulating multiple parallel machines to solve the decomposed subproblems, the decentralized approaches take a shorter computational time than the centralized approach as the size of the problem increases. In addition, it can be observed that the bilevel approach outperforms the multilevel one in terms of both solution quality and the computational time.

Acknowledgment

This work was supported by the Office of Naval Research (ONR) (Grant No. N000141310160). Such support does not constitute an endorsement by the funding agency of the opinions expressed in the paper. We would also like to thank all reviewers for their constructive and helpful comments and suggestions.

Nomenclature

e=

norm of difference of two vectors

f=

dynamic equations of the system

g=

inequality constraints for physical plant design variables

h=

equality constraints for physical plant design variables

p=

vector of control/path constraints

q=

vector of initial or final conditions of state variables

s=

time difference between two grid points

u=

vector of control variables

x=

state variables

E=

consistency constraint matrix

H=

Hamiltonian of control subproblems

L=

Lagrangian of physical plant design optimization subproblems

M=

number of subintervals in discretization

N=

number of subsystems

V=

set of pairs of connected subsystems

X=

combined optimizing variable vector

Z=

system-level objective function

tf=

final time

yi=

local physical plant design variable vector in the ith subsystem

y¯s=

vector of combined shared physical plant design variables

ysi,j=

shared physical plant design variable vector in ith and jth subsystems

wC, wP=

weighting coefficient associated with control and physical plant design objective function

zC, zP=

control and physical plant design objective function

( · )i=

symbol ( · ) with subscript i: variables/parameters in the ith subsystem

( · )k=

symbol ( · ) with subscript k: variables/parameters at grid point k

β=

step size constants

γ,η=

Lagrange multipliers of physical plant design constraints

ϵD,ϵZ,ϵR,ϵys=

preset tolerance value of convergence for optimizing variables, objective function, consistency constraint residual, and shared plant design variables

λ=

co-state variables

μ=

dual variable vectors associated with physical plant design consistency constraints

ν=

dual variable vectors in bilevel approach

ρ=

Lagrange multipliers associated with converted complicating constraints

σ=

common value vector of shared physical plant design variables

ϕ=

converted complicating constraints at collocation points

ψ=

scaling factors

Appendix: Derivation of the Decentralized Direct Collocation Method

In this paper, the discretization steps described in Ref. [17] are modified to fit a decentralized formulation.

In the direct collocation method, as shown in Fig. 11, the value of control variables between every two grid points are approximated as a linear interpolating function, and the values of state variables are approximated by a cubic polynomial function, known as cubic collocation at Lobatto points [17]. Hence, for tkttk+1,
ui(t)=ui,k+ttks(ui,k+1ui,k)
(A1)
xi(t)=l=03ci,k(l)(ttks)l
(A2)
where ci,k(l)(l=0,1,2,3) are coefficients of the cubic polynomial,
ci,k(0)=xi,kci,k(1)=sfi,kci,k(2)=3xi,k2sfi,k+3xi,k+1sfi,k+1ci,k(3)=2xi,k+sfi,k2xi,k+1+sfi,k+1fi,k:=fi(xi,k,ui,k,xj,k,yi,ysi,j)
(A3)
It is then proved [17] that by finding the coefficients ci,k(l), the dynamic equations are satisfied at all grid points. By matching the slopes of the cubic polynomials at the midpoints (or collocation points) at all discretized intervals, the approximating functions are regarded as the solution to the optimal control problem. In the subinterval tkttk+1, the variables at the collocation point are denoted with the subscript (k + 1/2) and are evaluated as [17] (k = 0, …, M − 1),
ui,k+1/2=12(ui,k+ui,k+1)xi,k+1/2=12(xi,k+xi,k+1)+s8(fi,kfi,k+1)x˙i,k+1/2=3(xi,k+1xi,k)2s14(fi,k+fi,k+1)
(A4)
Hence, the continuous-time problem is transformed to a nonlinear optimization problem, whose constraints include (1) the slopes of approximating polynomials matching with the value of dynamic equations at the midpoint of each discretized interval, (2) the control constraints for variables at all grid points, and (3) the boundary conditions at t0 = 0 and tM = tf. The constraints are formulated as follows:
x˙i,k+1/2=fi(xi,k+1/2,ui,k+1/2,xj,k+1/2)
(A5)
pi(xi,k,ui,k)0,qi(xi,1,xi,M,tM)=0
(A6)
A defect function ϕi,k for interval k in subsystem i is defined as follows:
ϕi,k=x˙i,k+1/2fi(xi,k+1/2,ui,k+1/2,xj,k+1/2)
(A7)
Therefore, the constraints of Eq. (A5) can be rewritten as ϕi,k=0. Furthermore, by plugging in the variables in Eq. (A4), ϕi,k can be expressed as follows:
ϕi,k=3(xi,k+1xi,k)2s14(fi,k+fi,k+1)fi(xi,k+1/2,ui,k+1/2,xj,k+1/2)=0
(A8)
Fig. 11
Direct collocation method discretization: (a) control—piecewise linear and (b) state—piecewise cubic polynomial
Fig. 11
Direct collocation method discretization: (a) control—piecewise linear and (b) state—piecewise cubic polynomial
Close modal

Note that ϕi,k is related to the variable values of SSj at the grid and the collocation points, i.e., xj,k, xj,k+1/2, and xj,k+1.

References

1.
Fathy
,
H. K.
,
Reyer
,
J. A.
,
Papalambros
,
P. Y.
, and
Ulsoy
,
A. G.
,
2001
, “
On the Coupling Between the Plant and Controller Optimization Problems
,”
Proceedings of the 2001 American Control Conference
, Vol.
3
,
Arlington, VA
,
June 25–27
, pp.
1864
1869
.
2.
Peters
,
D. L.
,
Papalambros
,
P.
, and
Ulsoy
,
A.
,
2011
, “
Control Proxy Functions for Sequential Design and Control Optimization
,”
ASME J. Mech. Des.
,
133
(
9
), p.
091007
. 10.1115/1.4004792
3.
Allison
,
J. T.
, and
Nazari
,
S.
,
2010
,
“Combined Plant and Controller Design Using Decomposition-Based Design Optimization and the Minimum Principle
,” Proceedings of the ASME 2010 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, ASME Paper No.
DETC2010-28887
. 10.1115/detc2010-28887
4.
Allison
,
J. T.
,
Guo
,
T.
, and
Han
,
Z.
,
2014
, “
Co-Design of an Active Suspension Using Simultaneous Dynamic Optimization
,”
ASME J. Mech. Des.
,
136
(
8
), p.
081003
. 10.1115/1.4027335
5.
Jiang
,
Y.
,
Wang
,
Y.
,
Bortoff
,
S. A.
, and
Jiang
,
Z.-P.
,
2015
, “
Optimal Codesign of Nonlinear Control Systems Based on a Modified Policy Iteration Method
,”
IEEE Trans. Neural Networks Learn. Syst.
,
26
(
2
), pp.
409
414
. 10.1109/TNNLS.2014.2382338
6.
Li
,
H.
,
Yan
,
W.
, and
Shi
,
Y.
,
2018
, “
Triggering and Control Codesign in Self-Triggered Model Predictive Control of Constrained Systems: With Guaranteed Performance
,”
IEEE Trans. Autom. Control
,
63
(
11
), pp.
4008
4015
. 10.1109/TAC.2018.2810514
7.
Xu
,
Y.
,
Cervin
,
A.
, and
Arzen
,
K. E.
,
2018
, “
Jitter-Robust LQG Control and Real-Time Scheduling Co-Design
,”
Proceedings of the American Control Conference
,
Milwaukee, WI
,
June 27–29
, pp.
3189
3196
http://dx.doi.org/10.23919/ACC.2018.8430953.
8.
Liu
,
T.
,
Azarm
,
S.
, and
Chopra
,
N.
,
2017
, “
On Decentralized Optimization for a Class of Multisubsystem Codesign Problems
,”
ASME J. Mech. Des.
,
139
(
12
), p.
121404
. 10.1115/1.4037893
9.
Behtash
,
M.
, and
Alexander-Ramos
,
M. J.
,
2018
, “
Decomposition-Based MDSDO for Co-Design of Large-Scale Dynamic Systems
,”
ASME 2018 International Design Engineering Technical Conferences
,
Quebec City, Quebec, Canada
,
Aug. 26–29
,
ASME Paper No. DETC2018-85935
http://doi.org/10.1115/DETC2018-85935.
10.
Betts
,
J. T.
,
1998
, “
Survey of Numerical Methods for Trajectory Optimization
,”
J. Guidance Control Dyn.
,
21
(
2
), pp.
193
207
. 10.2514/2.4231
11.
Rao
,
A. V.
,
2009
, “
A Survey of Numerical Methods for Optimal Control
,”
Adv. Astronaut. Sci.
,
135
(
1
), pp.
497
528
.
12.
Gerdts
,
M.
,
2003
, “
Direct Shooting Method for the Numerical Solution of Higher-Index DAE Optimal Control Problems
,”
J. Optim. Theory Appl.
,
117
(
2
), pp.
267
294
. 10.1023/A:1023679622905
13.
Diehl
,
M.
,
Bock
,
H. G.
,
Diedam
,
H.
, and
Wieber
,
P.-B.
,
2006
, “Fast Direct Multiple Shooting Algorithms for Optimal Robot Control,”
Fast Motions in Biomechanics and Robotics; Optimization and Feedback Control
,
M.
Diehl
and
K.
Mombaur
, eds.,
Springer
,
New York
, pp.
65
93
.
14.
Fahroo
,
F.
, and
Ross
,
I. M.
,
2002
, “
Direct Trajectory Optimization by a Chebyshev Pseudospectral Method
,”
J. Guidance Control Dyn.
,
25
(
1
), pp.
160
166
. 10.2514/2.4862
15.
Benson
,
D.
,
2005
, “
A Gauss Pseudospectral Transcription for Optimal Control
,”
Ph.D. dissertation
,
Massachusetts Institute of Technology
,
Cambridge, MA
.
16.
Francolin
,
C.
, and
Rao
,
A
,
2011
, “
Direct Trajectory Optimization and Costate Estimation of State Inequality Path-Constrained Optimal Control Problems Using a Radau Pseudospectral Method
,”
2012 Guidance, Navigation, and Control Conference
,
AIAA, Minneapolis, MN
,
Aug. 13–16
, p.
4528
.
17.
Von Stryk
,
D. M. O.
,
1993
, “Numerical Solution of Optimal Control Problems by Direct Collocation,”
Optimal Control: Calculus of Variations, Optimal Control Theory and Numerical Methods
,
R.
Bulirsch
,
A.
Miele
,
J.
Stoer
, and
K.
Well
, eds.,
Birkhäuser
,
Basel
, pp.
129
143
.
18.
Herber
,
D. R.
,
2017
, “
Advances in Combined Architecture, Plant, and Control Design
,”
Ph.D. dissertation
,
University of Illinois at Urbana-Champaign
,
Champaign, IL
.
19.
Chilan
,
C. M.
,
Herber
,
D. R.
,
Nakka
,
Y. K.
,
Chung
,
S. J.
,
Allison
,
J. T.
,
Aldrich
,
J. B.
, and
Alvarez-Salazar
,
O. S.
,
2016
, “
Co-Design of Strain Actuated Solar Arrays for Precision Pointing and Jitter Reduction
,”
AIAA J.
,
55
(
9
), pp.
3180
3195
. 10.2514/1.J055748
20.
Spielberg
,
A.
,
Araki
,
B.
,
Sung
,
C.
,
Tedrake
,
R.
, and
Rus
,
D.
,
2017
, “
Functional Co-Optimization of Articulated Robots
,”
IEEE International Conference on Robotics and Automation (ICRA)
,
Singapore
,
May 29–June 3
, pp.
5035
5042
http://dx.doi.org/10.1109/ICRA.2017.7989587.
21.
Geiger
,
B.
,
Horn
,
J.
,
DeLullo
,
A.
,
Niessner
,
A.
, and
Long
,
L.
,
2006
, “
Optimal Path Planning of UAVs Using Direct Collocation With Nonlinear Programming
,”
2006 Guidance, Navigation, and Control Conference
,
AIAA, Keystone, CO
,
Aug. 21–24
, p.
6199
.
22.
Stipanovic
,
D. M.
,
Inalhan
,
G.
,
Teo
,
R.
, and
Tomlin
,
C. J.
,
2004
, “
Decentralized Overlapping Control of a Formation of Unmanned Aerial Vehicles
,”
Automatica
,
40
(
8
), pp.
1285
1296
. 10.1016/j.automatica.2004.02.017
23.
Ghali
,
S.
,
Aousgi
,
I.
, and
Elloumi
,
S.
,
2017
, “
Decentralized Optimal Control of Interconnected Systems Using a Direct Approach
,”
2017 International Conference on Advanced Systems and Electrical Technologies
,
Hammamet, Tunisia
,
Jan. 14–17
, pp.
79
84
http://dx.doi.org/10.1109/ASET.2017.7983670.
24.
Mehraeen
,
S.
, and
Jagannathan
,
S.
,
2011
, “
Decentralized Optimal Control of a Class of Interconnected Nonlinear Discrete-Time Systems by Using Online Hamilton-Jacobi-Bellman Formulation
,”
IEEE Trans. Neural Networks
,
22
(
11
), pp.
1757
1769
. 10.1109/TNN.2011.2160968
25.
Boyd
,
S.
,
Xiao
,
L.
,
Mutapcic
,
A.
, and
Mattingley
,
J.
,
2007
,
Notes on Decomposition Methods
,
Stanford University
,
Stanford, CA
, https://web.stanford.edu/class/ee364b/lectures/decompositionnotes.pdf, Accessed Dec. 25, 2019.
26.
Peters
,
D. L.
,
Papalambros
,
P. Y.
, and
Ulsoy
,
A. G.
,
2009
, “
On Measures of Coupling Between the Artifact and Controller Optimal Design Problems
,” ASME 2009 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, ASME Paper No.
DETC2009-86868
. 10.1115/detc2009-86868
27.
Conejo
,
A. J.
,
Nogales
,
F. J.
, and
Prieto
,
F. J.
,
2002
, “
A Decomposition Procedure Based on Approximate Newton Directions
,”
Math. Program.
,
93
(
3
), pp.
495
515
. 10.1007/s10107-002-0304-3
28.
Moin
,
P.
,
2010
,
Fundamentals of Engineering Numerical Analysis
,
Cambridge University Press
,
Cambridge, UK
.
29.
Rutquist
,
P. E.
, and
Edvall
,
M. M.
,
2010
,
Propt-MATLAB Optimal Control Software
,
TOMLAB Optimization, Inc.
,
Pullman, WA
.
30.
Azarm
,
S.
, and
Papalambros
,
P.
,
1982
, “
An Interactive Design Procedure for Optimization of Helical Compression Springs
,”
University of Michigan
,
MI
,
Technical Report No. UMMEAM-82-7
.