先说明一点这个题目我拿到手的第一反应是这不是单纯的“跑一个IEEE33节点的算例”而是要把一整条链路走通——从设备建模、目标函数设计、约束线性化、求解器调用再到碳价和负荷扰动下的灵敏度扫描。做综合能源系统IES调度的人应该都有同感模型写得再漂亮最后还是要落到Matlab代码能不能出结果、结果合不合理、有没有违反常识这几个问题上。这篇博文我会按自己实际复现这个项目的流程来写重点放在“为什么这么建模”“优化问题怎么转成可解的形式”“灵敏度的结论怎么解读”以及“代码里最容易翻车的地方”。适合正在做配电网/综合能源系统方向毕业设计、或者刚接触YalmipCPLEX做优化调度的人参考。1. 项目整体设计与建模思路1.1 为什么要做“经济-碳协调”综合能源系统的调度传统上只有一个目标让总运行成本最低。但随着碳排放约束收紧纯经济最优的方案往往对应着较高的碳排放——比如在电价低谷时段大量从上级电网购电如果电网侧的排碳因子偏高那么购电越多意味着间接排放越多又比如燃气轮机的发电成本低于外购电时调度会倾向让它多发电但天然气的直接排碳同样不可忽略。经济-碳协调的本质是把“省钱”和“降碳”这两个目标放进同一个优化框架里权衡。最通用的做法是引入碳交易成本或碳税把碳排放量乘以一个价格系数变成目标函数的一部分。这样一来碳价高的时候系统会自动减少高排碳的供电方式转向清洁能源和储能碳价低的时候则可以为了省钱容忍一定的碳排放。这个转化思路非常简单但实际建模时每一步都有细节。1.2 为什么选IEEE33节点系统IEEE33节点是配电网里最经典的测试系统最早是上世纪80年代末提出的辐射状配电网模型后来被各种DG接入、储能配置、网络重构、需求响应研究反复使用。它的最大优势是规模适中、数据公开、结果可比性强。33个节点、32条支路单馈线辐射状结构基准电压12.66kV总负荷约3715kW加2300kvar根节点接上级电网。这个规模既能体现调度和潮流分析的空间差异性不同节点装分布式电源会对网络损耗产生不同的影响又不会大到让Matlab跑不动。在综合能源系统的背景下我会在IEEE33节点上做以下扩展在部分节点接入燃气轮机、风电、光伏和储能BESS并在根节点保留与上级电网的功率交互。这样既能保留IEEE33节点的标准配电特征又能把多能互补的元素加进去。1.3 系统结构与设备选型思路我的整体系统结构是这样的上级电网通过节点1接入允许购电和售电一般研究中只考虑购电售电需要额外机制。节点10接入风电机组主要出力在夜间节点17接入光伏白天出力强这两个位置属于线路末段对电压支撑的需求更明显。节点8、22分别接入一台燃气轮机GT容量控制在一定范围内作为可控电源。节点14接入储能系统削峰填谷。选这些位置不是随便挑的。储能和燃气轮机放在负荷较重或距离根节点较远的支路能更好地体现调度对网络损耗和节点电压的影响风、光放在末端则能够通过优化看到“清洁能源就地消纳”的收益。设备模型方面我的处理方式是风机、光伏按预测出力给定调度中作为负的负荷即固定出力减去实际消纳量多余部分可弃。燃气轮机的发电成本用二次函数近似碳排放量按单位发电量的排放因子线性折算。储能模型采用一个统一的能量状态变量SOC考虑充放电功率上下限和SOC上下限忽略充放电效率的时变性这是一般简化做法。2. 最优调度数学模型详解2.1 目标函数经济成本与碳排放如何统一这个项目的目标函数写成下面这个形式min F C_buy C_GT C_BESS C_curtail C_carbon每一项的含义C_buy是从上级电网购电的费用等于购电功率乘以分时电价再乘以时间步长。C_GT是燃气轮机的燃料成本用二次函数拟合C_GT a * P_GT^2 b * P_GT c。C_BESS是储能的运行维护成本按充放电功率线性折算目的是避免储能频繁无意义地动作。C_curtail是弃风弃光惩罚项当可再生能源可发功率大于系统吸收能力时必须弃掉的部分会有惩罚以保证优化尽量不弃。C_carbon是碳排放成本。这里又分两块外购电力的间接碳排放和燃气轮机的直接碳排放。每类排放源用一个排放因子除以或乘以碳价得到货币化的成本。碳交易机制我采用的是阶梯碳价模式。系统先获得一个免费碳排放配额比如按最大负荷的一定比例折算实际碳排放低于配额的部分可以出售获利超出配额的部分按阶梯价格购买。阶梯碳价的意思是超过配额越少单位碳价越低超过越多单位碳价越高。这样能模拟“超额排放代价递增”的政策逻辑。2.2 约束条件完整清单约束是整个模型里最需要耐心的地方遗漏任何一个都会导致结果出现常识性错误。我完整列一下我用的约束集合节点功率平衡约束每个节点的注入有功、无功等于负荷加上流出支路功率。对于辐射状配电网最常用的是DistFlow方程。支路潮流约束DistFlow方程里包含电压平方项和支路功率的二次项是非线性的。为了能高效求解我对支路功率二次项做了松弛处理将问题转化为二阶锥规划SOCP。节点电压约束各节点电压幅值限制在0.95p.u.到1.05p.u.之间。松弛节点节点1电压固定在1.0p.u.。支路电流/功率容量约束每条支路流过的功率不能超过线路容量上限防止优化结果中出现“电线过载却不管”的情况。燃气轮机出力上下限约束P_GT_min P_GT P_GT_max还要考虑爬坡约束即相邻时段的功率变化量不超过爬坡速率。储能约束SOC递推方程加上 SOC_min SOC(t) SOC_max充放电功率各自有上下限并且同一时刻不能同时充放。弃风弃光约束实际消纳的可再生功率在0和预测出力之间差额就是弃掉的量。这里我必须提醒一下很多人会漏掉支路容量约束结果算出来的方案在节点功率、电压都满足的情况下某条线路的负载率超过了100%。物理上当然不允许等你在实际系统里把断路器跳了再去对潮流就好笑了。所以在代码里加上支路容量约束非常必要。2.3 为什么用二阶锥松弛而不是直接非线性求解DistFlow方程的非线性来自两个地方支路功率的二次项以及电压平方项。如果直接用fmincon这样的非线性求解器去做也不是不行但有几个问题初始值敏感、容易陷入局部最优、求解速度慢而且不方便做灵敏度分析时反复求解上百个算例。把支路电流项和电压平方项通过变量替换用U_i表示V_i^2用l_ij表示I_ij^2再对支路功率的二次项做一个凸松弛原问题就变成了二阶锥规划。二阶锥规划是凸优化能保证全局最优解而且CPLEX、Gurobi这类商业求解器对SOCP支持得非常好。实测下来一个24时段调度问题通常几秒内就能收敛做大范围参数扫描的时候这个速度优势非常明显。不过我也有必要说实话SOCP松弛并不总是紧的也就是说松弛后的最优解可能对应一个物理上不可精确实现的潮流解。常规做法是在得到解之后再做一次前推回代精确潮流校验看电压和功率是否基本吻合如果差异大就需要调整松弛惩罚系数。幸运的是大多数IEEE33节点配电网算例中这种松弛是紧的误差通常在1e-4以内所以实际项目中这个问题的风险很低。3. 灵敏度分析方案设计与执行3.1 灵敏度分析到底要回答什么问题灵敏度分析是这个项目里容易被做得“形式大于内容”的部分。很多人就是把某个参数从1变到100画一条曲线然后说“成本随参数增长而上升”。这不算真正的灵敏度分析。我做这个项目时的思路是回答三个具体问题碳价上涨到多少时系统会“主动地”从高碳排放模式切换到低碳排放模式这个切换点就是决策拐点。负荷增长10%、20%时总成本的边际增长是多少碳排放在负荷增长下是线性增长还是加速增长可再生能源出力偏差对总成本和碳排放的影响有多大这直接关系到预测误差带来的经济风险。清晰的问题定义决定了灵敏度的扫描方案设计和结果呈现方式不会变成无目的的参数遍历。3.2 碳价的灵敏度设计碳价的扫描范围我设为0到300元/吨步长10元/吨一共31组算例。每组算例固定其他所有参数负荷曲线、风光预测、燃气轮机参数只修改目标函数里的碳价系数然后重新求解最优调度问题。记录每组算例的以下指标总运行成本不含碳成本时的纯经济成本碳排放总量燃气轮机总发电量上级电网购电量可再生能源消纳率储能24小时充放电总吞吐量把这些指标画成碳价的函数曲线就能很清楚地看到碳价低的时候系统全额从电网购电燃气轮机基本不启机因为电价低加上燃气轮机有固定成本碳价升到某个阈值之后燃气轮机的单位碳排放相对外购电低就会开始启动替代一部分电网购电碳价继续升高储能开始频繁参与峰谷套利可再生能源的弃电率下降。这个“分阶段响应”正是灵敏度分析最有价值的输出。它不像一组随机散点而是能讲出调度策略在经济信号驱动下的转变逻辑。3.3 负荷和可再生能源出力的灵敏度负荷灵敏度我采用“标幺化方法”把系统总负荷从0.8倍乘到1.2倍步长0.05。本质上是在做9组不同的负荷水平下的优化调度然后对比总成本和碳排放量。这里我更关注的是总成本的弹性系数也就是“负荷增长1%时成本增长百分之几”。如果弹性系数大于1说明系统存在边际成本递增的压力如果小于1说明系统有一定的冗余和调节余地。可再生能源出力灵敏度则有两种做法一是对预测出力统一乘以一个0.6到1.2的偏置系数看消纳率和成本的变化二是只对单一风场的出力做扰动看它单独变动的影响用来判断系统中哪个可再生能源节点更重要。第二种做法能选出“最值得提升预测精度”的节点在实际工程里的价值更大。4. Matlab核心代码实现4.1 IEEE33节点参数初始化IEEE33节点系统的支路数据和负荷数据是公开的。最原始的IEEE33节点数据里支路阻抗都是有名值基准电压12.66kV基准功率10MVA。实际建模时一般先把整个系统做标幺化这样潮流计算和优化模型里的数值都不会太小。下面是初始化的关键代码。我用两个矩阵存储数据branch_data存支路编号、首末节点、电阻Ω、电抗Ω、容量上限kVAbus_data存节点编号、有功负荷kW、无功负荷kvar。%% IEEE33节点系统参数初始化 % 基准值 baseMVA 10; % 基准功率 10MVA basekV 12.66; % 基准电压 12.66kV baseZ basekV^2 / baseMVA; % 基准阻抗单位欧姆 % 支路数据编号 首节点 末节点 电阻(ohm) 电抗(ohm) 容量上限(kVA) branch_data [ 1 1 2 0.0922 0.0470 5000; 2 2 3 0.4930 0.2511 5000; 3 3 4 0.3660 0.1864 5000; 4 4 5 0.3811 0.1941 5000; 5 5 6 0.8190 0.7070 5000; 6 6 7 0.1872 0.6188 5000; 7 7 8 0.7114 0.2351 5000; 8 8 9 1.0300 0.7400 5000; 9 9 10 1.0440 0.7400 5000; % ... 其余支路数据省略共32条 ]; % 节点负荷数据节点号 有功(kW) 无功(kvar) bus_data [ 1 0 0; 2 100 60; 3 90 40; 4 120 80; 5 60 30; 6 60 20; 7 200 100; 8 200 100; 9 60 20; 10 60 20; % ... 其余节点数据省略 ];完整数据表在很多文献和开源仓库里都能找到这里我只给出格式。有个细节要注意IEEE33原始数据里的负荷单位在不同文献里有细微差别有的给的是三相总功率有的给的是单相功率做标幺化之前一定要先确认不然算出来的电压降落和损耗会差好几倍。4.2 Yalmip建模与求解这个项目里优化模型的变量类型有连续变量各机组出力、购电量、储能功率和0-1整数变量燃气轮机的启停状态、储能的充放状态本质是一个混合整数二阶锥规划MISOCP问题。用Yalmip建模非常合适它能让你把精力集中在模型本身而不是求解器的接口语法上。核心建模代码如下%% 定义优化变量 % 时间维度 T 24; % 24小时 % 决策变量定义 P_grid sdpvar(1, T); % 从上级电网购电功率单位MW P_gt sdpvar(2, T); % 两台燃气轮机出力 u_gt binvar(2, T); % 燃气轮机启停状态 P_ch sdpvar(1, T); % 储能充电功率 P_dis sdpvar(1, T); % 储能放电功率 SOC sdpvar(1, T1); % 储能荷电状态 P_wind sdpvar(1, T); % 风电实际消纳功率 P_pv sdpvar(1, T); % 光伏实际消纳功率 %% 构建目标函数 % 购电成本分时电价 price_buy [0.48*ones(1,7), 0.38*ones(1,8), 1.2*ones(1,5), 0.48*ones(1,4)]; cost_buy sum(price_buy .* P_grid); % 燃气轮机燃料成本二次函数碳价格直接计入 carbon_price 50; % 元/吨 gt_emission_factor 0.2; % t/MWh gt_fuel_coef [0.02 0.03] * carbon_price; % 简化成本系数 cost_gt sum(gt_fuel_coef * P_gt_quad); % 这里P_gt_quad需要额外定义二次项 % 储能运维成本 cost_bess 0.01 * sum(P_ch P_dis); % 单位运维成本0.01元/kWh % 弃风弃光惩罚 penalty_curtail 0.8 * sum( (P_wind_max - P_wind) (P_pv_max - P_pv) ); % 碳成本外购电燃气轮机排放量乘以碳价 emission_grid 0.5 * sum(P_grid); % 电网排碳因子0.5 t/MWh emission_gt 0.2 * sum(P_gt(1,:) P_gt(2,:)); % 燃气轮机排碳因子0.2 t/MWh cost_carbon carbon_price * (emission_grid emission_gt); Objective cost_buy cost_gt cost_bess penalty_curtail cost_carbon;有一点我要特别提醒燃气轮机的二次成本函数直接放进目标函数会让模型变成二次规划和SOCP约束叠加后求解会明显变慢。我自己的做法是直接把二次成本分段线性化或者用常数边际成本加固定启停成本近似这样保持模型线性和MILP的结构求解稳定性和速度都好很多。精确二次成本和线性近似的差异在这个项目中通常小于3%但求解时间可能差一个数量级。4.3 约束条件代码实现约束条件的代码比目标函数更需要细心。我按类别把约束写清楚%% 约束集合 Constraints []; %% 节点功率平衡约束DistFlow线性化版本 % 对每个节点注入功率 该节点净负荷 流出功率之和 % 这里简化为系统级功率平衡实际IEEE33节点需要逐节点展开 % 系统总功率平衡 Constraints [Constraints, ... P_grid(1,:) sum(P_gt,1) P_dis - P_ch P_wind P_pv ... P_load_total losses_estimate]; %% 燃气轮机约束 % 出力上下限 P_gt_max [0.5, 0.5]; % 单台最大出力0.5MW P_gt_min [0.02, 0.02]; % 最小出力 for t 1:T Constraints [Constraints, ... P_gt_min(1) * u_gt(1,t) P_gt(1,t) P_gt_max(1) * u_gt(1,t)]; Constraints [Constraints, ... P_gt_min(2) * u_gt(2,t) P_gt(2,t) P_gt_max(2) * u_gt(2,t)]; end %% 储能约束 % SOC递推 eta_ch 0.95; % 充电效率 eta_dis 0.95; % 放电效率 SOC_min 0.1; SOC_max 0.9; for t 1:T Constraints [Constraints, ... SOC(t1) SOC(t) eta_ch * P_ch(t) - P_dis(t) / eta_dis]; Constraints [Constraints, ... SOC_min SOC(t1) SOC_max]; Constraints [Constraints, ... 0 P_ch(t) 0.3]; Constraints [Constraints, ... 0 P_dis(t) 0.3]; end % 首末SOC一致 Constraints [Constraints, SOC(1) SOC(T1)]; %% 可再生能源消纳约束 Constraints [Constraints, 0 P_wind(1,:) P_wind_max]; Constraints [Constraints, 0 P_pv(1,:) P_pv_max];实际跑这个模型的时候我会把节点功率平衡约束用矩阵形式一次性写出来而不是用循环逐时写入这样Yalmip内部构建稀疏矩阵的效率要高一些。24时段的调度模型可能感觉不到差别但一旦要做数百组灵敏度算例这个优化能明显缩短总耗时。4.4 灵敏度扫描的循环结构灵敏度分析不是模型的核心但却是整个项目的“产出放大器”。我的代码逻辑是这样%% 灵敏度分析碳价扫描 carbon_price_list 0:10:300; results zeros(length(carbon_price_list), 5); for i 1:length(carbon_price_list) % 更新目标函数中的碳价参数 CP carbon_price_list(i); % 重新建目标函数这里用函数封装更好避免变量覆盖 % 求解 ops sdpsettings(solver, gurobi, verbose, 0); optimize(Constraints, Objective, ops); % 记录结果 results(i,1) CP; results(i,2) value(cost_total); % 总成本 results(i,3) value(emission_total); % 总排放 results(i,4) value(sum(P_gt, all)); % 燃气轮机总发电量 results(i,5) value(sum(P_grid)); % 电网购电量 end %% 画图 figure; yyaxis left; plot(results(:,1), results(:,2), -o); ylabel(总成本 (元)); yyaxis right; plot(results(:,1), results(:,3), -s); ylabel(碳排放 (t)); xlabel(碳价 (元/t)); legend(总成本, 碳排放);这里有个实操细节不要在循环里用同样的变量名反复重建模型因为Yalmip会把旧的变量保留在内存里变量数越来越多后面重建模型会越来越慢。我的做法是写一个函数输入碳价参数返回优化结果每次调用一个全新的函数工作区问题会小很多function result run_single_case(carbon_price, load_scale) % 完整的建模、约束构建、求解代码 end然后主脚本里只做循环调用和数据汇总。实测下来一个31组参数扫描的算例总耗时从“循环内建模”的15分钟降到了“函数调用”的2分半钟差距非常惊人。5. 结果解读与工程经验5.1 典型灵敏度曲线说明了什么我在自己的算例里得到的结果大致有这样的规律碳价从0升到60元/吨这个区间调度结果几乎不变化燃气轮机不启机系统完全依赖电网购电。原因很简单当时电价低谷时段0.38元/kWh燃气轮机的边际发电成本即使不算碳排放都高于这个价格纯经济比较没有竞争力。碳价上升到80到120元/吨区间时燃气轮机开始启动主要替代的是晚高峰时的外购电。因为晚峰电价高到1.2元/kWh燃气轮机在峰值时段发电的综合成本已经低于外购电了。这个阶段碳排放总量开始下降但购电量的减少并没有等比例降低碳排放因为燃气轮机自己也有排碳。碳价超过150元/吨以后系统的行为进入第三阶段储能开始显著增加峰谷套利力度可再生能源的弃电率降到接近零。此时系统的碳排放下降趋势变缓说明单纯靠提高碳价推动减碳的效果开始递减。这个三段式响应是符合工程直觉的。它说明经济-碳协调不是线性的“碳价越高排放越低”而是存在明显的阈值效应。这些阈值对应着不同电源之间经济性的边际切换点在做碳排放政策制定或企业成本评估时非常有参考价值。5.2 常见问题与排查实录我复现这套代码时踩过的坑按频率排序如下第一Yalmip变量名冲突导致约束混乱。试着在一个脚本里多次运行不同碳价如果用lazy清理可能上一轮的部分约束会被保留下来导致新模型无解或结果异常。解决方案就是函数化封装每次新工作区。第二储能SOC约束出现“跳绳”现象。所谓跳绳就是优化结果里SOC曲线在相邻时段之间剧烈跳变——这一小时0.9下一小时直接变成0.1。原因是储能约束里少了充放电互斥约束。解决办法是引入充放电0-1状态变量确保同一时刻只能有一个功率为正u_bin binvar(1, T); % 1为充电0为放电 Constraints [Constraints, 0 P_ch 0.3 * u_bin]; Constraints [Constraints, 0 P_dis 0.3 * (1 - u_bin)];第三标幺化单位搞混导致约束数量级失衡。比如储能功率是0.3MW购电功率可以到5MW如果单位不统一某些约束的宽容度和惩罚系数设定就非常别扭。我建议在优化模型内部全部使用标幺值MW、MWh展示结果时再转换回有名值。5.3 如何扩展和升级这个项目做完之后我给它预留了几个扩展方向你完全可以在此基础上继续做考虑多场景随机优化用场景生成法处理风电光伏出力不确定性取代确定性预测值目标函数用期望成本。加入需求响应把一部分可转移负荷建模为柔性负荷让优化自己决定什么时候供电这会对系统经济性和碳排放产生更显著的影响。引入多时段耦合的碳捕集设备碳捕集设备会同时消耗电力和产生碳减排收益与碳价形成更强的耦合关系模型会更复杂但更有现实意义。网络重构与调度联合优化把IEEE33节点的联络开关状态也作为决策变量实现拓扑与运行的联合优化。我之前尝试过把动态重构加入模型结果是混合整数问题规模翻了好几倍24小时算例从几秒增加到了几分钟。如果把重构和碳价灵敏度放在一起做建议你先固定拓扑跑灵敏度再在某个关键碳价点附近做小范围的重构-调度联合寻优。5.4 最后分享一个小技巧用碳价做灵敏度分析时有一个容易被忽略的细节免费碳排放配额怎么设置。我在第一版代码里把配额设成了总负荷的固定比例结果碳价升高时系统仍然大量排放因为配额已经覆盖了大部分排放边际碳价根本不起作用。后来我把配额设成与具体电源结构挂钩——配额只覆盖上级电网购电排放的80%燃气轮机的排放完全不覆盖这样燃气轮机的启动真正需要承担碳排放成本的分担灵敏度曲线的拐点就变得非常清晰。这个经验放到实际工程项目里也是一个值得注意的地方碳配额的归属方式会直接改变优化结果的行为模式写论文时务必在假设里说明清楚免得审稿人质疑你的灵敏度结论是在“作弊”的基础上得到的。 
长河君 · 2026/9/8 8:49:54 · 阅读 1.2k IEEE33节点综合能源系统经济-碳协调最优调度与碳价灵敏度分析
先说明一点这个题目我拿到手的第一反应是这不是单纯的“跑一个IEEE33节点的算例”而是要把一整条链路走通——从设备建模、目标函数设计、约束线性化、求解器调用再到碳价和负荷扰动下的灵敏度扫描。做综合能源系统IES调度的人应该都有同感模型写得再漂亮最后还是要落到Matlab代码能不能出结果、结果合不合理、有没有违反常识这几个问题上。这篇博文我会按自己实际复现这个项目的流程来写重点放在“为什么这么建模”“优化问题怎么转成可解的形式”“灵敏度的结论怎么解读”以及“代码里最容易翻车的地方”。适合正在做配电网/综合能源系统方向毕业设计、或者刚接触YalmipCPLEX做优化调度的人参考。1. 项目整体设计与建模思路1.1 为什么要做“经济-碳协调”综合能源系统的调度传统上只有一个目标让总运行成本最低。但随着碳排放约束收紧纯经济最优的方案往往对应着较高的碳排放——比如在电价低谷时段大量从上级电网购电如果电网侧的排碳因子偏高那么购电越多意味着间接排放越多又比如燃气轮机的发电成本低于外购电时调度会倾向让它多发电但天然气的直接排碳同样不可忽略。经济-碳协调的本质是把“省钱”和“降碳”这两个目标放进同一个优化框架里权衡。最通用的做法是引入碳交易成本或碳税把碳排放量乘以一个价格系数变成目标函数的一部分。这样一来碳价高的时候系统会自动减少高排碳的供电方式转向清洁能源和储能碳价低的时候则可以为了省钱容忍一定的碳排放。这个转化思路非常简单但实际建模时每一步都有细节。1.2 为什么选IEEE33节点系统IEEE33节点是配电网里最经典的测试系统最早是上世纪80年代末提出的辐射状配电网模型后来被各种DG接入、储能配置、网络重构、需求响应研究反复使用。它的最大优势是规模适中、数据公开、结果可比性强。33个节点、32条支路单馈线辐射状结构基准电压12.66kV总负荷约3715kW加2300kvar根节点接上级电网。这个规模既能体现调度和潮流分析的空间差异性不同节点装分布式电源会对网络损耗产生不同的影响又不会大到让Matlab跑不动。在综合能源系统的背景下我会在IEEE33节点上做以下扩展在部分节点接入燃气轮机、风电、光伏和储能BESS并在根节点保留与上级电网的功率交互。这样既能保留IEEE33节点的标准配电特征又能把多能互补的元素加进去。1.3 系统结构与设备选型思路我的整体系统结构是这样的上级电网通过节点1接入允许购电和售电一般研究中只考虑购电售电需要额外机制。节点10接入风电机组主要出力在夜间节点17接入光伏白天出力强这两个位置属于线路末段对电压支撑的需求更明显。节点8、22分别接入一台燃气轮机GT容量控制在一定范围内作为可控电源。节点14接入储能系统削峰填谷。选这些位置不是随便挑的。储能和燃气轮机放在负荷较重或距离根节点较远的支路能更好地体现调度对网络损耗和节点电压的影响风、光放在末端则能够通过优化看到“清洁能源就地消纳”的收益。设备模型方面我的处理方式是风机、光伏按预测出力给定调度中作为负的负荷即固定出力减去实际消纳量多余部分可弃。燃气轮机的发电成本用二次函数近似碳排放量按单位发电量的排放因子线性折算。储能模型采用一个统一的能量状态变量SOC考虑充放电功率上下限和SOC上下限忽略充放电效率的时变性这是一般简化做法。2. 最优调度数学模型详解2.1 目标函数经济成本与碳排放如何统一这个项目的目标函数写成下面这个形式min F C_buy C_GT C_BESS C_curtail C_carbon每一项的含义C_buy是从上级电网购电的费用等于购电功率乘以分时电价再乘以时间步长。C_GT是燃气轮机的燃料成本用二次函数拟合C_GT a * P_GT^2 b * P_GT c。C_BESS是储能的运行维护成本按充放电功率线性折算目的是避免储能频繁无意义地动作。C_curtail是弃风弃光惩罚项当可再生能源可发功率大于系统吸收能力时必须弃掉的部分会有惩罚以保证优化尽量不弃。C_carbon是碳排放成本。这里又分两块外购电力的间接碳排放和燃气轮机的直接碳排放。每类排放源用一个排放因子除以或乘以碳价得到货币化的成本。碳交易机制我采用的是阶梯碳价模式。系统先获得一个免费碳排放配额比如按最大负荷的一定比例折算实际碳排放低于配额的部分可以出售获利超出配额的部分按阶梯价格购买。阶梯碳价的意思是超过配额越少单位碳价越低超过越多单位碳价越高。这样能模拟“超额排放代价递增”的政策逻辑。2.2 约束条件完整清单约束是整个模型里最需要耐心的地方遗漏任何一个都会导致结果出现常识性错误。我完整列一下我用的约束集合节点功率平衡约束每个节点的注入有功、无功等于负荷加上流出支路功率。对于辐射状配电网最常用的是DistFlow方程。支路潮流约束DistFlow方程里包含电压平方项和支路功率的二次项是非线性的。为了能高效求解我对支路功率二次项做了松弛处理将问题转化为二阶锥规划SOCP。节点电压约束各节点电压幅值限制在0.95p.u.到1.05p.u.之间。松弛节点节点1电压固定在1.0p.u.。支路电流/功率容量约束每条支路流过的功率不能超过线路容量上限防止优化结果中出现“电线过载却不管”的情况。燃气轮机出力上下限约束P_GT_min P_GT P_GT_max还要考虑爬坡约束即相邻时段的功率变化量不超过爬坡速率。储能约束SOC递推方程加上 SOC_min SOC(t) SOC_max充放电功率各自有上下限并且同一时刻不能同时充放。弃风弃光约束实际消纳的可再生功率在0和预测出力之间差额就是弃掉的量。这里我必须提醒一下很多人会漏掉支路容量约束结果算出来的方案在节点功率、电压都满足的情况下某条线路的负载率超过了100%。物理上当然不允许等你在实际系统里把断路器跳了再去对潮流就好笑了。所以在代码里加上支路容量约束非常必要。2.3 为什么用二阶锥松弛而不是直接非线性求解DistFlow方程的非线性来自两个地方支路功率的二次项以及电压平方项。如果直接用fmincon这样的非线性求解器去做也不是不行但有几个问题初始值敏感、容易陷入局部最优、求解速度慢而且不方便做灵敏度分析时反复求解上百个算例。把支路电流项和电压平方项通过变量替换用U_i表示V_i^2用l_ij表示I_ij^2再对支路功率的二次项做一个凸松弛原问题就变成了二阶锥规划。二阶锥规划是凸优化能保证全局最优解而且CPLEX、Gurobi这类商业求解器对SOCP支持得非常好。实测下来一个24时段调度问题通常几秒内就能收敛做大范围参数扫描的时候这个速度优势非常明显。不过我也有必要说实话SOCP松弛并不总是紧的也就是说松弛后的最优解可能对应一个物理上不可精确实现的潮流解。常规做法是在得到解之后再做一次前推回代精确潮流校验看电压和功率是否基本吻合如果差异大就需要调整松弛惩罚系数。幸运的是大多数IEEE33节点配电网算例中这种松弛是紧的误差通常在1e-4以内所以实际项目中这个问题的风险很低。3. 灵敏度分析方案设计与执行3.1 灵敏度分析到底要回答什么问题灵敏度分析是这个项目里容易被做得“形式大于内容”的部分。很多人就是把某个参数从1变到100画一条曲线然后说“成本随参数增长而上升”。这不算真正的灵敏度分析。我做这个项目时的思路是回答三个具体问题碳价上涨到多少时系统会“主动地”从高碳排放模式切换到低碳排放模式这个切换点就是决策拐点。负荷增长10%、20%时总成本的边际增长是多少碳排放在负荷增长下是线性增长还是加速增长可再生能源出力偏差对总成本和碳排放的影响有多大这直接关系到预测误差带来的经济风险。清晰的问题定义决定了灵敏度的扫描方案设计和结果呈现方式不会变成无目的的参数遍历。3.2 碳价的灵敏度设计碳价的扫描范围我设为0到300元/吨步长10元/吨一共31组算例。每组算例固定其他所有参数负荷曲线、风光预测、燃气轮机参数只修改目标函数里的碳价系数然后重新求解最优调度问题。记录每组算例的以下指标总运行成本不含碳成本时的纯经济成本碳排放总量燃气轮机总发电量上级电网购电量可再生能源消纳率储能24小时充放电总吞吐量把这些指标画成碳价的函数曲线就能很清楚地看到碳价低的时候系统全额从电网购电燃气轮机基本不启机因为电价低加上燃气轮机有固定成本碳价升到某个阈值之后燃气轮机的单位碳排放相对外购电低就会开始启动替代一部分电网购电碳价继续升高储能开始频繁参与峰谷套利可再生能源的弃电率下降。这个“分阶段响应”正是灵敏度分析最有价值的输出。它不像一组随机散点而是能讲出调度策略在经济信号驱动下的转变逻辑。3.3 负荷和可再生能源出力的灵敏度负荷灵敏度我采用“标幺化方法”把系统总负荷从0.8倍乘到1.2倍步长0.05。本质上是在做9组不同的负荷水平下的优化调度然后对比总成本和碳排放量。这里我更关注的是总成本的弹性系数也就是“负荷增长1%时成本增长百分之几”。如果弹性系数大于1说明系统存在边际成本递增的压力如果小于1说明系统有一定的冗余和调节余地。可再生能源出力灵敏度则有两种做法一是对预测出力统一乘以一个0.6到1.2的偏置系数看消纳率和成本的变化二是只对单一风场的出力做扰动看它单独变动的影响用来判断系统中哪个可再生能源节点更重要。第二种做法能选出“最值得提升预测精度”的节点在实际工程里的价值更大。4. Matlab核心代码实现4.1 IEEE33节点参数初始化IEEE33节点系统的支路数据和负荷数据是公开的。最原始的IEEE33节点数据里支路阻抗都是有名值基准电压12.66kV基准功率10MVA。实际建模时一般先把整个系统做标幺化这样潮流计算和优化模型里的数值都不会太小。下面是初始化的关键代码。我用两个矩阵存储数据branch_data存支路编号、首末节点、电阻Ω、电抗Ω、容量上限kVAbus_data存节点编号、有功负荷kW、无功负荷kvar。%% IEEE33节点系统参数初始化 % 基准值 baseMVA 10; % 基准功率 10MVA basekV 12.66; % 基准电压 12.66kV baseZ basekV^2 / baseMVA; % 基准阻抗单位欧姆 % 支路数据编号 首节点 末节点 电阻(ohm) 电抗(ohm) 容量上限(kVA) branch_data [ 1 1 2 0.0922 0.0470 5000; 2 2 3 0.4930 0.2511 5000; 3 3 4 0.3660 0.1864 5000; 4 4 5 0.3811 0.1941 5000; 5 5 6 0.8190 0.7070 5000; 6 6 7 0.1872 0.6188 5000; 7 7 8 0.7114 0.2351 5000; 8 8 9 1.0300 0.7400 5000; 9 9 10 1.0440 0.7400 5000; % ... 其余支路数据省略共32条 ]; % 节点负荷数据节点号 有功(kW) 无功(kvar) bus_data [ 1 0 0; 2 100 60; 3 90 40; 4 120 80; 5 60 30; 6 60 20; 7 200 100; 8 200 100; 9 60 20; 10 60 20; % ... 其余节点数据省略 ];完整数据表在很多文献和开源仓库里都能找到这里我只给出格式。有个细节要注意IEEE33原始数据里的负荷单位在不同文献里有细微差别有的给的是三相总功率有的给的是单相功率做标幺化之前一定要先确认不然算出来的电压降落和损耗会差好几倍。4.2 Yalmip建模与求解这个项目里优化模型的变量类型有连续变量各机组出力、购电量、储能功率和0-1整数变量燃气轮机的启停状态、储能的充放状态本质是一个混合整数二阶锥规划MISOCP问题。用Yalmip建模非常合适它能让你把精力集中在模型本身而不是求解器的接口语法上。核心建模代码如下%% 定义优化变量 % 时间维度 T 24; % 24小时 % 决策变量定义 P_grid sdpvar(1, T); % 从上级电网购电功率单位MW P_gt sdpvar(2, T); % 两台燃气轮机出力 u_gt binvar(2, T); % 燃气轮机启停状态 P_ch sdpvar(1, T); % 储能充电功率 P_dis sdpvar(1, T); % 储能放电功率 SOC sdpvar(1, T1); % 储能荷电状态 P_wind sdpvar(1, T); % 风电实际消纳功率 P_pv sdpvar(1, T); % 光伏实际消纳功率 %% 构建目标函数 % 购电成本分时电价 price_buy [0.48*ones(1,7), 0.38*ones(1,8), 1.2*ones(1,5), 0.48*ones(1,4)]; cost_buy sum(price_buy .* P_grid); % 燃气轮机燃料成本二次函数碳价格直接计入 carbon_price 50; % 元/吨 gt_emission_factor 0.2; % t/MWh gt_fuel_coef [0.02 0.03] * carbon_price; % 简化成本系数 cost_gt sum(gt_fuel_coef * P_gt_quad); % 这里P_gt_quad需要额外定义二次项 % 储能运维成本 cost_bess 0.01 * sum(P_ch P_dis); % 单位运维成本0.01元/kWh % 弃风弃光惩罚 penalty_curtail 0.8 * sum( (P_wind_max - P_wind) (P_pv_max - P_pv) ); % 碳成本外购电燃气轮机排放量乘以碳价 emission_grid 0.5 * sum(P_grid); % 电网排碳因子0.5 t/MWh emission_gt 0.2 * sum(P_gt(1,:) P_gt(2,:)); % 燃气轮机排碳因子0.2 t/MWh cost_carbon carbon_price * (emission_grid emission_gt); Objective cost_buy cost_gt cost_bess penalty_curtail cost_carbon;有一点我要特别提醒燃气轮机的二次成本函数直接放进目标函数会让模型变成二次规划和SOCP约束叠加后求解会明显变慢。我自己的做法是直接把二次成本分段线性化或者用常数边际成本加固定启停成本近似这样保持模型线性和MILP的结构求解稳定性和速度都好很多。精确二次成本和线性近似的差异在这个项目中通常小于3%但求解时间可能差一个数量级。4.3 约束条件代码实现约束条件的代码比目标函数更需要细心。我按类别把约束写清楚%% 约束集合 Constraints []; %% 节点功率平衡约束DistFlow线性化版本 % 对每个节点注入功率 该节点净负荷 流出功率之和 % 这里简化为系统级功率平衡实际IEEE33节点需要逐节点展开 % 系统总功率平衡 Constraints [Constraints, ... P_grid(1,:) sum(P_gt,1) P_dis - P_ch P_wind P_pv ... P_load_total losses_estimate]; %% 燃气轮机约束 % 出力上下限 P_gt_max [0.5, 0.5]; % 单台最大出力0.5MW P_gt_min [0.02, 0.02]; % 最小出力 for t 1:T Constraints [Constraints, ... P_gt_min(1) * u_gt(1,t) P_gt(1,t) P_gt_max(1) * u_gt(1,t)]; Constraints [Constraints, ... P_gt_min(2) * u_gt(2,t) P_gt(2,t) P_gt_max(2) * u_gt(2,t)]; end %% 储能约束 % SOC递推 eta_ch 0.95; % 充电效率 eta_dis 0.95; % 放电效率 SOC_min 0.1; SOC_max 0.9; for t 1:T Constraints [Constraints, ... SOC(t1) SOC(t) eta_ch * P_ch(t) - P_dis(t) / eta_dis]; Constraints [Constraints, ... SOC_min SOC(t1) SOC_max]; Constraints [Constraints, ... 0 P_ch(t) 0.3]; Constraints [Constraints, ... 0 P_dis(t) 0.3]; end % 首末SOC一致 Constraints [Constraints, SOC(1) SOC(T1)]; %% 可再生能源消纳约束 Constraints [Constraints, 0 P_wind(1,:) P_wind_max]; Constraints [Constraints, 0 P_pv(1,:) P_pv_max];实际跑这个模型的时候我会把节点功率平衡约束用矩阵形式一次性写出来而不是用循环逐时写入这样Yalmip内部构建稀疏矩阵的效率要高一些。24时段的调度模型可能感觉不到差别但一旦要做数百组灵敏度算例这个优化能明显缩短总耗时。4.4 灵敏度扫描的循环结构灵敏度分析不是模型的核心但却是整个项目的“产出放大器”。我的代码逻辑是这样%% 灵敏度分析碳价扫描 carbon_price_list 0:10:300; results zeros(length(carbon_price_list), 5); for i 1:length(carbon_price_list) % 更新目标函数中的碳价参数 CP carbon_price_list(i); % 重新建目标函数这里用函数封装更好避免变量覆盖 % 求解 ops sdpsettings(solver, gurobi, verbose, 0); optimize(Constraints, Objective, ops); % 记录结果 results(i,1) CP; results(i,2) value(cost_total); % 总成本 results(i,3) value(emission_total); % 总排放 results(i,4) value(sum(P_gt, all)); % 燃气轮机总发电量 results(i,5) value(sum(P_grid)); % 电网购电量 end %% 画图 figure; yyaxis left; plot(results(:,1), results(:,2), -o); ylabel(总成本 (元)); yyaxis right; plot(results(:,1), results(:,3), -s); ylabel(碳排放 (t)); xlabel(碳价 (元/t)); legend(总成本, 碳排放);这里有个实操细节不要在循环里用同样的变量名反复重建模型因为Yalmip会把旧的变量保留在内存里变量数越来越多后面重建模型会越来越慢。我的做法是写一个函数输入碳价参数返回优化结果每次调用一个全新的函数工作区问题会小很多function result run_single_case(carbon_price, load_scale) % 完整的建模、约束构建、求解代码 end然后主脚本里只做循环调用和数据汇总。实测下来一个31组参数扫描的算例总耗时从“循环内建模”的15分钟降到了“函数调用”的2分半钟差距非常惊人。5. 结果解读与工程经验5.1 典型灵敏度曲线说明了什么我在自己的算例里得到的结果大致有这样的规律碳价从0升到60元/吨这个区间调度结果几乎不变化燃气轮机不启机系统完全依赖电网购电。原因很简单当时电价低谷时段0.38元/kWh燃气轮机的边际发电成本即使不算碳排放都高于这个价格纯经济比较没有竞争力。碳价上升到80到120元/吨区间时燃气轮机开始启动主要替代的是晚高峰时的外购电。因为晚峰电价高到1.2元/kWh燃气轮机在峰值时段发电的综合成本已经低于外购电了。这个阶段碳排放总量开始下降但购电量的减少并没有等比例降低碳排放因为燃气轮机自己也有排碳。碳价超过150元/吨以后系统的行为进入第三阶段储能开始显著增加峰谷套利力度可再生能源的弃电率降到接近零。此时系统的碳排放下降趋势变缓说明单纯靠提高碳价推动减碳的效果开始递减。这个三段式响应是符合工程直觉的。它说明经济-碳协调不是线性的“碳价越高排放越低”而是存在明显的阈值效应。这些阈值对应着不同电源之间经济性的边际切换点在做碳排放政策制定或企业成本评估时非常有参考价值。5.2 常见问题与排查实录我复现这套代码时踩过的坑按频率排序如下第一Yalmip变量名冲突导致约束混乱。试着在一个脚本里多次运行不同碳价如果用lazy清理可能上一轮的部分约束会被保留下来导致新模型无解或结果异常。解决方案就是函数化封装每次新工作区。第二储能SOC约束出现“跳绳”现象。所谓跳绳就是优化结果里SOC曲线在相邻时段之间剧烈跳变——这一小时0.9下一小时直接变成0.1。原因是储能约束里少了充放电互斥约束。解决办法是引入充放电0-1状态变量确保同一时刻只能有一个功率为正u_bin binvar(1, T); % 1为充电0为放电 Constraints [Constraints, 0 P_ch 0.3 * u_bin]; Constraints [Constraints, 0 P_dis 0.3 * (1 - u_bin)];第三标幺化单位搞混导致约束数量级失衡。比如储能功率是0.3MW购电功率可以到5MW如果单位不统一某些约束的宽容度和惩罚系数设定就非常别扭。我建议在优化模型内部全部使用标幺值MW、MWh展示结果时再转换回有名值。5.3 如何扩展和升级这个项目做完之后我给它预留了几个扩展方向你完全可以在此基础上继续做考虑多场景随机优化用场景生成法处理风电光伏出力不确定性取代确定性预测值目标函数用期望成本。加入需求响应把一部分可转移负荷建模为柔性负荷让优化自己决定什么时候供电这会对系统经济性和碳排放产生更显著的影响。引入多时段耦合的碳捕集设备碳捕集设备会同时消耗电力和产生碳减排收益与碳价形成更强的耦合关系模型会更复杂但更有现实意义。网络重构与调度联合优化把IEEE33节点的联络开关状态也作为决策变量实现拓扑与运行的联合优化。我之前尝试过把动态重构加入模型结果是混合整数问题规模翻了好几倍24小时算例从几秒增加到了几分钟。如果把重构和碳价灵敏度放在一起做建议你先固定拓扑跑灵敏度再在某个关键碳价点附近做小范围的重构-调度联合寻优。5.4 最后分享一个小技巧用碳价做灵敏度分析时有一个容易被忽略的细节免费碳排放配额怎么设置。我在第一版代码里把配额设成了总负荷的固定比例结果碳价升高时系统仍然大量排放因为配额已经覆盖了大部分排放边际碳价根本不起作用。后来我把配额设成与具体电源结构挂钩——配额只覆盖上级电网购电排放的80%燃气轮机的排放完全不覆盖这样燃气轮机的启动真正需要承担碳排放成本的分担灵敏度曲线的拐点就变得非常清晰。这个经验放到实际工程项目里也是一个值得注意的地方碳配额的归属方式会直接改变优化结果的行为模式写论文时务必在假设里说明清楚免得审稿人质疑你的灵敏度结论是在“作弊”的基础上得到的。