Computing the \(H_\infty\)-norm of transfer functions that do not depend on parameters

We denote the algebra of real rational functions on the imaginary axis \(i \mathbb{R}\) which are proper and have no poles on \(i \mathbb{R}\) by \(RL_\infty := \Bigg\{ \frac{n(i \omega)}{d(i \omega)} \mid n,d \in \mathbb{R}[i \omega], \ \gcd(n,d) = 1, \\ \deg_\omega(n) \leq \deg_\omega(d), V_{i \mathbb{R}}(\langle d \rangle) = \emptyset \Bigg\}.\)

For a transfer function \(G \in \mathbb{R}(s)^{u \times v}\) such that \(G_{i \mathbb{R}} \in RL^u_{\infty}\), let \(\gamma > 0\) and \(\Phi_\gamma(s) := \gamma^2 I_v - G^T(-s)G(s),\) where \(G^T\) is the transpose of \(G\).

We then can write \(\det(\Phi_\gamma(i \omega)) = \frac{n(\omega, \gamma)}{d(\omega)},\) where \(n \in \mathbb{R}[\omega, \gamma]\) and \(d \in \mathbb{R}[\omega]\) are coprime. Since \(G\) has no poles on the imaginary axis, \(d(\omega)\) does not vanish on \(\mathbb{R}\).

We denote by \(\tilde{n} \in \mathbb{R}[\omega, y]\) the square-free part of \(n\). The \(H_\infty\)-norm of \(G\) is then given by:

$$\|G\|_\infty = \max \left\{ \pi_y \left( V_{\mathbb{R}} \left( \tilde{n}, \frac{\partial \tilde{n}}{\partial \omega} \right) \right) \cup V_{\mathbb{R}} \left( \langle L_\omega(\tilde{n}) \rangle \right) \right\},$$

where \(\pi_y\) gives the \(y\)-projection.

Below, we outline the, so-called, Sturm-Habicht method to compute \(\|G\|_\infty\) as presented in Yacine Bouzidi, Alban Quadrat, Fabrice Rouillier, and Grace Younes, "Computation of the \(\mathcal{L}_{\infty}\)-norm of Finite-Dimensional Linear Systems," in Mathematical Software – ICMS 2021, Springer, 2021.

Algorithm: Sturm-Habicht Method

Input: The polynomial \(P \in (\mathbb{Z}[y])[x]\) such that the curve \(\{(x, y) \in \mathbb{R}^2 \mid P(x, y) = 0\}\) is bounded in the \(y\)-direction.

Output: Isolating interval of \(\max \left\{ \pi_y \big( V_\mathbb{R} \big( \langle P, \frac{\partial P}{\partial x} \rangle \big) \big) \cup V_\mathbb{R}(\langle L_c(P) \rangle) \right\}.\)

Steps:

  1. Compute \(\{\mathrm{Sres}_j(P, \frac{\partial P}{\partial x}, x)\}_{j=0,\ldots,\deg_x(P)}\).

  2. Compute \(y_1 < \dots < y_m\), the real roots of \(\mathrm{Sres}_0\).

  3. For \(i\) from \(1\) to \(m\):

    • If \(y_{1-i+m} \in V_\mathbb{R}(\langle L_c(P) \rangle)\), return the isolating interval of \(y_{1-i+m}\).

    • Else if \(\mathrm{Sign\_Variations}(\{\mathrm{sign}(\mathrm{stha}_d(x)(y_{1-i+m})), \dots, \mathrm{sign}(\mathrm{stha}_1(y_{1-i+m}))\}) > 0\) return the isolating interval of \(y_{1-i+m}\).

Adding and loading the necessary packages

begin
    using Pkg
    Pkg.activate();
    using  Nemo, AbstractAlgebra, PaceControl
end

Two ways to call the Hinf_StHa function:

1. With the transfer matrix \(G(s)\)

begin
    S1, (y, s) = rational_function_field(AbstractAlgebra.QQ, [:y, :s])
    G = [s//(s+1) -s//(s+1); -s//(s+1) 1//(s+1)]
end
2×2 Matrix{AbstractAlgebra.Generic.RationalFunctionFieldElem{Rational{BigInt}, AbstractAlgebra.Generic.MPoly{Rational{BigInt}}}}:
 s//(s + 1)   -s//(s + 1)
 -s//(s + 1)  1//(s + 1)
Hinf_StHa(G)
[1.59, 1.62]
Hinf_StHa(G,starting_precision=Int32(10))
[1.61792, 1.61804]

2. With the polynomial \(n(\gamma,\omega)\)

begin
    S, γ = polynomial_ring(Nemo.QQ, "γ");
    S2, ω = polynomial_ring(S, "ω");
    P = (4*γ − 1)*(4*γ + 1)*ω^4 − 4*(2*γ − 1)*(2*γ + 1) *ω^2 + 16*(γ − 1)*(γ + 1)
end
(16*γ^2 - 1)*ω^4 + (-16*γ^2 + 4)*ω^2 + 16*γ^2 - 16 
Hinf_StHa(P)
[1.09, 1.12]
Hinf_StHa(P,starting_precision=Int32(10))
[1.09436, 1.09448]