Plotting parameters and variables on the same plot for Differential Equation solution

General Usage
, ,
1

Hi Everyone,
I am using DifferentialEquations.jl to solve an ODE. After solving it, I plot the different variables.
The problem is:
I want to add a parameter (p) to the same plot. Is it possible? I only see different ways of plotting the variables, but no way to add parameters to the plot.
The reason I want to see the parameter, is because it is changed a few times using iterative callbacks (DiffEqCallbacks).

2

Could you explain what you mean by “see the parameter”? How would it appear on the plot - is it a point on the plot, a line, some text somewhere?

3

Do you want to display the parameter as a text ( annotation) on your plot? Are you plotting a solution coordinates, t\to x_i(t), for each i, or the trajectory t\to (x_1(t),\ldots ,x_n(t)), n=2 or 3?
A minimal example and a more precise formulation of your intention to display the parameter for some interval of time, t, would be useful.

4

My plot is two variable plotted over time (which is the solution of the differential equation). While the parameter is a constant value given as an input to the differential equation.
The parameter would be a line on the plot. It would have been a straight horizontal line, but since the callback function changes the value at specified times, it should look like steps when plotted over time.

Here is my code, if it might help:

using Pkg, ModelingToolkit, BlockSystems, OrdinaryDiffEq, Plots, Interpolations, Symbolics

function interp(t) 
    time =   [0.05, 1.0,   2.0, 3.0, 6.0, 7.0, 8.0]
    P_elec = [1.0, 1.0,   0.6, 1.0, 1.0, 1.6,  1.6]
    interpoleer = linear_interpolation(time, P_elec, extrapolation_bc=Line())
    return interpoleer(t)
end

@parameters t M D P_m(t) P_e(t)
@variables ω(t)

dt = Differential(t)
swing = IOBlock([dt(ω) ~ 1/M * (P_m - D*ω - P_e)], # equations
                 [P_m],                       # inputs
                 [ω])                              # outputs

swing = set_p(swing, :D=>1, :M=>1)

@variables P_m(t)
pfix = IOBlock([P_m ~ 1], [], [P_m])
swing_with_P_m_P_e = IOSystem([pfix.P_m => swing.P_m], [pfix, swing]) |> connect_system

gen = generate_io_function(swing_with_P_m_P_e)

u0 = zeros(length(gen.states))
odef = ODEFunction((du, u, p, t) -> gen.f_ip(du, u, nothing, p, t); mass_matrix=gen.massm, syms=Symbol.(gen.states))

tspan = (0.0, 10.0)
p = [0.99] #P_e as parameter
prob = ODEProblem(odef, u0, tspan, p)

time_choice(integrator) = 0.05 + integrator.t
affect!(integrator) = integrator.p .= [interp(integrator.t)] 
cb = IterativeCallback(time_choice, affect!)
sol = solve(prob, Rodas4(autodiff=false), callback = cb) 

size = (500, 300)
linewidth=2
palette=:tab10
grid=false
xlabel = "Time"

plot(sol; linewidth, legend=:right, palette, grid, xlabel, ylabel="Frequency / Power", size, title="fixed P_m")
5

@empet, I added my code to try better explain this in a reply to @nilshg . Hope it helps?

6

A horizontal line is:

hline!([parameter_value], label = "parameter")

but if you want P_elec at different points in time it’s:

plot!(time, P_elec, label = "parameter")

In general the differential equations ecosystem uses Plots.jl recipes so all the plots once created can just be used like any other Plots.jl plot, i.e. they can be modified by the relevant commands with a bang ! added.

7

Thanks.

Though it didn’t really solve the problem. I am trying to plot “p” over “t”, not “P_elec” over “time”.
The code: affect!(integrator) = integrator.p .= [interp(integrator.t)], calls values of P_elec at different times, but plotting P_elec over Time will not give me the same as “p” over “t”.

8

You’d have to save it yourself. p isn’t logged.

1 Like
9

Thanks @ChrisRackauckas.
I’ll go and try that. Do you have any references/links where I can learn more about logging parameters?

10

Just add a push in the callback

11

Thanks. It works!

For others that might want to see the push!():

function simulate(sys)
    gen = generate_io_function(sys)

    u0 = zeros(length(gen.states))
    odef = ODEFunction((du, u, p, t) -> gen.f_ip(du, u, nothing, p, t); mass_matrix=gen.massm, syms=Symbol.(gen.states))

    #@info "Gen state" gen.states
    tspan = (0, 5)
    p = [0.8] #P_e as parameter >>>Is this only P_e, or is this a combination of the D, M, and P_e parameters?
    #@info p
    prob = ODEProblem(odef, u0, tspan, p)

    
    # Logging callback
    logged_p = []
    logged_t = []
    time_choice(integrator) = 0.001 + integrator.t
    affect!(integrator) = begin
        integrator.p .= [interp(integrator.t)]
        push!(logged_p, integrator.p[1])
        push!(logged_t, integrator.t)
    end
    cb = IterativeCallback(time_choice, affect!)

    sol = solve(prob, Rodas4(autodiff=false), callback=cb)
    return sol, logged_p, logged_t
end