Download Modeling and Simulation in Scilab/Scicos

148
5 Examples
5.4.1 Mathematical Model
Many applications use polynomial models. In some cases these are argued for on first
principles. In other applications they occur because if linear equations do not suffice, then
the next natural thing to try is a polynomial equation. If a linear approximation is not
adequate, then the next level of Taylor approximation uses quadratic functions. The DAE
we consider arises as the model of a batch reactor and is sometimes used as a test problem
in the chemical engineering literature [17, 31]. The equations are
ẏ1 = −p3 y2 y8 ,
ẏ2 = −p1 y2 y6 + p2 y10 − p3 y2 y8 ,
(5.59a)
(5.59b)
ẏ3 = p3 y2 y8 + p4 y4 y6 − p5 y9 ,
ẏ4 = −p4 y4 y6 + p5 y9 ,
(5.59c)
(5.59d)
ẏ5 = p1 y2 y6 − p2 y10 ,
ẏ6 = −p1 y2 y6 − p4 y4 y6 + p2 y10 + p5 y9
(5.59e)
(5.59f)
0 = −0.0131 + y6 + y8 + y9 + y10 − y7 ,
0 = p7 y1 − y8 (p7 + y7 ),
(5.59g)
(5.59h)
0 = p8 y3 − y9 (p8 + y7 ),
0 = p6 y5 − y10 (p6 + y7 ).
(5.59i)
(5.59j)
In its original form this problem is known to be very stiff and is used for a test problem
for parameter sensitivity software. We shall take a slightly modified set of initial conditions.
5.4.2 Scilab Implementation
To illustrate the procedure we shall proceed as follows. Our observational data will be
given by a simulation followed by the addition of random perturbations. We begin with a
fixed value of p, which we denote by ptrue. Assuming that p is given in the current Scilab
environment, the system is coded as follows:
function [r,ires]=reactor(t,y,ydot)
ires=0;
r(1)=ydot(1)+p(3)*y(2)*y(8);
r(2)=ydot(2)+p(1)*y(2)*y(6)-p(2)*y(10)+p(3)*y(2)*y(8);
r(3)=ydot(3)-p(3)*y(2)*y(8)-p(4)*y(4)*y(6)+p(5)*y(9);
r(4)=ydot(4)+p(4)*y(4)*y(6)-p(5)*y(9);
r(5)=ydot(5)-p(1)*y(2)*y(6)+p(2)*y(10);
r(6)=ydot(6)+p(1)*y(2)*y(6)+p(4)*y(4)*y(6)-p(2)*y(10)-p(5)*y(9);
r(7)=-0.0131+y(6)+y(8)+y(9)+y(10)-y(7);
r(8)=p(7)*y(1)-y(8)*(p(7)+y(7));
r(9)=p(8)*y(3)-y(9)*(p(8)+y(7));
r(10)=p(6)*y(5)-y(10)*(p(6)+y(7));
endfunction
Since the model (5.59) is a DAE, only some initial conditions will be consistent. We
shall assume that {p6 , p7 , p8 } are known and fixed at the values ptrue(6:8). In this case
the space of consistent initial conditions is independent of the remaining values of pi and
we do not have to worry about finding new consistent initial conditions as the parameters
are varied.