Non linear greybox modelling: iddata - input arguments

Warning and error: SIP Warning: There are more output channels in the data object than samples. Check if the output data matrix should be transposed. > In warning at 26 In @iddata\private\datachk at 35 In iddata.iddata at 192 In SIP at 10 Warning: There are more output channels in the data object than samples. Check if the output data matrix should be transposed. > In warning at 26 In @iddata\private\datachk at 35 In iddata.pvset at 138 In iddata.set at 143 In iddata.iddata at 227 In SIP at 10 ipfd Error using iddata (line 229) The number of elements in the "OutputUnit" property value must be equal to the number of output channels. Use a cell array of 1000000 string(s).

Error in SIP (line 10) sway = iddata(O, [], 0.001, 'Name', 'Sway', 'OutputUnit', 'rad');

Error using ipfd (line 34) Not enough input arguments.

    %%An inverted pendulum model directly applied to normal human gait
    %Model predicts values of ground reaction forces and velocity of human gait.
    %An inverted pendulum is defined basing on the average center-of-pressure 
    %to the instantaneous center-of-mass and equations of motion during single 
    %support that allowed a telescoping action were derived.

    %%Data

    load('O');
    sway = iddata(O, [], 0.001, 'Name', 'Sway', 'OutputUnit', 'rad');
    set(sway, 'Tstart', 0, 'TimeUnit', 's');

    load('r');
    leg = iddata(r, [], 0.001, 'Name', 'Leg length', 'OutputUnit', 'm');
    set(leg, 'Tstart', 0, 'TimeUnit', 's');

    stab = load('stab');
    %stab = [g, w0, dt, m]
    g = stab(1);
    set(g, 'OutputName', 'Gravity constant');
    set(g, 'OutputUnit', 'm/s^2');
    w0 = stab(2);
    set(w0, 'OutputName', 'Initial angular velocity');
    set(w0, 'OutputUnit', 'rad/s');
    dt = stab(3);
    set(dt, 'OutputName', 'Integration step');
    set(dt, 'OutputUnit', 's');
    m = stab(4);
    set(m, 'OutputName', 'Weight');
    set(m, 'OutputUnit', 'kg');
    %%Pendulum Modeling
    % This paragraph presents a model structure for the pendulum system.
    % The forward and inverse dynamics of it was studied by the same authors 
    %before. The structure is as follows:
    %
    %a(i) = (-g*cos(O(i))/r
    %w(i+1) = w(i) + a(i)*dt
    %O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
    %Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
    %Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
    %v = w*r*sin(O)
    %
    % having parameters (or constants):
    %    a - angular acceleration [rad/s^2]
    %    g - the gravity constant [m/s^2]
    %    O - sway [rad]
    %    r - the length of the leg [m]
    %    w - angular velocity [rad/s]
    %    dt - integration steps [s]
    %    m - the weight [kg]
    %    Fh - horizontal ground reaction force [
    %    Fv - vertical ground reaction force
    %    v - velocity

    %We enter this information into the MATLAB file ipfd.m.

The next step is to create a generic IDNLGREY object for describing the pendulum. We also enter information about the names and the units of the inputs, states, outputs and parameters. Owing to the physical reality all model parameters must be positive and this is imposed by setting the 'Minimum' property of all parameters to the smallest recognizable positive value in MATLAB®, eps(0).

    FileName      = 'ipfd';        % File describing the model structure.
    Order         = [1 0 2];             % Model orders [ny nu nx]. ???
    Parameters    = [g; w0; d; m];   % Initial parameters.
    InitialStates = [1; 0];              % Initial states.
    Ts            = 0;                   % Time-continuous system.
    nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts);
    %Output [Fx, Fy, v, x, y]
    set(nlgr, 'OutputName', {'Horizontal ground reaction force' 'Vertical ground reaction force' 'velocity' 'Horizontal position' 'Vertical position'}, ...
    'OutputUnit', {'N' 'N' 'm/s' 'm' 'm'},            ...
    'TimeUnit', 's'); %output
    nlgr = setinit(nlgr, 'Name', {'Sway' 'Leg length'});
    nlgr = setinit(nlgr, 'Unit', {'rad' 'm'});
    nlgr = setpar(nlgr, 'Name', {'Gravity constant' 'Integration step' ...
    'Initial angular velocity' 'Mass'}); %parameters
    nlgr = setpar(nlgr, 'Unit', {'m/s^2' 's' 'rad/s' 'kg'}); %units
    nlgr = setpar(nlgr, 'Minimum', {eps(0) eps(0) eps(0) eps(0)});   % All parameters > 0.

    %%Performance of Pendulum Model
    % Before trying to estimate any parameter we simulate the system with the
    % guessed parameter values. The solver used is Runge-Kutta 45 with adaptive 
    %step length(ode45). The Runge-Kutta 45 (ode45) solver often returns high 
    %quality results relatively fast, and is therefore chosen as the default 
    %differential equation solver in IDNLGREY.

    % C. Model computed with adaptive Runge-Kutta 45 ODE solver.
    nlgrrk45 = nlgr;
    nlgrrk45.Algorithm.SimulationOptions.Solver = 'ode45';     % Runge-Kutta 45.

    %%Parameter Estimation
    %Some parameters fixed, other not for simulation
    nlgrrk45 = setpar(nlgrrk45, 'Fixed', {true true false true});

    trk45 = clock;
    nlgrrk45 = pem(z, nlgrrk45, 'Display', 'Full');   % Perform parameter estimation.
    trk45 = etime(clock, trk45);
    disp('           Estimation time   Estimated value');
    fprintf('   ode45:  %3.1f               %1.3f\n', trk45, nlgrrk45.Parameters(3).Value);
    %the Final Prediction Error (FPE) criterion
    fpe(nlgrrk45);

    figure
    t = linspace(0, 1000, 1000000);
    subplot(3,1,1);
    plot(t,Fh)

    subplot(3,1,2);
    plot(t,Fv)

    subplot(3,1,3);
    plot(t,v)

[Walter Roberson]

This code makes use of material from the paper,

Performance of an inverted pendulum model directly applied to normal human gait. Buczek FL, Cooney KM, Walker MR, Rainbow MJ, Concha MC, Sanders JO. Performance of an inverted pendulum model directly applied to normal human gait. Clin Biomech (Bristol, Avon). 2006 Mar; 21(3):288-96.

There does not appear to a free version of the paper available.

Answer 1

[Walter Roberson]

Each of the places you have a single OutputUnit, such as

sway = iddata(O, [], 0.001, 'Name', 'Sway', 'OutputUnit', 'rad');

change that to be a cell array containing the string:

sway = iddata(O, [], 0.001, 'Name', 'Sway', 'OutputUnit', {'rad'});

It appears that the content of O is a row vector. You need to pass in a column vector instead of a row vector.

sway = iddata(O.', [], 0.001, 'Name', 'Sway', 'OutputUnit', {'rad'});

You might encounter the same difficulty in the other datasets. If you know that they are each vectors but you do not know the orientation, then you can force column vector by

sway = iddata(O(:), [], 0.001, 'Name', 'Sway', 'OutputUnit', {'rad'})

[Katarzyna Wieciorek]

Could you please tell me how to correct OutputName and Name, is it question of comas or semicolons?

Error: Error using idnlmodel/set (line 89) The value of the "OutputName" property must be a string vector with as many entries as outputs. Type "help DynamicSystem.OutputName" for more information.

Error in SIP (line 71) set(nlgr, 'OutputName', {'Horizontal ground reaction force' 'Vertical ground reaction force' 'velocity' 'Horizontal position' 'Vertical position'}, ...

%Output [Fx, Fy, v, x, y]
set(nlgr, 'OutputName', {'Horizontal ground reaction force' 'Vertical ground reaction force' 'velocity' 'Horizontal position' 'Vertical position'}, ...
          'OutputUnit', {'N' 'N' 'm/s' 'm' 'm'},            ...
          'TimeUnit', 's'); %output
nlgr = setinit(nlgr, 'Name', {'Sway' 'Leg length'});
nlgr = setinit(nlgr, 'Unit', {'rad' 'm'});
nlgr = setpar(nlgr, 'Name', {'Gravity constant' 'Integration step' ...
                      'Initial angular velocity' 'Mass'}); %parameters
nlgr = setpar(nlgr, 'Unit', {'m/s^2' 's' 'rad/s' 'kg'}); %units
nlgr = setpar(nlgr, 'Minimum', {eps(0) eps(0) eps(0) eps(0)});   % All parameters > 0.

[Walter Roberson]

When you use Order of [1 0 2], that defines there as being only 1 output. You need exactly one OutputName entry for each output, but you are trying to create 5 OutputNames instead of the 1 mandated by Order of 1. If you want 5 outputs you need to have Order start with 5. Note that your ipfd.m file must be a function returning two outputs, the second of which must be a column vector with the same number of components as you have defined outputs.

[Katarzyna Wieciorek]

I attached ipfd.m to let you see assumptions of the model. I have 5 outputs and 6 inputs, 2 of them vary in time.

"Number of outputs, inputs, and states of the model Order = [Ny Nu Nx]. For time series, Nu is set to 0, and for static model structures, Nx is set to 0."

My model have 2 time series and 4 constant variables as input, 5 variables as output, what is the order then?

Thank you in advance.

[Walter Roberson]

You can't return 5 MATLAB variables from the function you name as a idnlgrey model. Look at the documentation on how it has to be structured: http://www.mathworks.com/help/ident/ug/estimating-nonlinear-grey-box-models.html#bq5_fno

The output variables are:

    dx — Represents the right side(s) of the state-space equation(s). A column vector with Nx entries. For static models, dx=[].

    For discrete-time models. dx is the value of the states at the next time step x(t+Ts).

    For continuous-time models. dx is the state derivatives at time t, or dxdt.

    y — Represents the right side(s) of the output equation(s). A column vector with Ny entries.

You have 2 states, so your Nx is 2. You have 5 model outputs so your Ny is 5. You have a time series so your Nu is set to 0. Your Order should therefore be [5 0 2]. And your function needs to be altered as follows (omitting the comments). Note the 4 lines I marked CHANGED

function [dx, Y] = ipfd(g, O, r, w0, dt, m)  %CHANGED
%State equations
a = (-g*cos(O))/r; %angular acceleration, not output
w = w0 + a*dt; %Angular velocity
O = O + w*dt + 1/2*a*(dt.^2);  %CHANGED: not O(0), cannot index at 0

%Output equations
Fx = m.*((diff(r, 2)-r.*w.^2)*cos(O)-(r.*a+2.*diff(r, 1)*w)*sin(O));
Fy = m.*((diff(r, 2)-r.*w.^2)*sin(O)-(r.*a+2.*diff(r, 1)*w)*cos(O)) + m.*g;
v = w*r*sin(O); 
x = r*cos(O);
y = r*sin(O);

dx = [w;O];               %CHANGED
Y = [Fx; Fy; v; x; y];    %CHANGED

[Katarzyna Wieciorek]

Thank you very much.

I have one problem more, what variable should I use in this line to not get error while solving ivp? for Y it doesn't work, for w and w0 too.

nlgrrk45 = pem(?, nlgrrk45, 'Display', 'Full');

more:

%%Parameter Estimation
%Some parameters fixed, other not for simulation
nlgrrk45 = setpar(nlgrrk45, 'Fixed', {true false true true});

trk45 = clock;
nlgrrk45 = pem(Y, nlgrrk45, 'Display', 'Full');   % Perform parameter estimation.
trk45 = etime(clock, trk45);
disp('           Estimation time   Estimated value');
fprintf('   ode45:  %3.1f               %1.3f\n', trk45, nlgrrk45.Parameters(3).Value);
%the Final Prediction Error (FPE) criterion
fpe(nlgrrk45);

Thank you in advance.

[Walter Roberson]

What you pass in to pem() in the first parameter needs to be something created by iddata(). In your case, that would be either sway or leg.

I suspect this might still fail in your situation, but I need to keep reading so I will post this now. If my suspicion is right then I need to know that size(O) is the same as size(r) and that both are vectors.

[Walter Roberson]

I think you are going to need to construct your iddata differently. I think you are going to need:

    load('O');
    load('r');

    swayleg = iddata([O(:),r(:)], [], 0.001, 'Name', 'Sway/Leg', 'OutputUnit', {'m', 'rad'});
    set(swayleg, 'Tstart', 0, 'TimeUnit', 's');

And then I think you are going to use

 nlgrrk45 = pem(swayleg, nlgrrk45, 'Display', 'Full');

This will fail if your O and r are not the same length as each other. If they are not vectors then that would be important to know for interpretation purposes.

[Paulo Costa]

Now the error is this:

      >> main
      Error using idnlgrey/pem (line 127)
      The number of inputs and outputs of the model must match that of the data.

      Error in pem (line 54)
         [sys, varargout{1:nargout-1}] = pem(varargin{:});

      Error in main (line 94)
      nlgrrk45 = pem(swayleg, nlgrrk45, 'Display', 'Full'); % nlgrrk45 = pem(Y, nlgrrk45,
      'Display', 'Full');   % Perform parameter estimation.

And yes, O and r are vectors.

What could we do?

[Walter Roberson]

Sorry this is all new to me.

What I can see is that to do parameter estimation, you need a set of measured data outputs, one column per model output, and one row per time step. You have 5 outputs from your model so you would need 5 columns. But what you have instead is 2 columns, and those 2 columns appear to correspond to state rather than to outputs, but I am not sure of that yet. I think you might be needing to use those as InputData rather than as OutputData, I am still checking on that.

I also see now that there are errors in your ifpd. I think I am close to being able to fix those.

[Walter Roberson]

Revised. There was a parameter "u" missing, and the state was not being decoded properly.

function [dx, Y] = ipfd(g, wO, u, r, w0, dt, m)
%State equations
%    O - sway [rad]
%    w - angular velocity [rad/s]
w = wO(1);
O = wO(2);
a = (-g*cos(O))/r; %angular acceleration, not output
w = w + a*dt; %Angular velocity
O = O + w*dt + 1/2*a*(dt.^2);

Everything above here changed. The rest stayed the same

%Output equations
Fx = m.*((diff(r, 2)-r.*w.^2)*cos(O)-(r.*a+2.*diff(r, 1)*w)*sin(O));
Fy = m.*((diff(r, 2)-r.*w.^2)*sin(O)-(r.*a+2.*diff(r, 1)*w)*cos(O)) + m.*g;
v = w*r*sin(O); 
x = r*cos(O);
y = r*sin(O);

dx = [w;O];
Y = [Fx; Fy; v; x; y];

[Walter Roberson]

Dang, the above is going to have to be fixed again. I'm working on it.

[Walter Roberson]

When you do the

stab = load('stab');

what does

class(stab)

indicate?

Also, please post your current version of your code. The more I read the more it looks to me like you could not have gotten as far as you did with the code that you posted.

[Paulo Costa]

Stab does nothing, we erased it.

This is the new code of "ipfd":

function [dx, Y] = ipfd(g, O, r, w0, dt, m)
%A pendulum system

%Model:
%a(i) = (-g*cos(O(i))/r
%w(i+1) = w(i) + a(i)*dt
%O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
%Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
%Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
%v = w*r*sin(O)

%State equations
a = (-g*cos(O))/r; %angular acceleration, not output
w = w0 + a*dt; %Angular velocity
O = O + w*dt + 1/2*a*(dt.^2);  %cannot index at 0
%Output equations
Fx = m.*((diff(r, 2)-r.*w.^2)*cos(O)-(r.*a+2.*diff(r, 1)*w)*sin(O));
Fy = m.*((diff(r, 2)-r.*w.^2)*sin(O)-(r.*a+2.*diff(r, 1)*w)*cos(O)) + m.*g;
v = w*r*sin(O); 
x = r*cos(O);
y = r*sin(O);
dx = [O; r];               %input
Y = [Fx; Fy; v; x; y];    %output
end

[Paulo Costa]

And this is the newest code for "main":

   %%An inverted pendulum model directly applied to normal human gait
  %Model predicts values of ground reaction forces and velocity of human gait.
  %An inverted pendulum is defined basing on the average center-of-pressure 
  %to the instantaneous center-of-mass and equations of motion during single 
  %support that allowed a telescoping action were derived.

%%Data

% load('O');
% O = iddata(O.', [], 0.001, 'Name', 'Sway', 'OutputUnit', {'rad'});
% set(sway, 'Tstart', 0, 'TimeUnit', 's');
% 
% load('r');
% r = iddata(r.', [], 0.001, 'Name', 'Leg length', 'OutputUnit', {'m'});
% set(leg, 'Tstart', 0, 'TimeUnit', 's');

load('O');
load('r');
swayleg = iddata([O(:),r(:)], [], 0.001, 'Name', 'Sway/Leg', 'OutputUnit', {'m', 'rad'});
set(swayleg, 'Tstart', 0, 'TimeUnit', 's');

%%Pendulum Modeling
% This paragraph presents a model structure for the pendulum system.
% The forward and inverse dynamics of it was studied by the same authors 
%before. The structure is as follows:
%
%a(i) = (-g*cos(O(i))/r
%w(i+1) = w(i) + a(i)*dt
%O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
%Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
%Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
%v = w*r*sin(O)
%
% having parameters (or constants):
%    a - angular acceleration [rad/s^2]
%    g - the gravity constant [m/s^2]
%    O - sway [rad]
%    r - the length of the leg [m]
%    w - angular velocity [rad/s]
%    dt - integration steps [s]
%    m - the weight [kg]
%    Fh - horizontal ground reaction force [
%    Fv - vertical ground reaction force
%    v - velocity

%We enter this information into the MATLAB file ipfd.m.

The next step is to create a generic IDNLGREY object for describing the pendulum. We also enter information about the names and the units of the inputs, states, outputs and parameters. Owing to the physical reality all model parameters must be positive and this is imposed by setting the 'Minimum' property of all parameters to the smallest recognizable positive value in MATLAB®, eps(0).

FileName      = 'ipfd';        % File describing the model structure.
Order         = [5 0 2];             % Model orders [ny nu nx]. 
%Vector [Ny Nu Nx], specifying the number of model outputs Ny, inputs Nu, 
%and states Nx.
Parameters    = [9.81 5.86 0.001 39.3];   % Initial parameters.
InitialStates = [1; 0];              % Initial states.
Ts            = 0;                   % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts);
%Output [Fx, Fy, v, x, y]
set(nlgr, 'OutputName', {'Horizontal ground reaction force' 'Vertical ground reaction force' 'velocity' 'Horizontal position' 'Vertical position'}, ...
          'OutputUnit', {'N' 'N' 'm/s' 'm' 'm'},            ...
          'TimeUnit', 's'); %output
nlgr = setinit(nlgr, 'Name', {'Sway' 'Leg length'});
nlgr = setinit(nlgr, 'Unit', {'rad' 'm'});
nlgr = setpar(nlgr, 'Name', {'Gravity constant' 'Integration step' ...
                      'Initial angular velocity' 'Mass'}); %parameters
nlgr = setpar(nlgr, 'Unit', {'m/s^2' 's' 'rad/s' 'kg'}); %units
nlgr = setpar(nlgr, 'Minimum', {eps(0) eps(0) eps(0) eps(0)});   % All parameters > 0.

%%Simulation
%Simulation before estimating parameters
data = [O; r];
% sim(nlgr, data);
%%Performance of Pendulum Model
% Before trying to estimate any parameter we simulate the system with the
% guessed parameter values. The solver used is Runge-Kutta 45 with adaptive 
%step length(ode45). The Runge-Kutta 45 (ode45) solver often returns high 
%quality results relatively fast, and is therefore chosen as the default 
%differential equation solver in IDNLGREY.

% C. Model computed with adaptive Runge-Kutta 45 ODE solver.
nlgrrk45 = nlgr;
nlgrrk45.Algorithm.SimulationOptions.Solver = 'ode45';     % Runge-Kutta 45.

%%Parameter Estimation
%Some parameters fixed, other not for simulation
nlgrrk45 = setpar(nlgrrk45, 'Fixed', {true false true true});

trk45 = clock;
nlgrrk45 = pem(swayleg, nlgrrk45, 'Display', 'Full'); % nlgrrk45 = pem(Y, nlgrrk45, 'Display', 'Full');   % Perform parameter estimation.
trk45 = etime(clock, trk45);
disp('           Estimation time   Estimated value');
fprintf('   ode45:  %3.1f               %1.3f\n', trk45, nlgrrk45.Parameters(3).Value);
%the Final Prediction Error (FPE) criterion
fpe(nlgrrk45);

figure
t = linspace(0, 1000, 1000000);
subplot(3,1,1);
plot(t,Fh)

subplot(3,1,2);
plot(t,Fv)

subplot(3,1,3);
plot(t,v)

[Paulo Costa]

and this is for "fakedata":

%fake data maker

%every time varying variable needs to be saved to separate .mat file
%input variables: g, O, r, w0, dt, m
O = zeros(1,1000000);
O1 = linspace(1.05, 1.92, 500);
for i = 1:1000
    b = (2*i-1)*500;
    a = b-499;
    r(1,a:b) = O1(:);
end
for i = 1:1000
    b = i*1000;
    a = b-499;
    r(1,a:b) = fliplr(O1(:));
end

r = zeros(1,1000000); 
%pendulum lenght [m], height = 1.46 [m], leg is half of height = 0.73
r1 = linspace(0.73, 0.65, 500);
for i = 1:1000
    b = (2*i-1)*500;
    a = b-499;
    r(1,a:b) = r1(:);
end
for i = 1:1000
    b = i*1000;
    a = b-499;
    r(1,a:b) = fliplr(r1(:));
end

%stab = [g, w0, dt, m]
%stab = [9.81, 0.3, 0.001, 32]; 

%Saving
save('O.mat','O');
save('r.mat','r');
%save('stab.mat','stab');

[Walter Roberson]

When you do the pem() parameter estimation, are you intending that the initial angular velocity be varied (Fixed = false) or that the integration step be varied (Fixed = false) ? The values that you pass in as your initial parameters are in the order [g, w0, dt, m], but you have

nlgr = setpar(nlgr, 'Name', {'Gravity constant' 'Integration step' ...
  'Initial angular velocity' 'Mass'}); %parameters

which means the order should be [g, dt, w0, m] if we are to match that.

Your model needs r, which the comments in the code list as a constant, which means that it should go into the parameters. And your initial angular velocity should not be a parameter: it should be how you initialize your states: where you currently have

InitialStates = [1; 0];

it should be

InitialStates = [w0; 0];

The equations in the comments do not need the initial angular velocity: it only forms the initial condition for the state. But if you move it out of parameters then there is no way it could be varied within pem(), which is why it becomes particularly important whether you had intended w0 to be the same parameter that you were varying, why it is important whether the parameter order is [g, w0, dt, m] or [g, dt, w0, m] when you are setting the second parameter to be varied by pem().

The order of Parameters needs to match the order of the variables inside the ipfd function, so when you add in 'r' as a parameter, that function needs to be modified too. (There are other reasons why that file as it is is broken even after the corrections I posted above.)

The iddata() objects you have are for Sway and Leg Length, neither of which are outputs created by your model.

Sway appears to be the same as Sway, O(i) in the model, which is a state in the model, not a measurement.

Leg Length appears to be the same as r in the model, which for the purposes of the equations is a constant, not a state and not a measurement. Notice that the iddata() object "leg" that was constructed from the r vector is not made part of the idnlgrey object in any way -- it is not passed in to the initial object and it is not set*() to any of the objects.

So what's going on here? It kind of looks like you want Sway and Leg Length to be inputs O and r. Since O is currently a state using Sway as input would be physically correspond to "forcing". If you do want to use Sway as an input then it would need to be moved to the "u" parameter, as in Nu.

I would have thought that if you use Sway as an input that you would have to switch to discrete time rather than continuous time because each Sway value would correspond to an input at a specific time step, and that you would have to adjust the Ts parameter to equal the time step of the Sway measurements, but I see an example on the idnlgrey() page for twotanks_c which uses an input and which is continuous time. Unfortunately I do not have access to the demo source to study how that was done or what it means.

I also see from that demo example that it is valid to have states (Nx = 2 in the example) while having inputs (Nu = 1 in the example), which implies that the model would not be considered either a time series (because Nu must be 0 for a time series) nor a static model (because Nx must be 0 for a static model). Maybe that's okay; I am still new to this.

If you want to use the input data "r" in the model, it too would have to become an input, added to u, Nu > 0. It would seem to physically correspond to people shortening their leg while walking, which does not seem physically realistic unless it is accounting for the leg effectively becoming slightly longer as people go up on to their toes as they move faster. Yes, I could see that having an effect, but I would have though it would be a lot more important to incorporate a "spring" effect, that as the swing reaches further the pull against the elastic muscles dampen the motion. Whatever; failure of the model to take some things into effect is your problem, not mine :=p

[Walter Roberson]

In your fakedata routine you are leaving O set as all 0.

Change

O = zeros(1,1000000);
O1 = linspace(1.05, 1.92, 500);
for i = 1:1000
    b = (2*i-1)*500;
    a = b-499;
    r(1,a:b) = O1(:);
end
for i = 1:1000
    b = i*1000;
    a = b-499;
    r(1,a:b) = fliplr(O1(:));
end

to

O = zeros(1,1000000);
O1 = linspace(1.05, 1.92, 500);
for i = 1:1000
    b = (2*i-1)*500;
    a = b-499;
    O(1,a:b) = O1(:);
end
for i = 1:1000
    b = i*1000;
    a = b-499;
    O(1,a:b) = fliplr(O1(:));
end

[Walter Roberson]

Your ipfd is definitely wrong and the earlier corrections I posted will not be enough to fix it. But before I can post a revised version I will need to know the order of parameters, whether it is w0 followed by dt or if it is dt followed by w0. I will also need to know whether r is to be treated as a parameter (a constant) that you did not add to the list yet or if it is to be treated as an input. I could continue to treat w0 as a parameter that is useless to the model... but if you are wanting to vary it for parameter estimation then things get more complicated.

[Paulo Costa]

From what I understood, the order should be dt followed by w0 and supposedly r should be treated as an input. w0 is the initial angular velocity, it should be a parameter.

[Walter Roberson]

Okay, and what about Sway? Should it be a state or an input?

[Paulo Costa]

From what I've seen sway should be a state.

[Walter Roberson]

Doesn't seem to be much point in generating sway as fake data if it is a state. But anyhow...

[Walter Roberson]

function [dx, Y] = ipfd(t, w_O, r, g, dt, m)   %CHANGED
    %A pendulum system

    %Model:
    %a(i) = (-g*cos(O(i))/r
    %w(i+1) = w(i) + a(i)*dt
    %O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
    %Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
    %Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
    %v = w*r*sin(O)

    %State equations
    w = w_O(1);         %CHANGED
    O = w_O(2);         %CHANGED
    a = (-g*cos(O))/r; %angular acceleration, not output
    w = w0 + a*dt; %Angular velocity
    O = O + w*dt + 1/2*a*(dt.^2);  %cannot index at 0
    %Output equations
    Fx = m.*((diff(r, 2)-r.*w.^2)*cos(O)-(r.*a+2.*diff(r, 1)*w)*sin(O));
    Fy = m.*((diff(r, 2)-r.*w.^2)*sin(O)-(r.*a+2.*diff(r, 1)*w)*cos(O)) + m.*g;
    v = w*r*sin(O); 
    x = r*cos(O);
    y = r*sin(O);
    dx = [w;O];               %CHANGED
    Y = [Fx; Fy; v; x; y];    %output
  end

[Paulo Costa]

It says that

>> ipfd
Error using ipfd (line 15)
Not enough input arguments.

line 15 is

   w = w_O(1);

[Walter Roberson]

I cut out the parameter estimation because I am getting tired and confused and I'm not sure that it is meaningful. It might be possible to take the output of the simulation as input for parameter estimation, but I think it would require adding the leg to the iddata that results from sim(), which should be possible.

I'm working on this without access to the toolbox itself and I've never looked at any of this before you asked your question. My mind is full for the moment so I will release this back to you.

I believe everything up to and including the sim() should work. I am not as confident that the plotting right after that will work. If you run sim by itself without assigning the output to Y then it should plot by itself.

     %%An inverted pendulum model directly applied to normal human gait
    %Model predicts values of ground reaction forces and velocity of human gait.
    %An inverted pendulum is defined basing on the average center-of-pressure 
    %to the instantaneous center-of-mass and equations of motion during single 
    %support that allowed a telescoping action were derived.

%%Data

% load('O');
% O = iddata(O.', [], 0.001, 'Name', 'Sway', 'OutputUnit', {'rad'});
% set(sway, 'Tstart', 0, 'TimeUnit', 's');

load('r');
leg = iddata([], r.', 0.001, 'Name', 'Leg length', 'OutputUnit', {'m'});
set(leg, 'Tstart', 0, 'TimeUnit', 's');

%%Pendulum Modeling
% This paragraph presents a model structure for the pendulum system.
% The forward and inverse dynamics of it was studied by the same authors 
%before. The structure is as follows:
%
%a(i) = (-g*cos(O(i))/r
%w(i+1) = w(i) + a(i)*dt
%O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
%Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
%Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
%v = w*r*sin(O)
%
% having parameters (or constants):
%    a - angular acceleration [rad/s^2]
%    g - the gravity constant [m/s^2]
%    O - sway [rad]
%    r - the length of the leg [m]
%    w - angular velocity [rad/s]
%    dt - integration steps [s]
%    m - the weight [kg]
%    Fh - horizontal ground reaction force [
%    Fv - vertical ground reaction force
%    v - velocity

%We enter this information into the MATLAB file ipfd.m.

The next step is to create a generic IDNLGREY object for describing the pendulum. We also enter information about the names and the units of the inputs, states, outputs and parameters. Owing to the physical reality all model parameters must be positive and this is imposed by setting the 'Minimum' property of all parameters to the smallest recognizable positive value in MATLAB®, eps(0).

FileName      = 'ipfd';        % File describing the model structure.
Order         = [5 1 2];             % Model orders [ny nu nx]. 
%Vector [Ny Nu Nx], specifying the number of model outputs Ny, inputs Nu, 
%and states Nx.

w0 = 5.86;
g = 9.81;
dt = 0.001;
m = 39.3;
Parameters    = [g, dt, m];   % Initial parameters.
InitialStates = [w0; 0];              % Initial states.
Ts            = 0;                   % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts);
%Output [Fx, Fy, v, x, y]
set(nlgr, 'OutputName', {'Horizontal ground reaction force' 'Vertical ground reaction force' 'velocity' 'Horizontal position' 'Vertical position'}, ...
          'OutputUnit', {'N' 'N' 'm/s' 'm' 'm'},            ...
          'TimeUnit', 's'); %output
nlgr = setinit(nlgr, 'Name', {'Angular Velocity', 'Sway'});
nlgr = setinit(nlgr, 'Unit', {'rad/s' 'rad'});
nlgr = setpar(nlgr, 'Name', {'Gravity constant', 'Integration step', 'Mass'}); %parameters
nlgr = setpar(nlgr, 'Unit', {'m/s^2', 's', 'kg'}); %units
nlgr = setpar(nlgr, 'Minimum', {eps(0) eps(0) eps(0)});   % All parameters > 0.

%%Simulation
%Simulation before estimating parameters

%the output of sim() should be an iddata object with as many
%output channels as there are output parameters. The
%documentation implies that it would have no input channels

Y = sim(nlgr, leg);

figure
plot(Y(:,1,1));  %plot Fh
title('plot of Fh');

figure
plot(Y(:,1,2));  %plot Fv
title('plot of Fv');

[Walter Roberson]

You can't just run ipfd by itself; it has to be passed input parameters. It is the work routine called inside the greybox estimating.

[Paulo Costa]

When we are running the code you posted, we get this message:

Error using iddata (line 229)
The number of elements in the "OutputUnit"
property value must be equal to the number
of output channels. Use a cell array of 0
string(s).

Error in main (line 122)
leg = iddata([], r.', 0.001, 'Name', 'Leg
length', 'OutputUnit', {'m'});

The program won't let the code run because of the error, what should we do?

[Katarzyna Wieciorek]

We also consider simplification of the model but we would like to ask you first: what do you think, may it make anything easier when we are at this stage?

We have the same error as in the very beginning, r.' or r(:) give the same. I am confused.

In simplified model we could make r another constant parameter and have only one time-varying component, O. Will it be much easier and doable if current work seems impossible to finish? We are just little students with no more ideas what to do :( and we are very sorry for making you tired and confused. We can live without rkf45 but we can't without working script.

[Walter Roberson]

For leg, change OutputUnit to InputUnit

[Paulo Costa]

When we changed now it gives us this error

Error using idnlgrey/isvalid (line 93)
  Error or mismatch in the specified ODE
  file "ipfd".
  The error message is: "Too many input
  arguments."

    Error in idnlgrey/sim (line 102)
    error(isvalid(nlsys, 'OnlyFileName',
    data));

    Error in main (line 182)
    sim(nlgr, leg);% Y = sim(nlgr, leg);

[Walter Roberson]

You currently have Sway (O) and Angular Velocity (w) as states. In terms of the modeling, a state is something about the model that gets initialized and varies in time as the model runs, in response to conditions of the model, with the new values being determined based upon any other inputs and the old values. In your model, the Sway and Angular Velocity of one step depend upon the Sway and Angular Velocity of the preceding step. Anything that is computed based upon current conditions and is needed to know to compute the next conditions is "state".

Any external force or change being imposed upon the system is "input". You have indicated that Leg Length, r, should be input.

The pendulum equations say that the distance traveled increases as the leg length increases, and that makes some sense: clearly a person with very short legs would have to take more steps to keep up with someone with much longer legs.

People can only increase their leg length by wearing taller shoes or by lifting up on to the ball of their foot. A quick estimate is that the length gained by lifting up their arch might be on the order of 15%. But one must keep in mind that when people wear high-heels (which lift the arch and make the leg longer), they tend to walk more slowly not faster. Therefore leg length alone should not be considered as sufficient input for determining walking (or running) speed. It might be good enough for a demonstration model but not for serious work.

For serious gait analysis work, you would need to model the effects of muscles, possibly as springs; including that as people walk (and especially as they run) the mechanics of motion stores energy in the bones and muscles of the forward foot and springs back against the ground; see http://www.sciencedaily.com/releases/2015/04/150401132742.htm

But I digress.

[Katarzyna Wieciorek]

Walter,

I am afraid you didn't understand us well :(

We model human (exactly kids) gait basing on paper where all equations are given. You can see them in our script in comments %%Pendulum modelling. A leg is treated as pendulum and its length changes with time. O and r are time varying and input. We load r, not w.

My question was: if we make r parameter not dependent on time, having only time-varying O will it be better?

Every model is simplification, our paper doesn't care about muscles. Maybe we could simplify even more and get rid of time-varying r?

[Walter Roberson]

Your code has Sway (O) as state, not as input. Look at your comments in the model:

    %a(i) = (-g*cos(O(i))/r
    %w(i+1) = w(i) + a(i)*dt
    %O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2

The current Sway, O(i) is examined and an acceleration is calculated from it. The acceleration times the timestep is added to the current angular velocity to calculate the updated angular velocity, w(i+1). The acceleration and current angular velocity are used to calculate an updated Sway, O(i+1). Then forces and velocities are calculated based upon those. The updated Sway and Angular velocity are output as state in the dx output variable. Those two values, the updated O and updated w, then become what is passed in the next time the routine is called: they are remembered, and form part of the state.

You can use leg-length as an input, with the physical meaning being questionable if length varies more than about 15% (going up on the ball of the foot.)

If you use Sway as an input. However, if the input Sway at time point i, O(i), varies too much from the natural state calculated out of time point i-1, then you could have the effect of a whole series of pushes in a very jarring manner. But that's something for you to examine.

I will provide two versions of ipfd, one for the situation where O is not an Input (because it is a State), and the other for the case where it is an Input. I need to do some housekeeping now, though, so it will have to be later tonight.

[Walter Roberson]

Here is the ipfd version for having r as an Input, but Sway still as State. It is fairly similar to what I posted before, but when I looked back I saw that I made two small mistakes in what I had then. Also, the documentation specifies an additional parameter of uncertain meaning. I also renamed the input variables to match the template documentation, to make it easier to compare the file against the template.

function [dx, Y] = ipfd(t, x, u, p1, p2, p3, varargin)
    %A pendulum system

    %Model:
    %a(i) = (-g*cos(O(i))/r
    %w(i+1) = w(i) + a(i)*dt
    %O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
    %Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
    %Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
    %v = w*r*sin(O)

    %states
    wi = x(1);
    Oi = x(2);

    %inputs
    r = u(1);

    %parameters
    g = p1;
    dt = p2;
    m = p3;

    %State equations
    a = (-g*cos(Oi))/r; %angular acceleration, not output
    w = wi + a*dt; %Angular velocity  %Fixed
    O = Oi + wi*dt + 1/2*a*(dt.^2);   %Fixed

    %Output equations
    Fx = m.*((diff(r, 2)-r.*w.^2)*cos(O)-(r.*a+2.*diff(r, 1)*w)*sin(O));
    Fy = m.*((diff(r, 2)-r.*w.^2)*sin(O)-(r.*a+2.*diff(r, 1)*w)*cos(O)) + m.*g;
    v = w*r*sin(O); 
    x = r*cos(O);
    y = r*sin(O);
    dx = [w;O];
    Y = [Fx; Fy; v; x; y];    %output
  end

This goes together with this driver. I fixed the InputUnit problem. As well I had the Parameter as the wrong vector shape, which might have caused the "too many parameter" problem, perhaps (but I don't really know what caused that.)

     %%An inverted pendulum model directly applied to normal human gait
    %Model predicts values of ground reaction forces and velocity of human gait.
    %An inverted pendulum is defined basing on the average center-of-pressure 
    %to the instantaneous center-of-mass and equations of motion during single 
    %support that allowed a telescoping action were derived.

%%Data

% load('O');
% O = iddata(O.', [], 0.001, 'Name', 'Loaded O', 'InputName', {'Sway'}, 'InputUnit', {'rad'});
% set(sway, 'Tstart', 0, 'TimeUnit', 's');

load('r');
leg = iddata([], r.', 0.001, 'Name', 'Loaded r', 'InputName', {'Leg Length'}, 'InputUnit', {'m'});
set(leg, 'Tstart', 0, 'TimeUnit', 's');

%%Pendulum Modeling
% This paragraph presents a model structure for the pendulum system.
% The forward and inverse dynamics of it was studied by the same authors 
%before. The structure is as follows:
%
%a(i) = (-g*cos(O(i))/r
%w(i+1) = w(i) + a(i)*dt
%O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
%Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
%Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
%v = w*r*sin(O)
%
% having parameters (or constants):
%    a - angular acceleration [rad/s^2]
%    g - the gravity constant [m/s^2]
%    O - sway [rad]
%    r - the length of the leg [m]
%    w - angular velocity [rad/s]
%    dt - integration steps [s]
%    m - the weight [kg]
%    Fh - horizontal ground reaction force [
%    Fv - vertical ground reaction force
%    v - velocity

%We enter this information into the MATLAB file ipfd.m.

The next step is to create a generic IDNLGREY object for describing the pendulum. We also enter information about the names and the units of the inputs, states, outputs and parameters. Owing to the physical reality all model parameters must be positive and this is imposed by setting the 'Minimum' property of all parameters to the smallest recognizable positive value in MATLAB®, eps(0).

FileName      = 'ipfd';        % File describing the model structure.
Order         = [5 1 2];             % Model orders [ny nu nx]. 
%Vector [Ny Nu Nx], specifying the number of model outputs Ny, inputs Nu, 
%and states Nx.

w0 = 5.86;
g = 9.81;
dt = 0.001;
m = 39.3;
Parameters    = [g; dt; m];   % Initial parameters.  %CHANGED
InitialStates = [w0; 0];              % Initial states.
Ts            = 0;                   % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts);
%Output [Fx, Fy, v, x, y]
set(nlgr, 'OutputName', {'Horizontal ground reaction force' 'Vertical ground reaction force' 'velocity' 'Horizontal position' 'Vertical position'}, ...
          'OutputUnit', {'N' 'N' 'm/s' 'm' 'm'},            ...
          'TimeUnit', 's'); %output
nlgr = setinit(nlgr, 'Name', {'Angular Velocity', 'Sway'});
nlgr = setinit(nlgr, 'Unit', {'rad/s' 'rad'});
nlgr = setpar(nlgr, 'Name', {'Gravity constant', 'Integration step', 'Mass'}); %parameters
nlgr = setpar(nlgr, 'Unit', {'m/s^2', 's', 'kg'}); %units
nlgr = setpar(nlgr, 'Minimum', {eps(0) eps(0) eps(0)});   % All parameters > 0.

%%Simulation
%Simulation before estimating parameters

%the output of sim() should be an iddata object with as many
%output channels as there are output parameters. The
%documentation implies that it would have no input channels

%Y = sim(nlgr, leg);
sim(nlgr, leg);

%{
 %block commented out
  figure
  plot(Y(:,1,1));  %plot Fh
  title('plot of Fh');

    figure
    plot(Y(:,1,2));  %plot Fv
    title('plot of Fv');
  %}

[Walter Roberson]

Here is the ipfd variation having Sway and r as a Input, with Sway no longer a State.

function [dx, Y] = ipfd(t, x, u, p1, p2, p3, varargin)
    %A pendulum system

    %Model:
    %a(i) = (-g*cos(O(i))/r
    %w(i+1) = w(i) + a(i)*dt
    %O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
    %Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
    %Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
    %v = w*r*sin(O)

    %states
    wi = x(1);

    %inputs
    r = u(1);
    Oi = u(2);  %CHANGES

    %parameters
    g = p1;
    dt = p2;
    m = p3;

    %State equations
    a = (-g*cos(Oi))/r; %angular acceleration, not output
    w = wi + a*dt; %Angular velocity  %Fixed
    O = Oi + wi*dt + 1/2*a*(dt.^2);   %Fixed

    %Output equations
    Fx = m.*((diff(r, 2)-r.*w.^2)*cos(O)-(r.*a+2.*diff(r, 1)*w)*sin(O));
    Fy = m.*((diff(r, 2)-r.*w.^2)*sin(O)-(r.*a+2.*diff(r, 1)*w)*cos(O)) + m.*g;
    v = w*r*sin(O); 
    x = r*cos(O);
    y = r*sin(O);
    dx = [w];                 %changes, O not a state
    Y = [Fx; Fy; v; x; y];    %output
  end

This goes together with the below driver.

     %%An inverted pendulum model directly applied to normal human gait
    %Model predicts values of ground reaction forces and velocity of human gait.
    %An inverted pendulum is defined basing on the average center-of-pressure 
    %to the instantaneous center-of-mass and equations of motion during single 
    %support that allowed a telescoping action were derived.

%%Data

load('O');
% O = iddata(O.', [], 0.001, 'Name', 'Loaded O', 'InputName', {'Sway'}, 'InputUnit', {'rad'});
% set(sway, 'Tstart', 0, 'TimeUnit', 's');

load('r');
% leg = iddata([], r.', 0.001, 'Name', 'Loaded r', 'InputName', {'Leg Length'}, 'InputUnit', {'m'});  %CHANGED
% set(leg, 'Tstart', 0, 'TimeUnit', 's');

leg_sway = iddata([], [r(:), O(:)], 0.001, 'Name', 'Loaded r O', 'InputName', {'Leg Length', 'Sway'}, 'InputUnit', {'m', 'rad'});

%%Pendulum Modeling
% This paragraph presents a model structure for the pendulum system.
% The forward and inverse dynamics of it was studied by the same authors 
%before. The structure is as follows:
%
%a(i) = (-g*cos(O(i))/r
%w(i+1) = w(i) + a(i)*dt
%O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
%Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
%Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
%v = w*r*sin(O)
%
% having parameters (or constants):
%    a - angular acceleration [rad/s^2]
%    g - the gravity constant [m/s^2]
%    O - sway [rad]
%    r - the length of the leg [m]
%    w - angular velocity [rad/s]
%    dt - integration steps [s]
%    m - the weight [kg]
%    Fh - horizontal ground reaction force [
%    Fv - vertical ground reaction force
%    v - velocity

%We enter this information into the MATLAB file ipfd.m.

The next step is to create a generic IDNLGREY object for describing the pendulum. We also enter information about the names and the units of the inputs, states, outputs and parameters. Owing to the physical reality all model parameters must be positive and this is imposed by setting the 'Minimum' property of all parameters to the smallest recognizable positive value in MATLAB®, eps(0).

FileName      = 'ipfd';        % File describing the model structure.
Order         = [5 2 1];             % Model orders [ny nu nx]. 
%Vector [Ny Nu Nx], specifying the number of model outputs Ny, inputs Nu, 
%and states Nx.

w0 = 5.86;
g = 9.81;
dt = 0.001;
m = 39.3;
Parameters    = [g; dt; m];   % Initial parameters.  %CHANGED
InitialStates = [w0];              % Initial states.
Ts            = 0;                   % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts);
%Output [Fx, Fy, v, x, y]
set(nlgr, 'OutputName', {'Horizontal ground reaction force' 'Vertical ground reaction force' 'velocity' 'Horizontal position' 'Vertical position'}, ...
          'OutputUnit', {'N' 'N' 'm/s' 'm' 'm'},            ...
          'TimeUnit', 's'); %output
nlgr = setinit(nlgr, 'Name', {'Angular Velocity'});
nlgr = setinit(nlgr, 'Unit', {'rad/s'});
nlgr = setpar(nlgr, 'Name', {'Gravity constant', 'Integration step', 'Mass'}); %parameters
nlgr = setpar(nlgr, 'Unit', {'m/s^2', 's', 'kg'}); %units
nlgr = setpar(nlgr, 'Minimum', {eps(0) eps(0) eps(0)});   % All parameters > 0.

%%Simulation
%Simulation before estimating parameters

%the output of sim() should be an iddata object with as many
%output channels as there are output parameters. The
%documentation implies that it would have no input channels

%Y = sim(nlgr, leg);
sim(nlgr, leg);

%{
 %block commented out
  figure
  plot(Y(:,1,1));  %plot Fh
  title('plot of Fh');

    figure
    plot(Y(:,1,2));  %plot Fv
    title('plot of Fv');
  %}

[Paulo Costa]

After trying the first one (Sway as a State):

>> main
Error using idnlgrey/isvalid (line 93)
The second output argument of the ODE file "ipfd" must contain 5 elements. Type "help
idnlgrey" for more information.

Error in idnlgrey/sim (line 102)
error(isvalid(nlsys, 'OnlyFileName', data));

Error in main (line 260)
sim(nlgr, leg);

After trying the second one (Sway as input):

Error using idnlgrey/sim (line 100)
Simulation data and model sizes are incompatible. For simulation without noise, data should
have 2 inputs. For simulation with noise, data should have 7 inputs.

Error in main (line 342)
sim(nlgr, leg);

[Walter Roberson]

For the second one, the call to sim() should be

sim(nlgr, leg_sway);

and the second one should not even have a variable "leg" defined so if it made it that far you must have changed the code.

[Walter Roberson]

I see what is happening in the first situation; I understand why it is saying that there needs to be 5 outputs.

The code there calculates some output values in terms of the first and second derivatives of r, r' and r''. The code attempts to implement those by using diff(r,1) and diff(r,2). However, diff() applied to a numeric object is not derivative: it is the difference function.

diff(x) = [x(2)-x(1), x(3)-x(2), x(4)-x(3), ... x(n)-x(n-1)]]

When applied to a scalar, there is no second element to take the difference against, so the result is the empty vector. When r is a scalar, diff(r,1) is [] and diff(r,2) is [] as well. So the variables calculated in terms of those expressions come out empty, and that makes the overall list Y come out with fewer than 5 elements.

I think I can see a way to deal with it, but it requires that I impose an artificial boundary condition such as initializing r' and r'' to 0 at the beginning. I need more rest before I code that.

[Paulo Costa]

In version 1 it still gives us the same problem basicly:

Error using idnlgrey/isvalid (line 93)
The second output argument of the ODE file
"ipfd_1" must contain 5 elements. Type "help
idnlgrey" for more information.

Error in idnlgrey/sim (line 102)
error(isvalid(nlsys, 'OnlyFileName', data));

Error in main_1 (line 71)
sim(nlgr, leg);

In version 2 I saw that it should be leg_sway, but it still gives us the 5 element error

Error using idnlgrey/isvalid (line 93)
The second output argument of the ODE file
"ipfd_2" must contain 5 elements. Type "help
idnlgrey" for more information.

Error in idnlgrey/sim (line 102)
error(isvalid(nlsys, 'OnlyFileName', data));

Error in main_2 (line 72)
sim(nlgr, leg_sway);

[Walter Roberson]

The error about needing to contain 5 elements is due to the issue with diff(r,1) and diff(r,2) resulting in empty matrices because diff() applied to numeric values is not the differential operator. diff() is only the differential operator when diff() is being applied to a symbolic object.

I think I know how to deal with it, but I am still waking up.

[Paulo Costa]

Walter, we saw what you said and we changed it. Instead of usig de diff() function we used the ' or D() wich is the operator.

I've used version 1 of your code and it simulates what we want, but when we try to plot it we get:

Error using iddata/subsref (line 50)
Input index exceeds number of inputs.

Error in main_1 (line 75)
  plot(Y(:,1,1));  %plot Fh

What do you think it is now?

[Walter Roberson]

No, ' is just complex conjugate transpose, and D() is not defined.

[Walter Roberson]

Sway as state version. I modified it to treat handle r' and r'' through states

function [dx, Y] = ipfd(t, x, u, p1, p2, p3, varargin)
    %A pendulum system

    %Model:
    %a(i) = (-g*cos(O(i))/r
    %w(i+1) = w(i) + a(i)*dt
    %O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
    %Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
    %Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
    %v = w*r*sin(O)

    %states
    wi = x(1);
    Oi = x(2);
    old_r = x(3); %previous r
    old_dr = x(4); %previous r'

    %inputs
    r = u(1);

    %parameters
    g = p1;
    dt = p2;
    m = p3;

    %State equations
    a = (-g .* cos(Oi)) ./ r; %angular acceleration, not output
    w = wi + a .* dt; %Angular velocity
    O = Oi + wi .* dt + 1/2 .* a .* (dt.^2);

    dr = (r - old_r) ./ dt;
    ddr = (dr - old_dr) ./ dt;

    %Output equations

    mrr = (ddr - r .* w.^2);
    rarw = (r .* a + 2 .* dr .* w);
    Fx = m .* (mrr .* cos(O) - rarw .* sin(O));
    Fy = m .* (mrr .* sin(O) - rarw .* cos(O)) + m .* g;
    v = w*r*sin(O);
    x = r*cos(O);
    y = r*sin(O);

    dx = [w; O; r; dr];
    Y = [Fx; Fy; v; x; y];
  end

Corresponding driver:

%%An inverted pendulum model directly applied to normal human gait
%Model predicts values of ground reaction forces and velocity of human gait.
%An inverted pendulum is defined basing on the average center-of-pressure 
%to the instantaneous center-of-mass and equations of motion during single 
%support that allowed a telescoping action were derived.

%%Data

% load('O');
% O = iddata(O.', [], 0.001, 'Name', 'Loaded O', 'InputName', {'Sway'}, 'InputUnit', {'rad'});
% set(sway, 'Tstart', 0, 'TimeUnit', 's');

load('r');
leg = iddata([], r.', 0.001, 'Name', 'Loaded r', 'InputName', {'Leg Length'}, 'InputUnit', {'m'});
set(leg, 'Tstart', 0, 'TimeUnit', 's');

%%Pendulum Modeling
% This paragraph presents a model structure for the pendulum system.
% The forward and inverse dynamics of it was studied by the same authors 
%before. The structure is as follows:
%
%a(i) = (-g*cos(O(i))/r
%w(i+1) = w(i) + a(i)*dt
%O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
%Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
%Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
%v = w*r*sin(O)
%
% having parameters (or constants):
%    a - angular acceleration [rad/s^2]
%    g - the gravity constant [m/s^2]
%    O - sway [rad]
%    r - the length of the leg [m]
%    w - angular velocity [rad/s]
%    dt - integration steps [s]
%    m - the weight [kg]
%    Fh - horizontal ground reaction force [
%    Fv - vertical ground reaction force
%    v - velocity

%We enter this information into the MATLAB file ipfd.m.

The next step is to create a generic IDNLGREY object for describing the pendulum. We also enter information about the names and the units of the inputs, states, outputs and parameters. Owing to the physical reality all model parameters must be positive and this is imposed by setting the 'Minimum' property of all parameters to the smallest recognizable positive value in MATLAB®, eps(0).

FileName      = 'ipfd';        % File describing the model structure.
Order         = [5 1 4];             % Model orders [ny nu nx]. 
%Vector [Ny Nu Nx], specifying the number of model outputs Ny, inputs Nu, 
%and states Nx.

w0 = 5.86;
g = 9.81;
dt = 0.001;
m = 39.3;
Parameters    = [g; dt; m];   % Initial parameters.  %CHANGED
InitialStates = [w0; 0, r(1), 0];              % Initial states.
Ts            = 0;                   % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts);
%Output [Fx, Fy, v, x, y]
set(nlgr, 'OutputName', {'Horizontal ground reaction force' 'Vertical ground reaction force' 'velocity' 'Horizontal position' 'Vertical position'}, ...
          'OutputUnit', {'N' 'N' 'm/s' 'm' 'm'},            ...
          'TimeUnit', 's'); %output
nlgr = setinit(nlgr, 'Name', {'Angular Velocity', 'Sway', 'Old r', 'Old dr'});
nlgr = setinit(nlgr, 'Unit', {'rad/s' 'rad', 'm', 'm/s'});
nlgr = setpar(nlgr, 'Name', {'Gravity constant', 'Integration step', 'Mass'}); %parameters
nlgr = setpar(nlgr, 'Unit', {'m/s^2', 's', 'kg'}); %units
nlgr = setpar(nlgr, 'Minimum', {eps(0) eps(0) eps(0)});   % All parameters > 0.

%%Simulation
%Simulation before estimating parameters

%the output of sim() should be an iddata object with as many
%output channels as there are output parameters. The
%documentation implies that it would have no input channels

%Y = sim(nlgr, leg);
sim(nlgr, leg);

%{
%block commented out

figure
plot(Y(:,1,1));  %plot Fh
title('plot of Fh');

    figure
    plot(Y(:,1,2));  %plot Fv
    title('plot of Fv');
  %}

[Walter Roberson]

Sway as input version.

function [dx, Y] = ipfd(t, x, u, p1, p2, p3, varargin)
    %A pendulum system

    %Model:
    %a(i) = (-g*cos(O(i))/r
    %w(i+1) = w(i) + a(i)*dt
    %O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
    %Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
    %Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
    %v = w*r*sin(O)

    %states
    wi = x(1);
    old_r = x(2); %previous r
    old_dr = x(3); %previous r'

    %inputs
    r = u(1);
    Oi = x(2);

    %parameters
    g = p1;
    dt = p2;
    m = p3;

    %State equations
    a = (-g .* cos(Oi)) ./ r; %angular acceleration, not output
    w = wi + a .* dt; %Angular velocity
    O = Oi + wi .* dt + 1/2 .* a .* (dt.^2);

    dr = (r - old_r) ./ dt;
    ddr = (dr - old_dr) ./ dt;

    %Output equations

    mrr = (ddr - r .* w.^2);
    rarw = (r .* a + 2 .* dr .* w);
    Fx = m .* (mrr .* cos(O) - rarw .* sin(O));
    Fy = m .* (mrr .* sin(O) - rarw .* cos(O)) + m .* g;
    v = w*r*sin(O);
    x = r*cos(O);
    y = r*sin(O);

    dx = [w; r; dr];
    Y = [Fx; Fy; v; x; y];
  end

Corresponding driver:

%%An inverted pendulum model directly applied to normal human gait
%Model predicts values of ground reaction forces and velocity of human gait.
%An inverted pendulum is defined basing on the average center-of-pressure 
%to the instantaneous center-of-mass and equations of motion during single 
%support that allowed a telescoping action were derived.

%%Data

load('O');
% O = iddata(O.', [], 0.001, 'Name', 'Loaded O', 'InputName', {'Sway'}, 'InputUnit', {'rad'});
% set(sway, 'Tstart', 0, 'TimeUnit', 's');

load('r');
% leg = iddata([], r.', 0.001, 'Name', 'Loaded r', 'InputName', {'Leg Length'}, 'InputUnit', {'m'});
% set(leg, 'Tstart', 0, 'TimeUnit', 's');

   leg_sway = iddata([], [r(:), O(:)], 0.001, 'Name', 'Loaded r O', 'InputName', {'Leg Length', 'Sway'}, 'InputUnit', {'m', 'rad'});
  set(leg_sway, 'Tstart', 0, 'TimeUnit', 's');

%%Pendulum Modeling
% This paragraph presents a model structure for the pendulum system.
% The forward and inverse dynamics of it was studied by the same authors 
%before. The structure is as follows:
%
%a(i) = (-g*cos(O(i))/r
%w(i+1) = w(i) + a(i)*dt
%O(i+1) = O(i) + w(i)*dt + 1/2*a(i)*dt^2
%Fh = m*[r''-r*w^2)*cos(O) - (r*a+2*r'*w)*sin(O))]
%Fv = m*[r''-r*w^2)*sin(O) - (r*a+2*r'*w)*cos(O))] + m*g
%v = w*r*sin(O)
%
% having parameters (or constants):
%    a - angular acceleration [rad/s^2]
%    g - the gravity constant [m/s^2]
%    O - sway [rad]
%    r - the length of the leg [m]
%    w - angular velocity [rad/s]
%    dt - integration steps [s]
%    m - the weight [kg]
%    Fh - horizontal ground reaction force [
%    Fv - vertical ground reaction force
%    v - velocity

%We enter this information into the MATLAB file ipfd.m.

The next step is to create a generic IDNLGREY object for describing the pendulum. We also enter information about the names and the units of the inputs, states, outputs and parameters. Owing to the physical reality all model parameters must be positive and this is imposed by setting the 'Minimum' property of all parameters to the smallest recognizable positive value in MATLAB®, eps(0).

FileName      = 'ipfd';        % File describing the model structure.
Order         = [5 2 3];             % Model orders [ny nu nx]. 
%Vector [Ny Nu Nx], specifying the number of model outputs Ny, inputs Nu, 
%and states Nx.

w0 = 5.86;
g = 9.81;
dt = 0.001;
m = 39.3;
Parameters    = [g; dt; m];   % Initial parameters.  %CHANGED
InitialStates = [w0; r(1), 0];              % Initial states.
Ts            = 0;                   % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts);
%Output [Fx, Fy, v, x, y]
set(nlgr, 'OutputName', {'Horizontal ground reaction force' 'Vertical ground reaction force' 'velocity' 'Horizontal position' 'Vertical position'}, ...
          'OutputUnit', {'N' 'N' 'm/s' 'm' 'm'},            ...
          'TimeUnit', 's'); %output
nlgr = setinit(nlgr, 'Name', {'Angular Velocity', 'Old r', 'Old dr'});
nlgr = setinit(nlgr, 'Unit', {'rad/s', 'm', 'm/s'});
nlgr = setpar(nlgr, 'Name', {'Gravity constant', 'Integration step', 'Mass'}); %parameters
nlgr = setpar(nlgr, 'Unit', {'m/s^2', 's', 'kg'}); %units
nlgr = setpar(nlgr, 'Minimum', {eps(0) eps(0) eps(0)});   % All parameters > 0.

%%Simulation
%Simulation before estimating parameters

%the output of sim() should be an iddata object with as many
%output channels as there are output parameters. The
%documentation implies that it would have no input channels

%Y = sim(nlgr, leg_sway);
sim(nlgr, leg_sway);

%{
%block commented out

figure
plot(Y(:,1,1));  %plot Fh
title('plot of Fh');

    figure
    plot(Y(:,1,2));  %plot Fv
    title('plot of Fv');
  %}

[Walter Roberson]

I have not figured out plotting. That's why I left it as sim() with no outputs (which causes plotting) and commented out the rest.