Documentation |
Estimate and validate linear models from multiple-input/single-output (MISO) data to find the one that best describes the system dynamics.
After completing this tutorial, you will be able to accomplish the following tasks using the command line:
Create data objects to represent data.
Plot the data.
Process data by removing offsets from the input and output signals.
Estimate and validate linear models from the data.
Simulate and predict model output.
Note: This tutorial uses time-domain data to demonstrate how you can estimate linear models. The same workflow applies to fitting frequency-domain data. |
This tutorial uses the data file co2data.mat, which contains two experiments of two-input and single-output (MISO) time-domain data from a steady-state that the operator perturbed from equilibrium values.
In the first experiment, the operator introduced a pulse wave to both inputs. In the second experiment, the operator introduced a pulse wave to the first input and a step signal to the second input.
Load the data.
load co2data;
This command loads the following five variables into the MATLAB Workspace:
Input_exp1 and Output_exp1 are the input and output data from the first experiment, respectively.
Input_exp2 and Output_exp2 are the input and output data from the second experiment, respectively.
Time is the time vector from 0 to 1000 minutes, increasing in equal increments of 0.5 min.
For both experiments, the input data consists of two columns of values. The first column of values is the rate of chemical consumption (in kilograms per minute), and the second column of values is the current (in amperes). The output data is a single column of the rate of carbon-dioxide production (in milligrams per minute).
Plot the input and output data from both experiments.
plot(Time,Input_exp1,Time,Output_exp1) legend('Input 1','Input 2','Output 1') figure plot(Time,Input_exp2,Time,Output_exp2) legend('Input 1','Input 2','Output 1')
The first plot shows the first experiment, where the operator applies a pulse wave to each input to perturb it from its steady-state equilibrium.
The second plot shows the second experiment, where the operator applies a pulse wave to the first input and a step signal to the second input.
Plotting the data, as described in Plotting the Input/Output Data, shows that the inputs and the outputs have nonzero equilibrium values. In this portion of the tutorial, you subtract equilibrium values from the data.
The reason you subtract the mean values from each signal is because, typically, you build linear models that describe the responses for deviations from a physical equilibrium. With steady-state data, it is reasonable to assume that the mean levels of the signals correspond to such an equilibrium. Thus, you can seek models around zero without modeling the absolute equilibrium levels in physical units.
Zoom in on the plots to see that the earliest moment when the operator applies a disturbance to the inputs occurs after 25 minutes of steady-state conditions (or after the first 50 samples). Thus, the average value of the first 50 samples represents the equilibrium conditions.
Use the following commands to remove the equilibrium values from inputs and outputs in both experiments:
Input_exp1 = Input_exp1-... ones(size(Input_exp1,1),1)*mean(Input_exp1(1:50,:)); Output_exp1 = Output_exp1-... mean(Output_exp1(1:50,:)); Input_exp2 = Input_exp2-... ones(size(Input_exp2,1),1)*mean(Input_exp2(1:50,:)); Output_exp2 = Output_exp2-... mean(Output_exp2(1:50,:));
The System Identification Toolbox™ data objects, iddata and idfrd, encapsulate data values and data properties into a single entity. You can use the System Identification Toolbox commands to conveniently manipulate these data objects as single entities.
In this portion of the tutorial, you create two iddata objects, one for each of the two experiments. You use the data from Experiment 1 for model estimation, and the data from Experiment 2 for model validation. You work with two independent data sets because you use one data set for model estimation and the other for model validation.
Note: When two independent data sets are not available, you can split one data set into two parts, assuming that each part contains enough information to adequately represent the system dynamics. |
You must have already loaded the sample data into the MATLAB^{®} workspace, as described in Loading Data into the MATLAB Workspace.
Use these commands to create two data objects, ze and zv :
Ts = 0.5; % Sampling interval is 0.5 min
ze = iddata(Output_exp1,Input_exp1,Ts);
zv = iddata(Output_exp2,Input_exp2,Ts);
ze contains data from Experiment 1 and zv contains data from Experiment 2. Ts is the sampling interval.
The iddata constructor requires three arguments for time-domain data and has the following syntax:
data_obj = iddata(output,input,sampling_interval);
To view the properties of an iddata object, use the get command. For example, type this command to get the properties of the estimation data:
get(ze)
ans = Domain: 'Time' Name: '' OutputData: [2001x1 double] y: 'Same as OutputData' OutputName: {'y1'} OutputUnit: {''} InputData: [2001x2 double] u: 'Same as InputData' InputName: {2x1 cell} InputUnit: {2x1 cell} Period: [2x1 double] InterSample: {2x1 cell} Ts: 0.5000 Tstart: [] SamplingInstants: [2001x0 double] TimeUnit: 'seconds' ExperimentName: 'Exp1' Notes: {} UserData: []
To learn more about the properties of this data object, see the iddata reference page.
To modify data properties, you can use dot notation or the set command. For example, to assign channel names and units that label plot axes, type the following syntax in the MATLAB Command Window:
% Set time units to minutes ze.TimeUnit = 'min'; % Set names of input channels ze.InputName = {'ConsumptionRate','Current'}; % Set units for input variables ze.InputUnit = {'kg/min','A'}; % Set name of output channel ze.OutputName = 'Production'; % Set unit of output channel ze.OutputUnit = 'mg/min'; % Set validation data properties zv.TimeUnit = 'min'; zv.InputName = {'ConsumptionRate','Current'}; zv.InputUnit = {'kg/min','A'}; zv.OutputName = 'Production'; zv.OutputUnit = 'mg/min';
You can verify that the InputName property of ze is changed, or "index" into this property:
ze.inputname;
Tip Property names, such as InputUnit, are not case sensitive. You can also abbreviate property names that start with Input or Output by substituting u for Input and y for Output in the property name. For example, OutputUnit is equivalent to yunit. |
You can plot iddata objects using the plot command.
Plot the estimation data.
plot(ze)
The bottom axes show inputs ConsumptionRate and Current, and the top axes show the output ProductionRate .
Plot the validation data in a new figure window.
figure % Open a new MATLAB Figure window plot(zv) % Plot the validation data
Zoom in on the plots to see that the experiment process amplifies the first input ( ConsumptionRate ) by a factor of 2, and amplifies the second input (Current ) by a factor of 10.
Before you begin, create a new data set that contains only the first 1000 samples of the original estimation and validation data sets to speed up the calculations.
Ze1 = ze(1:1000); Zv1 = zv(1:1000);
For more information about indexing into iddata objects, see the corresponding reference page.
Frequency-response and step-response are nonparametric models that can help you understand the dynamic characteristics of your system. These models are not represented by a compact mathematical formula with adjustable parameters. Instead, they consist of data tables.
In this portion of the tutorial, you estimate these models using the data set ze. You must have already created ze, as described in Creating iddata Objects.
The response plots from these models show the following characteristics of the system:
The response from the first input to the output might be a second-order function.
The response from the second input to the output might be a first-order or an overdamped function.
The System Identification Toolbox product provides three functions for estimating the frequency response:
etfe computes the empirical transfer function using Fourier analysis.
spa estimates the transfer function using spectral analysis for a fixed frequency resolution.
spafdr lets you specify a variable frequency resolution for estimating the frequency response.
Use spa to estimate the frequency response.
Ge=spa(ze);
Plot the frequency response as a Bode plot.
bode(Ge)
The amplitude peaks at the frequency of about 0.7 rad/s, which suggests a possible resonant behavior (complex poles) for the first input-to-output combination - ConsumptionRate to Production .
In both plots, the phase rolls off rapidly, which suggests a time delay for both input/output combinations.
To estimate the step response from the data, first estimate a non-parametric impulse response model (FIR filter) from data and then plot its step response.
% model estimation Mimp = impulseest(Ze1,60); % step response step(Mimp)
The step response for the first input/output combination suggests an overshoot, which indicates the presence of an underdamped mode (complex poles) in the physical system.
The step response from the second input to the output shows no overshoot, which indicates either a first-order response or a higher-order response with real poles (overdamped response).
The step-response plot indicates a nonzero delay in the system, which is consistent with the rapid phase roll-off you got in the Bode plot you created in Estimating the Empirical Step Response.
To identify parametric black-box models, you must specify the input/output delay as part of the model order.
If you do not know the input/output delays for your system from the experiment, you can use the System Identification Toolbox software to estimate the delay.
In the case of single-input systems, you can read the delay on the impulse-response plot. However, in the case of multiple-input systems, such as the one in this tutorial, you might be unable to tell which input caused the initial change in the output and you can use the delayest command instead.
The command estimates the time delay in a dynamic system by estimating a low-order, discrete-time ARX model with a range of delays, and then choosing the delay that corresponding to the best fit.
The ARX model structure is one of the simplest black-box parametric structures. In discrete-time, the ARX structure is a difference equation with the following form:
$$\begin{array}{l}y(t)+{a}_{1}y(t-1)+\dots +{a}_{na}y(t-{n}_{a})=\\ {\text{b}}_{\text{1}}u(t-{n}_{k})+\dots +{b}_{nb}u(t-{n}_{k}-{n}_{b}+1)+e(t)\end{array}$$
y(t) represents the output at time t, u(t) represents the input at time t, n_{a} is the number of poles, n_{b} is the number of b parameters (equal to the number of zeros plus 1), n_{k} is the number of samples before the input affects output of the system (called the delay or dead time of the model), and e(t) is the white-noise disturbance.
delayest assumes that n_{a}=n_{b}=2 and that the noise e is white or insignificant, and estimates n_{k}.
To estimate the delay in this system, type:
delayest(ze)
ans = 5 10
This result includes two numbers because there are two inputs: the estimated delay for the first input is 5 data samples, and the estimated delay for the second input is 10 data samples. Because the sampling interval for the experiments is 0.5 min, this corresponds to a 2.5 -min delay before the first input affects the output, and a 5.0 -min delay before the second input affects the output.
There are two alternative methods for estimating the time delay in the system:
Plot the time plot of the input and output data and read the time difference between the first change in the input and the first change in the output. This method is practical only for single-input/single-output system; in the case of multiple-input systems, you might be unable to tell which input caused the initial change in the output.
Plot the impulse response of the data with a 1-standard-deviation confidence region. You can estimate the time delay using the time when the impulse response is first outside the confidence region.
Model order is one or more integers that define the complexity of the model. In general, model order is related to the number of poles, the number of zeros, and the response delay (time in terms of the number of samples before the output responds to the input). The specific meaning of model order depends on the model structure.
To compute parametric black-box models, you must provide the model order as an input. If you do not know the order of your system, you can estimate it.
After completing the steps in this section, you get the following results:
For the first input/output combination: n_{a}=2, n_{b}=2, and the delay n_{k}=5.
For the second input/output combination: n_{a}=1, n_{b}=1, and the delay n_{k}=10.
Later, you explore different model structures by specifying model-order values that are slight variations around these initial estimate.
In this portion of the tutorial, you use struc, arxstruc, and selstruc to estimate and compare simple polynomial (ARX) models for a range of model orders and delays, and select the best orders based on the quality of the model.
The following list describes the results of using each command:
struc creates a matrix of model-order combinations for a specified range of n_{a}, n_{b}, and n_{k} values.
arxstruc takes the output from struc, systematically estimates an ARX model for each model order, and compares the model output to the measured output. arxstruc returns the loss function for each model, which is the normalized sum of squared prediction errors.
selstruc takes the output from arxstruc and opens the ARX Model Structure Selection window, which resembles the following figure, to help you choose the model order.
You use this plot to select the best-fit model.
The horizontal axis is the total number of parameters — n_{a} + n_{b}.
The vertical axis, called Unexplained output variance (in %), is the portion of the output not explained by the model—the ARX model prediction error for the number of parameters shown on the horizontal axis.
The prediction error is the sum of the squares of the differences between the validation data output and the model one-step-ahead predicted output.
n_{k} is the delay.
Three rectangles are highlighted on the plot in green, blue, and red. Each color indicates a type of best-fit criterion, as follows:
Red — Best fit minimizes the sum of the squares of the difference between the validation data output and the model output. This rectangle indicates the overall best fit.
Green — Best fit minimizes Rissanen MDL criterion.
Blue — Best fit minimizes Akaike AIC criterion.
In this tutorial, the Unexplained output variance (in %) value remains approximately constant for the combined number of parameters from 4 to 20. Such constancy indicates that model performance does not improve at higher orders. Thus, low-order models might fit the data equally well.
Note: When you use the same data set for estimation and validation, use the MDL and AIC criteria to select model orders. These criteria compensate for overfitting that results from using too many parameters. For more information about these criteria, see the selstruc reference page. |
In this tutorial, there are two inputs to the system and one output and you estimate model orders for each input/output combination independently. You can either estimate the delays from the two inputs simultaneously or one input at a time.
It makes sense to try the following order combinations for the first input/output combination:
n_{a}=2:5
n_{b}=1:5
n_{k}=5
This is because the nonparametric models you estimated in Estimating Impulse Response Models show that the response for the first input/output combination might have a second-order response. Similarly, in Estimating Delays in the Multiple-Input System, the delay for this input/output combination was estimated to be 5.
To estimate model order for the first input/output combination:
Use struc to create a matrix of possible model orders.
NN1 = struc(2:5,1:5,5);
Use selstruc to compute the loss functions for the ARX models with the orders in NN1.
selstruc(arxstruc(ze(:,:,1),zv(:,:,1),NN1))
This command opens the interactive ARX Model Structure Selection window.
Note: The Rissanen MDL and Akaike AIC criteria produces equivalent results and are both indicated by a blue rectangle on the plot. |
The red rectangle represents the best overall fit, which occurs for n_{a}=2, n_{b}=3, and n_{k}=5. The height difference between the red and blue rectangles is insignificant. Therefore, you can choose the parameter combination that corresponds to the lowest model order and the simplest model.
Click the blue rectangle, and then click Select to choose that combination of orders:
n_{a}=2
n_{b}=2
n_{k}=5
To continue, press any key while in the MATLAB Command Window.
It makes sense to try the following order combinations for the second input/output combination:
n_{a}=1:3
n_{b}=1:3
n_{k}=10
This is because the nonparametric models you estimated in Estimating Impulse Response Models show that the response for the second input/output combination might have a first-order response. Similarly, in Estimating Delays in the Multiple-Input System, the delay for this input/output combination was estimated to be 10.
To estimate model order for the second input/output combination:
Use struc to create a matrix of possible model orders.
NN2 = struc(1:3,1:3,10);
Use selstruc to compute the loss functions for the ARX models with the orders in NN2.
selstruc(arxstruc(ze(:,:,2),zv(:,:,2),NN2))
This command opens the interactive ARX Model Structure Selection window.
Note: The Akaike AIC and the overall best fit criteria produces equivalent results. Both are indicated by the same red rectangle on the plot. |
The height difference between all of the rectangles is insignificant and all of these model orders result in similar model performance. Therefore, you can choose the parameter combination that corresponds to the lowest model order and the simplest model.
Click the yellow rectangle on the far left, and then click Select to choose the lowest order: n_{a}=1, n_{b}=1, and n_{k}=10.
To continue, press any key while in the MATLAB Command Window.
In this portion of the tutorial, you estimate a continuous-time transfer function. You must have already prepared your data, as described in Preparing Data. You can use the following results of estimated model orders to specify the orders of the model:
For the first input/output combination, use:
Two poles, corresponding to na=2 in the ARX model.
Delay of 5, corresponding to nk=5 samples (or 2.5 minutes) in the ARX model.
For the second input/output combination, use:
One pole, corresponding to na=1 in the ARX model
Delay of 10, corresponding to nk=10 samples (or 5 minutes) in the ARX model.
You can estimate a transfer function of these orders using the tfest command. For the estimation, you can also choose to view a progress report by setting the Display option to on in the option set created by the tfestOptions command.
Opt = tfestOptions('Display', 'on');
Collect the model orders and delays into variables to pass to tfest.
np = [2 1]; ioDelay = [2.5 5];
Estimate the transfer function.
mtf = tfest(Ze1, np, [], ioDelay, Opt);
View the model's coefficients.
mtf
mtf = From input "ConsumptionRate" to output "Production": 91.75 s - 1.272 exp(-2.5*s) * --------------------- s^2 + 37.76 s + 1.009 From input "Current" to output "Production": 5.186 exp(-5*s) * ---------- s + 0.5066 Continuous-time identified transfer function. Parameterization: Number of poles: [2 1] Number of zeros: [1 0] Number of free coefficients: 6 Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using TFEST on time domain data "Ze1". Fit to estimation data: 85.89% (simulation focus) FPE: 6.424, MSE: 6.259
The model's display shows more than 85% fit to estimation data.
In this portion of the tutorial, you create a plot that compares the actual output and the model output using the compare command.
compare(Zv1, mtf)
The comparison shows about 60% fit.
Use the resid command to evaluate the characteristics of the residuals.
resid(Zv1, mtf)
The residuals show high degree of auto-correlation. This is not unexpected since the model mtf does not have any components to describe the nature of the noise separately. To model both the measured input-output dynamics and the disturbance dynamics you will need to use a model structure that contains elements to describe the noise component. You can use bj, ssest and procest commands, which create models of polynomial, state-space and process structures respectively. These structures, among others, contain elements to capture the noise behavior.
In this portion of the tutorial, you estimate a low-order, continuous-time transfer function (process model). the System Identification Toolbox product supports continuous-time models with at most three poles (which might contain underdamped poles), one zero, a delay element, and an integrator.
You must have already prepared your data, as described in Preparing Data.
You can use the following results of estimated model orders to specify the orders of the model:
For the first input/output combination, use:
Two poles, corresponding to n_{a}=2 in the ARX model.
Delay of 5, corresponding to n_{k}=5 samples (or 2.5 minutes) in the ARX model.
For the second input/output combination, use:
One pole, corresponding to n_{a}=1 in the ARX model.
Delay of 10, corresponding to n_{k}=10 samples (or 5 minutes) in the ARX model.
Note: Because there is no relationship between the number of zeros estimated by the discrete-time ARX model and its continuous-time counterpart, you do not have an estimate for the number of zeros. In this tutorial, you can specify one zero for the first input/output combination, and no zero for the second-output combination. |
Use the idproc command to create two model structures, one for each input/output combination:
midproc0 = idproc({'P2ZUD','P1D'}, 'TimeUnit', 'minutes');
idproc accepts a cell array that contains two strings that specify the model structure for each input/output combination:
'P2ZUD' represents a transfer function with two poles ( P2 ), one zero ( Z ), underdamped (complex-conjugate) poles ( U ) and a delay ( D ).
'P1D' represents a transfer function with one pole ( P1 ) and a delay ( D ).
The example also uses the TimeUnit parameter to specify the unit of time used.
View the two resulting models.
midproc0
midproc0 = Process model with 2 inputs: y = G11(s)u1 + G12(s)u2 From input 1 to output 1: 1+Tz*s G11(s) = Kp * ---------------------- * exp(-Td*s) 1+2*Zeta*Tw*s+(Tw*s)^2 Kp = NaN Tw = NaN Zeta = NaN Td = NaN Tz = NaN From input 2 to output 1: Kp G12(s) = ---------- * exp(-Td*s) 1+Tp1*s Kp = NaN Tp1 = NaN Td = NaN Parameterization: 'P2DUZ' 'P1D' Number of free coefficients: 8 Use "getpvec", "getcov" for parameters and their uncertainties. Status: Created by direct construction or transformation. Not estimated.
The parameter values are set to NaN because they are not yet estimated.
Set the time delay property of the model object to 2.5 min and 5 min for each input/output pair as initial guesses. Also, set an upper limit on the delay because good initial guesses are available.
midproc0.Structure(1,1).Td.Value=2.5; midproc0.Structure(1,2).Td.Value=5; midproc0.Structure(1,1).Td.Maximum=3; midproc0.Structure(1,2).Td.Maximum=7;
Note: When setting the delay constraints, you must specify the delays in terms of actual time units (minutes, in this case) and not the number of samples. |
procest is an iterative estimator of process models, which means that it uses an iterative nonlinear least-squares algorithm to minimize a cost function. The cost function is the weighted sum of the squares of the errors.
Depending on its arguments, procest estimates different black-box polynomial models. You can use procest, for example, to estimate parameters for linear continuous-time transfer-function, state-space, or polynomial model structures. To use procest, you must provide a model structure with unknown parameters and the estimation data as input arguments.
For this portion of the tutorial, you must have already defined the model structure, as described in Specifying the Structure of the Process Model. Use midproc0 as the model structure and Ze1 as the estimation data:
midproc = procest(Ze1, midproc0); present(midproc)
midproc = Process model with 2 inputs: y = G11(s)u1 + G12(s)u2 From input "ConsumptionRate" to output "Production": 1+Tz*s G11(s) = Kp * ---------------------- * exp(-Td*s) 1+2*Zeta*Tw*s+(Tw*s)^2 Kp = -1.1804 +/- 0.29989 Tw = 1.6566 +/- 629.69 Zeta = 15.917 +/- 6039.2 Td = 2.425 +/- 56.972 Tz = -109.2 +/- 57.374 From input "Current" to output "Production": Kp G12(s) = ---------- * exp(-Td*s) 1+Tp1*s Kp = 10.264 +/- 0.048403 Tp1 = 2.0488 +/- 0.054899 Td = 4.9175 +/- 0.034373 Parameterization: 'P2DUZ' 'P1D' Number of free coefficients: 8 Use "getpvec", "getcov" for parameters and their uncertainties. Status: Termination condition: Maximum number of iterations reached. Number of iterations: 20, Number of function evaluations: 276 Estimated using PROCEST on time domain data "Ze1". Fit to estimation data: 86.2% FPE: 6.08, MSE: 5.984 More information in model's "Report" property.
Unlike discrete-time polynomial models, continuous-time models let you estimate the delays. In this case, the estimated delay values are different from the initial guesses you specified of 2.5 and 5 , respectively. The large uncertainties in the estimated values of the parameters of G_1(s) indicate that the dynamics from input 1 to the output are not captured well.
To learn more about estimating models, see Process Models.
In this portion of the tutorial, you create a plot that compares the actual output and the model output.
compare(Zv1,midproc)
The plot shows that the model output reasonably agrees with the measured output: there is an agreement of 60% between the model and the validation data.
Use resid to perform residual analysis:
resid(Zv1,midproc)
Because the sample system has two inputs, there are two cross-correlation plots of the residuals with each input, as shown in the following figure.
Autocorrelation and Cross-Correlations of Residuals with the First Input
The cross-correlation between the second input and residual errors is significant.
After MATLAB displays the first plot, press Enter to view the cross-correlation with the second input, as shown in the following figure.
Cross-Correlations of Residuals with the Second Input
In the preceding figure, the autocorrelation plot shows values outside the confidence region and indicates that the residuals are correlated.
Change the algorithm for iterative parameter estimation to Levenberg-Marquardt.
Opt = procestOptions;
Opt.SearchMethod = 'lm';
midproc1 = procest(Ze1, midproc, Opt);
Tweaking the algorithm properties or specifying initial parameter guesses and rerunning the estimation may improve the simulation results. Adding a noise model may improve prediction results but not necessarily the simulation results.
This portion of the tutorial shows how to estimate a process model and include a noise model in the estimation. Including a noise model typically improves model prediction results but not necessarily the simulation results.
Use the following commands to specify a first-order ARMA noise:
Opt = procestOptions;
Opt.DisturbanceModel = 'ARMA1';
midproc2 = procest(Ze1, midproc0, Opt);
You can type 'dist' instead of 'DisturbanceModel'. Property names are not case sensitive, and you only need to include the portion of the name that uniquely identifies the property.
Compare the resulting model to the measured data.
compare(Zv1,midproc2)
The plot shows that the model output maintains reasonable agreement with the validation-data output.
Perform residual analysis.
figure resid(Zv1,midproc2)
Press Enter to view the cross-correlation of the residuals with the second input.
The next plot shows that adding a noise model produces uncorrelated residuals: the top set of axes show that the autocorrelation values are inside the confidence bounds. This indicates a more accurate model.
In this portion of the tutorial, you estimate several different types of black-box, input-output polynomial models.
You must have already prepared your data, as described in Preparing Data.
You can use the following previous results of estimated model orders to specify the orders of the polynomial model:
For the first input/output combination, use:
Two poles, corresponding to n_{a}=2 in the ARX model.
One zero, corresponding to n_{b}=2 in the ARX model.
Delay of 5, corresponding to n_{k}=5 samples (or 2.5 minutes) in the ARX model.
For the second input/output combination, use:
One pole, corresponding to n_{a}=1 in the ARX model.
No zeros, corresponding to n_{b}=1 in the ARX model.
Delay of 10, corresponding to n_{k}=10 samples (or 5 minutes) in the ARX model.
About ARX Models. For a single-input/single-output system (SISO), the ARX model structure is:
$$\begin{array}{l}y(t)+{a}_{1}y(t-1)+\dots +{a}_{na}y(t-{n}_{a})=\\ {\text{b}}_{\text{1}}u(t-{n}_{k})+\dots +{b}_{nb}u(t-{n}_{k}-{n}_{b}+1)+e(t)\end{array}$$
y(t) represents the output at time t, u(t) represents the input at time t, n_{a} is the number of poles, n_{b} is the number of zeros plus 1, n_{k} is the number of samples before the input affects output of the system (called the delay or dead time of the model), and e(t) is the white-noise disturbance.
The ARX model structure does not distinguish between the poles for individual input/output paths: dividing the ARX equation by A, which contains the poles, shows that A appears in the denominator for both inputs. Therefore, you can set n_{a} to the sum of the poles for each input/output pair, which is equal to 3 in this case.
The System Identification Toolbox product estimates the parameters $${a}_{1}\dots {a}_{n}$$ and $${b}_{1}\dots {b}_{n}$$ using the data and the model orders you specify.
Estimating ARX Models Using arx. Use arx to compute the polynomial coefficients using the fast, noniterative method arx:
marx = arx(Ze1,'na',3,'nb',[2 1],'nk',[5 10]); present(marx) % Displays model parameters % with uncertainty information
marx = Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t) A(z) = 1 - 1.027 (+/- 0.02917) z^-1 + 0.1675 (+/- 0.04214) z^-2 + 0.01307 (+/- 0.02591) z^-3 B1(z) = 1.86 (+/- 0.1896) z^-5 - 1.608 (+/- 0.1894) z^-6 B2(z) = 1.612 (+/- 0.07417) z^-10 Sample time: 0.5 minutes Parameterization: Polynomial orders: na=3 nb=[2 1] nk=[5 10] Number of free coefficients: 6 Use "polydata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using ARX on time domain data "Ze1". Fit to estimation data: 90.7% (prediction focus) FPE: 2.75, MSE: 2.719 More information in model's "Report" property.
MATLAB estimates the polynomials A , B1 , and B2. The uncertainty for each of the model parameters is computed to 1 standard deviation and appears in parentheses next to each parameter value.
Alternatively, you can use the following shorthand syntax and specify model orders as a single vector:
marx = arx(Ze1,[3 2 1 5 10]);
Accessing Model Data. The model you estimated, marx, is a discrete-time idpoly object. To get the properties of this model object, you can use the get function:
get(marx)
a: [1 -1.0266 0.1675 0.0131] b: {[0 0 0 0 0 1.8600 -1.6085] [0 0 0 0 0 0 0 0 0 0 1.6117]} c: 1 d: 1 f: {[1] [1]} IntegrateNoise: 0 Variable: 'z^-1' ioDelay: [0 0] Structure: [1x1 pmodel.polynomial] NoiseVariance: 2.7611 Report: [1x1 idresults.arx] InputDelay: [2x1 double] OutputDelay: 0 Ts: 0.5000 TimeUnit: 'minutes' InputName: {2x1 cell} InputUnit: {2x1 cell} InputGroup: [1x1 struct] OutputName: {'Production'} OutputUnit: {'mg/min'} OutputGroup: [1x1 struct] Name: '' Notes: {} UserData: [] SamplingGrid: [1x1 struct]
You can access the information stored by these properties using dot notation. For example, you can compute the discrete poles of the model by finding the roots of the A polynomial.
marx_poles=roots(marx.a)
marx_poles = 0.7953 0.2883 -0.0570
In this case, you access the A polynomial using marx.a.
The model marx describes system dynamics using three discrete poles.
Tip You can also use pole to compute the poles of a model directly. |
Learn More. To learn more about estimating polynomial models, see Input-Output Polynomial Models.
For more information about accessing model data, see Data Extraction.
About State-Space Models. The general state-space model structure is:
$$\begin{array}{l}\frac{dx}{dt}=Ax(t)+Bu(t)+Ke(t)\\ y(t)=Cx(t)+Du(t)+e(t)\end{array}$$
y(t) represents the output at time t, u(t) represents the input at time t, x(t) is the state vector at time t, and e(t) is the white-noise disturbance.
You must specify a single integer as the model order (dimension of the state vector) to estimate a state-space model. By default, the delay equals 1.
The System Identification Toolbox product estimates the state-space matrices A, B, C, D, and K using the model order and the data you specify.
The state-space model structure is a good choice for quick estimation because it contains only two parameters: n is the number of poles (the size of the A matrix) and nk is the delay.
Estimating State-Space Models Using n4sid. Use the n4sid command to specify a range of model orders and evaluate the performance of several state-space models (orders 2 to 8):
mn4sid = n4sid(Ze1,2:8,'InputDelay',[4 9]);
This command uses the fast, noniterative (subspace) method and opens the following plot. You use this plot to decide which states provide a significant relative contribution to the input/output behavior, and which states provide the smallest contribution.
The vertical axis is a relative measure of how much each state contributes to the input/output behavior of the model (log of singular values of the covariance matrix). The horizontal axis corresponds to the model order n. This plot recommends n=4, indicated by a red rectangle.
To select this model order, select 4 from the Model Order drop-down list and click Apply.
By default, n4sid uses a free parameterization of the state-space form. To estimate a canonical form instead, set the value of the SSParameterization property to 'Canonical' . You can also specify the input-to-output delay (in samples) using the 'InputDelay' property.
mCanonical = n4sid(Ze1, 3,'SSParameterization', 'canonical','InputDelay', [4 9]); present(mCanonical); % Display model properties
mCanonical = Discrete-time identified state-space model: x(t+Ts) = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x3 x1 0 1 0 x2 0 0 1 x3 0.0738 +/- 0.05922 -0.6083 +/- 0.1624 1.445 +/- 0.1284 B = ConsumptionR Current x1 1.844 +/- 0.1756 0.5631 +/- 0.1225 x2 1.064 +/- 0.1679 2.309 +/- 0.1226 x3 0.2769 +/- 0.0966 1.878 +/- 0.1061 C = x1 x2 x3 Production 1 0 0 D = ConsumptionR Current Production 0 0 K = Production x1 0.8676 +/- 0.03183 x2 0.6864 +/- 0.04166 x3 0.5113 +/- 0.04379 Input delays (sampling periods): 4 9 Sample time: 0.5 minutes Parameterization: CANONICAL form with indices: 3. Feedthrough: none Disturbance component: estimate Number of free coefficients: 12 Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using N4SID on time domain data "Ze1". Fit to estimation data: 91.37% (prediction focus) FPE: 2.456, MSE: 2.339 More information in model's "Report" property.
Note: mn4sid and mCanonical are discrete-time models. To estimate a continuous-time model, set the 'Ts' property to 0 in the estimation command, or use the ssest command: mCT1 = n4sid(Ze1, 3, 'Ts', 0, 'InputDelay', [2.5 5]) mCT2 = ssest(Ze1, 3,'InputDelay', [2.5 5]) |
Learn More. To learn more about estimating state-space models, see State-Space Models.
About Box-Jenkins Models. The general Box-Jenkins (BJ) structure is:
$$y(t)={\displaystyle \sum _{i=1}^{nu}\frac{{B}_{i}(q)}{{F}_{i}(q)}{u}_{i}\left(t-n{k}_{i}\right)}+\frac{C(q)}{D(q)}e(t)$$
To estimate a BJ model, you need to specify the parameters n_{b}, n_{f}, n_{c}, n_{d}, and n_{k}.
Whereas the ARX model structure does not distinguish between the poles for individual input/output paths, the BJ model provides more flexibility in modeling the poles and zeros of the disturbance separately from the poles and zeros of the system dynamics.
Estimating a BJ Model Using pem. You can use pem to estimate the BJ model. pem is an iterative method and has the following general syntax:
pem(data,'na',na,'nb',nb,'nc',nc,'nd',nd,'nf',nf,'nk',nk)
To estimate the BJ model, type:
na = 0; nb = [ 2 1 ]; nc = 1; nd = 1; nf = [ 2 1 ]; nk = [ 5 10]; mbj = polyest(Ze1,[na nb nc nd nf nk]);
This command specifies nf=2 , nb=2 , nk=5 for the first input, and nf=nb=1 and nk=10 for the second input.
Display the model information.
present(mbj)
mbj = Discrete-time BJ model: y(t) = [B(z)/F(z)]u(t) + [C(z)/D(z)]e(t) B1(z) = 1.823 (+/- 0.1792) z^-5 - 1.315 (+/- 0.2368) z^-6 B2(z) = 1.791 (+/- 0.06435) z^-10 C(z) = 1 + 0.1066 (+/- 0.04015) z^-1 D(z) = 1 - 0.7453 (+/- 0.027) z^-1 F1(z) = 1 - 1.321 (+/- 0.06939) z^-1 + 0.5911 (+/- 0.05516) z^-2 F2(z) = 1 - 0.8314 (+/- 0.006445) z^-1 Sample time: 0.5 minutes Parameterization: Polynomial orders: nb=[2 1] nc=1 nd=1 nf=[2 1] nk=[5 10] Number of free coefficients: 8 Use "polydata", "getpvec", "getcov" for parameters and their uncertainties. Status: Termination condition: Near (local) minimum, (norm(g) < tol). Number of iterations: 6, Number of function evaluations: 13 Estimated using POLYEST on time domain data "Ze1". Fit to estimation data: 90.75% (prediction focus) FPE: 2.757, MSE: 2.689 More information in model's "Report" property.
The uncertainty for each of the model parameters is computed to 1 standard deviation and appears in parentheses next to each parameter value.
The polynomials C and D give the numerator and the denominator of the noise model, respectively.
Tip Alternatively, you can use the following shorthand syntax that specifies the orders as a single vector: mbj = bj(Ze1,[2 1 1 1 2 1 5 10]); bj is a version of pem that specifically estimates the BJ model structure. |
Learn More. To learn more about identifying input-output polynomial models, such as BJ, see Input-Output Polynomial Models.
Compare the output of the ARX, state-space, and Box-Jenkins models to the measured output.
compare(Zv1,marx,mn4sid,mbj)
compare plots the measured output in the validation data set against the simulated output from the models. The input data from the validation data set serves as input to the models.
To perform residual analysis on the ARX model, type the following command:
resid(Zv1,marx)
Because the sample system has two inputs, there are two plots showing the cross-correlation of the residuals with each input. Press Enter to view the cross-correlation with the second input.
To perform residual analysis on the state-space model, type the following command:
resid(Zv1,mn4sid)
Finally, to perform residual analysis on the BJ model, type the following command:
resid(Zv1,mbj)
All three models simulate the output equally well and have uncorrelated residuals. Therefore, choose the ARX model because it is the simplest of the three input-output polynomial models and adequately captures the process dynamics.
In this portion of the tutorial, you simulate the model output. You must have already created the continuous-time model midproc2, as described in Estimating Process Models.
Simulating the model output requires the following information:
Input values to the model
Initial conditions for the simulation (also called initial states)
For example, the following commands use the iddata and idinput commands to construct an input data set, and use sim to simulate the model output:
% Create input for simulation U = iddata([],idinput([200 2]),'Ts',0.5,'TimeUnit','min'); % Simulate the response setting initial conditions equal to zero ysim_1 = sim(midproc2,U);
To maximize the fit between the simulated response of a model to the measured output for the same input, you can compute the initial conditions corresponding to the measured data. The best fitting initial conditions can be obtained by using findstates(idParametric) on the state-space version of the estimated model. The following commands estimate the initial states X0est from the data set Zv1:
% State-space version of the model needed for computing initial states
midproc2_ss = idss(midproc2);
X0est = findstates(midproc2_ss, Zv1);
Next, simulate the model using the initial states estimated from the data.
% Simulation input Usim = Zv1(:,[],:); Opt = simOptions('InitialCondition',X0est); ysim_2 = sim(midproc2_ss, Usim, Opt);
Compare the simulated and the measured output on a plot.
figure plot([ysim_2.y, Zv1.y]) legend({'model output','measured'}) xlabel('time'), ylabel('Output')
Many control-design applications require you to predict the future outputs of a dynamic system using the past input/output data.
For example, use predict to predict the model response five steps ahead:
predict(midproc2,Ze1,5)
Compare the predicted output values with the measured output values. The third argument of compare specifies a five-step-ahead prediction. When you do not specify a third argument, compare assumes an infinite prediction horizon and simulates the model output instead.
compare(Ze1,midproc2,5)
Use pe to compute the prediction error Err between the predicted output of midproc2 and the measured output. Then, plot the error spectrum using the spectrum command.
[Err] = pe(midproc2,Zv1); spectrum(spa(Err,[],logspace(-2,2,200)))
The prediction errors are plotted with a 1-standard-deviation confidence interval. The errors are greater at high frequencies because of the high-frequency nature of the disturbance.