# -*- coding: utf-8 -*-
"""
CUMCM 2023 A —— 定日镜场的优化设计（The Heliostat Field Optimization）

确定性真源模型（纯标准库，零随机噪声；SEED 仅作版本锚定）。
建模对象：塔式太阳能光热电站的定日镜场。
  Q1 单台定日镜几何/成本优化：在反射面成本 + 支撑结构成本（随面积超线性增长）
      与光学年发电之间权衡，求使单位发电成本(LCOE)最低的单镜尺寸与长宽比。
  Q2 镜场布置优化：以塔为中心做径向交错(radial-staggered)布局，逐台计算余弦效率、
      大气衰减、阴影/遮挡，求总年发电与镜场 LCOE；并对比三种布局策略。
  Q3 经济性：CAPEX(镜场+塔/吸热器+土地) + OPEX，折算 LCOE，做参数敏感性。

输出：
  - assets/problems/data/cumcm2023a.csv  （参数 + 单镜优化 + 各环指标 + 布局对比）
  - 返回 D（dict），供范文附录独立复现全部权威数字。

用法：cd tools/ && python3 gen_cumcm2023a.py
"""
import os
import csv
import math

SEED = 2023  # 版本锚定，模型本身确定性、无随机

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

# ---------------- 物理 / 经济常数 ----------------
H_TOWER = 120.0       # 定日镜塔高 (m)
RHO = 0.92            # 镜面反射率
DNI = 2100.0          # 年法向直射辐照 (kWh/m^2/yr, 典型沙漠)
ETA_TE = 0.40         # 热->电转换效率 (power block)
ATT_DELTA = 1500.0    # 大气衰减尺度 (m)，eta_att = exp(-d/ATT_DELTA)
SPILL = 0.97          # 溢出/截断损失系数
HOURS = 8760.0        # 年小时数（仅用于说明，DNI 已为年量）

# 单镜成本模型: C(S) = C_M*S + C_S*S^1.3 + C_B
C_M = 150.0           # 镜面成本 ($/m^2)
C_S = 8.0             # 支撑结构成本系数（随面积超线性）
C_B = 2000.0          # 基座/传动基础成本 ($)
CRF = 0.08            # 资本回收因子 (annualization)
OPEX_RATE = 0.02      # 年运维占 CAPEX 比例

# 镜场布局
N_HELIO = 600
INNER_R = 60.0        # 最内圈半径 (m)
RING_DR = 18.0        # 环间距 (m)
SPACING = 14.0        # 同环相邻镜中心距 (m)
FIELD_MAX_R = 600.0
TOWER_COST = 30.0e6   # 塔+吸热器 CAPEX ($)
LAND_UNIT = 3.0       # 土地单价 ($/m^2)


# ---------------- Q1: 单台定日镜几何/成本优化 ----------------
def single_opt():
    """求使单镜 LCOE 最低的面积 S* 与方镜边长。

    解析推导：年发电 E(S)=DNI*S*eta0*ETA_TE (eta0 为参考位置光学效率，与 S 无关)，
    年成本 C_a(S)=(CRF+OPEX)*(C_M*S + C_S*S^1.3 + C_B)。
    LCOE(S)=C_a/(E0*S)，其中 E0=DNI*eta0*ETA_TE。
    令 dLCOE/dS=0 => (C_M + 1.3*C_S*S^0.3)*S = C_M*S + C_S*S^1.3 + C_B
                 => 0.3*C_S*S^1.3 = C_B
                 => S* = (C_B/(0.3*C_S))**(1/1.3)
    长宽比固定为 1（正方形在给定面积下周长最小、支撑最省，且光学效率与形状无关）。
    """
    S_star = (C_B / (0.3 * C_S)) ** (1.0 / 1.3)
    w = h = math.sqrt(S_star)
    # 参考位置（单镜最优布设点，取 r=0 即塔基正下方远处投影的代表值用 r=120 近似）
    r0 = 120.0
    d0 = math.sqrt(r0 ** 2 + H_TOWER ** 2)
    eta_cos0 = 0.90 * (H_TOWER / d0) ** 0.5
    eta_att0 = math.exp(-d0 / ATT_DELTA)
    eta_sb0 = 1.0 - 0.10 * (INNER_R / r0)
    eta0 = RHO * eta_cos0 * eta_att0 * eta_sb0 * SPILL
    E_single = DNI * S_star * eta0 * ETA_TE            # kWh/yr
    C_single = C_M * S_star + C_S * (S_star ** 1.3) + C_B
    C_a = (CRF + OPEX_RATE) * C_single                 # 年化成本
    LCOE_single = C_a / E_single * 1000.0              # $/MWh
    return dict(S_star=S_star, w=w, h=h, aspect=w / h, eta0=eta0,
                E_single=E_single, C_single=C_single, C_a=C_a,
                LCOE_single=LCOE_single)


# ---------------- Q2: 镜场布置 ----------------
def build_field(inner_r=INNER_R, ring_dr=RING_DR, n=N_HELIO, spacing=SPACING,
                max_r=FIELD_MAX_R):
    """径向交错布局：逐环半径 r_k=inner_r+k*ring_dr，每环 M_k=floor(2π r_k/spacing) 台。
    返回 helios 列表，每项为 dict(r, d, x, y, eta_cos, eta_att, eta_sb, eta, E)。"""
    helios = []
    k = 0
    while len(helios) < n and inner_r + k * ring_dr <= max_r:
        r = inner_r + k * ring_dr
        M = max(1, int(2 * math.pi * r / spacing))
        # 交错：奇数环整体错开半间距
        phase = (math.pi / M) if (k % 2 == 1) else 0.0
        for j in range(M):
            if len(helios) >= n:
                break
            ang = phase + 2 * math.pi * j / M
            x = r * math.cos(ang)
            y = r * math.sin(ang)
            d = math.sqrt(r ** 2 + H_TOWER ** 2)
            eta_cos = 0.90 * (H_TOWER / d) ** 0.5
            eta_att = math.exp(-d / ATT_DELTA)
            eta_sb = 1.0 - 0.10 * (INNER_R / r)
            eta = RHO * eta_cos * eta_att * eta_sb * SPILL
            E = DNI * 0.0  # placeholder, filled per-area below
            helios.append(dict(r=r, d=d, x=x, y=y, eta_cos=eta_cos,
                               eta_att=eta_att, eta_sb=eta_sb, eta=eta, E=E))
        k += 1
    return helios


def field_metrics(helios, S):
    """给定单镜面积 S，计算各镜年发电与镜场汇总。"""
    total_E = 0.0
    cos_sum = att_sum = sb_sum = eta_sum = 0.0
    rings = {}
    for hh in helios:
        E = DNI * S * hh["eta"] * ETA_TE
        hh["E"] = E
        total_E += E
        cos_sum += hh["eta_cos"]; att_sum += hh["eta_att"]
        sb_sum += hh["eta_sb"]; eta_sum += hh["eta"]
        rk = round(hh["r"])
        rings.setdefault(rk, dict(n=0, E=0.0, cos=0.0, att=0.0, sb=0.0))
        rings[rk]["n"] += 1
        rings[rk]["E"] += E
        rings[rk]["cos"] += hh["eta_cos"]
        rings[rk]["att"] += hh["eta_att"]
        rings[rk]["sb"] += hh["eta_sb"]
    n = len(helios)
    rs = sorted(rings.keys())
    ring_rows = []
    R_max = max(hh["r"] for hh in helios)
    for rk in rs:
        g = rings[rk]
        ring_rows.append(dict(r=rk, n=g["n"], E=g["E"] / 1000.0,
                               cos=g["cos"] / g["n"], att=g["att"] / g["n"],
                               sb=g["sb"] / g["n"]))
    return dict(n=n, total_E=total_E / 1000.0,  # MWh/yr
                avg_cos=cos_sum / n, avg_att=att_sum / n,
                avg_sb=sb_sum / n, avg_eta=eta_sum / n,
                R_max=R_max, ring_rows=ring_rows,
                land_area=math.pi * R_max ** 2)


def layout_comparison(S):
    """三种布局策略对比：返回 (name, helios, metrics, CAPEX, LCOE)。"""
    configs = [
        ("A-径向交错(默认)", dict(inner_r=60, ring_dr=18, spacing=14)),
        ("B-同心均布", dict(inner_r=60, ring_dr=18, spacing=18)),
        ("C-内密外疏", dict(inner_r=45, ring_dr=22, spacing=12)),
    ]
    out = []
    for name, kw in configs:
        hs = build_field(**kw)
        m = field_metrics(hs, S)
        land = m["land_area"] * LAND_UNIT
        capex = hs.__len__() * (C_M * S + C_S * (S ** 1.3) + C_B) + TOWER_COST + land
        annual = capex * CRF + capex * OPEX_RATE
        lcoe = annual / m["total_E"]  # $/MWh (total_E 已为 MWh)
        out.append(dict(name=name, n=m["n"], total_E=m["total_E"],
                        R_max=m["R_max"], land_cost=land, capex=capex,
                        lcoe=lcoe))
    return out


# ---------------- 主入口 ----------------
def gen_cumcm2023a(seed=SEED, out_csv="cumcm2023a.csv"):
    so = single_opt()
    S = so["S_star"]
    hs_default = build_field()
    m = field_metrics(hs_default, S)
    comp = layout_comparison(S)
    best = min(comp, key=lambda c: c["lcoe"])

    # 经济性（以默认布局计）
    land = m["land_area"] * LAND_UNIT
    capex = m["n"] * (C_M * S + C_S * (S ** 1.3) + C_B) + TOWER_COST + land
    annual = capex * CRF + capex * OPEX_RATE
    lcoe_field = annual / m["total_E"]
    # 敏感性：DNI、CRF、ETA_TE 各 ±10%
    sens = []
    for label, fac in [("DNI+10%", 1.10), ("DNI-10%", 0.90),
                       ("CRF+10%", 1.10), ("CRF-10%", 0.90),
                       ("ETA_TE+10%", 1.10), ("ETA_TE-10%", 0.90)]:
        if "DNI" in label:
            E2 = m["total_E"] * fac
        elif "ETA" in label:
            E2 = m["total_E"] * fac
        else:
            capex2 = capex * fac
            E2 = m["total_E"]
            lcoe2 = (capex2 * CRF + capex2 * OPEX_RATE) / E2
            sens.append((label, round(lcoe2 - lcoe_field, 2)))
            continue
        lcoe2 = annual / E2
        sens.append((label, round(lcoe2 - lcoe_field, 2)))

    # 写 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(["单镜最优面积S*", round(S, 2), "m2", "使单镜LCOE最低"])
        w.writerow(["单镜边长(方镜)", round(so["w"], 2), "m", "长宽比=1(正方形)"])
        w.writerow(["单镜反射效率eta0", round(so["eta0"], 4), "-", "RHO*cos*att*sb*spill(参考位)"])
        w.writerow(["单镜年发电", round(so["E_single"] / 1000.0, 2), "MWh/yr", "DNI*S*eta0*ETA_TE"])
        w.writerow(["单镜CAPEX", round(so["C_single"], 0), "$", "C_M*S+C_S*S^1.3+C_B"])
        w.writerow(["单镜LCOE", round(so["LCOE_single"], 2), "$/MWh", "年化成本/年发电"])
        w.writerow(["镜场定日镜数N", m["n"], "台", "径向交错布局"])
        w.writerow(["镜场最外圈半径", round(m["R_max"], 1), "m", "土地半径"])
        w.writerow(["镜场年发电", round(m["total_E"], 1), "MWh/yr", "Σ各镜"])
        w.writerow(["平均光学效率", round(m["avg_eta"], 4), "-", "RHO*cos*att*sb*spill"])
        w.writerow(["镜场CAPEX", round(capex / 1e6, 2), "M$", "镜场+塔/吸热器+土地"])
        w.writerow(["镜场LCOE", round(lcoe_field, 2), "$/MWh", "年化/年发电"])
        w.writerow(["最优布局", best["name"], "-", "LCOE最低策略"])
        w.writerow([])
        w.writerow(["环号(半径m)", "镜数", "单环年发电(MWh)", "平均余弦效率", "平均大气衰减", "平均遮挡效率"])
        for rr in m["ring_rows"]:
            w.writerow([rr["r"], rr["n"], round(rr["E"], 2),
                        round(rr["cos"], 4), round(rr["att"], 4), round(rr["sb"], 4)])
        w.writerow([])
        w.writerow(["布局策略", "镜数", "年发电(MWh)", "最外半径(m)", "土地成本(M$)", "CAPEX(M$)", "LCOE($/MWh)"])
        for c in comp:
            w.writerow([c["name"], c["n"], round(c["total_E"], 1),
                        round(c["R_max"], 1), round(c["land_cost"] / 1e6, 3),
                        round(c["capex"] / 1e6, 2), round(c["lcoe"], 2)])

    return dict(
        single=so, S=S, field=m, layouts=comp, best=best,
        capex=capex, lcoe_field=lcoe_field, sensitivity=sens,
        csv_path=path,
        summary="CUMCM2023A 定日镜场优化：单镜方镜S*=%g m2、镜场N=%d、年发电%.1f MWh、LCOE=%.2f $/MWh"
                 % (S, m["n"], m["total_E"], lcoe_field),
    )


if __name__ == "__main__":
    D = gen_cumcm2023a()
    print(D["summary"])
    print("单镜 S* = %.2f m2 (%.2f x %.2f, 长宽比=%.3f)"
          % (D["S"], D["single"]["w"], D["single"]["h"], D["single"]["aspect"]))
    print("单镜 LCOE = %.2f $/MWh" % D["single"]["LCOE_single"])
    print("镜场 N = %d, 年发电 = %.1f MWh, LCOE = %.2f $/MWh"
          % (D["field"]["n"], D["field"]["total_E"], D["lcoe_field"]))
    print("平均光学效率 = %.4f (cos %.4f / att %.4f / sb %.4f)"
          % (D["field"]["avg_eta"], D["field"]["avg_cos"],
             D["field"]["avg_att"], D["field"]["avg_sb"]))
    print("布局对比 (按LCOE):")
    for c in sorted(D["layouts"], key=lambda x: x["lcoe"]):
        print("  %s: N=%d E=%.1f MWh LCOE=%.2f" % (c["name"], c["n"], c["total_E"], c["lcoe"]))
    print("敏感性(对LCOE的$/MWh变化):", D["sensitivity"])