去年做新能源电力系统调度项目时甲方给我提了一个特别具体的要求调度方案里的每个发电计划都必须扛得住任意两台设备同时跳闸也就是N-k安全校验。一开始我觉得这不难等真正把风电、光伏、光热三种电源塞进同一个优化模型再把N-k约束逐条加进去才发现模型规模和求解难度完全是两回事。这篇文章是我把“含风电-光伏-光热电站的N-k安全优化调度模型”从公式一路写成Matlab代码的完整记录包括三类电源的建模差异、预想事故集的筛选逻辑、YALMIP下的约束生成以及实测里踩过的性能坑希望对做电力系统调度方向的同学有点参考价值。1. 新能源高占比场景下N-k安全为什么成为刚需1.1 从单重故障到多重故障N-1准则不够用了传统电力系统运行调度最常用的安全准则是N-1任意一条线路、一台变压器或一台发电机退出运行后系统仍能维持稳定供电。这个准则在传统电源为主的电网里够用了因为同步机组惯量大、调节能力强单重故障发生后备用容量和AGC能在较短时间内把系统拉回安全区间。但新能源渗透率上来以后情况发生了变化。风电、光伏出力随天气波动系统等效惯量下降功频调节能力变弱。在这种运行条件下单重故障叠加新能源出力的快速波动很可能在短时间内演变成多重故障。比如大风天气下风电满发多条线路负载率本来就高此时再发生一条线路跳闸潮流会转移到相邻线路如果转移潮流超过相邻线路的热稳定极限就会引发连锁跳闸。N-k准则正是针对这类场景提出的系统在任意k个元件同时故障退出运行后依然能满足安全运行要求。N-1到N-k看似只是故障数量从1变成k背后复杂度是指数级上升的。以118节点系统为例N-1只需要检查一百多条线路的单重故障而N-2要检查的组合数是C(189,2)大概有17766种。如果直接把这近两万种场景的潮流约束全部塞进优化模型即便是商用求解器也未必能在可接受时间内给出结果。1.2 风光间歇性叠加故障场景风险被放大风电光伏的出力本质上是一个随机过程。风速变化、云层遮挡都会让实际出力偏离预测值。在做日前调度时我们拿到的是预测曲线但真实运行时的出力跟预测可能偏差很大。这就带来一个隐患调度方案是基于预测出力求解的如果实际风电、光伏出力比预测值低常规机组就要顶上去如果比预测值高多余功率又要在电网里找出路。N-k安全约束与风光不确定性放在一起会产生一个叠加效应。比如调度模型在某个中等风电出力场景下满足了N-1安全但实际运行中风速突增风电出力大幅上升某些线路的潮流已经逼近极限此时再叠加一个元件故障线路潮流越限几乎是必然的。所以处理高比例新能源系统的N-k安全不能只看基态安全还要把不确定性纳入调度模型这正是这个模型要做的事。1.3 光热电站的角色从新能源到可调度电源在这个模型里光热电站不是普通的“新能源电源”它是作为可调度电源出现的。光热电站由集热场、储热系统、发电系统三部分组成。太阳集热场收集热量把热量储存到储热罐里需要发电时再释放热量驱动汽轮机。储热系统的存在让光热电站的出力可以在时间上“平移”中午太阳能多的时候少发电多储热晚上光伏出力为零的时候放热发电。这意味着光热电站既包含随机性集热量取决于太阳辐射又具有可控性储热罐可以缓冲出力计划能够调节。在N-k安全调度中它是一个天然的灵活性资源。风电光伏出力波动的风险可以由光热电站的储热来对冲故障场景下缺失的功率也可以由光热电站的快速调节来补足。所以含光热的系统做N-k安全调度比单纯的风光系统要好解得多前提是光热的储热模型建得够细。2. 三类电站的出力建模抓住随机性与可调度性的本质差异2.1 风电风速概率分布与出力场景生成风电出力建模的第一步是描述风速。经典做法是用两参数Weibull分布拟合风速的历史数据然后通过风机的功率曲线把风速转换成出力。功率曲线一般分成三段切入风速以下出力为零风速在切入和额定之间时出力近似线性上升额定风速以上保持额定功率切出风速以上为了保护风机而停机。在调度模型里直接采用风速分布非常麻烦工程上更常见的是场景法。我习惯用历史风速数据做蒙特卡洛抽样生成大量风速序列再用同步回代削减法把场景压缩到10~20个每个场景附带一个概率。这样风电出力在模型里既不是单点预测也不是完全确定的负荷而是以场景集合的形式参与求解结果对风速波动具备一定的鲁棒性。如果不想引入场景也可以把风电出力设为一个可调区间做区间优化或者鲁棒优化。区间取值太小调度结果在真实运行时很可能扛不住N-k故障取得太大系统要预留过多备用经济性又变差。实际经验是先用预测误差的历史分位数来确定区间边界比如取90%置信水平的预测区间再结合故障场景校核能在一个可接受的计算代价内保证安全性。2.2 光伏光照强度模型与日间出力曲线光伏出力与光照强度直接相关光照强度一般用Beta分布描述。光伏出力可以简单表示为额定容量乘以光照强度转换系数再考虑温度修正。光伏的特点在于它的日间规律非常明显夜间出力为零正午前后达到峰值晴天和阴天的差异是造成不确定性的主要原因。在调度模型中光伏的随机性处理方式与风电类似可以用场景法也可以视为波动区间。不过相比风电光伏的时间相关性更强24小时的出力曲线形状相对固定这给调度带来一点便利光伏的爬坡集中在日出和日落时段只要在这两个时段给常规机组留足爬坡空间大部分情况下系统响应是跟得上的。另外要提醒一点模型里光伏的“出力上限”是可以主动弃光的。很多时候最优解并不是让光伏满发而是在某些时段通过弃光来避免潮流越限这样整体调度成本反而更低。如果你在建模型时把光伏出力设置成固定注入功率那N-k安全约束很可能直接不可行。2.3 光热集热-储热-发电三层能量流建模光热电站的建模和风电光伏有本质区别它是一套多能量流系统涉及太阳能到热能、热能到电能的两次转换。调度时不能只盯着发电功率还要跟踪储热罐里的热量状态。集热环节的输入是太阳直接辐射DNI集热量等于集热效率乘以集热面积再乘以DNI。储热环节的动态约束是核心每个时段的储热增量等于充热功率乘以充热效率减去放热功率再减去储热损失同时储热量不能超过储热罐容量充热和放热不能同时进行。发电环节则比较简单发电功率等于放热功率乘以发电效率并受汽轮机最大功率限制。我建议这样定义变量来避免单位混乱P_ch表示从集热场进入储热罐的热功率P_dis表示从储热罐放出进入发电系统的热功率储热量递推关系为E_tes(t) E_tes(t-1) η_store * P_ch(t) - P_dis(t)发电功率P_csp(t) η_pb * P_dis(t)。这样每个变量的物理意义都很清楚写代码时不容易搞混。这个模型里最关键的约束是储热罐的时序约束它把各时段耦合在一起。比如中午的太阳能不能全部上网可以先存进储热罐晚上再释放发电。很多初学的人容易漏掉“充放热不能同时进行”这个约束导致模型在某些时段出现同时充放热的无效解白白浪费储热资源。实际建模中可以用二元变量处理也可以采用工程简化近似——设置一个很小的惩罚项同时充放热在经济上不划算从而在最优解里自然规避。2.4 为什么说光热让调度模型“活”了起来把三类电源放在一起看风电光伏是不可控的随机电源光热是可控但受热量约束的调节电源。调度的本质就是在风电光伏出力波动的情况下利用常规机组和光热储热来保证功率平衡与N-k安全。光热的加入让目标函数的调整空间变大了。比如夜间负荷高峰期光伏为零如果没有光热常规机组必须承担全部压力有了光热白天储存的热量可以在夜间释放系统不再那么依赖备用机组。在N-k故障场景下某条线路跳闸后潮流重新分配光热电站可以快速改变出力来消除越限这种调节能力是风电光伏不具备的。我把三类电源的建模维度整理了一下方便对照理解电源类型不确定性来源可控性建模核心风电风速随机波动基本不可控只能弃风功率曲线 场景生成光伏光照强度随机波动基本不可控只能弃光Beta分布 日间出力曲线光热太阳辐射波动储热罐提供可调度性集热-储热-发电三层能量流3. N-k安全约束的建模从组合爆炸到可解问题3.1 预想事故集N-k场景怎么筛直接枚举所有N-k故障组合是不现实的。118节点系统189条线路N-2有17766种N-3直接到一百万以上就算生成约束不花时间求解器也会被撑爆。所以第一步必须做故障筛选。我做故障筛选的思路是两层。第一层是电气距离筛选先分析基态潮流把负载率高的线路、以及与重载线路电气距离近的线路挑出来构成候选故障集。第二层是故障影响分析对候选故障集中的每个组合计算故障后的潮流分布以“是否引起其他线路越限”为判据保留那些真正会造成越限的严重故障场景剔除不影响安全校验的场景。筛选参数可以按项目要求调整。如果要求偏保守就保留更多场景如果计算资源紧张就把筛选阈值调严一些。实际项目里我会先做一次“筛选-校验-再筛选”的迭代先用初始筛选集求解模型再用求解结果重新评估故障影响把漏掉的严重故障场景补进事故集最多迭代两三次就能收敛。3.2 故障态潮流与直流潮流的快速计算N-k安全约束的计算核心是故障态潮流。如果用交流潮流每个故障场景都要重新迭代求解模型根本没法用优化器直接处理。工程上几乎都用直流潮流近似在这个模型里故障态潮流可以用线路开断分布因子LODF快速计算。LODF的基本思想是线路k断开后线路l上潮流的变化量与线路k断开前的潮流成正比。用公式表示就是 f_l_after f_l_0 LODF_{l,k} * f_k_0其中 f_l_0 是基态潮流LODF_{l,k} 是开断分布因子由网络阻抗参数决定不随调度方案变化。计算一次LODF矩阵之后每个故障场景的潮流约束都是线性组合求解效率非常高。使用LODF时要注意两个细节。一是LODF计算依赖于基态潮流的参考方向线路潮流为负时开断转移方向也可能反转代码里必须处理好正负号。二是当故障涉及多条线路同时断开时不能用单条线路的LODF简单叠加需要根据实际并列断开的线路重新计算或者使用广义开断分布因子GDDF。我的代码里保留的是单条和多条两类接口按需调用。3.3 直接枚举还是分解求解我的取舍方案处理N-k约束有两种主流思路。第一种是把所有筛选后的故障场景约束全部显式写入优化模型用求解器一次性求解这种方式实现简单适用于事故集规模不大的情况。第二种是Benders分解主问题只考虑基态约束子问题专门做安全校验发现越限场景后生成割约束回传给主问题迭代求解。我的做法是两者结合。对于N-1和筛选后规模较小的N-2事故集直接显式写入约束利用YALMIP的批量约束生成能力代码很简洁。如果项目里要求覆盖更广的N-2甚至N-3场景主问题就直接求解不动了我会在外面套一层Benders循环主问题求解一次、安全校验一次、生成割约束一次迭代10~20轮就能收敛比一次性写完整个模型要快得多。这也是我建议你在实际项目中优先考虑的方向。4. Matlab代码结构从数据读入到结果输出的完整链路4.1 环境搭建与文件组织我用的是Matlab YALMIP CPLEXYALMIP负责建模CPLEX负责求解。YALMIP的安装很简单把工具箱路径加进Matlab就行关键是求解器接口要装对。CPLEX的Matlab接口需要单独安装对应版本的插件版本不匹配时YALMIP会一直报错“no solver found”。代码按模块拆分文件我习惯这样组织main.m主程序定义系统参数和调度时段data_case.m节点系统数据包含负荷、线路、机组参数build_scenarios.m生成风光出力场景build_contingency.m预想事故集生成与筛选build_model.m构建优化模型决策变量、目标函数、约束solve_model.m调用求解器并处理结果plot_results.m出力曲线和储热状态可视化小系统验证用IEEE 30节点就够了跑通逻辑后想体现模型价值建议换IEEE 118节点系统N-k安全约束的规模效应会更明显。4.2 决策变量与目标函数的代码化模型里的决策变量包括常规机组出力、风电场出力、光伏出力、光热电站出力、储热罐充放热功率和储热量。用YALMIP定义非常直接代码如下P_g sdpvar(n_G, T, full); % 常规机组出力 P_w sdpvar(n_W, T, full); % 风电场出力 P_pv sdpvar(n_PV, T, full); % 光伏出力 P_csp sdpvar(n_CSP, T, full); % 光热发电功率 E_tes sdpvar(n_CSP, T, full); % 储热量 P_ch sdpvar(n_CSP, T, full); % 充热功率 P_dis sdpvar(n_CSP, T, full); % 放热功率目标函数我一般取全时段发电成本最小加弃风弃光惩罚。常规机组用二次成本函数风电光伏和光热的运行成本设为定值弃风弃光加一个较大的惩罚系数保证模型在满足安全约束时尽量不要弃新能源。这里有个经验惩罚系数不是越大越好太大会让求解器在目标函数数值上失去灵敏度一般取常规机组最大边际成本的1.5~2倍就够。4.3 N-k安全约束生成模块的实现细节安全约束生成是整个代码里最需要小心的地方。基态潮流约束是每个时段每条线路都要满足传输极限Constraints [Constraints, f_base B_net * theta(:, t)]; Constraints [Constraints, -f_max f_base f_max];然后把筛选好的事故集传给约束生成函数对每个场景构造故障态潮流约束。这里的关键是用预先算好的LODF矩阵把故障态潮流表示成基态潮流的线性组合for k 1:n_cont f_after f_base LODF_cell{k} * f_base; Constraints [Constraints, -f_max f_after f_max]; end这里我踩过一个坑LODF矩阵的行列对应关系必须和线路编号严格一致有一次调换线路序号后结果完全乱掉排查了很久才发现是数据对齐问题。建议在生成LODF矩阵之后用两三条已知线路的故障转移值做一次单元测试确认矩阵正确再往下走。光热储热约束和爬坡约束也需要在这个模块里一并建立% 储热动态递推 for t 2:T Constraints [Constraints, E_tes(:,t) E_tes(:,t-1) eta_store * P_ch(:,t) - P_dis(:,t)]; end % 储热容量与充放热上下限 Constraints [Constraints, E_min E_tes E_max]; Constraints [Constraints, 0 P_ch P_ch_max]; Constraints [Constraints, 0 P_dis P_dis_max];4.4 求解后的校验与可视化求解完成后不要直接拿结果走人一定要做一次独立校验把调度方案代入故障态潮流计算程序逐条检查所有故障场景是否越限。这一步虽然和模型里的约束是同一套逻辑但用独立脚本复核一遍能发现建模阶段漏掉的问题比如某些故障场景压根没被筛选进来。可视化方面我习惯画三张图常规机组和光热电站的出力堆叠图、风电场和光伏电站的出力与弃电情况、光热储热罐的储热量曲线。储热曲线特别能说明问题——如果储热量曲线频繁出现锯齿状震荡说明充放热约束或目标函数里可能有问题正常调度下储热量应该是平滑的充放循环。5. 实测性能与调试经验那些文档里不会写的坑5.1 模型规模急剧膨胀后怎么控制求解时间以IEEE 118节点系统、24个时段、20个风电场景、40个N-2故障场景为例模型的变量和约束数量会突破10万量级CPLEX求解时间可能从几分钟涨到几十分钟。我实测下来最有效的加速手段有三个。第一是给MIP设置合适的相对gap比如0.01而不是默认的0.0001求解时间能下降一半以上调度方案的工程误差完全可以接受。第二是合理给出变量初值用YALMIP的assign和initial_guess功能把基态调度的解作为热启动值二次求解和迭代求解时速度提升非常明显。第三是检查冗余约束比如某些故障场景约束在所有时段都非活跃可以直接删除我用了一段检测脚本自动剥离非活跃约束后模型规模能缩减两到三成。5.2 故障场景下模型不可行的排查链路调度模型在加入某个N-k场景后报infeasible是调试阶段最常遇到的问题。我总结了一套固定排查思路按顺序做能省很多时间。排查步骤操作可能结论第一步去掉所有N-k约束只求解基态调度基态不可行则是机组/功率平衡/储热约束问题第二步逐个加入单条线路故障约束找出引发不可行的具体故障线路第三步分析该线路故障后的潮流转移定位越限线路检查附近可调资源第四步检查风光场景区间与储热初始状态可能是不确定性范围过窄或储热初值不合理还有一个容易忽略的点是光热储热的初始状态。模型第一个时段的储热量如果初始值设得太小而前几个时段辐照又不足夜间需要放热时储热罐可能已经空了整个调度就不可行。我通常把初始储热量设为容量的一半以上并在模型末尾增加一个储热量终端约束要求结束时的储热量不低于初始值这样模型不会把储热罐“用完”。5.3 风光出力场景的粒度怎么选风光场景的粒度直接影响N-k安全结果的可信度。场景太少模型会漏掉那些对安全威胁最大的极端出力场景太多求解时间又会失控。我的经验是先用历史数据统计风电光伏出力的极端区间比如5%概率的大风低光照场景和95%概率的小风场景把这些极端场景强制加入场景集再配上10~15个常规场景。如果你做的是日前调度风光的预测窗口是24小时我建议把预测误差按4~6小时分段处理这样场景之间的相关性比简单整段随机抽样更合理。分段后每个时段内的不确定性更均匀N-k校核的结论也更稳定。5.4 把光热模型移植到电池储能、抽水蓄能的思路这套Matlab代码的可扩展性挺好。光热电站的储热模型本质上是一个储能模型把E_tes、P_ch、P_dis这三个变量改一下参数就能变成电池储能模型储能容量改成电池容量充放热效率改成充放电效率汽轮机约束改成逆变器功率上限即可。抽水蓄能也是类似区别是上下库容和抽水、发电效率要分别处理。如果要做N-1-1相继故障第一个故障后还没来得及恢复又发生第二个故障LODF部分需要改成顺序开断的迭代计算约束生成逻辑会复杂一些但整体框架不用大改。扩展的时候注意一点储能类约束的时序耦合特性决定了模型会越来越难解遇到求解时间失控时优先检查是否可以用场景削减或Benders分解来缓解不要盲目堆硬件资源。最后说一句实在话。N-k安全调度模型这东西理论框架就那么多真正的壁垒在细节数据对不对齐、场景筛得准不准、求解器参数调没调好任何一个环节都会让模型不可用。我的建议是先从IEEE 30节点小系统跑通再一步步往大系统上扩展把每个模块的输出打印出来肉眼检查一遍比任何复杂的算法技巧都管用。这套Matlab代码作为起点你完全可以在此基础上改出自己的版本遇到具体问题也欢迎一起交流讨论。