Pros and Cons IIT Bombay In each case, we get highly nonlinear - - PowerPoint PPT Presentation

Automation Lab IIT Bombay

Pros and Cons

In each case, we get highly nonlinear optimization

problem

Important to give good initial guess to obtain

reliable model parameters: not an easy task

It is assumed that ‘structure’ of the linearized

state space model / q-TF matrix is know a-priori, which may not be the case

Effect of unmeasured inputs (disturbances) on

system dynamics is not captured by the model

Possible Remedy: Black box modeling 4/10/2012 System Identification 63

Automation Lab IIT Bombay

Summary

If a reliable mechanistic / grey box model

is available, then a linear perturbation model can be developed using the Taylor series expansion in the neighborhood of the operating point

If structure of linear perturbation model

matrices is known, then the parameters can

• f the model can be estimated using input
• utput data.

4/10/2012 System Identification 64

Automation Lab IIT Bombay

Stability of Linear Dynamic Model

4/10/2012 System Identification 65

Model

• n

Perturbati Linear Time Continuous SISO Consider Results Time Continuous

• f

Review Quick

[ ]

solution the

• f

behavior asymptotic the determine matrix

• f

es eigen valu i.e. det ) ( polynomial r denominato

• f

Roots A A I − = s s d

[ ]

( )

[ ] [ ] [ ]

A I B A I C B A I C − − = − = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ =

−

s s adj s d s n s u s s u s d s n s y det ) ( ) ( ) ( ) ( ) ( ) ( ) ( Function Transfer

1

) ( = = + = x Cx B Ax x y u dt d

Automation Lab IIT Bombay

Stability of Linear Dynamic Model

4/10/2012 System Identification 66

[ ]

Stable Marginaly is System Dynamic LHP in rest and axis imaginary

• n the

are roots the

• f

some If Stable ally Asymptotic is System Dynamic plane

• s

complex

• f

Plane Half Left strictly in are det ) ( polynomial r denominato

• f

roots the All : Case ⇒ ⇒ − = A I s s d

[ ]

Unstable is System Dynamic plane

• s

complex

• f

Plane Half Right in are det ) ( polynomial r denominato

• f

root nay Is : Case ⇒ − = A I s s d

Automation Lab IIT Bombay

Stability: Discrete Time Systems

4/10/2012 System Identification 67

Consider SISO discrete time linear perturbation model

[ ]

( )

[ ] ( ) [ ]

Φ I Γ Φ I C Γ Φ I C

• z

det

• z

) ( ) ( ) (

• z

) ( ) ( (z)

1

adj z d z n z u z d z n y = = =

−

x Cx Γ Φx x = = + = + ) ( ) ( ) ( ) ( ) ( ) 1 ( k k y k u k k

[ ]

solution the

• f

behavior asymptotic the determine matrix

• f

s eigenvalue i.e. det ) ( polynomial r denominato

• f

Roots case time continuous the to Analogous Φ Φ − = I z z d

Automation Lab IIT Bombay

Stability: Discrete Time Systems

4/10/2012 System Identification 68

) ( ) ( .... .......... ) ( ) 1 ( ) 2 ( ) ( ) 1 ( u(k) and zero

• non

is ) ( scenario Consider

2

x Φ x x Φ Φx x Φx x x

k

k = = = = =

1 1 2 2 1

able diagonaliz is where case special Consider

− − −

Ψ ΨΛ = Φ ⇒ Ψ ΨΛ = Φ ⇒ ΨΛΨ = Φ Φ

k k

( ) ( ) ( ) ⎥

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = Λ

k n k k

λ λ λ ... ... ... ... ... ... ... ...

2 1

Automation Lab IIT Bombay

Stability of Linear Dynamic Model

4/10/2012 System Identification 69

[ ]

( )

[ ] [ ]

Stable ally Asymptotic is System Dynamic ) ( ) ( k as Then n ... 1,2, i for i.e. plane

• z

complex in circle unit the inside strictly in are

• f

s eigenvalue i.e. det ) ( polynomial r denominato

• f

roots the All : Case

k k i i

⇒ → Φ = ⇒ → Ψ ΨΛ = Φ ⇒ → Λ ⇒ ∞ → → = < Φ Φ − =

−

1

1

x x I

k k k

k z z d λ λ

Automation Lab IIT Bombay

Stability of Linear Dynamic Model

4/10/2012 System Identification 70

[ ]

( )

Stable Marginally is System Dynamic k as bounded remains ) ( ) ( k as bounded remains k as bounded remains Then n ... 1,2, i for i.e. plane

• z

complex in circle unit the

• n
• r

inside in are

• f

s eigenvalue i.e. det ) ( polynomial r denominato

• f

roots the All : Case

k i i

⇒ ∞ → Φ = ⇒ ∞ → Ψ ΨΛ = Φ ⇒ ∞ → = ≤ Φ Φ − =

−

1

1

x x I

k k k

k z z d λ λ

Automation Lab IIT Bombay

Stability of Linear Dynamic Model

4/10/2012 System Identification 71

[ ]

( )

[ ] [ ]

Unstable is System Dynamic Linear ) ( ) ( k as Then j some for i.e. plane

• z

complex

• f

circle unit the inside

• utside

are

• f

s eigenvalue i.e. det ) ( polynomial r denominato

• f

roots more

• r
• ne

any If : Case

k k j j

⇒ ∞ → Φ = ⇒ ∞ → Ψ ΨΛ = Φ ⇒ ∞ → Λ ⇒ ∞ → ∞ → > Φ Φ − =

−

1

1

x x I

k k k

k z z d λ λ

Automation Lab IIT Bombay

Poles in s-plan and z-plane

4/10/2012 System Identification 72

[ ]

Φ Φ Ψ Ψ = = Φ ⇒ ΨΛΨ =

Λ

• f

eigenvalue is e then A

• f

eigenvalue is if Also, rs eigenvecto identical have Matrix and Matrix that implies This

T

• 1
• 1

λ

λ A A

A T T

e e

[ ] [ ]

T

e by determined is plane

• z

in

• f

Location ) sin( ) cos( ) sin( ) cos(

• f

eigenvalue an represent Let

α λ α β α λ

β β β β β α λ

T T T j T T

e T j T T j T e e e j ⇒ = + + = = + =

+

1 A

Automation Lab IIT Bombay

Poles in s-plan and z-plane

4/10/2012 System Identification 73

plane

• z

in

• rigin

the around circle unit a inside map plane

• s
• f

half left in poles All e then If

T

⇒ < < 1

α

α plane

• z

in circle unit the

• f

boundary the

• n

map plane

• s

in axis imaginary the

• n

poles All e then If

T

⇒ = = 1

α

α plane

• z

in circle unit the e

• utsidesid

map plane

• s
• f

half left in poles All e then If

T

⇒ > > 1

α

α

Automation Lab IIT Bombay

Poles in s-plan and z-plane

4/10/2012 System Identification 74

Conformal map z = exp(sT)

Automation Lab IIT Bombay

Poles in s-plan and z-plane

4/10/2012 System Identification 75

possible. always be not may reverse the However, specified. is terval samplingin the

• nce

system time discrete equivalent an find can we system linear time continuous every For : Note

5 . ) exp(-T/ such that

• f

value real with ) ( ) ( dt dy form the

• f

system time continuous no is There 0.5

• z

: pole domain 5 1 Example − = = + = − + = + τ τ κ τ t u t y z u(k) y(k) .

• )

y(k

Automation Lab IIT Bombay

4/10/2012 System Identification 76

Reference Textbooks

Astrom, K. J., and B. Wittenmark, Computer

Controlled Systems, Prentice Hall India (1994).

Franklin, G. F., Powell, J. D., and M. L. Workman,

Digital Control Systems, Addison Wesley, 1990.

Ljung, L., Glad, T., Modeling of Dynamic Systems,

Prentice Hall, N. J., 1994.

Development of Control Relevant Linear Perturbation Models

Part II: Development of Black Box Models

Sachin C. Patwardhan

• Dept. of Chemical Engineering, IIT Bombay

Email: sachinp@iitb.ac.in

Automation Lab IIT Bombay

System Identification 2 4/10/2012 System Identification 2

Data Driven Models

Development of linear state space/transfer models starting from first principles/gray box models is impractical proposition in many situations. Practical Approach

• Conduct experiments by perturbing process

around chosen operating point

• Collect input-output data
• Fit a differential equation or difference

equation model using optimization Difficulties

• Measurements are inaccurate
• Process is influenced by unknown disturbances
• Models are approximate

Automation Lab IIT Bombay

System Identification 3 4/10/2012 System Identification 3

Discrete Model Development

2 4 6 8 10 12 14 16 18 20 2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 Sampling Instant Manipulated Input

Excite plant around the desired operating point by injecting input perturbations Process

5 10 15 20 1.8 2 2.2 2.4 2.6 2.8 3 3.2

Sampling Instant Measured Output

Input excitation for model identification Unmeasured Disturbances Measured output response Measurement Noise

Automation Lab IIT Bombay

System Identification 4

Black Box Models

Black Box Models: Static maps (correlations) and

• r

dynamic models (difference equations) developed directly from historical input-output data

Valid over limited operating range covered by the

training data set

Provide NO insight into internal working of system

under consideration

Development process: much less time consuming

and comparatively easy. Does not need services of a domain expert

4/10/2012 System Identification 4

Automation Lab IIT Bombay

System Identification 5

Two Non-Interacting Tanks Setup

4/10/2012 System Identification 5

Tank 1 Tank 2 Control Valve 1 Control Valve 2

Automation Lab IIT Bombay

System Identification 6 4/10/2012 System Identification 6

Two Non-Interacting Tanks Setup

Manipulated Input

Pump1 Tank 2 Tank 1

Non Interacting Tank Level Control setup

Sump Pump2 CV 2 CV 1

LT

Disturbance Input SISO System Controlled Output: Level in Tank 2 Sampling Interval T = 5 sec

Automation Lab IIT Bombay

System Identification 7 4/10/2012 System Identification 7

Input Output Data

50 100 150 200 250 4.5 5 5.5 6 6.5

Output (mA) Raw Input and output signals

50 100 150 200 250 10 12 14

Time Input (mA)

Automation Lab IIT Bombay

System Identification 8 4/10/2012 System Identification 8

Perturbation Data for Identification

50 100 150 200 250

• 1

1

Output Input and output signals

50 100 150 200 250

• 2

2

Time Input

Mean values removed from Input and Output data

Automation Lab IIT Bombay

System Identification 9 4/10/2012 System Identification 9

Black Box Models

Dynamic Models: Given observed data

[ ] [ ]

) ( ... ) 2 ( ) 1 ( : Outputs Measured ) ( ... ) 2 ( ) 1 ( : Inputs past

• f

Set

) ( ) (

k y y y Y k u u u U

k k

= =

we are looking for relationship

( )

) ( , , ) (

) 1 ( ) 1 (

k e Y U k

k k

+ Ω =

− −

θ y

such that noise (residuals) e(k) are as small as possible

vector parameter represents

p

R ∈ θ

{ } { }

) ( ),....., 1 ( ), ( U : Set Input d Manipulate ) ( ),....., 1 ( ), ( Y : Set Outputs Measured ) variables (deviation Collected Data al Experiment

N N

N N u u u y y y ≡ ≡

Automation Lab IIT Bombay

System Identification 10 4/10/2012 System Identification 10

Linear Difference Equation Model

Output Error (OE) Model

) ( ) ( ) ( ) ( ) ( ) ( k v k x k y k u q G k x + = =

Deterministic component Residue: unmeasured dist.+ measurement noise

N Nq

h q h q h q G q a q a q a q b q b q b q G q a q a q b q b q G

− − − − − − − − − − − − − − − −

+ + + = + + + + + = + + + = .... ) ( 1 ) ( 1 ) ( ) G(q for Choices Possible

2 2 1 1 1 3 3 2 2 1 1 3 3 2 2 1 1 1 2 2 1 1 2 2 1 1 1

• 1

Finite Impulse Response (FIR) Model

Automation Lab IIT Bombay

System Identification 11 4/10/2012 System Identification 11

Parameterized OE model

) ( ) ( ) ( ) 3 ( b ) 2 ( ) 2 ( ) 1 (

• a

(k) to equivalent is which ) ( q a + q a + 1 q b q b ) ( period sampling 1 to equal found was time) (dead delay Time

2 1 2 1 1 2 2 1

• 1

2

• 2

1

• 1

k v k x k y k- u k- u b k- x a k- x x k u q k x + = + + − = + =

− −

Proposed Model: System is expected to have second

• rder dynamics if we neglect valve and sensor dynamics

) ( 1 ) ( model time discrete

• rder

nd 2' to equivalent is which ) ( ) 1 )( 1 ( ) (

2 2 1 1 2 2 1 1 2 1

k u q a q a q b q b k x s u s s k s x

p − − − −

+ + + = + + = τ τ

Automation Lab IIT Bombay

System Identification 12 4/10/2012 System Identification 12

Model Parameter Estimation

Difficulty with signals in model x(k) : contribution to dynamics due to u(k) alone (cannot be measured) v(k) : measurement noise / disturbance (cannot be measured) Consequence Parameter estimation problem has to be solved using nonlinear optimization problem

Automation Lab IIT Bombay

System Identification 13 4/10/2012 System Identification 13

Model Parameter Estimation

) ( ) 1 ( ) 1 ( ) 2 ( ) 3 ( ˆ follows as ) ( estimate y recursivel can we )) 2 ( ), 1 ( , , , , ( for guess Given

2 1 2 1 2 1 2 1

u b u b x a x a x k x x x b b a a + + − − =

) 3 ( ) 2 ( ) 2 ( ˆ ) 1 ( ˆ ) ( ˆ

2 1 2 1

− + − + − − − − = N u b N u b N x a N x a N x

) 1 ( ) 2 ( ) 2 ( ) 3 ( ˆ ) 4 ( ˆ

2 1 2 1

u b u b x a x a x + + − − =

.......... ) 2 ( ) 3 ( ) 3 ( ˆ ) 4 ( ˆ ) 5 ( ˆ

2 1 2 1

u b u b x a x a x + + − − =

N k for k x k y k v ,.... 4 , 3 ) ( ˆ ) ( ) ( ˆ = − =

Automation Lab IIT Bombay

System Identification 14 4/10/2012 System Identification 14

Parameter Estimation

[ ] [ ]

) 3 ( ) 2 ( ) 2 ( ˆ ) 1 ( ˆ ) ( ˆ ) ( ˆ ) ( ) ( ˆ )) 2 ( ), 1 ( , , , , ( respect to with minimized is (k) v ˆ ) ( v ˆ ),... 2 ( v ˆ such that )) 2 ( ), 1 ( , , , , ( Estimate

2 1 2 1 2 1 2 1 3 2 2 1 2 1

− + − + − − − − = − = = Ψ

∑

=

k u b k u b k x a k x a k x k x k y k v x x b b a a k x x b b a a

N k

Possible Simplification : Choose x(0) = x(1) = 0

Nonlinear Optimization Problem

q 0.6841 q 1.653

• 1

) A(q q 0.01269 q 006

• 4.567e

) B(q ) Toolbox tion Identifica System MATLAB (using Parameters Model Identified

2

• 1
• 1
• 3
• 2
• 1
• +

= + =

Automation Lab IIT Bombay

System Identification 15 4/10/2012 System Identification 15

OE Model for Two Tank System

50 100 150 200

• 0.5

0.5

y(k) OE(2,2,2): Measured and Simulated Outputs

50 100 150 200

• 0.05

0.05 0.1 0.15 v(k)

Automation Lab IIT Bombay

System Identification 16 4/10/2012 System Identification 16

Unmeasured Disturbance Modeling

Measured output y(k) contains contributions from

Measurement errors (noise) Unmeasured disturbances

In additions modeling (approximation) errors arise while developing linear perturbation models If process is open loop unstable, OE model cannot be estimated using the prediction error approach

Thus, in order to extract true model parameters from the data, we need to carry out modeling of unmeasured disturbances (or noise) Noise is modeled as a stochastic process (sequence of random variables, which are typically correlated in time)

Automation Lab IIT Bombay

System Identification 17

Random Variable

Ref.: http://www.rslabntu.net/Random_Processes/Chapter1.pdf

A random variable is a mapping function which assigns

• utcomes of a random experiment to real numbers

Occurrence of the outcome follows certain probability

• distribution. Therefore, a random variable is completely

characterized by its probability density function (PDF).

Automation Lab IIT Bombay

System Identification 18 4/10/2012 System Identification 18

Stochastic Process

variable. random a is ,.) v(k , k k fixed a For process. stochastic the

• f

n' realizatio ' called function time

• rdinary

an is ) v(., function the , fixed a For = = ξ ξ ξ

.....} 1, 0, 1,

• ...

k : {v(k) variables stochastic

• f

family a as regarded be can process stochastic discrete A =

). v(k, variables two

• f

function as considered be may process random A ξ systems. control

• f

context in instants sampling represents k Index

Foundations of theory laid by Kolmogorov around 1930

Automation Lab IIT Bombay

Stochastic Process

System Identification 19 4/10/2012

Realizations of a stochastic process x

Automation Lab IIT Bombay

System Identification 20 4/10/2012 System Identification 20

Stochastic Process

Ref.: http://www.rslabntu.net/Random_Processes/Chapter1.pdf

Automation Lab IIT Bombay

System Identification 21 4/10/2012 System Identification 21

Stochastic Process

Examples of processes modeled as stochastic time series 1. Stock market and exchange rate fluctuations 2. Signals such as speech, audio and video, 3. Medical data such as a patient's EKG, EEG, blood pressure or temperature 4. Random movement such as Brownian motion or random walks. Heater Mixer Experimental Setup: Cold Water Temperature

Automation Lab IIT Bombay

Example: Global Annual Mean Surface Air Temperature Change

System Identification 22

Line plot of global mean land-ocean temperature index, 1880 to present, with the base period 1951-1980. The dotted black line is the annual mean and the solid red line is the five- year moving

• average. The

green bars show uncertainty estimates. Ref.: http://data.giss.na sa.gov/gistemp/gra phs_v3/

Automation Lab IIT Bombay

Example: Speech Recording

System Identification 23

200 400 600 800 1000 500 1000 1500 2000 2500 3000 3500 4000 4500 Sample (k) Amplitude Speech Recording of the syllable aaa...hhh

Sampled at 10000 Points per second with N =1020 Analysis of this signal is

• f interest

in computer recording of speech

Automation Lab IIT Bombay

System Identification 24 4/10/2012 System Identification 24

Discrete Time White Noise

Consider a random process {e(k): k= 1, 2, …} such that e(k) are independent and equally distributed random variables

2

(k)} variance{e mean{e(k)} that such variables random ed uncorrelat

• f

Collection σ = =

Such sequence is called discrete time white noise Simplest Random Process

variables random normal

• r

Gaussian (iid) d distrubute y identicall and nt independed are e(k) Noise White Gaussian

Automation Lab IIT Bombay

Example: Gaussian White Noise

System Identification 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

• 2
• 1

1 2 3 Time instant (k) Value Three realizations of zero mean white noise process with std = 1

Automation Lab IIT Bombay

Moving Average Process

System Identification 26

Consider a Gaussian white noise process {e(k): k= 1, 2, …} We can introduce a moving average that smoothes the series v(k) = (1/3) { e(k-1) + e(k) + e(k+1) } {v(k): k= 1, 2, …)} called Moving Average Process We can introduce an auto-regressive series x(k) = x(k-1) – 0.9 x(k-2) + e(k) {x(k): k= 1, 2, …)} called Auto Regressive Process Capable of capturing Quasi-Periodic Behavior in Data

Automation Lab IIT Bombay

Example: Moving Average Process

System Identification 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

• 4
• 3
• 2
• 1

1 2 Time instant (k) e(k) and v(k) White noise e(k) and Moving average process v(k) White Noise e(k) Moving average v(k)

Automation Lab IIT Bombay

Example: Auto-Regressive Process

System Identification 28

5 10 15 20

• 3
• 2
• 1

1 2 3 4 Time instant (k) e(k) and x(k) White Noise e(k) and Auto Regressive process x(k) White Noise e(k) Auto Regressive Process x(k)

Automation Lab IIT Bombay

PDF of Stochastic Processes

4/10/2012 System Identification 29

Let {v(k) : k=0,1,2,…} be a stochastic process

variables random l domensiona

• n

are points time distinct n at process random a

• f

values The

{ }

process random the

• f

function

• n

distributi l dimensiona finite the called is ies probabilit denotes P where ) ( ,....., ) v(k P ) ,..k k , k ; ,..., F( function The

1 n 2 1 n 1 n n

k v ξ ξ ξ ξ ≤ ≤ =

1

Automation Lab IIT Bombay

Gaussian or Normal Process

4/10/2012 System Identification 30

subsets possible all for and n all for Gaussian jointly is points time finite at ) ( )..... ( ), v(k subset any

• f
• n

distributi joint Then,

1 n

k v k v

2

Each variable v(k) has a Gaussian distribution for k=0,1,2,…

{ } ( )(

) { }

j j i i j i

k v k v k k v μ μ μ − − = = ) ( ) ( ) , ) ( E r(k variances its and E mean its by zed characteri completely is process Gaussian A

i i

Automation Lab IIT Bombay

Mean of a Stochastic Process

4/10/2012 System Identification 31

{ }

process stochastic the

• f

mean varying Time : (k) space) (sample ..... 2, , 0,1 k for ) , ( ) ( ) ( μ ξ ξ ξ μ Ω ∈ = = =

∫

Ω

k dF k k v E

Expected value

Consider discrete time stochastic process Y(k) = cos(2πk/25) + e(k) Example 1: Signal plus noise Let {e(k)} be a discrete time zero mean white noise process with standard deviation σ

Automation Lab IIT Bombay

Example: Mean

System Identification 32

E{y(k)} = E{ cos(2πk/25) } + E{ e(k) } = cos(2πk/25) Example 2: Moving Average process v(k) = (1/3) { e(k-1) + e(k) + e(k+1) } E{v(k)) = (1/3) E{ e(k-1) } + E {e(k) }+ E {e(k+1) } = 0 Example 3: Auto Regressive process x(k) = x(k-1) – 0.9 x(k-2) + e(k) E{x(k)} = E{x(k-1)} – 0.9 E{x(k-2)} + E{e(k)} Let μ = E{x(k)} = E{x(k-1)} = E{x(k-2)}

• 0.9 μ = E{e(k)} = 0 μ = 0

Automation Lab IIT Bombay

Stochastic Processes

4/10/2012 System Identification 33

Consider the stochastic process { v(k) } Do random two random variable at two different time instants, say v(k) and v(t), have any relation? Second order moments of the stochastic processes can be used to answer this question. Let {w(k) : k=0,1,2,…} be another stochastic process. Do random two random variable at two different time instants, say v(k) and w(t), have any relation?

Automation Lab IIT Bombay

Second Order Moments

System Identification 34

[ ] [ ][ ]

{ }

[ ][ ]

) , ; , ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ), ( ) , ( is process stochastic a

• f

covariance

• Auto

t s dF k k t t v k k v E t v k v Cov t k r

v 2 1 2 1

ξ ξ μ ξ μ ξ μ μ

∫∫

− − = − − = =

4/10/2012 System Identification 34

v(t) varoiable random with v(k) variable random

• f

dependence quantifies function covariance

• Auto

[ ] [ ]

{ }

2

) ( ) ( ) ( ), ( ) , ( k k v E k v k v Cov k k r

v

μ − = = = variance to reduces covariance auto t, k For

Automation Lab IIT Bombay

System Identification 35 4/10/2012 System Identification 35

Example: Discrete Time White Noise

Consider discrete time zero mean white noise process {e(k): k= 1, 2, …} with standard deviation σ Since e(k) and e(j) are uncorrelated when j ‡ k, it follows that

⎩ ⎨ ⎧ ± ± = = = ,.... , ) ( 2 1

2

τ τ σ τ for for r

e

[ ]

k j when ≠ = 0 ) ( ) ( j e k e E

Thus, auto-covariance of white noise process is

Automation Lab IIT Bombay

Example: MA process

System Identification 36

Consider the Moving Average Process v(k) = (1/3) { e(k-1) + e(k) + e(k+1) }

[ ] [ ][ ] { } [ ][ ] { }

) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ), ( ) , ( 1 1 1 1 9 1 + + + − + + + − = − − = = t e t e t e k e k e k e E t v k v E t v k v Cov t k r

v

[ ][ ] { }

{ } { } { } [ ]

2 2 2 2

3 1 1 1 9 1 1 1 1 1 9 1 σ = + + + − = + + + − + + + − = ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , ( k e E k e E k e E k e k e k e k e k e k e E k k r

vv

For k = t

Automation Lab IIT Bombay

Example: MA process

System Identification 37

[ ][ ] { }

{ } { } [ ]

2 2 2

9 2 1 9 1 1 1 2 1 9 1 1 σ = + + = + + + − + + + + = + ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , ( k e E k e E k e k e k e k e k e k e E k k r

vv

⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ≥ = = = = 3 2 9 1 9 2 3

2 2 2

t

• k

t

• k

/ t

• k

/ t k / ) , ( Summary σ σ σ t k r

vv

Other auto-covariance coefficients can be derived in a similar manner Note: Auto-covariance For the MA process depends only on the time separation (k-t) or lag and not n the absolute location of time points in series

Automation Lab IIT Bombay

Auto-correlation function

System Identification 38

) , ( ) , ( ) , ( ) , ( t t r k k r t k r t k

v v v v

= ρ

A normalized measure of association with values between ± 1 Measures predictability of the series at time t > k, i.e. v(t), using only the value v(k)

1 1 ≤ ≤ − ) , ( that follows it , inequality Schawtz

• Cauchy

From

v

t k ρ

Automation Lab IIT Bombay

Interpretation of Correlation Function

System Identification 39 4/10/2012

n correlatio strong means

• ne

to close

• f

Values

v

) , ( τ ρ + k k n correlatio NO means zero to equal

• f

Value

v

) , ( τ ρ + k k n correlatio negative means

• f

Values

v

< + ) , ( τ ρ k k

Investigation of shape of correlation function indicates temporal dependencies of the process

then ) ( v(t) relation linear a through v(t) predict can we If k v

1

β β + = ⎪ ⎩ ⎪ ⎨ ⎧ < − = > = when 1 when when 1

1 1 1

β β β ρ ) , ( t k

v

Automation Lab IIT Bombay

System Identification 40 4/10/2012 System Identification 40

Stationary Stochastic Process

{ }

v(k) E constant) (i.e. time

• f

t independen is Mean Process Random Stationary Weakly

v

μ = An ideal stochastic process often used in time series Modeling Examples: White noise, Auto Regressive Process, Moving average process

) ( ) ( ) , ( τ τ

v v v

r t k r t k r = − = = t)

• (k

difference

• f

function solely is time in dependence function ation Autocorrel

Automation Lab IIT Bombay

Cross-Covariance Function

System Identification 41

{ } { } [ ][ ] { }

) ( ) ( ) ( ) ( ) , ( as defined is w(k) and v(k) processes stochastic

• f

function covariance

• Cross

k t w k k v E t k r

w v vw

μ μ − − =

4/10/2012 System Identification 41

Measure of predictability of the series w(t) from the series v(k)

functions ation autocorrel are ) , ( and ) , ( where ) , ( ) , ( ) , ( ) , ( Function n Correlatio Cross t t r k k r t t r k k r t k r t k

w v w v vw v

= ρ

Automation Lab IIT Bombay

Example: AR Process

System Identification 42

Consider Auto Regressive (AR) series x(k) = x(k-1) – 0.9 x(k-2) + e(k)

[ ] [ ][ ] { } [ ][ ] { }

) ( ) ( ) ( . ) ( ) ( ) ( ) ( ), ( ) , ( k e k e k x k x E k e k x E k e k x Cov k k r

xe

+ − − − = − − = = 2 9 1

[ ][ ] { } [ ][ ] { } [ ][ ] { }

2

2 1 σ = = = − = − ) ( ) ( ) , ( ) ( ) ( ) ( ) ( Since k e k e E k k r k e k x E k e k x E

xe

[ ] [ ][ ] { } [ ][ ] { }

2

• 0.9

σ = − − − = − + − − − = − = − ) ( ) ( . ) ( ) ( ) ( . ) ( ) ( ), ( ) , ( 2 2 9 2 2 9 1 2 2 k e k x E k e k e k x k x E k e k x Cov k k r

xe

Automation Lab IIT Bombay

System Identification 43 4/10/2012 System Identification 43

Computing Statistics

{ } ∫

Ω

= = ⇒ = ≡ ) ( ) ( ) ( ) ( ) ( ) ( ) ( Function Density y Probabilit Cumulative ) ( ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ dP f g E d p dP f d dP P

{ } ∫

Ω

= Ω ∈ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ d f g g E g ) ( ) ( ) ( ), g(

• f

value expected Then,

• f

function continuos some be ) ( and ) f( p.d.f. with variable random a be Let

{ }

∑ ∫ ∑

= Ω =

≈ ≈ =

N i i N i i i N

g N P g dP g g E

1 1 2

1 ) ( ) ( ) ( ) ( ) ( ) ( ) f( from drawn samples (i.i.d) d Distribute y Identicall and t Independen represent ,...... , Let

1

ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ

Automation Lab IIT Bombay

System Identification 44 4/10/2012 System Identification 44

Estimating Auto-covariance From Data [ ][ ]

∑

− =

− + − =

τ

μ τ μ τ

N k v v v

k v k v N r

1

1 ˆ ) ( ˆ ) ( ) ( ˆ as estimated be can covariance

• auto

Sample

∑

=

=

N k

k v N

1 v

) ( 1 ˆ as estimated be can mean Sample time in signal

• f

values using estimated be can statistics sample process, random stationary a For μ

[ ][ ]

∑

− =

− + − =

τ

μ τ μ τ

N k w v vw

k w k v N r

1

1 ˆ ) ( ˆ ) ( ) ( ˆ as estimated be can covariance

• cross

Sample

Automation Lab IIT Bombay

Sample Correlation Functions

System Identification 45 4/10/2012 System Identification 45

) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ

v v v v v v

r r r r r τ τ τ ρ = =

[ ] [ ]

N E

v v

1 = ≠ = ) ( ˆ and all for ) ( ˆ then process noise white a is v(k) If v(k) process stationary a for ACF

• f
• n

distributi sample Large τ ρ σ τ τ ρ

When we estimate sample autocorrelation for a Gaussian white noise process, it turns out that

when ) ( ˆ > ≠ τ τ ρ

v

neglected? be to enough" small " is ) , ( ˆ whether judge we do How t k

v

ρ

Automation Lab IIT Bombay

Sample Correlation Functions

System Identification 46 4/10/2012 System Identification 46

0.95

• r

0.99

• f

value Typical ) P( with deviation normal standard

• f

value denotes where interval the

• utside

is peak

• bserved

the whether g determinin by t significan are ) ( ˆ in peaks whether assessing

• f

method rough a

• btain

we result, above

• n

Based Interval Confidence

/ / /

= = > ≠ α α ζ ζ ζ ζ ζ τ ρ

α α α 2 2 2

N

zero" to close "

• r

ant insignific considered be can . ) ( ˆ . i.e. 200 2.5758 ) ( ˆ 200 2.5758

• )

( ˆ

• f

values 0.99, and 200 N For : Example

/2 /2

• 1821

1821 < < − = < < = = = τ ρ ζ τ ρ ζ τ ρ α

α α

Automation Lab IIT Bombay

Discrete Time White Noise

System Identification 47

Autocorrelation function of 200 samples of a zero mean Gaussian White Process with σ =1

5 10 15 20

• 0.4
• 0.2

0.2 0.4 0.6 0.8 1 1.2 Autocorrelation of 200 samples of Gaussian White Process lags 99 % Confidence Interval (i.e. alfa = 0.99)

Automation Lab IIT Bombay

System Identification 48 4/10/2012 System Identification 48

Example: White Noise

Mean = 27.33 0C Variance = 0.633 0C

20 40 60 80 100 26 26.5 27 27.5 28 28.5 29

Sampling Instant Temperature (deg C)

Measurement Mean value

Experimental Data Temperature sensor kept in a beaker containing water at the room temperature

Automation Lab IIT Bombay

System Identification 49 4/10/2012 System Identification 49

Measurement Errors: Histogram

• 1.5
• 1
• 0.5

0.5 1 1.5 5 10 15 20 25 30 35 40 45

Measurement Error

• No. of Samples

Histogram of Measurement Errors

Experimental Data Estimated e(k)=T(k)-Tmean

Is this a white noise?

Automation Lab IIT Bombay

Estimated Autocorrelation Function

System Identification 50

Measurement errors can be modeled as a Gaussian white noise!

5 10 15 20

• 0.2

0.2 0.4 0.6 0.8 1 1.2 Autocorrelation Function: Temperature Data lags

Automation Lab IIT Bombay

Measurement Errors: Estimated pdf

System Identification 51

• 1.5
• 1
• 0.5

0.5 1 1.5 0.1 0.2 0.3 0.4 0.5 0.6 0.7

Deviation Temperature

Probability Density

Gaussian Density Function estimated fromTemperature data Estimated mean = -2.3e-15 Estimated std. dev. = 0.633

Automation Lab IIT Bombay

Example: Autocorrelation

System Identification 52

Strong positive Correlation with temperature deviations of previous 30 years

10 20 30 40 50

• 0.4
• 0.2

0.2 0.4 0.6 0.8 1 1.2 Autocorrelation Function : Yearly average global temperature deviations lags

Automation Lab IIT Bombay

Example: Speech Data

System Identification 53

50 100 150 200 250

• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 1 Autocorrelation Function of Speech Data lags

Original Signal: Sequence of repeating short signals ACF confirms This behavior

Automation Lab IIT Bombay

Sample Correlation Functions

System Identification 54 4/10/2012 System Identification 54

) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ

w v vw vw

r r r τ τ ρ =

[ ]

N

vw

1 = ) ( ˆ and mean zero with normal is ) ( ˆ

• f
• n

distributi sample large then white are w(k) and v(k) processes stationary If : fact following

• n

based intervals confidence Compute

vw

τ ρ σ τ ρ

1 Note ≤ ≤ − ) ( ˆ τ ρvw 1

neglected? be to enough" small " is ) ( ˆ whether judge we do How τ ρvw

Automation Lab IIT Bombay

El Nino and Fish Population

SOI: measure of changes in air pressure,

related to sea surface temperatures in central pacific (Central pacific region warms every 3 to 7 years due to the El Nino effect)

Both series exhibit repetitive behavior, with

regularly repeating cycles

Does fish population some how depend on

variation of SOI?

Use Cross correlation analysis to uncover

possible relation

System Identification 55

Automation Lab IIT Bombay

El Nino and Fish Population

System Identification 56

1950 1955 1960 1965 1970 1975 1980 1985

• 1

1 Southern Oscillation Index SOI 1950 1955 1960 1965 1970 1975 1980 1985 50 100 Estimated (Scaled) Fish Polulation Data Years

Automation Lab IIT Bombay

El Nino and Fish Population

System Identification 57

Strong positive correlation for values separated by 12 months Observations separated by 6 months are negatively correlated shows +ve excursions are related to –ve excursions 6 months removed.

Automation Lab IIT Bombay

El Nino and Fish Population

System Identification 58

• 50

50

• 0.8
• 0.6
• 0.4
• 0.2

0.2 0.4 Estimated Cross Correlation bet. SOI and Fish Polulation Lags

Cross correlation Peaks at lag of -6 Showing that SOI Measured at k-6 is associated with Fist population at time k

Automation Lab IIT Bombay

System Identification 59 4/10/2012 System Identification 59

OE Model for Two Tank System

50 100 150 200

• 0.5

0.5

y(k) OE(2,2,2): Measured and Simulated Outputs

50 100 150 200

• 0.05

0.05 0.1 0.15 v(k)

Automation Lab IIT Bombay

System Identification 60 4/10/2012 System Identification 60

OE Model : Correlation Functions

5 10 15 20 25

• 0.5

0.5 1 1.5 Output Error Model OE(2,2,2): Correlation function of residuals v(k) lag

• 30
• 20
• 10

10 20 30

• 0.4
• 0.2

0.2 0.4 Cross corr. function between input u(k) and residuals v(k) lag

Is the residual sequence auto- correlated? Is the residual sequence correlated with manipulated input?

Automation Lab IIT Bombay

System Identification 61 4/10/2012 System Identification 61

Signal Spectrum or spectral density

{ }

) ( frequency

• f

function real a always is ) ( exists. sum infinite the provided ) ( ) ( v(k) signal

• f

spectrum power define We ω ω τ ω

ωτ τ s j v v

e r Φ = Φ

− ∞ −∞ =

∑

{ } { }

) ( ) ( ) ( have we ) ( ) ( Defining τ τ − = ∞ → =

∑

=

k v k v E r k v N N Lim k v E

v N k 1

1

Note: Spectrum of signal v(k) represents Fourier transform of auto-covariance function

∫

−

=

π π ω

ω ω φ d e k r

v ik v

) ( ) (

Inverse Transform:

Automation Lab IIT Bombay

System Identification 62 4/10/2012 System Identification 62

Power Spectrum Density

Physical Interpretation

) , ( band in signal

• f

power Represents ) (

1 2

2 1

ω ω ω ω φ

ω ω

∫

d

Area under spectral density band thus represents the signal power in certain frequency band Total area under the curve is proportional to Variance of the signal

Automation Lab IIT Bombay

White Noise Process

System Identification 63

Power spectrum density of white noise π σ ω φ 2

2

= ) ( Power spectrum density of white noise is constant. Analogy with spectral properties of white light explains name given to the process All colored stochastic signals are generated by filtering white noise through a transfer function

Automation Lab IIT Bombay

System Identification 64 4/10/2012 System Identification 64

Spectrum of White Noise

0.2 0.4 0.6 0.8 1 10

• 0 .3

10

• 0 .2

10

• 0 .1

10 10

0 .1

P xx - X P ower S pectral D ensity F requency E stim ated

Theoretical

1

2 =

σ with noise white Gaussian mean zero

• f

Spetrum

White Noise: Uniform power spectrum density at all frequencies

) / ( frequency Nyquist where ) / ( Frequency Normalized T

N N

π ω ω ω = ≡ ≡

Automation Lab IIT Bombay

System Identification 65 4/10/2012 System Identification 65

Spectrum: Colored Noise

0.2 0.4 0.6 0.8 1 1 2 3 4 5 Pxx - X Power Spectral Density Frequency Theoretical Estimated from data

) ( ) ( . ) ( k e k v k v + − = 1 5 Process Regressive Auto

1

2 =

σ with mean zero is )} ( { : k e Case

Automation Lab IIT Bombay

Power Spectrum: Speech Signal

System Identification 66

0.2 0.4 0.6 0.8 1 10

1

10

2

10

3

10

4

10

5

10

6

10

7

Speech Signal: Power Spectral Density Normalized Frequency

This is certainly a colored noise and spectrum looks qualitatively similar to that of

• f AR process!

So, is there some AR process that governs this behavior and can we find it using data?

Automation Lab IIT Bombay

System Identification 67 4/10/2012 System Identification 67

Noise Modeling

Information about unmeasured disturbances in the past is contained in the output measurement record indirectly through their influence on past outputs ) ( ) ( ) ( ) ( k v k u q G k y + = Deterministic component Residue: unmeasured disturbances + measurement noise

[ ]

) ( p)

• y(k

....., 1),

• y(k

m),

• u(k

1),...,

• u(k

) ( k e f k y + =

Thus, an obvious choice of model structure is

Automation Lab IIT Bombay

System Identification 68 4/10/2012 System Identification 68

Equation Error Model

A discrete linear model, which captures the effect of past unmeasured disturbances, can be proposed as follows

time dead delay / Time : d ) ( ) ( ... ) 1 ( ) ( ... ) 1 ( ) (

1 1

k e m d k u b d k u b n k y a k y a k y

m n

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

How many past outputs do we include in the model?

• We can choose n such that error sequence {e(k)} is a

white noise sequence i.e. it contains no information about the disturbances in the past

Automation Lab IIT Bombay

System Identification 69 4/10/2012 System Identification 69

ARX Model Development

) ( ) 3 ( ) 2 ( ) 2 ( ) 1 ( ) ( ˆ .......... ) 4 ( ) 1 ( ) 2 ( ) 2 ( ) 3 ( ) 4 ( ˆ

2 1 2 1 2 1 2 1

N e N u b N u b N y a N y a N y e u b u b y a y a y + − + − + − − − − = + + + − − =

) ( ) 3 ( ) 2 ( ) 2 ( ) 1 ( ) ( 1 d with model ARX

• rder

nd 2' consider Thus,

2 1 2 1

k e k u b k u b k y a k y a k y + − + − + − − − − = =

Consider data obtained from two tank system and let us try developing an ARX model with n=2

) 3 ( ) ( ) 1 ( ) 1 ( ) 2 ( ) 3 ( have we equation, model the

• f

use recursive With

2 1 2 1

e u b u b y a y a y + + + − − =

Automation Lab IIT Bombay

System Identification 70 4/10/2012 System Identification 70

ARX : Parameter Identification

{

4 3 4 2 1 4 4 4 4 4 4 4 4 4 3 4 4 4 4 4 4 4 4 4 2 1 4 3 4 2 1 e θ Y ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ Ω ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − − − − − − − − = ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ e(N) ... ) e( ) e( b b a a N u N u N y N y u u y y u u y y y(N) ... ) y( ) y( 4 3 ) 3 ( ) 2 ( ) 2 ( ) 1 ( .. .. .. .. ) 1 ( ) 2 ( ) 2 ( ) 3 ( ) ( ) 1 ( ) 1 ( ) 2 ( 4 3

2 1 2 1

form matrix in Arranging

“Linear in Parameter” Model Advantage: Least square parameter estimation problem can be solved analytically

e θ Y + Ω =

Automation Lab IIT Bombay

System Identification 71 4/10/2012 System Identification 71

ARX : Parameter Identification

[ ] [ ]

θ θ Y θ Y e e = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = ∂ ∂ Ω − Ω − = =

T T T

b b a a

2 1 2 1

• ptimality

for condition Necessary : Function Objective ϕ ϕ ϕ ϕ ϕ ϕ

[ ]

e e θ θ

T

min ) ( ) , , , ( min ˆ estimation parameter square Least

2 2 2 1 2 1

= =

∑

= N k

k e b b a a

[ ] [ ] [ ]

θ Y θ θ Y θ Y θ = Ω − Ω − = ∂ Ω − Ω − ∂ = ∂ ∂

T T

ϕ

Automation Lab IIT Bombay

System Identification 72 4/10/2012 System Identification 72

ARX : Parameter Identification

[ ]

Y θ

T T LS

Ω Ω Ω = ⇒

−1

ˆ

definite positive is matrix Hessian minimum be to

• ptimum

for condition ufficient

2 2

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ θ ϕ S

ϕ ϕ

• f

minimum global be to happens it fact, In minimum a is ˆ definite ve always is which

2 2 LS T

θ θ ⇒ + Ω Ω = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂

q 0.01191 q 0.002506 B(q) q 0.5876 q 1.558

• 1

) A(q e(k) )u(k) B(q )y(k) A(q

3

• 2
• 2
• 1
• 1
• 1
• 1

+ = + = + =

Least square Estimates for Two tank system

Automation Lab IIT Bombay

System Identification 73

ARX Model: Output and Residuals

50 100 150 200 250

• 1
• 0.5

0.5 1 y(k)

ARX Model ARX(2,2,2): Comparison of Predicted and Measured Outputs

50 100 150 200 250

• 0.1
• 0.05

0.05 0.1 e(k) Time (samples)