Computationally Efficient Nonlinear Model Predictive Control Using the L1 Cost-Function

Model Predictive Control (MPC) algorithms typically use the classical L2 cost function, which minimises squared differences of predicted control errors. Such an approach has good numerical properties, but the L1 norm that measures absolute values of the control errors gives better control quality. If a nonlinear model is used for prediction, the L1 norm leads to a difficult, nonlinear, possibly non-differentiable cost function. A computationally efficient alternative is discussed in this work. The solution used consists of two concepts: (a) a neural approximator is used in place of the non-differentiable absolute value function; (b) an advanced trajectory linearisation is performed on-line. As a result, an easy-to-solve quadratic optimisation task is obtained in place of the nonlinear one. Advantages of the presented solution are discussed for a simulated neutralisation benchmark. It is shown that the obtained trajectories are very similar, practically the same, as those possible in the reference scheme with nonlinear optimisation. Furthermore, the L1 norm even gives better performance than the classical L2 one in terms of the classical control performance indicator that measures squared control errors.


Introduction
In Model Predictive Control (MPC), a dynamical model of the process is used online to repeatedly make predictions of the future values of the controlled variables and to optimise the current and future control policy [1,2]. Due to such a formulation, very good closed-loop control accuracy is obtained and all necessary constraints may be easily imposed on process variables. MPC algorithms are utilised for industrial process control, e.g., chemical reactors [3] and distillation columns [4]. In addition to that, MPC algorithms are also used for fast dynamical systems, e.g., electromagnetic mills [5], electromechanical systems [6], servomotors [7], quadrotors [8], autonomous vehicles [9][10][11], unmanned aerial vehicles [12] and stochastic systems [13]. Four conditions must be fulfilled to obtain good control quality: good measurements, a precise model, an adequate choice of the MPC cost function and a fast MPC computation scheme carried out on-line.
MPC algorithms depend heavily on precise measurements of process variables that are provided by sensors. In some versions of automatic control systems, e.g., those described in [12,14,15], it is stressed that all necessary variables must be measured because otherwise, a significant loss in control performance is unavoidable. If the measurements are not available, the typical approach is to perform on-line estimation using Kalman or Extended Kalman filters [16]. Furthermore, there are other alternative approaches to solve this problem to some extent. They are usually developed for particular applications. For example, when the measurement of lateral vehicle speed is not available, it can be estimated from other measurable parameters [10]. The study presented in [17] shows another approach of dealing with faulty sensors for linear systems by changing the representation of the process model. The work in [18] presents a real vehicle that uses an external camera to detect obstacles and lanes on the road as well as external rear-corner radars to detect objects coming from the rear. An interesting example is presented in [19], where an anemometer is used to measure external disturbances such as wind force and direction. In [20], a depth sensor is installed for sea ship depth measurement. Moreover, the heave speed is obtained by calculating the derivative from depth sensor data. Finally, an example outside the vehicle field where fault-tolerant control is handled using MPC is stiction in control valves [21].
Typically, the minimised cost function used in MPC measures the sum of squared control errors predicted over some time horizon (the MPC-L 2 algorithm). Such a formulation is computationally simple and has good numerical properties. The use of the sum of absolute values of the predicted control errors (the MPC-L 1 algorithm) is much less popular. However, it may be easily verified that the second approach leads to better control quality [22][23][24][25]; the result is independent of the indicator used [23]. Furthermore, the classical L 2 cost function is not always suitable from the point of stability, as pointed out in [26]. It must be also stressed that in regression (model identification), the classical L 2 approach yields solutions fragile to outliers, while the L 1 norm leads to robustness [27].
It is important to emphasise that the numerical difficulty associated with the MPC-L 1 optimisation task depends on the type of dynamical model used for prediction. Provided that a linear model is used, one obtains a linear optimisation task [2,22], which can be solved very efficiently using the classical simplex method. Unfortunately, it is not possible when a nonlinear model is used. In such a case, a nonlinear optimisation problem with constraints must be solved at each sampling instant on-line, which may be computationally too demanding or completely impossible. Examples of nonlinear MPC-L 1 algorithms in which the Sequential Quadratic Programming algorithm is used are described in [23,24]. A trust-region sequential quadratic programming method is used in [28].
This work describes computationally efficient approaches to the MPC-L 1 algorithms for nonlinear processes. In spite of the fact that a nonlinear model is used for prediction, a computationally simple quadratic optimisation task is solved at each sampling instant on-line, and nonlinear optimisation is not used. The proposed solution is based on two concepts: a neural approximator and advanced on-line linearisation. A practical approach to nonlinear MPC is to perform model or trajectory linearisation on-line. As a result, an easyto-solve quadratic optimisation problem is obtained. Details of a few such algorithms may be found in [29]. An alternative is to use the fuzzy approach, in which a combination of multiple linear models approximates the nonlinear process, which also results in quadratic optimisation, e.g., [30,31]. Unfortunately, in research carried out so far, only the classical MPC-L 2 cost function has been used. The main difficulty is the fact that the cost function in the rudimentary MPC-L 1 optimisation task is not differentiable. To make it possible, in the presented approach, the classical MPC-L 1 cost function is replaced by its differentiable representation. Such a representation is obtained by means of the neural approximator. Neural networks of the Multi Layer Perceptron (MLP) structure with two layers are used for two reasons. Firstly, they are universal approximators, capable of approximating nonlinear relations with great accuracy. Secondly, the obtained neural approximation is differentiable and poses no numerical problems for the presented computational scheme. The efficiency of the presented approaches is demonstrated for a simulated neutralisation process. In particular, the discussed MPC-L 1 algorithms are compared with the classical MPC-L 2 ones.
It is worth stressing that different neural structures have been used in MPC for different purposes. First of all, they may be used for prediction: (a) Neural networks are most often used as black-box models of dynamical processes.
Various structures are used: the classical MLP networks [29,32,33], Radial Basis Function (RBF) networks [34][35][36], Long Short-Term Memory (LSTM) [37][38][39] and Gated Recurrent Unit (GRU) [39] structures. Typically, the input-output neural models are used. The state-space neural models [29,40,41] are used when the state-space process description is necessary, although such an approach is significantly less popular. (b) Neural networks are also used in block-oriented cascade dynamical models, which consist of neural static blocks and linear dynamic blocks connected in series. Hammerstein [42][43][44] and Wiener [45][46][47] models are the most frequently used cascade models. (c) Quasi-linear neural models [48,49]. In this approach, the dynamical model has the classical linear form, but its parameters depend on the operating point of the process and their values are determined on-line by neural networks. (d) Neural step response models [50]. In this approach, time-varying coefficients of the model are computed on-line by a neural network. (e) Neural multi-models [51,52]. In this approach, separate networks calculate the predictions for the consecutive sampling instants over the prediction horizon. As a result, the neural model is not used recurrently, which significantly simplifies training. Additionally, prediction errors are not propagated. (f) Hybrid neural models [53]. In this approach, neural networks are used to calculate the parameters of the first-principle models. (g) Neural networks may be used for modelling and MPC of distributed parameter systems [54][55][56].
Additionally, neural networks may be utilised to accelerate and simplify on-line calculations in MPC: (a) Neural inverse static models are used to try to cancel process nonlinearity. In particular, such a method is frequently used when Wiener cascade models are considered [47,57]. As a result, a quadratic optimisation task is obtained in place of a nonlinear one. (b) A neural approximator may be used to find the initial solution of the MPC optimisation problem, which speeds up calculations [58,59]. (c) Neural networks are able to approximate the MPC control law [60][61][62]. For training, sufficiently rich data sets are necessary, obtained for different operating points. (d) Specialised recurrent neural networks may be used to solve the MPC optimisation task on-line [63,64]. As a result, numerical optimisation is not necessary.
The article is organised in the following way. Section 2 recalls the general MPC-L 1 and MPC-L 2 optimisation tasks. The main part of the article, given in Section 3, details the computationally efficient nonlinear MPC schemes in which the L 1 cost function is used. Section 4 thoroughly discusses simulation results for a neutralisation reactor. In particular, MPC-L 1 and MPC-L 2 algorithms are compared. Finally, Section 5 summarises the whole article.

Problem Formulation
In this work, MPC of a Single-Input Single-Output processes is considered. The input of the process, i.e., the manipulated variable, is denoted by u. The output of the process, i.e., the controlled variable, is denoted by y. The vector of decision variables calculated on-line at each sampling instant of MPC is defined as a set of N u increments of the manipulated variable: . . .
where N u is named the control horizon. The decision variables (1) are computed from a constrained MPC optimisation task. The general form of the rudimentary MPC optimisation task considered in this work is: In this work, two types of constraints are considered: the magnitude constraints imposed on the manipulated variable, defined by u min and u max , and the constraints imposed on the increments of that variable, defined by u min and u max . These constraints are considered over the control horizon, N u . Having calculated the decision vector (1) from the MPC optimisation problem (2), the first element of the obtained sequence is applied to the process and the whole computational scheme is repeated at the consecutive sampling instant. The classical minimised cost function used in MPC (the MPC-L 2 algorithm) has the form: The first part of the cost function measures the sum of squared control errors predicted over the prediction horizon N; the symbols y sp (k + p|k) andŷ(k + p|k) denote the set-point and the predicted value of the controlled variable, respectively, both for the future sampling instant k + p known (it refers to the set-point) or calculated (it refers to the predictions) at the current sampling instant k. The predictions,ŷ(k + p|k), are calculated using a dynamical model of the controlled process. The second part of the cost function is used to minimise excessive changes of the manipulated variable; λ is a weighting coefficient. Additionally, it provides good numerical properties.
In the the MPC-L 1 algorithm, the sum of absolute values of the predicted control errors over the prediction horizon are taken into account rather than the sum of squared errors. Hence, the minimised cost function is: To penalise significant changes of the manipulated variable and to obtain good numerical properties, the second part of the cost function is the same as in the MPC-L 2 formulation [23].

Computationally Efficient Nonlinear MPC Using the L 1 Cost-Function
Let us formulate the MPC-L 1 optimisation task. Taking into account the general MPC problem (2) and the MPC-L 1 cost function (4), we obtain: The use of a nonlinear model for prediction leads to two computational difficulties. Firstly, predictionsŷ(k + p|k) are nonlinear functions of the calculated decision vector (1), which means that the cost function is nonlinear. As a result, we obtain a nonlinear optimisation task that must be repeated at each sampling instant. Secondly, the absolute value function is not differentiable, which means that the classical gradient-based optimisation method cannot be used. In spite of the second of the mentioned computational difficulties, in the literature, is it possible to find applications of gradient-based nonlinear optimisation methods to solve the MPC-L 1 optimisation task (5) [23,24,28]. The objective of this work is to derive a much more computationally simple approach to the MPC-L 1 problem in which nonlinear optimisation is not used. To achieve this result, two concepts are used: (a) The first part of the non-differentiable cost function (4) is replaced by its differentiable representation. For this purpose, a neural network approximation of the absolute value function is used. (b) The J 1 (k) cost function with a neural approximator is differentiable but nonlinear in terms of the computed control moves (1). To simplify the calculation scheme, an advanced trajectory linearisation method is used. As a result, a simple-to-solve quadratic optimisation task is obtained in place of the nonlinear one. Quadratic optimisation problems, for λ > 0, have only one minimum, which is the global one.

Neural Approximation of the MPC-L 1 Cost-Function
Only the first part of the cost function J 1 (k) defined by Equation (4) is not differentiable. Hence, a differentiable approximator of the first part is only necessary. At first, let us define the predicted control error for the future sampling instant k + p computed at the current instant k: where p = 1, . . . , N. Hence, the cost function J 1 (k) may be rewritten compactly: In order to use quadratic optimisation in MPC, it is postulated to use the general form of the differentiable approximation of the cost function (7): where α(e(k + p|k)) describes some nonlinear function of the predicted control error. Comparing the cost functions (7) and (8), it is clear that the following approximation of the absolute value function is required: for all p = 1, . . . , N. In order to obtain a differentiable approximation of the absolute value function, the neural network of a Multi Layer Perceptron (MLP) type is used [65]. The network has two layers, the first of which (the hidden one) is nonlinear, and a linear output. It is defined by the following equation: where K denotes the number of hidden neurons and ϕ is a nonlinear activation function. The weights of the first layer of the network are denoted by w 1 i,0 and w 1 i,1 for i = 1, . . . , K. The weights of the second layer are denoted by w 2 0 and w 2 i for i = 1, . . . , K. Provided that a differentiable activation function ϕ is used, the function α(e(k + p|k)) is also differentiable. It is of course true for the tanh function, which is used in this work.

Advanced Trajectory Linearisation of the MPC-L 1 Cost-Function
Thanks to the use of neural approximation (10) of the absolute value function, the cost function J 1 (k) defined by Equation (8) is differentiable, but still nonlinear. Additionally, the model of the controlled process is nonlinear. As a result, the cost function J 1 (k) is nonlinear in terms of the calculated future control increments (1). To simplify the problem, we use the advanced on-line trajectory linearisation approach. Linearisation is not carried out in a simple way, for the current or past operating point of the process, but for some assumed future trajectory of the manipulated variable, defined over the control horizon: . . .
Using the Taylor series expansion formula, the linear approximation of the multivariable function α(e(k + p|k)) that has N u arguments, defined by the vector u traj (k), is: for p = 1, . . . , N. In order to find the partial derivatives in the right side of Equation (12), the neural approximator defined by Equation (10) is differentiated with respect to the assumed future trajectory of the manipulated variable (11), which gives: for all p = 1, . . . , N, r = 0, . . . , N u − 1. The inputs of the hidden nodes of the neural network are: where i = 1, . . . , K, p = 1, . . . , N. If the tanh function is used in the hidden layer of the neural network: Finally, from Equations (13)-(15), we obtain: for all p = 1, . . . , N, r = 0, . . . , N u − 1. The partial derivatives in the right side of Equation (16) are derived for a particular model structure used for prediction. Differentiability of the model is required. Let us stress that in Equation (12), independent linear approximations are obtained for the consecutive sampling instants over the prediction horizon, i.e., for p = 1, . . . , N. In MPC, we need approximations of the absolute value of the predicted error over the whole prediction horizon. In order to simplify derivations, a compact vector matrix notation is used. The predicted trajectory of the function α, i.e., the vector form of Equation (12), is the following: where vectors of length N have the forms: and the vector of length N u , corresponding to the vector of increments (1), is: 20) and the N × N u matrix of partial derivatives has the structure: The matrix (21), where its entries are defined by Equation (16), is calculated for the specific neural approximator and the nonlinear model used. Because in MPC we calculate not the future values of the manipulated variables (Equation (20)) but the corresponding increments (1), we have to express the vector equation of the linearised trajectory (Equation (17)) in terms of the vector of control increments u(k). We obtain: where the auxiliary matrix of dimensionality N u × N u is: and the vector of length N u has the structure:

Formulation of the Computationally Simple MPC-L 1 Quadratic Optimisation Task
In order to obtain a quadratic optimisation MPC-L 1 optimisation problem, we take into account the general nonlinear MPC-L 1 optimisation task defined by Equation (5) in which the first part of the minimised cost function is approximated by Equation (22). As a result, we obtain: Thanks to linearisation, the minimised cost function is quadratic in terms of the decision vector, u(k), and all the constraints are linear with respect to the vector u(k). The auxiliary vectors of length N u are: The obtained optimisation problem (25) can now be transformed into the standard form, typical of quadratic optimisation tasks: where the inequality constraints are defined by the matrix: and the vector: while the box constraints imposed on the decision vector are: By differentiating the cost function J 1 (k) used in the quadratic optimisation problem (25) with respect to the decision variables, u(k), we obtain: Hence, the Hessian matrix necessary in the classical form of the quadratic optimisation task (27) is: while the vector f QP is: Calculation of the matrix of derivatives, as well as the trajectory, are explained for the specific form of the model used in simulations presented in Section 4. In general, two versions of the presented algorithm are possible. Firstly, in the MPC algorithm with nonlinear prediction and linearisation along the trajectory (MPC-NPLT-L 1 ), trajectory linearisation may be executed once at each sampling instant for the assumed trajectory u traj (k). It means that only one quadratic optimisation problem is solved at every sampling instant. Alternatively, in the MPC algorithm with nonlinear prediction and linearisation along the predicted trajectory (MPC-NPLPT-L 1 ), a few internal iterations are possible at each sampling instant. The first internal iteration is the same as in the MPC-NPLT-L 1 scheme. In the consecutive ones, the optimal solution obtained in the previous internal iteration is used for linearisation.

The Neutralisation Reactor
In this work, the neutralisation reactor is used as a benchmark to evaluate and compare all considered MPC methods. The process has one manipulated variable, which is a base (NaOH) stream q 1 (mL/s) and one controlled variable, which is the value of pH of the product. The detailed fundamental model of the process is given in [66]. The process is nonlinear, since its static and dynamic properties depend on the operating point. Hence, it is frequently used as a good benchmark to evaluate model identification algorithms and advanced nonlinear control methods.

Neutralisation Reactor Modelling for MPC
In this work, a Wiener model [67] of the neutralisation reactor is used in MPC for prediction. Such a block-orientated model is composed of two parts: a linear dynamic block followed by a nonlinear static one. The first part of the Wiener model is described by second-order dynamics: The auxiliary variable between two model blocks is denoted by v. The second part of the model is represented by a differentiable function: Combining Equations (34) and (35), we obtain the output of the Wiener model for the sampling instant k: The relations between the process input and output variables, i.e., q 1 and pH, respectively, and the input and output of the model, i.e., u and y, respectively, are given by the following relations: where in the nominal operating point, q 1,0 = 15.5 and pH 0 = 7. In this study, a sigmoidlike neural network is used to represent the nonlinear static part of the model. Details of model training, validation and selection are given in [68]. An alternative, a Support Vector Machine (SVM) Wiener model, is presented in [69]. The sampling time of the Wiener model is 10 s.

Calculation of the Predicted Trajectories for the Wiener Model of the Neutralisation Reactor
The trajectory α(e traj (k)), defined by Equation (19), is calculated using the neural approximator from Equation (10). The predicted control errors, e(k + p|k), are computed from Equation (6). The predicted values of the controlled variable,ŷ traj (k + p|k), are calculated from the Wiener model in the following way. From the linear block of the model (Equation (34)), we have: From the nonlinear block of the model (Equation (36)), the output predictions are: for any sampling time instant k + p|k, where p = 1, . . . , N. The disturbance estimate is calculated using Equation (36), which yields:

Calculation of the Matrices of Derivatives for the Wiener Model of the Neutralisation Reactor
The elements of the matrix dα(e traj (k)) du traj (k) defined by Equation (21) can be found using Equation (16). The derivatives ∂ŷ traj (k+p|k) ∂u traj (k+r|k) for all p = 1, . . . , N and r = 0, . . . , N u − 1 are calculated for the Wiener model used in the following way. By differentiating Equation (41), the following formula is obtained: The first derivative in Equation (43) depends on the type of nonlinear static function. The second is calculated recursively for k + 1|k, . . . , k + N|k. For the first sampling time instant, k + 1, calculated at the current time instant, k; by differentiating the formula (38), the following formula is acquired: Predictionŷ(k + 1|k) does not depend on control signals u(k + 1|k), u(k + 2|k), . . .. With that in mind, the following formula is obtained: For the next time instant, k + 2, by differentiating Equation (39), the following is acquired: where the derivatives ∂u traj (k+p|k) ∂u traj (k+r|k) can only take the values 0 or 1: For the time instant k + p|k, where p = 3, . . . , N, by differentiating Equation (40), we obtain: Analogously, using Equation (45), a general regularity can be observed:

1.
For the Wiener model, the disturbance estimate is calculated from Equation (42).

2.
The trajectory of the manipulated variable, u traj (k), that defines the linearisation point (Equation (11)), is formed. Three possible choices are discussed in the next section.
In the case of the MPC-NPLT-L 1 algorithm, the first element of the obtained decision vector, u opt (k), is applied to the process, i.e., u(k) = u opt (k|k) + u(k − 1). 7.
In the case of the MPC-NPLPT-L 1 algorithm, steps 2-5 are repeated a few times (in this work, maximally five times). The trajectory used for linearisation is defined as u traj (k) = J u opt (k) + u(k − 1), where the matrix J and the vector u(k − 1) are defined by Equations (23) and (24), respectively, and u opt (k) denotes the optimal solution calculated in the previous internal iteration (for the current sampling instant k). When the internal iterations are terminated, the first element of the decision vector computed in the last internal iteration is applied to the process.

Comparison of MPC-L 1 and MPC-L 2 Algorithms for the Neutralisation Reactor
At first, let us verify the usefulness of the MLP neural network to serve as an approximator of the absolute value function of the predicted control error (Equation (10)). The neural network with K = 10 hidden nodes of the tanh type is used. Figure 1 compares the non-differentiable absolute value (abs) function and its differentiable neural approximation. The range of the control error, e, is adequate for further use of the neural approximator in MPC for the considered neutralisation reactor. For the chosen neural network structure and the number of hidden nodes, the obtained approximation accuracy is very good. In the following part of the article, the following MPC algorithms are considered: • MPC-NO-L 1 : the MPC algorithm with nonlinear optimisation with the L 1 norm used in the first part of the minimised cost function defined by Equation (4). The resulting nonlinear optimisation task is given by Equation (5). Two versions of the MPC-NO-L 1 are considered: the non-differentiable absolute value function or its differentiable neural approximation may be used. • MPC-NPLT1-L 1 : the discussed MPC algorithm with nonlinear prediction and linearisation along the trajectory. The neural network is used to approximate the nondifferentiable absolute value function. Moreover, a linear approximation of the nonlinear trajectory of the predicted control errors over the prediction horizon is used in the cost function. The resulting quadratic optimisation task is given by Equation (25). The trajectory used for linearisation, i.e., u traj (k) (Equation (11)) is constant; all its elements are equal to the value of the manipulated variable calculated at the previous sampling instant, i.e., u(k − 1), and applied to the process. • MPC-NPLT2-L 1 : the trajectory used for linearisation is defined by the last N u − 1 elements of the optimal solution u(k) computed at the previous sampling instant.
Only the first element of this sequence is actually used for control. • MPC-NPLT3-L 1 : the trajectory used for linearisation is constant, all its elements are equal to the value of the process input corresponding to the current output set-point. For this purpose, the inverse static model of the process is used: u sp (k) =g(y sp (k)).
In this work, a neural network of the MLP type with two layers serves as the inverse model (the first nonlinear layer contains 10 hidden nodes of the tanh type). • MPC-NPLPT-L 1 : the discussed MPC algorithm with nonlinear prediction and linearisation along the predicted trajectory. In this case, trajectory linearisation and quadratic optimisation are repeated maximally five times at each sampling instant. The trajectory used for linearisation is taken from the previous internal iteration of the algorithm. In the first internal iteration, for linearisation, the trajectory obtained from the inverse static model for the current set-point is used, exactly as it is done in the MPC-NPLT3-L 1 scheme.
Similarly, the following MPC algorithms with the L 2 cost function (3) are considered: • MPC-NO-L 2 : the MPC algorithm with nonlinear optimisation with the L 2 norm used in two parts of the minimised cost function. The resulting nonlinear optimisation task is given by Equation (2).
• MPC-NPLT1-L 2 , MPC-NPLT2-L 2 and MPC-NPLT3-L 2 : the MPC algorithm with nonlinear prediction and linearisation along the trajectory [29]. Trajectory linearisation and quadratic optimisation are performed once at each sampling instant. • MPC-NPLPT-L 2 : the MPC algorithm with Nonlinear Prediction and Linearisation along the Predicted Trajectory [29]. Trajectory linearisation and quadratic optimisation are repeated maximally five times at each sampling instant.
In all simulations presented next, the same parameters of MPC are used: N = 10, N u = 3, λ = 1. The constraints are imposed on the magnitude of the manipulated variable: q min 1 = 0, q max 1 = 30. All simulations are performed in MATLAB. The fmincon and quadprog functions are used for nonlinear and quadratic optimisation, respectively; default parameters are used in both cases.
Firstly, let us evaluate the efficiency of two versions of the MPC-NO-L 1 control scheme. In both of them, nonlinear optimisation is used at each sampling instant on-line, but the objective of the comparison is to only check the efficiency of the neural approximator of the absolute value function. In the first case, the non-differentiable absolute value function is used in the first part of the minimised cost function, while in the second case, a neural approximator is considered. Figure 2 depicts the obtained results for a few changes of the set-point. In general, the non-differentiable absolute value function results in some numerical problems for the optimisation routine. The influence of such problems on the resulting trajectories of the process is best visible for the second step of the set-point. For better comparison, two bottom panels of Figure 2 show an enlarged fragment of the obtained results. For the non-differentiable absolute value function, the sign of the real control error changes a few times, while the neural approximator smooths the trajectory of the process output.  Figure 3 compares the simulation results of three computationally efficient MPC algorithms with one repetition of trajectory linearisation and quadratic optimisation at each sampling instant, but the trajectory used for linearisation, i.e., u traj (k), is chosen in different ways. In all algorithms, the L 1 norm cost function is used. The first two initialisation methods, used in the MPC-NPLT1-L 1 and MPC-NPLT2-L 1 algorithms, lead to some problems for the second set-point step. The explanation is the following: when the set-point changes abruptly, the linear approximation of the predicted trajectory performed using the past operating point (MPC-NPLT1-L 1 ) or the part of the optimal trajectory for the past operating point (MPC-NPLT2-L 1 ) differs significantly from the true nonlinear trajectory, performed in the MPC-NO-L 1 algorithm without any linearisation ( Figure 2). Of course, increasing the value of the coefficient λ would solve the problem, but it leads to slower control. The best results are obtained in the MPC-NPLT3-L 1 algorithm, in which the trajectory used for linearisation is constant, and all its elements are equal to the value of the process input corresponding to the current output set-point. Such an input value is calculated using an inverse neural static model of the process u sp (k) =g(y sp (k)).
Let us remind ourselves that our objective is to obtain an MPC algorithm that uses the L 1 cost function but it should be computationally efficient, i.e., quadratic optimisation should be used in place of nonlinear optimisation. On the other hand, the "ideal" MPC-NO-L 1 algorithm is treated as the reference. In different words, we develop a method that is computationally much simpler, but we hope to obtain control accuracy similar to that possible in the reference MPC-NO-L 1 algorithm. Figure 4 compares the best computationally efficient MPC algorithms with one repetition of trajectory linearisation and quadratic optimisation at each sampling instant, i.e., the MPC-NPLT3-L 1 scheme, with the reference MPC-NO-L 1 scheme. We can see that the MPC-NPLT3-L 1 algorithm discussed in this work gives very good trajectories, very similar to those obtained in the MPC-NO-L 1 control approach. However, there are some small discrepancies between the obtained trajectories. The maximal allowed number of such repetitions is five. The stopping criteria (continuation and termination of the internal iterations) are defined by an additional parameter δ [29]. The internal iterations are continued if the sum of squared differences between the prediction of the model output and the required set-point is greater than δ. The internal iterations are terminated if the sum of squared differences between the optimal solution obtained in two consecutive internal iterations is lower than δ. In general, the lower the parameter δ, the better control quality we expect to obtain, but the number of repetitions of the internal iterations increases. Figure 5 shows the obtained trajectories of the MPC-NPLPT-L 1 algorithm for three values of the parameter δ. For the considered neutralisation process, the obtained differences are small, but the best results are obtained for δ = 10 −1 , as it gives the best trajectory for the second set-point step. The considered neutralisation process is nonlinear. The process gain and its dynamics heavily depend on the current operating points. Hence, for different operating points, different overshoot and settling times are observed. Figure 6 depicts the number of internal iterations in the consecutive sampling instants of the MPC-NPLPT-L 1 algorithm for three values of the parameter δ. More than one internal iteration is necessary when the set-point changes; when the process is close to the required operating point, one internal iteration is sufficient. It is clear that the lower the parameter δ, the more internal iterations are necessary. This is true provided that a perfect model is used in MPC. In reality, the model is an approximation of the process. In our case, the unavoidable process-model mismatch results in the lowest overshoot for δ = 10 −1 . Remembering our objective the most important issue is to compare the trajectories of the reference MPC-NO-L 1 algorithm and the computationally efficient MPC-NPLPT-L 1 scheme in the latter one, δ = 10 −1 . Figure 7 shows the obtained trajectories. It is very important to note that the recommended MPC-NPLPT-L 1 algorithm with on-line trajectory linearisation and quadratic optimisation gives very similar, practically the same, trajectories as the algorithm with nonlinear optimisation. It is necessary to recall Figure 4, which compares the performance of the MPC-NPLT3-L 1 algorithm in which linearisation and optimisation are performed only once at each sampling instant. Comparing Figures 4 and 7, it is clear that multiple linearisation and optimisation possible at each sampling instant in the MPC-NPLPT-L 1 approach is really beneficial for the considered process. In Figure 8, the recommended MPC-NPLPT-L 1 algorithm, which uses the L 1 norm, is compared with its counterpart, MPC-NPLPT-L 2 , in which the classical L 2 norm is used. The use of the L 1 norm leads to a lower overshoot and shorter setting time. Having visually compared the obtained trajectories of the discussed MPC algorithms, it is interesting to analyse their performance and differences using some numerical control quality indicators [70]. We consider the following indices: • The sum of absolute values of control errors for the whole simulation horizon defined as: where y(k) denotes the real value of the process output obtained in simulation. • The sum of absolute values of differences between the output of the process when it is controlled by the "ideal" MPC-NO-L 1 algorithm (y MPC-NO-L 1 (k)) and the output of the process when it is controlled by a compared MPC scheme (y(k)). These differences are considered for the whole simulation horizon: • The sum of squared control errors for the whole simulation horizon defined as: • The sum of squared differences between the output of the process when it is controlled by the "ideal" MPC-NO-L 1 algorithm and the output of the process when it is controlled by a compared MPC scheme. These differences are considered for the whole simulation horizon: It is important to stress the fact that the performance criterion E 1 is not directly minimised in MPC algorithms that use the L 1 cost function J 1 (Equation (4)). This is because the actually minimised cost function J 1 takes into account the absolute value of control errors predicted over the prediction horizon. The horizon is shifted at each sampling instant. Moreover, it takes into account the second part of the cost function in which the sum of squared moves of the manipulated variable of the control horizon is penalised. Similarly, the performance criterion E 2 is not directly minimised in MPC algorithms that use the L 2 cost function J 2 (Equation (3)). Let us state our expectations and objectives. When the norm L 1 is used in MPC, we expect to obtain the lowest value of the indices E 1 and E MPC-NO-L 1 1 . The more similar the obtained trajectories are to those obtained in the MPC-NO-L 1 scheme, the closer the index E MPC-NO-L 1 1 is to 0. Table 1 presents numerical values of the control quality indices (50)- (53). The following algorithms are considered: MPC-NPLT1, MPC-NPLT2 and MPC-NPLT3 algorithms, three versions of the MPC-NPLPT scheme with different stopping criterion defined by the parameter δ and the MPC-NO approach. The results are divided into two parts: in the first one, the norm L 1 is used in MPC; in the second one, the norm L 2 is considered; the cost function type minimised in MPC is denoted in colour. In addition to the control quality indices, Table 1 also specifies the relative calculation time of all algorithms; the results are given in percentages in such a way that the calculation time for the most demanding algorithm, i.e., MPC-NO-L 1 , is treated as 100%. Considering the obtained numerical results, we are able to formulate the following observations concerned with control quality: 1.
Comparing the MPC algorithms with the norm L 1 , in which one trajectory linearisation and quadratic optimisation are executed at each sampling instant, the best results are obtained in the MPC-NPLT3-L 1 scheme, in which the trajectory linearisation is performed using an inverse static model of the process. That algorithm gives the lowest values of the performance indices E 1 and E MPC-NO-L 1 1 . It confirms the comparison given in Figure 3 are lower. It confirms the comparison given in Figure 7. Moreover, the lower the parameter δ, the more similar the obtained trajectory is to that possible in the reference MPC-NO-L 1 scheme.

3.
Bearing in mind our expectations and objectives, the classical MPC algorithms that use the L 2 norm give a worse performance. This confirms the comparison given in Figure 8. , which is natural, but also makes it possible to reduce the indices E 2 and E MPC-NO-L 1 2 . For all pairs of algorithms (with L 1 and L 2 norms), the E 2 index is slightly lower when the L 1 norm is used. This difference is even more clear when we consider the E MPC-NO-L 1 2 index. It confirms the comparison given in Figure 8. Furthermore, it is also possible to put forward the following observations concerned with computation time:

1.
In general, all MPC algorithms with the norm L 1 are more computationally demanding than their counterparts that use the norm L 2 . This is because, in the first case, in all calculations, i.e., in prediction, linearisation and optimisation, the neural approximator determines the absolute values of the control errors over the prediction horizon whereas, in the second case, no approximator is used, the predictions and control errors are used directly in all calculations.

2.
All MPC algorithms with linearisation and quadratic optimisation are less computationally demanding than the reference "ideal" MPC-NO algorithm.

3.
The more complicated the trajectory linearisation, the longer the calculation time. The lowest calculation time is observed in the MPC-NPLT1, MPC-NPLT2 and MPC-NPLT3 algorithms, with one repetition of linearisation and quadratic optimisation at each sampling instant. The calculation time becomes longer in the MPC-NPLPT scheme, with a few repetitions of linearisation and optimisation at each instant; the lower the parameter δ, the longer the calculation time.
Finally, let us compare the results of the MPC-NPLPT-L 1 and MPC-NPLPT-L 2 algorithms when the penalty coefficient is increased. In our study, the value λ = 5 is considered. Figure 9 presents the obtained trajectories. Two top panels show the results for the whole simulation scenario, whereas the bottom ones show enlarged fragments for three set-point changes. It is interesting to see an additional advantage of the recommended L 1 norm because, in such a case, the overshoot is much lower and the setting time is much shorter when compared with the trajectories possible when the classical L 2 norm is used. In the considered comparison, for λ = 5, the benefits of using the L 1 norm are even more clear than in the case of the default value of the penalty coefficient, i.e., λ = 1, as presented in Figure 8.

Conclusions
The presented MPC-L 1 algorithms have three essential advantages. Firstly, thanks to using a neural differentiable approximator of the non-differentiable absolute value function and on-line advanced trajectory linearisation, computationally simple quadratic optimisation is used in place of demanding nonlinear optimisation. Secondly, the obtained trajectories are very similar, practically the same, as those possible in the reference scheme with nonlinear optimisation. Thirdly, it is shown that the use of the recommended L 1 norm gives better results defined using different control quality criteria, such as the sum of absolute values or squared control errors, overshoot and setting time. Furthermore, it must be stressed that the L 1 norm even gives better results than the classical L 2 one in terms of the classical control performance indicator that measures squared control errors. The discussed MPC is very universal, it may be used in industrial process control applications and in fast embedded systems. The only condition is that the model used for prediction and the neural approximator must be differentiable.
As future works, the following issues are worth investigating: taking into account not only constraints imposed on the manipulated variable, but also on the predicted controlled one, the development of nonlinear MPC-L 1 algorithms for multivariable processes with multiple inputs and multiple outputs, considering state-space process description in place of the input-output one and considering alternative types of the cost function, not only the discussed L 1 norm. Funding: This research was financed by Warsaw University of Technology in the framework of the project for the scientific discipline automatic control, electronics and electrical engineering.