Nonlinear Programming Properties for Stable and Robust NMPC

Preprints, 5th IFAC Conference on Nonlinear Model Predictive Preprints, Preprints, 5th IFAC IFAC Conference Conference on on Nonlinear Nonlinear Model Model Predictive Predictive Control 5th Preprints, 5th IFAC Conference on Nonlinear Model Predictive Control Available online at www.sciencedirect.com Control September 17-20, 2015. Seville, Spain Control September September 17-20, 17-20, 2015. 2015. Seville, Seville, Spain Spain September 17-20, 2015. Seville, Spain

ScienceDirect

IFAC-PapersOnLine 48-23 (2015) 388–397

Nonlinear Programming Properties for Stable and Nonlinear Programming Properties for Stable and Nonlinear Programming Properties for Stable Nonlinear Programming Properties for Stable and and Robust NMPC Robust NMPC Robust NMPC Robust NMPC

Xue Yang, Devin W. Griffith, Lorenz T. Biegler Xue Yang, Devin Devin W. Griffith, Griffith, Lorenz T. Biegler Xue Xue Yang, Yang, Devin W. W. Griffith, Lorenz Lorenz T. T. Biegler Biegler Chemical Engineering Department, Carnegie Mellon University, Pittsburgh, Chemical Engineering Department, Carnegie Mellon University, Pittsburgh, Chemical Department, Carnegie Mellon PA 15213, USA (e-mail: [email protected]) Chemical Engineering Engineering Department, Carnegie Mellon University, University, Pittsburgh, Pittsburgh, PA 15213, USA (e-mail: [email protected]) PA 15213, USA (e-mail: [email protected]) PA 15213, USA (e-mail: [email protected]) Abstract: We explore properties of nonlinear programming problems (NLPs) that arise in the formulation Abstract: We explore properties of nonlinear programming problems (NLPs) that arise in the formulation Abstract: explore nonlinear programming (NLPs) that arise the of NMPCWe subproblems and showof influence on stabilityproblems and robustness NMPC. that satisfy Abstract: We explore properties properties oftheir nonlinear programming problems (NLPs)of that arise in inNLPs the formulation formulation of NMPC subproblems and show their influence on stability and robustness of NMPC. NLPs that satisfy satisfy of NMPC subproblems and show their influence on stability and robustness of NMPC. NLPs that linear independence constraint qualification (LICQ), second order sufficient conditions (SOSC) andsatisfy strict of NMPC subproblems and show their influence on stability and robustness of NMPC. NLPs that linear independence constraint qualification (LICQ), second order sufficient conditions (SOSC) and strict linear independence constraint qualification second order sufficient conditions (SOSC) and complementarity (SC), have solutions that (LICQ), are continuous and differentiable with perturbations ofstrict the linear independence constraint qualification (LICQ), second order sufficient conditions (SOSC) and strict complementarity (SC), have solutions that are continuous and differentiable with perturbations of the complementarity (SC), have solutions that are continuous and differentiable with perturbations of problem data. As a result, they are important prerequisites for nominal and ISS stability of NMPC complementarity (SC), havethey solutions that are continuous andfordifferentiable with perturbations of the the problem data. As aa result, are important prerequisites nominal and ISS stability of NMPC problem data. As result, they are important prerequisites for nominal and ISS stability of NMPC controllers. Moreover, we show thatimportant ensuring prerequisites these properties possible reformulation of problem data. As a result, they are for is nominal andthrough ISS stability of NMPC controllers. Moreover, we show that ensuring these properties is possible through reformulation of controllers. Moreover, we show that ensuring these properties is possible through reformulation of the NLP subproblem for NMPC, through the addition of 1 penalty and barrier terms. We show how controllers. Moreover, we show that ensuring these properties is possible through reformulation of penalty and barrier terms. We show how the NLP subproblem for NMPC, through the addition of  1 penalty and barrier terms. We show how the NLP subproblem for NMPC, through the addition of  1 these properties also establish ISS of related sensitivity-based NMPC controllers, such as asNMPC and penalty and barrier terms. We show how the NLP subproblem for NMPC, through the addition of  1 these properties also establish ISS of related sensitivity-based NMPC controllers, such as asNMPC and these properties also establish ISS related sensitivity-based NMPC controllers, such as and amsNMPC. Finally, demonstrate impact of our reformulated NLPs on several examples that have these properties alsowe establish ISS of ofthe related sensitivity-based NMPC controllers, such as asNMPC asNMPC and amsNMPC. Finally, we demonstrate the impact of our reformulated NLPs on several examples that have amsNMPC. Finally, we demonstrate the impact of our reformulated NLPs on several examples that have shown nonrobust performance on earlier NMPCofstrategies. amsNMPC. Finally, we demonstrate the impact our reformulated NLPs on several examples that have shown nonrobust performance on earlier NMPC strategies. shown nonrobust performance on NMPC strategies. shown performance on earlier earlier strategies. © 2015,nonrobust IFAC (International Federation of NMPC Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Nonlinear Programming, Constraint Qualifications, Optimality Conditions, Penalty Keywords: Nonlinear Programming, Constraint Qualifications, Optimality Conditions, Penalty Keywords: Nonlinear Functions, NMPC Keywords: Nonlinear Programming, Programming, Constraint Constraint Qualifications, Qualifications, Optimality Optimality Conditions, Conditions, Penalty Penalty Functions, NMPC Functions, Functions, NMPC NMPC

1. INTRODUCTION 1. INTRODUCTION INTRODUCTION 1. 1. INTRODUCTION Model Predictive Control (MPC) is widely accepted in the Model Predictive Predictive Control Control (MPC) (MPC) is is widely widely accepted in in the the Model process industries Control as a generic multivariable controllerinwith Model Predictive (MPC) is widely accepted accepted the process industries as a generic multivariable controller with process industries aa generic multivariable with constraint handling.as More recently, MPC has controller been extended process industries as generic multivariable controller with constraint handling. handling. More More recently, recently, MPC MPC has has been been extended extended constraint to Nonlinear Model More Predictive Control in extended order to constraint handling. recently, MPC(NMPC) has been to Nonlinear Model Predictive Predictive Control (NMPC) in order order to to to Nonlinear Model Control (NMPC) in realize high-performance control Control of highly(NMPC) nonlinearinprocesses. to Nonlinear Model Predictive order to realize high-performance control of highly nonlinear processes. realize high-performance control of nonlinear processes. In particular, NMPC allows incorporation first principle prorealize high-performance control of highly highlyof nonlinear processes. In particular, particular, NMPC allows allows incorporation of first principle principle proIn NMPC incorporation of first process models and integrates with on-line optimization strategies In particular, NMPC allowswith incorporation of first principle process models and and integrates on-line optimization optimization strategies cess models integrates with on-line strategies consistent tasks, including scheduling and cess modelswith and higher-level integrates with on-line optimization strategies consistent with higher-level tasks, including scheduling and consistent with higher-level including scheduling and planning. NMPC for trackingtasks, and so-called “economic” stage consistent with higher-level tasks, including scheduling and planning. NMPC for tracking and so-called “economic” stage planning. NMPC for tracking so-called “economic” stage costs, as well as associated stateand estimation tasks, are reviewed, planning. NMPC for tracking and so-called “economic” stage costs, as as well well as as associated associated state state estimation estimation tasks, tasks, are are reviewed, reviewed, costs, formulated and in state considerable detail in Rawlings and costs, as well asanalyzed associated estimation tasks, are reviewed, formulated and analyzed in considerable considerable detail in Rawlings Rawlings and formulated and analyzed in detail in and Mayne (2009); Mayne et in al.considerable (2000). Due detail to advances described formulated and analyzed in Rawlings and Mayne (2009); (2009); Mayne Mayne et et al. al. (2000). (2000). Due Due to to advances advances described described Mayne in Chen(2009); and Allg¨ oweret(1998); Magni and Scattolinidescribed (2007), Mayne Mayne al. (2000). Due to advances in Chen Chen and Allg¨ Allg¨oower wer (1998); (1998); Magni Magni and and Scattolini Scattolini (2007), (2007), in fundamental robustness of NMPC are in Chen and and stability Allg¨owerand (1998); Magniproperties and Scattolini (2007), fundamental stability and robustness properties of NMPC are fundamental stability and robustness properties of NMPC are well-known, and many of the key issues related to the applicafundamental and stability robustness properties of the NMPC are well-known, manyand of the the key issues issues related to to applicawell-known, and many of key related the applicability and relevance of NMPC areissues well understood. well-known, and many of the key related to theNevertheapplicability and relevance of NMPC are well understood. Neverthebility and relevance of well Nevertheless, an is theare cost andunderstood. reliability of on-line bility andimportant relevancehurdle of NMPC NMPC well Nevertheless, an an important hurdle is the theare cost andunderstood. reliability of of on-line less, important hurdle is cost and reliability on-line computation; lengthy and unreliable optimization calculations less, an important hurdle is the cost and reliability of on-line computation; lengthy and unreliable optimization calculations computation; lengthy and optimization lead to unsuccessful performance. computation; lengthycontroller and unreliable unreliable optimization calculations calculations lead to unsuccessful controller performance. lead unsuccessful controller performance. lead to tostudy unsuccessful performance. This addressescontroller this issue by examining existence and This study addresses this issue by examining examining existence existence and and This study addresses this issue uniqueness properties of nonlinear (NLPs), and and the This study properties addresses this issue by byprograms examining existence uniqueness of nonlinear nonlinear programs (NLPs), and the the uniqueness properties of programs (NLPs), and stability andproperties sensitivity of their solutions. With the and former, uniqueness of nonlinear programs (NLPs), the stability and sensitivity of their solutions. With the former, stability and sensitivity solutions. the former, nominal trajectories can of be their determined for With NMPC subprobstability and sensitivity of their solutions. With the former, nominal trajectories trajectories can can be be determined determined for for NMPC NMPC subprobsubprobnominal lems, andtrajectories robustnesscan of the NMPC controller (through ISS) nominal be determined for NMPC subproblems, and robustness of the NMPC controller (through ISS) lems, and robustness of the NMPC controller (through ISS) are guided by the latter. We also show that the formulation lems,guided and robustness of the controller ISS) are by the the latter. latter. We NMPC also show show that the the(through formulation are guided by We also that formulation of the NMPC subproblem a key on formulation these NLP are guided by the latter. Wehas also showimpact that the of the NMPC subproblem has a key impact on these NLP of the subproblem has key on these properties. Finally, the existence NLPimpact solutions areNLP difof the NMPC NMPC subproblem has aaof on that these properties. Finally, the existence existence of key NLPimpact solutions that areNLP difproperties. Finally, the of NLP solutions that are ferentiable with respect to problem data leads to the developproperties. the existence of NLP thatdevelopare difdifferentiable Finally, with respect respect to problem problem data solutions leads to to the the ferentiable with to data leads development of sensitivity-based NMPC, which reduces on-line ferentiable with respect to problem datagreatly leads to the development of sensitivity-based NMPC, which greatly reduces on-line ment of NMPC, which computation and computational ment of sensitivity-based sensitivity-based NMPC,delay. which greatly greatly reduces reduces on-line on-line computation and computational computational delay. computation and delay. computation and computational delay. The remainder of this section introduces NLP-based strategies The remainder remainder of of this this section section introduces introduces NLP-based NLP-based strategies strategies The for NMPC to set the section stage for the analysis. Section 2 then The remainder of this introduces NLP-based strategies for NMPC to set the stage for the analysis. Section then for NMPC to set the stage for the analysis. Section 222 then provides some background properties for the NLP subproblem for NMPC to set the stage for the analysis. Section then provides some background properties for the NLP subproblem provides provides some some background background properties properties for for the the NLP NLP subproblem subproblem

of NMPC, including sufficient conditions for optimality and of NMPC, NMPC, including including sufficient sufficient conditions conditions for for optimality optimality and and of constraint These conditions properties for areoptimality then usedand in of NMPC, qualifications. including sufficient constraint qualifications. These properties are then used in constraint qualifications. properties are then used in Section 3 to reformulate These the NMPC subproblem in order to constraint qualifications. These properties are then used in Section 3 to reformulate the NMPC subproblem in order to Section 33stability to reformulate the NMPC in order to promote and differentiability ofsubproblem NLP solutions. Section Section to reformulate the NMPC subproblem in order to promote stability stability and and differentiability differentiability of of NLP NLP solutions. solutions. Section Section promote 4promote builds stability on theseand sensitivity properties order to present and differentiability of in NLP solutions. Section 4 builds on these sensitivity properties in order to present and 44 builds on sensitivity properties in to and analyze asNMPC strategy with clipping. Section builds the on these these sensitivity properties in order orderFinally, to present present and analyze the asNMPC strategy with clipping. clipping. Finally, Section the asNMPC strategy with Finally, Section 5analyze demonstrates these results on three well-known examples, analyze the asNMPC strategy with clipping. Finally, Section demonstrates these these results results on on three three well-known well-known examples, examples, 555 demonstrates where previous NMPC strategies to be nonrobust. demonstrates these results on were three shown well-known examples, where previous NMPC NMPC strategies were shown to be be nonrobust. nonrobust. where previous strategies were shown to where previous NMPC strategies were shown to be nonrobust. 1.1 NLP Strategies for NMPC 1.1 NLP NLP Strategies for for NMPC 1.1 1.1 NLP Strategies Strategies for NMPC NMPC Consider the following discrete-time nonlinear dynamic model Consider the the following following discrete-time discrete-time nonlinear nonlinear dynamic model model Consider of the plant uncertainties: Consider thewith following discrete-time nonlinear dynamic dynamic model of the the plant plant with uncertainties: of with uncertainties: of the plant x(kwith + 1) uncertainties: = fˆˆ(x(k), u(k), w(k)) x(k + 1) 1) = = ffˆ(x(k), (x(k), u(k), u(k), w(k)) w(k)) x(k + x(k + 1) = = ffˆ(x(k), (x(k), u(k)) u(k), w(k)) + d(x(k), u(k), w(k)) (1) = ff (x(k), (x(k), u(k)) u(k)) + + d(x(k), d(x(k), u(k), u(k), w(k)) w(k)) (1) = (1) n n n u and f (x(k), u(k)) + d(x(k), (1) where x(k) ∈ ℜnxx=, u(k) ∈ℜ w(k) ∈u(k), ℜnwww(k)) are the plant n u nu and nw are where x(k) x(k) ∈ ∈ℜ ℜnnx ,, u(k) ∈ ℜ w(k) ∈ ℜ the plant where u(k) ∈ ℜ and w(k) ∈ ℜ are the plant n n x u w states, controls and disturbance signals, respectively, defined at where x(k) ∈ ℜand, disturbance u(k) ∈ ℜ signals, and w(k) ∈ ℜ are defined the plant states, controls respectively, at nx +nu → states, controls and disturbance respectively, defined at time steps tk with integers k > signals, 0. The mapping f :ℜ n +n states, controls and disturbance signals, respectively, defined at x +nu → n time steps t with integers k > 0. The mapping f : ℜ x u k with nx with time steps t integers k > 0. The mapping f : ℜ →  n +n k x u ℜ f (0, 0) = 0 represents the nominal model, while the n time steps t with integers k > 0. The mapping f : ℜ →  x k nx with ℜ f (0, 00 represents the nominal model, while the x +n0) u +n= w → ℜ 0) = while the nx with term d : ℜffnn(0, ℜnnxx is usedthe to nominal describe model, modeling errors, ℜ with (0, 0) = 0 represents represents the nominal model, while the x +n u +n w → n +n +n n term d : ℜ ℜ is used to describe modeling errors, x u w x term d : ℜ →  ℜ is used to describe modeling errors, n +n +n n x u w x estimation errors and disturbances. We assume that f (·, ·) and term d : ℜ errors and → disturbances. ℜ is used toWe describe modeling errors, estimation assume that f (·, ·) and estimation and We ·) d(·, ·, ·) areerrors Lipschitz continuous, and that the that noiseff (·, is estimation errors and disturbances. disturbances. We assume assume that (·,w(k) ·) and and d(·, ·, ·) are Lipschitz continuous, and that the noise w(k) is d(·, ·, ·) are Lipschitz continuous, and that the noise w(k) is drawn from a bounded set W . d(·, ·, ·) are Lipschitz continuous, and that the noise w(k) is drawn from a bounded set W . drawn from a bounded set W . drawn from a bounded set W . With this model description, we compute an estimate of the With this model description, we compute an estimate of the With this model compute estimate of current state x(k)description, that can bewe used for ouran With this model description, we compute annonlinear estimate modelof the the current state x(k) that can be used for our nonlinear modelcurrent state x(k) that can be used nonlinear modelbased controller (NMPC), defined byfor theour following nonlinear current state x(k) that can be used for our nonlinear modelbased controller controller (NMPC), (NMPC), defined defined by by the following following nonlinear based programming problem (NLP): based controller (NMPC), defined by the the following nonlinear nonlinear programming problem (NLP): programming problem (NLP): N−1 programming problem (NLP): N−1 ψ(zl , vl ) (2a) JN (p0 ) := min Ψ(zN ) + N−1 ∑ N−1 Ψ(z )+ ψ(z (2a) JJN (p zl ,vl N 0 ) := min ll ,, v ∑ Ψ(z + ψ(z vvll ))) (2a) l=0 N) ∑ zzl ,v l (p00 )) := := min min Ψ(z ) + ψ(z , (2a) JNN (p N l l ,v ∑ l=0 l l l ,vl = f (z , v ) l=0 s.t. zzl+1 l = 0, . . . N − 1 (2b) l l l=0 s.t. zl+1 = f (zl ,, vvl )) ll = 0, .. .. .. N −1 (2b) s.t. = (2b) l , vl ) l = s.t. zzl+1 = pff (z (z = 0, 0, . . . N N− −1 1 (2b) l+1 (2c) z0 = 0 l l p (2c) zz0 = 0 = p (2c) 0 0 p0 vl ∈ U, zN ∈ X f . (2c) zz0l = ∈ X, (2d) X, v ∈ U, z ∈ X . (2d) zzl ∈ N f l ∈ X, v ∈ U, z ∈ X . (2d) N f l l (2d) zl ∈ X, vl ∈ U, zN ∈ X f .

Copyright 2015 IFAC 389 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2015, IFAC (International Federation of Automatic Control) Copyright © 2015 IFAC 389 Copyright 2015 responsibility IFAC 389Control. Peer review© of International Federation of Automatic Copyright ©under 2015 IFAC 389 10.1016/j.ifacol.2015.11.310

2015 IFAC NMPC September 17-20, 2015. Seville, Spain

Xue Yang et al. / IFAC-PapersOnLine 48-23 (2015) 388–397

Here p0 = x(k) is a fixed parameter in the NLP determined by the actual or estimated plant state. We assume that the states and controls are restricted to the domains X and U, respectively. X f is the terminal set with X f ⊂ X. We also assume that N is sufficiently long and Ψ(zN ) is sufficiently large so that zN ∈ X f is always true for the solution of (2). As shown in Pannocchia et al. (2011) and Gr¨une (2013), N sufficiently long and Ψ(zN ) sufficiently large allows X f to be omitted in (2). The set U is compact and contains the origin; the sets X and X f are closed and contain the origin in their interiors. The stage cost is given by ψ(·, ·) : ℜnx +nu → ℜ, while the terminal cost is denoted by Ψ(·) : ℜnx → ℜ; both are assumed to be Lipschitz continuous.

After solution of (2) the control action is extracted from the optimal trajectory {z∗0 , ..., z∗N , v∗0 , ..., v∗N−1 } as u(k) = v∗0 . At the next time step, the plant evolves as in (1) and we shift the time sequence one step forward, k = k + 1, obtain the next state estimate x(k + 1) and solve the next NMPC problem (2). This recursive strategy gives rise to the feedback law, u(k) = κ(x(k)) with κ(·) : ℜnx → ℜnu and system (1) becomes: x(k + 1) = fˆ(x(k), κ(x(k)), w(k))

(3)

2. EXISTENCE PROPERTIES FOR NLP SOLUTIONS We rewrite Problem (2), with x = (z0 , ..., zN , v0 , ..., vN−1 ) and p = x(k) as: min F(x, p), s.t. c(x, p) = 0, g(x, p) ≤ 0. (4) x

An important characterization of the solution of (4) is the concept of a KKT point, which satisfies the Karush-KuhnTucker conditions for (4): Definition 1. (KKT, see Nocedal and Wright (2006)) KKT conditions for Problem (4) are given by: ∗

c(x ) = 0,

∇g j (x∗ , p)T q ≤ 0,

(5)

0 ≤ ν ⊥ g(x ) ≤ 0 ∗

for some multipliers (λ , ν), where x∗ is a KKT point. We also define L = F(x) + c(x)T λ + g(x)T ν as the Lagrange function of (4). A constraint qualification (CQ) is required so that a KKT point is necessary for a local minimizer of (4) (Nocedal and Wright (2006)). For problem (4) the following CQ is widely invoked. Definition 2. (LICQ, Nocedal and Wright (2006)) The linear independence constraint qualification (LICQ) holds at x∗ when the gradient vectors ∇c(x∗ , p) and ∇g j (x∗ , p); j ∈ J where J = { j|g j (x∗ , p) = 0} (6) are linearly independent. LICQ also implies that the multipliers λ , ν are unique. 390

for ν j = 0, j ∈ J.

A stronger second order condition (which is also easier to evaluate) requires positive definiteness of the Hessian in additional directions, given as follows. Definition 4. (SSOSC, Robinson (1980)) The KKT point with multipliers λ and ν is a strict local optimum if the following strong second order sufficient conditions (SSOSC) hold at x∗ : qT ∇xx L(x∗ , λ , ν, p)q > 0 such that

∇ci (x∗ , p)T q = 0, T

∇g j (x , p) q = 0,

Hence, we replace d(x(k), u(k), w(k)) with d(x(k), w(k)) since u(k) = κ(x(k)). We refer to the above strategy as ideal NMPC (iNMPC), where the on-line calculation time is neglected. Correspondingly we denote the feedback law of iNMPC as uid (k) = κ id (x(k)). iNMPC has well-known stability properties (see Magni and Scattolini (2007)). In this study we show conditions on (2) under which its solution, particularly u(k), remains continuous and differentiable with respect to changes in x(k).

∇F(x∗ ) + ∇c(x∗ )λ + ∇g(x∗ )ν = 0

Definition 3. (SOSC, Fiacco (1983)) The KKT point with multipliers λ and ν is a strict local optimum if the following second order sufficient conditions (SOSC) hold at x∗ : qT ∇xx L(x∗ , λ , ν, p)q > 0 for all q = 0 (7) such that ∇ci (x∗ , p)T q = 0, i = 1, .., nc (8) ∇g j (x∗ , p)T q = 0, for ν j > 0, j ∈ J

∗

= f (x(k), κ(x(k))) + d(x(k), κ(x(k)), w(k)) = f (x(k), κ(x(k))) + d(x(k), w(k))

389

for all q = 0

i = 1, .., nc for ν j > 0, j ∈ J.

(9) (10)

Note that by adding

2 (11) x − x∗ W to the objective in (4), the KKT conditions of (4) are unchanged, (x∗ , λ , ν) remains a KKT point. Moreover, by defining matrix Z as a basis of the nullspace of active constraint gradients (in Definition 4) and choosing any positive semi-definite W , with sufficiently large eigenvalues for Z T W Z, SOSC, SSOSC and GSSOSC (given below) can always be satisfied at x∗ .

For solutions of (4) that satisfy LICQ and SSOSC, Lipschitz continuity of x∗ (p) can be guaranteed by the strong regularity theorem due to Robinson (1980), who treated the KKT conditions as generalized equations. In addition, to provide differentiability of x∗ (p) with respect to p, we need to consider the following definition. Definition 5. (Strict Complementarity, Fiacco (1983)) At a KKT point of (4) (x∗ , λ , ν), the strict complementarity condition (SC) holds only if ν j − g j (x∗ , p) > 0 for each j ∈ J. We now state the key property at the solution of (4), which provides the foundation for sensitivity-based NMPC. Theorem 6. (Implicit function theorem applied to (5), Fiacco (1983)) Let x∗ (p) be a KKT point that satisfies (5), and assume that SC, LICQ and SOSC hold at x∗ . Further let the functions F, c, g be at least k + 1 times differentiable in x and k times differentiable in p. Then • x∗ is an isolated minimizer, and the associated multipliers λ and ν are unique. • for p in a neighborhood of p0 the set of active constraints remains unchanged, • for p in a neighborhood of p0 there exists a k times  differentiable function s(p) = x∗ (p)T , λ (p)T , ν(p)T , that corresponds to a locally unique minimum for (4).

More general results on uniform continuity of the solution of (17) can be derived under the following assumptions. Definition 7. (MFCQ, Nocedal and Wright (2006)) For Problem (4), the Mangasarian-Fromovitz constraint qualification (MFCQ) holds at the optimal point x∗ (p) if and only if • ∇c(x∗ , p) is linearly independent.

2015 IFAC NMPC 390 September 17-20, 2015. Seville, Spain

Xue Yang et al. / IFAC-PapersOnLine 48-23 (2015) 388–397

• There exists a vector q such that T

T

∇c(x , p) q = 0, ∇g j (x , p) q < 0 ∗

∗

j ∈ J.

(12)

MFCQ implies that the set of KKT multipliers is a closed convex polytope (Gauvin (1977)). Another useful constraint qualification is given as: Definition 8. (CRCQ, Janin (1984)) For Problem (17), the constant rank constraint qualification (CRCQ) holds at (x∗ , p0 ), when for all subsets of the active constraint gradients ∇g j (x, p) j ∈ J and ∇c(x, p) (13) retain constant rank for all points in a neighborhood of (x∗ , p0 ). Finally, if MFCQ holds at a KKT point but not LICQ, the multipliers λ , ν are no longer unique, and we need a more general second order condition. Definition 9. (GSSOSC, Ralph and Dempe (1995)) The generalized strong second order sufficient condition (GSSOSC) holds at x∗ when the SSOSC holds for all KKT multipliers λ , ν that satisfy the KKT conditions of (4). MFCQ and GSSOSC at KKT points are the weakest conditions under which the perturbed solution of (4) is locally unique (Kojima (1980)). For these more general conditions we cannot expect x∗ (p) to be differentiable (because active sets are nonunique). However, with these conditions and CRCQ, directional derivatives for x∗ (p) can be calculated with a particular QP formulation (Ralph and Dempe (1995)), and this is sufficient to obtain sensitivity updates in an NMPC context.

through addition of a sufficiently large quadratic regularization term (11). These terms are quite compatible with tracking stage costs as well as economic stage costs (J¨aschke et al. (2014)). With MFCQ, CRCQ and GSSOSC satisfied, selecting ρ larger ¯ will drive ξl to zero, where ρ¯ is than a finite threshold, ρ > ρ, the dual norm of the multipliers at the solution of problem (2). If ξl = 0, then the solution of (14) is identical to the solution of problem (2). Therefore, nominal stability properties of (14) are identical to those of (2). Since a solution with ξl > 0 for arbitrarily large values of ρ implies that problem (2) is locally infeasible, we assume that a finite ρ¯ can be found as long as problem (2) is well-posed. This corresponds to the common assumption that there exists a feasible input sequence, which steers the system to the terminal region, i.e., the horizon N is long enough to satisfy the terminal conditions. By eliminating the ξ variables in (14) it is straightforward to show that problem (14) is equivalent to the following. min Ψ(zN ) + ρg+ (zN )1 + zl ,vl

s.t. zl+1 = f (zl , vl ), z0 = x(k), vl ∈ U, l = 1, . . . , N − 1 ( j)

where g+ (zl ) = max(0, g( j) (zl )). Through the redefinition: ψ(zl , vl ) := ψ(zl , vl ) + ρg+ (zl ) and Ψ(zl , vl ) := Ψ(zN ) + ρg+ (zN ), (15) can be rewritten as: JN (x(k)) = min Ψ(zN ) + vl ,zl

To develop a robust problem formulation we relax X and X f with 1 penalty terms. We assume, without loss of generality that X, X f and U can be represented by simple upper and lower bounds, and we write the bounds from X and X f as inequalities g(zl ) ≤ 0 and g(zN ) ≤ 0, respectively. This leads to the following reformulation of (2): zl ,vl ,ξl

(14)

l=0

s.t. zl+1 = f (zl , vl ), z0 = x(k) g(zl ) ≤ ξl , ξl ≥ 0;

vl ∈ U,

∑ (ψ(zl , vl ) + ρξlT 1)

N−1

∑ ψ(zl , vl )

(16)

l=0

s. t. zl+1 = f (zl , vl ), z0 = x(k),

In order to exploit the properties of NLP solutions and sensitivity, we strengthen our assumptions on the dynamic model and objective function, and assume that f (x, u), ψ(x, u) and Ψ(x, u) have Lipschitz continuous second derivatives with respect to their arguments. We also use the arguments from Pannocchia et al. (2011) and Gr¨une (2013), and assume N sufficiently long and Ψ(zN ) sufficiently large, so that zN ∈ X f is always satisfied in the nominal case.

N−1

∑ (ψ(zl , vl ) + ρg+ (zl )1 )

l=0

(15)

3. NMPC PROBLEM REFORMULATION

min Ψ(zN ) + ρξNT 1 +

N−1

l = 0, ..., N − 1

vl ∈ U, l = 1, . . . , N − 1

which replaces problem (2). In the subsequent development we refer to problem (16) as PN (x(k)). Note that the redefined objective function in (15) is no longer differentiable everywhere, but still Lipschitz continuous, which is sufficient for the stability analysis in Section 4.1. Moreover, a number of smoothing strategies have been applied to the redefined objective function, which allow the properties of SSOSC and LICQ to hold for problem (16). These include smoothed max operators (Oliveira and Biegler (1994)), augmented Lagrangian functions (Zavala and Anitescu (2010)) and barrier functions (Wills and Heath (2004)). Consequently, u(k) = κ(x(k)) and JN (x(k)), either from the NLP (14) or the smoothed version of (16), are always Lipschitz continuous in x(k). In the next subsection we present the NLP sensitivity background to show how these properties lead to continuity and differentiability of the solution of PN (x(k)) with respect to disturbances.

l = 0, ..., N

l = 0, ..., N − 1

3.1 Interior-Point NLP Solvers [1, 1, . . . , 1]T .

where ξl is a penalty variable vector and 1 = It is easy to see that the gradients of the equality constraints contain a nonsingular basis matrix, and are linearly independent. Moreover, it is straightforward to show that the MFCQ always holds at the solution of (14) (see J¨aschke et al. (2014)), and because the inequalities are simple bounds, CRCQ is also satisfied. Under these conditions the multipliers of (14) are bounded. As noted above, GSSOSC is straightforward to satisfy 391

To explore continuity and sensitivity properties of problem (14), we rewrite (4) slightly as: min F(x, p), s.t. c(x, p) = 0, g(x, p) + y = 0, y ≥ 0 (17) x

In interior-point solvers, the inequality constraints of problem (17) are handled implicitly by adding barrier terms to the objective function,

2015 IFAC NMPC September 17-20, 2015. Seville, Spain

min F(x, p) − µ

Xue Yang et al. / IFAC-PapersOnLine 48-23 (2015) 388–397

ny

∑ ln(y( j) ),

(18a)

j=1

s.t. c(x, p) = 0, g(x, p) + y = 0

(18b)

where y( j) denotes the jth component of vector y. Solving (18) for the sequence of µ l → 0, with l = 0, 1, 2, . . . , ∞, leads to solution of (17). As shown in Forsgren et al. (2002), convergence of solutions of (18) to (17) has been proved under GSSOSC and MFCQ. For a given barrier parameter value µ, IPOPT (W¨achter and Biegler (2006)) solves the primal-dual optimality conditions of barrier problems (18) directly, ∇F(x, p) + ∇c(x, p)λ + ∇g(x, p)ν = 0 c(x, p) = 0

g(x, p) + y = 0 Yν = µ1,

(19a) (19b) (19c)

where Y = diag(y). To solve this system of nonlinear equations, IPOPT uses an exact Newton method; at the ith Newton iteration, the search direction ∆si = si+1 − si is computed by linearization of the KKT conditions (19), 

Hi Ac i  Ac i T 0   Ag i T I 0 Vi

Ag i 0 0 0

 I  ∆xi  ∇F(xi , p)+Ac i λi +Ag i νi c(xi , p) ∆yi   0    = −  (20) g(xi , p) + y 0  ∆λi Yi ν − µe Yi ∆νi

where Ac i := ∇c(xi , p), Ag i := ∇g(xi , p) and Hi ∈ ℜn×n is the Hessian of the Lagrange function. After solving a sequence of barrier problems for µ → 0, the solver returns the solution s∗,T (p) = [x∗,T y∗,T , λ ∗,T ν ∗,T ] for problem (17). For sensitivity of the NLP solution, we note that problem (17) is parametric in the data p and the optimal primal and dual variables can be treated as implicit functions of p. For a sufficiently small µ l , the KKT conditions (19) of the barrier problem (18) can be expressed as the equations ϕ(s∗ (p), p) = 0 and we denote K∗ (p0 ) as the KKT matrix in (20). To compute approximate neighboring solutions around an already available nominal solution s∗ (p0 ), we invoke the following classical results. Since (18) also satisfies LICQ and SSOSC at the solution of (18), Theorem 6 is also true for (18). Therefore, these results allow the application of the implicit function theorem to (19) at s∗ (p0 ) to yield:  ∂ s∗ ∂ ϕ(s(p), p)  ∗ K (p0 ) =− (21) ∗ ∂p ∂p s (p0 ),p0

K∗ (p

∗ where 0 ) is the KKT matrix (20) evaluated at s (p0 ), and first-order estimates of neighboring solutions are obtained from: ∂ s∗ (p − p0 ) + O(|p − p0 |2 ) (22) s∗ (p) = s∗ (p0 ) + ∂p ∂ s∗ s(p) ˜ = s∗ (p0 ) + (23) (p − p0 ) ∂p where s(p) ˜ is an approximate solution of s∗ (p). From continuity and differentiability of the optimal solution vector, there exists a positive Lipschitz constant Lq such that,

|s(p) ˜ − s∗ (p)| ≤ Lq |p − p0 |2 .

(24) 392

391

and a first order difference approximation ∆s(p) = s(p) ˜ − s∗ (p0 ) can be obtained from: K∗ (p0 )∆s(p) = −ϕ(s∗ (p0 ), p). (25) Also, for the step length τ ∈ [0, 1], the quantity (26) sˆ = s∗ (p0 ) + τ∆s(p) has the property: ∂ s∗ (27) |sˆ − s∗ (p)| ≤ (1 − τ)| | |p − p0 | + Lq |p − p0 |2 . ∂p This result will be used in the next section to analyze the clipping strategy for asNMPC. 3.2 Active Set Changes When the perturbation p − p0 is large enough to induce a change in the active constraint set at the solution of (17), yˆ( j) < 0 or νˆ ( j) < 0 with τ = 1, which violates the KKT conditions. Tracking s∗ (p) becomes nonsmooth, the linearization (25) is no longer valid and Theorem 6 no longer holds. Nevertheless, GSSOSC, MFCQ and CRCQ still lead to Lipschitz continuity of s∗ (p) as well as the directional derivative of x∗ (p) for p− p0 . As developed in Ralph and Dempe (1995), this directional derivative is determined by solving a specialized QP that can be viewed as an extension of (25), with linearized inequality constraints included and strongly active inequalities (with positive multipliers) expressed as equality constraints. Because multipliers at the active set transition are nonunique, a linear program (LP) in dual space must first be solved to maximize a linearized prediction of the Lagrange function. Both Zavala and Anitescu (2010) and J¨aschke et al. (2014) develop pathfollowing algorithms to track x∗ (p) with respect to p. Detailed presentation of these path-following approaches and general underlying concepts that govern NLP sensitivity are described in both studies. On the other hand, while path-following leads to a rigorous treatment of NLP sensitivity, it is more expensive than the simple update (25). A cheaper alternative called “clipping in the first interval” perturbs the solution only up to the active set change, using (26) so that Theorem 6 still holds. Here, feasibility of the soft constrained formulation (16) only requires vl ∈ U, and clipping ensures that the perturbed control variable value, u(k) remains within its bounds. Because the clipping strategy requires no additional computational cost beyond (25), we incorporate this approach within the advanced step NMPC developed and analyzed in the following section. Additional discussion of these strategies can be found in J¨aschke et al. (2014) and Yang and Biegler (2013). 4. ADVANCED-STEP NMPC STRATEGIES Treatment of NMPC problem PN (p0 ) with the above NLP and sensitivity tools corresponds to the off-line and on-line components, respectively, of advanced step strategies. At time tk we use the current estimate x(k) and control u(k) to predict the future state at tk+1 , but we also begin execution of problem (16) at tk with the prediction of the state at tk+1 as the initial condition. At the solution of (16) we obtain the KKT matrix K∗ . Once the state estimate x(k +1) becomes available, we compute a fast backsolve with K∗ to obtain the control action u(k +

2015 IFAC NMPC 392 September 17-20, 2015. Seville, Spain

Xue Yang et al. / IFAC-PapersOnLine 48-23 (2015) 388–397

1). Consequently, the proposed framework for the advancedstep NMPC (asNMPC) strategy, with clipping to ensure u(k + 1) ∈ U, can be summarized as follows. In background, between tk and tk+1 : (1) Use x(k) and u(k) to predict the future state through ¯ + 1) and solve the x(k ¯ + 1) ≡ f (x(k), u(k)). Set p0 = x(k predicted problem PN (x(k ¯ + 1)). (2) At the solution s∗ (p0 ), retain factors of K∗ (p0 ). On-line, at tk+1 : (1) Obtain the state estimate x(k + 1) and set p = x(k + 1). Compute the sensitivity step ∆s(p) from (25). Extract the appropriate element of ∆s(p), called ∆v0 which corresponds to the perturbed control variable at l = 0. • Calculate v˜0 = v∗0 + ∆v0 . Note that elements of v∗0 that are strongly active have zero perturbations. • If v˜0 ∈ U, set τ = 1. • Else, find the largest value of τ ∈ [0, 1] such that v∗0 + τ∆v0 ∈ U. (2) Set u(k +1) = v∗0 +τ∆v0 ∈ U, and return to the background phase with k := k + 1. Stability and robustness properties of the advanced step strategy have been analyzed in Zavala and Biegler (2009) and Biegler et al. (2015). Here we restate and slightly extend these results. 4.1 Nominal and ISS Stability of asNMPC We begin with well-known results on nominal and robust stability for the iNMPC controller uid = κ id (x) (Magni and Scattolini (2007); Keerthi and Gilbert (1988)), which follow from the solution of (16). Assumption 10. (Nominal Stability Assumptions of iNMPC) • The terminal penalty Ψ(·), satisfies Ψ(z) > 0, ∀z ∈ X f \{0}, • There exists a local control law u = κ f (z) defined on X f , such that f (z, κ f (z)) ∈ X f , ∀z ∈ X f , and Ψ( f (z, κ f (z))) − Ψ(z) ≤ −ψ(z, κ f (z)), ∀z ∈ X f . • The optimal stage cost ψ(x, u) = ψ(x, κ(x)) satisfies α p (|x|) ≤ ψ(x, u) ≤ αq (|x|) where α p (·) and αq (·) are K functions.

• The system (1) is ISS in X if there exists a K L function β , and a K function γ such that for all w in the bounded set W , |x(k)| ≤ β (|x(0)|, k) + γ(|w|), ∀ k ≥ 0, ∀x(0) ∈ X (28) • A function V (·) is called an ISS-Lyapunov function for system (1) if there exist a set X, K functions α1 , α2 , α3 and σ such that, ∀x ∈ X and ∀w ∈ W , we have α1 (|x|) ≤ V (x) ≤ α2 (|x|) and V ( fˆ(x, u, w)) − V (x) ≤ −α3 (|x|) + σ (|w|) Moreover, if X is a robustly invariant set for system (1) and V (·) is an ISS-Lyapunov function for this system, then the resulting system is ISS in X (Chen and Allg¨ower (1998); Magni and Scattolini (2007)). Note that for PN (x), X = ℜnx . To analyze the robustness of the asNMPC controller, we recognize that with x(k ¯ + 1) = f (x(k), u(k)), there exists a future mismatch x(k + 1) − x(k ¯ + 1) = d(x(k), w(k)) at the next time ¯ + 1)) and step, giving rise to two different problems PN (x(k PN (x(k + 1)), with optimal costs JN (x(k ¯ + 1)) and JN (x(k + 1)), respectively. To account for this, we define the mismatch term, ε(x(k + 1)) := JN (x(k + 1)) − JN (x(k ¯ + 1)). From Theorem 6, there exists a local positive Lipschitz constant LJ such that ∀x ∈ X, ε(x(k + 1)) ≤ LJ |d(x(k), w(k))|. We make the following assumptions and establish robust stability of the iNMPC controller from the following theorem. Assumption 13. (Robust Stability Assumptions)

• The solution of PN , given by s∗ (p) in problem (17) satisfies LICQ, SOSC and SC. From Theorem 6, the objective function and s∗ (p) are therefore continuous and differentiable with respect to p and the resulting feedback law, derived from s∗ (p), can be represented as u = κ(x). • d(x, w) is Lipschitz with respect to its arguments with Lipschitz constant Ld and |d(x, 0)| ≤ α0 (|x|), where α0 (|x|) is a K function. Theorem 14. (Robust ISS Stability of iNMPC (Theorem 2 in Magni and Scattolini (2007), see also Jiang and Wang (2001)) Under Assumptions 10 and 13 with α0 (|x|) ≤ LηJ α p (|x|) and η ∈ (0, 1), the cost function JN (x) obtained from the solution of PN (x) (16) is an ISS-Lyapunov function and the resulting closed-loop system is ISS stable.

For the analysis of the robustness properties of the asNMPC controller we need to consider the effect of NLP sensitivity errors. The forward simulation x(k ¯ +1) = f (x(k), u(k)) predicts Nominal stability of iNMPC can be paraphrased by the follow- the future state at tk+1 but the plant will evolve with uncertain dynamics generating x(k +1). Moreover, we need to distinguish ing theorem. iNMPC using uid (k) and asNMPC, which generates Theorem 11. (Nominal Stability of iNMPC Rawlings and Mayne between as as (2009)) Consider the moving horizon problem (16) and associ- u (k) = κ (x(k)). To interpret this difference we consider an id extended problem PN+1 (x(k), u(k)): ˆ ated control law u = u , that satisfies Assumption 10. Then, N−1 JN (x) from PN (x) is a Lyapunov function and the closed-loop ˆ ˆ + ∑ ψ(zl , vl ) J(x(k), u(k)) ˆ := min Ψ(zN ) + ψ(x(k), u(k)) system is asymptotically stable. zl ,vl

In the nominal case both asNMPC and iNMPC controllers produce identical control actions uid (k + 1) = uas (k + 1) (because d(x, w) = 0 and ∆s(p) = 0). Nominal stability of asNMPC then follows directly from Theorem 11 (Zavala and Biegler (2009)). For the analysis of robust stability properties of the iNMPC we consider Input-to-State Stability (ISS) (Magni and Scattolini (2007); Jiang and Wang (2001)). Definition 12. (Input-to-State Stability) 393

s.t. zl+1 = f (zl , vl )

l=0

l = 1, . . . N − 1

ˆ vl ∈ U z0 = f (x(k), u(k)),

This problem has an equivalent solution to problem PN (z0 ) ˆ u) and we consider J(x, ˆ as our candidate ISS Lyapunov ˆ function. We define J as (x(k)) := J(x(k), uas (k)), J id (x(k)) := id id ˆ ˆ J(x(k), u (k)), and also J (x(k)) ¯ := J(x(k), ¯ u¯id (k)), where

2015 IFAC NMPC September 17-20, 2015. Seville, Spain

Xue Yang et al. / IFAC-PapersOnLine 48-23 (2015) 388–397

uid (k) and u¯id (k) are determined as u(k) ˆ ∈ U in PN+1 . For the next time step, we define the following residuals as: εs (x(k + 1)) := J id (x(k + 1)) − J id (x(k ¯ + 1)) as

id

εas (x(k + 1)) := J (x(k + 1)) − J (x(k + 1))

(29a) (29b)

where εs accounts for the model mismatch while εas accounts for approximation errors introduced by NLP sensitivity. From (27) and Theorem 6 we have positive Lipschitz constants LJ , Lu , Lsv and Lsq such that ∀x ∈ X, εs (x(k + 1)) ≤ LJ |x(k + 1) − x(k ¯ + 1)| ≤ LJ |d(x(k), w(k))|

393

problems (see reviews in Biegler et al. (2015) and J¨aschke et al. (2014)). To demonstrate the effectiveness of the NLP reformulations for robust NMPC and provide further insight into the NLP properties that lead to ISS stability, we consider three examples from Grimm et al. (2004). All of these examples exhibit nominal stability with NMPC, but NMPC is not robust with respect to data perturbations in these dynamic models. Here we replace (2) by (14), leading to a robust NMPC controller applied to these examples. Because these examples exhibit some nonsmooth elements, we also applied appropriate smooth approximations so that NLP solvers can be applied.

For these examples we show how the NLP (2) becomes nonrobust, because relevant CQs are violated. As seen below, the modified NMPC based on (14) satisfies the relevant CQs and εas (x(k + 1)) ≤ Lu (|uas (k + 1) − uid (k + 1)|) therefore demonstrates robust stability. Moreover, robustness is ≤ LJ ((1 − τ)Lsv + Lsq |d(x(k), w(k))|) · |d(x(k), w(k))|also observed for the sensitivity-based asNMPC controller. By comparing the successive costs J as (x(k)) and J as (x(k + 1)), 5.1 Example 1 - Artstein’s Circles we arrive at an ISS property, slightly extended from the result in Biegler et al. (2015). The first example has two states, one control and the following Theorem 15. (Robust Stability of asNMPC). Given Assumptions definition for x(k + 1) = f (x(k), u(k)) with: 10 and 13 with K functions α0 (|x|) and α p (|x|) that satisfy: LJ (1 + Lsv + Lsq (α0 (|x|) + 2Ld wmax )α0 (|x|) ≤ ηα p (|x|), (30)  T with |w(k)| ≤ wmax and η ∈ (0, 1), then the cost function J as (x) −|x|u + x(1) x(2) f (x, u) = (32) , obtained from the solution of PN+1 with u = uas is an ISS1 + |x|2 u2 − 2x(1) u 1 + |x|2 u2 − 2x(1) u Lyapunov function and the resulting closed-loop system is ISS stable. U = [−1, 1], X = {x : x ∈ R2 , x(1) ≤ c} with c ∈ (0, 1), and X f = εB2 such that X f ⊂ X. For our calculations, we set Proof : We compare the costs J as (x(k)), J as (x(k + 1)) and use c = 0.25 and ε = 0.1. The states of the nominal system always the above mismatch terms in to obtain, lie on the same circle of the family shown in Figure 1. Above the horizontal axis, negative u generates a clockwise movement, J as (x(k + 1)) − J as (x(k)) while positive u generates counter-clockwise movement. The id as id id ¯ + 1)) − J (x(k)) + J (x(k + 1)) − J (x(k ¯ + 1)) = J (x(k directions are reversed below the horizontal axis. Note that larger magnitudes produce larger movements in u, and magni+J as (x(k + 1)) − J id (x(k + 1)) tude of u diminishes asymptotically as the origin is approached. as ≤ −ψ(x(k), u (k)) + εs (x(k + 1)) + εas (x(k + 1)) In fact, the origin is not reachable for the nominal system in any (31) number of time steps. Also note that the solid circle in Figure 1 ≤ −α p (|x(k)|) + εs (x(k + 1)) + εas (x(k + 1)). is a feasible trajectory if the state can be driven from point b’ to The last two inequalities follow by noting that the solution point b” in one time step, since constraints are only checked at of PN+1 at k provides a feasible solution to PN+1 at k + 1, discrete time steps. and from Assumption 10. Substituting the bounds for the error terms leads to: Let ψ(x, u) be the length of the shorter arc between x and the εs (x(k + 1)) + εas (x(k + 1)) ≤ LJ (1 + (1 − τ)Lsv closest reachable point to the origin in one step. Let Ψ(x) be the length of the shorter arc between the state and the origin. +Lsq (|d(x(k), 0)| + Ld |w(k)|))(|d(x(k), 0)| + Ld |w(k)|) Formulating these expressions into a usable equation presents a ≤ LJ (1 + Lsv + Lsq (|d(x(k), 0)| + 2Ld wmax ))|d(x(k), 0)| challenge. Since the system is symmetric about both the x and y axis, we will first reflect the initial condition to the first quadrant +Ld LJ (1 + Lsv + Lsq Ld |w(k)|)|w(k)| before solving an NLP, and use the following cost definitions: ≤ LJ (1 + Lsv + Lsq (α0 (|x|) + 2Ld wmax ))α0 (|x|) ψ(x, u) = +Ld LJ (1 + Lsv + Lsq Ld |w(k)|)|w(k)| |x|(cos−1 δ )x(1) f1 (x, −1) + (x(2) − |x|)( f2 (x, −1) − |x|) ≤ ηα p (|x(k)|) + σ |w(k)| ([(x(1) )2 + (x(2) − |x|)2 ][ f1 (x, −1)2 + ( f2 (x, −1) − |x|)2 ])1/2 where the second inequality follows from τ ≥ 0, and the last (33) inequality follows from (30) and σ (|w|) = Ld LJ (1 + Lsv + Lsq Ld |w(k)|)|w(k)|. The theorem is proved by substituting this −|x|2 (cos−1 δ )(x(2) − |x|) result into (31) to yield: (34) Ψ(x) = |x|((x(1) )2 + (x(2) − |x|)2 )1/2 J as (x(k + 1)) − J as (x(k)) ≤ (η − 1)α p (|x(k)|) + σ |w(k)|.  where δ = 0.9999, for example, to keep the expression from having an unbounded derivative. 5. REFORMULATING NONROBUST NMPC: THREE CASE STUDY EXAMPLES Source of nonrobustness: The nonrobustness in this problem The reformulated NLP (14) has been applied to a number of is caused by the hard state constraint, x(1) ≤ c. It is obvious large-scale process systems, including distillation and reaction that it is possible for the system to be taken outside of X. 394

2015 IFAC NMPC 394 September 17-20, 2015. Seville, Spain

Xue Yang et al. / IFAC-PapersOnLine 48-23 (2015) 388–397

this NLP is infeasible for any staring point in the first quadrant outside the critical circle. After soft constraints are applied, we obtain the optimal solutions shown in blue in Fig. 2(a), where the state goes beyond x1 = c in one step and then converges to the origin.

−0.91

0.7

soft+NMPC soft+asNMPC

soft+NMPC soft+asNMPC −0.92 0.6 −0.93 0.5 −0.94 0.4

u

X2

−0.95

−0.96

0.3

−0.97 0.2 −0.98 0.1 −0.99

0 0

0.05

0.1

0.15

0.2

0.25

X

0.3

0.35

−1 0

5

10

(a) Trajectory of x

Figure 1. Constrained circles (Grimm et al. (2004)) However, it is also interesting to note nonrobustness even when xk ∈ X ∀ k > 0. As explained in Grimm et al. (2004), there exists a critical circle (e.g., the solid circle in Figure √ 1, given by (x(1) )2 + (x(2) − rc )2 = rc2 , where rc := c/ 1 − c2 ), for which no initial condition above the horizontal axis outside the critical circle has a nominal clock-wise trajectory that does not violate the state constraint. So, if the initial condition is perturbed from point a to point e in Figure 1, or if a subsequent point is perturbed from point b to point d, then the c.w. solution becomes infeasible. Therefore, the optimal value function is discontinuous at the critical circle, which violates a key assumption in the robust stability analysis. Furthermore, we can also pinpoint the cause of nonrobustness from an NLP perspective. Starting at an initial point on the critical circle, e.g. x0 = [.0905, .5]T which corresponds to point a in Figure 1, the NLP (2) suffers unbounded multipliers, which implies a violation of MFCQ. This is a violation of the minimum condition for continuity of the value function and therefore leads to nonrobustness. Reformulation: The example has both hard state and terminal constraints. Here, we soften the state constraints and retain the terminal constraint. As a result, the state constraint becomes (1)

xl ≤ c + ξl , l = 0, . . . , N − 1 and we add 1 penalties to the objective function: N−1

N−1

l=0

l=0

∑ ψ(xl , ul ) + Ψ(xN ) + ρ ∑ ξl

(35)

(36)

The controller does not violate the state constraint in the nominal case. However, in the robust case any violation will be minimized to some degree due to the objective function penalty. Furthermore, the multipliers of the reformulated NLP remain bounded, MFCQ is satisfied, and robustness is obtained. Results: We introduce additive disturbances to the state equations (32) with Gaussian noise with zero mean and standard deviation of 0.05. Horizon length is N = 10 and simulation time is 30. We choose c = 0.25, ε = 0.1 and weight on slack variables are ρ = 105 . The system starts from (0.055, 0.51) such that it hits x1 = 0.25 on its trajectory to the origin. We first attempt to solve NMPC problem (2) without soft constraints and are unable to get a feasible solution. It turns out 395

15

20

25

30

time step

1

(b) Trajectory of u

Figure 2. Example 1. Trajectory of x and u. Reformulated NMPC is shown in blue; reformulated asNMPC is shown in red. We also apply advanced-step NMPC with clipping to the reformulated NLP (14). The state and control profiles are shown in red in Fig. 2. We observe that asNMPC leads to nearly the same results as iNMPC, and robustness is preserved. We also observe in Fig. 2(b) that asNMPC takes larger control moves. Note that after time = 15 there are greater differences between the control profiles generated by NMPC and asNMPC. However, since |x| is close to zero, the effect of different controls could be neglected, according to the model (32).

5.2 Example 2 Example 2 from Grimm et al. (2004), also has a two state system, x+ = f (x, u) with one control where: f (x, u) = [x(1) (1 − u), |x|u ]T x := [x(1) , x(2) ]T

(37)

∈ ℜ2

= X. We choose where u ∈ [0, 1] = U and ψ(x, u) = |x|2 , Ψ = 0, the terminal constraint X f = {0} and a length N = 2. Source of nonrobustness: An interesting characteristic of this system is that from any initial condition away from the origin, there exists a unique input sequence, (u1 , u2 ) = (1, 0), that steers the state to the origin in two steps. This gives the intermediate state (0, |x0 |). Therefore, the optimal control calculated at x1 = 0 is zero, and the optimal control in every other case is 1. Thus, the control law and optimal value function may be discontinuous, which leads to nonrobustness. Specifically, LICQ is violated at the point x0 = (0, |x2,0 |T ). To see this, consider the constraints of the NLP:   (1) (1) x1 − x0 (1 − u0 )   (2)   (38) c(x, u) =  x1(1) − u0 |x0 |  = 0.  x (1 − u1 )  1 u1 |x1 | To provide bounded derivatives we use a small positive constant ε to smooth the square root, i.e., (|x1 |2 + ε)1/2 . The Jacobian (2) evaluated at x0 = [0, x0 ]T is:

2015 IFAC NMPC September 17-20, 2015. Seville, Spain

x1 |x1 |

(2)

u(1) x1 |x1 |

 (1) x0 0 −|x0 | 0   (1)  0 −x1   0 |x1 |

also work well when we increase the standard deviation of the noise to 0.5. 0.5

0.4

(39)   10 0 0 0 1 −|x(2) | 0  0 = 1 0 0 0 00 0 0 Note that the first and third constraints are dependent, and LICQ is violated. Reformulation: For this problem, it is possible to either extend the controller horizon or soften the terminal constraint. We will compare the results for both cases. Note that if we add the slack variable ξ and a soft inequality constraint, then the terminal constraint becomes: |xN |2 ≤ ξ (40) with ξ ≥ 0 and an 1 penalty (ρξ ) added to the objective function. The augmented equation system is:   (1) (1) x1 − x0 (1 − u0 )   (2)  x1 − u0 |x0 |    (1) (41) c(x, u, ξ ) = x − x(1) (1 − u1 ) = 0   2 (2) 1  x − u |x |  1

2

1 0  0 1   − 1 0 u  1 ∇c(x, u, ξ )T =   u1 x1(1) u1 x1(2) − |x1 | − |x1 |  0 0  1 0 0 0  0 1 −|x0(2) | 0  = −1 0 0 0  0 0 0 0 0 0 0 0

1

(1)

x0 0 −|x0 | 0 (2) 0 x1

0 0 1

0 0 0

−|x1 |

0 0 0 1 0 0

0 0 0 1 0

0  0 0  0  0 1

0 (1)

u x − 1|x 2| 2

0.6 0.25

0.2

0.4 0.15

0.1

0.2

0.05

0 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0 0

0.5

X1

(a) Trajectory of x, noise=0.05 or 0.5

20

30

40

50

60

70

80

90

100

(b) Trajectory of u, noise=0.05

Figure 3. Example 2. Trajectory of x and u. N = 2 without soft constraint is shown in blue; N = 10 without soft constraint is shown in red; N = 2 with soft constraint is shown in green.

In order to show that applying the advanced step strategy will not deteriorate robustness, we apply asNMPC to the soft constrained problem with N = 2 and original problem with N = 10. With standard deviation of the noise at 0.05, robustness is obtained either with soft constraints or a longer horizon with N = 10. We increase the standard deviation of the noise from 0.05 to 0.5 and observe that the states still converge to the origin, as shown in Fig. 4. Then we apply asNMPC and observe that robustness is not lost when the advanced-step strategy is applied. It can be observed that when the noise is small (0.05), there is no significant difference between performance of iN 0 MPC and asNMPC; however, with a larger level of noise (0.5), the state converges slower with the advanced-step strategy. 0   0   0  1 0.5

N=2,noise=0.05 N=2,soft,noise=0.05 N=2,soft+asNMPC,noise=0.05 N=2,soft,noise=0.5 N=2,soft+asNMPC,noise=0.5

0.4

0.35

0.3

0.3

0.25

0.25

0.2

0.2

0.15

0.15

0.1

0.1

0.05

0.05

0 0

N=2,noise=0.05 N=10,noise=0.05 N=10,asNMPC,noise=0.05 N=10,noise=0.5 N=10,asNMPC,noise=0.5

0.45

0.35

(1)

10

time step

0.4

u x − 1|x 2| 2

N=2,soft

0.8

0.5

1

N=10

0.3

0.45

0

N=2 1

0.35

X2



ξ − |x2 |2

N=2,noise=0.05 N=10,noise=0.05 N=2,soft,noise=0.05 N=10,noise=0.5 N=2,soft,noise=0.5

0.45

u

u

0 1 0

395

X2

1  0  ∇c(x, u)T =  1 − u(1)  (1) (1)

X2



Xue Yang et al. / IFAC-PapersOnLine 48-23 (2015) 388–397

0.05

0.1

0.15

0.2

0.25

X1

0.3

0.35

0.4

0.45

0.5

0 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

X1

(a) N = 2, asNMPC, with soft con- (b) N = 10, asNMPC, without soft straints constraints

(42) This matrix has full row rank, and LICQ is now satisfied. As noted in Grimm et al. (2004), this example can also be solved with an extended control horizon.

Figure 4. Example 2. Trajectory of x. N = 2 without soft constraint is shown in blue; noise standard deviation=0.05 is shown in green or brown; noise standard deviation=0.5 is shown in red.

Results: We introduce additive disturbance vectors to the state equations, with Gaussian noise with zero mean and standard deviation of 0.05. The horizon length is N = 2, simulation time is 100, ρ = 10 and the system starts from (0.5, 0). We first consider the original NLP (2) with N = 2 and without soft constraints. We plot the trajectory of x within the simulation period in Fig. 3(a). Although we use a long simulation time (100), x still does not converge to (0, 0). The corresponding control profile u in Fig. 3(b) stays close to 1, from x2 (k + 1) = |x(k)|, and does not decay to zero. After increasing N to 10, x converges to the origin within only a few steps. On the other hand, with N = 2 and the soft constraint (40), we observe that x also converges within a few steps from a different trajectory. Moreover, from Fig. 3(b) we see that u(k) ∈ [0, 1] with a longer horizon N = 10 or with a soft constraint, and both x1 and x2 could decrease to zero. The reformulation schemes 396

Figure 5. Plot of γ(θx ) in Example 3 (Grimm et al. (2004)) 5.3 Example 3 Example 3 from Grimm et al. (2004) consists of the following state equations:

2015 IFAC NMPC 396 September 17-20, 2015. Seville, Spain

f (x, u) = (1 + |x(1) |)



Xue Yang et al. / IFAC-PapersOnLine 48-23 (2015) 388–397

 sin(u(1) ) sin(u(2) ) + γ(θx )x π cos(u(2) ) 8 cos(u(1) ) (43a)

X = R2 , X f = B2 (43b) π π (43c) U = [0, ] × [0, ] N = 2 4 2 where γ(θx ) is given in Figure 5, x π is the state rotated π8 8 radians clockwise, and θx is the angle between the state and the positive vertical axis, increasing clockwise. Here

We first solve the (smoothed) original NLP problem with N = 2 and without soft constraints. We plot the trajectories of |x| and x within the simulation period. As shown in Fig. 6(a), |x| stalls at |x| = 1 and as shown in Fig. 6(b), x oscillates within a very small region around (0, 1). We then increase N to 10 and observe that x eventually converges to the origin. We also observe that |x| decreases to zero with a different trajectory with N = 2 and soft constraint (44). Finally, we then increase the the standard deviation of the noise to 0.25, and observe convergence to the origin, but with a slower rate and larger values of |x|. 2.5

2.5

N=2,noise=0.01 N=10,noise=0.01 N=2,soft,noise=0.01 N=10,noise=0.25 N=2,soft,noise=0.25

   2   γ(θx ) = min  (θx − π) , 1 π       (2) (2) −x x −x(1) θx (x) = max , 0 + cos−1 2cos−1 |x| |x| |x(1) |

N=2,noise=0.01 N=10,noise=0.01 N=2,soft,noise=0.01 N=10,noise=0.25 N=2,soft,noise=0.25

2

2 1.5

1.5

X2

|x|

1

0.5

1

0

0.5 −0.5

On the other hand, from an NLP perspective we note that KKT conditions do not apply to this problem because of the nondifferentiable absolute value, max and min operators in the model. These terms must be smoothed, e.g., by |x(1) | ≈ ((x(1) )2 + ε)1/2 ,

min(η, ζ ) ≈ 1/2(η + ζ + ((ζ − η)2 + ε)1/2 ),

max(η, ζ ) ≈ 1/2(η + ζ + ((η − ζ )2 + ε)1/2 )

with ε > 0 small. However, for any ε > 0 the terminal region is no longer reachable in two steps. Thus for the smoothed model, softening the terminal constraint and a longer horizon N are necessary.

5

10

15

20

25

−1 −0.2

30

0

0.2

0.4

0.6

time step

0.8

X1

1

1.2

1.4

1.6

1.8

(a) Trajectory of |x|,noise=0.01 or (b) Trajectory of x,noise=0.01 or 0.25 0.25

Figure 6. Example 3. Trajectory of |x| and x. N = 2 without soft constraint is shown in blue; N = 10 without soft constraint is shown in red; N = 2 with soft constraint is shown in green. To show that applying the advanced step strategy will not deteriorate robustness, we apply asNMPC with clipping to the soft constrained problem (14) with N = 2 and the original problem (2) with N = 10. In the first case robustness is obtained by use of soft constraints. With an increased standard deviation from 0.01 to 0.25 we observe that the states still converge to the origin, shown in Fig. 7. When we apply asNMPC, we observe that robustness is not lost, although with the larger level of noise, the convergence rate is slower. Moreover, for the hard constrained case and N = 10, robustness is obtained by use of a long enough horizon. Increasing the standard deviation from 0.01 to 0.25 allows the states to still converge to the origin, as shown in Fig. 8. With asNMPC we also observe robustness in Figs. 8(a) and 8(b). However, with a larger level of noise (0.25), asNMPC converges with a slower rate and larger values of |x|. 2.5

2.5 N=2,noise=0.01 N=2,soft,noise=0.01 N=2,soft+asNMPC,noise=0.01 N=2,soft,noise=0.25 N=2,soft+asNMPC,noise=0.25

N=2,noise=0.01 N=2,soft,noise=0.01 N=2,soft+asNMPC,noise=0.01 N=2,soft,noise=0.25 N=2,soft+asNMPC,noise=0.25

2

2

1.5 1.5

X2

Source of nonrobustness: This example behaves similarly to Example 2. Choosing an initial condition of x0 = [0, 1]T we note that a two step solution can be generated. On the other (1) (2) hand, for any initial condition x0 = [x0 , x0 ]T outside of the unit disk, the first state must be zeroed out in the subsequent (1) step, leaving the intermediate state [0, 1+|x0 |]T so that the unit disk is reached in two steps. When the state is perturbed outside the unit disk, the controller will simply repeat the process of zeroing out the first state at each time step.

0 0

|x|

The terminal cost is Ψ(x) = |x| and we choose ψ(x, u) = |x|2 + |u|2 , a unit circle terminal region X f = B2 and horizon length N = 2.

1

1

0.5

Reformulation: As in Example 2, robustness may be obtained either by lengthening the control horizon or by softening the terminal constraint as: |xN | ≤ 1 + ξN (44) and adding an 1 penalty to the objective function: N−1

∑ ψ(xl , ul ) + Ψ(xN ) + ρξN

0.5

0

0 0

5

10

15

20

25

30

−0.5 −0.5

0

0.5

time step

(a) Trajectory of |x|

1

X

1.5

2

2.5

1

(b) Trajectory of x

Figure 7. Example 3. N = 2, asNMPC, with soft constraints, noise=0.01 or 0.25.

(45)

l=0

This reformulation leads to satisfaction of the KKT conditions along with SC and LICQ. Results: We introduce additive state disturbances as Gaussian noise with zero mean and standard deviation of 0.01. Horizon length is N = 2, simulation time is 30, ρ = 1, and the system starts from (1, 1.5). 397

ACKNOWLEDGEMENTS This work is partially supported by the National Science Foundation Graduate Research Fellowship Program Grant No. DGE1252522. The first author thanks the Department of Chemical Engineering and the second author thanks the Pittsburgh chapter of the ARCS Foundation and the Choctaw Nation of Oklahoma for generous support.

2015 IFAC NMPC September 17-20, 2015. Seville, Spain

2.5

Xue Yang et al. / IFAC-PapersOnLine 48-23 (2015) 388–397

2.5 N=2,noise=0.01 N=10,noise=0.01 N=10,asNMPC,noise=0.01 N=10,noise=0.25 N=10,asNMPC,noise=0.25

N=2,noise=0.01 N=10,noise=0.01 N=10,asNMPC,noise=0.01 N=10,noise=0.25 N=10,asNMPC,noise=0.25

2

2 1.5

1.5

|x|

X2

1

0.5 1

0 0.5 −0.5

0 0

5

10

15

20

25

30

−1 −0.5

0

time step

0.5

1

X

1.5

2

2.5

1

(a) Trajectory of |x|

(b) Trajectory of x

Figure 8. Example 3. N = 10, asNMPC, without soft constraints, noise=0.01 or 0.25. REFERENCES Biegler, L.T., Yang, X., and Fischer, G.G. (2015). Advances in sensitivity-based nonlinear model predictive control and dynamic real-time optimization. Journal of Process Control, accepted for publication. Chen, H. and Allg¨ower, F. (1998). A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. Automatica, 34, 1205–1217. Fiacco, A. (1983). Introduction to Sensitivity and Stability Analysis in Nonlinear Programming. Academic Press, New York. Forsgren, A., Gill, P., and Wright, M. (2002). Interior methods for nonlinear optimization. SIAM Review, 44/4, 525–597. Gauvin, J. (1977). A necessary and sufficient regularity condition to have bounded multipliers in nonconvex programming. Mathematical Programming, 12(1), 136–138. Grimm, G., Messina, M.J., Tuna, S., and Teel, A. (2004). Examples when nonlinear model predictive control is nonrobust. Automatica, 40, 523–533. Gr¨une, L. (2013). Economic receding horizon control without terminal constraints. Automatica, 49, 725–734. Janin, R. (1984). Directional derivative of the marginal function in nonlinear programming. In A. Fiacco (ed.), Sensitivity, Stability and Parametric Analysis, volume 21 of Mathematical Programming Studies, 110–126. Springer Berlin Heidelberg. J¨aschke, J., Yang, X., and Biegler, L.T. (2014). Fast economic model predictive control based on nlp-sensitivities. Journal of Process Control, 24, 1260–1272. Jiang, Z. and Wang, Y. (2001). Input-to-state stability for discrete-time nonlinear systems. Automatica, 37, 857–869. Keerthi, S. and Gilbert, E. (1988). Optimal infinite-horizon feedback laws for general class of constrained discretetime systems: Stability and moving-horizon approximations. IEEE Trans. Auto. Cont., 57, 265–293. Kojima, M. (1980). Strongly stable stationary solutions in nonlinear programming. In S.M. Robinson (ed.), Analysis and Computation of Fixed points. Academic Press, New York. Magni, L. and Scattolini, R. (2007). Robustness and robut design of mpc for nonlinear discrete-time systems. In R. Findeisen, F. Allg¨ower, and L. Biegler (eds.), Assessment and Future Directions of Nonlinear Model Predictive Control, 239–254. Springer, Berlin. Mayne, D., Rawlings, J., Rao, C., and Scokaert, P. (2000). Constrained model predictive control: stability and optimality. Automatica, 36, 789–814. Nocedal, J. and Wright, S. (2006). Numerical Optimization. Operations Research and Financial Engineering. Springer, New York, 2nd edition. 398

397

Oliveira, N.M.C. and Biegler, L.T. (1994). Constraint handling and stability properties of model predictive control. AIChE J., 40, 1138 – 1155. Pannocchia, G., Rawlings, J., and Wright, S. (2011). Conditions under which supboptimal nonlinear mpc is inherently robust. Systems & Control Letters, 60, 747–755. Ralph, D. and Dempe, S. (1995). Directional derivatives of the solution of a parametric nonlinear program. Mathematical Programming, 70(1-3), 159–172. Rawlings, J.B. and Mayne, D.Q. (2009). Model Predictive Control: Theory and Design. Nob Hill Publishing, LLC. Robinson, S.M. (1980). Strongly regular generalized equations. Math. Oper. Res., 5, 43–62. W¨achter, A. and Biegler, L.T. (2006). On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming. Mathematical Programming, 106(1), 25–57. Wills, A.G. and Heath, W.P. (2004). Barrier function based model predictive control. Automatica, 40, 1415 – 1422. Yang, X. and Biegler, L.T. (2013). Advanced-multi-step nonlinear model predictive control. J. Process Control, 23, 1116– 1128. Zavala, V. and Anitescu, M. (2010). Real-time nonlinear optimization as a generalized equation. SIAM J. Control Optim., 48, 5444–5467. Zavala, V. and Biegler, L. (2009). The advanced step nmpc controller: Optimality, stability and robustness. Automatica, 45, 86–93.