diff --git a/docs/src/man/numerical.md b/docs/src/man/numerical.md
index 8bef8b97c..8e4284b82 100644
--- 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
index 303610731..db2dcbf49 100644
--- 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
index d28300391..42c86ed93 100644
--- 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
index 146469585..d37fa60a2 100644
--- a/lib/ControlSystemsBase/src/analysis.jl
+++ b/lib/ControlSystemsBase/src/analysis.jl
@@ -150,7 +150,7 @@ function det(sys::Matrix{S}) where {S<:SisoZpk}
 end
 
 """
-    dcgain(sys, ϵ=0)
+    dcgain(sys, ϵ=0; balance=true)
 
 Compute the dcgain of system `sys`.
 
@@ -158,8 +158,11 @@ 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
index 8477a4f8c..d1758388d 100644
--- 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)
 
diff --git a/lib/ControlSystemsBase/src/freqresp.jl b/lib/ControlSystemsBase/src/freqresp.jl
index 9ec2f9b3d..11d289d53 100644
--- a/lib/ControlSystemsBase/src/freqresp.jl
+++ b/lib/ControlSystemsBase/src/freqresp.jl
@@ -31,7 +31,8 @@ end
 
 BodemagWorkspace(sys::LTISystem, ω::AbstractVector) = BodemagWorkspace(sys, length(ω))
 
-function freqresp(sys::LTISystem, w::Real)
+function freqresp(sys::LTISystem, w::Real; balance=true)
+    sys = _balance(sys, balance)
     # Create imaginary freq vector s
     if iscontinuous(sys)
         s = im*w
@@ -41,10 +42,10 @@ function freqresp(sys::LTISystem, w::Real)
     evalfr(sys, s)
 end
 
-freqresp(G::Union{UniformScaling, AbstractMatrix, Number}, w::Real) = G
+freqresp(G::Union{UniformScaling, AbstractMatrix, Number}, w::Real; balance=true) = G
 
 """
-    sys_fr = freqresp(sys, w)
+    sys_fr = freqresp(sys, w; balance=true)
 
 Evaluate the frequency response of a linear system.
 
@@ -54,6 +55,7 @@ For discrete systems, computes `G(e^{jωT}) = C(e^{jωT}I - A)^{-1}B + D`
 # Arguments
 - `sys::LTISystem`: The system to analyze
 - `w::AbstractVector{<:Real}`: Frequency vector (rad/s)
+- `balance`: Call [`balance_statespace`](@ref) on the system before computing the frequency response. Balancing improves the numerical accuracy for poorly scaled systems and does not change the input-output map. Pass `balance=false` to avoid the (small) overhead.
 
 # Returns
 - `sys_fr`: Complex frequency response array of size `(ny, nu, length(w))`
@@ -68,20 +70,24 @@ w = exp10.(LinRange(-2, 2, 200))
 resp = freqresp(sys, w)
 ```
 """
-@autovec () function freqresp(sys::LTISystem, w_vec::AbstractVector{W}) where W <: Real
+@autovec () function freqresp(sys::LTISystem, w_vec::AbstractVector{W}; balance=true) where W <: Real
     te = timeevol(sys)
     ny,nu = noutputs(sys), ninputs(sys)
     T = promote_type(Complex{real(numeric_type(sys))}, Complex{W})
     R = Array{T, 3}(undef, ny, nu, length(w_vec))
-    freqresp!(R, sys, w_vec)
+    freqresp!(R, sys, w_vec; balance)
 end
 
 """
-    freqresp!(R::Array{T, 3}, sys::LTISystem, w_vec::AbstractVector{<:Real})
+    freqresp!(R::Array{T, 3}, sys::LTISystem, w_vec::AbstractVector{<:Real}; balance=true)
 
 In-place version of [`freqresp`](@ref) that takes a pre-allocated array `R` of size (ny, nu, nw)`
+
+Note that `balance=true` allocates a balanced copy of `sys`. Pass `balance=false` to keep the
+call free of allocations that scale with the state dimension.
 """
-function freqresp!(R::Array{T,3}, sys::LTISystem, w_vec::AbstractVector{<:Real}) where T
+function freqresp!(R::Array{T,3}, sys::LTISystem, w_vec::AbstractVector{<:Real}; balance=true) where T
+    sys = _balance(sys, balance)
     te = sys.timeevol
     ny,nu = noutputs(sys), ninputs(sys)
     @boundscheck size(R) == (ny,nu,length(w_vec))
@@ -91,7 +97,7 @@ function freqresp!(R::Array{T,3}, sys::LTISystem, w_vec::AbstractVector{<:Real})
     R
 end
 
-function freqresp!(R::Array{T,3}, sys::TransferFunction, w_vec::AbstractVector{<:Real}) where T
+function freqresp!(R::Array{T,3}, sys::TransferFunction, w_vec::AbstractVector{<:Real}; balance=true) where T
     te = sys.timeevol
     ny,nu = noutputs(sys), ninputs(sys)
     @boundscheck size(R) == (ny,nu,length(w_vec))
@@ -101,18 +107,18 @@ function freqresp!(R::Array{T,3}, sys::TransferFunction, w_vec::AbstractVector{<
     R
 end
 
-@autovec () function freqresp(G::AbstractMatrix, w_vec::AbstractVector{<:Real})
+@autovec () function freqresp(G::AbstractMatrix, w_vec::AbstractVector{<:Real}; balance=true)
     repeat(G, 1, 1, length(w_vec))
 end
 
-@autovec () function freqresp(G::Number, w_vec::AbstractVector{<:Real})
+@autovec () function freqresp(G::Number, w_vec::AbstractVector{<:Real}; balance=true)
     fill(G, 1, 1, length(w_vec))
 end
 
 _freq(w, ::Continuous) = complex(0, w)
 _freq(w, te::Discrete) = cis(w*te.Ts)
 
-@autovec () function freqresp!(R::Array{T,3}, sys::AbstractStateSpace, w_vec::AbstractVector{W}) where {T, W <: Real}
+@autovec () function freqresp!(R::Array{T,3}, sys::AbstractStateSpace, w_vec::AbstractVector{W}; balance=true) where {T, W <: Real}
     ny, nu = size(sys)
     @boundscheck size(R) == (ny,nu,length(w_vec))
     if sys.nx == 0 # Only D-matrix
@@ -121,13 +127,14 @@ _freq(w, te::Discrete) = cis(w*te.Ts)
         end
         return R
     end
+    sys = _balance(sys, balance)
     local F, Q
     try
         F = hessenberg(sys.A)
         Q = Matrix(F.Q)
     catch e
         # For matrix types that do not have a hessenberg implementation, we call the standard version of freqresp.
-        (e isa @static VERSION < v"1.12" ? Union{MethodError, ErrorException} : Union{MethodError, ErrorException, FieldError}) && return freqresp_nohess!(R, sys, w_vec)
+        (e isa @static VERSION < v"1.12" ? Union{MethodError, ErrorException} : Union{MethodError, ErrorException, FieldError}) && return freqresp_nohess!(R, sys, w_vec; balance=false) # sys is already balanced above
         # ErrorException appears if we try to access Q on a type which does not have Q as a field or property, notably HessenbergFactorization from GenericLinearAlgebra, on julia v1.12, this is instead a FieldError
         rethrow()
     end
@@ -205,21 +212,23 @@ function ldiv2!(u, cs, F::UpperHessenberg, B::AbstractVecOrMat; shift::Number=fa
     return X
 end
 
-function freqresp_nohess(sys::AbstractStateSpace, w_vec::AbstractVector{W}) where W <: Real
+function freqresp_nohess(sys::AbstractStateSpace, w_vec::AbstractVector{W}; balance=true) where W <: Real
     ny, nu = size(sys)
     T = promote_type(Complex{real(eltype(sys.A))}, Complex{W})
     R = Array{T, 3}(undef, ny, nu, length(w_vec))
-    freqresp_nohess!(R, sys, w_vec)
+    freqresp_nohess!(R, sys, w_vec; balance)
 end
 
 """
-    freqresp_nohess(sys::AbstractStateSpace, w_vec::AbstractVector{<:Real})
+    freqresp_nohess(sys::AbstractStateSpace, w_vec::AbstractVector{<:Real}; balance=true)
 
 Compute the frequency response of `sys` without forming a Hessenberg factorization.
 This function is called automatically if the Hessenberg factorization fails.
+
+- `balance`: Call [`balance_statespace`](@ref) on the system before computing the frequency response.
 """
 freqresp_nohess
-@autovec () function freqresp_nohess!(R::Array{T,3}, sys::AbstractStateSpace, w_vec::AbstractVector{W}) where {T, W <: Real}
+@autovec () function freqresp_nohess!(R::Array{T,3}, sys::AbstractStateSpace, w_vec::AbstractVector{W}; balance=true) where {T, W <: Real}
     ny, nu = size(sys)
     @boundscheck size(R) == (ny,nu,length(w_vec))
     nx = sys.nx
@@ -229,6 +238,7 @@ freqresp_nohess
         end
         return R
     end
+    sys = _balance(sys, balance)
     A,B,C0,D = ssdata(sys)
     C = complex.(C0) # We make C complex in order to not incur allocations in mul! below
     te = sys.timeevol
@@ -263,7 +273,9 @@ end
 Evaluate the transfer function of the LTI system sys
 at the complex number s=x (continuous-time) or z=x (discrete-time).
 
-For many values of `x`, use `freqresp` instead.
+For many values of `x`, use `freqresp` instead. Unlike [`freqresp`](@ref), this function does
+not balance the realization of `sys`, call [`balance_statespace`](@ref) manually if `sys` is
+poorly scaled.
 """
 function evalfr(sys::AbstractStateSpace, s::Number)
     T = _evalfr_return_type(sys, s)
@@ -312,7 +324,7 @@ function (sys::TransferFunction)(z_or_omegas::AbstractVector, map_to_unit_circle
 end
 
 """
-    mag, phase, w = bode(sys[, w]; unwrap=true)
+    mag, phase, w = bode(sys[, w]; unwrap=true, balance=true)
 
 Compute the magnitude and phase parts of the frequency response of system `sys`
 at frequencies `w`. The frequency response is evaluated as `G(jω)` for continuous
@@ -322,6 +334,7 @@ systems and `G(e^{jωT})` for discrete systems.
 - `sys::LTISystem`: The system to analyze
 - `w::AbstractVector`: Frequency vector (rad/s). If omitted, a default frequency range is used.
 - `unwrap::Bool`: If true (default), apply phase unwrapping to avoid discontinuities
+- `balance`: Call [`balance_statespace`](@ref) on the system before computing the frequency response, see [`freqresp`](@ref).
 
 # Returns
 - `mag`: Magnitude of frequency response, size `(ny, nu, length(w))`
@@ -339,14 +352,14 @@ sys = tf(1, [1, 1])
 mag, phase, w = bode(sys)
 ```
 """ 
-@autovec (1, 2) function bode(sys::LTISystem, w::AbstractVector; unwrap=true)
-    resp = freqresp(sys, w)
+@autovec (1, 2) function bode(sys::LTISystem, w::AbstractVector; unwrap=true, balance=true)
+    resp = freqresp(sys, w; balance)
     angles = angle.(resp)
     unwrap && unwrap!(angles,3)
     @. angles = rad2deg(angles)
     return abs.(resp), angles, w
 end
-@autovec (1, 2) bode(sys::LTISystem) = bode(sys, _default_freq_vector(sys, Val{:bode}()))
+@autovec (1, 2) bode(sys::LTISystem; unwrap=true, balance=true) = bode(sys, _default_freq_vector(sys, Val{:bode}()); unwrap, balance)
 
 # Performance difference between bode and bodemag for tf. Note how expensive the phase unwrapping is.
 # using ControlSystemsBase
@@ -361,7 +374,7 @@ end
 # # 2.991 ms (1 allocation: 64 bytes)
 
 """
-    mag = bodemag!(ws::BodemagWorkspace, sys::LTISystem, w::AbstractVector)
+    mag = bodemag!(ws::BodemagWorkspace, sys::LTISystem, w::AbstractVector; balance=true)
 
 Compute the Bode magnitude operating in-place on an instance of [`BodemagWorkspace`](@ref).
 
@@ -369,6 +382,7 @@ Compute the Bode magnitude operating in-place on an instance of [`BodemagWorkspa
 - `ws::BodemagWorkspace`: Pre-allocated workspace created with [`BodemagWorkspace`](@ref)
 - `sys::LTISystem`: The system to analyze
 - `w::AbstractVector`: Frequency vector (rad/s)
+- `balance`: Call [`balance_statespace`](@ref) on the system before computing the frequency response. This allocates a balanced copy of `sys`, pass `balance=false` to keep the call free of allocations that scale with the state dimension.
 
 # Returns
 - `mag`: Magnitude of frequency response, size `(ny, nu, length(w))`. 
@@ -392,20 +406,20 @@ ws = BodemagWorkspace(sys, w)
 mag = bodemag!(ws, sys, w)
 ```
 """
-function bodemag!(ws::BodemagWorkspace, sys::LTISystem, w::AbstractVector)
-    freqresp!(ws.R, sys, w)
+function bodemag!(ws::BodemagWorkspace, sys::LTISystem, w::AbstractVector; balance=true)
+    freqresp!(ws.R, sys, w; balance)
     @. ws.mag = abs(ws.R)
     ws.mag
 end
 
-function bodemag_nohess!(ws::BodemagWorkspace, sys::LTISystem, w::AbstractVector)
-    freqresp_nohess!(ws.R, sys, w)
+function bodemag_nohess!(ws::BodemagWorkspace, sys::LTISystem, w::AbstractVector; balance=true)
+    freqresp_nohess!(ws.R, sys, w; balance)
     @. ws.mag = abs(ws.R)
     ws.mag
 end
 
 """
-    re, img, w = nyquist(sys[, w])
+    re, img, w = nyquist(sys[, w]; balance=true)
 
 Compute the real and imaginary parts of the frequency response of system `sys`
 at frequencies `w`. The frequency response is evaluated as `G(jω)` for continuous
@@ -414,6 +428,7 @@ systems and `G(e^{jωT})` for discrete systems.
 # Arguments
 - `sys::LTISystem`: The system to analyze
 - `w::AbstractVector`: Frequency vector (rad/s). If omitted, a default frequency range is used.
+- `balance`: Call [`balance_statespace`](@ref) on the system before computing the frequency response, see [`freqresp`](@ref).
 
 # Returns
 - `re`: Real part of frequency response, size `(ny, nu, length(w))`
@@ -430,14 +445,14 @@ w = logspace(-2, 2, 100)
 re, img, w = nyquist(sys, w)
 ```
 """ 
-@autovec (1, 2) function nyquist(sys::LTISystem, w::AbstractVector)
-    resp = freqresp(sys, w)
+@autovec (1, 2) function nyquist(sys::LTISystem, w::AbstractVector; balance=true)
+    resp = freqresp(sys, w; balance)
     return real(resp), imag(resp), w
 end
-@autovec (1, 2) nyquist(sys::LTISystem) = nyquist(sys, _default_freq_vector(sys, Val{:nyquist}()))
+@autovec (1, 2) nyquist(sys::LTISystem; balance=true) = nyquist(sys, _default_freq_vector(sys, Val{:nyquist}()); balance)
 
 """
-    sv, w = sigma(sys[, w])
+    sv, w = sigma(sys[, w]; balance=true)
 
 Compute the singular values of the frequency response of system `sys` at
 frequencies `w`. The frequency response is evaluated as `G(jω)` for continuous
@@ -446,6 +461,7 @@ systems and `G(e^{jωT})` for discrete systems.
 # Arguments
 - `sys::LTISystem`: The system to analyze
 - `w::AbstractVector`: Frequency vector (rad/s). If omitted, a default frequency range is used.
+- `balance`: Call [`balance_statespace`](@ref) on the system before computing the frequency response, see [`freqresp`](@ref).
 
 # Returns
 - `sv`: Singular values of frequency response, size `(min(ny, nu), length(w))`
@@ -463,8 +479,8 @@ sys = ss([-1 0; 0 -2], [1 0; 0 1], [1 1; 0 1], 0)
 sv, w = sigma(sys)
 ```
 """ 
-@autovec (1,) function sigma(sys::LTISystem, w::AbstractVector)
-    resp = freqresp(sys, w)
+@autovec (1,) function sigma(sys::LTISystem, w::AbstractVector; balance=true)
+    resp = freqresp(sys, w; balance)
     ny, nu = size(sys)
     if ny == 1 || nu == 1 # Shortcut available
         sv = Matrix{real(eltype(resp))}(undef, 1, length(w))
@@ -476,7 +492,7 @@ sv, w = sigma(sys)
     end
     return sv, w
 end
-@autovec (1,) sigma(sys::LTISystem) = sigma(sys, _default_freq_vector(sys, Val{:sigma}()))
+@autovec (1,) sigma(sys::LTISystem; balance=true) = sigma(sys, _default_freq_vector(sys, Val{:sigma}()); balance)
 
 function _default_freq_vector(systems::Vector{<:LTISystem}, plot; adaptive=false)
     if adaptive
diff --git a/lib/ControlSystemsBase/src/hammerstein_weiner.jl b/lib/ControlSystemsBase/src/hammerstein_weiner.jl
index fd6e74572..c3b955d2d 100644
--- a/lib/ControlSystemsBase/src/hammerstein_weiner.jl
+++ b/lib/ControlSystemsBase/src/hammerstein_weiner.jl
@@ -1,5 +1,5 @@
-function freqresp!(R::Array{T,3}, sys::HammersteinWienerSystem, ω::AbstractVector{W}) where {T, W <: Real}
-    all(f isa Offset for f in sys.f) && return freqresp!(R, lft(sys.P.P, ss(I(length(sys.f)), timeevol(sys))), ω)
+function freqresp!(R::Array{T,3}, sys::HammersteinWienerSystem, ω::AbstractVector{W}; balance=true) where {T, W <: Real}
+    all(f isa Offset for f in sys.f) && return freqresp!(R, lft(sys.P.P, ss(I(length(sys.f)), timeevol(sys))), ω; balance)
     throw(ArgumentError("Frequency response is not defined for HammersteinWienerSystem with nonlinearities. Call linearize to obtain a linearized system"))
 end
 
diff --git a/lib/ControlSystemsBase/src/matrix_comps.jl b/lib/ControlSystemsBase/src/matrix_comps.jl
index 0f81b8655..03514ee65 100644
--- a/lib/ControlSystemsBase/src/matrix_comps.jl
+++ b/lib/ControlSystemsBase/src/matrix_comps.jl
@@ -273,7 +273,7 @@ covar(sys::TransferFunction, W) = covar(ss(sys), W)
 # Note: the H∞ norm computation is probably not as accurate as with SLICOT,
 # but this seems to be still reasonably ok as a first step
 """
-    norm(sys, p=2; tol=1e-6)
+    norm(sys, p=2; tol=1e-6, balance=true)
 
 `norm(sys)` or `norm(sys,2)` computes the H2 norm of the LTI system `sys`.
 
@@ -287,17 +287,19 @@ It represents the desired relative accuracy for the computed L∞ norm
 (this is not an absolute certificate however).
 
 `sys` is first converted to a `StateSpace` model if needed.
+
+`balance`: Call [`balance_statespace`](@ref) on the system before computing the norm. The norm depends on the system only through its input-output map, and balancing improves the numerical conditioning for poorly scaled systems.
 """
-function LinearAlgebra.norm(sys::AbstractStateSpace, p::Real=2; tol=1e-6)
+function LinearAlgebra.norm(sys::AbstractStateSpace, p::Real=2; tol=1e-6, balance=true)
     if p == 2
-        return sqrt(max(0,tr(covar(sys, I))))
+        return sqrt(max(0,tr(covar(_balance(sys, balance), I))))
     elseif p == Inf
-        return hinfnorm(sys; tol=tol)[1]
+        return hinfnorm(sys; tol, balance)[1]
     else
         error("`p` must be either `2` or `Inf`")
     end
 end
-LinearAlgebra.norm(sys::TransferFunction, p::Real=2; tol=1e-6) = norm(ss(sys), p, tol=tol)
+LinearAlgebra.norm(sys::TransferFunction, p::Real=2; tol=1e-6, balance=true) = norm(ss(sys), p; tol, balance)
 
 
 """
@@ -320,7 +322,7 @@ function schur_form(sys)
 end
 
 """
-    Ninf, ω_peak = hinfnorm(sys; tol=1e-6, resid_tol=nothing)
+    Ninf, ω_peak = hinfnorm(sys; tol=1e-6, resid_tol=nothing, balance=true)
 
 Compute the H∞ norm `Ninf` of the LTI system `sys`, together with a frequency
 `ω_peak` at which the gain Ninf is achieved.
@@ -342,6 +344,8 @@ boundary as genuine.
 
 `sys` is first converted to a state space model if needed.
 
+`balance`: Call [`balance_statespace`](@ref) on the system before computing the norm. The norm depends on the system only through its input-output map, and balancing improves the numerical conditioning for poorly scaled systems.
+
 The continuous-time L∞ norm computation implements the 'two-step algorithm' in:\\
 **N.A. Bruinsma and M. Steinbuch**, 'A fast algorithm to compute the H∞-norm of
 a transfer function matrix', Systems and Control Letters (1990), pp. 287-293.
@@ -352,12 +356,12 @@ state space systems in continuous and discrete time', American Control Conferenc
 
 See also [`linfnorm`](@ref).
 """
-hinfnorm(sys::AbstractStateSpace{<:Continuous}; tol=1e-6, resid_tol=nothing) = _infnorm_two_steps_ct(schur_form(sys)[1], :hinf, tol; resid_tol)
-hinfnorm(sys::AbstractStateSpace{<:Discrete}; tol=1e-6, resid_tol=nothing) = _infnorm_two_steps_dt(schur_form(sys)[1], :hinf, tol; resid_tol)
-hinfnorm(sys::TransferFunction; tol=1e-6, resid_tol=nothing) = hinfnorm(ss(sys); tol, resid_tol)
+hinfnorm(sys::AbstractStateSpace{<:Continuous}; tol=1e-6, resid_tol=nothing, balance=true) = _infnorm_two_steps_ct(schur_form(_balance(sys, balance))[1], :hinf, tol; resid_tol)
+hinfnorm(sys::AbstractStateSpace{<:Discrete}; tol=1e-6, resid_tol=nothing, balance=true) = _infnorm_two_steps_dt(schur_form(_balance(sys, balance))[1], :hinf, tol; resid_tol)
+hinfnorm(sys::TransferFunction; tol=1e-6, resid_tol=nothing, balance=true) = hinfnorm(ss(sys); tol, resid_tol, balance)
 
 """
-    Ninf, ω_peak = linfnorm(sys; tol=1e-6, resid_tol=nothing)
+    Ninf, ω_peak = linfnorm(sys; tol=1e-6, resid_tol=nothing, balance=true)
 
 Compute the L∞ norm `Ninf` of the LTI system `sys`, together with a frequency
 `ω_peak` at which the gain `Ninf` is achieved.
@@ -372,6 +376,8 @@ the computed L∞ norm (this is not an absolute certificate however).
 
 `sys` is first converted to a state space model if needed.
 
+`balance`: Call [`balance_statespace`](@ref) on the system before computing the norm. The norm depends on the system only through its input-output map, and balancing improves the numerical conditioning for poorly scaled systems.
+
 The continuous-time L∞ norm computation implements the 'two-step algorithm' in:\\
 **N.A. Bruinsma and M. Steinbuch**, 'A fast algorithm to compute the H∞-norm of
 a transfer function matrix', Systems and Control Letters (1990), pp. 287-293.
@@ -382,15 +388,15 @@ state space systems in continuous and discrete time', American Control Conferenc
 
 See also [`hinfnorm`](@ref).
 """
-function linfnorm(sys::AbstractStateSpace; tol=1e-6, resid_tol=nothing)
-    sys2, _ = schur_form(sys)
+function linfnorm(sys::AbstractStateSpace; tol=1e-6, resid_tol=nothing, balance=true)
+    sys2, _ = schur_form(_balance(sys, balance))
     if iscontinuous(sys2)
         return _infnorm_two_steps_ct(sys2, :linf, tol; resid_tol)
     else
         return _infnorm_two_steps_dt(sys2, :linf, tol; resid_tol)
     end
 end
-linfnorm(sys::TransferFunction; tol=1e-6, resid_tol=nothing) = linfnorm(ss(sys); tol, resid_tol)
+linfnorm(sys::TransferFunction; tol=1e-6, resid_tol=nothing, balance=true) = linfnorm(ss(sys); tol, resid_tol, balance)
 
 """
     _modal_residues(sys, suspects)
diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl
index 10730917d..f270279a2 100644
--- a/lib/ControlSystemsBase/src/plotting.jl
+++ b/lib/ControlSystemsBase/src/plotting.jl
@@ -293,15 +293,10 @@ _span(vec) = -(reverse(extrema(vec))...)
     grid   --> true
 
     for (si,s) = enumerate(systems)
-        if balance
-            sbal = balance_statespace(s)[1]
-        else
-            sbal = s
-        end
-        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 _PlotScale == "dB" # Set by setPlotScale(str) globally
             mag = 20*log10.(mag)
         elseif 0 ∈ mag
@@ -344,7 +339,7 @@ _span(vec) = -(reverse(extrema(vec))...)
                 end
                 plotphase || continue
 
-                if adjust_phase_start == true && isrational(sbal)
+                if adjust_phase_start == true && isrational(s)
                     if intexcess != 0
                         # Snap phase so that it starts at -90*intexcess
                         nineties = round(Int, phasedata[1] / 90)
@@ -437,10 +432,7 @@ 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]
+        re_resp, im_resp = nyquist(s, w; balance)[1:2]
         for j=1:nu
             for i=1:ny
                 redata = re_resp[i, j, :]
@@ -564,6 +556,8 @@ fontsize = 10
 
 `val` ∈ [0,1] determines the brightness of the gain lines
 
+`balance`: Call [`balance_statespace`](@ref) on the system before computing the frequency response, see [`freqresp`](@ref).
+
 Additional keyword arguments are sent to the function plotting the systems and can be
 used to specify colors, line styles etc. using regular RecipesBase.jl syntax
 
@@ -579,7 +573,8 @@ nicholsplot
     pInc     = 30,
     sat      = 0.4,
     val      = 0.85,
-    fontsize = 10)
+    fontsize = 10,
+    balance  = true)
 
     plots_id = Base.PkgId(UUID("91a5bcdd-55d7-5caf-9e0b-520d859cae80"), "Plots")
     haskey(Base.loaded_modules, plots_id) || error("Call using Plots before calling this function")
@@ -612,7 +607,7 @@ nicholsplot
     Ni_La(ϕ)        = @. 0.090*10^(ϕ/60)
     getColor(mdb)   = convert(Colors.RGB,Colors.HSV(360*((mdb-minimum(Gains))/(maximum(Gains)-minimum(Gains)))^1.5,sat,val))
 
-    megaangles      = vcat(map(s -> 180/π*angle(vec(freqresp(s, w))), systems)...)
+    megaangles      = vcat(map(s -> 180/π*angle(vec(freqresp(s, w; balance))), systems)...)
     filter!(x-> !isnan(x), megaangles)
     extremeangles = extrema(megaangles)
     extremeangles = floor(extremeangles[1]/180)*180, ceil(extremeangles[2]/180)*180
@@ -694,7 +689,7 @@ nicholsplot
     extremas = extrema(Gains)
     # colors = [:blue, :cyan, :green, :yellow, :orange, :red, :magenta]
     for (sysi,s) = enumerate(systems)
-        ℜresp, ℑresp        = nyquist(s, w)[1:2]
+        ℜresp, ℑresp        = nyquist(s, w; balance)[1:2]
         ℜdata               = dropdims(ℜresp, dims=(1,2))
         ℑdata               = dropdims(ℑresp, dims=(1,2))
         mag                 = 20*log10.(sqrt.(ℜdata.^2 + ℑdata.^2))
@@ -736,10 +731,7 @@ 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]'
+        sv = sigma(s, w; balance)[1]'
         if extrema && size(sv, 2) > 2
             sv = sv[:, [1, end]]
         end
@@ -794,14 +786,11 @@ marginplot
     layout --> (2ny, nu)
     label --> ""
     for (si, s) in enumerate(systems)
-        if balance
-            s = balance_statespace(s)[1]
-        end
-        bmag, bphase = bode(s, w)
+        bmag, bphase = bode(s, w; balance)
 
         for j=1:nu
             for i=1:ny
-                wgm, gm, wpm, pm, fullPhase, phasedata = sisomargin(s[i,j],w; full=true, allMargins=true, adjust_phase_start)
+                wgm, gm, wpm, pm, fullPhase, phasedata = sisomargin(s[i,j], w; full=true, allMargins=true, adjust_phase_start, balance)
                 if length(gm) > 5
                     @warn "Only showing smallest 5 out of $(length(gm)) gain margins"
                     idx = sortperm(gm)
@@ -1038,10 +1027,7 @@ 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))
+        sv = abs.(relative_gain_array(s, w; balance))
         for j in 1:size(sv, 1)
             for i in 1:size(sv, 2)
                 @series begin
diff --git a/lib/ControlSystemsBase/src/types/conversion.jl b/lib/ControlSystemsBase/src/types/conversion.jl
index b01990997..4c9c75ce1 100644
--- a/lib/ControlSystemsBase/src/types/conversion.jl
+++ b/lib/ControlSystemsBase/src/types/conversion.jl
@@ -237,6 +237,18 @@ end
 
 balance_statespace(sys, args...) = sys, I # For system types that do not have an implementation
 
+"""
+    _balance(sys, balance::Bool)
+
+Return `balance_statespace(sys)[1]` if `balance` is `true`, otherwise `sys` unchanged.
+
+This helper is used by functions that depend on `sys` only through its input-output map,
+for which balancing improves the numerical conditioning without changing the result. It is
+a no-op for system representations that have no state-space realization to balance, such as
+transfer functions.
+"""
+_balance(sys, balance::Bool) = balance ? balance_statespace(sys)[1] : sys
+
 """
 `T = balance_transform{R}(A::AbstractArray, B::AbstractArray, C::AbstractArray, perm::Bool=false)`
 
@@ -252,7 +264,7 @@ See also `balance_statespace`, `balance`
 function balance_transform(A::AbstractArray, B::AbstractArray, C::AbstractArray, perm::Bool=false)
     nx = size(A, 1)
     # Compute a scaling of the system matrix M
-    R = promote_type(eltype(A), eltype(B), eltype(C), Float32) # Make sure we get at least BlasFloat
+    R = promote_type(float(promote_type(eltype(A), eltype(B), eltype(C))), Float32) # Make sure we get at least BlasFloat. The inner `float` prevents an integer or rational system from being demoted to `Float32`.
     T = R[A B; C zeros(R, size(C*B))]
 
     size(T,1) < size(T,2) && (T = [T; zeros(R, size(T,2)-size(T,1),size(T,2))])
diff --git a/lib/ControlSystemsBase/src/types/staticsystems.jl b/lib/ControlSystemsBase/src/types/staticsystems.jl
index ef2445f42..1fabdb454 100644
--- a/lib/ControlSystemsBase/src/types/staticsystems.jl
+++ b/lib/ControlSystemsBase/src/types/staticsystems.jl
@@ -42,6 +42,24 @@ to_sized(sys::AbstractStateSpace) = HeteroStateSpace(sys, to_sized)
 
 StaticStateSpace(G::TransferFunction) = StaticStateSpace(ss(G))
 
+# The generic method for `AbstractStateSpace` returns a `StateSpace`, which would turn a
+# system with static matrices into one with dense matrices. The balancing transformation
+# itself is computed with dense matrices since it relies on LAPACK, the matrix container
+# type is restored afterwards. The element type is allowed to change, balancing a system
+# with integer matrices produces floating-point ones.
+function balance_statespace(sys::HeteroStateSpace, perm::Bool=false; kwargs...)
+    A, B, C, T = balance_statespace(_dense(sys.A), _dense(sys.B), _dense(sys.C), perm; kwargs...)
+    f = _matrix_container(sys.A)
+    HeteroStateSpace(f(A), f(B), f(C), sys.D, sys.timeevol), T
+end
+
+_dense(A::AbstractMatrix) = A
+_dense(A::Union{SArray, SizedArray}) = Matrix(A)
+
+_matrix_container(::AbstractMatrix) = identity
+_matrix_container(::SArray) = to_static
+_matrix_container(::SizedArray) = to_sized
+
 # function to_static(sys::DelayLtiSystem)
 #     innerP = to_static(sys.P.P)
 #     partP = PartitionedStateSpace(innerP, sys.P.nu1, sys.P.ny1)
@@ -201,7 +219,7 @@ end
 
 
 
-@autovec () function freqresp_nohess!(R::Array{T,3}, sys::StaticStateSpace, w_vec::AbstractVector{W}) where {T, W <: Real}
+@autovec () function freqresp_nohess!(R::Array{T,3}, sys::StaticStateSpace, w_vec::AbstractVector{W}; balance=true) where {T, W <: Real}
     ny, nu = size(sys)
     @boundscheck size(R) == (ny,nu,length(w_vec))
     nx = sys.nx
@@ -211,6 +229,7 @@ end
         end
         return R
     end
+    sys = _balance(sys, balance)
     A,B,C0,D = ssdata(sys)
     C = complex.(C0) # Still important when using ForwardDiff
     te = sys.timeevol
diff --git a/lib/ControlSystemsBase/test/test_freqresp.jl b/lib/ControlSystemsBase/test/test_freqresp.jl
index dda4a8402..f421a3e96 100644
--- a/lib/ControlSystemsBase/test/test_freqresp.jl
+++ b/lib/ControlSystemsBase/test/test_freqresp.jl
@@ -177,6 +177,63 @@ mag, mag, ws2 = bode(sys2)
 
 @test margin(tf(1, [1, -1], 0.01)).gm == [2;;]
 
+## Balancing ##
+# The functions that depend on the system only through its input-output map balance the
+# realization by default. The results must agree with the unbalanced computation for a
+# well-conditioned system, and the keyword must be accepted by every entry point.
+@testset "balance keyword" begin
+    Random.seed!(0)
+    Gbal = ssrand(2, 2, 4)
+    Gsiso = ssrand(1, 1, 4)
+    wbal = exp10.(LinRange(-2, 2, 50))
+
+    @test freqresp(Gbal, wbal) ≈ freqresp(Gbal, wbal; balance=false)
+    @test freqresp(Gbal, 1.0) ≈ freqresp(Gbal, 1.0; balance=false)
+    @test bode(Gbal, wbal)[1] ≈ bode(Gbal, wbal; balance=false)[1]
+    @test bode(Gbal)[1] ≈ bode(Gbal; balance=false)[1]
+    @test nyquist(Gbal, wbal)[1] ≈ nyquist(Gbal, wbal; balance=false)[1]
+    @test sigma(Gbal, wbal)[1] ≈ sigma(Gbal, wbal; balance=false)[1]
+    @test dcgain(Gbal) ≈ dcgain(Gbal; balance=false)
+    @test relative_gain_array(Gbal, wbal) ≈ relative_gain_array(Gbal, wbal; balance=false)
+    @test relative_gain_array(Gbal, 1.0) ≈ relative_gain_array(Gbal, 1.0; balance=false)
+    @test margin(Gsiso, wbal).gm ≈ margin(Gsiso, wbal; balance=false).gm
+    @test delaymargin(Gsiso) ≈ delaymargin(Gsiso; balance=false)
+    @test hinfnorm(Gbal)[1] ≈ hinfnorm(Gbal; balance=false)[1] rtol=1e-6
+    @test linfnorm(Gbal)[1] ≈ linfnorm(Gbal; balance=false)[1] rtol=1e-6
+    @test norm(Gbal) ≈ norm(Gbal; balance=false) rtol=1e-8
+    @test norm(Gbal, Inf) ≈ norm(Gbal, Inf; balance=false) rtol=1e-6
+
+    ws_bal = ControlSystemsBase.BodemagWorkspace(Gbal, wbal)
+    @test bodemag!(ws_bal, Gbal, wbal) ≈ bode(Gbal, wbal)[1]
+    @test bodemag!(ws_bal, Gbal, wbal; balance=false) ≈ bode(Gbal, wbal; balance=false)[1]
+
+    @test ControlSystemsBase.freqresp_nohess(Gbal, wbal) ≈ freqresp(Gbal, wbal)
+    @test ControlSystemsBase.freqresp_nohess(Gbal, wbal; balance=false) ≈ freqresp(Gbal, wbal)
+
+    # The autovec versions forward the keyword as well
+    @test bodev(Gsiso, wbal; balance=false)[1] ≈ bodev(Gsiso, wbal)[1]
+    @test nyquistv(Gsiso; balance=false)[1] ≈ nyquistv(Gsiso)[1]
+    @test sigmav(Gsiso; balance=false)[1] ≈ sigmav(Gsiso)[1]
+    @test freqrespv(Gsiso, wbal; balance=false) ≈ freqrespv(Gsiso, wbal)
+
+    # Representations without a state-space realization accept the keyword and ignore it
+    Ptf = tf(1, [1, 1])
+    @test freqresp(Ptf, wbal; balance=false) ≈ freqresp(Ptf, wbal)
+    @test freqresp(ss(2.0), wbal; balance=false) ≈ freqresp(ss(2.0), wbal)
+
+    # A poorly scaled realization is where balancing actually matters
+    Aill = [-1e6 1e5; 1e-5 -1.0]
+    Gill = ss(Aill, [1e-4; 1.0;;], [1.0 1e4], 0)
+    @test ControlSystemsBase.balance_statespace(Gill)[1].A != Gill.A
+    @test freqresp(Gill, wbal) ≈ freqresp(Gill, wbal; balance=false) rtol=1e-6
+
+    # Balancing keeps an integer realization at full precision, see issue with
+    # `balance_transform` demoting integer systems to Float32
+    Gint = ss([0 1; -1 -1], [0; 1], [1 0], 0)
+    @test eltype(ControlSystemsBase.balance_statespace(Gint)[2]) === Float64
+    @test hinfnorm(Gint, tol=1e-10)[1] ≈ 2/sqrt(3) rtol=1e-10
+end
+
 end
 
 
diff --git a/lib/ControlSystemsBase/test/test_staticsystems.jl b/lib/ControlSystemsBase/test/test_staticsystems.jl
index 6e4b7f0b9..94e3604b1 100644
--- a/lib/ControlSystemsBase/test/test_staticsystems.jl
+++ b/lib/ControlSystemsBase/test/test_staticsystems.jl
@@ -156,6 +156,20 @@ freqresp_nohess!(R, syss, w) # precompile
 D = StaticStateSpace(randn(2,2))
 @test freqresp_nohess(D, w) ≈ freqresp_nohess(ss(D), w)
 
+# Balancing a system with static matrices must not turn it into a dense `StateSpace`
+@test ControlSystemsBase.balance_statespace(syss)[1] isa StaticStateSpace
+@test freqresp(syss, w) ≈ freqresp(ss(syss), w)
+@test freqresp(syss, w; balance=false) ≈ freqresp(ss(syss), w)
+sys_ill = ss([-1e6 1e5; 1e-5 -1.0], [1e-4; 1.0;;], [1.0 1e4], 0)
+syss_ill = StaticStateSpace(sys_ill)
+@test ControlSystemsBase.balance_statespace(syss_ill)[1] isa StaticStateSpace
+@test ControlSystemsBase.balance_statespace(syss_ill)[1].A != syss_ill.A
+# Balancing an integer `HeteroStateSpace` promotes the element type rather than erroring
+sysh_int = HeteroStateSpace([-5 1; 0 -2], [2; 1;;], [3 0], [0;;], Continuous())
+@test ControlSystemsBase.balance_statespace(sysh_int)[1] isa HeteroStateSpace
+@test freqresp(sysh_int, w) ≈ freqresp(ss(sysh_int), w)
+@test ControlSystemsBase.balance_statespace(to_sized(sys_ill))[1].A isa SizedArray
+
 # Benchmarks with length(w) = 20000
 # @btime freqresp_nohess!(R, syss, w);
 # 29.660 ms (240001 allocations: 409.55 MiB)