%
% SIMPLE CARBON REMOVAL PLANT MODEL
%
% Described in Appendix B
%
%
% Sensitivity Calculations
%
% Bob Newell, January 1996
%
global vt
%
vt(1) = 50 ;           % feed flowrate (Ml/d)
vt(2) = 200 ;          % feed substrate (mg/l)
vt(3) = 25 ;           % recycle flowrate (Ml/d)
vt(4) = 0.6 ;          % wastage flowrate (Ml/d)
vt(5) = 5.36 ;         % air flowrate (m3/min)
vt(6) = 4 ;            % tank/effluent substrate (mg/l)
vt(7) = 3 ;            % tank/effluent oxygen (mg/l)
vt(8) = 1449 ;         % tank heterotrophs (mg/l)
vt(9) = 12.5 ;         % tank volume (Ml)
vt(10) = 10 ;          % saturated oxygen (mg/l)
vt(11) = 4.9 ;         % maximum specific growth rate (1/d)
vt(12) = 30 ;          % substrate saturation constant (mg/l)
vt(13) = 0.2 ;         % oxygen saturation constant (mg/l)
vt(14) = 0.4 ;         % biomass decay rate (1/d)
vt(15) = 0.6 ;         % yield coefficient (-)
vt(16) = 0.1 ;         % fraction inerts on decay (-)
vt(17) = 12 ;          % KLa-airflow regression slope (min/m3/d)
vt(18) = 12 ;          % KLa-airflow regression intercept (1/d)
%
% dynamic sensitivity calculations
%
sm = [] ;
h1 = '                  Plus 1 percent        ' ;
h2 = '        Minus 1 percent      ' ;
disp( [ h1 h2 ] ) ;
h1 = '     No        Ss        So        Xh   ' ;
h2 = '     Ss        So        Xh   ' ;
disp( [ h1 h2 ] ) ,
for i = 1:18 ,
   smp = [ i ] ;
   vtsave = vt(i) ;
   vt(i) = 0.99 * vtsave ;
   xo = vt(6:8) ;
   smp = [ smp State_Derivatives( 0, xo ) ./ xo ] ;
   vt(i) = 1.01 * vtsave ;
   xo = vt(6:8) ;
   smp = [ smp State_Derivatives( 0, xo ) ./ xo ] ;
   vt(i) = vtsave ;
   sm = [ sm; smp ] ;
end ;
disp( sm )
%
disp( '            Sensitivity      Nonlinearity' ) ;
disp( '    No    Ss    So    Xh    Ss    So    Xh' ) ;
sm = abs( sm ) ;
smmax = max( sm ) ;
mag = zeros( size( sm ) ) ;
for i = 1:18 ,
   mag(i,1) = floor( sm(i,1) ) ;
   for j = 2:4 ,
      mag(i,j) = floor( 3 * sm(i,j) / smmax(j) ) ;
   end ;  
   for j = 5:7 ,
      if mag(i,j-3) > 0,
         mag(i,j) = floor( 15 * abs( sm(i,j) - sm(i,j-3) ) ...
	         / sm(i,j-3) ) ;
      end ;
   end ;  
end ;
disp( mag )
%
%
% steady state sensitivity calculations
%
%
% initial state vector
%
xo = vt(6:8) ;
%
x = zeros( 1, 3 );
global x
%
h1 = '                  Plus 1 percent        ' ;
h2 = '        Minus 1 percent      ' ;
disp( [ h1 h2 ] ) ;
h1 = '     No        Ss        So        Xh   ' ;
h2 = '     Ss        So        Xh   ' ;
disp( [ h1 h2 ] ) ,
sm = [] ;
for i = 1:8 ,
   vtsave = vt(i+10) ;
   vt(i+10) = 1.01 * vtsave ;
	c1 = vt(14) + vt(4) * ( vt(1) + vt(3) ) ...
      / vt(9) / ( vt(4) + vt(3) ) ;
   tau = vt(9) / vt(1) ;
   xhmax1 = ( vt(2) - 0.1 ) / tau / ( c1 / vt(15) ...
      - ( 1 - vt(16) ) * vt(14) ) ;   
   kla = vt(17) * vt(5) + vt(18) ;
   xhmax2 = ( 0.01 / tau - kla * ( vt(10) - 0.01 ) ) / ...
      ( ( vt(15) - 1 ) * c1 / vt(15) ) ;
   xhmax = min( [ xhmax1 xhmax2 ] ) ;
   fmin( 'ap2func3', 500, xhmax ) ;
   smp = [ i+10 ( x - xo ) ./ xo ] ;
   vt(i+10) = 0.99 * vtsave ;
	c1 = vt(14) + vt(4) * ( vt(1) + vt(3) ) ...
      / vt(9) / ( vt(4) + vt(3) ) ;
   tau = vt(9) / vt(1) ;
   xhmax1 = ( vt(2) - 0.1 ) / tau / ( c1 / vt(15) ...
      - ( 1 - vt(16) ) * vt(14) ) ;   
   kla = vt(17) * vt(5) + vt(18) ;
   xhmax2 = ( 0.01 / tau - kla * ( vt(10) - 0.01 ) ) / ...
      ( ( vt(15) - 1 ) * c1 / vt(15) ) ;
   xhmax = min( [ xhmax1 xhmax2 ] ) ;
   fmin( 'ap2func3', 500, xhmax ) ;
   smp = [ smp ( x - xo ) ./ xo ] ;
   vt(i+10) = vtsave ;
   sm = [ sm; smp ] ;
end ;
disp( sm ) ;
%
disp( '            Sensitivity      Nonlinearity' ) ;
disp( '    No    Ss    So    Xh    Ss    So    Xh' ) ;
sm = abs( sm ) ;
smmax = max( sm ) ;
mag = zeros( size( sm ) ) ;
for i = 1:8 ,
   mag(i,1) = floor( sm(i,1) ) ;
   for j = 2:4 ,
      mag(i,j) = floor( 3 * sm(i,j) / smmax(j) ) ;
   end ;  
   for j = 5:7 ,
      if mag(i,j-3) > 0,
         mag(i,j) = floor( 15 * abs( sm(i,j) - sm(i,j-3) ) ...
	         / sm(i,j-3) ) ;
	   end ;
   end ;  
end ;
disp( mag )
%
%
