Nonlinear Grey Box Modeling code - really slow or not working?
10 views (last 30 days)
Show older comments
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');
0 Comments
Accepted Answer
Rajiv Singh
on 28 Jul 2011
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.
More Answers (1)
Giovanny Mateus
on 11 Jan 2018
Edited: Giovanny Mateus
on 11 Jan 2018
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
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!