03. Practical Case Study: Mass-Spring-Damper 2nd-Order Oscillator & Root Locus PID Tuning¶

Mechanical Resonance, Root Locus Parameter Sweeps, and Target Damping Ratio Synthesis

This practical case study analyzes the dynamics of a canonical 2nd-order mechanical Mass-Spring-Damper oscillator, investigates the root locus trajectories under parameter variations, and synthesizes a full PID controller to achieve exact transient performance specifications.

In [ ]:

import matplotlib.pyplot as plt
import numpy as np

import ctrlpy as cp
from ctrlpy.controllers import pid_parallel
from ctrlpy.plotting_plotly import plot_root_locus_plotly

%matplotlib inline

import matplotlib.pyplot as plt import numpy as np import ctrlpy as cp from ctrlpy.controllers import pid_parallel from ctrlpy.plotting_plotly import plot_root_locus_plotly %matplotlib inline

1. Physics & Dynamics of 2nd-Order Mechanical Oscillators¶

A mass $m$ sliding on a surface is attached to a linear spring of stiffness $k$ and a viscous damper with damping coefficient $c$. An external control force $F(t)$ is applied to position the mass:

Equation of Motion (Newton's 2nd Law):¶

$$m\, \ddot{x}(t) + c\, \dot{x}(t) + k\, x(t) = F(t)$$

Taking the Laplace transform with zero initial conditions: $$(m s^2 + c s + k) X(s) = F(s)$$

The force-to-displacement plant transfer function is: $$G(s) = \frac{X(s)}{F(s)} = \frac{1}{m s^2 + c s + k} = \frac{1/m}{s^2 + 2\zeta \omega_n s + \omega_n^2}$$

where:

• Natural Frequency: $\omega_n = \sqrt{\frac{k}{m}}$
• Damping Ratio: $\zeta = \frac{c}{2\sqrt{m k}}$

2. Open-Loop Dynamics & Underdamped Regime¶

Let's consider an underdamped mechanical system with the following parameters:

• Mass: $m = 1.0\text{ kg}$
• Viscous damping: $c = 0.4\text{ N}\cdot\text{s}/\text{m}$
• Spring stiffness: $k = 4.0\text{ N}/\text{m}$

In [ ]:

# Physical parameters
m = 1.0  # kg
c = 0.4  # N*s/m
k = 4.0  # N/m

# Analytical parameters
wn = np.sqrt(k / m)
zeta = c / (2.0 * np.sqrt(m * k))

print(f"Natural Frequency wn : {wn:.4f} rad/s")
print(f"Damping Ratio zeta   : {zeta:.4f} (Underdamped: 0 < zeta < 1)")

# Transfer Function: G(s) = 1 / (m*s^2 + c*s + k)
G_msd = cp.tf([1.0], [m, c, k])

print("\nMass-Spring-Damper Transfer Function:")
print(G_msd)
print(f"Open-Loop Poles: {G_msd.poles()}")

# Physical parameters m = 1.0 # kg c = 0.4 # N*s/m k = 4.0 # N/m # Analytical parameters wn = np.sqrt(k / m) zeta = c / (2.0 * np.sqrt(m * k)) print(f"Natural Frequency wn : {wn:.4f} rad/s") print(f"Damping Ratio zeta : {zeta:.4f} (Underdamped: 0 < zeta < 1)") # Transfer Function: G(s) = 1 / (m*s^2 + c*s + k) G_msd = cp.tf([1.0], [m, c, k]) print("\nMass-Spring-Damper Transfer Function:") print(G_msd) print(f"Open-Loop Poles: {G_msd.poles()}")

In [ ]:

# Simulate open-loop step response
ol_step = cp.step_response(G_msd, T=25.0)

print("--- Open-Loop Step Response ---")
print(f"Steady-State Position: {ol_step.steady_state_value():.4f} m (1/k = {1 / k:.4f})")
print(f"Percent Overshoot    : {ol_step.overshoot():.2f} %")
print(f"Settling Time (2%)   : {ol_step.settling_time(tolerance=0.02):.4f} s")
print(f"Peak Time tp         : {ol_step.peak_time():.4f} s")

fig, ax = cp.plot_step(G_msd, T=25.0)
ax.set_title("Open-Loop Step Response (Lightly Damped Oscillator)")
plt.show()

# Simulate open-loop step response ol_step = cp.step_response(G_msd, T=25.0) print("--- Open-Loop Step Response ---") print(f"Steady-State Position: {ol_step.steady_state_value():.4f} m (1/k = {1 / k:.4f})") print(f"Percent Overshoot : {ol_step.overshoot():.2f} %") print(f"Settling Time (2%) : {ol_step.settling_time(tolerance=0.02):.4f} s") print(f"Peak Time tp : {ol_step.peak_time():.4f} s") fig, ax = cp.plot_step(G_msd, T=25.0) ax.set_title("Open-Loop Step Response (Lightly Damped Oscillator)") plt.show()

3. Root Locus Parameter Sweeps¶

Let's study how the closed-loop pole locations migrate as we vary the proportional feedback gain $K$ in a unity feedback loop:

In [ ]:

# Plot interactive Root Locus for G_msd(s)
fig_rl = plot_root_locus_plotly(G_msd)
fig_rl.show()

# Plot interactive Root Locus for G_msd(s) fig_rl = plot_root_locus_plotly(G_msd) fig_rl.show()

Notice that with pure proportional control, the closed-loop poles travel vertically along the imaginary direction: the natural frequency increases, but the real part remains fixed at $-\zeta\omega_n = -0.2\text{ s}^{-1}$. Thus, pure proportional gain cannot improve the settling time ($t_s \approx \frac{4}{\zeta\omega_n} \approx 20\text{ s}$).

To improve damping and accelerate settling, we need Derivative action ($K_d$) to pull the poles into the left-half plane, and Integral action ($K_i$) to achieve unit tracking.

4. PID Controller Synthesis¶

We design a parallel PID controller: $$C(s) = K_p + \frac{K_i}{s} + K_d s = \frac{K_d s^2 + K_p s + K_i}{s}$$

Target Design Specifications:¶

1. Settling time: $t_s \le 1.5\text{ s}$ (requires $\mathrm{Re}(p) \le -2.67$)
2. Damping ratio: $\zeta \ge 0.707$ (percent overshoot $< 5\%$)
3. Steady-state error: $e_{ss} = 0$ for step position commands.

In [ ]:

# PID controller gains
Kp = 30.0
Ki = 15.0
Kd = 8.0

C_pid = pid_parallel(Kp=Kp, Ki=Ki, Kd=Kd)
print("PID Controller Transfer Function C(s):")
print(C_pid)

# Form closed loop system: T_pid = feedback(C * G, 1)
L_pid = cp.series(C_pid, G_msd)
T_pid = cp.feedback(L_pid, 1.0)

print("\nClosed-Loop Transfer Function T_pid(s):")
print(T_pid)
print(f"Closed-Loop Poles: {T_pid.poles()}")

# PID controller gains Kp = 30.0 Ki = 15.0 Kd = 8.0 C_pid = pid_parallel(Kp=Kp, Ki=Ki, Kd=Kd) print("PID Controller Transfer Function C(s):") print(C_pid) # Form closed loop system: T_pid = feedback(C * G, 1) L_pid = cp.series(C_pid, G_msd) T_pid = cp.feedback(L_pid, 1.0) print("\nClosed-Loop Transfer Function T_pid(s):") print(T_pid) print(f"Closed-Loop Poles: {T_pid.poles()}")

In [ ]:

# Also create P and PD controllers for progressive comparison
C_p = cp.tf([Kp], [1.0])
T_p = cp.feedback(cp.series(C_p, G_msd), 1.0)

C_pd = cp.tf([Kd, Kp], [1.0])
T_pd = cp.feedback(cp.series(C_pd, G_msd), 1.0)

# Simulate responses
t_sim = np.linspace(0.0, 5.0, 1000)
resp_ol = cp.step_response(G_msd * k, T=t_sim)  # Scaled to unit steady-state for comparison
resp_p = cp.step_response(T_p, T=t_sim)
resp_pd = cp.step_response(T_pd, T=t_sim)
resp_pid = cp.step_response(T_pid, T=t_sim)

print("=== PID Closed-Loop Performance ===")
print(f"Steady-State Value : {resp_pid.steady_state_value():.5f} m")
print(f"Rise Time (10%-90%): {resp_pid.rise_time():.4f} s")
print(f"Settling Time (2%) : {resp_pid.settling_time(tolerance=0.02):.4f} s")
print(f"Percent Overshoot  : {resp_pid.overshoot():.2f} %")

# Also create P and PD controllers for progressive comparison C_p = cp.tf([Kp], [1.0]) T_p = cp.feedback(cp.series(C_p, G_msd), 1.0) C_pd = cp.tf([Kd, Kp], [1.0]) T_pd = cp.feedback(cp.series(C_pd, G_msd), 1.0) # Simulate responses t_sim = np.linspace(0.0, 5.0, 1000) resp_ol = cp.step_response(G_msd * k, T=t_sim) # Scaled to unit steady-state for comparison resp_p = cp.step_response(T_p, T=t_sim) resp_pd = cp.step_response(T_pd, T=t_sim) resp_pid = cp.step_response(T_pid, T=t_sim) print("=== PID Closed-Loop Performance ===") print(f"Steady-State Value : {resp_pid.steady_state_value():.5f} m") print(f"Rise Time (10%-90%): {resp_pid.rise_time():.4f} s") print(f"Settling Time (2%) : {resp_pid.settling_time(tolerance=0.02):.4f} s") print(f"Percent Overshoot : {resp_pid.overshoot():.2f} %")

In [ ]:

# Matplotlib Comparison Plot
plt.figure(figsize=(10, 5))
plt.plot(t_sim, resp_ol.y, "k--", label="Open-Loop (Scaled)")
plt.plot(t_sim, resp_p.y, "r-.", label="P-Only (High Oscillation, Off-Target)")
plt.plot(t_sim, resp_pd.y, "g--", label="PD (Well-Damped, but Non-Zero Error)")
plt.plot(t_sim, resp_pid.y, "b-", lw=2.5, label="PID (Fast, Damped, Zero Error)")
plt.axhline(1.0, color="k", linestyle=":", label="Target Position (1.0 m)")
plt.title("Mass-Spring-Damper: Controller Evolution (Open-Loop vs P vs PD vs PID)")
plt.xlabel("Time [s]")
plt.ylabel("Position $x(t)$ [m]")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# Matplotlib Comparison Plot plt.figure(figsize=(10, 5)) plt.plot(t_sim, resp_ol.y, "k--", label="Open-Loop (Scaled)") plt.plot(t_sim, resp_p.y, "r-.", label="P-Only (High Oscillation, Off-Target)") plt.plot(t_sim, resp_pd.y, "g--", label="PD (Well-Damped, but Non-Zero Error)") plt.plot(t_sim, resp_pid.y, "b-", lw=2.5, label="PID (Fast, Damped, Zero Error)") plt.axhline(1.0, color="k", linestyle=":", label="Target Position (1.0 m)") plt.title("Mass-Spring-Damper: Controller Evolution (Open-Loop vs P vs PD vs PID)") plt.xlabel("Time [s]") plt.ylabel("Position $x(t)$ [m]") plt.grid(True) plt.legend() plt.tight_layout() plt.show()

5. State-Space Representation & State Trajectories¶

We can convert the closed-loop system into state-space form to inspect both position $x_1(t)$ and velocity $x_2(t) = \dot{x}(t)$:

In [ ]:

# Convert closed-loop system to StateSpace
T_ss = T_pid.to_ss()

# Simulate step response capturing state vector x(t)
step_ss = cp.step_response(T_ss, T=5.0)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 6), sharex=True)

ax1.plot(step_ss.t, step_ss.y, "b-", lw=2)
ax1.axhline(1.0, color="r", linestyle="--")
ax1.set_ylabel("Position $x(t)$ [m]")
ax1.set_title("Closed-Loop State Trajectories")
ax1.grid(True)

if step_ss.x is not None and step_ss.x.shape[1] >= 2:
    vel = np.gradient(step_ss.y, step_ss.t)
    ax2.plot(step_ss.t, vel, "m-", lw=2)
    ax2.set_ylabel(r"Velocity $\dot{x}(t)$ [m/s]")
    ax2.set_xlabel("Time [s]")
    ax2.grid(True)

plt.tight_layout()
plt.show()

# Convert closed-loop system to StateSpace T_ss = T_pid.to_ss() # Simulate step response capturing state vector x(t) step_ss = cp.step_response(T_ss, T=5.0) fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 6), sharex=True) ax1.plot(step_ss.t, step_ss.y, "b-", lw=2) ax1.axhline(1.0, color="r", linestyle="--") ax1.set_ylabel("Position $x(t)$ [m]") ax1.set_title("Closed-Loop State Trajectories") ax1.grid(True) if step_ss.x is not None and step_ss.x.shape[1] >= 2: vel = np.gradient(step_ss.y, step_ss.t) ax2.plot(step_ss.t, vel, "m-", lw=2) ax2.set_ylabel(r"Velocity $\dot{x}(t)$ [m/s]") ax2.set_xlabel("Time [s]") ax2.grid(True) plt.tight_layout() plt.show()

---

Key Takeaways¶

1. Physical mechanical oscillators can be parameterized by natural frequency $\omega_n$ and damping ratio $\zeta$.
2. Root locus diagrams reveal the limitations of single-gain feedback.
3. Derivative action provides synthetic damping, enabling orders-of-magnitude reduction in settling time.
4. Integral action eliminates steady-state position error against spring restoring forces.

Next, explore 04_frequency_domain_deep_dive.ipynb for advanced Bode margins and Nyquist contour analysis!