{ "cells": [ { "cell_type": "markdown", "id": "50f2d44d", "metadata": {}, "source": [ "\n", "****\n", "# **Code for submitted paper**:\n", "## **The economic value of load-shifting and its viability as a substitute for grid expansion**\n", "\n", "### Doing the Modelling and Forecasting\n", " - Forecasting\n", " - Calculation of incentives\n", " - Plots of forecasting results\n", " - Calculation of incentives under the perfect fit (Forecast = Actual)\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "eef9be35", "metadata": {}, "outputs": [], "source": [ "%reset\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib as mpl\n", "from sklearn.metrics import mean_pinball_loss\n", "import plotly.graph_objs as go\n", "from plotly.subplots import make_subplots\n", "import plotly.graph_objects as go\n", "from sklearn.ensemble import HistGradientBoostingRegressor\n", "from sklearn.model_selection import RandomizedSearchCV, TimeSeriesSplit\n", "from sklearn.metrics import make_scorer, mean_pinball_loss\n", "from pathlib import Path\n", "\n", "plt.style.use(\"seaborn-v0_8-whitegrid\")\n", "mpl.rcParams['axes.linewidth'] = 1\n", "np.random.seed(42)" ] }, { "cell_type": "code", "execution_count": 12, "id": "d083851c", "metadata": {}, "outputs": [], "source": [ "# ======================================================================\n", "# Reading-In the created dataframes\n", "# ======================================================================\n", "\n", "DATA_DIR = next(\n", " p for p in [\n", " Path(\"Data\"),\n", " Path(\"../Data\"),\n", " Path(\"ReadHandIn/Data\"),\n", " ]\n", " if p.is_dir()\n", ")\n", "\n", "df_EV20_PV25_5 = pd.read_csv(DATA_DIR / \"df_EV20_PV25_5.csv\", parse_dates=[\"clock_local\"])\n", "df_EV40_PV41_5 = pd.read_csv(DATA_DIR / \"df_EV40_PV41_5.csv\", parse_dates=[\"clock_local\"])\n", "df_EV75_PV75 = pd.read_csv(DATA_DIR / \"df_EV75_PV75.csv\", parse_dates=[\"clock_local\"])\n", "\n", "df_EV20_PV25_5['clock_local'] = pd.to_datetime(df_EV20_PV25_5['clock_local'])\n", "df_EV40_PV41_5['clock_local'] = pd.to_datetime(df_EV40_PV41_5['clock_local'])\n", "df_EV75_PV75['clock_local'] = pd.to_datetime(df_EV75_PV75['clock_local'])" ] }, { "cell_type": "code", "execution_count": 13, "id": "3efa14d1", "metadata": {}, "outputs": [], "source": [ "# ======================================================================\n", "# Function to get economic incentives of load-shifting (Following the approach from Appendix D)\n", "# ======================================================================\n", "\n", "def make_incentives (forecast_results, scenario):\n", "\n", " res = forecast_results[['clock_local','Adjusted_constumption_to_model','prediction']]\n", " res = res.rename(columns={\"Adjusted_constumption_to_model\": \"ytrue_cont\", \"prediction\": \"ypred_cont\"})\n", "\n", "\n", " # Thresholds --------------------------------------------------------------------------------------------------------------------------\n", " q85 = res[res['ytrue_cont'] > 0]['ytrue_cont'].quantile(0.85)\n", " q95 = res[res['ytrue_cont'] > 0]['ytrue_cont'].quantile(0.95)\n", " B0 = 0\n", "\n", " # \"Predicted Exceedances to shift\" regardless of true/false positive ------------------------------------------------------------------\n", " q85_shift_kw = (res[res['ypred_cont'] > q85]['ypred_cont'] - q85).sum()\n", " q95_shift_kw = (res[res['ypred_cont'] > q95]['ypred_cont'] - q95).sum()\n", " B0_shift_kw = (res[res['ypred_cont'] < B0]['ypred_cont'] - B0).sum() * (-1)\n", "\n", " #Saving kW--------------------------------------------------------------------------------------------------------------------------------\n", " missed_q85 = (res['ytrue_cont'].max()) - (res[res['ypred_cont'] < q85]['ytrue_cont'].max())\n", " missed_q95 = (res['ytrue_cont'].max()) - (res[res['ypred_cont'] < q95]['ytrue_cont'].max())\n", " missed_B0 = ((res['ytrue_cont'].min()) - (res[res['ypred_cont'] > B0]['ytrue_cont'].min())) * (-1)\n", "\n", " underp_q85 = (res['ytrue_cont'].max()) - (((res[res['ypred_cont'] >= q85]['ytrue_cont']) - res[res['ypred_cont'] >= q85]['ypred_cont']).max() + q85)\n", " underp_q95 = (res['ytrue_cont'].max()) - (((res[res['ypred_cont'] >= q95]['ytrue_cont']) - res[res['ypred_cont'] >= q95]['ypred_cont']).max() + q95)\n", " underp_B0 = ((res['ytrue_cont'].min()) - (((res[res['ypred_cont'] <= B0]['ytrue_cont']) - res[res['ypred_cont'] <= B0]['ypred_cont']).min())) * (-1)\n", "\n", " Sav_CAP_q85 = min(missed_q85, underp_q85)\n", " Sav_CAP_q95 = min(missed_q95, underp_q95)\n", " Sav_CAP_B0 = min(missed_B0, underp_B0)\n", " print(\"Underp\", underp_B0)\n", " print(\"Missed\", missed_B0)\n", " print(\"SAc_CAP\",Sav_CAP_B0)\n", " \n", " #grid expansion costs €/kw--------------------------------------------------------------------------------------------------------------------------------\n", " grid_exp_costs_ait_2030 = 442.1769\n", " grid_exp_costs_ait_2040 = 384.9372\n", "\n", " grid_exp_cost = 0\n", " if scenario == \"low\":\n", " grid_exp_cost = grid_exp_costs_ait_2030\n", " elif scenario == \"medium\":\n", " grid_exp_cost = grid_exp_costs_ait_2040\n", " elif scenario == \"high\":\n", " grid_exp_cost = grid_exp_costs_ait_2040\n", "\n", " #budget----------------------------------------------------------------------------------------------------------------------------------------------------\n", " budget_q85 = Sav_CAP_q85 * grid_exp_cost\n", " budget_q95 = Sav_CAP_q95 * grid_exp_cost\n", " budget_B0 = Sav_CAP_B0 * grid_exp_cost\n", "\n", " # Results -------------------------------------------------------------------------------------------------------------------------------------------\n", " print(\"Old Max\", res['ytrue_cont'].max())\n", " print(\"Old Min\", res['ytrue_cont'].min())\n", " print(\"\")\n", " print(\"New Max q85\", res['ytrue_cont'].max() - Sav_CAP_q85)\n", " print(\"New Max q95\", res['ytrue_cont'].max() - Sav_CAP_q95)\n", " print(\"New Min B0\", res['ytrue_cont'].min() + Sav_CAP_B0)\n", " print(\"\")\n", " print(\"Sum KW to shift q85\", q85_shift_kw)\n", " print(\"Sum KW to shift q95\", q95_shift_kw)\n", " print(\"Sum KW to shift B0\", B0_shift_kw)\n", " print(\"\")\n", " print(\"Sum Budget from Cap Saving q85\", budget_q85)\n", " print(\"Sum Budget from Cap Saving q95\", budget_q95)\n", " print(\"Sum Budget from Cap Saving B0\", budget_B0)\n", " print(\"\")\n", "\n", " costs = budget_B0 + budget_q85\n", " annual_budget_bo_q85 = costs * ((0.0611)/(1-(1+0.0611)**(-20)))\n", " total_exc = B0_shift_kw + q85_shift_kw\n", " annual_kw_incentive_bo_q85 = annual_budget_bo_q85 / total_exc\n", " print('ALL (q85 + B0) Incentive per household € per year: ', (annual_budget_bo_q85 / 1292))\n", " print('ALL (q85 + B0) Incentive per €/kw: ', annual_kw_incentive_bo_q85)\n", " print(\"annual_budget_bo_q85\", annual_budget_bo_q85)\n", "\n", " costs = budget_B0 + budget_q95\n", " annual_budget_bo_q95 = costs * ((0.0611)/(1-(1+0.0611)**(-20)))\n", " total_exc = B0_shift_kw + q95_shift_kw\n", " annual_kw_incentive_bo_q95 = annual_budget_bo_q95 / total_exc\n", " print(\"\")\n", " print('ALL (q95 + B0) Incentive per household € per year: ', (annual_budget_bo_q95 / 1292))\n", " print('ALL (q95 + B0) Incentive per €/kw: ', annual_kw_incentive_bo_q95)\n", " print(\"annual_budget_bo_q95\", annual_budget_bo_q95)\n", "\n", " print(\"\")\n", " print(\"False Positive freq q85: \", len(res[(res['ypred_cont'] > q85) & (res['ytrue_cont'] < q85)]) / 35040)\n", " print(\"False Positive freq q95: \", len(res[(res['ypred_cont'] > q95) & (res['ytrue_cont'] < q95)]) / 35040)\n", " print(\"False Positive freq B0: \", len(res[(res['ypred_cont'] < B0) & (res['ytrue_cont'] > B0)]) / 35040)\n" ] }, { "cell_type": "code", "execution_count": 23, "id": "dc4015d6", "metadata": {}, "outputs": [], "source": [ "# ======================================================================\n", "# Function to Plot the results of the forecast if wanted (the calls to it are currently commented out in the code at the end)\n", "# ======================================================================\n", "\n", "\n", "def plot_fcst(forecast_results, title=\"EV X, PV Y\"):\n", "\n", " plot_df = forecast_results.copy()\n", " time_col = \"clock_local\"\n", " observed_col = \"Adjusted_constumption_to_model\"\n", "\n", " plot_df[time_col] = pd.to_datetime(plot_df[time_col])\n", "\n", " q_85 = plot_df[plot_df['Adjusted_constumption_to_model'] > 0]['Adjusted_constumption_to_model'].quantile(0.85)\n", " q_95 = plot_df[plot_df['Adjusted_constumption_to_model'] > 0]['Adjusted_constumption_to_model'].quantile(0.95)\n", "\n", "\n", " plot_df[\"actual_below_0\"] = (plot_df[observed_col] < 0).astype(int)\n", " plot_df[\"actual_above_q_85\"] = (plot_df[observed_col] > q_85).astype(int)\n", " plot_df[\"actual_above_q_95\"] = (plot_df[observed_col] > q_95).astype(int)\n", "\n", " fig = make_subplots(rows=4, cols=1, shared_xaxes=True, vertical_spacing=0.06, row_heights=[0.45, 0.18, 0.18, 0.19], subplot_titles=(f\"{title}: Observed vs. predicted\", \"P(y < 0)\", \"P(y > q85)\", \"P(y > q95)\"))\n", " \n", " fig.add_trace(go.Scatter(x=plot_df[time_col], y=plot_df[observed_col], mode=\"lines\", name=\"observed\", line=dict(color=\"black\", width=2)), row=1, col=1)\n", " fig.add_trace(go.Scatter(x=plot_df[time_col], y=plot_df[\"prediction\"], mode=\"lines\", name=\"prediction\", line=dict(color=\"royalblue\", width=2)), row=1, col=1)\n", " fig.add_hline(y=0, line_dash=\"dash\", line_color=\"red\", annotation_text=\"0\", row=1, col=1)\n", " fig.add_hline(y=q_85, line_dash=\"dash\", line_color=\"orange\", annotation_text=\"q85\", row=1, col=1)\n", " fig.add_hline(y=q_95, line_dash=\"dash\", line_color=\"purple\", annotation_text=\"q95\", row=1, col=1)\n", "\n", " fig.add_trace(go.Scatter(x=plot_df[time_col], y=plot_df[\"prob_below_0\"], mode=\"lines\", name=\"P(y < 0)\", line=dict(color=\"red\", width=2), fill=\"tozeroy\", fillcolor=\"rgba(255, 0, 0, 0.15)\"), row=2, col=1)\n", " fig.add_trace(go.Scatter(x=plot_df.loc[plot_df[\"actual_below_0\"] == 1, time_col], y=plot_df.loc[plot_df[\"actual_below_0\"] == 1, \"actual_below_0\"], mode=\"markers\", name=\"actual y < 0\", marker=dict(color=\"black\", size=5),), row=2, col=1)\n", " fig.add_trace(go.Scatter(x=plot_df[time_col], y=plot_df[\"actual_above_q_85\"], mode=\"lines\", name=\"P(y > q85)\", line=dict(color=\"orange\", width=2), fill=\"tozeroy\", fillcolor=\"rgba(255, 165, 0, 0.18)\"), row=3, col=1)\n", " fig.add_trace(go.Scatter(x=plot_df.loc[plot_df[\"actual_above_q_85\"] == 1, time_col], y=plot_df.loc[plot_df[\"actual_above_q_85\"] == 1, \"actual_above_q_85\"], mode=\"markers\", name=\"actual y > q_85\", marker=dict(color=\"black\", size=5),), row=3, col=1)\n", " fig.add_trace(go.Scatter(x=plot_df[time_col], y=plot_df[\"prob_above_q_95\"], mode=\"lines\", name=\"P(y > q95)\", line=dict(color=\"purple\", width=2), fill=\"tozeroy\", fillcolor=\"rgba(128, 0, 128, 0.15)\"), row=4, col=1)\n", "\n", " fig.add_trace(go.Scatter(x=plot_df.loc[plot_df[\"actual_above_q_95\"] == 1, time_col],y=plot_df.loc[plot_df[\"actual_above_q_95\"] == 1, \"actual_above_q_95\"], mode=\"markers\", name=\"actual y > q95\", marker=dict(color=\"black\", size=5),), row=4, col=1)\n", "\n", " for r in [2, 3, 4]:\n", " fig.add_hline(y=0.5, line_dash=\"dot\", line_color=\"gray\", annotation_text=\"50%\", row=r, col=1)\n", "\n", " fig.update_yaxes(title_text=\"Value\", fixedrange=True, row=1, col=1)\n", " fig.update_yaxes(title_text=\"Probability\", range=[0, 1], tickformat=\".0%\", fixedrange=True, row=2, col=1)\n", " fig.update_yaxes(title_text=\"Probability\", range=[0, 1], tickformat=\".0%\", fixedrange=True, row=3, col=1)\n", " fig.update_yaxes(title_text=\"Probability\", range=[0, 1], tickformat=\".0%\", fixedrange=True, row=4, col=1)\n", " fig.update_xaxes(fixedrange=False, row=1, col=1)\n", " fig.update_xaxes(fixedrange=False, row=2, col=1)\n", " fig.update_xaxes(fixedrange=False, row=3, col=1)\n", " fig.update_xaxes(title_text=\"Time\", rangeslider=dict(visible=True, thickness=0.035), fixedrange=False, row=4, col=1)\n", " fig.update_layout(height=850, width=1200, hovermode=\"x unified\", template=\"plotly_white\", title=dict(text=title, x=0.5), legend=dict(orientation=\"h\", yanchor=\"bottom\", y=1.04, xanchor=\"left\", x=0), dragmode=\"zoom\")\n", "\n", " return fig" ] }, { "cell_type": "code", "execution_count": null, "id": "d01691a3", "metadata": {}, "outputs": [], "source": [ "# ======================================================================\n", "# Function to create the Load Forecast (Following the implementation described in section 4.2 of the paper)\n", "# ======================================================================\n", "\n", "\n", "def make_fcst_quantile_no_nan(df):\n", "\n", "\n", "\n", " # 1. Base setup ----------------------------------------------------------------------------------------------\n", " dd = df.copy()\n", " dd[\"clock_local\"] = pd.to_datetime(dd[\"clock_local\"])\n", "\n", " y_col = \"Adjusted_constumption_to_model\"\n", "\n", " dd.loc[dd[\"clock_local\"].dt.minute != 0, \"Prg.Globalstrahlung\"] = np.nan\n", "\n", " dd[\"Prg.Globalstrahlung\"] = (dd[\"Prg.Globalstrahlung\"].interpolate(method=\"cubic\", limit_direction=\"both\")).clip(lower=0)\n", " dd[\"Prg.Globalstrahlung\"] = (dd[\"Prg.Globalstrahlung\"].interpolate(method=\"linear\", limit_direction=\"both\").ffill().bfill())\n", "\n", " df_features = dd.copy()\n", "\n", "\n", "\n", " # 2. Feature Engineeriing ----------------------------------------------------------------------------------------------\n", " raw_hour = df_features[\"clock_local\"].dt.hour\n", " df_features[\"hour_sin\"] = np.sin(2 * np.pi * raw_hour / 24.0)\n", " df_features[\"hour_cos\"] = np.cos(2 * np.pi * raw_hour / 24.0)\n", " df_features[\"day_sin\"] = np.sin(2 * np.pi * (raw_hour * 60 + df_features[\"clock_local\"].dt.minute) / 1440.0)\n", " df_features[\"day_cos\"] = np.cos(2 * np.pi * (raw_hour * 60 + df_features[\"clock_local\"].dt.minute) / 1440.0)\n", " day_of_week = df_features[\"clock_local\"].dt.dayofweek\n", " df_features[\"week_sin\"] = np.sin(2 * np.pi * day_of_week / 7.0)\n", " df_features[\"week_cos\"] = np.cos(2 * np.pi * day_of_week / 7.0)\n", " df_features[\"is_weekend\"] = day_of_week.isin([5, 6]).astype(int)\n", "\n", " df_features[\"lag_96\"] = df_features[y_col].shift(96)\n", " df_features[\"lag_192\"] = df_features[y_col].shift(192)\n", " historical_window = df_features[y_col].shift(96)\n", " df_features[\"lag_window_mean\"] = (historical_window.rolling(window=96, min_periods=1).mean())\n", " df_features[\"lag_window_std\"] = (historical_window.rolling(window=96, min_periods=2).std().fillna(0.0))\n", " df_features[\"lag_window_max\"] = (historical_window.rolling(window=96, min_periods=1).max())\n", " df_features[\"lag_window_min\"] = (historical_window.rolling(window=96, min_periods=1).min())\n", " df_features[\"lag_window_min_dip\"] = (historical_window.rolling(window=96, min_periods=1).min())\n", " df_features[\"prev_week_negative_fraction\"] = ((historical_window < 0).astype(int).rolling(window=168, min_periods=1).mean())\n", "\n", " df_features[\"temp_lag_96\"] = df_features[\"Prg.Temperatur\"].shift(96)\n", " df_features[\"radiation_lag_96\"] = df_features[\"Prg.Globalstrahlung\"].shift(96)\n", " temp_denom = (df_features[\"Prg.Temperatur\"] + 40).clip(lower=1)\n", " df_features[\"radiation_temp_ratio\"] = (df_features[\"Prg.Globalstrahlung\"] / temp_denom)\n", " df_features[\"radiation_squared\"] = (df_features[\"Prg.Globalstrahlung\"] ** 2)\n", " df_features[\"is_sunny_hours\"] = (df_features[\"Prg.Globalstrahlung\"] > 10).astype(int)\n", " df_features[\"is_below_zero\"] = (df_features[\"lag_96\"] < 0).astype(int)\n", " df_features[\"midday_solar_intensity\"] = (df_features[\"Prg.Globalstrahlung\"] * df_features[\"day_sin\"].clip(lower=0))\n", "\n", " df_features[\"lag_96_momentum\"] = (df_features[y_col].shift(96) - df_features[y_col].shift(100))\n", " df_features[\"lag_window_negative_severity\"] = (historical_window.clip(upper=0).rolling(window=96, min_periods=1).mean())\n", " df_features[\"solar_dip_potential\"] = (df_features[\"Prg.Globalstrahlung\"] ** 1.5) / (df_features[\"Prg.Temperatur\"] + 30).clip(lower=1)\n", " df_features[\"midday_dip_vulnerability\"] = (df_features[\"midday_solar_intensity\"] * df_features[\"day_cos\"].clip(upper=0).abs())\n", " df_features[\"negative_streak_rolling\"] = ((historical_window < 0).astype(int).rolling(window=96, min_periods=1).sum())\n", "\n", " df_features = df_features.replace([np.inf, -np.inf], np.nan)\n", "\n", "\n", "\n", " # 3. Feature list ------------------------------------------------------------------\n", " lag_features = [\"lag_96\", \"lag_192\", \"lag_window_mean\", \"lag_window_std\", \"lag_window_max\", \"lag_window_min\", \"lag_window_min_dip\", \"prev_week_negative_fraction\",]\n", " time_features = [\"hour_sin\", \"hour_cos\", \"day_sin\", \"day_cos\", \"week_sin\", \"week_cos\", \"is_weekend\",]\n", " weather_features = [\"Prg.Temperatur\", \"Prg.Globalstrahlung\", \"temp_lag_96\", \"radiation_lag_96\", \"radiation_temp_ratio\", \"radiation_squared\", \"is_sunny_hours\", \"midday_solar_intensity\",]\n", " negative_dip_features = [\"lag_96_momentum\", \"lag_window_negative_severity\", \"solar_dip_potential\", \"midday_dip_vulnerability\", \"negative_streak_rolling\",]\n", " candidate_features = (lag_features + time_features + weather_features + negative_dip_features)\n", "\n", "\n", "\n", " # 4. Output columns ------------------------------------------------------------------\n", " quantiles = np.array([0.01, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99])\n", " quantile_cols = [f\"q_{int(q * 100):02d}\" for q in quantiles]\n", "\n", " dd[\"prediction\"] = np.nan\n", "\n", " for col in quantile_cols:\n", " dd[col] = np.nan\n", "\n", " dd[\"pi_50_lower\"] = np.nan\n", " dd[\"pi_50_upper\"] = np.nan\n", "\n", " dd[\"pi_80_lower\"] = np.nan\n", " dd[\"pi_80_upper\"] = np.nan\n", "\n", " dd[\"pi_90_lower\"] = np.nan\n", " dd[\"pi_90_upper\"] = np.nan\n", "\n", " dd[\"pi_98_lower\"] = np.nan\n", " dd[\"pi_98_upper\"] = np.nan\n", "\n", " dd[\"prob_below_0\"] = np.nan\n", " dd[\"prob_above_q_85\"] = np.nan\n", " dd[\"prob_above_q_95\"] = np.nan\n", "\n", "\n", "\n", " # 5. Forecast start ------------------------------------------------------------------\n", " start_test_date = dd[\"clock_local\"].iloc[8832]\n", " test_data = df_features[df_features[\"clock_local\"] >= start_test_date]\n", " positive_test_values = test_data.loc[test_data[y_col] > 0, y_col]\n", "\n", " q_85_threshold = positive_test_values.quantile(0.85)\n", " q_95_threshold = positive_test_values.quantile(0.95)\n", "\n", " test_months = test_data[\"clock_local\"].dt.to_period(\"M\").unique()\n", "\n", " print(f\"Starting predictions from: {start_test_date}\")\n", "\n", "\n", "\n", " # 6. Hyperparameter setup ------------------------------------------------------------------\n", " param_distributions = {\n", " \"learning_rate\": [0.005, 0.01, 0.015, 0.02, 0.03],\n", " \"max_iter\": [1200, 1400, 1600],\n", " \"max_depth\": [5, 7, 9, 11, None],\n", " \"min_samples_leaf\": [10, 20, 30, 50],\n", " \"l2_regularization\": [0.0, 0.1, 0.2, 0.3, 0.5],\n", " }\n", "\n", " n_iter_search = 10\n", " median_pinball_scorer = make_scorer(mean_pinball_loss, alpha=0.50, greater_is_better=False )\n", "\n", " \n", "\n", " # 7. Monthly walk-forward forecast ------------------------------------------------------------------\n", " for target_month in test_months:\n", "\n", " train_month = target_month - 1\n", "\n", " print(f\"\\nTraining for {target_month} \")\n", "\n", " test_mask = (df_features[\"clock_local\"].dt.to_period(\"M\") == target_month)\n", " test_mask = (test_mask & (df_features[\"clock_local\"] >= start_test_date))\n", " test_chunk = df_features[test_mask].copy()\n", "\n", " train_mask = (df_features[\"clock_local\"].dt.to_period(\"M\") == train_month)\n", " train_chunk = df_features[train_mask].copy()\n", " train_chunk = train_chunk.dropna(subset=[y_col])\n", "\n", " valid_features = []\n", "\n", " for feature in candidate_features:\n", "\n", " train_ok = train_chunk[feature].notna().all()\n", " test_ok = test_chunk[feature].notna().all()\n", "\n", " train_finite = np.isfinite(train_chunk[feature]).all()\n", " test_finite = np.isfinite(test_chunk[feature]).all()\n", "\n", " if train_ok and test_ok and train_finite and test_finite:\n", " valid_features.append(feature)\n", "\n", " X_train = train_chunk[valid_features]\n", " y_train = train_chunk[y_col]\n", " X_test = test_chunk[valid_features]\n", "\n", " print(f\" -> train rows: {len(train_chunk)}\")\n", " print(f\" -> test rows: {len(test_chunk)}\")\n", " print(f\" -> usable features: {len(valid_features)} / {len(candidate_features)}\")\n", "\n", "\n", " n_splits = min(3, max(2, len(train_chunk) // 500))\n", "\n", " tscv = TimeSeriesSplit(n_splits=n_splits)\n", "\n", " base_median_model = HistGradientBoostingRegressor(\n", " loss=\"quantile\",\n", " quantile=0.50,\n", " random_state=42\n", " )\n", "\n", " tuner = RandomizedSearchCV(\n", " estimator=base_median_model,\n", " param_distributions=param_distributions,\n", " n_iter=n_iter_search,\n", " scoring=median_pinball_scorer,\n", " cv=tscv,\n", " random_state=42,\n", " n_jobs=1,\n", " refit=True\n", " )\n", "\n", " tuner.fit(X_train, y_train)\n", " best_params = tuner.best_params_\n", " print(f\" -> best parameters: {best_params}\")\n", "\n", " q_preds = pd.DataFrame(index=test_chunk.index)\n", "\n", " for q in quantiles:\n", "\n", " q_col = f\"q_{int(q * 100):02d}\"\n", "\n", " q_model = HistGradientBoostingRegressor(\n", " loss=\"quantile\",\n", " quantile=q,\n", " learning_rate=best_params[\"learning_rate\"],\n", " max_iter=best_params[\"max_iter\"],\n", " max_depth=best_params[\"max_depth\"],\n", " min_samples_leaf=best_params[\"min_samples_leaf\"],\n", " l2_regularization=best_params[\"l2_regularization\"],\n", " random_state=42\n", " )\n", "\n", " q_model.fit(X_train, y_train)\n", " q_preds[q_col] = q_model.predict(X_test)\n", "\n", " # Enforce non-crossing quantiles\n", " q_array = np.maximum.accumulate(\n", " q_preds[quantile_cols].to_numpy(),\n", " axis=1\n", " )\n", "\n", " for i, q_col in enumerate(quantile_cols):\n", " dd.loc[test_chunk.index, q_col] = q_array[:, i]\n", "\n", " # q50 as point forecast\n", " q50_idx = quantile_cols.index(\"q_50\")\n", " dd.loc[test_chunk.index, \"prediction\"] = q_array[:, q50_idx]\n", "\n", "\n", " dd.loc[test_chunk.index, \"pi_50_lower\"] = q_array[:, quantile_cols.index(\"q_25\")]\n", " dd.loc[test_chunk.index, \"pi_50_upper\"] = q_array[:, quantile_cols.index(\"q_75\")]\n", "\n", " dd.loc[test_chunk.index, \"pi_80_lower\"] = q_array[:, quantile_cols.index(\"q_10\")]\n", " dd.loc[test_chunk.index, \"pi_80_upper\"] = q_array[:, quantile_cols.index(\"q_90\")]\n", "\n", " dd.loc[test_chunk.index, \"pi_90_lower\"] = q_array[:, quantile_cols.index(\"q_05\")]\n", " dd.loc[test_chunk.index, \"pi_90_upper\"] = q_array[:, quantile_cols.index(\"q_95\")]\n", "\n", " dd.loc[test_chunk.index, \"pi_98_lower\"] = q_array[:, quantile_cols.index(\"q_01\")]\n", " dd.loc[test_chunk.index, \"pi_98_upper\"] = q_array[:, quantile_cols.index(\"q_99\")]\n", "\n", " # Threshold probabilities\n", " dd.loc[test_chunk.index, \"prob_below_0\"] = [\n", " np.interp(0.0, row, quantiles, left=0.0, right=1.0)\n", " for row in q_array\n", " ]\n", "\n", " dd.loc[test_chunk.index, \"prob_above_q_85\"] = [\n", " 1.0 - np.interp(q_85_threshold, row, quantiles, left=0.0, right=1.0)\n", " for row in q_array\n", " ]\n", "\n", " dd.loc[test_chunk.index, \"prob_above_q_95\"] = [\n", " 1.0 - np.interp(q_95_threshold, row, quantiles, left=0.0, right=1.0)\n", " for row in q_array\n", " ]\n", "\n", "\n", "\n", " # 8. Final result ------------------------------------------------------------------\n", " forecast_results = dd[\n", " dd[\"clock_local\"] >= start_test_date\n", " ][\n", " [\n", " \"clock_local\",\n", " y_col,\n", " \"prediction\",\n", "\n", " \"q_01\",\n", " \"q_05\",\n", " \"q_10\",\n", " \"q_25\",\n", " \"q_50\",\n", " \"q_75\",\n", " \"q_90\",\n", " \"q_95\",\n", " \"q_99\",\n", "\n", " \"pi_50_lower\",\n", " \"pi_50_upper\",\n", " \"pi_80_lower\",\n", " \"pi_80_upper\",\n", " \"pi_90_lower\",\n", " \"pi_90_upper\",\n", " \"pi_98_lower\",\n", " \"pi_98_upper\",\n", "\n", " \"prob_below_0\",\n", " \"prob_above_q_85\",\n", " \"prob_above_q_95\",\n", " ]\n", " ].reset_index(drop=True)\n", "\n", " return forecast_results" ] }, { "cell_type": "code", "execution_count": null, "id": "99a65405", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "23e3ceb3", "metadata": {}, "outputs": [], "source": [ "# ======================================================================\n", "# Running for scenario: EV 20%, PV 25.5%\n", "# ======================================================================\n", "\n", "result_EV20_PV25_5_CONS = make_fcst_quantile_no_nan(df_EV20_PV25_5)\n", "incentives_EV20_PV25_5_CONS = make_incentives(result_EV20_PV25_5_CONS, \"low\")" ] }, { "cell_type": "code", "execution_count": null, "id": "15399422", "metadata": {}, "outputs": [], "source": [ "# ======================================================================\n", "# Running for scenario: EV 20%, PV 25.5% - PERFECT FIT\n", "# ======================================================================\n", "\n", "perfect_fit_EV20_PV25_5 = result_EV20_PV25_5_CONS.copy()\n", "perfect_fit_EV20_PV25_5['prediction'] = perfect_fit_EV20_PV25_5['Adjusted_constumption_to_model']\n", "make_incentives(perfect_fit_EV20_PV25_5, \"low\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "919c1725", "metadata": {}, "outputs": [], "source": [ "# ======================================================================\n", "# Running for scenario: EV 40%, PV 41.5%\n", "# ======================================================================\n", "\n", "result_EV40_PV41_5_CONS = make_fcst_quantile_no_nan(df_EV40_PV41_5)\n", "incentives_EV40_PV41_5_CONS = make_incentives(result_EV40_PV41_5_CONS, \"low\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "38ec7c17", "metadata": {}, "outputs": [], "source": [ "# ======================================================================\n", "# Running for scenario: EV 40%, PV 41.5% - PERFECT FIT\n", "# ======================================================================\n", "\n", "perfect_fit_EV40_PV41_5 = result_EV40_PV41_5_CONS.copy()\n", "perfect_fit_EV40_PV41_5['prediction'] = perfect_fit_EV40_PV41_5['Adjusted_constumption_to_model']\n", "make_incentives(perfect_fit_EV40_PV41_5, \"low\")\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# ======================================================================\n", "# Running for scenario: EV 75%, PV 75%\n", "# ======================================================================\n", "\n", "result_EV75_PV75_CONS = make_fcst_quantile_no_nan(df_EV75_PV75)\n", "incentives_EV75_PV75_CONS = make_incentives(result_EV75_PV75_CONS, \"low\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "069a5336", "metadata": {}, "outputs": [], "source": [ "# ======================================================================\n", "# Running for scenario: EV 75%, PV 75% - PERFECT FIT\n", "# ======================================================================\n", "\n", "perfect_fit_EV75_PV75 = result_EV75_PV75_CONS.copy()\n", "perfect_fit_EV75_PV75['prediction'] = perfect_fit_EV75_PV75['Adjusted_constumption_to_model']\n", "make_incentives(perfect_fit_EV75_PV75, \"low\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "43beaf37", "metadata": {}, "outputs": [], "source": [ "# ======================================================================\n", "# Exporting\n", "# ======================================================================\n", "\n", "DATA_DIR.mkdir(parents=True, exist_ok=True)\n", "result_EV20_PV25_5_CONS.to_csv(DATA_DIR / \"forecast_EV20_PV25_5_day_ahead.csv\", index=False)\n", "result_EV40_PV41_5_CONS.to_csv(DATA_DIR / \"forecast_EV40_PV41_5_day_ahead.csv\", index=False)\n", "result_EV75_PV75_CONS.to_csv(DATA_DIR / \"forecast_EV75_PV75_day_ahead.csv\", index=False)\n", "\n", "print(f\"Forecasts exported to: {DATA_DIR.resolve()}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "932d1cb9", "metadata": {}, "outputs": [], "source": [ "plot_EV20_PV25_5_CONS = plot_fcst(result_EV20_PV25_5_CONS, \"EV 20%, PV 25.5%\")\n", "plot_EV40_PV41_5_CONS = plot_fcst(result_EV40_PV41_5_CONS, \"EV 40%, PV 41.5%\")\n", "plot_EV75_PV75_CONS = plot_fcst(result_EV20_PV25_5_CONS, \"EV 75%, PV 75%\")\n", "\n", "plot_EV20_PV25_5_CONS\n", "plot_EV40_PV41_5_CONS\n", "plot_EV75_PV75_CONS\n" ] } ], "metadata": { "kernelspec": { "display_name": "gdal_env", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.14" } }, "nbformat": 4, "nbformat_minor": 5 }