Multivariable constrained process control via Lyapunov R-functions

Multivariable constrained process control via Lyapunov R-functions Aldo Balestrino a , Andrea Caiti a , Sergio Grammatico a,? , a Department

of Energy and Systems Engineering, University of Pisa, Largo Lucio Lazzarino 1, 56122 Pisa, Italy.

Abstract This paper proposes a smooth control Lyapunov function with maximal controlled domain of attraction for the state-feedback control of constrained uncertain linear systems describing the dynamics of multivariable processes. Constructive algorithms are shown to build such a control Lyapunov function within a merging procedure due to the so called R-functions. Robustness with respect to polytopic model uncertainties and disturbance rejection are also discussed. The controlled dynamics of two chemical reactions are simulated to show the benefits of the proposed strategy. Key words: Generated Lyapunov functions, Uncertain linear systems, Process control.

1

Introduction

Chemical processes are continuously faced with the requirements of becoming safer, more reliable, and more economical in operation. The design of effective chemical process control systems, inherently Multi-Input Multi-Output (MIMO) and nonlinear, needs to be both rigorous and practical [12]. Moreover, the unavoidable presence of physical constraints on the process variables and in the capacity of control actuators not only limit the nominal performance of the controlled system, but can also affect the stability of the overall system. As a consequence, the stabilization of such processes is one of the most attractive research areas for the chemical and control engineering community [11]. ? Corresponding author Sergio Grammatico, e-mail: [email protected].

Preprint submitted to Elsevier

4 June 2012

Model Predictive Control (MPC) [22], [21], also known as Receding Horizon Control (RHC) [31], [20], can handle both state and control input constraints within an optimal control setting [34], [24]. These approaches, as well as the Explicit MPC [7], [19], can be quite computationally-demanding. Hence a large literature has been developed for fast computation of sub-optimal (robust) MPC solutions, see [28], [29], [37] among others in recent literature.

In [24], an interesting Lyapunov-based MPC approach has been proposed for the control of an exothermic chemical reaction, taking place in a Continuous Stirred Tank Reactor (CSTR). In particular, a quadratic Control Lyapunov Function (CLF) is used together with an horizon-1 MPC. However, an ellipsoidal set can not accurately fit the polyhedral state constraints describing the limits on the admissible concentration of the chemical reactant and on the reactor temperature. Therefore a Quadratic CLF (QCLF) can not guarantee stability over the whole controllable invariant set [2]. On the other side, the estimate of the controlled invariant state-space region can be enlarged via the synthesis of Polyhedral CLFs (PCLFs) [8], composite-quadratic CLFs [17], [18], smoothed PCLFs [10], Truncated Ellipsoid (TE) Control Lyapunov Functions [26], [36], and smoothed TE CLFs [1], [5].

This paper considers the control of constrained CSTRs and develops a technique based on the so called Control Lyapunov R-Functions (CLRFs) for its solution. This approach has been proposed by the authors in [2] and [4] and it is here made more general. The proposed CLF has inner level curves that can be made, independently from the external one, as-close-as-desired to any choice of smooth ones. The main contribution of this paper is: a) the tuning of a free design parameter via constructive algorithms, in order to trade-off between optimality requirements and the size of the guaranteed domain of attraction; additional contributions are: b) the discussion about robustness to model uncertainties and disturbance rejection; c) the possibility to handle asymmetric domain of attractions. Unlike switching control strategies [14], within this novel approach both constraints and optimality arguments can be handled by a unique smooth CLF, at least for uncertain linear systems, together with a continuous control [2], [3], [4].

The paper is organized as follows. The state-feedback control of a constrained CSTR is presented as motivating example in the next section. Sections 3, 4 present the main theoretical tools here used, while two constructive procedures for the synthesis of suitable CLRFs are presented in Section 5. Section 6 applies the results to the CSTR case studies. In last section we conclude the paper and outline some future work. The proofs are given in Appendix. 2

1.1

Notation

In denotes the n × n identity matrix. The closed k-level set of a continuous function V : X ⊆ Rn → R, i.e. {x ∈ X : V (x) ≤ k}, is denoted by N [V, k]. A set S ⊆ Rn is called C-set if it is a convex and compact set including the origin in its interior [9]. Ir denotes {n ∈ Z+ : n ≤ r}.

2

A constrained uncertain Continuous Stirred Tank Reactor as motivating example k

Consider an irreversible, exothermic first-order reaction of the form A → B, taking place in a CSTR. The inlet stream consists of pure A at flow rate F , concentration CA0 and temperature TA0 . The dynamic model of the process is of the form E F CA C˙A = (CA0 − CA ) − k0 exp − V RTR   F ∆H E Q ˙ TR = (TA0 − TR ) − , k0 exp − CA + V ρcp RTR ρcp V 



(1)

where CA denotes the concentration of the species A, TR denotes the temperature of the reactor, Q is the heat input to the reactor, V is the volume of the reactor, k0 , E, ∆H are, respectively, the pre-exponential constant, the activation energy and the enthalpy of the reaction, cp and ρ are, respectively, the heat capacity and the fluid density in the reactor. The strategy of controlling the temperature to an equilibrium point may be unusual from a practical point of view. However, this model structure has been widely addressed in the literature [11–14,23,24] as a benchmark case study. The numerical values of the process parameters, taken from [24], are shown in Table 1. The open-loop unstable equilibrium point x¯ = [CAs , TAs ]> = [0.57 kmol/m3 , 395.3 ◦ K]> , u¯ = 0, of the dynamical model (1) has to be stabilized according to the state constraints |CA − CAs | ≤ 0.16 kmol/m3 , |TR − TRs | ≤ 3 ◦ K.

(2)

The control variables υ are the variation of the inlet concentration of species A, υ1 = ∆CA0 = CA0 − CA0s , and the heat input to the reactor υ2 = Q. These manipulated control inputs are constrained as follows: |∆CA0 | ≤ 1 kmol/m3 , |Q| ≤ 1 kJ/h. 3

(3)

Table 1 Process parameters and steady state values V

0.1

R

8.314

CA0s

1

TA0s

310

Qs ∆H

0 −4.78 · 104

m3 kJ/(kmol ◦ K) kmol/m3 ◦K

kJ/min kJ/kmol

k0

72 · 109

E

8.314 · 104

cp

0.239

kJ/(kg ◦ K)

ρ

1000

kg/m3

F

0.1

TRs

395.3

CA s

0.57

1/min kJ/kmol

m3 /min ◦K

kmol/m3

It is indeed desired to find a continuous control law υ(t), satisfying the control constraints, that drives the state ξ(t) = [CA (t), TR (t)]> to x¯, also in accordance to the state constraints. In fact, discontinuous and/or chattering control laws, such as the ones usually obtained by switching controllers, are actually not well implementable on real actuators [2]. As we are focusing on constrained uncertain linear systems, let us consider system (1) linearized in a neighborhood of the equilibrium point (¯ x, u¯), namely 







−1.7428 −0.0271  0  1  x˙ =   x +   u, 148.5626 4.4191 0 0.0418

(4)

where x = ξ − x¯ and u = υ − u¯.

2.1

Static state-feedback based on the Riccati-optimal quadratic control Lyapunov function

An optimization problem is usually considered in many control approaches for multivariable chemical processes [13], [23], [24], and, typically, a (piecewise) QCLF is designed. This particular choice is motivated by the fact that, given an unconstrained linear system x˙ = Ax + Bu, the gradient-based control u(x) = −R−1 B > P ∗ x, being P ∗  0 the (unique) solution of the Algebraic 4

398.3

397.3

TR (°K)

396.3

395.3

394.3

393.3

392.3 0.42

0.47

0.52

0.57 CA (kmol/m3)

0.62

0.67

0.72

Fig. 1. Level sets of the Riccati-optimal quadratic control Lyapunov function. The box in solid line indicates the state constraints.

Riccati Equation (ARE) A> P +P A+Q−P BR−1 B > P = 0, asymptotically stabilizes the unconstrained linear system itself, while minimizing the quadratic performance cost J(x, u) =

Z +∞ 



x(τ )> Qx(τ ) + u(τ )> Ru(τ ) dτ.

(5)

0

In the case of constrained systems, both linear and nonlinear ones, the particular choice of the CLF is a critical point in the control design, since the largest (indeed non conservative) estimate for the controllable state space set should be provided. Considering the weight matrices Q = R = I2 , the candidate QCLF deriving from the solution of the corresponding ARE for system (4) is x> P ∗ x, with   P∗ =  

20.3640 1.6312  1.6312

0.1979

.

(6)

As shown in Figure 1, since a quadratic function can not fit well the polyhedral state constraints, only a shrunk QCLF can be used for any static control design, in order to guarantee the fulfillment of the state constraints. As a result, the estimate of the largest invariant set is reduced (at most) to the largest ellipsoidal set included in the maximal controllable set.

2.2

Enlarging the controlled invariant state-space region

In the case of constrained (uncertain) linear systems, the maximal controlled set can be computed via PCLFs, for instance by using the procedure proposed in [25] for discrete-time systems. Another way to obtain a larger estimate of 5

the controllable set is the use of a Truncated Ellipsoid (TE) [26] [36] CLF. However, since the standard TE CLF is not differentiable and thus (optimal) nonlinear gradient-based controllers can lead to a non-continuous control signal [30], a smoothing technique has been proposed in [5] via the framework of CLRFs. Note that in order to handle both constraints and local optimality, the aim here is to design a CLF with external level sets in accordance to the maximal controllable set, provided by a PCLF, and with internal level sets close to the quadratic optimal ones, provided by the solution of the ARE. Next section presents the tools of R-functions in order to obtain such kind of smooth CLF having non-homothetic level sets of the chosen shape.

3

On the use of R-functions

The framework of R-functions has been introduced in [32] for geometric applications of logic algebra. The interested reader is referred to [33] for an intensive overview of the theory on the matter. The use of R-functions for state-feedback stabilization has been firstly proposed in [1], [5]. In the following, after recalling the basic notions on R-functions, the useful composition rule introduced in [2] is reported as the basis of our developments. Definition 1 A function r : Fn ⊆ Rn → R is an R-function w.r.t. the arguments ri : Fni ⊆ Rn → R, i = 1, 2, ..., n, if there exists a Boolean function R : Bn → B, where B = {0, 1}, such that h (r) = R (h (r1 ) , h (r2 ) , . . . , h (rn )) ,

(7)

where h(·) is the standard Heaviside step function. Informally, a real function r is an R-function if it can change its sign only when some of its arguments change their sign [5]. In the following, we are concerned with (7) having R equal to the Boolean and and with n = 2. The Boolean operator and can be “recovered” for two real functions r1 , r2 : Rn → R according to the following R-function r∧α : Rn → R. φ



r∧α (x) = ρ(α, φ) φr1 (x) + r2 (x) −

q

φ

(φr1

(x))2

+ r2

(x)2



− 2αφr1 (x)r2 (x) , (8)

+ where α ∈ [0, 1] ⊂ R+ 0 , φ ∈ R and



ρ(α, φ) = φ + 1 −

q

φ2

6

+ 1 − 2αφ

−1

∈ R+

(9)

is the normalization factor such that, for any k ∈ R+ ∀α ∈ [0, 1], ∀φ ∈ R+ .

r1 (x) = r2 (x) = k ⇒ r∧α (x) = k φ

(10)

α

The interpretation of the and composition ∧φ is that the composed function r∧α is positive when evaluated in x if and only if both functions r1 (x) and r2 (x) φ are positive in x. The result is obtained by exploiting the triangle inequality and the law of cosines, according to the following Lemma. Lemma 1 Consider the R-function r∧α (8). Then, ∀x ∈ Rn φ

r∧α (x) > 0 ⇔ {r1 (x) > 0 ∧ r2 (x) > 0} φ

∀α ∈ [0, 1], ∀φ ∈ R+ .

(11)

Remark 1 The composition rule (8) generalizes the original formulation of R-functions [32], that is recovered for φ ≡ 1. The additional parameter φ preserves function r∧α to be a suitable R-function, because sign (ri ) = sign (φri ) φ

∀φ ∈ R+ . Remark 2 Let α = 1. Then r∧α

1

φ

= r1 ∧φ r2 = min {φr1 , r2 } .

(12)

Therefore parameter α links the R-composition and the non-differentiable min operator, used in the synthesis of PCLFs and TE CLFs. 3.1

Using R-functions in the motivating example revisited

Considering the case of study of Section 2 with an uncertainty of ±25% on the input flow parameter F , therefore, for w ∈ [−0.25, 0.25] ∈ R, 

 −1.7428(1 + w)

A = A(w) = 

−4.4191(1 + w)

148.5626 

1 + w

B = B(w) = 

−0.0271

0

0 0.042

  ,

(13)

  .

(14)

The maximal robust controlled polyhedral set is of the kind X = {x ∈ R2 : kGxk∞ ≤ 1} and the procedure of [25] on the Euler Approxi7

398.3

397.3

TR (°K)

396.3

395.3

394.3

393.3

392.3 0.42

0.47

0.52

0.57 CA (kmol/m3)

0.62

0.67

0.72

Fig. 2. Level sets of the polyhedral control Lyapunov function with maximal controlled set.

mating System (EAS) with time-step δt = 0.1, yields the polyhedral function of Figure 2. The provided matrix G has got 33 rows. Two rows of G, say G1 and G2 , describe the state constraints, namely



> > {|x1 | ≤ 0.16, |x2 | ≤ 3} ⇔

[G> 1 , G2 ] x

∞

≤ 1.

(15)

The smoothed polyhedral function [10] (of the second order ) V1 (x) = kGxk22p , p ∈ Z+ , can be composed with the quadratic function V2 (x) = x> P ∗ x, respectively in their 1-level sets N [V1 , 1] and N [V2 , 1], within the framework of R-functions. Define functions R1 (x) = 1 − V1 (x) and R2 (x) = 1 − V2 (x). Without loss of generality, these functions are normalized so that their maximum α α value is 1. Then compute the R-intersection (and rule ∧φ ) R∧α = R1 ∧φ R2 , φ

according to (8), for arbitrary values of α ∈ [0, 1], φ ∈ R+ . The composed function R∧α is the (smoothed) intersection between the smoothed polyheφ dral function and the quadratic one in the sense that R∧α is positive inside φ

the geometric intersection region N [V1 , 1] ∩ N [V2 , 1], it is zero on the boundary, negative outside, and its maximum value is 1 at the origin. Therefore V∧α = 1 − R∧α is the positive definite function associated to R∧α . φ

φ

φ

Since N [V∧α , 1] = N [V1 , 1] ∩ N [V2 , 1], in order to obtain the maximal conφ

trolled set N [V1 , 1], function V2 can be enlarged by scaling matrix P and consequently also the performance cost J, such that N [V1 , 1] ⊆ N [V2 , 1], therefore N [V∧α , 1] = N [V1 , 1] ∩ N [V2 , 1] = N [V1 , 1]. φ

The sublevel sets of the function V∧α are shown in Figure 3, with parameters φ p = 10, α = 0.1, φ = 20. 8

398.3

397.3

TR (°K)

396.3

395.3

394.3

393.3

392.3 0.42

0.47

0.52

0.57 CA (kmol/m3)

0.62

0.67

0.72

Fig. 3. Level sets of the control Lyapunov R-function having maximal controlled set and “close-to-optimal” internal sublevel sets.

Remark 3 The “intersection” of a polyhedral function with an ellipsoidal one has been used as a candidate TE CLF in [26], [36]. Within the framework of R-functions, the TE is recovered as a special case (α ≡ 1, φ ≡ 1). 3.2

Asymmetric controlled set

The framework of R-functions can be also used to shape the sublevel sets of a positive definite function with asymmetric domain. It is in fact sufficient to merge an asymmetric polyhedral function V1 of the kind 2



V1 (x) = max{Gi x} i∈Ir

,

(16)

where Gi is the ith row of matrix G ∈ Rr×n and the linear inequalities Gx ≤ 1 describe the asymmetric polyhedron. For instance, in the motivating example, by merging the PCLF corresponding to the asymmetric state constraints 



I2 , −I2 x ≤

>



0.16 3 −0.08 −3

(17)

and the Riccati-optimal QCLF V2 , one obtains the CLF shown in Figure 4. Remark 4 Asymmetric controlled domain of attractions have been recently investigated in [35] via the use of “barrier CLFs”. The benefit of using the Rcomposition (8) is that the inner sublevel sets can be shaped as desired while preserving the asymmetric domain in the state space.

9

398.3

397.3

TR (°K)

396.3

395.3

394.3

393.3

392.3 0.42

0.47

0.52

0.57 CA (kmol/m3)

0.62

0.67

0.72

Fig. 4. Level sets of the control Lyapunov R-function having asymmetric maximal controlled set and “close-to-optimal” internal sublevel sets.

3.3

Some technical properties of the R-composition

The external shape of the overall region N [V∧α , 1] is not affected by the choice φ

of the parameters α, φ [2]. Function V∧α is differentiable in the set N [V∧α , 1] φ

φ

with non-homothetic level sets for α ∈ [0, 1). Parameter α affects the smoothness of the inner sublevel sets of the composed function, while parameter φ can make the shape of the sublevel sets of the composed function closer to one of the two generating functions, as shown in Figure 3.

Proposition 1 Consider the set S = int ({x ∈ Rn : r1 (x) > 0, r2 (x) > 0}). The composed function r∧α converges pointwise to r2 (r1 ) as parameter φ apφ

proaches infinity (zero). lim r α (x) = r2 (x) ∀α ∈ [0, 1], ∀x ∈ S

(18)

lim r α (x) = r1 (x) ∀α ∈ [0, 1], ∀x ∈ S.

(19)

φ→+∞ ∧φ

φ→0+ ∧φ

Remark 5 If N [V1 , k] ⊆ N [V2 , k] and φ = 0, then R∧α (x) = R1 (x) i.e. φ

V∧α (x) = V1 (x) ∀x ∈ N [V∧α , k] = N [V1 , k]. φ

φ

Proposition 2 The gradient of the composed R-function r∧α (8) can be comφ

10

puted as ∇r∧α (x) = ρ(α, φ) [c1 (x, α, φ)φ∇r1 (x) + c2 (x, α, φ)∇r2 (x)] ,

(20)

φ

where c1 (x, α, φ) = 1 + q c2 (x, α, φ) = 1 + q

3.4

−φr1 (x) + αr2 (x) (φr1 (x))2 + r2 (x)2 − 2αφr1 (x)r2 (x) −r2 (x) + αφr1 (x) (φr1 (x))2 + r2 (x)2 − 2αφr1 (x)r2 (x)

,

(21)

.

(22)

Composing two Lyapunov functions

Theorem 1 Assume that functions Vi : N [Vi , k] → R+ 0 , i = 1, 2, are two ˙ Lyapunov functions with time derivatives Vi (x) ≤ −ηVi (x), i = 1, 2, respectively in their domains N [Vi , k]. Then the R-composed function V∧α is φ

a Lyapunov function with decreasing rate η ∈ R+ in the intersection set N [V∧α , k] = N [V1 , k] ∩ N [V2 , k], ∀α ∈ [0, 1], ∀φ ∈ R+ . φ

Remark 6 Theorem 1 extends the result in [3], because it holds ∀α ∈ [0, 1] and ∀φ ∈ R+ . For ease of notation, in the rest of the paper, the notation ∧ is used in place α of ∧φ .

4

4.1

Control Lyapunov R-functions for uncertain linear systems

Problem statement

Consider constrained uncertain linear systems of the kind x˙ = A(w)x + B(w)u,

(23)

where w ∈ W, A(w) ∈ A = conv {A1 , A2 , ..., As }

(24)

B(w) ∈ B = conv {B1 , B2 , ..., Bs } , subject to the constraints x ∈ X = {y ∈ Rn : kGyk∞ ≤ 1} , 11

(25)

being G ∈ Rr×n full column rank, and u ∈ U = {v ∈ Rm : kvk∞ ≤ 1} .

(26)

The objective is the state-feedback stabilization of (23) via a continuous control, static or horizon-1 optimal, guaranteeing both a large controlled set and locally-optimal performance w.r.t. to the quadratic cost J (5).

4.2

Conditions for stabilizability via control Lyapunov R-functions

In this section, the R-composed function V∧ , built-up by merging a smoothed PCLF V1 and a QCLF V2 , is considered as candidate CLF for solving the control problem of Section 4.1. According to the R-composition procedure shown in Section 3: Ri = k − Vi , i = 1, 2; R∧ = R1 ∧ R2 ; V∧ = k − R∧ .

(27)

Assumption 1 Function V1 (x) = kGxk22p , for p ∈ Z+ sufficiently large, is the smoothed robust PCLF that shapes the desired controlled set for (23), namely N [V1 , k]. Assumption 2 System (23) is quadrically stabilizable and the Riccati-optimal V2 (x) = x> P x, being P a solution of the ARE up to a scaling factor, is a robust QCLF in the set N [V2 , k]. A necessary and sufficient Petersen-like [30] condition, used in [5] in presence of control constraints of the kind (26), can be extended to the case of systems (23) with certain B as shown in Proposition 3. The case of uncertain B can be addressed under the additional “matching” assumption Bi = B(I + ∆i ) [6], by following the approach of [10] (Section V). Therefore, from now on, the notation B stands for the “nominal” input matrix. Proposition 3 Function V∧ : Rn → R+ 0 is a robust control Lyapunov function in S = N [V∧ , k] if and only if there exists η ∈ R+ such that 1 (

max

x∈S\0

max ∇V∧ (x)A(w)x − w∈W

m X

)

| (∇V∧ (x)B)i | + η

kxk22

≤ 0.

(28)

i=1

1

for a nonlinear system x˙ = f (x, w) + g(x)u, V∧ is a CLF if and only if ∃η > 0 such that: maxx∈S\0 maxw∈W ∇V∧ (x)(f (x, w) + g(x)u) + η kxk ≤ 0.

12

Corollary 1 Consider Assumptions 1, 2. Assume there exist a common continuous state feedback control law κ(x) and a decreasing rate η ∈ R+ s.t. ∀i ∈ I2 max ∇Vi (x) (A(w)x + B(w)κ(x)) ≤ −ηVi (x) ∀x ∈ N [Vi , k]. w∈W

(29)

Then the R-composed function V∧ is a control Lyapunov function with decreasing rate η in the set N [V∧ , k] = N [V1 , k]∩N [V2 , k] for system (23), ∀α ∈ [0, 1], ∀φ ∈ R+ . Corollary 2 [(BMI relaxation)] Assume that there exist K ∈ Rm×n , P ∈ Rn×n , P  0, η ∈ R+ and γijk ∈ R+ 0 , for i = 1, ..., s , j, h = 1, ..., r , such that > (Ai + Bi K)> G> h Gh + Gh Gh (Ai + Bi K) 4

− 2ηG> h Gh +

r X



> γijh G> j Gj − Gh Gh



j=1 >

(Ai + Bi K) P + P (Ai + Bi K) 4 −2ηP

∀i ∈ Is , ∀h ∈ Ir (30)

−¯1m ≤ Kv (l) ≤ ¯1m ∀l, h

(31) i

where v (l) are the vertices of the polyhedron L kGxk2∞ , k . Then the smoothed polyhedral CLRF iV∧ is ha CLF, with decreasing rate η, for system (23) in the h i 2 > set N kGxk∞ , k ∩ N x P x, k , ∀α ∈ [0, 1], ∀φ ∈ R+ . Note that, while Proposition 3 is a necessary and sufficient condition for the robust stabilizability via a given candidate V∧ , Corollaries 1, 2 are only sufficient.

4.3

Lyapunov-based nonlinear control

As the R-composed function V∧ can guarantee both constraints fulfillment and locally-optimal shape, the following horizon-1 control can be derived from an approximate Hamilton–Jacobi–Bellman (HJB) approach 2 . u(x) = arg min

υ∈U¯(x)

n

o

∇V∧ (x) (A(0)x + Bυ) + x> Qx + υ > Rυ ,

2

(32)

for a nonlinear system x˙ = f (x, w) + g(x)u, given a robust CLF V∧ , an approximate HJB control is u(x) = arg minυ∈U˜(x) ∇V∧ (x)(f (x, 0) + g(x)u) + L(x, u), where ˜ U(x) = {υ ∈ U : maxw∈W ∇V∧ (x)(f (x, w) + g(x)υ) + η kxk ≤ 0} and L is the Lagrangian of the cost the be minimized.

13

where 



¯ U(x) = υ ∈ U : max ∇V∧ (x) (A(w)x + Bυ) ≤ −η kxk22 . w∈W

(33)

Control (32) is continuous [3] for the nominal case w ≡ 0 because, after the change of variable u 7→ R1/2 u + 21 R−1/2 B > ∇V∧ (x), it derives from a minimal selection control [15] that is known to be continuous. 4.4

Perturbed uncertain linear systems

Consider the presence of a bounded time-variant disturbance d(t) ∈ D ⊂ Rq , namely x˙ = A(w)x + B(w)u + Ed. (34) In the perturbed case, the control objective is to keep the state vector in a large (robust) controlled invariant set until the disturbance is acting, while optimal stabilization is desired if the disturbance vanishes. The proposed CLRF is also suited to solve this kind of control problem for perturbed uncertain systems (34). In fact, on one hand, the maximal controlled invariant set can be approximately computed by passing to the EAS as suggested in [25], in other to guarantee persistent disturbance rejection and robustness to uncertainties. Moreover, on the other hand, in the case of vanishing disturbances, a locally-optimal shape is obtained as the inner sublevel sets of the CLRF are close to the Riccati-optimal quadratic ones.

5

Constructive algorithms

In this section, two constructive algorithms are presented to build-up a suitable CLRF V∧ that is again the merging of CLFs V1 and V2 satisfying Assumptions 1, 2. The two algorithms represent the trade-off between the volume of the achieved controlled invariant set versus optimality. Roughly speaking, Algorithm 1 looks for the largest domain of attraction with an a-priori-fixed “level of optimality”; Algorithm 2 looks for the best “level of optimality” after a-priori-fixing a desired domain of attraction. In both algorithms, CLF V1 , shaping the maximal controlled set N [V1 , k], is fixed. The parameters used are the desired decreasing rate η ∈ R+ , the R+ functions parameters α ∈ [0, 1), φ ∈ R+ 0 , the step tolerance  ∈ R and the scaling factor δ for V2 . 14

5.1

Algorithm 1: priority to optimality

In Algorithm 1, a close-to-optimal shape is a-priori fixed for V∧ . Then, the maximal domain of attraction for V∧ (with such a fixed shape) is searched. At the end, the obtained domain of attraction is at least the one of the optimal CLF V2 . Technically, Algorithm 1 starts from the merging of V1 and V2 with φ  1. Function V2 is initially scaled such that N [V1 , k] ⊂ N [V2 , k], so that N [V∧ , k] = N [V1 , k] ∩ N [V2 , k] = N [V1 , k]. Then N [V2 , k] is progressively reduced, until V∧ is a valid CLF (in the progressively reduced domain N [V∧ , k]). Algorithm 1. (priority to optimality) 1. Initialization Define parameters η, α, φ  1, ; i = 0; δ (0) s.t. N [V1 , k] ⊂ N [δ (0) V2 , k]. 2. Iteration 

(i)

α

(i)



V α = k − (k − V1 ) ∧φ (k − δ V2 ) , ∧φ

(35)

α

according to the R-composition ∧φ (8). 3. Feasibility test (i) if condition (28) is true for V α (x) ∧φ

then Stop: (i) (i) V α is a valid CLF in N [V α , k], with decreasing rate η. ∧φ

∧φ

else δ (i+1) = δ (i) + ; i 7→ i + 1 Go to Step 2.  (i)

(i+1)

Clearly N [V∧ , k] ⊇ N [V∧ , k], therefore, in the worst-case, Algorithm 1 ends when N [V∧ , k] = N [V2 , k] that is a robust controlled set in light of Assumption 2.

5.2

Algorithm 2: priority to the size of the controlled set

In Algorithm 2, the domain of attraction for V∧ is a-priori fixed to be the one of V1 , i.e. N [V∧ , k] = N [V1 , k] (the latter is ad-hoc built in order to be “large”, if not the largest possible). Then, the inner level curves of V∧ are 15

made as-close-as-possible to the optimal ones of V2 , still guaranteeing that V∧ is a valid CLF in N [V∧ , k] = N [V1 , k]. Therefore we start from the merging of V1 and V2 , again initially scaled such that N [V1 , k] ⊂ N [V2 , k], but with shape parameter φ = 0. Therefore, as initial guess V∧ = V1 according to Proposition 1. Then φ is increased so that the level curves of the CLRF V∧ become closer to the optimal ones of V2 , as long as V∧ remains a valid CLF, while still guaranteeing the desired domain of attraction. Algorithm 2. (priority to the size of the controlled set) 1. Initialization Define parameters η, α, ; δ s.t. N [V1 , k] ⊂ N [δV2 , k]; i = 0; φ(0) = 0. 2. Iteration (i)



α



V α = k − (k − V1 ) ∧φ(i) (k − δV2 ) , ∧φ

(36)

α

according to the R-composition ∧φ (8). 3. Feasibility test (i) if condition (28) is false for V α (x) ∧φ

then Stop: (i−1) (i−1) Vα is a valid CLF in N [V α , k], with decreasing rate η. ∧φ

∧φ

else φ(i+1) = φ(i) + ; i 7→ i + 1; Go to Step 2. 

6

Simulations

In this section, two cases of study, both modeling controlled chemical reactions in a CSTR, are simulated with both nominal conditions and randomly-taken model uncertainties w and external disturbances d. Five control algorithms are tested: the standard LQR and the linear RHC (with a 100-steps prediction and control horizon) both based on the nominal system; the horizon-1 Lyapunov-based (32) control with three different CLFs, namely the Riccati-optimal QCLF for the nominal system, the robust smoothed PCLF with maximal controlled invariant set and the CLRF that merges the above two CLFs. 16

398.3 397.3

TR (°K)

396.3 395.3 394.3 393.3 392.2 0.42

0.47

0.52

0.57 CA (kmol/m3)

0.62

0.67

0.72

Fig. 5. Constraint violation induced by the use of the Riccati-optimal QCLF, even in the case of certain dynamics. The solid line between the ellipse and the external box of the state constraints denotes the maximal controllable set.

In the synthesis of the CLRF V∧ , function V1 shaping the maximal robust controlled set N [V1 , k] (computed according to [25]) is merged with the Riccatioptimal QCLF scaled, together with the quadratic performance cost J (5), such that N [V2 , k] ⊃ N [V1 , k]. Parameters η = 10−4 (that makes the left-hand side of (28) strictly negative), α = 0.1, p = 10 are always used in the merging procedures. Based on our numerical experience, the closed-loop performances do not change much as a function of α and p, as far as α < 1 and p is finite. The choice of p trades-off between the smoothness of the external level set of V∧ and the size of the guaranteed domain of attraction N [V∧ , k] (as p < ∞ “smooths the corners”). What mainly affects the closed-loop performances is the inner shape of the CLF. This motivates a fine tuning of the parameter φ only. The LQR and the Lyapunov-based control associated to the Riccati QCLF are control strategies focused on optimality (of the nominal model). The simulation results of these control strategies are not used for comparisons since in both cases the constraints are often violated: this is the case of the nonadmissible controlled state trajectory in Figure 5. On the contrary, the use of the PCLF is suited for robustness and disturbance-rejection arguments. Simulating the nominal dynamics, the RHC (requiring the highest computational effort to be performed in “real-time”) provides the best performances, while in the perturbed case, the presence of a model uncertainty w(t) and an (asymmetric) external disturbance d(t) yields to constraint violations. Therefore a more complicated, and hence computationally demanding, MPC should be used in order to take in account model uncertainties and disturbance rejection. 17

Table 2 Example 1. Average simulation results starting from 100 randomly-taken initial states, normalized with respect to the RHC. Row J stands for the performance cost, while row T stands for the required computational time. RHC

PCLF

CLRF

J

1

1.516

1.013

T

1

0.172

0.086

As expected, only the Lyapunov-based controllers associated to the smoothed PCLF and to the CLRF yield to the constraints fulfillment also in the perturbed cases. The main benefit of the proposed approach is that the horizon-1 CLRF-based controller (32) is characterized by the maximal robust controlled invariant set and by close-to-optimal performances, see Tables 2, 3.

6.1

The motivating example revisited

Algorithm 2 is used to tune parameter φ, still preserving the maximal controlled set. Setting a step  = 1, the algorithm returns φ = 51. In fact, for φ = 52 the state x∗ = (−0.0469, −0.9438)> leads to a violation of the feasibility condition in (28). Table 2 shows the comparison between the above-mentioned control algorithms for the motivating example described in Section 2 in the nominal case. Averaging, the performance of the horizon-1 CLRF approach is only 1% worse than the RHC performance. Moreover, unlike the nominal RHC, in the case of model uncertainty, the horizon-1 CLRF approach still guarantee the constraints fulfilment in the maximal robust controlled set. Figure 6 shows some controlled state trajectories, of the uncertain dynamics, starting from initial states close to the boundary of the polyhedron that shapes the maximal controlled set. Figure 7 compares the control inputs of the three considered control strategies. From Table 2, one can see that the closed-loop performance of the proposed horizon-1 CLRF approach are very close to the ones of a horizon-100 RHC. In fact, Figure 7 shows that the control values taken by the CLRF-based strategy are quite close to the ones of the RHC.

18

398.3 397.3

TR (°K)

396.3 395.3 394.3 393.3 392.3 0.42

0.47

0.52

0.57 3 CA (kmol/m )

0.62

0.67

0.72

Fig. 6. Example 1 (uncertain dynamics). Some controlled state trajectories converging to the origin according to the level sets of the control Lyapunov R-function.

u1 (kmol/m3)

1 CLRF RHC PCLF

0.5

0

−0.5

−1 0

0.2

0.4

t (min)

0.6

0.8

1

1 CLRF RHC PCLF

u2 (kJ/h)

0.5

0

−0.5

−1 0

0.2

0.4

t (min)

0.6

0.8

1

Fig. 7. Example 1 (uncertain dynamics, random initial state). Control inputs trajectories for the three compared feedback strategies. The horizon-1 CLRF-based control is quite close to the standard horizon-100 RHC.

19

6.2

A perturbed three-dimensional Continuous Stirred Tank Reactor

In this subsection we consider a three-dimensional irreversible, first-order reaction taking place in a CSTR. This example is not a direct extension of the model of Section 6.1, but it is taken from [16,27]. The nonlinear state space model is E F0 (c0 − c) − k0 c exp − c˙ = 2 πr h RT   E Uh F0 (T0 − T ) −∆H ˙ + k0 c exp − +2 (Tc − T ) T = 2 πr h ρcp RT rρcp F0 − F h˙ = , πr2 



(37)

where the parameters are available in [27] (Table I). The perturbed linearized model is x˙ = Ax + Bu + Ed, being d ∈ [−8, 8] the external disturbance corresponding to a variation of the inlet flow rate [27], where 



A=

 −0.7489    11.0600  

−0.0034 −0.0007  −0.6704

0



0

−3  −5.426 · 10

B=

    

  , −2.5450   

0

15 · 1.53 · 10−5

1.297

15 · 0.1218

0

15 · (−6.592 · 10−2 )

    ,  





 

 

−5  −1.762 · 10 

E =  7.784 · 10−2  . (38)   

6.592 · 10−2



The state variables x = (c − cs , T − Ts , h − hs )>

(39)

respectively consist of the deviation from the nominal molar concentration c, the reactor temperature T and the tank level h, while the control variables are u = (Tc − Tc,s , F − Fs )> . (40) For for state and control, box constraints are used, namely |x1 | ≤ 0.8, |x2 | ≤ 6, |x3 | ≤ 0.6; kuk∞ ≤ 1. 20

(41)

Table 3 Example 2. Average simulation results, normalized with respect to the RHC. Row J stands for the performance cost, while row T stands for the required computational time. RHC

PCLF

CLRF

J

1

14.937

1.010

T

1

0.095

0.094

Moreover, it is supposed that the element (2, 2) of matrix A is subject to an uncertainty of ±50%, therefore A(w) ∈ conv {A1 , A2 }, where 







 −0.7489 −0.0034 −0.0007 

 −0.7489 −0.0034 −0.0007   





    11.0600 −0.3352 −2.5450  , A2 =  11.0600 −1.0056 −2.5450  . A1 =          0 0 0 0 0 0 (42)

From [27], the weight matrices of the quadratic performance cost J (5) are Q=

h 1 0 0

0 0 0

0 0 1

i

, R = 0.1I2 .

Adopting Algorithm 1 with fixed φ = 100, we give priority to optimality. Condition (28) is satisfied for N [V∧ , 1] = N [V1 , 1]. Table 3 shows the averaged simulation results for the nominal dynamics. The averaged performance of the horizon-1 CLRF strategy is about 1% worse than the horizon-100 RHC performance.

In the perturbed case, since the disturbance d(t), randomly takes values in the interval [−8, 8], the RHC strategy based on the nominal model is not admissible because the hard constraints are not always satisfied. Moreover, if the disturbance is asymmetric (yielding an asymmetric maximal polyhedral domain), for instance in the interval [0, 8], then the RHC is not very efficient even from the performance point of view. We show the simulation results in the case of model uncertainty w(t) and external perturbation d(t), the latter vanishing after 2 minutes. Figures 8, 9 corresponds to a simulation with randomly taken initial condition x0 very close to the boundary of the hard constraints. For this reason, Figure 9 makes clear how, in the first time instants, the proposed CLRF-based strategy has to be “more aggressive” then the standard RHC, in order to be robust to the worst-case uncertainties and disturbances. On the contrary, the standard RHC does not guarantee that the state variable will remain in the admissible region. 21

0.06 CLRF RHC PCLF

∆ c (mol/l)

0.05 0.04 0.03 0.02 0.01 0 0

1

2

3

4 t (min)

5

6

7

8

1 0 CLRF RHC PCLF

∆ T (°K)

−1 −2 −3 −4 −5 −6 0

1

2

3

4 t (min)

5

6

7

8

0.1 0 CLRF RHC PCLF

∆ h (m)

−0.1 −0.2 −0.3 −0.4 −0.5 0

1

2

3

4 t (min)

5

6

7

8

Fig. 8. Example 2 (uncertain and perturbed dynamics, random initial state close to the boundary). State trajectories of the proposed horizon-1 CLRF-based control and of the standard horizon-100 RHC.

Finally, we simulate the case of persistent disturbances. We consider the case in which d(t) is piecewise-constant and it jumps from its maximum value dmax to dmax /8 at the time instant t = 2 min. Clearly, the state variable can not be stabilized (as the disturbance values are unknown to the controllers), but it reaches a neighborhood of the origin. Figure 10 shows that the control inputs still maintain a qualitatively smooth behavior.

22

1 CLRF RHC u1 (°K)

0.5

0

−0.5

−1 0

0.5

1

1.5 t (min)

2

2.5

3

1 CLRF RHC u2 (l/min)

0.5

0

−0.5

−1 0

0.5

1

1.5 t (min)

2

2.5

3

Fig. 9. Example 2 (uncertain and perturbed dynamics, random initial state close to the boundary). Control input trajectories of the proposed horizon-1 CLRF-based control and of the standard horizon-100 RHC. The control values are very close as Algorithm 1 gives priority to optimality. The control inputs of the PCLF-based strategy are not reported because they are significantly, badly, chattering.

7

Conclusion and future work

This paper considers the state-feedback stabilization of constrained uncertain linear systems, with application to the control of benchmark chemical reactions taking place in a Continuous Stirred Tank Reactor. It has been shown that the class of control Lyapunov R-functions is particularly useful because, unlike switching control strategies, both robustness and optimality arguments can be addresses by designing a continuous control associate to a unique control Lyapunov function with maximal controlled set and inner sublevel sets with locally-optimal shape. Moreover, also asymmetric domain of attractions can be handled as well. For these reasons, the use of control Lyapunov R-functions is particularly effective for the feedback control of multivariable chemical processes, where uncertainties, constraints and optimality arguments are to be faced simultaneously. 23

0.2 CLRF RHC

0.15

u1 (°K)

0.1 0.05 0 −0.05 −0.1 0

0.5

1

1.5

2 t (min)

2.5

3

3.5

4

0.7 CLRF RHC

u2 (l/min)

0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.5

1

1.5

2 t (min)

2.5

3

3.5

4

Fig. 10. Example 2 (uncertain dynamics and setpoint disturbance, random initial state close to the boundary). Control input trajectories of the proposed horizon-1 CLRF-based control and of the standard horizon-100 RHC.

As future line of research, it could be interesting to formulate the Hamilton– Jacobi–Bellman control, together with the proposed control Lyapunov Rfunction, in a discrete-time setting. Further investigations will be focused on the control of the nonlinear plant dynamics, by incorporating the nonlinear process model into the controller structure.

Appendix For ease of notation, in the proofs the dependence of functions R˙ 1 (x), R˙ 2 (x), R1 (x), R2 (x), ρ(α, φ), c1 (x, α, φ), c2 (x, α, φ) on their arguments is not explicitly shown.

Proof of Proposition 1

The proof is similar to the one of Proposition 2 in [3], whose result is here extended to all α ∈ [0, 1], besides all φ ∈ R+ . 24

Proof of Proposition 2

Trivial.

Proof of Theorem 1

Passing to the R-functions Ri = k − Vi , the assumption is equivalent to R˙ i ≥ η(k − Ri ), i = 1, 2, ∀x ∈ N [V∧ , k]. Moreover, the thesis is equivalent to R˙ ∧ ≥ η(k − R∧ ) ∀x ∈ N [V∧ , k]. From Proposition 2, h

R˙ ∧ = ρ φc1 R˙ 1 + c2 R˙ 2

i

≥ ηρ [φc1 (k − R1 ) + c2 (k − R2 )] = ηρ [k(φc1 + c2 ) − (φc1 R1 + c2 R2 )] . (43)

First, note that  

ρ(φc1 R1 + c2 R2 ) = ρ φ 1 + q  1 + q



−φR1 + αR2 (φR1

)2

+

R22

− 2αφR1 R2

−R2 + αφR1

(φR1 )2 + R22 − 2αφR1 R2





 R2 



 R1 +

= 

φR1 )2 + R22 − 2αφR1 R2  ρ φR1 + R2 − q = R∧ . (44) (φR1 )2 + R22 − 2αφR1 R2

Therefore, (43) is equivalent to R˙ ∧ ≥ η [k · ρ(φc1 + c2 ) − R∧ ] and it has to be proved that ρ(φc1 + c2 ) ≥ 1. 25

(45)

ρ(φc1 + c2 ) ≥ 1 ⇔ 

1 −φR1 + αR2 φ + φ q √ 2 + φ + 1 − φ + 1 − 2αφ (φR1 )2 + R22 − 2αφR1 R2 

−R2 + αφR1

≥1 ⇔ (φR1 )2 + R22 − 2αφR1 R2 q φ2 R1 + R2 − αφ(R1 + R2 ) φ2 + 1 − 2αφ ≥ q . (46) (φR1 )2 + R22 − 2αφR1 R2

1+ q

By direct substitution in (46), the latter inequality is true for both the extreme cases {R1 = 0, R2 = k} and {R1 = k, R2 = 0}. Therefore, by dividing by R2 and defining S = R1 /R2 , the right-hand-side of (46) can be written as φ2 S + 1 − αφ(1 + S) √ 2 2 . φ S + 1 − 2αφS

(47)

Then q

φ2 + 1 − 2αφ ≥

φ2 S + 1 − αφ(1 + S) √ 2 2 ⇔ φ S + 1 − 2αφS 

(φ2 + 1 − 2αφ)(φ2 S 2 + 1 − 2αφS) − φ2 S + 1 − αφ(1 + S) 2

2

2

≥ 0⇔ 2

φ (1 − α )(S − 1) ≥ 0, (48) that is always true. This concludes the proof, as ρ(φc1 + c2 ) ≥ 1 ⇒ R˙ ∧ ≥ η (k − R∧ ).

Proof of Proposition 3

Condition 3 is the robust version for uncertain linear systems presented in [2]. It is, in fact, directly obtained by considering the constrained control that   > > ˙ minimize V∧ (x) pointwise, i.e. u(x) = −sign B ∇V∧ (x) .

Proof of Corollary 1

The proof is the same of Corollary 1 in [2] as Theorem 1 is here extended to all α ∈ [0, 1], besides all φ ∈ R+ . 26

Proof of Corollary 2

The proof is the same of Theorem 2 in [3] as Theorem 1 is here extended to all α ∈ [0, 1], besides all φ ∈ R+ .

Acknowledgements

The authors would like to thank prof. Franco Blanchini and prof. Gabriele Pannocchia for useful discussions on Lyapunov-based and predictive control.

References

[1] A. Balestrino, A. Caiti, E. Crisostomi, and S. Grammatico. Stabilizability of linear differential inclusions via R-functions. Proc. of the IFAC Symposium on Nonlinear Control Systems, Bologna (Italy), 2010. [2] A. Balestrino, A. Caiti, and S. Grammatico. Constrained stabilization of a continuous stirred tank reactor via smooth control Lyapunov R-functions. Proc. of the IEEE Conference on Decision and Control, Orlando (Florida, USA), 2011. [3] A. Balestrino, A. Caiti, and S. Grammatico. A new class of control Lyapunov functions for constrained stabilization of linear systems. submitted to Automatica (under review), 2011. [4] A. Balestrino, A. Caiti, and S. Grammatico. Stabilizability of constrained uncertain linear systems via smooth control Lyapunov R-functions. Proc. of the IEEE Conference on Decision and Control, Orlando (Florida, USA), 2011. [5] A. Balestrino, E. Crisostomi, S. Grammatico, and A. Caiti. Stabilization of constrained linear systems via smoothed truncated ellipsoids. Proc. of the IFAC World Congress, Milan (Italy), 2011. [6] B.R. Barmish, M. Corless, and G. Leitmann. A new class of stabilizing controllers for uncertain dynamically systems. SIAM Journal on Control and Optimization, 21(2), 1983. [7] A. Bemporad, M. Morari, V. Dua, and E.N. Pistikopoulos. The explicit linear quadratic regulator for constrained systems. Automatica, 38:3–20, 2002. [8] F. Blanchini. Nonquadratic Lyapunov functions for robust control. Automatica, 31(3):451–461, 1995. [9] F. Blanchini and S. Miani. Constrained stabilization via smooth Lyapunov functions. Systems & Control Letters, 35:155–163, 1998.

27

[10] F. Blanchini and S. Miani. A new class of universal Lyapunov functions for the control of uncertain linear systems. IEEE Trans. on Automatic Control, 44(3):641–647, 1999. [11] N.H. El-Farra and P.D. Christofides. Integrating robustness, optimality, and constraints in control of nonlinear processes. Chemical Engineering Science, 56:1841–1868, 2001. [12] N.H. El-Farra and P.D. Christofides. Bounded robust control of constrained multivariabe nonlinear processes. Chemical Engineering Science, 58:3025–3047, 2003. [13] N.H. El-Farra, P. Mhaskar, and P.D. Christofides. Hybrid predictive control of nonlinear systems: method and applications to chemical processes. International Journal of Robust and Nonlinear Control, 14:199–225, 2004. [14] N.H. El-Farra, P. Mhaskar, and P.D. Christofides. Output feedback control of switched nonlinear systems using multiple Lyapunov functions. Systems & Control Letters, 54(12):1163–1182, 2005. [15] R.A. Freeman and P.V. Kokotovi´c. Robust nonlinear control design. Birkh¨auser, Boston, 1996. [16] M.A. Henson and D.E. Seborg. Nonlinear process control. Prentice Hall, Upper Saddle River, 1997. [17] T. Hu and F. Blanchini. Non-conservative matrix inequality conditions for stability/stabilizability of linear differential inclusions. Automatica, 46:190–196, 2010. [18] T. Hu and Z. Lin. Constrained quadratic Lyapunov functions for constrained control systems. IEEE Trans. on Automatic Control, 48(3):440–450, 2003. [19] A. Kojima and M. Morari. LQ control for constrained continuous-time systems. Automatica, 40:1143–1155, 2004. [20] W.H. Kwon and S. Han. Receding Horizon Control. Springer, London, 2005. [21] D. Q. Mayne, M.M. Seron, and S. Rakovi´c. Robust model predictive control of constrained linear systems with bounded disturbances. Automatica, 41:219–224, 2005. [22] D.Q. Mayne, J.B. Rawlings, C.V. Rao, and P.O.M. Scokaert. Constrained model predictive control: stability and optimality. Automatica, 36:789–814, 2000. [23] P. Mhaskar, N.H. El-Farra, and P.D. Christofides. Robust hybrid predictive control of nonlinear systems. Automatica, 41:209–217, 2005. [24] P. Mhaskar, N.H. El-Farra, and P.D. Christofides. Stabilization on nonlinear systems with state and control constraints using Lyapunov-based predictive control. Systems & Control Letters, 55(8):650–659, 2006.

28

[25] S. Miani and C. Savorgnan. MAXIS-G: a software package for computing polyhedral invariant sets for constrained LPV systems. Proc. of the IEEE Conference on Decision and Control, Seville (Spain), 2005. [26] B. O’Dell and E.A. Misawa. Semi-ellipsoidal controlled invariant sets for constrained linear systems. Journal of Dynamic Systems, Measurement and Control, 124:98–103, 2002. [27] G. Pannocchia and J.B. Rawlings. Disturbance models for offset-free modelpredictive control. AIChE Journal, 49(2):426–437, 2003. [28] G. Pannocchia, J.B. Rawlings, and S.J. Wright. Fast, large–scale model predictive control by partial enumeration. Automatica, 43(5):852–860, 2007. [29] G. Pannocchia, S.J. Wright, and J.B. Rawlings. Partial enumeration MPC: Robust stability results and application to an unstable cstr. Journal of Process Control, 21(10):1459–1466, 2011. [30] I. R. Petersen and B. R. Barmish. Control effort considerations in the stabilization of uncertain dynamical systems. Systems and Control Letters, 9:417–422, 1987. [31] J.A. Primbs, V. Nevistic, and J.C. Doyle. A receding horizon generalization of pointwise min-norm controllers. IEEE Trans. on Automatic Control, 45:898– 909, 2000. [32] V.L. Rvachev. Geometric applications of logic algebra (in Russian). Naukova Dumka, 1967. [33] V. Shapiro. Theory of R-functions and applications: a primer. CPA Technical Report CPA88-3, Cornell Programmable Automation, 1988. [34] M. Sznaier, R. Suarez, and J. Cloutier. Suboptimal control of constrained nonlinear systems via receding horizon constrained control lyapunov functions. International Journal of Robust and Nonlinear Control, 13:247–259, 2003. [35] K.P. Tee, S.S. Ge, and E.H. Tay. Barrier Lyapunov functions for the control of output–constrained nonlinear systems. Automatica, 45(4):918–927, 2009. [36] T. Thibodeau, W. Tong, and T. Hu. Set invariance and performance analysis of linear systems via truncated ellipsoids. Automatica, 45:2046–2051, 2009. [37] M.N Zeilinger, C.N. Jones, and M. Morari. Real-time suboptimal model predictive control using a combination of explicit MPC and online optimization. IEEE Trans. on Automatic Control, 56(7):1524–1534, 2011.

29