Skip to content

Commit 494df16

Browse files
BenLauwensBenLauwens
committed
direct qss solver implemented
1 parent cc44986 commit 494df16

7 files changed

Lines changed: 291 additions & 73 deletions

File tree

‎src/SimJulia.jl‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -28,8 +28,8 @@ module SimJulia
2828
export Coroutine, @coroutine
2929
export Container, Resource, Store
3030
export Put, Get, Request, Release, cancel, capacity, request, @request
31-
export Continuous, Variable, ZeroCrossing
32-
export @model, @continuous
31+
export Model, Continuous, Variable, ZeroCrossing
32+
export @model, @continuous, @trigger, @zerocrossing
3333
export QSS
3434
export non_stiff, stiff
3535

@@ -50,5 +50,6 @@ module SimJulia
5050
include("odes/base.jl")
5151
include("odes/macros.jl")
5252
include("continuous.jl")
53+
include("odes/roots.jl")
5354
include("odes/QSS.jl")
5455
end

‎src/continuous.jl‎

Lines changed: 32 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,15 @@
11
mutable struct Variable <: AbstractEvent
22
bev :: BaseEvent
33
id :: UInt
4-
f :: Function
54
x :: Taylor1
65
t :: Float64
7-
function Variable(env::Environment, id::Int)
8-
new(BaseEvent(env), UInt(id))
6+
function Variable(env::Environment, id::Int, x::Taylor1, t::Float64)
7+
new(BaseEvent(env), id, x, t)
98
end
109
end
1110

1211
function advance_time(var::Variable, t::Float64)
13-
x = var.x
14-
Δt = t - var.t
15-
var.x = evaluate(x, Δt + Taylor1(x.order))
12+
var.x = evaluate(var.x, t - var.t + Taylor1(var.x.order))
1613
var.t = t
1714
var.x.coeffs[1]
1815
end
@@ -21,23 +18,35 @@ struct ZeroCrossing <: AbstractEvent
2118

2219
end
2320

21+
macro zerocrossing(expr::Expr)
22+
expr.head != :call && error("Expression is not a function call!")
23+
nothing
24+
end
25+
2426
struct Continuous <: AbstractProcess
2527
bev :: BaseEvent
28+
integrator :: Integrator
2629
vars :: Vector{Variable}
27-
p :: Vector{Float64}
28-
zcs :: Vector{ZeroCrossing}
29-
function Continuous(env::Environment, p::Vector{Float64})
30-
new(BaseEvent(env), Vector{Variable}(), p, Vector{ZeroCrossing}())
30+
function Continuous(env::Environment, integrator::Integrator)
31+
cont = new(BaseEvent(env), integrator, Vector{Variable}())
32+
t = now(env)
33+
x = initial_values(integrator, t)
34+
for (i, x₀) in enumerate(x)
35+
push!(cont.vars, Variable(env, i, x₀, t))
36+
end
37+
for var in cont.vars
38+
schedule(var, cont, integrator)
39+
end
40+
cont
3141
end
3242
end
3343

34-
function Continuous{I<:Integrator}(model::Function, ::Type{I}, env::Environment, x₀::Vector{Float64}, p::Vector{Float64}=Float64[]; args...)
35-
cont = Continuous(env, p)
36-
for i in 1:length(x₀)
37-
push!(cont.vars, Variable(env, i))
44+
function Continuous{I<:Integrator}(model::Model, env::Environment, ::Type{I}, x₀::Vector{Float64}, p::Vector{Float64}=Float64[]; args...)
45+
for p₀ in p
46+
push!(model.p, p₀)
3847
end
39-
I(model, cont, now(env), x₀; args...)
40-
cont
48+
integrator = I(model, now(env), x₀; args...)
49+
Continuous(env, integrator)
4150
end
4251

4352
macro continuous(expr::Expr)
@@ -54,5 +63,11 @@ macro continuous(expr::Expr)
5463
push!(args, expr.args[i])
5564
end
5665
end
57-
esc(:(Continuous($(func), $(args...); $(params...))))
66+
esc(:(Continuous($func(), $(args...); $(params...))))
67+
end
68+
69+
function schedule(var::Variable, cont::Continuous, integrator::Integrator, Δt::Float64=0.0)
70+
var.bev = BaseEvent(environment(var))
71+
@callback step(var, cont, integrator)
72+
schedule(var, Δt)
5873
end

‎src/odes/QSS.jl‎

Lines changed: 79 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -2,70 +2,102 @@
22

33
struct QSS{T} <: Integrator
44
order :: UInt8
5-
deps :: Matrix{Bool}
5+
model :: Model
66
q :: Vector{Taylor1}
77
t :: Vector{Float64}
8-
function QSS{T}(model::Function, cont::Continuous, t::Float64, vec_x₀::Vector{Float64}; order::Number=4) where T
9-
f, deps = model()
10-
integrator = new(UInt8(order), deps, Vector{Taylor1}(), Vector{Float64}())
11-
zero_taylor = 0*Taylor1(Float64, integrator.order+1)
12-
for x₀ in vec_x₀
13-
push!(integrator.t, t)
14-
push!(integrator.q, x₀ + zero_taylor)
8+
Δrel :: Float64
9+
Δabs :: Float64
10+
function QSS{T}(model::Model, t::Float64, x₀::Vector{Float64};
11+
order::Number=4, Δrel::Float64=1e-6, Δabs::Float64=1e-6) where T
12+
qss = new(UInt8(order), model, Vector{Taylor1}(), Vector{Float64}(), Δrel, Δabs)
13+
t₀ = t + Taylor1(Float64, qss.order+1)
14+
for q₀ in x₀
15+
push!(qss.t, t)
16+
push!(qss.q, q₀ + Taylor1(zeros(Float64, qss.order+1)))
1517
end
16-
t₀ = t + Taylor1(Float64, integrator.order+1)
17-
vec_x = Vector{Taylor1}()
18-
for (i, x₀) in enumerate(vec_x₀)
19-
push!(vec_x, integrate(f[i](t₀, integrator.q, cont.p), x₀))
20-
end
21-
for i in 1:integrator.order-1
22-
for (j, x) in enumerate(vec_x)
23-
integrator.q[j] = copy(x)
24-
end
25-
for (j, x₀) in enumerate(vec_x₀)
26-
vec_x[j] = integrate(f[j](t₀, integrator.q, cont.p), x₀)
18+
for i in 1:qss.order-1
19+
for (j, q₀) in enumerate(x₀)
20+
qss.q[j] = integrate(model.f[j](t₀, qss.q, model.p), q₀)
2721
end
2822
end
29-
for (i, x) in enumerate(vec_x)
30-
var = cont.vars[i]
31-
var.f = f[i]
32-
var.x = x
33-
var.t = t
34-
@callback step(var, cont, integrator)
35-
schedule(var)
36-
end
37-
integrator
23+
qss
3824
end
3925
end
4026

41-
function Continuous(model::Function, env::Environment, x₀::Vector{Float64}, p::Vector{Float64}=Float64[]; order::Number=4)
42-
Continuous(model, QSS{non_stiff}, env, x₀, p; order=order)
27+
function Continuous(model::Model, env::Environment, x₀::Vector{Float64}, p::Vector{Float64}=Float64[];
28+
stiff::Bool=false, order::Number=4, Δrel::Float64=1e-6, Δabs::Float64=1e-6)
29+
stiff ? Continuous(model, env, QSS{SimJulia.stiff}, x₀, p; order=order, Δrel=Δrel, Δabs=Δabs) :
30+
Continuous(model, env, QSS{SimJulia.non_stiff}, x₀, p; order=order, Δrel=Δrel, Δabs=Δabs)
31+
end
32+
33+
function initial_values(qss::QSS, t::Float64)
34+
t₀ = t + Taylor1(Float64, qss.order+1)
35+
x₀ = Vector{Taylor1}()
36+
for (i, f) in enumerate(qss.model.f)
37+
push!(x₀, integrate(f(t₀, qss.q, qss.model.p), qss.q[i].coeffs[1]))
38+
end
39+
x₀
4340
end
4441

45-
function step(var::Variable, cont::Continuous, integrator::QSS)
42+
function step(var::Variable, cont::Continuous, qss::QSS)
43+
t = now(environment(var))
44+
n = length(qss.model.f)
4645
i = var.id
47-
env = environment(var)
48-
t = now(env)
46+
t₀ = t + Taylor1(Float64, qss.order+1)
4947
x₀ = advance_time(var, t)
50-
update_quantized_state(cont, integrator, i, t)
48+
println(t, " ", i, " ", x₀)
49+
update_quantized_state(qss, cont.vars, i, t)
50+
Δt = compute_next_time(var.x, max(qss.Δrel*x₀, qss.Δabs))
51+
schedule(var, cont, qss, Δt)
52+
for j in filter(j->qss.model.deps[j,i], 1:n)
53+
dep = cont.vars[j]
54+
x₀ = evaluate(dep.x, t - dep.t)
55+
dep.t = t
56+
for k in filter(k->qss.model.deps[j,k], 1:n)
57+
advance_time(qss, k, t)
58+
end
59+
dep.x = integrate(qss.model.f[j](t₀, qss.q, qss.model.p), x₀)
60+
Δt = recompute_next_time(qss, var.x, qss.q[j], max(qss.Δrel*x₀, qss.Δabs))
61+
schedule(dep, cont, qss, Δt)
62+
end
63+
end
64+
65+
function advance_time(qss::QSS, i::Int, t::Float64)
66+
qss.q[i] = evaluate(qss.q[i], t - qss.t[i] + Taylor1(qss.q[i].order))
67+
qss.t[i] = t
68+
qss.q[i].coeffs[1]
69+
end
70+
71+
function update_quantized_state(qss::QSS{non_stiff}, vars::Vector{Variable}, i::UInt, t::Float64)
72+
qss.q[i] = copy(vars[i].x)
73+
qss.q[i].coeffs[end] = 0.0
74+
qss.t[i] = t
75+
end
76+
77+
function update_quantized_state(qss::QSS{stiff}, vars::Vector{Variable}, i::UInt, t::Float64)
78+
for (j, istrue) in enumerate(qss.model.deps[i, :])
79+
istrue && advance_time(qss, j, t)
80+
end
81+
q₋ = deepcopy(qss.q)
5182
end
5283

53-
function advance_time(integrator::QSS, i::Int, t::Float64)
54-
q = integrator.q[i]
55-
Δt = t - integrator.t[i]
56-
integrator.q[i] = evaluate(q, Δt + Taylor1(q.order))
57-
integrator.t[i] = t
58-
integrator.q[i].coeffs[1]
84+
function compute_next_time(x::Taylor1, Δq::Float64)
85+
(abs(Δq/x.coeffs[end]))^(1.0/x.order)
5986
end
6087

61-
function update_quantized_state(cont::Continuous, integrator::QSS{non_stiff}, i::UInt, t::Float64)
62-
integrator.q[i] = Taylor1(cont.vars[i].x.coeffs[1:integrator.order])
63-
integrator.t[i] = t
88+
function recompute_next_time(::QSS{non_stiff}, x::Taylor1, q::Taylor1, Δq::Float64)
89+
p = x.coeffs - q.coeffs
90+
p[1] -= Δq
91+
neg = roots(p)
92+
p[1] += 2Δq
93+
pos = roots(p)
94+
select_root(neg, pos)
6495
end
6596

66-
function update_quantized_state(cont::Continuous, integrator::QSS{stiff}, i::UInt, t::Float64)
67-
for (j, istrue) in enumerate(integrator.deps[i, :])
68-
istrue && advance_time(integrator, j, t)
97+
function select_root(neg::Vector{Complex{Float64}}, pos::Vector{Complex{Float64}})
98+
selected = Inf
99+
for v in filter(v->abs(imag(v)) < 1e-15 && real(v) >= 0, [neg..., pos...])
100+
selected = real(v) < selected ? real(v) : selected
69101
end
70-
q₋ = deepcopy(integrator.q)
102+
selected
71103
end

‎src/odes/base.jl‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1 +1,10 @@
11
abstract type Integrator end
2+
3+
struct Model
4+
f :: Vector{Function}
5+
p :: Vector{Float64}
6+
deps :: Matrix{Bool}
7+
function Model(f::Vector{Function}, deps::Matrix{Bool})
8+
new(f, Vector{Float64}(), deps)
9+
end
10+
end

‎src/odes/macros.jl‎

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -24,8 +24,13 @@ macro model(expr::Expr)
2424
ex isa Expr && ex.head == Symbol("=") && check_dependencies(ex.args[2], deps, ex.args[1].args[2], expr.args[1].args[3])
2525
end
2626
esc(:(function $func_name()
27-
f = Array{Function}(length($f_vec))
28-
$((:(f[$(f_vec[i].args[1].args[2])] = (t::TaylorSeries.Taylor1, q::Vector{TaylorSeries.Taylor1}, p::Vector{Float64})->$(f_vec[i].args[2])) for i in 1:length(:($f_vec)))...)
29-
f, $deps
27+
f = Array{Function}($n)
28+
$((:(f[$(f_vec[i].args[1].args[2])] = ($(expr.args[1].args[2])::TaylorSeries.Taylor1, $(expr.args[1].args[3])::Vector{TaylorSeries.Taylor1}, $(expr.args[1].args[4])::Vector{Float64})->$(f_vec[i].args[2])) for i in 1:length(:($f_vec)))...)
29+
Model(f, $deps)
3030
end))
3131
end
32+
33+
macro trigger(expr::Expr)
34+
expr.head != :function && error("Expression is not a function definition!")
35+
nothing
36+
end

0 commit comments

Comments
 (0)