简介:该资源是一份用牛拉法(Newton-Raphson法)实现含风电场电力系统潮流计算的程序文档,面向电力系统方向的研究人员、工程师及电气类课程学习者,用于处理风电并网后输出随机、不确定条件下的电网电压与功率分布分析。资源包共1个docx文件,压缩包约40KB,内容以MATLAB程序代码与算法说明为主,包含导纳矩阵与雅可比矩阵构建、节点类型判定、迭代收敛判断等核心模块。程序通过节点数、支路数、平衡节点号与误差精度四个输入量组织计算,并在节点类型为4时按风电场参数(有功功率、定子电抗、转子漏抗、转子电阻、励磁电抗)单独处理其电压与电流特性,区分平衡节点、PQ节点、PV节点与风电场节点四类求解逻辑。已有105人学习浏览,适合希望理解牛拉法迭代流程、雅可比矩阵装配与高斯消去求解,并据此开展风电并网潮流仿真与电网性能优化的读者参考。
1. 风电场接入后,潮流计算为什么不能再直接套用常规程序
一个 100 MW 风电场通过 220 kV 升压站并入电网,有人图省事,直接把-100 MW当作恒功率负荷灌进常规潮流计算程序,算出来并网点电压 0.94 p.u.,而现场录波是 0.99 p.u.。差的这 5% 不是程序写错了,是风电场本身不是"只出有功"的恒功率源:恒速异步机组在发出有功的同时必须从系统吸收无功建立磁场,双馈机组在变流器容量内可以支撑机端电压,直驱机组无功范围更宽。三种机型在潮流方程里的处理方式完全不同,节点类型也会随控制策略切换。这个标题真正要解决的,是把风速、风机等值电路和节点类型一起塞进牛顿-拉夫逊或者 PQ 分解法的迭代框架,让含风电场的电网算得出、算得准,而不是把风电场简化成一条负的负荷线。适合手上已经有一份潮流计算程序、正准备接入风电场的电气工程师,也适合做新能源并网仿真的研究生。
2. 风电场在潮流方程里怎么摆:节点类型选择与三种风机的建模路径
2.1 从 PQ、PV、平衡节点说起,风电场卡在哪一步
常规电力系统潮流计算把每个节点归成三类:PQ 节点给定有功、无功,待求电压幅值与相角;PV 节点给定有功和电压幅值,待求无功与相角;平衡节点给定幅值与相角,承担全网功率差额。负荷节点是 PQ,火电、水电机组母线是 PV,这种分工几十年没变过。
风电场卡在这一步:有功由风速决定,本身就是一个随时间波动的量;无功取决于机端电压和控制策略,既不能像负荷那样事先给定,也不能像同步机那样在一个宽范围内自由调节。如果按"负的负荷"处理,风速一变、机端电压一变,无功误差会直接传导到整个区域的电压分布上,潮流计算的收敛性和准确性一起垮掉。
工程上比较可靠的做法,是把风电场节点在 PQ 与 PV 之间按控制方式动态切换,对恒速机组则引入一个把转差率作为内层变量的 RX 模型。节点类型不是拍脑袋定的,而是由机组类型和变流器控制方式推出来的。
2.2 恒速恒频异步风机:RX 模型与转差率内层迭代
固定转速的鼠笼式异步发电机没有励磁控制,必须从电网吸收无功建立磁场,所以它的无功 Q 不是给定值,而是机端电压 U 和输出有功 P 的函数。给定风速下的机械功率之后,机端无功与电压耦合,必须用内层迭代解出来。
等值电路的结构是:定子电阻 R1 与定子漏抗 X1 串联,后面接励磁电抗 Xm 与转子支路(R2/s + jX2)的并联,其中 s 是转差率。给定机端相电压 Uph,输入阻抗写成
Zin = R1 + jX1 + (jXm)(R2/s + jX2) / (jXm + R2/s + jX2)电流 I = Uph / Zin,三相视在功率 S = 3 · Uph · conj(I)。输入有功 Pe(s) 在正常转差范围内随 s 单调上升,所以求满足 Pe(s) = Pg 的 s 用二分法最省事:下限取 1e-5,上限取 0.3,迭代五六十次精度足够。求出 s 之后,Q = imag(S) 就是这台机组从系统吸收的无功,作为该风电节点的无功负荷进入潮流方程。
注意:二分的单调区间只在转差率小于临界转差率时成立。如果风速给的有功超过机组最大出力,二分法会一路把 s 推到上界,这时候要先把有功限幅,否则每次潮流迭代都会得到同一个错误无功。
三个参数最敏感:Xm 决定励磁无功的量级,R2 决定转差率对有功的灵敏度,R1+X1 影响机端电压跌落。这三个值如果拿不到厂家数据,可以从机组铭牌的额定功率、额定电压、额定功率因数和最大转矩倍数反推,误差能控制在可接受范围。
2.3 双馈与直驱风机:PV 节点和恒功率因数控制的差别
双馈机组通过转子侧变流器实现有功无功解耦,无功调节能力受转子电流限制;直驱永磁机组通过全功率变流器并网,无功范围更宽,通常按电网导则要求参与电压调节。潮流计算里对应两种处理:
- 恒功率因数控制模式:按 PQ 节点处理,Q = P · tan(arccos φ)。风速变化引起 P 变化时 Q 要跟着重算,建议在每轮潮流迭代前更新一次风机出力,而不是整段仿真用一组固定值。
- 恒电压控制模式:按 PV 节点处理,Q 待求但有上下限。每轮迭代后检查 Q 是否越限,越限就把该节点转成 PQ,Q 固定在限值上再继续迭代,这就是无功越限检查,和同步机 PV-PQ 转换是同一套逻辑。
风电场关口一般还配有 SVG 或者并联电容器组,这部分补偿容量在潮流计算里并进等值注入,不要漏掉,否则并网点电压会偏低零点几到一 p.u. 的量级。
2.4 风电场内部等值与并网点:箱变和集电线路怎么处理
多台风机经箱变升压,再通过集电线路汇到主变接入并网点。如果只关心并网点电压水平,可以把整个风电场等值为一台机组加一条等值阻抗,计算量下降一个量级;如果要评估集电线路损耗或者机端电压分布,就得保留每台机或每一回集电线路作为独立节点。
| 机型 | 有功给定方式 | 无功特性 | 潮流节点类型 | 内层变量 |
|---|---|---|---|---|
| 恒速异步(鼠笼) | 风速查功率曲线 | 与电压、转差率强耦合,吸收无功 | PQ + RX 内层迭代 | 转差率 s |
| 双馈(DFIG) | 风速查功率曲线 | 变流器可调,常按恒功率因数 | PQ,越限后固定 Q | 无功越限标志 |
| 直驱永磁(PMSG) | 风速查功率曲线 | 全功率变流器,无功范围宽 | PV 或 PQ | Qmin / Qmax |
选择哪种处理方式,取决于研究目的和能拿到的数据。算区域电网电压分布,等值成一台机组够了;算风电场内部损耗和机端过电压,必须保留详细模型。
3. 用 MATLAB 搭一套含风电场的潮流计算程序
3.1 数据组织:节点表、支路表和风机参数表
写程序之前先把数据表定死,后面改参数只动数据不动算法。节点表每行至少包含编号、节点类型、基准电压、有功负荷、无功负荷、电压初值、相角初值;支路表每行包含首末端、电阻、电抗、对地电纳和变比;风机参数单独一张表,存单机容量、功率曲线关键点和等值电路参数。
%% data_prep.m —— 组装节点表与支路表 % 节点表 bus: [编号 类型 基准电压kV Pd Qd V0 th0] % 类型: 1=PQ, 2=PV, 3=平衡节点 % 支路表 branch: [首端 末端 R X 对地电纳/2 变比] baseMVA = 100; % 系统基准容量,风电场参数必须换算到同一基准 bus = [ 1 3 220 0 0 1.00 0 % 平衡节点 2 1 220 80 30 1.00 0 % 地区负荷 3 1 35 20 8 1.00 0 % 低压侧负荷 4 2 220 0 0 1.01 0]; % 风电场并网点,先按 PV 处理 branch = [ 1 2 0.010 0.080 0.020 1.00 1 4 0.015 0.100 0.025 1.00 3 2 0.020 0.120 0.015 0.95];节点 4 先按 PV 给定电压 1.01 p.u.,如果算完发现所需无功超过风电场实际能提供的范围,再改成 PQ 并给定 Q,这是标准的无功越限处理流程。基准容量一定要统一,风机厂家给的等值电路参数往往是折算到机组自身容量下的标幺值,直接拿去用会导致整个雅可比矩阵比例失衡、迭代发散。
3.2 节点导纳矩阵的生成
function Y = build_ybus(nb, br) % nb: 节点数; br: 支路表 [from to R X B/2 k] Y = zeros(nb); for t = 1:size(br,1) i = br(t,1); j = br(t,2); y = 1/(br(t,3) + 1i*br(t,4)); % 串联导纳 b = 1i*br(t,5); % 线路对地电纳 k = br(t,6); % 变比,1 表示无变压器 Y(i,i) = Y(i,i) + y/(k*k) + b; % 非标准变比按 π 型等值折算 Y(j,j) = Y(j,j) + y + b; Y(i,j) = Y(i,j) - y/k; Y(j,i) = Y(j,i) - y/k; endbuild_ybus的逻辑很简单:每条支路把自己的串联导纳和对地电纳累加到对应的对角和非对角元素上。变比不为 1 时,首端一侧的自导纳要除以 k 的平方,互导纳除以 k。这段代码的常见坑是把对地电纳两边都按 B/2 加了却忘了变压器励磁支路,或者把变比方向搞反,导致并网点电压偏高。
3.3 牛顿-拉夫逊主循环与数值雅可比
含风电场之后,风电节点的无功是电压的函数,解析雅可比要多推好几项,写错一项就发散。用数值差分求雅可比是更稳的选择:对每个待求变量加一个小扰动,观察不平衡量的变化率。
function [V, th, iter] = nr_pf(Y, Pspec, Qspec, type, V0, th0, tol, maxit) % type: 1=PQ, 2=PV, 3=平衡节点 idxP = find(type == 1 | type == 2); % 有功方程:PQ 与 PV 节点 idxQ = find(type == 1); % 无功方程:仅 PQ 节点 V = V0(:); th = th0(:); for iter = 1:maxit mis = mismatch(Y, Pspec, Qspec, idxP, idxQ, th, V); if max(abs(mis)) < tol, break; end h = 1e-6; nx = numel(idxP) + numel(idxQ); J = zeros(nx, nx); for k = 1:numel(idxP) % 对相角求偏导 t2 = th; t2(idxP(k)) = t2(idxP(k)) + h; J(:,k) = (mismatch(Y,Pspec,Qspec,idxP,idxQ,t2,V) - mis) / h; end for k = 1:numel(idxQ) % 对电压幅值求偏导 v2 = V; v2(idxQ(k)) = v2(idxQ(k)) + h; J(:,numel(idxP)+k) = (mismatch(Y,Pspec,Qspec,idxP,idxQ,th,v2) - mis) / h; end dx = J \ mis; th(idxP) = th(idxP) + dx(1:numel(idxP)); V(idxQ) = V(idxQ) + dx(numel(idxP)+1:end); end end function mis = mismatch(Y, Pspec, Qspec, idxP, idxQ, th, V) Vc = V .* exp(1i*th); S = Vc .* conj(Y * Vc); mis = [Pspec(idxP) - real(S(idxP)); Qspec(idxQ) - imag(S(idxQ))]; endtol一般取 1e-8 到 1e-6,maxit取 20 到 30。数值雅可比的好处是加进风机模型后不用重推公式,代价是每轮要多算2n次不平衡量,节点数上千时速度会明显慢下来,这时候再换解析雅可比。差分步长h取 1e-6 附近比较合适,取太大会把非线性误差引进来,取太小会被浮点精度吃掉。
3.4 风电场节点的无功修正嵌进每轮迭代
关键一步是在每轮迭代前,用当前机端电压更新风电机组的无功注入,让功率不平衡量和雅可比都反映最新的电压水平。
% 每轮潮流迭代前更新风电场无功注入 for w = 1:nw U = abs(V(busW(w))) * Vbase; % 机端线电压,kV [Qw, s] = induction_Q(U, Pw(w), mach); % RX 内层迭代求无功 Qspec(busW(w)) = -Qw / baseMVA; % 吸收无功记为负注入 endinduction_Q内部用二分法求转差率,返回机组吸收的无功。如果风电场按恒功率因数控制,这一段直接换成Qspec = -Pw*tan(acos(pf))/baseMVA即可,程序结构不用动。注意 Q 的符号约定要和节点表一致,同一套程序里混用两种符号约定是排查半天找不到原因的经典事故。
3.5 结果输出与并网点指标
收敛之后建议输出三组量:各节点电压幅值和相角、各支路首末端功率和损耗、风电场的实际无功出力与功率因数。并网点电压是判断风电场是否满足并网要求的核心指标,网损率用来验证结果是否合理,风电场实际功率因数和预设值对比可以验证控制模式有没有实现正确。第一版程序跑通一个 demo 算例就够,先确认潮流收敛、结果能对上手算,再往里面加规模和时序。
4. 风速、功率曲线与多机等值的参数设定
4.1 从风速到风机有功:Weibull 采样与功率曲线插值
风电场出力来自风速,风速本身是随机变量,常用两参数 Weibull 分布描述。做时序潮流或者概率潮流时,先用逆变换法生成风速序列,再查功率曲线得到有功。
% 风速采样与功率曲线查表 k = 2.0; c = 9.0; % Weibull 形状参数与尺度参数 u = rand(8760,1); v = c * (-log(1-u)).^(1/k); % 逆变换采样,生成全年逐时风速 % 2 MW 机组功率曲线关键点:风速(m/s) — 出力(MW) pc = [0 0; 3 0; 4 0.10; 6 0.60; 9 1.50; 12 2.00; 25 2.00]; P = interp1(pc(:,1), pc(:,2), v, 'linear', 0); P(v < 3 | v > 25) = 0; % 切出段置零k影响风速分布的形状,内陆风场通常取 1.8 到 2.2;c是尺度参数,量级和年平均风速接近,9 m/s 对应年平均风速约 8 m/s 的地区。interp1的第四个参数0表示超出表格范围时外插值取 0,防止高风速段出现负出力这种荒唐结果。P(v < 3 | v > 25) = 0这行显式处理切入和切出风速,别指望插值函数自动帮你做这件事。
4.2 多台风机的等值聚合:容量加权与等值风速
一个 50 MW 风电场可能装了 25 台 2 MW 机组,逐台建模会让节点数膨胀。常见等值方法有三种,选哪种看你要算什么。
| 等值方法 | 适用场景 | 关键参数 | 误差来源 |
|---|---|---|---|
| 单机倍乘 | 同型号、布置集中 | 台数 N | 尾流导致的出力差异 |
| 容量加权等值 | 不同型号混装 | 容量权重系数 | 参数分散性 |
| 分群等值 | 风速差异明显的大风场 | 分群数与质心风速 | 分群边界选取 |
等值阻抗常按容量加权求倒数和:Zeq = 1 / sum(1./Zi)。如果所有机组型号相同,直接倍乘最简单,等值机容量取全场容量,等值阻抗取单机阻抗除以台数,等值风速取各机风速的容量加权平均。尾流效应明显时,下游机组风速要按尾流模型衰减,否则算出来的风电场出力会系统性偏高。
4.3 一份能直接抄的参数表
| 参数 | 推荐取值 | 说明 |
|---|---|---|
| 基准容量 | 100 MVA | 与全网一致,风机标幺值必须换算 |
| 收敛精度 tol | 1e-6 p.u. | 小于 1e-8 收益有限且易触发浮点噪声 |
| 最大迭代次数 | 20~30 | 超过 30 次还没收敛基本是模型或初值问题 |
| 转差率上限 | 0.3 | 超过说明有功给大了,要先限幅 |
| RX 内层迭代次数 | 60 | 二分法,60 次精度远超潮流需求 |
| 数值雅可比步长 | 1e-6 | 太大引入非线性误差,太小被精度吃掉 |
| 风速 Weibull k | 1.8~2.2 | 内陆取小值,沿海取大值 |
| 风速 Weibull c | 年平均风速 × 1.1 | 粗略估计,有条件用实测数据拟合 |
| 切入 / 额定 / 切出风速 | 3 / 12 / 25 m/s | 按机型手册填,不要凭经验猜 |
参数最容易被忽视的是基准容量。风机等值电路参数通常以机组容量为基准,如果不去折算直接塞进 100 MVA 基准的网络,导纳矩阵元素会差几十倍,潮流注定发散,而且报错信息只会告诉你矩阵奇异,不会告诉你基准不一致。
5. 不收敛、程序报错与结果校核的排查技巧
5.1 潮流发散先查这四处
风机无功模型接进来之后发散,九成问题出在这几个地方。第一,节点类型冲突:把风电场同时设成 PV 又给它固定 Q,或者平衡节点选在了风电场上。第二,初值太差:重载线路用平启动,第一轮电压修正量过大直接跳出收敛域,改成用上一次断面结果热启动通常就好。第三,RX 内层迭代不收敛:有功超过机组最大出力,转差率顶到上界,无功不再变化但也不正确,先检查功率曲线和机组容量是否匹配。第四,Q 的符号约定:吸收无功记为负注入还是正负荷,程序里必须统一,混用会让不平衡量差出两倍。
5.2 MATLAB 常见报错与"无法定位程序输入点"
Matrix is singular to working precision一般意味着雅可比矩阵奇异,常见原因是某个孤岛节点没有和主网连通,或者 PV 节点转 PQ 之后方程数量对不上。Index exceeds matrix dimensions大多是节点编号和数组下标相差 1,MATLAB 从 1 开始编号,而很多数据文件从 0 开始。如果把程序用 MATLAB Compiler 打包成独立 exe 发给现场运行,可能出现"无法定位程序输入点",这是打包版本对应的运行时库和现场安装版本不匹配导致的,重新分发对应版本的运行库即可,和潮流算法本身无关。调试时把每轮迭代的最大不平衡量打印出来,收敛过程应该是单调下降,出现振荡或者突然上升,基本能定位到具体节点。
5.3 用功率因数和网损率做结果校核
算完不要只看程序有没有报错,至少做三项校核。并网点实际功率因数应该落在风电场设计范围内,恒功率因数控制模式下算出来的功率因数如果偏离设定值超过 0.02,说明无功注入没更新。全网网损率拿来做横向对比,含风电场的算例网损率通常比无风电场时略高或略低,但不会出现数量级差异,差太多要检查等值阻抗和变压器变比。最后把平衡节点的出力差额和风电场出力相加,和总负荷加网损对一下,功率平衡差在 1e-6 p.u. 量级才算收敛到了正确解,而不是收敛到了一个数学上成立但物理上不合理的点。
本文还有配套的精品资源,点击获取