Public API

hurwitz_stable(coefs::Vector) -> Bool

True if every root of the polynomial of coefficients coefs, highest degree first, has a negative real part: PACE.hurwitz_stable, the Liénard-Chipart criterion on subresultants computed by LACE, which also handles interval coefficients (Kharitonov's theorem). A root on the imaginary axis gives false.

hurwitz_stable([1, 3, 3, 1])     # (s + 1)^3: true
hurwitz_stable([1, 1, -1, -1])   # (s - 1)(s + 1)^2: false
hurwitz_stable([1, 0, 1])        # s^2 + 1: false

stability_parametric(RD, CD, D1, D2) -> (stable, unstable)

The same computation as stability_parametric(D), on its intermediate polynomials: RD, CD in $\mathbb{Q}[x, y, u_1, u_2]$ (the conditions on the torus) and D1, D2 (the conditions on the discs, univariate with coefficients in $\mathbb{Q}[u_1, u_2]$). Returns the indices, in the list of sample points of the partial CAD, of the stable and of the unstable cells.

stability_parametric(D) -> (stable, unstable)

Structural stability of the 2-D discrete system of denominator D(z1, z2, u1, u2), for the parameters (u1, u2): the system is stable when D has no zero in the closed unit polydisc $|z_1| \le 1, |z_2| \le 1$.

D is a polynomial in four variables, z1, z2, u1, u2 in this order, with integer coefficients (its ring may be over Nemo.ComplexField()). Returns the rational sample points [u1, u2] of the cells of a partial CAD of the parameter plane where the system is stable, and those where it is not.

The parameter plane is decomposed along the discriminant variety of the conditions on the torus and the Liénard-Chipart conditions on D(z, 1) and D(1, z); at one sample point per cell, the conditions on the torus are solved by the bivariate route of LACE (PACE.biv_real_solutions).

using PaceControl, Nemo
R, (z1, z2, u1, u2) = polynomial_ring(Nemo.ComplexField(), ["z1", "z2", "u1", "u2"])
D = u2*z1*z2 - u1*z2 - u2*z1 + z1*z2 - z1 + 1
stable, unstable = stability_parametric(D)    # 1 stable cell, 11 unstable ones

Hinf_StHa(P; starting_precision::Int32 = Int32(2))

An isolating interval of the largest $\gamma$ such that the curve $P(\omega, \gamma) = 0$ has a real critical point above it, by the Sturm-Habicht method; for the polynomial $n(\omega, \gamma)$ built from a transfer matrix $G$ (see the method below), this is $\|G\|_\infty = \sup_{\omega \in \mathbb{R}} \sigma_{\max}(G(i\omega))$.

P is a polynomial in $\omega$ with coefficients in $\mathbb{Q}[\gamma]$ (a Nemo polynomial over a univariate polynomial ring), whose real curve is bounded in the $\gamma$ direction. The result is an interval of MPFI at the precision starting_precision (in bits) or finer.

The candidates are the real roots $y_1 < \dots < y_m$ of the last signed subresultant $\mathrm{Sres}_0(P, \partial P / \partial \omega)$, scanned from the largest: the first $y_k$ that is a root of the leading coefficient $L_\omega(P)$, or at which the sequence of the principal Sturm-Habicht coefficients has a positive number of sign variations, is returned. The principal Sturm-Habicht coefficients are those of the pair $(P, \partial P / \partial \omega)$, computed by PACE.signed_principal_subresultant_coefs (multi-modular, LACE); their signs at the roots are the exact ones of PACE.signs_at.

using PaceControl, Nemo
S, γ = polynomial_ring(QQ, "γ")
S2, ω = polynomial_ring(S, "ω")
P = (4γ - 1)*(4γ + 1)*ω^4 - 4*(2γ - 1)*(2γ + 1)*ω^2 + 16*(γ - 1)*(γ + 1)
Hinf_StHa(P)                                  # encloses sqrt((5 + sqrt(21))/8)
Hinf_StHa(P, starting_precision = Int32(10))

Reference: Y. Bouzidi, A. Quadrat, F. Rouillier, G. Younes, Computation of the $L_\infty$-norm of finite-dimensional linear systems, Maple in Mathematics Education and Research, Springer, 2021.

Hinf_StHa(G::Matrix; starting_precision::Int32 = Int32(2))

An isolating interval of $\|G\|_\infty = \sup_{\omega \in \mathbb{R}} \sigma_{\max}(G(i\omega))$, by the Sturm-Habicht method, for a transfer matrix $G(s) \in \mathbb{Q}(s)^{u \times v}$ without pole on the imaginary axis.

The entries of G are rational functions of an AbstractAlgebra rational function field in two variables: the first one stands for $\gamma$, the second one is $s$. The numerator $n(\omega, \gamma)$ of $\det(\gamma^2 I_v - G^T(-i\omega) G(i\omega))$ is formed, then Hinf_StHa(n) is called.

using PaceControl, AbstractAlgebra
S, (y, s) = rational_function_field(AbstractAlgebra.QQ, [:y, :s])
G = [s//(s+1) -s//(s+1); -s//(s+1) 1//(s+1)]
Hinf_StHa(G)                                  # encloses (1 + sqrt(5))/2
Hinf_StHa(G, starting_precision = Int32(10))

hinf_peak(M::Dict; precfin = 64, maxrounds = 12, trace = false, maxbez = 0, threads) -> NormResult

The peak gain $\gamma = \sup_{\omega \in \mathbb{R}} \sigma_{\max}(G(i\omega))$ of $G(s) = C_1 (sI - A)^{-1} B_1 + D_{11}$, the realisation M being read by read_matrices. The result encloses $\Gamma = \gamma^2$ between two rationals, r.lo <= gamma^2 <= r.hi, or has r.infinite set when a pole lies on the imaginary axis. When G is stable (HinfNorm.stability), $\gamma = \|G\|_\infty$.

The peak is the largest Gamma on the curve c(Gamma, Omega) = 0 with Omega >= 0 (HinfNorm.curve): at Omega = 0, in the limit Omega -> infinity, or at a critical point in Omega; the candidates are isolated exactly and refined until the largest one is separated from the others.

Keywords:

• precfin: the enclosure of the peak has relative width at most 2^-precfin;
• maxrounds: maximal number of precision rounds to separate the candidates;
• trace: print each step and its time;
• maxbez: when positive, give up the interior critical points when their Bezout bound exceeds it;
• threads: number of threads of LACE for this call (default HinfNorm.get_threads()).

M = read_matrices("test/data/AC1_matrices.txt")
r = hinf_peak(M; threads = 4)
r.infinite || (sqrt(Float64(r.lo)), sqrt(Float64(r.hi)))

hinf_peak(Num, den; kw...) -> NormResult

The peak gain of a system given by its transfer matrix H(s) = Num(s) / den(s). Same keywords, same result. Num is a matrix of polynomials (vectors of coefficients, constant term first, Nemo polynomials or numbers), or a single polynomial for a one-input one-output system; see HinfNorm.shared_parts. read_transfer reads them from a file.

hinf_peak([0, 1], [1, 1])                       # s/(s+1): gamma^2 = 1, at omega -> infinity
hinf_peak(reshape([[0, 1], [0, -1], [0, -1], [1]], 2, 2), [1, 1])
hinf_peak([1], [0, 0, 1]).infinite              # 1/s^2: pole on the axis

hinf_peak(A, B1, C1, D11; kw...) -> NormResult

The peak gain of the realisation given by four matrices with exact entries (integers or rationals, not floats); hinf_peak(read_matrices(path)) reads them from a file.

hinf_peak(P0::HinfNorm.Parts; kw...) -> NormResult

The peak gain from parts already built by HinfNorm.shared_parts.

NormResult

The result of hinf_peak.

• infinite: a pole lies on the imaginary axis (src = :axe), the norm is infinite;
• lo, hi: rationals with lo <= gamma^2 <= hi, gamma the peak gain;
• src: where the peak is attained: :omega_0 (at omega = 0), :infini_omega (limit omega -> infinity, gamma = sigma_max(D)) or :critique (at an interior critical point);
• olo, ohi: enclosure of Omega = omega^2 at the peak (0 at omega = 0, negative for omega -> infinity);
• nbound, nint: number of boundary and interior candidates; prec: final precision in bits; msg: diagnostic message;
• degG, degO: degrees of the curve c(Gamma, Omega); dred: degree in Omega of the reduced curve; mode, ncomp, shear: how the interior system was solved (see HinfNorm.InterPrep);
• t: seconds spent in the axis test, the curve, the boundary, the interior system and the isolation.

read_matrices(path) -> Dict{String, Matrix{Rational{BigInt}}}

The realisation A, B1, C1, D11 read from a file written in Maple syntax, A := Matrix(r, c, [[...], ...]):, as exact rational matrices (decimals are read as the rationals they spell). $G(s) = C_1 (sI - A)^{-1} B_1 + D_{11}$.

read_transfer(path) -> (Num, den)

A transfer matrix read from a file, ready for hinf_peak(Num, den).

The format is one polynomial per line, coefficients separated by spaces, constant term first; # starts a comment and blank lines are ignored.

# H(s) = Num(s) / den(s)
den: 2 3 1              # s^2 + 3s + 2
num 1 1: 1 0 1          # s^2 + 1        (output 1, input 1)
num 1 2: 0 1            # s
num 2 1: 3/4            # a rational is written a/b
num 2 2: 1.25           # a decimal denotes the rational it spells

An entry that is not named is zero, and the size of the matrix is the largest pair of indices that appears. Coefficients are read exactly: 1.25 is 5/4, never the nearest double.

HinfNorm

Exact peak gain $\gamma = \sup_{\omega} \sigma_{\max}(G(i\omega))$ of a continuous-time linear system, its $H_\infty$ norm when the system is stable. The system is given by a realisation $G(s) = C (sI - A)^{-1} B + D$ with exact rational entries, or by a transfer matrix $G(s) = N(s) / d(s)$.

hinf_peak (exported by PaceControl) gives an enclosure of $\Gamma = \gamma^2$ between two rationals; the other exported names are its steps: minimal realisation (shared_parts), poles on the imaginary axis (axis_verdict), stability (stability), the curve $c(\Gamma, \Omega)$ of the squared singular values in $\Omega = \omega^2$ (curve), its boundary candidates (boundary_candidates) and interior critical points (interior_prepare, interior_boxes).

M = read_matrices("AC1_matrices.txt")
r = hinf_peak(M)          # r.lo <= gamma^2 <= r.hi

Parts

The transfer matrix of the minimal realisation, G(s) = Num(u) / (kq * dmin(u)) with u = q*s (see transfer_parts); deg_d is the degree of the denominator before the common factors were removed. Built once by shared_parts and passed to axis_verdict, stability, curve and hinf_peak.

shared_parts(A, B, C, D) -> Parts
shared_parts(M::Dict) -> Parts

The minimal transfer matrix of the realisation (A, B, C, D) (exact rational matrices), or of M["A"], M["B1"], M["C1"], M["D11"] as read by read_matrices.

shared_parts(Num, den) -> Parts

The system given by its transfer matrix H(s) = Num(s) / den(s) rather than by a realisation.

Num is a matrix of polynomials in s – Num[i, j] the numerator of the transfer from input j to output i – and den the denominator they share. A polynomial is a vector of coefficients with the constant term first, a polynomial of Nemo, or a number. Coefficients may be integers or rationals, exactly; a decimal read from a file must be turned into the rational it spells (dec2rat_pub) and not into the nearest double.

A single polynomial instead of a matrix is the one-input one-output case.

shared_parts([1, 0, 1], [2, 3, 1])         # (s^2 + 1) / (s^2 + 3s + 2), SISO
shared_parts([[1, 0, 1]], [2, 3, 1])       # the same, written as one row
shared_parts([[1], [0, 1]], [2, 3, 1])     # one input, two outputs: a row

A vector of numbers is a single polynomial; a vector of polynomials is one row.

The denominator is shared by the whole matrix: a system whose entries have different denominators is put over their common one first.

transfer_parts(A, B, C, D, S) -> (p, Num, kq, q)

p and the entries of Num in S = ZZ[u], u = q*s, such that, exactly,

G(s) = Num(u) / (kq * p(u)),    u = q*s

with p(u) = det(uI - qA). The constant kq does not matter for the poles, but it scales Gamma in the curve.

reduced_transfer(A, B, C, D, S) -> (p, Num, kq, q)

transfer_parts with the numerators and the denominator divided by their common gcd, which gives the minimal realisation. A non-minimal mode at s = 0 would give the curve a factor Omega and make c(Gamma, 0) identically zero.

minimal_denominator(p, Num)

p divided by the gcd of p and every entry of Num: modes that are neither controllable nor observable drop out. Returns (dmin, deg p, deg dmin).

dec2rat_pub(s) -> Rational{BigInt}

A decimal string as the exact rational it spells: "1.25" is 5//4.

axis_verdict(A, B, C, D) -> AxisResult
axis_verdict(M::Dict), axis_verdict(P0::Parts), axis_verdict(Num, den)

Whether G has a pole on the imaginary axis, hence an infinite norm.

dmin(i y) = P(y) + i Q(y): a purely imaginary root of dmin is a real root common to P and Q, hence a real root of their gcd. The scaling u = q s with q > 0 does not move the axis.

The axis verdict of a system given by its transfer matrix.

AxisResult

Result of axis_verdict: on_axis is true when a pole lies on the imaginary axis; deg_d and deg_dmin are the degrees of the denominator before and after the removal of the non-minimal modes; nroots is the number of distinct real roots y with dmin(iy) = 0 (-1 when dmin vanishes on the whole axis).

n_real_roots(f)

Number of distinct real roots of f, exactly: isolated by PACE.isolation, which takes the square-free part itself.

stability(A, B, C, D; precs = (64, 128, 256, 512), balls = true) -> StabResult
stability(M::Dict; kw...)

Asymptotic stability of G, decided on the denominator of the minimal realisation.

Tests, in this order: no pole; a pole on the imaginary axis (unstable, infinite norm); coefficients of different signs (unstable); the roots enclosed by ball arithmetic at the precisions precs (when balls); the Routh-Hurwitz criterion of PACE (routh_stable).

stability(P0::Parts; precs, balls) -> StabResult

The same verdict, from the minimal denominator P0.dmin of parts already built.

StabResult

Result of stability: stable; axis, true when a pole lies on the imaginary axis; method, what decided: :aucun_pole (no pole), :axe_imaginaire (pole on the axis), :coefficient_non_positif (coefficients of different signs), :boules (ball arithmetic, at precision prec) or :routh (integer Routh table).

routh_stable(co::Vector{ZZRingElem}) -> Bool

True if the polynomial of coefficients co, from degree 0 to degree n, is Hurwitz stable: PACE.hurwitz_stable on the coefficients highest degree first, exact over the integers. A root on the imaginary axis, or a zero leading coefficient, gives false.

curve(A, B, C, D) -> CurveResult
curve(M::Dict), curve(P0::Parts), curve(Num, den)

The curve c(Gamma, Omega): numerator of det(Gamma I - G~(s) G(s)), with G~(s) = G(-s)^T, after s^2 = -Omega, formed on the smaller side; its real points with Omega >= 0 are the squared singular values Gamma of G(i omega), Omega = omega^2.

M = N~ N (or N N~) and e = d(u) d(-u) are polynomial, hence so is det(Gamma e I - M) = chi_M(Gamma e): no rational function appears. The curve is formed on the minimal realisation (see reduced_transfer), divided by its content in ZZ and by its content in Gamma (strip_content_G).

The curve of a system given by its transfer matrix. shared_parts(Num, den) says what the two arguments may be.

CurveResult

Result of curve: c in ZZ[Gamma, Omega], its degrees degG and degO, the side side on which the Gram matrix was formed (:z for G G~, :w for G~ G) and its size r.

eval_curve(cr, Om, Gam) -> (value, scale)

The curve evaluated exactly at the rationals (Omega, Gamma), in QQ, with scale the sum of the moduli of its terms: |value| / scale measures how close the point is to the curve, including at a multiple root.

strip_content_G(c, T)

c divided by its content seen as a polynomial in Gamma, that is, by its largest factor depending on Omega alone. Such a factor comes from e(Omega) = |d(i omega)|^2 (multiple poles), not from a singular value.

coeff_in(p, ivar, k, S)

The coefficient of x_ivar^k in the bivariate p, as a polynomial of S in the other variable.

coeff_O(p, k, S)

The coefficient of Omega^k in p(Gamma, Omega), a polynomial of S in Gamma.

embed(u, T, ivar)

The univariate u as a polynomial of the bivariate ring T in its variable ivar.

Candidate

A positive real root, PACE.Root refined in place by LACE, and the family src it comes from (:omega_0, :infini_omega).

isolate_pos(f, src; prec = 64) -> Vector{Candidate}

The positive real roots of f, increasing, of width below 2^-prec, tagged src. The root 0 is discarded: it is a maximum only for the zero system, handled by the caller.

largest_root(cs; maxrounds = 24, target = 10^-40) -> Candidate

The largest of the roots cs, refining (PACE.refine!) the leader and the candidates that overlap it until they separate. Candidates that still overlap once all are narrower than the relative width target are taken as equal.

boundary_candidates(cr::CurveResult) -> BoundResult

The candidates of the univariate families: positive roots of c(Gamma, 0) (tagged :omega_0) and of the leading coefficient in Omega (limit Omega -> infinity, tagged :infini_omega). The content ccst of c in Omega divides c(Gamma, 0), so its roots are already candidates.

BoundResult

Result of boundary_candidates: the candidates cands and the polynomials in Gamma c0 = c(Gamma, 0), cinf (leading coefficient in Omega) and ccst (content in Omega).

peak_boundary(M::Dict; cr = nothing) -> PeakResult

The largest boundary candidate of the curve cr (by default curve(M)), enclosed in Gamma = gamma^2: a lower bound on the peak, which hinf_peak compares with the interior critical points.

PeakResult

Result of peak_boundary: [lo, hi] encloses the largest boundary candidate, of family src (:nul when there is none), among ncand candidates.

interior_reduce(cr::CurveResult; trace = false)

The reduced curve c2 whose critical points in Omega are computed, or nothing when there is no interior system. Three successive divisions:

gcd(c, dc/dOmega) repeated factors; no extremum is lost content in Gamma factors free of Gamma, which are not singular values but degeneracies of e content in Omega factors free of Omega, which would give a one-dimensional component; their roots are already among the boundary candidates

interior_prepare(cr::CurveResult; trace = false, maxbez = 0) -> InterPrep

A parametrisation of the critical points of the reduced curve {c2 = 0, dc2/dOmega = 0} (interior_reduce).

In degree 1 and 2 in Gamma the elimination is explicit (modes :deg1, :deg2). Otherwise the triangular decomposition biv_decomposition of LACE is used (:decomp): each generic component reads A(Omega) Gamma + B(Omega) = 0 above the roots of a polynomial in Omega, whence Gamma = -B/A; a non-generic one of degree one is solved by rational_fibre. The sheared RUR (biv_parameterization, mode :rur) is used only for the other non-generic components: the curve has zeros at infinity in both directions, which the decomposition puts aside. With maxbez > 0, the computation is given up when the Bezout bound exceeds maxbez.

InterPrep

Result of interior_prepare: ok; mode (:deg1, :deg2, :decomp, :rur, or :aucun on failure); the components comps, each giving Gamma and Omega as rational functions at the roots of a univariate eliminant; the shear coefficient a (0 without shear); the total degree nsol of the eliminants; the degree degO of the reduced curve in Omega; and a message msg.

interior_boxes(pr::InterPrep, prec) -> (boxes, ndrop)

The critical points with Gamma > 0 and Omega > 0, enclosed in Gamma and Omega, and the number ndrop of real roots of the eliminants where the parametrisation reads 0/0. The real roots of each eliminant are isolated by LACE, and Gamma and Omega evaluated there with prec significant bits (PACE.rur_coords). All real roots are kept: under a shear the parametrising variable is not Omega.

set_threads!(n) -> Int

Set the number of threads LACE may use from here on (1, the default, is sequential). Returns the previous value.

set_threads!(Sys.CPU_THREADS)

get_threads() -> Int

The number of threads LACE may use.

biv_time() -> (decomposition, rur)

Seconds the last hinf_peak spent in the bivariate decomposition of LACE, and in the sheared RUR it falls back on.