Analyse de production photovoltaïque avec pvlib¶

Ce notebook présente l'utilisation de la bibliothèque pvlib pour modéliser et analyser des systèmes photovoltaïques.

Introduction à pvlib¶

pvlib est une bibliothèque Python qui fournit un ensemble d'outils pour la modélisation des performances des systèmes photovoltaïques. Elle permet de simuler la production d'énergie en fonction de différents paramètres comme l'emplacement géographique, l'orientation des panneaux, les conditions météorologiques, etc.

In [1]:

# Importation des bibliothèques nécessaires
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pvlib
from pvlib import location
from pvlib import irradiance
from pvlib import atmosphere
from pvlib import pvsystem
from pvlib import temperature
from pvlib import tracking
from pvlib import modelchain
from datetime import datetime, timedelta
import pytz

# Configuration de l'affichage
plt.style.use('ggplot')
%matplotlib inline

# Importation des bibliothèques nécessaires import numpy as np import pandas as pd import matplotlib.pyplot as plt import pvlib from pvlib import location from pvlib import irradiance from pvlib import atmosphere from pvlib import pvsystem from pvlib import temperature from pvlib import tracking from pvlib import modelchain from datetime import datetime, timedelta import pytz # Configuration de l'affichage plt.style.use('ggplot') %matplotlib inline

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 5
      3 import pandas as pd
      4 import matplotlib.pyplot as plt
----> 5 import pvlib
      6 from pvlib import location
      7 from pvlib import irradiance

ModuleNotFoundError: No module named 'pvlib'

Définition de l'emplacement¶

Nous allons définir un emplacement en France pour notre système photovoltaïque.

In [ ]:

# Définition de l'emplacement (Paris, France)
latitude = 48.8566
longitude = 2.3522
altitude = 35  # mètres
tz = 'Europe/Paris'

# Création de l'objet Location
site = pvlib.location.Location(latitude, longitude, tz, altitude)

print(f"Emplacement: {site.name}")
print(f"Latitude: {site.latitude:.4f}°")
print(f"Longitude: {site.longitude:.4f}°")
print(f"Altitude: {site.altitude} m")
print(f"Fuseau horaire: {site.tz}")

# Définition de l'emplacement (Paris, France) latitude = 48.8566 longitude = 2.3522 altitude = 35 # mètres tz = 'Europe/Paris' # Création de l'objet Location site = pvlib.location.Location(latitude, longitude, tz, altitude) print(f"Emplacement: {site.name}") print(f"Latitude: {site.latitude:.4f}°") print(f"Longitude: {site.longitude:.4f}°") print(f"Altitude: {site.altitude} m") print(f"Fuseau horaire: {site.tz}")

Calcul de la position du soleil¶

Nous allons calculer la position du soleil pour une journée donnée.

In [ ]:

# Définition de la période (une journée d'été)
start = pd.Timestamp('2023-06-21 00:00:00', tz=tz)
end = pd.Timestamp('2023-06-21 23:59:59', tz=tz)
times = pd.date_range(start=start, end=end, freq='10min')

# Calcul de la position solaire
solpos = site.get_solarposition(times)

# Affichage des résultats
plt.figure(figsize=(12, 6))
plt.plot(times, solpos['elevation'], label='Élévation')
plt.plot(times, solpos['azimuth'], label='Azimut')
plt.xlabel('Heure')
plt.ylabel('Angle (degrés)')
plt.title('Position du soleil le 21 juin 2023 à Paris')
plt.legend()
plt.grid(True)
plt.show()

# Affichage des heures de lever et coucher du soleil
sunrise = solpos[solpos['elevation'] > 0].index[0]
sunset = solpos[solpos['elevation'] > 0].index[-1]
print(f"Lever du soleil: {sunrise.strftime('%H:%M')}")
print(f"Coucher du soleil: {sunset.strftime('%H:%M')}")
print(f"Durée du jour: {sunset - sunrise}")

# Définition de la période (une journée d'été) start = pd.Timestamp('2023-06-21 00:00:00', tz=tz) end = pd.Timestamp('2023-06-21 23:59:59', tz=tz) times = pd.date_range(start=start, end=end, freq='10min') # Calcul de la position solaire solpos = site.get_solarposition(times) # Affichage des résultats plt.figure(figsize=(12, 6)) plt.plot(times, solpos['elevation'], label='Élévation') plt.plot(times, solpos['azimuth'], label='Azimut') plt.xlabel('Heure') plt.ylabel('Angle (degrés)') plt.title('Position du soleil le 21 juin 2023 à Paris') plt.legend() plt.grid(True) plt.show() # Affichage des heures de lever et coucher du soleil sunrise = solpos[solpos['elevation'] > 0].index[0] sunset = solpos[solpos['elevation'] > 0].index[-1] print(f"Lever du soleil: {sunrise.strftime('%H:%M')}") print(f"Coucher du soleil: {sunset.strftime('%H:%M')}") print(f"Durée du jour: {sunset - sunrise}")

Calcul de l'irradiance solaire¶

Nous allons calculer l'irradiance solaire pour notre emplacement.

In [ ]:

# Calcul de l'irradiance extraterrestre
irradiance_et = pvlib.irradiance.get_extra_radiation(times)

# Calcul de l'airmass
airmass = pvlib.atmosphere.get_relative_airmass(solpos['apparent_zenith'])
airmass_abs = pvlib.atmosphere.get_absolute_airmass(airmass)

# Calcul de l'irradiance par ciel clair (modèle INEICHEN)
irradiance_clear = pvlib.irradiance.ineichen(
    solpos['apparent_zenith'], 
    airmass_abs, 
    altitude=altitude,
    dni_extra=irradiance_et
)

# Affichage des résultats
plt.figure(figsize=(12, 6))
plt.plot(times, irradiance_et, label='Irradiance extraterrestre')
plt.plot(times, irradiance_clear['ghi'], label='GHI (ciel clair)')
plt.plot(times, irradiance_clear['dni'], label='DNI (ciel clair)')
plt.plot(times, irradiance_clear['dhi'], label='DHI (ciel clair)')
plt.xlabel('Heure')
plt.ylabel('Irradiance (W/m²)')
plt.title('Irradiance solaire le 21 juin 2023 à Paris')
plt.legend()
plt.grid(True)
plt.show()

# Affichage des valeurs maximales
print(f"Irradiance extraterrestre maximale: {irradiance_et.max():.2f} W/m²")
print(f"GHI maximale (ciel clair): {irradiance_clear['ghi'].max():.2f} W/m²")
print(f"DNI maximale (ciel clair): {irradiance_clear['dni'].max():.2f} W/m²")
print(f"DHI maximale (ciel clair): {irradiance_clear['dhi'].max():.2f} W/m²")

# Calcul de l'irradiance extraterrestre irradiance_et = pvlib.irradiance.get_extra_radiation(times) # Calcul de l'airmass airmass = pvlib.atmosphere.get_relative_airmass(solpos['apparent_zenith']) airmass_abs = pvlib.atmosphere.get_absolute_airmass(airmass) # Calcul de l'irradiance par ciel clair (modèle INEICHEN) irradiance_clear = pvlib.irradiance.ineichen( solpos['apparent_zenith'], airmass_abs, altitude=altitude, dni_extra=irradiance_et ) # Affichage des résultats plt.figure(figsize=(12, 6)) plt.plot(times, irradiance_et, label='Irradiance extraterrestre') plt.plot(times, irradiance_clear['ghi'], label='GHI (ciel clair)') plt.plot(times, irradiance_clear['dni'], label='DNI (ciel clair)') plt.plot(times, irradiance_clear['dhi'], label='DHI (ciel clair)') plt.xlabel('Heure') plt.ylabel('Irradiance (W/m²)') plt.title('Irradiance solaire le 21 juin 2023 à Paris') plt.legend() plt.grid(True) plt.show() # Affichage des valeurs maximales print(f"Irradiance extraterrestre maximale: {irradiance_et.max():.2f} W/m²") print(f"GHI maximale (ciel clair): {irradiance_clear['ghi'].max():.2f} W/m²") print(f"DNI maximale (ciel clair): {irradiance_clear['dni'].max():.2f} W/m²") print(f"DHI maximale (ciel clair): {irradiance_clear['dhi'].max():.2f} W/m²")

Modélisation d'un système photovoltaïque¶

Nous allons maintenant modéliser un système photovoltaïque simple et calculer sa production.

In [ ]:

# Définition des paramètres du système PV
surface_tilt = 30  # inclinaison des panneaux en degrés
surface_azimuth = 180  # orientation des panneaux (sud = 180°)
albedo = 0.2  # albédo du sol

# Calcul de l'irradiance sur le plan incliné
poa_irradiance = pvlib.irradiance.get_total_irradiance(
    surface_tilt=surface_tilt,
    surface_azimuth=surface_azimuth,
    dni=irradiance_clear['dni'],
    ghi=irradiance_clear['ghi'],
    dhi=irradiance_clear['dhi'],
    solar_zenith=solpos['apparent_zenith'],
    solar_azimuth=solpos['azimuth'],
    model='haydavies',
    albedo=albedo
)

# Affichage des résultats
plt.figure(figsize=(12, 6))
plt.plot(times, poa_irradiance['poa_global'], label='Irradiance globale sur plan incliné')
plt.plot(times, irradiance_clear['ghi'], label='GHI (plan horizontal)', linestyle='--')
plt.xlabel('Heure')
plt.ylabel('Irradiance (W/m²)')
plt.title(f'Irradiance sur un plan incliné à {surface_tilt}° orienté au sud')
plt.legend()
plt.grid(True)
plt.show()

# Calcul du gain d'irradiance
ghi_daily = irradiance_clear['ghi'].sum() * 10/60  # kWh/m²/jour (10 min = 1/6 heure)
poa_daily = poa_irradiance['poa_global'].sum() * 10/60  # kWh/m²/jour
gain = (poa_daily / ghi_daily - 1) * 100

print(f"Irradiance journalière sur plan horizontal: {ghi_daily:.2f} kWh/m²/jour")
print(f"Irradiance journalière sur plan incliné: {poa_daily:.2f} kWh/m²/jour")
print(f"Gain d'irradiance: {gain:.1f}%")

# Définition des paramètres du système PV surface_tilt = 30 # inclinaison des panneaux en degrés surface_azimuth = 180 # orientation des panneaux (sud = 180°) albedo = 0.2 # albédo du sol # Calcul de l'irradiance sur le plan incliné poa_irradiance = pvlib.irradiance.get_total_irradiance( surface_tilt=surface_tilt, surface_azimuth=surface_azimuth, dni=irradiance_clear['dni'], ghi=irradiance_clear['ghi'], dhi=irradiance_clear['dhi'], solar_zenith=solpos['apparent_zenith'], solar_azimuth=solpos['azimuth'], model='haydavies', albedo=albedo ) # Affichage des résultats plt.figure(figsize=(12, 6)) plt.plot(times, poa_irradiance['poa_global'], label='Irradiance globale sur plan incliné') plt.plot(times, irradiance_clear['ghi'], label='GHI (plan horizontal)', linestyle='--') plt.xlabel('Heure') plt.ylabel('Irradiance (W/m²)') plt.title(f'Irradiance sur un plan incliné à {surface_tilt}° orienté au sud') plt.legend() plt.grid(True) plt.show() # Calcul du gain d'irradiance ghi_daily = irradiance_clear['ghi'].sum() * 10/60 # kWh/m²/jour (10 min = 1/6 heure) poa_daily = poa_irradiance['poa_global'].sum() * 10/60 # kWh/m²/jour gain = (poa_daily / ghi_daily - 1) * 100 print(f"Irradiance journalière sur plan horizontal: {ghi_daily:.2f} kWh/m²/jour") print(f"Irradiance journalière sur plan incliné: {poa_daily:.2f} kWh/m²/jour") print(f"Gain d'irradiance: {gain:.1f}%")

Modélisation d'un module photovoltaïque¶

Nous allons utiliser un module photovoltaïque standard de la base de données CEC pour calculer la production d'électricité.

In [ ]:

# Récupération des données du module PV depuis la base CEC
sandia_modules = pvlib.pvsystem.retrieve_sam('CECMod')
module_name = 'Canadian_Solar_Inc__CS6X_300M'
module = sandia_modules[module_name]

# Affichage des caractéristiques du module
print(f"Module: {module_name}")
print(f"Puissance nominale: {module['STC']} W")

# Affichage des paramètres clés du module
print(f"Technologie: {module['Technology']}")
print(f"Surface: {module['A_c']:.2f} m²")
print(f"Coefficient de température: {module['PTC'] / module['STC']:.3f} %/°C")

# Récupération des données du module PV depuis la base CEC sandia_modules = pvlib.pvsystem.retrieve_sam('CECMod') module_name = 'Canadian_Solar_Inc__CS6X_300M' module = sandia_modules[module_name] # Affichage des caractéristiques du module print(f"Module: {module_name}") print(f"Puissance nominale: {module['STC']} W") # Affichage des paramètres clés du module print(f"Technologie: {module['Technology']}") print(f"Surface: {module['A_c']:.2f} m²") print(f"Coefficient de température: {module['PTC'] / module['STC']:.3f} %/°C")

In [ ]:

# Calcul de la température du module
temp_air = pd.Series(25, index=times)  # température ambiante constante de 25°C
wind_speed = pd.Series(1, index=times)  # vitesse du vent constante de 1 m/s

# Modèle SAPM pour la température du module
temp_model_params = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS['sapm']['open_rack_glass_glass']
temp_cell = pvlib.temperature.sapm_cell(poa_irradiance['poa_global'], temp_air, wind_speed, **temp_model_params)

# Calcul de la puissance DC
effective_irradiance = poa_irradiance['poa_global']
dc_parameters = {
    'pdc0': module['STC'],
    'gamma_pdc': module['PTC'] / module['STC'] - 1,  # approximation du coefficient de température
    'temp_ref': 25
}

# Utilisation du modèle PVWatts pour calculer la puissance DC
dc_power = pvlib.pvsystem.pvwatts_dc(effective_irradiance, temp_cell, **dc_parameters)

# Calcul de la puissance AC (avec un onduleur de rendement 96%)
ac_power = dc_power * 0.96

# Affichage des résultats
plt.figure(figsize=(12, 6))
plt.plot(times, dc_power, label='Puissance DC')
plt.plot(times, ac_power, label='Puissance AC')
plt.xlabel('Heure')
plt.ylabel('Puissance (W)')
plt.title(f'Production d\'un module {module_name} de {module["STC"]} W')
plt.legend()
plt.grid(True)
plt.show()

# Calcul de la production journalière
dc_energy = dc_power.sum() * 10/60 / 1000  # kWh/jour
ac_energy = ac_power.sum() * 10/60 / 1000  # kWh/jour

print(f"Production DC journalière: {dc_energy:.2f} kWh/jour")
print(f"Production AC journalière: {ac_energy:.2f} kWh/jour")
print(f"Rendement global: {ac_energy/(poa_daily * module['A_c']):.1%}")

# Calcul de la température du module temp_air = pd.Series(25, index=times) # température ambiante constante de 25°C wind_speed = pd.Series(1, index=times) # vitesse du vent constante de 1 m/s # Modèle SAPM pour la température du module temp_model_params = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS['sapm']['open_rack_glass_glass'] temp_cell = pvlib.temperature.sapm_cell(poa_irradiance['poa_global'], temp_air, wind_speed, **temp_model_params) # Calcul de la puissance DC effective_irradiance = poa_irradiance['poa_global'] dc_parameters = { 'pdc0': module['STC'], 'gamma_pdc': module['PTC'] / module['STC'] - 1, # approximation du coefficient de température 'temp_ref': 25 } # Utilisation du modèle PVWatts pour calculer la puissance DC dc_power = pvlib.pvsystem.pvwatts_dc(effective_irradiance, temp_cell, **dc_parameters) # Calcul de la puissance AC (avec un onduleur de rendement 96%) ac_power = dc_power * 0.96 # Affichage des résultats plt.figure(figsize=(12, 6)) plt.plot(times, dc_power, label='Puissance DC') plt.plot(times, ac_power, label='Puissance AC') plt.xlabel('Heure') plt.ylabel('Puissance (W)') plt.title(f'Production d\'un module {module_name} de {module["STC"]} W') plt.legend() plt.grid(True) plt.show() # Calcul de la production journalière dc_energy = dc_power.sum() * 10/60 / 1000 # kWh/jour ac_energy = ac_power.sum() * 10/60 / 1000 # kWh/jour print(f"Production DC journalière: {dc_energy:.2f} kWh/jour") print(f"Production AC journalière: {ac_energy:.2f} kWh/jour") print(f"Rendement global: {ac_energy/(poa_daily * module['A_c']):.1%}")

Analyse de sensibilité¶

Étudions l'impact de l'inclinaison et de l'orientation des panneaux sur la production.

In [ ]:

# Fonction pour calculer la production en fonction de l'inclinaison et de l'orientation
def calculate_energy(tilt, azimuth):
    # Calcul de l'irradiance sur le plan incliné
    poa = pvlib.irradiance.get_total_irradiance(
        surface_tilt=tilt,
        surface_azimuth=azimuth,
        dni=irradiance_clear['dni'],
        ghi=irradiance_clear['ghi'],
        dhi=irradiance_clear['dhi'],
        solar_zenith=solpos['apparent_zenith'],
        solar_azimuth=solpos['azimuth'],
        model='haydavies',
        albedo=albedo
    )
    
    # Calcul de la température du module
    temp_cell = pvlib.temperature.sapm_cell(poa['poa_global'], temp_air, wind_speed, **temp_model_params)
    
    # Calcul de la puissance DC
    dc_power = pvlib.pvsystem.pvwatts_dc(poa['poa_global'], temp_cell, **dc_parameters)
    
    # Calcul de la puissance AC
    ac_power = dc_power * 0.96
    
    # Calcul de l'énergie journalière
    ac_energy = ac_power.sum() * 10/60 / 1000  # kWh/jour
    
    return ac_energy

# Création d'une grille d'inclinaisons et d'orientations
tilts = np.arange(0, 91, 10)
azimuths = np.arange(90, 271, 10)  # 90=Est, 180=Sud, 270=Ouest

# Calcul de la production pour chaque combinaison
energy_matrix = np.zeros((len(tilts), len(azimuths)))

for i, tilt in enumerate(tilts):
    for j, azimuth in enumerate(azimuths):
        energy_matrix[i, j] = calculate_energy(tilt, azimuth)

# Normalisation par rapport à la valeur maximale
energy_matrix_norm = energy_matrix / np.max(energy_matrix)

# Affichage des résultats
plt.figure(figsize=(12, 8))
contour = plt.contourf(azimuths, tilts, energy_matrix_norm, 20, cmap='viridis')
plt.colorbar(contour, label='Production relative')
plt.xlabel('Azimut (degrés)')
plt.ylabel('Inclinaison (degrés)')
plt.title('Production relative en fonction de l\'inclinaison et de l\'orientation')
plt.grid(True)

# Marquer le point optimal
max_idx = np.unravel_index(np.argmax(energy_matrix), energy_matrix.shape)
opt_tilt = tilts[max_idx[0]]
opt_azimuth = azimuths[max_idx[1]]
plt.plot(opt_azimuth, opt_tilt, 'ro', markersize=10)
plt.annotate(f'Optimal: {opt_tilt}°, {opt_azimuth}°', 
             xy=(opt_azimuth, opt_tilt), xytext=(opt_azimuth+20, opt_tilt+10),
             arrowprops=dict(facecolor='black', shrink=0.05))

plt.show()

print(f"Inclinaison optimale: {opt_tilt}°")
print(f"Orientation optimale: {opt_azimuth}° (180° = Sud)")
print(f"Production maximale: {energy_matrix[max_idx]:.2f} kWh/jour")

# Fonction pour calculer la production en fonction de l'inclinaison et de l'orientation def calculate_energy(tilt, azimuth): # Calcul de l'irradiance sur le plan incliné poa = pvlib.irradiance.get_total_irradiance( surface_tilt=tilt, surface_azimuth=azimuth, dni=irradiance_clear['dni'], ghi=irradiance_clear['ghi'], dhi=irradiance_clear['dhi'], solar_zenith=solpos['apparent_zenith'], solar_azimuth=solpos['azimuth'], model='haydavies', albedo=albedo ) # Calcul de la température du module temp_cell = pvlib.temperature.sapm_cell(poa['poa_global'], temp_air, wind_speed, **temp_model_params) # Calcul de la puissance DC dc_power = pvlib.pvsystem.pvwatts_dc(poa['poa_global'], temp_cell, **dc_parameters) # Calcul de la puissance AC ac_power = dc_power * 0.96 # Calcul de l'énergie journalière ac_energy = ac_power.sum() * 10/60 / 1000 # kWh/jour return ac_energy # Création d'une grille d'inclinaisons et d'orientations tilts = np.arange(0, 91, 10) azimuths = np.arange(90, 271, 10) # 90=Est, 180=Sud, 270=Ouest # Calcul de la production pour chaque combinaison energy_matrix = np.zeros((len(tilts), len(azimuths))) for i, tilt in enumerate(tilts): for j, azimuth in enumerate(azimuths): energy_matrix[i, j] = calculate_energy(tilt, azimuth) # Normalisation par rapport à la valeur maximale energy_matrix_norm = energy_matrix / np.max(energy_matrix) # Affichage des résultats plt.figure(figsize=(12, 8)) contour = plt.contourf(azimuths, tilts, energy_matrix_norm, 20, cmap='viridis') plt.colorbar(contour, label='Production relative') plt.xlabel('Azimut (degrés)') plt.ylabel('Inclinaison (degrés)') plt.title('Production relative en fonction de l\'inclinaison et de l\'orientation') plt.grid(True) # Marquer le point optimal max_idx = np.unravel_index(np.argmax(energy_matrix), energy_matrix.shape) opt_tilt = tilts[max_idx[0]] opt_azimuth = azimuths[max_idx[1]] plt.plot(opt_azimuth, opt_tilt, 'ro', markersize=10) plt.annotate(f'Optimal: {opt_tilt}°, {opt_azimuth}°', xy=(opt_azimuth, opt_tilt), xytext=(opt_azimuth+20, opt_tilt+10), arrowprops=dict(facecolor='black', shrink=0.05)) plt.show() print(f"Inclinaison optimale: {opt_tilt}°") print(f"Orientation optimale: {opt_azimuth}° (180° = Sud)") print(f"Production maximale: {energy_matrix[max_idx]:.2f} kWh/jour")

Conclusion¶

Ce notebook a présenté l'utilisation de pvlib pour modéliser et analyser un système photovoltaïque. Nous avons vu comment :

1. Calculer la position du soleil pour un emplacement donné
2. Estimer l'irradiance solaire par ciel clair
3. Calculer l'irradiance sur un plan incliné
4. Modéliser la production d'un module photovoltaïque
5. Analyser l'impact de l'inclinaison et de l'orientation sur la production

Ces outils permettent d'optimiser la conception des systèmes photovoltaïques et d'estimer leur production d'énergie.