Estimate a process model and compare its response with the measured output.

Load the input/output data, which is stored in an iddata object z1.

load iddata1 z1

Estimate a first-order process model sys that contains one pole and no zeroes or delays. This model structure has type P1.

sys = procest(z1,'P1');

Compare the simulated model response with the measured output.

compare(z1,sys)

The fit percentage for the model is low. Add a delay to the model and compare the simulated and measured outputs.

sys = procest(z1,'P1D');
compare(z1,sys)

The fit percentage has improved, but is still below 50%. The plot shows that the model output peaks do not attain the height of the measured output peaks, which indicates that the model needs to include more dynamics.

Create a second-order process model with complex (underdamped) poles.

sys = procest(z1,'P2U');
compare(z1,sys)

The fit now exceeds 70%.

You can view more information about the estimation by exploring the idproc property sys.Report.

sys.Report

ans = 
              Status: 'Estimated using PROCEST'
              Method: 'PROCEST'
    InitialCondition: 'zero'
                 Fit: [1x1 struct]
          Parameters: [1x1 struct]
         OptionsUsed: [1x1 idoptions.procest]
           RandState: []
            DataUsed: [1x1 struct]
         Termination: [1x1 struct]

View the estimated gain Kp.

Kp = sys.Kp

Kp = 7.6818

Estimate a process model after specifying initial guesses for parameter values and bounding them.

Obtain input/output data.

data = idfrd(idtf([10 2],[1 1.3 1.2],'iod',0.45),logspace(-2,2,256));

Specify the parameters of the estimation initialization model.

type = 'P2UZD';
init_sys = idproc(type);

init_sys.Structure.Kp.Value = 1; 
init_sys.Structure.Tw.Value = 2; 
init_sys.Structure.Zeta.Value = 0.1; 
init_sys.Structure.Td.Value = 0; 
init_sys.Structure.Tz.Value = 1; 
init_sys.Structure.Kp.Minimum = 0.1;
init_sys.Structure.Kp.Maximum = 10;
init_sys.Structure.Td.Maximum = 1;
init_sys.Structure.Tz.Maximum = 10;

Specify the estimation options.

opt = procestOptions('Display','full','InitialCondition','Zero');
opt.SearchMethod = 'lm'; 
opt.SearchOptions.MaxIterations = 100;

Estimate the process model.

sys = procest(data,init_sys,opt);

Since the 'Display' option is specified as 'full', the estimation progress is displayed in a separate Plant Identification Progress window.

Compare the data to the estimated model.

compare(data,sys);

Load the data.

load iddata1ic z1i

Estimate a first-order plus dead time process model sys and return the initial conditions in ic. First specify 'estimate' for 'InitialCondition' to force the software to estimate ic. The default 'auto' setting uses the 'estimate' method only when the influence of the initial conditions on the overall model error exceed a threshold. When the initial conditions have a negligible effect on the overall estimation-error minimization process, the 'auto' setting uses 'zero'.

opt = procestOptions('InitialCondition','estimate');
[sys,offset,ic] = procest(z1i,'P1D',opt);
ic

ic = 
  initialCondition with properties:

     A: -3.8997
    X0: -1.0871
     C: 4.5652
    Ts: 0

ic is an initialCondition object that encapsulates the free response of sys, in state-space form, to the initial state vector in X0. You can incorporate ic when you simulate sys with the z1i input signal and compare the response with the z1i output signal.

Obtain input/output data.

load iddata1 z1
load iddata2 z2
data = [z1 z2(1:300)];

data is a data set with 2 inputs and 2 outputs. The first input affects only the first output. Similarly, the second input affects only the second output.

In the estimated process model, the cross terms, which model the effect of the first input on the second output and vice versa, should be negligible. If the estimation process instead assigns higher orders to the cross dynamics, the degrees of estimation uncertainty for those terms should be high.

Estimate the process model.

type = 'P2UZ'; 
sys = procest(data,type);

The type variable denotes a model with complex-conjugate pair of poles, a zero, and a delay.

To evaluate the uncertainties, plot the frequency response.

w = linspace(0,20*pi,100);
h = bodeplot(sys,w);
showConfidence(h);

The responses from the cross pairs show larger uncertainty, indicating that using a single type for each input/output pair results in too much energy in the cross pairs.

Use regularization to estimate parameters of an overparameterized process model.

Load the data.

load iddata1 z1;

Construct an initial system sysi by specifying parameter values for a model that includes three poles, one zero, and underdamped modes. Assume that gain Kp is known with a higher degree of confidence than the other model parameters.

sysi = idproc('P3UZ','Kp',7.5,'Tw',0.25,'Zeta',0.3,'Tp3',20,'Tz',0.02);

Estimate an unregularized process model sys1 using sysi to initialize the estimation model.

sys1 = procest(z1,sysi);

Estimate a regularized process model sys2 from sysi. Because K has a higher level of confidence, set the regularization constant R higher than for the other model parameters. This setting causes the estimation process to place more emphasis on maintaining the initial value of K.

opt = procestOptions;
opt.Regularization.Nominal = 'model';
opt.Regularization.R = [100;1;1;1;1];
opt.Regularization.Lambda = 0.1;
sys2 = procest(z1,sysi,opt);

Compare the model outputs with data.

compare(z1,sys1,sys2);

Regularization helps steer the estimation process towards the correct parameter values, as the better fit for sys2 shows.

Compare the estimated gain values for sys1 and sys2.

g1 = sys1.Kp

g1 = -0.2320

g2 = sys2.Kp

g2 = 6.6236

The Kp value for the regularized system is much closer to the initial value than for the unregularized system.

Obtain the measured input-output data.

load iddemo_heatexchanger_data;
data = iddata(pt,ct,Ts);  
data.InputName  = '\Delta CTemp';
data.InputUnit  = 'C';
data.OutputName = '\Delta PTemp';
data.OutputUnit = 'C';
data.TimeUnit   = 'minutes';

Estimate a first-order plus dead time process model.

type = 'P1D';
sysP1D = procest(data,type);

Compare the model with the data.

compare(data,sysP1D)

Plot the model residuals.

figure
resid(data,sysP1D);

The figure shows that the residuals are correlated. To account for that, add a first order ARMA disturbance component to the process model.

opt = procestOptions('DisturbanceModel','ARMA1');
sysP1D_noise = procest(data,'p1d',opt);

Compare the models.

compare(data,sysP1D,sysP1D_noise)

Plot the model residuals.

figure
resid(data,sysP1D_noise);

The residues of sysP1D_noise are uncorrelated.