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"
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"])
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/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
+```
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/docs/src/man/numerical.md b/docs/src/man/numerical.md
index 906631aa6..5ca66768b 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 (but not the input-output mapping). If balancing is performed before simulation, the output will correspond to the output of the original system, but the state trajectory will not.
+2. Allocates some memory.
+
+Balancing is also automatically performed when a transfer function is converted to a statespace system using `ss(G)`, to convert without balancing, call `convert(StateSpace, G, balance=false)`.
+
+Intuitively (and simplified), balancing may be beneficial when the magnitude of the elements of the ``B`` matrix are vastly different from the magnitudes of the element of the ``C`` matrix, or when the ``A`` matrix contains very large coefficients. An example that exhibits all of these traits is the following
+```@example BALANCE
+using ControlSystemsBase, LinearAlgebra
+A = [-6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9 0.0; 0.0 -6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9; 2.0163803998106024e8 2.0163803998106024e8 -0.006223894167415392 -1.551620418759878e8 0.002358202548321148 0.002358202548321148; 0.0 0.0 5.063545034365582e-9 -0.4479539754649166 0.0 0.0; -2.824060629317756e8 2.0198389074625736e8 -0.006234569427701143 -1.5542817673286995e8 -0.7305736722226711 0.0023622473513548576; 2.0198389074625736e8 -2.824060629317756e8 -0.006234569427701143 -1.5542817673286995e8 0.0023622473513548576 -0.7305736722226711]
+B = [0.004019511633336128; 0.004019511633336128; 0.0; 0.0; 297809.51426114445; 297809.51426114445]
+C = [0.0 0.0 0.0 1.0 0.0 0.0]
+D = [0.0]
+linsys = ss(A,B,C,D)
+norm(linsys.A, Inf), norm(linsys.B, Inf), norm(linsys.C, Inf)
+```
+which after balancing becomes
+```@example BALANCE
+bsys, T = balance_statespace(linsys)
+norm(bsys.A, Inf), norm(bsys.B, Inf), norm(bsys.C, Inf)
+```
+If you plot the frequency-response of the two systems using [`bodeplot`](@ref), you'll see that they differ significantly (the balanced one is correct).
+
 ## Frequency-response calculation
 For small systems (small number of inputs, outputs and states), evaluating the frequency-response of a transfer function is reasonably accurate and very fast.
 
@@ -82,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/Project.toml b/lib/ControlSystemsBase/Project.toml
index fc6b7ac9a..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.1"
+version = "1.5.0"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
@@ -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..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
@@ -233,4 +235,6 @@ function __init__()
     end
 end
 
+include("precompilation.jl")
+
 end
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/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/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
 
 """
diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl
index 3886f602d..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"
@@ -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)`
+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.
+- If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- `balance`: Call [`balance_statespace`](@ref) on the system before plotting.
 
 `kwargs` is sent as argument to RecipesBase.plot.
 """
@@ -273,7 +274,7 @@ function _get_plotlabel(s, i, j)
     end
 end
 
-@recipe function bodeplot(p::Bodeplot; plotphase=true, ylimsphase=(), unwrap=true, hz=false)
+@recipe function bodeplot(p::Bodeplot; plotphase=true, ylimsphase=(), unwrap=true, hz=false, balance=true)
     systems, w = _processfreqplot(Val{:bode}(), p.args...)
     ws = (hz ? 1/(2π) : 1) .* w
     ny, nu = size(systems[1])
@@ -286,6 +287,9 @@ end
     link --> :x
 
     for (si,s) = enumerate(systems)
+        if balance
+            s = balance_statespace(s)[1]
+        end
         mag, phase = bode(s, w)[1:2]
         if _PlotScale == "dB" # Set by setPlotScale(str) globally
             mag = 20*log10.(mag)
@@ -378,13 +382,13 @@ optionally provided.
 - `Ms_circles`: draw circles corresponding to given levels of sensitivity (circles around -1 with  radii `1/Ms`). `Ms_circles` can be supplied as a number or a vector of numbers. A design staying outside such a circle has a phase margin of at least `2asin(1/(2Ms))` rad and a gain margin of at least `Ms/(Ms-1)`.
 - `Mt_circles`: draw circles corresponding to given levels of complementary sensitivity. `Mt_circles` can be supplied as a number or a vector of numbers.
 - `critical_point`: point on real axis to mark as critical for encirclements
-
-If `hz=true`, the hover information will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- If `hz=true`, the hover information will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- `balance`: Call [`balance_statespace`](@ref) on the system before plotting.
 
 `kwargs` is sent as argument to plot.
 """
 nyquistplot
-@recipe function nyquistplot(p::Nyquistplot; Ms_circles=Float64[], Mt_circles=Float64[], unit_circle=false, hz=false, critical_point=-1)
+@recipe function nyquistplot(p::Nyquistplot; Ms_circles=Float64[], Mt_circles=Float64[], unit_circle=false, hz=false, critical_point=-1, balance=true)
     systems, w = _processfreqplot(Val{:nyquist}(), p.args...)
     ny, nu = size(systems[1])
     nw = length(w)
@@ -394,6 +398,9 @@ nyquistplot
     θ = range(0, stop=2π, length=100)
     S, C = sin.(θ), cos.(θ)
     for (si,s) = enumerate(systems)
+        if balance
+            s = balance_statespace(s)[1]
+        end
         re_resp, im_resp = nyquist(s, w)[1:2]
         for j=1:nu
             for i=1:ny
@@ -645,17 +652,18 @@ end
 @userplot Sigmaplot
 """
     sigmaplot(sys, args...; hz=false)
-    sigmaplot(LTISystem[sys1, sys2...], args...; hz=false)
+    sigmaplot(LTISystem[sys1, sys2...], args...; hz=false, balance=true)
 
 Plot the singular values of the frequency response of the `LTISystem`(s). A
 frequency vector `w` can be optionally provided.
 
-If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- `balance`: Call [`balance_statespace`](@ref) on the system before plotting.
 
 `kwargs` is sent as argument to Plots.plot.
 """
 sigmaplot
-@recipe function sigmaplot(p::Sigmaplot; hz=false)
+@recipe function sigmaplot(p::Sigmaplot; hz=false, balance=true)
     systems, w = _processfreqplot(Val{:sigma}(), p.args...)
     ws = (hz ? 1/(2π) : 1) .* w
     ny, nu = size(systems[1])
@@ -664,6 +672,9 @@ sigmaplot
     xguide --> (hz ? "Frequency [Hz]" : "Frequency [rad/s]")
     yguide --> "Singular Values $_PlotScaleStr"
     for (si, s) in enumerate(systems)
+        if balance
+            s = balance_statespace(s)[1]
+        end
         sv = sigma(s, w)[1]'
         if _PlotScale == "dB"
             sv = 20*log10.(sv)
@@ -689,15 +700,17 @@ _to1series(y) = _to1series(1:size(y,3),y)
 
 @userplot Marginplot
 """
-    fig = marginplot(sys::LTISystem [,w::AbstractVector];  kwargs...)
-    marginplot(sys::Vector{LTISystem}, w::AbstractVector;  kwargs...)
+    fig = marginplot(sys::LTISystem [,w::AbstractVector];  balance=true, kwargs...)
+    marginplot(sys::Vector{LTISystem}, w::AbstractVector;  balance=true, kwargs...)
 
 Plot all the amplitude and phase margins of the system(s) `sys`.
-A frequency vector `w` can be optionally provided.
+
+- A frequency vector `w` can be optionally provided.
+- `balance`: Call [`balance_statespace`](@ref) on the system before plotting.
 
 `kwargs` is sent as argument to RecipesBase.plot.
 """
-@recipe function marginplot(p::Marginplot; plotphase=true)
+@recipe function marginplot(p::Marginplot; plotphase=true, hz=false, balance=true)
     systems, w = _processfreqplot(Val{:bode}(), p.args...)
     ny, nu = size(systems[1])
     s2i(i,j) = LinearIndices((nu,(plotphase ? 2 : 1)*ny))[j,i]
@@ -710,6 +723,9 @@ A frequency vector `w` can be optionally provided.
     layout --> (2ny, nu)
     label --> ""
     for (si, s) in enumerate(systems)
+        if balance
+            s = balance_statespace(s)[1]
+        end
         bmag, bphase = bode(s, w)
         for j=1:nu
             for i=1:ny
@@ -909,17 +925,18 @@ end
 @userplot Rgaplot
 """
     rgaplot(sys, args...; hz=false)
-    rgaplot(LTISystem[sys1, sys2...], args...; hz=false)
+    rgaplot(LTISystem[sys1, sys2...], args...; hz=false, balance=true)
 
 Plot the relative-gain array entries of the `LTISystem`(s). A
 frequency vector `w` can be optionally provided.
 
-If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.
+- `balance`: Call [`balance_statespace`](@ref) on the system before plotting.
 
 `kwargs` is sent as argument to Plots.plot.
 """
 rgaplot
-@recipe function rgaplot(p::Rgaplot; hz=false)
+@recipe function rgaplot(p::Rgaplot; hz=false, balance=true)
     systems, w = _processfreqplot(Val{:sigma}(), p.args...)
     ws = (hz ? 1/(2π) : 1) .* w
     ny, nu = size(systems[1])
@@ -928,6 +945,9 @@ rgaplot
     xguide --> (hz ? "Frequency [Hz]" : "Frequency [rad/s]")
     yguide --> "Element magnitudes"
     for (si, s) in enumerate(systems)
+        if balance
+            s = balance_statespace(s)[1]
+        end
         sv = abs.(relative_gain_array(s, w))
         for j in 1:size(sv, 1)
             for i in 1:size(sv, 2)
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
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/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/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
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 c97071f70..5350bebb1 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)
@@ -187,9 +187,13 @@ end
 #     balance_statespace(A2, B2, C2, perm)
 # end
 
-function balance_statespace(sys::StateSpace, 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 ss(A,B,C,sys.D,sys.timeevol), T
+    if hasfield(S, :sys)
+        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
 end
 
 # Method that might fail for some exotic types, such as TrackedArrays
@@ -212,6 +216,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)`
 
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_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)
diff --git a/lib/ControlSystemsBase/test/test_conversion.jl b/lib/ControlSystemsBase/test/test_conversion.jl
index a53aff512..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]
@@ -176,4 +180,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
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)
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/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))
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