-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathRST_ROBUST.m
192 lines (162 loc) · 6.82 KB
/
RST_ROBUST.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
192
%% INÍCIO DA ROTINA
%% LUÍS AUGUSTO MESQUITA DE CASTRO (01/04/2017)
% Instituto Federal do Pará (IFPA)
% Universidade Federal do Pará (UFPA)
% Controle Digital de Sistemas (Mestrado em Engenharia Elétrica - UFPA)
% Controle Preditivo e Estocástico (Mestrado em Engenharia Elétrica - UFPA)
%% DaqDuino Data Acquisition device.
% DAQ-Duino, 2013-2016
% Author: Prof. Dr. Antonio Silveira (asilveira@ufpa.br)
% Laboratory of Control and Systems (LACOS), UFPA (www.ufpa.br)
%% Limpar todas as variáveis do workspace
% clear; close all; clc
%% Obter realização em função de transferência discreto do modelo identificado
disp('ANÁLISE DE ROBUSTEZ DA MALHA DE CONTROLE');
%% Análise de Sensibilidade
% Az = [1 a1 a2]; Bz = [b0 b1]; Rz = [r0 r1]; Sz = [s0 s1 s2 s3]; Tz = [t0 t1 t2]; Dz = [1 -1];
Tmf = tf(conv(Bz,Tz),conv(conv(Dz,Az),Rz)+conv(Bz,[Sz 0]),Ts); % Função de Sensibilidade Complementar
Si = tf(conv(Bz,conv(Dz,Rz)),conv(conv(Dz,Az),Rz)+conv(Bz,[Sz 0]),Ts); % Função de Sensibilidade de entrada
So = tf(conv(conv(Dz,Az),Rz),conv(conv(Dz,Az),Rz)+conv(Bz,[Sz 0]),Ts); % Função de Sensibilidade de saída
% % Priorizando R(z^-1)
% Az = [1 a1 a2]; Bz = [b0 b1]; Rz = [r0r r1r]; Sz = [s0r s1r s2r s3r]; Tz = [t0 t1 t2]; Dz = [1 -1];
% Tmfr = tf(conv(Bz,Tz),conv(conv(Dz,Az),Rz)+conv(Bz,Sz),Ts); % Função de Sensibilidade Complementar
% Sir = tf(conv(Bz,conv(Dz,Rz)),conv(conv(Dz,Az),Rz)+conv(Bz,Sz),Ts); % Função de Sensibilidade de entrada
% Sor = tf(conv(conv(Dz,Az),Rz),conv(conv(Dz,Az),Rz)+conv(Bz,Sz),Ts); % Função de Sensibilidade de saída
% % Priorizando S(z^-1)
% Az = [1 a1 a2]; Bz = [b0 b1]; Rz = [r0s r1s]; Sz = [s0s s1s s2s s3s]; Tz = [t0 t1 t2]; Dz = [1 -1];
% Tmfs = tf(conv(Bz,Tz),conv(conv(Dz,Az),Rz)+conv(Bz,Sz),Ts); % Função de Sensibilidade Complementar
% Sis = tf(conv(Bz,conv(Dz,Rz)),conv(conv(Dz,Az),Rz)+conv(Bz,Sz),Ts); % Função de Sensibilidade de entrada
% Sos = tf(conv(conv(Dz,Az),Rz),conv(conv(Dz,Az),Rz)+conv(Bz,Sz),Ts); % Função de Sensibilidade de saída
%% Diagrama de Bode do sistema em malha aberta e em malha fechada
% Resposta em frequência do sistema em malha aberta
[a,b,c] = bode(Gz);
%% Inicializar vetores
w1 = zeros(1,length(a)); % Vetor de frequências
pha1 = zeros(1,length(a)); % Vetor de margem de fase
mag1 = zeros(1,length(a)); % Vetor de margem de ganho
for k = 1:length(w1)
w1(1,k) = c(k,1);
pha1(1,k) = b(1,1,k);
mag1(1,k) = a(1,1,k);
end
magdB1 = mag2db(mag1); % Converter valores para decibels (dB)
% Resposta em frequência do sistema em malha fechada
[a,b,c] = bode(Tmf,w1);
%% Inicializar vetores
w2 = zeros(1,length(a)); % Vetor de frequências
pha2 = zeros(1,length(a)); % Vetor de margem de fase
mag2 = zeros(1,length(a)); % Vetor de margem de ganho
for k = 1:length(w2)
w2(1,k) = c(k,1);
pha2(1,k) = b(1,1,k);
mag2(1,k) = a(1,1,k);
end
magdB2 = mag2db(mag2); % Converter valores para decibels (dB)
%% Resultados
figure(1); % Figura 1
subplot(211);
semilogx(w1,magdB1,'b','linewidth',2); hold on; grid on
semilogx(w2,magdB2,'r','linewidth',2); hold on; grid on
set(gca,'fontsize',14);
set(gca,'linewidth',1);
set(gca,'xscale','log');
ylim([min([magdB1 magdB2]) max([magdB1 magdB2])+5]);
title('Bode diagram: open loop and closed loop system');
ylabel('magnitude (dB)');
legend('Open loop','Closed loop','Location','southeast');
subplot(212);
semilogx(w1,pha1,'b','linewidth',2); hold on; grid on
semilogx(w2,pha2,'r','linewidth',2); hold on; grid on
set(gca,'fontsize',14);
set(gca,'linewidth',1);
set(gca,'xscale','log');
ylim([min([pha1 pha2]) max([pha1 pha2])+50]);
xlabel('frequency (rad/s)'); ylabel('phase (degrees)');
legend('Open loop','Closed loop','Location','southeast');
%% Funções de sensibilidade
disp('FUNÇÕES DE SENSIBILIDADE:');
Mt = norm(Tmf,Inf); % Norma infinita (valor absoluto - taxa de amplificação)
Msi = norm(Si,Inf); % Norma infinita (valor absoluto - taxa de amplificação)
Mso = norm(So,Inf); % Norma infinita (valor absoluto - taxa de amplificação)
disp('Valor máximo de T(z):');
disp(Mt);
disp('Valor máximo de Si(z):');
disp(Msi);
disp('Valor máximo de So(z):');
disp(Mso);
% Margem de ganho
MGT = 1+(1/Mt); % Valor absoluto
MGTdB = mag2db(MGT); % Valor em dB
MGSi = (Msi/(Msi-1)); % Valor absoluto
MGSidB = mag2db(MGSi); % Valor em dB
MGSo = (Mso/(Mso-1)); % Valor absoluto
MGSodB = mag2db(MGSo); % Valor em dB
disp('Margem de ganho da função de sensibilidade complementar T(z) em dB:');
display(MGTdB);
disp('Margem de ganho da função de sensibilidade de entrada Si(z) em dB:');
display(MGSidB);
disp('Margem de ganho da função de sensibilidade de saída So(z) em dB:');
display(MGSodB);
% Margem de fase
MFT = 2*asin(1/(2*Mt))*(180/pi);
MFSi = 2*asin(1/(2*Msi))*(180/pi);
MFSo = 2*asin(1/(2*Mso))*(180/pi);
disp('Margem de fase da função de sensibilidade complementar T(z) em graus:');
display(MFT);
disp('Margem de fase da função de sensibilidade de entrada Si(z) em graus:');
display(MFSi);
disp('Margem de fase da função de sensibilidade de saída So(z) em graus:');
display(MFSo);
%% Diagrama de Bode de Si(z)
% Resposta em frequência da função de sensibilidade de entrada
[a,b,c] = bode(Si,w1);
%% Inicializar vetores
w1 = zeros(1,length(a)); % Vetor de frequências
pha1 = zeros(1,length(a)); % Vetor de margem de fase
mag1 = zeros(1,length(a)); % Vetor de margem de ganho
for k = 1:length(w1)
w1(1,k) = c(k,1);
pha1(1,k) = b(1,1,k);
mag1(1,k) = a(1,1,k);
end
magdB1 = mag2db(mag1); % Converter valores para decibels (dB)
%% Resultados
figure(2); % Figura 2
% subplot(211);
semilogx(w1,magdB1,'r','linewidth',2); hold on; grid on
set(gca,'fontsize',14);
set(gca,'linewidth',1);
set(gca,'xscale','log');
ylim([min(magdB1) max(magdB1)+5]);
title('Input sensitivity function S_i(z^{-1})');
ylabel('singular values (dB)');
xlabel('frequency (rad/s)');
legend('S_i(z^{-1})','Location','northwest');
%% Diagrama de Bode de So(z)
% Resposta em frequência da função de sensibilidade de saída
[a,b,c] = bode(So,w1);
%% Inicializar vetores
w1 = zeros(1,length(a)); % Vetor de frequências
pha1 = zeros(1,length(a)); % Vetor de margem de fase
mag1 = zeros(1,length(a)); % Vetor de margem de ganho
for k = 1:length(w1)
w1(1,k) = c(k,1);
pha1(1,k) = b(1,1,k);
mag1(1,k) = a(1,1,k);
end
magdB1 = mag2db(mag1); % Converter valores para decibels (dB)
%% Resultados
figure(3); % Figura 3
% subplot(211);
semilogx(w1,magdB1,'r','linewidth',2); hold on; grid on
set(gca,'fontsize',14);
set(gca,'linewidth',1);
set(gca,'xscale','log');
ylim([min(magdB1) max(magdB1)+5]);
title('Output sensitivity function S_o(z^{-1})');
ylabel('singular values (dB)');
xlabel('frequency (rad/s)');
legend('S_o(z^{-1})','Location','northwest');
disp('FIM DA ANÁLISE DE ROBUSTEZ DA MALHA DE CONTROLE');
%% FIM DA ROTINA
%% LUÍS AUGUSTO MESQUITA DE CASTRO (01/04/2017)