计及源荷不确定性的综合能源系统容量配置与运行调度优化

做综合能源系统优化的人应该都遇到过这种憋屈时刻费了一个月把容量配置模型搭好了用确定性预测曲线跑出来的方案投资成本降得挺漂亮结果把真实的风电、光伏和负荷曲线塞回去一验证系统动不动就弃风、切负荷甚至储能直接干到0%不满足约束。问题就出在题目里那六个字——计及源荷不确定性。综合能源生产单元的运行调度和容量配置本质上是两个时间尺度的耦合决策如果不提前把不确定性建模进去你得到的就是一套看起来很优、实际很脆的方案。这篇博文我想完整梳理一下在Matlab里实现计及源荷不确定性的综合能源生产单元运行调度与容量配置优化的整个思路包括不确定性怎么刻画、双层模型怎么搭、Yalmip配合求解器怎么落地以及我自己调模型时踩过的那些坑。同样做这个方向的朋友特别是刚入门的硕士研究生可以把这篇文章当成一份从零到一的项目模板来参考。1. 问题拆解为什么要把运行调度和容量配置放在一个框架里1.1 两个时间尺度在模型里的矛盾容量配置解决的是建多大的问题运行调度解决的是怎么开的问题。前者是年尺度甚至全生命周期尺度的投资决策后者是小时级、分钟级的操作决策。很多人习惯两步走先按典型日确定容量再在固定容量下分析运行策略。听起来没毛病但仔细想想就发现里面藏着矛盾——容量配置需要知道设备在未来各种工况下的利用率而这个利用率恰恰取决于运行调度策略反过来运行调度的可行域又完全被设备容量框住了。两个问题互为前置条件标准解法就是把它们放进同一个优化框架里也就是双层优化。不过这并不意味着一定要用严格的双层优化数学形式。实际工程里常用的大概有三类处理方式第一类把容量变量和运行变量全部拉平做成一个大混合整数线性规划(MILP)约束包含所有场景的运行约束这叫全场景扩展。第二类采用两级迭代外层用启发式或群智能算法搜索容量方案内层对每个方案求解运行调度子问题返回运行成本给外层。第三类利用KKT条件把下层问题转化为上层的约束形成单层数学规划(MPEC)再线性化处理。在Matlab实现时我个人的经验是如果算例规模不大比如单站多能互补系统设备类型不超过十种场景数压缩到十几个第一类全场景扩展最省事直接用Yalmip建模交给Gurobi/CPLEX求解稳定性最高。如果要做园区级甚至区域级网络节点数上来了或者想在优化里同时考虑网络拓扑那么第二类内外迭代框架更灵活。1.2 不确定性要是被忽略配置结果会偏成什么样在动手建模之前先想清楚一个问题源荷不确定性到底会让结果往哪个方向偏。风电光伏的出力预测误差是典型的低谷乐观——预测曲线看着午后光伏有几十千瓦出力实际可能一阵云飘过来直接归零。负荷侧也有类似问题尤其是冷热电联供系统里的热负荷用户行为随机性强尖峰往往比预测值晚到或者持续时间更长。如果用确定性模型做容量配置它会觉得光伏出力很稳定、负荷曲线很光滑于是倾向于减少储能和可调机组的容量因为反正预测的光伏够用。等到实际运行光伏波动要靠燃气轮机扛负荷尖峰要靠储能顶设备容量不足系统只能切负荷。切负荷在模型里往往被惩罚成天价成本但你做容量配置的人不会看到这个因为你根本没把不确定性场景放进约束。所以确定性模型的容量配置结果是在牺牲可靠性换取低成本而且这个可靠性损失在模型里完全不可见直到实测验证才会爆炸。反过来如果过度保守把所有极端场景都按最恶劣情况约束储能和机组容量会配得特别大投资成本高到项目没法落地。所以计及源荷不确定性的本质是在投资成本、运行成本和可靠性之间找平衡点而平衡的标尺就是你如何定义不确定性集合或场景集。2. 源荷不确定性的数学刻画与场景削减2.1 用场景集替代单一预测曲线蒙特卡洛生成与聚类削减在做随机规划时最常见的方式是场景法。思路很直白既然预测曲线不准那就把每种可能的偏差都变成一个场景每个场景就是一条完整的风电/光伏/负荷时序。场景怎么来标准做法是假设预测误差服从某种概率分布然后用蒙特卡洛抽样。以光伏为例实际出力可以表示为预测值加上误差项P_pv_actual(t) P_pv_forecast(t) delta_pv(t)其中误差项delta_pv(t)往往假设为正态分布但要注意两个问题一是光伏出力在接近零的时候不可能抽成负数所以要做截断处理二是同一时刻不同时间点的误差之间存在相关性不能用独立抽样否则生成出来的场景特别毛糙缺乏连续天气过程的那种持续性。我在Matlab里一般用Cholesky分解或者基于历史误差协方差矩阵生成多变量正态分布随机数这样能保留时间相关性% 根据历史预测误差构造协方差矩阵 Sigma % 对 Sigma 做 Cholesky 分解 L chol(Sigma, lower) % 生成白噪声向量 z再通过 mu L*z 得到相关误差序列 errors mu L * randn(nTime, nScenarios);负荷侧类似但负荷还有明显的日内周期性。直接在原始序列上叠加正态误差容易破坏周期性更好的做法是先把负荷分解为典型日标幺曲线 × 日峰值然后在日峰值和标幺曲线两个层面分别加扰动。这样生成的场景既保留负荷形状又能体现峰值波动。场景生成只是第一步。蒙特卡洛抽个三五千个场景直接扔进优化模型基本跑不动。几万个时间点乘上几千场景约束规模瞬间爆炸。所以必须做场景削减。主流的削减方法有快速前向选择、后向消除和k-means聚类。我自己在Matlab里用得最多的是k-means因为它不仅能把场景数量压到十个以内还能顺便做场景合并给每个聚类中心一个权重代表这一簇场景的概率。% 将每个场景展开成行向量比如72小时的新能源出力曲线 scenarioMatrix reshape(allScenarios, nScenarios, []); % 用 kmeans 聚类为 nCluster 个典型场景 [idx, centers] kmeans(scenarioMatrix, nCluster, MaxIter, 500); % 每个聚类中心的权重 该簇样本数 / 总样本数 weights histcounts(idx, size(centers,1)) / nScenarios;聚类削减后的场景虽然数量少了但基本保留了原始场景集的统计特征。这里有个细节削减的对象是整个时序场景而不是每个时间点独立聚类。如果你只对每个时刻单独聚类会破坏时序相关性跑出来的结果可能出现这小时光伏出力很高下一小时突然很低这种不真实组合。2.2 鲁棒优化下的不确定集合设计如果说场景法是从概率视角描述不确定性那鲁棒优化就是从最坏情况视角描述不确定性它不需要知道精确分布只需给不确定性一个范围集合。集合怎么设计直接决定模型的保守程度。最简单的盒式集合U { P_wind(t) in [P_wind_forecast(t) - delta_wmax(t), P_wind_forecast(t) delta_wmax(t)] }盒式集合约束所有时间点同时达到最坏情况结果会异常保守。更合理的是带预算的集合通过引入预算参数Gamma限制整个调度周期内不确定性实际偏离预测值的总次数或总幅度sum_t |delta_w(t)| / delta_wmax(t) Gamma_wGamma_w越大越保守等于时间周期长度时退化为盒式集合。在Matlab里用Yalmip配合鲁棒优化建模可以直接用鲁棒约束的方法定义不确定变量但很多老版本或者不装额外工具包的环境下更通用的做法是把它转写成对偶形式把不确定变量消掉得到等价的确定性鲁棒对等模型。对偶转化的核心思想对于线性约束ax bu c其中u属于一个锥多面体不确定集合可以把这个约束转化为有限个确定性约束。过程不复杂但对初学者来说容易算错对偶变量方向。我自己实操时如果不是论文要求必须给出纯鲁棒对等模型通常更倾向于用场景法加一个最恶劣场景作为鲁棒性的近似校验。这个最恶劣场景可以从历史观测中挑也可以人为构建一个极端组合——比如最大负荷日叠加光伏零出力、风电最大爬坡等。跑一遍调度模型看看各设备是否还有裕度。这种做法工程上完全够用且比纯鲁棒模型容易调试得多。3. 双层优化模型在Matlab里的落地思路3.1 上层容量配置投资决策与设备选型综合能源生产单元一般包含光伏、风电如果有资源条件、燃气轮机、电储能、热储能、电锅炉、吸收式制冷机等。容量配置的目标是决定每种设备的额定容量。设备投资成本通常可以近似为额定容量的线性函数再加上固定安装费所以上层目标函数大致为min invCost sum_i (c_i * S_i f_i * binary_i)其中S_i是连续变量表示额定容量binary_i是是否配置该设备c_i是单位容量投资成本f_i是固定投资成本。加binary变量是因为很多设备有最小配置门槛比如一台燃气轮机要么不装要装就至少装0.5MW这种非凸约束必须靠整数变量表达。上层约束还包括最大可用土地面积、最大总投资预算、设备容量上下限、能源品类多样化要求比如可再生能源渗透率不低于一定比例等。这些约束一般比较简单关键是它们会传递到下层运行调度模型里作为设备容量的参数。3.2 下层运行调度线性化机组模型与经济调度下层模型是在给定设备容量后求解典型场景下的最优日运行策略目标通常是运行成本最小包含购电费用、燃气费用、弃电惩罚、负荷失配惩罚以及设备启停成本。下层模型的决策变量包括每时刻每台机组的出力、储能充放电功率及荷电状态、电网交互功率、各能流的分配比例等。约束包括功率平衡约束、设备出力上下限、爬坡约束、储能能量平衡约束、网络或输能容量约束、以及设备运行特性线性化约束。燃气轮机的发电效率不是一个常数而是随负载率变化的曲线。如果直接引入二次效率曲线模型就变成MIQP混合整数二次规划求解难度会上一个台阶。实操中我通常做分段线性化把负载率区间分成3到5段每段对应不同的发电效率效率乘以燃料热值得到燃料消耗量用分段线性表达式替代二次曲线。这样精确度足够模型仍然是MILP求解器处理起来非常顺手。储能的约束是典型线性约束但有个易错点充放电功率往往乘一个效率因子导致同一个时刻既不能同时充电和放电。只用功率约束不表达互斥关系时求解器可能给出一边充一边放的荒唐结果。虽然在有损耗的情况下这种荒唐解是次优的但数值上有时会出现。稳妥做法是引入两个0-1变量分别表示充放电状态并加约束u_charge u_discharge 1。3.3 上下层迭代求解的两种工程方案方案A双层嵌套。外层可以用MATLAB自带的遗传算法ga函数或者粒子群算法particleswarm搜索容量方案内层针对每个方案调用YalmipGurobi求解MILP运行子问题。这个方案实现起来最容易理解代码结构清晰缺点是每个外层个体都要多次求解内层模型计算量大。我跑典型案例时如果场景数为10种群规模100迭代50次意味着要解50000个MILP单个MILP秒级求解整体也要十几个小时。所以一定要控制规模和求解精度。方案B全场景单层MILP。把N个场景的运行变量全部并列容量变量共享目标函数是投资成本加各场景运行成本的期望值按场景权重加权。约束把所有场景的约束一并列出共享的容量变量作为场景间耦合约束。这种做法的优点是模型理论上能求到全局最优解不用像启发式那样担心收敛到局部最优。缺点是场景数和系统规模变大后变量数量线性增长可能内存爆炸。方案CBenders分解。把问题分解为主问题和子问题主问题决定容量子问题检验运行可行性并生成最优割。Benders分解本质上就是为这类问题设计的但实现复杂度高而且需要处理整数变量的可行性割如果不是论文需要完全自定义算法不太建议从零手写。网格环境或者教学可以把Benders作为扩展方向但实际项目我更推荐先用方案B或方案A。作为一个折中我自己常用的套路是先用方案B跑一个场景数为8到10的小规模模型得到一组可作为参考的最优解然后把这组解作为初始可行解喂给方案A的粒子群算法让它在大搜索空间里微调。这样既保证可靠初值又能探索更多设备组合可能效率比直接跑大场景单层模型高很多。4. 代码实现中的关键细节与求解技巧4.1 用Yalmip定义不确定变量和场景在Matlab里Yalmip是建模层Gurobi或CPLEX是底层求解器。用Yalmip定义多场景模型时核心技巧是下标里加场景维度。比如光伏出力变量不用P_pv(t)这种只有时间维度的定义而是用P_pv(t, k)表示第k个场景t时刻的光伏出力。Yalmip支持这种矩阵化定义约束也用循环或矢量化写。一个典型的容量-运行联合MILP模型在Yalmip里的骨架长这样% 定义容量变量 S_pv sdpvar(1, 1); % 光伏容量 S_wind sdpvar(1, 1); % 风电容量 S_gas sdpvar(1, 1); % 燃气轮机容量 S_es sdpvar(1, 1); % 电储能容量 % 定义运行变量场景集合K时间集T P_pv sdpvar(nTime, nScenarios); % 实际消纳的光伏功率 P_gas sdpvar(nTime, nScenarios); % 燃气轮机出力 P_ch sdpvar(nTime, nScenarios); % 储能充电 P_dis sdpvar(nTime, nScenarios); % 储能放电 SOC sdpvar(nTime1, nScenarios); % 荷电状态 % 约束定义示例 Constraints []; for k 1:nScenarios for t 1:nTime Constraints [Constraints, P_pv(t,k) S_pv * availableFactor(t,k)]; Constraints [Constraints, P_gas(t,k) S_gas]; Constraints [Constraints, P_ch(t,k) P_dis(t,k) S_es]; % 简化 end end % 目标函数 Objective investCost sum(weights .* operateCost); % 求解 optimize(Constraints, Objective, sdpsettings(solver,gurobi,verbose,2));关键点是所有容量变量不能带场景下标所有运行变量必须带场景下标。Yalmip自动识别这种维度差异约束写起来很顺畅。但要注意别混淆矩阵乘法和逐元素乘法Yalmip对sdpvar矩阵的号规则和Matlab一致不是逐元素乘一定要用.。4.2 从求解器报错看模型的病态来源整数变量、双线性项与Big-M初学阶段最常见的报错是Non-convex quadratic或者Unknown solver这类问题大多源于引入了两个连续变量相乘的双线性项。比如储能充放电互斥约束如果直接用P_ch * P_dis 0这就不是线性约束。正确做法是用二进制变量加Big-M线性化u_ch binvar(nTime, nScenarios); u_dis binvar(nTime, nScenarios); Constraints [Constraints, P_ch P_ch_max * u_ch]; Constraints [Constraints, P_dis P_dis_max * u_dis]; Constraints [Constraints, u_ch u_dis 1];Big-M的取值也有讲究。M值必须大于变量理论上能达到的最大值但又不能太大否则会引入数值病态。比如储能最大充电功率不超过储能容量那M取S_es就行了别随便写个1e6。数值病态的典型表现是求解器报告Numerical trouble、迭代停滞或者给出违反约束的解。如果你遇到这类报错先系统检查一遍所有Big-M约束把M压到理论最小值的1.1到1.5倍。另一类常见报错是Unbounded objective或Infeasible problem。无界问题通常是因为某个变量缺少上下限约束。比如电网交互功率P_grid如果没有定义上下限目标函数里购电成本为正理论上P_grid会趋向负无穷售电除非你在目标里加了售电收入否则这个解就是无界的。所以每个有成本系数的变量一定要检查它是否有对应的约束把它框住。4.3 计算效率优化热启动、约束削减与并行化当模型规模上来了求解时间不是线性增长而是指数增长这时候效率优化就变得重要。我的经验有以下几条第一热启动。Gurobi支持传MIP start也就是把已知的可行解作为初值。怎么获取已知可行解可以用之前说的确定性模型跑一遍把容量和运行变量的值作为start值传进去assign(S_pv, 1.2); % 给变量赋值 assign(P_gas, initialGuessMatrix); optimize(Constraints, Objective, sdpsettings(solver,gurobi,mipstart,1));热启动能大幅减少求解器无谓的搜索特别是你的模型里有大量0-1变量时。第二约束削减。很多约束在大多数场景下是不紧的。比如储能的上下限约束在低光照低负荷的场景下储能根本不会满充满放这些约束对最优解没影响。你可以先把约束全部写上求解一轮看哪些设备的影子价格一直为零再把那些场景或时段对应的约束删掉。这种操作在工程上叫Lazy Constraints虽然Gurobi也有lazy的API但在Yalmip里用起来不方便我更倾向于手动做一轮两轮削减。第三并行化。如果用的是方案A外层启发式那Matlab的parfor非常香。把粒子群或遗传算法的适应度函数改成parfor循环几十个workers一起解MILP时间能从十个小时压到一两个小时。但要注意parfor里调用Yalmipsolver需要每个worker都有相应的工具箱和求解器并且Yalmip的变量传递不太方便建议把内层模型封装成一个独立的m函数参数通过结构体传入求解结果通过结构体传出。5. 算例测试一个微电网-多能互补系统的结果对照5.1 算例参数与场景设置为了验证模型我搭建了一个典型的综合能源生产单元包含光伏50kW可扩展、燃气轮机热电联产可扩展、电储能可扩展、蓄热罐可扩展系统与上级电网有连接同时向园区供应电和热。设备参数如下设备单位投资成本运行效率/性能参数备注光伏8500元/kW标况效率20%温度修正系数-0.4%/°C年利用小时数约1100h燃气轮机4200元/kW发电效率40%热电比1.2最小出力30%启停成本500元/次电储能1800元/kWh充放电效率95%SOC范围15%~95%循环寿命按6000次折算蓄热罐900元/kWh散热损失5%/h最大容量不超过日峰值热负荷3倍负荷和光伏数据来自某园区历史数据时间分辨率取1小时典型日选取春夏秋冬各一个。不确定性建模时对光伏出力生成50个场景负荷生成50个场景两两组合得到2500个原始场景再用kmeans聚类成10个典型场景每个场景有对应的权重。对比模型包括确定性模型直接用预测曲线、随机规划模型10场景、鲁棒优化模型盒式集合Gamma6。5.2 不同不确定性处理方式下的配置与运行结果确定性模型跑下来的最优配置是光伏35kW燃气轮机30kW电储能80kWh蓄热罐120kWh。总年化投资成本约11.6万元运行成本约9.8万元总成本21.4万元。随机规划模型的结果是光伏32kW燃气轮机45kW电储能120kWh蓄热罐180kWh。投资成本上升到13.2万元但因为储能和蓄热罐多了可以在电价低谷多充电、热价低谷多蓄热运行成本降到了8.4万元总成本21.6万元——和确定性模型几乎一样但可靠性大幅提高。鲁棒优化Gamma6的结果则是光伏28kW燃气轮机55kW电储能150kWh蓄热罐230kWh。投资成本15.8万元运行成本7.6万元总成本23.4万元。显然鲁棒模型更保守容量更大可靠性最高的代价是总成本增加了9%。这三种结果的对比最有意思的点在于如果只看总成本确定性模型和随机规划模型差不多但如果你把实际光伏出力偏差场景带回去验证确定性模型的失负荷率是12%而随机规划模型只有0.5%。这就是为什么我说确定性模型的低成本是一种虚假的低成本因为你没把不可靠带来的损失算进去。5.3 结果能说明什么成本与可靠性的帕累托边界容量配置和运行调度联合优化的本质是在成本与可靠性之间找帕累托前沿。把鲁棒优化的Gamma参数从0变到24对应全天每个时刻都不允许偏差可以得到一串配置结果Gamma越小越接近确定性方案Gamma越大容量越往保守方向走。把每种方案的期望失负荷率EENS和总成本画在二维坐标里就是一条单调下降的曲线。项目决策者可以根据对可靠性的容忍度在这条曲线上选点。在Matlab里这一步不难。你只需要在模型里设置Gamma为一个变量循环跑多个场景记录下不同Gamma下的总成本和EENSgammaValues 0:2:24; for idx 1:length(gammaValues) % 更新不确定预算参数 % 求解模型 % 保存总成本和EENS end % 用plot或scatter绘制帕累托前沿不过还有一个值得注意的点储能配置在帕累托曲线上会出现平台期。当储能容量从80kWh增加到120kWh时系统可靠性提升明显但超过150kWh以后可靠性提升非常缓慢。原因在于储能主要解决的是短时尺度的功率波动而长时尺度的能量缺额需要靠燃气轮机和电网交互光堆储能没法覆盖连续多日的极端天气。所以做容量配置时别盲目追求储能容量要看看缺的是什么类型的灵活性。6. 调试中踩过的坑与避坑建议6.1 蒙特卡洛场景数选多少才够生成场景时很多人会问蒙特卡洛抽样到底抽多少个合适我的经验是先抽500个聚类成10个场景跑一版模型再抽2000个聚类成10个场景再跑一版比较两次优化结果。如果容量和总成本差异在3%以内说明场景集已经收敛。如果差异很大就继续增加原始场景数。另一个判断指标是看聚类中心权重的稳定性——如果重复多次聚类得到的聚类中心序列变化不大说明场景数量足够反之则需要增加抽样数或改用更好的聚类方法。我见过有人直接把10000个场景扔进模型跑了两天没出来这是完全没必要的。场景削减的目的不是保留尽可能多的信息而是在保留统计特征的前提下尽量压缩计算规模。10个场景和20个场景的结果往往差不了太多但计算时间可能要翻好几倍。6.2 无可行解的排查思路模型报infeasible是最折磨人的。我的排查习惯是从最简单的模型开始逐步加复杂度。比如先跑一个确定性单场景模型确认能出可行解然后加第二个场景跑一次看是否可行如果第二个场景加进去就无可行解那基本可以断定是新增场景里的某个约束冲突了。具体到综合能源系统最常见的原因是储能的初值SOC设置与日内能量平衡矛盾。比如你给定了SOC(0)50%但模型要求一天结束时SOC(T)50%这是为了保证日循环如果系统被迫把储能放光来满足某一时段的负荷而夜间没有足够的充电空间就可能无解。这时候要么放宽SOC终值范围要么在目标函数里给SOC偏离初值加惩罚项允许储能跨日调度。这种方法不仅仅是为了求解可行性也更好地贴近实际因为储能本来就是可以跨日调度的。6.3 储能荷电状态初值的敏感性储能SOC初值对运行调度结果影响非常大。很多人习惯设为50%但在多场景随机规划中每个场景的SOC初值应该保持一致否则模型会把初始储能当成免费的灵活性资源导致某些场景过度放电。正确做法是让SOC(1,k)对所有场景k都是一个共享变量SOC_initial sdpvar(1,1); for k 1:nScenarios Constraints [Constraints, SOC(1,k) SOC_initial]; Constraints [Constraints, SOC(end,k) SOC_initial]; % 日循环约束 end这样模型会自己优化出每个典型日的初始荷电状态而不是人为指定一个拍脑袋的数值。我实测过固定SOC50%和SOC为优化变量两种方案燃气轮机的容量配置差异能达到20%以上因为固定的初值限制了储能在全天运行时段的灵活性。另一个相关的小陷阱是SOC上限设置。如果储能容量非常大但SOC上下限设成15%~95%那实际可用容量只有额定容量的80%而投资成本是按额定容量算的。想在某个算例里体现10%的社会余量没问题但论文或者项目汇报时要明确说明否则别人会误以为你的储能利用率只有80%。最后说一个我用Matlab跑这类模型多年的习惯不管模型多小从来不用默认求解设置跑MILP而是会手动把MIPGap设置到0.1%或更小Matlab里Gurobi通过sdpsettings的gurobi.MIPGapAbs等字段设置。初始求解时可以先放宽到5%拿到一个可行解然后再收紧MIPGap用热启动继续优化。这个两段式求解法在综合能源系统这种动辄上万个约束、几千个变量的模型里比一次性用严格精度求解快非常多而且你可以在调试阶段快速地拿到一个近似解来检查约束的合理性——等你确认模型逻辑没毛病再挂机等严格最优解省下的时间足够多刷好几篇论文。