# Stability regions in a parameter plane # # Run from this folder: maple -q 05_stability_regions.mpl # (the drawing is written next to the script, as 05_stability_regions.png). # # A 2-D discrete system has a transfer function N(z1,z2)/D(z1,z2); it is structurally stable # when its denominator D has no zero in the closed unit polydisc |z1| <= 1, |z2| <= 1. # When D also depends on two parameters u1, u2, the question becomes: for which (u1, u2) is # the system structurally stable? The answer is a union of regions of the parameter plane, # bounded by algebraic curves. We compute these curves exactly, take one rational point in # each region, and decide stability there with the non-parametric test of worksheet 3. # # Example from the literature (Li, Hou, 2021): # D = u2 z1 z2 - u1 z2 - u2 z1 + z1 z2 - z1 + 1 . # # Everything below is exact (rational arithmetic), except the drawing of section 4. with(LACE): t0 := time[real](): interface(screenwidth = 100): D0 := u2*z1*z2 - u1*z2 - u2*z1 + z1*z2 - z1 + 1: printf("D = %a\n", D0); # ## 1. The polynomial conditions # # A classical theorem (Huang 1972, refined by DeCarlo, Murray, Saeks 1977) reduces the test # on the polydisc to three tests: D(z, 1) has no zero in |z| <= 1, D(1, z) has no zero in # |z| <= 1, and D has no zero on the torus |z1| = |z2| = 1. # # As in worksheet 3, with the parameters kept symbolic: # - the torus. The Moebius map z = (x - i)/(x + i) sends the real line onto the unit circle. # After it, the numerator of D splits into a real part R and an imaginary part I, in # Q[x, y, u1, u2]; D vanishes on the torus iff {R = I = 0} has a real solution (x, y). # - the discs. The map z = (1 - x)/(1 + x) sends the closed unit disc onto the closed right # half-plane; D(z, 1) has no zero in |z| <= 1 iff the numerator D1(x) of D((1-x)/(1+x), 1) # is Hurwitz (all its roots have a negative real part). Same for D2 from D(1, z). # # The Julia worksheet calls prepare_input_stability (PaceControl); here is the same in Maple. # D is first made homogeneous in (z1, Z1) and in (z2, Z2), so that the numerators are # obtained by substituting z and Z separately. prepare_input_stability := proc(D, z1, z2, x, y) local Z1, Z2, d1, d2, H, T, RD, CD, D1, D2; d1 := degree(D, z1); d2 := degree(D, z2); H := expand(Z1^d1*Z2^d2*subs(z1 = z1/Z1, z2 = z2/Z2, D)); # the torus: z1 = (x - i)/(x + i), z2 = (y - i)/(y + i); evalc takes x, y, u real T := expand(subs(z1 = x - I, Z1 = x + I, z2 = y - I, Z2 = y + I, H)); RD := evalc(Re(T)); CD := evalc(Im(T)); # the discs: z = (1 - x)/(1 + x), the other variable set to 1 D1 := collect(expand(subs(z2 = 1, Z2 = 1, z1 = 1 - x, Z1 = 1 + x, H)), x); D2 := collect(expand(subs(z1 = 1, Z1 = 1, z2 = 1 - x, Z2 = 1 + x, H)), x); [RD, CD, D1, D2] end proc: RD, CD, D1, D2 := op(prepare_input_stability(D0, z1, z2, x, y)): printf("R = %a\nI = %a\nD1 = %a\nD2 = %a\n", RD, CD, D1, D2); # ## 2. Where can the answer change? # # On the torus. For (u1, u2) fixed, {R = I = 0} has finitely many solutions in (x, y). The # discriminant variety of this system with respect to (u1, u2) is a curve of the parameter # plane outside which the number of real solutions is locally constant (Lazard, Rouillier # 2007). Maple computes it with RootFinding:-Parametric:-DiscriminantVariety: dv := RootFinding:-Parametric:-DiscriminantVariety([RD, CD], [], [x, y]): printf("discriminant variety: %a\n", dv); # On the discs. The Lienard-Chipart criterion (1914) decides whether a real polynomial # p(x) = F(x^2) + x G(x^2), of degree n, is Hurwitz: its coefficients and the signed # principal subresultant coefficients of (x G, F) when n is odd, of (F, G) otherwise, are all # nonzero and of the sign of its leading coefficient. The signed principal subresultant # coefficients (Sturm-Habicht coefficients) are determinants built from the coefficients of # the two polynomials (Basu, Pollack, Roy, Algorithms in Real Algebraic Geometry, ch. 8). # With u symbolic, the criterion gives polynomials in u whose signs decide the question. # # The Julia worksheet uses lienard_chipart of PACE; the same procedure in Maple follows. The # coefficients of LACE:-biv_psc are the same, but it takes only one parameter. # sRes_p, ..., sRes_0 of (A, B), deg A = p > deg B = q, by the Sylvester-Habicht determinants signed_psc := proc(A, B, x) local p, q, j, k, i, M, s, eps; eps := k -> (-1)^(k*(k - 1)/2); p := degree(A, x); s := Array(0 .. p, fill = 0); s[p] := lcoeff(A, x); if B = 0 then return [seq(s[p - j], j = 0 .. p)] end if; q := degree(B, x); s[q] := eps(p - q)*lcoeff(B, x)^(p - q); for j from 0 to q - 1 do # rows: x^(q-j-1) A, ..., A, x^(p-j-1) B, ..., B; columns: x^(p+q-j-1), ..., x^j M := Matrix(p + q - 2*j, p + q - 2*j, (r, c) -> `if`(r <= q - j, coeff(expand(x^(q - j - r)*A), x, p + q - j - c), coeff(expand(x^(p - j - (r - q + j))*B), x, p + q - j - c))); s[j] := eps(p - j)*LinearAlgebra:-Determinant(M); end do; [seq(expand(s[p - j]), j = 0 .. p)] end proc: # [coefficients of p, constant term first], [signed principal subresultant coefficients] lienard_chipart := proc(p, x) local n, c, F, G, k; n := degree(p, x); c := [seq(coeff(p, x, k), k = 0 .. n)]; if n < 1 then return c, [] end if; F := add(c[2*k + 1]*x^k, k = 0 .. iquo(n, 2)); G := add(c[2*k + 2]*x^k, k = 0 .. iquo(n - 1, 2)); if type(n, odd) then c, signed_psc(expand(x*G), F, x) else c, signed_psc(F, G, x) end if end proc: lc1, sr1 := lienard_chipart(D1, x): lc2, sr2 := lienard_chipart(D2, x): lcs := [op(lc1), op(sr1), op(lc2), op(sr2)]: printf("Lienard-Chipart polynomials:\n"); for f in lcs do printf(" %a\n", f) end do; # A check with Maple's own tool: DynamicSystems:-RouthTable accepts symbolic coefficients and # returns the condition under which a polynomial is Hurwitz (the Routh-Hurwitz criterion). printf("D1 Hurwitz iff %a\n", DynamicSystems:-RouthTable(D1, x, 'stablecondition' = true)); printf("D2 Hurwitz iff %a\n", DynamicSystems:-RouthTable(D2, x, 'stablecondition' = true)); # The boundary of the regions is contained in the zero set of the irreducible factors of all # these polynomials: polys := remove(type, [op(map(op, dv)), op(lcs)], constant): boundary := [op({seq(seq(g[1], g in factors(f)[2]), f in polys)})]: # normalize each factor to a positive coefficient of u1, as in the Julia output boundary := map(g -> `if`(signum(coeff(g, u1, 1)) < 0, -g, g), boundary): boundary := [op({op(boundary)})]: # the order of the Julia output: by decreasing constant term, then increasing u2 coefficient cst := q -> eval(q, [u1 = 0, u2 = 0]): boundary := sort(boundary, (a, b) -> evalb(cst(a) > cst(b) or (cst(a) = cst(b) and coeff(a, u2, 1) < coeff(b, u2, 1)))): printf("boundary: %a\n", boundary); # Three lines. # ## 3. One point per region # # Partial CAD. A cylindrical algebraic decomposition (CAD) of the plane adapted to these # curves cuts it into cells: vertical strips u1 in (a, b) above which the curves do not cross, # each cut into the pieces between consecutive curves. Every polynomial of the list has a # constant sign on each cell, so one rational point per two-dimensional cell suffices. We # use the CAD of RegularChains (Chen, Moreno Maza), with u1 as the first variable as in Julia. with(RegularChains): with(SemiAlgebraicSetTools): Ru := PolynomialRing([u2, u1]): cad := CylindricalAlgebraicDecompose(boundary, Ru, output = list): # two-dimensional cells: odd index in u1 and in u2; their sample points are rational cells := select(c -> type(c[1][1], odd) and type(c[1][2], odd), cad): # a cell: [index, [regular chain, [interval of u1, interval of u2]]] samples := map(c -> [c[2][2][1][1], c[2][2][2][1]], cells): if not andmap(c -> c[2][2][1][1] = c[2][2][1][2] and c[2][2][2][1] = c[2][2][2][2], cells) then error "a sample point is not rational" end if; # The test at a sample point: the Lienard-Chipart signs, then the real solutions of # {R = I = 0} with the parameters replaced by their values. The bivariate solver of LACE # computes a rational univariate representation (RUR) of the system from subresultants, with # no Groebner basis, then isolates its real solutions with certified intervals. When the # resultant is a constant, there is no solution: LACE raises an error, caught here. biv_real_solutions := proc(a, b) local comps, c; try comps := LACE:-biv_parameterization([a, b], [x, y], verbose = false); catch "LACE (biv_parameterization)": if lastexception[3] <> "the resultant is a constant" then error end if; comps := []; end try; [seq(op(LACE:-real_solutions(c[1], c[2], [c[4], c[3]], prec = 30)), c in comps)] end proc: # signum, not sign: Maple's sign(0) is 1 hurwitz_at := proc(coefs, sres, p) local s, f, pt; pt := [u1 = p[1], u2 = p[2]]; s := signum(eval(coefs[-1], pt)); s <> 0 and andmap(f -> signum(eval(f, pt)) = s, [op(coefs), op(sres)]) end proc: stable_at := proc(p) local pt; pt := [u1 = p[1], u2 = p[2]]; if not (hurwitz_at(lc1, sr1, p) and hurwitz_at(lc2, sr2, p)) then return false end if; evalb(biv_real_solutions(eval(RD, pt), eval(CD, pt)) = []) end proc: verdict_cad := map(stable_at, samples): stable := [seq(`if`(verdict_cad[i], samples[i], NULL), i = 1 .. nops(samples))]: printf("CAD cells: %d, stable: %d, unstable: %d\n", nops(samples), nops(stable), nops(samples) - nops(stable)); printf("stable cells: %a\n", stable); # The Julia worksheet finds 12 cells: the projection of PACE also keeps the constant terms of # the polynomials in u2, which adds the superfluous value u1 = -3. The CAD of RegularChains # needs only u1 = -1 and u1 = 1, hence 3 strips of 3 cells. RootFinding:-Parametric has # its own decomposition, CellDecomposition, of the torus conditions only; it takes u2 as the # first variable and gives other cells, which cover the plane just as well: cd := RootFinding:-Parametric:-CellDecomposition([RD = 0, CD = 0], [x, y]): samples_cd := map(s -> eval([u1, u2], s), cd:-SamplePoints): verdict_cd := map(stable_at, samples_cd): printf("CellDecomposition: %d cells, stable: %a\n", nops(samples_cd), [seq(`if`(verdict_cd[i], samples_cd[i], NULL), i = 1 .. nops(samples_cd))]); # Connected components. A CAD cuts the regions into cylindrical cells, more than needed. # The Julia worksheet gets the connected components of the complement of the curve # prod q = 0 from its certified topology (CurveTopology), which has no equivalent in Maple. # Here the curve is made of three lines, and then the components are known exactly: each # one is an intersection of open half-planes, so it is the set of points with a given sign # vector (sign(q1), sign(q2), sign(q3)). Grouping the cells by their sign vector gives the # regions. This argument holds for lines, not for general curves. signvec := p -> map(q -> signum(eval(q, [u1 = p[1], u2 = p[2]])), boundary): regions := [op({op(map(signvec, samples))})]: printf("regions (sign vectors of the boundary lines): %d\n", nops(regions)); for sv in regions do i := ListTools:-Search(sv, map(signvec, samples)); printf(" signs %a sample %a stable: %a\n", sv, samples[i], verdict_cad[i]); end do: # The same test at the sample points of the components found by CurveTopology in Julia: K := [[-2, -3], [-2, -1], [-2, 1], [0, -1], [2, -9], [2, -2], [2, 5]]: verdict := map(stable_at, K): printf("Julia samples: %a\nverdict: %a\n", K, verdict); # Seven regions, nine CAD cells; one stable region in both cases. # ## 4. The stable region # # Its description by signs of the boundary polynomials: p := K[ListTools:-Search(true, verdict)]: for q in boundary do printf(" %a : %d\n", q, signum(eval(q, [u1 = p[1], u2 = p[2]]))) end do; # The system is stable exactly in the open triangle u1 < 1, u1 - 2 u2 > 1, u1 + 2 u2 > -3, # of vertices (1, 0), (1, -2), (-1, -1). The vertices, intersections of the lines: printf("vertices: %a\n", [seq(seq(solve({boundary[i], boundary[j]}, {u1, u2}), j = i + 1 .. 3), i = 1 .. 2)]); # The drawing: the three lines and the stable region shaded (it contains (0, -1)). # This drawing is NOT certified: plots:-implicitplot and plots:-inequal sample a grid and # draw what they see at its points; they can miss a small region or a crossing. The exact # result is the list of regions and verdicts above, not the picture. A certified drawing, # as written by CurveTopology in the Julia worksheet, has no equivalent in Maple. pl := plots:-display( plots:-inequal({boundary[1] > 0, boundary[2] > 0, boundary[3] < 0}, u1 = -4 .. 3, u2 = -4 .. 3, optionsfeasible = [color = "LightGreen"], optionsexcluded = [color = "White"]), seq(plots:-implicitplot(boundary[k], u1 = -4 .. 3, u2 = -4 .. 3, color = ["Red", "Blue", "Black"][k], thickness = 2), k = 1 .. 3), plots:-pointplot([[0, -1]], symbol = solidcircle, symbolsize = 12), plots:-textplot([[-3.5, 0.2, "u1+2u2=-3"], [2, 1.2, "u1-2u2=1"], [1.6, -3.5, "u1=1"]], align = right), title = "Regions of the (u1, u2) plane (not certified)", labels = ["u1", "u2"], axes = boxed, scaling = constrained): pl := subsindets(pl, specfunc(LEGEND), () -> NULL): # no legend plottools:-exportplot("05_stability_regions.png", pl): printf("drawing written to 05_stability_regions.png (not certified)\n"); # ## Summary # # | Step | Tool | # |-----------------------------------------------|----------------------------------------------| # | conditions on the torus and on the discs | Moebius maps, exact (procedure above) | # | where the number of real solutions can change | RootFinding:-Parametric:-DiscriminantVariety | # | where the Hurwitz conditions can change | Lienard-Chipart (procedure above), RouthTable| # | one point per region | CAD (RegularChains), sign vectors of lines | # | the answer in each region | RUR and certified isolation (LACE) | # | drawing | implicitplot, inequal: NOT certified | # # The result is exact: a finite list of regions, each certified stable or unstable. # A grid of the parameter plane would only give a picture, blind to regions smaller than its # step; the drawing above is such a picture. printf("time: %.2f s\n", time[real]() - t0); LACE:-finish():