diff --git a/data_base_influx.py b/data_base_influx.py index ad0adfb..c69ddb8 100644 --- a/data_base_influx.py +++ b/data_base_influx.py @@ -1,5 +1,7 @@ from influxdb_client import InfluxDBClient, Point, WritePrecision from datetime import datetime +import datetime as dt +import pandas as pd class DataBaseInflux: def __init__(self, url: str, token: str, org: str, bucket: str): @@ -25,4 +27,22 @@ class DataBaseInflux: # Punkt in InfluxDB schreiben self.write_api.write(bucket=self.bucket, org=self.org, record=point) + def store_forecasts(self, forecast_name: str, data: pd.Series): + + measurement = forecast_name + run_tag = dt.datetime.now(dt.timezone.utc).replace(second=0, microsecond=0).isoformat(timespec="minutes") + + pts = [] + + series = pd.to_numeric(data, errors="coerce").dropna() + + for ts, val in series.items(): + pts.append( + Point(measurement) + .tag("run", run_tag) + .field("value", float(val)) + .time(ts.to_pydatetime(), WritePrecision.S) + ) + + self.write_api.write(bucket=self.bucket, org=self.org, record=pts) diff --git a/forecaster/weather_forecaster.py b/forecaster/weather_forecaster.py new file mode 100644 index 0000000..a95d5d9 --- /dev/null +++ b/forecaster/weather_forecaster.py @@ -0,0 +1,61 @@ +#!/usr/bin/env python3 +import time +import datetime as dt +import requests +from zoneinfo import ZoneInfo +from matplotlib import pyplot as plt +import pandas as pd + +TZ = "Europe/Berlin" +DAYS = 2 + +OPEN_METEO_URL = "https://api.open-meteo.com/v1/forecast" + +class WeatherForecaster: + def __init__(self, latitude, longitude): + self.lat = latitude + self.lon = longitude + + def get_hourly_forecast(self, start_hour, days): + start_hour_local = start_hour + end_hour_local = start_hour_local + dt.timedelta(days=days) + + params = { + "latitude": self.lat, + "longitude": self.lon, + "hourly": ["temperature_2m", "shortwave_radiation", "wind_speed_10m"], + "timezone": TZ, + "start_hour": start_hour_local.strftime("%Y-%m-%dT%H:%M"), + "end_hour": end_hour_local.strftime("%Y-%m-%dT%H:%M") + } + + h = requests.get(OPEN_METEO_URL, params=params).json()["hourly"] + + time_stamps = h["time"] + time_stamps = [ + dt.datetime.fromisoformat(t).replace(tzinfo=ZoneInfo(TZ)) + for t in time_stamps + ] + + weather = pd.DataFrame(index=time_stamps) + weather["ghi"] = h["shortwave_radiation"] + weather["temp_air"] = h["temperature_2m"] + weather["wind_speed"] = h["wind_speed_10m"] + + return weather + + + +if __name__=='__main__': + + weather_forecast = WeatherForecaster(latitude=48.041, longitude=7.862) + while True: + now = dt.datetime.now() + secs = 60 - now.second #(60 - now.minute) * 60 - now.second # Sekunden bis volle Stunde + time.sleep(secs) + + now_local = dt.datetime.now() + start_hour_local = (now_local + dt.timedelta(hours=1)).replace(minute=0, second=0, microsecond=0) + time_stamps, temps, ghi, wind_speed = weather_forecast.get_hourly_forecast(start_hour_local, DAYS) + plt.plot(time_stamps, temps) + plt.show() \ No newline at end of file diff --git a/main.py b/main.py index e77b395..2edccf5 100644 --- a/main.py +++ b/main.py @@ -1,12 +1,16 @@ import time from datetime import datetime from data_base_influx import DataBaseInflux +from forecaster.weather_forecaster import WeatherForecaster from heat_pump import HeatPump from pv_inverter import PvInverter +from simulators.pv_plant_simulator import PvWattsSubarrayConfig, PvWattsPlant from solaredge_meter import SolaredgeMeter from shelly_pro_3m import ShellyPro3m from energysystem import EnergySystem from sg_ready_controller import SgReadyController +from pvlib.location import Location +import datetime as dt # For dev-System run in terminal: ssh -N -L 127.0.0.1:8111:10.0.0.10:502 pi@192.168.1.146 # For productive-System change IP-adress in heatpump to '10.0.0.10' and port to 502 @@ -31,6 +35,22 @@ meter = SolaredgeMeter(device_name='solaredge_meter', ip_address='192.168.1.112' es.add_components(hp_master, hp_slave, shelly, wr, meter) controller = SgReadyController(es) +# FORECASTING +latitude = 48.041 +longitude = 7.862 +TZ = "Europe/Berlin" +HORIZON_DAYS = 2 +weather_forecaster = WeatherForecaster(latitude=latitude, longitude=longitude) +site = Location(latitude=latitude, longitude=longitude, altitude=35, tz=TZ, name="Gundelfingen") + +p_module = 435 +upper_roof_north = PvWattsSubarrayConfig(name="north", pdc0_w=(29+29+21)*p_module, tilt_deg=10, azimuth_deg=20, dc_loss=0.02, ac_loss=0.01) +upper_roof_south = PvWattsSubarrayConfig(name="south", pdc0_w=(29+21+20)*p_module, tilt_deg=10, azimuth_deg=200, dc_loss=0.02, ac_loss=0.01) +upper_roof_east = PvWattsSubarrayConfig(name="east", pdc0_w=7*p_module, tilt_deg=10, azimuth_deg=110, dc_loss=0.02, ac_loss=0.01) +upper_roof_west = PvWattsSubarrayConfig(name="west", pdc0_w=7*p_module, tilt_deg=10, azimuth_deg=290, dc_loss=0.02, ac_loss=0.01) +cfgs = [upper_roof_north, upper_roof_south, upper_roof_east, upper_roof_west] +pv_plant = PvWattsPlant(site, cfgs) + while True: now = datetime.now() if now.second % interval_seconds == 0 and now.microsecond < 100_000: @@ -42,5 +62,11 @@ while True: else: mode_as_binary = 1 db.store_data('sg_ready', {'mode': mode_as_binary}) + + if now.minute == 0 and now.second == 0: + start_hour_local = (now + dt.timedelta(hours=1)).replace(minute=0, second=0, microsecond=0) + weather = weather_forecaster.get_hourly_forecast(start_hour_local, HORIZON_DAYS) + total = pv_plant.get_power(weather) + db.store_forecasts('pv_forecast', total) time.sleep(0.1) diff --git a/requirements.txt b/requirements.txt index 08e3db0..c52e3d0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,5 @@ pymodbus~=3.8.6 pandas openpyxl -sshtunnel \ No newline at end of file +sshtunnel +pvlib \ No newline at end of file diff --git a/simulators/pv_plant_simulator.py b/simulators/pv_plant_simulator.py new file mode 100644 index 0000000..37c03a9 --- /dev/null +++ b/simulators/pv_plant_simulator.py @@ -0,0 +1,210 @@ +from __future__ import annotations +from dataclasses import dataclass +from typing import Optional, Dict, List, Literal, Tuple, Union + +import numpy as np +import pandas as pd +import pvlib +import matplotlib.pyplot as plt +from pvlib.location import Location +from pvlib.pvsystem import PVSystem +from pvlib.modelchain import ModelChain + +SeriesOrArray = Union[pd.Series, np.ndarray] + +# ----------------------------- Konfiguration ----------------------------- + +@dataclass +class PvWattsSubarrayConfig: + name: str + pdc0_w: float # STC-DC-Leistung [W] + tilt_deg: float # Neigung (0=horizontal) + azimuth_deg: float # Azimut (180=Süd) + gamma_pdc: float = -0.004 # Tempkoeff. [1/K] + eta_inv_nom: float = 0.96 # WR-Wirkungsgrad (nominal) + albedo: float = 0.2 # Bodenreflexion + + # Pauschale Verluste (PVWatts-Losses) + dc_loss: float = 0.0 + ac_loss: float = 0.0 + soiling: float = 0.0 + + # Modell + transposition_model: Literal["perez","haydavies","isotropic","klucher","reindl"] = "perez" + + +# ------------------------------ Subarray --------------------------------- + +class PvWattsSubarray: + """ + Ein Subarray mit pvlib.ModelChain (PVWatts). + Berechnet automatisch DNI/DHI aus GHI (ERBS-Methode) + und nutzt ein SAPM-Temperaturmodell. + """ + def __init__(self, cfg: PvWattsSubarrayConfig, location: Location): + self.cfg = cfg + self.location = location + self._mc: Optional[ModelChain] = None + + # --------------------------------------------------------------------- + def _create_modelchain(self) -> ModelChain: + """Erzeuge eine pvlib.ModelChain-Instanz mit PVWatts-Parametern.""" + temp_params = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS["sapm"]["open_rack_glass_polymer"] + + system = PVSystem( + surface_tilt=self.cfg.tilt_deg, + surface_azimuth=self.cfg.azimuth_deg, + module_parameters={"pdc0": self.cfg.pdc0_w, "gamma_pdc": self.cfg.gamma_pdc}, + inverter_parameters={"pdc0": self.cfg.pdc0_w, "eta_inv_nom": self.cfg.eta_inv_nom}, + albedo=self.cfg.albedo, + temperature_model_parameters=temp_params, + module_type="glass_polymer", + racking_model="open_rack", + ) + + mc = ModelChain( + system, self.location, + transposition_model=self.cfg.transposition_model, + solar_position_method="nrel_numpy", + airmass_model="kastenyoung1989", + dc_model="pvwatts", + ac_model="pvwatts", + aoi_model="physical", + spectral_model=None, + losses_model="pvwatts", + temperature_model="sapm", + ) + + mc.losses_parameters = { + "dc_loss": float(self.cfg.dc_loss), + "ac_loss": float(self.cfg.ac_loss), + "soiling": float(self.cfg.soiling), + } + + self._mc = mc + return mc + + # --------------------------------------------------------------------- + def calc_dni_and_dhi(self, weather: pd.DataFrame) -> pd.DataFrame: + """ + Berechnet DNI & DHI aus GHI über die ERBS-Methode. + Gibt ein neues DataFrame mit 'ghi', 'dni', 'dhi' zurück. + """ + if "ghi" not in weather: + raise ValueError("Wetterdaten benötigen mindestens 'ghi'.") + # Sonnenstand bestimmen + sp = self.location.get_solarposition(weather.index) + erbs = pvlib.irradiance.erbs(weather["ghi"], sp["zenith"], weather.index) + out = weather.copy() + out["dni"] = erbs["dni"].clip(lower=0) + out["dhi"] = erbs["dhi"].clip(lower=0) + return out + + # --------------------------------------------------------------------- + def _prepare_weather(self, weather: pd.DataFrame) -> pd.DataFrame: + """Sichert vollständige Spalten (ghi, dni, dhi, temp_air, wind_speed).""" + if "ghi" not in weather or "temp_air" not in weather: + raise ValueError("weather benötigt Spalten: 'ghi' und 'temp_air'.") + + w = weather.copy() + + # Zeitzone prüfen + if w.index.tz is None: + w.index = w.index.tz_localize(self.location.tz) + else: + if str(w.index.tz) != str(self.location.tz): + w = w.tz_convert(self.location.tz) + + # Wind default + if "wind_speed" not in w: + w["wind_speed"] = 1.0 + + # DNI/DHI ergänzen (immer mit ERBS) + if "dni" not in w or "dhi" not in w: + w = self.calc_dni_and_dhi(w) + + return w + + # --------------------------------------------------------------------- + def get_power(self, weather: pd.DataFrame) -> pd.Series: + """ + Berechnet AC-Leistung aus Wetterdaten. + """ + w = self._prepare_weather(weather) + mc = self._create_modelchain() + mc.run_model(weather=w) + return mc.results.ac.rename(self.cfg.name) + + +# ------------------------------- Anlage ---------------------------------- + +class PvWattsPlant: + """ + Eine PV-Anlage mit mehreren Subarrays, die ein gemeinsames Wetter-DataFrame nutzt. + """ + def __init__(self, site: Location, subarray_cfgs: List[PvWattsSubarrayConfig]): + self.site = site + self.subs: Dict[str, PvWattsSubarray] = {c.name: PvWattsSubarray(c, site) for c in subarray_cfgs} + + def get_power( + self, + weather: pd.DataFrame, + *, + return_breakdown: bool = False + ) -> pd.Series | Tuple[pd.Series, Dict[str, pd.Series]]: + """Berechne Gesamtleistung und optional Einzel-Subarrays.""" + parts: Dict[str, pd.Series] = {name: sub.get_power(weather) for name, sub in self.subs.items()} + + # gemeinsamen Index bilden + idx = list(parts.values())[0].index + for s in parts.values(): + idx = idx.intersection(s.index) + parts = {k: v.reindex(idx).fillna(0.0) for k, v in parts.items()} + + total = sum(parts.values()) + total.name = "total_ac" + + if return_breakdown: + return total, parts + return total + + +# --------------------------- Beispielnutzung ----------------------------- +if __name__ == "__main__": + # Standort + site = Location(latitude=52.52, longitude=13.405, altitude=35, tz="Europe/Berlin", name="Berlin") + + # Zeitachse: 1 Tag, 15-minütig + times = pd.date_range("2025-06-21 00:00", "2025-06-21 23:45", freq="15min", tz=site.tz) + + # Dummy-Wetter + ghi = 1000 * np.clip(np.sin(np.linspace(0, np.pi, len(times)))**1.2, 0, None) + temp_air = 16 + 8 * np.clip(np.sin(np.linspace(-np.pi/2, np.pi/2, len(times))), 0, None) + wind = np.full(len(times), 1.0) + weather = pd.DataFrame(index=times) + weather["ghi"] = ghi + weather["temp_air"] = temp_air + weather["wind_speed"] = wind + + # Zwei Subarrays + cfgs = [ + PvWattsSubarrayConfig(name="Sued_30", pdc0_w=6000, tilt_deg=30, azimuth_deg=180, dc_loss=0.02, ac_loss=0.01), + PvWattsSubarrayConfig(name="West_20", pdc0_w=4000, tilt_deg=20, azimuth_deg=270, soiling=0.02), + ] + plant = PvWattsPlant(site, cfgs) + + # Simulation + total, parts = plant.get_power(weather, return_breakdown=True) + + # Plot + plt.figure(figsize=(10, 6)) + plt.plot(total.index, total / 1000, label="Gesamtleistung (AC)", linewidth=2, color="black") + for name, s in parts.items(): + plt.plot(s.index, s / 1000, label=name) + plt.title("PV-Leistung (PVWatts, ERBS-Methode für DNI/DHI)") + plt.ylabel("Leistung [kW]") + plt.xlabel("Zeit") + plt.legend() + plt.grid(True, linestyle="--", alpha=0.5) + plt.tight_layout() + plt.show()