-
Notifications
You must be signed in to change notification settings - Fork 17
/
solve.jl
412 lines (329 loc) · 13.9 KB
/
solve.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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
# ======================================
# Solve function for continuous systems
# ======================================
const AbstractContinuousSystem_ = AbstractContinuousSystem # to trick the linter
"""
solve(ivp::IVP{<:AbstractContinuousSystem}, tspan, alg; kwargs...)
Solves the initial-value problem defined by `ivp` over the time span `tspan`,
using the algorithm `alg`. If no algorithm is given, a default algorithm is chosen.
### Input
- `ivp` -- initial-value problem
- `tspan` -- time span for this initial-value problem
- `alg` -- reachability algorithm
Additional options are passed as arguments or keyword arguments; see the notes
below for details. See the online documentation for examples.
### Output
The solution of a reachability problem, as an instance of a `ReachSolution`.
### Notes
- Use the `alg`, `algorithm` or `opC` keyword arguments to specify the algorithm
to solve the initial-value problem. Algorithm-specific options should be passed
to the algorithm constructor as well.
- Use the `tspan` keyword argument to specify the time span; it can be:
- a tuple,
- an interval, or
- a vector with two components.
- Use the `T` keyword argument to specify the time horizon; the initial time is
then assumed to be zero.
- Use the `static` keyword argument to force conversion to static arrays in the
algorithm (should be supported by the algorithm).
- Use the `NSTEPS` keyword argument to specify the number of discrete steps
solved in the set-based recurrence.
- Use the `threading` option to use multi-threading parallelism. This option applies
for initial-value problems whose initial condition is a vector of sets.
"""
function solve(ivp::IVP{<:AbstractContinuousSystem_}, args...; kwargs...)
# preliminary checks
_check_dim(ivp)
# get time span (or the emptyset if NSTEPS was specified)
tspan = _get_tspan(args...; kwargs...)
# get the continuous post or find a default one
cpost = _get_cpost(ivp, args...; kwargs...)
if isnothing(cpost)
cpost = _default_cpost(ivp, tspan; kwargs...)
end
# run the continuous-post operator
F = post(cpost, ivp, tspan; kwargs...)
# optionally compute ensemble simulations
got_ensemble = get(kwargs, :ensemble, false)
dict = Dict{Symbol,Any}(:ensemble => nothing)
if got_ensemble
@requires OrdinaryDiffEq
ensemble_sol = _solve_ensemble(ivp, args...; kwargs...)
dict[:ensemble] = ensemble_sol
end
return _solve_return(ivp, F, cpost, dict, args...; kwargs...)
end
# wrap the flowpipe and algorithm in a solution structure
function _solve_return(ivp::IVP{<:AbstractContinuousSystem}, F, cpost, dict, args...; kwargs...)
return sol = ReachSolution(F, cpost, dict)
end
# solve for distributed initial conditions; uses multi-threaded implementation by default
function solve(ivp::IVP{AT,VT}, args...;
kwargs...) where {AT<:AbstractContinuousSystem,
ST<:Union{LazySet,IntervalBox},VT<:AbstractVector{ST}}
_check_dim(ivp)
tspan = _get_tspan(args...; kwargs...)
cpost = _get_cpost(ivp, args...; kwargs...)
if isnothing(cpost)
cpost = _default_cpost(ivp, tspan; kwargs...)
end
X0 = initial_state(ivp)
S = system(ivp)
threading = haskey(kwargs, :threading) ? kwargs[:threading] : (Threads.nthreads() > 1)
F = _solve_distributed(cpost, S, X0, tspan, Val(threading); kwargs...)
return ReachSolution(MixedFlowpipe(F), cpost)
end
function _solve_distributed(cpost, S, X0, tspan, threading::Val{false}; kwargs...)
return [post(cpost, IVP(S, X0i), tspan; kwargs...) for X0i in X0]
end
function _solve_distributed(cpost, S, X0, tspan, threading::Val{true}; kwargs...)
nsets = length(X0)
FT = Flowpipe{numtype(cpost),rsetrep(cpost)} # TODO add third parameter
sol_tot = Vector{FT}(undef, nsets)
Threads.@threads for i in eachindex(X0)
sol_i = ReachabilityAnalysis.post(cpost, IVP(S, X0[i]), tspan; kwargs...)
sol_tot[i] = sol_i
end
return sol_tot
end
# ====================================
# Continuous post operator interface
# ====================================
"""
AbstractPost
Abstract supertype of all post operator types.
"""
abstract type AbstractPost end
"""
AbstractContinuousPost
Abstract supertype of all continuous post operators.
"""
abstract type AbstractContinuousPost <: AbstractPost end
setrep(::InitialValueProblem{HS,ST}) where {HS,ST} = ST
# ==================
# Argument handling
# ==================
# extend dimension for common IntervalArithmetic types
LazySets.dim(::IA.Interval) = 1
LazySets.dim(::IA.IntervalBox{D,N}) where {D,N} = D
# lazy sets (or sets that behave as such)
_dim(X::Union{<:LazySet,<:IA.Interval,<:IA.IntervalBox}) = dim(X)
# singleton elements
_dim(X::Number) = 1
_dim(X::AbstractVector{N}) where {N<:Number} = length(X)
# vector of sets
function _dim(X::AbstractVector{UT}) where {UT<:Union{<:LazySet,
<:IA.Interval,<:IA.IntervalBox}}
n = _dim(first(X))
all(X -> _dim(X) == n, X) || throw(ArgumentError("dimension mismatch between " *
"the initial sets in this array; expected only sets of dimension $n"))
return n
end
_dim(X) = throw(ArgumentError("the type of the initial condition, $(typeof(X)), cannot be handled"))
function _check_dim(S, X0; throw_error::Bool=true)
n = statedim(S)
d = _dim(X0)
n == d && return true
if throw_error
throw(ArgumentError("the state-space dimension should match " *
"the dimension of the initial state, but they are " *
"$n and $d respectively"))
end
return false
end
_dim(X0::Tuple{VT,VT}) where {N,VT<:AbstractVector{N}} = length(X0[1]) + length(X0[2])
_dim(X0::Tuple{<:LazySet{N},<:LazySet{N}}) where {N} = dim(X0[1]) + dim(X0[2])
function _check_dim(S::Union{SecondOrderLinearContinuousSystem,
SecondOrderAffineContinuousSystem,
SecondOrderConstrainedLinearControlContinuousSystem,
SecondOrderConstrainedAffineControlContinuousSystem},
X0; throw_error::Bool=true)
n = statedim(S)
d = _dim(X0)
d == 2 * n && return true
if throw_error
throw(ArgumentError("the dimension of the initial state should be " *
"twice the state-space dimension, but they are " *
"$d and $n respectively"))
end
return false
end
function _check_dim(ivp::InitialValueProblem{<:AbstractContinuousSystem}; throw_error::Bool=true)
S = system(ivp)
X0 = initial_state(ivp)
return _check_dim(S, X0; throw_error=throw_error)
end
# promotion from tuple-like arguments
@inline _promote_tspan((t1, t2)::Tuple{T,T}) where {T} = TimeInterval(t1, t2)
@inline _promote_tspan((t1, t2)::Tuple{T,S}) where {T,S} = TimeInterval(promote(t1, t2))
# no-op, corresponds to (inf(tspan), sup(tspan))
@inline _promote_tspan(tspan::IA.Interval) = tspan
# no-op, takes interval wrapped data; corresponds to (min(tspan), max(tspan))
@inline _promote_tspan(tspan::Interval) = tspan.dat
# number T defaults to a time interval of the form [0, T]
@inline _promote_tspan(tspan::Number) = TimeInterval(zero(tspan), tspan)
# catch-all array like tspan
@inline function _promote_tspan(tspan::AbstractArray)
@assert length(tspan) == 2 throw(ArgumentError("the length of `tspan` must be two, but " *
"it is of length $(length(tspan))"))
return TimeInterval(tspan[1], tspan[2])
end
function _get_tspan(args...; kwargs...)
got_tspan = haskey(kwargs, :tspan)
got_T = haskey(kwargs, :T)
got_NSTEPS = haskey(kwargs, :NSTEPS)
# throw error for repeated arguments
if got_tspan && got_T
throw(ArgumentError("cannot parse the time horizon `T` and the " *
"time span `tspan` simultaneously; use one or the other"))
elseif got_tspan && got_NSTEPS
throw(ArgumentError("cannot parse the time span `tspan` and the " *
"number of steps `NSTEPS` simultaneously; use one or the other"))
elseif got_T && got_NSTEPS
throw(ArgumentError("cannot parse the time horizon `T` and the " *
"number of steps `NSTEPS` simultaneously; use one or the other"))
end
# parse time span
if got_tspan
tspan = _promote_tspan(kwargs[:tspan])
elseif got_T
tspan = _promote_tspan(kwargs[:T])
elseif got_NSTEPS
tspan = emptyinterval() # defined a posteriori
elseif !isempty(args) && applicable(_promote_tspan, args[1])
# applies to tuples, vectors, LazySets.Interval, IA.Interval, and numbers
tspan = _promote_tspan(args[1])
#=
(args[1] isa Tuple{Float64, Float64} || # time span given as tuple
args[1] isa Vector{Float64} || # time span given as vector
args[1] isa Interval || # time span given as LazySets.Interval
args[1] isa IA.Interval || # time span given as IA.Interval
args[1] isa Real) # got time horizon as first argument
=#
elseif haskey(kwargs, :max_jumps)
# got maximum number of jumps, applicable to hybrid systems
tspan = emptyinterval()
else
# couldn't find time span => error
throw(ArgumentError("the time span has not been specified, but is required " *
"for `solve`; you should specify either the time horizon " *
"`T=...`, the time span `tspan=...`, or the number of steps, `NSTEPS=...`"))
end
return tspan
end
# return the time horizon given a time span
# the check_positive flag is used for algorithms that do not support negative
# times
function _get_T(tspan::TimeInterval; check_zero::Bool=true, check_positive::Bool=true)
t0 = inf(tspan)
if check_zero
@assert iszero(t0) "this algorithm can only handle zero initial time"
end
T = sup(tspan)
if check_positive
@assert T > 0 "the time horizon should be positive"
end
return T
end
tstart(Δt::TimeInterval) = inf(Δt)
tend(Δt::TimeInterval) = sup(Δt)
tspan(Δt::TimeInterval) = Δt
function _get_cpost(ivp, args...; kwargs...)
got_alg = haskey(kwargs, :alg)
got_algorithm = haskey(kwargs, :algorithm)
got_opC = haskey(kwargs, :opC)
no_args = isempty(args) || args[1] === nothing
# continuous post was specified
if got_alg
cpost = kwargs[:alg]
elseif got_algorithm
cpost = kwargs[:algorithm]
elseif got_opC
cpost = kwargs[:opC]
# check if either args[1] or args[2] are the post
elseif !no_args && typeof(args[1]) <: AbstractContinuousPost
cpost = args[1]
elseif !no_args && length(args) == 2 && typeof(args[2]) <: AbstractContinuousPost
cpost = args[2]
elseif !no_args && length(args) > 2
throw(ArgumentError("the number of arguments, $(length(args)), is not valid"))
# no algorithm found => compute default
else
cpost = nothing
end
return cpost
end
# ===================
# Algorithm defaults
# ===================
# number of discrete steps used if only the time horizon is given
const DEFAULT_NSTEPS = 100
function _default_cpost(S::AbstractContinuousSystem, X0, ishybrid, tspan; kwargs...)
if isaffine(S)
if haskey(kwargs, :δ)
δ = kwargs[:δ]
elseif haskey(kwargs, :N)
N = kwargs[:N]
δ = diam(tspan) / N
elseif haskey(kwargs, :NSTEPS)
NSTEPS = kwargs[:NSTEPS]
δ = diam(tspan) / NSTEPS
elseif haskey(kwargs, :num_steps)
num_steps = kwargs[:num_steps]
δ = diam(tspan) / num_steps
else
δ = diam(tspan) / DEFAULT_NSTEPS
end
if statedim(S) == 1 && !is_second_order(S)
opC = INT(; δ=δ)
else
static = haskey(kwargs, :static) ? kwargs[:static] : false
if ishybrid || !(X0 isa AbstractZonotope)
opC = GLGM06(; δ=δ, static=static, approx_model=Forward())
else
opC = GLGM06(; δ=δ, static=static)
end
end
else
# check additional kwargs options if they exist, allowing some aliases
orderQ = get(kwargs, :orderQ, DEFAULT_ORDER_Q_TMJETS)
orderT = get(kwargs, :orderT, DEFAULT_ORDER_T_TMJETS)
if haskey(kwargs, :max_steps)
maxsteps = kwargs[:max_steps]
elseif haskey(kwargs, :maxsteps)
maxsteps = kwargs[:maxsteps]
else
maxsteps = DEFAULT_MAX_STEPS_TMJETS
end
if haskey(kwargs, :abs_tol)
abstol = kwargs[:abs_tol]
elseif haskey(kwargs, :abstol)
abstol = kwargs[:abstol]
else
abstol = DEFAULT_ABS_TOL_TMJETS
end
opC = TMJets(; orderQ=orderQ, orderT=orderT, maxsteps=maxsteps, abstol=abstol)
end
return opC
end
"""
_default_cpost(ivp::IVP{<:AbstractContinuousSystem}, tspan; kwargs...)
### Input
- `ivp` -- initial-value problem
- `tspan` -- time-span
### Output
A continuous post-operator that can handle the given initial-value problem.
### Notes
If the system is affine, then:
- If it is one-dimensional, algorithm `INT` is used, otherwise,
- Algorithm `GLGM06` is used.
If the system is not affine, then the algorithm `TMJets` is used.
"""
function _default_cpost(ivp::IVP{<:AbstractContinuousSystem}, tspan; kwargs...)
return _default_cpost(system(ivp), initial_state(ivp), false, tspan; kwargs...)
end
# return a default algorithm, assuming that the first mode is representative
function _default_cpost(ivp::IVP{<:AbstractHybridSystem}, tspan; kwargs...)
first_mode = mode(system(ivp), 1)
return _default_cpost(first_mode, initial_state(ivp), true, tspan; kwargs...)
end