# Convergence Guaranteed Nonlinear Constraint Model Predictive Control via I/O Linearization

Authors: Xiaobing Kong; Xiangjie Liu; Xiuming Yao
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Constituting reliable optimal solution is a key issue for the nonlinear constrained model predictive control. Input-output feedback linearization is a popular method in nonlinear control. By using an input-output feedback linearizing controller, the original linear input constraints will change to nonlinear constraints and sometimes the constraints are state dependent. This paper presents an iterative quadratic program (IQP) routine on the continuous-time system. To guarantee its convergence, another iterative approach is incorporated. The proposed algorithm can reach a feasible solution over the entire prediction horizon. Simulation results on both a numerical example and the continuous stirred tank reactors (CSTR) demonstrate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

Model predictive control (MPC) is a popular algorithm in process control which online solves an optimization problem at each time step [1]. MPC considers process input, output, and state constraints directly in the control variable calculation. Under the linear model, with a quadratic objective function, MPC utilizes the convex quadratic program (QP), which can easily find the optimal solution [2].Industrial processes are generally nonlinear, due to the frequent changes of the operating point right across the whole operation range. In general, the nonlinear model predictive control (NMPC) also online solves an optimization problem, by using the sequential quadratic program (SQP). The resulting nonlinear programming problems are usually nonconvex, and the online computational burden is generally large for most complex systems. A general way to solve this nonlinear optimal problem is to use approximate approach. For example, the first control move is calculated exactly, since it is actually implemented. The rest control moves can be approximated, since they are not implemented [3]. Thus the number of decision variables in the online optimization problem is equal to the number of inputs, instead of the number of inputs multiplied by the control horizon for a conventional NMPC algorithm. Paper [4] extended the first prediction of the linear constraints to the whole control horizon. Paper [5] studies the stability and region of attraction properties of a family of nonlinear MPC systems. Paper [6] presents nonlinear multivariable predictive control using neurofuzzy networks.“Jacobian linearization” and “the input-output feedback linearization (IOFL)” are the two popular approaches in nonlinear control area. While the former can reflect the nonlinear model only at some certain point, the later represents the nonlinear system over a much wider operating range. Thus, IOFL is utilized in this paper for constituting nonlinear MPC, since the IOFL can offer a linear dynamic system so that the total optimal problem can be solved using the QP routine. Nevertheless, this can make the constraints to be nonlinear and state dependent. The approximate method developed in paper [7] can guarantee a feasible solution over the entire prediction horizon. The neural network is utilized to model the nonlinear discrete-time system. Paper [8] also introduces a technique with an affine transformation of the feasible region using feedback linearization scheme for handling input constraints. The neural network is also utilized to model the nonlinear system.For digital MPC controller design, since the real-time control systems are in discrete-time forms, most NMPC research adopts the discrete-time feedback linearization. It should be noticed that the IOFL in differential geometry has been well developed for continuous-time system. Some famous methods in continuous-time system [9], for example, the extended system method, are not suitable to feedback linearise a general discrete-time system. Direct application of feedback linearization for a general discrete-time system may either be impossible or involve the time-consuming search algorithm. In [10], the authors proposed a linear model predictive control strategy via input/output linearization to the nonlinear process, by using one-step constraint algorithm. This paper aims to make full use of the advantage of continuous-time system IOFL and then try to reach the convergence to a feasible solution over the entire prediction horizon within the available time. Since the input constraints are transferred into nonlinear constraints from initial linear constraints, the convergent algorithm is constituted to guarantee a feasible solution without constraints violation. Simulation results on both a numerical example and the continuous stirred tank reactors (CSTR) demonstrate the effectiveness of the NMPC method.


## 2. The Linear Control Structure via Input-Output Feedback Linearization

The IOFL is to transform the original nonlinear system into linear input-output relationship, generally by using a static state feedback control law [11]. Consider the SISO affine state-space model as follows:
(1)x-˙=f(x-)+g(x-)u,y=h(x-),
where x-∈Rn is state variables, u is the manipulated input variable, y is the controlled output variable, and f, g, and h are smooth functions in a domain D⊂Rn.Definition 1 (relative degree).
The nonlinear system (1) is said to have relative degree γ, 1≤γ≤n, in a region D0⊂D if
(2)LgLfih(x-)≡0,i=0,…,γ-2;LgLfγ-1h(x-)≠0,∀x-∈D0.
If the relative degreeγ=n, then for every x0∈D, a neighborhood N of x0 exists such that the map Φ:x-→x=[hLfh⋯Lfn-1h]T restricted to N is a diffeomorphism on N.The mappingΦ transfers the system (1) to a new system as follows:
(3)x˙1=x2,x˙2=x3,x˙2⋮x˙n=b(x)+a(x)u,y=x1,
where a(x-)=LgLfn-1h(x-), b(x-)=Lfnh(x-), which can also be expressed as a(x), b(x) in the new coordinates.The feedback law is(4)u=1a(x)[-b(x)+v],
where v is the transformed input variable. The new linear system is
(5)x˙=Ax+Bv,y=Cx,
where
(6)A=[010⋱⋱⋱100],B=[00⋮1],C=[10⋯0].The resulting linear state-space system (5) could be used for constituting “standard” linear MPC [1]. Discretization (5) can result in
(7)x(k+1)=Adx(k)+Bdv(k),y(k)=Cdx(k).
The matrices Ad, Bd, and Cd are expressed as Ad=eAT, Bd=∫0TeAtdt·B, Cd=C, where T is sampling time.DefineΔx(k+1)=x(k+1)-x(k); Δx(k)=x(k)-x(k-1); Δv(k)=v(k)-v(k-1). Choose a new group of state variables xu(k)=[Δx(k)Ty(k)]T; the augmented system is
(8)[Δx(k+1)y(k+1)]︷xu(k+1)=[Ad0dTCdAd1]︷Au[Δx(k)y(k)]︷xu(k)+[BdCdBd]︷BuΔv(k),y(k)=[0d1]︷Cu[Δx(k)y(k)],
where 0d=[00⋯0]︷n. Based on the state-space model (Au,Bu,Cu), the future outputs along the predictive horizon Np are calculated sequentially as
(9)y(k+1∣k)=CuAuxu(k)+CuBuΔv(k),y(k+2∣k)=CuAu2xu(k)+CuAuBuΔv(k)+CuBuΔv(k+1)⋮(10)y(k+Np∣k)=CuAuNpxu(k)+CuAuNp-1BuΔv(k)+CuAuNp-2BuΔv(k+1)+⋯+CuAuNp-NcBuΔv(k+Nc-1).
Define vectors
(11)Y=[y(k+1∣k)⋯y(k+Np∣k)]T,ΔV=[Δv(k)Δv(k+1∣k)⋯Δv(k+Nc-1∣k)]T,
where Nc is the control horizon.Collect (9) and (10) together in a compact matrix form as
(12)Y=Fxu(k)+ΨΔV,
where(13)F=[CuAuCuAu2CuAu3⋮CuAuNp],Ψ=[CuBu00⋯0CuAuBuCuBu0⋯0CuAu2BuCuAuBuCuBu⋯0⋮CuAuNp-1BuCuAuNp-2BuCuAuNp-3Bu⋯CuAuNp-NcBu].Similarly, we can have the future state variables along the predictive horizonNc based on the model of (Ad,Bd,Cd)(14)X=A~ΔV+γ,
where
(15)X=[x(k)x(k+1∣k)⋯x(k+Nc-1∣k)]T,A~=[00⋯0Bd0⋯0(Ad+I)BdBd⋯0⋮⋮⋮⋮(∑i=1Nc-1Adi-1)Bd(∑i=1Nc-2Adi-1)Bd⋯0],γ=[x(k)x(k)+AdΔx(k)x(k)+(∑i=12Adi)Δx(k)⋮x(k)+(∑i=1Nc-1Adi)Δx(k)].
The object function to be minimised is a quadratic criterion on ΔV as
(16)J=(Rs-Y)T(Rs-Y)+ΔVTR-ΔV.
By using (12), the optimization problem can be attributed to
(17)Jmin=minΔV{12ΔVTHΔV+ΓTΔV}subjecttoΔV_≤ΔV≤ΔV¯,
where H=ΨTΨ+R-, Γ=ΨT(-Rs+Fxu(k)), ΔV_ and ΔV¯ indicates the minimal and maximum value of ΔV.


## 3. Convergence Algorithm for Constraint Optimal MPC

### 3.1. The Nonlinear Constraint Handling

In practice, the process inputs are frequently subjected to the following level inequality constraints:(18)U_≤U≤U-,
where U=[uk⋯uk+Nc-1]T represents the control inputs, and U_ and U¯ indicate the minimal and maximum value of U.After the realization of the IOFL, the input vectorU is transformed into a new one V. With this nonlinear and state dependent input constraint, traditional QP cannot be used to calculate the optimal control sequence. Though the sequential quadratic programming (SQP) technique is generally used on nonlinear constraint MPC, it is often nonconvex and can cause large computation burden. This property will be demonstrated later in the simulation example of CSTR.In order to solve this nonlinear constraints problem, an iterative QP is first used to try to get the optimal solution. The key issue is to establish the nonlinear relationship between vectorsV and U and then get the explicit expression of U from V.To establish the nonlinear relationship between vectorsV and U, expand the feedback law (4) over the control horizon Nc,
(19)v(k)=a(x(k))u(k)+b(x(k)),v(k+1)=a(x(k+1))u(k+1)+b(x(k+1)),⋮v(k+Nc-1)=a(x(k+Nc-1))u(k+Nc-1)+b(x(k+Nc-1)).Considering the expression (14), (19) can be rearranged in the matrix form
(20)V=Gu[U,ΔV].
The above inequality can be rewritten in the following form:
(21)G[U,ΔV]=0
due to
(22)V=[v(k)v(k+1)⋮v(k+Nc-1)]=[1011111⋮⋱⋱⋱11⋯11]·[Δv(k)Δv(k+1)⋮Δv(k+Nc-1)]+v(k-1).
In general, v(k+i-1) can be rewritten as
(23)v(k+i-1)=v(k-1)+∑j=1iΔv(k+j-1)(i=1,…,Nc),
with the constraints
(24)v_≤v≤v¯,
where v_, v¯ indicate the minimal and maximum value of v, which are state dependent.Thus(25)v_k+i-1(x(k+i-1))-v(k-1)≤∑j=1iΔv(k+j-1)≤v-k+i-1(x(k+i-1))-v(k-1).
The above inequality can be rewritten over the entire horizon,
(26)ΛTΔV(k)≤c(X(ΔV(k)))T,
where(27)Λ=[LT-LT],L=[1011⋮11⋯1],c=[v-k-v(k-1)⋯v-k+Nc-1-v(k-1)v(k-1)-v_k⋯v(k-1)-v_k+Nc-1].From (26), the above equation can be rewritten as
(28)ΛTΔV(k)≤c(X(ΔV(U(k))))T.This iterative QP algorithm is presented below:Step  1: Initializingc(X(ΔV(U0(k)))) within the input constraint.Step  2: Solve QP forΔVi(k) subject to ΛTΔVi(k)≤c(X(ΔV(Ui-1(k))))TStep  3: Calculatec(X(ΔV(Ui(k))))Step  4: Test ifΛTΔVi(k)≤c(X(ΔV(Ui(k))))T, end; otherwise, i=i+1, go to Step  2.The above iterative QP is effective for solving the problem if the initial condition is properly chosen. Nevertheless, it has no guaranteed convergence to a feasible solution. To overcome this problem, another iterative algorithm is needed to guarantee a feasible solution over the entire prediction horizon is incorporated. This idea is previously presented by paper [7] on a input/output model. Extending it to the state-space equation is straightforward.In order to solve the convergence problem, linear approximation is necessary. Now consider the approximation by linearization through Taylor’s expansion aboutU0 of the nonlinear constraints (21) as
(29)ΔV=ΔV0+g[U0](h~)+r(h~),
where U0 is the chosen initial operating trajectory, ΔV0 is obtained from (21) when U=U0, the (Nc×Nc) dimensional matrix g[U0] represents the Jacobian matrix ∂ΔV/∂U in (21) at the operating trajectory, and r(h~) corresponds to the higher order terms of the approximation given by
(30)r(h~)=12h~Tg[U0+θh~]h~,
where θ∈[0,1], h~=U-U0.Neglect the higher order terms as(31)ΔV=ΔV0+g[U0](U-U0).
Equation (31) can lead to an explicitly linear function relationship between the control sequence U and ΔV:
(32)U=g-1[U0]ΔV+U0-g-1[U0]ΔV0.
Now the optimization problem is defined as
(33)Jmin=minΔV{12ΔVTHΔV+ΓTΔV}.
subject to the convex set of approximate linear constraints
(34)U_≤MΔV+m0≤U-,
with M=g-1[U0], m0=U0-g-1[U0]×ΔV0.In order for the algorithm to reach the convergence to a feasible solution, the optimal problem can be rewritten as(35)ΔV*=argminΔV{12ΔVTHΔV+ΓTΔV}(36)subjecttoU_≤MαΔV+m0α≤U-,
where 0<α<1 is decreasing parameter with initial value α=1. Define Mα and m0α, respectively, by
(37)Mα=(g-1[U0])α(I-(g-1[U0])α)-1,m0α=U0-(g-1[U0])α(I-(g-1[U0])α)-1ΔV0.
Define Ul*=MαΔV*+m0α to be the linear control sequence and Unl* to be the nonlinear control sequence obtained from G[U,ΔV*]=0. Then the major goal is to keep
(38)U_≤Unl*≤U-.
Thus the optimization problem with the approximate constraints ΔV* is a feasible solution, which can be solved iteratively using QP.This approximate constraints algorithm is presented follows.Step  1: assure feasibility of the new operating trajectory:U_≤Ui≤U-. Solve for ΔV0, M,m0Step  2:αi-1=1Step  3: actualize the squeezing factorα:
(39)α=α×λ.Step  4: new linearized constraints atU0: U_≤MαΔV+m0α≤U-.Step  5: Solve QP for the new linearized constraints:(40)ΔVi+1*=argminΔVJ(ΔV).Step  6: test ifUnl* violate the original constraint. If: U_≤Unl*≤U-⇒End.Step  7: leti=i+1, go to 3.


### 3.2. The Convergence Property

With the object function (35), assume the practical input sequence constraints
(41)|U|≤1.
Assume a given initial vector representing the operating trajectory, U0, verifying the conditions
(42)|U0|≤1-ε1,
where ε1 is small positive scalar values, ε1>0.For any given scalar0<α<1, the solution of the optimization problem (35) is subjected to
(43)|U0+(g-1[U0])α(I-(g-1[U0])α)(ΔV*-ΔV0)|≤1.
which can be further written as
(44)|ΔV*|≤[I-(g-1[U0])-α](g[U0])α(1-U0)+ΔV0.
In the limit, as α→0, (g-1[U0])-α→I, (g[U0])α→I, the solution of the optimization problem is given by
(45)limα→0|ΔV*|≤ΔV0.
Due to the smooth function G[U,ΔV]=0 and definition (42), this last result leads to
(46)|U*|=|U0|≤1-ε1,
signifying that the resulting practical nonlinear control sequence is always within the bounds, and so a feasible solution is always guaranteed.


### 3.3. The Algorithm for Constraint Optimal MPC

The MPC algorithm, which can guarantee convergence to a feasible solution, has to deal with the constraints as the following steps:Step  1:i=0.Step  2: if(N-i<0), go to 8.Step  3: assure feasibility of the new operating trajectory:U_≤Ui≤U-.Step  4: solve QP atUi:
(47)ΔVi+1*=argminΔVJ(ΔV)(48)subjecttoΛTΔVi+1*≤c(X(ΔV(Ui)))T.Step  5: test ifΔVnl* violates the constraints: if ΛTΔV≤c(X(ΔV(Ui+1)))T⇒END
whereUi+1=Unl* is nonlinear control sequence, obtained from G[Unl*,ΔVi+1*]=0.Step  6:i=i+1. Go to 2Step  7:U0=Ui, ΔV0=ΔVN*. Solve for M, m0.Step  8:α=1.Step  9: actualize the squeezing factorα:
(49)α=α×λ.Step  10: new linearized constraints atU0: U_≤MαΔV+m0α≤U-Step  11: Solve QP for the new linearized constraints:(50)ΔVi+1*=argminΔVJ(ΔV).Step  12: test ifUnl* violate the original constraint. If: U_≤Unl*≤U-⇒EndStep  13: leti=i+1. Go to 9.The total algorithm performance is quite related to the tuning of two parameters:N: the number of iterations of the QP algorithm;λ: the squeezing parameter of the convergence guaranteed algorithm.The convergence speed of this iterative algorithm can be set using the decreasing rateα, defined here through parameter 0<λ<1.Choosing the initial valueU0 is a key issue in solving the optimal problem, since it is quite related to the resulting computation burden. Suppose that the final input sequence at instant k can be expressed as Unl*=[u(0)u(1)⋯u(Nc-2)]T. Then, the initial value U0 at instant k+1 should be chosen as  U0=[u(1)⋯u(Nc-2)0]T; for example, the first Nc-1 factors are taken directly from the prediction of Unl*.


## 3.1. The Nonlinear Constraint Handling

In practice, the process inputs are frequently subjected to the following level inequality constraints:(18)U_≤U≤U-,
where U=[uk⋯uk+Nc-1]T represents the control inputs, and U_ and U¯ indicate the minimal and maximum value of U.After the realization of the IOFL, the input vectorU is transformed into a new one V. With this nonlinear and state dependent input constraint, traditional QP cannot be used to calculate the optimal control sequence. Though the sequential quadratic programming (SQP) technique is generally used on nonlinear constraint MPC, it is often nonconvex and can cause large computation burden. This property will be demonstrated later in the simulation example of CSTR.In order to solve this nonlinear constraints problem, an iterative QP is first used to try to get the optimal solution. The key issue is to establish the nonlinear relationship between vectorsV and U and then get the explicit expression of U from V.To establish the nonlinear relationship between vectorsV and U, expand the feedback law (4) over the control horizon Nc,
(19)v(k)=a(x(k))u(k)+b(x(k)),v(k+1)=a(x(k+1))u(k+1)+b(x(k+1)),⋮v(k+Nc-1)=a(x(k+Nc-1))u(k+Nc-1)+b(x(k+Nc-1)).Considering the expression (14), (19) can be rearranged in the matrix form
(20)V=Gu[U,ΔV].
The above inequality can be rewritten in the following form:
(21)G[U,ΔV]=0
due to
(22)V=[v(k)v(k+1)⋮v(k+Nc-1)]=[1011111⋮⋱⋱⋱11⋯11]·[Δv(k)Δv(k+1)⋮Δv(k+Nc-1)]+v(k-1).
In general, v(k+i-1) can be rewritten as
(23)v(k+i-1)=v(k-1)+∑j=1iΔv(k+j-1)(i=1,…,Nc),
with the constraints
(24)v_≤v≤v¯,
where v_, v¯ indicate the minimal and maximum value of v, which are state dependent.Thus(25)v_k+i-1(x(k+i-1))-v(k-1)≤∑j=1iΔv(k+j-1)≤v-k+i-1(x(k+i-1))-v(k-1).
The above inequality can be rewritten over the entire horizon,
(26)ΛTΔV(k)≤c(X(ΔV(k)))T,
where(27)Λ=[LT-LT],L=[1011⋮11⋯1],c=[v-k-v(k-1)⋯v-k+Nc-1-v(k-1)v(k-1)-v_k⋯v(k-1)-v_k+Nc-1].From (26), the above equation can be rewritten as
(28)ΛTΔV(k)≤c(X(ΔV(U(k))))T.This iterative QP algorithm is presented below:Step  1: Initializingc(X(ΔV(U0(k)))) within the input constraint.Step  2: Solve QP forΔVi(k) subject to ΛTΔVi(k)≤c(X(ΔV(Ui-1(k))))TStep  3: Calculatec(X(ΔV(Ui(k))))Step  4: Test ifΛTΔVi(k)≤c(X(ΔV(Ui(k))))T, end; otherwise, i=i+1, go to Step  2.The above iterative QP is effective for solving the problem if the initial condition is properly chosen. Nevertheless, it has no guaranteed convergence to a feasible solution. To overcome this problem, another iterative algorithm is needed to guarantee a feasible solution over the entire prediction horizon is incorporated. This idea is previously presented by paper [7] on a input/output model. Extending it to the state-space equation is straightforward.In order to solve the convergence problem, linear approximation is necessary. Now consider the approximation by linearization through Taylor’s expansion aboutU0 of the nonlinear constraints (21) as
(29)ΔV=ΔV0+g[U0](h~)+r(h~),
where U0 is the chosen initial operating trajectory, ΔV0 is obtained from (21) when U=U0, the (Nc×Nc) dimensional matrix g[U0] represents the Jacobian matrix ∂ΔV/∂U in (21) at the operating trajectory, and r(h~) corresponds to the higher order terms of the approximation given by
(30)r(h~)=12h~Tg[U0+θh~]h~,
where θ∈[0,1], h~=U-U0.Neglect the higher order terms as(31)ΔV=ΔV0+g[U0](U-U0).
Equation (31) can lead to an explicitly linear function relationship between the control sequence U and ΔV:
(32)U=g-1[U0]ΔV+U0-g-1[U0]ΔV0.
Now the optimization problem is defined as
(33)Jmin=minΔV{12ΔVTHΔV+ΓTΔV}.
subject to the convex set of approximate linear constraints
(34)U_≤MΔV+m0≤U-,
with M=g-1[U0], m0=U0-g-1[U0]×ΔV0.In order for the algorithm to reach the convergence to a feasible solution, the optimal problem can be rewritten as(35)ΔV*=argminΔV{12ΔVTHΔV+ΓTΔV}(36)subjecttoU_≤MαΔV+m0α≤U-,
where 0<α<1 is decreasing parameter with initial value α=1. Define Mα and m0α, respectively, by
(37)Mα=(g-1[U0])α(I-(g-1[U0])α)-1,m0α=U0-(g-1[U0])α(I-(g-1[U0])α)-1ΔV0.
Define Ul*=MαΔV*+m0α to be the linear control sequence and Unl* to be the nonlinear control sequence obtained from G[U,ΔV*]=0. Then the major goal is to keep
(38)U_≤Unl*≤U-.
Thus the optimization problem with the approximate constraints ΔV* is a feasible solution, which can be solved iteratively using QP.This approximate constraints algorithm is presented follows.Step  1: assure feasibility of the new operating trajectory:U_≤Ui≤U-. Solve for ΔV0, M,m0Step  2:αi-1=1Step  3: actualize the squeezing factorα:
(39)α=α×λ.Step  4: new linearized constraints atU0: U_≤MαΔV+m0α≤U-.Step  5: Solve QP for the new linearized constraints:(40)ΔVi+1*=argminΔVJ(ΔV).Step  6: test ifUnl* violate the original constraint. If: U_≤Unl*≤U-⇒End.Step  7: leti=i+1, go to 3.


## 3.2. The Convergence Property

With the object function (35), assume the practical input sequence constraints
(41)|U|≤1.
Assume a given initial vector representing the operating trajectory, U0, verifying the conditions
(42)|U0|≤1-ε1,
where ε1 is small positive scalar values, ε1>0.For any given scalar0<α<1, the solution of the optimization problem (35) is subjected to
(43)|U0+(g-1[U0])α(I-(g-1[U0])α)(ΔV*-ΔV0)|≤1.
which can be further written as
(44)|ΔV*|≤[I-(g-1[U0])-α](g[U0])α(1-U0)+ΔV0.
In the limit, as α→0, (g-1[U0])-α→I, (g[U0])α→I, the solution of the optimization problem is given by
(45)limα→0|ΔV*|≤ΔV0.
Due to the smooth function G[U,ΔV]=0 and definition (42), this last result leads to
(46)|U*|=|U0|≤1-ε1,
signifying that the resulting practical nonlinear control sequence is always within the bounds, and so a feasible solution is always guaranteed.


## 3.3. The Algorithm for Constraint Optimal MPC

The MPC algorithm, which can guarantee convergence to a feasible solution, has to deal with the constraints as the following steps:Step  1:i=0.Step  2: if(N-i<0), go to 8.Step  3: assure feasibility of the new operating trajectory:U_≤Ui≤U-.Step  4: solve QP atUi:
(47)ΔVi+1*=argminΔVJ(ΔV)(48)subjecttoΛTΔVi+1*≤c(X(ΔV(Ui)))T.Step  5: test ifΔVnl* violates the constraints: if ΛTΔV≤c(X(ΔV(Ui+1)))T⇒END
whereUi+1=Unl* is nonlinear control sequence, obtained from G[Unl*,ΔVi+1*]=0.Step  6:i=i+1. Go to 2Step  7:U0=Ui, ΔV0=ΔVN*. Solve for M, m0.Step  8:α=1.Step  9: actualize the squeezing factorα:
(49)α=α×λ.Step  10: new linearized constraints atU0: U_≤MαΔV+m0α≤U-Step  11: Solve QP for the new linearized constraints:(50)ΔVi+1*=argminΔVJ(ΔV).Step  12: test ifUnl* violate the original constraint. If: U_≤Unl*≤U-⇒EndStep  13: leti=i+1. Go to 9.The total algorithm performance is quite related to the tuning of two parameters:N: the number of iterations of the QP algorithm;λ: the squeezing parameter of the convergence guaranteed algorithm.The convergence speed of this iterative algorithm can be set using the decreasing rateα, defined here through parameter 0<λ<1.Choosing the initial valueU0 is a key issue in solving the optimal problem, since it is quite related to the resulting computation burden. Suppose that the final input sequence at instant k can be expressed as Unl*=[u(0)u(1)⋯u(Nc-2)]T. Then, the initial value U0 at instant k+1 should be chosen as  U0=[u(1)⋯u(Nc-2)0]T; for example, the first Nc-1 factors are taken directly from the prediction of Unl*.


## 4. Case Studies

Two examples are presented to illustrate the implementation and the performance of the proposed nonlinear MPC. In example 1, a numerical state-space equation is used, and in example 2, the control problem of CSTR is considered.


### 4.1. Numerical Example

Assume the state-space equations to be(51)x-˙1=x-2,x-˙2=-3x-12x-2-x-13sin(x-2)+u,y=x-1.The input constraint is-0.5≤u≤0.5.Using (4), the static feedback law: v=u+(-3x12x2-x13sinx2)The resulting feedback linearization linear system is(52)[x˙1x˙2]=[0100][x1x2]+[01]v,y=x1y=x1.Choose the predictive horizonNp=8, the control horizon Nc=5, and the sampling time τ=1 s. Figure 1 shows the system response for a step signal using the proposed IOFL MPC method. The algorithm guarantee’s that the process input does not violate the constraints. The system reaches stability in 3 seconds without overshoot. Figure 2 shows the number of iterations at each optimization step. From Figure 2, the process needs 65 iterations at the beginning but soon reduces the iteration and finally reaches the convergence using only one iteration.Figure 1
The closed-loop response for a step signal under the proposed IOFL MPC.Figure 2
The number of iterations at each optimization step.Simulations were then repeated under different predictive horizons. Table1 lists the comparison of the control performance. From Table 1, the computing burden increases when the prediction horizon increases, meanwhile the tracking property improves with the increment of the prediction horizon.Table 1
The comparison of control performance at different prediction horizon.
Prediction horizon
Computing burden(CPU time)
Sum squared error
8
1.0469
11.4042
10
1.1875
11.4036
12
1.2813
11.4032


### 4.2. The Control of a First-Order CSTR

CSTRs are typical chemical processes representing a variety of complex industry systems. The CSTR problem discussed here represents a first-order, irreversible, exothermic kinetics reaction, which can be described by the following equation [12]:
(53)dx1dt=-x1+Dα·exp(x21+x2/γ),dx2dt=-x2(1+β)+H·Dα·(1-x1)·exp(x21+x2/γ)+βu,y=x1,
where the two state variables x1 and x2 are the reactor temperature and the normalized reactant concentration. The control variable u is the normalized cooling water temperature. Dα, β, and γ are the constants of the system.The CSTR control objectives contain two aspects: set point tracking and regulation of reactant concentration under the perturbation of feed water temperature.With the static feedback law,(54)v=au+b,
where
(55)a=β(ξ1+ξ2)·(1-(1r)ln((ξ1+ξ2)Dα))2,b=ξ1-(ξ1+ξ2)+(ξ1+ξ2)·[ln((ξ1+ξ2)/Dα)1-(1/r)ln((ξ1+ξ2)/Dα)(1+β)+H·(1-ξ1)(ξ1+ξ2)ln((ξ1+ξ2)/Dα)1-(1/r)ln((ξ1+ξ2)/Dα)].
The resulting feedback linearization linear system is
(56)[ξ˙1ξ˙2]=[0100][ξ1ξ2]+[01]v,y=[10][ξ1ξ2].Suppose the constraint to be-3≤u≤3. Let r,Dα,H=1, β=0.01, the predictive horizon Np=10, the control horizon Nc=5, and the sampling time τ=1 s. Figure 3 shows the system output and the control under step changes in the reference signal.Figure 3
The system output and the control under step changes in the reference signal.From the simulation result, the manipulated variableu does not violate the defined constraint. The system reaches stability in 30 steps without obvious overshoot.As mentioned before, there exist two general approaches for NMPC. One is to use the SQP method. Another is to approximate the first constraints as the constraints over the entire control horizon. The proposed method is then compared with these two methods, considering this CSTR plant. The simulation results using the three predictive control schemes forNc=5 are given in Figure 4, adopting the best N and λ for each Np in the proposed IOFL NMPC.Figure 4
The comparison of the three responses.Computation burden is one of the most concerning problems in CSTR real-time control. In using MPC, though choosing larger predictive horizon can improve system performance, it can, meantime, result in the increment of computation burden, especially in nonlinear system. For further comparison, the controller performance is analysed with respect to two variables, for example, sum squared output error (SSE) and the relative optimization time for the simulation. Table2 shows the performance of each controller configuration for various prediction horizons. The three methods, namely, the one-step constraints, the MPC with IOFL and the nonlinear MPC, can all guarantee a feasible control solution over the complete prediction horizon. The comparison of the three shows a similar trend performance in terms of the system tracking error, although the one-step constraint is less computational time demanding, as shown in Figure 5. With the increase of the prediction horizon, the time increase in nonlinear MPC is obvious, while the time increase in one-step constraints and MPC with IOFL optimization is quite minor. This is due to the efficient iterative process of the proposed IOFL MPC. The reduced computing burden for the one-step constraints method is obvious, since only the first time-step-ahead inputs are constrained, leaving the other Np-1 inputs unconstrained. The shortcoming is that the resulting optimum may not be feasible. From Figure 5, the total trade-off between the computational demands and the gain in optimality is obviously favourable to the proposed MPC with IOFL.Table 2
Left: CPU time of the entire simulation. Right: sum squared output error (SSE) in the simulation.
Predictive horizon
CPU time
SSE
One-step constraints
IOFL MPC
SQP
One-step constraints
IOFL MPC
SQP
10
0.1563
0.2011
0.3054
11.4181
11.4045
11.4042
12
0.1719
0.2407
0.3602
11.4180
11.4034
11.4031
14
0.2031
0.2795
0.4415
11.4177
11.4026
11.4018
16
0.2614
0.3306
0.5728
11.4130
11.4019
11.4003Figure 5
The comparison of control performance.


## 4.1. Numerical Example

Assume the state-space equations to be(51)x-˙1=x-2,x-˙2=-3x-12x-2-x-13sin(x-2)+u,y=x-1.The input constraint is-0.5≤u≤0.5.Using (4), the static feedback law: v=u+(-3x12x2-x13sinx2)The resulting feedback linearization linear system is(52)[x˙1x˙2]=[0100][x1x2]+[01]v,y=x1y=x1.Choose the predictive horizonNp=8, the control horizon Nc=5, and the sampling time τ=1 s. Figure 1 shows the system response for a step signal using the proposed IOFL MPC method. The algorithm guarantee’s that the process input does not violate the constraints. The system reaches stability in 3 seconds without overshoot. Figure 2 shows the number of iterations at each optimization step. From Figure 2, the process needs 65 iterations at the beginning but soon reduces the iteration and finally reaches the convergence using only one iteration.Figure 1
The closed-loop response for a step signal under the proposed IOFL MPC.Figure 2
The number of iterations at each optimization step.Simulations were then repeated under different predictive horizons. Table1 lists the comparison of the control performance. From Table 1, the computing burden increases when the prediction horizon increases, meanwhile the tracking property improves with the increment of the prediction horizon.Table 1
The comparison of control performance at different prediction horizon.
Prediction horizon
Computing burden(CPU time)
Sum squared error
8
1.0469
11.4042
10
1.1875
11.4036
12
1.2813
11.4032


## 4.2. The Control of a First-Order CSTR

CSTRs are typical chemical processes representing a variety of complex industry systems. The CSTR problem discussed here represents a first-order, irreversible, exothermic kinetics reaction, which can be described by the following equation [12]:
(53)dx1dt=-x1+Dα·exp(x21+x2/γ),dx2dt=-x2(1+β)+H·Dα·(1-x1)·exp(x21+x2/γ)+βu,y=x1,
where the two state variables x1 and x2 are the reactor temperature and the normalized reactant concentration. The control variable u is the normalized cooling water temperature. Dα, β, and γ are the constants of the system.The CSTR control objectives contain two aspects: set point tracking and regulation of reactant concentration under the perturbation of feed water temperature.With the static feedback law,(54)v=au+b,
where
(55)a=β(ξ1+ξ2)·(1-(1r)ln((ξ1+ξ2)Dα))2,b=ξ1-(ξ1+ξ2)+(ξ1+ξ2)·[ln((ξ1+ξ2)/Dα)1-(1/r)ln((ξ1+ξ2)/Dα)(1+β)+H·(1-ξ1)(ξ1+ξ2)ln((ξ1+ξ2)/Dα)1-(1/r)ln((ξ1+ξ2)/Dα)].
The resulting feedback linearization linear system is
(56)[ξ˙1ξ˙2]=[0100][ξ1ξ2]+[01]v,y=[10][ξ1ξ2].Suppose the constraint to be-3≤u≤3. Let r,Dα,H=1, β=0.01, the predictive horizon Np=10, the control horizon Nc=5, and the sampling time τ=1 s. Figure 3 shows the system output and the control under step changes in the reference signal.Figure 3
The system output and the control under step changes in the reference signal.From the simulation result, the manipulated variableu does not violate the defined constraint. The system reaches stability in 30 steps without obvious overshoot.As mentioned before, there exist two general approaches for NMPC. One is to use the SQP method. Another is to approximate the first constraints as the constraints over the entire control horizon. The proposed method is then compared with these two methods, considering this CSTR plant. The simulation results using the three predictive control schemes forNc=5 are given in Figure 4, adopting the best N and λ for each Np in the proposed IOFL NMPC.Figure 4
The comparison of the three responses.Computation burden is one of the most concerning problems in CSTR real-time control. In using MPC, though choosing larger predictive horizon can improve system performance, it can, meantime, result in the increment of computation burden, especially in nonlinear system. For further comparison, the controller performance is analysed with respect to two variables, for example, sum squared output error (SSE) and the relative optimization time for the simulation. Table2 shows the performance of each controller configuration for various prediction horizons. The three methods, namely, the one-step constraints, the MPC with IOFL and the nonlinear MPC, can all guarantee a feasible control solution over the complete prediction horizon. The comparison of the three shows a similar trend performance in terms of the system tracking error, although the one-step constraint is less computational time demanding, as shown in Figure 5. With the increase of the prediction horizon, the time increase in nonlinear MPC is obvious, while the time increase in one-step constraints and MPC with IOFL optimization is quite minor. This is due to the efficient iterative process of the proposed IOFL MPC. The reduced computing burden for the one-step constraints method is obvious, since only the first time-step-ahead inputs are constrained, leaving the other Np-1 inputs unconstrained. The shortcoming is that the resulting optimum may not be feasible. From Figure 5, the total trade-off between the computational demands and the gain in optimality is obviously favourable to the proposed MPC with IOFL.Table 2
Left: CPU time of the entire simulation. Right: sum squared output error (SSE) in the simulation.
Predictive horizon
CPU time
SSE
One-step constraints
IOFL MPC
SQP
One-step constraints
IOFL MPC
SQP
10
0.1563
0.2011
0.3054
11.4181
11.4045
11.4042
12
0.1719
0.2407
0.3602
11.4180
11.4034
11.4031
14
0.2031
0.2795
0.4415
11.4177
11.4026
11.4018
16
0.2614
0.3306
0.5728
11.4130
11.4019
11.4003Figure 5
The comparison of control performance.


## 5. Conclusions

MPC has been widely used in industrial process. One limitation of MPC is that it is mostly based on a linear model. The performance will deteriorate when the current output is relatively far away from the operating point at which the linear control model was generated, especially in the presence of constrains. The approach presented in this paper aims at combining the nonlinear differential geometry technique with existing MPC technique to avoid nonlinear optimization when using a nonlinear model of the process. Considering the continuous-time system, the detailed routine for reaching a feasible solution through iterative QP optimization is presented in this paper. The initial value is properly chosen to greatly reduce the computation burden. Comparing results on the CSTR considering, both the SSE and the relative optimization time needed for completing the simulation have also been addressed in detail. The advantage is that it avoids a direct nonlinear optimization on a nonlinear model of the process and incorporates the linear optimization technique in the controller design.



---

*Source: 476367-2013-12-05.xml*