做热电联供微网优化有一段时间了从最初的确定性调度一步步做到考虑源荷随机特征的随机优化踩了不少坑也积累了一些能直接用的经验。这篇就围绕“考虑源荷随机特征的热电联供微网优化研究Matlab代码实现”这个方向把整个建模到代码落地的链路完整拆开讲一遍。适合正在做微网优化调度、综合能源系统方向课题的学生也适合刚入行做园区能源规划、想了解随机优化怎么落地的工程师。内容会涉及模型思路、数学公式怎么转成代码、场景生成与削减的实操技巧以及我自己在跑仿真时遇到的一堆典型问题。1. 热电联供微网优化到底在解决什么问题1.1 先搞清楚系统的“源”和“荷”长什么样热电联供微网简单说就是一个麻雀虽小五脏俱全的小型能源系统。它通常包含发电设备微型燃气轮机、内燃机、供热设备燃气锅炉、余热锅炉、可再生能源光伏、风电、储能装置蓄电池、蓄热罐以及电负荷和热负荷两类用户需求。这个系统的核心价值就在于用一套设备同时满足电和热两种需求实现能源的梯级利用。很多人一开始容易把“优化”想复杂。其实本质上就是回答一个问题在未来24小时里每一台设备在每一个时刻该发多少电、产多少热、充放多少储能才能让总运行成本最低如果光伏、风电出力是确定已知的负荷也是固定数值这是个标准的混合整数线性规划MILP问题直接丢给求解器就能出结果。但现实没那么理想。光伏靠天吃饭云一遮就断崖式下跌风电更是看风速脸色用户的电热负荷也受天气、作息、生产计划影响波动很明显。如果调度方案完全基于预测的确定性数值来制定实际运行时很容易出现功率缺额或者能源浪费。这就引出了标题里的核心关键词——源荷随机特征。1.2 随机特征怎么进入优化框架处理源荷随机特征工程和学术上常用的路线有三条。第一条是随机规划Stochastic Programming用多个离散场景描述可能的出力与负荷情况让优化结果在所有场景下综合表现最优。第二条是鲁棒优化考虑最坏情况追求“不管怎么波动都不会出事”。第三条是分布鲁棒优化介于两者之间考虑概率分布的不确定性。我这次实现用的是随机规划中的场景法也是目前最直观、最容易在Matlab里落地的一种。核心思路分三步先对光伏出力、风电出力、电负荷、热负荷的预测误差分别建模用蒙特卡洛抽样生成大量场景然后用场景削减技术去掉冗余场景保留少数几个有代表性的场景最后把这几个场景一起塞进优化模型目标函数变成“所有场景下运行成本的概率加权求和”约束条件则要保证每个场景下系统都能安全运行。这样做的好处很明显计算复杂度可控求解器处理起来不费劲同时能让调度决策在随机环境下保持可行性不为单一预测值“赌运气”。缺点也有就是结果对场景质量比较敏感。场景生成得不好优化结果就会偏离实际这一块后面会详细讲。2. 目标函数与约束条件的建模细节2.1 目标函数不是只算发电成本做热电联供微网优化第一步就是把目标函数写清楚。注意这里的“成本”不是单一的电费而是一整套运行费用的总和。我常用的目标函数形式如下min F Σs ωs × Σt [ CF(t) CM(t) CB(t) CG(t) CP(t) ]其中s是场景索引ωs是场景s的概率权重CF(t)是t时刻的燃料成本包括燃气轮机的天然气消耗和燃气锅炉的天然气消耗CM(t)是设备维护成本通常按出力线性折算CB(t)是储能设备劣化成本用于避免蓄电池和蓄热罐频繁充放CG(t)是向大电网购电的成本如果允许向电网售电还需要扣掉售电收益CP(t)是弃风弃光惩罚项用来保证优化不会把可再生出力全砍掉来省成本。很多初学者忽略了一个细节热电联产机组的发电和发热不是独立的两者之间存在强耦合关系。同一个机组发电出力定了可回收的余热量就有个上限不是想发多少热就发多少热。这个耦合关系通常用一个“热电比”系数来描述具体实现时要查机组的技术参数手册或者用厂家提供的性能曲线拟合出可行域不等式。还有一个点容易被忽略——机组启停状态。燃气轮机频繁启停不仅损耗设备寿命还会抬升成本。所以模型里要引入二进制变量表示机组启停状态并且配合最小启停时间约束。这就把模型从LP变成了MILP求解时间会明显增加。后面在调试时MILP求解时间是大家问得最多的问题之一。2.2 约束条件的完整清单与物理含义约束条件是优化模型的骨架每一条都对应着实际的物理限制。我在写代码时会把约束分成四大类方便检查和调试。第一类是功率平衡约束。电功率平衡要求任意场景、任意时刻燃气轮机发电、光伏出力、风电出力、蓄电池放电、电网购电之和恰好等于电负荷加上蓄电池充电和电转热设备如果有的话的耗电。热功率平衡要求燃气轮机余热、燃气锅炉产热、蓄热罐放热之和等于热负荷加上蓄热罐充热。要注意电平衡里电网交互功率是变量既可以是正值购电也可以是负值售电很多论文会用正负两个变量来线性化处理避免目标函数里出现绝对值项。第二类是设备自身约束。比如燃气轮机出力的上下限约束、爬坡约束相邻时刻出力变化不能超过额定值的一定比例、最小启停时间约束。光伏和风电的出力上限就是该时刻场景提供的最大可用功率这个值随场景和时刻变化。燃气锅炉同样有出力上下限和爬坡约束。蓄电池和蓄热罐需要重点处理——不仅要限制充放电功率上限还要保证存储能量的动态递推关系当前时刻能量等于上一时刻能量加上充入量减去放出量并考虑充放效率损失。第三类是储能状态约束。储能容量要在安全范围内不能过充过放。我一般还会加一个调度周期始末能量一致的约束也就是日始日终的SOC荷电状态要相等保证调度方案可以循环执行。这个约束在带储能的微网优化里几乎是标配。第四类是随机性相关约束。在场景法中每个场景都要满足前述所有平衡和运行约束。也就是说不管哪个场景发生了系统都必须能安全运行。这就是所谓的“非预期性约束”和“场景可行性约束”。有的论文还会加机会约束比如“热负荷供应不足的概率不超过5%”这种形式在Matlab里实现会更复杂通常需要做转化处理这次先不展开。把上述约束全部堆进模型后整个优化问题的规模会达到几千个变量、几千条约束。对MatlabYalmip求解器来说不算大事但你得保证建模型的时候矩阵别搞错不然一跑起来全是不可行解。3. 源荷随机场景的生成与削减实操3.1 用蒙特卡洛模拟生成原始场景随机优化效果好不好场景生成是第一步。我的做法是先对光伏、风电、负荷的预测值加上带统计特征的误差项。误差的分布不是随便拍脑袋定的得基于历史数据和天气预报水平。光伏误差一般用Beta分布描述风电误差常用Weibull分布负荷预测误差在很多文献里被近似为均值为零的正态分布也有人用t分布来捕捉厚尾特征。具体操作时我一般先对每个时段独立生成误差样本按以下方式叠加到预测值上P_scenario(t, s) P_forecast(t) ε(t, s)这里 ε(t, s) 就是从误差分布中抽出来的随机扰动。但这里有个隐患直接逐时段独立抽样会让生成的场景波动非常剧烈相邻时段忽高忽低完全不符合实际出力曲线的时间连续性。解决这个问题可以在误差抽样时引入时间相关性比如用一阶自回归模型AR(1)来生成平滑的误差序列或者用Cholesky分解对协方差矩阵做相关性处理。我第一次做的时候图省事直接独立抽样结果场景削减后保留的那些“代表场景”全是锯齿状的调度结果惨不忍睹。后来老老实实加上时间相关性效果立刻就不一样了。蒙特卡洛抽样的数量我建议至少在500到1000个场景起步。场景太少覆盖不了真实的极端情况太多后续场景削减的计算量会明显上升。我自己常用的组合是初始生成1000个场景削减到10到20个。这个数量对MILP求解来说是完全可以接受的。3.2 场景削减不是简单聚类要保留概率信息场景削减听起来像聚类但目标不太一样。随机规划中的场景削减要在保证削减后的场景集合在概率分布意义下最接近原始场景集合的前提下尽量压缩场景数量。换句话说场景削减结束后留下的每个场景还得带一个权重这个权重就是它在新概率体系中的出现概率。主流的场景削减算法有快速前向选择法Fast Forward Selection、同步回代消除法Simultaneous Backward Reduction以及基于K-means聚类的简化做法。快速前向选择的原理是从初始场景集中逐次挑出K个最能代表整体分布的场景让原始场景与替代场景之间的Kantorovich距离最小。同步回代消除法反过来每次删除一个对分布影响最小的场景并把被删场景的概率叠加到离它最近的保留场景上。在Matlab里实现同步回代消除法并不复杂。核心循环是先算所有场景两两之间的距离矩阵然后找总删除代价最小的那个场景删掉同时更新概率。不过这个算法朴素实现有个毛病每一次删除都要重新计算距离和代价循环K次就要算K×N×N次距离当初始场景数达到1000以上时跑起来有点慢。我实际用的方法是先对原始场景做一次K-means聚类把1000个场景粗分为50个簇然后把每个簇的中心作为初始场景集合再在这个小集合上用同步回代消除法精准削减到15个场景。这样组合拳的效率比单纯用同步回代消除法高很多而且削减后的场景效果差距不大。如果你不想自己写距离计算也可以用Matlab的pdist函数直接算欧氏距离但要注意不同变量光伏、风电、负荷量纲不一样最好先做归一化再算距离否则数值大的变量会主导距离指标。还有一个操作细节场景削减时应该把光伏、风电、电负荷、热负荷的多维场景合并成一个向量一起削减而不是分别削减再随机组合。因为原始随机变量之间存在相关性分开削减会破坏这个相关性结构。3.3 场景质量的检验方法场景生成完别急着丢进优化模型。我会先画图做直观检查——把保留的场景曲线画在同一个坐标系里看看它们是否包络住了预测曲线是否覆盖了原始场景的上下界范围。这个目测很重要如果削减后的场景全挤在预测线附近说明极端情况都被削没了如果场景之间相差太大说明削减得不够。更严格一点的检验是计算场景削减前后的概率距离。用Kantorovich距离来衡量削减前后距离越小说明效果越好。还可以算削减后场景的均值与方差和原始样本的均值方差对比偏差在5%以内我就认为可以接受。这个检验代码很短十几行就能搞定属于性价比很高的步骤。4. Matlab代码架构与核心实现4.1 代码整体框架从数据到求解器这部分的代码实现思路我按模块拆分来讲。整个项目我习惯分成五个脚本模块数据定义模块、场景生成模块、场景削减模块、优化建模模块、结果分析与绘图模块。用脚本而不是一个巨大的函数调试时更方便哪一步出问题直接定位。数据定义模块里包括设备参数额定容量、效率、爬坡率、储能参数、预测出力曲线光伏、风电、电负荷、热负荷24小时数据、电价信息分时电价是分时还实时这直接影响调度结果。场景生成模块的伪代码思路如下%% 场景生成模块 ns 1000; % 初始场景数 T 24; % 调度周期 pv_err zeros(ns, T); % 光伏误差场景矩阵 % 假设光伏预测误差服从Beta分布经验参数 for t 1:T a 2; b 2; % 根据历史数据拟合 pv_err(:, t) betarnd(a, b, ns, 1) - a/(ab); end % 同理生成风电、电负荷、热负荷误差 % 对误差序列做AR(1)时间平滑 phi 0.8; % 相关性系数 for s 1:ns for t 2:T pv_err(s,t) phi * pv_err(s,t-1) (1-phi) * pv_err(s,t); end end % 叠加到预测值上 pv_scenarios pv_forecast pv_err * pv_capacity;4.2 场景削减核心代码示例同步回代消除法的核心逻辑可以浓缩成一个循环。我直接把可运行的简化代码放出来大家改改数据和参数就能跑。%% 场景削减同步回代消除法 % scenarios: ns维矩阵, 每行为一个场景向量 % ns: 初始场景数 % rs: 目标保留场景数 % probs: 初始概率, 一般为1/ns function [rs_scen, rs_prob] sceneReduction(scenarios, probs, rs) ns size(scenarios, 1); % 计算距离矩阵 D squareform(pdist(scenarios)); D(D0) inf; % 排除自身 while ns rs % 计算每个场景的删除代价: 到最近其他场景的距离 * 概率 cost_vec zeros(1, ns); for i 1:ns cost_vec(i) probs(i) * min(D(i, :)); end % 找到删除代价最小的场景 [~, del_idx] min(cost_vec); % 找到被删场景的最近场景 [~, nearest_idx] min(D(del_idx, :)); % 概率叠加到最近场景 probs(nearest_idx) probs(nearest_idx) probs(del_idx); % 删除场景和对应距离 scenarios(del_idx, :) []; probs(del_idx) []; D(del_idx, :) []; D(:, del_idx) []; ns ns - 1; end rs_scen scenarios; rs_prob probs; end这段代码的核心逻辑就是反复执行“找到最不重要的场景、删掉它、把它的概率转给离它最近的场景”这三个动作。注意每删一个场景都要重新算一次距离矩阵这是正确性优先的写法。如果初始场景数特别大比如5000以上建议先K-means粗削减。4.3 优化模型在Yalmip中的建模要点搭建优化模型我强烈建议用Yalmip工具箱。它能把几十行约束用接近数学公式的语法写出来省去手推矩阵的噩梦。贴上核心建模片段%% 优化模型构建 (Yalmip) % 决策变量 P_gt sdpvar(T, 1); % 燃气轮机出力 H_gt sdpvar(T, 1); % 燃气轮机余热回收量 P_gb sdpvar(T, 1); % 燃气锅炉出力 P_pv sdpvar(T, 1); % 光伏实际使用量 P_wt sdpvar(T, 1); % 风电实际使用量 P_grid sdpvar(T, 1); % 购电量正表示购电负表示售电 SOC_es sdpvar(T1, 1); % 蓄电池SOC P_ch sdpvar(T, 1); % 充电功率 P_dis sdpvar(T, 1); % 放电功率 on_gt binvar(T, 1); % 燃气轮机启停状态 % 目标函数各场景概率加权成本之和 F 0; for s 1:ns_keep % 燃料费用燃气轮机燃气锅炉 fuel_cost c_gas * sum(P_gt / eta_gt P_gb / eta_gb); % 购电费用 grid_cost sum(price .* P_grid); % 储能劣化费用 battery_cost c_battery * sum(P_ch P_dis); % 弃光弃风惩罚 curtail_cost c_curtail * sum((P_pv_avail(:,s) - P_pv) (P_wt_avail(:,s) - P_wt)); % 场景s对应的总运行成本 F_s fuel_cost grid_cost battery_cost curtail_cost; % 概率加权求和 F F probs_keep(s) * F_s; end注意这里有个很重要的建模细节判断决策变量是否跨场景复用。燃气轮机的出力、启停状态、购电计划这类“可调整决策”在随机规划中属于先行动作通常在未知场景实现前就要敲定所以跨场景复用同一个变量而光伏实际使用量、储能充放电量属于可调整变量可以在场景实现后再决定增减。这部分如果不做区分就建模随机优化的意义就大打折扣了。约束条件按之前说的四类添加到模型中。加上约束后用以下代码一行求解%% 求解 ops sdpsettings(solver, gurobi, verbose, 2, mipgap, 0.001); sol optimize(cons, F, ops);如果不装外部求解器也可以指定solver, cplex或者用Matlab自带的intlinprog。但要提醒一句模型规模上来后自带求解器的效率比Gurobi、Cplex差不少。用Gurobi求解这类带随机场景约束的MILP问题通常几秒到几十秒就能出结果intlinprog可能要跑好几分钟还不一定能找到高质量整数解。求解完成后记得检查sol.problem是否为0如果不是0该看约束可行性检查的排查步骤这部分放到第5节讲。4.4 结果绘图与小技巧优化结果出来后绘图是让方案可视化的重要一环。我一般画四张图第一张电功率平衡图。把燃气轮机、光伏、风电、购电、储能的充放电量堆叠成柱状图叠起来的高度应该正好等于电负荷。第二张热功率平衡图。展示余热回收、燃气锅炉、蓄热罐以及热负荷的平衡情况。第三张SOC曲线。看蓄电池和蓄热罐一天之内的能量变化确认SOC在安全区间内且始末一致。第四张不同场景下的运行成本对比。把消减前后的成本画出来直观体现随机优化相对确定性调度的优势。画图时有个小建议用面积图而不是多条折线叠加否则设备出力互相堆叠时很难看出总量关系。Matlab里用area函数画堆叠面积图非常方便比自己手动累加折线省事得多。5. 仿真调试中常见的坑与排查技巧5.1 求解器报“不可行”的排查套路我刚做这个项目时最头疼的就是solve返回 infeasible。这玩意一旦出现意味着你构造的约束条件根本不存在可行解。排查思路是有套路的按以下顺序检查大概率能解决问题第一步检查约束是否矛盾。最典型的矛盾是把电平衡等式写成了硬约束同时又把光伏出力上限设成等于预测值而预测值本身又低于必须满足的电负荷。这时候要不就是弃负荷、要不就得提高购电上限。看等式约束两边可调节变量的余量这是最快定位冲突的方式。第二步检查储能约束。始末SOC相等这个约束经常和储能容量上下限打架。比如蓄电池容量很小但始末必须都是80%一刻钟里又要支撑一个很大的负荷峰值那肯定不可行。解决办法是放宽末时刻SOC范围比如允许在60%到90%之间或者把储能容量参数调大。第三步检查设备爬坡约束。燃气轮机爬坡率设得太小加上负荷陡增时也会导致无解。把爬坡约束先注释掉看模型是否可解。如果注释后可行小步调试爬坡率就行了。第四步检查场景约束写法。如果在场景循环里误用了同一个变量名或者漏了下标也可能出现约束错乱。建议加一段约束数量校验代码把Yalmip自动生成的约束个数打印出来跟你预期的数量对比能快速发现漏加或者重复加了约束。5.2 随机优化结果反而“更差”我遇到过几次很困惑的情况随机优化的期望成本比确定性调度的成本还要高领导就说你这随机优化有意义吗其实这是正常现象。随机优化不是要让“期望成本最低”这件事绝对优于确定性方案而是在不确定性环境下给出一个稳健的、在多个可能场景下都能跑的方案。它牺牲了一部分在最优场景下的经济性换来的是在最恶劣场景下依然可行且成本可控。具体来说随机优化的成本通常比“事后诸葛亮”的确定性调度高但比“预测错了硬扛”的确定性调度低很多。所以对比的时候一定要拉开差距一个是按实际出力场景回放确定性优化得到的决策结果一个是随机优化得到的决策结果然后在同样一批真实场景上回放比较这样才有意义。有的论文里管这叫“模拟回放检验”这是评价随机优化价值的金标准。5.3 场景削减后结果波动大是为什么场景削减数量太少比如只留3到5个会导致随机优化的结果对场景集合特别敏感换一批场景结果就差很大。我的经验是热电联供微网这种场景保留场景数至少要在10个以上最好在15到20个之间。太少随机分布的细节捕捉不够太多求解时间指数级上升收益递减。另外场景削减前一定要做归一化。我踩过最深的坑就在这里。原始场景向量里电负荷的数值在几百千瓦量级而光伏在几十千瓦量级pdist算距离时电负荷一方主导了整个距离矩阵导致削减出来的场景全部在优化电负荷的拟合光伏的多样性直接丢失。后来把所有变量先除以其最大值做归一化再削减问题迎刃而解。5.4 Matlab求解时间太长怎么办MILP求解时间过长最常见的瓶颈是整数变量太多以及场景数太多导致约束规模膨胀。我的优化策略有这几个第一把不必要的整数变量实数化。比如储能充放电状态可以用两个连续变量充和放加上一个很小的惩罚项代替二进制变量只要惩罚项不影响最优解太多就好。第二矩阵化约束写写避免循环里逐条addConstraintYalmip虽然方便但循环10万次加约束也会很慢。尽量用向量化表达式一次性构造约束。第三合理设MIP Gap。我一般设成0.001也就是允许千分之一的误差这能大幅加速收敛。如果只是做趋势分析设0.01也够用求解速度能快好几倍。5.5 设备参数设置不合理导致的“假最优”还有一个隐蔽陷阱当光伏渗透率很高、弃光惩罚设置过低时优化结果会倾向于少用光伏、多用燃气轮机因为边际成本低还在燃气轮机出力范围内。看起来结果很漂亮成本也不高但实际上没把问题体现在结果里。我的做法是把弃光惩罚设置为燃气轮机发电边际成本的1.2到1.5倍这样优化器才会优先消纳光伏而不是“花钱买省事”。6. 一点个人经验与后续扩展思路跑完整个项目我最大的感受是随机优化这部分风险不在数学模型的推导而在“从公式到Matlab代码”这段路上。公式里的场景概率、约束下标在代码里稍有偏差就导致结果完全不可信。所以强烈建议每写一个模块就用一组小规模测试数据验证一下——单独测试场景生成模块、单独测试场景削减模块、单独测试不含随机场景的确定性模型确认确定性模型跑得没问题再把场景迭进去。这样做定位bug的时间能缩短一半以上。后
做热电联供微网优化有一段时间了从最初的确定性调度一步步做到考虑源荷随机特征的随机优化踩了不少坑也积累了一些能直接用的经验。这篇就围绕“考虑源荷随机特征的热电联供微网优化研究Matlab代码实现”这个方向把整个建模到代码落地的链路完整拆开讲一遍。适合正在做微网优化调度、综合能源系统方向课题的学生也适合刚入行做园区能源规划、想了解随机优化怎么落地的工程师。内容会涉及模型思路、数学公式怎么转成代码、场景生成与削减的实操技巧以及我自己在跑仿真时遇到的一堆典型问题。1. 热电联供微网优化到底在解决什么问题1.1 先搞清楚系统的“源”和“荷”长什么样热电联供微网简单说就是一个麻雀虽小五脏俱全的小型能源系统。它通常包含发电设备微型燃气轮机、内燃机、供热设备燃气锅炉、余热锅炉、可再生能源光伏、风电、储能装置蓄电池、蓄热罐以及电负荷和热负荷两类用户需求。这个系统的核心价值就在于用一套设备同时满足电和热两种需求实现能源的梯级利用。很多人一开始容易把“优化”想复杂。其实本质上就是回答一个问题在未来24小时里每一台设备在每一个时刻该发多少电、产多少热、充放多少储能才能让总运行成本最低如果光伏、风电出力是确定已知的负荷也是固定数值这是个标准的混合整数线性规划MILP问题直接丢给求解器就能出结果。但现实没那么理想。光伏靠天吃饭云一遮就断崖式下跌风电更是看风速脸色用户的电热负荷也受天气、作息、生产计划影响波动很明显。如果调度方案完全基于预测的确定性数值来制定实际运行时很容易出现功率缺额或者能源浪费。这就引出了标题里的核心关键词——源荷随机特征。1.2 随机特征怎么进入优化框架处理源荷随机特征工程和学术上常用的路线有三条。第一条是随机规划Stochastic Programming用多个离散场景描述可能的出力与负荷情况让优化结果在所有场景下综合表现最优。第二条是鲁棒优化考虑最坏情况追求“不管怎么波动都不会出事”。第三条是分布鲁棒优化介于两者之间考虑概率分布的不确定性。我这次实现用的是随机规划中的场景法也是目前最直观、最容易在Matlab里落地的一种。核心思路分三步先对光伏出力、风电出力、电负荷、热负荷的预测误差分别建模用蒙特卡洛抽样生成大量场景然后用场景削减技术去掉冗余场景保留少数几个有代表性的场景最后把这几个场景一起塞进优化模型目标函数变成“所有场景下运行成本的概率加权求和”约束条件则要保证每个场景下系统都能安全运行。这样做的好处很明显计算复杂度可控求解器处理起来不费劲同时能让调度决策在随机环境下保持可行性不为单一预测值“赌运气”。缺点也有就是结果对场景质量比较敏感。场景生成得不好优化结果就会偏离实际这一块后面会详细讲。2. 目标函数与约束条件的建模细节2.1 目标函数不是只算发电成本做热电联供微网优化第一步就是把目标函数写清楚。注意这里的“成本”不是单一的电费而是一整套运行费用的总和。我常用的目标函数形式如下min F Σs ωs × Σt [ CF(t) CM(t) CB(t) CG(t) CP(t) ]其中s是场景索引ωs是场景s的概率权重CF(t)是t时刻的燃料成本包括燃气轮机的天然气消耗和燃气锅炉的天然气消耗CM(t)是设备维护成本通常按出力线性折算CB(t)是储能设备劣化成本用于避免蓄电池和蓄热罐频繁充放CG(t)是向大电网购电的成本如果允许向电网售电还需要扣掉售电收益CP(t)是弃风弃光惩罚项用来保证优化不会把可再生出力全砍掉来省成本。很多初学者忽略了一个细节热电联产机组的发电和发热不是独立的两者之间存在强耦合关系。同一个机组发电出力定了可回收的余热量就有个上限不是想发多少热就发多少热。这个耦合关系通常用一个“热电比”系数来描述具体实现时要查机组的技术参数手册或者用厂家提供的性能曲线拟合出可行域不等式。还有一个点容易被忽略——机组启停状态。燃气轮机频繁启停不仅损耗设备寿命还会抬升成本。所以模型里要引入二进制变量表示机组启停状态并且配合最小启停时间约束。这就把模型从LP变成了MILP求解时间会明显增加。后面在调试时MILP求解时间是大家问得最多的问题之一。2.2 约束条件的完整清单与物理含义约束条件是优化模型的骨架每一条都对应着实际的物理限制。我在写代码时会把约束分成四大类方便检查和调试。第一类是功率平衡约束。电功率平衡要求任意场景、任意时刻燃气轮机发电、光伏出力、风电出力、蓄电池放电、电网购电之和恰好等于电负荷加上蓄电池充电和电转热设备如果有的话的耗电。热功率平衡要求燃气轮机余热、燃气锅炉产热、蓄热罐放热之和等于热负荷加上蓄热罐充热。要注意电平衡里电网交互功率是变量既可以是正值购电也可以是负值售电很多论文会用正负两个变量来线性化处理避免目标函数里出现绝对值项。第二类是设备自身约束。比如燃气轮机出力的上下限约束、爬坡约束相邻时刻出力变化不能超过额定值的一定比例、最小启停时间约束。光伏和风电的出力上限就是该时刻场景提供的最大可用功率这个值随场景和时刻变化。燃气锅炉同样有出力上下限和爬坡约束。蓄电池和蓄热罐需要重点处理——不仅要限制充放电功率上限还要保证存储能量的动态递推关系当前时刻能量等于上一时刻能量加上充入量减去放出量并考虑充放效率损失。第三类是储能状态约束。储能容量要在安全范围内不能过充过放。我一般还会加一个调度周期始末能量一致的约束也就是日始日终的SOC荷电状态要相等保证调度方案可以循环执行。这个约束在带储能的微网优化里几乎是标配。第四类是随机性相关约束。在场景法中每个场景都要满足前述所有平衡和运行约束。也就是说不管哪个场景发生了系统都必须能安全运行。这就是所谓的“非预期性约束”和“场景可行性约束”。有的论文还会加机会约束比如“热负荷供应不足的概率不超过5%”这种形式在Matlab里实现会更复杂通常需要做转化处理这次先不展开。把上述约束全部堆进模型后整个优化问题的规模会达到几千个变量、几千条约束。对MatlabYalmip求解器来说不算大事但你得保证建模型的时候矩阵别搞错不然一跑起来全是不可行解。3. 源荷随机场景的生成与削减实操3.1 用蒙特卡洛模拟生成原始场景随机优化效果好不好场景生成是第一步。我的做法是先对光伏、风电、负荷的预测值加上带统计特征的误差项。误差的分布不是随便拍脑袋定的得基于历史数据和天气预报水平。光伏误差一般用Beta分布描述风电误差常用Weibull分布负荷预测误差在很多文献里被近似为均值为零的正态分布也有人用t分布来捕捉厚尾特征。具体操作时我一般先对每个时段独立生成误差样本按以下方式叠加到预测值上P_scenario(t, s) P_forecast(t) ε(t, s)这里 ε(t, s) 就是从误差分布中抽出来的随机扰动。但这里有个隐患直接逐时段独立抽样会让生成的场景波动非常剧烈相邻时段忽高忽低完全不符合实际出力曲线的时间连续性。解决这个问题可以在误差抽样时引入时间相关性比如用一阶自回归模型AR(1)来生成平滑的误差序列或者用Cholesky分解对协方差矩阵做相关性处理。我第一次做的时候图省事直接独立抽样结果场景削减后保留的那些“代表场景”全是锯齿状的调度结果惨不忍睹。后来老老实实加上时间相关性效果立刻就不一样了。蒙特卡洛抽样的数量我建议至少在500到1000个场景起步。场景太少覆盖不了真实的极端情况太多后续场景削减的计算量会明显上升。我自己常用的组合是初始生成1000个场景削减到10到20个。这个数量对MILP求解来说是完全可以接受的。3.2 场景削减不是简单聚类要保留概率信息场景削减听起来像聚类但目标不太一样。随机规划中的场景削减要在保证削减后的场景集合在概率分布意义下最接近原始场景集合的前提下尽量压缩场景数量。换句话说场景削减结束后留下的每个场景还得带一个权重这个权重就是它在新概率体系中的出现概率。主流的场景削减算法有快速前向选择法Fast Forward Selection、同步回代消除法Simultaneous Backward Reduction以及基于K-means聚类的简化做法。快速前向选择的原理是从初始场景集中逐次挑出K个最能代表整体分布的场景让原始场景与替代场景之间的Kantorovich距离最小。同步回代消除法反过来每次删除一个对分布影响最小的场景并把被删场景的概率叠加到离它最近的保留场景上。在Matlab里实现同步回代消除法并不复杂。核心循环是先算所有场景两两之间的距离矩阵然后找总删除代价最小的那个场景删掉同时更新概率。不过这个算法朴素实现有个毛病每一次删除都要重新计算距离和代价循环K次就要算K×N×N次距离当初始场景数达到1000以上时跑起来有点慢。我实际用的方法是先对原始场景做一次K-means聚类把1000个场景粗分为50个簇然后把每个簇的中心作为初始场景集合再在这个小集合上用同步回代消除法精准削减到15个场景。这样组合拳的效率比单纯用同步回代消除法高很多而且削减后的场景效果差距不大。如果你不想自己写距离计算也可以用Matlab的pdist函数直接算欧氏距离但要注意不同变量光伏、风电、负荷量纲不一样最好先做归一化再算距离否则数值大的变量会主导距离指标。还有一个操作细节场景削减时应该把光伏、风电、电负荷、热负荷的多维场景合并成一个向量一起削减而不是分别削减再随机组合。因为原始随机变量之间存在相关性分开削减会破坏这个相关性结构。3.3 场景质量的检验方法场景生成完别急着丢进优化模型。我会先画图做直观检查——把保留的场景曲线画在同一个坐标系里看看它们是否包络住了预测曲线是否覆盖了原始场景的上下界范围。这个目测很重要如果削减后的场景全挤在预测线附近说明极端情况都被削没了如果场景之间相差太大说明削减得不够。更严格一点的检验是计算场景削减前后的概率距离。用Kantorovich距离来衡量削减前后距离越小说明效果越好。还可以算削减后场景的均值与方差和原始样本的均值方差对比偏差在5%以内我就认为可以接受。这个检验代码很短十几行就能搞定属于性价比很高的步骤。4. Matlab代码架构与核心实现4.1 代码整体框架从数据到求解器这部分的代码实现思路我按模块拆分来讲。整个项目我习惯分成五个脚本模块数据定义模块、场景生成模块、场景削减模块、优化建模模块、结果分析与绘图模块。用脚本而不是一个巨大的函数调试时更方便哪一步出问题直接定位。数据定义模块里包括设备参数额定容量、效率、爬坡率、储能参数、预测出力曲线光伏、风电、电负荷、热负荷24小时数据、电价信息分时电价是分时还实时这直接影响调度结果。场景生成模块的伪代码思路如下%% 场景生成模块 ns 1000; % 初始场景数 T 24; % 调度周期 pv_err zeros(ns, T); % 光伏误差场景矩阵 % 假设光伏预测误差服从Beta分布经验参数 for t 1:T a 2; b 2; % 根据历史数据拟合 pv_err(:, t) betarnd(a, b, ns, 1) - a/(ab); end % 同理生成风电、电负荷、热负荷误差 % 对误差序列做AR(1)时间平滑 phi 0.8; % 相关性系数 for s 1:ns for t 2:T pv_err(s,t) phi * pv_err(s,t-1) (1-phi) * pv_err(s,t); end end % 叠加到预测值上 pv_scenarios pv_forecast pv_err * pv_capacity;4.2 场景削减核心代码示例同步回代消除法的核心逻辑可以浓缩成一个循环。我直接把可运行的简化代码放出来大家改改数据和参数就能跑。%% 场景削减同步回代消除法 % scenarios: ns维矩阵, 每行为一个场景向量 % ns: 初始场景数 % rs: 目标保留场景数 % probs: 初始概率, 一般为1/ns function [rs_scen, rs_prob] sceneReduction(scenarios, probs, rs) ns size(scenarios, 1); % 计算距离矩阵 D squareform(pdist(scenarios)); D(D0) inf; % 排除自身 while ns rs % 计算每个场景的删除代价: 到最近其他场景的距离 * 概率 cost_vec zeros(1, ns); for i 1:ns cost_vec(i) probs(i) * min(D(i, :)); end % 找到删除代价最小的场景 [~, del_idx] min(cost_vec); % 找到被删场景的最近场景 [~, nearest_idx] min(D(del_idx, :)); % 概率叠加到最近场景 probs(nearest_idx) probs(nearest_idx) probs(del_idx); % 删除场景和对应距离 scenarios(del_idx, :) []; probs(del_idx) []; D(del_idx, :) []; D(:, del_idx) []; ns ns - 1; end rs_scen scenarios; rs_prob probs; end这段代码的核心逻辑就是反复执行“找到最不重要的场景、删掉它、把它的概率转给离它最近的场景”这三个动作。注意每删一个场景都要重新算一次距离矩阵这是正确性优先的写法。如果初始场景数特别大比如5000以上建议先K-means粗削减。4.3 优化模型在Yalmip中的建模要点搭建优化模型我强烈建议用Yalmip工具箱。它能把几十行约束用接近数学公式的语法写出来省去手推矩阵的噩梦。贴上核心建模片段%% 优化模型构建 (Yalmip) % 决策变量 P_gt sdpvar(T, 1); % 燃气轮机出力 H_gt sdpvar(T, 1); % 燃气轮机余热回收量 P_gb sdpvar(T, 1); % 燃气锅炉出力 P_pv sdpvar(T, 1); % 光伏实际使用量 P_wt sdpvar(T, 1); % 风电实际使用量 P_grid sdpvar(T, 1); % 购电量正表示购电负表示售电 SOC_es sdpvar(T1, 1); % 蓄电池SOC P_ch sdpvar(T, 1); % 充电功率 P_dis sdpvar(T, 1); % 放电功率 on_gt binvar(T, 1); % 燃气轮机启停状态 % 目标函数各场景概率加权成本之和 F 0; for s 1:ns_keep % 燃料费用燃气轮机燃气锅炉 fuel_cost c_gas * sum(P_gt / eta_gt P_gb / eta_gb); % 购电费用 grid_cost sum(price .* P_grid); % 储能劣化费用 battery_cost c_battery * sum(P_ch P_dis); % 弃光弃风惩罚 curtail_cost c_curtail * sum((P_pv_avail(:,s) - P_pv) (P_wt_avail(:,s) - P_wt)); % 场景s对应的总运行成本 F_s fuel_cost grid_cost battery_cost curtail_cost; % 概率加权求和 F F probs_keep(s) * F_s; end注意这里有个很重要的建模细节判断决策变量是否跨场景复用。燃气轮机的出力、启停状态、购电计划这类“可调整决策”在随机规划中属于先行动作通常在未知场景实现前就要敲定所以跨场景复用同一个变量而光伏实际使用量、储能充放电量属于可调整变量可以在场景实现后再决定增减。这部分如果不做区分就建模随机优化的意义就大打折扣了。约束条件按之前说的四类添加到模型中。加上约束后用以下代码一行求解%% 求解 ops sdpsettings(solver, gurobi, verbose, 2, mipgap, 0.001); sol optimize(cons, F, ops);如果不装外部求解器也可以指定solver, cplex或者用Matlab自带的intlinprog。但要提醒一句模型规模上来后自带求解器的效率比Gurobi、Cplex差不少。用Gurobi求解这类带随机场景约束的MILP问题通常几秒到几十秒就能出结果intlinprog可能要跑好几分钟还不一定能找到高质量整数解。求解完成后记得检查sol.problem是否为0如果不是0该看约束可行性检查的排查步骤这部分放到第5节讲。4.4 结果绘图与小技巧优化结果出来后绘图是让方案可视化的重要一环。我一般画四张图第一张电功率平衡图。把燃气轮机、光伏、风电、购电、储能的充放电量堆叠成柱状图叠起来的高度应该正好等于电负荷。第二张热功率平衡图。展示余热回收、燃气锅炉、蓄热罐以及热负荷的平衡情况。第三张SOC曲线。看蓄电池和蓄热罐一天之内的能量变化确认SOC在安全区间内且始末一致。第四张不同场景下的运行成本对比。把消减前后的成本画出来直观体现随机优化相对确定性调度的优势。画图时有个小建议用面积图而不是多条折线叠加否则设备出力互相堆叠时很难看出总量关系。Matlab里用area函数画堆叠面积图非常方便比自己手动累加折线省事得多。5. 仿真调试中常见的坑与排查技巧5.1 求解器报“不可行”的排查套路我刚做这个项目时最头疼的就是solve返回 infeasible。这玩意一旦出现意味着你构造的约束条件根本不存在可行解。排查思路是有套路的按以下顺序检查大概率能解决问题第一步检查约束是否矛盾。最典型的矛盾是把电平衡等式写成了硬约束同时又把光伏出力上限设成等于预测值而预测值本身又低于必须满足的电负荷。这时候要不就是弃负荷、要不就得提高购电上限。看等式约束两边可调节变量的余量这是最快定位冲突的方式。第二步检查储能约束。始末SOC相等这个约束经常和储能容量上下限打架。比如蓄电池容量很小但始末必须都是80%一刻钟里又要支撑一个很大的负荷峰值那肯定不可行。解决办法是放宽末时刻SOC范围比如允许在60%到90%之间或者把储能容量参数调大。第三步检查设备爬坡约束。燃气轮机爬坡率设得太小加上负荷陡增时也会导致无解。把爬坡约束先注释掉看模型是否可解。如果注释后可行小步调试爬坡率就行了。第四步检查场景约束写法。如果在场景循环里误用了同一个变量名或者漏了下标也可能出现约束错乱。建议加一段约束数量校验代码把Yalmip自动生成的约束个数打印出来跟你预期的数量对比能快速发现漏加或者重复加了约束。5.2 随机优化结果反而“更差”我遇到过几次很困惑的情况随机优化的期望成本比确定性调度的成本还要高领导就说你这随机优化有意义吗其实这是正常现象。随机优化不是要让“期望成本最低”这件事绝对优于确定性方案而是在不确定性环境下给出一个稳健的、在多个可能场景下都能跑的方案。它牺牲了一部分在最优场景下的经济性换来的是在最恶劣场景下依然可行且成本可控。具体来说随机优化的成本通常比“事后诸葛亮”的确定性调度高但比“预测错了硬扛”的确定性调度低很多。所以对比的时候一定要拉开差距一个是按实际出力场景回放确定性优化得到的决策结果一个是随机优化得到的决策结果然后在同样一批真实场景上回放比较这样才有意义。有的论文里管这叫“模拟回放检验”这是评价随机优化价值的金标准。5.3 场景削减后结果波动大是为什么场景削减数量太少比如只留3到5个会导致随机优化的结果对场景集合特别敏感换一批场景结果就差很大。我的经验是热电联供微网这种场景保留场景数至少要在10个以上最好在15到20个之间。太少随机分布的细节捕捉不够太多求解时间指数级上升收益递减。另外场景削减前一定要做归一化。我踩过最深的坑就在这里。原始场景向量里电负荷的数值在几百千瓦量级而光伏在几十千瓦量级pdist算距离时电负荷一方主导了整个距离矩阵导致削减出来的场景全部在优化电负荷的拟合光伏的多样性直接丢失。后来把所有变量先除以其最大值做归一化再削减问题迎刃而解。5.4 Matlab求解时间太长怎么办MILP求解时间过长最常见的瓶颈是整数变量太多以及场景数太多导致约束规模膨胀。我的优化策略有这几个第一把不必要的整数变量实数化。比如储能充放电状态可以用两个连续变量充和放加上一个很小的惩罚项代替二进制变量只要惩罚项不影响最优解太多就好。第二矩阵化约束写写避免循环里逐条addConstraintYalmip虽然方便但循环10万次加约束也会很慢。尽量用向量化表达式一次性构造约束。第三合理设MIP Gap。我一般设成0.001也就是允许千分之一的误差这能大幅加速收敛。如果只是做趋势分析设0.01也够用求解速度能快好几倍。5.5 设备参数设置不合理导致的“假最优”还有一个隐蔽陷阱当光伏渗透率很高、弃光惩罚设置过低时优化结果会倾向于少用光伏、多用燃气轮机因为边际成本低还在燃气轮机出力范围内。看起来结果很漂亮成本也不高但实际上没把问题体现在结果里。我的做法是把弃光惩罚设置为燃气轮机发电边际成本的1.2到1.5倍这样优化器才会优先消纳光伏而不是“花钱买省事”。6. 一点个人经验与后续扩展思路跑完整个项目我最大的感受是随机优化这部分风险不在数学模型的推导而在“从公式到Matlab代码”这段路上。公式里的场景概率、约束下标在代码里稍有偏差就导致结果完全不可信。所以强烈建议每写一个模块就用一组小规模测试数据验证一下——单独测试场景生成模块、单独测试场景削减模块、单独测试不含随机场景的确定性模型确认确定性模型跑得没问题再把场景迭进去。这样做定位bug的时间能缩短一半以上。后
