多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

风电场并网潮流计算:等值模型、MATPOWER与雅可比修正

风电场并网潮流计算:等值模型、MATPOWER与雅可比修正 简介这是一份面向电力系统研究人员、电气工程专业学生及并网工程技术人员的潮流计算程序文档聚焦含风电场接入场景下的电网稳态分析。资源以MATLAB代码讲解为主围绕牛拉法Newton-Raphson法求解非线性潮流方程组展开帮助读者理解风电场输出功率随风速波动的随机性对电压与功率分布的影响。压缩包仅含1个docx文件约40KB属于轻量级技术文档便于随取随读。内容从导纳矩阵的构建讲起说明支路参数矩阵与节点参数矩阵的数据组织方式并区分平衡节点、PQ节点、PV节点与风电场节点的不同处理逻辑尤其对节点类型为4时的风电场电压电流特性计算、雅可比矩阵元素推导及高斯消去法求解修正量等环节有较为完整的体现。该文档已有105人学习适合作为课程设计、毕业设计或科研入门阶段理解风电并网潮流建模与迭代求解思路的参考材料也可用于对照代码结构梳理算法流程与调试要点。1. 风电场接进来以后潮流算不准的根子在无功常规潮流里节点只有三种PQ 节点功率固定、PV 节点电压和有功固定、平衡节点兜底。风电场哪一种都不太像——它的有功随风速走更麻烦的是定速异步机吸收的无功与机端电压互相耦合电压掉下去它吸的无功反而变多节点无功不再是个常数。把风电场简单当成一个负的负荷去算并网点电压和实际值差两三个百分点是常事重载时连线路潮流方向都会判错。这篇把三件事讲透风电场怎么等值成潮流里的节点、在 MATLAB/MATPOWER 里怎么改数据、交替迭代不收敛时该动雅可比矩阵的哪一项。做新能源接入评审、风电场接入系统设计或者课程里要交一份含风电场的潮流计算程序都能直接照着跑。2. 含风电场的节点等值从风速序列到 P、Q 注入2.1 风速到有功功率曲线、尾流折减与场用电风电场送出的有功不是直接给定的得先从风速推出来。单机功率曲线是一条折线切入风速 3 m/s 以下不出力3 m/s 到 12 m/s 之间按风速立方关系爬升12 m/s 到 25 m/s 维持额定25 m/s 以上切出。风电场总出力不等于单机出力乘台数中间要乘两个系数——尾流折减系数通常取 0.900.95场用电率通常取 1%2%两个都乘完才是并网点注入。function P wind_power(v, vci, vr, vco, Pr) % v 当前风速 m/s % vci 切入风速vr 额定风速vco 切出风速单位均为 m/s % Pr 单机额定有功MW if v vci || v vco P 0; % 停机区间 elseif v vr P Pr * (v^3 - vci^3) / (vr^3 - vci^3); % 立方段近似风功率密度关系 else P Pr; % 额定段限功率运行 end end这段立方关系来自风功率密度与风速三次方成正比做规划阶段够用。但真正跑接入系统设计时我一般直接把厂家给的功率曲线表读进来用interp1(v_table, p_table, v, pchip)插值比套公式准也避免了不同机型曲线形状不一样带来的偏差。另外注意风速序列不要只取一个平均值——风速的威布尔分布在出力上不平均用平均风速算出来的出力会系统性偏大后面第五章会讲怎么按场景批量处理。2.2 风电场在潮流里的三种节点模型与选型把风电场折算成潮流节点工程上有四条常见路线选哪条决定了后面雅可比矩阵要不要动、程序好不好收。等值模型节点类型无功来源适用机组需要改雅可比吗恒功率因数 PQPQQ P·tanφφ 由并网点给定双馈、直驱自带无功控制不需要P-Q(V) / RXPQ由滑差方程解出 Q f(U)定速鼠笼异步机需要对角块加一项等值 PVPVSVG/STATCOM 维持电压配无功补偿装置的风电场需要加 Q 限值越限转 PQ恒阻抗 ZPQQ U²/X粗估网损、规划初筛不需要选型的判断很直接如果风机是双馈或直驱机侧变流器本身能调无功并网点功率因数可以做到 0.95 超前到 0.95 滞后用恒功率因数 PQ 模型误差在可接受范围这也是工程上八成的场景先用的模型。如果风电场装的是定速异步机或者要研究电压稳定裕度就必须上 P-Q(V) 模型因为无功—电压的负反馈特性正是电压失稳的诱因。等值 PV 模型只在风电场配了容量足够的无功补偿装置时才成立而且要在程序里加 Q 越限判断不然并网点电压会被强行定住算出来的无功需求超出补偿容量也不报警。等值参数计算里最容易写错的一步是阻抗折算。n 台同型机组并联等值成一台等值阻抗是被除不是被乘Z_eq Z_machine / n容量 S_eq n·S_machine。集电线路35 kV 电缆为主再按各自的长度和单位阻抗加到等值机与并网点之间。如果场内机型或风速分布差异大按风速分群做多机等值每群一台等值机加一段等值电缆通常 35 群就能把场内电压分布误差压到 1% 以内。2.3 定速机的 P-Q(V) 模型滑差方程怎么数值求解定速异步机的无功不能像 PQ 节点那样一次给定。它的稳态等值电路是定子阻抗 Rs jXs并联励磁支路 jXm再串转子支路 Rr/s jXr滑差 s 在发电机工况下为负值。给定机端电压 U 和有功出力 P滑差由有功平衡条件唯一确定无功随之算出。这个方程有解析解但形式很丑写成数值求解更省事、也更不容易在参数代入时出错。function Q q_induction(U, P, p) % U 机端电压(pu)P 有功出力(pu)p 为参数结构体 % 返回该机吸收的无功(pu)正号表示吸收 function F residual(s) Zin p.Rs 1j*p.Xs 1/(1/(1j*p.Xm) 1/(p.Rr/s 1j*p.Xr)); S U^2 * conj(1/Zin); % 从系统视角看机组的复功率 F abs(real(S)) - P; % 有功出力必须与给定值匹配 end s fzero(residual, [-0.05, -1e-6]); % 滑差为负发电机工况 Zin p.Rs 1j*p.Xs 1/(1/(1j*p.Xm) 1/(p.Rr/s 1j*p.Xr)); Q imag(U^2 * conj(1/Zin)); endfzero的搜索区间要卡在稳定运行段内左边界不能越过最大转矩对应的滑差一般在 0.030.05 之间越过去解出来的是不稳定工作点。有一组典型参数 Rs 0.005、Xs 0.10、Xm 2.5、Rr 0.008、Xr 0.10均为标幺值在 U 1.0、P 1.0 时解出 s ≈ -0.009Q ≈ 0.61对应并网点功率因数约 0.85。这个数字很关键定速机不带补偿时每发 1 MW 有功大约要从电网吸 0.6 Mvar 无功这部分无功得穿过集电电缆送进来电缆上的电压损耗会进一步压低机端电压形成正反馈。做风电场接入计算时看到功率因数低于 0.9基本就是这个原因。3. 用 MATPOWER 跑通含风电场潮流计算的最小流程3.1 算例准备与并网点选择先用现成算例把流程跑通比一上来就自己写牛顿法划算得多。以 MATPOWER 自带的case14为例14 个节点里 1、2、3、6、8 号是发电机节点其余是纯负荷节点。风电场要挂在 PQ 节点上选 9 号母线变压器低压侧或者 14 号母线末端都能看到电压变化这里用 9 号。选并网点时有个容易忽略的点不要把风电场挂在 PV 节点上。PV 节点的电压由同步机调压器撑着挂上去以后风电注入的无功引起的那点电压波动全被同步机吸收掉算出来的并网点电压几乎不动你以为模型是对的其实什么都没验证到。3.2 恒功率因数模型的写法负负荷注入MATPOWER 的bus矩阵里第 3 列是 Pd有功负荷MW第 4 列是 Qd无功负荷Mvar。风电场向外送有功、向内吸无功用负负荷表达最直接不用改gen矩阵也不用担心 Qmax/Qmin 的约束逻辑。%% wind_pf_pq.m —— 恒功率因数模型 mpc loadcase(case14); wf_bus 9; % 风电场并网点母线号 Pwf 60; % 风电场有功出力MW pf 0.98; % 并网点功率因数滞后风电场吸收无功 Qwf Pwf * tan(acos(pf)); % 折算成吸收的无功Mvar idx find(mpc.bus(:,1) wf_bus); mpc.bus(idx, 3) mpc.bus(idx, 3) - Pwf; % Pd 减去风电有功 mpc.bus(idx, 4) mpc.bus(idx, 4) Qwf; % Qd 加上风电吸收的无功 res runpf(mpc, mpoption(verbose, 0, out.all, 0)); fprintf(并网点 %d 电压 %.4f pu无功注入 %.2f Mvar\n, ... wf_bus, res.bus(idx, 8), Pwf*tan(acos(pf)));mpoption里verbose设 0 是为了批量跑的时候不刷屏单次调试时把它设成 2能看到每次迭代的节点电压和最大不平衡量排查收敛问题全靠这个输出。判断是否收敛不要只看命令窗口读res.success这个字段返回 1 才是收敛不收敛时res.bus里装的是最后一次迭代的结果直接拿去用会得到完全错误的结论。3.3 P-Q(V) 模型的交替迭代外循环恒功率因数是把无功钉死在初值上如果并网点电压偏离 1.0 比较多这个假设就崩了。P-Q(V) 模型的解法是交替迭代先用一个假设电压算出无功注入跑一次潮流得到新电压再用新电压重算无功如此往复。因为无功对电压的灵敏度有限这个外循环一般 48 次就稳定下来。%% wind_pf_pqv.m —— P-Q(V) 模型交替迭代 mpc0 loadcase(case14); wf_bus 9; Pwf 60; N 12; % 12 台等值机 p struct(Rs,0.005,Xs,0.10,Xm,2.5,Rr,0.008,Xr,0.10); Sb 5; % 单机容量 MVA U 1.0; tol 1e-5; for k 1:20 Qunit q_induction(U, Pwf/(N*Sb), p); % 单机吸收无功(pu) Qwf N * Qunit * Sb; % 风电场总无功(Mvar) mpc mpc0; idx find(mpc.bus(:,1) wf_bus); mpc.bus(idx,3) mpc.bus(idx,3) - Pwf; mpc.bus(idx,4) mpc.bus(idx,4) Qwf; r runpf(mpc, mpoption(verbose,0,out.all,0)); Unew r.bus(idx,8); fprintf(iter %2d U %.5f Qwf %8.2f Mvar\n, k, Unew, Qwf); if abs(Unew - U) tol, break; end U Unew; end这段代码里有三个参数决定了结果对不对。Sb是单机容量基准Pwf/(N*Sb)把风电场总有功折算成单机标幺值这一步搞错的话q_induction里解出来的滑差会跑到不稳定区。N是等值台数等值阻抗随 N 变化但在这个模型里体现为无功—电压曲线的斜率台数越多曲线越平。tol取 1e-5 是电压收敛判据比潮流本身 1e-8 的收敛精度松一个量级这样外循环不会因为内层残差噪声来回震荡。3.4 结果校验三个必看的量跑完以后不要只看并网点电压合不合格至少核对三个量。第一是并网点电压正常情况下应落在 0.971.07 之间偏低说明无功补偿不够。第二是风电场注入电网的无功把Qwf加上补偿装置出力得到的是并网点净无功这个值要和接入系统设计给的风电场功率因数范围对照。第三是网损变化把风电出力投运前后的系统总网损相减新能源接入后网损下降是常态但如果反而上升往往说明潮流反向穿越了主变这时候要检查风电出力是不是超过了本地负荷。4. 含风电场潮流不收敛时解法与雅可比怎么调4.1 快速解耦法在风电场面前为什么会失效很多人习惯用快速解耦法XB/BX 版本跑高压网潮流速度快、内存小。但一旦把风电场加进来这个方法经常报错或者迭代几十次不收敛原因在 R/X 比上。快速解耦法成立的前提是线路电抗远大于电阻高压架空线 R/X 通常在 0.1 以下这个假设没问题而风电场内部集电线路以 35 kV 电缆为主截面大、距离短R/X 可以到 1 甚至更高PQ 分解时那个有功主要受相角影响、无功主要受电压影响的解耦假设直接不成立。经验做法是主网部分继续用快速解耦没问题但只要算例里含风电场等值网络就换成完整牛顿—拉夫逊法。MATPOWER 里通过mpoption(pf.alg, NR)指定代价是每次迭代计算完整雅可比、耗时增加但收敛可靠得多。如果算例规模大又必须提速可以考虑把风电场内部的电缆网络先做等值合并把节点数压下来再跑。4.2 Q 随 V 变化时雅可比对角块要补哪一项这是含风电场潮流计算里最核心的一处代码改动。标准牛顿法的修正方程写成矩阵形式是[ ΔP ] [ H N ] [ Δθ ] [ ΔQ ] [ M L ] [ ΔV/V ]其中不平衡量定义为 ΔP P_spec − P_calc、ΔQ Q_spec − Q_calc。常规 PQ 节点的 Q_spec 是常数对电压求导为零所以 L 块只包含 −∂Q_calc/∂V 这一部分。但风电场节点的 Q_spec f(V) 是电压的函数线性化时必须多出一项L(i,i) −V_i · ∂Q_calc,i/∂V_i V_i · dQ_spec,i/dV_i也就是在原来 L 块的对角元上加上 V·dQ/dV 这一项。系数用中心差分算就行dV 1e-4; dQdV (q_induction(UdV, Ppu, p) - q_induction(U-dV, Ppu, p)) / (2*dV); L_diag_extra U * dQdV * (N*Sb) / baseMVA; % 折算到系统基准典型参数下这个导数在 1.52.5 之间乘上电压就是 2 左右的量级和 L 块原有的对角元同阶绝不是可以忽略的小量。漏掉它的后果很具体解本身还是对的因为每次迭代都会用更新后的电压重算 Q_spec不动点没变但收敛性质从二次退化成线性迭代次数从 5 次左右涨到 1520 次重载或者无功缺额大的场景直接发散。很多人调这个程序时反复怀疑是初值问题其实根子在这里。4.3 不收敛时的五个排查点排查顺序检查项典型现象处理方式1等值阻抗是否除了台数无功需求异常大Z_eq Z/n容量 S_eq n·S2滑差求解区间是否越过最大转矩点解出的 Q 是负值或数量级离谱收缩 fzero 区间到 [-0.04, -1e-6]3外循环是否震荡电压在两组值之间来回跳加阻尼U 0.5·U_old 0.5·U_new4是否用了快速解耦法报 maximum iteration exceeded换成pf.alg NR5是否有孤立节点或零阻抗支路雅可比奇异矩阵求逆失败检查集电电缆合并后的支路参数第 3 条值得展开说。交替迭代在无功—电压曲线斜率很陡的时候会震荡也就是 dQ/dV 特别大的工况比如风电场无功补偿容量刚好卡在临界点。这时候不要急着换算法先给外循环加个 0.5 的松弛因子通常两三次迭代就稳住了。如果加了松弛还震荡那说明物理上确实没有稳定解得回头看无功补偿容量配置。5. 多风速场景批量计算与并网点电压校验单个风速点算出来的结果没有多少说服力接入系统评审要的是全年或典型日的风速序列下的电压包络。批量的关键不是循环本身而是两件事场景怎么取、初值怎么给。场景不要按风速均匀取按风速的威布尔分布分区间取概率中值比如把 025 m/s 分成 10 段每段取该段的期望风速作为代表再按各段的年小时数加权统计结果。这样 10 个场景就能覆盖全年比均匀取 25 个点还有代表性。初值给法直接影响批量跑的时间。每个场景都从平启动电压全 1.0、相角全 0开始20 个场景的耗时是单场景的 20 倍把上一个场景收敛后的电压和相角直接作为下一个场景的初值因为相邻场景的风速接近、工况连续通常两次迭代就收敛整体耗时能降到三成左右。MATPOWER 里通过修改mpc.bus(:, 8)电压幅值和mpc.bus(:, 9)相角来实现热启动。%% wind_pf_scan.m —— 风速序列批量潮流与电压校验 v_set [5 7 9 11 13 15 18 22]; % 代表性风速 m/s mpc0 loadcase(case14); wf_bus 9; Pr 60; U_prev ones(size(mpc0.bus,1),1); % 热启动初值 A_prev zeros(size(mpc0.bus,1),1); tab zeros(numel(v_set), 3); for i 1:numel(v_set) Pwf wind_power(v_set(i), 3, 12, 25, 1.0) * Pr; mpc mpc0; idx find(mpc.bus(:,1) wf_bus); mpc.bus(:,8) U_prev; mpc.bus(:,9) A_prev; mpc.bus(idx,3) mpc.bus(idx,3) - Pwf; mpc.bus(idx,4) mpc.bus(idx,4) Pwf*tan(acos(0.98)); r runpf(mpc, mpoption(verbose,0,out.all,0,pf.alg,NR)); if ~r.success, fprintf(风速 %g m/s 不收敛\n, v_set(i)); continue; end U_prev r.bus(:,8); A_prev r.bus(:,9); tab(i,:) [v_set(i), Pwf, r.bus(idx,8)]; end disp(array2table(tab, VariableNames, {风速,出力MW,并网点电压pu}));跑完拿到这张表判据就三条并网点电压在所有场景下都落在 0.971.07 之间网损随风速的变化单调或接近单调出现非单调拐点说明潮流发生了反向风电场出力最大的那几个场景并网点功率因数不低于接入系统设计给定的下限。三条里任意一条越限回到第三章调无功补偿容量或者风电场等值参数重算不要在收敛判据上找补——那只是把问题藏起来。本文还有配套的精品资源点击获取
返回列表