城市交通碳排放系统动力学模型：从因果回路到Python代码

简介面向城市交通能源消费与碳排放研究的学术资料以系统动力学模型为核心结合镇江市案例构建城市客运交通仿真框架完整呈现从模型设计、实证分析到政策建议的研究链条。适合交通规划、能源经济与环境政策方向的研究者、高校师生以及从事低碳交通工作的从业者参考。文档共1个PDF文件压缩包大小514KB内容为2017年中国城市交通规划年会论文全文兼顾方法综述与实例演算文中梳理了LEAP、终端能源消费模型等主流预测方法并对比宏观与微观模型的局限可辅助读者快速梳理系统动力学在城市交通能耗预测中的建模思路与关键参数。已有249人学习浏览相关论文写作或城市交通碳减排量化分析可从中获得可复用的模型框架、情景设置与政策参考兼具理论与应用价值。1. 系统动力学模型在城市交通能源与碳排放研究中的位置城市交通碳排放不是一条随GDP线性抬升的曲线而是由机动车保有量、路网容量、出行需求、燃料结构和政策反馈共同塑造的累积过程。传统回归模型把参数固定在历史均值上模型一旦外推超过五年就开始失真而系统动力学模型恰好擅长描述这种“存量累积—反馈延迟—非线性波动”的机制。它可以回答“如果新能源车渗透率在2030年达到40%而公交出行分担率提升10%城市交通碳排放会在什么时候达峰”这类跨期政策问题。本文以城市交通能源消费和碳排放为研究对象不用图表软件直接用Python差分方程搭建一套可运行的系统动力学模型从因果回路、变量边界、参数标定到情景模拟和敏感性分析一步步落到可复现的代码上。适合双碳规划、交通政策预研以及想在业务上把宏观问题翻译成微分方程的技术人。2. 先立理论交通能源碳排放系统的存量、流量与反馈结构2.1 为什么是系统动力学而不是回归或LEAP交通能源预测领域有三种常用路径统计回归、部门叠加、系统动力学。统计回归用GDP、人口、机动车保有量拟合历史油耗上手快但解释不了政策干预为什么会出现“先降后升”的反弹。LEAP一类自下而上模型会把每种车型的行驶里程、单位能耗、排放因子乘起来结果细致可是情景之间彼此隔离油价上涨、拥堵恶化、公交分担率提升这些因素无法在模型内部形成反馈。城市交通能源消费和碳排放恰好是一个反馈网络车辆越多路越堵车速下降导致单位里程能耗上升能源消费增加促使政府调整限购和公交政策而政策又反过来改变车辆增长率。这类跨年份的动态响应只有系统动力学的存量-流量-反馈结构能在同一时间轴上模拟出来。2.2 模型边界与变量分类先回答算的是谁的交通和谁的能源建模第一步不是写方程而是划边界。我一般把边界定为城市内客运交通包含私人小客车、公共汽电车和轨道交通不包含货运和城际出行。货运的驱动因子是货物周转量和物流效率反馈回路与客运完全不同强塞进来会让模型参数失去解释力。边界划清后核心存量只有两个机动车保有量和累计碳排放。其余都是流量或辅助变量例如年新增车辆、年报废车辆、出行需求、平均车速、单位能耗、新能源渗透率。下表给出最小变量集单位统一使用公制避免后续换算混乱。变量类型单位初值/标定方式机动车保有量存量辆车辆登记数据累计碳排放存量万吨CO20年新增车辆流量辆/年购买率×保有量年报废车辆流量辆/年报废率×保有量城市人口外生万人统计年鉴日均出行需求辅助万公里/日人均出行强度×人口拥堵指数辅助无量纲由保有量与道路容量计算平均车速辅助km/h由拥堵指数计算单位能耗辅助L/km或kWh/km车辆技术水平新能源渗透率辅助无量纲政策情景输入这个表是建模的锚点。当一个人拿到统计口径混乱的数据时先按这个表把变量归位后面推导方程才不会跑偏。2.3 把因果回路翻译成一组差分方程系统动力学的核心是因果回路。我抽两条主要的反馈链正向回路是城市人口增长和经济发展推动出行需求上升私家车保有量增加路网拥堵加剧平均车速下降单位能耗上升最终能源消费和碳排放增长负向回路是拥堵加剧后出行时间成本上升部分出行者转向公共交通或取消出行出行需求回落机动车净增速度被抑制。把这两条回路写成方程时存量更新和辅助变量计算有严格顺序。常用的差分形式是P_{t1} P_t (P_t × r_purchase - P_t × r_scrap) × ΔtD_t D_0 × (1 g)^t × (1 - η × C_t)C_t (P_t / R)^αV_t V_max / (1 e^{β(C_t - 1)})E_t a_0 / V_t a_1 × V_t^2其中P表示机动车保有量D表示出行需求C表示拥堵指数V表示平均车速E表示单位里程能耗R为等效道路容量α、β为弹性系数。这里最关键的是顺序先根据上一期的速率更新存量再用本期的存量计算辅助变量最后用辅助变量求流量并累积到累计碳排放。如果把辅助变量和存量在同一时刻互推方程就变成了代数环模拟会发散。3. 用Python搭建可运行的交通能源碳排放SD模型3.1 欧拉积分与循环结构系统动力学代码的最小骨架把上一章的差分方程变成代码我一般不直接上Vensim而是先用numpy做一阶欧拉积分。原因是模型探索阶段需要频繁修改方程结构脚本能最快暴露错误等结构稳定后再迁移到专业工具也不迟。欧拉积分的写法很简单用循环扫描时间步每个步长内先计算存量的净速率再更新存量。时间步长取1年但如果模型中存在数值较大的增长率要么把dt缩小到0.5年要么改用龙格库塔法否则会出现存量振荡甚至负值。3.2 30行Python代码从2020年到2060年的模拟下面是一个最小可运行模型只保留机动车保有量和累计碳排放两个存量其余都是辅助变量。代码使用pandas输出结果方便直接查看。import numpy as np import pandas as pd t_start, t_end, dt 2020, 2060, 1.0 years np.arange(t_start, t_end 1, dt) n len(years) P np.zeros(n) # 机动车保有量辆 CUM np.zeros(n) # 累计碳排放万吨CO2 P[0] 250e4 # 初始250万辆 # 核心参数 r_purchase 0.08 # 年度新车购买率 r_scrap 0.05 # 年度报废率 R 300e4 # 等效道路容量辆 alpha 0.8 # 拥堵指数弹性 v_max 60.0 # 自由流车速km/h beta 2.0 # 拥堵敏感系数 a0, a1 0.3, 0.002 # 能耗拟合系数 vmt_base 12000.0 # 基准年均行驶里程km charge 0.2 # 拥堵对出行需求的抑制系数 new_share np.linspace(0.05, 0.8, n) # 新能源渗透率线性上升 cf_fuel 2.31 # 燃油碳强度kg CO2/L cf_ele 0.6 # 电力碳强度kg CO2/kWh fe 8.0 # 燃油车百公里油耗L/100km for i in range(1, n): # 1. 先更新存量基于上一期的净速率 P[i] P[i-1] (P[i-1] * r_purchase - P[i-1] * r_scrap) * dt # 2. 再算辅助变量拥堵指数、平均车速、行驶里程 C (P[i] / R) ** alpha V v_max / (1 np.exp(beta * (C - 1.0))) M vmt_base * (1 - charge * C) # 3. 最后算流量和累计排放 fuel_use P[i] * (1 - new_share[i]) * M * fe / 100 nev_use P[i] * new_share[i] * M * 0.15 total_co2 (fuel_use * cf_fuel nev_use * cf_ele) / 1e7 CUM[i] CUM[i-1] total_co2 * dt df pd.DataFrame({年份: years, 保有量_辆: P, 累计碳排放_万吨: CUM}) print(df.tail())这段代码的关键在循环顺序先更新保有量存量再计算拥堵指数C和平均车速V之后才计算燃油使用与电力使用。new_share是政策输入向量在这个示例中被设置为从5%线性升到80%代表新能源车替代速度。注意单位换算fuel_use是升cf_fuel是kg/L两者相乘后除以1e7把单位统一成万吨CO2。如果直接用P[i]的万辆单位计算排放量级会差一万倍这是最常见的单位陷阱。3.3 参数表与业务含义哪些需要标定哪些是政策输入模型中并非所有参数都需要精确标定。有的参数来自统计数据有的是政策目标有的是反向拟合值。下面这张表列出了一组默认参数和它们的业务来源方便其他人接手时快速理解。参数默认值单位业务来源r_purchase0.081/年新车销售/保有量比值r_scrap0.051/年车辆报废年限倒数R300万辆道路里程折算容量alpha0.8无量纲拥堵弹性标定v_max60km/h城市限速与自由流速度beta2.0无量纲车速-拥堵曲线拟合charge0.2无量纲出行需求对拥堵弹性new_share线性无量纲双碳政策目标标定参数时r_purchase和r_scrap可以从车辆销售与注销数据中获得alpha和charge如果拿不到调研数据可以先设置一组合理区间做再一次敏感性分析不等得出精确实数。模型有六个以上待标定参数时不要相信人工调参必须用优化器。4. 情景模拟与敏感性分析让政策参数可调、可比较4.1 把情景定义为参数矩阵三种常见政策包实际研究中只跑一条基准曲线没有价值必须把政策组合成情景。常见做法是把run_sd封装成函数参数用dict传入这样同一个模型可以批量跑多个情景例如基准情景、低碳情景、强化政策情景。下面代码用run_sd返回累计碳排放便于比较。def run_sd(params): P np.zeros(n) CUM np.zeros(n) P[0] params[p0] for i in range(1, n): P[i] P[i-1] (P[i-1] * params[r_purchase] - P[i-1] * params[r_scrap]) * dt C (P[i] / params[R]) ** alpha V v_max / (1 np.exp(beta * (C - 1.0))) M vmt_base * (1 - charge * C) fuel_use P[i] * (1 - new_share[i]) * M * fe / 100 nev_use P[i] * new_share[i] * M * 0.15 total_co2 (fuel_use * cf_fuel nev_use * cf_ele) / 1e7 CUM[i] CUM[i-1] total_co2 * dt return P, CUM base {p0: 250e4, r_purchase: 0.08, r_scrap: 0.05, R: 300e4} low {p0: 250e4, r_purchase: 0.06, r_scrap: 0.06, R: 320e4} strong {p0: 250e4, r_purchase: 0.04, r_scrap: 0.07, R: 340e4} for name, params in [(基准, base), (低碳, low), (强化, strong)]: P, CUM run_sd(params) print(name, 2050年累计排放: %.2f 万吨 % CUM[np.where(years2050)[0][0]])这里的三个参数包分别代表基准情景维持现有购车和报废节奏低碳情景略微收紧新车购买并提高报废率强化情景再增加道路容量建设降低拥堵弹性。注意new_share仍然是全局变量也就是说新能源替代速度三个情景共用一条目标曲线。如果你想更细可以把new_share也放进params里按年份分段定义。4.2 敏感性分析用Morris抽样而不是靠猜系统动力学模型的参数多单变量摄动法无法捕捉参数之间的交互。Morris方法是一种全局敏感性分析方法适合在模型跑一次耗时几秒的探索阶段使用。它按“一次一步”的方式抽样用少量样本估算每个参数对输出的影响。SALib库实现了现成的Morris采样与分析方法可以在需要时用pip install SALib安装。from SALib.sample import morris from SALib.analyze import morris as morris_analyze problem { num_vars: 3, names: [r_purchase, alpha, charge], bounds: [[0.02, 0.1], [0.2, 1.5], [0.0, 0.5]] } param_values morris.sample(problem, N50, seed42) outputs [] for row in param_values: params base.copy() params[r_purchase], params[alpha], params[charge] row[0], row[1], row[2] P, CUM run_sd(params) outputs.append(CUM[-1]) Si morris_analyze(problem, param_values, np.array(outputs)) for name, mu_star, sigma in zip(problem[names], Si[mu_star], Si[sigma]): print(f{name}: μ* {mu_star:.3f}, σ {sigma:.3f})mu_star越大说明该参数对累计碳排放的影响越强烈sigma大说明该参数与其他参数存在明显交互。这个结果可以指导建模者把精力放在高影响参数上。例如如果r_purchase的mu_star远大于charge那说明控制车辆增长的政策比限行更有杠杆效应而不是反向。4.3 三处必查的坑反馈极性、步长与单位模型跑出结果后第一件事是检查反馈极性是否错误。例如拥堵指数必须随保有量上升而上升如果代码里写成C (R / P[i]) ** alpha模拟结果会出现排放先升后降的假反转而且不容易被抽查数据发现。第二坑是时间步长过大。增长率超过0.5时欧拉法在dt1年下会产生振荡。判断方法很简单把dt改成0.25重新跑一次如果累计排放曲线与dt1年的结果相差超过5%说明步长不够。第三个坑是单位不一致尤其是P以“万辆”为单位、M以“公里”为单位、cf以“kg/L”为单位时最终CO2需要做多次换算。我一般在写代码时把单位写进变量名比如P_veh、M_km_per_year即使变量名长一些也能减少换算出错。5. 用历史数据校准和蒙特卡洛区间扩展模型5.1 用最小二乘把模型拉回历史曲线任何系统动力学模型在做预测前都要先做历史校验。r_purchase、alpha、charge这些参数很多时候不能直接观测我一般用历史保有量和碳排放数据做参数校准。scipy.optimize.least_squares可以完成这个任务。先定义一个残差函数返回模型输出与统计数据的差再让优化器调整参数。from scipy.optimize import least_squares def run_sd_full(params): P, CUM run_sd(params) return pd.DataFrame({保有量: P, 累计碳排放: CUM}) def residual(theta): params base.copy() params[r_purchase], params[alpha] theta[0], theta[1] sim run_sd_full(params) return sim[保有量].to_numpy()[:len(hist)] - hist.to_numpy() hist np.array([250, 270, 291, 314, 340]) # 示例历史数据 res least_squares(residual, x0[0.08, 0.8], bounds([0.01, 0.1], [0.2, 2.0])) print(校准后 r_purchase%.3f, alpha%.3f % (res.x[0], res.x[1]))残差函数里sim[保有量]前五年要和历史数据长度一致。如果优化过后残差仍然很大先不要继续调参应该回看模型结构是不是缺少限购、出租车电动化、停车费弹性这些机制。参数补偿结构缺陷是拟合中最大的误区。5.2 蒙特卡洛区间把模型预测变成可决策的概率带点估计预测没有传递参数不确定性所以我在方案比选时会跑蒙特卡洛模拟。用拉丁超立方或普通均匀抽样生成200组参数运行模型后取5%和95%分位数作为区间。下面的代码展示了最简单的区间构造方法rng np.random.default_rng(42) samples rng.uniform( low[0.04, 0.5, 0.05], high[0.10, 1.2, 0.35], size(200, 3) ) results [] for row in samples: params base.copy() params[r_purchase], params[alpha], params[charge] row P, CUM run_sd(params) results.append(CUM[-1]) lo, hi np.percentile(results, [5, 95]) print(2050年累计碳排放区间: %.2f ~ %.2f 万吨 % (lo, hi))这种区间比单条曲线更能支持决策。比较基准情景和强化政策情景时不要只看两条均值线是否分开而要比较两个概率区间的重叠程度如果重叠面积大说明政策效果在参数不确定性下并不显著。用这一招可以很快筛掉那些“看起来很美”但鲁棒性差的政策方案。本文还有配套的精品资源点击获取