engineeringkit.control_systems#

Domain hub · History · Examples · API reference

Classical stability, compensator design, PID with anti-windup, nonlinear loops, state feedback, LQR and estimation.

engineeringkit.control_systems: feedback, stability, frequency response, state space and optimal control.

Black’s feedback amplifier and Minorsky’s PID controller with Ziegler-Nichols tuning; Maxwell’s characteristic roots, the Routh and Hurwitz tests and Lyapunov’s direct method; Nyquist’s encirclement criterion, Bode plots, stability margins, the Nichols chart and Bode’s sensitivity functions; lead, lag and lead-lag compensators, system type and error constants; Evans’s root locus; Hazen’s second-order servo specifications; Padé dead-time approximation and Smith’s predictor; zero-order-hold and Tustin discretisation and Jury’s test; PID anti-windup; describing functions, limit cycles and the phase plane; Kalman’s controllability and observability, pole placement and Luenberger observers; and the LQR, discrete LQR and Kalman-Bucy Riccati designs. Models are continuous-time and SISO unless stated; SciPy supplies the LTI objects, Riccati/Lyapunov solvers, pole placement and step responses.

class engineeringkit.control_systems.CompensatorResult(kind, compensator, loop, compensated_loop, phase_margin_deg, gain_cross_frequency, gain_margin)[source]#

Bases: object

A lead, lag or lead-lag compensator designed for a loop, with the margins it achieves.

Parameters:

kind (str)
compensator (Any)
loop (Any)
compensated_loop (Any)
phase_margin_deg (float)
gain_cross_frequency (float)
gain_margin (float)

compensated_loop: Any#

Compensated loop \(C(s)L(s)\).

Type:: scipy.signal.TransferFunction

compensator: Any#

Compensator \(C(s)\), unit gain at \(s = 0\).

Type:: scipy.signal.TransferFunction

gain_cross_frequency: float#

Gain-crossover frequency of the compensated loop [rad/s].

Type:: float

gain_margin: float#

Gain margin of the compensated loop [-].

Type:: float

kind: str#

"lead", "lag" or "lead-lag".

Type:: str

loop: Any#

Uncompensated loop \(L(s)\).

Type:: scipy.signal.TransferFunction

phase_margin_deg: float#

Phase margin of the compensated loop [deg].

Type:: float

plot(w=None, axes=None)[source]#: Bode plots of the loop before and after compensation; return (fig, axs).

summary()[source]#

Achieved margins, crossover and the compensator’s high-frequency gain.

Return type:: dict

class engineeringkit.control_systems.PID(kp, ki=0.0, kd=0.0)[source]#

Bases: object

Ideal parallel PID controller \(C(s) = K_p + K_i/s + K_d s\).

Derivative action is ideal (no roll-off filter). The linear methods have no actuator limit; simulate_saturated() adds one, with optional anti-windup. With ki == 0 no integrator is included.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> r = PID(1, 1).simulate(transfer_function([1], [1, 1]), np.linspace(0, 5, 51))
>>> round(r.summary()["final"], 4)
0.9933

Parameters:

kp (float)
ki (float)
kd (float)

closed_loop(plant)[source]#

Unity-feedback closed loop \(CG/(1 + CG)\) from reference to output.

Parameters:: plant (scipy.signal.TransferFunction or tuple) – Plant \(G(s)\), as a transfer function or (num, den) polynomial coefficients in descending powers of \(s\).
Returns:: Closed-loop transfer function.
Return type:: scipy.signal.TransferFunction

kd: float = 0.0#

Derivative gain [input s/error].

Type:: float

ki: float = 0.0#

Integral gain [input/(error s)].

Type:: float

kp: float#

Proportional gain [input/error].

Type:: float

simulate(plant, T=None)[source]#

Unit-step reference response of the closed loop.

Parameters:

plant (scipy.signal.TransferFunction or tuple) – Plant \(G(s)\), as in closed_loop().
T (array_like, optional) – Sample times [s]; chosen by scipy.signal.step() if omitted.

Returns:

Time [s] and output for a unit step in the reference.

Return type:

ResponseResult

simulate_saturated(plant, T, u_min, u_max, anti_windup=None, tracking_time=None, reference=1.0)[source]#

Step response with a saturating actuator, without or with anti-windup.

The controller computes \(v = K_p e + z - K_d\dot y\), the actuator delivers \(u = \min(\max(v, u_{min}), u_{max})\), and the integrator state \(z\) obeys

None: \(\dot z = K_i e\), which winds up while \(u \ne v\);
"clamping" (conditional integration): \(\dot z = 0\) while the actuator is saturated and \(K_i e\) would drive it further, else \(K_i e\);
"back_calculation": \(\dot z = K_i e + (u - v)/T_t\), bleeding the integrator off at the tracking time constant \(T_t\).

Parameters:

plant (scipy.signal.TransferFunction or tuple) – Strictly proper plant \(G(s)\); with kd != 0 its relative degree must be at least 2.
T (array_like) – At least two increasing sample times [s], fine enough for a fixed-step RK4 integration.
u_min (float) – Actuator limits [input units], u_min < u_max.
u_max (float) – Actuator limits [input units], u_min < u_max.
anti_windup ({None, "clamping", "back_calculation"}, optional) – Integrator protection.
tracking_time (float, optional) – \(T_t\) [s] for back-calculation; defaults to \(\sqrt{T_iT_d}\) with \(T_i = K_p/K_i\), \(T_d = K_d/K_p\), or \(T_i\) without derivative action.
reference (float, optional) – Height of the reference step [output units].

Returns:

Time [s], output and actuator output u.

Return type:

ResponseResult

Raises:

ValueError – For an improper or biproper plant, a derivative on a plant of relative degree one, bad limits or an unknown anti-windup mode.

Notes

Derivative action acts on the measurement, so the reference step gives no derivative kick; within the limits the response then differs from simulate() by that kick only. Without an integral gain the integrator is absent and anti_windup has no effect. Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

References

K. J. Åström and T. Hägglund, PID Controllers: Theory, Design, and Tuning, 2nd ed. (ISA, 1995), sec. 3.5.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> G = transfer_function([1], [1, 0])
>>> r = PID(2, 1).simulate_saturated(G, np.linspace(0, 20, 2001), -0.5, 0.5, "back_calculation")
>>> round(r.summary()["final"], 3), r.summary()["input_max"]
(1.0, 0.5)

transfer_function()[source]#: Controller transfer function \(C(s)\); all-zero gains have none.

class engineeringkit.control_systems.PhasePortraitResult(t, x, v, acceleration)[source]#

Bases: object

Trajectory of a second-order system \(\ddot x = f(x, \dot x)\) in its phase plane.

Parameters:

t (ndarray)
x (ndarray)
v (ndarray)
acceleration (Callable)

acceleration: Callable#

acceleration(x, v), the right-hand side \(f\).

Type:: callable

plot(ax=None, field=True)[source]#: Plot the trajectory over the direction field; return (fig, ax).

summary()[source]#

Final state and the largest position excursion.

Return type:: dict

t: ndarray#

Time samples [s].

Type:: ndarray

v: ndarray#

Velocity at each sample, in position units per second.

Type:: ndarray

x: ndarray#

Position at each sample, in the model’s units.

Type:: ndarray

class engineeringkit.control_systems.ResponseResult(t, y, u=None)[source]#

Bases: object

Time response of a simulated control loop.

Parameters:

t (ndarray)
y (ndarray)
u (ndarray | None)

plot(ax=None)[source]#: Plot the response against time; return (fig, ax).

summary()[source]#

Final value and peak of the response, and the input range when recorded.

Return type:: dict

t: ndarray#

Time samples [s].

Type:: ndarray

u: ndarray | None = None#

Actuator output at each sample, when the simulation records it.

Type:: ndarray or None

y: ndarray#

Output at each sample, in output units per unit reference.

Type:: ndarray

engineeringkit.control_systems.as_polynomial(coeffs, name='coeffs', min_degree=0)[source]#

Validate a descending-power polynomial and strip leading zeros.

Parameters:

coeffs (array_like) – Real coefficients, highest power first.
name (str, optional) – Argument name used in error messages.
min_degree (int, optional) – Smallest admissible degree after stripping leading zeros.

Returns:

Coefficients with a nonzero leading term.

Return type:

ndarray

Raises:

ValueError – If the coefficients are not a finite real vector, are all zero, or the degree is below min_degree.

Notes

Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> as_polynomial([0, 1, 3, 2])
array([1., 3., 2.])

engineeringkit.control_systems.bode(sys, w=None)[source]#

Bode magnitude and phase of a continuous-time model.

Parameters:

sys (scipy.signal.lti) – Continuous-time model.
w (array_like, optional) – Angular frequencies [rad/s]; chosen automatically if omitted.

Returns:

(w, magnitude, phase): angular frequency [rad/s], gain [dB] and phase [deg].

Return type:

tuple of ndarray

Notes

Wraps scipy.signal.bode(). See Control systems: matrix contracts and numerical stability.

References

H. W. Bode, Network Analysis and Feedback Amplifier Design (Van Nostrand, 1945).

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> w, mag, phase = bode(transfer_function([1], [1, 1]), [1.0])
>>> round(float(mag[0]), 3), round(float(phase[0]), 1)
(-3.01, -45.0)

engineeringkit.control_systems.complementary_sensitivity(loop)[source]#

Complementary sensitivity \(T = L/(1 + L)\) of a unity-feedback loop.

Parameters:: loop (scipy.signal.lti) – SISO open-loop model \(L(s) = N(s)/D(s)\).
Returns:: \(N/(D + N)\), from the reference to the output and, with a sign change, from sensor noise to the output.
Return type:: scipy.signal.TransferFunction

Notes

\(S + T = 1\) identically. See Control systems: matrix contracts and numerical stability.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> T = complementary_sensitivity(transfer_function([1], [1, 0]))
>>> T.num.tolist(), T.den.tolist()
([1.0], [1.0, 1.0])

engineeringkit.control_systems.controllability_matrix(A, B)[source]#

Kalman controllability matrix \([B,\ AB,\ \dots,\ A^{n-1}B]\).

Parameters:

A (array_like) – State matrix, shape (n, n).
B (array_like) – Input matrix, shape (n, m).

Returns:

Shape (n, n m).

Return type:

ndarray

Notes

See Control systems: matrix contracts and numerical stability.

References

R. E. Kalman, “On the general theory of control systems,” Proc. 1st IFAC Congress (Moscow, 1960) 481-492.

Examples

>>> controllability_matrix([[0, 1], [0, 0]], [[0], [1]]).tolist()
[[0.0, 1.0], [1.0, 0.0]]

engineeringkit.control_systems.dc_gain(sys)[source]#

Steady-state gain \(G(0)\) for a constant input.

Parameters:: sys (scipy.signal.lti) – SISO continuous-time model.
Returns:: Output per unit input at zero frequency.
Return type:: float
Raises:: ValueError – If the model has a pole at \(s = 0\) (an integrator), so the gain is unbounded.

Notes

The steady state is reached only if the system is stable. See Control systems: matrix contracts and numerical stability.

Examples

>>> dc_gain(transfer_function([2, 2], [1, 3, 2]))
1.0

engineeringkit.control_systems.dead_zone_describing_function(amplitude, width, slope=1.0)[source]#

Describing function of a dead zone: zero output for inputs below width, then gain slope.

Parameters:

amplitude (float or array_like) – Input amplitude \(A\) [input units], positive.
width (float) – Half-width \(\delta\) of the dead zone [input units], positive.
slope (float, optional) – Gain \(k\) outside the dead zone [output units per input unit], positive.

Returns:

\(N = 0\) for \(A \le \delta\), else \(k\) minus the saturation describing function with limit \(\delta\).

Return type:

float or ndarray

Notes

A dead zone is a linear gain minus a saturation, and the describing function is linear in the nonlinearity. Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> float(dead_zone_describing_function(0.5, 1.0))
0.0

engineeringkit.control_systems.describing_function(nonlinearity, amplitude, samples=4096)[source]#

Describing function of a memoryless nonlinearity, from the first Fourier harmonic of its output.

Parameters:

nonlinearity (callable) – nonlinearity(x) evaluated elementwise on an array of inputs.
amplitude (float or array_like) – Input amplitude \(A\) [input units], positive.
samples (int, optional) – Points per period of the quadrature, at least 16.

Returns:

\(N(A) = (b_1 + j a_1)/A\), where \(b_1\) and \(a_1\) are the sine and cosine coefficients of the output’s fundamental [output units per input unit].

Return type:

complex or ndarray of complex

Notes

The rectangle rule on a uniform grid is spectrally accurate for smooth nonlinearities and first-order accurate across a jump. An odd, single-valued nonlinearity has a real \(N\). Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> round(float(describing_function(lambda x: x**3, 2.0).real), 6)
3.0

engineeringkit.control_systems.design_lag(loop, phase_margin_deg, extra_deg=5.0, zero_ratio=10.0, w=None)[source]#

Lag compensator that lowers a loop’s crossover to where it has the required phase margin.

Find the frequency \(\omega_c\) where \(\angle L = -180^\circ + \gamma + \delta\), the allowance \(\delta\) paying for the lag’s residual phase there.
Choose \(\beta = |L(j\omega_c)|\), so that the lag’s high-frequency gain \(1/\beta\) makes \(\omega_c\) the new crossover.
Put the lag zero at \(\omega_c/\) zero_ratio.

Parameters:

loop (scipy.signal.lti) – Open-loop model \(L(s)\), with its gain already set for the required error constant.
phase_margin_deg (float) – Required phase margin \(\gamma\) [deg], strictly between 0 and 90.
extra_deg (float, optional) – Allowance \(\delta\) for the lag’s phase at \(\omega_c\) [deg], between 0 and 45.
zero_ratio (float, optional) – Ratio of \(\omega_c\) to the lag zero [-], greater than 1.
w (array_like, optional) – Positive angular frequencies [rad/s] for reading the response.

Returns:

The compensator, both loops and the margins actually achieved.

Return type:

CompensatorResult

Raises:

ValueError – If the phase never reaches the target, or the gain there is already below one (a lag cannot raise it).

Notes

A zero a decade below \(\omega_c\) costs about \(\arctan 0.1 \approx 5.7^\circ\) there, which is why \(\delta = 5^\circ\) is the usual allowance. The price is a slower response: the crossover, and so the bandwidth, goes down. Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> r = design_lag(transfer_function([5], [0.5, 1.5, 1, 0]), 40.0)
>>> round(r.phase_margin_deg, 1), round(r.gain_cross_frequency, 3)
(39.9, 0.564)

engineeringkit.control_systems.design_lead(loop, phase_margin_deg, extra_deg=5.0, w=None)[source]#

Lead compensator that raises a loop’s phase margin, by Ogata’s frequency-response procedure.

Read the present phase margin \(\gamma_0\) of \(L\).
The lead must supply \(\phi_m = \gamma - \gamma_0 + \delta\), where \(\delta\) (extra_deg) pays for the extra lag at the crossover, which moves up.
Put the peak lead at the frequency where \(|L| = \sqrt\alpha\), so that the compensator’s gain \(1/\sqrt\alpha\) makes it the new crossover.

Parameters:

loop (scipy.signal.lti) – Open-loop model \(L(s)\), with its gain already set for the required error constant.
phase_margin_deg (float) – Required phase margin \(\gamma\) [deg], strictly between 0 and 90.
extra_deg (float, optional) – Safety allowance \(\delta\) [deg], between 0 and 45.
w (array_like, optional) – Positive angular frequencies [rad/s] for reading the response; defaults to 20000 points from \(10^{-4}\) to \(10^4\).

Returns:

The compensator, both loops and the margins actually achieved.

Return type:

CompensatorResult

Raises:

ValueError – If the loop has no gain crossover, already meets the margin, or needs 90° or more of lead from one stage.

Notes

A single lead stage gives less than 90° and, in practice, about 60°; beyond that use two stages or design_lead_lag(). The achieved margin differs from \(\gamma\) by however far \(\delta\) misjudges the extra lag. Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> r = design_lead(transfer_function([40], [1, 2, 0]), 50.0)
>>> round(r.phase_margin_deg, 1), round(r.gain_cross_frequency, 2)
(49.8, 8.85)

engineeringkit.control_systems.design_lead_lag(loop, phase_margin_deg, crossover_frequency, extra_deg=5.0, zero_ratio=10.0, w=None)[source]#

Lead-lag compensator giving a chosen crossover frequency and phase margin.

The lead puts its peak at the chosen crossover \(\omega_c\) and supplies \(\phi_m = \gamma + \delta - 180^\circ - \angle L(j\omega_c)\); the lag then sets \(\beta = |L(j\omega_c)|/\sqrt\alpha\) so that the total gain there is one.

Parameters:

loop (scipy.signal.lti) – Open-loop model \(L(s)\), with its gain already set for the required error constant.
phase_margin_deg (float) – Required phase margin \(\gamma\) [deg], strictly between 0 and 90.
crossover_frequency (float) – Required gain-crossover frequency \(\omega_c\) [rad/s], positive.
extra_deg (float, optional) – Allowance \(\delta\) for the lag’s phase at \(\omega_c\) [deg], between 0 and 45.
zero_ratio (float, optional) – Ratio of \(\omega_c\) to the lag zero [-], greater than 1.
w (array_like, optional) – Positive angular frequencies [rad/s] for reading the achieved margins.

Returns:

The compensator, both loops and the margins actually achieved.

Return type:

CompensatorResult

Raises:

ValueError – If the required lead is not between 0 and 90°, or the loop gain at \(\omega_c\) is too low for a lag to be needed (then use design_lead()).

Notes

The lead buys phase at the crossover, the lag gives back the excess gain the lead adds, so the error constants are unchanged and the crossover lands where it was asked. The phase of \(L\) is read in \([-360^\circ, 0^\circ)\). Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> r = design_lead_lag(transfer_function([20], [1, 3, 2, 0]), 45.0, 1.0)
>>> round(r.phase_margin_deg, 1), round(r.gain_cross_frequency, 3)
(44.6, 1.004)

engineeringkit.control_systems.discretize_tustin(sys, dt, prewarp_frequency=None)[source]#

Tustin’s bilinear discretisation of a transfer function, optionally prewarped.

Every \(s\) in \(G(s)\) is replaced by \(c\,(z - 1)/(z + 1)\) with \(c = 2/T\), the trapezoidal rule for integration, or \(c = \omega_0/\tan(\omega_0 T/2)\) to make the discrete response match the continuous one exactly at \(\omega_0\).

Parameters:

sys (scipy.signal.lti) – Proper SISO continuous-time model.
dt (float) – Sampling period \(T\) [s], positive.
prewarp_frequency (float, optional) – Frequency \(\omega_0\) [rad/s] to match, positive and below the Nyquist frequency \(\pi/T\).

Returns:

Discrete-time model in descending powers of \(z\) with dt set, normalised to a monic denominator.

Return type:

scipy.signal.TransferFunction

Raises:

ValueError – If the model is improper or the prewarping frequency is not below \(\pi/T\).

Notes

A pole \(p\) maps to \(z = (c + p)/(c - p)\), so a stable model stays stable and the imaginary axis lands on the unit circle. The whole frequency axis is squeezed into \([0, \pi/T)\): the continuous frequency \(\omega\) appears at \(\omega_d = (2/T)\arctan(\omega/c)\), which is why prewarping is needed near \(\pi/T\). Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

References

A. Tustin, “A method of analysing the behaviour of linear systems in terms of time series,” J. IEE Part IIA 94 (1947) 130-142.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> Gd = discretize_tustin(transfer_function([1], [1, 0]), 0.1)
>>> Gd.num.tolist(), Gd.den.tolist()
([0.05, 0.05], [1.0, -1.0])

engineeringkit.control_systems.discretize_zoh(A, B, dt)[source]#

Exact zero-order-hold discretisation \((A_d, B_d)\) of \(\dot x = Ax + Bu\).

Parameters:

A (array_like) – Continuous state matrix, shape (n, n) [1/s].
B (array_like) – Continuous input matrix, shape (n, m).
dt (float) – Sampling period [s], positive.

Returns:

(Ad, Bd) for \(x_{k+1} = A_d x_k + B_d u_k\).

Return type:

tuple of ndarray

Notes

Both matrices come from one matrix exponential of the augmented matrix \(\begin{bmatrix}A & B\\ 0 & 0\end{bmatrix} T\), which is valid even when \(A\) is singular. See Control systems: matrix contracts and numerical stability.

References

J. R. Ragazzini and L. A. Zadeh, “The analysis of sampled-data systems,” Trans. AIEE 71 (1952) 225-234.

Examples

>>> Ad, Bd = discretize_zoh([[0.0]], [[1.0]], 0.5)
>>> Ad.tolist(), Bd.tolist()
([[1.0]], [[0.5]])

engineeringkit.control_systems.dlqr(A, B, Q, R)[source]#

Discrete-time infinite-horizon linear-quadratic regulator.

Minimises \(\sum_{k \ge 0} (x_k^\mathsf{T} Q x_k + u_k^\mathsf{T} R u_k)\) over \(u_k = -K x_k\) for \(x_{k+1} = A x_k + B u_k\); the steady state of Bellman’s backward dynamic-programming recursion.

Parameters:

A (array_like) – Discrete state matrix, shape (n, n).
B (array_like) – Discrete input matrix, shape (n, m).
Q (array_like) – Symmetric positive semidefinite state weight, shape (n, n).
R (array_like) – Symmetric positive definite input weight, shape (m, m).

Returns:

(K, P, poles): gain \(K = (R + B^\mathsf{T}PB)^{-1}B^\mathsf{T}PA\), the discrete Riccati solution \(P\) and the closed-loop eigenvalues (inside the unit circle when stabilising).

Return type:

tuple

Notes

See Control systems: matrix contracts and numerical stability.

References

Bellman, Dynamic Programming (Princeton University Press, 1957).

Examples

>>> K, P, poles = dlqr([[1.0]], [[1.0]], [[1.0]], [[1.0]])
>>> round(float(P[0, 0]), 6), round(float(K[0, 0]), 6)
(1.618034, 0.618034)

engineeringkit.control_systems.error_constants(loop)[source]#

Static error constants and steady-state errors of a stable unity-feedback loop.

Parameters:: loop (scipy.signal.lti) – SISO open-loop model \(L(s)\); the closed loop \(L/(1+L)\) must be stable.
Returns:: type (int); Kp [-], Kv [1/s] and Ka [1/s²]; and the steady-state errors step_error \(1/(1+K_p)\), ramp_error \(1/K_v\) [s] and parabola_error \(1/K_a\) [s²], per unit step, unit-slope ramp and \(t^2/2\) input. Constants are inf above the system type and 0 below it; errors are then 0 and inf.
Return type:: dict
Raises:: ValueError – If the closed loop is not asymptotically stable, so that no steady state exists.

Notes

A type-\(N\) loop follows a polynomial input of degree below \(N\) with zero error, one of degree \(N\) with a finite error, and loses higher ones. Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> e = error_constants(transfer_function([10], [1, 5, 0]))
>>> e["type"], e["Kv"], e["ramp_error"]
(1, 2.0, 0.5)

engineeringkit.control_systems.feedback(G, H=1.0, sign=-1)[source]#

Closed-loop transfer function of a feedback loop, \(G/(1 - \sigma GH)\).

Parameters:

G (scipy.signal.lti) – Forward-path model.
H (scipy.signal.lti or float, optional) – Feedback-path model or constant gain (unity feedback by default).
sign ({-1, 1}, optional) – Feedback sign \(\sigma\); -1 is negative feedback.

Returns:

Closed-loop model from reference to output.

Return type:

scipy.signal.TransferFunction

Notes

Pole-zero cancellations are not removed. See Control systems: matrix contracts and numerical stability.

References

H. S. Black, “Stabilized feedback amplifiers,” Bell Syst. Tech. J. 13 (1934) 1-18.

Examples

>>> T = feedback(transfer_function([1000], [1]), 0.01)
>>> round(dc_gain(T), 3)
90.909

engineeringkit.control_systems.feedforward_gain(A, B, C, K)[source]#

Reference precompensator \(N\) giving unit DC gain under SISO state feedback.

With \(u = -Kx + Nr\), the closed loop \(\dot x = (A - BK)x + BNr\) tracks a constant reference when \(N = -1/\left(C (A - BK)^{-1} B\right)\).

Parameters:

A (array_like) – State matrix, shape (n, n).
B (array_like) – Input matrix, shape (n, 1).
C (array_like) – Output matrix, shape (1, n).
K (array_like) – State-feedback gain, shape (1, n).

Returns:

Precompensator gain [input/reference units].

Return type:

float

Raises:

ValueError – If the system is not SISO, or the closed loop has a zero at the origin so no finite \(N\) exists.

Notes

Exact only for a perfect model; integral action is needed for robust tracking. See Control systems: matrix contracts and numerical stability.

Examples

>>> feedforward_gain([[-1]], [[1]], [[1]], [[2]])
3.0

engineeringkit.control_systems.hurwitz_determinants(coeffs)[source]#

Leading principal minors \(\Delta_1, \dots, \Delta_n\) of the Hurwitz matrix.

For \(a_0 s^n + a_1 s^{n-1} + \dots + a_n\) with \(a_0 > 0\), all roots lie in the open left half-plane if and only if every \(\Delta_k > 0\). The Hurwitz matrix has entries \(H_{ij} = a_{2j-i}\) (1-based, \(a_k = 0\) outside \(0 \le k \le n\)).

Parameters:: coeffs (array_like) – Real coefficients in descending powers of \(s\), degree at least one. They are normalised by the leading coefficient.
Returns:: The \(n\) determinants.
Return type:: ndarray

Notes

On the stability boundary \(\Delta_{n-1} = 0\) (a pair of imaginary roots) or \(\Delta_n = 0\) (a root at the origin). See Control systems: matrix contracts and numerical stability.

References

A. Hurwitz, “Ueber die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt,” Math. Ann. 46 (1895) 273-284.

Examples

>>> hurwitz_determinants([1, 3, 3, 1]).round(6).tolist()
[3.0, 8.0, 8.0]

engineeringkit.control_systems.impulse_response(sys, T=None)[source]#

Response to a unit impulse input from rest.

Parameters:

sys (scipy.signal.lti) – Continuous-time model.
T (array_like, optional) – Time samples [s]; chosen automatically if omitted.

Returns:

(t, h): time [s] and impulse response [output/(input s)].

Return type:

tuple of ndarray

Notes

Wraps scipy.signal.impulse(). See Control systems: matrix contracts and numerical stability.

Examples

>>> t, h = impulse_response(transfer_function([1], [1, 1]), [0.0, 1.0])
>>> round(float(h[-1]), 4)
0.3679

engineeringkit.control_systems.is_controllable(A, B)[source]#

Whether the controllability matrix has full state rank.

Parameters:

A (array_like) – State matrix, shape (n, n).
B (array_like) – Input matrix, shape (n, m).

Returns:

True if every state can be steered by the inputs.

Return type:

bool

Notes

Rank is computed numerically with SVD tolerances. See Control systems: matrix contracts and numerical stability.

Examples

>>> is_controllable([[-1, 0], [0, -1]], [[1], [1]])
False

engineeringkit.control_systems.is_observable(A, C)[source]#

Whether the observability matrix has full state rank.

Parameters:

A (array_like) – State matrix, shape (n, n).
C (array_like) – Output matrix, shape (p, n).

Returns:

True if every state can be reconstructed from the outputs.

Return type:

bool

Notes

See Control systems: matrix contracts and numerical stability.

Examples

>>> is_observable([[0, 1], [0, 0]], [[1, 0]])
True

engineeringkit.control_systems.jury_stable(coeffs)[source]#

Whether every root of a polynomial in \(z\) lies strictly inside the unit circle.

Jury’s table reduces \(P(z) = a_0 z^n + a_1 z^{n-1} + \dots + a_n\) one degree at a time. Each row is \(a_0 P(z) - a_n z^n P(1/z)\), divided by \(z\), and the polynomial is stable if and only if \(|a_n| < |a_0|\) holds at every step down to degree zero.

Parameters:: coeffs (array_like) – Real coefficients in descending powers of \(z\), degree at least one.
Returns:: True if all roots satisfy \(|z| < 1\); roots on the unit circle count as unstable.
Return type:: bool

Notes

This is the discrete counterpart of Routh’s test: no roots are computed. The quick necessary conditions \(P(1) > 0\) and \((-1)^n P(-1) > 0\) (for \(a_0 > 0\)) follow from it. Each row is normalised to unit leading coefficient, and a relative tolerance of \(10^{-12}\) decides the boundary cases. See Control systems: matrix contracts and numerical stability.

References

E. I. Jury and J. Blanchard, “A stability test for linear discrete systems in table form,” Proc. IRE 49 (1961) 1947-1948.

Examples

>>> jury_stable([1, -1.0, 0.5]), jury_stable([1, -2.5, 1.0])
(True, False)

engineeringkit.control_systems.kalman_gain(A, C, Q, R)[source]#

Steady-state Kalman-Bucy estimator gain for continuous-time measurements.

For \(\dot x = Ax + w\), \(y = Cx + v\) with white process and measurement noise of spectral densities \(Q\) and \(R\), the estimator \(\dot{\hat x} = A\hat x + L(y - C\hat x)\) with \(L = PC^\mathsf{T}R^{-1}\) minimises the steady error covariance \(P\).

Parameters:

A (array_like) – State matrix, shape (n, n) [1/s].
C (array_like) – Output matrix, shape (p, n).
Q (array_like) – Process-noise spectral density, symmetric positive semidefinite, shape (n, n).
R (array_like) – Measurement-noise spectral density, symmetric positive definite, shape (p, p).

Returns:

(L, P, poles): observer gain, error covariance and estimator poles [1/s].

Return type:

tuple

Notes

\(P\) solves \(AP + PA^\mathsf{T} - PC^\mathsf{T}R^{-1}CP + Q = 0\). \(Q\) and \(R\) are spectral densities, not sampled covariances. See Control systems: matrix contracts and numerical stability.

References

R. E. Kalman and R. S. Bucy, “New results in linear filtering and prediction theory,” J. Basic Eng. 83 (1961) 95-108.

Examples

>>> L, P, poles = kalman_gain([[0]], [[1]], [[4]], [[1]])
>>> L.round(6).tolist(), P.round(6).tolist()
([[2.0]], [[2.0]])

engineeringkit.control_systems.lag_compensator(beta, zero_frequency)[source]#

Phase-lag compensator \((Ts + 1)/(\beta Ts + 1)\) that cuts the high-frequency gain by \(\beta\).

Parameters:

beta (float) – Pole-zero ratio \(\beta\) [-], greater than 1.
zero_frequency (float) – Zero \(1/T\) [rad/s], positive; the pole is at \(1/(\beta T)\).

Returns:

\(C(s)\) with unit gain at DC and \(1/\beta\) at high frequency.

Return type:

scipy.signal.TransferFunction

Notes

Placed a decade or so below the crossover, the lag lowers the loop gain there by \(\beta\) while adding only a few degrees of lag. Equivalently, it raises the low-frequency gain, and so the error constants, by \(\beta\) relative to the crossover gain. Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> C = lag_compensator(10.0, 0.1)
>>> C.num.tolist(), C.den.tolist()
([0.1, 0.01], [1.0, 0.01])

engineeringkit.control_systems.lead_compensator(phase_lead_deg, frequency)[source]#

Phase-lead compensator \((Ts + 1)/(\alpha Ts + 1)\) with its peak lead at a given frequency.

Parameters:

phase_lead_deg (float) – Maximum phase lead \(\phi_m\) [deg], strictly between 0 and 90.
frequency (float) – Frequency of the maximum lead, \(\omega_m\) [rad/s], positive.

Returns:

\(C(s)\) with \(\alpha = (1 - \sin\phi_m)/(1 + \sin\phi_m)\) and \(T = 1/(\omega_m\sqrt\alpha)\): unit gain at DC, gain \(1/\sqrt\alpha\) at \(\omega_m\) and \(1/\alpha\) at high frequency.

Return type:

scipy.signal.TransferFunction

Notes

The zero \(1/T\) and pole \(1/(\alpha T)\) sit symmetrically about \(\omega_m\) on a log scale. Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> C = lead_compensator(30.0, 2.0)
>>> [round(float(c), 4) for c in C.num / C.den[-1]], [round(float(c), 4) for c in C.den / C.den[-1]]
([0.866, 1.0], [0.2887, 1.0])

engineeringkit.control_systems.lead_lag_compensator(phase_lead_deg, lead_frequency, beta, lag_zero_frequency)[source]#

Series lead-lag compensator, the product of lead_compensator() and lag_compensator().

Parameters:

phase_lead_deg (float) – Maximum phase lead of the lead part [deg], strictly between 0 and 90.
lead_frequency (float) – Frequency of that maximum lead [rad/s], positive.
beta (float) – Pole-zero ratio of the lag part [-], greater than 1.
lag_zero_frequency (float) – Zero of the lag part [rad/s], positive.

Returns:

\(C(s)\) with unit gain at DC.

Return type:

scipy.signal.TransferFunction

Notes

Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> C = lead_lag_compensator(30.0, 2.0, 10.0, 0.1)
>>> len(C.den) - 1
2

engineeringkit.control_systems.limit_cycle(loop, describing, amplitude_bounds, w=None, samples=400)[source]#

Limit cycles predicted by the describing function for a loop closed through a real nonlinearity.

With a real \(N(A)\), \(L(j\omega) = -1/N(A)\) needs \(\angle L = -180^\circ\) (a phase crossover \(\omega_{180}\)) and \(N(A) = 1/|L(j\omega_{180})|\), the gain margin there.

Parameters:

loop (scipy.signal.lti) – Linear part \(L(s)\) of the loop, open-loop stable.
describing (callable) – describing(A) returning the real describing function, vectorised over amplitudes.
amplitude_bounds (tuple of float) – (A_min, A_max) [input units], the positive amplitude range searched.
w (array_like, optional) – Positive angular frequencies [rad/s] for finding the phase crossovers, as in stability_margins().
samples (int, optional) – Log-spaced amplitudes scanned for sign changes, at least 2.

Returns:

One dict per predicted limit cycle, in increasing frequency then amplitude: amplitude [input units], frequency [rad/s], period [s] and stable (bool). Empty if none is predicted.

Return type:

list of dict

Raises:

ValueError – If the bounds are not increasing and positive, or the describing function is not real.

Notes

By Loeb’s criterion, a cycle is stable when \(N\) falls through the gain margin as \(A\) grows: a larger oscillation then sees a loop gain below the critical one and decays back. The method assumes \(L\) filters out the harmonics, so it gives estimates, not proofs. Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

References

R. J. Kochenburger, “A frequency response method for analyzing and synthesizing contactor servomechanisms,” Trans. AIEE 69 (1950) 270-284.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> relay = lambda A: relay_describing_function(A, 1.0)
>>> (cycle,) = limit_cycle(transfer_function([1], [1, 3, 3, 1]), relay, (0.01, 10))
>>> round(cycle["amplitude"], 4), round(cycle["frequency"], 3)
(0.1592, 1.732)

engineeringkit.control_systems.lqr(A, B, Q, R)[source]#

Continuous-time linear-quadratic regulator.

Minimises \(J = \int_0^\infty (x^\mathsf{T} Q x + u^\mathsf{T} R u)\,dt\) over \(u = -Kx\) for \(\dot x = Ax + Bu\).

Parameters:

A (array_like) – State matrix, shape (n, n) [1/s].
B (array_like) – Input matrix, shape (n, m).
Q (array_like) – Symmetric positive semidefinite state weight, shape (n, n).
R (array_like) – Symmetric positive definite input weight, shape (m, m).

Returns:

(K, P, poles): gain \(K = R^{-1} B^\mathsf{T} P\), the Riccati solution \(P\) (so the optimal cost from \(x_0\) is \(x_0^\mathsf{T} P x_0\)) and the closed-loop poles [1/s].

Return type:

tuple

Notes

\(P\) solves \(A^\mathsf{T}P + PA - PBR^{-1}B^\mathsf{T}P + Q = 0\). A stabilising solution needs \((A, B)\) stabilisable and \((Q^{1/2}, A)\) detectable; SciPy errors are not masked. See Control systems: matrix contracts and numerical stability.

References

R. E. Kalman, “Contributions to the theory of optimal control,” Bol. Soc. Mat. Mexicana 5 (1960) 102-119.

Examples

>>> K, P, poles = lqr([[0]], [[1]], [[4]], [[1]])
>>> K.round(6).tolist(), poles.real.round(6).tolist()
([[2.0]], [-2.0])

engineeringkit.control_systems.lyapunov_equation(A, Q=None)[source]#

Solve the continuous Lyapunov equation \(A^\mathsf{T} P + P A = -Q\).

If \(A\) is stable, \(V(x) = x^\mathsf{T} P x\) is a Lyapunov function: it is positive and decreases at the rate \(\dot V = -x^\mathsf{T} Q x\) along every trajectory of \(\dot x = Ax\).

Parameters:

A (array_like) – State matrix, shape (n, n) [1/s].
Q (array_like, optional) – Symmetric positive semidefinite decay-rate matrix, shape (n, n); the identity by default.

Returns:

Symmetric solution \(P\).

Return type:

ndarray

Raises:

ValueError – If two eigenvalues of \(A\) sum to zero, so the equation has no unique solution.

Notes

See Control systems: matrix contracts and numerical stability.

References

A. M. Lyapunov, The General Problem of the Stability of Motion (Kharkov, 1892; English translation, Taylor & Francis, 1992).

Examples

>>> lyapunov_equation([[-1.0]]).tolist()
[[0.5]]

engineeringkit.control_systems.lyapunov_stable(A)[source]#

Asymptotic stability of \(\dot x = Ax\) by Lyapunov’s direct method.

Parameters:: A (array_like) – State matrix, shape (n, n) [1/s].
Returns:: True if \(A^\mathsf{T} P + PA = -I\) has a positive definite solution \(P\).
Return type:: bool

Notes

A singular equation (eigenvalues summing to zero, e.g. an undamped oscillator) means the origin is not asymptotically stable. See Control systems: matrix contracts and numerical stability.

Examples

>>> lyapunov_stable([[0, 1], [-2, -3]]), lyapunov_stable([[0, 1], [-1, 0]])
(True, False)

engineeringkit.control_systems.nichols(sys, w=None)[source]#

Open-loop Nichols locus of a loop, with the closed-loop response and resonant peak it implies.

Parameters:

sys (scipy.signal.lti) – SISO open-loop model \(L(s)\).
w (array_like, optional) – Positive angular frequencies [rad/s]; defaults to 4000 points from \(10^{-3}\) to \(10^3\).

Returns:

w [rad/s]; the open-loop gain_db [dB] and continuous phase_deg [deg]; the closed-loop closed_gain_db [dB] and closed_phase_deg [deg] of \(L/(1 + L)\); and the resonant peak resonant_peak_db [dB] with its resonant_frequency [rad/s], the largest closed-loop gain on the grid.

Return type:

dict

Notes

The resonant peak is the highest M-contour the locus touches. The open-loop phase is unwrapped from the lowest frequency, so it can leave \([-360^\circ, 0^\circ]\). Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

References

H. M. James, N. B. Nichols and R. S. Phillips, Theory of Servomechanisms, MIT Radiation Laboratory Series 25 (McGraw-Hill, 1947).

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> n = nichols(transfer_function([1], [1, 1, 0]))
>>> round(n["resonant_peak_db"], 2), round(n["resonant_frequency"], 3)
(1.25, 0.707)

engineeringkit.control_systems.nichols_closed_loop(gain_db, phase_deg)[source]#

Closed-loop gain and phase of a unity-feedback loop from its open-loop Nichols coordinates.

Parameters:

gain_db (float or array_like) – Open-loop gain \(20\log_{10}|L|\) [dB].
phase_deg (float or array_like) – Open-loop phase \(\angle L\) [deg].

Returns:

(closed_gain_db, closed_phase_deg) of \(T = L/(1 + L)\), the phase in \((-180^\circ, 180^\circ]\). The gain grows without bound towards the critical point \(L = -1\).

Return type:

tuple of ndarray

Notes

The contours of constant \(|T|\) drawn in these coordinates are the M-contours of the Nichols chart; they repeat every \(360^\circ\) of open-loop phase. Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

References

H. M. James, N. B. Nichols and R. S. Phillips, Theory of Servomechanisms, MIT Radiation Laboratory Series 25 (McGraw-Hill, 1947).

Examples

>>> m, p = nichols_closed_loop(0.0, -90.0)
>>> round(float(m), 3), round(float(p), 1)
(-3.01, -45.0)

engineeringkit.control_systems.num_den(sys)[source]#

Numerator and denominator polynomials of a SISO continuous-time model.

Parameters:: sys (scipy.signal.lti) – A TransferFunction, or any SciPy LTI model that converts to one (state space, zeros-poles-gain).
Returns:: (num, den) in descending powers of \(s\).
Return type:: tuple of ndarray
Raises:: ValueError – If the model is not single-input single-output.

Notes

Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> num_den(signal.TransferFunction([2, 2], [1, 3, 2]))
(array([2., 2.]), array([1., 3., 2.]))

engineeringkit.control_systems.nyquist(sys, w)[source]#

Complex frequency response \(G(j\omega)\).

Parameters:

sys (scipy.signal.lti) – SISO continuous-time model.
w (float or array_like) – Angular frequencies [rad/s]; negative values give the mirror image.

Returns:

Output/input gain at each frequency.

Return type:

complex or ndarray of complex

Notes

See Control systems: matrix contracts and numerical stability.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> complex(np.round(nyquist(transfer_function([2, 2], [1, 3, 2]), 2.0), 12))
(0.5-0.5j)

engineeringkit.control_systems.nyquist_encirclements(sys, w=None)[source]#

Net clockwise encirclements of \(-1\) by the Nyquist plot of a loop.

Nyquist’s criterion: the closed loop \(L/(1+L)\) has \(Z = N + P\) right-half-plane poles, where \(P\) is the number of open-loop right-half-plane poles and \(N\) the value returned here.

Parameters:

sys (scipy.signal.lti) – Proper open-loop model \(L(s)\) with no poles on the imaginary axis.
w (array_like, optional) – Positive angular frequencies [rad/s] spanning all features of \(L\); defaults to 40001 points from \(10^{-4}\) to \(10^4\).

Returns:

Clockwise encirclements of \(-1\) (negative if counter-clockwise).

Return type:

int

Raises:

ValueError – If \(L\) is improper, has imaginary-axis poles, or passes through \(-1\) on the grid.

Notes

The plot is traced for \(\omega\) from \(-\infty\) to \(\infty\); for a proper loop the infinite semicircle maps to a single point. A grid that misses a fast phase swing can miscount. See Control systems: matrix contracts and numerical stability.

References

H. Nyquist, “Regeneration theory,” Bell Syst. Tech. J. 11 (1932) 126-147.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> nyquist_encirclements(transfer_function([12], [1, 3, 3, 1]))
2

engineeringkit.control_systems.observability_matrix(A, C)[source]#

Kalman observability matrix \([C;\ CA;\ \dots;\ CA^{n-1}]\).

Parameters:

A (array_like) – State matrix, shape (n, n).
C (array_like) – Output matrix, shape (p, n).

Returns:

Shape (n p, n).

Return type:

ndarray

Notes

See Control systems: matrix contracts and numerical stability.

Examples

>>> observability_matrix([[0, 1], [0, 0]], [[1, 0]]).tolist()
[[1.0, 0.0], [0.0, 1.0]]

engineeringkit.control_systems.observer_gain(A, C, desired_poles)[source]#

Luenberger observer gain \(L\) placing the eigenvalues of \(A - LC\).

The observer \(\dot{\hat x} = A\hat x + Bu + L(y - C\hat x)\) makes the estimation error obey \(\dot e = (A - LC)e\).

Parameters:

A (array_like) – State matrix, shape (n, n).
C (array_like) – Output matrix, shape (p, n); (A, C) must be observable.
desired_poles (array_like) – n error-dynamics poles [1/s], conjugate pairs for complex values.

Returns:

Gain \(L\), shape (n, p).

Return type:

ndarray

Notes

Computed by duality, placing the poles of \(A^\mathsf{T} - C^\mathsf{T} L^\mathsf{T}\), so with one output a repeated pole is placed by Ackermann’s formula, as in state_feedback_poles(). See Control systems: matrix contracts and numerical stability.

References

D. G. Luenberger, “Observing the state of a linear system,” IEEE Trans. Mil. Electron. 8 (1964) 74-80.

Examples

>>> observer_gain([[0, 1], [0, 0]], [[1, 0]], [-1, -2]).round(6).tolist()
[[3.0], [2.0]]

engineeringkit.control_systems.pade_delay(delay, order=1)[source]#

Diagonal Padé approximant of a time delay \(e^{-s\theta}\).

Parameters:

delay (float) – Dead time \(\theta\) [s], positive.
order (int, optional) – Degree \(n \ge 1\) of numerator and denominator.

Returns:

\(D(-s)/D(s)\) with \(D(s) = \sum_{k=0}^{n} \frac{(2n-k)!\,n!}{(2n)!\,k!\,(n-k)!}(s\theta)^k\).

Return type:

scipy.signal.TransferFunction

Notes

The approximant is all-pass, \(|D(-j\omega)/D(j\omega)| = 1\), and matches the Taylor series of \(e^{-s\theta}\) to order \(2n\), so its phase is accurate for \(\omega\theta\) up to about \(n\). It has \(n\) right-half-plane zeros, so its step response starts with an undershoot that the true delay does not have. See Control systems: matrix contracts and numerical stability.

References

H. Padé, “Sur la représentation approchée d’une fonction par des fractions rationnelles,” Ann. Sci. École Norm. Sup. 9 (1892) 3-93.

Examples

>>> D = pade_delay(2.0)
>>> D.num.tolist(), D.den.tolist()
([-1.0, 1.0], [1.0, 1.0])

engineeringkit.control_systems.phase_portrait(acceleration, x0, v0, t)[source]#

Trajectory of \(\ddot x = f(x, \dot x)\) in the phase plane.

Parameters:

acceleration (callable) – acceleration(x, v) returning \(\ddot x\) [position units/s²]; it is also evaluated on arrays to draw the direction field.
x0 (float) – Initial position [position units] and velocity [position units/s].
v0 (float) – Initial position [position units] and velocity [position units/s].
t (array_like) – At least two strictly increasing sample times [s], fine enough for a fixed-step RK4 integration.

Returns:

Time [s], position and velocity, and the right-hand side for plotting the direction field.

Return type:

PhasePortraitResult

Notes

Fixed RK4 steps tolerate a discontinuous \(f\), such as a relay, at first-order accuracy across each switch. Time is eliminated in the plane: the slope of a trajectory is \(dv/dx = f(x, v)/v\). Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

References

A. A. Andronov and S. E. Khaikin, Theory of Oscillations (ONTI, Moscow, 1937; with A. A. Vitt, 2nd ed., 1959).

Examples

>>> r = phase_portrait(lambda x, v: -x, 1.0, 0.0, np.linspace(0, np.pi, 2001))
>>> round(r.summary()["final_x"], 6)
-1.0

engineeringkit.control_systems.poles(sys)[source]#

Poles of a SISO model, the roots of its denominator.

Parameters:: sys (scipy.signal.lti) – SISO continuous-time model.
Returns:: Poles [1/s].
Return type:: ndarray of complex

Notes

The system is asymptotically stable when every pole has a negative real part. See Control systems: matrix contracts and numerical stability.

Examples

>>> np.sort(poles(transfer_function([1], [1, 3, 2]))).real.tolist()
[-2.0, -1.0]

engineeringkit.control_systems.relay_describing_function(amplitude, output_level, hysteresis=0.0)[source]#

Describing function of a relay, ideal or with hysteresis.

Parameters:

amplitude (float or array_like) – Input amplitude \(A\) [input units], greater than hysteresis.
output_level (float) – Relay output \(\pm M\) [output units], positive.
hysteresis (float, optional) – Switching threshold \(h\) [input units], nonnegative: the relay switches up when the input rises past \(h\) and down when it falls past \(-h\).

Returns:

\(N(A) = \frac{4M}{\pi A}e^{-j\arcsin(h/A)}\); real when \(h = 0\) [output units per input unit].

Return type:

float, complex or ndarray

Notes

Hysteresis makes the relay lag the input by \(\arcsin(h/A)\). Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

References

R. J. Kochenburger, “A frequency response method for analyzing and synthesizing contactor servomechanisms,” Trans. AIEE 69 (1950) 270-284.

Examples

>>> round(float(relay_describing_function(1.0, np.pi / 4)), 12)
1.0

engineeringkit.control_systems.root_locus(num, den, gains)[source]#

Closed-loop characteristic roots of \(D(s) + K N(s)\) for each gain.

Parameters:

num (array_like) – Open-loop numerator and denominator, descending powers of \(s\). The open loop must be strictly proper so the number of roots is constant.
den (array_like) – Open-loop numerator and denominator, descending powers of \(s\). The open loop must be strictly proper so the number of roots is constant.
gains (array_like) – Nonempty vector of loop gains \(K\) [-].

Returns:

Shape (len(gains), len(den) - 1); row i holds the roots for gains[i]. Roots within a row are not tracked along branches.

Return type:

ndarray of complex

Notes

See Control systems: matrix contracts and numerical stability.

References

W. R. Evans, “Graphical analysis of control systems,” Trans. AIEE 67 (1948) 547-551.

Examples

>>> np.sort(root_locus([1], [1, 3, 2], [0.25])[0].real).tolist()
[-1.5, -1.5]

engineeringkit.control_systems.routh_first_column(coeffs)[source]#

First column of the Routh array of a real polynomial.

The number of sign changes down the column equals the number of roots in the open right half-plane.

Parameters:: coeffs (array_like) – Real coefficients in descending powers of \(s\), with a nonzero leading coefficient and degree at least one.
Returns:: First-column entries for rows \(s^n\) down to \(s^0\), normalised so the leading coefficient is one.
Return type:: ndarray

Notes

A row of zeros is replaced by the derivative of the auxiliary polynomial formed from the row above. A zero pivot is replaced by a small positive \(\varepsilon\), so this is a numerical teaching method rather than a symbolic proof for roots on or near the imaginary axis. See Control systems: matrix contracts and numerical stability.

References

E. J. Routh, A Treatise on the Stability of a Given State of Motion (Macmillan, 1877).

Examples

>>> routh_first_column([1, 3, 3, 1]).round(4).tolist()
[1.0, 3.0, 2.6667, 1.0]

engineeringkit.control_systems.routh_stable(coeffs)[source]#

Whether every root of a polynomial lies strictly in the left half-plane.

Parameters:: coeffs (array_like) – Real coefficients in descending powers of \(s\).
Returns:: True if the Routh column has no sign change and no root lies on the imaginary axis; marginal stability counts as unstable.
Return type:: bool

Notes

Routh signs alone cannot distinguish the auxiliary-row marginal case, so root locations confirm it with a relative tolerance. See Control systems: matrix contracts and numerical stability.

Examples

>>> routh_stable([1, 3, 3, 1]), routh_stable([1, 0, 1])
(True, False)

engineeringkit.control_systems.saturation_describing_function(amplitude, limit, slope=1.0)[source]#

Describing function of a saturation: linear with gain slope for inputs up to limit.

Parameters:

amplitude (float or array_like) – Input amplitude \(A\) [input units], positive.
limit (float) – Input \(a\) at which the output saturates [input units], positive; the output saturates at \(\pm k a\).
slope (float, optional) – Linear gain \(k\) [output units per input unit], positive.

Returns:

\(N = k\) for \(A \le a\), else \(\frac{2k}{\pi}\left[\arcsin\frac aA + \frac aA \sqrt{1 - (a/A)^2}\right]\).

Return type:

float or ndarray

Notes

\(N\) falls from \(k\) towards \(4ka/(\pi A)\), the relay value, as the amplitude grows. Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> float(saturation_describing_function(0.5, 1.0, 2.0))
2.0

engineeringkit.control_systems.second_order_metrics(zeta, omega_n)[source]#

Peak time, percent overshoot and 2 % settling time of an underdamped second-order step response.

Parameters:

zeta (float or array_like) – Damping ratio [-], strictly between 0 and 1.
omega_n (float or array_like) – Undamped natural frequency [rad/s], positive.

Returns:

peak_time \(\pi/\omega_d\) [s], percent_overshoot \(100\,e^{-\pi\zeta/\sqrt{1-\zeta^2}}\) [%] and settling_time_2pct \(4/(\zeta\omega_n)\) [s], with \(\omega_d = \omega_n\sqrt{1-\zeta^2}\).

Return type:

dict

Notes

Peak time and overshoot are exact; the settling time is the usual envelope estimate, not the exact 2 % crossing. See Control systems: matrix contracts and numerical stability.

References

H. L. Hazen, “Theory of servo-mechanisms,” J. Franklin Inst. 218 (1934) 279-331.

Examples

>>> m = second_order_metrics(0.5, 2.0)
>>> round(float(m["percent_overshoot"]), 2), float(m["settling_time_2pct"])
(16.3, 4.0)

engineeringkit.control_systems.sensitivity(loop)[source]#

Sensitivity function \(S = 1/(1 + L)\) of a unity-feedback loop.

Parameters:: loop (scipy.signal.lti) – SISO open-loop model \(L(s) = N(s)/D(s)\).
Returns:: \(D/(D + N)\), from an output disturbance to the output and from the reference to the error.
Return type:: scipy.signal.TransferFunction

Notes

Pole-zero cancellations are not removed. See Control systems: matrix contracts and numerical stability.

References

H. W. Bode, Network Analysis and Feedback Amplifier Design (Van Nostrand, 1945).

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> S = sensitivity(transfer_function([1], [1, 0]))
>>> S.num.tolist(), S.den.tolist()
([1.0, 0.0], [1.0, 1.0])

engineeringkit.control_systems.sensitivity_peak(loop, w=None)[source]#

Peaks of \(|S|\) and \(|T|\) and the stability margins they guarantee.

Parameters:

loop (scipy.signal.lti) – SISO open-loop model \(L(s)\) whose closed loop is stable.
w (array_like, optional) – Positive angular frequencies [rad/s]; defaults to 20000 points from \(10^{-4}\) to \(10^4\).

Returns:

Ms and Mt [-], the peaks of \(|S(j\omega)|\) and \(|T(j\omega)|\); Ms_frequency and Mt_frequency [rad/s], where they occur; gain_margin_bound \(M_s/(M_s - 1)\) [-] and phase_margin_bound_deg \(2\arcsin(1/(2M_s))\) [deg], lower bounds on the margins.

Return type:

dict

Notes

\(L(j\omega)\) stays outside the circle of radius \(1/M_s\) about \(-1\), which gives both bounds; \(M_s\) between 1.3 and 2 is the usual design range. Peaks are read on the grid, and they mean nothing for an unstable closed loop. Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> p = sensitivity_peak(transfer_function([1], [1, 0]))
>>> round(p["Ms"], 6), round(p["Mt"], 6)
(1.0, 1.0)

engineeringkit.control_systems.smith_predictor(controller, plant, delay, order=3)[source]#

Equivalent controller of a Smith predictor around a delayed plant.

The plant is \(G(s)e^{-s\theta}\). The predictor feeds back the output of a delay-free model \(G\) and corrects it by the difference between the measured output and the delayed model output. The result is the single controller

\[C_{\mathrm{SP}} = \frac{C}{1 + C G\,(1 - e^{-s\theta})},\]

which makes the closed loop \(CGe^{-s\theta}/(1 + CG)\): the delay-free design, delayed.

Parameters:

controller (scipy.signal.lti) – Primary controller \(C(s)\), designed for the delay-free plant.
plant (scipy.signal.lti) – Delay-free plant model \(G(s)\).
delay (float) – Dead time \(\theta\) [s], positive.
order (int, optional) – Padé order used for the delay inside the predictor.

Returns:

\(C_{\mathrm{SP}}(s)\) with \(e^{-s\theta}\) replaced by pade_delay() of the given order.

Return type:

scipy.signal.TransferFunction

Notes

The predictor relies on a perfect model: mismatch in \(G\) or \(\theta\) reappears in the characteristic equation. It cannot stabilise an unstable plant, whose mode stays in the open loop. In a computer the delay is a buffer of past commands, not a Padé factor. See Control systems: matrix contracts and numerical stability.

References

O. J. M. Smith, “Closer control of loops with dead time,” Chem. Eng. Prog. 53 (1957) 217-219.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> Csp = smith_predictor(transfer_function([1], [1]), transfer_function([1], [1, 1]), 1.0, order=1)
>>> Csp.num.tolist(), Csp.den.tolist()
([1.0, 3.0, 2.0], [1.0, 5.0, 2.0])

engineeringkit.control_systems.stability_margins(sys, w=None)[source]#

Worst gain and phase margins over all crossovers of a loop’s frequency response.

Parameters:

sys (scipy.signal.lti) – SISO open-loop model \(L(s)\).
w (array_like, optional) – Positive angular frequencies [rad/s]; defaults to 20000 points from \(10^{-4}\) to \(10^4\).

Returns:

gain_margin [-], the margin at the phase crossover nearest instability (inf without a phase crossover); phase_margin_deg [deg], the margin of smallest magnitude over the gain crossovers (nan without one); gain_cross_frequency and phase_cross_frequency [rad/s], the crossovers at which those worst margins are read (nan if absent); and the arrays phase_cross_frequencies [rad/s] with their gain_margins [-], and gain_cross_frequencies [rad/s] with their phase_margins_deg [deg], over every crossover in increasing frequency.

Return type:

dict

Notes

A phase crossover is any frequency where the phase passes \(-180^\circ + 360^\circ k\). Its gain margin \(1/|L(j\omega)|\) is below one when reducing the gain by that factor reaches instability, as in a conditionally stable loop; the worst gain margin is the one with the smallest \(|\log \mathrm{GM}|\). Phase margins are wrapped into \([-180^\circ, 180^\circ)\); a negative one means the plot crosses the unit circle on the far side of \(-1\), and a phase lead of that size reaches instability. With one crossover and a stable open loop that signals an unstable closed loop, but an open-loop unstable loop that needs counter-clockwise encirclements can cross on both sides, so the worst phase margin is the one of smallest magnitude. Crossings are interpolated linearly in log-frequency between grid points, so a grid that misses a crossover misses its margin. Time delays need a rational approximation such as pade_delay(). See Control systems: matrix contracts and numerical stability.

References

H. W. Bode, “Relations between attenuation and phase in feedback amplifier design,” Bell Syst. Tech. J. 19 (1940) 421-454.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> m = stability_margins(transfer_function([2], [1, 3, 2, 0]))
>>> round(m["gain_margin"], 2), round(m["phase_cross_frequency"], 3)
(3.0, 1.414)

engineeringkit.control_systems.state_feedback_poles(A, B, desired_poles)[source]#

State-feedback gain \(K\) placing the eigenvalues of \(A - BK\).

Parameters:

A (array_like) – State matrix, shape (n, n).
B (array_like) – Input matrix, shape (n, m); (A, B) must be controllable.
desired_poles (array_like) – n closed-loop poles [1/s]; complex poles must come in conjugate pairs.

Returns:

Gain \(K\), shape (m, n), for the control law \(u = -Kx\).

Return type:

ndarray

Notes

Wraps scipy.signal.place_poles(), which cannot repeat a pole more often than the rank of \(B\). With one input the gain is unique, so a repeated pole is placed by Ackermann’s formula \(K = [0\ \cdots\ 0\ 1]\,\mathcal{C}^{-1}\phi(A)\) instead, where \(\mathcal{C}\) is the controllability matrix and \(\phi\) the desired characteristic polynomial. See Control systems: matrix contracts and numerical stability.

References

W. M. Wonham, “On pole assignment in multi-input controllable linear systems,” IEEE Trans. Autom. Control 12 (1967) 660-665.

Examples

>>> state_feedback_poles([[0, 1], [0, 0]], [[0], [1]], [-1, -2]).round(6).tolist()
[[2.0, 3.0]]

engineeringkit.control_systems.state_matrices(A, B=None, C=None)[source]#

Validate and convert state-space matrices to float arrays.

Parameters:

A (array_like) – Square state matrix, shape (n, n) [1/s].
B (array_like, optional) – Input matrix, shape (n, m).
C (array_like, optional) – Output matrix, shape (p, n).

Returns:

(A, B, C), with None for matrices not supplied.

Return type:

tuple of ndarray

Raises:

ValueError – If a matrix is nonfinite, empty or has an inconsistent shape.

Notes

Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> A, B, C = state_matrices([[0, 1], [-2, -3]], B=[[0], [1]])
>>> A.shape, B.shape, C
((2, 2), (2, 1), None)

engineeringkit.control_systems.state_space(A, B, C, D)[source]#

Continuous-time state-space model \(\dot x = Ax + Bu,\ y = Cx + Du\).

Parameters:

A (array_like) – State matrix, shape (n, n) [1/s].
B (array_like) – Input matrix, shape (n, m).
C (array_like) – Output matrix, shape (p, n).
D (array_like) – Feedthrough matrix, shape (p, m).

Returns:

Continuous-time model.

Return type:

scipy.signal.StateSpace

Notes

Matrix units follow the chosen state, input and output units; see Control systems: matrix contracts and numerical stability.

Examples

>>> sys = state_space([[-2]], [[2]], [[1]], [[0]])
>>> float((-sys.C @ np.linalg.solve(sys.A, sys.B) + sys.D).item())
1.0

engineeringkit.control_systems.step_response(sys, T=None)[source]#

Response to a unit step input from rest.

Parameters:

sys (scipy.signal.lti) – Continuous-time model.
T (array_like, optional) – Time samples [s]; chosen automatically if omitted.

Returns:

(t, y): time [s] and output per unit step.

Return type:

tuple of ndarray

Notes

Wraps scipy.signal.step(). See Control systems: matrix contracts and numerical stability.

Examples

>>> t, y = step_response(transfer_function([1], [1, 1]), [0.0, 1.0])
>>> round(float(y[-1]), 4)
0.6321

engineeringkit.control_systems.system_type(loop)[source]#

Number of free integrators (poles at \(s = 0\)) in an open loop.

Parameters:: loop (scipy.signal.lti) – SISO open-loop model \(L(s)\).
Returns:: Poles at the origin minus zeros at the origin; negative when the loop differentiates.
Return type:: int

Notes

Only exact zero coefficients count, so a pole at \(10^{-9}\) is not an integrator. See Control systems: matrix contracts and numerical stability.

Examples

>>> from engineeringkit.control_systems import transfer_function
>>> system_type(transfer_function([1, 2], [1, 3, 0, 0]))
2

engineeringkit.control_systems.transfer_function(num, den)[source]#

Continuous-time transfer function \(N(s)/D(s)\).

Parameters:

num (array_like) – Numerator and denominator coefficients in descending powers of \(s\). The denominator must have degree at least zero and must not be identically zero.
den (array_like) – Numerator and denominator coefficients in descending powers of \(s\). The denominator must have degree at least zero and must not be identically zero.

Returns:

Continuous-time SISO model.

Return type:

scipy.signal.TransferFunction

Notes

Units of \(G\) are output units per input unit. See Control systems: matrix contracts and numerical stability.

Examples

>>> G = transfer_function([1], [1, 1])
>>> G.den.tolist()
[1.0, 1.0]

engineeringkit.control_systems.weight_matrix(matrix, size, name, strict=False)[source]#

Validate a symmetric positive (semi)definite weighting or covariance matrix.

Parameters:

matrix (array_like) – Candidate matrix, shape (size, size).
size (int) – Required dimension.
name (str) – Argument name used in error messages.
strict (bool, optional) – Require positive definite (True) rather than semidefinite.

Returns:

The matrix as a float array.

Return type:

ndarray

Raises:

ValueError – If the matrix has the wrong shape, is not symmetric, or is indefinite.

Notes

Inputs are raw SI numbers with no dimensional metadata; see Control systems: matrix contracts and numerical stability.

Examples

>>> weight_matrix([[2.0]], 1, "R", strict=True)
array([[2.]])

engineeringkit.control_systems.zeros(sys)[source]#

Zeros of a SISO model, the roots of its numerator.

Parameters:: sys (scipy.signal.lti) – SISO continuous-time model.
Returns:: Zeros [1/s].
Return type:: ndarray of complex

Notes

See Control systems: matrix contracts and numerical stability.

Examples

>>> zeros(transfer_function([2, 2], [1, 3, 2])).real.tolist()
[-1.0]

engineeringkit.control_systems.ziegler_nichols(ku, tu, controller='PID')[source]#

Ziegler-Nichols ultimate-sensitivity tuning.

Raise a proportional gain until the loop oscillates steadily: that gain is the ultimate gain \(K_u\) and the oscillation period is \(T_u\). The classic rules then set

type	\(K_p\)	\(T_i\)	\(T_d\)
P	\(0.5K_u\)
PI	\(0.45K_u\)	\(T_u/1.2\)
PID	\(0.6K_u\)	\(T_u/2\)	\(T_u/8\)

with \(K_i = K_p/T_i\) and \(K_d = K_p T_d\).

Parameters:

ku (float) – Ultimate proportional gain [input/error], positive.
tu (float) – Ultimate oscillation period [s], positive.
controller ({"P", "PI", "PID"}, optional) – Controller structure.

Returns:

Tuned controller.

Return type:

PID

Notes

The rules aim for roughly quarter-amplitude decay and are only a starting point; they give large overshoot on many plants. See Control systems: matrix contracts and numerical stability.

References

J. G. Ziegler and N. B. Nichols, “Optimum settings for automatic controllers,” Trans. ASME 64 (1942) 759-768.

Examples

>>> ziegler_nichols(8.0, 4.0)
PID(kp=4.8, ki=2.4, kd=2.4)

Visualizers#

Matplotlib visualizers for engineeringkit.control_systems.

engineeringkit.control_systems.visualizers.plot_bode(sys, w=None, axes=None, label=None)[source]#

Bode magnitude and phase plots on a shared logarithmic frequency axis.

Parameters:

sys (scipy.signal.lti) – Continuous-time model.
w (array_like, optional) – Angular frequencies [rad/s].
axes (sequence of two matplotlib.axes.Axes, optional) – Magnitude and phase axes; a new figure is created if omitted.
label (str, optional) – Legend label.

Returns:

(magnitude_axes, phase_axes).

Return type:

tuple of matplotlib.axes.Axes

engineeringkit.control_systems.visualizers.plot_nichols(sys, w=None, ax=None, label=None, contours=(-20.0, -12.0, -6.0, -3.0, -1.0, 0.0, 0.5, 1.0, 3.0, 6.0))[source]#

Nichols chart: open-loop gain against phase over constant closed-loop gain (M) contours.

Parameters:

sys (scipy.signal.lti) – SISO open-loop model.
w (array_like, optional) – Positive angular frequencies [rad/s], as in nichols().
ax (matplotlib.axes.Axes, optional) – Axes to draw on; a new figure is created if omitted.
label (str, optional) – Legend label.
contours (sequence of float or None, optional) – Closed-loop gains [dB] to draw as contours; None draws none.

Return type:

matplotlib.axes.Axes

engineeringkit.control_systems.visualizers.plot_nyquist(sys, w=None, ax=None, label=None)[source]#

Nyquist plot of a loop for positive and negative frequencies, marking -1.

Parameters:

sys (scipy.signal.lti) – SISO open-loop model.
w (array_like, optional) – Positive angular frequencies [rad/s]; defaults to a log grid from \(10^{-3}\) to \(10^3\).
ax (matplotlib.axes.Axes, optional) – Axes to draw on; a new figure is created if omitted.
label (str, optional) – Legend label.

Return type:

matplotlib.axes.Axes

engineeringkit.control_systems.visualizers.plot_phase_portrait(x, v, acceleration=None, ax=None, label=None)[source]#

Phase-plane trajectory, optionally over the direction field of \(\ddot x = f(x, \dot x)\).

Parameters:

x (array_like) – Matching position and velocity samples.
v (array_like) – Matching position and velocity samples.
acceleration (callable, optional) – acceleration(x, v), evaluated on arrays, to draw normalised direction arrows over the plotted region.
ax (matplotlib.axes.Axes, optional) – Axes to draw on; a new figure is created if omitted.
label (str, optional) – Legend label.

Return type:

matplotlib.axes.Axes

engineeringkit.control_systems.visualizers.plot_response(t, y, ax=None, label=None)[source]#

Plot a time response.

Parameters:

t (array_like) – Matching time samples [s] and output values.
y (array_like) – Matching time samples [s] and output values.
ax (matplotlib.axes.Axes, optional) – Axes to draw on; a new figure is created if omitted.
label (str, optional) – Legend label.

Return type:

matplotlib.axes.Axes

engineeringkit.control_systems.visualizers.plot_root_locus(num, den, gains, ax=None)[source]#

Plot closed-loop root loci with open-loop poles (x) and zeros (o).

Parameters:

num (array_like) – Strictly proper open-loop numerator and denominator.
den (array_like) – Strictly proper open-loop numerator and denominator.
gains (array_like) – Loop gains to evaluate.
ax (matplotlib.axes.Axes, optional) – Axes to draw on; a new figure is created if omitted.

Return type:

matplotlib.axes.Axes