"""
Photovoltaics Module
====================
Photovoltaic power generation modeling based on EU PVGIS methodology.
:author: Dipl.-Ing. (FH) Jonas Pfeiffer
"""
import numpy as np
import pandas as pd
from districtheatingsim.heat_generators.solar_radiation import calculate_solar_radiation
from districtheatingsim.utilities.test_reference_year import import_TRY
# Constant for degree-radian conversion
DEG_TO_RAD = np.pi / 180
[docs]
def Calculate_PV(
TRY_data: str,
Gross_area: float,
Longitude: float,
STD_Longitude: float,
Latitude: float,
Albedo: float,
East_West_collector_azimuth_angle: float,
Collector_tilt_angle: float,
) -> tuple[float, float, np.ndarray]:
"""
Calculate photovoltaic power output based on EU PVGIS methodology.
:param TRY_data: Path to Test Reference Year data
:type TRY_data: str
:param Gross_area: PV system area [m²]
:type Gross_area: float
:param Longitude: Geographic longitude [degrees]
:type Longitude: float
:param STD_Longitude: Standard longitude for time zone [degrees]
:type STD_Longitude: float
:param Latitude: Geographic latitude [degrees]
:type Latitude: float
:param Albedo: Ground reflection coefficient [-]
:type Albedo: float
:param East_West_collector_azimuth_angle: Azimuth angle [degrees]
:type East_West_collector_azimuth_angle: float
:param Collector_tilt_angle: Tilt angle from horizontal [degrees]
:type Collector_tilt_angle: float
:return: (yield_kWh, P_max, P_L)
:rtype: tuple
"""
# Import Test Reference Year meteorological data
Ta_L, W_L, D_L, G_L, _ = import_TRY(TRY_data)
# Define photovoltaic system constants based on EU PVGIS methodology
eff_nom = 0.199 # Nominal PV module efficiency under STC conditions [-]
sys_loss = 0.14 # Total system losses including inverter, cabling, soiling [-]
U0 = 26.9 # Thermal loss coefficient constant [W/(°C·m²)]
U1 = 6.2 # Thermal loss coefficient wind-dependent [W·s/(°C·m³)]
# EU PVGIS efficiency model coefficients for crystalline silicon
k1, k2, k3, k4, k5, k6 = -0.017237, -0.040465, -0.004702, 0.000149, 0.000170, 0.000005
# Generate annual hourly time series (8760 hours)
start_date = np.datetime64("2024-01-01T00:00")
time_steps = start_date + np.arange(8760) * np.timedelta64(1, "h")
# Calculate solar irradiation on tilted PV surface
GT_L, _, _, _ = calculate_solar_radiation(
time_steps,
G_L,
D_L,
Longitude,
STD_Longitude,
Latitude,
Albedo,
East_West_collector_azimuth_angle,
Collector_tilt_angle,
)
# Convert irradiation from W/m² to kW/m² for efficiency calculations
G1 = GT_L / 1000
# Calculate PV module temperature using thermal model
# Accounts for ambient temperature, solar heating, and wind cooling
Tm = Ta_L + GT_L / (U0 + U1 * W_L)
T1m = Tm - 25 # Temperature difference from Standard Test Conditions
# Calculate relative efficiency based on irradiation and temperature
eff_rel = np.ones_like(G1)
non_zero_mask = G1 != 0
# Apply EU PVGIS efficiency model for non-zero irradiation periods
eff_rel[non_zero_mask] = (
1
+ k1 * np.log(G1[non_zero_mask])
+ k2 * np.log(G1[non_zero_mask]) ** 2
+ k3 * T1m[non_zero_mask]
+ k4 * T1m[non_zero_mask] * np.log(G1[non_zero_mask])
+ k5 * Tm[non_zero_mask] * np.log(G1[non_zero_mask]) ** 2
+ k6 * Tm[non_zero_mask] ** 2
)
# Set efficiency to zero for no irradiation periods
eff_rel[~non_zero_mask] = 0
eff_rel = np.nan_to_num(eff_rel, nan=0)
# Calculate instantaneous PV power output [kW]
# P = Irradiation × Area × Nominal_Efficiency × Relative_Efficiency × (1 - System_Losses)
P_L = G1 * Gross_area * eff_nom * eff_rel * (1 - sys_loss)
# Calculate performance metrics
P_max = np.max(P_L) # Maximum instantaneous power [kW]
E = np.sum(P_L) # Total annual energy [kWh]
# Round results for practical use
yield_kWh = round(E, 2)
P_max = round(P_max, 2)
return yield_kWh, P_max, P_L
[docs]
def azimuth_angle(direction: str) -> float | None:
"""
Convert cardinal direction to azimuth angle.
:param direction: Cardinal direction (N, S, E, W, O, NE, etc.)
:type direction: str
:return: Azimuth angle [degrees] or None
:rtype: float or None
"""
# Azimuth angle mapping for PV system orientations
# Following solar energy conventions: 0° = South, 90° = West
azimuths = {
"N": 180, # North - minimal irradiation
"W": 90, # West - afternoon energy production
"S": 0, # South - optimal for Northern Hemisphere
"O": 270, # East ("Ost" in German) - morning energy production
"NO": 225, # Northeast ("Nordost") - morning-biased production
"SO": 315, # Southeast ("Südost") - morning-optimal production
"SW": 45, # Southwest - afternoon-optimal production
"NW": 135, # Northwest - afternoon-biased production
}
return azimuths.get(direction.upper(), None)
[docs]
def calculate_building(TRY_data: str, building_data: str, output_filename: str) -> None:
"""
Calculate photovoltaic yield for multiple buildings.
:param TRY_data: Path to Test Reference Year data
:type TRY_data: str
:param building_data: Path to building CSV (building_id, roof_area, orientation)
:type building_data: str
:param output_filename: Path for results CSV output
:type output_filename: str
"""
# Load building data from CSV file with robust error handling
try:
gdata = np.genfromtxt(building_data, delimiter=";", skip_header=1, dtype=None, encoding="utf-8")
except Exception as e:
raise ValueError(f"Error loading building data from {building_data}: {e}") from e
# Location parameters for PV calculation (customizable)
# Default: Munich, Germany - Central European location
Longitude = 11.576 # Munich longitude [degrees]
STD_Longitude = 15.0 # Central European Time standard longitude [degrees]
Latitude = 48.137 # Munich latitude [degrees]
# PV system parameters
Albedo = 0.2 # Typical ground reflection coefficient [-]
Collector_tilt_angle = 36 # Optimal tilt angle for latitude [degrees]
# Initialize time series for annual hourly data
Annual_hours = np.arange(1, 8761) # Hours 1-8760 for full year
# Initialize results DataFrame with time column
df = pd.DataFrame()
df["Annual Hours"] = Annual_hours
print("Calculating PV yield for buildings...")
print(f"Processing {len(gdata)} building systems...")
# Process each building in the input data
for _idx, (building, area, direction) in enumerate(gdata):
try:
# Convert numpy types to Python native types
building = str(building)
area = float(area)
direction = str(direction)
print(f"Processing {building}: {area:.1f} m² {direction}-facing")
# Get azimuth angle for the specified direction
current_azimuth = azimuth_angle(direction)
# Handle East-West (OW) configuration with split installation
if current_azimuth is None and direction.upper() == "OW":
# Split area equally between East and West orientations
area_split = area / 2
directions = ["O", "W"] # East and West in German notation
print(f" → East-West configuration: {area_split:.1f} m² each direction")
else:
# Single orientation configuration
directions = [direction]
# Calculate PV performance for each orientation
for orientation in directions:
current_azimuth = azimuth_angle(orientation)
if current_azimuth is not None:
# Use split area for EW systems, full area for single orientation
system_area = area_split if direction.upper() == "OW" else area
# Calculate PV performance using EU PVGIS methodology
yield_kWh, max_power, P_L = Calculate_PV(
TRY_data,
system_area,
Longitude,
STD_Longitude,
Latitude,
Albedo,
current_azimuth,
Collector_tilt_angle,
)
# Generate appropriate column suffix for EW systems
suffix = f" {orientation}" if direction.upper() == "OW" else ""
# Display calculation results
print(f" → PV yield {building}{suffix}: {yield_kWh / 1000:.1f} MWh")
print(f" → Maximum power {building}{suffix}: {max_power:.1f} kW")
# Add results to DataFrame with descriptive column name
column_name = f"{building}{suffix} {system_area:.1f} m² [kW]"
df[column_name] = P_L
else:
print(f" → Warning: Invalid orientation '{orientation}' for {building}")
except Exception as e:
print(f" → Error processing {building}: {e}")
continue
# Save comprehensive results to CSV file
try:
df.to_csv(output_filename, index=False, sep=";", encoding="utf-8-sig")
print(f"\nResults successfully saved to: {output_filename}")
# Calculate and display summary statistics
pv_columns = [col for col in df.columns if "[kW]" in col]
total_systems = len(pv_columns)
total_capacity = df[pv_columns].max().sum()
annual_yield = df[pv_columns].sum().sum() / 1000 # Convert to MWh
print("Summary:")
print(f" → Total systems processed: {total_systems}")
print(f" → Total PV capacity: {total_capacity:.1f} kW")
print(f" → Annual district yield: {annual_yield:.1f} MWh")
print(f" → Average capacity factor: {annual_yield * 1000 / (total_capacity * 8760):.2f}")
except Exception as e:
raise OSError(f"Error saving results to {output_filename}: {e}") from e