# Filled Function Method for Nonlinear Model Predictive Control

Authors: Hajer Degachi; Bechir Naffeti; Wassila Chagra; Moufida Ksouri
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

A new method is used to solve the nonconvex optimization problem of the nonlinear model predictive control (NMPC) for Hammerstein model. Using nonlinear models in MPC leads to a nonlinear and nonconvex optimization problem. Since control performances depend essentially on the results of the optimization method, in this work, we propose to use the filled function as a global optimization method to solve the nonconvex optimization problem. Using this method, the control law can be obtained through two steps. The first step consists of determining a local minimum of the objective function. In the second step, a new function is constructed using the local minimum of the objective function found in the first step. The new function is called the filled function; the new constructed function allows us to obtain an initialization near the global minimum. Once this initialization is determined, we can use a local optimization method to determine the global control sequence. The efficiency of the proposed method is proved firstly through benchmark functions and then through the ball and beam system described by Hammerstein model. The results obtained by the presented method are compared with those of the genetic algorithm (GA) and the particle swarm optimization (PSO).


---

## Body

## 1. Introduction

Model predictive control (MPC) is a powerful control tool. This control strategy is formulated as the repeated solution of an open loop control problem [1]. At each sampling time, a control sequence is determined where only the first component is applied to the system and the optimization task is repeated at the next sampling instant.The MPC algorithm presents the major advantage to efficiently handle constraints on input and output [2–4]. MPC is also able to control a wide variety of processes starting from systems that present a simple behavior like linear process [5] as well as those that exhibit more complex behavior like nonlinear [6, 7] and multivariable process [8].The MPC theory was successfully applied in many areas of application such as chemical, petrochemical, pulp and paper, aerospace and defense, and food processing [9].Firstly, MPC was well developed for linear models. Since almost real processes have a nonlinear behavior, this fact motivates the development of the NMPC strategy [4].The NMPC is considered as a purely optimization-based algorithm. In fact, the minimization problem of the NMPC problem is nonlinear and nonconvex due to the nonlinear nature of the model. The determination of the control sequence, solution of the NMPC optimization problem, should be done using a global optimization algorithm to ensure good control performances. For this, a variety of solutions were proposed in literature to solve this kind of problem. Many works were focused on online linearization of the nonlinear model [10–12]. A second solution was proposed in some other works that consists of reversing the nonlinear block to remove nonlinearity. This solution is only addressed to block-oriented models [13] where the nonlinear block is described by a polynomial. The two mentioned solutions allow obtaining a quadratic cost function at each sampling time where the global solution can be easily determined.Since solving a nonlinear optimization problem at each sampling time is a hard-computational task, [14] has proposed to describe the nonlinear process by a set of uncertain linear models to overcome the online computational complexity of the NMPC. Since the used model in the NMPC is nonlinear and since real process must operate under rigorous conditions, the nonlinear characteristic of the model should be kept without simplification during the online computation of the control sequence to get performances with sufficient accuracy.For this, genetic algorithm (GA) [15], particle swarm optimization (PSO) [16], and neural network (NN) [17] were used as global optimization methods to solve the nonlinear optimization problem. Also, deterministic method such as Generalized Geometric Programming (GGP) was used as global optimization method for the NMPC [18, 19].In this work, we will use the filled function as a global deterministic optimization method to solve the NMPC optimization problem. This method is based on the construction of a new function starting from the first local minimum found of the original function to be minimized and carry out the search of the global minimum using a local optimization method. This function allows determining the global minimum by finding an appropriate initialization.In this work, the NMPC strategy is based on Hammerstein model. The Hammerstein model belongs to block-oriented models that is described by a nonlinear static block followed by a linear dynamic one. The special structure of the Hammerstein model was widely exploited to describe a variety of processes such as heat exchanger [17], SOFC [20], and pH neutralization process [21].The rest of the paper is organized as follows. In Section2, we describe the Hammerstein model. The NMPC based on Hammerstein model is presented in Section 3. In Section 4, we outline the optimization method used to solve the NMPC minimization problem. Simulation results on benchmark functions based on the global optimization method are given in Section 5. The efficiency of the proposed method for NMPC is illustrated in Section 6.


## 2. Hammerstein Model

The Hammerstein model is composed of a static nonlinear block followed by a linear dynamic one as depicted in Figure1.Figure 1
Hammerstein model.The outputs of the two blocks are given by the following equation:(1)S:sk=fukAq-1yk=Bq-1skThe polynomialsAq-1 and Bq-1 are defined as(2)Aq-1=1+a1q-1+a2q-2+…+anaq-na(3)Bq-1=b1q-1+b2q-2+…+bnbq-nbwhereq-1 is the unit delay operator and na and nb define the order of polynomials Aq-1 and Bq-1.Various types of nonlinearity can be used to describe the nonlinear block of the Hammerstein model. In this work, the polynomial form is chosen, and the functionf is defined as(4)fuk=c1uk+c2uk2+⋯+cr1ukncwherenc is the order of the polynomial.


## 3. Nonlinear Model Predictive Control Design

The aim of the NMPC strategy is to compute the control sequenceU=uk,uk+1,…,uk+Nu-1 by minimizing the criterion J defined as the difference between the desired trajectory ysp(k+j) and the model output y^(k+j) over the prediction horizon Np and penalizing the control increments over the control horizon Nu(5)J=∑i=1Npyspk+i-y^k+i2+λ∑i=0Nu-1Δuk+i2(6)Δumin≤Δuk+j≤Δumax,j=0,…,Nu-1umin≤uk+j≤umax,j=0,…,Nu-1whereΔuk+i=uk+i-uk+i-1 and λ is a positive weighting coefficient.Based on the receding horizon principle, only the first component will be applied to the system. Then the horizon will be shifted one step forward and the whole procedure will be repeated.The predicted output of the Hammerstein model is defined as(7)y^k+j=-∑i=1naaiy^k-i+j+∑i=1nbbis^k-i+js ^ ( k + j ) is the output of the nonlinear block defined as(8)s^k+j=∑l=1nccluk+jlBy using (7) and (8), the criterion expression (5) becomes(9)Jk=∑j=1Npyspk+j+∑i=1naaiy^k-i+j-∑i=1nbbi∑l=1nccluk-i+jl2+λ∑j=0Nu-1Δuk+j∣k2Due to the nonlinear nature of Hammerstein model, the obtained objective function defined by (9) of the NMPC problem is nonconvex.Solving the presented minimization problem using standard optimization methods will lead to suboptimal results in control. Deterministic optimization methods have the major benefit to offer the global solution. For this, we propose in this work to use the filled function method to solve the NMPC optimization problem. The present method was found to be easier to implement and offer an accurate solution but suffer from the high computation time.


## 4. Optimization Method

In this work, we will use the filled function method to solve the NMPC optimization problem.The concept of this global optimization method consists of determining a good initialization to obtain the global minimum. Therefore, the key of this method is to construct a function “P’’ that allows moving from the local minimum of the objective function to obtain a more promising initialization point. This point presents a stationary point of the function P. The new initialization will allow finding easily the global minimum.Thus, we consider a functionf(x) to be minimized. The principle of the filled function method can be summarized as follows [22].Step 1.
Use a local optimization method to find a local minimumx1 of f(x).Step 2.
Construct the filled functionP(x) at the local minimum x1 of f(x). Minimize the new function P using a local optimization method and so the minimum, t1, of the filled function is found. The minimum t1of the filled function is used as initialization to minimize again the functionf, where a new better minimum will be found. Many filled functions are proposed in literature. Most of them are described using some parameters that can be difficult to be adjusted [23]. Other works minimize this number of parameters and reduce the filled function into a mono adjustable parameter function.
In this work, we will use the nonparametric filled function method proposed by [22] and defined as(10)Pt,x1=-∫x1xft-fx1dtx≥x1-∫xx1ft-fx1dtx≤x1
The nonparametric filled functionP should satisfy the following properties.Theorem 1.
Ifx1is an isolated minimizer of f:Ω⊂R→R, then x1must be a maximizer of P(x,x1).Theorem 2.
Ifx1is an isolated minimizer of f:Ω⊂R→R, then P(x,x1) does not have any stationary point for fx>fx1.Theorem 3.
Ifx1is an isolated minimizer of f:Ω⊂R→R, then Px,x1 must have a stationary point for fx≤fx1. The flowchart of the considered method is given in Figure 2 that details the different steps of the algorithm. The variable α in the algorithm presents a small step within the interval [0,1].Figure 2
Flowchart of the filled function algorithm.


## 5. Simulation Results on Benchmark Function

The efficiency of the filled function is demonstrated through benchmark functions defined in Table1. We consider 3 multimodal functions as depicted in Figures 3, 4, and 5. The function f1 has a global minimum atx = 17.0392 and two local minima atx = 5.3622 andx = 10.4535. The function f2 has a global minimum atx = 5.1994 and two local minima atx = 3.4392 andx = 7.0678. The function f3 has a global minimum atx = -1 and a local onex = 2. The minimization is carried out using the filled function and the local method illustrated by the “fmincon” function of the MATLAB environment.Table 1
Benchmark functions.
Function Min Research interval f 1 x = sin ⁡ x + sin ⁡ 2 x 3 17.03 [ 0,20 ]  f 2 x = sin ⁡ x + sin ⁡ 10 x 3 + ln ⁡ x - 0.84 x 5.36 [ 2.7,7.5 ]  f 3 x = x 4 - 3 x 3 - 1.5 x 2 + 10 x -1 [ - 2,3 ]Figure 3
Evolution of the functionf1.Figure 4
Evolution of the functionf2.Figure 5
Evolution of the functionf3.Consider the first functionf1.Iteration 1.Step 1.
The first step in the filled function optimization method consists of using a local optimization method to determine a local minimum off. We choose as initialization x0=2.
In this work, we use the gradient descent method to find a local minimumx1 of f. Using x0=2 as initialization and using the gradient we can get a local minimum, x1 =5.3622, as indicated in Figure 6Figure 6
Evolution off1 and P in the first iteration.Step 2.
Usingx1 we construct the filled function P(t,x1) using (10) as illustrated in Figure 6.Step 3.
Once we haveP, we apply a local optimization method to determine a local minimum ofP, and we set initialization as t0=x1+α, where α is a small step real value. The obtained minimum of P is t1Iteration 2.
Ift1 belongs to the interval of variation of the manipulated variable, it will be used in the second iteration as initialization forf.Step 1.
Since this condition is true in our case,t1 presents a new initialization for f as depicted in Figure 6.
The new initialization presents a promising point to get a best minimum than that found in the first step.
The new minimum,x1, of f found with initialization t1 is shown in Figure 7.Figure 7
Evolution off1 and P in the second iteration.Step 2.
Using the new found minimumx1 indicated in Figure 7 we construct the filled function P(t,x1) using (10) as illustrated in Figure 7. We can conclude from Figure 7 that P does not have a stationary point in the considered interval. So, the global minimum of f1 is x1=17.0328.The minimization results for the considered function,f1, f2, and f3, are illustrated, respectively, in Tables 1, 2, and 3Table 2
Minimization off1using global and local optimization method.
Initialization Reached min(proposed method) Reached min(local method) 3 17.0328 5.3622  7 17.0328 5.3622  9 17.0328 10.4535  14 17.0328 17.0328Table 3
Minimization off2using global and local optimization method.
Initialization Reached min(proposed method) Reached min(local method) 4.5 5.1995 3.4392  6 5.2006 7.0678  6.5 5.1994 5.1998  7 5.1994 7.0697Based on results given by Tables2, 3, and 4, we can note that, for different initialization values, the convergence to the global minimum using the local optimization method is not always guaranteed. Contrariwise, we can conclude that whatever the value of the initialization the minimization of the considered function using the proposed method is always guaranteed.Table 4
Minimization off3using global and local optimization method.
Initialization Reached min(proposed method) Reached min(local method) 2.5 -1 -1  1.5 -1 2  -1.5 -1 -1


## 6. Simulation Results

The system under consideration is a ball and beam which is composed of ball, ball position sensor, a center pivoted beam on which the ball rolls on, beam angle sensor, and the servomotor. The control objective is to fix the ball on a desired position. The ball and beam system is described by a Hammerstein model defined as [24](11)S:yk=-a1yk-1+b1sk-1sk=c1uk+c2uk2The different parameters of the linear and nonlinear blocks are given bya1=-0.62, b1=1, c1=0.396, and c2=-0.036.The prediction and the control horizons are fixed toNp = 2, Nu = 1 and constraints on the manipulated variable u are fixed as -5≤u(k)≤5.First, the NMPC algorithm for the ball and beam system is implemented using the filled function method to solve the minimization problem. The desired trajectory is fixed as -5 and 5 represented by three step changes. The NMPC algorithm based on the proposed method is compared to a local optimization method illustrated by “fmincon” function of the MATLAB environment. The simulation results are depicted in Figure8.Figure 8
Output tracking using the proposed method and a local method forλ=0.1.We can note from Figure8 that the NMPC performances based on the filled function method are significantly better than those of the NMPC based on a local optimization method. The good ability of the filled function method to solve the optimization problem and offer the optimal control quality is justified by the fact that the NMPC based on the proposed method ensures the setpoint tracking with zero steady state error. Also, we can note that the response of the system based on the local optimization method is characterized by a high overshoot compared to the results given by the response of the system based on the filled function method.In the sequel, the NMPC based on the proposed method is compared with the NMPC based on GA and PSO algorithms.To compare the performances of the different algorithms we consider two performance indexes: the SSE and the SCV.The SSE is defined as(12)SSE=1N∑i=1Nyspk-yk2Since our goal in control is to reach the desired trajectory and ensure the good output tracking by determining the optimal control sequence that gives the minimum dissipated energy, we will consider the second performance index SCV defined as(13)SCV=1N∑k=1Nuk2We can note from Figure9 that the three algorithms are able to reach the desired trajectory and ensure good output tracking as proved by the SSE value. The significant difference lays in the control quality. Based on Figure 9 we can note than each algorithm reaches the desired trajectory with a different control sequence. Table 5 shows that the best control quality is given by the proposed method since it offers the small SCV. Also, we can remark that the filled function method offers the less overshoot as illustrated by Table 5. We can conclude that for the proposed method the setpoint is reached with the best control quality and the less overshoot. Although the results of setpoint tracking obtained by the GA, PSO, and the proposed method are quite similar, it is very clear that the control quality of GA and PSO is inferior to the filled function method. In fact, setting parameters of evolutionary algorithms as coding of individuals, population size, crossover and mutation rates, and stopping criterion has an important role in the determination of an accurate solution [25].Table 5
The performance of the NMPC based on the three algorithms forλ=0.1.
Method Overshoot SCV SSE Proposed method 0.0142 1.635 3.114 GA 0.0178 10.729 3.120 PSO 0.0178 5.990 3.112Figure 9
Output tracking using the GA, PSO, and the proposed method forλ=0.1.To show the effectiveness of our proposed method with different values ofλ we suggest varying λ and carry out simulation for λ=0.2.We can notice from Figure10 and Tables 5 and 6 that the increase of the weighting factor λ results in a high overshoot value and slow system response. Based on Table 4 we can conclude that the SCV offered by the filled function method is better than that given by GA and PSO. This confirms the good ability of the filled function method to give the best control quality.Table 6
The performance of the NMPC based on the three algorithms forλ=0.2.
Method Overshoot SCV SSE proposed method 0.0494 1.646 3.114 GA 0.0594 5.994 3.123 PSO 0.0594 10.681 3.124Figure 10
Output tracking using the GA, PSO, and the proposed method forλ=0.2.Tables5 and 6 give the performance of the different algorithms for λ=0.1 and λ=0.2. We can note that for different values of λ the less overshoot and SCV are given by the filled function method.


## 7. Conclusion

In this work, a nonlinear model predictive control for Hammerstein model is presented. Using nonlinear model in the MPC strategy complicates the solution of the optimization problem and leads to computational difficulties. Consequently, an efficient optimization algorithm should be used to overcome these difficulties and ensure good control results. The filled function, as a global optimization method, is used in this work to solve the nonconvex minimization problem of the NMPC. The proposed method proves its efficiency to give good control performances through nonlinear process described by a ball and beam system. It offers a good output tracking and a less overshoot. The filled function method is then compared with two heuristic global optimization algorithms: GA and PSO. We can remark that the GA and the PSO algorithms are quite similar. We can also conclude that the NMPC with the proposed method is much better than the GA and PSO since it offers less overshoot and gives the best control sequence at each sampling time.



---

*Source: 9497618-2018-06-04.xml*