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, yi, which exist only in the ith subsystem; and shared physical variables, $ysi,j$, which are shared between the ith and jth subsystems. The state and control variables of the ith subsystem are represented by xi(t) and ui(t). When shared physical variables appear in more than two subsystems, for example, among ith, jth, and kth 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 ith and jth subsystems is determined by the coupling matrices Aij in the dynamic equations. The set V is defined such that if Aij is a nonzero matrix, then the integer pair (i, j) ∈ V; otherwise, if Aij = [0], then $(i,j)∉V$.

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.

For a system consisting of N subsystems, the codesign problem for the ith subsystem is formulated as
$minui(t),yi,ysi,jZi=wPizPi(yi,ysi,j)+12wCi∫0τ(xiT(t)Qixi(t)+uiT(t)Riui(t)+2xiT(t)Siui(t))dtsubject to the constraints fi(yi,ysi,j)≤0 (physical inequality constraints)gi(yi,ysi,j)=0 (physical equality constraints)ẋi(t)=Ai(yi,ysi,j)xi(t)+Bi(yi,ysi,j)ui(t) +∑(i,j)∈VAij(yj,ysi,j)xj(t) (system dynamics)xi(0)=xi0 (initial conditions)$
(1)

The weights associated with the physical and control objective functions of the ith subsystem are $wPi$ and $wCi$, where $wPi∈[0,1], wCi∈[0,1]$, and $wPi+wCi=1$ for all i ∈ {1,…, N}. In the control objective function, Qi and Ri are real symmetric positive definite weight matrices. The matrices Ai, Bi, and Aij representing the system dynamics can be linearly or nonlinearly parameterized in terms of the physical variables.

The overall codesign problem is formulated as
$minui(t),yi,ysi,jZ=∑i=1NwPizPi(yi,ysi,j)+12∑i=1NwCi∫0τ(xiT(t)Qixi(t)+uiT(t)Riui(t)+2xiT(t)Siui(t))dtsubject to the constraintsfi(yi,ysi,j)≤0, gi(yi,ysi,j)=0ysi,j=ysj,i, ∀(i,j)∈Vẋi(t)=Ai(yi,ysi,j)xi(t)+Bi(yi,ysi,j)ui(t) +∑(i,j)∈VAij(yj,ysi,j)xj(t)xi(0)=xi0(i,j=1,…,N,j≠i)$
(2)

In Sec. 3, a multilevel decentralized approach is developed to obtain optimal solution to the N-subsystem problem described by Eq. (2).

## 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 (SSi, 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 Lagrangian and the Hamiltonian can be written as
$L=∑i=1N(wPizPi(yi,ysi,j)+ηiTfi(yi,ysi,j)+γiTgi(yi,ysi,j))+μT(ysi,j−ysj,i)$
(3)

$H=∑i=1N[12(xiT(t)Qixi(t)+uiT(t)Riui(t)+2xiT(t)Siui(t)) +piT(t)(Ai(y)xi(t)+Biui(t)+∑(i,j)∈VAij(y)xj(t))]$
(4)
The necessary optimality conditions for the overall codesign problem (Eq. (2)) are
$fi(yi,ysi,j)≤0, gi(yi,ysi,j)=0$
(5)

$ηi≥0, ηi°fi(yi,ysi,j)=0$
(6)

$ysi,j=ysj,i, ∀(i,j)∈V$
(7)

$∂H∂ui=0, ∂H∂xi=−ṗi(t), ∂H∂pi=−ẋi(t)$
(8)

$xi(0)=xi0$
(9)

$∂L∂yi+∂H∂yi=0$
(10)

$∂L∂ysi,j+∂H∂ysi,j=0(i,j=1,…,N,j≠i)$
(11)

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

To devise an iterative approach between the physical and control parts, the necessary optimality conditions are considered with respect to the physical variables that involve the Hamiltonian from the control part, namely, Eqs. (10) and (11).

### 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 nc shared physical variables. A vector $σ∈Rnc$ is introduced to represent the common values of the shared physical variables. Let $yŝ$ 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̂=Eσ$, where E is a matrix representing the set of consistency constraints for shared physical variables, with the matrix entries Eij = 1 when the physical variables are shared between the subsystems SSi and SSj, and Eij = 0 otherwise.

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 Ai, Bi, Aij, Qi, Ri, and Si, and initial conditions xi0 for all i = 1, 2,…, N, the solution steps are detailed as in the following:

1. (S1)

Select the initial values of yi0 and $ysi,j0$ for the physical variables of the ith subsystem. Set the current physical variables' values, yic and $ysi,j−c$, equal to yi0 and $ysi,j0$ for all i.

2. (S2)
Insert yic and $ysi,j−c$ in the dynamic equations of the problem (Eq. (2)) and find Ai, Aij, and Bi matrices. For the ith subsystem (i = 1,…, N), solve the independent matrix Riccati equation and obtain Pi(t)
$P˙i(t)=−Pi(t)(Ai−BiRi−1SiT)−(AiT−SiRi−1BiT)Pi(t) +Pi(t)BiRi−1BiTPi(t)−Qi+SiRi−1SiT$
(12)
with Pi(τ) = [0]. Next, select the initial values of $αi(0)(t)$, and compute $vi(0)(t)$ as
$vi(0)(t)≡∑(i,j)∈VAijxj0, ∀ t∈[0,τ]$
(13)
where the superscript (⋅)(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 lth inner loop.
3. (S3)
For the lth control inner loop, use $αi(l−1)(t)$ and $vi(l−1)(t)$ to solve $hi(l)(t)$ and $xi(l)(t)$ with the boundary condition $hi(l)(τ)=0$ and the initial condition $xi(l)(0)=xi0$ as
$h˙i(l)(t)=−(AiT−Pi(t)BiRi−1BiT−SiRi−1BiT)hi(l)(t) −Pi(t)vi(l−1)(t)+∑(i,j)∈VAjiαj(l−1)(t)$
(14)

$x˙i(l)(t)=(Ai−BiRi−1BiTPi(t)−BiRi−1SiT)xi(l)(t) −BiRi−1BiThi(l)(t)+vi(l−1)(t)$
(15)
Compute the costate variables $λi(l)(t)$
$λi(l)(t)=Pi(t)xi(l)(t)+hi(l)(t)$
(16)
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.
4. (S4)
Compute the overall interaction error $eC(l)$
$eC(l)=∑i=1N∫0τ(‖vi(l−1)(t)−∑(i,j)∈VAijxj(l)(t)‖22)dt$
(17)
If $eC(l)>ε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 $αi(l)(t)$ and $vi(l)(t)$
$αi(l)(t)=−λi(l)(t)$
(18)

$vi(l)(t)=∑(i,j)∈VAijxj(l)(t)$
(19)
and repeat from the step (S3). Otherwise, continue to the step (S5). Let xi(t), λi(t), and ui(t) denote the time-variant values of the state, costate and control variables in the last inner loop for the ith subsystem, where
$ui(t)=−Ri−1(BiTPi(t)+SiT)xi(t)−Ri−1BiThi(t)$
(20)
5. (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 yic and $ysi,j−c$ values
$∂H∂yi|yi−c,ysi,j−c=∑i=1N[wCi∫0τλiT(t)∂∂yi(Ai(yi,ysi,j)xi(t)+Bi(yi,ysi,j)ui(t)+∑(i,j)∈VAij(yi,ysi,j)xj(t))dt]$
(21)

$∂H∂ysi,j|yi−c,ysi,j−c=∑i=1N[wCi∫0τλiT(t)∂∂ysi,j(Ai(yi,ysi,j)xi(t)+Bi(yi,ysi,j)ui(t)+∑(i,j)∈VAij(yi,ysi,j)xj(t))dt]$
(22)
6. (S6)
Set the initial dual variable vector μ(0) = 0. For the ith subsystem, the Lagrangian of the physical part is
$Li=wPizPi(yi,ysi,j)+ηiTfi(yi,ysi,j)+γiTgi(yi,ysi,j)+μTysi,j$
(23)
Solve the physical variables $yi,ysi,j$ by the necessary optimality conditions in the ith subsystem
$∂Li∂yi+∂H∂yi|yi−c,ysi,j−c=0∂Li∂ysi,j+∂H∂ysi,j|yi−c,ysi,j−c=0∂Li∂γi=0fi(yi,ysi,j)≤0ηi≥0, ηi°fi(yi,ysi,j)=0$
(24)
7. (S7)
For the lth plant inner loop, let $yŝ$ be the vector of all $ysi,j$ obtained in all subsystems, compute the average of the shared physical variables by
$σ=(ETE)−1ETyŝ$
(25)
and evaluate the consistency constraint's residual
$eR(l)=‖yŝ−Eσ‖22$
(26)
If $eR(l)>εR$, then update the dual variable vector via a subgradient method with a fixed step size β
$μ(l)=μ(l−1)+βyŝ$
(27)
and repeat from Eq. (23). Otherwise, results in the plant's inner loop are obtained as yi–new and $ysi,j−new$, and continue to step (S8).
8. (S8)

If $max{‖yi−new−yi−c‖,‖ysi,j−new−ysi,j−c‖}>εP$ for i = 1,…, N, update $yi−c, ysi,j−c$ with the new physical design values $yi−new, ysi,j−new$, 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.

The numerical example is modified from an example in Ref. [26], which is originally an optimal control problem. In order to make it as a codesign problem, five physical design variables (m1 to m5) 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:
$minm1,…,m5,u1(t),u2(t)Z=wP((m1−1)2+(m2−2)2+(m3−1)2 +(m4−2)2+(m5−3)2)+12wC∫01(xT(t)Qx(t) +uT(t)Ru(t)+2xT(t)Su(t))dtsubject to the constraintsm12+m22+m32+m42−8≤0m3+m4+2m5−8=0x.(t)=Ax(t)+Bu(t)x0=[−1, 0.1, 1,−0.5]T$
(28)
where
$wP=wC=0.5A=[20.1m130.1m400.2m1−m20.1−0.50.50.15m50.50−0.2m3−0.25−m52], B=[100.1000.500.25]Q=diag(2, 1, 1, 2), R=diag(1,2)S=[11000001]T$
In order to solve the problem in a decentralized manner, the problem is decomposed into two subsystems (SS1 and SS2). The weights for the objective functions in each subsystem are $wP1=wP2=wC1=wC2=0.5$. Local physical variables in the two subsystems are y1 = [m1, m2]T and y2 = m5. The shared physical variables between the subsystems are $ys1,2=ys2,1=[m3,m4]T$. The subsystems are then formulated as
$SS1miny1,ys1,2,u1(t)Z1=wP1((m1−1)2+(m2−2)2+12(m3−1)2+12(m4−2)2)+12wC1∫01(x1T(t)Q1x1(t)+R1u12(t)+2x1T(t)S1u1(t))dtsubject to the constraintsf1(y1,ys1,2):=m12+m22+m32+m42−8≤0ẋ1(t)=A1x1(t)+B1u1(t)+A12x2(t)x10=[−1,0.1]TA1=[20.1m130.2m1−m2], A12=[0.1m400.1−0.5]B1=[1, 0.1]T, Q1=diag(2,1), R1=1, S1=[1, 1]T$
(29)

$SS2miny2,ys2,1,u2(t)Z2=wP2(12(m3−1)2+12(m4−2)2+(m5−3)2) +12wC2∫01(x2T(t)Q2x2(t)+R2u22(t)+2x2T(t)S2u2(t))dtsubject to the constraintsg2(y2,ys2,1):=m3+m4+2m5−8=0ẋ2(t)=A2x2(t)+B2u2(t)+A21x1(t)x20=[1,−0.5]TA2=[m50.5−0.25−m52], A21=[0.50.150−0.2m3]B2=[0.5, 0.25]T, Q2=diag(1,2), R2=2, S2=[0, 1]T$
(30)

The tolerance values εC = 0.01, εR = 0.01, and εP = 0.001 are, respectively, used to check whether the control inner loop, plant inner loop, and outer loop have converged.

The problem is solved below following the steps of the proposed approach.

#### Approach: Solution Steps

1. (S1)

Select the initial values of the physical variables y10 = [1,1]T, y20 = 1, and $ys1,2-0=[1,1]T$ as the initial values for the physical variables. So, $y1−c=[1,1]T, y2−c=1$, and $ys1,2−c=[1,1]T$.

2. (S2)
Obtain A1 and A2 matrices using the initial values of physical design variables
$A1=[20.10.2−1], A2=[10.5−0.25−1]A12=[0.100.1−0.5], A21=[0.50.150−0.2]$
(31)
Then, solve Eq. (12) for i = 1, 2 with boundary conditions P1(1) = P2(1) = [0] to obtain P1(t) and P2(t). Next, set initial values of the coordinator variables of α(0)(t) and v(0)(t) for t ∈ [0, 1]
$α1(0)(t)≡[0.5,0.5]Tα2(0)(t)≡[0.75,0.75]Tv1(0)(t)≡A12x20≡[0.1,0.35]Tv2(0)(t)≡A21x10≡[−0.45,−0.02]T$
(32)
3. (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=[−1,0.1]T, x20=[1,−0.5]T$ by Eq. (15).

4. (S4)
Check interaction error $eC(1)$ in Eq. (17)
$eC(1)=∫01(‖v1(0)(t)−A12x2(1)(t)‖22+‖v2(0)(t)−A21x1(1)(t)‖22)dt=0.64$
(33)
Since $eC(1)>εC$, update the coordinator variables to $αi(1)(t)$ and $vi(1)(t)$ by Eqs. (18) and (19)
$α1(1)(t)=−λ1(1)(t)=−P1(t)x1(1)(t)−h1(1)(t)α2(1)(t)=−λ2(1)(t)=−P2(t)x2(1)(t)−h2(1)(t)v1(1)(t)=A12x2(1)(t)v2(1)(t)=A21x1(1)(t)$
(34)
Repeat steps (S3) to (S4) until $eC(l)<εC$ after the lth iteration of the control's inner loop.
5. (S5)
After $eC(l)$ becomes sufficiently small, using the values of λi(t)s and xi(t)s, in the last inner loop, compute the gradient of the Hamiltonian of the control part with respect to the physical variables
$∂H∂y1|y1−c,ys1,2−c=0.5∫01λ1T(t)∂∂y1(A1x1(t))dt=(−0.074−0.048)$
(35)

$∂H∂y2|y2−c,ys1,2−c=0.5∫01λ2T(t)∂∂y2(A2x2(t))dt=0.40$
(36)

$∂H∂ys1,2|y1−c,ys1,2−c=∑i=120.5∫01λiT(t)∂∂ys1,2(Aixi(t) +∑(i,j)∈VAijxj(t))dt=(0.0035−0.047)$
(37)
6. (S6)
The Lagrangians of the physical parts in SS1 and SS2 are
${L1=wP1zP1(y1,ys1,2)+η1f1(y1,ys1,2)+μTys1,2L2=wP2zP2(y2,ys2,1)+γ2g2(y2,ys2,1)−μTys2,1$
(38)
Set the initial value of the dual variable vector μ = 0. Solve for y1 and $ys1,2$, and y2 and $ys2,1$ separately in the two subsystems using the necessary optimality conditions, as in Eq. (24)
$SS1:∂L1∂y1+∂H∂y1|y1−c,ys1,2−c=0∂L1∂ys1,2+∂H∂ys1,2|y1−c,ys1,2−c=0f1(y1,ys1,2)≤0, η1≥0, η1f1(y1,ys1,2)=0$
(39)

$SS2:∂L2∂y2+∂H∂y2|y2−c,ys1,2−c=0, ∂L2∂ys2,1+∂H∂ys1,2=0∂L2∂γ2=0, g2=0$
(40)
which yield
$y1=[0.97, 1.85]T, ys1,2=[0.82, 1.72]T$
(41)

$y2=2.53, ys2,1=[0.92, 2.02]T$
(42)
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
$yŝ=[0.82, 1.72, 0.92, 2.02]T$
(43)
7. (S7)
Since the consistency constraint is $ys1,2=ys2,1$, the matrix E is therefore
$E=[10100101]T$
Now, calculate σ by Eq. (25) and evaluate the consistency constraint residual by Eq. (26)
$σ=(ETE)−1ETyŝ=[0.87, 1.87]T$
(44)

$eR(1)=0.22$
(45)
Since $eR(1)>εR, μ$ is updated by Eq. (27) with β = 0.05
$μ=[−0.0026,−0.0075, 0.0026, 0.0075]T$
(46)
and repeat from step (S6) until the residual is small enough. The new values for the physical variables are
$y1−new=[0.92, 1.76]Ty2−new=2.65ys1,2−new=[0.87, 1.82]T$
(47)
8. (S8)

Since , update $y1−c, y2−c, ys1,2−c$ with the values of $y1−new, y2−new, ys1,2−new$, 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.

Figure 4 shows the trajectories of the optimal controller u2(t) of SS2. As shown, u2(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 ith subsystem, the state variable xi(t) = [xi1(t), xi2(t)]T denotes position and velocity of the ith mass at time t, and the control variable ui(t) is the input applied to the mass. The control objective function of each subsystem is a quadratic cost function, with the weighting matrices Qi = diag(1, 1) and Ri = 1 for 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 εC = 0.1, εR = 0.01, and εP = 0.001.

In the physical design part, it is assumed that the damping coefficients and the masses are known as $c1=10 kg·s−1, c2=15 kg·s−1, c3=20 kg·s−1$, and m1 = m2 = m3 = 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., yi = di (i = 1, 2, 3). The physical objective function in the ith 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:
$min zDi=(10.24G)(π12.8Fu)2C−2yi4subject to the constraintsC−1−κ1C−1yi−1≤1(inside diameter)κ2yi−1≤1(lower bound on wire diameter)0.8FuDsGC3yi−2≤1(clash allowance)$
(48)
where G is shear modulus (G = 30 MPa, for chrome silicon), Fu is maximum allowable force (N), C is spring index (dimensionless), κ1 is the minimum allowable inside diameter (m), κ2 is the lower bound on physical variables (m), and Ds is a clearance constant (dimensionless). The spring constants can be calculated as [4]
$ki=di4G8Di3Na(1+12C2)$
(49)

where Di = Cdi is coil diameter, and Na is the number of active coils of the spring.

In this example, the parameters are set as C = 8, Fu = 6, κ1 = 0.05, κ2 = 0.1, Ds = 0.25, and Na = 100. The problem is formulated as the following:
$minyi,ui(t) Z=∑i=13wPi(10.24G)(π12.8Fu)2C−2yi4 +12∑i=13wCi∫05(xiT(t)Qixi(t)+Riui2(t))dtsubject to the constraints(C−1−1)yi−κ1C−1≤0, i=1, 2, 3κ2−yi≤0, i=1, 2, 30.8FuDsGC3yi−2≤1, i=1, 2, 3ẋ1(t)=[01−k1+k2m1−c1+c2m1]x1(t)+[00k2m1c2m1]x2(t) +[01m1]u1(t)ẋ2(t)=[01−k2+k3m2−c2+c3m2]x2(t)+[00k2m2c2m2]x1(t) +[00k3m2c3m2]x3(t)+[01m2]u2(t)ẋ3(t)=[01−k3m3−c3m3]x3(t)+[00k3m3c3m3]x2(t) +[01m3]u3(t)xi0=[1,1]T, i=1, 2, 3$
(50)

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 ui(t) and velocities of all masses xi2(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.

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

•
• eC =

interaction error in the control's inner loop

•
• eR =

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

•
• nc =

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

•
• wC =

weighting coefficient associated with the control objective function

•
• wP =

weighting coefficient associated with the physical design objective function

•
• x =

state variables

•
• x0 =

initial conditions of state variables

•
• yi =

local physical design variables in ith subsystem

•
• $ysi,j$ =

shared physical design variables in ith and jth subsystems

•
• $ŷs$ =

combined shared physical design variables

•
• Z =

system-level objective function

•
• zC =

control objective function

•
• zP =

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 ith subsystem

## References

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
,”
American Control Conference
(
ACC
), Arlington, VA, June 25–27, pp.
1864
1869
.
2.
Allison
,
J. T.
, and
Herber
,
D. R.
,
2014
, “
Multidisciplinary Design Optimization of Dynamic Engineering Systems
,”
AIAA J.
,
52
(
4
), pp.
691
710
.
3.
Allison
,
J. T.
, and
Nazari
,
S.
,
2010
, “
Combined Plant and Controller Design Using Decomposition-Based Design Optimization and the Minimum Principle
,”
ASME
Paper No. DETC2010-28887.
4.
Allison
,
J. T.
, and
Han
,
Z.
,
2011
, “
Co-Design of an Active Suspension Using Simultaneous Dynamic Optimization
,”
ASME
Paper No. DETC2011-48521.
5.
Allison
,
J. T.
,
Guo
,
T.
, and
Han
,
Z.
,
2015
, “
Co-Design of an Active Suspension Using Simultaneous Dynamic Optimization
,”
ASME J. Mech. Des.
,
136
(
8
), p.
081003
.
6.
Allison
,
J. T.
,
Herber
,
D. R.
, and
Deshmukh
,
A. P.
,
2015
, “
Integrated Design of Dynamic Sustainable Energy Systems
,”
20th International Conference on Engineering Design
(
ICED
) Milan, Italy, July 27–30, pp. 299–308.
7.
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
.
8.
Peters
,
D. L.
,
Papalambros
,
P. Y.
, and
Ulsoy
,
A. G.
,
2013
, “
Sequential Co-Design of an Artifact and Its Controller Via Control Proxy Functions
,”
Mechatronics
,
23
(
4
), pp.
409
418
.
9.
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.
Ravichandran
,
T.
,
Wang
,
D.
, and
Heppler
,
G.
,
2006
, “
Simultaneous Plant-Controller Design Optimization of a Two-Link Planar Manipulator
,”
Mechatronics
,
16
(
3
), pp.
233
242
.
11.
Sandoval
,
L. R.
,
Budman
,
H. M.
, and
Douglas
,
P. L.
,
2008
, “
Simultaneous Design and Control of Processes Under Uncertainty: A Robust Modelling Approach
,”
J. Process Control
,
18
(
7
), pp.
735
752
.
12.
Alyaqout
,
S. F.
,
Papalambros
,
P. Y.
, and
Ulsoy
,
A. G.
,
2011
, “
Combined Robust Design and Robust Control of an Electric DC Motor
,”
IEEE/ASME Trans. Mechatronics
,
16
(
3
), pp.
574
582
.
13.
Reyer
,
J. A.
, and
Papalambros
,
P. Y.
, 2002, “Combined Optimal Design and Control With Application to an Electric DC Motor,”
ASME. J. Mech. Des.
, 124(2), pp. 183–191.
14.
Deshmukh
,
A. P.
, and
Allison
,
J. T.
,
2013
, “
Design of Nonlinear Dynamic Systems Using Surrogate Models of Derivative Functions
,”
ASME
Paper No. DETC2013-12262.
15.
Zhang
,
W.
,
Li
,
Q.
, and
Guo
,
L.
,
1999
, “
Integrated Design of Mechanical Structure and Control Algorithm for a Programmable Four-Bar Linkage
,”
IEEE/ASME Trans. Mechatronics
,
4
(
4
), pp.
354
362
.
16.
Yan
,
H.-S.
, and
Yan
,
G.-J.
,
2009
, “
Integrated Control and Mechanism Design for the Variable Input-Speed Servo Four-Bar Linkages
,”
Mechatronics
,
19
(
2
), pp.
274
285
.
17.
Sridharan
,
S.
,
Echols
,
J. A.
,
Rodriguez
,
A. A.
, and
Mondal
,
K.
,
2014
, “
Integrated Design and Control of Hypersonic Vehicles
,”
American Control Conference
, (
ACC
) Portland, OR, June 4–6, pp.
1371
1376
.
18.
Kim
,
M.-J.
, and
Peng
,
H.
,
2006
, “
Combined Control/Plant Optimization of Fuel Cell Hybrid Vehicles
,”
American Control Conference
(
ACC
), Minneapolis, MN, June 14–16, pp. 1–6.
19.
Geoffrion
,
A. M.
,
1972
, “
Generalized Benders Decomposition
,”
J. Optim. Theory Appl.
,
10
(
4
), pp.
237
260
.
20.
Tosserams
,
S.
,
Etman
,
L. F. P.
, and
Rooda
,
J. E.
,
2007
, “
An Augmented Lagrangian Decomposition Method for Quasi-Separable Problems in MDO
,”
Struct. Multidiscip. Optim.
,
34
(
3
), pp.
211
227
.
21.
Tosserams
,
S.
,
Etman
,
L. F. P.
, and
Rooda
,
J. E.
,
2008
, “
,”
Int. J. Numer. Methods Eng.
,
73
(
13
), pp.
1885
1910
.
22.
Geromel
,
J. C.
, and
Bernussou
,
J.
,
1982
, “
Optimal Decentralized Control of Dynamic Systems
,”
Automatica
,
18
(
5
), pp.
545
557
.
23.
Nedic
,
A.
, and
Ozdaglar
,
A.
,
2009
, “
Distributed Subgradient Methods for Multi-Agent Optimization
,”
IEEE Trans. Autom. Control
,
54
(
1
), pp.
48
61
.
24.
Nedic
,
A.
,
Ozdaglar
,
A.
, and
Parrilo
,
P. A.
,
2010
, “
Constrained Consensus and Optimization in Multi-Agent Networks
,”
IEEE Trans. Autom. Control
,
55
(
4
), pp.
922
938
.
25.
Maasoumy
,
M.
,
Pinto
,
A.
, and
Sangiovanni-Vincentelli
,
A.
,
2011
, “
Model-Based Hierarchical Optimal Control Design for HVAC Systems
,”
ASME
Paper No. DSCC2011-6078.
26.
Jamshidi
,
M.
,
1996
, “
A Multi-Level Structure for a Class of Dynamical Optimization Problems
,” Master's thesis, Case Western Reserve University, Cleveland, OH.
27.
Jamshidi, M., 1997,
Large-Scale Systems: Modeling, Control, and Fuzzy Logic
,
Prentice Hall
28.
Cohen
,
G.
, and
Joalland
,
G.
,
1976
, “
Coordination Methods by the Prediction Principle in Large Dynamic Constrained Optimization Problems
,”
IFAC Proc. Vol.
,
9
(
3
), pp.
539
547
.
29.
Smith
,
N. J.
, and
Sage
,
A. P.
,
1973
, “
An Introduction to Hierarchical Systems Theory
,”
Comput. Electr. Eng.
,
1
(
1
), pp.
55
71
.
30.
,
N.
, and
,
A.
,
2006
, “
Optimal Control of Robot Manipulators With a New Two-Level Gradient-Based Approach
,”
Electr. Eng.
,
88
(
5
), pp.
383
393
.
31.
,
A. J.
, and
Bell
,
D. J.
,
1986
, “
A Simplified Algorithm for Optimization of Large-Scale Gas Networks
,”
Optim. Control Appl. Methods
,
7
(
1
), pp.
95
104
.
32.
Joalland
,
G.
, and
Cohen
,
G.
,
1980
, “
Optimal Control of a Water Distribution Network by Two Multilevel Methods
,”
Automatica
,
16
(
1
), pp.
83
88
.
33.
Boyd
,
S.
,
Xiao
,
L.
,
Mutapcic
,
A.
, and
Mattingley
,
J.
,
2007
, “
Notes on Decomposition Methods Notes for EE364B
,” Stanford University, Stanford, CA, accessed Sept. 26, 2017, hthttps://stanford.edu/class/ee364b/lectures/decomposition_notes.pdf
34.
Peters
,
D. L.
,
Papalambros
,
P. Y.
, and
Ulsoy
,
A. G.
,
2009
, “
On Measures of Coupling Between the Artifact and Controller Optimal Design Problems
,”
ASME
Paper No. DETC2009-86868.
35.
Rutquist
,
P. E.
, and
Edvall
,
M. M.
,
2010
, “
Propt-MATLAB Optimal Control Software
,” TOMLAB Optimization, Pullman, WA.
36.
Azarm
,
S.
, and
Papalambros
,
P.
,
1982
, “
An Interactive Design Procedure for Optimization of Helical Compression Springs
,” University of Michigan, Ann Arbor, MI, Technical Report No.
UM-MEAM-82-7
.