Nonlinear Grey Box Modeling code - really slow or not working?
10 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
Rajiv Singh
le 28 Juil 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.
Plus de réponses (1)
Giovanny Mateus
le 11 Jan 2018
Modifié(e) : Giovanny Mateus
le 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 commentaires
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!