微电网日前日内多时间尺度源荷储协调调度的Matlab+yalmip实现

最近帮一个做微电网方向的朋友把“日前日内多阶段多时间尺度源荷储协调调度”的整套Matlab程序从头到尾跑了一遍从模型搭建、yalmip建模到求解器配置、结果出图算是完整捋顺了。这套东西在电力系统方向的论文里出镜率极高几乎可以用作研究生课题的标准开局风光不确定性要处理、储能充放电要约束、机组启停要决策、还得兼顾经济性和新能源消纳一个模型不复杂但细节非常多。这个项目解决的痛点很明确单靠人工经验排计划根本扛不住光伏、风机和负荷的波动储能也会被频繁充放搞到寿命骤减。而“日前日内多阶段多时间尺度”的思路就是先把“明天怎么跑”的大计划定下来再在运行当天根据最新的预测数据滚动修正相当于开车先看地图规划路线再靠实时路况调整方向。对于刚接触电力系统优化的同学或者正在写微电网、综合能源、园区能量管理方向论文的工程师这套Matlab代码够你直接用也够你改造成自己的模型。我把自己在复现、扩展过程中碰到的问题和修改经验都整理在下面了包括模型怎么设计、代码怎么写、求解器怎么配以及最常见的报错和排查思路。1. 项目整体设计与思路拆解1.1 源荷储协调调度到底在解决什么问题先聊清楚“源荷储”这三个字。源是发电侧包括光伏、风机、微型燃气轮机、柴油发电机延伸一点还能算上从大电网买的电荷是负荷侧包括居民用电、工商业负载以及可以平移的柔性负荷储就是储能电池有时候还包括蓄冷蓄热。所谓“协调调度”就是在满足负荷需求的前提下把源、荷、储三部分的运行计划一起优化出来让整个系统在一天内运行成本最低、弃风弃光最小、或者碳排放最少。这里有个很关键的认识源、荷、储不是独立优化而是耦合在一起的。光伏大发的时候要么给储能充电要么卖给电网负荷高峰的时候储能放电还是机组加出力要看电价和燃料成本谁更划算。所以本质上这是一个多变量、多约束、带时间耦合的优化问题最顺手的方式就是用Matlab加yalmip工具箱建模把目标函数和约束条件写出来丢给求解器去解。做这个项目的时候我遇到最多的问题是把模型理解成“算一道静态题目”。其实调度问题最大的难点在于时间耦合储能电池的SOC状态会从第一个时段一路传递到最后一个时段机组的爬坡约束也要考虑相邻时段的变化。所以模型的维度天生是二维的变量矩阵的规模是“设备数量乘时段数量”这也是后面写代码时容易踩坑的地方。1.2 为什么必须分“日前—日内”两个阶段很多初学者会问既然要优化为什么不直接拿未来24小时的所有数据一次算完非得分“日前”和“日内”这就要说到预测精度和不确定性了。提前一天做的日前调度依赖的是次日负荷和风光出力的预测值这个预测精度通常比较粗误差可能到10%以上。如果直接拿这个粗糙的计划去运行白天实际光照突然下降、或者负荷突然飙升系统就会出现功率不平衡。所以在运行当天要做日内滚动优化每隔一段时间比如4小时重新读一次最新预测数据把未来几个小时的计划重新算一遍只执行第一步的结果等下一轮滚动再更新。这就是“日前计划做粗调、日内滚动做细调”的思路和导航软件先规划全局路线、再根据实时路况调整很像。在代码设计上“多阶段”体现在两层优化之间有信息传递日前调度算出的机组开停机状态、储能整体充放电策略、联络线购电计划会作为日内滚动优化时的基准值。日内层在目标函数里加入了对“偏离日前计划”的惩罚项修正幅度尽量小这样既利用了最新预测信息又不会让机组忽开忽停。表日前调度与日内滚动的典型特征对比对比维度日前调度层日内滚动层优化周期24小时4~6小时滚动窗口控制步长1小时15分钟或30分钟数据来源次日预测值超短期预测修正值决策重点机组组合、储能总计划、购电计划跟踪日前计划、微调出力求解规模较大含启停整数变量较小窗口短鲁棒性要求中高这两层合在一起就是标题里“多时间尺度”的含义。如果进一步扩展到实时控制层还可以把秒级的功率分配做进去但目前大部分论文和项目做到日前加日内两层就已经足够撑起一个完整的研究框架了。1.3 整体框架与数据流设计整体代码框架我习惯分成三大块数据输入、模型求解、结果分析。数据输入包括负荷曲线、光伏出力和风机出力的预测序列以及设备参数额定功率、效率、爬坡速率、储能容量、电价曲线等。模型求解分两层进行先解日前再把日前结果传给日内滚动层。结果分析就是画图、统计成本和消纳率。这个框架里有一个容易被忽略的设计点日前和日内的步长不一定相同。日前调度步长是1小时一天24个时段日内滚动步长可以压缩到15分钟未来4小时的窗口就有16个时段。所以从日前传递到日内的时候需要把日前的1小时结果“插值”或者“保持”成日内对应的基准值序列否则维度对不上代码直接报错。我在实际代码里是这样处理的日内滚动的每个阶段先用repmat函数把日前结果扩展到15分钟精度的序列作为基准值然后只对当前滚动窗口内的变量做优化。这样两层之间的衔接就非常干净也不容易出现索引错位。2. 核心模型拆解与数学表达2.1 日前调度模型的目标函数与约束先说目标函数。日前调度的目标一般是综合运行成本最小我这里考虑了以下几项机组燃料成本用二次函数近似或者分段线性化、机组启停成本、向电网购电费用、储能充放电的损耗成本用来抑制频繁充放、以及弃风弃光的惩罚项。目标函数写出来大概是Objective sum(GasCost) sum(StartCost) sum(GridBuyCost) ... sum(BattCost) sum(CurtailPenalty);其中天然气发电成本每时段按出力P的二次函数计算但求解器处理非线性比较慢我一般用分段线性化或者直接用一次项近似。优化目标是minimize所以惩罚系数可以调大一些让模型优先消纳新能源。等式约束最关键的是功率平衡约束每一时段都要满足P_PV(t) P_WT(t) P_MT(t) P_DE(t) P_GridBuy(t) P_Dis(t) ... P_Load(t) P_GridSell(t) P_Ch(t) P_Curtail(t);注意储能充电、卖给电网、弃电都在等式右边相当于“被消费掉了”。这个等式写反或者漏项是整个模型最容易出错的地方后面我会专门讲怎么排查。不等式约束包括各机组出力的上下限、爬坡约束、储能充放电功率限制和SOC范围、与电网交换功率的上限以及不允许同时买电和卖电可以通过二元变量强制。开机组还要考虑最小开机时间和最小停机时间不过这个在初版代码里我一般先不加等基础模型跑通再说因为最小启停时间会引入大量整数变量求解速度会明显下降。2.2 日内滚动优化模型日内滚动层的模型和日前层结构类似但有两个显著区别。第一它的目标函数里多了一项“跟踪日前计划”也就是说日内优化的核心不是重新制定一个最优计划而是在不偏离日前计划太多的情况下做修正。这在实际运行中很重要因为机组频繁变动不仅影响效率还可能导致计划外调整。所以在日内层目标函数里面加上了到日前计划值的偏差惩罚项变量多出来的这一块用二次项表示。第二日内层不需要重新决策机组启停。日前层已经定好了哪些机组在今天哪些时段开机日内层只在开机状态下调整出力大小这样就大大减少了整数变量求解速度很快。下面是日内滚动循环的基本结构for k 1 : N_windows % 截取当前滚动窗口的预测数据和日前基准值 idx k : k window_len - 1; P_PV_win P_PV_forecast(idx); P_Load_win P_Load_forecast(idx); % 定义窗口内的连续变量以15分钟为步长 P_MT_win sdpvar(1, length(idx)); P_Dis_win sdpvar(1, length(idx)); P_Ch_win sdpvar(1, length(idx)); SOC_win sdpvar(1, length(idx)); % 约束与目标函数... % 求解 optimize(Cons, Obj, options); % 只保存窗口第一个时段的决策结果 P_MT_result(idx(1)) value(P_MT_win(1)); SOC_refresh value(SOC_win(1)); end注意这个“只取第一步结果”的操作这是滚动优化的核心。窗口虽然预测了未来4个小时但真正下发执行的只有当下这一个时段下一轮循环会重新滚动、重新预测这样才叫“滚动”。2.3 储能、机组和可平移负荷的建模细节储能是调度模型里最重要的元件因为它的变量是时间耦合的。SOC的递推公式要写进约束里否则储能就成了一个没有记忆的“充电宝”毫无意义。SOC递推写成for t 1 : T-1 Cons [Cons, SOC(t1) SOC(t) - P_Dis(t)/eta_dis/ESS_cap ... P_Ch(t)*eta_ch/ESS_cap]; end如果不考虑充放电切换时的能量损失这个式子就是动能守恒。但要真正控制储能不出现“又充又放”的浪费行为要么加二元变量要么靠目标函数里的损耗成本自然抑制。初版代码我建议先加损耗项简单有效速度也快。可平移负荷也是多阶段模型里常见的元件比如洗衣机、工业生产线它们的特征是可以从高峰时段挪到低谷时段但必须连续运行H个时段。这个建模比普通负荷麻烦一些需要用0-1变量表示“启动时刻”再加连续运行约束。如果初版代码想跑通我建议先不加可平移负荷把普通负荷作为已知量处理否则整数变量爆炸调试难度直线上升。机组建模里爬坡约束最容易写错。爬坡是“相邻时段出力之差”受限公式是Cons [Cons, P_MT(t1) - P_MT(t) RampUp * dt, P_MT(t) - P_MT(t1) RampDown * dt];这里dt如果是小时爬坡速率单位就是MW/h如果换成15分钟步长爬坡速率要除以4。很多人改完步长忘了改爬坡参数结果日前能解、日内就报无解很大概率是这个原因。3. Matlab代码实现与实操要点3.1 代码整体结构与模块划分完整跑通的Matlab代码我不建议所有内容堆在一个脚本里。我自己习惯分四个文件main.m负责流程调度、load_data.m负责读入和生成数据、build_model.m负责建立yalmip变量和约束、plot_result.m负责画图。如果模型复杂再拆出case_data.m专门存放设备参数矩阵。main.m的主流程大概这样clc; clear; close all; addpath(genpath(./yalmip)); % 1. 数据准备 data load_data(); % 2. 日前调度 [schedule_da, result_da] solve_day_ahead(data); % 3. 日内滚动调度将日前结果作为基准 result_intra solve_intraday(data, schedule_da); % 4. 结果后处理 plot_result(result_da, result_intra);分段写的好处是日后如果要换数据、换模型、换成多目标或者鲁棒优化只需要动对应模块不需要把整个文件重写一遍。这也是很多开源代码的标准结构看代码的人能快速定位问题。3.2 yalmip建模关键步骤用yalmip建优化问题步骤非常固定定义变量写目标函数写约束调optimize求解提取value。最容易出错的地方就是变量维度的定义统一。我建议所有时段的变量都定义成行向量比如P_MT sdpvar(1, 24)这样和约束矩阵相乘的时候不容易转置出错。如果要加机组数量维度第一个维度放机组编号第二个放时段。比如P_mt sdpvar(3, 24)表示3台机组、24个时段然后用sum(P_mt, 1)求总出力。目标函数里要善用sum函数。比如购电成本是分时电价向量Price乘以购电功率向量P_buy直接写成sum(Price .* P_buy)就行Matlab的向量化计算会让代码简洁很多也能避免循环里重复复制变量导致速度慢。一个我反复踩过的坑yalmip里定义约束时sdpvar变量和约束集合要一行一行往下拼接比如Cons []; Cons [Cons, P_MT P_MT_min, P_MT P_MT_max]; Cons [Cons, SOC(2:end) SOC(1:end-1) - P_Dis(2:end)/eta_dis/cap P_Ch(2:end)*eta_ch/cap];拼接的时候要小心逗号分隔的多个约束会合并成一个约束集合但如果有维度不匹配报错会非常让人抓狂提示“Unable to perform assignment”或者“Dimensions of arrays being concatenated are not consistent”。这时候我一般用size(P_MT)挨个检查变量维度。3.3 求解器选择与参数配置这个模型带机组启停和储能充放电方向的整数变量本质是MILP问题。yalmip本身不求解优化问题它只是一个建模语言背后要接求解器。最常用的商业求解器是CPLEX和GurobiMatlab自带的intlinprog也能解适合小规模快速验证。求解器参数设置在yalmip里很简单options sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 0.001, gurobi.TimeLimit, 300);MIPGap控制在0.1%左右足够了取到0.0001会很慢。TimeLimit一定要设尤其数据集很大的时候否则求解器可能跑十几分钟不结束你都不知道该等还是该停。如果没有Gurobi或者CPLEX授权可以用Matlab内置的intlinprog。yalmip里直接写options sdpsettings(solver, intlinprog, verbose, 1);实测下来纯连续变量的优化问题intlinprog表现尚可但一旦整数变量超过几十个速度明显下滑。所以如果打算做实打实的研究还是建议申请Gurobi的学术免费license个人学术用途基本免费网上搜一下就能找到申请入口。3.4 调度曲线可视化与结果导出调度结果不画图等于白做。我常用的画图方式是4个子图放在一张图里第一个子图画电源出力堆叠图第二个画储能SOC和充放电功率第三个画负荷平衡情况或者叫功率不平衡量第四个画与大电网的交换功率和对应电价。这样一张图基本把调度结果的关键信息全涵盖了。画图的时候有几个小技巧。第一堆叠图用area函数而不是plot因为面积图能直观看出出力结构占比。第二负荷平衡曲线要在图上叠加看等式约束是否严格满足如果不满足大概率约束写错了。第三Excel导出用writematrix保存调度计划给后续分析模块用。figure; subplot(2,2,1); area(t, [P_PV; P_WT; P_MT; P_DE]); legend(光伏, 风机, 燃气轮机, 柴油机); xlabel(时间/h); ylabel(出力/MW); subplot(2,2,2); yyaxis left; plot(t, SOC); ylabel(SOC); yyaxis right; stairs(t, P_Dis, r); hold on; stairs(t, -P_Ch, b); legend(SOC, 放电, 充电);日期轴上如果显示“1,2,3,4”而不是具体时刻我一般先把时间轴数据做成分钟或小时数plot完之后再设置xticklabel。4. 常见问题与排查技巧实录4.1 求解无解与结果异常的排查思路这是频率最高的报错yalmip返回的info是“Infeasible problem”意思是当前约束条件下找不到任何可行解。这不代表模型错了但大概率是约束之间互相矛盾。排查顺序我一般固定如下。先看功率平衡等式把每时段的可用出力上限加起来看看能不能覆盖负荷。比如负荷峰值为100MW但光伏风机机组储能最大出力只有80MW那必然无解。然后看储能SOC初值设太低或者容量太小也会导致充放电空间不够把SOC初始值改成50%或者调大容量试试。再看爬坡约束如果机组爬坡速率限制得太死高峰时段的负荷可能跟不上。最后看与大电网的交换功率上限如果买的电不够负荷用同样无解。一个实用的小技巧在约束集合里临时删除一部分约束比如先注释掉爬坡约束看问题是否变可行。通过二分法定位是哪个约束组导致的冲突。这个办法我屡试不爽。4.2 求解速度太慢的优化方法刚开始跑通的时候如果日前的24个时段全部用整数变量比如3台机组加储能方向可能就有72个0-1变量MILP求解时间会从几秒涨到几分钟。优化速度的办法也很多。第一减少整数变量的数量。如果不需要精确的启停成本可以把机组启停固定为日前第一轮求解的结果日内层就全是连续变量速度飞起。第二把SOC的递推约束写得尽量紧凑不要写成三重循环yalmip对向量化约束的支持很好直接用矩阵运算。第三合理设置MIPGap求解器没必要证到“绝对最优”0.5%的gap在很多工程场景下完全够用。第四检查是否有冗余约束比如SOC大于0小于1这种边界条件如果模型本身已经通过递推自然满足加不加意义不大。实测下来同一个模型在Gurobi和intlinprog之间速度差距能到10倍以上。如果就是跑着玩、验证逻辑intlinprog够用如果要做蒙特卡洛或者参数敏感性分析必须上商用求解器。4.3 环境配置类问题yalmip、求解器、中文注释乱码Matlab环境问题是很多新手过不去的坎。yalmip安装其实非常简单从GitHub上下载一个压缩包解压后addpath(genpath(yalmip目录))就行。但如果Matlab版本比较老建议找对应时期的yalmip版本新版yalmip对老版本Matlab的支持可能有问题提示无法加载coder或者找不到classdef。求解器集成也有讲究。Gurobi和CPLEX的安装包里会自带Matlab接口安装完后要把对应的接口路径也加到Matlab路径里然后重启Matlab。验证是否成功可以在命令窗口输入yalmiptest会列出检测到的所有求解器。中文注释乱码这事也是高频搜索词。Matlab 2023a及以后版本默认注解编码可能和旧代码不一致打开旧代码全是乱码。解决办法有两个一是把脚本文件编码从GBK转成UTF-8用VS Code或者记事本另存为UTF-8二是在Matlab的预设选项里把“代码文件编码”设置成默认UTF-8。注意转换编码之前先备份我以前转换中文字符串导致字符串被截断折腾了很久才发现是编码抖动的问题。表常见报错与排查速查表报错/症状可能原因排查和解决办法Infeasible problem约束矛盾负荷超过总可控出力上限删除部分约束定位冲突检查平衡等式和爬坡限制Solver not foundyalmip路径没加好或求解器未安装运行yalmiptest查看可用求解器列表NaN或Inf出现在结果里某个变量没有被约束到超出定义域检查所有sdpvar是否都为0索引是否越界维度不匹配报错变量向量方向不一致1x24 vs 24x1统一向量方向用size打印变量维度中文注释乱码GBK/UTF-8编码不一致用编辑器转换编码或修改Matlab预设编码求解时间过长整数变量多MIPGap设置过小固定机组启停、放宽MIPGap、换Gurobi结果出图时曲线断层日前和日内步长不一致导致衔接错位日前结果先按15分钟粒度扩展再传递储能SOC长时间满或空充放电效率参数设置不合理或惩罚项不足调整效率增加损耗成本系数4.4 结果合理性的人工校验求解器给出结果不代表结果就一定合理这一点我吃过很多亏。所以每次跑完我都习惯人工验算几个关键指标功率平衡是否每个时段都严格满足、SOC首尾是否一致如果设定了调度周期末SOC要回到初始值看约束有没有生效、机组出力有没有抄过额定值或频繁启停。另一个很隐蔽的问题如果目标函数里弃风弃光惩罚系数设得太低模型会为了省燃料成本故意弃掉便宜的风光导致新能源消纳率很低。如果惩罚系数设得太高模型又可能为了让储能多充电买高价电。这种目标函数系数怎么调没有绝对标准只能根据你的调度偏好来定。我的经验是弃电惩罚系数通常设置成购电最高电价的1.5~2倍这样模型会优先消纳新能源而不是贪便宜买电。5. 扩展方向与下一步建议5.1 从单目标走向多目标调度基础版本是单目标即成本最小。实际论文里往往要求多目标比如同时考虑运行成本和碳排放。改成多目标有两条路一是用加权和把碳排放折算成碳成本加进目标函数二是用NSGA-II这类启发式算法求Pareto前沿。前者在yalmip框架里几乎是零改动后者需要把优化问题打包成一个函数再套用遗传算法框架。如果时间紧我建议先走加权和这条路。把碳排放系数乘以机组出力得到一个“碳成本”项和燃料成本加权求和。权重系数慢慢调可以画出一条Pareto曲线论文图和故事性都够了。5.2 引入不确定性建模日前调度的输入是预测值但预测不准确怎么办最常见的做法是场景法或者鲁棒优化。场景法就是把光伏、负荷的可能情况生成多个场景每个场景对应一组约束让解对所有场景都可行鲁棒优化则是把不确定量控制在某个区间里求解最坏情况下的最优解。这些方法在Matlab里都可实现但计算量会显著上升。我给的建议是先把确定性模型跑熟再逐步加不确定性。否则遇到无解的时候你根本分不清是场景生成问题还是原始模型的问题。5.3 与深度学习预测结合还有一个很好的扩展方向用Matlab的深度学习工具箱训练光伏和负荷预测模型把预测结果作为日前调度的输入。这是因为调度结果的优劣很大程度上取决于预测精度目前学术界和工程界都很看重“预测优化”的联合框架。比如用BP神经网络或者LSTM做完预测再把预测序列直接喂给前面这套调度程序整个代码流程已经天然打通了。Matlab里用trainNetwork训练LSTM非常方便预测结果存在矩阵里后续直接索引取用即可。这样一套“预测—调度”联动的代码框架无论是毕业设计还是实际项目都很有竞争力。5.4 从离线优化走向实时闭环如果程序已经跑通还希望落到实际系统中可以考虑把调度结果做成查询表或者在线应用。实时控制层接收实际测量数据后在短时间内重新求解一个简化的功率分配问题。这个实时层不需要含启停整数变量只做连续功率分配求解时间可以压到几十毫秒级别和PLC或SCADA联动是完全可行的。我自己在测试这套方案时把实时层的优化频率设置为每分钟调用一次滚动窗口设定为未来5分钟、步长1分钟。结果和离线优化的差距很小但应对突发波动的能力明显提升。这块算是从“仿真程序”走向“实际系统”的一道分水岭。最后再分享一个小技巧。多阶段调度程序的调试一定要从单层开始。先把日前调度跑通画出正确的出力和SOC曲线再加日内滚动。否则两层耦合在一起出问题的话你根本找不到是日前结果不对、还是滚动逻辑有误。另外数据文件的命名和列顺序一定要固定下来没人愿意在写完模型之后还要回头猜测第5列到底是负荷还是光伏。