Add more powerful solver (non-linear capability)

[nhoening]

Next to Cbc, we want to ship support for a "stronger" solver in FlexMeasures, which nonlinear constraint support. The first candidate is ipopt.

With Cbc, there is a limit to the workarounds we can think of to get our storage scheduler to function well under various conditions. The main issue we're running into is that we cannot currently forbid our scheduler to make energy disappear. The only thing that is preventing that from happening is by introducing unfavourable conditions to do so. In other words, negative prices are a problem, and workarounds can reasonably be implemented only in relatively non-complex systems.

Because we plan to do optimization on larger (neighbourhood) systems, with possibly more than one energy contract, we'll really need to include a nonlinear constraint to solve this issue once and for all, for which cbc is no longer adequate.

We know by experience that getting a solver to run on any machine is using up developer time, so hopefully we can use ipopt within the FlexMeasures Docker image (like we could with Cbc). One alternative is to start using a web service like https://neos-server.org/neos/, but that would come with adding a major extra dependency (and also an external queueing system).

[ewinnington]

HiGHS might be an interesting solver to add support for, it is developing nicely for energy optimizations. https://github.com/ERGO-Code/HiGHS
Would be a simple and easy one to add to the docker image.

Also, could put CBC/HiGHS in a separate docker with an API on top to make the opt more scalable and then if you need to have a more powerful solver, an this API can interface to CPLEX or Gurobi. This could be put in place for companies that have a license for those solvers. Supporting Gurobi compute servers would allow to move the optimization to a separate server cluster, typical in an Energy company's setup.

As for Energy "disappearing" during an LP/MILP optimization, this can be dealt with by having the correct lp model with slack variables.

Happy to help discuss both these topics.

[Flix6x]

Thank you for pointing us to this, and thinking along. I'm just starting to read up on this solver, but I already found that work has been done on interfacing to the solver via Pyomo. I guess next steps could be to see:

• if the pyomo interface works
• if it can solve our LP problems
• if it can handle the constraint I have in mind to solve the issue of disappearing energy under negative prices

[ewinnington]

Do you have the equations for the https://github.com/FlexMeasures/flexmeasures/blob/main/flexmeasures/data/models/planning/linear_optimization.py documented somewhere as Latex or markdown+MathJax? If so I could take a look. Or you can just export the final equation set prior to solving a problem with pprint. I can take a look at what the issue is with your formulation. CBC is super capable and should not be causing any issues with negative prices.

I uncommented the lines

# model.pprint()
# model.display()

in the docker image of FlexMeasures, but was unable to find where the model would have been logged or printed (didn't seem to print to console the model). If I had found it I wouldn't bother you :)

[nhoening]

docker logs flexmeasures-worker-1 might be helpful here, as well (given the container name in he tutorial is likely to be "flexmeasures-worker-1").

[ewinnington]

The problem is that to stay in the linear world, you can't specify x*y = 0 as a constraint.

You have to either go with using a binary variable Z (and become a MILP problem) with:
$$x <= \infty Z ;\qquad \ y <= \infty (1-Z) ;\qquad \ Z \in \lbrace 0,1 \rbrace \in \mathbb{Z}$$

And as we know, Binary variables expensive and we would like to avoid them. But CBC can solve many binaries, so it is worth it to have the mode because it's easy to setup.

But If I understand the base problem you are defining, it is possible to restrict the state change of the battery to the max of the charge or discharge, so that we can get rid of some of the issues, namely the "double spending".

• $t$ is time point
• battery state at $b^t$ limited by $b_{max}$
• battery has a max charge and discharge rate that can be calculated once or per timepoint $r^t = \max({charge}, {discharge})$
• electricity price $e^t$
• charge efficiency $c_{eff}$ is in range {0..1}
• discharge efficiency $d_{eff}$ is 1/discharge - due to the way I use it to compute the output power
• Power input into battery $p^t$ with $p^t_{max}$ being the maximum power intake
• Power output from battery $o^t$ with $o^t_{max}$ being the maximum power outtake

If I was trying to maximize the value, the Objective would be: $\max \sum e^t * (p^t - o^t)$

Subject to:

• (0) Battery limits $0 <= b^t <= b_{max}$
• (1a) Power limits $p^t \leq p^t_{max}; \qquad o^t \leq o^t_{max}$
• (2) Battery balance equation $b^{t+1} = b^t + c_{eff} * p^t - d_{eff} * o^t$

To impose the additional restriction and staying in the linear world, we can decide to place limits the state charge of the battery by the sum of the

• (3a) charge/discharge rate limit $c_{eff} * p^t + d_{eff} * o^t \leq r^t$.

This would represent running the battery in a proportional time to the ratio of charge and discharge in a sub-timestep action and is actually a valid solution if you would allow the battery to switch between charge and discharge in a sub-timestep way.

I hope this helps. Sadly I'm not used to the way the pyomo sets up the equations, so I haven't understood your objective function yet, but the constraints explained here should be valid regardless of your objective function.

(edit - after posting I thought a little more) Instead of 3a, which causes problems if the charge rate and discharge rates are very different , we can formalize the fraction $f^t$ which actually represents the fraction of the timestep spent in the charging (respectively discharging) state. This replaces equations 1a and 3a to reduce to only 1b.

• fraction of timestep spent in potential charging state $f^t \in [0..1] \in \mathbb{R}$
• (1b) Power limits $p^t \leq f^t * p^t_{max}; \qquad o^t \leq (1-f^t) * o^t_{max}$

Again, still allows the battery to take sub timestep actions but formalises it and limits the in and out flow of power to the correct values, while only adding one linear variable. And if you want to make it a MILP problem, you simply change $f^t \in \mathbb{Z}$ as a binary => done!

[Flix6x]

Good stuff. I see no flaws in your formulation, but I do have practical issues.

1. The one extra binary variable per time slot amounts to quite a large combinatory problem for a scheduling horizon of even a few steps. I believe it scales with $2^n$ with $n$ the number of steps, which makes the problem intractable to solve in practice for scheduling horizons we typically encounter (which is from 24 up to 960 steps). (See my test results below.)
2. The fraction $f^t$ contains vital information for interpreting the resulting schedule, so it would need to be communicated together with the schedule, essentially making the schedule a time series of tuple values. Of course, this is all possible (we could store the fraction as sensor data, and communicate it together with the power values whenever a schedule is fetched), but, if we take a step back, this work would all just be for essentially creating a way to be wasteful with energy.

Test results

I ran a test with negative price values, with and without including a binary variable for each time slot, and recorded the time pytest took to complete (which has considerable overhead).

n steps	seconds with binary	seconds without binary
8	14	17
12	18	17
16	21	17
20	34	17
24	136	17
28	664	17
32	>180	18

[ewinnington]

So we need to distinguish two cases:

1. Is the output of the model correct or "better" with the new constraint?

1a) "It's not possible to both charge and discharge at some rate for the full duration of a single timeslot (these modes cannot be in effect simultaneously)" <- this should be fixed by the Binary constraint and mitigated by the linear constraint to represent only feasible operation plans where charging and discharging can only sum up to the timestep and not exceed it. Is this the case it is fixed?

1b) "Second, our scheduler returns the resulting schedule as the sum of the charging schedule and the discharging schedule. In the above example, that would equate to a consumption value of 19 kWh. The scheduler "knowns" that this energy represents an energy loss and doesn't end up increasing the battery's SoC, but any other code outside of the scheduler interprets this as energy that is scheduled to be consumed by the battery, resulting in 17.1 kWh being added to the SoC, and 1.9 kWh lost in the conversion. Downstream, when the schedule is used to compute how the SoC changes over time, this also leads to SoC's that surpass 100% of the storage capacity." This should be checked if the constraint change does resolve this problem and if not, what leads to this issue. Also, reading this again indicates to me that the the state of the battery is not captured from the output variable of the optimization, where it is available - but re-calculated afterwards.

2. Performance of a model can always be improved by multiple techniques once we know we have a correct model:

2A) In the first instance, we already have identified the specific case where the situation can happen: Negative prices. Therefore, you can generally make the model Linear and only create a binary variable if the price of the hour is negative. Since in LP, the prices have perfect foresight, this is a valid solution. Even in the worst case in the lowly connected grid of West Texas, only 1/3 of the hours of the Day Ahead Market were negative so it is a huge saving to only impose the condition when needed. It is easy to setup the variables to support that in code.

2B) MILP variables can be saved and reloaded from a previous run, or from a heuristic algorithm. So you can either solve the problem once and cache the results and pre-load them at the next run, or find some initial heuristic to pre-load binary variables that will be close to the solution vector. IE: If the prices are negative, you can heuristically set the binary variable to 1 as in we expect to charge, if the prices are positive we set the binary variable as a 0 (expect to discharge). The opt can still change these binary variables, they are just initial preloaded vars. CBC and other solvers support this.

2C) Depending on the length of the problem, you can cut the problem into smaller parts and define problem boundaries. This can then allow the solver to be run in parallel. This requires some more "thought" but is a common method used in large scales (european wide LP models - style of Energy Exemplar Plexos for example).

I can list many other techniques for both the modelling and performance, but we need to know which one to address.

3. if we take a step back, this work would all just be for essentially creating a way to be wasteful with energy. which is as specified in the problem currently a valid solution to our equations and in reality something that batteries can do to dissipate excess energy. Leads back to 1) change the model if that is not a desired behaviour - or identify a variable to penalize it.

[Flix6x]

Thank you, these are excellent suggestions. I need some time to test the impact of 2B and 2C. Hopefully, this will already bring the computation time down to a practical range. It sounds promising.

• 1A: I already verified that adding the binary constraint indeed solves the problem. 👍
• 1B: Indeed, we are not yet storing the SoC as computed by the scheduler. In fact, we aren't storing the SoC information at all yet. But that doesn't matter here, I think: in FlexMeasures, we made a modelling choice to store (and communicate) prosumption schedules as time series of power flows (positive values indicating consumption and negative values indicating production), rather than a tuple of time series with distinct consumption schedules and production schedules. This modelling choice does not work well when the scheduler can still decide to both consume and produce during the same time slot.

[nhoening]

@ewinnington here is an update of sorts:

We have now created a nice documentation of our scheduler, with formulas: https://flexmeasures.readthedocs.io/en/latest/concepts/device_scheduler.html

We also checked that the HiGHs solver is compatible, and it will be supported in the Docker image out of the box. Thanks for the hint! It seems that computation-time-wise, it might have a positive effect only on larger problems than we are dealing with at the moment (see a visual profile of where time is spent on the PR)

@Flix6x when you're back from vacation, it would be great to report back on 2B and 2C. We can also make issues out of them. For 2B, could we store MILP variables in Redis? Or in a blob on the sensor/asset we are scheduling? For 2C, I guess that is an option when we actually compute solutions for largers sites (with more flexible assets).

[Flix6x]

Looks like HiGHS is doing a much better job with adding a binary constraint for each time slot! Based on the test results below, I suggest we move to have FlexMeasures use HiGHS by default, and add the binary constraint for each time slot.

Test results with HiGHS

I reran my previous test with negative price values, this time including results for the HiGHS solver.

n steps	seconds with CBC and with binary	seconds with HiGHS and with binary	seconds with cbc and without binary	seconds with HiGHS and without binary
8	14	14	12	9
12	13	14	12	15
16	17	15	13	16
20	30	14	16	14
24	118	16	15	12
28	🐌	15	16	12
32	🐌	16	14	14
...	...	...	...	...
96	🐌	19	12	13