cola_paper.m

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

% CHECK:
% check stability by computing eigenvalues along the trajectory
Mwl=30;
Mw0 = Mwl*0.99 + 40*0.01; LTW = 2.70629*Mw0;   % Mass reflux constant
Vs=3.20629; Fs=1.0; zFs=0.5;
[A,B,C,D]=cola_linearize('colalv4mass_lin',Xinit',[LTW Vs Fs zFs]);
P = eig(A); Ps = sort(P);
[A2,B2,C2,D2]=cola_linearize('colalv4mass_lin',Xinit',[LTW Vs Fs zFs*0.98]);
P2 = eig(A2); P2s = sort(P2);

UstPo = [];
EUstPo = [];
PT = [];
for i=1:size(x,1)
  disp(['Point number: ', int2str(i)] );
  [Ah,Bh,Ch,Dh]=cola_linearize('colalv4mass_lin',...
           x(i,:),[LTW Vs Fs zFs*0.98]);
  Ph = eig(Ah); mPh = max(Ph);
  if mPh >= 0
    disp(['Point ', int2str(i),' is unstable: ' ])
    UstPo  = [UstPo; i];
    EUstPo = [EUstPo; Ph.']
  end
  Phs = sort(Ph);
  PT = [PT; Phs(1:10).'];
end

% save EigsOnTra UstPo EUstPo PT 
%
% NOTE: DOES INDEED BECOME UNSTABLE ALONG THER PATH FOR ML=30.
% END CHECK

%---------------------------------------------
% Figure 9: Nonlinear Simulation with other configurations (close level loops)
%---------------------------------------------

% -------------------------------------------------------------
% Simulate a change in feed rate with level loops open 
% -------------------------------------------------------------
% FIRST GO INTO THE FILE cola.m and change F=1.00 to F=1.01. 
% SAVE THE FILE as cola_F1.m_and simulate:

clear all
load Xinit

[t,x]=ode15s('cola4_F1',[0 500],Xinit); 
t0 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting

% Plot reboiler holdup (Apparent delay of about 1.26 min)
plot(t0,M1); axis([0 3 0.5 0.52]); title('MB: reboiler holdup')
plot(t0,xB); title('xB: reboiler composition')
plot(t0,yD); title('yD: distillate composition')

% LV-configuration
%FIRST GO INTO THE FILE cola_lv.m and change F to 1.01.  Then use 
[t,x]=ode15s('cola_lv_F1',[0 500],Xinit); 
tlv=t; M1lv=x(:,42); xBlv = x(:,1); yDlv = x(:,41); %save data for plot 
% Do the same for other configurations (REMEMBER TO CHANGE F!!!)
% LB-configuration
[t,x]=ode15s('cola_lb_F1',[0 500],Xinit); 
tlb=t; M1lb=x(:,42); xBlb = x(:,1); yDlb = x(:,41);
% Double ratio
[t,x]=ode15s('cola_rr_F1',[0 500],Xinit);
trr=t; M1rr=x(:,42); xBrr = x(:,1); yDrr = x(:,41); 
% DB
[t,x]=ode15s('cola_db_F1',[0 500],Xinit);
tdb=t; M1db=x(:,42); xBdb = x(:,1); yDdb = x(:,41); 
taul=0.063; LB = x(:,43)/taul;
plot(tdb,LB);

% Plot them together
plot(t0,M1,'-',tlv,M1lv,':',tlb,M1lb,'-.',trr,M1rr,'--'); 
title('MB: reboiler (bottom) holdup [mol]');
set(gca,'FontSize',14); xlabel('Time')
%hl = legend('OL','LV','LB','RR')
%set(hl,'Visible','off','Position', [0.15 0.68 0.15 0.25] )

plot(t0,xB,'-',tlv,xBlv,':',tlb,xBlb,'--',trr,xBrr,'--',tdb,xBdb,'--'); 
set(gca,'FontSize',14); xlabel('Time'); %ylabel('xB')
text(350,0.016,'OL','VerticalAlignment','top');
text(225,0.0151,'LV','VerticalAlignment','bottom',...
                    'HorizontalAlignment','right');
text(350,0.01,'RR','VerticalAlignment','bottom');
text(350,0.00695,'LB','VerticalAlignment','bottom');
text(350,0.0032,'DB','VerticalAlignment','bottom');
%print -deps fig5xb
%!mv fig5xb.eps  ~skoge/latex/work/fig/fig5xb.eps
%!mv fig5xb.eps ../latex/fig

plot(t0,yD,'-',tlv,yDlv,':',tlb,yDlb,'--',trr,yDrr,'--',tdb,yDdb,'--'); 
set(gca,'FontSize',14); xlabel('Time'); %ylabel('yD') 
text(350,0.99665,'DB','VerticalAlignment','top',...
                     'HorizontalAlignment','left');
text(350,0.9924,'OL','VerticalAlignment','bottom');
text(150,0.9918,'LV','VerticalAlignment','top');
text(350,0.9900,'RR','VerticalAlignment','top');
text(350,0.9849,'LB','VerticalAlignment','bottom');

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

%----------------------------------------------------------------
% Figure 12: Comparisons with the perfect operator. One-point control
%----------------------------------------------------------------

% 1. Feedback control one point. No meas. delay
clear all
L0 = 2.70629;
V0 = 3.20629;
TL = [0    L0;
      1000 L0];
TV = [0    V0;
      9.9  V0;
      10   V0*1.01
      1000 V0*1.01];
colas_PItop;
%----------------------------------------------------------
% IMPORTANT: YOU MUST NOW Start simulation in simulink
%----------------------------------------------------------
t1cr=t; M1cr = y4; xBcr = y2-0.01; yDcr = y1-0.99; 

% 2. Perfect operator at steady-state, use the value of L from 
% feedback control save as cola_lv_op.m
load Xinit
tsim = 500;
[t0,x]=ode15s('cola_lv_op',[0 tsim],Xinit); 
t0op = t0; M1op=x(:,42); xBop = x(:,1)-0.01; yDop = x(:,41)-0.99; 

close all

%PLOT TOP COMPOSITION
plot(t1cr,xBcr,'--',t0op,xBop,'-',t1cr,yDcr,'--',t0op,yDop,'-')
set(gca,'FontSize',14); xlabel('Time'); %ylabel('Compsitions')
hl = legend('--','LC','-','OP')
set(hl,'Visible','off','Position', [0.2 0.45 0.15 0.17] )
text(400,-0.79e-3,'xB')
text(400,1.2e-4,'yD')
%print -deps fig6c
%!mv fig6c.eps ~skoge/latex/work/fig/fig6.eps
%!mv fig6c.eps fig6c.eps ../latex/fig

% PLOT REFLUX:  u1 is L
DL = TL(1,2)+u1(size(u1,1),1)
plot(t1cr,u1+TL(1,2),'--',[0 9.9 10 tsim],[TL(1,2) TL(1,2) DL DL],'-')
set(gca,'FontSize',14); xlabel('Time'); %ylabel('Control value')
hl = legend('--','LC','-','OP')
set(hl,'Visible','off','Position', [0.2 0.45 0.15 0.17] )
%axis([0 tsim 0 0.05])

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


%-------------------------------------------------------------
% Figure 13:  Same with two-point control 
%-------------------------------------------------------------

% 1. First two-point control
clear all 
close all
L0 = 2.70629;
V0 = 3.20629;
TL = [0    L0;
      1000 L0];
TV = [0    V0;
      1000 V0];
rxB = [0    0.01;
       1000 0.01];
ryD = [0    0.99;
       9.9  0.99;
       10   0.995;
       1000 0.995];
colas_PItu1_nodelay   % with no time delay in measurements
% Simulation start, the simulation takes abot 2-4 min.
% Then:
t1cr=t; M1cr = y4; xBcr = y2-0.01; yDcr = y1-0.99; 

% 2. Now do the "open-loop change (Perfect operator)
load Xinit
tsim = 500;
[t,x]=ode15s('cola_lv_op2',[0 tsim],Xinit); 
t0op = t; M1op=x(:,42); xBop = x(:,1)-0.01; yDop = x(:,41)-0.99; 

close all
plot(t1cr,xBcr,'--',t1cr,yDcr,'--',t0op,yDop,'-',...
       t0op,xBop,'-',[0 tsim],[5e-3 5e-3],':');
set(gca,'FontSize',14); xlabel('Time'); %ylabel('Compsitions')
hl = legend('--','LC','-','OP')
set(hl,'Visible','off','Position', [0.3 0.53 0.1 0.17] )
text(400,-0.7e-3,'xB')
text(400,4e-3,'yD')
%print -deps fig7c
%!mv fig7c.eps  ~skoge/latex/work/fig/fig7c.eps

% To plot the inputs (L and V),

DL = TL(1,2)+u1(size(u1,1),1)
DL = TL(1,2)+u1(size(u1,1),1)
DV = TV(1,2)+u2(size(u1,1),1)
plot(t1cr,u1+TL(1,2),'--',t1cr,u2+TV(1,2),'--',...
     [0 9.9 10 tsim],[TL(1,2) TL(1,2) DL DL],'-',...
     [0 9.9 10 tsim],[TV(1,2) TV(1,2) DV DV],'-')
set(gca,'FontSize',14); xlabel('Time'); %ylabel('Control value')
hl = legend('--','LC','-','OP')
set(hl,'Visible','off','Position', [0.2 0.53 0.1 0.17] )
%axis([0 tsim 0 0.05])
text(400,DL+0.04,'L');
text(400,DV-0.02,'V','VerticalAlignment','top');

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


%------------------------------------------------------------
% Simulate same with time delay (not much effect; see also Fig. 18)
%------------------------------------------------------------

clear all 
colas_PI_dist
L0 = 2.70629;
V0 = 3.20629;
TL = [0    L0;
      1000 L0];
TV = [0    V0;
      1000 V0];
TF = [0    1;
      1000 1];
TzF= [0     0.5;
      1000  0.5];
rxB = [0    0.01;
%       99.9, 0.01;
%       100  0.005;
%       1000 0.005];
       1000 0.01];
ryD = [0    0.99;
       9.9  0.99;
       10   0.995;
       1000 0.995];
% Simulation start, the simulation takes abot 2-4 min.
% Then:
t1=t; M1 = y4; xB1 = y2-0.01; yD1 = y1-0.99; 
L1 = u1; V1 = u2;
close all
tsim=200;
figure(1)
plot(t1,xB1,'--',t1,yD1,'-',[0 tsim],[5e-3 5e-3],':',...
                           [0 tsim],[0 0],':');
text(150,0.00001,'xB','VerticalAlignment','bottom')
text(150,0.005-0.00001,'yD','VerticalAlignment','top');
set(gca,'FontSize',14); xlabel('Time'); ylabel('Compsitions')
axis([0 tsim -0.006 0.006])
%print -deps fig9tu1
%!mv fig9tu1.eps  ~skoge/latex/work/fig/fig9tu1.eps

%------------------------------------------------------------
% Figure 18: Simulation with disturbances and time delay
%------------------------------------------------------------
clear all 
L0 = 2.70629;
V0 = 3.20629;
TL = [0    L0;
      1000 L0];
TV = [0    V0;
      1000 V0];
rxB = [0    0.01;
       1000 0.01];
ryD = [0    0.99;
       199.9  0.99;
       200   0.995;
       1000  0.995];
TF = [0    1;
      9.9  1;
      10   1.2
      1000 1.2];

TzF= [0     0.5;
      99.9  0.5
      100   0.6
      1000  0.6];

colas_PI_dist
% Simulation start, the simulation takes abot 2-4 min.
% Then:
t1=t; M1 = y4; xB = y2-0.01; yD = y1-0.99; 
L = u1; V = u2;

subplot(2,1,1)
plot(t1,yD,t1,xB,'--')
set(gca,'FontSize',14); xlabel('Time'); ylabel('Compsitions')
text(30,2.2e-3,'xB','VerticalAlignment','bottom')
text(30,-3.2e-3,'xD','VerticalAlignment','top')
axis([0 300 -0.006 0.006])
%print -deps dist5yn
%!mv dist5yn.eps  ~skoge/latex/work/fig/dist5yn.eps



% -----------------------------------------------------
% Figure 10: Effect of slow level tuning for DV-configuration 
% -----------------------------------------------------

clear all 
load Xinit

tsim = 400;
% Increase V by 1 % (from 0.50 to 0.495)
% First for condenser level gain K=10
[t,x]=ode15s('coladv8_10',[0 tsim],Xinit);
t8 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD8 = yD(:) - 0.99; dxB8 = xB(:)-0.01;
plot(t8,dyD8,t8,dxB8,'--');
hold on

% Condenser level gain K=1
[t,x]=ode15s('coladv8_1',[0 tsim],Xinit);
t8 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD8 = yD(:) - 0.99; dxB8 = xB(:)-0.01;
plot(t8,dyD8,t8,dxB8,'--');

% Condenser level gain K=0.1
[t,x]=ode15s('coladv8_01',[0 tsim],Xinit);
t8 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
dyD8 = yD(:) - 0.99; dxB8 = xB(:)-0.01;