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