搞微网优化的人尤其是做热电联供CHP方向的研究生和工程师基本都遇到过同一个窘境模型建得挺漂亮约束也写得齐全但一碰到实际项目就翻车。为什么因为真实场景里的光伏出力和负荷需求从来就不是一条平滑曲线它天生带着随机性和波动性。如果你把源荷都用确定性数值固定下来那优化结果拿到现场基本没法用——要么配出来的调度方案过于乐观要么在极端天气下直接崩溃。这个项目做的东西正是把源荷随机特征纳入热电联供微网的优化模型里用Matlab实现一套可复现的求解框架。这篇文章不是来给你念教材的而是从实际做的角度把整套东西拆开讲清楚源荷随机性怎么建模、场景怎么生成和削减、优化目标与约束怎么列、Matlab代码结构怎么搭、以及我踩过的那些坑。适合正在写微网优化方向论文的学生、需要做园区综合能源调度的工程师以及对这类问题有兴趣的自学者。你拿到的思路和框架可以直接套进自己的场景里去改。1. 整体设计思路为什么必须考虑源荷随机特征1.1 确定性优化模型的问题出在哪先看一个最简单的例子。假设你的微网里有一台燃气轮机、一组光伏板、一个电储能和一个蓄热罐负荷和光伏出力都用预测的确定性数值带入优化模型。这样可以快速求出一组调度方案但问题很明显光伏的预测出力在正午可能偏差很大如果实际出力比预测低20%燃气轮机就得在短时间内快速爬坡补电甚至被迫切除部分可中断负荷。如果实际热负荷比预测高蓄热罐又不够用就得启动备用锅炉运行成本直接飙升。换句话说确定性优化本质上是在“用一个确定的值去代替一个不确定的分布”。这种方式在规划阶段还能接受但在运行调度层面它忽略了两个关键事实第一光伏和风电出力本质上是随机过程受天气、云层、风速影响第二负荷也不是一成不变的用户行为、气温变化都会造成预测误差。不考虑这些得到的所谓“最优”调度方案往往是不经济的甚至是不可行的。1.2 随机优化让方案在多个场景下都“过得去”这个项目采用的核心思路是随机优化中的期望值模型。简单说就是生成足够多的源荷随机场景每个场景代表一种可能发生的“天气负荷”组合每个场景都有一定的概率然后优化目标变成让“所有场景下的期望运行成本”最低。举个例子一个晴天场景光伏出力高、负荷一般另一个阴天场景光伏出力低、负荷偏高。如果只按晴天场景调度阴天场景下就需要额外购电成本很高如果只按阴天场景调度晴天场景下燃气轮机就会过度运行浪费燃料。随机优化的好处在于它找到的是一个折中方案既保证大多数场景下成本可控又不会让极端场景全面崩盘。这里还要区分一下随机优化和鲁棒优化。鲁棒优化追求的是“最坏情况下也安全”结果往往非常保守经济性差随机优化则用概率加权的方式求期望最优在经济性和可靠性之间取得更合理的平衡。对于园区微网这类场景随机优化更符合实际需求这也是这个项目选择它的主要原因。1.3 求解框架从场景生成到最优调度整个项目的技术链路可以分成四段不确定性建模、场景生成与削减、建立优化模型、求解与分析。任何一段做不好结果都会出问题。不确定性建模是基础你需要用概率分布描述光伏、风电和负荷的随机特性。场景生成是把概率分布转化成一组“看得见、能计算”的数值样本比如用蒙特卡洛采样或拉丁超立方采样生成2000个场景。场景削减是为了降低计算量——2000个场景塞进一个MILP模型里求解时间很可能从分钟级变成小时级这在科研和工程上都不太能接受。削减到2030个代表性场景精度损失通常可以控制在误差范围内。之后才是建优化模型、列约束、调求解器。2. 核心细节源荷随机特征建模与场景削减实操2.1 风速、光照和负荷的概率分布怎么选这个项目的源荷随机建模部分核心是确定三个变量的概率分布参数。风速通常用两参数Weibull分布描述。形状参数k决定曲线的偏斜程度尺度参数c决定平均风速大小。以我手里一个常见案例为例某园区所在地区平均风速约6.5m/s我取k2.3、c8.5m/s这样生成的样本风速分布比较贴合实际记录。光伏出力受光照强度影响也受云层随机遮挡影响。光照强度常用Beta分布描述两个形状参数α和β需要根据历史辐照度的均值和方差估计。比如把光照强度归一化到01区间后取α3.5、β2.2可以模拟出晴天多、阴天少的辐照特征。太阳位置和光伏板倾角也要做成全天24小时曲线。负荷的不确定性相对温和一般用正态分布来描述预测误差。以预测值为均值标准差取预测值的5%10%。要注意电负荷和热负荷应该分别建模因为它们的波动特性不同——电负荷在一天内会有明显的早高峰和晚高峰热负荷则更受室外温度影响变化更平缓。这些参数直接用Matlab的概率分布对象去拟合历史数据用fitdist对历史风速、光照、负荷做参数估计比单纯拍脑袋定参数靠谱得多。2.2 场景生成蒙特卡洛采样和它的两个隐藏问题确定概率分布后就可以采样生成场景了。最基础的做法是蒙特卡洛采样对每个时段24小时的每个随机变量独立采样拼成一条24小时的“场景曲线”重复2000次得到一个2000×3×24的三维场景数组。但这里有两个问题我在实验中踩过必须提醒你。第一个是采样维度爆炸。24个时段乘以3个随机源意味着每次采样的是一个72维的向量。纯随机采样在高维空间里样本会非常稀疏生成的场景可能“看起来”不太像真实的源荷曲线——比如光伏出力在夜间突然出现非零值。解决办法是加限制条件采样时强制夜间光伏为0或者在采样后直接把这些不物理的样本剔除。第二个问题是场景相关性问题。光照和风速其实不是完全独立的——大风天往往云多、光照弱。如果忽略这一点生成的场景组合会高估某些情况的发生概率。比较好的做法是用Copula函数或者Cholesky分解引入相关系数矩阵把风速和光照的负相关性塞进采样过程。Matlab里直接用copularnd([Gaussian, rho, N])就能生成带相关性的样本很方便。2.3 场景削减从2000到20精度还不掉2000个场景什么都好就是求解太慢。一个带二进制变量的热电联供优化模型场景数翻倍求解时间可能翻好几倍2000个场景基本让你等到怀疑人生。场景削减算法的思路很简单从原始场景集合中找到一个数量更少的子集使得删除的场景与保留场景之间的“概率距离”最小。最经典的实现是同步回代消除法。通俗地讲就是反复比较所有场景对之间的欧氏距离每次把与其他场景最相似距离最近的那个场景删掉把它本来承担的概率转移到最相似的邻居身上直到场景数量降到目标值。我试过用K-means聚类做削减也能用但结果不如同步回代消除法稳定。聚类方法在某些场景下会把两个分布完全不同的场景归为一类导致削减后的场景集合失真。实验中我用同步回代消除法把2000个场景削减到20个再对比削减前后的期望目标函数值误差控制在3%以内效果非常理想。这段在Matlab里写起来也不复杂。核心步骤包括计算场景两两距离矩阵、循环查找最小距离、更新概率权重。用向量化操作可以大幅提速千万别在循环里一个个算距离2000×2000的距离矩阵也就几秒而已。3. 优化模型核心要点目标函数、约束与线性化处理3.1 设备建模和决策变量怎么定义热电联供微网里每个设备都得有对应的数学模型。燃气轮机是关键设备它同时产电和产热所以模型中它的出力包括电功率P_MT(t)和热功率H_MT(t)二者通过热电比η_CHP关联。烟气余热回收还有一个换热效率这些参数都从设备手册拿。模型里我还加入了光伏阵列、风力发电机、电储能、蓄热罐、燃气锅炉和电锅炉。每个设备的出力上下限、爬坡率限制、效率曲线都要体现在约束里。以电储能的SOC为例SOC(t1) SOC(t) (P_ch(t)·η_ch - P_dis(t)/η_dis)·Δt/Capacity这个状态递推方程是整个储能模型的核心。蓄热罐同理只是把电量换成热量把充电效率换成蓄热效率。决策变量包括各设备时段出力、储能和蓄热罐的充放状态、购售电功率、备用容量还有燃气轮机的启停状态。启停状态需要整数变量表示0/1二进制变量因此整个问题是个混合整数线性规划问题MILP。3.2 目标函数不只算燃料成本目标函数设计得合理不合理直接决定优化结果是否贴合实际。这个项目里总运行成本包括六项燃气轮机燃料成本通常按耗量特性曲线折算成二次函数或分段线性函数。购电成本跟电网分时电价有关峰谷电价差越大储能套利空间越大。启停成本如果燃气轮机频繁启停这部分成本会非常高模型里必须有惩罚项。弃风弃光惩罚为了防止优化结果为了保证燃气轮机利用率而白白扔掉清洁能源给弃电加一个惩罚系数。储能折旧成本每充放一次电池寿命就短一点这部分也要折算进去。最后是失负荷惩罚极端场景下切负荷的代价设为正常电价的几十倍防止模型为了省钱随便切负载。将六个成本加总求期望最小就是随机优化目标。在Matlab里用Yalmip写目标函数非常顺手支持分段函数、逻辑约束和二次项配合cplex或gurobi求解器效率很高。3.3 约束清单功率平衡之外容易漏掉的那些约束条件是一个优化模型的“边界”边界画错结果全错。电气功率平衡约束是刚性的电网购电功率 光伏出力 风电出力 燃气轮机发电 储能放电 电负荷 储能充电 电锅炉耗电。这个等式必须在每个场景每个时段都成立。热功率平衡约束同样重要燃气轮机余热 燃气锅炉产热 蓄热罐放热 电锅炉产热 热负荷 蓄热罐充热。我一开始把电锅炉放到了电负荷侧结算结果热平衡漏了一项调试时找到半夜才发现。设备运行约束别漏爬坡率。燃气轮机一分钟能爬多少兆瓦设备手册写得一清二楚但很多模型为了简化就省了这条。省掉的后果是优化出来的调度方案在现实中根本执行不了因为机组爬坡速度跟不上指令。备用容量约束是随机优化的重要组成部分。为了保证在源荷偏离预测时系统不至于崩溃模型会要求燃气轮机和储能预留一定的向上、向下备用容量。备用容量约束写成线性不等式的形式分别对应各个时段的旋转备用需求。小技巧我可以提前把设备参数、场景数据、电价参数都写到Excel表格里主程序统一读取而不是在代码里硬编码。这样以后换场景、换参数都不用改代码直接改表格就行迭代效率高很多。3.4 非线性项的线性化处理优化模型里很容易出现非线性项比如燃气轮机效率随出力变化导致的非线性成本、储能设备的充放电效率系数都是非线性的。MILP求解器处理不了非线性必须做线性化。常用手段有三个。第一是分段线性化把设备的耗量特性曲线分成若干段每段用斜率不同的线性函数逼近。第二是引入辅助二进制变量比如逻辑约束的转化、带最小启停时间约束的机组组合问题都需要这个技巧。第三是绝对值线性化处理类似储能SOC偏差的绝对值惩罚项时用两不等式等价转化即可。我遇到过的最考验功夫的线性化细节是处理储能同时充放电的问题。优化结果里可能会出现储能一边充电一边放电的“漏洞”白白浪费能量。解决方法是加入二进制变量y_storage(t)把充电功率和放电功率分别限定在y(t)·P_max和(1-y(t))·P_max以内二者互斥确保同一时刻只会充或只会放。这个约束虽然加上了二进制变量、增加了求解难度但必不可少。4. 实操过程Matlab代码框架与关键实现4.1 代码整体结构这个项目的Matlab代码我建议按模块化思路组织不要把所有逻辑堆在一个脚本里。我的目录结构大致是main.m入口脚本、数据读取模块、场景生成和削减模块、优化建模模块用Yalmip、结果展示模块、数据目录。main.m的主流程是读取数据 → 生成场景 → 场景削减 → 构建优化模型 → 求解 → 输出并绘图。每部分一个函数调试时定位问题非常方便。%% main.m - 热电联供微网随机优化主程序 % 清理环境 clear; clc; close all; % 初始化随机数种子保证实验可复现 rng(2024); % 读取基础数据设备参数、负荷曲线、分时电价等 data load_microgrid_data(data_microgrid.xlsx); % Step 1: 生成原始随机场景 2000个 scenarios_orig generate_scenarios(data, 2000); % Step 2: 场景削减保留20个典型场景 scenarios_red reduce_scenarios(scenarios_orig, 20); % Step 3: 构建并求解随机优化模型 result solve_stochastic_optimization(data, scenarios_red); % Step 4: 绘图与分析 plot_dispatch_result(result, data);4.2 场景生成与削减的关键代码逻辑场景生成部分核心是生成一天24小时的风电、光伏和负荷曲线。为了效率用向量化操作批量生成。function scenarios generate_scenarios(data, N) % 模拟N个随机场景 % data中包含Wekaill分布参数、Beta分布参数、负荷正态参数 T 24; % 时段数 scenarios.wind zeros(N, T); scenarios.pv zeros(N, T); scenarios.load zeros(N, T); % 仓位风速采样用wblrnd wind_speed wblrnd(data.wind_k, data.wind_c, N, T); % 转成风电出力简化V-I特性曲线 scenarios.wind wind_to_power(wind_speed, data.wind_rated, data.wind_cutin, data.wind_cutout); % 光伏光照强度满足Beta分布 midday最大值1 radiation betarnd(data.pv_alpha, data.pv_beta, N, T); radiation radiation .* repmat(data.pv_shape, N, 1); % 叠加日辐照形状 scenarios.pv radiation * data.pv_capacity * data.pv_eff; % 负荷预测值 正态偏差 scenarios.load repmat(data.load_predict, N, 1) .* (1 0.05 * randn(N, T)); end这里有个小坑需要注意直接用randn生成负荷误差可能出现负负荷得在生成后做一次截断处理。光伏出力的Beta分布在夜间时段也没有自动归零需要强行把夜间时段乘0。这两个小点不处理场景数据看起来就会很荒谬。场景削减部分的代码关键在同步回代消除的循环实现。核心逻辑如下function scenarios_red reduce_scenarios(scenarios_orig, target_count) % 基于同步回代消除法的场景削减 % 输入原始场景结构输出削减后的场景 N length(scenarios_orig.prob); dist_matrix compute_scenario_distance(scenarios_orig); % 未删除集合 active_set 1:N; while length(active_set) target_count % 寻找最小距离的场景对(i, j)其中i是要删除的候选 min_dist inf; del_idx 0; neighbor_idx 0; for i 1:length(active_set) for j 1:length(active_set) if j ~ i dist_matrix(active_set(i), active_set(j)) min_dist min_dist dist_matrix(active_set(i), active_set(j)); del_idx active_set(i); neighbor_idx active_set(j); end end end % 更新邻居概率 scenarios_orig.prob(neighbor_idx) scenarios_orig.prob(neighbor_idx) scenarios_orig.prob(del_idx); % 删除场景 active_set(active_set del_idx) []; end end双重循环写成这样很直观但削减2000个到20个时这个双重循环会非常慢每轮都要全量扫描一次距离矩阵。实际运行时建议用矩阵运算一次找出最小距离对我在本地实测中通过预计算距离矩阵和向量化查找降到了传统写法的1/10时间。4.3 求解Yalmip建模与求解器配置我用Yalmip做建模和求解。Yalmip的最大好处是建模和求解器分离你可以轻松切换cplex、gurobi或内置的sdpam来试验不同求解器的性能。对MILP来说我强烈建议用gurobi或cplex内置的线性求解器在大规模场景下会慢到让你怀疑人生。function result solve_stochastic_optimization(data, scenarios) % 定义变量 T 24; S length(scenarios.prob); P_mt sdpvar(S, T, full); % 燃气轮机发电 H_mt sdpvar(S, T, full); % 燃气轮机产热 P_pv_use sdpvar(S, T, full); % 光伏实际利用 P_wt_use sdpvar(S, T, full); % 风电实际利用 P_grid sdpvar(S, T, full); % 电网购电(正)/售电(负分开) P_buy sdpvar(S, T, full); P_sell sdpvar(S, T, full); P_es sdpvar(S, T, full); % 储能充放正充负放 SOC sdpvar(S, T1, full); ... % 目标函数期望运行成本 Objective 0; for s 1:S prob_s scenarios.prob(s); Objective Objective prob_s * (... sum(fuel_cost(P_mt(s,:))) ... sum(price_buy(s,:) .* P_buy(s,:)) - ... sum(price_sell(s,:) .* P_sell(s,:)) ... ... ); end optimize(constraints, Objective, opts); end约束部分最需要注意的变量维度要一致。我第一版代码里SOC设成了(S, T1)而充放电功率是(S, T)结果Yalmip报了一堆维度不匹配的报错。后来我统一约定所有时段变量用(S, T)状态递推单独在约束里处理SOC(t1) SOC(t) ...把SOC多出的T1列和P_es的T列对齐写约束时都要手动对齐时段坐标。求解器配置方面要设置求解时间上限和gap容忍度。默认的不设置会一直求解到全局最优但大模型下这个耗时不可控。我一般设置相对gap 0.01时间上限600秒对论文级应用完全够用。求解完成后顺手检查solvertime和gap信息判断求解质量。4.4 结果展示把结果可视化出来几乎是每个导师第一眼要看的东西。我用两张图一张是电功率平衡堆叠图另一张是热功率平衡堆叠图。X轴是24小时Y轴是功率堆叠面积图可以很直观地看出燃气轮机、光伏、风电、储能的出力构成。另一张重要的图是原始场景和削减后场景的对比图。把2000条半透明灰色场景曲线画成背景再把削减后的20条高亮色曲线叠加在上面一眼就能确认削减后的场景保留了原始分布特征不丢主要信息。这个图写论文的时候几乎是必放图。%% plot_dispatch_result.m - 结果可视化 plot(1:24, mean(scenarios_orig.load), k--, LineWidth, 2); hold on; for s 1:length(scenarios_red.prob) plot(1:24, scenarios_red.load(s,:), LineWidth, 1.2, ... Color, [0.2 0.4 0.8 0.6]); end xlabel(时段); ylabel(负荷/kW); legend(原始场景均值, 削减后场景);5. 常见问题与排查技巧实录5.1 为什么求解器给出的结果是0新手最常见的问题模型求解完所有决策变量全是0但目标函数值又不是0。这种情况十有八九是约束条件存在不可行求解器返回的是一个“可行的虚拟解”。排查技巧是先检查约束中最可能与数据冲突的部分。我的排查路径是先把所有设备上下限约束注释掉求解看是否正常出结果如果出来了就逐个添加约束每加一条求解一次哪个约束加进去后结果变0就是它的问题。实际操作中我遇到过两次一次是燃气轮机的热电比约束和总热功率平衡约束在某个场景下联立无解另一次是蓄热罐的SOC初始值和终值约束起了冲突。用二分法逐个约束排查很快就能定位。5.2 Yalmip的sdpvar维度不对错误信息一般类似“Unable to perform assignment because size of left side is...”。这个问题十有八九是你把变量当成了矩阵去切片或者某个约束里变量维度对不上。我的经验是定义变量时就用全维度的sdpvar(S, T, full)所有操作都保证在(S, T)这个二维平面上进行。如果某个约束里需要某一列直接用sdpvar(:, t)取列不要用squeeze转来转去那样容易产生维度静默变化的问题极其隐蔽难查。5.3 求解时间太长如果你的场景数量从20个加到100个就能明显觉得卡顿那大概率哪里写得不合理。比较常见的原因是用了过多的二进制变量或者有约束在循环里写了太多条。Yalmip里把矩阵约束一次性写进去远比循环写约束高效得多。以储能同时充放电约束为例你应该用矩阵形式的二进制变量y_storage一次性加上(1-y)相关的约束而不是在for循环里逐条添加。代码量少是次要的重要的是求解器拿到的模型更紧凑预求解阶段的工作量小很多计算效率会有明显的提升。场景削减参数也是一个平衡点。目标数量20个和50个求解时间可能差出35倍但目标函数值差异很小。我建议你做个简单的敏感性测试分别用10、20、30、50个场景求解画出一条目标值随场景数的变化曲线然后选那个“拐点”上面的场景数性价比最高。5.4 极端场景导致数值溢出有几次出现NaN或Inf的结果检查后发现是场景生成时随机采样到了极端参数。比如风速采样到了一个超过切出风速的值对应的风电出力函数除零或者负荷误差的高斯尾巴生成了负值再开方。应对措施是给所有随机变量加上物理下限和上限的截断。在采样后统一做一次clip将风速范围限定在[0, 40]m/s把负荷数值限制在非负范围内绝对不会影响概率分布特征但能彻底避免数值问题。这种“脏数据防护”在科研代码里价值极高虽然论文里不会写但实际操作中少了它会反复浪费时间。5.5 一天到底怎么分时段别全用1小时我在初始版本里用的是24个1小时时段后来项目扩展到对比不同调度周期时发现1小时粒度对于储能和蓄热罐的动态响应有些粗糙。改用15分钟粒度96个时段后储能的调度行为明显更精细高频波动也能捕捉到但场景生成和求解时间也会飙升。权衡下来如果你做的是规划层优化24点足够如果你做运行层优化建议用96点。这个选择要提前做不要模型建完之后再改否则所有维度都要重调浪费大量时间。6. 后续扩展与个人体会6.1 从期望值模型到鲁棒优化你的场景如果极端事件非常多比如台风多发区域的园区期望值模型在极端场景下会捉襟见肘。这时可以和鲁棒优化结合采用两阶段鲁棒优化第一阶段的决策先定下来第二阶段看最坏场景下的调整成本。这个扩展在Matlab里也比较好实现难点主要在分解算法上——一般用列与约束生成算法CCG处理两阶段问题。对论文产出来说这个点是很好的创新方向。6.2 把模型迁移到多微网联合调度热电联供单微网做扎实后下一步值得尝试的是多微网互联调度不同微网之间通过联络线交换功率形成共享备用和互济。源荷随机性在多微网场景下会被空间平滑效应部分抵消——这片区域的云来了隔壁区域的阳光可能还好。这个方向随机特征建模会更复杂但发CVPR级别的论文都有不少先例。6.3 我做这个项目的几点体会项目做到后面最深的体会是数学模型占三分工程实现占七分。源荷随机特征建模无论你写得多么精细最终能不能支撑调度决策还得看代码写得好不好、求解器调参到不到位、数据清洗不走样。很多学生卡在论文的仿真环节出不来其实都是代码工程能力拖了后腿。另外千万别迷信求解器默认参数。gurobi和cplex都有大量调参空间比如MIPFocus、Threads、TimeLimit、MIPGap针对不同模型规模做一些简单配置求解速度可能翻倍。我在实验中把Threads设成实际物理核数后大场景求解时间压缩了一半多。最后还要留意一个数学细节随机优化目标里各场景概率权重需要归一化有些用了削减前归一化的结果期望值偏差得很隐蔽。这个坑非常深你会发现在结果分析中目标函数值偏大或偏小怎么调约束都调不对根源就在概率和不是1。做场景削减之后记得重新检查一遍所有场景概率之和必须严格等于1。这个小检查花你十秒钟但能省掉至少半天的迷惑排查时间。