做了好几年电力系统调度相关的代码我最大的体会是——安全约束这东西平时看不见可一旦真出事它是唯一能兜底的关口。今天想聊的这套“计及N-k安全约束的含光热电站电力系统优化调度模型”简单说就是在常规经济调度的基础上强制要求调度计划在任意k个关键元件线路、变压器、机组同时停运时仍然可用同时系统里还接入了光热电站这种带储热特性的可调度可再生能源。它适合两类人看一类是正在复现相关论文、准备写毕业设计或期刊代码的研究生另一类是想把N-k约束补进自己调度程序里的工程师。这篇文章不会只给结论我会把整条链路拆开讲N-k约束到底在约束什么、光热电站的储热方程怎么进优化器、IEEE14节点和IEEE118节点算例该怎么设计最后给出一套可以在Matlab里跑起来的代码骨架以及我调试时踩过的几个坑。读完你至少能照着把模型搭起来并且知道求解器报错时该先去查哪里。1. 为什么要把N-k安全约束写进调度模型从“可用性”到“生存性”的约束升级1.1 N-1/N-k到底在查什么电力系统里说的“N”指的是系统关键元件的总数比如线路、变压器、发电机组都算。传统的经济调度只考虑正常情况下所有元件都在线的状态目标是最小化发电成本输出一组机组出力和潮流分布。但电网运行里最麻烦的恰恰是“正常状态”随时可能被打破一条线路遭雷击跳闸、一台机组因故障解列这些都是分钟级甚至秒级的事。N-1安全准则的意思是系统失去任意一个元件之后仍然不能出现负荷丢失、线路过载或电压越限。这是绝大多数实际电网调度运行的最低安全门槛。而N-k是它的推广形式k可以取2、3甚至更高代表的是“任意k个元件同时停运”这种更极端但仍然可能发生的场景比如同塔双回线同时遭外力破坏、一座变电站内多条出线同时跳闸。我要提醒一个容易混淆的点N-k不是“所有N个元件里必须坏k个”而是“在任何一组k个元件停运的组合下系统都必须保持安全”。所以当k增大时需要考虑的场景组合数量是组合数C(N, k)这个数字涨得非常快。以一个20条线路的系统为例N-1只有20个场景N-2就有190个N-3直接到1140个。这也是为什么很多调度模型在论文里写着“考虑N-k”实际代码里却只跑了N-1因为场景规模会迅速压垮求解器。1.2 为什么用了光热电站之后N-k约束反而更难处理光热电站和光伏最大的区别在于它自带储热系统白天可以把多余的热量存在熔盐罐里晚上再释放出来发电。听起来很灵活但这恰恰让N-k约束变得更复杂。传统火电或者水电的N-k校验主要关心机组停运后剩下机组能否顶上、线路潮流是否越限。光热电站接入后它既是电源又是储能运行状态受一天之内太阳辐照度变化的影响还要满足储热罐的能量平衡。这意味着在做N-k场景校验时不能简单地假设光热电站“满发”而要看它在当前调度时刻是否有足够的热能支撑出力承诺。换句话说光热电站的N-k时段可用容量是动态变化的它取决于此前若干时段的集热和储热状态。我在实现时遇到的一个典型现象是如果只把光热电站当作一个普通可调电源设定它的最大出力为装机容量那么调度结果会显得很乐观N-k场景校验时光热也能顶上去。但一旦把储热罐容量和能量平衡写进约束就会发现在连续阴天或傍晚时段光热的实际可发电力远低于装机值原本认为安全的调度计划在N-k校验中会暴露出大量的负荷丢失。这就是“物理模型太粗糙导致安全评估失真”的典型例子。1.3 预防性安全与纠正性安全两种约束写法差距很大在把N-k约束写进优化模型时先要想清楚写的是“预防性安全”还是“纠正性安全”。预防性安全的写法是当前正常状态下的调度计划直接放到任意故障场景下都必须满足约束故障发生后不需要再调整任何机组出力。这种写法实现最简单把每个故障场景的潮流约束都加到同一个优化问题里就行但代价是正常运行成本会比较高因为某些机组必须预留更多出力空间或保持更高开机容量来应对极端场景。纠正性安全则允许故障发生后短时间内对可用机组进行重新调度比如把某些机组出力提高、把光热电站储热释放出来。典型做法是给每个故障场景单独定义一套机组出力调整变量并限制调整幅度在某个范围比如±20%额定出力、调整时间10分钟。这么写更贴近实际调度模型也更经济但代码复杂度和求解难度都会上升一个台阶。我在自己的实现里采用了一个折中方案正常调度计划使用预防性安全约束保证基态结果能在所有预想故障场景下维持不失负荷但对每台机组添加一个“短时过载能力”参数在故障场景校验时允许机组出力在短时间内越过正常运行上限。这样既不会把模型搞成完全每个场景一套变量的庞大结构又能反映实际机组具备的短时调节能力。2. 光热电站不是“带电池的光伏”CSP的调度模型该怎么建2.1 光热电站的物理链条与三个子系统光热电站Concentrating Solar Power简称CSP在整个调度模型里需要拆成三个部分来看集热场、储热罐、汽轮发电机组。太阳直射辐射DNI照射到集热器上被转化为热能热能一部分直接送入蒸汽发生器发电另一部分存入熔盐罐需要发电时再从储热罐取出热量驱动汽轮机。在数学上这套能量链条可以写成三层关系。第一层是集热场的热功率上限Q_sf(t) η_sf × DNI(t) × A_sf其中A_sf是集热场有效采光面积η_sf是光热转换效率DNI是当前时段的直射辐照度。这一层决定了“太阳能进来多少”。第二层是储热罐的能量平衡E_tes(t) E_tes(t-1) (Q_sf(t) - Q_pb(t) - Q_spill(t)) × Δt这里Q_pb是送入汽轮机发电的热功率Q_spill是集热场超出需求后被舍弃的热功率用大白话说就是“收进来的热发掉一部分、存一部分、实在多到装不下就扔一部分”。第三层是热电转换P_csp(t) η_pb × Q_pb(t)η_pb是汽轮发电机组的热电转换效率。整条链路的物理意义可以类比成一个移动电源太阳能是充电器熔盐罐是电池汽轮机是用电器。充电速度受阳光限制电池容量有限放电速度还要受汽轮机自身特性的约束。调度模型要做的就是在这个链条的每一个环节上都写上对应约束。2.2 储热罐的时序耦合是最容易丢的约束我见过很多第一次建CSP模型的人把光热电站简单处理成三个不等式出力有上下限、有爬坡约束、有电量限制然后就没有然后了。这么做最直接的后果是储热罐的跨时段能量状态丢失了优化器可以“凭空”在某个时段给出很高的出力而根本不管之前有没有存够热。实际上储热罐是一个强时序耦合的状态量第t时段的储热水平由第t-1时段的状态和本时段的热输入输出共同决定。这正是CSP调度里最核心、也最容易被忽略的部分。在设计模型时一定要把E_tes作为决策变量放进优化问题并显式写出相邻时段的递推关系。同时还要加两条硬性约束第一条是储热罐容量上下界0 ≤ E_tes(t) ≤ C_tes其中C_tes由熔盐罐的物理容量决定。第二条是调度周期末段的储热水平约束通常写成E_tes(T) ≥ λ × C_tesλ取0.1到0.3之间。我最早跑24小时调度时没有加这条末期约束结果优化器把储热罐在夜里全部放空第二天早上系统又面临负荷高峰光热电站无法提供支撑出力曲线看起来就像“为了省弃热而透支了明天的能力”。加上末期约束之后结果立刻合理很多。2.3 CSP机组的运行区间与爬坡限制不是想发就发光热电站的汽轮机组本质上是一台热力机组它有最小技术出力有爬坡速率限制也有启停状态约束。很多简化模型把P_csp当连续变量用0到额定容量直接约束这在大多数情况下是过度乐观的。实际汽轮机在低负荷区域运行效率极差甚至可能无法稳定运行。我通常会给光热电站定义一个正常运行区间比如20%到100%额定出力低于20%时必须停机。同时爬坡速率也要写进去熔盐储热虽然能平抑太阳波动但汽轮机自身的进气阀门、缸体热应力决定了它不能像锂电池放电那样瞬间改变出力。实测下来一台50MW的CSP机组爬坡速率一般取每分钟2%到5%额定出力折合到小时级调度模型里就是每时段24到60MW/h的变化上限。如果你用的是小时级调度模型爬坡约束可能没有那么突出但启停状态一定要考虑。我的做法是用一个0/1变量u_csp(t)表示汽轮机是否在线然后把P_csp(t)写成u_csp(t) × P_csp_min ≤ P_csp(t) ≤ u_csp(t) × P_csp_max同时加上最小开机时间和最小停机时间约束否则优化器会为了某个时段的负荷波动而让汽轮机频繁启停这在物理上根本实现不了。这一步做完CSP才真正像一个可以参与调度调节的电源而不是一个“随叫随到的电热水壶”。3. 模型总装目标函数、运行约束与N-k故障场景的表达式3.1 目标函数怎么设置才能兼顾经济性和N-k惩罚我用的目标函数是系统总运行成本最小化包含三个层次min ∑ₜ [ ∑₉ (a_g × P_g(t)² b_g × P_g(t) c_g × u_g(t)) c_csp × P_csp(t) ] ∑ₛ π_s × ∑ₜ c_curtail × P_curtail(s,t)第一项是常规火电机组的发电成本通常用二次函数拟合煤耗第二项是光热电站的运行维护成本取值很小比如2到5美元/MWh主要让优化器不会在两种可行方案间随意抖动第三项是故障场景下的失负荷惩罚c_curtail要取得非常大远高于正常发电成本保证优化器只有在极端场景下才会选择切负荷。这里有一个值得注意的设计细节N-k安全约束本身并不直接进目标函数而是通过约束条件来体现。目标函数里的失负荷惩罚是为了在“故障场景确实无法完全满足负荷需求”时给优化器一条软性退路。比如一个故障场景把系统分裂成孤岛孤岛内发电容量不足这时硬性要求功率平衡且不失负荷优化问题就会直接无解。加入失负荷变量并配上高惩罚系数后问题永远有解同时优化器会尽量避免失负荷。3.2 基础运行约束的完整集合在写N-k场景约束之前先把正常状态下的基础约束列齐。我在Matlab实现中通常包含以下五组一是节点有功平衡约束对每个节点i和时段tΣ_g∈Gi P_g(t) P_csp(t) - P_load_i(t) Σ_l B_il × (θ_i(t) - θ_l(t))这里用的是直流潮流模型忽略无功和网损把节点注入有功和支路有功潮流用节点相角表示。二是机组出力上下界约束P_g_min × u_g(t) ≤ P_g(t) ≤ P_g_max × u_g(t)同时要有爬坡约束P_g(t) - P_g(t-1) ≤ RU_gP_g(t-1) - P_g(t) ≤ RD_g。三是机组最小启停时间约束这个需要一组二进制变量表示启停事件并用专门的不等式组实现。很多入门代码会漏掉这一条导致优化结果里同一台机组在相邻时段反复启停。四是线路潮流限值约束对每条支路l|B_l × (θ_from(l) - θ_to(l))| ≤ P_l_max。五是CSP子系统约束包括上一节写的集热场热功率约束、储热罐能量平衡、容量限额、末期储热约束以及汽轮机的出力区间和启停约束。3.3 N-k故障场景的建模预想故障集、孤岛处理与失负荷变量N-k约束的加入方式核心是对每个预想故障场景s重复一遍潮流安全约束只是系统拓扑变成了“移除故障元件之后的拓扑”。具体做法是预先枚举出预想故障集F对每个场景s∈F生成对应的节点电纳矩阵B_s。如果场景s涉及线路断线就把该线路从支路列表中剔除重新组装B矩阵如果涉及机组故障就限制该机组出力强制为0。然后对每个场景s添加以下约束Σ_g∈G_avail(s) P_g(t) P_csp(t) P_curtail(s,t) P_load(t)以及该场景下的支路潮流限额约束-B_l_max ≤ B_s × (θ_i(s,t) - θ_j(s,t)) ≤ B_l_max这里每个故障场景其实都有一套自己的相角变量θ(s,t)不能直接沿用正常状态的θ因为系统拓扑变了之后相角分布必须重新求解。我开始实现时图省事想用同一套θ变量同时校验所有场景结果发现一条线路断线后的相角物理上就和正常状态完全不同这样约束本身就没有意义了。正确的做法是给每个场景单独定义一组相角变量把拓扑变化的影响真正放进优化问题。失负荷变量P_curtail(s,t)是保证可行性的关键。故障场景s中如果可用发电容量确实小于负荷优化器会切掉一部分负荷但目标函数里对应的高惩罚会让它尽量不这么做。实际结果解读时可以统计每个场景的最大失负荷量以此判断系统当前调度计划的N-k安全裕度。关于预想故障集怎么选我的建议是分三层第一层是全部N-1场景数量少且物理意义明确第二层是高风险N-2场景比如同塔双回线、同母线多出线、共用通道的线路组第三层是根据基态潮流结果筛选的高负载率线路组合。不要一上来就把所有C(N,k)组合全塞进模型那样求解器会先把你自己的耐心耗光。4. 算例设计IEEE14节点做机理验证IEEE118节点做规模压力测试4.1 两套标准算例系统的基本差异算例选型上我同时跑IEEE14节点和IEEE118节点两套系统目的完全不同。IEEE14节点系统规模小14条母线、约5台机组、20条线路左右非常适合做机制验证。在这个系统上N-1场景全枚举只有20个左右N-2线路场景也不过190个全部塞进优化模型也能在几十秒内求解。我习惯用它来检查约束逻辑是否正确比如光热电站的热能平衡是否合理、故障场景下线路潮流是否真的没有越限、失负荷惩罚是否真的只在极端场景才触发。IEEE118节点系统规模大得多118条母线、54台机组左右、线路约186条是检验算法效率和可扩展性的“压力测试机”。在这个系统上如果仍然全枚举所有N-2场景就是C(186,2)约1.7万个场景直接建一个包含所有场景约束的MILP模型规模会膨胀到求解器难以接受的程度。我在118节点上只保留了全部N-1场景和约50个高风险N-2场景才把求解时间控制在可接受的分钟级范围。下表是我实现时关注的关键差异指标 | IEEE14节点 | IEEE118节点 母线数 | 14 | 118 机组数 | 约5 | 约54 支路数 | 约20 | 约186 N-1场景规模 | 几十个 | 两百多个 N-2全枚举难度 | 容易可直接枚举 | 组合爆炸必须筛选 主要用途 | 逻辑验证、教学演示 | 大规模算法评测、性能测试需要注意IEEE标准算例的不同版本在机组数、线路参数上会有差异我从Matpower的case14和case118数据文件读入数据时就发现机组报价参数很多版本里根本没提供需要自己根据煤耗特性补一套。算例数据的版本一致性会在结果对比时非常折腾建议在代码开头固定数据来源版本并写进注释。4.2 N-2场景选择策略与故障集设计IEEE118节点上N-2全枚举不现实我的筛选策略是按“故障后果严重度”来取舍。先用正常状态的最优调度结果做一次潮流计算列出负载率最高的前20条线路然后在这些线路之间以及它们周围电气距离较近的线路之间生成N-2组合。同时手动添加经典的高风险场景同一走廊的双回线、同一厂站出线、主变与其高压侧线路组合。这里有个容易被忽略的点N-k约束不是场景越多越好。把不相关的低风险场景加进模型不仅无助于提高安全性还会白白增加求解时间。我在118节点上做过对比把N-2场景从全部1.7万个缩减到前50个高负载场景后求解时间从数小时量级降到几分钟量级而最严重场景的失负荷指标几乎没有变化因为低风险场景本来就不会造成约束起作用。4.3 什么结果才算“调度模型真正起到了作用”算例跑完之后一定要看几个关键指标来判断模型有没有真正发挥作用。第一个指标是带N-k约束和不带N-k约束两种情况下的正常运行成本差。如果两种情况下成本几乎一样说明预想故障集选得太温和安全约束没有真正绷紧。正常情况下考虑N-k安全约束后运行成本会上升幅度一般在2%到10%这部分成本就是“买安全保险”的代价。第二个指标是故障场景下最大失负荷量。如果没有N-k约束某些高风险场景下可能直接出现数十MW的失负荷加上约束之后这个量应该显著下降或降为零。我在IEEE14节点上验证时最严重的N-2故障场景从失负荷约30MW降到0MW代价是正常状态总成本上升约4%这个权衡关系非常清楚。第三个指标是求解时间和MIP gap。IEEE118节点上如果gap设置过松比如1%优化器会提前终止此时看似有一种可行解但可能距离真正的最优解很远而且安全约束虽然满足运行成本却虚高不少。我的经验是把gap控制在0.1%以内求解时间允许放宽到10分钟。5. MatlabYALMIP代码结构从数据读入到Gurobi求解的完整骨架5.1 主程序的七步骨架与数据组织我整个代码按照下面这个流程组织每个子函数各自独立方便替换算例数据和故障集。以Matlab配合YALMIP工具箱为框架求解器用Gurobi或CPLEX。七步骨架是读数据、定义变量、收集正常态约束、收集N-k场景约束、设置目标函数、求解、后处理校验。主函数开头长这样function [result] csp_nk_schedule(casefile, K) mpc loadcase(casefile); % 读取IEEE标准算例 T 24; % 调度时段数 dt 1; % 时段时长小时 faultSet build_fault_set(mpc, K); % 构造预想故障集 %% 定义决策变量YALMIP Pg sdpvar(mpc.ng, T, full); % 常规机组出力 ug binvar(mpc.ng, T, full); % 机组启停状态 theta sdpvar(mpc.nb, T, full); % 正常态节点相角 theta_s cell(length(faultSet), 1); % 各故障场景相角 Pcsp sdpvar(1, T, full); % 光热电站出力 Etes sdpvar(1, T, full); % 储热罐能量 Pcurtail sdpvar(length(faultSet), T, full); % 场景失负荷 %% 收集约束 C []; C [C, BaseConstraints(...)]; for s 1:length(faultSet) theta_s{s} sdpvar(mpc.nb, T, full); C [C, FaultScenarioConstraints(...)]; end %% 目标函数与求解 objective ObjectiveFunc(Pg, ug, Pcsp, Pcurtail); ops sdpsettings(solver, gurobi, mipgap, 1e-4, verbose, 1); optimize(C, objective, ops); %% 回读并校验结果 result extract_result(...); end这段代码里有一个关键点每个故障场景都单独定义了theta_s{s}这一套相角变量绝不能复用正常状态的theta。我第一次实现时为了减少变量数偷懒复用theta结果约束虽然形式上写了实际却完全没有反映故障拓扑算出来的结果在N-k校验里一塌糊涂。5.2 YALMIP变量声明、约束收集与求解器配置YALMIP里最坑的一个细节是矩阵变量的默认对称性。sdpvar(n, T)默认声明一个n×T的对称矩阵这对调度问题来说是完全错误的。机组出力矩阵需要每一行是一台机组、每一列是一个时段所以必须写成sdpvar(ng, T, full)加上第三个参数full告诉YALMIP这是完整矩阵。我一开始就是漏了这个full导致模型维度莫名其妙扩大求解速度极慢还报出一堆“变量越界”的错误。约束收集时另一个值得注意的点是不要在循环里反复用方括号拼接超大的约束数组。场景数多、变量维数大时YALMIP内部每次拼接都会复制整个约束集内存开销非常离谱。我的做法是为每个故障场景先生成一个单独的约束元胞数组subC{s}全部循环结束后再用方括号一次性展开C_all [C_base, subC{:}];求解器配置上我比较常用的设置是ops sdpsettings(solver, gurobi, mipgap, 1e-4, ... mipfocus, 1, timelimit, 600, ... numericalemphasis, 1, verbose, 2);mipgap设成1e-4代表相对最优间隙到0.01%才停保证结果可信。numericalemphasis设为1让求解器在遇到数值病态问题时更积极处理对带多种约束尺度差异的调度问题很有帮助。5.3 一段可以直接参考的CSP储热约束代码下面这段代码是我在Matlab中实现光热电站核心约束的片段可以直接嵌到模型里。% 输入已知量 DNI load_dni_curve(); % 24小时DNI曲线W/m2 A_sf 2.5e6; % 集热场面积m2 eta_sf 0.4; % 光热转换效率 eta_pb 0.38; % 热电转换效率 C_tes 300; % 储热罐容量MWh P_csp_max 50; % 光热额定出力MW P_csp_min 10; % 光热最小技术出力MW lambda 0.2; % 末期储热比例 % 决策变量 % Q_sf_total 是集热场热功率Q_spill 是弃热 % Q_pb 是送入汽轮机的热功率 % 集热场热功率上限 C [C, Q_sf_total eta_sf * A_sf * DNI(t) / 1e6]; % 单位MW % 储热罐能量平衡 C [C, Etes(t) Etes(t-1) dt * (Q_sf_total - Q_pb - Q_spill)]; % 储热罐容量与末期约束 C [C, 0 Etes C_tes]; C [C, Etes(T) lambda * C_tes]; % 电热转换与机组运行区间 C [C, Pcsp(t) eta_pb * Q_pb(t)]; C [C, u_csp(t) * P_csp_min Pcsp(t) u_csp(t) * P_csp_max];如果你在实际项目中看到“光热出力曲线很漂亮但储热罐能量毫无变化”的结果那大概率就是储能平衡约束没写进模型或者在YALMIP变量声明时把Etes误声明成了对称矩阵。6. 求解性能优化与我最常遇到的几个报错陷阱6.1 组合爆炸场景筛选比硬解全部组合靠谱N-k约束的最大敌人永远是组合爆炸。我自己的一个深刻教训是在IEEE118节点上第一次尝试全枚举N-2场景代码写完后一跑发现YALMIP光是模型编译阶段就吃掉了近20GB内存Gurobi在半小时内一条可行解都没出来。后来改用“基态潮流预筛选”的方法只挑出线路负载率排行前20的线路及其组合模型规模瞬间缩小了两个数量级。场景筛选的落地做法是先用不带N-k约束的基本调度模型求一版结果然后对这个结果做全场景N-k潮流扫描把实际出现线路越限或母线电压问题的故障场景提取出来作为约束集加入优化模型再重新求解迭代2到3轮。这种“预调度—故障筛选—再调度”的外循环方式在工程上远比一次性枚举所有场景高效而且结果基本等价。不少商业调度软件内部也是这么做的。6.2 直流潮流导致的孤岛失负荷问题DC潮流模型把系统当作纯有功网络线路停运后如果系统被分成两个电气孤岛每个孤岛内部必须自己保持功率平衡。有些预想故障场景会产生完全没有电源的孤岛负荷比如某条线断开后一小片负荷区域与主网失去连接而该区域又没有本地机组。这种“孤儿负荷”场景放进模型后如果没有失负荷变量优化问题直接无解。我的处理方式是给所有负荷节点都加上故障场景失负荷变量P_curtail(s,t)并允许其在极端场景下切掉部分或全部负荷。目标函数里的失负荷惩罚取一个极大值正常场景下优化器不可能主动切负荷真遇到孤岛无电源的场景时切负荷变量兜住可行性模型仍然有解同时结果报告里会把这个孤岛失负荷量明确标出来。这一步做完N-k调度在DC潮流框架下的数值稳定性会好很多。6.3 结果不合理时的排查顺序跑出来的结果不对劲时我习惯按下面这个顺序排查效率高很多。如果你手头的模型也出现了类似问题可以直接照这个思路走症状 | 可能原因 | 排查方向 求解时间爆炸 | 故障场景全枚举过多 | 改用迭代场景筛选限制N-2组合范围 正常态和故障态结果几乎一致 | 故障场景约束没有真正生效 | 检查场景s的相角变量是否独立B_s是否重新组装 光热电站全天满发 | 储热末期约束缺失 | 加末期最低储热水平约束 机组频繁启停 | 最小启停时间缺失或gap过松 | 补启停时间约束收紧mipgap 无解且提示infeasible | 孤岛负荷无电源且无失负荷变量 | 给故障场景加失负荷变量和高惩罚系数 线路潮流普遍偏大但安全校验通过 | B矩阵组装错误 | 对比正常态B_s与去掉线路后的B_s差异最后再分享一个我自己调试时的高频操作任何时候怀疑N-k约束没生效最简单粗暴的验证方法是人为设置一个必会造成越限的故障场景比如把系统中负载最重的线路单独断开。如果约束正确模型要么自动调整机组出力消除越限要么在结果中给出该场景失负荷量如果两种迹象都没有说明某个环节的约束根本没有写进优化器。这一招在IEEE14节点上验证逻辑特别方便几乎一测一个准。从14节点搬到118节点的过程中我最大的感慨是让N-k安全约束真正落地难的从来不是约束公式本身而是场景管理和求解器的脾气。希望这篇笔记能帮你少踩几个我已经替你踩过的坑。
计及N-k安全约束的光热电站优化调度：模型、算例与Matlab实现
做了好几年电力系统调度相关的代码我最大的体会是——安全约束这东西平时看不见可一旦真出事它是唯一能兜底的关口。今天想聊的这套“计及N-k安全约束的含光热电站电力系统优化调度模型”简单说就是在常规经济调度的基础上强制要求调度计划在任意k个关键元件线路、变压器、机组同时停运时仍然可用同时系统里还接入了光热电站这种带储热特性的可调度可再生能源。它适合两类人看一类是正在复现相关论文、准备写毕业设计或期刊代码的研究生另一类是想把N-k约束补进自己调度程序里的工程师。这篇文章不会只给结论我会把整条链路拆开讲N-k约束到底在约束什么、光热电站的储热方程怎么进优化器、IEEE14节点和IEEE118节点算例该怎么设计最后给出一套可以在Matlab里跑起来的代码骨架以及我调试时踩过的几个坑。读完你至少能照着把模型搭起来并且知道求解器报错时该先去查哪里。1. 为什么要把N-k安全约束写进调度模型从“可用性”到“生存性”的约束升级1.1 N-1/N-k到底在查什么电力系统里说的“N”指的是系统关键元件的总数比如线路、变压器、发电机组都算。传统的经济调度只考虑正常情况下所有元件都在线的状态目标是最小化发电成本输出一组机组出力和潮流分布。但电网运行里最麻烦的恰恰是“正常状态”随时可能被打破一条线路遭雷击跳闸、一台机组因故障解列这些都是分钟级甚至秒级的事。N-1安全准则的意思是系统失去任意一个元件之后仍然不能出现负荷丢失、线路过载或电压越限。这是绝大多数实际电网调度运行的最低安全门槛。而N-k是它的推广形式k可以取2、3甚至更高代表的是“任意k个元件同时停运”这种更极端但仍然可能发生的场景比如同塔双回线同时遭外力破坏、一座变电站内多条出线同时跳闸。我要提醒一个容易混淆的点N-k不是“所有N个元件里必须坏k个”而是“在任何一组k个元件停运的组合下系统都必须保持安全”。所以当k增大时需要考虑的场景组合数量是组合数C(N, k)这个数字涨得非常快。以一个20条线路的系统为例N-1只有20个场景N-2就有190个N-3直接到1140个。这也是为什么很多调度模型在论文里写着“考虑N-k”实际代码里却只跑了N-1因为场景规模会迅速压垮求解器。1.2 为什么用了光热电站之后N-k约束反而更难处理光热电站和光伏最大的区别在于它自带储热系统白天可以把多余的热量存在熔盐罐里晚上再释放出来发电。听起来很灵活但这恰恰让N-k约束变得更复杂。传统火电或者水电的N-k校验主要关心机组停运后剩下机组能否顶上、线路潮流是否越限。光热电站接入后它既是电源又是储能运行状态受一天之内太阳辐照度变化的影响还要满足储热罐的能量平衡。这意味着在做N-k场景校验时不能简单地假设光热电站“满发”而要看它在当前调度时刻是否有足够的热能支撑出力承诺。换句话说光热电站的N-k时段可用容量是动态变化的它取决于此前若干时段的集热和储热状态。我在实现时遇到的一个典型现象是如果只把光热电站当作一个普通可调电源设定它的最大出力为装机容量那么调度结果会显得很乐观N-k场景校验时光热也能顶上去。但一旦把储热罐容量和能量平衡写进约束就会发现在连续阴天或傍晚时段光热的实际可发电力远低于装机值原本认为安全的调度计划在N-k校验中会暴露出大量的负荷丢失。这就是“物理模型太粗糙导致安全评估失真”的典型例子。1.3 预防性安全与纠正性安全两种约束写法差距很大在把N-k约束写进优化模型时先要想清楚写的是“预防性安全”还是“纠正性安全”。预防性安全的写法是当前正常状态下的调度计划直接放到任意故障场景下都必须满足约束故障发生后不需要再调整任何机组出力。这种写法实现最简单把每个故障场景的潮流约束都加到同一个优化问题里就行但代价是正常运行成本会比较高因为某些机组必须预留更多出力空间或保持更高开机容量来应对极端场景。纠正性安全则允许故障发生后短时间内对可用机组进行重新调度比如把某些机组出力提高、把光热电站储热释放出来。典型做法是给每个故障场景单独定义一套机组出力调整变量并限制调整幅度在某个范围比如±20%额定出力、调整时间10分钟。这么写更贴近实际调度模型也更经济但代码复杂度和求解难度都会上升一个台阶。我在自己的实现里采用了一个折中方案正常调度计划使用预防性安全约束保证基态结果能在所有预想故障场景下维持不失负荷但对每台机组添加一个“短时过载能力”参数在故障场景校验时允许机组出力在短时间内越过正常运行上限。这样既不会把模型搞成完全每个场景一套变量的庞大结构又能反映实际机组具备的短时调节能力。2. 光热电站不是“带电池的光伏”CSP的调度模型该怎么建2.1 光热电站的物理链条与三个子系统光热电站Concentrating Solar Power简称CSP在整个调度模型里需要拆成三个部分来看集热场、储热罐、汽轮发电机组。太阳直射辐射DNI照射到集热器上被转化为热能热能一部分直接送入蒸汽发生器发电另一部分存入熔盐罐需要发电时再从储热罐取出热量驱动汽轮机。在数学上这套能量链条可以写成三层关系。第一层是集热场的热功率上限Q_sf(t) η_sf × DNI(t) × A_sf其中A_sf是集热场有效采光面积η_sf是光热转换效率DNI是当前时段的直射辐照度。这一层决定了“太阳能进来多少”。第二层是储热罐的能量平衡E_tes(t) E_tes(t-1) (Q_sf(t) - Q_pb(t) - Q_spill(t)) × Δt这里Q_pb是送入汽轮机发电的热功率Q_spill是集热场超出需求后被舍弃的热功率用大白话说就是“收进来的热发掉一部分、存一部分、实在多到装不下就扔一部分”。第三层是热电转换P_csp(t) η_pb × Q_pb(t)η_pb是汽轮发电机组的热电转换效率。整条链路的物理意义可以类比成一个移动电源太阳能是充电器熔盐罐是电池汽轮机是用电器。充电速度受阳光限制电池容量有限放电速度还要受汽轮机自身特性的约束。调度模型要做的就是在这个链条的每一个环节上都写上对应约束。2.2 储热罐的时序耦合是最容易丢的约束我见过很多第一次建CSP模型的人把光热电站简单处理成三个不等式出力有上下限、有爬坡约束、有电量限制然后就没有然后了。这么做最直接的后果是储热罐的跨时段能量状态丢失了优化器可以“凭空”在某个时段给出很高的出力而根本不管之前有没有存够热。实际上储热罐是一个强时序耦合的状态量第t时段的储热水平由第t-1时段的状态和本时段的热输入输出共同决定。这正是CSP调度里最核心、也最容易被忽略的部分。在设计模型时一定要把E_tes作为决策变量放进优化问题并显式写出相邻时段的递推关系。同时还要加两条硬性约束第一条是储热罐容量上下界0 ≤ E_tes(t) ≤ C_tes其中C_tes由熔盐罐的物理容量决定。第二条是调度周期末段的储热水平约束通常写成E_tes(T) ≥ λ × C_tesλ取0.1到0.3之间。我最早跑24小时调度时没有加这条末期约束结果优化器把储热罐在夜里全部放空第二天早上系统又面临负荷高峰光热电站无法提供支撑出力曲线看起来就像“为了省弃热而透支了明天的能力”。加上末期约束之后结果立刻合理很多。2.3 CSP机组的运行区间与爬坡限制不是想发就发光热电站的汽轮机组本质上是一台热力机组它有最小技术出力有爬坡速率限制也有启停状态约束。很多简化模型把P_csp当连续变量用0到额定容量直接约束这在大多数情况下是过度乐观的。实际汽轮机在低负荷区域运行效率极差甚至可能无法稳定运行。我通常会给光热电站定义一个正常运行区间比如20%到100%额定出力低于20%时必须停机。同时爬坡速率也要写进去熔盐储热虽然能平抑太阳波动但汽轮机自身的进气阀门、缸体热应力决定了它不能像锂电池放电那样瞬间改变出力。实测下来一台50MW的CSP机组爬坡速率一般取每分钟2%到5%额定出力折合到小时级调度模型里就是每时段24到60MW/h的变化上限。如果你用的是小时级调度模型爬坡约束可能没有那么突出但启停状态一定要考虑。我的做法是用一个0/1变量u_csp(t)表示汽轮机是否在线然后把P_csp(t)写成u_csp(t) × P_csp_min ≤ P_csp(t) ≤ u_csp(t) × P_csp_max同时加上最小开机时间和最小停机时间约束否则优化器会为了某个时段的负荷波动而让汽轮机频繁启停这在物理上根本实现不了。这一步做完CSP才真正像一个可以参与调度调节的电源而不是一个“随叫随到的电热水壶”。3. 模型总装目标函数、运行约束与N-k故障场景的表达式3.1 目标函数怎么设置才能兼顾经济性和N-k惩罚我用的目标函数是系统总运行成本最小化包含三个层次min ∑ₜ [ ∑₉ (a_g × P_g(t)² b_g × P_g(t) c_g × u_g(t)) c_csp × P_csp(t) ] ∑ₛ π_s × ∑ₜ c_curtail × P_curtail(s,t)第一项是常规火电机组的发电成本通常用二次函数拟合煤耗第二项是光热电站的运行维护成本取值很小比如2到5美元/MWh主要让优化器不会在两种可行方案间随意抖动第三项是故障场景下的失负荷惩罚c_curtail要取得非常大远高于正常发电成本保证优化器只有在极端场景下才会选择切负荷。这里有一个值得注意的设计细节N-k安全约束本身并不直接进目标函数而是通过约束条件来体现。目标函数里的失负荷惩罚是为了在“故障场景确实无法完全满足负荷需求”时给优化器一条软性退路。比如一个故障场景把系统分裂成孤岛孤岛内发电容量不足这时硬性要求功率平衡且不失负荷优化问题就会直接无解。加入失负荷变量并配上高惩罚系数后问题永远有解同时优化器会尽量避免失负荷。3.2 基础运行约束的完整集合在写N-k场景约束之前先把正常状态下的基础约束列齐。我在Matlab实现中通常包含以下五组一是节点有功平衡约束对每个节点i和时段tΣ_g∈Gi P_g(t) P_csp(t) - P_load_i(t) Σ_l B_il × (θ_i(t) - θ_l(t))这里用的是直流潮流模型忽略无功和网损把节点注入有功和支路有功潮流用节点相角表示。二是机组出力上下界约束P_g_min × u_g(t) ≤ P_g(t) ≤ P_g_max × u_g(t)同时要有爬坡约束P_g(t) - P_g(t-1) ≤ RU_gP_g(t-1) - P_g(t) ≤ RD_g。三是机组最小启停时间约束这个需要一组二进制变量表示启停事件并用专门的不等式组实现。很多入门代码会漏掉这一条导致优化结果里同一台机组在相邻时段反复启停。四是线路潮流限值约束对每条支路l|B_l × (θ_from(l) - θ_to(l))| ≤ P_l_max。五是CSP子系统约束包括上一节写的集热场热功率约束、储热罐能量平衡、容量限额、末期储热约束以及汽轮机的出力区间和启停约束。3.3 N-k故障场景的建模预想故障集、孤岛处理与失负荷变量N-k约束的加入方式核心是对每个预想故障场景s重复一遍潮流安全约束只是系统拓扑变成了“移除故障元件之后的拓扑”。具体做法是预先枚举出预想故障集F对每个场景s∈F生成对应的节点电纳矩阵B_s。如果场景s涉及线路断线就把该线路从支路列表中剔除重新组装B矩阵如果涉及机组故障就限制该机组出力强制为0。然后对每个场景s添加以下约束Σ_g∈G_avail(s) P_g(t) P_csp(t) P_curtail(s,t) P_load(t)以及该场景下的支路潮流限额约束-B_l_max ≤ B_s × (θ_i(s,t) - θ_j(s,t)) ≤ B_l_max这里每个故障场景其实都有一套自己的相角变量θ(s,t)不能直接沿用正常状态的θ因为系统拓扑变了之后相角分布必须重新求解。我开始实现时图省事想用同一套θ变量同时校验所有场景结果发现一条线路断线后的相角物理上就和正常状态完全不同这样约束本身就没有意义了。正确的做法是给每个场景单独定义一组相角变量把拓扑变化的影响真正放进优化问题。失负荷变量P_curtail(s,t)是保证可行性的关键。故障场景s中如果可用发电容量确实小于负荷优化器会切掉一部分负荷但目标函数里对应的高惩罚会让它尽量不这么做。实际结果解读时可以统计每个场景的最大失负荷量以此判断系统当前调度计划的N-k安全裕度。关于预想故障集怎么选我的建议是分三层第一层是全部N-1场景数量少且物理意义明确第二层是高风险N-2场景比如同塔双回线、同母线多出线、共用通道的线路组第三层是根据基态潮流结果筛选的高负载率线路组合。不要一上来就把所有C(N,k)组合全塞进模型那样求解器会先把你自己的耐心耗光。4. 算例设计IEEE14节点做机理验证IEEE118节点做规模压力测试4.1 两套标准算例系统的基本差异算例选型上我同时跑IEEE14节点和IEEE118节点两套系统目的完全不同。IEEE14节点系统规模小14条母线、约5台机组、20条线路左右非常适合做机制验证。在这个系统上N-1场景全枚举只有20个左右N-2线路场景也不过190个全部塞进优化模型也能在几十秒内求解。我习惯用它来检查约束逻辑是否正确比如光热电站的热能平衡是否合理、故障场景下线路潮流是否真的没有越限、失负荷惩罚是否真的只在极端场景才触发。IEEE118节点系统规模大得多118条母线、54台机组左右、线路约186条是检验算法效率和可扩展性的“压力测试机”。在这个系统上如果仍然全枚举所有N-2场景就是C(186,2)约1.7万个场景直接建一个包含所有场景约束的MILP模型规模会膨胀到求解器难以接受的程度。我在118节点上只保留了全部N-1场景和约50个高风险N-2场景才把求解时间控制在可接受的分钟级范围。下表是我实现时关注的关键差异指标 | IEEE14节点 | IEEE118节点 母线数 | 14 | 118 机组数 | 约5 | 约54 支路数 | 约20 | 约186 N-1场景规模 | 几十个 | 两百多个 N-2全枚举难度 | 容易可直接枚举 | 组合爆炸必须筛选 主要用途 | 逻辑验证、教学演示 | 大规模算法评测、性能测试需要注意IEEE标准算例的不同版本在机组数、线路参数上会有差异我从Matpower的case14和case118数据文件读入数据时就发现机组报价参数很多版本里根本没提供需要自己根据煤耗特性补一套。算例数据的版本一致性会在结果对比时非常折腾建议在代码开头固定数据来源版本并写进注释。4.2 N-2场景选择策略与故障集设计IEEE118节点上N-2全枚举不现实我的筛选策略是按“故障后果严重度”来取舍。先用正常状态的最优调度结果做一次潮流计算列出负载率最高的前20条线路然后在这些线路之间以及它们周围电气距离较近的线路之间生成N-2组合。同时手动添加经典的高风险场景同一走廊的双回线、同一厂站出线、主变与其高压侧线路组合。这里有个容易被忽略的点N-k约束不是场景越多越好。把不相关的低风险场景加进模型不仅无助于提高安全性还会白白增加求解时间。我在118节点上做过对比把N-2场景从全部1.7万个缩减到前50个高负载场景后求解时间从数小时量级降到几分钟量级而最严重场景的失负荷指标几乎没有变化因为低风险场景本来就不会造成约束起作用。4.3 什么结果才算“调度模型真正起到了作用”算例跑完之后一定要看几个关键指标来判断模型有没有真正发挥作用。第一个指标是带N-k约束和不带N-k约束两种情况下的正常运行成本差。如果两种情况下成本几乎一样说明预想故障集选得太温和安全约束没有真正绷紧。正常情况下考虑N-k安全约束后运行成本会上升幅度一般在2%到10%这部分成本就是“买安全保险”的代价。第二个指标是故障场景下最大失负荷量。如果没有N-k约束某些高风险场景下可能直接出现数十MW的失负荷加上约束之后这个量应该显著下降或降为零。我在IEEE14节点上验证时最严重的N-2故障场景从失负荷约30MW降到0MW代价是正常状态总成本上升约4%这个权衡关系非常清楚。第三个指标是求解时间和MIP gap。IEEE118节点上如果gap设置过松比如1%优化器会提前终止此时看似有一种可行解但可能距离真正的最优解很远而且安全约束虽然满足运行成本却虚高不少。我的经验是把gap控制在0.1%以内求解时间允许放宽到10分钟。5. MatlabYALMIP代码结构从数据读入到Gurobi求解的完整骨架5.1 主程序的七步骨架与数据组织我整个代码按照下面这个流程组织每个子函数各自独立方便替换算例数据和故障集。以Matlab配合YALMIP工具箱为框架求解器用Gurobi或CPLEX。七步骨架是读数据、定义变量、收集正常态约束、收集N-k场景约束、设置目标函数、求解、后处理校验。主函数开头长这样function [result] csp_nk_schedule(casefile, K) mpc loadcase(casefile); % 读取IEEE标准算例 T 24; % 调度时段数 dt 1; % 时段时长小时 faultSet build_fault_set(mpc, K); % 构造预想故障集 %% 定义决策变量YALMIP Pg sdpvar(mpc.ng, T, full); % 常规机组出力 ug binvar(mpc.ng, T, full); % 机组启停状态 theta sdpvar(mpc.nb, T, full); % 正常态节点相角 theta_s cell(length(faultSet), 1); % 各故障场景相角 Pcsp sdpvar(1, T, full); % 光热电站出力 Etes sdpvar(1, T, full); % 储热罐能量 Pcurtail sdpvar(length(faultSet), T, full); % 场景失负荷 %% 收集约束 C []; C [C, BaseConstraints(...)]; for s 1:length(faultSet) theta_s{s} sdpvar(mpc.nb, T, full); C [C, FaultScenarioConstraints(...)]; end %% 目标函数与求解 objective ObjectiveFunc(Pg, ug, Pcsp, Pcurtail); ops sdpsettings(solver, gurobi, mipgap, 1e-4, verbose, 1); optimize(C, objective, ops); %% 回读并校验结果 result extract_result(...); end这段代码里有一个关键点每个故障场景都单独定义了theta_s{s}这一套相角变量绝不能复用正常状态的theta。我第一次实现时为了减少变量数偷懒复用theta结果约束虽然形式上写了实际却完全没有反映故障拓扑算出来的结果在N-k校验里一塌糊涂。5.2 YALMIP变量声明、约束收集与求解器配置YALMIP里最坑的一个细节是矩阵变量的默认对称性。sdpvar(n, T)默认声明一个n×T的对称矩阵这对调度问题来说是完全错误的。机组出力矩阵需要每一行是一台机组、每一列是一个时段所以必须写成sdpvar(ng, T, full)加上第三个参数full告诉YALMIP这是完整矩阵。我一开始就是漏了这个full导致模型维度莫名其妙扩大求解速度极慢还报出一堆“变量越界”的错误。约束收集时另一个值得注意的点是不要在循环里反复用方括号拼接超大的约束数组。场景数多、变量维数大时YALMIP内部每次拼接都会复制整个约束集内存开销非常离谱。我的做法是为每个故障场景先生成一个单独的约束元胞数组subC{s}全部循环结束后再用方括号一次性展开C_all [C_base, subC{:}];求解器配置上我比较常用的设置是ops sdpsettings(solver, gurobi, mipgap, 1e-4, ... mipfocus, 1, timelimit, 600, ... numericalemphasis, 1, verbose, 2);mipgap设成1e-4代表相对最优间隙到0.01%才停保证结果可信。numericalemphasis设为1让求解器在遇到数值病态问题时更积极处理对带多种约束尺度差异的调度问题很有帮助。5.3 一段可以直接参考的CSP储热约束代码下面这段代码是我在Matlab中实现光热电站核心约束的片段可以直接嵌到模型里。% 输入已知量 DNI load_dni_curve(); % 24小时DNI曲线W/m2 A_sf 2.5e6; % 集热场面积m2 eta_sf 0.4; % 光热转换效率 eta_pb 0.38; % 热电转换效率 C_tes 300; % 储热罐容量MWh P_csp_max 50; % 光热额定出力MW P_csp_min 10; % 光热最小技术出力MW lambda 0.2; % 末期储热比例 % 决策变量 % Q_sf_total 是集热场热功率Q_spill 是弃热 % Q_pb 是送入汽轮机的热功率 % 集热场热功率上限 C [C, Q_sf_total eta_sf * A_sf * DNI(t) / 1e6]; % 单位MW % 储热罐能量平衡 C [C, Etes(t) Etes(t-1) dt * (Q_sf_total - Q_pb - Q_spill)]; % 储热罐容量与末期约束 C [C, 0 Etes C_tes]; C [C, Etes(T) lambda * C_tes]; % 电热转换与机组运行区间 C [C, Pcsp(t) eta_pb * Q_pb(t)]; C [C, u_csp(t) * P_csp_min Pcsp(t) u_csp(t) * P_csp_max];如果你在实际项目中看到“光热出力曲线很漂亮但储热罐能量毫无变化”的结果那大概率就是储能平衡约束没写进模型或者在YALMIP变量声明时把Etes误声明成了对称矩阵。6. 求解性能优化与我最常遇到的几个报错陷阱6.1 组合爆炸场景筛选比硬解全部组合靠谱N-k约束的最大敌人永远是组合爆炸。我自己的一个深刻教训是在IEEE118节点上第一次尝试全枚举N-2场景代码写完后一跑发现YALMIP光是模型编译阶段就吃掉了近20GB内存Gurobi在半小时内一条可行解都没出来。后来改用“基态潮流预筛选”的方法只挑出线路负载率排行前20的线路及其组合模型规模瞬间缩小了两个数量级。场景筛选的落地做法是先用不带N-k约束的基本调度模型求一版结果然后对这个结果做全场景N-k潮流扫描把实际出现线路越限或母线电压问题的故障场景提取出来作为约束集加入优化模型再重新求解迭代2到3轮。这种“预调度—故障筛选—再调度”的外循环方式在工程上远比一次性枚举所有场景高效而且结果基本等价。不少商业调度软件内部也是这么做的。6.2 直流潮流导致的孤岛失负荷问题DC潮流模型把系统当作纯有功网络线路停运后如果系统被分成两个电气孤岛每个孤岛内部必须自己保持功率平衡。有些预想故障场景会产生完全没有电源的孤岛负荷比如某条线断开后一小片负荷区域与主网失去连接而该区域又没有本地机组。这种“孤儿负荷”场景放进模型后如果没有失负荷变量优化问题直接无解。我的处理方式是给所有负荷节点都加上故障场景失负荷变量P_curtail(s,t)并允许其在极端场景下切掉部分或全部负荷。目标函数里的失负荷惩罚取一个极大值正常场景下优化器不可能主动切负荷真遇到孤岛无电源的场景时切负荷变量兜住可行性模型仍然有解同时结果报告里会把这个孤岛失负荷量明确标出来。这一步做完N-k调度在DC潮流框架下的数值稳定性会好很多。6.3 结果不合理时的排查顺序跑出来的结果不对劲时我习惯按下面这个顺序排查效率高很多。如果你手头的模型也出现了类似问题可以直接照这个思路走症状 | 可能原因 | 排查方向 求解时间爆炸 | 故障场景全枚举过多 | 改用迭代场景筛选限制N-2组合范围 正常态和故障态结果几乎一致 | 故障场景约束没有真正生效 | 检查场景s的相角变量是否独立B_s是否重新组装 光热电站全天满发 | 储热末期约束缺失 | 加末期最低储热水平约束 机组频繁启停 | 最小启停时间缺失或gap过松 | 补启停时间约束收紧mipgap 无解且提示infeasible | 孤岛负荷无电源且无失负荷变量 | 给故障场景加失负荷变量和高惩罚系数 线路潮流普遍偏大但安全校验通过 | B矩阵组装错误 | 对比正常态B_s与去掉线路后的B_s差异最后再分享一个我自己调试时的高频操作任何时候怀疑N-k约束没生效最简单粗暴的验证方法是人为设置一个必会造成越限的故障场景比如把系统中负载最重的线路单独断开。如果约束正确模型要么自动调整机组出力消除越限要么在结果中给出该场景失负荷量如果两种迹象都没有说明某个环节的约束根本没有写进优化器。这一招在IEEE14节点上验证逻辑特别方便几乎一测一个准。从14节点搬到118节点的过程中我最大的感慨是让N-k安全约束真正落地难的从来不是约束公式本身而是场景管理和求解器的脾气。希望这篇笔记能帮你少踩几个我已经替你踩过的坑。
