Optimizing a step function

[cstjean]

Hi! I'm evaluating InfiniteOpt to model our chemical process. I've got an ODE problem where the control signal is a step function (i.e. changing the operating conditions requires physical intervention and is done once a week), and I'd like to figure out the optimal value of that control signal at each steps (i.e. MPC). That sounds like something that should be easy to do in InfiniteOpt, but I'm new to JuMP so it's not immediately obvious. Should one use @finite_parameter for such a thing?

Another question: if I'm reading the MathOptAI documentation right, one can now use machine learning models inside JuMP/InfiniteOpt, but the training still has to be done outside the model, right? It's not possible to do the Scientific ML thing of fitting an ML model embedded inside a set of ODEs, minimizing the RMSE between the prediction of the ODE and the data points?

[pulsipher]

Hi @cstjean,

To better help you with the step functions, can you provide some more details? What are the time horizon and step size you have in mind? How many times will you be able to update the control variable? I suspect, we can model this using domain restrictions, but I can provide more specific instructions with some more info (a small toy example formulation would be even better).

Yes, we can embed ML models in InfiniteOpt via MathOptAI. And we can train them using https://lanl-ansi.github.io/MathOptAI.jl/stable/api/#MathOptAI.replace_weights_with_variables. Just keep in mind that training these directly can be computationally expensive. We are currently developing some efficient ways to handle this though.

[cstjean]

We have an electrochemical reactor. Current density cannot be changed very easily. The control signal (current density) Startup + normal operation + shutdown looks like this:

Except that the "normal operation" part should be 100X longer (startup takes 12hours, normal operation takes months).

The main question to answer for us is: what is the optimal stable operation current (the final current density).

Other modeled quantities (temperature, concentration etc.) are a function of the current density, and rapidly reach a steady point.

From my limited understanding of collocation methods, we could have a discretization that puts many points around step changes, then space out those points during stable operation. If that could be achieved with InfiniteOpt, that'd be a good step forward for us.

[pulsipher]

Thanks for the info, states with significantly different time scales typically presents an interesting challenge. Strategically placing collocation points would certainly help. Some ways to handle the step inputs would include:

• Defining finite elements that align with times to adjust the current density, adding collocation points to more accurately integrate the states, and enforcing that the control signal be kept constant over each finite element via https://infiniteopt.github.io/InfiniteOpt.jl/stable/manual/variable/#InfiniteOpt.constant_over_collocation-Tuple{InfiniteVariableRef,%20GeneralVariableRef}.
• Using domain restrictions to keep the inputs constants over predefined ranges
• Segment the time horizon into coupled time stages

[cstjean]

I tried domain restrictions today, and everything went smoothly, I can build a step function. Looking forward to writing our model POC in April!

Terminology question: does the kwarg supports: Sets the support points for the parameters. refer to collocation points? Or "nodes"? I still don't have the hang of the numerical details.

[pulsipher]

Supports refer to the set of discretization points applied to an infinite parameter. For instance, let's define t with 3 support points:

@infinite_parameter(model, t in [0, 1], supports = [0, 0.5, 1])

This defines $t\in[0,1]$ that will discretized using $\{0, 0.5, 1\}$. Note that by default these supports will be used with backward finite difference to approximate derivatives. To use orthogonal collocation over finite elements, we have to specify the numerical method we would like to apply to the derivatives via derivative_method:

@infinite_parameter(model, t in [0, 1], supports = [0, 0.5, 1], derivative_method = OrthogonalCollocation(3))

This will apply collocation where $\{0,0.5,1\}$ constitute the finite elements where each use 3 Lobatto collocation points (i.e., number of nodes). Hence, the total set of discretization points becomes $\{0, 0.25, 0.5, 0.75, 1\}$. To learn more, see https://infiniteopt.github.io/InfiniteOpt.jl/stable/guide/derivative/#Derivative-Methods and https://infiniteopt.github.io/InfiniteOpt.jl/stable/manual/derivative/#InfiniteOpt.OrthogonalCollocation.

[cstjean]

UPDATE: I studied more and realized my misunderstandings. I'll post some update later today.

I've been making good progress. I have nearly all of my model basics translated over, and things usually work (in 100 Ipopt iterations or so), but for some seemingly-random values of the inputs, Ipopt fails to find a solution (hits maxiter=3000).

EDIT: to make it clearer, currently the optimized input signal on successful runs looks like this; we want to find the optimal value of each plateau:

I suspect that I'm doing something wrong. Here's how the input signal I is done:

""" Inside half-open-interval. Left-side is open EXCEPT at 0, where it is closed, to cover the space.
"""
ino(interval) = in(Interval{minimum(interval) == 0 ? :closed : :open, :closed}(interval))

function model()
    t_step_change = 0:20:200   # the boundaries of the plateaus
    ϵ = 0.1
    supportΔs = vcat(-ϵ, 0, +ϵ, 0.2:0.4:5)   # add some support after the plateau boundaries, to properly model the transient
    supports = filter(x->0<=x<=200,
                      sort(reduce(vcat, [step .+ supportΔs for step in t_step_change])))    

    @infinite_parameter(m, t in [0, max_t]; supports=supports)
                        #derivative_method=OrthogonalCollocation(3))

    @variable m 0 <= I_plateau[1:N_electro, 1:N_steps] <= 14  # having the <=14 constraint on I
                                                              # instead works out better for some
                                                              # problems, worse for others
    @variable m I[1:N_electro] Infinite(t)
    #constant_over_collocation.(I, t)

    # Inputs
    for E in 1:N_electro
        for i in 1:N_steps
            @constraint(m, I[E] == I_plateau[E, i],
                        DomainRestriction(ino(t_step_change[i]..t_step_change[i+1]), t))
        end
    end
    ...
end

OrthogonalCollocation(3) works, but usually takes longer than the default (or hits maxiter). When I use both OrthogonalCollocation and the commented out constant_over_collocation call, I get TOO_FEW_DEGREES_OF_FREEDOM, and I don't understand why.

Another issue apparently comes from the timescale mismatch: I'd like to have maybe 5 hours of dense points after each plateau start, and then months of stable operation covered by only a few points (t is in hours). But apparently, the longer my plateaus are, the bigger I need to make ϵ. Naively, this was unexpected.

Sorry that I couldn't share the full example. Maybe I can produce an MWE later this week if needed.

Segment the time horizon into coupled time stages

Not sure what you're referring to here.

InfiniteOpt has been a pleasure to use so far 👍

[cstjean]

I had more time to study today and review the documentation.

Starting from your suggestions:

• A) Defining finite elements that align with times to adjust the current density, adding collocation points to more accurately integrate the states, and enforcing that the control signal be kept constant over each finite element via https://infiniteopt.github.io/InfiniteOpt.jl/stable/manual/variable/#InfiniteOpt.constant_over_collocation-Tuple{InfiniteVariableRef,%20GeneralVariableRef}.
• B) Using domain restrictions to keep the inputs constants over predefined ranges
• C) Segment the time horizon into coupled time stages

My previous question:

When I use both OrthogonalCollocation and the commented out constant_over_collocation call, I get TOO_FEW_DEGREES_OF_FREEDOM, and I don't understand why.

I can now answer: I was trying to implement both your A) and B) suggestions at the same time, but that's not a good idea. A and B are conflicting, hence the TOO_FEW_DEGREES_OF_FREEDOM outcome.

Domain restriction approach

I tried to make the domain restriction approach work better today (see the code from my previous post), as it already works about 50% of the time for small problems. One discovery is that while the final solution's input signal correctly looks like this:

If I interrupt the optimization early on, it looks like this:

I had naively hoped that the equality constraints over the first plateau's $I_1 = I_2 = I_3 = ... = I_{numsupportperplateau} = Iplateau_1$

    for E in 1:N_electro
        for i in 1:N_steps
            @constraint(m, I[E] == I_plateau[E, i],
                        DomainRestriction(ino(t_step_change[i]..t_step_change[i+1]), t))
        end
    end

would be recognized as such by the solver, and that it would have somehow eliminated the redundant variables as a preprocessing step, optimizing just one variable per plateau instead of N variables. Unfortunately, that is not the case, so the problem becomes exponentially more difficult for the solver. Is there any way that I can fix that? Should I use a different solver?

constant_over_collocation approach

I have made it work on some small problems, but in general it did not perform well for me. And I have two issues with it:

• If I understand well, each element should be as big as a single plateau. If I do that, then is there any way to "place the nodes" closer to the start of the plateau? I guess I could split each plateau into two elements, where the first one is smaller to accurately model the transient
• How can I tell what's going on inside the elements? With the domain restriction approach, I get some nice plausible curves for the integrated variables , but if I have only one or two elements per plateau, I get only one or two points. Related to Automatically choosing the "right" interpolation type #419

Incidentally, these plots' plateaus are very short for us. The process never gets into steady state, which is precisely the problem I'm trying to solve. The longer I make the plateaus, the more the solver struggles.

[pulsipher]

Looking at the Time Stage Approach, it feels rather similar to the domain restrictions approach. Wouldn't the transcription be very similar?

It is similar, but the time stage approach is a little more fool proof to implement.

In MTK, simple equations like cost_of_electricity = 10 * electricity_usage are basically free, because the solver will treat them as observables and substitute them in where needed. In InfiniteOpt, it looks like I should avoid introducing such extra variables?

You can accomplish this in InfiniteOpt and JuMP via @expression.

In our case, num_supports=low_value would be a good easy problem.

You can set the start values based on a previous solve via set_start_values. Though, I suspect using an integrator with fixed input would be more robust.

If this transcribes to a 10^6 difference in coefficients, then according to JuMP it's fundamentally problematic. It's a bummer. Would there be a way to automatically rescale the transcribed variables to make the coefficients reasonable?

Scaling can make a big difference. InfiniteOpt and JuMP do not have a way to automate this though, it needs to be done manually. I don't quite have enough context to describe how to do so in your case. However, having minute resolution over a year long horizon seems like an unnecessarily long time horizon. Roughly what is the time constant of the system? Have you considered a rolling horizon (MPC-type) approach?

[cstjean]

We made more good progress today by tweaking the start values and simplifying our non-linear operator. We're now able to solve for 4 years. Another trick was that I originally used (approximately) t+20minutes, t+40minutes, t+60minutes, t+6months for the support set (where t is the start of each plateau), but based on your advice, t+20minutes, t+40minutes, t+60minutes, t+1week, t+2week, ... t+26weeks should reduce the coefficient magnitude disparity, without increasing the support set size too much.

Thank you for all the help.

[pulsipher]

Ah I see, you are using a diffusing horizon, that makes sense.

[cstjean]

Yeah, I guess in our case it's a bit different. Basically we'll use a "fine grid" close to operating condition change, and a "coarse grid" further off (to model the steady state operation). I'd love to be able to use something like an exponentially increasing time horizon (like adaptive time-stepping does), but IIUC it wouldn't help with the constraint that the smallest timestep and the biggest timestep should be within 10^6 of each other, to avoid conditioning problems...?

In any case, given that 1month / 1minute = 4.e4, such a coarse / fine grid ratio seems fine for us both with regards to conditioning and inference speed.

[pulsipher]

Fair points. One thing I would point out is that typically the difference in scale is only critical between adjacent time steps (e.g., t_i and t_i+1) as are needed for the numerical methods that integrate the ODEs. Elsewhere, the time steps only act as indices for the variables unless time is explicitly part of an algebraic expression in a constraint. So, I would say that exponentially increasing the time step does not necessarily cause ill-conditioning. One could even use a discretization based on Gauss Laguerre quadrature which tailored for infinite time horizons.