run.py 8.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225
  1. import pandas as pd
  2. import numpy as np
  3. from matplotlib import pyplot as plt
  4. import matplotlib.dates as mdates
  5. import argparse
  6. import yaml
  7. from mix_simul.scenarios import Scenario
  8. parser = argparse.ArgumentParser()
  9. parser.add_argument("--begin", help="begin date (YYYY-MM-DD), between 1985-01-01 and 2015-01-01", required=True)
  10. parser.add_argument("--end", help="end date (YYYY-MM-DD), between 1985-01-01 and 2015-01-01", required=True)
  11. parser.add_argument("--flexibility", help="enable load flexibility modeling", action="store_true")
  12. parser.add_argument("--scenarios", help="path to scenarios parameters yml file", default="scenarios/rte.yml")
  13. args = parser.parse_args()
  14. with open(args.scenarios, 'r') as f:
  15. scenarios = yaml.load(f, Loader=yaml.FullLoader)
  16. potential = pd.read_parquet("data/potential.parquet")
  17. potential = potential.loc['1985-01-01 00:00:00':'2015-01-01 00:00:00', :]
  18. potential.fillna(0, inplace=True)
  19. begin = args.begin
  20. end = args.end
  21. flexibility = args.flexibility
  22. potential = potential.loc[(
  23. slice(f'{begin} 00:00:00', f'{end} 00:00:00'), 'FR'), :]
  24. # intermittent sources potential
  25. p = potential[["onshore", "offshore", "solar"]].to_xarray().to_array()
  26. p = np.insert(p, 3, 0.7, axis=0) # nuclear power-like
  27. times = potential.index.get_level_values(0)
  28. n_scenarios = len(scenarios)
  29. fig, axes = plt.subplots(nrows=n_scenarios, ncols=2, sharex="col", sharey=True)
  30. w, h = fig.get_size_inches()
  31. fig.set_figwidth(w*1.5)
  32. fig.set_figheight(h*1.5)
  33. fig_storage, axes_storage = plt.subplots(nrows=n_scenarios, ncols=2, sharex="col", sharey=True)
  34. fig_storage.set_figwidth(w*1.5)
  35. fig_storage.set_figheight(h*1.5)
  36. fig_dispatch, axes_dispatch = plt.subplots(nrows=n_scenarios, ncols=2, sharex="col", sharey=True)
  37. fig_dispatch.set_figwidth(w*1.5)
  38. fig_dispatch.set_figheight(h*1.5)
  39. fig_gap_distribution, axes_gap_distribution = plt.subplots(nrows=int(np.ceil(n_scenarios/2)), ncols=2, sharex="col", sharey=True)
  40. fig_gap_distribution.set_figwidth(w*1.5)
  41. fig_gap_distribution.set_figheight(h*1.5)
  42. months = [
  43. "Février",
  44. "Juin"
  45. ]
  46. labels = [
  47. "load (GW)",
  48. "production (GW)",
  49. "available power (production-storage) (GW)",
  50. "power deficit"
  51. ]
  52. labels_storage = [
  53. "Batteries (TWh)",
  54. "STEP (TWh)",
  55. "P2G (TWh)"
  56. ]
  57. labels_dispatch = [
  58. "Hydro (GW)",
  59. "Biomass (GW)",
  60. "Thermal (GW)"
  61. ]
  62. date_fmt = mdates.DateFormatter('%d/%m')
  63. # for step in np.linspace(start, stop, 2050-2022, True)[::-1]:
  64. row = 0
  65. for scenario in scenarios:
  66. if not flexibility:
  67. scenarios[scenario]["flexibility_power"] = 0
  68. scenario_model = Scenario(**scenarios[scenario])
  69. S, load, production, gap, storage, dp = scenario_model.run(times, p)
  70. print(f"{scenario}:", S, gap.max(), np.quantile(gap, 0.95))
  71. print(f"exports: {np.minimum(np.maximum(-gap, 0), 39).sum()/1000} TWh; imports: {np.minimum(np.maximum(gap, 0), 39).sum()/1000} TWh")
  72. print(f"dispatchable: " + ", ".join([f"{dp[i].sum()/1000:.2f} TWh" for i in range(dp.shape[0])]))
  73. potential['load'] = load
  74. potential['production'] = production
  75. potential['available'] = production-np.diff(storage.sum(axis=0), append=0)
  76. for i in range(3):
  77. potential[f"storage_{i}"] = np.diff(storage[i,:], append=0)#storage[i,:]/1000
  78. potential[f"storage_{i}"] = storage[i,:]/1000
  79. for i in range(dp.shape[0]):
  80. potential[f"dispatch_{i}"] = dp[i,:]
  81. potential["dispatch"] = dp.sum(axis=0)
  82. data = [
  83. potential.loc[(slice('2013-02-01 00:00:00',
  84. '2013-03-01 00:00:00'), 'FR'), :],
  85. potential.loc[(slice('2013-06-01 00:00:00',
  86. '2013-07-01 00:00:00'), 'FR'), :]
  87. ]
  88. for col in range(2):
  89. ax = axes[row, col] if axes.ndim > 1 else axes[col]
  90. ax.plot(data[col].index.get_level_values(0), data[col]
  91. ["load"], label="adjusted load (GW)", lw=1)
  92. ax.plot(data[col].index.get_level_values(0), data[col]
  93. ["production"], label="production (GW)", ls="dotted", lw=1)
  94. ax.plot(data[col].index.get_level_values(0), data[col]["available"],
  95. label="available power (production-d(storage)/dt) (GW)", lw=1)
  96. ax.fill_between(
  97. data[col].index.get_level_values(0),
  98. data[col]["available"],
  99. data[col]["load"],
  100. where=data[col]["load"] > data[col]["available"],
  101. color='red',
  102. alpha=0.15
  103. )
  104. ax.xaxis.set_major_formatter(date_fmt)
  105. ax.text(
  106. 0.5, 0.87, f"Scénario {scenario} ({months[col]})", ha='center', transform=ax.transAxes)
  107. ax.set_ylim(10, 210)
  108. ax = axes_storage[row, col] if axes.ndim > 1 else axes_storage[col]
  109. for i in np.arange(3):
  110. if i == 2:
  111. base = 0
  112. else:
  113. base = np.sum([data[col][f"storage_{j}"] for j in np.arange(i+1,3)], axis=0)
  114. ax.fill_between(data[col].index.get_level_values(0), base, base+data[col]
  115. [f"storage_{i}"], label=f"storage {i}", alpha=0.5)
  116. ax.plot(data[col].index.get_level_values(0), base+data[col]
  117. [f"storage_{i}"], label=f"storage {i}", lw=0.25)
  118. ax.xaxis.set_major_formatter(date_fmt)
  119. ax.text(
  120. 0.5, 0.87, f"Scénario {scenario} ({months[col]})", ha='center', transform=ax.transAxes)
  121. ax = axes_dispatch[row, col] if axes.ndim > 1 else axes_dispatch[col]
  122. for i in range(dp.shape[0]):
  123. if i == 0:
  124. base = 0
  125. else:
  126. base = np.sum([data[col][f"dispatch_{j}"] for j in np.arange(i)], axis=0)
  127. ax.fill_between(data[col].index.get_level_values(0), base, base+data[col]
  128. [f"dispatch_{i}"], label=f"dispatch {i}", alpha=0.5)
  129. ax.plot(data[col].index.get_level_values(0), base+data[col]
  130. [f"dispatch_{i}"], label=f"dispatch {i}", lw=0.25)
  131. ax.xaxis.set_major_formatter(date_fmt)
  132. ax.text(
  133. 0.5, 0.87, f"Scénario {scenario} ({months[col]})", ha='center', transform=ax.transAxes)
  134. div = int(np.ceil(n_scenarios/2))
  135. ax = axes_gap_distribution[row%div, row//div] if axes_gap_distribution.ndim > 1 else axes_gap_distribution[row%div]
  136. hist, bin_edges = np.histogram(-gap, bins=1000)
  137. hist = np.cumsum(hist)
  138. hist = 100*(hist - hist.min()) / hist.ptp()
  139. keep = np.abs(bin_edges[:-1]) < 50
  140. ax.plot(bin_edges[:-1][keep], hist[keep], label="power gap")
  141. ax.text(
  142. 0.5, 0.87, f"Scénario {scenario}", ha='center', transform=ax.transAxes)
  143. row += 1
  144. for axs in [axes, axes_dispatch, axes_storage]:
  145. if axes.ndim > 1:
  146. for label in axs[-1, 0].get_xmajorticklabels() + axs[-1, 1].get_xmajorticklabels():
  147. label.set_rotation(30)
  148. label.set_horizontalalignment("right")
  149. else:
  150. for label in axs[0].get_xmajorticklabels() + axs[-1].get_xmajorticklabels():
  151. label.set_rotation(30)
  152. label.set_horizontalalignment("right")
  153. flex = "With" if flexibility else "Without"
  154. plt.subplots_adjust(wspace=0, hspace=0)
  155. fig.suptitle(f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)")
  156. fig.text(1, 0, 'Lucas Gautheron', ha="right")
  157. fig.legend(labels, loc='lower right', bbox_to_anchor=(1, -0.1),
  158. ncol=len(labels), bbox_transform=fig.transFigure)
  159. fig.savefig("output/load_supply.png", bbox_inches="tight", dpi=200)
  160. fig_storage.suptitle(f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)")
  161. fig_storage.text(1, 0, 'Lucas Gautheron', ha="right")
  162. fig_storage.legend(labels_storage, loc='lower right', bbox_to_anchor=(1, -0.1),
  163. ncol=len(labels_storage), bbox_transform=fig_storage.transFigure)
  164. fig_storage.savefig("output/storage.png", bbox_inches="tight", dpi=200)
  165. fig_dispatch.suptitle(f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)")
  166. fig_dispatch.text(1, 0, 'Lucas Gautheron', ha="right")
  167. fig_dispatch.legend(labels_dispatch, loc='lower right', bbox_to_anchor=(1, -0.1),
  168. ncol=len(labels_dispatch), bbox_transform=fig_dispatch.transFigure)
  169. fig_dispatch.savefig("output/dispatch.png", bbox_inches="tight", dpi=200)
  170. fig_gap_distribution.suptitle(f"Power gap cumulative distribution (%)\nSimulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)")
  171. fig_gap_distribution.legend(["Power gap (available-load) (GW)"], loc='lower right', bbox_to_anchor=(1, -0.1),
  172. ncol=1, bbox_transform=fig_dispatch.transFigure)
  173. fig_gap_distribution.text(1, 0, 'Lucas Gautheron', ha="right")
  174. fig_gap_distribution.savefig("output/gap_distribution.png", bbox_inches="tight", dpi=200)
  175. plt.show()