Solar PV Forecasting with ML

Forecast solar power output using SARIMAX, Random Forest, and XGBoost.

Introduction

30 Second Summary

Solar panels on rooftops generate power based on weather conditions that change hour to hour. Grid operators need accurate forecasts of this production to balance supply and demand without blackouts or wasted energy.

In this project, you will build a machine learning pipeline that predicts hourly solar power output for a rooftop installation in Lyon, France. You will compare three forecasting approaches and discover why gradient boosting dominates energy forecasting.

What You'll Build

You run your notebook, feed it a year of real weather data, and watch three models compete to predict when the sun delivers power and when clouds steal it away.

By the end of this project, you'll have:

  • A simulated PV production dataset for a 5 kWc system driven by real historical weather data from Open-Meteo, showing clear daily and seasonal patterns plus cloud-induced dips you can visualize.
  • A persistence baseline that exposes why naive "same as yesterday" forecasts fail spectacularly on days following weather transitions.
  • A side-by-side model comparison of SARIMAX, Random Forest, and XGBoost with MAE, RMSE, and skill score metrics proving which model wins and a feature importance plot showing why.
  • Secret Mission: Generate a 48-hour forward forecast and identify surplus grid injection hours where predicted PV output exceeds household consumption.

Are there any prerequisites?

You should be comfortable with Python fundamentals and have a basic understanding of pandas and NumPy. Familiarity with machine learning concepts like train/test splits and regression will help, but is not strictly required.

Before We Start

Before We Start...

Before we dive in, let's lock in what you're about to build and why it matters. Solar power is intermittent, and grid operators need accurate forecasts to balance supply and demand without blackouts or waste.

Set Up Your Environment

This project combines three different Python ecosystems: PV physics modeling, classical statistics, and machine learning. All of these libraries need to coexist in one environment before you can build anything.

In this step, you will install Python 3.12+, grab all the required packages, and launch a Jupyter Notebook where the forecasting work happens.

In this step, get ready to:
  • Install Python 3.12+ and verify the version.
  • Install PV modeling, statistical, and ML libraries.
  • Launch Jupyter Notebook and confirm all imports work.
Install Python 3.12+

This project requires Python 3.12 or higher because XGBoost 3.3.0 (the gradient boosting library you will use) only supports Python 3.12+.

  • Open your terminal by pressing Cmd+Space (macOS) or the Windows key (Windows), typing Terminal, and pressing Enter.
  • Check your Python version by running this command:
python --version

What do you see after running the command?

✔️ I see version 3.12 or higher

Your Python version is compatible. Move on to the next substep.

ⓧ I see an older version

Your system has an older Python installed that does not meet the 3.12+ requirement. On macOS, the system Python is often an older version that shadows newer installs.

  • Download the latest Python 3.12 installer from python.org/downloads.
  • Run the installer and follow the prompts.
  • On macOS, the new install may not override the system python command. Try running python3.12 --version to confirm the new install is accessible.
  • If python3.12 works, use python3.12 and pip3.12 in place of python and pip for the rest of this project.

Re-run the version check to confirm you now see 3.12 or higher.

ⓧ Command not found

Python is not installed on your machine yet.

  • Go to python.org/downloads and download the Python 3.12 installer for your operating system.
  • Run the installer. On Windows, check the box that says Add Python to PATH before clicking Install Now.
  • Close and reopen your terminal after installation.
  • Run python --version again to confirm you see 3.12 or higher.

macOS alternative command

On macOS, the command may be python3 instead of python. Try python3 --version if the plain python command is not found.

Install all project dependencies

You need libraries spanning three domains: pvlib for solar physics simulation, scikit-learn and xgboost for machine learning, and statsmodels for classical time series modeling.

  • Install all required packages by running this command:
pip install pvlib==0.15.2 scikit-learn==1.9.0 xgboost==3.3.0 statsmodels==0.14.6 matplotlib pandas numpy requests jupyter

What does this command do?

  • pvlib==0.15.2 provides physics-based PV system modeling (simulating solar panel output from weather data).
  • scikit-learn==1.9.0 provides the Random Forest regressor and evaluation metrics.
  • xgboost==3.3.0 provides the gradient boosting model that will outperform other approaches.
  • statsmodels==0.14.6 provides the SARIMAX time series model.
  • matplotlib, pandas, numpy, and requests handle data manipulation, visualization, and API calls.
  • jupyter gives you the interactive notebook environment for running code cell by cell.

Before you run the next command, what do you think will happen if you try to import all three ML libraries in one line?

  • Verify the installation succeeded by running this command:
python -c "import pvlib; import xgboost; import statsmodels; print('All imports OK')"

You should see All imports OK printed with no errors.

Seeing an import error?

  • If you see a ModuleNotFoundError, the pip you used may not match the python you are running. Try pip3 install ... or python -m pip install ... instead.
  • On macOS, if you installed Python 3.12 separately, use python3.12 -m pip install ... to ensure packages install for the correct version.

Help me troubleshoot a Python import error after pip install

Launch Jupyter and create your notebook

Jupyter Notebook lets you run code in small cells and see results immediately. This is ideal for data science workflows where you build up analysis piece by piece.

  • Move to your Desktop by running this command:
cd ~/Desktop
  • Create a project folder and move into it by running these commands:
mkdir pv-forecasting
cd pv-forecasting
  • Launch Jupyter Notebook by running:
jupyter notebook

What happens when you run this?

Jupyter starts a local web server and opens your browser to http://localhost:8888. Your terminal will show server logs and remain occupied. Leave it running.

If your browser does not open automatically, look for a URL in the terminal output that starts with http://localhost:8888 and paste it into your browser.

  • In the Jupyter interface, click New in the top right corner.
  • Select Python 3 from the dropdown. A new notebook opens in a separate tab.
  • Click the title at the top of the page (it says Untitled) and rename it to pv_forecasting_lyon.

Your notebook is now saved as pv_forecasting_lyon.ipynb inside the pv-forecasting folder on your Desktop.

Write and run the imports cell

The first cell of every data science notebook imports all the libraries and sets global configuration. This confirms your environment is correctly wired.

  • In the first empty cell of your notebook, paste this import block:
# Cell 1: Imports and Configuration
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import requests
from datetime import datetime, timedelta

import pvlib
from pvlib.pvsystem import PVSystem
from pvlib.location import Location
from pvlib.modelchain import ModelChain
from pvlib.temperature import TEMPERATURE_MODEL_PARAMETERS

from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error

import xgboost as xgb
from statsmodels.tsa.statespace.sarimax import SARIMAX

import warnings
warnings.filterwarnings('ignore')

What do these imports cover?

  • pvlib classes (PVSystem, Location, ModelChain) simulate solar panel output from weather data.
  • sklearn provides the Random Forest model and evaluation metrics (MAE, RMSE).
  • xgb is the XGBoost gradient boosting library.
  • SARIMAX is the classical time series model that handles both autocorrelation and exogenous weather inputs.
  • warnings.filterwarnings('ignore') suppresses deprecation warnings that would clutter your notebook output.
  • Still in the same cell, add this configuration block below the imports:
# Configuration
LATITUDE = 45.76
LONGITUDE = 4.83
TIMEZONE = 'Europe/Paris'
SYSTEM_CAPACITY_KW = 5.0

print("All imports OK")
print(f"pvlib version: {pvlib.__version__}")
print(f"XGBoost version: {xgb.__version__}")

What do these constants represent?

These define the PV system you will forecast for: a 5 kWc rooftop installation in Lyon, France (latitude 45.76, longitude 4.83). Lyon has a good solar resource with significant cloud variability, making it ideal for testing forecast models.

Before you run this cell, what do you expect to see printed if all libraries installed correctly?

  • Press Shift+Enter to run the cell.

You should see three lines of output: All imports OK, then the pvlib version (0.15.2), then the XGBoost version (3.3.0).

Seeing a ModuleNotFoundError in the notebook?

  • This means Jupyter is using a different Python than the one you installed packages for. Close Jupyter (Ctrl+C in the terminal), then relaunch it using python -m jupyter notebook to force the correct Python.
  • If you installed Python 3.12 separately on macOS, try python3.12 -m jupyter notebook instead.

Help me fix a ModuleNotFoundError inside Jupyter Notebook

✔️ Awesome, I've got everything!

Great. Confirm you have saved the notebook (Cmd+S or Ctrl+S).

ⓧ I'd like to double check the full code

Your Cell 1 should contain two sections that you pasted above, in this order:

  • The imports block (starting with # Cell 1: Imports and Configuration through warnings.filterwarnings('ignore')).
  • The configuration and print block (starting with # Configuration through the final print statement).

Both sections go in the same cell. Scroll up and confirm each block is present and in order.

Your environment is ready with all the PV modeling, statistics, and ML libraries working together. Next up, you will pull real weather data from Open-Meteo and simulate a full year of solar power production.

Fetch Weather Data and Simulate PV Output

Your libraries are installed and your notebook is running. Now you need the raw material for any forecasting model: data.

Machine learning models learn patterns from historical examples. For PV forecasting, that means pairing real weather conditions (irradiance, temperature, clouds) with the power a solar panel actually produced under those conditions. In this step, you fetch one year of hourly weather for Lyon from Open-Meteo, simulate "ground truth" PV output using a physics model, and visualize the result.

In this step, get ready to:
  • Fetch one year of hourly weather data (GHI, DNI, DHI, temperature, wind speed, cloud cover) from Open-Meteo for Lyon.
  • Simulate hourly AC power output for a 5 kWc rooftop PV system using pvlib.
  • Plot the full-year production curve and a weekly zoom showing cloud impact.
Fetch historical weather data from Open-Meteo

Open-Meteo provides free historical weather data based on ERA5 reanalysis at hourly resolution. You need irradiance components (GHI, DNI, DHI) to drive the PV model, plus temperature and wind speed which affect panel efficiency.

  • In your notebook, create a new cell below Cell 1. Paste the following code to set up the API request:
# Cell 2: Fetch Historical Weather Data from Open-Meteo
end_date = "2025-12-31"
start_date = "2025-01-01"

url = "https://archive-api.open-meteo.com/v1/archive"
params = {
    "latitude": LATITUDE,
    "longitude": LONGITUDE,
    "start_date": start_date,
    "end_date": end_date,
    "hourly": [
        "shortwave_radiation",
        "direct_radiation",
        "diffuse_radiation",
        "direct_normal_irradiance",
        "temperature_2m",
        "wind_speed_10m",
        "cloud_cover"
    ],
    "timezone": TIMEZONE
}

response = requests.get(url, params=params)
data = response.json()

What does this code do?

This builds and sends a GET request to the Open-Meteo Historical Weather API. The key pieces:

  • `shortwave_radiation` is the total Global Horizontal Irradiance (GHI) in W/m2. This is the primary driver of PV output.
  • `direct_normal_irradiance` (DNI) and `diffuse_radiation` (DHI) split the total irradiance into its beam and scattered components. pvlib needs all three to model a tilted panel.
  • `temperature_2m` and `wind_speed_10m` affect panel cell temperature and thus efficiency.
  • Directly below that code in the same cell, paste this to parse the JSON response into a DataFrame:
hourly = data["hourly"]
weather_df = pd.DataFrame({
    "datetime": pd.to_datetime(hourly["time"]),
    "ghi": hourly["shortwave_radiation"],
    "dhi": hourly["diffuse_radiation"],
    "dni": hourly["direct_normal_irradiance"],
    "temp_air": hourly["temperature_2m"],
    "wind_speed": hourly["wind_speed_10m"],
    "cloud_cover": hourly["cloud_cover"]
})
weather_df = weather_df.set_index("datetime")
weather_df.index = weather_df.index.tz_localize(TIMEZONE)

print(f"Weather data shape: {weather_df.shape}")
print(f"Date range: {weather_df.index[0]} to {weather_df.index[-1]}")
weather_df.head()

What does this parsing do?

The API returns JSON with a nested `hourly` dictionary. Each key holds an array of 8,760 values (one per hour of the year). You reshape these arrays into a pandas DataFrame with a timezone-aware datetime index, which pvlib requires for solar position calculations.

The column rename from `shortwave_radiation` to `ghi` makes the DataFrame match standard solar energy terminology.

  • Run Cell 2 by pressing **Shift+Enter**.

You should see output showing the DataFrame shape (8760, 6) and a date range covering all of 2025, plus a preview of the first few rows.

Seeing an error from the API call?

  • If you see a `ConnectionError`, check your internet connection. Open-Meteo requires no API key, but you need network access.
  • If the response contains `"error": true`, the API rejected your parameters. Double-check that `start_date` and `end_date` use the `YYYY-MM-DD` format exactly as shown.
  • If you see a `KeyError` on `hourly`, print `data` to inspect the raw response. A misspelled variable name in the `hourly` list causes it to be silently ignored.

Help me debug my Open-Meteo API call in Jupyter Notebook.

Simulate PV output with pvlib

Weather data alone is not a training target. You need the power a real panel would produce under those conditions. The pvlib library models the full physics chain from sunlight hitting a tilted surface through to AC power out of the inverter.

  • Create a new cell below Cell 2. Paste this code to configure and run the PV simulation:
# Cell 3: PV System Simulation with pvlib
location = Location(
    latitude=LATITUDE,
    longitude=LONGITUDE,
    tz=TIMEZONE,
    altitude=170
)

# Define a 5 kWc system: fixed tilt 30 deg, south-facing (azimuth=180)
temperature_params = TEMPERATURE_MODEL_PARAMETERS["sapm"]["open_rack_glass_glass"]

system = PVSystem(
    surface_tilt=30,
    surface_azimuth=180,
    module_parameters={"pdc0": SYSTEM_CAPACITY_KW * 1000, "gamma_pdc": -0.004},
    inverter_parameters={"pdc0": SYSTEM_CAPACITY_KW * 1000 * 1.1, "eta_inv_nom": 0.96},
    temperature_model_parameters=temperature_params
)

mc = ModelChain.with_pvwatts(system, location)

What is this PV system?

You are modeling a typical French rooftop installation: 5 kWc (kilowatt-peak) of panels at 30 degrees tilt, facing due south (azimuth 180). The PVWatts model simplifies the electrical chain into two parameters:

  • `pdc0` = 5000 W nameplate DC capacity.
  • `gamma_pdc` = -0.004 temperature coefficient (power drops 0.4% per degree above 25 C).
  • `eta_inv_nom` = 0.96 means the inverter converts 96% of DC to AC.

ModelChain wires together irradiance transposition, cell temperature modeling, and power conversion in one call.

  • Directly below that code in the same cell, paste this to run the simulation and store results:
# Prepare weather for pvlib
pvlib_weather = weather_df[["ghi", "dhi", "dni", "temp_air", "wind_speed"]].copy()
pvlib_weather.columns = ["ghi", "dhi", "dni", "temp_air", "wind_speed"]

mc.run_model(pvlib_weather)

# Extract AC power in kW
weather_df["pv_power_kw"] = mc.results.ac / 1000.0
weather_df["pv_power_kw"] = weather_df["pv_power_kw"].clip(lower=0)

print(f"Max PV power: {weather_df['pv_power_kw'].max():.2f} kW")
print(f"Annual production: {weather_df['pv_power_kw'].sum():.0f} kWh")

Why clip to zero?

The PVWatts model can produce tiny negative values at night due to floating-point arithmetic. Clipping to zero ensures the output is physically meaningful. The annual production total (around 5,500-6,500 kWh) is consistent with a real 5 kWc system in Lyon.

  • Run Cell 3 by pressing **Shift+Enter**.

You should see the max PV power (close to 5 kW on the clearest summer day) and the annual production total printed below the cell.

Seeing a ModelChain error?

  • If you see a `KeyError` about missing columns, check that `weather_df` from Cell 2 has columns named exactly `ghi`, `dhi`, `dni`, `temp_air`, `wind_speed`. Re-run Cell 2 if needed.
  • If you see `TypeError` about timezone, confirm Cell 2's `tz_localize` line ran successfully. pvlib requires timezone-aware indices.

Help me fix my pvlib ModelChain simulation error.

Visualize the production curve

Numbers alone do not tell you whether the simulation is realistic. A plot reveals the expected seasonal arc (high summer, low winter) and the stochastic cloud dips that make forecasting difficult.

  • Create a new cell below Cell 3. Paste this visualization code:
# Cell 4: Visualize PV Production
fig, axes = plt.subplots(2, 1, figsize=(14, 8))

# Full year
axes[0].plot(weather_df.index, weather_df["pv_power_kw"], linewidth=0.3, color="orange")
axes[0].set_title("Simulated PV Production - Lyon 5 kWc - Full Year 2025")
axes[0].set_ylabel("Power (kW)")
axes[0].set_xlabel("Date")

# One week zoom (June)
week_start = "2025-06-09"
week_end = "2025-06-15"
week_data = weather_df.loc[week_start:week_end]
axes[1].plot(week_data.index, week_data["pv_power_kw"], color="orange", label="PV Power")
axes[1].fill_between(week_data.index, week_data["pv_power_kw"], alpha=0.3, color="orange")
axes[1].set_title("Weekly Zoom - June 2025 (notice cloud impact)")
axes[1].set_ylabel("Power (kW)")
axes[1].set_xlabel("Hour")
axes[1].legend()

plt.tight_layout()
plt.show()

What do the two panels show?

  • The top panel shows the full year. You can see the seasonal envelope: production peaks in June/July and drops in December/January.
  • The bottom panel zooms into one June week. Clear days show smooth bell curves reaching nearly 5 kW at solar noon. Cloudy days show irregular dips and reduced peaks. This variability is exactly what your forecasting models need to predict.

Before you run this cell, what do you expect the summer peak power to be, given a 5 kWc system?

  • Run Cell 4 by pressing **Shift+Enter**.

You should see a two-panel matplotlib plot. The top shows the full-year seasonal arc with dense orange traces. The bottom shows a June week with clear bell-curve days and cloud-disrupted days side by side.

Plot not appearing?

  • If you see a blank output, make sure you ran Cells 1-3 in order first. Cell 4 depends on `weather_df` having the `pv_power_kw` column.
  • If the plot appears but looks empty, check that `weather_df['pv_power_kw'].max()` is greater than 0 (run it in a new cell to confirm).

Help me debug why my matplotlib plot is not showing PV production data.

✔️ Awesome, I've got everything!

Great. Make sure all four cells are saved in your notebook.

ⓧ I'd like to double check the full code

Scroll through your notebook and confirm you have these cells in order, each producing output:

  • Cell 1: Imports and Configuration (prints "All imports OK" plus version numbers).
  • Cell 2: Fetch Historical Weather Data from Open-Meteo (prints DataFrame shape and date range, displays first 5 rows).
  • Cell 3: PV System Simulation with pvlib (prints max PV power and annual production).
  • Cell 4: Visualize PV Production (displays a two-panel matplotlib plot).

You now have a complete training dataset: real weather inputs paired with physics-simulated PV output. Next up, you will engineer domain-specific features and discover why a naive persistence forecast falls short on cloudy days.

Engineer Features and Build a Persistence Baseline

You have a full year of simulated PV production. But if you tried to forecast power output using only raw weather columns like GHI and temperature, you'd miss the predictable patterns of solar geometry that dominate production.

In this step, you'll engineer domain-specific features that separate the predictable solar cycle from unpredictable weather effects. Then you'll test the simplest possible forecast approach to see exactly why smarter models are needed.

In this step, get ready to:
  • Engineer domain-specific solar features including clear-sky index, solar zenith angle, and cyclical time encodings.
  • Split data into train and test sets preserving temporal order.
  • Build a persistence baseline and visualize where it fails on weather transitions.
Engineer solar features

The clear-sky index is the single most important feature in PV forecasting. It's the ratio of actual GHI to what GHI would be on a perfectly clear day. A value near 1 means clear skies. A value near 0 means heavy cloud cover.

This separates the predictable solar geometry (which pvlib already models) from the unpredictable weather effect you actually need to forecast.

  • Create a new cell in your notebook (Cell 5) and add the following feature engineering code:
# Cell 5: Feature Engineering
# Clear-sky model for clear-sky index
clearsky = location.get_clearsky(weather_df.index)
weather_df["ghi_clearsky"] = clearsky["ghi"]

# Clear-sky index (kt): ratio of actual to clear-sky GHI
weather_df["clearsky_index"] = np.where(
    weather_df["ghi_clearsky"] > 10,
    weather_df["ghi"] / weather_df["ghi_clearsky"],
    0
)
weather_df["clearsky_index"] = weather_df["clearsky_index"].clip(0, 1.5)

# Solar position
solar_pos = location.get_solarposition(weather_df.index)
weather_df["solar_zenith"] = solar_pos["zenith"]

# Cyclical time features
weather_df["hour_sin"] = np.sin(2 * np.pi * weather_df.index.hour / 24)
weather_df["hour_cos"] = np.cos(2 * np.pi * weather_df.index.hour / 24)
weather_df["month_sin"] = np.sin(2 * np.pi * weather_df.index.month / 12)
weather_df["month_cos"] = np.cos(2 * np.pi * weather_df.index.month / 12)

print("Features engineered:")
print(weather_df[["clearsky_index", "solar_zenith", "hour_sin", "month_sin"]].describe())

What does this code do?

This cell engineers four types of domain features from the raw data:

  • clearsky_index normalizes GHI against the theoretical clear-sky maximum. This isolates the weather effect from predictable solar geometry.
  • solar_zenith tells the model where the sun is in the sky at each hour.
  • Cyclical sin/cos encodings for hour and month let the model understand that hour 23 is close to hour 0, not far away. Plain integers would break this circular relationship.
  • Run the cell.

You should see a describe() table showing statistics for your new features. The clearsky_index mean should be well below 1.0, reflecting Lyon's mix of sunny and cloudy days.

Seeing a KeyError or empty DataFrame?

  • Make sure you ran Cells 1 through 4 first. The location object and weather_df must exist before this cell can run.
  • Check that your weather_df index is timezone-aware. The get_clearsky() method requires a tz-aware DatetimeIndex.

Help me debug my feature engineering cell

Build a persistence baseline and split the data

The simplest possible PV forecast is persistence: assume tomorrow's production at each hour equals today's production at the same hour. Every more complex model must beat this baseline to justify its existence.

Before you run the next cell, consider: do you think repeating yesterday's production will produce an accurate forecast? What happens on days when the weather changes suddenly?

  • Create a new cell (Cell 6) and add the persistence baseline code:
# Cell 6: Persistence Baseline
# Persistence: forecast = same hour yesterday
weather_df["persistence"] = weather_df["pv_power_kw"].shift(24)

# Drop NaN rows
df_clean = weather_df.dropna().copy()

# Filter daytime only (solar zenith < 85) for meaningful evaluation
df_day = df_clean[df_clean["solar_zenith"] < 85].copy()

# Train/test split: first 10 months train, last 2 months test
split_date = "2025-11-01"
train = df_day.loc[:split_date]
test = df_day.loc[split_date:]

# Persistence metrics on test
persistence_mae = mean_absolute_error(test["pv_power_kw"], test["persistence"])
persistence_rmse = np.sqrt(mean_squared_error(test["pv_power_kw"], test["persistence"]))

print(f"Persistence Baseline (daytime hours):")
print(f"  MAE:  {persistence_mae:.3f} kW")
print(f"  RMSE: {persistence_rmse:.3f} kW")
print(f"  Test samples: {len(test)}")

What does this code do?

  • shift(24) copies yesterday's power value at the same hour as the "forecast" for today.
  • Filtering to solar_zenith < 85 keeps only daytime hours. Evaluating on nighttime (when power is always zero) would inflate accuracy artificially.
  • The train/test split uses the first 10 months for training and the last 2 months (November-December) for testing. Temporal order is preserved so no future data leaks into the training set.
  • Run the cell.

You should see a MAE above 0.5 kW. That means on average, the persistence forecast is off by more than half a kilowatt every hour. For a 5 kW system, that's a significant error.

What are MAE and RMSE?

MAE (Mean Absolute Error) is the average of all absolute differences between predicted and actual values. It gives you a straightforward "how far off am I on average" in kW.

RMSE (Root Mean Squared Error) penalizes large errors more heavily. If RMSE is much larger than MAE, it means your model occasionally makes very large mistakes.

Getting NaN or zero for the metrics?

  • Confirm Cell 5 ran successfully. The solar_zenith column must exist for the daytime filter to work.
  • Check that pv_power_kw exists in weather_df. This was created in Cell 3.

Help me debug my persistence baseline metrics

Visualize the persistence failure

The numbers tell you persistence is bad. A plot shows you exactly why it's bad.

  • Create a new cell (Cell 7) and add the visualization code:
# Cell 7: Visualize Persistence Failure
fig, ax = plt.subplots(figsize=(14, 4))
plot_start = "2025-11-03"
plot_end = "2025-11-07"
plot_data = test.loc[plot_start:plot_end]

ax.plot(plot_data.index, plot_data["pv_power_kw"], label="Actual", color="orange", linewidth=2)
ax.plot(plot_data.index, plot_data["persistence"], label="Persistence (yesterday)", color="gray", linestyle="--")
ax.set_title("Persistence Fails on Weather Transitions")
ax.set_ylabel("Power (kW)")
ax.legend()
plt.tight_layout()
plt.show()

What does this plot show?

This zooms into a 4-day window in November. The orange line is actual PV production. The gray dashed line is the persistence forecast (yesterday's values shifted forward 24 hours).

On days where the weather is similar to the previous day, persistence looks fine. But on days following a weather transition (sunny to cloudy, or cloudy to sunny), the gray line lags behind badly.

  • Run the cell.

You should see the plot with the persistence forecast (gray dashed line) clearly lagging behind the actual production (orange) on days where the weather changed.

This is the gap you'll close in the next steps. Persistence ignores both the temporal autocorrelation in the data and the weather forecast information you have available. A model that uses these should do much better.

Why is persistence the designed obstacle?

In solar forecasting, persistence is the standard "no-skill" reference. Any model you build must beat it, or it's adding complexity for nothing. The skill score metric you'll use later is defined as 1 minus (model error / persistence error).

A positive skill score means the model adds real value. This is how energy companies evaluate forecast providers.

Plot looks flat or empty?

  • Make sure Cell 6 ran without errors. The test DataFrame must exist with both pv_power_kw and persistence columns.
  • If the plot only shows nighttime zeros, check that your test data includes daytime hours. The filter solar_zenith < 85 should have kept only hours with the sun above the horizon.

Help me fix my persistence failure visualization

✔️ Awesome, I've got everything!

Great. Double check you've saved your notebook with all three new cells (5, 6, and 7) present and run.

ⓧ I'd like to double check the full code

Your notebook should now have Cells 1 through 7. Scroll up through your notebook and confirm these three cells are present in order, matching the code you added above:

  • Cell 5: Feature Engineering with clear-sky index, solar zenith, and cyclical time features.
  • Cell 6: Persistence Baseline with the shift(24), daytime filter, train/test split, and MAE/RMSE output.
  • Cell 7: Visualize Persistence Failure with the 4-day comparison plot.

You now have a clear picture of the problem: persistence can't handle weather transitions, and it sets the bar your ML models must beat. Next up, you'll build a SARIMAX model that uses both temporal autocorrelation and exogenous weather inputs to close this gap.

Forecast with SARIMAX

The persistence baseline proved that "just use yesterday" fails whenever the weather changes between days. It ignores two critical signals: the temporal autocorrelation within a day's power curve, and the exogenous weather variables that drive production.

SARIMAX addresses both. The "AR" and "MA" components model how each hour relates to prior hours. The seasonal component captures the repeating 24-hour diurnal cycle. The "X" means it can accept exogenous inputs like GHI and temperature. In this step, you fit a SARIMAX model and measure how much it improves over persistence.

In this step, get ready to:
  • Fit a SARIMAX model on downsampled training data with weather inputs.
  • Generate predictions on the test period and compare MAE and RMSE to persistence.
  • Evaluate where SARIMAX still struggles with rapid intra-day variability.
Fit the SARIMAX model on downsampled training data

SARIMAX fitting is computationally expensive on large hourly datasets. To keep the fit under 2 minutes in your notebook, you downsample the training data to every 3rd row while still retaining enough signal to learn the temporal patterns.

  • Create a new cell below Cell 7 and add the following code to fit the SARIMAX model:
# Cell 8: SARIMAX Model
train_sarimax = train.iloc[::3][["pv_power_kw", "ghi", "temp_air"]].copy()

exog_train = train_sarimax[["ghi", "temp_air"]]
endog_train = train_sarimax["pv_power_kw"]

print("Fitting SARIMAX (this may take 1-2 minutes)...")
model = SARIMAX(
    endog_train,
    exog=exog_train,
    order=(1, 0, 1),
    seasonal_order=(1, 0, 1, 24),
    enforce_stationarity=False,
    enforce_invertibility=False
)
results = model.fit(disp=False, maxiter=100)
print("SARIMAX fitted.")
print(f"AIC: {results.aic:.0f}")

What does this code do?

train.iloc[::3] takes every 3rd row from the training set, reducing the dataset size so SARIMAX can fit in a reasonable time.

order=(1, 0, 1) sets one autoregressive lag, no differencing (PV production is seasonal but stationary within a year), and one moving-average term. seasonal_order=(1, 0, 1, 24) adds a seasonal layer with period 24, capturing the repeating daily power curve.

The exog parameter feeds GHI and temperature as external weather inputs. This is the "X" in SARIMAX, and it is what lets the model react to actual weather conditions instead of relying purely on past power values.

  • Run the cell. You will see "Fitting SARIMAX (this may take 1-2 minutes)..." printed while the model trains.
  • Once complete, you should see "SARIMAX fitted." followed by an AIC value (a model quality score where lower is better).

Model taking more than 5 minutes?

Make sure you used train.iloc[::3] to downsample. Without this, SARIMAX attempts to fit thousands of hourly observations and can run for 10+ minutes.

If you see a convergence warning, that is expected with maxiter=100. The model still produces usable predictions.

Help me debug my SARIMAX fitting issue

Generate predictions and compare to persistence

With the model fitted, you can now forecast on the test period using the actual weather observations as exogenous inputs. This simulates a "perfect weather forecast" scenario, isolating how well SARIMAX converts weather knowledge into power predictions.

  • Create a new cell below Cell 8 and add the following code to generate SARIMAX predictions and compute error metrics:
# Cell 9: SARIMAX Prediction on Test Set
exog_test = test[["ghi", "temp_air"]]

forecast_sarimax = results.get_forecast(
    steps=len(test),
    exog=exog_test.values
)
test_copy = test.copy()
test_copy["sarimax_pred"] = forecast_sarimax.predicted_mean.clip(lower=0).values

sarimax_mae = mean_absolute_error(test_copy["pv_power_kw"], test_copy["sarimax_pred"])
sarimax_rmse = np.sqrt(mean_squared_error(test_copy["pv_power_kw"], test_copy["sarimax_pred"]))

print(f"SARIMAX Results:")
print(f"  MAE:  {sarimax_mae:.3f} kW")
print(f"  RMSE: {sarimax_rmse:.3f} kW")
print(f"  Improvement over persistence: {(1 - sarimax_mae/persistence_mae)*100:.1f}%")

What does this code do?

results.get_forecast(steps=len(test), exog=exog_test.values) asks the fitted model to predict forward for as many steps as the test set has rows, feeding in the real weather data as the exogenous input.

.clip(lower=0) ensures predictions never go negative. A solar panel cannot produce negative power, but the model's linear components can occasionally predict below zero at night or during heavy cloud.

The improvement percentage uses the formula 1 - (model_MAE / persistence_MAE). This is the skill score concept you will use for the final comparison. A positive percentage means SARIMAX is beating persistence.

Before you run this cell, do you think SARIMAX will improve on persistence by more or less than 20%? Consider that persistence had no weather awareness at all.

  • Run the cell.

You should see the SARIMAX MAE printed below the persistence MAE value, confirming it has dropped. The improvement percentage shows how much closer to zero the error moved.

Where does SARIMAX still struggle?

SARIMAX captures linear relationships between weather and power output. It handles the broad daily shape and responds to GHI changes. But rapid intra-day variability, like a cloud passing over for 30 minutes mid-afternoon, creates non-linear interactions that a linear model cannot fully track.

This is why tree-based models (the next step) outperform SARIMAX in PV forecasting. They can learn the complex, non-linear mapping between multiple weather features and power output simultaneously.

SARIMAX MAE is higher than persistence?

Check that exog_test uses the same columns (ghi and temp_air) in the same order as exog_train. Mismatched column order causes the model to apply GHI coefficients to temperature values.

Verify that Cell 8 completed with "SARIMAX fitted." printed. If it was interrupted mid-fit, the model parameters are unreliable.

Help me understand why my SARIMAX results are worse than persistence

✔️ Awesome, I've got everything!

Great. Double check you have both Cell 8 and Cell 9 saved in your notebook.

ⓧ I'd like to double check the full code

Here is Cell 8:

# Cell 8: SARIMAX Model
train_sarimax = train.iloc[::3][["pv_power_kw", "ghi", "temp_air"]].copy()

exog_train = train_sarimax[["ghi", "temp_air"]]
endog_train = train_sarimax["pv_power_kw"]

print("Fitting SARIMAX (this may take 1-2 minutes)...")
model = SARIMAX(
    endog_train,
    exog=exog_train,
    order=(1, 0, 1),
    seasonal_order=(1, 0, 1, 24),
    enforce_stationarity=False,
    enforce_invertibility=False
)
results = model.fit(disp=False, maxiter=100)
print("SARIMAX fitted.")
print(f"AIC: {results.aic:.0f}")

Here is Cell 9:

# Cell 9: SARIMAX Prediction on Test Set
exog_test = test[["ghi", "temp_air"]]

forecast_sarimax = results.get_forecast(
    steps=len(test),
    exog=exog_test.values
)
test_copy = test.copy()
test_copy["sarimax_pred"] = forecast_sarimax.predicted_mean.clip(lower=0).values

sarimax_mae = mean_absolute_error(test_copy["pv_power_kw"], test_copy["sarimax_pred"])
sarimax_rmse = np.sqrt(mean_squared_error(test_copy["pv_power_kw"], test_copy["sarimax_pred"]))

print(f"SARIMAX Results:")
print(f"  MAE:  {sarimax_mae:.3f} kW")
print(f"  RMSE: {sarimax_rmse:.3f} kW")
print(f"  Improvement over persistence: {(1 - sarimax_mae/persistence_mae)*100:.1f}%")

SARIMAX proves that incorporating weather information and temporal structure beats naive persistence. But its linear assumptions limit how much further it can go. Next, you will unleash tree-based models that handle non-linear feature interactions natively and use all your engineered features at once.

Train Tree-Based Models and Compare All Forecasters

SARIMAX improved over persistence by modeling temporal autocorrelation with weather inputs. But it is fundamentally a linear model. It cannot capture the complex non-linear interactions between irradiance, temperature, cloud cover, and PV output that happen simultaneously.

Tree-based models like Random Forest and XGBoost handle these non-linear interactions naturally and can use all 12 engineered features at once. In this step, you'll train both and pit all four approaches against each other in a final comparison.

In this step, get ready to:
  • Train Random Forest and XGBoost regressors on all engineered features.
  • Compare all 4 models with MAE, RMSE, and skill score metrics.
  • Visualize XGBoost feature importance to understand what drives PV forecasting accuracy.
Train Random Forest and XGBoost regressors

Unlike SARIMAX which only used GHI and temperature, tree-based models can ingest every feature you engineered. They learn decision boundaries that naturally capture interactions like "high GHI + low cloud cover + noon = peak production."

  • Create a new cell in your notebook and add the following code to define the feature set and train a Random Forest:
# Cell 10: Tree-Based Models
# Feature columns for ML models
feature_cols = [
    "ghi", "dni", "dhi", "temp_air", "wind_speed", "cloud_cover",
    "solar_zenith", "clearsky_index",
    "hour_sin", "hour_cos", "month_sin", "month_cos"
]

X_train = train[feature_cols]
y_train = train["pv_power_kw"]
X_test = test[feature_cols]
y_test = test["pv_power_kw"]

# Random Forest
rf = RandomForestRegressor(n_estimators=100, random_state=42, n_jobs=-1)
rf.fit(X_train, y_train)
rf_pred = rf.predict(X_test).clip(min=0)

rf_mae = mean_absolute_error(y_test, rf_pred)
rf_rmse = np.sqrt(mean_squared_error(y_test, rf_pred))

What does this code do?

  • feature_cols lists all 12 features: raw weather variables, solar geometry, clear-sky index, and cyclical time encodings.
  • X_train and X_test extract only those columns for the models, keeping temporal order intact.
  • The Random Forest uses 100 decision trees (n_estimators=100) that each vote on the prediction. .clip(min=0) ensures no negative power predictions.
  • In the same cell, add the XGBoost regressor below the Random Forest code:
# XGBoost
xgb_model = xgb.XGBRegressor(
    n_estimators=200,
    max_depth=6,
    learning_rate=0.1,
    random_state=42,
    n_jobs=-1
)
xgb_model.fit(X_train, y_train)
xgb_pred = xgb_model.predict(X_test).clip(min=0)

xgb_mae = mean_absolute_error(y_test, xgb_pred)
xgb_rmse = np.sqrt(mean_squared_error(y_test, xgb_pred))

print(f"Random Forest:  MAE={rf_mae:.3f} kW, RMSE={rf_rmse:.3f} kW")
print(f"XGBoost:        MAE={xgb_mae:.3f} kW, RMSE={xgb_rmse:.3f} kW")

What does this code do?

  • XGBoost uses gradient boosting: each of its 200 trees corrects the errors of the previous tree, learning progressively.
  • max_depth=6 limits tree depth to prevent overfitting. learning_rate=0.1 controls how aggressively each tree corrects errors.
  • The print statements show both models' MAE and RMSE so you can immediately compare them to SARIMAX.
  • Run the cell.

You should see both models' MAE printed, with XGBoost showing a lower error than Random Forest, and both significantly lower than SARIMAX.

Seeing a MemoryError or slow training?

  • If your machine is low on RAM, change n_jobs=-1 to n_jobs=1 for both models. This uses a single CPU core and reduces memory usage.
  • Ensure you ran the feature engineering cell (Cell 5) and the train/test split cell (Cell 6) before this cell.
  • help me debug my tree-based model training

Why do tree-based models outperform SARIMAX here?

SARIMAX models linear relationships between inputs and output. PV production has strongly non-linear dependencies. For example, power drops sharply when cloud cover crosses a threshold, but the exact threshold depends on solar angle.

Decision trees capture these interaction effects automatically through their branching structure. They don't need you to specify which features interact.

Build the final 4-model comparison

Individual MAE numbers are hard to interpret in isolation. A skill score normalizes each model's error against persistence, giving you a single number that says "how much better is this than the naive approach."

  • Create a new cell and add this comparison code:
# Cell 11: Final Model Comparison
# Skill score: 1 - (model_MAE / persistence_MAE)
results_df = pd.DataFrame({
    "Model": ["Persistence", "SARIMAX", "Random Forest", "XGBoost"],
    "MAE (kW)": [persistence_mae, sarimax_mae, rf_mae, xgb_mae],
    "RMSE (kW)": [persistence_rmse, sarimax_rmse, rf_rmse, xgb_rmse],
    "Skill Score": [
        0,
        1 - sarimax_mae / persistence_mae,
        1 - rf_mae / persistence_mae,
        1 - xgb_mae / persistence_mae
    ]
})

print("\n=== MODEL COMPARISON ===")
print(results_df.to_string(index=False))

What is a skill score?

Skill score is defined as 1 - (model_MAE / persistence_MAE). A score of 0 means "no better than persistence." A score of 0.7 means the model removes 70% of the persistence error. This is the standard metric in solar forecasting research.

  • In the same cell, add the bar chart visualization below the print statements:
# Bar chart comparison
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

colors = ["gray", "steelblue", "forestgreen", "darkorange"]
axes[0].bar(results_df["Model"], results_df["MAE (kW)"], color=colors)
axes[0].set_title("MAE Comparison (lower is better)")
axes[0].set_ylabel("MAE (kW)")

axes[1].bar(results_df["Model"], results_df["Skill Score"], color=colors)
axes[1].set_title("Skill Score vs Persistence (higher is better)")
axes[1].set_ylabel("Skill Score")
axes[1].axhline(y=0, color="black", linestyle="-", linewidth=0.5)

plt.tight_layout()
plt.show()

What does this chart show?

  • The left panel shows raw MAE in kW. Lower bars mean more accurate forecasts.
  • The right panel shows skill score. Persistence is the zero baseline. Every bar above zero represents improvement over the naive approach.
  • XGBoost (orange) should have the tallest skill score bar and the shortest MAE bar.
  • Run the cell.

You should see a printed comparison table followed by a side-by-side bar chart. XGBoost achieves the lowest MAE and the highest skill score of all four models.

Seeing a NameError for one of the metric variables?

  • If you see a NameError for sarimax_mae or rf_mae, re-run the cells that computed them (Cell 9 for SARIMAX, Cell 10 for tree models).
  • Jupyter notebooks lose variables if you restart the kernel. If you restarted, re-run all cells from Cell 1 onward.
  • help me fix my model comparison cell
Visualize XGBoost feature importance

Knowing which features matter most tells you where to focus data quality efforts in production. If clear-sky index dominates, then accurate clear-sky modeling is critical to your forecast quality.

  • Create a new cell and add this feature importance visualization:
# Cell 12: XGBoost Feature Importance
importance = pd.Series(
    xgb_model.feature_importances_,
    index=feature_cols
).sort_values(ascending=True)

fig, ax = plt.subplots(figsize=(8, 5))
importance.plot(kind="barh", ax=ax, color="darkorange")
ax.set_title("XGBoost Feature Importance for PV Forecasting")
ax.set_xlabel("Importance")
plt.tight_layout()
plt.show()

What does this code do?

  • feature_importances_ extracts how much each feature contributed to reducing prediction error across all 200 trees.
  • Sorting ascending puts the most important features at the top of the horizontal bar chart, making it easy to read.

Before you run this cell, predict: which feature do you think XGBoost relies on most for PV forecasting?

  • Run the cell.

You should see a horizontal bar chart with clear-sky index and GHI as the top two predictors. This makes physical sense. Clear-sky index isolates the weather effect from predictable solar geometry, giving the model the most informative single signal about cloud impact.

Feature importance chart looks empty?

  • Make sure you ran Cell 10 (the tree-based models cell) first. The xgb_model variable must exist before you can access its feature_importances_ attribute.
  • Check that feature_cols is defined in Cell 10. If you split the cells differently, this variable might be missing.
  • help me debug my feature importance plot

✔️ Awesome, I've got everything!

Great. Double check you've saved your notebook before continuing.

ⓧ I'd like to double check the full code

Your notebook should have Cells 10, 11, and 12 added in this step. Scroll up and confirm each is present in order:

  • Cell 10: defines feature_cols, creates train/test arrays, trains Random Forest and XGBoost, prints both MAE/RMSE values.
  • Cell 11: builds results_df with all 4 models, prints the comparison table, and displays the side-by-side bar chart.
  • Cell 12: extracts xgb_model.feature_importances_ and displays the horizontal bar chart.

You now have quantitative proof that gradient boosting dominates PV forecasting. Your XGBoost model, powered by the clear-sky index feature you engineered, delivers the most accurate predictions. Next up in the Secret Mission, you'll put this model to work generating a live 48-hour forecast and identifying grid injection surplus hours.

Secret mission

48-Hour Forward Forecast and Grid Injection Analysis

Your models can predict historical PV output, but a grid operator needs to know what happens tomorrow. In this extension, you fetch live weather forecast data, run your best model on it, and compute an operational metric: how many hours will your rooftop system inject surplus power into the grid?

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally and uses the free Open-Meteo API. There are no ongoing costs or cloud resources to worry about.

Decide whether to keep your resources running, pause them to come back later, or delete them entirely.

Resources you used:

  • Jupyter notebook server (running process in your terminal).
  • Notebook file: pv_forecasting_lyon.ipynb.
  • Python packages: pvlib, scikit-learn, xgboost, statsmodels, matplotlib, pandas, numpy, requests, jupyter.

Keep everything running

No action needed. Choose this if you're still actively building or plan to extend the notebook with new models or locations.

  • Your Jupyter server stays running in the terminal where you launched it.
  • Your notebook and all installed packages remain available for immediate use.
  • The Open-Meteo API requires no key or account, so there is nothing to maintain on that side.

Pause - I'll come back to this later

Shut down the Jupyter server to free up the terminal and memory. Your notebook file and packages stay on disk for next time.

  • Switch to the terminal where Jupyter is running and press Ctrl+C to stop the server.
  • Confirm the shutdown when prompted.
  • To resume later, navigate to the folder containing your notebook and run jupyter notebook again.

Delete - I don't want to use this again

Remove the notebook file and uninstall the project packages to start fresh.

  • Stop the Jupyter server by pressing Ctrl+C in the terminal where it is running.
  • Delete the notebook file by running:
rm ~/Desktop/pv_forecasting_lyon.ipynb

On Windows, run:

del ~\Desktop\pv_forecasting_lyon.ipynb
  • Uninstall the project packages by running:
pip uninstall pvlib scikit-learn xgboost statsmodels matplotlib pandas numpy requests jupyter -y
  • Confirm the packages are removed by running pip list and checking that pvlib, xgboost, and statsmodels no longer appear.

Nice Work!

Nice Work!

You just built a production-grade solar PV forecasting pipeline that takes raw weather data all the way to hourly power predictions for a 5 kWc rooftop system in Lyon.

You've learned how to:

  • Build a complete PV forecasting pipeline from Open-Meteo weather data through pvlib simulation to ML-powered predictions for a real rooftop installation.
  • Engineer domain-specific features like clear-sky index, solar zenith angle, and cyclical time encodings that transform raw weather into high-signal inputs for energy models.
  • Implement and compare SARIMAX, Random Forest, and XGBoost against a persistence baseline, proving quantitatively why gradient boosting leads in energy forecasting.
  • Secret Mission: Produce an operational 48-hour forward forecast with grid injection surplus analysis, the exact output a grid operator needs for dispatch decisions.

Ready to quiz yourself?