Control Systems Analysis¶

This notebook demonstrates state-space control systems analysis using Maxima's symbolic capabilities combined with numerical computation from numerics. We analyse linear time-invariant systems in state-space form:

$$\frac{d\mathbf{x}}{dt} = A\mathbf{x} + B\mathbf{u}, \quad \mathbf{y} = C\mathbf{x}$$

Starting from circuit equations, we derive the state-space model symbolically, then use numerical tools for eigenvalue analysis, step response simulation, and frequency response (Bode) plots.

In [13]:
load("numerics")$
load("ax-plots")$

RLC Circuit as State-Space System [S+N]¶

Consider a series RLC circuit driven by an input voltage $V_{\text{in}}$. The two energy-storage elements (inductor $L$ and capacitor $C$) give rise to two state variables: the inductor current $i_L$ and the capacitor voltage $v_C$.

From Kirchhoff's voltage law:

  • $L \frac{di_L}{dt} = -R\,i_L - v_C + V_{\text{in}}$
  • $C \frac{dv_C}{dt} = i_L$

Writing in state-space form with $\mathbf{x} = [i_L,\; v_C]^T$ and $u = V_{\text{in}}$:

In [14]:
/* RLC series circuit: L*di/dt = -R*i - v_C + V_in */
/*                     C*dv_C/dt = i                 */
/* State: x = [i_L, v_C], Input: u = V_in            */
/* dx/dt = A*x + B*u                                 */
A_sym : matrix([-R/L, -1/L], [1/C_cap, 0]);
B_sym : matrix([1/L], [0]);
C_out : matrix([0, 1]);

matrix([-(R/L),-(1/L)],[1/C_cap,0])
matrix([1/L],[0])
Out[14]:
In [15]:
/* Substitute numeric values: R=1 Ohm, L=0.5 H, C=1 F */
params : [R=1, L=1/2, C_cap=1]$
A_num : subst(params, A_sym)$
print("A (numeric) =", A_num)$

A : ndarray(float(A_num))$
B : ndarray(float(subst(params, B_sym)))$
print("A ndarray:", A)$
print("B ndarray:", B)$

A (numeric) = matrix([-2,-2],[1,0])
A ndarray: ndarray(\#\<ndarray\ 1\ DOUBLE\-FLOAT\ \(2\ 2\)\>)
B ndarray: ndarray(\#\<ndarray\ 2\ DOUBLE\-FLOAT\ \(2\ 1\)\>)

Eigenvalue Analysis — Stability [S+N]¶

The stability of a linear system is determined by its eigenvalues. If all eigenvalues have strictly negative real parts, the system is asymptotically stable. The eigenvalues also reveal the natural frequency $\omega_n$ and damping ratio $\zeta$.

For a 2nd-order system with characteristic polynomial $s^2 + 2\zeta\omega_n s + \omega_n^2 = 0$:

  • $\omega_n = \sqrt{\det(A)}$
  • $2\zeta\omega_n = -\text{tr}(A)$
In [16]:
/* Symbolic eigenvalue analysis */
char_poly : charpoly(A_sym, s);
print("Characteristic polynomial:", char_poly)$

/* Solve for eigenvalues symbolically */
evals_sym : solve(char_poly, s);
print("Symbolic eigenvalues:", evals_sym)$

/* Numeric eigenvalue computation via np_eig */
[evals_nd, evecs_nd] : np_eig(A)$
print("Numeric eigenvalues:", evals_nd)$

1/(C_cap*L)-(-s-R/L)*s
Characteristic polynomial: 1/(C_cap*L)-(-s-R/L)*s
[s = -((sqrt(C_cap^2*R^2-4*C_cap*L)+C_cap*R)/(2*C_cap*L)),
 s = (sqrt(C_cap^2*R^2-4*C_cap*L)-C_cap*R)/(2*C_cap*L)]
Symbolic eigenvalues:
                     [s = -((sqrt(C_cap^2*R^2-4*C_cap*L)+C_cap*R)
                        /(2*C_cap*L)),
                      s = (sqrt(C_cap^2*R^2-4*C_cap*L)-C_cap*R)/(2*C_cap*L)]
Numeric eigenvalues: ndarray(\#\<ndarray\ 3\ COMPLEX\-DOUBLE\-FLOAT\ \(2\)\>)
In [17]:
/* Natural frequency and damping ratio */
omega_n : sqrt(subst(params, 1/(L*C_cap)));
zeta : subst(params, R / (2*sqrt(L/C_cap)));

print("Natural frequency omega_n =", omega_n)$
print("Damping ratio zeta =", zeta)$

if zeta < 1 then print("System is underdamped")
elseif zeta = 1 then print("System is critically damped")
else print("System is overdamped")$

/* Check stability: all eigenvalues have negative real parts */
print("All eigenvalues have Re < 0 (stable):",
  is(subst(params, realpart(rhs(evals_sym[1]))) < 0)
  and is(subst(params, realpart(rhs(evals_sym[2]))) < 0))$

sqrt(2)
1/sqrt(2)
Natural frequency omega_n = sqrt(2)
Damping ratio zeta = 1/sqrt(2)
System is underdamped
All eigenvalues have Re < 0 (stable): true

Step Response via Matrix Exponential¶

For a step input $u = 1$ applied at $t = 0$ with zero initial conditions, the state evolves as:

$$\mathbf{x}(t + \Delta t) = e^{A \Delta t}\,\mathbf{x}(t) + \Delta t \cdot B\,u$$

This simple Euler-like time-stepping with the matrix exponential $e^{A\Delta t}$ (computed once via np_expm) gives the step response trajectories for both state variables.

In [18]:
/* Time-stepping for step response */
dt : 0.01$
nt : 500$
t_vals : makelist(i * dt, i, 0, nt - 1)$

/* Precompute matrix exponential for one time step */
eAdt : np_expm(np_scale(dt, A))$

/* Step input as a 1x1 matrix */
u_step : ndarray(matrix([1.0]))$
Bu : np_scale(dt, np_matmul(B, u_step))$

/* Simulate */
x : np_zeros([2, 1])$  /* initial state: [i_L=0, v_C=0] */
i_out : []$
v_out : []$
for k : 1 thru nt do (
  i_out : endcons(np_ref(x, 0, 0), i_out),
  v_out : endcons(np_ref(x, 1, 0), v_out),
  x : np_add(np_matmul(eAdt, x), Bu)
)$

print("Simulation complete:", nt, "time steps")$
print("Final i_L =", last(i_out), ", Final v_C =", last(v_out))$

Simulation complete: 500 time steps
Final i_L = -0.0031410220404779685 , Final v_C = 1.0047288668145948
In [19]:
/* Plot step response: current i_L(t) and voltage v_C(t) */
ax_draw2d(
  color="blue", line_width=2, name="i_L(t) — current",
  lines(t_vals, i_out),
  color="red", line_width=2, name="v_C(t) — voltage",
  lines(t_vals, v_out),
  title="Step Response of RLC Circuit",
  xlabel="Time (s)", ylabel="Amplitude",
  showlegend=true, grid=true
)$
No description has been provided for this image

Frequency Response — Bode Plot [S+N]¶

The transfer function from input $V_{\text{in}}$ to output $v_C$ is:

$$H(s) = C (sI - A)^{-1} B$$

We derive this symbolically, then evaluate $|H(j\omega)|$ and $\angle H(j\omega)$ numerically across a range of frequencies to produce Bode magnitude and phase plots.

In [20]:
/* Derive transfer function symbolically.
   For series RLC with output v_C: impedance divider gives
   H(s) = (1/Cs) / (Ls + R + 1/Cs) = 1 / (LCs^2 + RCs + 1) */
H_sym : 1 / (L * C_cap * s^2 + R * C_cap * s + 1);
print("Transfer function H(s) =", H_sym)$

/* Substitute numeric parameters */
H_num : ratsimp(subst(params, H_sym));
print("H(s) with numeric params =", H_num)$

1/(C_cap*L*s^2+C_cap*R*s+1)
Transfer function H(s) = 1/(C_cap*L*s^2+C_cap*R*s+1)
2/(s^2+2*s+2)
H(s) with numeric params = 2/(s^2+2*s+2)
In [21]:
/* Generate Bode plot data */
omega_vals : np_logspace(-1, 2, 200)$

/* Separate numerator and denominator to avoid symbolic complex division */
H_n : num(H_num)$
H_d : denom(H_num)$

mag_list : []$
phase_list : []$
omega_list : []$
for i : 0 thru np_size(omega_vals) - 1 do block(
  [w, n_val, d_val, d_re, d_im, d2, re, im, mag_db, phase_deg],
  w : np_ref(omega_vals, i),
  n_val : float(subst(s=%i*w, H_n)),
  d_val : subst(s=%i*w, H_d),
  d_re : float(realpart(d_val)),
  d_im : float(imagpart(d_val)),
  d2 : d_re^2 + d_im^2,
  re : n_val * d_re / d2,
  im : -n_val * d_im / d2,
  mag_db : float(10 * log(re^2 + im^2) / log(10)),
  phase_deg : float(atan2(im, re) * 180 / %pi),
  omega_list : endcons(w, omega_list),
  mag_list : endcons(mag_db, mag_list),
  phase_list : endcons(phase_deg, phase_list)
)$

print("Bode data computed for", length(omega_list), "frequencies")$

Bode data computed for 200 frequencies
In [22]:
/* Bode Magnitude Plot */
ax_draw2d(
  color="blue", line_width=2,
  lines(omega_list, mag_list),
  title="Bode Magnitude Plot",
  xlabel="Frequency (rad/s)", ylabel="Magnitude (dB)",
  logx=true, grid=true
)$
No description has been provided for this image
In [23]:
/* Bode Phase Plot */
ax_draw2d(
  color="red", line_width=2,
  lines(omega_list, phase_list),
  title="Bode Phase Plot",
  xlabel="Frequency (rad/s)", ylabel="Phase (degrees)",
  logx=true, grid=true
)$
No description has been provided for this image

Poles and Zeros¶

The poles of the transfer function are the roots of the denominator — they coincide with the eigenvalues of $A$. For stability, all poles must lie in the left half of the complex plane ($\text{Re}(s) < 0$).

In [24]:
/* Find poles: roots of the denominator of H(s) */
H_denom : denom(H_sym)$
poles_sym : solve(subst(params, H_denom), s);
print("Poles:", poles_sym)$

/* Extract real and imaginary parts for plotting */
pole_re : makelist(float(realpart(rhs(p))), p, poles_sym)$
pole_im : makelist(float(imagpart(rhs(p))), p, poles_sym)$

print("Pole locations (Re, Im):")$
for i : 1 thru length(poles_sym) do
  print("  s_", i, "=", pole_re[i], "+", pole_im[i], "j")$

/* Plot poles in the complex plane */
ax_draw2d(
  color="red", marker_size=8, name="Poles",
  points(pole_re, pole_im),
  color="gray", line_width=1,
  explicit(0, x, -3, 1),
  title="Pole-Zero Map",
  xlabel="Real", ylabel="Imaginary",
  grid=true
)$

[s = -%i-1,s = %i-1]
Poles: [s = -%i-1,s = %i-1]
Pole locations (Re, Im):
  s_ 1 = -1.0 + -1.0 j
  s_ 2 = -1.0 + 1.0 j
No description has been provided for this image

Summary¶

The state-space representation provides a systematic framework for analysing linear control systems:

  1. Symbolic derivation from physical equations (RLC circuit) reveals the structure of $A$, $B$, $C$ in terms of component values.
  2. Eigenvalue analysis determines stability, natural frequency, and damping ratio.
  3. Step response via matrix exponential time-stepping shows the transient behaviour.
  4. Bode plots characterise the frequency response, showing resonance and roll-off.
  5. Pole-zero maps give a geometric view of stability in the complex plane.

The symbolic-numeric bridge in Maxima lets us derive general results with parameters and then instantiate them for specific numerical evaluation — combining the best of both worlds.