%
% Feedback Control Simulation
%
% Process: dissolved oxygen
% Controller: PID
%
% Bob Newell, February 1996
%
%----------------------------------------------------------
%
fig = 15;
%
% 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)
ssbar = [ 1.1 4.4 13.2 ];
sobar = 3;                % oxygen concentration (mg/l)
xhbar = 851.5;            % biomass concentration (mg/l)
vabar = [ 8.99 31.5 95.5 ];  % 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 13.10
% open loop test imc values ( K=1.47 tau1=5, tau2=20 )
%        filter constant = 1, 5, 12
%kp = [ 17.0  3.40  1.42  ];
%ki = [ 0.680 0.136 0.057 ];
%kd = [ 68.0  13.6  5.67  ];
% book figure 13.11 & 13.12
%kp = [ 13.0  8.50  3.40  ];
%ki = [ 3.25  0.340 0.136 ];
%kd = [ 13.0  34.0  0     ];
% book figure 17.13
if fig == 13,
  kp = [ 3.40  3.40  3.40  ];
  ki = [ 0.136 0.136 0.136  ];
  kd = [ 0     0     0     ];
end;
% book figure 17.15
if fig == 15,
  kp = [ 3.40/2.0127  3.40  3.40/0.2845  ];
  ki = [ 0.136/2.0127 0.136 0.136/0.2845  ];
  kd = [ 0     0     0     ];
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
ydstep = 0.05;
xhstep = 0;               % size of step biomass change
ssstep = 0;               % size of step substrate change
%ssstep = -0.4;
pastep = 0;               % size of step pressure change
%----------------------------------------------------------
%
% determine simulations to run
%
n = length( kp );         % number of controllers
m = length( ssbar );      % number of processes
l = 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(j) + 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:l,
%
%        controller
         e = yd - y;
         esum = esum + e * dt;
         dedt = ( e - laste ) / dt;
         laste = e;
         u = uo(j) + kp(j) * e + ki(j) * esum + kd(j) * dedt ;
         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;
plot( dt*(1:k)', yyy, '-' );
set( gca, 'fontsize', 16 );
xlabel( 'Time in minutes' );
ylabel( 'Sensor signals' );
% book figures use
axis( [ 0 60 0.50 0.58 ] );
if fig == 13,
  arrow( [ 30 0.524 ], [ 25 0.536 ] );
  h = text( 28, 0.52, 'high' );
  set( h, 'fontsize', 16 );
  arrow( [ 8 0.572 ], [ 13 0.557 ] );
  h = text( 3, 0.575, 'medium' );
  set( h, 'fontsize', 16 );
  arrow( [ 29 0.57 ], [ 18 0.563 ] );
  h = text( 30, 0.57, 'low' );
  set( h, 'fontsize', 16 );
  print -f1 -deps2 fig17-13.eps
end;
if fig == 15,
  arrow( [ 30 0.523 ], [ 25 0.545 ] );
  h = text( 28, 0.52, 'high' );
  set( h, 'fontsize', 16 );
  arrow( [ 8 0.572 ], [ 14 0.557 ] );
  h = text( 3, 0.575, 'medium' );
  set( h, 'fontsize', 16 );
  arrow( [ 39 0.57 ], [ 28 0.559 ] );
  h = text( 40, 0.57, 'low' );
  set( h, 'fontsize', 16 );
  print -f1 -deps2 fig17-15.eps
end;
%
%----------------------------------------------------------
% the end





