10385477682

Committed 14 Aug 2024 10:00AM UTC coverage: 67.393% (-0.03%) from 67.423%

Build # 10385477682

Build Type

Pull #1144

github

Committed by

victorgarcia98

Commit Message

Merge remote-tracking branch 'origin/feature/planning/multi-commitments' into feature/planning/multi-commitments

Signed-off-by: Victor Garcia Reolid <victor@seita.nl>

Pull Request Pull Request #1144: Feature/planning/multi-commitments

Run Details

6112 of 9891 branches covered (61.79%)

Branch coverage included in aggregate %.

59 of 92 new or added lines in 3 files covered. (64.13%)

48 existing lines in 4 files now uncovered.

7684 of 10580 relevant lines covered (72.63%)

5.75 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

85.09

/flexmeasures/data/models/planning/linear_optimization.py

from __future__ import annotations

from flask import current_app
import pandas as pd
import numpy as np
from pandas.tseries.frequencies import to_offset
from pyomo.core import (
    ConcreteModel,
    Var,
    RangeSet,
    Set,
    Param,
    Reals,
    NonNegativeReals,
    NonPositiveReals,
    Binary,
    Constraint,
    Objective,
    minimize,
)
from pyomo.environ import UnknownSolver  # noqa F401
from pyomo.environ import value
from pyomo.opt import SolverFactory, SolverResults

from flexmeasures.data.models.planning.utils import initialize_series, initialize_df
from flexmeasures.utils.calculations import apply_stock_changes_and_losses

infinity = float("inf")


def device_scheduler(  # noqa C901
    device_constraints: list[pd.DataFrame],
    ems_constraints: pd.DataFrame,
    commitment_quantities: list[pd.Series] | None = None,
    commitment_downwards_deviation_price: list[pd.Series] | list[float] | None = None,
    commitment_upwards_deviation_price: list[pd.Series] | list[float] | None = None,
    commitments: list[pd.DataFrame] | None = None,
    initial_stock: float = 0,
) -> tuple[list[pd.Series], float, SolverResults, ConcreteModel]:
    """This generic device scheduler is able to handle an EMS with multiple devices,
    with various types of constraints on the EMS level and on the device level,
    and with multiple market commitments on the EMS level.
    A typical example is a house with many devices.
    The commitments are assumed to be with regard to the flow of energy to the device (positive for consumption,
    negative for production). The solver minimises the costs of deviating from the commitments.

    :param device_constraints:   Device constraints are on a device level. Handled constraints (listed by column name):
        max: maximum stock assuming an initial stock of zero (e.g. in MWh or boxes)
        min: minimum stock assuming an initial stock of zero
        equal: exact amount of stock (we do this by clamping min and max)
        efficiency: amount of stock left at the next datetime (the rest is lost)
        derivative max: maximum flow (e.g. in MW or boxes/h)
        derivative min: minimum flow
        derivative equals: exact amount of flow (we do this by clamping derivative min and derivative max)
        derivative down efficiency: conversion efficiency of flow out of a device (flow out : stock decrease)
        derivative up efficiency: conversion efficiency of flow into a device (stock increase : flow in)
        stock delta: predefined stock delta to apply to the storage device. Positive values cause an increase and negative values a decrease
    :param ems_constraints:     EMS constraints are on an EMS level. Handled constraints (listed by column name):
        derivative max: maximum flow
        derivative min: minimum flow
    :param commitments:         Commitments are on an EMS level. Handled parameter (listed by column name):
        quantity:                   for example, 5.5
        downwards deviation price:  10.1
        upwards deviation price:    10.2
        group:                      1 (defaults to the enumerate time step j)

    Potentially deprecated arguments:
        commitment_quantities: amounts of flow specified in commitments (both previously ordered and newly requested)
            - e.g. in MW or boxes/h
        commitment_downwards_deviation_price: penalty for downwards deviations of the flow
            - e.g. in EUR/MW or EUR/(boxes/h)
            - either a single value (same value for each flow value) or a Series (different value for each flow value)
        commitment_upwards_deviation_price: penalty for upwards deviations of the flow

    All Series and DataFrames should have the same resolution.

    For now, we pass in the various constraints and prices as separate variables, from which we make a MultiIndex
    DataFrame. Later we could pass in a MultiIndex DataFrame directly.
    """

    model = ConcreteModel()

    # If the EMS has no devices, don't bother
    if len(device_constraints) == 0:
        return [], 0, SolverResults(), model

    # Get timing from first device
    start = device_constraints[0].index.to_pydatetime()[0]
    # Workaround for https://github.com/pandas-dev/pandas/issues/53643. Was: resolution = pd.to_timedelta(device_constraints[0].index.freq)
    resolution = pd.to_timedelta(device_constraints[0].index.freq).to_pytimedelta()
    end = device_constraints[0].index.to_pydatetime()[-1] + resolution

    # Move commitments from old structure to new
    if commitments is None:
        commitments = []
    if commitment_quantities is not None:
        for quantity, down, up in zip(
            commitment_quantities,
            commitment_downwards_deviation_price,
            commitment_upwards_deviation_price,
        ):

            # Turn prices per commitment into prices per commitment flow
            if all(isinstance(price, float) for price in down) or isinstance(
                down, float
            ):
                down = initialize_series(down, start, end, resolution)
            if all(isinstance(price, float) for price in up) or isinstance(up, float):
                up = initialize_series(up, start, end, resolution)

            group = initialize_series(list(range(len(down))), start, end, resolution)
            df = initialize_df(
                ["quantity", "downwards deviation price", "upwards deviation price"],
                start,
                end,
                resolution,
            )
            df["quantity"] = quantity
            df["downwards deviation price"] = down
            df["upwards deviation price"] = up
            df["group"] = group
            commitments.append(df)

    # Check if commitments have the same time window and resolution as the constraints
    for commitment in commitments:
        start_c = commitment.index.to_pydatetime()[0]
        resolution_c = pd.to_timedelta(commitment.index.freq)
        end_c = commitment.index.to_pydatetime()[-1] + resolution
        if not (start_c == start and end_c == end):
            raise Exception(
                "Not implemented for different time windows.\n(%s,%s)\n(%s,%s)"
                % (start, end, start_c, end_c)
            )
        if resolution_c != resolution:
            raise Exception(
                "Not implemented for different resolutions.\n%s\n%s"
                % (resolution, resolution_c)
            )

    def convert_commitments_to_subcommitments(
        dfs: list[pd.DataFrame],
    ) -> list[pd.DataFrame]:
        """Transform commitments, each specifying a group for each time step, to sub-commitments, one per group.

        We also enumerate the time steps in a new column "j".

        For example, given contracts A and B (represented by 2 DataFrames), each with 3 groups,
        we return (sub)commitments A1, A2, A3, B1, B2 and B3,
        where A,B,C is the enumerated contract and 1,2,3 is the enumerated group.
        """
        all_groups = []
        for df in dfs:
            df["j"] = range(len(df.index))
            groups = list(df["group"].unique())
            for group in groups:
                sub_commitment = df[df["group"] == group].drop(columns=["group"])
                all_groups.append(sub_commitment)
                # todo: raise if sub_commitment has non-unique prices or quantities
        return all_groups

    commitments: list[pd.DataFrame] = convert_commitments_to_subcommitments(commitments)

    bigM_columns = ["derivative max", "derivative min", "derivative equals"]
    # Compute a good value for M
    M = np.nanmax([np.nanmax(d[bigM_columns].abs()) for d in device_constraints])

    # M has to be 1 MW, at least
    M = max(M, 1)

    for d in range(len(device_constraints)):
        if "stock delta" not in device_constraints[d].columns:
            device_constraints[d]["stock delta"] = 0
        else:
            device_constraints[d]["stock delta"] = device_constraints[d][
                "stock delta"
            ].fillna(0)

    # Add indices for devices (d), datetimes (j) and commitments (c)
    model.d = RangeSet(0, len(device_constraints) - 1, doc="Set of devices")
    model.j = RangeSet(
        0, len(device_constraints[0].index.to_pydatetime()) - 1, doc="Set of datetimes"
    )
    model.c = RangeSet(0, len(commitments) - 1, doc="Set of commitments")

    # Add 2D indices for commitment datetimes (cj)

    def commitments_init(m):
        return ((c, j) for c in m.c for j in commitments[c]["j"])

    model.cj = Set(dimen=2, initialize=commitments_init)

    # Add parameters
    def price_down_select(m, c):
        return commitments[c]["downwards deviation price"].iloc[0]

    def price_up_select(m, c):
        return commitments[c]["upwards deviation price"].iloc[0]

    def commitment_quantity_select(m, c):
        return commitments[c]["quantity"].iloc[0]

    def device_max_select(m, d, j):
        min_v = device_constraints[d]["min"].iloc[j]
        max_v = device_constraints[d]["max"].iloc[j]
        equal_v = device_constraints[d]["equals"].iloc[j]
        if np.isnan(max_v) and np.isnan(equal_v):
            return infinity
        else:
            if not np.isnan(equal_v):
                # make min_v < equal_v
                equal_v = np.nanmax([equal_v, min_v])

            return np.nanmin([max_v, equal_v])

    def device_min_select(m, d, j):
        min_v = device_constraints[d]["min"].iloc[j]
        max_v = device_constraints[d]["max"].iloc[j]
        equal_v = device_constraints[d]["equals"].iloc[j]
        if np.isnan(min_v) and np.isnan(equal_v):
            return -infinity
        else:
            if not np.isnan(equal_v):
                # make equal_v <= max_v
                equal_v = np.nanmin([equal_v, max_v])

            return np.nanmax([min_v, equal_v])

    def device_derivative_max_select(m, d, j):
        max_v = device_constraints[d]["derivative max"].iloc[j]
        equal_v = device_constraints[d]["derivative equals"].iloc[j]
        if np.isnan(max_v) and np.isnan(equal_v):
            return infinity
        else:
            return np.nanmin([max_v, equal_v])

    def device_derivative_min_select(m, d, j):
        min_v = device_constraints[d]["derivative min"].iloc[j]
        equal_v = device_constraints[d]["derivative equals"].iloc[j]
        if np.isnan(min_v) and np.isnan(equal_v):
            return -infinity
        else:
            return np.nanmax([min_v, equal_v])

    def ems_derivative_max_select(m, j):
        v = ems_constraints["derivative max"].iloc[j]
        if np.isnan(v):
            return infinity
        else:
            return v

    def ems_derivative_min_select(m, j):
        v = ems_constraints["derivative min"].iloc[j]
        if np.isnan(v):
            return -infinity
        else:
            return v

    def device_efficiency(m, d, j):
        """Assume perfect efficiency if no efficiency information is available."""
        try:
            eff = device_constraints[d]["efficiency"].iloc[j]
        except KeyError:
            return 1
        if np.isnan(eff):
            return 1
        return eff

    def device_derivative_down_efficiency(m, d, j):
        """Assume perfect efficiency if no efficiency information is available."""
        try:
            eff = device_constraints[d]["derivative down efficiency"].iloc[j]
        except KeyError:
            return 1
        if np.isnan(eff):
            return 1
        return eff

    def device_derivative_up_efficiency(m, d, j):
        """Assume perfect efficiency if no efficiency information is available."""
        try:
            eff = device_constraints[d]["derivative up efficiency"].iloc[j]
        except KeyError:
            return 1
        if np.isnan(eff):
            return 1
        return eff

    def device_stock_delta(m, d, j):
        return device_constraints[d]["stock delta"].iloc[j]

    model.up_price = Param(model.c, initialize=price_up_select)
    model.down_price = Param(model.c, initialize=price_down_select)
    model.commitment_quantity = Param(model.c, initialize=commitment_quantity_select)
    model.device_max = Param(model.d, model.j, initialize=device_max_select)
    model.device_min = Param(model.d, model.j, initialize=device_min_select)
    model.device_derivative_max = Param(
        model.d, model.j, initialize=device_derivative_max_select
    )
    model.device_derivative_min = Param(
        model.d, model.j, initialize=device_derivative_min_select
    )
    model.ems_derivative_max = Param(model.j, initialize=ems_derivative_max_select)
    model.ems_derivative_min = Param(model.j, initialize=ems_derivative_min_select)
    model.device_efficiency = Param(model.d, model.j, initialize=device_efficiency)
    model.device_derivative_down_efficiency = Param(
        model.d, model.j, initialize=device_derivative_down_efficiency
    )
    model.device_derivative_up_efficiency = Param(
        model.d, model.j, initialize=device_derivative_up_efficiency
    )
    model.stock_delta = Param(model.d, model.j, initialize=device_stock_delta)

    # Add variables
    model.ems_power = Var(model.d, model.j, domain=Reals, initialize=0)
    model.device_power_down = Var(
        model.d, model.j, domain=NonPositiveReals, initialize=0
    )
    model.device_power_up = Var(model.d, model.j, domain=NonNegativeReals, initialize=0)
    model.device_power_sign = Var(model.d, model.j, domain=Binary, initialize=0)
    model.commitment_downwards_deviation = Var(
        model.c, domain=NonPositiveReals, initialize=0
    )
    model.commitment_upwards_deviation = Var(
        model.c, domain=NonNegativeReals, initialize=0
    )

    # Add constraints as a tuple of (lower bound, value, upper bound)
    def device_bounds(m, d, j):
        """Apply conversion efficiencies to conversion from flow to stock change and vice versa,
        and apply storage efficiencies to stock levels from one datetime to the next."""
        stock_changes = [
            (
                m.device_power_down[d, k] / m.device_derivative_down_efficiency[d, k]
                + m.device_power_up[d, k] * m.device_derivative_up_efficiency[d, k]
                + m.stock_delta[d, k]
            )
            for k in range(0, j + 1)
        ]
        efficiencies = [m.device_efficiency[d, k] for k in range(0, j + 1)]
        return (
            m.device_min[d, j],
            [
                stock - initial_stock
                for stock in apply_stock_changes_and_losses(
                    initial_stock, stock_changes, efficiencies
                )
            ][-1],
            m.device_max[d, j],
        )

    def device_derivative_bounds(m, d, j):
        return (
            m.device_derivative_min[d, j],
            m.device_power_down[d, j] + m.device_power_up[d, j],
            m.device_derivative_max[d, j],
        )

    def device_down_derivative_bounds(m, d, j):
        """Strictly non-positive."""
        return (
            min(m.device_derivative_min[d, j], 0),
            m.device_power_down[d, j],
            0,
        )

    def device_up_derivative_bounds(m, d, j):
        """Strictly non-negative."""
        return (
            0,
            m.device_power_up[d, j],
            max(0, m.device_derivative_max[d, j]),
        )

    def device_up_derivative_sign(m, d, j):
        """Derivative up if sign points up, derivative not up if sign points down."""
        return m.device_power_up[d, j] <= M * m.device_power_sign[d, j]

    def device_down_derivative_sign(m, d, j):
        """Derivative down if sign points down, derivative not down if sign points up."""
        return -m.device_power_down[d, j] <= M * (1 - m.device_power_sign[d, j])

    def ems_derivative_bounds(m, j):
        return m.ems_derivative_min[j], sum(m.ems_power[:, j]), m.ems_derivative_max[j]

    def ems_flow_commitment_equalities(m, c, j):
        """Couple EMS flows (sum over devices) to each commitment."""
        return (
            0,
            m.commitment_quantity[c]
            + m.commitment_downwards_deviation[c]
            + m.commitment_upwards_deviation[c]
            - sum(m.ems_power[:, j]),
            # take max deviation within commitment group
            0 if len(commitments[c]) == 1 else None,
        )

    def ems_flow_commitment_up_equalities(m, c, j):
        """Couple EMS flows (sum over devices) to each up commitment."""
        return (
            0,
            # sum(m.ems_power[:, j]) - m.commitment_quantity[c] - m.commitment_upwards_deviation[c],
            m.commitment_quantity[c]
            - sum(m.ems_power[:, j])
            + m.commitment_upwards_deviation[c],
            None,
            # sum(m.ems_power[:, j]) - m.commitment_quantity[c] < m.commitment_upwards_deviation[c]
        )

    def ems_flow_commitment_down_equalities(m, c, j):
        """Couple EMS flows (sum over devices) to each down commitment."""
        return (
            None,
            # sum(m.ems_power[:, j]) - m.commitment_quantity[c] - m.commitment_downwards_deviation[c],
            m.commitment_quantity[c]
            - sum(m.ems_power[:, j])
            + m.commitment_downwards_deviation[c],
            0,
            # m.commitment_downwards_deviation[c]  < sum(m.ems_power[:, j]) - m.commitment_quantity[c]
        )

    def device_derivative_equalities(m, d, j):
        """Couple device flows to EMS flows per device."""
        return (
            0,
            m.device_power_up[d, j] + m.device_power_down[d, j] - m.ems_power[d, j],
            0,
        )

    model.device_energy_bounds = Constraint(model.d, model.j, rule=device_bounds)
    model.device_power_bounds = Constraint(
        model.d, model.j, rule=device_derivative_bounds
    )
    model.device_power_down_bounds = Constraint(
        model.d, model.j, rule=device_down_derivative_bounds
    )
    model.device_power_up_bounds = Constraint(
        model.d, model.j, rule=device_up_derivative_bounds
    )
    model.device_power_up_sign = Constraint(
        model.d, model.j, rule=device_up_derivative_sign
    )
    model.device_power_down_sign = Constraint(
        model.d, model.j, rule=device_down_derivative_sign
    )
    model.ems_power_bounds = Constraint(model.j, rule=ems_derivative_bounds)
    model.ems_power_commitment_equalities = Constraint(
        model.cj, rule=ems_flow_commitment_equalities
    )
    # model.ems_power_commitment_down_equalities = Constraint(
    #     model.cj, rule=ems_flow_commitment_down_equalities
    # )
    # model.ems_power_commitment_up_equalities = Constraint(
    #     model.cj, rule=ems_flow_commitment_up_equalities
    # )
    model.device_power_equalities = Constraint(
        model.d, model.j, rule=device_derivative_equalities
    )

    # Add objective
    def cost_function(m):
        costs = 0
        for c in m.c:
            costs += m.commitment_downwards_deviation[c] * m.down_price[c]
            costs += m.commitment_upwards_deviation[c] * m.up_price[c]
        return costs

    model.costs = Objective(rule=cost_function, sense=minimize)

    # Solve
    solver_name = current_app.config.get("FLEXMEASURES_LP_SOLVER")

    solver = SolverFactory(solver_name)

    # disable logs for the HiGHS solver in case that LOGGING_LEVEL is INFO
    if current_app.config["LOGGING_LEVEL"] == "INFO" and (
        "highs" in solver_name.lower()
    ):
        solver.options["output_flag"] = "false"

    # load_solutions=False to avoid a RuntimeError exception in appsi solvers when solving an infeasible problem.
    results = solver.solve(model, load_solutions=False)

    # load the results only if a feasible solution has been found
    if len(results.solution) > 0:
        model.solutions.load_from(results)

    planned_costs = value(model.costs)
    planned_power_per_device = []
    for d in model.d:
        planned_device_power = [model.ems_power[d, j].value for j in model.j]
        planned_power_per_device.append(
            initialize_series(
                data=planned_device_power,
                start=start,
                end=end,
                resolution=to_offset(resolution),
            )
        )

    model.pprint()
    model.display()
    print(results.solver.termination_condition)
    print(planned_costs)
    return planned_power_per_device, planned_costs, results, model

1	from __future__ import annotations	8✔
2
3	from flask import current_app	8✔
4	import pandas as pd	8✔
5	import numpy as np	8✔
6	from pandas.tseries.frequencies import to_offset	8✔
7	from pyomo.core import (	8✔
8	ConcreteModel,
9	Var,
10	RangeSet,
11	Set,
12	Param,
13	Reals,
14	NonNegativeReals,
15	NonPositiveReals,
16	Binary,
17	Constraint,
18	Objective,
19	minimize,
20	)
21	from pyomo.environ import UnknownSolver # noqa F401	8✔
22	from pyomo.environ import value	8✔
23	from pyomo.opt import SolverFactory, SolverResults	8✔
24
25	from flexmeasures.data.models.planning.utils import initialize_series, initialize_df	8✔
26	from flexmeasures.utils.calculations import apply_stock_changes_and_losses	8✔
27
28	infinity = float("inf")	8✔
29
30
31	def device_scheduler( # noqa C901	8✔
32	device_constraints: list[pd.DataFrame],
33	ems_constraints: pd.DataFrame,
34	commitment_quantities: list[pd.Series] \| None = None,
35	commitment_downwards_deviation_price: list[pd.Series] \| list[float] \| None = None,
36	commitment_upwards_deviation_price: list[pd.Series] \| list[float] \| None = None,
37	commitments: list[pd.DataFrame] \| None = None,
38	initial_stock: float = 0,
39	) -> tuple[list[pd.Series], float, SolverResults, ConcreteModel]:
40	"""This generic device scheduler is able to handle an EMS with multiple devices,
41	with various types of constraints on the EMS level and on the device level,
42	and with multiple market commitments on the EMS level.
43	A typical example is a house with many devices.
44	The commitments are assumed to be with regard to the flow of energy to the device (positive for consumption,
45	negative for production). The solver minimises the costs of deviating from the commitments.
46
47	:param device_constraints: Device constraints are on a device level. Handled constraints (listed by column name):
48	max: maximum stock assuming an initial stock of zero (e.g. in MWh or boxes)
49	min: minimum stock assuming an initial stock of zero
50	equal: exact amount of stock (we do this by clamping min and max)
51	efficiency: amount of stock left at the next datetime (the rest is lost)
52	derivative max: maximum flow (e.g. in MW or boxes/h)
53	derivative min: minimum flow
54	derivative equals: exact amount of flow (we do this by clamping derivative min and derivative max)
55	derivative down efficiency: conversion efficiency of flow out of a device (flow out : stock decrease)
56	derivative up efficiency: conversion efficiency of flow into a device (stock increase : flow in)
57	stock delta: predefined stock delta to apply to the storage device. Positive values cause an increase and negative values a decrease
58	:param ems_constraints: EMS constraints are on an EMS level. Handled constraints (listed by column name):
59	derivative max: maximum flow
60	derivative min: minimum flow
61	:param commitments: Commitments are on an EMS level. Handled parameter (listed by column name):
62	quantity: for example, 5.5
63	downwards deviation price: 10.1
64	upwards deviation price: 10.2
65	group: 1 (defaults to the enumerate time step j)
66
67	Potentially deprecated arguments:
68	commitment_quantities: amounts of flow specified in commitments (both previously ordered and newly requested)
69	- e.g. in MW or boxes/h
70	commitment_downwards_deviation_price: penalty for downwards deviations of the flow
71	- e.g. in EUR/MW or EUR/(boxes/h)
72	- either a single value (same value for each flow value) or a Series (different value for each flow value)
73	commitment_upwards_deviation_price: penalty for upwards deviations of the flow
74
75	All Series and DataFrames should have the same resolution.
76
77	For now, we pass in the various constraints and prices as separate variables, from which we make a MultiIndex
78	DataFrame. Later we could pass in a MultiIndex DataFrame directly.
79	"""
80
81	model = ConcreteModel()	8✔
82
83	# If the EMS has no devices, don't bother
84	if len(device_constraints) == 0:	8✔
85	return [], 0, SolverResults(), model	×
86
87	# Get timing from first device
88	start = device_constraints[0].index.to_pydatetime()[0]	8✔
89	# Workaround for https://github.com/pandas-dev/pandas/issues/53643. Was: resolution = pd.to_timedelta(device_constraints[0].index.freq)
90	resolution = pd.to_timedelta(device_constraints[0].index.freq).to_pytimedelta()	8✔
91	end = device_constraints[0].index.to_pydatetime()[-1] + resolution	8✔
92
93	# Move commitments from old structure to new
94	if commitments is None:	8✔
NEW 95	commitments = []	×
96	if commitment_quantities is not None:	8✔
NEW 97	for quantity, down, up in zip(	×
98	commitment_quantities,
99	commitment_downwards_deviation_price,
100	commitment_upwards_deviation_price,
101	):
102
103	# Turn prices per commitment into prices per commitment flow
NEW 104	if all(isinstance(price, float) for price in down) or isinstance(	×
105	down, float
106	):
NEW 107	down = initialize_series(down, start, end, resolution)	×
NEW 108	if all(isinstance(price, float) for price in up) or isinstance(up, float):	×
NEW 109	up = initialize_series(up, start, end, resolution)	×
110
NEW 111	group = initialize_series(list(range(len(down))), start, end, resolution)	×
NEW 112	df = initialize_df(	×
113	["quantity", "downwards deviation price", "upwards deviation price"],
114	start,
115	end,
116	resolution,
117	)
NEW 118	df["quantity"] = quantity	×
NEW 119	df["downwards deviation price"] = down	×
NEW 120	df["upwards deviation price"] = up	×
NEW 121	df["group"] = group	×
NEW 122	commitments.append(df)	×
123
124	# Check if commitments have the same time window and resolution as the constraints
125	for commitment in commitments:	8!
126	start_c = commitment.index.to_pydatetime()[0]	8✔
127	resolution_c = pd.to_timedelta(commitment.index.freq)	8✔
128	end_c = commitment.index.to_pydatetime()[-1] + resolution	8✔
129	if not (start_c == start and end_c == end):	8✔
130	raise Exception(	×
131	"Not implemented for different time windows.\n(%s,%s)\n(%s,%s)"
132	% (start, end, start_c, end_c)
133	)
134	if resolution_c != resolution:	8✔
135	raise Exception(	×
136	"Not implemented for different resolutions.\n%s\n%s"
137	% (resolution, resolution_c)
138	)
139
140	def convert_commitments_to_subcommitments(	8✔
141	dfs: list[pd.DataFrame],
142	) -> list[pd.DataFrame]:
143	"""Transform commitments, each specifying a group for each time step, to sub-commitments, one per group.
144
145	We also enumerate the time steps in a new column "j".
146
147	For example, given contracts A and B (represented by 2 DataFrames), each with 3 groups,
148	we return (sub)commitments A1, A2, A3, B1, B2 and B3,
149	where A,B,C is the enumerated contract and 1,2,3 is the enumerated group.
150	"""
151	all_groups = []	8✔
152	for df in dfs:	8✔
153	df["j"] = range(len(df.index))	8✔
154	groups = list(df["group"].unique())	8✔
155	for group in groups:	8✔
156	sub_commitment = df[df["group"] == group].drop(columns=["group"])	8✔
157	all_groups.append(sub_commitment)	8✔
158	# todo: raise if sub_commitment has non-unique prices or quantities
159	return all_groups	8✔
160
161	commitments: list[pd.DataFrame] = convert_commitments_to_subcommitments(commitments)	8✔
162
163	bigM_columns = ["derivative max", "derivative min", "derivative equals"]	8✔
164	# Compute a good value for M
165	M = np.nanmax([np.nanmax(d[bigM_columns].abs()) for d in device_constraints])	8✔
166
167	# M has to be 1 MW, at least
168	M = max(M, 1)	8✔
169
170	for d in range(len(device_constraints)):	8✔
171	if "stock delta" not in device_constraints[d].columns:	8✔
172	device_constraints[d]["stock delta"] = 0	×
173	else:
174	device_constraints[d]["stock delta"] = device_constraints[d][	8✔
175	"stock delta"
176	].fillna(0)
177
178	# Add indices for devices (d), datetimes (j) and commitments (c)
179	model.d = RangeSet(0, len(device_constraints) - 1, doc="Set of devices")	8✔
180	model.j = RangeSet(	8✔
181	0, len(device_constraints[0].index.to_pydatetime()) - 1, doc="Set of datetimes"
182	)
183	model.c = RangeSet(0, len(commitments) - 1, doc="Set of commitments")	8✔
184
185	# Add 2D indices for commitment datetimes (cj)
186
187	def commitments_init(m):	8✔
188	return ((c, j) for c in m.c for j in commitments[c]["j"])	8✔
189
190	model.cj = Set(dimen=2, initialize=commitments_init)	8✔
191
192	# Add parameters
193	def price_down_select(m, c):	8✔
194	return commitments[c]["downwards deviation price"].iloc[0]	8✔
195
196	def price_up_select(m, c):	8✔
197	return commitments[c]["upwards deviation price"].iloc[0]	8✔
198
199	def commitment_quantity_select(m, c):	8✔
200	return commitments[c]["quantity"].iloc[0]	8✔
201
202	def device_max_select(m, d, j):	8✔
203	min_v = device_constraints[d]["min"].iloc[j]	8✔
204	max_v = device_constraints[d]["max"].iloc[j]	8✔
205	equal_v = device_constraints[d]["equals"].iloc[j]	8✔
206	if np.isnan(max_v) and np.isnan(equal_v):	8✔
207	return infinity	8✔
208	else:
209	if not np.isnan(equal_v):	8✔
210	# make min_v < equal_v
211	equal_v = np.nanmax([equal_v, min_v])	8✔
212
213	return np.nanmin([max_v, equal_v])	8✔
214
215	def device_min_select(m, d, j):	8✔
216	min_v = device_constraints[d]["min"].iloc[j]	8✔
217	max_v = device_constraints[d]["max"].iloc[j]	8✔
218	equal_v = device_constraints[d]["equals"].iloc[j]	8✔
219	if np.isnan(min_v) and np.isnan(equal_v):	8✔
220	return -infinity	8✔
221	else:
222	if not np.isnan(equal_v):	8✔
223	# make equal_v <= max_v
224	equal_v = np.nanmin([equal_v, max_v])	8✔
225
226	return np.nanmax([min_v, equal_v])	8✔
227
228	def device_derivative_max_select(m, d, j):	8✔
229	max_v = device_constraints[d]["derivative max"].iloc[j]	8✔
230	equal_v = device_constraints[d]["derivative equals"].iloc[j]	8✔
231	if np.isnan(max_v) and np.isnan(equal_v):	8✔
UNCOV 232	return infinity	×
233	else:
234	return np.nanmin([max_v, equal_v])	8✔
235
236	def device_derivative_min_select(m, d, j):	8✔
237	min_v = device_constraints[d]["derivative min"].iloc[j]	8✔
238	equal_v = device_constraints[d]["derivative equals"].iloc[j]	8✔
239	if np.isnan(min_v) and np.isnan(equal_v):	8✔
UNCOV 240	return -infinity	×
241	else:
242	return np.nanmax([min_v, equal_v])	8✔
243
244	def ems_derivative_max_select(m, j):	8✔
245	v = ems_constraints["derivative max"].iloc[j]	8✔
246	if np.isnan(v):	8✔
247	return infinity	8✔
248	else:
249	return v	8✔
250
251	def ems_derivative_min_select(m, j):	8✔
252	v = ems_constraints["derivative min"].iloc[j]	8✔
253	if np.isnan(v):	8✔
254	return -infinity	8✔
255	else:
256	return v	8✔
257
258	def device_efficiency(m, d, j):	8✔
259	"""Assume perfect efficiency if no efficiency information is available."""
260	try:	8✔
261	eff = device_constraints[d]["efficiency"].iloc[j]	8✔
UNCOV 262	except KeyError:	×
UNCOV 263	return 1	×
264	if np.isnan(eff):	8✔
265	return 1	8✔
266	return eff	8✔
267
268	def device_derivative_down_efficiency(m, d, j):	8✔
269	"""Assume perfect efficiency if no efficiency information is available."""
270	try:	8✔
271	eff = device_constraints[d]["derivative down efficiency"].iloc[j]	8✔
UNCOV 272	except KeyError:	×
UNCOV 273	return 1	×
274	if np.isnan(eff):	8✔
275	return 1	8✔
276	return eff	8✔
277
278	def device_derivative_up_efficiency(m, d, j):	8✔
279	"""Assume perfect efficiency if no efficiency information is available."""
280	try:	8✔
281	eff = device_constraints[d]["derivative up efficiency"].iloc[j]	8✔
UNCOV 282	except KeyError:	×
283	return 1	×
284	if np.isnan(eff):	8✔
285	return 1	8✔
286	return eff	8✔
287
288	def device_stock_delta(m, d, j):	8✔
289	return device_constraints[d]["stock delta"].iloc[j]	8✔
290
291	model.up_price = Param(model.c, initialize=price_up_select)	8✔
292	model.down_price = Param(model.c, initialize=price_down_select)	8✔
293	model.commitment_quantity = Param(model.c, initialize=commitment_quantity_select)	8✔
294	model.device_max = Param(model.d, model.j, initialize=device_max_select)	8✔
295	model.device_min = Param(model.d, model.j, initialize=device_min_select)	8✔
296	model.device_derivative_max = Param(	8✔
297	model.d, model.j, initialize=device_derivative_max_select
298	)
299	model.device_derivative_min = Param(	8✔
300	model.d, model.j, initialize=device_derivative_min_select
301	)
302	model.ems_derivative_max = Param(model.j, initialize=ems_derivative_max_select)	8✔
303	model.ems_derivative_min = Param(model.j, initialize=ems_derivative_min_select)	8✔
304	model.device_efficiency = Param(model.d, model.j, initialize=device_efficiency)	8✔
305	model.device_derivative_down_efficiency = Param(	8✔
306	model.d, model.j, initialize=device_derivative_down_efficiency
307	)
308	model.device_derivative_up_efficiency = Param(	8✔
309	model.d, model.j, initialize=device_derivative_up_efficiency
310	)
311	model.stock_delta = Param(model.d, model.j, initialize=device_stock_delta)	8✔
312
313	# Add variables
314	model.ems_power = Var(model.d, model.j, domain=Reals, initialize=0)	8✔
315	model.device_power_down = Var(	8✔
316	model.d, model.j, domain=NonPositiveReals, initialize=0
317	)
318	model.device_power_up = Var(model.d, model.j, domain=NonNegativeReals, initialize=0)	8✔
319	model.device_power_sign = Var(model.d, model.j, domain=Binary, initialize=0)	8✔
320	model.commitment_downwards_deviation = Var(	8✔
321	model.c, domain=NonPositiveReals, initialize=0
322	)
323	model.commitment_upwards_deviation = Var(	8✔
324	model.c, domain=NonNegativeReals, initialize=0
325	)
326
327	# Add constraints as a tuple of (lower bound, value, upper bound)
328	def device_bounds(m, d, j):	8✔
329	"""Apply conversion efficiencies to conversion from flow to stock change and vice versa,
330	and apply storage efficiencies to stock levels from one datetime to the next."""
331	stock_changes = [	8✔
332	(
333	m.device_power_down[d, k] / m.device_derivative_down_efficiency[d, k]
334	+ m.device_power_up[d, k] * m.device_derivative_up_efficiency[d, k]
335	+ m.stock_delta[d, k]
336	)
337	for k in range(0, j + 1)
338	]
339	efficiencies = [m.device_efficiency[d, k] for k in range(0, j + 1)]	8✔
340	return (	8✔
341	m.device_min[d, j],
342	[
343	stock - initial_stock
344	for stock in apply_stock_changes_and_losses(
345	initial_stock, stock_changes, efficiencies
346	)
347	][-1],
348	m.device_max[d, j],
349	)
350
351	def device_derivative_bounds(m, d, j):	8✔
352	return (	8✔
353	m.device_derivative_min[d, j],
354	m.device_power_down[d, j] + m.device_power_up[d, j],
355	m.device_derivative_max[d, j],
356	)
357
358	def device_down_derivative_bounds(m, d, j):	8✔
359	"""Strictly non-positive."""
360	return (	8✔
361	min(m.device_derivative_min[d, j], 0),
362	m.device_power_down[d, j],
363	0,
364	)
365
366	def device_up_derivative_bounds(m, d, j):	8✔
367	"""Strictly non-negative."""
368	return (	8✔
369	0,
370	m.device_power_up[d, j],
371	max(0, m.device_derivative_max[d, j]),
372	)
373
374	def device_up_derivative_sign(m, d, j):	8✔
375	"""Derivative up if sign points up, derivative not up if sign points down."""
376	return m.device_power_up[d, j] <= M * m.device_power_sign[d, j]	8✔
377
378	def device_down_derivative_sign(m, d, j):	8✔
379	"""Derivative down if sign points down, derivative not down if sign points up."""
380	return -m.device_power_down[d, j] <= M * (1 - m.device_power_sign[d, j])	8✔
381
382	def ems_derivative_bounds(m, j):	8✔
383	return m.ems_derivative_min[j], sum(m.ems_power[:, j]), m.ems_derivative_max[j]	8✔
384
385	def ems_flow_commitment_equalities(m, c, j):	8✔
386	"""Couple EMS flows (sum over devices) to each commitment."""
387	return (	8✔
388	0,
389	m.commitment_quantity[c]
390	+ m.commitment_downwards_deviation[c]
391	+ m.commitment_upwards_deviation[c]
392	- sum(m.ems_power[:, j]),
393	# take max deviation within commitment group
394	0 if len(commitments[c]) == 1 else None,
395	)
396
397	def ems_flow_commitment_up_equalities(m, c, j):	8✔
398	"""Couple EMS flows (sum over devices) to each up commitment."""
NEW UNCOV 399	return (	×
400	0,
401	# sum(m.ems_power[:, j]) - m.commitment_quantity[c] - m.commitment_upwards_deviation[c],
402	m.commitment_quantity[c]
403	- sum(m.ems_power[:, j])
404	+ m.commitment_upwards_deviation[c],
405	None,
406	# sum(m.ems_power[:, j]) - m.commitment_quantity[c] < m.commitment_upwards_deviation[c]
407	)
408
409	def ems_flow_commitment_down_equalities(m, c, j):	8✔
410	"""Couple EMS flows (sum over devices) to each down commitment."""
NEW UNCOV 411	return (	×
412	None,
413	# sum(m.ems_power[:, j]) - m.commitment_quantity[c] - m.commitment_downwards_deviation[c],
414	m.commitment_quantity[c]
415	- sum(m.ems_power[:, j])
416	+ m.commitment_downwards_deviation[c],
417	0,
418	# m.commitment_downwards_deviation[c] < sum(m.ems_power[:, j]) - m.commitment_quantity[c]
419	)
420
421	def device_derivative_equalities(m, d, j):	8✔
422	"""Couple device flows to EMS flows per device."""
423	return (	8✔
424	0,
425	m.device_power_up[d, j] + m.device_power_down[d, j] - m.ems_power[d, j],
426	0,
427	)
428
429	model.device_energy_bounds = Constraint(model.d, model.j, rule=device_bounds)	8✔
430	model.device_power_bounds = Constraint(	8✔
431	model.d, model.j, rule=device_derivative_bounds
432	)
433	model.device_power_down_bounds = Constraint(	8✔
434	model.d, model.j, rule=device_down_derivative_bounds
435	)
436	model.device_power_up_bounds = Constraint(	8✔
437	model.d, model.j, rule=device_up_derivative_bounds
438	)
439	model.device_power_up_sign = Constraint(	8✔
440	model.d, model.j, rule=device_up_derivative_sign
441	)
442	model.device_power_down_sign = Constraint(	8✔
443	model.d, model.j, rule=device_down_derivative_sign
444	)
445	model.ems_power_bounds = Constraint(model.j, rule=ems_derivative_bounds)	8✔
446	model.ems_power_commitment_equalities = Constraint(	8✔
447	model.cj, rule=ems_flow_commitment_equalities
448	)
449	# model.ems_power_commitment_down_equalities = Constraint(
450	# model.cj, rule=ems_flow_commitment_down_equalities
451	# )
452	# model.ems_power_commitment_up_equalities = Constraint(
453	# model.cj, rule=ems_flow_commitment_up_equalities
454	# )
455	model.device_power_equalities = Constraint(	8✔
456	model.d, model.j, rule=device_derivative_equalities
457	)
458
459	# Add objective
460	def cost_function(m):	8✔
461	costs = 0	8✔
462	for c in m.c:	8✔
463	costs += m.commitment_downwards_deviation[c] * m.down_price[c]	8✔
464	costs += m.commitment_upwards_deviation[c] * m.up_price[c]	8✔
465	return costs	8✔
466
467	model.costs = Objective(rule=cost_function, sense=minimize)	8✔
468
469	# Solve
470	solver_name = current_app.config.get("FLEXMEASURES_LP_SOLVER")	8✔
471
472	solver = SolverFactory(solver_name)	8✔
473
474	# disable logs for the HiGHS solver in case that LOGGING_LEVEL is INFO
475	if current_app.config["LOGGING_LEVEL"] == "INFO" and (	8✔
476	"highs" in solver_name.lower()
477	):
UNCOV 478	solver.options["output_flag"] = "false"	×
479
480	# load_solutions=False to avoid a RuntimeError exception in appsi solvers when solving an infeasible problem.
481	results = solver.solve(model, load_solutions=False)	8✔
482
483	# load the results only if a feasible solution has been found
484	if len(results.solution) > 0:	8✔
485	model.solutions.load_from(results)	8✔
486
487	planned_costs = value(model.costs)	8✔
488	planned_power_per_device = []	8✔
489	for d in model.d:	8✔
490	planned_device_power = [model.ems_power[d, j].value for j in model.j]	8✔
491	planned_power_per_device.append(	8✔
492	initialize_series(
493	data=planned_device_power,
494	start=start,
495	end=end,
496	resolution=to_offset(resolution),
497	)
498	)
499
500	model.pprint()	8✔
501	model.display()	8✔
502	print(results.solver.termination_condition)	8✔
503	print(planned_costs)	8✔
504	return planned_power_per_device, planned_costs, results, model	8✔

FlexMeasures / flexmeasures / 10385477682

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous