-
Notifications
You must be signed in to change notification settings - Fork 17
/
TMJets.jl
108 lines (91 loc) · 3.83 KB
/
TMJets.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
@testset "TMJets algorithm (TMJets21b)" begin
prob, tspan = vanderpol()
# default algorithm for nonlinear systems
sol = solve(prob; tspan=tspan)
@test sol.alg isa TMJets
# pass the algorithm explicitly
sol = solve(prob; tspan=tspan, alg=TMJets())
@test sol.alg isa TMJets
# pass options outside algorithm
sol = solve(prob; T=0.2, maxsteps=2100)
@test sol.alg isa TMJets && sol.alg.maxsteps == 2100
sol = solve(prob; T=0.2, maxsteps=2100) # alias
@test sol.alg isa TMJets && sol.alg.maxsteps == 2100
sol = solve(prob; T=0.2, orderT=7, orderQ=1)
@test sol.alg isa TMJets && sol.alg.orderT == 7 && sol.alg.orderQ == 1
sol = solve(prob; T=0.2, abstol=1e-11)
@test sol.alg isa TMJets && sol.alg.abstol == 1e-11
sol = solve(prob; T=0.2, abstol=1e-11)
@test sol.alg isa TMJets && sol.alg.abstol == 1e-11
@test ReachabilityAnalysis.tspan(shift(sol, 1.0)) == ReachabilityAnalysis.tspan(sol) + 1.0
# split initial conditions
X0, S = initial_state(prob), system(prob)
X0s = split(X0, [2, 1]) # split along direction x
sols = solve(IVP(S, X0s); T=0.1, threading=false) # no threading
@test flowpipe(sols) isa MixedFlowpipe
sols = solve(IVP(S, X0s); T=0.1, threading=true) # with threading (default)
@test flowpipe(sols) isa MixedFlowpipe
end
@testset "TMJets algorithm (TMJets21b): linear IVPs" begin
prob, dt = exponential_1d()
sol = solve(prob; tspan=dt, alg=TMJets())
@test sol.alg isa TMJets
# getter functions for a taylor model reach-set
R = sol[1]
@test domain(R) == tspan(R)
@test diam(remainder(R)[1]) < 1e-9
@test get_order(R) == [8]
@test polynomial(R) isa Vector{Taylor1{TaylorN{Float64}}}
@test expansion_point(R) ≈ [IntervalArithmetic.interval(0.0)]
# test intersection with invariant
prob, dt = exponential_1d(; invariant=HalfSpace([-1.0], -0.3)) # x >= 0.3
sol_inv = solve(prob; tspan=dt, alg=TMJets())
@test [0.3] ∈ overapproximate(sol_inv[end], Zonotope)
m = length(sol_inv)
# check that the following reach-set escapes the invariant
@test [0.3] ∉ overapproximate(sol[m + 1], Zonotope)
# TODO test higher order system
# prob, tspan = linear5D_homog()
# sol = solve(prob, tspan=tspan, alg=TMJets())
# @test sol.alg isa TMJets
# TODO test linear system with input
# prob, tspan = linear5D()
# sol = solve(prob, tspan=tspan, alg=TMJets())
# @test sol.alg isa TMJets
end
@testset "TMJets algorithm (TMJets21a): linear IVPs" begin
prob, dt = exponential_1d()
sol = solve(prob; tspan=dt, alg=TMJets21a())
@test sol.alg isa TMJets21a
# getter functions for a taylor model reach-set
R = sol[1]
@test domain(R) == tspan(R)
@test diam(remainder(R)[1]) < 1e-13
@test get_order(R) == [8]
@test polynomial(R) isa Vector{Taylor1{TaylorN{Float64}}}
@test expansion_point(R) ≈ [IntervalArithmetic.interval(0.0)]
# test intersection with invariant
prob, dt = exponential_1d(; invariant=HalfSpace([-1.0], -0.3)) # x >= 0.3
sol_inv = solve(prob; tspan=dt, alg=TMJets21a())
@test [0.3] ∈ overapproximate(sol_inv[end], Zonotope)
m = length(sol_inv)
# check that the following reach-set escapes the invariant
@test [0.3] ∉ overapproximate(sol[m + 1], Zonotope)
end
@testset "1D Burgers equation (TMJets21b)" begin
L0 = 1.0 # domain length
U0 = 1.0 # Re = 20.
x = range(-0.5 * L0, 0.5 * L0; length=4)
# Initial velocity
X0 = Singleton(-U0 * sin.(2 * π / L0 * x))
# IVP definition
prob = @ivp(x' = burgers!(x), dim = 4, x(0) ∈ X0)
sol = solve(prob; tspan=(0.0, 1.0), alg=TMJets21b())
@test dim(sol) == 4
end
#=
alg = TMJets(abstol=1e-10, orderT=10, orderQ=2)
# reach mode
sol = solve(P, T=7.0, alg=alg)
@test set(sol[1]) isa Hyperrectangle # check default set representation
=#