528 lines
22 KiB
Python
528 lines
22 KiB
Python
# -*- coding: utf-8 -*-
|
|
"""Trayectoria optima (X, Y) de la celda por correlacion G_POA ~ Isc.
|
|
|
|
Metodo (ver conversacion 2026-10-18):
|
|
1. Base de datos: todos los logs de logs_nuevos/ con dia (G > 20 W/m2),
|
|
ventana robusta P0.5-P99.5 para descartar outliers de la trayectoria.
|
|
2. Normalizacion solar: rho = Isc / Isc_max(dia). Validacion del proxy:
|
|
Pearson r y Spearman rho de Isc vs G por dia (r >= 0.7 => G es proxy).
|
|
3. Cielo claro: indice de claridad CI = integral(G) / integral(G_clear),
|
|
donde G_clear es la envolvente superior (max por bin de 10 min de hora
|
|
solar) entre todos los dias. Definen la trayectoria los dias con
|
|
CI >= 0.85*CI_max (top 3 por CI); con envolvente patchwork de pocos
|
|
dias el CI integral no alcanza 0.8 absoluto.
|
|
4. Tabla operativa T: h -> (X*, Y*). Por bin de 15 min de hora solar,
|
|
mediana de X e Y de los dias claros; relleno de huecos + suavizado.
|
|
5. Ruta espacial: Dijkstra sobre la rejilla (costo = -score, 8 vecinos)
|
|
desde el primer punto de la tabla hasta el ultimo.
|
|
6. Backtest por dia (incluye nublados): E_real, E_open (tabla en bucle
|
|
abierto, muestreo bilineal de la rejilla del dia) y E_hc (tabla +
|
|
hill-climbing de 1 paso = max del bloque 3x3). Para los dias claros la
|
|
tabla se construye SIN el dia evaluado (leave-one-out); los nublados
|
|
se evaluan contra la tabla global (referencia). Metrica: dE %.
|
|
7. Salida: figura de 3 paneles (mapa + rutas, barras dE %, curva de la
|
|
tabla) y CSV: tabla_trayectoria.csv (10 min) y resumen_backtest.csv.
|
|
|
|
Mismo mapeo y calibracion que heatmap_extremos_p80.py (adc_manager.cpp:
|
|
C=canal1 Isc, P=canal2 piranometro on-board, V=canal3 Voc; ADS122C04
|
|
24-bit, Vref 2.048 V, canal P en complemento a 2).
|
|
|
|
Uso: py -3.12 trayectoria_optima.py
|
|
"""
|
|
|
|
import glob
|
|
import os
|
|
import re
|
|
from datetime import datetime
|
|
|
|
import numpy as np
|
|
import matplotlib
|
|
matplotlib.use("Agg")
|
|
import matplotlib.pyplot as plt
|
|
from matplotlib.colors import Normalize
|
|
from matplotlib.patches import FancyArrowPatch
|
|
|
|
ROOT = os.path.dirname(os.path.abspath(__file__))
|
|
LOGS_DIR = os.path.join(ROOT, "logs_nuevos")
|
|
OUT_DIR = os.path.join(ROOT, "figures", "trayectoria_optima")
|
|
|
|
ADC_MAX = 8388607.0
|
|
GAIN_SHUNT = 20.0
|
|
R_SHUNT = 0.075
|
|
DEG_PER_USTEP = 80.0 / 53000.0
|
|
|
|
DAY_G_MIN = 20.0
|
|
PAD_PCT = 0.5
|
|
T0, T1 = 7.5, 17.5 # ventana de hora solar de la tabla
|
|
BIN_TAB = 0.25 # bin de la tabla (h)
|
|
TAB_STEP = 0.1 # paso del CSV (h)
|
|
NX, NY = 61, 46 # rejilla comun
|
|
|
|
RE = re.compile(
|
|
r"LOG:\[(?P<t>[^\]]+)\]\s*"
|
|
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+)"
|
|
)
|
|
|
|
|
|
def trapz(y, x):
|
|
f = getattr(np, "trapezoid", None)
|
|
if f is None:
|
|
f = np.trapz
|
|
return f(y, x)
|
|
|
|
|
|
def sign24(raw):
|
|
raw &= 0xFFFFFF
|
|
return raw - 0x1000000 if raw & 0x800000 else raw
|
|
|
|
|
|
def isc_ma(c_raw):
|
|
vsh_mV = c_raw / ADC_MAX * 2048.0 / GAIN_SHUNT
|
|
return (vsh_mV / 1000.0) / R_SHUNT * 1000.0
|
|
|
|
|
|
def robust_limits(arr, lo=PAD_PCT, hi=100.0 - PAD_PCT):
|
|
a = np.asarray(arr, dtype=float)
|
|
return float(np.percentile(a, lo)), float(np.percentile(a, hi))
|
|
|
|
|
|
def load_log(path):
|
|
ts, xs, ys, cs, ps = [], [], [], [], []
|
|
with open(path, "r", encoding="utf-8", errors="replace") as f:
|
|
for line in f:
|
|
m = RE.search(line)
|
|
if m:
|
|
ts.append(m.group("t"))
|
|
xs.append(int(m.group("x")))
|
|
ys.append(int(m.group("y")))
|
|
cs.append(int(m.group("c")))
|
|
ps.append(sign24(int(m.group("p"))))
|
|
return (np.array(ts), np.array(xs, dtype=float), np.array(ys, dtype=float),
|
|
np.array(cs, dtype=float), np.array(ps, dtype=float))
|
|
|
|
|
|
def parse_t(t_strs):
|
|
out = np.empty(len(t_strs), dtype=float)
|
|
for i, s in enumerate(t_strs):
|
|
dt = datetime.strptime(s, "%Y-%m-%d %H:%M:%S")
|
|
out[i] = dt.hour + dt.minute / 60.0 + dt.second / 3600.0
|
|
return out
|
|
|
|
|
|
def bin_stat(x, y, vals, xbins, ybins):
|
|
xi = np.clip(np.searchsorted(xbins, x, side="right") - 1, 0, len(xbins) - 2)
|
|
yi = np.clip(np.searchsorted(ybins, y, side="right") - 1, 0, len(ybins) - 2)
|
|
shape = (len(ybins) - 1, len(xbins) - 1)
|
|
flat_sum = np.zeros(shape[0] * shape[1])
|
|
flat_cnt = np.zeros(shape[0] * shape[1])
|
|
for a, b, v in zip(yi, xi, vals):
|
|
k = a * shape[1] + b
|
|
flat_sum[k] += v
|
|
flat_cnt[k] += 1.0
|
|
with np.errstate(invalid="ignore"):
|
|
flat = flat_sum / np.where(flat_cnt == 0, np.nan, flat_cnt)
|
|
return flat.reshape(shape)
|
|
|
|
|
|
def bilinear(grid, x, y, xbins, ybins):
|
|
"""Muestreo bilineal de la rejilla en posiciones (x, y) arbitraras."""
|
|
nx, ny = grid.shape[1], grid.shape[0]
|
|
fx = (np.asarray(x, dtype=float) - xbins[0]) / (xbins[-1] - xbins[0]) * nx
|
|
fy = (np.asarray(y, dtype=float) - ybins[0]) / (ybins[-1] - ybins[0]) * ny
|
|
fx = np.clip(fx, 0, nx - 1)
|
|
fy = np.clip(fy, 0, ny - 1)
|
|
j0 = np.clip(fx.astype(int), 0, nx - 2)
|
|
i0 = np.clip(fy.astype(int), 0, ny - 2)
|
|
wx = fx - j0
|
|
wy = fy - i0
|
|
v = (grid[i0, j0] * (1 - wx) * (1 - wy) + grid[i0, j0 + 1] * wx * (1 - wy)
|
|
+ grid[i0 + 1, j0] * (1 - wx) * wy + grid[i0 + 1, j0 + 1] * wx * wy)
|
|
return v
|
|
|
|
|
|
def nan_max_pool3(grid):
|
|
"""Maximo del bloque 3x3 (nan-aware): aproxima un hill-climbing de 1 paso."""
|
|
import warnings
|
|
p = np.pad(grid, 1, mode="edge")
|
|
n, m = grid.shape
|
|
with np.errstate(invalid="ignore"), warnings.catch_warnings():
|
|
warnings.simplefilter("ignore", RuntimeWarning)
|
|
out = np.nanmax(
|
|
np.stack([p[di:di + n, dj:dj + m] for di in range(3) for dj in range(3)]),
|
|
axis=0)
|
|
return out
|
|
|
|
|
|
def dijkstra(grid, i0, j0, i1, j1):
|
|
"""Ruta de maximo score (costo = -score) con 8 vecinos en la rejilla."""
|
|
import heapq
|
|
n, m = grid.shape
|
|
INF = np.inf
|
|
best = np.full((n, m), INF)
|
|
prev = np.full((n, m), -1, dtype=int)
|
|
best[i0, j0] = 0.0
|
|
pq = [(0.0, i0, j0)]
|
|
nbrs = [(-1, -1, 1.0), (-1, 0, 1.414), (-1, 1, 1.414), (0, -1, 1.0),
|
|
(0, 1, 1.0), (1, -1, 1.414), (1, 0, 1.0), (1, 1, 1.414)]
|
|
while pq:
|
|
d, i, j = heapq.heappop(pq)
|
|
if d > best[i, j] + 1e-12:
|
|
continue
|
|
if i == i1 and j == j1:
|
|
break
|
|
for di, dj, w in nbrs:
|
|
ni, nj = i + di, j + dj
|
|
if ni < 0 or ni >= n or nj < 0 or nj >= m:
|
|
continue
|
|
if not np.isfinite(grid[ni, nj]):
|
|
continue
|
|
nd = d + w * (0.0 if grid[ni, nj] >= 1 else 1.0 - grid[ni, nj])
|
|
if nd < best[ni, nj] - 1e-12:
|
|
best[ni, nj] = nd
|
|
prev[ni, nj] = i * m + j
|
|
heapq.heappush(pq, (nd, ni, nj))
|
|
if not np.isfinite(best[i1, j1]):
|
|
return None
|
|
path = []
|
|
cur = i1 * m + j1
|
|
while cur != -1:
|
|
i, j = divmod(cur, m)
|
|
path.append((i, j))
|
|
cur = prev[i, j]
|
|
path.reverse()
|
|
return path
|
|
|
|
|
|
def build_table(days_subset, tb):
|
|
"""Tabla T: centro de cada bin de tb -> (X*, Y*). Mediana por bin +
|
|
relleno de huecos + suavizado. Devuelve None si hay menos de 4 bins
|
|
con datos."""
|
|
tab_x = np.full(len(tb) - 1, np.nan)
|
|
tab_y = np.full(len(tb) - 1, np.nan)
|
|
for k in range(len(tb) - 1):
|
|
pts_x, pts_y = [], []
|
|
for d in days_subset:
|
|
m = (d["t_solar"] >= tb[k]) & (d["t_solar"] < tb[k + 1])
|
|
pts_x.extend(d["x"][m].tolist())
|
|
pts_y.extend(d["y"][m].tolist())
|
|
if len(pts_x) >= 3:
|
|
tab_x[k] = np.median(pts_x)
|
|
tab_y[k] = np.median(pts_y)
|
|
tf = (tb[:-1] + tb[1:]) / 2
|
|
fin = np.isfinite(tab_x)
|
|
if fin.sum() < 4:
|
|
return None
|
|
tab_x = np.interp(tf, tf[fin], tab_x[fin])
|
|
tab_y = np.interp(tf, tf[fin], tab_y[fin])
|
|
# suavizado con nucleo [1,2,3,2,1]/9
|
|
ker = np.array([1, 2, 3, 2, 1], dtype=float) / 9
|
|
tab_x = np.convolve(tab_x, ker, mode="same")
|
|
tab_y = np.convolve(tab_y, ker, mode="same")
|
|
return tf, tab_x, tab_y
|
|
|
|
|
|
def cell_index(x, y, xbins, ybins):
|
|
nx, ny = len(xbins) - 1, len(ybins) - 1
|
|
j = int(np.clip((x - xbins[0]) / (xbins[-1] - xbins[0]) * nx, 0, nx - 1))
|
|
i = int(np.clip((y - ybins[0]) / (ybins[-1] - ybins[0]) * ny, 0, ny - 1))
|
|
return i, j
|
|
|
|
|
|
def main():
|
|
os.makedirs(OUT_DIR, exist_ok=True)
|
|
paths = sorted(glob.glob(os.path.join(LOGS_DIR, "log_*.txt")),
|
|
key=lambda p: re.search(r"(\d{4}-\d{2}-\d{2})", p).group(1))
|
|
print("Logs disponibles:")
|
|
days = []
|
|
for p in paths:
|
|
ts, xs, ys, cs, ps = load_log(p)
|
|
g = np.maximum(0.0, ps) / ADC_MAX * 2048.0 / 2.0
|
|
isc = isc_ma(cs)
|
|
day = g > DAY_G_MIN
|
|
if day.sum() < 50:
|
|
print(" %-28s sin dia suficiente, se omite" % os.path.basename(p))
|
|
continue
|
|
x0, x1 = robust_limits(xs[day])
|
|
y0, y1 = robust_limits(ys[day])
|
|
ins = day & (xs >= x0) & (xs <= x1) & (ys >= y0) & (ys <= y1)
|
|
if ins.sum() < 50:
|
|
print(" %-28s se omite" % os.path.basename(p))
|
|
continue
|
|
t = parse_t(ts[ins])
|
|
# mediodia solar: maximo de G con media movil de 1 h
|
|
k = max(3, int(np.round(1.0 / max(1e-6, np.median(np.diff(t))))))
|
|
kwin = 2 * k + 1
|
|
kern = np.ones(kwin) / kwin
|
|
gg_pad = np.pad(g[ins], k, mode="edge")
|
|
gg_s = np.convolve(gg_pad, kern, mode="valid")[:len(g[ins])]
|
|
t_noon = t[int(np.argmax(gg_s))]
|
|
t_solar = t - t_noon + 12.0
|
|
days.append(dict(
|
|
name=os.path.splitext(os.path.basename(p))[0],
|
|
x=xs[ins], y=ys[ins], c=isc[ins], g=g[ins],
|
|
t=t, t_solar=t_solar,
|
|
x0=x0, x1=x1, y0=y0, y1=y1,
|
|
))
|
|
print(" %-28s n_dia=%d" % (os.path.basename(p), int(ins.sum())))
|
|
if len(days) < 2:
|
|
print("No hay suficientes dias de datos.")
|
|
return
|
|
|
|
# rejilla comun sobre la union de las ventanas de los dias
|
|
xlo, xhi = min(d["x0"] for d in days), max(d["x1"] for d in days)
|
|
ylo, yhi = min(d["y0"] for d in days), max(d["y1"] for d in days)
|
|
xbins = np.linspace(xlo, xhi, NX)
|
|
ybins = np.linspace(ylo, yhi, NY)
|
|
|
|
# --- paso 1: metricas por dia y rejillas -------------------------------
|
|
bins_clear = np.arange(7.0, 18.0 + 1e-9, 1.0 / 6.0) # bins de 10 min
|
|
clear_vals = np.full(len(bins_clear) - 1, np.nan)
|
|
for d in days:
|
|
d["rho"] = d["c"] / d["c"].max()
|
|
d["grid_rho"] = bin_stat(d["x"], d["y"], d["rho"], xbins, ybins)
|
|
d["grid_g"] = bin_stat(d["x"], d["y"], d["g"], xbins, ybins)
|
|
sel = (d["t_solar"] >= 7.0) & (d["t_solar"] <= 18.0)
|
|
for k in range(len(bins_clear) - 1):
|
|
m = sel & (d["t_solar"] >= bins_clear[k]) & (d["t_solar"] < bins_clear[k + 1])
|
|
if m.any():
|
|
clear_vals[k] = np.nanmax([clear_vals[k], d["g"][m].max()])
|
|
# interpolacion de huecos de la envolvente
|
|
fin = np.isfinite(clear_vals)
|
|
clear_vals = np.interp(np.arange(len(clear_vals)), np.where(fin)[0], clear_vals[fin])
|
|
tc = (bins_clear[:-1] + bins_clear[1:]) / 2
|
|
for d in days:
|
|
g_clear_d = np.interp(d["t_solar"], tc, clear_vals)
|
|
# CI por integral: dias con cielo claro capturan mas energia total
|
|
with np.errstate(invalid="ignore"):
|
|
d["CI"] = float(trapz(d["g"], d["t_solar"]) /
|
|
max(1e-9, trapz(g_clear_d, d["t_solar"])))
|
|
if d["c"].std() > 0 and d["g"].std() > 0:
|
|
d["r"] = float(np.corrcoef(d["c"], d["g"])[0, 1])
|
|
else:
|
|
d["r"] = np.nan
|
|
cr = np.argsort(np.argsort(d["c"]))
|
|
gr = np.argsort(np.argsort(d["g"]))
|
|
d["rho_s"] = float(np.corrcoef(cr, gr)[0, 1])
|
|
ci_max = max(d["CI"] for d in days)
|
|
# umbral relativo: la envolvente es un patchwork con pocos dias, el CI
|
|
# integral nunca llega a 0.8 absoluto; el dia mas claro queda ~0.7-0.9
|
|
clear = sorted([d for d in days if d["CI"] >= 0.85 * ci_max],
|
|
key=lambda d: d["CI"], reverse=True)[:3]
|
|
print("\nDias que definen la trayectoria (CI >= 0.85*CI_max = %.2f): %s"
|
|
% (0.85 * ci_max,
|
|
", ".join("%s(%.2f)" % (d["name"], d["CI"]) for d in clear)))
|
|
if not clear:
|
|
print("Ningun dia claro; se usa el de mayor CI.")
|
|
clear = [max(days, key=lambda d: d["CI"])]
|
|
|
|
# --- paso 2: tabla T: h -> (X*, Y*) ------------------------------------
|
|
tb = np.arange(T0, T1 + 1e-9, BIN_TAB)
|
|
res = build_table(clear, tb)
|
|
if res is None:
|
|
print("Tabla sin datos suficientes.")
|
|
return
|
|
tab_t, tab_x, tab_y = res
|
|
|
|
# tabla a paso de 10 min para el CSV
|
|
csv_t = np.arange(tab_t[0], tab_t[-1] + 1e-9, TAB_STEP)
|
|
csv_x = np.interp(csv_t, tab_t, tab_x)
|
|
csv_y = np.interp(csv_t, tab_t, tab_y)
|
|
csv_path = os.path.join(OUT_DIR, "tabla_trayectoria.csv")
|
|
with open(csv_path, "w", encoding="utf-8") as f:
|
|
f.write("hora_solar,X_ref_upaso,Y_ref_upaso,X_ref_deg,Y_ref_deg\n")
|
|
for a, b, c in zip(csv_t, csv_x, csv_y):
|
|
f.write("%.1f,%.0f,%.0f,%.2f,%.2f\n"
|
|
% (a, b, c, b * DEG_PER_USTEP, c * DEG_PER_USTEP))
|
|
|
|
# --- paso 3: rejilla pool y Dijkstra ------------------------------------
|
|
import warnings
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore", RuntimeWarning)
|
|
pool_rho = np.nanmean(np.stack([d["grid_rho"] for d in clear]), axis=0)
|
|
pool_g = np.nanmean(np.stack([d["grid_g"] for d in clear]), axis=0)
|
|
pool_rho_n = pool_rho / np.nanmax(pool_rho)
|
|
i0, j0 = cell_index(tab_x[0], tab_y[0], xbins, ybins)
|
|
i1, j1 = cell_index(tab_x[-1], tab_y[-1], xbins, ybins)
|
|
path = dijkstra(pool_rho_n, i0, j0, i1, j1)
|
|
if path:
|
|
pcx = (xbins[np.array([j for i, j in path])] +
|
|
xbins[np.array([j for i, j in path]) + 1]) / 2
|
|
pcy = (ybins[np.array([i for i, j in path])] +
|
|
ybins[np.array([i for i, j in path]) + 1]) / 2
|
|
else:
|
|
pcx, pcy = np.array([]), np.array([])
|
|
|
|
# --- paso 4: backtest por dia -------------------------------------------
|
|
rows = []
|
|
for d in days:
|
|
is_clear = any(x is d for x in clear)
|
|
g_hc = nan_max_pool3(d["grid_rho"])
|
|
m = (d["t_solar"] >= tab_t[0]) & (d["t_solar"] <= tab_t[-1])
|
|
tt = d["t_solar"][m]
|
|
e_real = float(trapz(d["rho"][m], tt)) if len(tt) > 1 else 0.0
|
|
if is_clear and len(clear) >= 2:
|
|
# leave-one-out: tabla construida sin el dia evaluado
|
|
r_loo = build_table([x for x in clear if x is not d], tb)
|
|
if r_loo is None:
|
|
r_loo = (tab_t, tab_x, tab_y)
|
|
tipo = "self-fit"
|
|
else:
|
|
tipo = "LOO"
|
|
elif is_clear:
|
|
r_loo = (tab_t, tab_x, tab_y)
|
|
tipo = "self-fit"
|
|
else:
|
|
r_loo = (tab_t, tab_x, tab_y)
|
|
tipo = "referencia"
|
|
u_t, u_x, u_y = r_loo
|
|
lo = max(tt.min(), u_t[0]) if len(tt) > 0 else None
|
|
hi = min(tt.max(), u_t[-1]) if len(tt) > 0 else None
|
|
if len(tt) < 2 or lo is None or lo >= hi:
|
|
e_open = 0.0
|
|
e_hc = 0.0
|
|
else:
|
|
tx_t = u_t[(u_t >= lo) & (u_t <= hi)]
|
|
tx_x = np.interp(tx_t, u_t, u_x)
|
|
tx_y = np.interp(tx_t, u_t, u_y)
|
|
g_num = np.nan_to_num(d["grid_rho"], nan=0.0)
|
|
g_hc_num = np.nan_to_num(g_hc, nan=0.0)
|
|
if len(tx_t) < 2:
|
|
e_open = 0.0
|
|
e_hc = 0.0
|
|
else:
|
|
e_open = float(trapz(bilinear(g_num, tx_x, tx_y, xbins, ybins), tx_t))
|
|
e_hc = float(trapz(bilinear(g_hc_num, tx_x, tx_y, xbins, ybins), tx_t))
|
|
de_open = (e_open - e_real) / e_real * 100.0 if e_real > 0 else np.nan
|
|
de_hc = (e_hc - e_real) / e_real * 100.0 if e_real > 0 else np.nan
|
|
rows.append((d["name"], d["CI"], d["r"], d["rho_s"], e_real,
|
|
de_open, de_hc, "claro" if is_clear else "", tipo))
|
|
print("%-22s CI=%.2f r=%.3f rho_s=%.3f dE_open=%+.1f%% dE_hc=%+.1f%% [%s] %s"
|
|
% (d["name"].replace("log_", ""), d["CI"], d["r"], d["rho_s"],
|
|
de_open, de_hc, tipo, "<-- claro" if is_clear else ""))
|
|
res_path = os.path.join(OUT_DIR, "resumen_backtest.csv")
|
|
with open(res_path, "w", encoding="utf-8") as f:
|
|
f.write("log,CI,r_isc_g,spearman,E_real,dE_open_pct,dE_hc_pct,rol,tipo\n")
|
|
for r in rows:
|
|
f.write("%s,%.3f,%.3f,%.3f,%.4f,%.1f,%.1f,%s,%s\n"
|
|
% (r[0], r[1], r[2], r[3], r[4], r[5], r[6], r[7], r[8]))
|
|
|
|
# --- paso 5: figura de 3 paneles -----------------------------------------
|
|
fig = plt.figure(figsize=(18.5, 11.0))
|
|
gs = fig.add_gridspec(2, 2, width_ratios=[1.7, 1.0],
|
|
left=0.05, right=0.985, top=0.845, bottom=0.09,
|
|
wspace=0.22, hspace=0.30)
|
|
axA = fig.add_subplot(gs[:, 0])
|
|
axB = fig.add_subplot(gs[0, 1])
|
|
axC = fig.add_subplot(gs[1, 1])
|
|
|
|
im = axA.imshow(
|
|
np.ma.masked_invalid(pool_rho),
|
|
extent=[xbins[0], xbins[-1], ybins[0], ybins[-1]],
|
|
origin="lower", aspect="auto", cmap="inferno",
|
|
norm=Normalize(vmin=0, vmax=np.nanmax(pool_rho)),
|
|
interpolation="nearest")
|
|
cb = fig.colorbar(im, ax=axA, pad=0.02)
|
|
cb.set_label("Isc media por celda, normalizada por dia (pool de dias claros)",
|
|
fontsize=9)
|
|
cb.ax.tick_params(labelsize=8)
|
|
cont_g = axA.contour((xbins[:-1] + xbins[1:]) / 2, (ybins[:-1] + ybins[1:]) / 2,
|
|
pool_g, levels=5, colors="white", linewidths=0.7, alpha=0.6)
|
|
axA.clabel(cont_g, inline=True, fontsize=6.5, fmt="%.0f")
|
|
|
|
cmap = plt.get_cmap("tab10")
|
|
for k, d in enumerate(days):
|
|
axA.plot(d["x"], d["y"], ".", ms=1.4, color=cmap(k % 10), alpha=0.22,
|
|
zorder=2, label=d["name"].replace("log_", "10-"))
|
|
|
|
axA.plot(tab_x, tab_y, "-", color="white", lw=2.6, zorder=5,
|
|
label="Tabla operativa T: h -> (X*, Y*)")
|
|
step = max(1, len(tab_x) // 12)
|
|
for k in range(0, len(tab_x) - 1, step):
|
|
axA.add_patch(FancyArrowPatch((tab_x[k], tab_y[k]), (tab_x[k + 1], tab_y[k + 1]),
|
|
arrowstyle="-|>", mutation_scale=14,
|
|
color="white", lw=1.2, alpha=0.9, zorder=5))
|
|
for h in range(8, 18):
|
|
if h < tab_t[0] or h > tab_t[-1]:
|
|
continue
|
|
k = int(np.argmin(np.abs(tab_t - h)))
|
|
axA.plot(tab_x[k], tab_y[k], "o", ms=7, mfc="gold", mec="black",
|
|
mew=0.8, zorder=7)
|
|
axA.annotate(str(h), (tab_x[k], tab_y[k]), textcoords="offset points",
|
|
xytext=(5, 5), fontsize=7.5, color="gold", zorder=7)
|
|
if len(pcx):
|
|
axA.plot(pcx, pcy, "--", color="yellow", lw=1.8, zorder=4,
|
|
label="Ruta de maximo score (Dijkstra)")
|
|
axA.set_title("A) Trayectoria optima sobre el mapa de Isc (pool de dias claros)",
|
|
fontsize=11)
|
|
axA.set_xlabel("X (micro-pasos)")
|
|
axA.set_ylabel("Y (micro-pasos)")
|
|
axA.legend(loc="upper center", bbox_to_anchor=(0.5, -0.12), ncol=2,
|
|
fontsize=7.5, framealpha=0.95, borderpad=0.4, handletextpad=0.5,
|
|
labelspacing=0.4, columnspacing=1.0)
|
|
tx = axA.secondary_xaxis("top", functions=(lambda s: s * DEG_PER_USTEP,
|
|
lambda d: d / DEG_PER_USTEP))
|
|
tx.set_xlabel("X [deg]", fontsize=8, labelpad=2)
|
|
tx.tick_params(labelsize=7)
|
|
ty = axA.secondary_yaxis("right", functions=(lambda s: s * DEG_PER_USTEP,
|
|
lambda d: d / DEG_PER_USTEP))
|
|
ty.set_ylabel("Y [deg]", fontsize=8, labelpad=2, color="dimgray")
|
|
ty.tick_params(labelsize=7, colors="dimgray")
|
|
|
|
names = [r[0].replace("log_", "") for r in rows]
|
|
xx = np.arange(len(names))
|
|
w = 0.38
|
|
axB.bar(xx - w / 2, [r[5] for r in rows], w, color="steelblue",
|
|
label="tabla, bucle abierto")
|
|
axB.bar(xx + w / 2, [r[6] for r in rows], w, color="darkorange",
|
|
label="tabla + hill-climbing 1 paso")
|
|
axB.axhline(0, color="black", lw=0.9)
|
|
for k, r in enumerate(rows):
|
|
axB.annotate("CI=%.2f r=%.2f" % (r[1], r[2]),
|
|
(k, max(2.0, max(r[5], r[6])) + 2.0),
|
|
ha="center", fontsize=7)
|
|
axB.set_xticks(xx)
|
|
axB.set_xticklabels(names, rotation=30, ha="right", fontsize=8)
|
|
axB.set_ylabel("dE vs registrado [%]")
|
|
axB.set_title("B) Backtest por dia (positivo = la tabla captura mas)\n"
|
|
"dias claros evaluados con leave-one-out",
|
|
fontsize=11)
|
|
axB.grid(axis="y", ls="--", alpha=0.35)
|
|
axB.legend(fontsize=8, loc="upper left")
|
|
|
|
axC.plot(csv_t, csv_x, "-", color="steelblue", lw=2.0, label="X*(h) [upaso]")
|
|
axC.set_xlabel("hora solar")
|
|
axC.set_ylabel("X* [upaso]", color="steelblue")
|
|
axC.tick_params(axis="y", labelcolor="steelblue")
|
|
axC2 = axC.twinx()
|
|
axC2.plot(csv_t, csv_y, "--", color="darkorange", lw=2.0, label="Y*(h) [upaso]")
|
|
axC2.set_ylabel("Y* [upaso]", color="darkorange")
|
|
axC2.tick_params(axis="y", labelcolor="darkorange")
|
|
vxp = float(np.abs(np.gradient(csv_x)).max()) * 60.0
|
|
vyp = float(np.abs(np.gradient(csv_y)).max()) * 60.0
|
|
axC.set_title("C) Tabla operativa (para firmware, 10 min)\n"
|
|
"velocidad maxima: X %.0f Y %.0f upaso/min" % (vxp, vyp),
|
|
fontsize=11)
|
|
axC.grid(True, ls="--", alpha=0.35)
|
|
h1, l1 = axC.get_legend_handles_labels()
|
|
h2, l2 = axC2.get_legend_handles_labels()
|
|
axC.legend(h1 + h2, l1 + l2, fontsize=8, loc="upper left")
|
|
|
|
fig.suptitle(
|
|
"Trayectoria optima de la celda por correlacion G_POA ~ Isc - "
|
|
"%d logs, %d dias claros (CI >= 0.85*CI_max)\n"
|
|
"logica hibrida por hora del dia: setpoint de la tabla T + hill-climbing "
|
|
"sobre Isc solo con cielo claro; sin busqueda local con nubes"
|
|
% (len(days), len(clear)),
|
|
fontsize=13, y=0.99)
|
|
out = os.path.join(OUT_DIR, "trayectoria_optima_resumen.png")
|
|
fig.savefig(out, dpi=150, bbox_inches="tight")
|
|
plt.close(fig)
|
|
print("\nTabla -> %s" % os.path.basename(csv_path))
|
|
print("Figura -> %s" % os.path.basename(out))
|
|
print("Backtest -> %s" % os.path.basename(res_path))
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|