MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn moreOpportunities for recent engineering grads.

Apply Today
Asked by Hill on 16 Sep 2013

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 by Tom Lane on 17 Sep 2013

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

## 0 Comments