Skip to content

feat: value forecast and departure risk in the schedule #182

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
andig wants to merge 1 commit into from
Closed

feat: value forecast and departure risk in the schedule #182

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
14 changes: 14 additions & 0 deletions client/client.gen.go
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

25 changes: 25 additions & 0 deletions openapi.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
20 changes: 17 additions & 3 deletions src/optimizer/app.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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', {
Expand All @@ -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', {
Expand Down Expand Up @@ -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
Expand All @@ -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,
Expand All @@ -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:
Expand Down
109 changes: 68 additions & 41 deletions src/optimizer/optimizer.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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:
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down
Loading
Loading