%
% Averaging Level Control Simulation
%
% Process: tank with level sensor and valve
% Controller: PID
%
% Bob Newell, April 1996
%
%----------------------------------------------------------
%
% Process Parameters
%
area = 100;               % tank cross-sectional area (m2)
vtau = 3;                 % valve time constant (min)
cvmax = 14.142;           % valve size coefficient
linear = 1;               % linear or equal percentage
stau = 2;                 % sensor time constant (min)
sgain = 0.5 / 2.5;        % sensor gain
%
fbar = [ 50 ];            % flowrate
hbar = 2.0;               % initial height
pbar = 100;               % pump delivery pressure (kPa)
%
%----------------------------------------------------------
%
% Controller Parameters
%     u = uo + kp*e + ki*integral(e) + kd*derivative(e)
%
% open loop - no control
%kp = [ 0 ];
%ki = [ 0 ];
%kd = [ 0 ]; 
% 
kp = [ -1.00  -0.1   ];
ki = [ -0.05  -0.005 ];
kd = [ -2.00  -0.2   ];
%
dt = 0.5;                 % sampling interval
if linear,
   uo = fbar ./ cvmax ./ ...
        sqrt( pbar - 50 );  % initial controller output
else,
   uo = 1 + log( fbar ./ cvmax ./ 2.236 ./ ...
        sqrt( pbar - 50 ) ) ./ ...
	    log( 20 );            % initial controller output
end;
yo = hbar * sgain;            % initial process output
%
%----------------------------------------------------------
%
% Define the disturbance
%
fmag = 40;                % size of flow changes
flen = 10;                % length of pulses
fdur = 50;                % time between pulses
%----------------------------------------------------------
%
% determine simulations to run
%
n = length( kp );         % number of controllers
l = 1000;                  % length of simulation
%
yyy = [];
uuu = [];
%
for i = 1:n,
%
% initialise simulation   
      y = yo;             % initial output
      yd = 0.4;           % setpoint
      u = uo;             % initial input
	  uptime = 10;
	  downtime = flen + 10;
	  h = hbar;
	  fi = fbar;
	  fo = fbar;
%
      esum = 0;           % initial integral error sum
      laste = 0;          % initial last error
%
      yy = [];
      uu = [];
%
      for k = 1:l,
%
%        controller
         e = yd - y;
         esum = esum + e * dt;
         dedt = ( e - laste ) / dt;
         laste = e;
         u = uo + kp(i) * e + ki(i) * esum + kd(i) * dedt ;
         uu = [ uu; [ u fi/fbar/2 ] ];
%
%        process
         if linear,
            f = cvmax * u * sqrt( pbar - 50 );
         else,
		    f = cvmax * ( 2.236 * 20 ^ ( u - 1 ) ) ...
			                        * sqrt( pbar - 50 );
		 end;
         dfodt = - fo / vtau + f / vtau;
		 if (k*dt) > uptime,
		    fi = fbar + fmag;
			uptime = uptime + fdur;
		 end;
		 if (k*dt) > downtime,
		    fi = fbar;
			downtime = downtime + fdur;
		 end;
         dhdt = ( fi - fo ) / area;
         dydt = - y / stau + h * sgain / stau;
%
         fo = fo + dfodt * dt;
         h = h + dhdt * dt;
         y = y + dydt * dt;
%
         yy = [ yy; y ];
%
      end;
%
      yyy = [ yyy yy ];
	  uuu = [ uuu uu ];
%
end;
%
h = figure;
subplot( 2, 1, 1 );
plot( dt*(600:1000)', uuu(600:1000,:), '-' );
set( gca, 'fontsize', 16 );
xlabel( 'Time in minutes' );
ylabel( 'Process inputs' );
% book figures use
axis( [ 300 500 0.40 1.1 ] );
arrow( [ 340 0.95 ], [ 326 0.75 ] );
h = text( 325, 1, 'conventional' );
set( h, 'fontsize', 16 );
arrow( [ 390 0.95 ], [ 372 0.75 ] );
h = text( 375, 1, 'disturbance' );
set( h, 'fontsize', 16 );
arrow( [ 440 0.95 ], [ 435 0.61 ] );
h = text( 425, 1, 'averaging' );
set( h, 'fontsize', 16 );
%
subplot( 2, 1, 2 );
plot( dt*(600:1000)', yyy(600:1000,:), '-' );
set( gca, 'fontsize', 16 );
xlabel( 'Time in minutes' );
ylabel( 'Process outputs' );
% book figures use
axis( [ 300 500 0.0 1.0 ] );
arrow( [ 340 0.83 ], [ 355 0.40 ] );
h = text( 325, 0.9, 'conventional' );
set( h, 'fontsize', 16 );
arrow( [ 440 0.83 ], [ 435 0.58 ] );
h = text( 425, 0.9, 'averaging' );
set( h, 'fontsize', 16 );
%
print -f1 -deps2 fig17-12.eps
%
%----------------------------------------------------------
% the end
