Control in Octave

A control loop reads a sensor, compares it with a target, and drives a motor to close the gap. Before touching hardware you can predict whether that loop settles smoothly, overshoots, or tears itself apart. Octave's control package does this in a few lines, and the maths underneath, polynomials and their roots, is simple enough to check with NumPy.

Transfer functions

A transfer function G(s) describes how a system turns an input into an output, as a ratio of polynomials in s. Think of s as "take a derivative" and 1/s as "integrate". For a DC motor driving a joint, voltage makes the speed rise with a time constant of 0.5 s to 2 rad/s per volt, and the angle is the integral of the speed. Together:

G(s) = 4 / (s * (s + 2)) = 4 / (s^2 + 2*s)

In Octave you list the polynomial coefficients, highest power first: tf(4, [1 2 0]). The denominator [1 2 0] means 1*s^2 + 2*s + 0.

Poles and stability

The poles are the roots of the denominator. Each pole p contributes a term that behaves like exp(p * t) in the response. A negative real part means the term decays. A positive real part means it grows without limit. An imaginary part means it oscillates. So the rule is short: a system is stable if and only if every pole has a negative real part.

The motor alone has poles at 0 and -2. The pole at 0 is the integrator: apply a constant voltage and the angle never settles. Now close the loop: measure the angle, subtract it from the command and feed the error to the motor, which is unity negative feedback. The closed loop is T = G / (1 + G), here 4 / (s^2 + 2*s + 4). Feedback moved the poles, which is the point of control.

pkg load control
G = tf(4, [1 2 0]);            % open loop: 4 / (s^2 + 2 s)
T = feedback(G, 1);            % unity negative feedback
p = pole(T);
wn = abs(p(1));
zeta = -real(p(1)) / wn;
printf("closed-loop poles: %.4f +/- %.4fi\n", real(p(1)), abs(imag(p(1))));
printf("wn = %.3f rad/s, zeta = %.3f\n", wn, zeta);

closed-loop poles: -1.0000 +/- 1.7321i
wn = 2.000 rad/s, zeta = 0.500

Compare with the standard form s^2 + 2*zeta*wn*s + wn^2: wn^2 = 4 gives wn = 2, and 2*zeta*wn = 2 gives zeta = 0.5. The poles are -1 plus or minus 1.7321j. The real part sets how fast the response decays, the imaginary part sets the ringing. The damping ratio zeta predicts the overshoot, exp(-pi*zeta / sqrt(1 - zeta^2)), which is 16.3 percent for zeta = 0.5.

Step response

Apply a unit step to the command. step(T) with no outputs opens a plot; to get numbers, pass a time vector and capture the result.

t = 0:0.01:8;
y = step(T, t);
[ymax, k] = max(y);
printf("overshoot %.1f %%, peak at t = %.1f s\n", (ymax - 1) * 100, t(k));
printf("value at t = 8 s: %.3f\n", y(end));

overshoot 16.3 %, peak at t = 1.8 s
value at t = 8 s: 1.000

The figure from step(T) starts at 0 with zero slope, rises through 0.85 at 1 s, peaks near 1.16 at 1.8 s, dips to about 0.97 at 3.5 s and settles within 2 percent of 1.0 after about 4 s, which is the textbook estimate 4 / (zeta*wn). The final value of 1 means the loop tracks the command exactly, because the integrator in G forces zero steady-state error.

When more gain breaks the loop

Real loops have extra lag from drivers, filters and sensors. Suppose the plant is K / (s*(s+1)*(s+2)). With unity feedback the closed-loop poles are the roots of s^3 + 3*s^2 + 2*s + K, and a larger gain K pushes them toward the right half of the plane.

for K = [4 8]
  p = roots([1 3 2 K]);        % characteristic polynomial
  worst = max(real(p));
  if worst < 0
    verdict = "stable";
  else
    verdict = "unstable";
  end
  printf("K = %d: worst real part %+.4f, %s\n", K, worst, verdict);
end

K = 4: worst real part -0.1018, stable
K = 8: worst real part +0.0832, unstable

At K = 4 the response rings but decays; at K = 8 the oscillation grows every cycle. The boundary is K = 6, where 3 times 2 equals K and the poles sit on the imaginary axis. The command rlocus(G) draws how the poles move as K increases, the standard tool for choosing a safe gain.

The same maths in NumPy

No control package is needed to see the same numbers: poles are just np.roots, and a step response is a differential equation you can integrate. The equation behind T is y'' + 2*y' + 4*y = 4*u. We split it into position y and velocity v and use Euler's method: take a tiny time step dt and move each state by its rate times dt.

numpy_version.py

ResetRun

import numpy as np # Closed loop of G(s) = 4 / (s^2 + 2 s) with unity feedback: T(s) = 4 / (s^2 + 2 s + 4) poles = sorted(np.roots([1.0, 2.0, 4.0]), key=lambda p: -p.imag) print("poles: %.4f %+.4fi and %.4f %+.4fi" % (poles[0].real, poles[0].imag, poles[1].real, poles[1].imag)) wn = abs(poles[0]) zeta = -poles[0].real / wn print("wn = %.3f rad/s, zeta = %.3f" % (wn, zeta)) print("predicted overshoot: %.1f %%" % (100 * np.exp(-np.pi * zeta / np.sqrt(1 - zeta ** 2)))) # Euler integration of y'' + 2 y' + 4 y = 4 u, with a unit step u dt = 0.001 y, v = 0.0, 0.0 wd = np.sqrt(3.0) print(" t euler exact") for n in range(6001): t = n * dt if n % 500 == 0: exact = 1 - np.exp(-t) * (np.cos(wd * t) + np.sin(wd * t) / wd) print("%.1f %.3f %.3f" % (t, y, exact)) a = 4.0 - 2.0 * v - 4.0 * y y += dt * v v += dt * a # Gain sweep: characteristic polynomial s^3 + 3 s^2 + 2 s + K for K in (4, 8): worst = max(np.roots([1, 3, 2, K]).real) print("K = %d: worst real part %+.4f, %s" % (K, worst, "stable" if worst < 0 else "unstable"))

The table compares Euler with the exact formula: they differ by at most 0.001. The 0.5 s sampling misses the true peak of 1.163 at 1.81 s, so the largest printed value is 1.154 at 2.0 s. Euler is first order and crude; a smaller dt, or Octave's ode45, is the cure when accuracy matters. The last two lines repeat the gain sweep with np.roots alone. Change K to 6 and worst becomes essentially zero: a loop on the edge of oscillating forever.