An interpretative and adaptive MPC for nonlinear systems

Most industrial approaches are based on linear models obtained from the linearization technique in solving MPC problems for nonlinear systems. A nonlinear plant generally admits a locally-linearized model when considering regulating a particular operating point. In most cases, a linear MPC yield adequate performance, especially for chemical process industries. This approach with one linear model has an advantage in designing and deploying the offline explicit MPC solutions based on multi-parametric quadratic programming [7] and the online MPC solutions based on fast convex optimization algorithms for embedded platforms to satisfy real-time requirements. For tracking problems that operate over a wide range of operating conditions, multiple linear model based MPC addresses the non-adequate accuracy issue in one linear model based MPC. Multiple linear models are also called linear parameter varying (LPV), often used in the context of gain-scheduling, namely an MPC controller scheduled by linear models [8]. One approach for obtaining an LPV model is the online successive linearization at the current states based on a first-principles model [9].

One improved approach is the linear time-varying (LTV) model, obtained from linearizing the nonlinear model at each prediction horizon [10]. However, the LTV-MPC approach comes at a higher computational cost than the LPV-MPC approach. Since the states and inputs during the prediction horizon are unknown, the typical way utilizes the shifted optimal inputs sequence of the previous MPC solution as the inputs, and the states are calculated by integrating the nonlinear model. Another successful and broadly used approach is the real-time iteration (RTI) scheme. The RTI approach assumes that the MPC solution at the current time is very similar to the solution obtained at the previous time. Under this assumption, a full Newton step can be taken, providing an excellent approximation of the fully converged NMPC solution. Not only that, but the RTI algorithm also divides the whole computation into preparation phase and feedback phase phases, to achieve a shorter feedback delay [11]. The RTI-NMPC approach has been implemented into the open-source software ACADO Toolkit, which allows exporting the optimized C-code for deployment [12, 13]. Another approximate scheme in the literature is the Continuation/GMRES method [14]. In [14, 15, 16], the authors also provide its full implementation within the tool AutoGenU to support the C-code generation.

In addition to the above approximating approaches, incorporating nonlinearity directly into the MPC problem can provide a systematic way of dealing with systems with nonlinear dynamics, constraints, and objectives but at the cost of heavy online computation. Solving the resulting non-convex optimization problem relies on an efficient and reliable nonlinear programming algorithm, which has been alleviated to some extent with the advent of tailored and professional software tools, such as the CasADi [17] and FORCES NLP [18].

To further reduce the online computation cost of NMPC, the tremendous advances in machine learning and deep learning inspired many works, such as learning a globally-linear model or learning an efficient representation of approximated MPC laws via deep neural networks. Herein we only name a few works related to MPC for nonlinear systems. In [19], the author utilized a lifting operator, called the Koopman operator, to lift the nonlinear dynamics into a higher dimensional space where its evolution is approximately linear. Then an MPC controller based on the lifted linear model can replace the NMPC. The lifted linear model can be learned by a user-defined dictionary library or deep neural networks from data [20, 21]. Another approach that directly learns the approximated MPC law to reduce the online computation cost has been explored in the literature. In [22], the authors show that a neural network can represent exactly the MPC law described by the piecewise affine function in the case of linear MPC. The authors further extended to the case of nonlinear systems [23], in which the deep neural network is utilized to learn the robust nonlinear MPC law.

In the practice of the MPC technique, in addition to numerical optimization algorithms needed to solve the MPC problem, another critical ingredient is estimation algorithms. Generally, the nonlinear dynamic model can be an interpretative first-principle state-space model or a black-box model built from experimental data via system identification (or modern machine learning) in an NMPC setting. A nonlinear first-principle state-space model can provide interpretability that is preferred in practice, but it often suffers the issues like unavailable full-state information, unmeasured disturbance, or unmodelled time-varying terms. Therefore, these cases require a state estimation algorithm or a joint state parameter estimation algorithm in which unmeasured disturbances and unmodeled time-varying terms are considered as parameters to be estimated. Many state estimation methods for nonlinear systems exist. The well-known extended Kalman Filter (EKF) is undoubtedly the predominant state estimation technique [24], and the EKF has also been applied in many joint state-parameter estimation applications [25, 26].

The identified nonlinear black-box model used in NMPC can be divided into two categories: state-space or input-output formulation, such as state-space based recurrent neural networks (RNNs) or input-output based neural-network autoregressive model with exogenous inputs (NNARX). One obvious advantage of input-output based models over state-space based models is that input-output based models don’t require an state estimation algorithm to estimate the generally higher-dimensional states. Despite this, input-output based models still require an online estimation algorithm to update the model parameters when discrepancies happen between the model output predictions and the streaming output measurements. The linear input-output ARX models are well known for their adaptability allowing efficient online recursive estimation algorithms, like recursive least-squares (RLS) or Kalman Filter [27]. In fact, nonlinear ARX models also have good adaptability, updating their model parameters online via the EKF algorithm [28].

As discussed above, the EKF algorithm is an essential component in NMPC practice, both for the first-principle based and black-box models. The EKF is based on the first-order linearization of nonlinear dynamics around the previously estimated vectors. Although the first-order linear approximation used by the EKF is sometimes not accurate enough, other advanced techniques can provide better estimation accuracy, but this comes at the cost of complicated implementation and expensive computational costs. The ease of implementation and computation burden are of great importance, especially in the cooperation of NMPC and EKF. Despite the widely practical usefulness of the EKF, its convergence and stability guarantee requires strict conditions, such as satisfying the nonlinear observability rank condition and having sufficiently small initial estimation error as well as disturbing noise term [29, 30]. In this paper we assume that the EKF algorithm will not suffer from divergence problems in its applications.

It is worth noting that the widely used SL-MPC approach and the EKF algorithm are both based on linearization techniques. Their combination, namely the SL-MPC with the EKF method presented in [31], has become a common choice to handle the nonlinear first-principle based model in the industrial practice of NMPC technique.In its EKF part, the states and unknown disturbances are recursively updated from the streaming input-output data based on the linearization of the discrete-time nonlinear model. In its SL-MPC part, the control input is calculated based on the estimated states and disturbances at the current time and the linearized state-space model, which is also updated online based on the linearization of a first-principle model at the estimated states and disturbances. Thus, the SL-MPC with the EKF method, in a sense, is utilizing the EKF algorithm and the streaming input-output data to continuously update online a linear state-space model, which is under the representation of states defined by the first-principle-based model. And a linear input-output plant could have infinite equivalent linear state-space realizations via coordinate transformation. Our previous work [32] also illustrated that a linear state-space model can be equivalently transformed into an input-output ARX model, which is well known for its adaptibility. In this paper, the SL-MPC with EKF method is the starting point to derive our interpretative and adaptive MPC (IA-MPC) method, which is illustrated in detail in Section II.

In our IA-MPC method, a linear state-space model is first obtained by performing the linearization of a first-principle based model at the initial point, then transformed into an equivalent ARX model; all of these are performed offline. This offline ARX model acquisition not only makes the black-box ARX model inherit the interpretability of the first-principle based model but also exploits its adaptive nature; this is why we call our method Interpretive and Adaptive MPC (IA-MPC). In the online closed-loop control, the EKF algorithm is used to update the ARX model parameters, and our previous construction-free CDAL-ARX algorithm [33] is used to compute the MPC control input.

The advantage of our IA-MPC method are sumarized as follows:

1. 1.

   In addition to the gained interpretability and adaptability, the acquisition of the offline ARX model also makes it possible to have a sufficiently small initial estimation error in ARX model parameters. All the EKF algorithm needs to do is the system tracking (or feedback correction), not the system identification.

2. 2.

   Compared to the EKF algorithm in the SL-MPC method, the EKF algorithm in our IA-MPC method can be decoupled when updating the ARX model parameters, thanks to the decoupling feature of ARX model. The decoupled EKF computation avoids matrix inversion and only involves scalar division, even allowing parallel calculation according to the number of outputs.

3. 3.

   The corresponding ARX-based MPC problem is solved by our previous construction-free, matrix-free, and library-free CDAL-ARX algorithm. Thus, our IA-MPC method would significantly reduce the difficulty in deploying nonlinear MPC on embedded platforms.

The extended Kalman Filter exploits the idea that performs the linearization at each sampling time to approximate the nonlinear system as a time-varying system and apply the linear filtering theory. In digital controller design, $u$ and $d$ are assumed to a constant value between the sampling time. Thus, a discrete-time model of (1) can be formulated as follows:

where $F_{t_{s}}(x_{k-1},u_{k-1},d_{k-1})$ represents the the terminal state vector obtained by integrating the ordinary differential equation (1a) for one sample interval $t_{s}$ with the initial condition of $x_{k-1}$ and constant inputs of $u_{k-1}$ and $d_{k-1}$. For a better state estimation, the unmeasured disturbance should also be estimated simultaneously. In general, $d$ can be modeled as the following stochastic difference equation

where $w_{k-1}$ is white noise sequence. Combining the Eqn (2) with (3), the following augmented model is obtained

where the measurement noise sequence $v_{k}$ is often to appear in the measured output $y_{k}$. By linearizing (4) at $\{\hat{x}_{k-1},\hat{d}_{k-1}\}$, the EKF computes the new estimates $[\hat{x}_{k}^{\prime},\hat{d}_{k}^{\prime}]^{\prime}$ from the feedback measurement $y_{k}$ and the model prediction (4e)

where

$\Lambda_{k-1}$,$\Lambda_{k-1}^{d}$, $H_{k-1}$ and $H_{k-1}^{d}$ are calculated using the following formula

where $\partial F_{t_{s}}/\partial x$ and $\partial F_{t_{s}}/\partial d$ represent Jacobian matrices of $F_{t_{s}}$ with respect to $x$ and $d$, respectively. And $\partial g/\partial x$ and $\partial g/\partial d$ represent Jacobian matrices of $g$ with respect to $x$ and $d$, respectively. Thus, we can conclude the EKF computation procedure which has been well presented and analysed, see Tab. I.

By recursively handling the streaming input-output data $(u_{k-1},y_{k-1})$, the EKF algorithm could calculate the new estimated states and disturbances $(\hat{x}_{k},\hat{d}_{k})$. The $(\hat{x}_{k},\hat{d}_{k})$ is not only used as the nominal values to linearize the nonlinear first-principle model but also used as the initial conditions of the linearized state-space model in the SL-MPC method. Note that the continuous-time model (1) has to be transformed into the discrete-time model for digital MPC design. Two approaches, first discretise then linearise and first linearise then discretise, have been reported. The EKF algorithm presented in Section II-A is based on the first discretise then linearise approach. In the context of SL-MPC, it is more common to use the first linearise then discretise. Here the continuous-time first-principle-based model (1) is linearized at current points $\{\hat{x}_{k},u_{k-1},\hat{d}_{k}\}$,

where $A_{c}=\frac{\partial f}{\partial x}|_{\hat{x}_{k},u_{k-1},\hat{d}_{k}}$, $B_{c}=\frac{\partial f}{\partial u}|_{\hat{x}_{k},u_{k-1},\hat{d}_{k}}$, $e_{c}=f(\hat{x}_{k},u_{k-1},\hat{d}_{k})$, $C=\frac{\partial g}{\partial x}|_{\hat{x}_{k},\hat{d}_{k}}$ and $h_{c}=g(\hat{x}_{k},\hat{d}_{k})$. Here the differential eqn (8a) has to be discretised for obtaining a discrete-time model. The discretization method includes the exact and the approximated discretization. One common discretization method is the the one-step Euler’s approximate method, then a discrete-time model is obtained as follows

where $e=t_{s}\left(e_{c}-A_{c}\hat{x}_{k}-B_{c}u_{k-1}\right)$, $A=I+t_{s}A_{c}$, $B=t_{s}B_{c}$, $h=h_{c}-C\hat{x}_{k}$. Then, a MPC tracking problem formulation is listed as follows

The MPC tracking problem (II-B) can be solved by many quadratic programming (QP) algorithms, such as the active-set, the interior-point and the ADMM-based solver. Note that the model parameters $\{A,B,e,C,h\}$ in (9) are related to the recursively updated estimates $\{\hat{x}_{k},\hat{d}_{k}\}$ and the last control input $u_{k-1}$, it means that the MPC tracking problem (II-B) has to be reformulated at each sampling time. Most QP solvers requires an explicit condensing or sparse MPC-to-QP construction. In such a SL-MPC scenario, the MPC-to-QP construction has to be performed online at each sampling time, as does solving the QP problem. Indeed, often the online MPC-to-QP construction has a comparable computational cost to solving the QP problem itself online, especially when warm-starting strategies are employed. Our previous work [34], a construction-free CDAL solver, can avoid explicit MPC-to-QP construction in state-space MPC problem.

In our preivous work [32], we conclude three SS-to-ARX transformations, namely the Cayley-Hamilton-based, Observer-Theory-based and Kalman-Filter-based transformations. The Cayley-Hamilton-based tranformation can derive the unique and equivalent ARX model but with sensitivity to noises. The Observer-Theory-based and Kalman-Filter-based transformations has the same structure in the transformed ARX model, and the Kalman-Filter-based transformation requires a larger order than the Observer-Theory-based transformation. In fact, the order determines the noises robustness of the ARX model, which can easily tuned by the pole placement method. Clearly, larger order brings better noises robustness but also leads to more computational cost in ARX model update and ARX-based MPC. Herein we mainly present how to transform the adaptive linearized state-space model (9) with the help of the the Observer-Theory-based transformation.

By utilizing a gain matrix $L$ in the SS model (9) to begin the derivation of the SS-to-ARX transformation as follows

Based on the above evolution equation (11), calculating from previous time $t-p$ to current time $t$ can lead to the following equation

If a gain matrix $L$ exists such that $(A-LC)^{p}$ vanishes to zero, i.e.,

From linear system theory, the existence of such an gain matrix $L$ is assured as long as the system is observable. Thus, the term $(A-LC)^{p}x_{t-p}$ in (12) is zero for $k\geq p$, then multiplying the matrix $C$ on both sides of the equation (12) can derive the following ARX model,

where $\Psi_{i}=C(A-LC)^{i-1}L$, $\Omega_{i}=C(A-LC)^{i-1}B$ and $\zeta=\sum_{i=1}^{k}C(A-LC)^{i-1}e$

Next, we provide the analysis to tell why the gain matrix $L$ can be viewed as an observer gain. The original state-space model (9) has an observer gain $L$ of the following form

where $\hat{x}_{k}$ is the estimated state. The state estimation error can be denoted as $\epsilon_{k}=x_{k}-\hat{x}_{k}$, then $\epsilon_{k}$ follows the dynamic equation,

The estimated state $\hat{x}_{k}$ will converge the actual value $x_{k}$ as $k$ tends to infinity if the matrix $A-LC$ is asymptotically stable, namely the condition (13), which can be satisfied by using the pole placement method.

By viewing the time-varying ARX model parameters as system states, the ARX model (14) is re-formulateed as the following state-space form

where $\{w_{k}\}$ is a sequence of independent random vectors and $\{v_{k}\}$ is a output noise sequence. Although the variances of $\{w_{k}\}$ and $\{v_{k}\}$ are unknown in real applications and the actual parameters $\{\theta_{k}\}$ may change differently from the above random walk model (17a), the EKF algorithm can still work very well [36].

Let’s write the Eqn (17b) as the following separated formulation

where $n_{y}$ denotes the dimensions of $y$, $y_{k}(j)$ denotes the $j$-th element of $y_{k}$. $\theta_{k}$ and $\varphi_{k}$ have the following relationship according to (14)

The advantage of the above separated ARX formulation is that it allows the EKF algorithm to run in a parallel way and avoids the matrix-inverse operation, only involving scalar division, compared to the EKF in the SL-MPC (see Tab. I). We list the EKF computation scheme of our IA-MPC framework in Tab. III.

In our IA-MPC framework, the ARX model is updated by the EKF algorithm at each sampling time. Thus it leads to the corresponding time-changing ARX-based MPC tracking problem as follows.,

In such a a time-changing setting, the computation time spent in constructing its optimization problem and solving itself should be together considered. Indeed, often the online MPC-to-QP construction has a comparable computational cost to solving the MPC problem itself online, especially when warm-starting strategies are employed. Our previous construction-free CDAL-ARX algorithm [33] can avoid the explicit MPC-to-QP construction, thus allowing it to be suitable for this scenario.

In each iteration of our CDAL-ARX algorithm, an augmented Lagrangian (AL) subproblem is solved by the coordinate descent (CD) method, and the accelerated Nesterov’s scheme is used to update the Lagrangian dual variables. In particular, an efficient coupling scheme between CD and AL method is proposed to exploit the structure of the ARX-MPC problem, which can reduce the computation cost of each inner iteration. The detailed description of algorithm implementation can be found in [33]. In addition to the notable construction-free feature, our CDAL-ARX is also matrix-free and library-free which makes it practically useful in embedded deployment.

1. 1.

   Two tank problem: we consider the cascaded two tanks system, which is a fluid level control system. The input signal controls the water pump that pumps the water from a reservoir into the upper water tank. The water of the upper water tank also flows through a small opening into the lower water tank. The water of the lower water tank flows through a small opening into the reservoir. Herein without considering the overflow effect, we adopt the Bernoulli’s principle and conservation of mass to derive the following first-principle model:

   	$\displaystyle\dot{x}_{1}$	$\displaystyle=-k_{1}\sqrt{x_{1}}+k_{2}u$		(21)
   	$\displaystyle\dot{x}_{2}$	$\displaystyle=k_{1}\sqrt{x_{1}}-k_{3}\sqrt{x_{2}}$
   	$\displaystyle y$	$\displaystyle=x_{2}$

   where $x_{1}$ and $x_{2}$ are the water level of the upper and lower water tank, respectively. The full states cannot be measured only the $x_{2}$ as the measured output $y$. $u$ is the input signal, and $k_{1}$ , $k_{2}$ and $k_{3}$ are constants depending on the system properties, herein we adopt value $0.5$ for all of them. The sampling time of discrete digital control is $0.2$ s and the control goal is to make the output $y$ track the given reference signal subject to the input constraints $0\leq u\leq 2$ and the input increment constraints $-0.5\leq\Delta u\leq 0.5$.

2. 2.

   Bilinear motor problem: one common nonlinear control benchmark is a bilinear DC motor plant, whose equation is described as follows,

   	$\displaystyle\dot{x}_{1}$	$\displaystyle=-\left(R_{a}/L_{a}\right)x_{1}-\left(k_{m}/L_{a}\right)x_{2}u+u_{a}/L_{a}$		(22)
   	$\displaystyle\dot{x}_{2}$	$\displaystyle=-(B/J)x_{2}+\left(k_{m}/J\right)x_{1}u-\tau_{l}/J$
   	$\displaystyle y$	$\displaystyle=x_{2}$

   where $x_{1}$ and $x_{2}$ are the rotor current and angular velocity, respectively. The full states cannot be measured only the $x_{2}$ as the measured output $y$. The control input $u$ is the stator current. The system parameters are $L_{a}=0.314$, $R_{a}=12.345$, $k_{m}=0.253$, $J=0.00441$, $B=0.00732$, $\tau_{l}=1.47$, $u_{a}=60$. The bilinearity of (22) appears between the state and the control input. The sampling time of discrete digital control is $0.01$ s and the control goal is to make the output $y$ track the given reference signal subject to the input constraints $0\leq u\leq 2$ and the input increment constraints $-1\leq\Delta u\leq 1$.

3. 3.

   CSTR problem: one typical nonlinear process control benchmark is a continuous stirred tank reactor (CSTR) problem, which is described by the following continuous-time nonlinear model,

   	$\displaystyle\dot{C_{A}}$	$\displaystyle=C_{A,i}-C_{A}-k_{0}e^{\frac{-EaR}{T}}C_{A}$		(23)
   	$\displaystyle\dot{T}$	$\displaystyle=T_{i}+0.3T_{c}-1.3T+11.92k_{0}e^{\frac{-EaR}{T}}C_{A}$
   	$\displaystyle y$	$\displaystyle=T$

   where $C_{A}$ is the concentration of reagent A, $T$ is the temperature of the reactor, $C_{A,i}$ is the inlet feed stream concentration, which is assumed to have the constant value $10.0$ kgmol/m³. The unmeasured disturbance comes from the inlet feed stream temperature $T_{i}$, which has slow fluctuations represented by $T_{i}=298.15+5\sin(0.05t)$ $K$. The manipulated variable is the coolant temperature $T_{c}$. The constants $k_{0}=34930800$ and $EaR=-5963.6$ (in MKS units). The sampling time of discrete digital control is $0.5$ s and the control goal is to manipulates the coolant temperature $T_{c}$ to track a higher temperature of the reactor (equals a higher conversion rate) as well as reject the unmeasured disturbance $T_{i}$. The physical constraints comes from the input increment constraints $-1\leq\Delta T_{c}\leq 1$. Note that the unmeasured disturbance $T_{i}$ is estimated by the EKF algorithm in SL-MPC and NMPC methods.

4. 4.

   Van der Pol oscillator problem: we consider the nonlinear Van der Pol oscillator with a time varying parameter, and its dynamics are given by,

   	$\displaystyle\dot{x}$	$\displaystyle=v$		(24)
   	$\displaystyle\dot{v}$	$\displaystyle=\mu(t)(1-x^{2})v-x+u$

   where $x$ is the position, $v$ is the velocity, $u$ is the control input and $\mu(t)$ is a piecewise function defined as

   $$\mu(t)=\begin{cases}1&\text{ if }t\leq 50\\ 3&\text{ if }t>50\end{cases}$$

   The sampling time of discrete digital control is $0.2$ s and the control goal is to make the output $y$ track the given reference signal subject to the input constraints $-10\leq u\leq 10$ and the input increment constraints $-10\leq\Delta u\leq 10$. Note that the unmeasured time-varying parameter $\mu(t)$ is estimated by the EKF algorithm in SL-MPC and NMPC methods.

The above four problems share the same MPC and EKF settings: the MPC prediction horizon $N_{p}=10$, the MPC output weight $W_{y}=10$ and the MPC input increment weight $W_{\Delta u}=0.1$, the initial error covariance matrix of EKF $P_{0}=10I$, the process noise covariance matrix of EKF $Q=0.01I$ and the measurement noise covariance of EKF $R=0.01$. We consider two simulation scenarios, that is without and with process noises, respectively.

1. 1.

   Two tank problem: its process noise scenario is given a random noise $0.05\times$rand(2,1), and the initial conditions are $x_{1}=1$, $x_{2}=1$ and $u=1$ for obtaining the linearized state-space model, which is then transformed into the ARX model with order $p=3$ by using poles $[0.01,0.02]$. Its tracking signal is randomly selected every $20$ s in the range $[1,3]$.

2. 2.

   Bilinear motor problem: its process noise scenario is given a random noise $1\times$rand(2,1), and the initial conditions are $x_{1}=5.2542$, $x_{2}=-19.2205$ and $u=1$ for obtaining the linearized state-space model, which is then transformed into the ARX model with order $p=5$ by using poles $[0.05,0.1]$. Its tracking signal is randomly selected every $0.4$ s in the range $[-10,10]$.

3. 3.

   CSTR problem: its process noise scenario is given a random noise $0.1\times$rand(2,1), and the initial state of the reactor is at a low conversion rate, with $C_{A}=8.57$ kgmol/m³, $T=311$ K, and then the linearized state-space model around the initial conditions is transformed into the ARX model with order $p=3$ by using poles $[0.01,0.02]$. Its tracking signal gradually changes from $311.2639$ K to $370$ K during $50$ $100$ s and then holds constant.

4. 4.

   Van der Pol oscillator problem: its process noise scenario is given a random noise $1\times$rand(2,1), and the initial conditions are $x_{1}=0$, $x_{2}=0$ and $u=0$ for obtaining the linearized state-space model, which is then transformed into the ARX model with order $p=3$ by using poles $[0.005,0.01]$. Its tracking signal switches between $0$ and $1$ every $10$ s.

Their simulation results are plotted in Fig. 3, 4, 5 and 6, respectively. The Two tank and Bilinear motor problem have no unknown disturbances and the NMPC, SL-MPC and our IA-MPC method generate the same offset-free tracking performances subject to the input constraints in their noise-free scenario, which are plotted in Fig. 3(a) and 4(a). The CSTR and Van der Pol oscillator problems have the unknown time-varying disturbance, it’s found that our IA-MPC method generate the better offset-free tracking performances than the NMPC and SL-MPC method under the given EKF setting in their noise-free scenario, which are plotted in Fig. 5(a) and 6(a). Under the process noises, our IA-MPC method generate the better noise robustness than the NMPC and SL-MPC method in the four problems, which are plotted in Fig. 3(b), 4(b), 5(b) and 6(b). As expected, the online computation time of the exact NMPC method is much longer than the approximate SL-MPC and IA-MPC method, and the online computation time of the SL-MPC and IA-MPC method are almost the same since they both utilized our developed construction-free CDAL method, having no online construction time.