做微网优化调度的朋友应该都有这种体会纯电系统还好一旦带上热负荷模型复杂度直接翻倍。要是在这个基础上再考虑源荷不确定性很多人的第一反应就是“先按期望值算个确定性模型算了”。但这样算出来的调度方案在实际运行中往往会被风电出力的突然下跌或者热负荷的意外爬升打得措手不及。今天要拆解的这套“源荷不确定性下考虑随机特征的热电联供微网优化模型及Matlab程序实现”要解决的就是这个问题——同时考虑电、热两种负荷和风光出力的随机波动用场景法把不确定性引入优化模型最终给出一个经济性最优、且能应对一定随机波动的调度方案。这篇文章我会从建模思路讲到Matlab代码实现全程避坑适合正在做微网优化、综合能源系统方向的研究生以及刚接触随机优化的工程师参考。1. 模型整体设计与优化思路拆解1.1 为什么要考虑源荷不确定性先想清楚一个事确定性模型和随机优化模型之间的差别到底在哪。传统的确定性调度是把风电出力、光伏出力、电负荷、热负荷都当成已知参数输入一个固定的值然后求解出一个确定的机组启停和出力计划。它的问题是——风电光伏的预测误差通常在10%到30%之间负荷预测也有一定的偏差这些误差叠加起来足够让调度方案在实际执行的时候出现功率不平衡或者备用不足。而随机优化模型的核心思路是把这些不确定参数用一组带概率的随机场景来描述。每个场景就是一组可能的源荷实现值比如“风电低出力、负荷高峰”算一个场景“风电高出力、负荷低谷”算另一个场景。优化模型在决策的时候让最终方案在所有场景下都能满足运行约束同时最小化期望成本。这样做出来的方案不再是最优的一个点而是在概率意义上整体最优的一个区间策略。1.2 整体技术路线场景生成-场景削减-优化求解这套模型的技术路线可以用三步来概括场景生成、场景削减和优化求解。场景生成就是根据源荷预测值和误差分布抽样生成大量可能的源荷场景。这一步最常用的是蒙特卡洛抽样也可以做拉丁超立方采样后者在小样本下分布覆盖更均匀。场景削减则是把规模庞大的初始场景压缩成少量代表性场景避免优化模型因为场景数太多而解不动。最后用YALMIP工具箱建模调用Cplex或者Gurobi求解混合整数线性规划问题。在建模层我比较推荐用随机期望值模型Expected Value Model的框架目标函数是各个场景下运行成本的期望值最小化约束条件细分成每个场景都要满足的“场景约束”和与场景无关的“第一阶段决策约束”。这样写的好处是逻辑清晰后面用Matlab代码实现的时候约束矩阵的组装也方便。2. 源荷不确定性建模方法与随机特征刻画2.1 风电、光伏和负荷的概率分布选择不确定性的随机特征主要体现在三块风电出力、光伏出力和负荷波动。这三类不确定性的统计特性其实差得挺远建模方式也要分别处理。风电出力预测误差通常用正态分布近似。简单做法是认为实际出力等于预测值加上一个零均值、标准差为预测值一定比例的正态扰动项。更精细一点的做法是用Weibull分布描述风速再通过风机功率曲线转换为出力分布但这个过程中涉及非线性变换场景生成后分布形态会偏移工程上用得不多。光伏出力受光照强度影响理论上是Beta分布。标准做法是对光照强度做归一化然后用Beta分布拟合。但由于Beta分布抽样需要先估计形状参数alpha和beta操作起来比正态分布麻烦所以很多论文也简化为正态分布来处理。实际上如果样本数据足够直接用历史出力数据的经验分布做自助法Bootstrap抽样效果往往比强行套理论分布更好。电负荷和热负荷的不确定性一般也假设为预测值加正态扰动。不过要注意热负荷和电负荷往往有一定的相关性——室外温度升高时电负荷因为空调制冷上升热负荷反而下降。如果想建模得更精细可以考虑用Copula或者协方差矩阵来刻画这种相关性。在做毕设或者论文时相关性这一块做到位评审老师会认为你的模型“考虑随机特征”不只是噱头。2.2 蒙特卡洛抽样与拉丁超立方采样对比场景生成最常用的两个方法是蒙特卡洛抽样和拉丁超立方采样两者我都试过说下感受。蒙特卡洛抽样实现极其简单几行代码就能搞定就是对着每个不确定变量的分布函数直接randn或者random函数抽样。但它的缺点是会有聚类效应——抽样次数不够多的时候某些区域点很密某些区域点很稀疏场景的代表性不稳。拉丁超立方采样LHS的思路是把每个变量的分布区间均匀分层每层保证抽到一个样本点这样抽出来的场景分布更均匀。对于随机优化这种需要靠有限场景代表真实分布的场景LHS在同样抽样次数下通常比蒙特卡洛更稳。实际写Matlab代码时可以用lhsdesign函数或者手动分层抽样。用lhsdesign要注意它默认生成的是[0,1]区间的均匀分层样本需要再用icdf函数转成目标分布的分位数操作起来也不复杂。2.3 场景削减同步回代消除法的Matlab实现场景生成一万个很容易但优化模型里一万个场景根本解不动——变量规模直接膨胀一万倍。所以场景削减是必须的。削减的原则是去掉冗余场景保留的少量场景在概率分布上尽量接近原场景集合。最经典的方法是同步回代消除法Fast Forward Selection / Fast Backward Reduction。核心思想是贪心地、迭代地去找一个与剩余场景集合“距离最近”的场景用它来代替某几个被删除的场景同时累加概率距离直到场景数削减到预设值。Matlab实现这段逻辑其实不复杂核心步骤如下function [scen_keep, prob_keep] backwardReduction(scen, prob, keepNum) % 输入: scen - Ns x nVar 的原始场景矩阵 % prob - Ns x 1 的原始概率向量 % keepNum - 保留场景个数 Ns size(scen, 1); active true(Ns, 1); distMat zeros(Ns, Ns); for i 1:Ns for j 1:Ns if i ~ j distMat(i, j) norm(scen(i, :) - scen(j, :), 2); end end end while sum(active) keepNum minDistForEach inf(Ns, 1); for i 1:Ns if active(i) distToActive distMat(i, active); distToActive(i) []; % 排除自身 minDistForEach(i) min(distToActive); end end % 找删除代价最小的场景 [~, delIdx] min(prob .* active .* minDistForEach); % 找与之最近的活跃场景 distVec distMat(delIdx, :); distVec(~active) inf; distVec(delIdx) inf; [~, nearestIdx] min(distVec); % 概率转移 prob(nearestIdx) prob(nearestIdx) prob(delIdx); active(delIdx) false; end scen_keep scen(active, :); prob_keep prob(active); prob_keep prob_keep / sum(prob_keep); end这里有几个坑提醒一下。第一个是距离矩阵场景维度多、场景数又大时双层循环非常慢可以用pdist2函数一次性计算距离矩阵来加速。第二个是删除准则的权重经典算法是用“概率乘以到最近场景的距离”也就是说概率大但孤立的场景更不该删因为删了它代价大。第三削减完成之后被删除场景的概率要转移到离它最近的那个保留场景上这样总概率仍为1且保留场景的概率代表了原始场景集合的概率密度分布。3. 热电联供微网各设备建模与约束解析3.1 热电联产机组电热运行区间与可行域热电联产机组是这个微网的核心供能设备。它的特点就是电和热的生产是耦合的——发多少电就带出多少热。这个电热耦合关系决定了整个模型的可行域形状。燃气轮机型号不同运行区间也不同。常见的是背压式机组它的电热关系近似线性热出力等于电出力乘以热电比。还有抽凝式机组灵活性更高热出力可以在一个区间内调节对应的可行域是一个类似平行四边形的区域。用数学语言描述典型的CHP机组建模有一套不等式约束P_chp_min P_chp(t) P_chp_max H_chp(t) eta_recovery * P_chp(t) % 背压式简化模型这里P_chp是电出力H_chp是热出力eta_recovery是余热回收系数。对于抽凝式机组约束变为P_chp(t) - P_chp_min 0 P_chp_max - P_chp(t) 0 H_chp(t) c1 * P_chp(t) c2 H_chp(t) c3 * P_chp(t) c4后一组不等式围成的区域就形成了电热功率可行域。建这个可行域的目的是让优化器在寻优时知道电出力和热出力不能随便乱取必须落在机组物理特性允许的范围内。3.2 储能设备建模电储能与蓄热罐储能是微网灵活性的关键来源。电储能方面我一般用能量型模型把荷电状态SOC作为连续变量。SOC(t1) SOC(t) (P_ch_eff * P_ch(t) - P_dis(t) / P_dis_eff) * dt SOC_min SOC(t) SOC_max 0 P_ch(t) P_ch_max * u_ch(t) 0 P_dis(t) P_dis_max * u_dis(t) u_ch(t) u_dis(t) 1最后一条约束是防止同时充电和放电物理上这没有意义而且会让目标函数寻优钻空子。这里用了一个二进制变量u来保证储能不能同时充放。充放电效率两个参数值通常在0.95左右但注意要在模型中分开乘在充电功率和放电功率上别搞混。蓄热罐的建模逻辑和电储能几乎一样只是把SOC换成储热量充放电功率换成吸放热功率。由于热惯性的存在蓄热罐的时间常数往往比电储能大所以调度周期内可以考虑更长的状态保持时间。温度分层蓄热罐需要更复杂的模型但线性优化框架下简化成“能量罐”模型是最实用的做法。3.3 电网交互与购售电约束微网和外网之间有一个联络线功率交换有上下限这个约束一定要加否则优化器会可能让你完全不发电全部从电网买电这样的方案虽然理论上“最省”但不满足微网作为自治系统的建模初衷。切网模式下不涉及但如果并网运行就需要考虑。如果模型的调度周期是24小时、时间步长1小时那购电价格和售电价格会随峰谷时段变化。经典做法是给定分时电价向量目标函数中购电成本按购电价格计算售电收益按售电价格计算。注意售电价格通常低于购电价格这意味着你要在模型里用两个非负变量分别表示购电量和售电量而不是用一个有正有负的净交换功率变量。3.4 目标函数与约束全集构建目标函数是典型的经济调度目标最小化总运行成本期望值min sum_scenario prob(s) * [ sum_t (C_buy * P_buy(s,t) - C_sell * P_sell(s,t)) sum_t C_gas * F_gas(s,t) ]其中F_gas是燃气轮机消耗的燃料量C_gas是燃料单价C_buy和C_sell对应分时购售电价。储能和蓄热罐的运行成本可以设置一个很小的单位充放损耗系数主要目的是防止设备频繁无意义动作数值上设成0.001这样的小量就行。约束全集包括电功率平衡约束、热功率平衡约束、CHP机组运行约束、储能设备动态约束、储热罐动态约束、购售电功率上下限约束、爬坡约束。每一类约束都必须在每个场景、每个时段下成立。整个模型最终可以整理成混合整数线性规划MILP用求解器可以直接求解。4. Matlab程序实现从建模到求解的完整流程4.1 求解环境搭建与工具选择程序实现的第一步是准备求解环境。我用的组合是Matlab YALMIP Cplex或Gurobi。YALMIP是一个Matlab上的建模工具箱它最大的价值是你不用手动组装矩阵直接用符号化的约束表达式写模型特别适合做教学和原型验证。Cplex和Gurobi是商业求解器MILP求解效率远高于Matlab内置的intlinprog。如果你只有Matlab内置求解器小规模算例可以用intlinprog跑通但场景数一旦上百求解时间会让人崩溃。Cplex安装时有个小坑Cplex自带的Matlab接口目录要添加到Matlab路径中否则YALMIP调用时会报“找不到cplex”的错误。Gurobi也一样要配置环境变量并用gurobi_setup把接口加进Matlab路径。% 在Matlab中添加YALMIP路径 addpath(genpath(D:\Program Files\MATLAB\yalmip-master)); savepath; % 添加Cplex路径 addpath(C:\Program Files\IBM\ILOG\CPLEX_Studio_Community\cplex\matlab\x64_win64); savepath;路径换成你自己机器上的安装位置即可重点是savepath否则下次启动又要重新添加。4.2 数据结构设计与主程序框架程序的数据结构直接决定后面写约束的时候顺不顺畅。我的习惯是建一个结构体Params统一存放所有参数运行时段T24场景数Ns等作为全局控制变量。主程序流程是读取/生成参数 - 生成不确定场景 - 场景削减 - 定义决策变量 - 编写目标函数与约束 - 求解 - 结果可视化。决策变量的定义用YALMIP比较方便P_chp sdpvar(T, Ns, full); % CHP电出力24时段 x Ns场景 H_chp sdpvar(T, Ns, full); % CHP热出力 SOC sdpvar(T1, Ns, full); % 电储能SOC注意T1因为SOC(1)是初始值 P_buy sdpvar(T, Ns, full); % 从电网购电功率 P_sell sdpvar(T, Ns, full); % 向电网售电功率 u_chp binvar(T, Ns); % CHP启停变量如需要 u_ch binvar(T, Ns); % 储能充电状态 u_dis binvar(T, Ns); % 储能放电状态这里有一个关键点第一类决策变量机组启停状态通常是与场景无关的也就是说在调度开始前就要确定下来不能随场景变化。在建模上如果决策变量对所有场景共用同一组值就不需要加s标签直接定义成T x 1的变量即可。但如果是目前这类简单的经济调度模型机组始终在运行状态不做启停规划就可以直接简化掉启停变量。4.3 约束批量添加与表达式优化技巧写约束时避免用for循环一个一个加YALMIP支持矩阵化约束表达式。比如电功率平衡约束可以写成Constraints []; % 电功率平衡CHP出力 光伏 风电 购电 电负荷 充电 售电 Constraints [Constraints, P_chp P_pv P_wt P_buy - P_ch - P_sell P_load]; % 热功率平衡CHP热出力 蓄热罐放热 热负荷 蓄热罐吸热 Constraints [Constraints, H_chp H_dis - H_ch H_load];注意如果P_chp、P_ch、P_load等变量都是T x Ns维度那么等式约束会自动批量生成T*Ns个约束。这样写的好处是代码简洁而且求解器性能更好。储能SOC递推约束只能逐时段写循环因为SOC(t1)依赖SOC(t)这没法向量化。不过循还次数只有T次不影响性能。SOC(:, 1) SOC_init; % 初始SOC for t 1:T Constraints [Constraints, SOC(t1, :) SOC(t, :) ... (eta_ch * P_ch(t, :) / Cap_E - P_dis(t, :) / (eta_dis * Cap_E)) * dt]; end这里SOC除以电容量的原因是为了把SOC归一化到0到1之间。4.4 求解配置与结果导出模型建好后调用求解器的设置也要注意。Cplex的求解时间受MIP gap控制默认gap是1e-4对这类模型来说有点过严。实际上设到1e-3就够了求解速度提升明显且方案在实际工程中没有肉眼可见的差别。options sdpsettings(solver, cplex, verbose, 2, cplex.mip.tolerances.mipgap, 1e-3); sol optimize(Constraints, Objective, options); if sol.problem 0 % 求解成功 P_chp_opt value(P_chp); SOC_opt value(SOC); else disp(求解失败); disp(sol.info); end结果导出时建议把每个场景下的电功率平衡情况、热功率平衡情况分别算出来画成堆叠面积图检查是否有负值或者不合理的突变。我见过很多模型求解成功但结果明显不符合物理直觉的情况基本都是约束漏写了或者符号写错了。5. 常见问题与排查技巧实录5.1 求解器报错与YALMIP调试经验求解失败是新手最常见的痛点。我个人经验是分三步排查先看sol.info给出的错误信息再检查是否有变量维度不匹配最后用“松弛约束法”定位问题约束。YALMIP报“Infeasible problem”最常见的原因有三个约束之间自相矛盾、变量的上下界互相冲突、某个参数漏写导致非法约束。排查时把约束注释掉一部分逐步松开很快能找到问题。比如先把储能SOC约束注释掉再跑如果可行了就说明问题出在SOC约束的参数上。5.2 场景数量与求解时间的平衡场景数的选择直接决定求解规模。场景削减到5个时模型很快就能解出来但结果可能受少数场景的代表性影响较大场景数调到30个求解时间会增加数倍但结果的稳定性明显提高。我一般建议削减到10到20个场景之间。特别提醒场景削减之后保留场景的概率分布才是最优的。如果直接用削减前生成的场景去做期望成本计算得到的结果会偏大因为冗余场景的高概率峰谷都被均匀稀释掉了。5.3 常见问题速查表问题现象可能原因排查思路求解器报Infeasible约束矛盾或参数越界逐步注释约束定位冲突约束结果出现剧烈波动缺少爬坡约束或目标函数缺项检查是否有爬坡限制储能SOC长期不变化充放电价格无差异检查分时电价是否真的形成套利空间场景削减后概率之和不等于1削减后未归一化对概率向量重新归一化Cplex求解时间过长二进制变量太多增大MIP GAP容差或减少场景数热功率不平衡蓄热罐约束遗漏检查每个时段的热功率等式约束5.4 数值稳定性与参数灵敏度心得最后分享几个我踩过坑之后的收获帮助你们少走弯路。第一单位统一很重要。功率用kW能量用kWh时间用h价格用元/kWh所有的参数都严格按这个单位体系来模型就不会出现数量级混乱的问题。第二储能初始SOC对第一个时段的调度影响非常大。如果初始SOC设置得不合理比如让小容量储能从满电开始调度优化器可能会在第一时段疯狂放电造成结果失真。常规做法是把初始SOC设成50%左右或者直接把它当作决策变量在优化中一起确定。第三爬坡约束一定要加。很多模型不加CHP机组爬坡约束结果求出来的出力曲线锯齿状剧烈波动现场根本无法执行。热电联产机组的爬坡率通常在每分钟0.02到0.05倍额定功率之间按60分钟折算成每时段爬坡限制去设置。我在实际跑模型的过程中最大的体会是这类随机优化项目建模占了60%的精力代码调试占了30%真正求解只用了10%的时间。很多人一上来就直接写代码结果后面反复折腾。正确的方法是先在纸上把目标函数、约束、变量维度全部列清楚再做参数初始化、再写代码一晚上就能跑通整个流程。按照这个思路去做你的热电联供微网优化模型不仅能顺利求解还能经得起评审老师和数据验证的推敲。