import json, sys, os import numpy as np from scipy.stats import chi2 sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from senl import FeederPV, NLEstimator # noqa FEEDERS = os.path.abspath(os.environ.get("G4W_FEEDERS", "feeders")) CASES = { "ieee13": (FEEDERS + "/13Bus/IEEE13Nodeckt.dss", "650632", ("650",), "SOURCEBUS"), "ieee34": (FEEDERS + "/34Bus/ieee34Mod1.dss", "L1", ("800",), "SOURCEBUS"), "ieee123": (FEEDERS + "/123Bus/IEEE123Master.dss", "L115", ("149",), "150"), } # parâmetros dos dados de carga (cenário sintético inspirado nos problemas descritos para ADMS) P = dict(ami_cov=0.6, gran=((0.03, 0.4), (0.06, 0.3), (0.15, 0.3)), # (erro relativo, fração) para 5, 15, 60 min miss=0.10, gross=0.03, typical_sig=0.5, naive_sig=0.05, pv_irr=(0.8, 0.1)) def profile(fd, rng): """Atributos fixos de cada carga: tem medição, granularidade.""" n = len(fd.loads) has = rng.random(n) < P["ami_cov"] g = rng.choice(len(P["gran"]), n, p=[f for _, f in P["gran"]]) return has, np.array([P["gran"][k][0] for k in g]) def pseudo(fd, est, kw, kvar, pvkw_irr, has, gerr, rng, mode): """Gera P,Q (kW) por carga e sigma (kVA) segundo o tratamento de dados ('ingenuo' ou 'tratado').""" n = len(kw); P_net = kw - pvkw_irr; Q = kvar.copy() nomP, nomQ = np.array([L["kw"] for L in fd.loads]), np.array([L["kvar"] for L in fd.loads]) pv_exp = fd.pv_kw * P["pv_irr"][0] estP = np.empty(n); estQ = np.empty(n); sig = np.empty(n) miss = rng.random(n) < P["miss"]; gross = rng.random(n) < P["gross"] for i in range(n): typP = nomP[i] - (pv_exp[i] if mode == "tratado" else 0.0) # curva típica (tratado desconta GD cadastrada) if not has[i]: estP[i], estQ[i], sig[i] = typP, nomQ[i], P["typical_sig"] * abs(nomP[i] + 1j * nomQ[i]); continue e = 1 + gerr[i] * rng.standard_normal() mP, mQ = P_net[i] * e, Q[i] * e prosumer = fd.pv_kw[i] > 0 if mode == "ingenuo": if miss[i]: mP, mQ = 0.0, 0.0 # intervalo faltante vira zero if prosumer: mP = abs(mP) # sem sentido de fluxo: só o módulo if gross[i]: mP, mQ = 4 * mP, 4 * mQ # kWh de 15 min lido como kW estP[i], estQ[i], sig[i] = mP, mQ, P["naive_sig"] * max(abs(mP + 1j * mQ), 0.1 * abs(nomP[i] + 1j * nomQ[i])) else: bad = miss[i] or (gross[i] and abs(4 * P_net[i]) > 3 * abs(nomP[i]) + 1e-9) if prosumer and not bad: # sentido desconhecido: usa sinal pelo cadastro de GD s = np.sign(typP) if typP != 0 else 1.0 mP = s * abs(mP); sg = 0.3 * abs(nomP[i] + 1j * nomQ[i]) elif gross[i] and not bad: mP, mQ = 4 * mP, 4 * mQ; sg = gerr[i] * abs(mP + 1j * mQ) * 4 # erro grosseiro não detectado else: sg = gerr[i] * abs(mP + 1j * mQ) + 0.02 * abs(nomP[i] + 1j * nomQ[i]) if bad: mP, mQ, sg = typP, nomQ[i], P["typical_sig"] * abs(nomP[i] + 1j * nomQ[i]) estP[i], estQ[i], sig[i] = mP, mQ, sg S = -(est.M @ (estP + 1j * estQ)) # injeção (sai da rede = consumo) por nó, VA sigN = np.sqrt((np.abs(est.M) ** 2) @ sig ** 2) + 50.0 # M já converte kVA -> VA return S, sigN def run(name, N=200, seed=2026, out=".", order=None): p, sl, sv, root = CASES[name] fd = FeederPV(p, sub_meter_line=sl, sub_v_buses=sv, source_buses=(root,)) rng = np.random.default_rng(seed) has, gerr = profile(fd, rng) nL = len(fd.loads) kwN = np.array([L["kw"] for L in fd.loads]); kvN = np.array([L["kvar"] for L in fd.loads]) sensors_all = order or [] variants = [(m, K) for K in (0, 2, 4) for m in ("ingenuo", "auto", "tratado")] ests = {K: NLEstimator(fd, sensors_all[:K]) for K in (0, 2, 4)} skipped = 0 res = {v[0] + f"_K{v[1]}": {"status": [], "err_v": [], "chi_rej": [], "it": [], "err_p": [], "removed": []} for v in variants} for s in range(N): k = np.clip(1 + 0.3 * rng.standard_normal(nL), 0.2, 2.0) * np.clip(1 + 0.15 * rng.standard_normal(), 0.6, 1.4) irr = float(np.clip(rng.normal(*P["pv_irr"]), 0.3, 1.0)) Vt = fd.truth(k, irr) if Vt is None: skipped += 1; continue for mode, K in variants: est = ests[K]; r = res[f"{mode}_K{K}"] S, sg = pseudo(fd, est, kwN * k, kvN * k, fd.pv_kw * irr, has, gerr, np.random.default_rng(1000 + s), "ingenuo" if mode == "auto" else mode) z, sig = est.build_z(Vt, S, sg, np.random.default_rng(5000 + s)) if mode == "auto": nz = len(est.zero_nodes); nl = len(est.load_nodes) x, it, st, J, rem, sig = est.solve_clean(z, sig, list(range(2 * nz, 2 * nz + 2 * nl))) r["removed"].append(rem) else: x, it, st, J = est.solve(z, sig) r["status"].append(st); r["it"].append(it) if st == "ok": n = est.n; V = (x[:n] + 1j * x[n:]) * est.vb nodes = [i for i in range(n) if i not in fd.source] r["err_v"].append(float(np.max(np.abs(np.abs(V[nodes]) - np.abs(Vt[nodes])) / est.vb[nodes]) * 100)) dof = len(z) - 2 * n r["chi_rej"].append(bool(J > chi2.ppf(0.99, max(dof, 1)))) L = est.load_nodes St = Vt[L] * np.conj(fd.Y[L] @ Vt); Se = V[L] * np.conj(fd.Y[L] @ V) r["err_p"].append(float(np.sum(np.abs(Se.real - St.real)) / np.sum(np.abs(St.real)) * 100)) fd.reset() summ = {} for key, r in res.items(): st = np.array(r["status"]); ok = st == "ok" summ[key] = {"conv": float(ok.mean() * 100), "div": float((st != "ok").mean() * 100), "rej": float(np.mean(r["chi_rej"]) * 100) if r["chi_rej"] else None, "reliable": float(np.sum(ok & np.isin(np.arange(len(st)), np.where(ok)[0][~np.array(r["chi_rej"], bool)] if r["chi_rej"] else [])) / len(st) * 100), "err_v_med": float(np.median(r["err_v"])) if r["err_v"] else None, "err_v_p90": float(np.percentile(r["err_v"], 90)) if r["err_v"] else None, "err_p_med": float(np.median(r["err_p"])) if r["err_p"] else None, "it_med": float(np.median(r["it"])), "removed_med": float(np.mean(r["removed"])) if r["removed"] else 0.0, "status_counts": {u: int((st == u).sum()) for u in set(st)}} print(name, key, summ[key]) json.dump({"feeder": name, "params": P, "sensors": sensors_all[:4], "summary": summ, "raw": res, "loads": nL, "skipped_truth": skipped, "pv_loads": int((fd.pv_kw > 0).sum()), "ami": int(has.sum())}, open(os.path.join(out, f"p2_{name}.json"), "w"), indent=1, default=float) if __name__ == "__main__": out = sys.argv[1]; N = int(os.environ.get("N", "200")) p1 = os.environ.get("P1_RESULTS", "") for nm in sys.argv[2:]: order = json.load(open(f"{p1}/results_{nm}.json"))["modes"]["VI"]["greedy_order"] if p1 else None run(nm, N=N, out=out, order=order)