|
|
@@ -0,0 +1,596 @@
|
|
|
+import pandas as pd
|
|
|
+import numpy as np
|
|
|
+import numba
|
|
|
+
|
|
|
+from numpy.random import dirichlet
|
|
|
+from scipy.optimize import linprog
|
|
|
+from scipy.signal import square
|
|
|
+from statsmodels.formula.api import ols
|
|
|
+
|
|
|
+import cvxpy as cp
|
|
|
+
|
|
|
+from matplotlib import pyplot as plt
|
|
|
+import matplotlib.dates as mdates
|
|
|
+
|
|
|
+import argparse
|
|
|
+
|
|
|
+parser = argparse.ArgumentParser()
|
|
|
+parser.add_argument("--begin")
|
|
|
+parser.add_argument("--end")
|
|
|
+parser.add_argument("--flexibility", action="store_true")
|
|
|
+args = parser.parse_args()
|
|
|
+
|
|
|
+
|
|
|
+def achieved_load_factors(sources=[]):
|
|
|
+ # installed capacity
|
|
|
+ infra = pd.read_csv(
|
|
|
+ "registre-national-installation-production-stockage-electricite-agrege.csv", sep=";")
|
|
|
+ infra["Year"] = infra["dateMiseEnService"].str[-4:].astype(str)
|
|
|
+ infra.dropna(subset=["Year", "puisMaxInstallee"], inplace=True)
|
|
|
+ infra = infra[infra["Year"] != 'nan']
|
|
|
+ infra = infra.groupby(["filiere", "Year"]).agg(
|
|
|
+ pi=('puisMaxInstallee', 'sum')
|
|
|
+ ).reset_index().sort_values(["filiere", "Year"])
|
|
|
+
|
|
|
+ infra['pi'] = infra.groupby(
|
|
|
+ 'filiere')['pi'].transform(pd.Series.cumsum)/1000
|
|
|
+ infra["Year"] = infra["Year"].astype(int)
|
|
|
+ infra.set_index(["filiere", "Year"], inplace=True)
|
|
|
+ infra = infra.reindex(
|
|
|
+ index=[
|
|
|
+ (filiere, year) for filiere in set(infra.index.get_level_values(0)) for year in np.arange(2000, 2022)
|
|
|
+ ]
|
|
|
+ )
|
|
|
+ infra = infra.sort_index(level=0).ffill().reindex(infra.index)
|
|
|
+ infra.reset_index(inplace=True)
|
|
|
+ infra = infra.pivot(index="Year", columns="filiere", values="pi")
|
|
|
+
|
|
|
+ # production
|
|
|
+ prod = pd.read_csv("eco2mix-national-cons-def.csv", sep=";")
|
|
|
+ prod["time"] = pd.to_datetime(
|
|
|
+ prod["Date et Heure"].str.replace(":", ""), format="%Y-%m-%dT%H%M%S%z", utc=True
|
|
|
+ )
|
|
|
+ prod["Year"] = prod["time"].dt.year
|
|
|
+ prod = prod.merge(infra, left_on="Year", right_on="Year", how="left").sort_values("time").set_index("time")
|
|
|
+
|
|
|
+ # consumption
|
|
|
+ consumption = pd.read_csv("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)
|
|
|
+
|
|
|
+ print(consumption["Consommation brute électricité (MW) - RTE"])
|
|
|
+
|
|
|
+ for source in sources:
|
|
|
+ fig, ax1 = plt.subplots()
|
|
|
+ ax2 = ax1.twinx()
|
|
|
+ prod[f"{source} (load factor)"] = prod[f"{source} (MW)"]/prod[source]
|
|
|
+ print(prod[f"{source} (load factor)"])
|
|
|
+ ax1.scatter(prod.index, prod[f"{source} (load factor)"], s=0.5, color="red")
|
|
|
+ ax2.plot(consumption.index, consumption["Consommation brute électricité (MW) - RTE"], lw=0.5, color="blue")
|
|
|
+ plt.show()
|
|
|
+ fig.savefig(f"test_{source}.png", dpi=200)
|
|
|
+ plt.clf()
|
|
|
+ plt.cla()
|
|
|
+
|
|
|
+ # plt.show()
|
|
|
+
|
|
|
+
|
|
|
+# achieved_load_factors(["Nucléaire", "Hydraulique"])
|
|
|
+
|
|
|
+
|
|
|
+def generate_load_curve(times, yearly_total=645*1000):
|
|
|
+ # load observed load curve
|
|
|
+ consumption = pd.read_csv("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()
|
|
|
+ hourly['h'] = ((hourly.index - hourly.index.min()
|
|
|
+ ).total_seconds()/3600).astype(int)
|
|
|
+
|
|
|
+ hourly['Y'] = hourly.index.year-hourly.index.year.min()
|
|
|
+
|
|
|
+ # generate fourier components
|
|
|
+ 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(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)
|
|
|
+ results = model.fit()
|
|
|
+
|
|
|
+ # normalize according to the desired total yearly consumption
|
|
|
+ intercept = results.params[0]+results.params[1]
|
|
|
+ results.params *= yearly_total/(intercept*365.25*24)
|
|
|
+
|
|
|
+ # 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 - hourly.index.min()
|
|
|
+ ).total_seconds()/3600).astype(int)
|
|
|
+ times['Y'] = 1
|
|
|
+ for i, f in enumerate(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 = results.predict(times)
|
|
|
+ return np.array(curve)
|
|
|
+
|
|
|
+
|
|
|
+@numba.njit
|
|
|
+def storage_iterate(dE, capacity, efficiency, n):
|
|
|
+ storage = np.zeros(n)
|
|
|
+
|
|
|
+ for i in np.arange(1, n):
|
|
|
+ if dE[i] >= 0:
|
|
|
+ dE[i] *= efficiency
|
|
|
+
|
|
|
+ storage[i] = np.maximum(0, np.minimum(capacity, storage[i-1]+dE[i-1]))
|
|
|
+
|
|
|
+ return storage
|
|
|
+
|
|
|
+
|
|
|
+def objective_with_storages(
|
|
|
+ potential,
|
|
|
+ units_per_region,
|
|
|
+ dispatchable,
|
|
|
+ p_load=2664*1000/(365*24),
|
|
|
+ storage_capacities=5*250,
|
|
|
+ storage_max_loads=5*30,
|
|
|
+ storage_max_deliveries=5*30,
|
|
|
+ storage_efficiencies=0.67
|
|
|
+):
|
|
|
+
|
|
|
+ power = np.einsum('ijk,ik', potential, units_per_region)
|
|
|
+ power_delta = power-p_load
|
|
|
+
|
|
|
+ T = len(power)
|
|
|
+
|
|
|
+ available_power = np.array(power_delta)
|
|
|
+ excess_power = np.maximum(0, available_power)
|
|
|
+ deficit_power = np.maximum(0, -available_power)
|
|
|
+
|
|
|
+ n_storages = len(storage_capacities)
|
|
|
+ storage_try_load = np.zeros((n_storages, T))
|
|
|
+ storage_try_delivery = np.zeros((n_storages, T))
|
|
|
+
|
|
|
+ storage = np.zeros((n_storages, T))
|
|
|
+ storage_impact = np.zeros((n_storages, T))
|
|
|
+ dE_storage = np.zeros((n_storages, T))
|
|
|
+
|
|
|
+ for i in range(n_storages):
|
|
|
+ storage_try_load[i] = np.minimum(excess_power, storage_max_loads[i])
|
|
|
+ storage_try_delivery[i] = np.minimum(deficit_power, storage_max_deliveries[i])
|
|
|
+
|
|
|
+ dE_storage[i] = storage_try_load[i]-storage_try_delivery[i]
|
|
|
+
|
|
|
+ storage[i] = storage_iterate(
|
|
|
+ dE_storage[i], storage_capacities[i], storage_efficiencies[i], T
|
|
|
+ )
|
|
|
+
|
|
|
+ # 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/storage_efficiencies[i], 1)
|
|
|
+ )
|
|
|
+
|
|
|
+ available_power += storage_impact[i]
|
|
|
+ excess_power = np.maximum(0, available_power)
|
|
|
+ deficit_power = np.maximum(0, -available_power)
|
|
|
+
|
|
|
+ gap = p_load-power-storage_impact.sum(axis=0)
|
|
|
+
|
|
|
+ # optimize dispatch
|
|
|
+ n_dispatchable_sources = dispatchable.shape[0]
|
|
|
+ dispatch_power = cp.Variable((n_dispatchable_sources, T))
|
|
|
+
|
|
|
+ constraints = [
|
|
|
+ dispatch_power >= 0,
|
|
|
+ cp.sum(dispatch_power, axis=0) <= np.maximum(gap, 0)
|
|
|
+ ] + [
|
|
|
+ dispatch_power[i,:] <= dispatchable[i,:,0].sum()
|
|
|
+ for i in range(n_dispatchable_sources)
|
|
|
+ ] + [
|
|
|
+ cp.sum(dispatch_power[i]) <= dispatchable[i,:,1].sum()
|
|
|
+ for i in range(n_dispatchable_sources)
|
|
|
+ ]
|
|
|
+
|
|
|
+ prob = cp.Problem(
|
|
|
+ cp.Minimize(
|
|
|
+ cp.sum(cp.pos(gap-cp.sum(dispatch_power, axis=0))) + cp.max(gap-cp.sum(dispatch_power, axis=0))
|
|
|
+ ),
|
|
|
+ constraints
|
|
|
+ )
|
|
|
+
|
|
|
+ prob.solve(solver=cp.ECOS, max_iters=300)
|
|
|
+ dp = dispatch_power.value
|
|
|
+
|
|
|
+ gap -= dp.sum(axis=0)
|
|
|
+ S = np.maximum(gap, 0).mean()
|
|
|
+
|
|
|
+ return S, power, gap, storage, dp
|
|
|
+
|
|
|
+def objective_with_storage_and_flexibility(
|
|
|
+ potential,
|
|
|
+ units_per_region,
|
|
|
+ dispatchable,
|
|
|
+ p_load=2664*1000/(365*24),
|
|
|
+ storage_capacities=5*250,
|
|
|
+ storage_max_loads=5*30,
|
|
|
+ storage_max_deliveries=5*30,
|
|
|
+ storage_efficiencies=0.67,
|
|
|
+ flexibility_power=17,
|
|
|
+ flexibility_time=8
|
|
|
+):
|
|
|
+
|
|
|
+ power = np.einsum('ijk,ik', potential, units_per_region)
|
|
|
+
|
|
|
+ T = len(power)
|
|
|
+ tau = flexibility_time
|
|
|
+
|
|
|
+ h = cp.Variable((T, tau+1))
|
|
|
+
|
|
|
+ constraints = [
|
|
|
+ h >= 0,
|
|
|
+ h <= 1,
|
|
|
+ h[:, 0] >= 1-flexibility_power/p_load,
|
|
|
+ cp.multiply(p_load, cp.sum(h, axis=1))-p_load <= flexibility_power,
|
|
|
+ ] + [
|
|
|
+ h[t, 0]+h[t-1, 1]+h[t-2, 2]+h[t-3, 3]+h[t-4, 4] +
|
|
|
+ h[t-5, 5]+h[t-6, 6]+h[t-7, 7]+h[t-8, 8] == 1
|
|
|
+ for t in np.arange(flexibility_time, T)
|
|
|
+ ]
|
|
|
+
|
|
|
+ prob = cp.Problem(
|
|
|
+ cp.Minimize(
|
|
|
+ cp.sum(cp.pos(cp.multiply(p_load, cp.sum(h, axis=1))-power))
|
|
|
+ ),
|
|
|
+ constraints
|
|
|
+ )
|
|
|
+
|
|
|
+ prob.solve(verbose=True, solver=cp.ECOS, max_iters=300)
|
|
|
+
|
|
|
+ hb = np.array(h.value)
|
|
|
+ p_load = p_load*np.sum(hb, axis=1)
|
|
|
+
|
|
|
+ power_delta = power-p_load
|
|
|
+ available_power = np.array(power_delta)
|
|
|
+ excess_power = np.maximum(0, available_power)
|
|
|
+ deficit_power = np.maximum(0, -available_power)
|
|
|
+
|
|
|
+ n_storages = len(storage_capacities)
|
|
|
+ storage_try_load = np.zeros((n_storages, T))
|
|
|
+ storage_try_delivery = np.zeros((n_storages, T))
|
|
|
+
|
|
|
+ storage = np.zeros((n_storages, T))
|
|
|
+ storage_impact = np.zeros((n_storages, T))
|
|
|
+ dE_storage = np.zeros((n_storages, T))
|
|
|
+
|
|
|
+ for i in range(n_storages):
|
|
|
+ storage_try_load[i] = np.minimum(excess_power, storage_max_loads[i])
|
|
|
+ storage_try_delivery[i] = np.minimum(deficit_power, storage_max_deliveries[i])
|
|
|
+
|
|
|
+ dE_storage[i] = storage_try_load[i]-storage_try_delivery[i]
|
|
|
+
|
|
|
+ storage[i] = storage_iterate(
|
|
|
+ dE_storage[i], storage_capacities[i], storage_efficiencies[i], T
|
|
|
+ )
|
|
|
+
|
|
|
+ # 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/storage_efficiencies[i], 1)
|
|
|
+ )
|
|
|
+
|
|
|
+ available_power += storage_impact[i]
|
|
|
+ excess_power = np.maximum(0, available_power)
|
|
|
+ deficit_power = np.maximum(0, -available_power)
|
|
|
+
|
|
|
+ # power missing to meet demand
|
|
|
+ gap = p_load-power-storage_impact.sum(axis=0)
|
|
|
+
|
|
|
+ # optimize dispatch
|
|
|
+ n_dispatchable_sources = dispatchable.shape[0]
|
|
|
+ dispatch_power = cp.Variable((n_dispatchable_sources, T))
|
|
|
+
|
|
|
+ constraints = [
|
|
|
+ dispatch_power >= 0,
|
|
|
+ cp.sum(dispatch_power, axis=0) <= np.maximum(gap, 0)
|
|
|
+ ] + [
|
|
|
+ dispatch_power[i,:] <= dispatchable[i,:,0].sum()
|
|
|
+ for i in range(n_dispatchable_sources)
|
|
|
+ ] + [
|
|
|
+ cp.sum(dispatch_power[i]) <= dispatchable[i,:,1].sum()
|
|
|
+ for i in range(n_dispatchable_sources)
|
|
|
+ ]
|
|
|
+
|
|
|
+ prob = cp.Problem(
|
|
|
+ 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
|
|
|
+
|
|
|
+ gap -= dp.sum(axis=0)
|
|
|
+ S = np.maximum(gap, 0).mean()
|
|
|
+
|
|
|
+ return S, power, gap, storage, p_load, dp
|
|
|
+
|
|
|
+
|
|
|
+potential = pd.read_parquet("potential.parquet")
|
|
|
+potential = potential.loc['1985-01-01 00:00:00':'2015-01-01 00:00:00', :]
|
|
|
+potential.fillna(0, inplace=True)
|
|
|
+
|
|
|
+potential = potential.loc[(
|
|
|
+ slice(f'{args.begin} 00:00:00', f'{args.end} 00:00:00'), 'FR'), :]
|
|
|
+
|
|
|
+load = generate_load_curve(potential.index.get_level_values(0))
|
|
|
+
|
|
|
+p = potential[["onshore", "offshore", "solar"]].to_xarray().to_array()
|
|
|
+p = np.insert(p, 3, 0.68, axis=0) # nuclear power-like
|
|
|
+# p[-1,:,0] = 0.68+0.25*np.cos(2*np.pi*np.arange(len(potential))/(365.25*24)) # nuclear seasonality ersatz
|
|
|
+
|
|
|
+dispatchable = np.zeros((3, 1, 2))
|
|
|
+
|
|
|
+# hydro power
|
|
|
+dispatchable[0][0][0] = 22
|
|
|
+dispatchable[0][0][1] = 63*1000*len(potential)/(365.25*24)
|
|
|
+
|
|
|
+# biomass
|
|
|
+dispatchable[1][0][0] = 2
|
|
|
+dispatchable[1][0][1] = 12*1000*len(potential)/(365.25*24)
|
|
|
+
|
|
|
+# thermique tradi
|
|
|
+dispatchable[2][0][0] = 0.5
|
|
|
+dispatchable[2][0][1] = 1e8
|
|
|
+
|
|
|
+n_sources = p.shape[0]
|
|
|
+n_dt = p.shape[1]
|
|
|
+n_regions = p.shape[2]
|
|
|
+
|
|
|
+# start = np.array([18, 0, 10, 61.4, 3]) # current mix
|
|
|
+# # stop = np.array([65, 40, 143, 0, 3]) # negawatt like mix
|
|
|
+# stop = np.array([74, 62, 208, 0, 3]) # RTE 100% EnR
|
|
|
+
|
|
|
+scenarios = {
|
|
|
+ 'M0': np.array([74, 62, 208, 0]),
|
|
|
+ 'M1': np.array([59, 45, 214, 16]),
|
|
|
+ 'M23': np.array([72, 60, 125, 16]),
|
|
|
+ 'N1': np.array([58, 45, 118, 16+13]),
|
|
|
+ 'N2': np.array([52, 36, 90, 16+23]),
|
|
|
+ 'N03': np.array([43, 22, 70, 24+27]),
|
|
|
+}
|
|
|
+
|
|
|
+BATTERY_CAPACITY=4 # 4 hour capacity
|
|
|
+STEP_CAPACITY=24
|
|
|
+P2G_CAPACITY=24*7*2
|
|
|
+
|
|
|
+battery_power = {
|
|
|
+ 'M0': 26,
|
|
|
+ 'M1': 21,
|
|
|
+ 'M23': 13,
|
|
|
+ 'N1': 9,
|
|
|
+ 'N2': 2,
|
|
|
+ 'N03': 1
|
|
|
+}
|
|
|
+
|
|
|
+p2g_power = {
|
|
|
+ 'M0': 29,
|
|
|
+ 'M1': 20,
|
|
|
+ 'M23': 20,
|
|
|
+ 'N1': 11,
|
|
|
+ 'N2': 5,
|
|
|
+ 'N03': 0
|
|
|
+}
|
|
|
+
|
|
|
+step_power = 8
|
|
|
+
|
|
|
+fig, axes = plt.subplots(nrows=6, ncols=2, sharex="col", sharey=True)
|
|
|
+w, h = fig.get_size_inches()
|
|
|
+fig.set_figwidth(w*1.5)
|
|
|
+fig.set_figheight(h*1.5)
|
|
|
+
|
|
|
+fig_storage, axes_storage = plt.subplots(nrows=6, ncols=2, sharex="col", sharey=True)
|
|
|
+fig_storage.set_figwidth(w*1.5)
|
|
|
+fig_storage.set_figheight(h*1.5)
|
|
|
+
|
|
|
+fig_dispatch, axes_dispatch = plt.subplots(nrows=6, ncols=2, sharex="col", sharey=True)
|
|
|
+fig_dispatch.set_figwidth(w*1.5)
|
|
|
+fig_dispatch.set_figheight(h*1.5)
|
|
|
+
|
|
|
+
|
|
|
+# for step in np.linspace(start, stop, 2050-2022, True)[::-1]:
|
|
|
+row = 0
|
|
|
+for scenario in scenarios:
|
|
|
+ units = np.transpose([scenarios[scenario]])
|
|
|
+
|
|
|
+ # 8 GW STEP capacity common to all scenarios
|
|
|
+ storage_power = battery_power[scenario] + p2g_power[scenario] + step_power
|
|
|
+ storage_max_load = battery_power[scenario]*BATTERY_CAPACITY + p2g_power[scenario]*24*7 + step_power*24
|
|
|
+
|
|
|
+ # nuclear_load_model(p[:,:3], units, load, 20)
|
|
|
+
|
|
|
+ if args.flexibility:
|
|
|
+ S, production, gap, storage, adjusted_load, dp = objective_with_storage_and_flexibility(
|
|
|
+ p, units,
|
|
|
+ dispatchable,
|
|
|
+ p_load=load,
|
|
|
+ storage_capacities=[battery_power[scenario]*BATTERY_CAPACITY, step_power*STEP_CAPACITY, p2g_power[scenario]*P2G_CAPACITY],
|
|
|
+ storage_max_deliveries=[battery_power[scenario], step_power, p2g_power[scenario]],
|
|
|
+ storage_max_loads=[battery_power[scenario], step_power, p2g_power[scenario]],
|
|
|
+ storage_efficiencies=[0.8, 0.8, 0.3]
|
|
|
+ )
|
|
|
+
|
|
|
+ print(f"{scenario}:", S, gap.max(), np.quantile(gap, 0.95))
|
|
|
+ else:
|
|
|
+ S, production, gap, storage, dp = objective_with_storages(
|
|
|
+ p,
|
|
|
+ units,
|
|
|
+ dispatchable,
|
|
|
+ p_load = load,
|
|
|
+ storage_capacities=[battery_power[scenario]*BATTERY_CAPACITY, step_power*STEP_CAPACITY, p2g_power[scenario]*P2G_CAPACITY],
|
|
|
+ storage_max_deliveries=[battery_power[scenario], step_power, p2g_power[scenario]],
|
|
|
+ storage_max_loads=[battery_power[scenario], step_power, p2g_power[scenario]],
|
|
|
+ storage_efficiencies=[0.8, 0.8, 0.3]
|
|
|
+ )
|
|
|
+ adjusted_load = load
|
|
|
+ print(f"{scenario} w/o flexibility:", S, gap.max(), np.quantile(gap, 0.95))
|
|
|
+
|
|
|
+ print(f"exports: {np.minimum(np.maximum(-gap, 0), 39).sum()/1000} TWh; imports: {np.minimum(np.maximum(gap, 0), 39).sum()/1000} TWh")
|
|
|
+ print(f"dispatchable: " + ", ".join([f"{dp[i].sum()/1000:.2f} TWh" for i in range(dp.shape[0])]))
|
|
|
+
|
|
|
+ potential['adjusted_load'] = adjusted_load
|
|
|
+ potential['production'] = production
|
|
|
+ potential['available'] = production-np.diff(storage.sum(axis=0), append=0)
|
|
|
+
|
|
|
+ for i in range(3):
|
|
|
+ potential[f"storage_{i}"] = np.diff(storage[i,:], append=0)#storage[i,:]/1000
|
|
|
+ potential[f"storage_{i}"] = storage[i,:]/1000
|
|
|
+
|
|
|
+ for i in range(dp.shape[0]):
|
|
|
+ potential[f"dispatch_{i}"] = dp[i,:]
|
|
|
+
|
|
|
+ potential["dispatch"] = dp.sum(axis=0)
|
|
|
+
|
|
|
+ data = [
|
|
|
+ potential.loc[(slice('2013-02-01 00:00:00',
|
|
|
+ '2013-03-01 00:00:00'), 'FR'), :],
|
|
|
+ potential.loc[(slice('2013-06-01 00:00:00',
|
|
|
+ '2013-07-01 00:00:00'), 'FR'), :]
|
|
|
+ ]
|
|
|
+
|
|
|
+ months = [
|
|
|
+ "Février",
|
|
|
+ "Juin"
|
|
|
+ ]
|
|
|
+
|
|
|
+ labels = [
|
|
|
+ "adjusted load (GW)",
|
|
|
+ "production (GW)",
|
|
|
+ "available power (production-storage) (GW)",
|
|
|
+ "power deficit"
|
|
|
+ ]
|
|
|
+
|
|
|
+ labels_storage = [
|
|
|
+ "Batterie (TWh)",
|
|
|
+ "STEP (TWh)",
|
|
|
+ "P2G (TWh)"
|
|
|
+ ]
|
|
|
+
|
|
|
+ labels_dispatch = [
|
|
|
+ "Hydraulique (GW)",
|
|
|
+ "Biomasse (GW)",
|
|
|
+ "Thermique gaz (GW)"
|
|
|
+ ]
|
|
|
+
|
|
|
+ for col in range(2):
|
|
|
+ ax = axes[row, col]
|
|
|
+ ax.plot(data[col].index.get_level_values(0), data[col]
|
|
|
+ ["adjusted_load"], label="adjusted load (GW)", lw=1)
|
|
|
+ ax.plot(data[col].index.get_level_values(0), data[col]
|
|
|
+ ["production"], label="production (GW)", ls="dotted", lw=1)
|
|
|
+ ax.plot(data[col].index.get_level_values(0), data[col]["available"],
|
|
|
+ label="available power (production-d(storage)/dt) (GW)", lw=1)
|
|
|
+
|
|
|
+ ax.fill_between(
|
|
|
+ data[col].index.get_level_values(0),
|
|
|
+ data[col]["available"],
|
|
|
+ data[col]["adjusted_load"],
|
|
|
+ where=data[col]["adjusted_load"] > data[col]["available"],
|
|
|
+ color='red',
|
|
|
+ alpha=0.15
|
|
|
+ )
|
|
|
+
|
|
|
+ fmt = mdates.DateFormatter('%d/%m')
|
|
|
+ ax.xaxis.set_major_formatter(fmt)
|
|
|
+ ax.text(
|
|
|
+ 0.5, 0.87, f"Scénario {scenario} ({months[col]})", ha='center', transform=ax.transAxes)
|
|
|
+
|
|
|
+ ax.set_ylim(25, 225)
|
|
|
+
|
|
|
+ ax = axes_storage[row, col]
|
|
|
+ for i in np.arange(3):
|
|
|
+ if i == 2:
|
|
|
+ base = 0
|
|
|
+ else:
|
|
|
+ base = np.sum([data[col][f"storage_{j}"] for j in np.arange(i+1,3)], axis=0)
|
|
|
+
|
|
|
+ ax.fill_between(data[col].index.get_level_values(0), base, base+data[col]
|
|
|
+ [f"storage_{i}"], label=f"storage {i}", alpha=0.5)
|
|
|
+
|
|
|
+ ax.plot(data[col].index.get_level_values(0), base+data[col]
|
|
|
+ [f"storage_{i}"], label=f"storage {i}", lw=0.25)
|
|
|
+
|
|
|
+ ax.xaxis.set_major_formatter(fmt)
|
|
|
+ ax.text(
|
|
|
+ 0.5, 0.87, f"Scénario {scenario} ({months[col]})", ha='center', transform=ax.transAxes)
|
|
|
+
|
|
|
+ ax = axes_dispatch[row, col]
|
|
|
+ for i in range(dp.shape[0]):
|
|
|
+ if i == 0:
|
|
|
+ base = 0
|
|
|
+ else:
|
|
|
+ base = np.sum([data[col][f"dispatch_{j}"] for j in np.arange(i)], axis=0)
|
|
|
+
|
|
|
+ ax.fill_between(data[col].index.get_level_values(0), base, base+data[col]
|
|
|
+ [f"dispatch_{i}"], label=f"dispatch {i}", alpha=0.5)
|
|
|
+
|
|
|
+ ax.plot(data[col].index.get_level_values(0), base+data[col]
|
|
|
+ [f"dispatch_{i}"], label=f"dispatch {i}", lw=0.25)
|
|
|
+
|
|
|
+ ax.xaxis.set_major_formatter(fmt)
|
|
|
+ ax.text(
|
|
|
+ 0.5, 0.87, f"Scénario {scenario} ({months[col]})", ha='center', transform=ax.transAxes)
|
|
|
+
|
|
|
+ row += 1
|
|
|
+
|
|
|
+for label in axes[-1, 0].get_xmajorticklabels() + axes[-1, 1].get_xmajorticklabels():
|
|
|
+ label.set_rotation(30)
|
|
|
+ label.set_horizontalalignment("right")
|
|
|
+
|
|
|
+for label in axes_storage[-1, 0].get_xmajorticklabels() + axes_storage[-1, 1].get_xmajorticklabels():
|
|
|
+ label.set_rotation(30)
|
|
|
+ label.set_horizontalalignment("right")
|
|
|
+
|
|
|
+for label in axes_dispatch[-1, 0].get_xmajorticklabels() + axes_dispatch[-1, 1].get_xmajorticklabels():
|
|
|
+ label.set_rotation(30)
|
|
|
+ label.set_horizontalalignment("right")
|
|
|
+
|
|
|
+flex = "With" if args.flexibility else "Without"
|
|
|
+
|
|
|
+plt.subplots_adjust(wspace=0, hspace=0)
|
|
|
+fig.suptitle(f"Simulations based on {args.begin}--{args.end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)")
|
|
|
+fig.text(1, 0, 'Lucas Gautheron', ha="right")
|
|
|
+fig.legend(labels, loc='lower right', bbox_to_anchor=(1, -0.1),
|
|
|
+ ncol=len(labels), bbox_transform=fig.transFigure)
|
|
|
+fig.savefig("output.png", bbox_inches="tight", dpi=200)
|
|
|
+
|
|
|
+fig_storage.suptitle(f"Simulations based on {args.begin}--{args.end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)")
|
|
|
+fig_storage.text(1, 0, 'Lucas Gautheron', ha="right")
|
|
|
+fig_storage.legend(labels_storage, loc='lower right', bbox_to_anchor=(1, -0.1),
|
|
|
+ ncol=len(labels_storage), bbox_transform=fig_storage.transFigure)
|
|
|
+fig_storage.savefig("output_storage.png", bbox_inches="tight", dpi=200)
|
|
|
+
|
|
|
+fig_dispatch.suptitle(f"Simulations based on {args.begin}--{args.end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)")
|
|
|
+fig_dispatch.text(1, 0, 'Lucas Gautheron', ha="right")
|
|
|
+fig_dispatch.legend(labels_dispatch, loc='lower right', bbox_to_anchor=(1, -0.1),
|
|
|
+ ncol=len(labels_dispatch), bbox_transform=fig_dispatch.transFigure)
|
|
|
+fig_dispatch.savefig("output_dispatch.png", bbox_inches="tight", dpi=200)
|
|
|
+plt.show()
|
