using constraints in ode45

Question

Alex 2014-12-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/167527-using-constraints-in-ode45

编辑： Matt J 2014-12-22

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2014-12-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/167527-using-constraints-in-ode45#answer_163058

编辑：Matt J 2014-12-22

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Alex 2014-12-22

编辑：Alex 2014-12-22

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 2014-12-22

编辑：Matt J 2014-12-22

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?

请先登录，再进行评论。