%
% Feedback Control Simulation
%
% Process: denitrification
% Controller: PID
%
% Bob Newell, November 1997
%
%----------------------------------------------------------
%
% Process Parameters
%
tau = 120;                % residence time (mins)
yh = 0.67;                % yield coefficient
muhat = 6.0 / 1440;       % max growth rate (1/mins)
knox = 0.5;               % nitrate satn constant
koa = 20;                 % substrate satn constant
cdtau = 5;                % dosage control time constant (min)
cdgain = 0.480555 / 0.5;  % dosage control gain
cdmax = 1;                % maximum dosage capacity
stau = 15;                % sensor time constant (min)
sgain = 0.5 / 1;          % sensor gain
%
snoxinbar = [ 8 ];        % influent substrate concentration (mg/l)
snoxbar = 1.0;            % effluent substrate concentration (mg/l)
ssinbar = 5;              % cod influent concentration (mg/l)
ssbar = 2;                % cod effluent concentration (mg/l)
xhbar = 1341.34;          % 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 = [ -1 ];
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
xhstep = 0;               % size of step biomass change
snoxstep = 0;             % size of step nitrate change
ssstep = 0;               % size of step substrate change
%----------------------------------------------------------
%
% determine simulations to run
%
n = length( kp );         % number of controllers
len = 480;                % length of simulation
%
yyy = [];
uuu = [];
%
for i = 1:n,
%
% initialise simulation   
  y = yo;             % initial output
  yd = yo + ydstep;   % setpoint
  u = uo;             % initial input
  cd = cdgain * uo;   % initial dosage
  xh = xhbar + xhstep;  % initial biomass change
  snox = snoxbar;     % initial nitrate
  snoxin = snoxinbar + snoxstep;  % initial substrate change
  ss = ssbar;     % initial nitrate
  ssin = ssinbar + ssstep;  % 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 > cdmax / cdgain, u = cdmax / cdgain; end;
    end;
    uu = [ uu; u ];
%
%   process
%
%   a) carbon dosage control dynamics (and simple euler integrator)
    dcddt = - cd / cdtau + u * cdgain / cdtau;
    cd = cd + dcddt * dt;
%   b) denitrification biomass respiration rate
    dur = muhat * xh * ...
    	   ( snox / ( knox + snox ) ) * ( ss / ( ks + ss ) );
%   c) nitrate mass balance (and simple euler integrator)
    dsnoxdt = snoxin / tau - snox / tau + dur * ( yh - 1 ) / yh / 2.86;
    snox = snox + dsnoxdt * dt;
%   d) carbon mass balance (and simple euler integrator)
    dssdt = ssin / tau - ss / tau - dur / yh + cd;
    ss = ss + dssdt * dt;
%   e) sensor dynamics (and simple euler integrator)
    dydt = - y / stau + snox * 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












