Newer
Older
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
"""
`freqamp` frequency perturbation amplitude
`n` number of perturbations
Example
≡≡≡≡≡≡≡≡
julia>
J = 1 # only node 1
J = [1,2,3] # set of nodes
"""
function stability_problem(p, J; freqamp = 1.0, t_final=120.0, tolerance = 0.1)
# default ode problem
u_sync = synchronous_state(p)
ode_prob = ODEProblem(swing_dynamics!, u_sync, (0, t_final), p)
# low discrepancy sequence for unit inverval
N = nv(p.grid)
m = length(J)
P = primes(4*m^2+1)[1:2m]
H = Halton.(P)
# perturb intial condition
function _bstab_init(prob, i, repeat)
u0 = u_sync[:]
Z = getindex.(H, i) .* 2 .- 1
u0[J] .= Z[1:m] .* pi
u0[J .+ N] .= Z[1+m:2m] .* freqamp
prob.u0 .= u0
prob
end
# check convergence to sync
function _bstab_out(sol, i)
u = sol.u[end];
Δu = maximum( abs.( u[1+N:end] .- u_sync[1+N:end] ) )
S = (Δu < tolerance)
return (S, false)
end
return prob = EnsembleProblem(ode_prob;
output_func = _bstab_out,
prob_func = _bstab_init
)
end