Nonlinear Grey Box Modeling code - really slow or not working?

Hi all, At a Matlab user's great suggestion (see here: http://mathworks.com/matlabcentral/answers/12469-estimating-parameters-for-model-of-3-differential-equations) I've tried to start using nonlinear grey box modeling to estimate model parameters. I've read a lot about it and done the tutorial and come up with code that doesn't produce errors and may be working, but it's been running for about 15 hours and I'm not sure if it's still computing or it's not really going to work. I'm also not totally sure where the output will be, or if I've done this correctly. The data that I'm comparing the model to is 20 days of data on algal biomass, algal quota, and concentration of available resource. Any tips would be much appreciated! Thank you!

___Model file (Droop_batch_m.m)______

function [dx, y] = droop_batch_m(t, x, u, p_max, K_p, u_Amax, Q_Amin, varagin)
%script for grey box modeling
     % Output equations.
   y = [x(1);                           ... % Algal biomass.
        x(2);                           ... % Quota.
        x(3)                            ... % Resource.
        ];

     % State equations.
     dx = [u_Amax*x(1)-((Q_Amin/x(2))*(u_Amax*x(1)));... % Algal biomass.
          ((p_max * x(3)) / (K_p + x(3)))-(x(2)*u_Amax)-(Q_Amin*u_Amax);... % Quota.
          x(1)*((p_max * x(3)) / (K_p + x(3)))  ... % Resource.
          ];

____Script to run the estimator (based off of tutorial idnylgreydemo3.m____

%%estimator

FileName      = 'droop_batch_m';              % File describing the model structure.
Order         = [3 0 3];                % Model orders [ny nu nx].
Parameters    = struct('Name',    {'Maximum intake rate of nutrient by algae' 'Half saturation constant of nutrient intake by algae'   ...
                                   'Apparent max growth rate of algae' 'Subsistence quota of algae'}, ...
                       'Unit',    {'micromole-P/(mg-C*day)' 'micromole-P/L' '1/day' 'micromole-P/mg-C'},                    ...
                       'Value',   {1 0.001 0.1 0.01},                                        ...
                       'Minimum', {0 0 0 0},                                    ...
                       'Maximum', {Inf Inf Inf Inf},                                        ...
                       'Fixed',   {false false false false});  
InitialStates = struct('Name',    {'Algal biomass' 'Quota' 'Resource'}, ...
                       'Unit',    {'mg-C/L' 'micromole-P/mg-C' 'micromole-P/L'},     ...
                       'Value',   {0.02 0.25 50},                                         ...
                       'Minimum', {0 0 0},                                             ...
                       'Maximum', {Inf Inf Inf},                                         ...
                       'Fixed',   {true true true});   
Ts            = 0;                      % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts,                  ...
                'Name', 'Droop batch model',               ...
                'OutputName', {'Algal biomass' 'Quota' 'Resource'}, ...
                'OutputUnit', {'mg-C/L' 'micromole-P/mg-C' 'micromole-P/L'},     ...
                'TimeUnit', 'day');

%%check model

present(nlgr);

%%Load data

load 'batch_model_data.csv' %load data for culture
z = iddata(batch_model_data, [], 0.01, 'Name', 'Droop batch model');
set(z, 'OutputName', {'Algal biomass' 'Quota' 'Resource'}, ...
       'Tstart', 0, 'TimeUnit', 'Day');

%%parameter estimation
nlgr = pem(z, nlgr, 'Display', 'Full');

Answer 1

[Rajiv Singh]

A common reason for slowness is that the system is slow to simulate. You can quickly assess this by running a simulation command on the initial model before starting estimations. For example, check how long "step(nlgr, z.Samp(end))" takes on your system. If you find that the simulation is slow, your system could be stiff. In that case setting the solver to 'ode15s' or 'ode23s' might help (see nlgr.Algorithm.SimulationOptions.Solver)

Usually it helps to study your ODEs using MATLAB's ODE solvers (does the simulation run without errors, how long it takes etc) before wrapping them into grey box framework for estimation.

If choosing stiff solvers does not help, try writing the ODE function as a MEX file. The demos would illustrate how to do this with little knowledge of MEX files. Using a MEX function for the ODE can significantly improve the performance. But before you do that, make sure that that your equations are well-posed and simulate successfully using the chosen initial parameter values and estimation input signal.

[LS]

Thank you again!!

Answer 2

[Giovanny Mateus]

Hi,

I have a similar Problem. My system can be solved with ODE15s in 30 seconds, but with function compare (from the non-lineargreybox toolbox) it could´t be solved apparently after more than 10 hours of running (Matlab status "Busy") .. I appreciate any suggestion