From 13bad7684ff769d0b714c8e85321fd499bc3a4c5 Mon Sep 17 00:00:00 2001 From: Claude Opus 5 Date: Fri, 31 Jul 2026 10:20:46 +0200 Subject: [PATCH] feat(forecast): project irradiance onto each PV array's plane Radiation-bearing providers previously produced one flat estimate, rated x (W/m2 / 1000), which ignored panel orientation entirely: a south-facing 35 deg roof and a flat one got the same forecast. When per-plane geometry is configured (Weather tab: tilt/azimuth/kWp), global horizontal irradiance is now projected onto each plane via the existing sunpos physics and summed. Providers publishing only GHI get an Erbs correlation to split direct from diffuse first. Sites with no arrays keep the previous behaviour, and Forecast.Solar - which already returns site-calibrated watts - is deliberately left untouched so its numbers are not scaled twice. Split out of the original #718 per review: this half is user-visible and improves every radiation-bearing provider today. The STRANG client and performance scoring that shared that branch move to a follow-up PR. Co-authored-by: HuggeK <48095810+HuggeK@users.noreply.github.com> Signed-off-by: Hugo Karlsson <48095810+HuggeK@users.noreply.github.com> --- .changeset/strang-poa-pv-forecast.md | 15 ++++ go/internal/forecast/forecast.go | 49 ++++++++++- go/internal/forecast/forecast_test.go | 94 ++++++++++++++++++++ go/internal/sunpos/sunpos.go | 118 ++++++++++++++++++++++---- go/internal/sunpos/sunpos_test.go | 87 +++++++++++++++++++ 5 files changed, 341 insertions(+), 22 deletions(-) create mode 100644 .changeset/strang-poa-pv-forecast.md diff --git a/.changeset/strang-poa-pv-forecast.md b/.changeset/strang-poa-pv-forecast.md new file mode 100644 index 00000000..687c8a7e --- /dev/null +++ b/.changeset/strang-poa-pv-forecast.md @@ -0,0 +1,15 @@ +--- +"ftw": minor +--- + +PV forecasts are now orientation-aware. When per-plane geometry is configured +(the Weather tab's PV arrays: tilt/azimuth/kWp), a radiation-bearing forecast +provider's global horizontal irradiance is projected onto each panel plane via +the physics `sunpos` model and summed, instead of the previous flat +`rated × (W/m² / 1000)` estimate that ignored panel orientation. Sites with no +arrays configured keep the existing behaviour, and providers that already return +site-calibrated watts (Forecast.Solar) are left untouched. + +Providers that publish only global horizontal irradiance get an Erbs correlation +to split it into direct and diffuse components before projection, so a +south-facing 35° roof and a flat one no longer receive the same forecast. diff --git a/go/internal/forecast/forecast.go b/go/internal/forecast/forecast.go index 08a55efb..6222871b 100644 --- a/go/internal/forecast/forecast.go +++ b/go/internal/forecast/forecast.go @@ -29,6 +29,7 @@ import ( "github.com/srcfl/ftw/go/internal/config" "github.com/srcfl/ftw/go/internal/state" + "github.com/srcfl/ftw/go/internal/sunpos" ) // Provider is implemented by each weather source. @@ -243,10 +244,16 @@ func EstimatePVW(lat, lon float64, t time.Time, cloudPct *float64, ratedW float6 // Service wraps a provider + store + scheduler for forecasts. type Service struct { - Provider Provider - Store *state.Store - Lat, Lon float64 - RatedPVW float64 // total rated PV across all arrays (used for estimate) + Provider Provider + Store *state.Store + Lat, Lon float64 + RatedPVW float64 // total rated PV across all arrays (used for estimate) + + // Arrays holds per-plane geometry (tilt/azimuth/kWp) mirrored from the + // weather config. When set, a radiation-bearing provider's horizontal + // GHI is projected onto each plane via sunpos and summed, instead of the + // orientation-blind flat rated×(W/m²/1000) estimate. Empty → flat estimate. + Arrays []Array stop chan struct{} done chan struct{} @@ -286,10 +293,21 @@ func FromConfig(cfg *config.Weather, ratedPVW float64, st *state.Store, userAgen default: return nil } + // Mirror per-plane geometry for the POA path. Shared across all + // providers: forecast_solar already applies geometry server-side (so + // this stays unused there), but open_meteo / STRÅNG-style GHI providers + // use it to project irradiance onto each plane in fetchAndStore. + var arrays []Array + for _, a := range cfg.PVArrays { + if a.KWp > 0 { + arrays = append(arrays, Array{TiltDeg: a.TiltDeg, AzimuthDeg: a.AzimuthDeg, KWp: a.KWp}) + } + } return &Service{ Provider: p, Store: st, Lat: cfg.Latitude, Lon: cfg.Longitude, RatedPVW: ratedPVW, + Arrays: arrays, stop: make(chan struct{}), done: make(chan struct{}), } @@ -341,7 +359,11 @@ func (s *Service) fetchAndStore(ctx context.Context) { switch { case r.PVWEstimated != nil: pvW = *r.PVWEstimated + case r.SolarWm2 != nil && len(s.Arrays) > 0: + // Orientation-aware: project GHI onto each configured plane and sum. + pvW = poaPVWattsFromGHI(s.Lat, s.Lon, r.HourStart, *r.SolarWm2, s.Arrays) case r.SolarWm2 != nil && s.RatedPVW > 0: + // No geometry configured — flat, orientation-blind estimate. pvW = s.RatedPVW * (*r.SolarWm2) / 1000.0 default: pvW = EstimatePVW(s.Lat, s.Lon, r.HourStart, r.CloudCoverPct, s.RatedPVW) @@ -365,6 +387,25 @@ func (s *Service) fetchAndStore(ctx context.Context) { slog.Info("forecast fetched", "count", len(points), "provider", s.Provider.Name()) } +// poaPVWattsFromGHI converts a global-horizontal irradiance (W/m², positive) +// into expected DC PV output (W, positive) by projecting it onto each +// configured array's plane via sunpos and scaling by nameplate. This is the +// orientation-aware replacement for the flat rated×(W/m²/1000) estimate; it +// is used whenever the provider supplies GHI and the site has per-plane +// geometry. Returns 0 when the sun is down or no arrays produce output. +func poaPVWattsFromGHI(lat, lon float64, t time.Time, ghiWm2 float64, arrays []Array) float64 { + var total float64 + for _, a := range arrays { + if a.KWp <= 0 { + continue + } + poa := sunpos.POAFromGHI(t, lat, lon, ghiWm2, a.TiltDeg, a.AzimuthDeg) + // kWp×1000 = nameplate W at STC (1000 W/m²); scale by POA/1000. + total += a.KWp * 1000.0 * (poa / 1000.0) + } + return total +} + // Load returns forecasts in [sinceMs, untilMs]. func (s *Service) Load(sinceMs, untilMs int64) ([]state.ForecastPoint, error) { return s.Store.LoadForecasts(sinceMs, untilMs) diff --git a/go/internal/forecast/forecast_test.go b/go/internal/forecast/forecast_test.go index 73450fa0..e0cb5526 100644 --- a/go/internal/forecast/forecast_test.go +++ b/go/internal/forecast/forecast_test.go @@ -247,3 +247,97 @@ func TestFromConfigBuildsMetNo(t *testing.T) { if s.Lat != 59 { t.Errorf("lat: %f", s.Lat) } if s.RatedPVW != 10000 { t.Errorf("rated: %f", s.RatedPVW) } } + +func TestFromConfigPopulatesArrays(t *testing.T) { + st, _ := state.Open(filepath.Join(t.TempDir(), "t.db")) + defer st.Close() + cfg := &config.Weather{ + Provider: "open_meteo", Latitude: 59, Longitude: 18, + PVArrays: []config.PVArray{ + {Name: "south", KWp: 6, TiltDeg: 35, AzimuthDeg: 180}, + {Name: "east", KWp: 4, TiltDeg: 30, AzimuthDeg: 90}, + {Name: "empty", KWp: 0, TiltDeg: 10, AzimuthDeg: 200}, // skipped (kWp 0) + }, + } + s := FromConfig(cfg, 10000, st, "ua") + if s == nil { t.Fatal("expected service") } + if len(s.Arrays) != 2 { + t.Fatalf("expected 2 arrays (kWp>0 only), got %d", len(s.Arrays)) + } + if s.Arrays[0].KWp != 6 || s.Arrays[1].AzimuthDeg != 90 { + t.Errorf("array geometry mismatch: %+v", s.Arrays) + } +} + +// ---- POA-per-array (orientation-aware) estimate ---- + +func TestPOAPVWattsSumsArrays(t *testing.T) { + tt := time.Date(2026, 6, 21, 11, 0, 0, 0, time.UTC) + one := poaPVWattsFromGHI(59.3293, 18.0686, tt, 700, []Array{{TiltDeg: 35, AzimuthDeg: 180, KWp: 5}}) + two := poaPVWattsFromGHI(59.3293, 18.0686, tt, 700, []Array{ + {TiltDeg: 35, AzimuthDeg: 180, KWp: 5}, + {TiltDeg: 35, AzimuthDeg: 180, KWp: 5}, + }) + if one <= 0 { + t.Fatalf("expected positive POA watts, got %.1f", one) + } + if math.Abs(two-2*one) > 1e-6 { + t.Errorf("two identical arrays should double output: one=%.2f two=%.2f", one, two) + } +} + +func TestPOAPVWattsZeroAtNight(t *testing.T) { + tt := time.Date(2026, 12, 21, 23, 0, 0, 0, time.UTC) + w := poaPVWattsFromGHI(59.3293, 18.0686, tt, 500, []Array{{TiltDeg: 35, AzimuthDeg: 180, KWp: 10}}) + if w != 0 { + t.Errorf("night POA watts should be 0, got %.2f", w) + } +} + +// End-to-end: with per-plane geometry, a GHI-bearing provider's stored PV +// estimate comes from the POA path and differs from the flat rated×GHI/1000. +func TestServicePOAPathDiffersFromFlat(t *testing.T) { + handler := http.HandlerFunc(func(w http.ResponseWriter, r *http.Request) { + resp := map[string]any{ + "hourly": map[string]any{ + "time": []string{"2026-06-21T11:00"}, + "shortwave_radiation": []float64{700}, + "cloud_cover": []float64{5}, + "temperature_2m": []float64{20}, + }, + } + _ = json.NewEncoder(w).Encode(resp) + }) + srv := httptest.NewServer(handler) + defer srv.Close() + + st, _ := state.Open(filepath.Join(t.TempDir(), "state.db")) + defer st.Close() + + p := NewOpenMeteo() + p.BaseURL = srv.URL + s := &Service{ + Provider: p, Store: st, Lat: 59.3293, Lon: 18.0686, RatedPVW: 10000, + Arrays: []Array{{TiltDeg: 35, AzimuthDeg: 180, KWp: 10}}, + } + s.fetchAndStore(context.Background()) + + tt := time.Date(2026, 6, 21, 11, 0, 0, 0, time.UTC) + rows, err := st.LoadForecasts(tt.UnixMilli(), tt.Add(time.Hour).UnixMilli()) + if err != nil { + t.Fatal(err) + } + if len(rows) != 1 || rows[0].PVWEstimated == nil { + t.Fatalf("expected 1 forecast with PV estimate, got %+v", rows) + } + got := *rows[0].PVWEstimated + flat := 10000 * 700.0 / 1000.0 // orientation-blind estimate = 7000 W + want := poaPVWattsFromGHI(59.3293, 18.0686, tt, 700, s.Arrays) + if math.Abs(got-want) > 1.0 { + t.Errorf("service should use POA path: got %.1f want %.1f", got, want) + } + if math.Abs(got-flat) < 1.0 { + t.Errorf("POA estimate should differ from flat %.0f, got %.1f", flat, got) + } + t.Logf("POA-per-array estimate %.0fW vs flat %.0fW", got, flat) +} diff --git a/go/internal/sunpos/sunpos.go b/go/internal/sunpos/sunpos.go index 8ce95cc3..787e7c27 100644 --- a/go/internal/sunpos/sunpos.go +++ b/go/internal/sunpos/sunpos.go @@ -125,35 +125,117 @@ func ClearSkyW(t time.Time, lat, lon float64) float64 { } // POA estimates plane-of-array irradiance for one tilted panel using the -// isotropic-sky model. Splits clear-sky horizontal irradiance into beam -// (DNI) and diffuse (DHI) components via a simple Erbs correlation, then -// projects each onto the panel. +// isotropic-sky model, driven by the package's own clear-sky prior. It splits +// that clear-sky horizontal irradiance into beam (DNI) and diffuse (DHI) +// components with a fixed 20% diffuse fraction, then projects each onto the +// panel. Used as the prior signal for the PV twin when no measured irradiance +// is available. // // Returns W/m² on the panel surface; clamped to ≥ 0. +// +// When a data source supplies measured irradiance, prefer the two variants +// below: POAFromComponents (GHI + DHI both known, e.g. SMHI STRÅNG params +// 117 + 122) or POAFromGHI (only GHI known, e.g. Open-Meteo shortwave). func POA(t time.Time, lat, lon, panelTiltDeg, panelAzDeg float64) float64 { sun := At(t, lat, lon) - if sun.ZenithDeg >= 90 { - return 0 - } ghi := ClearSkyW(t, lat, lon) - if ghi <= 0 { + // No measured diffuse component from the clear-sky prior, so keep the + // historical fixed 20% diffuse fraction for this variant. + return POAFromComponents(sun, ghi, 0.2*ghi, panelTiltDeg, panelAzDeg) +} + +// POAFromComponents projects measured global (GHI) and diffuse (DHI) +// horizontal irradiance onto a tilted panel using the isotropic-sky model, +// reusing AOI for the beam projection. All irradiances in W/m²; returns the +// plane-of-array irradiance in W/m², clamped ≥ 0. +// +// Use this when a source gives both GHI and DHI directly (e.g. SMHI STRÅNG +// parameters 117 + 122). When only GHI is available use POAFromGHI, which +// estimates the diffuse split via the Erbs correlation first. +func POAFromComponents(sun Position, ghi, dhi, panelTiltDeg, panelAzDeg float64) float64 { + if sun.ZenithDeg >= 90 || ghi <= 0 { return 0 } - // Erbs et al. (1982) clearness-based diffuse fraction. We don't have - // real GHI measurements, so kt comes from our own clear-sky model → - // always ~0.7-0.75 → diffuse ratio ~0.2. Conservative. - kd := 0.2 - dhi := ghi * kd - dni := (ghi - dhi) / math.Cos(sun.ZenithDeg*math.Pi/180) - + if dhi < 0 { + dhi = 0 + } + if dhi > ghi { + dhi = ghi + } + tiltR := panelTiltDeg * math.Pi / 180 + diffusePOA := dhi * (1 + math.Cos(tiltR)) / 2 // isotropic sky dome aoi := AOI(sun, panelTiltDeg, panelAzDeg) if aoi > 90 { - // Sun behind panel — only diffuse counts. - return dhi * (1 + math.Cos(panelTiltDeg*math.Pi/180)) / 2 + // Sun behind the panel — only diffuse reaches the surface. + return diffusePOA } + cosZ := math.Cos(sun.ZenithDeg * math.Pi / 180) + if cosZ < 0.01 { + // Sun on the horizon: beam projection is ill-conditioned + // (divide-by-~0) and diffuse dominates anyway. + return diffusePOA + } + dni := (ghi - dhi) / cosZ beamPOA := dni * math.Cos(aoi*math.Pi/180) - diffusePOA := dhi * (1 + math.Cos(panelTiltDeg*math.Pi/180)) / 2 out := beamPOA + diffusePOA - if out < 0 { out = 0 } + if out < 0 { + out = 0 + } return out } + +// POAFromGHI projects a measured/forecast global horizontal irradiance (GHI, +// W/m²) onto a tilted panel when no diffuse component is available. It +// estimates the diffuse fraction from the hourly clearness index via the Erbs +// et al. (1982) correlation, then delegates to POAFromComponents. +// +// Use this for radiation providers that expose shortwave/GHI but not diffuse +// (e.g. Open-Meteo shortwave_radiation, or SMHI STRÅNG global-only windows). +func POAFromGHI(t time.Time, lat, lon, ghi, panelTiltDeg, panelAzDeg float64) float64 { + sun := At(t, lat, lon) + if sun.ZenithDeg >= 90 || ghi <= 0 { + return 0 + } + cosZ := math.Cos(sun.ZenithDeg * math.Pi / 180) + i0h := extraterrestrialHorizontalW(t, cosZ) + kt := 0.0 + if i0h > 0 { + kt = ghi / i0h + } + dhi := ghi * ErbsDiffuseFraction(kt) + return POAFromComponents(sun, ghi, dhi, panelTiltDeg, panelAzDeg) +} + +// ErbsDiffuseFraction returns the diffuse fraction (DHI/GHI) for an hourly +// clearness index kt, per Erbs, Klein & Duffie (1982). Result is in +// [0.165, 1]: overcast skies (low kt) are almost entirely diffuse, clear +// skies (high kt) settle near 16.5% diffuse. +func ErbsDiffuseFraction(kt float64) float64 { + switch { + case kt <= 0: + return 1 + case kt <= 0.22: + return 1 - 0.09*kt + case kt <= 0.80: + return 0.9511 - 0.1604*kt + 4.388*kt*kt - 16.638*kt*kt*kt + 12.336*kt*kt*kt*kt + default: + return 0.165 + } +} + +// extraterrestrialHorizontalW returns top-of-atmosphere irradiance on a +// horizontal surface (W/m²) at time t for a solar cosine-zenith cosZ. Used as +// the denominator of the clearness index. Returns 0 when the sun is at/below +// the horizon. +func extraterrestrialHorizontalW(t time.Time, cosZ float64) float64 { + if cosZ <= 0 { + return 0 + } + doy := float64(t.UTC().YearDay()) + gamma := 2 * math.Pi * (doy - 1) / 365 + e0 := 1.000110 + + 0.034221*math.Cos(gamma) + 0.001280*math.Sin(gamma) + + 0.000719*math.Cos(2*gamma) + 0.000077*math.Sin(2*gamma) + const i0 = 1361.0 // solar constant W/m² + return i0 * e0 * cosZ +} diff --git a/go/internal/sunpos/sunpos_test.go b/go/internal/sunpos/sunpos_test.go index efe085fa..fb48f8d5 100644 --- a/go/internal/sunpos/sunpos_test.go +++ b/go/internal/sunpos/sunpos_test.go @@ -66,3 +66,90 @@ func TestPOAFlatEqualsGHI(t *testing.T) { t.Errorf("flat POA should equal GHI: ghi=%.1f poa=%.1f", ghi, poa) } } + +// A flat panel receives all of GHI regardless of the diffuse split, because +// beam-on-horizontal + diffuse-on-horizontal reconstructs GHI exactly. +func TestPOAFromComponentsFlatEqualsGHI(t *testing.T) { + tt := time.Date(2026, 6, 21, 11, 0, 0, 0, time.UTC) + sun := At(tt, 59.33, 18.07) + const ghi = 700.0 + for _, dhi := range []float64{0, 140, 350, 700} { + poa := POAFromComponents(sun, ghi, dhi, 0, 180) + if math.Abs(poa-ghi) > 0.5 { + t.Errorf("flat POA should equal GHI for dhi=%.0f: got %.2f want %.0f", dhi, poa, ghi) + } + } +} + +// A south-tilted panel at Stockholm winter noon (low sun) collects more than a +// flat one — the whole point of projecting onto the plane. +func TestPOAFromComponentsSouthTiltBeatsFlatInWinter(t *testing.T) { + tt := time.Date(2026, 12, 21, 11, 0, 0, 0, time.UTC) // ~noon local, low sun + sun := At(tt, 59.33, 18.07) + if sun.ZenithDeg >= 90 { + t.Skip("sun below horizon") + } + const ghi, dhi = 200.0, 60.0 + flat := POAFromComponents(sun, ghi, dhi, 0, 180) + tilt := POAFromComponents(sun, ghi, dhi, 45, 180) // 45° south + if tilt <= flat { + t.Errorf("south-tilted winter POA (%.1f) should beat flat (%.1f)", tilt, flat) + } +} + +// Sun behind the panel → only the diffuse dome contributes (no negative beam). +func TestPOAFromComponentsSunBehindPanelIsDiffuseOnly(t *testing.T) { + tt := time.Date(2026, 6, 21, 5, 0, 0, 0, time.UTC) // morning, sun in the east + sun := At(tt, 59.33, 18.07) + if sun.ZenithDeg >= 90 { + t.Skip("sun below horizon") + } + const ghi, dhi = 300.0, 90.0 + // A steep west-facing wall can't see the eastern morning sun's beam. + poa := POAFromComponents(sun, ghi, dhi, 90, 270) + diffuseOnly := dhi * (1 + math.Cos(90*math.Pi/180)) / 2 + if math.Abs(poa-diffuseOnly) > 0.5 { + t.Errorf("sun-behind panel should be diffuse-only %.2f, got %.2f", diffuseOnly, poa) + } +} + +// Erbs diffuse fraction: overcast → ~all diffuse, clear → floor at 0.165, +// monotonically non-increasing across the mid range. +func TestErbsDiffuseFraction(t *testing.T) { + if f := ErbsDiffuseFraction(0); f != 1 { + t.Errorf("kt=0 → fraction 1, got %.3f", f) + } + if f := ErbsDiffuseFraction(1.0); math.Abs(f-0.165) > 1e-9 { + t.Errorf("kt>0.8 → 0.165, got %.3f", f) + } + clear := ErbsDiffuseFraction(0.75) + murky := ErbsDiffuseFraction(0.30) + if !(clear < murky) { + t.Errorf("clearer sky should have less diffuse: clear=%.3f murky=%.3f", clear, murky) + } + if clear < 0.165 || clear > 1 { + t.Errorf("fraction out of [0.165,1]: %.3f", clear) + } +} + +// POAFromGHI (Erbs split) on a flat panel still reconstructs ~GHI. +func TestPOAFromGHIFlatApproxGHI(t *testing.T) { + tt := time.Date(2026, 6, 21, 11, 0, 0, 0, time.UTC) + const ghi = 600.0 + poa := POAFromGHI(tt, 59.33, 18.07, ghi, 0, 180) + if math.Abs(poa-ghi) > 0.5 { + t.Errorf("flat POAFromGHI should equal GHI: got %.2f want %.0f", poa, ghi) + } +} + +// Night → zero from both measured-irradiance variants regardless of input. +func TestPOAVariantsZeroAtNight(t *testing.T) { + tt := time.Date(2026, 12, 21, 23, 0, 0, 0, time.UTC) + sun := At(tt, 59.33, 18.07) + if p := POAFromComponents(sun, 500, 100, 35, 180); p != 0 { + t.Errorf("night POAFromComponents should be 0, got %.2f", p) + } + if p := POAFromGHI(tt, 59.33, 18.07, 500, 35, 180); p != 0 { + t.Errorf("night POAFromGHI should be 0, got %.2f", p) + } +}