using constraints in ode45

7 views (last 30 days)
Alex
Alex on 22 Dec 2014
Edited: Matt J on 22 Dec 2014
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);

Answers (1)

Matt J
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
Alex
Alex on 22 Dec 2014
Edited: Alex on 22 Dec 2014
no, p is the vector of parameters being optimized. the concentrations are used in the objective function (least squares minimization), and also calculated as unretained intermediate values in the call to ode45
Matt J
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?

Sign in to comment.

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!