%
% Self-Tuning Control Simulation
%
% Process: dissolved oxygen
% Controller: PID by Cameron and Seborg
% Estimator: RLS by Bierman
%
% Bob Newell, February 1996
%
%----------------------------------------------------------
%
% Process Parameters
%
tau = 360;                % residence time (mins)
yh = 0.6;                 % yield coefficient
muhat = 4.9 / 1440;       % max growth rate (1/mins)
ks = 30;                  % substrate satn constant
koh = 0.2;                % oxygen satn constant
klac = 2.655 / 1440;      % KLa square constant
sosat = 10;               % satd oxygen concentration
cvmax = 1.516;            % 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 = [ 1 4 7 ];        % substrate concentration (mg/l)
sobar = 3;                % oxygen concentration (mg/l)
xhbar = 1449;             % biomass concentration (mg/l)
% vabar = [ 2.89 5.36 6.77 ]; % air flowrate (m3/min)
vabar = [ 2.89 ]; % air flowrate (m3/min)
pabar = 100;              % air pressure (kpa)
%
%----------------------------------------------------------
%
% Controller Parameters
%    u = uo + kp*e + ki*integ(e) + kd*deriv(e)
%
kp = [ 0.2 ];     % proportional gain (initial) 
ki = [ 0.025 ];   % integral gain (initial)
kd = [ 0 ];       % derivative gain (initial)
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.2;             % 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
%----------------------------------------------------------
%
% determine simulations to run
%
n = length( vabar );      % number of processes
l = 800;                  % length of simulation
%
yyy = [];
uuu = [];
ppp = [];
ccc = [];
%
for i = 1:n,
%
% initialise simulation   
   y = yo;             % initial output
   yd = yo;            % initial setpoint
   u = uo(i);          % initial input
   so = sobar;         % initial do
   xh = xhbar + xhstep;  % initial biomass change
   ss = ssbar(i) + ssstep;  % initial substrate change
   pa = pabar + pastep;  % initial pressure change
%
   esum = 0;           % initial integral error sum
   laste = 0;          % initial last error
%
% estimator initialisation
%
   pd = 0;
   nu = 1;
   noise = 0.5;
   update = 300;
%
   yf = 0;
   phi = [ 0 0 0 ];
   p = [ 0; 0; 0 ];
   dmat = [ 1 1 1 ];
   umat = [ 0 0 0 ];
%
   yy = [];
   uu = [];
   pp = [];
   cc = [];
%
   for k = 1:l,
      if k > 400,
	     yd = yo + ydstep;   % setpoint
      end;
%
%     controller
      e = yd - y + randn * noise;
      esum = esum + e * dt;
      dedt = ( e - laste ) / dt;
      laste = e;
      u = uo(i) + kp * e + ki * esum + kd * dedt ;
      cc = [ cc; [ kp 10*ki kd ] ];
      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 = - so / tau + ...
           ( ( yh - 1 ) / yh ) * muhat * xh * ...
    	   ( ss / ( ks + ss ) ) * ( so / ( koh + so ) ) + ...
      	   ( sosat - so ) * klac * va * va;
      so = so + dsodt * dt;
      dydt = - y / dotau + so * dogain / dotau;
      y = y + dydt * dt;
%
      yy = [ yy; y ];
%
% parameter estimation
%
      if k > 1,
%     phi update
	     phi = [ yf; phi(1); ( uu(k) - uu(k-1) ) / nu ];
%     error filter
         yf = e - pd * yf;
%     rls estimator
         [ p, dmat, umat ] = rls( yf, phi, ...
		                   p, dmat, umat, 3, 0.95 );
%     controller update
         if k > update,
		    alpha = 1 + pd;
			kp = nu * ( ( 2 * p(2) + p(1) ) / p(3) ) / alpha;
			ki = nu * ( ( 1 - p(1) - p(2) ) / p(3) ) / alpha;
			kd = nu * ( p(2) / p(3) ) / alpha;
		 end;
%
	  end;
      pp = [ pp; p' ];
%
   end;
%
   yyy = [ yyy yy ];
   uuu = [ uuu uu ];
   ppp = [ ppp pp ];
   ccc = [ ccc cc ];
%
   disp( [ 'Final parameters = ' ...
               num2str( mean( ppp(k-50:k,1) ) ) '  ' ...
               num2str( mean( ppp(k-50:k,2) ) ) '  ' ...
               num2str( mean( ppp(k-50:k,3) ) ) ] );
   disp( [ 'Final kp = ' num2str( kp ) ] );
   disp( [ 'Final ki = ' num2str( ki ) ] );
%
end;
%
ls = 'legend( ';
for j = 1:n,
   if j == 1,
      ls = [ ls sprintf( '''%5.2f''', vabar(j) ) ];
   else,
      ls = [ ls sprintf( ', ''%5.2f''', vabar(j) ) ];
   end;
end;
ls = [ ls ')' ];
%
figure
plot( dt*(1:k)', uuu )
xlabel( 'Time' )
ylabel( 'Inputs' )
%
figure
plot( dt*(1:k)', yyy )
xlabel( 'Time' )
ylabel( 'Outputs' )
%eval( ls )
%
figure
plot( dt*(1:k)', ppp )
xlabel( 'Time' )
ylabel( 'Parameters' )
%eval( ls )
%
figure
plot( dt*(1:k)', ccc )
xlabel( 'Time' )
ylabel( 'Constants' )
%eval( ls )
%
%----------------------------------------------------------
% the end
