做虚拟电厂调度的朋友应该都有同一个感受所谓“最优调度方案”本质上是在跟预测较劲。风电场中午报的出力曲线还是平缓的结果上午十点风向一变实际出力掉了一截光伏更是明显昨天看今天的曲线还是平滑钟形当天飘过来一片云直接凹进去一个坑。这时候你昨天辛辛苦苦算出来的那个24小时全局最优计划放到今天真实运行环境里可能就不是最优了甚至严格来说根本没法执行。这个问题在行业里已经有了标准应对手段——把调度拆成日前调度和日内调度两个时间尺度来做。这也是虚拟电厂多时间尺度调度优化研究的核心思路很多顶级SCI论文的模型框架就是这么搭的。这篇文章我把这块内容掰开揉碎讲清楚为什么非得分两个尺度、两个尺度各自的数学模型怎么建、在Matlab里怎么用YalmipCplex把模型跑起来以及复现过程中最常见的坑在哪。适合正在做虚拟电厂、微电网优化调度方向研究或者准备复现相关SCI论文做毕设/横向项目的同学参考。1. 为什么要把调度拆成“日前日内”两个时间尺度预测误差是根源1.1 虚拟电厂调度的底层逻辑在不确定性中做决策虚拟电厂本身不是一座物理电厂而是把分布式光伏、风电、储能、小型燃气轮机、可控负荷这些分散资源通过聚合管控的方式对外统一表现为一个可调的“发电单元”。这个“聚合”的难点在于里面既有可控性很强的微燃机、储能又有完全看天吃饭的光伏和风电。调度的任务就是在每个决策周期内决定各个设备产出多少功率、储能是充还是放、要不要向电网买电卖电使得整个虚拟电厂在满足负荷需求的前提下经济性最好同时满足各类设备自身的运行约束。问题是所有这些决策必须依赖对未来一段时间风光出力、负荷大小、市场电价的预测。而预测这个东西天然是时间尺度越长、误差越大。1.2 预测精度与时域的必然关系电力系统里预测有个经验规律时间尺度越长预测不确定性的边界越大。日前预测提前24小时的分布式光伏出力均方根误差大概在20%—30%的水平超短期预测未来4小时以内的误差通常能压到5%—10%。风功率预测基本也是这个趋势只是幅度更大。预测类型提前时间典型误差水平用途日前预测24小时20%—30%日前计划、机组组合日内超短期预测4小时以内5%—10%滚动修正、实时跟踪这就产生了一个矛盾如果只做日前调度计划做得很“全局最优”但执行起来因为预测误差大概率走样如果只做日内调度虽然预测准了但那些响应慢、启停成本高的设备比如微燃机、需要预先安排的可中断负荷根本来不及反应。所以行业共识是分层协调日前调度负责“定组合”日内调度负责“调出力”。1.3 两个尺度各管什么定组合与调出力的分工逻辑日前调度提前一天运行以1小时或15分钟为时段粒度覆盖未来24小时。它的核心任务是确定机组的启停状态、储能是否参与充放电调峰、以及各时段与大电网的购售电计划。这些问题包含整数决策变量属于典型的**混合整数线性规划MILP**问题。日内调度当天运行以15分钟或5分钟为粒度预测精度更高、系统状态信息更全。它不再重新决定“哪台机组开、哪台机组停”而是在日前给出的启停状态和计划轨迹基础上用滚动优化的方式调整各设备的出力水平、储能充放电功率和购售电计划。这个世界上已经没有了整数变量前提是日内不调整启停求解速度快很多属于**线性规划LP**问题。用生活类比解释就是你计划一个月后去外地参加一个论坛提前订好机票和酒店这是日前计划当天早上起来查实时路况决定走哪条高速、要不要中途充电这是日内修正。不是非此即彼而是配合着来。2. 日前调度混合整数线性规划建模的完整拆解2.1 决策变量整数变量与连续变量的分工先搭一个典型的虚拟电厂骨架。假设里面包含2台微型燃气轮机MT1、MT2单机容量200kW1个风电场WT额定容量200kW1个光伏电站PV额定容量150kW1个储能系统ESS容量600kWh最大充放电功率150kW1个可中断负荷IL最大可削减量50kW与上级电网的联络线购售电功率上限均为300kW日前调度的决策变量可以分成两类第一类是整数变量典型代表是微燃机的启停状态u(i,t)1表示运行、0表示停机。为什么必须要整数因为启停是一个“是/否”的离散决策——你不能开出0.7台燃气轮机来。这个整数的引入直接让日前调度从线性规划“升级”成了混合整数线性规划求解复杂度也上了一个台阶。第二类是连续变量包括机组出力P(i,t)、储能充电功率P_ch(t)、放电功率P_dis(t)、购电功率P_buy(t)、售电功率P_sell(t)、可中断负荷削减量P_IL(t)。这些变量的定义域是实数区间比如机组出力必须在最小稳定出力50kW到额定200kW之间连续调节。2.2 目标函数五项成本怎么加才完整日前调度最常见的优化目标是整个调度周期内的运行成本最小化。假设时段粒度为1小时调度周期为24小时目标函数写成$$ \min \sum_{t1}^{24} \left[ \sum_{i1}^{2} C_{MT,i}(P_i(t)) \sum_{i1}^{2} C_{su,i}(t) C_{grid}(t) C_{IL}(t) C_{ESS}(t) \right] $$其中每一项的含义C_MT是微燃机的燃料成本通常用二次函数描述a*P^2 b*P c。在MATLAB里为了方便求解会把二次项做分段线性化处理或者直接用二次规划求解器。C_su是启停成本。燃气轮机冷启动一次代价不小论文里通常设置一个固定值比如每次启动120元通过启停状态变量的变化量来触发。C_grid是购售电成本对应分时电价下的购电支出减去售电收益。这里有个细节购电和售电价格通常不对称购电价执行的是用户侧峰谷电价售电价执行的是上网电价后者一般低于前者。C_IL是可中断负荷的补偿成本。削减用户负荷需要付给用户补偿费补偿价格通常高于正常电价所以只有在系统实在平衡不了的时候才会削减。C_ESS是储能的运行损耗成本可以用充放电过程中折算的损耗费用来近似也可以直接写为储能的日均折旧成本乘充放电量。2.3 约束条件功率平衡、储能SOC、机组爬坡这些硬约束建模中最容易出问题的地方在约束条件。基本盘至少要有下面这些功率平衡约束是等式约束也是模型的核心骨架所有时段必须满足负荷功率 储能充电功率 风电出力 光伏出力 微燃机出力 储能放电功率 购电功率 - 售电功率 - 可中断负荷削减量注意储能充电要放在等号左边作为“用电设备”放电放在右边作为“电源”这个方向搞反了模型就废了。机组出力上下限约束当机组运行时出力必须在最小、最大稳定出力之间当机组停机时出力必须为0。这个逻辑需要用大M法或者直接写成另一种形式u(i,t) * P_min(i) P(i,t) u(i,t) * P_max(i)。爬坡约束机组在两个相邻时段之间的出力变化量不能超过爬坡速率限制。燃气轮机升负荷速率和降负荷速率可能不同虽然这样会多写两个不等式但在复现SCI论文时需要留个心眼很多论文是直接统一成同一个爬坡速率。储能SOC递推约束SOC(t) SOC(t-1) (η_ch * P_ch(t) - P_dis(t) / η_dis) * Δt / E_max其中 η_ch 和 η_dis 分别是充、放电效率一般取0.95左右。这里有个容易忽略的点P_ch * Δt的单位是kWh要先除以储能额定容量 E_max 再换算成SOC的百分比增量。如果直接用kW的功率数值去加SOC量纲就错了。储能容量约束SOC必须落在[SOC_min, SOC_max]内论文里一般取10%—90%。充放电功率也要有上限而且还要加一个“不能同时充电和放电”的逻辑约束否则求解器可能会给你算出一个又充电又放电的荒唐解来“刷平衡”。这个逻辑约束可以用两个互补的0-1变量控制也可以用P_ch * P_dis 0这类非线性约束但不推荐会增加求解难度。可中断负荷约束削减量不超过上报的最大可削减量且一天内总削减时段数有限制体现负荷方“可以配合但不会无限配合”的实际约束。备用容量约束系统要留出一部分旋转备用应对预测误差。一般是要求运行中机组最大出力之和减去当前出力加上储能可放电容量大于预测误差的某个比例。这些约束写全了日前模型才算完整。缺了备用约束论文评审时容易被质疑系统可靠性缺了储能同时充放电约束算出来的结果在工程上是不可执行的。2.4 日前调度结果的输出形态日前调度求解完成后输出的核心结果包括各微燃机各时段启停状态矩阵0-1矩阵也是日内调度的固定输入各微燃机各时段出力计划储能各时段充放电功率和SOC变化轨迹各时段购电/售电功率曲线可中断负荷削减计划这些结果构成了一组“基准计划曲线”交给电网调度机构和电力市场。但注意这个基准计划在真实运行时大概率会被修正——这就轮到日内调度出场了。3. 日内调度模型预测控制思想下的滚动修正3.1 日内调度的决策空间与日前有什么不同日内调度不再是重新求解一个独立的24小时优化问题而是采用滚动时域优化的方式在每个调度时刻对未来的一个预测窗口建模求解然后只执行窗口内第一步的指令。窗口长度通常取未来4小时16个15分钟时段也有文献取2小时或更长核心权衡在于窗口太短前瞻性不够窗口太长超短期预测的优势发挥不出来计算量还大。日内调度的关键区别在于日前确定的机组启停状态已经固化作为常数参数传入日内模型不再作为整数变量参与优化。这意味着原本的整数变量被“降级”成了常数模型从MILP退化为LP或者二次规划求解过程几乎瞬间完成。日内调度的决策变量变成了微燃机各时段出力调整量、储能充放电功率调整量、购售电调整量、可中断负荷调整量。这些调整量是连续变量实际出力 日前计划出力 调整量。3.2 滚动时域优化的“只执行第一步”原则滚动优化是模型预测控制MPC思想在调度场景中的工程化应用核心操作可以概括为三个步骤在当前时刻 t0读取最新超短期预测数据未来H个时段的风光出力、负荷、电价求解包含未来H个时段决策变量的优化模型只执行 t0 时刻的决策指令然后等待下一时刻到来刷新数据重新求解为什么要“只执行第一步”因为预测是不断更新的。你在 t0 算出来的“未来4小时最优计划”只是基于当前信息的局部最优判断。等时间走到 t015分钟你获得了最新的实际运行数据和新的超短期预测之前的决策可能已经不是最优了所以必须滚动更新。这一步逻辑特别像导航软件的路径规划——你不会出门前按一次路径规划就一路死磕到底而是每隔几十秒根据实时路况重新计算路径。调度也一样只是刷新周期是15分钟而不是几十秒。在Matlab里滚动优化的主循环框架是这样的horizon 16; % 4小时15分钟粒度 total_slots 96; % 全天96个15分钟时段 for step 1:total_slots % 1. 更新超短期预测数据从预测数据文件中读取当前step之后的horizon个点 pred_wind wind_intraday(step:stephorizon-1); pred_pv pv_intraday(step:stephorizon-1); pred_load load_intraday(step:stephorizon-1); % 2. 构建并求解日内滚动优化模型 [P_MT, P_ESS, P_grid, SOC] solve_intraday_model(...); % 3. 只执行第一个时段的指令更新实际状态 P_exec(step) value(P_MT(1)); % 只取第一步 SOC_actual(step1) update_soc(SOC_actual(step), value(P_ESS(1))); end3.3 与日前计划的衔接偏差约束设计日内调度不是完全推倒日前计划重来的。如果每次滚动都完全按最新预测自由优化日前的机组组合和经济性安排就白做了而且可能频繁出现大范围调整对设备寿命和电网调度都不友好。因此日内模型必须“尊重”日前计划通常通过两种方式实现硬约束方式日内出力不能偏离日前计划超过一定比例比如|P_intra - P_dayahead| delta * P_dayahead。这个方式逻辑简单但缺点也很明显当预测偏差确实很大时硬约束可能使可行域为空模型直接无解。软惩罚方式在目标函数中加入偏差惩罚项例如$$ \min \sum_t C_{adj}(t) \lambda \sum_t |P_{intra}(t) - P_{dayahead}(t)| $$其中C_adj是调整成本lambda是偏差惩罚系数。这个方式更灵活求解器可以自动权衡“跟着日前走”和“调整到更优运行点”之间的经济性。工程上我更推荐软惩罚因为它的可行域始终非空极端情况下多花钱也能保证有解。日内调度的目标函数也可以看作两部分之和调整成本最小化 对日前计划的偏差最小化。调整成本包括储能额外充放电损耗、购售电价格波动带来的额外购电成本、可中断负荷补偿。偏差惩罚的系数 lambda 在不同论文里取值方式不同我复现的时候一般取单位电量电价的一个比例比如0.3—0.5倍的电价。4. MatlabYalmipCplex代码实现从建模到求解的完整链路4.1 环境准备与求解器配置要在Matlab里跑通这套模型我建议的软件组合是MATLAB R2020b及以上 Yalmip工具箱 CPLEX或Gurobi。Yalmip是建模语言层的工具箱它本身不求解问题而是把优化模型转换成底层的求解器能识别的格式。用Yalmip的好处是你不用手动处理矩阵拼装模型写出来跟数学公式几乎一一对应调试效率高很多。安装配置里最容易踩的坑是路径问题。下载Yalmip和Cplex的压缩包解压后必须把整个文件夹用addpath(genpath(...))加进MATLAB路径而且要注意Cplex的Matlab接口文件夹在解压根目录的matlab子目录下不是根目录。安装完成后用一条命令验证yalmiptest看到Cplex那一行显示found说明求解器对接成功。如果显示not found优先检查路径是否添加完整其次检查MATLAB版本和Cplex版本是否兼容。官方支持矩阵里一般会写清楚哪些版本组合能正常工作我实测下来MATLAB 2021b配Cplex 12.10是稳定的组合。4.2 代码架构设计分模块开发的好处这种两时间尺度模型代码不要全部挤在一个脚本里建议按模块拆vpp_scheduling/ main_day_ahead.m % 日前调度主程序 main_intraday.m % 日内调度主程序 data/ % 数据文件 load_forecast_DA.m wind_forecast_DA.m pv_forecast_DA.m price_DA.m intraday_forecast.m models/ build_day_ahead_model.m build_intraday_model.m utils/ yalmip_solve.m plot_vpp_results.m calculate_metrics.m这样拆的好处很明显数据、模型、求解、绘图各管各的日后改参数、换数据、加约束不需要理解整个工程只改对应模块就行。我的习惯是主程序保持非常短几十行顶天了真正的逻辑全部封装在函数里。这跟写论文一样有清晰的层次别人拿到代码也好复核这对学术复现极其重要。4.3 关键约束的Yalmip写法与易错点Yalmip里变量定义很直观二进制变量用binvar连续变量用sdpvar。日前调度模型的核心变量定义如下T 24; % 日前时段数 n_MT 2; % 燃机台数 u binvar(n_MT, T, full); % 启停状态 P sdpvar(n_MT, T, full); % 燃机出力 P_ch sdpvar(1, T); % 储能充电 P_dis sdpvar(1, T); % 储能放电 P_buy sdpvar(1, T); % 购电 P_sell sdpvar(1, T); % 售电 P_IL sdpvar(1, T); % 可中断负荷 SOC sdpvar(1, T1); % SOC增加一个初始点这里有一个我复现论文时遇到的坑SOC变量的长度要设成T1因为递推关系需要用到SOC(1)作为初始状态。如果粗心只定义T个变量最后写SOC约束时索引会越界。功率平衡约束的写法Constraints []; for t 1:T Constraints [Constraints, ... P_load(t) P_ch(t) ... P_wind(t) P_pv(t) sum(P(:,t)) P_dis(t) P_buy(t) - P_sell(t) - P_IL(t)]; end储能SOC递推和容量约束SOC_init 0.2; % 初始SOC 20% SOC_min 0.1; SOC_max 0.9; eta_ch 0.95; eta_dis 0.95; E_max 600; % kWh Constraints [Constraints, SOC(1) SOC_init]; for t 1:T % SOC递推 Constraints [Constraints, ... SOC(t1) SOC(t) (eta_ch * P_ch(t) - P_dis(t)/eta_dis) / E_max * 1]; % 容量限制 Constraints [Constraints, SOC_min SOC(t1) SOC_max]; Constraints [Constraints, 0 P_ch(t) 150, 0 P_dis(t) 150]; end注意递推式里最后乘的那个1是时间步长小时。如果时段粒度是15分钟这里就要乘0.25。单位换算是这个模型里最常见的Bug来源没有之一。有一个快速验证SOC递推对不对的方法让储能只充电一小时算一下SOC增量是不是等于eta_ch * P_ch / E_max跟手算结果对不上就说明单位出问题了。防止储能同时充放电的约束用0-1变量控制z binvar(1, T, full); % 1表示充电状态 M 300; % 大M常数 for t 1:T Constraints [Constraints, 0 P_ch(t) z(t)*M]; Constraints [Constraints, 0 P_dis(t) (1-z(t))*M]; end大M常数取一个足够大的数就行但别取太大否则会引入数值问题。用150充放电功率上限做M就够了——如果有了z(t)1充电上限就是储能额定功率150放电同理。设定求解器参数并求解ops sdpsettings(solver, cplex, verbose, 1, debug, 1); ops.cplex.mip.tolerances.mipgap 0.001; % MIP间隙容差 ops.cplex.timelimit 300; % 求解时间上限 result optimize(Constraints, Objective, ops); if result.problem 0 u_opt value(u); P_opt value(P); SOC_opt value(SOC); else disp(求解失败); disp(result.info); % 开启Yalmip debug后会输出更详细的诊断信息 end4.4 求解结果可视化与指标计算模型跑完之后别急着写论文先画图检查结果合不合理。我最常画的四张图有功功率平衡图负荷、风电、光伏、燃机、储能、购售电的堆叠面积图一眼能看出功率是否平衡、弃风弃光是否异常储能SOC曲线检查SOC轨迹是否在上下限内以及是否呈现“谷充峰放”的合理规律微燃机启停状态计特图用stairs画阶梯图检查机组是否有违反最小启停时间的逻辑问题购售电曲线和分时电价叠加图验证购电行为是否跟随电价信号绘制SOC轨迹的示例代码figure; stairs(0:T, SOC_opt, LineWidth, 1.5); hold on; yline(0.9, r--, SOC上限); yline(0.1, r--, SOC下限); xlabel(时段(h)); ylabel(SOC); title(储能SOC轨迹); grid on;看到SOC曲线在每个循环周期内“涨得上去、放得下来”且没有长时间顶在边界上基本说明模型逻辑是对的。5. 复现中的常见问题与调试经验5.1 不可行解的定位方法做优化建模的同学一定都遇到过求解器直接报“infeasible”的情况。我最开始复现的时候也经常被这个问题卡住后来总结出一套逐步排查的套路。第一步打开Yalmip的debug模式sdpsettings(debug, 1)。Yalmip会尝试定位是哪条约束造成不可行。如果它给出的信息不够明确就用“约束逐步启停法”先把所有约束注释掉验证目标函数和变量定义没问题、能正常求解然后从功率平衡这种核心约束开始一条一条加回来每加一条就求解一次直到某一步突然报infeasible就能锁定问题出在哪条约束上。这个方法虽然笨但非常有效。我遇到最多的infeasible根源是两个SOC递推约束的时间步长单位写错导致储能电量“凭空消失”怎么都平衡不了以及爬坡约束与机组最小出力约束冲突比如两台机组都停机时突然要求一台在下一时段必须满发以满足负荷但爬坡速率限制让它在物理上做不到。5.2 数值尺度问题从kW到p.u.这是一个非常容易被忽略但杀伤力很大的问题。如果一个模型里同时出现kW量级比如3500和小数量级比如0.95的效率、0.1的SOC下限Cplex在数值处理上可能出现精度问题表现为求解器告警“Numerical difficulties”或者解的质量很差。解决办法有两个方向。第一个是统一单位把功率全部折算成MW或p.u.。第二个是调整量级让目标函数值的数量级保持在1000—10000之间不要动不动冲到1e10以上。我在复现论文时倾向于用p.u.建模——以100kVA为基准容量光伏出力、负荷功率、储能功率全部除以基准值SOC本来就在0—1之间这样整个模型的所有变量都在相近的数量级数值稳定性明显改善。5.3 求解时间过长怎么优化日前调度如果只有2台燃机、24个时段Cplex基本几十秒就能给出最优解。但如果扩展到10台机组、96个时段MILP的求解时间可能从一个数量级跳到另一个数量级。这也是SCI论文里做大规模算例时常用的处理手段之一。我的实用建议是设置mipgap工程上最优性间隙0.1%已经足够好没必要非要证明到0.0001%。ops.cplex.mip.tolerances.mipgap 0.001能显著缩短求解时间减少整数变量不是所有设备都需要0-1变量。比如储能充放电状态的互补约束如果只是为了防止同时充放电可以用其他方式替代不一定非要引入额外的0-1变量合理设置初始可行解如果知道某个启停方案是可行的可以先固定这个方案算一个可行解作为MIP的初始下界帮助求解器裁剪搜索树5.4 滚动衔接时的时序索引偏移日内调度模块里最容易出bug的地方是时间索引的管理。因为你有一个“当前时刻step”有一个预测窗口长度horizon还有一个全天总时段数total_slots三套索引一旦没对齐结果就是灾难性的。具体场景举例假设当前是第20个15分钟时段也就是上午9点超短期预测给出的数据是从第20时段到第35时段的风光负荷曲线。在Matlab里如果你直接wind(20:35)去取值会取到第20到35个原始数据点——如果原始数据是从全天开头存的这么取是对的但如果窗口模型内部的索引t1对应的是“当前时刻”而预测数据数组索引是从窗口起点开始存的那映射关系就要做一次偏移。这里的关键点是窗口内t1的预测数据在全局数组中的索引是step。我建议把所有全局数据都带时间戳管理结构体里的.time字段留着写代码时对照时间戳而不是单纯的数组下标能少踩很多坑。另一个我反复遇到的问题日内滚动更新SOC时用的是实际执行的充放电功率而不是优化值。如果直接用value(P_ch)去更新忽略了实际执行过程中的效率损耗和偏差那么跑了几个周期之后模型的SOC状态会和真实状态严重偏离后续滚动优化的起点就错了整个控制逻辑失效。正确的做法是在主循环里单独维护一个SOC_actual变量用实际功率值去更新作为下一轮优化的初始状态。5.5 验证模型正确的三个自检习惯最后分享三个我复现这类模型时的自检习惯看起来简单但真能在关键时刻救命。第一个习惯是“跑一半”。建完模型先别跑全天96个时段只跑前4个时段对应1小时用日粒度数据验证一遍。如果模型逻辑有错小规模问题跑完之后人工验算一遍结果就能发现如果直接跑大规模算例错误结果藏在海量数据里特别难找。第二个习惯是“查平衡”。取任意一个时段手工计算功率平衡等式的左右两边残差看是否在1e-6量级。如果残差偏大大概率是某个约束漏掉了别急着画图先把残差问题解决掉。第三个习惯是“看SOC闭环”。一个调度周期结束后检查SOC是否回到了起点附近或是否满足周期边界条件。如果SOC一路向下或者一路向上说明储能约束方向写反了或者充放电效率的作用方向反了。这类问题用肉眼看SOC图就能一眼识破。有了这套模型搭建和代码实现的经验后再回去读那些顶级SCI里的虚拟电厂多时间尺度调度框架你会发现核心逻辑其实很清晰日前通过MILP解决带整数变量的组合优化日内通过滚动MPC解决带反馈的实时修正两个时间尺度通过机组启停状态和偏差惩罚实现协调。把这个主脉络抓住剩下的就是往框架里填充更复杂的细节——比如多场景随机规划、鲁棒优化界、分布式求解ADMM这些都是在同一个底座上的进阶玩法。我用这套思路复现了几个不同算例系统从单VPP到多VPP互联核心的日前日内衔接机制始终没变过建议大家也照这个路径往下走踩坑会少很多。
做虚拟电厂调度的朋友应该都有同一个感受所谓“最优调度方案”本质上是在跟预测较劲。风电场中午报的出力曲线还是平缓的结果上午十点风向一变实际出力掉了一截光伏更是明显昨天看今天的曲线还是平滑钟形当天飘过来一片云直接凹进去一个坑。这时候你昨天辛辛苦苦算出来的那个24小时全局最优计划放到今天真实运行环境里可能就不是最优了甚至严格来说根本没法执行。这个问题在行业里已经有了标准应对手段——把调度拆成日前调度和日内调度两个时间尺度来做。这也是虚拟电厂多时间尺度调度优化研究的核心思路很多顶级SCI论文的模型框架就是这么搭的。这篇文章我把这块内容掰开揉碎讲清楚为什么非得分两个尺度、两个尺度各自的数学模型怎么建、在Matlab里怎么用YalmipCplex把模型跑起来以及复现过程中最常见的坑在哪。适合正在做虚拟电厂、微电网优化调度方向研究或者准备复现相关SCI论文做毕设/横向项目的同学参考。1. 为什么要把调度拆成“日前日内”两个时间尺度预测误差是根源1.1 虚拟电厂调度的底层逻辑在不确定性中做决策虚拟电厂本身不是一座物理电厂而是把分布式光伏、风电、储能、小型燃气轮机、可控负荷这些分散资源通过聚合管控的方式对外统一表现为一个可调的“发电单元”。这个“聚合”的难点在于里面既有可控性很强的微燃机、储能又有完全看天吃饭的光伏和风电。调度的任务就是在每个决策周期内决定各个设备产出多少功率、储能是充还是放、要不要向电网买电卖电使得整个虚拟电厂在满足负荷需求的前提下经济性最好同时满足各类设备自身的运行约束。问题是所有这些决策必须依赖对未来一段时间风光出力、负荷大小、市场电价的预测。而预测这个东西天然是时间尺度越长、误差越大。1.2 预测精度与时域的必然关系电力系统里预测有个经验规律时间尺度越长预测不确定性的边界越大。日前预测提前24小时的分布式光伏出力均方根误差大概在20%—30%的水平超短期预测未来4小时以内的误差通常能压到5%—10%。风功率预测基本也是这个趋势只是幅度更大。预测类型提前时间典型误差水平用途日前预测24小时20%—30%日前计划、机组组合日内超短期预测4小时以内5%—10%滚动修正、实时跟踪这就产生了一个矛盾如果只做日前调度计划做得很“全局最优”但执行起来因为预测误差大概率走样如果只做日内调度虽然预测准了但那些响应慢、启停成本高的设备比如微燃机、需要预先安排的可中断负荷根本来不及反应。所以行业共识是分层协调日前调度负责“定组合”日内调度负责“调出力”。1.3 两个尺度各管什么定组合与调出力的分工逻辑日前调度提前一天运行以1小时或15分钟为时段粒度覆盖未来24小时。它的核心任务是确定机组的启停状态、储能是否参与充放电调峰、以及各时段与大电网的购售电计划。这些问题包含整数决策变量属于典型的**混合整数线性规划MILP**问题。日内调度当天运行以15分钟或5分钟为粒度预测精度更高、系统状态信息更全。它不再重新决定“哪台机组开、哪台机组停”而是在日前给出的启停状态和计划轨迹基础上用滚动优化的方式调整各设备的出力水平、储能充放电功率和购售电计划。这个世界上已经没有了整数变量前提是日内不调整启停求解速度快很多属于**线性规划LP**问题。用生活类比解释就是你计划一个月后去外地参加一个论坛提前订好机票和酒店这是日前计划当天早上起来查实时路况决定走哪条高速、要不要中途充电这是日内修正。不是非此即彼而是配合着来。2. 日前调度混合整数线性规划建模的完整拆解2.1 决策变量整数变量与连续变量的分工先搭一个典型的虚拟电厂骨架。假设里面包含2台微型燃气轮机MT1、MT2单机容量200kW1个风电场WT额定容量200kW1个光伏电站PV额定容量150kW1个储能系统ESS容量600kWh最大充放电功率150kW1个可中断负荷IL最大可削减量50kW与上级电网的联络线购售电功率上限均为300kW日前调度的决策变量可以分成两类第一类是整数变量典型代表是微燃机的启停状态u(i,t)1表示运行、0表示停机。为什么必须要整数因为启停是一个“是/否”的离散决策——你不能开出0.7台燃气轮机来。这个整数的引入直接让日前调度从线性规划“升级”成了混合整数线性规划求解复杂度也上了一个台阶。第二类是连续变量包括机组出力P(i,t)、储能充电功率P_ch(t)、放电功率P_dis(t)、购电功率P_buy(t)、售电功率P_sell(t)、可中断负荷削减量P_IL(t)。这些变量的定义域是实数区间比如机组出力必须在最小稳定出力50kW到额定200kW之间连续调节。2.2 目标函数五项成本怎么加才完整日前调度最常见的优化目标是整个调度周期内的运行成本最小化。假设时段粒度为1小时调度周期为24小时目标函数写成$$ \min \sum_{t1}^{24} \left[ \sum_{i1}^{2} C_{MT,i}(P_i(t)) \sum_{i1}^{2} C_{su,i}(t) C_{grid}(t) C_{IL}(t) C_{ESS}(t) \right] $$其中每一项的含义C_MT是微燃机的燃料成本通常用二次函数描述a*P^2 b*P c。在MATLAB里为了方便求解会把二次项做分段线性化处理或者直接用二次规划求解器。C_su是启停成本。燃气轮机冷启动一次代价不小论文里通常设置一个固定值比如每次启动120元通过启停状态变量的变化量来触发。C_grid是购售电成本对应分时电价下的购电支出减去售电收益。这里有个细节购电和售电价格通常不对称购电价执行的是用户侧峰谷电价售电价执行的是上网电价后者一般低于前者。C_IL是可中断负荷的补偿成本。削减用户负荷需要付给用户补偿费补偿价格通常高于正常电价所以只有在系统实在平衡不了的时候才会削减。C_ESS是储能的运行损耗成本可以用充放电过程中折算的损耗费用来近似也可以直接写为储能的日均折旧成本乘充放电量。2.3 约束条件功率平衡、储能SOC、机组爬坡这些硬约束建模中最容易出问题的地方在约束条件。基本盘至少要有下面这些功率平衡约束是等式约束也是模型的核心骨架所有时段必须满足负荷功率 储能充电功率 风电出力 光伏出力 微燃机出力 储能放电功率 购电功率 - 售电功率 - 可中断负荷削减量注意储能充电要放在等号左边作为“用电设备”放电放在右边作为“电源”这个方向搞反了模型就废了。机组出力上下限约束当机组运行时出力必须在最小、最大稳定出力之间当机组停机时出力必须为0。这个逻辑需要用大M法或者直接写成另一种形式u(i,t) * P_min(i) P(i,t) u(i,t) * P_max(i)。爬坡约束机组在两个相邻时段之间的出力变化量不能超过爬坡速率限制。燃气轮机升负荷速率和降负荷速率可能不同虽然这样会多写两个不等式但在复现SCI论文时需要留个心眼很多论文是直接统一成同一个爬坡速率。储能SOC递推约束SOC(t) SOC(t-1) (η_ch * P_ch(t) - P_dis(t) / η_dis) * Δt / E_max其中 η_ch 和 η_dis 分别是充、放电效率一般取0.95左右。这里有个容易忽略的点P_ch * Δt的单位是kWh要先除以储能额定容量 E_max 再换算成SOC的百分比增量。如果直接用kW的功率数值去加SOC量纲就错了。储能容量约束SOC必须落在[SOC_min, SOC_max]内论文里一般取10%—90%。充放电功率也要有上限而且还要加一个“不能同时充电和放电”的逻辑约束否则求解器可能会给你算出一个又充电又放电的荒唐解来“刷平衡”。这个逻辑约束可以用两个互补的0-1变量控制也可以用P_ch * P_dis 0这类非线性约束但不推荐会增加求解难度。可中断负荷约束削减量不超过上报的最大可削减量且一天内总削减时段数有限制体现负荷方“可以配合但不会无限配合”的实际约束。备用容量约束系统要留出一部分旋转备用应对预测误差。一般是要求运行中机组最大出力之和减去当前出力加上储能可放电容量大于预测误差的某个比例。这些约束写全了日前模型才算完整。缺了备用约束论文评审时容易被质疑系统可靠性缺了储能同时充放电约束算出来的结果在工程上是不可执行的。2.4 日前调度结果的输出形态日前调度求解完成后输出的核心结果包括各微燃机各时段启停状态矩阵0-1矩阵也是日内调度的固定输入各微燃机各时段出力计划储能各时段充放电功率和SOC变化轨迹各时段购电/售电功率曲线可中断负荷削减计划这些结果构成了一组“基准计划曲线”交给电网调度机构和电力市场。但注意这个基准计划在真实运行时大概率会被修正——这就轮到日内调度出场了。3. 日内调度模型预测控制思想下的滚动修正3.1 日内调度的决策空间与日前有什么不同日内调度不再是重新求解一个独立的24小时优化问题而是采用滚动时域优化的方式在每个调度时刻对未来的一个预测窗口建模求解然后只执行窗口内第一步的指令。窗口长度通常取未来4小时16个15分钟时段也有文献取2小时或更长核心权衡在于窗口太短前瞻性不够窗口太长超短期预测的优势发挥不出来计算量还大。日内调度的关键区别在于日前确定的机组启停状态已经固化作为常数参数传入日内模型不再作为整数变量参与优化。这意味着原本的整数变量被“降级”成了常数模型从MILP退化为LP或者二次规划求解过程几乎瞬间完成。日内调度的决策变量变成了微燃机各时段出力调整量、储能充放电功率调整量、购售电调整量、可中断负荷调整量。这些调整量是连续变量实际出力 日前计划出力 调整量。3.2 滚动时域优化的“只执行第一步”原则滚动优化是模型预测控制MPC思想在调度场景中的工程化应用核心操作可以概括为三个步骤在当前时刻 t0读取最新超短期预测数据未来H个时段的风光出力、负荷、电价求解包含未来H个时段决策变量的优化模型只执行 t0 时刻的决策指令然后等待下一时刻到来刷新数据重新求解为什么要“只执行第一步”因为预测是不断更新的。你在 t0 算出来的“未来4小时最优计划”只是基于当前信息的局部最优判断。等时间走到 t015分钟你获得了最新的实际运行数据和新的超短期预测之前的决策可能已经不是最优了所以必须滚动更新。这一步逻辑特别像导航软件的路径规划——你不会出门前按一次路径规划就一路死磕到底而是每隔几十秒根据实时路况重新计算路径。调度也一样只是刷新周期是15分钟而不是几十秒。在Matlab里滚动优化的主循环框架是这样的horizon 16; % 4小时15分钟粒度 total_slots 96; % 全天96个15分钟时段 for step 1:total_slots % 1. 更新超短期预测数据从预测数据文件中读取当前step之后的horizon个点 pred_wind wind_intraday(step:stephorizon-1); pred_pv pv_intraday(step:stephorizon-1); pred_load load_intraday(step:stephorizon-1); % 2. 构建并求解日内滚动优化模型 [P_MT, P_ESS, P_grid, SOC] solve_intraday_model(...); % 3. 只执行第一个时段的指令更新实际状态 P_exec(step) value(P_MT(1)); % 只取第一步 SOC_actual(step1) update_soc(SOC_actual(step), value(P_ESS(1))); end3.3 与日前计划的衔接偏差约束设计日内调度不是完全推倒日前计划重来的。如果每次滚动都完全按最新预测自由优化日前的机组组合和经济性安排就白做了而且可能频繁出现大范围调整对设备寿命和电网调度都不友好。因此日内模型必须“尊重”日前计划通常通过两种方式实现硬约束方式日内出力不能偏离日前计划超过一定比例比如|P_intra - P_dayahead| delta * P_dayahead。这个方式逻辑简单但缺点也很明显当预测偏差确实很大时硬约束可能使可行域为空模型直接无解。软惩罚方式在目标函数中加入偏差惩罚项例如$$ \min \sum_t C_{adj}(t) \lambda \sum_t |P_{intra}(t) - P_{dayahead}(t)| $$其中C_adj是调整成本lambda是偏差惩罚系数。这个方式更灵活求解器可以自动权衡“跟着日前走”和“调整到更优运行点”之间的经济性。工程上我更推荐软惩罚因为它的可行域始终非空极端情况下多花钱也能保证有解。日内调度的目标函数也可以看作两部分之和调整成本最小化 对日前计划的偏差最小化。调整成本包括储能额外充放电损耗、购售电价格波动带来的额外购电成本、可中断负荷补偿。偏差惩罚的系数 lambda 在不同论文里取值方式不同我复现的时候一般取单位电量电价的一个比例比如0.3—0.5倍的电价。4. MatlabYalmipCplex代码实现从建模到求解的完整链路4.1 环境准备与求解器配置要在Matlab里跑通这套模型我建议的软件组合是MATLAB R2020b及以上 Yalmip工具箱 CPLEX或Gurobi。Yalmip是建模语言层的工具箱它本身不求解问题而是把优化模型转换成底层的求解器能识别的格式。用Yalmip的好处是你不用手动处理矩阵拼装模型写出来跟数学公式几乎一一对应调试效率高很多。安装配置里最容易踩的坑是路径问题。下载Yalmip和Cplex的压缩包解压后必须把整个文件夹用addpath(genpath(...))加进MATLAB路径而且要注意Cplex的Matlab接口文件夹在解压根目录的matlab子目录下不是根目录。安装完成后用一条命令验证yalmiptest看到Cplex那一行显示found说明求解器对接成功。如果显示not found优先检查路径是否添加完整其次检查MATLAB版本和Cplex版本是否兼容。官方支持矩阵里一般会写清楚哪些版本组合能正常工作我实测下来MATLAB 2021b配Cplex 12.10是稳定的组合。4.2 代码架构设计分模块开发的好处这种两时间尺度模型代码不要全部挤在一个脚本里建议按模块拆vpp_scheduling/ main_day_ahead.m % 日前调度主程序 main_intraday.m % 日内调度主程序 data/ % 数据文件 load_forecast_DA.m wind_forecast_DA.m pv_forecast_DA.m price_DA.m intraday_forecast.m models/ build_day_ahead_model.m build_intraday_model.m utils/ yalmip_solve.m plot_vpp_results.m calculate_metrics.m这样拆的好处很明显数据、模型、求解、绘图各管各的日后改参数、换数据、加约束不需要理解整个工程只改对应模块就行。我的习惯是主程序保持非常短几十行顶天了真正的逻辑全部封装在函数里。这跟写论文一样有清晰的层次别人拿到代码也好复核这对学术复现极其重要。4.3 关键约束的Yalmip写法与易错点Yalmip里变量定义很直观二进制变量用binvar连续变量用sdpvar。日前调度模型的核心变量定义如下T 24; % 日前时段数 n_MT 2; % 燃机台数 u binvar(n_MT, T, full); % 启停状态 P sdpvar(n_MT, T, full); % 燃机出力 P_ch sdpvar(1, T); % 储能充电 P_dis sdpvar(1, T); % 储能放电 P_buy sdpvar(1, T); % 购电 P_sell sdpvar(1, T); % 售电 P_IL sdpvar(1, T); % 可中断负荷 SOC sdpvar(1, T1); % SOC增加一个初始点这里有一个我复现论文时遇到的坑SOC变量的长度要设成T1因为递推关系需要用到SOC(1)作为初始状态。如果粗心只定义T个变量最后写SOC约束时索引会越界。功率平衡约束的写法Constraints []; for t 1:T Constraints [Constraints, ... P_load(t) P_ch(t) ... P_wind(t) P_pv(t) sum(P(:,t)) P_dis(t) P_buy(t) - P_sell(t) - P_IL(t)]; end储能SOC递推和容量约束SOC_init 0.2; % 初始SOC 20% SOC_min 0.1; SOC_max 0.9; eta_ch 0.95; eta_dis 0.95; E_max 600; % kWh Constraints [Constraints, SOC(1) SOC_init]; for t 1:T % SOC递推 Constraints [Constraints, ... SOC(t1) SOC(t) (eta_ch * P_ch(t) - P_dis(t)/eta_dis) / E_max * 1]; % 容量限制 Constraints [Constraints, SOC_min SOC(t1) SOC_max]; Constraints [Constraints, 0 P_ch(t) 150, 0 P_dis(t) 150]; end注意递推式里最后乘的那个1是时间步长小时。如果时段粒度是15分钟这里就要乘0.25。单位换算是这个模型里最常见的Bug来源没有之一。有一个快速验证SOC递推对不对的方法让储能只充电一小时算一下SOC增量是不是等于eta_ch * P_ch / E_max跟手算结果对不上就说明单位出问题了。防止储能同时充放电的约束用0-1变量控制z binvar(1, T, full); % 1表示充电状态 M 300; % 大M常数 for t 1:T Constraints [Constraints, 0 P_ch(t) z(t)*M]; Constraints [Constraints, 0 P_dis(t) (1-z(t))*M]; end大M常数取一个足够大的数就行但别取太大否则会引入数值问题。用150充放电功率上限做M就够了——如果有了z(t)1充电上限就是储能额定功率150放电同理。设定求解器参数并求解ops sdpsettings(solver, cplex, verbose, 1, debug, 1); ops.cplex.mip.tolerances.mipgap 0.001; % MIP间隙容差 ops.cplex.timelimit 300; % 求解时间上限 result optimize(Constraints, Objective, ops); if result.problem 0 u_opt value(u); P_opt value(P); SOC_opt value(SOC); else disp(求解失败); disp(result.info); % 开启Yalmip debug后会输出更详细的诊断信息 end4.4 求解结果可视化与指标计算模型跑完之后别急着写论文先画图检查结果合不合理。我最常画的四张图有功功率平衡图负荷、风电、光伏、燃机、储能、购售电的堆叠面积图一眼能看出功率是否平衡、弃风弃光是否异常储能SOC曲线检查SOC轨迹是否在上下限内以及是否呈现“谷充峰放”的合理规律微燃机启停状态计特图用stairs画阶梯图检查机组是否有违反最小启停时间的逻辑问题购售电曲线和分时电价叠加图验证购电行为是否跟随电价信号绘制SOC轨迹的示例代码figure; stairs(0:T, SOC_opt, LineWidth, 1.5); hold on; yline(0.9, r--, SOC上限); yline(0.1, r--, SOC下限); xlabel(时段(h)); ylabel(SOC); title(储能SOC轨迹); grid on;看到SOC曲线在每个循环周期内“涨得上去、放得下来”且没有长时间顶在边界上基本说明模型逻辑是对的。5. 复现中的常见问题与调试经验5.1 不可行解的定位方法做优化建模的同学一定都遇到过求解器直接报“infeasible”的情况。我最开始复现的时候也经常被这个问题卡住后来总结出一套逐步排查的套路。第一步打开Yalmip的debug模式sdpsettings(debug, 1)。Yalmip会尝试定位是哪条约束造成不可行。如果它给出的信息不够明确就用“约束逐步启停法”先把所有约束注释掉验证目标函数和变量定义没问题、能正常求解然后从功率平衡这种核心约束开始一条一条加回来每加一条就求解一次直到某一步突然报infeasible就能锁定问题出在哪条约束上。这个方法虽然笨但非常有效。我遇到最多的infeasible根源是两个SOC递推约束的时间步长单位写错导致储能电量“凭空消失”怎么都平衡不了以及爬坡约束与机组最小出力约束冲突比如两台机组都停机时突然要求一台在下一时段必须满发以满足负荷但爬坡速率限制让它在物理上做不到。5.2 数值尺度问题从kW到p.u.这是一个非常容易被忽略但杀伤力很大的问题。如果一个模型里同时出现kW量级比如3500和小数量级比如0.95的效率、0.1的SOC下限Cplex在数值处理上可能出现精度问题表现为求解器告警“Numerical difficulties”或者解的质量很差。解决办法有两个方向。第一个是统一单位把功率全部折算成MW或p.u.。第二个是调整量级让目标函数值的数量级保持在1000—10000之间不要动不动冲到1e10以上。我在复现论文时倾向于用p.u.建模——以100kVA为基准容量光伏出力、负荷功率、储能功率全部除以基准值SOC本来就在0—1之间这样整个模型的所有变量都在相近的数量级数值稳定性明显改善。5.3 求解时间过长怎么优化日前调度如果只有2台燃机、24个时段Cplex基本几十秒就能给出最优解。但如果扩展到10台机组、96个时段MILP的求解时间可能从一个数量级跳到另一个数量级。这也是SCI论文里做大规模算例时常用的处理手段之一。我的实用建议是设置mipgap工程上最优性间隙0.1%已经足够好没必要非要证明到0.0001%。ops.cplex.mip.tolerances.mipgap 0.001能显著缩短求解时间减少整数变量不是所有设备都需要0-1变量。比如储能充放电状态的互补约束如果只是为了防止同时充放电可以用其他方式替代不一定非要引入额外的0-1变量合理设置初始可行解如果知道某个启停方案是可行的可以先固定这个方案算一个可行解作为MIP的初始下界帮助求解器裁剪搜索树5.4 滚动衔接时的时序索引偏移日内调度模块里最容易出bug的地方是时间索引的管理。因为你有一个“当前时刻step”有一个预测窗口长度horizon还有一个全天总时段数total_slots三套索引一旦没对齐结果就是灾难性的。具体场景举例假设当前是第20个15分钟时段也就是上午9点超短期预测给出的数据是从第20时段到第35时段的风光负荷曲线。在Matlab里如果你直接wind(20:35)去取值会取到第20到35个原始数据点——如果原始数据是从全天开头存的这么取是对的但如果窗口模型内部的索引t1对应的是“当前时刻”而预测数据数组索引是从窗口起点开始存的那映射关系就要做一次偏移。这里的关键点是窗口内t1的预测数据在全局数组中的索引是step。我建议把所有全局数据都带时间戳管理结构体里的.time字段留着写代码时对照时间戳而不是单纯的数组下标能少踩很多坑。另一个我反复遇到的问题日内滚动更新SOC时用的是实际执行的充放电功率而不是优化值。如果直接用value(P_ch)去更新忽略了实际执行过程中的效率损耗和偏差那么跑了几个周期之后模型的SOC状态会和真实状态严重偏离后续滚动优化的起点就错了整个控制逻辑失效。正确的做法是在主循环里单独维护一个SOC_actual变量用实际功率值去更新作为下一轮优化的初始状态。5.5 验证模型正确的三个自检习惯最后分享三个我复现这类模型时的自检习惯看起来简单但真能在关键时刻救命。第一个习惯是“跑一半”。建完模型先别跑全天96个时段只跑前4个时段对应1小时用日粒度数据验证一遍。如果模型逻辑有错小规模问题跑完之后人工验算一遍结果就能发现如果直接跑大规模算例错误结果藏在海量数据里特别难找。第二个习惯是“查平衡”。取任意一个时段手工计算功率平衡等式的左右两边残差看是否在1e-6量级。如果残差偏大大概率是某个约束漏掉了别急着画图先把残差问题解决掉。第三个习惯是“看SOC闭环”。一个调度周期结束后检查SOC是否回到了起点附近或是否满足周期边界条件。如果SOC一路向下或者一路向上说明储能约束方向写反了或者充放电效率的作用方向反了。这类问题用肉眼看SOC图就能一眼识破。有了这套模型搭建和代码实现的经验后再回去读那些顶级SCI里的虚拟电厂多时间尺度调度框架你会发现核心逻辑其实很清晰日前通过MILP解决带整数变量的组合优化日内通过滚动MPC解决带反馈的实时修正两个时间尺度通过机组启停状态和偏差惩罚实现协调。把这个主脉络抓住剩下的就是往框架里填充更复杂的细节——比如多场景随机规划、鲁棒优化界、分布式求解ADMM这些都是在同一个底座上的进阶玩法。我用这套思路复现了几个不同算例系统从单VPP到多VPP互联核心的日前日内衔接机制始终没变过建议大家也照这个路径往下走踩坑会少很多。
