Fit a clear-sky 'ideal produce' PV curve from Home Assistant recorder history and emit a template sensor (companion to the harizanov.com clear-sky line post)

fit_clearsky.py

#!/usr/bin/env python3
"""Fit a clear-sky ("ideal produce") PV curve to your own Home Assistant history.

The idea: clouds only ever subtract. Bin a few years of hourly PV output by sun
position and the top ~2-4% of every bin is what your array does under a clear sky.
Fit a small physical model to that envelope and you get an "ideal produce" curve
you can overlay on actual production - no panel datasheet or roof measuring needed.
At the end this script prints a ready-to-paste Home Assistant template sensor.

Companion post: https://harizanov.com/  (search "clear-sky line")

Usage
-----
1. Find your PV power sensor in the recorder DB (default HAOS path shown):

     sqlite3 -readonly /homeassistant/home-assistant_v2.db \
       "SELECT id, statistic_id FROM statistics_meta WHERE statistic_id LIKE '%photovolta%';"

2. Export hourly long-term statistics for it (and optionally an OUTDOOR
   temperature sensor - genuinely outdoor, check it reads winter-cold in January):

     sqlite3 -readonly -csv /homeassistant/home-assistant_v2.db \
       "SELECT start_ts, mean FROM statistics WHERE metadata_id=<PV_ID> ORDER BY start_ts;" > pv_hourly.csv
     sqlite3 -readonly -csv /homeassistant/home-assistant_v2.db \
       "SELECT start_ts, mean FROM statistics WHERE metadata_id=<TEMP_ID> ORDER BY start_ts;" > temp_hourly.csv

3. Fill in SITE CONFIG below, `pip install numpy scipy`, run the script.

Notes
-----
- The more full years of data the better; 2+ years covers all seasons.
- Export/feed-in limits poison the envelope: if your inverter curtails at a flat
  ceiling (mine did at 5.0 kW for two summers), those hours read as "clear sky"
  lower than reality. Pool years where the limit was off, or expect the fit to
  underestimate midday summer. Check with a histogram of near-noon values.
- Timestamps are unix UTC throughout; sun position is computed in UTC, so DST
  never enters the picture.
- Written for a single-orientation array in the northern hemisphere; widen the
  panel-azimuth bounds in FIT if yours faces far from south.
"""
import numpy as np
from scipy.optimize import least_squares

# ---------------- SITE CONFIG ----------------
LAT = None          # decimal degrees, e.g. 42.75  (keep your exact coords private!)
LON = None          # decimal degrees east
ALT_M = None        # site altitude in meters, e.g. 720
PV_CSV = 'pv_hourly.csv'
TEMP_CSV = 'temp_hourly.csv'          # set to None if you have no outdoor temp history
TEMP_ENTITY = 'sensor.your_outdoor_temperature'  # used in the emitted HA template
# ----------------------------------------------

assert None not in (LAT, LON, ALT_M), "fill in LAT / LON / ALT_M first"
PRESS = np.exp(-ALT_M / 8435.0)  # broadband pressure correction for altitude


def load_csv(path, ncols):
    rows = []
    with open(path) as f:
        for line in f:
            parts = line.strip().split(',')
            if len(parts) < ncols:
                continue
            try:
                rows.append([float(p) if p != '' else np.nan for p in parts[:ncols]])
            except ValueError:
                pass
    return np.array(rows)


def solar_position(ts):
    """NOAA solar position; ts = unix seconds (array). Returns apparent elevation
    and azimuth (degrees, azimuth 0=N clockwise), incl. refraction correction."""
    ts = np.asarray(ts, dtype=float)
    jd = ts / 86400.0 + 2440587.5
    T = (jd - 2451545.0) / 36525.0
    L0 = np.mod(280.46646 + T * (36000.76983 + T * 0.0003032), 360)
    M = 357.52911 + T * (35999.05029 - 0.0001537 * T)
    Mr = np.radians(M)
    ecc = 0.016708634 - T * (0.000042037 + 0.0000001267 * T)
    C = (np.sin(Mr) * (1.914602 - T * (0.004817 + 0.000014 * T))
         + np.sin(2 * Mr) * (0.019993 - 0.000101 * T) + np.sin(3 * Mr) * 0.000289)
    omega = 125.04 - 1934.136 * T
    lam = L0 + C - 0.00569 - 0.00478 * np.sin(np.radians(omega))
    eps = (23 + (26 + (21.448 - T * (46.8150 + T * (0.00059 - T * 0.001813))) / 60) / 60
           + 0.00256 * np.cos(np.radians(omega)))
    epsr, lamr = np.radians(eps), np.radians(lam)
    decl = np.arcsin(np.sin(epsr) * np.sin(lamr))
    y = np.tan(epsr / 2) ** 2
    L0r = np.radians(L0)
    eqtime = 4 * np.degrees(y * np.sin(2 * L0r) - 2 * ecc * np.sin(Mr)
                            + 4 * ecc * y * np.sin(Mr) * np.cos(2 * L0r)
                            - 0.5 * y * y * np.sin(4 * L0r) - 1.25 * ecc * ecc * np.sin(2 * Mr))
    tst = np.mod(ts, 86400.0) / 60.0 + eqtime + 4 * LON
    ha = np.mod(tst / 4.0 - 180.0 + 180, 360) - 180
    har, latr = np.radians(ha), np.radians(LAT)
    cos_zen = np.clip(np.sin(latr) * np.sin(decl) + np.cos(latr) * np.cos(decl) * np.cos(har), -1, 1)
    el = 90 - np.degrees(np.arccos(cos_zen))
    elr = np.radians(np.clip(el, -1, 89))
    refr = np.where(el > 85, 0,
                    np.where(el > 5, 58.1 / np.tan(elr) - 0.07 / np.tan(elr) ** 3 + 0.000086 / np.tan(elr) ** 5,
                             np.where(el > -0.575,
                                      1735 + el * (-518.2 + el * (103.4 + el * (-12.79 + el * 0.711))),
                                      -20.772 / np.tan(elr)))) / 3600.0
    sin_zen = np.sin(np.arccos(cos_zen)) + 1e-12
    sin_az = -np.cos(decl) * np.sin(har) / sin_zen
    cos_az = (np.sin(decl) - np.sin(latr) * cos_zen) / (np.cos(latr) * sin_zen)
    return el + refr, np.mod(np.degrees(np.arctan2(sin_az, cos_az)), 360)


def doy_frac(ts):
    dt64 = np.asarray(ts, dtype=float).astype('datetime64[s]')
    return (dt64 - dt64.astype('datetime64[Y]')).astype('timedelta64[s]').astype(float) / 86400.0 + 1


# ---------------- load & scrub ----------------
pv = load_csv(PV_CSV, 2)
ts, p = pv[:, 0], pv[:, 1]
ok = ~np.isnan(p)
ts, p = ts[ok], p[ok]
# stuck-sensor flatlines (same value >=3h while producing) and nonzero night rows
stuck = np.zeros(len(p), bool)
for i in range(2, len(p)):
    if p[i] > 100 and p[i] == p[i - 1] == p[i - 2]:
        stuck[i - 2:i + 1] = True
el0, _ = solar_position(ts + 1800)
drop = stuck | ((el0 < -3) & (p > 50))
print(f"scrubbed {drop.sum()} suspect rows")
ts, p = ts[~drop], p[~drop]

# ambient temperature: hourly history if provided, else seasonal sinusoid fit
if TEMP_CSV:
    tmp = load_csv(TEMP_CSV, 2)
    tok = ~np.isnan(tmp[:, 1])
    tts, tvs = tmp[tok, 0], tmp[tok, 1]
    tdoy, thr = doy_frac(tts), np.mod(tts, 86400) / 3600.0
    A = np.column_stack([np.ones_like(tdoy),
                         np.cos(2 * np.pi * tdoy / 365.25), np.sin(2 * np.pi * tdoy / 365.25),
                         np.cos(2 * np.pi * thr / 24), np.sin(2 * np.pi * thr / 24)])
    tc, *_ = np.linalg.lstsq(A, tvs, rcond=None)
    temp_map = dict(zip(tmp[:, 0], tmp[:, 1]))

    def get_temp(tsx):
        d, h = doy_frac(tsx), np.mod(tsx, 86400) / 3600.0
        out = (tc[0] + tc[1] * np.cos(2 * np.pi * d / 365.25) + tc[2] * np.sin(2 * np.pi * d / 365.25)
               + tc[3] * np.cos(2 * np.pi * h / 24) + tc[4] * np.sin(2 * np.pi * h / 24))
        for i, t in enumerate(tsx):
            v = temp_map.get(t - np.mod(t, 3600), np.nan)
            if not np.isnan(v):
                out[i] = v
        return out
else:
    def get_temp(tsx):
        return np.full(len(tsx), 15.0)  # crude; thermal derate folds into kwp

# ---------------- model ----------------
# params: tilt, panel_az, kwp_eff, tau0, tau_cos, tau_sin, am_exp, kdiff, horiz_e, horiz_w, clip_w
def model_power(params, el, az, doy, tamb):
    tilt, azp, kwp, tau0, tauc, taus, bexp, kd, hde, hdw, pclip = params
    elr = np.radians(np.maximum(el, 0.0))
    sinel = np.sin(elr)
    up = el > 0.5
    E0 = 1361.0 * (1 + 0.033 * np.cos(2 * np.pi * (doy - 3) / 365.25))
    am = PRESS / (sinel + 0.50572 * (np.maximum(el, 0) + 6.07995) ** -1.6364)
    ang = 2 * np.pi * doy / 365.25
    tau = tau0 * (1 + tauc * np.cos(ang) + taus * np.sin(ang))  # seasonal turbidity
    dni = E0 * np.exp(-tau * am ** bexp) * up
    hor = np.where(az < 180, hde, hdw)          # east/west horizon obstruction
    dni = dni / (1 + np.exp(-(el - hor)))
    dhi = kd * E0 * np.maximum(sinel, 0) ** 0.6 * up
    tr = np.radians(tilt)
    cosaoi = sinel * np.cos(tr) + np.cos(elr) * np.sin(tr) * np.cos(np.radians(az - azp))
    poa = (dni * np.maximum(cosaoi, 0) + dhi * (1 + np.cos(tr)) / 2
           + 0.2 * (dni * sinel + dhi) * (1 - np.cos(tr)) / 2)
    pdc = kwp * poa * (1 - 0.0038 * (tamb + 0.032 * poa - 25.0))
    return np.minimum(np.maximum(pdc, 0), pclip)


# hourly stats are means, so average the model over 6 sub-samples per hour
SUB = np.arange(300, 3600, 600.0)
tsub = (ts[:, None] + SUB[None, :]).ravel()
el_sub, az_sub = solar_position(tsub)
el_sub, az_sub = el_sub.reshape(len(ts), 6), az_sub.reshape(len(ts), 6)
el_mid, az_mid = el_sub[:, 2:4].mean(1), az_sub[:, 2:4].mean(1)
tamb = get_temp(ts)
doy = doy_frac(ts)


def model_hourly(params):
    out = np.zeros(len(ts))
    for j in range(6):
        out += model_power(params, el_sub[:, j], az_sub[:, j], doy, tamb)
    return out / 6.0


# ---------------- envelope selection ----------------
day = el_mid > 1
binid = np.clip((el_mid / 2).astype(int), 0, 40) * 100 + np.clip(((az_mid - 40) / 5).astype(int), 0, 60)
env_idx = []
for b in np.unique(binid[day]):
    idx = np.where(day & (binid == b))[0]
    if len(idx) < 8:
        continue
    lo, hi = np.percentile(p[idx], [95.5, 99.8])  # top slice, extreme outliers off
    env_idx.extend(idx[(p[idx] >= lo) & (p[idx] <= hi)])
env_idx = np.array(sorted(env_idx))
print(f"envelope points: {len(env_idx)} of {day.sum()} day hours")

# ---------------- FIT ----------------
pmax = np.nanmax(p)
x0 = np.array([25.0, 180.0, pmax / 1100, 0.28, -0.1, 0.0, 0.65, 0.10, 2.0, 2.0, pmax])
lb = np.array([5.0, 120.0, 0.3 * x0[2], 0.08, -0.6, -0.6, 0.3, 0.02, -2.0, -2.0, 0.85 * pmax])
ub = np.array([60.0, 240.0, 2.0 * x0[2], 0.80, 0.6, 0.6, 2.5, 0.30, 12.0, 12.0, 1.25 * pmax])


def resid(params):
    r = model_hourly(params)[env_idx] - p[env_idx]
    return r * np.where(r < 0, 2.0, 1.0) / 100.0  # asymmetric: model should cap the envelope


fit = least_squares(resid, x0, bounds=(lb, ub))
pars = fit.x
names = ['tilt', 'panel_az', 'kwp_eff', 'tau0', 'tau_cos', 'tau_sin', 'am_exp', 'kdiff', 'horiz_e', 'horiz_w', 'clip_w']
print("\nfitted parameters:")
for n, v in zip(names, pars):
    print(f"  {n:9s} = {v:.4f}")

# ---------------- validation ----------------
m = model_hourly(pars)
r = p[env_idx] - m[env_idx]
print(f"\nenvelope fit: rmse={np.sqrt(np.mean(r ** 2)):.0f} W  bias={np.mean(r):.0f} W")
exc = day & (p > m)
print(f"hours actual>model: {100 * exc.sum() / day.sum():.1f}%  (expect ~3-4% for an envelope)")
ldays = ((ts + 7200) // 86400).astype(int)
ratios = []
for d in np.unique(ldays):
    msk = (ldays == d) & day
    if msk.sum() >= 4 and m[msk].sum() > 1000:
        ratios.append(p[msk].sum() / m[msk].sum())
ratios = np.array(ratios)
print(f"daily actual/ideal ratio: p50={np.percentile(ratios, 50):.2f}  p98={np.percentile(ratios, 98):.2f} "
      f"(p98 should sit just under 1.0; clear days ~0.95+)")

# ---------------- emit HA template sensor ----------------
tilt, azp, kwp, tau0, tauc, taus, bexp, kd, hde, hdw, pclip = pars
print(f"""
# ---- paste into a package / configuration.yaml ----
template:
  - sensor:
      - name: "Solar ideal produce"
        unique_id: solar_ideal_produce_clearsky
        unit_of_measurement: "W"
        device_class: power
        state_class: measurement
        state: >-
          {{% set el = state_attr('sun.sun','elevation') | float(-90) %}}
          {{% if el <= 0.5 %}}
            0
          {{% else %}}
            {{% set az = state_attr('sun.sun','azimuth') | float(180) %}}
            {{% set doy = now().timetuple().tm_yday %}}
            {{% set tamb = states('{TEMP_ENTITY}') | float(15) %}}
            {{% set w = 2 * pi * doy / 365.25 %}}
            {{% set rad = pi / 180 %}}
            {{% set elr = el * rad %}}
            {{% set sinel = sin(elr) %}}
            {{% set e0 = 1361 * (1 + 0.033 * cos(2 * pi * (doy - 3) / 365.25)) %}}
            {{% set am = {PRESS:.4f} / (sinel + 0.50572 * (el + 6.07995) ** (-1.6364)) %}}
            {{% set tau = {tau0:.4f} * (1 + {tauc:.4f} * cos(w) + {taus:.4f} * sin(w)) %}}
            {{% set dni = e0 * e ** (-tau * am ** {bexp:.4f}) %}}
            {{% set hor = {hde:.4f} if az < 180 else {hdw:.4f} %}}
            {{% set dni = dni / (1 + e ** (-(el - hor))) %}}
            {{% set dhi = {kd:.4f} * e0 * sinel ** 0.6 %}}
            {{% set tr = {tilt:.4f} * rad %}}
            {{% set cosaoi = [sinel * cos(tr) + cos(elr) * sin(tr) * cos((az - {azp:.4f}) * rad), 0] | max %}}
            {{% set poa = dni * cosaoi + dhi * (1 + cos(tr)) / 2 + 0.2 * (dni * sinel + dhi) * (1 - cos(tr)) / 2 %}}
            {{% set pdc = {kwp:.4f} * poa * (1 - 0.0038 * (tamb + 0.032 * poa - 25)) %}}
            {{{{ [[pdc, 0] | max, {pclip:.1f}] | min | round(0) }}}}
          {{% endif %}}
""")