Merge pull request #800 from JuliaControl/highorder

warn more against high-order transfer functions

Author: baggepinnen

docs/src/man/numerical.md

@@ -3,6 +3,21 @@
 ## Numerical accuracy
 Transfer functions, and indeed polynomials in general, are infamous for having poor numerical properties. Consider the simple polynomial $ax^n - 1$ which, due to rounding of the polynomial coefficients, is represented as $(a+\epsilon)x^n - 1$ where $\epsilon$ is on the order of `eps(a)`. The roots of this polynomial have a much larger $\epsilon$, due to the n:th root in the expression $\dfrac{1}{\sqrt[n]{(a + \epsilon)}}$. For this reason, it's ill-advised to use high-order transfer functions. Orders as low as 6 may already be considered high. When a transfer function is converted to a state-space representation using `ss(G)`, balancing is automatically performed in an attempt at making the numerical properties of the model better.
 
+This problem is illustrated below, where we first create a statespace system ``G`` and convert this to a transfer function ``G_1``. We then perturb a *single element* of the dynamics matrix ``A`` by adding the machine epsilon for `Float64` (`eps() = 2.22044e-16`), and convert this perturbed statespace system to a transfer function ``G_2``. The difference between the two transfer functions is enormous, the norm of the difference in their denominator coefficient vectors is on the order of ``10^{96}``.
+
+```julia
+sys = ssrand(1,1,100);
+G1 = tf(sys);
+sys.A[1,1] += eps();
+G2 = tf(sys);
+norm(denvec(G1)[] - denvec(G2)[])
+6.270683106765845e96
+```
+If we plot the poles of the two systems, they are also very different
+```julia
+scatter(poles(G1)); scatter!(poles(G2))
+```
+![Noisy poles](https://user-images.githubusercontent.com/3797491/215962177-38447944-6cca-4070-95ea-7f3829efee2e.png))
 
 
 ## Frequency-response calculation

lib/ControlSystemsBase/src/types/SisoTfTypes/SisoRational.jl

@@ -9,6 +9,7 @@ struct SisoRational{T} <: SisoTf{T}
             # The numerator is zero, make the denominator 1
             den = one(den)
         end
+        T <: AbstractFloat && length(den) > 20 && eps(T) >= eps(Float64) && @warn "High-order transfer functions are highly sensitive to numerical errors. The result may be inaccurate. Consider making use of statespace systems instead" maxlog=1
         new{T}(num, den)
     end
 end

lib/ControlSystemsBase/src/types/conversion.jl

@@ -252,6 +252,7 @@ function convert(::Type{TransferFunction{TE,SisoRational{T}}}, sys::AbstractStat
     matrix = Matrix{SisoRational{T}}(undef, size(sys))
 
     A, B, C, D = ssdata(sys)
+    T <: AbstractFloat && size(A, 1) > 20 && eps(T) >= eps(Float64) && @warn "High-order transfer functions are highly sensitive to numerical errors. The result may be inaccurate." maxlog=1
 
     # The following follows from the matrix inversion lemma:
     # det(X + uᵀv) = det(X)(1 + vᵀX⁻¹u), or