From 74740e8fd74055f0d075f1b60842939943b67afb Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Wed, 19 Apr 2023 08:06:11 +0200
Subject: [PATCH 01/24] handle breaking change in Polynomials

closes #828
---
 lib/ControlSystemsBase/Project.toml            | 2 +-
 lib/ControlSystemsBase/src/types/conversion.jl | 2 +-
 lib/ControlSystemsBase/test/test_conversion.jl | 5 +++++
 3 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml
index fc6b7ac9a..eef6e6e47 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/ControlSystemsBase.jl.git"
-version = "1.4.1"
+version = "1.4.2"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
diff --git a/lib/ControlSystemsBase/src/types/conversion.jl b/lib/ControlSystemsBase/src/types/conversion.jl
index c97071f70..6486dc5f9 100644
--- a/lib/ControlSystemsBase/src/types/conversion.jl
+++ b/lib/ControlSystemsBase/src/types/conversion.jl
@@ -137,7 +137,7 @@ function siso_tf_to_ss(T::Type, f::SisoRational)
     N = length(den) - 1 # The order of the rational function f
 
     # Get numerator coefficient of the same order as the denominator
-    bN = length(num) == N+1 ? num[1] : zero(num[1])
+    bN = length(num) == N+1 ? num[1] : zero(eltype(num))
 
     @views if N == 0 #|| num == zero(Polynomial{T})
         A = zeros(T, 0, 0)
diff --git a/lib/ControlSystemsBase/test/test_conversion.jl b/lib/ControlSystemsBase/test/test_conversion.jl
index a53aff512..25b80ef06 100644
--- a/lib/ControlSystemsBase/test/test_conversion.jl
+++ b/lib/ControlSystemsBase/test/test_conversion.jl
@@ -176,4 +176,9 @@ s3,s4 = promote(s1,s2)
 @test s3 == s4 == s2
 @test typeof(s3) == typeof(s4) == typeof(s2)
 
+# https://github.com/JuliaControl/ControlSystems.jl/issues/828
+mytf = [tf(1,[1,1]) tf(0); tf(0) tf(1,[1,1])]
+myss = ss(mytf)
+@test myss == ss([-1 0; 0 -1],[1 0; 0 1],[1 0; 0 1],[0 0; 0 0])
+
 end

From 604a9abbeac2a8ef9f314bb992d76b80e15fb519 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Wed, 26 Apr 2023 07:16:52 +0200
Subject: [PATCH 02/24] customizable default names

---
 lib/ControlSystemsBase/src/types/Lti.jl | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/lib/ControlSystemsBase/src/types/Lti.jl b/lib/ControlSystemsBase/src/types/Lti.jl
index fcd83adac..12fbed210 100644
--- a/lib/ControlSystemsBase/src/types/Lti.jl
+++ b/lib/ControlSystemsBase/src/types/Lti.jl
@@ -106,10 +106,10 @@ system_name(::LTISystem) = ""
 
 Get a vector of strings with the names of the inputs of `P`, or the `i`:th name if and index is given.
 """
-input_names(P::LTISystem) = [input_names(P, i) for i in 1:ninputs(P)]
-function input_names(P::LTISystem, i)
+input_names(P::LTISystem; kwargs...) = [input_names(P, i; kwargs...) for i in 1:ninputs(P)]
+function input_names(P::LTISystem, i; default = "u")
     i <= ninputs(P) || throw(BoundsError(P, i))
-    ninputs(P) == 1 && (return "u")
+    ninputs(P) == 1 && (return default)
     "u($i)"
 end
 
@@ -119,10 +119,10 @@ end
 
 Get a vector of strings with the names of the outputs of `P`, or the `i`:th name if and index is given.
 """
-output_names(P::LTISystem) = [output_names(P, i) for i in 1:noutputs(P)]
-function output_names(P::LTISystem, i)
+output_names(P::LTISystem; kwargs...) = [output_names(P, i; kwargs...) for i in 1:noutputs(P)]
+function output_names(P::LTISystem, i; default = "y")
     i <= noutputs(P) || throw(BoundsError(P, i))
-    noutputs(P) == 1 && (return "y")
+    noutputs(P) == 1 && (return default)
     "y($i)"
 end
 
@@ -132,9 +132,9 @@ end
 
 Get a vector of strings with the names of the states of `P`, or the `i`:th name if and index is given.
 """
-state_names(P::LTISystem) = [state_names(P, i) for i in 1:nstates(P)]
-function state_names(P::LTISystem, i)
+state_names(P::LTISystem; kwargs...) = [state_names(P, i; kwargs...) for i in 1:nstates(P)]
+function state_names(P::LTISystem, i; default = "x")
     i <= nstates(P) || throw(BoundsError(P, i))
-    nstates(P) == 1 && (return "x")
+    nstates(P) == 1 && (return default)
     "x($i)"
 end
\ No newline at end of file

From 1e1512e819ad064626abd82d732360b7bd27ea57 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Wed, 26 Apr 2023 09:18:17 +0200
Subject: [PATCH 03/24] better error message for improper delay conversion
 (#831)

* better error message for improper delay conversion

* create s var
---
 lib/ControlSystemsBase/src/types/DelayLtiSystem.jl | 4 ++++
 lib/ControlSystemsBase/test/test_conversion.jl     | 4 ++++
 2 files changed, 8 insertions(+)

diff --git a/lib/ControlSystemsBase/src/types/DelayLtiSystem.jl b/lib/ControlSystemsBase/src/types/DelayLtiSystem.jl
index e404b4b49..f9e94e151 100644
--- a/lib/ControlSystemsBase/src/types/DelayLtiSystem.jl
+++ b/lib/ControlSystemsBase/src/types/DelayLtiSystem.jl
@@ -59,6 +59,10 @@ function Base.convert(::Type{DelayLtiSystem{T1,S}}, d::T2) where {T1,T2 <: Abstr
 end
 
 function Base.convert(::Type{DelayLtiSystem{T,S}}, sys::TransferFunction{TE}) where {T,S,TE}
+    if issiso(sys) && length(numvec(sys.matrix[1,1])) > length(denvec(sys.matrix[1,1]))
+        error("The transfer function is not proper and can not be converted to a DelayLtiSystem type. If you tried to form the system `exp(sL) * B / A` where `B / A` is proper, add parenthesis to make it `exp(sL) * (B / A)`.")
+    end
+       
     DelayLtiSystem{T,S}(convert(StateSpace{TE, T}, sys))
 end
 # Catch convertsion between T
diff --git a/lib/ControlSystemsBase/test/test_conversion.jl b/lib/ControlSystemsBase/test/test_conversion.jl
index 25b80ef06..806c1f06d 100644
--- a/lib/ControlSystemsBase/test/test_conversion.jl
+++ b/lib/ControlSystemsBase/test/test_conversion.jl
@@ -62,6 +62,10 @@ H3 = zpk([(-1+im*sqrt(79))/8, (-1-im*sqrt(79))/8], [-3/2], 2)
 @test zpk(G3) ≈ H3 rtol=1e-15
 @test zpk(H3) == H3
 
+s = tf("s")
+@test_throws ErrorException delay(5)*((s+1))/((s+2)*(s+0.5))
+@test_throws ErrorException delay(5)*zpk((s+1))/zpk((s+2)*(s+0.5))
+
 # Test complex 1
 A = [1.0 + im 1; 0 -2-3im]
 B = [0;2]

From b2dff5de647e16b4e39702a72d3f34991f41be97 Mon Sep 17 00:00:00 2001
From: Albin Heimerson <albin.heimerson@control.lth.se>
Date: Fri, 28 Apr 2023 09:02:20 +0200
Subject: [PATCH 04/24] remove instance with h where it should be Ts

Looking through some old issues I saw #447 which seemed mostly resolved, though there was this one random occurrence of `h` which was not even referenced in the function body, so it felt like it should be removed since it could accidentally lead to bugs with silently accepting the `h` and returning something other than intended.
---
 lib/ControlSystemsBase/src/pid_design.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/lib/ControlSystemsBase/src/pid_design.jl b/lib/ControlSystemsBase/src/pid_design.jl
index 7696fb87b..49fbd181d 100644
--- a/lib/ControlSystemsBase/src/pid_design.jl
+++ b/lib/ControlSystemsBase/src/pid_design.jl
@@ -196,7 +196,7 @@ KN \\dfrac{s + b}{s + bN} = K \\dfrac{1 + s/b}{1 + s/(bN)}
 
 See also `leadlinkat` `laglink`
 """
-function leadlink(b, N, K=1; h=nothing, Ts=nothing)
+function leadlink(b, N, K=1; Ts=nothing)
     Ts !== nothing && (Ts ≥ 0 || throw(ArgumentError("Negative `Ts` is not supported.")))
     N > 1 || @warn "N should be ≥ 1 for the link to be phase advancing."
     numerator = [1/b, 1]
@@ -204,7 +204,6 @@ function leadlink(b, N, K=1; h=nothing, Ts=nothing)
     gain = K
     G = tf(gain*numerator,denominator)
     return  isnothing(Ts) ? G : c2d(G,Ts)
-
 end
 
 """

From dd4c9109ccf2880f688455c0fae6623d8d504a98 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 12 May 2023 11:58:57 +0200
Subject: [PATCH 05/24] apply numerical balancing before plotting

mitigates #835
---
 lib/ControlSystemsBase/src/plotting.jl | 50 ++++++++++++++++++--------
 1 file changed, 35 insertions(+), 15 deletions(-)

diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl
index 3886f602d..9bb23869e 100644
--- a/lib/ControlSystemsBase/src/plotting.jl
+++ b/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,11 +700,13 @@ _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.
 """
@@ -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)

From 5afca01df142bb2dbad28307ab5a09604b84af57 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 12 May 2023 12:21:01 +0200
Subject: [PATCH 06/24] add default fallback for balancing

---
 lib/ControlSystemsBase/src/types/conversion.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/lib/ControlSystemsBase/src/types/conversion.jl b/lib/ControlSystemsBase/src/types/conversion.jl
index 6486dc5f9..ee52021e9 100644
--- a/lib/ControlSystemsBase/src/types/conversion.jl
+++ b/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)`
 

From 3342dea7f9d9861e98494afba0b30a0f460f84e9 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 12 May 2023 12:36:17 +0200
Subject: [PATCH 07/24] add some documentation help on numerical balancing

---
 docs/src/man/numerical.md | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/docs/src/man/numerical.md b/docs/src/man/numerical.md
index 906631aa6..59bd596aa 100644
--- a/docs/src/man/numerical.md
+++ b/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
+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.
 

From 441a7a1dac7d7544087be052186246f6e340ef9d Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 12 May 2023 13:35:18 +0200
Subject: [PATCH 08/24] test fixes

---
 docs/src/man/numerical.md                           | 4 ++--
 lib/ControlSystemsBase/src/plotting.jl              | 2 +-
 lib/ControlSystemsBase/test/test_delayed_systems.jl | 2 +-
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/src/man/numerical.md b/docs/src/man/numerical.md
index 59bd596aa..6fd34d26b 100644
--- a/docs/src/man/numerical.md
+++ b/docs/src/man/numerical.md
@@ -30,8 +30,8 @@ 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
-2. Allocates some memory
+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)`.
 
diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl
index 9bb23869e..d8b6b8cac 100644
--- a/lib/ControlSystemsBase/src/plotting.jl
+++ b/lib/ControlSystemsBase/src/plotting.jl
@@ -710,7 +710,7 @@ Plot all the amplitude and phase margins of the system(s) `sys`.
 
 `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]
diff --git a/lib/ControlSystemsBase/test/test_delayed_systems.jl b/lib/ControlSystemsBase/test/test_delayed_systems.jl
index 4a69e825d..ceacd5872 100644
--- a/lib/ControlSystemsBase/test/test_delayed_systems.jl
+++ b/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)

From ba2713dd8142ee8e4c3ebfb48ee515f0dc7be6ce Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 12 May 2023 14:54:26 +0200
Subject: [PATCH 09/24] Update Project.toml

---
 lib/ControlSystemsBase/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml
index eef6e6e47..8eb570ad0 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/ControlSystemsBase.jl.git"
-version = "1.4.2"
+version = "1.4.3"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"

From d3c75c6ea8f4b4a4ef3ed0f421827b7bb065285f Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 12 May 2023 15:05:25 +0200
Subject: [PATCH 10/24] make `balance_statespace` work for more system types

closes #835
---
 lib/ControlSystemsBase/src/types/conversion.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/lib/ControlSystemsBase/src/types/conversion.jl b/lib/ControlSystemsBase/src/types/conversion.jl
index ee52021e9..bcdc6d9e2 100644
--- a/lib/ControlSystemsBase/src/types/conversion.jl
+++ b/lib/ControlSystemsBase/src/types/conversion.jl
@@ -187,9 +187,10 @@ end
 #     balance_statespace(A2, B2, C2, perm)
 # end
 
-function balance_statespace(sys::StateSpace, perm::Bool=false)
+function balance_statespace(sys::AbstractStateSpace, perm::Bool=false)
     A, B, C, T = balance_statespace(sys.A,sys.B,sys.C, perm)
-    return ss(A,B,C,sys.D,sys.timeevol), T
+
+    return basetype(sys)(A,B,C,sys.D, ntuple(i->getfield(sys, i+4), fieldcount(T)-4)...), T
 end
 
 # Method that might fail for some exotic types, such as TrackedArrays

From 5e121f57748f4e0e503ed3f354aa72d91fbdd144 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 12 May 2023 15:26:53 +0200
Subject: [PATCH 11/24] mention balancing in simulation section as well

---
 docs/src/man/numerical.md                      | 2 ++
 lib/ControlSystemsBase/src/types/conversion.jl | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/docs/src/man/numerical.md b/docs/src/man/numerical.md
index 6fd34d26b..5ca66768b 100644
--- a/docs/src/man/numerical.md
+++ b/docs/src/man/numerical.md
@@ -114,3 +114,5 @@ Linear systems with zero-order-hold inputs can be exactly simulated in discrete
 
 For discrete-time systems, the function [`lsim!`](@ref) accepts a pre-allocated workspace objects that can be used to avoid allocations for repeated simulations.
 
+### Numerical balancing
+If you are only interested in the simulated outputs, not the state trajectories, you may consider applying balancing to the statespace model using [`balance_statespace`](@ref) before simulating, see the section on [State-space balancing](@ref) above. If the state trajectories are of interest, balancing can still be performed before simulation, and the inverse transformation applied to the state trajectories after simulation.
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/src/types/conversion.jl b/lib/ControlSystemsBase/src/types/conversion.jl
index bcdc6d9e2..329dc0b33 100644
--- a/lib/ControlSystemsBase/src/types/conversion.jl
+++ b/lib/ControlSystemsBase/src/types/conversion.jl
@@ -190,7 +190,7 @@ end
 function balance_statespace(sys::AbstractStateSpace, perm::Bool=false)
     A, B, C, T = balance_statespace(sys.A,sys.B,sys.C, perm)
 
-    return basetype(sys)(A,B,C,sys.D, ntuple(i->getfield(sys, i+4), fieldcount(T)-4)...), T
+    return basetype(sys)(A,B,C,sys.D, ntuple(i->getfield(sys, i+4), fieldcount(sys)-4)...), T
 end
 
 # Method that might fail for some exotic types, such as TrackedArrays

From 3c7632f5a316abc320410e863978e11906103de4 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Sat, 13 May 2023 07:46:49 +0200
Subject: [PATCH 12/24] Update conversion.jl

---
 lib/ControlSystemsBase/src/types/conversion.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/ControlSystemsBase/src/types/conversion.jl b/lib/ControlSystemsBase/src/types/conversion.jl
index 329dc0b33..932fa2875 100644
--- a/lib/ControlSystemsBase/src/types/conversion.jl
+++ b/lib/ControlSystemsBase/src/types/conversion.jl
@@ -190,7 +190,7 @@ end
 function balance_statespace(sys::AbstractStateSpace, perm::Bool=false)
     A, B, C, T = balance_statespace(sys.A,sys.B,sys.C, perm)
 
-    return basetype(sys)(A,B,C,sys.D, ntuple(i->getfield(sys, i+4), fieldcount(sys)-4)...), T
+    return basetype(sys)(A,B,C,sys.D, ntuple(i->getfield(sys, i+4), fieldcount(typeof(sys))-4)...), T
 end
 
 # Method that might fail for some exotic types, such as TrackedArrays

From 3cf5dd730af1230658607b75f98b1e11ad48e4b7 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Sat, 13 May 2023 10:18:58 +0200
Subject: [PATCH 13/24] Update Project.toml

---
 lib/ControlSystemsBase/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml
index 8eb570ad0..3c6f894b7 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/ControlSystemsBase.jl.git"
-version = "1.4.3"
+version = "1.4.4"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"

From 642f35139871793d16e031eeb65741f875b6e916 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 15 May 2023 09:35:09 +0200
Subject: [PATCH 14/24] add examples to plot docs

---
 docs/src/lib/plotting.md | 68 ++++++++++++++++++++++++++++++++++------
 1 file changed, 58 insertions(+), 10 deletions(-)

diff --git a/docs/src/lib/plotting.md b/docs/src/lib/plotting.md
index 09deed54b..3ef5d386b 100644
--- a/docs/src/lib/plotting.md
+++ b/docs/src/lib/plotting.md
@@ -22,31 +22,59 @@ Private = false
 ### Bode plot
 
 ![bode](https://github.com/JuliaControl/ControlExamplePlots.jl/blob/master/src/figures/bode.png?raw=true)
-
+```julia
+tf1 = tf([1],[1,1])
+tf2 = tf([1/5,2],[1,1,1])
+sys = [tf1 tf2]
+ws = exp10.(range(-2,stop=2,length=200))
+bodeplot(sys, ws)
+```
 ### Sigma plot
 
 ![sigma](https://github.com/JuliaControl/ControlExamplePlots.jl/blob/master/src/figures/sigma.png?raw=true)
-
+```julia
+sys = ss([-1 2; 0 1], [1 0; 1 1], [1 0; 0 1], [0.1 0; 0 -0.2])
+sigmaplot(sys)
+```
 ### Margin
 
 ![margin](https://github.com/JuliaControl/ControlExamplePlots.jl/blob/master/src/figures/margin.png?raw=true)
-
+```julia
+tf1 = tf([1],[1,1])
+tf2 = tf([1/5,2],[1,1,1])
+ws = exp10.(range(-2,stop=2,length=200))
+marginplot([tf1, tf2], ws)
+```
 ### Gangoffour plot
 
 ![gangoffour](https://github.com/JuliaControl/ControlExamplePlots.jl/blob/master/src/figures/gangoffour.png?raw=true)
-
+```julia
+tf1 = tf([1.0],[1,1])
+gangoffourplot(tf1, [tf(1), tf(5)])
+```
 ### Nyquist plot
 
 ![nyquist](https://github.com/JuliaControl/ControlExamplePlots.jl/blob/master/src/figures/nyquist.png?raw=true)
-
+```julia
+sys = ss([-1 2; 0 1], [1 0; 1 1], [1 0; 0 1], [0.1 0; 0 -0.2])
+ws = exp10.(range(-2,stop=2,length=200))
+nyquistplot(sys, ws, Ms_circles=1.2, Mt_circles=1.2)
+```
 ### Nichols plot
 
 ![nichols](https://github.com/JuliaControl/ControlExamplePlots.jl/blob/master/src/figures/nichols.png?raw=true)
-
+```julia
+tf1 = tf([1],[1,1])
+ws = exp10.(range(-2,stop=2,length=200))
+nicholsplot(tf1,ws)
+```
 ### Pole-zero plot
 
 ![pzmap](https://github.com/JuliaControl/ControlExamplePlots.jl/blob/master/src/figures/pzmap.png?raw=true)
-
+```julia
+tf2 = tf([1/5,2],[1,1,1])
+pzmap(c2d(tf2, 0.1))
+```
 ### Rlocus plot
 
 ![rlocus](https://github.com/JuliaControl/ControlExamplePlots.jl/blob/master/src/figures/rlocus.png?raw=true)
@@ -54,12 +82,32 @@ Private = false
 ### Lsim response plot
 
 ![lsim](https://github.com/JuliaControl/ControlExamplePlots.jl/blob/master/src/figures/lsim.png?raw=true)
-
+```julia
+sys = ss([-1 2; 0 1], [1 0; 1 1], [1 0; 0 1], [0.1 0; 0 -0.2])
+sysd = c2d(sys, 0.01)
+L = lqr(sysd, [1 0; 0 1], [1 0; 0 1])
+ts = 0:0.01:5
+plot(lsim(sysd, (x,i)->-L*x, ts; x0=[1;2]), plotu=true)
+```
 ### Impulse response plot
 
 ![impulse](https://github.com/JuliaControl/ControlExamplePlots.jl/blob/master/src/figures/impulse.png?raw=true)
-
+```julia
+tf1 = tf([1],[1,1])
+tf2 = tf([1/5,2],[1,1,1])
+sys = [tf1 tf2]
+sysd = c2d(ss(sys), 0.01)
+plot(impulse(sysd, 5), l=:blue)
+```
 ### Step response plot
 
 ![step](https://github.com/JuliaControl/ControlExamplePlots.jl/blob/master/src/figures/step.png?raw=true)
-
+```julia
+tf1 = tf([1],[1,1])
+tf2 = tf([1/5,2],[1,1,1])
+sys = [tf1 tf2]
+sysd = c2d(ss(sys), 0.01)
+res = step(sysd, 5)
+plot(res, l=(:dash, 4))
+# plot!(stepinfo(step(sysd[1,1], 5))) # adds extra info to the plot
+```

From a00e241b31a080f5f6090c0f770a86f5a141a684 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 15 May 2023 09:50:55 +0200
Subject: [PATCH 15/24] add precompilation workload (#838)

* add precompilation workload

* add dep and compat
---
 lib/ControlSystemsBase/Project.toml           |  2 ++
 .../src/ControlSystemsBase.jl                 |  2 ++
 lib/ControlSystemsBase/src/precompilation.jl  | 31 +++++++++++++++++++
 3 files changed, 35 insertions(+)
 create mode 100644 lib/ControlSystemsBase/src/precompilation.jl

diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml
index 3c6f894b7..395112cca 100644
--- a/lib/ControlSystemsBase/Project.toml
+++ b/lib/ControlSystemsBase/Project.toml
@@ -15,6 +15,7 @@ MatrixEquations = "99c1a7ee-ab34-5fd5-8076-27c950a045f4"
 MatrixPencils = "48965c70-4690-11ea-1f13-43a2532b2fa8"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
@@ -33,6 +34,7 @@ MatrixEquations = "1, 2.1"
 MatrixPencils = "1.6"
 Polyester = "0.6, 0.7"
 Polynomials = "1.1.10, 2.0, 3.0"
+PrecompileTools = "1"
 RecipesBase = "1"
 StaticArrays = "1"
 julia = "1.6"
diff --git a/lib/ControlSystemsBase/src/ControlSystemsBase.jl b/lib/ControlSystemsBase/src/ControlSystemsBase.jl
index b1da3b541..d777244d8 100644
--- a/lib/ControlSystemsBase/src/ControlSystemsBase.jl
+++ b/lib/ControlSystemsBase/src/ControlSystemsBase.jl
@@ -233,4 +233,6 @@ function __init__()
     end
 end
 
+include("precompilation.jl")
+
 end
diff --git a/lib/ControlSystemsBase/src/precompilation.jl b/lib/ControlSystemsBase/src/precompilation.jl
new file mode 100644
index 000000000..5aa8ad6d0
--- /dev/null
+++ b/lib/ControlSystemsBase/src/precompilation.jl
@@ -0,0 +1,31 @@
+using PrecompileTools
+
+
+PrecompileTools.@setup_workload begin
+    P = ssrand(1,1,2, proper=true)
+    PrecompileTools.@compile_workload begin
+        Pd = c2d(P, 0.1)
+        for P in (P, Pd)
+            bode(P)
+            # nyquist(P)
+            step(P)
+            feedback(P)
+            feedback(P, P)
+            minreal(P)
+            balance_statespace(P)
+            P*P
+            P+P
+            2.0*P
+            # hinfnorm(P)
+            poles(P)
+            tzeros(P)
+        end
+        G = tf(1.0, [1.0, 1])
+        ss(G)
+        G = tf(1.0, [1.0, 1], 1)
+        ss(G)
+
+        # Pdel = P*delay(1.0)
+        # pade(Pdel, 2)
+    end
+end
\ No newline at end of file

From d39a2134c4c92242fd7ce9e03151d07cbd6ca7dd Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 15 May 2023 17:08:31 +0200
Subject: [PATCH 16/24] improve docs for magnitude scale

---
 lib/ControlSystemsBase/src/plotting.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl
index d8b6b8cac..c510360f7 100644
--- a/lib/ControlSystemsBase/src/plotting.jl
+++ b/lib/ControlSystemsBase/src/plotting.jl
@@ -10,7 +10,7 @@ _PlotScaleStr = ""
     setPlotScale(str)
 
 Set the default scale of magnitude in `bodeplot` and `sigmaplot`.
-`str` should be either `"dB"` or `"log10"`.
+`str` should be either `"dB"` or `"log10"`. The default scale if none is chosen is `"log10"`.
 """
 function setPlotScale(str::AbstractString)
     if str == "dB"
@@ -249,7 +249,7 @@ end
     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)`
+optionally provided. To change the Magnitude scale see [`setPlotScale`](@ref). The default magnitude scale is "log10" (absolute scale).
                                             
 - 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.

From f62c25c1859803d45deb2376de9adbcfc825c6db Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Wed, 17 May 2023 11:24:26 +0200
Subject: [PATCH 17/24] add test for MIMO extended_gangoffour

---
 .../src/sensitivity_functions.jl                  |  6 +++---
 lib/ControlSystemsBase/test/test_connections.jl   | 15 +++++++++++++++
 2 files changed, 18 insertions(+), 3 deletions(-)

diff --git a/lib/ControlSystemsBase/src/sensitivity_functions.jl b/lib/ControlSystemsBase/src/sensitivity_functions.jl
index fa393ced2..184b5d8bb 100644
--- a/lib/ControlSystemsBase/src/sensitivity_functions.jl
+++ b/lib/ControlSystemsBase/src/sensitivity_functions.jl
@@ -148,9 +148,9 @@ T  = G[2, 2]
 
 # For MIMO P
 S  = G[1:P.ny,     1:P.nu]
-PS = G[1:P.ny,     P.nu+1:end]
-CS = G[P.ny+1:end, 1:P.nu]
-T  = G[P.ny+1:end, P.nu+1:end]
+PS = G[1:P.ny,     P.ny+1:end]
+CS = G[P.ny+1:end, 1:P.ny]
+T  = G[P.ny+1:end, P.ny+1:end] # Input complimentary sensitivity function
 ```
 
 The gang of four can be plotted like so
diff --git a/lib/ControlSystemsBase/test/test_connections.jl b/lib/ControlSystemsBase/test/test_connections.jl
index 216040e6d..ce024bf10 100644
--- a/lib/ControlSystemsBase/test/test_connections.jl
+++ b/lib/ControlSystemsBase/test/test_connections.jl
@@ -393,6 +393,21 @@ G = extended_gangoffour(P, K, true)
 @test tf(G[2,2]) ≈ tf(-input_comp_sensitivity(P, K))
 
 
+## MIMO
+P = ssrand(3,2,1,proper=false)
+K = ssrand(2,3,1,proper=false)
+G = extended_gangoffour(P, K, false)
+S  = G[1:P.ny,     1:P.nu]
+PS = G[1:P.ny,     P.ny+1:end]
+CS = G[P.ny+1:end, 1:P.ny]
+T  = G[P.ny+1:end, P.ny+1:end]
+
+@test tf(S) ≈ tf(sensitivity(P, K))
+@test tf(CS) ≈ tf(G_CS(P, K))
+@test tf(PS) ≈ tf(G_PS(P, K))
+@test tf(T) ≈ tf(input_comp_sensitivity(P, K))
+
+
 # https://github.com/JuliaControl/ControlSystems.jl/issues/815
 function feedback_ctrl(G, K)
     ny,nu = size(G)

From 0bae63a1f164f2f62feba128207d340728384ead Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Wed, 17 May 2023 13:19:36 +0200
Subject: [PATCH 18/24] optimize performance of `sigma`

---
 lib/ControlSystemsBase/src/freqresp.jl       | 12 ++++++++++--
 lib/ControlSystemsBase/test/test_freqresp.jl | 18 ++++++++++++++++++
 2 files changed, 28 insertions(+), 2 deletions(-)

diff --git a/lib/ControlSystemsBase/src/freqresp.jl b/lib/ControlSystemsBase/src/freqresp.jl
index dbfe6ff15..a1820936a 100644
--- a/lib/ControlSystemsBase/src/freqresp.jl
+++ b/lib/ControlSystemsBase/src/freqresp.jl
@@ -358,10 +358,18 @@ end
 Compute the singular values `sv` of the frequency response of system `sys` at
 frequencies `w`. See also [`sigmaplot`](@ref)
 
-`sv` has size `(max(ny, nu), length(w))`""" 
+`sv` has size `(min(ny, nu), length(w))`""" 
 @autovec (1) function sigma(sys::LTISystem, w::AbstractVector)
     resp = freqresp(sys, w)
-    sv = dropdims(mapslices(svdvals, resp, dims=(1,2)),dims=2)
+    ny, nu = size(sys)
+    if ny == 1 || nu == 1 # Shortcut available
+        sv = Matrix{real(eltype(resp))}(undef, 1, length(w))
+        for i = eachindex(w)
+            @views sv[1, i] = norm(resp[:,:,i])
+        end
+    else
+        sv = dropdims(mapslices(svdvals, resp, dims=(1,2)),dims=2)::Matrix{real(eltype(resp))}
+    end
     return sv, w
 end
 @autovec (1) sigma(sys::LTISystem) = sigma(sys, _default_freq_vector(sys, Val{:sigma}()))
diff --git a/lib/ControlSystemsBase/test/test_freqresp.jl b/lib/ControlSystemsBase/test/test_freqresp.jl
index 43f0c4c48..89a385e75 100644
--- a/lib/ControlSystemsBase/test/test_freqresp.jl
+++ b/lib/ControlSystemsBase/test/test_freqresp.jl
@@ -128,6 +128,24 @@ end
 @test sigma(sys, ws)[1] ≈ sigs
 
 
+## Test sigma when either ny = 1 or nu = 1
+sys = ssrand(1,2,1)
+resp = freqresp(sys, ws)
+sigs = Array{Float64}(undef, 1, 50)
+for i in eachindex(ws)
+    sigs[:, i] =  svdvals(resp[:,:,i])
+end
+@test sigma(sys, ws)[1] ≈ sigs
+
+sys = ssrand(2,1,1)
+resp = freqresp(sys, ws)
+sigs = Array{Float64}(undef, 1, 50)
+for i in eachindex(ws)
+    sigs[:, i] =  svdvals(resp[:,:,i])
+end
+@test sigma(sys, ws)[1] ≈ sigs
+
+
 # test unwrap
 P = tf(1, [1, 1, 0]) * tf(1, [1, 1])
 w = exp10.(LinRange(-2, 4, 200))

From 863ec58b9d63b070908a553920a9bfcbf6542f7c Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Wed, 17 May 2023 13:57:56 +0200
Subject: [PATCH 19/24] try harder to return the correct type from matrix comps

---
 lib/ControlSystemsBase/src/types/StateSpace.jl | 6 +++++-
 lib/ControlSystemsBase/src/types/conversion.jl | 9 ++++++---
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/lib/ControlSystemsBase/src/types/StateSpace.jl b/lib/ControlSystemsBase/src/types/StateSpace.jl
index c0757f40c..0bc42662c 100644
--- a/lib/ControlSystemsBase/src/types/StateSpace.jl
+++ b/lib/ControlSystemsBase/src/types/StateSpace.jl
@@ -486,7 +486,11 @@ function minreal(sys::T, tol=nothing; fast=false, atol=0.0, kwargs...) where T <
         atol = tol
     end
     Ar, Br, Cr = MatrixPencils.lsminreal(A,B,C; atol, fast, kwargs...)
-    basetype(T)(Ar,Br,Cr,D, ntuple(i->getfield(sys, i+4), fieldcount(T)-4)...)
+    if hasfield(T, :sys)
+        basetype(T)(ss(Ar,Br,Cr,D), ntuple(i->getfield(sys, i+1), fieldcount(T)-1)...)
+    else
+        basetype(T)(Ar,Br,Cr,D, ntuple(i->getfield(sys, i+4), fieldcount(T)-4)...)
+    end
 end
 
 
diff --git a/lib/ControlSystemsBase/src/types/conversion.jl b/lib/ControlSystemsBase/src/types/conversion.jl
index 932fa2875..528923349 100644
--- a/lib/ControlSystemsBase/src/types/conversion.jl
+++ b/lib/ControlSystemsBase/src/types/conversion.jl
@@ -187,10 +187,13 @@ end
 #     balance_statespace(A2, B2, C2, perm)
 # end
 
-function balance_statespace(sys::AbstractStateSpace, perm::Bool=false)
+function balance_statespace(sys::S, perm::Bool=false) where S <: AbstractStateSpace
     A, B, C, T = balance_statespace(sys.A,sys.B,sys.C, perm)
-
-    return basetype(sys)(A,B,C,sys.D, ntuple(i->getfield(sys, i+4), fieldcount(typeof(sys))-4)...), T
+    if hasfield(S, :sys)
+        basetype(S)(ss(A,B,C,D), ntuple(i->getfield(sys, i+1), fieldcount(S)-1)...), T
+    else
+        basetype(S)(A,B,C,sys.D, ntuple(i->getfield(sys, i+4), fieldcount(S)-4)...), T
+    end
 end
 
 # Method that might fail for some exotic types, such as TrackedArrays

From 74855be8b13f48f180ddd6cedc36c8754f16fae6 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 19 May 2023 11:10:40 +0200
Subject: [PATCH 20/24] add AD to lqr opt

---
 docs/Project.toml                             |   3 +
 .../src/examples/automatic_differentiation.md | 146 ++++++++++++++----
 2 files changed, 117 insertions(+), 32 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 965226a0b..53e9fdfcf 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -7,9 +7,12 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 ForwardDiffChainRules = "c9556dd2-1aed-4cfe-8560-1557cf593001"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
+ImplicitDifferentiation = "57b37032-215b-411a-8a7c-41a003a55207"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationGCMAES = "6f0a0517-dbc2-4a7a-8a20-99ae7f27e911"
+OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
diff --git a/docs/src/examples/automatic_differentiation.md b/docs/src/examples/automatic_differentiation.md
index 70c1be86d..a1d08f11b 100644
--- a/docs/src/examples/automatic_differentiation.md
+++ b/docs/src/examples/automatic_differentiation.md
@@ -96,7 +96,7 @@ function plot_optimized(P, params, res)
     for (i,params) = enumerate((params, res))
         ls = (i == 1 ? :dash : :solid)
         lab = (i==1 ? "Initial" : "Optimized")
-        C, G, r1, r2 = systems(P, params)
+        C, G, r1, r2 = systems(params, P)
         mag = reshape(bode(G, Ω)[1], 4, :)'[:, [1, 2, 4]]
         plot!([r1, r2]; title="Time response", subplot=1,
             lab = lab .* [" \$r → e\$" " \$d → e\$"], legend=:bottomright, ls,
@@ -110,7 +110,7 @@ function plot_optimized(P, params, res)
     fig
 end
 
-function systems(P, params)
+function systems(params, P)
     kp,ki,kd,Tf = params # We optimize parameters in 
     C    = pid(kp, ki, kd; form=:parallel, Tf, state_space=true)
     G    = extended_gangoffour(P, C) # [S PS; CS T]
@@ -122,12 +122,12 @@ end
 
 σ(x) = 1/(1 + exp(-x)) # Sigmoid function used to obtain a smooth constraint on the peak of the sensitivity function
 
-@views function cost(P, params::AbstractVector{T}) where T
-    C, G, res1, res2 = systems(P, params)
+@views function cost(params::AbstractVector{T}, P) where T
+    C, G, res1, res2 = systems(params, P)
     R,_ = bode(G, Ω, unwrap=false)
     S = sum(σ.(100 .* (R[1, 1, :] .- Msc))) # max sensitivity
     CS = sum(Ω .* R[2, 1, :])               # frequency-weighted noise sensitivity
-    perf = mean(abs2, res1.y .* res1.t) + mean(abs2, res2.y .* res2.t)
+    perf = mean(abs2, res1.y .* res1.t') + mean(abs2, res2.y .* res2.t')
     return perf + Sweight*S + CSweight*CS # Blend all objectives together
 end
 
@@ -139,18 +139,62 @@ params  = [kp, ki, kd, 0.01] # Initial guess for parameters
 using Optimization
 using OptimizationGCMAES
 
-fopt = OptimizationFunction((x, _)->cost(P, x), Optimization.AutoForwardDiff())
-prob = OptimizationProblem(fopt, params, lb=zeros(length(params)), ub = 10ones(length(params)))
+fopt = OptimizationFunction(cost, Optimization.AutoForwardDiff())
+prob = OptimizationProblem(fopt, params, P, lb=zeros(length(params)), ub = 10ones(length(params)))
 solver = GCMAESOpt()
-res = solve(prob, solver; maxiters=1000); res.objective
+res = solve(prob, solver; maxiters=1000)
 plot_optimized(P, params, res.u)
 ```
 
 The optimized controller achieves more or less the same low peak in the sensitivity function, but does this while *both* making the step responses significantly faster *and* using much less controller gain for large frequencies (the orange sensitivity function), an altogether better tuning. The only potentially negative effect of this tuning is that the overshoot in response to a reference step increased slightly, indicated also by the slightly higher peak in the complimentary sensitivity function (green). However, the response to reference steps can (and most often should) be additionally shaped by reference pre-filtering (sometimes referred to as "feedforward" or "reference shaping"), by introducing an additional filter appearing in the feedforward path only, thus allowing elimination of the overshoot without affecting the closed-loop properties.
 
 ### LQG controller
-We could attempt a similar automatic tuning of an LQG controller. This time, we choose to optimize the weight matrices of the LQR problem and the state covariance matrix of the noise. Since this is a SISO system, we do not need to tune the control-input matrix or the measurement covariance matrix, since any non-unit weight assigned to those can be associated with the state matrices instead. Since these matrices are supposed to be positive semi-definite, we optimize Cholesky factors rather than the full matrices.
+We could attempt a similar automatic tuning of an LQG controller. This time, we choose to optimize the weight matrices of the LQR problem and the state covariance matrix of the noise. The synthesis of an LQR controller involves the solution of a Ricatti equation, which in turn involves performing a Schur decomposition. These steps hard hard to differentiate through in a conventional way, but we can make use of implicit differentiation using the implicit function theorem. To do so, we load the package `ImplicitDifferentiation`, and defines the conditions that hold at the solution of the Ricatti equaiton:
+```math
+A^TX + XA - XBR^{-1}B^T X + Q = 0
+```
+
+We then define differentiable versions of [`lqr`](@ref) and [`kalman`](@ref) that make use of the "implicit function".
+
+```@example autodiff
+using ImplicitDifferentiation
+
+function forward_are(Q0; A,B,R)
+    Q = reshape(Q0, size(A))
+    ControlSystemsBase.are(Continuous, A, B, Q, R), 0
+end
+
+function conditions_are(Q, X, noneed; A,B,R)
+    XB = X*B
+    T = promote_type(eltype(Q), eltype(X))
+    C = T.(Q)
+    mul!(C,A',X,1,1)
+    mul!(C,X,A,1,1)
+    mul!(C, XB, R\XB', -1, 1)
+end
+
+const implicit_are = ImplicitFunction(forward_are, conditions_are)
+
+function lqrdiff(P::AbstractStateSpace{Continuous}, Q, R)
+    (; A, B) = P
+    X0, _ = implicit_are(Q; A, B, R)
+    X = X0 isa AbstractMatrix ? X0 : reshape(X0, size(A))
+    R\(B'X)
+end
+
+function kalmandiff(P::AbstractStateSpace{Continuous}, Q, R)
+    A = P.A
+    B = P.C'
+    X0, _ = implicit_are(Q; A, B, R)
+    X = X0 isa AbstractMatrix ? X0 : reshape(X0, size(A))
+    (R\(B'X))'
+end
+```
+
+Since this is a SISO system, we do not need to tune the control-input matrix or the measurement covariance matrix, any non-unit weight assigned to those can be associated with the state matrices instead. Since these matrices are supposed to be positive semi-definite, we optimize Cholesky factors rather than the full matrices. We will also reformulate the problem slightly and use proper constraints to limit the sensitivity peak.
 ```@example autodiff
+using Ipopt, OptimizationMOI
+MOI = OptimizationMOI.MOI
 function triangular(x)
     m = length(x)
     n = round(Int, sqrt(2m-1))
@@ -164,36 +208,74 @@ function triangular(x)
 end
 invtriangular(T) = [T[i,j] for i = 1:size(T,1) for j = i:size(T,1)]
 
-function systems(P, params)
-    Qchol = triangular(params[1:10])
-    Rchol = triangular(params[11:20])
-    Q    = Qchol'Qchol
-    R    = Rchol'Rchol
-    L    = lqr(P, Q, I(1))
-    K    = kalman(P, R, I(1))
-    C    = observer_controller(P, L, K)
-    G    = extended_gangoffour(P, C) # [S PS; CS T]
-    Gd   = c2d(G, 0.1)   # Discretize the system
-    res1 = step(Gd[1,1], 0:0.1:15) # Simulate S
-    res2 = step(Gd[1,2], 0:0.1:15) # Simulate PS
-    C, G, res1, res2
+function systemslqr(params, P)
+    n2 = length(params) ÷ 2
+    Qchol = triangular(params[1:n2])
+    Rchol = triangular(params[n2+1:2n2])
+    Q = Qchol'Qchol
+    R = Rchol'Rchol
+    L = lqrdiff(P, Q, I(1))
+    K = kalmandiff(P, R, I(1))
+    C = observer_controller(P, L, K)
+    G = extended_gangoffour(P, C) # [S PS; CS T]
+    C, G
+end
+
+function systems(params, P)
+    C, G = systemslqr(params, P)
+    C, G, sim(G)...
 end
 
-Q0 = diagm([1,1,1,1]) # Initial guess LQR state penalty
-R0 = diagm([1,1,1,1]) # Initial guess Kalman state covariance
+function sim(G)
+    Gd = c2d(G, 0.1, :zoh)   # Discretize the system
+    res1 = step(Gd[1, 1], 0:0.1:15) # Simulate S
+    res2 = step(Gd[1, 2], 0:0.1:15) # Simulate PS
+    res1, res2
+end
+
+@views function cost(params::AbstractVector{T}, P) where T
+    CSweight = 0.001 # Noise amplification penalty
+    C, G = systemslqr(params, P)
+    res1, res2 = sim(G)
+    R, _ = bodev(G[2, 1], Ω; unwrap=false)
+    CS = sum(R .*= Ω)               # frequency-weighted noise sensitivity
+    perf = mean(abs2, res1.y .*= res1.t') + mean(abs2, res2.y .*= res2.t')
+    return perf + CSweight * CS # Blend all objectives together
+end
+
+@views function constraints(res, params::AbstractVector{T}, P) where T
+    C, G = systemslqr(params, P)
+    S, _ = bodev(G[1, 1], Ω; unwrap=false)
+    res .= maximum(S) # max sensitivity
+    nothing
+end
+
+Msc = 1.3   # Constraint on Ms
+
+Q0 = diagm([1.0, 1, 1, 1]) # Initial guess LQR state penalty
+R0 = diagm([1.0, 1, 1, 1]) # Initial guess Kalman state covariance
 params = [invtriangular(cholesky(Q0).U); invtriangular(cholesky(R0).U)]
 
-fopt = OptimizationFunction((x, _)->cost(P, x))
-prob = OptimizationProblem(fopt, params, lb=fill(-10, length(params)), ub = fill(10, length(params)))
-solver = GCMAESOpt()
-res = solve(prob, solver; maxiters=1000); res.objective
+fopt = OptimizationFunction(cost, Optimization.AutoForwardDiff(); cons=constraints)
+prob = OptimizationProblem(fopt, params, P;
+    lb    = fill(-10.0, length(params)),
+    ub    = fill(10.0, length(params)),
+    ucons = fill(Msc, 1),
+    lcons = fill(-Inf, 1),
+)
+solver = Ipopt.Optimizer()
+MOI.set(solver, MOI.RawOptimizerAttribute("print_level"), 0)
+MOI.set(solver, MOI.RawOptimizerAttribute("acceptable_tol"), 1e-1)
+MOI.set(solver, MOI.RawOptimizerAttribute("acceptable_constr_viol_tol"), 1e-2)
+MOI.set(solver, MOI.RawOptimizerAttribute("acceptable_iter"), 5)
+MOI.set(solver, MOI.RawOptimizerAttribute("hessian_approximation"), "limited-memory")
+
+res = solve(prob, solver)
 plot_optimized(P, params, res.u)
 ```
 
-This controller should perform better than the PID controller, which is known to be incapable of properly damping the resonance in a double-mass system. 
+This controller should perform better than the PID controller, which is known to be incapable of properly damping the resonance in a double-mass system. However, we did not include any integral action in the LQG controller, which has implication for the disturbance response, as indicated by the green step response in the simulation above.
 
-!!! note "No automatic differentiation"
-    This example did not use automatic differentiation like we did when optimizing the PID controller. The problematic functions are the ones that solve the Riccati equations, these make use of the Schur factorization which does not have differentiation rules defined.
 
 #### Robustness analysis
 To check the robustness of the designed LQG controller w.r.t. parametric uncertainty in the plant, we load the package [`MonteCarloMeasurements`](https://github.com/baggepinnen/MonteCarloMeasurements.jl) and recreate the plant model with 20% uncertainty in the spring coefficient.
@@ -202,7 +284,7 @@ using MonteCarloMeasurements
 Pu = DemoSystems.double_mass_model(k = Particles(32, Uniform(80, 120))) # Create a model with uncertainty in spring stiffness k ~ U(80, 120)
 unsafe_comparisons(true) # For the Bode plot to work
 
-C,_ = systems(P, res.u)             # Get the controller assuming P without uncertainty
+C,_ = systemslqr(res.u, P)             # Get the controller assuming P without uncertainty
 Gu = extended_gangoffour(Pu, C)     # Form the gang-of-four with uncertainty
 w = exp10.(LinRange(-1.5, 2, 500))
 bodeplot(Gu, w, plotphase=false, ri=false, N=32, ylims=(1e-1, 30), layout=1, sp=1, c=[1 2 4 3], lab=["S" "CS" "PS" "T"])

From cc5dcda88e003f652d2e13da6b13e9c4b8682aa5 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 19 May 2023 14:34:32 +0200
Subject: [PATCH 21/24] fix balance statespace

---
 lib/ControlSystemsBase/Project.toml            | 2 +-
 lib/ControlSystemsBase/src/types/conversion.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml
index 395112cca..0822e7bfc 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/ControlSystemsBase.jl.git"
-version = "1.4.4"
+version = "1.4.5"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
diff --git a/lib/ControlSystemsBase/src/types/conversion.jl b/lib/ControlSystemsBase/src/types/conversion.jl
index 528923349..5350bebb1 100644
--- a/lib/ControlSystemsBase/src/types/conversion.jl
+++ b/lib/ControlSystemsBase/src/types/conversion.jl
@@ -190,7 +190,7 @@ end
 function balance_statespace(sys::S, perm::Bool=false) where S <: AbstractStateSpace
     A, B, C, T = balance_statespace(sys.A,sys.B,sys.C, perm)
     if hasfield(S, :sys)
-        basetype(S)(ss(A,B,C,D), ntuple(i->getfield(sys, i+1), fieldcount(S)-1)...), T
+        basetype(S)(ss(A,B,C,sys.D), ntuple(i->getfield(sys, i+1), fieldcount(S)-1)...), T
     else
         basetype(S)(A,B,C,sys.D, ntuple(i->getfield(sys, i+4), fieldcount(S)-4)...), T
     end

From df73a7600c440d6d6e033bd2e866345bc8bb92e0 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 22 May 2023 15:50:04 +0200
Subject: [PATCH 22/24] Add interface to DSP (#843)

* add interface to dsp

* add interface to dsp

* add additional DSP methods

* add docs on filter design
---
 docs/src/lib/constructors.md                  |  1 +
 docs/src/man/creating_systems.md              | 17 +++++++-
 .../src/ControlSystemsBase.jl                 |  2 +
 lib/ControlSystemsBase/src/dsp.jl             | 39 +++++++++++++++++++
 lib/ControlSystemsBase/test/runtests.jl       |  3 +-
 lib/ControlSystemsBase/test/test_dsp.jl       | 30 ++++++++++++++
 src/dsp.jl                                    | 31 +++++++++++++++
 test/test_dsp.jl                              | 10 +++++
 8 files changed, 131 insertions(+), 2 deletions(-)
 create mode 100644 lib/ControlSystemsBase/src/dsp.jl
 create mode 100644 lib/ControlSystemsBase/test/test_dsp.jl
 create mode 100644 src/dsp.jl
 create mode 100644 test/test_dsp.jl

diff --git a/docs/src/lib/constructors.md b/docs/src/lib/constructors.md
index 444586128..57f244253 100644
--- a/docs/src/lib/constructors.md
+++ b/docs/src/lib/constructors.md
@@ -20,4 +20,5 @@ zpk
 delay
 pade
 ssdata
+ControlSystemsBase.seriesform
 ```
diff --git a/docs/src/man/creating_systems.md b/docs/src/man/creating_systems.md
index e2cfed794..836cb6191 100644
--- a/docs/src/man/creating_systems.md
+++ b/docs/src/man/creating_systems.md
@@ -433,4 +433,19 @@ uF  │       │ │  │       ├──────►         │ yC │  uP
                              └───────────────────────────────────┘
 ```
 
-See code example [complicated_feedback.jl](https://github.com/JuliaControl/RobustAndOptimalControl.jl/blob/master/examples/complicated_feedback.jl).
\ No newline at end of file
+See code example [complicated_feedback.jl](https://github.com/JuliaControl/RobustAndOptimalControl.jl/blob/master/examples/complicated_feedback.jl).
+
+## Filter design
+Filters can be designed using [DSP.jl](https://docs.juliadsp.org/stable/filters/). This results in filter objects with types from the DSP package, which can be converted to transfer functions using [`tf`](@ref) from ControlSystemsBase.
+
+```@example FilterDesign
+using DSP, ControlSystemsBase, Plots
+
+fs = 100
+df = digitalfilter(Bandpass(5, 10; fs), Butterworth(2))
+G = tf(df, 1/fs) # Sample time must be provided in the conversion to get the correct frequency scale in the Bode plot
+bodeplot(G, xscale=:identity, yscale=:identity, hz=true)
+```
+
+See also
+- [`ControlSystemsBase.seriesform`](@ref)
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/src/ControlSystemsBase.jl b/lib/ControlSystemsBase/src/ControlSystemsBase.jl
index d777244d8..6f3445549 100644
--- a/lib/ControlSystemsBase/src/ControlSystemsBase.jl
+++ b/lib/ControlSystemsBase/src/ControlSystemsBase.jl
@@ -129,6 +129,7 @@ import LinearAlgebra: BlasFloat
 export lyap # Make sure LinearAlgebra.lyap is available
 import Printf
 import Printf: @printf, @sprintf
+import DSP
 import DSP: conv
 using ForwardDiff
 import MatrixPencils
@@ -198,6 +199,7 @@ include("nonlinear_components.jl")
 include("types/staticsystems.jl")
 
 include("plotting.jl")
+include("dsp.jl")
 
 @deprecate pole poles
 @deprecate tzero tzeros
diff --git a/lib/ControlSystemsBase/src/dsp.jl b/lib/ControlSystemsBase/src/dsp.jl
new file mode 100644
index 000000000..c580926bc
--- /dev/null
+++ b/lib/ControlSystemsBase/src/dsp.jl
@@ -0,0 +1,39 @@
+
+tf(p::DSP.PolynomialRatio{:z}, h::Real = 1) = tf(DSP.coefb(p), DSP.coefa(p), h)
+tf(p::DSP.PolynomialRatio{:s}) = tf(DSP.coefb(p), DSP.coefa(p))
+tf(p::DSP.ZeroPoleGain{:z}, h::Real = 1) = tf(DSP.PolynomialRatio(p), h)
+tf(p::DSP.ZeroPoleGain{:s}) = tf(DSP.PolynomialRatio(p))
+
+function DSP.PolynomialRatio(G::TransferFunction{<:Discrete})
+    DSP.PolynomialRatio(
+        numvec(G)[],
+        denvec(G)[],
+    )
+end
+
+function TransferFunction(b::DSP.Biquad, h::Real = 1)
+    b0, b1, b2, a1, a2 = b.b0, b.b1, b.b2, b.a1, b.a2
+    tf([b0, b1, b2], [1, a1, a2], h)
+end
+
+
+"""
+    Gs, k = seriesform(G::TransferFunction{Discrete})
+
+Convert a transfer function `G` to a vector of second-order transfer functions and a scalar gain `k`, the product of which equals `G`.
+"""
+function seriesform(G::TransferFunction{<:Discrete})
+    Gs = DSP.SecondOrderSections(DSP.PolynomialRatio(G))
+    bqs = TransferFunction.(Gs.biquads, G.Ts)
+    bqs, Gs.g
+end
+
+
+function DSP.SecondOrderSections(G::TransferFunction{<:Discrete})
+    DSP.SecondOrderSections(DSP.PolynomialRatio(G))
+end
+
+function zpk(p::DSP.ZeroPoleGain, h::Real)
+    z,p,k = p.z, p.p, p.k
+    zpk(z, p, k, h)
+end
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/test/runtests.jl b/lib/ControlSystemsBase/test/runtests.jl
index 4f6dca27d..60bd2a2fe 100644
--- a/lib/ControlSystemsBase/test/runtests.jl
+++ b/lib/ControlSystemsBase/test/runtests.jl
@@ -39,7 +39,8 @@ my_tests = [
             "test_hammerstein_wiener",
             "test_demo_systems",
             "test_autovec",
-            "test_plots"
+            "test_plots",
+            "test_dsp",
             ]
 
 @testset "All Tests" begin
diff --git a/lib/ControlSystemsBase/test/test_dsp.jl b/lib/ControlSystemsBase/test/test_dsp.jl
new file mode 100644
index 000000000..532a460bb
--- /dev/null
+++ b/lib/ControlSystemsBase/test/test_dsp.jl
@@ -0,0 +1,30 @@
+@testset "DSP interoperability" begin
+    @info "Testing DSP interoperability"
+    import DSP 
+    G = DemoSystems.resonant()*DemoSystems.resonant(ω0=2) |> tf
+    Gd0 = c2d(DemoSystems.resonant(), 0.1) |> tf
+    Gd = c2d(G, 0.1)
+    Gs, k = ControlSystemsBase.seriesform(Gd)
+    @test k*prod(Gs) ≈ Gd atol=1e-3
+    Gds = DSP.SecondOrderSections(Gd)
+
+    u = randn(100)
+    uf = DSP.filt(Gds, u, zeros(2,2))
+    uls = lsim(Gd, u').y'
+    @test uf[1:end-1] ≈ uls[2:end]
+
+    z,p,k = randn(3), randn(3), randn()
+
+    f = DSP.ZeroPoleGain(z,p,k)
+    fcs = zpk(f, 1)
+
+    @test fcs.matrix[1].z ≈ z
+    @test fcs.matrix[1].p ≈ p
+    @test fcs.matrix[1].k ≈ k
+
+    u = randn(10)
+    uf = DSP.filt(f, u)
+    uls = lsim(fcs, u').y'
+    @test uf ≈ uls
+    
+end
diff --git a/src/dsp.jl b/src/dsp.jl
new file mode 100644
index 000000000..a301f413c
--- /dev/null
+++ b/src/dsp.jl
@@ -0,0 +1,31 @@
+ControlSystems.tf(p::DSP.PolynomialRatio, h::Real = 1) = tf(DSP.coefb(p), DSP.coefa(p), h)
+ControlSystems.tf(p::DSP.ZeroPoleGain, h::Real = 1) = tf(DSP.PolynomialRatio(p), h)
+
+DSP.PolynomialRatio(G::TransferFunction{<:Discrete}) = DSP.PolynomialRatio(numpoly(G)[], denpoly(G)[])
+
+function ControlSystems.TransferFunction(b::DSP.Biquad, h::Real = 1)
+    b0, b1, b2, a1, a2 = b.b0, b.b1, b.b2, b.a1, b.a2
+    tf([b0, b1, b2], [1, a1, a2], h)
+end
+
+
+"""
+    Gs, k = seriesform(G::TransferFunction)
+
+Convert a transfer function `G` to a vector of second-order transfer functions, the product of which equals `G`. `k` is a scalar gain.
+"""
+function seriesform(G::TransferFunction{<:Discrete})
+    Gs = DSP.SecondOrderSections(DSP.PolynomialRatio(G))
+    bqs = TransferFunction.(Gs.biquads, G.Ts)
+    bqs, Gs.g
+end
+
+
+function DSP.SecondOrderSections(G::TransferFunction{<:Discrete})
+    DSP.SecondOrderSections(DSP.PolynomialRatio(G))
+end
+
+function ControlSystems.zpk(p::DSP.ZeroPoleGain, h::Real)
+    z,p,k = p.z, p.p, p.z
+    zpk(z,p,k,h)
+end
\ No newline at end of file
diff --git a/test/test_dsp.jl b/test/test_dsp.jl
new file mode 100644
index 000000000..ae658cd56
--- /dev/null
+++ b/test/test_dsp.jl
@@ -0,0 +1,10 @@
+@testset "DSP interoperability" begin
+    @info "Testing DSP interoperability"
+    import DSP 
+    G = DemoSystems.resonant()*DemoSystems.resonant(ω0=2) |> tf
+    Gd = c2d(G, 0.1)
+    Gs, k = seriesform(Gd)
+    @test k*prod(Gs) ≈ Gd
+    Gds = DSP.SecondOrderSections(Gd)
+    
+end
\ No newline at end of file

From 32d76b45ebc54010a9a612facea1f29885a96341 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 22 May 2023 15:50:36 +0200
Subject: [PATCH 23/24] Update Project.toml

---
 lib/ControlSystemsBase/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml
index 0822e7bfc..5c65f35a4 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/ControlSystemsBase.jl.git"
-version = "1.4.5"
+version = "1.5.0"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"

From 544a3ff0b9e078fc7e507b6d40f88115e67f13eb Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 22 May 2023 16:38:25 +0200
Subject: [PATCH 24/24] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 5e84af8e4..a3f828e26 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,7 +2,7 @@ name = "ControlSystems"
 uuid = "a6e380b2-a6ca-5380-bf3e-84a91bcd477e"
 authors = ["Dept. Automatic Control, Lund University"]
 repo = "https://github.com/JuliaControl/ControlSystems.jl.git"
-version = "1.7.2"
+version = "1.7.3"
 
 [deps]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"