Reference
Category : GAMS NOA library
Mainfile : protein.gms
$onText
Optimal production of secreted protein in a fed-batch reactor.
Park, S., Ramirez, W.F., Optimal production of secreted protein in fed-batch
reactors. A.I.Ch.E. Journal, 34, 1988, pp.1550-1558.
$offText
$if set n $set nh %n%
$if not set nh $set nh 500
set nh Number of subintervals / 0*%nh% /;
alias (nh,k);
Scalar tf final time / 10 /
x1_0 initial value for x1 / 1.0 /
x2_0 initial value for x2 / 5.0 /
x3_0 initial value for x3 / 0.0 /
x4_0 initial value for x4 / 0.0 /
x5_0 initial value for x5 / 1.0 /
h ;
h=tf/%nh%;
Variables x1(nh)
x2(nh)
x3(nh)
x4(nh)
x5(nh)
u(nh) control variable
a1(nh)
a2(nh)
a3(nh)
obj criterion ;
Equations eobj criterion definition
state1(nh) state equation 1
state2(nh) state equation 2
state3(nh) state equation 3
state4(nh) state equation 4
state5(nh) state equation 5
ea1
ea2
ea3;
eobj..
obj =e= x4['%nh%']*x5['%nh%'] ;
state1(nh(k+1))..
x1[k+1] =e= x1(k)+
(h/2)*( a1(k)*x1(k) - u(k)*x1(k)/x5(k) +
a1(k+1)*x1(k+1) - u(k+1)*x1(k+1)/x5(k+1) ) ;
state2(nh(k+1))..
x2[k+1] =e= x2(k)+
(h/2)*( -7.3*a1(k)*x1(k) - u(k)*(x2(k)-20)/x5(k)
-7.3*a1(k+1)*x1(k+1) - u(k+1)*(x2(k+1)-20)/x5(k+1) );
state3(nh(k+1))..
x3[k+1] =e= x3(k)+
(h/2)*( a2(k)*x1(k) - u(k)*x3(k)/x5(k) +
a2(k+1)*x1(k+1) - u(k+1)*x3(k+1)/x5(k+1) );
state4(nh(k+1))..
x4[k+1] =e= x4(k)+
(h/2)*( a3(k)*(x3(k)-x4(k)) - u(k)*x4(k)/x5(k) +
a3(k+1)*(x3(k+1)-x4(k+1)) - u(k+1)*x4(k+1)/x5(k+1) );
state5(nh(k+1))..
x5[k+1] =e= x5(k) + (h/2)*( u(k) + u(k+1) );
ea1(nh(k)).. a1(k) =e= 21.87*x2(k)/((x2(k)+0.4)*(x2(k)+62.5));
ea2(nh(k)).. a2(k) =e= (x2(k)*exp(-5*x2(k)))/(0.1+x2(k));
ea3(nh(k)).. a3(k) =e= 4.75*a1(k)/(0.12+a1(k));
*Initial point
x1.l[nh]=1.0;
x2.l[nh]=5.0;
x3.l[nh]=0.0;
x4.l(nh)=0.0;
x5.l(nh)=1.0;
u.l(nh) =0.0;
x1.fx ['0'] = x1_0;
x2.fx ['0'] = x2_0;
x3.fx ['0'] = x3_0;
x4.fx ['0'] = x4_0;
x5.fx ['0'] = x5_0;
*Bounds
u.lo(nh) = 0.0; u.up(nh) = 5;
Model protein /all/;
option reslim=60000;
protein.iterlim=80000;
Solve protein maximizing obj using nlp;
$ifThenI x%mode%==xbook
file stat1 /protein1.dat/;
file stat2 /protein2.dat/;
file stat3 /protein3.dat/;
file stat4 /protein4.dat/;
file stat5 /protein5.dat/;
file cont /protein.dat/;
put stat1;
loop(nh, put x1.l(nh):10:5,',', put/)
put stat2;
loop(nh, put x2.l(nh):10:5,',', put/)
put stat3;
loop(nh, put x3.l(nh):10:5,',', put/)
put stat4;
loop(nh, put x4.l(nh):10:5,',', put/)
put stat5;
loop(nh, put x5.l(nh):10:5,',', put/)
put cont;
loop(nh, put u.l(nh):10:5,',', put/)
$endIf
* End of protein