Check out InfiniteOpt.jl: Home · InfiniteOpt.jl.
using InfiniteOpt
import Ipopt
import Plots
τ = 1.0
σ = 49.0
E0 = 100.0 # Note I made this smaller
fmax = 10.0
T = 20.0
model = InfiniteModel(Ipopt.Optimizer)
@infinite_parameter(model, t in [0, T], num_supports = 1_000)
@variables(model, begin
e >= 0, Infinite(t)
v, Infinite(t)
0 <= f <= fmax, Infinite(t)
end)
@objective(model, Max, integral(v, t))
@constraints(model, begin
deriv(v, t) == -v / τ + f
v(0) == 0
deriv(e, t) == σ - f * v
e(0) == E0
end)
optimize!(model)
@assert termination_status(model) == LOCALLY_SOLVED
ts, e_opt, v_opt, f_opt = supports(t), value(e), value(v), value(f)
Plots.plot(
Plots.plot(ts, e_opt, ylabel = "e(t)"),
Plots.plot(ts, v_opt, ylabel = "v(t)"),
Plots.plot(ts, f_opt, ylabel = "f(t)", xlabel = "t"),
layout = (3, 1),
legend = false,
)
cc @pulsipher