cola_paper.m

国外经典书籍MULTIVARIABLE FEEDBACK CONTROL-多变量反馈控制 的源码

% This file is used to generate the examples and figures in:
% S. Skogestad: Dynamics and control of distillation columns -
% A tutorial introduction. Presented at Distillationa and Absorbtion
% 97, Maastricht, Netherlands, 8-10 Sept. 1997

%---------------------------------------------
% First find steady-state 
%---------------------------------------------

% Do this by simulating 20000 min with stabilized LV-model:
[t,x]=ode15s('cola_lv',[0 20000],0.5*ones(1,82)'); 
Xinit = x(sel(size(x),1,1),:)'; % This data has been saved in cola_init.mat
                               % and can be retrieved using the command
                               % load cola_init
  
% Perform nonlinear simulation
[t,x]=ode15s('cola4',[0 500],Xinit);   
plot(t,x);  %Nothing happens so we are at steady-state...
%save Xinit Xinit
    
% ---------------------------------------------
% Figure 2: Composition profiles
% -------------------------------------------

xprofile0 = Xinit(1:41);
stages = [1:1:41]';

% Now after change in external flows (increase L by 0.02 with V constant)
[t,x]=ode15s('colalv1large',[0 20000],Xinit); 
xdum = x(sel(size(x),1,1),:)';
xprofile1 = xdum(1:41);

% Now after change in internal flows (increase both L and V by 1)
[t,x]=ode15s('colalv2large',[0 20000],Xinit); 
xdum = x(sel(size(x),1,1),:)';
xprofile2 = xdum(1:41);

plot(xprofile0,stages); hold
plot(xprofile1,stages,'-.'); 
plot(xprofile2,stages,'--'); hold off;
axis([0 1 1 41])
ylabel('Stage no.')
xlabel('Composition')
hl = legend('Nom','Ext','Int')
set(hl,'Visible','off','Position', [0.15 0.7 0.2 0.2] )
%print -deps profile
%!mv profile.eps ~skoge/latex/work/fig/profile.eps

% Logarithmic compositions
Xprofile0 = log(xprofile0(:)./(1-xprofile0(:)) );
Xprofile1 = log(xprofile1(:)./(1-xprofile1(:)) );
Xprofile2 = log(xprofile2(:)./(1-xprofile2(:)) );

plot(Xprofile0,stages); hold;
plot(Xprofile1,stages,'-.'); 
plot(Xprofile2,stages,'--'); hold off
axis([-6 6 1 41 ])
ylabel('Stage no.')
xlabel('Logaritmic composition')
hl = legend('Nom','Ext','Int')
set(hl,'Visible','off','Position', [0.15 0.7 0.2 0.2] )
%print -deps profilelog
%!mv profilelog.eps ~skoge/latex/work/fig/profilelog.eps
hold off

% ------------------------------------------------
% Figure 4: Respons to bottom liquid flow L_B = L_2
% ------------------------------------------------
clear all
load Xinit
tsim=7;
[t2,x]=ode23s('colalv2',[0:0.01:tsim],Xinit,[1e-5]); 
[t2n,xn]=ode15s('colalv2n',[0:0.01:tsim],Xinit); 
    taul=0.063;     	% time constant for liquid dynamics (min)

LB = x(:,43)/taul; % This is L - L0i + M0i/taul; assumes lambda=0
LT=2.70629*1.1;                          % Reflux

plot(t2,LB)

t2 = [-1; t2];
L = [LB(1);LB];
LTv = [ 0;0;LT-LT/1.1;LT-LT/1.1];
tv  = [-1;0; 0;tsim];

plot(t2,L-L(1),tv,LTv,'--');
xlabel('Time'); ylabel('Change in L');
hl = legend('LB','LT')
set(hl,'Visible','off','Position', [0.7 0.475 0.11 0.2] )
set(gca,'FontSize',14)
axis([-1 7 -0.05 0.3])
%print -deps figlv2a
%!mv figlv2a.eps ~skoge/latex/work/fig/figlv2a.eps
%!mv figlv2a.eps ../latex/fig

%-----------------------------------------------------------------
% (NO FIGURE): K2 effect (parameter lambda). Incerease V by 0.1
% ---------------------------------------------------------------a
tsim=7;
!cp colamodk2_0.m colamodk2.m
[t0,x]=ode15s('cola4k2',[0 tsim],Xinit); % lambda=0
M10=x(:,42); xb0 = x(:,1);

!cp colamodk2_05.m colamodk2.m
[t1,x]=ode15s('cola4k2',[0 tsim],Xinit); % lambda=0.5
M11=x(:,42); xb1=  x(:,1);

!cp colamodk2_1.m colamodk2.m
[t2,x]=ode15s('cola4k2',[0 tsim],Xinit); % lambda=1
M12=x(:,42); xb2 =  x(:,1);

!cp colamodk2_15.m colamodk2.m
[t3,x]=ode15s('cola4k2',[0 tsim],Xinit); % lambda=1.5
M13=x(:,42); xb3 =  x(:,1);

!cp colamodk2_m5.m colamodk2.m
[t4,x]=ode15s('cola4k2',[0 tsim],Xinit); % lambda=-0.5
M14=x(:,42); xb4 =  x(:,1);

plot(t0,M10,t1,M11,t2,M12,t3,M13,t4,M14); % Inverse response for lambda>1
plot(t0,xb0,t1,xb1,t2,xb2,t3,xb3,t4,xb4); % nothing strange for lambda=0.5

t2 = [-1; t2];
L = [LB(1);LB];
LTv = [ 0;0;LT-LT/1.1;LT-LT/1.1];
tv  = [-1;0; 0;tsim];


% -------------------------------------------------------------
% Figure 5: LV-configuration Simulate a change in reflux rate with V constant 
% -------------------------------------------------------------
% FIRST GO INTO THE FILE colalv.m and increase L by 0.1%
% SAVE THE FILE as colalv1.m simulate:
clear all
load Xinit
tsim=500;
[t,x]=ode15s('colalv1',[0 tsim],Xinit); 
t1 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD = yD(:) - 0.99; dxB = xB(:)-0.01;
% 1n. same without flow dynamics
[t,x]=ode15s('colalv1n',[0 tsim],Xinit);
t1n = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD1n = yD(:) - 0.99; dxB1n = xB(:)-0.01;
% Now do the same change in boilup
[t,x]=ode15s('colalv1v',[0 tsim],Xinit); 
t1v = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyDv = yD(:) - 0.99; dxBv = xB(:)-0.01;
% 1n. same without flow dynamics
[t,x]=ode15s('colalv1vn',[0 tsim],Xinit);
t1vn = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD1vn = yD(:) - 0.99; dxB1vn = xB(:)-0.01;

plot(t1,dxB,'--',t1,dyD,t1n,dxB1n,':',t1n,dyD1n,':');
hold on
plot(t1v,dxBv,'--',t1v,dyDv,t1vn,dxB1vn,':',t1vn,dyD1vn,':'); 
%title('0.1% increase in external flows'); xlabel('Time (min)')
set(gca,'Fontsize',14)
xlabel('Time')
ylabel('Composition')
hl = legend('xB','yD','NoH')
set(hl,'Visible','off','Position', [0.15 0.7 0.2 0.2] )
text(200,1e-3,'Increase in L with V constant');
text(200,-1e-3,'Increase in V with L constant');
%print -deps figlv1
%!mv figlv1.eps ~skoge/latex/work/fig/figlv1.eps
hold off


% -------------------------------------------------------------
% Figure 6: LV-configuration Simulate a change in  internal flows
% -------------------------------------------------------------
% Increse L by 10% and increase V by same amount som D is constant
% SAVE THE FILE as colalv2.m simulate:
clear all
load Xinit
tsim=500;
[t,x]=ode15s('colalv2',[0 tsim],Xinit); 
t2 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD = yD(:) - 0.99; dxB = xB(:)-0.01;


%plot(t2,dxB,t2,dyD); title('10% increase in internal flows') 
% 2n. Do the same without flow dynamics
[t,x]=ode15s('colalv2n',[0 tsim],Xinit); 
t2n = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD2n = yD(:) - 0.99; dxB2n = xB(:)-0.01;

plot(t2,dxB,'--',t2,dyD,t2n,dxB2n,':',t2n,dyD2n,':'); 
%title('10% increase in internal flows'); xlabel('Time (min)') 
set(gca,'Fontsize',14)
xlabel('Time')
ylabel('Composition')
hl = legend('xB','yD','NoH')
set(hl,'Visible','off','Position', [0.15 0.475 0.2 0.2] )
%print -deps figlv2
%!mv figlv2.eps ~skoge/latex/work/fig/figlv2.eps

%-----------------------------------------------
% Linearized model for LV- configuration
%------------------------------------------------

clear all
load Xinit
%The linear model for the LV-configurtation (4 inputs, 2 outputs)
Ls=2.70629; Vs=3.20629; Fs=1.0; zFs=0.5;
[A,B,C,D]=cola_linearize('cola_lv_lin',Xinit',[Ls Vs Fs zFs]);
Glvu =  pck(A,B,C,D);
eig(A)
% compare with model with no flow dynamics
[A,B,C,D]=cola_linearize('cola_lvn_lin',Xinit',[Ls Vs Fs zFs]);
eig(A)

%-----------------------------------------------
% Figure 7: Compare linear and nonlinear simulations
%------------------------------------------------
tsim=30;

% Nonlinear of magnitude 0.1%
[t,x]=ode15s('colalv1',[0 tsim],Xinit);
t3 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD3 = (yD(:) - 0.99)/(2.70629*0.001);
YD = log(yD./(1-yD)); dYD3 = (YD(:) - YD(1))/(2.70629*0.001);

% simulate linear step in L of magnitude 1
t3l = linspace(0,tsim,10*tsim); zero = zeros(10*tsim,1);
[Y,X] = step(A,B,C,D,1,t3l);
dyDl = Y(:,1);

% Nonlinear of magnitude 1%
[t,x]=ode15s('colalv1_1',[0 tsim],Xinit);
t3_1 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD3_1 = (yD(:) - 0.99)/(2.70629*0.01);
YD = log(yD./(1-yD)); dYD3_1 = (YD(:) - YD(1))/(2.70629*0.01);

% Nonlinear of magnitude 10%
[t,x]=ode15s('colalv1_10',[0 tsim],Xinit);
t3_10 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD3_10 = (yD(:) - 0.99)/(2.70629*0.1);
YD = log(yD./(1-yD)); dYD3_10 = (YD(:) - YD(1))/(2.70629*0.1);

% Nonlinear of magnitude 50%
[t,x]=ode15s('colalv1_50',[0 tsim],Xinit);
t3_50 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD3_50 = (yD(:) - 0.99)/(2.70629*0.5);
YD = log(yD./(1-yD)); dYD3_50 = (YD(:) - YD(1))/(2.70629*0.5);

% Plotting

plot(t3,dyD3); hold on
plot(t3l,dyDl,'--'); 
plot(t3l,zero,':');
plot(t3_1,dyD3_1); 
plot(t3_10,dyD3_10); 
plot(t3_50,dyD3_50); 
hold off
text(25,0.112,'lin','HorizontalAlignment','left',...
                    'VerticalAlignment','bottom');
text(25,0.104,'+0.1','HorizontalAlignment','left','VerticalAlignment','top');
text(25,0.09,'+1','HorizontalAlignment','left','VerticalAlignment','top');
text(25,0.034,'+10','HorizontalAlignment','left',...
                    'VerticalAlignment','bottom');
text(25,0.007,'+50','HorizontalAlignment','left',...
                    'VerticalAlignment','bottom');

%title('Nonlinearity xD for increase in L');  xlabel('Time (min)')
set(gca,'FontSize',14,'Position',[0.13 0.11 0.775 0.8]); 
set(gcf,'Position',[296 420 560 420],'PaperPosition',[0.25 2.5 8 6])
xlabel('Time'); ylabel('Distillate composition')
%print -deps figlv3
%!mv figlv3.eps ~skoge/latex/work/fig/figlv3.eps

plot(t3,dYD3); hold;
plot(t3_1,dYD3_1);
plot(t3_10,dYD3_10);
plot(t3_50,dYD3_50);
text(25,10.9,'+0.1','HorizontalAlignment','left',...
                     'VerticalAlignment','bottom');
text(25,10.5,'+1','HorizontalAlignment','left','VerticalAlignment','top');
text(25,9.1,'+10','HorizontalAlignment','left',...
                    'VerticalAlignment','top');
text(25,6.8,'+50','HorizontalAlignment','left',...
                    'VerticalAlignment','bottom');
%title('Nonlinearity with log compositions');  xlabel('Time (min)')
set(gca,'FontSize',14,'Position',[0.13 0.11 0.775 0.8]); 
set(gcf,'Position',[296 420 560 420],'PaperPosition',[0.25 2.5 8 6])
xlabel('Time'); ylabel('Distillate composition')
%print -deps figlv3log
%!mv figlv3log.eps  ~skoge/latex/work/fig/figlv3log.eps
hold;

%---------------------------------------------
% Figure 8: Effect of mass flows
%---------------------------------------------
clear all
load Xinit
tsim = 1000;

% Decrease zf by 1 % (from 0.50 to 0.495)
[t,x]=ode15s('colalv4',[0:1:tsim],Xinit);
t4 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD4 = yD(:) - 0.99; dxB4 = xB(:)-0.01;
plot(t4,dyD4,t4,dxB4); 
hold

% Same with MASS reflux constant
[t,x]=ode15s('colalv4mass',[0:1:tsim],Xinit);
%[t,x]=ode15s('dum',[0:1:tsim],Xinit);
%[t,x]=ode23s('colalv4mass',[0:1:tsim],Xinit);
%[t,x]=ode45('colalv4mass',0, 1000, Xinit);
t4mass = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD4mass = yD(:) - 0.99; dxB4mass = xB(:)-0.01;

plot(t4mass,dyD4mass,'--',t4mass,dxB4mass,'--'); 
set(gca,'FontSize',14); xlabel('Time'); ylabel('Composition')
hl = legend('--','MR','-','MF')
set(hl,'Visible','off','Position', [0.15 0.12 0.15 0.18] )
text(400,-0.003,'xB'); text(400,-0.02,'yD');
hold off;
%title('Effect of mass reflux');  xlabel('Time (min)')
axis([0 1000 -0.05 0] )

%print -deps figlv4
%!mv figlv4.eps  ~skoge/latex/work/fig/figlv4.eps