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

Chemical and Biological Engineering

Model Predictive Control: Background

• B. Wayne Bequette
• Concise Review of Undergraduate Process

Control

• Introduction to Discrete-time Systems
• Digital PID
• Discrete Internal Model Control (IMC)

January 2015

Automation and Control

feedforward feedback

Control Algorithm

l Compares measured process output with desired setpoint,

and calculates a manipulated input (often a flowrate)

l Proportional-integral-derivative (PID)

Ø Long history, “workhorse”, lower-level control loops

l Model predictive control (MPC)

Ø Most widely applied “advanced” control algorithm Ø Constraints, multivariable systems

Common notation

r = setpoint (desired value of process output) y = measured process output e = error (setpoint – output, r – y) u = manipulated input (often a flowrate)

Proportional-Integral-Derivative (PID) Control

( ) ( ) ( ) ( ) ( ) ( ) ( )

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + + + = − =

∫

dt t de dt t e t e k u t u t y t r t e

D t I c

τ τ 1

Manipulated Input Error = setpoint - measured output Proportional gain Integral time Derivative time

l

Single-input, single-output design

Ø

Problems with one controller can impact another controller

l

Constraints

Ø

Can cause “windup” problems

l

Does not explicitly require a process model

Chemical and Biological Engineering

Continuous Linear Models State Space and Transfer Function

• B. Wayne Bequette
• Linearization
• State Space Form
• Transfer Function
• Step Responses
• MV Properties

Nonlinear ODE Models

l General “Lumped Parameter” Form

( )

p u x f x , , =  y = g x,u,p

( )

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ =

n

x x x x 

2 1

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ =

m

u u u u 

2 1

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ =

r

y y y y 

2 1

differential state equations algebraic output equations ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ =

q

p p p p 

2 1

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = dt dx dt dx dt dx x

n

 

2 1

states state derivatives (from accumulation term) inputs

• utputs

parameters

Linearization

s s s s s s s s

u x j i ij u x j i ij u x j i ij u x j i ij

u g D x g C u f B x f A

, , , ,

∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = = = = u D x C y u B x A x ʹ″ + ʹ″ = ʹ″ ʹ″ + ʹ″ = ʹ″ 

( )

p u x f x , , =  y = g x,u,p

( )

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − = ʹ″

ns n s s

x x x x x x x 

2 2 1 1

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − = ʹ″

ms m s s

u u u u u u u 

2 2 1 1

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − = ʹ″

rs r s s

y y y y y y y 

2 2 1 1

‘perturbation’ or ‘deviation’ variables

eqn i state j eqn i input j

steady-state solution, us, xs, ys Where:

Example: Van de vuuse Reaction

( )

B A B B A A A Af A

C k C k C V F dt dC C k C k C C V F dt dC

2 1 2 3 1

− + − = − − − =

[ ]

Bs B Afs Af s Bs B As A

C C x y C C V F V F u C C C C x x x − = = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =

2 2 1

, ,

[ ] [ ]

, 1 , , 2

2 1 3 1

= = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − − = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − − − − − = D C C V F C C B k V F k C k k V F A

Bs s As Afs s As s

= f1 = f2

etc C A u x

As s s

C f x f A

, 1 , 1 1 11

∂ ∂ ∂ ∂ = =

“hidden slide” provided for additional background

Continuous State Space Model

Du Cx y Bu Ax x + = + =  ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡

m rm r m n rn r n r m nm n m n nn n n n

u u d d d d x x c c c c y y u u b b b b x x a a a a x x                            

1 1 1 11 1 1 1 11 1 1 1 1 11 1 1 1 11 1

n by n n by m r by n r by m Assumes deviation variable form

Stability

l Eigenvalues of A

Ø n x n A matrix = n eigenvalues = n states Ø If all eigenvalues have a negative real portion – stable Ø If any eigenvalue has a positive real portion – unstable Ø Complex – generally ‘oscillatory’

( )

det = − A I λ

l Characteristic Polynomial (n roots):

Ø Must have at least 2 states to oscillate (be complex)

eigenvalues are the roots of the equation Identity matrix

Example: CSTR at 2 Operating Points

Operating condition 1 Operating condition 2

A

1 =

−1.1680 −0.0886 2.0030 −0.2443 ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ A2 = −1.8124 −0.2324 9.6837 1.4697 ⎡ ⎣ ⎢ ⎤ ⎦ ⎥

( ) ( )( ) ( )( )

4628 . 4123 . 1 003 . 2 0886 . 2443 . 1680 . 1 det 2443 . 0030 . 2 0886 . 1680 . 1

2 1 1

= + + = = − − + + = − ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + − + = − λ λ λ λ λ λ λ λ A I A I

λ =−0.8955 hr

• 1 and λ = −0.5168 hr
• 1

λ =−0.8366 hr

• 1 and λ = 0.4939 hr
• 1

Operating Point 2 = Unstable Operating Point 1 = Stable

Can show that the Eigenvalues for Operating Point 2 are:

“hidden slide” provided for additional background

Laplace Transform

l Convert Differential Equations to Algebraic Equations

Ø Design controllers using algebra rather than differential equations

l Easy Analysis of “Block Diagrams”

L f t

( )

[ ]= F s

( ) =

f t

( )e

−stdt ∞

∫

L 1

[ ]= 1

s

( ) ( ) ( )

f s sF dt t df L − = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ Definition of Laplace Transform Unit step Derivative

Laplace Transforms, cont’d

l Most Undergraduate Courses

Ø Much (perhaps too much) time is spent:

n Taking Laplace transform of process differential equation n Laplace transform of “forcing function” (typically a step) n Multiply, then perform partial fraction expansion n Invert each term back to the time domain for an analytical expression

l In Practice

Ø Step response behavior for process understanding Ø Main use of transfer function is for “controller synthesis” Ø Can easily convert from differential equations (“state space”)

to transfer function form using MATLAB, etc.

Ø Closed-loop block diagram analysis

Very painful, many nightmares?

Block Diagram Analysis

Setpoint Error Manipulated input Measured output

g c e(s) controller (s) u(s) g (s)

p

y(s) r(s) + − process

( ) ( ) ( ) ( ) ( ) ( )

s r s g s g s g s g s y

p c p c

⋅ + = 1

Closed-loop characteristic equation (roots determine stability) Routh Array, etc.

State Space to Transfer Function Form

Du Cx y Bu Ax x + = + = 

( ) ( ) ( ) ( ) ( ) ( )

s Du s Cx s y s Bu s Ax s sx + = + =

( ) ( ) ( ) ( ) ( )

[ ] ( )

( ) ( ) ( )

s u s G s y s u D B A sI C s y s Bu A sI s x

1 1

= + − = − =

− −

Take the Laplace Transform (s = transform variable)

Transfer function matrix

usually, D = 0

( ) ( )

[ ]

B A sI C s G

1 −

− =

Matrix Transfer Function Form

( ) ( ) ( )

s u s G s y =

y1 s

( )

 yr s

( )

! " # # # # $ % & & & & = g11 s

( )

 g1m s

( )

 gr1 s

( )

 grm s

( )

! " # # # # $ % & & & & u1 s

( )

 um s

( )

! " # # # # $ % & & & &

yi s

( ) = gij s ( )uj s ( )

r outputs m inputs

• utput i

input j r rows by m columns Transfer Function Matrix First subscript = output Second subscript = input

Dynamic Behavior: SISO

( )

1 1 1 1 1 1

a s a s a s a b s b s b s b s g

n m n n m m m m p

+ + + + + + + + =

− − − −

 

( ) ( )( ) ( ) ( )( ) ( )

n m pz p

p s p s p s z s z s z s k s g − − − − − − =  

2 1 2 1

( ) ( )( ) ( )

( )( ) ( )

1 1 1 1 1 1

2 1 2 1

+ + + + + + = s s s s s s k s g

pn p p nm n n p p

τ τ τ τ τ τ   Polynomial Pole-zero Gain-time constant Zeros = roots of numerator Poles = roots of denominator = determine stability

1 1

1 p

p

− = τ

Large time constant = small, negative, pole

Relative order = n-m

Example of Pole-Zero Cancellation

l Number of states = Number of poles (order of numerator

• f transfer function), except when some poles and zeros

‘cancel’

[ ] [ ]

1 8255 . 3 5301 . 1 5640 . 2 75 . 9056 . = = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡− = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − = D C B A

( ) ( ) ( )( )

2640 . 2 3 . 3 . 5302 . 1 6792 . 564 . 2 4590 . 5302 . 1

2

+ + + − = + + − − = s s s s s s s g p gp s

( ) =

−0.6758 0.4417s + 1

Bioreactor model

“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)

Real axis Imaginary axis x x more oscillatory faster response

• inverse response zeros

x unstable poles

x = pole

• = zero

“hidden slide” provided for additional background

First-Order ( ) ( )

s u s k s y

p p

1 + = τ τ p dy dt + y = kpu dy dt = − 1 τ p y + kp τ p u

• r

gain time constant

2 4 6 8 10 12 14 16 18 20 5 10 15 20 t, minutes y, deg C 2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 t, minutes u, kW

Example: gain = 1.43 oC/kW time constant = 5 minutes

• utput

input

First Order + Time-delay

( ) ( )

s u s e k s y

p s p

1 + =

−

τ

θ

( )

θ τ − = + t u k y dt dy

p p

gain time constant

time- delay time constant y u long-term change in process output change in process input 63.2 %

• f change

in process

• utput

time-delay

Second-order Underdamped

time period of oscillation a b c rise time decay ratio = b/a

• vershoot ratio = a/c

time to first peak

y s

( ) =

k τ

2s 2 + 2ζτs +1 ⋅u s

( )

Im Re 1

2 , 1 2 2 , 1

j p p ± = − ± − = τ ζ τ ζ

1 < ζ

step response

Numerator Dynamics

10 20 30 40 50 60

• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 1 1.2 time

• 15
• 10
• 3

3 10 15 20 various values of numerator time constant y

( ) ( ) ( )( )

1 15 1 3 1 + + + = s s s s g

n p

τ

n

z τ 1 − = > → < z

n

τ

“inverse response” (right-half-plane zeros) [challenging control problem]

response faster large → =

n

τ

Steady-state Gain

l For stable systems, the steady-state gain is found from the

long-term response

( ) ( )( ) ( )

( )( ) ( )

s u s s s s s s k s y

pn p p nm n n p

Δ ⋅ + + + + + + = 1 1 1 1 1 1

2 1 2 1

τ τ τ τ τ τ  

( ) ( ) ( )( ) ( )

( )( ) ( )

( ) ( )

u k s u k s s sy t y s u s s s s s s k s s sy t y

p p s t pn p p nm n n p s t

Δ = Δ ⋅ ⋅ = = Δ ⋅ + + + + + + ⋅ = =

→ ∞ → → ∞ → 2 1 2 1

lim lim 1 1 1 1 1 1 lim lim τ τ τ τ τ τ  

( )

u y u t y k

t p

Δ Δ = Δ =

∞ →

lim

Step input Final value theorem

“hidden slide” provided for additional background

Process Gain: Controller Implications

Long-term behavior from steady-state information Controller really serves as “inverse” of process u k y

pΔ

= Δ y k u

p

Δ = Δ 1

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

( ) ( )

s y s g s u

p

1 ) ( =

Inverse of gp(s) is unstable if gp(s) has a right-half-plane zero

( ) ( )

s u s g s y

p

= ) (

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

Matrix-vector Notation

( ) ( )

s u G s y

p

=

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )⎥

⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ s u s u s g s g s g s g s y s y

2 1 22 21 12 11 2 1

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

s u s g s u s g s y s u s g s u s g s y

2 22 1 21 2 2 12 1 11 1

+ = + = 2 input – 2 output Example Similar to SISO, controller serves as “process inverse”

( ) ( )

s y G s u

p 1 −

=

Steady-state Implications

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = 1 1 1 1

1 p

K

2 1 2 2 1 1

1 1 1 1 u u y u u y + − = + =

System 1 Which system will be more difficult to control? System 2

2 1 2 2 1 1

1 1 1 1 u u y u u y + = + = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 1 1 1 1

2 p

K

y K u

p 1 −

=

Steady-state Implications

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = 1 1 1 1

1 p

K

2 1 2 2 1 1

1 1 1 1 u u y u u y + − = + =

System 1 Inverse does not exist (Gain Matrix is singular, rank = 1) Cannot independently control both outputs System 2

2 1 2 2 1 1

1 1 1 1 u u y u u y + = + = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 1 1 1 1

2 p

K

y K u

p 1 −

=

Another Example

2 1 2 2 1 1

1 1 95 . 1 u u y u u y + = + = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 1 1 95 . 1

3 p

K Inverse exists – No longer singular (“gut feeling” that there may still be a problem…)

Directional Sensitivity: Singular Value Decomposition (SVD) G = UΣV

T

T V U G

                          ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − − ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − − = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡

Σ

72 . 70 . 70 . 72 . 0253 . 98 . 1 70 . 72 . 72 . 70 . 1 1 95 . 1 Gain matrix

Left singular Vector matrix Right singular Vector matrix Matrix of singular values

Example

Strongest input direction Strongest

• utput

direction Max singular value

Example, cont’d

T V U G

                          ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − − ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − − = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡

Σ

72 . 70 . 70 . 72 . 0253 . 98 . 1 70 . 72 . 72 . 70 . 1 1 95 . 1 ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ 41 . 1 38 . 1 70 . 72 . 1 1 95 . 1

2 1

y y ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ 70 . 72 .

2 1

u u ⎥ ⎦ ⎤ ⎢ ⎣ ⎡− = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ 72 . 70 .

2 1

u u ⎥ ⎦ ⎤ ⎢ ⎣ ⎡− = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡− ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ 018 . 018 . 72 . 70 . 1 1 95 . 1

2 1

y y

Input in Strongest Direction Output in Strongest Direction Input in Weakest Direction Output in Weakest Direction

Same magnitude input = very different magnitude output

y = Kp u

Example, cont’d

T V U G

                          ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − − ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − − = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡

Σ

72 . 70 . 70 . 72 . 0253 . 98 . 1 70 . 72 . 72 . 70 . 1 1 95 . 1 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

G1 s

( ) =

2.6 62s + 1 1.5 23s +1

( ) 62s +1 ( )

1.4 30s +1

( ) 90s +1 ( )

2.8 90s +1

( )

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

G2 s

( ) =

1.5 63s +1 2.5 39s +1

( ) 63s +1 ( )

2.5 56s +1

( ) 91s +1 ( )

1.6 91s +1

( )

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

z = −0.057 and + 0.013 sec

• 1

z = −0.060 and − 0.018 sec

• 1

Operating Point 1 – Minimum Phase Operating Point 2 – Nonminimum Phase (RHPT zero) Matrix inverse is unstable “Transmission zeros” are both negative

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

• f 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)

Chemical and Biological Engineering

Discrete Linear Models

• B. Wayne Bequette
• Sampling Rules/Assumptions
• Continuous to Discrete
• Z-transform
• Dynamic Properties of Discrete Systems

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

k k k k k k

Du Cx y u x x + = Γ + Φ =

+1

State Space Input-Output (ARX)

m k m k k k n k n k k k

u b u b u b u b y a y a y a y

− − − − − −

+ + + + + − − − − =  

2 2 1 1 2 2 1 1

( )

[ ]

k

y Z z y =

Z-transform

[ ]

( )

z y z y Z

k 1 1 − − =

Backwards shift operator

( ) ( ) ( ) ( )

z u z b z b z b b z y z a z a z a

m m m m n n n n − + − − − − + − − −

+ + + + = + + + +

1 1 1 1 1 1 1 1

... ... 1

Solve for y(z) usually b0 = 0

Some texts/papers have different sign conventions

Discrete Models

Discrete Transfer Function Multiply by

( ) ( )

z u z a z a z a z b z b z b b z y

n n n n m m m m − + − − − − + − − −

+ + + + + + + + =

1 1 1 1 1 1 1 1

... 1 ... gp z

( ) = b0 + b

1z −1 +... + b m−1z − m+1 + bmz − m

1 + a

1z−1 + ... + an−1z− n+1 + anz− n

gp z

( ) =

b0z

n + b 1z n−1 + ... + bmz n−m

zn + a

1z−1 + ... + an −1z −n+1 + an n n

z z

zeros = roots of numerator polynomial poles = roots of denominator polynomial ( )

( )( ) ( ) ( )( ) ( )

n m pz p

p z p z p z z z z z z z K z g − − − − − − ⋅ =  

2 1 2 1

Stability

Discrete Continuous

Stable Stable

poles in LHP are stable poles inside unit circle are stable

Example

y k

( ) = −a1y k −1 ( )+ b

1u k −1

( )

gp z

( ) =

b

1z −1

1 + a1z −1 = b

1

z + a1

p = −a1 a1 = 0.5,−0.5,1.5,−1.5 p = −0.5,0.5,−1.5,1.5

a1 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 behavior O s c i l l a t o r y, stable M o n o t o n i c , stable O s c i l l a t o r y, unstable M o n o t o n i c , unstable

Let y(0) = 1, u(k) = 0

Value of the pole

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

g p s

( )=

1 s +1

( )

3

( ) ( )( ) ( )

3

3679 . 1238 . 7990 . 1 0803 . − + + = z z z z g p

Sample time = 1

poles = 0.3679 (multiplicity 3) zeros = -1.7990 & -0.1238 Unstable inverse Relative order = 1 Relative order = 3 poles = -1 (multiplicity 3)

Final & Initial Value Theorems

lim

n→∞ y nΔt

( ) = lim

z→1 1− z −1

( )y z

( )

y z

( ) = gp z ( )u z ( ) = gp z ( )⋅

1 1− z−1 lim

t→ ∞ y t

( )= lim

z →1 1 − z −1

( )

y z

( ) = lim

z →1 1 − z −1

( )

gp z

( )

1 1− z

−1 = lim z →1 gp z

( )

lim

t→0 y t

( ) = lim

z →∞ 1− z −1

( )y z

( )

Initial value theorem Final value theorem Long-term step response unit step So, simply set z = 1 in gp(z) for long-term unit step response

“hidden slide” provided for additional background

Final Value Theorems

( ) ( )

3

1 1 + = s s g p

( ) ( )( ) ( )

3

3679 . 1238 . 7990 . 1 0803 . − + + = z z z z g p

( ) ( )( ) ( )

1 2525 . 2525 . 3679 . 1 1238 . 1 7990 . 1 1 0803 . 1

3

= = − + + = → z g p

( ) ( )

1 1 1

3 =

+ = → s g p

Long-term outputs for unit step inputs continuous discrete

“hidden slide” provided for additional background

State Space Models

Du Cx y Bu Ax x + = + = 

k k k k k k

Du Cx y u x x + = Γ + Φ =

+1

Continuous ? Discrete Finite differences approximation for derivative

( ) ( ) ( )

[ ]

k k k k k k k k k

tu B x t A I x Bu Ax t x x t x x t t k x t k x x Δ + Δ + = + ≈ Δ − Δ − = Δ Δ − Δ + ≈

+ + + 1 1 1

1  Φ Γ

How good are the approximations?

Exact Discretization

( ) ( ) ( )

( )

k t A k t A k t t t k A t A k t A k

Bu A I e x e x t Bu d e e t x e t t x

k k

1 1 − Δ Δ + Δ + − Δ Δ

− + = + = Δ +

∫

σ

σ

( )

B A I e e

t A t A 1 − Δ Δ

− = Γ = Φ

[ ]

t B t A I Δ = Γ Δ + = Φ Exact Approximate

Example Discretization

( )

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − − Φ = Γ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − = = Φ

− Δ

0157 . 2592 . 1 . 04 . 04 . 1 . 8869 . 0974 . 7408 . 12 . 12 . 3 . exp

1

I e

t A

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = Γ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = Φ 3 . 88 . 12 . 7 . 12 . 12 . 3 . 1 1 Exact Approximate

[ ]

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡

2 1 2 1 2 1

1 1 . 04 . 04 . 1 . x x y u x x x x  

( ) ( )( )

1 25 1 10 1 + + = s s s g p

3 = Δt

Example, Continued

Exact Approximate

( ) ( )

D zI C z g p + Γ Φ − =

−1

( )

2 1 2 1 2

65702 . 6277 . 1 1 0136 . 0157 . 65702 . 6277 . 1 0136 . 0157 .

− − − −

+ − + = + − + = z z z z z z z z g p

( )

2 1 2 2

616 . 58 . 1 1 036 . 616 . 58 . 1 036 .

− − −

+ − = + − = z z z z z z g p

Poles/eigenvalues = 0.7408 & 0.8869 Poles/eigenvalues = 0.7 & 0.88

Discrete Transfer Function

“hidden slide” provided for additional background

• 2

2 4 6 8 0.5 1 1.5 2 temp, deg C discrete-time step

• 2

2 4 6 8 0.2 0.4 0.6 0.8 1 discrete-time step psig

s

1

s s s s

2 4 5 3 Etc.

S = s1 s2 s3 s4 s5  sN ! " # $

T

Example Step Response Model

Step Response Model

∞ − ∞ + − − + − − ∞ = −

Δ + + Δ + Δ + + Δ = Δ =∑

k N N k N N k N k i i k i k

u s u s u s u s u s y  

1 1 1 1 1

( )

. ,

1 1 1 1 1 1 1 1 1 1

∑

− = − − ∞ − − + − − − ∞ − − + − − −

Δ + = Δ + + Δ + Δ + + Δ = Δ + + Δ + Δ + + Δ =

−

N i i k i N k N k u k N k N N k N k k N N k N N k N k k

u s u s y u u s u s u s u s u s u s u s y

N k

            

• 2

2 4 6 8 10 0.2 0.4 0.6 0.8 1 discrete-time step psig

• 2

2 4 6 8 10

• 0.2

0.2 0.4 0.6 0.8 1 discrete-time step temp, deg C

h h h 1 2 3

hi = si − si−1 si = hj

j=1 i

∑

Example Impulse Response Model

Impulse and step response coefficients are related

Parameter Estimation, ARX Models

2 2 1 1 2 2 1 1 2 1 1 2 1 1 2 1 2 1 1 2 1 1

ˆ ˆ ˆ

− − − − − −

+ + − − = + + − − = + + − − =

N N N N N

u b u b y a y a y u b u b y a y a y u b u b y a y a y  ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡

− − − − − − 2 1 2 1 2 1 2 1 1 1 1

. . . . . . . . ˆ . . ˆ b b a a u u y y u u y y y y

N N N N N

Θ Φ = ˆ Y

model output measured output known (measured) input time step index N pts could change sign on columns 1 & 2 of Φ

Optimization/Parameter Solution

yi − ˆ yi

( )

2 i=1 N

∑

= Y − ˆ Y

( )

T Y − ˆ

Y

( ) = Y −ΦΘ

( )

T Y −ΦΘ

( )

( )

∑

=

−

N i i i b b a a

y y

1 2 , , ,

ˆ min

2 1 2 1

Θ = Φ

TΦ

( )

−1Φ TY

Objective Function Solution

parameter estimate vector measured

• utput

vector

Example

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

• 0.2

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

• 1
• 0.5

0.5 1 time, min u

PRBS input

Result

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − − − = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − − − − − − − − − − − − − − − − − − − − − − − − − − − = Φ 0502 . 1134 . 1162 . 1513 . 0172 . 3846 . 5007 . 2003 . 0463 . 1123 . 1879 . 0485 . 1558 . 2171 . 3446 . 1771 . 4618 . 1564 . 0137 . 1 1 1162 . 1134 . 1 1 1513 . 1162 . 1 1 0172 . 1513 . 1 1 3846 . 0172 . 1 1 5007 . 3846 . 1 1 2003 . 5007 . 1 1 0463 . 2003 . 1 1 1123 . 0463 . 1 1 1879 . 1123 . 1 1 0485 . 1879 . 1 1 1558 . 0485 . 1 1 2171 . 1558 . 1 1 3446 . 2171 . 1 1 1771 . 3446 . 1 1 4618 . 1771 . 1 1 1564 . 4618 . 1 1 0137 . 1564 . 1 1 0889 . 0137 . 1 1 0889 . Y

Θ = Φ

TΦ

( )

−1Φ TY

Θ = 1.1196 −0.3133 −0.0889 0.2021 ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ = −a1 −a2 b1 b2 ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥

( ) ( ) ( )( )

5481 . 5716 . 274 . 2 0889 . 3133 . 1196 . 1 1 2021 . 0889 . 3133 . 1196 . 1 2021 . 0889 .

2 1 2 1 2 2 1 2 2 1

− − − − = + − + − = + − + − = + + + =

− − − −

z z z z z z z z z z a z a z b z b z g p

Identified Model

( ) ( ) ( )

3133 . 1196 . 1 1 2021 . 0889 . 3133 . 1196 . 1 2021 . 0889 .

2 1 2 1 2 2 1 2 2 1

z u z z z z z y z z z a z a z b z b z g p ⋅ + − + − = + − + − = + + + =

− − − −

2021 . 0889 . 3133 . 1196 . 1

2 1 2 1 − − − −

+ − − =

k k k k k

u u y y y

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

Chemical and Biological Engineering

Discrete (Digital) Control

• B. Wayne Bequette
• Review of Digital PID
• Review of Model-based Digital Control
• Discrete Internal Model Control
• Examples
• Summary

Proportional-Integral-Derivative (PID) Control

( ) ( ) ( ) ( ) ( ) ( ) ( )

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + + + = − =

∫

dt t de dt t e t e k u t u t y t r t e

D t I c

τ τ 1

Manipulated Input Error = setpoint – measured output Proportional gain Integral time Derivative time

e k

( )= e tk

( )

( ) ( ) ( ) ( ) ( ) ( )

∑ ∑ ∫

= =

Δ = Δ ≈ Δ ⋅ + + Δ ⋅ + Δ ⋅ ≈

k i k i i k t

t i e t t e t t e t t e t t e dt t e

k

1 1 2 1

 de tk

( )

dt ≈ e k

( )− e k −1 ( )

Δt

Substituting each term, we find the discrete controller

( ) ( ) ( ) ( ) ( )⎥

⎦ ⎤ ⎢ ⎣ ⎡ − − − Δ + Δ + − + = −

∑

− = 1

2 1 ) ( 1 1

k i D I c

k e k e t i e t k e k u k u τ τ

u k

( ) = u

0 + kc e k

( )+ Δt

τI e(i) + τ D Δt e k

( )− e k −1 ( )

( )

i= 0 k

∑

⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ u k

( ) = u k −1 ( )+ kc

1+ Δt τ I + τ D Δt ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ e k

( )+ −1− 2τ D

Δt ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ e k −1

( )+ τ D

Δt e k − 2

( )

⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥

It is convenient to work with the “velocity form.” First, generate the equation for step k-1 Then subtract u(k-1) from u(k) to find

u k

( ) = u k −1 ( )+ b2e k ( )+ b

1e k −1

( )+ b0e k − 2 ( )

b2 = kc 1+ Δt τ I + τ D Δt " # $ % & ', b

1 = −kc 1+ 2τ D

Δt " # $ % & ', b0 = kcτ D Δt

u k

( ) = u k −1 ( )+ kc

1+ Δt τ I + τ D Δt ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ e k

( )+ −1− 2τ D

Δt ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ e k −1

( )+ τ D

Δt e k − 2

( )

⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥

( )

[ ]

( ) ( )

[ ]

( ) ( )

[ ]

( ).

2 , 1 ,

2 1

z e z k e z e z k e z e k e

− −

= − Ζ = − Ζ = Ζ

( ) ( ) ( ) ( ) ( )

2 1 1

1 2

− + − + = − − k e b k e b k e b k u k u ( )

[ ]

( ) ( )

[ ]

( ) ( )

[ ]

( ).

2 , 1 ,

2 1

z u z k u z u z k u z u k u

− −

= − Ζ = − Ζ = Ζ

( ) ( ) ( ) ( )

z e z b z b b z u z

2 1 1 2 1

1

− − −

+ + = −

( ) (

) ( )

z e z z b z b b z u

1 2 1 1 2

1

− − −

− + + =

( ) (

) ( )

( ) ( )

z e z g z e z z b z b z b z u

c

= − + + =

2 1 2 2

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

l Continuous Model l PID Parameters

( )

( )

. 2 , 5 . , 5 . 5 . θ τ θ τ τ θ τ τ θ λ θ τ + = + = + + =

p p D p I p p c

k k  gp(s) = kpe−θs τ ps+1

λ is the tuning parameter – related to the closed-loop time constant

Digital Control: Block Diagram

• r(z)
• e(z)
• u(z)
• gc(z)
• gp(z)
• d(z)
• y(z)
• +
• +
• +
• _

l Stability analysis l Poles of CLTF must be

inside the unit circle

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

z r z g z y z r z g z g z g z g z y

CL c p c p

1 = ⋅ + =

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

r(z) e(z) u(z) gc(z) gp(z) d(z) y(z) + + + _

l Digital Controller Design

Ø Response specification-based

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

z r z g z y z r z g z g z g z g z y

CL c p c p

1 = ⋅ + =

Desired response

Controller Design

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

z g z g z g z g z g z g z g z g z g

CL CL p c c p c p CL

− ⋅ = + = 1 1 1 Desired response solve for controller

Internal Model Form

u(z) q(z) r(z) d(z) y(z) + + +

• + -

~ g (z)

p

d(z) ~ r(z) ~ g (z)

p

y(z) ~ Process Process model IMC controller

For controller design, consider perfect model and no disturbances

IMC Design

r(z) u(z) q(z) gp(z) y(z)

( ) ( ) ( ) ( ) ( ) ( ) ( )

z r z g z y z r z q z g z y

CL p

= =

Desired response

( ) ( ) ( )

z g z g z q

p CL

=

solve for controller Really based on the model

( ) ( ) ( ) ( ) ( )

z g z g z g z g z q

p CL p CL 1

~ ~

−

= =

Can implement IMC in Classic Feedback Form

( ) ( ) ( ) ( )

z g z q z q z g

p c

~ 1− =

Digital Controller Design

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

John Ragazzini

Deadbeat

l Achieve setpoint in the minimum time (if no time-delay,

then one time step)

k k

r y =

+1

( ) ( )

z r z z y

1 −

=

Backwards shift operator

( ) ( )

z g z z q

p 1 1~− −

=

( )

1 −

= z z gCL IMC Form

( ) ( )

1 1 1

1 ~

− − −

− ⋅ = z z z g z g

p c

Classic Feedback Form

Example 0, FO Process

( )

1 10 1 ~ + = s s g p

( )

1 1

9048 . 1 0952 . ~

− −

− = z z z g p

Δt = 1

( )

( )

0952 . 9048 . 0952 . 9048 . 1

1 1 1

+ − = − =

− − −

z z z z z z q

( ) ( ) ( ) ( ) ( )

( )

( )

1 1 1 1 1 1 1 2 1 1

1 504 . 9 1 504 . 10 1 504 . 9 504 . 10 1 0952 . 9048 . step unit 1 1

− − − − − − − − − −

− − + − = − − = − − = = − = = z z z z u z z z z z z z u z z r let z r z q z u Large control action up, then down Control action setpoint change

Classic Feedback Form for Deadbeat Design ( ) ( ) ( )

1 1 1 1 1 1 1 1 1

1 9048 . 1 5 . 10 1 9048 . 1 0952 . 1 1 9048 . 0952 . 1 1 1 0952 . 9048 . 0952 . 9048 . 0952 . 9048 . 1 ~ 9048 . 0952 . 9048 . 1 0952 . ~

− − − − − − − − −

− − ⋅ = − − ⋅ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = − − ⋅ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = − ⋅ − = − = − = − = − = z z z z z z z z z g z z z z g z z z z g

c p p

( ) ( ) ( )

1 1 ~ 1 ~

1 1 1 1

− ⋅ = − ⋅ =

− − − −

z z g z z z g z g

p p c

Classic Feedback Implementation ( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( ) ( )

1 1 1 1 1 1 1 1

9048 . 5 . 10 9048 . 5 . 10 9048 . 1 5 . 10 1 1 9048 . 1 5 . 10

− − − − − − − −

− ⋅ + = − ⋅ = − − ⋅ = − = − − ⋅ =

k k k k k k k k c c

e e u u e e u u z e z z u z z e z g z u z z z g

Note that this is a PI controller

uk = uk−1 +10.5⋅ ek − 0.9048ek−1

( )

u k

( ) = u k −1 ( )+ b2e k ( )+ b

1e k −1

( )+ b0e k − 2 ( )

b2 = kc 1+ Δt τ I + τ D Δt " # $ % & ' =10.5 b

1 = −kc 1+ 2τ D

Δt " # $ % & ' = −10.5*0.9048 = −9.5 b0 = kcτ D Δt = 0 kc = 9.5 τ I = 9.5 The PI parameters are then

IMC Implementation

u(z) q(z) r(z) d(z) y(z) + + +

• + -

~ g (z)

p

d(z) ~ r(z) ~ g (z)

p

y(z) ~ Process Process model IMC controller

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

z r z q z u z d z r z r z y z y z d z u z g z y

p

~ ~ ~ ~ ~ ~ ~ = − = − = =

( ) ( )

( ) ( )

( ) ( )

( )

( ) ( )

( )

1 1 1 1 1 1 1 1 1

~ 9048 . ~ 0952 . 1 ~ 9048 . 1 0952 . 1 ~ 9048 . 1 0952 . ~ ~ ~ ~ 0952 . ~ 9048 . ~ 0952 . ~ 9048 . 1

− − − − − − − − −

− ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = − = − = − = + = = −

k k k k k k k k k k k k

r r u z r z z u z r z z z u z d r r y y d u y y z u z z y z

Discussion ( )

1 1 1

~ 9048 . ~ 0952 . 1 ~ ~ ~ ~ 0952 . ~ 9048 . ~

− − −

− ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = − = − = + =

k k k k k k k k k k k k

r r u d r r y y d u y y

( )

1 1

9048 . 5 . 10

− −

− ⋅ + = − =

k k k k k k k

e e u u y r e

Standard feedback control clearly has the current control action as a function of the previous control action. Why doesn’t IMC? Standard Feedback IMC

Example 0 (First-order), Deadbeat Design

• 2

2 4 6 8 0.5 1 time y deadbeat, FO example

• 2

2 4 6 8 5 10 time u

Example 1, Second-Order Process

( ) ( )( )

1 25 1 10 1 ~ + + = s s s g p

( )

2 1 2 1

65702 . 6277 . 1 1 0136 . 0157 . ~

− − − −

+ − + = z z z z z g p

Δt = 3

( )

0136 . 0157 . 65702 . 6277 . 1 0136 . 0157 . 65702 . 6277 . 1

2 2 2 1 3 2 1

+ + + − = + + − =

− − − − −

z z z z z z z z z z q

“ringing” behavior as shown next zero = -0.0136/0.0157 = -0.8694 controller has pole = -0.8694, which is stable, but oscillates

Example 1 (Second-order): Deadbeat Design

• 5

5 10 15 20 25 30 0.5 1 1.5

time y deadbeat, SO example

• 5

5 10 15 20 25 30

• 100
• 50

50 100

time u

Perfect control, at sample times “ringing” behavior (intersample)

Dahlin’s Controller

l Desired first-order + time-delay response to setpoint

change

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

z g z z g z z z g z g z q

p p p CL 1 1 1 1

~ 1 ~ 1 1 ~

− − − −

⋅ − − = ⋅ − − = = α α α α

( ) ( )

1 1

1 1

− −

− − = z z z gCL α α

where α = exp(-Δt/λ) For no time-delay:

( ) ( ) ( ) ( )

0136 . 0157 . 65702 . 6277 . 1 1 0136 . 0157 . 65702 . 6277 . 1 1 1 1

2 2 1 2 1 1 1

+ + − ⋅ − − = + + − ⋅ − − =

− − − − − −

z z z z z z z z z z z q α α α α For the second-order example:

Second-order Process: Dahlin’s Controller

Still “ringing”, but more damped than deadbeat For λ = 3, α = 0.3679

• 5

5 10 15 20 25 30 0.5 1 1.5 time y Dahlin, SO example

• 5

5 10 15 20 25 30

• 100
• 50

50 100 time u

First-order response, at sample times

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 gc(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

( )

0293 . 65702 . 6277 . 1 0293 . 65702 . 6277 . 1 1 ~

2 2 2 1 1

+ + + − = + − = =

− − − −

z z z z z z g z q

p

( )

0136 . 0157 . 0136 . 0157 . 65702 . 6277 . 1 1 0136 . 0157 . 65702 . 6277 . 1 1 0136 . 0157 . ~

2 1 2 1 2 1 2 1

+ + ⋅ + − + = + − + =

− − − − − − − −

z z z z z z z z z g p

( )

z g p− ~

( )

z g p+ ~ second-order example

Second-order Example: State Deadbeat Design

• 5

5 10 15 20 25 30 0.5 1 1.5 time y State Deadbeat, SO example

• 5

5 10 15 20 25 30

• 100
• 50

50 100 time u

Controller Forms

( ) ( )

( ) ( )( ) ( )

( ) ( )

n n k k pz p

p z p z a z a z a z a z K z G z g − − − − − − = =

+ − + + − −

  

1 1 1 1

~

Zafiriou & Morari notation Assumes stable poles

+ i

a

− i

a

zero inside unit circle = = zero outside unit circle

( ) ( ) ( )

( ) ( ) n

k pz n SD

z a a K p z p z z q

− −

− − − − = 1 1

1 1

 

( ) ( ) ( ) ( )

( ) ( )

( )

1 1 1 1

1 1 1 1

− − − −

− − ⋅ − − − − = = z z z a a K p z p z z F q z q

n k pz n SD VE

α α   State Deadbeat: State Deadbeat with Filter (Vogel-Edgar):

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”)

( ) ( ) ( )

z g z g z g

p p p + −

= ~ ~ ~

• zeros outside unit circle
• zeros inside unit circle that are negative
• “all-pass” by including pole at 1/zi for each

positive zi outside the unit circle (but not the negative ones)

• numerator one order

less than denominator

( ) ( ) ( )

z F z g z q

p 1

~−

−

=

l Controller: Invert “good stuff” part of model

• “Filter” for desired response, often first-order

Example 2: Third-Order (3-tank)

( ) ( )

3

1 1 ~ + = s s g p

( ) ( )( ) ( )

3

3679 . 1238 . 7990 . 1 0803 . ~ − + + = z z z z g p

( ) ( )( ) ( ) ( )( ) ( )( ) 2

3 2

1238 . 1 7990 . 2 1238 . 7990 . 1 3679 . 1238 . 1 7990 . 2 0803 . ~ z z z z z z g p + + ⋅ − =

( ) ( ) ( )

z g z g z g

p p p + −

= ~ ~ ~

Sample time = 1

( ) ( ) ( ) ( ) ( ) ( ) ( )

α α α α − − ⋅ − = − − ⋅ − = =

− − − −

z z z z z z z z F z g z q

p

1 2526 . 3679 . 1 1 2526 . 3679 . ~

2 3 1 1 2 3 1

Note gain is 1 (value at z = 1) where α = exp(-Δt/λ) zeros at -1.7990 (outside unit circle)

• 0.1238 (inside, but negative)

Example 2 (Third-Order) IMC Design

• 2

2 4 6 8 10 0.5 1 1.5 time y Discrete IMC, TO example

• 2

2 4 6 8 10

• 1

1 2 3 4 5 time u

For λ = 1, α = 0.135

Third-Order, Inverse Response (Ex. 3, Z&M, 1985)

( ) ( ) ( )( )( )

5 . 2 2 1 5 . 1 333 . 3 ~ + + + + − = s s s s s g p

( ) ( )( ) ( )( )( )

779 . 819 . 905 . 792 . 162 . 1 01316 . ~ − − − + − − = z z z z z z g p

( ) ( )( )( ) ( )( )( ) ( )( ) ( )( )

792 . 1 162 . 162 . 1 1 162 . 1 1 1 792 . 162 . 1 779 . 819 . 905 . 162 . 1 1 1 162 . 1 1 792 . 1 162 . 01316 . ~ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − + − ⋅ − − − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − − − = z z z z z z z z z z g p

Δt = 0.1

( ) ( ) ( ) ( )( )( ) ( ) ( )

α α − − ⋅ − − − − = =

− −

z z z z z z z F z g z q

p

1 8606 . 02741 . 779 . 819 . 905 . ~ 1

where α = exp(-Δt/λ) zeros at 1.162 (outside unit circle)

• 0.792 (inside, but negative)

( ) ( ) ( )( )( )

779 . 819 . 905 . 8606 . 02741 . ~ − − − − =

−

z z z z z z g p

Response (Ex. 3, Z&M ’85)

• 1

1 2 3 4 5

• 1
• 0.5

0.5 1 1.5 time y Discrete IMC, RHPZ - Example 3 Z&M

• 1

1 2 3 4 5

• 2

2 4 6 8 10 time u

For λ = 0.5, α = 0.8187

“Real-World” Discussion

l Have assumed a “perfect model” for these simulations. In

practice, the real-world “plant” is not perfectly modeled (indeed, there are usually nonlinearities involved).

l Can approximate real-world challenges by having the

“plant” be different than the model used for controller

• design. Also, it is important to incorporate constraints

and noise in the simulations.

Summary

l Digital Control Techniques

Ø Deadbeat Ø Dahlin’s Ø State Deadbeat, State Deadbeat w/filter (Vogel-Edgar) Ø Modified Dahlin’s (mis-developed by Dahlin!)

l Internal Model Control

Ø Factorization of zeros outside unit circle and negative zeros inside the

circle

Ø Form “all-pass”(reflection of positive zeros outside the unit circle) Ø Invert “good stuff” and use first-order filter

Suggestion

l Work through Workshop on Discrete IMC