-
Notifications
You must be signed in to change notification settings - Fork 4
/
rhpz_sensitivity.m
125 lines (106 loc) · 3.42 KB
/
rhpz_sensitivity.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
% Sensitivity Limits Imposed by RHP Zero:
% Sensitivity bounds based on RHPZ and Weighting function
% SISO LTI plant with P-K classic feedback structure
% Stable plant
% -----------------------------------------------------------
clear;
s = tf('s');
% Available Bandwidth
wp = 1e5;
% RHP zero
z = 10;
% % % Generic Weighting FUnction:
% W = ((s+ws*nthroot(M,k2))^k2 * (s+wp/nthroot(M,k3))^k3)/ ...
% ((s+ws*nthroot(epsilon,k1))^k1 * (s+ws)^(k2-k1) * (s+wp)^k3);
% Order of transfer function in a given frequency range
% First slope:
k1 = 1;
% Second slope:
k2 = 1;
% Third slope:
k3 = 1;
% Performance Bandwidth Vector
ws_vec = logspace(-1,1,1000);
% Peak Sensitivity Value Vector
M_vec = logspace(0,2,1000);
% Epsilon (magnitude of sensitivity at low freq.)
% Assume that epsilon is 0 (or 0+ to be precise)
% --------------------------------------------------
% % Solve for M for different value of ws
% Expression involving M and ws
% z^k1 * (z+ws)^(k2-k1) * (z+wp)^k3 ...
% - (z+ws*nthroot(M,k2))^k2 * (z+wp/nthroot(M,k3))^k3 == 0;
% Preallocate
Msoln_vec = NaN*ones(numel(ws_vec),1);
validM_vec = NaN*ones(numel(ws_vec),1);
for ws_ind = 1:numel(ws_vec)
ws = ws_vec(ws_ind);
M_param_func = @(M,ws,wp,z,k1,k2,k3) ...
(z+ws*nthroot(M,k2))^k2 * (z+wp/nthroot(M,k3))^k3 ...
- z^k1 * (z+ws)^(k2-k1) * (z+wp)^k3;
% "Single" variable function:
M_singlevar_fun = @(M) M_param_func(M,ws,wp,z,k1,k2,k3);
% Initial Point:
M_init = 1;
% Solve for M (root of nonlinear expression):
try
M_soln = fzero(M_singlevar_fun,M_init);
catch
M_soln = NaN; %
end
% Store the solution:
Msoln_vec(ws_ind) = M_soln;
% For the assumed upper bound on sensitivity, the relations are valid
% under certain assumption. This assumption is based on relation
% between ws and wp.
validM_vec(ws_ind) = (wp/ws)^((k2*k3)/(k2+k3));
end
% % Plot
figure;
semilogx(ws_vec,mag2db(Msoln_vec),'-b');
grid on; hold on;
semilogx(ws_vec,mag2db(validM_vec),'-r');
title('LB on M vs \omega_s');
ylabel('LB on M (dB)');
xlabel('\omega_s (rad/s)');
plot_axis;
ylim([-10 60])
% -------------------------------------------------------------
% Solve for ws for different value of M
% Expression involving ws:
% ws <= (pi*p + k3/nthroot(M,k3)*wp - k3*wp)/(-k1 - k2*nthroot(M,k2) + k2)
% Preallocate
ws_vec = NaN*ones(numel(ws_vec),1);
validws_vec = NaN*ones(numel(ws_vec),1);
for M_ind = 1:numel(M_vec)
M = M_vec(M_ind);
ws_param_func = @(ws,M,wp,z,k1,k2,k3) ...
(z+ws*nthroot(M,k2))^k2 * (z+wp/nthroot(M,k3))^k3 ...
- z^k1 * (z+ws)^(k2-k1) * (z+wp)^k3;
% "Single" variable function:
ws_singlevar_fun = @(ws) ws_param_func(ws,M,wp,z,k1,k2,k3);
% Initial Point:
ws_init = 1;
% Solve for ws (root of nonlinear expression):
try
ws_soln = fzero(ws_singlevar_fun,ws_init);
catch
ws_soln = NaN; %
end
% Store ws
ws_vec(M_ind) = ws_soln;
% For the assumed upper bound on sensitivity, the relations are valid
% under certain assumption. This assumption is based on relation
% between ws and wp.
validws_vec(M_ind) = wp/(nthroot(M,k2)*nthroot(M,k3));
end
% % Plot
figure;
semilogy(mag2db(M_vec),ws_vec,'-b');
grid on; hold on;
semilogy(mag2db(M_vec),validws_vec,'-r');
title('UP on \omega_s vs M');
xlabel('M (dB)');
ylabel('UB on Perf. Bandwidth \omega_s (rad/s)');
plot_axis;
ylim([1e-2 1e1])