Model predictive control for ARMAX processes with additive outlier noise

Abstract

The Autoregressive Moving Average (ARMAX) model with exogenous input is a widely used discrete time series model, but its special structure allows outliers of its process to affect multiple output data items, thereby significantly affecting the output. In this paper, a regularized model predictive control (MPC) is proposed for an ARMAX process affected by outlier noise. The outlier noise is modeled as an auxiliary variable in the ARMAX model, and the MPC cost function is reconstructed to reduce the influence of outlier noise on multiple data items. The stability of the proposed method and the convergence of output/input and state are guaranteed. The degree to which regularization affects the system can be adjusted by an optional parameter. This paper provides some helpful insights on how to choose this optional parameter in the cost function. The effectiveness of the proposed method is demonstrated by the results of 200 repeated simulations.

Introduction

The autoregressive moving average (ARMAX) model ¹ with exogenous input is a discrete time series model that is widely used in system description and control design.^2–4 The output of the ARMAX process is an autoregressive structure influenced by a sequence of historical inputs, outputs and white noise, so the model can accurately describe a large number of real processes of general complexity with a good balance between complexity and performance, ⁵ but this structure allows outliers of the ARMAX process to affect multiple data items and thus significantly affect the output, which makes estimating the parameters of the ARMAX process challenging, which has become a research hotspot in processing estimation problems. ⁶

Model predictive control (MPC) is usually used to solve the optimization problem of the ARMAX process.^7,8 At each sampling instant, the open-loop constrained optimal control problem of MPC is solved in a finite horizon, and the control input sequence is obtained. The first data from the obtained optimal control sequence is applied to the system and the problem is repeated at the next sampling instant according to the new measurements.^9,10 In recent years, the MPC algorithm based on the neural network model^11,12 has also received attention. General MPC methods are accurate and undisturbed for system models, but they are inadequate to deal with the various uncertainties of the actual systems, such as model inaccuracy, time variations, and unknown disturbance. Robust model predictive control (RMPC)^13,14 is a model-based control technology, which is suitable for uncertain or nonlinear systems subject to physical constraints.^15,16 By considering the influence of model uncertainties in the prediction horizon and solving the open-loop optimal problem at each sampling time, this method ensures the closed-loop stability and makes the uncertainties meet the constraints.^17,18 Generally when studying control algorithms for models with additional uncertain disturbances, a high-priority option is to improve the performance of RMPC.^19,20

In machine learning and control engineering, outliers are generally regarded as anomalies, inconsistencies, and deviations in the sequence in data statistics studies.^21,22 Although outliers rarely appear, they will affect the actual project analysis and cause deviations in the results. ²³ Therefore, they cannot be ignored in the interference input process. Outliers may occur in data sets due to sensor failures, manual entry errors, or abnormal events. ²⁴ Regularization is an effective solution to the problem that the output is affected by outliers in the ARMAX process. The authors Forero et al. ²⁵ proposed a block coordinate descent method, and obtained an iterative algorithm with guaranteed convergence and less computational complexity. The authors Mateos and Giannakis ²⁶ developed a robust principal component analysis (PCA) approach that encourages sparsity in the outlier matrix by explicitly modeling outliers. The authors Farahmand et al. ²⁷ developed robust smoothers for dynamical processes that are polluted by outliers in both measurements and state dynamics. $l_1$ -norm and $l_2$ -norm are used in Yin et al. ²⁸ and Yin and Gao ²⁹ respectively, to estimate the state of outliers in the ARMAX process. Most of these works introduce additional regularization terms in the optimization function, which can be seen as adding soft constraints to the state. The above literatures have proposed to effectively solve the problem of outlier noise pollution, but the penalty for outlier vectors is relatively conservative, and the value of the regularization parameter is relatively vague.

In order to solve the problem that the ARMAX process is polluted by outlier noise and obtain a clear parameter selection decision. In this paper, we propose an MPC method based on $l_2$ -norm regularization to solve the control problem of ARMAX processes in which the system state is polluted by outliers. We focus on dealing with uncertainties in the state and explicitly modeling outliers as auxiliary variables. Similar to the practice in many previous studies, we modify the cost function by adding extra regularization terms to it. The final cost function is realized by minimizing a quadratic cost function composed of four terms. The first one tracks the state affected by the outliers; the second one is an optimization weighting item for the decision optimization variable; the third one is a terminal penalty term to ensure the stability of the system; the fourth one is the regularization term, which penalizes the auxiliary variable that simulates the outlier state and contains a scalar parameter that allows reducing noise pollution by adjusting the regularization quantity. The state error of the improved MPC converges to an asymptotically stable steady state, while ensuring the stability and unbiasedness of the system. Simulation cases show that this method can effectively reduce the impact of external abnormal noise on the system without polluting other data items.

The rest of this paper is organized as follows. Section 2 describes the structure of the ARMAX process and presents some fundamentals of MPC. Section 3 details how we develop the MPC for the ARMAX process with outlier modeling. Section 4 discusses and analyzes the simulation results. The conclusion is drawn in Section 5.

Preliminaries

The ARMAX process

Consider the following common ARMAX process:

$$A_i(q^{-1})z_k=B_i(q^{-1})u_k+C_i(q^{-1}){\omega }_k$$

(1)

$$y_k=z_k$$

(2)

where

$$\begin{array}{c}{A}_i(q^{-1})=1+a_1{q}^{-1}+\cdots +a_n{q}^{-n}\hfill \\ \multicolumn{1}{c}{B_i(q^{-1})=b_1{q}^{-1}+b_2{q}^{-2}+\cdots +b_n{q}^{-n}}\\ \multicolumn{1}{c}{C_i(q^{-1})=1+c_1{q}^{-1}+\cdots +c_n{q}^{-n}}\end{array}$$

are matrix polynomials in backward-shift operator $q^{-1}$ , $q^{-1}{z}_k=z_{k-1}$ , where $k$ is the step size, and $a_i\in {\mathbb{R}}^{m\times m}$ , $b_i\in {\mathbb{R}}^{m\times l}$ , $c_i\in {\mathbb{R}}^{m\times m}$ . $u_k\in {\mathbb{R}}^{n\times 1}$ is the input of process, $w_k\in {\mathbb{R}}^{l\times 1}$ is a bounded disturbance, $y_k$ is the actual observed output, and $z_k$ is the process output given by (2). In this paper, we focus on the abnormal noise added on the state, so here we set $y_k$ and $z_k$ to be equal.

The state space of the ARMAX process in (1) can be given as follows:

$$x_{k+1}=A{x}_k+B{u}_k+G{w}_k$$

(3)

$$z_k=H{x}_k+w_k$$

(4)

where

$$\begin{array}{cc}A\hfill & =\left[\begin{array}{ccccc}-a_1\hfill & I\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{-a_2}& 0\hfill & I\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{-a_{n-1}}& 0\hfill & 0\hfill & \cdots \hfill & I\hfill \\ \multicolumn{1}{c}{-a_n}& 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right],G=\left[\begin{array}{c}{c}_1-a_1\hfill \\ \multicolumn{1}{c}{c_2-a_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{c_n-a_n}\end{array}\right]\hfill \\ \multicolumn{1}{c}{B}& ={\left[\begin{array}{cccc}{b}_1\hfill & b_2\hfill & \cdots \hfill & b_n\hfill \end{array}\right]}^T,H=\left[\begin{array}{cccc}I\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right]\hfill \end{array}$$

MPC for the ARMAX process

To thoroughly analyze the problem, we first derive a standard MPC for the ARMAX process with general disturbances in this section. By expanding the derivation equation (3), the following $N+1$ equations can be obtained:

$$\begin{array}{l}{x}_{k+1}=A{x}_k+B{u}_k+G{w}_k\\ x_{k+2}=A^2{x}_k+AB{u}_k+B{u}_{k+1}+AG{w}_k+G{w}_{k+1}\\ ⋮\\ x_{k+N}=A^N{x}_k+\sum _{i=1}^N{A}^{i-1}B{u}_{k-i+1}\\ +\sum _{i=1}^N{A}^{i-1}G{w}_{k-i+1}\end{array}$$

(5)

where $u_k$ is the input at time $k$ of the ARMAX process. The input follows the feedback control rate and is determined by state feedback and process decision variables. It is given as follows:

$$u_{k+i}=K{x}_{k+i}+v_{k+i}$$

(6)

where $K$ is the linear system feedback gain or state feedback matrix. To facilitate the modification of the cost function later in this study, we rewrite the input quantity in the state prediction equation as the decision variable. For notational convenience, the vectors are defined as

$$\begin{array}{c}\multicolumn{1}{c}{X_k=[x_k^T,x_{k+1}^T,\dots ,x_{k+N}^T{]}^T}\\ \multicolumn{1}{c}{V_k=[v_k^T,v_{k+1}^T,\dots ,v_{k+N-1}^T{]}^T}\\ \multicolumn{1}{c}{W_k=[w_k^T,w_{k+1}^T,\dots ,w_{k+N}^T{]}^T}\end{array}$$

and the matrices are defined as

$$\begin{array}{c}{S}_x=\left[\begin{array}{c}(A+\mathrm{BK})A^T\hfill \\ \multicolumn{1}{c}{(A+\mathrm{BK})(A^2{)}^T}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{(A+\mathrm{BK})(A^N{)}^T{]}^T}\end{array}\right]\hfill \\ \multicolumn{1}{c}{S_v=\left[\begin{array}{ccccc}0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{B}& 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{\mathrm{AB}}& B\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{A^{N-1}B}& A^{N-2}B\hfill & A^{N-3}B\hfill & \cdots \hfill & 0\hfill \end{array}\right]}\\ \multicolumn{1}{c}{S_w=\left[\begin{array}{ccccc}0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{G}& 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{\mathrm{AG}}& G\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{A^{N-1}G}& A^{N-2}G\hfill & A^{N-3}G\hfill & \cdots \hfill & 0\hfill \end{array}\right]}\end{array}$$

where $X_k$ is the state trajectory predicted in the prediction horizon $N$ at time $k$ , $V_k$ is a sequence of decision variables, which is the objective solution of the optimization problem, and $W_k$ is the disturbance sequence acting on the system at time $k$ .

Equation (5) can be written in a more compact form as

$$X_k=S_x{x}_k+S_v{V}_k+S_w{W}_k$$

(7)

For the equation for predicting the state of the MPC system, it is given leveraging the system’s special structure and collected it in one place using a uniform notation. Next, we will improve the method described here to obtain a more general method. The goal of MPC is usually presented in the form of quadratic programing. At time $k$ , the objective function can be expressed as follows:

$$\begin{array}{c}\multicolumn{1}{c}{J=\sum _{i=0}^{N-1}[x_{k+i|k}^{T}Q{x}_{k+i|k}+u_{k+i|k}^{T}R{u}_{k+i|k}]}\\ \multicolumn{1}{c}{+x_{k+N|k}^T{P}_N{x}_{k+N|k}}\end{array}$$

(8)

where $x_{k+i|k}$ is a state vector that predicts $k+i$ at time $k$ and is extracted from the prediction model equation (7), $u_{k+i|k}$ is the input at time $k$ to predict time $k+1$ , $x_{k+N|k}$ is the terminal state, $P_N$ is the terminal penalty term given by solving the positive definite solution of the Lyapunov equation, $Q$ and $R$ are the weighted matrices for state quantities and control quantities, respectively, and $N$ is the prediction horizon.

Equation (6) is substituted and converted into a matrix equivalent cost function:

$$J=X_k^T\mathrm{Qxx}{X}_k+X_k^T\mathrm{Qxv}{V}_k+V_k^T\mathrm{Qvv}{V}_k$$

(9)

where

$$\begin{array}{c}{Q}_{\mathrm{xx}}=\left[\begin{array}{cccc}Q+K^T\mathrm{RK}\hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & Q+K^T\mathrm{RK}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & 0\hfill & P_N\hfill \end{array}\right]\hfill \\ \multicolumn{1}{c}{Q_{\mathrm{xv}}=\left[\begin{array}{ccccc}2K^{T}R\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 2K^{T}R\hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & 2K^{T}R\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & 0\hfill & 0\hfill \end{array}\right]}\\ \multicolumn{1}{c}{Q_{\mathrm{vv}}=\left[\begin{array}{cccc}R\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& R\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & R\hfill \end{array}\right]}\end{array}$$

Substituting equation (7) into equation (9), the deformation of the cost function $J$ can be given as follows:

$$J=V_k^T{\mathrm{\Phi }}_a{V}_k+x_k^T{\mathrm{\Theta }}_a{V}_k+x_k^T{\mathrm{\Gamma }}_a{x}_k$$

(10)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{\Phi }}_a=S_v^T(Q_{\mathrm{xx}}{S}_v+Q_{\mathrm{xv}})+Q_{\mathrm{vv}}}\\ \multicolumn{1}{c}{{\mathrm{\Theta }}_a=S_x^T(2Q_{\mathrm{xx}}{S}_v+Q_{\mathrm{xv}})}\\ \multicolumn{1}{c}{{\mathrm{\Gamma }}_a=S_x^T{Q}_{\mathrm{xx}}{S}_x}\end{array}$$

To find the minimum value, suppose that the gradient of $V_k$ with respect to the cost function $J_k$ is zero, and the analytic solution of $V_k$ can be obtained as follows:

$$V_k=(-0.5x_k^T\mathrm{\Theta }{\mathrm{\Phi }}^{-1}{)}^T$$

(11)

Here we find the solution of $V_k$ in the general MPC. In the following sections, we will compare the solution of the modified cost function with the target solution in the ordinary form to draw some conclusions.

Regularized MPC

In this part, we explicitly model outliers by introducing auxiliary vectors $O_k$ and further describe the proposed MPC. The objective of the controller is to stably produce the output of the MPC process, denoted by $y_k$ , taking into account abnormally contaminated state values, denoted by $x_k$ .

Outlier modeling

We define the auxiliary variable $O_k$ and introduce equation (3) to explicitly model the outlier noise, and propose an MPC for the ARMAX process with additive abnormal noise. Before proceeding further, we first define the vector containing time $k$ for the receding horizon of the latest abnormal noise batch as follows:

$$O_k=[o_k^T,o_{k+1}^T,\dots ,o_{k+N}^T{]}^T$$

(12)

The state space of the MPC is remodeled so that the state quantity is affected by outliers. The new process still retains the features of the old one. The new process model is given as follows:

$$\begin{array}{c}{x}_{k+1}=A{x}_k+B{u}_k+G{\omega }_k+o_k\hfill \\ \multicolumn{1}{c}{y_k=H{x}_k+{\omega }_k}\end{array}$$

(13)

The prediction equation for $X_k$ is rewritten according to the new equation of state. Equation (12) is added to equation (7) to obtain a new prediction equation of state as follows:

$$X_{o,k}=S_x{x}_k+S_v{V}_k+S_w{W}_k+O_k$$

(14)

It can be seen from the equation that data vector $X_k$ is affected by two terms related to abnormal noise $o_k$ and process noise ${\omega }_k$ . It is also affected by the historical state and noise, so the occurrence of a single error can affect multiple data items as well as the performance of controllers that do not explicitly consider abnormal noise.

• Remark 1. The difficulty with introducing $O_k$ is that calculating $u_k$ and $O_k$ based only on $x_k$ would be an undetermined problem (with fewer equations than unknowns), and it usually has an infinite number of solutions. Therefore, it is necessary to explicitly model the outliers to contribute a unique optimal solution to the problem. In practical applications, outlier sequences are usually related to large errors caused by unexpected events such as large perturbations and machine failures. We define outliers as continuous Gaussian noise in the time series, and there is a large number of outliers at a certain time point, which allows us to solve subsequent problems more easily.

The $l_2-\mathrm{norm}$ regularized MPC

The handling of outliers and data mutations is a major topic in the field of science. Outliers in measurements are usually considered as accidental large errors. ³⁰ In this specific case, $l_1$ -norm regularization can be used to restore sparsity and detect outliers. ³¹ In this paper, we consider a more general case of a zero-mean random variable with outliers. Considering that $l_1$ -norm regularization has low computational efficiency in the non-sparse case, we introduce $l_2$ -norm regularization, which has simple computationally efficient analytical solutions, into $O_k$ .

Based on the points discussed above, we introduce a regularization term into the cost function equation (9) of MPC. The regularization term consists of an optional parameter and the quadratic of $l_2$ -norm of the abnormal vector. The optimization problem after processing is equivalent to minimizing the following cost function:

$$\begin{array}{c}{J}_k=X_{o,k}^T\mathrm{Qxx}{X}_{o,k}+X_{o,k}^T\mathrm{Qxv}{V}_k+V_k^T\mathrm{Qvv}{V}_k\hfill \\ \multicolumn{1}{c}{+\beta {\left|\left|O_k\right|\right|}^2}\end{array}$$

(15)

where $\beta$ is a positive constant tuning parameter, which affects the degree of regularization, and $Q_{\mathrm{xx}}$ , $Q_{\mathrm{xv}}$ , and $Q_{\mathrm{vv}}$ are given by function (9). The objective problem can be described as the follows:

$$\begin{array}{c}[V_k,O_k]=\arg \min J_k\hfill \\ \multicolumn{1}{c}{s.t.x_{\min }\le x_k\le x_{\max }}\\ \multicolumn{1}{c}{u_{\min }\le u_k\le x_{\max }}\end{array}$$

where $x_{\min },x_{\max }\in \mathbb{R}$ and $u_{\min },u_{\max }\in \mathbb{R}$ . The disturbance $w$ is assumed to be bounded, which is represented as $w\in W$ , where $W$ is the compact set containing the origin.

Substituting equation (14) into equation (15) and based on $l_2$ -norm regularization, the optimization problem can be considered as solving the optimal solution of the following minimization function:

$$J_k=V_k^T\mathrm{\Phi }{V}_k+x_k^T\mathrm{\Theta }{V}_k+O_k^T\mathrm{\Gamma }{V}_k+O_k^T\mathrm{\Pi }{O}_k+x_k^T{T}_N{O}_k$$

(16)

where $\mathrm{\Phi }$ , $\mathrm{\Theta }$ , $\mathrm{\Gamma }$ , $\mathrm{\Pi }$ , $T_N$ are given by:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Phi }=S_v^T(Q_{\mathrm{xx}}{S}_v+Q_{\mathrm{xv}})+Q_{\mathrm{vv}}}\\ \multicolumn{1}{c}{\mathrm{\Theta }=S_x^T(2Q_{\mathrm{xx}}{S}_v+Q_{\mathrm{xv}})}\\ \multicolumn{1}{c}{\mathrm{\Gamma }=2Q_{\mathrm{xx}}{S}_v+Q_{\mathrm{xv}}}\\ \multicolumn{1}{c}{\mathrm{\Pi }=B_N+Q_{\mathrm{xx}}}\\ \multicolumn{1}{c}{T_N=2S_x^T{Q}_{\mathrm{xx}}}\\ \multicolumn{1}{c}{B_N=\mathrm{blkdiag}(\beta ,\dots ,\beta )}\end{array}$$

Gradient descent is an effective algorithm to solve the regularization of $l_2$ -norm. It works by step-by-step iterative work to obtain the minimum cost function and variable parameters of the model. The necessary condition for minimizing the cost function (16) using this algorithm is given by:

$$\begin{array}{c}{\nabla }_{V_k}{J}_k=0\hfill \\ \multicolumn{1}{c}{{\nabla }_{O_k}{J}_k=0}\end{array}$$

(17)

where

$$\begin{array}{c}{\nabla }_{V_k}{J}_k=2V_k^T\mathrm{\Phi }+x_k^T\mathrm{\Theta }+O_k^T\mathrm{\Gamma }\hfill \\ \multicolumn{1}{c}{{\nabla }_{O_k}{J}_k=V_k^T{\mathrm{\Gamma }}^T+2\mathrm{\Pi }{O}_k+x_k^T{T}_N}\end{array}$$

(18)

Combining equations (17) and (18), we can obtain

$$V_k=(-0.5(x_k^T\mathrm{\Theta }+O^T\mathrm{\Gamma }){\mathrm{\Phi }}^{-1}{)}^T$$

(19)

$$O_k=(-0.5(V_k^T{\mathrm{\Gamma }}^T+x_k^T{T}_N){\mathrm{\Pi }}^{-1}{)}^T$$

(20)

The optimization variable of the proposed MPC problem is the decision variable $V_k$ , so the ultimate goal of the algorithm is to obtain the value of $V_k$ when the function is optimal. Following this idea, we substitute equation (19) into equation (20) and apply the matrix inversion lemma to obtain

$$V_k=(-0.5(x_k^T\mathrm{\Theta }-0.5x_k^T{T}_N{\mathrm{\Pi }}^{-1}\mathrm{\Gamma }){\mathrm{\Phi }}^{-1}{E}^{-1}{)}^T$$

(21)

where

$$E=I-0.25{\mathrm{\Gamma }}^T{\mathrm{\Pi }}^{-1}\mathrm{\Gamma }{\mathrm{\Phi }}^{-1}$$

(22)

and $\mathrm{\Theta }$ , $T_N$ , $\mathrm{\Pi }$ and $\mathrm{\Gamma }$ are all given by function (16). Clearly, $E$ can be computed offline. In this method, the computational difficulty of the algorithm is reduced and its computational speed is improved.

• Remark 2. Regularized MPC has a similar form to normal MPC, but is actually more common because the former can be adjusted and selected to be the same form as the latter. Compared to equation (11), equation (21) has an additional term to $V_k$ , as shown in equation (22), where the effect of outlier noise is taken into account. Control becomes less reliable in the presence of outlier noise, so the term is used to adjust the relative trust in control behavior. If the control quantity is subject to a large number of uncertainties, the system needs to make a larger penalty for outliers, it is preferable to choose a smaller value for $\beta$ . Conversely, if anomalous noise is negligible and does not rely on the algorithm to reduce the output pollution rate, a larger $\beta$ is preferred. Actually, when $\beta \to \infty$ , $E\to 0$ and the regularized controller in equation (21) reduces to the ordinary controller in equation (11).

Examples

In this section, some examples are given to verify the effectiveness of the proposed MPC, and its performance is tested with simulation and experimental data.

Algorithm The $l_2$ -regularized MPC.
1. Give the initial values for $x_0$ , set predictive horizon $N$ and select a suitable value for $\beta$ .
2. Solve the Lyakonov equation, solve for $P_N$ , and calculate for $K$ .
3. Obtain the coefficient matrix $S_x$ , $S_v$ of the prediction equation, and calculate the weighted matrices $\mathrm{\Phi }$ , $\mathrm{\Theta }$ , $\mathrm{\Gamma }$ , $\mathrm{\Pi }$ , $T_N$ and $E$ offline.
4. Solve the optimization problem in (16) by (21) and get $V_k$ .
5. Calculate the system input $u_k$ and get the current state at the moment, so as to get the output value $y_k$ by (13).
6. Set $k=k+1$ and go back to step 4.

Example 1: Outlier

Consider the ARMAX process in (1) where

$$\begin{array}{c}\multicolumn{1}{c}{A_i(q^{-1})=I+a_1{q}^{-1}+a_2{q}^{-2}}\\ \multicolumn{1}{c}{B_i(q^{-1})=b_1{q}^{-1}+b_2{q}^{-2}}\\ \multicolumn{1}{c}{C_i(q^{-1})=I+c_1{q}^{-1}+c_2{q}^{-2}}\end{array}$$

(23)

where

$$\begin{array}{cc}{a}_1=-1,\hfill & a_2=-2\hfill \\ \multicolumn{1}{c}{b_1=0,}& b_2=1\hfill \\ \multicolumn{1}{c}{c_1=2,}& c_2=1\hfill \end{array}$$

and the measurement noise ${\omega }_k$ is $N(0,0.1)$ . The state and control inputs satisfy the constraints with $x_{\min }=[-3,-4]$ , $x_{\max }=[6,3]$ and $[u_{\min },u_{\max }]=[-5,5]$ . For the MPC, we choose a predictive horizon $N=6$ and a separable quadratic objective with $Q=I$ and $R=I$ . The outlier $O_k$ is modeled as the Gaussian noise in the polluted state. At $k=100$ , the outlier is given by

$$\begin{array}{c}{O}_k=\{\begin{array}{cc}-2,\hfill & k=100\hfill \\ \multicolumn{1}{c}{N(0,0.1),}& \mathrm{other}.\hfill \end{array}\hfill \end{array}$$

(24)

The process is again initialized with $x_0=[0,0]$ and $0.1$ is selected as the suitable variable $\beta$ . For simplicity, the reference signal is set to 0.

The selection of optional parameters $\beta$ determines the degree of regularization in the system. To evaluate how the selection of optional parameters $\beta$ on the system, different $\beta$ values are tested.

Figure 1 depicts the output of the proposed MPC when $\beta =4$ , 10 and 100, respectively. It can be seen that the amount of data at the polluted moment is affected by the value of beta when Gaussian noise has outliers. A smaller $\beta$ has a smaller penalty for outlier noise, but the output is greatly affected by the outlier. It is also clear that a smaller value of $\beta$ will make the regularization more conservative. The number of polluted items decreases when the value of $\beta$ increases, which greatly reduces the adverse effect of outliers on the system. When the value of $\beta$ continues to increase, the output curve will become much less obvious. This shows that the appropriate choice of $\beta$ can better handle outliers. Compared with the method of state estimation, the algorithm in this paper is more intuitive and reduces the complexity of the algorithm under the premise of achieving good performance.

Figure 1.

The output for the proposed MPC with different values of $\beta$ .

Example 2: Performance comparison between MPC and PID for ARMAX processes

Consider the ARMAX process described by equation (23) where the noise ${\omega }_k$ is generated in the range $(0,0.1)$ and the outlier noise $o_k$ is given by equation (24). The PID controller can be expressed in the following discrete form:

$$u_k=K_p{e}_k+K_i\sum _{j=1}^k{e}_j+K_d(e_k-e_{k-1})$$

(25)

where $e_k$ represents the difference between the reference signal $R_k$ and the output signal $y_k$ , which is given by

$$e_k=R_k-y_k$$

(26)

For simplicity, the reference signal is set to 0, and the PID controller parameters are set to $K_p=1$ , $K_i=0.5$ and $K_d=0$ , respectively. The simulation is again initialized with $x_0=[3,3]$ . Figures 2 and 3 compare the output and input signals of the PID controller and the regularized MPC, respectively. The results of the algorithm are simulated for 200 runs, the average of which is represented in a white curve. It can be seen from Figure 3 that the control deviation of the proposed MPC is smaller than that of the PID controller when the system is affected by outliers. Compared with the regularized MPC, the PID controller has a more obvious overshoot when returning to the reference trajectory. Hence, the proposed MPC delivers more precise control performance than the PID controller. This example shows that the proposed MPC can effectively prevent multiple data entries from being polluted by the outlier, and its overall performance is better than that of the PID controller.

Figure 2.

Output signals of PID and proposed MPC. The white curve represents the average of 200 runs: (a) PID and (b) regularized MPC.

Figure 3.

Input signals of PID and proposed MPC. The white curve represents the average of 200 runs: (a) PID and (b) regularized MPC.

Conclusions

In this paper, we present an MPC method for ARMAX processes with outlier noise. The outliers are modeled as auxiliary variables, and the cost function of the MPC is reconstructed by means of $l_2$ -norm regularization. The explicit modeling of outlier noise can help to prevent or mitigate the influence of outlier state on multiple instants of time, which leads to improved control accuracy. The regularization term is introduced into the cost function of MPC, which makes the cost function easier to parse due to the explicit modeling of noise. After deducing the analytical solution for the MPC, we recommend that different adjustment parameters should be chosen for different environments. Since the adjustment parameters are optional, the MPC is well suited for performing optimization behaviors in the presence of outlier noise. In addition, the calculation example shows that the proposed MPC can deal well with outliers to ensure excellent control performance, thereby freeing the control system from performance limitations and excessive errors.

The outlier noise introduced in this paper is universal to the mathematical model, and the proposed regularized MPC should be applicable to more scenarios. In future work, the applicability of this algorithm will be further explained and improved.