★ FEATURED ARTICLE 简介这份源程序围绕虚拟电厂/微网单元的日前随机优化调度展开面向电气工程、能源系统方向的毕业设计及科研学习者重点解决光伏出力与负荷功率双重不确定性下的调度决策问题。资源包以RAR压缩形式提供整体约3.66MB已有87人学习下载。程序采用蒙特卡洛算法对预测的光伏和负荷曲线进行多场景生成再利用快概率距离快速消除法将场景削减至5个构建了计及源荷双重不确定性的虚拟电厂随机优化调度模型并通过随机规划法对多场景下的调度策略进行优化实现效果良好。代码覆盖从场景生成、削减到优化求解的完整流程可帮助读者快速复现论文算例理解随机规划在能源系统中的应用逻辑也可作为毕业设计或课题研究的改造基础。1. 日前随机优化调度把预测会错本身写进模型而不是事后补救微网日前调度里有一道分水岭把风电、光伏和负荷的预测值当成确定数去优化还是把预测会错这件事本身写进模型。前者的结果在数值预测稍有偏差时就会失效后者的代表做法就是标题里《计及源荷双重不确定性的虚拟电厂微网日前随机优化调度》这套论文加源程序。源程序解决的具体问题可以一句话讲清用场景法把光伏、风机出力与负荷的不确定性转成一组组可能的未来曲线再在期望成本最小化的目标下求日前计划让计划在所有场景下都可行又不至于保守到浪费资源。适合正在做虚拟电厂、微网日前调度课题的同学和工程师既得看懂随机优化建模的套路还得能把源程序跑起来改掉参数后能复现出论文表里的结果。下文先讲建模骨架再给可抄的MATLAB实现最后是排查记录。2. 源荷双重不确定性的场景建模从点预测到1000条可能的未来2.1 为什么点预测会骗人不确定性要靠场景来描述常规的确定性日前调度是把第二天的光伏、负荷预测曲线直接塞进优化模型。真实曲线不可能是预测值那条光滑线当光照被云遮住、负荷峰谷偏移功率平衡约束可能直接被击穿这时候只能临时依赖实时市场补电或紧急切负荷成本完全不可控。把不确定性写进模型的第一步是用概率分布描述预测误差对光伏来说预测误差的标准差通常取预测值的10%~20%用正态或截断正态分布负荷预测相对靠谱标准差取3%~5%就够风机出力误差往往偏态明显用Beta分布拟合更合适。场景法的做法就是按这些分布反复采样生成例如1000条完整的24点出力曲线每条曲线配一个等概率。纯蒙特卡洛实现最简单但收敛慢拉丁超立方采样把区间分层低位收敛明显更快。如果光伏和负荷误差存在负相关或同区域风电场出力受地形牵连可以用Cholesky分解把协方差矩阵折算进采样过程。下面是生成场景的最小脚本我习惯单独存成gen_scenarios.m% gen_scenarios.m % 生成光伏出力与负荷的日前场景T24点 % 输出: P_pv_scen [N_scen, T], P_load_scen [N_scen, T] N_scen 1000; T 24; sigma_pv 0.15; % 光伏预测误差标准差比例 sigma_load 0.05; % 负荷预测误差标准差比例 % P_pv_forecast / P_load_forecast: 预测曲线1xT由外部数据读入 P_pv_scen zeros(N_scen, T); P_load_scen zeros(N_scen, T); for s 1:N_scen for t 1:T % 截断到[0, 装机容量]避免出现负出力 P_pv_scen(s,t) max(0, P_pv_forecast(t) ... sigma_pv * P_pv_forecast(t) * randn()); P_load_scen(s,t) max(0, P_load_forecast(t) ... sigma_load * P_load_forecast(t) * randn()); end end这段代码把所有随机性压在两行randn上。sigma_pv取0.15的意思是预测100MW的光伏出力场景的标准差是15MW24个点各自独立扰动。实际工程里误差在时间上是相关的比如连续两朵云不会让第1小时偏大、第2小时立刻偏小。想更真实可以对误差做一阶自回归AR(1)e_t rho*e_{t-1} sqrt(1-rho^2)*sigma*randnrho取0.5~0.8。负荷的sigma取0.05已经偏保守再大会出现负负荷这种脱离物理的场景。注意场景数和点数决定了后面场景削减和MIP求解的规模1000×24只是入门量级别一上来就上10000个场景后面削减代码的O(N²)距离矩阵会先把你卡住。2.2 随机优化的目标函数期望成本不是一次成本确定性模型是一个场景下的总成本随机模型是所有场景下成本的期望。数学形式可以写成最小化第一阶段成本加上所有场景第二阶段成本的概率加权和。第一阶段成本在日前就要定死不随场景变第二阶段成本是每个场景下调整量带来的代价。常见的成本项大概分这几类成本项归属阶段说明上级电网购电成本一按日前电价 × 计划购电量机组燃料成本一线性或分段线性近似机组启停成本一与0-1变量绑定弃风弃光惩罚二每弃1MWh给惩罚通常取几百元/MWh切负荷惩罚二要设得很大否则模型会主动切负荷省钱调整量惩罚二场景内实际出力偏离计划出力的代价这里有一个很容易翻车的点切负荷惩罚系数。很多初版模型把切负荷惩罚设为500元/MWh结果随机优化求解出来的计划在稍微不利的场景就果断切负荷总成本反而比确定性模型漂亮很多。这不是模型变强了是模型在作弊。常见做法是参考失负荷价值量级至少取到5000~10000元/MWh让切负荷只在极端场景下出现。弃风弃光的惩罚也要高于机组发电边际成本不然模型宁可弃掉也不发。2.3 两阶段变量分工哪些变量不带场景下标随机优化的核心语法不在求解器而在变量维度的声明。第一阶段变量不带场景下标包括机组启停0-1变量u、机组计划出力P_g、储能日前充放电计划、与上级电网的购电计划。第二阶段变量带场景下标包括机组出力调整量、弃风弃光量、切负荷量、储能的实时修正。在YALMIP里我习惯用cell数组给每个场景单独建变量避免三维sdpvar在部分版本上的兼容问题% 第一阶段不随场景变化 u binvar(N_gen, T, full); % 机组启停0-1变量 P_g sdpvar(N_gen, T, full); % 机组计划出力 P_buy sdpvar(1, T, full); % 向电网购电 % 第二阶段随场景变化用cell数组存放 delta_g cell(N_scen, 1); % 机组调整量 P_curt cell(N_scen, 1); % 弃风弃光量 P_shed cell(N_scen, 1); % 切负荷量 P_dis cell(N_scen, 1); % 储能放电 P_ch cell(N_scen, 1); % 储能充电 for s 1:N_scen delta_g{s} sdpvar(N_gen, T, full); P_curt{s} sdpvar(1, T, full); P_shed{s} sdpvar(1, T, full); P_dis{s} sdpvar(1, T, full); P_ch{s} sdpvar(1, T, full); endbinvar和sdpvar的第一个参数都是矩阵维度full告诉YALMIP这个矩阵不是对称的。cell数组的做法在循环里组装约束时非常直观。需要特别提醒储能充放电两个变量之间要加互斥约束否则优化器会让它们同时为正、在目标函数里互相抵消白赚效率损耗——这是最隐蔽的建模错误之一后面排查章会单独提。2.4 约束骨架功率平衡、爬坡、SOC与旋转备用这节约束是照着源荷双重不确定性逐个落地的。对每个场景s、每个时段t功率平衡必须严格成立Σ P_g(i,t) P_pv_scen(t,s) P_dis{s}(t) - P_ch{s}(t) P_buy(t) P_load_scen(t,s) P_curt{s}(t) - P_shed{s}(t)注意P_buy不带场景下标因为购电计划在日前就确定弃风和切负荷都放在等式右侧作为松弛手段但必须被惩罚系数套住。旋转备用约束是随机模型和确定性模型的另一个分水岭Σ (P_g_max - P_g(i,t)) ≥ r_load × P_load_scen(t,s) r_pv × P_pv_scen(t,s)这条约束必须带场景下标用来保证即使某场景风光全无也能用热备用顶起来。r_load取5%~10%r_pv取10%~15%是常见起步值。储能SOC递推约束写成SOC(t1,s) SOC(t,s) η_ch × P_ch{s}(t) - P_dis{s}(t)/η_dis同时加上SOC(T,s) SOC(1,s)保证一个调度周期内储能不能凭空生电。很多初版模型会漏掉SOC终值约束储能会在最后几小时疯狂放电把能量用光看起来成本很低实际不可运行。把这几组约束写全模型骨架就立住了。3. 场景削减与两阶段随机优化把1000个场景压到80个还不丢精度3.1 为什么1000个场景不能直接喂给混合整数规划把1000个场景、每个24时段放进MIP约束数量约为每个场景每个时段一条功率平衡、每条机组爬坡、SOC递推等轻松到5万~10万条0-1变量有N_gen×T个不随场景数变。真正膨胀的是连续变量1000场景×24时段×5类调节变量就是12万个变量。对大MIP来说12万变量未必不能解但配上0-1变量形成的分支定界树求解时间会从几分钟膨胀到几小时甚至内存溢出。这就是场景削减存在的理由用几十个有代表性的场景去近似上千个场景代表的概率分布代价是引入一点近似误差收益是求解时间从小时级降到分钟级。3.2 可抄的场景削减代码同步回代消除SCR同步回代消除是最常出现在这类论文里的削减方法。思路是迭代合并每次找到距离最近的两个场景把概率小的那个场景的概率累加到概率大的场景上然后删掉概率小的场景直到剩下N_keep个。它不是简单随机抽而是尽量保留概率质量大和位置极端的场景。function [idx_keep, p_keep] scr_reduction(P_scen, N_keep) % 同步回代消除场景削减 % 输入: P_scen [N_scen, T*2] 每行一个场景列是24点光伏24点负荷 % N_keep 削减后保留的场景数 % 输出: idx_keep 保留的场景编号, p_keep 归一化后的概率 N size(P_scen, 1); p ones(N, 1) / N; % 等概率场景若原始场景带权重则替换此取值 idx_act (1:N); % 当前还活着的场景编号 while length(idx_act) N_keep nact length(idx_act); % 计算活场景之间的欧氏距离矩阵 D zeros(nact, nact); for i 1:nact for j i1:nact d norm(P_scen(idx_act(i),:) - P_scen(idx_act(j),:)); D(i,j) d; D(j,i) d; end end D(D0) inf; % 跳过自身 % 找距离最近的一对 [i_min, j_min] find(D min(D(:)), 1); % 概率小的场景并到概率大的场景上 if p(idx_act(i_min)) p(idx_act(j_min)) p(idx_act(i_min)) p(idx_act(i_min)) p(idx_act(j_min)); idx_act(j_min) []; else p(idx_act(j_min)) p(idx_act(j_min)) p(idx_act(i_min)); idx_act(i_min) []; end end idx_keep idx_act; p_keep p(idx_keep) / sum(p(idx_keep)); % 重新归一化 end这个函数最需要注意的坑是索引D是活场景之间的矩阵i_min和j_min是D的下标要经过idx_act转回全局场景编号直接拿着i_min去删P_scen的行是错的。另一个坑是while循环里每次都对活场景重算距离矩阵整体复杂度接近O(N²×T)N超过2000时建议先用k-means粗削减到500再进SCR。等概率场景直接p1/N如果原始场景来自带权重的采样把p向量换成对应权重即可。调用方式是把光伏和负荷场景拼接成一列方便算距离P_all [P_pv_scen, P_load_scen]; % 拼成 [N_scen, 48] [idx_keep, p_keep] scr_reduction(P_all, 80); P_pv_red P_pv_scen(idx_keep, :); P_load_red P_load_scen(idx_keep, :);N_keep的经验值是50~100。取20时结果往往只剩几个代表遇到强非线性约束会失真后面代价就大了。3.3 削减质量的检验削减前后期望成本偏差要小于2%削减完别急着跑优化。先做一个便宜的自检把削减前的场景集合和削减后的场景集合各算一次期望成本两个值的相对偏差超过2%就说明N_keep选小了。更严格的做法是计算两个概率测度之间的Wasserstein距离但工程上期望成本偏差更直观cost_before evaluate_expected_cost(P_scen); % 削减前期望成本 cost_after evaluate_expected_cost(P_red); % 削减后期望成本 dev abs(cost_before - cost_after) / cost_before * 100; assert(dev 2, 场景削减偏差过大请增大N_keep);如果削减后期望成本明显偏低往往不是好事被删掉的场景里可能有决定备用方案上限的关键场景比如光伏接近零的阴雨天被合并掉了。建议削减后人工检查保留场景里是否覆盖了最坏情况光伏接近装机10%以下、负荷接近峰值的场景至少各留1~2个。evaluate_expected_cost可以是简化版先对每个场景求确定性成本再求加权平均不需要跑完整MIP。3.4 求解器配置与时长管控YALMIP调用CPLEX/Gurobi的参数模型和场景都准备好后求解配置直接决定你能不能下班。我常用的配置如下ops sdpsettings(solver, cplex, ... verbose, 2, ... savesolveroutput, 1, ... cplex.mip.tolerances.mipgap, 0.001, ... cplex.mip.tolerances.integrality, 1e-6, ... cplex.threads, 8, ... cplex.timelimit, 3600); sol optimize(Constraints, Objective, ops); if sol.problem 0 fprintf(最优目标值: %.2f\n, value(Objective)); else disp(sol.info); endmipgap0.001意思是当上下界差距小于0.1%就停止这是把求解时间从几小时压到几分钟的核心参数。timelimit给3600秒兜底防止它周末还在转。verbose2可以看到分支定界的gap曲线卡住时能分辨是收敛慢还是不可行。integrality保持1e-6就好别追求变量精确到1e-9而增大计算量。求解器授权问题经常卡住新手优先检查solver字段拼写、确认求解器安装目录在系统PATH里。sol.problem返回0表示求解成功其余情况打印sol.info看具体原因。4. 随机优化模型的常见问题排查5个踩坑记录与解决路径4.1 现象模型直接报infeasible根本没有可行解原因随机模型里功率平衡约束在极端场景下无解最常见是傍晚高峰负荷叠加零光伏的时段机组出力和储能全放也不够。也可能是P_buy被写成定值没给购电留上下限调节空间。解决先给功率平衡加松弛变量把等式拆成和各加非负松弛松弛量在目标函数里乘以很大的惩罚系数比如10000元/MWh保证平时不松弛、只在极端场景松。调试阶段用温和场景测试——把场景数临时削减到50看是否场景太极端导致无解。还有一招把SOC终值约束暂时去掉如果立即变可行说明储能初值设定有问题。4.2 现象求解器报license错误或求解到一半内存溢出原因license错误通常是求解器安装和系统环境变量没配对YALMIP找不到求解器接口内存溢出则几乎都是场景削减没做、变量维度过大尤其是SCR函数里距离矩阵是N×N2000场景就是几百万个浮点数。解决先单独测solver参数确认YALMIP能识别求解器内存问题先砍N_keep到100再把SCR里的norm改成向量化计算最后检查约束里是否误用了矩阵平方这类非线性算子——YALMIP会把它们转成二阶锥约束维度直接翻倍内存和耗时都起飞。4.3 现象场景削减到20个成本低得可疑原因N_keep太小SCR把高成本场景合并进低成本场景概率被摊薄期望成本系统性偏低。解决对N_keep做扫描例如[20, 50, 80, 100, 150]画出期望成本随N_keep的曲线找拐点同时检查保留场景中光伏最低值是否被删掉必要时强制把最坏场景加入保留集。我见过不少把N_keep直接设成20就交差的版本复现时一对比论文数值就露馅。4.4 现象储能SOC一天下来不回起点充放电效率被白算原因漏了SOC循环约束优化器有动机在最后时段把储能放空用把昨天的能量今天用掉来压成本另一个是充放互斥约束缺失充电和放电同时为正目标函数里一加一减互相抵消但不会多花钱模型会利用这个漏洞。解决两条约束都补上——SOC循环约束写进Constraints充放互斥用P_ch P_dis 某上限 × binvar的方式或者直接加大惩罚。调试时打印value(SOC)看是否从0开始、最后回到0这是最直观的检查手段。4.5 现象随机优化成本比确定性模型还低解释不通原因十有八九是切负荷惩罚太低或者旋转备用约束没带场景下标。第一次搭随机模型时容易犯这个错切负荷惩罚设成800元/MWh随机模型在不利场景果断切负荷总成本反而比不敢切负荷的确定性模型低。解决把切负荷惩罚提到5000元以上同时给旋转备用加最低功率爬坡约束再用同一组场景跑确定性模型做对照两个模型约束完全一致只差不确定性处理方式对比才有意义。5. 源程序到手后先做这四步验证复现数值、压场景数、做对照、灵敏度第一步从文献库把论文正文和附录下下来先对着算例参数表把源程序里的装机容量、负荷曲线、电价、sigma值一项项对齐。跑一遍原算例记录目标函数值和论文表格里的数值对比偏差在3%以内说明手里的版本和发表版一致偏差过大先查场景生成是否用了同一随机种子再查成本项系数有没有漏项。第二步把N_keep从200、150、100、50、20往下扫观察期望成本变化。如果成本在100以后基本平住就在100附近定参如果一直往下掉说明削减偏差一直存在要回头查SCR实现。这个实验顺便能生成一条成本-场景数曲线写进论文里就是很好的敏感性分析图。第三步用同一套参数跑确定性模型再把随机优化的解拿到确定性模型里评估实际成本。两模型成本差距越大说明考虑双重不确定性的价值越大如果差距很小反过来说明这个算例里不确定性本身不强选题需要换场景。第四步改旋转备用系数和切负荷惩罚两个敏感参数各算一遍观察成本变化方向和幅度用来支撑结论里的鲁棒性论述。A同学在复现某同类源程序时曾把场景数从100压到10成本低了8%一度以为是随机优化的功劳后来查出来是备用约束少写了一个场景下标——这类源程序的坑大多不在算法本身而在约束细节尤其是带不带场景下标这一条。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
