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

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).

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?

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.

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")



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

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.

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â€ť.

Youâ€™d have to save it yourself. p isnâ€™t logged.

1 Like

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

Just add a push in the callback

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