using constraints in ode45

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);

Answer 1

[Matt J]

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);

[Alex]

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]

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?