Structural stability of a 2-D system: from a polydisc to real solutions

Download: Jupyter notebook, Julia script, Maple version (same conclusions). To install PACE and run them, see Running the worksheets.

A discrete 2-D system with transfer function $N/D$ is structurally stable when $D(z_1,z_2) \neq 0$ on the closed unit bidisc $\overline{\mathbb{U}}^2$. This is a question about infinitely many points. We turn it, exactly, into the question does a system of two polynomial equations have a real solution?, and answer it with certified tools.

Running example: $D(z_1,z_2) = (2z_1^2+10z_1+12) + (z_1^2+5z_1+6)\,z_2 = (z_1+2)(z_1+3)(z_2+2),$ stable, which lets us check each step. The algorithms do not factor $D$.

using AbstractAlgebra, PACE, PaceControl

S, (z1, z2) = polynomial_ring(QQ, ["z1", "z2"])
D = (2z1^2 + 10z1 + 12) + (z1^2 + 5z1 + 6) * z2

Step 1: two 1-D tests (DeCarlo, Murray, Saeks)

\[D\]

is stable if and only if $D(z_1,1)$ and $D(1,z_2)$ have no root in the closed unit disc and $D$ has no zero on the torus $\mathbb{T}^2$.

A univariate $p$ of degree $d$ has all its roots outside the closed disc if and only if $(1+s)^d\,p\bigl(\tfrac{1-s}{1+s}\bigr)$ is Hurwitz: the map $s\mapsto\frac{1-s}{1+s}$ sends the open left half-plane onto the exterior of the disc. The exact Hurwitz test of the first worksheet then decides.

U, s = polynomial_ring(QQ, "s")

function no_root_in_disc(p)
    d = degree(p)
    h = sum(coeff(p, j) * (1 - s)^j * (1 + s)^(d - j) for j in 0:d)
    hurwitz_stable(Rational{BigInt}[coeff(h, k) for k in degree(h):-1:0])
end

p1 = AbstractAlgebra.evaluate(D, [gen(U), one(U)])      ## D(z1, 1)
p2 = AbstractAlgebra.evaluate(D, [one(U), gen(U)])      ## D(1, z2)
(p1, no_root_in_disc(p1)), (p2, no_root_in_disc(p2))

((3*s^2 + 15*s + 18, true), (12*s + 24, true))

Step 2: from the torus to real points (Möbius)

\[z = \frac{x-i}{x+i}\]

maps the real line onto $\mathbb{T}\setminus\{1\}$. Substituting it for $z_1$ and $z_2$ and clearing denominators, the numerator of $D$ is $\mathcal{R}(x_1,x_2) + i\,\mathcal{I}(x_1,x_2)$ with rational $\mathcal{R}$, $\mathcal{I}$, and $D$ has no zero on $(\mathbb{T}\setminus\{1\})^2$ if and only if $\mathcal{R}=\mathcal{I}=0$ has no real solution. (The points with a coordinate equal to $1$ are covered by step 1.)

Complex polynomials are handled as pairs (real part, imaginary part):

R, (x1, x2) = polynomial_ring(QQ, ["x1", "x2"])
cmul((a, b), (c, d)) = (a * c - b * d, a * d + b * c)
cpow(p, k) = k == 0 ? (one(R), zero(R)) : cmul(p, cpow(p, k - 1))

function moebius(D)
    d1, d2 = degree(D, 1), degree(D, 2)
    re, im = zero(R), zero(R)
    for (c, (j, k)) in zip(coefficients(D), exponent_vectors(D))
        t = cmul(cmul(cpow((x1, -one(R)), j), cpow((x1, one(R)), d1 - j)),
                 cmul(cpow((x2, -one(R)), k), cpow((x2, one(R)), d2 - k)))
        re += c * t[1]
        im += c * t[2]
    end
    re, im
end

Re, Im = moebius(D)

(36*x1^2*x2 - 10*x1 - 6*x2, 12*x1^2 + 30*x1*x2 - 2)

Step 3: real solutions of a bivariate system

\[\mathcal{R}\]

and $\mathcal{I}$ have no common factor, so the system has finitely many complex solutions. For two equations in two unknowns, LACE has a dedicated solver, which uses subresultants over the integers, computed modulo many primes, and no Gröbner basis. It returns a rational univariate representation (RUR): a polynomial rext in a new variable $T$ (here $T = x_2 - x_1$), and two rational functions of $T$ that give $x_1$ and $x_2$ at each of its roots.

gcd(Re, Im)

comps = biv_parameterization([Re, Im]; verbose = false)
comps[1].rext

5-element Vector{BigInt}:
 100
   0
 369
   0
 324

This is $324T^4 + 369T^2 + 100 = (36T^2+25)(9T^2+4)$: no real root. The certified isolation of the real solutions finds none:

The bivariate RUR of LACE (no Gröbner basis), then the certified real solutions. A constant resultant means no solution, reported as a PACEError of code -1.

function biv_real_solutions(a, b)
    comps = try
        biv_parameterization([a, b]; verbose = false)
    catch e
        e isa PACEError && e.code == -1 || rethrow()
        []
    end
    vcat([real_solutions(c.rext, c.rden, [c.num_Y, c.num_X]; prec = 30) for c in comps]...)
end

real_torus_zeros(a, b) = length(biv_real_solutions(a, b))
real_torus_zeros(Re, Im)

0

\[D\]

is structurally stable. This is a proof: every step is exact or certified.

Near the boundary: a tolerance cannot decide

Take $D_\varepsilon = 2 + z_1 + z_2 + \varepsilon$. For $\varepsilon > 0$ it is stable ($|z_1+z_2| \le 2$ on the bidisc); for $\varepsilon < 0$ it vanishes at $z_1 = z_2 = -1-\varepsilon/2$, inside the bidisc. With $\varepsilon = \pm 10^{-12}$, a floating-point test samples $|D_\varepsilon|$ on the torus:

function min_on_torus(ε, N)
    m = Inf
    for j in 0:N-1, k in 0:N-1
        w1, w2 = cis(2π * j / N), cis(2π * k / N)
        m = min(m, abs(2 + w1 + w2 + ε))
    end
    m
end

[(ε, min_on_torus(ε, 1000)) for ε in (1e-12, -1e-12)]

2-element Vector{Tuple{Float64, Float64}}:
 (1.0e-12, 1.0000000299951953e-12)
 (-1.0e-12, 1.0000000299951953e-12)

Both minima are positive and equal (the grid contains $z_1=z_2=-1$): a sampling (or any test with a tolerance above $10^{-12}$) cannot tell the stable case from the unstable one. The exact test does:

for ε in (1 // 10^12, -1 // 10^12)
    Dε = 2 + z1 + z2 + ε
    a, b = moebius(Dε)
    q1 = AbstractAlgebra.evaluate(Dε, [gen(U), one(U)])
    q2 = AbstractAlgebra.evaluate(Dε, [one(U), gen(U)])
    println("ε = ", ε, ": disc tests ", (no_root_in_disc(q1), no_root_in_disc(q2)),
            ", real zeros on the torus: ", real_torus_zeros(a, b))
end

ε = 1//1000000000000: disc tests (true, true), real zeros on the torus: 0
ε = -1//1000000000000: disc tests (true, true), real zeros on the torus: 2

For $\varepsilon = -10^{-12}$, the two real solutions are the two zeros of $D_\varepsilon$ on the torus, at distance about $10^{-6}$ from $z_1 = z_2 = -1$: the system is not stable.

Summary

Sampling the torus gives no proof, and near the boundary no tolerance separates the two answers.