L. T. Biegler Carnegie Mellon University and Zhejiang University August, 2013

1. Advanced Algorithms for Nonlinear Model Predictive Control and Dynamic Real-time Optimization L. T. Biegler Carnegie Mellon University and Zhejiang University August, 2013

2. Complex dynamic process models are available for powerful off-line optimization strategies. Can these dynamic models be applied directly for on-line control and real-time optimization?

3. Overview • On-line Dynamic Process Optimization • Three Key Concepts • Representation of Dynamic Problem • Newton Barrier NLP • NLP Sensitivity • Advanced Model Predictive Control • NMPC  asNMPC • MHE  asMHE • RTO  D-RTO • Where to NMPC and NLP? • Conclusions

4. RTO – State of the Art w Plant APC y u RTO c(x, u, p) = 0 DR-PE c(x, u, p) = 0 p • Data Reconciliation & Parameter Identification • Estimation problem formulations • Steady state model • Maximum likelihood objective • functions considered to get • parameters (p) • MinpF(x, y, p, w) • s.t. c(x, u, p, w) = 0 • x  X, p  P • Real-time optimization • Steady state model for states (x) • Supply setpoints (u) to APC • (control system) • Model mismatch, measured and • unmeasured disturbances (w) • Minu F(x, u, w) • s.t. c(x, u, p, w) = 0 • x  X, u  U 9

5. RTO Characteristics w Plant APC y u RTO c(x, u, p) = 0 DR-PE c(x, u, p) = 0 p • Data reconciliation – identify gross errors and consistency in data • Periodic update of process model identification • Usually requires APC loops (MPC, DMC, etc.) • RTO/APC interactions: Assume decomposition of time scales • APC to handle disturbances and fast dynamics • RTO to handle static operations • Typical cycle: 1-2 hours, closed loop • State of art operation - widely implemented in refineries and petrochemical • plants (Marlin and Hrymak, 1997; Perkins, 1998) 10

6. Nonlinear Optimization for Advanced Process Control NMPC Estimation and Control • Why NMPC? • Track a profile – evolve from linear dynamic models (MPC) • Severe nonlinear dynamics (e.g, sign changes in gains) • Operate process over wide range (e.g., startup and shutdown) z : differential states y : algebraic states Process d : disturbances u : manipulated variables NMPC Controller ysp : set points Model Updater NMPC Subproblem

7. Nonlinear Optimization for Advanced Process Control

8. w Plant APC y u RTO c(x, u, p) = 0 DR-PE c(x, u, p) = 0 p Dynamic On-line Optimization: • Goal: Integrate On-line Optimization with APC • Consistent, first-principle dynamic models • Consistent, feed-forward optimization • Increase in computational complexity • Time-critical calculations • Essential for: • Feed changes • Nonstandard operations • Optimal disturbance rejection

9. w Plant m PC u y D-RTO c(x, x’, u, p) = 0 DR-PE c(x, x’, u, p) = 0 p Dynamic On-line Optimization: • Goal: Integrate On-line Optimization with APC • Consistent, first-principle dynamic models • Consistent, feed-forward optimization • Increase in computational complexity • Time-critical calculations • Essential for: • Feed changes • Nonstandard operations • Optimal disturbance rejection

10. Overview • On-line Dynamic Process Optimization • Three Key Concepts • Representation of Dynamic Problem • Newton Barrier NLP • NLP Sensitivity • Advanced Model Predictive Control • NMPC  asNMPC • MHE  asMHE • RTO  D-RTO • Where to NMPC and NLP? • Conclusions

11. Dynamic Optimization - Underlying Problem min F(z(tf)) s.t. dz(t)/dt = f(z(t), y(t), u(t), t, p), z(0) = z0 0 = g(z(t), y(t), u(t), t, p) zl ≤ z(t) ≤ zu yl ≤ y(t) ≤ yu ul ≤ u(t) ≤ uu pl ≤ p ≤ pu t, time z, differential variables y, algebraic variables tf, final time u, control variables p, time independent parameters

12. DAE Models in Process Engineering • Differential Equations • Conservation Laws (Mass, Energy, Momentum) • Algebraic Equations • Constitutive Equations, Equilibrium (physical properties, hydraulics, rate laws) • Semi-explicit form • Assume to be index one (i.e., algebraic variables can be solved uniquely by algebraic equations) • If not, DAE can be reformulated to index one (see Ascher and Petzold) • Characteristics • Large-scale models – not easily scaled • Sparse but no regular structure • Direct linear solvers widely used • Coarse-grained decomposition of linear algebra

13. DAE Optimization Problem Apply a NLP solver Efficient forconstrained problems Multiple Shooting Embeds DAE Solvers/Sensitivity Dynamic Optimization Approaches Pontryagin(1962) Indirect/Variational Discretize controls Hasdorff (1977), Sullivan (1977), Vassiliadis (1994)… Single Shooting Bock and coworkers Discretize state, control variables Simultaneous Approach Larger NLP Handles instabilities • Take Full Advantage of Open Structure • Many Degrees of Freedom • Periodic Boundary Conditions • Multi-stage Formulations • … Simultaneous Collocation (Direct Transcription) Large/Sparse NLP - Betts; B…

14. RadauCollocationonFiniteElements Polynomials · · · · · · · · · · · · · ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ ´ tf ti to Collocation points hi Piecewise Continuous Algebraic and Control Variables Mesh points Finite element, i ´ Continuous Differential Variables ´ ´ ´ ´ ´ ´ element i q = 2 ´ ´ ´ q = 1 ´ ´ Monotonic polynomials for u(t) (e.g., piecewise linear)

15. Collocation on finite Elements Nonlinear Programming Problem (NLP) Discretized variables NonlinearProgrammingFormulation Nonlinear Dynamic Optimization Problem (Piecewise) Continuous profiles

16. Nonlinear Programming Problem s.t.

17. How Process Optimization?Optimization Formulations and Models Open First and Second Derivatives NLP Barrier Exact First Derivatives Compute Efficiency rSQP Finite Differences SQP DFO Black Box Closed 102 104 106 100 Variables/Constraints

18. Barrier Methods for Large-Scale Nonlinear Programming Can generalize for Original Formulation a ≤ x ≤ b Barrier Approach • As mè 0, x*(m) è x* Fiacco and McCormick (1968)

19. Barrier Problem Solution • Newton Directions (KKT System) • Reducing the System • Solve IPOPT Code – www.coin-or.org

20. IPOPT Features(Wächter, Laird, B., 2002-2009) • Hessian Calculation • - BFGS (full/LM and reduced space) • - SR1 (full/LM and reduced space) • - Exact full Hessian (direct) • - Exact reduced Hessian (direct) • - Preconditioned CG • Widely Available • Eclipse License and COIN-OR distribution: http://www.coin-or.org • Solved on many thousands of test problems and applications • Broad, growing user community Newton-Based Barrier Method • Globally, superlinearly convergent (Wächter and B., 2005) • Easily tailored to different problem structures Line Search Globalization • l2exact penalty merit function • augmented Lagrangian merit function • Filter method (extended from Fletcher and Leyffer)

21. Case Study: High Purity Distillation Column(Scalable Process Model - trays, comps, complexity) Liquid / vapor Product stream V L D F • Methanol, n-Propanol Separation • Detailed dynamic MESH model • Use Reflux, Reboiler duty to maintain purity • Infer composition with temperature measurements • Setpoints on Trays 14 and 28 V’ L’ B

22. Case Study: High Purity Column (Diehl, 2001) • MESH Equations for Distillation Column Assumption: Vapor holdups are negligible. Ideal vapor phases. Well mixed entering streams. Constant pressure drops Constant tray efficiency Li-1 Vi Fi Mi Li Vi+1 Mass balance: Comp. balance: Energy balance: Phase equilibrium: Index 2 DAE System • Reformulated to index 1 Summation: Hydrodynamics:

23. Distillation Model Structure • Block Tridiagonal Structure • (2 Nc + 3) eqns + and thermo fcns per stage • Nt rows • Scaleable – Nt: 10 – 200, Nc: 10 – 1000 • Nonlinear extensions in stage eqns. • azeotropic behavior (singular points) • mass transfer extension (MERQ eqns.) • tray hydraulics …

24. Dynamic Optimization Distillation Results (Setpoint Change) Tray 14 Temperature Tray 28 Temperature Reboiler duty Reflux rate • 40 trays, 252 DAEs  19814 vars/ 19794 eqns • NLP Solution Time = 9.4 CPUs, KKT tolerance = 10-8

25. Other Dynamic Optimization Case Studies • Off-line Case Studies • Dynamic Bioprocess Optimization • Parameter Estimation of Batch Data • Synthesis of Reactor Networks • Crystallization Temperature Profiles • Optimal Batch Distillation Operation • Satellite Trajectories in Astronautics • Batch Process Integration • Simulated Moving Bed Optimization • Optimal Pressure Swing Adsorption • Optimization of Polymerization Processes • Profiles for Continuous Distillation • On-line Case Studies • NMPC of Tenn. Eastman Process • Source Detection of Water Networks • Cross-directional Sheet-forming Processes • Thermo-mech. Pulping NMPC • D-RTO for Gas Pipelines • Air Traffic Conflict Resolution • NMPC for Refinery Distillation • Startup for Combined Cycle Power Generation • Cyclic Operation for LDPE Process • Ramping for Air Separation Columns

26. Overview • On-line Dynamic Process Optimization • Three Key Concepts • Representation of Dynamic Problem • Newton Barrier NLP • NLP Sensitivity • Advanced Model Predictive Control • NMPC  asNMPC • MHE  asMHE • RTO  D-RTO • Where to NMPC and NLP? • Conclusions

27. On-line Issues: Model Predictive Control (NMPC) NMPC Estimation and Control • Why NMPC? • Track a profile – evolve from linear dynamic models (MPC) • Severe nonlinear dynamics (e.g, sign changes in gains) • Operate process over wide range (e.g., startup and shutdown) z : differential states y : algebraic states Process d : disturbances u : manipulated variables NMPC Controller ysp : set points Model Updater NMPC Subproblem

28. MPC - Background Motivate: embed dynamic model in moving horizon framework to drive process to desired state (Rawlings and Mayne, 2009) Generic MIMO controller Direct handling of input and output constraints Slow process time-scales – consistent with dynamic operating policies Different types Linear Models: Step Response (DMC) and State-space Empirical Models: Neural Nets, Volterra Series Hybrid Models: (QP/MIQP…), apply parametric programming and Explicit MPC First Principle Models – direct link to off-line planning NMPC Pros and Cons + Operate process over wide range (e.g., startup and shutdown) + Vehicle for Dynamic Real-time Optimization - Need Fast NLP Solver for time-critical, on-line optimization - Computational Delay from On-line Optimization degrades performance

29. What about Fast NMPC? • Fast NMPC is not just NMPC with a fast solver (Engell, 2007) • Computational delay – between receipt of process measurement and injection of control, determined by cost of dynamic optimization • Leads to loss of performance and stability(see Rawlings and Mayne, 2009; Findeisenand Allgöwer, 2004; Santos et al., 2001) • Can computational delay be overcome? • Fast Newton-based NMPC • Cheap NLP Sensitivity

30. What about Fast NMPC? • Fast NMPC is not just NMPC with a fast solver (Engell, 2007) • Computational delay – between receipt of process measurement and injection of control, determined by cost of dynamic optimization • Leads to loss of performance and stability(see Rawlings and Mayne, 2009; Findeisenand Allgöwer, 2004; Santos et al., 2001) • Can computational delay be overcome? • Fast Newton-based NMPC • Cheap NLP Sensitivity

31. Can we avoid on-line optimization? • Divide Dynamic Optimization Problem: • preparation, feedback response and transition stages • solve complete NLP in background (‘between’ sampling times) as part of preparation and transition stages • solve perturbed problem on-line based on NLP sensitivity • > two orders of magnitude reduction in on-line computation • Based on NLP sensitivity of z0 for dynamic systems • Extended to Collocation approach – Zavala et al. (2008, 2009) • Similar approach for Moving Horizon Estimation – Zavala et al. (2008) • Stability Properties (Zavala et al., 2009) • Nominal stability – no disturbances nor model mismatch • Lyapunov-based analysis for NMPC • Robust stability – some degree of mismatch • Input to State Stability (ISS) from Magni et al. (2005) • Extension to economic objective functions

32. Solution Triplet Optimality Conditions NLP Sensitivity Parametric Program NLP Sensitivity  Rely upon Existence and Differentiability of Path  Main Idea: Obtain and find by Taylor Series Expansion

33. NLP Sensitivity with IPOPT (Pirnay, Lopez Negrete, B., 2011) Apply Implicit Function Theorem to around Obtaining Optimality Conditions of KKT Matrix IPOPT  Already Factored at Solution  Sensitivity Calculation from Single Backsolve  Approximate Solution Retains Active Set

34. Advanced Step Nonlinear MPC (Zavala, B., 2009) Solve NLP in background (between steps, not on-line) Update using sensitivity on-line x(k) xk+1|k u(k) tk tk+1 tk+2 tk+N Solve NLP(k) in background (between tk and tk+1)

35. Advanced Step Nonlinear MPC (Zavala, B., 2009) Solve NLP in background (between steps, not on-line) Update using sensitivity on-line xk+1|k x(k) x(k+1) u(k+1) u(k) tk tk+1 tk+2 tk+N Solve NLP(k) in background (between tk and tk+1) Sensitivity to update problem on-line to get (u(k+1))

36. Advanced Step Nonlinear MPC (Zavala, B., 2009) Solve NLP in background (between steps, not on-line) Update using sensitivity on-line xk+2|k+1 x(k) x(k+1) u(k+1) u(k) tk tk+1 tk+2 tk+N Solve NLP(k) in background (between tk and tk+1) Sensitivity to update problem on-line to get (u(k+1)) Solve NLP(k+1) in background (between tk+1 and tk+2)

37. NMPC for High Purity Distillation • Air Separation Unit in IGCC-based Power Plants • Need for high purity O2 • Respond quickly to changes in process demand • Large, highly nonlinear dynamic separation (MESH) models • Methanol distillation (Diehl, Bock et al., 2005) • 40 trays, 210 DAEs, 19746 discretized equations • Argon Recovery Column • 50 trays, 260 DAEs, 21306 discretized equations • Double Column ASU Case Study • 80 trays, 1520 DAEs, 116,900 discretized equations

38. Nonlinear Model Predictive Control Air Separation Unit (Huang, B., 2009) Objective: force the production rates to follow the set-points, while main their purities. 4 manipulated variables. 4 output variables. Horizon: 100 minutes in 20 finite elements. Sampling time: 5 minutes. Reformulated index 1 system 320 ODEs, 1200 AEs.

39. ASU Nonlinear MPC Comparison with Linear MPC t = 30-60 min, product rates are ramped down by 30%. t =1000-1030 min, they are ramped back. Output Variables Manipulated Variables Set-points, linear MPC, AS-NMPC profile. No noise nor model mismatch All tuning parameters determined for linear MPC

40. NMPC of Air Separation Unit (w/noise) Comparision of NMPC with asNMPC At t = 30-60 min, product rates are ramped down by 40%. At t =1000-1030 min, they are ramped back. 5% disturbance is added to Mi. N = 20, K = 3 320 ODEs, 1200 AEs. Variables: 117,140 Constraints: 116,900 400 NLPs solved Background: 200 CPUs, 6 itn. Online: 1 CPUs Computational Delay Reduced from 200  1 CPUs. Blue dashed lines are ideal NMPC profile Red lines are AS-NMPC profile. Linear MPC is unstable for this example

41. Combining MHE & NMPC (Huang, Patwardhan, B., 2009) • Offset-free Formulation • Apply MHE results as state and output corrections for NMPC problem • Modify with an advanced step approach  as-MHE z : differential states y : algebraic states Process d : disturbances u : manipulated variables NMPC Controller ysp : set points Model Estimator

42. Combining MHE & NMPC (Huang, Patwardhan, B., 2009) • Offset-free Formulation • Apply MHE results as state and output corrections for NMPC problem • Modify with an advanced step approach  as-MHE

43. Advanced-step MHE (Zavala, Lopez Negrete, B. 2009 - 2011) Measured outputs Estimated states Background: At tk, having xkanduk, approximatexk+1and yk+1. Solve the extended MHE problem from k-N to k+1. Let p0 = approximate yk +1.

44. Advanced-step MHE (Zavala, Lopez Negrete, B. 2009 - 2011) Measured outputs Estimated states Background: At tk, having xkanduk, approximatexk+1and yk+1. Solve the extended MHE problem from k-N to k+1. Let p0 = approximate yk +1. On-line update:At tk+1, obtain yk+1. Set p = yk+1and use NLP sensitivity to get fast update xk+1.

45. Advanced-step MHE (Zavala, Lopez Negrete, B. 2009 - 2011) Measured outputs Estimated states Background: At tk, having xkanduk, approximatexk+1and yk+1. Solve the extended MHE problem from k-N to k+1. Let p0 = approximate yk +1. On-line update:At tk+1, obtain yk+1. Set p = yk+1and use NLP sensitivity to get fast update xk+1. Iterate:Set k = k+1 and go to background. NLP Sensitivity used for State Approximation and Covariance Updates

46. Combined MHE and NMPC for ASU (Huang, B., 2010) • Change in model mismatch • (hydraulic parameter) • Setpoint • MHE/NMPC with correction • NMPC without correction • AMPL/IPOPT (DuoCore 2.4 GHz) • NMPC background (6 iter/200 CPUs) • MHE background (15 iter/90 CPUs) • N=5, K=3 • 29,285 constraints • 30,885 variables • Online NMPC+MHE ~ 2 CPUs

47. Recent asNMPC Studies • NMPC for CSTRs (Zavala, B., 2009) • NMPC and D-RTO for LDPE reactors (Zavala, B., 2009) • NMPC and D-RTO for air separation units (Huang, Zavala, B., 2009; Huang, B., 2010) • NMPC for Thermo-mechanical pulping (Harinath, B., Dumont, 2010) • NMPC for C3 splitters (Fischer, B., 2011) • NMPC for power generation cycles (d’Amato, Kumar, Lopez-Negrete, B., 2011) • Extension of asNMPC over multiple time steps (Yang, B., 2013)

48. NMPC with Economic ObjectivesBeyond RTO and MPC Regulation  D-RTO • Benefits of combining RTO with NMPC? • Direct, dynamic production maximization • Remove artificial setpointobjective • Remove artificial steady state problem • Overcomes neglect of dynamic uncertainty • Leads to significant improvements (up to • 10%) over steady state RTO w Plant APC y u RTO c(x, u, p) = 0 DR-PE c(x, u, p) = 0 Plant m p PC u y D-RTO c(x, x’, u, p) = 0 DR-PE c(x, x’, u, p) = 0 p

49. Challenges with D-RTO Replace regulation objective with economic objective in NMPC? Bartusiak, Young et al. (2007) Chachuat et al. (2008), Dadhe and Engell (2008), Engell(2007, 2009) Busch, Kadam Marquardt et al. (2008) Odloak, Zanin, Tvrzska de Gouvea (2002) Zanin, Tvrzska de GouveaOdloak (2000) Diehl, Amrit and Rawlings (2010) Angeli and Rawlings (2010) Angeli, Amrit and Rawlings (2011) Robust Stability of Lyapunov function must be K∞ function (e.g., strong convexity of stage cost) Many open Stability/RobustnessQuestions Still Remain - does optimum go to a steady state or not? - how do we enforce optimal steady state? - how to consider cyclic problems? Remedy:Regularize economic objective with K∞ function for stage cost?

50. Challenges with D-RTO Replace regulation objective with economic objective in NMPC? Bartusiak, Young et al. (2007) Chachuat et al. (2008), Dadhe and Engell (2008), Engell(2007, 2009) Busch, Kadam Marquardt et al. (2008) Odloak, Zanin, Tvrzska de Gouvea (2002) Zanin, Tvrzska de GouveaOdloak (2000) Diehl, Amrit and Rawlings (2010) Angeli and Rawlings (2010) Angeli, Amrit and Rawlings (2011) Robust Stability of Lyapunov function must be K∞ function (e.g., strong convexity of stage cost) Many open Stability/RobustnessQuestions Still Remain - does optimum go to a steady state or not? - how do we enforce optimal steady state? - how to consider cyclic problems? Remedy:Regularize economic objective with K∞ function for stage cost?