Decomposition and Projection Methods for Distributed Robustness Analysis of Interconnected Uncertain Systems

(1)

Technical report from Automatic Control at Linköpings universitet

Decomposition and Projection Methods

for Distributed Robustness Analysis of

Interconnected Uncertain Systems

Sina Khoshfetrat Pakazad, Martin S. Andersen,

Anders Hansson, Anders Rantzer

Division of Automatic Control

E-mail: sina.kh.pa@isy.liu.se, martin.andersen@isy.liu.se,

hansson@isy.liu.se, anders.rantzer@control.lth.se

2nd December 2011

Report no.: LiTH-ISY-R-3033

Submitted to Technical Report

Address:

Department of Electrical Engineering Linköpings universitet

SE-581 83 Linköping, Sweden

WWW: http://www.control.isy.liu.se

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

Technical reports from the Automatic Control group in Linköping are available from http://www.control.isy.liu.se/publications.

(2)

We consider a class of convex feasibility problems where the constraints

thatdescribethefeasiblesetareloosely coupled. These problems arise in

robuststabilityanalysisoflarge, weaklyinterconnectedsystems. To

facili-tatedistributedimplementationofrobuststabilityanalysisofsuchsystems,

weproposetwoalgorithmsbasedon decompositionandsimultaneous

pro-jections. The rst algorithm is a nonlinear variant of Cimmino's mean

projection algorithm, but by taking the structure of the constraints into

account,wecanobtainafasterrateofconvergence. Thesecondalgorithm

is devisedby applyingthealternatingdirection methodof multiplierstoa

convex minimization reformulation of theconvex feasibility problem. We

usenumericalresultstoshowthatbothalgorithmsrequirefarlessiterations

thantheacceleratednonlinearCimminoalgorithm.

(3)

Distributed Robustness Analysis of

Interconnected Uncertain Systems

Sina Khoshfetrat Pakazad, Martin S. Andersen, Anders Hansson

and Anders Rantzer

2011-12-02

Abstract

Weconsideraclassofconvexfeasibilityproblemswheretheconstraints

thatdescribe the feasibleset are loosely coupled. These problems arise

inrobuststability analysis of large,weaklyinterconnected systems. To

facilitate distributedimplementationof robuststability analysis of such

systems,weproposetwoalgorithmsbasedondecompositionand

simulta-neousprojections. TherstalgorithmisanonlinearvariantofCimmino's

meanprojectionalgorithm,butbytakingthestructureoftheconstraints

into account, we can obtain a faster rate of convergence. The second

algorithmisdevisedbyapplyingthealternatingdirectionmethodof

mul-tipliers toa convexminimization reformulation ofthe convexfeasibility

problem. Weusenumericalresultstoshowthatbothalgorithmsrequire

(4)

(5)

1

Decomposition and Projection Methods for

Distributed Robustness Analysis of

Interconnected Uncertain Systems

Sina Khoshfetrat Pakazad, Martin S. Andersen, Anders Hansson, Anders Rantzer,

Abstract—We consider a class of convex feasibility problems

where the constraints that describe the feasible set are loosely coupled. These problems arise in robust stability analysis of large, weakly interconnected systems. To facilitate distributed implementation of robust stability analysis of such systems, we propose two algorithms based on decomposition and simulta-neous projections. The first algorithm is a nonlinear variant of Cimmino’s mean projection algorithm, but by taking the structure of the constraints into account, we can obtain a faster rate of convergence. The second algorithm is devised by applying the alternating direction method of multipliers to a convex minimization reformulation of the convex feasibility problem. We use numerical results to show that both algorithms require far less iterations than the accelerated nonlinear Cimmino algorithm.

Index Terms—Robust stability analysis, convex feasibility

prob-lems, projection algorithms, decomposition, distributed comput-ing.

I. INTRODUCTION

D

ISTRIBUTED and parallel solutions for control prob-lems are of great importance. Prioritizing such ap-proaches over the centralized solutions is mainly motivated by structural or computational reasons, [1]–[5]. Such solutions are common in the control of distributed interconnected systems, which appear in many practical applications, [2], [5]–[10]. Considering the fact that many of the control design methods are model based, they are vulnerable to model uncertainty. Some of these design methods consider the model uncertainty, but many, specially the ones employed in industry, neglect it. As a result, for such methods it is important to address the stability of the closed loop system with respect to the model uncertainty.

Methods for robust stability analysis of uncertain systems has been studied thoroughly, e.g., see [11]–[14]. One of popular methods for robust stability analysis is the so called

µ-analysis, an approximation of which solves a Convex

Fea-sibility Problem, CFP, for each frequency on a finite grid of frequencies. For networks of interconnected uncertain systems, each CFP is a global feasibility problem that involves the entire network. In other words, each CFP depends on all the local models, i.e., the models associated with the subsystems in the network, as well as the network topology. In networks

S. K. Pakazad, M. S. Andersen, and A. Hansson are with the Division of Automatic Control, Department of Electrical Engineering, Link¨oping University, Sweden.

A. Rantzer is with The Department of Automatic Control, Lund University, Sweden.

where all the local models are available to a central unit and the size of the problem is not computationally prohibitive, each feasibility problem can be solved efficiently in a centralized manner. However distributed or parallel methods are desirable in networks, for instance, where the interconnected system di-mension is large or where the local models are either private or available to only a small number of subsystems in the network, [15]–[18]. In [15]–[17], for specific network structures, the system matrix is diagonalized and decoupled using decoupling transformations. This decomposes the analysis and design problem and provides the possibility to produce distributed and parallel solutions. An overview of such algorithms for analysis of systems governed by Partial Differential Equations, PDEs, is given in [19].

In this paper, we use decomposition methods to facili-tate distributed robust stability analysis of large networks of weakly interconnected uncertain systems. Decomposition allows us to reformulate a single global feasibility constraint

x ∈ C ⊂ Rn_involving_{m uncertain systems as a set of N < m} loosely coupled constraints

x ∈ Ci ≡ {x ∈ Rn| fi(x) K 0}, i = 1, . . . , N (1) wherefi(x) Ki0 is a linear conic inequality with respect to

a convex coneKi. We will assume thatfi depends only on a small subset of variables{xk| k ∈ Ji}, Ji ⊆ {1, . . . , n}. The number of constraintsN is less than the number of systems m

in the network, and hence the functionsfi generally involve problem data from more than one subsystem in the network. We will henceforth assume that eachfiis described in terms of only a small number of local models such that the Euclidean projection PCi(x) of x onto Ci involves just a small group

of systems in the network. This assumption is generally only valid if the network is weakly interconnected.

One algorithm that is suited for distributed solution of the CFP is the nonlinear Cimmino algorithm [20], [21] which is also known as the Mean Projection Algorithm. This algorithm is a fixed-point iteration which takes as the next iteratex(k+1) a convex combination of the projectionsPCi(x

(k)_{), i.e.,} x(k+1)_:= N X i=1 αiPCi(x (k)₎ ₍₂₎

where PN_i=1αi = 1 and α1, . . . , αN > 0. Notice that each iteration consists of two steps: a parallel projection step which is followed by a consensus step that can be solved by means of distributed weighted averaging, e.g., [22]–[26]. Iusem &

(6)

2

De Pierro [27] have proposed an accelerated variant of the nonlinear Cimmino algorithm that takes as the next iterate a convex combination of the projections of x(k) _{on only the} sets for which x(k) _{∈ C}_/

i. This generally improves the rate of convergence when only a few constraints are violated. However, the nonlinear Cimmino algorithm and its accelerated variant may take unnecessarily conservative steps when the sets Ci are loosely coupled. We will consider two algorithms that can exploit this type of structure, and both algorithms are closely related to the nonlinear Cimmino in that each iteration consists of a parallel projection step and a consensus step.

The first algorithm that we consider is equivalent to von Neumann Alternating Projection (AP) algorithm [28], [29] in a product space E = R|J1|_{× · · · × R}|JN| _{of dimension}

PN

i=1|Ji|. As a consequence, this algorithm converges at a linear rate under mild conditions, and its behavior is well-understood also when the CFP is infeasible [30]. Using the ideas from [31], [32], we also show how this algorithm can be implemented in a distributed manner.

A CFP can also be formulated as a convex minimization problem which can be solved with distributed optimization algorithms; see e.g. [23], [31]–[33]. The second algorithm that we consider is the Alternating Direction Method of Multipliers (ADMM) [31], [32], [34], [35], applied to a convex minimization formulation of the CFP. Unlike AP, ADMM also makes use of dual variables, and when applied to the CFP, it is equivalent to Dykstra’s alternating projection method [30], [36] in the product spaceE. Although there exist problems for

which Dykstra’s method is much slower than the classical AP algorithm, it generally outperforms the latter in practice. Outline

The paper is organized as follows. In Section II, we present a product space formulation of CFPs with loosely coupled sets, and we propose an algorithm based on von Neumann AP method. In Section III, we consider a convex minimization reformulation of the CFP, and we describe an algorithm based on the ADMM. We discuss distributed implementation of both algorithms in Section IV, and in Section V, we consider distributed solution of the robust stability analysis problem. We present numerical results in Section VI, and we conclude the paper in Section VII.

Notation

We denote with Npthe set of positive integers{1, 2, . . . , p}. Given a setJ ⊂ {1, 2, . . . , n}, the matrix EJ ∈ R|J|×nis the 0-1 matrix that is obtained from an identity matrix of order n by deleting the rows indexed by Nn\ J. This means that EJx is a vector with the components of x that correspond to the elements in J. The distance from a point x ∈ Rn _{to a set} C ⊆ Rn _{is denoted as} _{dist(x, C), and it is defined as}

dist(x, C) = inf

y∈Ckx − yk2. (3)

Similarly, the distance between two sets C1, C2 ⊆ Rn is defined as

dist(C1, C2) = inf y∈C1,x∈C2

kx − yk₂. (4)

The relative interior of a setC is denoted relint(C), and D = diag(a1, . . . , an) is a diagonal matrix of order n with diagonal entriesDii = ai. The column space of a matrixA is denoted C_(A).

II. DECOMPOSITION AND PROJECTION METHODS

A. Decomposition and convex feasibility in product space Given N loosely coupled sets C1, . . . , CN, as defined in (1), we defineN lower-dimensional sets

¯

Ci= {si∈ R|Ji|| EJTisi∈ Ci}, i = 1, . . . , N (5)

such thatsi∈ ¯Ci impliesEJTisi∈ Ci. With this notation, the

standard form CFP can be rewritten as find s1, s2, . . . , sN, x subject to si∈ ¯Ci, i = 1, . . . , N

si= EJix, i = 1, . . . , N

(6)

where the equality constraints are the so called coupling con-straints that ensure that the variabless1, . . . , sN are consistent with one another. In other words, if the setsCiandCj(i 6= j) depend on xk, then the kth component of EJTisi and E

T Jjsj

must be equal.

The formulation (6) decomposes the global variablex into N coupled variables s1, . . . , sN. This allows us to rewrite the problem as a CFP with two sets

find S subject to S ∈ C, S ∈ D (7) where S = (s1, . . . , sN) ∈ R|J1|× · · · × R|JN| C = ¯C1× · · · × ¯CN D = { ¯Ex | x ∈ Rn_} ¯ E =ET J1 · · · E T JN T .

The formulation (7) can be thought of as a “compressed” product space formulation of a CFP with the constraints in (1), and it is similar to the consensus optimization problems described in [31, Sec. 7.2].

B. Von Neumann’s alternating projection in product space The problem (7) can be solved using the von Neumann’s AP method. Given a CFP with two sets

find x

subject to x ∈ A, x ∈ B (8)

and a starting point y(0) _{= x}

0, von Neumann’s AP method computes two sequences [37]

z(k+1)= PA y(k) (9a) y(k+1)= PB z(k+1). (9b)

If the CFP is feasible, i.e.,A ∩ B 6= 0, then both sequences

converge to a point inA ∩ B. We discuss infeasible problems

(7)

3

the von Neumann AP algorithm to the CFP (7) results in the update rule S(k+1)= PC X(k) (10a) =P_C¯₁ EJ1x (k)_{, . . . , P} ¯ CN EJNx (k) X(k+1)= ¯E ¯ET_E_¯−1_E_¯T_S(k+1) | {z } x(k+1) , (10b)

where (10a) is the projection ontoC, and (10b) is the projection

onto the column space of ¯E. The projection onto the set C can

be computed in parallel by N computing agents, i.e., agent i computes s(k)_i = P_C¯_i(EJix

(k)_{), and the second projection} can be interpreted as a consensus step that can be solved via distributed averaging.

C. Convergence

Gubin et al. [38] showed that a sequence of cyclic projec-tions converges linearly to somex∗_{∈ C = ∩}

iCi6= ∅ if the set {C1, . . . , CN} is boundedly linearly regular, i.e., there exists θB > 0 for every bounded set B in Rn such that

dist(x, C) ≤ θB max

1≤i≤N{dist(x, Ci)}, ∀x ∈ B. (11) It was later shown by Bauschke et al. [39] and Beck & Teboulle [40] that Slater’s condition for the CFP implies bounded linear regularity, i.e., (11) holds and hence the cyclic projection algorithm converges linearly if

k \ i=1 Ci ! \ \N i=k+1 relint(Ci) ! 6= ∅ (12)

where C1, . . . , Ck are polyhedral sets andCk+1, . . . , CN are general closed, convex sets. For two sets, cyclic projections is equivalent to alternating projections, and when {A, B} is

boundedly linearly regular, it follows from [40, Thm. 2.2] that

dist(x(k+1), C) ≤ γBdist(x(k), C) with γB= s 1 − 1 θ2 B (13) where C = A ∩ B and θB a positive constant that depends of the starting point x(0) _{and that satisfies (11) with} _{B =} {x | kx − zk ≤ kx(0)_{− zk} for any z ∈ C.}

For infeasible problems (i.e.,A∩B = ∅), the iterates satisfy y(k)− z(k), y(k)− z(k+1)→ v, kvk2= dist(A, B). (14) Moreover, the sequencesy(k)_and_z(k)_{converge if the distance} dist(A, B) is attained, and otherwise ky(k)_k

2 → ∞ and

kz(k)_{k → ∞ [30], [37]. This means that we can use the} difference v(k) _{= y}(k) _{− z}(k) _{to detect infeasibility. If the} set C in the CFP in (7) is a closed set, it is possible to use

the statement in (14) for detecting the infeasibility of the CFP in (6). Note that this is the case for the problem described in Section V.

III. CONVEX MINIMIZATION REFORMULATION

The CFP (6) is equivalent to a convex minimization problem minimize PN_i=1gi(si)

subject to si= EJix, i = 1, . . . , N,

(15) with variablesS and x, where

gi(si) = (

∞ si6∈ ¯Ci 0 si∈ ¯Ci

(16) is the indicator function for the set ¯Ci. In the following, we apply the ADMM to the problem (15).

A. Solution via Alternating Direction Method of Multipliers Consider the following equality constrained convex opti-mization problem

minimize F (x)

subject to Ax = b (17)

with variable x ∈ Rn_{, and where} _{A ∈ R}m×n_,_{B ∈ R}m_{, and} F : Rn_{7→ R. The augmented Lagrangian for this problem can} be expressed as Lρ(x, λ) = F (x) + λT(Ax − b) + ρ 2kAx − bk 2 2, (18) where ρ > 0 and λ ∈ Rm_{. This can also be rewritten in} normalized form as ¯ Lρ(x, ¯λ) = F (x) + ρ 2 Ax − b + ¯λ 2₂−ρ 2 ¯λ 2₂ (19) where ¯Lρ(x, ¯λ) is referred to as the normalized augmented Lagrangian, and ¯λ = λ/ρ is the so-called normalized Lagrange

variable.

A general convex optimization problem of the form minimize F (S)

subject to S ∈ A (20)

can be cast as a problem of the form (17) by means of the indicator functiong(S) for the set A, i.e.,

minimize F (S) + g(q)

subject to S − q = 0 (21)

whereq is an auxiliary variable. The normalized augmented

Lagrangian for this problem is given by

¯ Lρ(S, q, ¯λ) = F (S) + g(q) +ρ 2 S − q − ¯λ 2 2− 1 2ρ ¯λ 2 2. (22) The problem in (21) can be solved by applying the ADMM to (22); see e.g. [31] and [32, Sec. 3.4]. This results in the following iterative algorithm

S(k+1)= argmin S F (S) +ρ 2 S − q(k)− ¯λ(k) 2 2 q(k+1)= argmin q g(q) +ρ 2 S(k+1)− q − ¯λ(k) 2 2 = PA(S(k+1)+ ¯λ(k)), ¯ λ(k+1)= ¯λ(k)+S(k+1)− q(k+1).

(8)

4

The problem (15) can be solved using a similar approach. The normalized augmented Lagrangian for this problem is given by Lρ(S, x, ¯λ) = N X i=1 n gi(si) + ρ 2 si− EJix − ¯λi 2 2− ρ 2 ¯λi 2 2 o (23) where ¯λi ∈ R|Ji| and ¯λ = ¯λ1, . . . , ¯λN. If we apply the ADMM to (23), we obtain the following iterative method

S(k+1)= argmin S (N X i=1 gi(si) + ρ 2 si− EJ_ix(k)₋¯λ(k)i 2 2 ) x(k+1)= argmin x (N X i=1 ρ 2 s (k+1) i − EJ_ix − ¯λ(k)i 2 2 ) ¯ λ(k+1)i = ¯λ (k) i + (s (k+1) i − EJ_ix(k+1)). (24)

The update ofS can be decomposed into N subproblems s(k+1)_i = argmin si gi(si) + ρ 2 si− EJix (k)_{− ¯}_λ(k) i 2 2 = P_C¯_i(EJix (k)_{− ¯}_λ(k) i ), (25) and the update ofx can be computed as

x(k+1)= ¯ETE¯−1E¯TS(k+1)+ ¯λ(k). (26) Using the definition X = ¯Ex, we can write the update of X

as

X(k+1)= ¯E ¯ETE¯−1E¯TS(k+1)+ ¯λ(k). (27) Furthermore, given the starting pointx(0) _{= x}

0 and ¯λ(0)= 0, we obtain ¯ λ(k)₌ k X i=1 S(i)_{− X}(i)_. ₍₂₈₎

The matrix ¯E ¯ET_E_¯−1_E_¯T _{in (27) defines an orthogonal} projection onto the column space of the matrix ¯E, and hence

¯

E ¯ET_E_¯−1_E_¯T¯_λ(k) _{= 0 for all k ≥ 0. This allows us to} simplify (27) as

X(k+1)= ¯E ¯ET_E_¯−1_E_¯T_S(k+1)_. ₍₂₉₎ To summarize, the update rules for the ADMM, applied to (15), are given by

S(k+1)= PC(X(k)− ¯λ(k)) (30a) X(k+1)= ¯E ¯ET_E_¯−1_E_¯T_S(k+1) _(30b)

¯

λ(k+1)= ¯λ(k)+ (S(k+1)− X(k+1)). (30c) We remark that this algorithm is a special case of the algorithm developed in [31, Sce. 7.2] for consensus optimization. We also note that this algorithm is equivalent to Dykstra’s alter-nating projection method for two sets where one of the sets is affine [30]. Moreover, it is possible to detect infeasibility in the same way that we can detect infeasibility for the AP method, i.e., X(k)_{− S}(k)_{→ v where kvk}

2= dist(C, D) [30].

Algorithm 1 Alternating projection algorithm in product space 1: Given x(0) _{and for all} _{i = 1, . . . , N , each agent i should}

2: repeat

3: s(k+1)_i ← P_C¯_i

x(k)_J_i .

4: Communicate with all agentsr belonging to Ne(i).

5: for allj ∈ Ji do 6: x(k+1)_j =_|I1 j| P q∈Ij ET Jqs (k+1) q j. 7: end for 8: k ← k + 1 9: until forever

Note that, unlike the AP method, the iteration (30) does not necessarily converge with a linear rate when the feasible sets satisfy (11).

IV. DISTRIBUTED IMPLEMENTATION

In this section, we describe how the algorithms can be implemented in a distributed manner using techniques that are similar to those described in [31] and [32, Sec. 3.4]. Specifically, the parallel projection steps in the algorithms expressed by the update rules in (10) and (30) are amenable to distributed implementation. We will henceforth assume that a network of N computing agents is available.

LetIi = {k | i ∈ Jk} ⊆ Nn denote the set of constraints that depend onxi. Then it is easy to verify that

¯

ET_{E = diag(|I}_¯

1|, . . . , |In|) (31) and consequently, thejth component of x in the update rules

(10b) and (30b) is of the form

x(k+1)_j = 1 |Ij| X q∈Ij EJTqs (k+1) q j. (32)

In other words, the agents in the setIj must solve a distributed averaging problem to computex(k+1)_j . LetxJi = EJix. Since

the set ¯Ciassociated with agenti involves the variables index byJi, agenti should update xJi by performing the update in

(32) for allj ∈ Ji. This requires agenti to communicate with all agents in the set

Ne(i) =nj | Ji \

Jj6= ∅ o

, (33)

which we call the neighbors of agenti. Each agent j ∈ Ne(i)

shares one or more variables with agent i. Distributed

vari-ants of the algorithms presented in the previous sections are summarized in Algorithms 1 and 2.

A. Feasibility Detection

For strictly feasible problems, the algorithms 1 and 2 converge to a feasible solution. We now discuss different techniques for checking feasibility.

1) Global Feasibility Test: Perhaps the easiest way to check feasibility is to directly check the feasibility of x(k) _with respect to all the constraints. This can be accomplished by explicitly forming x(k)_{. This requires a centralized unit that} recieves the local variables from the individual agents.

(9)

5

Algorithm 2 ADMM based algorithm

1: Given x(0)_{, ¯}_λ(0) _{= 0, for all i = 1, . . . , N , each agent i} should

2: repeat

3: s(k+1)_i ← P_C¯_i

x(k)_J_i − ¯λ(k)_i .

4: Communicate with all agents r belonging to Ne(i).

5: for allj ∈ Ji do 6: x(k+1)_j = _|I1 j| P q∈Ij ET Jqs (k+1) q j. 7: end for 8: λ¯(k+1)_i ← ¯λ(k)_i + s(k+1)_i − x(k+1)_J_i 9: k ← k + 1 10: until forever

2) Local Feasibility Test: It is also possible to chech feasibility locally. Instead of sending the local variables to a central unit, each agent declares its feasibility status with respect to its local constraint. This type of feasibility detection method is based on the following Lemmas.

Lemma 1: If x(k)_J_i ∈ ¯Ci, for all i = 1, · · · , N then using these vectors a feasible solution,x, can be constructed for the

original problem, i.e.,x ∈TN_i=1Ci.

Proof: The update rule for each of the components of the variablex is merely the average of all the corresponding

local iterates from other agents. Hence, all iterates of each of the components of the x variable are equal. As a result,

based on the definition of the ¯Ci, by havingx(k)_J_i ∈ ¯Ci for all i = 1, · · · , N , it is possible to construct a feasible vector x(k) from the iterates x(k)_J_i .

Lemma 2: If kX(k+1) _{− X}(k)_k

2 = 0 and kS(k+1) − X(k+1)_k

2= 0, then a feasible solution, x, can be constructed for the original problem, i.e.,x ∈TN_i=1Ci.

Proof: Consider the update rules for Algorithm 1. The conditions stated in the lemma imply that s(k+1)_i = x(k)_J_i . As a result for alli = 1, · · · , N

s(k+1)_i = P_C¯_i

x(k)_J_i = x(k)_J_i , (34) which implies thatx(k)_J_i ∈ ¯Ci. Therefore, similar to Lemma 1, it is possible to generate a feasible solution out of vectorsx(k)_J_i . Consider the update rules in Algorithm 2 for the ADMM based algorithm. Then the conditions stated above imply that

S(k+1)= PC X(k)− ¯λ(k)= X(k), (35a) X(k+1)= ¯E ¯ET_E_¯−1_E_¯T_X(k+1)_{= X}(k) _(35b) ¯ λ(k+1)= ¯λ(k)+S(k+1)− X(k+1)= ¯λ(k). (35c) Hence any solution that satisfies the conditions in the lemma, constitutes an equilibrium point for the update rule in (30), and as a result

P_C¯_i

x(k)_J_i = x(k)_J_i , (36) which in turn provides the possibility to construct a feasible solution.

Remark 1: For Algorithm 1, conditions in lemmas 1 and 2 are equivalent. However, this is not the case for Algorithm 2. For this algorithm, satisfaction of the conditions in Lemma 2 implies the conditions in Lemma 1.

Remark 2: Feasibility detection using Lemma 1, requires additional computations for Algorithm 2. These computations include local feasibility check of the iterates xJi. Note that

this check does not incur any additional cost for Algorithm 1.

B. Infeasibility Detection

Recall from Section II that if the CFP is infeasible, then the sequencekS(k)_{− X}(k)_k

2will converge to a nonzero constant kvk2 = dist(C, D). Therefore, in practice, it is possible to detect infeasible problems by limiting the number of iterations, i.e., if kS(k) _{− X}(k)_k

2 is not sufficiently small after the maximum number of iterations, the problem is considered to be infeasible.

V. ROBUSTSTABILITYANALYSIS

Robust stability of uncertain large scale weakly intercon-nected systems with structured uncertainty can be analyzed through the so-called µ-analysis framework [41]. This leads

to CFP which is equivalent to a semidefinite programming (SDP) problem involving the system matrix.

Consider the following system description

Y (s) = M (s)U (s), (37)

whereM (s) is a m × m transfer function matrix, and let

U (s) = ∆Y (s), (38)

where∆ = diag(δi), with δi ∈ R, |δi| ≤ 1 for i = 1, · · · , m, represent the uncertainties in the system. The system is said to be robustly stable if there exists a diagonal positive definite

X(ω) and 0 < µ < 1 such that

M (jω)∗X(ω)M (jω) − µ2X(ω) ≺ 0, (39) for allω ∈ R. Note that this problem is infinite dimensional,

and in practice, it is often solved approximately by discretizing the frequency variable.

In the following, we investigate only a single frequency point such that the dependence on the frequency is dropped. Moreover, for the sake of simplicity, we assume thatM is

real-valued. The extension to complex valuedM is straightforward.

As a result, feasibility of the following CFP is a sufficient condition for robust stability of the system

find X

subject to M′_{XM − X −ǫI} xi≥ ǫ, i = 1, . . . , m

(40)

forǫ > 0 and where xi are the diagonal elements ofX. A large scale network of weakly interconnected uncertain systems can also be represented in the form (37)-(38). In the case of weakly interconnected system, the system matrix

(10)

6 1 x1 2 x2 3 x3 m − 2 xm₋₂ m − 1 xm₋₁ m xm · · · w1 y1 z1 w2 y2 z2 w1 y1 z1 w2 y2 z2 wm₋₃ ym₋₃ zm₋₃ wm₋₄ ym₋₄ zm₋₄ wm₋₃ ym₋₃ zm₋₃ wm₋₄ ym₋₄ zm₋₄

Fig. 1. A chain of m uncertain systems with problem variables x1, . . . , xm and auxiliary (coupling) variables wi, yi, zi, i = 1, . . . , m − 3. Each box

corresponds to an uncertain system, and the dashed lines separate the m − 2 groups of systems that are associated with the setsC¯1, . . . , ¯Cm₋₂.

M (s) ∆

Y (s) U (s)

Fig. 2. An uncertain system.

example, we consider a chain of systems which leads to a tri-diagonal system matrix

M =         g1 h1 0 0 0 f2 g2 h2 0 0 0 f3 . .. . .. 0 0 0 . .. gm−1 hm−1 0 0 0 fm gm         . (41)

This system matrix is obtained if the input-output relation for the underlying systems are given by

p1= g1q1+ h1z2 z1= q1 q1= δ1p1, (42) pm= gmqm+ fmzm−1 zm= qm qm= δ1pm, (43) and pi= giqi+ fizi−1+ hizi+1 zi= qi qi= δipi, (44)

fori = 2, . . . , m − 1. The tri-diagonal structure in the system

matrix implies that the LMI defined in (40) becomes banded. This is a special case of a so-called chordal sparsity pattern, and these have been exploited in SDP by several authors; see [42]–[45].

Positive semidefinite matrices with chordal sparsity patterns can be decomposed into a sum of matrices that are positive semidefinite [44], [46, Sec. 5.1]. For example, a positive semidefinite band matrix can be decomposed into a sum of positive semidefinite matrices, as illustrated in Figure 4. Note that the nonzero blocks marked in the matrices on the right-hand side in Figure 4 are structurally equivalent to the block

M1(s) M2(s) Mm(s) δ1 δ2 δm p1 q1 z1 z2 p2 q2 z2 z3 pm qm zm−1 zm · · ·

Fig. 3. An example of a weakly interconnected uncertain system.

= + · · · +

Fig. 4. A positive semidefinite band matrix of order n with half-bandwidth p can be decomposed as a sum of n − p positive semidefinite matrices.

in the left-hand side, but the numerical values are generally different. Using this decomposition technique, the problem in (40) can be reformulated as in (7).

By introducing auxiliary variables w, y, z ∈ Rm−3 _and lettingq = (x, w, y, z), the five-diagonal matrix M′_{XM − X} can be decomposed as f1(q) = _g 1 f2 h1g2 0 h2 x1 0 0 x2 g1 f2 h1 g2 0 h2 ′ − _x 1 0 0 0 x2−w1−y1 0 −y1 −z1 (45a) fm−1(q) = _h m−1 0 gm−1 fm hm−1gm h xm−1 0 0 xm i fm−1 0 gm−1 fm hm−1gm ′ (45b) − _w m−3 ym−3 0 ym−3 wm−3+xm−1 0 0 0 xm and fi(q) = _f i+1 gi+1 hi+1 [xi+1] _f i+1 gi+1 hi+1 ′ (45c) − _w i−1 yi−1 0

yi−1 zi−1+xi+1−wi−yi

0 −yi −zi

for i = 2, . . . , m − 3. Notice that (45a) and (45b) depend

on data from two subsystems whereas (45c) depends on data from only one subsystem. This dependence is also illustrated in Figure 1. The right-hand side of the LMI in (40) can be

(11)

7 0 5 10 15 20 25 0 5 10 15 20 25 30 35 40 45 50 No. Iterations

No. Violated Constraints

Fig. 5. Algorithm 1: Number of violated constraints in the CFP after each update of the global variables. This figure illustrates the results for 50 randomly generated problems.

decomposed in a similar manner

F1= h_ǫ_{0 0} 0 ǫ 0 0 0 0 i , Fm−2 = h_{0 0 0} 0 ǫ 0 0 0 ǫ i , Fi= h_{0 0 0} 0 ǫ 0 0 0 0 i , i = 2, . . . , m − 3.

With this decomposition, we can reformulate the LMI in (40) as a set of m − 2 LMIs

q ∈ Ci= {q | fi(q) Fi}, i = 1, . . . , m − 2, or equivalently, as the constraints

si∈ ¯Ci= {EJiq | fi(q) Fi} ⊆ R

|Ji|_{, s}

i= EJiq, (46)

fori = 1, . . . , m−2. Here Jiis the set of indices of the entries ofq that are required to evaluate fi. Notice that (46) is of the form (6). The CFP defined by the constraints in (46) can be solved over a network of m − 2 agents. The connectivity of

this network and the coupling variables among the agents are illustrated in Figure 1.

Remark 3: As was mentioned above, the CFP in (40) is often solved over a frequency grid. If these frequency points are chosen sufficiently close, it is probable that the CFPs for adjacent frequencies are similar. As a result, it is very likely that a solution to a CFP for one frequency is either a solution to or close to a solution to the CFP for the adjacent frequencies. In this case, warm starting the projection-based algorithms using solutions to CFPs at adjacent frequency points may significantly reduce the computational cost.

VI. NUMERICAL RESULTS

In this section, we apply Algorithms 1 and 2 to a family of random problems involving a chain of uncertain systems. These problems have the same band structure as described in Section V. We use the local feasibility detection method introduced in Section IV-A. Note that in order to avoid numerical problems and unnecessarily many projections, the projections are performed for slightly tighter bounds than those

0 5 10 15 20 25 10−8 10−6 10−4 10−2 k X k− X k− 1k 2 2 0 5 10 15 20 25 10−7 10−6 10−5 10−4 k S k− X kk 2 2 No. Iterations Fig. 6. Algorithm 1: X(k)− X(k−1) 2 2 and S(k)− X(k) 2 2 with

respect to the iteration number. This figure illustrates the results for 50 randomly generated problems.

0 5 10 15 20 25 30 0 5 10 15 20 25 30 35 40 45 50 No. Iterations

Fig. 7. Algorithm 2: Number of violated constraints in the CFP after each update of the global variables. This figure illustrates the results for 50 randomly generated problems.

used for feasibility checks of the local iterates. This setup is the same for all the experiments presented in this section.

Figures 5 and 6 show the behavior of Algorithm 1 for 50 randomly generated problems withm = 52. These problems

decompose into 50 subproblems which 50 agents solve col-laboratively. Figure 5 shows the number of agents that are locally infeasible at each iteration. This can be used to detect convergence to a feasible solution by checking the satisfaction of the conditions in Lemma 1. The figure demonstrates that the global variables satisfy the constraints of the original problem after at most 22 iterations. Considering Remark 1 the same performance should also be observed by checking the conditions in Lemma 2. This is confirmed by the results illustrated in Figure 6.

We performed the same experiment with Algorithm 2, and the results are illustrated in Figures 7 and 8. As can be seen from these figures, the feasibility detection based on Lemmas 1 and 2 requires at most 7 and 27 iterations to

(12)

8 0 5 10 15 20 25 30 10−20 10−10 100 k X k− X k− 1k 2 2 0 5 10 15 20 25 30 10−15 10−10 10−5 100 k S k− X kk 2 2 No. Iterations Fig. 8. Algorithm 2: X(k)− X(k−1) 2 2 and S(k)− X(k) 2 2 with

respect to the iteration number. This figure illustrates the results for 50 randomly generated problems.

0 5 10 15 20 25 30 0 50 100 150 200 250 300 350 400 No. Iterations

Fig. 9. Algorithm 2: Number of violated constraints in the CFP after each update of the global variables. This figure illustrates the results for problems with different number of constraints.

detect convergence to a feasible solution, respectively. In our experiments, Algorithm 2 is faster when feasibility detection is based on the conditions in Lemma 1, and Algorithm 1 is faster when the condition in Lemma 2 is used for feasibility detection. It is also worth mentioning that with the accelerated nonlinear Cimmino algorithm, more than 1656 iterations were needed to obtain a feasible solution.

In the next experiment, we investigate the performance of the algorithms as a function of the number of systems in the chain. The results in Figures9 and 10 indicate that the problem parameterm does not affect the performance of Algorithm 1

much. The number of iterations required to reach consensus is only increased slightly. Figures 11 and 12 verifies that Algorithm 2 behaves in a similar manner.

VII. CONCLUSION

In this paper, we have shown that it is possible to solve CPFs with loosely coupled constraints efficiently in a distributed

0 5 10 15 20 25 30 10−8 10−6 10−4 10−2 k X k− X k− 1k 2 2 0 5 10 15 20 25 30 10−8 10−6 10−4 10−2 k S k− X kk 2 2 No. Iterations Fig. 10. Algorithm 1: X(k)− X(k−1) 2 2 and S(k)− X(k) 2 2 with

respect to the iteration number. This figure illustrates the results for problems with different number of constraints.

0 5 10 15 20 25 30 0 50 100 150 200 250 300 350 400 No. Iterations

Fig. 11. Algorithm 2: Number of violated constraints in the CFP after each update of the global variables. This figure illustrates the results for problems with different number of constraints.

manner. We have proposed two algorithms. One is based on von Neumann’s AP method, and hence it enjoys the linear convergence rate that charaterizes this algorithms when applied to strictly feasible problems. The other method is based on the ADMM, and it generally outperforms the AP method in prac-tice in terms of the number of iterations required to obtain a feasible solution. Both methods can detect infeasible problems. For structured problems that arise in robust stability analysis of a large-scale weakly interconnected uncertain systems, our numerical results show that both algorithms outperform the classical projection-based algorithms.

ACKNOWLEDGMENT

This research is supported by the Swedish department of education within the ELLIIT research program.

REFERENCES

[1] D. Axehill and A. Hansson, “Towards parallel implementation of hybrid MPC – A survey and directions for future research,” in Distributed

(13)

9 0 5 10 15 20 25 30 10−20 10−10 100 k X k− X k− 1k 2 2 0 5 10 15 20 25 30 10−15 10−10 10−5 100 k S k− X kk 2 2 No. Iterations Fig. 12. Algorithm 2: X(k)− X(k−1) 2 2 and S(k)− X(k) 2 2 with

respect to the iteration number. This figure illustrates the results for problems with different number of constraints.

Decision Making and Control, ser. Lecture Notes in Control and

Information Sciences, R. Johansson and A. Rantzer, Eds. Springer Verlag, 2011.

[2] C. Langbort, R. Chandra, and R. D’Andrea, “Distributed control design for systems interconnected over an arbitrary graph,” IEEE Transactions

on Automatic Control, vol. 49, no. 9, pp. 1502 – 1519, Sep. 2004.

[3] A. Venkat, I. Hiskens, J. Rawlings, and S. Wright, “Distributed MPC strategies with application to power system automatic generation con-trol,” IEEE Transactions on Control Systems Technology, vol. 16, no. 6, pp. 1192–1206, 2008.

[4] J. Sandell, N., P. Varaiya, M. Athans, and M. Safonov, “Survey of de-centralized control methods for large scale systems,” Automatic Control,

IEEE Transactions on, vol. 23, no. 2, pp. 108 – 128, Apr. 1978.

[5] D. Dochain, G. Dumont, D. Gorinevsky, and T. Ogunnaike, “IEEE transactions on control systems technology special issue on control of industrial spatially distributed processes,” IEEE Transactions on Control

Systems Technology, vol. 11, no. 5, pp. 609 – 611, Sep. 2003.

[6] A. Bayen, R. Raffard, and C. Tomlin, “Adjoint-based constrained control of eulerian transportation networks: Application to air traffic control,” vol. 6, 2004, pp. 5539–5545.

[7] M. Jovanovic, M. Arcak, and E. Sontag, “Remarks on the stability of spatially distributed systems with a cyclic interconnection structure,” 2007, pp. 2696–2701.

[8] A. Hansen, P. Sorensen, F. Iov, and F. Blaabjerg, “Centralized power control of wind farm with doubly fed induction generators,” Renewable

Energy, vol. 31, no. 7, pp. 935–951, 2006.

[9] M. Cantoni, E. Weyer, Y. Li, S. Ooi, I. Mareels, and M. Ryan, “Control of large-scale irrigation networks,” Proceedings of the IEEE, vol. 95, no. 1, pp. 75–91, 2007.

[10] G. Stewart, D. Gorinevsky, and G. Dumont, “Feedback controller design for a spatially distributed system: The paper machine problem,” IEEE

Transactions on Control Systems Technology, vol. 11, no. 5, pp. 612–

628, 2003.

[11] K. Zhou, J. C. Doyle, and K. Glover, Robust and optimal control. Prentice Hall, 1997.

[12] K. Zhou and J. C. Doyle, Essentials of robust control. Prentice Hall, 1998.

[13] S. Skogestad and I. Postlethwaite, Multivariable feedback control. Wi-ley, 2007.

[14] A. Megretski and A. Rantzer, “System analysis via integral quadratic constraints,” IEEE Transactions on Automatic Control, vol. 42, no. 6, pp. 819 –830, Jun 1997.

[15] H. Fang and P. Antsaklis, “Distributed control with integral quadratic constraints,” vol. 17, no. 1, 2008.

[16] J. Rice and M. Verhaegen, “Robust control of distributed systems with sequentially semi-separable structure,” 2009.

[17] J. G. VanAntwerp, A. P. Featherstone, and R. D. Braatz, “Robust cross-directional control of large scale sheet and film processes,” Journal of

Process Control, vol. 11, no. 2, pp. 149 – 177, 2001.

[18] L. Grasedyck, W. Hackbusch, and B. Khoromskij, “Solution of large scale algebraic matrix Riccati equations by use of hierarchical matrices,”

Computing Vienna/New York, vol. 70, no. 2, pp. 121–165, 2003.

[19] J. Rice and M. Verhaegen, “Distributed control: A sequentially semi-separable approach for spatially heterogeneous linear systems,” IEEE

Transactions on Automatic Control, vol. 54, no. 6, pp. 1270 –1283,

Jun. 2009.

[20] Y. Censor and T. Elfving, “A nonlinear cimmino algorithm,” Department of Mathematics, University of Haifa, Tech. Rep. Report 33, 1981. [21] G. Cimmino, “Calcolo approssimato per le soluzioni dei sistemi di

equazioni lineari,” La Ricera Scientifica Roma, vol. 1, pp. 326–333. [22] A. Nedic and A. Ozdaglar, “Distributed subgradient methods for

multi-agent optimization,” IEEE Trans. Autom. Control, vol. 54, no. 1, pp. 48–61, 2009.

[23] J. N. Tsistsiklis, “Problems in decentralized decision making and com-putation,” Ph.D. dissertation, MIT, 1984.

[24] J. Tsitsiklis, D. Bertsekas, and M. Athans, “Distributed asynchronous de-terministic and stochastic gradient optimization algorithms,” Automatic

Control, IEEE Transactions on, vol. 31, no. 9, pp. 803 – 812, Sep 1986.

[25] X. Lin and S. Boyd, “Fast linear iterations for distributed averaging,” in Proceedings 42nd IEEE Conference on Decision and Control, 2003., vol. 5, Dec. 2003, pp. 4997 – 5002 Vol.5.

[26] S. Rajagopalan and D. Shah, “Distributed averaging in dynamic net-works,” IEEE Journal of Selected Topics in Signal Processing, vol. 5, no. 4, pp. 845 –854, Aug. 2011.

[27] A. N. Iusem and A. R. De Pierro, “Convergence results for an accelerated nonlinear Cimmino algorithm,” Numerische Mathematik, vol. 49, pp. 367–378, 1986.

[28] J. von Neumann, The Geometry of Orthogonal Spaces, ser. Functional Operators. Princeton University Press, 1950, vol. 2.

[29] L. Bregman, “The method of successive projection for finding a common point of convex sets,” Soviet Math. Dokl., vol. 162, pp. 487–490, 1965. [30] H. H. Bauschke and J. M. Borwein, “Dykstra’s alternating projection algorithm for two sets,” Journal of Approximation Theory, vol. 79, no. 3, pp. 418–443, 1994.

[31] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends in Machine Learning, vol. 3, no. 1, pp. 1–122, 2011.

[32] D. P. Bertsekas and J. N. Tsitsiklis, Parallel and Distributed

Computa-tion: Numerical Methods. Athena Scientific, 1997.

[33] A. Nedic, A. Ozdaglar, and P. Parrilo, “Constrained consensus and optimization in multi-agent networks,” IEEE Trans. Autom. Control, vol. 55, no. 4, pp. 922–938, April 2010.

[34] R. Glowinski and A. Marroco, “Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualité d’une classe de problèmes de dirichlet non linéaires,” Revue française d’automatique,

informatique, recherche op´erationnelle, vol. 9, no. 2, pp. 41–76, 1975.

[35] D. Gabay and B. Mercier, “A dual algorithm for the solution of nonlinear variational problems via finite element approximation,” Computers &

Mathematics with Applications, vol. 2, no. 1, pp. 17–40, 1976.

[36] R. L. Dykstra, “An algorithm for restricted least squares regression,”

Journal of the American Statistical Association, vol. 78, no. 384, pp.

837–842, 1983.

[37] H. H. Bauschke and J. M. Borwein, “On the convergence of Von Neumann’s alternating projection algorithm for two sets,” Set-Valued

Analysis, vol. 1, pp. 185–212, 1993.

[38] L. G. Gubin, B. T. Polyak, and E. V. Raik, “The method of projections for finding the common point of convex sets,” USSR Computational

Mathematics and Mathematical Physics, vol. 7, no. 6, pp. 1–24, 1967.

[39] H. H. Bauschke, J. M. Borwein, and W. Li, “Strong conical hull intersection property, bounded linear regularity, jameson’s property (g), and error bounds in convex optimization,” Mathematical Programming, vol. 86, no. 1, pp. 135–160, 1999.

[40] A. Beck and M. Teboulle, “Convergence rate analysis and error bounds for projection algorithms in convex feasibility problems,” Optimization

Methods and Software, vol. 18, no. 4, pp. 377–394, 2003.

[41] M. K. H. Fan, A. L. Tits, and J. C. Doyle, “Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics,” IEEE Trans.

Autom. Control, vol. 36, no. 1, pp. 25 –38, jan. 1991.

[42] M. S. Andersen, J. Dahl, and L. Vandenberghe, “Implementation of nonsymmetric interior-point methods for linear optimization over sparse matrix cones,” Mathematical Programming Computation, vol. 2, pp. 167–201, 2010.

[43] K. Fujisawa, S. Kim, M. Kojima, Y. Okamoto, and M. Yamashita, “User’s manual for SparseCoLo: Conversion methods for SPARSE COnnic-form Linear Optimization,” Department of Mathematical and

(14)

10

Computing Sciences, Tokyo Institute of Technology, Tokyo, Tech. Rep., 2009.

[44] S. Kim, M. Kojima, M. Mevissen, and M. Yamashita, Mathematical

Programming, pp. 1–36, 2010.

[45] M. Fukuda, M. Kojima, M. Kojima, K. Murota, and K. Nakata, “Ex-ploiting sparsity in semidefinite programming via matrix completion i: General framework,” SIAM Journal on Optimization, vol. 11, pp. 647– 674.

[46] N. Kakimura, “A direct proof for the matrix decomposition of chordal-structured positive semidefinite matrices,” Linear Algebra and its

(15)

Division,Department

DivisionofAutomaticControl

DepartmentofElectricalEngineering

Date 2011-12-02 Språk Language Svenska/Swedish Engelska/English Rapporttyp Reportcategory Licentiatavhandling Examensarbete C-uppsats D-uppsats Övrigrapport

URLförelektroniskversion

http://www.control.isy. liu .se

ISBN

ISRN

Serietitelochserienummer

Titleofseries,numbering

ISSN

1400-3902

LiTH-ISY-R-3033

Titel

Title

Decompositionand Projection Methods for Distributed Robustness Analysis of

Intercon-nectedUncertain Systems

Författare

Author

SinaKhoshfetratPakazad,MartinS.Andersen,AndersHansson,AndersRantzer

Sammanfattning

Abstract

Weconsider aclassofconvex feasibilityproblemswhere the constraintsthat describethe

feasibleset areloosely coupled. These problemsarise inrobust stabilityanalysisoflarge,

weaklyinterconnectedsystems. Tofacilitatedistributedimplementationofrobuststability

analysisofsuchsystems,wepropose twoalgorithmsbasedon decompositionand

simulta-neousprojections. TherstalgorithmisanonlinearvariantofCimmino'smeanprojection

algorithm,butbytakingthestructureoftheconstraintsintoaccount,wecanobtainafaster

rateofconvergence. Thesecondalgorithmisdevisedbyapplyingthe alternatingdirection

methodofmultiplierstoaconvexminimizationreformulationoftheconvexfeasibility

prob-lem. Weusenumericalresultstoshowthatbothalgorithmsrequirefarlessiterationsthan

theacceleratednonlinearCimminoalgorithm.

Nyckelord