Game-Theoretic Mixed $H_2/H_{infty}$ Control with Sparsity Constraint for Multi-agent Networked Control Systems

Game-Theoretic Mix ed H 2 /H ∞ Control with Sparsity Constraint for Multi-agent Network ed Control Systems Feier Lian, Aranya Chakrabortty Senior Member , IEEE , and Alexandra Duel-Hallen F ellow , IEEE Abstract Multi-agent networked control systems (NCSs) are often subject to model uncertainty and are limited by large communication cost, associated with feedback of data between the system nodes. T o provide robustness against model uncertainty and to reduce the communication cost, this paper in vestigates the mixed H 2 /H ∞ control problem for NCS under the sparsity constraint. First, proximal alternating linearized minimization (P ALM) is employed to solve the centralized social optimization where the agents hav e the same optimization objecti ve. Next, we in vestig ate a sparsity-constrained noncooperative game, which accommodates different control-performance criteria of different agents, and propose a best-response dynamics algorithm based on P ALM that con ver ges to an approximate Generalized Nash Equilibrium (GNE) of this game. A special case of this game, where the agents hav e the same H 2 objectiv e, produces a partially-distributed social optimization solution. W e validate the proposed algorithms using a network with unstable node dynamics and demonstrate the superiority of the proposed P ALM-based method to a previously in vestigated sparsity-constrained mixed H 2 /H ∞ controller . Index T erms sparse controller , H 2 and H ∞ control, Linear Matrix Inequality (LMI), model uncertainty , nonconv ex nonsmooth optimization, game theory The authors are with the Department of Electrical and Computer Engineering, North Carolina State University , Raleigh, NC 27695 USA (e-mail: flian2@ncsu.edu; aranya.chakrabortty@ncsu.edu; sasha@ncsu.edu) Financial support from the NSF Grant EECS 1544871 is gratefully acknowledged. 1 Game-Theoretic Mix ed H 2 /H ∞ Control with Sparsity Constraint for Multi-agent Network ed Control Systems I . I N T RO D U C T I O N Recent research on multi-agent networks has proposed various methods for reducing communication cost for control by using sparse H 2 control designs [1]–[6]. Howe ver , most of this work ignores the effects of model uncertainties, which are bound to arise in most practical large-scale systems since the network operating conditions and topology change frequently ov er time. Even if the topology is fixed, the exact model parameters are not always av ailable to the designer . The sparse control design in such cases must also be robust against the uncertainties. Robust designs have been reported in sev eral recent papers such as [7]–[10] using H ∞ control, which is suitable for handling norm-bounded uncertainties in the system dynamics. In particular , [8], [9] employ both H 2 and H ∞ control, thus balancing the H 2 performance of the nominal system and robustness objectiv e. Howe ver , optimizing the H 2 performance under the H ∞ and sparsity constraints has not been in vestigated by other researchers. Moreover , while a global control performance cost was optimized, this metric did not address the individual objecti ves of multiple agents under uncertainty . Control of NCS under uncertainties has received signifi- cant attention recently in various domains, such as wide-area control of power systems [2], multi-robot coordination, multi- access broadcast channel, vehicle formation, and wireless sensor networks [11], [12]. In these problems, game theory becomes a po werful tool, with dif ferent control inputs modeled as game players, where each player aims to optimize its individual objective using an associated control policy . Differ- ential games hav e been in vestigated for uncertain multi-agent systems, and algorithms for finding an equilibrium point have been proposed based on solving a set of coupled optimization problems. The works [13], [14] extend Nash-type differential game in [15] by finding robust Nash strategies while either considering polytopic uncertainty or formulating uncertain external disturbance as a fictitious player . The works [16]–[19] model the uncertainty using stochastic differential equations, and the Nash strategies are found by solving cross-coupled matrix equations through necessary optimality conditions or Karush-Kuhn-T ucker (KKT) conditions. Recently , reinforce- ment learning has been applied to multi-agent control Nash games when the system parameters are completely or partially unknown [20]–[23]. Howe ver , these reported game-theoretic designs did not address any sparsity constraint. In this paper , we in vestig ate controller designs for multi- agent systems that aim to reduce H 2 cost under H ∞ and sparsity constraints. Both social optimization and a noncon vex game, where the H 2 -objectiv es of the agents are same and dif- ferent, respectiv ely , are in vestigated. W e model our uncertainty as a norm-bounded parameteric uncertainty that translates into an H ∞ constraint as in [9], [24]. W e employ the proximal alternating linearized minimization (P ALM) [5], which has been shown to be effecti ve for optimization for noncon ve x nonsmooth problems [25] and was utilized in [5] in a sparsity- constrained output-feedback co-design problem. First, a cen- tralized sparsity-constrained mixed H 2 /H ∞ controller , which represents the social optimization, is addressed. Preliminary results on this topic were recently reported in our conference paper [24], where we dev eloped a centralized controller under the sparsity and H ∞ constraints using a greedy gradient sup- port pursuit (GraSP) method [26]. Howe ver , the algorithm in [24] requires the kno wledge of an initial stabilizing feedback gain that satisfies the sparsity constraint. W e eliminate this requirement and sho w the advantages of the P ALM-based controller in this paper ov er that in [24] in terms of the quadratic H 2 -cost. Next, we extend the proposed design to the multi-agent scenario where each agent designs its own part of the feedback matrix, subject to a shared global H ∞ -norm and sparsity constraints. Since each agent has different individual cost, the control design is modeled as a noncooperativ e game with shared constraints. W e dev elop a numerical algorithm to find the Generalized Nash Equilibrium (GNE) [27] of this game. The proposed algorithm has partially-distributed computation, i.e., in the first stage, each player computes its own feedback matrix while in the second stage, the sparse links are chosen globally based on the results from the first stage. Third, assuming all players of the game hav e the same H 2 -optimization objectiv e, we develop a potential game that yields a partially-distributed implementation of the social op- timization. W e perform numerical simulations to demonstrate the performance of the proposed algorithms and discuss their con ver gence properties. The main contributions of the paper can, therefore, be summarized as: • Development and analysis of a centralized, sparsity- constrained mix ed H 2 /H ∞ controller for social optimiza- tion of multi-agent systems with norm-bounded uncer- tainty . • Development of game-theoretic, partially-distributed al- gorithms that aim to minimize the H 2 -norms of the agents’ transfer functions under shared sparsity and H ∞ - norm constraints. 2 T ABLE I: Notation T erm Definition M  0(  0) Matrix M is positive definite (semidefinite) M ≺ 0(  0) M is negati ve definite (semidefinite) σ max ( M ) maximum singular value of M || K || F Frobenius norm of the matrix K , defined by p trace( K T K ) . card( K ) Cardinality of matrix K , defined by the number of nonzero elements in K . ∇ K J ( K ) The gradient of the scalar function J ( K ) with respect to the matrix K . Assuming K ∈ R m × n , ∇ K J ( K ) is gi ven by a m × n matrix with the elements [ ∇ K J ( K )] ij = ∂ J /∂K ij . [ K ] s The matrix obtained by preserving only the s largest- magnitude entries of the matrix K and setting all other entries to zero. H 2 norm The system H 2 norm in the time domain is || G || 2 =  R ∞ 0  trace( H ( t ) T H ( t ))  dt  1 / 2 , where H ( t ) is the impulse response of the system. H ∞ norm The system H ∞ norm || G || ∞ = sup ω σ max ( G ( j ω )) , where G ( j ω ) is the system transfer matrix. The rest of the paper is organized as follows. Section II presents the system model with parametric uncertainty and dev elops a centralized P ALM algorithm for social optimization using sparsity-constrained mixed H 2 /H ∞ control. Section III describes a multi-agent system with parametric uncertainty , proposes a noncooperative game with shared sparsity and H ∞ constraints, and de velops partially-distributed numerical algo- rithms for this game and for social optimization. Section IV demonstrates effecti veness of the proposed algorithms using numerical simulations, and discusses their complexities and con ver gence properties. Section V discusses future directions and concludes the paper . Throughout the paper , matrices are denoted with boldface capital letters. If M is a symmetric matrix, the upper block matrices are sometimes denoted by ∗ to save space. Some notation used is summarized in T able I. I I . P A L M A L G O R I T H M F O R C E N T R A L I Z E D S PAR S I T Y - C O N S T R A I N E D M I X E D H 2 /H ∞ C O N T RO L A. System model and mixed H 2 /H ∞ contr ol Consider the following linear time-in variant system with model uncertainty: ˙ x ( t ) = ( A + ∆ A ) | {z } ˆ A x ( t ) + ( B + ∆ B ) | {z } ˆ B u ( t ) + B 2 w 2 ( t ) z 2 ( t ) = C 2 x ( t ) + D 2 u ( t ) + D 22 w 2 ( t ) y ( t ) = Cx ( t ) , (1) where x ( t ) ∈ R n × 1 is the state vector , u ( t ) ∈ R m × 1 is the control input vector , w 2 ( t ) ∈ R m 2 × 1 is the exogenous input, z 2 ( t ) ∈ R p 2 × 1 is the performance output, and y ( t ) ∈ R p × 1 is the measured output. The matrices A and B are the nominal values of the state and input matrices, respectiv ely , while ∆ A and ∆ B model the respectiv e uncertainties. W e make the following assumptions: Assumption 1: (i) The pair ( A , B ) is controllable, ( C , A ) is observable, ( C 2 , A ) is observable. (ii) ∆ A and ∆ B hav e the form [9] [∆ A ∆ B ] = B 1 ∆ δ [ C 1 D 1 ] , (2) where B 1 ∈ R n × m 1 , C 1 ∈ R p 1 × n , D 1 ∈ R p 1 × m are known ma- trices, and ∆ δ ∈ R m 1 × p 1 is an unknown matrix which is norm- bounded, satisfying ∆ δ T ∆ δ  1 /γ 2 I for any scalar γ > 0 . Assumption 2: Matrices C 2 and D 2 hav e the following form: C 2 =  C 1 2 , 0  , D 2 =  0 D 2 2  , C T 2 D 2 = 0 . (3) Using assumptions 1 and 2, the system (1) can be expressed as the feedback interconnection of the following two subsys- tems: Σ :          ˙ x ( t ) = Ax ( t ) + Bu ( t ) + B 1 w 1 ( t ) + B 2 w 2 ( t ) z 1 ( t ) = C 1 x ( t ) + D 1 u ( t ) z 2 ( t ) = C 2 x ( t ) + D 2 u ( t ) + D 22 w 2 ( t ) y ( t ) = Cx ( t ) (4) Σ K : n w 1 ( t ) = ∆ δ z 1 ( t ) , (5) where z 1 ( t ) ∈ R p 1 × 1 , w 1 ( t ) ∈ R m 1 × 1 . Our goal in this section is to find a linear static output- feedback controller u ( t ) = − Ky ( t ) that stabilizes (1), i.e., guarantees || T z 1 w 1 ( K ) || ∞ < γ . Follo wing [28], this H ∞ -norm constraint can be transformed into an Linear Matrix Inequality (LMI) condition as stated in the following theorem. Theorem 1: The H ∞ -norm constraint || T z 1 w 1 ( K ) || ∞ < γ holds if and only if there exists an X = X T that satisfies the LMI  A cl ( K ) X + XA cl ( K ) T + B 1 B T 1 X C cl 1 ( K ) T C cl 1 ( K ) X − γ 2 I  ≺ 0 X  0 (6) where A cl ( K ) = A − BKC , C cli ( K ) = C i − D i K C for i = 1 , 2 . The mixed H 2 /H ∞ control problem can then be stated as: Gi ven an achiev able H ∞ -norm bound γ , find a feedback controller K ∈ R m × p that solves Minimize K || T z 2 w 2 ( K ) || 2 , (7) s.t u ( t ) = − Ky ( t )) , equation (4) holds, and (8) || T z 1 w 1 ( K ) || ∞ < γ ( or equiv alently , (6) holds ) . (9) T z 2 w 2 is the closed-loop transfer function from w 2 to z 2 . For simplicitly , and without loss of generality , we set D 22 = 0 in (1). Follo wing standard rob ust control results, such as in [28], it can be shown that the squared H 2 norm from w 2 to z 2 for the system (4) is || T z 2 w 2 ( K ) || 2 2 = J ( K ) := trace( B T 2 PB 2 ) (10) where P is the solution of the L yapunov equation P A cl ( K ) + A cl ( K ) T P + C cl 2 ( K ) T C cl 2 ( K ) = 0 . (11) W e can define Q = ( C 1 2 ) T C 1 2  0 , R = ( D 2 2 ) T D 2 2  0 (12) in which case the objectiv e J ( K ) in (10) can be written as J ( K ) = Z ∞ t =0  x T ( t ) Qx ( t ) + u T ( t ) Ru ( t )  dt. (13) B. Sparsity-constrained mixed H 2 /H ∞ contr ol The solution K in problem (7)-(9) is usually a dense matrix, meaning that ev ery sensor must send the sensed outputs to ev ery controller . This can result in a large communication cost. 3 T o reduce this cost, we impose a sparsity constraint on the feedback matrix [2], [24], resulting in the following sparsity- constrained mixed H 2 /H ∞ problem: min K || T z 2 w 2 ( K ) || 2 , s.t. || T z 1 w 1 ( K ) || ∞ < γ , card( K ) ≤ s, (14) with the plant model satisfying (4). For simplicity , we define each nonzero entry in K as one communication link. Alterna- tiv e definitions of sparsity and their effects on the actual cost of communication are discussed in [2]. C. Overview of the centralized P ALM algorithm The sparsity-constrained mixed H 2 /H ∞ problem (14) was solved in our recent paper [24] using the GraSP algorithm, assuming that for any gi ven value of s we can find an initial guess for K that satisfies the s -sparse structure. Depending on the plant model and the uncertainty in (4), finding such a feasible initial guess in reality , ho wev er , can be quite difficult. In this section, we eliminate this requirement by introducing a sparsity-constrained optimization algorithm based on P ALM. For this, we first transform (14) into a problem with two optimization v ariables, K and F , where K is defined in (8), and F represents the sparse feedback matrix that satisfies the cardinality constraint. The problem (14) can then be reformulated as follows: min K , F J ( K ) + ρ 2 || K − F || 2 F , s.t. || T z 1 w 1 ( K ) || ∞ < γ , card( F ) ≤ s, (15) where the penalty term ρ/ 2 || K − F | 2 F is used to regularize the difference between K and F . When the parameter ρ is chosen large enough, this term can be reduced sufficiently . There are two constrained v ariables in (15). W e next transform (15) to an unconstrained optimization problem by defining the follo wing indicator functions: g ( K ) =  0 , T ∞ ( K ) < γ + ∞ , O .W . (16) f ( F ) =  0 , card( F ) ≤ s + ∞ , O .W . (17) Using (16) and (17) the problem (15) can be written as min K , F Φ( K , F ) , (18) where Φ( K , F ) = J ( K ) + g ( K ) + f ( F ) + H ( K , F ) (19) with H ( K , F ) = ρ 2 || K − F || 2 F (20) being the coupling function between K and F . The P ALM algorithm proceeds by alternating the minimiza- tion on the variables ( K , F ) through separate subproblems [25], which simplifies solving (18), as described belo w . When K is fixed, the optimization (18) reduces to minimizing the sum of a nonsmooth function f ( F ) and a smooth function H ( K , F ) of F . From the result of proximal forward-backward splitting algorithm [25], minimizing f + H can be relaxed as iterativ ely upper bounding the objecti ve and minimizing the upper bound [29]. The iteration can be written as F k +1 = (21) arg min F  h F − F k , ∇ F H ( F , K ) i + t 2 || F − F k || 2 F + f ( F )  , where h , i denotes inner product. Minimizing the first two terms is equiv alent to minimizing the first order (linear) approximation of H ( K , F ) at F = F k , regularized by a trust-region penalty near F k . When t ∈ ( L, ∞ ) and L is the Lipschitz constant (see S.II [30] and Appendix B.3 [31]) for ∇ F H ( K , F ) , the regularized linear approximation provides an upper bound on H ( K , F ) [29]. Eq. (21) can be rewritten compactly using the definition of a proximal map as F k +1 ∈ prox f t  F k − 1 /t ∇ F H ( K , F k )  , (22) where for σ : R d → ( ∞ , ∞ ] , a proper and lower semicontinu- ous function, x ∈ R d and t > 0 , the proximal map associated with σ at point x is pro x σ t ( x ) = arg min u ∈ R d  σ ( u ) + t 2 || u − x || 2  . (23) Similar analysis can be carried out for the minimization of (18) when F is fixed. In summary , the P ALM algorithm minimizes (18) by alternativ ely finding the proximal maps: F k +1 ∈ pro x f a k  F k − 1 /a k ∇ F H ( K k , F k )  (24) K k +1 ∈ pro x J + g b k  K k − 1 /b k ∇ K H ( K k , F k +1 )  (25) where a k and b k are positiv e constants that are greater than the Lipschitz constants L 1 ( K k ) and L 2 ( F k +1 ) of ∇ F H ( K k , F ) and ∇ K H ( K , F k +1 ) , respecti vely . D. Algorithm description W e summarize the P ALM algorithm for sparsity-constrained mixed H 2 /H ∞ control in Algorithm 1. In Steps 2 and 3 of Algorithm 1, F -minimization (24) and K -minimization (25) are performed, respecti vely . In Step 2, we perform iterative F -minimization (24), which can be re written from (23) as: F k +1 = arg min F n f ( F ) + a k 2 || F − Z k || 2 F o (26) where Z k is the point within the parenthesis in (24), found in Step 2.2. It is easy to see that the partial gradients of H ( K , F ) , Algorithm 1 P ALM algorithm for the mixed H 2 /H ∞ control algorithm with sparsity constraint Given s : sparsity constraint, γ : H ∞ -norm bound. 1. Initialization : K 0 : any stabilizing feedback gain with T ∞ ( K 0 ) < γ . F 0 : any stabilizing feedback gain F 0 . Compute a := γ 1 ρ , b := γ 2 ρ . for k = 1 , 2 , ...k max until || K k +1 − K k || F <  1 or || F k +1 − F k || F <  2 do // 2. F -minimization step 2.1 Compute Z k := F k − 1 a ∇ F H ( K k , F k ) 2.2 Prune Z k : F k +1 := [ Z k ] s . // 3. K -minimization step 3.1 Compute X k := K k − 1 b ∇ K H ( K k , F k +1 ) . 3.2 Update K k +1 : K k +1 := K P R OX O P ( K k , X k , b ) . end f or 4 defined in (20) are ∇ K H ( K , F ) = ρ ( K − F ) ∇ F H ( K , F ) = ρ ( F − K ) . (27) From (27), the Lipschitz constant L 1 ( K k ) = ρ , and thus the constant a k in (24) and (26) is defined as a k = a = γ 1 ρ (28) with γ 1 > 1 . In Step 3 of Algorithm 1, we perform iterativ e K - minimization: K k +1 = arg min K  J ( K ) + g ( K ) + b k 2 || K − X k || 2 F  , (29) which is equiv alent to (25), and b k is chosen as b k = b = γ 2 ρ , with γ 2 > 1 . In the following, we present the solutions for Eq. (26) and (29) used in Steps 2 and 3 of Algorithm 1. 1) F -minimization: Applying the proximal operator (24) of function f is equiv alent to minimizing a regularized version of f . In (26), f is an indicator function of the set X = { F | card( K ) ≤ s } (17), so the proximal operator in (26) (Step 2.3 in Algorithm 1) is equi valent to the projection of Z k onto the set X . Therefore, we can re write (26) as F k +1 = arg min F || F − Z k || 2 F s.t. card( F ) ≤ s. (30) As sho wn in [5], [25], the solution to (30) is [ Z k ] s (see T able I), which is Step 2.3 of Algorithm 1. Algorithm 2 K P RO X O P : Subroutine to solve (31) 1: procedure K P RO X O P ( K cur , X k , b ) 2: while T rue do 3: K prev := K cur 4: if ||∇ K h ( K cur ) || F <  3 then 5: // Stationary point in the interior of the H ∞ - constraint set. 6: break 7: end if 8: // T ake a gradient-descent step in the interior of the H ∞ -constraint set. 9: K cur := K prev − d ∇ K h ( K prev ) , where step size d > 0 is chosen by backtracking line search [32] s.t. T ∞ ( K cur ) < γ 10: if d <  2 then 11: // K prev is near the boundary of the H ∞ - constraint set 12: Solve for K in using (36). Let ∆ K cur := K in − K prev . 13: K cur := K prev + d 0 ∆ K cur , where d 0 is deter- mined by backtracking line search [32]. 14: end if 15: if || K cur − K prev || F <  1 then 16: break 17: end if 18: end while 19: end procedur e 2) K -minimization: Next we focus on the proximal opera- tor for (29), which is equi valent to min K h ( K ) s.t. T ∞ ( K ) < γ (31) where h ( K ) ,  J ( K ) + b k 2 || K − X k || 2 F  . (32) W e propose to solve (31) using a feasible direction method in the search space of K , summarized in Algorithm 2. Starting from an interior point of the feasible region of the problem (In step 3.2 of Algorithm 1, K k always satisfies T ∞ ( K k ) < γ ), the algorithm first descends along the gradient of h ( K ) until the solution reaches a stationary point in the interior (line 6) or on the boundary of the constraint set. When the solution is in the interior of the feasible region, a gradient-descent update step is used (line 9). When the current solution is at the boundary and the gradient-descent direction violates the H ∞ -norm constraint, we seek a direction that reduces the minimization objective and simultaneously moves the solution away from the boundary of the H ∞ -norm constraint set to its interior (lines 12-13 in Algorithm 2). In lines 12-13, we find the improving feasible direction for (31) when the solution is at the boundary of the feasible re gion. W e recall the Zoutendijk’ s method [33] as the foundation for general feasible direction methods, which requires ev aluating the gradients of both the objectiv e and the constraint functions, i.e., ∇ K h ( K ) and ∇ K T ∞ ( K ) . In the original formulation of Zoutendijk’ s method (S.I [30], Appendix B.2 in [31]), the gradient of the constraint function is ev aluated to form conditions for the improving feasible direction. Howe v er , due to the dif ficulty in ev aluating the gradient of an H ∞ norm, we utilize the concept of lev el sets as in [34], as well as their LMI interpretation, to dev elop an alternati ve condition. In each step of the Zoutendijk’ s method, a linear program- ming subproblem is solved to find the improving feasible direction. The inequality trace[( ∇ K h ( K )) T · ∆ K ] < 0 guar- antees that an update direction ∆ K decreases h ( K ) in (31). Moreov er , the inequality which in volv es gradient of the H ∞ norm, i.e., trace[( ∇ K || T z 1 w 1 ( K ) || ∞ ) T · ∆ K ] < 0 can be used to check if ∆ K mo ves away from the H ∞ bound. In the gain space of K ∈ R m × n , the set of all stabilizing K which satisfy (6), i.e., with H ∞ norm smaller than γ , is a le vel set K ( γ ) := { K | || T z 1 w 1 ( K ) || ∞ < γ } . (33) Giv en a stabilizing gain K 0 , the algorithm in [34] proceeds by first finding a sufficiently small γ 0 such that K 0 ∈ K ( γ 0 ) . Next, a con vex subset ˆ K ( γ 0 ) of K ( γ 0 ) , which also contains K 0 near the boundary , can be formed using an LMI sufficient condition. Then an inner point K in of ˆ K ( γ 0 ) is found using the follo wing suf ficient LMI condition. For the above ( K 0 , γ 0 ), K in is an inner point of ˆ K ( γ 0 ) if the matrix function G ( K in ; K 0 ) is positive definite, i.e., G ( K in ; K 0 )  0 ⇒ K in ∈ ˆ K ( γ 0 ) . (34) The details of computing G ( K in ; K 0 ) are provided in [24]. W e combine the LMI condition (34) and the gradient of h ( K ) in (31) to form the iterati ve algorithm for solving (31). 5 The gradient of h ( K ) is ∇ K h ( K ) = 2( RK C − B T P ) LC T + b k ( K − X k ) . (35) Thus, given a current solution K cur near the boundary of the H ∞ -norm constraint set, an improving feasible point K in can be found by solving the follo wing linear matrix inequality: max z , K in z s.t. trace[( ∇ K h ( K cur )) T ( K in − K cur )] + z ≤ 0 G ( K in ; K cur ) − θ z · I  0 . (36) The parameter θ ≥ 0 is a predetermined factor that controls how far K moves away from the H ∞ -norm boundary . The value of θ determines the speed of reduction of the H ∞ norm, with a small value of θ resulting in a less aggressive shrinkage of the H ∞ norm. If the solution z ∗ in (36) is positiv e, then K in − K cur is an improving feasible direction; otherwise, an improving feasible direction cannot be found. Giv en the current solution K cur and the inner point K in solved in (36), the update rule is giv en in lines 12–13 of Algorithm 2, where d 0 ≤ 1 is the step size found by a backtracking line search using the Armijo condition [33]. I I I . S PA R SI T Y - C O N S T R A I N E D N O N C O O P E R A T I V E G A M E S F O R M U LT I - AG E N T C O N T RO L A. Multi-agent model and generalized Nash equilibrium Next, we extend the optimization in Section II to the case when the agents have dif ferent optimization objectives. T o accommodate this scenario, we consider the following multi- agent system with model uncertainty in A and B matrices. Consider a network of N agents, where agent i employs its control strategy u i ( t ) ∈ R q i × 1 , i = 1 , ..., N . Thus, (1) becomes ˙ x ( t ) = ( A +∆ A ) x ( t ) + N X i =1  B ( i ) +∆ B ( i )  u i ( t ) + B 2 w 2 ( t ) y ( t ) = Cx ( t ) u i ( t ) = − K i y ( t ) , i = 1 , ..., N , (37) where A ∈ R n × n , B ( i ) ∈ R n × q i represent the nominal values of the state and control matrix, respectiv ely , for the i -th control input. W e assume all agents know A and B ( i ) for i = 1 , ..., N , and the uncertain matrices ∆ A ∈ R n × n and ∆ B , [∆ B (1) , ∆ B (2) , ..., ∆ B ( N ) ] satisfy the norm-bounded assumption (2), where ∆ B ( i ) ∈ R q i × n . Note that B ( i ) is the column block of B in (1), with P N i =1 B ( i ) u i = Bu , and K i ∈ R q i × p is the row block of K in (8) associated with the rows corresponding to the control inputs for agent i . Thus, the multi-agent system can be expressed in the form (4 – 5), with the first equation in (4) replaced by ˙ x ( t ) = Ax ( t ) + N X i =1 B ( i ) u i ( t ) + B 1 w 1 ( t ) + B 2 w 2 ( t ) . (38) W e introduce the following notation. Let K − i denote the set of strategies j 6 = i, j = 1 , ..., N . When agent i chooses its strategy K i in (37) gi ven K − i , we refer to the resulting feedback gain matrix K as { K i ; K − i } . In the multi-agent system (38), the single performance output z 2 in (1) is replaced by N indi vidual performance outputs of the agents z 2 , ( i ) . Assuming that each performance output z 2 , ( i ) = C 2 , ( i ) x + D 2 , ( i ) u i has a form that satisfies (3), the H 2 -cost from w 2 to agent i ’ s performance output can equiv alently be defined as the indi vidual LQR cost of agent i : J i ( K ) = Z ∞ t =0  x T ( t ) Q i x ( t ) + u T i ( t ) R i u i ( t )  dt s.t. w 1 ( t ) = 0 , w 2 ( t ) = δ ( t ) (39) where Q i ∈ R n × n  0 and R i ∈ R q i × q i  0 are weight matrices for state and control input of agent i , respectively , and w 2 ( t ) is an impulse disturbance. Similar to the centralized case, for stabilization of (37), the joint control strategy K needs to satisfy (6). In addition, we are interested in implementing a sparse con- troller subject to a global sparsity constraint. In the following, we develop a noncooperati ve game where agent i is modeled as a game player , with its strategy giv en by the control policy represented by K i . The joint strategies { K 1 , K 2 , ..., K N } must guarantee stability of the uncertain system (37) with at most s communication links in total. Thus, the set of admissible strategies { K 1 , K 2 , ..., K N } must satisfy || T z 1 w 1 ( { K 1 , K 2 , ..., K N } ) || ∞ < γ , (40) card( { K 1 , K 2 , ..., K N } ) ≤ s, (41) and the set of feasible strategies for player i , giv en other players’ strategies K − i , must satisfy G i ( K − i ) = { K i | card( { K i ; K − i } ) ≤ s, || T z 1 w 1 ( { K i ; K − i } ) || ∞ < γ . (42) Giv en K − i , player i solves the following optimization: min K i J i ( { K i ; K − i } ) s . t . K i ∈ G i ( K − i ) . (43) Follo wing [27], we can say that the set of strategies ( K ∗ 1 , K ∗ 2 , ..., K ∗ N ) is a Generalized Nash Equilibrium (GNE) if J i ( { K ∗ i ; K ∗ − i } ) ≤ J i ( { K i ; K ∗ − i } ) , ∀ K i ∈ G i ( K ∗ − i ) , i = 1 , ..., N . (44) In GNE, no user can unitarily deviate from the equilibrium to improv e his utility giv en that the strategy satisfies the global constraint [27]. A GNE differs from Nash equilibrium (NE) due to the presence of global constraints. B. P ALM algorithm for GNE W e propose to solve the generalized Nash strategies (44) using the best-response dynamic (43) where each player takes its turn to maximize its payoff based on other players’ strate- gies. The steps are listed in Algorithm 3. Recall Algorithm 1, where the tuple K , F was iterati vely optimized to solve the penalized optimization (18). Similarly , given K − i , player i ’ s optimization (43) can be written in the penalized form using indicator functions as min K i , F Φ i ( K i , F ; K − i ) (45) with Φ i ( K i , F ; K − i ) , J i ( { K i ; K − i } ) + h ( { K i ; K − i } ) + f ( F ) + H ( { K i ; K − i } , F ) , (46) where the indicator functions h ( · ) and f ( · ) are giv en by (16,17), and the matrix { K i ; K − i } is defined after (38) . In 6 this optimization, K i ∈ R q i × n is viewed as the feedback gain of agent i that satisfies || T z 1 w 1 ( { K i ; K − i } ) || ∞ < γ , and F ∈ R m × n represents the system-wide sparse feedback gain matrix that satisfies the global sparsity constraint. In the func- tion Φ( K , F ) (18), the variables K and F were of the same size, and they represented the same global sparsity-constrained feedback gain. Howe ver , when minimizing Φ i ( K i , F ; K − i ) (45), the variable K i is the robust feedback gain for player i while F represents the global feedback gain that satisfies the sparsity constraint. In the best-response dynamic, in each round the players take turns to minimize their own respective Φ i functions over K i and F . The equilibrium point is achieved when no player can improve its Φ i using K i and F while K j is fixed for j 6 = i . Note that in the initial best response update steps, gi ven non-sparse K − i , the minimization objective (45) cannot drive the coupling function H ( { K i ; K − i } , F ) in (20) to zero. This is because P i 6 = j || K i − ( F ) i || 2 F ≈ 0 only when K − i approaches the desired lev el of sparsity , where ( F ) i ∈ R q i × n denotes the ro w block of F that corresponds to the feedback gain of the i th player . The minimization of (45) is similar to the minimization of (18). Thus, modified Algorithm 1 is used in line 9 of Algorithm 3 to solve (45). Given its K l i , F l at iteration l , the following proximal operators are performed by player i in the minimization of line 9 of Algorithm 3. 1) F -minimization: : Compute the proximal point Z k for F k : Z k = F k − 1 a ∇ F H ( K k , F k ) = F k − γ a ( F k − { K k i ; K − i } ) (47) Solve the proximal operator: F k +1 = arg min F a 2 || F − Z k || 2 F s . t . card( F ) ≤ s, (48) and get F k +1 = [ Z k ] s , similar to Step 2 of Algorithm 1. 2) K -minimization: : Algorithm 3 P ALM algorithm for computing GNE (44) 1: Given s : global sparsity constraint, γ : H ∞ -norm bound. 2: Initialization : 3: K 0 : any stabilizing feedback gain with T ∞ ( K 0 ) < γ . 4: F 0 : any stabilizing feedback gain F 0 . 5: for l = 1 ...l max until || F 1 − F l − 1 || F <  3 do 6: K l := K l − 1 , F l := F l − 1 7: for i = 1 ...N do 8: // Solve using (47 – 50) with K l i , F l as the initial values: 9: ˆ K i , ˆ F = arg min K i , F Φ i ( K i , F ; K l − i ) 10: // Update K l and F l : 11: K l = { ˆ K i ; K l − i } 12: F l = ˆ F 13: end f or 14: end for 15: Output: K GNE ( s ) := F l . Compute proximal point X k i for K i : X k i = K k i − 1 b ∇ K i H ( K k i − ( F k +1 ) i ) = K k i − ρ b ( K k i − ( F k +1 ) i ) . (49) Solve the proximal operator: K k +1 i = arg min K i  J i ( { K i ; K − i } ) + b 2 || K i − X k i || 2 F  s . t . || T z 1 w 1 ( { K i ; K − i } ) || ∞ < γ . (50) Solving (50) is similar to solving (31). In (50), player i aims to update its control strategy K i giv en other players’ strategies K − i . Algorithm 2 is applied to solve (50) with sev eral modifications. The minimization cost in (50) is defined as h i ( K i , K − i ) , J i ( { K i ; K − i } ) + b 2 || K i − X k i || 2 F , and the gradient with respect to K i in line 4 of Algorithm 2 is replaced by ∇ K i h i ( K i , K − i ) = (51) 2( R i K i C − B T ( i ) P ( K i , K − i )) L ( K i , K − i ) C T + b ( K i − X k i ) where L ( K i , K − i ) and P ( K i , K − i ) are the solution of the following set of equations: ( ¯ A cl ( K i , K − i )) T P ( K i , K − i ) + P ( K i , K − i ) ¯ A cl ( K i , K − i ) + ¯ Q i ( K i , K − i ) = 0 ¯ A cl ( K i , K − i ) L ( K i , K − i ) + L ( K i , K − i )( ¯ A cl ( K i , K − i )) T + B 2 B T 2 = 0 (52) and ¯ A cl ( K i , K − i ) = A − X j 6 = i B ( j ) K j C − B ( i ) K i C (53) ¯ Q i ( K i , K − i ) = Q i + C T ( X j 6 = i ( K j ) T R j K j + K T i R i K i ) C . Q i and R i are defined in (39). Similar to lines 10–14 in Algo- rithm 2, when player i ’ s strategy K cur i is near the boundary of the H ∞ -norm constraint given other players’ strategies K − i , an improving feasible direction for K i can be found by solving an LMI such as (36) for scalar z and K in i ∈ R q i × n : Maximize z , K in i z s.t. trace[( ∇ K i h i ( K cur i , K − i ) T ( K in i − K cur i )] + z ≤ 0 G  { K in i ; K − i } ; { K cur i ; K − i }  − θ z · I  0 . (54) Here, θ is the factor to control the speed of reduction of H ∞ norm, and G is defined as in (34). Then the update direction of K i can be formed as ∆ K i = K in i − K cur i . W e note that in K -minimization step, each player updates its own strategy K i , while in F -minimization step, the strategies of all the players are jointly updated. Thus, Algorithm 3 has partially distributed computation. Finally , we note that the centralized problem (14) can be represented as a potential game by modifying the noncoop- erativ e game (44). A game { N , {A i } , { J i }} with N players, action set {A i } N i =1 and utilities { J i } N i =1 , is an exact potential game [35] if there exists a global function Φ , such that for ev ery player i ∈ N , a − i ∈ A − i and a 0 i , a 00 i ∈ A i , J i ( a 0 i , a − i ) − J i ( a 00 i , a − i ) = Φ( a 0 i , a − i ) − Φ( a 00 i , a − i ) . (55) W e employ a common assumption that the input penalty of each user is uncorrelated, i.e., u T Ru = P N i =1 u T i R ( i ) u i , where R ( i ) is the submatrix of R that represents the weight 7 ˙ x i = A ii x i + X j 6 = i e ˆ l ( i,j ) x j + B i u i player 1 player 2 Fig. 1: The 5-node open-loop unstable network matrix for u i ( t ) . Thus, the objectiv e J ( K ) in (13) can be expressed as J ( K ) = J ( { K i , K − i } ) = (56) Z ∞ 0   x T   Q + C T ( X j 6 = i K T j R ( j ) K j ) C   x + u T i R ( i ) u i   dt, with u i = − K i y . Thus, the minimization objectiv es J i ( { K i ; K − i } ) of all players in (43) are replaced by the global LQR cost (56). T o con vert the game in (43) into an exact potential game, we set Q i = Q + C T ( P j 6 = i K T j R ( j ) K j ) C , R i = R ( i ) in the individual cost (39), which is consistent with (56). The GNE strategies K ∗ 1 , K ∗ 2 , ..., K ∗ N for the potential game can be written as J ( { K ∗ i ; K ∗ − i } ) ≤ J ( { K i ; K ∗ − i } ) , ∀ K i ∈ G i ( K ∗ − i ) , i = 1 , ..., N . (57) W e employ Algorithm 3 to compute (57). Players update their control strategies in the K -minimization step distributi v ely , and jointly update their strategies in the F -minimization step, thereby obtaining a partially distributed implementation of the centralized sparsity-constrained problem (14). I V . N U M E R I C A L R E S U L T S A N D C O N V E R G E N C E A N A L Y S I S A. Network model W e consider an example of an uncertain network model from [36]. The network consists of N connected nodes dis- tributed randomly on a L unit by L unit square area. Each node is an unstable second-order system coupled with other nodes through an exponentially decaying function of the Euclidean distance ˆ l ( i, j ) [1], [36]. The state-space representation of node i is given as:  ˙ x 1 i ˙ x 2 i  = ˆ A ii  x 1 i x 2 i  + X j 6 = i e − ˆ l ( i,j )  x 1 j x 2 j  +  0 1  ( d i + u i ) . (58) In the abov e state-space representation, the state matrix ˆ A ii , i = 1 , ..., N and the Euclidean distances ˆ l ( i, j ) are not known exactly . In particular, ˆ A ii = A ii + A ii   θ 11 θ 12 θ 21 θ 22  ˆ l ( i, j ) = l ( i, j ) · (1 + δ i,j ) , (59) where A ii and l ( i, j ) are the nominal v alues, and δ ij and θ ij are independent random perturbations, uniformly distributed in the range ± 20% . The operator  denotes element-wise multiplication. As in (1), A denotes the nominal value of the state matrix of this N -node unstable system, and ˆ A denotes one realization of the perturbed state matrix. The uncertain matrix ∆ A = ˆ A − A in (1). The control input matrix is assumed to be known for this example, so that ˆ B = B = 1 N ⊗ B ii , where B ii =  0 1  T , and ⊗ denotes the Kronecker product [37]. In this simulation study , we collected 200 random samples of ˆ A . T o guarantee closed-loop stability of (58), we numerically compute the worst-case ˆ A as, ˆ A worst = arg max ˆ A σ max ( ˆ A − A ) . (60) Using the singular value decomposition, we obtain U S V T = ˆ A worst − A . Normalizing S by σ max ( S ) , we set B 1 = p σ max ( S ) U , C 1 = p σ max ( S ) V T in (2). Due to this normalization, γ = 1 . The following parameters are employed in the simulations. W e set L = 2 and N = 5 , thus A ∈ R 10 × 10 , B ∈ R 10 × 5 . The output matrix C = I 10 . The dense feedback matrix K has card( K ) = 50 . When the feedback controller is completely decentralized, i.e., feedback links only exist between states and controllers within the same node, card( K ) = 10 . The performance index for the LQR cost employs Q = 100 · I and R = I in (12) for the centralized problem (14). For the noncooperativ e game (44), we consider a two-player game as shown in Figure 1, where player 1 is in charge of the control inputs in nodes 1 and 3 and player 2 is in charge of the control inputs in nodes 2, 4, 5. The performance index matrices Q i , R i , i = 1 , 2 for the LQR cost in (39) satisfy: x T Q 1 x + u T 1 R 1 u 1 = 100[( x 11 − x 13 ) 2 +( x 21 − x 23 ) 2 ]+ u 2 1 + u 2 3 x T Q 2 x + u T 2 R 2 u 2 = 100 X j =2 , 4 , 5 ( x 2 1 j + x 2 2 j ) + X j =2 , 4 , 5 u 2 j . (61) W e solve all the LMIs using the CVX package [38]. B. Social optimization First, we present simulation results for the problem (14) applied to the system in (58) with γ = 1 in (14) ov er a range of s -values. W e implement Algorithm 1, with the resulting feedback matrix denoted as K ∗ palm ( s ) . For the same problem (14), we also use Algorithm 3 applied to the potential game (57), with the solution denoted by K ∗ P ALMPG ( s ) , giv en the sparsity constraint s . For comparison, we also run the GraSP algorithm that was used in [24], with the resulting feedback denoted by K ∗ GraSP ( s ) , initialized by a stabilizing decentralized controller K dec with card( K dec ) = 10 . In general, GraSP needs to be initialized by a K 0 that satisfies card( K 0 ) ≤ s and T ∞ ( K 0 ) < γ , which in reality might be difficult to find. In contrast, the P ALM-based Algorithm 1 of this paper does not rely on any such sparse initialization. Finally , we show performance of the dense mixed H 2 /H ∞ controller using the simple gradient method in [28]. Figure 2 illustrates the optimal LQR cost J in problem (14) and the associated H ∞ norm vs. sparsity constraint s . For 15 ≤ s ≤ 50 , the centralized Algorithm 1 and the potential game using Algorithm 3 both conv erge to a solution with sufficiently small coupling function in (18), which indicates F ≈ K . From Figure 2(a), we observe that the H 2 norms of all sparsity-constrained methods decrease as s is relaxed, and approach to that of the dense controller [4]. Howe ver , the 8 15 20 25 30 35 40 45 50 s 73.5 74 74.5 75 75.5 76 76.5 77 J PALM centralized (Alg. 1) PALM-Potential game (Alg. 3) GraSP [10] Gradient method [29] (a) J vs. sparsity constraint s . 15 20 25 30 35 40 45 50 s 0.8 0.85 0.9 0.95 1 1.05 H ∞ γ =1 PALM centralized (Alg.1) PALM-Potential game (Alg.3) GraSP [10] Gradient method [29] (b) H ∞ norm v .s. sparsity constraint s . Fig. 2: The LQR cost J and H ∞ norm vs. sparsity constraint s . P ALM-based methods have similar LQR costs and outperform significantly the greedy GraSP algorithm in [24]. In GraSP , the choice of activ e coordinates only depends on the gradient information of the function J . At con ver gence, the solution of the mixed H 2 /H ∞ problem has the sparsity structure gi ven by the greedy selection step. For the P ALM algorithm, since we iteratively compute the proximal map on X k and Z k , the support is chosen based on the information on both the LQR cost J ( K ) and the H ∞ -norm constraint T ∞ ( K ) . Thus, at conv ergence, the P ALM method finds a critical point of problem (14) while GraSP does not necessarily achiev e it. Figure 2(b) sho ws the H ∞ norms of K ∗ P ALM ( s ) , K ∗ P ALMPG ( s ) and K ∗ GraSP ( s ) . W e observe that for both GraSP and P ALM methods, the solution is found in the interior of the H ∞ -norm constraint for s ≥ 30 , and on the boundary for s ≤ 25 , which indicates that when the sparsity constraint becomes more stringent, satisfying the sparsity and H ∞ -norm constraints simultaneously becomes challenging. Both Algorithm 1 (the social optimization) and Algorithm 3 (the potential game) are found to con ver ge for all s -values for this system. Figure 3 shows the error in consecutiv e steps for variable K at the end of step 3 of Algorithm 1 as a function of iteration step, for different s -values. W e found that ∆ F k has a similar trend to ∆ K k . The errors in consecuti ve steps are defined as ∆ K k , K k − K k − 1 and ∆ F k , F k − F k − 1 . W e note that the error conv erges faster for larger s -values, which might be explained by the fact that that for s > 25 , the minima are found in the interior of the H ∞ -norm constraint set (see Figure 2). For Algorithm 3 (potential game), the penalized cost function Φ i and || K − F || 2 F (line 9) hav e similar trends to those for Algorithm 1. Moreover , it is demonstrated in Fig 4 that although Algorithm 1 con ver ges to a critical point of Φ ( K , F ) , the coupling function H ( K , F ) > 0 for s < 15 . As a result, when Algorithm 1 conv erges for these s -values, K 6 = F , so a sparse feedback solution that satisfies (14) cannot be found. Thus, in Figure 2, we only sho w the LQR cost and H ∞ -norm for 15 ≤ s ≤ 50 . C. The noncooperative game W e in vestigate performance of Algorithm 3 for the noncoop- erativ e game with different indi vidual costs (61) for the system (58). W e use K GNE ( s ) = { K GNE 1 ( s ) , K GNE 2 ( s ) } to denote the two players’ feedback produced by Algorithm 3 when the sparsity constraint is gi ven by s . Figure 5 shows the errors in consecutiv e steps of player i ’ s strategic v ariables K i , F i for i = 1 , 2 vs iteration round l in Algorithm 3. W e observe that 0 50 100 150 200 250 300 PALM iteration index k 10 -4 10 -2 10 0 10 2 || ∆ K k || s=5 s=15 s=25 s=29 s=31 s=40 Fig. 3: The error in K vs. iteration k in P ALM Algorithm 1 Step 2 and 3 for different s -values. 0 10 20 30 40 PALM iteration index k 0 100 200 300 400 ||K k - F k || 2 s=5 s=10 0 20 40 60 80 100 PALM iteration index k 10 -4 10 -2 10 0 10 2 ||K k - F k || 2 s=15 s=25 s=29 s=30 s=35 s=45 0 20 40 60 80 100 PALM iteration index k 10 1 10 2 10 3 10 4 Φ (K k , F k ) s=15 s=25 s=29 s=30 s=35 s=45 0 10 20 30 40 PALM iteration index k 0 0.5 1 1.5 2 Φ (K k , F k ) × 10 5 s=5 s=10 (a) (b) (c) (d) Fig. 4: The penalized cost function Φ( K k , F k ) and the coupling function || K k − F k || 2 F vs iteration k in the end of Step 3 of Algorithm 1 for multiple s -values. both || ∆ K i || F and || ∆ F i || F decrease significantly within the first 10 iterations and then saturate to small values as l grows, resulting in the saturation of the penalized cost function Φ i in line 9 of Algorithm 3, which corresponds to an approximate equilibrium point as discussed in section IV -D. The normalized coupling function 1 ρ H ( K l , F l ) = || K l − F l || 2 F (20) decreases with iteration l , following the trend in Figure 4. For s > 20 , the square error || K l − F l || 2 F reaches a sufficiently small value ( < 10 − 4 ) at the equilibrium point, while for s ≤ 20 , the square error is larger , causing T ∞ ( K GNE ( s )) > T ∞ ( K l ) , which results in T ∞ ( K GNE ( s )) > 1 during conv ergence. One way to a void this discrepency and still guarantee closed-loop stability is to replace γ in (16) with γ −  , and pro vide a margin that compensates for the square error between K and F . For this example, we set  = 0 . 01 , so that T ∞ ( K GNE ( s )) < 0 . 99 . Figure 6 illustrates the individual LQR costs J i as in (39), and the global H ∞ norm when K GNE ( s ) is implemented. W e observe that in Figure 6(a), for each player i , the LQR cost achieved at the equilibrium point J i ( K GNE ( s )) tends to decrease with s , which indicates that there is a trade-off between the selfish LQR cost and the global shared sparsity constraint. Figure 6(b) shows that T ∞ ( K GNE ( s )) < 1 for 15 ≤ s ≤ 45 , indicating that the Nash strategies in K GNE ( s ) are guaranteed to stabilize the uncertain system (58). 9 0 20 40 60 80 100 l 10 -4 10 -2 10 0 10 2 Error in consecutive steps s=15 || ∆ K 1 l || || ∆ K 2 l || || ∆ F 1 l || || ∆ F 2 l || (a) (b) 0 50 100 150 200 250 300 l 10 -4 10 -2 10 0 10 2 Error in consecutive steps s=25 || ∆ K 1 l || || ∆ K 2 l || || ∆ F 1 l || || ∆ F 2 l || Fig. 5: Errors in consecutive steps of K l i and F l i for players i = 1 , 2 vs. step l in Algorithm 3 (the noncooperative game). 15 20 25 30 35 40 45 50 s 26 28 30 J 1 (K GNE (s)) 15 20 25 30 35 40 45 50 s 44 45 46 J 2 (K GNE (s)) (a) J i ( K GNE ( s )) vs. sparsity constraint s at GNE for i = 1 , 2 . 15 20 25 30 35 40 45 50 s 0.95 0.96 0.97 0.98 0.99 1 T ∞ (K GNE (s)) 0 2 4 6 8 10 nz = 45 0 2 4 6 s=45 0 2 4 6 8 10 nz = 35 0 2 4 6 s=35 0 2 4 6 8 10 nz = 15 0 2 4 6 s=15 0 2 4 6 8 10 nz = 30 0 2 4 6 s=30 (b) T ∞ ( K GNE ( s )) vs. the sparsity constraint s at GNE. Fig. 6: The indi vidual LQR cost and global H ∞ norm of K GNE ( s ) vs sparsity constraint s at GNE and the sparsity pattern of K GNE ( s ) for dif ferent s values. D. Algorithm Con ver gence and Complexity W e close this section by providing some final comments about the con vergence properties and complexity of the pro- posed algorithms. Global conv ergence of P ALM for noncon- ve x nonsmooth functions was studied in [25], while that of P ALM-based output feedback co-design under block-sparsity constraints was established in [5]. Results in [5], [25] are extended to analyze the conv ergence properties of Algorithm 1 in [30], [31]. It has been prov en in [Bolte et al., 2014] that if Lemma 2 – 4 in [30] hold, then the sequence generated by P ALM algorithm globally conv erges. In addition, if Lemma 5 of [30] holds, the sequence con ver ges to a critical point [Bolte et al., 2014] of Φ . This confirms conv ergence of Algorithm 1 to a sparsity-constrained mixed H 2 /H ∞ controller , which corresponds to a critical point of Φ under mild assumptions on the functions J and g . Next, we briefly discuss the conv ergence properties of Algorithm 3. Suppose a GNE (44) is given by ( K ∗ 1 , ..., K ∗ N ) . Then the following condition holds for each player i [39]: ∇ G i ( K ∗ − i ) ,η J i ( { K ∗ i ; K ∗ − i } ) = 0 , (62) where ∇ G i ( K ∗ − i ) ,η J i ( { K i ; K ∗ − i } ) is the projected gradient of cost J i (39) onto the constraint set G i (42) for player i , i = 1 , ..., N . Again, from [39], we can write ∇ G i ( K ∗ − i ) ,η J i ( { K i ; K ∗ − i } ) , 1 η ( K i − Π G i ( K ∗ − i ) [ K i − η ∇ K i J i ( { K i , K ∗ − i } )]) (63) where η > 0 and the operator Π K ( · ) denotes projection onto the set K . Similarly , in line 9 of Algorithm 3, a necessary condition for Φ i to achiev e its minimum is that the projected gradient ∇ G i ( K l − i ) ,η J i ( { K i ; K l − i } ) = 0 . Instead of seeking an exact equilibrium point as GNE, we assume conv ergence of Al- gorithm 3 when this projected gradient is suf ficiently small, which is a necessary condition for an approximate local equi- librium [39]. At iteration l , the ˆ K i in line 9 can be vie wed as an approximation of Π G i ( K l − i ) [ K l − 1 i − η ∇ K i J i ( { K l − 1 i , K l − i ) } ] . Thus, the norm of the projected gradient is proportional to || K l i − K l − 1 i || , implying that small values of ∆ K l i and ∆ F l i indicate con ver gence of Algorithm 3. This is illustrated in Figure 5. Moreover , we note that there is no theoretical guarantee for the existence of GNE for the game in (43). If a GNE exists for the potential game (57) then this GNE satisfies the necessary condition for the minimizer of (14). The main numerical complexity of Algorithms 1 and 3 is dominated by the K -minimization step (Step 3 of Algorithm 1 and line 9 of Algorithm 3), which has polynomial complexity on the number of variables in the feedback matrix [40]. V . C O N C L U S I O N The P ALM method was exploited to solve the sparsity- constrained mixed H 2 /H ∞ control problem for multi-agent systems. First, a centralized social-optimization algorithm was in vestigated. Second, we dev eloped noncooperative and potential games that hav e partially-distributed computation. The proposed algorithms were validated using an open-loop unstable network dynamic system. It was demonstrated that the centralized P ALM method outperforms the GraSP-based method for most sparsity lev els, and conv erges both theo- retically as well as in simulation results. Moreover , a best- response dynamics algorithm for the proposed games con- ver ges to an approximate GNE point. The performance of the potential game for social optimization closely approximates that of the centralized algorithm. 10 R E F E R E N C E S [1] F . Lin, M. Fardad, and M. R. Jovanovi ´ c, “Design of optimal sparse feedback gains via the alternating direction method of multipliers, ” IEEE T rans. on Aut. Ctrl. , vol. 58, no. 9, pp. 2426–2431, 2013. [2] F . Lian, A. Chakrabortty , and A. Duel-Hallen, “Game-theoretic multi- agent control and network cost allocation under communication con- straints, ” IEEE Journal on Selected Areas in Communications , vol. 35, no. 2, pp. 330–340, 2017. [3] N. Monshizadeh, H. L. T rentelman, and M. K. Camlibel, “Projection- based model reduction of multi-agent systems using graph partitions, ” IEEE Tr ans. on Contr ol of Network Systems , vol. 1, no. 2, pp. 145–154, 2014. [4] F . D ¨ orfler , M. R. Jov anovi ´ c, M. Chertkov , and F . Bullo, “Sparsity- promoting optimal wide-area control of power networks, ” IEEE T r ans. on P ower Systems , vol. 29, no. 5, pp. 2281–2291, 2014. [5] F . Lin and V . Adetola, “Co-design of sparse output feedback and row/column-sparse output matrix, ” in American Control Confer ence (ACC), 2017 . IEEE, 2017, pp. 4359–4364. [6] N. Matni and V . Chandrasekaran, “Regularization for design, ” IEEE T ransactions on Automatic Contr ol , vol. 61, no. 12, pp. 3991–4006, 2016. [7] C. Lidstr ¨ om and A. Rantzer , “Optimal H ∞ state feedback for systems with symmetric and hurwitz state matrix, ” in American Contr ol Confer - ence , 2016, pp. 3366–3371. [8] R. Arastoo, M. Baha varnia, M. V . K othare, and N. Motee, “Closed-loop feedback sparsification under parametric uncertainties, ” in IEEE 55th Confer ence on Decision and Control , 2016, pp. 123–128. [9] M. Bahav arnia and N. Motee, “Sparse memoryless LQR design for uncertain linear time-delay systems, ” IF A C-P apersOnLine , vol. 50, no. 1, pp. 10 395–10 400, 2017. [10] M. Bahav arnia, C. Somarakis, and N. Motee, “State feedback controller sparsification via a notion of non-fragility , ” in Decision and Contr ol (CDC), 2017 IEEE 56th Annual Conference on . IEEE, 2017, pp. 4205– 4210. [11] S. Seuken and S. Zilberstein, “Formal models and algorithms for decentralized decision making under uncertainty , ” Autonomous Agents and Multi-Agent Systems , vol. 17, no. 2, pp. 190–250, 2008. [12] P . Ogren, E. Fiorelli, and N. E. Leonard, “Cooperativ e control of mobile sensor networks: Adaptiv e gradient climbing in a distributed en vironment, ” IEEE Tr ansactions on Automatic contr ol , vol. 49, no. 8, pp. 1292–1302, 2004. [13] M. Jungers, E. B. Castelan, E. R. De Pieri, and H. Abou-Kandil, “Bounded nash type controls for uncertain linear systems, ” Automatica , vol. 44, no. 7, pp. 1874–1879, 2008. [14] N. de la Cruz and M. Jimenez-Lizarraga, “Finite time robust feed- back nash equilibrium for linear quadratic games, ” IF AC-P apersOnLine , vol. 50, no. 1, pp. 11 794–11 799, 2017. [15] T . Bas ¸ar and G. J. Olster, Dynamic noncooperative game theory . SIAM, 1995, vol. 200. [16] H. Mukaidani, “Robust guaranteed cost control for uncertain stochastic systems with multiple decision makers, ” Automatica , vol. 45, no. 7, pp. 1758–1764, 2009. [17] ——, “H 2 /H ∞ control problem for stochastic delay systems with multiple decision makers, ” in Decision and Control (CDC), 2014 IEEE 53r d Annual Conference on . IEEE, 2014, pp. 2648–2653. [18] H. Mukaidani and H. Xu, “Stackelberg strategies for stochastic systems with multiple followers, ” Automatica , vol. 53, pp. 53–59, 2015. [19] H. Mukaidani, H. Xu, and V . Dragan, “Dynamic games for markov jump stochastic delay systems, ” in Recent Results on Time-Delay Systems . Springer , 2016, pp. 207–227. [20] D. Vrabie and F . Lewis, “Integral reinforcement learning for online computation of feedback nash strategies of nonzero-sum differential games, ” in Decision and Control (CDC), 2010 49th IEEE Conference on . IEEE, 2010, pp. 3066–3071. [21] K. G. V amvoudakis, “Non-zero sum Nash Q-learning for unkno wn deterministic continuous-time linear systems, ” Automatica , vol. 61, pp. 274–281, 2015. [22] R. Song, F . L. Lewis, and Q. W ei, “Off-policy integral reinforcement learning method to solv e nonlinear continuous-time multiplayer nonzero- sum games, ” IEEE transactions on neural networks and learning sys- tems , vol. 28, no. 3, pp. 704–713, 2017. [23] K. G. V amvoudakis, H. Modares, B. Kiumarsi, and F . L. Lewis, “Game theory-based control system algorithms with real-time reinforcement learning: How to solve multiplayer games online, ” IEEE Contr ol Sys- tems , vol. 37, no. 1, pp. 33–52, 2017. [24] F . Lian, A. Chakrabortty , F . W u, and A. Duel-Hallen, “Sparsity- constrained mixed H 2 /H ∞ control, ” in American Control Conference (ACC), 2018 , 2018. [25] J. Bolte, S. Sabach, and M. T eboulle, “Proximal alternating linearized minimization or nonconv ex and nonsmooth problems, ” Mathematical Pr ogramming , vol. 146, no. 1-2, pp. 459–494, 2014. [26] S. Bahmani, B. Raj, and P . T . Boufounos, “Greedy sparsity-constrained optimization, ” The Journal of Mac hine Learning Resear ch , vol. 14, no. 1, pp. 807–841, 2013. [27] D. Paccagnan, B. Gentile, F . Parise, M. Kamgarpour, and J. L ygeros, “Distributed computation of generalized nash equilibria in quadratic aggregati ve games with affine coupling constraints, ” in Decision and Contr ol (CDC), 2016 IEEE 55th Confer ence on . IEEE, 2016, pp. 6123–6128. [28] Y . Kami and E. Nobuyama, “ A gradient method for the static output feedback mixed H 2 /H ∞ control, ” IF AC Pr oceedings V olumes , vol. 41, no. 2, pp. 7838–7842, 2008. [29] N. Parikh, S. Boyd et al. , “Proximal algorithms, ” F oundations and T r ends R  in Optimization , vol. 1, no. 3, pp. 127–239, 2014. [30] F . Lian, A. Chakrabortty , and A. Duel-Hallen, “Supplementary materials for ‘Game-theoretic mixed h 2 /h ∞ control with sparsity constraint for multi-agent networked control systems’. ” [31] F . Lian, “Communication-cost-constrained algorithms and games for multi-agent control systems, ” Ph.D. dissertation, North Carolina State Univ ersity , Raleigh, NC, USA, 2019. [32] D. G. Luenberger and Y . Y e, Linear and nonlinear pro gramming . Springer , 1984, vol. 2. [33] M. S. Bazaraa, H. D. Sherali, and C. M. Shetty , Nonlinear pr ogramming: theory and algorithms . John Wile y & Sons, 2013. [34] M. Saeki, “Static output feedback design for H ∞ control by descent method, ” in 45th IEEE Confer ence on Decision and Control , 2006, pp. 5156–5161. [35] N. Li and J. R. Marden, “Designing games for distributed optimization, ” IEEE J ournal of Selected T opics in Signal Processing , vol. 7, no. 2, pp. 230–242, 2013. [36] N. Motee and A. Jadbabaie, “Optimal control of spatially distributed systems, ” IEEE T rans. on Aut. Ctrl. , vol. 53, no. 7, pp. 1616–1629, 2008. [37] C. D. Meyer , Matrix analysis and applied linear algebra . Siam, 2000, vol. 71. [38] M. Grant and S. Boyd, “CVX: Matlab software for disciplined con ve x programming, version 2.1, ” http://cvxr .com/cvx, Mar. 2014. [39] E. Hazan, K. Singh, and C. Zhang, “Efficient regret minimization in non-con vex games, ” arXiv preprint , 2017. [40] P . Gahinet, A. Nemirovskii, A. J. Laub, and M. Chilali, “The lmi control toolbox, ” in Pr oceedings of 1994 33r d IEEE Confer ence on Decision and Contr ol , vol. 3. IEEE, 1994, pp. 2038–2041. 11 Supplemental Materials for “Game-Theoretic Mix ed H 2 /H ∞ Control with Sparsity Constraint for Multi-agent Network ed Control Systems” by Feier Lian, Aran ya Chakrabortty , and Ale xandra Duel-Hallen S . I : O V E RV I E W O F Z O U T E N D I J K ’ S M E T H O D The Zoutendijk’ s method [Bazaraa et al., 2013] is an ap- proach to constrained optimization, where an improving fea- sible direction is generated by solving a subproblem, usually a linear program. W e briefly overvie w Zoutendijk’ s method for the case of nonlinear inequality constraints. Consider the following constrained optimization problem: Minimize f ( x ) s.t. g i ( x ) ≤ 0 , i = 1 , ..., m , (S1) where x ∈ R n × 1 and f ( x ) and g i ( x ) are differentiable at x . At point x , I is the set of active constraint I = { i | g i ( x ) = 0 } . An improving feasible direction d can be found by the following linear programming problem [Bazaraa et al., 2013]: Maximize z , d z s. t. ∇ f ( x ) T d + z ≤ 0 , ∇ g i ( x ) T d + z ≤ 0 ∀ i ∈ I , − 1 ≤ d j ≤ 1 , ∀ j = 1 , ..., n, (S2) where the third normalizing constraint prevents the optimal z from approaching ∞ . It was shown [Bazaraa et al., 2013] that if the optimal v alue of z , denoted as z ∗ , satisfies z ∗ > 0 , then d is an improving direction since d satisfies ∇ f ( x ) T d < 0 and ∇ g i ( x ) T d < 0 ∀ i ∈ I . Otherwise if z ∗ = 0 , then the current x is a Fritz John point [Bazaraa et al., 2013], which satisfies the necessary condition for the local minimum of (S1). S . I I : D E FI N I T I O N S O F T E R M S I N S E C T I O N I I - C Definition 1. (Lipschitz constant) A function f : R d → R with the gradient function ∇ f is Lipschitz continuous with Lipschitz constant L on S ∈ R d if ||∇ f ( x ) − ∇ f ( y ) || ≤ L || x − y || for all x , y ∈ S [Bazaraa et al., 2013]. Definition 2. (Proper) The function σ : S → R is a proper function if σ ( x ) > −∞ for all x ∈ S , and σ ( x ) < ∞ for at least one point x ∈ S . Definition 3. (Lower semicontinuous) The function σ : S → R is lo wer semicontinuous at ¯ x ∈ S if for all  > 0 there e xists a δ such that x ∈ S and || x − ˆ x || < δ imply σ ( x ) − σ ( ¯ x ) > −  . S . I I I : N O T A T I O N U S E D I N K U R DY K A - Ł O JA S I E W I C Z ( K L ) P RO P E RT Y , E M P L OY E D I N C O N V E R G E N C E A N A L Y S I S O F A L G O R I T H M 1 Definition 4. (Distance.) For any subset S ⊂ R d and any point x ∈ R d , the distance from x to S is defined and denoted by dist( x , S ) := inf {|| y − x || : y ∈ S } . (S3) When S = ∅ , we hav e dist( x , S ) = ∞ for all x . Let η ∈ [0 , ∞ ] . W e denote by Φ η the class of all concave and continuous functions ϕ : [0 , η ) → R + which satisfy the following conditions: (i) ϕ (0) = 0 . (ii) ϕ has first-order continous deriv ati ve on (0 , η ) and conti- nous at 0 ; (iii) for all s ∈ (0 , η ) : ϕ 0 ( s ) > 0 . For proper and lower semicontinous functions, the subdiffer- entials are defined below [Bolte et al., 2014]: Definition 5. (Subdifferentials) Let σ : R d → ( −∞ , ∞ ] be a proper and lower semicontinous function. (i) For a giv en x ∈ dom σ , the Fr ´ echet subdifferential of σ at x , written as ˆ ∂ σ ( x ) , is the set of all vectors u ∈ R d which satisfy lim y 6 = x inf y → x σ ( y ) − σ ( x ) − h u , y − x i || y − x || ≥ 0 . (S4) When x / ∈ dom σ , we set ˆ ∂ σ ( x ) = ∅ . (ii) The limiting subdifferential, or subdifferential, of σ at x ∈ R d , written ∂ σ ( x ) , is defined as ∂ σ ( x ) :=  u ∈ R d : ∃ x k → x , σ ( x k ) → σ ( x ) (S5) and u k ∈ ˆ ∂ σ ( x k ) → u as k → ∞ o . (S6) Points whose subdifferentials contains 0 are called (limiting- )critical points. Definition 6. (Kurdyka-Łojasie wicz (KL) Property) Let σ : R d → ( −∞ , + ∞ ] be proper and lower semicontinuous. (i) The function σ is said to have the Kurdyka-Łojasie wicz (KL) Property at ¯ u ∈ dom ∂ σ := { u ∈ R d : ∂ σ ( u ) 6 = ∅} if there exist η ∈ (0 , ∞ ] , a neighorhood U of ¯ u and a function ϕ ∈ Φ η , such that for all u ∈ U ∩ [ σ ( ¯ u ) < σ ( u ) < σ ( ¯ u ) + η ] , (S7) 12 the follo wing inequality holds ϕ 0 ( σ ( u ) − σ ( ¯ u )) dist(0 , ∂ σ ( u )) ≥ 1 . (S8) (ii) If σ satisfies the KL property at each point of dom ∂ σ , then σ is called a KL function. It is shown in [Bolte et al., 2014] that KL functions arise in many applications for optimization, in particular , semi- algebraic functions are KL functions. The definitions for semi- algebraic function is giv en as follo ws. Definition 7. (Semi-algebraic sets and functions). (i) A subset S ∈ R d is real semi-algebraic set if there exists a finite number of real polynomial functions g ij , h ij : R d → R such that S = ∪ p j =1 ∩ q i =1  u ∈ R d : g ij ( u ) = 0 and h ij ( u ) < 0  . (S9) (ii) A function h : R d → ( −∞ , + ∞ ] is called semi-algebraic if its graph  ( u , t ) ∈ R d +1 : h ( u ) = t  (S10) is a semi-algebraic subset of R d +1 . S . I V: P R O O F O F G L O B A L C O N V E R G E N C E O F A L G O R I T H M 1 In this section, we employ results from [Lin and Adetola, 2017], [Bolte et al., 2014] to analyze conv ergence of Algo- rithm 1. T o simplify notation, we define ˜ g ( K ) , J ( K ) + g ( K ) , where J ( K ) is the performance index of H 2 cost (10) and g ( K ) is the indicator function for the H ∞ constraint (16). Lemma 2: ˜ g : R m × p → ( −∞ , ∞ ] and f : R m × p → ( −∞ , ∞ ] are proper and lo wer semicontinuous functions. Pr oof. In problem (18) the function J ( K ) is the LQR cost when using the feedback gain K . Clearly f ( K ) > −∞ , and J ( K ) < + ∞ if K is stabilizing, thus function J is proper . In addition, J ( K ) is continuous in K [Rautert and Sachs, 1997], and therefore lower semicontinuous. The function g ( K ) in (16) is the indicator function for the lev el set K ( γ ) in (33), and thus can take either 0 or + ∞ , with g ( K ) = 0 whenev er K ∈ K ( γ ) . Thus, ˜ g ( K ) is proper . In addition, g ( K ) is an indicator function of an open set, thus it is lower semicontinuous. Giv en J and g are both proper and lower semicontinuous, the summation ˜ g = J + g is proper and lower semicontinuous. Similarly , f ( F ) in (17) is a proper function. Moreover , it is sho wn in [Bolte et al., 2014] that it is lo wer semicontinuous. Lemma 3: H : R m × p × R m × p → R is a continuously differentiable function, i.e., H ∈ C 1 . Pr oof. The gradient of H ( K , F ) (27) is countinous in K , F . Thus, H ∈ C 1 . Lemma 4: (i) inf R m × p , R m × p Φ > −∞ , inf R m × p f > −∞ , and inf R m × p ˜ g > −∞ , where Φ is gi ven by (19). (ii) The partial gradient ∇ K H ( K , F ) is globally Lipschitz with moduli L 1 ( F ) , that is [Bolte et al., 2014], ||∇ K H ( K 1 , F ) − ∇ K H ( K 2 , F ) || ≤ L 1 ( F ) || K 1 − K 2 || . Like wise, the partial gradient ∇ F H ( K , F ) is globally Lips- chitz with moduli L 2 ( K ) . (iii) There exist bounds λ − i , λ + i , i = 1 , 2 such that inf { L 1 ( F k ) : k ∈ N } ≥ λ − 1 , inf { L 2 ( K k ) : k ∈ N } ≥ λ − 2 sup { L 1 ( F k ) : k ∈ N } ≤ λ + 1 , sup { L 2 ( K k ) : k ∈ N } ≤ λ + 2 (S11) (iv) ∇ H , ( ∇ K H , ∇ F H ) is Lipschitz continuous [Luen- berger and Y e, 1984] on bounded subsets of R m × p × R m × p . That is, for each bounded subset B 1 × B 2 of R m × p × R m × p there exists M > 0 , such that for all ( K i , F i ) ∈ ( B 1 , B 2 ) , ||∇ K H ( K 1 , F 1 ) − ∇ K H ( K 2 , F 2 ) || 2 F + ||∇ F H ( K 1 , F 1 ) − ∇ F H ( K 2 , F 2 ) || 2 F ≤ M ( || K 1 − K 2 || 2 F + || F 1 − F 2 || 2 F ) (S12) Pr oof. (i)–(iv) are stated as assumptions in [Bolte et al., 2014]. W e show that these assumptions hold for our sparsity- constrained mixed H 2 /H ∞ problem. It is easy to see that (i) holds since f and g are indicator functions. Since J is the LQR performance index, J ( K ) > 0 . Thus ˜ g ( K ) > −∞ . In (20), H ( K , F ) ≥ 0 , so Φ( K , F ) > −∞ . Properties (ii) and (iii) require the partial gradient of H to be globally Lipschitz, and the Lipschitz constant be upper and lo wer bounded, which is easy to verify since L 1 ( F k ) = L 2 ( K k ) = ρ (27). Property (iv) holds since the left-hand side of (S12) can be expressed as: LH S = 2 ρ 2 || ( K 1 − K 2 ) − ( F 1 − F 2 ) || 2 F ≤ 4 ρ 2 ( || K 1 − K 2 || 2 F + || F 1 − F 2 || 2 F ) . Assumption 1: Function J is a semi-algebraic function [Bolte et al., 2014]. Lemma 5: The objecti ve function Φ of (18) is a Kurdyka- Łojasiewicz (KL) function [Bolte et al., 2014]. Remark 1. A broad class of functions satisfy the semi- algebraic property , including polynomial functions, ` 0 -norm function and indicator function of positiv e semidefinite cones [Bolte et al., 2014]. The function f ( F ) is the indicator function for the semi-algebraic set { F | card( F ) ≤ s } . Thus, function f is semi-algebraic [Lin and Adetola, 2017], [Bolte et al., 2014]. The function g ( K ) is the indicator function for the le vel set K ( γ ) , which is approximated by the conv ex lev el set ˆ K ( γ 0 ) , represented by the LMI (34), and ˆ K ( γ 0 ) is a semi-algebraic set [Netzer, 2016]. The coupling function H is polynomial, so it is semi-algebraic [Bolte et al., 2014]. Moreov er , J is a semi-algebraic function by Assumption 1. Thus, each term of Φ is semi-algebraic, and since a finite sum of semi-algebraic functions is also semi-algebraic, Φ is semi-algebraic. It is shown in Theorem 5.1 in [Bolte et al., 2014] that a semi- algebraic function satisfies the KL property at any point in its domain. Thus, Φ is KL. It has been proved in [Bolte et al., 2014] that if Lemma 2 – 4 hold, then the sequence generated by P ALM algorithm globally con ver ges. In addition, if Lemma 5 holds, the se- quence conv erges to a critical point [Bolte et al., 2014] of Φ . This confirms con ver gence of Algorithm 1 to a sparsity- constrained mixed H 2 /H ∞ controller , which corresponds to a critical point of Φ under mild assumptions on the functions J and g . 13 R E F E R E N C E [Bazaraa et al., 2013] Bazaraa, M. S., Sherali, H. D., and Shetty , C. M. (2013). Nonlinear progr amming: theory and algorithms . John W iley & Sons. [Bolte et al., 2014] Bolte, J., Sabach, S., and T eboulle, M. (2014). Proximal alternating linearized minimization or noncon vex and nonsmooth problems. Mathematical Pro gramming , 146(1-2):459–494. [Lin and Adetola, 2017] Lin, F . and Adetola, V . (2017). Co-design of sparse output feedback and row/column-sparse output matrix. In American Contr ol Conference (ACC), 2017 , pages 4359–4364. IEEE. [Luenberger and Y e, 1984] Luenberger , D. G. and Y e, Y . (1984). Linear and nonlinear pro gramming , volume 2. Springer. [Netzer , 2016] Netzer, T . (2016). Real algebraic geometry and its applica- tions. arXiv pr eprint arXiv:1606.07284 . [Rautert and Sachs, 1997] Rautert, T . and Sachs, E. W . (1997). Computa- tional design of optimal output feedback controllers. SIAM Journal on Optimization , 7(3):837–852.