Lecture
The System Identification Toolbox package contains tools for building mathematical models of dynamic systems based on observed input and output data. It has a flexible graphical interface that helps organize data and build models [2, 12, 29, 30]. The identification methods included in the package are applicable to a wide range of problems: from the design of control systems and signal processing to the analysis of time series and vibration. The main features of the package are:
• a simple and flexible interface;
• preliminary data processing, including pre-filtering, removal of trends and offsets; selection of the data range for analysis;
• autoregression methods;
• response analysis in the time and frequency domains;
• display of the zeros and poles of the system's transfer function;
• residual analysis when testing the model;
• construction of complex diagrams, such as the Nyquist diagram,
and others.
The graphical interface simplifies data preprocessing as well as the interactive process of model identification. It is also possible to work with the package in command mode and using the Simulink extension. Operations for loading and saving data, selecting a range, and removing offsets and trends require minimal effort and are located in the main menu.
The presentation of data and identified models is organized graphically in such a way that during the interactive identification process the user can easily return to a previous step. For beginners there is an option to preview
the following possible steps: the graphical tools make it possible to find any of the previously obtained models and assess its quality in comparison with other models.
Starting from output and input measurements, one can build a parametric model of the system that describes its dynamic behavior. The package supports all traditional model structures, including autoregression, the Box–Jenkins structure, and others. It supports linear state-space models, which can be defined in either discrete or continuous space. These models can include an arbitrary number of inputs and outputs. The package contains functions that can be used as test data for the identified models. Identification of linear models is widely used in control system design, when a model of the plant needs to be created. In signal processing problems, models can be used for adaptive signal processing. Identification methods are also successfully applied to financial applications.
2.3.1. Types of models in the System Identification Toolbox package
The discrete-time models used in the System Identification Toolbox package, which account for the effect of observation noise, can be represented as follows:
• the AR (AutoRegressive) model:
A(z) y(t) = e(t),
where e(t) is discrete white noise; A(z) =1+ a1z−1 + a2 z−2 +K+ ana z−na ;
• the ARX model (AutoRegressive with external input):
A(z) y(t) = B(z)u(t) +e(t),
where B(z) =1+b1z−1 +b2 z−2 +K+bnb z−nb ;
• the ARMAX model (AutoRegressive-Moving Average with external input) — the autoregressive moving-average model:
y(t) +a1 y(t −1) +K+ ana y(t −n) =
= b1u(t) +b2u(t −1) +K+bnb y(t −m) +e(t), A(z) y(t) = B(z)u(t −nk) +C(z)e(t),
81
where nk is the delay (lag) value;
C(z) =1+c1z−1 +c2 z−2 +K+cnc z−nc ;
• the Output–Error model:
y(t) =
B(z)
u(t −nk) +e(t) ,
F(z)
where F(z) =1+ f z−1
+ f
2
z−2 +K+ f
nf
z−nf ;
1
• the Box–Jenkins (BJ) model:
y(t) =
B(z)
u(t
−nk) +
C(z)e(t) ,
F(z)
D(z)
where D(z) =1+ d1z−1 + d2 z−2 +K+dnd z−nd .
These models can be regarded as special cases of the generalized parametric linear structure:
A(z) y(t) = FB((zz))u(t −nk) + CD((zz))e(t).
At the same time, all of them allow extension to multivariable plants (having several inputs and outputs);
• the state-space model: x(t +1) = Ax(t) + Bu(t),
y(t) = C(x(t) + Du(t) +v(t),
where A, B, C, D are the coefficient, input, output, and feedthrough matrices of the system; v(t) is correlated observation white noise. Another (canonical) form of representing this model is also possible:
x(t +1) = Ax(t) + Bu(t) + Ke(t),
y(t) = C(x(t) + Du(t) +e(t),
where K is some weighting matrix; e(t) is discrete white noise. In its operation, the System Identification Toolbox package uses three
internal types of matrix model representation, which, using the package's operators and functions, are converted into all of the aforementioned types of plant models:
82
• theta format (for time-domain models);
• frequency format (for frequency-domain models);
• zero-pole format.
2.3.2. Using the System Identification Toolbox package
to build a system model
The main operators and functions of the System Identification Toolbox package are typed at the MatLab command line or can be used when writing programs for m-files.
A list of the main operators and functions and their descriptions is given in Appendix 1.
The rules for working with these blocks are similar to the rules for working with other Simulink blocks. The resulting model is displayed
in the main MatLab window.
As an example of using the System Identification Toolbox package for identifying technical systems, a spray dryer is considered [12].
It implements a certain technological process, in which
the input action on the system is the gas flow rate, m3/h, and the output controlled parameter is temperature, °C. The identification process includes the following stages:
• a priori analysis of the plant in order to select the model structure;
• conducting a preliminary study of the plant in order to refine the estimate of the model structure;
• developing a methodology for the main experimental study of the plant, drawing up an experiment plan;
• conducting the main experimental study to obtain an array of data (ui, yi);
• mathematical processing of the data array (using the System Identification Toolbox package) in order to determine the model parameters and its adequacy, the confidence bounds of the parameters, and the output coordinate of the model.
At the same time, in the process of studying the system it is necessary to accept
certain assumptions that allow applying the well-developed apparatus for analyzing stationary, linear plants:
• the technical plant is a system with lumped parameters;
83
• the technical plant is stationary, i.e., its static and dynamic properties are invariant over time;
• the equations of the models of the system under study are linearized locally, i.e., for small deviations ±∆yi from the chosen "operating" point
(operating mode) of the plant.
The data array ui, yi is formed as a result of the labor-intensive operation of
decoding the recorded charts from the measuring system's instruments. However, the wide development of microprocessor and computing technology and its introduction into production technological processes have made it possible to significantly improve the technical support for system identification. Processing the data array using the System Identification Toolbox package involves the following stages:
• processing and converting data in order to create a file;
• analysis of experimental data in order to preliminarily determine the main characteristics of the plant;
• parametric estimation of data in order to create various types of models in theta format;
• specifying the model structure;
• modifying and refining the model structure (if necessary);
• checking the adequacy of and comparing various models in order to select the best one;
• converting the model from theta format into a form convenient for further use in the analysis and synthesis of the control system.
At each stage of identification it is possible to graphically display the modeling results and extract the necessary information about the plant.
2.3.2.1. Data processing when building a plant model
Building an analytical model of a technical plant begins with a description of the operating principle of the plant — the dryer. Slurry is fed into the dryer, where it is sprayed. Injection burners create a high temperature in the material spraying zone. The sprayed particles, losing moisture, are already collected as powder at the bottom of the dryer, from where they go directly to the hopper above the presses.
84
The main task is to build the model in the form of a transfer function. It is also necessary to take into account the presence of random disturbances (i.e., factors unavoidable in a real setting), the effect of which must be accounted for in the modeling. In particular, such disturbances may include random changes in the intensity of slurry feed into the dryer or random changes in the quality of fuel in the burners, etc. These changes must be taken into account when building the dryer model. Theoretically, the problem reduces to studying the behavior of the system taking into account the effect of random factors on it.
An analytical model of complex systems can be built based on an array of input and output data obtained as a result of a physical experiment conducted on the system. In this case, the effect of random factors must be modeled.
The statistical data obtained while studying the thermal plant are placed in a separate file (named Project), which includes a data array consisting of 1000 values of the input parameter — gas flow rate, m3/h — and 1000 values of the output parameter — temperature in the plant, °C.
To load the file into the MatLab Workspace, the following command must be executed in command-line mode:
>>load Project.
As a result of executing the command, an array of input variables u2 and an array of output variables
y2 will appear in the Workspace. The sampling interval (T0 = 0.08) must be specified additionally at the command line:
>> ts = 0.08.
For further use in the System Identification Toolbox package, this data array must be combined into a single file containing the necessary information about the input and output parameters of the plant, their values, and units of measurement.
To combine the source data into a single file dan.m, the command used is
>> dan=iddata(y2,u2,ts).
85
The result of executing the command is commented on by the following
MatLab phrase:
Time domain data set with 1000 samples. Sampling interval: 0.08 Outputs Unit (if specified) у1 Inputs Unit (if specified) u1.
The generated file dan.m indicates that it contains the results of 1000 measurements with a sampling interval of 0.08 s. The input variable is the array u1, and the output parameter is y1.
For clarity of the generated file, it is necessary to introduce into its structure designations for the input and output data, indicating the parameter units:
>>dan.outputn='temperature'; >>dan.input='gas flow rate'; >>dan.inputUnit='m3/h', >>dan.outputUnit='deg.C 100'.
As a result, the generated data file dan.m has the following
form:
Time domain data set with 1000 samples. Sampling interval: 0.08 Outputs Unit (if specified) temperature deg.C 100 Inputs Unit (if specified) gas flow rate m3/h.
Complete information about the file dan.m can be obtained using the command
>>get(dan).
To graphically display the data, the command plot(dan) or the command idplot(datta) is used; however, in the latter case the graphs will not contain information about the variable names and their units. The source data using the command >>plot(dan) is shown in fig. 2.2.
For further use of the obtained source data, it is necessary to carry out their preliminary processing in order to remove the trend (constant component) from the data set and, if necessary, filter the data using the tools available in the System Identification Toolbox package.
To remove the trend, the function dtrend is used:
>> zdan=dtrend(dan).
86
Fig. 2.2. Source data for identification of the control plant
Executing the function will produce the following entry at the command line
Time domain data set with 1000 samples. Sampling interval: 0.08 Outputs Unit (if specified) temperature deg.C 100 Inputs Unit (if specified) gas flow rate m3/h.
The file zdan will appear in the system's Workspace. A new file zdan.m has been obtained, in which the constant
component of the signals is absent (fig. 2.3). The figure is obtained after executing the command
>> plot(zdan)).
The file will subsequently be used to build models of the control plant.
In addition to the trend removal command mentioned, the System Identification Toolbox package has other experimental data processing functions, which are given in the description of the System Identification Toolbox package. These functions are applied in cases where a preliminary analysis of the identification plant has been carried out and possible disturbances or certain other dynamic characteristics
87
have been determined, or when it becomes necessary to change the sampling interval in the case of increased error in the representation of the system model during parametric estimation.
Fig. 2.3. Source data for plant identification in the absence of a trend
The next stage of identification is determining the statistical and frequency characteristics of the arrays of source data.
2.3.2.2. Estimating the statistical
and frequency characteristics of the source data
As already noted above, when forming the array of source data using a physical experiment on the technical control plant, the input signal acting on the plant was represented as a random process with zero mathematical expectation (i.e., centered after removing the trend). We will consider the process to be ergodic, which is necessary for practical applications of the theory of random processes, since it makes it possible to judge the statistical characteristics of a random process from a single sufficiently long realization of it. In accordance with the properties of a stationary ergodic process, any statistical characteristic obtained by averaging over the ensemble of possible
88
realizations can, with a probability arbitrarily close to unity, be obtained by averaging over a sufficiently long time interval from a single realization of the random process. Therefore, any realization of the source data can be used to obtain the statistical characteristics of the arrays of source data, since during the planning and conducting of the experiment it is impossible to say in advance which realization the process will follow.
To characterize the relationship between the values of a random process at different points in time, the correlation function Rx(τ) and the spectral density of the random process Sx(ω) are used.
The statistical properties of the relationship between two random processes X(t) and G(t) can be characterized by the mutual correlation function
Rxg(τ).
The System Identification Toolbox package has four functions — cra, etfe, conf, and spa — for nonparametric estimation of the set of experimental data. The function cra performs the calculation of auto- and cross-correlation functions, and estimation of the impulse response by the correlation analysis method for a single-input single-output plant based on the array of experimental data. This function is written as follows:
cra(z);
[ir,R,cl] = cra(z, М, na, plot), cra(R).
Here z is the matrix of experimental data of the form z = [y2 u2]; y2 is the column vector corresponding to the output data; u2 is the column vector corresponding to the input data; М is the maximum value of the discrete argument for which the impulse response estimate is calculated (by default М = 20); na is the order of the autoregressive model (polynomial order) used to calculate the parameters of the whitening filter (by default na = 10). When na = 0, the untransformed input signal is used as the identifying signal; plot indicates whether a graph is produced:
if plot = 0, there is no graph,
if plot = 1, a graph of the obtained impulse response estimate is plotted together with the 99% confidence interval,
89
if plot = 2, graphs of all correlation functions are output. Return values:
ir — the estimate (vector of values) of the impulse response;
R — a matrix whose first-column elements are the values of the discrete argument, second-column elements are the values of the estimated autocorrelation function of the output signal, third-column elements are the values of the estimated autocorrelation function of the input signal, and fourth-column elements are the values of the estimated cross-correlation function.
For example, for the slurry dryer these values are as follows.
М and na are taken by default [], [].
>>[ir,R,cI]=cra(zdan,[], [],2).
Fig. 2.4 shows the results of calculating the autocorrelation function of the output signal (Covf for filtered y); the autocorrelation function of the input signal (Covf for prewhitened u); the cross-correlation function (Correlation from u to у); the impulse response
estimate (Impulse response estimate).
A more detailed graph of the impulse response can be obtained by running the function cra with a single argument zdan
(fig. 2.5):
>> cra(zdan).
It should be noted that on the graphs the abscissa shows time intervals τ = ti −ti−1, and the ordinate shows the values of the
correlation functions for the input u2 and output y2 signals; the values of the cross-correlation function and the impulse response.
From the characteristics obtained, it follows that as τ increases, a sharp drop in the correlation dependence of the input signal is observed, which indicates a weak relationship between sections of the process corresponding to arbitrary points in time (the process is closer to white noise, and the autocorrelation function is closer to a delta function). The output quantity, on the contrary, changes states from one moment to another more smoothly, and therefore the relationship between the previous and subsequent values of the output signal is closer than that of the input.
a b
c d
Fig. 2.4. Function graphs: a — autocorrelation function of the output signal; b — autocorrelation function of the input signal; c — cross-correlation function; d — impulse response
Fig. 2.5. Impulse response
To obtain the frequency characteristics of the experimental data, we will use the frequency-characteristic estimation functions. The function spa returns the frequency characteristics of a single-input single-output plant and estimates of the spectral density of its signals for the generalized linear plant model:
91
[g, phiv]= spa(z);
[g, phiv, z_spe]=spa(z, М, maxsize, Т);
here z is the matrix of source data; М is the width of the time window (by default М = min(30, length(z)/10), where length(z) is the number of rows of the matrix
z); Т is the sampling interval; maxsize is a parameter that determines the maximum size of the matrices created during the computation.
Return values:
g — the estimate of W(e jωT) in frequency format;
phiv — the estimate of the spectral density of the noise v(t);
z_spe — the matrix of spectral densities of the input and output signals.
The Bode diagram (magnitude and phase response) is built using the functions spa and bodeplot and the data obtained while studying the plant and contained in the file Project:
>>load Project;
>>х=[у2 u2];
>>S=spa(z);
>>bodeplot(g).
The modeling results without a confidence interval are shown in fig. 2.6.
The dependencies obtained confirm the high-frequency component of the values of the input and output signals. The limits of frequency variation on the graphs are set by default.
To obtain the frequency characteristics together with a confidence interval of width three standard deviations, the System Identification Toolbox package in MatLab offers the following options:
• setting the limits of frequency variation using the command
>>w=logspace(w1, w2, N),
where w1 is the lower bound of the frequency range (10w1); w2 is the upper bound of the frequency range (10w2); N is the number of points on the graph;
92
Fig. 2.6. Frequency characteristics of the control plant
• plotting the magnitude response, phase response, and S(ω) — the spectral density function of the noise e(t);
• calculating the g-estimate of the magnitude and phase response in frequency format and the phiv-estimate of the noise spectral density using the command
>>[g, phiv]=spa(z, [], w).
The graphs of the magnitude response, phase response, and S(ω) are plotted with a confidence interval of three standard deviations using the command
>> bodeplot([g р], 'sd',3,'fill'),
where 'sd' indicates a solid line for the confidence interval (by default this line is dashed); 3 is the width of the confidence band in standard deviations; 'fill' is the way the confidence band is filled (in yellow).
Plotting the magnitude and phase response using the functions spa, bodeplot, logspace, and the data obtained in the file Project with the corresponding confidence band is done as follows:
>> w =logspace(-2,pi,128); >>[jg, phiv]=spa(z, [], w);
93
>> bodeplot([g, phiv], 3, 'fill').
The modeling results are shown in fig. 2.7.
Fig. 2.7. Estimates of the magnitude and phase response together with the confidence interval
To plot the graph of the noise spectral density estimate
with a confidence interval, the following command is executed:
>>bodeplot([phiv], 'sd', 3, 'fill').
The modeling results are shown in fig. 2.8.
Fig. 2.8. Estimates of S(ω) together with the confidence interval
94
The resulting graph of the noise spectral density estimate together with the confidence interval indicates a uniform distribution of signal power over the frequency spectrum, followed by a drop in power at frequencies above 1.1 rad/s.
Next, parametric estimation must be performed.
2.3.3.3. Parametric estimation of data
Parametric estimation of experimental data is carried out in order to determine the parameters of a model of a given structure by minimizing a chosen quality criterion of the model (most often, the mean square of the mismatch between the plant outputs and its postulated model).
To perform parametric estimation, the array of experimental data must be conditionally divided into two parts (not necessarily equal):
>>zdanv=zdan(I:500);
>>zdane=zdan(501:1000).
The first part of the data array will be used for parametric estimation and building the system model. The second part will be needed for verification (quality checking) of the model, assessing the adequacy of the resulting model, and determining identification errors. It should be noted that parametric identification in the System Identification Toolbox package is performed in discrete form, and the resulting models are discrete.
The System Identification Toolbox package considers various types of models, which describe the identification plant with varying degrees of accuracy. To select the most suitable structure and type of model during parametric estimation of experimental data, the System Identification Toolbox package has special functions for:
• parametric estimation;
• specifying the model structure;
• modifying and refining the model structure and selecting the model structure.
The estimation function ar determines the autoregression parameters:
95
A(z) y(t) = e(t),
where A(z) =1+ a z−1
+ a z−2
+Ka
na
z−na , i.e., the coefficients of the polynomial
1
2
A(z) when modeling scalar time sequences. The function has the syntax
th = ar(у, n)
or another form of writing that allows changing the modeling parameters:
[th, refl]=ar(y, n, approach, win, maxsize, Т),
where у is the column vector of data containing N elements; n is the model order (number of estimated coefficients); approach is an argument (string variable) determining the estimation method:
• '1s' — least squares method;
• 'yw' — Yule-Walker method;
• 'burg' — Burg's method (a combination of the least squares method
with minimization of the harmonic mean);
• 'gl' — a method using the harmonic mean.
If any of these values ends with a zero (for example, burg0), the calculation is accompanied by estimation of the correlation functions.
win — (string variable) used in the case of missing data:
• win='now' — only the available data is used (used by default, except when approach = 'yw');
• win='prw' — missing initial data is replaced with zeros, so that summation starts from the zero time instant;
• win='prw' — missing subsequent data is replaced with zeros, so that summation extends to the time instant N+n;
• win='ppw' — both initial and subsequent missing data are replaced with zeros (used in the Yule-Walker algorithm);
maxsize determines the maximum dimensionality of the problem; Т is the sampling interval.
Return values:
th — the resulting autoregression model in theta format (the internal matrix format for representing parametric models of the
System Identification package);
96
relf — information about the coefficients and the loss function.
To use the parametric estimation function ar, it is necessary to extract the output variable y from the array of experimental data recorded in the file dan, using the command
>>y=dan.у;
which is equivalent to the command
>>y=get(dan,'у'); >>th=ar(y,4);
Discrete time IDPOLY model: A(q)y(t)=e(t)
A(q)=1 -2.195q^-1+1.656q^-2-0.4167q^-3-0.04415q^- 4 Estimated using AR ('fb'/'now') from data set y
Loss function 0.0104048 and FPE 0.0104884 Sampling interval: 1.
Complete information about the autoregression model th can be obtained using the command
>> present(th).
Discrete-time IDPOLY model: A(q)y(t)=e(t) A(q)=1-2.195(±0.03165 q^-1+1.656(±0.07533 q^-2- -0.4167(±0.07542)q^-3-0.04415(±0.03174)q^- 4 Estimated using AR ('fb'/'now') from data set y Loss function 0.0104048 and FPE 0.0104884 Sampling interval: 1
Created: 07-Jan-2012 08:47:27
Last modified: 07-Jan-2012 08:47:27.
The information shows that the model is discrete and that a forward-backward method (a variant of the least squares method) is used to estimate its parameters, as indicated by the string variable 'fb' (used by default); only the available data y is required to build the model, as indicated by the string variable 'now' (used by default); the following are determined: the loss function as the residual sum of squared errors, and the so-called theoretical
Akaike Information Criterion (AIC), FPE; and the
sampling interval.
97
The next function, arx, estimates the parameters of the AR and ARX model, represented by the relationship
A(z) y(t) = B(z)u(t) +e(t)
or in expanded form:
y(t) + a1 y(t −1) +K+ ana y(t −n) =
= b1u(t) +b2u(t −1) +K+bnb y(t −m) +e(t).
Here and below, e(t) is discrete white noise.
B(z) =1+b1z−1 +b2 z−2 +K+bnb z−nb .
The function has the following syntax: dar=arx(z,nn)
or another form of writing that allows changing the modeling parameters:
dar=arx(z,nn,maxsize,Т),
where z is the experimental data; nn is the set of specified model parameters (the argument nn contains three parameters: na — the order (number of coefficients) of the polynomial A(z); nb — the order of the polynomial B(z); nk — the delay value); maxsize is the maximum dimensionality of the problem; Т is the sampling interval.
When choosing a polynomial, the question of its degree arises. It is known that as the order of the polynomials increases, the degree of adequacy of the model to the real plant improves. However, this results in cumbersome expressions and increases the modeling time. Therefore, to find the optimal polynomial order, the model structure selection functions can be used.
The function arxstruc computes the loss functions for a number of different competing single-output ARX models:
v=arxstruc(ze,zv,NN)
or
v=arxstruc(ze,zv,NN,maxsize);
where ze, zv are, respectively, the matrices of experimental data for model estimation and verification; NN is the matrix specifying the competing structures, with rows of the form nn=[na nb nk]; maxsize is the maximum dimensionality of the problem.
98
The returned value v is a matrix whose upper elements in each column (except the last) are the loss function values for the ARX models, the structure of which is indicated by the subsequent elements of the columns (i.e., each column corresponds to one model). The first element of the last column is the number of experimental data values used for model verification.
The function selstruc selects the best model structure from a range of possible options:
[nn,vmod]=selstruc(v), [nn,vmod]=selstruc(v,с),
where v is the matrix returned by the function arxstruc; с is a string variable determining whether a graph is output or the criterion for selecting the best structure:
• с='plot' — a graph is output of the loss function versus the number of estimated model coefficients;
• с='log' — a graph of the logarithm of the loss function is output;
• с='aic' — no graph is output, but the structure minimizing the theoretical Akaike Information Criterion is returned;
• с='mdl' — the structure providing the minimum of Rissanen's Minimum Description Length criterion is returned;
• с=a — the structure that minimizes the value of the loss function is selected:
vmod = v(1+a(d / N)),
where v is the value of the loss function; d is the number of estimated model coefficients; N is the size of the sample of experimental data used for estimation.
Return values: nn — the selected structure;
vmod — the value of the corresponding criterion. For example, for the Project data one can:
• set the limits for varying the model order: >>NN=struc(1:10,1:10,1);
• calculate the loss functions: >>v=arxstruc(zdane,zdаnv,NN);
99
• select the best structure of polynomial orders: >>[nn,vmod]=selstruc(v,'plot'),
where 'plot' is a string variable determining the output of a graph of the loss function versus the number of estimated coefficients of the
model (fig. 2.9).
Fig. 2.9. Model structure selection window
In the window that appears, the bars indicate the value of the loss function. When the cursor is moved to the corresponding bar, the values of the polynomial orders na, nb, nk are displayed in the right-hand field of the window. Recommendations for choosing the bar color appear in the graph field.
Instead of the string variable 'plot', the following options are possible:
• 'log' — a graph of the logarithm of the loss function is output;
• 'aic' — no graph is output, but the structure minimizing the theoretical Akaike Information Criterion
(AIC) is returned, FPE vmod:
vmod =
v
,
+2(d
))
(1
N
where v is the value of the loss function; d is the number of estimated model coefficients; N is the amount of experimental data used for estimation;
• 'mdl' — the structure providing the minimum of Rissanen's Minimum Description Length criterion is returned:
vmod =
v
;
+log(N)(d
))
(1
N
• with the string variable set to some numerical value a, the structure is selected that minimizes:
vmod =
v
.
+a(d
))
(1
N
The selection of the best structure of polynomial orders can also be done using a simpler command:
>>nn=selstruc(v,0)
MatLab returns
nn=7 9 1.
Taking into account the selected model structure, let us determine the form of the ARX model by executing the function arx:
>>darx=arx(zdanv,nn).
A matrix of 100 columns and 4 rows is returned with the values of various criteria:
vmod=…………. .
A discrete model is returned, presented in theta format (the internal form of matrix models).
Discrete-time IDPOLY model: А(q)у(t)=B(q)u(t)+e(t)
A(q)=1- 0.9601q^-1+0.03384q^-2-0.1035 q^-3+0.089q^-4- 0.02827^q-5+
+0.1383q^-6-0.06615q^-7 B(q)=0.0009708q^-1+0.003863q^-2+0.06428q^-3+0.06245q^-4+ +0.01997q^-5-0.007896q^-6-0.01869q^-7-0.01931 q^-8-0.007944q^-9 Estimated using ARX from data set zdanv
Loss function 0.00139284 and FPE 0.00148493 Sampling interval: 0.08.
101
The function armax estimates the parameters of the ARMAX model: >>darmax=armax(zdanv,[2 2 2 1]).
Arguments of the function:
zdanv — the vector of experimental data;
[nа nb nc nk] — the polynomial degrees and the delay value. A discrete model is returned in theta format:
Discrete-time IDPOLY model: A(q)y(t=B(q)u(t)+C(q)e(t) A(q=1-1.648q^-1+0.7054q^'-2 B(q)=0.005017q^-1+0.04556q^-2 C(q)=1-0.278q^-1-0.4333q^-2
Estimated using ARMAX from data set zdanv Loss function 0.00723841 and FPE 0.00741677 Sampling interval: 0.08.
The function oe estimates the parameters of the OE model: >>zoe=oe(zdanv,[2 2 I]).
A discrete model is returned, presented in theta format:
Discrete time IDPOLY model: y(t)=[B(q)/F(q)]u(t)+e(t) B(q)=-0.04799q^-1+0.1169q^-2 F(q)=1-1.513q^-1+0.5894^q-2
Estimated using OE from data set zdanv Loss function 0.0180258 and FPE 0.0183195 Sampling interval: 0.08.
The function bj estimates the parameters of the Box-Jenkins model: >>zbj=bj(zdanv,[2 2 2 2 1]).
A discrete model is returned, presented in theta format:
Discrete time IDPOLY model: y(t)=[B(q)F(q)]u(t)+[C(q)/D(q)]e(t) B(q)=0.01242q^-1+0.03574q^-2 C(q)=1+0.5362q^-1+0.1415q^-2
D(q)=1-0.9029q^-1+ 0.2424q^-2 F(q)=1-1.657q^-1-0.7119q^-2 Estimated using BJ from data set zdanv
Loss function 0.00672765 and FPE 0.00695098 Sampling interval: 0.08.
The function n4sid is used to estimate the parameters of state-space models in canonical form for an arbitrary number of inputs and outputs:
102
[zn4s,АO] = n4sid(z,order,ny,auxord),
where z is the matrix of experimental data; order specifies the model order. If this argument is entered as a row vector, preliminary calculations are performed for models of all specified orders (by default from the first to the tenth), with a graph output that allows choosing the optimal order. If order = 'best' (the default), the model of the best order is selected; ny is the number of outputs (by default ny = 1); auxord is the auxiliary order used by the estimation algorithm. This order must be greater than the order specified by the parameter order (by default the auxiliary order equals (1.2*order+3)). If this argument is entered as a row vector, the model of the best order is selected.
For the example under consideration, Project, we have: >>zn4s= n4sid(zdanv,[1:10],11:10]),
where the first square brackets specify the interval of model orders, order; preliminary calculations are performed for models of all specified orders from 1 to 10, with a graph output that allows choosing the optimal order. After this, the order must be typed at the MatLab command line and calculation of the model coefficients continued by pressing Enter (fig. 2.10).
Fig. 2.10. Graph for selecting the optimal model order
103
The second square brackets specify the auxiliary order used by the estimation algorithm (by default the auxiliary order equals (1.2*order+3)). In this case the optimal order is chosen without outputting the corresponding graph.
The result of executing the command is the output of the estimation procedure:
Warning: Input arguments must be scalar. >In n4sid>transf at 1027
In n4sid at 134
Select model order:('Return' gives default).
When clicking In n4sid>transf at 1027, In n4sid at 134, the M-file program editor window will appear.
When pressing Enter, "Order chosen to 3" appears. State-space model: x(t +Ts) = Ax(t) + Bu(t) + Ke(t),
y(t) = Cx(t) + Du(t) +e(t).
A =
x1
x2
x3
x1
0.96014
−0.21598
0.062944
x2
0.24873
0.66536
0.2574
x3
−0.036067 −0.64612
0.14772
B =
gas flow rate
x1
−0.00029117
x2
−0.012463
x3
−0.032588
C =
temperature
x1
x2
x3
18.093 0.078884 −0.1539
D =
gas flow rate−temperature
0
104
К = temperature
x1
0.026294
x2
−0.00982
x3
−0.032588
x(0) =
x1
−0.0073071
x2
−0.012254
x3
−0.05598
Estimated using N4SID from data set zdanv Loss function 0.0020154 and FPE 0.00211451 Sampling interval: 0.08.
The function pem estimates the parameters of a generalized multivariable linear model:
>>zpem=pem(zdanv).
State-space model:
x(t +Ts) = Ax(t) + Bu(t) + Ke(t),
y(t) = Cx(t) + Du(t) +e(t).
A =
x1
x2
x3
x1
0.9525
−0.20939
0.049221
x2
0.25496
0.65832
0.2538
x3
−0.038767
−0.64818
0.14259
B =
gas flow rate x1 −0.0002445 x2 −0.011953 x3 −0.041018
C =
temperature
x1
x2
x3
18.093
0.078583
−0.15352
105
D =
gas flow rate−temperature
0
К = temperature x1 0.05191 x2 −0.006843 x3 −0.067671 x(0) =
x1 0 x2 0 x3 0
Estimated using PEM from data set zdanv Loss function 0.00157651 and FPE 0.0016343 Sampling interval: 0.08.
2.3.3.4. Model conversion functions
For further use of the obtained models in the analysis and synthesis of systems, the System Identification Toolbox package has special functions that allow converting these models from theta format (the internal, discrete form of matrix models) into other forms, and in particular, from a discrete model into a continuous model in the form of a transfer function.
The function th2arx converts a theta-format model into an ARX model. The function has the syntax
>> [А,B]=th2arx(darx),
where darx is the theta-format model:
А=
Columns 1 through 7
1.0000 -0.9601 0.0338 -0.1035 0.0890 -0.0283 0.1383 Column 8
-0.0661
В=
106
Columns 1 through 7
0 0.0010 0.0039 0.0643 0.0625 0.0200 -0.0079 Columns 8 through 10
-0.0187 -0.0193 -0.0079
The th2ff function calculates the frequency characteristics and the corresponding standard deviations for a model in theta format. Any of the previously discussed models, for example darx, can serve as the argument of this function:
>>[g,phiv]=th2ff(darx)
IDFRD model g.
Contains Frequency Response Data for 1 output and 1 input and SpectrumData for disturbances at 1
output at 140 frequency points, ranging from 0.1 rad/s to 39.27 rad/s. Output Channels: temperature
Input Channels: gas flow rate Sampling time: 0.08
Estimated from data set zdanv using ARX. IDFRD model phiv.
Contains SpectrumData for 1 signal
at 105 frequency points, ranging from 0.1 rad/s to 39.27 rad/s. Output Channels: temperature
Sampling time: 0.08
Estimated from data set zdanv using ARX.
The th2poly function converts a model matrix in theta format into the matrices of a generalized (multivariable) linear model:
>>[A,B,C,D,K,lan,T]=th2poly(zpem)
A= 1.0000 -1.7534 1.0765 -0.2489 B= 0 0.0009 0.0012 0.0688 C= 1.0000 -0.8043 0.2732 0.1086 D=1
K= l
lan= .0016 T= 0.0800.
Here the parameter lan determines the intensity of the observation noise.
107
The th2ss function converts a theta model into a state-variable model. Any of the previously discussed models, for example darmax, can serve as the argument of this function:
>>[A,B,C,D,K,x0]=th2ss(zn4s)
A= 1.6483 1.0000 -0.7054 0 B= 0.0050 0.0456
C= l 0 D=0
K= 1.3703 -1.1387 x0= 0 0.
The th2tf function converts a theta-format model of a multivariable object into a vector of transfer functions associated with the selected input:
>>[num,den]= th2tf(zn4s)
num= 0 -0.0012 0.0069 0.0632 den= 1.0000 -1.7732 1.1013 -0.2554,
The tf command is used to represent a transfer function as a ratio
>>zzn4s=tf(num,den,0.08)
Transfer function: -0.001236z^2+0.006901z+0.06318
------------------------------------------
z^3-1.773z^2+1.101z-0.2554 Sampling time: 0.08.
The thd2thc function converts a discrete model into a continuous one. For example: to convert the discrete theta-format model zn4s (a state-variable model in canonical form for an arbitrary number of inputs and outputs) into a continuous model and to represent it as a transfer function. To do this, one must first run the thd2thc function, which converts the discrete model into a continuous one, then run the th2tf function, which converts the theta-format model of a multivariable object into a vector of transfer functions associated with the selected input, and then the tf command to represent
the transfer function as a ratio:
108
>>sn4s=thd2thc(zn4s);
>>[num,den]=th2tf(zn4s); >>sysn4s=tf(num,den).
Transfer function: 0.7165s^2-22.12s+264.4
-------------------------------------
s^3+17.06s^2+134.4s+279.1.
For the inverse conversion of a continuous model into a discrete model, there is a function thc2thd.
The th2zp function calculates the zeros, poles and static transfer coefficients (gains) of the theta-format model zn4s of a multivariable object:
>>[zepo,k]=th2zp(zn4s)
zepo =
1.000
61.0000
21.0000
81.0000
10.4622
8.7202
0.7842
0.0085
-4.8834
3.3276
0.4945+0.2848i
0.0095+0.0176i
Inf
Inf
0.4945-0.2848i
0.6312
k=
1.0000
0.9476
0.0145.
Information about the zeros and poles of the zn4s model can be obtained using the command
>>[zero,polus]=getzp(zepo)
zero = 10.4662 -4.8834
polus =0.7842 0.4945+0.2848i 0.4945 – 0.2848i.
Using the zpplot command, one can plot the zeros and poles of the zn4s model:
>>zpplot(zpform(zepo)).
Fig. 2.11 shows a plot of the zeros (marked with a circle) and poles (marked with a cross) of the zn4s model, obtained using the zpplot command.
109
Fig. 2.11. Plots of the zeros and poles of the zn4s model
The plot data show that the model is stable: the poles of the model lie inside a circle of radius 1 passing through the point with coordinates (–1; j0).
2.3.3.5. Checking the adequacy of the model
One of the important stages in the identification of automation objects is checking the quality of the model against a chosen criterion of closeness between the model output and the object output, i.e., checking its adequacy. In the System Identification Toolbox, the fit adequacy estimate is adopted as such a criterion, which is calculated using the formula
fit = norm( yh) / N ,
where norm is the norm of the vector; yh and y are the outputs of the model and the object, respectively; N is the number of elements in the data array.
To check the adequacy of the previously obtained models, we use the function
>>compare(zdane,zn4s,zpem,zoe,zbj,darx,darmax),
where zdane is the object output; zn4s, zpem, zoe, zbj, darx, darmax are the outputs of the models zn4s, zpem, zoe, zbj, darx, darmax.
The result of executing the command is a plot of the object output together with the outputs of the constructed models (Fig. 2.12).
Fig. 2.12. Plots of the object output and the model outputs
The plot shows, using colored lines, the outputs of the obtained models and the values of the adequacy criterion, expressed as percentages. The best results are shown by the models darx, zn4s and zpem.
To check the adequacy of the zn4s model, we use the function
>>compare(zdane,zn4s).
The result of executing the command is a plot of the object shown in Fig. 2.13.
Fig. 2.13. Plots of the outputs of the object and the zn4s model
111
The System Identification Toolbox also makes it possible to predict the modeling error for a given input action u(t) and a known output coordinate of the object y(t). The estimation is performed using the Predictive Error Method, abbreviated PEM, which works as follows. Suppose the model of the object under study has the form of a generalized linear model:
y(t) =W (z)u(t) +v(t),
where W (z) is the discrete transfer function of any of the previously discussed models. Here the noise v(t) can be represented as v(t) = H (z)e(z),
where e(z) is discrete white noise, which in fact characterizes the model error; H (z) is a certain polynomial in z that converts discrete white noise into the actual disturbances present when measuring the output parameters of the object.
From these expressions it follows that
e(t) = H −1(z)[ y(t) −W (z)u(t)].
The resid function calculates the residual error e for a given model, as well as r – the matrix of values of the autocorrelation function of the process e(t) and the values of the cross-correlation function between the residual errors e(t) and the outputs of the automation object y(t), together with the corresponding 99% confidence intervals.
In addition to the stated values, plots of these functions are also output. As an example, let us compare the residual errors and the corresponding correlation functions for the previously obtained models darx and zbj, which have the maximum and minimum adequacy estimates, using the commands
>>[e,r]=resid(zdan,darx); >>[e1,r1]=resid(zdan,zbj).
The plots shown (Fig. 2.14, a, b) characterize a uniform distribution of residual errors over the entire range of variation of the time intervals τ, with the residual error values for the darx model being almost twice as large as for the zbj model. To display the plots, the command resid(r) must be executed.
112
a
b
Fig. 2.16. Plots of the autocorrelation and cross-correlation functions for the model: a – zbj; b – darx
113
After executing the function
[e,r]=resid(zdan,darx)
MatLab returns:
Time domain data set with 1000 samples. Sampling interval: 0.08
Outputs Unit (if specified)
e@temperature deg.C 100 Inputs Unit (if specified)
u1
r= 1.0e+003 *.
After executing the command >>resid(r), the plot of the autocorrelation and cross-correlation functions for the model is displayed.
Thus, in the course of evaluating the adequacy of the various object models, the models darx, zn4s and zpem were identified as having the maximum adequacy criterion values and, consequently, can be used further in the analysis and synthesis of systems.
2.3.3.6. Analysis of the model of the technical object
For the analysis of the system model, the zn4s model is considered, which has one of the best adequacy indicators. Previously, various forms of this model were obtained:
• zn4s – discrete model in theta format (LTI format):
A =
x1
x2
x3
x1
0.96014
−0.21598
0.062944
x2
0.24873
0.66536
0.2574
x3
−0.036067
−0.64612
0.14772
B =
gas flow rate x1 −0.00029117 x2 −0.012463 x3 −0.032588
114
C =
temperature
x1
x2
x3
18.093
0.078884
−0.1539
D =
gas flow rate−temperature
0
K = temperature
x1
0.026294
x2
−0.00982
x3
−0.032588
x(0) =
x1
−0.0073071
x2
−0.012254
x3
−0.05598
Estimated using N4SID from data set zdanv Loss function 0.0020154 and FPE 0.00211451 Sampling interval: 0.08.
• sn4s – continuous model in theta format (LTI format) Statespace model:
dx(t) = Ax(t) + Bu(t) + Ke(t)
dt
y(t) = Cx(t) + Du(t) +e(t)
A =
x1
x2
x3
x1
−0.24803
−2.0709
2.222
x2
3.2929
−1.3184
6.3825
x3
1.8425
−16.783
−15.496
115
B =
gas flow rate x1 0.032775 x2 −0.028698 x3 −0.81694
C =
temperature
x1
x2
x3
18.093
0.078583
−0.1539
D =
gas flow rate−temperature
0
K = temperature x1 0.43342 x2 0.13241
x3 −1.7679 x(0) =
x1 −0.0073071 x2 −0.012254 x3 −0.05598
Estimated using N4SID (later modified)
Loss function 0.00213802 and FPE 0.0022164
As can be seen in the presented models, the values of the coefficients of the matrices A, B, C, D, K are different. This is explained by the fact that for the continuous model a Z-transform was performed in order to obtain the discrete model;
• zzn4s – discrete model in the form of a transfer function
W (z) =
−0.001236z2
+0.006901z +0.06318
;
z3 −1.773z2
+1.101z −0.2554
116
• sysn4s – continuous model in the form of a transfer function
W (s) =
0.7165s2
+ 22.12s + 264.4
.
s3 +17.06s2
+134.4s +279.1
The presented forms are one and the same model, written in different forms and formats.
The analysis of the system is carried out on the basis of studying its dynamic and frequency characteristics.
The step function calculates and plots the step response of the system:
>>step(sys)
>>step(sys,t) >>step(sys1,sys2,..,sysN,t)
>>step(sys1,'PlotStyle1',…,sysN, 'PlotStyleN') >>[y,t,x]=step(sys).
For models given in state space, the initial conditions are taken to be zero. The arguments of the function are as follows:
•sys,sys1,sys2,…,sysN – names of the models for which the step responses are constructed;
•t – argument specifying the moment at which the simulation ends, either
in the form t =Tfinal (in seconds), or in the form t = 0:dt: Tfinal. For discrete models, the value of dt must be equal to the sampling interval; for continuous models it must be sufficiently small to account for the fastest changes in the transient process;
•'PlotStyle1',..,'PlotStyleN' – string variables specifying the styles (line types) when displaying several plots at once.
Return values.
•plots of the transient processes;
•y, x, t – respectively, vectors containing the values of the transient process, the state variables and the time instants (when these values are returned, the transient process plot is not displayed).
Plots of the step response of the system, represented by the discrete model zzn4s and the continuous model sysn4s, constructed using the step function:
>>step(zzn4s,sysn4s).
117
After executing the step command, MatLab returns plots of the transient process (Fig. 2.15). Clicking the left mouse button anywhere on the transient process plot causes a pop-up information tooltip to appear, showing the current numerical value of the transient process and the corresponding time instant.
Clicking the right mouse button anywhere on the transient process plot causes a pop-up menu to appear for editing the pop-up information tooltip window.
Fig. 2.15. Plots of the transient processes of the models zzn4s and sysn4s
On the transient-process plots, the step-shaped line represents the transient process of the discrete model, while the solid line represents that of the continuous model. In addition, the plot area shows the main characteristics of the transient process:
•settling time – 1.9 s for both models;
•steady-state value of the output coordinate 0.948 for both models.
To plot the impulse response of the models, use the command
>>impulse(zzn4s,sysn4s).
After executing the impulse command, MatLab returns the plots
(Fig. 2.16).
118
Fig. 2.16. Plots of the impulse responses of the models zzn4s and sysn4s
The main characteristics of the system model when a unit impulse action is applied at the input are:
•peak amplitude – for the discrete model it is 0.14, and for the continuous model – 1.77;
•settling time for both the discrete and continuous models is 1.85 s.
To determine the static gain of the system model, one can use the degain command:
>>k=degain(sysn4s); >>k=degain(zzn4s).
For the continuous and discrete system, the gain equals 0.9476.
To determine the frequency response of the models, the bode command is used. The frequency responses of the discrete model zzn4s and the continuous model sysn4s are shown in Fig. 2.17.
On the frequency-response plots, the values of the gain margin (Gain Margin) are indicated, which for the discrete model amount to 10.9 dB, and for the continuous model – 9.34 dB.
119
Fig. 2.17. Frequency responses of the models
The stability margin values can also be determined in MatLab command-line mode using the commands
>>[Gm,Pm,Wcg,Wcp]=margin(sysn4s); >>[Gm1,Pm1,Wcg1,Wcp1]=margin(zzn4s).
For the continuous model: Gm=3.3642
Pm=Inf
Wcg=7.7471
Wcp=NaN.
For the discrete model: Gm1=2.9158
Pm1=Inf
Wcg1=6.8677
Wcp1=NaN,
where Gm is the gain margin in natural units at frequency Wcg, Pm is the phase margin at frequency Wcp.
To determine the stability margins on a logarithmic scale, the following operations must be performed:
>>Gmlog=20*log10(Gm1) – for the discrete model:
Gmlog=9.2951, >>Gmlog=20*log10(Gm) – for the continuous model:
Gmlog = 10.5377.
As can be seen, determining the stability margins by this last method allows these values to be calculated more precisely than from the frequency-response plots. Analysis of the frequency responses shows that the models zzn4s and sysn4s are stable, with the corresponding gain margins. The phase margin is infinite.
This conclusion is also confirmed by the complex amplitude-phase characteristic APC (Fig. 2.18), since the APC locus does not cross the point of the complex plane with coordinates (–1, j0).
To construct the APC, use the command
>>nyquist(zzn4s,sysn4s).
The stability of the models can also be determined using the zero-pole map, based on the position of the model's zeros relative to the unit-radius circle in the complex plane, as shown in Fig. 2.18.
Fig. 2.18. APC locus for the continuous and discrete models
121
The zero-pole map of the models can also be constructed using the
commands pzmap(zzn4s,sysn4s) or pzmap(zn4s,sn4s).
It should be recalled that the output parameter of the automation object (in our case, temperature) y(t) consists of two components: the output of the object y(t) and additive noise e(t), caused by the influence of external factors on the automation object. The latter component is a purely random quantity and is characterized by statistical parameters. In the course of identification, dependencies were obtained that determine the component e(t) and the correlation dependencies between e(t) and the output of the object y(t):
e(t) = H1(z)[ y(t) −W (z)u(t)].
Using the plot(e) command, a plot of the change in e(t) is constructed
and the main statistical characteristics of the noise are determined
(Fig. 2.19).
Fig. 2.19. Plot of the additive model e(t)
To obtain the statistical characteristics, in the plot's menu bar under Tools, select the Data statistics option. As a result of executing this command, a window will appear showing the main statistical characteristics of the random process of change of e(t) over time (Fig. 2.20):
•min and max – the minimum and maximum values of the noise;
•mean – the arithmetic mean value;
•median – the median of the process;
•std – the standard deviation;
•range
продолжение следует...
Часть 1 2.3. Software Tools for System Identification and Modeling
Часть 2 - 2.3. Software Tools for System Identification and Modeling
Comments