Reference
Neculai Andrei, Nonlinear Optimization Applications Using the GAMS Technology,Springer Optimization and Its Applications, Model Protein (12.23) in chapter Optimal Control , 2013
$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