Project 12  ·  AME 451

Linear Control Systems — MATLAB Portfolio

AME 451 · Linear Control Systems · University of Southern California

MATLAB Root Locus Bode Plots Nyquist Criterion Routh-Hurwitz Lead-Lag Compensation PID Tuning Nichols Chart Frequency Response Transfer Functions

Course

AME 451 — Linear Control Systems

Software

MATLAB / Symbolic Math Toolbox

Topics

Root Locus, Bode, Nyquist, Lead-Lag

Scripts

40+ .m files & .mlx live scripts

Key Tool

Custom Routh-Hurwitz function

Contents

  1. Routh-Hurwitz Stability Criterion
  2. Root Locus Analysis & Design
  3. Bode Plots & Frequency Response
  4. Nyquist, Nichols & Multi-Domain Stability
  5. Lead-Lag Compensator Design
  6. Step Response & Gain Tuning
Stability Criterion

Routh-Hurwitz Stability Analysis

A custom routh_hurwitz() function was written to automate Routh array construction for both numerical and symbolic polynomials. It handles both special cases — a zero in the first column (replaced with epsilon, then the limit taken as ε→0⁺) and a row of zeros (replaced using the auxiliary polynomial method). The function returns the full Routh array M and the simplified first column L; sign changes in L equal the number of roots in the right-half plane.

routh_hurwitz.m  ·  Custom Stability Function 
function [M, L] = routh_hurwitz(P, N)
% Returns Routh array M and simplified first column L.
% Handles special cases: zero in first column (epsilon method),
% and row of zeros (auxiliary polynomial method).
% P — numerical or symbolic coefficient array [a_n, a_{n-1}, ..., a_0]
% N — precision digits (default 10); values < 10^-N treated as zero.

if nargin < 2, N = 10; end
n = length(P);
P = sym(P);
syms epsilon

nfil = n;  ncol = ceil(n/2);
M = sym(zeros(nfil, ncol));

% Fill first two rows from coefficient array
k = 1;
for j = 1:ncol
    for i = 1:2
        M(i,j) = P(k);  k = k+1;
    end
end

% Build remaining rows via 2x2 determinants
for i = 3:nfil
    % Special case: first element of row is zero
    if abs(double(M(i-1,1))) < 10^(-N)
        M(i-1,1) = epsilon;                          % epsilon replacement
    % Special case: entire row is zero
    elseif sum(abs(double(M(i-1,:)))) < 10^(-N)
        for j = 1:ncol                               % use auxiliary polynomial
            M(i-1,j) = M(i-2,j) * (n-(i-2)-2*(j-1));
        end
    end
    for j = 1:ncol-1
        M(i,j) = det([M(i-2,1) M(i-2,j+1); ...
                      M(i-1,1) M(i-1,j+1)]) / (-M(i-1,1));
    end
end

M = vpa(M);
M1 = limit(M, epsilon, 0, 'right');     % resolve epsilon → 0+
L  = vpa(zeros(n,1));
for i = 1:n
    L(i,1) = M1(i,1);                  % sign changes in L = # RHP roots
end

Used across homework and final exam to determine stability margins symbolically — solving for the critical gain K at which the first column element crosses zero:

FinalQ3.m  ·  Symbolic K margin via Routh array 
syms s K
den_ol = expand((s+1)*(3*s+1)*(0.4*s+1));   % open-loop denominator
num_ol = sym(9);                              % gain k = 9

% Characteristic equation: 1 + K·G(s) = 0
CharEq     = simplify(den_ol + K*num_ol);
CharCoeff  = fliplr(coeffs(CharEq, s));      % [a_n ... a_0]
RouthArray = simplify(routh_hurwitz(CharCoeff));

% Solve for K where first-column element → 0 (stability boundary)
K_critical = solve(RouthArray(3,1))          % K ≈ 15.86
Key Result
Characteristic equation: 1.2s³ + 4.6s² + 4.4s + 1 + K·0.4·... = 0. Routh stability boundary: K_critical ≈ 15.86 — the maximum proportional gain before the closed-loop system becomes unstable. For the operating gain K = 9 all first-column elements are positive → system is stable with positive phase and gain margins.
Root Locus

Root Locus Analysis & Compensator Design

Root locus techniques were applied to visualize how closed-loop poles migrate as gain K is varied. Performance specifications (max overshoot Mp, natural frequency ωn) were overlaid as constraint curves using sgrid(), enabling graphical compensator placement. The open-loop transfer function for the final exam root locus question included a lead compensator k = 4.35(s+1) designed to pull loci into the acceptable damping/frequency region.

FinalQ1.m  ·  Root locus with ζ = 0.4, ωn = 2 constraints 
clear; clc; close all
syms s k

% Lead compensator + plant
k      = 4.35*(s+1);
num_ol = sym(40*s + 40);                         % k·(s+1) factor
den_ol = expand((s+10) * (10*s^2) * (0.08*s+1));

g_ol = tf(sym2poly(num_ol), sym2poly(den_ol));

% Root locus with performance constraints overlaid
rlocus(g_ol)
xlim([-3, 1]);  ylim([-5, 5])
xline(-1, ':')                                   % Re(s) = -1 boundary
yline([1, -1], ':')                              % Im(s) = ±1 lines
sgrid(0.4, 2)                                    % ζ = 0.4, ωn = 2 rad/s

Breakaway and break-in points were found symbolically using the backward asymptote equation d/ds[−P(s)/Z(s)] = 0, and Routh-Hurwitz was applied to find exact cross-over gain:

FinalQ2A.m  ·  Breakaway points & stability margin 
clear; clc; close all
syms s k

num = expand((s+5)*2);          % G(s) = 2(s+5) / [(s+1)(s+3)]
den = expand((s+1)*(s+3));

% Breakaway equation: d/ds[-den/num] = 0
bkwyEq = simplify(diff(-den/num));
bkwyPt = double(solve(bkwyEq))  % real breakaway point(s)

% Characteristic equation for Routh analysis
CharEq     = simplify(den + k*num);
CharCoeff  = fliplr(coeffs(CharEq, s));
RouthArray = simplify(routh_hurwitz(CharCoeff));

g = tf(sym2poly(num), sym2poly(den));
rlocus(g);  grid
Result
Plant G(s) = 2(s+5)/[(s+1)(s+3)] has zeros at s = −5 and poles at s = −1, −3. Breakaway point computed symbolically. Root locus exits the real axis and the system remains stable for all finite K > 0 (minimum-phase, more zeros than poles in RHP → none).
Frequency Response

Bode Plots & Asymptotic Approximation

Frequency response was analyzed by plotting the exact Bode magnitude and phase alongside manually-constructed asymptotic (straight-line) approximations using a custom asymp() helper. A separate BodePaper() utility generates blank logarithmic graph paper as a MATLAB figure — useful for overlaying exact vs. approximate curves.

HW8Q1A.m  ·  Exact Bode vs. asymptotic approximation 
clear; clc; close all

% G(s) = 12 / [(s+4)(s+1/3)]
%   Two real poles at s = -4, -1/3
%   Low-frequency gain = 12 / (4 × 1/3) = 9  ≈ +19.1 dB
num = [12];
den = conv([1 4], [1 1/3]);            % (s+4)(s+1/3)
G = tf(num, den);

% Blank Bode paper (freq: 0.01–100, mag: -60 to +20 dB, phase: -225° to 0°)
BodePaper(0.01, 100, -60, 20, -225, 0);
hold on

% Exact Bode
bode(G);  grid on

% Asymptotic (straight-line) approximation overlaid
asymp(G, 0.01, 100, -60, 20);
Key Concept
Corner frequencies at ω = 4 rad/s and ω = 1/3 rad/s. Low-frequency asymptote: |G(j0)| = 9 ≈ +19.1 dB. Above each corner the magnitude slope drops −20 dB/decade; at high frequencies the total slope is −40 dB/decade. Phase transitions from 0° to −180°.
Multi-Domain Stability

Nyquist, Nichols & Closed-Loop Frequency Response

For Final Q3, a complete multi-domain stability analysis was performed on a third-order system. The same plant was analyzed using four complementary representations: Nyquist plot (encirclement criterion), open-loop Bode plot (gain/phase margins), Nichols chart (closed-loop M-circles), and closed-loop Bode (bandwidth and peak M). Together they paint a full picture of stability robustness.

FinalQ3.m  ·  Nyquist · Bode · Nichols · Closed-loop Bode 
clear; clc; close all
syms s K
k = 9;

% Transfer functions — open-loop and closed-loop
num_ol = sym(k);
den_ol = expand((s+1)*(3*s+1)*(0.4*s+1));       % three real poles
g_ol   = tf(sym2poly(num_ol), sym2poly(den_ol));

num_cl = expand(k*(0.4*s+1));
den_cl = expand(den_ol + k);                     % 1 + G(s) in denominator
t_cl   = tf(sym2poly(num_cl), sym2poly(den_cl));

% ── 1. Nyquist plot ──────────────────────────────
figure
nyquist(g_ol);  grid on
xlim([-1.2, 0.4]);  ylim([-0.5, 0.5])           % zoom near -1+0j point

% ── 2. Open-loop Bode (gain & phase margins) ────
figure
bode(g_ol);  grid on

% ── 3. Nichols chart ─────────────────────────────
figure
nichols(t_cl);  grid on                          % M-circles & N-circles

% ── 4. Closed-loop Bode (bandwidth, peak M) ──────
figure
bode(t_cl);  grid on

% ── Routh-Hurwitz critical gain ───────────────────
CharEq    = simplify(den_ol + K*num_ol);
CharCoeff = fliplr(coeffs(CharEq, s));
Routh     = simplify(routh_hurwitz(CharCoeff));
K_crit    = solve(Routh(3,1))                   % K_critical ≈ 15.86
Stability Summary — K = 9
Nyquist: locus does not encircle −1+0j → no RHP closed-loop poles, system stable. Gain margin from Bode: K_GM = 15.86 / 9 ≈ 4.9 dB of margin before instability. Routh-Hurwitz independently confirms critical gain at K ≈ 15.86. All four methods in agreement — robust stability at the design operating point.
Compensator Design

Lead-Lag Compensator Design

HW10 demonstrates the full compensator synthesis workflow: translate time-domain specifications (Mp, tp, Kv) into frequency-domain requirements (phase margin, bandwidth, DC gain), then design lead and lag networks to meet them. The lead compensator adds phase at the crossover frequency; the lag compensator boosts low-frequency gain to satisfy the velocity constant without reducing phase margin.

Plant: G(s) = K / [s(s+8)(s+30)]  ·  Specs: Mp = 10 %, tp = 0.6 s, Kv = 10 s⁻¹

HW10.m  ·  Full lead-lag synthesis from specs 
clear; clc; close all
syms s

G_plant_sym = 1/(s*(s+8)*(s+30));

%% ── Step 1: Translate specs to frequency domain ──────────────────────────
M_p = 0.10;   t_p = 0.6;   K_v = 10;

zeta    = abs(log(M_p) / sqrt(pi^2 + log(M_p)^2))          % damping ratio

% Bandwidth from peak-time requirement
omega_B = (pi/(t_p*sqrt(1-zeta^2))) * ...
          sqrt((1-2*zeta^2) + sqrt(4*zeta^4 - 4*zeta^2 + 2))

% DC gain K from velocity constant requirement: Kv = lim s→0 [s·K·G]
K = K_v / limit(s*G_plant_sym, s, 0)

%% ── Step 2: Lead compensator parameters ──────────────────────────────────
phi_max  = atand(2*zeta / sqrt(-2*zeta^2 + sqrt(1 + 4*zeta^4)))
omega_C  = 0.8 * omega_B                           % design crossover frequency
x        = phi_max - (180 - (138.36 + 5))          % required phase addition
alpha    = (1 - sind(x)) / (1 + sind(x))           % pole/zero separation ratio
beta     = 1/alpha

% Lead zero and pole
T_1      = 1 / (omega_C * sqrt(alpha))
z_lead   = 1/T_1
p_lead   = z_lead / alpha

G_lead   = (1/beta) * ((s + z_lead)/(s + p_lead));

%% ── Step 3: Lag compensator parameters ───────────────────────────────────
z_lag    = omega_C / 10                            % decade below crossover
p_lag    = z_lag / beta

G_lag    = beta * ((s + z_lag)/(s + p_lag));

%% ── Step 4: Build compensated open-loop TF ───────────────────────────────
num_plant = K;
den_plant = s*(s+8)*(s+30);

num_comp  = K * (s + z_lag) * (s + z_lead);
den_comp  = (s + p_lag) * (s + p_lead) * den_plant;

G_comp    = tf(sym2poly(num_comp), sym2poly(den_comp));

figure; bode(G_comp); grid on; title('Compensated Open-Loop Bode')

%% ── Step 5: Compare step responses ───────────────────────────────────────
den_uncomp_cl = den_plant + num_plant;
den_comp_cl   = den_comp  + num_comp;

G_unc_cl = tf(sym2poly(num_plant),    sym2poly(den_uncomp_cl));
G_com_cl = tf(sym2poly(num_comp),     sym2poly(den_comp_cl));

figure; hold on
step(G_unc_cl);  uncompensated = stepinfo(G_unc_cl)
step(G_com_cl);  compensated   = stepinfo(G_com_cl)
legend('Uncompensated','Compensated', 'Interpreter','latex')
Design Results
Computed: ζ = 0.591 (from Mp = 10%), ωB ≈ 8.5 rad/s, Kv = 10 sets DC gain K. Lead compensator adds phase at ωC; lag boosts low-frequency gain by factor β. Step response post-compensation shows reduced overshoot, faster settling, and near-zero steady-state error vs. the uncompensated plant — all three specs met simultaneously.
Step Response

Proportional Gain Tuning & Step Response Comparison

A simple but illuminating exercise: comparing two proportional gains on a first-order plant to observe the classic speed vs. overshoot trade-off. The closed-loop transfer function is derived analytically, then step() responses for Kp = 56 and Kp = 36 are overlaid. Higher proportional gain accelerates the response but increases overshoot; lower gain is smoother but sluggish — a fundamental trade-off resolved properly only by adding derivative or integral action.

HW4Q3.m  ·  Proportional gain comparison loop 
clear; clc; close all

% Plant: G(s) = 0.5 / (6s + 2)   (first-order, time constant τ = 3 s)
% Closed-loop with proportional controller K_p:
%   T(s) = 0.5·K_p / (6s + 2 + 0.5·K_p) = K_eff / (τ_eff·s + 1)

K_p = [56, 36];
hold on

for n = 1:length(K_p)
    Kp  = K_p(n);
    num = 0.5*Kp / (2 + 0.5*Kp);          % closed-loop DC gain
    den = [6/(2 + 0.5*Kp),  1];            % [τ_cl, 1]
    T_s = tf(num, den);
    step(T_s);  grid on
end

legend('$K_p = 56$', '$K_p = 36$', ...
       'Interpreter', 'latex', ...
       'Location',    'south', ...
       'Orientation', 'horizontal')
Observation
K_p = 56: τ_cl ≈ 0.18 s, DC gain → 0.93 — fast but higher steady-state offset.  |  K_p = 36: τ_cl ≈ 0.32 s, DC gain → 0.90 — slower, marginally lower gain. Both remain first-order (no overshoot) but neither eliminates steady-state error — integral action (Ki) is required to achieve zero offset.