%
% Example 18.1: Feedforward Control Simulation
%
% Process: dissolved oxygen
% Controller: lead-lag + PID
%
% Bob Newell, March 1996
%
%----------------------------------------------------------
%
fig = 1;
%
% Process Parameters
%
vfbar = 50 / 1440;        % feed flowrate (Ml/min)
ssf = 151.7;              % feed concentration (mg/l)
vol = 12.5;               % volume (Ml) 
yh = 0.6;                 % yield coefficient
muhat = 3.733 / 1440;     % max growth rate (1/mins)
ks = 20;                  % substrate satn constant
koh = 0.2;                % oxygen satn constant
fp = 0.1;                 % fraction inerts
bh = 0.4 / 1440;          % biomass decay rate
klaa = 166 / 1440;        % KLa constants
klab = 80;
sosat = 10;               % satd oxygen concentration
cvmax = 8.909;            % air valve size coefficient
linear = 1;               % linear or equal percentage
dotau = 1.5;              % do sensor time constant (min)
dogain = 0.5 / 3;         % do sensor gain
%
ssbar = [ 4.4 ];          % substrate concentration (mg/l)
sobar = 3;                % oxygen concentration (mg/l)
xhbar = 851.5;            % biomass concentration (mg/l)
vabar = [ 31.5 ];         % air flowrate (m3/min)
%vabar = [ 29 31.5 34 ];
pabar = 100;              % air pressure (kpa)
%
%----------------------------------------------------------
%
% Controller Parameters
%     u = uo + kp*e + ki*integral(e) + kd*derivative(e)
%
% open loop - no controllers at all
if fig == 0,
  kp = [ 0 ];     % proportional gain (vector)
  ki = [ 0 ];     % integral gain (vector)
  kd = [ 0 ];     % derivative gain (vector)
  ksff = [ 0 ];   % static feedforward gain (vector)
  au = [ 0 ];     % dynamic feedforward u coeff (vector)
  ad = [ 0 ];     % dynamic feedforward d coeff (vector)
  bdobu = [ 0 ];  % dynamic feedforward gain (vector)
end;
% feedforward only control - no feedback
% book figure 18.4 left plot
if fig == 1,
  kp = [ 0 0 0 0 ];
  ki = [ 0 0 0 0 ];
  kd = [ 0 0 0 0 ];
  ksff = [ 0 11.7 11.7 11.14 ];
  au = [ 0 0 -0.9785 0 ];
  ad = [ 0 0 -0.9840 -0.942 ];
  bdobu = [ 0 1 0.7442 0.058 ];
end;
% feedforward and feedback control
% book figure 18.4 right plot
if fig == 2,
  kp = [ 3.40  3.40  3.40  3.40 ];
  ki = [ 0.136 0.136 0.136 0.136 ];
  kd = [ 13.6  13.6  13.6  13.6 ];
  ksff = [ 0 11.7 11.7 11.14 ];
  au = [ 0 0 -0.9785 0 ];
  ad = [ 0 0 -0.9840 -0.942 ];
  bdobu = [ 0 1 0.7442 0.058 ];
end;
%
dt = 0.5;                 % sampling interval
if linear,
   uo = vabar ./ cvmax ./ ...
        sqrt( pabar - 50 );  % initial controller output
else,
   uo = 1 + log( vabar ./ cvmax ./ 2.236 ./ ...
        sqrt( pabar - 50 ) ) ./ ...
	    log( 20 );            % initial controller output
end;
yo = 0.5;                 % initial process output
%
%----------------------------------------------------------
%
% Define the disturbance
%
ydstep = 0;               % size of setpoint change
xhstep = 0;               % size of step biomass change
ssstep = 0;               % size of step substrate change
pastep = 0;               % size of step pressure change
vfstep = -5 / 1440;       % size of feed flowrate change
%----------------------------------------------------------
%
% determine simulations to run
%
n = length( kp );         % number of controllers
m = length( vabar );      % number of processes
l = 400;                  % length of simulation
%
yyy = [];
uuu = [];
%
for i = 1:n,
  for j = 1:m,
%
% initialise simulation   
    y = yo;             % initial output
    yd = yo + ydstep;   % setpoint
    u = uo(j);          % initial input
    so = sobar;         % initial do
    xh = xhbar + xhstep;  % initial biomass change
    ss = ssbar + ssstep;  % initial substrate change
    pa = pabar + pastep;  % initial pressure change
    vf = vfbar + vfstep;  % initial influent flowrate
%
    esum = 0;           % initial integral error sum
    laste = 0;          % initial last error
    lastu = 0;
    lastd = 0;
%
    yy = [];
    uu = [];
%
    for k = 1:l,
%
%        feedback controller
      e = yd - y;
      esum = esum + e * dt;
      dedt = ( e - laste ) / dt;
      laste = e;
      u = uo(j) + kp(i) * e + ki(i) * esum + kd(i) * dedt;
%
%        feedforward compensator
      dist = ( vf - vfbar ) ;
      uff = - ad(i) * lastu + bdobu(i) * ...
                        ( dist + au(i) * lastd );
      lastu = uff;
      lastd = dist;
      u = u + ksff(i) * uff;
      uu = [ uu; u ];
%
%        process
      if linear,
        va = cvmax * u * sqrt( pa - 50 );
      else,
        va = cvmax * ( 2.236 * 20 ^ ( u - 1 ) ) ...
                       * sqrt( pa - 50 );
      end;
      dssdt = ( ssf - ss ) * vf / vol - ...
           ( 1 / yh ) * muhat * xh * ...
    	   ( ss / ( ks + ss ) ) * ( so / ( koh + so ) ) + ...
      	   ( 1 - fp ) * bh * xh;
      dsodt = vf / vol - 2 * so * vf / vol + ...
           ( ( yh - 1 ) / yh ) * muhat * xh * ...
    	   ( ss / ( ks + ss ) ) * ( so / ( koh + so ) ) + ...
      	   ( sosat - so ) * klaa * ( 1 - exp ( - va / klab ) );
      ss = ss + dssdt * dt;
      so = so + dsodt * dt;
      dydt = - y / dotau + so * dogain / dotau;
      y = y + dydt * dt;
%
      yy = [ yy; [ y ] ];
%
    end;
%
    yyy = [ yyy yy ];
    uuu = [ uuu uu ];
%
  end;
end;
%
figure;
plot( dt*(1:k)', yyy, '-' );
set( gca, 'fontsize', 24 );
xlabel( 'Time in minutes' );
ylabel( 'Sensor signals' );
% for the book figures
if fig == 1,
  h = text( 70, 0.56, 'open loop' );
  set( h, 'fontsize', 24 );
  arrow( [ 40 0.53 ], [ 12 0.494 ] );
  h = text( 35, 0.535, 'Dynamic FF (STD)' );
  set( h, 'fontsize', 24 );
  arrow( [ 78 0.518 ], [ 32 0.503 ] );
  h = text( 80, 0.518, 'Dynamic FF (MBD)' );
  set( h, 'fontsize', 24 );
  h = text( 50, 0.492, 'static FF' );
  set( h, 'fontsize', 24 );
  print -f1 -deps2 fig18-4a.eps
end;
if fig == 2,
  axis( [ 0 200 0.49 0.51 ] );
  arrow( [ 43 0.5085 ], [ 25 0.5067 ] );
  h = text( 45, 0.5085, 'FB only' );
  set( h, 'fontsize', 24 );
  arrow( [ 48 0.506 ], [ 12 0.5005 ] );
  h = text( 50, 0.506, 'FB + Dynamic FF (MBD)' );
  set( h, 'fontsize', 24 );
  arrow( [ 48 0.496 ], [ 15 0.499 ] );
  h = text( 50, 0.496, 'FB + Dynamic FF (STD)' );
  set( h, 'fontsize', 24 );
  arrow( [ 33 0.493 ], [ 8 0.496 ] );
  h = text( 35, 0.493, 'FB + static FF' );
  set( h, 'fontsize', 24 );
  print -f1 -deps2 fig18-4b.eps
end;
%
%----------------------------------------------------------
% the end
