% 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: braatz@uiuc.edu % 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)') ;