% L18_1dfan.m - ducted fan example for L18.1
% RMM, 24 Feb 03

%% Ducted fan dynamics 
J = 0.0475;				% inertia around pitch axis
m = 1.5;				% mass of fan
r = 0.25;				% distance to flaps
d = 0.2;				% damping factor (estimated)
g = 10;					% gravitational constant
l = 0.05;				% offset of center of mass
dfan = tf([r], [J, d, m*g*l]);

%% Lead compensator
a = 25; b = 300; k = 15;		% controller parameters
ctrl = tf([k, k*a], [1/b, 1]);

%% Utility transfer functions
P = dfan; C = ctrl;
S = 1/(1+P*C);				% sensitivity function
T = P*C/(1+P*C);			% complementary senstivity function

%% Frequency vectors for evaluating performance
Wv = logspace(0, 3);
Sv = reshape(freqresp(S, Wv), 1, length(Wv));
Tv = reshape(freqresp(T, Wv), 1, length(Wv));

%%
%% Nominal performance
%%
wp = 12;				% performance weight parameters
W1 = tf([20], [1/wp^2, 2/wp, 1]);
W1v = reshape(freqresp(W1, Wv), 1, length(Wv));
bode(W1*S);				% check for nominal performance

%%
%% Perturbation #1: 20% uncertainty in r
%%
W2 = tf([0.2], [1]);			% multiplicative perturbation
W2v = reshape(freqresp(W2, Wv), 1, length(Wv));
bode(W2*T);				% check for robust stability

RPv = abs(W1v.*Sv) + abs(W2v.*Tv);	% robust performance magnitude
loglog(Wv, RPv, Wv, ones(length(Wv)));	% plot to see if it is OK

%%
%% Perturbation #2: 50% uncertainty in l
%%
W2 = tf([0.25*m*g/r], [1]);		% feedback perturbatoin
W2v = reshape(freqresp(W2, Wv), 1, length(Wv));
bode(W2*P*S);				% check for robust stability

%%
%% Perturbation #3: motor dynamics (pure time delay)
%%
[num, den] = pade(0.1, 2);		% second order pade approximation
W2 = tf(num, den) - 1;			% multiplicative perturbation
W2v = reshape(freqresp(W2, Wv), 1, length(Wv));

bode(W2*T);				% check for robust stability
RSv = abs(W2v .* Tv);			% more careful check
loglog(Wv, RSv, Wv, ones(length(Wv)));	

RPv = abs(W1v.*Sv) + abs(W2v.*Tv);	% robust performance magnitude
loglog(Wv, RPv, Wv, ones(length(Wv)));	% plot to see if it is OK