List of the model parameters

## 1. Introduction

Model predictive control (MPC) has made a significant impact on control engineering. It has been applied in almost all of industrial fields such as petrochemical, biotechnical, electrical and mechanical processes. MPC is one of the most applicable control algorithms which refer to a class of control algorithms in which a dynamic process model is used to predict and optimize process performance. Linear model predictive control (LMPC) has been successfully used for years in numerous advanced industrial applications. It is mainly because they can handle multivariable control problems with inequality constraints both on process inputs and outputs.

Because properties of many processes are nonlinear and linear models are often inadequate to describe highly nonlinear processes and moderately nonlinear processes which have large operating regimes, different nonlinear model predictive control (NMPC) approaches have been developed and attracted increasing attention over the past decade [1-5].

On the other hand, since the incorporation of nonlinear dynamic model into the MPC formulation, a non-convex nonlinear optimal control problem (NOCP) with the initial state must be solved at each sampling instant. At the result only the first element of the control policy is usually applied to the process. Then the NOCP is solved again with a new initial value coming from the process. Due the demand of an on-line solution of the NOCP, the computation time is a bottleneck of its application to large-scale complex processes and NMPC has been applied almost only to slow systems. For fast systems where the sampling time is considerably small, the existing NMPC algorithms cannot be used. Therefore, solving such a nonlinear optimization problem efficiently and fast has attracted strong research interest in recent years [6-11].

To solve NOCP, the control sequence will be parameterized, while the state sequence can be handled with two approaches: sequential or simultaneous approach. In the sequential approach, the state vector is handled implicitly with the control vector and initial value vector. Thus the degree of freedom of the NLP problem is only composed of the control parameters. The direct single shooting method is an example of the sequential method. In the simultaneous approach, state trajectories are treated as optimization variable. Equality constraints are added to the NLP and the degree of freedom of the NLP problem is composed of both the control and state parameters. The most well-known simultaneous method is based on collocation on finite elements and multiple shooting.

Both single shooting method and multiple shooting based optimization approaches can then be solved by a nonlinear programming (NLP) solver. The conventional iterative optimization method，such as sequential quadratic programming (SQP) has been applied to NMPC. As a form of the gradient-based optimization method, SQP performs well in local search problems. But it cannot assure that the calculated control values are global optimal because of its relatively weak global search ability. Moreover, the performance of SQP greatly depends on the choice of some initialization values. Improper initial values will lead to local optima or even infeasible solutions.

Genetic Algorithms (GAs) are a stochastic search technique that applies the concept of process of the biological evolution to find an optimal solution in a search space. The conceptual development of the technique is inspired by the ability of natural systems for adaptation. The increasing application of the algorithm has been proved to be efficient in solving complicated nonlinear optimization problems, because of their ability to search efficiently in complicated nonlinear constrained and non-convex optimization problem, which makes them more robust with respect to the complexity of the optimization problem compared to the more conventional optimization techniques.

Compared with SQP, GAs can reduce the dimension of search space efficiently. Indeed, in SQP the state sequence is treated as additional optimization variables; as such, the number of decision variables is the sum of the lengths of both the state sequence and the control sequence. In contrast, in GAs, state equations can be included in the objective function, thus the number of decision variables is only the length of control sequence. Furthermore, the search range of the input variable constraints can be the search space of GA during optimization, which makes it easier to handle the input constraint problem than other descent-based methods.

However, a few applications of GAs to nonlinear MPC [12][13] can partially be explained by the numerical complexity of the GAs, which make the suitable only for processes with slow dynamic. Moreover, the computational burden is much heavier and increases exponentially when the horizon length of NMPC increases. As a result, the implementation of NMPC tends to be difficult and even impossible.

In this paper an improved NMPC algorithm based on GA is proposed to reduce the severe computational burden of conventional GA-based NMPC algorithms. A conventional NMPC algorithm seeks the exact global solution of nonlinear programming, which requires the global solution be implemented online at every sampling time. Unfortunately, finding the global solution of nonlinear programming is in general computationally impossible, not mention under the stringent real-time constraint. We propose to solve a suboptimal descent control sequence which satisfies the control, state and stability constraints in the paper. The solution does not need to minimize the objective function either globally or locally, but only needs to decrease the cost function in an effective manner. The suboptimal method has relatively less computational demands without deteriorating much to the control performance.

The rest of the paper is organized as follows. Section 2 briefly reviews nonlinear model predictive control. Section 3 describes the basics of GAs, followed by a new GA-based computationally efficient NMPC algorithm. Section 4 analyses the stability property of our nonlinear model predictive control scheme for closed-loop systems. Section 5 demonstrates examples of the proposed control approach applied to a coupled-tank system and CSTR. Finally we draw conclusions and give some directions for future research.

## 2. Nonlinear model predictive control

### 2.1. System

Consider the following time-invariant, discrete-time system with integer *k* representing the current discrete time event:

In the above,

### 2.2. Objective function

The objective function in the NMPC is a sum over all stage costs plus an additional final state penalty term [14], and has the form:

where *x*(*k*+*j|k*) and *u*(*k*+*j|k*) are predicted values at time *k* of *x*(*k*+*j*)*,u*(*k*+*j*). *P* is the prediction horizon. In general, *Q*≥*0*,*R>0*.

### 2.3. General form of NMPC

The general form of NMPC law corresponding to (1) and (2) is then defined by the solution at each sampling instant of the following problem:

where *X* _{F }is a terminal stability constraint, and u(*k*)=[*u*(*k*|*k*),…,*u*(*k*+*P*-1|*k*)] is the control sequence to be optimized over.

The following assumptions A1 - A4 are made:

A1: *X* _{F }⊂ *X*, *X* _{F }closed, 0∈ *X* _{F}

A2: the local controller *κ* _{F }(x)∈ *U*, *x*∈ *X* _{F}

A3: *f*(*x*, *κ* _{F }(x)) ∈ *X* _{F }, *x*∈ *X* _{F}

A4: *F*(*f*(*x*, *κ* _{F }(x)))-*F*(*x*)+*l*(*x*, *κ* _{F }(*x*)) ≤ 0, *x*∈ *X* _{F}

Based on the formulation in (3), model predictive control is generally carried out by solving online a finite horizon open-loop optimal control problem, subject to system dynamics and constraints involving states and controls. At the sampling time *k*, a NMPC algorithm attempts to calculate the control sequence u(*k*) by optimizing the performance index (3) under constraints (4) (5) and terminal stability constraint (6). The first input *u*(*k*|*k*) is then sent into the plant, and the entire calculation is repeated at the subsequent control interval *k*+1.

## 3. NMPC algorithm based on genetic algorithm

### 3.1. Handling constraints

An important characteristic of process control problems is the presence of constraints on input and state variables. Input constraints arise due to actuator limitations such as saturation and rate-of-change restrictions. Such constraints take the form:

State constraints usually are associated with operational limitations such as equipment specifications and safety considerations. System state constraints are defined as follows:

where ∆u(*k*)=[*u*(*k*|*k*)-*u* _{k }-1,…, *u*(*k* +*P*-1|*k*)- *u*(*k* +*P*-2|*k*)], *x*(*k*)=[*x*(*k*+1|*k*),…,*x*(*k*+*P*|*k*)].

The constraints (7) and (8) can be written as an equivalent inequality:

### 3.2. Genetic algorithm

GA is known to have more chances of finding a global optimal solution than descent-based nonlinear programming methods and the operation of the GA used in the paper is explained as follows.

#### 3.2.1. Coding

Select the elements in the control sequence u(*k*) as decision variables. Each decision variable is coded in real value and *n* _{u }**P* decision variables construct the *n* _{u }**P* -dimensional search space.

#### 3.2.2. Initial population

Generate initial control value in the constraint space described in (10). Calculate the corresponding state value sequence x(*k*) from (4). If the individual (composing control value and state value) satisfies the state constraints (9) and terminal constraints (6), select it into the initial population. Repeat the steps above until PopNum individuals are selected.

#### 3.2.3. Fitness value

Set the fitness value of each individual as 1/(*J+*1).

#### 3.2.4. Genetic operators

Use roulette method to select individuals into the crossover and mutation operator to produce the children. Punish the children which disobey the state constraints (9) and terminal constraints (6) with death penalty. Select the best PopNum individuals from the current parent and children as the next generation.

#### 3.2.5. Termination condition

Repeat the above step under certain termination condition is satisfied, such as evolution time or convergence accuracy.

### 3.3. Improved NMPC algorithm based on GA

In recent years, the genetic algorithms have been successfully applied in a variety of fields where optimization in the presence of complicated objective functions and constraints. The reasons of widely used GAs are its global search ability and independence of initial value. In this paper GAs are adopted in NMPC applications to calculate the control sequence. If the computation time is adequate, GAs can obtain the global optimal solution. However, it needs to solve on line a non-convex optimization problem involving a total number of *n* _{u }**P* decision variables at each sampling time. To obtain adequate performance, the prediction horizon should be chosen to be reasonably large, which results in a large search space and an exponentially-growing computational demand. Consequently, when a control system requires fast sampling or a large prediction horizon for accurate performance, it becomes computationally infeasible to obtain the optimal control sequence via the conventional GA approach. There is thus a strong need for fast algorithms that reduce the computational demand of GA.

The traditional MPC approach requires the global solution of a nonlinear optimization problem. This is in practice not achievable within finite computing time. An improved NMPC algorithm based on GA does not necessarily depend on a global or even local minimum. The optimizer provides a feasible decedent solution, instead of finding a global or local minimum. The feasible solution decreases the cost function instead of minimizing the cost function. Judicious selection of the termination criterions of GA is the key factor in reducing the computation burden in the design of the suboptimal NMPC algorithm. To this end, the following two strategies at the (*k*+1)-th step are proposed.

The control sequence output at the

*k*-th control interval in the genetic algorithm is always selected as one of the initial populations at the (*k*+1)-th control interval. Furthermore, some of the best individuals at the*k*-th control interval are also selected into the current initial population. Most of all, the elite-preservation strategy is adopted. Figure 1 shows one of the choices of the initial population per iteration. This strategy guarantees the quality of current population and the stability of the NMPC algorithm.

Stopping criterions of GA are the key factor of decreasing the computation burden. GA is used to compute the control sequence. The objective value

*J*(*k*+1) at the (*k*+1)-th control interval is computed and compared with*J*(*k*) that stored at the*k*-th control interval. If*J*(*k*+1) is smaller than*J*(*k*), that at the*k*-th control interval, then the control sequence u(*k*+1) is retained as a good feasible solution, and its first element*u*(*k*+1|*k*+1) is sent to the plant. Otherwise, if there does not exist a feasible value for u(*k*+1) to yield*J*(*k*+1) <*J*(*k*), then the best u(*k*+1) is chosen to decrease the objective function the most.

The traditional MPC approach requires the global solution of a nonlinear optimization problem. This is in practice not achievable within finite computing time. An improved NMPC algorithm based on GA does not necessarily depend on a global or even local minimum. The optimizer provides a feasible decedent solution, instead of finding a global or local minimum. The feasible solution decreases the cost function instead of minimizing the cost function. Judicious selection of the termination criterions of GA is the key factor in reducing the computation burden in the design of the suboptimal NMPC algorithm. To this end, the following two strategies at the (*k*+1)-th step are proposed.

With the above two strategies, the computational complexity of the control calculation is substantially reduced. Summarizing, our proposed improved NMPC algorithm performs the following iterative steps:

[Initialization]:

choose parameters *P*, *X* _{F }, *Q*, *R*, *Q*’ and model *x*(*k*+1) = *f*(*x*(*k*),*u*(*k*)); initialize the state and control variables at *k* = 0; compute and store *J*(0).

[modified Iteration]:

at the

*k*-th control interval, determine the control sequence u(*k*) using GA satisfies constraints (4) (5), terminal stability constraint (6) and*J*(*k*) <*J*(*k*-1). The first input*u*(*k*|*k*) is then sent into the plant.store

*J*(*k*) and set*k*=*k*+1;if there does not exist a feasible value for u(

*k*) to yield*J*(*k*) <*J*(*k*-1), then the best u(*k*) is chosen to decrease the objective function the most.[Termination]

The entire calculation is repeated at subsequent control interval *k*+1 and goes to Step 2.

Though the proposed method does not seek a globally or locally optimal solution within each iteration step, it may cause little performance degradation to the original GA due to its iterative nature, which is known to be capable of improving suboptimal solution step by step until reaching near-optimal performance at the final stage. Besides its near-optimal performance, the proposed algorithm possesses salient feature; it guarantees overall system stability and, most of all, leads to considerable reduction in the online computation burden. Finding a control sequence that satisfies a set of constraints is significantly easier than solving a global optimization problem. Here it is possible to obtain the suboptimal control sequence via GA for practical systems with very demanding computation load, that is, systems with a small sampling time or a large prediction horizon.

## 4. Stability of nonlinear model predictive control system

The closed-loop system controlled by the improved NMPC based on GA is proved to be stable.

Theorem 1: For a system expressed in (1) and satisfying the assumption A1-A4, the closed-loop system is stable under the improved NMPC framework.

Proof: Suppose there are an admissible control sequence u(*k*) and a state sequence x(*k*) that satisfy the input, state and terminal stability constraints at the sampling time *k*.

At the sampling time *k*, the performance index, which is related to u(*k*) and x(*k*), is described as

In the closed-loop system controlled by the improved NMPC, define the feasible input and state sequences for the successive state are *x* ^{+}=*f*(*x*, *u*(*k*|*k*)).

The resulting objective function of u(*k*+1) and x(*k*+1) at the (*k*+1)-th sampling time is

If the optimal solution is found, it follows that

From A4, the following inequality holds

Or using the improved NMPC algorithm, it follows that

Thus the sequence *P* time indices decreases. As such and given the fact that the cost function

Also, because the sequence

## 5. Simulation and experiment results

### 5.1. Simulation results to a continuous stirred tank reactor plant

#### 5.1.1. Model of continuous stirred tank reactor plant

Consider the highly nonlinear model of a chemical plant (continuous stirred tank reactor-CSTR). Assuming a constant liquid volume, the CSTR for an exothermic, irreversible reaction, A→B, is described by

where *C* _{A }is the concentration of *A* in the reactor, *T* is the reactor temperature and *T* _{c }is the temperature of the coolant stream. The parameters are listed in the Table 1.

Variables | Meaning | Value | Unit |

q | the inlet flow | 100 | l/min |

V | the reactor liquid volume | 100 | l |

CAf | the concentration of inlet flow | 1 | mol/l |

k _{0} | reaction frequency factor | 7.2*10^{10} | min-1 |

E/R | 8750 | K | |

E | activation energy | ||

R | gas constant | 8.3196*10^{3} | J/(mol K) |

Tf | the temperature of inlet flow | 350 | K |

∆H | the heat of reaction | -5*10^{4} | J/mol |

ρ | the density | 1000 | g/l |

CP | the specific heat capacity of the fluid | 0.239 | J/(g K) |

UA | 5*10^{4} | J/(min K) | |

U | the overall heat transfer coefficient | ||

A | the heat transfer area |

#### 5.1.2. Simulation results

The paper present CSTR simulated examples to confirm the main ideas of the paper. The nominal conditions, *C* _{A }= 0.5mol/l, *T* = 350K, *T* _{c }= 300K, correspond to an unstable operating point. The manipulated input and controlled output are the coolant temperature (*T* _{c }) and reactor temperature (*T*). And the following state and input constraints must be enforced:

The simulation platform is MATLAB and simulation time is 120 sampling time. The sampling time is *Ts* = 0.05s, mutation probability is *P* _{c }=0.1, population size is 100, maximum generation is 100, and the fitness value is 1/(*J*+1).

method | settling time, min | percent overshoot,% |

NMPC algorithm based on GA | 0.75 | 1 |

Suboptimal NMPC algorithm based on GA | 2.5 | 0 |

For comparison, the same simulation setup is used to test both the conventional NMPC algorithm based GA and the suboptimal NMPC algorithm. The resulting control values are depicted in Fig 2. Table 2 compares the performance of the two algorithms using the metrics settling time and percent of overshoot. The conventional NMPC algorithm has a faster transient phase and a smaller percentage of overshoots.

When the population sizes or the maximum generation is relatively large, the time consumption of the two methods is compared in Fig 3.

From Fig 2 and Table 2, it is apparent that the control performance of the two methods is almost same. But from Fig 3, it is evident that the suboptimal NMPC algorithm based on GA has a considerably reduced demand on computational complexity.

### 5.2. Simulation results to a coupled-tank system

#### 5.2.1. Model of coupled-tank system

The apparatus [15], see Fig.4, consists of two tanks *T* _{1} and *T* _{2}, a reservoir, a baffle valve *V* _{1} and an outlet valve *V* _{2}. *T* _{1} has an inlet commanded through a variable pump based on *PMW* and *T* _{2} has an outlet that can be adjusted through a manually controlled valve only. The outlets communicate to a reservoir from which the pumps extract the water to deliver it to the tank. The two tanks are connected through the baffle valve, which again can only be adjusted manually. The objective of the control problem is to adjust the inlet flow so as to maintain the water level of the second tank close to a desired setpoint.

The water levels *h* _{1} and *h* _{2}, which are translated through the pressure transducer into a *DC* voltage ranging from 0*V* to 5*V*, are sent to *PC* port via *A/D* transition. The tank pump control, which is computed by the controller in *PC* with the information of the water level *h* _{1} and *h* _{2}, is a current level in the range 4*mA* to 20*mA*, where these correspond to the pump not operating at all, and full power respectively.

The dynamics of the system are modeled by the state-space model equations:

where the flows obey Bernoulli’s equation^{[16]}, i.e.

*z*.

The output equation for the system is

The cross-section areas, i.e. *A* and *S,* are determined from the diameter of the tanks and pipes. The flow coefficients, *μ* _{1} and *μ* _{2, }have experimentally (from steady-state measurements) been determined. Table 3 is the meanings and values of all the parameters in Eqn.15

Signal | Physics Meaning | Value |

A | Cross-section area of tank | 6.3585×10^{-3}m^{2} |

S | Cross-section area of pipe | 6.3585×10^{-5}m^{2} |

g | acceleration of gravity | 9.806m/s^{2} |

μ _{1} | flow coefficient 1 | 0.3343 |

μ _{2} | flow coefficient 2 | 0.2751 |

Several constraints have to be considered. Limited pump capacity implies that values of *cm* ^{3 }*/s*. The limits for the two tank levels, *h* _{1} and *h* _{2}, are from 0 to 50*cm*.

#### 5.2.2. Simulation results

The goal of the couple-tank system is to control the level of Tank 2 to setpoint. The initial levels of the two tanks, *h* _{1} *, h* _{2}, are 0*cm*. The objectives and limits of the tank system: Input constraint is 0≤*u*≤100%; State objectives are 0≤*h* _{1} *, h* _{2}≤0.5*m,* and the setpoint of Tank 2 is 0.1m.

The simulation platform is MATLAB and simulation time is 80 sample time. In NMPC, select prediction horizon *P=*10, weighting parameters*R=*1, sample time *Ts=5s*, mutation probability *Pc=*0.1, population size is 200, maximum generation is 100, the fitness value is 1/(*J+*1).

For the purpose of comparison the same simulation is carried out with the conventional NMPC algorithm based GA and the fast NMPC algorithm. The result is shown in Fig 5. The performance indexes of the two algorithms are shown in Table 4.

method | Settling time, s | percent overshoot, % |

NMPC algorithm based on GA | 70 | 3 |

Fast NMPC algorithm based on GA | 100 | 7 |

When the population sizes or the maximum generation is relatively larger, the time consumptions of the two method is shown in Fig 6.

From Fig 5 and Table 4, it is apparent that the control performance of the two methods is almost same. But from Fig 6, the computation demand reduces significant when the fast NMPC algorithm based on GA is brought into the system.

#### 5.2.3. Experiment results

The objectives and limits of the system: Input constraint is 0*≤u≤*100%; State objectives are *0≤h* _{1} *,h* _{2} *≤0.5m,* and the setpoint of Tank 2 is 0.1m. Select prediction horizon *P=*10, weighting parameters*R=*1, sample time Ts=5s, mutation probability *Pc=*0.1, population size is 200, maximum generation is 100.

The tank apparatus is controlled with the NMPC algorithm based on conventional GA, the experimental curve is shown in Fig 7 and performance index is shown in Table 5.

Time, s | 0-400 | 401-800 | 801-1200 |

Setpoint, m | 0.1 | 0.15 | 0.07 |

Settling time, s | 159 | 190 | 210 |

Percent overshoot | None | None | None |

The same experiment is carried out with fast NMPC algorithm based on GA. The result is shown in Fig 8 and performance index is shown in Table 5.

Time, s | 0-400 | 401-800 | 801-1200 |

Setpoint, m | 0.1 | 0.15 | 0.07 |

Settling time, s | 190 | 190 | 260 |

Percent overshoot | None | None | None |

## 6. Conclusions

In this paper an improved NMPC algorithm based on GA has been proposed. The aim is to reduce the computational burden without much deterioration to the control performance. Compared with traditional NMPC controller, our approach has much lower computational burden, which makes it practical to operate in systems with a small sampling time or a large prediction horizon.

The proposed approach has been tested in CSTR and a real-time tank system. Both computer simulations and experimental testing confirm that the suboptimal NMPC based on GA resulted in a controller with less computation time.

## Acknowledgments

This work is supported by National Natural Science Foundation of China (Youth Foundation, No. 61004082) and Special Foundation for Ph. D. of Hefei University of Technology (No. 2010HGBZ0616).