Plot nominal dynamics for Van der Pol system

Author: Zinoex

examples/Project.toml

@@ -8,6 +8,7 @@ IntervalSySCoRe = "6214ee8e-79dc-449b-ae00-4c9ff82a02a9"
 LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

examples/systems/van_der_pol.jl

@@ -1,6 +1,6 @@
 using Revise
 
-using LinearAlgebra, LazySets
+using LinearAlgebra, LazySets, Plots
 using IntervalMDP, IntervalSySCoRe
 
 
@@ -37,6 +37,31 @@ function van_der_pol_decoupled(; state_split=(50, 50), input_split=10)
     return mdp, reach, avoid
 end
 
+function van_der_pol_plot_nominal()
+    sys = van_der_pol_sys()
+
+    X = Hyperrectangle(low=[-4.0, -4.0], high=[4.0, 4.0])
+    state_abs = StateUniformGridSplit(X, (50, 50))
+
+    U = Hyperrectangle(low=[-1.0], high=[1.0])
+    input_abs = InputLinRange(U, 10)
+
+    R = IntervalSySCoRe.regions(state_abs)[837]
+    u = IntervalSySCoRe.inputs(input_abs)[3]
+
+    Y = nominal(dynamics(sys), R, u)
+
+    xs = sample(R, 10)
+    ys = [nominal(dynamics(sys), x, element(u)) for x in xs]
+
+    p = plot(R, color=:blue, label="X")
+    scatter!(p, Plots.unzip(Tuple.(xs)), color=:blue, label="Xhat", markershape=:xcross)
+    plot!(p, Y, color=:red, label="Y")
+    scatter!(p, Plots.unzip(Tuple.(ys)), color=:red, label="Yhat", markershape=:xcross)
+
+    display(p)
+end
+
 function main()
     @time "abstraction" mdp, reach, avoid = van_der_pol_decoupled(; state_split=(200, 200), input_split=21)
 

src/dynamics/NonlinearAdditiveNoiseDynamics.jl

@@ -90,7 +90,9 @@ function nominal(dyn::NonlinearAdditiveNoiseDynamics, X::Hyperrectangle, U::Hype
     return Yconv
 end
 
-function nominal(dyn::NonlinearAdditiveNoiseDynamics, X::Hyperrectangle, U::Singleton)
+nominal(dyn::NonlinearAdditiveNoiseDynamics, X::Hyperrectangle, U::Singleton) = nominal(dyn, X, element(U))
+
+function nominal(dyn::NonlinearAdditiveNoiseDynamics, X::Hyperrectangle, u::AbstractVector)
     # Use the Taylor model to over-approximate the reachable set
     order = 1
 
@@ -102,7 +104,6 @@ function nominal(dyn::NonlinearAdditiveNoiseDynamics, X::Hyperrectangle, U::Sing
     set_variables(Float64, "x"; order=2order, numvars=dimstate(dyn))
 
     x = [TaylorModelN(i, order, IntervalBox(x0), dom) for i in 1:dimstate(dyn)]
-    u = element(U)
 
     # Perform the Taylor expansion
     y = dyn.f(x, u)
@@ -127,7 +128,17 @@ function nominal(dyn::NonlinearAdditiveNoiseDynamics, X::Hyperrectangle, U::Sing
     return Yconv
 end
 
+function nominal(dyn::NonlinearAdditiveNoiseDynamics, X::Singleton, U::Singleton)
+    x = element(X)
+    u = element(U)
+
+    y = dyn.f(x, u)
+
+    return Singleton(y)
+end
+
 
+nominal(dyn::NonlinearAdditiveNoiseDynamics, x::AbstractVector, u::AbstractVector) = dyn.f(x, u)
 
 noise(dyn::NonlinearAdditiveNoiseDynamics) = dyn.w
 dimstate(dyn::NonlinearAdditiveNoiseDynamics) = dyn.nstate