Merge pull request #836 from JuliaControl/balanceplot

apply numerical balancing before plotting

docs/src/man/numerical.md

@@ -20,6 +20,38 @@ scatter(poles(G1)); scatter!(poles(G2))
 ![Noisy poles](https://user-images.githubusercontent.com/3797491/215962177-38447944-6cca-4070-95ea-7f3829efee2e.png))
 
 
+#### State-space balancing
+
+The function [`balance_statespace`](@ref) can be used to compute a balancing transformation ``T`` that attempts to scale the system so that the row and column norms of
+```math
+\begin{bmatrix}
+TAT^{-1} & TB\\
+CT^{-1} & 0
+\end{bmatrix}
+```
+are approximately equal. This typically improves the numerical performance of several algorithms, including frequency-response calculations and continuous-time simulations. When frequency-responses are plotted using any of the built-in functions, such as [`bodeplot`](@ref) or [`nyquistplot`](@ref), this balancing is performed automatically. However, when calling [`bode`](@ref) and [`nyquist`](@ref) directly, the user is responsible for performing the balancing. The balancing is a relatively cheap operation, but it
+1. Changes the state representations of the system (but not the input-output mapping). If balancing is performed before simulation, the output will correspond to the output of the original system, but the state trajectory will not.
+2. Allocates some memory.
+
+Balancing is also automatically performed when a transfer function is converted to a statespace system using `ss(G)`, to convert without balancing, call `convert(StateSpace, G, balance=false)`.
+
+Intuitively (and simplified), balancing may be beneficial when the magnitude of the elements of the ``B`` matrix are vastly different from the magnitudes of the element of the ``C`` matrix, or when the ``A`` matrix contains very large coefficients. An example that exhibits all of these traits is the following
+```@example BALANCE
+using ControlSystemsBase, LinearAlgebra
+A = [-6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9 0.0; 0.0 -6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9; 2.0163803998106024e8 2.0163803998106024e8 -0.006223894167415392 -1.551620418759878e8 0.002358202548321148 0.002358202548321148; 0.0 0.0 5.063545034365582e-9 -0.4479539754649166 0.0 0.0; -2.824060629317756e8 2.0198389074625736e8 -0.006234569427701143 -1.5542817673286995e8 -0.7305736722226711 0.0023622473513548576; 2.0198389074625736e8 -2.824060629317756e8 -0.006234569427701143 -1.5542817673286995e8 0.0023622473513548576 -0.7305736722226711]
+B = [0.004019511633336128; 0.004019511633336128; 0.0; 0.0; 297809.51426114445; 297809.51426114445]
+C = [0.0 0.0 0.0 1.0 0.0 0.0]
+D = [0.0]
+linsys = ss(A,B,C,D)
+norm(linsys.A, Inf), norm(linsys.B, Inf), norm(linsys.C, Inf)
+```
+which after balancing becomes
+```@example BALANCE
+bsys, T = balance_statespace(linsys)
+norm(bsys.A, Inf), norm(bsys.B, Inf), norm(bsys.C, Inf)
+```
+If you plot the frequency-response of the two systems using [`bodeplot`](@ref), you'll see that they differ significantly (the balanced one is correct).
+
 ## Frequency-response calculation
 For small systems (small number of inputs, outputs and states), evaluating the frequency-response of a transfer function is reasonably accurate and very fast.
 

lib/ControlSystemsBase/src/plotting.jl

@@ -246,12 +246,13 @@ end
 ## FREQUENCY PLOTS ##
 """
     fig = bodeplot(sys, args...)
-    bodeplot(LTISystem[sys1, sys2...], args...; plotphase=true, kwargs...)
+    bodeplot(LTISystem[sys1, sys2...], args...; plotphase=true, balance = true, kwargs...)
 
 Create a Bode plot of the `LTISystem`(s). A frequency vector `w` can be
 optionally provided. To change the Magnitude scale see `setPlotScale(str)`
                                             
-If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- `balance`: Call [`balance_statespace`](@ref) on the system before plotting.
 
 `kwargs` is sent as argument to RecipesBase.plot.
 """
@@ -273,7 +274,7 @@ function _get_plotlabel(s, i, j)
     end
 end
 
-@recipe function bodeplot(p::Bodeplot; plotphase=true, ylimsphase=(), unwrap=true, hz=false)
+@recipe function bodeplot(p::Bodeplot; plotphase=true, ylimsphase=(), unwrap=true, hz=false, balance=true)
     systems, w = _processfreqplot(Val{:bode}(), p.args...)
     ws = (hz ? 1/(2π) : 1) .* w
     ny, nu = size(systems[1])
@@ -286,6 +287,9 @@ end
     link --> :x
 
     for (si,s) = enumerate(systems)
+        if balance
+            s = balance_statespace(s)[1]
+        end
         mag, phase = bode(s, w)[1:2]
         if _PlotScale == "dB" # Set by setPlotScale(str) globally
             mag = 20*log10.(mag)
@@ -378,13 +382,13 @@ optionally provided.
 - `Ms_circles`: draw circles corresponding to given levels of sensitivity (circles around -1 with  radii `1/Ms`). `Ms_circles` can be supplied as a number or a vector of numbers. A design staying outside such a circle has a phase margin of at least `2asin(1/(2Ms))` rad and a gain margin of at least `Ms/(Ms-1)`.
 - `Mt_circles`: draw circles corresponding to given levels of complementary sensitivity. `Mt_circles` can be supplied as a number or a vector of numbers.
 - `critical_point`: point on real axis to mark as critical for encirclements
-
-If `hz=true`, the hover information will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- If `hz=true`, the hover information will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- `balance`: Call [`balance_statespace`](@ref) on the system before plotting.
 
 `kwargs` is sent as argument to plot.
 """
 nyquistplot
-@recipe function nyquistplot(p::Nyquistplot; Ms_circles=Float64[], Mt_circles=Float64[], unit_circle=false, hz=false, critical_point=-1)
+@recipe function nyquistplot(p::Nyquistplot; Ms_circles=Float64[], Mt_circles=Float64[], unit_circle=false, hz=false, critical_point=-1, balance=true)
     systems, w = _processfreqplot(Val{:nyquist}(), p.args...)
     ny, nu = size(systems[1])
     nw = length(w)
@@ -394,6 +398,9 @@ nyquistplot
     θ = range(0, stop=2π, length=100)
     S, C = sin.(θ), cos.(θ)
     for (si,s) = enumerate(systems)
+        if balance
+            s = balance_statespace(s)[1]
+        end
         re_resp, im_resp = nyquist(s, w)[1:2]
         for j=1:nu
             for i=1:ny
@@ -645,17 +652,18 @@ end
 @userplot Sigmaplot
 """
     sigmaplot(sys, args...; hz=false)
-    sigmaplot(LTISystem[sys1, sys2...], args...; hz=false)
+    sigmaplot(LTISystem[sys1, sys2...], args...; hz=false, balance=true)
 
 Plot the singular values of the frequency response of the `LTISystem`(s). A
 frequency vector `w` can be optionally provided.
 
-If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- `balance`: Call [`balance_statespace`](@ref) on the system before plotting.
 
 `kwargs` is sent as argument to Plots.plot.
 """
 sigmaplot
-@recipe function sigmaplot(p::Sigmaplot; hz=false)
+@recipe function sigmaplot(p::Sigmaplot; hz=false, balance=true)
     systems, w = _processfreqplot(Val{:sigma}(), p.args...)
     ws = (hz ? 1/(2π) : 1) .* w
     ny, nu = size(systems[1])
@@ -664,6 +672,9 @@ sigmaplot
     xguide --> (hz ? "Frequency [Hz]" : "Frequency [rad/s]")
     yguide --> "Singular Values $_PlotScaleStr"
     for (si, s) in enumerate(systems)
+        if balance
+            s = balance_statespace(s)[1]
+        end
         sv = sigma(s, w)[1]'
         if _PlotScale == "dB"
             sv = 20*log10.(sv)
@@ -689,15 +700,17 @@ _to1series(y) = _to1series(1:size(y,3),y)
 
 @userplot Marginplot
 """
-    fig = marginplot(sys::LTISystem [,w::AbstractVector];  kwargs...)
-    marginplot(sys::Vector{LTISystem}, w::AbstractVector;  kwargs...)
+    fig = marginplot(sys::LTISystem [,w::AbstractVector];  balance=true, kwargs...)
+    marginplot(sys::Vector{LTISystem}, w::AbstractVector;  balance=true, kwargs...)
 
 Plot all the amplitude and phase margins of the system(s) `sys`.
-A frequency vector `w` can be optionally provided.
+
+- A frequency vector `w` can be optionally provided.
+- `balance`: Call [`balance_statespace`](@ref) on the system before plotting.
 
 `kwargs` is sent as argument to RecipesBase.plot.
 """
-@recipe function marginplot(p::Marginplot; plotphase=true)
+@recipe function marginplot(p::Marginplot; plotphase=true, hz=false, balance=true)
     systems, w = _processfreqplot(Val{:bode}(), p.args...)
     ny, nu = size(systems[1])
     s2i(i,j) = LinearIndices((nu,(plotphase ? 2 : 1)*ny))[j,i]
@@ -710,6 +723,9 @@ A frequency vector `w` can be optionally provided.
     layout --> (2ny, nu)
     label --> ""
     for (si, s) in enumerate(systems)
+        if balance
+            s = balance_statespace(s)[1]
+        end
         bmag, bphase = bode(s, w)
         for j=1:nu
             for i=1:ny
@@ -909,17 +925,18 @@ end
 @userplot Rgaplot
 """
     rgaplot(sys, args...; hz=false)
-    rgaplot(LTISystem[sys1, sys2...], args...; hz=false)
+    rgaplot(LTISystem[sys1, sys2...], args...; hz=false, balance=true)
 
 Plot the relative-gain array entries of the `LTISystem`(s). A
 frequency vector `w` can be optionally provided.
 
-If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- `balance`: Call [`balance_statespace`](@ref) on the system before plotting.
 
 `kwargs` is sent as argument to Plots.plot.
 """
 rgaplot
-@recipe function rgaplot(p::Rgaplot; hz=false)
+@recipe function rgaplot(p::Rgaplot; hz=false, balance=true)
     systems, w = _processfreqplot(Val{:sigma}(), p.args...)
     ws = (hz ? 1/(2π) : 1) .* w
     ny, nu = size(systems[1])
@@ -928,6 +945,9 @@ rgaplot
     xguide --> (hz ? "Frequency [Hz]" : "Frequency [rad/s]")
     yguide --> "Element magnitudes"
     for (si, s) in enumerate(systems)
+        if balance
+            s = balance_statespace(s)[1]
+        end
         sv = abs.(relative_gain_array(s, w))
         for j in 1:size(sv, 1)
             for i in 1:size(sv, 2)

lib/ControlSystemsBase/src/types/conversion.jl

@@ -212,6 +212,8 @@ function _balance_statespace(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMa
     return A, B, C, T
 end
 
+balance_statespace(sys, args...) = sys, I # For system types that do not have an implementation
+
 """
 `T = balance_transform{R}(A::AbstractArray, B::AbstractArray, C::AbstractArray, perm::Bool=false)`
 

lib/ControlSystemsBase/test/test_delayed_systems.jl

@@ -30,7 +30,7 @@ using ControlSystemsBase
 # Stritcly proper system
 P1 = DelayLtiSystem(ss(-1.0, 1, 1, 0))
 P1_fr = 1 ./ (im*ω .+ 1)
-@test freqresp(P1, ω)[:] == P1_fr
+@test freqresp(P1, ω)[:] ≈ P1_fr
 
 # Not stritcly proper system
 P2 = DelayLtiSystem(ss(-2.0, -1, 1, 1)) # (s+1)/(s+2)