-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathplot_comparison_results.m
191 lines (159 loc) · 6.31 KB
/
plot_comparison_results.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
%% plot the resilient obserevrs comparison result
% Yu Zheng, 2021, Jan,31
%% Extracting values
time_vec = out.logsout.getElement('x').Values.Time;
% State vectors
x = out.logsout.getElement('x').Values.Data; % true states
error_LO = out.logsout.getElement('error_LO').Values.Data; % error of luenburger
error_L1O = out.logsout.getElement('error_L1O').Values.Data; % error of L1 observer
error_ELO = out.logsout.getElement('error_ELO').Values.Data; % error of event-triggered luenburger
error_MMO = out.logsout.getElement('error_MMO').Values.Data; % error of multi-model observer
error_WL1P = out.logsout.getElement('error_WL1P').Values.Data; % error of pruning observer
% error_H2O = out.logsout.getElement('error_H2O').Values.Data; % error of robust H2 observer
x_hat_LO = out.logsout.getElement('x_hat_LO').Values.Data; % Observed states (Luenberger)
x_hat_L1O = out.logsout.getElement('x_hat_L1O').Values.Data; % Observed states (Unconstrained l1)
x_hat_MMO = out.logsout.getElement('x_hat_MMO').Values.Data; % Observed states (muti-model)
x_hat_ELO = out.logsout.getElement('x_hat_ELO').Values.Data; % Observed states (event-triggered luenburger)
x_hat_WL1P = out.logsout.getElement('x_hat_WL1P').Values.Data; % Observed states (pruning)
% x_hat_H2O = out.logsout.getElement('x_hat_H2O').Values.Data; % Observed states (robust H2 observer)
% FDIA performance
aux_model_res = out.logsout.getElement('y_obsv').Values.Data; % measured y for observer (y_obsv = y-D*u)
BDD_res = out.logsout.getElement('BDD_res').Values.Data; % Bad data residue
% prior performance
q = out.logsout.getElement('q').Values.Data;
q_hat = out.logsout.getElement('q_hat').Values.Data;
%% Ploting/Visualization/Metric tables
% figure(1) % errors
figure(2) % states
% figure(3) % BDD
LW = 1.3; % linewidth
FS = 15; % font size
% % error visulization
% figure(1)
% subplot(5,1,1)
% plot(time_vec,error_LO,'LineWidth',LW)
% ylabel('error','FontWeight','bold','Fontsize',FS)
% title('Luenburger Observer','Fontsize',FS)
% set(gca,'fontweight','bold','fontsize',FS)
% set(gca,'LineWidth',1)
%
% figure(1)
% subplot(5,1,2)
% plot(time_vec,error_L1O,'LineWidth',LW)
% ylabel('error','FontWeight','bold','Fontsize',FS)
% title('Unconstrained L_1 Observer','Fontsize',FS)
% set(gca,'fontweight','bold','fontsize',FS)
% set(gca,'LineWidth',1)
%
%
% figure(1)
% subplot(5,1,3)
% plot(time_vec,error_MMO,'LineWidth',LW)
% ylabel('error','FontWeight','bold','Fontsize',FS)
% title('Multi-model Observer','Fontsize',FS)
% set(gca,'fontweight','bold','fontsize',FS)
% set(gca,'LineWidth',1)
%
%
% figure(1)
% subplot(5,1,4)
% plot(time_vec,error_ELO,'LineWidth',LW)
% ylabel('error','FontWeight','bold','Fontsize',FS)
% title('Event-triggered Luenburger Observer','Fontsize',FS)
% set(gca,'fontweight','bold','fontsize',FS)
% set(gca,'LineWidth',1)
%
%
% figure(1)
% subplot(5,1,5)
% plot(time_vec,error_WL1P,'LineWidth',LW)
% ylabel('error','FontWeight','bold','Fontsize',FS)
% title('Resilient Pruning Observer','Fontsize',FS)
% set(gca,'fontweight','bold','fontsize',FS)
% set(gca,'LineWidth',1)
% states estimation compare
n_comp = 5; % number of observes in comparison
n_pres = 5; % number of states presented
figure (2)
for iter=1:n_pres
subplot(n_pres,n_comp,1+(iter-1)*n_comp)
plot(time_vec,error_LO(:,iter),'k','LineWidth',LW)
if(iter == 1)
title('LO','Fontsize',FS)
end
ylabel(['$$\delta_{' num2str(iter) '}-\hat{\delta_{' num2str(iter) '}}$$'],'FontWeight','bold','Fontsize',FS,'interpreter','latex')
set(gca,'fontweight','bold','fontsize',12)
set(gca,'LineWidth',1)
axis([0 8 -50e-4 50e-4])
subplot(n_pres,n_comp,2+(iter-1)*n_comp)
plot(time_vec,error_L1O(:,iter),'k','LineWidth',LW)
if(iter == 1)
title('UL1O','Fontsize',FS)
end
set(gca,'fontweight','bold','fontsize',12)
set(gca,'LineWidth',1)
axis([0 8 -50e-4 50e-4])
subplot(n_pres,n_comp,3+(iter-1)*n_comp)
plot(time_vec,error_ELO(:,iter),'k','LineWidth',LW)
if(iter == 1)
title('ETLO','Fontsize',FS)
end
set(gca,'fontweight','bold','fontsize',12)
set(gca,'LineWidth',1)
axis([0 8 -50e-4 50e-4])
subplot(n_pres,n_comp,4+(iter-1)*n_comp)
plot(time_vec,error_MMO(:,iter),'k','LineWidth',LW)
if(iter == 1)
title('MMO','Fontsize',FS)
end
set(gca,'fontweight','bold','fontsize',12)
set(gca,'LineWidth',1)
axis([0 8 -50e-4 50e-4])
subplot(n_pres,n_comp,5+(iter-1)*n_comp)
plot(time_vec,error_WL1P(:,iter),'k','LineWidth',LW)
if(iter == 1)
title('RPO','Fontsize',FS)
end
set(gca,'fontweight','bold','fontsize',12)
set(gca,'LineWidth',1)
axis([0 8 -50e-4 50e-4])
end
% BDD_res
figure(3)
plot(time_vec,BDD_res,'k','LineWidth',LW), hold on
plot(time_vec,BDD_thresh*ones(length(time_vec),1),'r--','LineWidth',2*LW)
ylabel('Bad Data Detection Residual','FontWeight','bold')
% error metric table (for rotor angle)
n_delta = size(x,2)/2; % number of generator angles
error_LO_delta = error_LO(:,1:n_delta).';
error_L1O_delta = error_L1O(:,1:n_delta).';
error_ELO_delta = error_ELO(:,1:n_delta).';
error_MMO_delta = error_MMO(:,1:n_delta).';
error_WL1P_delta = error_WL1P(:,1:n_delta).';
metric_rms = @(x,dim) sqrt(sum(x.^2,dim)/size(x,dim));
metric_max = @(x,dim) max(abs(x),[],dim);
error_table = [metric_rms(error_LO_delta,2) metric_rms(error_L1O_delta,2) metric_rms(error_ELO_delta,2) metric_rms(error_MMO_delta,2) metric_rms(error_WL1P_delta,2) ...
metric_max(error_LO_delta,2) metric_max(error_L1O_delta,2) metric_max(error_ELO_delta,2) metric_max(error_MMO_delta,2) metric_max(error_WL1P_delta,2)];
%% prior evaluation
tot = size(time_vec,1);
PPV = zeros(tot,1);
for j = 1:tot
safe_hat = find(q_hat(j,:));
safe = find(q(j,:));
TP = numel(intersect(safe_hat, safe));
risk = find(~q(j,:));
FP = numel(intersect(safe_hat, risk));
PPV(j) = TP/(TP+FP);
end
% plot
figure(4)
plot(time_vec,PPV,'k','LineWidth',LW)
ylabel('PPV','FontWeight','bold','Fontsize',12)
xlabel('time','FontWeight','bold','Fontsize',12)
% title('Precision of prior','Fontsize',FS)
set(gca,'fontweight','bold','fontsize',12)
ylabel('PPV','FontWeight','bold','Fontsize',FS)
xlabel('time','FontWeight','bold','Fontsize',FS)
title('Precision of prior','Fontsize',FS)
set(gca,'fontweight','bold','fontsize',FS)
set(gca,'LineWidth',1)