bump.m

"Teaching antiwindup, bumpless transfer, and split-range control"一文的仿真程序实现

% Copyright c 1998-2002 The Board of Trustees of the University of Illinois
% 		  All rights reserved.
% Developed by:	Large Scale Systems Research Laboratory
%		Professor Richard Braatz, Director
%               Department of Chemical Engineering
%       	University of Illinois
%		http://brahms.scs.uiuc.edu
% 
% Permission hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to
% deal with the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the 
% Software is furnished to do so, subject to the following conditions:
% 		1. Redistributions of source code must retain the above copyright
%		   notice, this list of conditions and the following disclaimers.
%		2. Redistributions in binary form must reproduce the above 
%		   copyright notice, this list of conditions and the following 
%		   disclaimers in the documentation and/or other materials 
%		   provided with the distribution.
%		3. Neither the names of Large Scale Research Systems Laboratory,
%		   University of Illinois, nor the names of its contributors may
%		   be used to endorse or promote products derived from this 
%		   Software without specific prior written permission.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 
% THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
% OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 
% ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
% DEALINGS IN THE SOFTWARE.
%
%
%	Benchmark for studying antiwindup and bumpless transfer
%
%	Based on experimental apparatus used for collecting
%	diffusion data described in Phillip Drake's PhD thesis in chemistry.
%	The control algorithm is based on two IMC-PID controllers
%	implemented in velocity form described in the manuscript:
%
%	S.H. Chung and R.D. Braatz. Teaching antiwindup, bumpless
%	transfer, and split-range control. Chemical Engineering
%	Education, 32:222-223, 1998.
%
%	See the paper for the precise details of the problem statement.
%	Users of this code should cite the above paper.
%
%	Authors:  	      Serena H. Chung and Richard D. Braatz
%	Date last modified:   January 2, 2002
%	www:		      http://brahms.scs.uiuc.edu/lssrl/software
%	e-mail:		      [e-mail]
%

clear
clear global
format long
format compact
clf

delt =  0.1 ;                           % initialize step size
tau = [9.5 1.7] ;                       % process time constants
kp = [1.0 0.068] ;			% process gains
delay = [24 14] ;                       % set this to an integer value
theta = [delay(1)*delt delay(2)*delt] ; % process time delays
finalk =  7000 ;                        % set this to an integer value
T = finalk ;				% final simulation time
finaltime = finalk * delt ;
beta = 0.013;                       	% set measurement noise level

y = zeros(finalk+2,1) ;			% initialization
e = zeros(finalk+1,1) ;			% setting vector dimensions
d = zeros(finalk+2,1) ;			% in this way speeds MATLAB
u1 = zeros(1,finalk+1) ;		% calculations
u2 = zeros(1,finalk+1) ;	
h1 = zeros(finalk+T+25,1) ;		
h2 = zeros(finalk+T+15,1) ;
a1 = zeros(finalk+T+1,1) ;
a2 = zeros(finalk+T+1,1) ;

d(1) = beta*randn(1) ;			% disturbance
for i = 2:(finalk)+2,
	d(i) = d(i-1) + beta*randn(1) ;
end

disp('  ')
disp('Hit enter to use the IMC-PID parameters in Handout B.')
whichpidparam = input('Enter 1 to enter your own IMC-PID parameters: ') ;
disp('  ')
if isempty(whichpidparam),			% use default IMC-PID controller
   disp('   ')
   disp('Either enter the IMC filter time constant (lambda) or')
   lambda = input('hit return to use the default (lambda=1):') ;
   disp('  ')
   if isempty(lambda), lambda = 1.0 ; end
   for i = 1:2,                       	
        kc(i) = (2*tau(i)+theta(i))/2/kp(i)/(lambda+theta(i)) ;
        taui(i) = tau(i) + theta(i)/2 ;
        taud(i) = tau(i)*theta(i)/(2*tau(i)+theta(i)) ;
	tauf(i) = lambda*theta(i)/2/(lambda+theta(i)) ;
   end
else
   for i = 1:2, 
        kc(i) = input([ 'Enter the controller gain (kc) for P', ...
		   	   num2str(i), ': ']) ;
        taui(i) = input([ 'Enter the integral time constant (tauI) for P', ...
		   	   num2str(i), ': ']) ;
        taud(i) = input([ 'Enter the derivative time constant (tauD) for P', ...
		   	   num2str(i), ': ']) ;
        tauf(i) = input([ 'Enter the filter time constant (tauF) for P', ...
		   	   num2str(i), ': ']) ;
   end
end

alpha1 = tauf./delt./(1+tauf./delt) ;
alpha2 = kc./(1+tauf./delt) ;
alpha3 = taud./delt ;
alpha4 = delt./taui ;

y0 = 21.0 ;   			    	% set initial temperature
setpoint = ones(finalk+T+1,1) ;		% set desired temperature

for jj=1:1301+700,
	setpoint(jj) = 120.0 - y0 ;
end
for jj=1302+700:5018+700,
	setpoint(jj) = -40.0/(2858.0-1287.0)*(jj-2858.0-700.0)+80.0-y0 ;
end
for jj=5019+700:9001+700,
	setpoint(jj) = 25.0 - y0;
end

			% calculate coefficients for impulse response model
for ii = 1:(finalk+T+1),
	h1(delay(1)+ii) = kp(1)*(exp(-(ii-1) ...
               		* delt/tau(1))-exp(-ii*delt/tau(1))) ;
	h2(delay(2)+ii) = kp(2)*(exp(-(ii-1) ...
               		* delt/tau(2))-exp(-ii*delt/tau(2))) ;
end

a1(1) = h1(1) ;	        % calculate coefficients for step response model
a2(1) = h2(1) ;

for ii = 2:(finalk+T+1),
	a1(ii) = h1(ii) + a1(ii-1) ;
	a2(ii) = h2(ii) + a2(ii-1) ;
end

% p os a flag that indicates the controller and process to use
% u1 is the manipulated variable for the process with heat sink 1
% u2 is the manipulated variable for the process with heat sink 2
% y1 is the temperature response to the u1 control moves 
% y2 is the temperature response to the u2 control moves

if (setpoint(1)+y0) < 30,
	u1(1) = 0.0 ; p = 2 ;
	u2(1) = kc(p)/(1+tauf(p)/delt)*(setpoint(1) - y(1) ...
	      + delt/taui(p)*(setpoint(1)-y(1))) ;
	u1(2) = 0.0 ;  
	u2(2) = u2(1) + tauf(p)/delt/(1+tauf(p)/delt)*u2(1) ...
       	      + kc(p)/(1+tauf(p)/delt)*(setpoint(2) - y(2) ...
	      - setpoint(1) + y(1) + delt/taui(p)*(setpoint(2) - y(2))) ;
else
	u2(1) = 0.0 ; p = 1 ;
	u1(1) = kc(p)/(1+tauf(p)/delt)*(setpoint(1) - y(1) ...
	      + delt/taui(p)*(setpoint(1)-y(1))) ;
	u2(2) = 0.0 ; 
	u1(2) = u1(1) + tauf(p)/delt/(1+tauf(p)/delt)*u1(1) ...
       	      + kc(p)/(1+tauf(p)/delt)*(setpoint(2) - y(2) ...
	      - setpoint(1) + y(1) + delt/taui(p)*(setpoint(2) - y(2))) ;
end

if u1(1) > 100.0, u1(1) = 100.0 ; end
if u1(2) > 100.0, u1(2) = 100.0 ; end
if u2(1) > 100.0, u2(1) = 100.0 ; end
if u2(2) > 100.0, u2(2) = 100.0 ; end
if u1(1) < 0.0, u1(1) = 0.0 ; end
if u1(2) < 0.0, u1(2) = 0.0 ; end
if u2(1) < 0.0, u2(1) = 0.0 ; end
if u2(2) < 0.0, u2(2) = 0.0 ; end

for kk = 1:finalk,
    if ((round(kk/100)-kk/100) == 0),
	disp(['time = ' num2str(round(kk/10))]) 
    end
    e(kk+1) = setpoint(kk+1) - y(kk+1) ;
    if kk > 1,
	if (setpoint(kk)+y0) < 30,
		u1(kk+1) = 0.0 ; p = 2 ; 	
		u2(kk+1) = u2(kk) + alpha1(p)*(u2(kk)-u2(kk-1)) ...
            		 + alpha2(p)*(e(kk+1)-e(kk) ...
			 - alpha3(p)*(y(kk+1)-2*y(kk)+y(kk-1)) ...
			 + alpha4(p)*e(kk+1)) ;
	else
		u2(kk+1) = 0.0 ; p = 1 ; 	
		u1(kk+1) = u1(kk) + alpha1(p)*(u1(kk)-u1(kk-1)) ...
            		 + alpha2(p)*(e(kk+1)-e(kk) ...
			 - alpha3(p)*(y(kk+1)-2*y(kk)+y(kk-1)) ...
			 + alpha4(p)*e(kk+1)) ;
	end

	if u1(kk+1) > 100.0, u1(kk+1) = 100.0; end
	if u1(kk+1) < 0.0, u1(kk+1) = 0.0; end
	if u2(kk+1) > 100.0, u2(kk+1) = 100.0; end
	if u2(kk+1) < 0.0, u2(kk+1) = 0.0; end
    end
    sum = (u1(kk+1:-1:2) - u1(kk:-1:1))*a1(1:kk) ...
	+ (u2(kk+1:-1:2) - u2(kk:-1:1))*a2(1:kk) ...
	+ a1(kk+1)*u1(1) + a2(kk+1)*u2(1) ;

    y(kk+2) = sum + d(kk+2) ;		% add noise
end

Y = y + y0 ;
transition = 30.0*ones(1,finalk) ;

clf						% plot the output
subplot(221)
plot(delt*[1:finalk]', setpoint(1:finalk,1)+y0, '-')
	xlabel('Time (minutes)')
	ylabel('Temperature (Centigrade)')
	axis([0 700 0 140])
	text(330,-37,'(a)') ;
subplot(222)
plot(delt*[1:finalk]', Y(1:finalk,1), '-', ...
	delt*[1:finalk]', transition(1:finalk), '--')
	xlabel('Time (minutes)')
	ylabel('Temperature (Centigrade)')
	axis([0 700 0 140])
	text(330,-37,'(b)') ;
subplot(223)
plot(delt*[1:finalk]', u1(1,1:finalk), '--', ...
     	delt*[1:finalk]', u2(1,1:finalk), '-')
	xlabel('Time (minutes)')
	ylabel('Percent power')
	axis([0 700 0 100])
	text(330,-25,'(c)') ;
subplot(224)
plot(delt*[700:finalk]',setpoint(700:finalk,1)+y0'-Y(700:finalk,1),'-')
	xlabel('Time (minutes)')
	ylabel('Temperature (Centigrade)')
	axis([0 700 -1.5 1.0])
	text(330,-2.125,'(d)') ;