From 9afd1e74e94a36a3ece210fbbc640523460f4e2b Mon Sep 17 00:00:00 2001 From: andig Date: Sun, 4 Oct 2026 16:42:21 +0200 Subject: [PATCH] feat: value forecast and departure risk in the schedule (closes #163) p_departure per battery weighs p_a by the probability the device is still there, ft_err per step prices the grid at the mean over ft +- ft_err with the schedule shared and only the grid split per scenario. Both optional, absent they reproduce the previous model. --- README.md | 4 +- client/client.gen.go | 14 +++ openapi.yaml | 25 ++++ src/optimizer/app.py | 20 +++- src/optimizer/optimizer.py | 109 ++++++++++------- .../029-departure-risk-charges-early.json | 109 +++++++++++++++++ .../030-forecast-band-charges-early.json | 109 +++++++++++++++++ tests/test_bad_request.py | 2 +- tests/test_continuity_api.py | 4 +- tests/test_risk.py | 85 +++++++++++++ tools/bench_risk.py | 113 ++++++++++++++++++ 11 files changed, 546 insertions(+), 48 deletions(-) create mode 100644 test_cases/029-departure-risk-charges-early.json create mode 100644 test_cases/030-forecast-band-charges-early.json create mode 100644 tests/test_risk.py create mode 100644 tools/bench_risk.py diff --git a/README.md b/README.md index 8e3230a31..1ce2d18c6 100644 --- a/README.md +++ b/README.md @@ -92,13 +92,15 @@ A schedule the solver stopped on at the limit is reported as `Feasible` rather t "c_min": 1380, "c_max": 11000, "d_max": 0, // W "s_goal": [0, 0, 0, 0, 0, 0, 0, 40000, 0, 0, 0, 0], // Wh per step "charge_from_grid": true, - "p_a": 0.00022 // value of stored energy, per Wh + "p_a": 0.00022, // value of stored energy, per Wh + "p_departure": [0, 0, 0, 0, 0, 0, 0, 0.1, 0.3, 0.4, 0, 0] // chance of leaving, per step } ], "time_series": { "dt": [3600, 3600, 3600, 3600, 3600, 3600, 3600, 3600, 3600, 3600, 3600, 3600], "gt": [230, 210, 205, 205, 240, 380, 780, 920, 640, 480, 430, 460], // demand, Wh "ft": [0, 0, 0, 0, 0, 60, 320, 850, 1600, 2500, 3300, 3900], // PV forecast, Wh + "ft_err": [0, 0, 0, 0, 0, 30, 160, 430, 800, 1250, 1650, 1950], // its standard deviation, Wh "p_N": [0.00026, 0.00024, 0.00023, 0.00023, 0.00025, 0.00030, 0.00036, 0.00041, 0.00038, 0.00032, 0.00027, 0.00022], // import, per Wh "p_E": [0.00008, 0.00008, 0.00008, 0.00008, 0.00008, 0.00008, diff --git a/client/client.gen.go b/client/client.gen.go index bce40aa71..abeeb67bf 100644 --- a/client/client.gen.go +++ b/client/client.gen.go @@ -83,6 +83,13 @@ type BatteryConfig struct { // PDemand Minimum charge demand per time step (Wh) PDemand []float32 `json:"p_demand,omitempty"` + // PDeparture Probability that the device leaves during each time step, summing to at most 1. The remainder + // is the probability it stays past the horizon. p_a then values the state of charge the device + // is expected to leave with instead of the state at the end of the horizon, so a charge that may + // be cut short by the departure is worth less than one that is already in. Without it the device + // is assumed to stay. + PDeparture []float32 `json:"p_departure,omitempty"` + // SCapacity The capacity at 100% SOC in Wh. If not specified s_capacity will be set to s_max. // s_initial must be less or equal s_capacity, otherwise the optimization will return an error. SCapacity float32 `json:"s_capacity,omitempty"` @@ -250,6 +257,13 @@ type TimeSeries struct { // Ft Forecasted energy generation (e.g., solar PV) at each time step (Wh) Ft []float32 `json:"ft"` + // FtErr Standard deviation of the generation forecast at each time step (Wh). Where given, the grid + // cost of the step is the mean over the generation coming in at ft - ft_err and at ft + ft_err + // instead of the cost at ft. The schedule is the same in both, the deviation lands on the grid: + // relying on forecasted generation that may fall short is charged the import it then needs. + // Steps where export pays more than import are evaluated at ft. + FtErr []float32 `json:"ft_err,omitempty"` + // Gt Household energy demand at each time step (Wh) Gt []float32 `json:"gt"` diff --git a/openapi.yaml b/openapi.yaml index 61d12b131..522238a59 100644 --- a/openapi.yaml +++ b/openapi.yaml @@ -280,6 +280,19 @@ components: maximum: 2 default: 0 description: Charging and discharging priority 0..2 compared to other batteries. 2 = highest priority. + p_departure: + type: array + items: + type: number + minimum: 0 + maximum: 1 + description: | + Probability that the device leaves during each time step, summing to at most 1. The remainder + is the probability it stays past the horizon. p_a then values the state of charge the device + is expected to leave with instead of the state at the end of the horizon, so a charge that may + be cut short by the departure is worth less than one that is already in. Without it the device + is assumed to stay. + example: [0, 0, 0.1, 0.3, 0.4, 0] TimeSeries: type: object @@ -325,6 +338,18 @@ components: minimum: 0 description: Grid export remuneration per Wh at each time step (currency units/Wh) example: [0.15, 0.12, 0.10, 0.11, 0.14, 0.16] + ft_err: + type: array + items: + type: number + minimum: 0 + description: | + Standard deviation of the generation forecast at each time step (Wh). Where given, the grid + cost of the step is the mean over the generation coming in at ft - ft_err and at ft + ft_err + instead of the cost at ft. The schedule is the same in both, the deviation lands on the grid: + relying on forecasted generation that may fall short is charged the import it then needs. + Steps where export pays more than import are evaluated at ft. + example: [0, 600, 1600, 2100, 1600, 500] OptimizationInput: type: object diff --git a/src/optimizer/app.py b/src/optimizer/app.py index ab1603045..1e8c0d639 100644 --- a/src/optimizer/app.py +++ b/src/optimizer/app.py @@ -144,7 +144,9 @@ def handle_validation_error(error): 'c_active': fields.Boolean(required=False, description='Whether the device is charging at the start of the time horizon.'), 'd_max': fields.Float(required=True, description='Maximum discharge power (W)'), 'p_a': fields.Float(required=True, description='Monetary value per Wh at end of the optimization horizon'), - 'c_priority': fields.Integer(required=False, description='Charging and discharging priority compared to other batteries. 2 = highest priority.') + 'c_priority': fields.Integer(required=False, description='Charging and discharging priority compared to other batteries. 2 = highest priority.'), + 'p_departure': fields.List(fields.Float, required=False, + description='Probability that the device leaves during each time step. The remainder stays past the horizon.'), }) time_series_model = api.model('TimeSeries', { @@ -153,6 +155,8 @@ def handle_validation_error(error): 'ft': fields.List(fields.Float, required=True, description='Forecasted solar generation at each time step (Wh)'), 'p_N': fields.List(fields.Float, required=True, description='Price per Wh taken from grid at each time step'), 'p_E': fields.List(fields.Float, required=True, description='Remuneration per Wh fed into grid at each time step'), + 'ft_err': fields.List(fields.Float, required=False, + description='Standard deviation of the solar forecast at each time step (Wh). Grid cost is taken as the mean over ft +- ft_err.'), }) optimization_input_model = api.model('OptimizationInput', { @@ -236,6 +240,7 @@ def post(self): d_max=bat_data['d_max'], p_a=bat_data['p_a'], c_priority=bat_data.get('c_priority', 0), + p_departure=bat_data.get('p_departure'), )) # Parse time series data @@ -245,6 +250,7 @@ def post(self): ft=data['time_series']['ft'], p_N=data['time_series']['p_N'], p_E=data['time_series']['p_E'], + ft_err=data['time_series'].get('ft_err'), ) # Validate time series lengths. dt included: the model indexes every series by it, @@ -254,12 +260,20 @@ def post(self): 'p_N': len(time_series.p_N), 'p_E': len(time_series.p_E), 'p_demand': [len(bat.p_demand) for bat in batteries if bat.p_demand is not None], 's_goal': [len(bat.s_goal) for bat in batteries if bat.s_goal is not None], + 'p_departure': [len(bat.p_departure) for bat in batteries if bat.p_departure is not None], } + if time_series.ft_err is not None: + lengths['ft_err'] = len(time_series.ft_err) - if len({*[v for k, v in lengths.items() if k not in ('p_demand', 's_goal')], - *lengths['p_demand'], *lengths['s_goal']}) > 1: + per_battery = ('p_demand', 's_goal', 'p_departure') + if len({*[v for k, v in lengths.items() if k not in per_battery], + *[n for k in per_battery for n in lengths[k]]}) > 1: api.abort(400, "All time series must have the same length", lengths=lengths) + for bat in batteries: + if bat.p_departure is not None and (min(bat.p_departure) < 0 or sum(bat.p_departure) > 1 + 1e-6): + api.abort(400, "p_departure must be probabilities summing to at most 1") + except BadRequest: raise except Exception as e: diff --git a/src/optimizer/optimizer.py b/src/optimizer/optimizer.py index 7bb16a8da..07e3d2ab1 100644 --- a/src/optimizer/optimizer.py +++ b/src/optimizer/optimizer.py @@ -207,6 +207,7 @@ class BatteryConfig: s_goal: Optional[List[float]] = None # Goal state of charge (Wh) c_priority: int = 0 c_active: bool = False # Whether the device is charging at the start of the horizon + p_departure: Optional[List[float]] = None # Probability the device leaves during each step @dataclass @@ -216,6 +217,7 @@ class TimeSeriesData: ft: List[float] # Forecasted production [Wh] p_N: List[float] # Import prices [currency unit/Wh] p_E: List[float] # Export prices [currency unit/Wh] + ft_err: Optional[List[float]] = None # Standard deviation of the production forecast [Wh] class Optimizer: @@ -314,6 +316,46 @@ def __init__(self, strategy: OptimizationStrategy, grid: GridConfig, batteries: # grid sides leveled by the active peak attenuation strategy, empty for all other strategies self.peak_sides = PEAK_STRATEGY_SIDES.get(strategy.charging_strategy, ()) + # steps whose grid cost is taken as the mean over the forecast band ft +- ft_err instead of + # at ft. Skipped where export pays more than import: the band's import and export are not + # locked against each other by a direction binary, so there the pair could pocket the + # difference on paper. + self.band_steps = set() + if time_series.ft_err is not None: + self.band_steps = {t for t in self.time_steps + if time_series.ft_err[t] > 0 and time_series.p_N[t] >= time_series.p_E[t]} + + def _soc_value_weights(self, bat: BatteryConfig) -> np.ndarray: + """ + Probability that s[t] is the state the battery is left in, which p_a then values. Without + departure probabilities that is the last step alone. + """ + weights = np.zeros(self.T) + if bat.p_departure is not None: + weights[:] = bat.p_departure + if self.T: + weights[-1] += 1 - weights.sum() + return weights + + def _grid_value(self, value=lambda var: var): + """ + Export revenue less import cost over the horizon. value maps a variable to what enters the + sum: the variable itself for the objective, pulp.value for the result. Import beyond + p_max_imp is priced like the rest, the limit itself is enforced by its penalty. + """ + total = 0 + for t in self.time_steps: + p_N, p_E = self.time_series.p_N[t], self.time_series.p_E[t] + if t in self.band_steps: + total += sum(0.5 * (value(self.variables[f'e_{k}'][t]) * p_E - value(self.variables[f'n_{k}'][t]) * p_N) + for k in ('lo', 'hi')) + continue + e_imp = value(self.variables['n'][t]) + if self.grid.p_max_imp is not None: + e_imp += value(self.variables['e_imp_lim_exc'][t]) + total += value(self.variables['e'][t]) * p_E - e_imp * p_N + return total + def create_model(self): """ Create and initialize the MILP model @@ -387,6 +429,11 @@ def _setup_variables(self): self.variables['n'] = [pulp.LpVariable(f"n_{t}", lowBound=0) for t in self.time_steps] self.variables['e'] = [pulp.LpVariable(f"e_{t}", lowBound=0) for t in self.time_steps] + # grid import/export with the production one ft_err below (lo) and above (hi) the forecast [Wh] + for k in ('lo', 'hi'): + self.variables[f'n_{k}'] = {t: pulp.LpVariable(f"n_{k}_{t}", lowBound=0) for t in self.band_steps} + self.variables[f'e_{k}'] = {t: pulp.LpVariable(f"e_{k}_{t}", lowBound=0) for t in self.band_steps} + # penalty variables for exceeding grid power limits (W) # for grid import if self.grid.p_max_imp is not None: @@ -449,30 +496,14 @@ def _setup_target_function(self): ############################################################################ # actual cost & benefit elements - # Grid import cost (negative because we want to minimize cost) [currency unit] - for t in self.time_steps: - # if a demand rate beyond p_max_imp is applied, both portions have to be considered - # for energy cost. If only an import limit is given, there should never be power - # import beyond p_max_imp, however, if the limit gets violated, we account for - # the energy cost as well to stay consistent. - if self.grid.p_max_imp is not None: - objective -= ( - # grid import up to the demand rate threshold - self.variables['n'][t] - # import beyond the threshold - + self.variables['e_imp_lim_exc'][t] - ) * self.time_series.p_N[t] - else: - # standard case - objective -= self.variables['n'][t] * self.time_series.p_N[t] - - # Grid export revenue [currency unit] - for t in self.time_steps: - objective += self.variables['e'][t] * self.time_series.p_E[t] + # Grid export revenue less import cost [currency unit] + objective += self._grid_value() - # Final state of charge value [currency unit] + # Value of the state of charge the battery is left in [currency unit] for i, bat in enumerate(self.batteries): - objective += self.variables['s'][i][-1] * bat.p_a + for t, weight in enumerate(self._soc_value_weights(bat)): + if weight: + objective += self.variables['s'][i][t] * weight * bat.p_a # charge for import power demand rate. The demand rate is applied to the maximum # power draw beyond the threshold within the time horizon. @@ -603,6 +634,16 @@ def _add_energy_balance_constraints(self): == e_grid_exp + self.time_series.gt[t]) + # the same balance with the production off by ft_err either way. The schedule stays, + # the deviation lands on the grid: a shortfall cannot raise the export, a surplus + # cannot raise the import. + if t in self.band_steps: + err = self.time_series.ft_err[t] + self.problem += self.variables['n_lo'][t] - self.variables['e_lo'][t] == e_grid_imp - e_grid_exp + err + self.problem += self.variables['n_hi'][t] - self.variables['e_hi'][t] == e_grid_imp - e_grid_exp - err + self.problem += self.variables['e_lo'][t] <= e_grid_exp + self.problem += self.variables['n_hi'][t] <= e_grid_imp + # Constraints (4)-(5): Grid flow direction. Export/import have natural # per-step caps (export: solar plus discharge-to-grid capacity; import: # demand plus grid-charge capacity) that tighten the LP relaxation vs @@ -1277,30 +1318,16 @@ def get_clean_objective_value(self): ''' recalculate the objective value without penalties and strategy icentives ''' - clean_objective = 0 - # Grid import cost (negative because we want to minimize cost) [currency unit] - for t in self.time_steps: - if self.grid.p_max_imp is not None: - clean_objective -= ( - # grid import up to the demand rate threshold - pulp.value(self.variables['n'][t]) - # import beyond the threshold - + pulp.value(self.variables['e_imp_lim_exc'][t]) - ) * self.time_series.p_N[t] - else: - # standard case - clean_objective -= pulp.value(self.variables['n'][t]) * self.time_series.p_N[t] - - # Grid export revenue [currency unit] - for t in self.time_steps: - clean_objective += pulp.value(self.variables['e'][t]) * self.time_series.p_E[t] + # Grid export revenue less import cost [currency unit] + clean_objective = self._grid_value(pulp.value) # Value of the energy the horizon added to the batteries [currency unit]. The reference is # bat.s_initial and not s[0]: s[0] is the state after the first time step and already # carries that step's charging, so subtracting it dropped the first step from the result. for i, bat in enumerate(self.batteries): - clean_objective += (pulp.value(self.variables['s'][i][self.T-1]) - - bat.s_initial) * bat.p_a + for t, weight in enumerate(self._soc_value_weights(bat)): + if weight: + clean_objective += (pulp.value(self.variables['s'][i][t]) - bat.s_initial) * weight * bat.p_a # charge for import power demand rate. The demand rate is applied to the maximum # power draw beyond the threshold within the time horizon. diff --git a/test_cases/029-departure-risk-charges-early.json b/test_cases/029-departure-risk-charges-early.json new file mode 100644 index 000000000..23c0d8130 --- /dev/null +++ b/test_cases/029-departure-risk-charges-early.json @@ -0,0 +1,109 @@ +{ + "strict": true, + "request": { + "batteries": [ + { + "charge_from_grid": true, + "discharge_to_grid": false, + "s_capacity": 5000, + "s_min": 0, + "s_max": 5000, + "s_initial": 0, + "c_min": 0, + "c_max": 5000, + "d_max": 0, + "p_a": 0.00035, + "p_departure": [ + 0.3, + 0, + 0, + 0 + ] + } + ], + "time_series": { + "dt": [ + 3600, + 3600, + 3600, + 3600 + ], + "gt": [ + 0, + 0, + 0, + 0 + ], + "ft": [ + 0, + 0, + 0, + 0 + ], + "p_N": [ + 0.0003, + 0.000305, + 0.00029, + 0.00028 + ], + "p_E": [ + 0.0001, + 0.0001, + 0.0001, + 0.0001 + ] + }, + "eta_c": 0.95, + "eta_d": 0.95 + }, + "expected_response": { + "status": "Optimal", + "objective_value": 0.1625000000000003, + "limit_violations": { + "grid_import_limit_exceeded": false, + "grid_export_limit_hit": false + }, + "batteries": [ + { + "charging_power": [ + 5000.0, + 0.0, + 0.0, + 0.0 + ], + "discharging_power": [ + 0.0, + 0.0, + 0.0, + 0.0 + ], + "state_of_charge": [ + 4750.0, + 4750.0, + 4750.0, + 4750.0 + ] + } + ], + "grid_import": [ + 5000.0, + 0.0, + 0.0, + 0.0 + ], + "grid_export": [ + 0.0, + 0.0, + 0.0, + 0.0 + ], + "flow_direction": [ + 0, + 0, + 0, + 0 + ], + "grid_import_overshoot": [], + "grid_export_overshoot": [] + } +} diff --git a/test_cases/030-forecast-band-charges-early.json b/test_cases/030-forecast-band-charges-early.json new file mode 100644 index 000000000..be3243b01 --- /dev/null +++ b/test_cases/030-forecast-band-charges-early.json @@ -0,0 +1,109 @@ +{ + "strict": true, + "request": { + "batteries": [ + { + "charge_from_grid": false, + "discharge_to_grid": false, + "s_capacity": 2000, + "s_min": 0, + "s_max": 2000, + "s_initial": 0, + "c_min": 0, + "c_max": 2000, + "d_max": 0, + "p_a": 0.0002 + } + ], + "time_series": { + "dt": [ + 3600, + 3600, + 3600, + 3600 + ], + "gt": [ + 0, + 0, + 0, + 0 + ], + "ft": [ + 0, + 2000, + 0, + 2000 + ], + "ft_err": [ + 0, + 200, + 0, + 1000 + ], + "p_N": [ + 0.0003, + 0.0003, + 0.0003, + 0.0003 + ], + "p_E": [ + 0.0001, + 0.00011, + 0.0001, + 0.0001 + ] + }, + "eta_c": 0.95, + "eta_d": 0.95 + }, + "expected_response": { + "status": "Optimal", + "objective_value": 0.5984210502, + "limit_violations": { + "grid_import_limit_exceeded": false, + "grid_export_limit_hit": false + }, + "batteries": [ + { + "charging_power": [ + 0.0, + 1105.2632, + 0.0, + 1000.0 + ], + "discharging_power": [ + 0.0, + 0.0, + 0.0, + 0.0 + ], + "state_of_charge": [ + 0.0, + 1050.0, + 1050.0, + 2000.0 + ] + } + ], + "grid_import": [ + 0.0, + 0.0, + 0.0, + 0.0 + ], + "grid_export": [ + 0.0, + 894.73684, + 0.0, + 1000.0 + ], + "flow_direction": [ + 1, + 1, + 1, + 1 + ], + "grid_import_overshoot": [], + "grid_export_overshoot": [] + } +} diff --git a/tests/test_bad_request.py b/tests/test_bad_request.py index e9efd8b55..fb877dd43 100644 --- a/tests/test_bad_request.py +++ b/tests/test_bad_request.py @@ -86,7 +86,7 @@ def test_length_mismatch_names_the_series(payload, series, capsys): response = app.test_client().post('/optimize/charge-schedule', json=payload) assert response.status_code == 400 - lengths = {'dt': 2, 'gt': 2, 'ft': 2, 'p_N': 2, 'p_E': 2, 'p_demand': [], 's_goal': [2]} + lengths = {'dt': 2, 'gt': 2, 'ft': 2, 'p_N': 2, 'p_E': 2, 'p_demand': [], 's_goal': [2], 'p_departure': []} lengths[series] = 1 assert response.json == {'message': 'All time series must have the same length', 'lengths': lengths} assert json.loads(capsys.readouterr().out)['bad_request'] == { diff --git a/tests/test_continuity_api.py b/tests/test_continuity_api.py index 30e65afaa..591bb8f90 100644 --- a/tests/test_continuity_api.py +++ b/tests/test_continuity_api.py @@ -15,7 +15,7 @@ def test_api_returns_continuous_equal_price_sessions(second_c_min: float | None) model.batteries.append(replace(model.batteries[0], c_min=second_c_min)) request = { 'batteries': [{key: value for key, value in asdict(battery).items() if value is not None} for battery in model.batteries], - 'time_series': asdict(model.time_series), + 'time_series': {key: value for key, value in asdict(model.time_series).items() if value is not None}, 'eta_c': 1, 'eta_d': 1, } @@ -38,7 +38,7 @@ def test_api_keeps_a_running_session_charging(): model.batteries[0].c_active = True request = { 'batteries': [{key: value for key, value in asdict(model.batteries[0]).items() if value is not None}], - 'time_series': asdict(model.time_series), + 'time_series': {key: value for key, value in asdict(model.time_series).items() if value is not None}, 'eta_c': 1, 'eta_d': 1, } diff --git a/tests/test_risk.py b/tests/test_risk.py new file mode 100644 index 000000000..e4cecdccb --- /dev/null +++ b/tests/test_risk.py @@ -0,0 +1,85 @@ +import copy +import json +import pathlib + +import numpy +import pytest + +from optimizer.app import app + +# the two stored risk cases, each against the same request without the risk field. The stored +# case pins the early schedule, these pin that it is the field that moves it there. +DEPARTURE = json.loads(pathlib.Path('test_cases/029-departure-risk-charges-early.json').read_text())['request'] +BAND = json.loads(pathlib.Path('test_cases/030-forecast-band-charges-early.json').read_text())['request'] + + +def solve(request): + response = app.test_client().post('/optimize/charge-schedule', json=request) + assert response.status_code == 200, response.json + assert response.json['status'] == 'Optimal' + return response.json + + +def charging(result): + return numpy.array(result['batteries'][0]['charging_power']) + + +def test_departure_probability_moves_the_charge_ahead_of_the_departure(): + without = copy.deepcopy(DEPARTURE) + del without['batteries'][0]['p_departure'] + + # the cheapest hour is the last one, so without the risk the charge waits for it + assert charging(solve(without)).argmax() == 3 + # a 30 % chance of leaving in the first hour outweighs the 0.1 the wait would save + assert charging(solve(DEPARTURE)).argmax() == 0 + + +def test_departure_value_is_the_expected_state_of_charge(): + result = solve(DEPARTURE) + bat = DEPARTURE['batteries'][0] + soc = numpy.array(result['batteries'][0]['state_of_charge']) + weights = numpy.array(bat['p_departure'], float) + weights[-1] += 1 - weights.sum() + expected = ((soc - bat['s_initial']) * weights).sum() * bat['p_a'] \ + - (numpy.array(result['grid_import']) * numpy.array(DEPARTURE['time_series']['p_N'])).sum() + assert numpy.isclose(result['objective_value'], expected) + + +def test_forecast_band_moves_the_charge_to_the_surer_surplus(): + without = copy.deepcopy(BAND) + del without['time_series']['ft_err'] + + # the first surplus pays a little more exported, so at the forecast the battery waits for the second + assert charging(solve(without)).argmax() == 3 + # the second surplus may fall 1000 Wh short, which the battery would then draw from the grid + assert charging(solve(BAND)).argmax() == 1 + + +def test_forecast_band_cost_is_the_mean_over_the_band(): + result = solve(BAND) + series = BAND['time_series'] + net = numpy.array(result['grid_import']) - numpy.array(result['grid_export']) + value = 0 + for err in (-1, 1): + shifted = net - err * numpy.array(series['ft_err']) + value += 0.5 * (numpy.minimum(shifted, 0) * -numpy.array(series['p_E']) + - numpy.maximum(shifted, 0) * numpy.array(series['p_N'])).sum() + value += result['batteries'][0]['state_of_charge'][-1] * BAND['batteries'][0]['p_a'] + assert numpy.isclose(result['objective_value'], value) + + +@pytest.mark.parametrize('p_departure', [[0.5, 0.6, 0, 0], [-0.1, 0.5, 0, 0]]) +def test_departure_probabilities_are_validated(p_departure): + request = copy.deepcopy(DEPARTURE) + request['batteries'][0]['p_departure'] = p_departure + response = app.test_client().post('/optimize/charge-schedule', json=request) + assert response.status_code == 400 + assert 'p_departure' in response.json['message'] + + +def test_short_ft_err_is_a_length_mismatch(): + request = copy.deepcopy(BAND) + request['time_series']['ft_err'] = [0, 200] + response = app.test_client().post('/optimize/charge-schedule', json=request) + assert response.status_code == 400 + assert response.json['lengths']['ft_err'] == 2 diff --git a/tools/bench_risk.py b/tools/bench_risk.py new file mode 100644 index 000000000..90e92602f --- /dev/null +++ b/tools/bench_risk.py @@ -0,0 +1,113 @@ +"""Benchmark the forecast band and the departure probabilities against a reference revision. + +Runs every stored test case four ways: through the reference revision, through the working tree as +stored, with a forecast band of 20 percent of ft growing with the square root of the lead time, and +with a 30 percent departure chance spread over hours 6 to 12 on every battery that cannot discharge. +Reports solve time, model size and whether the schedule moved. + + uv run python tools/bench_risk.py [--rev main] [--repeat 3] +""" +import argparse +import copy +import json +import pathlib +import sys +import time + +import numpy as np + +ROOT = pathlib.Path(__file__).resolve().parents[1] +sys.path.insert(0, str(ROOT / 'src')) +sys.path.insert(0, str(ROOT / 'tools')) + +from bench_peak_leveling import build, load_reference # noqa: E402 + +import optimizer.optimizer as new # noqa: E402 the working tree under test, needs the path above + + +def with_band(req): + req = copy.deepcopy(req) + hours = np.cumsum(req['time_series']['dt']) / 3600 + ft = np.array(req['time_series']['ft'], float) + req['time_series']['ft_err'] = (0.2 * np.sqrt(hours / 24) * ft).tolist() + return req + + +def with_departure(req): + req = copy.deepcopy(req) + hours = np.cumsum(req['time_series']['dt']) / 3600 + window = (hours > 6) & (hours <= 12) + if not window.any(): + return None + p = np.where(window, 0.3 / window.sum(), 0.0) + for bat in req['batteries']: + if bat['d_max'] == 0: + bat['p_departure'] = p.tolist() + return req + + +def without_risk(req): + """the request as the reference revision understands it""" + req = copy.deepcopy(req) + req['time_series'].pop('ft_err', None) + for bat in req['batteries']: + bat.pop('p_departure', None) + return req + + +def build_new(req): + opt = build(new, req) + for bat, cfg in zip(req['batteries'], opt.batteries): + cfg.p_departure = bat.get('p_departure') + opt.time_series.ft_err = req['time_series'].get('ft_err') + return new.Optimizer(opt.strategy, opt.grid, opt.batteries, opt.time_series, opt.eta_c, opt.eta_d, opt.M) + + +def run(make, req, repeat): + best, res, opt = None, None, None + for _ in range(repeat): + o = make(req) + started = time.perf_counter() + r = o.solve() + best = min(best or 1e9, time.perf_counter() - started) + res, opt = r, o + return {'t': best, 'res': res, 'rows': len(opt.problem.constraints), 'cols': len(opt.problem.variables())} + + +def charging(res): + return np.array([b['charging_power'] for b in res['batteries']]) + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument('--rev', default='main', help='revision to compare against') + ap.add_argument('--repeat', type=int, default=3, help='solves per case, fastest counts') + args = ap.parse_args() + ref = load_reference(args.rev) + + print(f'{"case":46s} {"t_ref":>6s} {"t_new":>6s} {"t_band":>6s} {"t_dep":>6s} {"rows":>12s} {"cols":>10s} band departure') + totals = np.zeros(4) + for path in sorted((ROOT / 'test_cases').glob('*.json')): + req = json.loads(path.read_text())['request'] + a = run(lambda r: build(ref, without_risk(r)), req, args.repeat) + b = run(build_new, req, args.repeat) + band = run(build_new, with_band(req), args.repeat) + dep_req = with_departure(req) + dep = run(build_new, dep_req, args.repeat) if dep_req else None + times = [a['t'], b['t'], band['t'], dep['t'] if dep else 0] + totals += times + # a stored case with a risk field is a different model for the reference by design + same = (np.allclose(charging(a['res']), charging(b['res']), atol=1) + or np.isclose(a['res']['objective_value'], b['res']['objective_value'], atol=1e-4) + or without_risk(req) != req) + band_moved = not np.allclose(charging(b['res']), charging(band['res']), atol=1) + dep_moved = dep and not np.allclose(charging(b['res']), charging(dep['res']), atol=1) + print(f'{path.stem:46s} {a["t"]:6.3f} {b["t"]:6.3f} {band["t"]:6.3f} {dep["t"] if dep else 0:6.3f} ' + f'{a["rows"]:5d}->{band["rows"]:5d} {a["cols"]:4d}->{band["cols"]:4d} ' + f'{"moved" if band_moved else "same ":5s} {"moved" if dep_moved else "same" if dep else "-"}' + f'{"" if same else " REF DIFFERS"}') + print(f'{"total":46s} {totals[0]:6.3f} {totals[1]:6.3f} {totals[2]:6.3f} {totals[3]:6.3f}') + + +if __name__ == '__main__': + main()