123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159 |
- import pandas as pd
- import numpy as np
- import cvxpy as cp
- from statsmodels.formula.api import ols
- from typing import Union
- class ConsumptionModel:
- def __init__(self):
- pass
- def get(self):
- pass
- class ThermoModel(ConsumptionModel):
- def __init__(self):
- pass
- def get(temperatures: np.ndarray):
- pass
- class FittedConsumptionModel(ConsumptionModel):
- def __init__(self, yearly_total: float):
- """Consumption Model fitted against observations of French electricity consumption (2012-2021).
- The model captures seasonal, weekly and diurnal variations. It has only 1 parameter, the total annual consumption.
- :param yearly_total: Total annual consumption in GWh.
- :type yearly_total: float
- """
- self.yearly_total = yearly_total
- self.fit()
- def get(self, times: Union[pd.Series, np.ndarray]) -> np.ndarray:
- """Retrieve the consumption for each timestamp from the input array.
- :param times: 1D array containing the input timestamps
- :type times: Union[pd.Series, np.ndarray]
- :return: 1D array of floats (consumption in GW) with the same length as the input.
- :rtype: np.ndarray
- """
- # compute the load for the desired timestamps
- times = pd.DataFrame({"time": times}).set_index("time")
- times.index = pd.to_datetime(times.index, utc=True)
- times["h"] = (
- (times.index - self.time_reference).total_seconds() / 3600
- ).astype(int)
- times["Y"] = 1
- for i, f in enumerate(self.frequencies):
- times[f"c_{i+1}"] = (times["h"] * f * 2 * np.pi).apply(np.cos)
- times[f"s_{i+1}"] = (times["h"] * f * 2 * np.pi).apply(np.sin)
- curve = self.fit_results.predict(times).values
- return curve
- def fit(self):
- consumption = pd.read_csv("data/consommation-quotidienne-brute.csv", sep=";")
- consumption["time"] = pd.to_datetime(
- consumption["Date - Heure"].str.replace(":", ""),
- format="%Y-%m-%dT%H%M%S%z",
- utc=True,
- )
- consumption.set_index("time", inplace=True)
- hourly = consumption.resample("1H").mean()
- self.time_reference = hourly.index.min()
- hourly["h"] = (
- (hourly.index - self.time_reference).total_seconds() / 3600
- ).astype(int)
- hourly["Y"] = hourly.index.year - hourly.index.year.min()
- # generate fourier components
- self.frequencies = (
- list((1 / (365.25 * 24)) * np.arange(1, 13))
- + list((1 / (24 * 7)) * np.arange(1, 7))
- + list((1 / 24) * np.arange(1, 12))
- )
- components = []
- for i, f in enumerate(self.frequencies):
- hourly[f"c_{i+1}"] = (hourly["h"] * f * 2 * np.pi).apply(np.cos)
- hourly[f"s_{i+1}"] = (hourly["h"] * f * 2 * np.pi).apply(np.sin)
- components += [f"c_{i+1}", f"s_{i+1}"]
- hourly.rename(
- columns={"Consommation brute électricité (MW) - RTE": "conso"}, inplace=True
- )
- hourly["conso"] /= 1000
- # fit load curve to fourier components
- model = ols("conso ~ " + " + ".join(components) + " + C(Y)", data=hourly)
- self.fit_results = model.fit()
- # normalize according to the desired total yearly consumption
- intercept = self.fit_results.params[0] + self.fit_results.params[1]
- self.fit_results.params *= self.yearly_total / (intercept * 365.25 * 24)
- class ConsumptionFlexibilityModel:
- def __init__(self, flexibility_power: float, flexibility_time: int):
- """Consumption flexibility model.
- Adapts consumption to supply under constraints.
- :param flexibility_power: Maximum power that can be postponed at any time, in GW
- :type flexibility_power: float
- :param flexibility_time: Maximum consumption delay in hours
- :type flexibility_time: int
- """
- self.flexibility_power = flexibility_power
- self.flexibility_time = flexibility_time
- def run(self, load: np.ndarray, supply: np.ndarray) -> np.ndarray:
- """Runs the model.
- Given initial load and supply, the model returns an adjusted load optimized under constraints.
- :param load: 1D array containing the initial load at each timestep, in GW
- :type load: np.ndarray
- :param supply: 1D array containing the power supply at each timestep, in GW
- :type supply: np.ndarray
- :return: 1D array containing the adjusted load at each timestep, in GW
- :rtype: np.ndarray
- """
- from functools import reduce
- T = len(supply)
- tau = self.flexibility_time
- h = cp.Variable((T, tau + 1))
- constraints = [
- h >= 0,
- h <= 1,
- # bound on the fraction of the demand at any time that can be postponed
- h[:, 0] >= 1 - self.flexibility_power / load,
- cp.multiply(load, cp.sum(h, axis=1)) - load <= self.flexibility_power,
- ] + [
- # total demand conservation
- reduce(lambda x, y: x + y, [h[t - l, l] for l in range(tau)]) == 1
- for t in np.arange(tau, T)
- ]
- prob = cp.Problem(
- cp.Minimize(cp.sum(cp.pos(cp.multiply(load, cp.sum(h, axis=1)) - supply))),
- constraints,
- )
- prob.solve(verbose=True, solver=cp.ECOS, max_iters=300)
- hb = np.array(h.value)
- return load * np.sum(hb, axis=1)
|