# Fast Model Predictive Control Based on Adaptive Alternating Direction Method of Multipliers

Authors: Yu Li; Qiming Zou; Xiaoru Ji; Chanyuan Zhang; Ke Lu
Journal: Journal of Chemistry
 (2019)

---

## Abstract

Model Predictive Control (MPC) can effectively handle control problem with disturbances, multicontrol variables, and complex constraints and is widely used in various control systems. In MPC, the control input at each time step is obtained by solving an online optimization problem, which will cause a time delay in real time on embedded computers with limited computational resources. In this paper, we utilize adaptive Alternating Direction Method of Multipliers (a-ADMM) to accelerate the solution of MPC. This method adaptively adjusts penalty parameter to balance the value of primal residual and dual residual. The performance of this approach is profiled via the control of a quadcopter with 12 states and 4 controls and prediction horizon ranging from 10 to 40. The simulation results demonstrate that the MPC based on a-ADMM has a significant improvement in real-time and convergence performance and thus is more suitable for solving large-scale optimal control problems.


---

## Body

## 1. Introduction

Model Predictive Control (MPC) [1, 2], a.k.a. receding horizon control, is a classical optimization-based control strategy for multivariable constrained systems. Since its appearance in the 1970s, MPC has been widely applied in a variety of areas, such as the transportation system [3, 4], grid system [5], and water distribution system [6]. MPC only utilizes the first element of the control input sequence obtained by solving an optimization problem. On the one hand, although receding horizon way makes MPC robust to model bias and external noises, MPC requires solving the optimization problem in each control cycle, which may cause a large time delay. On the other hand, low-complexity optimization algorithms have been motivated by the fact that several of these recent applications exploit microcontroller that has limited computing power [7–9].For linear systems with a quadratic cost function and multiple linear constraints, MPC can be easily transformed into a standard quadratic programming (QP) form [10]. Therefore, off-the-shelf QP solvers can be used in this case. Most of them can be categorized into one of the two main classes, i.e., active-set and interior-point methods. Interior-point methods have become common choice for implementation because of its good convergence performance [11, 12]. However, it tends to demand significant online computation effort to compute the input. Recently, some modified methods have been proposed to reduce the computational complexity [13]. For example, Domahidi et al. decreased the computational burden by introducing a low-rank matrix forward substitution scheme and extended the interior-point methods to applications on low-cost embedded hardware. Compared to interior-point methods, active-set methods need on average substantially more iterations, but each iteration is computationally much cheaper and can be ward- or hot-start. For example, a parametric active-set algorithm for quadratic programming (qpOASES) is another commonly used QP solver which is applicable when numerous related QPs are solved sequentially, by exploiting the geometric property of state space [14]. Unfortunately, since qpOASES is an active-set-based method, it is limited to small- to medium-scale QP problem and sensitive to the external disturbance [14].Aforementioned methods are centralized optimization algorithms, which are not appropriate for real-time large-scale applications on low-cost embedded computers. On the contrary, due to its success in solving large-scale optimization problems and pretty low computational complexity, distributed optimization algorithms have drawn increasing attention [15]. As a distributed optimization algorithm, Alternating Direction Method of Multipliers (ADMM) is a strong algorithm for solving convex optimization problems [16, 17]. ADMM was first introduced in the mid-1970s for the numerical solution of partial differential equations [18, 19] and particularly useful for solving optimization problems that are too large to be handled by generic optimization solvers. Because of its fast processing speed and good convergence performance, the method has been used in a large number of areas, such as compressed sensing [20], neural network [21], and image processing [22]. More details can be found in [17]. This wide range of applications has triggered a great interest in developing a better and more comprehensive understanding of the theoretical properties of ADMM.ADMM has been proved to be an excellent fit for QP problems, since it can achieve linear convergence with strong convexity assumptions [23]. However, a drawback of ADMM is that the number of iterations required to converge is highly dependent on the step-size parameters. Unfortunately, analytic results for the optimal choice of penalty parameter are not available except for very special problems [24, 25]. In the special case of a strictly convex quadratic problem, criteria for choosing an optimal constant penalty have been recently proposed in [24, 26]. Although the optimal penalty parameter can significantly improve practical convergence, these methods make strong assumptions about the objective and constraint matrix. By introducing relaxation parameter, a method called relaxed ADMM uses past iterations to compute the current one. Instead of the fixed penalty parameter, Aybat and Iyengar [27] suggested using a penalty parameter sequence. Another different approach to the above methods is to accelerate ADMM by adaptively modifying the penalty parameter. As one of the adaptive ADMM methods, residual balancing (RB) [28] tunes penalty parameter to keep both residuals of similar magnitude. We borrow the idea from this method to further improve the convergence of ADMM combined with MPC.In this paper, we utilize a-ADMM to efficiently solve the MPC problem, and the performance is verified by multiple simulation experiments. To be more specific, we transform the standard MPC problem into a distributed formulation and then take advantage of a-ADMM to solve this problem. We verified the method on an Unmanned Aerial Vehicle (UAV) control problem, and this mechanical control system includes 12 states and 4 control inputs [29]. We tested our method on different settings, i.e., the prediction horizon N ranges from 10 to 40. For each setting, MPC requires the solution of a QP with 12+4×N variables and 12+4×2×N inequality constraints. In this complex UAV control problem, a-ADMM algorithm achieves an absolute improvement of real-time response ability over conventional interior-point method and convergence performance over other ADMM variations.The paper is organized as follows. Section2 briefly introduces MPC, ADMM, and some ADMM variations. Section 3 describes the main procedures of the transformation from MPC to the distributed optimization problem and the application of a-ADMM in this specific problem. In Section 4, we compared the method with some baselines in the UAV control task. Finally, the conclusion is drawn and our future work is declared in Section 5. Appendix A describes how to formulate the QP problem of MPC in some details.


## 2. Backgrounds

### 2.1. Model Predictive Control

MPC is an optimal control approach which can handle multiple constraints on states and control inputs. Using a receding horizon control framework, MPC is robust to the model bias and external noises, which makes MPC fit well for highly complex multivariable industry processes.In this paper, the following discrete-time linear time-invariant system is considered:(1)Xt+1=AXt+BUt,where Xt∈ℝn and Ut∈ℝm denote the state and control input, respectively, at time step t. A∈ℝn×n and B∈ℝn×m are the system and control matrices, respectively. For the regulation problem, at each time step t, MPC solves the finite horizon optimal control problem:(2a)minUXNTPXN+∑t=1N−1XtTQXt+UtTRUt,(2b)s.t.Xt+k+1t=AXt+kt+BUt+kt,k=0,…,N−1,(2c)Xmin≤Xt+kt≤Xmax,k=1,…,N,(2d)Umin≤Ut+kt≤Umax,k=1,…,N,(2e)X0=Xt,where N is the prediction horizon. Xt+kt denotes the predicted state vector at time t+k, obtained by applying the input sequence ut,…,ut+k−1. Xmax∈ℝn and Xmin∈ℝn define constraints on state variable X, Umax∈ℝm and Umin∈ℝm define constraints on input variable U. P∈ℝn×n and Q∈ℝn×n are stage and terminal matrices which penalize the state deviation at the end of the prediction horizon N and over the entire horizon, respectively. R∈ℝm×m is the cost matrix for the control inputs. Typically, P≽0, Q≽0, and R≻0.In MPC, control inputs are obtained by solving open-loop optimal control problem (2a)–(2e) [10]. To be more specific, at each time step t, the optimization problem (2a)–(2e) is solved by taking the current state Xt as initial state of the problem and generates a sequence of optimal control input. However, according to receding horizon framework of MPC, only the first element of the control input sequence is applied to the plant. The remaining solutions are culled, and optimal control problem with new initial state Xt+1 will be solved again at the next time step t+1.


### 2.2. Alternating Direction Method of Multipliers

ADMM is a strong algorithm which can efficiently solve convex optimization problems, and it can decompose the global problem into several smaller and easy-to-solve local subproblems [30]. Because of its fast processing speed and good convergence performance, ADMM has received tremendous interest and approval for solving numerous problems in machine learning, statistics, and signal processing [17, 31–33].ADMM can efficiently solve the optimization problems of the following standard form:(3a)minUfU,(3b)s.t.U∈C,where f and C are convex. This problem can be rewritten in another form as(4a)minU,ZfU+gZ,(4b)s.t.U−Z=0,where g is the indicator function of C, U is primal variable, and Z is separate variable. The augmented Lagrangian is defined as(5)LρU,Z,c=fU+gZ+ρ2U−Z+c22,ρ>0,where ρ is the step size (or penalty parameter) and c is the scaled dual variable. The general augmented Lagrange multiplier method is to minimize Lρ, w.r.t. U and Z at the same time. This is in difficulty and does not exploit the fact that the objective function is separable. To remedy this issue, ADMM decomposes the global problem (4) into two subproblems and minimizes U and Z, respectively, with the following three steps:(6a)Uk+1≔argminfU+ρ2U−Zk+ck22,(6b)Zk+1≔argmingZ+ρ2Uk+1−Z+ck22,(6c)ck+1≔ck+Uk+1−Zk+1.Under rather mild conditions, ADMM converges for any constant penalty parameter [34]. The convergence of ADMM is based on the primal and dual residuals:(7)rk=Uk−Zk,sk=−ρZk−Zk−1.It has been observed that the algorithm approaches to the optimal solution when the primal and dual residuals approach to zero. Generally, the iterative process is terminated if(8a)rk2≤pεabs+εrelmaxUk2,Zk2,(8b)sk2≤nεabs+εrelcρk2,where εabs>0 is absolute tolerance and εrel>0 is relative tolerance. However, ρ has an immediate effect on the convergence element of the algorithm, and inappropriate tuning of this algorithm parameter may render the method moderative.One way of accelerating the convergence properties of the algorithm is to utilize previous iterates when computing the following ones. This practical method is called relaxed Alternating Direction Method of Multipliers (r-ADMM) [35] that replaces Uk+1 with αkUk+1−1−αkZk in the Z− and c− updates (6a)–(6c), yielding(9a)Uk+1≔argminfU+ρ2U−Zk+ck22,(9b)U˜k+1≔αkUk+1−1−αkZk,(9c)Zk+1≔argmingZ+ρ2Uk+1−Z+ck22,(9d)ck+1≔ck+Uk+1−Zk+1.The parameterαk∈0,2 is called the relaxation parameter. Note that letting αk=1 for all k recovers the original ADMM iterations. Empirical studies have illustrated that overrelaxation, i.e., letting αk>1, is often beneficial and the guideline αk∈1.5,1.8 has been proposed [36].Another approach for improving the convergence properties of the algorithm is to utilize different penalty parameters for each iteration. Adaptive Alternating Direction Method of Multipliers (a-ADMM) [28] is based on the following observation: increasing ρk strengthens the penalty term and leads to smaller primal residuals but larger dual ones; on the contrary, decreasing ρk results in smaller dual residuals and larger primal residuals. As both residuals must be small at convergence, it requires adjusting ρ to keep both primal residuals and dual residuals of the same magnitude. A simple and effective scheme for this purpose is as follows:(10)ρk+1≔τincrρk,ifrk2>μsk2,ρkτdecr,ifsk2>μrk2,ρk,otherwise,where μ > 1, τincr>1, and τdecr>1. Using different penalty parameter ρk for each iteration k, a-ADMM can converge quickly and be insensitive to the initial penalty parameter. However, without a careful choice of μ, τincr, and τdecr, this algorithm may fail to converge [28].


## 2.1. Model Predictive Control

MPC is an optimal control approach which can handle multiple constraints on states and control inputs. Using a receding horizon control framework, MPC is robust to the model bias and external noises, which makes MPC fit well for highly complex multivariable industry processes.In this paper, the following discrete-time linear time-invariant system is considered:(1)Xt+1=AXt+BUt,where Xt∈ℝn and Ut∈ℝm denote the state and control input, respectively, at time step t. A∈ℝn×n and B∈ℝn×m are the system and control matrices, respectively. For the regulation problem, at each time step t, MPC solves the finite horizon optimal control problem:(2a)minUXNTPXN+∑t=1N−1XtTQXt+UtTRUt,(2b)s.t.Xt+k+1t=AXt+kt+BUt+kt,k=0,…,N−1,(2c)Xmin≤Xt+kt≤Xmax,k=1,…,N,(2d)Umin≤Ut+kt≤Umax,k=1,…,N,(2e)X0=Xt,where N is the prediction horizon. Xt+kt denotes the predicted state vector at time t+k, obtained by applying the input sequence ut,…,ut+k−1. Xmax∈ℝn and Xmin∈ℝn define constraints on state variable X, Umax∈ℝm and Umin∈ℝm define constraints on input variable U. P∈ℝn×n and Q∈ℝn×n are stage and terminal matrices which penalize the state deviation at the end of the prediction horizon N and over the entire horizon, respectively. R∈ℝm×m is the cost matrix for the control inputs. Typically, P≽0, Q≽0, and R≻0.In MPC, control inputs are obtained by solving open-loop optimal control problem (2a)–(2e) [10]. To be more specific, at each time step t, the optimization problem (2a)–(2e) is solved by taking the current state Xt as initial state of the problem and generates a sequence of optimal control input. However, according to receding horizon framework of MPC, only the first element of the control input sequence is applied to the plant. The remaining solutions are culled, and optimal control problem with new initial state Xt+1 will be solved again at the next time step t+1.


## 2.2. Alternating Direction Method of Multipliers

ADMM is a strong algorithm which can efficiently solve convex optimization problems, and it can decompose the global problem into several smaller and easy-to-solve local subproblems [30]. Because of its fast processing speed and good convergence performance, ADMM has received tremendous interest and approval for solving numerous problems in machine learning, statistics, and signal processing [17, 31–33].ADMM can efficiently solve the optimization problems of the following standard form:(3a)minUfU,(3b)s.t.U∈C,where f and C are convex. This problem can be rewritten in another form as(4a)minU,ZfU+gZ,(4b)s.t.U−Z=0,where g is the indicator function of C, U is primal variable, and Z is separate variable. The augmented Lagrangian is defined as(5)LρU,Z,c=fU+gZ+ρ2U−Z+c22,ρ>0,where ρ is the step size (or penalty parameter) and c is the scaled dual variable. The general augmented Lagrange multiplier method is to minimize Lρ, w.r.t. U and Z at the same time. This is in difficulty and does not exploit the fact that the objective function is separable. To remedy this issue, ADMM decomposes the global problem (4) into two subproblems and minimizes U and Z, respectively, with the following three steps:(6a)Uk+1≔argminfU+ρ2U−Zk+ck22,(6b)Zk+1≔argmingZ+ρ2Uk+1−Z+ck22,(6c)ck+1≔ck+Uk+1−Zk+1.Under rather mild conditions, ADMM converges for any constant penalty parameter [34]. The convergence of ADMM is based on the primal and dual residuals:(7)rk=Uk−Zk,sk=−ρZk−Zk−1.It has been observed that the algorithm approaches to the optimal solution when the primal and dual residuals approach to zero. Generally, the iterative process is terminated if(8a)rk2≤pεabs+εrelmaxUk2,Zk2,(8b)sk2≤nεabs+εrelcρk2,where εabs>0 is absolute tolerance and εrel>0 is relative tolerance. However, ρ has an immediate effect on the convergence element of the algorithm, and inappropriate tuning of this algorithm parameter may render the method moderative.One way of accelerating the convergence properties of the algorithm is to utilize previous iterates when computing the following ones. This practical method is called relaxed Alternating Direction Method of Multipliers (r-ADMM) [35] that replaces Uk+1 with αkUk+1−1−αkZk in the Z− and c− updates (6a)–(6c), yielding(9a)Uk+1≔argminfU+ρ2U−Zk+ck22,(9b)U˜k+1≔αkUk+1−1−αkZk,(9c)Zk+1≔argmingZ+ρ2Uk+1−Z+ck22,(9d)ck+1≔ck+Uk+1−Zk+1.The parameterαk∈0,2 is called the relaxation parameter. Note that letting αk=1 for all k recovers the original ADMM iterations. Empirical studies have illustrated that overrelaxation, i.e., letting αk>1, is often beneficial and the guideline αk∈1.5,1.8 has been proposed [36].Another approach for improving the convergence properties of the algorithm is to utilize different penalty parameters for each iteration. Adaptive Alternating Direction Method of Multipliers (a-ADMM) [28] is based on the following observation: increasing ρk strengthens the penalty term and leads to smaller primal residuals but larger dual ones; on the contrary, decreasing ρk results in smaller dual residuals and larger primal residuals. As both residuals must be small at convergence, it requires adjusting ρ to keep both primal residuals and dual residuals of the same magnitude. A simple and effective scheme for this purpose is as follows:(10)ρk+1≔τincrρk,ifrk2>μsk2,ρkτdecr,ifsk2>μrk2,ρk,otherwise,where μ > 1, τincr>1, and τdecr>1. Using different penalty parameter ρk for each iteration k, a-ADMM can converge quickly and be insensitive to the initial penalty parameter. However, without a careful choice of μ, τincr, and τdecr, this algorithm may fail to converge [28].


## 3. Fast MPC Based on Adaptive ADMM

### 3.1. Adaptive ADMM MPC

We focus on improving MPC by solving its QP problem with a-ADMM. In this work, the standard MPC problem is transformed into an equivalent QP problem which could be easily decomposed into several smaller local subproblems. And a-ADMM could be finally adopted to efficiently solve the resulting problem.MPC problem ((2a)–(2e)) can be expressed as a QP problem when the system dynamic and constraints are linear and the cost function is quadratic [10]. Specifically, by substituting Xt+kt=AkXt+∑j=0k−1AjBUt+k−1−j into the cost function, the optimization problem equations (2a)–(2e) can be formulated as(11a)minUUTHU+XtTCTU+XtTYXt,(11b)s.t.GU≤W+SXt,where U=UtT,Ut+1T,…,Ut+N−1TT is the optimization vector and H, C, Y, G, W, and S are constant matrices that can be easily obtained offline. Note that H has to be positive definite for the convexity of this problem. The detailed process of this transformation is shown in Appendix A.As the problem is determined by the current stateXt, the implementation of MPC requires the online solution of problem ((11a) and (11b)) at each time step. Although off-the-shelf QP solvers, e.g., active-set methods [37] and interior-point methods [11], are available, it may demand significant computation effort in order to compute the input Ut online.Fortunately, the QP problem (11a) and (11b) can be further transformed into ADMM standard form by introducing a split vector Z named slack vector (see [17]):(12a)minU,ZfU+gZ,(12b)s.t.U−Z=0,where fU=UTHU+XtTCTU+XtTYXt is the initial objective function with qualified domain and gZ=IZ is defined as the indicator function:(13)IZ≔0,ifGZ≤W+SXt,+∞,otherwise.The augmented Lagrangian is defined as(14)LρU,Z,c=fU+gZ+ρ2U−Z+c22,ρ>0,where ρ is the penalty parameter and c is the dual parameter. We propose to solve problem (11a) and (11b) with a-ADMM as follows:(15a)Uk+1≔argminUTHU+XtTCTU+ρ2U−Zk+Ck22,(15b)Zk+1≔argminIZ+ρ2Uk+1−Z+Ck22,(15c)Ck+1≔Ck+Uk+1−Zk+1,where the penalty parameter ρ is a piecewise linear function as equation (10). As the optimizer U is only needed, the term involving Y is usually removed from problem (11a) and (11b). Since both residuals must be small at convergence, ρ is tuned to keep both residuals of similar magnitude.


### 3.2. Convergence Analysis of Adaptive ADMM-MPC

He et al. [28] proved that convergence is guaranteed for a-ADMM when either of the two following conditions is satisfied:Condition 1 (bounded increasing):(16)∑k=1∞ηk2<∞,whereηk=maxρk/ρk−1,1−1.Condition 2 (bounded decreasing):(17)∑k=1∞θk2<∞,where θk=maxρk−1/ρk,1−1.Condition 1 (Condition 2) suggests that the increasement (decreasement) of adaptive penalty parameter is bounded.Based on the above analysis, the framework of a-ADMM algorithm is portrayed in Algorithm1.Algorithm 1:a-ADMM algorithm framework.
a-ADMM algorithm
Initialization: optimization variables U0, Z0, C0
While not converge by stopping criteria ((8a) and (8b)) and k¡maxiter do
Uk+1≔argminUTHU+XtTCTU+ρ/2U−Zk+Ck22
Zk+1≔argminIZ+ρ/2Uk+1−Z+Ck22
Ck+1≔Ck+Uk+1−Zk+1
if rk2>sk2
ρk+1=τincrρk
else if sk2>rk2
ρk+1=ρk/τdecr
else ρk+1=ρk
end if
k=k+1
end while
Output: optimal control input U∗


## 3.1. Adaptive ADMM MPC

We focus on improving MPC by solving its QP problem with a-ADMM. In this work, the standard MPC problem is transformed into an equivalent QP problem which could be easily decomposed into several smaller local subproblems. And a-ADMM could be finally adopted to efficiently solve the resulting problem.MPC problem ((2a)–(2e)) can be expressed as a QP problem when the system dynamic and constraints are linear and the cost function is quadratic [10]. Specifically, by substituting Xt+kt=AkXt+∑j=0k−1AjBUt+k−1−j into the cost function, the optimization problem equations (2a)–(2e) can be formulated as(11a)minUUTHU+XtTCTU+XtTYXt,(11b)s.t.GU≤W+SXt,where U=UtT,Ut+1T,…,Ut+N−1TT is the optimization vector and H, C, Y, G, W, and S are constant matrices that can be easily obtained offline. Note that H has to be positive definite for the convexity of this problem. The detailed process of this transformation is shown in Appendix A.As the problem is determined by the current stateXt, the implementation of MPC requires the online solution of problem ((11a) and (11b)) at each time step. Although off-the-shelf QP solvers, e.g., active-set methods [37] and interior-point methods [11], are available, it may demand significant computation effort in order to compute the input Ut online.Fortunately, the QP problem (11a) and (11b) can be further transformed into ADMM standard form by introducing a split vector Z named slack vector (see [17]):(12a)minU,ZfU+gZ,(12b)s.t.U−Z=0,where fU=UTHU+XtTCTU+XtTYXt is the initial objective function with qualified domain and gZ=IZ is defined as the indicator function:(13)IZ≔0,ifGZ≤W+SXt,+∞,otherwise.The augmented Lagrangian is defined as(14)LρU,Z,c=fU+gZ+ρ2U−Z+c22,ρ>0,where ρ is the penalty parameter and c is the dual parameter. We propose to solve problem (11a) and (11b) with a-ADMM as follows:(15a)Uk+1≔argminUTHU+XtTCTU+ρ2U−Zk+Ck22,(15b)Zk+1≔argminIZ+ρ2Uk+1−Z+Ck22,(15c)Ck+1≔Ck+Uk+1−Zk+1,where the penalty parameter ρ is a piecewise linear function as equation (10). As the optimizer U is only needed, the term involving Y is usually removed from problem (11a) and (11b). Since both residuals must be small at convergence, ρ is tuned to keep both residuals of similar magnitude.


## 3.2. Convergence Analysis of Adaptive ADMM-MPC

He et al. [28] proved that convergence is guaranteed for a-ADMM when either of the two following conditions is satisfied:Condition 1 (bounded increasing):(16)∑k=1∞ηk2<∞,whereηk=maxρk/ρk−1,1−1.Condition 2 (bounded decreasing):(17)∑k=1∞θk2<∞,where θk=maxρk−1/ρk,1−1.Condition 1 (Condition 2) suggests that the increasement (decreasement) of adaptive penalty parameter is bounded.Based on the above analysis, the framework of a-ADMM algorithm is portrayed in Algorithm1.Algorithm 1:a-ADMM algorithm framework.
a-ADMM algorithm
Initialization: optimization variables U0, Z0, C0
While not converge by stopping criteria ((8a) and (8b)) and k¡maxiter do
Uk+1≔argminUTHU+XtTCTU+ρ/2U−Zk+Ck22
Zk+1≔argminIZ+ρ/2Uk+1−Z+Ck22
Ck+1≔Ck+Uk+1−Zk+1
if rk2>sk2
ρk+1=τincrρk
else if sk2>rk2
ρk+1=ρk/τdecr
else ρk+1=ρk
end if
k=k+1
end while
Output: optimal control input U∗


## 4. Simulation Experiment

We evaluated our approach on an UAV control task. The dynamic model of UAV used in this experiment is from [29] and can be expressed as follows:(18a)mx¨=−usinθ−βx˙,(18b)my¨=ucosθsinϕ−βy˙,(18c)mz¨=ucosθcosϕ−mg−βz˙,(18d)θ¨=t˜θ,(18e)ϕ¨=t˜ϕ,(18f)ψ¨=t˜ψ,where x, y, and z denote the position of the UAV and θ, ϕ, and ψ denote the rotation angle of the UAV around the Cartesian coordinate axis. The damping coefficient β takes into account the actual friction effect. This dynamic model has 12 states X∈ℝ12 and 4 control inputs U∈ℝ4. For the purpose of designing a appropriate MPC controller for UAV, the nonlinear dynamic model ((18a)–(18f)) is linearized at the equilibrium point (see [29]):(19)Xt+1=AXt+BUt,Yt=CXt,where Xt=θ,ϕ,ψ,x,y,z,θ˙,ϕ˙,ψ˙,x˙,y˙,z˙′∈ℝ12 is the state vector, Ut=u,t˜θ,t˜ϕ,t˜ψ∈ℝ4 is the input vector, and Yt∈ℝ12 is the output vector. Assuming that the UAV system is completely measurable, and the observation matrix C is a 12 × 12 identity matrix. The coefficient matrices A and B are as follows:(20)A=1000000.10000001000000.10000001000000.10000.0488001000.0016000.0992000−0.048800100−0.0016000.09920000001000000.09920000001000000000000100000000000010000.9734000000.0488000.9846000−0.973400000−0.0488000.98460000000000000.9846,B=0−0.072600.0726−0.072600.07260−0.01520.0152−0.01520.01520−0.0006−0.00000.00060.00060−0.000600.01060.01060.01060.01060−1.451201.4512−1.451201.45120−0.30490.3049−0.30490.30490−0.023600.02360.02360−0.023600.21070.21070.21070.2107.The baseline methods we compared in our experiments are as follows.


### 4.1. Interior-Point Method (IPM)

IPM is the mostly used QP solver modeling the problem constraints as parametrized penalty functions which also referred to as barrier functions. At each iteration, an unconstrained optimization problem is solved for varying barrier function until the optimum is achieved.


### 4.2. qpOASES: a Parametric Active-Set Algorithm for Quadratic Programming

qpOASES is a very used QP solver which is applicable when numerous related QPs are solved sequentially, by exploiting the geometric property of state space [14].


### 4.3. Alternating Direction Method of Multipliers (ADMM)

ADMM is a strong algorithm for solving convex optimization problems [16, 17] which decomposes the global problem into several smaller and easy-to-solve local subproblems [30].


### 4.4. Relaxed Alternating Direction Method of Multipliers (r-ADMM)

r-ADMM is one of the variants of ADMM which introduces relaxation parameterα to relax the update of the optimal solution. Empirical studies have illustrated that overrelaxation is beneficial to the convergence [36].We first fixed the prediction horizonN=10 to find the optimal setting of ADMM and r-ADMM: the penalty parameter ρ and the relaxation parameter α. Figure 1 gives the results of different values of parameters. As can be seen, the average number of iterations depends significantly on the penalty parameter ρ and the relaxation parameter α. Based on this result, we simply set ρ=0.7 and α=1.5 in the following experiments.Figure 1
Effect of penalty parameterρ of ADMM (a) and relaxation parameter α of r-ADMM (b) in terms of average number of iterations taken to converge. For each setting, algorithms are performed 50 times on the UAV control task (N = 10) with random initial states.
(a)
(b)Similarly, we used the same control task to find the optimal setting of a-ADMM. According to the theory of a-ADMM, the relaxation parameterα is fixed and can be set α = 1.5 based on previous results. However, different from ADMM, the penalty parameter ρ of a-ADMM is a piecewise linear function as equation (10) influenced by three other parameters, i.e., μ, τincr, and τdecr. Figure 2 gives the results of different values of μ, τincr, and τdecr on the UAV control task when the prediction horizon N=10. Based on the lowest average iterations of each setting, we set μ=10, τincr=2, and τdecr=10 in the following experiments, and the optimization variables U0, Z0, and C0 are all initialized to a zero vector, so that the initial strategy of optimization variables does not affect the performance comparison of these algorithms (IPM, qpOASES, ADMM, r-ADMM, and a-ADMM). The parameters of these algorithms are listed in Table 1.Figure 2
Effect ofμ (a), τincr, (b) and τdecr (c) in terms of average number of iterations taken to converge. For each setting, a-ADMM is performed 50 times on the UAV control task (N = 10) with random initial states.
(a)
(b)
(c)Table 1
List of algorithmic parameters.
Parameters
IPM
qpOASES
ADMM
r-ADMM
a-ADMM
ρ
—
—
0.7
0.7
Adaptive
U
0
0
0
0
0
0
Z
0
0
0
0
0
0
C
0
0
0
0
0
0
α
—
—
1
1.5
1
ε
abs
—
—
1e − 2
1e − 2
1e − 2
ε
rel
—
—
1e − 4
1e − 4
1e − 4
μ
—
—
—
—
10
τ
incr
—
—
—
—
2
τ
decr
—
—
—
—
10
Maxiter
1000
1000
1000
1000
1000The feedback loop controls the quadcopter to land on a platform located at coordinate position0,0,0. As shown in Figure 3, the closed-loop performance (including the feasibility and stability) of ADMM and a-ADMM is exactly the same. The reason is that, with the same absolute tolerance εabs and relative tolerance εrel, ADMM and a-ADMM converge to the same optimal solution leading to the same closed-loop performance.Figure 3
Closed-loop performance (running time ony-axis) versus number of time steps (time steps on x-axis) of ADMM and a-ADMM on UAV control task with N = 10.
(a)
(b)
(c)
(d)Figure4 gives our real-time performance of each method on the UAV control task. This figure shows the cumulative running time of our method and other baselines with the prediction horizon N ranging from 10 to 40. From the simulation results, it can be seen that, combined with MPC, ADMM and its variants significantly outperform custom QP solver, IPM, and qpOASES. With the increment of the task complexity, the running time of IPM and qpOASES increases obviously, while that of ADMM, r-ADMM, and a-ADMM basically remains unchanged.Figure 4
Real-time performance (running time ony-axis) versus number of time steps (time steps on x-axis) of IPM, qpOASES, ADMM, r-ADMM, and a-ADMM on UAV control task with N ranging from 10 to 40: (a) N = 10, (b) N = 20, (c) N = 30, and (d) N = 40.
(a)
(b)
(c)
(d)As shown in Figure4, there is not much difference in the cumulative running time of several variants of the ADMM algorithm. In this part, we aim to further investigate the difference in the cumulative number of iterations between them. Figure 5 portrays the cumulative iteration number of each method when the prediction horizon N ranging from 10 to 40. Apparently, although ADMM and its variations have the similar real-time performance, a-ADMM significantly outperforms all other baselines in terms of convergence performance under all the task settings. qpOASES proposes to move on a straight line in the state space when transitioning from the “old” QP problem to the “new” one. The reason why the cumulative iteration number of qpOASES stays stable after several time steps is mainly because the state of quadcopter has converged to the objective region, and the adjacent QP problem is quite similar.Figure 5
Convergence performance (cumulative number of iterations ony-axis) versus number of time steps (time steps on x-axis) of IPM, qpOASES, ADMM, r-ADMM, and a-ADMM on UAV control task with N ranging from 10 to 40: (a) N = 10, (b) N = 20, (c) N = 30, and (d) N = 40.
(a)
(b)
(c)
(d)


## 4.1. Interior-Point Method (IPM)

IPM is the mostly used QP solver modeling the problem constraints as parametrized penalty functions which also referred to as barrier functions. At each iteration, an unconstrained optimization problem is solved for varying barrier function until the optimum is achieved.


## 4.2. qpOASES: a Parametric Active-Set Algorithm for Quadratic Programming

qpOASES is a very used QP solver which is applicable when numerous related QPs are solved sequentially, by exploiting the geometric property of state space [14].


## 4.3. Alternating Direction Method of Multipliers (ADMM)

ADMM is a strong algorithm for solving convex optimization problems [16, 17] which decomposes the global problem into several smaller and easy-to-solve local subproblems [30].


## 4.4. Relaxed Alternating Direction Method of Multipliers (r-ADMM)

r-ADMM is one of the variants of ADMM which introduces relaxation parameterα to relax the update of the optimal solution. Empirical studies have illustrated that overrelaxation is beneficial to the convergence [36].We first fixed the prediction horizonN=10 to find the optimal setting of ADMM and r-ADMM: the penalty parameter ρ and the relaxation parameter α. Figure 1 gives the results of different values of parameters. As can be seen, the average number of iterations depends significantly on the penalty parameter ρ and the relaxation parameter α. Based on this result, we simply set ρ=0.7 and α=1.5 in the following experiments.Figure 1
Effect of penalty parameterρ of ADMM (a) and relaxation parameter α of r-ADMM (b) in terms of average number of iterations taken to converge. For each setting, algorithms are performed 50 times on the UAV control task (N = 10) with random initial states.
(a)
(b)Similarly, we used the same control task to find the optimal setting of a-ADMM. According to the theory of a-ADMM, the relaxation parameterα is fixed and can be set α = 1.5 based on previous results. However, different from ADMM, the penalty parameter ρ of a-ADMM is a piecewise linear function as equation (10) influenced by three other parameters, i.e., μ, τincr, and τdecr. Figure 2 gives the results of different values of μ, τincr, and τdecr on the UAV control task when the prediction horizon N=10. Based on the lowest average iterations of each setting, we set μ=10, τincr=2, and τdecr=10 in the following experiments, and the optimization variables U0, Z0, and C0 are all initialized to a zero vector, so that the initial strategy of optimization variables does not affect the performance comparison of these algorithms (IPM, qpOASES, ADMM, r-ADMM, and a-ADMM). The parameters of these algorithms are listed in Table 1.Figure 2
Effect ofμ (a), τincr, (b) and τdecr (c) in terms of average number of iterations taken to converge. For each setting, a-ADMM is performed 50 times on the UAV control task (N = 10) with random initial states.
(a)
(b)
(c)Table 1
List of algorithmic parameters.
Parameters
IPM
qpOASES
ADMM
r-ADMM
a-ADMM
ρ
—
—
0.7
0.7
Adaptive
U
0
0
0
0
0
0
Z
0
0
0
0
0
0
C
0
0
0
0
0
0
α
—
—
1
1.5
1
ε
abs
—
—
1e − 2
1e − 2
1e − 2
ε
rel
—
—
1e − 4
1e − 4
1e − 4
μ
—
—
—
—
10
τ
incr
—
—
—
—
2
τ
decr
—
—
—
—
10
Maxiter
1000
1000
1000
1000
1000The feedback loop controls the quadcopter to land on a platform located at coordinate position0,0,0. As shown in Figure 3, the closed-loop performance (including the feasibility and stability) of ADMM and a-ADMM is exactly the same. The reason is that, with the same absolute tolerance εabs and relative tolerance εrel, ADMM and a-ADMM converge to the same optimal solution leading to the same closed-loop performance.Figure 3
Closed-loop performance (running time ony-axis) versus number of time steps (time steps on x-axis) of ADMM and a-ADMM on UAV control task with N = 10.
(a)
(b)
(c)
(d)Figure4 gives our real-time performance of each method on the UAV control task. This figure shows the cumulative running time of our method and other baselines with the prediction horizon N ranging from 10 to 40. From the simulation results, it can be seen that, combined with MPC, ADMM and its variants significantly outperform custom QP solver, IPM, and qpOASES. With the increment of the task complexity, the running time of IPM and qpOASES increases obviously, while that of ADMM, r-ADMM, and a-ADMM basically remains unchanged.Figure 4
Real-time performance (running time ony-axis) versus number of time steps (time steps on x-axis) of IPM, qpOASES, ADMM, r-ADMM, and a-ADMM on UAV control task with N ranging from 10 to 40: (a) N = 10, (b) N = 20, (c) N = 30, and (d) N = 40.
(a)
(b)
(c)
(d)As shown in Figure4, there is not much difference in the cumulative running time of several variants of the ADMM algorithm. In this part, we aim to further investigate the difference in the cumulative number of iterations between them. Figure 5 portrays the cumulative iteration number of each method when the prediction horizon N ranging from 10 to 40. Apparently, although ADMM and its variations have the similar real-time performance, a-ADMM significantly outperforms all other baselines in terms of convergence performance under all the task settings. qpOASES proposes to move on a straight line in the state space when transitioning from the “old” QP problem to the “new” one. The reason why the cumulative iteration number of qpOASES stays stable after several time steps is mainly because the state of quadcopter has converged to the objective region, and the adjacent QP problem is quite similar.Figure 5
Convergence performance (cumulative number of iterations ony-axis) versus number of time steps (time steps on x-axis) of IPM, qpOASES, ADMM, r-ADMM, and a-ADMM on UAV control task with N ranging from 10 to 40: (a) N = 10, (b) N = 20, (c) N = 30, and (d) N = 40.
(a)
(b)
(c)
(d)


## 5. Conclusions

In this paper, we utilize a-ADMM to accelerate the solution of MPC. We firstly transform the standard MPC problem into a distributed formulation and then take advantage of a-ADMM to solve this problem. We verified this method on an UAV control problem and compared a-ADMM with IPM, ADMM, and r-ADMM under different control settings, i.e., the prediction horizonN ranging from 10 to 40. The experimental results have shown that a-ADMM can greatly improve the effectiveness of MPC. As a future and noteworthy research direction, one might provide better scheme for adaptively adjusting penalty parameter considering more complex constraints.



---

*Source: 8035204-2019-09-30.xml*