[FEATURE] Semidefinite integral

[cstjean]

Now that our core model is written, I was pleasantly surprised to see that I could add extra @expression's at no additional cost. Unfortunately, while I can compute this definite integral

@expression m total_cost integral(hourly_cost, t, 0, 100)

What I'd really like would be

@expression m cumulative_cost integral(hourly_cost, t, 0, t)

that computes the semi-definite integral

$$ \int_{0}^{t} hourlycost(\tau), d\tau $$

Is there any way to get that in InfiniteOpt? We have a lot of similar expressions, and adding them all to our model as variable + constraint makes me worry about runtime.

[pulsipher]

InfiniteOpt does not support using infinite parameters for integral bounds. Can you better clarify your use case? Do you want to compute these as a post-solve calculation or use them to pose your optimization problem? If you want to use them to pose your optimization problem, could you better clarify how?

If you are trying trying to define an accumulation quantity for the optimization problem, than the easiest thing to do would be to reformulate the integral into an ODE. For instance, let $c(t)$ be the accumulated cost at time $t$ and $ci(t)$ be the instantaneous cost at time $t$ such that we have

$$c(t) = \int_0^t ci(\tau) d\tau$$

Differentiating we obtain the equivalent ODE

$$\frac{d c(t)}{dt} = ci(t), \quad c(0) = 0$$

If you are trying to do a post-solve calculation, then the best thing to do would be to extract the optimal trajectories along with the support values and then numerically integrate the quantities you want using a dedicated package like QuadGK or NumericalIntegration, or even just use trapezoid rule directly.

[cstjean]

Thank you for your answer.

Do you want to compute these as a post-solve calculation or use them to pose your optimization problem?

Both. Let me explain our use case.

From my perspective, our mathematical problem can be divided into

• A) Core variables that must be jointly solved: say, relating pressure and temperature inside our reactor
• B) Peripheral variables that our user might "compute as a post-solve calculation, or use them to pose an optimization problem", like the cumulative electricity cost (an integral), total cost of ownership, etc.
• C) Additional constraints, like "today we'd like to do a shutdown between 3 and 5PM"
• D) Objective, like "minimize electricity cost", "minimize total cost of ownership", "minimize risk"

So far I've implemented A) and everything works well. Now comes B). Ideally I'd like to be able to put all the peripheral variables and equations as well into the model, and allow our users to pick their objectives, additional constraints etc. For instance, "Minimize electricity cost, then plot the total cost of ownership", or the reverse: "Minimize total cost of ownership, then plot the cumulative electricity cost". If all variables in B) are part of the model, and I do as you say:

the easiest thing to do would be to reformulate the integral into an ODE

then the code flows very naturally, as there is no difference between A) and B) quantities, and I can achieve the desired flexibility. Unfortunately, our model will be much bigger than just solving for A) quantities, and the solver will take longer the more peripheral variables we add. This is a worry.

@expression m hourly_electricity_cost I * power_price is great because

• It registers m[:hourly_electricity_cost]
• It incurs no additional solver cost if it is unused

So a user can put m[:hourly_electricity_cost] inside either the objective, or a plot call, and from my understanding it does not incur extra cost to the solver (beyond the objective evaluation); no new variable.

If semi-definite integrals were supported as @expression, then we could have the best of both worlds even for our quantities like cumulative_electricity_cost = integral(hourly_electricity_cost, t).

Perhaps this feature request is out-of-scope for InfiniteOpt. I believe that I can support this by creating my own wrapping model type, but of course having this logic natively supported would make it that much smoother.

[cstjean]

Basically, right now we have:

@expression m electricity_cost, power * hourly_rate
@constraint m electricity_cost(10) <= 5
...
plot(m[:time], m[:electricity_cost])

I'd like to have the same for cumulative_electricity_cost:

@expression m cumulative_electricity_cost, integral(electricity_cost, t, 0, t)
@constraint m cumulative_electricity_cost(10) <= 10
...
plot(m[:time], m[:cumulative_electricity_cost])

As you pointed out, all the tools are already there to do it myself. For the constraint, I can use another integral with a shorter integration bound. For the plot, I can use QuadGK. But it's very clunky compared to the elegance of the electricity_cost code.

Incidentally, given m[:cumulative_electricity_cost], is there a way to recover that it's the (definite) integral of electricity_cost? In JuMP, I can do

julia> @expression m ex v*sin(v)
v * sin(v)

julia> ex.args
2-element Vector{Any}:
 v
 sin(v)

to inspect the object, but the GeneralVariableRef that I get for cumulative_electricity_cost (when defined as a definite integral) is opaque.

[pulsipher]

It is important to point out that while expressions can be used for convenience in post analysis, this has not been the intent behind them, so the codebase is not optimized for this use case.

Now with regard to semi-definite integrals, I would recommend using if statements switch between use cases. If you want to use it with the solver, then pose it as an ODE (this will be better behaved numerically then a naive semi-definite integral anyways). Otherwise, if you want to post calculate with it then use QuadGK based on the optimal values. Since this is being used as a code base, I would recommend following the suggestions provided at https://jump.dev/JuMP.jl/stable/tutorials/getting_started/design_patterns_for_larger_models/.

In the longer term, we can add support for semi-definite integrals, but this will require some breaking changes to the integral API and as such will need to wait for the next major breaking release of InfiniteOpt.

To more intuitively inspect InfiniteOpt variables, use dispatch_variable_ref(vref::GeneralVariableRef) to get a more informative type. See https://infiniteopt.github.io/InfiniteOpt.jl/stable/guide/expression/#Variable-Hierarchy.

[cstjean]

Now with regard to semi-definite integrals, I would recommend using if statements switch between use cases. If you want to use it with the solver, then pose it as an ODE (this will be better behaved numerically then a naive semi-definite integral anyways)

Curious: are you suggesting that

@variable m integral_of_x Infinite(t)
@constraint m ∂(integral_of_x, t) == x
@constraint m intregal_of_x(0) == 0
@objective integral_of_x(10)

is recommended over

@objective m integral(x, t, 0, 10)

?

In the longer term, we can add support for semi-definite integrals, but this will require some breaking changes to the integral API and as such will need to wait for the next major breaking release of InfiniteOpt.

Thank you for the clear feedback. I believe that there are enough options in there to satisfy our needs in the short term.

[pulsipher]

Curious: are you suggesting that

@variable m integral_of_x Infinite(t)
@constraint m ∂(integral_of_x, t) == x
@constraint m intregal_of_x(0) == 0
@objective integral_of_x(10)
is recommended over

@objective m integral(x, t, 0, 10)
?

My observation is specific to semi-definite integrals that results in a very different transcribed algebraic structure than definite integrals. For definite integrals, I would recommend the use of integral directly. Though there could be certain use cases where an ODE is preferable. For instance, constrained integrals can result in KKT systems with a more challenging sparsity structure than using an equivalent ODE.