%
% Example 18.2: Feedforward Control Simulation
%
% Process: dissolved oxygen
% Controller: lead-lag from influent + PID
%
% Bob Newell, March 1996
%
%----------------------------------------------------------
%
% Process Parameters
%
vfbar = 50 / 1440;        % feed flowrate (Ml/min)
ssfbar = 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
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)
pabar = 100;              % air pressure (kpa)
%
%----------------------------------------------------------
%
% Controller Parameters
%     u = uo + kp*e + ki*integral(e) + kd*derivative(e)
%
% FB, FB + flow FF, FB + flow&COD FF
% book figure 18.6
kp = [ 0 3.40  3.40  3.40 ];
ki = [ 0 0.136 0.136 0.136 ];
kd = [ 0 13.6  13.6  13.6 ];
ksff = [ 0 0 11.14 18.85 ];
au = [ 0 0 0 0 ];
ad = [ 0 0 -0.888 -0.888 ];
bdobu = [ 0 0 0.112 0.112 ];
%
dt = 1;                 % sampling interval
uo = vabar ./ cvmax ./ ...
        sqrt( pabar - 50 );  % initial controller output
yo = 0.5;                 % initial process output
%
%----------------------------------------------------------
%
% determine simulations to run
%
n = length( kp );              % number of controllers
l = 1440 * 4;                  % length of simulation
%
%----------------------------------------------------------
%
% Define the disturbance
%
load D18_2;
t = dt * ( 1:l )';
% dd(:,1) is modelled flowrate
dd = vfbar * ( 1 - 0.35 * ( sin( pi * t / 720 + pi ) - ...
                 0.9 * sin( pi * t / 360 + 0.15 * pi ) ) );
% dd(:,2) is modelled cod
dd = [ dd ssfbar * ( 1 - 0.246 * ( sin( pi * t / 720 + pi ) - ...
                 0.9 * sin( pi * t / 360 + 0.15 * pi ) ) ) ];
% dd(:,3) is actual flowrate
dp = [];
for i = 1:47,
  dp = [ dp; d(i,1)+(d(i+1,1)-d(i,1))/120*(0:119)' ];
end;
dd = [ dd [ dp; d(48,1)*ones(120,1) ]/1440 ];
% dd(:,4) is actual cod
dp = [];
for i = 1:47,
  dp = [ dp; d(i,2)+(d(i+1,2)-d(i,2))/120*(0:119)' ];
end;
dd = [ dd [ dp; d(48,2)*ones(120,1) ] ];
%
%----------------------------------------------------------
%
%figure
%plot( t, dd(:,1)/vfbar, t, dd(:,2)/ssfbar, ...
%         t, dd(:,3)/vfbar, t, dd(:,4)/ssfbar )
%pause
%
yyy = [];
uuu = [];
%
for i = 1:n,
%
% initialise simulation   
  y = yo;             % initial output
  yd = yo;            % setpoint
  u = uo(j);          % initial input
  so = sobar;         % initial do
  xh = xhbar;         % initial biomass change
  ss = ssbar;         % initial substrate change
  pa = pabar;         % initial pressure change
  vf = vfbar;         % initial influent flowrate
%
  esum = 0;           % initial integral error sum
  laste = 0;          % initial last error
  lastu = 0;
  lastd = 0;
%
  yy = [];
  uu = [];
%
  for k = 1:l,
%
%        disturbance
    vf = dd( k, 3 );
    ssf = dd( k, 4 );
%
%        feedback controller
    e = yd - y;
    esum = esum + e * dt;
    dedt = ( e - laste ) / dt;
    laste = e;
    u = uo + kp(i) * e + ki(i) * esum + kd(i) * dedt;
%
%        feedforward compensator
    dist = ( dd( k, 3 ) - vfbar ) ;
    uff = - ad(i) * lastu + bdobu(i) * ...
                   ( dist + au(i) * lastd );
    lastu = uff;
    lastd = dist;
    u = u + ksff(i) * uff;
    uu = [ uu; u ];
%
%        process
    va = cvmax * u * sqrt( pa - 50 );
    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;
%
h = figure;
plot( dt*(1:2880)'/60, 6*yyy(1:2880,[2 3 4]), '-' );
set( gca, 'fontsize', 16 );
xlabel( 'Time in hours' );
ylabel( 'DO in mg/l' );
%
arrow( [ 36 3.35 ], [ 28.5 3.2 ] );
h = text( 37, 3.35, 'Type 1' );
set( h, 'fontsize', 16 );
arrow( [ 36 3.25 ], [ 28.6 3.1 ] );
h = text( 37, 3.25, 'Type 2' );
set( h, 'fontsize', 16 );
arrow( [ 36 3.15 ], [ 29.5 3.02 ] );
h = text( 37, 3.15, 'Type 3' );
set( h, 'fontsize', 16 );
%
print -f1 -deps2 fig18-5.eps
%
%----------------------------------------------------------
% the end
