PID Controller Tuning with L-BFGS¶

Optimize PID gains $(K_p, K_i, K_d)$ for a second-order plant by minimizing a performance metric. We derive the closed-loop state-space symbolically, simulate the step response numerically via matrix exponential, and tune the gains with np_minimize.

Since the cost function involves a simulation loop, we use central finite-difference gradient approximation rather than analytic gradients.

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

/* Finite-difference gradient helper */
fd_gradient(f, x, eps) := block([nd, g],
  nd : np_size(x), g : np_zeros([nd]),
  for j : 0 thru nd - 1 do block([xp, xm],
    xp : np_copy(x), np_set(xp, j, np_ref(x, j) + eps),
    xm : np_copy(x), np_set(xm, j, np_ref(x, j) - eps),
    np_set(g, j, (f(xp) - f(xm)) / (2 * eps))),
  g)$

The Plant: Second-Order System [S+N]¶

Consider a second-order plant (e.g. motor speed control): $$G(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}$$

With $\omega_n = \sqrt{5}$ and $\zeta = 1/\sqrt{5}$ (underdamped), this gives $G(s) = \frac{5}{s^2 + 2s + 5}$.

In [2]:
/* Plant transfer function */
G_sym : omega_n^2 / (s^2 + 2*zeta*omega_n*s + omega_n^2)$
plant_params : [omega_n = sqrt(5), zeta = 1/sqrt(5)]$
G_num : ratsimp(subst(plant_params, G_sym));

/* State-space: controllable canonical form */
Ap_sym : matrix([0, 1], [-omega_n^2, -2*zeta*omega_n])$
Bp_sym : matrix([0], [omega_n^2])$
Cp_sym : matrix([1, 0])$

Ap : ndarray(float(subst(plant_params, Ap_sym)))$
Bp : ndarray(float(subst(plant_params, Bp_sym)))$
Cp : ndarray(float(subst(plant_params, Cp_sym)))$
print("Plant A =", subst(plant_params, Ap_sym))$

5/(s^2+2*s+5)
Plant A = matrix([0,1],[-5,-2])

Open-Loop Step Response¶

Before adding a controller, let us see the plant's natural step response using $e^{A \Delta t}$ time-stepping.

In [3]:
/* Simulate step response */
dt : 0.01$  nt : 500$
t_sim : makelist(i * dt, i, 0, nt - 1)$
eAdt_p : np_expm(np_scale(dt, Ap))$
Bu_p : np_scale(dt, Bp)$

x_p : np_zeros([2, 1])$
y_open : []$
for k : 1 thru nt do (
  y_open : endcons(np_ref(np_matmul(Cp, x_p), 0, 0), y_open),
  x_p : np_add(np_matmul(eAdt_p, x_p), Bu_p))$

ax_draw2d(
  color="blue", line_width=2, name="Open-loop",
  lines(t_sim, y_open),
  color="black", dash="dash", explicit(1, t, 0, 5),
  title="Open-Loop Step Response",
  xlabel="Time (s)", ylabel="Output y(t)",
  grid=true
)$
No description has been provided for this image

PID Controller: Closed-Loop State Space [S+N]¶

A PID controller $C(s) = K_p + K_i/s + K_d s$ is handled in state-space by augmenting the plant state with an integrator state $x_I$ (integral of err):

$$\dot{x}_I = e = r - y, \quad u = K_p e + K_i x_I + K_d(-\dot{y})$$

The augmented state is $[x_1, x_2, x_I]^T$ and the closed-loop system matrices depend on the gains.

In [4]:
/* Build closed-loop state-space as function of gains */
build_cl_matrices(Kp, Ki, Kd) := block(
  [A_cl, B_cl],
  /* Augmented state: [x1, x2, x_I]
     e = r - x1,  u = Kp*e + Ki*x_I + Kd*(-x2)
     x1' = x2
     x2' = -5*x1 - 2*x2 + 5*(Kp*(r-x1) + Ki*x_I + Kd*(-x2))
          = -(5+5*Kp)*x1 - (2+5*Kd)*x2 + 5*Ki*x_I + 5*Kp*r
     x_I' = r - x1 */
  A_cl : ndarray(float(matrix(
    [0,            1,       0],
    [-(5+5*Kp),    -(2+5*Kd), 5*Ki],
    [-1,           0,       0]))),
  B_cl : ndarray(float(matrix([0], [5*Kp], [1]))),
  [A_cl, B_cl])$
In [5]:
/* Step response simulation for given gains */
simulate_cl(gains) := block(
  [Kp, Ki, Kd, A_cl, B_cl, eAdt, Bu, x_cl, y_out],
  Kp : np_ref(gains, 0),
  Ki : np_ref(gains, 1),
  Kd : np_ref(gains, 2),
  [A_cl, B_cl] : build_cl_matrices(Kp, Ki, Kd),
  eAdt : np_expm(np_scale(dt, A_cl)),
  Bu : np_scale(dt, B_cl),
  x_cl : np_zeros([3, 1]),
  y_out : np_zeros([nt]),
  for k : 0 thru nt - 1 do (
    np_set(y_out, k, np_ref(x_cl, 0, 0)),
    x_cl : np_add(np_matmul(eAdt, x_cl), Bu)),
  y_out)$

/* ISE cost: integral of (1 - y(t))^2 */
ise_cost(gains) := block([y_out, err],
  y_out : simulate_cl(gains),
  err : np_sub(1.0, y_out),
  np_sum(np_pow(err, 2)) * dt)$

ise_grad(gains) := fd_gradient(ise_cost, gains, 1e-4)$

Optimizing PID Gains¶

We start with a moderate initial guess and let L-BFGS find gains that minimize ISE (Integral Squared Error).

In [6]:
gains0 : ndarray([1.0, 0.5, 0.1])$
print("Initial ISE:", ise_cost(gains0))$

[gains_opt, ise_opt, conv] : np_minimize(
  ise_cost, ise_grad, gains0, 1e-6, 100)$
print("Optimized gains: Kp =", np_ref(gains_opt, 0),
      ", Ki =", np_ref(gains_opt, 1),
      ", Kd =", np_ref(gains_opt, 2))$
print("ISE:", ise_opt, " Converged:", conv)$

Initial ISE: 0.6720801307270194
*************************************************
  N=    3   NUMBER OF CORRECTIONS= 3
       INITIAL VALUES
 F=  6.720801307270194D-01   GNORM=  7.783820447813677D-01
*************************************************
   I  NFN     FUNC                    GNORM                   STEPLENGTH
   1    2     4.292460873046933D-01   6.370509763731969D-01   1.284716170811572D+00
  17   18     4.135824936053501D-02   2.085490695328005D-04   1.000000000000000D+00
 THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.
 IFLAG = 0
Optimized gains: Kp = 217.10935222673945 , Ki = 83.7917154842175 , Kd =
                     5.9994785945444296
ISE: 0.04135824936053501  Converged: true
In [7]:
/* Before/after step response comparison */
y_before : simulate_cl(gains0)$
y_after : simulate_cl(gains_opt)$

ax_draw2d(
  color="red", dash="dash", line_width=2, name="Initial gains",
  lines(t_sim, np_to_list(y_before)),
  color="blue", line_width=2, name="Optimized gains",
  lines(t_sim, np_to_list(y_after)),
  color="black", dash="dot", name="Reference (r=1)",
  explicit(1, t, 0, 5),
  title="Step Response: Before vs After Tuning",
  xlabel="Time (s)", ylabel="Output y(t)",
  grid=true, showlegend=true
)$
No description has been provided for this image

Performance Metrics¶

Compare key step-response metrics before and after tuning.

In [8]:
/* Compute overshoot and settling time */
compute_metrics(y_list) := block([ymax, overshoot, settled],
  ymax : lmax(y_list),
  overshoot : (ymax - 1.0) * 100,
  /* Settling: last time |y-1| > 0.02 */
  settled : 0,
  for i : length(y_list) step -1 thru 1 do
    if abs(y_list[i] - 1.0) > 0.02 and settled = 0 then
      settled : i * dt,
  [overshoot, settled])$

[os_b, ts_b] : compute_metrics(np_to_list(y_before))$
[os_a, ts_a] : compute_metrics(np_to_list(y_after))$
print("Before: overshoot =", os_b, "%, settling =", ts_b, "s")$
print("After:  overshoot =", os_a, "%, settling =", ts_a, "s")$

Before: overshoot = -13.262478000353362 %o164, settling = 5.0 s
After:  overshoot = 17.608009138448022 %o164, settling = 0.27 s

ITAE Cost — Alternative Metric¶

The Integral of Time-weighted Absolute Error penalizes late errs more: $$J_{\text{ITAE}} = \int_0^T t\,|e(t)|\,dt$$

Optimizing ITAE typically produces a different tradeoff: less overshoot but possibly slower rise time.

In [9]:
itae_cost(gains) := block([y_out, err, t_nd, te],
  y_out : simulate_cl(gains),
  err : np_abs(np_sub(1.0, y_out)),
  t_nd : np_scale(dt, np_arange(nt)),
  te : np_mul(t_nd, err),
  np_sum(te) * dt)$

itae_grad(gains) := fd_gradient(itae_cost, gains, 1e-4)$

[gains_itae, _, conv_itae] : np_minimize(
  itae_cost, itae_grad, gains0, 1e-6, 100)$
print("ITAE-optimal: Kp =", np_ref(gains_itae, 0),
      ", Ki =", np_ref(gains_itae, 1),
      ", Kd =", np_ref(gains_itae, 2))$

*************************************************
  N=    3   NUMBER OF CORRECTIONS= 3
       INITIAL VALUES
 F=  2.674550592459656D+00   GNORM=  4.556348577799761D+00
*************************************************
   I  NFN     FUNC                    GNORM                   STEPLENGTH
   1    2     5.689865927611260D-01   1.325966209741216D+00   2.194739895170389D-01
  42   66     1.432549548773363D-02   5.937319339606287D-05   1.000000000000000D+00
 THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.
 IFLAG = 0
ITAE-optimal: Kp = 62.3998574830572 , Ki = 13.458982005066446 , Kd =
                  5.654847790622001
In [10]:
y_itae : simulate_cl(gains_itae)$
ax_draw2d(
  color="red", line_width=2, name="ISE-optimal",
  lines(t_sim, np_to_list(y_after)),
  color="green", line_width=2, name="ITAE-optimal",
  lines(t_sim, np_to_list(y_itae)),
  color="black", dash="dot", name="Reference",
  explicit(1, t, 0, 5),
  title="ISE vs ITAE Tuning",
  xlabel="Time (s)", ylabel="Output y(t)",
  grid=true, showlegend=true
)$
No description has been provided for this image

Stability Check [S+N]¶

We verify the optimized closed-loop system is stable by computing eigenvalues of $A_{cl}$. All eigenvalues must have negative real parts.

In [11]:
[A_opt, _] : build_cl_matrices(
  np_ref(gains_opt, 0), np_ref(gains_opt, 1), np_ref(gains_opt, 2))$
[evals_cl, _] : np_eig(A_opt)$

/* Extract real and imaginary parts */
evals_list : np_to_list(evals_cl)$
re_parts : map(realpart, evals_list)$
im_parts : map(imagpart, evals_list)$
print("Closed-loop eigenvalues:", evals_list)$
print("All stable:", every(lambda([r], r < 0), re_parts))$

ax_draw2d(
  color="blue", marker_size=8, name="CL poles",
  points(re_parts, im_parts),
  color="gray", dash="dash",
  explicit(0, u, min(apply(min, re_parts) - 1, -1),
              max(apply(max, re_parts) + 1, 1)),
  title="Pole Map (Optimized Controller)",
  xlabel="Re", ylabel="Im", grid=true
)$

Closed-loop eigenvalues:
                        [28.783422765935292*%i-15.804422117024941,
                         -(28.783422765935292*%i)-15.804422117024941,
                         -0.3885487386722597]
All stable: true
No description has been provided for this image

Summary¶

Metric Initial Gains ISE-Optimal ITAE-Optimal
Cost Manual guess Minimized Minimized
Gradient — Finite differences Finite differences

The finite-difference gradient approach is practical when the cost involves a simulation loop. Each L-BFGS iteration requires $2p+1$ cost evaluations (here $p=3$ gains), each involving a 3$\times$3 matrix exponential and 500 time steps — fast enough for interactive use.