using constraints in ode45
7 views (last 30 days)
Show older comments
I'm solving a system of ODEs for parameter values that best fit a given data set (reaction kinetics). The concentrations are given as fractions and thus should always sum to one, but given the extra degree of separation between fmincon, I'm not sure how to implement that constraint. Because C is only calculated within ode45, there just might not be a simple way around this, and the parameter estimates this code arrives at are the best I can do
function [K EXITFLAG] = parameterEstimation(p0,t,C,LB,UB)
%function to be called to solve problem. p0 is the initial guess at the vector of parameters to be optimized, p. t and C are data vectors, which will be passed to the nested error function: this is what will be passed to fmincon
function error = objfun(p)
[t est] = ode45(@(t,C) ode(t,C,p),t,[1 0 0]);
error=sum((est(:,1) - C(:,1).^2) + (est(:,2) - C(:,2)).^2 + (est(:,3)...
- C(:,3)).^2);
end
options = optimset('Algorithm','interior-point','TolFun',10e-10);
[X, FVAL, EXITFLAG] = fmincon(@objfun,p0,[],[],[],[],LB,UB,[],options);
K = X;
end
function dCdt = ode(t,conc,p)
%system of ode's
dCdt = zeros(3,1);
dCdt(1) = -(conc(1) - conc(2)/p(1)) - 100*(conc(1) - conc(3)/p(2));
dCdt(2) = (conc(1) - conc(2)/p(1));
dCdt(3) = 100*(conc(1) - conc(3)/p(2));
end
t = linspace(0,15,16);
Ca = [1 0.248 0.180 0.152 0.113 0.113 0.102 0.0937 0.0799 0.0862 0.0825 ...
0.089 0.0832 0.0854 0.0796 0.0791]';
Cd = [0 0.316 0.509 0.622 0.697 0.735 0.775 0.767 0.773 0.779 0.781 0.794 ...
0.796 0.794 0.808 0.804]';
Cu = [0 0.384 0.274 0.222 0.174 0.159 0.14 0.147 0.14 0.114 0.129 0.116 ...
0.119 0.122 0.131 0.118]';
C = [Ca Cd Cu]; p0 = [100 1]; LB = 0; UB = 10^4;
[p, EXITFLAG] = parameterEstimation(p0,t,C,LB,UB);
0 Comments
Answers (1)
Matt J
on 22 Dec 2014
Edited: Matt J
on 22 Dec 2014
Are we to understand that the unknown vector 'p' are the concentrations? If so, it is a linear equality constraint. Use the Aeq and beq arguments,
[X, FVAL, EXITFLAG] = fmincon(@objfun,p0,[],[],Aeq,beq,LB,UB,[],options);
2 Comments
Matt J
on 22 Dec 2014
Edited: Matt J
on 22 Dec 2014
Because sum(dCdt)=0, the constraint should be satisfied automatically as long as you supply initial conditions satisfying sum(C)=1, which you appear to be doing. Have you actually checked the output of ode45 in isolation, without combining it with fmincon? Does sum(C,2) really deviate from 1 by more than what floating point errors would account for?
See Also
Categories
Find more on Optimization Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!