cola_paper.m

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

plot(t8,dyD8,t8,dxB8,'--');
hold off

%title('Effect of level tuning for DV');  
xlabel('Time')
axis([0 tsim -20e-4 6e-4])
text(50,3.3e-4,'K=10')
text(50,-5.4e-4,'K=10')
text(150,2.7e-4,'K=1','VerticalAlignment','top')
text(150,-6.1e-4,'K=1','VerticalAlignment','top')
text(250,-0.8e-4,'K=0.1','VerticalAlignment','top')
text(250,-10e-4,'K=0.1','VerticalAlignment','top')
hl = legend('yD','xB')
set(hl,'Visible','off','Position', [0.75 0.1 0.12 0.2] )

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

%------------------------------------------------------------
%  Linear analysis (Freq.dep. RGA etc.)
% ---------------------------------------------------------

clear all
load Xinit

% LV-configurations
Ls=2.706; Vs=3.206; 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);  % Model has 4 inputs and 2 outputs
% Scale the model:
Si = diag([1 1 .2 .1]); So = diag([1/0.01 1/0.01]); %scalings
Glv = mmult (So, Glvu, Si);           

%For the double ratio (L/D-V/B) configuration (with 4 inputs and 2 outputs):
R1s = 2.70629/0.5; R2s=3.20629/0.5; Fs=1.0; zFs=0.5;
[A,B,C,D]=cola_linearize('cola_rr_lin',Xinit',[R1s R2s Fs zFs]);
Grru =  pck(A,B,C,D);
Si = diag([1 1 .2 .1]); So = diag([1/0.01 1/0.01]); %scalings
Grr = mmult (So, Grru, Si);           

% We could alternatively start directly from the linearized 4x4 model
% I'll do that for the rest
Ls=2.706; Vs=3.206; Ds=0.5; Bs=0.5; Fs=1.0; zFs=0.5;
[A,B,C,D]=cola_linearize('cola4_lin',Xinit',[Ls Vs Ds Bs Fs zFs]);
G4u =  pck(A,B,C,D); 
Si = diag([1 1 1 1 .2 .1]); So = diag([1/0.01 1/0.01 1 1]); %scalings
G4 = mmult (So, G4u, Si);               % scaled 82 state model

% Use sysic to generate DV-configuration
Kl = 10; Kb = 10; % P-controllers for levels (bandwidth = 10 rad/min)
systemnames = 'G4 Kl Kb';
inputvar = '[D(1); V(1); d(2)]'; outputvar = '[G4(1);G4(2)]';
input_to_G4 = '[Kl; V; D; Kb; d ]';
input_to_Kl = '[G4(3)]'; input_to_Kb = '[G4(4)]';
sysoutname ='Gdv'; cleanupsysic='yes'; sysic;

% DB-configuration
Kl = 10; Kv = 10;
systemnames = 'G4 Kl Kv';
inputvar = '[D(1); B(1); d(2)]'; outputvar = '[G4(1);G4(2)]';
input_to_G4 = '[Kl; Kv; D; B; d ]';
input_to_Kl = '[G4(3)]'; input_to_Kv = '[G4(4)]';
sysoutname ='Gdb'; cleanupsysic='yes'; sysic;

    % Now do the same as a function of freq
w = sort([logspace(-4,0,161) 1.1:0.1:10 2.02:0.02:3]);     

% For each configuration do the following
G = sel(Glv,':',[1 2]); Gd= sel(Glv,':',[3 4]);
Glvf = frsp(G,w); Gdlvf=frsp(Gd,w);  % Frequency repsonse of G and Gd
G = sel(Gdv,':',[1 2]); Gd= sel(Gdv,':',[3 4]);
Gdvf = frsp(G,w); Gddvf=frsp(Gd,w);  % Frequency repsonse of G and Gd
G = sel(Gdb,':',[1 2]); Gd= sel(Gdb,':',[3 4]);
Gdbf = frsp(G,w); Gddbf=frsp(Gd,w);  % Frequency repsonse of G and Gd
G = sel(Grr,':',[1 2]); Gd= sel(Grr,':',[3 4]);
Grrf = frsp(G,w); Gdrrf=frsp(Gd,w);  % Frequency repsonse of G and Gd

%-------------------------
% Figure 11: Diagonal RGA-elements
%-------------------------
rlv = sel(vrga(Glvf),1,1);
rdv = sel(vrga(Gdvf),1,1);
rdb = sel(vrga(Gdbf),1,1);
rrr = sel(vrga(Grrf),1,1);

vplot('liv,lm',rlv,rdv,'--',rdb,'--',rrr,'--',1,':');  
axis([0.0001,10,.1,100])
text(.00015,35,'LV','VerticalAlignment','bottom')
text(.00015,17.2,'DB','VerticalAlignment','bottom')
text(.00015,0.44,'DV','VerticalAlignment','top')
text(.00015,3.25,'L/D V/B','VerticalAlignment','top')
xlabel('Frequency [rad/min]')
ylabel('Magnitude of diagonal RGA-element')
%print -deps figrga
%!mv figrga.eps  ~skoge/latex/work/fig/figrga.eps
%!mv figrga.eps ../latex/fig

%--------------------------------
% Figure 14: Open-loop effect of disturbances
%--------------------------------


% Consider effect of feed rate on top composition
dlv = sel(Gdlvf,1,[1 ]);
ddv = sel(Gddvf,1,[1 ]);
ddb = sel(Gddbf,1,[1 ]);
drr = sel(Gdrrf,1,[1 ]);

vplot('liv,lm',dlv,ddv,'-.',ddb,'--',drr,1,':');
axis([0.001,1,.01,100])
text(.002,65,'DB')
text(.002,9,'LV, DV')
text(.004,0.016,'L/D V/B')
xlabel('Frequency [rad/min]')
ylabel('Magnitude')
%print -deps figdf0
%!mv figdf0.eps  ~skoge/latex/work/fig/figdf0.eps
%!mv figdf0.eps  ../latex/fig

% Consider effect of feed composition on top composition
dlv = sel(Gdlvf,1,[2 ]);
ddv = sel(Gddvf,1,[2 ]);
ddb = sel(Gddbf,1,[2 ]);
drr = sel(Gdrrf,1,[2 ]);

vplot('liv,lm',dlv,ddv,'-.',ddb,'--',drr,1,':');
axis([0.001,1,.01,100])
text(.0015,8,'LV, DV, DB, L/D V/B (all configurations)',...
      'VerticalAlignment', 'bottom')
xlabel('Frequency [rad/min]')
ylabel('Magnitude')
%print -deps figdz0
%!mv figdz0.eps  ~skoge/latex/work/fig/figdz0.eps
%!mv figdz0.eps  ../latex/fig

%------------------------------------------------
% Figure 15: Close bottom composition loop (partial control)
%------------------------------------------------
gd1 = sel(Gdlvf,1,[1 2]); gd2 = sel(Gdlvf,2,[1 2]); 
g12 = sel(Glvf,1,2);  g22 = sel(Glvf,2,2);
Pdlv = msub(gd1,mmult(g12,minv(g22),gd2));
gd1 = sel(Gddvf,1,[1 2]); gd2 = sel(Gddvf,2,[1 2]); 
g12 = sel(Gdvf,1,2);  g22 = sel(Gdvf,2,2);
Pddv = msub(gd1,mmult(g12,minv(g22),gd2));
gd1 = sel(Gddbf,1,[1 2]); gd2 = sel(Gddbf,2,[1 2]); 
g12 = sel(Gdbf,1,2);  g22 = sel(Gdbf,2,2);
Pddb = msub(gd1,mmult(g12,minv(g22),gd2));
gd1 = sel(Gdrrf,1,[1 2]); gd2 = sel(Gdrrf,2,[1 2]); 
g12 = sel(Grrf,1,2);  g22 = sel(Grrf,2,2);
Pdrr = msub(gd1,mmult(g12,minv(g22),gd2));

% Feed rate disturbance
plv = sel(Pdlv,1,1);
pdv = sel(Pddv,1,1);
pdb = sel(Pddb,1,1);
prr = sel(Pdrr,1,1);

vplot('liv,lm',plv,pdv,'-.',pdb,'--',prr,1,':');
axis([0.001,1,.01,100])
text(.0015,18,'DV, DB','VerticalAlignment','bottom')
text(.0015,1.3,'LV','VerticalAlignment','bottom')
text(.004,0.016,'L/D V/B','VerticalAlignment','bottom')
xlabel('Frequency [rad/min]')
ylabel('Magnitude')
%print -deps figdf1
%!mv figdf1.eps  ~skoge/latex/work/fig/figdf1.eps
%!mv figdf1.eps  ../latex/fig

% Feed composition disturbance
plv = sel(Pdlv,1,2);
pdv = sel(Pddv,1,2);
pdb = sel(Pddb,1,2);
prr = sel(Pdrr,1,2);

vplot('liv,lm',plv,pdv,'-.',pdb,'--',prr,1,':');
axis([0.001,1,.01,100])
text(.0015,0.0325,'LV','VerticalAlignment','bottom')
text(.0015,19,'DV, DB','VerticalAlignment','bottom')
text(.0015,1.3,'L/D V/B','VerticalAlignment','bottom')
xlabel('Frequency [rad/min]')
ylabel('Magnitude')
%print -deps figdz1
%!mv figdz1.eps  ~skoge/latex/work/fig/figdz1.eps
%!mv figdz1.eps  ../latex/fig

%-------------------------------------------
% Figure 16: Two-point control (CLDG)
%-------------------------------------------

gdiag=vdiag(vdiag(Glvf)); prga = mmult(gdiag,minv(Glvf)); 
cldglv=mmult(prga,Gdlvf);
gdiag=vdiag(vdiag(Gdvf)); prga = mmult(gdiag,minv(Gdvf)); 
cldgdv=mmult(prga,Gddvf);
gdiag=vdiag(vdiag(Gdbf)); prga = mmult(gdiag,minv(Gdbf)); 
cldgdb=mmult(prga,Gddbf);
gdiag=vdiag(vdiag(Grrf)); prga = mmult(gdiag,minv(Grrf)); 
cldgrr=mmult(prga,Gdrrf);

% Feed rate disturbance on xD
clv = sel(cldglv,1,1);
cdv = sel(cldgdv,1,1);
cdb = sel(cldgdb,1,1);
crr = sel(cldgrr,1,1);

vplot('liv,lm',clv,cdv,'-.',cdb,'--',crr,1,':');
axis([0.001,1,.01,100])
text(.0015,45,'LV','VerticalAlignment','bottom')
text(.004,15,'DB','VerticalAlignment','bottom')
text(.0015,8.5,'DV','VerticalAlignment','bottom')
text(.0015,0.0185,'L/D V/B','VerticalAlignment','bottom')
xlabel('Frequency [rad/min]')
ylabel('Magnitude')
%print -deps figdf2
%!mv figdf2.eps  ~skoge/latex/work/fig/figdf2.eps
%!mv figdf2.eps  ../latex/fig

% Feed composition disturbance on xD
clv = sel(cldglv,1,2);
cdv = sel(cldgdv,1,2);
cdb = sel(cldgdb,1,2);
crr = sel(cldgrr,1,2);

vplot('liv,lm',clv,cdv,'-.',cdb,'--',crr,1,':');
axis([0.001,1,.01,100])
text(.0015,1.15,'LV','VerticalAlignment','bottom')
text(.0015,35.5,'DB','VerticalAlignment','bottom')
text(.0015,8.5,'DV','VerticalAlignment','bottom')
text(.0015,4,'L/D V/B','VerticalAlignment','top')
xlabel('Frequency [rad/min]')
ylabel('Magnitude')
%print -deps figdz2
%!mv figdz2.eps  ~skoge/latex/work/fig/figdz2.eps
%!mv figdz2.eps  ../latex/fig/figdz2.eps

% Now evaluate the CLDGs for the BOTTOM:
% Feed rate disturbance on xB
clv = sel(cldglv,2,1);
cdv = sel(cldgdv,2,1);
cdb = sel(cldgdb,2,1);
crr = sel(cldgrr,2,1);

vplot('liv,lm',clv,cdv,'-.',cdb,'--',crr,1,':');
axis([0.001,1,.01,100])
xlabel('Frequency [rad/min]')
ylabel('Magnitude')
text(.0065,43,'LV','VerticalAlignment','bottom')
text(.0065,10.5,'DB','VerticalAlignment','bottom')
text(.0065,1.15,'DV','VerticalAlignment','bottom')
text(.002,0.046,'L/D V/B','VerticalAlignment','bottom')
%print -deps figdf2b
%!mv figdf2b.eps  ~skoge/latex/work/fig/figdf2b.eps
%!mv figdf2b.eps  ../latex/fig/figdf2b.eps

% Feed composition disturbance on xB
clv = sel(cldglv,2,2);
cdv = sel(cldgdv,2,2);
cdb = sel(cldgdb,2,2);
crr = sel(cldgrr,2,2);

vplot('liv,lm',clv,cdv,'-.',cdb,'--',crr,1,':');
axis([0.001,1,.01,100])
text(.0015,9.0,'LV','VerticalAlignment','bottom')
text(.0015,36,'DB','VerticalAlignment','bottom')
text(.0015,0.125,'DV','VerticalAlignment','bottom')
text(.0015,6.2,'L/D V/B','VerticalAlignment','top')
xlabel('Frequency [rad/min]')
ylabel('Magnitude of xB (scaled)')
%print -deps figdz2b
%!mv figdz2b.eps  ~skoge/latex/work/fig/figdz2b.eps
%!mv figdz2b.eps  ../latex/fig/figdz2b.eps


%------------------------------------------------------------
% Figure 17: Two point control LV configuration, with given controllers
%------------------------------------------------------------

k1 = nd2sys([3.76 1],[3.76 1e-8],26.1/100);
k2 = nd2sys([3.31 1],[3.31 1e-8],-37.5/100);

Kdiag = daug(k1,k2);

Kdiagjw = frsp(Kdiag,w);
Llvjw = mmult(Glvf,Kdiagjw);
gdiaglv=vdiag(vdiag(Glvf)); prgalv = mmult(gdiaglv,minv(Glvf)); 
cldglv=mmult(prgalv,Gdlvf);

figure(1)
vplot('liv,lm', sel(Llvjw,1,1),'--',sel(prgalv,1,1),sel(prgalv,1,2),...
                sel(cldglv,1,1),sel(cldglv,1,2),1,':')
xlabel('Frequency [rad/min]')
ylabel('Magnitude')
axis([1e-3 10 1e-2 1e3])
text(1e-2,3e2,'l1','VerticalAlignment','bottom');
text(2.5,3.4,'p11','VerticalAlignment','bottom');
text(3.5,0.28,'p12','VerticalAlignment','top');
text(1e-2,20,'c11','VerticalAlignment','bottom');
text(0.1,0.14,'c12','VerticalAlignment','bottom');
%print -deps fPIfl1
%!mv fPIfl1.eps  ~skoge/latex/work/fig/fPIfl1.eps


figure(2)
vplot('liv,lm', sel(Llvjw,2,2),'--',sel(prgalv,2,1),sel(prgalv,2,2),...
                sel(cldglv,2,1),sel(cldglv,2,2), 1,':')
xlabel('Frequency [rad/min]')
ylabel('Magnitude')
axis([1e-3 10 1e-2 1e3])
text(2e-2,1.5e2,'l2','VerticalAlignment','bottom');
text(0.6e-2,1.5e2,'l2d','VerticalAlignment','bottom');
text(2e-3,62,'c21','VerticalAlignment','bottom');
text(2.5,5.7,'p21','VerticalAlignment','bottom');
text(1.7,0.6,'p22','VerticalAlignment','top');
text(2e-3,9,'c22','VerticalAlignment','bottom');
%print -deps fPIfl2
%!mv fPIfl2.eps  ~skoge/latex/work/fig/fPIfl2.eps


% fjerne figurtekst: close all
% Etterpaa figurer: ps2frag *.eps