alkahest.controls

(gain_margin, phase_margin, phase_crossover, gain_crossover).

Same order and conventions as control.margin: gain_margin is an absolute ratio (inf when the Nyquist plot never crosses the negative real axis), phase_margin in degrees (inf when |G| never crosses 1); phase_crossover is where the gain margin is measured (∠G = −180°) and gain_crossover where |G(jω)| = 1 (nan when absent). gain_margin_db is 20·log₁₀ of the ratio.

gain_margin_db

LTI system ẋ = A x + B u, y = C x + D u (MIMO).

A (n×n), B (n×m), C (p×n), D (p×m, default zero) are Matrix or nested lists; entries may be symbolic. dt marks a discrete-time system x[k+1] = A x[k] + B u[k].

nstates: int

ninputs: int

noutputs: int

is_siso() → bool

characteristic_polynomial(var)

to_tf(var)

poles(var=None) → list

zeros(var=None) → list

dc_gain()

controllability_matrix()

observability_matrix()

controllability_rank() → int

observability_rank() → int

is_controllable() → bool

is_observable() → bool

series(other)

parallel(other)

feedback(other=None, sign=-1)

impulse_response(t)

step_response(t)

discretize(T, method='zoh')

d2c(method='zoh')

minreal(var=None, method='auto')

subs(mapping)

SISO transfer function G(var) = num(var) / den(var).

numExpr | int

Numerator polynomial in var — or, when den is omitted, any rational function of var (it is then split by cancel()).

denExpr | int, optional

Denominator polynomial in var.

varExpr

The transform variable (s, or z for a discrete system).

dtExpr | float, optional

Sampling period of a discrete-time system; None for continuous.

num and den are kept as given — common factors are not cancelled behind the caller's back, because a pole–zero cancellation can hide an unstable mode. minreal() cancels them and states the hypothesis. Coefficients may carry symbolic parameters.

>>> G = TransferFunction(1, m*s**2 + b*s + k, s)
>>> G.to_ss().A        # controllable canonical form (scipy's tf2ss)
[[-b/m, -k/m], [1, 0]]

num

den

num_coeffs: list

den_coeffs: list

expr

is_discrete: bool

from_coeffs(num, den, var, *, dt=None)

degree() → tuple

is_proper() → bool

is_strictly_proper() → bool

__call__(value)

subs(mapping)

freqresp(w) → list

poles() → list

zeros() → list

dc_gain()

series(other)

parallel(other)

feedback(other=1, sign=-1)

equals(other) → bool

minreal()

to_ss()

impulse_response(t)

step_response(t)

discretize(T, method='tustin', var=None)

d2c(method='tustin', var=None)

bode(w)

margins()

A controls-layer operation refused (E-CTRL-001 … E-CTRL-006).

001 malformed input; 002 an improper transfer function where a proper one is needed; 003 a numeric routine on data with free parameters; 004 the requested object does not exist (singular Lyapunov operator, no stabilising Riccati solution, infinite DC gain); 005 out of the exact engine's reach — use the numeric route; 006 numpy is required and missing.

Bode data of a SISO system: (magnitude_dB, phase_deg) at frequencies w.

Matches scipy.signal.bode(): magnitude 20·log₁₀|G(jω)| and the unwrapped phase in degrees. Numeric: the system must have no free parameters (E-CTRL-003; use G.subs({...})). Returns numpy arrays when numpy is available, else lists.

Stabilising solution of the continuous algebraic Riccati equation.

Aᵀ X + X A − X B R⁻¹ Bᵀ X + Q = 0 (scipy's solve_continuous_are convention). The stabilising solution is the one with A − B R⁻¹ Bᵀ X Hurwitz; it is unique when it exists, and real.

method="exact" (n ≤ 2, rational or symbolic data): solves the polynomial system for the entries of a symmetric X with solve(), computes the Routh–Hurwitz conditions of each candidate's closed loop, and returns the candidate whose conditions hold at sample — a {symbol: number} point inside the parameter region of interest (default: every parameter set to a distinct value in (1, 2)). The Hurwitz conditions of the returned branch are stated in controls_side_conditions() (e.g. "… r … > 0"): the result is the stabilising solution exactly where they hold, whichever sample was used. No stabilising candidate refuses with E-CTRL-004.

method="exact" for rational data of any size (and as the fallback when the parametric system above does not solve): the stable invariant subspace of the Hamiltonian H = [[A, −BR⁻¹Bᵀ], [−Q, −Aᵀ]] is computed exactly. The characteristic polynomial of H is factored over ℚ; every irreducible factor is classified by an exact Routh test as having all its roots in the open left half plane, all in the right one, or neither. The kernels of f(H)^m for the stable factors span the subspace [U₁; U₂] and X = U₂U₁⁻¹ is rational. X is verified: symmetric with an exactly zero Riccati residual. When an irreducible factor has roots on both sides (the typical generic case: eigenvalues ±√c) the stable subspace is not defined over ℚ; "exact" refuses with E-CTRL-005 and "auto" falls back to the numeric route and states that in controls_side_conditions(). Roots on the imaginary axis, or a singular U₁, mean no stabilising solution (E-CTRL-004).

method="numeric": the Hamiltonian [[A, −BR⁻¹Bᵀ], [−Q, −Aᵀ]]'s stable invariant subspace (ordered Schur via scipy when available, else eigenvectors), X = U₂U₁⁻¹; returns a numpy.ndarray. Eigenvalues on the imaginary axis (no stabilising solution) refuse with E-CTRL-004.

"auto": exact when every entry is exact and n ≤ 2; for n > 2 exact data the exact Hamiltonian route, falling back to numeric (which needs numeric data — E-CTRL-003 if parameters remain) when the stable subspace is not rational. Floating-point data is always numeric.

>>> X = care([[0, 1], [0, 0]], [[0], [1]], [[1, 0], [0, 1]], [[r]])  # double integrator
>>> controls_side_conditions()   # r ≠ 0 and the closed-loop Hurwitz conditions (r > 0)

Hypotheses the most recent controls call on this thread rests on.

One string per condition, e.g. "m ≠ 0" (a leading coefficient that was divided by) or "(s + a) ≠ 0" (a cancelled pole–zero pair). An empty list means the last result holds for every value of its parameters. Each public entry point of controls resets the list.

Stabilising solution of the discrete algebraic Riccati equation (numeric).

Aᵀ X A − X − (Aᵀ X B)(R + Bᵀ X B)⁻¹(Bᵀ X A) + Q = 0 (scipy.linalg.solve_discrete_are()'s convention); the stabilising solution makes A − B (R + BᵀXB)⁻¹ BᵀXA Schur (all eigenvalues inside the unit disc). Returns a numpy.ndarray.

The solution is the stable invariant subspace of the symplectic matrix M = [[A + G A⁻ᵀ Q, −G A⁻ᵀ], [−A⁻ᵀ Q, A⁻ᵀ]] (G = B R⁻¹ Bᵀ, A invertible): ordered Schur via scipy when available, else eigenvectors; X = U₂U₁⁻¹. A singular A uses scipy's generalised-pencil solver (E-CTRL-006 when scipy is missing). The residual of the returned X is checked; eigenvalues of M on the unit circle mean no stabilising solution (E-CTRL-004). The data must be numeric (E-CTRL-003 if parameters remain); method="exact" is refused with E-CTRL-005 — the exact engine covers the continuous equation (care()).

Discrete LQR: (K, X) with K = (R + BᵀXB)⁻¹ BᵀXA and X from dare().

u[k] = −K x[k] minimises Σ xᵀQx + uᵀRu. Numeric (numpy.ndarray).

Solve the discrete Lyapunov (Stein) equation A X Aᵀ − X + Q = 0.

(python-control's dlyap(A, Q) and scipy.linalg.solve_discrete_lyapunov() convention.) The operator is invertible iff no two eigenvalues satisfy λᵢ λⱼ = 1; a proven singular one refuses with E-CTRL-004 (in particular an eigenvalue −1 or a reciprocal pair).

method="exact": for rational data of any size the equation is mapped by the Cayley transform F = (A − I)(A + I)⁻¹ to the continuous Lyapunov equation F X + X Fᵀ + ½(I − F) Q (I − F)ᵀ = 0, which lyapunov() solves with its structured exact solver, and the result is verified against the Stein equation (residual exactly zero). Symbolic data (up to 3×3) vectorises (I − A ⊗ A) vec X = vec Q with an exact inverse; its determinant is a side condition. E-CTRL-005 beyond. "numeric" solves the vectorised system with numpy; "auto" is exact unless an entry is a float.

>>> dlyap(Matrix([[1, 0], [0, 1]]) / 2, Matrix([[1, 0], [0, 1]]))   # (4/3)·I

Kronecker product A ⊗ B.

For A of shape (m, n) and B of shape (p, q) returns the (m·p, n·q) block matrix [[a₁₁ B, …, a₁ₙ B], …]. Entries may be symbolic. Satisfies vec(A X B) = (Bᵀ ⊗ A) vec(X) with column-major vec, which is how lyapunov() builds its linear system.

>>> kronecker(Matrix([[a, 0], [0, b]]), Matrix([[1, 1]]))
[[a, a, 0, 0], [0, 0, b, b]]

Continuous LQR: (K, X) with K = R⁻¹ Bᵀ X and X from care().

u = −K x minimises ∫ xᵀQx + uᵀRu dt. Exact or numeric as care(); side conditions as care() states them.

Solve the continuous Lyapunov equation A X + X Aᵀ + Q = 0.

(python-control's lyap(A, Q) convention; scipy.linalg.solve_continuous_lyapunov() solves A X + X Aᴴ = Q, i.e. pass -Q there.)

method="exact": the operator X ↦ A X + X Aᵀ is invertible iff no two eigenvalues of A satisfy λᵢ + λⱼ = 0; a proven singular one refuses with E-CTRL-004, a symbolic determinant is a side condition (controls_side_conditions). For rational A and Q of any size the solve is structured: Cayley–Hamilton turns the Sylvester operator into a polynomial in A and Aᵀ (O(n) matrix products over ℚ, no n²×n² system), and the result is verified — the residual is exactly zero. Symbolic A (up to 3×3) vectorises: (I ⊗ A + A ⊗ I) vec X = −vec Q with an exact inverse, then verifies the residual cancels. Other non-rational data (radicals, RootOf) up to 8×8; E-CTRL-005 beyond.

method="numeric" solves the same system with numpy and returns a numpy.ndarray. "auto" picks exact unless an entry is a float.

Matrix | numpy.ndarray

>>> lyapunov(Matrix([[0, 1], [-k, -b]]), Matrix([[1, 0], [0, 1]]))  # X for AX+XAᵀ+I=0

Gain and phase margins of a continuous SISO loop transfer function.

The crossover frequencies are roots of real polynomials in ω: Im[N(jω)·D(−jω)] = 0 (phase crossover; kept where G(jω) ≤ 0) and |N(jω)|² − |D(jω)|² = 0 (gain crossover). For rational coefficients the roots are isolated exactly (real_roots()) and refined by exact bisection; otherwise numpy's polynomial roots are used. The reported margin is the one control.margin reports: the gain margin closest to 1 in log scale and the smallest-magnitude phase margin.

Discrete time (dt set): the loop is evaluated on the unit circle z = e^{jθ}, θ = ωT ∈ [0, π] (frequencies up to the Nyquist frequency π/T). With c = cos θ both crossover conditions are real polynomials in c over ℚ — the phase condition Im[N(z)·D(z̄)] / sin θ = Σ eₖ Uₖ₋₁(c) (Chebyshev U; θ = 0, π are always on the real axis) and the gain condition |N|² − |D|² = Σ ... Tₖ(c) — and their roots in [−1, 1] are isolated exactly like the continuous ones. phase_crossover / gain_crossover are reported in rad per time unit (θ / T).

Refuses systems with free parameters (E-CTRL-003).