751 lines
34 KiB
Python
751 lines
34 KiB
Python
#!/usr/bin/env python
|
|
# -*- coding: utf-8 -*-
|
|
"""
|
|
Parser inteligente de logs del WiT-CPV
|
|
(webpage_dynamictrack_autotune_log_psram.ino)
|
|
===============================================================================
|
|
|
|
Formato de lineas soportado:
|
|
|
|
[AAAA-MM-DD HH:MM:SS] X:<x>,Y:<y>,C:<c>,V:<v>,P:<p>,S:<s>,er:<er> -> muestra periodica (~2 s)
|
|
[ts] FineTune: dC=.., dynE=.., e=.., near=0/1, Cmax=.. -> ciclo de fine-tune
|
|
[ts] SunPos: Az=.., Alt=.., Az_corr=.., Alt_corr=.., posX=.., posY=.. -> prediccion astronomica
|
|
[ts] Ejecutando comando: N -> cambio de estado
|
|
[ts] ErrMax=.. / MinErrMax=.. / CicleMax=.. / Az_deviation=.. -> parametros por web
|
|
[ts] AT: rec small=.., abrupt=.., reset=.., peak_prox=.. -> auto-tune
|
|
[ts] demas eventos (Timeinfo, Debug, Recibi, Moviendo, ...)
|
|
|
|
Calibracion (identica a metrica1.py), ADS122C04 de 23 bit, Vref = 2.048 V:
|
|
|
|
G_POA [W/m2] = P_raw/8388607*2048/2 (piranometro, divisor /2)
|
|
Voc [V] = V_raw/8388607*2.048*2 (divisor 10k/10k)
|
|
Isc [A] = (C_raw/8388607*2048/20)/1e3/0.075 (AD8418 20x, shunt 75 mOhm)
|
|
P_CPV [W] = Voc * Isc * FF con FF = 0.82
|
|
|
|
Geometria del soporte (trackSunByTime en el firmware):
|
|
|
|
x = (90 - alt_corr)/80*53000 + 53000 + Az_deviation/90*y
|
|
y = (az_corr - 180)/80*53000 + 53000
|
|
=> cada eje recorre 80 grados con 53000 micro-pasos (1 ustep = 80/53000 deg)
|
|
|
|
Salida: articulo/resultados_<tag>/{samples.csv, events.csv, metrics.json}
|
|
articulo/figs/fig_{power,position,finetune,recovery}_<tag>.pdf
|
|
|
|
Uso: py -3.12 parse_log.py 21092025.txt
|
|
"""
|
|
|
|
import json
|
|
import re
|
|
import sys
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
import pandas as pd
|
|
|
|
try:
|
|
import matplotlib
|
|
matplotlib.use("Agg")
|
|
import matplotlib.dates as mdates
|
|
import matplotlib.pyplot as plt
|
|
HAS_MPL = True
|
|
except Exception:
|
|
HAS_MPL = False
|
|
|
|
# ----------------------------------------------------------------------------
|
|
# Calibracion (identica a metrica1.py)
|
|
# ----------------------------------------------------------------------------
|
|
ADC_MAX = 8388607.0 # 2^23 - 1
|
|
VREF_MV = 2048.0 # Vref interno del ADS122C04 (mV)
|
|
GAIN_SHUNT = 20.0 # ganancia del amplificador de instrumentacion
|
|
R_SHUNT = 0.075 # Ohm
|
|
FF = 0.82 # fill factor caracterizado
|
|
|
|
# Geometria del soporte
|
|
STEPS_PER_80DEG = 53000.0 # micro-pasos por eje para 80 grados
|
|
DEG_PER_USTEP = 80.0 / STEPS_PER_80DEG
|
|
|
|
ROOT = Path(__file__).resolve().parent
|
|
ARTICULO = ROOT / "articulo"
|
|
|
|
|
|
def _trapz(y, x):
|
|
"""Integracion trapezoidal compatible con numpy >= 2 (np.trapezoid)."""
|
|
if hasattr(np, "trapezoid"):
|
|
return float(np.trapezoid(y, x=x))
|
|
return float(np.trapz(y, x=x))
|
|
|
|
|
|
# ----------------------------------------------------------------------------
|
|
# Expresiones regulares de los eventos del log
|
|
# ----------------------------------------------------------------------------
|
|
RE_SAMPLE = re.compile(
|
|
r"^X:(?P<x>-?\d+),Y:(?P<y>-?\d+),C:(?P<c>\d+),V:(?P<v>\d+),"
|
|
r"P:(?P<p>\d+),S:(?P<s>\d+),er:(?P<er>\d+)\s*$"
|
|
)
|
|
RE_FINE = re.compile(
|
|
r"FineTune: dC=(?P<dc>\d+), dynE=(?P<dynE>\d+), e=(?P<e>\d+), "
|
|
r"near=(?P<near>\d), Cmax=(?P<cmax>\d+)"
|
|
)
|
|
RE_SUNPOS_NEW = re.compile(
|
|
r"SunPos: Az=(?P<az>[\d.]+), Alt=(?P<alt>[\d.]+), "
|
|
r"Az_corr=(?P<azc>[\d.]+), Alt_corr=(?P<altc>[\d.]+), "
|
|
r"posX=(?P<px>\d+), posY=(?P<py>\d+)"
|
|
)
|
|
RE_SUNPOS_OLD = re.compile(
|
|
r"SunPos: Az=(?P<az>[\d.]+), Zen=(?P<zen>[\d.]+) \(raw=[\d.]+\), "
|
|
r"posX=(?P<px>\d+), posY=(?P<py>\d+)"
|
|
)
|
|
RE_CMD = re.compile(r"Ejecutando comando: (\d+)")
|
|
RE_PARAM = re.compile(r"^(ErrMax|MinErrMax|CicleMax|Az_deviation|TiltAngle)=(\S+)")
|
|
RE_AT_REC = re.compile(
|
|
r"AT: rec small=(?P<small>\d+), abrupt=(?P<abrupt>\d+), "
|
|
r"reset=(?P<reset>\d+), peak_prox=(?P<peak>[\d.]+)"
|
|
)
|
|
RE_AT_MEAN = re.compile(r"AT: meanI=([\d.]+)")
|
|
|
|
|
|
def parse_log(file_path):
|
|
"""Lee el log y devuelve (samples, events) como DataFrames."""
|
|
samples = []
|
|
events = []
|
|
with open(file_path, "r", encoding="utf-8", errors="replace") as f:
|
|
for line in f:
|
|
line = line.rstrip("\n")
|
|
# El firmware antiguo escribe cada linea con prefijo "LOG:" (antes
|
|
# del timestamp). Retirarlo aqui (es un no-op en el formato nuevo)
|
|
# permite que el mismo parser lea ambos formatos de log.
|
|
if line.startswith("LOG:"):
|
|
line = line[4:]
|
|
m = re.match(r"^\[(\d{4}-\d{2}-\d{2} \d{2}:\d{2}:\d{2})\] (.*)$", line)
|
|
if not m:
|
|
continue
|
|
ts, body = m.group(1), m.group(2)
|
|
|
|
ms = RE_SAMPLE.match(body)
|
|
if ms:
|
|
samples.append(
|
|
{
|
|
"timestamp": ts,
|
|
"x": int(ms["x"]),
|
|
"y": int(ms["y"]),
|
|
"c_raw": int(ms["c"]),
|
|
"v_raw": int(ms["v"]),
|
|
"p_raw": int(ms["p"]),
|
|
"s_uptime": int(ms["s"]),
|
|
"er": int(ms["er"]),
|
|
}
|
|
)
|
|
continue
|
|
|
|
mf = RE_FINE.search(body)
|
|
if mf:
|
|
events.append(
|
|
{"timestamp": ts, "kind": "fine",
|
|
"dC": int(mf["dc"]), "dynE": int(mf["dynE"]),
|
|
"e": int(mf["e"]), "near": int(mf["near"]),
|
|
"cmax": int(mf["cmax"])}
|
|
)
|
|
continue
|
|
|
|
mp = RE_SUNPOS_NEW.search(body)
|
|
if mp is None:
|
|
mp = RE_SUNPOS_OLD.search(body)
|
|
if mp:
|
|
d = mp.groupdict()
|
|
if "alt" in d and d["alt"] is not None:
|
|
alt, altc = float(d["alt"]), float(d["altc"])
|
|
else:
|
|
zen = float(d["zen"])
|
|
alt = 90.0 - zen
|
|
altc = 90.0 - zen
|
|
events.append(
|
|
{"timestamp": ts, "kind": "sunpos", "az": float(d["az"]),
|
|
"alt": alt, "azc": float(d.get("azc", d["az"])),
|
|
"altc": altc, "px": int(d["px"]), "py": int(d["py"])}
|
|
)
|
|
continue
|
|
|
|
mc = RE_CMD.search(body)
|
|
if mc:
|
|
events.append({"timestamp": ts, "kind": "cmd", "cmd": int(mc.group(1))})
|
|
continue
|
|
|
|
mp = RE_PARAM.match(body)
|
|
if mp:
|
|
events.append({"timestamp": ts, "kind": "param",
|
|
"name": mp.group(1), "value": float(mp.group(2))})
|
|
continue
|
|
|
|
mat = RE_AT_REC.search(body)
|
|
if mat:
|
|
events.append({"timestamp": ts, "kind": "at_rec",
|
|
**{k: float(mat[k])
|
|
for k in ("small", "abrupt", "reset", "peak")}})
|
|
continue
|
|
|
|
if body.startswith("AT:"):
|
|
events.append({"timestamp": ts, "kind": "at_info", "text": body})
|
|
continue
|
|
|
|
# eventos menores (contexto)
|
|
matched = False
|
|
for prefix in ("Timeinfo:", "UT_hour:", "Debug:", "Recibi", "Moviendo"):
|
|
if body.startswith(prefix):
|
|
events.append({"timestamp": ts, "kind": "info", "text": body})
|
|
matched = True
|
|
break
|
|
if not matched:
|
|
events.append({"timestamp": ts, "kind": "other", "text": body})
|
|
|
|
df_s = pd.DataFrame(samples)
|
|
if df_s.empty:
|
|
raise ValueError("Sin muestras X,Y,C,V,P en el log.")
|
|
df_s["timestamp"] = pd.to_datetime(df_s["timestamp"])
|
|
df_e = pd.DataFrame(events)
|
|
if not df_e.empty:
|
|
df_e["timestamp"] = pd.to_datetime(df_e["timestamp"])
|
|
return df_s, df_e
|
|
|
|
|
|
|
|
# ----------------------------------------------------------------------------
|
|
# Unidades fisicas + maquina de estados + referencia P_ref
|
|
# ----------------------------------------------------------------------------
|
|
CMD_NAMES = {2: "Track solar", 3: "Fine-tune", 4: "Homing",
|
|
5: "Recentering", 6: "Set position", 10: "Driver status", 11: "Demo"}
|
|
|
|
|
|
def build_physical(df_s, df_e):
|
|
"""Añade unidades físicas, P_ref, estado/fase y posicion objetivo."""
|
|
# 1) unidades fisicas (calibracion de metrica1.py)
|
|
df_s["G_POA"] = df_s["p_raw"] / ADC_MAX * (VREF_MV / 2.0) # W/m2
|
|
df_s["Voc"] = df_s["v_raw"] / ADC_MAX * 2.048 * 2.0 # V
|
|
df_s["Isc"] = (df_s["c_raw"] / ADC_MAX * (VREF_MV / GAIN_SHUNT)
|
|
/ 1000.0) / R_SHUNT # A
|
|
df_s["P_CPV"] = df_s["Voc"] * df_s["Isc"] * FF # W
|
|
|
|
# 2) maquina de estados a partir de los comandos recibidos
|
|
status = 0
|
|
phase = "Startup"
|
|
phase_marks = []
|
|
prev_ts = df_s["timestamp"].iloc[0]
|
|
cmds = (df_e[df_e["kind"] == "cmd"] if not df_e.empty
|
|
else pd.DataFrame(columns=["timestamp", "cmd"]))
|
|
for ts_c, cmd in zip(cmds["timestamp"], cmds["cmd"]):
|
|
phase_marks.append((prev_ts, ts_c, phase, status))
|
|
if cmd in (2, 3):
|
|
status = 3
|
|
phase = "Solar tracking"
|
|
elif cmd == 4:
|
|
phase = "Homing"
|
|
elif cmd == 5:
|
|
phase = "Recentering"
|
|
elif cmd == 11:
|
|
phase = "Demo"
|
|
else:
|
|
phase = CMD_NAMES.get(cmd, f"Cmd {cmd}")
|
|
prev_ts = ts_c
|
|
phase_marks.append((prev_ts, df_s["timestamp"].iloc[-1], phase, status))
|
|
|
|
# 3) asignar estado/fase a cada muestra
|
|
ph = pd.DataFrame(phase_marks, columns=["t0", "t1", "phase", "status"])
|
|
idx = ph["t0"].values.searchsorted(df_s["timestamp"].values, side="right") - 1
|
|
idx = np.clip(idx, 0, len(ph) - 1)
|
|
df_s["status"] = ph["status"].values[idx]
|
|
df_s["phase"] = ph["phase"].values[idx]
|
|
|
|
# 4) objetivo astronomico mas reciente (fwd-fill) + Ang_dev aplicado
|
|
if not df_e.empty:
|
|
sp = df_e[df_e["kind"] == "sunpos"].copy()
|
|
if not sp.empty:
|
|
tgt = pd.DataFrame({"timestamp": sp["timestamp"],
|
|
"tgt_x": sp["px"], "tgt_y": sp["py"],
|
|
"az": sp["az"], "alt": sp["alt"],
|
|
"azc": sp["azc"], "altc": sp["altc"]})
|
|
merged = pd.merge_asof(df_s[["timestamp"]], tgt,
|
|
on="timestamp", direction="backward")
|
|
for col in ("tgt_x", "tgt_y", "az", "alt", "azc", "altc"):
|
|
df_s[col] = merged[col]
|
|
azdev = 0.0
|
|
par = df_e[df_e["kind"] == "param"]
|
|
if not par.empty:
|
|
azp = par[par["name"] == "Az_deviation"]
|
|
if not azp.empty:
|
|
azser = pd.Series(azp["value"].values, index=azp["timestamp"])
|
|
ff = azser.reindex(df_s["timestamp"], method="ffill").dropna()
|
|
azdev = float(ff.iloc[0]) if len(ff) else 0.0
|
|
# angulos reconstruidos invirtiendo el mapeo del firmware
|
|
x_base = df_s["x"] - azdev / 90.0 * df_s["y"]
|
|
df_s["alt_est"] = 90.0 - (x_base - STEPS_PER_80DEG) / STEPS_PER_80DEG * 80.0
|
|
df_s["az_est"] = 180.0 + (df_s["y"] - STEPS_PER_80DEG) / STEPS_PER_80DEG * 80.0
|
|
df_s["azdev_used"] = azdev
|
|
else:
|
|
for col in ("tgt_x", "tgt_y", "az", "alt", "azc", "altc", "azdev_used"):
|
|
df_s[col] = np.nan
|
|
df_s["azdev_used"] = 0.0
|
|
|
|
# 5) ciclos de fine-tune alineados por timestamp (merge asof)
|
|
if not df_e.empty:
|
|
ft = df_e[df_e["kind"] == "fine"].copy()
|
|
if not ft.empty:
|
|
df_s = pd.merge_asof(df_s, ft[["timestamp", "e", "dC", "dynE",
|
|
"near", "cmax"]].sort_values("timestamp"),
|
|
on="timestamp", direction="backward")
|
|
else:
|
|
for col in ("e", "dC", "dynE", "near", "cmax"):
|
|
df_s[col] = np.nan
|
|
|
|
# 6) referencia optima: P_ref(t) = R_opt,max * G(t) * Voc(t) * FF
|
|
# R_opt,max se estima bajo cielo directo (G >= 600 W/m2) para que los
|
|
# transitorios nubosos (respuesta angular del piranometro vs optica
|
|
# concentradora) no inflen la responsividad de referencia.
|
|
lit = df_s[df_s["G_POA"] > 600.0]
|
|
R_opt = float((lit["Isc"] / lit["G_POA"]).max()) if len(lit) else np.nan
|
|
df_s["R_opt"] = R_opt
|
|
df_s["Isc_ref"] = R_opt * df_s["G_POA"]
|
|
df_s["P_ref"] = R_opt * df_s["G_POA"] * df_s["Voc"] * FF
|
|
|
|
# 7) actividad mecanica e intervalos temporales
|
|
df_s["dx"] = df_s["x"].diff().abs()
|
|
df_s["dy"] = df_s["y"].diff().abs()
|
|
df_s["dt_s"] = df_s["timestamp"].diff().dt.total_seconds().fillna(0.0)
|
|
|
|
return df_s, ph
|
|
|
|
|
|
|
|
# ----------------------------------------------------------------------------
|
|
# Metricas (seccion 4 del articulo)
|
|
# ----------------------------------------------------------------------------
|
|
def window_metrics(df, label):
|
|
"""Metricas de una ventana del log (df ya enriquecido con P_ref)."""
|
|
w = df[df["G_POA"] > 50.0] if len(df) else df
|
|
if w.empty:
|
|
return {"label": label, "n_samples": 0}
|
|
# segundos desde el inicio de la ventana (independiente de la resolucion
|
|
# de datetime de pandas 3.x)
|
|
ts = w["timestamp"].values
|
|
x = (ts - ts[0]) / np.timedelta64(1, "s")
|
|
E = _trapz(w["P_CPV"].values, x) / 3600.0 # Wh
|
|
Pref = _trapz(w["P_ref"].values, x) # W.s
|
|
H = _trapz(w["G_POA"].values, x) / 3600.0 / 1000.0 # kWh/m2
|
|
eta = 100.0 * E * 3600.0 / Pref if Pref > 0 else np.nan
|
|
S_el, S_az = int(w["dx"].sum()), int(w["dy"].sum())
|
|
return {
|
|
"label": label,
|
|
"n_samples": int(len(w)),
|
|
"t_inicio": str(w["timestamp"].iloc[0]),
|
|
"t_fin": str(w["timestamp"].iloc[-1]),
|
|
"E_Wh": round(E, 3),
|
|
"H_kWh_m2": round(H, 4),
|
|
"E_norm_Wh_per_kWhm2": round(E / H, 3) if H > 0 else None,
|
|
"eta_tr_pct": round(eta, 2) if eta == eta else None,
|
|
"MAE_I_mA": round(1000.0 * np.mean(np.abs(w["Isc_ref"] - w["Isc"])), 2),
|
|
"MAE_P_mW": round(1000.0 * np.mean(np.abs(w["P_ref"] - w["P_CPV"])), 2),
|
|
"P_CPV_mean_mW": round(1000.0 * np.mean(w["P_CPV"]), 1),
|
|
"P_ref_mean_mW": round(1000.0 * np.mean(w["P_ref"]), 1),
|
|
"S_el_usteps": S_el,
|
|
"S_az_usteps": S_az,
|
|
"S_tot_usteps": S_el + S_az,
|
|
"S_tot_deg": round((S_el + S_az) * DEG_PER_USTEP, 1),
|
|
}
|
|
|
|
|
|
def recovery_t95(df, df_e, horizon_min=20):
|
|
"""t95 tras cada comando de re-apuntado/re-centrado.
|
|
|
|
Se define la recuperacion como el primer instante, en una ventana de
|
|
`horizon_min` minutos despues del trigger, en que P_CPV alcanza de nuevo
|
|
el 95% de P_ref. Si no se alcanza dentro de la ventana se reporta None
|
|
(y, como dato informativo, el instante en que si se alcanza).
|
|
"""
|
|
out = []
|
|
cmds = df_e[df_e["kind"] == "cmd"] if not df_e.empty else pd.DataFrame()
|
|
if cmds.empty:
|
|
return pd.DataFrame(out)
|
|
for _, row in cmds.iterrows():
|
|
n, t0 = int(row["cmd"]), row["timestamp"]
|
|
if n not in (2, 5):
|
|
continue
|
|
w = df[(df["timestamp"] >= t0) &
|
|
(df["timestamp"] <= t0 + pd.Timedelta(minutes=horizon_min))].reset_index(drop=True)
|
|
if w.empty:
|
|
continue
|
|
thr = 0.95 * w["P_ref"].values
|
|
P = w["P_CPV"].values
|
|
below = np.where(P < thr)[0]
|
|
if below.size == 0:
|
|
out.append({"trigger": CMD_NAMES.get(n, n), "t_trigger": str(t0),
|
|
"t95_s": 0.0, "dropped": False, "P_min_mW": None})
|
|
continue
|
|
i_drop = int(below[0])
|
|
P_min = float(P[i_drop:].min())
|
|
t95 = None
|
|
above = np.where(P[i_drop:] >= thr)[0]
|
|
if above.size:
|
|
t95 = (w["timestamp"].iloc[i_drop + int(above[0])] - t0).total_seconds()
|
|
# cruce posterior (fuera de la ventana), solo informativo
|
|
t_late = None
|
|
rest = df[df["timestamp"] > t0 + pd.Timedelta(minutes=horizon_min)]
|
|
if t95 is None and not rest.empty:
|
|
above_late = np.where(rest["P_CPV"].values >= 0.95 * rest["P_ref"].values)[0]
|
|
if above_late.size:
|
|
t_late = (rest["timestamp"].iloc[int(above_late[0])] - t0).total_seconds()
|
|
out.append({"trigger": CMD_NAMES.get(n, n), "t_trigger": str(t0),
|
|
"t95_s": round(t95, 1) if t95 is not None else None,
|
|
"t95_late_s": round(t_late, 1) if t_late is not None else None,
|
|
"dropped": True, "P_min_mW": round(P_min * 1000, 1)})
|
|
return pd.DataFrame(out)
|
|
|
|
|
|
def cloud_events(df):
|
|
"""Transitorios nubosos: G < 0.85 * mediana movil 10 min (ventanas >= 120 s)."""
|
|
if df.empty:
|
|
return []
|
|
d = df.set_index("timestamp").sort_index()
|
|
g = d["G_POA"]
|
|
if g.max() < 300:
|
|
return []
|
|
med = g.rolling("10min", center=True, min_periods=10).median()
|
|
mask = (g < 0.85 * med) & (g > 30)
|
|
events = []
|
|
t = d.index
|
|
blocks, start = [], None
|
|
for i in range(len(mask)):
|
|
if mask.iloc[i] and start is None:
|
|
start = i
|
|
elif not mask.iloc[i] and start is not None:
|
|
blocks.append((start, i - 1))
|
|
start = None
|
|
if start is not None:
|
|
blocks.append((start, len(mask) - 1))
|
|
for i0, i1 in blocks:
|
|
dur = (t[i1] - t[i0]).total_seconds()
|
|
if dur < 120:
|
|
continue
|
|
g_seg = g.iloc[i0:i1 + 1]
|
|
base = float(g.iloc[max(0, i0 - 5):i0].mean())
|
|
g_min = float(g_seg.min())
|
|
rec_idx = np.where(g.iloc[i1:].values >= 0.95 * base)[0]
|
|
t_rec = (t[i1 + int(rec_idx[0])] - t[i0]).total_seconds() if rec_idx.size else None
|
|
events.append({
|
|
"t_start": str(t[i0]), "dur_s": int(dur),
|
|
"G_base_Wm2": round(base, 0), "G_min_Wm2": round(g_min, 0),
|
|
"depth_pct": round(100.0 * (1 - g_min / base), 1) if base > 0 else None,
|
|
"t95_G_s": round(t_rec, 1) if t_rec is not None else None,
|
|
})
|
|
return events
|
|
|
|
|
|
def position_errors(df, df_e):
|
|
"""Error de posicion vs objetivo astronomico tras cada reticle.
|
|
|
|
'inicial': primera muestra despues de la orden (fin del recorrido).
|
|
'convergido': mediana del offset durante los ultimos 5 min de los
|
|
primeros 10 min tras el reticle (periodo de convergencia del fine-tune;
|
|
despues, el objetivo geometrico queda obsoleto porque el sol sigue
|
|
moviendose y el controlador persigue el pico electrico, no el target).
|
|
"""
|
|
out = []
|
|
if df_e.empty:
|
|
return pd.DataFrame(out)
|
|
sp = df_e[df_e["kind"] == "sunpos"]
|
|
for i, (_, row) in enumerate(sp.iterrows()):
|
|
t0, tx, ty = row["timestamp"], row["px"], row["py"]
|
|
seg = df[(df["timestamp"] >= t0) & (df["status"] == 3)]
|
|
if seg.empty:
|
|
continue
|
|
first = seg.iloc[0]
|
|
conv = seg[(seg["timestamp"] >= t0 + pd.Timedelta(minutes=5)) &
|
|
(seg["timestamp"] <= t0 + pd.Timedelta(minutes=10))]
|
|
if conv.empty:
|
|
conv = seg.tail(50)
|
|
ex0, ey0 = abs(first["x"] - tx) * DEG_PER_USTEP, abs(first["y"] - ty) * DEG_PER_USTEP
|
|
exs = float(np.median([abs(r["x"] - tx) * DEG_PER_USTEP for _, r in conv.iterrows()]))
|
|
eys = float(np.median([abs(r["y"] - ty) * DEG_PER_USTEP for _, r in conv.iterrows()]))
|
|
out.append({
|
|
"t_reticle": str(t0), "az_corr": row["azc"], "alt_corr": row["altc"],
|
|
"errX0_deg": round(ex0, 2), "errY0_deg": round(ey0, 2),
|
|
"errX0_usteps": int(abs(first["x"] - tx)), "errY0_usteps": int(abs(first["y"] - ty)),
|
|
"errX_ss_deg": round(exs, 3), "errY_ss_deg": round(eys, 3),
|
|
})
|
|
return pd.DataFrame(out)
|
|
|
|
|
|
|
|
# ----------------------------------------------------------------------------
|
|
# Figuras (PDF vectoriales para el articulo)
|
|
# ----------------------------------------------------------------------------
|
|
CMD_COLORS = {2: "tab:green", 3: "tab:purple", 4: "tab:blue", 5: "tab:brown"}
|
|
|
|
|
|
def make_figures(df, df_e, outdir, tag):
|
|
if not HAS_MPL:
|
|
print("[aviso] matplotlib no disponible; se omiten figuras.")
|
|
return
|
|
figs_dir = ARTICULO / "figs"
|
|
figs_dir.mkdir(parents=True, exist_ok=True)
|
|
t = df["timestamp"]
|
|
lit_mask = (df["G_POA"] > 50.0).values
|
|
cmds = df_e[df_e["kind"] == "cmd"] if not df_e.empty else pd.DataFrame()
|
|
|
|
# --- Fig 1: potencia generada vs referencia + irradiancia
|
|
fig, ax1 = plt.subplots(figsize=(7.2, 3.4), dpi=150)
|
|
ax1.plot(t, df["P_CPV"] * 1000, color="tab:red", lw=1.0, label="$P_{CPV}$")
|
|
ax1.plot(t, df["P_ref"] * 1000, color="k", ls="--", lw=1.0, label="$P_{ref}$")
|
|
ax1.fill_between(t, df["P_CPV"] * 1000, df["P_ref"] * 1000,
|
|
where=lit_mask, color="tab:red", alpha=0.12,
|
|
label="Loss vs. $P_{ref}$")
|
|
ax1.set_ylabel("Power (mW)", color="tab:red")
|
|
ax1.tick_params(axis="y", labelcolor="tab:red")
|
|
ax2 = ax1.twinx()
|
|
ax2.plot(t, df["G_POA"], color="tab:orange", lw=0.8, alpha=0.85, label="$G_{POA}$")
|
|
ax2.set_ylabel("Irradiance (W/m$^2$)", color="tab:orange")
|
|
ax2.tick_params(axis="y", labelcolor="tab:orange")
|
|
labeled = False
|
|
for _, r in cmds.iterrows():
|
|
if int(r["cmd"]) != 2:
|
|
continue # solo se anota "Track solar": los demas comandos
|
|
# estan todos dentro de ~75 s entre si y sus etiquetas se solaparian
|
|
ax1.axvline(r["timestamp"], color=CMD_COLORS.get(r["cmd"], "gray"),
|
|
ls=":", lw=1.0, alpha=0.8)
|
|
if not labeled: # etiqueta unica: ambos Track-solar estan a < 2 min
|
|
ax1.text(r["timestamp"], ax1.get_ylim()[1] * 0.95,
|
|
f"cmd {int(r['cmd'])}\n({CMD_NAMES.get(int(r['cmd']), '')})",
|
|
fontsize=6.5, color=CMD_COLORS.get(r["cmd"], "gray"))
|
|
labeled = True
|
|
h1, l1 = ax1.get_legend_handles_labels()
|
|
h2, l2 = ax2.get_legend_handles_labels()
|
|
ax1.legend(h1 + h2, l1 + l2, loc="upper right", fontsize=7, framealpha=0.6)
|
|
ax1.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
|
|
ax1.set_xlabel("Local time (21 Sep 2025, UTC-6)")
|
|
ax1.set_title(f"WiT-CPV field test - session {tag}", fontsize=9)
|
|
fig.tight_layout()
|
|
fig.savefig(figs_dir / f"fig_power_{tag}.pdf")
|
|
plt.close(fig)
|
|
|
|
# --- Fig 2: posicion real vs objetivo astronomico
|
|
fig, axes = plt.subplots(2, 1, figsize=(7.2, 4.2), dpi=150, sharex=True)
|
|
for ax, col, name, color in ((axes[0], "x", "$X$ (zenithal, $\\mu$steps)", "tab:blue"),
|
|
(axes[1], "y", "$Y$ (azimuthal, $\\mu$steps)", "tab:olive")):
|
|
ax.plot(t, df[col], color=color, lw=0.9, label="Actual position")
|
|
tc = df["tgt_x" if col == "x" else "tgt_y"]
|
|
ax.step(t, tc, where="post", color="k", ls="--", lw=1.0,
|
|
label="Predicted sun position")
|
|
ax.set_ylabel(name)
|
|
if not cmds.empty:
|
|
for _, r in cmds.iterrows():
|
|
ax.axvline(r["timestamp"], color=CMD_COLORS.get(r["cmd"], "gray"),
|
|
ls=":", lw=0.8, alpha=0.7)
|
|
ax.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
|
|
ax.tick_params(labelsize=7)
|
|
axes[1].set_xlabel("Local time (21 Sep 2025, UTC-6)")
|
|
fig.tight_layout()
|
|
fig.savefig(figs_dir / f"fig_position_{tag}.pdf")
|
|
plt.close(fig)
|
|
|
|
|
|
|
|
# --- Fig 3: fine-tune adaptativo (e_dyn, |dC|)
|
|
fig, ax1 = plt.subplots(figsize=(7.2, 3.2), dpi=150)
|
|
ax1.plot(t, df["er"], color="tab:purple", lw=0.7, alpha=0.9,
|
|
label="Dynamic threshold $e_{dyn}$")
|
|
ax2 = ax1.twinx()
|
|
ax2.semilogy(t, df["dC"].clip(lower=1), color="tab:red", lw=0.5, alpha=0.7,
|
|
label="$|\\Delta I_{sc}|$")
|
|
ax2.set_ylabel("$|\\Delta I_{sc}|$ (ADC counts, log)", color="tab:red")
|
|
ax2.tick_params(axis="y", labelcolor="tab:red")
|
|
ax1.set_ylabel("$e_{dyn}$ ($\\mu$steps)")
|
|
at = df_e[df_e["kind"] == "at_rec"] if not df_e.empty else pd.DataFrame()
|
|
for _, r in at.iterrows():
|
|
ax1.axvline(r["timestamp"], color="k", ls=":", lw=1.0)
|
|
ax1.text(r["timestamp"], ax1.get_ylim()[1] * 0.9, "auto-tune\nrecommendation",
|
|
fontsize=6.5, ha="left")
|
|
h1, l1 = ax1.get_legend_handles_labels()
|
|
h2, l2 = ax2.get_legend_handles_labels()
|
|
ax1.legend(h1 + h2, l1 + l2, loc="upper right", fontsize=7)
|
|
ax1.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
|
|
ax1.set_xlabel("Local time (21 Sep 2025, UTC-6)")
|
|
fig.tight_layout()
|
|
fig.savefig(figs_dir / f"fig_finetune_{tag}.pdf")
|
|
plt.close(fig)
|
|
|
|
# --- Fig 4: recuperacion t95 alrededor de los triggers
|
|
rec = [r for _, r in cmds.iterrows() if int(r["cmd"]) in (2, 5)]
|
|
if rec:
|
|
fig, axes = plt.subplots(1, min(2, len(rec)), figsize=(7.2, 2.9), dpi=150)
|
|
axes = np.atleast_1d(axes)
|
|
for ax, r in zip(axes, rec[:2]):
|
|
t0 = r["timestamp"]
|
|
w = df[(df["timestamp"] >= t0 - pd.Timedelta(minutes=1)) &
|
|
(df["timestamp"] <= t0 + pd.Timedelta(minutes=20))]
|
|
ax.plot(w["timestamp"], w["P_CPV"] * 1000, color="tab:red", lw=1.0)
|
|
ax.plot(w["timestamp"], w["P_ref"] * 1000, color="k", ls="--", lw=1.0)
|
|
ax.plot(w["timestamp"], 0.95 * w["P_ref"] * 1000, color="tab:green",
|
|
ls=":", lw=1.0, label="0.95 $P_{ref}$")
|
|
ax.axvline(t0, color="tab:blue", ls=":", lw=1.0)
|
|
thr = 0.95 * w["P_ref"].values
|
|
P = w["P_CPV"].values
|
|
below = np.where(P < thr)[0]
|
|
if below.size:
|
|
i = int(below[0])
|
|
above = np.where(P[i:] >= thr)[0]
|
|
if above.size:
|
|
t95s = (w["timestamp"].iloc[i + int(above[0])] -
|
|
w["timestamp"].iloc[i]).total_seconds()
|
|
ax.annotate(f"$t_{{95}}$ = {t95s:.0f} s",
|
|
xy=(w["timestamp"].iloc[i + int(above[0])],
|
|
P[i + int(above[0])] * 1000),
|
|
fontsize=8, color="tab:green")
|
|
ax.set_title(f"cmd {int(r['cmd'])} ({CMD_NAMES.get(int(r['cmd']), '')}) "
|
|
f"@ {str(r['timestamp'])[11:16]}", fontsize=8)
|
|
ax.set_ylabel("Power (mW)")
|
|
ax.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M:%S"))
|
|
ax.tick_params(labelsize=7)
|
|
ax.legend(fontsize=7)
|
|
fig.tight_layout()
|
|
fig.savefig(figs_dir / f"fig_recovery_{tag}.pdf")
|
|
plt.close(fig)
|
|
print(f"[figs] guardadas en {figs_dir}")
|
|
|
|
|
|
|
|
# ----------------------------------------------------------------------------
|
|
# Main
|
|
# ----------------------------------------------------------------------------
|
|
def main():
|
|
args = sys.argv[1:]
|
|
log_file = args[0] if args else "21092025.txt"
|
|
tag = Path(log_file).stem
|
|
path = ROOT / log_file
|
|
outdir = ARTICULO / f"resultados_{tag}"
|
|
outdir.mkdir(parents=True, exist_ok=True)
|
|
|
|
print(f"[1/5] Parseando {path} ...")
|
|
df_s, df_e = parse_log(path)
|
|
print(f" muestras={len(df_s)} eventos={len(df_e)} "
|
|
f"({df_s['timestamp'].iloc[0]} -> {df_s['timestamp'].iloc[-1]})")
|
|
|
|
print("[2/5] Unidades fisicas + estados + P_ref ...")
|
|
df, ph = build_physical(df_s, df_e)
|
|
lit = df[df["G_POA"] > 200.0]
|
|
R_opt = df["R_opt"].iloc[0]
|
|
print(f" R_opt,max = {1000 * R_opt:.3f} mA/(W/m2) "
|
|
f"(ventana iluminada: {lit['timestamp'].iloc[0]} -> {lit['timestamp'].iloc[-1]})")
|
|
|
|
print("[3/5] Metricas ...")
|
|
phases = ph[ph["status"] == 3]
|
|
metrics = {"tag": tag,
|
|
"log": path.name,
|
|
"session_start": str(df["timestamp"].iloc[0]),
|
|
"session_end": str(df["timestamp"].iloc[-1]),
|
|
"R_opt_max_mA_per_Wm2": round(1000 * R_opt, 3),
|
|
"FF_used": FF}
|
|
metrics["full_session"] = window_metrics(df, "Full session (G > 50 W/m2)")
|
|
metrics["daylight"] = window_metrics(lit, "Daylight (G > 200 W/m2)")
|
|
metrics["tracking_active"] = window_metrics(df[df["status"] == 3],
|
|
"Tracking active (status=3)")
|
|
clear = df[df["G_POA"] > 600.0]
|
|
if len(clear):
|
|
m_clear = window_metrics(clear, "Clear-sky (G > 600 W/m2)")
|
|
m_clear["frac_P_above_095Pref"] = round(
|
|
float((clear["P_CPV"] >= 0.95 * clear["P_ref"]).mean()), 3)
|
|
m_clear["ratio_P_Pref_median"] = round(
|
|
float((clear["P_CPV"] / clear["P_ref"]).median()), 3)
|
|
metrics["clear_sky"] = m_clear
|
|
else:
|
|
metrics["clear_sky"] = None
|
|
metrics["phases"] = [window_metrics(df[df["phase"] == p], f"Phase: {p}")
|
|
for p in phases["phase"].drop_duplicates()]
|
|
|
|
rec = recovery_t95(df, df_e)
|
|
metrics["t95_recovery"] = rec.to_dict("records") if not rec.empty else []
|
|
metrics["cloud_events"] = cloud_events(df)
|
|
pe = position_errors(df, df_e)
|
|
metrics["position_errors"] = pe.to_dict("records") if not pe.empty else []
|
|
er_seq = df["er"].dropna().astype(int)
|
|
changes = er_seq[er_seq != er_seq.shift(1)]
|
|
changes = changes.iloc[1:] if changes.index[0] == df.index[0] else changes
|
|
metrics["er_progression"] = [
|
|
{"t": str(ts), "e": int(ev)}
|
|
for ts, ev in zip(df["timestamp"].values[changes.index], changes.values)]
|
|
if not df_e.empty:
|
|
at = df_e[df_e["kind"] == "at_rec"]
|
|
cols = [c for c in ("timestamp", "small", "abrupt", "reset", "peak")
|
|
if c in at.columns]
|
|
metrics["at_recommendations"] = at[cols].to_dict("records")
|
|
else:
|
|
metrics["at_recommendations"] = []
|
|
|
|
print("[4/5] Guardando salida ...")
|
|
keep = ["timestamp", "x", "y", "c_raw", "v_raw", "p_raw", "s_uptime", "er",
|
|
"G_POA", "Voc", "Isc", "P_CPV", "Isc_ref", "P_ref", "status", "phase",
|
|
"az_est", "alt_est", "tgt_x", "tgt_y", "az", "alt", "azc", "altc",
|
|
"azdev_used", "e", "dC", "dynE", "near", "cmax"]
|
|
df[keep].to_csv(outdir / "samples.csv", index=False)
|
|
df_e.to_csv(outdir / "events.csv", index=False, na_rep="")
|
|
with open(outdir / "metrics.json", "w", encoding="utf-8") as f:
|
|
json.dump(metrics, f, indent=2, ensure_ascii=False, default=str)
|
|
|
|
print("[5/5] Figuras ...")
|
|
make_figures(df, df_e, outdir, tag)
|
|
|
|
# -------- resumen en consola --------
|
|
print("\n" + "=" * 68)
|
|
print(f" RESUMEN DE METRICAS (log {path.name})")
|
|
print("=" * 68)
|
|
for key in ("full_session", "daylight", "tracking_active", "clear_sky"):
|
|
m = metrics[key]
|
|
if m is None:
|
|
continue
|
|
extra = (f" | %P>=0.95Pref = {m['frac_P_above_095Pref']:.0%} | "
|
|
f"med(P/Pref) = {m['ratio_P_Pref_median']}"
|
|
if "frac_P_above_095Pref" in m else "")
|
|
print(f"\n--- {m['label']} ({m['t_inicio'][11:19]} -> {m['t_fin'][11:19]}, "
|
|
f"n={m['n_samples']})")
|
|
print(f" E = {m['E_Wh']:.2f} Wh | H = {m['H_kWh_m2']:.4f} kWh/m2 | "
|
|
f"E/H = {m['E_norm_Wh_per_kWhm2']} Wh/(kWh/m2)")
|
|
print(f" eta_tr = {m['eta_tr_pct']} % | MAE_I = {m['MAE_I_mA']} mA | "
|
|
f"MAE_P = {m['MAE_P_mW']} mW{extra}")
|
|
print(f" S_el = {m['S_el_usteps']:,} | S_az = {m['S_az_usteps']:,} | "
|
|
f"S_tot = {m['S_tot_usteps']:,} usteps ({m['S_tot_deg']} deg)")
|
|
for m in metrics["phases"]:
|
|
print(f"\n [{m['label']}] {m['t_inicio'][11:16]} -> {m['t_fin'][11:16]} | "
|
|
f"eta={m['eta_tr_pct']}% E={m['E_Wh']}Wh H={m['H_kWh_m2']}kWh/m2 "
|
|
f"S={m['S_tot_usteps']:,}")
|
|
if metrics["t95_recovery"]:
|
|
print("\n--- t95 (recuperacion a 95% de P_ref tras trigger, horizonte 20 min)")
|
|
for r in metrics["t95_recovery"]:
|
|
extra = (f" (P_min = {r['P_min_mW']} mW)"
|
|
if r.get("P_min_mW") is not None else "")
|
|
late = (f" [primer cruce posterior: {r['t95_late_s']} s despues]"
|
|
if r.get("t95_late_s") else "")
|
|
print(f" {r['trigger']:>12} @ {r['t_trigger'][11:19]} -> t95 = "
|
|
f"{r['t95_s']} s{extra}{late}")
|
|
if metrics["cloud_events"]:
|
|
print(f"\n--- {len(metrics['cloud_events'])} transitorios nubosos detectados")
|
|
for c in metrics["cloud_events"][:8]:
|
|
print(f" {c['t_start'][11:19]} G: {c['G_base_Wm2']:.0f} -> "
|
|
f"{c['G_min_Wm2']:.0f} W/m2 ({c['depth_pct']}%) "
|
|
f"dur={c['dur_s']}s t95_G={c['t95_G_s']}s")
|
|
if metrics["position_errors"]:
|
|
print("\n--- Error de posicion vs objetivo astronomico")
|
|
for p in metrics["position_errors"]:
|
|
print(f" {p['t_reticle'][11:19]} inicial: {p['errX0_deg']}/"
|
|
f"{p['errY0_deg']} deg ({p['errX0_usteps']}/{p['errY0_usteps']} "
|
|
f"usteps) | regimen: {p['errX_ss_deg']}/{p['errY_ss_deg']} deg")
|
|
if metrics.get("er_progression"):
|
|
print("\n--- Progresion del umbral adaptativo e_dyn (cambios)")
|
|
print(" " + " -> ".join(f"{r['t'][11:19]}:{r['e']}"
|
|
for r in metrics["er_progression"]))
|
|
if metrics.get("at_recommendations"):
|
|
for a in metrics["at_recommendations"]:
|
|
print(f"\n--- Auto-tune @ {str(a['timestamp'])[11:19]}: "
|
|
f"small={a['small']:.0f}, abrupt={a['abrupt']:.0f}, "
|
|
f"reset={a['reset']:.0f}, peak_prox={a['peak']}")
|
|
print(f"\n[ok] salida en {outdir}\n")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|