piecewise curve fit of a biexponential function

Here is the scenario, a stirred tank reactor that has constant input and output rates that will change at some time t to different input and output rates that will mimic the following curve:

clear;clc;
A=24;
B=6;
alpha=1.5;
beta=0.25;
t=0:.5:30;
Cp=A*exp(-alpha*t)+B*exp(-beta*t); %standard 2 compartment pk curve
semilogy(t,Cp,'o')

using these equation that can be optimized using least squares regression:

%from time 0 to t  
Cplin1=Ao*((ri-ro)*t+Vo).^(-ri./(ri-ro))*Vo.^(ro./(ri-ro)); 
%from time t to end
Cplin2=Ao1*((ri1-ro1)*t+Vo1).^(-ri1./(ri1-ro1))*Vo1.^(ro1./(ri1-ro1));

where: Ao and Vo in this example are 300 and 10, respectively; Ao1 and Vo1 are determined at some time t from Cplin1; the t in Cplin2 most likely will require a phase shift; and using a gridded search (because that is what I'm most comfortable doing, but I'd be fine with another method if you wrote the code) to find ro, ri, ro1, ri1 and t.

here is what I have, but it doesn't run:

clear;clc;
Ao=300;%amount of drug
Vo=10;%volume drug is in
t12=3;%half life of drug
k=log(2)/t12;
A=24;
B=6
alpha=1.5
beta=0.25
ri=1;%flow rate into stirred tank
ro=.5;%flow rate out of stirred tank
t=0:.5:30; %time
phaseTime=t;
Cp=A*exp(-alpha*t)+B*exp(-beta*t); %standard 2 compartment pk curve
semilogy(t,Cp,'o')
SSBest=10000 %dummy
Vo1=Vo %dummy?
Ao1=Ao %dummy?
for ri=10:1:20; %flow rate in to a stirred tank
  for ro=10:1:20; %flow rate out of a stirred tank
  for ri1=10:1:20;        
  for ro1=10:1:20;
      %for n=1:length(t);
      if (ri*t-ro*t+Vo)>0 %AND ri~=ro %if ri=ro then i'll be dividing by 0
      if (ri1*t-ro1*t+Vo1)>0  %&&ri1~=ro1
          Cplin1=Ao*((ri-ro)*phaseTime+Vo).^(-ri./(ri-ro))*Vo.^(ro./(ri-ro)); %this is the right equation, to look at concentration
          Ao1=Ao*((ri-ro)*phaseTime+Vo).^(-ri./(ri-ro))*Vo.^(ro./(ri-ro))*((ri-ro)*phaseTime+Vo);
          Vo1=((ri-ro)*phaseTime+Vo);
          Cplin2=Ao1*((ri1-ro1)*phaseTime+Vo1).^(-ri1./(ri1-ro1))*Vo1.^(ro1./(ri1-ro1));   
       Cplin=[Cplin1(1:4) Cplin2(5:end)]; 
       SS=sum((Cp-Cplin).^2); %i'm fitting the curve using least squares
         if (SS<SSBest) %keeps only the ri and ro that minimize SSE
         RinBest=ri;
         RoutBest=ro;
         Rin1Best=ri1;
         Rout1Best=ro1;
         Ao1Best=Ao1;
         Vo1Best=Vo1;
         CplinBest=Cplin;
         disp([ri,ro,ri1,ro1,SS])
         SSBest=SS;
         %end 
      end   
   end
      end
  end
  end
  end
end
Cplin=Ao*(RinBest*t-RoutBest*t+Vo).^(-RinBest./(RinBest-RoutBest))*Vo.^(RoutBest./(RinBest-RoutBest)); %generate the curve
hold on
plot(t,CplinBest) 
hold off

Answer 1

[Tom Lane]

You have a number of places where you are multiplying two vectors. Most likely you want element-wise multiplication (".*") rather than matrix multiplication ("*"). I noticed these two lines, at least:

Ao1=Ao*((ri-ro)*phaseTime+Vo).^(-ri./(ri-ro))*Vo.^(ro./(ri-ro)).*((ri-ro)*phaseTime+Vo);

Cplin2=Ao1.*((ri1-ro1)*phaseTime+Vo1).^(-ri1./(ri1-ro1)).*Vo1.^(ro1./(ri1-ro1));