Efficient parallel coordinate descent algorithm for convex optimization problems with separable constraints: application to distributed MPC

Introduction

Model predictive control (MPC) has become a popular advanced control technology implemented in network systems due to its ability to handle hard input and state constraints . Network systems are usually modeled by a graph whose nodes represent subsystems and whose arcs indicate dynamic couplings. These types of systems are complex and large in dimension, whose structures may be hierarchical and they have multiple decision-makers (e.g. process control , traffic and power systems , flight formation ).

Decomposition methods represent a very powerful tool for solving distributed MPC problems in network systems. The basic idea of these methods is to decompose the original large optimization problem into smaller subproblems. Decomposition methods can be divided in two main classes: primal and dual decomposition methods. In primal decomposition the optimization problem is solved using the original formulation and variables via methods such as interior-point, feasible directions, Gauss-Jacobi type and others . In dual decomposition the original problem is rewritten using Lagrangian relaxation for the coupling constraints and the dual problem is solved with a Newton or (sub)gradient algorithm . In cooperative based distributed MPC algorithms are proposed based on Gauss-Jacobi iterations, where asymptotic convergence to the centralized solution and feasibility for their iterates is proved. In non-cooperative algorithms are derived for distributed MPC problems, where communication takes place only between neighbors. In a distributed algorithm based on interior-point methods is proposed whose iterates converge to the centralized solution. In dual distributed gradient algorithms based on Lagrange relaxation of the coupling constraints are presented for solving MPC problems, algorithms which usually produce feasible and optimal primal solutions in the limit. While much research has focused on a dual approach, our work develops a primal method that ensures constraint feasibility, has low iteration complexity and provides estimates on suboptimality.

Further, MPC schemes tend to be quite costly computation-wise compared with classical control methods, e.g. PID controllers, so that for these advanced schemes we need hardware with a reasonable amount of computational power that is embedded on the subsystems. Therefore, research for distributed and embedded MPC has gained momentum in the past few years. The concept behind embedded MPC is designing a control scheme that can be implemented on autonomous electronic hardware, e.g programmable logic controllers (PLC) or field-programmable gate arrays (FPGAs) . Such devices vary widely in both computational power and memory storage capabilities as well as cost. As a result, there has been a growing focus on making MPC schemes faster by reducing problem size and improving the computational efficiency through decentralization , moving block strategies (e.g. by using latent variables or Laguerre functions ) and other procedures, allowing these schemes to be implemented on cheaper hardware with little computational power.

The main contribution of this paper is the development of a parallel coordinate descent algorithm for smooth convex optimization problems with separable constraints that is computationally efficient and thus suitable for MPC schemes that need to be implemented distributively or in hardware with limited computational power. This algorithm employs parallel block-coordinate updates for the optimization variables and has similarities to the optimization algorithm proposed in , but with simpler implementation, lower iteration complexity and guaranteed rate of convergence. We derive (sub)linear rate of convergence for the new algorithm whose proof relies on the Lipschitz property of the gradient of the objective function. The new parallel algorithm is used for solving MPC problems for general linear network systems in a distributed fashion using local information. For ensuring stability of the MPC scheme, we use a terminal cost formulation derived from a distributed synthesis and we eliminate the need for a terminal state constraint. Compared with the existing approaches based on an end point constraint, we reduce the conservatism by combining the underlying structure of the system with distributed optimization . Because the MPC optimization problem is usually terminated before convergence, our MPC controller is a form of suboptimal control. However, using the theory of suboptimal control we can still guarantee feasibility and stability.

This paper is organized as follows. In Section 2 we derive our parallel coordinate descent optimization algorithm and prove the convergence rate for it. In Sections 3.1-3.2 we introduce the model for general network systems, present the MPC problem with a terminal cost formulation and provide the means for which this terminal cost can be synthesized distributively. In Sections 3.3-3.4 we employ our algorithm for distributively solving MPC problems arising from network systems and discuss details regarding its implementation. In Section 4 we compare its performance with other algorithms and test it on a real application - a quadruple water tank process.

A parallel coordinate descent algorithm for smooth convex problems with separable constraints

In this section we propose a parallel coordinate descent based algorithm for efficiently solving the general convex optimization problem of the following form:

Let us partition the identity matrix in accordance with the structure of the decision variable $\mathbf{u}$:

Due to the assumption that $f$ is coordinate-wise Lipschitz continuous, it can be easily deduced that :

which will prove useful for estimating the rate of convergence for our algorithm. Additionally, if function $f$ is smooth and strongly convex with regards to $\left\|\cdot\right\|_{1}$ with a parameter $\sigma_{1}$, then :

Note that if $f$ is strongly convex w.r.t the standard Euclidean norm $\left\|\cdot\right\|$ with a parameter $\sigma_{0}$, then $\sigma_{0}\geq\sigma_{1}L^{i}_{\max}$, where $L^{i}_{\max}=\displaystyle\max_{i}L_{i}$. By taking $\mathbf{w}=\mathbf{v}+E^{i}h_{i}$ and $\mathbf{v}=\mathbf{u}$ in (5) we also get:

and combining with (3) we also deduce that $\sigma_{1}\leq 1$.

We now define the constrained coordinate update for our algorithm:

The optimality conditions for the previous optimization problem are:

Taking $\mathbf{v}^{i}=\mathbf{u}^{i}$ in the previous inequality and combining with (3) we obtain the following decrease in the objective function:

We now present our Parallel Coordinate Descent Method, that resembles the method in but with simpler implementation, lower iteration complexity and guaranteed rate of convergence, and is a parallel version of the coordinate descent method from :

Algorithm PCDM Choose $\mathbf{u}_{0}^{i}\in\mathbf{U}^{i}$ for all $i=1,\dots,M$. For $k\geq 0$: 1. Compute in parallel $\mathbf{\bar{v}}^{i}(\mathbf{u}_{k}),\;i=1,\dots,M.$ 2. Update in parallel: $\mathbf{u}_{k+1}^{i}=\frac{1}{M}\mathbf{\bar{v}}^{i}(\mathbf{u}_{k})+\frac{M-1}{M}\mathbf{u}_{k}^{i},\;i=1,\dots,M.$

From (6)-(7), convexity of $f$ and $\mathbf{u}_{k+1}=\sum_{i}\frac{1}{M}\bar{\mathbf{u}}^{i}(\mathbf{u}_{k})$ we see immediately that method PCDM decreases strictly the objective function at each iteration, provided that $\mathbf{u}_{k}\neq\mathbf{u}_{*}$, where $\mathbf{u}_{*}$ is the optimal solution of (1), i.e.:

Let $f^{*}$ be the optimal value in optimization problem (1). The following theorem derives convergence rate of Algorithm PCDM and employs standard techniques for proving rate of convergence of the gradient method :

If function $f$ in optimization problem (1) has a coordinate-wise Lipschitz continuous gradient with constants $L_{i}$ as given in (2), then Algorithm PCDM has the following sublinear rate of convergence:

where $r_{0}=\left\|\mathbf{u}_{0}\!-\!\mathbf{u}_{*}\right\|_{1}$.

where $\mathbf{u}_{*}$ is the optimal solution of (1) and $\mathbf{u}^{i}_{*}=(E^{i})^{T}\mathbf{u}_{*}$. Then, using similar derivations as in , we have:

Taking into account that our algorithm is a descent algorithm, i.e. $f(\mathbf{u}_{j})\geq f(\mathbf{u}_{k+1})$ for all $j\leq k$ and by the previous inequality the proof is complete. ∎

Now, we derive linear convergence rate for Algorithm PCDM, provided that $f$ is additionally strongly convex:

Under the assumptions of Theorem 1 and if we further assume that $f$ is strongly convex with regards to $\left\|\cdot\right\|_{1}$ with a constant $\sigma_{1}$ as given in (5), then the following linear rate of convergence is achieved for Algorithm PCDM:

We take $\mathbf{w}=\mathbf{u}^{*}$ and $\mathbf{v}=\mathbf{u}_{k}$ in (5) and through (9) we get:

From the strong convexity of $f$ in (5) we also get:

We now define $\gamma=\frac{2\sigma_{1}}{1+\sigma_{1}}\in$ and using the previous inequality we obtain the following result:

Applying this inequality iteratively, we obtain the following for $k\geq 0$:

and by replacing $\gamma=\frac{2\sigma_{1}}{1+\sigma_{1}}$ we complete the proof. ∎

The following properties follow immediately for our Algorithm PCDM.

For the optimization problem (1), with the assumptions of Theorem 2, we have the following statements: (i) Given any feasible initial guess $\mathbf{u}_{0}$, the iterates of the Algorithm PCDM are feasible at each iteration, i.e. $\mathbf{u}_{k}^{i}\in\mathbf{U}^{i}$ for all $k\geq 0$. (ii) The function $f$ is nonincreasing, i.e. $f(\mathbf{u}_{k+1})\leq f(\mathbf{u}_{k})$ according to (8). (iii) The sub(linear) rate of convergence of Algorithm PCDM is given in Theorem 1 (Theorem 2).

Application of Algorithm PCDM to distributed suboptimal MPC

The Algorithm PCDM can be used to solve distributively input constrained MPC problems for network systems after state elimination. In this section we show that the MPC scheme obtained by solving approximately the corresponding optimization problem with Algorithm PCDM is stable and distributed.

In this paper we consider discrete-time network systems, which are usually modeled by a graph whose nodes represent subsystems and whose arcs indicate dynamic couplings, defined by the following linear state equations :

We can now formulate the MPC problem for system (11) over a prediction horizon of length $N$ and a given initial state $x$ as :

Further, we denote the approximate solution produced by Algorithm PCDM for problem (14) after certain number of iterations with $\mathbf{u}^{\text{CD}}$. We also consider that at each MPC step the Algorithm PCDM is initialized (warm start) with the shifted sequence of controllers obtained at the previous step and the feedback controller $\kappa(\cdot)$ computed in Section 3.2 below. The suboptimal MPC scheme corresponding to (14) would now be:

2 Distributed synthesis for a terminal cost

The task of finding suitable $P^{i}$, $F^{i}$ and $W^{\mathcal{N}^{i}}$ matrices is now reduced to the following optimization problem:

It can be easily observed that if the optimal value $\delta^{*}\leq 0$, consequently $W\leq 0$ and (16) holds. This optimization problem, in its current nonconvex form, cannot be solved efficiently. However, it can be recast as a sparse SDP if we can reformulate (18) as an LMI. We need now to make the assumption that all the subsystems have the same dimension for the states, i.e. $n_{i}=n_{j}$ for all $i,j$. Subsequently, we introduce the well-known linearizations: $P^{i}=(S^{i})^{-1},\;F^{i}=Y^{i}G^{-1}$ and a series of matrices that will be of aid in formulating the LMIs:

where the $0$ blocks are of appropriate dimensions By $I_{n}$ we denote the identity matrix of size $n\times n$, by $\otimes$ we denote the standard Kronecker product and by $\left|{\mathcal{N}^{i}}\right|$ the cardinality of the set $\mathcal{N}^{i}$..

has an optimal value $\delta^{*}\leq 0$, then (16) holds By $*$ we denote the transpose of the symmetric block of the matrix..

From (21) we observe that $S^{\mathcal{N}^{i}}\succ 0$, so that $(S^{\mathcal{N}^{i}}-G^{\mathcal{N}^{i}})^{T}(S^{\mathcal{N}^{i}})^{-1}(S^{\mathcal{N}^{i}}\\ -G^{\mathcal{N}^{i}})\succeq 0$, which in turn implies

If we apply the Schur complement to (21), we obtain:

There exist in literature many optimization algorithms (see e.g. ) for solving distributively sparse SDP problems in the form (20).

3 Stability of the MPC scheme

satisfies (16). Then, using Theorem 3 from we have that our MPC controller stabilizes asymptotically the system for all initial states $x\in X_{N}$, where

4 Distributed implementation of the MPC scheme based on Algorithm PCDM

In this section we discuss some technical aspects for the distributed implementation of the MPC scheme derived above when using Algorithm PCDM to solve the control problem (14). Usually, in the linear MPC framework, the local stage and final cost are taken of the following quadratic form:

where $\mathbf{Q}$ is positive definite due to the assumption that all $R^{i}$ are positive definite. Usually, for the dynamics (11) the corresponding matrices $\mathbf{Q}$ and $\mathbf{W}$ obtained after eliminating the states are dense and despite the fact that Algorithm PCDM can perform parallel computations (i.e. each subsystem needs to solve small local problems) we need all to all communication between subsystems. However, for the dynamics (12) the corresponding matrices $\mathbf{Q}$ and $\mathbf{W}$ are sparse and in this case in our Algorithm PCDM we can perform distributed computations (i.e. the subsystems solve small local problems in parallel and they need to communicate only with their neighborhood subsystems as detailed below). Indeed, if the dynamics of the system are given by (12), then

and thus the matrices $\mathbf{Q}$ and $\mathbf{W}$ have a sparse structure (see also ). Let us define the neighborhood subsystems of a certain subsystem $i$ as $\hat{\mathcal{N}}^{i}={\mathcal{N}}^{i}\cup\{l:\;l\in{\mathcal{N}}^{j},j\in\bar{\mathcal{N}}^{i}\}$, where $\bar{\mathcal{N}}^{i}=\{j:\;i\in{\mathcal{N}}^{j}\}$, then the matrix $\mathbf{Q}$ has all the $(i,j)$ block matrices $\mathbf{Q}^{ij}=0$ for all $j\notin\hat{\mathcal{N}}^{i}$ and the matrix $\mathbf{W}$ has all the block matrices $\mathbf{W}^{ij}=0$ for all $j\notin\bar{\mathcal{N}}^{i}$, for any given subsystem $i$. As a result, we can express the objective function of problem (1) as a sum of local functions with sparse structure:

Thus, the $i$th block components of $\nabla f$ can be computed using only local information:

Note that in Algorithm PCDM the only parameters that we need to compute are the Lipschitz constants $L_{i}$. However, in the MPC problem, $L_{i}$ does not depend on the initial state $x$ and can be computed locally by each subsystem as: $L_{i}=\lambda_{\max}(\mathbf{Q}^{ii})$. From the previous discussion it follows immediately that the iterations of Algorithm PCDM can be performed in parallel using distributed computations (see (24)).

Further, our Algorithm PCDM has a simpler implementation of the iterates than the algorithm from : in Algorithm PCDM the main step consists of computing local projections on the sets $\mathbf{U}^{i}$ (in the context of MPC usually these sets are simple and the projections can be computed in closed form); while in the algorithm from this step is replaced with solving local dense QP problems with the feasible set given by $\mathbf{U}^{i}$ (even in the context of MPC this local QP problems cannot be solved in closed form and an additional QP solver needs to be used). Finally, the number of iterations for finding an approximate solution can be easily predicted in our Algorithm (see Theorems 1 and 2), while in the algorithm from the authors prove only asymptotic converge.

Numerical Results

Since our Algorithm PCDM has similarities with the algorithm from , in this section we compare these two algorithms on controlling a laboratory setup with DMPC (4 tank process) and on MPC problems for random network systems of varying dimension.

To demonstrate the applicability of our Algorithm PCDM, we apply this newly developed method for solving the optimization problems arising from the MPC problem for a process consisting of four interconnected water tanks, see Fig. 1 for the process diagram, whose objective is to control the level of water in each of the four tanks. For this plant, there are two types of system inputs that can be considered: the pump flows, when the ratios of the three way valves are considered fixed, or the ratios of the three way valves, whilst having fixed flows from the pumps. In this paper, we consider the latter option, with the valve ratios denoted by $\gamma_{a}$ and $\gamma_{b}$, such that tanks 1 and 3 have inflows $\gamma_{a}q_{a}$ and $(1-\gamma_{a})q_{a}$, while tanks 2 and 4 have inflows $\gamma_{b}q_{b}$ and $(1-\gamma_{b})q_{b}$. The simplified continuous nonlinear model of the plant is well known . We use the following notation: $h_{i}$ are the levels and $a_{i}$ are the discharge constants of tank $i$, $S$ is the cross section of the tanks, $\gamma_{a}$, $\gamma_{b}$ are the three-way valve ratios, both in $$, while$q_{a}$and$q_{b}$ are the pump flows.

The discharge constants $a_{i}$, with $i=1,\dots,4$ and the other parameters of the model are determined experimentally from our laboratory setup (see Table 1). We can obtain a linear continuous state-space model by linearizing the nonlinear model at an operating point given by $h_{i}^{0}$, $\gamma_{a}^{0}$, $\gamma_{b}^{0}$, and the maximum inflows from the pumps, with the deviation variables $x^{i}=h_{i}-h_{i}^{0}$, $u^{1}=\gamma_{a}-\gamma_{a}^{0}$, $u^{2}=\gamma_{b}-\gamma_{b}^{0}$:

where $\tau_{i}=\frac{S}{a_{i}}\sqrt{\frac{2h_{i}^{0}}{g}}$, $i=1,\dots,4$, is the time constant for tank $i$.

Using zero-order hold method with a sampling time of $5$ seconds we obtain the discrete time model of type (12), with the partition $x^{1}\leftarrow\left[x^{1}~x^{4}\right]^{T}$ and $x^{2}\leftarrow\left[x^{2}~x^{3}\right]^{T}$. For the input constraints of the MPC scheme we consider the practical constraints of the ratios of the three way valves for our plant, i.e $u^{i}\in[0.15,\ 0.8]-\gamma^{i}_{0}$, where $\gamma^{i}_{0}$ is the linearization input. Due to the fact that our plant has overflow sensors fitted to the tanks and an emergency shutoff program, we do not introduce constraints for the states. For the stage cost we have taken the weighting matrices to be $Q^{i}=I_{n_{i}}$ and $R^{i}=0.01I_{m_{i}}$.

2 Implementation of the MPC scheme using MPI

In this section we underline the benefits of Algorithm PCDM when it is implemented in an appropriate fashion for the quadruple tank MPC scheme. We implemented for comparison, Algorithm PCDM and that of . Both algorithms were implemented in C programming language, with parallelization ensured via MPI and linear algebra operations done with CLAPACK. Algorithm requires solving, at each step, $2$ QP problems in parallel, problems which cannot be solved in closed form. For solving these QP problems, we use the qpip routine of the QPC toolbox . The algorithms were implemented on a PC, with 2 Intel Xeon E5310 CPUs at 1.60 GHz and 4Gb of RAM. For the MPC problem in this subsection we control the plant such that the levels and inputs will reach those of the steady state linearization values $h^{0}$ and $\gamma^{0}$.

Figure 2 outlines a comparison of the two algorithms for solving this quadruple tank MPC problem, considering a prediction time of $150$ seconds, for different prediction horizons $N$ and sampling time $\tau$, such that $\tau N=150$ seconds. The bar values represent the total sum $\displaystyle\sum_{t=1}^{50}V_{N}(x_{t},\mathbf{u})$ for $50$ MPC steps, where $\mathbf{u}$ is calculated either with PCDM or with the algorithm from . For the same $50$ simulation steps, we outline in Table 2 a comparison of the average number of iterations achieved by both algorithms and the performance loss, i.e. a percentile difference between the suboptimal cost achieved in Figure 2 ($\displaystyle\sum_{t=1}^{50}V_{N}(x_{t},\mathbf{u})$) and the optimal costs that were precalculated with Matlab’s quadprog ($\displaystyle\sum_{t=1}^{50}V_{N}^{*}(x_{t})$), both for the time $\tau$ and prediction horizon $N$. Note that our total cost is usually better than that of when the available time is short ($\tau<2$) and for $\tau\geq 2$ both algorithms solve the corresponding optimization problem exactly. Also note that, due to its low complexity iteration, our algorithm performs more than ten times the amount of iterations than the algorithm from .

3 Implementation of the MPC scheme using Siemens S7-1200 PLC

Due to this cycle time, the limited size of the S7-1200’s work memory and its processing speed, the number of iterations of the Algorithm PCDM that can be computed are also limited. In Table 3 the number of iterations available per prediction horizon, included in the $5$ seconds cycle time, and the memory requirements for these prediction horizons are presented. Although the numbers of computed iterations seem small, we have found in practice that the suboptimal MPC scheme still stabilizes the quadruple tank process and ensures set point tracking.

The results of the control process are presented in Fig. 3 for a prediction horizon $N=20$: the continuous lines represent the evolution of water levels in each of the four tanks, while the dashed lines are their respective set points. We choose two set points. We first let the plant get near its first setpoint, after which we choose a new set point which is an equilibrium point for the plant. As it can be observed from the figure, the MPC scheme still steers the process to the respective set points.

4 Implementation of MPC scheme for random network systems

We now wish to outline a comparison of results between algorithm PCDM and that of when solving QP problems arising from MPC for random network systems. Both algorithms were implemented in the same manner as described in Section 4.2. We considered random network systems with dynamics (11) generated as follows: the entries of system matrices $A^{ij}$ and $B^{ij}$ are taken from a normal distribution with zero mean and unit variance. Matrices $A^{ij}$ are then scaled, so that they become neutrally stable. Matrices $Q^{i}\succeq 0$ and $R^{i}\succ 0$ are random. The input variables are constrained to lie in box sets whose boundaries are generated randomly. The terminal cost matrices $P^{i}$ are taken to be the solution of the SDP problem given in Lemma 2. For each subsystem the number of inputs is taken $m_{i}=5$ or $m_{i}=10$. We let the prediction horizon range between $N=6$ to $N=120$. The subsystems are arranged in a ring, i.e. $\mathcal{N}^{i}=\{i-1,i,i+1\}$. We first considered $M=8$ subsystems, matching the number of cores on our PC. Parallel implementation was also carried out for $M=16$ subsystems, with each core of the PC running two processes. The resulting random QP problems have $p=MNm_{i}$ variables. The stopping criterion for each algorithm is $f(\mathbf{u}_{k})-f^{*}\leq 0.001$, with $f^{*}$ being precomputed for each problem using Matlab’s quadprog. For each prediction horizon, $10$ simulations were run, starting from different random initial states.

Table 4 presents the average CPU time in seconds for the execution of each algorithm. It illustrates that Algorithm PCDM, with its design for distributed computations and simple iterations, usually performs better than that in , where the assumption is that for each iteration, a QP problem of size $p/M$ needs to be solved. The entries with $*$ denote that the algorithm would have taken over $5$ hours to complete. Also note that our implementation of the algorithm from , for problems of larger dimensions, i.e starting with $p=3200$, takes less time for it to complete if the problem is divided between $M=16$ subsystems than $M=8$. This is due to the fact that the solver qpip takes much more time to solve problems of size $600$ in the case of $p=4800$ and $M=8$ than problems of size $300$ for $p=4800$ and $M=16$. Also, the transmission delays between subsystems are negligible in comparison with these qpip times. We have also implemented Algorithm PCDM in a centralized manner, i.e. without using MPI and, as can be seen from the table, we gain speedups of computation when the algorithm is parallelized. Algorithm PCDM is outperformed by Matlab’s quadprog, but do note that quadprog is not designed for distributed implementation and there are no transmission delays between processes.

Conclusions

In this paper we have proposed a parallel optimization algorithm for solving smooth convex problems with separable constraints that may arise e.g in MPC for general linear systems comprised of interconnected subsystems. The new optimization algorithm is based on the block coordinate descent framework but with very simple iteration complexity and using local information. We have shown that for strongly convex objective functions it has linear convergence rate. An MPC scheme based on this optimization algorithm was derived, for which every subsystem in the network can compute feasible and stabilizing control inputs using distributed computations. An analysis for obtaining local terminal costs from a distributed viewpoint was made which guarantees stability of the closed-loop interconnected system. Preliminary numerical tests show that this algorithm is suitable for MPC applications, especially those with hardware that has low computational power.

References