%
% Feedback Control Simulation
%
% Process: dissolved oxygen
% Controller: on-off
%
% Bob Newell, October 1997
%
%----------------------------------------------------------
%
% 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
dotau = 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
%
somax = 3.5;
somin = 2.5;
%
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 = dogain * somin ;                 % initial process output
%
%----------------------------------------------------------
%
%----------------------------------------------------------
%
% determine simulations to run
%
n = 1;         % number of controllers
m = 1;         % number of processes
l = 120;                  % length of simulation
%
yyy = [];
uuu = [];
%
for i = 1:n,
   for j = 1:m,
%
% initialise simulation   
      y = yo;             % initial output
      u = 1;              % initial input
      so = somin;         % initial do
      xh = xhbar;         % initial biomass change
      ss = ssbar;         % initial substrate change
      pa = pabar;         % initial pressure change
%
      yy = [];
      uu = [];
%
      for k = 1:l,
%
%        controller
         if y > dogain * somax,
           u = 0;
         end;
         if y < dogain * somin,
           u = 1;
         end; 
         uu = [ uu; u ];
%
%        process
         if linear,
            va = cvmax * u * sqrt( pa - 50 );
         else,
            va = cvmax * ( 2.236 * 20 ^ ( u - 1 ) ) ...
	                        * sqrt( pa - 50 );
	 end;
         dsodt = 1 / tau - 2 * so / tau + ...
           ( ( yh - 1 ) / yh ) * muhat * xh * ...
    	   ( ss / ( ks + ss ) ) * ( so / ( koh + so ) ) + ...
      	   ( sosat - so ) * klaa * ( 1 - exp ( - va / klab ) );
         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;
%
h = figure;
hax = plotyy( dt*(1:k)', yyy, dt*(1:k)', uuu );
axesyy( hax, [ 0 60 0.15 0.65 -1 5 ], ...
   'Time in minutes', ...
   'DO sensor signal', ...
   'Valve signal', 2, ...
   [ 0.4 0.5 0.6 ], [ 0 1 ] );
%
h = line( [ 0 60 ], [ dogain*somax dogain*somax ] );
set( h, 'color', 'k' );
h = line( [ 0 60 ], [ dogain*somin dogain*somin ] );
set( h, 'color', 'k' );
%
arrow( [ 16 dogain*somin-0.05 ], [ 16 dogain*somin ] );
arrow( [ 16 dogain*somax+0.05 ], [ 16 dogain*somax ] );
h = text( 1, 0.37, 'deadband' );
set( h, 'fontsize', 16 );
%
printsto -f1 -deps2 fig17-6.eps
%
%----------------------------------------------------------
% the end
