Lucas Gautheron 1 vuosi sitten
vanhempi
commit
3b42f89a02
4 muutettua tiedostoa jossa 90 lisäystä ja 77 poistoa
  1. 44 33
      mix_simul/consumption.py
  2. 18 23
      mix_simul/production.py
  3. 9 9
      mix_simul/scenarios.py
  4. 19 12
      mix_simul/storage.py

+ 44 - 33
mix_simul/consumption.py

@@ -6,6 +6,7 @@ from statsmodels.formula.api import ols
 
 from typing import Union
 
+
 class ConsumptionModel:
     def __init__(self):
         pass
@@ -13,6 +14,7 @@ class ConsumptionModel:
     def get(self):
         pass
 
+
 class ThermoModel(ConsumptionModel):
     def __init__(self):
         pass
@@ -20,6 +22,7 @@ class ThermoModel(ConsumptionModel):
     def get(temperatures: np.ndarray):
         pass
 
+
 class FittedConsumptionModel(ConsumptionModel):
     def __init__(self, yearly_total: float):
         self.yearly_total = yearly_total
@@ -27,58 +30,66 @@ class FittedConsumptionModel(ConsumptionModel):
 
     def get(self, times: Union[pd.Series, np.ndarray]):
         # compute the load for the desired timestamps
-        times = pd.DataFrame({'time': times}).set_index("time")
+        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
+        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)
+            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["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()
+        hourly = consumption.resample("1H").mean()
 
         self.time_reference = hourly.index.min()
 
-        hourly['h'] = ((hourly.index - self.time_reference
-                        ).total_seconds()/3600).astype(int)
+        hourly["h"] = (
+            (hourly.index - self.time_reference).total_seconds() / 3600
+        ).astype(int)
 
-        hourly['Y'] = hourly.index.year-hourly.index.year.min()
+        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))
+        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)
+            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}']
+            components += [f"c_{i+1}", f"s_{i+1}"]
 
         hourly.rename(
-            columns={'Consommation brute électricité (MW) - RTE': 'conso'}, inplace=True)
+            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)
+        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)
+        intercept = self.fit_results.params[0] + self.fit_results.params[1]
+        self.fit_results.params *= self.yearly_total / (intercept * 365.25 * 24)
 
-class ConsumptionFlexibilityModel:
 
+class ConsumptionFlexibilityModel:
     def __init__(self, flexibility_power: float, flexibility_time: int):
         self.flexibility_power = flexibility_power
         self.flexibility_time = flexibility_time
@@ -87,28 +98,28 @@ class ConsumptionFlexibilityModel:
         from functools import reduce
 
         T = len(supply)
-        tau = self.flexibility_time    
+        tau = self.flexibility_time
 
-        h = cp.Variable((T, tau+1))
+        h = cp.Variable((T, tau + 1))
 
         constraints = [
             h >= 0,
             h <= 1,
-            h[:, 0] >= 1-self.flexibility_power/load, # bound on the fraction of the demand at any time that can be postponed
-            cp.multiply(load, cp.sum(h, axis=1))-load <= self.flexibility_power, 
+            # 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,
         ] + [
-            reduce(lambda x, y: x+y, [h[t-l, l] for l in range(tau)]) == 1 # total demand conservation
+            # 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
+            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)
+        return load * np.sum(hb, axis=1)

+ 18 - 23
mix_simul/production.py

@@ -1,6 +1,7 @@
 import numpy as np
 import cvxpy as cp
 
+
 class PowerSupply:
     def __init__(self):
         pass
@@ -10,22 +11,17 @@ class PowerSupply:
 
 
 class IntermittentArray(PowerSupply):
-    def __init__(self,
-        potential: np.ndarray,
-        units_per_region: np.ndarray
-    ):
+    def __init__(self, potential: np.ndarray, units_per_region: np.ndarray):
 
         self.potential = potential
         self.units_per_region = units_per_region
 
     def power(self):
-        return np.einsum('ijk,ik', self.potential, self.units_per_region)    
+        return np.einsum("ijk,ik", self.potential, self.units_per_region)
 
 
 class DispatchableArray(PowerSupply):
-    def __init__(self,
-        dispatchable: np.ndarray
-    ):
+    def __init__(self, dispatchable: np.ndarray):
         self.dispatchable = np.array(dispatchable)
         self.n_dispatchable_sources = self.dispatchable.shape[0]
 
@@ -34,25 +30,24 @@ class DispatchableArray(PowerSupply):
         T = len(gap)
         dispatch_power = cp.Variable((self.n_dispatchable_sources, T))
 
-        constraints = [
-            dispatch_power >= 0,
-            cp.sum(dispatch_power, axis=0) <= np.maximum(gap, 0)
-        ] + [
-            dispatch_power[i,:] <= self.dispatchable[i,0]
-            for i in range(self.n_dispatchable_sources)
-        ] + [
-            cp.sum(dispatch_power[i]) <= self.dispatchable[i,1]*T/(365.25*24)
-            for i in range(self.n_dispatchable_sources)
-        ]
+        constraints = (
+            [dispatch_power >= 0, cp.sum(dispatch_power, axis=0) <= np.maximum(gap, 0)]
+            + [
+                dispatch_power[i, :] <= self.dispatchable[i, 0]
+                for i in range(self.n_dispatchable_sources)
+            ]
+            + [
+                cp.sum(dispatch_power[i]) <= self.dispatchable[i, 1] * T / (365.25 * 24)
+                for i in range(self.n_dispatchable_sources)
+            ]
+        )
 
         prob = cp.Problem(
-            cp.Minimize(
-                cp.sum(cp.pos(gap-cp.sum(dispatch_power, axis=0)))
-            ),
-            constraints
+            cp.Minimize(cp.sum(cp.pos(gap - cp.sum(dispatch_power, axis=0)))),
+            constraints,
         )
 
         prob.solve(solver=cp.ECOS, max_iters=300)
         dp = dispatch_power.value
 
-        return dp
+        return dp

+ 9 - 9
mix_simul/scenarios.py

@@ -6,12 +6,13 @@ import yaml
 
 
 class Scenario:
-    def __init__(self,
-        yearly_total = 645*1000,
+    def __init__(
+        self,
+        yearly_total=645 * 1000,
         sources: dict = {},
         multistorage: dict = {},
         flexibility_power=0,
-        flexibility_time=8
+        flexibility_time=8,
     ):
         self.yearly_total = yearly_total
         self.sources = sources
@@ -27,8 +28,7 @@ class Scenario:
 
         # intermittent sources (or sources with fixed load factors)
         intermittent_array = IntermittentArray(
-            intermittent_load_factors,
-            np.transpose([self.sources["intermittent"]])
+            intermittent_load_factors, np.transpose([self.sources["intermittent"]])
         )
         power = intermittent_array.power()
 
@@ -38,18 +38,18 @@ class Scenario:
             )
             load = flexibility_model.run(load, power)
 
-        power_delta = power-load
+        power_delta = power - load
 
         # adjust power to load with storage
         storage_model = MultiStorageModel(
             self.multistorage["capacity"],
             self.multistorage["power"],
             self.multistorage["power"],
-            self.multistorage["efficiency"]
+            self.multistorage["efficiency"],
         )
 
         storage, storage_impact = storage_model.run(power_delta)
-        gap = load-power-storage_impact.sum(axis=0)
+        gap = load - power - storage_impact.sum(axis=0)
 
         # further adjust power to load with dispatchable power sources
         dispatchable_model = DispatchableArray(self.sources["dispatchable"])
@@ -58,4 +58,4 @@ class Scenario:
         gap -= dp.sum(axis=0)
         S = np.maximum(gap, 0).mean()
 
-        return S, load, power, gap, storage, dp
+        return S, load, power, gap, storage, dp

+ 19 - 12
mix_simul/storage.py

@@ -1,7 +1,8 @@
-import numpy as np 
+import numpy as np
 import numba
 from typing import Tuple
 
+
 @numba.njit
 def storage_iterate(dE, capacity, efficiency, n):
     storage = np.zeros(n)
@@ -10,26 +11,28 @@ def storage_iterate(dE, capacity, efficiency, n):
         if dE[i] >= 0:
             dE[i] *= efficiency
 
-        storage[i] = np.maximum(0, np.minimum(capacity, storage[i-1]+dE[i-1]))
+        storage[i] = np.maximum(0, np.minimum(capacity, storage[i - 1] + dE[i - 1]))
 
     return storage
 
+
 class StorageModel:
     def __int__(self):
         pass
 
 
 class MultiStorageModel(StorageModel):
-    def __init__(self,
+    def __init__(
+        self,
         storage_capacities: np.ndarray,
         storage_max_loads=np.ndarray,
         storage_max_deliveries=np.ndarray,
-        storage_efficiencies=np.ndarray
+        storage_efficiencies=np.ndarray,
     ):
 
         self.storage_max_loads = np.array(storage_max_loads)
         self.storage_max_deliveries = np.array(storage_max_deliveries)
-        self.storage_capacities = np.array(storage_capacities)*self.storage_max_loads
+        self.storage_capacities = np.array(storage_capacities) * self.storage_max_loads
         self.storage_efficiencies = np.array(storage_efficiencies)
 
         self.n_storages = len(self.storage_capacities)
@@ -38,7 +41,6 @@ class MultiStorageModel(StorageModel):
         assert len(storage_max_deliveries) == self.n_storages
         assert len(storage_efficiencies) == self.n_storages
 
-    
     def run(self, power_delta: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
         T = len(power_delta)
 
@@ -55,23 +57,28 @@ class MultiStorageModel(StorageModel):
 
         for i in range(self.n_storages):
             storage_try_load[i] = np.minimum(excess_power, self.storage_max_loads[i])
-            storage_try_delivery[i] = np.minimum(deficit_power, self.storage_max_deliveries[i])
+            storage_try_delivery[i] = np.minimum(
+                deficit_power, self.storage_max_deliveries[i]
+            )
 
-            dE_storage[i] = storage_try_load[i]-storage_try_delivery[i]
+            dE_storage[i] = storage_try_load[i] - storage_try_delivery[i]
 
             storage[i] = storage_iterate(
-                dE_storage[i], self.storage_capacities[i], self.storage_efficiencies[i], T
+                dE_storage[i],
+                self.storage_capacities[i],
+                self.storage_efficiencies[i],
+                T,
             )
 
-            # impact of storage on the available power        
+            # impact of storage on the available power
             storage_impact[i] = -np.diff(storage[i], append=0)
             storage_impact[i] = np.multiply(
                 storage_impact[i],
-                np.where(storage_impact[i] < 0, 1/self.storage_efficiencies[i], 1)
+                np.where(storage_impact[i] < 0, 1 / self.storage_efficiencies[i], 1),
             )
 
             available_power += storage_impact[i]
             excess_power = np.maximum(0, available_power)
             deficit_power = np.maximum(0, -available_power)
 
-        return storage, storage_impact
+        return storage, storage_impact