% IDDEMOPR: Demo of process models echo off % Lennart Ljung % Copyright 1986-2003 The MathWorks, Inc. % $Revision: 1.7.4.1 $ $Date: 2004/04/10 23:16:08 $ if isempty(findstr(matlabpath,'simulink')) error('This demo requires SIMULINK') end echo on % This demo illustrates how to build simple process models often used in % process industry. % They are of the basic type 'Static Gain + Time Constant + Time Delay'. % pause % Strike a key to continue % % The models covered are of the general types % 1 + Tz*s % P(s) = Kp * exp(-Td*s) * ---------------------- % (1 + Tp1*s)(1 + Tp2*s) % % or an integrating process % 1 + Tz*s % P(s) = Kp * exp(-Td*s) * ------------------------ % s(1 + Tp1*s)(1 + Tp2*s) % % where the user can determine the number of real poles (0, 1, 2 or 3), as well % as the presence of a zero in the numerator, the presence of an integration % and the presence of a time delay (Td). Underdamped (complex) poles can % also be allowed. % % The particular type of model is defined by the letters P (for process model), % D (for time delay) Z (for a zero) and I (for integration). An integer will denote % the number of poles. The models are generated by the command IDPROC. % % This should be clear from the following examples. pause % Strike a key to continue idproc('P1') pause % Strike a key to continue idproc('P2DIZ') pause % Strike a key to continue idproc('P0ID') pause % Strike a key to continue % Consider the system described by the following SIMULINK scheme. open_system('iddempr1') set_param('iddempr1/Random Number','seed','0') % The red part is the system, the blue part is the controller and % the reference signal is a swept sinusoid (a chirp signal). % The data sampling time is set to 0.5 seconds. pause %Strike a key to continue % The system can in terms of SITB objects be described by m0 = idpoly(1,0.1,1,1,[1 0.5],'Ts',0,'InputDelay',1.57,'NoiseVariance',0.01); % or as an idproc model m0p = idproc('p1d','Kp',0.2,'Tp1',2,'Td',1.57); pause %Strike a key to continue m0p % To estimate the model, first simulate the system pause, sim('iddempr1') % Strike a key to simulate dat1e = iddata(y,u,0.5); % The IDDATA object % Let us look at the data: pause, figure(1), plot(dat1e), % Strike a key to continue % Then estimate the model by applying the PEM command to the data with % a model structure determined by IDPROC: pause %Strike a key to continue m1 = pem(dat1e,'p1d'); % As you see, it is necessary to give an upper bound for the % delay. That can be done as in % m1 = pem(dat1e,'p1d','Td',{'max',2}); % % Check the result: pause %Strike a key to continue m1 pause %Strike a key to continue % To get information about uncertainties, use present(m1) % The resulting model M is an IDPROC model object to which all of % the toolbox's model commands can be applied: pause % Strike a key to continue step(m1,m0) pause bode(m1,m0,'sd',5,'fill') pause compare(dat1e,m0,m1) % etc. pause % Strike a key to continue % It may be important (at least for slow sampling) to consider the % intersample behavior of the input data. To illustrate this, let % us study the same system as before, but without the sample and % hold circuit: % open_system('iddempr5') pause % Strike a key to continue % Simulate this system with the same sampling interval: pause, sim('iddempr5') % Strike a key to simulate dat1f = iddata(y,u,0.5); % The IDDATA object % and estimate the model: m2 = pem(dat1f,'p1d','Td',{'max',2}); pause % strike a key to continue m2 % This model has a slightly less precise estimate of the delay than % the previous one, m1: [m0p.td.value, m1.td.value, m2.td.value] step(m0,m1,m2) pause % strike a key to continue % However, by telling the estimation process that the intersample % behavior is first-order-hold (an approximation to the true % continuous) input, we do better: dat1f.InterSample = 'foh'; m3 = pem(dat1f,'p1d','Td',{'max',2}) % Compare the four models m0 (true) % m1 (obtained from zoh input) % m2 (obtained for continuous input, with zoh assumption) and % m3 (obtained for the same input,but with foh assumption): pause % strike a key to continue [m0p.td.value, m1.td.value, m2.td.value, m3.td.value] compare(dat1e,m0,m1,m2,m3) pause % strike a key to continue step(m0,m1,m2,m3) pause % strike a key to continue % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close_system('iddempr1'),close_system('iddempr5') % Let us now consider a more complex process, with integration, that is % operated in closed loop: open_system('iddempr2') % The true system can be represented by m0 = idproc('p2zdi','Kp',1,'Tp1',1,'Tp2',5,'Tz',3,'Td',2.2); % The process is controlled by a PD regulator with limited input amplitude % and a zero order hold device. The sampling interval is 1 second set_param('iddempr2/Random Number','seed','0') pause, sim('iddempr2') % Strike a key to simulate dat2 = iddata(y,u,1); % The IDDATA object % % Two different simulations are made, the second one for % validation purposes. % set_param('iddempr2/Random Number','seed','13') sim('iddempr2') dat2v = iddata(y,u,1); % Let us look at the data: % pause, plot(dat2), % Strike a key to continue % To estimate the model use pause, m2 = pem(dat2,'p2zdi','td',{'max',5}) % Strike a key to continue pause % Strike a key to continue compare(dat2v,m2,m0) % Gives very good agreement with data pause % Strike a key to continue bode(m2,m0) pause % Strike a key to continue % Compare also with the parameters of the true system present(m2) pause, m0 % Strike a key to continue [m0.par m2.par] % A word of caution. Identification of several real time constants may sometimes % be an ill-conditioned problem, especially if the data are collected in closed loop. % % To illustrate this, let us estimate a model based on the validation data: % pause, % Strike a key to continue m2v = pem(dat2v,'p2diz','td',{'max',5}) % This model has much worse parameter values. On the other hand, it performs nearly % identically to the true system m0 when tested on the other data set dat2: pause % Strike a key to continue compare(dat2,m0,m2,m2v) pause % Strike a key to continue % Suppose we know from other sources that one time constant is 1: m2v.Tp1 = 1; m2v.Tp1.status='fix'; % We can fix this value, while estimating the other parameters by % ('trace','full' means full info about the iterations to the screen) % pause, m2v = pem(dat2v,m2v,'trace','full') % Strike a key to continue % pause % Strike a key to continue % This also indicates that simple approximation should do well on the data: m1 = pem(dat2v,'p2d') pause % Strike a key to continue % compare(dat2,m0,m2,m2v,m1) % However, m1 does not contain any integration, so the open loop long % time range behaviour will be quite different: % pause % Strike a key to continue step(m0,m2,m2v,m1) pause % Strike a key to finish demo bdclose('iddempr2'), close(1) echo off