add example with uncertain delay (#54)

* add example with uncertain delay

* add more examples and tests

docs/src/uncertainty.md

@@ -339,3 +339,36 @@ TODO
 [^Skogestad]: Skogestad, "Multivariable Feedback Control: Analysis and Design"
 
 [^Doyle91]: Doyle, Packard, Zhou, "Review of LFTs, LMIs and μ". [`https://www.researchgate.net/publication/257200344_Review_of_LFTs_LMIs_and_mu`](https://www.researchgate.net/publication/257200344_Review_of_LFTs_LMIs_and_mu)
+
+
+## Uncertain time delays
+
+Modeling uncertain time delays can be done in several ways, one approach is to make use of a multiplicative uncertainty weight created using [`neglected_delay`](@ref) multiplied by an uncertain element created using [`δr`](@ref) or [`δc`](@ref), example:
+```@example uncertain_delay
+using RobustAndOptimalControl, ControlSystems, MonteCarloMeasurements
+a  = 10
+P  = ss([0 a; -a 0], I(2), [1 a; -a 1], 0) # Plant
+W0 = neglected_delay(0.005) |> ss # Weight
+W  = I(2) + W0*I(2) * uss([δr(), δr()]) # Create a diagonal real uncertainty weighted in frequency by W0
+Ps = P*W # Uncertain plant
+Psamples = rand(Ps, 100) # Sample the uncertain plant for plotting
+w = exp10.(LinRange(-1, 3, 300)) # Frequency vector
+bodeplot(Psamples, w)
+```
+Note how this approximation approach imparts some uncertainty also in the gain.
+
+The other alternative is to use use sampled uncertain delays. The next example shows how we can create a system with an uncertain delay, where we know that the delay is an integer number of milliseconds between 1ms and 4ms.
+```@example uncertain_delay
+using RobustAndOptimalControl, ControlSystems, MonteCarloMeasurements
+unsafe_comparisons(true)
+L = Particles(collect((1:4) ./ 1000)) # Uncertain time delay, an integer number of milliseconds between 1ms and 4ms
+P = delay(L)*tf(1, [0.01, 1])
+C = pid(kp=2, ki=1, series=true)
+w = exp10.(-1:0.01:4)
+plot(
+     bodeplot(P, exp10.(-1:0.001:3)),
+     plot(step(feedback(P, C), 0:0.0001:0.05), lab="L = " .* string.(P.Tau[].particles'), title="Disturbance response"),
+     nyquistplot(P*C, w[1:10:end], points=true, xlims=(-3.5, 2.5), ylims=(-5, 1.5), Ms_circles=[1.5, 2], alpha=1) # Note, the nyquistplot with uncertain coefficients requires manual selection of plot limits
+)
+```
+Notice how the gain is completely certain, while the phase starts becoming very uncertain for high frequencies. The phase unwrapping for the `bodeplot` does not handle the uncertainty correctly, so interpret the phase curve with care. 

src/uncertainty_interface.jl

@@ -464,6 +464,25 @@ function sys_from_particles(P, i)
     ss(vecindex(A, i), vecindex(B, i), vecindex(C, i), vecindex(D, i))
 end
 
+function sys_from_particles(P::DelayLtiSystem{T,S}, i) where {T, S<:AbstractParticles}
+    A,B,C,D = ssdata(ss(P.P.P))
+    ControlSystems.DelayLtiSystem(
+        ControlSystems.PartitionedStateSpace(ss(vecindex(A, i), vecindex(B, i), vecindex(C, i), vecindex(D, i)), P.P.nu1, P.P.ny1),
+        vecindex(P.Tau, i)
+    )
+end
+
+function sys_from_particles(P)
+    map(1:nparticles(ControlSystems.numeric_type(P))) do i
+        sys_from_particles(P, i)
+    end
+end
+
+
+function ControlSystems.lsim(sys::DelayLtiSystem{T,S}, u, t::AbstractArray{<:Real}, args...; x0=fill(zero(T), nstates(sys)), kwargs...) where {T, S<:AbstractParticles}
+    syss = sys_from_particles(sys)
+    [lsim(sysi, u, t, args...; x0 = vecindex(x0, i), kwargs...) for (i, sysi) in enumerate(syss)]
+end
 
 # function any0det(D::Matrix{<:Complex{<:AbstractParticles}})
 #     D0 = similar(D, ComplexF64)

test/test_uncertainty.jl

@@ -330,4 +330,15 @@ dmm = argmin(dm->dm.margin, dm.simultaneous_input)
 @test dmm.margin ≈ 0.5335 atol = 0.001
 
 
+## Uncertain delays
+using RobustAndOptimalControl, ControlSystems, MonteCarloMeasurements
+unsafe_comparisons(true)
+L = Particles(collect((1:4) ./ 1000)) # Uncertain time delay, an integer number of milliseconds between 1ms and 4ms
+P = delay(L)*tf(1, [0.01, 1])
+C = pid(kp=2, ki=1, series=true)
+w = exp10.(-1:0.01:4)
+
+bodeplot(P, exp10.(-1:0.001:3))
+plot(step(feedback(P, C), 0:0.0001:0.05), lab="L = " .* string.(P.Tau[].particles'), title="Disturbance response")
+nyquistplot(P*C, w[1:10:end], points=true, xlims=(-3.5, 2.5), ylims=(-5, 1.5), Ms_circles=[1.5, 2], alpha=1) # Note, the nyquistplot with uncertain coefficients requires manual selection of plot limits