含新能源的N-k安全约束经济调度：建模与MATLAB实现

做电力系统调度优化的这些年N-k安全约束是我觉得最难跟外行讲清楚、也是实际工程里最见真章的一个概念。随着风电、光伏、光热在电网里的渗透率越来越高调度模型再拿N-1当唯一安全标准说实话已经不太够用了。我近期在复现和研究一个带风电-光伏-光热电站的N-k安全经济调度模型用MATLAB实现了整套代码把建模思路、故障集生成方式、求解器配合以及我踩过的几个大坑一次性写出来给准备做这个方向或者正在被论文代码折磨的师弟师妹做个参考。这类模型不是那种“跑通就行”的玩具代码。你要真把它用到IEEE 30节点、IEEE 118节点甚至更大系统里就会遇到组合爆炸、无可行解、SOC漂移一整套问题。这篇文章我会从模型拆解、三类电源数学建模、N-k场景怎么生成、MATLAB实现框架一直讲到结果分析和排查思路全程按我实际调试的经验来讲不是教科书式的抄公式。1. 先拆清楚这个模型到底在算什么1.1 N-k安全约束到底是什么先说N-k最简单的理解。N代表系统里元件的总数一条线路算一个元件、一台变压器算一个元件、一台发电机也算一个元件。N-1是指任意一台元件或线路退出运行后系统仍然能安全稳定运行。N-k则是更狠的条件同时退出k个元件系统也不能出问题。传统的调度模型里安全校核做到N-1已经算不错了。但实际电网中极端天气、连锁故障、保护拒动这些情况往往不是单个元件故障那么简单。两个、三个元件同时出问题在运行层面是真实存在的。N-k安全约束就是把这种“多重故障”直接放进了优化模型里要求调度方案在正常态下做出来之后任何一个预想的k重故障场景下都不会出现线路越限、失负荷这类问题。把这个约束写进优化调度模型本质上是让决策者在“成本经济性”和“运行安全性”之间做一个数学意义上的折中。追求高安全必然要留更多备用、调整机组出力成本会上升完全不管安全成本是低了但故障一来就是连锁反应。这也是这类模型最有研究价值的地方。1.2 风电、光伏和光热在调度里的定位差异很多人第一次看到“含风电-光伏-光热电站”这个组合会以为光热跟光伏差不多都是看天吃饭。这个理解必须纠正否则后面建模全乱。风电和光伏本质上是被动电源。风来了、阳光来了能发就发不能发就瞪眼。它们参与调度的方式主要是“弃风弃光”和“预测出力约束”——也就是说模型里风电光伏的实际出力不能超过当前时刻预测的可用功率你可以少用它们但不可能让它们多发电。光热电站完全不一样。光热电站的核心结构是太阳岛集热场、储热系统熔盐罐和动力岛汽轮发电机组。太阳能先被集热场收集转化成热能加热熔盐热能可以立即用来产生蒸汽发电也可以先存在储热罐里留着晚上或者阴天再发。这相当于给可再生能源装了一个“充电宝”而且是热形式的充电宝。所以在优化模型里光热电站的出力不是完全由当前光照决定的它有一个腾挪空间光照充足时多收集热量、把热量存起来光照不足或者系统需要的时候再放热发电。它既像新能源又具备类似火电的可调度性。这个特性在N-k安全约束下尤其有价值——故障发生后光热电站可以快速响应承担一部分功率支撑。1.3 优化调度模型的目标函数这类模型的目标函数通常不是单一的“发电成本最小”而是一个包含多项成本的综合目标。以我的实现为例目标函数包括四块第一是常规火电机组的煤耗成本通常用二次函数表示也就是出力P的二次项加一次项再加固定成本。第二是弃风弃光惩罚风电光伏的预测可用功率没有被充分利用时要付出惩罚成本这个惩罚系数要设置得比火电边际成本高才能让模型优先消纳新能源。第三是光热电站的运行成本包括集热场维护、熔盐泵耗电、启停费用等这部分成本相对低但必须加上否则模型会让光热电站疯狂出力。第四是N-k安全约束带来的隐性成本——比如为了满足故障态约束正常态下某些火电必须保底运行、机组出力必须压低这些都会反映在总成本里。这里要强调一点如果只做单目标优化N-k安全约束是以“硬约束”形式存在的。你可以通过调整故障集的大小来观察它对成本的影响而不是把安全性写进目标函数做多目标优化。这样求解相对简单结果也更容易解释。2. 三类电源的数学建模2.1 风电、光伏的出力约束风电和光伏在调度模型里属于“可用出力约束型”电源。关键变量有两个预测可用功率和实际调度出力。预测可用功率来自日前预测比如风电在某个时段预测可用功率是200MW但实际上电网可能只需要它发150MW那50MW就是弃风。数学上写出来就是[ 0 \le P_{w,t} \le P_{w,t}^{avail} ] [ 0 \le P_{pv,t} \le P_{pv,t}^{avail} ]其中[P_{w,t}^{avail}]和[P_{pv,t}^{avail}]是预测值是已知参数。如果你的模型是多时段的还要考虑风电光伏出力的时序耦合但我这里先讲单时段版本——也就是典型的经济调度切片后面扩展多时段时再叠加爬坡约束和储热SOC的时间递推关系。有人会问要不要用场景法描述风电不确定性当然可以。比较完整的做法是用蒙特卡洛生成大量风电出力场景再通过K-means或者同步回代缩减出典型场景把确定性模型变成多场景随机优化。这样做当然更严谨但组合N-k故障场景一起用问题规模会爆炸。我的经验是第一次实现N-k模型先别急着做随机优化用确定性预测值把安全约束的逻辑跑通再回头升级不确定性模块。这样每步的坑都能单独排查。2.2 光热电站的储热细节光热电站的建模是这个模型里最需要仔细处理的部分。完整的CSP电站模型有三块集热场、储热系统、发电系统。集热场的产热功率由太阳法向直接辐射DNI决定近似写成[ Q_{sf,t} \eta_{sf} \cdot DNI_t \cdot A_{sf} ]这个热量有两条出路直接送去蒸汽发生器发电或者送入储热罐存起来。储热罐的能量平衡关系是[ E_{t1} E_t \eta_{ch} \cdot P_{ch,t} - \frac{P_{dis,t}}{\eta_{dis}} ][ E_{min} \le E_t \le E_{max} ]其中[E_t]是储热罐的当前储热量SOC[P_{ch,t}]是储热功率[P_{dis,t}]是放热功率[η_ch]和[η_dis]分别是充放热效率。发电系统的功率输出取决于输入蒸汽发生器的热量加上从储热罐释放的热量再乘以热电转换效率。所以光热的净出力约束可以写成[ P_{csp,t} \eta_{sg} \cdot (Q_{sf,t} P_{dis,t} - P_{ch,t}) ][ 0 \le P_{csp,t} \le P_{csp}^{max} ]这几个式子放在一起光热电站的调度自由度就体现出来了某时刻DNI很高产热很多但如果此时电价不高或者系统不缺电可以选择少发电、多充储热把热量留到晚上系统负荷高峰再放出来。在N-k安全约束下这种“时间搬移”能力非常有用故障时段如果有额外功率需求光热可以释放储热顶上。注意一点储热罐的充放热不能同时进行实际运行中熔盐泵和换热器的工作模式决定了充放互斥。严格建模要加一个0-1变量[ P_{ch,t} \le M \cdot z_t, \quad P_{dis,t} \le M \cdot (1-z_t), \quad z_t \in {0,1} ]如果不加这个约束模型会利用充放同时进行的漏洞白赚效率结果就是SOC一直虚高实际根本做不到。这个细节我在第一次实现时踩了坑后面在常见问题部分再详细说。2.3 网络潮流约束DC潮流和PTDF安全约束的核心是线路潮流。交流潮流是非线性模型放进N-k优化里基本没法大规模求解所以工程和学术上都用直流潮流近似。DC潮流假设电压幅值接近1、线路两端相角差很小、线路电阻远小于电抗从而把非线性潮流简化成线性方程[ P_{line} PTDF \cdot (P_{injection} - P_{demand}) ]PTDF全称是功率传输分布因子它表示某个节点注入单位功率变化时各条线路上潮流的变化量。给定网络拓扑和线路参数PTDF矩阵可以一次性求出来后面所有故障场景只需要变化拓扑再重算PTDF或者用线路开断分布因子LODF来修正。正常状态下线路潮流直接由PTDF矩阵和节点注入功率决定[ P_l^0 PTDF_0 \cdot (P_G P_{wind} P_{pv} P_{csp} - P_D) ]每个节点要满足功率平衡也就是注入功率等于负荷加上对外送电。在数学上节点注入功率等于该节点连接的发电机出力总和减去该节点的负荷。模型里不会再单独写每个节点的功率平衡因为PTDF已经隐含了系统级平衡。但你需要在变量定义层面保证系统总发电等于总负荷否则PTDF算出来的是不平衡解没有意义。3. N-k预想故障集与安全约束写法3.1 故障集生成枚举与筛选N-k模型第一步就是生成预想故障集。给定系统N个元件通常取线路和发电机要选k个同时故障组合数是C(N,k)。以IEEE 30节点系统为例支路大约41条加上6台发电机一共47个元件C(47,2)是1081C(47,3)就接近一万六。看起来还能算但如果系统换到IEEE 118节点支路186条加54台机组C(240,3)超过227万直接枚举再逐个加约束模型就没法解了。所以实际做法分几个层级。最简单的用Matlab的nchoosek直接枚举适合k1、k2的小系统或者教学演示。稍微复杂一点先根据电气距离、断面负荷水平、历史故障概率把故障组合筛选出“高风险集合”只对高风险集合施加安全约束。再进一步用优化-检验迭代法先不加N-k约束求解把结果代入故障校验模块找出越限最严重的故障场景加入约束再重新求解反复迭代直到没有新增越限场景。我在代码里实现的是枚举加约束生成混合策略。k2时直接枚举所有线路两两组合k3时只枚举包含关键断面的组合。这样既能保证论文复现的完整性又能控制求解时间。3.2 故障状态下线路潮流怎么算得到故障集之后核心问题变成某个故障场景下各线路潮流是多少这里有一个关键选择故障后机组出力变不变。如果采用严格的预防性调度认为调度指令已经在正常态下固定故障后短时间内机组不重新调度那故障态潮流只取决于网络拓扑变化和节点注入变化。这种假设下用LODF可以直接算[ P_l^s P_l^0 \sum_{m \in F} LODF_{l,m} \cdot P_m^0 ]LODF是线路开断分布因子它表示线路m开断后潮流转移到线路l的比例。这个公式只处理线路故障不处理发电机故障。如果故障集里包含发电机就需要先修正节点注入故障机组出力变为0其功率缺额按参与因子分配到其他机组然后再用故障后的PTDF矩阵计算潮流。我实际用的是更稳妥也更好调试的方法对每个故障场景直接用故障后的拓扑重新计算PTDF矩阵再用故障后的节点注入向量乘PTDF。虽然每次求逆耗时多一点但逻辑清晰不容易出错。MATLAB里这一步可以直接用Matpower的makePTDF函数传入修改后的branch矩阵。3.3 预防性安全约束的完整写法在所有变量和潮流表达准备好之后N-k安全约束就写得很直白了。对故障集里的每一个场景s要求[ -P_{l,max} \le P_l^s \le P_{l,max}, \quad \forall l ]其中[P_l^s]是场景s下线路l的潮流由场景相关PTDF和调度变量线性表示。每个场景加两组不等式即可。这里要注意的是“正常态变量和故障态变量是不是同一套”。我在代码里采用的是同一套调度变量也就是预防性调度。这样可以保证正常态做的决策在故障态下不需要重新调度就能满足安全要求。这也是N-k安全约束的强形式工程上最常用的解释是调度员“预料到可能发生的故障并且提前安排了一个即使故障发生也不用动的运行点”。如果你想把模型扩展成“故障后允许校正调度”那就需要引入第二层决策变量变成两个层级的鲁棒优化或者MPC滚动优化模型复杂度完全上一个量级。我的建议是先把预防性模型吃透校正式模型是在这基础上做后备、调频、切负荷的扩展。4. MATLAB代码实现解析4.1 整体框架选择我的实现思路是Matpower负责读数据和计算PTDFYalmip负责建模Gurobi负责求解。这套组合在电力系统优化里已经算标准配置原因是分工清晰Matpower的case格式是公开通用的Yalmip的约束写法接近数学表达式不容易出错Gurobi对大规模线性规划和二次规划的性能足够强。如果你没有Gurobi用Cplex或者开源求解器SCIP也可以。Yalmip对这三者都有接口。唯一要提醒的是N-k场景混合整数变量一旦多了用免费求解器很容易出现“几分钟不出结果”的情况。我的习惯是先把k1跑通验证模型正确性再逐步加k2场景这个过程里不要一下子把所有故障集都塞进去。4.2 数据准备与PTDF计算数据准备阶段我建议手动构造一个3节点或者5节点的简单系统来验证代码然后再切到IEEE 30节点。原因是N-k多故障场景下手算都不容易直接上大系统一旦出错你根本分辨不清是潮流算错、约束写错还是求解器的问题。下面是我在小系统上验证PTDF计算的一段核心代码mpc loadcase(case30); % 读取IEEE 30节点数据 baseMVA mpc.baseMVA; bus mpc.bus; branch mpc.branch; gen mpc.gen; [PTDF0] makePTDF(baseMVA, bus, branch); % 正常态PTDF n_line size(branch, 1); % 线路条数 % 遍历所有N-1支路开断 PTDF_cell cell(n_line, 1); for s 1:n_line branch_s branch; branch_s(s, :) []; % 断开第s条线路 PTDF_cell{s} makePTDF(baseMVA, bus, branch_s); end这里有个细节makePTDF要求传入的bus和branch矩阵是Matpower格式断开支路时直接删行但保留节点编号不变Matpower会自动按剩余支路重新编号内部导纳矩阵。多故障场景同理删除多行支路再调makePTDF。4.3 Yalmip建模关键代码变量定义部分我用Yalmip的sdpvar来建连续变量用binvar建0-1变量。核心片段Pg sdpvar(n_gen, 1); % 火电出力 Pwind sdpvar(n_wind, 1); % 风电实际出力 Ppv sdpvar(n_pv, 1); % 光伏实际出力 Pcsp sdpvar(n_csp, 1); % 光热电站净出力 Pch sdpvar(n_csp, 1); % 储热充电功率耗热 Pdis sdpvar(n_csp, 1); % 储热放电功率 SOC sdpvar(n_csp, 1); % 储热罐SOC % 目标煤耗成本 弃风弃光惩罚 光热运行成本 fuel_cost sum(a .* Pg.^2 b .* Pg c); curtail_cost 50 * sum(Pwind_avail - Pwind) 40 * sum(Ppv_avail - Ppv); csp_cost 5 * sum(Pcsp); Objective fuel_cost curtail_cost csp_cost;约束组装也直接逐条写Constraints [sum(Pg) sum(Pwind) sum(Ppv) sum(Pcsp) total_demand]; % 火电约束 Constraints [Constraints, Pgmin Pg Pgmax]; % 出力上下限 % 风电光伏出力约束 Constraints [Constraints, 0 Pwind Pwind_avail]; Constraints [Constraints, 0 Ppv Ppv_avail]; % 光热约束 Constraints [Constraints, 0 Pdis Pdis_max]; Constraints [Constraints, 0 Pch Pch_max]; Constraints [Constraints, SOCmin SOC SOCmax]; Constraints [Constraints, SOC SOC0 eta_ch*Pch - Pdis/eta_dis];光热净出力与热量的关系我在代码里写成% Q_sf 是集热场产热来自DNI数据 Constraints [Constraints, Pcsp eta_sg * (Q_sf Pdis - Pch)]; Constraints [Constraints, 0 Pcsp Pcsp_max];接下来是N-k安全约束的组装这是全代码最费内存的地方for s 1:n_scenario PTDF_s PTDF_cell{s}; % 节点注入向量发电机节点出力 - 负荷节点功率 Injection_s gen_to_bus * Pg wind_to_bus * Pwind ... pv_to_bus * Ppv csp_to_bus * Pcsp - Pd; Pflow_s PTDF_s * Injection_s; Constraints [Constraints, -line_max Pflow_s line_max]; end这里的gen_to_bus是发电机到节点的关联矩阵维度是n_bus乘n_gen。每台发电机实际接入哪个节点从mpc.gen的第一列GEN_BUS读取后转成稀疏矩阵。这个关联矩阵一定要仔细检查我见过太多人在这里犯迷糊——PTDF的维度是线路数乘节点数注入向量的维度是节点数乘1两者维度对不上整个模型第一跑就崩。4.4 求解与结果输出求解部分用Yalmip的optimize函数参数设置最关键的是gap和time limitoptions sdpsettings(solver,gurobi, ... gurobi.MIPGap, 1e-4, ... gurobi.TimeLimit, 3600); diagnosis optimize(Constraints, Objective, options);跑完后用value()取出各变量值再用正常态PTDF算一遍潮流人工检查一下有没有线路达到上限。我曾经遇到过求解器返回成功但线路潮流约束实际没有被满足的情况后来发现是PTDF矩阵的单位错了——baseMVA没有统一导致潮流结果小了100倍。所以我说结果校验这步绝不能省。5. 结果怎么分析才有说服力5.1 成本与安全之间的权衡模型跑通之后第一个要做对比实验就是设置不同的k值k0、k1、k2分别求最优调度成本。你会看到成本随k值上升而上升这是必然的。关键是上升的幅度和结构。以IEEE 30节点系统为例如果成本上升不到2%说明系统本身裕度充足如果上升超过10%甚至出现无可行解说明系统在特定断面上确实存在安全短板。这种分析用一张表展示最直观方案故障集规模总成本弃风弃光量光热SOC均值k00基准基准基准k1线路数略增可能下降明显提高k2C(线路,2)明显上升下降明显进一步上升光热SOC均值上升这个指标值得注意。当系统加了N-k安全约束之后光热电站会倾向于提前多存热量因为故障随时可能发生储热就是安全裕度。量化这个变化就相当于证明了光热储热在安全调度中的价值。5.2 光热“有储热”和“无储热”的对照实验要证明光热电站的价值最干净的做法是做一个对照组一个版本光热电站带储热罐一个版本去掉储热罐光热只能即时生产即时发电。两个版本加同样的N-k约束比较总成本和安全性。我在实际跑下来发现带储热的光热电站能显著降低N-k约束带来的成本增量。原因是它能把正午时段的过剩热能转移到傍晚负荷高峰而且故障需要功率支撑时储热罐是随时可用的备用容量。无储热版本的光热下午光照一弱就完全帮不上忙系统只能让火电顶着成本自然高。这种对照实验写论文出来特别有说服力建议大家都做一下。5.3 灵敏度分析除了k值还可以做两类重要的灵敏度分析。第一类是故障集里的元件类型选择只考虑线路故障和同时考虑线路加发电机故障结果会有明显差别后者约束更强成本更高。第二类是光热储热容量把储热罐的容量从4小时逐步增加到8小时、12小时观察总成本变化曲线。容量增加到一定程度后成本下降会趋于平缓这时候就能找到“储热容量最经济区间”。风功率预测偏差的灵敏度也很值得做。我在模型里把预测出力乘以不同的误差比例看调度方案成本的变化。这个分析能让你的模型更有实际意义因为N-k安全约束假设了预测值准确而现实中预测不可能准确。6. 我踩过的坑和排查思路6.1 模型报“无可行解”怎么排查N-k模型一上来就无可行解是频率最高的问题。我的排查顺序是第一步关掉所有N-k约束只跑正常态调度如果没有可行解说明基础约束就错了优先检查功率平衡等式和机组上下限。第二步加N-1约束看是否可行如果无解大概率是某条线路故障后潮流越限无法消除可以把最严重的故障线路打印出来。第三步加N-2约束如果无解往往不是全部N-2都无解而是个别断面组合太极端。实际工程里如果某些故障场景确实无法满足有一个变通方法允许这些场景下的负荷削减并在目标函数里加上失负荷惩罚。这相当于把“硬安全约束”降级为“带惩罚的软约束”在极端场景下牺牲一点负荷换取模型整体可行。论文里可以写成切负荷变量和切负荷成本。6.2 求解时间爆炸怎么办N-k模型最常见的性能瓶颈不是变量数而是故障场景数。每个场景加一组线路潮流约束场景上了几百个之后Yalmip内部展开的约束规模会非常大矩阵组装耗内存、求解器预处理也慢。我的解决办法有三个。第一尽可能把共享项提出来——比如正常态潮流约束单独写故障态约束只加故障后变化量部分减少重复计算。第二采用约束生成迭代而非一次性加全部约束。第三把故障场景先做筛选只保留初始解下越限最严重的Top-N场景求解后再迭代校验。实测下来这个方法能把求解时间缩短一个数量级以上。另外Yalmip建模时尽量避免在循环里使用string变量名拼接变量要用cell数组存储变量引用。MATLAB的字符串处理会拖慢循环性能这个在场景数上千时差距非常明显。6.3 光热电站SOC漂移问题前面提到的充放同时进行问题我实际跑的时候确实出现过模型为了满足N-k约束会在某些时段同时让储热罐充电和放电结果SOC始终保持在边界看着很好看但实际物理上做不到。解决方法是加入充放互斥的0-1约束。代价是引入整数变量求解变慢。如果系统规模很大可以做一个折中不要求全程严格互斥而是只对储热SOC接近上限或下限的时段做互斥约束其余时段放开效果也不错。另一个SOC相关的坑是SOC的初始值设置。模型里光热电站的SOC初始值定了末端SOC要不要强制回到初始值如果不约束模型会把储热当成免费资源结尾时段疯狂放热把SOC用完。实际运行中调度周期末端需要保留一定的储热。我的做法是加上末端SOC不低于初始SOC的80%这样的约束具体比例可以根据你的实际场景调整。6.4 PTDF符号和单位错误PTDF矩阵的正负号取决于你如何定义潮流正方向和节点注入正方向。我习惯上规定注入为正、负荷为负线路潮流正方向是支路首端到末端。如果你从Matpower拿到的数据里bus矩阵的PD列是负荷值构建注入向量时要记得取负号。单位方面Matpower的makePTDF默认用标幺值如果你把出力单位设成MW但没有除以baseMVA潮流结果会被放大或者缩小100倍这类错误非常隐蔽查起来最耗时间。6.5 关于后续扩展的一点体会模型跑通、结果分析做完之后再往后走有几个很自然的扩展方向。比如把确定性风电出力改成场景树加入风功率预测误差的机会约束比如把单时段调度扩展为多时段机组组合增加爬坡约束和启停变量再比如在N-k安全约束基础上叠加考虑储能系统的快速响应支撑。我个人在实际调试这个模型时的最大体会是N-k安全约束的核心难点从来不在某个公式上而在于故障场景和调度变量之间那一大坨线性关系的组织。只要把PTDF计算、节点注入关联矩阵、故障集枚举这三件事理清楚剩下的就是把约束按数学表达式翻译成Yalmip代码。希望这篇文章能让你少走几个我走过的弯路真的把代码跑起来、把结果做出来我写这些的目的就达到了。