Skip to content

Commit 02dbef4

Browse files
dylewskydylewsky
authored
Add files via upload
1 parent 42b3941 commit 02dbef4

File tree

78 files changed

+8537
-0
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

78 files changed

+8537
-0
lines changed

‎CH09/BPOD.m‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
sysBPOD = BPOD(sysFull,sysAdj,r)
2+
3+
[yFull,t,xFull] = impulse(sysFull,0:1:(r*5)+1);
4+
sysAdj = ss(sysFull.A',sysFull.C',sysFull.B',sysFull.D',-1);
5+
[yAdj,t,xAdj] = impulse(sysAdj,0:1:(r*5)+1);
6+
% Not the fastest way to compute, but illustrative
7+
% Both xAdj and xFull are size m x n x 2
8+
HankelOC = []; % Compute Hankel matrix H=OC
9+
for i=2:size(xAdj,1) % Start at 2 to avoid the D matrix
10+
Hrow = [];
11+
for j=2:size(xFull,1)
12+
Ystar = permute(squeeze(xAdj(i,:,:)),[2 1]);
13+
MarkovParameter = Ystar*squeeze(xFull(j,:,:));
14+
Hrow = [Hrow MarkovParameter];
15+
end
16+
HankelOC = [HankelOC; Hrow];
17+
end
18+
[U,Sig,V] = svd(HankelOC);
19+
Xdata = [];
20+
Ydata = [];
21+
for i=2:size(xFull,1) % Start at 2 to avoid the D matrix
22+
Xdata = [Xdata squeeze(xFull(i,:,:))];
23+
Ydata = [Ydata squeeze(xAdj(i,:,:))];
24+
end
25+
Phi = Xdata*V*Sig^(-1/2);
26+
Psi = Ydata*U*Sig^(-1/2);
27+
Ar = Psi(:,1:r)'*sysFull.a*Phi(:,1:r);
28+
Br = Psi(:,1:r)'*sysFull.b;
29+
Cr = sysFull.c*Phi(:,1:r);
30+
Dr = sysFull.d;
31+
sysBPOD = ss(Ar,Br,Cr,Dr,-1);

‎CH09/CH09_SEC02_1_GramianPlot.m‎

Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
clear all, close all, clc
2+
3+
A = [-.75 1; -.3 -.75];
4+
B = [2; 1];
5+
C = [1 2];
6+
D = 0;
7+
8+
sys = ss(A,B,C,D);
9+
10+
Wc = gram(sys,'c'); % Controllability Gramian
11+
Wo = gram(sys,'o'); % Observability Gramian
12+
13+
[sysb,g,Ti,T] = balreal(sys); % Balance the system
14+
15+
BWc = gram(sysb,'c') % Balanced Gramians
16+
BWo = gram(sysb,'o')
17+
18+
%% Plot Gramians
19+
theta = 0:.01:2*pi;
20+
xc = cos(theta);
21+
yc = sin(theta);
22+
CIRC = [xc; yc];
23+
24+
ELLIPb = Ti*sqrt(BWc)*T*CIRC;
25+
ELLIPc = sqrt(Wc)*CIRC;
26+
ELLIPo = sqrt(Wo)*CIRC;
27+
28+
plot(xc,yc,'k--','LineWidth',2)
29+
hold on
30+
31+
% Draw controllability Gramian (unbalanced)
32+
plot(ELLIPc(1,:),ELLIPc(2,:),'r','LineWidth',2)
33+
patch(ELLIPc(1,:),ELLIPc(2,:),'r','FaceAlpha',.75)
34+
35+
% Draw observability Gramian (unbalanced)
36+
plot(ELLIPo(1,:),ELLIPo(2,:),'b','LineWidth',2)
37+
patch(ELLIPo(1,:),ELLIPo(2,:),'b','FaceAlpha',.75)
38+
39+
% Draw balanced Gramians
40+
patch(ELLIPb(1,:),ELLIPb(2,:),'k','FaceColor',[.5 0 .5],'FaceAlpha',.25)
41+
plot(ELLIPb(1,:),ELLIPb(2,:),'Color',[.35 0 .35],'LineWidth',2)
42+
43+
%% Formatting
44+
axis equal, grid on
45+
set(gcf,'Position',[1 1 550 400])
46+
set(gcf,'PaperPositionMode','auto')
47+
% print('-depsc2', '-loose', '../figures/FIG_BT_GRAM');
48+
49+
%% Manually compute scaled balancing transformation
50+
[Tu,D] = eig(Wc*Wo); % Tu are unscaled e-vecs
51+
52+
Atu = inv(Tu)*A*Tu;
53+
Btu = inv(Tu)*B;
54+
Ctu = C*Tu;
55+
Dtu = 0;
56+
syst = ss(Atu,Btu,Ctu,Dtu);
57+
58+
Sigmac = gram(syst,'c') % Diagonal Gramians
59+
Sigmao = gram(syst,'o') % (but not equal)
60+
Sigmas = diag(Sigmac)./diag(Sigmao);
61+
62+
% Scaled balancing transformation
63+
T = Tu*diag(Sigmas.^(1/4));
64+
65+
% Permute columns of T to order Sigma
66+
Sigma = diag(Sigmac).^(1/2).*diag(Sigmao).^(1/2);
67+
[sigsort,permind] = sort(Sigma,'descend');
68+
T = T(:,permind); % Hierarchical
69+
70+
% Compute balanced system
71+
At = inv(T)*A*T;
72+
Bt = inv(T)*B;
73+
Ct = C*T;
74+
Dt = 0;
75+
sysBal = ss(At,Bt,Ct,Dt);
76+
77+
BWc = gram(sysBal,'c') % Balanced Gramians
78+
BWo = gram(sysBal,'o')
79+
80+
81+
82+
%%
83+
% ELLIPb = sqrt(BWc)*CIRC;
84+
85+
% A = [-1 2; 0 -1];
86+
% B = [1; 1];
87+
% C = [1 1];
88+
% D = 0;
89+
90+
91+
% A = [-1 2; -.5 -1];
92+
% B = [1; 1];
93+
% C = [1 1];
94+
% D = 0;

‎CH09/CH09_SEC02_2_BalancedTruncation.m‎

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
clear all, close all, clc
2+
load testSys; % previously saved system
3+
4+
q = 2; % Number of inputs
5+
p = 2; % Number of outputs
6+
n = 100; % State dimension
7+
sysFull = drss(n,p,q); % Discrete random system
8+
r = 10; % Reduced model order
9+
10+
%% Plot Hankel singular values
11+
hsvs = hsvd(sysFull); % Hankel singular values
12+
figure
13+
subplot(1,2,1)
14+
semilogy(hsvs,'k','LineWidth',2)
15+
hold on, grid on
16+
semilogy(r,hsvs(r),'ro','LineWidth',2)
17+
ylim([10^(-15) 100])
18+
subplot(1,2,2)
19+
plot(0:length(hsvs),[0; cumsum(hsvs)/sum(hsvs)],'k','LineWidth',2)
20+
hold on, grid on
21+
plot(r,sum(hsvs(1:r))/sum(hsvs),'ro','LineWidth',2)
22+
set(gcf,'Position',[1 1 550 200])
23+
set(gcf,'PaperPositionMode','auto')
24+
% print('-depsc2', '-loose', '../figures/FIG_BT_HSVS');
25+
26+
%% Exact balanced truncation
27+
sysBT = balred(sysFull,r); % Balanced truncation
28+
29+
%% Compute BPOD
30+
[yFull,t,xFull] = impulse(sysFull,0:1:(r*5)+1);
31+
sysAdj = ss(sysFull.A',sysFull.C',sysFull.B',sysFull.D',-1);
32+
[yAdj,t,xAdj] = impulse(sysAdj,0:1:(r*5)+1);
33+
% Not the fastest way to compute, but illustrative
34+
% Bot xAdj and xFull are size m x n x 2
35+
HankelOC = []; % Compute Hankel matrix H=OC
36+
for i=2:size(xAdj,1) % start at 2 to avoid the D matrix
37+
Hrow = [];
38+
for j=2:size(xFull,1)
39+
Ystar = permute(squeeze(xAdj(i,:,:)),[2 1]);
40+
MarkovParameter = Ystar*squeeze(xFull(j,:,:));
41+
Hrow = [Hrow MarkovParameter];
42+
end
43+
HankelOC = [HankelOC; Hrow];
44+
end
45+
[U,Sig,V] = svd(HankelOC);
46+
Xdata = [];
47+
Ydata = [];
48+
for i=2:size(xFull,1) % start at 2 to avoid the D matrix
49+
Xdata = [Xdata squeeze(xFull(i,:,:))];
50+
Ydata = [Ydata squeeze(xAdj(i,:,:))];
51+
end
52+
Phi = Xdata*V*Sig^(-1/2);
53+
Psi = Ydata*U*Sig^(-1/2);
54+
Ar = Psi(:,1:r)'*sysFull.a*Phi(:,1:r);
55+
Br = Psi(:,1:r)'*sysFull.b;
56+
Cr = sysFull.c*Phi(:,1:r);
57+
Dr = sysFull.d;
58+
sysBPOD = ss(Ar,Br,Cr,Dr,-1);
59+
60+
%% Plot impulse responses for all methods
61+
figure
62+
impulse(sysFull,0:1:60), hold on;
63+
impulse(sysBT,0:1:60)
64+
impulse(sysBPOD,0:1:60)
65+
legend('Full model, n=100','Balanced truncation, r=10','Balanced POD, r=10')
66+
67+
%% Plot impulse responses for all methods
68+
figure
69+
[y1,t1] = impulse(sysFull,0:1:200);
70+
[y2,t2] = impulse(sysBT,0:1:100)
71+
[y5,t5] = impulse(sysBPOD,0:1:100)
72+
subplot(2,2,1)
73+
stairs(y1(:,1,1),'LineWidth',2);
74+
hold on
75+
stairs(y2(:,1,1),'LineWidth',1.2);
76+
stairs(y5(:,1,1),'LineWidth',1.);
77+
ylabel('y_1')
78+
title('u_1')
79+
set(gca,'XLim',[0 60]);
80+
grid on
81+
subplot(2,2,2)
82+
stairs(y1(:,1,2),'LineWidth',2);
83+
hold on
84+
stairs(y2(:,1,2),'LineWidth',1.2);
85+
stairs(y5(:,1,2),'LineWidth',1.);
86+
title('u_2')
87+
set(gca,'XLim',[0 60]);
88+
grid on
89+
subplot(2,2,3)
90+
stairs(y1(:,2,1),'LineWidth',2);
91+
hold on
92+
stairs(y2(:,2,1),'LineWidth',1.2);
93+
stairs(y5(:,2,1),'LineWidth',1.);
94+
xlabel('t')
95+
ylabel('y_2')
96+
set(gca,'XLim',[0 60]);
97+
grid on
98+
subplot(2,2,4)
99+
stairs(y1(:,2,2),'LineWidth',2);
100+
hold on
101+
stairs(y2(:,2,2),'LineWidth',1.2);
102+
stairs(y5(:,2,2),'LineWidth',1.);
103+
xlabel('t')
104+
set(gca,'XLim',[0 60]);
105+
grid on
106+
subplot(2,2,2)
107+
legend('Full model, n=100',['Balanced truncation, r=',num2str(r)],['Balanced POD, r=',num2str(r)])
108+
set(gcf,'Position',[100 100 550 350])
109+
set(gcf,'PaperPositionMode','auto')
110+
% print('-depsc2', '-loose', '../figures/FIG_BT_IMPULSE');

‎CH09/CH09_SEC03_ERA_OKID.m‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
clear all
2+
close all
3+
load testSys
4+
5+
%% Obtain impulse response of full system
6+
[yFull,t] = impulse(sysFull,0:1:(r*5)+1);
7+
YY = permute(yFull,[2 3 1]); % Reorder to be size p x q x m
8+
% (default is m x p x q)
9+
10+
%% Compute ERA from impulse response
11+
mco = floor((length(yFull)-1)/2); % m_c = m_o = (m-1)/2
12+
[Ar,Br,Cr,Dr,HSVs] = ERA(YY,mco,mco,numInputs,numOutputs,r);
13+
sysERA = ss(Ar,Br,Cr,Dr,-1);
14+
15+
%% Compute random input simulation for OKID
16+
uRandom = randn(numInputs,200); % Random forcing input
17+
yRandom = lsim(sysFull,uRandom,1:200)'; % Output
18+
19+
%% Compute OKID and then ERA
20+
H = OKID(yRandom,uRandom,r);
21+
mco = floor((length(H)-1)/2); % m_c = m_o
22+
[Ar,Br,Cr,Dr,HSVs] = ERA(H,mco,mco,numInputs,numOutputs,r);
23+
sysERAOKID = ss(Ar,Br,Cr,Dr,-1);
24+
25+
26+
27+
%% Plot impulse responses for all methods
28+
figure
29+
[y1,t1] = impulse(sysFull,0:1:200);
30+
[y2,t2] = impulse(sysERA,0:1:100);
31+
[y3,t3] = impulse(sysERAOKID,0:1:100);
32+
subplot(2,2,1)
33+
stairs(y1(:,1,1),'LineWidth',2);
34+
hold on
35+
% stairs(y2(:,1,1),'r','LineWidth',1.2);
36+
stairs(y2(:,1,1),'LineWidth',1.2);
37+
stairs(y3(:,1,1),'LineWidth',1.);
38+
% stairs(y5(:,1,1),'m--','LineWidth',1.);
39+
set(gca,'XLim',[0 60]);
40+
grid on
41+
ylabel('y_1')
42+
title('u_1')
43+
subplot(2,2,2)
44+
stairs(y1(:,1,2),'LineWidth',2);
45+
hold on
46+
% stairs(y2(:,1,2),'r','LineWidth',1.2);
47+
stairs(y2(:,1,2),'LineWidth',1.2);
48+
stairs(y3(:,1,2),'LineWidth',1.);
49+
% stairs(y5(:,1,2),'m--','LineWidth',1.);
50+
set(gca,'XLim',[0 60]);
51+
grid on
52+
title('u_2')
53+
subplot(2,2,3)
54+
stairs(y1(:,2,1),'LineWidth',2);
55+
hold on
56+
% stairs(y2(:,2,1),'r','LineWidth',1.2);
57+
stairs(y2(:,2,1),'LineWidth',1.2);
58+
stairs(y3(:,2,1),'LineWidth',1.);
59+
% stairs(y5(:,2,1),'m--','LineWidth',1.);
60+
set(gca,'XLim',[0 60]);
61+
grid on
62+
xlabel('t')
63+
ylabel('y_2')
64+
subplot(2,2,4)
65+
stairs(y1(:,2,2),'LineWidth',2);
66+
hold on
67+
% stairs(y2(:,2,2),'r','LineWidth',1.2);
68+
stairs(y2(:,2,2),'LineWidth',1.2);
69+
stairs(y3(:,2,2),'LineWidth',1.);
70+
% stairs(y5(:,2,2),'m--','LineWidth',1.);
71+
set(gca,'XLim',[0 60]);
72+
grid on
73+
xlabel('t')
74+
legend('Full model, n=100',['ERA, r=',num2str(r)],['ERA/OKID, r=',num2str(r)])
75+
set(gcf,'Position',[100 100 550 350])
76+
set(gcf,'PaperPositionMode','auto')
77+
print('-depsc2', '-loose', '../figures/FIG_ERA_IMPULSE');
78+
79+
%%
80+
%% plot input/output pair for OKID
81+
figure
82+
subplot(1,2,1)
83+
title('Inputs')
84+
stairs(uRandom(1,:))
85+
hold on
86+
stairs(uRandom(2,:))
87+
legend('u_1','u_2')
88+
% xlabel('Time')
89+
% ylabel('Amplitude')
90+
xlabel('t')
91+
ylabel('u')
92+
subplot(1,2,2)
93+
title('Outputs')
94+
stairs(yRandom(1,:))
95+
hold on
96+
stairs(yRandom(2,:))
97+
xlabel('t')
98+
ylabel('y')
99+
legend('y_1','y_2')
100+
% xlabel('Time')
101+
% ylabel('Amplitude')
102+
set(gcf,'Position',[100 100 550 200])
103+
set(gcf,'PaperPositionMode','auto')
104+
print('-depsc2', '-loose', '../figures/FIG_ERA_IODATA');

0 commit comments

Comments
 (0)