前阵子为了把一篇题为《可再生能源发电与电动汽车的协同调度策略研究》的硕士论文从头到尾复现一遍我把自己关在电脑前整整折腾了两周。这篇论文的主体是一个混合整数线性规划模型用Matlab调用商业求解器求解目标是让风电、光伏和电动汽车的充放电行为在一个调度周期内互相配合从而降低系统运行成本、提高新能源消纳率。如果你正在复现同方向的论文或者想搞明白这类调度模型从数学公式到代码实现之间到底隔了哪些坑这篇文章应该能帮你省下不少时间。我把整个复现过程拆成五个部分来讲模型的边界条件怎么理解、数学公式里的每一项对应什么物理含义、Matlab代码怎么搭骨架、求解器参数怎么调、以及我自己踩过的几个典型问题。文章里出现的所有代码片段都是从可运行工程里抽出来的简化版参数和时段都做了删减但核心逻辑保持一致。1. 从模型假设到代码骨架这篇论文的边界条件1.1 论文到底在解决什么问题先说背景。风电和光伏的出力受天气影响具有明显的波动性和不确定性而传统调度模型通常把机组出力当成平滑可调变量来处理两者一叠加系统很容易出现“该消纳风光时火电压不下来需要顶出力时风光又不够”的局面。这不是单纯增加装机就能解决的问题而是需要在时间尺度上重新安排各类电源和负荷的配合。电动汽车恰好提供一个抓手大多数车辆一天中大部分时间停在车位充电时间可以错峰甚至可以在高峰时段向电网放电。于是“可再生能源发电 电动汽车协同调度”就变成一个很自然的优化命题把EV当作可平移、可反向输出的储能资源与火电、风电、光伏放在同一个优化框架里。复现这种论文最大的价值不是跑出一张漂亮的功率平衡图而是理解作者在建立模型时做了哪些取舍。比如为什么要把EV聚合成一个等效集群为什么不考虑三相潮流为什么使用混合整数线性规划而不是启发式算法。这些取舍决定了结果能回答什么问题也决定了代码能不能在合理时间内求解。我自己的经验是先花两天读模型再动手写代码比一边读一边写要快得多。1.2 系统框架与调度时间尺度这篇论文采用的是一个典型的单节点日前调度框架。系统母线连接三类常规电源火电机组、可再生能源风电场和光伏电站、一个电动汽车集群以及基础负荷。调度周期为24小时时间分辨率1小时共24个时段。所有资源共用一个功率平衡方程不考虑线路潮流和电压约束也不考虑网损。这种设定对硕士论文阶段的研究来说是合理的研究重心在于多资源之间的时间耦合关系而不是配电网的物理分布特性。时间分辨率的选择也需要刻意关注。有些论文会用15分钟间隔模型规模相当于1小时间隔的4倍求解时间膨胀可能不止4倍。如果原论文是1小时步长复现时尽量不要擅自改成15分钟因为预测精度、爬坡约束尺度、EV接入时间窗的语义都会变化。步长一变结果没有可比性。1.3 关键假设与代码映射我复现时把作者在原文中列出的假设整理成了四条每条都直接映射到代码里的一个步骤单节点模型所有电源和负荷在同一个母线上功率平衡用等式表达。EV聚合模型不逐辆跟踪每辆车而是把车队当作一个“大电池”用总容量、总充放电功率、总SOC来建模。风光出力给定风电和光伏的可用功率作为场景参数直接输入模型只能决定消纳多少、弃掉多少。电池寿命损耗线性化把充放电对电池寿命的影响近似成与充放电电量成正比的惩罚成本。这些假设不是凭空来的而是为了让MILP模型在可接受时间内求得全局最优解。如果一开始就加入非线性电池模型、配电网潮流和逐车SOC那问题会变成难以求解的MINLP。复现的第一步就是忠实地把作者假设的模型实现出来之后再谈扩展。2. 目标函数、EV约束与功率平衡模型的拼装逻辑2.1 目标函数成本项的写法与取舍目标函数是整个模型的指挥棒写错一项后面所有结果都会被带偏。我采用的总目标是最小化系统在调度周期内的总成本主要包含六项内容火电机组燃料成本机组启动成本弃风、弃光惩罚成本EV电池损耗成本向主网购电的成本如果存在向主网售电则这一项为负收益。可以用一个简化表达式描述min F Σ_t [ Σ_g (a_g*P_g(t)^2 b_g*P_g(t) c_g) Σ_g SU_g(t) λ_w*(W_avail(t)-W_use(t)) λ_s*(S_avail(t)-S_use(t)) λ_bat*(P_ch(t)P_dis(t)) π_buy(t)*P_buy(t) - π_sell(t)*P_sell(t) ]这里有几个决策点需要反复斟酌。第一个决策点是燃料成本里的二次项。如果保留二次项问题属于MIQPGurobi这类求解器也能处理但论文复现场景下我通常建议先按原论文的写法来。如果原论文用的是分段线性化就把二次项拆成几段线性成本保持MILP结构如果原论文直接保留了二次项那么我调用Gurobi时会把P_g(t)手动预置为非负连续变量让求解器走MIQP路径。两种方式的数值往往有细微差异论文里通常只报告一种你最好在复现说明里写清楚自己用的是哪一种。第二个决策点是惩罚系数的大小。弃风弃光惩罚系数必须显著高于机组单位发电成本否则优化结果会为了省一点调节成本而主动放弃清洁能源这不符合协同调度的初衷。但系数也不能高到离谱否则任何约束冲突都会通过弃电惩罚来硬解导致机组的启停模式失真。我一般先按论文给出的值复现再做一个灵敏度分析观察弃风率随系数变化的趋势确认系数选择在合理区间。2.2 风光出力与弃电惩罚风电和光伏在模型中不是决策变量而是作为可用功率参数出现。真正被决策的是每一时段的实际消纳功率它必须介于0和可用功率之间0 ≤ W_use(t) ≤ W_avail(t) 0 ≤ S_use(t) ≤ S_avail(t)弃风量就是W_avail(t) - W_use(t)。这里有一个容易忽略的细节如果风光出力不是作为参数而是被写成变量并加上“实际等于预测”的约束那模型会变成随机规划求解方式完全不同。所以复现前必须确认论文里是确定性参数还是场景集合。我遇到的论文多数采用多场景随机规划生成若干组风电光伏场景每组场景对应一个概率目标函数是所有场景下的期望成本。如果导出的代码只有一组场景那就相当于只取了期望值场景。想验证不确定性建模是否正确可以先跑单场景再逐步增加场景数量观察成本变化是否收敛。场景数量从10加到50如果成本变化不超过1%说明场景规模已经够用了。2.3 EV聚合模型与SOC递推EV集群的建模是这篇论文最有意思的部分。把车队等效成一个电池后需要描述它的电量变化规律。每时段的SOC递推可以写成Soc(t1) Soc(t) (η_ch * P_ch(t) - P_dis(t)/η_dis) * Δt / Cap其中Cap是整个集群的等效总容量单位是MWh而P_ch和P_dis单位是MWΔt单位是h。很多刚上手的人会忘记除以Cap导致SOC变成一个很大的数值模型仍然觉得“合理”。这属于单位换算错误在结果里往往表现为EV的总电量变化和功率曲线对不上。除了SOC递推还需要考虑三条边界约束。第一SOC上下限通常取0.1到0.9预留安全裕度第二充放电功率不能超过聚合后允许的总功率上限第三一个时段内不能同时充电和放电需要引入二进制变量P_ch(t) ≤ M * y(t) P_dis(t) ≤ M * (1 - y(t))M取值不能太小否则会限制功率上限也不能盲目取特别大以免数值病态。我习惯取总充放电功率上限的两倍刚好满足大M约束的常规要求。电池损耗惩罚项的本质是“让EV别被白嫖”。如果不加这项调度模型会把EV电池当作免费储能可能让它在几个时段内反复满充满放优化结果非常激进。加上一个与充放电电量成正比的损耗成本后EV才会在真正有价值的时候比如高电价时段放电参与调节。这个系数的量纲通常是元/MWh我复现时通过调节它来匹配论文里EV的利用频率。2.4 功率平衡和备用约束的隐性耦合功率平衡约束是硬性的不允许任何偏差Σ_g P_g(t) W_use(t) S_use(t) P_dis(t) P_buy(t) P_load(t) P_ch(t) P_sell(t)如果EV只能从电网充电放电功率为零这个等式就退化成传统电力平衡。加入放电后等式的右边相当于“净负荷 原负荷 充电功率 - 放电功率”这也是很多论文在画图时把EV净充电曲线画成虚线的原因。除了功率平衡还要加系统旋转备用约束因为新能源预测不可能绝对准确。典型写法是正备用容量不小于负荷的5%加风光预测出力的10%负备用类似。备用约束会占用机组可调出力空间也会影响EV的充放电行为——在备用紧张时段EV不能无限制地充电因为充电会推高净负荷消耗更多上备用。这是协同调度中容易忽略的隐性耦合。网架简化的核心意义是把非线性潮流从问题中剥离。如果后续要考虑具体配电网需要把单节点功率平衡替换为DistFlow方程再通过二阶锥松弛转成凸问题。那是另一个话题不是这篇论文的研究重点我建议在基础版本完全跑通后再碰。3. Matlab数据流与Yalmip建模从输入到求解的落地细节3.1 数据目录与预处理脚本好的代码结构比好听的变量名更重要。我在复现时把所有输入数据集中到几个初始化脚本里方便后续替换论文中的不同算例。目录结构大致这样main.m主入口负责读取数据、建模、求解、画图data/load_profile.m负荷曲线data/wind_solar.m风光出力data/ev_config.mEV数量、容量、接入时间窗、初始SOCdata/unit_config.m火电机组成本系数、上下限、爬坡率utils/yalmip_model.m构建约束和目标函数utils/plot_results.m结果可视化。这样做的好处是想调整某个算例只需要修改一个数据脚本不需要在模型代码里到处找参数。我第一次复现时把所有参数写在一个大脚本里结果为了对比三组数据复制出三个版本改了一个参数忘了同步最后对不上结果白忙了一整天。后来才痛定思痛把数据和模型彻底分开。一个小细节是固定随机数种子。如果EV初始SOC、场景生成或者负荷曲线包含随机因素在脚本开头用rng(42)固定种子否则每次运行结果都会有细微差异论文里没法交代。3.2 Yalmip变量与约束的写法我习惯用Yalmip作为建模层它的好处是语法简洁底层可以自由切换求解器。以一个三台火电机组、24时段的小算例为例变量声明可以这么写T 24; n_gen 3; Pg sdpvar(n_gen, T, full); % 机组出力 u binvar(n_gen, T, full); % 机组启停状态 v binvar(n_gen, T, full); % 机组启动指示 P_ch sdpvar(1, T, full); % EV集群充电功率 P_dis sdpvar(1, T, full); % EV集群放电功率 Soc sdpvar(1, T1, full); % EV集群SOC y binvar(1, T, full); % 充放电互斥标志注意Soc的维度是T1对应0到24共25个点。凡是涉及Soc(t1)和Soc(t)的约束循环范围应当是t1:T。这个维度错位是新手最容易踩的坑后面我会专门展开。约束组装我倾向于用循环因为可读性好排查问题方便。比如功率平衡Constraints []; for t 1:T Constraints [Constraints, ... sum(Pg(:,t)) W_use(t) S_use(t) P_dis(t) ... P_load(t) P_ch(t)]; end爬坡约束同样写进循环注意首时段没有上一时段循环从2开始Ramp_up 30; Ramp_down 30; for t 2:T Constraints [Constraints, Pg(:,t) - Pg(:,t-1) Ramp_up]; Constraints [Constraints, Pg(:,t-1) - Pg(:,t) Ramp_down]; end目标函数在Yalmip中直接写成可加和的表达式Objective sum(sum(a .* Pg.^2 b .* Pg c)) ... sum(SU_cost .* v) ... wind_pen * sum(W_avail - W_use) ... solar_pen * sum(S_avail - S_use) ... bat_pen * sum(P_ch P_dis) ... sum(price .* (P_buy P_sell));这里有一个非常重要的细节a .* Pg.^2会让问题变成MIQP。如果论文明确要对燃料成本做分段线性化就不能这样写。分段线性化的标准做法是用一组辅助变量表示每段功率然后加约束% 假设每台机组分为3段 P_seg sdpvar(n_gen, T, 3, full); seg_lim [0.4, 0.3, 0.3]; % 每段占最大出力比例 seg_cost [15, 18, 22]; % 每段边际成本然后保证P_seg的每一段非负且不超过上限总出力等于各段之和。这样得到的目标函数是线性的模型就是MILP。我建议在代码注释里写清楚自己复现的是MIQP还是MILP版本因为求解器日志中的objective值不同容易造成“对不上论文”的错觉。3.3 求解器选择和调参求解器是MILP的引擎。我测试过的几个选择如下求解器授权方式求解速度适用场景Gurobi学术免费/商业付费非常快中大规模MILP/MIQPCplex学术免费/商业付费非常快同上习惯用哪个都行CBC开源免费较慢验证模型用大规模不太建议MATLAB intlinprog随Matlab中等不想装外部求解器时备选在Yalmip里切换求解器非常方便只需要在sdpsettings里指定名称。以Gurobi为例ops sdpsettings(solver, gurobi); ops.verbose 1; ops.gurobi.MIPGap 1e-4; ops.gurobi.TimeLimit 300; sol optimize(Constraints, Objective, ops);MIPGap是求解器判定最优的阈值设成1e-4意味着允许与最优值有0.01%的偏差。对调度问题来说这个精度足够能大幅缩短求解时间。TimeLimit是硬保护防止模型特别难解时一直占用机器。如果sol.problem返回0表示求解成功返回1通常表示不可行返回2表示数值问题。我建议拿到结果后马上判断这个返回值再继续。调参经验是先不加EV用纯火电风光跑通再加EV充电最后加EV放电和电池损耗惩罚。每一步都对比求解器日志里的约束数量和变量数量确认模型规模在预期范围内。数据量特别大时可以用ops.cachesol1让Yalmip缓存上次解方便热启动。3.4 结果保存与可视化求解完成后的第一件事不是画图而是把结果保存下来result.Pg value(Pg); result.P_ch value(P_ch); result.P_dis value(P_dis); result.Soc value(Soc); result.Objective value(Objective); save(result.mat, result);然后再画图。我最常用的是一张堆叠面积图表示功率平衡figure; area([result.Pg, W_use, S_use, result.P_dis]); hold on; plot(P_load result.P_ch, k-, LineWidth, 2); xlabel(时段/h); ylabel(功率/MW); legend(火电,风电,光伏,EV放电,净负荷);还有一张EV的SOC时序曲线figure; plot(0:T, result.Soc, LineWidth, 2); hold on; yline(0.1, r--); yline(0.9, r--); xlabel(时段/h); ylabel(SOC);这两张图基本能回答所有“结果对不对”的直觉问题SOC是否在上下限内、功率平衡曲线是否闭合、EV在高电价时段是否有放电尖峰。如果一眼看去不满足常识后面再做定量验证也不迟。4. 不可行、求解爆炸与数据错位三类典型坑4.1 模型不可行的完整定位链路我在第一次跑通Yalmip模型时optimize返回的信息是infeasible。这很打击人但也是绝大多数新手会遇到的第一道坎。我整理出完整的排查链路。第一步把目标函数改成常数0Objective 0; ops.verbose 2; sol optimize(Constraints, Objective, ops);如果还是不可行说明可行域本身就空问题不在目标函数。第二步从约束里去除EV相关约束。具体做法是在循环里只保留火电和风光看能不能跑出可行解。如果能问题出在EV约束如果还是不可行问题在火电爬坡或启停逻辑。第三步定位到EV约束后优先检查Soc的索引。Yalmip的check函数可以直接检查每条约束是否满足residuals check(Constraints);residuals大于一个阈值表示约束越限负值或接近于0表示满足。配合assign给变量赋一组可行初值你就能秒看出是哪一行约束把模型搞矛盾的。这种问题的根子往往是“维度对不上”。SOC有T1个点而充放电功率和负荷只有T个点写循环时差一个下标整个可行域就全乱了。所以我在代码里统一了习惯所有和SOC递推有关的循环变量都从1到T右侧的SOC索引写t1左侧写t凡是用0到T初始SOC的图形单独保存。4.2 EV规模扩大后求解时间爆炸论文里如果EV数量是几百辆逐车建模可能还能跑一旦变成几千辆或者几万辆0-1变量数量会完全失控。我做过一个对比实验同样的24时段问题100辆车逐车建模时0-1变量约2400个5000辆车就是12万个求解时间从十几秒涨到十几分钟再加充放电互斥和机组启停直接几小时不出结果。解决办法是采用聚合模型。聚合成一个集群后不管车有多少辆变量数量固定是24个时段乘以几个决策变量。具体实现时把总容量设为所有车辆电池容量之和把最大充放电功率设为所有车辆可同时充放电时功率之和。如果车辆接入时间窗不一样可以按接入时段分成两到三个集群每个集群有独立的SOC和充放电功率但变量数量仍然有限。这个选择背后的理由是论文研究的是系统级调度策略不是车与车之间的电量均衡。EV用户个体的差异通过聚合参数来反映比如平均初始SOC、平均最大充放电功率、可调度时段概率分布。如果你非要逐车建模那就要考虑用递推聚合算法或者分布式求解复杂度会高很多。4.3 风光数据单位不匹配我复现时发现一个奇怪现象程序跑出来的弃风率特别高高到不像是“协同调度”的结果。查了很久发现原因不在模型而在数据。我用的风电数据文件里存的是标幺值就是按额定容量归一化后的0到1而我读进代码后没有乘以额定容量就直接当成MW用了等于把300MW的风电场当成了最多1MW的电源。这当然不可能消纳多少风电。排查这类问题我会先做一个“数据报表”disp([风电总可用电量(MWh): , num2str(sum(W_avail))]); disp([光伏总可用电量(MWh): , num2str(sum(S_avail))]); disp([负荷总电量(MWh): , num2str(sum(P_load))]);如果风电总可用电量只有几十MWh而负荷是几千MWh那肯定有问题。还要检查时间粒度和单位。15分钟数据重采样到1小时功率用15分钟的平均值而不是累加值如果直接累加会把功率当成能量结果翻4倍。另一方面如果论文用的是“典型日”曲线可能已经把数据做了平滑处理不适合直接比较绝对值。处理办法是写一个独立的数据预处理脚本输入原始数据输出统一格式的struct并在脚本里加assert判断关键量级assert(max(W_avail) WindCap * 1.05); assert(all(P_load 0));这样能在跑模型之前就发现数据异常而不是等结果出来再猜。4.4 结果对不上论文时的排查顺序当我终于跑出一个“看起来合理”的解却发现和论文里的成本数值、弃风率对不上这时候要做的不是立刻调参数而是按下面这套顺序排查。先做极端测试把风光出力全设为0EV也不参与放电目标函数应该等于纯火电满足负荷的燃料成本。这一步可以用手算大概验证比如负荷电量已知、单位煤耗成本已知算出总成本量级。如果这一步都对不上问题一定在机组成本系数、单位换算或者负荷数据。再做功率平衡自检求解完成后把所有功率信号代回平衡方程求最大误差balance sum(result.Pg, 1) W_use S_use result.P_dis ... - P_load - result.P_ch; assert(max(abs(balance)) 1e-4);如果这个断言通过至少说明模型内部自洽。然后比较趋势而不要急着比数值论文里EV应该表现出“夜间充电、白天放电”的模式你的结果如果也这样参数差异只是量级问题。最后再依据灵敏度分析表逐一核验系数。偏差现象常见原因检查思路弃风率明显偏高弃风惩罚太小或风电单位错误增大惩罚系数或核对标幺值EV几乎没有充放电动作电池损耗惩罚过大降低损耗系数成本与论文差一个数量级燃料成本单位或负荷单位错误核对元/MWh与MW/MWh关系启停变量异常频繁启动成本过低或爬坡约束过松调高启动成本检查爬坡率5. 验证模型正确性后还能往哪里扩展5.1 我给自己定的几条硬性检查规则模型正确性比结果“好看”重要得多。我给自己定了几条硬性检查规则。第一每个约束必须能用常识解释出力最大不超过上限爬坡速度不超过额定值SOC在0到1之间功率平衡处处成立。第二求解完成后把Soc(t)和充放电功率画到同一张图上观察电量变化是否等于充放电功率积分。第三做简单的灵敏度分析负荷曲线整体上调10%总成本应当明显上升电价低谷时段EV充电量应当增加。只要有一条不符合就说明模型里还有没找到的逻辑问题。check(Constraints)返回值是很好的辅助工具但它只能告诉我们“约束是否满足”不能告诉我们“约束是否合理”。所以我还经常把最优解的变量值打印出来逐时段对着功率平衡表用计算器加一遍。虽然累但能发现很多隐蔽的单位错误。5.2 从基础模型到扩展模型的方向基础版本稳定后如果还想做点自己的东西我建议按这样的优先级扩展。第一把EV从单节点扩展到多节点配电网用DistFlow方程替代单一功率平衡但要注意模型从线性变成非线性需要做二阶锥松弛。第二考虑EV用户响应不确定性和充电桩接入的不确定性用机会约束或分布鲁棒优化替代确定性场景。第三把集中式调度改成“调度中心-聚合商”分层结构用交替方向乘子法之类的分布式算法求解这样能回答“EV聚合商是否愿意透露内部参数”这类问题。第四把日前调度和实时滚动优化结合起来把MILP解出来的结果作为MPC的参考轨迹。需要注意每次扩展都会带来新的求解复杂度不要一上来就都加进去。我先用单节点的基础模型把算法框架调稳再逐步增加模块每增加一个模块就做一次对照实验记录对目标函数、求解时间、解质量的影响。最后再分享一个自己的小习惯每次跑完模型我都把result.mat、求解器日志和数据脚本版本一起存档文件名里带上日期和参数标记。改参数后如果发现回不到原来的结果能快速定位是哪一步覆盖了数据。复现论文本质上是把别人的想法翻译成自己的代码耐心比聪明重要得多尤其是当你盯着infeasible发呆的时候先按思路拆解再动手改代码往往比反复重启Matlab更有效。
