可再生能源与电动汽车协同调度策略论文复现：建模、求解与代码实现

复现过这篇论文的朋友应该都有同感：题目里“可再生能源发电”“电动汽车”“协同调度策略”每一个词都是热点，组合在一起却是个硬骨头。新能源出力的随机性怎么刻画，EV集群的充放电行为怎么建模，双边的“协同”到底协同什么，目标函数怎么写才能让审稿人挑不出毛病——光是这几个问题就能劝退不少人。更别说还要用Matlab和Python两套语言各写一遍，工作量直接翻倍。

我这篇文章就把这个方向完整拆开讲透，从数学模型到代码架构，从场景削减到求解器配置，把复现过程中踩过的坑和验证过的技巧一并整理出来。适合正在做电力系统优化调度方向的硕士生、刚入门论文复现的研一同学，以及想从Matlab迁移到Python做科研的工程师。你不需要提前学过很多知识，但如果你用过Yalmip或者pulp里的任何一个，上手会快很多。

1. 论文复现的整体思路与系统建模

1.1 从标题读懂论文的三层结构

凡是标题里带“协同调度策略”的论文，本质上都在讲一件事：让不同类型的灵活性资源互相配合，在满足系统安全运行的前提下，实现经济性或低碳性最优。放到这个题目里，灵活性资源一头是可再生能源（风电、光伏），另一头是电动汽车。

先把这三层结构拆开看。第一层是新能源发电侧建模：风速服从威布尔分布，光照强度服从贝塔分布，风机出力和光伏出力都有各自的非线性转换关系，这是整个模型的随机性来源。第二层是负荷侧或者说EV侧建模：每辆EV有接入时间、离开时间、起始SOC、目标SOC、充电功率上限，成百上千辆EV聚合在一起形成集群，集群的总充电能力随时段变化。第三层才是“协同”——在同一个优化框架里同时调度火电机组出力、可再生能源实际出力和EV的充放电功率，让系统总成本最小。

这里有一个多数初学者容易理解偏的地方：“协同”不是先算风电预测、再算EV充电计划，而是把两个问题放进同一个目标函数里联立求解。风电大发时，如果系统调峰能力不足，就会出现弃风，所以要让EV尽量在这个时段多充电；EV大规模集中充电时，又会让负荷曲线出现新的尖峰，所以要让充电时间避开原始负荷高峰。这个双向耦合关系，就是协同调度的核心价值所在。

我在复现这类论文时习惯先画一张信息流图（用纸笔就行）：哪些量是已知参数，哪些量是决策变量，哪些约束把不同模块串在一起。画完这张图，论文的模型就清晰了一大半。

1.2 核心数学模型：目标函数与约束条件

这个方向的论文，目标函数再怎么变，都逃不出一个基本框架：火电机组的燃料成本和启停成本 + 弃风弃光惩罚 + EV调度费用（如果允许V2G，就是充放电收益）。写成数学形式就是一个最小化问题：

[ \min \sum_{t=1}^{T} \sum_{g=1}^{N_g} \left( a_g P_{g,t}^2 + b_g P_{g,t} + c_g u_{g,t} + SUC_g \cdot y_{g,t} \right) + \sum_{t=1}^{T} \lambda_w \cdot P_{w,curt,t} + \sum_{t=1}^{T} \lambda_{ev} \cdot \left( P_{ev,ch,t} + P_{ev,dis,t} \right) ]

这里面每个符号都有实际工程含义：( P_{g,t} )是火电机组g在t时段的出力，( u_{g,t} )是启停状态（0/1变量），( y_{g,t} )是启动动作指示变量，( P_{w,curt,t} )是弃风量，( \lambda_w )是弃风惩罚系数。二次项( a_g P^2 )描述煤耗特性，但MILP求解器不直接吃二次函数，所以复现时一般用分段线性化处理，这个细节后面专门讲。

约束条件分四类。功率平衡约束是最硬的一条：所有发电出力之和必须等于负荷加EV充电负荷。火电约束包括出力上下限、爬坡速率限制和最小启停时间——爬坡约束是调度模型里最容易把系统搞到无解的约束之一。新能源约束是实际出力不能超过预测可用出力，弃风量就是这个差值。EV约束包括集群总充电功率上限、每辆EV的SOC动态方程、离网时SOC要达到用户设定值。

这些约束看起来多，但写成代码其实就几行的事情。真正的难点在于量纲和单位——风电出力单位是MW，EV充电功率单位是kW，如果图省事直接相加，结果一定是错的。复现的时候我建议全部统一到MW，EV集群聚合之后用MW表示，精度完全够用。

1.3 模型为什么落地为MILP而不是启发式算法

很多同学会问：既然目标函数里有二次项，为什么不用粒子群或者遗传算法？这里要说清楚一个关键问题：论文复现和实际工程出结果，首要要求是“结果可以解释、误差可控”。MILP模型有严格的数学最优性证明，求解器会给你一个最优间隙（gap），你知道当前解离全局最优还有多远。启发式算法跑十次得到十个不同的解，你没法判断哪个更优，审稿人也没法接受这种不确定的结果。

我复现过一篇用粒子群做协同调度的论文，跑了20次，最好解和最差解相差8%以上，而同一模型用Cplex只需要几十秒就能给出gap小于0.1%的全局最优解。这就是MILP的核心优势。

当然，MILP也有代价：求解时间随着整数变量数量增加而急剧上升。火电启停变量是0/1变量，机组数量×时段数就是整数变量的规模。10台机组、24个时段就是240个整数变量，对现代求解器来说是小菜一碟。但如果把每辆EV是否充电也做成0/1变量，几千辆EV就直接把模型压垮了。解决方法是EV聚合——把同一接入时段、同一离开时段的EV归为一个集群，集群层面只约束总功率和总能量，这样就把几千个整数变量压缩到十几个连续变量，模型规模天差地别。

2. 新能源出力不确定性与场景削减

2.1 风光出力模型的概率描述

风速概率分布通常用两参数威布尔分布描述，形状参数k和尺度参数c决定了分布的整体形态。实际风电场的数据拟合出来，k通常在1.8到2.3之间，c在6到9之间。风机出力与风速的关系是分段函数：切入风速以下不出力，额定风速到切出风速之间满发，中间段近似线性。

光伏出力模型相对简单，核心输入是辐照度( G )，标准测试条件（STC）下辐照度是1000W/m²，实际出力按比例缩放，再考虑温度修正系数。但光伏模型有个容易忽略的点：云层遮挡导致的辐照度突变在时序上具有很强的自相关性——前一时刻是多云，后一时刻大概率还是多云。用独立抽样生成辐照度序列是不对的，必须考虑时间相关性。

处理这个问题的标准做法是用一阶马尔可夫链：把辐照度分成若干状态（比如晴天、多云、阴天），状态转移矩阵描述前后时刻的演变规律，这样生成的场景序列在时间维度上就是连续的、符合物理规律的。我复现时用了一阶马尔可夫链加核密度估计来生成风速和辐照度场景，效果比直接蒙特卡洛独立抽样好很多。

2.2 蒙特卡洛抽样与K-means场景削减

随机优化的标准流程是：先根据概率分布抽样生成大量场景（比如1000个），再用场景削减算法筛出少数有代表性的场景（比如10个），每个场景带一个概率权重，最后把确定性问题变成带场景权重的期望值优化问题。

场景削减最常用的是K-means聚类。直接把每一天的风速曲线和光照曲线拼成一个向量，然后做K-means聚类，把1000条曲线聚成10类，每类的质心就是一个代表场景，该类中场景数量除以总数就是该场景的概率。这里有个实操技巧：聚类之前一定要做归一化，不然风速数值大（个位数到十几），光照数值也大，但两个向量的量纲不同会影响距离计算。归一化之后，每个维度的贡献才均匀。

另外一个更精细的方法是快速前向选择法（Fast Forward Selection），它基于 Kantorovich 距离逐次删除对概率分布影响最小的场景。相比K-means，前向选择法在电力系统场景削减里理论依据更扎实，算法实现也不复杂，几十行代码就能搞定。我在复现时两种方法都试了，统一到10个场景的情况下，前向选择法生成的场景多样性略好，K-means生成的场景更“典型”。如果论文里没有特别指定，用K-means完全够用，而且好向别人解释。

使用场景法的时候，我把优化问题改写成：

[ \min \sum_{s=1}^{S} \pi_s \cdot f(x, \xi_s) ]

其中( \pi_s )是场景s的概率权重，( \xi_s )是该场景下的风电光伏预测出力，( x )是决策变量。注意：第一阶段的决策变量（机组启停、EV调度计划）在所有场景下是同一个值，这是随机优化的核心约束——你必须在看到真实风电出力之前做出决策。这个“非预期性约束”在代码里不显式存在，但建模时一定要想清楚哪些变量是场景独立的，哪些变量是场景相关的。

2.3 随机优化与鲁棒优化的选型考虑

复现这个方向时，你还会看到很多论文用鲁棒优化替代随机优化。鲁棒优化的思路是：不假设场景的概率分布，只给定不确定量的变化区间，然后求解在最坏情况下的最优解。好处是不需要概率信息、求解结果非常保守可靠，坏处是过于保守——为了抵御最坏情况，系统会预留大量旋转备用，经济性明显下降。

在实际复现中，我给的选型建议是：如果论文的核心贡献在协同调度策略本身，用场景法（随机优化）就够了，模型简单、求解快、结果容易解释；如果论文想强调应对极端天气的能力，再上鲁棒优化，但要做好经济性变差的准备。也有论文做两阶段鲁棒优化，第一阶段确定机组启停，第二阶段在不确定集内做经济调度，这种模型用列与约束生成算法（C&CG）求解，复杂度明显上升——如果不是为了冲顶刊，不建议新手一上来就啃这个。

如果你用的是Matlab+Yalmip，鲁棒优化可以通过 Yalmip 的不确定变量描述（uncertain）配合求解器做鲁棒对等转换，但实际效果依赖求解器支持程度，我实测下来还是自己手动写对等模型更可控。

3. Matlab代码实现：从数据到结果的全流程

3.1 代码架构设计：四段式结构

Matlab复现时我建议把代码拆成四个文件或者四个清晰分隔的区域：数据准备区、模型构建区、求解配置区、结果分析区。很多新手喜欢把数据计算和模型约束混在一起写，看起来每一步都连着，实际上调试起来欲哭无泪——想改一个约束条件，得在一大段代码里反复寻找。

第一个区域负责生成或读入所有已知参数：机组参数（容量、煤耗系数、爬坡率、最小启停时间）、负荷曲线（24小时）、风速和光照场景（从2.2节的场景削减结果读入）、EV集群参数（各时段可用EV数量、充电功率上限、SOC约束）。第一区域跑完，工作区里全是清晰的变量名。

第二个区域用Yalmip定义决策变量和目标函数，拼接约束条件。第三区域配置求解器参数并调用optimize。第四区域提取并可视化结果。这种四段式结构最大的好处是：你改了一个参数（比如上调弃风惩罚系数），只需要重新运行第一区域到第三区域，结果分析区域完全不用动。

3.2 Yalmip求解MILP的关键写法

用Yalmip建模MILP的标准写法是：

%% 决策变量定义 P_g = sdpvar(n_gen, T); % 火电出力，连续变量 u = binvar(n_gen, T); % 启停状态，0/1变量 P_w = sdpvar(n_w, T); % 风电实际出力 P_ev = sdpvar(n_ev_cluster, T); % EV集群充放电功率，正为充电 %% 目标函数（二次煤耗函数使用分段线性近似） objective = 0; for t = 1:T for g = 1:n_gen % 这里用分段线性系数矩阵 c_seg 替代二次项 objective = objective + c_seg * P_g(g,t) + startup_cost * y(g,t); end end %% 约束拼接 Constraints = []; for t = 1:T Constraints = [Constraints, sum(P_g(:,t)) + sum(P_w(:,t)) ... == P_load(t) + sum(P_ev(:,t))]; % 功率平衡 Constraints = [Constraints, ... P_g(:,t) >= P_gmin .* u(:,t), ... P_g(:,t) <= P_gmax .* u(:,t)]; % 出力上下限 end %% 求解 ops = sdpsettings('solver', 'cplex', 'verbose', 2, 'savesolveroutput', 1); result = optimize(Constraints, objective, ops);

这段代码里有几个关键细节要特别注意。

第一，sdpvar定义连续变量，binvar定义0/1变量，两种变量拼接同一个约束时没有问题，但千万不能把binvar写成sdpvar再加整数约束——性能会差很多。Yalmip会自动识别变量类型并交给求解器的对应分支处理。

第二，约束的拼接方式[Constraints, constraint1, constraint2]是Yalmip的标准用法，效率没问题。但如果你有成千上万条约束，建议用cell数组先存起来最后一次性拼接，能显著减少Yalmip内部的符号计算开销。

第三，目标函数里的二次项。Yalmip本身就支持二次目标函数，Cplex也支持二次约束规划（MIQP），但求解速度比MILP慢很多。复现时建议把煤耗的二次曲线做分段线性化处理，切成4到5段就足够精确——这是我从英文论文里学到的经验，用4段线性化拟合二次函数，误差不到0.5%，但求解速度提升了好几倍。

3.3 非线性约束的处理技巧

EV的SOC动态方程本身是线性的：( SOC_{t+1} = SOC_t + P_{ch,t} \cdot \eta \cdot \Delta t / C_{bat} - P_{dis,t} \cdot \Delta t / (C_{bat} \cdot \eta) )，这个直接写成线性约束没问题。真正麻烦的是EV充电功率是双向的——一辆EV既可能充电也可能放电（V2G），如果你允许V2G，充电功率( P_{ev} )就是一个可正可负的连续变量，这时候目标函数里对( P_{ev} )取绝对值（充电费用/放电收益）就变成非线性的了。

处理办法是引入两个非负变量( P_{ch} )和( P_{dis} )，约束( P_{ev} = P_{ch} - P_{dis} )和( P_{ch}, P_{dis} \geq 0 )，再把目标函数里的( |P_{ev}| )替换成( P_{ch} + P_{dis} )。这样既保证了线性性，又通过目标函数的系数自然避免“同时充电又放电”的无意义解——因为同时充放会让目标函数多付一倍的调度费用，最优解一定会选择只充或者只放。

3.4 结果可视化与灵敏度分析

解出来之后，可视化决定了你的论文图能不能让导师眼前一亮。我习惯画三张图：

第一张是24小时功率平衡图：火电、风电、光伏的出力堆叠图，最上面叠上负荷曲线，EV充电功率用虚线画在负荷上方。第二张是EV集群的SOC变化曲线，展示调度策略下EV的充放电行为是否合理。第三张是弃风率对比柱状图——无EV调度、有序充电、V2G三种模式下的弃风率，直观体现协同调度的效果。

画图用Matlab基础绘图函数就能实现，但要注意堆叠图用area函数，柱状图用bar函数，配色尽量用色盲友好的方案（如cmocean或wong配色）。顺便说一句，论文图的字体大小至少10pt，线宽至少1.5pt，否则放到Word里会显得很糊。

灵敏度分析一般做两个维度：一是EV渗透率从10%到100%变化，观察系统总成本和弃风率的变化趋势；二是惩罚系数( \lambda_w )从低到高变化，观察调度策略从“接受弃风”到“不惜成本消纳风电”的过渡。这两组分析做完，论文的“讨论”章节素材就齐了。

4. Python代码实现：从Matlab迁移的完整思路

4.1 Python复现的核心动机

既然Matlab版本已经能跑通全部模型，为什么还要用Python再写一遍？我个人的动机有两点：一是Python生态完全开源，在Linux服务器上部署不需要买license，而且后续如果要把模型接到机器学习代理模型（用神经网络替代场景抽样）上，Python无缝衔接；二是很多审稿人现在更偏好Python代码，因为审稿和复现的门槛更低。

我用Python复现这套调度模型的体验是：如果你的核心诉求是“快速复现论文结果”，Matlab+Yalmip依然是最快的路径；但如果要做二次开发和扩展，Python的工程友好性明显更强——数据预处理用pandas，场景削减用scikit-learn，求解用pulp或Python-MIP，可视化用matplotlib，一整套流程都能在一个生态里完成。

4.2 建模库选型对比

Python里做MILP建模有几个选择：pulp最简单，适合教学和验证；Python-MIP（mip包）性能和可扩展性更好，内置CBC求解器，也能无缝切换Gurobi；cvxpy适合做凸优化和锥规划，对MILP的支持也比较完善。我实测下来的选型建议是：如果场景规模在10个以内、机组不超过20台，pulp完全够用，代码简洁易读；如果要跑大规模场景或者做两阶段鲁棒优化，用mip包会更稳妥，CBC不行就换Gurobi，接口不变。

pulp的建模方式比Yalmip更“显式”，每个变量都要声明上下界和类型：

import pulp # 创建问题实例 prob = pulp.LpProblem("Cooperative_Dispatch", pulp.LpMinimize) # 决策变量: 火电出力 (连续), 启停 (0/1), EV充放电功率 (连续) P_g = pulp.LpVariable.dicts("P_g", (range(n_gen), range(T)), lowBound=0) u = pulp.LpVariable.dicts("u", (range(n_gen), range(T)), cat="Binary") P_ev = pulp.LpVariable.dicts("P_ev", (range(n_cluster), range(T)), lowBound=-20, upBound=20) # 目标函数（以线性煤耗简化示例） objective = 0 for g in range(n_gen): for t in range(T): objective += cost_coeff[g] * P_g[g][t] + startup_cost[g] * u[g][t] for t in range(T): objective += 50 * (P_ev_sum_forecast[t] - P_ev_actual[t]) # 弃风惩罚等 prob += objective # 约束条件 for t in range(T): prob += pulp.lpSum(P_g[g][t] for g in range(n_gen)) + \ pulp.lpSum(P_w[n][t] for n in range(n_w)) == \ load[t] + pulp.lpSum(P_ev[c][t] for c in range(n_cluster)) for g in range(n_gen): prob += P_g[g][t] <= P_gmax[g] * u[g][t] prob += P_g[g][t] >= P_gmin[g] * u[g][t] # 求解 status = prob.solve() print("目标值: ", pulp.value(prob.objective))

这段代码里有几个容易踩坑的地方。第一，pulp的LpVariable.dicts创建多维变量时，嵌套字典的索引要小心写错——我建议用(g, t)元组作为单一索引，比双重嵌套清晰很多。第二，pulp默认的CBC求解器是开源的，性能比不上Cplex/Gurobi，但多数调度模型规模不大，CBC完全能hold住。第三，约束里用pulp.lpSum而不是Python内置的sum，pulp的lpSum是懒加载的，构建大规模模型时效率更高。

4.3 从Matlab迁移到Python的注意点

迁移过程不是简单的逐行翻译，有几个地方必须重新设计。

第一，Yalmip的binvar和sdpvar在pulp里分别对应cat="Binary"和不指定cat（默认连续）。但Yalmip对变量维度是隐含的（矩阵），pulp是显式的，每个变量都要单独声明。维度多了以后，建议写一个辅助函数自动生成变量字典，避免每行都手写。

第二，Yalmip的约束拼接是[C, C1, C2]，pulp的约束是prob += constraint逐条添加。批量生成约束时，注意pulp的约束表达式要能在add到problem之后再修改——不能在执行prob.solve()之后再改约束，否则求解器会报异常。

第三，Python场景削减可以直接调scikit-learn的KMeans，比Matlab手写K-means方便得多。但有一点要注意：KMeans聚类的输入是二维数组（样本数×特征数），你的场景维度是24小时风速曲线，特征数是24，1000个场景就是1000×24的矩阵。聚类结果打乱顺序，但每个样本的聚类标签能帮你统计每个场景的概率。

第四，如果你需要对比Matlab和Python的求解结果，务必保证两边用的求解器相同。我用Cplex两边各跑一次，偏差在1e-6以内；但如果一边用Gurobi一边用CBC，数值差异可能会到1%以上，这种差异来自求解器内部的算法选择，不是你的代码有bug。沟通实验对比结果时，这一点要提前说明，不然容易被质疑复现有问题。

5. 常见问题与排查技巧实录

5.1 求解器报“infeasible problem”怎么办

这是复现调度模型时最常遇到的错误，而且大多出现在你信心满满以为模型没问题的时候。排查顺序我建议按下面几步走：

先检查功率平衡约束里是否混入了符号错误。利用风电场-负荷-充电功率的数值代入，手算一下各时段是否满足sum(发电) == sum(负荷)。很多情况是EV放电功率符号定义反了——你定义P_ev为负表示放电，但功率平衡约束里写成了加号，结果等式永远不成立。

再检查机组爬坡约束和最小启停约束是否相互矛盾。爬坡约束要求机组出力变化不能超过爬坡率，如果你设定的爬坡率特别小，而负荷曲线又剧烈波动，系统可能真的找不到可行解。最小启停时间约束更是严格：机组的启动动作触发后，必须保持开启至少( T_{up} )小时。如果负荷低谷很短，这个约束会和机组最小出力约束打架。解决方法是把最小启停时间数值调小甚至置零，先验证模型框架能用，再逐步恢复完整约束。

最后检查场景数据是否越界。风电场景里有负风速、光伏场景里有大于1的归一化辐照度，这通常是不小心把概率分布的随机抽样序列用了错误的参数。所有场景数据的取值范围在建模前做一次数值校验，能省掉大量排查时间。

5.2 运行时间过长怎么优化

模型规模一大，求解时间从几十秒涨到几十分钟一点都不奇怪。我常用的优化手段有三个：

第一，给所有0/1变量提供一个好的初始解。用启发式算法（比如先不考虑EV约束的机组组合问题）跑出一个启停方案，通过Yalmip的assign或pulp的setInitialValue给求解器当热启动，Cplex/Gurobi的求解时间能减少50%以上。

第二，削减整数变量的规模。前面提到的EV集群聚合就是一个例子——用集群替代单车，整数变量数量不变，但连续变量大幅减少，整体模型求解速度更快。还可以对机组做“必开/必停”预处理：如果某台机组在当前负荷水平下无论如何都不可能开机，就直接把它从候选机组里剔除。

第三，给求解器设定一个可接受的MIP gap上限，比如1%。对工程应用来说，1%的次优解与最优解的成本差异几乎可以忽略，但求解时间可能缩短一个数量级。实测下来，用ops = sdpsettings('solver', 'cplex', 'mipgap', 0.01)把gap从默认的1e-4放宽到1%，24时段、10机组的模型求解时间从大约200秒降到30秒以内，完美满足快速迭代的需求。

5.3 结果不符合预期的排查路径

求解成功、目标值也正常，但画出来的图看起来不合理——比如EV充电时段完全没避开负荷高峰，或者弃风率比无序充电还高。这种情况最可能的原因是目标函数中相应部分的权重系数设置不当。

EV调度费用系数如果远小于火电煤耗成本，优化器就会放任EV在任何时段充电，因为“不划算”去专门调整充电时段。弃风惩罚系数如果低于EV充电成本，系统就会更倾向于弃风而不是让EV多充电消纳——这从优化角度是“正确”的，但不符合你想要的“协同”效果。所以复现论文时，除了复现公式，更要复现论文中的参数具体数值。如果论文没给全，根据结果反推参数是非常实用的技巧——先设一组合理初值，看结果趋势，再调整参数直到曲线形态和论文图基本一致。

还有一类问题是EV的SOC约束太松或者太紧。SOC约束设在0.1到0.9之间，充电需求设为离网时达到0.8，如果可用充电时间很短，约束会让充电功率在可用时段内冲到上限——这种“强制充电”行为是物理约束导致的，不是bug。但如果SOC上限设成1.0而电池模型又带了充电效率，你会看到SOC在调度末期才勉强达到目标值，中间全程贴着上限跑，这就不太正常了。

5.4 实测问题速查表

最后聊两句我在复现过程中比较深的体会。这个题目说难也难，说简单也简单——核心数学模型并不复杂，难的是把每个环节都做扎实：场景削减是否有理论依据、EV聚合是否物理合理、目标函数的权重设置是否符合论文意图。很多同学复现论文失败，不是模型算不出来，而是结果图跟论文对不上，于是开始怀疑自己的代码有bug，实际上往往是参数设置的问题。我的经验是先复现目标函数的量级和趋势，再逐步逼近论文的数值结果——一步到位复现出完全一样的数字几乎不可能，但复现出图线趋势一致、数量级正确的结果，已经足够支撑论文的对比实验了。如果你正在做这个方向，建议先把最简单的“无EV调度”场景跑通做baseline，再一步步加入EV有序充电和V2G，每一步都有对照，出问题才好定位。