System Identification for epidemic modeling

I am trying to estimate parameters using the system identification toolbox. I used the example of three ecological systems provided in the documentation, and worked on it o build my own model.

Here is my model:

function [dx, y] = shigella(t, u, x,a,b,q,m,r,varargin)

y = [x(1); x(2); x(3)];

dx = [q - m*x(1) - b*x(1)*x(2) + a*x(3); (b*x(1)) - (m+r)*x(2); (r*x(2)) - (m+a)*x(3)];

-----------------------------------------------------------------------------------------------

And the system Id code:

FileName      = 'shigella';              % File describing the model structure.
Order         = [3 0 3];                % Model orders [ny nu nx].
Parameters    = struct('Name',    {'Loss of immunity', 'infectivity rate', 'population renewal', 'death rate', 'recovery rate'}, ...
   'Unit',    {'1/week' '1/week' '1/week' '1/week', '1/week'},                    ...
   'Value',   {0.25 0.002 10 0.012 0.14},                                        ...
   'Minimum', {-Inf -Inf -Inf -Inf -Inf},                                    ...
   'Maximum', {Inf Inf Inf Inf Inf},                                        ...
   'Fixed',   {false false false false false});   % Estimate all 4 parameters.
InitialStates = struct('Name', {'Susceptible', 'Infected', 'Recovered'}, ...
   'Unit',    {'Size (in thousands)' 'Size (in thousands)', 'Size(inthousands)'},     ...
   'Value',   {20 5 10},                                         ...
   'Minimum', {0 0 0},                                             ...
   'Maximum', {Inf Inf Inf},                                         ...
   'Fixed',   {false false false});   % Estimate both initial states.
Ts            = 0;                      % Time-continuous system.
nlgr = idnlgrey(FileName,Order, Parameters, InitialStates, Ts,                  ...
   'Name', 'Population Variation',               ...
   'OutputName', {'Susceptible', 'Infected', 'Recovered'}, ...
   'OutputUnit', {'Size (in thousands)' 'Size (in thousands)', 'Size (in thousands)'},     ...
   'TimeUnit', 'week');
present(nlgr)

w = [1000 5 0];
u = [80 15 3];
v = [65 24 10];
f = [90 27 8];
z = [87 20 10];
g = [w; u; v; f; z];

z = iddata(g,[], 1, 'Name', 'Population Variation');
set(z, 'OutputName', {'Susceptible', 'Infected', 'Recovered'}, ...
   'Tstart', 0, 'TimeUnit', 'Week');

compare(z, nlgr, 1)

opt = nlgreyestOptions;
opt.Display = 'on';
opt.SearchOption.MaxIter = 50;
nlgr = nlgreyest(nlgr, opt);

disp('   True      Estimated parameter vector');
ptrue = [0.25 0.002 10 0.012 0.14];
fprintf('   %6.3f    %6.3f\n', [ptrue'; getpvec(nlgr)']);
disp(' ');
disp('   True      Estimated initial states');
x0true = [20 5 10];
fprintf('   %6.3f    %6.3f\n', [x0true'; cell2mat(getinit(nlgr, 'Value'))']);

--------------------------------------------------------------------------------------------

When I run the code part by part, i get the following errors:

??? Error using ==> idnlgrey.isvalid at 165 Error or mismatch in the specified ODE file 'shigella'. The error message is: 'Attempted to access x(1); index out of bounds because numel(x)=0.'

I tried to get rid of the compare part and move on, but i got this error

??? Undefined function or variable 'nlgreyestOptions'.

---------------------------------------------------------------------------

I would really appreciate any help I can get on this.

[Ojaswita]

Any possible help, please?

[Star Strider]

Just looking at your code, I find it difficult to understand your ‘shigella’ ODE. Before you use it with the System Identification Toolbox, be certain your ODE works with ode45 (or whatever solver is appropriate for it). See the documentation for ode45 to use the correct syntax for your ‘shigella’ function.

[Ojaswita]

Thanks alot! I also thought that the DEs are not being solved. But I was trying to follow the syntax that was given in the documentation. Since I was not sure about some steps, i tried to follow step by step as it said. If I need to change it, how do i solve it? Should I Solve it using ode45 in the same file?

Pls let me know

[Ojaswita]

Is there a possibility of that the functions work on later versions of MATLAB?