Rounded values linked from subproblem violate constraints in next subproblem

[nmgeek]

#395 pointed out how solvers may have different rounding strategies and this can break tests. This issue is different in that a) only the CBC solver is used and b) the rounding leads to infeasible termination during the model run.

Specifically, I am running a week-by-week optimization (53 subproblems) and the max_footroom_energy_rule constraint is violated by approximately .00001 MWh. This constraint is:

mod.Stor_Downward_Reserves_MW[s, tmp] \
        * mod.hrs_in_tmp[tmp] \
        * mod.stor_charging_efficiency[s] \
        <= mod.Energy_Capacity_MWh[s, mod.period[tmp]] \
        * mod.Availability_Derate[s, tmp] - \
        mod.Stor_Starting_Energy_in_Storage_MWh[s, tmp] \
        - mod.Stor_Charge_MW[s, tmp] \
        * mod.hrs_in_tmp[tmp] \
        * mod.stor_charging_efficiency[s] \
        + mod.Stor_Discharge_MW[s, tmp] \
        * mod.hrs_in_tmp[tmp] \
        / mod.stor_discharging_efficiency[s]

Reserves and derate are set to defaults and the violation occurs during charging so it simplifies to:

0 <= Energy_Capacity_MWh - Stor_Starting_Energy_in_Storage_MWh - Stor_Charge_MW * stor_charging_efficiency

Stor_Charge_MW, as stored in the stor_linked_timepoints.tab, is rounded up in the 4th digit so the expression on the right evaluates to slightly less than zero, leading to infeasibility.

Note that in a healthy power system storage will be near full-capacity at midnight, when optimizations with weekly or daily subproblems would roll over to the next subproblem. So it is actually quite common for the capacity to be reaching the full mark in the last timepoint of a subproblem.

As a test, I added floor rounding in the code that saves stor_linked_timepoints.tab in stor.py. That change worked around this bug. The delta is below.

I feel like this has to be an issue with the output of CBC rather than rounding in python. Since python has infinite precision there should be no rounding in python.

My fix below is a suggestion for how this might be resolved. Perhaps there is a feature in CBC which would allow you to specify how rounding is done? I have not dug into that. Or, would it work to make the charge and discharge variables integers?

If we use this fix there should also be a "mirror fix" for discharge power.

             for (p, tmp) in sorted(mod.STOR_OPR_TMPS):
                 if tmp in tmps_to_link:
+                    # Ensure that rounded values do not violate max energy storage constraints.
+                    starting = max(value(mod.Stor_Starting_Energy_in_Storage_MWh[p, tmp]), 0)
+                    charge = max(value(mod.Stor_Charge_MW[p, tmp]), 0)
+                    print('TTT rounded', charge, starting)
+                    charge = math.floor(100000 * min(charge,
+                                                     (value(mod.Energy_Capacity_MWh[p, mod.period[tmp]])
+                                                      - starting) / value(mod.stor_charging_efficiency[p]))) / 100000.
                     writer.writerow([
                         p,
                         tmp_linked_tmp_dict[tmp],
-                        max(value(mod.Stor_Starting_Energy_in_Storage_MWh[
-                                      p, tmp]),
-                            0
-                            ),
+                        starting,
                         max(value(mod.Stor_Discharge_MW[p, tmp]), 0),
-                        max(value(mod.Stor_Charge_MW[p, tmp]), 0)
+                        charge,
                     ])

[anamileva]

I'd have to think about a generic fix as I wouldn't want to force rounding to a particular number of sig figs either.

[nmgeek]

Personally, I would not worry if it rounded to the nearest Watt. I can't think of any application of GridPath that would need finer resolution. Regardless, any kind of rounding in Gridpath seems messy: it would be nice to find a more elegant solution.

Another solution would be to "clamp" the Charge_MW (and Discharging_MW) parameter value when loading the data from stor_linked_timepoints.tab. You would add code like shown below to stor.load_model_data() after the parameter values are loaded from the file. I prefer this approach because it protects GP from out-of-bounds input data.

Let me know if you would like a PR.

(In this case "clamping" is distinct from rounding.)

mod.stor_linked_charge.value = min(value(mod.stor_linked_charge), (value(mod.Energy_Capacity_MWh) - value(mod.stor_linked_starting_energy_in_storage)) / value(mod.stor_charging_efficiency))