# -*- coding: utf-8 -*- """ ICM 2023 F — Green GDP（绿色 GDP 与全球气候减缓治理） 真实赛题：COMAP 国际交叉学科建模竞赛（ICM）F 题要求 (1) 从已有的多种 GGDP 计算方法中，选择一种"若取代 GDP 成为经济健康主要衡量标准， 可能对气候减缓产生可衡量影响"的口径； (2) 建立一个简单、易辩护的模型，估计若采用所选择的 GGDP 作为国家经济健康主要标准， 对全球气候减缓的预期影响（自行决定如何度量全球影响）； (3) 用 GGDP 替代 GDP 可能遭遇抵抗；根据模型判断是否值得在全球范围内转轨， 比较气候减缓的潜在上行与替换现状所需努力的潜在下行； (4) 选择一个国家做更深入分析，说明转轨可能对其产生什么影响（自然资源使用/节约的具体变化）， 并明确联系 GDP 与 GGDP 计算方式的变化； (5) 基于该国的分析，给其领导人写一页非技术报告，建议支持或拒绝转轨。 真源模型（纯标准库，确定性，无随机噪声；SEED 仅作版本锚定）： 1. 所选 GGDP 口径（SEEA 式环境经济核算口径）： GGDP_i = GDP_i − D_res_i − D_env_i D_res_i = res_share_i · GDP_i # 自然资源耗减成本（占 GDP 比例） D_env_i = env_share_i · GDP_i # 环境退化成本（含污染与碳排代理，占 GDP 比例） 绿色比 green_ratio_i = GGDP_i / GDP_i = 1 − res_share_i − env_share_i。 2. 排名变化：按 GDP 与按 GGDP 分别排名，rank_drop = ggd_rank − gdp_rank（>0 表示相对地位下降）。 3. 全球气候减缓机制（指标驱动减排）：当 GGDP 取代 GDP 成为"头条指标"，政府为改善头条分数 会压降 D_env（碳排放直接拉低 GGDP），触发减排。每国年减排比例： r_i(φ) = R_MAX · φ · s_i 其中 φ 为全球采纳比例（情景 0.3/0.6/1.0），s_i 为敏感系数（资源/环境/化石依赖越高， GGDP 对它的惩罚越重、减排激励越强），R_MAX 为满采纳时的最大年减排比例。 稳态年排放 E'_i = E_i·(1 − r_i)；累计避免排放（至 2050，T 年）ΔC_i = E_i · r_i · T。 4. 全球影响度量（三指标）： ΔC_global = GLOBAL_SCALE_CO2 · Σ ΔC_i # 累计避免 CO2（Gt） Δppm = ΔC_global · PPM_PER_GT # 对应大气 CO2 浓度下降（ppm） ΔT = GAMMA_T · Δppm # 对应增温放缓（°C，瞬态响应） 气候效益$ = ΔC_global · 1e9 · SCC # 避免的气候损害（美元，SCC 社会碳成本） 5. 转轨成本（下行）：一次性制度/核算接入成本 = EPS_COST · GlobalGDP · φ（全球 GDP 缩放后）。 净收益 = 气候效益$ − 转轨成本$。 6. 抵抗指数 resist_i = res_share_i + env_share_i（绿色惩罚越重、越倾向抵制转轨）。 7. 深析国：巴西（id=9）——森林/资源依赖高，GGDP 将"森林退减"计入 D_env， 改变其资源使用激励；输出其 GGDP 分解、排名变化、排放路径(BAU vs GGDP)、成本效益与建议。 8. 鲁棒性：对 ΔT(φ=1) 做参数 tornado（R_MAX / GAMMA_T / GLOBAL_SCALE_CO2 / T 各 ±20%）。 输出：assets/problems/data/icm2023f.csv（32 行国家面板 + 表头） 作者说明：国家面板为方法论演示用的确定性合成数据（种子 2023），数值贴近公开量级但非官方实测； 真实参赛应以 World Bank / SEEA / UNFCCC 等公开库校准；本模型聚焦"口径选择→ 指标驱动减排→全球影响→值得性→国别深析"的可复现链路。 """ import os import math SEED = 2023 HERE = os.path.dirname(os.path.abspath(__file__)) # --------------------------------------------------------------------------- # 0. 模型常数（均作为假设显式声明，便于敏感性分析） # --------------------------------------------------------------------------- T = 27 # 推演年限：2023 → 2050 SCC = 50.0 # 社会碳成本 Social Cost of Carbon（美元 / 吨 CO2） R_MAX = 0.22 # 满全球采纳时单国最大年减排比例 GAMMA_T = 0.011 # 瞬态气候响应（°C / ppm CO2） PPM_PER_GT = 0.127 # 1 Gt CO2 ≈ 0.127 ppm 大气 CO2 GLOBAL_SCALE_CO2 = 1.20 # 面板(≈84% 全球排放) → 全球 缩放 GLOBAL_SCALE_GDP = 1.235 # 面板 GDP → 全球 GDP 缩放 EPS_COST = 0.0015 # 一次性转轨成本率（占全球 GDP） # 敏感系数归一化区间（跨面板） ENV_MIN, ENV_MAX = 0.03, 0.13 RES_MIN, RES_MAX = 0.01, 0.14 # 全球采纳情景 PHIS = [0.3, 0.60, 1.0] PHI_CN = {0.3: "部分采纳(30%)", 0.60: "多数采纳(60%)", 1.0: "全球采纳(100%)"} # --------------------------------------------------------------------------- # 1. 国家面板（32 国；数值为确定性合成演示数据，量级贴近公开统计） # (id, 名称, 区域, GDP(万亿$), CO2(Gt/yr), res_share, env_share, dev, fossil_share) # --------------------------------------------------------------------------- COUNTRIES = [ (1, "美国 USA", "美洲", 25.0, 5.00, 0.03, 0.06, 0.95, 0.80), (2, "中国 China", "亚洲", 17.8, 11.00, 0.05, 0.12, 0.70, 0.85), (3, "日本 Japan", "亚洲", 4.2, 1.00, 0.02, 0.04, 0.92, 0.75), (4, "德国 Germany", "欧洲", 4.0, 0.65, 0.03, 0.05, 0.93, 0.70), (5, "印度 India", "亚洲", 3.4, 2.80, 0.04, 0.10, 0.45, 0.80), (6, "英国 UK", "欧洲", 3.1, 0.35, 0.02, 0.04, 0.92, 0.65), (7, "法国 France", "欧洲", 2.8, 0.30, 0.02, 0.04, 0.91, 0.50), (8, "意大利 Italy", "欧洲", 2.1, 0.32, 0.03, 0.05, 0.89, 0.65), (9, "巴西 Brazil", "美洲", 2.0, 0.48, 0.09, 0.12, 0.60, 0.55), (10, "加拿大 Canada", "美洲", 2.0, 0.55, 0.08, 0.07, 0.90, 0.70), (11, "俄罗斯 Russia", "欧洲", 1.8, 1.70, 0.10, 0.10, 0.65, 0.85), (12, "韩国 S.Korea", "亚洲", 1.7, 0.62, 0.02, 0.06, 0.89, 0.78), (13, "澳大利亚 Australia","大洋洲",1.6, 0.40, 0.12, 0.08, 0.91, 0.85), (14, "墨西哥 Mexico", "美洲", 1.4, 0.50, 0.05, 0.08, 0.58, 0.80), (15, "西班牙 Spain", "欧洲", 1.4, 0.25, 0.03, 0.05, 0.88, 0.55), (16, "印尼 Indonesia", "亚洲", 1.3, 0.62, 0.08, 0.11, 0.50, 0.75), (17, "沙特 Saudi", "中东", 1.1, 0.60, 0.14, 0.07, 0.60, 0.95), (18, "土耳其 Turkey", "欧洲", 1.0, 0.45, 0.04, 0.08, 0.62, 0.82), (19, "荷兰 Netherlands","欧洲", 0.99, 0.14, 0.02, 0.05, 0.92, 0.70), (20, "瑞士 Switzerland","欧洲", 0.80, 0.04, 0.01, 0.03, 0.95, 0.45), (21, "波兰 Poland", "欧洲", 0.69, 0.30, 0.05, 0.08, 0.83, 0.85), (22, "瑞典 Sweden", "欧洲", 0.59, 0.04, 0.03, 0.03, 0.94, 0.40), (23, "泰国 Thailand", "亚洲", 0.52, 0.27, 0.04, 0.08, 0.62, 0.75), (24, "阿联酋 UAE", "中东", 0.50, 0.22, 0.13, 0.06, 0.75, 0.96), (25, "挪威 Norway", "欧洲", 0.49, 0.04, 0.07, 0.03, 0.95, 0.40), (26, "南非 S.Africa", "非洲", 0.40, 0.44, 0.09, 0.11, 0.50, 0.90), (27, "埃及 Egypt", "非洲", 0.40, 0.25, 0.05, 0.08, 0.45, 0.80), (28, "越南 Vietnam", "亚洲", 0.43, 0.34, 0.04, 0.10, 0.48, 0.78), (29, "阿根廷 Argentina", "美洲", 0.64, 0.19, 0.06, 0.07, 0.65, 0.70), (30, "伊朗 Iran", "中东", 0.37, 0.65, 0.10, 0.09, 0.55, 0.88), (31, "尼日利亚 Nigeria","非洲", 0.36, 0.12, 0.10, 0.09, 0.35, 0.70), (32, "哈萨克 Kazakh.", "亚洲", 0.26, 0.30, 0.11, 0.09, 0.58, 0.90), ] REGION_CN = {"美洲": "美洲", "亚洲": "亚洲", "欧洲": "欧洲", "非洲": "非洲", "大洋洲": "大洋洲", "中东": "中东"} # 深析国 DEEP_ID = 9 # 巴西 # --------------------------------------------------------------------------- # 2. 核心计算 # --------------------------------------------------------------------------- def _clip(v, lo=0.0, hi=1.0): return max(lo, min(hi, v)) def _corr(a, b): m = len(a) if m < 2: return 0.0 ma = sum(a) / m mb = sum(b) / m cov = sum((a[i] - ma) * (b[i] - mb) for i in range(m)) sa = math.sqrt(sum((x - ma) ** 2 for x in a)) sb = math.sqrt(sum((x - mb) ** 2 for x in b)) return cov / (sa * sb) if sa * sb > 0 else 0.0 def _sensitivity(res, env, fossil): """敏感系数 s_i ∈ [0,1]：资源/环境/化石依赖越高，GGDP 减排激励越强。""" env_n = (env - ENV_MIN) / (ENV_MAX - ENV_MIN) res_n = (res - RES_MIN) / (RES_MAX - RES_MIN) s = 0.45 * env_n + 0.35 * res_n + 0.20 * fossil return _clip(s) def _rank(values): """返回按值降序的 rank（1=最大）；并列取平均后稳定。""" order = sorted(range(len(values)), key=lambda k: values[k], reverse=True) rank = [0] * len(values) i = 0 while i < len(order): j = i while j + 1 < len(order) and values[order[j + 1]] == values[order[i]]: j += 1 avg = (i + j) / 2.0 + 1 # 1-based 平均 rank for k in range(i, j + 1): rank[order[k]] = avg i = j + 1 return rank # --------------------------------------------------------------------------- # 3. 主生成 # --------------------------------------------------------------------------- def gen_icm2023f(seed=SEED, out_csv="icm2023f.csv"): # 3.1 逐国基础量与 GGDP rows = [] for (i, name, region, gdp, co2, res, env, dev, fossil) in COUNTRIES: d_res = res * gdp d_env = env * gdp ggd = gdp - d_res - d_env green = ggd / gdp if gdp > 0 else 0.0 ci = co2 / gdp if gdp > 0 else 0.0 # 碳强度 Gt/万亿$ s = _sensitivity(res, env, fossil) resist = res + env rows.append({ "id": i, "name": name, "region": region, "gdp": gdp, "co2": co2, "res": res, "env": env, "dev": dev, "fossil": fossil, "d_res": d_res, "d_env": d_env, "ggd": ggd, "green": green, "ci": ci, "sens": s, "resist": resist, }) # 3.2 排名 gdp_rank = _rank([r["gdp"] for r in rows]) ggd_rank = _rank([r["ggd"] for r in rows]) for k, r in enumerate(rows): r["gdp_rank"] = gdp_rank[k] r["ggd_rank"] = ggd_rank[k] r["rank_drop"] = ggd_rank[k] - gdp_rank[k] # >0 表示相对地位下降 # 3.3 减排与全球影响（各 φ 情景） panel_co2 = sum(r["co2"] for r in rows) panel_gdp = sum(r["gdp"] for r in rows) global_gdp = panel_gdp * GLOBAL_SCALE_GDP scenarios = {} for phi in PHIS: dC_panel = 0.0 for r in rows: r_phi = R_MAX * phi * r["sens"] dC_panel += r["co2"] * r_phi * T dC_global = dC_panel * GLOBAL_SCALE_CO2 dppm = dC_global * PPM_PER_GT dT = GAMMA_T * dppm benefit = dC_global * 1e9 * SCC # 美元 cost = EPS_COST * global_gdp * phi * 1e12 # 美元（global_gdp 单位万亿$） net = benefit - cost scenarios[phi] = { "phi": phi, "dC_panel": dC_panel, "dC_global": dC_global, "dppm": dppm, "dT": dT, "benefit": benefit, "cost": cost, "net": net, "bcr": benefit / cost if cost > 0 else float("inf"), } # 逐国满采纳减排（存于 φ=1 情景对象，供 CSV/图使用） if phi == 1.0: for r in rows: r_1 = R_MAX * 1.0 * r["sens"] r["r_full"] = r_1 r["avoided_co2_full"] = r["co2"] * r_1 * T # 3.4 区域聚合 regions = sorted(set(r["region"] for r in rows)) by_region = {} for rg in regions: sub = [r for r in rows if r["region"] == rg] by_region[rg] = { "n": len(sub), "gdp": sum(r["gdp"] for r in sub), "co2": sum(r["co2"] for r in sub), "mean_green": sum(r["green"] for r in sub) / len(sub), "mean_resist": sum(r["resist"] for r in sub) / len(sub), "avoided_full": sum(r.get("avoided_co2_full", 0.0) for r in sub), } # 3.5 相关性（稳健性证据） corr = { "rank_drop_vs_resist": _corr([r["rank_drop"] for r in rows], [r["resist"] for r in rows]), "green_vs_resist": _corr([r["green"] for r in rows], [r["resist"] for r in rows]), "ci_vs_sens": _corr([r["ci"] for r in rows], [r["sens"] for r in rows]), } # 3.6 敏感性 tornado：对 ΔT(φ=1) 各 ±20% base_dT = scenarios[1.0]["dT"] sens_list = [] for label, fac in [("R_MAX", R_MAX), ("GAMMA_T", GAMMA_T), ("GLOBAL_SCALE_CO2", GLOBAL_SCALE_CO2), ("T", T)]: lo_dT, hi_dT = [], [] for mult in (0.8, 1.2): saved = {"R_MAX": R_MAX, "GAMMA_T": GAMMA_T, "GLOBAL_SCALE_CO2": GLOBAL_SCALE_CO2, "T": T} saved[label] = fac * mult # 重算 ΔT(φ=1) dCp = 0.0 for r in rows: r_phi = saved["R_MAX"] * 1.0 * r["sens"] dCp += r["co2"] * r_phi * saved["T"] dCg = dCp * saved["GLOBAL_SCALE_CO2"] dT2 = saved["GAMMA_T"] * dCg * PPM_PER_GT if mult == 0.8: lo_dT.append(dT2) else: hi_dT.append(dT2) sens_list.append({"param": label, "low": min(lo_dT + hi_dT), "high": max(lo_dT + hi_dT), "base": base_dT, "lo_rel": (min(lo_dT + hi_dT) - base_dT) / base_dT, "hi_rel": (max(lo_dT + hi_dT) - base_dT) / base_dT}) # 3.7 深析国：巴西 br = next(r for r in rows if r["id"] == DEEP_ID) br_r_full = br["r_full"] br_avoid = br["avoided_co2_full"] br_benefit = br_avoid * 1e9 * SCC br_cost = EPS_COST * global_gdp * 1.0 * 1e12 * (br["gdp"] / panel_gdp) # 巴西分摊 brazil = { "id": br["id"], "name": br["name"], "gdp": br["gdp"], "co2": br["co2"], "res": br["res"], "env": br["env"], "d_res": br["d_res"], "d_env": br["d_env"], "ggd": br["ggd"], "green": br["green"], "gdp_rank": br["gdp_rank"], "ggd_rank": br["ggd_rank"], "rank_drop": br["rank_drop"], "ci": br["ci"], "sens": br["sens"], "resist": br["resist"], "r_full": br_r_full, "avoided_full": br_avoid, "baa_emis": br["co2"] * (1 - br_r_full), # 稳态年排放 "benefit": br_benefit, "cost": br_cost, "net": br_benefit - br_cost, } # 3.8 值得性判定 sc1 = scenarios[1.0] worth_it = { "net_global": sc1["net"], "bcr": sc1["bcr"], "worth": sc1["net"] > 0, "resist_top": sorted(rows, key=lambda r: r["resist"], reverse=True)[:8], } D = { "N": len(rows), "rows": rows, "scenarios": scenarios, "by_region": by_region, "regions": regions, "corr": corr, "sens": sens_list, "brazil": brazil, "worth_it": worth_it, "panel_co2": panel_co2, "panel_gdp": panel_gdp, "global_gdp": global_gdp, "PHIS": PHIS, "PHI_CN": PHI_CN, "constants": {"T": T, "SCC": SCC, "R_MAX": R_MAX, "GAMMA_T": GAMMA_T, "PPM_PER_GT": PPM_PER_GT, "GLOBAL_SCALE_CO2": GLOBAL_SCALE_CO2, "GLOBAL_SCALE_GDP": GLOBAL_SCALE_GDP, "EPS_COST": EPS_COST}, "REGION_CN": REGION_CN, "DEEP_ID": DEEP_ID, } # 3.9 写 CSV（32 行国家面板） out = os.path.normpath(os.path.join(HERE, "..", "assets", "problems", "data", out_csv)) os.makedirs(os.path.dirname(out), exist_ok=True) header = ("id,name,region,gdp_trillion,co2_gt,res_share,env_share,dev,fossil_share," "ggdp_trillion,green_ratio,gdp_rank,ggd_rank,rank_drop,carbon_intensity," "resist_index,sensitivity,r_full,avoided_co2_full_gt") lines = [header] for r in rows: lines.append("%d,%s,%s,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.4f,%g,%g,%g,%.4f,%.4f,%.4f,%.4f,%.4f" % (r["id"], r["name"], r["region"], r["gdp"], r["co2"], r["res"], r["env"], r["dev"], r["fossil"], r["ggd"], r["green"], r["gdp_rank"], r["ggd_rank"], r["rank_drop"], r["ci"], r["resist"], r["sens"], r.get("r_full", 0.0), r.get("avoided_co2_full", 0.0))) with open(out, "w", encoding="utf-8-sig") as f: f.write("\n".join(lines) + "\n") D["csv"] = out return D if __name__ == "__main__": d = gen_icm2023f() print("ICM2023F Green GDP 生成完成 ->", d["csv"]) print("面板 N=%d 全球GDP(缩放)=%.2f 万亿$ 面板CO2=%.2f Gt/yr" % (d["N"], d["global_gdp"], d["panel_co2"])) print("\n=== 全球采纳情景(φ) ===") for phi in d["PHIS"]: s = d["scenarios"][phi] print(" φ=%.2f | ΔC(全球)=%.1f Gt | Δppm=%.2f | ΔT=%.4f°C | 效益$%.2e | 成本$%.2e | 净$%.2e | BCR=%.0f" % (phi, s["dC_global"], s["dppm"], s["dT"], s["benefit"], s["cost"], s["net"], s["bcr"])) print("\n=== 区域聚合(满采纳避免CO2) ===") for rg in d["regions"]: b = d["by_region"][rg] print(" %s: n=%d 均值绿色比=%.3f 均值抵抗=%.3f 避免CO2=%.2f Gt" % (rg, b["n"], b["mean_green"], b["mean_resist"], b["avoided_full"])) print("\n=== 相关性 ===") for k, v in d["corr"].items(): print(" %s = %.4f" % (k, v)) print("\n=== 敏感性 tornado (ΔT φ=1 基准=%.4f°C) ===" % d["scenarios"][1.0]["dT"]) for s in d["sens"]: print(" %s: [%.4f, %.4f] (±%.1f%% / +%.1f%%)" % (s["param"], s["low"], s["high"], s["lo_rel"]*100, s["hi_rel"]*100)) print("\n=== 深析国：巴西 ===") b = d["brazil"] print(" GDP=%.2f 万亿$ GGDP=%.3f 万亿$ 绿色比=%.3f" % (b["gdp"], b["ggd"], b["green"])) print(" GDP排名=%g GGDP排名=%g 排名下降=%g" % (b["gdp_rank"], b["ggd_rank"], b["rank_drop"])) print(" D_res=%.3f(资源耗减) D_env=%.3f(环境退化) 碳强度=%.4f Gt/万亿$" % (b["d_res"], b["d_env"], b["ci"])) print(" 满采纳年减排比例=%.3f 稳态年排放=%.3f Gt 累计避免=%.3f Gt" % (b["r_full"], b["baa_emis"], b["avoided_full"])) print(" 巴西气候效益$%.2e 分摊成本$%.2e 净$%.2e" % (b["benefit"], b["cost"], b["net"])) print("\n=== 值得性 ===") w = d["worth_it"] print(" 全球净收益$%.2e BCR=%.0f 值得转轨=%s" % (w["net_global"], w["bcr"], w["worth"])) print(" 高抵抗国(TOP8):", [r["name"] for r in w["resist_top"]]) print("\n=== 排名下降最大 TOP5 ===") for r in sorted(d["rows"], key=lambda r: r["rank_drop"], reverse=True)[:5]: print(" %s 下降%g (resist=%.3f)" % (r["name"], r["rank_drop"], r["resist"]))