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

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.

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) where 𝑥∈𝑅^𝑛 is state variables, 𝑢 is the manipulated input variable, 𝑦 is the controlled output variable, and 𝑓, 𝑔, and ℎ are smooth functions in a domain 𝐷⊂𝑅^𝑛.

Definition 1 (relative degree). The nonlinear system (1) is said to have relative degree 𝛾, 1≤𝛾≤𝑛, in a region 𝐷₀⊂𝐷 if 𝐿_𝑔𝐿^𝑖_𝑓ℎ(𝑥)𝐿≡0,𝑖=0,…,𝛾−2;_𝑔𝐿_𝑓^𝛾−1ℎ(𝑥)≠0,∀𝑥∈𝐷₀. (2)
If the relative degree 𝛾=𝑛, then for every 𝑥₀∈𝐷, a neighborhood 𝑁 of 𝑥₀ exists such that the map Φ∶𝑥→𝑥=[ℎ𝐿_𝑓ℎ⋯𝐿_𝑓^𝑛−1ℎ]^𝑇 restricted to 𝑁 is a diffeomorphism on 𝑁.

The mapping Φ transfers the system (1) to a new system as follows: 𝑥̇₁=𝑥₂,𝑥̇₂=𝑥₃,⋮𝑥̇_𝑛=𝑏(𝑥)+𝑎(𝑥)𝑢,𝑦=𝑥₁, (3) where 𝑎(𝑥)=𝐿_𝑔𝐿_𝑓^𝑛−1ℎ(𝑥), 𝑏(𝑥)=𝐿^𝑛_𝑓ℎ(𝑥), which can also be expressed as 𝑎(𝑥), 𝑏(𝑥) in the new coordinates.

The feedback law is 1𝑢=𝑎(𝑥)[−𝑏(𝑥)+𝑣], (4) where 𝑣 is the transformed input variable. The new linear system is 𝑥̇=𝐴𝑥+𝐵𝑣,𝑦=𝐶𝑥, (5) where 00⋮1𝐴=010⋱⋱⋱100,𝐵=,𝐶=10⋯0. (6)

The resulting linear state-space system (5) could be used for constituting “standard” linear MPC [1]. Discretization (5) can result in 𝑥(𝑘+1)=𝐴_𝑑𝑥(𝑘)+𝐵_𝑑𝑣(𝑘),𝑦(𝑘)=𝐶_𝑑𝑥(𝑘). (7) The matrices 𝐴_𝑑, 𝐵_𝑑, and 𝐶_𝑑 are expressed as 𝐴_𝑑=𝑒^𝐴𝑇, 𝐵_𝑑=∫^𝑇₀𝑒^𝐴𝑡𝑑𝑡⋅𝐵, 𝐶_𝑑=𝐶, where 𝑇 is sampling time.

Define Δ𝑥(𝑘+1)=𝑥(𝑘+1)−𝑥(𝑘); Δ𝑥(𝑘)=𝑥(𝑘)−𝑥(𝑘−1); Δ𝑣(𝑘)=𝑣(𝑘)−𝑣(𝑘−1). Choose a new group of state variables 𝑥_𝑢(𝑘)=[Δ𝑥(𝑘)^𝑇𝑦(𝑘)]^𝑇; the augmented system is ^𝑥𝑢(^()𝑘+1)𝑦()=Δ𝑥𝑘+1𝑘+1^𝐴𝑢𝐴_𝑑0^𝑇_𝑑𝐶_𝑑𝐴_𝑑1^{𝑥𝑢(𝑘)}(𝑘)𝑦(𝑘)+Δ𝑥𝐵^𝐵𝑢_𝑑𝐶_𝑑𝐵_𝑑(𝑘),𝑦(𝑘)=Δ𝑣_𝐶𝑢0_𝑑1(𝑘)𝑦(𝑘),Δ𝑥 (8) where 0_𝑑=[00^𝑛⋯0]. Based on the state-space model (𝐴_𝑢,𝐵_𝑢,𝐶_𝑢), the future outputs along the predictive horizon 𝑁_𝑝 are calculated sequentially as 𝑦(𝑘+1∣𝑘)=𝐶_𝑢𝐴_𝑢𝑥_𝑢(𝑘)+𝐶_𝑢𝐵_𝑢(𝑘),𝑦()Δ𝑣𝑘+2∣𝑘=𝐶_𝑢𝐴²_𝑢𝑥_𝑢(𝑘)+𝐶_𝑢𝐴_𝑢𝐵_𝑢Δ𝑣(𝑘)+𝐶_𝑢𝐵_𝑢()⋮𝑦Δ𝑣𝑘+1𝑘+𝑁_𝑝∣𝑘=𝐶_𝑢𝐴^𝑁𝑝_𝑢𝑥_𝑢(𝑘)+𝐶_𝑢𝐴^𝑁𝑝_𝑢⁻¹𝐵_𝑢Δ𝑣(𝑘)+𝐶_𝑢𝐴^𝑁𝑝_𝑢⁻²𝐵_𝑢Δ𝑣(𝑘+1)+⋯+𝐶_𝑢𝐴_𝑢^{𝑁𝑝−𝑁𝑐}𝐵_𝑢Δ𝑣𝑘+𝑁_𝑐−1. (9)(10) Define vectors 𝑌=𝑦(𝑘+1∣𝑘)⋯𝑦𝑘+𝑁_𝑝∣𝑘^𝑇,Δ𝑉=Δ𝑣(𝑘)Δ𝑣(𝑘+1∣𝑘)⋯Δ𝑣𝑘+𝑁_𝑐−1∣𝑘^𝑇, (11) where 𝑁_𝑐 is the control horizon.

Collect (9) and (10) together in a compact matrix form as 𝑌=𝐹𝑥_𝑢(𝑘)+ΨΔ𝑉, (12) where𝐶𝐹=^𝑢𝐴_𝑢𝐶_𝑢𝐴_𝑢²𝐶_𝑢𝐴_𝑢³⋮𝐶_𝑢𝐴_𝑢^𝑁𝑝𝐶,Ψ=^𝑢𝐵_𝑢𝐶00⋯0_𝑢𝐴_𝑢𝐵_𝑢𝐶_𝑢𝐵_𝑢𝐶0⋯0_𝑢𝐴_𝑢²𝐵_𝑢𝐶_𝑢𝐴_𝑢𝐵_𝑢𝐶_𝑢𝐵_𝑢⋮𝐶⋯0_𝑢𝐴_𝑢^𝑁𝑝−1𝐵_𝑢𝐶_𝑢𝐴_𝑢^𝑁𝑝−2𝐵_𝑢𝐶_𝑢𝐴_𝑢^𝑁𝑝−3𝐵_𝑢⋯𝐶_𝑢𝐴_𝑢^{𝑁𝑝−𝑁𝑐}𝐵_𝑢. (13)

Similarly, we can have the future state variables along the predictive horizon 𝑁_𝑐 based on the model of (𝐴_𝑑,𝐵_𝑑,𝐶_𝑑)𝑋=𝐴Δ𝑉+𝛾, (14) where 𝑋=𝑥(𝑘)𝑥(𝑘+1∣𝑘)⋯𝑥𝑘+𝑁_𝑐−1∣𝑘^𝑇,𝐵𝐀=00⋯0_𝑑𝐴0⋯0_𝑑+𝐼𝐵_𝑑𝐵_𝑑⋯0⋮⋮⋮⋮_{𝑁𝑐−1𝑖=1}𝐴_𝑑^𝑖−1𝐵_𝑑^𝑁𝑐−2_𝑖=1𝐴_𝑑^𝑖−1𝐵_𝑑,𝑥(𝑘)𝑥(𝑘)⋯0𝛾=+𝐴_𝑑(𝑘)𝑥(𝑘)+Δ𝑥_2𝑖=1𝐴_𝑑^𝑖(𝑘)⋮𝑥(𝑘)+Δ𝑥_{𝑁𝑐−1𝑖=1}𝐴_𝑑^𝑖(𝑘).Δ𝑥 (15) The object function to be minimised is a quadratic criterion on Δ𝑉 as 𝐽=𝑅_𝑠−𝑌^𝑇𝑅_𝑠−𝑌+Δ𝑉^𝑇𝑅Δ𝑉. (16) By using (12), the optimization problem can be attributed to 𝐽_min=min_Δ𝑉12Δ𝑉^𝑇𝐻Δ𝑉+Γ^𝑇Δ𝑉subjecttoΔ𝑉≤Δ𝑉≤Δ𝑉, (17) where 𝐻=Ψ^𝑇Ψ+𝑅, Γ=Ψ^𝑇(−𝑅_𝑠+𝐹𝑥_𝑢(𝑘)), Δ𝑉 and Δ𝑉 indicates the minimal and maximum value of Δ𝑉.

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) where 𝑈=[𝑢_𝑘⋯𝑢_{𝑘+𝑁𝑐−1}]^𝑇 represents the control inputs, and 𝑈 and 𝑈 indicate the minimal and maximum value of 𝑈.

After the realization of the IOFL, the input vector 𝑈 is transformed into a new one 𝑉. 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 vectors 𝑉 and 𝑈 and then get the explicit expression of 𝑈 from 𝑉.

To establish the nonlinear relationship between vectors 𝑉 and 𝑈, expand the feedback law (4) over the control horizon 𝑁_𝑐, 𝑣(𝑘)(𝑥(𝑘𝑢(𝑘)(𝑥(𝑘,𝑣()(𝑥(𝑢()(𝑥(,⋮𝑣=𝑎))+𝑏))𝑘+1=𝑎𝑘+1))𝑘+1+𝑏𝑘+1))𝑘+𝑁_𝑐𝑥−1=𝑎𝑘+𝑁_𝑐−1𝑢𝑘+𝑁_𝑐−1+𝑏𝑥𝑘+𝑁_𝑐−1. (19)

Considering the expression (14), (19) can be rearranged in the matrix form 𝑉=𝐺_𝑢[𝑈,Δ𝑉]. (20) The above inequality can be rewritten in the following form: 𝐺[𝑈,Δ𝑉]=0 (21) due to 𝑣(𝑘)𝑣()⋮𝑣𝑉=𝑘+1𝑘+𝑁_𝑐=⋅(𝑘)()⋮−11011111⋮⋱⋱⋱11⋯11Δ𝑣Δ𝑣𝑘+1Δ𝑣𝑘+𝑁_𝑐().−1+𝑣𝑘−1 (22) In general, 𝑣(𝑘+𝑖−1) can be rewritten as 𝑣()()+𝑘+𝑖−1=𝑣𝑘−1^𝑖_𝑗=1Δ𝑣𝑘+𝑗−1𝑖=1,…,𝑁_𝑐, (23) with the constraints 𝑣≤𝑣≤𝑣, (24) where 𝑣, 𝑣 indicate the minimal and maximum value of 𝑣, which are state dependent.

Thus 𝑣_{𝑘+𝑖−1}(𝑥(()≤𝑘+𝑖−1))−𝑣𝑘−1_𝑖𝑗=1≤Δ𝑣𝑘+𝑗−1𝑣_{𝑘+𝑖−1}(𝑥(𝑘+𝑖−1))−𝑣(𝑘−1). (25) The above inequality can be rewritten over the entire horizon, Λ^𝑇Δ𝑉(𝑘)≤𝑐(𝑋(Δ𝑉(𝑘)))^𝑇, (26) whereΛ=𝐿^𝑇−𝐿^𝑇⋮,𝐿=,101111⋯1𝑐=𝑣_𝑘−𝑣(𝑘−1)⋯𝑣_{𝑘+𝑁𝑐−1}−𝑣(𝑘−1)𝑣(𝑘−1)−𝑣_𝑘⋯𝑣(𝑘−1)−𝑣_{𝑘+𝑁𝑐−1}. (27)

From (26), the above equation can be rewritten as Λ^𝑇Δ𝑉(𝑘)≤𝑐(𝑋(Δ𝑉(𝑈(𝑘))))^𝑇. (28)

This iterative QP algorithm is presented below: Step  1: Initializing 𝑐(𝑋(Δ𝑉(𝑈⁰(𝑘)))) within the input constraint. Step  2: Solve QP for Δ𝑉^𝑖(𝑘) subject to Λ^𝑇Δ𝑉^𝑖(𝑘)≤𝑐(𝑋(Δ𝑉(𝑈^𝑖−1(𝑘))))^𝑇 Step  3: Calculate 𝑐(𝑋(Δ𝑉(𝑈^𝑖(𝑘)))) Step  4: Test if Λ^𝑇Δ𝑉^𝑖(𝑘)≤𝑐(𝑋(Δ𝑉(𝑈^𝑖(𝑘))))^𝑇, end; otherwise, 𝑖=𝑖+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 about 𝑈₀ of the nonlinear constraints (21) as Δ𝑉=Δ𝑉₀+𝑔𝑈₀ℎ+𝑟ℎ, (29) where 𝑈₀ is the chosen initial operating trajectory, Δ𝑉₀ is obtained from (21) when 𝑈=𝑈₀, the (𝑁_𝑐×𝑁_𝑐) dimensional matrix 𝑔[𝑈₀] represents the Jacobian matrix 𝜕Δ𝑉/𝜕𝑈 in (21) at the operating trajectory, and 𝑟(ℎ) corresponds to the higher order terms of the approximation given by 𝑟ℎ=12ℎ^𝑇𝑔𝑈₀+𝜃ℎℎ, (30) where 𝜃∈[0,1], ℎ=𝑈−𝑈₀.

Neglect the higher order terms as Δ𝑉=Δ𝑉₀+𝑔𝑈₀𝑈−𝑈₀. (31) Equation (31) can lead to an explicitly linear function relationship between the control sequence 𝑈 and Δ𝑉: 𝑈=𝑔⁻¹𝑈₀Δ𝑉+𝑈₀−𝑔⁻¹𝑈₀Δ𝑉₀. (32) Now the optimization problem is defined as 𝐽_min=min_Δ𝑉12Δ𝑉^𝑇𝐻Δ𝑉+Γ^𝑇Δ𝑉. (33) subject to the convex set of approximate linear constraints 𝑈≤𝑀Δ𝑉+𝑚₀≤𝑈, (34) with 𝑀=𝑔⁻¹[𝑈₀], 𝑚₀=𝑈₀−𝑔⁻¹[𝑈₀]×Δ𝑉₀.

In order for the algorithm to reach the convergence to a feasible solution, the optimal problem can be rewritten as Δ𝑉^∗=argmin_Δ𝑉12Δ𝑉^𝑇𝐻Δ𝑉+Γ^𝑇Δ𝑉subjectto𝑈≤𝑀^𝛼Δ𝑉+𝑚₀^𝛼≤𝑈, (35)(36) where 0<𝛼<1 is decreasing parameter with initial value 𝛼=1. Define 𝑀^𝛼 and 𝑚^𝛼₀, respectively, by 𝑀^𝛼=𝑔⁻¹𝑈₀^𝛼𝐼−𝑔⁻¹𝑈₀^𝛼−1,𝑚^𝛼₀=𝑈₀−𝑔⁻¹𝑈₀^𝛼𝐼−𝑔⁻¹𝑈₀^𝛼−1Δ𝑉₀. (37) Define 𝑈^∗_𝑙=𝑀^𝛼Δ𝑉^∗+𝑚^𝛼₀ to be the linear control sequence and 𝑈^∗_𝑛𝑙 to be the nonlinear control sequence obtained from 𝐺[𝑈,Δ𝑉^∗]=0. Then the major goal is to keep 𝑈≤𝑈^∗_𝑛𝑙≤𝑈. (38) Thus the optimization problem with the approximate constraints Δ𝑉^∗ 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: 𝑈≤𝑈_𝑖≤𝑈. Solve for Δ𝑉₀, 𝑀,𝑚₀ Step  2: 𝛼_𝑖−1=1 Step  3: actualize the squeezing factor 𝛼: 𝛼=𝛼×𝜆. (39) Step  4: new linearized constraints at 𝑈₀: 𝑈≤𝑀^𝛼Δ𝑉+𝑚^𝛼₀≤𝑈. Step  5: Solve QP for the new linearized constraints: Δ𝑉^∗_𝑖+1=argmin_Δ𝑉𝐽(Δ𝑉). (40) Step  6: test if 𝑈^∗_𝑛𝑙 violate the original constraint. If: 𝑈≤𝑈^∗_𝑛𝑙≤𝑈⇒𝐸𝑛𝑑. Step  7: let 𝑖=𝑖+1, go to 3.

3.2. The Convergence Property

With the object function (35), assume the practical input sequence constraints |𝑈|≤1. (41) Assume a given initial vector representing the operating trajectory, 𝑈₀, verifying the conditions 𝑈₀≤1−𝜀₁, (42) where 𝜀₁ is small positive scalar values, 𝜀₁>0.

For any given scalar 0<𝛼<1, the solution of the optimization problem (35) is subjected to 𝑈₀+𝑔⁻¹𝑈₀^𝛼𝐼−𝑔⁻¹𝑈₀^𝛼Δ𝑉^∗−Δ𝑉₀≤1. (43) which can be further written as Δ𝑉^∗≤𝐼−𝑔⁻¹𝑈₀^−𝛼𝑔𝑈₀^𝛼1−𝑈₀+Δ𝑉₀. (44) In the limit, as 𝛼→0, (𝑔⁻¹[𝑈₀])^−𝛼→𝐼, (𝑔[𝑈₀])^𝛼→𝐼, the solution of the optimization problem is given by lim_𝛼→0Δ𝑉^∗≤Δ𝑉₀. (45) Due to the smooth function 𝐺[𝑈,Δ𝑉]=0 and definition (42), this last result leads to 𝑈^∗=𝑈₀≤1−𝜀₁, (46) 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: 𝑖=0. Step  2: if (𝑁−𝑖<0), go to 8. Step  3: assure feasibility of the new operating trajectory: 𝑈≤𝑈_𝑖≤𝑈. Step  4: solve QP at 𝑈_𝑖: Δ𝑉^∗_𝑖+1=argmin_Δ𝑉𝐽(Δ𝑉)subjecttoΛ^𝑇Δ𝑉^∗_𝑖+1≤𝑐𝑋Δ𝑉𝑈_𝑖^𝑇. (47)(48) Step  5: test if Δ𝑉^∗_𝑛𝑙 violates the constraints: if Λ^𝑇Δ𝑉≤𝑐(𝑋(Δ𝑉(𝑈_𝑖+1)))^𝑇⇒𝐸𝑁𝐷   where 𝑈_𝑖+1=𝑈^∗_𝑛𝑙 is nonlinear control sequence, obtained from 𝐺[𝑈^∗_𝑛𝑙,Δ𝑉^∗_𝑖+1]=0. Step  6: 𝑖=𝑖+1. Go to 2 Step  7: 𝑈₀=𝑈_𝑖, Δ𝑉₀=Δ𝑉^∗_𝑁. Solve for 𝑀, 𝑚₀. Step  8: 𝛼=1. Step  9: actualize the squeezing factor 𝛼: 𝛼=𝛼×𝜆. (49) Step  10: new linearized constraints at 𝑈₀: 𝑈≤𝑀^𝛼Δ𝑉+𝑚₀^𝛼≤𝑈 Step  11: Solve QP for the new linearized constraints: Δ𝑉^∗_𝑖+1=argmin_Δ𝑉𝐽(Δ𝑉). (50) Step  12: test if 𝑈^∗_𝑛𝑙 violate the original constraint. If: 𝑈≤𝑈^∗_𝑛𝑙≤𝑈⇒𝐸𝑛𝑑 Step  13: let 𝑖=𝑖+1. Go to 9.

The total algorithm performance is quite related to the tuning of two parameters: 𝑁: 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 value 𝑈₀ 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 𝑘 can be expressed as 𝑈^∗_𝑛𝑙=[𝑢(0)𝑢(1)⋯𝑢(𝑁_𝑐−2)]^𝑇. Then, the initial value 𝑈₀ at instant 𝑘+1 should be chosen as  𝑈₀=[𝑢(1)⋯𝑢(𝑁_𝑐−2)0]^𝑇; for example, the first 𝑁_𝑐−1 factors are taken directly from the prediction of 𝑈^∗_𝑛𝑙.

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 𝑥̇‌₁=𝑥₂,𝑥̇₂=−3𝑥²₁𝑥₂−𝑥³₁sin𝑥₂+𝑢,𝑦=𝑥₁. (51)

The input constraint is −0.5≤𝑢≤0.5.

Using (4), the static feedback law: 𝑣=𝑢+(−3𝑥²₁𝑥₂−𝑥³₁sin𝑥₂)

The resulting feedback linearization linear system is 𝑥̇₁𝑥̇₂=01𝑥00𝑥¹₂+01𝑣,𝑦=𝑥₁𝑦=𝑥₁. (52)

Choose the predictive horizon 𝑁_𝑝=8, the control horizon 𝑁_𝑐=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. Table 1 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.

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]: 𝑑𝑥₁𝑑𝑡=−𝑥₁+𝐷_𝛼𝑥⋅exp²1+𝑥₂,/𝛾𝑑𝑥₂𝑑𝑡=−𝑥₂1+𝛽+𝐻⋅𝐷_𝛼⋅1−𝑥₁𝑥⋅exp²1+𝑥₂/𝛾+𝛽𝑢,𝑦=𝑥₁, (53) where the two state variables 𝑥₁ and 𝑥₂ are the reactor temperature and the normalized reactant concentration. The control variable 𝑢 is the normalized cooling water temperature. 𝐷_𝛼, 𝛽, 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) where 𝑎=𝛽𝜉₁+𝜉₂⋅11−𝑟𝜉ln¹+𝜉₂𝐷_𝛼²,𝑏=𝜉₁−𝜉₁+𝜉₂+𝜉₁+𝜉₂⋅ln𝜉₁+𝜉₂/𝐷_𝛼1−(1/𝑟)ln𝜉₁+𝜉₂/𝐷_𝛼1+𝛽+𝐻⋅1−𝜉₁𝜉₁+𝜉₂. (55) The resulting feedback linearization linear system is 𝜉̇₁𝜉̇₂=01𝜉00𝜉¹₂+01𝜉𝑣,𝑦=10¹𝜉₂. (56)

Suppose the constraint to be −3≤𝑢≤3. Let 𝑟,𝐷_𝛼,𝐻=1, 𝛽=0.01, the predictive horizon 𝑁_𝑝=10, the control horizon 𝑁_𝑐=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 variable 𝑢 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 for 𝑁_𝑐=5 are given in Figure 4, adopting the best 𝑁 and 𝜆 for each 𝑁_𝑝 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. Table 2 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 𝑁_𝑝−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.

Figure 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.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grants 60974051, 61273144, and 61203041, the Natural Science Foundation of Beijing under Grant 4122071, the Chinese National Postdoctor Science Foundation under Grants 2011M500217 and 2012T50036, the Doctoral Fund of Ministry of Education of China under Grant 20120036120013, and the Fundamental Research Funds for the Central Universities 13XS09.