# -*- coding: utf-8 -*- """ ICM 2019 D — Time to leave the Louvre（离开卢浮宫的时间） 真实赛题：COMAP 国际交叉学科建模竞赛（ICM）D 题——卢浮宫应急疏散模型。 T1 开发适应型应急疏散模型，使博物馆能探索疏散游客的一系列选项，同时允许急救人员尽快进入。 T2 模型须适应各类潜在威胁（火灾/化学/活跃威胁等），可动态重配置可用出口。 T3 确定向出口移动人流的潜在瓶颈（bottleneck）。 T4 验证模型，讨论卢浮宫如何实施；提出应急管理的政策与程序建议。 T5 讨论如何将模型调整并应用到其他类似大型拥挤建筑。 真源模型（纯标准库，确定性；SEED=2019 仅作版本锚定，模型无随机噪声）： 将卢浮宫抽象为「多层多核」网络，采用容量比例分配（capacity-proportional）的 流水线疏散模型（避免纯最短路导致的 Braess 式悖论）： · 5 层（B2/B1/RDC/L1/L2），每层 8 区，每 2 区归 1 个竖向疏散核心（core），共 4 核， 分别与 4 个主出口方向对应（金字塔/黎塞留/卡鲁塞尔/狮门）。 · 边容量 C = 宽度 × 最大舒适密度 ρ_max × 步行速度 V。 · 每个核心 k 的负载 L_k 沿竖向核心下放到地面，再由该核「开放出口集合」 （自身主出口若开放 + 自身辅助出口若开放）按容量比例分配： flow_{k,e} = L_k · C_e / Σ_{e∈open(k)} C_e 实现各出口同步均衡清空，出口清除时间 CT_e = flow_{k,e} / C_e = L_k / ΣC_k。 · 竖向核心清除时间 V_k = L_k / C_vert（序列在出口之前）。 · 最不利分区行程偏移 offset_k = 区→核 + 最不利层竖向行程（取最深海下层）。 · 全局疏散时间 T = max_k( offset_k + V_k + G_k )，G_k = L_k / ΣC_k 为核内出口均衡清除时间。 · 威胁场景：关闭受影响主出口 → 该核负载回落到辅助出口 → 重算 T，瓶颈随重配置转移。 · 辅助出口阈值：仅主线时瓶颈严重度高于有辅助时，量化开辅助出口的增益。 · 急救人员进入（FREM）：自黎塞留员工入口逆行进入、抵达最不利分区（狮门翼 B2）， 与疏散人流在同边逆向冲突 → 等待黎塞留出口边与狮门竖向核心清除后方能通行，量化到达延迟。 · 敏感性：占用 +30% / 步行速度 −20% / 通道宽度 −20%。 输出：assets/problems/data/icm2019d.csv（8 行情景面板 + 表头） 作者说明：本模型为方法论演示用的确定性合成网络（种子 2019），几何/容量参数取自公开疏散工程量级 （步行 1.25 m/s、ρ_max=1.5 人/m²、主通道宽 5–8 m），数值非卢浮宫实测蓝图； 真实参赛应以建筑 BIM/消防图纸与现场演练标定；本模型聚焦「容量比例分配-瓶颈-自适应- 急救」的可复现链路，结论（辅助出口显著降险、竖向核心是普遍瓶颈、化学威胁最不利）稳健。 """ import os import math import csv SEED = 2019 HERE = os.path.dirname(os.path.abspath(__file__)) # --------------------------------------------------------------------------- # 0. 模型常数（均作为假设显式声明，便于敏感性分析） # --------------------------------------------------------------------------- V_WALK = 1.25 # 普通游客步行速度 m/s（≈4.5 km/h） V_RESP = 1.50 # 急救人员步行速度 m/s（训练有素） RHO_MAX = 1.50 # 通道最大舒适人流密度 人/m² OCC_VIS = 16000 # 高峰同时在馆游客（确定性合成） OCC_STAFF = 2000 # 员工/安保（确定性合成） OCC = OCC_VIS + OCC_STAFF FLOORS = 5 FLOOR_W = [0.18, 0.20, 0.25, 0.22, 0.15] # 各层人数占比（B2/B1/RDC/L1/L2） FLOOR_NAME = {0: "B2 地下二层", 1: "B1 地下一层", 2: "RDC 地面层", 3: "L1 二层", 4: "L2 三层"} ZONE_W = [0.15, 0.15, 0.15, 0.15, 0.10, 0.10, 0.10, 0.10] # 每层 8 区权重(和=1) CORE_SHORT = ["金字塔", "黎塞留", "卡鲁塞尔", "狮门"] # 主出口（地面层 f=2 的核心出口）：节点名 / 对应核 / 通道长 m / 通道宽 m MAIN = [("E_pyramid", 0, 40, 8), ("E_richelieu", 1, 90, 6), ("E_carrousel", 2, 70, 6), ("E_lions", 3, 110, 5)] MAIN_NAME = {m[0]: m[1] for m in MAIN} # 辅助出口（员工/秘密/紧急，仅员工与消防知情）：6 个，分配到各核 AUX_N = 6 AUX_CORE = [0, 0, 1, 2, 2, 3] AUX_L = 70.0 AUX_W = 3.0 # 几何参数 L_ZC = 50.0 # 分区→核心 通道长 m W_ZC = 5.0 L_VERT = 50.0 # 竖向核心 层间长 m（梯段+平台） W_VERT = 8.0 # 分区人数（基线） ZONE_POP = {} for _f in range(FLOORS): for _i in range(8): ZONE_POP[(_f, _i)] = FLOOR_W[_f] * ZONE_W[_i] * OCC # 各核负载（基线） CORE_LOAD = {} for _k in range(4): CORE_LOAD[_k] = sum(ZONE_POP[(_f, _i)] for _f in range(FLOORS) for _i in range(8) if _i // 2 == _k) # 各层占用 FLOOR_OCC = {_f: sum(ZONE_POP[(_f, _i)] for _i in range(8)) for _f in range(FLOORS)} # --------------------------------------------------------------------------- # 1. 容量与行程工具 # --------------------------------------------------------------------------- def _cap(width, V): return width * RHO_MAX * V def _offset(V): """最不利分区行程偏移：区→核 + 最深海下层(2 层)竖向行程。""" return L_ZC / V + 2.0 * L_VERT / V # --------------------------------------------------------------------------- # 2. 单情景求解（容量比例分配流水线） # --------------------------------------------------------------------------- def _run(open_main, open_aux, occ_mult=1.0, v_mult=1.0, w_mult=1.0): _V = V_WALK * v_mult _L = {_k: CORE_LOAD[_k] * occ_mult for _k in range(4)} _C_vert = _cap(W_VERT * w_mult, _V) # 每个核的开放出口集合（容量） _open_exits = {} # k -> [(label, C_e)] _exit_cap = {} for _name, _k, _Lm, _Wm in MAIN: _exit_cap[_name] = _cap(_Wm * w_mult, _V) if _name in open_main: _open_exits.setdefault(_k, []).append((_name, _exit_cap[_name])) for _a in range(AUX_N): _k = AUX_CORE[_a] _nm = "A%d" % _a _exit_cap[_nm] = _cap(AUX_W * w_mult, _V) if _nm in open_aux: _open_exits.setdefault(_k, []).append((_nm, _exit_cap[_nm])) _off = _offset(_V) _T = -1.0 _Tk = {} _exit_flow = {} # 出口标签 -> 流量 _exit_ct = {} # 出口标签 -> 清除时间 _binding = None _binding_ct = -1.0 _severe = -1.0 for _k in range(4): _exits = _open_exits.get(_k, []) _sumC = sum(_c for _, _c in _exits) _Lk = _L[_k] if _sumC <= 0 or not _exits: # 核无可用出口（理论上不会发生，因辅助出口常开） _Tk[_k] = 1e9 _T = max(_T, 1e9) continue _Gk = _Lk / _sumC # 核内出口均衡清除时间 _Vk = _Lk / _C_vert # 竖向核心清除时间 _Tk[_k] = _off + _Vk + _Gk _T = max(_T, _Tk[_k]) # 出口级流量与清除时间 for _nm, _c in _exits: _fl = _Lk * _c / _sumC _ct = _fl / _c # = Gk（均衡） _exit_flow[_nm] = _exit_flow.get(_nm, 0.0) + _fl _exit_ct[_nm] = _ct if _ct > _binding_ct: _binding_ct = _ct _binding = _nm # 竖向核心也是潜在瓶颈 if _Vk > _binding_ct: _binding_ct = _Vk _binding = "核心%d竖向井(%s)" % (_k, CORE_SHORT[_k]) # 严重度（出口超载比 = 该出口实际吞吐 / 容量）：均衡下 = Gk*ΣC/Lk... 用 demand/cap for _nm, _c in _exits: _demand = _exit_flow.get(_nm, 0.0) / _T if _T > 0 else 0.0 _ov = _demand / _c if _ov > _severe: _severe = _ov _frem = _frem_delay(_exit_flow, _exit_cap, _C_vert, _L, _V) _occ_eff = OCC * occ_mult return dict(T=_T, Tk=_Tk, off=_off, exit_flow=_exit_flow, exit_ct=_exit_ct, C_vert=_C_vert, binding=_binding, binding_ct=_binding_ct, overload=_severe, occ_eff=_occ_eff, exit_cap=_exit_cap, L=_L, V=_V, frem=_frem) # --------------------------------------------------------------------------- # 3. 急救人员进入（FREM）延迟 # --------------------------------------------------------------------------- def _frem_delay(exit_flow, exit_cap, C_vert, L, V): # 自黎塞留员工入口逆行进入，抵达最不利分区（狮门翼 B2：核3） # 自身行程（逆行，训练速度 V_RESP） _own_len = 90.0 + 60.0 + 60.0 + 50.0 + 50.0 + 50.0 # 黎塞留边+2段环廊+2段竖井+区核 _own = _own_len / V_RESP # 冲突延迟：与疏散人流在同边逆向 _delays = {} # (1) 黎塞留主出口边：疏散者经此外出 _er = exit_flow.get("E_richelieu", 0.0) _delays["黎塞留主出口边"] = _er / exit_cap["E_richelieu"] if "E_richelieu" in exit_cap else 0.0 # (2) 狮门核竖向核心：疏散者(核3)下行，急救者上行逆向 _delays["狮门核心竖向井"] = L[3] / C_vert return dict(arrival=_own + sum(_delays.values()), own=_own, delays=_delays) # --------------------------------------------------------------------------- # 4. 情景定义 # --------------------------------------------------------------------------- _ALL_MAIN = [m[0] for m in MAIN] _ALL_AUX = ["A%d" % a for a in range(AUX_N)] SCEN = [ ("S0_主线_only", _ALL_MAIN, []), ("S1_全开(推荐)", _ALL_MAIN, _ALL_AUX), ("S2_火灾_地面出口受损", ["E_richelieu", "E_lions"], _ALL_AUX), ("S3_化学泄漏_主出口停用", [], _ALL_AUX), ("S4_活跃威胁_黎塞留翼封锁", ["E_pyramid", "E_carrousel", "E_lions"], _ALL_AUX), ("S5_高峰+30%占用", _ALL_MAIN, _ALL_AUX, 1.30, 1.0, 1.0), ("S6_步行速度-20%", _ALL_MAIN, _ALL_AUX, 1.0, 0.80, 1.0), ("S7_通道宽度-20%", _ALL_MAIN, _ALL_AUX, 1.0, 1.0, 0.80), ] # --------------------------------------------------------------------------- # 5. 主入口：生成数据 + CSV # --------------------------------------------------------------------------- def gen_icm2019d(seed=SEED, out_csv="icm2019d.csv"): _results = {} for _row in SCEN: _name = _row[0]; _m = _row[1]; _a = _row[2] _om = _row[3] if len(_row) > 3 else 1.0 _vm = _row[4] if len(_row) > 4 else 1.0 _wm = _row[5] if len(_row) > 5 else 1.0 _res = _run(_m, _a, _om, _vm, _wm) _results[_name] = dict(res=_res, open_main=list(_m), open_aux=list(_a), occ_mult=_om, v_mult=_vm, w_mult=_wm) # ---- 写 CSV（8 行情景面板 + 表头）---- OUT_DIR = os.path.normpath(os.path.join(HERE, "..", "assets", "problems", "data")) os.makedirs(OUT_DIR, exist_ok=True) csv_path = os.path.join(OUT_DIR, out_csv) cols = ["scenario", "occ_eff", "open_exits", "T_s", "T_min", "binding_edge", "binding_CT_s", "worst_overload", "frem_delay_s", "demand_rate_pers"] with open(csv_path, "w", newline="", encoding="utf-8-sig") as _f: _w = csv.writer(_f) _w.writerow(cols) for _row in SCEN: _name = _row[0] _r = _results[_name]["res"] _fr = _r["frem"] _w.writerow([ _name, round(_r["occ_eff"], 0), len(_results[_name]["open_main"]) + len(_results[_name]["open_aux"]), round(_r["T"], 1), round(_r["T"] / 60.0, 2), _r["binding"], round(_r["binding_ct"], 1), round(_r["overload"], 3), round(_fr["arrival"], 1), round(_r["occ_eff"] / _r["T"], 2) if _r["T"] > 0 else "", ]) # ---- 高亮速览（四路一致核心数字）---- _hl = {} for _row in SCEN: _name = _row[0] _r = _results[_name]["res"] _hl[_name] = dict(T_s=round(_r["T"], 1), T_min=round(_r["T"] / 60.0, 2), binding=_r["binding"], binding_CT=round(_r["binding_ct"], 1), overload=round(_r["overload"], 3), frem=round(_r["frem"]["arrival"], 1)) _s1 = _results["S1_全开(推荐)"]["res"] _s0 = _results["S0_主线_only"]["res"] _s3 = _results["S3_化学泄漏_主出口停用"]["res"] D = dict( SEED=SEED, V_WALK=V_WALK, V_RESP=V_RESP, RHO_MAX=RHO_MAX, OCC=OCC, OCC_VIS=OCC_VIS, OCC_STAFF=OCC_STAFF, FLOORS=FLOORS, FLOOR_W=FLOOR_W, CORE_SHORT=CORE_SHORT, MAIN=MAIN, AUX_N=AUX_N, AUX_CORE=AUX_CORE, CORE_LOAD=CORE_LOAD, FLOOR_OCC=FLOOR_OCC, results=_results, highlights=_hl, scen_names=[s[0] for s in SCEN], csv_path=csv_path, summary=dict( base_T_s=round(_s1["T"], 1), base_T_min=round(_s1["T"] / 60.0, 2), main_only_T_s=round(_s0["T"], 1), aux_gain_pct=round((_s0["T"] - _s1["T"]) / _s0["T"] * 100, 1), worst_threat="S3_化学泄漏_主出口停用", worst_threat_T_s=round(_s3["T"], 1), worst_threat_T_min=round(_s3["T"] / 60.0, 2), frem_base_s=round(_s1["frem"]["arrival"], 1), frem_mainonly_s=round(_s0["frem"]["arrival"], 1), ), ) return D if __name__ == "__main__": D = gen_icm2019d() s = D["summary"] print("并发在馆人数 OCC =", D["OCC"], "(游客 %d + 员工 %d)" % (D["OCC_VIS"], D["OCC_STAFF"])) print("基线(全开)疏散时间 T = %.1f s ≈ %.2f min" % (s["base_T_s"], s["base_T_min"])) print("仅主线 T = %.1f s；开辅助出口增益 = %.1f%%" % (s["main_only_T_s"], s["aux_gain_pct"])) print("最不利威胁(S3 化学) T = %.1f s (≈%.2f min)" % (s["worst_threat_T_s"], s["worst_threat_T_min"])) print("急救到达延迟：全开=%.1f s，仅主线=%.1f s" % (s["frem_base_s"], s["frem_mainonly_s"])) print("\n情景面板：") for nm in D["scen_names"]: h = D["highlights"][nm] fr = D["results"][nm]["res"]["frem"] print(" %-28s T=%.1fs(%.2fmin) 瓶颈=%s CT=%.1fs 超载=%.2f 急救=%.1fs" % (nm, h["T_s"], h["T_min"], h["binding"], h["binding_CT"], h["overload"], fr["arrival"])) print("\n各层占用：", {FLOOR_NAME[k]: round(v) for k, v in D["FLOOR_OCC"].items()}) print("各核负载：", {CORE_SHORT[k]: round(v) for k, v in D["CORE_LOAD"].items()}) print("CSV ->", D["csv_path"])