# -*- coding: utf-8 -*- """ gen_dgcup2020a.py —— 2020 电工杯 A 题《高铁牵引供电系统运行数据分析及等值建模》 ================================================================================= 确定性合成练习数据生成器（纯标准库 + 固定随机种子 SEED=2020）。 说明：赛题官方原始数据（20kHz 采样的三相电压电流录波、全天功率曲线）无法公开下载， 本脚本按赛题给定的采样率、窗口位置与物理机理构造**可复现的合成练习数据**， 供范文的四问贯通演算使用。所有数字由本脚本一次算出，范文正文、图表、附录共用， 保证"图 / 文 / 附录 / 源码"四方一致。合成数据不等于官方原始数据。 赛题约定（与题面一致）： 采样频率 fs = 20 kHz，工频 f0 = 50 Hz； 空载窗口 数据点 [1, 4000] → 0.000 ~ 0.200 s 牵引窗口 数据点 [1600000, 1604000] → 80.000 ~ 80.200 s 制动窗口 数据点 [4000000, 4004000] → 200.000 ~ 200.200 s 全程录波长度约 200 s；全天功率数据 86400 点（1 s 采样）。 运行： python3 gen_dgcup2020a.py 输出： ../assets/problems/data/dgcup2020a.csv """ import os import csv import math import random SEED = 2020 HERE = os.path.dirname(os.path.abspath(__file__)) OUT_DIR = os.path.normpath(os.path.join(HERE, "..", "assets", "problems", "data")) # ----------------------------- 全局物理常量 ----------------------------- FS = 20000.0 # 采样频率 Hz F0 = 50.0 # 工频 Hz NWIN = 4000 # 每个工况窗口点数（0.2 s = 10 个工频周期） U_HV = 220.0 # 电源侧标称线电压 kV U_TR = 27.5 # 牵引侧标称电压 kV PRICE = 0.63 # 电价 元/kWh LINE_KM = 24.6 # 单侧供电臂长度 km REC_SEC = 201 # 全程录波时长 s（含 t=200 s 制动窗口） COND = ["空载", "牵引", "制动"] WIN = {"空载": (1, 4000), "牵引": (1600000, 1604000), "制动": (4000000, 4004000)} TWIN = {"空载": 0, "牵引": 80, "制动": 200} # 窗口对应时刻 s # 各工况谐波含有率标称值（%，相对基波），以奇次为主 HARM = { "空载": {3: 12.6, 5: 9.30, 7: 6.40, 9: 4.10, 11: 3.00, 13: 2.20, 15: 1.50, 17: 1.10, 19: 0.80, 21: 0.60, 23: 0.50, 25: 0.40}, "牵引": {3: 3.90, 5: 2.80, 7: 2.05, 9: 1.20, 11: 1.55, 13: 1.10, 15: 0.52, 17: 0.44, 19: 0.33, 21: 0.26, 23: 0.62, 25: 0.41}, "制动": {3: 5.60, 5: 4.30, 7: 3.10, 9: 1.80, 11: 2.40, 13: 1.70, 15: 0.72, 17: 0.60, 19: 0.48, 21: 0.36, 23: 0.90, 25: 0.58}, } KEEP = (1, 3, 5, 7, 11, 13) # 等值模型保留的谐波源次数 # 牵引侧单相电流有效值 A 与功率因数（负值表示能量回馈） I_TR = {"空载": 5.40, "牵引": 302.0, "制动": 108.0} PF_TR = {"空载": 0.086, "牵引": 0.962, "制动": -0.941} # 电源侧电流不平衡度目标值 %（V/v 接线固有不平衡） EPS_I = {"空载": 41.5, "牵引": 28.4, "制动": 32.6} I1_ANG = {"空载": -84.6, "牵引": -31.5, "制动": 146.2} I2_ANG = {"空载": 137.2, "牵引": 142.4, "制动": -58.7} # 电源侧三相电压有效值 kV（相电压）与相角（度） U_HV_MAG = {"空载": (127.05, 127.02, 127.06), "牵引": (125.90, 127.42, 127.66), "制动": (127.61, 127.12, 126.98)} U_HV_ANG = {"空载": (0.0, -120.02, 119.98), "牵引": (0.0, -119.10, 120.72), "制动": (0.0, -120.28, 119.82)} A120 = complex(math.cos(2 * math.pi / 3), math.sin(2 * math.pi / 3)) # ----------------------------- 基础工具 ----------------------------- def _pol(mag, ang_deg): a = math.radians(ang_deg) return complex(mag * math.cos(a), mag * math.sin(a)) def _wave(n, rms, ang_deg, harms, noise=0.0022): """按基波有效值 + 谐波含有率合成一段时域波形（长度 n）。""" a1 = rms * math.sqrt(2.0) p1 = math.radians(ang_deg) hs = [] for h in sorted(harms): r = harms[h] * (1.0 + random.gauss(0.0, 0.05)) / 100.0 hs.append((h, a1 * r, random.uniform(0.0, 2.0 * math.pi))) out = [] w = 2.0 * math.pi * F0 / FS for k in range(n): v = a1 * math.sin(w * k + p1) for h, a, p in hs: v += a * math.sin(w * h * k + p) v += a1 * random.gauss(0.0, noise) out.append(v) return out def _spectrum(x, hmax=25): """整周期 DFT 求 1..hmax 次谐波幅值相量（复数，幅值制）。""" n = len(x) sp = {} for h in range(1, hmax + 1): c = 0.0 s = 0.0 w = 2.0 * math.pi * h * F0 / FS for k in range(n): a = w * k c += x[k] * math.cos(a) s += x[k] * math.sin(a) sp[h] = complex(2.0 * s / n, 2.0 * c / n) return sp def _thd(sp): f = abs(sp[1]) if f <= 0: return 0.0 q = sum(abs(sp[h]) ** 2 for h in sp if h >= 2) return math.sqrt(q) / f * 100.0 def _rms(x): return math.sqrt(sum(v * v for v in x) / len(x)) def _seq(pa, pb, pc): """对称分量法：返回 (正序, 负序, 零序) 相量。""" p = (pa + A120 * pb + A120 ** 2 * pc) / 3.0 n = (pa + A120 ** 2 * pb + A120 * pc) / 3.0 z = (pa + pb + pc) / 3.0 return p, n, z def _inv_seq(i1, i2): """由正序、负序相量反求三相相量（零序为 0）。""" ia = i1 + i2 ib = A120 ** 2 * i1 + A120 * i2 ic = A120 * i1 + A120 ** 2 * i2 return ia, ib, ic def _solve(a, b): """高斯消元解 a x = b。""" n = len(b) m = [row[:] + [b[i]] for i, row in enumerate(a)] for i in range(n): piv = max(range(i, n), key=lambda r: abs(m[r][i])) m[i], m[piv] = m[piv], m[i] d = m[i][i] if abs(d) < 1e-14: continue for j in range(i, n + 1): m[i][j] /= d for r in range(n): if r == i: continue f = m[r][i] if f: for j in range(i, n + 1): m[r][j] -= f * m[i][j] return [m[i][n] for i in range(n)] def _lsq(X, y): """最小二乘：返回 (系数, R2, RMSE, MAPE, 拟合值)。""" m = len(X) k = len(X[0]) ata = [[sum(X[i][r] * X[i][c] for i in range(m)) for c in range(k)] for r in range(k)] atb = [sum(X[i][r] * y[i] for i in range(m)) for r in range(k)] beta = _solve(ata, atb) yh = [sum(beta[j] * X[i][j] for j in range(k)) for i in range(m)] ybar = sum(y) / m sst = sum((v - ybar) ** 2 for v in y) sse = sum((y[i] - yh[i]) ** 2 for i in range(m)) r2 = 1.0 - sse / sst if sst > 0 else 0.0 rmse = math.sqrt(sse / m) mape = sum(abs((y[i] - yh[i]) / y[i]) for i in range(m)) / m * 100.0 return beta, r2, rmse, mape, yh # ----------------------------- 问题一：电能质量 ----------------------------- def q1_power_quality(): """三工况录波合成 + 序分量 / 不平衡度 / 谐波 / 功率。""" res = {} for cd in COND: s_mva = U_TR * I_TR[cd] / 1000.0 # 视在容量 MVA i1t = s_mva * 1000.0 / (math.sqrt(3.0) * U_HV) # 电源侧正序电流 A i2t = i1t * EPS_I[cd] / 100.0 pa, pb, pc = _inv_seq(_pol(i1t, I1_ANG[cd]), _pol(i2t, I2_ANG[cd])) ph_i = [] ph_u = [] thd_i = [] for idx, pz in enumerate((pa, pb, pc)): xi = _wave(NWIN, abs(pz), math.degrees(math.atan2(pz.imag, pz.real)), HARM[cd]) spi = _spectrum(xi) ph_i.append(spi[1] / math.sqrt(2.0)) thd_i.append(_thd(spi)) xu = _wave(NWIN, U_HV_MAG[cd][idx], U_HV_ANG[cd][idx], {3: 0.42, 5: 0.61, 7: 0.38, 11: 0.22, 13: 0.15}, 0.0004) ph_u.append(_spectrum(xu, 13)[1] / math.sqrt(2.0)) i1, i2, i0 = _seq(*ph_i) u1, u2, u0 = _seq(*ph_u) # 牵引侧单相馈线电流谐波 xt = _wave(NWIN, I_TR[cd], 0.0, HARM[cd]) spt = _spectrum(xt) base = abs(spt[1]) hr = {h: abs(spt[h]) / base * 100.0 for h in range(2, 26)} top5 = sorted(hr.items(), key=lambda t: -t[1])[:5] p = s_mva * PF_TR[cd] q = s_mva * math.sqrt(max(0.0, 1.0 - PF_TR[cd] ** 2)) res[cd] = { "I_ph": [abs(v) for v in ph_i], "I1": abs(i1), "I2": abs(i2), "I0": abs(i0), "eps_I": abs(i2) / abs(i1) * 100.0, "U1": abs(u1), "U2": abs(u2), "U0": abs(u0), "eps_U": abs(u2) / abs(u1) * 100.0, "THD_I_ph": thd_i, "THD_I": sum(thd_i) / 3.0, "THD_TR": _thd(spt), "hr": hr, "top5": top5, "S": s_mva, "P": p, "Q": q, "pf": PF_TR[cd], "I_TR": I_TR[cd], "win": WIN[cd], "t": TWIN[cd], "drop": None, } return res # ----------------------------- 全程 200 s 录波 ----------------------------- def q1_record(pq): """全程 201 s 逐秒功率 / 电流 / 母线电压序列（牵引侧），与三窗口严格对齐。""" p_idle = pq["空载"]["P"] p_tr = pq["牵引"]["P"] p_bk = pq["制动"]["P"] k80 = math.sin(math.pi * (80 - 20) / 85.0) ** 0.5 # 归一化因子，保证 t=80s 精确对齐 ps = [] for t in range(REC_SEC): if t < 20: p = p_idle + random.gauss(0.0, 0.004) elif t < 105: p = p_tr / k80 * math.sin(math.pi * (t - 20) / 85.0) ** 0.5 + random.gauss(0.0, 0.05) elif t < 150: p = 0.92 + 0.28 * math.sin(2 * math.pi * (t - 105) / 22.0) + random.gauss(0.0, 0.03) else: p = p_bk * math.sin(math.pi * (t - 150) / 100.0) ** 0.8 + random.gauss(0.0, 0.03) ps.append(p) ps[80] = p_tr ps[200] = p_bk ii = [abs(v) * 1000.0 / (U_TR * 0.96) for v in ps] ii[80] = pq["牵引"]["I_TR"] ii[200] = pq["制动"]["I_TR"] uu = [U_TR * 1000.0 - 1.32 * v + random.gauss(0.0, 20.0) for v in ii] pos = [v for v in ps if v > 0] neg = [v for v in ps if v < 0] return { "P": ps, "I": ii, "U": [v / 1000.0 for v in uu], "U_V": uu, "Pmax": max(ps), "Pmin": min(ps), "E_tr": sum(pos) / 3.6 * 1000.0 / 1000.0, "E_bk": -sum(neg) / 3.6 * 1000.0 / 1000.0, "Imax": max(ii), "Umin": min(uu) / 1000.0, "t_pmax": ps.index(max(ps)), } # ----------------------------- 问题三：全天动态牵引负荷 ----------------------------- TYPES = [("CRH380B", 1.15), ("CRH3C", 1.00), ("CRH2A", 0.86)] P_BASE = 4.20 # 8 编组 CRH3C 下行基准峰值 MW K_G16 = 1.92 # 16 编组倍数 D_UP = 0.85 # 上行（上坡）附加功率 MW def q3_timetable(): """构造全天列车时刻表（合成）。""" trs = [] t = 6 * 3600 + 240 while t < 23 * 3600 + 1800: hh = t / 3600.0 head = 300 if (7.0 <= hh < 10.0 or 17.0 <= hh < 20.0) else (450 if hh < 22.0 else 660) g16 = 1 if random.random() < 0.42 else 0 up = 1 if random.random() < 0.5 else 0 tn, tc = TYPES[random.randrange(3)] pk = P_BASE * (K_G16 if g16 else 1.0) * tc + (D_UP if up else 0.0) + random.gauss(0.0, 0.11) dur = 176 + random.randrange(-13, 14) bdur = 96 + random.randrange(-8, 9) bk = (0.30 if up else 0.46) + random.gauss(0.0, 0.030) trs.append({"t": int(t), "g16": g16, "up": up, "type": tn, "tc": tc, "pk": pk, "dur": dur, "bdur": bdur, "bk": bk}) t += head + random.randrange(-30, 31) return trs def q3_daily(trs): """叠加成 86400 点全天功率曲线（MW）。""" cur = [0.013] * 86400 for tr in trs: t0 = tr["t"] for k in range(tr["dur"]): if t0 + k < 86400: cur[t0 + k] += tr["pk"] * math.sin(math.pi * (k + 0.5) / tr["dur"]) ** 0.6 b0 = t0 + tr["dur"] + 12 for k in range(tr["bdur"]): if b0 + k < 86400: cur[b0 + k] -= tr["bk"] * tr["pk"] * math.sin(math.pi * (k + 0.5) / tr["bdur"]) ** 0.8 return cur def q3_extract(cur, thr=1.50): """从全天曲线提取牵引脉冲（峰值 / 持续 / 电量 / 制动深度）。""" segs = [] i = 0 n = len(cur) while i < n: if cur[i] >= thr: j = i while j < n and cur[j] >= thr: j += 1 if j - i >= 40: seg = cur[i:j] tail = cur[j:min(n, j + 150)] segs.append({"t": i, "dur": j - i, "pk": max(seg), "E": sum(seg) / 3.6, "bpk": -min(tail) if tail else 0.0}) i = j else: i += 1 for s in segs: s["ratio"] = s["bpk"] / s["pk"] if s["pk"] > 0 else 0.0 return segs def q3_model(trs, segs): """脉冲—时刻表配对 → 分类辨识 + 多元回归。""" pairs = [] for sg in segs: best = None for tr in trs: d = abs(tr["t"] - sg["t"]) if best is None or d < best[0]: best = (d, tr) if best and best[0] <= 60: pairs.append((sg, best[1])) # 一级：编组辨识（峰值阈值） thr_g = 6.30 okg = sum(1 for sg, tr in pairs if (1 if sg["pk"] >= thr_g else 0) == tr["g16"]) acc_g = okg / len(pairs) * 100.0 # 二级：方向辨识（制动/牵引峰值比阈值） thr_d = 0.380 okd = sum(1 for sg, tr in pairs if (1 if sg["ratio"] < thr_d else 0) == tr["up"]) acc_d = okd / len(pairs) * 100.0 # 三级：六类（编组×方向）最近质心 keys = sorted(set((tr["g16"], tr["up"]) for _, tr in pairs)) cents = {} for k in keys: ss = [(sg["pk"], sg["ratio"]) for (sg, tr) in pairs if (tr["g16"], tr["up"]) == k] cents[k] = (sum(a for a, _ in ss) / len(ss), sum(b for _, b in ss) / len(ss)) ok4 = 0 for (sg, tr) in pairs: bk = min(cents, key=lambda k: ((sg["pk"] - cents[k][0]) / 2.4) ** 2 + ((sg["ratio"] - cents[k][1]) / 0.09) ** 2) if bk == (tr["g16"], tr["up"]): ok4 += 1 acc4 = ok4 / len(pairs) * 100.0 # 回归：以 CRH3C 为基准中心化车型系数，含编组×车型交互项 X = [[1.0, float(tr["g16"]), float(tr["up"]), tr["tc"] - 1.0, tr["g16"] * (tr["tc"] - 1.0)] for _, tr in pairs] y = [sg["pk"] for sg, _ in pairs] beta, r2, rmse, mape, yh = _lsq(X, y) return {"pairs": len(pairs), "beta": beta, "R2": r2, "RMSE": rmse, "MAPE": mape, "acc_g": acc_g, "acc_d": acc_d, "acc4": acc4, "thr_g": thr_g, "thr_d": thr_d, "cents": cents, "y": y, "yh": yh} def q3_predict(trs, beta): """按时刻表 + 辨识模型预测次日曲线，与"实测"对比。""" act = [] pre = [] for tr in trs: a = tr["pk"] + random.gauss(0.0, 0.13) p = (beta[0] + beta[1] * tr["g16"] + beta[2] * tr["up"] + beta[3] * (tr["tc"] - 1.0) + beta[4] * tr["g16"] * (tr["tc"] - 1.0)) act.append(a) pre.append(p) n = len(act) rmse = math.sqrt(sum((act[i] - pre[i]) ** 2 for i in range(n)) / n) mape = sum(abs((act[i] - pre[i]) / act[i]) for i in range(n)) / n * 100.0 hh_a = [0.0] * 24 hh_p = [0.0] * 24 for i, tr in enumerate(trs): h = min(23, tr["t"] // 3600) hh_a[h] += act[i] * tr["dur"] * 0.78 / 3.6 hh_p[h] += pre[i] * tr["dur"] * 0.78 / 3.6 er = [abs(hh_a[h] - hh_p[h]) / hh_a[h] * 100.0 for h in range(24) if hh_a[h] > 1.0] return {"RMSE": rmse, "MAPE": mape, "act": act, "pre": pre, "hh_a": hh_a, "hh_p": hh_p, "h_mape": sum(er) / len(er), "h_max": max(er)} # ----------------------------- 问题二：节能降耗 ----------------------------- def q2_energy(cur): """全天电量统计 + 再生制动能量回收两方案比选。""" e_tr = sum(v for v in cur if v > 0) / 3.6 e_bk = -sum(v for v in cur if v < 0) / 3.6 eta0 = 0.386 e_lost = e_bk * (1.0 - eta0) plans = { "A": {"name": "车载超级电容储能", "n_train": 24, "cap": 22.5, "unit": 26.5, "eta": 0.912, "share": 0.742}, "B": {"name": "变电所地面储能", "pcs": 4.0, "cap": 640.0, "unit": None, "eta": 0.884, "share": 0.836, "inv": 620.0}, } out = {"E_tr": e_tr, "E_bk": e_bk, "eta0": eta0, "E_lost": e_lost, "bk_ratio": e_bk / e_tr * 100.0, "bill_y": e_tr * 365.0 * PRICE / 10000.0} for tag in ("A", "B"): pl = plans[tag] inv = pl.get("inv") or pl["n_train"] * pl["unit"] rec = e_bk * pl["share"] * pl["eta"] save_d = rec * PRICE / 10000.0 save_y = save_d * 365.0 out[tag] = dict(pl, inv=inv, rec=rec, save_d=save_d, save_y=save_y, payback=inv / save_y, cut=rec / e_tr * 100.0) out["best"] = "B" if out["B"]["payback"] < out["A"]["payback"] else "A" return out # ----------------------------- 问题四：等值建模 ----------------------------- def q4_equivalent(rec, pq): """戴维南等值 + RLC 等效 + 谐波阻抗 + 模型校验。""" idx = [k for k, v in enumerate(rec["I"]) if v > 12.0] ii = [rec["I"][k] for k in idx] uu = [rec["U_V"][k] for k in idx] X = [[1.0, -v] for v in ii] beta, r2, rmse, mape, uh = _lsq(X, uu) e0 = beta[0] / 1000.0 zeq = beta[1] # RLC 等效：接触网分布参数 + 无功补偿电容器组 r_km, l_km, c_km = 0.1832, 1.3540, 0.01105 R = r_km * LINE_KM L = l_km * LINE_KM / 1000.0 C_line = c_km * LINE_KM * 1e-6 C_comp = 1.85e-6 C = C_line + C_comp w = 2.0 * math.pi * F0 zh = [] for h in range(1, 41): xl = h * w * L xc = 1.0 / (h * w * C) zh.append((h, math.sqrt(R * R + (xl - xc) ** 2))) h_res = min(zh, key=lambda t: t[1])[0] f_res = 1.0 / (2.0 * math.pi * math.sqrt(L * C)) # 模型校验：保留 KEEP 次谐波源 + 基波压降，误差由频谱与压降解析导出 chk = {} for cd in COND: hr = pq[cd]["hr"] drop = zeq * pq[cd]["I_TR"] / (U_TR * 1000.0) * 100.0 e_h = math.sqrt(sum((hr[h] / 100.0) ** 2 for h in hr if h not in KEEP)) e_f = math.sqrt(1.0 + sum((hr[h] / 100.0) ** 2 for h in hr)) eh = e_h / e_f * 100.0 err = math.sqrt(eh ** 2 + drop ** 2) meas = pq[cd]["I_TR"] chk[cd] = {"meas": meas, "model": meas * (1.0 - err / 100.0), "err_h": eh, "err_f": drop, "err": err} er = [chk[c]["err"] for c in COND] return {"E0": e0, "Zeq": zeq, "R2": r2, "RMSE_U": rmse, "n_fit": len(ii), "R": R, "L": L * 1000.0, "C": C * 1e6, "C_line": C_line * 1e6, "C_comp": C_comp * 1e6, "Zh": zh, "h_res": h_res, "f_res": f_res, "h_res_exact": f_res / F0, "chk": chk, "err_mean": sum(er) / 3.0, "err_max": max(er), "r_km": r_km, "l_km": l_km, "c_km": c_km, "Z50": math.sqrt(R * R + (w * L - 1.0 / (w * C)) ** 2), "ii": ii, "uu": uu, "uh": uh} # ----------------------------- 主生成函数 ----------------------------- def gen_dgcup2020a(seed=SEED, out_csv="dgcup2020a.csv"): """确定性生成全部权威数字，写出 CSV 并返回字典 D。""" random.seed(seed) pq = q1_power_quality() rec = q1_record(pq) trs = q3_timetable() cur = q3_daily(trs) segs = q3_extract(cur) mdl = q3_model(trs, segs) prd = q3_predict(trs, mdl["beta"]) en = q2_energy(cur) eq = q4_equivalent(rec, pq) hour_mean = [sum(cur[h * 3600:(h + 1) * 3600]) / 3600.0 for h in range(24)] hour_pk = [max(cur[h * 3600:(h + 1) * 3600]) for h in range(24)] n16 = sum(1 for t in trs if t["g16"]) nup = sum(1 for t in trs if t["up"]) tcount = {} for t in trs: tcount[t["type"]] = tcount.get(t["type"], 0) + 1 D = { "seed": seed, "fs": FS, "f0": F0, "nwin": NWIN, "price": PRICE, "U_HV": U_HV, "U_TR": U_TR, "line_km": LINE_KM, "rec_sec": REC_SEC - 1, "win": WIN, "twin": TWIN, "pq": pq, "rec": rec, "keep": KEEP, "trs": trs, "trs_n": len(trs), "n16": n16, "n8": len(trs) - n16, "nup": nup, "ndown": len(trs) - nup, "tcount": tcount, "segs": segs, "segs_n": len(segs), "extract_rate": len(segs) / len(trs) * 100.0, "mdl": mdl, "prd": prd, "en": en, "eq": eq, "hour_mean": hour_mean, "hour_pk": hour_pk, "P_day_max": max(cur), "P_day_min": min(cur), "P_day_mean": sum(cur) / 86400.0, "load_rate": (sum(cur) / 86400.0) / max(cur) * 100.0, "P_base": P_BASE, "K_G16": K_G16, "D_UP": D_UP, } if not os.path.isdir(OUT_DIR): os.makedirs(OUT_DIR) rows = [("参数", "数值", "单位", "说明")] rows.append(("采样频率fs", "%.0f" % FS, "Hz", "赛题给定")) rows.append(("窗口长度", "%d" % NWIN, "点", "0.2s=10个工频周期")) rows.append(("随机种子", "%d" % seed, "-", "合成练习数据可复现")) for cd in COND: q = pq[cd] rows.append(("%s-正序电流I1" % cd, "%.3f" % q["I1"], "A", "220kV侧对称分量")) rows.append(("%s-负序电流I2" % cd, "%.3f" % q["I2"], "A", "220kV侧对称分量")) rows.append(("%s-零序电流I0" % cd, "%.4f" % q["I0"], "A", "V/v接线无零序通路")) rows.append(("%s-电流不平衡度" % cd, "%.2f" % q["eps_I"], "%", "I2/I1")) rows.append(("%s-电压不平衡度" % cd, "%.3f" % q["eps_U"], "%", "U2/U1，限值2%")) rows.append(("%s-馈线电流THD" % cd, "%.2f" % q["THD_TR"], "%", "27.5kV馈线")) rows.append(("%s-有功功率P" % cd, "%.3f" % q["P"], "MW", "牵引侧")) rows.append(("%s-无功功率Q" % cd, "%.3f" % q["Q"], "Mvar", "牵引侧")) rows.append(("%s-功率因数" % cd, "%.3f" % q["pf"], "-", "负值表示回馈")) rows.append(("全程最大牵引功率", "%.3f" % rec["Pmax"], "MW", "%ds录波" % (REC_SEC - 1))) rows.append(("全程最大制动功率", "%.3f" % rec["Pmin"], "MW", "负值为回馈")) rows.append(("全程最大电流", "%.1f" % rec["Imax"], "A", "牵引侧馈线")) rows.append(("全程最低母线电压", "%.3f" % rec["Umin"], "kV", "标称27.5kV")) rows.append(("全天牵引电量", "%.1f" % en["E_tr"], "kWh", "86400点积分")) rows.append(("全天再生制动电量", "%.1f" % en["E_bk"], "kWh", "86400点积分")) rows.append(("再生占比", "%.2f" % en["bk_ratio"], "%", "制动/牵引")) rows.append(("现状邻车吸收率", "%.1f" % (en["eta0"] * 100.0), "%", "无储能时")) for tag in ("A", "B"): a = en[tag] rows.append(("方案%s-储能容量" % tag, "%.1f" % a["cap"], "kWh", a["name"])) rows.append(("方案%s-综合效率" % tag, "%.3f" % a["eta"], "-", "充放电往返")) rows.append(("方案%s-日回收电量" % tag, "%.1f" % a["rec"], "kWh", "可利用份额%.3f" % a["share"])) rows.append(("方案%s-投资" % tag, "%.1f" % a["inv"], "万元", a["name"])) rows.append(("方案%s-年节约电费" % tag, "%.1f" % a["save_y"], "万元/年", "电价0.63元/kWh")) rows.append(("方案%s-静态回收期" % tag, "%.2f" % a["payback"], "年", "投资/年节约")) rows.append(("方案%s-节电率" % tag, "%.2f" % a["cut"], "%", "回收/牵引电量")) rows.append(("推荐方案", en["best"], "-", "回收期最短")) rows.append(("全天列车数", "%d" % len(trs), "列", "合成时刻表")) rows.append(("16编组列车数", "%d" % n16, "列", "占比%.1f%%" % (n16 / len(trs) * 100.0))) rows.append(("上行列车数", "%d" % nup, "列", "上坡方向")) rows.append(("脉冲提取数", "%d" % len(segs), "个", "阈值1.5MW")) rows.append(("提取成功率", "%.2f" % D["extract_rate"], "%", "脉冲/列车")) rows.append(("编组辨识准确率", "%.2f" % mdl["acc_g"], "%", "峰值阈值%.2fMW" % mdl["thr_g"])) rows.append(("方向辨识准确率", "%.2f" % mdl["acc_d"], "%", "制动比阈值%.3f" % mdl["thr_d"])) rows.append(("四象限六类准确率", "%.2f" % mdl["acc4"], "%", "编组×方向最近质心")) rows.append(("回归R2", "%.4f" % mdl["R2"], "-", "峰值功率模型")) rows.append(("回归RMSE", "%.4f" % mdl["RMSE"], "MW", "峰值功率模型")) rows.append(("回归MAPE", "%.2f" % mdl["MAPE"], "%", "峰值功率模型")) rows.append(("次日预测RMSE", "%.4f" % prd["RMSE"], "MW", "逐车峰值")) rows.append(("次日预测MAPE", "%.2f" % prd["MAPE"], "%", "逐车峰值")) rows.append(("逐小时电量MAPE", "%.2f" % prd["h_mape"], "%", "24点")) rows.append(("全天最大功率", "%.3f" % D["P_day_max"], "MW", "86400点")) rows.append(("全天平均功率", "%.3f" % D["P_day_mean"], "MW", "86400点")) rows.append(("负荷率", "%.2f" % D["load_rate"], "%", "平均/最大")) rows.append(("戴维南电势E0", "%.4f" % eq["E0"], "kV", "最小二乘辨识")) rows.append(("等值阻抗Zeq", "%.4f" % eq["Zeq"], "Ω", "最小二乘辨识")) rows.append(("辨识R2", "%.4f" % eq["R2"], "-", "U-I线性拟合")) rows.append(("压降拟合RMSE", "%.2f" % eq["RMSE_U"], "V", "%d个样本" % eq["n_fit"])) rows.append(("等效电阻R", "%.4f" % eq["R"], "Ω", "%.1fkm供电臂" % LINE_KM)) rows.append(("等效电感L", "%.3f" % eq["L"], "mH", "%.1fkm供电臂" % LINE_KM)) rows.append(("等效电容C", "%.4f" % eq["C"], "μF", "接触网+补偿电容器组")) rows.append(("谐振频率", "%.1f" % eq["f_res"], "Hz", "串联谐振")) rows.append(("谐振次数", "%.2f" % eq["h_res_exact"], "次", "接近11/13次谐波")) for cd in COND: c = eq["chk"][cd] rows.append(("%s-等值模型误差" % cd, "%.2f" % c["err"], "%", "谐波截断%.2f%%+压降%.2f%%" % (c["err_h"], c["err_f"]))) rows.append(("等值模型平均误差", "%.2f" % eq["err_mean"], "%", "三工况")) rows.append(("等值模型最大误差", "%.2f" % eq["err_max"], "%", "三工况")) path = os.path.join(OUT_DIR, out_csv) with open(path, "w", encoding="utf-8-sig", newline="") as f: csv.writer(f).writerows(rows) D["csv"] = path D["csv_rows"] = len(rows) - 1 return D if __name__ == "__main__": D = gen_dgcup2020a() print("=" * 66) print("2020 电工杯 A 题 高铁牵引供电系统运行数据分析及等值建模") print("确定性合成练习数据（种子 %d，非官方原始数据）" % D["seed"]) print("=" * 66) print("[问题一] 三工况电能质量（fs=%.0fHz，窗口 %d 点）" % (D["fs"], D["nwin"])) for cd in COND: q = D["pq"][cd] print(" %s 窗口[%d,%d]@%ds I1=%.3fA I2=%.3fA I0=%.4fA εI=%.2f%% εU=%.3f%%" % (cd, q["win"][0], q["win"][1], q["t"], q["I1"], q["I2"], q["I0"], q["eps_I"], q["eps_U"])) print(" 馈线 I=%.1fA THD=%.2f%% P=%.3fMW Q=%.3fMvar pf=%.3f 前3谐波 %s" % (q["I_TR"], q["THD_TR"], q["P"], q["Q"], q["pf"], ", ".join("%d次%.2f%%" % (h, v) for h, v in q["top5"][:3]))) R = D["rec"] print("[全程录波] %ds 最大牵引 %.3fMW(t=%ds) 最大制动 %.3fMW 最大电流 %.1fA 最低母线 %.3fkV" % (D["rec_sec"], R["Pmax"], R["t_pmax"], R["Pmin"], R["Imax"], R["Umin"])) print(" 窗口内电量：牵引 %.2f kWh，再生 %.2f kWh" % (R["E_tr"], R["E_bk"])) E = D["en"] print("[问题二] 全天牵引电量 %.1f kWh，再生制动电量 %.1f kWh（占比 %.2f%%）" % (E["E_tr"], E["E_bk"], E["bk_ratio"])) print(" 现状邻车吸收率 %.1f%%，弃用 %.1f kWh/日；年电费 %.1f 万元" % (E["eta0"] * 100.0, E["E_lost"], E["bill_y"])) for tag in ("A", "B"): a = E[tag] print(" 方案%s %s：容量 %.1fkWh，效率 %.3f，回收 %.1f kWh/日 → " "%.1f 万元/年，投资 %.1f 万元，回收期 %.2f 年，节电率 %.2f%%" % (tag, a["name"], a["cap"], a["eta"], a["rec"], a["save_y"], a["inv"], a["payback"], a["cut"])) print(" 推荐方案：%s" % E["best"]) M = D["mdl"] P = D["prd"] print("[问题三] 全天列车 %d 列（16编组 %d / 8编组 %d；上行 %d / 下行 %d）" % (D["trs_n"], D["n16"], D["n8"], D["nup"], D["ndown"])) print(" 车型分布 %s" % D["tcount"]) print(" 脉冲提取 %d 个（成功率 %.2f%%），配对 %d 个" % (D["segs_n"], D["extract_rate"], M["pairs"])) print(" 编组辨识 %.2f%%，方向辨识 %.2f%%，六类最近质心 %.2f%%" % (M["acc_g"], M["acc_d"], M["acc4"])) print(" 峰值回归 P=%.4f%+.4f·G16%+.4f·UP%+.4f·(Kt-1)%+.4f·G16(Kt-1)" % (M["beta"][0], M["beta"][1], M["beta"][2], M["beta"][3], M["beta"][4])) print(" R2=%.4f RMSE=%.4fMW MAPE=%.2f%%；次日 RMSE=%.4fMW MAPE=%.2f%% 逐小时电量MAPE=%.2f%%" % (M["R2"], M["RMSE"], M["MAPE"], P["RMSE"], P["MAPE"], P["h_mape"])) print(" 全天最大 %.3fMW 平均 %.3fMW 负荷率 %.2f%%" % (D["P_day_max"], D["P_day_mean"], D["load_rate"])) Q = D["eq"] print("[问题四] 戴维南等值 E0=%.4fkV Zeq=%.4fΩ（R2=%.4f，压降RMSE=%.2fV，%d样本）" % (Q["E0"], Q["Zeq"], Q["R2"], Q["RMSE_U"], Q["n_fit"])) print(" RLC 等效（%.1fkm）R=%.4fΩ L=%.3fmH C=%.4fμF（网%.4f+补偿%.2f）" % (D["line_km"], Q["R"], Q["L"], Q["C"], Q["C_line"], Q["C_comp"])) print(" 串联谐振 %.1fHz = %.2f 次；工频阻抗模 %.3fΩ" % (Q["f_res"], Q["h_res_exact"], Q["Z50"])) for cd in COND: c = Q["chk"][cd] print(" %s 实测 %.2fA / 模型 %.2fA，误差 %.2f%%（谐波截断 %.2f%% + 压降 %.2f%%）" % (cd, c["meas"], c["model"], c["err"], c["err_h"], c["err_f"])) print(" 三工况平均误差 %.2f%%，最大误差 %.2f%%" % (Q["err_mean"], Q["err_max"])) print("[输出] %s（%d 行数据）" % (D["csv"], D["csv_rows"]))