A novel nonlinear model predictive control design based on a hybrid particle swarm optimization-sequential quadratic programming algorithm: application to an evaporator system

Abstract

This paper proposes a hybrid algorithm by combining particle swarm optimization (PSO) with sequential quadratic programming (SQP) to handle the optimization part of the nonlinear model predictive control (NMPC). In the proposed method, a nonlinear model of the plant is directly applied to predict the future behaviour of the system for a certain horizon. In each sampling time a constrained nonlinear optimization (CNO) problem has to be solved in order to find the optimizing input sequence with predefined length. These solutions minimize the errors between the predicted states and the desired ones of the system. Only the first optimizing input is applied to the system and the rest of the sequence is discarded. The whole prediction and optimization must be repeated for the next upcoming steps. The evolutionary and heuristic PSO algorithm in cooperation with the powerful local minimizing SQP algorithm can quickly find these optimizing solutions and also can satisfy the constrains. The evaporation process due to its nonlinear dynamics and the existence of disturbances affecting the evaporator is considered to evaluate the performance of the proposed method. The simulation results show that the controller designed via the proposed approach can make the system track the reference commands and satisfy the restrictions properly. Moreover, based on the simulation results it can be seen that this approach is more efficient in comparison with the NMPC method, in which only the SQP algorithm has been utilized for optimization, and also the linear model predictive control method.

Keywords

Linear model predictive control nonlinear model predictive control constrained nonlinear optimization particle swarm optimization sequential quadratic programming evaporator

Introduction

Model predictive control (MPC) is one of the few control techniques that is able to consider constraints, on both states and inputs of the system, in the design of the control law for linear, nonlinear or uncertain systems (Piche et al., 2000). MPC employs an explicit model to predict and optimize the future behaviour of a plant. At each sampling time, an open-loop optimal control sequence is obtained through solving an optimization problem, and then the first element of this sequence is applied to the plant, while the remaining elements are discarded. This optimization procedure is repeated in every step.

Although most industrial processes are nonlinear, applications of MPC are commonly based on linear dynamic models. In the standard MPC, usually a nonlinear model of the system is linearized around the operating point. Using linear model predictive control (LMPC) is usually recommended when the process operates at a single operating point and the primary use of the controller is to reject small disturbances. However, in practice it is often required to track a changing set point, and to minimize offset as much as possible in the presence of significant plant–model mismatches and disturbances. To address reference tracking problems, Alvardo et al. (2011) and Limon et al. (2008, 2010) have proposed to augment the linear model with a static model for the disturbance and to use estimators to estimate the states and disturbances. While these methods in the case of model mismatch are able to yield the offset free tracking for nonlinear plant, the transient behaviour and the control performance are not satisfactory, especially in cases in which there are changing disturbances of high magnitudes.

In cases that the exact model of the systems is unknown, a combination of MPC with intelligent methods such as neural network methods or fuzzy logic algorithms is a noteworthy option. For example, in Dalamagkidis et al. (2011) a combination of neural networks and nonlinear model predictive control (NMPC) has been proposed to control an unmanned helicopter. In Ding (2011) and Ding and Ping (2013), infinite MPC has been proposed to control systems modelled by the fuzzy logic method. In this method for each sampling interval, a linear Takagi–Sugeno model has been obtained based on the system behaviour and then it has been used in prediction control design. In Aswani et al. (2013), using two models, MPC has been applied to systems with uncertainties. The first model was an approximation of the model together with limited uncertainties and the second one updated each time based on the statistical method. Despite the advantages of such hybrid methods, they generally depend on input–output sample data and numerical solving of the optimization problem may be affected due to approximation in the predictive model. When a nonlinear plant model is available, it makes sense to use it for prediction in MPC. It decreases the prediction errors, resulting in better control decisions.

Morari and Maeder (2012) and Yang et al. (2012) proposed to augment the nonlinear model of the plant with a static model for the disturbance and, like the LMPC method, estimate the plant state and disturbance values. These methods are more efficient compared with the LMPC approaches in addressing nonlinear plants. However, these methods are not efficient when the disturbances affecting the process are changing frequently during the process. Furthermore, the calculation burden is very cumbersome because the use of the nonlinear model to predict and estimate is time consuming. This defect in on-line constrained nonlinear optimization (CNO), as the key part of MPC, is irreparable. This difficulty is cumbersome, especially in the case of large dimension and highly nonlinear processes.

Therefore, improvement in NMPC is of great importance, especially for general models of nonlinear systems and in the presence of disturbance. There are several methods to finding the solution to a CNO problem, for example gradient descents, iterative dynamic programming, sequential quadratic programming (SQP), the interior point algorithm and so on (Boggs and Tolle, 2000; Schweitzer, 1971; Sorensen, 2006). In order to decrease the calculation burden, in Almaro et al. (2014) and Yang and Biegler (2013), based on gradient and sensitivity of cost function, effective methods have been proposed. Although these methods are successful in NMPC, they are severely dependent on a suitable initial guess for a solution, especially in cases when the cost functions to be minimized contain different local minima. In NMPC, due to the use of the nonlinear model for prediction, the cost function to be minimized is usually non-convex and therefore it may have several local minima. To cope with this disadvantage, the paper proposes that an evolutionary algorithm, which is able to explore the admissible region of potential solutions, is used to find the best initial guess and then introduce it to the nonlinear optimization algorithms.

The proposed method is a combination of particle swarm optimization (PSO) as the evolutionary algorithm and SQP as the CNO approach to deal with the nonlinear programming (NLP) due to their accuracy and fast converging. SQP is a NLP method that starts from a single searching point and finds a solution using the gradient information (Stoer and Bulirsch, 2002). The SQP algorithm is described in Appendix A. In addition, the PSO is a heuristic optimization algorithm based on simulating the movement and flocking of birds (Kennedy and Eberhart, 1995). In a solution space, a randomly chosen flock of solutions, which are commonly called particles, are applied to initialize the search process. The position of each particle iteratively changes by its velocity. Each particle’s velocity is adjusted based on its performance and that of the whole swarm in finding the best position. Such iterative changes in the positions of the particles make the swarm move towards better solution regions. Summarized details of the PSO algorithm are included in Appendix B.

The PSO part of the hybrid algorithm searches the solution area, which is determined by the constraints of the plant. It starts searching by randomly initializing the position of a group of particles in that area. Each particle position is a potential solution, and the PSO algorithm tries to find the best one. After PSO finds a position close to the global best position in the solution space, the SQP algorithm then takes over and searches around that position in order to converge to the global optimum. That is because PSO usually can barely converge to the optimum after getting close to the best position. In this way, the hybrid algorithm may find an optimum more quickly and more accurately. Therefore, the paper goes beyond the previous study on SQP-based NMPC by improving the optimization using an evolutionary algorithm, and obtains much better optimization as an important part of NMPC.

An evaporator system as a test-bed is considered to evaluate the proposed approach. It is a nonlinear multi-input multi-output system. Compared with NMPC, in which a simple SQP is utilized for optimization, our method can considerably enhance the results for nonlinear systems. Simulation results obtained by applying the proposed scheme to the evaporator system show that the controller designed with the proposed approach can make the plant track the piecewise constant references despite the influence of the frequently changing disturbances affecting the plant.

The rest of the paper is organized as follows: an overview of NMPC is presented in the following section. The third section is dedicated to the proposed approach for NMPC problems. A case study to evaluate the proposed method is given in the fourth section. Simulation results are presented in the fifth section and, finally, in the sixth section concluding remarks are presented.

Nonlinear model predictive control

Problem statement

Consider the nonlinear system to be controlled is given by the following dynamic equation:

\begin{matrix} x_{p} (k + 1) = f_{p} (x_{p} (k), u (k), d (k)) \\ y_{p} (k) = g_{p} (x_{p} (k)) \end{matrix}

(1)

where $x_{p} \in x \subset R^{n}$ is the state variable vector, $u \in U \subset R^{m}$ is the input vector, $y_{p} \in Y \subset R^{ℓ}$ is the output vector and $d \in R^{q}$ is the disturbance vector of the process. X, U and Y are admissible sets due to physical constraints and design requirements. In addition, f_p and g_p are continuously differentiable nonlinear vector functions.

Consider a nominal model of the process as follows:

\begin{matrix} x (k + 1) = f (x (k), u (k)) \\ y (k) = g (x (k)) \end{matrix}

(2)

where $x \in X$ , $u \in U$ and $y \in Y$ are state, input and output variables of the model, respectively. f and g are continuously differentiable nonlinear vector functions.

We take $r \in Y \subset R^{ℓ}$ as the reference vector signal and assume it to be piecewise constant. It is desired to find a controller to make the process follow the reference signals as closely as possible while the constraints have to be satisfied at the same time.

Controller design

To achieve the control objective, LMPC as a popular predictive control method has been widely used in industry. In this method, it is convenient to extract a linear model around the operating point of the nonlinear plant. When disturbances affecting the process are constant, according to Limon et al. (2008, 2010) and Alvardo et al. (2011) the linear model is augmented with a static model for disturbances as shown in (3), and then it is applied to prediction phase of MPC:

\begin{matrix} \hat{x} (k + 1) = A \hat{x} (k) + B \hat{u} (k) + Γ d (k) \\ d (k + 1) = d (k) \\ \hat{y} (k) = C \hat{x} (k) + Ψ d (k) \end{matrix}

(3)

where $\hat{x} \in \hat{X} \subset R^{n}$ , $\hat{u} \in \hat{U} \subset R^{m}$ and $\hat{y} \in \hat{Y} \subset R^{ℓ}$ are respectively the state, input and output variables of the linearized model. These variables are deviations from the operating point. $\hat{X}$ , $\hat{U}$ and $\hat{Y}$ are the admissible sets for the deviations. A, B, Γ, C and Ψ are matrices of proper dimensions. Based on the linear augmented model, it is required that the $\hat{x} (k)$ and d(k) be estimated during each sampling time by an estimator.

In the case that the disturbances are varying frequently around their nominal value and moreover the varying is of a meaningful amount, LMPC is not able to reach the control objective. In this case, NMPC can effectively improve the system response since it uses the nonlinear model of the plant for prediction.

For the nonlinear model (2), in the current sample k, assuming the states variables of the process are accessible, consider the general optimal control problem as follows:

\begin{matrix} \min_{\underline{u}} & {(x_{N_{p}} - x_{e})}^{^{T}} ( \end{matrix} x_{N_{p}} - x_{e}) + \sum_{j = 1}^{N_{p} - 1} [{(x_{j} - x_{e})}^{^{T}} Q (x_{j} - x_{e}) + {(u_{j - 1} - u_{e})}^{^{T}} R (u_{j - 1} - u_{e})]

(4)

\begin{matrix} \begin{matrix} s . t . & x_{e} = f (x_{e}, u_{e}), r = g (x_{e}) \end{matrix} \\ x_{0} = x_{p} (k) \\ \begin{matrix} x_{i + 1} = f (x_{i}, u_{i}); & i = 0, 1, . . ., N_{p} - 1 \end{matrix} \\ \begin{matrix} y_{i} = g (x_{i}) \\ \begin{matrix} \begin{matrix} x_{i + 1} \in X & , x_{e} \in X \end{matrix} \\ u \begin{matrix} _{i} \in U & , u_{e} \in U \end{matrix} \end{matrix} \end{matrix} \end{matrix}

where $\underline{u} = (u_{0}, u_{1}, \dots, u_{N_{p} - 1})$ . x_j and u_j are the predicted state and input vectors j steps in the future, respectively, and (x_e , u_e ) is the equilibrium point of the nonlinear model corresponding to the reference r. In addition, Q and R are positive definite weighting matrices of proper sizes to penalize the deviations.

Considering the optimizer of (4) as $u^{*} = (u_{0}^{*}, u_{1}^{*}, \dots, u_{N_{p} - 1}^{*})$ , we apply the first element of u ^* as input to the system:

u (k) = u_{0}^{*}

(5)

To facilitate the calculation of (4), it is common to assume that the control decisions are the same after a horizon, the so-called the control horizon, that is, $\underline{u} = (u_{0}, u_{1}, \dots, u_{N_{c}}, u_{N_{c}}, \dots, u_{N_{c}})$ , where N_c is the control horizon.

The optimization problem (4) must be solved during every sampling time. Then the solution is applied to the process in a receding horizon manner.

Since the nonlinear model of the system is used for prediction, the optimization problem (4) becomes a CNO problem that is not convex and often has different local solutions. So, unlike in LMPC, the optimization cannot be solved by quadratic programming (QP) algorithms, as QP is the conventional approach in LMPC.

To find the solutions of (4) a nonlinear program is required. In order to reach this goal, there are several CNO methods that are able to handle non-convex optimization problems.

The proposed method

As mentioned in the last paragraphs of the second section, the CNO problem (4) requires that the NLP algorithm be solved numerically. In order to find the solution of (4), the SQP approach as a constrained nonlinear minimizing algorithm can be used because of its fast convergence and high accuracy around local minima. In this paper we consider a class of nonlinear autonomous discrete-time systems of the following equation form:

\begin{matrix} x_{p} (k + 1) = F_{p} (x_{p} (k), d (k)) + G_{p} (x_{p} (k), d (k)) u (k) \\ y_{p} (k) = g_{p} (x_{p} (k)) \end{matrix}

(6)

where x_p , y_p , u, d and g_p are the same as defined for (1). Also F_p , G_p are continuously differentiable nonlinear functions.

Therefore, the model used for predicting the future behaviours of the plant is as below:

\begin{matrix} x (k + 1) = F (x (k)) + G (x (k)) u (k) \\ y (k) = g (x (k)) \end{matrix}

(7)

where x, y and g are the same as defined for (2). F and G are continuously differentiable nonlinear functions.

Based on Equations (4) and (7), the optimization problem becomes

{\begin{matrix} \min_{\underline{u}} & J = [F (x_{N} {_{p}}_{-}_{1}) + G (x_{N} {_{p}}_{-}_{1}) u_{N} {_{p}}_{-}_{1} - x_{e}] \end{matrix}}^{T} Q [F (x_{N} {_{p}}_{-}_{1}) + G (x_{N} {_{p}}_{-}_{1}) u_{N} {_{p}}_{-}_{1} - x_{e}] + \sum_{j = 1}^{N_{p} - 1} {[{(F (x_{N} {_{p}}_{-}_{1}) + G (x_{N} {_{p}}_{-}_{1}) u_{N} {_{p}}_{-}_{1} - x_{e})}^{T} Q (F (x_{N} {_{p}}_{-}_{1}) + G (x_{N} {_{p}}_{-}_{1}) u_{N} {_{p}}_{-}_{1} - x_{e})]

(8)

\begin{matrix} \begin{matrix} s . t . & x_{e} = F (x_{e}) + G (x_{e}) u_{e}, r = g (x_{e}) \end{matrix} \\ x_{0} = x_{p} (k) \\ \begin{matrix} x_{i + 1} = F (x_{i}) + G (x_{i}) u_{i}; & i = 0, 1, . . ., N_{p} - 1 \end{matrix} \\ \begin{matrix} y_{i} = g (x_{i}) \\ \begin{matrix} \begin{matrix} x_{i + 1} \in X & , x_{e} \in X \end{matrix} \\ u \begin{matrix} _{i} \in U & , u_{e} \in U \end{matrix} \end{matrix} \end{matrix} \end{matrix}

Calculating Equation (8), we obtain a continuously differentiable nonlinear function of the input vector as follows:

\begin{matrix} J = J (\underline{u}) \\ \begin{matrix} s . t . & h (\underline{u}) = 0 \end{matrix} \\ Ω (\underline{u}) \leq 0 \end{matrix}

(9)

where $h (\underline{u})$ and $Ω (\underline{u})$ are derived from conditions of Equation (8). Therefore, the Lagrangian function, according to equation (A.2) in Appendix A, will be

L (\underline{u}, λ, μ) = J (\underline{u}) + λ^{T} h (\underline{u}) + μ^{T} Ω (\underline{u})

(10)

where the vectors $λ \in R^{m_{h}}$ and $μ \in R^{ℓ_{Ω}}$ are referred to as Lagrangian multipliers. Then, we can apply following sub-optimization problem to find an optimal solution:

\begin{matrix} {\begin{matrix} min_{\underline{Δ u}} & \nabla L ({\underline{u}}^{k}, λ^{k}, μ^{k}) \end{matrix}}^{T} \underline{Δ u} + \frac{1}{2} {\underline{Δ u}}^{T} H ({\underline{u}}^{k}, λ^{k}, μ^{k}) \underline{Δ u} \\ \begin{matrix} s . t . & h ({\underline{u}}^{k}) + \nabla h {({\underline{u}}^{k})}^{T} \underline{Δ u} = 0 \end{matrix} \\ Ω ({\underline{u}}^{k}) + \nabla Ω {({\underline{u}}^{k})}^{T} \underline{Δ u} \leq 0 \end{matrix}

(11)

where $\underline{Δ u} = \underline{u} - {\underline{u}}^{k}$ is the input growth vector and u^k corresponds to iterate $k . λ^{k}$ and $μ^{k}$ are the associated Lagrangian multipliers of ${\underline{u}}^{k}$ .

To find the optimum solution of problem (9) using SQP, it is required to initialize $\underline{u}$ with a proper point. The important defect of SQP (the same as for many other nonlinear optimizing algorithms, such as the interior point algorithm, iterative dynamic programming, etc.) is that it cannot escape local optima to reach the global optimum. Therefore, the result is severely dependent on a proper initial guess for a solution. Meeting this requirement is critical to avoid getting stuck in local optima. On the other hand, it is not easy to find such a proper initial guess. To cope with this problem, we propose the use of PSO as an evolutionary algorithm. The particles of swarm can explore a predefined search space determined by the conditions of Equation (9) and find the best initial position for $\underline{u}$ around the global optimum. Although there is no guarantee that PSO can converge to the global optimum, this evolutionary algorithm can pass more local optima based on the number of particles and on the number of iterations. This ability of PSO helps escape many local optima and reach a point that is better than those many local optima in minimizing the cost function. This ability can even result in a point in the valley where the global optimum is located.

In order to show how PSO can find a proper initial point for SQP, we consider the optimum problem for two typical cost functions with different local optima densities as follows.

In the case that the density of local optima in the admissible range is not high, PSO can easily converge to a proper initial point in the valley where the global optimum is located. Then SQP starts from this point and converges to the global optimum. Figure 1 shows a typical cost function $v = ξ (σ)$ that has few local optima. $σ_{L}$ and $σ_{U}$ are the lower and the upper bounds of the admissible set of $σ$ , respectively. In Figure 1 there are three valleys and $σ_{1}$ , $σ_{2}$ and $σ_{3}$ are the local minima of $ξ$ . $σ_{2}$ is also the global optimum. v ₁, v ₂ and v ₃ are values of $ξ$ in $σ_{1}$ , $σ_{2}$ and $σ_{3}$ , respectively. In Figure 1 we can see that PSO skipped $σ_{1}$ and $σ_{3}$ and converged to $σ_{init}$ . Starting from $σ_{init}$ , SQP will converge to global optimum ( $σ_{2}$ ), since these two points are in the same valley.

In the case that the density of local optima is high, PSO might converge to a point that is not in the global optimum’s valley. Therefore, starting from this point SQP will converge to an optimum that is not global. However, depending on the number of particles and the number of iterations in PSO, we can reach a local optimum that has acceptable performance in minimizing the cost function. For an optimization closer to that of the global optimum, we need to increase the number of particles and iterations. However, it must be noted that the increase in these numbers will take more time for calculation. This might even exceed the sample time, resulting in impractical control. For instance in Figure 2, the cost function $v = ξ (σ)$ has too many local optima. $σ_{1}$ , $σ_{2}$ ,…, $σ_{11}$ are the local minima and v ₁,v ₂,…,v ₁₁ are values of $ξ$ in local minima. In Figure 2, we can see that PSO converged to $σ_{init}$ in the $σ_{10}$ valley, then SQP will converge to $σ_{10}$ . It is obvious that $σ_{6}$ is the global optimum and $σ_{10}$ is the local one. However, y ₆ and y ₁₀ show that $σ_{10}$ has the closest performance to that of $σ_{6}$ in minimizing the cost function.

Therefore PSO can help to find the global optimum by delivering a much better point as an initial guess to SQP.

However, it must be noted that the time taken by the optimization process between two consecutive samples is required to be less than the sample time as much as possible. This is necessary in on-line optimization. Therefore, NMPC can be used with any learning method that is fast enough to converge to a proper point during the sample time. The selection of PSO as an evolutionally algorithm is motivated by the fact that it can converge to a proper initial point, quickly. Furthermore, unlike the other heuristic algorithms, such as the genetic algorithm (GA) and Differential Evolution (DE), it does not require many operators (Vesterstrom and Thomsen, 2004).

This combination of PSO and SQP constructs a hybrid optimizing algorithm with quick and accurate performance. This advantage of the proposed method can be remarkably useful in on-line optimization, which is a key part of NMPC. A block diagram of the proposal approach is shown in Figure 3.

The procedure for designing the NMPC controller using the proposed hybrid method can be summarized as follows.

Step1: Using the nonlinear model (2) and x_p (k) measured from plant, predict the future state variables over the prediction horizon N_p .

Step2: Place the predicted state variables into (4) to form a nonlinear cost function that is going to be minimized.

Step 3: Based on the constraints on the state and input variables of the plant, determine the admissible set of the solutions.

Step 4: Run the PSO algorithm to find the best solution of the resulting CNO problem. In PSO, the positions and the velocities of a group of particles must be initialized randomly within the range determined in the previous step. (Each particle represents a potential solution.)

Step 5: After PSO has terminated, apply the best solution obtained as the initial guess for the SQP algorithm.

Step 6: Run the SQP algorithm to search around the initial guess to find the optimum solution.

Step 7: When SQP has terminated, apply the first element of the optimal solution (optimal control horizon sequence) to the plant and discard the others and go back to step 1 for the next sample (k+1).

Case study: evaporator

In this section an industrial process, namely the evaporator, is introduced as a case study to implement the proposed approach. It possesses a highly nonlinear dynamic and can be affected by frequently varying disturbances during the process time. We assume that (i) the dynamics of the plant are nonlinear, (ii) the state variables can be measured and (iii) the states and inputs of the plant are subject to hard constraints. Figure 4 shows the diagram of the evaporator system, where F ₂, P ₁₀₀ and F ₂₀₀ are the manipulated variables. F ₁, X ₁, T ₁ and T ₂₀₀ are the input disturbances. L ₂, X ₂ and P ₂ are the state variables of the plant and are measurable. X ₂ and P ₂ are also the controlled variables. A description of the system variables and their values are given in Table 1.

Table 1.

Variables of the plant (Cao, 2005; Newell and Lee, 1989).

Variables	Descriptions	Value	Unit
F ₁	Feed flow rate	10	kg/min
F ₂	Product flow rate	2	kg/min
F ₃	Circulating flow rate	50	kg/min
F ₄	Vapour flow rate	8	kg/min
F ₅	Condensate flow rate	8	kg/min
X ₁	Feed composition	5	%
X ₂	Product composition	25	%
T ₁	Feed temperature	40	^°C
T ₂	Product temperature	84.6	^°C
T ₃	Vapour temperature	80.6	^°C
L ₂	Separator level	1	m
P ₂	Operating pressure	50.5	kPa
F ₁₀₀	Steam flow rate	9.3	kg/min
T ₁₀₀	Steam temperature	119.9	^°C
P ₁₀₀	Steam pressure	194.7	kPa
Q ₁₀₀	Heat duty	339	kW
F ₂₀₀	Cooling water flow rate	208	kg/min
T ₂₀₀	Inlet cooling water temperature	25	^°C
T ₂₀₁	Outlet cooling water temperature	46.1	^°C
Q ₂₀₀	Condenser duty	307.9	kW

The dynamic equations of the plant are as follows (Newell and Lee, 1989):

\begin{matrix} \frac{{dL}_{2}}{dt} = \frac{F_{1} - F_{4} - F_{2}}{20} \\ \frac{{dX}_{2}}{dt} = \frac{F_{1} X_{1} - F_{2} X_{2}}{20} \\ \frac{{dP}_{2}}{dt} = \frac{F_{4} - F_{5}}{4} \end{matrix}

(12)

where the equations to calculate F ₄ and F ₅ are given in Appendix C. We can see from Appendix C that F ₄ and F ₅ are nonlinear functions of state, input and disturbance variables. Therefore, F ₄, F ₅, F ₁ X ₁ and F ₂ X ₂ in Equation (12) are implying that the evaporator is a nonlinear system.

The control goal is to design a controller based on NMPC that can make controlled variables X ₂ and P ₂ track the stepwise, changing references s2 and s3, while disturbances as shown in Figure 5 are varying within ±20% of their nominal values. These nominal values are available in Table 1. This objective must be achieved under the conditions that the hard constraints (i) 0 ≤ F ₂ ≤ 4, (ii) 0 ≤ P ₁₀₀ ≤ 400, (iii) 0 ≤ F ₂₀₀ ≤ 400 must be satisfied and L ₂ must be kept as close as possible to its nominal value (s1 = 1 m).

Simulation results

To validate the NMPC controller designed based on the proposed hybrid method, the case study problem is simulated using MATLAB software. The sampling time T_s is assigned to 1 min, the predictive horizon N_p = 5 and the control horizon N_c = 3. In the cost function, the weighting matrices Q = 100 × I and R = I were considered. In the PSO algorithm the number of iterations and particles were both assigned 50, and the parameters w _max = 0.9, w _min = 0.4, c ₁ = c ₂ = 2. For the evaporator, the state variables L ₂, X ₂ and P ₂ initialized according to Table 1 at 1 (m), 25% and 50.5 (kPa), respectively.

Figure 6 shows a comparison between NMPC based on a simple SQP without utilizing PSO, proposed NMPC and LMPC. The simulation results show that despite the fact that reference inputs are changing stepwise, and the disturbances are changing frequently, the proposed NMPC controller is able to make the plant track the reference commands. LMPC, on the contrary, despite leading L ₂ towards impractical values, cannot perform as well as proposed method in tracking. For SQP-based NMPC, in each sampling interval, the optimization process is initialized with guesses of optimal inputs. These guesses are the average points chosen from a 10-time randomization in the admissible ranges. Figure 6 shows that SQP is not sufficient at all to handle the optimization without starting from proper points. It is shown in Figure 7 that the constraints on manipulated variables are satisfied by the hybrid approach. The achievement for the proposed controller is due to the use of a nonlinear model to predict the future behaviours of the plant and, of course, the aid of the hybrid PSO-SQP algorithm to solve the CNO problem quickly and accurately during each sampling time. Figure 8 shows the time taken by the proposed hybrid algorithm to converge to the optimal solution per step. It signifies the advantage of employing the proposed approach in NMPC in the case of industrial processes, which are usually not fast.

In order to show that utilizing the PSO algorithm can efficiently improve the performance compared with other intelligent optimization algorithms, Figure 9 shows a comparison between NMPC based on SQP-DE, SQP-GA and the proposed method. In simulations the population and iteration numbers in DE and GA were the same as in PSO. For DE we assigned the crossover constant 0.5 and the deferential variation factor 1. In GA we assigned the crossover rate 0.8 and the mutation rate 0.01. Further information on DE and GA and their parameters can be found in Cruz et al. (2003) and Schmitt (2004). As shown in Figure 9, the proposed method makes the outputs of the evaporator follow the reference inputs more precisely. Furthermore, Figure 10 shows that NMPC with PSO is faster than NMPC with DE or GA. Therefore, we can conclude from Figures 9 and 10 that NMPC is more effective with PSO.

Conclusion

In MPC it is very profitable to use the nonlinear model instead of the linearized model in nonlinear systems. However, it is not easy to solve the constrained optimization problem resulting in NMPC, because it usually contains different local optima. Moreover, the accuracy and fast converging to an optimum solution are very important in on-line optimization. Consequently, this drawback may result in non-optimized performance of the controller. This paper proposed a hybrid method combining PSO with SQP to cope with that shortage. The cooperation between PSO and SQP in NMPC yields a quick and accurate algorithm to find the solution to the optimization problem. Simulation results on a nonlinear evaporator system showed that the proposed method can be significantly useful to design controllers based on NMPC.

Footnotes

Appendix A: sequential quadratic programming

Consider a nonlinear optimization problem of the form

(A.1)

\min_{z} Φ (z)

\begin{matrix} s . t . h (z) = 0 \\ Ω (z) \leq 0 \end{matrix}

where $Φ : R^{n_{z}} \to R$ is the cost function to be minimized. $z \in F \subset R^{n_{z}}$ is the solution to the minimizing problem (A.1). $F = {z \in R^{n_{z}} | h (z) = 0, Ω (z) \leq 0}$ is a feasible set for z. $h : R^{n_{z}} \to R^{m_{h}}$ and $Ω : R^{n_{z}} \to R^{ℓ_{Ω}}$ contain the equality and inequality constraints.

This method is based on iterative formulation and on the solution of QP sub-problems (Nocedal and Wright, 2006). To construct the sub-problems, we first present the following definition.

Definition (Nocedal and Wright, 2006). Lagrangian functional associated with NLP:

The functional $L : R^{n_{z} \times m_{h} \times ℓ_{Ω}} \to R$ defined by means of

(A.2)

L (z, λ, μ) = Φ (z) + λ^{T} h (z) + μ^{T} Ω (z)

is called the Lagrangian functional of the NLP (A.1). The vectors $λ \in R^{m_{h}}$ and $μ \in R^{ℓ_{Ω}}$ are referred to as Lagrangian multipliers.

With this definition we can say now that the sub-problems are obtained by linearizing the constraints and approximating the Lagrangian function quadratically as follows:

(A.3)

\min_{d_{z}} \nabla L {(z^{k}, λ^{k}, μ^{k})}^{T} d_{z} + \frac{1}{2} d_{z}^{T} H (z^{k}, λ^{k}, μ^{k}) d_{z}

\begin{matrix} s . t . h (z^{k}) + \nabla h {(z^{k})}^{T} d_{z} = 0 \\ Ω (z^{k}) + \nabla Ω {(z^{k})}^{T} d_{z} \leq 0 \end{matrix}

where d_z = z–z^k and z^k is called iterate k, and λ^k, µ^k are its associated Lagrangian multipliers. The vector $\nabla L (z^{k}, λ^{k}, μ^{k})$ is the gradient of Lagrangian function at $(z^{k}, λ^{k}, μ^{k})$ , and the matrix $H (z^{k}, λ^{k}, μ^{k})$ is a positive definite approximation of the Hessian matrix of the Lagrangian function, which can be updated by any of the quasi-Newton methods (Shen et al., 2010).

This sub-problem is a QP sub-problem; its solution is used to form the next iterate:

(A.4)

z^{k + 1} = z^{k} + α_{k} d_{z}

The step length parameter α_k is determined by an appropriate line search procedure so that a sufficient decrease in a merit function is obtained (Nocedal and Wright, 2006). Although this optimizing method is less time consuming than the population-based search algorithms, it is highly dependent on the initial estimate of the solution.

Appendix B: particle swarm optimization

At each iteration, the velocity vector, which changes each particle’s position, is calculated as follows:

(B.1)

V_{i}^{t + 1} = {wV}_{i}^{t} + c_{1} r_{1} ({Pbest}_{i}^{t} - P_{i}^{t}) + c_{2} r_{2} ({Gbest}^{t} - P_{i}^{t})

where P_i ^t is the position of the particle i in the tth iteration, Pbest_i ^t is the best position of this particle up to iterate t (i.e. personal best), Gbest^t is the best position up to iterate t among all particles (i.e. global best), w is the inertia weight and c ₁ and c ₂ are acceleration coefficients. r ₁ and r ₂ are two random numbers in the range [0, 1]. The position of each particle can be updated for next iteration as follows:

(B.2)

P_{i}^{t + 1} = V_{i}^{t + 1} + P_{i}^{t}

The performance of PSO is highly reliant on the mechanism whereby the inertia weight factor varies during the search. A proper mechanism can improve PSO in exploration and convergence. In this paper we apply the following mechanism, which is proposed by Shi and Eberhart (1998):

(B.3)

w = w_{\min} + (\frac{t_{\max} - t}{t_{\max}}) . (w_{\max} - w_{\min})

where t _max is the maximum number of iterations and t is the current number of iterations.

Appendix C

(C.1)

T_{2} = 0.5616 P_{2} + 0.3126 X_{2} + 48.43

(C.2)

T_{3} = 0.507 P_{2} + 55.0

(C.3)

F_{4} = \frac{Q_{100} - 0.07 F_{1} (T_{2} - T_{1})}{38.5}

(C.4)

T_{100} = 0.1538 P_{100} + 90.0

(C.5)

Q_{100} = 0.16 (F_{1} + F_{3}) (T_{100} - T_{2})

(C.6)

F_{100} = \frac{Q_{100}}{36.6}

(C.7)

Q_{200} = \frac{0.9576 F_{200} (T_{3} - T_{200})}{0.14 F_{200} + 6.84}

(C.8)

T_{201} = T_{200} + \frac{13.68 F_{200} (T_{3} - T_{200})}{0.14 F_{200} + 6.84}

(C.9)

F_{5} = \frac{Q_{200}}{38.5}

Funding

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

References

Almaro

Ferramosca

Gonzalez

. (2014) A gradient-based strategy for the one-layer RTO+MPC controller. Journal of Process Control 24(4): 435–447.

Alvardo

Limon

Munaoz de la

. (2011) Comparative analysis of distributed MPC techniques applied to the HD-MPC four-tank benchmark. Journal of Process Control 21(5): 800–815.

Aswani

Gonzalez

Sastry

. (2013) Provably safe and robust learning-based model predictive control. Automatica 49(5): 1216–1226.

Boggs

Tolle

(2000) Sequential quadratic programming for large scale nonlinear optimization. Journal of Computational and Applied Mathematics 124(1-2): 123–137.

Cao

(2005) A formulation of nonlinear model predictive control using automatic differentiation. Journal of Process Control 15(8): 851–858.

Cruz

ILL

Willingenburg

LGV

Straten

(2003) Efficient Differential Evolution algorithms for multimodal optimal control problems. Applied Soft Computing 3(2): 97–122.

Dalamagkidis

Valavanis

Piegl

(2011) Nonlinear model predictive control with neural network optimization for autonomous autorotation of small unmanned helicopters. IEEE Transactions on Control Systems Technology 19(4): 1047–1056.

Ding

(2011) Dynamic output feedback predictive control for nonlinear systems represented by a Takagi–Sugeno Model. IEEE Transactions on Fuzzy Systems 19(5): 831–843.

Ding

Ping

(2013) Output feedback predictive control with one free control move for nonlinear systems represented by a Takagi–Sugeno Model. IEEE Transactions on Fuzzy Systems 21(5): 1–15.

10.

Kennedy

Eberhart

(1995) Particle swarm optimization. In: proceedings of the IEEE international conference on neural networks, pp.1942–1948.

11.

Limon

Alvarado

Alamo

. (2008) MPC for tracking of piece-wise constant references for constrained linear systems. Automatica 44(9): 2382–2387.

12.

Limon

Alvarado

Alamo

. (2010) Robust tube-based MPC for tracking of constrained linear systems with additive disturbances. Journal of Process Control 20(3): 248–260.

13.

Morari

Maeder

(2012) Nonlinear offset-free model predictive control. Automatica 48(9): 2059–2067.

14.

Newell

Lee

(1989) Applied Process Control: A Case Study. Englewood Cliffs, NJ: Prentice Hall.

15.

Nocedal

Wright

(2006) Numerical Optimization. 2nd ed. Series in Operations Research. New York: Springer-Verlag.

16.

Piche

Rodsari

Johnson

. (2000) Nonlinear model predictive control using neural networks. IEEE Control System Magazine 20(3): 53–62.

17.

Schmitt

(2004) Theory of Genetic Algorithms II: Models for genetic operators over the string-tensor representation of populations and convergence to global optima for arbitrary fitness function under scaling. Theoretical Computer Science 310(1–3): 181–231.

18.

Schweitzer

(1971) Iterative solution of the functional equations of undiscounted Markov renewal programming. Journal of Mathematical Analysis and Applications 34(3): 495–501.

19.

Shen

Xue

Chen

(2010) Global convergence of a robust filter SQP algorithm. European Journal of Operational Research 206(1): 34–45.

20.

Shi

Eberhart

(1998) Parameter selection in particle swarm optimization. In: Evolutionary Programming VII. Lecture Notes in Computer Science. San Diego, California: Springer, pp.591–600.

21.

Sorensen

(2006) Newton’s method with a model trust region modification. SIAM Journal on Numerical Analysis 19(2)(1–3): 409–426.

22.

Stoer

Bulirsch

(2002) Introduction to Numerical Analysis. 3rd ed. Berlin-Heidelberg-New York: Springer.

23.

Vesterstrom

Thomsen

(2004) A comparative study of differential evolution, particle swarm optimization, and evolutionary algorithms on numerical benchmark problems. In: proceedings of the IEEE international conference on evolutionary computation, Volume 2, pp.1980–1987.

24.

Yang

Qian

(2012) A genuine nonlinear approach for controller design of boiler-turbine system. ISA Transactions 51(3): 446–453.

25.

Yang

Biegler

(2013) Advanced-multi-step nonlinear model predictive control. Journal of Process Control 23(8): 1116–1128.