improve robustness of hinfnorm (#1015)

* improve robustness of hinfnorm

* bump version

Author: baggepinnen

lib/ControlSystemsBase/Project.toml

@@ -2,7 +2,7 @@ name = "ControlSystemsBase"
 uuid = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 authors = ["Dept. Automatic Control, Lund University"]
 repo = "https://github.com/JuliaControl/ControlSystems.jl.git"
-version = "1.18.0"
+version = "1.18.1"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

lib/ControlSystemsBase/src/matrix_comps.jl

@@ -300,6 +300,25 @@ end
 LinearAlgebra.norm(sys::TransferFunction, p::Real=2; tol=1e-6) = norm(ss(sys), p, tol=tol)
 
 
+"""
+    sysm, T, SF = schur_form(sys)
+
+Bring `sys` to Schur form.
+
+The Schur form is characterized by `A` being Schur with the real values of eigenvalues of `A` on the main diagonal. `T` is the similarity transform applied to the system such that 
+```julia
+sysm ≈ similarity_transform(sys, T)
+```
+`SF` is the Schur-factorization of `A`.
+"""
+function schur_form(sys)
+    SF = schur(sys.A)
+    A = SF.T
+    B = SF.Z'*sys.B
+    C = sys.C*SF.Z
+    ss(A,B,C,sys.D, sys.timeevol), SF.Z, SF
+end
+
 """
     Ninf, ω_peak = hinfnorm(sys; tol=1e-6)
 
@@ -324,8 +343,8 @@ state space systems in continuous and discrete time', American Control Conferenc
 
 See also [`linfnorm`](@ref).
 """
-hinfnorm(sys::AbstractStateSpace{<:Continuous}; tol=1e-6) = _infnorm_two_steps_ct(sys, :hinf, tol)
-hinfnorm(sys::AbstractStateSpace{<:Discrete}; tol=1e-6) = _infnorm_two_steps_dt(sys, :hinf, tol)
+hinfnorm(sys::AbstractStateSpace{<:Continuous}; tol=1e-6) = _infnorm_two_steps_ct(schur_form(sys)[1], :hinf, tol)
+hinfnorm(sys::AbstractStateSpace{<:Discrete}; tol=1e-6) = _infnorm_two_steps_dt(schur_form(sys)[1], :hinf, tol)
 hinfnorm(sys::TransferFunction; tol=1e-6) = hinfnorm(ss(sys); tol=tol)
 
 """
@@ -352,10 +371,11 @@ state space systems in continuous and discrete time', American Control Conferenc
 See also [`hinfnorm`](@ref).
 """
 function linfnorm(sys::AbstractStateSpace; tol=1e-6)
-    if iscontinuous(sys)
-        return _infnorm_two_steps_ct(sys, :linf, tol)
+    sys2, _ = schur_form(sys)
+    if iscontinuous(sys2)
+        return _infnorm_two_steps_ct(sys2, :linf, tol)
     else
-        return _infnorm_two_steps_dt(sys, :linf, tol)
+        return _infnorm_two_steps_dt(sys2, :linf, tol)
     end
 end
 linfnorm(sys::TransferFunction; tol=1e-6) = linfnorm(ss(sys); tol=tol)
@@ -443,7 +463,8 @@ function _infnorm_two_steps_ct(sys::AbstractStateSpace, normtype::Symbol, tol=1e
             end
         end
     end
-    error("In _infnorm_two_steps_dt: The computation of the H∞/L∞ norm did not converge in $maxIters iterations")
+    @error("In _infnorm_two_steps_dt: The computation of the H∞/L∞ norm did not converge in $maxIters iterations")
+    return T((1+tol)*lb), T(ω_peak)
 end
 
 function _infnorm_two_steps_dt(sys::AbstractStateSpace, normtype::Symbol, tol=1e-6, maxIters=250, approxcirc=1e-8)

lib/ControlSystemsBase/test/test_matrix_comps.jl

@@ -362,5 +362,42 @@ res = observability(sys)
 @test all(==(3), res.ranks)
 @test all(<(sqrt(eps())), res.sigma_min)
 
-end 
+
+## https://github.com/JuliaControl/ControlSystems.jl/issues/1014
+P_test = zpk(
+[-101.47795511977208 + 0.0im
+    -48.91219110762173 + 0.0im
+    -7.282563985219324 + 7.114985406231401im
+    -7.282563985219324 - 7.114985406231401im
+    -7.96322290641594 + 0.0im
+    -1.5268748507837735 + 1.2594070637611725im
+    -1.5268748507837735 - 1.2594070637611725im
+    -0.7114937019357614 + 0.0im
+],
+[ -101.32273797977184 + 0.0im
+ -49.510558929948274 + 0.0im
+               -20.0 + 0.0im
+  -11.82446898287247 + 0.0im
+ -10.604444850836952 + 0.0im
+  -5.297845964509693 + 6.146324852257861im
+  -5.297845964509693 - 6.146324852257861im
+ -1.4795349979133343 + 1.2376578249023653im
+ -1.4795349979133343 - 1.2376578249023653im
+ -0.6235091399063754 + 0.0im
+ -0.2743617810110765 + 0.0im
+],
+704.6392766532747
+);
+
+C_test = zpk(
+    [-0.5 + 0.0im],
+    [0.],
+    0.6
+);
+
+S_test = sensitivity(P_test, C_test);
+n,w = hinfnorm(S_test)
+@test n ≈ 1.3056118418593037 atol=1e-3
+@test w ≈ 5.687023116875403 atol=1e-3
+end