H∞ norm: exact peak gain

hinf_peak computes the peak gain of a continuous-time linear system

\[\gamma = \sup_{\omega \in \mathbb{R}} \sigma_{\max}(G(i\omega)),\]

which is $\|G\|_\infty$ when $G$ is stable, exactly: the result encloses $\Gamma = \gamma^2$ between two rationals, to a relative width $2^{-64}$ by default. It is infinite when $G$ has a pole on the imaginary axis. All the input data must be exact (integers, rationals, or decimals read as the rationals they spell).

This is a second implementation of the $H_\infty$ norm, independent of the Sturm-Habicht method: it takes realisations or transfer matrices with rational entries, works on the minimal realisation, and returns $\gamma^2$ rather than $\gamma$.

Input

A state-space realisation $G(s) = C_1 (sI - A)^{-1} B_1 + D_{11}$, read by read_matrices from a file in Maple syntax (the format of the COMPleib files of test/data):

using PaceControl

M = read_matrices("test/data/AC1_matrices.txt")   # Dict with keys "A", "B1", "C1", "D11"
r = hinf_peak(M)
r = hinf_peak(M["A"], M["B1"], M["C1"], M["D11"])  # the same, from the four matrices

A transfer matrix $H(s) = N(s) / d(s)$, the entries of $N$ and $d$ being given by their coefficients, constant term first, or as Nemo polynomials:

r = hinf_peak([0, 1], [1, 1])                       # s / (s + 1)
r = hinf_peak(reshape([[0, 1], [0, -1], [0, -1], [1]], 2, 2), [1, 1])
hinf_peak([1], [0, 0, 1]).infinite                  # 1 / s^2: true

read_transfer reads a transfer matrix from a file, one polynomial per line:

# 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
Num, den = read_transfer("H.txt")
r = hinf_peak(Num, den)

Result

The result is a NormResult:

r = hinf_peak([0, 1], [1, 1])
r.infinite            # false
r.lo <= 1 <= r.hi     # true: gamma^2 = 1
r.src                 # :infini_omega, the peak is the limit omega -> infinity

r.src tells where the peak is attained: at $\omega = 0$ (:omega_0), in the limit $\omega \to \infty$ (:infini_omega, $\gamma = \sigma_{\max}(D)$) or at an interior frequency (:critique), in which case r.olo, r.ohi enclose $\omega^2$.

Method

The steps are exported by the module HinfNorm:

  1. HinfNorm.shared_parts: the transfer matrix $N(u) / (k\, d(u))$, $u = qs$, with integer coefficients, by the Leverrier algorithm, and the minimal realisation obtained by dividing out the common factors;
  2. HinfNorm.axis_verdict: a pole on the imaginary axis makes the norm infinite;
  3. HinfNorm.curve: the curve $c(\Gamma, \Omega)$, numerator of $\det(\Gamma I - G^T(-s) G(s))$ with $s^2 = -\Omega$; its points with $\Omega \ge 0$ are the squared singular values of $G(i\omega)$, $\Omega = \omega^2$;
  4. HinfNorm.boundary_candidates: the roots of $c(\Gamma, 0)$ and of the leading coefficient of $c$ in $\Omega$;
  5. HinfNorm.interior_prepare, HinfNorm.interior_boxes: the critical points $\{c = 0, \partial c / \partial \Omega = 0\}$, by explicit elimination in degree at most 2 in $\Gamma$, otherwise by the triangular decomposition of LACE;
  6. the candidates are isolated and refined until the largest one is separated from the others.
using PaceControl: HinfNorm

P0 = HinfNorm.shared_parts(M)
HinfNorm.stability(P0).stable         # asymptotic stability
cr = HinfNorm.curve(P0)
cr.degG, cr.degO                      # degrees of c in Gamma and Omega
r = hinf_peak(P0)

Options

r = hinf_peak(M; precfin = 128)       # relative width 2^-128
r = hinf_peak(M; threads = 8)         # threads of LACE for this call
r = hinf_peak(M; trace = true)        # print the steps and their times
HinfNorm.set_threads!(Sys.CPU_THREADS)