Reference
Neculai Andrei, Nonlinear Optimization Applications Using the GAMS Technology,Springer Optimization and Its Applications, Model Ethanol (12.11) in chapter Optimal Control , 2013
$onText
Optimal control of fed-batch reactor for ethanol production.
A fed-batch bioreactor for the production of ethanol from the anaerobic
glucose fermentation by Saccharomyces cerevisiae.

Larrosa, J.A.E., New heuristics for global optimization of complex
bioprocesses. Ph.D. Thesis, Universidade de Vigo, Departamento de Enxeera
Qumica, Vigo, 2008.
$offText

$if     set n  $set nh %n%
$if not set nh $set nh 100

set nh Number of subintervals / 0*%nh% /;
alias (nh,k);

Scalar    tf   final time /54/
          y1_0 initial value for y1 / 20/
          y2_0 initial value for y2 / 150 /
          y3_0 initial value for y3 / 21 /
          y4_0 initial value for y4 / 10 /
          a    a parameter in objective function /0.0/
          h ;

          h=tf/%nh%;

Variables y1(nh)  microbial population concentration
          y2(nh)  substrate concentration
          y3(nh)  the product concentration
          y4(nh)  volume
          u(nh)   control variable
          p1(nh)
          p2(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
          ep1(nh)
          ep2(nh) ;

* Objective function to be maximized:
eobj..
 obj =e= y3['%nh%']*y4['%nh%'] +
         h*a*sum(nh(k+1), ((u(k+1)-u(k))*(u(k+1)-u(k))) );

* Constraints:
state1(nh(k+1))..
y1[k+1] =e= y1(k)+
   (h/2)*( p1(k)*y1(k)    - u(k)*y1(k)/y4(k) +
           p1(k+1)*y1(k+1)- u(k+1)*y1(k+1)/y4(k+1) ) ;

state2(nh(k+1))..
y2[k+1] =e= y2(k)+
   (h/2)*( (-10)*p1(k)*y1(k) + u(k)*((150-y2(k))/y4(k)) +
           (-10)*p1(k+1)*y1(k+1) + u(k+1)*((150-y2(k+1))/y4(k+1)) );

state3(nh(k+1))..
y3[k+1] =e= y3(k)+
   (h/2)*( p2(k)*y1(k) - u(k)*y3(k)/y4(k) +
           p2(k+1)*y1(k+1) - u(k+1)*y3(k+1)/y4(k+1) );

state4(nh(k+1))..
y4[k+1] =e= y4(k) + (h/2)*(u(k) + u(k+1));

ep1(nh(k)).. p1(k) =e= (0.408/(1+y3(k)/16))*(y2(k)/(0.22+y2(k)));
ep2(nh(k)).. p2(k) =e= (1/(1+y3(k)/71.5))*(y2(k)/(0.44+y2(k)));

*Initial point
y1.l[nh]=20;
y2.l[nh]=150;
y3.l[nh]=40;
y4.l[nh]=10;
u.l[nh] =10;

* Fixed values:
y1.fx ['0'] = y1_0;
y2.fx ['0'] = y2_0;
y3.fx ['0'] = y3_0;
y4.fx ['0'] = y4_0;
y4.fx ['%nh%'] = 200;

*Bounds
u.lo(nh) = 0;     u.up(nh) = 12;

y1.lo(nh) = 0;
y2.lo(nh) = 0;
y3.lo(nh) = 0;
y4.lo(nh) = 0;

Model ethanol /all/;

Solve ethanol maximizing obj using nlp;

$ifThenI x%mode%==xbook
file stat1 /eth1.dat/;
file stat2 /eth2.dat/;
file stat3 /eth3.dat/;
file stat4 /eth4.dat/;
file cont  /eth.dat/;

put stat1;
loop(nh, put y1.l(nh):10:5,',', put/)

put stat2;
loop(nh, put y2.l(nh):10:5,',', put/)

put stat3;
loop(nh, put y3.l(nh):10:5,',', put/)

put stat4;
loop(nh, put y4.l(nh):10:5,',', put/)

put cont;
loop(nh, put u.l(nh):10:5,',', put/)
$endIf

* End of ethanol