run.py 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338
  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(
  10. "--begin",
  11. help="begin date (YYYY-MM-DD), between 1985-01-01 and 2015-01-01",
  12. required=True,
  13. )
  14. parser.add_argument(
  15. "--end",
  16. help="end date (YYYY-MM-DD), between 1985-01-01 and 2015-01-01",
  17. required=True,
  18. )
  19. parser.add_argument(
  20. "--flexibility", help="enable load flexibility modeling", action="store_true"
  21. )
  22. parser.add_argument(
  23. "--scenarios",
  24. help="path to scenarios parameters yml file",
  25. default="scenarios/rte_2050.yml",
  26. )
  27. args = parser.parse_args()
  28. with open(args.scenarios, "r") as f:
  29. scenarios = yaml.load(f, Loader=yaml.FullLoader)
  30. potential = pd.read_parquet("data/potential.parquet")
  31. potential = potential.loc["1985-01-01 00:00:00":"2015-01-01 00:00:00", :]
  32. potential.fillna(0, inplace=True)
  33. begin = args.begin
  34. end = args.end
  35. flexibility = args.flexibility
  36. potential = potential.loc[(slice(f"{begin} 00:00:00", f"{end} 00:00:00"), "FR"), :]
  37. # intermittent sources potential
  38. p = potential[["onshore", "offshore", "solar"]].to_xarray().to_array()
  39. p = np.insert(p, 3, 0.7, axis=0) # nuclear power-like
  40. times = potential.index.get_level_values(0)
  41. n_scenarios = len(scenarios)
  42. fig, axes = plt.subplots(nrows=n_scenarios, ncols=2, sharex="col", sharey=True)
  43. w, h = fig.get_size_inches()
  44. fig.set_figwidth(w * 1.5)
  45. fig.set_figheight(h * 1.5)
  46. fig_storage, axes_storage = plt.subplots(
  47. nrows=n_scenarios, ncols=2, sharex="col", sharey=True
  48. )
  49. fig_storage.set_figwidth(w * 1.5)
  50. fig_storage.set_figheight(h * 1.5)
  51. fig_dispatch, axes_dispatch = plt.subplots(
  52. nrows=n_scenarios, ncols=2, sharex="col", sharey=True
  53. )
  54. fig_dispatch.set_figwidth(w * 1.5)
  55. fig_dispatch.set_figheight(h * 1.5)
  56. fig_gap_distribution, axes_gap_distribution = plt.subplots(
  57. nrows=int(np.ceil(n_scenarios / 2)), ncols=2, sharex="col", sharey=True
  58. )
  59. fig_gap_distribution.set_figwidth(w * 1.5)
  60. fig_gap_distribution.set_figheight(h * 1.5)
  61. months = ["Février", "Juin"]
  62. labels = [
  63. "load (GW)",
  64. "production (GW)",
  65. "available power (production-storage) (GW)",
  66. "power deficit",
  67. ]
  68. labels_storage = ["Batteries (TWh)", "STEP (TWh)", "P2G (TWh)"]
  69. labels_dispatch = ["Hydro (GW)", "Biomass (GW)", "Thermal (GW)"]
  70. date_fmt = mdates.DateFormatter("%d/%m")
  71. # for step in np.linspace(start, stop, 2050-2022, True)[::-1]:
  72. row = 0
  73. for scenario in scenarios:
  74. if not flexibility:
  75. scenarios[scenario]["flexibility_power"] = 0
  76. scenario_model = Scenario(**scenarios[scenario])
  77. S, load, production, gap, storage, dp = scenario_model.run(times, p)
  78. print(f"{scenario}:", S, gap.max(), np.quantile(gap, 0.95))
  79. print(
  80. f"exports: {np.minimum(np.maximum(-gap, 0), 39).sum()/1000} TWh; imports: {np.minimum(np.maximum(gap, 0), 39).sum()/1000} TWh"
  81. )
  82. print(
  83. f"dispatchable: "
  84. + ", ".join([f"{dp[i].sum()/1000:.2f} TWh" for i in range(dp.shape[0])])
  85. )
  86. potential["load"] = load
  87. potential["production"] = production
  88. potential["available"] = production - np.diff(storage.sum(axis=0), append=0)
  89. potential["gap"] = gap
  90. for i in range(3):
  91. potential[f"storage_{i}"] = np.diff(
  92. storage[i, :], append=0
  93. ) # storage[i,:]/1000
  94. potential[f"storage_{i}"] = storage[i, :] / 1000
  95. for i in range(dp.shape[0]):
  96. potential[f"dispatch_{i}"] = dp[i, :]
  97. potential["dispatch"] = dp.sum(axis=0)
  98. data = [
  99. potential.loc[(slice("2013-02-01 00:00:00", "2013-03-01 00:00:00"), "FR"), :],
  100. potential.loc[(slice("2013-06-01 00:00:00", "2013-07-01 00:00:00"), "FR"), :],
  101. ]
  102. for col in range(2):
  103. ax = axes[row, col] if axes.ndim > 1 else axes[col]
  104. ax.plot(
  105. data[col].index.get_level_values(0),
  106. data[col]["load"],
  107. label="adjusted load (GW)",
  108. lw=1,
  109. )
  110. ax.plot(
  111. data[col].index.get_level_values(0),
  112. data[col]["production"],
  113. label="production (GW)",
  114. ls="dotted",
  115. lw=1,
  116. )
  117. ax.plot(
  118. data[col].index.get_level_values(0),
  119. data[col]["available"],
  120. label="available power (production-d(storage)/dt) (GW)",
  121. lw=1,
  122. )
  123. ax.fill_between(
  124. data[col].index.get_level_values(0),
  125. data[col]["available"],
  126. data[col]["load"],
  127. where=data[col]["load"] > data[col]["available"],
  128. color="red",
  129. alpha=0.15,
  130. )
  131. ax.xaxis.set_major_formatter(date_fmt)
  132. ax.text(
  133. 0.5,
  134. 0.87,
  135. f"Scénario {scenario} ({months[col]})",
  136. ha="center",
  137. transform=ax.transAxes,
  138. )
  139. ax.set_ylim(10, 210)
  140. ax = axes_storage[row, col] if axes.ndim > 1 else axes_storage[col]
  141. for i in np.arange(3):
  142. if i == 2:
  143. base = 0
  144. else:
  145. base = np.sum(
  146. [data[col][f"storage_{j}"] for j in np.arange(i + 1, 3)], axis=0
  147. )
  148. ax.fill_between(
  149. data[col].index.get_level_values(0),
  150. base,
  151. base + data[col][f"storage_{i}"],
  152. label=f"storage {i}",
  153. alpha=0.5,
  154. )
  155. ax.plot(
  156. data[col].index.get_level_values(0),
  157. base + data[col][f"storage_{i}"],
  158. label=f"storage {i}",
  159. lw=0.25,
  160. )
  161. ax.xaxis.set_major_formatter(date_fmt)
  162. ax.text(
  163. 0.5,
  164. 0.87,
  165. f"Scénario {scenario} ({months[col]})",
  166. ha="center",
  167. transform=ax.transAxes,
  168. )
  169. ax = axes_dispatch[row, col] if axes.ndim > 1 else axes_dispatch[col]
  170. for i in range(dp.shape[0]):
  171. if i == 0:
  172. base = 0
  173. else:
  174. base = np.sum(
  175. [data[col][f"dispatch_{j}"] for j in np.arange(i)], axis=0
  176. )
  177. ax.fill_between(
  178. data[col].index.get_level_values(0),
  179. base,
  180. base + data[col][f"dispatch_{i}"],
  181. label=f"dispatch {i}",
  182. alpha=0.5,
  183. )
  184. ax.plot(
  185. data[col].index.get_level_values(0),
  186. base + data[col][f"dispatch_{i}"],
  187. label=f"dispatch {i}",
  188. lw=0.25,
  189. )
  190. ax.xaxis.set_major_formatter(date_fmt)
  191. ax.text(
  192. 0.5,
  193. 0.87,
  194. f"Scénario {scenario} ({months[col]})",
  195. ha="center",
  196. transform=ax.transAxes,
  197. )
  198. div = int(np.ceil(n_scenarios / 2))
  199. ax = (
  200. axes_gap_distribution[row % div, row // div]
  201. if axes_gap_distribution.ndim > 1
  202. else axes_gap_distribution[row % div]
  203. )
  204. hist, bin_edges = np.histogram(-gap, bins=1000)
  205. hist = np.cumsum(hist)
  206. hist = 100 * (hist - hist.min()) / hist.ptp()
  207. keep = np.abs(bin_edges[:-1]) < 50
  208. ax.plot(bin_edges[:-1][keep], hist[keep], lw=1, label="power gap", color="#ff7f00")
  209. years = pd.date_range(start=begin, end=end, freq="Y")
  210. for i in range(len(years) - 1):
  211. year_data = potential.loc[
  212. (slice(f"{years[i]} 00:00:00", f"{years[i+1]} 00:00:00"), "FR"), "gap"
  213. ]
  214. hist, bin_edges = np.histogram(-year_data, bins=1000)
  215. hist = np.cumsum(hist)
  216. hist = 100 * (hist - hist.min()) / hist.ptp()
  217. keep = np.abs(bin_edges[:-1]) < 50
  218. ax.plot(
  219. bin_edges[:-1][keep], hist[keep], lw=0.5, alpha=0.2, color="#ff7f00"
  220. )
  221. ax.text(0.5, 0.87, f"Scénario {scenario}", ha="center", transform=ax.transAxes)
  222. row += 1
  223. for axs in [axes, axes_dispatch, axes_storage]:
  224. if axes.ndim > 1:
  225. for label in (
  226. axs[-1, 0].get_xmajorticklabels() + axs[-1, 1].get_xmajorticklabels()
  227. ):
  228. label.set_rotation(30)
  229. label.set_horizontalalignment("right")
  230. else:
  231. for label in axs[0].get_xmajorticklabels() + axs[-1].get_xmajorticklabels():
  232. label.set_rotation(30)
  233. label.set_horizontalalignment("right")
  234. flex = "With" if flexibility else "Without"
  235. plt.subplots_adjust(wspace=0, hspace=0)
  236. fig.suptitle(
  237. f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)"
  238. )
  239. fig.text(1, 0, "Lucas Gautheron", ha="right")
  240. fig.legend(
  241. labels,
  242. loc="lower right",
  243. bbox_to_anchor=(1, -0.1),
  244. ncol=len(labels),
  245. bbox_transform=fig.transFigure,
  246. )
  247. fig.savefig("output/load_supply.png", bbox_inches="tight", dpi=200)
  248. fig_storage.suptitle(
  249. f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)"
  250. )
  251. fig_storage.text(1, 0, "Lucas Gautheron", ha="right")
  252. fig_storage.legend(
  253. labels_storage,
  254. loc="lower right",
  255. bbox_to_anchor=(1, -0.1),
  256. ncol=len(labels_storage),
  257. bbox_transform=fig_storage.transFigure,
  258. )
  259. fig_storage.savefig("output/storage.png", bbox_inches="tight", dpi=200)
  260. fig_dispatch.suptitle(
  261. f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)"
  262. )
  263. fig_dispatch.text(1, 0, "Lucas Gautheron", ha="right")
  264. fig_dispatch.legend(
  265. labels_dispatch,
  266. loc="lower right",
  267. bbox_to_anchor=(1, -0.1),
  268. ncol=len(labels_dispatch),
  269. bbox_transform=fig_dispatch.transFigure,
  270. )
  271. fig_dispatch.savefig("output/dispatch.png", bbox_inches="tight", dpi=200)
  272. fig_gap_distribution.suptitle(
  273. f"Power gap cumulative distribution (%)\nSimulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)"
  274. )
  275. fig_gap_distribution.legend(
  276. ["Power gap (available-load) (GW)"],
  277. loc="lower right",
  278. bbox_to_anchor=(1, -0.1),
  279. ncol=1,
  280. bbox_transform=fig_dispatch.transFigure,
  281. )
  282. fig_gap_distribution.text(1, 0, "Lucas Gautheron", ha="right")
  283. fig_gap_distribution.savefig(
  284. "output/gap_distribution.png", bbox_inches="tight", dpi=200
  285. )
  286. plt.show()