diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index c1860bf58..991821197 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -26,7 +26,7 @@ jobs:
           - ControlSystems
           - ControlSystemsBase
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
@@ -48,7 +48,7 @@ jobs:
         # continue-on-error: ${{ matrix.version == 'nightly' }} # Allow nightly to fail and workflow still count as completed
       - uses: julia-actions/julia-processcoverage@v1
         with:
-          directories: lib/ControlSystemsBase/src
+          directories: src,lib/ControlSystemsBase/src
         if: ${{ matrix.version == '1' }}
       - uses: codecov/codecov-action@v3
         if: ${{ matrix.version == '1' }}
diff --git a/Project.toml b/Project.toml
index e742a88ac..dbb7732c2 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.8.2"
+version = "1.9.0"
 
 [deps]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
diff --git a/docs/make.jl b/docs/make.jl
index a44c24450..f6f4571d9 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -39,6 +39,7 @@ makedocs(modules=[ControlSystems, ControlSystemsBase],
             "Properties of delay systems" => "examples/delay_systems.md",
             "Automatic differentiation" => "examples/automatic_differentiation.md",
             "Tune a controller using experimental data" => "examples/tuning_from_data.md",
+            "Analysis of hybrid continuous/discrete systems" => "examples/zoh.md",
         ],
         "Functions" => Any[
             "Constructors" => "lib/constructors.md",
diff --git a/docs/src/examples/zoh.md b/docs/src/examples/zoh.md
new file mode 100644
index 000000000..5356e00b7
--- /dev/null
+++ b/docs/src/examples/zoh.md
@@ -0,0 +1,152 @@
+# Analysis of hybrid continuous-discrete systems
+In this example, we analyze the effect of ZoH sampling in continuous time and compare it to the equivalent discrete-time system
+
+We will consider a simple second-order system ``P``, which we define in continuous time:
+```@example zoh
+using ControlSystems, Plots
+s = tf('s')
+P = tf(0.1, [1, 0.1, 0.1])
+```
+
+## Continuous to discrete
+
+Next, we discretize this system using the standard [`c2d`](@ref) function, which uses ZoH sampling by default. We compare the frequency response of the discretized system with the frequency response of the original continuous-time system multiplied by the transfer function of the ZoH operator
+```math
+Z(s) = \dfrac{1 - e^{-T_s s}}{s}
+```
+
+```@example zoh
+Ts = 1 # Sample interval
+Z = (1 - delay(Ts))/s # The transfer function of the ZoH operator
+Pd = c2d(P, Ts) # Discrete-time system obtained by ZoH sampling
+Pz = P*Z # The continuous-time version of the discrete-time system
+
+wd = exp10.(-2:0.01:log10(2*0.5))
+bodeplot(P, wd, lab="\$P(s)\$")
+bodeplot!(Pz, wd, lab="\$P(s)Z(s)\$")
+bodeplot!(Pd, wd, lab="\$P_d(z)\$ (ZoH sampling)", l=:dash)
+vline!([0.5 0.5], l=(:black, :dash), lab="Nyquist freq.", legend=:bottomleft)
+```
+The frequency response of `Pz` ``= P(s) Z(s)`` matches that of ``P_d(z)`` exactly, but these two differ from the frequency response of the original ``P(s)`` due to the ZoH operator.
+
+The step response of `Pz` ``= P(s) Z(s)`` matches the discrete output of ``P_d(z)`` delayed by half the sample time
+```@example zoh
+resP = step(P, 12)
+resPz = step(Pz, 12)
+resPd = step(Pd, 12)
+plot([resP, resPz, resPd], lab=["P" "Pz" "Pd"])
+
+t_shift = resPd.t .+ Ts / 2
+plot!(t_shift, resPd.y[:], lab="Pd shifted", m=:o, l=:dash)
+```
+With a piecewise constant input, even if it's not a step, we get the same result
+```@example zoh
+Tf = 20
+ufun = (x,t)->[sin(2pi*floor(t/Ts)*Ts/5)]
+resP = lsim(P, ufun, Tf)
+resPz = lsim(Pz, ufun, 0:0.01:Tf)
+resPd = lsim(Pd, ufun, Tf)
+plot([resP, resPz, resPd], lab=["P" "Pz" "Pd"])
+
+t_shift = resPd.t .+ Ts / 2
+plot!(t_shift, resPd.y[:], lab="Pd shifted", m=:o)
+```
+
+With a _continuous_ input signal, the result is different,
+after the initial transient, the output of `Pz` matches that of `Pd` exactly
+(try plotting with the plotly() backend and zoom in at the end)
+```@example zoh
+Tf = 150
+ufun = (x,t)->[sin(2pi*t/5)]
+resP = lsim(P, ufun, Tf)
+resPz = lsim(Pz, ufun, 0:0.01:Tf)
+resPd = lsim(Pd, ufun, Tf)
+plot([resP, resPz, resPd], lab=["P" "Pz" "Pd"]); lens!([Tf-10, Tf], [-0.1, 0.1], inset=(1, bbox(0.4, 0.02, 0.4, 0.3)))
+```
+
+
+## Discrete to continuous
+
+Discrete-time systems can be converted to continuous-time systems formulated with delay terms in such a way that the frequency-response of the two systems match exactly, using the substiturion ``z^{-1} = e^{-sT_s}``. To this end, the function [`d2d_exact`](@ref) is used. This is useful when analyzing hybrid systems with both continuous-time and discrete-time components and accurate frequency-domain calculations are required over the entire frequency range up to the Nyquist frequency.
+
+Below, we compare the frequency response of a discrete-time system with delays to
+two continuous-time systems, one translated using the exact method and one using the standard `d2c` method with inverse ZoH sampling.
+We extend the frequency axis for this example to extend beyond the Nyquist frequency.
+```@example zoh
+wd = exp10.(LinRange(-2, 1, 200))
+Pdc = d2c_exact(ss(Pd))
+Pc = d2c(Pd)
+bodeplot(Pd, wd, lab="\$P_d(z)\$")
+bodeplot!(Pdc, wd, lab="\$P_d(s)\$ (exact translation)", l=:dash)
+bodeplot!(Pc, wd, lab="\$P_d(s)\$ (inverse ZoH sampling)")
+vline!([0.5 0.5], l=(:black, :dash), lab="Nyquist freq.", legend=:bottomleft)
+```
+We see that the translation of the discrete-time system to continuous time using the standard inverse ZoH sampling (``d2c(Pd)``) is not accurate for frequencies close to and above the Nyquist frequency. The translation using exact method (``d2c_exact(Pd)``) matches the frequency response of the discrete-time system exactly over the entire frequency axis.
+
+## Time-domain simulation
+
+When analyzing hybrid systems, systems with both discrete-time and continuous-time components, it is often useful to simulate the system in time-domain. To assemble the complete hybrid system, we must either
+1. Convert continuous-time components to discrete time, or
+2. Convert discrete-time components to continuous time.
+
+If all inputs to the continuous-time components are piecewise constant, then the first option is the most natural. The ZoH discretization is exact for piece-wise constant inputs. This conversion can be performed using the function [`c2d`](@ref) with the (default) ZoH sampling method. If some of the inputs to the continuous-time components are continuously varying, then the second option may be used for maximum accuracy. This is particularly useful if the continuous-time inputs contain frequencies higher than the Nyquist frequency, or when the continuous-time plant contains significant dynamics (resonances etc.) at frequencies higher than the Nyquist frequency of the discrete-time part. This conversion can be performed using the function [`d2c_exact`](@ref) which has two modes of operation. The default method produces a causal system which can be simulated in the time domain, while the second method produces an acausal system of lower complexity, which is amenable to frequency-domain computations only. Since we are going to simulate in the time domain here, we will use the causal method (default).
+
+In this example, we will model a continuous-time system ``P(s)`` with resonances appearing at large frequencies, while the discrete-time control system is operating a significantly lower frequencies.
+```@example zoh
+A = [0.0    0.0    0.0   1.0     0.0     0.0
+   0.0    0.0    0.0   0.0     1.0     0.0
+   0.0    0.0    0.0   0.0     0.0     1.0
+ -10.0    0.0   10.0  -0.001   0.0     0.001
+  -0.0  -10.0   10.0  -0.0    -0.001   0.001
+  10.0   10.0  -20.0   0.001   0.001  -0.002]
+B = [0, 0, 0, 0.1, 0, 0]
+C = [0.0  0.0  0.0  0.0  0.0  1.0]
+P = minreal(ss(A,B,C,0))
+w = exp10.(LinRange(-2, 2, 1000))
+bodeplot(P, w, lab="\$P(s)\$", plotphase=false)
+```
+The discrete-time controller ``C(z)`` is a basic PI controller
+
+```@example zoh
+Ts = 1
+C = pid(0.01,10; Ts, Tf = 1/100, state_space=true)
+```
+If we choose option 1. in this case, we get a discretized plant that has a very poor frequency-response fit to the original continuous-time plant, making frequency-domain analysis difficult
+```@example zoh
+Pd = c2d(P, Ts)
+Ld = Pd*C
+wd = exp10.(-3:0.001:log10(1000*0.05))
+figb = bodeplot(P, wd, lab="\$P(s)\$")
+bodeplot!(Pd, wd, lab="\$P(z)\$ ZoH")
+bodeplot!(Ld, wd, lab="\$P(z)C(z)\$", legend=:bottomleft)
+```
+
+If we instead make use of the second method above, exactly translating the discrete-time controller to continuous time, we can more easily determine that the closed-loop system will be stable by inspecting the Bode plot
+```@example zoh
+Cc = d2c_exact(C)
+Lc = P*Cc
+
+bodeplot(P, wd, lab="\$P(s)\$")
+bodeplot!(Lc, wd, lab="\$P(s)C(s)\$", legend=:bottomleft, c=3)
+```
+
+If we form the closed-loop system from an input disturbance to the output
+```math
+PS = \dfrac{P(s)}{1 + P(s)C(s)}
+```
+we can simulate the response to a continuous-time input disturbance that contains frequencies higher than the Nyquist frequency of the discrete-time system:
+```@example zoh
+PS = feedback(P, Cc) # Use the continuous-time plant and continuous translation of the controller
+fd = 10 # Frequency of input disturbance (Hz) (Nyquist frequency is 0.5Hz)
+disturbance(x, t) = sin(2pi*fd*t) + 0.02*(t >= 5) # Continuous input disturbance
+plot(lsim(PS, disturbance, 0:1e-3:10), title="Continuous disturbance response")
+```
+If we had tried doing this with the discretized plant, we would have gotten a very poor result
+```@example zoh
+PSd = feedback(Pd, C) # Use the discretized plant and discrete controller
+plot(lsim(PSd, disturbance, 10), ylims=(0.0, 0.005), title="Discrete disturbance response")
+```
+The output in this case is zero before the step in the disturbance at ``t = 5``, since the frequency of the input disturbance is at an integer multiple of the Nyquist frequency.
+
+
+## Summary
diff --git a/docs/src/man/numerical.md b/docs/src/man/numerical.md
index 8e3b5cf76..29798bfcb 100644
--- a/docs/src/man/numerical.md
+++ b/docs/src/man/numerical.md
@@ -116,4 +116,7 @@ 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
+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.
+
+## Static arrays in StateSpace systems
+The special statespace system type [HeteroStateSapce](@ref) can be used to store statespace models with static arrays rather than the default matrix type `Matrix`. See [State-Space Systems](@ref) for more details.
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml
index f258434a9..760a0d57d 100644
--- a/lib/ControlSystemsBase/Project.toml
+++ b/lib/ControlSystemsBase/Project.toml
@@ -2,7 +2,7 @@ name = "ControlSystemsBase"
 uuid = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 authors = ["Dept. Automatic Control, Lund University"]
 repo = "https://github.com/JuliaControl/ControlSystems.jl.git"
-version = "1.9.5"
+version = "1.10.0"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
@@ -41,11 +41,16 @@ MacroTools = "0.5"
 MatrixEquations = "1, 2.1"
 MatrixPencils = "1.6"
 Polyester = "0.6, 0.7"
-Polynomials = "1.1.10, 2.0, 3.0"
+Polynomials = "3.0, 4.0"
 PrecompileTools = "1"
 RecipesBase = "1"
 StaticArraysCore = "1"
 julia = "1.6"
+LinearAlgebra = "<0.0.1, 1"
+Printf = "<0.0.1, 1"
+Random = "<0.0.1, 1"
+UUIDs = "<0.0.1, 1"
+SparseArrays = "<0.0.1, 1"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
diff --git a/lib/ControlSystemsBase/src/ControlSystemsBase.jl b/lib/ControlSystemsBase/src/ControlSystemsBase.jl
index 483cb1aee..a93e71a37 100644
--- a/lib/ControlSystemsBase/src/ControlSystemsBase.jl
+++ b/lib/ControlSystemsBase/src/ControlSystemsBase.jl
@@ -84,6 +84,7 @@ export  LTISystem,
         c2d,
         c2d_x0map,
         d2c,
+        d2c_exact,
         # Time Response
         step,
         impulse,
diff --git a/lib/ControlSystemsBase/src/discrete.jl b/lib/ControlSystemsBase/src/discrete.jl
index a2ebcac90..3d101c533 100644
--- a/lib/ControlSystemsBase/src/discrete.jl
+++ b/lib/ControlSystemsBase/src/discrete.jl
@@ -91,6 +91,8 @@ end
 Convert discrete-time system to a continuous time system, assuming that the discrete-time system was discretized using `method`. Available methods are `:zoh, :fwdeuler´.
 
 - `w_prewarp`: Frequency for pre-warping when using the Tustin method, has no effect for other methods.
+
+See also [`d2c_exact`](@ref).
 """
 function d2c(sys::AbstractStateSpace{<:Discrete}, method::Symbol=:zoh; w_prewarp=0)
     A, B, C, D = ssdata(sys)
@@ -125,6 +127,38 @@ end
 
 d2c(sys::TransferFunction{<:Discrete}, args...) = tf(d2c(ss(sys), args...))
 
+
+"""
+    d2c_exact(sys::AbstractStateSpace{<:Discrete}, method = :causal)
+
+Translate a discrete-time system to a continuous-time system by one of the substitutions
+- ``z^{-1} = e^{-sT_s}`` if `method = :causal` (default)
+- ``z = e^{sT_s}``  if `method = :acausal`
+The translation is exact in the frequency domain, i.e.,
+the frequency response of the resulting continuous-time system is identical to
+the frequency response of the discrete-time system.
+
+This method of translation is useful when analyzing hybrid continuous/discrete systems in the frequency domain and high accuracy is required.
+
+The resulting system will be be a static system in feedback with pure delays. When `method = :causal`, the delays will be positive, resulting in a causal system that can be simulated in the time domain. When `method = :acausal`, the delays will be negative, resulting in an acausal system that **can not** be simulated in the time domain. The acausal translation results in a smaller system with half as many delay elements in the feedback path.
+"""
+function d2c_exact(sys::AbstractStateSpace{<:Discrete}, method=:causal)
+    T = sys.Ts
+    A,B,C,D = ssdata(sys)
+    if method === :acausal
+        z = delay(-T)
+        LR = append([z for _ in 1:sys.nx]...) - ss(A + I)
+        C*feedback(I(sys.nx), LR)*B + D
+    elseif method === :causal
+        z1 = delay(T)
+        ZI1 = append([z1 for _ in 1:sys.nx]...)
+        LR = ZI1 * ss(-A)
+        C*ZI1*feedback(I(sys.nx), LR)*B + D
+    else
+        error("Unknown method: $method. Choose method = :acausal or :causal")
+    end
+end
+
 # c2d and d2c for covariance and cost matrices =================================
 """
     Qd     = c2d(sys::StateSpace{Continuous}, Qc::Matrix, Ts;             opt=:o)
diff --git a/lib/ControlSystemsBase/src/pid_design.jl b/lib/ControlSystemsBase/src/pid_design.jl
index 65e04a910..1641a80e6 100644
--- a/lib/ControlSystemsBase/src/pid_design.jl
+++ b/lib/ControlSystemsBase/src/pid_design.jl
@@ -10,7 +10,7 @@ The `form` can be chosen as one of the following
 * `:series` - `Kc*(1 + 1/(τi*s))*(τd*s + 1)`
 * `:parallel` - `Kp + Ki/s + Kd*s`
 
-If `state_space` is set to `true`, either `kd` has to be zero 
+If `state_space` is set to `true`, either `Kd` has to be zero
 or a positive `Tf` has to be provided for creating a filter on 
 the input to allow for a state space realization. 
 The filter used is `1 / (1 + s*Tf + (s*Tf)^2/2)`, where `Tf` can typically 
@@ -47,44 +47,42 @@ end
 
 @deprecate pid(; kp=0, ki=0, kd=0, series = false) pid(kp, ki, kd; form=series ? :series : :parallel)
 
-function pid_tf(param_p, param_i, param_d=zero(typeof(param_p)); form=:standard,  Tf=nothing)
-    Kp, Ti, Td = convert_pidparams_to_standard(param_p, param_i, param_d, form)
-    ia = Ti != Inf && Ti != 0 # integral action, 0 would result in division by zero, but typically indicates that the user wants no integral action
+function pid_tf(param_p, param_i, param_d=zero(typeof(param_p)); form=:standard, Tf=nothing)
+    Kp, Ki, Kd = convert_pidparams_to_parallel(param_p, param_i, param_d, form)
     if isnothing(Tf)
-        if ia
-            return tf([Kp * Td, Kp, Kp / Ti], [1, 0])
+        if Ki != 0
+            return tf([Kd, Kp, Ki], [1, 0])
         else
-            return tf([Kp * Td, Kp], [1])
+            return tf([Kd, Kp], [1])
         end
     else
-        if ia
-            return tf([Kp * Td, Kp, Kp / Ti], [Tf^2/2, Tf, 1, 0])
+        if Ki != 0
+            return tf([Kd, Kp, Ki], [Tf^2/2, Tf, 1, 0])
         else
-            return tf([Kp * Td, Kp], [Tf^2/2, Tf, 1])
+            return tf([Kd, Kp], [Tf^2/2, Tf, 1])
         end
     end
 end
 
-function pid_ss(param_p, param_i, param_d=zero(typeof(param_p)); form=:standard,  Tf=nothing)
-    Kp, Ti, Td = convert_pidparams_to_standard(param_p, param_i, param_d, form)
+function pid_ss(param_p, param_i, param_d=zero(typeof(param_p)); form=:standard, Tf=nothing)
+    Kp, Ki, Kd = convert_pidparams_to_parallel(param_p, param_i, param_d, form)
     TE = Continuous()
-    ia = Ti != Inf && Ti != 0 # integral action, 0 would result in division by zero, but typically indicates that the user wants no integral action
     if !isnothing(Tf)
-        if ia
+        if Ki != 0
             A = [0 1 0; 0 0 1; 0 -2/Tf^2 -2/Tf]
             B = [0; 0; 1]
-            C = 2 * Kp / Tf^2 * [1/Ti 1 Td]
+            C = 2 / Tf^2 * [Ki Kp Kd]
         else
             A = [0 1; -2/Tf^2 -2/Tf]
             B = [0; 1]
-            C = 2 * Kp / Tf^2 * [1 Td]
+            C = 2 / Tf^2 * [Kp Kd]
         end
         D = 0
-    elseif Td == 0
-        if ia
+    elseif Kd == 0
+        if Ki != 0
             A = 0
             B = 1
-            C = Kp / Ti # Ti == 0 would result in division by zero, but typically indicates that the user wants no integral action
+            C = Ki # Ti == 0 would result in division by zero, but typically indicates that the user wants no integral action
             D = Kp
         else
             return ss([Kp])
@@ -98,7 +96,7 @@ end
 """
     pidplots(P, args...; params_p, params_i, params_d=0, form=:standard, ω=0, grid=false, kwargs...)
 
-Plots interesting figures related to closing the loop around process `P` with a PID controller supplied in `params`
+Display the relevant plots related to closing the loop around process `P` with a PID controller supplied in `params`
 on one of the following forms:
 * `:standard` - `Kp*(1 + 1/(Ti*s) + Td*s)` 
 * `:series` - `Kc*(1 + 1/(τi*s))*(τd*s + 1)`
@@ -267,7 +265,8 @@ function stabregionPID(P, ω = _default_freq_vector(P,Val{:bode}()); kd=0, form=
     phi = angle.(Pv)
     kp  = @. -cos(phi)/r
     ki  = @. kd*ω^2 - ω*sin(phi)/r
-    kp, ki  = convert_pidparams_from_to(kp, ki, kd, :parallel, form)
+    K   = convert_pidparams_from_parallel.(kp, ki, kd, form)
+    kp, ki = getindex.(K, 1), getindex.(K, 2)
     fig = if doplot
         RecipesBase.plot(kp,ki,linewidth = 1.5, xlabel=L"k_p", ylabel=L"k_i", title="Stability region of P, k_d = $(round(kd, digits=4))")
     else 
@@ -283,7 +282,8 @@ function stabregionPID(P::Function, ω = exp10.(range(-3, stop=1, length=50)); k
     phi     = angle.(Pv)
     kp      = -cos.(phi)./r
     ki      = @. kd*ω^2 - ω*sin(phi)/r
-    kp, ki  = convert_pidparams_from_to(kp, ki, kd, :parallel, form)
+    K       = convert_pidparams_from_parallel.(kp, ki, kd, form)
+    kp, ki  = getindex.(K, 1), getindex.(K, 2)
     fig = if doplot
         RecipesBase.plot(kp,ki,linewidth = 1.5, xlabel=L"k_p", ylabel=L"k_i", title="Stability region of P, k_d = $(round(kd, digits=4))")
     else 
@@ -300,7 +300,7 @@ Selects the parameters of a PI-controller (on parallel form) such that the Nyqui
 
 The parameters can be returned as one of several common representations 
 chosen by `form`, the options are
-* `:standard` - ``K_p(1 + 1/(T_i s) + T_ds)``
+* `:standard` - ``K_p(1 + 1/(T_i s) + T_d s)``
 * `:series` - ``K_c(1 + 1/(τ_i s))(τ_d s + 1)``
 * `:parallel` - ``K_p + K_i/s + K_d s``
 
@@ -350,7 +350,7 @@ function loopshapingPI(P0, ωp; ϕl=0, rl=0, phasemargin=0, form::Symbol=:standa
     else
         nothing
     end
-    kp, ki = convert_pidparams_from_to(kp, ki, 0, :parallel, form)
+    kp, ki = convert_pidparams_from_parallel(kp, ki, 0, form)
     (; C, kp, ki, fig, CF)
 end
 
@@ -491,7 +491,7 @@ function loopshapingPID(P0, ω; Mt = 1.3, ϕt=75, form::Symbol = :standard, dopl
     verbose && ki < 0 && @warn "Calculated ki is negative, inspect the Nyquist plot generated with doplot = true and try adjusting ω or the angle ϕt"
     C = pid(kp, ki, kd, form=:parallel)
     any(real(p) > 0 for p in poles(C)) && @error "Calculated controller is unstable."
-    kp, ki, kd = ControlSystemsBase.convert_pidparams_from_to(kp, ki, kd, :parallel, form)
+    kp, ki, kd = convert_pidparams_from_parallel(kp, ki, kd, form)
     CF = C*F
     fig = if doplot
         w = exp10.(LinRange(log10(ω)-2, log10(ω)+2, 1000))
@@ -522,15 +522,42 @@ The `form` can be chosen as one of the following
 """
 function convert_pidparams_to_standard(param_p, param_i, param_d, form::Symbol)
     if form === :standard
-        return param_p, param_i, param_d
+        return (param_p, param_i, param_d)
     elseif form === :series
-        return @. (
+        return (
             param_p * (param_i + param_d) / param_i,
             param_i + param_d,
             param_i * param_d / (param_i + param_d)
         )
     elseif form === :parallel
-        return @. (param_p, param_p / param_i, param_d / param_p)
+        return (param_p, param_p / param_i, param_d / param_p)
+    else
+        throw(ArgumentError("form $(form) not supported."))
+    end
+end
+
+"""
+    Kp, Ti, Td = convert_pidparams_to_parallel(param_p, param_i, param_d, form)
+
+Convert parameters from form `form` to `:parallel` form.
+
+The `form` can be chosen as one of the following
+* `:standard` - ``K_p(1 + 1/(T_i s) + T_d s)``
+* `:series` - ``K_c(1 + 1/(τ_i s))(τ_d s + 1)``
+* `:parallel` - ``K_p + K_i/s + K_d s``
+"""
+function convert_pidparams_to_parallel(param_p, param_i, param_d, form::Symbol)
+    if form === :parallel
+        return (param_p, param_i, param_d)
+    elseif form === :series
+        # param_i = 0 would result in division by zero, but typically indicates that the user wants no integral action
+        param_i == 0 && return (param_p, 0, param_p * param_d)
+        return (param_p * (param_i + param_d) / param_i,
+                param_p / param_i,
+                param_p * param_d)
+    elseif form === :standard
+        param_i == 0 && return (param_p, 0, param_p * param_d)
+        return (param_p, param_p / param_i, param_p * param_d)
     else
         throw(ArgumentError("form $(form) not supported."))
     end
@@ -542,30 +569,62 @@ end
 Convert parameters to form `form` from `:standard` form. 
 
 The `form` can be chosen as one of the following
-* `:standard` - ``K_p(1 + 1/(T_i s) + T_ds)``
+* `:standard` - ``K_p(1 + 1/(T_i s) + T_d s)``
 * `:series` - ``K_c(1 + 1/(τ_i s))(τ_d s + 1)``
 * `:parallel` - ``K_p + K_i/s + K_d s``
 """
 function convert_pidparams_from_standard(Kp, Ti, Td, form::Symbol)
     if form === :standard
-        return Kp, Ti, Td
+        return (Kp, Ti, Td)
     elseif form === :series
-        return @. (
-            (Ti - sqrt(Ti * (Ti - 4 * Td))) / 2 * Kp / Ti,
-            (Ti - sqrt(Ti * (Ti - 4 * Td))) / 2,
-            (Ti + sqrt(Ti * (Ti - 4 * Td))) / 2,
-        )
+        Δ = Ti * (Ti - 4 * Td)
+        Δ < 0 && throw(DomainError("The condition Ti^2 ≥ 4Td*Ti is not satisfied: the PID parameters cannot be converted to series form"))
+        sqrtΔ = sqrt(Δ)
+        return ((Ti - sqrtΔ) / 2 * Kp / Ti,
+                (Ti - sqrtΔ) / 2,
+                (Ti + sqrtΔ) / 2)
     elseif form === :parallel
-        return @. (Kp, Kp/Ti, Td*Kp)
+        return (Kp, Kp/Ti, Td*Kp)
+    else
+        throw(ArgumentError("form $(form) not supported."))
+    end
+end
+
+
+"""
+    Kp, Ti, Td = convert_pidparams_from_parallel(Kp, Ki, Kd, to_form)
+
+Convert parameters from form `:parallel` to form `to_form`.
+
+The `form` can be chosen as one of the following
+* `:standard` - ``K_p(1 + 1/(T_i s) + T_d s)``
+* `:series` - ``K_c(1 + 1/(τ_i s))(τ_d s + 1)``
+* `:parallel` - ``K_p + K_i/s + K_d s``
+"""
+function convert_pidparams_from_parallel(Kp, Ki, Kd, to::Symbol)
+    if to === :parallel
+        return (Kp, Ki, Kd)
+    elseif to === :series
+        Ki == 0 && return (Kp, 0, Kp*Kd)
+        Δ = Kp^2-4Ki*Kd
+        Δ < 0 &&
+            throw(DomainError("The condition Kp^2 ≥ 4Ki*Kd is not satisfied: the PID parameters cannot be converted to series form"))
+        sqrtΔ = sqrt(Δ)
+        return ((Kp - sqrtΔ)/2, (Kp - sqrtΔ)/(2Ki), (Kp + sqrtΔ)/(2Ki))
+    elseif to === :standard
+        Kp == 0 && throw(DomainError("Cannot convert to standard form when Kp=0"))
+        Ki == 0 && return (Kp, Inf, Kd / Kp)
+        return (Kp, Kp / Ki, Kd / Kp)
     else
         throw(ArgumentError("form $(form) not supported."))
     end
 end
 
+
 """
     convert_pidparams_from_to(kp, ki, kd, from::Symbol, to::Symbol)
 """
 function convert_pidparams_from_to(kp, ki, kd, from::Symbol, to::Symbol)
-    kp, ki, kd = convert_pidparams_to_standard(kp, ki, kd, from)
-    convert_pidparams_from_standard(kp, ki, kd, to)
+    Kp, Ki, Kd = convert_pidparams_to_parallel(kp, ki, kd, from)
+    convert_pidparams_from_parallel(Kp, Ki, Kd, to)
 end
diff --git a/lib/ControlSystemsBase/src/types/DelayLtiSystem.jl b/lib/ControlSystemsBase/src/types/DelayLtiSystem.jl
index 94119960c..f178d8c00 100644
--- a/lib/ControlSystemsBase/src/types/DelayLtiSystem.jl
+++ b/lib/ControlSystemsBase/src/types/DelayLtiSystem.jl
@@ -94,9 +94,16 @@ function +(sys::DelayLtiSystem{T1,S}, n::T2) where {T1,T2<:Number,S}
 end
 
 
+"""
+    /(G1, G2::DelayLtiSystem)
+
+Compute ``G_1 * G_2^{-1} where ``G_2`` is a DelayLtiSystem.
+Throws a SingularException if ``G_2`` is not invertible.
+"""
 function /(anything, sys::DelayLtiSystem)
-    all(iszero, sys.Tau) || error("A delayed system can not be inverted. Consider use of the function `feedback`.")
-    /(anything, sys.P.P) # If all delays are zero, invert the inner system
+    ny,nu = size(sys)
+    ny == nu || error("The denominator system must be square")
+    return anything * feedback(I(nu), sys - I(nu))
 end
 
 for other_type in [:Number, :AbstractMatrix, :LTISystem]
diff --git a/lib/ControlSystemsBase/src/types/LFTSystem.jl b/lib/ControlSystemsBase/src/types/LFTSystem.jl
index 1975261f3..630e47d6a 100644
--- a/lib/ControlSystemsBase/src/types/LFTSystem.jl
+++ b/lib/ControlSystemsBase/src/types/LFTSystem.jl
@@ -1,6 +1,10 @@
 abstract type LFTSystem{TE, T} <: LTISystem{TE} end
 timeevol(sys::LFTSystem) = timeevol(sys.P)
 
+Base.zero(sys::LFTSystem) = ss(zero(sys.P.D), sys.P.timeevol)
+Base.zero(::Type{<:LFTSystem{Continuous, F}}) where {F} = ss([zero(F)], Continuous())
+
+
 function *(n::Number, sys::LFTSystem)
     new_C = [sys.P.C1*n; sys.P.C2]
     new_D = [sys.P.D11*n sys.P.D12*n; sys.P.D21 sys.P.D22]
diff --git a/lib/ControlSystemsBase/test/test_delayed_systems.jl b/lib/ControlSystemsBase/test/test_delayed_systems.jl
index b8038d47a..a848b8ff4 100644
--- a/lib/ControlSystemsBase/test/test_delayed_systems.jl
+++ b/lib/ControlSystemsBase/test/test_delayed_systems.jl
@@ -56,6 +56,13 @@ P2_fr = (im*ω .+ 1) ./ (im*ω .+ 2)
 @test freqresp(P2 * delay(1), ω)[:] ≈ P2_fr .* exp.(-im*ω) rtol=1e-15
 @test freqresp(delay(1) * P2, ω)[:] ≈ P2_fr .* exp.(-im*ω) rtol=1e-15
 
+## Division / feedback
+@test freqresp(1/(1+P1), ω) ≈ freqresp(feedback(I(size(P1, 1)), P1), ω) rtol=1e-15
+
+## Zero
+@test zero(typeof(P1)) == ss(0.0)
+@test zero(P1) == ss(0.0)
+
 ## append
 P12 = append(P1, P2)
 G12 = [P1 tf(0); tf(0) P2]
@@ -179,9 +186,10 @@ w = 10 .^ (-2:0.1:2)
 
 @test propertynames(delay(1.0)) == (:P, :Tau, :nu, :ny)
 
-@test_throws ErrorException 1/f2
-@test_throws ErrorException randn(2,2)/f2
-@test_throws ErrorException f2/f2
+@test_throws SingularException 1/f2
+@test_throws SingularException randn(2,2)/f2
+@test_throws SingularException f2/f2
+@test 1/(I(2)+f2) == feedback(I(2), f2)
 
 #FIXME: A lot more tests, including MIMO systems in particular
 
diff --git a/lib/ControlSystemsBase/test/test_discrete.jl b/lib/ControlSystemsBase/test/test_discrete.jl
index 13f050b17..37d509418 100644
--- a/lib/ControlSystemsBase/test/test_discrete.jl
+++ b/lib/ControlSystemsBase/test/test_discrete.jl
@@ -146,6 +146,24 @@ p = poles(Gcl)
 pd = c2d_poly2poly(A, 0.1)
 @test pd ≈ denvec(c2d(P, 0.1))[]
 
+
+@testset "d2c_exact" begin
+    @info "Testing d2c_exact"
+    sys = ssrand(2, 3, 2, Ts = 1)
+    sysc_causal = d2c_exact(sys, :causal)
+    sysc_acausal = d2c_exact(sys, :acausal)
+    @test_throws ErrorException d2c_exact(sys, :KentMorgan)
+    # bodeplot(sys, w, lab="Discrete SS")
+    # bodeplot!(sysc_causal, w, lab="Continuous SS with causal delays", l=:dash, size=(800, 800))
+    # bodeplot!(sysc_acausal, w, lab="Continuous SS with acausal (negative) delays", l=:dash, size=(800, 800))
+    w = exp10.(LinRange(-2, 2, 200))
+    @test freqresp(sys, w) ≈ freqresp(sysc_causal, w) atol = 1e-8
+    @test freqresp(sys, w) ≈ freqresp(sysc_acausal, w) atol = 1e-8
+
+    # Requires full ControlSystems but does pass
+    # rd = step(sys, 0:10)
+    # rc = step(sysc_causal, 0:10)
+    # @test rd.y ≈ rc.y atol = 1e-8
 end
 
 
@@ -233,3 +251,4 @@ Qd_van_load = [
 @test norm(Qd - Qd_van_load) < 1e-6
 
 
+end
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/test/test_pid_design.jl b/lib/ControlSystemsBase/test/test_pid_design.jl
index 6723884e9..c9df2fb9b 100644
--- a/lib/ControlSystemsBase/test/test_pid_design.jl
+++ b/lib/ControlSystemsBase/test/test_pid_design.jl
@@ -1,5 +1,7 @@
 @testset "test_pid_design" begin
 
+CSB = ControlSystemsBase
+
 # Test gof plot and loopshaping
 P = tf(1,[1,1])^4
 gangoffourplot(P,tf(1))
@@ -13,8 +15,18 @@ C, kp, ki = loopshapingPI(P, ωp, phasemargin=60, form=:parallel, doplot=true)
 # tf
 @test pid(1.0, 1, 1) == tf(1) + tf(1,[1,0]) + tf([1,0],[1])
 @test pid(1.0, Inf, 1) == tf(1) + tf([1, 0], [1])
+@test pid(1.0, 0, 1) == tf(1) + tf([1, 0], [1])
+@test pid(0.0, 1, 1; form=:parallel) == tf(0) + tf(1,[1,0]) + tf([1,0],[1])
+@test pid(1.0, 2, 3; Tf=2) == tf([3,1,0.5], [2,2,1,0])
+@test all(CSB.convert_pidparams_from_standard(CSB.convert_pidparams_from_parallel(1, 2, 3, :standard)...,
+                                                  :parallel) .≈ (1,2,3))
+@test_throws DomainError CSB.convert_pidparams_from_parallel(2, 3, 0.5, :series)
+@test_throws DomainError CSB.convert_pidparams_from_parallel(0, 3, 0.5, :standard)
+@test_throws DomainError CSB.convert_pidparams_from_standard(2, 1, 0.5, :series)
 # ss
 @test tf(pid(1.0, 1, 0; state_space=true)) == tf(1) + tf(1,[1,0])
+@test tf(pid(0.0, 2, 3; form=:parallel, state_space=true, Tf=2)) == tf([3,0,2], [2, 2, 1, 0])
+@test tf(pid(1.0, 2, 3; state_space=true, Tf=2)) == tf([3, 1, 0.5], [2, 2, 1, 0])
 
 # Discrete
 @test_throws ArgumentError pid(1.0, 1, 1, Ts=0.1)
@@ -72,6 +84,7 @@ C, Kp, Ti = placePI(P, 2, 0.7; form=:standard)
 @test [Kp, Ti] ≈ [9/5, 9/20]
 
 # Test internal functions convert_pidparams*
+# Standard
 params = (2, 3, 0.5)
 parallel_params = ControlSystemsBase.convert_pidparams_from_standard(params..., :parallel)
 @test parallel_params == (2, 2/3, 1)
@@ -80,6 +93,15 @@ series_params = ControlSystemsBase.convert_pidparams_from_standard(params..., :s
 @test series_params == ((3-sqrt(3))/3, (3-sqrt(3))/2, (3+sqrt(3))/2)
 @test ControlSystemsBase.convert_pidparams_to_standard(series_params..., :series) == params
 
+# Parallel
+params = (4, 3, 0.5)
+standard_params = ControlSystemsBase.convert_pidparams_from_parallel(params..., :standard)
+@test standard_params == (4, 4/3, 0.5/4)
+@test ControlSystemsBase.convert_pidparams_to_parallel(standard_params..., :standard) == params
+series_params = ControlSystemsBase.convert_pidparams_from_parallel(params..., :series)
+@test series_params == ((4-sqrt(10))/2, (4-sqrt(10))/6, (4+sqrt(10))/6)
+@test all(ControlSystemsBase.convert_pidparams_to_parallel(series_params..., :series) .≈ params)
+
 # lead lag link
 a = 1
 M = 10