add tutorial analyzing hybrid systems

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",

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