接手风电、光伏、光热打捆外送相关的课题时，我第一个要搭的就是“含风电-光伏-光热电站电力系统N-k安全优化调度模型”的Matlab代码框架。这个模型听起来唬人，拆开看其实就三件事：怎么建模含储热的光热电站，怎么把N-k安全约束塞进优化问题，以及怎么让求解器在可接受时间内跑出结果。这篇文章把这三件事一次讲透，代码骨架可以直接改到IEEE 30节点或IEEE 118节点系统上，适合电气工程研究生、新能源并网算法工程师和做电力系统优化仿真的同学参考。

1. 从N-1到N-k：为什么风电光伏光热打捆系统需要更强的安全约束

1.1 单一故障校验的局限性

传统经济调度里，N-1校验已经是一个“默认安全底线”：任意一条线路或一台机组退出运行，系统仍能稳定运行。但实际电力系统事故很少是孤立的。线路连锁跳闸、新能源场站因低电压穿越失败成片脱网、极端天气导致多回线同时跳闸，这些都是真实存在的情况。如果一个调度模型只保证单一故障下的安全，那么多个元件同时停运时，系统潮流可能严重越限，甚至引发大面积停电。

高比例风电、光伏接入后，问题的性质又变了。风电和光伏出力随机波动，同一个断面下的输电通道负载率随新能源出力大幅变化。调度人员很难事先预料哪个时刻、哪组故障最危险。如果只做N-1，系统在“单重意外”下虽安全，但在“多重意外”下可能完全没有调整空间，因为基态调度已经把线路和机组都压得比较满。

这就是N-k安全约束的意义：要求调度结果在任意k个元件同时退出（或按N-1-1方式相继退出）时，系统仍能通过快速调整恢复到安全运行状态。这里的“k”可以是线路、变压器、发电机组，甚至整个新能源场站。

1.2 风光热打捆外送场景中N-k的必要性

我现在手里的仿真系统是一个典型的风-光-热打捆外送算例：风电场、光伏电站和光热电站集中接入一个区域电网，通过两条联络线与主网相连。如果按传统N-1去做安全校核，单独断开任何一条联络线，另一条线都能勉强承受剩余功率；看起来结果是合格的。

但把故障条件升级为N-2，两条联络线同时断开，剩余通道的负载率会瞬间冲到140%以上，系统根本扛不住。更麻烦的是，这个故障场景恰恰是实际中有可能发生的——联络线通常架在同一走廊上，恶劣天气下同跳的概率并不低。如果调度模型没有把N-2约束写进去，最后给出的机组出力和新能源并网计划，实际上是一次“隐藏着越限风险”的调度方案。

光热电站的作用在这个场景里特别有意思。它带储热罐，夜间也能发电，而且爬坡速率和调节范围都比火电灵活。如果优化模型给储热罐预留了足够的放热能力，那么在极端断面断开时，光热可以短时增发，减轻通道压力。反过来，如果模型只是简单把光热当成普通火电，这个调节能力就会被白白浪费。这就是为什么含光热系统的安全调度模型一定要把储热状态变量和N-k安全约束放在同一个优化框架里一起求解。

2. 风电、光伏与光热电站建模：出力序列和储热罐时间耦合是关键

2.1 风电和光伏出力序列的简化处理

风电出力在调度模型中通常用两种方式表达：一种是直接用预测功率曲线，另一种是用威布尔分布抽样生成若干典型场景。实际写代码时，我更倾向于用“可用功率上限+弃电惩罚”的确定性等效模型。原因很简单：N-k安全约束本身就包含大量故障态约束，模型规模已经不小，如果再叠加随机场景，求解难度成倍增加。

具体来说，风电的可用功率记为 Pw_avail(t)，光伏的可用功率记为 Ppv_avail(t)，决策变量 Pw(t) 和 Ppv(t) 满足：

• 0 ≤ Pw(t) ≤ Pw_avail(t)
• 0 ≤ Ppv(t) ≤ Ppv_avail(t)

实际并网功率小于可用功率的部分，在目标函数里以弃风弃光惩罚项体现。这个做法虽然没有严格处理随机变量的完备分布，但在工程规划中足够实用，而且代码实现非常直观。

光伏出力可以用光照强度Beta分布抽样，再乘以光电转换效率。不过如果做的是长时间尺度调度（比如24小时、1小时一个点），直接用典型晴天辐照曲线乘以装机容量和效率会更省事。输入数据只要保证时序上有“白天发电、夜间为零”的趋势，优化模型就能正确反映光伏对系统的影响。

2.2 光热电站不是“能烧锅炉的火电”：储热罐约束详解

这是整篇代码里最容易出错的地方，必须专门拿出来说。很多初学者把光热电站（CSP）建模成一台“燃料免费的火电”，只约束出力上下限，完全没有储热的概念。实际上，光热电站由集热场、储热系统、发电岛三部分组成，电出力只是热力过程的最终结果。

我在代码里采用了如下一组热功率平衡约束：

设 Q_sf(t) 为集热场在t时段收集的热功率： Q_sf(t) = η_coll · DNI(t) · A_field

其中 η_coll 是集热效率，DNI(t) 是法向直接辐射照度，A_field 是集热场面积。夜间 DNI 为零，所以 Q_sf 也为零。

集热场得到的热功率有四个去向：直接进发电岛发电 Q_pb(t)，充入储热罐 Q_ch(t)，从储热罐放出 Q_dis(t)，以及多余的热量弃掉 Q_dump(t)。热功率平衡为：

Q_sf(t) + Q_dis(t) = Q_pb(t) + Q_ch(t) + Q_dump(t)

储热罐的蓄热量 E(t) 是状态变量，满足时间耦合约束：

E(t+1) = E(t) + (η_ch · Q_ch(t) - Q_dis(t) / η_dis) · Δt

其中 η_ch 和 η_dis 分别是充热和放热效率。发电岛电出力：

P_csp(t) = η_pb · Q_pb(t)

约束条件还包括：

• 0 ≤ E(t) ≤ E_bar
• 0 ≤ Q_ch(t) ≤ Q_ch_bar
• 0 ≤ Q_dis(t) ≤ Q_dis_bar
• Q_pb_min ≤ Q_pb(t) ≤ Q_pb_max
• E(0) = E(T)，也就是一天调度周期内储热罐要回到初始状态

这组约束的核心在于：储热罐让光热电站具备了“跨时段”调节能力，白天光伏大发时，光热可以少发电、多充热；晚上光伏归零，光热放热发电填补负荷缺口。如果忽略储热罐状态变量，这个电站就真的和火电没有区别了。

2.3 同时充放热是否需要互补约束

很多同学会问：要不要加 Q_ch(t) · Q_dis(t) = 0 来禁止同时充放热？我的经验是：如果充热效率和放热效率都不等于1，同时充放热必然造成能量损耗，目标函数中的成本或能量惩罚会自然抑制这种行为，不需要额外加非线性互补约束。

但在实际调参时要注意，如果效率参数设置得太接近1，优化器可能会在局部时段出现“边充边放”的退化解。遇到这种情况，我通常把充热效率调低到0.9左右，放热效率调高到0.95左右，同时给储热罐增加一个轻微的自损系数 σ，比如：

E(t+1) = (1 - σ) · E(t) + (η_ch · Q_ch(t) - Q_dis(t) / η_dis) · Δt

加入自损系数后，储热罐里存着的热量即使不用也会缓慢损失，这样模型就会主动避免无意义的充放循环。

3. N-k安全约束的数学化落地：故障集生成、PTDF/LODF校验与最坏故障辨识

3.1 安全约束如何写进经济调度模型

N-k安全约束的本质是：在基态机组组合和出力确定之后，还要保证对于预定义故障集 K 中的每一个故障 k，系统都能在校正控制时段内把潮流调整回安全范围。这里我们采用的是“预防控制为主+校正控制为辅”的工程化处理方式，不在故障态重新做机组组合，而是允许发电机出力在故障后一定范围内调整。

数学上，正常态直流潮流为：

B0 · θ0 = P_inj0

故障态 k 的潮流方程为：

B_k · θ_k = P_inj,k

其中 B_k 是故障后网络的直流潮流矩阵，P_inj,k 是故障态节点注入功率。故障态注入相比正常态的变化主要来自两部分：跳闸机组出力归零、以及非故障机组为恢复安全而进行的调整 ΔPG_k。

因此，对每个故障 k，需要添加下面这几组约束：

• 潮流方程：B_k · θ_k = P_inj0 + C_gen · ΔPG_k - P_loss，其中 P_loss 表示故障中损失的出力
• 支路容量：-Limit ≤ T_k · θ_k ≤ Limit
• 调整量限幅：-ΔPG_max ≤ ΔPG_k ≤ ΔPG_max

值得强调的是，ΔPG_max 的取值要尽量符合实际。我通常把它设为机组额定出力的10%-20%，对应故障后10分钟内的AGC或一次调频能力。如果这个值设得过大，安全约束就会形同虚设；设得过小，模型又可能直接不可行。

3.2 故障集生成策略：从全枚举到关键断面筛选

故障集 K 的生成方式直接决定模型规模和求解速度。k=1时很简单，每条线路、每台机组都算一个故障，线性枚举即可。k=2时要枚举线-线、线-机、机-机组合，规模变成 O(n²)；对IEEE 30节点改造系统来说完全能接受，但到了IEEE 118甚至更大系统，就必须要筛选。

我实际使用的筛选方法是“两步走”。第一步用基态潮流跑一遍，把所有线路按负载率降序排列，只保留负载率超过70%的线路，再加上关键断面上的联络线和几台大容量机组，组成候选元件集合。第二步在这个候选集合内做N-k组合枚举，对每个组合用PTDF/LODF快速算一遍越限量，保留越限量最大的前S个故障加入主优化模型。

这样做背后有一个工程经验：N-k最危险的故障往往不是随机组合，而是集中在重载线路、关键断面和大容量机组上。把候选集合控制住，约束数量能降一个数量级，而漏掉关键故障的概率很低。

下面是一个故障集预筛的伪代码流程：

1. 运行不含安全约束的基础经济调度 2. 提取线路负载率，确定候选元件集合 C 3. 生成 C 的所有 k 阶组合（k ≤ max_k） 4. 对每个组合，用 LODF 计算各支路越限量 5. 计算越限量的最大值，按降序排序 6. 保留越限量最大的前 S 个组合 7. 把这 S 个故障加入主优化模型

3.3 用PTDF和LODF快速计算故障后潮流

故障潮流不一定每次都要重新求解 B_kθ_k = P，尤其是在预筛阶段，用开断分布因子（LODF）可以大幅度提速。

直流潮流模型下，线路 m 开断后，线路 l 上的潮流变化量可以用LODF近似：

Δf_l = LODF_{l,m} · f_m0

其中 f_m0 是线路 m 在基态下的潮流。LODF 可以由PTDF矩阵直接派生。MATLAB里如果已经有MATPOWER，可以直接提取系统数据，用 makeBdc 函数生成直流潮流B矩阵，再计算PTDF。

计算单条线路开断LODF的核心代码思路如下：

function LODF = compute_lodf(PTDF, line_m, line_l) % line_m 开断线路，line_l 被影响线路 % 直流潮流近似 ptdf_m_m = PTDF(line_m, line_m); if abs(1 - ptdf_m_m) < 1e-6 LODF = 0; % 分母接近零时视为无转移 else LODF = PTDF(line_l, line_m) / (1 - ptdf_m_m); end end

对于同时开断多条线路的情况，可以用LODF叠加近似。预筛阶段对精度要求不高，近似计算完全可以承担“找最坏候选”的角色；真正优化求解时，再对筛选出的少数严重故障建立精确的故障态潮流约束即可。

4. Matlab程序实现：基于YALMIP的N-k调度模型装配与求解思路

4.1 程序整体结构与文件划分

我没有把这个模型写成一个几百行的单文件脚本，而是拆成几个职责清晰的函数。这样做的好处是，你换一个算例系统、改一种故障集生成策略，只需要替换对应模块，不用动主框架。

我的文件清单如下：

main_scuc_nk.m % 主程序，控制故障集生成与安全约束迭代 load_case_data.m % 加载改进的IEEE节点系统数据 input_profile.m % 生成风电、光伏、光热出力序列 build_fault_set.m % 生成候选N-k故障集 build_base_model.m % 用YALMIP建立基态优化模型 add_nk_constraints.m % 向模型追加N-k安全约束 check_result.m % 事后详细潮流校验与指标输出

主程序的核心逻辑是安全约束生成（SCUC）迭代，不是一次性把所有N-k约束全塞进去。流程如下：

初始化：故障集为空，建立基态优化模型 循环： 1. 求解当前优化模型 2. 对候选故障集做快速潮流校验 3. 若无任何故障越限，则退出 4. 否则，把越限最严重的故障k加入模型 5. 重新求解

这个迭代思路在工程界非常成熟，名字叫“security constraint generation”，也就是安全约束生成。它的效率比一次性全枚举高出好几倍，尤其适合N-k这类约束组合爆炸的问题。

4.2 IEEE 30节点算例的数据改造

我的测试算例是在IEEE 30节点系统基础上改造出来的。原系统保留火电骨架，额外接入风电场和光伏电站，并把一台常规机组替换为光热电站。下面的表是常用的一组改造参数：

风电和光伏的可用功率序列我一般用随机种子固定后生成一次，存成Excel或.mat文件，方便重复实验结果。风电采用威布尔分布抽样再加平滑，光伏用Beta分布抽样辐照度，光热DNI用带有云层扰动系数的正弦曲线近似。

4.3 YALMIP建模核心变量与约束

决策变量分为基态变量和故障态变量两大类。基态变量包括火电开机状态、火电出力、风电光伏并网功率、光热充放热功率和储热罐状态；故障态变量包括每个故障对应的发电机出力调整量和故障态相角。

YALMIP定义决策变量的代码片段如下：

u = binvar(nG, T, 'full'); % 火电开机状态，0/1 Pg = sdpvar(nG, T, 'full'); % 火电出力 Pw = sdpvar(nW, T, 'full'); % 风电并网功率 Ppv = sdpvar(nPV, T, 'full'); % 光伏并网功率 Ecsp = sdpvar(T, 1, 'full'); % 储热罐蓄热量 Qch = sdpvar(T, 1, 'full'); % 充热功率 Qdis = sdpvar(T, 1, 'full'); % 放热功率 Qpb = sdpvar(T, 1, 'full'); % 进入发电岛的热功率

目标函数我采用最小化总费用，包含火电燃料成本、启停成本、弃风弃光惩罚、失负荷惩罚和光热弃热惩罚：

% 火电燃料成本，二次成本按分段线性化处理 fuel_cost = sum(sum(repmat(cost_a', T, 1) .* Pg.^2 ... + repmat(cost_b', T, 1) .* Pg ... + repmat(cost_c', T, 1) .* u)); % 弃风弃光惩罚 renew_curtail = sum(sum(Pw_avail - Pw)) + sum(sum(Ppv_avail - Ppv)); % 光热弃热惩罚 csp_dump = sum(Qdump); obj = fuel_cost + start_cost + 200 * renew_curtail + 5000 * load_shed_total + 20 * csp_dump;

注意火电燃料成本里如果有二次项，YALMIP会把模型交给求解器的MIQP求解器。如果算例规模上去了，建议把二次成本改成三分段或四段线性化，把模型降成MILP，求解速度会有质的提升。

功率平衡约束、火电出力上下限、爬坡约束、备用约束这些都比较常规，不再逐条展开。光热储热罐约束是关键，前面已经推导过公式，YALMIP里直接按公式逐行写即可。

4.4 N-k约束追加的具体写法

在安全约束生成循环里，每轮迭代只会加入一个或少数几个严重故障，所以故障态变量的数量是可控的。YALMIP代码可以这样组织：

function cons = add_nk_constraints(cons, Bk, Ptrans, Pg, u, fault) % Bk: 故障态直流潮流矩阵 % Ptrans: 故障态节点注入修正 % fault: 故障描述结构体 [nG, T] = size(Pg); dPg = sdpvar(nG, T, 'full'); th_k = sdpvar(nBus-1, T, 'full'); cons = [cons, -deltaMax <= dPg <= deltaMax]; cons = [cons, Bk * th_k == Ptrans + C_gen * dPg]; cons = [cons, -LineCap <= Tline * th_k <= LineCap]; end

故障态变量 dPg 是连续变量，不涉及整数决策，所以即使故障集里有几十上百个故障，模型的总变量规模也不会失控。整数变量仍然只在基态里头，这是工程应用里很重要的一条经验：不要在故障态里重新安排机组启停，否则MIP求解难度会指数级上升。

4.5 结果校验与可视化

求解完成后，建议不要直接相信优化结果，要跑一次“事后校核”。我把最优解带入所有候选故障，重新计算精确的直流潮流，统计每个故障下的最大越限量。这一步是防污染的，因为YALMIP模型里我们可能用了预筛、近似LODF或分段线性化，最终潮流应该用原始数据复核。

可视化部分我会输出下面几张图：火电、风电、光伏、光热逐时段出力堆叠图，储热罐蓄热量曲线，以及各个故障下线路最大负载率柱状图。绘图用MATLAB内置的plot和bar就够。

5. 求解性能优化、参数调节与实测中的避坑经验

5.1 模型不可行的常见原因排查

我刚把N-k约束接上模型时，经常遇到的问题是“infeasible problem”。排错顺序很重要，我一般是按下面这个顺序查：

第一，故障后系统发电容量是否满足负荷。如果一个故障同时切掉了多台大机组，故障态的系统总可用容量低于负荷，那无论怎么调都不可能满足潮流约束。这时候要么增加备用约束，要么允许故障态少量切负荷，并在目标函数里给失负荷惩罚。

第二，储热罐终值约束是否过强。E(T) = E(0) 在夜间负荷高峰时可能锁死储热罐的放热空间。一个更稳妥的做法是把等式改成 E(T) ≥ 0.1 · E_bar，给模型一点回旋余地，代价是让一天内的充放策略略偏保守。

第三，线路容量数据本身是否合理。加入N-2故障后，某些线路会在单回线开断时承担双倍潮流，如果容量瓶颈不在模型而在原始数据，那应该调整线路容量或者改变断面结构，而不是硬调惩罚系数。

5.2 求解器选择与参数配置

YALMIP默认调用sedumi或linprog这类通用求解器，算IEEE 30节点、24时段、带少量N-2约束的模型还能跑，但一旦故障集扩大到几十个，速度就完全不能看了。我现在的标准配置是Gurobi或CPLEX，YALMIP只需要设置一行options。

典型求解参数如下：

options = sdpsettings('solver', 'gurobi', ... 'gurobi.MIPGap', 0.01, ... 'gurobi.TimeLimit', 600, ... 'verbose', 2);

MIPGap设成1%通常足够。如果你只是做论文里的算例对比，甚至可以放松到2%-3%，速度会快很多。但如果是做实际调度方案，建议严格控制在0.5%以内，因为安全约束对结果非常敏感。

如果同一套系统要跑几十组不同故障集算例，我还推荐用“热启动”技巧：先用N-1约束的模型求解一次，把机组组合结果作为N-k模型的初始解。这个技巧在Gurobi里直接用start值即可实现，能减少不少时间。

5.3 三个实测中得到的反直觉结论

第一个结论：最严重的N-2故障，往往不是负载率最高的两条线路同时断开。我的算例里有一条负载率92%的线路和一条负载率78%的线路，直觉告诉我这两条一起断开最危险。但实际最坏故障是两条负载率分别为75%和69%的并行断面线路同时断开，因为它们的断开导致第三条联络线转移了几乎全部断面功率，负载率冲到143%。这提醒我，构建候选故障集时一定要把“位于同一输电断面”的线路作为强候选，而不是单纯看单条线路负载率。

第二个结论：加入N-2安全约束后，光热电站的充放热策略会发生本质变化。没有N-2约束时，模型倾向于让光热在午间多充热、晚高峰放热，这符合直观。加上N-2约束后，模型开始主动在负荷低谷期额外充一些热，并在某些关键时段预留放热能力，相当于给系统增加了一个“隐性备用”。这个效果只有把储热罐状态变量和N-k约束同时放进模型才能体现出来，也是光热相对其他新能源最独特的价值。

第三个结论：弃风弃光惩罚系数不是越大越好。在N-k约束下，风电光伏出力越高，系统基态越接近运行边界，故障后调整空间越小。有时把弃电惩罚系数从200提到500，系统为了满足N-2安全约束，反而会主动降低新能源并网功率。所以只看弃电率指标是不全面的，还要同时关注安全裕度。我后来习惯在结果表里同时输出“最坏故障下的最大线路负载率”，这个指标比成本更能说明模型质量。

5.4 不同k值下的求解规模对比

下表是我在IEEE 30节点改造系统上的一组实测结果，故障集生成参数相同，只改变最大k值：

这里的最坏故障越限量是指在“只考虑N-1约束的调度结果”下，用N-2故障集去校核得到的越限量。可以看到，N-1调度结果在N-2故障下有62MW的越限，说明如果不做N-2约束，系统的安全裕度是明显不足的。

值得说明的是这张表是特定算例下的相对结果，不是标准值。不同系统、不同故障集筛选策略差异会很大，但它有助于理解：N-k约束不是简单翻倍增加计算量，而是非线性上涨的。所以控制候选故障集规模比堆求解器更重要。

6. 光热电站容量与储热时长选择的敏感性分析

6.1 储热时长对N-k可行的实际影响

既然前面已经强调了储热罐的重要性，这一步值得单独拿出来量化。我做了这样一组实验：光热电站装机保持50MW不变，把储热时长从2小时逐步调到6小时，分别跑N-2安全约束模型，记录系统总运行成本和最坏故障下的最大线路负载率。

最明显的现象是：储热时长从2小时增加到4小时，最坏故障负载率下降14个百分点，效果非常显著；但从4小时到6小时，改善幅度变小。这说明在本文算例里，4小时储热基本是性价比拐点。遇到实际项目，我建议都做一次这样的敏感性曲线，不要凭经验拍脑袋定储热参数。

6.2 光热参与故障态校正控制的调度效果

还有一个容易被忽视的问题：光热电站是否应该参与故障态出力调整。有些模型把光热当作风电场一样只给固定预测出力，不参与调整，这样会浪费它的快速调节能力。

我在代码里给光热也加了故障态调整变量 dPg_csp，限制在额定出力的30%以内。实测下来，含光热参与校正控制的模型，N-2最坏故障越限量从77MW降到了52MW；同时系统总成本没有明显上升，因为光热只是短时增发，并不改变基态的开机方式。

这也解释了为什么优化模型需要为故障态预留“可调整空间”。如果调度计划把所有机组都顶到满发，连光热的储热罐也用得干干净净，故障来临时就算控制器再聪明，也没有物理功率可以去支援。

7. 把模型改到更大系统的三点提醒

7.1 故障集生成必须从“全枚举”切换到“安全约束生成”

IEEE 118节点系统有186条支路，N-2组合数大约1.7万个。如果每个故障都建一组约束，哪怕每组只有几十行约束，YALMIP建模矩阵也会膨胀到几GB。这时候必须把安全约束生成循环作为主框架，每次迭代只加入越限最严重的少数故障。这个方法我在前面已经介绍过，但到了更大系统，它不再是“优化项”，而是“必须项”。

7.2 用直流潮流近似要小心弱环网

直流潮流在辐射状或弱环网上精度尚可，但在强环网、高阻抗比线路多的情况下误差会被放大。N-k故障校核阶段，我建议对事故筛选用的是PTDF/LODF快速近似，但对最终入选的最坏故障，要回到交流潮流或至少保留完整B矩阵的直流潮流重新校核一遍。这个过程虽然费时，却是论文审稿和工程验收都认的可靠做法。

7.3 计算时间不够时优先压缩时间粒度

如果你需要跑24小时以上的中长期调度，而N-k约束导致单时段模型已经很大，可以把时间粒度从1小时放宽到2小时或4小时，先验证故障集的逻辑是否合理，再逐步细化。反过来，如果只关注某几个关键时段，可以只对这几个时段加入详细N-k约束，其他时段用简化的N-1或纯基态约束，这样既能保留安全特征，又不至于让模型无法求解。

这套代码我从N-1调到N-2再尝试关键N-3，绕了不少弯，最大的体会是：模型复杂度可以慢慢加，但数据准备的坑必须一开始就避开。拿到一个新算例系统，我先跑一次纯经济调度，把潮流分布摸清楚，标出重载线路和关键断面，再上N-2约束。这样一旦模型不可行，排错范围会小很多。另外一个习惯是，每次调完参数都顺手保存一组完整结果到mat文件，包括目标函数值、故障集、线路负载率矩阵。做敏感性分析或者写报告时，这些历史结果能省掉大量重复计算。如果你的系统里也有光热电站，建议做一次不同储热时长和不同故障集规模下的对照实验，这两个维度的结果往往比单纯调惩罚系数更能说明模型的工程价值。