Terminal region characterization and stability analysis of discrete time quasi-infinite horizon nonlinear model predictive control

Journal of Process Control 83 (2019) 30–52

Contents lists available at ScienceDirect

Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont

Terminal region characterization and stability analysis of discrete time quasi-infinite horizon nonlinear model predictive control Chinmay Rajhans a , Devin W. Griffith b , Sachin C. Patwardhan c,∗ , Lorenz T. Biegler b , Harish K. Pillai a a

Department of Electrical Engineering, Indian Institute of Technology Bombay, Mumbai, India Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, USA c Department of Chemical Engineering, Indian Institute of Technology Bombay, Mumbai, India b

a r t i c l e

i n f o

Article history: Received 31 January 2019 Received in revised form 6 August 2019 Accepted 8 August 2019 Keywords: Nonlinear model predictive control Quasi-infinite horizon formulation Terminal region characterization Nominal stability Input to state stability

a b s t r a c t A terminal penalty term and terminal constraints were introduced in the nonlinear MPC (NMPC) formulation for establishing the nominal closed loop stability. A large terminal set gives rise to a large region of attraction for the closed loop system and can also help in reducing the on-line computation cost. In this work, two novel approaches have been developed for characterization of the terminal region for a discrete time quasi-infinite horizon NMPC formulation that employs quadratic stage cost. Similar to the continuous time QIH NMPC formulation by Chen and Allgöwer [1], a stabilizing linear controller developed using linearization at the origin is employed for characterization of the terminal set. The first approach permits use of an arbitrary linear controller while the second approach is based on LQR controller. Using a method of bounding only the higher order nonlinear effects of the system via simulations under linear controller, the proposed approach is further extended to handle terminal region computation for a large dimensional system. Unlike the approaches available in the literature, the proposed approaches provide sufficient degrees of freedom to shape the terminal set. The quadratic control Lyapunov function in the terminal set together with the quadratic stage cost is further used to establish exponential stability of the discrete time QIH NMPC formulation under nominal conditions. Additionally, it is shown that discrete time QIH-NMPC based on the nominal model is input to state stable in the face of bounded disturbances in the state dynamics. Efficacy of the proposed approaches for characterization of the terminal region is demonstrated using a two state example from Chen and Allgöwer [1] and a benchmark CSTR system. Moreover the extension to large dimensional systems is demonstrated using a system consisting of two distillation columns in series. Parametric studies reveal that the choice of sampling interval and the choices of the tuning matrices in both approaches have significant influence on the sizes of terminal sets. Moreover, the simulation results demonstrate that the systems under consideration converge asymptotically to the origin under the nominal conditions and are input to state stable in the face of bounded disturbance. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Nonlinear model predictive control (NMPC) schemes facilitate systematic handling of system nonlinearity, operating constraints and multivariable interactions. As a result, it is increasingly being employed in many industrial applications. Over the last three decades, a wide spectrum of issues associated with NMPC, ranging from fast real time computations to integration of economic objectives, have been investigated by the academic community [2,3]. In

∗ Corresponding author. E-mail address: [email protected] (S.C. Patwardhan). https://doi.org/10.1016/j.jprocont.2019.08.002 0959-1524/© 2019 Elsevier Ltd. All rights reserved.

particular, significant progress has been made on establishing the nominal and robust stability of nonlinear MPC schemes [3–5]. In a seminal contribution, Michalska and Mayne [6] developed a dual mode receding horizon formulation in which the concept of terminal set is introduced. The NMPC is expected to steer the system states into the terminal set over a finite time horizon in the future. The local linear controller would then move the system to the origin. The terminal set was characterized using a linear controller designed using local linearization at the origin. Taking motivation from this dual mode formulation, Chen and Allgöwer [1] developed aquasi-infinite horizon NMPC (QIH-NMPC) formulation. To ensure nominal stability with this approach, a terminal state constraint is introduced along with the terminal cost in the NMPC

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

formulation. The linear controller was used only for characterization of the terminal set and for construction of a control Lyapunov function in the terminal set. The set of all initial states from which QIH-NMPC can drive system states inside the terminal set determines the region of attraction for the closed loop system. The size of region of attraction is determined by choice of the prediction horizon and size of the terminal set. Subsequently, a significant fraction of the literature on NMPC has used inclusion of the terminal cost and terminal constraint for ensuring nominal stability [3–5]. If the system under consideration is open loop stable, then it is possible to do away with the terminal state constraint at the cost of significant increase in the prediction horizon [7]. This, however, can result in substantial increase in on-line computations. In a recent review article, Mayne [3] observes that “an explicit rather than an implicit terminal constraint yields a larger region of attraction”. A large terminal set gives rise to a large region of attraction for the closed loop system. Despite this attractive feature, the methods for characterization of the terminal region have not received much attention in the literature. Mayne [3] also observes that the state and terminal constraints are avoided (in industrial applications) possibly due to extra recovery measures needed to cope with violation of these constraints. Griffith et al. [8] have recently shown that it is possible to avoid such constraint transgressions by use of adaptive prediction horizon length. In fact, they show that it is also possible to achieve significant reduction in on-line computation time if their horizon length adaptation strategy is employed along with inclusion of the terminal constraint. The approach developed by Chen and Allgöwer [1] was for an NMPC formulation that employs a continuous time nonlinear dynamic model for prediction (i.e., a continuous time QIH-NMPC). Since real time implementation of any NMPC scheme involves digital controls discretized versions of prediction models have to be employed in NMPC realization. It is well established that discrete time control system requires distinctly different considerations for the closed loop stability analysis and controller synthesis [9]. The stability results that are established for the continuous time QIH-NMPC formulation may hold if the sampling time is chosen sufficiently small. However, this requirement on sampling time can lead to a large prediction horizon and, as a consequence, a significant increase in the online computation cost. Since it is advantageous to choose a relatively larger sampling time and smaller prediction horizon to make NMPC computations feasible in a real-time implementation, terminal region characterization methods that are tailored for discrete time QIH-NMPC formulation have been subsequently proposed in the literature. The approach presented by Rawlings et al. [10] makes use of tuning matrices in NMPC to develop a linear quadratic regulator using model linearization at the origin, and use this LQR for characterization of the terminal region. This scheme, however, does not provide any additional degrees of freedom to shape the region. The approach developed by Limon [11] allows use of any arbitrary stabilizing controller for characterization of the terminal set. A control Lyapunov equation is constructed using stabilized linearized system. He introduces an additional tuning parameter that can be varied to shape the terminal set. The approaches developed by Johansen [12] and recently by Yu et al. [13], on the other hand, can be viewed as a discrete time extension of the characterization methods developed by Chen and Allgöwer [1]. These approaches also provides only a single tuning parameter for shaping the terminal set. Also, these methods require finding a Lipschitz constant for the nonlinear part of the system or solving a series of nonconvex optimization problems to global optimality, either of which makes application to a large system cumbersome. Majority of available approaches for terminal region characterization make use of a linear controller in the neighborhood of the origin. Lucia et al. [14] have shown that the size of the terminal

31

region can be increased significantly if a nonlinear controller is used for characterizing the terminal set. The Taylor series expansions of the system dynamics and the stage cost in the neighborhood of the origin are used for constructing a nonlinear controller and characterize the terminal set. This approach, however, has been developed for continuous time NMPC and for a special class of nonlinear systems for which time derivatives are polynomial functions. Recently, Lazar and Tetteroo [15] have developed a terminal region characterization approach for discrete time NMPC that can use a linear as well as a nonlinear controller for terminal set characterization. A novel feature of their approach is use of periodically time varying terminal sets, terminal costs and terminal control laws, which facilitates enlargement of the domain of attraction of the closed loop system. Since characterization based on local linearization of the system dynamics can result in conservative estimates of the terminal sets, they propose to use second order Taylor series approximation of the system dynamics for enlarging the terminal sets. For discrete time nonlinear systems with input affine structure, they propose to enhance the region of attraction by employing a nonlinear terminal control law. It is important to note that this approach provides sufficient degrees of freedom for enlarging the domain of attraction. In the literature on NMPC formulations with guaranteed stability, work by Chen and Allgöwer [1] stands out because it establishes nominal stability for the most commonly used NMPC formulation (i.e. based on the quadratic stage cost) and it takes ‘constructive’ approach rather than merely assuming existence of a control Lyapunov function in the terminal set. The method developed by Chen and Allgöwer [1] for construction of the control Lyapunov function and explicit characterization of the terminal set enables them to prove stronger stability results. However, this approach for construction of the terminal region and its discrete time variants that have appeared in the literature subsequently provide only a single degree of freedom for shaping the terminal set. While designing NMPC for a large dimensional system, a single tuning parameter may prove inadequate for shaping the terminal set. Thus, it is necessary to develop methods that provide sufficient degrees of freedom to compute the terminal cost and to characterize the terminal region for discrete time QIH-NMPC. Moreover, the proposed approach needs to be scalable to large dimensional systems. Recently, Rajhans et al. [16] have proposed a terminal region characterization method for discrete time QIH-NMPC formulation, in which the tuning knob consists of a symmetric positive definite matrix, which provides significant flexibility to shape the terminal set. In this work, we built upon [16] and develop a new approach for terminal region characterization which is based on a linear quadratic regulator (LQR) designed in the neighborhood of the origin. Further, we propose a method of bounding only the higher order nonlinear effects of the system via simulations. This leads to a method of calculating terminal region for a large dimensional system. By assuming that the set of feasible states is compact, we establish exponential stability of the nominal discrete time QIHNMPC formulation that employs quadratic stage cost. We then proceed to show that the nominal controller is input to state stable (ISS) in the presence of bounded disturbances in the state dynamics. Efficacy of the proposed approaches for terminal region computation is demonstrated using a two state system from [1], a benchmark CSTR system [17] and a large scale system consisting of two distillation columns in series [18]. This paper is organized into six sections. Section 2 briefly describes the discrete time quasi-infinite horizon NMPC formulation. Section 3 presents approaches for the characterization of the terminal region and its computationally tractable extension to large dimensional systems. The nominal stability and robust stability of the discrete time QIH-NMPC formulation is established in Section 4. Three simulation case studies that demonstrate efficacy of the pro-

32

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

posed approach are discussed in Section 5 and main conclusions reached through analysis of the results are summarized in Section 6. 1.1. Notations and definitions √ For a vector z, the two norm is defined as |z|:= zT z. The induced norm of a matrix A is given as |A| = supz |Az| / 0. Weighted |z| for z = √ norm is defined as |z|A := zT Az. min (A) and max (A) denote minimum and maximum eigenvalues   of matrix A respectively. For a given sequence d(k) for k ≥ 0, d  sup|d(k)|.

A5 States {x(k) : k = 0, 1, 2, . . . } are perfectly measured. A6 There are no disturbancesin the system.  At instant k, let U[k,k+N−1] ≡ u(k), u(k + 1), . . ., u(k + N − 1) denote a set of future manipulated inputs. Then, for the system given by Eq. (2), a discrete time version of the quasi-infinite horizon nonlinear model predictive control (QIH-NMPC) [1] is formulated as follows: arg min J(x(k), U[k,k+N−1] ) where

k≥0



k+N−1

Consider an autonomous discrete time system x(k + 1) = ˜f(x(k))

(3)

U[k,k+N−1]

J(x(k), U[k,k+N−1] ) = (1)

l(z(i), u(i)) + Vf (z(k + N))

i=k



T



T

(4)

together with a set X ⊂ Rnx such that the origin is inside X. Further assume that ˜f( · ) is locally Lipschitz and x = 0 is an equilibrium point of the system. Definitions of asymptotically stable equilibrium point, class K function, class KL function are given as follows:

l(z(i), u(i)) = z(i) Wx z(i) + u(i) Wu u(i)

Definition 1.

z(i + 1) = F(z(i), u(i)) for k ≤ i ≤ k + N − 1

(7)

z(k) = x(k)

(8)

z(i) ∈ X for k ≤ i ≤ k + N

(9)

[5]

The equilibrium point, 0, of the autonomous discrete time system (1) is asymptotically stable if lim x(k) = 0 for all x(0) ∈ X with X being the region of attraction. Definition 2.

k→∞

[5,8]

A function (s) : R≥0 → R≥0 belongs to class K if it is continuous, and strictly increasing with s and (0) = 0. It belongs to class K∞ if it is class K and unbounded (i.e. (r)→ ∞ as r→ ∞). A function ˇ : R≥0 × Z≥0 → R≥0 is of class KL if, for each t ≥ 0, ˇ(· , t) is a class K function, and, for each s ≥ 0, ˇ(s, ·) is non-increasing and lim ˇ(s, t) = 0. t→∞

Definition 3.

[5]

A function V : Rnx → R≥0 is said to be a control Lyapunov function for the system x(k + 1) = ˜f(x(k)) if there exist class K∞ functions ˛1 (·) and ˛2 (·) and positive definite function ˛3 (·) such that ∀x(k) ∈ X, function V(x(k)) has the following properties: V (x(k)) ≥ ˛1 (|x(k)|)

Vf (z(i)) = z(i) Pz(i)

u(i) ∈ U for k ≤ i ≤ k + N − 1

(10)

z(k + N) ∈ 

(11)

where N is the prediction and control horizon length, l( · , · ) denotes the stage cost term, Vf (·) denotes the terminal cost function, Wx and Wu denote (nx × nx ) and (nu × nu ) symmetric positive definite weighting matrices, and P is the symmetric positive definite terminal weighting matrix  of the discrete time NMPC formulation. Here, z(i) : k < i ≤ N denotes the states predicted by the internal model used in the NMPC formulation and XN ⊆ X represents a set of all initial conditions {x(k) : k = 0, 1, . . . } such that z(k + N) ∈ . The set, , represents the terminal region in the neighborhood of the origin and is chosen such that the set is invariant for the nonlinear control system controlled by a fictitious local linear state feedback controller with gain matrix say, L. The quadratic terminal cost term is chosen such that the following holds ∀z(k + N) ∈ : Vf (z(k + N)) ≥

V (F(x(k))) ≤ V (x(k)) − ˛3 (|x(k)|)

∞ 

l(z(i), u(i))

(12)

i=k+N

u(i) = −Lz(i) ∈ U,

[5]

An autonomous discrete time system x(k + 1) = ˜f(x(k)) is locally asymptotically stable, if there exists a control Lyapunov function on region of attraction XN . 2. Discrete time quasi-infinite horizon NMPC formulation

It 

is

assumed

∗

∗

(k + N) ≤ i < ∞ that ∗

∗

(2)

where x(k) ∈ X ⊂ Rnx represents the system states and u(k) ∈ U ⊂ Rnu represent the vector of manipulated inputs, together with initial condition x0 . It is assumed that: A1 F : X × U → X is continuously differentiable over X ⊂ Rnx and U ⊆ Rnu . A2 Set U is compact and contains the origin in the interior. A3 The set X ⊂ Rnx is control positive invariant for F(· , ·), i.e. F(x(k), u(k)) ∈ X for all x(k) ∈ X and all u(k) ∈ U. Moreover, set X is compact and contains the origin in its interior. A4 Origin is an equilibrium state of the system, i.e. F(0, 0) = 0.

an

optimal



(13) solution,

∗ U[k,k+N−1]

≡

u (k), u (k + 1), . . ., u (k + N − 1) , exists and can be reached through a suitable numerical approach. The controller is implemented in the moving horizon framework, i.e. only the first optimal input move

u(k) = u (k)

Consider a discrete time nonlinear system as given by x(k + 1) = F(x(k), u(k))

(6)

subject to:

V (x(k)) ≤ ˛2 (|x(k)|)

Theorem 4.

(5)

T

(14)

is implemented in the system. At instant k + 1, the NMPC optimization is re-formulated over [k + 1, k + N + 1] and solved again. The term quasi-infinite horizon is employed because the cost function is formulated such that it bounds the infinite horizon cost function. The optimal manipulated input, u(k), is always computed by solving the discrete quasi-infinite horizon NMPC formulation, whether x(k) is inside or outside . The local feedback controller is called fictitious because it is used only for characterization of the terminal region, , and never implemented. To see how this is done, consider  theJacobian linearization of (2) in the neighborhood of the origin, 0, 0 , x(k + 1) = x(k) + u(k)

(15)

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

where







approach which is based on an LQR design and we further propose a method of bounding only the higher order nonlinear effects of the system via simulations under LQR control. This leads to a method of calculating terminal region for a large dimensional systems.



∂F ∂F and  = ∂x 0,0 ∂u 0,0

˚=

It is further assumed that A7 The linearized system (15) is controllable. Under this assumption, a stabilizing feedback control law, u(k) = − Lx(k), can be found using any of the available design approaches such that matrix c =  − L has eigenvalues strictly inside the unit circle. To characterize the terminal region, consider a scenario where u(k) = − Lx(k) whenever x(k) ∈ . The dynamics of the system (2), when restricted to set, , can be expressed as follows

3.1. Arbitrary stabilizing controller based approach Lemma 6. Suppose that Assumptions A1 to A7 hold and a stabilizing state feedback control law has been designed such that c = ( − L) has all eigenvalues inside the unit circle. Let matrix P denote the solution of the following Lyapunov equation ˚c P˚c − P = −(Q∗ + Q) T

(21)

where Q is a positive definite matrix and Q* is defined as follows

x(k + 1) = F(x(k), −Lx(k)) = (˚ − L) x(k) + (x(k))

∗

Q :=Wx + LT Wu L

where  (x(k)) = F(x(k), −Lx(k)) − (˚ − L) x(k)

(16)

The key step is to find the controller gain matrix, L, and the associated matrix, P > 0, such that inequality (12) holds in . It is easy to show that |(x(k))| can be bounded from above by a class K function in the set . We prove the following Lemma, which is subsequently used for characterization of  for large scale systems. There exists M, q ∈ R+ such that

Lemma 5.

∀x(k) ∈ 

| (x(k))| ≤ M|x(k)|q ,

T

x(k) ∇ 2 j ( x(k))x(k) d ,

∀j = 1. . .nx

0

where nx is the number of states. Note that, at x(k) = 0, j (x) = 0 and ∇ j (x(k))T = 0. Then



1 T

x(k) ∇ 2 j ( x(k))x(k) d

(18)

0

Since F is twice continuously differentiable, we define:

m :=

max

x(k) ∈ ,j ∈ {1...nx }

|∇ 2 j (x(k))|

(19)



T

x(k) ∈ Rnx : x(k) Px(k) ≤ ˛, −Lx(k) ∈ U

(23)

in the neighborhood of the origin such that  is an invariant set for the nonlinear system (2) controlled by u(k) = − Lx(k). Moreover, for any x(k + N) ∈ , the infinite horizon cost for the nonlinear system (2) starting from (k + N) and controlled by the local linear controller satisfies the inequality ∞ 

l(x(i), u(i))

(24)

Proof. Since the eigenvalues of c = ( − L) are inside the unit circle, from the solvability condition of the discrete time Lyapunov equation, it follows that a unique P > 0 can be found that solves Eq. (21). According to Assumption A1, the origin 0 ∈ Rnu is in the interior of the set U. Then, we can always find a constant and a set  in the neighborhood of the origin 0 ∈ Rnx such that  ≡





T

x(k) ∈ Rnx : x(k) Px(k) ≤ , −Lx(k) ∈ U

|j (x(k))| ≤

m |x(k)|2 2

(26)

where the term, (x), is defined by Eq. (16). Consider a candidate control Lyapunov function defined as

which gives |(x(k))|2 =

(25)

Now, let ˛ ∈ (0, ] specify a region of the form given by Eq. (23). Since the input constraints are satisfied in , the system dynamics can be viewed as unconstrained in . Thus, inside  the nonlinear system dynamics together with the LQR can be expressed as follows x(k + 1) = ˚c x(k) +  (x(k))

and thus

T

nx 

|j (x(k))|2 ≤ nx



m

2

|x(k)|2

Vf (x(k)) = x(k) Px(k)

2

where P is the solution of Eq. (21). Note that

therefore √

nx

m |x(k)|2 2

2



2

min (P) x(k) ≤ Vf (x(k)) ≤ max (P) x(k)

j=1

|(x(k))| ≤



i=k+N

1

1 j (x(k)) = 2

≡

Vf (x(k + N)) ≥

j (x(k)) = j (0) + ∇ j (0)T x



(22)

where matrices (Wx , Wu ) define the stage cost of QIH-NMPC. Then, there exists a constant, ˛ > 0, specifying an ellipsoid of the form

(17)

By Taylor’s Theorem we have

Proof.

1 + 2

33

(20)

䊐

Using Eqs. (21) and (26), the difference of Lyapunov candidate Vf (x(k)) can be expressed as follows Vf (x(k))

= Vf (x(k + 1)) − Vf (x(k)) T

T

T

= x(k) (˚c P˚c − P)x(k) + 2 (x(k)) P˚c x(k)

3. Two approaches for the characterization of the terminal region In this work, we develop two approaches for characterization of the terminal region, . The first approach permits use of an arbitrary stabilizing linear controller (of fixed gain) for the terminal region characterization. For a large dimensional multi input multi output (MIMO) system, it is relatively easy to design a stabilizing controller using LQR approach. Thus, we develop another tailored

(27)

(28)

T

+ (x(k)) P (x(k)) Combining Eqs. (21) with Eq. (28), we have Vf (x(k))

= −x(k) (Q∗ + Q)x(k) + 2 (x(k)) Pc x(k) T

T

T

+ (x(k)) P (x(k)) = −x(k) Q∗ x(k) − ϕ(x(k)) T

(29)

34

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

where T

T

T

ϕ(x(k)) = x(k) (Q)x(k) − 2 (x(k)) P˚c x(k) −  (x(k)) P(x(k)) (30) Defining induced norm of the operator (x(k)) as ˇ ≡ sup{

|(x(k))| / 0} |x(k) ∈ , x(k) = |x(k)|

(31)

we can derive the following inequalities for the components of ϕ (x(k)) T

 (x(k)) Pc x(k) ≤ ˇ |P˚c |.|x(k)|2

(32)

T

 (x(k)) P (x(k)) ≤ ˇ2 |P|.|x(k)|2 T

(34)

Combining these inequalities with Eq. (29), it follows that Vf (x(k)) ≤ −x(k) Q∗ x(k) − [min (Q) − 2ˇ |P˚c | − |P|ˇ2 ]|x(k)|2 T

(35) Thus, the discrete system (2) controlled by u(k) = − Lx(k) is asymptotically stable in the neighborhood of the origin for all x(k) ∈  if the following condition is met 2

min (Q) − 2ˇ |P˚c | − |P|ˇ ≥ 0

2

|P|ˇ + 2ˇ |P˚c | − min (Q) = 0

−|P˚c | +



|P˚c |2 + min (Q)|P| |P|

(37)

Since | (0)| = 0, and fact that |(x(k))| is bounded from above by a class K function (ref: Lemma 5), we can choose ˛ > 0 such that ˇ ≤ ˇ∗ is in the set ˛ defined by Eq. (23), and then Vf (x(k)) < 0 in . Thus, the scalar function, Vf (x(k)), is a Lyapunov function for the nonlinear system when x(k) ∈  and the nonlinear system together with u(k) = − Lx(k) is asymptotically stable at the origin. Moreover,  is an invariant set for the nonlinear system (2). Now, the satisfaction of inequality ϕ (x(k)) ≥ 0 in  implies that Vf (x(k + 1)) − Vf (x(k)) ≤ −x(k) (Q∗ )x(k) = l(x(k), −Lx(k)) T

(38)

when x(k) ∈ . Summing the inequality (38) on both the sides over the range, [k + N, ∞), it follows that |x(k + N)|2P ≥

∞  



JLQ

x˜ , u[k:∞)

=

∞   i=k+N

 2

2

|x(i)| + |u(i)| ˜x ˜u W W

(39)

such that x˜ = x(k + N) belongs to the terminal region to be determined. Here the weighting matrices in the LQR formulation are ˜ x, W ˜ u ) such that W ˜ x > Wx > 0, W ˜ u ≥ Wu > 0 where chosen as (W (Wx , Wu ) represent the weighting matrices in the optimization for˜ x, W ˜ u ), the mulation given by Eqs. (3)–(10). Given a choice of (W steady state solution of the discrete-time algebraic Riccati equation (DARE) equations [9], i.e. T

˜ u LLQ ) ˜ x + LT W ˚c PLQ ˚c − PLQ = −(W LQ T

˜ u +  PLQ  ) LLQ = (W

−1

(40)

T

 PLQ ˚

(41)

yields P(= PLQ ) and the controller gain matrix L(= LLQ ). Under the ˜ x, W ˜ u > 0, a unique positive definite solution assumption A6 and W to DARE exists and c = ( − LLQ ) is Hurwitz [5]. Taking advantage of the fact that the LQR design equation itself is the discrete time Lyapunov equation, we state a tailored result that facilitates characterization of the terminal region. Lemma 8.

Suppose that Assumptions A1 to A7 hold. Let (PLQ , LLQ )

˜ x > Wx > represent the solution of DARE (40) and (41) for chosen W ˜ 0, Wu ≥ Wu > 0 such that ( − LLQ ) is Hurwitz. Then, there exists a constant, ˛ > 0, specifying an ellipsoid of the form ≡





T

x(k) ∈ Rnx : x(k) PLQ x(k) ≤ ˛, −LLQ x(k) ∈ U

(42)

in neighborhood of the origin such that  is an invariant set for the nonlinear system (2) controlled by the LQR. Moreover, for any x(k + N) ∈  the infinite horizon cost (39) for the nonlinear system (2) starting from (k + N) and controlled by the LQR satisfies inequality |x(k + N)|2P

≥

∞  

|x(i)|2Wx + |u(i)|2Wu



(43)

i=k+N

Proof. Steps in the proof are similar to the proof of Lemma 6 with P = PLQ and L = LLQ except for the following modification: Eq. (28) is changed to

i=k+N

i.e. inequality (24) holds for any x(k + N) ∈ . 䊐 Remark 7. Limon [11] has presented an approach similar to that of the arbitrary controller based approach for the terminal region characterization of the finite horizon NMPC using modified discrete time Lyapunov equation. This approach proposes to use Q = ( − 1) · Q* with   1 in the Lyapunov equation (21). Thus, the approach presented by Limon [11] can be viewed as a subclass of the proposed arbitrary controller based approach with a single tuning parameter. Recently Yu et al. [13] have proposed an approach for characterization of the terminal set which can be viewed as an extension of approach in Chen and Allgöwer  [1]. This approach also

employs a single tuning parameter ∈



LQ

x(i) (Q∗ )x(i) T

Consider an LQR designed for the linearized system (15) that minimizes the following unconstrained cost function

(36)

To compute the limiting value of ˇ , consider the positive solution of the quadratic equation

ˇ∗ =

3.2. LQR based approach

(33)

x(k) Qx(k) ≥ min (Q)|x(k)|2

i.e.

to larger . However, for a system with  (˚c ) very close to 1 (possibly due to fast sampling), the tuning parameter can get restricted in a very narrow range. On the other hand, the proposed approach permits selection of any Q > 0 and provides nx (nx + 1)/2 degrees of freedom for shaping the terminal set.

1,

1 (˚c )

for determin-

ing size of the terminal region, where  (˚c ) denotes the spectral radius c . Their derivation indicates that larger value of will lead

Vf (x(k))

T

T

= x(k) (˚c PLQ ˚c − PLQ )x(k) T

T

(44)

+ 2 (x(k)) PLQ ˚c x(k) +  (x(k)) PLQ  (x(k)) Defining matrices Q* and Q as follows Q∗ :=Wx + LTLQ Wu LLQ

(45)

Q:=Wx + LTLQ Wu LLQ

(46)

˜ x − Wx Wx :=W

(47)

˜ u − Wu and Wu :=W

Eq. (40) can be expressed as follows ˚c PLQ ˚c − PLQ = −(Q∗ + Q)

(48)

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

35

and Eq. (44) reduces to Vf (x(k)) = −x(k) (Q∗ )x(k) − ϕ (x(k)) T

(49)

Where T

T

(50)

Arguments in the rest of the proof are identical to the proof of Lemma 6. 䊐 In Section 3.1, in the arbitrary stabilizing controller based approach for characterizing , has a linear controller gain matrix L that is computed. On the other hand, in the LQR based approach, ˜ x, W ˜ u ). Thus, the gain matrix LLQ as well as Q are functions of (W one approach does not reduce to the other even though Eqs. (21) and (40) are qualitatively similar. Remark 9. Johansen [12] have presented an approach similar to LQR based approach for the characterization of the terminal region for discrete time QIH-NMPC. Their approach requires the linear controller to be LQR, which is designed using (Wx , Wu ) as the weighting matrices. To find the terminal penalty matrix, they propose to solve the following modified Lyapunov equation ˚c P˚c − P = −(Q∗ + P) T

(51)

where > 0. Similar to approach by Limon [11], these approaches also provides only a single degree of freedom for shaping of the terminal region. On the other hand, the method presented by Rawlings and Mayne [5] is equivalent to selecting = 0, in Eq. (51), which leaves no degree for freedom for shaping of the terminal region. Remark 10. Recently Lazar and Tetteroo [15] have proposed inclusion of time varying terminal cost terms of the form (j)

Vf (x) = xT P(j) x,

P(j) > 0, j = 0, 1, . . ., M − 1; M ∈ N≥1

(52)

for guaranteeing nominal stability of the discrete time QIHNMPC with quadratic stage cost. This approach involves design of (j) a finite sequence linear terminal control  laws, u = − L x, and  associated ellipsoidal terminal sets (j) = xT P(j) x ≤˛, ˛ > 0 , which are required to satisfy of the following coupled inequalities

T

(j+1)

P

(j) ˚c

(j)

− (1 − j )P



+ Wx + L



(j) T

(j)

Wu L

(j)

≤0

(j)

for j = 0, 1, . . ., M − 1 where ˚c = ( −  L ) and 

(53) (j)

(x) =

(j)

F(x, −L(j) x) − ˚c x. The set of inequalities (53) are solved for finding (P(j) , L(j) ) by transforming them to LMIs and the terminal sets are characterized by solving for maxRj (x) ≤ 0

(54)

x ∈ j

where Rj (x) =

 

(j)

− j xT P(j) x

T (x)

P(j+1) 

(j)



(x) + 2 



T j

˚c

[Q∗ + Q] (˚c )

j

j=0

˚c x(k) −  (x(k)) PLQ  (x(k))

(j) ˚c

∞  

T

ϕ (x(k)) = x(k) (Q)x(k) − 2 (x(k)) PLQ



P=

(j)

T (x)

(j)

P(j+1) ˚c x (55)

and j = 0, 1, . . ., M − 1. Unlike most approaches available in the literature, this approach provides several tuning for  parameters  increasing the region of attraction, viz., M and L(j) , j for j = 0, 1, . . ., M − 1. To reduce conservatism of region of attraction, use of quasi-second order Taylor series approximation is proposed in place of the first order approximation while computing (j) (x). For discrete time nonlinear systems with input affine structure, they propose to enlarge the region of attraction by employing nonlinear terminal control law, which is similar to [14]. Remark 11. Since matrix c =  − L is asymptotically stable by construction

is the unique symmetric and positive definite solution of the discrete Lyapunov equation. Each term in the series solution has computational complexity of O(n3x ) where nx is the state dimension. Thus, for large dimensional systems, computing series solution can prove computationally efficient. However, this is a linear problem with N = nx (nx + 1)/2 unknowns and the complexity of solving the resulting problem is O(N3 ). Thus, for moderately large dimensional systems, exact solution of the discrete time Lyapunov equation can be computed relatively easily. 3.3. Computation of the terminal region Lemmas 6 and 8 give a method for explicit characterization of the terminal region, , as a set of all x(k) such that ϕ (x(k)) ≥ 0. Also, (x(k)) → 0 and |(x(k))|/|x(k)| → 0 as |x(k)| → 0. Instead of using norm based method, one can directly inequality ϕ (x(k)) ≥ 0 for the characterization of the terminal region. Accordingly, if we choose ˛ > 0 such that ϕ (x(k)) ≥ 0 in the set, , defined by Eq. (23), then Vf (x(k)) < 0 in . Thus,  can be characterized as the largest value of ˛ such that ≡



T

x(k) : x(k) ∈ Rnx , x(k) Px(k) ≤ ˛,



−Lx(k) ∈ U, ϕ (x(k)) ≥ 0

(56)

In this work, the terminal set is characterized as defined by (56). The steps involved in computation of  are as follows: Step 1: Computation of upper bound ␥: Compute the largest such that T

 ≡ {x(k) : x(k) ∈ Rnx , x(k) Px(k) ≤ ,

− Lx(k) ∈ U}

(57)

Step 2: Characterization of Terminal Region: Find the largest possible ˛ ∈ (0, ] such that



global min x(k) ∈ 



ϕ (x(k)) = 0



(58) T



where  = x(k) ∈ Rnx : x(k) Px(k) ≤ ˛ , and ϕ (x(k)) is defined by Eq. (30) or (50). This can be achieved by solving a series of optimization problems for different values of ˛ ≤ . When system dynamics exhibits severe nonlinearity, it is likely that a conventional NLP solver can get trapped in a local minimum during iterations. To avoid such situations, we have developed a method to characterize the terminal region based on random sampling. This method is used to generate good initial guesses for application of the conventional optimization based search. A combination of the random sampling based approach and a conventional NLP solver is used to locate the global minimum. Details of this method are discussed in Appendix A. To compare terminal regions computed using different approaches, it becomes necessary to quantify size of a terminal region. One of the ways of quantifying the terminal region is using the area or volume or hyper volume enclosed inside the closed region. Since P is symmetric and positive definite matrix, one can diagonalize it as P = T where the matrix  consists of eigenvalues of P and  is an orthonormal matrix consisting of eigenvectors of P as columns. Thus, defining a rotational transformation, y(k) = T x(k), one can write T

T

x(k) Px(k) = y(k) y(k) ≤ ˛

36

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

The hyper-volume of the region enclosed in y(k)T y(k) ≤ ˛ is given as [19] HV =

cnx ˛nx /2



1 . . .nx

=

cnx ˛nx /2



det(P)

where 1 , 2 , . . ., nx represent eigenvalues of P and cnx is a constant dependent only on the state dimension nx , and is given by Eq. (59). cnx =

nx /2



 nx 2



(59)

+1

where  is a gamma function [19]. Values of cnx include c2 = , c3 = 4/3, c4 = 2 /2, c5 = 82 /15, c6 = 3 /6. We also note that solving the globally optimization problem (58) and the RSG approximation are NP hard and can become very expensive for systems with large nx . For this reason, we consider a further approach for large systems. Remark 12. In the case of linear time invariant systems, (x(k)) = 0 for all x(k) and the descent property of the control Lyapunov function holds for any choice of Q. Thus, the terminal region can be characterized using only the input constraints, i.e.  =  . Shape of the terminal set, however, depends on P matrix, which, in turn, depend on choices of L and Q. While the descent property of the control Lyapunov function holds even for Q = [0], from the viewpoint of reducing on-line computations for systems with very fast dynamics, it is desirable to find a suitable combination of L and Q that maximizes the size of the terminal region [20]). For chosen L and Q, characterization of  this case reduces to solving a sequence of QP problems. 3.4. Extension for large scale systems As motivation for a further extension to this methodology, we note that the Lipschitz constant ˇ defined in (31) must be valid in the entirety of the terminal region , however if ˇ is checked in a conservatively large region to guarantee that (31) holds, then a conservatively large value of ˇ may be found, which can lead to a small . Thus, to avoid conservatively small terminal regions, cumbersome iterations (calculating ˇ , finding , and then checking to see if ˇ is valid in ) may be required. Furthermore, solving (36) to global optimality may be prohibitively difficult for a large system. Therefore, we propose an approach that scales more easily to large systems. From (20), we see that (17) is satisfied with √ M = nx 2m and q = 2. Note that, in general, the system dynamics (2) are obtained from an implicit discretization of a set of differential and algebraic equations (DAEs) so that  cannot be obtained explicitly, and actually quantifying (17) may be tedious. Instead, for practical purposes we will find M and q in (17) explicitly via simulations from a sampling of initial conditions in the state space, as will be shown in Section 5.3. Here the key advantages of this method are apparent, in that we only need to solve a series of one step simulations using the linear control, and do not need to iterate on regions for which a Lipschitz constant is valid. It is also possible, however, to obtain larger terminal regions, since a power law bound on a nonlinear term is expected to be a better fit than a linear bound. Given this form of bound, in place of (32) and (33) we have T

 (x(k)) Pc x(k) ≤ M|Pc ||x(k)|q+1 T

2

2q

 (x(k)) P(x(k)) ≤ M |P||x(k)|

(60) (61)

And in place of (36) we have: min (Q)|x(k)|2 − 2M|Pc ||x(k)|q+1 − M 2 |P||x(k)|2q ≥ 0

(62)

Then since (17) holds for all x ∈ , we can solve (62) for |x(k)|. Thus: M 2 |P||x(k)|2(q−1) + 2M|Pc ||x(k)|q−1 − min (Q) ≤ 0

(63)

and by the quadratic formula: √ 1/(q−1) −M|Pc | +  ) |x(k)| ≤ ( M 2 |P|

(64)

 = M 2 |Pc |2 + M 2 |P|min (Q) Since P =

∞

k=0

k˜ k ˜ ≡ (˚c ) Q(˚ c ) solves this equation, where Q

˜ x + LT W ˜ u L, we can write W |P| ≤

∞

T k˜ k  |(˚c ) Q(˚ c) | k=0

≤



 

˜ ␭max Q 1 − ˆ 2

=: P

(65)

˜ is the maximum eigenvalue of Q. ˜ Here, we assume where ␭max Q that the maximum singular value of c , ˆ ∈ [0, 1). Similarly, we have: ∞ 

T

|˚c P| ≤ |˚c |

k

T k |(˚c ) Q˜ (˚c ) | ≤  ˆ P

(66)

k=0

and thus from (64):

 |x(k)| ≤

− ˆ P+



2

( ˆ P ) + min (Q)P P M

1/(q−1) =: cf

(67)

which will be used to define . In contrast to the approach in Section 3.3, computing M and q and satisfying (62) tends to be less expensive and scales well for systems with large nx . On the other hand, a terminal region based on the inequality (62) leads to a more stringent requirement than (56). As we will observe in Section 5, the nonlinearity bound method can lead to smaller terminal regions. 4. Nominal and robust stability of discrete QIH-NMPC In NMPC formulation, the optimal control problem given by Eqs. (3)–(11) is solved at every sampling instant k over the horizon [k, ∗ k + N]. The input u(k) is chosen as u(k) = u (k) and implemented on the plant. To establish the nominal stability of the closed loop system, we first establish that feasibility of the discrete QIH-NMPC optimization problem at k = 0 implies feasibility at all k > 0. This result is later used to establish asymptotic stability of the closed loop system at the origin. 4.1. Recursive feasibility Feasibility of the optimization problem means that, for any k ≥ 0, there exists at least one input profile U[k,k+N−1] ∈ U such that the resulting trajectory using state dynamic equations (7) and (8) satisfies the terminal inequality constraint (11) with value of the objective functional given by Eq. (4) being bounded. We state the following lemma: Lemma 13. For the nominal system with no uncertainty in the measurement, no disturbance, and perfect state feedback, feasibility of the discrete time QIH-NMPC problem (3) with (4) subject to (7)–(11) at sampling instant k = 0 implies its feasibility for all k > 0. ∗

Proof. It is assumed that the optimal solution U[0,N−1] ∈ U of optimal control problem (3) with cost function defined by (4) subject to constraints given by (7)–(11), exists and is computed at k = 0. This finite horizon optimal input profile drives the system (7) from initial state x(0) into terminal region  along the optimal trajectory

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52 ∗

z∗ (j + 1) = F(z∗ (j), u (j)) where j = 0, 1, 2, . . ., N − 1. In NMPC implementation, the manip∗ ulated input at the sampling instant k is given as u(0) = u (0). Thus, for the nominal system under absence of disturbances, state at sampling instant (k = 1) is given as x(1) = z*(1). For solving the QIH-NMPC problem at sampling instant k = 1 with initial condition z(1) = x(1), following feasible input profile is chosen: ∗

U[1,N] = {U[1,N−1] , −Lz(N)}

(68)

where L is the state feedback gain used for stabilizing the linearized model at origin. According to the result from Lemma 6 or 8, under the linear state feedback u(k) = − Lx(k), set  is invariant for the nonlinear system (2). Thus, the choice of u(N) = −Lz(N) together with the fact z*(N) ∈  implies that z(N + 1) ∈ . Since the input constraints are satisfied in , the input profile given by (68) is feasible. Accordingly, feasibility of the open loop optimal control problem (3) with (4) subject to (7)–(10) at sampling instant k = 0 implies its feasibility at sampling instant k = 1. Thus, by induction, feasibility at sampling instant k implies its feasibility at all the future sampling instants k + 1 and so on. 䊐 4.2. Asymptotic stability To establish nominal asymptotic stability of discrete time QIHNMPC formulation, the following results from Rawlings and Mayne [5] are useful. Note that the objective function of QIH-NMPC is of the form



k+N−1

J(x(k), U[k,k+N−1] ) =

37

To begin with, it is shown that the optimal cost function is nonincreasing. We later invoke Lemma 14 and Proposition 15 to show that the optimal cost function is a Lyapunov function. Theorem 16. Let (a) Assumptions A1 to A7 hold and (b) the discrete time QIH-NMPC problem is feasible at k = 0. Then, for the nominal system (2) in the absence of disturbances and controlled with QIH-NMPC is exponentially stable at the origin. Let XN ⊆ X ⊂ Rnx denote the set of initial conditions such that condition (b) is satisfied. Then, XN is a region of attraction for the closed loop system. ∗

Proof. Given the optimal input profile, U[k,k+N−1] , and optimal predicted state sequence ∗ z∗ (k + j + 1) = F(z∗ (k + j), u (k + j)) z∗ (k) = x(k) for j = 0, 1, . . ., N − 1, we have



k+N−1

∗

J ∗ (x(k), U[k,k+N−1] ) =

∗

(|z∗ (i)|2Wx + |u (i)|2Wu ) + |z∗ (k + N)|2P

i=k

(71) ∗

Since NPMC implements u(k) = u (k) on the system (2), in absence of any disturbances, it follows that x(k + 1) = z*(k + 1). Now, at instant (k + 1), consider a feasible input profile defined as follows ∗

U[k+1,k+N] = {U[k+1,k+N−1] , −Lz(k + N)}

(72)

With this choice of feasible input profile, since x(k + 1) = z*(k + 1), it follows that the predictions generated using (7) over window [k + 1, k + N + 1] satisfy z(i) = z∗ (i) for (k + 1) ≤ i ≤ (k + N)

l(x(i), u(i)) + Vf (x(k + N))

The following results relate bounds on l(x(k), u(k)) and

Since z(k + N) ∈ , which is an invariant set for the nominal system, using inequality (38) over the interval [k + N, k + N + 1], we get

Lemma 14.

|z(k + N + 1)|2P ≤ |z∗ (k + N)|2P − |z∗ (k + N)|2Q∗

i=k

[5]

Suppose Assumptions A1 to A7 hold. In addition, the following conditions hold:

min

u(k) ∈ U

k+N

= 



i=k+1

Vf (x(k + 1)) + l(x(k), u(k))|x(k + 1) ∈ 

≤ Vf (x(k)),



|z (i) |2Wx + |u (i) |2Wu + |z(k + N + 1)|2P

(69)

+|z∗ (k + N)|2W + |u(k + N)|2W x

and J(x(k + 1), U[k+1,k+N] )



≤

∗

= J ∗ (x(k), U[k,k+N−1] ) − (|x(k)|2W + |u∗ (k)|2W ) + |z(k + N)|2Q ∗ x

u

+ |z(k + N + 1)|2P − |z∗ (k + N)|2P

∀x(k) ∈  (70)

Also, since u(k) ∈ U

(75)

|z∗ (k + N)|2Wx + |u(k + N)|2Wu = |z∗ (k + N)|2Q ∗



u

Since z*(k + N) ∈  and u(k + N) is chosen as − Lz*(k + N) it follows that,

ˇ( x(k) ) for all x(k) ∈ X. Note that from (38), we have

min

u

+|z(k + N + 1)|2P − |z∗ (k + N)|2P

[5]

Vf (F(x(k), −Lx(k))) + l(x(k), −Lx(k)) ≤ Vf (x(k)),

(74)

∗

x

∀x(k) ∈ 

Suppose that Assumptions A1 to A7 hold and additional conditions of Lemma 14 hold. Also,  contains origin in its interior, and that  ⊆ X where X is a compact set in Rnx . If there exists a class K∞ function ˛(·) such that J(x(k), U[k,k+N−1] ) ≤ ˛(|x(k)|) for all x(k) ∈ , then there exists another class K∞ function ˇ(·) such that J(x(k), U[k,k+N−1] ) ≤



= J ∗ (x(k), U[k,k+N−1] ) − (|x(k)|2W + |u∗ (k)|2W )

Then, J(x(k), U[k,k+N−1] ) ≤ Vf (x(k)) for all x(k) ∈ . Proposition 15.

Now, consider the cost function at instant (k + 1) for the chosen input profile given by Eq. (68) J(x(k + 1), U[k+1,k+N] )

(i) l(·) and Vf (·) are continuous, l(0, 0) = 0 and Vf (0) = 0 (ii)  is control invariant set for the nominal system, and



(73)



Vf (F(x(k), u(k))) + l(x(k), u(k))|x(k) ∈ , x(k + 1) ∈ 

Combining inequality (73) with (76), we arrive at J(x(k + 1), U[k+1,k+N] )

it follows that inequality (69) holds for all x(k) ∈ . We now state the main theorem. A candidate Lyapunov function is constructed using the optimal cost function of the QIH-NMPC.

∗

≤ J ∗ (x(k), U[k,k+N−1] ) − (|x(k)|2W + |u∗ (k)|2W )



Vf (F(x(k), −Lx(k))) + l(x(k), −Lx(k))|x(k) ∈ 

(76)

x

(77)

u

∗

Let U[k+1,k+N] denote the optimal solution of the discrete QIHNMPC. Then, the optimality of the solution implies that ∗

J ∗ (x(k + 1), U[k+1,k+N] ) ≤ J(x(k + 1), U[k+1,k+N] )

(78)

38

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

where c and e represent constants that depend on sets  and X. Thus, we have

Combining inequalities (77) and (78), we arrive at J ∗ (x(k

∗ + 1), U[k+1,k+N] ) − J ∗



∗ x(k), U[k,k+N−1]



(79)

≤ −(|x(k)|2W + |u∗ (k)|2W ) x

u

V (x(k + 1)) − V (x(k)) < −min (Wx )|x(k)|2

(|x(k)|2Wx + |u∗ (k)|2Wu ) > 0 except at the equilibrium point x(k) = 0 and u(k) = 0. Thus,   ∗ J ∗ x(k), U[k,k+N−1] is a strictly decreasing function along a trajectory starting from any x(0) ∈ XN . Now, when x(0) = 0, it follows from Assumption A2 that the optimal solution of the discrete time QIH-NMPC is ∗



0, 0, . . ., 0





∗





2



k+N−1

¯ ∗ (i)|2Wu ) (|z∗ (i)|2Wx + |u

+ Vf (x(k))

(80)

for all x(k) ∈ XN . Thus, there exists a class K∞ function (·) such that

∀x(k) ∈ 

with (x(k)) = min (Wx )|x(k)|2 . From Lemma 14, it follows that J(x(k), U[k,k+N−1] ) ≤ Vf (x(k)) ≤ max (P)|x(k)|2 ,

∀x(k) ∈ 

Since X is a compact set, from Proposition 15 and using arguments from the proof of Proposition 2.18 from Rawlings and Mayne [5], it follows that there exists a class K∞ function ˇ(·) such that J(x(k), U[k,k+N−1] ) ≤ ˇ(|x(k)|),

Remark 17. In the proposed approach, the terminal set is an ellipsoidal region. When P tends to be a dense matrix, handling terminal constraint can turn out to be a computationally difficult task. In such situations, a polytopic approximation of the terminal region,

4.3. Robust stability

i=k

(|x(k)|) ≤ V (x(k)) ,

Inequalities (82) and (83) imply that the nominal system is, in fact, exponentially stable at the origin (ref. [5]). 䊐

ˆ Since  ˆ is a proper subset of , can be changed to z(k + N) ∈ . the requirement that  is an invariant set is still satisfied and the stability proof holds. While such polytopic approximation can be more conservative, on-line handling of the terminal constraint can get considerably simplified by this modification.

From Eq. (71), it follows that min (Wx ) x (k) ≤ V (x(k)) =

(83)

ˆ ⊂  can be ˆ that is contained inside the terminal set, i.e.  say , constructed and the terminal constraint in discrete time QIH-NMPC

for all k > 0. We define a candidate control Lyapunov function V (x(k)) ≡ J ∗ x(k), U[k,k+N−1]

(82)

and

Now, since Wx > 0 and Wu > 0,

U[k,k+N−1] =

c 2 min (Wx )|x(k)| ≤ V (x(k)) ≤ ( )max (P)|x(k)|2 e

∀x(k) ∈ XN

In the previous section, the nominal stability under discrete time QIH-NMPC has been established in the absence of any disturbance input. This modeling assumption is relaxed in this sub-section and the nominal model based controller is shown to be input to state stable (ISS) when the system is subjected to bounded disturbances. Consider a scenario where the Assumption A6 is relaxed and the system under consideration evolves as x(k + 1) = F(x(k), u(k), d(k))

(84)

where d(k) ∈ D ⊂ Rnd represents unknown disturbance input. It is to be noted that, under the scenario d(k) = 0, the system dynamics reduces to (2), i.e. F(x(k), u(k), 0) ≡ F(x(k), u(k)). The QIH-NMPC formulation, on the other hand, is based on the nominal model, i.e.

Moreover, inequality (79) implies that

z(k + j + 1) = F(z(k + j), u(k + j), 0)

(85)

V (x(k + 1)) − V (x(k)) < − (|x(k)|)

z(k) = x(k)

(86)

and for any 0 ≤ l ≤ k we have

j = 0, 1, . . ., N − 1

(87)

k 

V (x(k)) − V (x(l)) < −

A definition of input to state stability is given as follows:  (|x(k)|)

Definition 18.

j=l

Let ␾(k, x(0), d) denote the system trajectory starting from x(0) with disturbance sequence d. System (84) is input to state stable if there exist a class KL function ı and class K function  such that for all  initial conditions x(0)  and disturbance sequences d, where d[k] = d(0), d(1), . . ., d(k) , Eq. (88) is satisfied:

Thus, it follows that 0<

∞ 

 (|x(k)|) < V (x(0)) − V (x(∞)) ≤ V (x(0))

j=0

i.e.

∞ j=0

 (|x(k)|) is bounded from above. The properties of V ( · )

 

/ 0, (ii) (|x(k)|) ≤ are: (i) V 0 = 0 and V (x(k)) > 0 when x(k) = V (x(k)) ≤ ˇ(|x(k)|), ∀x(k) ∈ XN where (.) and ˇ(·) are class K∞ functions and (iii) V (x(k)) is a strictly decreasing function along a trajectory starting from any x(0) ∈ XN . Thus, V (x(k)) is a Lyapunov function for (2) controlled with QIH-NMPC in the absence of disturbances. Then, from Theorem 4, the nominal system controlled with discrete time QIH-NMPC is asymptotically stable at the origin, i.e. |x(k)| → 0 as k→ ∞ for any trajectory starting from any x(0) ∈ XN and XN is a region of attraction for the closed loop system. From the proof of Proposition 15 presented in [5], it can be seen that c c ˇ(|x(k)|) = ( )˛(|x(k)|) = ( )max (P)|x(k)|2 e e

[21,22]

(81)





| (k, x(0), d)| ≤ ı(|x(0)|, k) + (d[k−1] ) Definition 19.

(88)

[22]

A function V : Rnx → R≥0 is said to be an ISS Lyapunov function for the system x(k + 1) = G(x(k), d(k)) if there exist class K∞ functions ˛1 (·), ˛2 (·), ˛3 (·) and a class K function (·) such that ∀x(k) ∈ X for all d ∈ D, function V(x(k)) has the following properties: V (x(k)) ≥ ˛1 (|x(k)|) V (x(k)) ≤ ˛2 (|x(k)|) V (G(x(k), d(k))) ≤ V (x(k)) − ˛3 (|x(k)|) + (|d|)

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

Theorem 20.

|F(x(k), u(k), d(k)) − F(x(k), u(k), 0)| ≤ Ld |d(k)|

[22,5]

An autonomous discrete time system x(k + 1) = G(x(k), d(k)) is ISS, if there exists an ISS Lyapunov function on region of attraction XN . To establish ISS property of the discrete time QIH-NMPC, we introduce following additional assumptions: A8 The uncertainties d(k) ∈ D are sufficiently small, the set D is compact and contains the origin. A9 The function F(x(k), u(k), d(k)) is continuously differentiable w.r.t. d(k). A10 The set X ⊂ Rnx is robustly positive invariant for F( · , · , · ), i.e. F(x(k), u(k), d(k)) ∈ X holds for all x(k) ∈ X and d(k) ∈ D for some controller function ␬(.) with u(k) = ␬(x(k)).  ∗ A11 The optimal value function J ∗ x(k), U[k,k+N−1] is a uniformly continuous function w.r.t. x(k) for all x(k) ∈ XN . Note that Assumption A10 is on the lines of Remark 4 in [21], which asserts existence of robustly positive invariant set when the uncertainties are sufficiently small (Assumption A8). From uniform continuity of J ∗ ( · , · ) w.r.t. x(k) and from compactness of the set X (i.e. Assumption A3), it follows that J ∗ ( · , · ) is Lipschitz continuous for all x(k) ∈ XN . This implies that optimization problem defined by (3)–(11) satisfies Mangasarian–Fromovitz constraint qualification (MFCQ) and strong second order sufficient conditions (SSOSC) [23]. If problem QIH-NMPC satisfies strong second order sufficient conditions (SSOSC) and the Mangasarian-Fromowitz constraint qualification (MFCQ), then the value function J(· , ·) is Lipschitz continuous with respect x(k). As shown in [8], these properties hold for QIH-NMPC if the state constraints (9) are replaced by soft constraints and exact penalty functions. This is a key result, as counter-examples in [8,24] show that NMPC is not robust if SSOSC and MFCQ are not satisfied. Theorem 21. Let Assumptions A1–A5, A7–A11 hold. Also, let the nominal NMPC be asymptotically stable. If the system ((1)) admits an ISS Lyapunov function, then the system is input to state stable under the nominal NMPC with bounded disturbance. Proof. Consider ISS  a candidate  Lyapunov function ∗ V (x(k)) ≡ J ∗ x(k), U[k,k+N−1] From the proof of Theorem 16, it follows that there exists K∞ functions (.) and ˇ(.)such that

∀x(k) ∈ XN

(|x(k)|) ≤ V (x(k)) ≤ ˇ(|x(k)|),

=

 ∗ J

+

= V (x(k + 1)) − V (x(k))



∗





∗

= J ∗ x(k + 1), U[k+1,k+N] − J ∗ x(k), U[k,k+N−1] ∗



∗



F(x(k), u(k), d(k)), U[k+1,k+N]

 ∗ J

F(x(k), u(k), 0), U[k+1,k+N]



∗

− J ∗ F(x(k), u(k), 0), U[k+1,k+N]



∗

− J ∗ x(k), U[k,k+N−1]

 (89)



(90)

∗

≤ −(|x(k)|) = −min (Wx )|x(k)|2

(91)

From uniform continuity of J*(· , ·) w.r.t. x(k) (i.e. Assumption A11), it follows that there exists a class K∞ function,  Jx (·) such that ∗





∗

J ∗ F(x(k), u(k), d(k)), U[k+1,k+N] − J ∗ F(x(k), u(k), 0), U[k+1,k+N]



(93)

Combining Eq. (93) into Eq. (92), we get





∗



∗

J ∗ F(x(k), u(k), d(k)), U[k+1,k+N] − J ∗ F(x(k), u(k), 0), U[k+1,k+N]

 

≤ Jx Ld sup

  d[k] 



(94)

Since set D is compact, there exists a bound d∞ = sup|d(k)|, J

 ∗

∗ F(x(k), u(k), d(k)), U[k+1,k+N]



∗

k≥0

 

− J ∗ F(x(k), u(k), 0), U[k+1,k+N] ≤ Jx [Ld d∞ ] and this leads to the following inequality V (x(k + 1)) − V (x(k)) ≤ −min (Wx )|x(k)|2 + Jx [Ld d∞ ]

(95)

Thus, V (x(k)) becomes an ISS Lyapunov function and using Theorem 20 QIH-NMPC based on the nominal model is input to state stable. 䊐 5. Simulation studies Efficacy of the proposed discrete time QIH-NMPC formulation is demonstrated using the following three benchmark systems: • Two state system given by Chen and Allgöwer [1]. • Two state CSTR system given by Hicks and Ray [25]. • Distillation column system given by Leer [18]. The first two examples are used to illustrate terminal region characterization using Lemmas 6 and 8. The distillation column is a moderately large dimensional system and is used to illustrate approach given in Section 3.4 for terminal region characterization for large scale systems. In all three case studies, the system dynamics are modeled as a set of ordinary differential equations (ODEs) of the form dx = f(x(t), u(t)) dt

(96)

Under the assumption that the manipulated inputs are piecewise constant,

where T represents the sampling interval and k = 0, 1, 2, . . ., a computer control oriented discrete time model of the form (2) is developed where



J ∗ (F(x(k), u(k), 0), U[k+1,k+N] ) − J ∗ (x(k), U[k,k+N−1] )



  d[k] 



Since the nominal NMPC is assumed to be asymptotically stable, it follows that ∗



≤ Ld sup

u(t) = u(k) for kT ≤ t < (k + 1)T

Now, consider V (x(k))

39



and x(k) ≡ x(kT). Unless indicated otherwise, F(x(k), u(k)) is approximated using variable step size Runge-Kutta solver ode45 of MATLAB with tolerance of 10−8 . For the nonlinearity bound method, orthogonal collocation on finite elements is used instead. Also, the partial derivatives of F(x(k), u(k)) at the origin are computed as follows:

 ˚

From Assumption A9 and compactness of set D (Assumption A8), it follows that there exists a Lipschitz constant Ld such that

=

 

(92)

f(x( ), u(k))d kT



≤ Jx |F(x(k), u(k), d(k)) − F(x(k), u(k), 0)|

(k+1)T

F(x(k), u(k)) ≡ x(k) +

=

⎡



∂F = exp ⎣T ∂x 0,0



∂F = ∂u 0,0

0

T



⎤



∂f ⎦ ∂x 0,0

⎡

⎤⎡ ⎤   ∂f ⎦ ⎣ ∂f   ⎦ d exp ⎣ ∂x 0,0 ∂u 0,0 



To illustrate applicability of Lemmas 6 and 8 for discrete time formulations for different choices of sampling intervals, it becomes

40

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

necessary to use equivalent controller weighting matrices to facilitate comparison. Thus, in the first two simulation examples, the controller weighting matrices are chosen as follows: Let (Wcx , Wcu ) denote weighting matrices in a continuous time QIH-NMPC formulation. Then, the equivalent weighting matrices for a discrete time QIH-NMPC formulation are computed according to the following conversion formulae given by Astrom and Wittenmark [9]: Wx = T.Wcx

Wu =

Wcu T

where T represents the sampling interval. Also, to illustrate the arbitrary controller based approach for characterization of , it becomes necessary to design a stabilizing controller. To simplify the controller design, we have used LQR designed using (Wx , Wu ) as the weighting matrices. To characterize  using the arbitrary controller based approach, matrix Q > 0 can be selected in an arbitrary manner. However, for the sake of simplifying analysis of results, Q it is parameterized as follows Q = x Wx + u LT Wu L where x > 0 and u ≥ 0 are positive scalars and (Wx , Wu ) are QIH-NMPC weighting matrices. Similarly, in the LQR based char˜ x, W ˜ u ) are selected as acterization of , the weighting matrices (W follows: ˜ x = x Wx W

5.1. Chen and Allgöwer system [1] Consider the two state system given by Chen and Allgöwer [1]. This example was used by many researchers because of its simple two state and single input structure. Graphical visualization is possible because of two state structure. Control problem is challenging as the open loop dynamics of the system are unstable, however stabilizable. Equations governing the system dynamics are given as follows: x˙ 1 = x2 + u(0 + (1 − 0 )x1 )

(97)

x˙ 2 = x1 + u(0 − 4(1 − 0 )x2 )

(98)

where 0 = 0.5. The steady state operating point under consideration corresponds to x = [0 0]T and u = 0. The open loop dynamics of the linearized model at the origin are unstable. Jacobian linearization at the origin yields the following matrices



1.0453

0.3045

0.3045

1.0453





and  =

0.1749



0.1749

and poles of the open loop system at the equilibrium point are (0.7408, 1.3499). The input constraints for the QIH-NMPC formulation are U = {u ∈ R| − 2 ≤ u ≤ 2}. The discrete time QIH-NMPC weighting matrices in the simulation study presented in this section are chosen as follows Wx = T.(10I2 ),

Wu =

0.05 T

Sampling time T

Gain L

Closed loop poles

0.1 s 0.2 s 0.3 s 0.4 s 0.5 s

[2.2534 2.2534] [2.6202 2.6202] [2.7186 2.7186] [2.5850 2.5850] [2.3605 2.3605]

(0.9048, 0.8682) (0.8187, 0.6413) (0.7408, 0.3987) (0.6703, 0.2205) (0.6065, 0.1174)

Table 2 Chen and Allgöwer example: values of ˛ for arbitrary controller based approach. x

u

T = 0.1 s

T = 0.2 s

T = 0.3 s

T = 0.4 s

T = 0.5 s

1 5 10 15 30 50 100 50 50 50 50 50

0 0 0 0 0 0 0 1.1 10 20 50 100

0.0167 0.6597 1.7195 2.7673 5.9940 10.3763 21.1493 11.5838 20.9694 29.3960 47.2042 64.7503

0.0356 0.7670 1.7220 2.6731 5.5245 9.3359 18.8438 9.7266 12.9735 16.5224 26.0301 38.4791

0.0698 0.9071 1.9323 2.9541 6.0192 10.0582 20.1213 10.2642 11.9468 13.8187 19.2238 27.3977

0.1244 1.1330 2.3393 3.5367 7.1350 11.8939 23.9157 12.0335 13.0133 14.1266 17.3774 22.5590

0.1776 1.3572 2.7764 4.1841 8.3950 14.0519 28.0984 14.1102 14.7016 15.4676 17.3834 20.6231

Table 3 Chen and Allgöwer example: values of ˛ for LQR based approach.

˜ u = u Wu and W

where x > 1 and u ≥ 1.

˚=

Table 1 Chen and Allgöwer example: arbitrary controller based characterization of : controller gain matrices, L, and the corresponding closed loop poles.

(99)

Aim of the simulation exercise to gain insight into effect of sampling time, T, and tuning parameters (x , u ) on the size (area) of the terminal set, . Thus, in each case, sampling time is varied as T = {0.1, 0.2, 0.3, 0.4, 0.5}. 5.1.1. Arbitrary controller based characterization of  In the arbitrary controller based approach, the first step is to design a linear controller. The controller gain matrices for different

x

u

T = 0.1 s

T = 0.2 s

T = 0.3 s

T = 0.4 s

T = 0.5 s

5 10 50 200 50 50 50 50 50

1 1 1 1 1.1 10 20 50 100

0.2284 0.6305 2.6745 8.4585 3.0380 19.9910 30.0170 46.4670 60.9550

0.3892 0.9965 5.3699 21.3974 5.4682 11.5326 16.1056 25.3854 35.6083

0.6587 1.5903 8.7468 35.5610 8.7881 11.2906 13.5543 18.8100 25.2285

0.9136 2.1382 11.7713 47.7064 11.7687 12.9082 14.1244 17.0736 21.0130

1.1269 2.5833 14.1061 57.3391 14.1164 14.8088 15.3269 17.0555 19.5053

values of T and the corresponding closed loop poles are reported in Table 1. To begin with, x was changed by keeping u = 0. Table 2 presents values of ˛. Fig. 1 presents variation of area of  as a function of x (with u = 0) for different values of sampling interval, T. Since the area was found to saturate around x = 50, parameter u was varied by keeping x = 50. The results are reported in Table 2. Fig. 2 presents variation of area of  as a function of u for different values of T (when x = 50). This figure indicates that term u LT Wu L has a significant influence on the size of the terminal region. Finally, Fig. 3 underscores the premise that the terminal sets obtained in the discrete time scenarios are significantly different than the terminal set obtained using Chen and Allgöwer’s method [1] for the continuous time QIH-NMPC. We note that the area of  obtained using Chen and Allgöwer’s approach is 0.1852, which is smaller than areas of  obtained using the proposed approach (Area = 0.2534 for T = 0.1, Area = 0.2138 for T = 0.3). 5.1.2. LQR based characterization of  Unlike the arbitrary controller based characterization approach, which uses a fixed controller for a fixed T, the controller gain matrix ˜ x, W ˜ u ). Thus, the conLLQ in this case depends on the choices of (W troller gain matrix LLQ (computed using dlqr function of MATLAB) is different for each choice of T, x and u . First, x was varied by keeping u = 1. Table 3 presents values of ˛. Fig. 4 presents variation of area of  as a function of x when u = 1. Similar to the arbitrary controller based approach, the area in this case increases with increase in the sampling interval.

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

41

Fig. 1. Chen and Allgower example: variation of area of  for arbitrary controller based approach as a function of x with u = 0.

Fig. 2. Chen and Allgower example: variation of area of  for arbitrary controller based approach as a function of u with x = 50.

Fig. 3. Chen and Allgöwer example: comparison of  for T = 0.1 and T = 0.3 obtained using arbitrary controller based approach (x = 50, u = 100) with  obtained by Chen and Allgower (1998).

We observe that the area does not increase significantly after x = 50. Thus, x was kept constant at 50 and parameter u was varied. The results are summarized in Table 3. Fig. 5 presents variation of area of  as a function of u when x = 50. The qualitative trend observed in this figure is similar to the trend observed in Fig. 2.

When x is increased while keeping u = 1, the LQR controllers become increasingly aggressive. In this case, for a fixed area x , area of  increases as sampling interval, T, increases. However, when u is changed with x fixed, the corresponding LQR controllers become less and less aggressive and this leads to significant increase in size

42

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

Fig. 4. Chen and Allgower example: variation of area of  for LQR based approach as a function of x and T with u = 1.

Fig. 5. Chen and Allgower example: variation of area of  for LQR based approach as a function of u and T with x = 50.

of . Also, for a fixed u , smaller sampling time gives larger area of . Fig. 6 compares  obtained using LQR based approach (T = 0.3 s, x = 50, u = 50) with



LLQ = 2.4279

2.4279





and

PLQ =

397.3556 64.9002

64.9002



397.3556 (100)

and  obtained using arbitrary controller based approach (T = 0.3 s, x = 50, u = 100) with



L = 2.7186 2.7186





and

P=

408.4463

69.3418

69.3418

408.4463



(101)

The terminal region obtained using the method of Chen and Allgöwer [1] is plotted in Figs. 3 and 6 as a reference. Note that the largest terminal sets obtained using two different approaches are almost identical. Now, we consider terminal region calculations for system given by Eqs. (97) and (98) using a nonlinearity bound according to Eq. (17), found via simulations and terminal region calculated via Eq. (67). We hold u = 1 constant, and so the nonlinearity bound must be found for each combination of sampling time T and x . We show this in detail for one case, T = 0.1 s and x = 1.1, in Table 7 . This is done by simulating the system one step forward under LQR control from a random sampling of 10,000 initial conditions, and then choosing M and q to bound the data. After this is completed for each case, the terminal region area may be found as cf2 , with results

shown in Fig. 8. The nonlinearity bound parameters M and q as well as cf are summarized in Table 4 for each case. In this example, it was observed that, in the case of nonlinearity bound method, the terminal region area increases with increase in the sampling interval. For smaller sampling intervals (i.e. T = 0.1, 0.2), areas of the terminal regions obtained using nonlinearity bound method is significantly smaller than the areas of the terminal regions obtained using LQR based and arbitrary controller based approaches. However, the terminal regions obtained using nonlinearity bound method become comparable to areas obtained using other two approaches for larger sampling interval (i.e. T = 0.4, 0.5). Table 5 summarizes results (largest area of ) obtained for different combinations of characterization approaches, tuning parameters, and sampling intervals. From this table, it is evident that, for a fixed sampling time, providing additional degree of freedom helps in increasing the size of . Also, when only a single degree of freedom is available, increasing sampling interval leads to increase in the size of . The largest terminal regions obtained using LQR based approach and Arbitrary Controller based approach are comparable to each other. Remark 22. Yu et al. [13] have demonstrated efficacy of their approach using a system obtained by discretizing model (97) and (98) using explicit Euler method with sampling interval equal to 0.1. Table 6 summarizes results obtained when the proposed approaches are used for terminal region characterization of the discrete time system in [1]. It is to be noted that the controller

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

43

Fig. 6. Chen and Allgöwer example: comparison of  for T = 0.3 obtained using arbitrary controller based approach (x = 50, u = 100) and LQR based approach (x = 50, u = 100).

Fig. 7. Chen and Allgöwer example: nonlinear bound for T = 0.1 s, x = 1.1.

Fig. 8. Chen and Allgöwer example: terminal region calculations with nonlinearity bound.

44

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

Table 4 Chen and Allgöwer example: terminal region results with nonlinearity bound (ln(M), q, cf ). T

x = 1.1

x = 5

x = 10

x = 50

x = 200

0.1 0.2 0.3 0.4 0.5

−0.5,2,1.18e−3 0.4,2.05,5.15e−3 0.9,2.1,1.15e−2 1.2,2.15,0.020 1.4,2.15,2.82e−2

−0.3,2,2.03e−2 0.7,2.05,6.77e−2 1,2.1,0.129 1.3,2.15,0.186 1.45,2.15,0.233

−0.1,2,2.78e−2 0.8,2.05,8.73e−2 1.1,2.1,0.156 1.35,2.15,0.214 1.5,2.15,0.262

0.3,2,4.22e−2 0.9,2.05,0.118 1.15,2.1,0.184 1.4,2.15,0.238 1.55,2.15,0.286

0.6,2,5.39e−2 1,2.05,0.129 1.2,2.1,0.190 1.45,2.15,0.243 1.6,2.15,0.290

Table 5 Chen and Allgöwer example: comparison of maximum terminal regions. Approach

Degrees of freedom

Area of  (T = 0.1 s)

Area of  (T = 0.3 s)

Area of  (T = 0.5 s)

Arbitrary controller Arbitrary controller LQR (inequality method) LQR (inequality method) LQR (nonlinearity bound)

x x and u x x and u x

0.1288 0.2534 0.0481 0.2502 0.009127

0.1262 0.2138 0.1211 0.2022 0.1134

0.1364 0.1671 0.1406 0.1614 0.2642

Table 6 Chen and Allgöwer example: comparison of terminal regions using LQR based approach, arbitrary controller based approach and approach by Yu et al. [13]. Approach

Parameter



˛

Area

Yu et al. (2007) LQR LQR Arbitrary controller Arbitrary controller

= 1/0.91 x = 100, u = 1 x = 100, u = 20 x = 100, u = 0 x = 100, u = 20

123.66 13.2126 380.2826 145.2788 541.4081

0.778 4.6069 40.8432 10.4295 38.8676

0.0271 0.0711 0.2411 0.1148 0.2211

Table 7 Initial conditions for the states for Chen and Allgöwer example.

Fig. 10 shows plots of manipulated input trajectories for different initial conditions. Notice that from initial conditions P1 − P6, the state trajectories reach inside the terminal region within 6 samples as these points belong to set X6 . When the terminal region is computed using the approach given by Chen and Allgöwer [1], the sampling interval has to be chosen small (T = 0.1 in [1]) to keep the digital implementation of QIH-NMPC close to the continuous time formulation. This requirement, however, results in a larger prediction horizon (N = 15 in [1]), which significantly increases the on-line computation time. 5.2. CSTR system [25]

Point

x1 value

x2 value

P1 P2 P3 P4 P5 P6

−0.683 −0.523 0.808 0.774 0.292 −0.080

−0.864 0.744 0.121 −0.222 0.228 −0.804

We now consider continuous stirred tank reactor (CSTR) presented by Hicks and Ray [25]. Dynamics of this system are governed by the following set of ODEs [17]:



dz c (1 − zc ) Ea = − k0 zc exp − m2 zT dt





P=

686.5122 407.9211

651.3010 375.4682

407.9211 686.5122

375.4682 651.3010





and PLQ =

and parameters of the terminal set are given in Table (6). From Table 6, it can be observed that the areas of largest terminal regions computed using Lemmas 6 and 8 are 8.89 and 8.16 times larger, respectively, when compared with the terminal region area obtained by Yu et al. [13]. 5.1.3. Discrete time QIH-NMPC simulation To demonstrate the effectiveness of the proposed discrete time QIH-NMPC, time simulations of the discrete time QIH-NMPC with terminal region obtained using arbitrary controller based approach with T = 0.3,

x = 50,

u = 100

are carried out starting from six initial conditions specified in [1] (cf. Table 7). The prediction and control horizons are set to N = 6. Fig. 9 shows a phase plane diagram of system state evolution from various initial conditions along with the terminal region.

(102)



f

used is identical to the controller given by Yu et al. [13] while controller based approach and LLQ = employing the arbitrary  2.3238 2.3238 . The resulting P matrices are as follows



(z − zT ) dz T Ea + k0 zc exp − = T m2 zT dt



− ˛0 m1 (zT − zTCW )

(103)

where zc and zT represent dimensionless concentration and dimensionless temperature, respectively. Manipulated inputs are cooling water flow rate m1 and inverse of the dilution rate m2 . Nominal values of the parameters are given below: Variable

Nominal value

zTCW f zT Ea ˛0 k0

0.38 0.395 5 1.95 × 10−4 300

To improve numerical stability, the inputs (m1 , m2 ) appearing in the system dynamics are scaled as u1 = m1 /600 and u2 = m2 /40. It is assumed that the system is operated at the following steady state operating point: Xs = [0.6416 0.5387]T

and

Us = [0.5833 0.5000]T

(104)

The input constraints are given as follows U = {u1 , u2 ∈ R| − 0.4167 ≤ u1 ≤ 0.4167; −0.4750 ≤ u2 ≤ 0.5} Jacobian linearization of the discrete time nonlinear system at (Xs , Us ) yields the following pair of system matrices



˚=

0.9207

−0.3272

0.0296

1.2051





and

 =

0.0029

−0.0367

−0.0202

0.0153



C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

45

Fig. 9. Chen and Allgöwer example: phase portrait of states under discrete time QIH-NMPC.

Fig. 10. Chen and Allgöwer example: manipulated input trajectories under discrete time QIH-NMPC.

Eigenvalues of the open loop discrete time dynamics are (0.9602, 1.1656). For the discrete time NMPC formulation, weighting matrices are chosen as,



 Wx = T.

10

0

0

2

and Wu =

1 0

0



0.5

(105)

T

Sampling interval of T = 1 unit is used. Linear gain matrix and terminal penalty matrix obtained using LQR based approach (T = 1, x = 50, u = 200) are given as follows:



LLQ =

−1.3839 −9.8163 0.0401 5

PLQ = 10 ×



10.1250 0.0642



and 0.1559

(106)

0.1559 1.1502

Linear gain matrix and terminal penalty matrix obtained using arbitrary controller based approach (T = 1, x = 50, u = 35) are given as follows:



L=

−1.6296 −9.9872 −1.8185

9.6754





and P = 104 ×

0.4791

0.5285



0.5285 2.0824 (107)

Fig. 11 the largest terminal regions obtained with the approaches in Sections 3.1 and 3.2. The maximum terminal region area obtained

Table 8 Terminal region results with LQR and nonlinearity bound, x = 50. u

M

q

cf

0.1 10 20 35 50 75 100

1.097e3 54.5982 28.50 25.79 23.34 21.11 20.09

2.7 2.3 2.25 2.2 2.18 2.16 2.15

0.00231 0.00472 0.00390 0.00204 0.00137 0.00083 0.00056

using the LQR based approach is approximately 2.8 times larger than the maximum terminal region area obtained using the arbitrary controller based approach. For the nonlinearity bound method, we consider that x = 50 is held constant while u varies along with the control gain L. Fig. 12 shows the nonlinear bound (M |x|q ) for the case that u = 0.1. Values of the bound parameters (M and q) as well as cf are summarized in Table 8. Terminal region area is shown as a function of u in Fig. 13. The nonlinearity bound approach is simulation based and the parameters M and q are (possibly loose) upper bounds on nonlinearity matrices. While this approach is less dependent on the number of states, this approach is more conservative than the first approach that solves Eq. (58). Also, the terminal region determined in (67) is strongly determined by the values of M and q. Related to ˜u the trends in Fig. 13, we observed that increasing u increases W

46

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

Fig. 11. Hicks and Ray example: comparison of largest  for T = 1 obtained using arbitrary controller based approach and LQR based approach.

Fig. 12. Hicks and Ray example: LQR, nonlinearity bound with x = 50 and u = 0.1. Table 9 Terminal region results with LQR and nonlinearity bound,  = x = u .  cf

.1 0.00027

10 0.00156

20 0.00162

35 0.00164

50 0.00165

75 0.00166

Table 10 CSTR system by Hicks and Ray: comparison of maximum terminal regions. 100 0.00167

and the terminal penalty matrix P, which has a nonlinear effect on the terminal region determined by (67) through P . On the other hand, Table 9 also shows that M in this example decreased with increasing u , leading to a competing effect, which can be observed through the product P M in (67). Consequently, cf in (67) has a nonlinear dependence on both P and M, which we believe is the source of the non-monotonic behavior in Fig. 13. We note that this particular trend depends strongly on the nonlinear plant dynamics of this example. Now, we consider the case that  = x = u . Results for the parameter cf are summarized in Table 9, and results for the terminal region area are shown in Fig. 14. In this example, however, it was observed that, areas of the terminal regions obtained using nonlinearity bound method are significantly smaller when com-

Approach Arbitrary controller Arbitrary controller LQR (inequality method) LQR (inequality method) LQR (nonlinearity bound)

Degrees of freedom x = 50 x = 50 and u = 35 x = 50 x = 50 and u = 200 x = 50 and u = 10

Area of 

˛ −3

3.781 × 10 8.541 × 10−3 1.712 × 10−3 2.444 × 10−2 6.999 × 10−5

1.6779 23.0408 0.4535 173.0460 cf = 0.00472

pared with areas of the terminal regions obtained using LQR based and arbitrary controller based approaches. Table 10 summarizes area of the maximum terminal regions obtained using both arbitrary controller based approach and LQR based approach for sampling interval of T = 1 unit for the CSTR system by Hicks and Ray [25]. As observed in the case of Chen and Allgöwer system, utilizing an additional degree of freedom helps in increasing the size of terminal region. To demonstrate the effectiveness of the proposed discrete time QIH-NMPC, time simulations of the discrete time QIH-NMPC with

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

47

Fig. 13. Hicks and Ray example: LQR, terminal region area with x = 50.

Fig. 14. Hicks and Ray example: LQR, terminal region area with  = x = u .

terminal region obtained using LQR based approach with x = 50, u = 35 are performed starting from five initial conditions as given in Table 11. The prediction and control horizons are set to N = 10. Noise is added to the both the manipulated inputs by modifying the system dynamics as follows x(k + 1) = F(x(k), u(k) + d(k)) where each di (k) for i = 1, 2 is generated from uniform distribution with ±5% of the steady states as the bounds. Fig. 15 shows phase portrait from various initial conditions along with the ter-

Table 11 Initial conditions for the states for Hicks and Ray example. Point

x1 value

x2 value

P1 P2 P3 P4 P5

0.100 −0.350 −0.100 0.100 0.180

0.052 0.000 −0.051 −0.071 −0.002

48

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

Fig. 15. Hicks and Ray example: phase portrait of states under discrete time QIH-NMPC case with noise.

Fig. 16. Hicks and Ray example: input trajectories under discrete time QIH-NMPC for with noise case.

minal region for the QIH-NMPC case with noise. Notice that from initial conditions P1 − P5 , the state trajectories reach inside the terminal region within 10 samples as these points belong to set X10 . Fig. 15 also presents a zoomed version around the origin, which shows that the trajectories continue to remain in a small neighborhood of the origin and the closed loop system is input to state stable. Fig. 16 plots the corresponding manipulated input trajectories.

5.3. System of distillation columns As a high dimensional example, we consider two distillation columns in series [18]. This process separates a mixture of components A, B, and C. The bottoms of the first column are fed to the second column. The distillate of the first column is to be 95% A, the distillate of the second column is to be 95% B, and the bottoms of the second column is to be 95% C. These are implemented as inequality constraints in the optimization problem. For vapor–liquid equilibrium, constant relative volatility is assumed with constants ˛A = 2 and ˛B = 1.5. For tray hydraulics, we use the Francis weir formula, with constant Kuf = 21.65 above the feed and constant Kbf = 29.65 below the feed. The weir height is 0.25, and the nominal liquid holdup is 0.5. The feed is saturated liquid, so that qF = 1. Each column has 41 equilibrium stages including the reboiler, giving a total of 246 states and 8 controls. The distillation column system (ref. Fig. 17) separates a mixture of components A, B, and C. The bot-

Table 12 Distillation example: terminal region results.  cf

10 0.3773

102 0.4241

103 0.4288

104 0.4293

105 0.4293

toms of the first column are fed to the second column. The system dynamic equations and the associated nominal model parameters can be found in [18,8]. For the stage costs we set Wx = 10I246 and Wu = I8 , and for the discretization step we set T = 1 min. For this method we will utilize the extension for large-scale systems with nonlinearity bound (17) and terminal regions characterized by (67). Furthermore, we consider tuning the terminal regions via  = x = u , so that L does not change, and therefore the nonlinearity bound does not need to be recomputed. The results of 10, 000 one-step simulations (performed in MATLAB with ode45) are shown in Fig. 18, and the nonlinearity bound parameters are found to be q = 1.8 and M = 0.0743. Then, cf is shown as a function of  in Fig. 19, and numerical values are shown in Table 12. As can be seen, there is little benefit in terminal region size to increasing  beyond 100. Finally, we show closed-loop simulations of the distillation system under QIH-NMPC for varying values of  with N = 10. To accommodate the optimal control problem, we discretize the DAE system using three point Radau collocation. The models are implemented in AMPL and solved with IPOPT on an Intel i7-4770 3.4 GHz CPU. The simulation results can be seen in Fig. (20). All cases are con-

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

49

Fig. 17. Distillation columns: schematic diagram.

Fig. 18. Distillation column: nonlinearity bound.

vergent to the setpoint, however setting very high values for  leads to degradation in the dynamic performance. This, combined with the observation that high values of  do not significantly increase the terminal region size, leads to the conclusion that moderate values of  are the most practical. 6. Conclusions Terminal costs and terminal constraints are essential in NMPC formulations in order to establish nominal stability. A large terminal set gives rise to a large region of attraction for the closed loop system and can also help in reducing the on-line computation cost. Despite these attractive features, the methods for characterization of the terminal region have not received much attention in the literature. In this work, two novel approaches have been devel-

oped for the characterization of the terminal region for a discrete time quasi-infinite horizon NMPC formulation. Similar to the continuous time QIH-NMPC formulation by Chen and Allgöwer [1], a stabilizing linear controller developed using linearization at the origin is employed for characterization of the terminal set. The first approach permits use of an arbitrary linear controller while the second approach is based on LQR controller. Unlike the approaches available in the literature, both the approaches provide sufficient degrees of freedom to shape the terminal set. Using a method of bounding only the higher order nonlinear effects of the system via simulations under linear controller, the proposed approach is further extended to handle terminal region computation for a large dimensional system. The quadratic control Lyapunov function in the terminal set together with the quadratic stage cost is further used to establish exponential stability of the discrete time QIH-

50

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

Fig. 19. Distillation example: terminal regions with  = x = u .

Fig. 20. Distillation example: dynamic trajectories.

NMPC formulation under the nominal conditions. Additionally, it is shown that discrete time QIH-NMPC based on the nominal model is input to state stable (ISS) in the face of bounded disturbances in the state dynamics. Efficacy of the proposed approaches for characterization of the terminal region is demonstrated using a two state example from [1] and a benchmark CSTR system (or Hicks’s CSTR system, [25]). In the case of the first example, comparison of the terminal sets obtained using the method given by Chen and Allgöwer [1] and the proposed two approaches reveals that shapes and sizes of the terminal sets obtained using discrete time controllers and discrete time model

are significantly different than the terminal set obtained using the continuous time model. Parametric studies reveal the choice of sampling interval and the choices of the tuning matrices in both the approaches have significant influence on the sizes of terminal sets. In particular, when the linear controllers become more aggressive, it was observed that area of the terminal region increases with increase in the sampling interval. While the arbitrary controller based method gives much larger terminal set in the case of a two state example by Chen and Allgöwer [1], the LQR based approach performs better in case of CSTR system by Hicks and Ray [25]. Efficacy of the extension for a large dimensional system is

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

demonstrated using terminal region characterization and closed loop simulation studies on a two distillation columns in series [18]. The simulation results demonstrate that the state trajectories converge asymptotically to the origin under the nominal conditions for all three systems. Closed loop simulations using the Hicks’s CSTR system [25] and distillation column systems in the presence of bounded disturbances show that the state trajectories hover in a bounded region in the neighborhood of the origin, which is in accordance with the ISS property of the discrete time QIH-NMPC formulation. Finally, economic NMPC (eNMPC) has generated significant interest in the control community. However, developing a terminal region approach is a significant extension beyond tracking NMPC, which is the scope of the current study. eNMPC usually has stage costs that are not class K functions, so the nominal and robust stability properties developed in this study do not apply. On the other hand, terminal regions can be computed using the proposed approaches to accommodate shifting system steady states due to some known parameter change. Some preliminary work in this direction has been reported in [26]. Also, a promising approach is the work of Limon et al. have considered NMPC with a penalty term for tracking objective to the eNMPC objective function and have shown stability results (ref. [27]). Nevertheless, for economic stage costs a major reformulation of our method would be needed, and we leave this as a topic for future work. Appendix A Terminal region  can be characterized as the largest value of ˛ such that



≡



T

x(k) ∈ Rnx : x(k) Px(k) ≤ ˛, −Lx(k) ∈ U, ϕ (x(k)) ≥ 0 (108)

where ˛ ≤ . A random sampling based approach for characterization of  is proposed here. The first step is to develop a box in Rnx using



 ≡



T

x(k) ∈ Rnx : x(k) Px(k) ≤ , −Lx(k) ∈ U

(109)

If the input constraints are of the form uL ≤ −Lx(k) ≤ uH then  can be characterized using the following quadratic programming problem x∗ (k) = arg minx(k) Px(k) T

x(k)

subject to constraints



−L



x(k)≤

L

uL



−uH



T

The value of  = x∗ (k) follows: Let





xi,max = ı

Px∗ (k) is used to create a box in Rnx as



x = −ı Pii i,min

Pii

 for i = 1, 2, . . ., nx

where ı > 1, denote the corner points of the box. Then, using the uniform distribution, a sufficiently large number of samples are generated inside the box such that B=



(j)

x(j) : j = 1, 2, . . ., M; xi,min ≤ xi ≤ xi,max

We then compute the set S=



x(j)

T



Px(j) , ϕ x(j)



: j = 1, 2, . . ., M

51



The smallest value of x(j)

T





Px(j) for which ϕ x(j) ≤ 0, yields a

(l)

point, say x˜ , which is very close to the boundary of the terminal set. This point can be used as an initial guess of the optimization problems discussed in Section 3.3. Alternatively, another box of (l)

smaller volume can be created in the neighborhood of x˜ and the sampling based procedure can be repeated to arrive at a reasonably accurate estimate of the terminal region. Application of the random sampling approach is used to compute the first estimate of ˛, say ˛0 . When ˛0 ≥ , the value of ˛ is set equal to . In the case when ˛0 < , the sequential optimization procedure discussed in Step 2 (Section 3.3) is invoked. A series of optimization problems is solved starting from ˛ = ˛1 where ˛0 < ˛1 < and ˛ is gradually reduced until ϕ (x(k)) ≥ 0. Depending on the initial guess, a degenerate scenario can occur when the conventional NLP solver approaches the origin and produces a trivial solution. When such problem is encountered, the random sampling method is re-invoked to obtain another initial guess for the conventional NLP solver. These steps are repeated until the NLP solver produces a solution close to ˛0 . References [1] H. Chen, F. Allgöwer, A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability, Automatica 34 (10) (1998) 1205–1217. [2] L. Biegler, X. Yang, G. Fischer, Advances in sensitivity-based nonlinear model predictive control and dynamic real-time optimization, J. Process Control 30 (2015) 104–116. [3] D.Q. Mayne, Model predictive control: recent developments and future promise, Automatica 50 (2014) 2967–2986. [4] D.Q. Mayne, J.B. Rawlings, C.V. Rao, P.O.M. Scokaert, Constrained model predictive control: stability and optimality, Automatica 36 (6) (2000) 789–814. [5] J. Rawlings, D. Mayne, Model Predictive Control: Theory, Computation, and Design, Nob Hill Publishing, 2009. [6] H. Michalska, D. Mayne, Robust receding horizon control of constrained nonlinear systems, IEEE Trans. Autom. Control 38 (11) (1993) 1623–1633. [7] H. Chen, F. Allgöwer, A computationally attractive nonlinear predictive control scheme with guaranteed stability for stable systems, J. Process Control 8 (5–6) (1998) 475–485. [8] D.W. Griffith, L.T. Biegler, S.C. Patwardhan, Robustly stable adaptive horizon nonlinear model predictive control, J. Process Control 70 (2018) 109–122. [9] K.J. Astrom, B. Wittenmark, Computer-Controlled Systems: Theory and Design, Prentice Hall, 1997. [10] J.B. Rawlings, D.Q. Mayne, M. Diehl, Model Predictive Control: Theory, Computation, and Design, 2nd ed., Nob Hill Publishing, LLC, 2017. [11] D. Limon, Control predictivo de sistemas no lineales con restricciones: estabilidad y robustez (Ph.D. thesis), University of Seville, 2002. [12] T.A. Johansen, Approximate explicit receding horizon control of constrained nonlinear systems, Automatica 40 (2) (2004) 293–300. [13] S. Yu, T. Qu, F. Xu, H. Chen, Y. Hu, Stability of finite horizon model predictive control with incremental input constraints, Automatica 79 (2017) 265–272. [14] S. Lucia, P. Rumschinski, A.J. Krener, R. Findeisen, Improved design of nonlinear model predictive controllers, IFAC-PapersOnLine 48 (23) (2015) 254–259. [15] M. Lazar, M. Tetteroo, Computation of terminal costs and sets for discrete time nonlinear MPC, IFAC-PapersOnLine 51 (20) (2018) 141–146. [16] C. Rajhans, S.C. Patwardhan, H. Pillai, Discrete time formulation of quasi infinite horizon nonlinear model predictive control scheme with guaranteed stability, IFAC-PapersOnLine 50 (1) (2017) 7181–7186. [17] R. Huang, S.C. Patwardhan, L.T. Biegler, Robust stability of nonlinear model predictive control based on extended Kalman filter, J. Process Control 22 (1) (2012) 82–89. [18] R. Leer, Self-Optimizing Control Structures for Active Constraint Regions of a Sequence of Distillation Columns (Master’s thesis), Norwegian University of Science and Technology, 2012. [19] D.M.Y. Sommerville, An Introduction to the Geometry of N Dimensions, E. P. Dutton and Company Incorporated, 1929. [20] X.-B. Hu, W.-H. Chen, Model predictive control: terminal region and terminal weighting matrix, Proc. Inst. Mech. Engrs. Part I: J. Syst. Control Eng. 222 (2) (2008) 69–79. ˜ J.M. Bravo, A. Ferramosca, [21] D. Limon, T. Alamo, D.M. Raimondo, D.M. de la Pena, E.F. Camacho, Input-to-state stability: a unifying framework for robust model predictive control, in: Nonlinear Model Predictive Control Towards New Challenging Applications, Springer Berlin Heidelberg, 2009, pp. 1–26. [22] Z.-P. Jiang, Y. Wang, Input-to-state stability for discrete-time nonlinear systems, Automatica 37 (6) (2001) 857–869.

52

C. Rajhans et al. / Journal of Process Control 83 (2019) 30–52

[23] L.T. Biegler, D.M. Thierry, Large-scale optimization formulations and strategies for nonlinear model predictive control, IFAC-PapersOnLine 51 (20) (2018) 1–15. [24] G. Grimm, M.J. Messina, S.E. Tuna, A.R. Teel, Examples when nonlinear model predictive control is nonrobust, Automatica 40 (2004) 1729–1738. [25] G.A. Hicks, W.H. Ray, Approximation methods for optimal control synthesis, Can. J. Chem. Eng. 49 (4) (1971) 522–528.

[26] D.W. Griffith, Advances in Nonlinear Model Predictive Control for Large-Scale Chemical Process Systems (Ph.D. thesis), Carnegie Mellon University, 2018. [27] D. Limón, A. Ferramosca, I. Alvarado, T. Alamo, Nonlinear MPC for tracking piece-wise constant reference signals, IEEE Trans. Autom. Control 63 (11) (2018) 3735–3750.