Codesign refers to the process of integrating the optimization of the physical plant design and control of a system. In this paper, a new class of codesign problems with a multisubsystem architecture in both design and control is formulated and solved. Our work here extends earlier research on models and solution approaches from single system to multisubsystem codesign. In this class, the optimization model for the physical design part in each subsystem is assumed to have a convex objective function with convex inequality and linear equality constraints. The optimization model for the control part of each subsystem belongs to a class of finite time-horizon linear quadratic regulator (LQR) feedback control. A new multilevel decentralized method is proposed that can obtain optimal or near-optimal solutions for this class of codesign problems. Details of the model and approach are presented and demonstrated by a numerical as well as a more complex spring–mass–damper system example. The proposed decentralized approach has been compared with a centralized approach. Using a scalable test problem, it is shown that as the size of the problem is increased, the computation effort for the decentralized approach increases linearly while that of the centralized approach increases nonlinearly.

## Introduction

Codesign refers to an integrated optimization of the physical plant and controller for an engineering system. One of the key challenges for solving codesign problems is that the optimization is simultaneously applied to both time-invariant variables (physical design variables) and time-variant variables (state and control variables), which are usually coupled via design constraints and dynamic equations.

Existing papers in codesign mainly focus on single-system problems, namely, physical, and control design subproblems are considered as parts of a single system. The traditional sequential method is to optimize the physical system first followed by that for the control system. However, such a method does not fully handle the coupling between physical and control design variables in the codesign problems, so it can lead to suboptimal solutions [1]. As a result, codesign methods have been developed and applied to engineering examples in order to optimize all variables via a single-system codesign structure [2]. For example, Allison and Nazari [3] developed a decomposition-based approach to solve single-system codesign problems. Allison et al. studied codesign of the active suspension systems [4,5] and dynamic sustainable energy systems [6] based on a direct transcription method. Peters et al. [7,8] derived control proxy functions for certain types of codesign problems so that the optimal solution of the modified sequential approach matches with the simultaneous solution. Jiang et al. [9] proposed an iterative method to solve codesign problems with nonlinear control systems. Ravichandran et al. [10] adopted a metaheuristic evolutionary algorithm to obtain the optimal physical design values and set-point control of a two-link planar manipulator for carrying different payloads. Sandoval et al. [11] investigated optimization of the parameters and input of the dynamics of a chemical process with uncertainty. Optimal codesign of many other engineering applications are reported as well, including direct current motors [12,13], wind turbines [14], four-bar mechanism [15,16], and vehicles [17,18]. In many of the reported papers, the codesign problems considered do not contain a physical design objective function when evaluating the system's performance. Furthermore, no literature was found that specifically proposed and formulated optimal codesign problems with a multisubsystem architecture in both design and control parts.

On the other hand, multisubsystem (or decentralized) optimization models and methods have been developed and utilized to solve several (physical) engineering design problems. Examples of such methods include Benders' decomposition [19] and augmented Lagrangian decomposition [20,21]. The basic idea behind these methods is to decompose the entire problem into multiple subproblems, optimize the subproblems, and coordinate the solutions among the subproblems. Decentralized approaches have also been applied to control systems. For example, Geromel and Bernussou [22] proposed a feasible direction method using a matrix gradient to find the decentralized optimal controllers. Nedic et al. [23,24] developed a decentralized subgradient method for multi-agent system optimization. Maasoumy et al. [25] derived and applied a hierarchical control method to a multiroom heating, ventilation, and air conditioning system model. One approach that has been extended in this paper is the interactive prediction method [26,27]. Cohen and Joalland [28] provided proofs of convergence of this method in the application of linear quadratic regulator (LQR) feedback control. Smith and Sage [29] developed the interaction prediction method to solve nonlinear control problems. Applications of the interactive prediction method include optimal control of robot manipulators [30] and the water distribution network [31,32]. In short, the literature on decentralized optimization has focused on the decomposition-based optimization models and methods for either the physical or the control design problems.

The main objective of this paper is to develop formulation and decentralized solution approach for a class of multisubsystem codesign problems. The proposed approach extends and integrates two previous hierarchical decomposition techniques, namely, the interaction prediction method [26] and the dual decomposition method [33]. The interaction prediction method is an indirect method, which efficiently solves the class of control problems considered. However, there can be limitations to extending the indirect method to more general problems. The dual decomposition method is used to solve the design part with convex inequality constraints. In the proposed approach, the necessary optimality conditions for the multisubsystem codesign problem are used, as motivated by the work of Fathy et al. [1], to establish the connection between designs of the physical and control parts of the subsystems, leading to system-level optimal or near-optimal solutions. Additionally, a weighted sum of the physical and control design objective functions is formulated, subject to the constraints, so as to evaluate the performance of each subsystem. Hence, multi-objective codesign problems are considered in the proposed work.

Compared with centralized solution methods (e.g., methods that solve all variables as in a single system), the proposed multilevel decentralized approach has two key advantages. The first advantage is that the optimization of one subsystem may not require information from all other subsystems. This property helps when optimization of each subsystem is preferred to be relatively independent. The second advantage of the decentralized approach is that the size of the individual subproblems to solve is smaller compared to that in a centralized approach. In particular, as shown in this paper by way of a scalable test problem, as the size of the problem is increased, the computational effort for the proposed decentralized approach increases linearly.

The rest of this paper is organized as follows: In Sec. 2, a description for a class of multisubsystem codesign problems, together with formulations and assumptions, is provided. In Sec. 3, a multilevel decentralized approach is proposed to solve the described problems, where necessary optimality conditions and dual decomposition are presented, and the solution steps are described including mathematical details of the hierarchical scheme. In Sec. 4, the results from two examples solved by the proposed decentralized approach are presented. The first numerical example demonstrates the approach step by step. The second example shows a model for a series of three spring–mass–damper subsystems. Optimal results obtained from the centralized and decentralized methods are compared and contrasted for both examples. Finally, a comparison of computational times for solving a scalable test problem using the decentralized and centralized approaches is presented. In Sec. 5, a summary for the paper is provided.

## Multisubsystem Codesign Problem Formulation

Consider a multisubsystem codesign problem with *N* coupled subsystems. As an example, Fig. 1 shows the structure of such a problem with three subsystems. The physical and control design variables are optimized in the problem. Two types of physical design variables are considered: local physical variables, *y _{i}*, which exist only in the

*i*th subsystem; and shared physical variables, $ysi,j$, which are shared between the

*i*th and

*j*th subsystems. The state and control variables of the

*i*th subsystem are represented by

*x*(

_{i}*t*) and

*u*(

_{i}*t*). When shared physical variables appear in more than two subsystems, for example, among

*i*th,

*j*th, and

*k*th subsystems, consistency constraints $ysi,j=ysj,k$, $ysi,k=ysj,k$ and $ysi,j=ysi,k$ are introduced.

The interaction between the state variables of the *i*th and *j*th subsystems is determined by the coupling matrices *A _{ij}* in the dynamic equations. The set

*V*is defined such that if

*A*is a nonzero matrix, then the integer pair (

_{ij}*i*,

*j*) ∈

*V*; otherwise, if

*A*= [0], then $(i,j)\u2209V$.

_{ij}The problem is formulated as a bi-objective optimization problem for each subsystem, which consists of one physical and one control design subproblem, and the objective function is the weighted sum of the two-part objective functions. The overall objective function of the problem is the summation of individual objective functions of all *N* subsystems.

One of the main challenges in the codesign problem is the coupling between the physical, state, and control design variables, as well as the physical design constraints and system dynamic equations. Another key challenge in the multisubsystem case, the coupling between state variables and physical variables of different subsystems, is also considered.

A solution to the codesign problems considered in this paper is developed under the following assumptions:

- (A1)
The overall physical and control design objective functions are additively separable into the individual subsystem physical and control objective functions. This assumption simplifies the decomposition of the problem.

- (A2)
The coupling is unidirectional. By unidirectional coupling, it is meant that the physical design objective functions and constraints are dependent only on physical variables, and the control objective functions and constraints are dependent on physical and control variables [34].

- (A3)
Physical design objective functions are convex, physical design inequality constraints are convex, and physical equality constraints are linear functions of physical variables. In addition, active inequality and equality constraints are assumed to be linearly independent. This assumption is made to guarantee uniqueness of the solution in each step of the approach.

- (A4)
The control part of each subsystem belongs to a class of continuous finite time-horizon LQR feedback control, where the time-horizon is within [0,

*τ*]. Cross terms between the state and control variables of the same subsystem are considered in the control objective function, while the cross terms between the state and control variables of different subsystems are not considered, for the sake of simplicity of the decomposition in the solution steps. - (A5)
The subsystems are controllable for all feasible values of the physical variables, which ensures the feasibility and existence of an optimal solution in the codesign problem.

*N*subsystems, the codesign problem for the

*i*th subsystem is formulated as

The weights associated with the physical and control objective functions of the *i*th subsystem are $wPi$ and $wCi$, where $wPi\u2208[0,1],\u2009wCi\u2208[0,1]$, and $wPi+wCi=1$ for all *i* ∈ {1,…, *N*}. In the control objective function, *Q _{i}* and

*R*are real symmetric positive definite weight matrices. The matrices

_{i}*A*,

_{i}*B*and

_{i},*A*representing the system dynamics can be linearly or nonlinearly parameterized in terms of the physical variables.

_{ij}## Multilevel Decentralized Approach

As an example, Fig. 2 shows an overview of the proposed multilevel decentralized approach for solving a three-subsystem codesign problem. The approach consists of the design of the physical and the control parts. Decomposition-based optimization techniques are applied in both parts.

First, the idea of interaction prediction method [26] is used in the control part (see Fig. 2, the “Control” box), as described next. The bottom level of the control part solves for the individual optimal controller of each subsystem (SS*i*, *i* = 1, 2, 3). At the top level, the overall interaction error is calculated by the state and coordinator variables from interconnected subsystems. The coordinator variables are used to estimate the coupling between the subsystems. The coordinator variables are updated when the error is large, namely, the estimation of the coupling is not accurate at the bottom level. When the error becomes sufficiently small, the decentralized solution of the control part is obtained.

For the physical part (see Fig. 2, the “Plant” box), the physical variables are optimized using a dual decomposition scheme via a projected subgradient method [33]. At the bottom level, the partial Lagrangian is formulated with an initial dual variable vector. The local physical variables and copies of shared physical variables are optimized in each subsystem. At the top level, a consistency constraint residual is checked, and the dual variables are updated if the residual is not small enough.

In order to obtain the optimal or near-optimal solution of the entire codesign problem, the connection between the physical and the control parts needs to be established. The two parts are connected via the gradients of the Hamiltonian of the control part with respect to the physical variables. These values are computed from the necessary optimality conditions for the entire codesign problem. Then the physical and the control parts are solved iteratively to obtain a solution of the codesign problem.

If the problem is to find an optimal codesign for only one system, the proposed approach becomes the same as the “nested” strategy in Ref. [1]. Since a multisubsystem codesign problem is considered in this paper, hierarchical optimization techniques in both physical and control parts are applied, and so the approach can be regarded as an extension of the nested strategy [1].

### Necessary Optimality Conditions.

In this section, the necessary optimality conditions for the proposed multisubsystem codesign are derived using the Karush–Kuhn–Tucker conditions and Pontryagin maximum principle. These conditions are used to solve for variable values later on in the solution steps of the proposed approach. The idea is motivated by the previous work from Fathy et al. [1], where necessary optimality conditions for single-system codesign problems were developed.

The notation “∘” in Eq. (6) refers to the Hadamard product, meaning an elementwise multiplication of two vectors or two matrices.

### Dual Decomposition.

A dual decomposition technique [33] is used to optimize the physical design variables of the plant (see Fig. 2). Since both the local and shared physical variables as well as the coupling constraints exist in the defined class of codesign problems considered here, the dual decomposition via projected subgradient method is used to solve the plant part in a decentralized way [33]. Each subsystem has the local physical variables, copies of the shared physical variables and the associated dual variables, which are locally optimized at the bottom level. At the top level, the dual variable vector is updated iteratively to ensure consistency of the copies of shared physical variables.

Suppose there are a total number of *n _{c}* shared physical variables. A vector $\sigma \u2208Rnc$ is introduced to represent the common values of the shared physical variables. Let $ys\u0302$ be the vector that includes all shared physical variables $ysi,j$ from all subsystems for all

*i*,

*j*. Then the consistency constraints can be written as $ys\u0302=E\sigma $, where

*E*is a matrix representing the set of consistency constraints for shared physical variables, with the matrix entries

*E*= 1 when the physical variables are shared between the subsystems SS

_{ij}*i*and SS

*j*, and

*E*= 0 otherwise.

_{ij}The details of the solution steps for the proposed approach are given in Sec. 3.3.

### Approach: Solution Steps.

The flowchart describing the approach is shown in Fig. 3. Step 1 (or S1) is to choose the initial values of the physical variables. The steps (S2)–(S4) are based on the interaction prediction method [26], which solves the control part in a decentralized way. The iterations to run the interaction prediction method are denoted as the “Control Inner Loop,” see Fig. 3. After the converged results are obtained in the control part, the gradients of the Hamiltonian with respect to the physical variables are computed, in the step (S5), which are the link between the plant and control parts. Steps (S6) and (S7) are denoted as the “Plant Inner Loop,” in which the plant part is optimized through the dual decomposition [33] via a projected subgradient method. The gradients of the Hamiltonian with respect to the physical variables are used at the bottom level to reveal the influence from the control part. At the top level, if the consistency constraint residual is not small enough, the dual variables are updated and the iteration is repeated from (S6). Otherwise, the process goes to the step (S8), checking the difference between the converged physical design variable values obtained in the step (S7) and the values in the step (S2). If the difference is small enough, the approach ends and the converged results are obtained. The iterations in the steps (S2)–(S8) are denoted as the “Outer Loop” of the proposed approach, see Fig. 3.

In the multisubsystem optimal codesign problem, given a set of fixed $wPi$'s and $wCi$'s, the parameter values in the matrices *A _{i}*,

*B*,

_{i}*A*,

_{ij}*Q*,

_{i}*R*and

_{i},*S*, and initial conditions

_{i}*x*

_{i}_{0}for all

*i*= 1, 2,…,

*N*, the solution steps are detailed as in the following:

- (S1)
Select the initial values of

*y*_{i}_{0}and $ysi,j0$ for the physical variables of the*i*th subsystem. Set the current physical variables' values,*y*_{i}_{–}and $ysi,j\u2212c$, equal to_{c}*y*_{i}_{0}and $ysi,j0$ for all*i*. - (S2)Insert
*y*_{i}_{–}and $ysi,j\u2212c$ in the dynamic equations of the problem (Eq. (2)) and find_{c}*A*,_{i}*A*, and_{ij}*B*matrices. For the_{i}*i*th subsystem (*i*= 1,…,*N*), solve the independent matrix Riccati equation and obtain*P*(_{i}*t*)with$P\u02d9i(t)=\u2212Pi(t)(Ai\u2212BiRi\u22121SiT)\u2212(AiT\u2212SiRi\u22121BiT)Pi(t)\u2003+Pi(t)BiRi\u22121BiTPi(t)\u2212Qi+SiRi\u22121SiT$(12)*P*(_{i}*τ*) = [0]. Next, select the initial values of $\alpha i(0)(t)$, and compute $vi(0)(t)$ aswhere the superscript (⋅)$vi(0)(t)\u2261\u2211(i,j)\u2208VAijxj0,\u2003\u2200\u2009t\u2208[0,\tau ]$(13)^{(0)}in this step refers to the initial value of a variable, and in the following steps, the superscript (⋅)^{(}^{l}^{)}refers to the value of a variable in the*l*th inner loop. - (S3)For the
*l*th control inner loop, use $\alpha i(l\u22121)(t)$ and $vi(l\u22121)(t)$ to solve $hi(l)(t)$ and $xi(l)(t)$ with the boundary condition $hi(l)(\tau )=0$ and the initial condition $xi(l)(0)=xi0$ as$h\u02d9i(l)(t)=\u2212(AiT\u2212Pi(t)BiRi\u22121BiT\u2212SiRi\u22121BiT)hi(l)(t)\u2003\u2212Pi(t)vi(l\u22121)(t)+\u2211(i,j)\u2208VAji\alpha j(l\u22121)(t)$(14)Compute the costate variables $\lambda i(l)(t)$$x\u02d9i(l)(t)=(Ai\u2212BiRi\u22121BiTPi(t)\u2212BiRi\u22121SiT)xi(l)(t)\u2003\u2212BiRi\u22121BiThi(l)(t)+vi(l\u22121)(t)$(15)Note that in (S2) and (S3), Eqs. (12)–(16) are calculated within each subsystem, so it is possible to compute these values in a parallel manner.$\lambda i(l)(t)=Pi(t)xi(l)(t)+hi(l)(t)$(16) - (S4)Compute the overall interaction error $eC(l)$If $eC(l)>\epsilon C$, namely, the error of the estimation of the interaction between subsystems in the step (S3) is not sufficiently small, then update the coordinator variables to $\alpha i(l)(t)$ and $vi(l)(t)$$eC(l)=\u2211i=1N\u222b0\tau (\Vert vi(l\u22121)(t)\u2212\u2211(i,j)\u2208VAijxj(l)(t)\Vert 22)dt$(17)$\alpha i(l)(t)=\u2212\lambda i(l)(t)$(18)and repeat from the step (S3). Otherwise, continue to the step (S5). Let$vi(l)(t)=\u2211(i,j)\u2208VAijxj(l)(t)$(19)
*x*(_{i}*t*),*λ*(_{i}*t*), and*u*(_{i}*t*) denote the time-variant values of the state, costate and control variables in the last inner loop for the*i*th subsystem, where$ui(t)=\u2212Ri\u22121(BiTPi(t)+SiT)xi(t)\u2212Ri\u22121BiThi(t)$(20) - (S5)Starting with this step, the approach moves to the physical part. Calculate the gradient of the Hamiltonian of the control part with respect to the physical variables at
*y*_{i}_{–}and $ysi,j\u2212c$ values_{c}$\u2202H\u2202yi|yi\u2212c,ysi,j\u2212c=\u2211i=1N[wCi\u222b0\tau \lambda iT(t)\u2202\u2202yi(Ai(yi,ysi,j)xi(t)+Bi(yi,ysi,j)ui(t)+\u2211(i,j)\u2208VAij(yi,ysi,j)xj(t))dt]$(21)$\u2202H\u2202ysi,j|yi\u2212c,ysi,j\u2212c=\u2211i=1N[wCi\u222b0\tau \lambda iT(t)\u2202\u2202ysi,j(Ai(yi,ysi,j)xi(t)+Bi(yi,ysi,j)ui(t)+\u2211(i,j)\u2208VAij(yi,ysi,j)xj(t))dt]$(22) - (S6)Set the initial dual variable vector
*μ*^{(0) }= 0. For the*i*th subsystem, the Lagrangian of the physical part isSolve the physical variables $yi,ysi,j$ by the necessary optimality conditions in the$Li=wPizPi(yi,ysi,j)+\eta iTfi(yi,ysi,j)+\gamma iTgi(yi,ysi,j)+\mu Tysi,j$(23)*i*th subsystem$\u2202Li\u2202yi+\u2202H\u2202yi|yi\u2212c,ysi,j\u2212c=0\u2202Li\u2202ysi,j+\u2202H\u2202ysi,j|yi\u2212c,ysi,j\u2212c=0\u2202Li\u2202\gamma i=0fi(yi,ysi,j)\u22640\eta i\u22650,\u2003\eta i\xb0fi(yi,ysi,j)=0$(24) - (S7)For the
*l*th plant inner loop, let $ys\u0302$ be the vector of all $ysi,j$ obtained in all subsystems, compute the average of the shared physical variables byand evaluate the consistency constraint's residual$\sigma =(ETE)\u22121ETys\u0302$(25)If $eR(l)>\epsilon R$, then update the dual variable vector via a subgradient method with a fixed step size$eR(l)=\Vert ys\u0302\u2212E\sigma \Vert 22$(26)*β*and repeat from Eq. (23). Otherwise, results in the plant's inner loop are obtained as$\mu (l)=\mu (l\u22121)+\beta ys\u0302$(27)*y*_{i}_{–new}and $ysi,j\u2212new$, and continue to step (S8). - (S8)
If $max{\Vert yi\u2212new\u2212yi\u2212c\Vert ,\Vert ysi,j\u2212new\u2212ysi,j\u2212c\Vert}>\epsilon P$ for

*i*= 1,…,*N*, update $yi\u2212c,\u2009ysi,j\u2212c$ with the new physical design values $yi\u2212new,\u2009ysi,j\u2212new$, and repeat from the step (S2). Otherwise, the approach has converged and the solution is obtained.

In Sec. 4, two examples are solved using the proposed approach.

## Examples

Two examples are given in this section. The first one is a numerical example, which is used to demonstrate the proposed approach of Sec. 3.3 step by step. The second example is an engineering model of a serial spring–mass–damper subsystems, in which each mass is separately controlled while the wire diameters of the springs are the physical design variables. Finally, a scalable test problem is implemented by increasing the number of the spring–mass–damper subsystems.

### Numerical Example.

*m*

_{1}to

*m*

_{5}) are introduced in the example. The physical design objective functions and constraints are also added following the assumptions in this paper. The overall problem is formulated as:

*y*

_{1}= [

*m*

_{1},

*m*

_{2}]

^{T}and

*y*

_{2}=

*m*

_{5}. The shared physical variables between the subsystems are $ys1,2=ys2,1=[m3,m4]T$. The subsystems are then formulated as

The tolerance values *ε _{C}* = 0.01,

*ε*= 0.01, and

_{R}*ε*= 0.001 are, respectively, used to check whether the control inner loop, plant inner loop, and outer loop have converged.

_{P}The problem is solved below following the steps of the proposed approach.

#### Approach: Solution Steps

- (S1)
Select the initial values of the physical variables

*y*_{10}= [1,1]^{T},*y*_{20}= 1, and $ys1,2-0=[1,1]T$ as the initial values for the physical variables. So, $y1\u2212c=[1,1]T,\u2009y2\u2212c=1$, and $ys1,2\u2212c=[1,1]T$. - (S2)Obtain
*A*_{1}and*A*_{2}matrices using the initial values of physical design variablesThen, solve Eq. (12) for$A1=[20.10.2\u22121],\u2003A2=[10.5\u22120.25\u22121]A12=[0.100.1\u22120.5],\u2003A21=[0.50.150\u22120.2]$(31)*i*= 1, 2 with boundary conditions*P*_{1}(1) =*P*_{2}(1) = [0] to obtain*P*_{1}(*t*) and*P*_{2}(*t*). Next, set initial values of the coordinator variables of*α*^{(0)}(*t*) and*v*^{(0)}(*t*) for*t*∈ [0, 1]$\alpha 1(0)(t)\u2261[0.5,0.5]T\alpha 2(0)(t)\u2261[0.75,0.75]Tv1(0)(t)\u2261A12x20\u2261[0.1,0.35]Tv2(0)(t)\u2261A21x10\u2261[\u22120.45,\u22120.02]T$(32) - (S3)
Compute $h1(1)(t)$ and $h2(1)(t)$ with boundary conditions $h1(1)(1)=h2(1)(1)=0$ by Eq. (14), and $x1(1)(t)$ and $x2(1)(t)$ with initial conditions $x10=[\u22121,0.1]T,\u2009x20=[1,\u22120.5]T$ by Eq. (15).

- (S4)Check interaction error $eC(1)$ in Eq. (17)Since $eC(1)>\epsilon C$, update the coordinator variables to $\alpha i(1)(t)$ and $vi(1)(t)$ by Eqs. (18) and (19)$eC(1)=\u222b01(\Vert v1(0)(t)\u2212A12x2(1)(t)\Vert 22+\Vert v2(0)(t)\u2212A21x1(1)(t)\Vert 22)dt=0.64$(33)Repeat steps (S3) to (S4) until $eC(l)<\epsilon C$ after the$\alpha 1(1)(t)=\u2212\lambda 1(1)(t)=\u2212P1(t)x1(1)(t)\u2212h1(1)(t)\alpha 2(1)(t)=\u2212\lambda 2(1)(t)=\u2212P2(t)x2(1)(t)\u2212h2(1)(t)v1(1)(t)=A12x2(1)(t)v2(1)(t)=A21x1(1)(t)$(34)
*l*th iteration of the control's inner loop. - (S5)After $eC(l)$ becomes sufficiently small, using the values of
*λ*(_{i}*t*)s and*x*(_{i}*t*)s, in the last inner loop, compute the gradient of the Hamiltonian of the control part with respect to the physical variables$\u2202H\u2202y1|y1\u2212c,ys1,2\u2212c=0.5\u222b01\lambda 1T(t)\u2202\u2202y1(A1x1(t))dt=(\u22120.074\u22120.048)$(35)$\u2202H\u2202y2|y2\u2212c,ys1,2\u2212c=0.5\u222b01\lambda 2T(t)\u2202\u2202y2(A2x2(t))dt=0.40$(36)$\u2202H\u2202ys1,2|y1\u2212c,ys1,2\u2212c=\u2211i=120.5\u222b01\lambda iT(t)\u2202\u2202ys1,2(Aixi(t)\u2003+\u2211(i,j)\u2208VAijxj(t))dt=(0.0035\u22120.047)$(37) - (S6)The Lagrangians of the physical parts in SS1 and SS2 areSet the initial value of the dual variable vector${L1=wP1zP1(y1,ys1,2)+\eta 1f1(y1,ys1,2)+\mu Tys1,2L2=wP2zP2(y2,ys2,1)+\gamma 2g2(y2,ys2,1)\u2212\mu Tys2,1$(38)
*μ*= 0. Solve for*y*_{1}and $ys1,2$, and*y*_{2}and $ys2,1$ separately in the two subsystems using the necessary optimality conditions, as in Eq. (24)$SS1:\u2202L1\u2202y1+\u2202H\u2202y1|y1\u2212c,ys1,2\u2212c=0\u2202L1\u2202ys1,2+\u2202H\u2202ys1,2|y1\u2212c,ys1,2\u2212c=0f1(y1,ys1,2)\u22640,\u2003\eta 1\u22650,\u2003\eta 1f1(y1,ys1,2)=0$(39)which yield$SS2:\u2202L2\u2202y2+\u2202H\u2202y2|y2\u2212c,ys1,2\u2212c=0,\u2003\u2202L2\u2202ys2,1+\u2202H\u2202ys1,2=0\u2202L2\u2202\gamma 2=0,\u2003g2=0$(40)$y1=[0.97,\u20091.85]T,\u2003ys1,2=[0.82,\u20091.72]T$(41)By the assumption A3, in this step, the group of equations (Eqs. (39) and (40)) have a unique solution. The combined shared physical variable vector is$y2=2.53,\u2003ys2,1=[0.92,\u20092.02]T$(42)$ys\u0302=[0.82,\u20091.72,\u20090.92,\u20092.02]T$(43) - (S7)Since the consistency constraint is $ys1,2=ys2,1$, the matrix
*E*is thereforeNow, calculate$E=[10100101]T$*σ*by Eq. (25) and evaluate the consistency constraint residual by Eq. (26)$\sigma =(ETE)\u22121ETys\u0302=[0.87,\u20091.87]T$(44)Since $eR(1)>\epsilon R,\u2009\mu $ is updated by Eq. (27) with$eR(1)=0.22$(45)*β*= 0.05and repeat from step (S6) until the residual is small enough. The new values for the physical variables are$\mu =[\u22120.0026,\u22120.0075,\u20090.0026,\u20090.0075]T$(46)$y1\u2212new=[0.92,\u20091.76]Ty2\u2212new=2.65ys1,2\u2212new=[0.87,\u20091.82]T$(47) - (S8)
Since $max{\Vert y1\u2212new\u2212y1\u2212c\Vert ,\Vert y2\u2212new\u2212y2\u2212c\Vert ,\Vert ys1,2\u2212new\u2212ys1,2\u2212c\Vert}\epsilon P$, update $y1\u2212c,\u2009y2\u2212c,\u2009ys1,2\u2212c$ with the values of $y1\u2212new,\u2009y2\u2212new,\u2009ys1,2\u2212new$, and repeat previous steps from step (S2). Stop when the maximal value of all the norms is less than

*ε*._{P}

#### Optimization Results and Discussion.

The final optimal solutions obtained by different approaches are shown in Table 1. The centralized solution is obtained by TOMLAB™ [35], which solves for all variables simultaneously. The results show that the optimal solution obtained by the decentralized method matches well with the centralized solution, when using the same initial conditions. For this example, the decentralized approach takes less than or equal to four iterations for all control inner loops, around 15 iterations for all plant inner loops, and a total of nine iterations for the outer loop to converge. As a comparison and as shown in Table 1, if the codesign problem is solved in a sequential manner, it obtains a solution, which is worse than the result for the centralized and decentralized approaches.

$y1\u2217$ | $y2\u2217$ | $ys1,2\u2217$ | Z^{*} | |
---|---|---|---|---|

Decentralized | [0.68, 1.29]^{T} | 2.31 | [1.33, 2.03]^{T} | 3.38 |

Centralized | [0.7064, 1.3217]^{T} | 2.3375 | [1.3263, 1.9987]^{T} | 3.3353 |

Sequential | [0.9313, 1.8626]^{T} | 2.7292 | [0.8052, 1.7365]^{T} | 4.1479 |

$y1\u2217$ | $y2\u2217$ | $ys1,2\u2217$ | Z^{*} | |
---|---|---|---|---|

Decentralized | [0.68, 1.29]^{T} | 2.31 | [1.33, 2.03]^{T} | 3.38 |

Centralized | [0.7064, 1.3217]^{T} | 2.3375 | [1.3263, 1.9987]^{T} | 3.3353 |

Sequential | [0.9313, 1.8626]^{T} | 2.7292 | [0.8052, 1.7365]^{T} | 4.1479 |

Figure 4 shows the trajectories of the optimal controller *u*_{2}(*t*) of SS2. As shown, *u*_{2}(*t*) obtained by the centralized and decentralized methods are nearly the same while the magnitudes of both are smaller than that obtained from the sequential approach, leading to a better optimal value of objective function.

Note that different values of initial points for physical design variables are tested in this example, including both feasible and infeasible initial points. The test results indicate that the converged solutions for the physical variables are nearly identical, which means that the proposed approach is robust to an initial point.

### Engineering Example: Spring–Mass–Damper System Model.

This example is composed of three spring–mass–damper subsystems, which are connected to each other in series, as shown in Fig. 5.

In the control part, for the *i*th subsystem, the state variable *x _{i}*(

*t*) = [

*x*

_{i}_{1}(

*t*),

*x*

_{i}_{2}(

*t*)]

^{T}denotes position and velocity of the

*i*th mass at time

*t*, and the control variable

*u*(

_{i}*t*) is the input applied to the mass. The control objective function of each subsystem is a quadratic cost function, with the weighting matrices

*Q*= diag(1, 1) and

_{i}*R*= 1 for

_{i}*i*= 1, 2, 3. The system-level objective function is a weighted sum of all subsystems' physical and control design objective functions, with the associated weights chosen as $wPi=wCi=0.5$ for all

*i*. The final time is:

*τ*= 5 s. The tolerance values are set as

*ε*= 0.1,

_{C}*ε*= 0.01, and

_{R}*ε*= 0.001.

_{P}*m*

_{1}=

*m*

_{2}=

*m*

_{3}= 5 kg. The physical design part of the model considered in this example is based on a spring design optimization formulation in Ref. [36]. The wire diameters of springs are selected as the physical variables, i.e.,

*y*=

_{i}*d*(

_{i}*i*= 1, 2, 3). The physical objective function in the

*i*th subsystem, $zDi$, is to maximize energy storage capacity. It is shown in Ref. [36] 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 (

*G*= 30 MPa, for chrome silicon),

*F*is maximum allowable force (N),

_{u}*C*is spring index (dimensionless),

*κ*

_{1}is the minimum allowable inside diameter (m),

*κ*

_{2}is the lower bound on physical variables (m), and

*D*is a clearance constant (dimensionless). The spring constants can be calculated as [4]

_{s}where *D _{i}* =

*Cd*is coil diameter, and

_{i}*N*is the number of active coils of the spring.

_{a}*C*= 8,

*F*= 6,

_{u}*κ*

_{1}= 0.05,

*κ*

_{2}= 0.1,

*D*= 0.25, and

_{s}*N*= 100. The problem is formulated as the following:

_{a}The proposed decentralized approach is applied to solve this problem. Note that all physical variables are shared variables since they exist in more than one subsystem. Thus, the copies of the shared physical variables are optimized in each subsystem.

#### Optimization Results and Discussion.

The problem is solved by the proposed decentralized approach and by the centralized approach using TOMLAB™ [35]. The solutions for the optimal damping coefficients are shown in Table 2. The profiles of optimal controllers *u _{i}*(

*t*) and velocities of all masses

*x*

_{i}_{2}(

*t*) obtained by both centralized and decentralized methods are shown in Figs. 6 and 7. As is observed, the state variable profiles of the decentralized approach are very close to those of the centralized approach.

$y1\u2217$ | $y2\u2217$ | $y3\u2217$ | Z^{*} | |
---|---|---|---|---|

Decentralized | 0.10 | 0.10 | 0.10 | 7.34 |

Centralized | 0.1166 | 0.1000 | 0.1000 | 6.8558 |

Sequential | 0.1000 | 0.1000 | 0.1000 | 6.9860 |

$y1\u2217$ | $y2\u2217$ | $y3\u2217$ | Z^{*} | |
---|---|---|---|---|

Decentralized | 0.10 | 0.10 | 0.10 | 7.34 |

Centralized | 0.1166 | 0.1000 | 0.1000 | 6.8558 |

Sequential | 0.1000 | 0.1000 | 0.1000 | 6.9860 |

From Table 2, it can be observed that the optimal value of the entire objective function of the decentralized approach is slightly higher than that of the centralized and the sequential methods. However, the limitation of a centralized or a sequential approach is that the optimization requires full information from all subsystems to solve all variables simultaneously. While in a decentralized manner, the full information is not needed, and smaller subproblems can be solved in the process. Take SS1 as an example; in the decentralized approach, SS1 can be optimized with only considering its coupling with SS2, and this subproblem has a smaller size than the overall problem.

#### Computational Cost.

In this section, it is shown that, by way of a scalable test problem, as the size of the codesign problem (e.g., number of subsystems) increases, the proposed decentralized approach performs better in terms of computational cost than the centralized approach. The test problem is a spring–mass–damper with *N* subsystems (with each subsystem composed of a simple spring–mass–damper subproblem), as shown in Fig. 8. With the same objective functions and constraints as in the engineering example, in each subproblem, there exist two state variables (displacement and velocity of the mass), one control variable (force on the mass), and one physical plant variable (spring's wire diameter). The number of subsystems considered for each test problem is: 5, 10, 20, 30, all the way to 100. Each of these 11 test problems is solved by the centralized approach as well as the proposed decentralized approach using the same tolerance value for convergence. As an example, when the number of subsystems is equal to 100, there are 200 state, 100 control, and 100 physical plant variables.

Figure 9 shows a comparison of the computational (CPU) time between the centralized and decentralized approach for the solution of this test problem as the number of subsystems is increased. The computational time of the centralized approach is based on the CPU time of the solutions obtained from TOMLAB™ [35], which employs an efficient approach for solving codesign problems by an all-in-one approach. Figure 9 clearly shows that as the number of variables increases, the computational time of the decentralized approach increases about linearly with respect to the number of subsystems. On the other hand, the computational cost of the centralized approach grows nonlinearly.

## Summary

This paper presents the formulation for a class of multisubsystem codesign problems. In this class of problems, the physical design part in each subsystem has a convex objective function with convex inequality and linear equality constraints. The control part of each subsystem belongs to a class of finite time-horizon LQR feedback control. While many papers on codesign optimization problems have been published in recent years, codesign problems with a multisubsystem structure that has both physical design and control domains in the same subsystem have not been specifically considered. The main challenge for this class of problems is that the optimal solution of the physical and control design variables within each subsystem and also across different subsystems can be dependent on each other.

The proposed decentralized approach is composed of two hierarchical decomposition techniques, one for the physical part and the other for the control part of the codesign subsystems. The codesign problems are solved in a decentralized manner via a hierarchical structure inside each part. The necessary optimality conditions are derived for the overall problem, as well as for the subsystem subproblems. Derived from the necessary optimality conditions, the gradients of the Hamiltonian of the control part with respect to the physical variables bridge the physical and control parts, providing the quantitative information of the coupling between the two parts.

The decentralized approach is applied to several examples. In the first numerical example, the approach is illustrated in a detailed step-by-step manner. In the second example, an engineering codesign model of a three spring–mass–damper system is solved. Finally, a scalable test problem is considered. For this problem, the test results show that the computational time of the proposed decentralized approach increases approximately linearly with respect to an increase in the number of subsystems while the computational cost of the centralized approach grows nonlinearly.

## Acknowledgment

We would also like to thank all reviewers for their constructive and helpful comments and suggestions.

## Funding Data

Naval Air Warfare Center, Aircraft Division (Agreement No. N004211720018).

Office of Naval Research (Grant No. N000141310160).

## Nomenclature

*A*,*B*=matrices in dynamic equations of control part

*E*=consistency constraint matrix

*e*=_{C}interaction error in the control's inner loop

*e*=_{R}consistency constraint residual in the plant's inner loop

*f*=physical inequality constraints

*g*=physical equality constraints

*h*=ancillary variables in the control's inner loop

*H*=Hamiltonian of the control's subproblems

*L*=Lagrangian of the physical design subproblems

*N*=number of subsystems

*n*=_{c}number of common shared physical variables

*p*=costate variables

*Q*,*R*,*S*=weighting matrices in objective function in linear quadratic regulator

*u*=control variables

*v*=coordinator variable in the control's inner loop

*V*=set of pairs of connected subsystems

*w*=_{C}weighting coefficient associated with the control objective function

*w*=_{P}weighting coefficient associated with the physical design objective function

*x*=state variables

*x*_{0}=initial conditions of state variables

*y*=_{i}local physical design variables in

*i*th subsystem- $ysi,j$ =
shared physical design variables in

*i*th and*j*th subsystems- $y\u0302s$ =
combined shared physical design variables

*Z*=system-level objective function

*z*=_{C}control objective function

*z*=_{P}design objective function

*α*=coordinator variables in the control's inner loop

*β*=step size in the plant's inner loop

*γ*,*η*=Lagrange multipliers of the physical design constraints

*ε*=_{C}preset tolerance value of convergence for state variables

*ε*=_{P}preset tolerance value of convergence for physical design variables

*ε*=_{R}preset tolerance value of convergence for consistency constraint residual

*λ*=ancillary costate variables

*μ*=dual variables for the physical design consistency constraints

*σ*=common value vector of shared physical variables

*τ*=final time

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

*i*: variables/parameters in the*i*th subsystem