# -*- coding: utf-8 -*-
"""
电工杯 2023 A —— 电采暖负荷参与电力系统功率调节的技术经济分析
（Electric Heating Loads Participating in Power-System Power Regulation: Techno-Economic Analysis）

确定性真源模型（纯标准库，固定 SEED=2023；仅初始温度分布的随机抽样使用 seed，其余为解析/确定性）。
建模对象：居民区电采暖（直热式/空气源热泵等效为额定 8kW 户用加热器）聚合参与电网功率调节。

真实赛题（已核验）包含六个部分：
  Q1 典型住户用电行为：单体两节点（室内—墙体）一阶热惯性模型稳态解性态 + 24h 室温/开关曲线。
  Q2 单户调节能力：给定室外温度，计算功率上/下调可持续时间，分析温度影响。
  Q3 多户（6 户）调节能力：总用电功率曲线 + 各时点可上/下调总功率，分析温度影响。
  Q4 住宅区（600 户）调节能力：24h 开关/总功率曲线 + 各时点可上/下调总功率曲线。
  Q5 削峰填谷收益：削峰时段最大持续下调功率、填谷时段最大持续上调功率、参与后温度可行性、
                  供暖期 180 天分温度档的收益/户均收益/节省供热成本百分比。
  Q6 展望（定性，省级 4000 万 m² 与南方空调）：本脚本仅给出定量量级，正文展开。

郑重声明：官方原赛题数据集（附件 A/B 的逐时负荷、补偿价表等）未公开可下载，本文件所生成数据为
「确定性合成数据集」，仅用于跑通建模方法与训练，NOT 官方原题数据。所有数字可由本脚本独立复现。

运行：cd tools/ && python3 gen_dgcup2023a.py
"""
import os
import csv
import math
import random

SEED = 2023  # 版本锚定 + 初始温度分布抽样的随机种子

HERE = os.path.dirname(os.path.abspath(__file__))
OUT_DIR = os.path.normpath(os.path.join(HERE, "..", "assets", "problems", "data"))

# ---------------- 物理 / 经济基本常数（显式声明，确定性） ----------------
P_HEAT = 8.0          # 户用制热额定功率 (kW)
T_LOW = 18.0          # 室内温控下界 (℃)
T_HIGH = 22.0         # 室内温控上界 (℃)
DEAD_ON = 19.0        # 恒温死区：低于此开
DEAD_OFF = 21.0       # 恒温死区：高于此关
ETA = 0.30            # 制热功率注入墙体的比例（其余进室内空气）
C_IN = 2.0            # 室内空气等效热容 (kWh/℃)
C_W = 6.0             # 墙体等效热容 (kWh/℃)
K1 = 0.20            # 室内—墙体 传热系数 (kW/℃)
K2 = 0.15            # 室内—室外 传热系数 (kW/℃)，含围护综合传热
K3 = 0.08            # 墙体—室外 传热系数 (kW/℃)
DT = 0.05            # 积分步长 (h)，24h 共 480 步

N_HOUSE = 600         # 住宅区住户数（赛题给定）
N_SAMPLE = 6          # Q3 抽样户数

# 经济参数（确定性合成，量级贴近典型峰谷分时与辅助服务补偿）
E_PEAK = 0.85         # 峰时段电价 (元/kWh)
E_FLAT = 0.55         # 平时段电价 (元/kWh)
E_VALLEY = 0.35       # 谷时段电价 (元/kWh)
C_DOWN = 0.60         # 削峰（向下调节）辅助服务补偿 (元/kWh)
C_UP = 0.50           # 填谷（向上调节）辅助服务补偿 (元/kWh)
PEAK_HOURS = 2.0      # 削峰时段长度 (h)：取 [19,21)
VALLEY_HOURS = 3.0    # 填谷时段长度 (h)：取 [1,4)
PEAK_START, PEAK_END = 19.0, 21.0
VALLEY_START, VALLEY_END = 1.0, 4.0

# 供暖期 180 天分温度档（合成 Table 2）：室外平均温度 ℃ -> 持续天数
TABLE2 = [(-15.0, 10), (-10.0, 30), (-5.0, 50), (0.0, 50), (5.0, 40)]


# ---------------- 两节点一阶热惯性模型 ----------------
def steady_state(Tout, u, k2=K2, cin=C_IN):
    """解析稳态：u=1 开、u=0 关。返回 (Tin, Tw)。可传 k2/cin 做敏感性。"""
    A = (1.0 - ETA) * P_HEAT * u
    B = ETA * P_HEAT * u
    # -(K1+k2) Tin + K1 Tw = -k2 Tout - A
    #  K1 Tin -(K1+K3) Tw = -K3 Tout - B
    a11, a12 = -(K1 + k2), K1
    a21, a22 = K1, -(K1 + K3)
    b1 = -k2 * Tout - A
    b2 = -K3 * Tout - B
    det = a11 * a22 - a12 * a21
    Tin = (b1 * a22 - a12 * b2) / det
    Tw = (a11 * b2 - b1 * a21) / det
    return Tin, Tw


def deriv(Tin, Tw, Tout, u, k2=K2, cin=C_IN):
    dTin = (K1 * (Tw - Tin) + k2 * (Tout - Tin) + (1.0 - ETA) * P_HEAT * u) / cin
    dTw = (K1 * (Tin - Tw) + K3 * (Tout - Tw) + ETA * P_HEAT * u) / C_W
    return dTin, dTw


def simulate(Tout, Tin0, Tw0, u0, hours=24.0):
    """基线恒温死区控制：低于 19℃ 开、高于 21℃ 关、否则保持。
    返回逐时（每步）数组 Tin, Tw, u 以及以 1h 抽样的序列。"""
    n = int(round(hours / DT))
    Tin, Tw, u = Tin0, Tw0, u0
    Tin_l, Tw_l, u_l = [], [], []
    for _ in range(n):
        dTin, dTw = deriv(Tin, Tw, Tout, u)
        Tin += dTin * DT
        Tw += dTw * DT
        # 死区控制
        if Tin < DEAD_ON:
            u = 1
        elif Tin > DEAD_OFF:
            u = 0
        Tin_l.append(Tin); Tw_l.append(Tw); u_l.append(u)
    return Tin_l, Tw_l, u_l


def sustain_time(Tin, Tw, Tout, u_forced, bound):
    """从 (Tin,Tw) 起以 u_forced 强制运行，直到 Tin 越过 bound（18 或 22）。
    返回可持续时间 (h)，上限 12h（视为可持续）。"""
    t = 0.0
    x, y = Tin, Tw
    for _ in range(int(12.0 / DT) + 1):
        dTin, dTw = deriv(x, y, Tout, u_forced)
        x += dTin * DT
        y += dTw * DT
        t += DT
        if u_forced == 0 and x <= bound:
            return t
        if u_forced == 1 and x >= bound:
            return t
    return 12.0


def capacity_curve(Tout, Tin0, Tw0, u0):
    """逐时（1h 分辨率）返回基线开关、总功率、可下调/上调功率与可持续时间。
    用于单户 / 多户 / 住宅区统一口径。"""
    Tin_l, Tw_l, u_l = simulate(Tout, Tin0, Tw0, u0)
    # 以 1h 为抽样的时刻点
    hours = list(range(24))
    res = []
    for h in hours:
        idx = int(round(h / DT)) - 1
        idx = max(0, min(len(u_l) - 1, idx))
        Tin, Tw, u = Tin_l[idx], Tw_l[idx], u_l[idx]
        # 下调（削峰）：当前开则关，直到 18℃
        if u == 1:
            down_p = P_HEAT
            down_t = sustain_time(Tin, Tw, Tout, 0, T_LOW)
        else:
            down_p = 0.0
            down_t = 0.0
        # 上调（填谷）：当前关则开，直到 22℃
        if u == 0:
            up_p = P_HEAT
            up_t = sustain_time(Tin, Tw, Tout, 1, T_HIGH)
        else:
            up_p = 0.0
            up_t = 0.0
        res.append(dict(h=h, Tin=Tin, Tw=Tw, u=u, power=u * P_HEAT,
                        down_p=down_p, up_p=up_p, down_t=down_t, up_t=up_t))
    return res


def aggregate(N, Tout, seed=SEED):
    """生成 N 户，初始温度在 [18,22] 均匀（seeded），模拟基线 24h，
    返回逐时：总功率、可下调总功率、可上调总功率、可下调户数、可上调户数。"""
    rng = random.Random(seed)
    households = []
    for _ in range(N):
        Tin0 = 18.0 + rng.random() * 4.0        # 均匀 [18,22]
        Tw0 = Tin0 - 2.0                          # 墙体略低
        u0 = 1 if Tin0 < 20.0 else 0
        households.append((Tin0, Tw0, u0))
    agg = [dict(h=h, total=0.0, down=0.0, up=0.0, n_on=0, n_off=0,
                down_t_min=12.0, up_t_min=12.0) for h in range(24)]
    per_house = []
    for (Tin0, Tw0, u0) in households:
        cap = capacity_curve(Tout, Tin0, Tw0, u0)
        per_house.append(cap)
        for h in range(24):
            c = cap[h]
            a = agg[h]
            a["total"] += c["power"]
            a["down"] += c["down_p"]
            a["up"] += c["up_p"]
            if c["u"] == 1:
                a["n_on"] += 1
            else:
                a["n_off"] += 1
            if c["u"] == 1:
                a["down_t_min"] = min(a["down_t_min"], c["down_t"])
            else:
                a["up_t_min"] = min(a["up_t_min"], c["up_t"])
    return dict(N=N, Tout=Tout, agg=agg, per_house=per_house,
                f_on=[a["n_on"] / N for a in agg])


# ---------------- Q1：单体稳态与 24h 曲线 ----------------
def q1():
    ss_on = steady_state(-15.0, 1)
    ss_off = steady_state(-15.0, 0)
    # 稳态随室外温度变化（开）
    ss_curve = []
    for Tout in [-20, -15, -10, -5, 0, 5]:
        ss_curve.append((Tout, *steady_state(Tout, 1)))
    # 参数敏感性对 ON 稳态室温（固定 Tout=-15）
    sens = []
    for name, k2b, cinb in [("K2+20%", 0.20 * 1.2, C_IN),
                            ("K2-20%", 0.20 * 0.8, C_IN),
                            ("C_IN+20%", K2, C_IN * 1.2),
                            ("C_IN-20%", K2, C_IN * 0.8)]:
        Tin, Tw = steady_state(-15.0, 1, k2=k2b, cin=cinb)
        sens.append((name, round(Tin, 3)))
    # 24h 曲线（基线，Tout=-5，典型供暖日）
    cap = capacity_curve(-5.0, 20.0, 18.0, 1)
    on_hours = sum(1 for c in cap if c["u"] == 1)
    energy = sum(c["power"] for c in cap)  # kW * 1h = kWh（逐时）
    return dict(ss_on=ss_on, ss_off=ss_off, ss_curve=ss_curve, sens=sens,
                q1b_Tout=-5.0, q1b_on_hours=on_hours, q1b_energy=energy,
                q1b_Tin=[round(c["Tin"], 2) for c in cap],
                q1b_u=[c["u"] for c in cap])


# ---------------- Q2：单户调节能力（Tout=-15, Tin0=20, u0=1） ----------------
def q2():
    Tout = -15.0
    cap = capacity_curve(Tout, 20.0, 18.0, 1)
    # 全天下调/上调可持续时间（取非零时刻平均与最大）
    down_times = [c["down_t"] for c in cap if c["u"] == 1]
    up_times = [c["up_t"] for c in cap if c["u"] == 0]
    avg_down = sum(down_times) / len(down_times) if down_times else 0.0
    max_down = max(down_times) if down_times else 0.0
    avg_up = sum(up_times) / len(up_times) if up_times else 0.0
    max_up = max(up_times) if up_times else 0.0
    # 不同室外温度对单户可持续时间的影响
    temp_scan = []
    for Tout2 in [-20, -15, -10, -5, 0, 5]:
        c2 = capacity_curve(Tout2, 20.0, 18.0, 1)
        dt2 = [x["down_t"] for x in c2 if x["u"] == 1]
        ut2 = [x["up_t"] for x in c2 if x["u"] == 0]
        temp_scan.append((Tout2,
                          round(sum(dt2) / len(dt2), 3) if dt2 else 0.0,
                          round(sum(ut2) / len(ut2), 3) if ut2 else 0.0))
    return dict(Tout=Tout, cap=cap, avg_down=avg_down, max_down=max_down,
                avg_up=avg_up, max_up=max_up, temp_scan=temp_scan)


# ---------------- Q3：6 户调节能力（Tout=-20） ----------------
def q3():
    rng = random.Random(SEED + 1)
    houses = []
    for i in range(N_SAMPLE):
        Tin0 = 18.0 + rng.random() * 4.0
        Tw0 = Tin0 - 2.0
        u0 = 1 if Tin0 < 20.0 else 0
        houses.append((Tin0, Tw0, u0))
    Tout = -20.0
    caps = [capacity_curve(Tout, Tin0, Tw0, u0) for (Tin0, Tw0, u0) in houses]
    total = [sum(c[h]["power"] for c in caps) for h in range(24)]
    down_total = [sum(c[h]["down_p"] for c in caps) for h in range(24)]
    up_total = [sum(c[h]["up_p"] for c in caps) for h in range(24)]
    # 不同室外温度
    temp_scan = []
    for Tout2 in [-20, -15, -10, -5, 0]:
        caps2 = [capacity_curve(Tout2, Tin0, Tw0, u0) for (Tin0, Tw0, u0) in houses]
        dt2 = sum(sum(c[h]["down_p"] for h in range(24)) for c in caps2) / 24.0
        ut2 = sum(sum(c[h]["up_p"] for h in range(24)) for c in caps2) / 24.0
        temp_scan.append((Tout2, round(dt2, 2), round(ut2, 2)))
    return dict(Tout=Tout, houses=houses, caps=caps, total=total,
                down_total=down_total, up_total=up_total, temp_scan=temp_scan)


# ---------------- Q4：住宅区 600 户调节能力 ----------------
def q4():
    aggs = {}
    for Tout in [-15.0, -5.0, 5.0]:
        aggs[Tout] = aggregate(N_HOUSE, Tout)
    # 默认展示 Tout=-15 的曲线
    a = aggs[-15.0]
    total_curve = [round(a["agg"][h]["total"], 1) for h in range(24)]
    down_curve = [round(a["agg"][h]["down"], 1) for h in range(24)]
    up_curve = [round(a["agg"][h]["up"], 1) for h in range(24)]
    return dict(aggs=aggs, total_curve=total_curve,
                down_curve=down_curve, up_curve=up_curve,
                agg_default_Tout=-15.0)


# ---------------- Q5：削峰填谷收益（180 天分温度档） ----------------
def q5():
    # 对每个温度档求住宅区可下调（削峰时段）/可上调（填谷时段）可持续总功率
    bucket = []
    total_heat_energy = 0.0
    for Tout, days in TABLE2:
        a = aggregate(N_HOUSE, Tout)
        # 削峰时段平均可下调功率（可持续，取削峰窗内最小）
        down_pk = [a["agg"][h]["down"] for h in range(int(PEAK_START), int(PEAK_END))]
        up_vl = [a["agg"][h]["up"] for h in range(int(VALLEY_START), int(VALLEY_END))]
        down_kW = min(down_pk) if down_pk else 0.0       # 持续最大下调
        up_kW = min(up_vl) if up_vl else 0.0             # 持续最大上调
        f_on = sum(a["agg"][h]["n_on"] for h in range(24)) / 24.0 / N_HOUSE
        daily_heat = f_on * N_HOUSE * P_HEAT * 24.0      # 社区单日采暖耗电 kWh
        bucket.append(dict(Tout=Tout, days=days, down_kW=down_kW, up_kW=up_kW,
                           f_on=f_on, daily_heat=daily_heat))
        total_heat_energy += daily_heat * days
    # 全年调节电量
    down_energy_annual = sum(b["down_kW"] * PEAK_HOURS * b["days"] for b in bucket)
    up_energy_annual = sum(b["up_kW"] * VALLEY_HOURS * b["days"] for b in bucket)
    # 收益
    rev_down = down_energy_annual * C_DOWN          # 辅助服务（削峰）补偿
    rev_up = up_energy_annual * C_UP                # 辅助服务（填谷）补偿
    rev_arb = down_energy_annual * (E_PEAK - E_VALLEY)  # 峰谷套利（削峰电量移至谷段）
    rev_service = rev_down + rev_up
    rev_total = rev_service + rev_arb
    base_cost = total_heat_energy * E_FLAT          # 基准供热成本（平价计）
    save_pct = rev_total / base_cost * 100.0
    per_house = rev_total / N_HOUSE
    return dict(bucket=bucket, down_energy_annual=down_energy_annual,
                up_energy_annual=up_energy_annual, rev_down=rev_down,
                rev_up=rev_up, rev_arb=rev_arb, rev_service=rev_service,
                rev_total=rev_total, base_cost=base_cost, save_pct=save_pct,
                per_house=per_house, total_heat_energy=total_heat_energy)


# ---------------- Q6：省级 / 南方空调展望（定量量级） ----------------
def q6():
    AREA = 4000e4        # 4000 万 m²
    LOAD_DENSITY = 8.0   # 单位采暖面积功率密度 W/m²
    total_kW = AREA * LOAD_DENSITY / 1000.0   # kW
    # 以住宅区 600 户（约等效面积）推算可调节功率比例
    a = aggregate(N_HOUSE, -5.0)
    down_kW_600 = min(a["agg"][h]["down"] for h in range(int(PEAK_START), int(PEAK_END)))
    equiv_houses = total_kW / (N_HOUSE * P_HEAT / N_HOUSE)  # = total_kW / 8 per house
    n_equiv = total_kW / P_HEAT
    province_down_kW = down_kW_600 * (n_equiv / N_HOUSE)
    return dict(area=AREA, load_density=LOAD_DENSITY, total_kW=total_kW,
                n_equiv=n_equiv, province_down_kW=province_down_kW)


# ---------------- 主入口 ----------------
def gen_dgcup2023a(seed=SEED, out_csv="dgcup2023a.csv"):
    r1 = q1()
    r2 = q2()
    r3 = q3()
    r4 = q4()
    r5 = q5()
    r6 = q6()

    # 写 CSV（参数 + 各节权威数字面板）
    os.makedirs(OUT_DIR, exist_ok=True)
    path = os.path.join(OUT_DIR, out_csv)
    with open(path, "w", newline="", encoding="utf-8-sig") as f:
        w = csv.writer(f)
        w.writerow(["模块", "参数/指标", "数值", "单位", "说明"])
        w.writerow(["常数", "户用制热功率 P", P_HEAT, "kW", "直热/热泵等效额定"])
        w.writerow(["常数", "温控区间", "%g~%g" % (T_LOW, T_HIGH), "℃", "舒适带"])
        w.writerow(["常数", "墙体注入比例 ETA", ETA, "-", "其余进室内空气"])
        w.writerow(["常数", "室内空气热容 C_IN", C_IN, "kWh/℃", "一阶热惯性"])
        w.writerow(["常数", "墙体热容 C_W", C_W, "kWh/℃", ""])
        w.writerow(["常数", "K1 室-墙", K1, "kW/℃", ""])
        w.writerow(["常数", "K2 室-外", K2, "kW/℃", "围护综合"])
        w.writerow(["常数", "K3 墙-外", K3, "kW/℃", ""])
        w.writerow([])
        w.writerow(["Q1", "开稳态室温(Tout=-15)", "%.3f" % r1["ss_on"][0], "℃", "解析"])
        w.writerow(["Q1", "开稳态墙温(Tout=-15)", round(r1["ss_on"][1], 3), "℃", "解析"])
        w.writerow(["Q1", "关稳态室温(Tout=-15)", round(r1["ss_off"][0], 3), "℃", "解析"])
        w.writerow(["Q1", "24h 开机时长(Tout=-5)", r1["q1b_on_hours"], "h", "死区控制"])
        w.writerow(["Q1", "24h 采暖耗电(Tout=-5)", round(r1["q1b_energy"], 1), "kWh", "逐时求和"])
        w.writerow([])
        w.writerow(["Q2", "单户平均下调可持续", round(r2["avg_down"], 3), "h", "Tout=-15"])
        w.writerow(["Q2", "单户最大下调可持续", "%.3f" % r2["max_down"], "h", "Tout=-15"])
        w.writerow(["Q2", "单户平均上调可持续", round(r2["avg_up"], 3), "h", "Tout=-15"])
        w.writerow(["Q2", "单户最大上调可持续", round(r2["max_up"], 3), "h", "Tout=-15"])
        w.writerow([])
        w.writerow(["Q3", "6户总功率均值(Tout=-20)", round(sum(r3["total"]) / 24.0, 1), "kW", ""])
        w.writerow(["Q3", "6户平均可下调总功率", round(sum(r3["down_total"]) / 24.0, 1), "kW", ""])
        w.writerow(["Q3", "6户平均可上调总功率", round(sum(r3["up_total"]) / 24.0, 1), "kW", ""])
        w.writerow([])
        a4 = r4["aggs"][-15.0]
        w.writerow(["Q4", "600户平均总功率(Tout=-15)", round(sum(a4["agg"][h]["total"] for h in range(24)) / 24.0, 1), "kW", ""])
        w.writerow(["Q4", "600户峰值总功率(Tout=-15)", round(max(a4["agg"][h]["total"] for h in range(24)), 1), "kW", ""])
        w.writerow(["Q4", "600户可下调功率均值(Tout=-15)", round(sum(a4["agg"][h]["down"] for h in range(24)) / 24.0, 1), "kW", ""])
        w.writerow(["Q4", "600户可上调功率均值(Tout=-15)", round(sum(a4["agg"][h]["up"] for h in range(24)) / 24.0, 1), "kW", ""])
        w.writerow([])
        w.writerow(["Q5", "全年下调电量", round(r5["down_energy_annual"], 0), "kWh", "180天·削峰2h"])
        w.writerow(["Q5", "全年上调电量", round(r5["up_energy_annual"], 0), "kWh", "180天·填谷3h"])
        w.writerow(["Q5", "辅助服务收益", round(r5["rev_service"], 0), "元", "补偿"])
        w.writerow(["Q5", "峰谷套利收益", round(r5["rev_arb"], 0), "元", ""])
        w.writerow(["Q5", "总收益", round(r5["rev_total"], 0), "元", "=辅助+套利"])
        w.writerow(["Q5", "基准供热成本", round(r5["base_cost"], 0), "元", "平价计"])
        w.writerow(["Q5", "节省供热成本%", "%.2f" % r5["save_pct"], "%", ""])
        w.writerow(["Q5", "平均每户收益", round(r5["per_house"], 0), "元", "总/600"])
        w.writerow([])
        w.writerow(["Q6", "省级采暖总功率(4000万m²)", round(r6["total_kW"], 0), "kW", ""])
        w.writerow(["Q6", "省级等效户数", round(r6["n_equiv"], 0), "户", "8kW/户"])
        w.writerow(["Q6", "省级可下调功率(估)", round(r6["province_down_kW"], 0), "kW", "线性外推"])

    return dict(q1=r1, q2=r2, q3=r3, q4=r4, q5=r5, q6=r6,
                csv_path=path,
                summary="电工杯2023A：P=8kW, 600户, 总收益%.0f元, 节省成本%.1f%%"
                        % (r5["rev_total"], r5["save_pct"]))


if __name__ == "__main__":
    D = gen_dgcup2023a()
    r1, r2, r3, r4, r5, r6 = D["q1"], D["q2"], D["q3"], D["q4"], D["q5"], D["q6"]
    print("=" * 60)
    print("电工杯 2023 A 真源模型 —— 确定性合成数据（SEED=%d）" % SEED)
    print("=" * 60)
    print("[Q1] 稳态(Tout=-15): 开(室%.3f,墙%.3f) 关(室%.3f)"
          % (r1["ss_on"][0], r1["ss_on"][1], r1["ss_off"][0]))
    print("     稳态随室外温度(开):", [(t, round(ti, 2), round(tw, 2)) for t, ti, tw in r1["ss_curve"]])
    print("     参数敏感性(开稳态室):", r1["sens"])
    print("     Q1b(Tout=-5): 开机%.0fh 耗电%.1f kWh"
          % (r1["q1b_on_hours"], r1["q1b_energy"]))
    print("[Q2] 单户(Tout=-15): 平均下调%.3fh 最大%.3fh | 平均上调%.3fh 最大%.3fh"
          % (r2["avg_down"], r2["max_down"], r2["avg_up"], r2["max_up"]))
    print("     温度扫描(室外,平均下调,平均上调):", r2["temp_scan"])
    print("[Q3] 6户(Tout=-20): 总功率均值%.1f 可下调均值%.1f 可上调均值%.1f"
          % (sum(r3["total"]) / 24.0, sum(r3["down_total"]) / 24.0, sum(r3["up_total"]) / 24.0))
    print("     温度扫描:", r3["temp_scan"])
    a4 = r4["aggs"][-15.0]
    print("[Q4] 600户(Tout=-15): 总功率均值%.1f 峰值%.1f | 可下调均值%.1f 可上调均值%.1f"
          % (sum(a4["agg"][h]["total"] for h in range(24)) / 24.0,
             max(a4["agg"][h]["total"] for h in range(24)),
             sum(a4["agg"][h]["down"] for h in range(24)) / 24.0,
             sum(a4["agg"][h]["up"] for h in range(24)) / 24.0))
    print("     可下调曲线:", r4["down_curve"])
    print("     可上调曲线:", r4["up_curve"])
    print("[Q5] 分温度档:", [(b["Tout"], round(b["down_kW"], 0), round(b["up_kW"], 0)) for b in r5["bucket"]])
    print("     全年下调电量%.0f kWh 上调电量%.0f kWh" % (r5["down_energy_annual"], r5["up_energy_annual"]))
    print("     辅助服务收益%.0f 套利%.0f 总收益%.0f 基准成本%.0f"
          % (r5["rev_service"], r5["rev_arb"], r5["rev_total"], r5["base_cost"]))
    print("     节省供热成本%.2f%% 户均收益%.0f 元" % (r5["save_pct"], r5["per_house"]))
    print("[Q6] 省级总功率%.0f kW 等效户数%.0f 可下调(估)%.0f kW"
          % (r6["total_kW"], r6["n_equiv"], r6["province_down_kW"]))
    print("CSV ->", D["csv_path"])