Storing values in an array from function file called by ode15s solver

1 view (last 30 days)
In program lumenslumens1,the creation of empty arrays such as errb,errd need to store new values.
In program bookbook1,ode15s solver is used which calls lumenslumens1 program in order to solver the differential equations and the ode15s sovler is decremented by for loop.
At n=20
At first time (t=1)), the errb and errd has values.
At the second time of integration(t=2),the errb and errd got new values.But the problem is,the values of errb and errd(at the first time of integration t=1),,did'nt store in that array for second time of integration.
I have to store all the values of errb and errd to the empty array of errb and errd,in order to perform integration(ie,erintb and erintd).
% function file consisting of differential equations
function res=lumenslumens1(t,z,n,xnp1,xnm1,s,x1,k,u,errb,errd,L,V,R)
lm=ones(12,1);
% % % %Equations based on Tray
if n<20 && n~=10
dxbdt=0;
dmfdt=0;
dxfdt=0;
elseif n==1
dmfdt=0;
dxfdt=0;
elseif n==10
dmdt=0;
dxdt=0;
dxbdt=0;
elseif n==20
dmdt=0;
dxdt=0;
dxbdt=0;
dmfdt=0;
dxfdt=0;
else
dmdt=0;
end
%given data
nf=10;
Ro=128.01;
R=ones(500,1);
V=178.01*ones(500,1);
Vo=178.01;
F=100;
beta=0.1;
alpha=2;
zf=0.55;
mo=10;
mbo=100;
mdo=100;
xbo=0.98;
xdo=0.02;
%initial value allocation to specific variables
m=z(1);
x=z(2);
xd=z(3);
xb=z(4);
mf=z(5);
xf=z(6);
mlast=z(7);
xlast=z(8);
%feedback controllers
kcb=1000;
kcd=1000;
taud=5;
taub=1.25;
%VLE for yten
ylast=(alpha*xlast)/(1+(alpha-1)*xlast);
if k==1
xlast=u(n);
end
dxddt=V(k,1)*(ylast-xd)/mdo;
%error values by set point
if k==1 || k>1
errb(k,1)=0.02-xbo;
errd(k,1)=0.98-xd;
elseif n==1
errb(k,1)=0.02-xb;
errd(k,1)=0.98-xd;
end
for i=1:k
i;
ybb=errb(i,1);
ydd=errd(i,1);
lm(i,1)=i;
lim=lm(i,1);
erintb=trapz(lim,ybb,1);
erintd=trapz(lim,ydd,1);
i=i+1;
end
V(k,1)=Vo-kcb*(errb(k,1)+erintb/taub);
R(k,1)=Ro+kcd*(errd(k,1)+erintd/taub);
%component continuity equation for condeser and reflux drum
% % dxddt=V(s,1)*(yten-xd)/mdo;
% dxddt=V(k,1)*(ylast-xd)/mdo;
%material balance on condeser and reflux drum
D=V(k,1)-R(k,1);
%material balance on reboiler
if n==1
B=L(1)-V(k,1);
end
%The flow rate at initial time,
Lo=Ro+F;
%liquid Hydraulic relationship
L(n)=Lo+((m-mo)/beta);
L(20)=Lo+((mlast-mo)/beta);
%VlE for y and ynf
y=(alpha*x)/(1+(alpha-1)*x);
ynf=(alpha*xf)/(1+(alpha-1)*xf);
L(n)=Lo+((m-10)/beta);
L(20)=Lo+((mlast-mo)/beta);
%Material balance on the last tray
dmlastdt=R(k,1)-L(20);
dxlastdt=((R(k,1)*xd)-(xlast*R(k,1)))/mlast;
if n==10
% Material balance of feed tray
ynm1=(alpha*xnm1)/(1+(alpha-1)*xnm1);
dmfdt=L(nf+1)-L(nf)+F;
dxfdt=((L(nf+1)*xnp1(k,1))-(L(nf)*xf)+(V(k,1)*(ynm1-ynf))+(F*zf)-(xf*(L(nf+1)-L(nf)+F)))/mf;
elseif n<20
ynm1=(alpha*xnm1)/(1+(alpha-1)*xnm1);
% Material balance on the general trays
dmdt=L(n+1)-L(n);
dxdt=((L(n+1)*xnp1(k,1))-(L(n)*x+V(k,1)*(ynm1-y))-(x*(L(n+1)-L(n))))/m;
elseif n==1
% Component Continuity Equation on reboiler
yb=(alpha*xb)/((1+(alpha-1)*xb));
dxbdt=(L(1)*x1(k,1)-B*xb-V(k,1)*yb)/mbo;
end
res=[dmdt;dxdt;dxddt;dxbdt;dmfdt;dxfdt;dmlastdt;dxlastdt];
end
% bookbook1 file-ode15s is used
s=[];
tp=[];
to=1;
ts=1;
tf=10;
k=[19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0];
xnp1=zeros(20,1);
u=[0.035,0.05719,0.08885,0.1318,0.18622,0.24951,0.31618,0.37498,0.43391,0.47688,0.51526,0.56295,0.61896,0.68052,0.74345,0.80319,0.85603,0.89995,0.93458,0.96079];
xnm1=zeros(20,1);
x1=zeros(20,1);
ze=zeros(2,8);
errb=zeros(15,1);
errd=zeros(15,1);
L=zeros(1,21);
R=ones(13,1);
V=178.01*ones(500,1);
% v=zeros(200,8);
w=zeros(500,8);
% L=zeros(1,20);
for n=20:-1:19
for k=1:1:13
if k==1
zo=[10,u(n),0.98,0.02,10,u(10),10,u(20)];
else
zo=[w(1,1),w(1,2),w(1,3),w(1,4),w(1,5),w(1,6),w(1,7),w(1,8)];
end
if n==19 || n<19
d=n+1;
xnm1=u(d);
end
tspan=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15];%[dmdt;dxdt;dxddt;dxbdt;dmfdt;dxfdt;dMtendt;dxtendt]
% tspan=[k:1:k+2];
options=odeset('Reltol',5.421011e5,'Abstol',5.421011e5);
[t,z]=ode15s(@(t,z) lumenslumens1(t,z,n,xnp1,xnm1,s,x1,k,u,errb,errd,L,V,R),tspan(k:k+1),zo,options);
for c=2:3
for e=1:8
ze(c-1,e)=z(c,e);
end
end
s=cat(1,s,ze);
tp=cat(1,tp,t);
for b=1:1:8%column%
w(1,b)=z(3,b);
end
end
for h=1:10 %row
if n==20
xnp1(h,1)=s(h,8);
else
%j=h*10+1;
xnp1(h,1)=s(h,2);
end
end
end

Answers (1)

Jan
Jan on 11 Apr 2021
The code is very confusing.
The relative tolerance of 5.421011e5 is extremely huge. Are you sure that this is meaningful in your case?
Avoid expressions like
'Reltol',5.421011e5
This looks strange:
for i=1:k
i; % Why?
...
i=i+1; % This does not influence the value of i?!
end
If FOR loops, the counter is controlled by the FOR already. So omit the two lines "i;" and "i=i+1".
If you want to store a variable persistently between calls of a function, use the "persistent" command.

Categories

Find more on Configure Simulation Conditions 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!