源荷不确定性在热电联供微网里是个绕不开的坎儿。以前做优化调度拿个确定性负荷曲线、固定风电出力就开算看着结果挺漂亮一放到实际运行里就露馅——风电一波动、热负荷一突变原计划的“最优”方案直接失效甚至出现弃风、切负荷或者爬坡跟不上这种尴尬局面。这个课题我前前后后折腾了快两个月从场景生成到随机优化求解再到Matlab代码逐行调通中间踩了不少坑也总结出了一套能直接套用的建模和程序实现流程。这篇就把整个项目的思路、模型推导、程序框架以及我在实际跑数据时遇到的典型问题和排查方法完整拆开讲清楚。1. 项目整体思路与建模方案设计1.1 热电联供微网的物理架构与能量流解析先把这个系统里有哪些设备、能量怎么流动搞清楚。一个典型的热电联供微网核心设备包括CHP机组也就是燃气轮机加余热锅炉的组合、燃气锅炉、电储能、热储能有些配置还会加风电和光伏。电网侧通过一个公共连接点跟大电网交互可以买电也可以卖电。热负荷由CHP的余热回收和燃气锅炉共同供应电负荷由风电、CHP发电、电储能以及网购电共同支撑。这个系统的关键特征就是“电-热耦合”。CHP机组发电的同时产生余热电出力越高热出力也跟着涨但它不能像燃气锅炉那样单独调热这就导致在调度时不能把电和热分开看必须放在一个模型里联合优化。我在建模时用的是可变热电比的CHP模型热出力在一定的电出力区间内可以调整这样比定热电比模型更接近实际运行特性但也让约束条件更复杂特别是可行域的刻画需要仔细处理。从系统运行的角度来看微网调度最核心的任务就是在满足电、热负荷平衡的前提下决定每一台设备在各时段的出力。注意这里说的是“各时段”因为优化模型是离散化的通常以1小时为间隔整个调度周期取24小时。每一时段内风电出力、电负荷、热负荷都有各自的不确定性这就是后面随机优化要入场的地方。1.2 为什么一定要把源荷不确定性写进模型有人可能会问直接按预测值做确定性优化留点备用容量不就行了我在做对比实验时发现问题没这么简单。风电出力的预测误差在接入比例较高时非常显著热负荷更是受天气、用户行为影响波动幅度很大。如果模型里不考虑这些随机性调度结果往往是“贴着约束边界走”——比如某个时段CHP出力刚好卡在最大爬坡速率附近遇上实际风电出力骤降机组来不及补功率就只能切负荷或者高价购电。把不确定性放进优化模型本质上是在回答一个问题面对多种可能发生的源荷场景怎么做一个“在各种情况下都不至于太差”的调度决策。这就涉及到随机规划。我在项目里采用的是基于场景的两阶段随机优化框架。第一阶段在不确定参数实现之前做决策对应的是机组启停、与主网交互的购售电计划等“在这里就要定下来”的变量第二阶段在不确定参数实现后做调整对应的是各台机组的实际出力修正和储能充放电的再调度。简单说就是先定一个大概率能走的方案再根据实际出现的场景去灵活调整。这种两阶段结构的价值在于它允许决策者考虑到未来信息的逐步释放而不是把整个调度周期所有操作都锁死。实际编程的时候这个结构对应的是带场景索引的变量集合代码实现上要特别注意变量维度的设计后面我会详细讲。1.3 技术路线选型为什么用场景法而不是鲁棒优化处理不确定性学术界和工程上大致有三条路随机规划、鲁棒优化、机会约束规划。我做这个项目时考虑过鲁棒优化它给出的解在最恶劣情况下也能满足约束但缺点是过于保守——为了防备那些概率极低、几乎不发生的极端场景正常运行成本会高出不少。我测试过一个配置鲁棒模型的调度成本比随机优化模型高出近7%在现实微网中这个代价不太容易接受。随机规划通过场景集来描述不确定性能够精确反映不同场景的发生概率在经济性和鲁棒性之间取得更好的平衡。而机会约束规划允许某些约束的小概率违反道理上跟随机规划相通但在求解上经常要额外做等价转化处理不当反而增加复杂度。考虑到项目本身要求“Matlab程序实现”和工程可复用的目标最终选择基于蒙特卡洛采样加场景削减的随机优化路线。这个方案对不确定性的刻画比较精细求解复杂度也在可控范围内而且Matlab下YALMIP加商用求解器的组合已经非常成熟调试起来效率高。2. 不确定性建模与随机特征量化2.1 风电出力的随机特征与概率分布模型把不确定性输入到优化模型里不可能直接丢一堆分散的数据点进去必须提炼出数学特征。风电出力的随机特性通常从预测误差入手。实测数据表明风电功率预测误差近似服从均值为0的正态分布但重尾特征明显直接用正态分布会低估极端情况。工程上常用Beta分布或者带偏度的分布来拟合我在这里用的是截断正态分布既能在均值附近保留较高的概率密度又能控制抽样范围避免产生负功率或者超过装机容量的无效场景。具体参数上误差标准差跟预测时段和风电装机容量有关我取的是预测功率的10%到20%这个范围在国内外文献里都有支撑。拟合完分布后用蒙特卡洛采样生成原始场景再配合时序相关性处理。如果完全独立地逐时段采样生成的风电出力场景会出现相邻时段大幅跳变这不符合风电出力的实际平滑特性。我用了一阶自回归模型来引入时间相关性核心思路是让当前时段的预测误差跟上一时段的误差相关相关系数取0.6到0.8之间这样生成的场景序列更贴近真实风电爬坡特性。2.2 电、热负荷的波动特征刻画电负荷和热负荷的处理方式类似但精细度要求更高。电负荷有典型的早晚高峰特征热负荷则跟室外温度强相关波动更平缓但更难精准预测。我的做法是用历史负荷数据统计出典型日各时段的均值再叠加一个随机偏差来描述预测不确定性。这个偏差也假设服从正态分布标准差取该时段平均负荷的5%左右。有个细节容易被忽略热负荷不确定性对系统运行的影响是“慢变量”不像电负荷那样需要秒级、分钟级的响应因此在1小时时间尺度的调度模型里热负荷波动对时段内功率平衡的影响相对有限但对储热罐的容量配置和CHP机组的调度方案影响很大。建模时这个特征不需要显式表达但在结果分析时要意识到热负荷场景的极值往往决定了储热系统的边界状态所以场景削减时不能把热负荷的特殊场景全消掉否则调度结果会偏乐观。2.3 场景生成、削减与概率量化原始场景数量直接决定了优化模型的规模。蒙特卡洛采样一万个场景丢进模型求解器的内存直接爆掉YALMIP建模速度也慢得无法忍受。所以场景削减是必不可少的环节。我采用的削减方法是同步回代消除法。简单说它的核心逻辑就是从原始场景集中逐步剔除那些“跟别的场景最像、被剔除后对整体概率分布影响最小”的场景然后把被剔除场景的概率累积到离它最近的保留场景上。这里的“最近”是用Kantorovich距离来度量的本质上是一种概率测度下的距离。削减前后要检查场景集的统计特征比如削减后的均值、标准差应该跟原始样本基本一致。我一般做两组对照一批削减前1万个场景削减后保留20个均衡场景另一批削减后保留50个用来观察场景数量对优化结果的影响。实际测试下来20个场景已经在成本和计算时间之间取得不错平衡50个场景的边际改进不到1%但求解时间增加近一倍。如果场景数量过多一台普通电脑跑YALMIP加Cplex很可能要等几个小时。2.4 从随机约束到确定性等价转换的技巧随机优化的约束里最棘手的是带随机参数的等式约束比如各场景下的功率平衡。处理办法很直接把这类约束写成“对每一个场景都必须满足”这样就把随机约束拆成了确定性约束的集合。目标函数里的期望项则写成各场景成本加权求和权重就是场景概率。这个过程叫确定性等价转化是程序实现中最核心的一步。但有一个坑必须提并不是所有约束都适合写成“逐场景满足”。比如旋转备用约束如果严格要求所有场景下都要满足备用需求结果会非常保守。更合理的做法是把它写成机会约束允许在低概率场景下备用容量不足。在Matlab里机会约束不能直接丢给Cplex求解需要结合场景法把它转化为带整数变量的约束引入一个0-1变量表示“该场景是否允许备用不足”再用一个约束把不足概率限制在设定值内。这种处理我在后面的代码框架里会给出示例。3. 优化模型构建目标函数与约束体系3.1 目标函数设计经济性与低碳性的平衡这个模型的目标函数包含两部分系统总运行成本和碳排放成本。前者包括燃料成本、购电成本、弃风惩罚成本、切负荷惩罚成本后者则是对CHP机组和燃气锅炉的碳排放施加费用。燃料成本的计算是一个非线性项因为CHP机组的燃料消耗跟电出力和热出力都相关。我在程序里用的是二次函数拟合但Cplex不直接支持二次目标跟混合整数变量混在一起求解所以做了一步处理在机组的运行区间里分段线性化。这个处理非常重要如果不做求解器很难收敛到全局最优甚至直接无法求解。分段数取3到4段就足够逼近原曲线了。弃风和切负荷的惩罚成本系数要合理设置。切负荷惩罚要比购电成本高一个数量级弃风惩罚略低于发电成本这样模型才不会出现“为了保证经济性而故意弃风”的怪现象。我是以系统内最高边际发电成本为基准弃风系数取0.3倍基准切负荷系数取5倍基准实测下来调度行为比较合理。3.2 设备运行约束详解从CHP到储能CHP机组的约束是这里最核心的部分。它的电-热可行域需要仔细刻画。我用的是一个凸多边形近似下包络线由最小电出力和对应的热出力决定上包络线受最大燃料输入量和最大热回收量限制。热出力还有一个独立的上限。这些约束写成线性不等式组每台CHP机组大概用十几行约束就能描述清楚可行域。电储能约束是标准的储能电量动态变化方程、充放电功率上下限、充放电效率、以及调度周期始末电量相等。这里有个编程细节——储能电量的连续变量要用上一个时段的电量、充电功率和放电功率联合表达充电效率和放电效率不相等所以不能简单用“净功率”替代否则电量的递推公式会出错。我一开始图省事用净功率建模结果电量越算越不对查了半天才发现是效率项漏了。热储能跟电储能结构类似但需要注意热储能的损耗比电储能小得多时间常数很长可以扛热负荷在几个小时内的波动。热储能的容量约束在场景削减比较彻底的时候尤其重要因为场景少意味着极端热负荷场景可能被削减掉热储能的运行边界就显得比较紧。3.3 旋转备用约束与系统可靠性保障旋转备用约束是整个模型里最能体现“随机特征”的地方。它要求系统在任何场景下都有足够的可调容量来应对预测误差。在场景框架下备用约束的量化方式要跟场景概率关联。我是用机会约束来处理的在95%的概率下系统总可用上调容量不小于该场景下负荷与风电差值预测的偏差需求。在Matlab代码里这个约束用一个大M法加整数变量来实现。常规写法是for s 1:Ns % 上调备用约束 sum_ramp_up(:,s) binary_up(:,s) * M reserve_required(:,s); end sum(binary_up) Ns * (1 - alpha); % alpha是置信水平这里的M取值不能太大太大会引起数值问题我在程序里取备用需求的5倍就够。这个约束写不好模型很容易出现“所有场景都不需要备用”的不现实解或者反过来全部场景都加备用导致成本虚高调M值的过程也算是这个项目里一个比较磨人的环节。4. Matlab程序实现与求解全过程4.1 求解器选型与YALMIP环境配置这个项目用的是YALMIP作为建模语言底层求解器选了Cplex。YALMIP的好处是语法接近数学表达调试方便而且能自动识别模型类型把混合整数二次规划交给合适的求解器。Cplex在求解这类中大规模的混合整数线性规划问题上非常稳处理几万行约束没问题。环境配置上注意Matlab版本跟YALMIP的兼容性。我用的Matlab R2023b配的是YALMIP R20230616版本Cplex 12.10。安装时有个常见的坑把求解器路径添加到Matlab之后还要在YALMIP里显式声明求解器否则它会默认选一个不适合的求解器来跑混合整数问题。options sdpsettings(solver,cplex,verbose,2); options.cplex.mip.tolerances.mipgap 0.0001; % 设置MIP间隙控制求解精度求解速度跟MIP间隙设置关系很大。如果追求快速得到一个工程可用的解把间隙放到0.01就够了如果要用于论文数据建议设到0.0001代价是求解时间可能翻倍。我做敏感性分析时就是用0.0001的精度来算跑一次两万多行约束的模型大约需要20分钟。4.2 关键代码模块拆解参数定义、场景生成与建模求解整段程序的架构分四个模块参数初始化模块、场景生成模块、优化建模模块、结果输出模块。场景生成模块单独抽出来做成函数输入是历史数据和分布参数输出是削减后的场景集和对应概率。优化建模模块用YALMIP语句描述目标函数和约束主体代码大约150行结构非常清晰。变量定义是程序里最需要小心的部分。两阶段随机优化中第一阶段变量没有场景索引第二阶段变量带场景索引。在YALMIP里带场景索引的变量要定义成矩阵例如P_chp sdpvar(N_gen, Ns, full)其中第一维是机组编号第二维是场景编号。这个维度设计一旦搞错后面所有约束写起来都会乱套尤其在处理“第一阶段变量对所有场景都是一样的”这个约束时要用repmat把决策变量扩展到每个场景再跟场景变量联立约束。核心建模代码节选% 第一阶段常规机组出力不随场景变化 P_chp sdpvar(1, 24, full); % 24时段CHP电出力 H_chp sdpvar(1, 24, full); % 24时段CHP热出力 % 第二阶段各场景下的调整量 delta_P sdpvar(1, 24, Ns, full); % 场景调整出力 P_buy sdpvar(1, 24, Ns, full); % 各场景购电功率 % 约束部分对每个时段、每个场景写平衡约束 for t 1:24 for s 1:Ns Constraints [Constraints, ... P_wind(t,s) P_chp(t) delta_P(t,s) P_buy(t,s) P_dis(t,s) ... P_load(t,s) P_chg(t,s)]; end end这些约束看起来不复杂但写成循环嵌套之后规模会迅速膨胀。24个时段乘20个场景就是480组平衡约束再算上机组爬坡、储能、网络、备用等约束总约束行数会到两万以上。YALMIP建模本身没问题但要注意别用中文变量名语法解析会有问题。我在调试时把所有约束保存到Constraints对象里然后在求解前用check(Constraints)检查可行性能快速定位是哪组约束出了问题。4.3 结果可视化与调度方案分析求解完成后结果的可视化对分析很有帮助。我通常画四类图第一类是各场景下风电与负荷的时序曲线叠加观察场景集是否覆盖了典型的波动区间第二类是CHP电、热出力的P-H图电热运行点图能看到调度解是否落在可行域边界上第三类是储能充放电功率和电量的时序曲线检查是否存在不合理的频繁充放第四类是各场景下的成本分布直方图看期望成本附近的情况。有一次我画完P-H图发现CHP机组很多时段的运行点都贴在可行域的下包络线上这说明热负荷不足正在限制电出力跟现实情况对上了——燃气轮机在低热负荷时电出力上不去系统只能靠购电补缺口。这类图上的直观发现反过来又能帮你校验模型约束是不是写得太紧要不要放宽热出力下限。5. 典型算例结果与方案分析5.1 确定性模型与随机优化模型的调度差异为了说明考虑不确定性带来的价值我用同一组基础数据跑了两个模型一个是不考虑随机性的确定性模型把风电和负荷都取预测期望值另一个是完整的随机优化模型。对比结果非常直观。确定性模型给出的24小时总运行成本在数值上更低但这个“低”是虚假的因为它假设预测完全准确不需要为不确定性预留额外的调整能力。如果拿确定性模型求出的调度方案去回代测试所有场景会出现大量时段功率不平衡也就是切负荷实际成本反而比随机优化模型高出约12%。随机优化模型给出的方案在各场景下都能保证功率平衡某些场景下虽然购电量多一点但这是在“用一点经济性换可靠性”从期望成本角度反而是更优的。这个结果也回应了很多人的疑问随机优化不是为了省钱而是为了在现实的不确定环境中获得更真实的调度策略。5.2 不确定性水平变化对运行成本的影响我对风电预测误差标准差做了从5%到30%的敏感性分析。结果是一条单调上升的成本曲线但斜率在不同区间差异明显。误差从5%升到15%这个区间成本上升相对平缓一旦超过20%成本急剧上升。原因在于误差越大备用约束就更容易达到临界点迫使系统更多依赖快速但昂贵的调节资源。这里有一个对工程实践有意义的观察与其不断追求更精准的预测从20%误差压到10%不如在调度模型里合理量化备用需求和场景覆盖同样能明显改善运行经济性。当然这俩不是二选一但很多项目在预测模型上花了大量精力却忽略了优化模型本身对不确定性的消化能力是很可惜的。5.3 场景数量与求解质量的平衡分析前面提到了场景削减后保留数量的选择。这里给出具体的测试数据。原始蒙特卡洛生成1万个场景分别削减到10、20、50、100个各跑一次完整的随机优化。保留场景数求解时间分钟期望成本万元相对偏差%102.519.621.85208.319.310.245021.719.280.0810056.419.260.03可以看到20个场景已经能把期望成本偏差控制在0.3%以内再往上加场景精度提升非常有限但计算时间几乎是线性增长。如果只是做原理验证或者教学演示10个场景也够用但做科研或者工程决策建议至少保留20个场景。场景削减算法本身的代码不算复杂但我在实现时有个心得削减后务必重新归一化场景概率让所有场景概率之和严格等于1否则模型解的每一个约束都会带一个微小的偏移累积起来对成本计算的影响虽然不大但对结果一致性检验很要命。6. 常见问题与排查技巧实录6.1 求解器无解或不可行先查哪个约束YALMIP返回infeasible的时候最有效的排查方法不是瞪着眼睛看代码而是用约束松弛法。把目标函数临时改成0并把约束拆成几组逐个测试可行性。如果去掉储能约束就能解那就是储能电量递推或者容量约束出了问题。最常见的原因有两个一个是储能初始电量跟末尾电量约束冲突24小时的充放电能力根本补不上电量差额另一个是CHP可行域约束写漏了一两条导致某些时段的电热组合实际不可达。我处理这类问题时的做法是把约束一条条打印出来检查每一行的系数矩阵。YALMIP里可以用Constraints(i)访问单条约束配合sdisplay查看其表达式能非常直观地看出有没有写反不等号方向、有没有把变量维度搞错。这个方法看起来笨但排查效率比瞎猜高得多。6.2 场景削减导致概率分布偏移有一种情况需要警惕场景削减后某些稀有但影响严重的事件比如极低风电加高峰负荷同时出现可能被完全消除导致优化结果过于乐观。排查方法是对比削减前后场景集中风电出力和负荷的极值、标准差。如果削减后的标准差显著小于原始样本多半是削减过度或者Kantorovich距离计算有误。我在程序里加了一个检验步骤削减之后自动计算场景集的均值和标准差跟原始样本对比偏差超过2%就报警提醒检查削减算法的参数设置。这个自动校验在批量跑实验时特别有用不然人很容易被一堆中间结果淹没。6.3 求解时间过长从模型规模与求解器参数入手模型规模上来之后求解时间的优化是绕不开的。我有三个实用技巧。第一尽量把整数变量的数量降下来。CHP机组的分段线性化会引入分段变量如果能证明某些机组在调度周期内大概率不会启停可以预先固定启停状态大幅减少整数变量数。第二设置合理的MIP gap。做批次实验时用0.01的gap先粗筛对结果有明显影响的方案再用0.0001精算能节省大量时间。第三Cplex里开启并行计算和预求解功能options.cplex.mip.strategy.probe 2这种参数组合对特定模型能提速30%以上。有一次我跑200个场景的完整随机优化内存占用一度超过20GB直接把电脑跑死了。后来发现是YALMIP在构建大规模稀疏矩阵时变量排序方式影响了效率。把变量定义从sdpvar(Ns, 24)改成sdpvar(24, Ns)求解时间缩短了将近一半。这个维度顺序的影响在YALMIP文档里没有明确说明但实际效果非常显著建议大家在定义大矩阵变量时多试几种维度排列。6.4 数据单位与数值尺度引起的隐性错误这个坑非常隐蔽。风电出力是兆瓦级而惩罚成本系数可能是元每兆瓦时数值动辄上百万两者在约束矩阵里的尺度差了好几个数量级。Cplex对数值尺度非常敏感太悬殊的系数范围会导致数值困难求解器给出的结果即使“收敛”也不可靠。解决方法是归一化。把功率统一折算成标幺值或者统一数量级比如成本系数全部换算成万元每兆瓦时让约束矩阵里的系数集中在0.01到100这个区间。我吃过这个亏数据没归一化之前同一个模型跑出来的成本结果在重复实验时会有微小波动归一化之后就稳定了而且求解速度也明显提升。7. 项目扩展方向与实际应用建议7.1 从随机优化到分布鲁棒优化的升级路径如果后续想在理论上再进一步可以尝试把随机优化升级为分布鲁棒优化。核心思路是不再假设不确定参数的分布是精确已知的而是用一个包含真实分布的不确定集来描述。这样做的好处是当天数据不足、分布估计不准时模型依然能保证一定的可靠性。代价是模型规模和求解复杂度都会上升需要在YALMIP里引入半定规划的相关约束求解器也要换成支持半定规划的版本。我做过一个小规模测试分布鲁棒模型的成本比随机优化模型高约3%但对分布估计偏差的容忍度大幅提升。如果你的项目里历史数据质量一般这个方向非常值得探索。7.2 把模型跟实时预测系统联动的工程建议模型本身是一个离线优化工具但如果要投入实际运行最好跟实时负荷和风电预测系统对接。具体做法是每15分钟更新一次预测数据重新生成场景集滚动求解未来4到6小时的调度计划。滚动优化能及时吸收最新信息弥补单次优化对远期预测误差无能为力的短板。我在程序框架里预留了滚动优化的接口核心就是把整个调度周期从24小时改成移动窗口每步推进1小时。Matlab里用循环控制就能实现但要注意每次滚动时初值要继承上一轮的解否则机组启停状态会产生剧烈震荡。这个细节在工程落地时非常重要。7.3 代码封装与复用的心得整个项目跑通之后我花了一些时间把程序整理成模块化结构。场景生成、模型构建、求解、结果分析分别封装成独立函数输入输出接口定义清楚。这样在更换历史数据、调整设备参数甚至修改目标函数时只需要改配置文件和对应模块函数不用重写整段逻辑。对于要长期维护的项目建议把所有可调参数集中放到一个结构体里比如params.wind.sigma、params.storage.capacity这种层级结构。Matlab的struct虽然简单但在代码可读性和后期调参效率上的提升非常明显。初学者往往把参数散落在程序各个角落后面找一个系数要翻半天代码这个习惯还是尽早改掉比较好。我在实际跑这个项目的过程中最大的感受是理论模型和程序实现之间的距离往往不在数学推导里而在那些不值一提的小细节——一个约束漏写了索引一个效率参数方向搞反一个场景没有归一化概率都可能让结果差之千里。但反过来只要把不确定性建模的路子走通把场景生成、约束表达、求解调试这些环节形成一套自己的标准流程这类热电联供微网优化问题基本就是换数据、换参数、重复利用框架的问题了。希望这篇文章里的思路和坑能帮你少走一些弯路。