Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Activity
janlisse
@janlisse
if i solve the 3 ODE's separatly it works fine, but then i have to assemble the solutions later on
that's why i wanted to use the Integrator interface in the first place. to save myself from combining 3 solutions. Runtime performance maybe another concern but isn't critical
BridgingBot
@GitterIRCbot
[slack] <SebastianM-C> I think this could also be modeled with a DiscreteCallback
you would need one (or more) parameter(s) that modify the equations when the failure occurs, and then modify that via the callback
(This would give you one solution automatically)
[slack] <chrisrackauckas> I still don't have that JSON though. Is it in the repo?
[slack] <chrisrackauckas> the path doesn't reference the package dir.
@SebastianM-C We also had a look at this. This could be another approach but is not easy for us since we have a lot of code that generates the systems of ODE's.
BridgingBot
@GitterIRCbot

[slack] <chrisrackauckas> ```
using Pkg
pkg"add PowerDynamics#PowerDrop"
using PowerDynamics: read_powergrid, Json, find_operationpoint, rhs, simulate, PowerDrop
using OrdinaryDiffEq
powergrid = read_powergrid(Pkg.dir("PowerDynamics")*"/ieee-14-minimal.json", Json)

operationpoint = find_operationpoint(powergrid)

function PowerDynamics.simulate(pd::PowerDrop, powergrid, x0, timespan)
@assert first(timespan) <= pd.tspan_fault[1] "fault cannot begin in the past"
@assert pd.tspan_fault[2] <= last(timespan) "fult cannot end in the future"

problem = ODEProblem{true}(rhs(powergrid), x0.vec, timespan)
integrator = init(problem, Rodas4(autodiff=false))

step!(integrator, pd.tspan_fault[1], true)

# update integrator with error
integrator.f = rhs(pd(powergrid))

step!(integrator, pd.tspan_fault[2], true)

# update integrator, clear error
integrator.f = rhs(powergrid)

solve!(integrator)

return PowerGridSolution(integrator.sol, powergrid)

end

result = simulate(PowerDrop(
fraction = 0.9,
node_number = 1,
tspan_fault = (2.,3.)),
powergrid, operationpoint, (0., 5.))
```

is the MWE. In the future it would be nice if the reproducible example was shared.

[slack] <chrisrackauckas> I'll dig in. If it's not a timing thing then :shrug:
janlisse
@janlisse
I wasn't aware that you could specify a branch in the ]add. good to know :)
and thanks @ChrisRackauckas !
BridgingBot
@GitterIRCbot
[slack] <SebastianM-C> Couldn't you use something like this for the problem function
if(!has_failed)
  rhs(powergrid)
else
  rhs(pd(powergrid))
end
janlisse
@janlisse
@SebastianM-C looks quite straightforward, i will give that a try!
BridgingBot
@GitterIRCbot
[slack] <Pooya Rashidi Ravari> It doesn't work with SRIW1()
[slack] <chrisrackauckas> yes, I said there is no bridging distribution
[slack] <chrisrackauckas> and no Z0
[slack] <chrisrackauckas> do w = TruncatedWienerProcess(0.0,zeros(5),zeros(5))
BridgingBot
@GitterIRCbot
[slack] <Pooya Rashidi Ravari> thanks
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> @janlisse I tried it and it seems like the error estimates are just high after you change the f, so it seems like an f issue. I assume there's a bunch of cache stored in the f?
[slack] <chrisrackauckas> yeah @janlisse, the issue is in your f there, something to do with what you're caching

[slack] <chrisrackauckas> ```
function PowerDynamics.simulate(pd::PowerDrop, powergrid, x0, timespan)
@assert first(timespan) <= pd.tspan_fault[1] "fault cannot begin in the past"
@assert pd.tspan_fault[2] <= last(timespan) "fult cannot end in the future"

problem = ODEProblem{true}(deepcopy(rhs(pd(powergrid))), x0.vec, timespan)
integrator = init(problem, Rodas4(autodiff=false))

@show pd.tspan_fault
OrdinaryDiffEq.step!(integrator, pd.tspan_fault[1], true)

# update integrator with error
integrator.f = rhs(pd(powergrid))
u_modified!(integrator,true)

OrdinaryDiffEq.step!(integrator, pd.tspan_fault[2], true)

# update integrator, clear error
integrator.f = rhs(powergrid)

solve!(integrator)

return PowerDynamics.PowerGridSolution(integrator.sol, powergrid)

end
```

[slack] <chrisrackauckas> that works, with the only change being the deepcopy
[slack] <chrisrackauckas> so basically that means that some values from the first part are entering the second part that you didn't intend?
janlisse
@janlisse
@ChrisRackauckas If i use your variant i get a solution indeed. But it starts the simulation with the fault, and the version that we need needs to start with the regular state.
Means i need:
deepcopy(rhs(powergrid))
instead of
deepcopy(rhs(pd(powergrid))
If i reverse that, than even with deepcopy() i get back to the infinite looping.
AFAIK there is no caching involved in the f.
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> is the model well-defined in that order?
janlisse
@janlisse
I'm not sure if i fully understand your question. When i solve 3 separate ODE's in that order it works out. Thus it should be the same for the Integrator version. Or is there a difference i don't see yet?
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> I'd have to dig into what you're doing in order to know what's up with your function. There are no arrays, nothing enclosed?
janlisse
@janlisse
@ChrisRackauckas There are arrays used internally. The powergrid model we are using basically consists of edges representing the lines and vertices representing electrical nodes. We then use a dedicated layer where we transform this graph structure into an ODEFunction that we could use for DiffEq. The most relevant code snippets for this transformation are:
function network_dynamics(vertices!::Array{ODEVertex,1}, edges!::Array{StaticEdge,1}, graph)
    @assert length(vertices!) == length(vertices(graph))
    @assert length(edges!) == length(edges(graph))

    dim_v = [v.dim for v in vertices!]
    dim_e = [e.dim for e in edges!]
    dim_nd = sum(dim_v)

    e_int = zeros(sum(dim_e))

    graph_stucture = create_graph_structure(graph, dim_v, dim_e, e_int)

    nd! = nd_ODE_Static(vertices!, edges!, graph, graph_stucture)

    # Construct mass matrix
    mass_matrix = construct_mass_matrix([v.mass_matrix for v in vertices!], dim_nd, graph_stucture)

    symbols = [Symbol(vertices![i].sym[j],"_",i) for i in 1:length(vertices!) for j in 1:dim_v[i]]

    ODEFunction(nd!,mass_matrix = mass_matrix,syms=symbols)
end
And:
function (d::nd_ODE_Static)(dx, x, p, t)
    gs = d.graph_stucture
    @views begin
    for i in 1:gs.num_e
        d.edges![i].f!(gs.e_int[gs.e_idx[i]], x[gs.s_idx[i]], x[gs.d_idx[i]], p, t)
    end
    for i in 1:gs.num_v
        d.vertices![i].f!(dx[gs.v_idx[i]], x[gs.v_idx[i]], gs.e_s[i], gs.e_d[i], p, t)
    end
    end # views
    nothing
end
Since for this case only the vertices contain dynamic behaviour and the lines/edges are static and not used in the ODE.
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> There's an issue somehow with how those get shared.
BridgingBot
@GitterIRCbot
[slack] <ashton.bradley> has anyone noticed any speed regressions today? I see a major speed regression once solve has completed. After that, even simple function calls are taking seconds, and the REPL feels cludgy, like it is out of memory. Julia claims to be only using about 3GB of 16. There was an atom update this morning, just testing on another machine to see if that could be it …
BridgingBot
@GitterIRCbot
[slack] <isaacsas> Did you just update to Julia 1.2?
[slack] <ashton.bradley> no, still on 1.1. Should I be updating yet?
[slack] <isaacsas> No, just noticed some people complaining about increased plotting times on 1.2 on discourse so was curious.
[slack] <ashton.bradley> It doesn’t seem to be anything to do with the Atom update
BridgingBot
@GitterIRCbot
[slack] <ashton.bradley> ok, looks like it could be an issue with my code...
BridgingBot
@GitterIRCbot
[slack] <tbeason> @chrisrackauckas in the DiffEqFinancial package you talk about "solving" those PDEs, but it seems when you say solve you mean solve for the index price as a function of time. Is that correct?
[slack] <tbeason> If so, do you have any examples where you solve the actual pricing equation, such as the very first equation presented here https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_equation
BridgingBot
@GitterIRCbot
[slack] <tbeason> I'm interested in contributing, but I am definitely a work-from-a-template kind of guy on these types of things
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> We just don't have the PDEs in there yet :shrug:
[slack] <chrisrackauckas> we need a higher level interface for defining PDEProblems
[slack] <chrisrackauckas> for now it would still be useful to directly define the ODE described by the PDE though.
[slack] <chrisrackauckas> The library is squarely about the SDEs right now though
Nicholas J. Matzke
@nmatzke
Hi all! I have been working with Julia's ODE tools and I find them amazing! I have a quick, possibly dumb question, however: Let's say my ODE has a rate parameter that should vary in time as a deterministic function of an external, empirical time series (e.g. a temperature curve from weather data). I think this would require something like a subfunction that (a) knows the absolute value of time-after-start-point, and (b) looks up the relevant temperature on the (smoothed) temperature curve. Is this kind of thing possible? Cheers and thanks!
BridgingBot
@GitterIRCbot

[slack] <yingbo_ma> Would

function foo(u, p, t)
    ... p(t)
end

work?

BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> @nmatzke mix @yingbo_ma’s answer with https://github.com/PumasAI/DataInterpolations.jl