## 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. [1–3]), 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 [14–16], and direct collocation method [17]. In co-design, the direct transcription method has been applied to single-system problems [4,18–20]. 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 *i*th 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.

Four types of variables are considered in this problem, as shown in Fig. 1: local plant design variables *y*_{i}, which exist only in the *i*th subsystem; shared plant design variables *y*_{s}_{i,j}, which exist in both the *i*th and *j*th subsystems; and state variables *x*_{i}(*t*) and control variables *u*_{i}(*t*) of the *i*th 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 *j*th subsystem exist in the dynamic equations of the *i*th subsystem, then the integer pair (*i*, *j*) ∈ *V*; otherwise, $(i,j)\u2209V$.

*i*th subsystem for a multisubsystem co-design problem is as follows:

*i*th 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\u2208[0,1]$, $wCi\u2208[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=\u2211i=1NZi$.

*g*

_{i}and

*h*

_{i}are physical plant constraints, $x\u02d9i(t)=fi$ represents the dynamic constraints, and

*p*

_{i}and

*q*

_{i}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,

*t*_{f}]. 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.

With the direct method, the entire time horizon [0, *t*_{f}] is discretized into *M* subintervals by grid points, 0 = *t*_{0} < *t*_{1} < · · · < *t*_{M} = *t*_{f}, for all subsystems. The time difference *s* between two subsequent grid points is considered to be constant, i.e., a uniform grid with *s* = *t*_{k} − *t*_{k−1} = *t*_{f} /*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 *t*_{k} are denoted as *x*_{i,k} and *u*_{i,k}.

*k*in

*SSi*is defined as follows:

*u*

_{i,k+1/2}and

*x*

_{i,k+1/2}of

*SSi*at

*t*

_{k+1/2}can be evaluated from

*x*

_{i,k}and

*u*

_{i,k}. The evaluation of values of

*x*

_{j,k+1/2}starts from the equations of variables at the collocation points:

*x*

_{j,k+1/2}can be evaluated by plugging Eq. (6) into Eq. (5). If

*X*

_{i}is defined as the vector for the state and control variables at all grid points for

*SSi*:

*X*

_{i}and

*X*

_{j}, i.e., $\varphi i,k(Xi,Xj)=0$. Hence, the control problem (Eq. (2)) can be converted as a nonlinear program:

*SSi*(

*i*= 1, …,

*N*) is formulated as follows:

_{i}are the Lagrange multipliers associated with the complicating (or coupling) constraints $\varphi i$, and the “^” (hat) notation stands for some fixed value of the variables.

In OCD, first the values $X^i$ and $\u03f1^i$ are initialized. At the bottom level, *X*_{i} 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 *X*_{i} and the previous values $X^i$ is calculated. If the difference is sufficiently small, the iteration stops. Otherwise, $X^i$ and $\u03f1^i$ are updated with *X*_{i} and ϱ_{i}, respectively, and the procedure continues to the next iteration.

## 4 Decentralized Approaches for Multisubsystem Co-Design Problems

### 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.

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

*y*_{i0}and*y*_{s}_{i,j0}of the*i*th subsystem. Set the current variable values, $y^i$ and $y^si,j$, equal to*y*_{i0}and*y*_{s}_{i,j0}for all*i*.- (S2) Apply the decentralized implementation of the direct collocation method in Sec. 3. Define
*X*_{i}for the control and state variables of*SSi*at all grid points:The variable and multiplier values $X^i$ and $\u03f1^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(10)$Xi=(ui,0T,\u2026,ui,MT,xi,0T,\u2026,xi,MT)T$*SSi*is then solved,to obtain $Xi\u2021$ and $\u03f1i\u2021$ for all(11)$minXiwCiz\xafCi(Xi)+\u2211(j,i)\u2208V\u2211k=0M\u03f1^j,kT\varphi j,k(Xi,X^j,y^j,y^sj,i)subjecttotheconstraints,\varphi i,k(Xi,X^j,y^i,y^si,j)=0pi(Xi)\u22640,qi(Xi)=0$*SSi*, where $Xi\u2021$ and $\u03f1i\u2021$ denote the optimal results of*X*_{i}and $\rho i$, respectively, in this step. At the top level, the difference between $Xi\u2021$ versus $X^i$ is evaluated,If(12)$e=\u2211i=1N(\Vert Xi\u2021\u2212X^i\Vert 2)$*e*is greater than some tolerance value $\u03f5D$, then $X^i$ and $\u03f1^i$ are updated with $Xi\u2021$ and $\u03f1i\u2021$, 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*i*th subsystem are denoted as $Xi\u2021$ and $\u03f1i\u2021$, 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:(13)$\u2202H\u2202yi|y^i,y^si,j=\u2211i=1NwCi[\u2211k=0M\u22121\lambda i,k+1/2\u2202fi,k+1/2(Xi\u2021)\u2202yi]$$\lambda i,k+1/2$ represents the approximating value of co-state variable at collocation points(14)$\u2202H\u2202ysi,j|y^i,y^si,j=\u2211i=1NwCi[\u2211k=0M\u22121\lambda i,k+1/2\u2202fi,k+1/2(Xi\u2021)\u2202ysi,j]$
*t*_{k+1/2},where $\psi i,k=z\xafCi,k+1/2/2s$ represents the scaling factor, which depends on the discretization [17].(15)$\lambda i,k+1/2=\u221232\psi i,k\rho i,k\u2021(k=0,1,\u2026,M\u22121)$ - (S4) The dual variable vector is initialized:
*μ*= 0. For the*i*th subsystem, the partial Lagrangian of the plant part is as follows:The plant design variables(16)$Li=wPizPi(yi,ysi,j)+\eta iTgi(yi,ysi,j)+\gamma iThi(yi,ysi,j)+\mu Tysi,j$*y*_{i},*y*_{s}_{i,j}are then solved by the necessary optimality conditions in the*i*th subsystem:where “$\u2218$” denotes the Hadamard product.(17)$\u2202Li\u2202yi+\u2202H\u2202yi|y^i,y^si,j=0,\u2202Li\u2202ysi,j+\u2202H\u2202ysi,j|y^i,y^si,j=0\u2202Li\u2202\gamma i=0gi(yi,ysi,j)\u22640,\eta i\u22650,\eta i\u2218gi(yi,ysi,j)=0$ - (S5) Let $y\xafs$ be the vector of all
*y*_{s}_{i,j}obtained in all subsystems. The consistency constraints can be represented as $y\xafs=E\sigma $ [25], where*E*is a matrix with binary entries.*E*_{ij}= 1 when a plant variable is shared between*SSi*and*SSj*, and 0 otherwise. For the*l*th iteration in the plant part, the approach computes the common values of shared plant design variables and evaluates the consistency constraint residual by(18)$\sigma =(ETE)\u22121ETy\xafs$If $e>\u03f5R$, it updates the dual variable vector via a subgradient method with step size $\beta $,(19)$e=\Vert y\xafs\u2212E\sigma \Vert 2$and repeat from Eq. (16). Otherwise, the results in the plant part are obtained as(20)$\mu =\mu +\beta y\xafs$*y*_{i−new}and*y*_{s}_{i,j−new}. Calculate the overall objective function value*Z*_{new}accordingly and continue to the next step. (S6) If $\Vert Znew\u2212Zcurrent\Vert 2>\u03f5Z$, update $y^i$, $y^si,j$ with the new plant design values

*y*_{i−new},*y*_{s}_{i,j−new}, 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.

#### 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
*X*_{i}for the optimizing variables in*SSi*as follows: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 $\u03f1^i,k$ for all(21)$Xi=(ui,0T,\u2026,ui,MT,xi,0T,\u2026,xi,MT,yiT)T$*SSi*. - (S2) By the decentralized implementation in Sec. 3, the formulation of
*SSi*can be expressed as follows:where $z\xafCi$ is the converted control objective function. $Xi\u2021$ and $\u03f1i,k\u2021$ are then solved, denoting the optimal optimizing variables and the Lagrange multipliers associated with the complicating constraints $\varphi i,k$ for all(22)$minXiwPizPi(Xi)+wCiz\xafCi(Xi)+\u2211(j,i)\u2208V\u2211k=0M\u03f1^j,kT\varphi j,k(Xi,X^j)subjecttotheconstraints,gi(Xi)\u22640,hi(Xi)=0\varphi i,k(Xi,X^j)=0pi(Xi)\u22640,qi(Xi)=0$*SSi*. - (S3) At the top level, check the difference of the previous optimizing variable values $X^i$ and the current optimized values $Xi\u2021$. Ifthe converged co-design results are obtained and the iteration stops. Otherwise, update $X^i$ and $\u03f1^i,k$ with $Xi\u2021$ and $\u03f1i,k\u2021$, and repeat from step ((23)$\u2211i=1N\Vert (X^i\u2212Xi\u2021)\Vert 2<\u03f5D$
**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
*X*_{i}for the variables in the*i*th subsystem as follows:In addition, the copy of the shared plant design variables is denoted as $ysi,j(i)$ in the(24)$Xi=(ui,0T,\u2026,ui,MT,xi,0T,\u2026,xi,MT,yiT)T$*SSi*. Initialize the vectors of optimizing variables $X^i$, the vectors of Lagrange multipliers $\u03f1^i,k$, and the dual variable vectors $\nu 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:where $\nu ~i,j=\nu i,j$ if(25)$minXi,ysi,j(i)wPizPi(Xi,ysi,j(i))+wCiz\xafCi(Xi)+\u2211(i,j)\u2208V\nu \xafi,jTysi,j(i)+\u2211(j,i)\u2208V\u2211k=0M\u03f1^j,kT\varphi j,k(Xi,X^j,y^si,j(j))subjecttotheconstraints,gi(Xi,ysi,j(i))\u22640,hi(Xi,ysi,j(i))=0\varphi i,k(Xi,X^j,ysi,j(i))=0pi(Xi)\u22640,qi(Xi)=0$*i*<*j*, and $\nu ~i,j=\u2212\nu i,j$ if*i*>*j*. Hence, $Xi\u2021$ and $ysi,j(i)\u2021$ for the*i*th subsystem can be solved. - (S3) At the top level, the consistency constraints $ysi,j(i)\u2021=ysi,j(j)\u2021$ are checked for the shared plant design variables. Ifthe shared plant variable values have converged, and the co-design results are obtained. Otherwise, update the dual variable vectors $\nu i,j$ by(26)$\u2211(i,j)\u2208V\Vert (ysi,j(i)\u2021\u2212ysi,j(j)\u2021)\Vert 2<\u03f5ys$where $\beta $ is selected as a fixed step size value. Also, $X^i$ and $\u03f1^i,k$ are updated with $Xi\u2021$ and $\u03f1i,k\u2021$, and repeat from step ((27)$\nu i,j=\nu i,j\u2212\beta (ysi,j(i)\u2021\u2212ysi,j(j)\u2021)$
**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 *SS*1 and *SS*2 are formulated as follows:

*SS*1:)

*SS*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 *y*_{1} = [*m*_{1}, *m*_{2}]^{T} and *y*_{2} = *m*_{5}, respectively. The shared design variable between the subsystems *y*_{s}_{1,2} = *y*_{s}_{2,1} = [*m*_{3}, *m*_{4}]^{T}. It can be verified that for all values of *m*_{1}–*m*_{5}, 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.

y_{1}* | y_{2}* | y_{s}_{1,2}* | Z* | |
---|---|---|---|---|

Centralized | [1.11, 1.80]^{T} | 2.75 | [0.79, 1.70]^{T} | 0.91 |

Decentralized [8] | [1.13, 1.78]^{T} | 2.75 | [0.79, 1.71]^{T} | 0.93 |

Decentralized (A1) | [1.14, 1.78]^{T} | 2.75 | [0.79, 1.71]^{T} | 0.93 |

Decentralized (A2) | [1.11, 1.80]^{T} | 2.75 | [0.79, 1.70]^{T} | 0.91 |

y_{1}* | y_{2}* | y_{s}_{1,2}* | Z* | |
---|---|---|---|---|

Centralized | [1.11, 1.80]^{T} | 2.75 | [0.79, 1.70]^{T} | 0.91 |

Decentralized [8] | [1.13, 1.78]^{T} | 2.75 | [0.79, 1.71]^{T} | 0.93 |

Decentralized (A1) | [1.14, 1.78]^{T} | 2.75 | [0.79, 1.71]^{T} | 0.93 |

Decentralized (A2) | [1.11, 1.80]^{T} | 2.75 | [0.79, 1.70]^{T} | 0.91 |

In Fig. 4, the optimal control variables of *SS*1 and *SS*2 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.

### 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 *x*_{i}(*t*) = [*x*_{i1}(*t*), *x*_{i2}(*t*)]^{T} denotes the displacement and velocity of the mass *m*_{i}, and the control variable *u*_{i}(*t*) is the force applied to the mass from 0 to 5 s.

*y*

_{i}=

*d*

_{i}for all

*i*. The objective function in the

*i*th 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:

*G*is shear modulus (MPa),

*F*

_{u}is maximum allowable force (N),

*C*is spring index (dimensionless), $\kappa 1$ is minimum allowable inside diameter (m), $\kappa 2$ is the lower bound on the plant design variables (m), and

*D*

_{s}is a clearance constant (dimensionless). The spring constants can be calculated as [4] follows:

*D*

_{i}=

*Cd*

_{i}is the coil diameter and

*N*

_{a}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

*i*th spring can be expressed as follows:

The parameters are selected as *c*_{1} = 10, *c*_{2} = 15, *c*_{3} = 20, *m*_{1} = *m*_{2} = *m*_{3} = 5, *G* = 30 (for chrome silicon), *C* = 8, *F*_{u} = 6, $\kappa 1=0.05$, $\kappa 2=0.1$, *D*_{s} = 0.5, and *N*_{a} = 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 *u*_{i}(*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.

#### 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 *m*_{i} = 5 and damping coefficients *c*_{i} = 10. The objective function of the plant design part of each subsystem is assumed to be (*y*_{i} − 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.

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.

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.

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.

## 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

*t*_{f}=final time

*y*_{i}=local physical plant design variable vector in the

*i*th subsystem- $y\xafs$=
vector of combined shared physical plant design variables

*y*_{s}_{i,j}=shared physical plant design variable vector in

*i*th and*j*th subsystems*w*_{C},*w*_{P}=weighting coefficient associated with control and physical plant design objective function

*z*_{C},*z*_{P}=control and physical plant design objective function

- ( · )
_{i}= symbol ( · ) with subscript

*i*: variables/parameters in the*i*th subsystem- ( · )
_{k}= symbol ( · ) with subscript

*k*: variables/parameters at grid point*k*- $\beta $=
step size constants

- $\gamma ,\eta $=
Lagrange multipliers of physical plant design constraints

- $\u03f5D,\u03f5Z,\u03f5R,\u03f5ys$=
preset tolerance value of convergence for optimizing variables, objective function, consistency constraint residual, and shared plant design variables

- $\lambda $=
co-state variables

*μ*=dual variable vectors associated with physical plant design consistency constraints

- $\nu $=
dual variable vectors in bilevel approach

- $\rho $=
Lagrange multipliers associated with converted complicating constraints

- $\sigma $=
common value vector of shared physical plant design variables

- $\varphi $=
converted complicating constraints at collocation points

- $\psi $=
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.

*t*

_{k}≤

*t*≤

*t*

_{k+1},

*t*

_{k}≤

*t*≤

*t*

_{k+1}, the variables at the collocation point are denoted with the subscript (

*k*+ 1/2) and are evaluated as [17] (

*k*= 0, …,

*M*− 1),

*t*

_{0}= 0 and

*t*

_{M}=

*t*

_{f}. The constraints are formulated as follows:

*k*in subsystem

*i*is defined as follows:

Note that $\varphi i,k$ is related to the variable values of *SSj* at the grid and the collocation points, i.e., *x*_{j,k}, *x*_{j,k+1/2}, and *x*_{j,k+1}.