%
% Feedback Control Simulation (plots in Figure 17-10)
%
% Process: dissolved oxygen
% Controller: PID
%
% Bob Newell, February 1996
%
%----------------------------------------------------------
%
fig = 10;
%
% Process Parameters
%
tau = 360;                % residence time (mins)
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
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
stau = 5;                 % sensor time constant (min)
sgain = 0.5 / 3;          % sensor gain
%
%ssbar = [ 1 4 7 ];
ssbar = [ 4.4 ];          % substrate concentration (mg/l)
sobar = 3;                % oxygen concentration (mg/l)
xhbar = 851.5;            % biomass concentration (mg/l)
%vabar = [ 31.2 31.5 31.8 ];  % open loop step tests
vabar = [ 31.5 ];         % air flowrate (m3/min)
pabar = 100;              % air pressure (kpa)
%
%----------------------------------------------------------
%
% Controller Parameters
%     u = uo + kp*e + ki*integral(e) + kd*derivative(e)
%
% open loop - no control
%kp = [ 0 ];
%ki = [ 0 ];
%kd = [ 0 ]; 
% Ziegler-Nichols open loop test ( L=2 T=32 R=0.046 )
%kp = [ 13.0 ];
%ki = [ 3.25 ];
%kd = [ 13.0 ];
% open loop test imc values ( K=1.47 tau=25 )
%        filter constant = 1, 5, 12
%kp = [ 8.50  3.40  1.42  ];
%ki = [ 0.340 0.136 0.057 ];
%kd = [ 0 0 0 ];
% book figure 17.10
% open loop test imc values ( K=1.47 tau1=5, tau2=20 )
%        filter constant = 1, 5, 12
if fig == 10,
  kp = [ 17.0  3.40  1.42  ];
  ki = [ 0.680 0.136 0.057 ];
  kd = [ 68.0  13.6  5.67  ];
end;
% book figure 17.11 & 17.12
if fig == 11 | fig == 12,
  kp = [ 13.0  8.50  3.40  ];
  ki = [ 3.25  0.340 0.136 ];
  kd = [ 13.0  34.0  0     ];
end;
%
dt = 0.5;                 % integrating interval
conint = 1;               % control interval as multiple of dt
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
%
if fig == 10 | fig == 11,
  ydstep = 0.05;
else
  ydstep = 0;               % size of setpoint change
end;
xhstep = 0;               % size of step biomass change
if fig == 10 | fig == 11,
  ssstep = 0;               % size of step substrate change
else
  ssstep = -0.4;
end;
pastep = 0;               % size of step pressure change
%----------------------------------------------------------
%
% determine simulations to run
%
n = length( kp );         % number of controllers
m = length( vabar );      % number of processes
len = 120;                  % 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
%
      esum = 0;           % initial integral error sum
      laste = 0;          % initial last error
%
      yy = [];
      uu = [];
%
      for k = 1:len,
%
%        controller
%
         e = yd - y;
         esum = esum + e * dt;
         dedt = ( e - laste ) / dt;
         laste = e;
         if floor(k/conint)*conint == k,
           u = uo(j) + kp(i) * e + ki(i) * esum + kd(i) * dedt;
         end;
         uu = [ uu; u ];
%
%        process
%
%        a) valve characteristics
         if linear,
            va = cvmax * u * sqrt( pa - 50 );
         else,
            va = cvmax * ( 2.236 * 20 ^ ( u - 1 ) ) ...
                                      * sqrt( pa - 50 );
         end;
%        b) biomass respiration rate
         our = ( ( yh - 1 ) / yh ) * muhat * xh * ...
    	   ( ss / ( ks + ss ) ) * ( so / ( koh + so ) );
%        c) oxygen mass balance (and simple euler integrator)
         dsodt = 1 / tau - 2 * so / tau + our + ...
      	   ( sosat - so ) * klaa * ( 1 - exp ( - va / klab ) );
         so = so + dsodt * dt;
%        c) sensor dynamics (and simple euler integrator)
         dydt = - y / stau + so * sgain / stau;
         y = y + dydt * dt;
%
         yy = [ yy; y ];
%
      end;
%
      yyy = [ yyy yy ];
	  uuu = [ uuu uu ];
%
   end;
end;
%
h = figure;
plot( dt*(1:k)', yyy, '-' );
set( gca, 'fontsize', 16 );
xlabel( 'Time in minutes' );
ylabel( 'Sensor signals' );
% book figures use
if fig == 10,
  axis( [ 0 60 0.50 0.58 ] );
  arrow( [ 7 0.566 ], [ 10 0.5515 ] );
  stext( 3, 0.57, '\16 \tau_m = 1' );
  arrow( [ 15 0.519 ], [ 8 0.532 ] );
  stext( 12, 0.515, '\16 \tau_m = 5' );
  arrow( [ 32 0.53 ], [ 19 0.538 ] );
  stext( 33, 0.53, '\16 \tau_m = 12' );
%  printsto -f1 -deps2 fig1710a.eps
end;
if fig == 11,
  axis( [ 0 60 0.50 0.58 ] );
  arrow( [ 23 0.53 ], [ 16 0.551 ] );
  stext( 20, 0.525, '\16 IMC PID \tau_m = 2' );
  arrow( [ 29 0.565 ], [ 16 0.558 ] );
  stext( 30, 0.565, '\16 IMC PI \tau_m = 5' );
  arrow( [ 19 0.575 ], [ 11 0.565 ] );
  stext( 20, 0.575, '\16 Z-N open loop test' );
%  printsto -f1 -deps2 fig1710b.eps
end;
if fig == 12,
  axis( [ 0 60 0.46 0.54 ] );
  arrow( [ 10 0.526 ], [ 16 0.505 ] );
  stext( 2, 0.53, '\16 IMC PID \tau_m = 2' );
  arrow( [ 30 0.517 ], [ 20 0.509 ] );
  stext( 30, 0.517, '\16 IMC PI \tau_m = 5' );
  arrow( [ 25 0.483 ], [ 16 0.498 ] );
  stext( 20, 0.48, '\16 Ziegler-Nichols' );
%  printsto -f1 -deps2 fig1711.eps
end;
%
%----------------------------------------------------------
% the end
