### Create a StudySoup account

#### Be part of our community, it's free to join!

Already have a StudySoup account? Login here

# 3RD WORLD POLITICS POL S 147

UCSB

GPA 3.65

### View Full Document

## 34

## 0

## Popular in Course

## Popular in Political Science

This 48 page Class Notes was uploaded by Johnathan McKenzie on Thursday October 22, 2015. The Class Notes belongs to POL S 147 at University of California Santa Barbara taught by Staff in Fall. Since its upload, it has received 34 views. For similar materials see /class/227000/pol-s-147-university-of-california-santa-barbara in Political Science at University of California Santa Barbara.

## Reviews for 3RD WORLD POLITICS

### What is Karma?

#### Karma is the currency of StudySoup.

#### You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!

Date Created: 10/22/15

UNDERGRADUATE LECTURE NOTES ON SYSTEM IDENTIFICATION Jo uo P Hespanha April 1 2007 Disclaimer This is a draft and probably contains several typos Comments and information about typos are welcome Please contact the author at he spanhaQ ece ucsb edul Copyright to Joao Hespanhal Please do not distribute this document Without the author s consentl Contents H O 00 g ComputerControlled Systems 3 111 Computer control i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 3 1 2 Continuoustime systems i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 4 11211 Steadystate response to sine waves i i i i i i i i i i i i i i i i i i i i 5 1 3 Discretetime systems i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 6 11311 Steadystate response to sine waves i i i i i i i i i i i i i i i i i i i i 7 1132 Step and impulse responses i i i i i i i i i i i i i i i i i i i i i i i i 8 1 4 Discrete vsi continuoustime transfer functions i i i i i i i i i i i i i i i i i 8 1 5 MATLAB hints i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 9 116 To probe further i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 11 Parametric identi cation using leastsquares 13 2 1 Parametric identi cation i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 13 2 2 Leastsquares line tting i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 14 23 Vector leastsquares i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 16 2 4 MATLAB hints i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 18 215 To probe further i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 18 216 Exercises i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 20 Parametric identi cation of an ARX model 23 3 1 ARX Model i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 23 32 Identi cation of an ARX model i i i i i i i i i i i i i i i i i i i i i i i i i i 24 3 3 Known parameters i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 26 3 4 MATLAB hints i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 26 315 Exercises i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 28 Practical considerations in parametric identi cation 29 4 1 Signal scaling i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 29 4 2 Choice of inputs i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 30 4 3 Choice of model order i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 31 4 4 Choice of sampling frequency i i i i i i i i i i i i i i i i i i i i i i i i i i i 32 415 Exercises i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 34 Joao P Hespanha 5 Nonparametric identi cation 37 Nonparametric methods i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 37 5 2 Timedomain identi cation i i i i i i i i i i i i i i i i i i i i i i i i i i i i 38 5 21 Impulse response method i i i i i i i i i i i i i i i i i i i i i i i i i 38 522 Step response i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 39 523 Other inputs i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 40 53 Frequency response identi cation i i i i i i i i i i i i i i i i i i i i i i i i i 40 5 31 Sinewave testing i i i i H 40 5 32 Correlation method 41 5 4 MATLAB hints i i i i i i i i i i i 42 55 To probe further i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 42 5 6 Exercises i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 43 Attention When a margin sidebar nishes With i i i 77 more information about the topic can be found at the end of the lecture in the To probe further77 sectioni Introduction to System Identi cation The goal of system identi cation is to utilize inputoutput experimental data to determine a system s mo eli For example one may apply to the system a speci c input u measure the observed discretetime output y an try to determine the system s transfer function cf Figure 1 Figure 1 System identi cation from inputoutput experimental data Prerequisites 1i Laplace transform continuoustime transfer functions stability and frequency re sponse brie y reviewed here i 2 ztransform discretetime transfer functions stability and impulsestepfrequency responses brie y reviewed here Computercontrol system design brie y reviewed here aw i Familiarity with basic vector and matrix operations 9quot Knowledge of MATLABSimulinki Jo o P Hespanha Lecture 1 ComputerControlled Systems Summary This lecture reviews basic concepts in computercontrolled systems 1 Inputoutput models for continuoustime systems time and Laplace domains 2 Inputoutput models for discretetime systems time and ztransform domains 3 Transfer functions for system with sample and hold Contents 11 Computer control i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 3 12 Continuoustime systems i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 4 13 Discretetime systems i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 6 14 Discrete vs continuoustime transfer functions i i i i i i i i i i i i i i i i i 8 15 MATLAB hints i i i i i i i i i i i i i i i i i i i i i i i i i i i i 9 16 To probe further i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 11 11 Computer control Figure 11 show the typical block diagram of a control system implemented using a digi tal computer Although the output yt of most physical systems vary continuously as a function of time it can only be measured at discrete time instants Typically it is sampled periodically as shown in the left plot of Figure 12 Moreover the control signal ut cannot be changed continuously and typically is held constant between sampling times as shown in the right plot of Figure 12 It is therefore sometimes convenient to regard the process to be controlled as a discretetime system whose inputs and outputs are the discretetime signals Mk Wst was ukTs k e 1 2 3 c c c 4 Joao P Hespanha 6 W W process l MW 1 samp e Figure 11 Computer control architecture mm W 3 T H 4 ud1 0105 Ts 2T53TS Ts 2T53TS Figure 12 Sampling left and holding right The identi cation techniques that we will study allow us to estimate the transfer function of this discretetime system We can then use standard techniques from digital control to approximately recover the underlying continuoustime transfer function This is discussed in Section 14 below 12 Continuoustime systems Figure 13 Continuoustime system Consider the continuoustime system in Figure 13 and assume that the input and output signals satisfy the following differential equation 2M 71yltrlgt z m oy amumgt am1um 1gt 01212 0111 you 11 De nition The We recall that given a signal 95t with Laplace transform Xs the Laplace transform of umlateml Laplace tmnsform of a signal zt is given by Xs ODo e Stztdt System Identi cation the th derivative of xt is given by s Xs 7 3471940 7 347254 0 7 7 Nshm In particular when x and all its derivatives are zero at time t 0 the Laplace transform of the th derivative avg t simpli es to s Xsi Taking Laplace transforms to both sides of L1 and assuming that both y and u are zero at time zero as well as all their derivatives we obtain S Ys n11snills 232Ys 13Ys 0Ys amsmUs am11sm 1Us a232Us a13Us a0Usi This leads to where 1753 YS HcSUS amsm amqsm l 0113 a0 3 nilsnil 513 0 12 is called the continuoustime transfer function of the systemi The roots of the denominator are called the poles of the system and the roots of the numerator are called the zeros of the systemi The system is BIBO stable if all poles have negative real part In this case for every input bounded signal ut the output yt remains boundedi 121 Steadystate response to sine waves Suppose that we apply a sinusoidal input with amplitude A and angular frequency w given by ut A coswt V1320 to the system 11 with transfer function 12 Assuming that the system is BlBO stable the corresponding output is of the form ya 9A cosw gtgt m steadyrstate V1320 6t V transient where the gain 9 and the phase lt1 are given by the following formulas 1 lHCOWla lt1 AHCOw and 6t is a transient signal that decays to zero at t A 00 The Bode plot of a transfer function 1753 depicts 13 14 Matlab hint 1 function with numerator and denominator speci ed by um Matlab hint 2 zpkzpk creates a continuousetime transfer function with zeros poles and gain speci ed by z p k Matlab hint 3 pole sys and zero sys compute the poles and zeros of a system respectively Sidebar 1 Since there are other notions of stability one should use the term BIBO stable to clarify that it refers to the property that Boundedelnputs lead to Boundedeoutputs Sidebar 2 The angular frequency o is measured in radian per second and is related to the regular frequency f by the formula o 27rf where f is given in Hertz or cycles per second Matlab hint 4 bode sys draws the Bode plot of the system sys Sidebar 3 Since the output cannot depend on future inputs we must have n 2 m for 15 to be physically realizable De nition The unilateral zimmform of a signal is given by Xz iz kz c 160 6 Joao P Hespanha 1 the norm of H5040 in decibels dB which is given by 2010g10 Hcjw and 2 the phase of H5040 both as a function of the angular frequency w in a logarithmic scale In View of 1 371i4 the Bode plot provides information about how much a sinusoidal signal is ampli ed and by how much its phase is changed 13 Discretetime systems We ucUcT a Mt WC yckT Hdz Figure 14 Discretetime system Consider now the discretetime system in Figure 13 and assume that the output can be computed recursively by the following difference equation yk n 7 n71yk n 7 1 7 7 zzk 2 7 1yk 1 7 oyU amu s m am71uk m 71 a2uk 2 a1uk 1 a0uki 15 To obtain an equation that mimics 11 we can rewrite 15 as M n nily s n 71 2yk 2 llk 1 ova amu s m am71uk m 71 a2uk 2 a1uk 1 a0uki 16 We recall that given a signal 950 with ztransform Xz the ztransform of 95 Z is given by 2 Xz 7 21950 7 2171950 7 7 295M 71 In particular when x is equal to zero before time k Z the ztransform of 95 2 simpli es to 2 Xzi Taking ztransforms to both sides of 16 and assuming that both y and u are zero before time n 2 m we obtain System Identi cation 2nYz n11zn 1Yz 222Yz 12Yz 0Yz amszz am11zm 1Uz a222Uz a12Uz a0Uzi This leads to where Y9 HdZUZa amzm aqum l 112 a0 Hdzi 2 n712n71 iz o 17 is called the discretetime transfer function of the systemi The roots of the denominator are called the poles of the system and the roots of the numerator are called the zeros of the systemi The system is BIBO stable if all poles have magnitude strictly smaller than one In this case for every input bounded signal Mk the output yk remains boundedi 131 Steadystate response to sine waves Suppose that we apply a sinusoidal input with amplitude A and discretetime angular fre quency Q E 0 7r g1ven by uUs A cosQk Vke012m to the system 11 with transfer function 12 Assuming that the system is BlBO stable the corresponding output is of the form 21 gAcos9k gt 602 a mi where the gain 9 and phase lt1 are given by steadyrstate 9 I lHdejnla V tran si ent Vke012m lt1 AHdem and 6t is a transient signal that decays to zero at k A 00 The Bode plot of a transfer function Ha z depicts 1 the norm of Hdejn in decibels dB which is given by 2Olog10 lHdejnl and 2 the phase of Hdejn as a function of the angular frequency 9 in a logarithmic scale In view of 1187119 the 18 19 Bode plot provides information about how much a sinusoidal signal is ampli ed magnitude and by how much its phase is changed Exercise 1 Suppose that you sample at lKHz a unitamplitude continuoustime sinusoid uct with frequency 10Hzi Write the expression for the corresponding discretetime sinusoid ud i Matlab hint 5 tf numdenTs creates a discreteetime transfer function with sampling time and numerator an denominator speci ed by mum Matlab hint 6 poles and gain speci ed by z p k Sidebar 4 Since there are other notions of stability one should use the term BIBO stable to clarify that it refers to the property that Boundedelnputs lead to Boundedeoutputs Sidebar 5 Discretetime angular frequencies only take values from D to 7r since for Q 7r we already get the fastest discretetime sig nal nk Acos7rk AH Attention For 1 argument to H is e and no i continuousetime systems This means that Hz will be evaluated over the unit circle instead of the imaginary axis Matlab hint 7 bode sys draws the Bode plot of the system sys 8 Joao P Hespanha 132 Step and impulse responses Suppose that we apply a discretetime impulse given by uk m 1 k 0 0 k y 0 to the system 11 With transfer function 12 The corresponding output Mk is called the impulse response and its ztransform is precisely the system s transfer function a H 0 Consider noW an input equal to a discretetime step With magnitude oz ie 0 klt0 oz 20 Mk 3k and let yk be the corresponding output Which is called the step response Since the the impulse 6k can be obtained from the step 3k by MPW Vk 012m the impulse response Mk Which is the output to 6k can be obtained from the output 21 t0 3 by WPW Vke012m This is a consequence of linearity 14 Discrete vs continuoustime transfer functions When the sampling time Ts is small one can easily go from continuous to discretetime transfer functions shown in Figure 15 To understand hoW this can be done consider a continuoustime integrator 39 u 110 With transfer function 7 WS 7 1 Hcs 7 m 7 111 Suppose that the signals y and u are sampled every Ts time units as shown in Figure 12 and that we de ne the discretetime signals Mk MTS was ukTs k e 1 2 3 i i i System Identi cation 9 we 7 3 wk Hdz Figure 15 Discrete vs continuoustime transfer functions Assuming that ut is approximately constant over an interval of length Ts we can take a nite differences approximation to the derivative in 110 which leads to W Ts W Ts In particular for t kTs we obtain ydk 1 ydk Ts udUs 112 Taking the ztransform we conclude that 2Ydz 7 Ydz 7 7 Ydz 7 Ts 7 Udz 42gt Ha z 7 12 7 2151 113 Comparing 111 with 113 we observe that we can go directly from a continuoustime transfer function to the discretetime one using the so called Euler transformations 2 7 1 Ts s gt 42gt 2 gt 1 3T5 A somewhat better approximation is obtained by taking the righthandside of 112 to be the average of ud s and ud s 1 This leads to the so called Tustin or bilinea r transformation 22 7 1 7 Z 2STS H H Tsz1 273T 15 MATLAB hints Matlab Hint 1 5 tf The command systftf numden assigns to sys tf a MAT LAB continuoustime transfer function mum is a vector with the coefficients of the numerator of the system s transfer function and den a vector with the coefficients of the denominator The last coefficient must always be the zeroorder one Eg to get 3373 you should use num 2 0 den1 O 3 Sidebar 6 In discretetime an integrator has a po e at z 1 whereas in continuous7ti1ne the pole is at s 0 Matlab hint 8 c2d and d2c convert transfer functions from continuous to discretetime and vice versa Sidebar 7 Why 10 Joao P Hespanha The command systftf numdenTs assigns to sys tf a MATLAB discretetime trans fer function with sampling time equal to Ts For transfer matrices mum and den are cell arrays Type help tf for examples Optionally one can specify the names of the inputs outputs and state to be used in subsequent plots as follows systftf numden InputName input1 input2 DutputName output1 output2 StateName input1 input2 The number of elements in the bracketed lists must match the number of inputsoutputs and state variables Matlab Hint 2 6 zpk The command systfzpkz p k assigns to sys tf a MATLAB continuoustime transfer function z is a vector with the zeros of the system p a vector with its poles and k the gain Eg to get MW you should use z0p1 3 k2 The command systfzpkz p k TS assigns to sys tf a MATLAB discretetime transfer function with sampling time equal to Ts For transfer matrices z and p are cell arrays and k a regular array Type help zpk for examp es Optionally one can specify the names of the inputs outputs and state to be used in subsequent plots as follows systfzpkzpk InputName input1 input2 DutputName output1 output2 StateName input1 input2 The number of elements in the bracketed lists must match the number of inputsoutputs and state variables D Matlab Hint 4 7 bode The command bodesys draws the Bode plot of the system sys To specify the system one can use 1 systf numden where mum is a vector with the coefficients of the numerator of the system s transfer function and den a vector with the coefficients of the denominator The last coefficient must always be the zeroorder one Eg to get 3 one should use mum 2 0 den1 O 3 2 syszpk z pk where z is a vector with the zeros of the system p a vector with its poles and k the gain Eg to get El n one should use z0p1 3 k2 3 sysssABCD where ABCD are a realization of the system B Matlab Hint 8 c2d and d2c The command c2dsyscTs tustin converts the continuoustime LTl model sysc to a discretetime model with sampling time Ts using the Tustin transformation To specify the continuoustime system sysc one can use System Identi cation 11 11 sysctf numden where mum is a vector with the coefficients of the numerator of the system s transfer function and den a vector with the coefficients of the denominator The last coefficient must always be the zeroorder one Big to get 337 you should use num 2 0 den1 O 3 1 sysczpkzpk where z is a vector with the zeros of the system p a vector with r r 2 1ts poles and k the ga1n1 Big to get W you should use z 0p 1 3 k 2 31 syscss ABCD where ABCD are a realization of the system The command d2c sysd tustin converts the discretetime LTl model sysd to a continuous time model using the Tustin transformation To specify the discretetime system sysd one can use 11 sysdtf numdenTs where Ts is the sampling time mum is a vector with the coef cients of the numerator of the system s transfer function and den a vector with the coefficients of the denominator The last coefficient must always be the zeroorder one Big to get 327 you should use num 2 0 den 1 O 3 1 sysdzpkzpkTs where Ts is the sampling time 2 is a vector with the zeros of the system p a vector with its poles and k the gain Big to get MW you should use z0p 1 3 k2 DJ 1 sysdss ABCDTS where Ts is the sampling time ABCD are a realization of the system B 16 To probe further Sidebar 7 Tustin Transformation Consider a continuoustime integrator in 1110 and assume that ut is approximately linear on the interval kTs k HTS iiei k 1 Ts it t7 k s ut udk Two 1t mu 1Ts1 Then NUTS k 1 Ts it t7 kTs ltudk in 1gtdt yum 1 ww Ts ICTs Ts Ts yak Twas udk 1 Taking the ztransform we conclude that 2Ydz 7 Ydz 1012 2Udz Ydz Ts1 2 H 7 7 Ts 2 a W UM 22 71 Comparing this with 1111 we observe that we can go directly from a continuoustime transfer function to the discretetime one using the so called Tustin or bilinea r transforma tion 22 7 1 2 3T5 H7 42gt 2H 1 273115 Tsz1 D Jo o P Hespanha Lecture 2 Parametric identi cation using leastsquares Summary This lecture presents the basic tools needed for parametric identi cation using the method of leastsquares 1i Formulation of the parametric transfer function identi cation problem 2 Line tting using the leastsquares method 3 Surface tting using the leastsquares method Contents 21 Parametric identi cation i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 13 2 2 Leastsquares line tting i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 14 23 Vector leastsquares i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 16 2 4 MATLAB hints i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 18 25 To probe further i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 18 20 2 6 Exercises i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 21 Parametric identi cation 1n parametric identi cation one attempts to determine a nite number of unknown param eters that characterize the model of the systemi The following is a typical problem in this class 14 Joao P Hespanha Problem 1 Parametric transfer function identi cation Determine the coefficients a and t of the system s rational transfer function HM 7 amzm amiizm l 112 a0 2 n712n71 iz 50 l The number of poles n and the number of zeros in are assumed known This can be done using the methods of leastsquares which we introduce next starting from a simple line tting problem and eventually building up to Problem 1 in Lecture 3 The methods of leastsquares is generally attributed to Gauss 3 but the American mathematician Adrain 1 seems to have arrived at the same results independently A detailed treatment of leastsquares estimation can be found eg in 5 Chapter 7 22 Leastsquares line tting Suppose that we have reasons to believe that two physical scalar variables x and y are approximately related by a linear equation of the form y ago b 21 but we do not know the value of the constants a and b The variables y and as could be eg a voltage and a current in a simple electric circuit or the height of a uid in a tank and the pressure at the bottom of the tank cf Figure 21 z current y voltage gt R i i i z height gt i i 3 pressure a Voltage vs current b Pressure vs height Figure 21 Linear models Numerical values for a and b can be determined by conducting several experiments and measuring the values obtained for x and y We will denote by 951 y the measurements obtained on the ith experiment of a total of N experiments As illustrated in Figure 22 it is unlikely that we have exactly ylaxlb Vi 12luN 22 because of measurement errors and also because often the model 21 is only approximate System Identi cation 15 a Measurement noise b Nonlinearity Figure 22 Main causes for discrepancy between experimental data and linear model Leastsquares identi cation attempts to nd values for a and b for which the left and the righthandsides of 22 di er by the smallest possible error More precisely the values of a and b leading to the smallest possible sum of squares for the errors over the N experiments This motivates the following problem Sidebar 8 Statistical interpretation 0 Problem 2 Basic LeastSquares Find values for a and b that minimize the following leasmquares39 quot Sumof Squares Error SSE N SSE Emmi b 7 31 11 Solution to Problem 2 This problem can be solved using standard calculus All we need to do is solve BSSE N all Z2xlaxlb7yl0 11 BSSE N T Z2axlb7yl0 11 for the unknowns a and b This is a simple task because it amounts to solving a system of linear equations N N N 2 Ema 2 2951 A 515115 N N A 29632 021951 ME z 2ltszgta2Nb2Zyz b 1n the above expressions we used f1 and g to denote the solution to the leastsquares mini mization These are called the leastsquares estimates of a and b respectively Sidebar 9 Why Notation 1 a 1 denotes the inner product of a and b Sidebar 10 The line tting problem is a special case of this one for z z 1 a 9 a 7 Attention We are now using to denote the ith experiment since the subscript in z39 is needed to denotes the jth entry of the vector 2 16 Joao P Hespanha Attention The above estimates for a and b are not valid when WE 96 i Z 962 0 z 1 because this would lead to a division by zero It turns out that 23 only holds when the as obtained in all experiments are exactly the same as shown in Figure 23 Such set of 23 Figure 23 Singularity in leastsquares line tting due to poor datal measurements does not provide enough information to estimate the parameter a and b that characterize the line de ned by 21 23 Vector leastsquares The leastsquare problem de ned in Section 22 can be extended to a more general linear model where y is a scalar 2 2 21 22 an mUector and these variables are related by a linear equation of the form 2m 24 where 9 2 91 92 to determine 19m is an mvector with parameters whose values we would like To determine the parameter vector 9 we conduct N experiments from which we obtain measurements i E 1 2 i i i N where each denotes one mvector and each a scalar The leastsquares problem is now formulated as follows Problem 3 Vector LeastSquares Find values for 9 91 92 irnize the following Sumof Squares Error SSE 19m that min N SSE 2 go 9 7 you Solution to Problem 5 This problem can also be solved using standard calculusi Essentially we have to solve BSSE 7 BSSE 7 7 BSSE 70 891 8192 quotquot 86m System Identi cation 17 for the unknown 91 The equation above can also be rewritten as Notation 2 szz denotes the gradient of 7 BSSE BSSE BSSE 7 z VesSE l 391 392 a em l 0 fl However to simplify the computation it is convenient to write the SSE in vector notation Sidebar 11 Vector De ning E to be an Nvector obtained by stacking all the errors 9 7 on top of Emma is convenient bOth or algebraic manipulations eaCh Other Lea and for e icient MATLAB computations 2 9 7 211 2 9 7 212 E 2 2W 9 MN NX1 We can write Notation 3 Given a matrix M we denote the SSE EE transpose obey Ml Moreover by viewing 9 as a column vector we can write E Z9 7 Y where Z denotes a N X m matrix obtained by stacking all the row vectors on top of each other and Y an mvector obtained by stacking all the on top of each other itet E13 it 2 y 2 Z Y Ngtltm Ngtlt1 Therefore SSE 9Z7YZ97Y9ZZ972YZ9YY 25 It is now straightforward to compute the gradient of the SSE and nd the solution to Sidebar 12 The gradient VesSE 0 of a quadratic function zQzczd is 7 given by VeSSE 29 Z Z 7 ZYZ 0 42gt 9 Y ZZ Z 1 26 WNW zQQgtc which yields the following leastsquares estimate for 9 Matlab him 9 ZY computes directly ZZ 1ZY A 7 in a very e 39icient way which 9 ZZ lzY 2 7 avoids inverting 22 by solving directly the linear Since Z Z E and Z Y E we can rewrite the above formula as system of equations in 26 This is desirable because N N this operation is often 9 Rilf R1 f 1117cond1tioned 11 11 Matlab hint 10 svdA provides a measure of how close A is to being singular 18 Joao P Hespanha One can gauge the quality of the t by computing the SSE achieved by the estimate ie the minimum achievable squarederror Replacing the estimate 27 in 25 we obtain SSE 7 miniH2 717 Y ZZ Z 1Z Y 7 7 Y Z 2 8 HYHQ HYHQ W W When this errortosignal ratio is much smaller than one the mismatch between the left and righthandsides of 24 has been made signi cantly smaller than the output D Attention The above estimate for 9 is not valid when the m X in matrix R is not invertible Singularity of R is a situation similar to 23 in line tting and it means that the are not suf ciently rich to estimate 9 In practice even when R is invertible poor estimates are obtained if R is close to a singu art D 24 MATLAB hints Matlab Hint 10 svd The command svdA computes the singular values of the matrix A The singular values of a matrix A are the squareroots of the eigenvalues of A A Singular values are always real and nonnegative 1 The largest singular value written omax provides a measure of how much a matrix can amplify a vector In particular HAxH amaxiAmzi F The largest singular value written 7min provides a measure of how little a matrix can amplify a vector In particular leH 2 amulAlHIll For singular matrices ominA 0 D5 i The ration between the largest and smallest singular values is called the condition number of a matrix UmaxA 0minA provides a measure of how close a matrix is to being singular It is equal to 0 for singular matrices lnverting matrices with large condition numbers is very difficult and often results in large numerical errors due to roundoff MA I 25 To probe further Notation 1 Inner product Given two mvectors a a1 a2 a1 112 am and a am a b denotes the inner product of a and b ie a b folly 11 System Identi cation 19 Notation 2 Gradient Given a scalar function ofn variables fx1 i i i 95m Vgcf denotes the gradient of f ie a a a Vscfxlax2 axm i amm D Sidebar 8 Statistical interpretation of leastsquares Suppose that the mismatch between left and the righthand sides of 22 is due to uncorrelated zeromean noise In particular that ylaxlbnl Vi l2iun where all the n are uncorrelated zeromean random variables with the same variance 72 ie Enl 0 72 Enl nj 0Vij y ii The GaussMarkov Theorem stated below justi es the wide use of leastsquares estimation Theorem 1 GaussMarkov The best linear unbiased estimator BLUE for the param eters a and b is the leastsquares estimator Some clari cation is needed to understand this theorem 1 The GaussMarkov Theorem 1 compares all linear estimators ie all estimators of t e orm ampa1y1anyn h 1yi nym where the 01 are coef cients that may depend on the 951 F It then asks what is the choice for the 01 t that lead to an estimator that is best in the sense that it is unbiased ie for which Em 7 1 BE 7 b and ii results in the smallest possible variance for the estimation error ie that minimizes Ema 7 12 8 7 b The conclusion is that the leastsquares estimator satis es these requirements Unbiasedness means that when we repeat the identi cation procedure many time although we may never get a a an the estimates obtained will be clustered around the true values Minimum variance means that the clusters so obtained will be as narrow as possible Sidebar 9 Singularity of line t The estimates for a and b are not valid when M as 7 Z a 7 0 20 Joao P Hespanha To understand when this can happen let us compute NEW 7 My N295 7 2NMle 7121 but N11 2195 so NEW 7 My N295 7 951 This means that we run into trouble when 2xz My 0 which can only occur when all the x1 are the same and equal to 111 D Sidebar 12 Gradient of a quadratic function Given a m X 711 matrix Q and a mrow vector c the gradient of the quadratic function W m m xQx cx Z qujxlxj Zeal 11 11 11 To determine the gradient of f we need to compute m m 811 m m 8951 m 8951 3995 7 007 c 89 2411qu 189mg 24qu J Back Z 189 However g2 and 3 are zero wheneverz39 y k and j y k respectively and 1 otherwise Therefore the summations above can be simpli ed to a I m m m 3 211111961 quxj Ck 201111 11111 Ck 11 11 11 Therefore mm 7 7355 735 7315 21110111 111961 C1 211 11m 111mm Cm xQQ c D 26 Exercises Exercise 2 Electrical circuit Consider the electrical circuit in Figure 214 where the voltage U across the source and the resistor R are unknown To determine the values of U and R you place several test resistors r1 between the terminals A B and measure the voltages V1 and currents 1 across themi System Identi cation 21 Figure 24 Electrical circuit ll Write a MATLAB script to compute the voltages V and currents I that would be obtained When U 10V R MS and the n are equally spaced in the interval 100910k li Add to the correct voltage and current values measurement noise Use for the currents and for the voltages Gaussian distributed noise With zero mean and standard deviation ilmA and lOmV respective y The script should take as input the number N of test resistorsi 2 Use leastsquares to determine the values of R and U from the measurements for N E 5 10 50 100 1000 Repeat each procedure 5 times and plot the average error that you obtained in your estimates versus Ni Use a logarithmic scale What do you conclude D Exercise 3 Nonlinear spring Consider a nonlinear spring with restoring force F 7z avg Where x denotes the spring displacementi Use leastsquares to determine linear models for the spring of the form F ax b Compute values for the parameters a and b for ll forces evenly distributed in the interval 71 ll 2 forces evenly distributed in the interval 75 i5 and 3 forces evenly distributed in the interval 71 1 For each case ll calculate the SSE 2i plot the actual and estimated forces vs at and 3 plot the estimation error vs 95 D Jo o P Hespanha Lecture 3 Parametric identi cation of an ARX model S ummary This lecture explains how the methods of leastsquares can be used to identify ARX models 1 AutoRegression model with eXogeneous inputs ARK 2 Identi cation of an ARX model 3 Partial identi cation known parameters Contents 31 ARX Model i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 23 32 Identi cation of an ARX model i i i i i i i i i i i i i i i i i i i i i i i i i i 24 3 3 Known parameters i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 26 3 4 MATLAB hints i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 26 3 5 Exercises i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 28 31 ARX Model Suppose that we want to determine the transfer function of the discretetime system in Figure Sill ln leastsquares identi cation one converts the problem of estimating into the vector leastsquares problem considered in Section 23 This is done using the ARX model that will be constructed below The ztransforms of the input and output of the system in Figure 31 are related by Yz amzmam71zm71a12a0 2 uilznil 12 0 23 24 Joao P Hespanha Figure 31 System identi cation from inputoutput experimental data Where the a and the t denote the coefficients of the numerator and denominator polyno n mials of Multiplying the numerator and denominator of by 2 we obtain the transfer function expressed in negative powers of 2 Yz 7 amz m am1z m 1 0112 1 0102 12 7 1 n71271 39 39 39 5124quot 502 and therefore Yz 1z 1Yz lzin1Yz ozi lz amzinerUQ am71zinm71Uz alzin1Uz aozi UQ Taking inverse ztransforms we obtain 21 n71yk 71 1yk n 1 oyk n amu s 7 n m am71uk 7 n m 71 a1uk 7 n1a0uk 7 31 This can be reWritten compactly as gas was6 32 Where the nm1vector 9 contains the coefficient of the transfer function and the vector k the past inputs and outputs ilel 67 n71 l o am am71 01 00 kyk71 yk7n uk7nm uk7n The vector k is called the regression vector and the equation 32 is called an ARX model a short form of AutoRegression model With eXogeneous inputsl 32 Identi cation of an ARX model We are noW ready to solve the system identi cation Problem 1 introduced in Lecture 2 by applying the leastsquares method to the ARX model Solution to Problem 1 1 Apply a probe input signal Mk k E 1 2 l l l N to the systeml 2 Measure the corresponding output yk k E 1 2 l l lNl System Identi cation 25 3 Determine the values for the parameter 9 that minimize the discrepancy between the left and the righthandsides of 32 in a leastsquares sense According to Section 23 the leastsquares estimate of 9 is given by where W1 W2 gt2 w N or equivalently R The quality of the t can be assessed by computing the errortosignal ratio When this quantity is much smaller than one the mismatch between the left and right handsides of 32 has been made signi cantly smaller than the output ltIgt ltIgt 1ltIgt Y y0 y17nu17nm ya min ult2inmgt M21 yN ngt uN4nmgt N R I lt1gt lt1gt ZWWWM k21 SSE 7 HWP 1 Y ltIgt YY Matlab hint 11 PHIY computes directly from the matrix PHI and the vector Y Y Attention When the system is at rest before 33 141 W k 1 uk 2116 0 39 2 70 2 We 3 0 Otherwise if the I 7 j past inputs are unknown uNLn MN one needs to remove from and Y all rows with unknown data Matlab hint 12 arx is very convenient to perform leastrsquares identi cation of ARX models because it does not require the need to explicitly construct the matrix and the vector N fYZy M k1 Matlab hint 1 3 compare D is useful to validate the quality of the t Attention One may combine several inputoutput pairs of signals to identify one ARX model by stacking all the data available in one matrix ltIgt and one vector Y Suppose for example that one wants to identify an ARX model with two poles and two zeros based on an inputoutput pair u1ky1k k 12 of length N1 and another input output pair WW 212k k 12 of length N2 One would construct 3112 3113 1 we 3123 mmiuwmim mM4mM4gt u N2 3111 13 12 3112 14 13 mm u1N71 3121 23 22 3122 24 23 Matlab hint 14 This arx command does not lend itself to this type of multieinput processing WM WNQ u11 3113 u12 3114 mmem mmgt u2lt1gt Y m3 u22 3124 mmim mm Sidebar 1 3 The transfer function Hz is proper only if Hz is strictly proper However this poses no di 39iculties for this method of system identi cation Sidebar 14 The new output Mk is not a causal function of the original output yk but this is of no consequence since we have the whole yk available we carry out identi cation 26 1050 P Hespanha Attention Two common causes for errors in leastsquares identi cations of ARX models are 1 incorrect construction of the matrix ltIgt andor vector Y 2 incorrect construction of the identi ed transfer function from the entries in the least squares estimate Both errors can be avoided using the MATLAB command arx 3 3 Known parameters Due to physical considerations one often knows one or more poleszeros of the process For examp e 1 one may know that the process has an integrator which corresponds to a continuous time pole at s 0 and consequently to a discretetime pole at 2 or vrr m that the process has a which cone pond to a quot V s 0 and consequently to a discretetime zero at 2 1 zero at In this case it suffices to identify the remaining poleszeros which can be done as follows Suppose that it is known that the transfer function has a pole at 2 ie Hz 1 H 2 Z 7 lt gt where corresponds to the unknown portion of the transfer function In this case the ztransforms of the input and output of the system are related by Yz 7 1 7 12 7 2 7 Hz and therefore m 7 Z W 7 H gt where 172 2 2 7 Yz 9 yk 1 7 yk This means that we can directly estimate by computing 9 prior to identi cation and then regarding this variable as the output instead of 34 MATLAB hints Matlab Hint 12 arx The command arx from the identi cation toolbox performs least squares identi cation of ARX models To use this command one must System Identi cation 1 Create a data object that encapsulates the inputoutput data using dataiddata y u Ts where u and y are vectors with the input and output data respectively and Ts is the sampling interval 2 Compute the estimated model using mode larx data nanb nkl where data is the object with inputoutput data na 11b nk are integers that de ne the degrees of the numerator and denominator of the transfer function see below and model a MATLAB object with the coefficients of the transfer function The coefficient are returned for negative powers of 2 In particular the estimated transfer function is given by Yz 7 2 r kmodelia1modelia22 1 modeliana 12 3 12 7 modelib1 modelib22 1 modelibnb 12 b The estimated discretetime transfer function can be recovered using the MATLAB com man sysd tfmodelMe asured The object model contains a lot more information about the estimated transfer function which includes 1 modelda and modeldb give the standard deviations associated with the estimates of the coef cient A large value for one of the standard deviations indicates that the matrix R ltIgtltIgt is close to singular and points to little con dence on the identi ed model 2 modelnoise provides a measure of the quality of t and is directly related to the However this measure is not normalized opposite to 28 You can type idprops idpoly at the MATLAB prompt to obtain more information about other information that is available in model The MATLAB command arx uses a more sophisticated algorithm than the one described in these notes so the results obtained using arx may not match exactly the ones obtained using the formula 33 Matlab Hint 13 compare The command compare from the identi cation toolbox allows one to compare experimental outputs with simulated values from an estimated model This is useful to validate model identi cation The command compare data model plots the measured output from the inputoutput data object data and the predictions obtained from the estimated model model See MATLAB hints 12 19 for details on these objects D Sidebar 15 MATLAB is an object oriented language and model a is the property a of the object model Sidebar 16 Typically one chooses nb equal to the expected system dimension nk equal to the expected input delay and uh 7 nk This corresponds to yk being a function of yk 7 1 yk 7 nb and of 1k 7 nk 1k 7 nb 28 Joao P Hespanha 35 Exercises Exercise 4 Known zero Suppose that the process is known to have a zero at z 4 ie that H2WHL where corresponds to the unknown portion of the transfer function How would you estimate What if you known that the process has both a zero at y and a pole at X D Exercise 5 Selected parameters The data provided was obtained from a system with transfer function 275 m eemeew where p is unknown Use the leastsquares method to determine the value of p D Lecture 4 Practical considerations in parametric identi cation Summary This lecture discusses several practical considerations that are crucial for successful identi cationsl 1 Signal scaling 2 Choice of inputs 3 Choice of model order 4 Choice of sampling frequency Contents 41 Signal scaling l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l 29 4 2 Choice of inputs l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l 30 4 3 Choice of model order l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l 31 4 4 Choice of sampling frequency 32 34 4 5 Exercises l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l 41 Signal scaling For the computation of the leastsquares estimate to be numerically well conditioned it is Attention After important that both inputs and outputs have roughly the same order of magnitude Because i ffnti c tion one m then the units of these variable are generally quite different this often requires scaling of these a JHSt t e sysiem gain to re ect the scaling performed variables It is good practice to scale both inputs and outputs so that all va39rzable take on the signals used for normal values in the same range eg the interval 71 Ill identi cation 29 Sidebar 1 7 WhyJ H1z 2z21 and H2z W have the same value for z ejW l Therefore they will lead to the same yk when uk is a sinusoid with frequency 7r4 This input will not allow us to distinguish between H1 and H2 30 Joao P Hespanha 42 Choice of inputs The quality of leastsquares estimates depends signi cantly on the input Mk used The choice of inputs should take into account the following requirements 1 9quot The input should be su iciently rich so that the matrix R is nonsingular and is well conditioned Otherwise estimation is not possible a Single sinusoids should be avoided because a sinusoid of frequency 9 will not allow us to distinguish between different transfer functions with exactly the same value at the point 2 519 b Good inputs include squarewaves the sum of many sinusoids or a chirp signal ie a signal with timevarying frequency The resulting output yk should be much larger than the measurement noise In particular large inputs are needed when one wants to probe the system at frequencies for which the noise is large However when the process is nonlinear but is being approximated by a linear model large inputs may be problematic because the linearization may not be valid As shown in Figure 41 there is usually an optimal input level that needs to be determined by trialanderror SSE H y H nonlinear D0155 behavior dominates optimal input magnitude Figure 41 Optimal choice of input magnitude The input should be representative of the class of inputs that are expected to appear in the feedback loop a The input should have strong components in the range of frequencies at which the system will operate b The magnitude of the input should be representative of the magnitudes that are expected in the closed loop c When a stabilizing controller C is available data should be collected in closed loop with a rich reference signal In practice one then identi es the closedloop transfer function PC T1PC 4 1 System Identi cation 31 and then recovers the openloop transfer function P by solving 41 for the process T P 7 Cl 7 T Attention ln closedloop identi cation if the controller is very good the output will follow the reference very closely and the identi ed closedloop transfer function will be almost the identity over the range of frequencies of interest regardless of the process transfer function This will generally lead to very poor identi cation For identi cation purposes a sluggish controller is desirable Validation The choice of input is perhaps the most critical aspect in good system identi cation After any identi cation procedure the following steps should be followed to make sure that the data collected is adequate 1 Check if the matrix R is far from singular If not a different richer input should be used Check the quality of t by computing If this quantity is not small one of the F following is probably occurring a There is too much noise and the input magnitude needs to be increased b The inputs are outside the range for which the process is approximately linear and the input magnitude needs to be decreased c The assumed degrees for the numerator and denominator are incorrect see Sec tion 43 9quot After an input signal Mk passed the checks above repeat the identi cation experiment with the input auk with oz 71 oz 5 If the process is in the linear region the measured outputs should roughly by equal to Gagk and the two additional experiments should approximately result in the same transfer function When one obtains larger gains in the experiment with oz 5 this typically means that the process is saturating and one needs to decrease the magnitude of the input 43 Choice of model order A signi cant difficulty in parametric identi cation of ARX models is that to construct the regression vector k one needs to know the degree n of the denominator In fact an incorrect choice for n will generally lead to dif culties 1 Selecting a value for 71 too small will lead to mismatch between the measured data and the model and ngg will be large Selecting a value for 71 too large is called overparametem39zation and it generally leads to to B being close to singular To understand why suppose we have a transfer function 1 271 Sidebar 18 When n is chosen correctly m too large will generally not lead to trouble so one may as well chose the largest one that is physically meaningful If it is known that there is a delay of at least 1 one should make m n 7 cf equation 31 Joao P Hespanha but for estimation purposes we assumed that n 2 and therefore attempted to deter mine constants ozl such that 1222 112 04 H Z 22 iz 50 If the model was perfect it should be possible to match the data with any 01 such that a2220112010 2 0120011010 Pa 22 iz o Zilxzip i1p 5010 where p can be my number This means that the data is not sufficient to determine the values of the parameters 10 g l which translates into B being singular When there is noise it will never be possible to perfectly explain the data and the smallest SSE will always be positive either with n 1 or n 2 However in general different values ofp will result in different values for SSE In this case leastsquares estimation will produce the speci c value ofp that is better at explaining the noise77 which is not physically meaningful When one is uncertain about which values to choose for m and n the following procedure should be followed 1 Select for n the smallest possible integer that is physically meaningful conceivably 0 2 Do leastsquare estimation with m n a If the SSE is large and R is far from singular increase 71 and repeat the identi cation b If R is close to singular or the SSE did not decrease signi cantly then 71 got too large go back to the previous 71 and stop One needs to exercise judgment in deciding when the SSE is large77 or when R is farclose to singular77 Physical knowledge about the model should play a major role in deciding model orders Moreover one should always be very concerned about identifying noise as opposed to actually identifying the process model 4 4 Choice of sampling frequency For a discretetime model to accurately capture the process continuoustime dynamics the sampling frequency should be as large as possible However this typically leads to difficulties in system identi cation due to a couple of factors 1 If the poles of the continuoustime system are equal to 31 32 1 Sn System Identi cation 33 then the poles of the discretetime system obtained using the Tustin transformation are given by 2 snTs 7 231Ts 7 232Ts 21 7 227 2731115 2732115 2 7 i l n 2 7 snTs When Ts is much smaller than the poles all the discretetime poles 2 become very close to one regardless of the continuoustime dynamics The same is true for the zeros Therefore when on identi es a discretetime system using a very high sampling fre quency it is common to obtain discretetime transfer functions with several poles and zeros close to one these poles correspond to the actual poles of the continuoustime process warped by the Tustin transformation as well as several other poles and zeros away from one that arise due to highfrequency noise Generally these need to be discarded F When the sampling frequency is very large the values of yk and yk 1 are very close to each other and the difference between them is very much affected by the quantization interval introduced by the analogue to digital conversion This leads to a very poor predictive power by ARX models that try to predict the value of yk based on the last few past values of the input and output With larger sampling times the timehorizon over which an ARX model attempts to make predictions is larger and consequently the outputs exhibit larger variations making them less sensitive to quantization noise In practice one wants to chose a sampling rate that is large enough so that the Tustin transformation is valid but sufficiently small so as to avoid the problems discussed above A good rule of thumb is to sample the system at a frequency 2040 times larger than the desired closedloop bandwidth for the system One should look at the range of values taken by yk over the horizon of the ARX model to make sure that quantization effect are not too noticeable Attention It sometimes happens that the hardware used to collect data allow us to sample the system at frequencies much larger than what is needed for identi cationioversampling In this case one generally needs to downsample the signals but this actually provides an opportunity to remove measurement noise from the signals Suppose that yk and Mk have been sampled with a period Tlow but one wants to obtain signals le and 110 sampled with period Thigh LTlow where L is some integer larger than one A simple method to obtain 90 and ri e consists of extracting only one out of each L samples of the original signals Mk 2AM 110 Lkl 42 However instead of discarding all the remaining samples one can achieve some noise re duction by averaging Eigi Matlab hint 15 A signal can by downesampled as in 42 using bary 1Lend ownesampling can also be achieved with the command barydownsample y L Matlab hint 16 A signal y can by downesampled as in 43 using baryy1Lendi2 y2Lend1 y3Lend3 Matlab hint 17 wnisampling with more sophisticated and longer tering can be achieved with the command resample eg baryresample y 1 L 1 34 Joao P Hespanha 90 W 43 Mk uLk71u kuLk1 44 or even longer averaging The downsampled signals obtained in this way exhibit lower noise than the original ones 4 5 Exercises Exercise 6 Input magnitude A Simulink block that models a nonlinear springmass damper system is provided I Use the Simulink block to generate the system s response to step inputs with amplitude 025 and 10 and no measurement noise F For each set of data use the leastsquares method to estimate the systems transfer function Try a few values for the degrees of the numerator and denominator polyno mials m and n Check the quality of the model by following the validation procedure outlines above Important write MATLAB scripts to automate these procedures These scripts should take as inputs the simulation data uk yk and the integers m n D5 Use the Simulink block to generate the system s response to step inputs and measure ment noise with intensity 001 For the best values of m and n determined above plot the SSE vs the step size Which stepsize leads to the best model B Exercise 7 Model order Use the data provided to identify the transfer function of the system Use the procedure outlined above to determine the order of the numerator and denominator polynomials Plot the largest and smallest singular value of R and the SSE as a function of 71 Important write MATLAB scripts to automate this procedure These scripts should take as inputs the simulation data uk yk and the integers m n D Exercise 8 Sampling frequency Consider a continuoustime system with transfer function 471392 P S 32 7rs 471392 H Build a Simulink model for this system and collect inputoutput data for an input square wave with frequency 25Hz and unit amplitude for two sampling frequencies Ts 25sec and Ts 0005sec F Identify the system s transfer function without downsampling System Identi cation 35 3 Identify the system s transfer function using the data collected With Ts i0005sec but downsampled Important Write MATLAB scripts to automate this procedures These scripts should take as inputs the simulation data Mk yk the integers m n and the downsampling period Li 4 Brie y comment on the results D Jo o P Hespanha Lecture 5 Nonparametric identi cation Summary This lecture presents basic techniques for nonparametric identi cation 1 Timedomain identi cation 2 Frequencydomain identi cation Contents 51 Nonparametric methods i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 37 5 2 Timedomain identi cation i i i i i i i i i i i i i i i i i i i i i i i i i i i i 38 53 Frequency response identi cation i i i i i i i i i i i i i i i i i i i i i i i i i 40 5 4 MATLAB hints i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 42 55 To probe further i i i i i i i 42 5 6 Exercises i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 43 51 Nonparametric methods Nonparametric identi cation attempts to directly determine the model of a system without assuming that its transfer function is rational and that we known the number of poles or zeros The following are typical problems in this class Problem 4 Nonparametric frequency response identi cation Determine the fre quency response Hem over a range of frequencies 9 E Qmim Qmaxli 01 Problem 5 Nonparametric time response identi cation Determine the impulse response Mk from time k 0 to k Ni Matlab hint 18 ffth computes fastrFourieretransform FFT of the signal hk k 6 12 N which provides the values of Hej0 for equally spaced frequencies 9 k e01N7 1 Sidebar 19 The main weakness of the impulseresponse method is that impulses are rarely representative of typical inputs that appear in closed loop so especially for nonlinear systems we estimate a regimen that ma be far from the dynamics that will appear in closed 38 Joao P Hespanha Problem 4 is useful for controller design based on loopshaping or the Nyquist criterion whereas Problem 5 is useful for controller design based on the ZieglerNichols rules to tune PID controllers For N large one can also recover the transfer function from the impulse response by taking the ztransform 00 N Ho Zhk Zhwyk s Emmi k0 k20 assuming that N is large enough so that Mk m 0 for k gt N Throughout this chapter we will assume that the system is BIBO stable ie that all poles of the transfer function H have magnitude smaller than one or equivalently that Mk converges to zero exponentially fast A detailed treatment of this subject can be found eg in 5 Chapter 7 52 Timedomain identi cation We start by discussing a few methods used to solve the time response identi cation Prob lem 5 521 Impulse response method The impulse response of a system can be determines directly by starting with the system at rest and applying an impulse at the input Mk oz k 0 0 k y 0 The output will be a scaled version of the impulse response possibly corrupted by some measurement noise n yk aMk Therefore 210 710 h k a a Hopefully the noise term nka is small when compared to Mk and one can use the following estimate for the impulse response Mk L0 0 ke12mN Attention The choice of oz is generally critical to obtain a good estimate for the impulse response System Identi cation 39 i lal should be large to make sure that nkoz is indeed negligible when compared to i ii lal should be small to make sure that the process does not leave the region where a linear model is valid and where it is safe to operate it openloop As one increases a a simple practical test to check if one is leaving the linear region of operation is to do identi cation both with some oz gt 0 and 7a lt 0 In the linear region of operation the identi ed impulseresponse does not change See also the discussion in Section 42 522 Step response Sidebar 20 In general The impulse response of a system can also be determined by starting with the system at steps are more representative rest and applying an step at the input gan1gtt ilsriepfgbtaf35 response is a better method a k 2 0 than the impulse response 0 k lt 0 The ztransform of the output will be given by a 172 1 m Hlt2gtUltzgt Na Hz Nltzgt where Nz denotes the ztransform of measurement noise Solving for we obtain oz oz Taking inverse ztransforms we conclude that W yltkgteyltke 1 7 nltkgtenltke 1 oz oz Hopefully the noise term w is small when compared to Mk and we can use following estimate for the impulse response MPW kg 12MN Attention The main weakness of the stepresponse method is that it can amplify noise because in the worst case W can be twice as large as 7 20 when nk 71 771k Although this may seem unlikely it is not unlikely to have all the nk independent and identically distributed with zero mean and standard deviation 7 In this case Sidebar 21 Why Sthev 771k 3 Sthev k k 1 7141quot oz oz oz oz which means that the noise is ampli ed by approximately 41 D Matlab hint 19 The identi cation toolbox command impulse performs impulse response estimation for arbitrary inputs Sidebar 22 To minimize the errors caused by noise the amplitude 0 should be large However it should still be su 39iciently small so that the process does not leave the linear regime 40 Joao P Hespanha 523 Other inputs One can determine the impulse response of a system using any input and not just a pulse or stepinput Take an arbitrary input Mk with ztransform The ztransform of the output will be given by Yz Nz where Nz denotes the ztransform of measurement noise Solving for we obtain Hltzgt U1ltzgtltYltzgt 7 M2 Taking inverse ztransforms we conclude that W Z llU 1zlyki mm where the denotes the convolution operation Hopefully the noise term is small when compared to Mk and we can use following estimate for the impulse response Mk zilwilw a yk k e 12 N 53 Frequency response identi cation We now consider a couple of methods used to solve the frequency response identi cation Problem 4 531 Sinewave testing Suppose one applies an sinusoidal input of the form ukozcos k Vk E 12N 51 Assuming that has no poles at 2 em the measured output is given by gas aAn m 4m saw mas 52 where 1 Ag lHejnl and 4 2 function at 2 e AHem are the magnitude and phase of the transfer 5 6k is a transient signal that converges to zero as fast as C yk where y is the magnitude of the pole of with largest magnitude and c some constant an 3 Mk corresponds to measurement noise When nk is much smaller than 1A9 and for k suf ciently large so that 6k is negligible we have yk m 1A9 cos 9k This allows us to recover both Ag lHejnl and 19 AHem from the magnitude and phase of Repeating this experiment for several inputs with different frequencies 9 one can estimate the frequency response H em over the range of frequencies of interest System Identi cation 41 532 Correlation method Especially when there is noise it may be difficult to determine the amplitude and phase of yk by inspection The correlation method aims at accomplishing this task Suppose that the input uUs in 51 was applied and that after some time No we can neglect the transient term 6k from the output equation 52 In the correlation method we compute M2 N Kg Ze mkyag 1c cosQk 7 jsinQkyk a H 1 N N 1A9 Z ltcosQkcosQk 7 j sinQkcosQk Ze jnkk s k21 k21 Using the trigonometric formulas cosacosb cosa 7 b cosa b and sinacosb sina 7 b sina b we conclude that N KQ gAQ Z cosdm cos29k we jsin gtg 7jsin29k Mn N Zemkek 1c 2H M2 N 3 m E 72jnk71 gtn 71m 2AQZe 2Agge e sknk k1 1 EAQEM EAnen39dm 572ij 2 2 M2 e mwe s a 1 For large N the rst term dominates over the remaining ones as long as ri e and ek do Sidebar 2 3 Why7 not contain pure sinusoidal terms of frequency 9 Therefore Kg m Age gtn Hemi This leads to the following estimate for Hem pgNE 0 Attention The main weaknesses of the frequency response methods are i To estimate Hem at a single frequency one still needs a long input ie N large In practice this leads to a long experimentation period to obtain the Bode plot over a wide range of frequencies ii It requires that we apply very speci c inputs sinusoids to the process This may not be possible or safe in certain applications Especially for processes that are openloop unstable 42 Joao P Hespanha 54 MATLAB hints Matlab Hint 19 impulse The command impulse from the identi cation toolbox per forms impulse response estimation from data collected With arbitrary inputs To use this command one must 1 Create a data object that encapsulates the inputoutput data using dataiddatayu Where u and y are vectors With the input and output data F Plot the estimated impulse response using impulse dat or compute the response using modelimpulse data Where In is a MATLAB object With an impulseresponse model The impulse response is given by model b 55 To probe further Sidebar 21 Noise in stepresponse method When nk and nk71 are independent ave Varnk 7 nk 7 1 Varnk Varnk 7 Therefore Sthev nk 720 71 Var 720 71 Varnk Varnk 7 1 a2 When both variables have the same standard deviation 7 we obtain 0 7 7 2 StheV M 3W2L2 amp g 0 OZ Sidebar 23 Correlation method In the correlation method one computes N V V N V N V Kg TQHQJQ gAQe N ZE QJWC Zeimkk s k21 k21 System Identi cation 43 The rst summation is given by 1 7 eanN N 4ij 7 E e 7 V 17 5 219 k1 and is therefore bounded as N A 0 as long as Q y 0 As N A 0 the second summation converges to wmw mas EM NM k21 Where and Nz denote the ztransforms of 5k and nk respectively Therefore this term is also bounded as long as nk and 5k do not contain pure sinusoidal terms of frequency 9 We therefore conclude that the term Hem dominates as N A 00 D 56 Exercises Exercise 9 Step response A Simulink block that models a pendulum with viscous friction is provided 1 Use the Simulink block Without noise to estimate the system s impulse response using both the impulse and step response methodsi What is the largest value oz of the input magnitude for Which the system still remains Within the approximately linear region of operation Hint to make sure that the estimate is good you can compare the impulses response from steps and impulsesi F Turn on the noise in the Simulink block and repeat the identi cation procedure above for a range of values for at For both methods plot the error Eh Hho s 7 Mk vs a Where ho s denotes the estimate obtained for an input With magnitude as What is your conclusion Take the estimate determined in l for the noisefree case as the ground truth has Exercise 10 Correlation method Modify the Simulink block provided for Exercise 9 to accept a sinusoidal input 1 Estimate the system s frequency response Hem at a representative set of log spaced frequencies 0 using the correlation methodi Select an input amplitude oz for Which the system remains Within the approximately linear region of operation 2 Turn on the noise in the Simulink block and repeat the identi cation procedure above for a range of values for at Plot the error 21 lIHaejnz 7 Hemz I vs a Where Joao P Hespanha Haejnz denotes the estimate obtained for an input With magnitude 01 What is your conclusion Take the estimate Hejn determined in 1 for the noisefree case as the ground truth Hemi i Compare your best frequency response estimate Hejn With the frequency responses that you would obtain by taking the Fourier transform of the best impulse response Mk that you obtained in Exercise 9 Hint Use the MATLAB command fft to compute the Fourier transformi This command gives you Hejn from Q 0 to Q 27r so you should discard the second half of the Fourier transform corresponding to 9 gt 7r

### BOOM! Enjoy Your Free Notes!

We've added these Notes to your profile, click here to view them now.

### You're already Subscribed!

Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'

## Why people love StudySoup

#### "There's no way I would have passed my Organic Chemistry class this semester without the notes and study guides I got from StudySoup."

#### "I used the money I made selling my notes & study guides to pay for spring break in Olympia, Washington...which was Sweet!"

#### "I was shooting for a perfect 4.0 GPA this semester. Having StudySoup as a study aid was critical to helping me achieve my goal...and I nailed it!"

#### "Their 'Elite Notetakers' are making over $1,200/month in sales by creating high quality content that helps their classmates in a time of need."

### Refund Policy

#### STUDYSOUP CANCELLATION POLICY

All subscriptions to StudySoup are paid in full at the time of subscribing. To change your credit card information or to cancel your subscription, go to "Edit Settings". All credit card information will be available there. If you should decide to cancel your subscription, it will continue to be valid until the next payment period, as all payments for the current period were made in advance. For special circumstances, please email support@studysoup.com

#### STUDYSOUP REFUND POLICY

StudySoup has more than 1 million course-specific study resources to help students study smarter. If you’re having trouble finding what you’re looking for, our customer support team can help you find what you need! Feel free to contact them here: support@studysoup.com

Recurring Subscriptions: If you have canceled your recurring subscription on the day of renewal and have not downloaded any documents, you may request a refund by submitting an email to support@studysoup.com

Satisfaction Guarantee: If you’re not satisfied with your subscription, you can contact us for further help. Contact must be made within 3 business days of your subscription purchase and your refund request will be subject for review.

Please Note: Refunds can never be provided more than 30 days after the initial purchase date regardless of your activity on the site.