Model Predictive Control: Background B. Wayne Bequette Concise - PowerPoint PPT Presentation

Dynamic Behavior: SISO Relative order = n-m m m 1 b s b s  b s b − + + + + Polynomial ( ) m m 1 1 0 g s − = p n n 1 a s a s  a s a − + + + + n m 1 1 0 − ( )( ) ( ) k s z s z  s z − − − pz 1 2 m Pole-zero ( ) g s = p ( )( ) ( ) s p s p  s p − − − 1 2 n ( )( ) ( ) k s 1 s 1  s 1 τ + τ + τ + Gain-time constant p n 1 n 2 nm ( ) g s = ( )( ) ( ) p s 1 s 1  s 1 τ + τ + τ + p 1 p 2 pn Zeros = roots of numerator Poles = roots of denominator = determine stability 1 p Large time constant = small, negative, pole τ = − p 1 1

Example of Pole-Zero Cancellation l Number of states = Number of poles (order of numerator of transfer function), except when some poles and zeros ‘cancel’ 0 0 . 9056 1 . 5301 ⎡− ⎡ ⎤ ⎤ A B Bioreactor = = ⎢ ⎥ ⎢ ⎥ 0 . 75 2 . 5640 3 . 8255 − − ⎣ ⎦ ⎣ ⎦ model [ ] [ ] C 1 0 D 0 = = ( ) 1 . 5302 s 0 . 4590 1 . 5302 s 0 . 3 − − − + ( ) g p s = = 2 ( )( ) s 2 . 564 s 0 . 6792 s 0 . 3 s 2 . 2640 + + + + − 0.6758 ( ) = g p s 0.4417 s + 1 “hidden slide” provided for additional background

Chemical Process Engineers l More familiar with gain-time constant form l Most chemical processes are stable Ø Exceptions: Exothermic or bioreactors, closed-loop systems (mistuned) Imaginary axis more oscillatory x x = pole inverse response zeros faster response Real axis o x o = zero unstable poles x “hidden slide” provided for additional background

First-Order y + k p dy dt = − 1 dy or u dt + y = k p u τ p τ p τ p gain k p ( ) ( ) y s u s = τ s 1 + 20 p y, deg C 15 10 time constant output 5 0 0 2 4 6 8 10 12 14 16 18 20 t, minutes 10 Example: gain = 1.43 o C/kW 8 u, kW 6 input time constant = 5 minutes 4 2 0 0 2 4 6 8 10 12 14 16 18 20 t, minutes

First Order + Time-delay gain time-delay s k e − θ dy p ( ) ( ) ( ) y s u s y k u t = τ + = − θ p p s 1 dt τ + p time constant long-term 63.2 % y change in of change process output in process output time- time delay constant change in u process input

Second-order Underdamped k 1 ( ) = ( ) y s 2 + 2 ζτ s + 1 ⋅ u s ζ < 2 s τ time to first peak 2 1 ζ − ζ overshoot ratio = a/c p = − ± a 1 , 2 τ τ decay ratio = b/a b p Re j Im = ± 1 , 2 period of oscillation rise time c step response time

Numerator Dynamics ( ) s 1 large faster response τ + τ = → ( ) g s n n = p ( )( ) 1.2 3 s 1 15 s 1 + + 20 1 z 1 = − τ 15 0.8 n 0.6 10 0.4 y 0.2 3 0 0 z 0 τ < → > -3 n -0.2 “inverse response” -10 various values of numerator time constant -0.4 -15 (right-half-plane zeros) -0.6 0 10 20 30 40 50 60 time [challenging control problem]

Steady-state Gain l For stable systems, the steady-state gain is found from the long-term response Step ( )( ) ( ) k s 1 s 1  s 1 u τ + τ + τ + Δ p n 1 n 2 nm ( ) y s = ⋅ input ( )( ) ( ) s 1 s 1  s 1 s τ + τ + τ + p 1 p 2 pn Final value theorem ( )( ) ( ) k s 1 s 1  s 1 u τ + τ + τ + Δ p n 1 n 2 nm ( ) ( ) lim y t lim sy s s = = ⋅ ⋅ ( )( ) ( )  s 1 s 1 s 1 s t s 0 τ + τ + τ + → ∞ → p 1 p 2 pn u Δ ( ) ( ) lim y t lim sy s s k k u = = ⋅ ⋅ = Δ p p s t s 0 → ∞ → ( ) lim y t y Δ k t → ∞ = = p u u Δ Δ “hidden slide” provided for additional background

Process Gain: Controller Implications Long-term behavior from steady-state information y k u Δ = p Δ Controller really serves as “inverse” of process 1 u y Δ = Δ k p Process with higher gain is generally easier to control, all else being equal …

Process Zero: Controller Implications For “tight” control, controller is “inverse” of process ( ) ( ) y ( s ) g s u s = p 1 ( ) ( ) u ( s ) y s = g s p Inverse of g p (s) is unstable if g p (s) has a right-half-plane zero

Transfer Function to State Space l An infinite number of state space models can yield a given transfer function model l Two different state space “realizations” are normally used Ø Controllable canonical form Ø Observable canonical form

Multivariable Systems: Properties 2 input – 2 output Example ( ) ( ) ( ) ( ) ( ) y s g s u s g s u s = + 1 11 1 12 2 ( ) ( ) ( ) ( ) ( ) y s g s u s g s u s = + 2 21 1 22 2 Matrix-vector Notation ( ) ( ) ( ) ( ) y s g s g s u s ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ 1 11 12 1 = ⎢ ⎥ ⎢ ⎥ ⎢ ( ) ⎥ ( ) ( ) ( ) y s g s g s u s ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ 2 21 22 2 ( ) ( ) y s G u s = p Similar to SISO, controller serves as “process inverse” 1 ( ) ( ) u s G y s − = p

Steady-state Implications System 1 System 2 y 1 u 1 u y 1 u 1 u = + = + 1 1 2 1 1 2 y 1 u 1 u y 1 u 1 u = − + = + 2 1 2 2 1 2 1 u K y − = p 1 1 1 1 ⎡ ⎤ ⎡ ⎤ K K = = ⎢ ⎥ ⎢ ⎥ p 1 p 2 1 1 1 1 − ⎣ ⎦ ⎣ ⎦ Which system will be more difficult to control?

Steady-state Implications System 1 System 2 y 1 u 1 u y 1 u 1 u = + = + 1 1 2 1 1 2 y 1 u 1 u y 1 u 1 u = − + = + 2 1 2 2 1 2 1 u K y − = p 1 1 1 1 ⎡ ⎤ ⎡ ⎤ K K = = ⎢ ⎥ ⎢ ⎥ p 1 p 2 1 1 1 1 − ⎣ ⎦ ⎣ ⎦ Inverse does not exist (Gain Matrix is singular, rank = 1) Cannot independently control both outputs

Another Example y 1 u 0 . 95 u = + 1 1 2 y 1 u 1 u = + 2 1 2 1 0 . 95 ⎡ ⎤ K = ⎢ ⎥ p 3 1 1 ⎣ ⎦ Inverse exists – No longer singular (“gut feeling” that there may still be a problem … )

Directional Sensitivity: Singular Value Decomposition (SVD) T G = U Σ V Left singular Right singular Matrix of Gain matrix singular Vector matrix Vector matrix values Example Max singular value T 1 0 . 95 0 . 70 0 . 72 1 . 98 0 0 . 72 0 . 70 − − − − ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 1 1 0 . 72 0 . 70 0 0 . 0253 0 . 70 0 . 72 − − ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦                           G U V Σ Strongest Strongest input direction output direction

Example, cont’d T 1 0 . 95 0 . 70 0 . 72 1 . 98 0 0 . 72 0 . 70 − − − − ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 1 1 0 . 72 0 . 70 0 0 . 0253 0 . 70 0 . 72 − − ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦                           G U V Σ Input in Strongest Direction Output in Strongest Direction u y 0 . 72 1 0 . 95 0 . 72 1 . 38 ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ 1 1 = = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ u 0 . 70 y 1 1 0 . 70 1 . 41 ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ 2 2 y = K p u Input in Weakest Direction Output in Weakest Direction u y 0 . 70 1 0 . 95 0 . 70 0 . 018 ⎡− ⎡− ⎡− ⎡ ⎤ ⎤ ⎡ ⎤ ⎡ ⎤ ⎤ ⎤ 1 1 = = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ u 0 . 72 y 1 1 0 . 72 0 . 018 ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ 2 2 Same magnitude input = very different magnitude output

Example, cont’d T 1 0 . 95 0 . 70 0 . 72 1 . 98 0 0 . 72 0 . 70 − − − − ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 1 1 0 . 72 0 . 70 0 0 . 0253 0 . 70 0 . 72 − − ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦                           G U V Σ High condition number indicates problems Condition number = σ 1 / σ 2 = 1.98/0.0253 = 78 Note : For SVD analysis it is important to properly scale the inputs and outputs

MV Dynamic Properties: Quadruple Tank Operating Point 1 – Minimum Phase “Transmission zeros” are both ⎡ ⎤ negative 2.6 1.5 ⎢ ⎥ ( ) 62 s + 1 ( ) -1 62 s + 1 23 s + 1 ⎢ ⎥ z = − 0.060 and − 0.018 sec ( ) = G 1 s ⎢ ⎥ 1.4 2.8 ⎢ ⎥ 30 s + 1 ) 90 s + 1 90 s + 1 ( ( ) ( ) ⎣ ⎦ Operating Point 2 – Nonminimum Phase (RHPT zero) ⎡ ⎤ 1.5 2.5 ⎢ ⎥ 63 s + 1 ( 39 s + 1 ) 63 s + 1 ( ) -1 ⎢ ⎥ z = − 0.057 and + 0.013 sec ( ) = G 2 s ⎢ 2.5 1.6 ⎥ ⎢ ⎥ ( 56 s + 1 ) 91 s + 1 ( ) ( 91 s + 1 ) ⎣ ⎦ Matrix inverse is unstable

Multivariable Systems l Can have right-half-plane “transmission zeros” even when no individual transfer function has a RHP zero l Can have individual RHP zeros yet not have a RHPT zero Ø Fine performance when constraints are not active Ø May fail when one constraint becomes active or a loop is “opened” l Can exhibit “directional sensitivity” – with some setpoint directions much easier to achieve than others l Some of these MV properties cause challenges independent of control strategy selected

Summary l Nonlinear Model Ø Solve for steady-state, then linearize (Taylor series expansion) l State Space (linear) Model Ø Deviation variable form (perturbations from steady-state) l Dynamics Ø Eigenvalues of A matrix = poles of transfer function matrix n Pole-zero cancellation may reduce number of poles Ø Right-half-plane zeros = inverse response Ø MV: Transmission Zeros l Long term behavior from steady-state gain l Singular Value Decomposition (SVD)

Suggestion l Work through Workshop on the Three Tank Model Derive state space model • MATLAB Commands • Eigenvalues (eig) • State Space model (ss), step input (step) • Integration using ode45 • Single integration • In a “loop,” integrating from time step to time step • Convert state space to transfer function (ss2tf) •

Discrete Linear Models B. Wayne Bequette • Sampling Rules/Assumptions • Continuous to Discrete • Z-transform • Dynamic Properties of Discrete Systems Chemical and Biological Engineering

Discrete Models Input held constant between sample times (zero-order hold) Sample time is constant Rule of thumb for discrete control – select sample time roughly 1/10 of dominant time constant

Discrete Models x x u = Φ + Γ k + 1 k k State Space y Cx Du = + k k k Some texts/papers have different sign conventions y a y a y  a y Input-Output = − − − − + k 1 k 1 2 k 2 n k n − − − b u b u b u  b u (ARX) + + + + 0 k 1 k 1 2 k 2 m k m − − − usually b 0 = 0 ( ) [ ] y z Z y Z-transform = k 1 [ ] ( ) Z y z y z − Backwards shift operator − = k 1 ( ) ( ) 1 n 1 n 1 a z ... a z a z y z − − + − + + + + = 1 n 1 n Solve for y(z) − ( ) ( ) 1 m 1 m b b z − ... b z − + b z − u z + + + + 0 1 m 1 m −

Discrete Models 1 m 1 m b b z ... b z b z − − + − + + + + ( ) 0 1 m 1 m ( ) y z u z − = 1 n 1 n 1 a z ... a z a z − − + − + + + + 1 n 1 n − − 1 + ... + b − m + 1 + b m z − m Discrete ( ) = b 0 + b 1 z m − 1 z g p z 1 z − 1 + ... + a n − 1 z − n + 1 + a n z − n Transfer Function 1 + a n z Multiply by n z n + b n − 1 + ... + b m z n − m b 0 z 1 z ( ) = g p z z n + a 1 z − 1 + ... + a n − 1 z − n + 1 + a n ( )( ) ( ) z z z z  z z − − − poles = roots of denominator polynomial ( ) zeros = roots of numerator polynomial g z K 1 2 m = ⋅ p pz ( )( ) ( ) z p z p  z p − − − 1 2 n

Stability Stable Stable Continuous Discrete poles in LHP are stable poles inside unit circle are stable

Example − 1 b 1 z b 1 ( ) = ( ) = − a 1 y k − 1 ( ) + b ( ) g p z 1 + a 1 z − 1 = y k 1 u k − 1 z + a 1 p = − a 1 Value of the pole a 1 = 0.5, − 0.5,1.5, − 1.5 p = − 0.5,0.5, − 1.5,1.5 Let y(0) = 1, u(k) = 0 a 1 0.5 -0.5 1.5 -1.5 p -0.5 0.5 -1.5 1.5 y (1) -0.5 0.5 -1.5 1.5 y (2) 0.25 0.25 2.25 2.25 y (3) -0.125 0.125 -3.375 3.375 y (4) 0.0625 0.0625 5.0625 5.0625 Characteristic O s c i l l a t o r y, M o n o t o n i c , O s c i l l a t o r y, M o n o t o n i c , behavior stable stable unstable unstable

Simple Example: Results l Poles inside circle Ø stable l Poles outside circle Ø unstable l Negative poles Ø oscillate l First-order discrete systems can oscillate Ø This cannot happen with continuous systems

Discrete zeros l Zeros outside unit circle Ø Inverse is unstable (not necessarily inverse response) l Any continuous system with relative order >2 will have an unstable inverse with a small enough sample time Astrom, KJ, P Hagander & J Sternby “Zeros of Sampled Systems,” Proceedings of the 1984 American Control Conference , 1077-1081 Relative order = 3 Relative order = 1 Sample time = 1 1 ( )( ) 0 . 0803 z 1 . 7990 z 0 . 1238 + + ( ) = ( ) g p s g p z = 3 ( ) 3 s + 1 ( ) z 0 . 3679 − poles = -1 (multiplicity 3) poles = 0.3679 (multiplicity 3) zeros = -1.7990 & -0.1238 Unstable inverse

Final & Initial Value Theorems ( ) y z − 1 ( ) = lim ( ) n →∞ y n Δ t lim z → 1 1 − z Final value theorem unit step 1 ( ) = g p z ( ) u z ( ) = g p z ( ) ⋅ y z Long-term step response 1 − z − 1 1 ( ) ( ) − 1 − 1 ( ) = lim ( ) = lim ( ) ( ) lim t → ∞ y t z → 1 1 − z y z z → 1 1 − z g p z − 1 = lim z → 1 g p z 1 − z So, simply set z = 1 in g p (z) for long-term unit step response ( ) y z − 1 ( ) = lim ( ) Initial value theorem lim t → 0 y t z →∞ 1 − z “hidden slide” provided for additional background

Final Value Theorems Long-term outputs for unit step inputs 1 ( ) continuous g p s = 3 ( ) s 1 + 1 ( ) g p s 0 1 → = 3 = ( ) 0 1 + ( )( ) 0 . 0803 z 1 . 7990 z 0 . 1238 + + ( ) g p z discrete = 3 ( ) z 0 . 3679 − ( )( ) 0 . 0803 1 1 . 7990 1 0 . 1238 0 . 2525 + + ( ) g p z 1 1 → = = = 3 ( ) 0 . 2525 1 0 . 3679 − “hidden slide” provided for additional background

State Space Models ? x x u  x Ax Bu = Φ + Γ = + k + 1 k k y Cx Du y Cx Du = + = + k k k Continuous Discrete Finite differences approximation for derivative ( ( ) ) ( ) x k 1 t x k t x x + Δ − Δ −  k 1 k x + ≈ = t t Δ Δ x x − k 1 k Ax Bu + ≈ + k k t Δ [ ] x I A t x B tu = + Δ + Δ k 1 k k + Φ Γ How good are the approximations?

Exact Discretization t t + Δ A t A t k A ( ) ( ) ( ) x t t e x t e e d Bu t Δ Δ − σ + Δ = + ∫ σ k k k t k ( ) A t A t 1 x e x e I A Bu Δ Δ − = + − k 1 k k + A t e Δ Φ = Exact ( ) A t 1 e I A B Δ − Γ = − [ ] I A t Φ = + Δ Approximate B t Γ = Δ

Example Discretization  x 0 . 1 0 x 0 . 1 − ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ 1 1 u = + 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥  x 0 . 04 0 . 04 x 0 − ( ) g p s = ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ 2 2 ( )( ) 10 s 1 25 s 1 + + x ⎡ ⎤ 1 [ ] y 0 1 = ⎢ ⎥ x ⎣ ⎦ 2 Δ t 3 0 . 3 0 0 . 7408 0 = − ⎡ ⎤ ⎡ ⎤ A t e exp Δ Φ = = = ⎢ ⎥ ⎢ ⎥ 0 . 12 0 . 12 0 . 0974 0 . 8869 − ⎣ ⎦ ⎣ ⎦ Exact 1 − 0 . 1 0 0 . 1 0 . 2592 − ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ( ) I Γ = Φ − = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 0 . 04 0 . 04 0 0 . 0157 − ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ 1 0 0 . 3 0 0 . 7 0 − ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ Φ = + = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 0 1 0 . 12 0 . 12 0 . 12 0 . 88 − Approximate ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ 0 . 3 ⎡ ⎤ Γ = ⎢ ⎥ 0 ⎣ ⎦

Example, Continued Discrete Transfer Function − 1 ( ) ( ) g p z C zI D = − Φ Γ + 1 2 0 . 0157 z 0 . 0136 0 . 0157 z 0 . 0136 z − − + + Exact ( ) g p z = = 2 1 2 z 1 . 6277 z 0 . 65702 1 1 . 6277 z 0 . 65702 z − − − + − + Poles/eigenvalues = 0.7408 & 0.8869 2 0 . 036 0 . 036 z − ( ) g p z = = 2 1 2 z 1 . 58 z 0 . 616 1 1 . 58 z 0 . 616 z − − − + − + Approximate Poles/eigenvalues = 0.7 & 0.88 “hidden slide” provided for additional background

Example Step Response Model 2 1.5 s s 5 Etc. 4 1 s 3 temp, deg C s 2 0.5 s 1 0 -2 0 2 4 6 8 discrete-time step 1 0.8 0.6 psig 0.4 0.2 0 -2 0 2 4 6 8 discrete-time step T ! #  S = s 1 s 2 s 3 s 4 s 5 s N " $

Step Response Model ∞ y s u = ∑ Δ k i k i − i 1 = s u  s u s u  s u = Δ + + Δ + Δ + + Δ 1 k 1 N k N N 1 k N 1 N k − − + − − + ∞ − ∞   y s u s u s u s u = Δ + + Δ + Δ + + Δ k 1 k 1 N 1 k N 1 N k N N k − − − + − − ∞  (  ) s u s u s u u , = Δ + + Δ + Δ + + Δ          1 k 1 N 1 k N 1 N k N k − − − + − − ∞ u k N − N 1 − y s u s u . ∑ = + Δ k N k N i k i − − i 1 =

Example Impulse Response Model 1 0.8 0.6 h 0.4 1 temp, deg C h 2 h 0.2 3 0 -0.2 -2 0 2 4 6 8 10 discrete-time step 1 0.8 0.6 psig 0.4 0.2 0 -2 0 2 4 6 8 10 discrete-time step h i = s i − s i − 1 Impulse and step response i coefficients are related ∑ s i = h j j = 1

Parameter Estimation, ARX Models model output measured output known (measured) input ˆ y a y a y b u b u = − − + + 1 1 0 2 1 1 0 2 1 − − ˆ y a y a y b u b u = − − + + 2 1 1 2 0 1 1 2 0 time step index  ˆ y a y a y b u b u = − − + + N 1 N 1 2 N 2 1 N 1 2 N 2 − − − − ˆ y y y u u a − ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ could change 1 0 1 0 1 1 − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ sign on columns . . . . . a N − 2 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 1 & 2 of Φ = pts . . . . . b ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ˆ y y y u u b ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ N N 1 N 2 N 1 N 2 2 − − − − ˆ Y = Φ Θ

Optimization/Parameter Solution N 2 ( ˆ ) Objective Function min y y ∑ − i i a , a , b , b 1 2 1 2 i 1 = T Y − ˆ N T Y −ΦΘ 2 = Y − ˆ ( ) ( ) = Y −ΦΘ ∑ ( y i − ˆ ) ( ) ( ) y i Y Y i = 1 − 1 Φ ( ) T Φ T Y Solution Θ = Φ measured parameter output estimate vector vector

Example PRBS estimation example 0.6 0.4 0.2 y 0 -0.2 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 time, min 1 0.5 PRBS 0 u input -0.5 -1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 time, min

Result 0 . 0889 0 1 1 0 . 0137 − ⎡ ⎤ ⎡ ⎤ − 1 Φ ( ) T Φ T Y ⎢ ⎥ ⎢ ⎥ 0 . 0137 0 . 0889 1 1 0 . 1564 Θ = Φ − ⎢ ⎥ ⎢ ⎥ 0 . 1564 0 . 0137 1 1 0 . 4618 ⎢ − ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 0 . 4618 0 . 1564 1 1 0 . 1771 − ⎢ ⎥ ⎢ ⎥ ⎡ ⎤ ⎡ ⎤ 1.1196 − a 1 ⎢ 0 . 1771 0 . 4618 1 1 ⎥ ⎢ 0 . 3446 ⎥ − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 0 . 3446 0 . 1771 1 1 0 . 2171 − − ⎢ ⎥ ⎢ ⎥ − 0.3133 − a 2 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 0 . 2171 0 . 3446 1 1 0 . 1558 − − ⎢ ⎥ ⎢ ⎥ Θ = = ⎢ ⎥ ⎢ ⎥ 0 . 1558 0 . 2171 1 1 0 . 0485 − − ⎢ ⎥ ⎢ ⎥ − 0.0889 b 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 0 . 0485 0 . 1558 1 1 0 . 1879 − − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 0.2021 b 2 0 . 1879 0 . 0485 1 1 Y 0 . 1123 ⎣ ⎦ ⎣ ⎦ Φ = ⎢ − ⎥ = ⎢ − ⎥ ⎢ ⎥ ⎢ ⎥ 0 . 1123 0 . 1879 1 1 0 . 0463 − − ⎢ ⎥ ⎢ ⎥ 0 . 0463 0 . 1123 1 1 0 . 2003 ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ 0 . 2003 0 . 0463 1 1 0 . 5007 − ⎢ ⎥ ⎢ ⎥ b z b 0 . 0889 z 0 . 2021 + − + ⎢ ⎥ ⎢ ⎥ 0 . 5007 0 . 2003 1 1 0 . 3846 g p ( ) z 1 2 − − = = ⎢ ⎥ ⎢ ⎥ 2 2 z a z a z 1 . 1196 z 0 . 3133 + + − + 0 . 3846 0 . 5007 1 1 0 . 0172 − − ⎢ ⎥ ⎢ ⎥ 1 2 ⎢ ⎥ ⎢ ⎥ 0 . 0172 0 . 3846 1 1 0 . 1513 1 2 0 . 0889 z 0 . 2021 z − − − − − + ⎢ ⎥ ⎢ ⎥ = 0 . 1513 0 . 0172 1 1 0 . 1162 − − − ⎢ ⎥ ⎢ ⎥ 1 2 1 1 . 1196 z − 0 . 3133 z − − + ⎢ ⎥ ⎢ ⎥ 0 . 1162 0 . 1513 1 1 0 . 1134 − − ( ) 0 . 0889 z 2 . 274 ⎢ ⎥ ⎢ ⎥ − − 0 . 1134 0 . 1162 1 1 0 . 0502 ⎢ ⎥ ⎢ ⎥ = − − − ⎣ ⎦ ⎣ ⎦ ( z 0 . 5716 )( z 0 . 5481 ) − −

Identified Model b z b 0 . 0889 z 0 . 2021 + − + ( ) g p z 1 2 = = 2 2 z a z a z 1 . 1196 z 0 . 3133 + + − + 1 2 1 2 0 . 0889 z 0 . 2021 z − − − + ( ) ( ) y z u z = ⋅ 1 2 1 1 . 1196 z 0 . 3133 z − − − + y 1 . 1196 y 0 . 3133 y 0 . 0889 u 0 . 2021 u = − − + k k 1 k 2 k 1 k 2 − − − −

Subspace Identification l Subspace ID can be used to develop discrete state space models from input-output data

Summary l Discrete models Ø State space, ARX, discrete transfer function l Zeros & poles Ø Poles inside unit circle are stable (negative = oscillate) Ø Zeros inside unit circle have stable inverses l Parameter estimation Ø Example with PRBS input l Step and impulse response models

Discrete (Digital) Control B. Wayne Bequette • Review of Digital PID • Review of Model-based Digital Control • Discrete Internal Model Control • Examples • Summary Chemical and Biological Engineering

Proportional-Integral-Derivative (PID) Control Error = setpoint – measured output ( ) ( ) ( ) e t r t y t = − ( ) 1 de t ⎡ ⎤ t ( ) ( ) ( ) u t u k e t e t dt = + + ∫ + τ ⎢ ⎥ 0 c D dt 0 τ ⎣ ⎦ I Manipulated Proportional Integral time Derivative time Input gain ( ) ( ) = e t k e k t k k k ( ) ( ) ( )  ( ) ( ) ( ) e t dt e t t e t t e t t e t t e i t ∑ ∑ ≈ ⋅ Δ + ⋅ Δ + + ⋅ Δ ≈ Δ = Δ ∫ 1 2 k i i 1 i 1 0 = = ( ) de t k ( ) − e k − 1 ( ) ≈ e k dt Δ t

Substituting each term, we find the discrete controller ⎡ ⎤ k ( ) + Δ t e ( i ) + τ D ( ) ⎢ ∑ ⎥ ( ) = u ( ) − e k − 1 ( ) u k 0 + k c e k Δ t e k ⎢ ⎥ τ I ⎣ ⎦ i = 0 It is convenient to work with the “ velocity form. ” First, generate the equation for step k-1 k 1 t − ⎡ ⎤ Δ τ ( ) ( ) ( ( ) ( ) ) ⎥ u k 1 u k e k 1 e ( i ) D e k 1 e k 2 ∑ − = + − + + − − − ⎢ 0 c t τ Δ ⎣ ⎦ i 0 I = Then subtract u(k-1) from u(k) to find ⎡ ⎤ ⎛ ⎞ ⎛ ⎞ 1 + Δ t ( ) + − 1 − 2 τ D + τ D ) + τ D ⎢ ⎥ ⎜ ⎟ ( ) = u k − 1 ( ) + k c ( ( ) u k ⎟ e k ⎟ e k − 1 Δ t e k − 2 ⎜ ⎜ ⎢ ⎝ ⎠ ⎥ Δ t Δ t τ I ⎝ ⎠ ⎣ ⎦

⎡ ⎤ ⎛ ⎞ ⎛ ⎞ 1 + Δ t ( ) + − 1 − 2 τ D + τ D ) + τ D ⎢ ⎥ ⎜ ⎟ ( ) = u k − 1 ( ) + k c ( ( ) u k ⎟ e k ⎟ e k − 1 Δ t e k − 2 ⎜ ⎜ ⎢ ⎝ ⎠ ⎥ Δ t Δ t τ I ⎝ ⎠ ⎣ ⎦ " % b 2 = k c 1 + Δ t + τ D ' , $ • Δ t τ I # & " % 1 = − k c 1 + 2 τ D b ' , $ # Δ t & b 0 = k c τ D Δ t ( ) = u k − 1 ( ) + b 2 e k ( ) + b ( ) + b 0 e k − 2 ( ) u k 1 e k − 1

( ) ( ) ( ) ( ) ( ) u k u k 1 b e k b e k 1 b e k 2 − − = + − + − 2 1 0 [ ] [ ] ( ) ( ) ( ) ( ) u k u z , e k e z , Ζ = Ζ = [ ] 1 [ ] 1 ( ) ( ) ( ) ( ) u k 1 z u z , e k 1 z e z , − − • Ζ − = Ζ − = [ ( ) ] 2 ( ) . [ ( ) ] 2 ( ) . u k 2 z u z e k 2 z e z − − Ζ − = Ζ − = ( ) ( ) ( ) ( ) 1 1 2 1 z u z b b z b z e z − − − − = + + 2 1 0 ( ) ( ) ( ) 1 2 b b z b z − − + + 2 1 0 u z e z = 1 1 z − − ( ) ( ) ( ) 2 b z b z b + + 2 1 0 ( ) ( ) u z e z g z e z = = c 2 z z −

PID Tuning l Usually based on continuous-time methods, including: l Ziegler-Nichols closed-loop oscillations l Cohen-Coon l IMC-based PID l Skogestad’s tuning method l Frequency response l SWAG (most common)

IMC-Based PID g p ( s ) = k p e − θ s l Continuous Model  τ p s + 1 l PID Parameters ( ) 0 . 5 τ + θ p k , = c ( ) k 0 . 5 λ + θ p 0 . 5 , τ = τ + θ I p τ θ p . τ = D 2 τ + θ p λ is the tuning parameter – related to the closed-loop time constant

Digital Control: Block Diagram • d(z) • r(z) • e(z) • u(z) • + • y(z) • + • g c (z) • g p (z) • + • _ l Stability analysis ( ) ( ) g z g z p c ( ) ( ) ( ) ( ) y z r z = ⋅ l Poles of CLTF must be 1 g z g z + p c inside the unit circle ( ) ( ) ( ) y z g z r z = CL Closed-loop transfer function

Suggestion l Work through Workshop on Digital PID Control Find discrete time model • Digital PID • Tune for Verge of Instability (continuous oscillations) • Find closed-loop poles • Increase gain for unstable closed-loop • Find closed-loop poles • Next Topic: Internal Model Control (following slides)

Manfred Morari

Digital Feedback Control d(z) r(z) e(z) u(z) + y(z) + g c (z) g p (z) + _ ( ) ( ) g z g z l Digital Controller Design p c ( ) ( ) ( ) ( ) y z r z = ⋅ Ø Response specification-based 1 g z g z + p c ( ) ( ) ( ) y z g z r z = CL Desired response

Controller Design ( ) ( ) g z g z p c ( ) g z Desired response = CL ( ) ( ) 1 g z g z solve for + p c controller ( ) 1 g z ( ) g z CL = ⋅ c ( ) ( ) g z 1 g z − p CL

Internal Model Form d ( z ) Process ~ u ( z ) + + r ( z ) r ( z ) q ( z ) g ( z ) y ( z ) p + - ~ IMC + - ~ y ( z ) controller g ( z ) p Process model ~ d ( z ) For controller design , consider perfect model and no disturbances

IMC Design r(z) u(z) y(z) q(z) g p (z) ( ) ( ) ( ) ( ) y z g z q z r z = p solve for ( ) ( ) ( ) y z g z r z = controller CL ( ) g z ( ) CL q z Desired response = ( ) g z p ( ) g z ~ 1 ( ) CL ( ) ( ) Really based on q z g z g − z = = ~ CL p ( ) g z the model p

Can implement IMC in Classic Feedback Form ( ) q z ( ) g z = ~ c ( ) ( ) 1 − q z g z p

Digital Controller Design l Deadbeat John Ragazzini l Dahlin’s Controller l State Deadbeat l State Deadbeat with Filter (Vogel-Edgar) l Modified Dahlin’s Controller l NEW (IMC) DESIGN

Deadbeat l Achieve setpoint in the minimum time (if no time-delay, then one time step) y r = k + 1 k Backwards shift operator 1 ( ) ( ) y z z r z − = 1 ( ) g CL z = z − 1 ~ − IMC Form 1 ( ) ( ) q z z g z − = p 1 z − ~ 1 g ( ) z g ( ) z − Classic Feedback Form = ⋅ c p 1 1 z − −

Example 0, FO Process 1 1 0 . 0952 z − ~ ~ ( ) g p s ( ) = g p z = 10 s 1 Δ t = 1 1 + 1 0 . 9048 z − − ( ) 1 1 z 1 0 . 9048 z z 0 . 9048 − − − − ( ) q z = = 1 0 . 0952 z 0 . 0952 z 0 − + ( ) ( ) ( ) u z q z r z = Control action 1 ( ) let r z unit step setpoint change = = 1 1 z − − 1 2 1 z 0 . 9048 z 10 . 504 9 . 504 z − − − − − ( ) u z = = ( ) 1 1 1 0 . 0952 z 1 z 1 z − − − − − 1 10 . 504 9 . 504 z − Large control action up, − ( ) u z = + 1 1 1 z − 1 z − then down − −

Classic Feedback Form for Deadbeat Design 1 z 1 − ~ ~ 1 1 ( ) ( ) ( ) g z g z g z − − = ⋅ = ⋅ c p p 1 1 z z 1 − − − 1 0 . 0952 z 0 . 0952 − ~ ( ) g z = = p 1 1 0 . 9048 z z 0 . 9048 − − − 1 1 0 . 9048 z z 0 . 9048 − − − ~ 1 ( ) g z − = = p 1 0 . 0952 z 0 . 0952 − z 0 . 9048 1 1 z 0 . 9048 − − ⎛ ⎞ ( ) g z = ⋅ = ⋅ ⎜ ⎟ c 0 . 0952 z 1 0 . 0952 z 1 − − ⎝ ⎠ 1 1 1 1 0 . 9048 z 1 0 . 9048 z − − − − ⎛ ⎞ 10 . 5 = ⋅ = ⋅ ⎜ ⎟ 1 1 0 . 0952 1 z 1 z − − − − ⎝ ⎠

Classic Feedback Implementation 1 1 0 . 9048 z − − ( ) g z 10 . 5 = ⋅ c 1 1 z − − ( ) ( ) ( ) u z g z e z = c ( ) ( ) ( ) ( ) 1 1 1 z u z 10 . 5 1 0 . 9048 z e z − − − = ⋅ − ( ) u u 10 . 5 e 0 . 9048 e − = ⋅ − k k 1 k k 1 − − ( ) u u 10 . 5 e 0 . 9048 e = + ⋅ − k k 1 k k 1 − −

Note that this is a PI controller ( ) u k = u k − 1 + 10.5 ⋅ e k − 0.9048 e k − 1 ( ) = u k − 1 ( ) + b 2 e k ( ) + b ( ) + b 0 e k − 2 ( ) u k 1 e k − 1 " % b 2 = k c 1 + Δ t + τ D ' = 10.5 $ Δ t τ I # & " % 1 = − k c 1 + 2 τ D b ' = − 10.5*0.9048 = − 9.5 $ # Δ t & b 0 = k c τ D Δ t = 0 τ I = 9.5 The PI parameters are then k c = 9.5

IMC Implementation d ( z ) Process ~ u ( z ) + r ( z ) + r ( z ) q ( z ) g ( z ) y ( z ) p + - ~ IMC + - ~ y ( z ) controller g ( z ) p Process model ~ d ( z ) ~ ( ) ( ) 1 1 ( ) 1 0 . 9048 z y z 0 . 0952 z u z − − − = ~ ~ ~ ~ ( ) ( ) ( ) y z g z u z = y 0 . 9048 y 0 . 0952 u = + p k k 1 k 1 − − ~ ~ ~ ~ ( ) ( ) ( ) d z y z y z d y y = − = − k k k ~ ~ ~ ~ ( ) ( ) ( ) r z r z d z r r d = − = − k k k ~ ~ ( ) ( ) ( ) ( ) ( ) u z q z r z 1 1 1 ( ) = 0 . 0952 z u z z 1 0 . 9048 z r z − − − = − 1 ~ ⎛ ⎞ ( ) ( ) 1 ( ) u z 1 0 . 9048 z − r z = − ⎜ ⎟ 0 . 0952 ⎝ ⎠ 1 ⎛ ⎞ ~ ~ ( ) u r 0 . 9048 r = − ⎜ ⎟ k k k 1 − 0 . 0952 ⎝ ⎠

Discussion Standard feedback control clearly has the current control action as a function of the previous control action. Why doesn’t IMC? e r y Standard Feedback = − k k k ( ) u u 10 . 5 e 0 . 9048 e = + ⋅ − k k 1 k k 1 − − ~ ~ y 0 . 9048 y 0 . 0952 u IMC = + k k 1 k 1 − − ~ ~ d y y = − k k k ~ ~ r r d = − k k k 1 ⎛ ⎞ ~ ~ ( ) u r 0 . 9048 r = − ⎜ ⎟ k k k 1 − 0 . 0952 ⎝ ⎠

Example 0 (First-order), Deadbeat Design deadbeat, FO example 1 0.5 y 0 -2 0 2 4 6 8 time 10 5 u 0 -2 0 2 4 6 8 time

Example 1, Second-Order Process 1 2 1 0 . 0157 z 0 . 0136 z − − + ~ ~ g p ( ) s ( ) g p z = = 1 2 ( )( ) 10 s 1 25 s 1 1 1 . 6277 z 0 . 65702 z Δ t = 3 − − + + − + zero = -0.0136/0.0157 = -0.8694 1 2 3 2 z 1 . 6277 z 0 . 65702 z z 1 . 6277 z 0 . 65702 − − − − + − + q ( ) z = = 1 2 2 0 . 0157 z 0 . 0136 z 0 . 0157 z 0 . 0136 z 0 − − + + + controller has pole = -0.8694, which is stable, but oscillates “ringing” behavior as shown next

Example 1 (Second-order): Deadbeat Design deadbeat, SO example 1.5 1 y 0.5 “ringing” behavior (intersample) Perfect control, at sample times 0 -5 0 5 10 15 20 25 30 time 100 50 0 u -50 -100 -5 0 5 10 15 20 25 30 time

Dahlin’s Controller l Desired first-order + time-delay response to setpoint change 1 ( ) 1 z − − α ( ) g CL z For no time-delay: where α = exp(- Δ t/ λ ) = 1 1 z − − α 1 ( ) ( ) ( ) g z 1 z 1 − − α − α ~ ~ ( ) 1 ( ) 1 ( ) q z CL g z g z − − = = ⋅ = ⋅ ~ p p ( ) 1 ( ) g z 1 z z − − α − α p For the second-order example: ( ) 1 1 2 1 z 1 1 . 6277 z 0 . 65702 z − − − − α − + ( ) q z = ⋅ 1 1 2 1 z − 0 . 0157 z − 0 . 0136 z − − α + 2 ( ) 1 z 1 . 6277 z 0 . 65702 − α − + = ⋅ ( ) z 0 . 0157 z 0 . 0136 − α +

Second-order Process: Dahlin’s Controller Dahlin, SO example 1.5 For λ = 3, α = 0.3679 1 y 0.5 0 First-order response, at sample times -5 0 5 10 15 20 25 30 time 100 50 0 u -50 -100 -5 0 5 10 15 20 25 30 time Still “ringing”, but more damped than deadbeat

Dahlin’s Modified Controller l Substitute zeros at origin for unstable (|zero|>1) or ringing (|zero|<1 but negative) zeros. Also, keep the gain the same l Problem: Dahlin applied this to g c (z), but should have applied it to the IMC form, q(z)! “hidden slide” provided for additional background

State Deadbeat Controller Design l Brings outputs & inputs to new steady-state in the minimum number of time steps l Does not invert zeros at all 1 2 0 . 0157 z 0 . 0136 z − − + ~ second-order example ( ) g p z = 1 2 1 1 . 6277 z 0 . 65702 z − − − + 1 2 0 . 0157 0 . 0136 0 . 0157 z 0 . 0136 z − − + + = ⋅ 1 2 1 1 . 6277 z 0 . 65702 z 0 . 0157 0 . 0136 − − − + + ~ ~ ( ) ( ) g p − z g p + z 1 2 2 1 1 . 6277 z 0 . 65702 z z 1 . 6277 z 0 . 65702 − − − + − + ~ 1 ( ) q z g − = = = p 2 − 0 . 0293 0 . 0293 z 0 z 0 + +

Second-order Example: State Deadbeat Design State Deadbeat, SO example 1.5 1 y 0.5 0 -5 0 5 10 15 20 25 30 time 100 50 0 u -50 -100 -5 0 5 10 15 20 25 30 time

Controller Forms a + = zero outside unit circle a − zero inside unit circle = i i ( ) ( )( ) ( ) K z a  z a z a  z a − − + + − − − − ~ pz 1 k k 1 n 1 ( ) ( ) g z G z + − = = p ( ) ( ) z p  z p − − 1 n Zafiriou & Morari notation Assumes stable poles State Deadbeat: ( ) ( ) z p  z p − − ( ) 1 n q z = ( ) ( ) n SD K 1 a  1 a z − − − − pz 1 k State Deadbeat with Filter (Vogel-Edgar): 1 ( ) ( ) ( ) z p  z p 1 z − − − − α ( ) ( ) q z q F z 1 n = = ⋅ ( ) ( ) VE SD n 1 K 1 a  1 a z 1 z − − − − − − α pz 1 k Notice that neither the state deadbeat nor Vogel-Edgar controllers try to invert zeros (even good ones!)

IMC Design Summary l Model Factorization (“good stuff” and “bad stuff”) ~ ~ ~ ( ) ( ) ( ) g z g z g z = p p p − + • numerator one order • zeros outside unit circle less than denominator • zeros inside unit circle that are negative • “all-pass” by including pole at 1/z i for each positive z i outside the unit circle (but not the negative ones) l Controller: Invert “good stuff” part of model ~ − 1 ( ) ( ) ( ) q z g z F z = p − • “Filter” for desired response, often first-order

Example 2: Third-Order (3-tank) Sample time = 1 1 ( )( ) 0 . 0803 z 1 . 7990 z 0 . 1238 ~ + + ~ ( ) g p s g p ( ) z = = 3 3 ( ) s 1 ( ) z 0 . 3679 + − zeros at -1.7990 (outside unit circle) -0.1238 (inside, but negative) ~ ~ ~ ( ) ( ) ( ) g z g z g z = p p p − + 2 ( )( ) ( )( ) 0 . 0803 2 . 7990 1 . 1238 z z 1 . 7990 z 0 . 1238 + + ~ ( ) g p z = ⋅ 3 ) 2 ( )( ( ) 2 . 7990 1 . 1238 z z 0 . 3679 − Note gain is 1 (value at z = 1) 3 1 ( ) ( ) z 0 . 3679 1 z − − − α ~ 1 ( ) ( ) ( ) q z g z F z − = = ⋅ p 2 1 − 0 . 2526 z 1 z − − α 3 ( ) ( ) z 0 . 3679 1 − − α where α = exp(- Δ t/ λ ) = ⋅ 2 0 . 2526 z z − α

Example 2 (Third-Order) IMC Design For λ = 1, Discrete IMC, TO example α = 0.135 1.5 1 y 0.5 0 -2 0 2 4 6 8 10 time 5 4 3 2 u 1 0 -1 -2 0 2 4 6 8 10 time