On Decentralized Optimization for a Class of Multisubsystem Codesign Problems

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.

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 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 x_i(t) and u_i(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 A_ij in the dynamic equations. The set V is defined such that if A_ij is a nonzero matrix, then the integer pair (i, j) ∈ V; otherwise, if A_ij = [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.

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, Q_i and R_i are real symmetric positive definite weight matrices. The matrices A_i, B_i, and A_ij representing the system dynamics can be linearly or nonlinearly parameterized in terms of the physical variables.

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

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

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.

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

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

The details of the solution steps for the proposed approach are given in Sec. 3.3.

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_i, and S_i, and initial conditions x_i0 for all i = 1, 2,…, N, the solution steps are detailed as in the following:

1. (S1)

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

2. (S2)
   Insert y_i–c and $ysi,j−c$ in the dynamic equations of the problem (Eq. (2)) and find A_i, A_ij, and B_i matrices. For the ith subsystem (i = 1,…, N), solve the independent matrix Riccati equation and obtain P_i(t)

   $P˙i(t)=−Pi(t)(Ai−BiRi−1SiT)−(AiT−SiRi−1BiT)Pi(t) +Pi(t)BiRi−1BiTPi(t)−Qi+SiRi−1SiT$

   (12)

   with P_i(τ) = [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 (⋅)⁽⁰⁾ 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 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 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 y_i–c 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. 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 $yŝ$ be the vector of all $ysi,j$ obtained in all subsystems, compute the average of the shared physical variables by

   $σ=(ETE)−1ETyŝ$

   (25)

   and evaluate the consistency constraint's residual

   $eR(l)=‖yŝ−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)+βyŝ$

   (27)

   and repeat from Eq. (23). Otherwise, results in the plant's inner loop are obtained as y_i–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.

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.

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.

1. (S1)

   Select the initial values of the physical variables y₁₀ = [1,1]^T, y₂₀ = 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 A₁ and A₂ 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 P₁(1) = P₂(1) = [0] to obtain P₁(t) and P₂(t). Next, set initial values of the coordinator variables of α⁽⁰⁾(t) and v⁽⁰⁾(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 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

   $∂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 y₁ and $ys1,2$⁠, and y₂ 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

   $yŝ=[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)−1ETyŝ=[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 $max{‖y1−new−y1−c‖, ‖y2−new−y2−c‖,‖ys1,2−new−ys1,2−c‖}>εP$⁠, 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.