generated from JuliaDynamics/DynamicalSystemsBase.jl
-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathoptimized_shooting.jl
99 lines (82 loc) · 3.12 KB
/
optimized_shooting.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
export periodic_orbit, periodic_orbits, OptimizedShooting
using NonlinearSolve
"""
OptimizedShooting(; kwargs...)
A shooting method [Dednam2014](@cite) that uses Levenberg-Marquardt optimization
to find periodic orbits of continuous-time dynamical systems.
## Keyword arguments
- `Δt::Float64 = 1e-6`: step between the points in the residual `R`. See below for details.
- `n::Int64 = 2`: `n*dimension(ds)` is the number of points in the residual `R`. See below
for details.
- `nonlinear_solve_kwargs = (reltol=1e-6, abstol=1e-6, maxiters=1000)`: keyword arguments
to pass to the `solve` function from
[`NonlinearSolve.jl`](https://github.com/SciML/NonlinearSolve.jl). For details on the
keywords see the respective package documentation.
## Description
Let us consider the following continuous-time dynamical system
```math
\\frac{dx}{dt} = f(x, p, t)
```
Dednam and Botha [Dednam2014](@cite) suggest to minimize the residual ``R`` defined as
```math
R = (x(T)-x(0), x(T+\\Delta t)-x(\\Delta t), \\dots,
x(T+(n-1)\\Delta t)-x((n-1)\\Delta t))
```
where ``T`` is unknown period of a periodic orbit and ``x(t)`` is a solution at time ``t``
given some unknown initial point. Initial guess of the period ``T`` and the initial point
is optimized by the Levenberg-Marquardt algorithm.
In our implementation, the keyword argument `n` corresponds to ``n`` in the residual ``R``.
The keyword argument `Δt` corresponds to ``\\Delta t`` in the residual ``R``.
"""
@kwdef struct OptimizedShooting{T} <: PeriodicOrbitFinder
Δt::Float64 = 1e-6
n::Int64 = 2
nonlinear_solve_kwargs::T = (reltol=1e-6, abstol=1e-6, maxiters=1000)
end
function periodic_orbit(ds::CoupledODEs, alg::OptimizedShooting, ig::InitialGuess)
D = dimension(ds)
f = (err, v, p) -> begin
if isinplace(ds)
u0 = @view v[1:D]
else
u0 = SVector{D}(v[1:D])
end
T = v[end]
bounds = zeros(eltype(v), alg.n*2)
for i in 0:alg.n-1
bounds[i+1] = i*alg.Δt
bounds[i+alg.n+1] = T + i*alg.Δt
end
tspan = (0.0, T + alg.n*alg.Δt)
sol = solve(SciMLBase.remake(ds.integ.sol.prob; u0=u0,
tspan=tspan); DynamicalSystemsBase.DEFAULT_DIFFEQ..., ds.diffeq..., saveat=bounds)
if (length(sol.u) == alg.n*2)
for i in 1:alg.n
err[D*i-(D-1):D*i] = (sol.u[i] - sol.u[i+alg.n])
end
else
fill!(err, Inf)
end
end
prob = NonlinearLeastSquaresProblem(
NonlinearFunction(f, resid_prototype = zeros(alg.n*dimension(ds))), [ig.u0..., ig.T])
sol = solve(prob, NonlinearSolve.LevenbergMarquardt(); alg.nonlinear_solve_kwargs...)
if sol.retcode == ReturnCode.Success
u0 = sol.u[1:end-1]
T = sol.u[end]
Δt = 0.1
return PeriodicOrbit(ds, u0, T, Δt)
else
return nothing
end
end
function periodic_orbits(ds::CoupledODEs, alg::OptimizedShooting, igs::Vector{<:InitialGuess})
pos = []
for ig in igs
res = periodic_orbit(ds, alg, ig)
if !isnothing(res)
push!(pos, res)
end
end
return pos
end