JuliaControl · baggepinnen · Sep 18, 2026 · Sep 18, 2026
diff --git a/docs/src/man/numerical.md b/docs/src/man/numerical.md
@@ -30,9 +30,11 @@ 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
+are approximately equal. This typically improves the numerical performance of several algorithms, including frequency-response calculations and continuous-time simulations.
+
+Balancing is performed automatically by the functions that depend on the system only through its input-output map. This includes [`freqresp`](@ref), [`bode`](@ref), [`nyquist`](@ref), [`sigma`](@ref), [`margin`](@ref), [`delaymargin`](@ref), [`relative_gain_array`](@ref), [`dcgain`](@ref), [`hinfnorm`](@ref), [`linfnorm`](@ref), `norm`, as well as the corresponding plot functions such as [`bodeplot`](@ref) and [`nyquistplot`](@ref). All of these accept the keyword argument `balance = false` to turn the balancing off. 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.
+2. Allocates some memory. This matters mostly for the in-place functions [`freqresp!`](@ref) and [`bodemag!`](@ref), which are otherwise free of allocations that scale with the state dimension. Pass `balance = false` to these when the realization is known to be well scaled.
 
 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)`.
 
@@ -51,7 +53,7 @@ which after balancing becomes
 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).
+If you plot the frequency-response of the two systems using `bodeplot(sys, balance=false)`, you'll see that they differ significantly (the balanced one is correct). With the default `balance = true`, both plots coincide since the balancing is then performed internally.
 
 ## Extended precision and exotic number types
 Most functions in ControlSystems.jl can operate on numbers of any type, for example, computations can be performed with increased precision using the `BigFloat` type. While the `BigFloat` type is part of julia base, some functionality such as matrix factorizations used in `c2d` and `minreal` etc. require the user to load external packages to work with `BigFloat` numbers. The list below indicates how to make a number of functions work with `BigFloat` (and other exotic number types):

diff --git a/lib/ControlSystemsBase/Project.toml b/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.21.0"
+version = "1.22.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

diff --git a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
@@ -14,7 +14,7 @@ using ControlSystemsBase: downsample, _processfreqplot, _default_freq_vector,
                           sisomargin, relative_gain_array, rlocus,
                           input_names, output_names, state_names, system_name,
                           iscontinuous, isdiscrete, issiso, isrational,
-                          integrator_excess, balance_statespace, LTISystem,
+                          integrator_excess, LTISystem,
                           _PlotScale, _PlotScaleFunc, _PlotScaleStr, _span  # Use existing plot scale settings
 
 # Helper function to get y-scale transform for Makie
@@ -166,14 +166,12 @@ function CSMakie.bodeplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTISy
 
     # Plot data for each system
     for (si, s) in enumerate(systems)
-        sbal = balance ? balance_statespace(s)[1] : s
-
         intexcess = 0
-        if plotphase && adjust_phase_start && isrational(sbal)
-            intexcess = integrator_excess(sbal)
+        if plotphase && adjust_phase_start && isrational(s)
+            intexcess = integrator_excess(s)
         end
 
-        mag, phase = bode(sbal, w; unwrap=false)
+        mag, phase = bode(s, w; unwrap=false, balance)
 
         if ControlSystemsBase._PlotScale == "dB"
             mag = 20*log10.(mag)
@@ -203,7 +201,7 @@ function CSMakie.bodeplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTISy
                 if plotphase
                     phasedata = vec(phase[i, j, :])
 
-                    if adjust_phase_start && isrational(sbal) && intexcess != 0
+                    if adjust_phase_start && isrational(s) && intexcess != 0
                         nineties = round(Int, phasedata[1] / 90)
                         phasedata .+= ((90*(-intexcess-nineties)) ÷ 360) * 360
                     end
@@ -286,8 +284,7 @@ function CSMakie.nyquistplot!(fig, systems::Union{LTISystem, AbstractVector{<:LT
     θ = range(0, 2π, length=100)
 
     for (si, s) in enumerate(systems)
-        sbal = balance ? balance_statespace(s)[1] : s
-        re_resp, im_resp = nyquist(sbal, w)[1:2]
+        re_resp, im_resp = nyquist(s, w; balance)[1:2]
 
         for j in 1:nu
             for i in 1:ny
@@ -361,8 +358,7 @@ function CSMakie.sigmaplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTIS
              ylabel = "Singular Values $(ControlSystemsBase._PlotScaleStr)")
 
     for (si, s) in enumerate(systems)
-        sbal = balance ? balance_statespace(s)[1] : s
-        sv = sigma(sbal, w)[1]'
+        sv = sigma(s, w; balance)[1]'
 
         if extrema && size(sv, 2) > 2
             sv = sv[:, [1, end]]
@@ -426,15 +422,15 @@ function CSMakie.marginplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTI
 
     # Plot data for each system
     for (si, s) in enumerate(systems)
-        sbal = balance ? balance_statespace(s)[1] : s
-        bmag, bphase = bode(sbal, w)
+        bmag, bphase = bode(s, w; balance)
 
         for j in 1:nu
             for i in 1:ny
-                wgm, gm, wpm, pm, fullPhase, phasedata = sisomargin(sbal[i,j], w; 
+                wgm, gm, wpm, pm, fullPhase, phasedata = sisomargin(s[i,j], w; 
                                                                      full=true, 
                                                                      allMargins=true, 
-                                                                     adjust_phase_start)
+                                                                     adjust_phase_start,
+                                                                     balance)
 
                 # Magnitude plot
                 ax_mag = axes_mag[i, j]
@@ -588,8 +584,7 @@ function CSMakie.rgaplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTISys
              ylabel = "Element magnitudes")
 
     for (si, s) in enumerate(systems)
-        sbal = balance ? balance_statespace(s)[1] : s
-        rga = abs.(relative_gain_array(sbal, w))
+        rga = abs.(relative_gain_array(s, w; balance))
 
         for j in 1:size(rga, 1)
             for i in 1:size(rga, 2)

diff --git a/lib/ControlSystemsBase/src/analysis.jl b/lib/ControlSystemsBase/src/analysis.jl
@@ -150,16 +150,19 @@ function det(sys::Matrix{S}) where {S<:SisoZpk}
 end
 
 """
-    dcgain(sys, ϵ=0)
+    dcgain(sys, ϵ=0; balance=true)
 
 Compute the dcgain of system `sys`.
 
 equal to G(0) for continuous-time systems and G(1) for discrete-time systems.
 
 `ϵ` can be provided to evaluate the dcgain with a small perturbation into
 the stability region of the complex plane.
+
+`balance`: Call [`balance_statespace`](@ref) on the system before evaluating the transfer function, see [`freqresp`](@ref).
 """
-function dcgain(sys::LTISystem, ϵ=0)
+function dcgain(sys::LTISystem, ϵ=0; balance=true)
+    sys = _balance(sys, balance)
     return iscontinuous(sys) ? evalfr(sys, -ϵ) : evalfr(sys, exp(-ϵ*sys.Ts))
 end
 dcgain(G::Union{UniformScaling, Number, AbstractMatrix}) = G
@@ -413,8 +416,8 @@ end
 # end
 
 """
-    relative_gain_array(G, w::AbstractVector)
-    relative_gain_array(G, w::Number)
+    relative_gain_array(G, w::AbstractVector; balance=true)
+    relative_gain_array(G, w::Number; balance=true)
 
 Calculate the relative gain array of `G` at frequencies `w`. 
 G(iω) .* pinv(tranpose(G(iω)))
@@ -439,11 +442,11 @@ Reference: "On the Relative Gain Array (RGA) with Singular and Rectangular Matri
 Jeffrey Uhlmann
 https://arxiv.org/pdf/1805.10312.pdf
 """
-function relative_gain_array(G, w::AbstractVector)
-    mapslices(relative_gain_array, freqresp(G, w), dims=(1,2))
+function relative_gain_array(G, w::AbstractVector; balance=true)
+    mapslices(relative_gain_array, freqresp(G, w; balance), dims=(1,2))
 end
 
-relative_gain_array(G, w::Number) = relative_gain_array(freqresp(G, w))
+relative_gain_array(G, w::Number; balance=true) = relative_gain_array(freqresp(G, w; balance))
 
 """
     relative_gain_array(A::AbstractMatrix; tol = 1.0e-15)
@@ -490,17 +493,18 @@ function relative_gain_array(A::AbstractMatrix; tol = 1e-15)
 end
 
 """
-    wgm, gm, wpm, pm = margin(sys::LTISystem, w::Vector; full=false, allMargins=false, adjust_phase_start=true)
+    wgm, gm, wpm, pm = margin(sys::LTISystem, w::Vector; full=false, allMargins=false, adjust_phase_start=true, balance=true)
 
 returns frequencies for gain margins, gain margins (magnitude), frequencies for phase margins, phase margins (degrees)
 
 - If `!allMargins`, return only the smallest margin
 - If `full` return also `fullPhase`
 - `adjust_phase_start`: If true, the phase will be adjusted so that it starts at -90*intexcess degrees, where `intexcess` is the integrator excess of the system.
+- `balance`: Call [`balance_statespace`](@ref) on each SISO channel before computing the frequency response, see [`freqresp`](@ref).
 
 See also [`delaymargin`](@ref) and [`RobustAndOptimalControl.diskmargin`](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/api/#RobustAndOptimalControl.diskmargin)
 """
-function margin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMargins=false, adjust_phase_start=true)
+function margin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMargins=false, adjust_phase_start=true, balance=true)
     ny, nu = size(sys)
 
     T = float(numeric_type(sys))
@@ -519,7 +523,7 @@ function margin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMargin
     end
     for j=1:nu
         for i=1:ny
-            wgm[i,j], gm[i,j], wpm[i,j], pm[i,j], fullPhase[i,j] = sisomargin(sys[i,j], w; full=true, allMargins, adjust_phase_start)
+            wgm[i,j], gm[i,j], wpm[i,j], pm[i,j], fullPhase[i,j] = sisomargin(sys[i,j], w; full=true, allMargins, adjust_phase_start, balance)
         end
     end
     if full
@@ -530,21 +534,24 @@ function margin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMargin
 end
 
 """
-    ωgm, gm, ωpm, pm = sisomargin(sys::LTISystem, w::Vector; full=false, allMargins=false, adjust_phase_start=true))
+    ωgm, gm, ωpm, pm = sisomargin(sys::LTISystem, w::Vector; full=false, allMargins=false, adjust_phase_start=true, balance=true))
 
 Return frequencies for gain margins, gain margins, frequencies for phase margins, phase margins. If `allMargins=false`, only the smallest margins are returned.
+
+`balance`: Call [`balance_statespace`](@ref) on the system before computing the frequency response, see [`freqresp`](@ref).
 """
-function sisomargin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMargins=false, adjust_phase_start=true)
+function sisomargin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMargins=false, adjust_phase_start=true, balance=true)
     ny, nu = size(sys)
     if ny !=1 || nu != 1
         error("System must be SISO, use `margin` instead")
     end
-    mag, phase, w = bode(sys, w)
+    sys = _balance(sys, balance) # All frequency-response evaluations below use the balanced realization
+    mag, phase, w = bode(sys, w; balance=false)
     wgm, = _allPhaseCrossings(w, phase)
     gm = similar(wgm)
     remove = Int[]
     for i = eachindex(wgm)
-        Giw = freqresp(sys,wgm[i])[1]
+        Giw = freqresp(sys, wgm[i]; balance=false)[1]
         if sign(w[1]) != sign(w[end]) && abs(Giw) > 1e6 && wgm[i] < 0.001
             # This tries to filter out extremely large gain margins that can arise when the Nyquist contour crosses the negative real axis at -Inf.
             # This is filter is in addition to the filter_th check in _findCrossings
@@ -560,7 +567,7 @@ function sisomargin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMa
     pm = similar(wpm)
     for i = eachindex(wpm)
         # We have to access the actual phase value from the `phase` array to get unwrapped phase. This value is not fully accurate since it is computed at a grid point, so we compute the more accurate phase at the interpolated frequency. This accurate value is not unwrapped, so we add an integer multiple of 360 to get the closest unwrapped phase.
-        φ_nom = rad2deg(angle(freqresp(sys,wpm[i])[1]))
+        φ_nom = rad2deg(angle(freqresp(sys, wpm[i]; balance=false)[1]))
         φ_rounded = phase[clamp(round(Int, fi[i]), 1, length(phase))] # fi is interpolated, so we round to the closest integer
         φ_int = φ_nom - 360 * round( (φ_nom - φ_rounded) / 360 )
 
@@ -655,17 +662,19 @@ function _findCrossings(w, n, res; filter_th=Inf)
 end
 
 """
-    dₘ = delaymargin(G::LTISystem)
+    dₘ = delaymargin(G::LTISystem; balance=true)
 
 Return the delay margin, dₘ. For discrete-time systems, the delay margin is normalized by the sample time, i.e., the value represents the margin in number of sample times. 
 Only supports SISO systems.
 
 The delay margin is computed as the phase margin in radians divided by the cross-over frequency in rad/s. The delay margin is the maximum time delay that can be added to the system before it becomes unstable.
+
+`balance`: Call [`balance_statespace`](@ref) on the system before computing the frequency response, see [`freqresp`](@ref).
 """
-function delaymargin(G::LTISystem)
+function delaymargin(G::LTISystem; balance=true)
     # Phase margin in radians divided by cross-over frequency in rad/s.
     issiso(G) || error("delaymargin only supports SISO systems")
-    ωgₘ, gₘ, ωϕₘa, ϕₘa = margin(G,allMargins=true)
+    ωgₘ, gₘ, ωϕₘa, ϕₘa = margin(G; allMargins=true, balance)
     isempty(ϕₘa[1]) && return Inf
     ϕₘ, i = findmin(sign.(ωϕₘa[1]) .* ϕₘa[1]) # flip sign of negative frequency margins 
     ϕₘ   *= π/180

diff --git a/lib/ControlSystemsBase/src/delay_systems.jl b/lib/ControlSystemsBase/src/delay_systems.jl
@@ -1,8 +1,8 @@
-function freqresp!(R::Array{T,3}, sys::DelayLtiSystem, ω::AbstractVector{W}) where {T, W <: Real}
+function freqresp!(R::Array{T,3}, sys::DelayLtiSystem, ω::AbstractVector{W}; balance=true) where {T, W <: Real}
     ny = noutputs(sys)
     nu = ninputs(sys)
     @boundscheck size(R) == (ny,nu,length(ω))
-    P_fr = freqresp(sys.P.P, ω)
+    P_fr = freqresp(sys.P.P, ω; balance)
 
     cache = cis.(ω[1].*sys.Tau)