%
% Feedback Control Simulation
%
% Process: denitrification
% Controller: Predictive
%
% Bob Newell, November 1997
%
%----------------------------------------------------------
function dxdt = dmc_function( t, x )
%
global dynmat kpred predvect delt nextt nextu lastu law
%
% 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
ks = 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
%
yo = 0.5;                 % initial process output
uo = 0.5;                 % initial process input
%
% open loop - no control
if law == 0,
  kp = 0;
  ki = 0;
  uo = 0.6;
end; 
% closed loop - P control test
if law == 1,
  tm = 50;
  kp = - 150 / 5 / tm;
  ki = - 1 / 5 / tm;
end; 
%
%----------------------------------------------------------
%
% Define the disturbance
%
ydstep = 0;               % size of setpoint change
if t > 250, ydstep = -0.25; end;
xhstep = 0;               % size of step biomass change
snoxstep = 2.5;           % size of step nitrate change
ssstep = 0;               % size of step substrate change
%----------------------------------------------------------
%
% initialise simulation   
%
yd = yo + ydstep;   % setpoint
xh = xhbar + xhstep;  % initial biomass change
snoxin = snoxinbar + snoxstep;  % initial substrate change
ssin = ssinbar + ssstep;  % initial substrate change
%
%   PI controller
%
e = yd - x(4);
if law < 2,
  u = uo + kp * e + ki * x(5);
  if u < 0, u = 0; end;
  if u > cdmax / cdgain, u = cdmax / cdgain; end;
end;
%
%   predictive controller
%
if law == 2,
  if t < nextt,
    u = lastu;
  elseif t > nextt + delt,
    lastu = nextu;
% update the prediction vector
    np = length( predvect );
    e = x(4) - predvect(1);
    predvect(1:np-1) = predvect(2:np);
    predvect = predvect + e;
% calculate the error vector
    err = yd*ones(np,1) - predvect;
% calculate optimal control action
    delu = kpred * err;
% calculate new prediction vector
    nu = length( delu );
    predvect = predvect + dynmat * [ delu(1); zeros(nu-1,1) ];
% calculate actual control action
    u = nextu + delu(1);
    nextu = u;
    nextt = nextt + delt;
  elseif t >= nextt,
    u = nextu;
  end;
end;
%
%   process
%
dxdt = zeros( 5, 1 );
%   a) carbon dosage control dynamics (and simple euler integrator)
dxdt(1) = - x(1) / cdtau + u * cdgain / cdtau;
%   b) denitrification biomass respiration rate
dur = muhat * xh * ...
    	   ( x(2) / ( knox + x(2) ) ) * ( x(3) / ( ks + x(3) ) );
%   c) nitrate mass balance (and simple euler integrator)
dxdt(2) = snoxin / tau - x(2) / tau + dur * ( yh - 1 ) / yh / 2.86;
%   d) carbon mass balance (and simple euler integrator)
dxdt(3) = ssin / tau - x(3) / tau - dur / yh + x(1);
%   e) sensor dynamics (and simple euler integrator)
dxdt(4) = - x(4) / stau + x(2) * sgain / stau;
%   f) controller integral of error
dxdt(5) = e;
%
%----------------------------------------------------------
% the end
