% L9_1_maglev.m - MATLAB code for lecture 9.1
%
% RMM, 24 Nov 03
%
% Files needed: none
%

%%
%% Magnetic Levitation System
%%
%%   z1 = position of ball from centerline (down is positive)
%%   z2 = velocity of ball
%%   Vy = output (position sensor) voltage
%%   Vu = input (electromagnet) voltage
%%
%% The dynamics are of the form xdot = f(x, u), y = h(x):
%%
%%        [          z2           ]
%%    f = [ g-km/m*(kA*Vu)^2/z1^2 ]
%%    
%%    h = kT*z1+v0
%%

% Parameter values for system
kT = 613.65;			% gain between position and voltage
v0 = -16.18;			% voltage offset at zero position
m = 0.2;			% mass of ball, kg
g = 9.81;			% gravitational constant
kA = 1;				% electromagnet conductance
km = 3.13e-3 * m/2 / kA^2;	% gain on magnetic attractive force

% Equilibrium point calculation
z10=-v0/kT;
z20=0;
Vu0=1/kA*sqrt(z10^2*m*g/km);

% Compute linearization around the equilibrium
A = [0, 1; 2*km/m*(kA*Vu0)^2/(z10^3), 0];
B = [0; -2*km/m*(kA^2)*Vu0/(z10^2)];
C = [kT 0];
D = 0;

% create a MATLAB LTI system from this representation
magP = ss(A, B, C, D);

%%
%% Maglev controller circuit
%%
%% Note: these values are only approximately equivalent to what exists
%% in the circuit.
%%

k1 = 0.5;				% gain set by gain pot
R1 = 22000;				% Internal resistor
R2 = 22000;				% Resistor plug-in
R = 2000; C = 1e-6;			% RC plug-in

% Controller based on analog circuit
magC1 = -tf([(R1+R)*C 1], [R*C 1])*k1*R2/R1;

magL1 = magP*magC1;			% loop transfer function
magS1 = feedback(1, magP*magC1);	% sensitivity function
magT1 = feedback(magP*magC1, 1);	% closed loop response
figure(1); nyquist(magP, magL1);	% Nyquist plots
figure(2); bode(magP, magC1, magL1);	% Open loop Bode plot

%%
%% Now try to improve the performance by increasing DC gain
%%

magC2 = magC1*5;			% increased gain
magL2 = magP*magC2;			% loop transfer function
magS2 = feedback(1, magP*magC2);	% sensitivity function
magT2 = feedback(magP*magC2, 1);	% closed loop response

magC3 = magC1*20;			% increased gain
magL3 = magP*magC3;			% loop transfer function
magS3 = feedback(1, magP*magC3);	% sensitivity function
magT3 = feedback(magP*magC3, 1);	% closed loop response

figure(3); step(magT1, magT2, magT3);	% Step response
figure(4); bode(magS1, magS2, magS3);	% Closed loop Bode plot
