Estimating Coefficients of ODEs to Fit Given Solution

This example shows how to estimate model parameters using linear and nonlinear grey-box modeling.

Use grey-box identification to estimate coefficients of ODEs that describe the model dynamics to fit a given response trajectory.

  • For linear dynamics, represent the model using a linear grey-box model (idgrey). Estimate the model coefficients using greyest.

  • For nonlinear dynamics, represent the model using a nonlinear grey-box model (idnlgrey). Estimate the model coefficients using pem.

In this example, you estimate the value of the friction coefficient of a simple pendulum using its oscillation data. The motion of a simple pendulum is described by:

ml2θ¨+bθ˙+mglsinθ=0

θ is the angular displacement of the pendulum relative to its state of rest. g is the gravitational acceleration constant. m is the mass of the pendulum and l is the length of the pendulum. b is the viscous friction coefficient whose value will be estimated to fit the given angular displacement data. There is no external driving force that is contributing to the pendulum's motion.

Load measured data.

load(fullfile(matlabroot, 'toolbox', 'ident', ...
     'iddemos', 'data', 'pendulumdata'));
data = iddata(y, [], 0.1, 'Name', 'Pendulum');
data.OutputName = 'Pendulum position';
data.OutputUnit = 'rad';
data.Tstart = 0;
data.TimeUnit = 's';

The measured angular displacement data is loaded and saved as data, an iddata object with a sampling time of 0.1 seconds. The set command is used to specify data attributes such as the output name, output unit and the start time and units of the time vector.

Perform linear grey-box estimation.

Assuming that the pendulum undergoes only small angular displacements, the equation describing the pendulum motion can be simplified:

ml2θ¨+bθ˙+mglθ=0

Using the angular displacement (θ) and the angular velocity (θ˙) as state variables, the simplified equation can be rewritten in the form:

X(t).=AX(t)+Bu(t)y(t)=CX(t)+Du(t)

Here,

X(t)=[θθ˙]A=[01glbml2]B=0C=[10]D=0

The B and D matrices are zero because there is no external driving force for the simple pendulum.

  1. Create an ODE file that relates the model coefficients to its state space representation.

    function [A,B,C,D] = LinearPendulum(m,g,l,b,Ts)
    % Function mapping ODE coefficients to state-space matrices.
    % Save this function as a file on your computer.
    A = [0 1; -g/l, -b/m/l^2];
    B = zeros(2,0);
    C = [1 0];
    D = zeros(1,0);
    end

    The function, LinearPendulum, returns the state space representation of the linear motion model of the simple pendulum using the model coefficients m, g, l and b. The sample time used is specified by Ts.

    Save this function as LinearPendulum.m.

  2. Create a linear grey-box model associated with the LinearPendulum function.

    m = 1;
    g = 9.81;
    l = 1;
    b = 0.2;
    linear_model = idgrey('LinearPendulum', {m,g,l,b},'c');

    m, g and l specify the values of the known model coefficients. b specifies the initial guess for the viscous friction coefficient.

    The function LinearPendulum must be on the MATLAB® path. Alternatively, you can specify the full path name for this function.

    The 'c' input argument in the call to idgrey specifies linear_model as a continuous-time system.

  3. Specify m, g and l as known parameters.

    linear_model.Structure.Parameters(1).Free = false;
    linear_model.Structure.Parameters(2).Free = false;
    linear_model.Structure.Parameters(3).Free = false;

    As defined in the previous step, m, g, and l are the first three parameters of linear_model. Using the Structure.Parameters.Free field for each of the parameters, m, g, and l are specified as fixed values.

  4. Create an estimation options set that specifies the initial state to be estimated and turns on the estimation progress display.

    opt = greyestOptions('InitialState','estimate','Display','on');
    opt.Focus = 'stability';
  5. Estimate the viscous friction coefficient.

    linear_model = greyest(data,linear_model,opt);

    The greyest command updates the parameter of linear_model.

    b_est = linear_model.Structure.Parameters(4).Value;
    [linear_b_est, dlinear_b_est] = getpvec(linear_model, 'free')

    getpvec returns, as dlinear_b_est, the 1 standard deviation uncertainty associated with b, the free estimation parameter of linear_model.

    The estimated value of b, the viscous friction coefficient, using linear grey-box estimation is 0.1178 with a standard deviation of 0.0088.

  6. Compare the response of the linear grey-box model to the measured data.

    compare(data,linear_model)

    The linear grey-box estimation model provides a 49.9% fit to measured data. This is not surprising given that the underlying assumption of the linear pendulum motion model is that the pendulum undergoes small angular displacements, whereas the measured data shows large oscillations.

Perform nonlinear grey-box estimation.

Nonlinear grey-box estimation requires that you express the differential equation as a set of first order equations.

Using the angular displacement (θ) and the angular velocity (θ˙) as state variables, the equation of motion can be rewritten as a set of first order nonlinear differential equations:

x1(t)=θ(t)x2(t)=θ˙(t)x˙1(t)=x2(t)x˙2(t)=glsin(x1(t))bml2x2(t)y(t)=x1(t)

  1. Create an ODE file that relates the model coefficients to its nonlinear representation.

    function [dx,y] = NonlinearPendulum(t,x,u,m,g,l,b,varargin)
    % Function that maps the ODE coefficients to state 
    % variable derivatives and output.
    % Save this function as a file on your computer.
    
    % Output equation.
    y = x(1); % Angular position.
    
    % State equations.
    dx = [x(2);                             ... % Angular position
          -(g/l)*sin(x(1))-b/(m*l^2)*x(2)   ... % Angular velocity
         ];
    end

    The function, NonlinearPendulum, returns the state derivatives and output of the nonlinear motion model of the pendulum using the model coefficients m, g, l and b.

    Save this function as NonlinearPendulum.m.

  2. Create a nonlinear grey-box model associated with the NonlinearPendulum function.

    m = 1; 
    g = 9.81; 
    l = 1; 
    b = 0.2;
    order         = [1 0 2];             
    parameters    = {m,g,l,b};           
    initial_states = [1; 0];             
    Ts            = 0;                   
    nonlinear_model = idnlgrey('NonlinearPendulum', order, ...
                      parameters, initial_states, Ts);
  3. Specify m, g and l as known parameters.

    setpar(nonlinear_model, 'Fixed', {true true true false});

    As defined in the previous step, m, g, and l are the first three parameters of nonlinear_model. Using the setpar command, m, g, and l are specified as fixed values and b is specified as a free estimation parameter.

  4. Estimate the viscous friction coefficient.

    nonlinear_model = pem(data,nonlinear_model,'Display','Full');

    The pem command updates the parameter of nonlinear_model.

    b_est = nonlinear_model.Parameters(4).Value;
    [nonlinear_b_est, dnonlinear_b_est] = ...
                                getpvec(nonlinear_model, 'free')

    getpvec returns, as dnonlinear_b_est, the 1 standard deviation uncertainty associated with b, the free estimation parameter of nonlinear_model.

    The estimated value of b, the viscous friction coefficient, using nonlinear grey-box estimation is 0.1 with a standard deviation of 0.0149.

  5. Compare the response of the linear and nonlinear grey-box models to the measured data.

    compare(data,linear_model,nonlinear_model)

    The nonlinear grey-box model estimation provides a closer fit (95%) to the measured data.

Was this topic helpful?