A Revised Mehrotra Predictor-Corrector algorithm for Model Predictive Control

TL;DR

This paper introduces a revised Mehrotra predictor-corrector interior-point algorithm tailored for quadratic programming in input-constrained model predictive control, significantly reducing computational time compared to standard solvers.

Contribution

A novel extension of Mehrotra's algorithm specifically designed for quadratic programming in MPC, enhancing computational efficiency for real-time applications.

Findings

The new algorithm is faster than MATLAB's solver.

It effectively handles input constraints in MPC.

Demonstrated improved computational performance in case studies.

Abstract

Input constrained Model predictive control (MPC) includes an optimization problem which should iteratively be solved at each time-instance. The well-known drawback of model predictive control is the computational cost of the optimization problem. This results in restriction of the application of MPC to systems with slow dynamics, e.g., process control systems and small-scale problems. Therefore, implementing fast numerical optimization algorithms has been a point of interest. Interior-point methods are proved to be appropriate algorithms, from computational cost point-of-vie, to solve input-constrained MPC. In this paper first a modified version of Mehrotra's predictor-corrector algorithm, a famous interior-point algorithm, is extended for quadratic programming problems and then is applied to the constrained model predictive control problems. Results show that as expected, the new algorithm is faster than Matlab solver's algorithm.

Tables

Table 1. Table 1: Mehrotra’s Predictor-Corrector Algorithm for Quadratic Programming

Mehrotra’s Predictor-Corrector Algorithm

for Quadratic Programming

Starting from an appropriate initial point $(x_0,y_0,{\lambda }_0)$

for $k=0,1,2,\dots$

Set $(x,y,\lambda )=(x_k,y_k,{\lambda }_k)$ and solve the following

linear system for $(\Delta x^{\text{aff}},\Delta y^{\text{aff}},\Delta {\lambda }^{\text{aff}})$

Calculate $\mu =y^{𝖳}\lambda /m$

Calculate

${\hat{\alpha }}_{\text{aff}}=\text{max}\left\{\alpha \in (0,1]\right|(y,\lambda )+\alpha (\Delta y^{\text{aff}},\Delta {\lambda }^{\text{aff}})\ge 0\}$

Calculate

${\mu }_{\text{aff}}=\frac{{\left(y+{\hat{\alpha }}^{\text{aff}}\Delta y^{\text{aff}}\right)}^{𝖳}}{m}{\left(\lambda +{\hat{\alpha }}^{\text{aff}}\Delta y^{\text{aff}}\right)}^{𝖳}\left(\lambda +{\hat{\alpha }}^{\text{aff}}\Delta {\lambda }^{\text{aff}}\right)$

set centering parameter to $\sigma ={\left({\mu }_{\text{aff}}/\mu \right)}^3$

solve the following linear system for $(\Delta x,\Delta y,\Delta \lambda )$

Choose ${\tau }_k\in (0,1)$ and set $\hat{\alpha }=\text{min}({\alpha }_{{\tau }_k}^{pri},{\alpha }_{{\tau }_k}^{dual})$

Set $(x_{k+1},y_{k+1},{\lambda }_{k+1})=(x_k,y_k,{\lambda }_k)+\hat{\alpha }(\Delta x,\Delta y,\Delta \lambda )$

end for

Table 2. Table 2: Revised Predictor-Corrector Algorithm for Quadratic Programming. ⋆ = − Λ 𝒴 e − Δ Λ aff Δ 𝒴 aff e + σ μ e \star=-\Lambda\mathcal{Y}e-\Delta\Lambda^{\text{aff}}\Delta\mathcal{Y}^{\text{aff}}e+\sigma\mu e and ∗ = − Λ 𝒴 e − α a Δ Λ aff Δ 𝒴 aff e + σ μ e \ast=-\Lambda\mathcal{Y}e-\alpha_{a}\Delta\Lambda^{\text{aff}}\Delta\mathcal{Y}^{\text{aff}}e+\sigma\mu e and Υ = − Λ ​ 𝒴 ​ e − α a ​ Δ ​ Λ aff ​ Δ ​ 𝒴 aff ​ e + σ ​ μ ​ e Υ Λ 𝒴 𝑒 subscript 𝛼 𝑎 Δ superscript Λ aff Δ superscript 𝒴 aff 𝑒 𝜎 𝜇 𝑒 \Upsilon=-\Lambda\mathcal{Y}e-\alpha_{a}\Delta\Lambda^{\text{aff}}\Delta\mathcal{Y}^{\text{aff}}e+\sigma\mu e .

Revised Predictor-Corrector Algorithm for

Quadratic Programming

Input: parameters $\gamma \in (0,1/4)$, $\epsilon >0$ and the

safeguard parameter $\beta \in [\gamma ,1/4)$;

Starting from an appropriate initial point

$(x_0,y_0,{\lambda }_0)\in {𝒩}_{\infty }^{-}\left(\gamma \right)$

begin

while $y^{𝖳}\lambda \ge \epsilon$ do

begin

Predictor Step

Solve the following linear system and compute

the allowable step size (${\alpha }_a$)

such that $(x\left({\alpha }_a\right),\lambda \left({\alpha }_a\right),y\left({\alpha }_a\right))\in ℱ$;

IF $\left(1-{\alpha }_a\right)y{\left({\alpha }_a\right)}^{𝖳}\lambda \left(\xi \right)\le \epsilon$ then

Let $y=y\left({\alpha }_a\right),\lambda =\lambda \left({\alpha }_a\right)$ and Stop.

end

end

begin

Corrector Step

IF ${\alpha }_a>\xi$, then let ${\alpha }_a=\xi$

end

IF ${\alpha }_a\ge 0.1$, then solve the following linear system

with $\mu ={\left(1-{\alpha }_a\right)}^3/{\mu }_g$ and compute the maximum

allowable ${\alpha }_c$ such that $(x\left({\alpha }_c\right),\lambda \left({\alpha }_c\right),y\left({\alpha }_c\right))\in {𝒩}_{\infty }^{-}$

end

IF , then solve the following linear system

with $\mu$ as defined before and compute the maximum

allowable step such that $(x\left({\alpha }_c\right),\lambda \left({\alpha }_c\right),y\left({\alpha }_c\right))\in {𝒩}_{\infty }^{-};$

end

IF , then solve the following linear system

with $\mu =\frac{\beta }{1-\beta }{\mu }_g$ and compute the maximum

allowable step size

such that $(x\left({\alpha }_c\right),\lambda \left({\alpha }_c\right),y\left({\alpha }_c\right))\in {𝒩}_{\infty }^{-}$

end

Set $(x,\lambda ,y)=(x\left({\alpha }_c\right),\lambda \left({\alpha }_c\right),y\left({\alpha }_c\right))$

end

end

Table 3. Table 3: Comparison between the speed of the proposed algorithm vs. MATLAB’s embedded algorithm. Algorithm 1: Proposed Algorithm, Algorithm 2: MATLAB’s Algorithm

$𝒜$	$ℬ$	Alg.	$N_u=3$	$N_u=10$	$N_u=10$
$\left[\begin{array}{c}1\\ -0.8\end{array}\right]$	$\left[\begin{array}{c}0.4\\ 0.6\end{array}\right]$	1	0.5472	0.6684	0.9543
		2	0.9242	0.963545	1.0217
$\left[\begin{array}{c}1\\ -1\\ -0.8\end{array}\right]$	$\left[\begin{array}{c}0.4\\ 0.6\end{array}\right]$	1	0.5358	0.6684	0.8579
		2	0.9425	0.9635	0.9785
$\left[\begin{array}{c}1\\ -1\\ -0.8\end{array}\right]$	$\left[\begin{array}{c}0.04\\ -6\end{array}\right]$	1	0.5358	0.7886	0.8051
		2	0.8961	0.9234	0.8456
$\left[\begin{array}{c}1\\ -1\\ 0.675\end{array}\right]$	$\left[\begin{array}{c}0.04\\ -6\end{array}\right]$	1	0.7489	0.8845	0.8152
		2	0.9124	0.9232	0.9490

Equations

$$\begin{gathered} J(N_1,N_2,N_u)\text{vbox}:=\sum _{j=N_1}^{N_2}\delta (j)[\hat{y}(t+j|t)-w(t+j){]}^2 \\ +\sum _{j=1}^{N_u}\eta (j)[\Delta u(t+j-1){]}^2 \end{gathered}$$

$$\mathcal{A}(z^{-1})y(t)=\mathcal{B}(z^{-1})z^{-d}u(t-1)+\mathcal{C}(z^{-1})\frac{e(t)}{\Delta }$$

$$\mathcal{A}(z^{-1})$$

$$\mathcal{B}(z^{-1})$$

$$\mathcal{C}(z^{-1})$$

$$\tilde{A}(z^{-1})$$

$$1=E_j(z^{-1})\tilde{A}(z^{-1})+z^{-j}{F}_j(z^{-1}),$$

$$J=\frac{1}{2}{u}^{\mathsf{T}}Gu+c^{\mathsf{T}}u+f_0$$

$$G=2({\Gamma }^{\mathsf{T}}\Gamma +\eta I)$$

$$c^{\mathsf{T}}=2(f-w{)}^{\mathsf{T}}\Gamma$$

$$f_0=(f-w{)}^{\mathsf{T}}(f-w)$$

$$y$$

$$f$$

$$F(z^{-1})$$

$${\Gamma }_j(z^{-1})$$

$$\Gamma$$

$$\Gamma^{{}^{\prime}}(z^{-1})=\\ \left[\begin{matrix}\Big{(}G_{d+1}(z^{-1})-g_{0}\Big{)}z\\ \Big{(}G_{d+2}(z^{-1})-g_{0}-g_{1}z^{-1}\Big{)}z^{2}\\ \vdots\\ \Big{(}G_{d+N}(z^{-1})-g_{0}-g_{1}z^{-1}-\dots-g_{N-1}z^{-(N-1)}\Big{)}z^{N}\end{matrix}\right]$$

$$u$$

$$y$$

$$w$$

$$\underset{x}{\min }$$

$$Ax\ge b$$

$$Gx-A^{\mathsf{T}}\lambda +c=0$$

$$Ax-y-b=0$$

$$y_i{\lambda }_i=0$$

$$(y,\lambda )\ge 0$$

$$\mu =\frac{y^{\mathsf{T}}\lambda }{m}$$

$$-\Lambda \mathcal{Y}e-\Delta {\Lambda }^{\text{aff}}\Delta {\mathcal{Y}}^{\text{aff}}e+\sigma \mu e$$

$${\alpha }_{\tau }^{\text{pri}}$$

$${\alpha }_{\tau }^{\text{dual}}$$

$$\Lambda$$

$$\mathcal{Y}$$

$$e$$

$${\mathcal{N}}_{\infty }^{-}(\gamma )=\{(x,\lambda ,y)\in {\mathcal{F}}^0|y_i{\lambda }_i\ge \gamma \mu \}$$

$$\gamma \in (0,1)$$

$$\mu =\frac{y^{\mathsf{T}}\lambda }{m}$$

$${\mathcal{F}}^0=\{(x,\lambda ,y)|Gx+c-A^{\mathsf{T}}\lambda =0,Ax-y=b,(\lambda ,y)\ge 0\}$$

$$\xi =1-{\left(\frac{2\gamma t}{1-\gamma }\right)}^{\frac{1}{3}}$$

$$t=\underset{i\in {\mathcal{I}}_{+}}{\max }\left(\frac{\Delta y_i^{\text{aff}}\Delta {\lambda }_i^{\text{aff}}}{y_i{\lambda }_i}\right)$$

Full text

1 Introduction

Model predictive control (MPC) is a technique where prediction of the future behavior of a dynamical system and the resulting optimal control action are determined by iteratively solving an online optimization problem at each time-step. At each time interval, model predictive control computes a sequence of future variable adjustments by solving a finite horizon open-loop optimal control problem and using the first element of the input sequence as the control action for the system. Then in the next iteration, all calculations are repeated. Nowadays MPC is applied to wide range of applications from aerospace to automotive and control of chemical processes [26].

MPC has many advantages, including its ability to be applied to multivariable systems and plants where offline computations are not feasible. Another desired property of MPC is its capability to include explicit expression of state and input constraints.

In absence of constraints, the infinite-horizon optimization problem of MPC is the LQR problem. A well-known drawback of MPC is its computational cost which restricts application of MPC to systems with slow dynamics [28, 22, 36] and makes MPC inappropriate for large-scale systems [27, 20, 3]. Even for such systems, computer technology required to perform the computations is sometimes costly and also causes software certification concerns [1]. For such systems, often the optimal control problem in MPC is a convex quadratic programming [28]. Different approaches have been proposed to solve the optimization problem of MPC numerically, based on linearity or nonlinearity of the system and constraints. For example, one of the simplest cases, linear-MPC, is when constraints and the dynamical system are all linear and the optimization problem’s objective function is a quadratic function [2]. Another simple expression of MPC that makes it solvable is to use a hybrid prediction model, and as a result categorizing the optimization problem as a mixed integer quadratic program (MIQP) or a mixed-integer linear program (MILP) [1].

While in general in solving the optimization problem of MPC feedback is computed online using iterative optimization algorithms, in special cases, e.g. explicit MPC [4, 5], feedback can be computed offline. In explicit MPC, a small but dense Hessian is resulted from removing states from the constraints and the objective function. If the horizon length is shown by $N$, number of constrained variables are described by $n$, and number of inputs is shown with $p$, the computational complexity for explicit MPC is $\mathcal{O}((n+p)p^{2}N^{3})$, while the computational complexity of implicit MPC is shown in [28] to be $\mathcal{O}((p+m)^{3}+l(p+m)^{2})N)$ where $m$ is the number of states [32]. Due to memory demand, explicit MPC cannot be used for large-scale applications [4, 1, 20, 16]. Even for small problems, solving the problem online is not always slower [36].

Increasing the speed of solving the optimization problem that we are facing with in MPC, has been a point of interest for many years. One approach is to increase the speed by parallel hardware implementation [17, 21, 35, 37]. Another approach is to increase the speed through improving numerical optimization algorithms to solve the system of equations. In the case of linear MPC, Gradient Method [14, 29, 25], interior-point methods, and active-set methods [15, 34, 6, 11, 10] are examples of numerical algorithms that are implemented to solve MPC.

Interior-point methods (IPM) [19] have been proved as an efficient way to solve large-scale quadratic programming (QP) problems which arise in MPC [34, 32]. The reason is that interior-point algorithms can be characterized by polynomial complexity and are believed to be effective in large-scale problems, and hence MPC. Many IPMs solve Karush-Kahn-Tucker (KKT) use Newton’s method [38] while inexact-IPMs do not need to solve KKT conditions [34].

From implementation point of view, interior-point methods and active-set methods are different. In contrast to active-set methods, interior-point methods cannot easily use warm-start strategies, although there are efforts to yield warm-start IPMs [9, 39, 13, 18]. See [33, 31, 32] for implementation of warm-start interior point methods in model predictive control context. Different interior-point methods are discussed in [38] and it has been resulted that the most practical IPM is Mehrotra’s Predictor-Corrector Primal-Dual interior point method [23, 38, 24].

2 Model Predictive Control

Although the methodology of all model predictive controllers share common strategy but they may be different in system’s model, objective function, and the procedure to obtain the control law. In Generalized Predictive Control (GPC) [7], the objective function is:

[TABLE]

where $N_{1}$ is the minimum costing horizon, $N_{2}$ is maximum costing horizon, $w$ is future reference trajectory, $\hat{y}(t+j|t)$ is prediction of system’s output $j$ step ahead, $\delta$ and $\eta$ are weighting sequences, $N_{u}$ is the control horizon and $u$ is the control sequence. Our goal is to make $w$ and $\hat{y}$ as close as possible by defining the appropriate value for $u$ sequence.

By CARIMA (Auto-Regressive Integrated Moving Average) model, system can be defined as:

[TABLE]

where $\Delta=1-z^{-1}$, and $\mathscr{A}$,$\mathscr{B}$, and $\mathscr{C}$ are

[TABLE]

By Diophantine equation, we have

[TABLE]

where $j$ is the number of steps ahead. The polynomials $E_{j}$ and $F_{j}$ can be obtained by dividing $1$ by $\tilde{A}(a^{-1})$ until the remainder can be described as $z^{-j}F_{j}(z^{-1})$. Having the definitions above, GPC results the following quadratic program

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

Here $u$, $y$, and $w$ are defined as

[TABLE]

3 Modified Predictor Corrector Algorithm

In [30] it is shown that Mehrotra’s algorithm [23] has drawbacks. By a numerical example it is shown that there are cases where, in practice, Mehrotra’s algorithm performs very small steps in order to simultaneously keep the resulting points in a specific neighborhood of the central path while also maintaining the computational complexity of the algorithm polynomial. In [30] a new algorithm with safeguard on step lengths is proposed for linear programming and it is shown that the computational complexity of the resulting algorithm is polynomial-time.

If we define the problem as

[TABLE]

where $G$ is symmetric and positive semidefinite, and matrix $A$ is $m$-by-$n$, then KKT conditions could be written as

[TABLE]

where $y$ is the vector of slack variables and $\lambda$ is the vector of Lagrange multipliers associated with inequality constraints. Complementarity measure is defined by

[TABLE]

The Mehrotra’s algorithm for quadratic programming problems [24] is given in Table 1.

In Table 1,

[TABLE]

and $\sigma\in[0,1]$. Matrices $\Lambda$ and $\mathcal{Y}$ are defined as

[TABLE]

where $\tau\in(0,1)$ and

[TABLE]

The algorithm proposed by Salahi and Terlaki[30] can be extended into the quadratic programming case and the algorithm in Table 2 is resulted. Here

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

4 Simulation

The algorithm in Table 2 is applied to the GPC problem. Different systems all including the inequality constraint ($-0.5\leq u\leq 1$ ) are taken into account and comparison is made between the proposed algorithm and MATLAB’s optimization toolbox version 5.0 embedded algorithm for solving quadratic programming problems using ”quadprog” command which uses a reflective Newton method [8, 12].

The simulation is implemented on a computer with 3.5 GB of RAM and Intel Core2 Due Processor E7300 with clock speed 2.66 GHz. MATLAB software version 7.10.0 was used for the simulation.

Results can be seen in Table 3. Four dynamical systems with different number of poles and zeros are taken into acount. In Figure 1, the output, control sequence, and $\Delta u$ for of one of the dynamical systems is given.

5 Conclusion

In this paper we have proposed a revised version of Mehrotra’s predictor-corrector algorithm for the optimization problem of constrained-linear model predictive control. The algorithm has been developed for quadratic programming since most model predictive control algorithms solve a quadratic program to obtain the control sequence. Numerical examples show that, at least for small problems, the proposed algorithm is faster than Matlab’s algorithm for solving model predictive control problems.