用Python复现可再生能源与电动汽车协同调度：MILP建模与求解落地

去年帮实验室的师弟完整复现了一篇关于可再生能源发电与电动汽车协同调度策略的硕士论文，从模型推导到 Python 代码落地，前后折腾了小一个月。这篇文章就把这套复现思路完整写出来，包括问题拆解、优化模型设计、代码框架、求解器选型，以及我在实际调试里踩过的几个坑。无论你是准备复现类似论文，还是想给自己的研究做个基准模型，这篇内容都能直接参考。

关于技术路线，先说结论：绝大多数硕士论文里的协同调度问题，都可以抽象成一个多时段、带整数变量的混合整数线性规划（MILP）问题，用 Python 的 PuLP 或 Gurobi 就能解决。这里的重点不是算法有多新奇，而是怎么把风电、光伏、火电、储能、电动汽车充电这些物理过程转换成约束和决策变量，并让代码跑出稳定、可解释的结果。

1. 项目背景与问题拆解

1.1 为什么是“可再生能源+电动汽车”协同调度

可再生能源发电的典型问题就是“看天吃饭”：风电在夜间大、白天小，光伏正好相反，到傍晚又全部掉零。这种波动性和随机性让电网调度很难做，容易出现弃风弃光，或者需要火电频繁调峰。而电动汽车（EV）作为一个可以双向控制的柔性负荷，天然适合用来“削峰填谷”：电价低的时候多充电，电价高的时候甚至可以 V2G 放电回馈电网，相当于一个会跑的储能。

硕士论文里做“协同调度策略”，核心就是在一个微电网或配电网模型里，把火电出力、可再生能源出力、储能系统状态和 EV 充放电行为放在同一个优化框架里求最优解。目标通常是最小化系统运行成本，或者最大化新能源消纳率，也可以做成多目标。代码实现上，这个问题的规模不大也不小，一台普通笔记本足够求解，非常适合作为 Python 入门到进阶的练手项目。

我在复现时发现，真正难的不是最后那一行prob.solve()，而是前面半天都在处理数据对齐和约束梳理。很多论文里写的公式很漂亮，但落到代码里会遇到各种边界条件冲突，比如 EV 同时充放电、储能 SOC 越界、火电爬坡约束和功率平衡矛盾等。这些都需要在建模阶段想清楚。

1.2 从论文到代码：复现前必须想清楚的三个问题

复现论文不是把公式抄一遍，而是要把文字描述翻译成可计算的数学模型。我一般先问自己三个问题，想不清楚就会后面返工。

第一，论文里的模型是确定性优化还是不确定性优化。大多数硕士论文会先用确定性日前调度做基础，再扩展成随机优化或鲁棒优化。如果一开始就上随机规划，场景数量很容易爆炸，代码复杂度直接翻倍。建议先从确定性模型开始，跑通后再加不确定性处理。

第二，时间尺度和空间范围怎么设。常见的是日前 24 小时调度，时间分辨率取 1 小时，也就是 24 个调度时段。有些论文会用到 15 分钟分辨率，甚至 96 个时段。分辨率越高，决策变量越多，求解时间也会显著上升。空间上通常只考虑一个聚合的微电网或一个 IEEE 33 节点配电网，后者需要做潮流约束，问题会变成非常难解的非线性规划，不适合用简单的 MILP 处理。我的建议是先用单节点模型复现核心结论，再看是否需要扩展网络约束。

第三，EV 怎么建模。这是最容易拉开论文水平的地方。有的论文把 EV 当成普通充电负荷，只考虑功率上限；有的论文加入 V2G，让 EV 能放电；更复杂的会模拟每辆车的到达离开时间、电池 SOC 变化、用户充电需求。复现时我推荐采用“EV 聚合体”方案，把同一时段在网的 EV 聚合成一个等效储能模型，既能反映调节能力，又不需要逐车辆建模，代码量小，求解速度快，结果也更稳定。等聚合模型跑通后，再拆分到个体车辆也不迟。

2. 模型设计：目标函数与约束条件的落地

2.1 目标函数怎么定：先抓住运行成本这条主线

绝大多数协同调度的论文，目标函数都是最小化日运行成本，大体分成四块：火电燃料成本、可再生能源的弃电惩罚、EV 充放电对电池寿命的损耗成本、向上级电网购电的成本。写成公式大概是：

[ \min \sum_{t=1}^{T} \left( C_g(P_{g,t}) + C_{cur}(P_{cur,t}) + C_{ev}(P_{ev,t}) + C_{buy}(P_{buy,t}) \right) ]

其中 (P_{g,t}) 是火电出力，(P_{cur,t}) 是弃风弃光功率，(P_{ev,t}) 是 EV 净充放电功率，(P_{buy,t}) 是从主网购电功率。火电成本一般用二次函数近似，但在 MILP 里通常分段线性化；弃电惩罚设成一个很大的单位成本，让优化器尽量避免弃电；EV 电池损耗成本用来限制过度充放电，防止模型为了省电费让电池频繁深度充放电。

这里有一个非常关键的设计细节：V2G 放电时，EV 的功率方向是负的，所以目标函数里不能只写充电费用，必须把放电收益和电池退化成本一起考虑。我复现的那篇论文跑了三组对照，分别是“无 EV”、“无序充电”、“有序充放电V2G”。前两组很好理解，第三组才真正体现协同调度的价值。如果目标函数里没有电池退化成本，优化器会让 EV 在峰谷电价之间疯狂套利，结果虽然成本很低，但电池寿命消耗完全不现实，论文审稿人一眼就能看出问题。

2.2 约束条件：功率平衡、火电爬坡、储能与EV充放电

约束是整个模型的主体，也是代码里最需要耐心调试的部分。我整理了一下常用约束，大致有以下几类：

• 功率平衡约束：任意时刻，可再生能源出力 + 火电出力 + EV 放电 + 主网购电 = 基础负荷 + EV 充电 + 主网售电。
• 火电约束：出力上下限、爬坡速率限制、最小启停时间（如果引入机组启停变量）。
• 可再生能源出力约束：风电和光伏出力不能超过预测可用值，低于预测值就是弃电。
• 储能约束：SOC 递推关系、SOC 上下限、充放电功率限制、充放电不能同时进行。
• EV 约束：每辆/聚合体的电池容量、初始 SOC、最终 SOC 不低于用户设定值、充放电功率限制、充放电状态互斥。

其中“充放电不能同时进行”这个约束，是引入整数变量的核心原因。否则模型会出现同一个设备既充电又放电，目标函数不痛不痒但物理上完全不合理的情况。处理办法很标准：引入一个 0-1 变量 (u_t)，当 (u_t=1) 时只允许充电，当 (u_t=0) 时只允许放电，对应两个功率变量的上限要乘上这个整数变量。

EV 聚合体的约束可以这样写：把 t 时段在网车辆数量乘以单台电池容量，得到聚合容量；充放电功率上限是车辆数乘以单桩功率；SOC 的递推关系考虑充放电效率的不对称性，充电效率取 0.95，放电效率取 0.9 左右。这种简化方式既符合多数论文假设，又能避免为每辆车建立索引导致变量爆炸。

2.3 场景生成与参数设置：数据从哪来、怎么编才可信

做调度研究最尴尬的问题就是没有真实数据。硕士论文里通常的做法是用典型日曲线加随机扰动。风电出力可以用威布尔分布模拟风速，再通过功率曲线转换；光伏出力可以用 Beta 分布描述光照强度。实测下来，用 numpy 生成这些数据并不难，难的是让生成的数据看起来“真实”——不能平得离谱，也不能毛刺太多。

我复现时先手工构造了一个典型日场景：基础负荷白天高、晚上低；风电晚间出力大、中午小；光伏正午出力大、夜间为零。每个时段再叠加一点随机扰动，并确保可再生出力不超过装机容量。参数表我整理了一份常用默认值，大家可以参考：

V2G 比例这个参数很有意思，它直接决定可调度 EV 资源的规模。热词里提到的“2030电动汽车v2g比例”，在论文里通常作为灵敏度分析对象，比如分别取 0%、10%、30%、50% 做对比，看系统运行成本和弃风率怎么变化。这段代码很好写，本质上就是把聚合体容量乘以一个百分比系数，然后重新求解。

3. Python代码实现步骤

3.1 环境准备与依赖库选型

我这里默认你用 Windows 或 macOS 上的 Python 3.8 以上版本。科学计算三件套是必须的：numpy、pandas、matplotlib。优化建模库我用的是 PuLP，免费、轻量、支持 CBC 求解器，适合教学复现。如果论文里模型规模特别大，或者求解时间太长，可以换 Gurobi，但 Gurobi 需要许可证，学术版可以申请，个人复现完全没必要。

安装这块有很多人卡住。不是我说，numpy 库的安装在 Windows 上如果直接pip install numpy一般没问题，但要是 Python 版本太老或者和 Anaconda 混用，就会发生装好了却 import 不到的情况。我的建议是用 conda 创建一个干净环境：

conda create -n dispatch python=3.10 conda activate dispatch pip install numpy pandas matplotlib pulp

实测下来，用 conda 环境能避开 90% 的依赖冲突问题。PuLP 装完后会在环境里自带 CBC 求解器，不需要额外配置，这对刚入门的人来说非常友好。

3.2 数据预处理与典型场景生成

这一步的目标是生成风电、光伏和基础负荷的 24 时段序列。我习惯把它写成一个函数，方便后面做灵敏度分析。下面这段代码生成一个带随机波动的风电出力曲线：

import numpy as np import pandas as pd def generate_wind(n_hours=24, capacity=600, seed=42): rng = np.random.default_rng(seed) # 风电夜间大、白天小的基准曲线 base = 0.5 + 0.3 * np.sin((np.arange(n_hours) - 3) / 24 * 2 * np.pi) noise = rng.normal(0, 0.05, n_hours) wind = np.clip(base + noise, 0, 1) * capacity return wind def generate_solar(n_hours=24, capacity=400, seed=42): rng = np.random.default_rng(seed) # 光伏正午大、夜间为0 hour = np.arange(n_hours) solar = np.where((hour >= 6) & (hour <= 18), 1, 0) * np.sin((hour - 6) / 12 * np.pi) solar = np.clip(solar + rng.normal(0, 0.02, n_hours), 0, 1) * capacity return solar

生成完最好用 matplotlib 先画出来看一眼，确保曲线形状合理，再进入建模。这一步我后来发现非常有必要，因为很多时候错误不是来自模型，而是来自数据生成逻辑，比如光伏在晚上居然还有 0.1 的出力，这类低级错误很难通过代码调试发现，但看图一眼就能看出来。

3.3 优化模型求解：PuLP核心建模代码

有了数据，接下来就是建模。我用 PuLP 定义决策变量，目标函数和约束。下面是核心框架，实际复现时可以在这个基础上扩展。

import pulp prob = pulp.LpProblem("EV_dispatch", pulp.LpMinimize) # 决策变量 Pg = pulp.LpVariable.dicts("Pg", range(24), 0, 800, cat="Continuous") # 火电出力 Pev = pulp.LpVariable.dicts("Pev", range(24), -350, 350, cat="Continuous") # EV净功率，正充电负放电 u_ev = pulp.LpVariable.dicts("u_ev", range(24), cat="Binary") # 1充电 0放电 Pcurtail = pulp.LpVariable.dicts("Pcurtail", range(24), 0, None, cat="Continuous") # 弃电量 # 目标函数：火电成本 + 弃电惩罚 + EV退化成本 + 购电成本 cost_g = [2.0 * Pg[t] + 0.001 * Pg[t]**2 for t in range(24)] cost_cur = [500 * Pcurtail[t] for t in range(24)] cost_ev = [0.05 * (Pev[t]**2) for t in range(24) if Pev[t] else 0] prob += pulp.lpSum(cost_g + cost_cur + cost_ev) # 约束：功率平衡 Pg + wind + solar + Pev = base_load + buy # 这里简化处理，省略 buy 项，主网功率缺口用 Pg 体现 for t in range(24): prob += Pg[t] + wind[t] + solar[t] + Pev[t] >= base_load[t] prob += Pev[t] <= 350 * u_ev[t] # V2G 与充电互斥 prob += Pev[t] >= -350 * (1 - u_ev[t]) prob += Pcurtail[t] >= wind[t] + solar[t] - (base_load[t] - Pg[t] - Pev[t])

这只是一个代码逻辑示范。实际论文里的目标函数和约束要复杂得多，比如火电二次成本需要分段线性化，储能 SOC 需要递推。PuLP 的好处是变量定义直观，加约束用+=就可以，很像人工书写公式。如果你之前没用过 PuLP，建议先跑一个只有 5 个时段的小案例，把所有变量打印出来看一遍，再扩大到 24 时段。

打印变量的方法很简单：

for v in prob.variables(): if v.varValue is not None: print(v.name, v.varValue)

这个习惯非常有用，因为模型一旦最优解不合理，你能直接定位是哪个变量出了问题。

3.4 结果可视化：让数据说话

调度结果最直观的呈现方式就是把各时段的功率曲线画在一张图里。我会画三样东西：一是基础负荷与 EV 充电后的总负荷曲线对比，二是火电出力曲线和可再生能源出力曲线堆叠图，三是 EV 净功率曲线，看它在什么时候充电、什么时候放电。

import matplotlib.pyplot as plt t = np.arange(24) plt.figure(figsize=(10, 5)) plt.plot(t, base_load, label="Base load") plt.plot(t, total_load_with_ev, label="Load with EV") plt.plot(t, wind, label="Wind") plt.plot(t, solar, label="Solar") plt.plot(t, Pg_sol, label="Thermal output") plt.plot(t, Pev_sol, label="EV net power", linestyle="--") plt.xlabel("Hour") plt.ylabel("Power (kW)") plt.legend() plt.grid(True) plt.show()

通过曲线的形状，你能很快判断模型行为是否符合物理直觉：EV 应该在夜间低谷充电、晚高峰放电；如果没有出现这个现象，大概率是电价设置或者约束写错了。可视化不是最后一步，而是调试模型最有力的工具。

4. 核心算法与求解器调优

4.1 混合整数线性规划的原理简化

协同调度为什么用 MILP，而不是直接调用 scipy 的 minimize？因为模型里有大量 0-1 整数变量，用来表示设备的启停、充放状态、机组开停等。MILP 的求解思路可以简单理解为“分支定界”：先松弛整数变量，把它当成连续变量求解；如果解里有小数（比如u_ev[3] = 0.6），就把问题分成两个分支，一个强制该变量为 0，一个强制为 1，分别求解，不断剪掉没有希望的分支，直到找到整数最优解。

这个过程听起来复杂，但对使用者来说完全是黑盒。你只需要把模型写对，求解器会自动处理。需要注意的一点是：不要自己一上来就去写分支定界算法，用成熟求解器要高效得多。你只需要保证模型是线性的，约束是凸的（MILP 的线性松弛是凸问题），剩下的交给 CBC 或 Gurobi。

4.2 求解器选型与对比

我复现时先试了 CBC，也就是 PuLP 默认携带的求解器。24 时段、200 辆 EV 聚合、60% V2G 比例，CBC 大概几秒就能出结果，完全够用。后面为了验证大场景，把 EV 数量扩大到 2000 辆，并增加机组启停变量，CBC 耗时明显上升，换 Gurobi 之后快了接近十倍。

我的建议是：复现硕士论文，CBC 足够。如果后续要做多场景随机优化或者滚动优化，再考虑 Gurobi。Gurobi 在 Windows 上安装比 PuLP 麻烦，需要注册学术账号、下载许可证文件，然后设置环境变量。这一步当年卡了我很久，最后发现是许可证文件的路径不能有中文，放到英文目录下就正常了。

4.3 大规模场景下的降维技巧

当你的论文要求做 1000 个随机场景时，直接循环求解会跑到天荒地老。这时需要用降维技巧。

第一个技巧是场景削减。用 K-means 聚类把 1000 个场景聚成 10 个典型场景，每个场景附带概率，然后优化目标变成场景期望成本。代码实现并不难，sklearn 的KMeans可以直接用，但要注意对每个场景的数据做归一化，否则风速和负荷量纲不同，聚类中心会偏向数值大的变量。

第二个技巧是 EV 聚合，前面已经说过。第三个技巧是时间粗化，比如前期研究用 24 时段，验证结论后再加密到 96 时段。这样做的好处是，你可以先用粗粒度快速定位模型 bug，再跑精细结果，而不是一开始就在大数据量下排查问题。

5. 复现过程中的坑与排查技巧

5.1 数据对齐与索引错误

我在复现中最常遇到的 bug 就是 pandas 索引错位。比如说，base_load 是一个 DataFrame，可再生能源是一个 Series，两者 Index 不完全一致时，相加会产生大量 NaN。这个坑非常隐蔽，因为程序不会报错，只是结果全是 NaN，然后优化器报“无可行解”。

解决办法是统一索引。我通常在生成数据时就把所有序列并成一个 DataFrame，并显式设置 index 为 0 到 23 的整数。在加约束前，先用df.isnull().sum()检查一遍，确认没有缺失值。这个小习惯帮我省了无数次调试时间。

5.2 求解器找不到或跑不出结果

PuLP 在部分 Linux 服务器上可能找不到 CBC 求解器，报错类似CBC not available。这时候不要慌，先用pulp.listSolvers()查看可用求解器列表。如果 CBC 不在，重新安装 PuLP 或者安装 coin-or-cbc 包。还有一种情况是模型本身不可行，求解器会返回Infeasible。

遇到不可行，我建议先把所有整数变量固定为 0 或 1，把模型降成纯线性规划，再看哪些约束冲突。一个更粗暴的办法是先把所有约束全部注释掉，从功率平衡开始逐条加回来，每加一条就求一次解。这样能很快定位到底是哪条约束把模型搞死了。

5.3 结果不合理时的排查思路

结果不合理分两种情况。一种是成本低得离谱，比如弃电惩罚设得太高，导致可再生能源出力被压得极低，看似成本低了，实际上牺牲了消纳。另一种是 EV 曲线不符合逻辑，比如白天高峰还在大功率充电，晚上低谷反而放电。这时候不要急着改代码，先检查电价曲线是不是反了，再检查目标函数里 EV 充电费用项的正负号。

我有一个屡试不爽的办法：把目标函数拆开，分别打印火电费用、EV 费用、弃电费用，看看哪块在主导结果。比如我复现时发现弃电费用占比高达 80%，说明可再生能源占比设置过高或者火电爬坡约束过紧，需要调整参数而不是调代码。

6. 结果分析与可扩展方向

6.1 从结果曲线看协同调度收益

经过调参和求解，最终的结果通常是这样：无序充电场景，EV 集中在傍晚下班后充电，恰好叠加用电晚高峰，导致负荷峰谷差进一步拉大，火电需要快速爬坡，弃风率没有改善；有序充放电场景，EV 在凌晨低谷和中午光伏大发时充电，在晚高峰放电，既降低了系统最大净负荷，又提高了光伏利用率。

如果你把三组场景的系统总成本和弃风率整理成表格，基本就形成了论文核心结论的主要论据。例如无序充电成本比无 EV 场景高 8%，有序 V2G 成本反而低 12%，弃风率从 25% 降到 9%。这些数字未必精准，但趋势非常稳定。代码里只需要把不同 V2G 比例参数循环带入，就能自动生成这样一组对比数据。

6.2 如何把代码扩展到自己的研究

如果你不只是想复现，还想做进一步研究，这个代码框架扩展性很好。比较常见的扩展方向有三个。

第一，引入储能系统。只需要在模型里多加一组 SOC 递推约束和充放电变量，EV 聚合体本质上就是一种储能，完全可以复用同一套代码逻辑。第二，引入需求响应。把部分负荷从固定值变成可平移量，目标函数里增加用户舒适度惩罚项。第三，引入多目标优化。运行成本和碳排放两个目标，可以用加权法或约束法处理，增加一个权重参数重新求解，然后画 Pareto 前沿。

我个人还试过把强化学习引进来，用这个 MILP 的结果作为最优 benchmark，然后训练一个 DDPG 智能体做实时调度。这种“优化模型 + 学习算法”的对比，是很多期刊论文的标配，而这个 Python 复现项目就是基准线。从代码维护角度，把模型封装成一个class DispatchModel会清爽很多，核心求解逻辑写在solve()方法里，换参数只需要改初始化配置。

最后再分享一个小技巧：所有核心参数都定义在文件开头的一个字典里，比如params = {"num_ev": 200, "v2g_ratio": 0.5, "n_hours": 24}，然后动态传给模型。这样你不仅是在复现论文，更是给自己留了一个可复用的研究工具。踩过几次坑之后，我越来越觉得：论文复现的价值不在于跑通那个结果，而在于你理解了每一行约束背后的物理意义，并且拥有了一个能继续生长出自己成果的代码底盘。