%   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