
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab script for stability analysis of feedback systems
%
% Vishal Saxena
% University of Idaho
% http://lumerink.com
% Copyright 2016
% PS: Verify the accuracy of results yourself when 
% this script using for HW/Projects 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; clc; close all;

%% Define Opamp parameters
A0 = 100e3; % Gain factor of the transfer function

% Pole locations, use positive values for LHP poles
wp1 = 0;
wp2 = 0;
wp3 = [];
wp =  [wp1 wp2 wp3];

% Zero locations, negative value indicates an LHP zero
%wz = [];
wz = -[10];

% Feedback factor (1/k)
k = 10;     % ideal closed loop gain

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read documentation for the zpk function before using this script.
% This function likes the values in the form:
%     A0(s-z1)/((s+p1)(s+p2)(s+p3))
% Modify A0 to make sure you plot the desired response
%  -> factor of prod(wp)/prod(wz) for non-zero poles to fix the response
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

WP = wp(find(wp)); % Product of non-zero Pole values
WZ = abs(wz(find(wz))); % Product of non-zero Zero values
Gain = A0*prod(WP)/prod(WZ);



%% Evaluate transfer functions:

% Forward Path response
G = zpk(wz, -wp, Gain)*(1/k);

% Loop response
L = G*(1/k)


% Unity-gain closed loop response
H_un = minreal(G/(1+G));
% Minreal cancels the common terms in the Nr and Dr

% Closed loop response
H_cl = minreal(k*L/(1+L)); % Minimum realization


% Screen Housekeeping
scrsz = get(0,'ScreenSize'); 
Xscr = scrsz(3); Yscr = scrsz(4); Width= 0.32*Xscr; Height= 0.25*Yscr;
% Usage: figure('Position',[left, bottom, width, height])


%% Plot Open Loop Frequency Response
figure('Position',[1 2*Yscr/3 Width Height])

Wmax = max(1e3, 10*max(wp)); % maximum freq for the plot

bode(L, {0.1, Wmax});
[Gm,Pm,Wg,Wp] = margin(L);
sprintf('Phase Margin = %g', Pm)
sprintf('Wun = %d', Wp)

hold on;
bode(G, {0.1, Wmax});
legend('L','G');
h = findobj(gcf,'type','line');
set(h,'linewidth',2);
grid on;
title('Open Loop Frequency Response, L(s)');

%% Plot Unity Gain Frequency Response
figure('Position',[1 Yscr/3 Width Height])
bode(H_un, 'r', {0.1, Wmax});
grid on;
h = findobj(gcf,'type','line');
set(h,'linewidth',2);
title('Unity Gain Response: G/(1+G)');

%% Plot Closed Loop Frequency Response
figure('Position',[1 1  Width Height])

bode(H_cl, 'r', {0.1, Wmax});
grid on;
h = findobj(gcf,'type','line');
set(h,'linewidth',2);
title('Closed Loop Frequency Response, k*L/(1+L)');

%% Plot Unity Gain Step Response
figure('Position',[Xscr/3 Yscr/3  Width Height])

step(H_un);
grid on;
h = findobj(gcf,'type','line');
set(h,'linewidth',2);
title('Unity Gain Step Response');

%% Plot Closed Loop Step Response
figure('Position',[Xscr/3 1  Width Height])

step(H_cl);
grid on;
h = findobj(gcf,'type','line');
set(h,'linewidth',2);
title('Closed Loop Step Response');


%% Plot Root Locus
figure('Position',[Xscr/3 2*Yscr/3 Width Height])
rlocus(G);
title('Root Locus');


%% Nyquist Plot
figure('Position',[2*Xscr/3 2*Yscr/3  Width Height])
nyquist(G);
grid on;
h = findobj(gcf,'type','line');
set(h,'linewidth',2);
title('Nyquist Stability Plot');

%% EOF