# The H-infinity norm: sampling, exact peak gain, and the curve behind it
#
# Maple 2026 version of the worksheet 04_hinf_norm.jl. Run it in batch, from this folder:
#
#     maple -q 04_hinf_norm.mpl
#
# It needs the LACE toolbox for Maple (with(LACE) gives 82 commands). The drawing of
# section 4 is written next to this file, as 04_hinf_norm.png.
#
# A stable linear system with transfer matrix G(s) answers an input at frequency omega
# with an amplification given by the matrix G(i omega). For a system with one input and
# one output, this amplification is the number |G(i omega)|; with several inputs and
# outputs, it is the largest singular value sigma_max(G(i omega)), the largest factor by
# which G(i omega) multiplies the length of a vector. The H-infinity norm is the largest
# amplification over all frequencies, the *peak gain*:
#
#     gamma = sup_{omega real} sigma_max(G(i omega)).
#
# Numerically, gamma is estimated by computing the amplification at a list of
# frequencies, or computed by the bisection method of Boyd, Balakrishnan and Kabamba
# (Maple's DynamicSystems:-NormHinf, Matlab's hinfnorm and getPeakGain), which needs a
# tolerance. Exactly, gamma^2 is the largest critical value of a plane algebraic curve
# (section 4), and it can be enclosed between two rationals.
#
# Two certified procedures, both written below on top of LACE, are used: HinfPeak, which
# encloses gamma^2 and tells at which frequency it is attained, and HinfStHa, which
# encloses gamma with Sturm-Habicht sequences (signed subresultants). They do in Maple
# what hinf_peak and Hinf_StHa of the Julia package PaceControl do.

with(LinearAlgebra):
with(LACE):
interface(prettyprint = 0, warnlevel = 0):
LACE:-set_verbose(false):
t_start := time[real]():

# ------------------------------------------------------------------------------------
# The procedures
#
# The squared singular values Gamma of G(i omega) are the eigenvalues of the Hermitian
# matrix G(i omega) G(i omega)^* (or G^* G, of the other size: the nonzero eigenvalues
# are the same; the smaller one is taken). For a real rational G, G(i omega)^* is
# G^T(-i omega). The determinant det(Gamma I - G(s) G^T(-s)) is even in s; with
# s^2 = -Omega, Omega = omega^2, its numerator is a polynomial c(Gamma, Omega) with
# integer coefficients: the real algebraic curve c = 0 carries the squared singular
# values as functions of Omega >= 0.

# Numerator of a rational function even in s, as a polynomial in Omega (s^2 = -Omega).
EvenToOmega := proc(f, s, Om)
    local p := expand(numer(normal(f))), k;
    if subs(s = -s, p) <> p then error "not even in %1", s end if;
    p := expand(add(coeff(p, s, 2*k) * (-Om)^k, k = 0 .. iquo(degree(p, s), 2)));
    return p / icontent(p);
end proc:

# The curve c(Gamma, Omega) of the transfer matrix G(s), primitive over the integers.
HinfCurve := proc(G::Matrix, s, Gam, Om)
    local m := RowDimension(G), n := ColumnDimension(G), Gm, H;
    Gm := map(e -> subs(s = -s, e), G);
    H := `if`(m <= n, G . Transpose(Gm), Transpose(Gm) . G);
    return EvenToOmega(Determinant(Gam * IdentityMatrix(min(m, n)) - H), s, Om);
end proc:

# Is there a pole of G on the imaginary axis? With d the least common denominator of
# the entries, d(i omega) = 0 for a real omega iff |d(i omega)|^2 = d(s) d(-s), seen as
# a polynomial e(Omega), has a root Omega >= 0. LACE:-isolate decides it exactly.
PoleOnAxis := proc(G::Matrix, s)
    local d := lcm(op(map(e -> denom(normal(e)), convert(G, list)))), Om;
    return nops(LACE:-isolate(EvenToOmega(d * subs(s = -s, d), s, Om), [0, infinity])) > 0;
end proc:

# The largest real root of the univariate f, isolated to a width below 2^(-prec):
# [lo, hi], or FAIL when f has no real root.
LargestRoot := proc(f, prec)
    local r;
    if type(f, constant) then return FAIL end if;
    r := LACE:-isolate(f, :-prec = prec);
    return `if`(r = [], FAIL, r[-1]);
end proc:

# HinfPeak(G, s): exact peak gain. The result is a record with
#   infinite  true when a pole lies on the imaginary axis (then gamma = infinity);
#   lo, hi    rationals with lo <= gamma^2 <= hi;
#   olo, ohi  an enclosure of Omega = omega^2 where gamma^2 is attained;
#   src       which kind of candidate gives the maximum (see below).
# The largest squared singular value is a continuous function of Omega on [0, infinity),
# so its supremum is taken at Omega = 0 (src = zero), as a limit at infinity, where
# G(i omega) tends to D = G(infinity) (src = infinity), on a branch Gamma = constant
# (src = constant), or at a real point of {c = 0, dc/dOmega = 0} with Omega > 0
# (src = critical). Every candidate is a squared singular value at some Omega >= 0,
# hence <= gamma^2, and gamma^2 is one of them: the largest one is gamma^2.
# The real points of the critical system are given by LACE:-real_solutions (a Rational
# Univariate Representation, then a certified isolation).
HinfPeak := proc(G::Matrix, s, {prec::posint := 128})
    local Gam, Om, c, f, cands, r, Dinf, g, c1, sols, b, best, k, p;
    if PoleOnAxis(G, s) then
        return Record('infinite' = true);
    end if;
    c := HinfCurve(G, s, Gam, Om);
    cands := [];
    r := LargestRoot(eval(c, Om = 0), prec);
    if r <> FAIL then cands := [op(cands), [r, [0, 0], zero]] end if;
    Dinf := map(e -> limit(e, s = infinity), G);
    r := LargestRoot(CharacteristicPolynomial(Transpose(Dinf) . Dinf, Gam), prec);
    if r <> FAIL then cands := [op(cands), [r, [infinity, infinity], infinity]] end if;
    # Factors of c in Gamma alone (constant singular values) and in Omega alone are
    # set apart; the other factors, without multiplicity, make c1.
    g := 1; c1 := 1;
    for f in factors(c)[2] do
        if not has(f[1], Om) then g := g * f[1];
        elif has(f[1], Gam) then c1 := c1 * f[1];
        end if;
    end do;
    r := LargestRoot(g, prec);
    if r <> FAIL then cands := [op(cands), [r, [0, infinity], constant]] end if;
    p := prec;
    do
        sols := LACE:-real_solutions([c1, diff(c1, Om)], [Gam, Om], :-prec = p);
        # a box of Omega straddling 0 (and not the point 0): refine
        if not ormap(b -> b[2][1] < 0 and b[2][2] > 0, sols) then break end if;
        p := 2 * p;
    end do;
    for b in sols do
        if b[2][1] > 0 then cands := [op(cands), [b[1], b[2], critical]] end if;
    end do;
    best := 1;
    for k from 2 to nops(cands) do
        if cands[k][1][2] > cands[best][1][2] then best := k end if;
    end do;
    return Record('infinite' = false,
                  'lo' = max(seq(b[1][1], b in cands)), 'hi' = cands[best][1][2],
                  'olo' = cands[best][2][1], 'ohi' = cands[best][2][2],
                  'src' = cands[best][3]);
end proc:

# HinfStHa(P, w, y): the Sturm-Habicht method of Bouzidi, Quadrat, Rouillier and Younes
# (2021). P(w, y) is a polynomial in w with coefficients in Q[y], whose real curve
# is bounded in the y direction. The result is an isolating interval [lo, hi] of the
# largest y above which the curve P = 0 has a real critical point in w.
# The candidates are the real roots y_1 < ... < y_m of stha_0 = Sres_0(P, dP/dw), the
# last principal Sturm-Habicht coefficient, scanned from the largest. The first y_k that
# is a root of the leading coefficient lc_w(P), or at which the principal Sturm-Habicht
# coefficients have PmV > 0 (permanences minus variations, which counts the real roots
# in w of P(w, y_k)), is returned. The signs are the exact ones of LACE:-signs_at.
HinfStHa := proc(P, w, y, {prec::posint := 80})
    local SH, f, bs, nz, polys, rts, sg, i, sv, j;
    SH := LACE:-biv_psc([P, diff(P, w)], [w, y], verbose = false);    # [stha_0, ..., stha_d]
    f := LACE:-squarefree_part(SH[1]);
    if not has(f, y) then return FAIL end if;
    bs := LACE:-isolate(f, :-prec = prec);
    nz := select(k -> SH[k] <> 0, [seq(1 .. nops(SH))]);
    polys := [lcoeff(P, w), seq(SH[k], k in nz)];
    rts, sg := op(LACE:-signs_at(f, polys, bs));
    for i from nops(rts) to 1 by -1 do
        if sg[i][1] = 0 then return rts[i] end if;
        sv := [0 $ nops(SH)];
        for j to nops(nz) do sv[nz[j]] := sg[i][j + 1] end do;
        if LACE:-pmv(ListTools:-Reverse(sv)) > 0 then return rts[i] end if;
    end do;
    return FAIL;
end proc:

# The polynomial n(omega, gamma) of HinfStHa for a transfer matrix: the curve with
# Gamma = gamma^2 and Omega = omega^2.
StHaInput := proc(G::Matrix, s, w, y)
    local Gam, Om;
    return subs(Gam = y^2, Om = w^2, HinfCurve(G, s, Gam, Om));
end proc:

# Floating point (hardware doubles), on the logarithmic grid of 200 frequencies
# 10^(-2), ..., 10^2, as a Bode plot would.
omegas := [seq(evalhf(10^(-2 + 4 * k / 199)), k = 0 .. 199)]:
SampleMax := proc(amp, ws) max(map(amp, ws)) end proc:

# ------------------------------------------------------------------------------------
# 1. A sharp resonance
#
# A lightly damped oscillator, with one input and one output:
# G(s) = 1/(s^2 + 2 zeta s + 1) with damping zeta = 10^(-3). Near omega = 1 the
# amplification has a very narrow and high peak. Its value is known in closed form:
# gamma = 1/(2 zeta sqrt(1 - zeta^2)) ~ 500.00025, at omega = sqrt(1 - 2 zeta^2).

printf("\n== 1. A sharp resonance, zeta = 1/1000 ==\n");
zeta := 1/1000:
G1 := Matrix([[1 / (s^2 + 2 * zeta * s + 1)]]):
gamma_exact := 1 / (2 * zeta * sqrt(1 - zeta^2)):
printf("closed form            gamma = %.16g\n", evalhf(gamma_exact));

# Sampling 200 frequencies, logarithmically spaced over four decades, as a Bode plot
# would: |G(i omega)| = 1/sqrt((1 - omega^2)^2 + (2 zeta omega)^2).

amp1 := proc(w) evalhf(1 / sqrt((1 - w^2)^2 + (2 * zeta * w)^2)) end proc:
printf("sampling 200 frequencies:      %.16g\n", SampleMax(amp1, omegas));

# The resonance is narrower than the spacing of the grid: the estimate is more than
# 20 times too small, and nothing in the result says so.
#
# DynamicSystems:-NormHinf is Maple's implementation of the bisection method, in
# floating point. It keeps a lower bound gmin, the largest amplification met so far, and
# tests whether g = (1 + 2 eps) gmin is above the norm: this is the case when a
# Hamiltonian matrix built from the realisation and g has no eigenvalue on the imaginary
# axis. It returns the middle of [gmin, g]. With the default relative tolerance
# eps = 10^(-6), the answer is 500.0005 at omega = 1: above the exact value by
# 5 10^(-7) in relative terms, inside the tolerance, which is all a floating-point method
# promises (the eigenvalue test itself uses a threshold, 10^(-10)). With
# eps = 10^(-10) it gets the closed form to about 15 digits.

sys1 := DynamicSystems:-TransferFunction(G1[1, 1]):
r := DynamicSystems:-NormHinf(sys1, output = [norm, peakfreq]):
printf("NormHinf (tolerance 1e-6):     %.16g at omega = %.16g\n", r[1], r[2]);
r := DynamicSystems:-NormHinf(sys1, 1e-10, output = [norm, peakfreq]):
printf("NormHinf (tolerance 1e-10):    %.16g at omega = %.16g\n", r[1], r[2]);

# HinfPeak encloses gamma^2 between two rationals; the exact value
# 1/(4 zeta^2 (1 - zeta^2)) must lie between them:

r1 := HinfPeak(G1, s):
Gam_exact := 1 / (4 * zeta^2 * (1 - zeta^2)):
printf("HinfPeak: %.16g <= gamma^2 <= %.16g, exact value inside: %a\n",
       evalf[16](r1:-lo), evalf[16](r1:-hi), evalb(r1:-lo <= Gam_exact and Gam_exact <= r1:-hi));
printf("          width of the enclosure: %.3g\n", evalf[3](r1:-hi - r1:-lo));
printf("          attained at omega^2 in [%.16g, %.16g] (%a); 1 - 2 zeta^2 = %.16g\n",
       evalf[16](r1:-olo), evalf[16](r1:-ohi), r1:-src, evalf[16](1 - 2 * zeta^2));

# The Sturm-Habicht method, on the transfer function itself, encloses gamma:

b1 := HinfStHa(StHaInput(G1, s, w, y), w, y):
printf("HinfStHa: [%s, %s]\n", convert(evalf[40](b1[1]), string),
       convert(evalf[40](b1[2]), string));
printf("          closed form inside: %a\n",
       evalb(b1[1]^2 <= Gam_exact and Gam_exact <= b1[2]^2 and b1[1] >= 0));

# Two independent certified methods, and the closed form, agree.

# ------------------------------------------------------------------------------------
# 2. A pole on the imaginary axis
#
# With zeta = 0 there is no damping: G = 1/(s^2 + 1) has poles at +-i, the amplification
# is infinite at omega = 1, and so is gamma. A sampling returns a finite number, large
# or not depending on how close a sample falls to omega = 1:

printf("\n== 2. A pole on the imaginary axis, zeta = 0 ==\n");
G0 := Matrix([[1 / (s^2 + 1)]]):
amp0 := proc(w) evalhf(abs(1 / (1 - w^2))) end proc:
printf("sampling 200 frequencies:      %.16g\n", SampleMax(amp0, omegas));

# NormHinf computes the eigenvalues of the state matrix of a realisation in floating
# point, and refuses when one of them has a real part below 10^(-10) in absolute value.
# Here the eigenvalues come out as +-1.0 i and the answer is right; but it is a
# floating-point test with a threshold, and it looks at the eigenvalues of the state
# matrix, which need not all be poles of G (see AC1 below).

try
    DynamicSystems:-NormHinf(DynamicSystems:-TransferFunction(G0[1, 1]));
catch:
    printf("NormHinf: %s\n", StringTools:-FormatMessage(lastexception[2 .. -1]));
end try;

# The exact method decides that a root of the denominator lies on the axis:

printf("HinfPeak: infinite = %a\n", HinfPeak(G0, s):-infinite);

# ------------------------------------------------------------------------------------
# 3. A model from COMPleib: AC1
#
# COMPleib (Leibfritz, 2004) is a library of benchmark models for control design. AC1 is
# an aircraft model: dx/dt = A x + B1 w, z = C1 x + D11 w, with 5 states, 3 inputs w and
# 2 outputs z, so that G(s) = C1 (sI - A)^(-1) B1 + D11 is a 2 x 3 matrix. Its data are
# decimal numbers, read as the rationals they spell (0.1 as 1/10). The matrices below
# are those of the file AC1_matrices.txt of PaceControl (test/data), which the Julia
# worksheet reads.

printf("\n== 3. COMPleib AC1 ==\n");
A := Matrix(5, 5, [
  [0, 0, 1.132, 0, -1],
  [0, -0.0538, -0.1712, 0, 0.0705],
  [0, 0, 0, 1, 0],
  [0, 0.0485, 0, -0.8556, -1.013],
  [0, -0.2909, 0, 1.0532, -0.6859]
]):
B1 := Matrix(5, 3, [
  [0.03593, 0, 0.01672],
  [0, 0.00989, 0],
  [0, -0.07548, 0],
  [0, 0, 0.05635],
  [0.00145, 0, 0.06743]
]):
C1 := Matrix(2, 5, [
  [0, 0.707106781186547, 0, 0, 0],
  [0, 0, 0.707106781186547, 0, 0]
]):
D11 := Matrix(2, 3, [
  [0, 0, 0],
  [0, 0, 0]
]):

# Every decimal as the rational it spells (convert(..., rational, exact)):
Aq, Bq, Cq, Dq := seq(map(convert, M, rational, exact), M in [A, B1, C1, D11]):
printf("sizes: A %a, B1 %a, C1 %a\n", [Dimension(Aq)], [Dimension(Bq)], [Dimension(Cq)]);

# Sampling sigma_max(G(i omega)), in hardware complex doubles:

Af, Bf, Cf, Df := seq(Matrix(map(evalhf, M), datatype = complex[8]), M in [Aq, Bq, Cq, Dq]):
sigmax := proc(w)
    local Gw;
    Digits := 16;    # else the scalar I*w is rounded to 10 digits before entering the matrix
    Gw := Cf . LinearSolve(Matrix(5, 5, shape = identity, datatype = complex[8])
                           * (I * w) - Af, Bf) + Df;
    return SingularValues(Gw)[1];
end proc:
printf("sampling 200 frequencies:      %.16g\n", SampleMax(sigmax, omegas));

# Here the sampling is close, 0.6 % below the true value, but only a certified
# computation tells how close. The exact peak gain, with the frequency where it is
# attained (G(s) is formed exactly, from the rational matrices):

GAC1 := map(normal, Cq . MatrixInverse(s * IdentityMatrix(5) - Aq) . Bq + Dq):
rA := HinfPeak(GAC1, s):
printf("HinfPeak: gamma^2 in [%.16g, %.16g]\n", evalf[16](rA:-lo), evalf[16](rA:-hi));
printf("          omega^2 in [%.16g, %.16g], where = %a\n",
       evalf[16](rA:-olo), evalf[16](rA:-ohi), rA:-src);
printf("          gamma ~ %.16g\n", evalhf(sqrt(rA:-lo)));

# The floating reference. NormHinf refuses the realisation of COMPleib: the state matrix
# A has the eigenvalue 0, on the axis. But 0 is not a pole of G: the first column of A
# and the first column of C1 are zero, so the first state is driven by the inputs but
# acts neither on the other states nor on the outputs. HinfPeak works on G itself, where
# this factor cancels. Without the first state the realisation gives the same G, and
# NormHinf answers when its data are the exact rationals; with the same data in
# floating point, its transfer matrix at s = 0 has infinite entries (a cancellation lost
# to rounding) and it fails.

shownorm := proc(label, sys)
    local r;
    try
        r := DynamicSystems:-NormHinf(sys, output = [norm, peakfreq]);
        printf("%s %.16g at omega^2 = %.16g\n", label, r[1], r[2]^2);
    catch:
        printf("%s error: %s\n", label, StringTools:-FormatMessage(lastexception[2 .. -1]));
    end try;
end proc:
shownorm("NormHinf, 5 states:           ", DynamicSystems:-StateSpace(A, B1, C1, D11));
shownorm("NormHinf, 4 states, rationals:",
         DynamicSystems:-StateSpace(Aq[2 .. 5, 2 .. 5], Bq[2 .. 5, ..], Cq[.., 2 .. 5], Dq));
shownorm("NormHinf, 4 states, floats:   ",
         DynamicSystems:-StateSpace(A[2 .. 5, 2 .. 5], B1[2 .. 5, ..], C1[.., 2 .. 5], D11));

# The peak is attained at a finite, nonzero frequency (where = critical: a critical
# point of the curve of section 4).

# ------------------------------------------------------------------------------------
# 4. The curve behind the norm
#
# The squared singular values Gamma of G(i omega), as functions of Omega = omega^2, lie
# on the real algebraic curve c(Gamma, Omega) = 0, the numerator of
# det(Gamma I - G(s) G^T(-s)) with s^2 = -Omega (G G^T(-s) is 2 x 2 here). Its degrees
# in Gamma and Omega:

printf("\n== 4. The curve c(Gamma, Omega) of AC1 ==\n");
cA := HinfCurve(GAC1, s, Gam, Om):
printf("degrees: Gamma %d, Omega %d\n", degree(cA, Gam), degree(cA, Om));

# Its coefficients are exact integers; the largest has this many digits:

printf("digits of the largest coefficient: %d\n",
       max(map(t -> length(abs(t)), [coeffs(cA, [Gam, Om])])));

# gamma^2 is the largest critical value of Gamma on the branch Omega >= 0 (or a limit at
# Omega = 0 or infinity). The curve is drawn below with Omega horizontal and Gamma
# vertical, around the peak, and the certified peak of HinfPeak is marked.
#
# This drawing is NOT certified. Maple's plots:-implicitplot evaluates c in floating
# point on a grid and connects the sign changes: it can miss a small oval, connect two
# close branches, or put a branch on the wrong side of a singular point
# (algcurves:-plot_real_curve gives no guarantee either). The Julia worksheet draws the
# same curve with CurveTopology, which first computes the certified topology of the
# curve; CurveTopology is not available in Maple. Only the marked point is certified.

pk := [(rA:-olo + rA:-ohi) / 2, (rA:-lo + rA:-hi) / 2]:
pc := plots:-implicitplot(cA, Om = 0 .. 4 * evalf(pk[1]), Gam = 0 .. 1.4 * evalf(pk[2]),
                          gridrefine = 4, crossingrefine = 4, color = "SteelBlue",
                          thickness = 2):
pp := plots:-pointplot([evalf(pk)], symbol = solidcircle, symbolsize = 14,
                       color = "Firebrick"):
plotfile := "04_hinf_norm.png":
plottools:-exportplot(plotfile, plots:-display(pc, pp, labels = ["Omega = omega^2", "Gamma"],
      title = "AC1: squared singular values Gamma against Omega (not certified)")):
printf("drawing (not certified) written to %s\n", plotfile);

# The top of the branch is the certified (Omega, Gamma) ~ (0.0334, 4.706).
printf("certified peak: (Omega, Gamma) ~ (%.4f, %.4f)\n", evalf(pk[1]), evalf(pk[2]));

# ------------------------------------------------------------------------------------
# Summary
#
# | Approach                    | Resonance zeta = 1e-3          | Pole on the axis | Guarantee |
# |-----------------------------|--------------------------------|------------------|-----------|
# | Sampling 200 frequencies    | more than 20 times too small   | finite           | none      |
# | NormHinf (bisection)        | within its tolerance           | refused (float)  | none      |
# | HinfPeak (exact peak gain)  | encloses gamma^2               | infinite         | proof     |
# | HinfStHa (Sturm-Habicht)    | encloses gamma                 |                  | proof     |
#
# The norm is a critical value of an algebraic curve: certified real solving decides
# which candidate is the maximum. The drawing of the curve is only an illustration in
# Maple: its topology is not certified.

printf("\ntotal time: %.1f s\n", time[real]() - t_start);
LACE:-finish():