Parallel and distributed optimization methods for estimation and control in networks

Introduction

In many application fields, the notion of networks has emerged as a central, unifying concept for solving different problems in systems and control theory such as analysis, process control and estimation. We live and operate in a networked world. We drive to work on networks of roads and communicate with each other using an elaborate set of devices such as phones or computers, that connect wirelessly and through the internet. Traditional networks include transportation networks (roads, rails) and networks of utilities (water, electricity, gas). But more recent examples of the increasing impact of networks include information technology networks (internet, mobile phones, acoustic networks, etc), information networks (co-author networks, bibliographic networks), social networks (collaborations, organizations), and biological and genetic networks.

These networks are often composed of multiple subsystems characterized by complex dynamics and mutual interactions such that local decisions have long-range effects throughout the entire network. Many problems associated to networked systems, such as state estimation and control, can be posed as coupled optimization problems (see e.g. , etc). Note that in these systems the interaction between subsystems gives rise to coupling in the cost or constraints, but with a specific algebraic structure, in particular sparse matrix representation that could be exploited in numerical algorithms. Therefore, in order to design an overall decision architecture for such complex networks we need to solve large coupled optimization problems but with specific structure. The major difficulty in these problems is that due to their size, communication restrictions, or requirements on robustness, often no central decisions can be taken; instead, the decisions have to be taken locally. In such a set-up, single units, or local agents, must solve local optimization subproblems and then they must negotiate their outcomes and requirements with their neighbors in order to achieve convergence to the global optimal solution. Basically, there are two general optimization approaches: (i) “Centralized” optimization algorithms: In this class the specific structure of the system is exploited, as it represents considerable sparsity in the optimization problem due to the local coupling between optimization variables (sometimes referred to as separable optimization problems). The sparsity of the problem, given by the influences between the subsystems, leads to coupling constraints represented by sparse matrices. Though parts of the algorithms will be parallelized, the parallelization in these algorithms is not restricted by e.g. limited communication between subsystems and is just for the sake of exploiting sparsity. In summary, “centralized” algorithms benefit from the sparsity induced by the networked system and solve the resulting optimization problems on a parallel computer architecture. Several standard parallel and distributed optimization methods can be found in the textbooks . Various survey papers also exist on optimization-based distributed control. In the 70’s Tamura and Mahmoud presented very comprehensive overviews. More recently, in the actual status of research in the field of coordinated optimization-based control is presented. Many different control topologies can be considered in distributed control, which have been reviewed recently in . When there is no need to solve the separable optimization problem on a parallel computer architecture, an alternative would be to solve the global optimization problem using sparse solvers that take into account the sparse structure of the problem at the linear algebra level of the optimization algorithm. In general, this choice could lead to faster algorithms in terms of CPU time than distributed or parallel algorithms. (ii) Distributed optimization algorithms (sometimes referred to as distributed multi-agent optimization algorithms): In contrast to the “centralized” algorithms, distributed algorithms on graphs have to satisfy an extra constraint, namely their computations shall be performed on all nodes in parallel, and the communication between nodes is restricted to the edges of the graph, i.e. such algorithms do not use all-to-all communication protocols. In many complex networked systems the desired behavior can be formulated as coupled optimization problems but with restrictions on communication due to the special network topology: e.g. estimation in sensor networks, consensus and rendezvous problems in multi-agent systems, resource allocation in computer networks . Some existing distributed methods that take into account explicitly information restrictions in the network combine consensus negotiations (as an efficient method for information fusion) with subgradient methods .

The goal of this paper is twofold: (i) to establish a relationship between estimation and control in networked systems and distributed optimization methods and demonstrate the effectiveness of utilizing optimization-theoretic approaches for controlling such complex systems; (ii) motivated by this connection, to build upon optimization based results to better accommodate a broader class of estimation and control problems. The core of this paper consists of Section 2, covering three applications of estimation and control that appear in the context of networked systems and then proving how we can reformulate them as coupled optimization problems. One of the key contributions of this paper is to provide an accessible, yet relatively comprehensive, overview of three classes of decomposition schemes from mathematical programming for solving distributively coupled optimization problems. We demonstrate how the decomposition schemes suggest network architectures and protocols with different properties in terms of convergence speed and coordination overhead. We also present new decomposition methods that are more efficient in terms of convergence speed than some classical decomposition schemes.

The paper is organized as follows. In Section 2 we introduce different estimation and control problems that appear in the context of complex systems with interacting subsystems dynamics and then we show how we can reformulate them as coupled optimization problems. In Section 3 we present several parallel and distributed methods for solving this type of structured optimization problems and analyze their performance. Section 3 thus serves both as a review of the necessary background and a summary of our new extensions on decomposition methods. For each of the applications, numerical experiments on different parallel and distributed algorithms are provided.

Estimation and control problems in networks

In this section we formulate different estimation and control problems for systems consisting of interconnected subsystems. In Subsection 2.1 we present a state estimation problem for a system, using a network of sensors which must exchange information in order to reach a consensus on the state estimated for the entire system. In Subsection 2.2 we will present the problem of optimal control for a large-scale system, whose subsystems are coupled with their neighbors but the objective function is decoupled. Finally, in Subsection 2.3 and 2.4 we will discuss the cooperative control problem for a group of systems (agents), which have decoupled or coupled dynamics but share a common goal.

In this section we formulate the distributed state estimation problem for systems using a sensor network based on the moving horizon estimation (MHE) approach . Sensor networks can be employed in many applications, such as monitoring, exploration, surveillance or tracking targets over specific regions. We consider the concept of MHE, as this framework offers multiple advantages: since a particular minimization problem must be solved on-line at each step, the observer is optimal with respect to the associated cost, and moreover, constraints on the state and on the noise can be taken into account .

The state estimation problem can be posed as follows. We assume that each sensor in the network measures some variables of a process, computes a local estimate of the entire state of the system, and exchanges the computed estimates with its neighbors. The solution to the estimation problem consists in finding a methodology which guarantees that all sensors asymptotically reach a reliable estimate of the overall state of the system. For the observed process we consider the following nonlinear dynamics:

For a given estimation horizon $N\geq 1$, at time $k$ given the past measurements $y_{k-N}^{i},\cdots,y_{k}^{i}$ provided by the $i$th sensor and the estimate $\hat{x}_{k-N}$, we formulate the moving horizon estimation (MHE) at $k$ as the solution to the following optimization problem :

where the matrix $\Pi_{k-N}$ is computed recursively from a Riccati difference equation in a centralized way . For the liner case, the distributed computation of this matrix can be done in many ways: e.g. using the steady-state MHE formulation (i.e. computing off-line $\Pi_{\infty}$, which is the solution of the corresponding algebraic Riccati equation) or updating $\Pi_{k-N}$ for all the sensors in the same way (using a common covariance matrix $R$ for all sensors in the Riccati difference equation update). For the nonlinear case, the update of $\Pi_{k-N}$ in a distributed fashion is still an open issue.

Note that $v_{t}^{i}=y_{t}^{i}-\theta^{i}(x_{t})$ and using the dynamics (.1), we can write $\sum_{t=k-N}^{k}||v_{t}^{i}||_{R_{i}^{-1}}^{2}$ as a function depending only on $(x_{k-N},w_{k-N},\cdots,w_{k-1})$. Therefore, by eliminating the states in (1) using the dynamics (.1) and introducing the notations:

the MHE problem (1) can be recast as an optimization problem with decoupled cost but a common decision variable x (DCx):

where the set $\textbf{X}=X\times W^{N}$.

We assume that the communication network among sensors is described by a graph $G=(V,E)$, where the nodes in $V=\{1,\cdots,M\}$ represent the sensors and the edge $(i,j)\in E\subseteq V\times V$ models that sensor $j$ sends information to sensor $i$. Then, the main challenge is to provide distributed algorithms for solving problem (1) or equivalently (DCx) which guarantee that all the sensors asymptotically reach a reliable estimate of the state variables using the information exchange model given by the graph $G$.

Example 2.1 In the particular case where the state and noise constraints $x\in X$ and $w\in W$ are described by linear inequalities (i.e. $X$ and $W$ are polyhedral sets) and the dynamics of the process and of the sensors are linear, i.e.

the MHE problem (1) can be recast as a separable convex quadratic program with decoupled cost but a common decision variable in the form (DCx):

where the matrices $H_{i}$ are positive definite and the constraint set X becomes in this case polyhedral (described only by linear inequalities).

2 Distributed optimal control problem

The application that we will discuss in this section is the distributed control of large-scale networked systems with interacting subsystem dynamics, which can be found in a broad spectrum of applications ranging from traffic networks, wind farms, to interconnected chemical plants. Distributed control is promising in applications for complex systems, since this framework allows us to design local subsystem-based controllers that take into account the interactions between different subsystems and physical constraints.

We consider discrete-time systems which can be decomposed into $M$ subsystems described by difference equations of the form:

The centralized optimal control problem over a prediction horizon $N$ reads:

where $x^{i}$ are the values of the initial state for subsystem $i$. Note that a similar formulation of distributed control for coupled subsystems with decoupled costs has been given in in the context of distributed model predictive control.

Now, we show that the optimization problem (5) can be recast as a separable optimization problem with a particular structure. To this purpose, we denote with $\textbf{X}^{i}=(X^{i})^{N}\times(U^{i})^{N}$ and

With these notations, problem (5) now reads as an optimization problem with decoupled cost and sparse coupled constraints (DCCC):

where the coupled constraints $h^{i}(\textbf{x}^{j};~j\in\mathcal{N}^{i})=0$ are obtained from the coupling between the subsystems, i.e. by stacking the constraints (.1) for a given $i$.

The centralized optimization problem (5) or (DCCC) becomes interesting if the computations can be distributed among the subsystems (agents), can be done in parallel and the amount of information that the agents must exchange is limited. In comparison with the centralized approach, a distributed strategy offers a series of advantages: first, the numerical effort is considerably smaller since we solve low dimension problems in parallel and secondly such a design is modular, i.e. adding or removing subsystems does not require any controller redesign.

Example 2.2 Many networked systems, e.g. wind farms , interconnected chemical processes , or urban traffic systems , can be decomposed into $M$ appropriate linear subsystems:

where the matrices $E_{i}$ are of appropriate dimensions and

with the matrices $A_{ij}^{-},B_{ij}^{-}$ being obtained from the matrices $A_{ij},B_{ij}$ by removing the rows with all entries equal to zero. We consider a quadratic performance index for each subsystem $i$ of the form:

where the matrices $Q_{i},R_{i}$ and $P_{i}$ are positive semidefinite. We also assume that the sets $X^{i}$ and $U^{i}$ that define the state and input constraints (4) are polyhedral. The centralized control problem over the prediction horizon $N$ for this application can be formulated as follows:

We can eliminate the state variables in the optimization problem (7) using the dynamics (.1). In this case we can define $\textbf{x}^{i}=[w_{0}^{iT}\cdots w_{N-1}^{iT}\ u_{0}^{iT}\cdots u_{N-1}^{iT}]^{T}$. Then, the control problem (7) can be recast as a separable convex quadratic program with decoupled cost and coupled constraints in the form (DCCC):

where the matrices $H^{i}$ are positive semidefinite, the local constraint sets $\textbf{X}^{i}$ are polyhedral and the coupled constraints $\sum_{i=1}^{M}G_{i}\textbf{x}^{i}=g$ are obtained from the coupling between the subsystems, i.e. by stacking the constraints (.2) for all $i,t$. Note that the number of rows of the matrices $G^{i}$ are equal to $N\sum_{i=1}^{M}p_{i}$.

3 Cooperative control problem of dynamically uncoupled systems

Cooperative control for dynamically uncoupled systems arises in a wide variety of applications like formation flying, mobile sensor networks, rendezvous problems or decentralized coordination. The cooperative control problem for dynamically uncoupled agents consists in controlling a group of independent subsystems (i.e. with decoupled dynamics), but sharing a common goal (see e.g. ).

We consider a set of $M$ identical subsystems, having the following state-space description:

The cooperative control problem over a finite horizon of length $N$, given the initial condition $x^{i}$ for each subsystem $i$, is formulated as follows:

and $\textbf{X}^{i}$ the constraint set defined by the state and input constraints (4) and by the $i$th subsystem dynamics $x_{t+1}^{i}=\phi(x^{i}_{t},u^{i}_{t})$ over the prediction horizon. Using these notations, the previous cooperative control problem can be recast as an optimization problem with coupled cost and decoupled constraints (CCDC):

We are interested in finding efficient parallel algorithms for solving problem (CCDC).

Example 2.3 We consider the formation flying for a group of satellites that are distributed along a circular orbit with independent dynamics but they have to maintain a constant distance with respect to the two nearest neighbors (see e.g. ). Using a discretized version of the linear Clohessy-Wiltshire equations of the $i$th satellite for a nominal circular trajectory :

where $x^{1,i}$, $x^{2,i}$, $x^{3,i}$ are the displacements in the radial, tangential and out-of-plane direction, $a^{1,i}$, $a^{2,i}$, $a^{3,i}$ represent the accelerations of the satellite $i$ due to propulsion or external disturbances and $\omega_{n}$ is the angular velocity at which the orbit is covered, we obtain a discrete-time linear system for the $i$th satellite of the form

Since the goal is to maintain a constant distance with respect to the two nearest neighbors, we choose the following stage cost at time $t$:

where the blocks of the positive semidefinite Hessian matrix $H=[H_{ij}]_{ij}$ satisfies $H_{ij}=0$ if $|i-j|>3$ for all $i,j$ and the sets $\textbf{X}^{i}$ are polyhedral.

4 Cooperative control problem of dynamically coupled systems

In this section we discuss the cooperation-based optimal control problem for a group of dynamically coupled subsystems . For the $i$th subsystem we consider the following linear dynamics:

Note that the dynamics described in (23) are a particular case of (6). We also assume local input constraints $u^{i}_{t}\in U^{i}$, where $U^{i}$ are convex sets.

In order to provide a cooperative behavior between subsystems we replace each local cost $f^{i}$ with one that represents the systemwide impact of local control actions. One choice is to employ a strong convex combination of local subsystems’ costs as the global objective function for the entire system. In these conditions, the cooperative control problem for coupled systems on a finite horizon $N$ will have the form:

where $\alpha_{i}>0$ and sum to 1. Note that in this form problem (26) is a particular case of problem (DCCC), where the variables associated to the $i$th subsystem are given by $[\overline{\textbf{x}}^{iT}\;\overline{\textbf{u}}^{iT}]^{T}$. However, by eliminating the states in (26) using the global dynamic model (26.1) we obtain a coupled objective function in the local variables $\textbf{x}^{i}=\overline{\textbf{u}}^{i}$ (i.e. in the local control actions) and decoupled constraints, which is a particular case of (CCDC) problem (see also Remark 2.3(ii)).

Parallel and distributed optimization algorithms for solving coupled optimization problems

In this section we present several parallel and distributed algorithms for solving the optimization problems arising in applications from estimation and control discussed in Section 2 and analyze their properties and performances, in particular we define conditions for which these algorithms converge For simplicity of the exposition, in this section we assume that all the functions are differentiable.. The presented algorithms can be classified, on the one hand in “centralized” algorithms (that in general take advantage of the sparsity of the problem and solve in parallel low dimension subproblems) and distributed algorithms (that take into account explicitly information restrictions in the network and combine consensus negotiations with optimization methods to solve distributively the problem) and on the other hand in primal and dual decomposition algorithms. The first class is based on decomposing the original optimization problem, while the second consists in decomposing the corresponding dual problem.

For a given problem representation there are often many choices of distributed algorithms, each with possible different characteristics: e.g. rate of convergence, tradeoff between local computation and global communication, and quantity of message passing. Which alternative is the best depends on the specifications of the application. However, for each algorithm we will discuss in details their main characteristics in terms of performance and properties.

In this section we study several distributed algorithms for solving separable optimization problems with decoupled cost and common decision variables in the form (DCx), that e.g. appear in the context of state estimation in sensor networks (see Section 2.1). We associate to the set of agents (e.g. sensors) a graph $G=(V,E)$ and then such distributed algorithms must satisfy the following constraint: the computations will be performed on all nodes in parallel, and the communication between nodes is restricted to the edges of the graph. Distributed optimization algorithms are mainly based on combining consensus negotiations (as an efficient method for information fusion) with optimization methods to solve distributively problems of type (DCx).

First we introduce the consensus problem for a group of $M$ agents that considers conditions under which using a certain message-passing protocol, the local variables of each agent will converge to the same value . There exist several results related to the convergence of local variables to a common value using various information exchange protocols among agents . One of the most used models for consensus is based on the following discrete-time iteration: to generate an estimate at iteration $k+1$, agent $i$ forms a convex combination of its estimate $\textbf{x}^{i}_{k}$ with the estimates received from other agents:

where $\gamma^{ij}_{k}$ represent nonnegative weights Naturally, an agent $i$ assigns zero weight to the estimates $\textbf{x}^{j}$ for those agents $j$ whose estimate information is not available at the update time. satisfying $\sum_{j}\gamma_{k}^{ij}=1$. At each iteration $k$ the information exchange among agents can be represented by a graph $(V,E_{k})$, where $E_{k}=\{(i,j):\gamma_{k}^{ij}>0\}$. We can also introduce the graph $(V,E_{\infty})$, where $E_{\infty}=\{(i,j):(i,j)\in E_{k}\;\text{for infinitely many}\;k\}$. The graphs $(V,E_{k})$ satisfy the bounded interconnection interval property if there exists an integer $\tau$ such that for any $(i,j)\in E_{\infty}$ agent $j$ sends its information to agent $i$ at least once every $\tau$ consecutive iterations. It has been proved in that under certain assumptions on the weights $\gamma^{ij}_{k}$ (e.g. stochasticity of the matrix $\Gamma_{k}=[\gamma^{ij}_{k}]_{ij}$, strong connectivity property of $(V,E_{\infty})$ and bounded interconnection interval property), the states $\textbf{x}^{i}_{k}$ of all agents converge to the same state $x^{*}$. Similar convergence results can be found in .

We return now to our optimization problem of type (DCx). In a distributed projected gradient algorithm is analyzed, which basically combines the consensus iteration presented above with a projected gradient update to generate the next estimate of the optimum. More specifically, an agent $i$ updates its estimate by combining the estimates received from its neighbors, then taking a gradient step to minimize its objective function $f^{i}$ and finally projecting on the set X:

where $\alpha_{k}$ is a common step size, $\nabla f^{i}$ denotes the gradient of the function $f^{i}$, and $[\cdot]_{\textbf{X}}$ denotes the Euclidian projection on the set X. The following convergence result holds for Algorithm dgp1 :

For the optimization problem (DCx) we assume that all the functions $f^{i}$ are convex and have bounded gradients, the set X is convex and the step size satisfies $\sum_{k}\alpha_{k}=\infty$ and $\sum_{k}\alpha_{k}^{2}<\infty$. Moreover, we assume that the weights $\gamma_{k}^{ij}$ satisfy the following properties: the matrices $\Gamma_{k}=[\gamma_{k}^{ij}]_{ij}$ are doubly stochastic, the graph $(V,E_{\infty})$ is connected and the bounded interconnection interval property holds. Then, the distributed projected gradient Algorithm dgp1 converges to an optimum of problem (DCx).

An interesting variant of a distributed gradient projected algorithm has been provided in . Compared to the previous distributed gradient Algorithm dgp1, in a fixed connected graph $(V,E)$ is taken over all iterations and the information exchange among the agents is represented by a doubly stochastic matrix $\Gamma=[\gamma^{ij}]_{ij}$ such that $\gamma_{ij}>0$ if $(i,j)\in E$. In this algorithm, first each agent implements the gradient update locally and then it runs a number $\mu$ of consensus iterations with its neighbors:

where $\Gamma^{\mu}_{ij}$ denotes the $(i,j)$ entry of the matrix $\Gamma^{\mu}$. Under similar assumptions as in Theorem 3.1, the authors in proved convergence of Algorithm dgp2 for a constant step size and for a sufficiently large $\mu$.

In the case when the set X is explicitly defined through a finite set of equalities and inequalities, an algorithm based on a penalty primal-dual approach has been recently proposed in . This algorithm allows the agents exchange information over networks with time-varying topologies and asymptotically agree on an optimal solution and the optimal value.

Another interesting approach for solving the optimization problem (DCx), but in a serial fashion, can be found in where an incremental gradient method is presented. Each step of the algorithm is a gradient iteration for a single component function $f^{i}$, and there is one step per component function. Thus, an iteration can be viewed as a cycle of $M$ subiterations, so that at $k+1$:

For convex problems, using an appropriate step size $\alpha_{k}$, the authors in show that this algorithm has much better practical rate of convergence than the classical gradient method.

Remark 3.2 (i) The convexity assumptions on the functions $f^{i}$ and the set X for convergence of the two Algorithms dgp1 and dgp2 are usually satisfied in many applications: see e.g. the state estimation problem for linear systems discussed in Example 2.1 which leads to the convex quadratic program (2). (ii) One of the main challenges when solving problems of type (DCx) is the time-dependent communication topology, as communication links can change due to changing distances, obstacles, or disturbances. While in a constant topology is assumed for Algorithm dgp2, the Algorithm dgp1 and the algorithm from are based on a changing topology, which makes them more suitable in practical applications. Moreover, the cyclical incremental algorithm can be implemented only when each agent identifies a suitable downstream and upstream neighbor. Note the existence of a cycle is a stronger assumption than connectivity. (iii) From simulations we have observed that the algorithms from are very sensitive to the choice of the weights that must be tuned, since they are considered as parameters in these methods. These algorithms do not provide a mathematical way of choosing the weights from the consensus protocol, which has a very strong influence on the convergence rate of these methods. Recently in , a distributed algorithm has been derived for solving particular cases of problems of type (DCx), where the nonnegative weights corresponding to the consensus process are interpreted as dual variables and thus they are updated using arguments from duality theory. Moreover, if the network is not densely connected (i.e. each sensor has a large number of neighbors), one can expect the performance of these algorithms from to be worse than that of the cyclic incremental gradient .

2 Decomposition algorithms for solving optimization problems (DCCC)

In this section we present several decomposition algorithms for solving separable optimization problems with decoupled cost but coupled constraints in the form (DCCC). Distributed control for complex processes with interacting subsystem dynamics usually leads to such optimization problems (see e.g. Section 2.2). We discuss two classes of decomposition principles: primal and dual. We use the terms primal and dual in their mathematical programming meaning: primal indicates that the optimization problems are solved using the original formulation and variables and dual indicates that the original problem has been rewritten using Lagrangian relaxation.

Compared to the general formulation of problem (DCCC), we focus in this section on decomposition methods for the particular case of separable convex problems with decoupled cost and coupled constraints For the nonconvex case of problem (DCCC) we can still obtain decomposition algorithms by combining sequential quadratic programming or sequential convex programming, in order to linearize the nonlinear coupled constraints, with decomposition methods that address the decomposable convex problems (see e.g. ).:

Each function $f^{i}$ is convex quadratic and $\textbf{X}^{i}$ are compact convex sets. Moreover, the Slater’s condition holds, i.e. there exist $\textbf{x}^{i}\in\text{int}(\textbf{X}^{i})$ such that $\sum_{i=1}^{M}G_{i}\textbf{x}^{i}=g$.

From Example 2.2 we have seen that centralized optimal control for interconnected linear systems leads to such a separable convex quadratic formulation, e.g. (8).

We begin with primal decomposition (see e.g. and the references therein). We can decompose the original problem (conv-DCCC) as follows: we introduce some auxiliary variables in order to separate the coupled linear equality constraints, i.e. we introduce the new variables $\textbf{t}^{1},\cdots,\textbf{t}^{M-1}$, and obtain $M$ subproblems:

for $i=1,\cdots,M-1$ and the $M$th subproblem

We now discuss dual decomposition . In dual decomposition methods we have the following economic interpretation: the master problem sets the prices for the resources to each subproblem which has to decide the amount of resources to be used depending on the price. The iteration continues until the best pricing strategy is obtained. Clearly, if the coupled constraints $\sum_{i}G_{i}\textbf{x}^{i}=g$ are absent, then the problem (conv-DCCC) can be decoupled. Therefore it makes sense to relax these coupled constraints using duality theory. We construct the partial augmented Lagrangian:

where $\mu>0$ and the functions $P_{\textbf{X}^{i}}$ associated to the sets $\textbf{X}^{i}$ (usually called prox functions) must have certain properties explained below. We also define the corresponding augmented dual function:

and from the structure of $L_{\mu}$ we obtain that (28) decouples in $M$ subproblems

We are interested in the properties of the family of augmented dual functions $\{d_{\mu}\}_{\mu>0}$. Note that $\lim_{\mu\to 0}d_{\mu}(\lambda)=d_{0}(\lambda)$, where $d_{0}(\lambda)=\min_{\textbf{x}^{i}\in\textbf{X}^{i}}L_{0}(\textbf{x},\lambda)$ is the standard dual function, whenever the prox functions $P_{\textbf{X}^{i}}$ are chosen to be continuous on the compact sets $\textbf{X}^{i}$ or are barrier functions associated to these sets (see ). The goal is to maximize the augmented dual function for $\mu$ sufficiently small:

in order to find an approximation of the optimal Lagrange multiplier $\lambda^{*}=\arg\max_{\lambda}d_{0}(\lambda)$ and then to recover an approximation of the corresponding optimal primal variables $\textbf{x}^{i*}$. We distinguish three algorithms, depending on the choice of the constant $\mu$ and of the prox functions $P_{\textbf{X}^{i}}$:

dual subgradient algorithm: $\mu=0$ and $P_{\textbf{X}^{i}}=0$

dual fast gradient algorithm: $\mu>0$ and $P_{\textbf{X}^{i}}$ are strongly convex functions

dual interior-point algorithm: $\mu>0$ and $P_{\textbf{X}^{i}}$ are barrier functions for the sets $\textbf{X}^{i}$.

The next theorem provides the main properties of the augmented dual function:

Under Assumption 3.3, the augmented dual function $d_{\mu}$ is characterized as follows: (I) For any $\mu\geq 0$ and convex functions $P_{\textbf{X}^{i}}$ a subgradient of $d_{\mu}$ at $\lambda$ is given by $\sum_{i}G_{i}\textbf{x}^{i}(\mu,\lambda)-g$. (II) For $\mu>0$ and strong convex functions $P_{\textbf{X}^{i}}$ the function $d_{\mu}$ has a Lipschitz continuous gradient. (III) For $\mu>0$ and barrier functions $P_{\textbf{X}^{i}}$ the function $d_{\mu}$ is self-concordant.

We denote $\textbf{x}_{k}^{i}=\textbf{x}^{i}(\mu_{k},\lambda_{k})$. The iterations of the three algorithms are:

where $\alpha_{k}$ is a step-size that can be chosen as in Remark 3.2 for algorithm (DS) or satisfying Armijo rule for algorithm (DIP), $L_{\mu}$ is the Lipschitz constant of the gradient $\nabla d_{\mu}$ and $\beta_{k}>0$ is defined iteratively as in . Moreover, in the dual interior-point algorithm (DIP) we have an outer iteration in $p$ where we decrease $\mu_{p}\to 0$ and an inner iteration in $k$ where we need to generate vectors close to the central path using Newton updates with $\nabla^{2}d_{\mu}(\lambda)$ representing the Hessian of the augmented dual function $d_{\mu}$ at $\lambda$ (see for more details).

The convergence of these three algorithms (DS), (DFG) and (DIP) can be established under suitable assumptions on problem (conv-DCCC) and on the prox functions $P_{\textbf{X}^{i}}$:

If Assumption 3.3 holds for the separable convex problem (conv-DCCC), then all three algorithms (DS), (DFG) and (DIP) are convergent under a suitable choice of the step-size. Moreover, the dual fast gradient algorithm (DFG) has complexity ${\mathcal{O}}(\frac{c_{1}}{\epsilon})$, while the dual interior-point algorithm (DIP) has complexity ${\mathcal{O}}\left(c_{2}\log(\frac{c_{3}}{\epsilon})\right)$, where $\epsilon$ is the accuracy of the approximation of the optimum for problem (conv-DCCC) and $c_{i}$ are some positive constants.

We should note that in the primal subgradient algorithm we maintain feasibility of the coupled constraints in the problem (conv-DCCC) at each iteration while for the dual algorithms feasibility holds only at convergence of these algorithms and not at the intermediate iterations. Since for control problems the coupled constraints represent the dynamics of the networked system over the prediction horizon, when using a dual algorithm these dynamics will be satisfied only at convergence. This is a major issue when we stop at an intermediate step of a dual based algorithm.

There are also other dual decomposition methods based on the concept of augmented Lagrangians: e.g. the alternating direction method , where a quadratic penalty term $\mu||\sum_{i}G_{i}\textbf{x}^{i}-g||^{2}$ is added to the standard Lagrangian $L_{0}$. A computational drawback of this scheme is that the quadratic penalty term is not separable in $\textbf{x}^{i}$. However, this is overcome by carrying out the minimization problem in a Gauss-Seidel fashion, followed by a steepest ascent update of the multipliers. In other dual decomposition methods, such as partial inverse method or proximal point method , for example a term of the form $\mu\sum_{i}||\textbf{x}^{i}-\textbf{x}^{i}_{k}||^{2}$ is added to the Lagrangian $L_{0}$. These schemes have been shown to be very sensitive to the value of the parameter $\mu$, with difficulties in practice to obtain the best convergence rate. Some heuristics for choosing $\mu$ can be found in the literature . However, these heuristics have not been formally analyzed from the viewpoint of efficiency estimates for the general case (linear convergence results have been obtained e.g. only for strongly convex functions).

The new decomposition methods called here “dual fast gradient” (DFG) and “dual interior-point” (DIP) obtained by smoothing the Lagrangian are more efficient in terms of number of iterations compared to the classical primal or dual subgradient algorithm (see also Table 2). We should note however, that algorithm (DFG) is more appropriate than the algorithm (DIP) when solving problems where the number of coupling constraints is large, since for (DIP) we need to invert at each iteration a square matrix of dimension $n_{\lambda}$, where $n_{\lambda}$ denotes the dimension of $\lambda$ (or equivalently the number of rows in the matrices $G_{i}$).

It is also clear that the update rules in algorithms (DS) and (DFG) are completely distributed, according to the communication graph between subsystems. Indeed, we recall that the coupling constraints $h^{i}(\textbf{x}^{j};~j\in\mathcal{N}^{i})=0$ in problem (conv-DCCC) are assumed to be linear, of type $G^{i}[\textbf{x}^{j}]_{j\in\mathcal{N}^{i}}=g_{i}$, i.e. we have $[G_{1}\cdots G_{M}]=[G^{1T}\cdots G^{MT}]^{T}$. Let $\lambda^{i}$ be the Lagrange multipliers for the constraints $G^{i}[\textbf{x}^{j}]_{j\in\mathcal{N}^{i}}=g_{i}$, and thus $\lambda=[\lambda^{1T}\cdots\lambda^{MT}]^{T}$. Then, the main update rules in Algorithms (DS) and (DFG) are distributed, each agent $i$ using information only from its neighbors, e.g.:

However, for the algorithm (DIP), the update of the Lagrange multiplier has to be done by a central agent, i.e. in this case we have a star-shaped topology for the communication among subsystems. Note that for this algorithm the sparsity of the graph will impose sparsity on the matrices $G_{i}$, which in turn will have a strong effect on the computation of the Hessian of the corresponding dual function (see for more details).

3 Parallel algorithms for solving optimization problems of type (CCDC)

In this section we study parallel algorithms for solving optimization problems with coupled cost but decoupled constraints in the form (CCDC), that e.g. appear in the context of cooperative control (see Sections 2.3 and 2.4). A well known parallel algorithm in linear algebra for solving systems of linear equations is the Jacobi algorithm that can be also used in the context of optimization . Applying Jacobi algorithm, we decompose our optimization problem of type (CCDC) into $M$ optimization subproblems of lower dimension. In this algorithm each agent updates its variable $\textbf{x}^{i}$ by solving a low dimension optimization problem where the values of the rest of variables are calculated at the previous iteration. An extension of the Jacobi algorithm is the Gauss-Seidel algorithm, where at each iteration each agent updates its variable by solving an optimization problem for which the rest of the variables are replaced with the most recent values computed.

It is clear that in the Jacobi algorithm the optimization subproblems can be solved in parallel at each iteration. The Gauss-Seidel algorithm can be also parallelized, providing that a coloring scheme can be applied (see for more details).

The convergence of these two algorithms can be established under suitable contraction assumptions on the mapping $\textbf{x}-\beta\Delta f(\textbf{x})$ with respect to the block-maximum norm $\|\textbf{x}\|=\max_{i}\|\textbf{x}^{i}\|/\zeta_{i}$ , where the $\zeta_{i}$’s are positive scalars and $\textbf{x}=[\textbf{x}^{1T}\cdots\textbf{x}^{MT}]^{T}$.

. For the optimization problem (CCDC) we assume that the objective function $f$ is differentiable and suppose that the mapping $\textbf{x}-\beta\Delta f(\textbf{x})$ is a contraction for some positive scalar $\beta$. Then, the Jacobi and Gauss-Seidel algorithms are well defined and the sequence $\{\textbf{x}_{k}\}_{k}$ converges to the minimum of (CCDC) linearly for both iterations.

For the Gauss-Seidel algorithm, the assumptions for convergence given in Theorem 3.7 can be relaxed, in particular the contraction assumption can be replaced with a convexity assumption on the objective function ($f$ needs to be differentiable and convex and, furthermore, the function $f$ needs to be strictly convex function of $\textbf{x}^{i}$ when the values of all the other components of x are held constant, for each $i$), see for more details. If $f$ is not differentiable, the Jacobi or Gauss-Seidel algorithm can fail to converge to the minimum of (CCDC) because it can stop at a non-optimal “corner” point at which $f$ is non-differentiable and from which $f$ cannot be reduced along any coordinate. The contraction assumption on the functions $f$ for convergence of these two algorithms is usually satisfied in many applications: see e.g. the cooperative control problem for satellite formation discussed in Example 2.3 which leads to the convex quadratic program (2.3) for which the Hessian satisfies the contraction assumption or the application from Section 2.4.

In the optimization problem (CCDC) has been solved using a coordinate descent method. The iteration $k+1$ of the algorithm has the following form:

where $i_{k}$ is chosen randomly based on a uniform distribution. Moreover, we assume componentwise Lipschitz continuity of the gradient of $f$ with the Lipschitz constant $L_{i}$, for all $i=1,\cdots,M$. In Nesterov proves ${\mathcal{O}}(\frac{1}{\epsilon})$ rate of convergence in probability for the coordinate descent algorithm.

For cooperative control problems of dynamically coupled systems (see Section 2.4), which also leads to optimization problems of the form (CCDC), various versions of Jacobi-based algorithms have been proposed in the literature. For example in the authors have proposed an algorithm of the following form:

where $\alpha_{i}$ are positive weights, summing to $1$. In the authors have shown that all the limit points of the sequence generated by the previous algorithm are optimal.

In the authors have proposed a decomposition of the problem (CCDC) into a set of local subproblems that are solved iteratively by a network of agents. Each subproblem ˆ is obtained from ˆ(CCDC) discarding from the objective $f$ the terms that do not depend on $\textbf{x}^{i}$ and with the constraint set $\textbf{X}^{i}$. A distributed algorithm based on the method of feasible directions has been proposed to generate the iterations of the agents:

where the local descent direction is $d_{k}^{i}=\hat{\textbf{x}}^{i}_{k}-\textbf{x}^{i}_{k}$, for ˆ$\hat{\textbf{x}}^{i}_{k}\in\textbf{X}^{i}$, and the step size $\alpha_{k}^{i}$ satisfies the Armijo rule . The local iterations require relatively low effort and arrive at a solution of (CCDC) at the expense of slower convergence and high communication among neighboring agents.

From the Tables 1, 2 and 3 we can observe that, in order to get an optimal solution, we need to perform a large number of iterations. Note however that in practical applications from control it is not always necessary to get an optimal solution, but we can also use a suboptimal solution that can still preserve some fundamental properties for the system such as robustness, stability, etc. Whenever a suboptimal solution is satisfactory we can stop the optimization algorithm at an intermediate iteration. Note that there exist many control strategies based on this principle of suboptimality (see e.g. ).

Conclusions

This paper has presented three applications from estimation and process control for networked systems that lead to coupled optimization problems with particular structure that can be exploited in decomposition algorithms. A systematic framework is then developed in the paper to explore several parallel and distributed algorithms for solving such structured optimization problems, each with a different tradeoff among convergence speed, message passing amount, and distributed computation architecture. For each application, numerical experiments on several parallel and distributed algorithms are provided.

References