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');
Best Answer