做电力系统调度的同学尤其是研究可再生能源和电动汽车入网方向的人大概率绕不开“协同调度”这四个字。我前阵子为了硕士论文复现把“可再生能源发电与电动汽车的协同调度策略”完整用Matlab走了一遍从场景生成、约束搭建到调用求解器踩了不少坑。这篇东西不是理论推导的复读而是告诉你实际建模时怎么定义变量、怎么列约束、怎么让求解器顺利跑出结果以及哪些错误是你大概率也会遇到的。这篇内容的适用对象很明确正在准备硕士论文里仿真章节的研究生、做微电网/配电网调度仿真的工程师还有那些装了Matlab但不知道Yalmip到底怎么用的人。如果你只是想要一份能直接跑的代码那这篇文章同样能帮你搞清楚代码每一步在做什么免得拿到开源代码后连报错都看不懂。1. 协同调度问题拆解定位核心模型1.1 问题边界与决策变量先说清楚模型边界。绝大多数硕士论文里的协同调度不会把电网全模型都搬进来而是把研究对象简化成一个微电网或者一个等效的配电网节点。我做的是微电网场景里面有风电机组、光伏阵列、传统燃煤机组、储能系统还有一定数量的电动汽车。电动汽车既可以作为普通负荷充电也可以在电价合适的时候向电网放电V2G模式。这个模型要回答的问题是在一天24小时或96个时段内每台机组出多少电、储能什么时候充什么时候放、电动汽车什么时候充电、允许多少车在某个时段放电才能让总成本最低同时把可再生能源尽量用掉。决策变量分为连续变量和二进制变量。连续变量包括常规机组出力、储能充放电功率、电动汽车充放电功率、向主网购电和售电功率。二进制变量则用于表示机组启停状态、电动汽车充放电状态。这里有一个容易忽略的地方电动汽车不能同时充电又放电所以至少需要两个二进制变量来分别代表充电状态和放电状态并且要加上互斥约束u_ch u_dis 1。很多初学者把模型定义得太复杂一上来就考虑每辆车的行程、用户行为、充电桩位置结果变量数量爆炸求解器根本算不动。硕士论文复现阶段我建议先把所有电动汽车聚合成一个等效储能池比如100辆车每辆车40kWh那么聚合容量就是4000kWh聚合最大充电功率就是700kW。这样做虽然忽略了个体差异但能先把调度逻辑跑通后续再扩展精细化模型。1.2 目标函数与成本建模目标函数看起来简单写起来容易乱因为成本项实在太多。常见的目标是最小化一天内的系统总运行成本可以拆成五个部分常规机组燃料成本、可再生能源发电的运维成本通常很小、弃风弃光的惩罚成本、储能和电动汽车的电池损耗成本以及从主网购电的成本减去向主网售电的收益。常规机组成本通常写成输出功率的二次函数C_G(P) a*P^2 b*P c。但二次函数会让模型变成MIQP求解器求解速度明显变慢。复现时我一般做分段线性化把功率范围切成几段每段用一个线性斜率表示这样一来模型变成MILPGurobi和Cplex跑起来快得多。如果不做线性化遇到大场景数求解时间可能从几分钟变成一小时。弃风弃光惩罚成本很关键它决定模型会优先消纳可再生能源。惩罚系数一般取5001000元/MWh远高于发电成本这样模型宁可从电网买电也不会轻易弃风。但注意惩罚系数不要设成1e6这种量级否则数值病态会让求解器报出各种奇怪的可行性问题。电动汽车V2G放电时电池会额外老化所以目标函数里要加放电惩罚项。基础做法是给放电功率乘一个单位损耗成本比如0.2元/kWh。如果省掉这一项模型会为了赚峰谷电价差疯狂放电产生不现实的调度结果。1.3 约束条件体系约束是调度模型的骨架少一条约束结果就可能是空中楼阁。最核心的是功率平衡约束所有电源常规机组、风电、光伏、储能放电、EV放电、购电的功率之和必须等于所有负荷基础负荷、储能充电、EV充电、售电的功率之和。这条约束在每个时段都要满足写成向量形式时别忘了用sum(x,1)来确保按时间维度求和。储能约束包括SOC状态转移方程、充放电功率上下限、SOC上下限以及一个很重要的“最终SOC等于初始SOC”约束否则模型可以把蓄电池的能量偷走造成成本偏低。SOC计算时注意单位一致性功率单位是kW储能容量单位是kWh时间步长如果是1小时那么SOC变化量就是功率除以容量如果时间步长是15分钟必须乘上0.25否则电量会差4倍。这个坑我见过不止一次。电动汽车约束比储能多一层“充电需求约束”。每辆车每天至少需要充入一定的电量否则用户没电用。比如聚合体一天需要充入2000kWh那么所有时段充电总和减去放电总和后的净充电量必须大于等于2000kWh。另外EV的SOC也要限定在10%到90%之间太深放电对电池不友好。最后是常规机组爬坡约束相邻两个时段的出力变化不能超过上限比如30MW/h。如果这一步漏了调度结果可能是机组出力剧烈波动看起来最优但实际没法执行。爬坡约束本身很简单就是abs(P_G(t)-P_G(t-1)) Ramp用Yalmip时写成P_G(:,2:T)-P_G(:,1:T-1) Ramp和P_G(:,1:T-1)-P_G(:,2:T) Ramp两条。2. 可再生能源不确定性处理从场景生成到鲁棒/随机优化2.1 场景生成方法光伏和风电出力都有随机性如果直接用期望值代替调度结果会很“天真”不考虑实际波动。处理随机性最常用的方法是场景法生成大量代表可能出力曲线的样本然后在样本上求期望成本。场景生成不是随便造数据。风电出力一般用威布尔分布描述风速再通过风功率曲线换算成出力光伏出力则用Beta分布描述光照强度。也可以用历史出力数据直接采样比如NREL或国内电网公开的出力数据集。我在复现时用Matlab自带的randn配合正态扰动生成200个初始场景每一步都加入时序相关性避免相邻时段出力突变。一个容易被忽略的点场景不仅包含可再生能源出力还应该包含负荷的波动尤其是电动汽车接入后的充电负荷。如果只考虑风电随机而忽略负荷随机论文审稿人一眼就看出来了。博士生不会太在意但硕士论文还是要把基础的不确定性框架搭完整。生成场景时要注意样本数量。500个场景做仿真求解时间不可接受20个场景则可能丢失尾部风险。我建议先做400个削减到15个这样兼顾速度和精度。2.2 场景削减算法场景削减的核心思想是用少数具有代表性的场景替代原场景集且保持概率分布特性。常见的算法有三种K-means聚类、快速前向选择FFS、同步回代消除SBR。K-means是最直观的把400个场景聚成15类取每个类的中心作为代表场景概率为该类样本比例。但K-means对初始聚类中心敏感聚类结果可能局部最优。FFS和SBR是经典的启发式场景削减基于场景两两间的距离通常用欧氏距离逐步合并。SBR的做法是每次找一对总概率最小的场景删除其中一个把它被删的概率加到另一个上直到目标数量。这个方法能保留原场景集的分布形状也比较容易手写实现。我用Matlab手写了一个SBR核心就是计算距离矩阵、更新概率。代码量大约60行比调用复杂工具箱更可控。削减后要验证计算削减前后场景集的均值和方差看差异是否在可接受范围内。差异太大说明目标场景数太少。下表是我常用的出力和负荷场景对比指标原始400场景削减后15场景误差风电平均出力/MW12.412.83.2%光伏平均出力/MW8.17.92.5%负荷峰值/MW21.721.50.9%风电标准差3.22.99.4%标准差误差稍高一点但对调度结果影响不大因为优化的是期望成本均值误差更关键。2.3 不确定性建模选择处理不确定性的范式有随机规划和鲁棒优化两种。随机规划相当于把多个场景都放进约束里目标函数是各场景成本的期望。实现方式就是在Yalmip里对每个场景分别建约束目标函数是sum(p_s*Objective_s)。这种方式简单直接求解器能处理但需要保证每个场景下的约束都可行经常会遇到某些极端场景导致问题不可行。鲁棒优化则是找一个最坏情况下的最优决策。典型表示为min max_s。解决鲁棒优化通常需要引入对偶变量把内层max问题转化为约束条件数学推导繁琐但所得方案在面对不确定时有更强的保障。硕士论文复现时如果原标题只写了“协同调度”大概率是随机优化没必要主动把难度升到鲁棒。我当初直接选随机优化把精力放在如何设计约束保证所有场景可行比如引入可调变量如负荷削减来处理极端场景下的失负荷。3. Matlab代码实现要点从建模到求解3.1 工具箱与求解器选择Matlab里做优化调度最推荐的方案是Yalmip商用求解器。Yalmip是一个建模工具箱它把模型翻译成求解器能理解的内部形式解完之后还能用value()把变量结果取出来。它支持Cplex、Gurobi、Mosek等。Gurobi学术版可以免费申请license速度非常快尤其求解MILP时比Matlab自带的intlinprog强一个量级。Cplex现在对学术用户也开放免费下载。如果不想申请也可以用免费的SCIP或CBC求解器但大规模问题会慢很多。安装过程常出问题很多同学把Yalmip下载后塞到当前路径但忘了把求解器路径加进Matlab。正确的做法是在startup.m里写入如下命令addpath(genpath(D:\yalmip)); addpath(genpath(D:\gurobi\win64\matlab)); gurobi_setup; savepath;然后运行yalmiptest看到所有求解器都显示OK才算装好。如果Gurobi后装的记得重新savepath不然Matlab重启后还得重新添加。求解器选择上如果模型是纯线性规划没有0-1变量可以用Gurobi的LP如果有0-1变量Gurobi用MILP算法默认会开并发割平面。设置选项时我习惯这样写solverOptions sdpsettings(solver,gurobi,verbose,2,gurobi.TimeLimit,300);设置时间限制很重要否则求解器可能陷入一个超大规模的MILP里跑几小时。我自己复现时遇到过一次半天都没跑完最后发现是二进制变量定义多了几百个把-SOC互斥约束拆成每个EV每个时段两个变量100辆车96时段就19200个二进制变量当然慢。后来聚合建模后变量降到几百个几秒就出结果。3.2 代码框架设计与数据结构不要把所有变量塞进一个大脚本。我用四个文件组织代码main.m负责总流程和调用data_define.m里放系统数据model_define.m里放Yalmip建模solve_and_plot.m负责求解和图形输出。如果后续要加多场景再加一个scenario_generate.m。变量定义的核心是用sdpvar声明连续变量、用binvar声明0-1变量。比如N_gen 3; T 24; P_G sdpvar(N_gen, T, full); % 常规机组出力 P_ev_ch sdpvar(N_ev, T, full); % EV充电功率 P_ev_dis sdpvar(N_ev, T, full); % EV放电功率 u_ch binvar(N_ev, T, full); % EV充电状态 u_dis binvar(N_ev, T, full); % EV放电状态这里full参数特别重要。如果不加Yalmip默认认为变量是个n*n方阵如果你传入的是1x24向量它会报维度错误。很多新手在这里卡了一晚上。时间轴统一为行向量功率矩阵统一为变量数 x T所有sum()操作明确指定维度sum(x,1)这样约束构建时不用反复转置。3.3 核心约束构建与目标函数先写功率平衡约束。这里涉及一个常见问题可再生能源出力在场景化之后是一个概率矩阵如果你用期望值那就退化成了确定性模型。我用削减后的场景集每个场景单独建约束所以模型里会有一系列带场景下标(s)的变量。但在确定性验证时可以先不考虑场景只用一个代表场景。model_define.m中的示意如下Constraints []; % 功率平衡约束 % sum(P_G,1)P_windP_pvP_St_dissum(P_ev_dis,1)P_buy ... % P_base P_St_ch sum(P_ev_ch,1) P_sell Constraints [Constraints, ... sum(P_G,1)P_wind_sceP_pv_sceP_St_dissum(P_ev_dis,1)P_buy ... P_base P_St_ch sum(P_ev_ch,1) P_sell];这里P_wind_sce、P_pv_sce是给定场景下的出力序列。如果允许弃风和弃光就把它们替换为P_wind_used并新增弃风变量P_wind_curtail约束为P_wind_used P_wind_curtail P_wind_sce目标函数中加入惩罚成本。EV充放电互斥约束Constraints [Constraints, u_ch u_dis 1]; Constraints [Constraints, P_ev_ch P_ev_ch_max * u_ch]; Constraints [Constraints, P_ev_dis P_ev_dis_max * u_dis];SOC约束要写成递推式% SOC_ev(t1) SOC_ev(t) eta_ch*P_ev_ch(t)/Cap_ev - P_ev_dis(t)/(eta_dis*Cap_ev) - E_drive(t)/Cap_ev for t 1:T-1 Constraints [Constraints, SOC_ev(t1) SOC_ev(t) ... eta_ch*P_ev_ch(:,t)/Cap_ev ... - P_ev_dis(:,t)/(eta_dis*Cap_ev) ... - E_drive(:,t)/Cap_ev]; end注意E_drive代表该时段电动汽车出行消耗的能量如果没考虑出行可以设为零。目标函数建议用分段线性机组成本的增量形式来实现。简单起见复现时可以直接用二次函数让Gurobi当作MIQP求解但一旦场景数超过100MIQP会慢到让人崩溃。所以我用增量成本分段线性化定义机组出力在第k段的变量P_seg(i,k,t)成本就是各段斜率之和。这样目标函数保持线性。3.4 求解与结果输出求解调用一行代码result optimize(Constraints, Objective, sdpsettings(solver,gurobi,verbose,2));求解完之后第一件事是检查状态if result.problem 0 disp(求解成功); else yalmiperror(result.problem) end不要只依赖result.problem为0有时候求解器返回1警告但依然可以接受。yalmiperror能给出比较明确的错误描述。结果提取用value()函数。例如P_G_opt value(P_G); P_ev_ch_opt value(P_ev_ch); P_ev_dis_opt value(P_ev_dis); SOC_ev_opt value(SOC_ev);画图时要注意时序对齐。MATLAB的stairs更适合画电力系统调度图因为功率在时段内保持恒定。使用subplot同时展示机组出力、EV充放电功率、储能SOC和弃风弃光量让调度策略一目了然。有经验的读者都会再做一个“可行性检查”把所有决策变量的最值打印出来看是否在合理区间内。比如min(P_G_opt(:))如果出现负值那一定是有约束写错了。4. 复现过程中常见问题与排查技巧4.1 求解失败与收敛慢最常遇到的情况是求解器提示“Infeasible”也就是约束之间互相矛盾。第一步检查可行性把约束拆分逐个添加进模型测试。比如先只加功率平衡看看能不能解再加机组上下限直到定位到哪条约束导致不可行。还要注意数值尺度。通常功率单位用kW电价单位用元/kWhSOC是小数。如果目标函数里某一部分是0.001量级另一部分是1e8量级求解器数值精度会出问题。解决方法是统一量纲或者把大数值项除以一个基准功率。求解速度慢时先看二进制变量数量。如果二进制变量上千MILP求解时间指数增长。另一个常用技巧是设置“mipgap”也就是最优间隙容忍度。调度问题允许1%的间隙结果差距不大但速度能快很多ops sdpsettings(solver,gurobi,gurobi.MIPGap,0.01);4.2 维度不匹配与变量定义错误Yalmip变量尺寸不一致是最容易报错的地方。比如你定义了P_G sdpvar(N_gen, T)这在Yalmip里默认是一个方阵不对sdpvar(N_gen, T)创建的是N_gen行T列的矩阵没有full也一样。但如果你传入了1作为第一个维度而T是24那么P_G是行向量后续约束里你可能以为它是列向量写P_G(:,1)会报错。最好的习惯是在变量定义前用assert(size(P_G,2)T)来校验。还有个常见坑是binvar默认是不带引用限制的。在不该用二进制的地方用了binvar就会把问题变成超大的MILP。比如表征“电动汽车只能在充电桩接入时充电”不需要用二进制变量可以用时间窗直接强制充电功率为0Constraints [Constraints, P_ev_ch(:, outside_window) 0];。4.3 场景削减后的结果偏差用场景法求期望成本结果往往比实际期望值偏乐观因为削减后的场景覆盖不到极端情况。解决办法是做回代检验把求出的最优解固定住代入原始的400个场景计算真实的期望成本和模型里的期望成本对比。如果偏差超过5%需要增加保留场景数或者换成更好的削减算法。另外有些论文里不是单独做场景削减而是直接用蒙特卡洛在约束里循环所有场景这种做法在小规模系统里可行但场景数一大就完蛋。如果坚持这样做至少要限制场景数在20以内。4.4 电动汽车参数设置的坑EV参数设置直接影响调度结果。常见问题聚合容量过大导致EV放电能力比电网功率还大模型会安排大量V2G把储能变成主力电源。实际中EV不可能全天候联网所以要设定接入时段比如只允许18:00到次日8:00充放电。SOC的初值和终值必须一致否则模型会从初始SOC里“偷能量”。如果设置终值比初值低那么模型相当于免费用掉了电池里的电成本会偏低。我建议加约束SOC_ev(1) SOC_ev(T1)强制一个运行周期内电量守恒。充电需求每天最少多需要仔细定义。如果设定所有EV一天要充入2000kWh而模型为了多放电把SOC压到很低然后集中在下半夜充电虽然满足约束但可能不现实。更稳妥的做法是设定每个EV的入网时间和离网时间离网时SOC要达到某个下限比如90%。但这样就引入了跟具体时间相关的复杂约束我建议在聚合模型里简化为总净充电量约束且最大放电功率不超过总负荷的30%确保V2G只是辅助手段。5. 扩展方向与个人踩坑记录5.1 从单节点到配电网潮流如果论文需要研究线路约束和电压分布那就不能只在单节点上做功率平衡。最简单的是用DistFlow线性化LDF或二阶锥松弛SOCP。DistFlow把功率、电压、损耗关系写成线性或近似线性的等式适合有载调压的设备SOCP则把潮流约束写成锥约束Yalmip能用cone或norm实现Gurobi支持二阶锥。我在复现时最初没有网络约束后来加上了33节点配电网约束数量增加但性能还好。但需要提醒加入网络约束后原本在单节点下可行的EV充放电策略在潮流模型下可能不可行尤其是节点电压越限。需要把EV充放电功率和网络边界联系起来。这一步如果做出来论文的含金量会上一个档次。5.2 从单目标到多目标扩展硕士论文常喜欢写“经济性环保性”协同优化也就是多目标优化。最简单的是加权求和法和ε-约束法。加权求和法把碳排放和运行成本加权成一个目标函数但权重选取很主观。ε-约束法则是把碳排放作为额外约束每取一个排放上限就求解一次得到帕累托前沿。在Matlab里用Yalmip配合循环就能轻松实现ε-约束法。我复现时先用加权法发现碳排放权重稍大调度就恨不得全部用EV放电成本剧增后来改用ε-约束法结果曲线更平滑。如果你也要写多目标建议用ε-约束法并用帕累托图展示。5.3 我的真实复现体验这个项目我前后花了约三周每天两小时。第一周卡在场景削减和变量维度上第二周模型跑通但结果不合常理比如EV疯狂放电导致SOC一直为零第三天突然发现SOC赋初值时用了0.8约束却限制SOC下限是0.9当然不可行。找到后真想拍桌子。个人最大的体会是一定要先从确定性模型开始把所有约束用单场景跑通再引入随机场景。不要一开始就上Monte Carlo加场景削减否则报错时你根本不知道是场景出了问题还是约束出了问题。我调试顺序是无EV、无储能、有储能、有EV一层层加。这样加一个部分就验证一个部分效率最高。最后再分享一个小技巧所有功率和电量的单位用kW和kWh所有成本单位用元电价用元/kWhSOC用0到1的小数。写代码前先把单位写在注释里因为当你盯着满屏数字找错误时单位不一致是最隐蔽的元凶。跑通了之后记得保存一份“运行正常”的备份再改下一步这是我和我身边人用血泪换来的经验。