%
% Feedback Control Simulation
%
% Process: nitrification
% Controller: PID
%
% Bob Newell, November 1997
%
%----------------------------------------------------------
%
% Process Parameters
%
tau = 360;                % residence time (mins)
ya = 0.24;                % yield coefficient
muhat = 0.8 / 1440;       % max growth rate (1/mins)
knh4 = 1.0;               % ammonia satn constant
koa = 0.4;                % oxygen satn constant
dotau = 5;                % do control time constant (min)
dogain = 2 / 0.5;         % do control gain
domax = 6;                % maximum do capacity
stau = 15;                % sensor time constant (min)
sgain = 0.5 / 1;          % sensor gain
%
snh4inbar = [ 20 ];       % influent substrate concentration (mg/l)
snh4bar = 1.0;            % effluent substrate concentration (mg/l)
sobar = 2;                % oxygen concentration (mg/l)
xabar = 54.72;            % biomass concentration (mg/l)
%
%----------------------------------------------------------
%
% Controller Parameters
%     u = uo + kp*e + ki*integral(e) + kd*derivative(e)
%
% open loop - no control
kp = [ 0 ];
ki = [ 0 ];
kd = [ 0 ]; 
% closed loop - P control test
kp = [ -15 ];
ki = [ 0 ];
kd = [ 0 ]; 
%
dt = 0.5;                 % integrating interval
conint = 1;               % control interval as multiple of dt
yo = 0.5;                 % initial process output
uo = 0.5;                 % initial process input
%
%----------------------------------------------------------
%
% Define the disturbance
%
ydstep = 0.05;            % size of setpoint change
xastep = 0;               % size of step biomass change
snh4step = 0;             % size of step substrate change
%----------------------------------------------------------
%
% determine simulations to run
%
n = length( kp );         % number of controllers
len = 720;                % length of simulation
%
yyy = [];
uuu = [];
%
for i = 1:n,
%
% initialise simulation   
  y = yo;             % initial output
  yd = yo + ydstep;   % setpoint
  u = uo;             % initial input
  so = dogain * uo;   % initial do
  snh4 = snh4bar;     % initial ammonia
  xa = xabar + xastep;  % initial biomass change
  snh4in = snh4inbar + snh4step;  % initial substrate 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 + kp(i) * e + ki(i) * esum + kd(i) * dedt;
      if u < 0, u = 0; end;
      if u > domax / dogain, u = domax / dogain; end;
    end;
    uu = [ uu; u ];
%
%   process
%
%   a) DO control dynamics (and simple euler integrator)
    dsodt = - so / dotau + u * dogain / dotau;
    so = so + dsodt * dt;
%   b) autotrophic biomass respiration rate
    nur = ( ( ya - 4.57 ) / ya ) * muhat * xa * ...
    	   ( snh4 / ( knh4 + snh4 ) ) * ( so / ( koa + so ) );
%   c) ammonia mass balance (and simple euler integrator)
    dsnh4dt = snh4in / tau - snh4 / tau - nur / ( ya - 4.57 );
    snh4 = snh4 + dsnh4dt * dt;
%   d) nitrate mass balance (and simple euler integrator)
%    dsnoxdt = snoxin / tau - snox / tau + ya * nur;
%    snox = snox + dsnoxdt * dt;
%   e) sensor dynamics (and simple euler integrator)
    dydt = - y / stau + snh4 * sgain / stau;
    y = y + dydt * dt;
%
    yy = [ yy; y ];
%
  end;
%
  yyy = [ yyy yy ];
  uuu = [ uuu uu ];
%
end;
%
h = figure;
subplot( 2, 1, 1 );
plot( dt*(1:k)', yyy );
set( gca, 'fontsize', 16 );
xlabel( 'Time in minutes' );
ylabel( 'Sensor signals' );
%
subplot( 2, 1, 2 );
plot( dt*(1:k)', uuu );
set( gca, 'fontsize', 16 );
xlabel( 'Time in minutes' );
ylabel( 'Valve signals' );
%
%----------------------------------------------------------
% the end












