A scalable controller synthesis method for the robust control of networked systems
aa r X i v : . [ m a t h . O C ] S e p A scalable controller synthesis method for therobust control of networked systems
Pieter Appeltans and Wim MichielsDepartment of Computer Science, KU Leuven { pieter,wim } . { appeltans,michiels } @cs.kuleuven.be Abstract
This manuscript discusses a scalable controller synthesismethod for networked systems with a large number of identical subsys-tems based on the H-infinity control framework. The dynamics of theindividual subsystems are described by identical linear time-invariant de-lay differential equations and the effect of transport and communicationdelay is explicitly taken into account. The presented method is based onthe result that, under a particular assumption on the graph describingthe interconnections between the subsystems, the H-infinity norm of theoverall system is upper bounded by the robust H-infinity norm of a singlesubsystem with an additional uncertainty. This work will therefore brieflydiscuss a recently developed method to compute this last quantity. Theresulting controller is then obtained by directly minimizing this upperbound in the controller parameters.
Keywords — Robust control, H-infinity norm, Distributed control, Networked systems,Spectral value set, Stability radius
AMS subject classifications — 93B36, 93B70, 93A14, 93C23, 93C73
This manuscript presents a controller synthesis method for networked systems.Such networked systems consist of a large number of smaller subsystems thatinteract over a network. The analysis and control of these networked systemsis challenging due to their large dimension and the presence of delays. Thesedelays originate from the time needed to transfer matter, energy and informa-tion between subsystems. In this context, the traditional approach of usingone global controller for the complete network is thus not feasible due to thehigh communication requirements and the poor scalability with respect to thenumber of subsystems. Furthermore, the assumption that all measurements arecentrally available, does often not hold for networked systems. These limita-tions inspired local control approaches, in which each subsystem has its ownlocal controller. Neighboring controllers can however communicate to improvecontrol performance.In this manuscript we consider networked systems in which the dynamics ofthe individual subsystems are described by identical linear time-invariant delaydifferential equations. The resulting local controllers are identical and minimize1n upper bound for the H ∞ -norm of the overall system. This H ∞ -norm is animportant performance measure in robust control theory, see [18].The computation cost of the standard algorithms for calculating the H ∞ -norm of dynamical systems with discrete delays, such as [7], scales cubically withrespect to the number of states (and hence the number of subsystems). However,for some networked systems this computation cost can be decreased significantlyusing the decoupling transformation presented in [11] and [4]. More specifically,if the subsystems are identical and the graph describing the interconnectionsbetween the subsystems fulfills a particular assumption, then the H ∞ -norm ofthe complete system is equal to the maximal H ∞ -norm of a single parametrizedsubsystem where the allowable values of the parameter correspond to the eigen-values of the adjacency matrix of the interconnection graph. Moreover, in [3]it was suggested to consider this parameter as an uncertainty bounded to aregion in the complex plane that comprises all these eigenvalues. As such, theworst-case H ∞ -norm of this uncertain subsystem gives an upper bound for the H ∞ -norm of the complete network. Furthermore, this worst-case H ∞ -norm isalso known as the robust H ∞ -norm and can be computed at a cost that only de-pends on the dimension of an individual subsystem using the method presentedin [1].For the controller synthesis, we will directly minimize the robust H ∞ -normof the uncertain subsystem in the controller parameters. Our method thus fits inthe frequency based, direct optimization framework, used in [7],[12],[3] and [14].This framework allows to easily incorporate constraints on the structure of thecontroller, such as PID or reduced order control. In contrast, H ∞ -controller de-sign methods based on Ricatti equations and linear matrix inequalities typicallygive rise to dense controllers with dimensions equal to that of the system. A no-table exception is [9], which allows to design reduced order controllers. Anotheradvantage of the direct optimization approach compared to methods based onRicatti equations and linear matrix inequalities, in particular for systems withdelays, is that the obtained results are less conservative. This comes however atthe cost of having to solve a non-convex and non-smooth optimization problem.The remainder of this work is structured as follows. First, Section 2 in-troduces the considered networked systems and details the aforementioned de-coupling transformation. Next, a recently developed method to compute therobust H ∞ -norm of uncertain linear time-invariant systems with discrete delaysis discussed in Section 3. Subsequently, the direct optimization approach tosynthesize the controller is outlined in Section 4. Finally, the resulting designmethodology is illustrated using an example problem in Section 5 and someconcluding remarks are given in Section 6. H ∞ -norm of networked sys-tems In Section 2.1 we introduce the considered networked systems and the controlobjective. Section 2.2 presents the decoupling transformation that allows tocompute an upper bound for the H ∞ -norm of the overall system at a computa-tion cost that does not depend on the number of subsystems.2 .1 System Description and Control Objective Here, we consider networked systems with N subsystems. The dynamics of theindividual subsystems are identical and described by a state-space representa-tion of the following form: ˙ x j ( t ) = K P k =0 A k x j ( t − τ k ) + B u u j ( t − τ u ) + B u n u nj ( t ) + B w w j ( t ) y j ( t ) = C y x j ( t ) y nj ( t ) = C y n x j ( t ) z j ( t ) = C z x j ( t ) for j = 1 , . . . , N (1)with x j ( t ) ∈ R n the state vector of subsystem j , u j ( t ) ∈ R m c its control input, y j ( t ) ∈ R p c its measured output, w j ( t ) ∈ R m its performance input, z j ( t ) ∈ R p its performance output, 0 = τ < τ < · · · < τ K discrete delays, τ u ≥ A k , B u , B u n , B w , C y , C y n and C z real-valued matrices of appropriatedimension. The input u nj ( t ) and the output y nj ( t ) model the interactions betweenthe subsystems: u nj ( t ) = N X i =1 P Nj,i y ni ( t − τ n ),with τ n ≥ P N = [ P Nj,i ] Nj,i =1 the adjacency matrixof the interconnection graph of the network. More specifically, an element P Nj,i is non-zero if and only if the dynamics of subsystem j are influenced bysubsystem i .Each subsystem is controlled using a local controller and these local con-trollers are identical. Here we will consider dynamic output feedback controllersof order n c : ( ˙ ξ j ( t ) = J p ξ j ( t ) + F p y j ( t ) + F n p u ncj ( t ) u j ( t ) = L p ξ j ( t ) + K p y j ( t ) + K n p u ncj ( t ) for j = 1 , . . . , N (2)with ξ j ( t ) ∈ R n c the controller state of the local controller associated with sub-system j . The matrices J p , F p , F n p , L p , K p and K n p are real-valued and ofappropriated dimension. The subscript p is used to indicate that these matri-ces depend on some tunable control parameters p . If F n p = 0 and/or K n p = 0,the local controllers can communicate their sensor measurements to neighboringsubsystems. It is however required that the adjacency matrix of the communi-cation graph is the equal to the adjacency matrix of the interaction graph: u ncj ( t ) = N P i =1 P Nj,i y i ( t − τ nc ) , with τ nc ≥ Remark 1.
In the remainder of this manuscript we restrict our attention tocontroller architectures where the local controllers can only share sensor mea-surements. Note however that the results can be extended to architectureswhere the local controllers can also share their internal state.3y eliminating the control and coupling variables, we find the followingstate-space description for the closed-loop of the complete networked system: ˙ x ( t ) = K P k =0 ( I N ⊗ A k ) x ( t − τ k ) + ( I N ⊗ B u K p C y ) x ( t − τ u ) +( P N ⊗ B u n C y n ) x ( t − τ n ) + (cid:0) I N ⊗ B u L p (cid:1) ξ ( t − τ u ) + (cid:0) P N ⊗ B u K n p C y (cid:1) x ( t − τ u − τ nc ) + ( I N ⊗ B w ) w ( t )˙ ξ ( t ) = ( I N ⊗ J p ) ξ ( t ) + ( I N ⊗ F p C y ) x ( t ) + (cid:0) P N ⊗ F n p C y (cid:1) x ( t − τ nc ) z ( t ) = ( I N ⊗ C z ) x ( t ) (3)with I N the identity matrix of size N , x ( t ) = [ x ( t ) T · · · x N ( t ) T ] T the com-bined state, ξ ( t ) = [ ξ ( t ) T · · · ξ N ( t ) T ] T the combined controller state, w ( t ) =[ w ( t ) T · · · w N ( t ) T ] T the combined performance input, z ( t ) = [ z ( t ) T · · · z N ( t ) T ] T the combined performance output and ⊗ the Kro-necker product. The corresponding transfer function from w to z is equal to: T ( s ; p , N ) = (cid:0) I N ⊗ (cid:2) C z (cid:3)(cid:1) (cid:0) I N ( n + n c ) s − I N ⊗ Q p ( s ) − P N ⊗ R p ( s ) (cid:1) − × (cid:16) I N ⊗ (cid:20) B w (cid:21) (cid:17) (4)with Q p ( s ) = (cid:20) A F p C y J p (cid:21) + K X k =1 (cid:20) A k
00 0 (cid:21) e − sτ k + (cid:20) B u K p C y B u L p (cid:21) e − sτ u and R p ( s ) = (cid:20) B u n C y n
00 0 (cid:21) e − sτ n + (cid:20) B u K n p C y
00 0 (cid:21) e − s ( τ u + τ nc ) + (cid:20) F n p C y (cid:21) e − sτ nc . If system (3) is exponentially stable, the H ∞ -norm of (4) equals: k T ( · ; p , N ) k H ∞ = max ω ∈ R + σ (cid:0) T ( ω ; p , N ) (cid:1) with σ ( · ) the largest singular value of its matrix argument [7]. Here we recallthat the H ∞ -norm is an important performance measure in robust control the-ory, used to asses the disturbance rejection of a dynamical system as it givesthe worst-case energy gain of the system with respect to energy-bounded noisesignals: k T ( · ; p , N ) k H ∞ = max w ∈ L m k z k L p k w k L m , with k w k L m = qR + ∞ k w ( t ) k dt , L m = { w : [0 , + ∞ ) R m such that k w k L m < + ∞} , k w ( t ) k the Euclidean norm, and k z k L p and L p defined analo-gously [18]. 4 .2 The robust H ∞ -norm of subsystem as upper bound forthe H ∞ -norm of the overall network In this subsection we show that under the following assumption on P N , thereexists a decoupling transformation that allows to compute an upper bound forthe H ∞ -norm of (4) at a computation cost that does not depend on the numberof subsystems. Assumption 1.
The matrix P N has real-valued eigenvalues confined to aninterval [ a, b ] and is diagonalizable by a unitary matrix V N , i.e. V N H P N V N = Λ N ,with Λ N = diag( λ , λ , . . . , λ N ) and λ j ∈ [ a, b ] for j = 1 , . . . , N .If we apply the following change of variables to the states, the controllerstates, the performance input and the performance output of system (3)¯ x ( t ) = ( V N H ⊗ I n ) x ( t )¯ ξ ( t ) = ( V N H ⊗ I n c ) ξ ( t )¯ w ( t ) = ( V N H ⊗ I m ) w ( t )¯ z ( t ) = ( V N H ⊗ I p ) z ( t )we obtain ˙¯ x ( t ) = K P k =0 ( I N ⊗ A k ) ¯ x ( t − τ k ) + ( I N ⊗ B u K p C y ) ¯ x ( t − τ u ) +(Λ N ⊗ B u n C y n ) ¯ x ( t − τ n ) + ( I N ⊗ B u L p ) ¯ ξ ( t − τ u )+(Λ N ⊗ B u K n p C y ) ¯ x ( t − τ u − τ nc ) + ( I N ⊗ B w ) ¯ w ( t )˙¯ ξ ( t ) = ( I N ⊗ J p ) ¯ ξ ( t ) + ( I N ⊗ F p C y ) ¯ x ( t )+ (cid:0) Λ N ⊗ F n p C y (cid:1) ¯ x ( t − τ nc )¯ z ( t ) = ( I N ⊗ C z ) ¯ x ( t ). (5)Notice that all matrices in (5) are block diagonal and hence the behavior of thistransformed system is fully characterized by its N independent subsystems.This leads to the following theorem. Theorem 2.
For a networked system of form (3) whose adjacency matrix fulfillsAssumption 1, it holds that k T ( · ; p , N ) k H ∞ = k ¯ T ¯ w ¯ z ( · ; p , N ) k H ∞ = max λ ∈{ λ ,...,λ N } k ˆ T ˆ w ˆ z ( · ; p , λ ) k H ∞ with ¯ T ¯ w ¯ z ( · ; p , N ) the transfer function from ¯ w to ¯ z of system (5) and ˆ T ˆ w ˆ z ( · ; p , λ ) the transfer function from ˆ w to ˆ z of the following system param-eterized in λ : ˙ˆ x ( t ) = K P k =0 A k ˆ x ( t − τ k ) + B u K p C y ˆ x ( t − τ u ) + λB u n C y n ˆ x ( t − τ n )+ λB u K n p C y ˆ x ( t − τ u − τ nc ) + B u L p ˆ ξ ( t − τ u ) + B w ˆ w ( t )˙ˆ ξ ( t ) = J p ˆ ξ ( t ) + F p C y ˆ x ( t ) + λF n p C y ˆ x ( t − τ nc )ˆ z ( t ) = C z ˆ x ( t ) . (6)5 roof. The provided proof is added for self-containedness and is similar to theones given in [11] and [4]. We refer to these papers for more details.The relation between T ( ω ; p , N ) and ¯ T ¯ w ¯ z ( ω ; p , N ) is given by T ( ω ; p , N ) = ( V N ⊗ I p ) ¯ T ¯ w ¯ z ( ω ; N ) (cid:0) V N H ⊗ I m (cid:1) .Because ( V N ⊗ I p ) and ( V N H ⊗ I m ) are unitary matrices if follows that, σ (cid:0) T ( ω ; p , N ) (cid:1) = σ (cid:16) ( V N ⊗ I p ) ¯ T ¯ w ¯ z ( ω ; p , N ) (cid:0) V N H ⊗ I m (cid:1) (cid:17) = σ (cid:0) ¯ T ¯ w ¯ z ( ω ; p , N ) (cid:1) . The second equality follows from the fact that¯ T ¯ w ¯ z ( ω ; p , N ) = blkdiag j =1 ,...,N (cid:0) ˆ T ˆ w ˆ z ( ω ; p , λ j ) (cid:1) ,and hence σ (cid:0) ¯ T ¯ w ¯ z ( ω ; p , N ) (cid:1) = max λ ∈{ λ ,...,λ N } σ (cid:0) ˆ T ˆ w ˆ z ( ω ; p , λ ) (cid:1) , which concludes the proof.Note that system (6) corresponds to a single subsystem in (3) where thenetwork connections are replaced by a parameter. By treating λ in (6) as anuncertainty confined to the interval [ a, b ], the robust H ∞ -norm associated with(6), which is defined as the maximal value of the H ∞ -norm over all instances ofthe uncertain parameter, k ˆ T ˆ w ˆ z ( · ; p , · ) k [ a, b ] H ∞ = max λ ∈ [ a, b ] k ˆ T ˆ w ˆ z ( · ; p , λ ) k H ∞ , (7)can be used as an upper bound for the H ∞ -norm of (3), as stated in the followingcorollary. Corollary 3.
The H ∞ -norm of a networked system of form (3) whose adja-cency matrix fulfills Assumption 1, is upper bounded by the robust H ∞ -norm of ˆ T ˆ w ˆ z ( · ; p , λ ) , with λ an uncertain parameter confined to [ a, b ] : k T ( · ; p , N ) k H ∞ ≤ k ˆ T ˆ w ˆ z ( · ; p , · ) k [ a, b ] H ∞ . Furthermore, if Assumption 1 holds with a and b independent of N , then k ˆ T ˆ w ˆ z ( · ; p , · ) k [ a, b ] H ∞ is also an upper bound for the supremum of k T ( · ; p , N ) k H ∞ over the number of subsystems: sup N =1 ,..., + ∞ k T ( · ; p , N ) k H ∞ ≤ k ˆ T ˆ w ˆ z ( · ; p , λ ) k [ a, b ] H ∞ . If, furthermore, the (Hausdorff ) distance between [ a, b ] and + ∞ S N =1 { λ ∈ C : det( I N λ − P N ) = 0 } goes to zero then sup N =1 ,..., + ∞ k T ( · ; p , N ) k H ∞ = k ˆ T ˆ w ˆ z ( · ; p , λ ) k [ a, b ] H ∞ . N − N N Figure 1: Bidirectional ring topol-ogy Figure 2: Bidirectional line topol-ogy
Example 4.
To illustrate the applicability of this result, we consider the fol-lowing adjacency matrices: P ring N = . . . . . .
5. . . . . . . . .0 . . . . and P line N = . . . . .
5. . . . . . . . .0 . . . . The first adjacency matrix, P ring N , corresponds to a bidirectional ring topology,see Figure 1; the second one corresponds to a bidirectional line topology, see Fig-ure 2. The eigenvalues of these adjacency matrices are n cos (cid:16) π ( j − N (cid:17)o ⌊ N +22 ⌋ j =1 and n cos (cid:16) jπN +1 (cid:17)o Nj =1 , respectively. The eigenvalues ofboth matrices are thus confined to the interval [ a, b ] = [ − ,
1] for all
N > H ∞ -norm of (4) that holds for all N > H ∞ -norm This section introduces a numerical algorithm to efficiently compute the robust H ∞ -norm of an uncertain linear time-invariant system with discrete delays: ˙ x ( t ) = R P r =0 ( H r + λ G r ) x ( t − τ r ) + B w w ( t ) z ( t ) = C z x ( t ) (8)with x ( t ) ∈ R n the state, w ( t ) ∈ R m performance input, z ( t ) ∈ R p the per-formance output, 0 = τ < τ < · · · < τ R discrete delays, H r , G r , B w and C z real-valued matrices of appropriate dimension and λ a real-valued, scalar7arameter addressed as an uncertainty confined to the interval [ a, b ]. Note thatsystem (6) fits this form.Under the assumption that system (8) is internally exponentially stable forall λ ∈ [ a, b ], its (asymptotic) input-output behavior for each allowable value of λ is described in the Laplace domain by the following transfer function: T ( s ; λ ) = C z (cid:16) I n s − R X r =0 ( H r + λG r ) e − sτ r (cid:17) − B w ,and the associated robust H ∞ -norm is equal to k T ( · ; · ) k [ a, b ] H ∞ = max λ ∈ [ a, b ] k T ( · ; λ ) k H ∞ = max λ ∈ [ a, b ] ω ∈ R + σ (cid:0) T ( ω ; λ ) (cid:1) . (9)In [1] a novel numerical algorithm to compute the robust H ∞ -norm of anuncertain time-delay system is presented. This algorithm is based on the relationbetween the robust H ∞ -norm and the robust stability radius of an “uncertain”characteristic matrix. This relation is illustrated in Section 3.1. The resultingalgorithm is given in Section 3.2. Consider the following “uncertain” characteristic matrix M ( s ; λ, ∆) := I n s − R P r =0 ( H r + λG r ) e − sτ r − B w ∆ C z , (10)with H r , G r , B w , C z , τ r and λ as defined above, I n the identity matrix of size n and ∆ ∈ C m × p a complex-valued uncertainty with k ∆ k ≤ ε and ε ≥ λ which is real-valued and bounded to the interval [ a, b ] and a m × p matrix ∆which is complex-valued and bounded in spectral norm by ε , or in other words∆ ∈ B C m × p k·k ≤ ε with B C m × p k·k ≤ ε := (cid:8) ∆ ∈ C m × p : k ∆ k ≤ ε (cid:9) . Next, we define three important concepts related to this uncertain charac-teristic matrix. The spectral value set of this uncertain characteristic matrix isdefined as Λ [ a, b ] ε := [ λ ∈ [ a, b ] [ ∆ ∈B C m × p k·k ≤ ε { s ∈ C : det ( M ( s ; λ, ∆)) = 0 } . Note that this spectral value set is symmetric with respect to the imaginaryaxis. The pseudo-spectral abscissa is defined as the real part of the right-mostpoint, i.e. the point with the largest real part, in this spectral value set, α [ a, b ] ε := max n ℜ ( s ) : s ∈ Λ [ a, b ] ε o .Finally, the robust stability radius is defined as the smallest ε for which thisthis pseudo-spectral abscissa becomes non-negative, r [ a, b ] := min n ε ∈ R ≥ : α [ a, b ] ε ≥ o . (11)8urthermore, because α [ a, b ] ε is a continuous function of ε (this can be shownusing a similar argument as in [2, Section IV]), the transition to a non-negativepseudo-spectral abscissa is characterized by an ε for which α [ a, b ] ε equals zero.This means that the robust stability radius can also be defined as the smallest ε for which the spectral value set touches the imaginary axis: r [ a, b ] = min n ε ∈ R ≥ : ∃ ω ∈ R + such that ω ∈ Λ [ a, b ] ε o . (12)The following example illustrates these three concepts in more detail. Example 5.
Consider the following uncertain characteristic matrix: (cid:20) (cid:21) s − (cid:18)(cid:20) − − (cid:21) + λ (cid:20) − − (cid:21)(cid:19) − (cid:18)(cid:20) − −
10 2 (cid:21) + λ (cid:20) − (cid:21)(cid:19) e − s − (cid:20) − (cid:21) ∆ (cid:2) (cid:3) . (13) Figure 3 shows the part of the spectral value set in the region [ − . , − . × [2 . , .
7] for [ a, b ] equal to [ − ,
1] and several values of ε . For ε = 0, onlythe real-valued uncertainty λ plays a role and the spectral value set is a curve.For nonzero ε , also the complex-valued uncertainty ∆ affects the characteristicmatrix and the spectral value set becomes a region in the complex plane whichgrows as ε increases. Figure 4 shows the pseudo-spectral abscissa α [ a, b ] ε infunction of ε . We find that r [ − , = 0 . − , . × [ − , ω = 0). -0.5 -0.45 -0.4 -0.35 -0.3 -0.252.12.22.32.42.52.62.7 Figure 3: The part of the spectral value set of (13) in the region [ − . , − . × [2 . , .
7] for [ a, b ] equal to [ − ,
1] and ε equal to 0 (full line), 0 .
05 (dashedline), 0 . / H ∞ -norm asso-ciated with uncertain system (8) and the robust stability radius of (10). Theorem 6.
If uncertain system (8) is internally exponentially stable for all λ ∈ [ a, b ] , its associated robust H ∞ -norm is equal to the reciprocal of the robuststability radius of (10) .Proof. The following proof is a simplification of the result in [1].A complex number ω lies in Λ [ a, b ] ε if and only if there exist λ ∈ [ a, b ] and∆ ∈ B C m × p k·k ≤ ε such thatdet (cid:0) M ( ω ; λ, ∆) (cid:1) = det (cid:16) I n ω − P Rr =0 ( H r + λ G r ) e − ωτ r − B w ∆ C z (cid:17) = 0.9 Figure 4: The pseudo-spectral ab-scissa of (13) for [ a, b ] equal to[ − ,
1] in function of ε . Figure 5: The part of the spec-tral value set of (13) in [ − , . × [ − ,
10] for [ a, b ] equal to [ − , ε = 0 . λ ∈ [ a, b ],this is equivalent withdet (cid:16) I n − (cid:0) I n ω − P Rr =0 ( H r + λG r ) e − ωτ r (cid:1) − B w ∆ C z (cid:17) = 0.By the Weinstein-Aronszajn identity, this last equality can be rewritten asdet (cid:16) I − C z (cid:0) Iω − P Rr =0 ( H r + λG r ) e − ωτ r (cid:1) − B w ∆ (cid:17) = det (cid:0) I − T ( ω ; λ )∆ (cid:1) = 0 . The characterization of the robust stability radius in (12) can thus be rewrittenas r [ a, b ] = min ω ∈ R + min λ ∈ [ a, b ] min ∆ ∈ C m × p n k ∆ k : det (cid:0) I − T ( ω ; λ )∆ (cid:1) = 0 o .Using min ∆ ∈ C m × p {k ∆ k : det ( I − M ∆) = 0 } = σ ( M ) − from [15] one findsthat r [ a, b ] = min λ ∈ [ a, b ] ω ∈ R + (cid:16) σ (cid:0) T ( ω ; λ ) (cid:1)(cid:17) − = (cid:16) max λ ∈ [ a, b ] ω ∈ R + σ (cid:0) T ( ω ; λ ) (cid:1)(cid:17) − = (cid:16) k T ( · ; · ) k [ a, b ] H ∞ (cid:17) − , which concludes the proof. Remark 7.
The presented relation can be generalized to systems with uncer-tainties on the delays, multiple uncertainties and other uncertainty structures,such as full block and diagonal uncertainties. Also systems with delays and un-certainties in the input, output and direct feed-through terms can be considered.For more information, see [1].
This subsection presents a numerical algorithm to compute the robust stabilityradius. Once this quantity is found, the robust H ∞ -norm associated with (8)follows immediately from Theorem 6. 10y (11), the robust stability radius is the zero-crossing of the function R + ∋ ε α [ a, b ] ε . This zero-crossing can be found using the Newton-Bisection method,see [16, Chapter 9.4] for a reference implementation. This root finding methodrequires the evaluation of both α [ a, b ] ε and its derivative with respect to ε forgiven ε (whenever this derivative exists). The quantity α [ a, b ] ε can be computedusing the method presented in [1], which notes that α [ a, b ] ε is the solution of thefollowing optimization problem:max s, λ, ∆ ℜ ( s ) , subject to det (cid:0) M ( s ; λ, ∆) (cid:1) = 0 ,λ ∈ [ a, b ] , ∆ ∈ B C m × p k·k ≤ ε . (14)Furthermore, the following proposition shows that there exists a ∆ of rank oneand norm ε associated with local optima of this optimization problem. Thisresult will allow us to reduce the search space for ∆ to the space of matrices ofrank one and spectral norm ε . We will denote this space as B C m × p k·k = ε, rank=1 . Lemma 8.
Let s ⋆ Λ [ a, b ]0 be a local right-most point of Λ [ a, b ] ε , then there exist λ ⋆ ∈ [ a, b ] and ∆ ⋆ = εuv H with u ∈ C m , v ∈ C p and k u k = k v k = 1 such that det ( M ( s ⋆ ; λ ⋆ , ∆ ⋆ )) = 0 .Proof. Firstly, using similar ideas as in the proof of Theorem 6 one can showthat Λ [ a, b ] ε = Λ [ a, b ]0 ∪ n s ∈ C \ Λ [ a, b ]0 : max λ ∈ [ a, b ] σ ( T ( s ; λ )) ≥ ε − o .Because s ⋆ Λ [ a, b ]0 and s ⋆ is a right-most point of of Λ [ a, b ] ε , s ⋆ must lie on theboundary of (cid:8) s ∈ C \ Λ [ a, b ]0 : max λ ∈ [ a, b ] σ ( T ( s ; λ )) ≥ ε − (cid:9) . Hence, it holdsthat there exists a λ ⋆ ∈ [ a, b ] such that σ ( T ( s ⋆ ; λ ⋆ )) = ε − .Secondly, it can easily be verified that det (cid:0) I − T ( s ⋆ ; λ ⋆ )( ενυ H ) (cid:1) = 0 for υ and ν the left and right normalized singular vectors associated with σ (cid:0) T ( s ⋆ ; λ ⋆ ) (cid:1) .Following the derivation in the first part of the proof of Theorem 6 in theopposite direction one finds det (cid:0) M (cid:0) s ⋆ ; λ ⋆ , ενυ H (cid:1)(cid:1) = 0.Constrained optimization problem (14) is solved using a projected gradientflow method. The idea of this approach is to define a flow in the space ofpermissible variables along which the objective function monotonically increasesand whose attractive stationary points are (local) optimizers of the optimizationproblem. These stationary points can be found by choosing initial parametersand discretizing the resulting path till convergence to a stationary point usingEuler’s forward method. The step size is chosen such that the objective functionmonotonically increases along the discretized path. For more details on theusage of these methods for the computation of extreme points of spectral valuesets we refer to [2],[5] and [6].In our case we are thus looking for a path [0 , + ∞ ) ∋ θ ( λ ( θ ) , ∆( θ )) ∈ [ a, b ] × B C m × p k·k = ε, rank=1 such that the function θ max {ℜ ( s ) : det (cid:0) M (cid:0) s ; λ ( θ ) , ∆( θ ) (cid:1)(cid:1) = 0 }
11s monotonically increasing and such that the (local) optimizers of (14) appearas attractive stationary points. Furthermore, to improve computational perfor-mance we employ an explicit decomposition of ∆( θ ) as εu ( θ ) v ( θ ) H inspired by [5]. Here we consider the following flow where ˙ u , ˙ v and ˙ λ denote the derivatives of u , v and λ with respect to θ and where thedependency of u ( θ ), v ( θ ) and λ ( θ ) on θ is omitted for notational convenience. ˙ u = εξ ( θ ) (cid:16)(cid:0) I − uu H (cid:1) B Tw ϕ ( θ ) ψ ( θ ) H C Tz v + ℑ (cid:0) u H B Tw ϕ ( θ ) ψ ( θ ) H C Tz v (cid:1) u (cid:17) ˙ v = εξ ( θ ) (cid:16)(cid:0) I − vv H (cid:1) C z ψ ( θ ) ϕ ( θ ) H B w u + ℑ (cid:0) v H C z ψ ( θ ) ϕ ( θ ) H B w u (cid:1) v (cid:17) ˙ λ = λ = b and (cid:16)P Rr =0 ℜ (cid:0) ϕ ( θ ) H G r ψ ( θ ) e − s ( θ ) τ k (cid:1)(cid:17) > λ = a and (cid:16)P Rr =0 ℜ (cid:0) ϕ ( θ ) H G r ψ ( θ ) e − s ( θ ) τ k (cid:1)(cid:17) < ξ ( θ ) R P r =0 ℜ (cid:0) ϕ ( θ ) H G r ψ ( θ ) e − s ( θ ) τ k (cid:1) otherwise (15)with s ( θ ) the right-most characteristic root of M (cid:0) s ; λ ( θ ) , εu ( θ ) v ( θ ) H (cid:1) , and ϕ ( θ )and ψ ( θ ) the associated left and right eigenvectors, normalized such that ξ ( θ ) = ϕ ( θ ) H (cid:16) I + P Rr =0 ( H r + λ ( θ ) G r ) τ r e − s ( θ ) τ r (cid:17) ψ ( θ ) > . These paths are a combination of the paths presented in [2] (computingthe pseudo-spectral abscissa for real-valued Frobenius norm bounded uncer-tainties) and [5] (computing the pseudo-spectral abscissa for complex-valuedspectral norm bounded uncertainties). The fulfillment of the constraints, themonotonous increase of the objective function and the (local) optimality of sta-tionary points follows from these works. To conclude, we give some intuitionbehind (15): the right-hand side can be interpreted as a projection of the gra-dient of the objective function on the tangent space of the feasible set. Thisprojection step is needed to ensure that the variables remain within the feasibleset.The algorithm for numerically solving (14) is summarized in Algorithm 1,where s R (cid:0) M ( s ; λ, εuv H ) (cid:1) gives the right-most characteristic root of M ( s ; λ, εuv H ) and ˙ u k , ˙ v k and ˙ λ k correspond to the right-hand side of (15)evaluated at u k , v k and λ k . Remark 9.
To compute the right-most characteristic root, we use the algorithmpresented in [17]. This algorithm exploits the relation between a non-linear de-lay eigenvalue problem and the linear eigenvalue problem corresponding to thesolution operator of the associated delay differential equation. More precisely,this method uses a spectral discretization of the solution operator to approxi-mate the characteristic roots. These roots are subsequently refined by applyingNewton corrections based on the original non-linear eigenvalue problem formu-lation. This methods is however restricted to small problems. For large sparsematrices one could use iterative methods such as [10] and [8] to compute theright-most characteristic roots.Once α [ a, b ] ε is computed, its derivative with respect to ε can be computedalmost everywhere as shown in the following proposition.12 lgorithm 1 Discretization algorithm to solve (14). k ← λ , u , v s − ← −∞ and s ← s R (cid:0) M (cid:0) s ; λ , εu v H (cid:1)(cid:1) while | s k − s k − | > tol · | s k + s k − | do Find h such that ℜ (cid:16) s R (cid:16) M (cid:0) s ; λ k + h ˙ λ k , ε ( u k + h ˙ u k )( v k + h ˙ v k ) H (cid:1)(cid:17) (cid:17) ≥ℜ (cid:16) s R (cid:0) M ( s ; λ k , εu k v Hk ) (cid:1)(cid:17) if No h > tol h is found then stop. else λ k +1 ← λ k + h ˙ λ k ; λ k +1 ← max { a, min { λ k +1 , b }} ; u k +1 ← u k + h ˙ u k ; u k +1 ← u k +1 k u k +1 k ; v k +1 ← v k + h ˙ v k ; v k +1 ← v k +1 k v k +1 k ; s k +1 ← s R (cid:0) M (cid:0) s ; λ k +1 , εu k +1 v Hk +1 (cid:1)(cid:1) ; k ← k + 1; end ifend whileProposition 10. Let λ ⋆ and ∆ ⋆ = εu ⋆ v ⋆H be the unique optimizers of (14) and assume that the right-most characteristic root of M ( s ; λ ⋆ , ∆ ⋆ ) is simple,then dα [ a, b ] ε dε = ℜ (cid:16) ϕ ⋆H B w u ⋆ v ⋆H C z ψ ⋆ (cid:17) ϕ ⋆H ( I + P Rr =0 ( A k + λ ⋆ G r ) τ r e − τ r s ⋆ ) ψ ⋆ ,with s ⋆ the right-most characteristic root of M ( s ; λ ⋆ , ∆ ⋆ ) and φ ⋆ and ψ ⋆ itscorresponding left and right eigenvectors, normalized such that the denominatoris real and positive.Proof. Similar as in [2, Theorem 2]. H ∞ -controller synthesis method In this section we will describe a controller synthesis method for networkedsystems of form (3) whose associated adjacency matrix fulfills Assumption 1.The idea behind the presented method is to find a suitable controller by directlyminimizing (7) in the controller parameters p . Or in other words, we look forcontroller parameters p ⋆ that fulfill p ⋆ ∈ arg min p k ˆ T ˆ w ˆ z ( · ; p , · ) k [ a, b ] H ∞ . (16)Note however that k ˆ T ˆ w ˆ z ( · ; p , · ) k [ a, b ] H ∞ is only an upper bound for the actual H ∞ -norm of (4). The resulting control parameters will therefore in most cases notminimize the actual H ∞ -norm of the system, but this methodology has theadvantage that its computation cost does not depend on the number of subsys-tems. Furthermore, if Assumption 1 remains valid with a and b independentof the number of subsystems, then the resulting controller guarantees a level ofdisturbance rejection even if the number of subsystems is unknown.The minimization of (16) is however not trivial, as the robust H ∞ -norm maybe a non-smooth and non-convex function of the controller parameters even if13he controller matrices are analytic functions of the controller parameters. Thisprecludes the usage of standard optimization methods. Therefore, we will useHANSO [13], which implements a combination of BFGS with weak Wolfe linesearch and gradient sampling. Furthermore, to decrease the chance of endingup at a local optimum we will restart the optimization algorithm from severalinitial points. The optimization procedure requires the evaluation of both theobjective function and its derivative with respect to the control parameterswhenever this derivative exists. To evaluate the objective function, we usethe method presented in Section 3. The derivative with respect to the controlparameters follows from the following proposition. Proposition 11.
If there exists a unique pair ( ω ⋆ , λ ⋆ ) ∈ argmax λ ∈ [ a, b ] , ω ∈ R + σ ( T ( ω ; p , λ )) and if σ ( T ( ω ⋆ ; p , λ ⋆ )) is simple, then the function p
7→ k ˆ T ˆ w ˆ z ( · ; p , · ) k [ a, b ] H ∞ isdifferentiable at p , with dd p k T ( · ; p , · ) k [ a, b ] H ∞ = ℜ (cid:18) u H (cid:18) ∂∂ p T ( ω ⋆ ; p , λ ⋆ ) (cid:19) v (cid:19) ,in which u and v are the normalized left and right singular vectors associatedwith σ (cid:0) T ( ω ⋆ ; p , λ ⋆ ) (cid:1) , respectively. Finally, to start the optimization process we need initial controller parame-ters p for which the system is internally exponentially stable for all λ ∈ [ a, b ].To find such a starting point, we use the method presented in [3]. In this example we consider a networked system that consists of N frictionlesscarts that are interconnected using identical springs and that each balance aninverted pendulum. Furthermore, the first and the last cart are connected to thewall using additional (but identical) springs. This set-up is illustrated in Figure 6for N equal to 3. The dynamics of an individual cart-pendulum subsystem aregoverned by the following non-linear delay differential equations ( M + m )¨ x j ( t ) + ml cos (cid:0) θ j ( t ) (cid:1) ¨ θ j ( t ) − ml sin (cid:0) θ j ( t ) (cid:1)(cid:0) ˙ θ j ( t ) (cid:1) + k (cid:0) x j ( t ) − x j +1 ( t ) (cid:1) + k (cid:0) x j ( t ) − x j − ( t ) (cid:1) − u j ( t − τ u ) − w j, ( t ) = 0 l ¨ θ j ( t ) − g sin (cid:0) θ j ( t ) (cid:1) + ¨ x j ( t ) cos (cid:0) θ j ( t ) (cid:1) − w j, ( t ) = 0for j = 1 , . . . , N and M the mass of the individual carts, m the mass of thependulum’s bob which is connected to the cart using a massless rod of length l , k the spring constant, u j a controllable force that acts on the cart with aninput delay τ u , w j, and w j, external disturbances, θ j the angular displace-ment of the pendulum of cart j , x j the horizontal displacement of the j th cart’scenter with respect to its equilibrium position and x = x N +1 = 0. By a lin-earization around the equilibrium point ( x j , ˙ x j , θ j , ˙ θ j ) = (0 , , ,
0) and choosing14 x j ( t ) θ j ( t ) (cid:3) T as both measured and performance output, we obtain the fol-lowing linear state-space model of form (1):˙ x j ( t ) = − kM − mgM
00 0 0 1 kMl gl + mgMl x j ( t ) + M − Ml u j ( t − τ u ) + kM − kMl u nj ( t )+ M − mM − Ml l + mMl w j ( t ) (17) y j ( t ) = (cid:20) (cid:21) x j ( t ) , y nj ( t ) = (cid:2) (cid:3) x j ( t ) , z j ( t ) = (cid:20) (cid:21) x j ( t )with, u nj ( t ) = N X i =1 P Nj,i y ni ( t ),and P N = P line N = . . .
5. . . . . . . . .0 . . . . The input signal u j is generated using a controller of form (2) with n c equalto 2. Furthermore, we will allow communication between the controllers ofneighboring subsystems and as control parameters p we choose the elements ofthe control matrices. k M m l θ x u k M m l θ x u k M m l θ x u k Figure 6: Schematic representation of the considered set-up for N = 3.As mentioned in Theorem 4, the eigenvalues of P line N are restricted to theinterval [ − ,
1] for all
N >
1. It thus follows from Theorem 3 that the robust H ∞ -norm of a single subsystem forms an upper bound for the H ∞ -norm of theoverall network that holds independent of N . Furthermore, this upper boundcan be computed efficiently using the method presented in Section 3.15or the following parameter values: M = 1 kg, m = 0 .
05 kg, k = 1 N / m, l = 1 m, g = 9 . / s , τ u = 0 . τ u nc = 0 . J p ⋆ = " − . − . − . − . F p ⋆ = " − . . F n p ⋆ = " − . − . − . − . L p ⋆ = h − . − . i K p ⋆ = h . i K n p ⋆ = h − . − . i , (18)in which the subscript p ⋆ is used to indicate that these matrices correspondto the minimizing control parameters. The corresponding upper bound for the H ∞ -norm equals 0.512249 (rounded to 6 digits). Table 1 gives the actual H ∞ -norm of the closed loop system with controller matrices (18) for several N . Weobserve that these values lie close to the computed upper bound. This tablealso gives the time required by the algorithm described in [7] to compute thesevalues. As expected, the computation time grows roughly cubically with respectto N . For comparison, the computation time of the robust H ∞ -norm of theassociated uncertain subsystem equals 75.70 s. For system with large numberof subsystems it is thus beneficiary to minimize this upper bound instead of theactual H ∞ -norm.Table 1: The H ∞ -norm (rounded to 6 digits accuracy) of the closed loop net-worked system for several N and the computation time required by the algo-rithm described in [7] to compute these values. N k T ( · ; p ⋆ , N ) k H ∞ N equal to 20. We simulate the system for t ∈ [0 ,
10] where both the state andthe controller state are equal to 0 for t <
0, and each performance input is lowpass filtered ( ω cutoff = 6 π ) Gaussian white noise which is scaled after filteringto have a root mean square energy of 0.1. Figure 7 shows the performanceinputs and outputs of the 10 th subsystem. We observe that the noise is wellattenuated. The root mean square energy of z , and z , are equal to 0.03610and 0.02094, respectively.As explained in Theorem 3, the robust H ∞ -norm of an associated uncertainsubsystem of form (6) was used as upper bound for the actual H ∞ -norm ofthe networked system. Figure 8 shows the “worst-case gain function” of thisassociated uncertain subsystem: ω max λ ∈ [ − , σ (cid:16) ˆ T ˆ w ˆ z ( ω ; p ⋆ , λ ) (cid:17) (19)for ω ∈ [10 − , ]. The robust H ∞ -norm, which is equal to the maximal valueof this worst-case gain function, is indicated with a dashed line. We observe that16 t (s) -0.5-0.4-0.3-0.2-0.100.10.20.3 t (s) -0.15-0.1-0.0500.050.1 Figure 7: Simulation of the performance inputs (left) and outputs (right) of the10 th subsystem for t ∈ [0 ,
10] of the closed loop system (17)-(18) for N = 20.The state and the controller state are equal to zero for t <
0. Each performanceinput is low pass filtered ( ω cutoff = 6 π ) Gaussian white noise, scaled afterfiltering to have a root mean square energy of 0.1.the worst-case gain function is flat and almost equal to the robust H ∞ -norm for ω ∈ [0 , -1 Figure 8: The worst-case gain function of the associated uncertain subsystem asdefined in (19) for ω ∈ [10 − , ]. The robust H ∞ -norm, k ˆ T ˆ w ˆ z ( · ; p ⋆ , · ) k [ − , H ∞ ,is indicated with a dashed line. In this manuscript we introduced a scalable controller synthesis method fornetworked systems with identical subsystems and an interconnection topologythat fulfills Assumption 1. Using the decoupling transformation of Section 2 wearrived at an upper bound for the H ∞ -norm which can be computed at a costthat does not dependent on the number of subsystems. This upper bound was17he robust H ∞ -norm of a single subsystem with an additional scalar uncertainty.Subsequently an algorithm to compute this robust H ∞ -norm was discussed inSection 3. Finally, Section 4 showed how a controller that minimizes this upperbound can be synthesized using the direct optimization framework.To conclude, we note that the presented method can be extended to moregeneral identical subsystems, e.g. subsystems with direct feed-through terms.For such systems, however, the H ∞ -norm is sensitive to infinitesimal delaychanges and the strong H ∞ -norm is a more appropriate performance measure[7]. Furthermore, as in Section 2, one can show that the strong H ∞ -norm ofthe overall system is upper bounded by the robust strong H ∞ -norm of a singleuncertain subsystem. This robust strong H ∞ -norm can be computed using themethod presented in [1]. Acknowledgment
This work was supported by the project C14/17/072 of the KU Leuven ResearchCouncil and by the project G0A5317N of the Research Foundation-Flanders(FWO - Vlaanderen).
References [1] Pieter Appeltans and Wim Michiels. A pseudo-spectra based characterisa-tion of the robust strong H-infinity norm of time-delay systems with real-valued and structured uncertainties.
Preprint online available on arXiv, arXiv:1909.07778 [math.NA], 2019.[2] Francesco Borgioli and Wim Michiels. A Novel Method to Computethe Structured Distance to Instability for Combined Uncertainties on De-lays and System Matrices.
IEEE Transactions on Automatic Control ,65(4):1747–1754, 2020.[3] Deesh Dileep, Francesco Borgioli, Laurentiu Hetel, Jean Pierre Richard,and Wim Michiels. A scalable design method for stabilising decentralisedcontrollers for networks of delay-coupled systems.
IFAC-PapersOnLine ,51(33):68–73, 2018.[4] Deesh Dileep, Ruben Van Parys, Goele Pipeleers, Laurentiu Hetel,Jean Pierre Richard, and Wim Michiels. Design of robust decentralisedcontrollers for MIMO plants with delays through network structure ex-ploitation.
International Journal of Control , 2018.[5] Nicola Guglielmi and Christian Lubich. Differential equations for roamingpseudospectra: paths to extremal points and boundary tracking.
SIAMJournal on Numerical Analysis , 49(3):1194–1209, 2011.[6] Nicola Guglielmi and Christian Lubich. Low-rank dynamics for computingextremal points of real pseudospectra.
SIAM Journal on Matrix Analysisand Applications , 34(1):40–66, 2013.187] Suat Gumussoy and Wim Michiels. Fixed-order H-Infinity control for in-terconnected systems using delay differential algebraic equations.
SIAMJournal on Control and Optimization , 49(5):2212–2238, 2011.[8] Stefan G¨uttel, Roel Van Beeumen, Karl Meerbergen, and Wim Michiels.NLEIGS: A Class of Fully Rational Krylov Methods for Nonlinear Eigen-value Problems.
SIAM Journal on Scientific Computing , 36(6):A2842–A2864, 2014.[9] Gijs Hilhorst, Goele Pipeleers, Wim Michiels, and Jan Swevers. SufficientLMI conditions for reduced-order multi-objective H 2 / H ∞ control of LTIsystems. European Journal of Control , 23:17–25, 2015.[10] Elias Jarlebring, Karl Meerbergen, and Wim Michiels. A Krylov Methodfor the Delay Eigenvalue Problem.
SIAM Journal on Scientific Computing ,32(6):3278–3300, 2010.[11] Paolo Massioni and Michel Verhaegen. Distributed control for identical dy-namically coupled systems: A decomposition approach.
IEEE Transactionson Automatic Control , 54(1):124–135, 2009.[12] Wim Michiels. Spectrum-based stability analysis and stabilisation of sys-tems described by delay differential algebraic equations.
IET Control The-ory & Applications , 5(16):1829–1842, 2011.[13] Michael L. Overton. HANSO: a hybrid algorithm for nonsmooth optimiza-tion.
Available online: https://cs.nyu.edu/overton/software/hanso/ ,2009.[14] S. Mert ¨Ozer and Altu Iftar. Decentralized controller design for time-delaysystems by optimization.
IFAC-PapersOnLine , 28(12):462–467, 2015.[15] Andrew Packard and John C. Doyle. The complex structured singularvalue.
Automatica , 29(1):71–109, 1993.[16] William H. Press, Saul A. Teukolsky, and William T. Vetterling.
Numer-ical recipes in Fortran 77 : the art of scientific computing . CambridgeUniversity press, Cambridge, 2nd edition, 1996.[17] Zhen Wu and Wim Michiels. Reliably computing all characteristic rootsof delay differential equations in a given right half plane using a spectralmethod.
Journal of Computational and Applied Mathematics , 236(9):2499–2514, 2012.[18] Kemin Zhou and John C. Doyle.