计及光伏并网的电力系统频率特性和网络方程电压可解性分析

2023-11-11王雪一郝文清李振垚甘德强倪静怡

电力自动化设备 2023年10期

王雪一，郝文清，李振垚，甘德强，毛荀，倪静怡

（1.浙江大学电气工程学院，浙江杭州 310027；2.国网安徽省电力有限公司电力科学研究院，安徽合肥 230601）

0 引言

中国电网呈现大规模交直流混联状态，送端新能源发电机组大量并网，受端直流大规模馈入。送受端耦合紧密，直流故障将导致功率大幅波动，对电网的影响由局部转为全局。随着常规机组被大量替换，系统转动惯量降低，电网承受有功冲击的能力下降，频率稳定问题日益突出［1］。同时，新能源机组由于耐频、耐压能力不足，在系统频率、电压大幅变化的情况下容易脱网，引发连锁故障［2］。为保障电力系统安全稳定运行，在传统稳控方案的基础上，需要研究紧急功率控制策略，使新能源参与一次调频、调压［3］。新能源逆变器具备毫秒级快速响应指令的潜力［4］，若其能够应用于快速实现功率回降或提升（送端电网以功率回降为主），则将减少切机量和切负荷量。将系统频率进行反馈后，参与新能源输出功率控制，可以有效缓减电网频率调节压力，文献［5］指出光伏逆变器快速频率响应能力优于常规火电机组一次调频，具有良好的应用前景。

目前，针对新能源附加频率控制的研究大多聚焦于控制策略的制定，根据仿真结果评估方案的有效性［6］，量化评估方法较少。衡量系统频率动态响应的基本特征量主要包括惯量响应阶段的频率变化率、最大频率偏差和稳态频率偏差。随着新能源大量并网，系统频率特性发生改变，亟需开展适用于高比例新能源电力系统的频率稳定量化评估方法研究。文献［7］利用单调控制理论，对低频减载的作用规律进行了论证。文献［8］对系统频率的共模分量进行了刻画。文献［9］基于系统全状态模型得到状态变量时域响应的解析解，进而可分析得到系统参数和系统频率动态响应间的关系。

新能源机组缺乏电压支撑能力已是业界共识，文献［10］指出规模化光伏在电网故障后，其“低电压大电流”特性会造成电网无功缺额，恶化电压稳定性。但上述文献所采用方法均基于仿真计算，缺乏理论分析，难以揭示光伏接入对短期电压稳定的影响机理。多机电力系统本质上可以用一组微分代数方程式来描述，电网稳定机理得益于内在的系统特性［11］，网络方程的可解性为电压稳定分析提供了量化指标。文献［12］针对只存在一个非线性负载的系统，给出了网络方程可解性的明确条件，但该方法无法应用于含多非线性负载的场景。文献［13］推广了上述结果，进一步将最小奇异值作为衡量短期电压稳定的指标，可应用于含多非线性负载的场景，但是所提指标为复数形式，仍无法解释实虚部的物理意义。文献［14］解析推导了电压-无功灵敏度雅可比矩阵的结构特点，可以从网络方程自身特点出发分析光伏接入对电压的影响。

本文首先介绍光伏发电系统的数学模型，指出其典型响应特性。其次，通过求解电力系统状态矩阵，分析系统频率相关特征量与控制参数的单调关系，指出系统频率响应主要由主导模式决定，给出主导特征值的解析计算式。然后，聚焦网络方程，基于雅可比矩阵奇异性，提出一种负定性指标来量化电压可解性，并将该指标推广到含光伏的电力系统中，借助图论分析方法，论证光伏渗透率提高将会恶化电压可解性。最后，在4 机11 节点经典算例和实际华东电网中，仿真分析了光伏并网对电力系统频率特性和网络方程可解性的影响，验证了频率、电压分析方法的有效性。

1 光伏发电系统机电暂态模型

1.1 计及低电压穿越控制

光伏机组通过逆变器并网，可以实现功率的快速回降或提升，具有和同步机不同的响应特性。光伏发电系统的一般结构如图1 所示，主要包括有功功率控制模块、无功功率控制模块、故障穿越状态判断模块等［15］。图中：|Imax|为光伏注入电网电流的最大幅值。

图1 光伏发电系统一般结构示意图Fig.1 General structure schematic diagram of photovoltaic power generation system

光伏控制模块一般采用双环解耦控制，由于电流内环响应迅速，在机电暂态时间尺度下可以忽略内环控制过程［16］，外环一般采用定有功功率控制和定无功功率控制。当电网受到大扰动冲击时（如三相短路），为避免低电压下逆变器过流而导致脱网，光伏机组将进入低电压穿越（low voltage ride through，LVRT），LVRT 控制模型见文献［17］。合理设置LVRT控制参数对系统稳定有重要作用，附录A图A1 举例说明了LVRT 恢复期间，增大有功电流恢复速度可以减少光伏注入交流电网有功功率缺额Aloss，降低故障对频率的冲击。本文暂不考虑光伏LVRT 控制过程，但本文所提方法可以进行拓展，分析LVRT 控制策略对电力系统稳定的影响，这是后续研究将展开的课题。

1.2 计及附加频率控制

为改善频率特性，本文进一步考虑光伏的附加频率控制功能，当系统频率偏离设定值时，光伏将快速调节输出有功功率，参与调频，控制框图如图2 所示。图中：Pref为光伏机组输出有功功率的初始参考值；ΔPref为光伏附加频率控制所调节的有功功率；为计及Pref和ΔPref后更新的有功功率参考值；Ppv为光伏机组输出有功功率的实际值；Id为光伏机组输出有功电流的参考值／实际值；f为电网惯性中心（center of inertia，COI）频率；Tf为时间常数；fref为频率的参考值；kf为光伏附加频率控制增益常数；PI表示比例积分控制器。

图2 光伏附加频率控制框图Fig.2 Block diagram of photovoltaic supplementary frequency control

图3 给出了光伏定无功功率的控制框图。图中：Qref为光伏机组输出无功功率的参考值；Qpv为光伏机组输出无功功率的实际值；Iq为光伏机组输出无功电流的参考值／实际值。

正常运行状态下，光伏数学模型如下：

式中：kpp、kpi为有功功率控制环节的PI 控制器参数；kqp、kqi为无功功率控制环节的PI 控制器参数；xp、xq、xf分别为光伏机组有功状态变量、无功状态变量和频率状态变量；Vpv为光伏机组机端电压，|Vpv|为其幅值。

2 电网频率特性量化分析

2.1 状态空间模型时域解

下面将分析光伏紧急功率控制对系统频率特性的影响。同步机采用经典模型，考虑光伏附加频率控制后的系统可以近似用如下在平衡点处线性化后的五阶模型表示：

式中：ωN为同步机角频率的基准值；“Δ”表示对应物理量与其平衡点数值的偏差；xp、xq、xf、Pref、Vpv分别为由xp、xq、xf、Pref、Vpv构成的向量；δ、ω分别为同步机转子角和角频率向量；M、Mraw分别为由同步机惯性时间常数形成的对角矩阵和行向量；D为由同步机阻尼系数形成的对角矩阵；Pm为注入同步机机械功率向量；Pe为同步机输出电磁功率向量；In×1为元素均为1 的n维向量，n为光伏节点数目；Ml为第l台同步机的惯性时间常数；C1—C5、M′为系数矩阵，具体表达式见附录B式（B1）。

为消去代数变量ΔPe和Δ|Vpv|，对网络方程做如下线性化处理：

式中：Yee、Yepv、Ypve、Ypvpv为节点导纳矩阵Y（已将无源节点并入节点导纳矩阵）中元素，下标e、pv 分别表示同步机、光伏，Yee为同步机节点和同步机节点之间的节点导纳子矩阵，其他类似；S为视在功率向量；上横线表示取共轭；Ppv、Qpv分别为由Ppv、Qpv构成的向量；E=E′q.*ejδ，E′q为同步机q轴暂态电势向量，.*表示按元素相乘；θpv为光伏电压相角向量；./表示按元素相除；H1—H4、N1、N2、J1、J2、L1为对应的偏导矩阵。规定电流注入为正方向，将光伏输出功率表达式式（2）代入式（4），可以得到Δ|Vpv|、Δθpv的表达式如下：

式中：α1—α4、β1—β4为对应的偏导矩阵。在暂态过程中，光伏输出无功功率几乎不变，可令Δxq=0，又因为N1、N2、J1、J2的数值很小，故α1、α2、α4数值很小，式（3）可以进行化简，最终得到系统的线性化雅可比矩阵J为：

式（6）中各子矩阵块的表达式见附录B式（B2）。下面令：

令一台同步机为参考机可消去零特征值，此时雅可比矩阵为J′。研究系统的频率响应特性等同于求解线性系统Δx˙=J′Δx+BΔu的零状态响应时域解，当输入变量为阶跃函数时，系统解满足如下表达形式：

式中：λi和vi、qi分别为矩阵J′第i个特征值及其对应的右特征向量、左特征向量。

2.2 光伏并网对频率特征量的影响

2.2.1 最大频率变化率

COI 频率ωcoi( )t在平衡点处线性化后的表达式为：

式中：ωl为第l台同步机的角频率。

结合式（8），将式（9）对时间求导，得到COI频率变化率为：

式中：Pω为由右特征值向量构成的矩阵P中与角频率相关的部分；Eλ为由eλit形成的对角矩阵；Q为由左特征向量构成的矩阵。令t=0，得到惯量响应环节的频率变化率为：

式（11）结果表明，零时刻频率变化率只和系统同步机总惯量相关（成反比），光伏附加频率控制不改变零时刻频率变化率。

2.2.2 稳态频率偏差

系统状态变量稳态值Δx(∞)计算公式如下：

式中：Λ为由矩阵J′的特征值λi构成的对角矩阵。式（12）表明状态变量稳态值可以用雅可比矩阵的逆来表示：

根据分块矩阵求逆公式，令：

式中：Kf为由kf构成的对角矩阵；D～为考虑光伏附加频率控制后，系统的等效阻尼矩阵。经过整理，得到各台机组稳态转速偏差Δω(∞)表达式为：

式中：I为单位矩阵；Kpp为由kpp构成的对角矩阵；-D～-1+D～-1L～′(I′D～-1L～′)-1I′D～-1为一个行相等矩阵，故所有发电机组最后收敛至同一转速。

在系统等效阻尼系数表达式中，H2近似为一个元素全负的实数矩阵，H4近似为一个M 矩阵，其逆为非负矩阵，所以H2H-14元素符号为负，Kf的作用相当于增加了系统阻尼，即D～ ≥D，从而减小了受扰动后系统频率的稳态误差。

定义参与因子pji，其表达式为pji=vjiqiBΔu，vji为右特征向量vi的第j个元素，参与因子量化了特征值λi对状态变量Δxj的贡献。对于实特征值σ和复特征值σ±jω，设其参与因子分别为pr和pr±jpi，则时域解分量fji(t)如下：

不同的特征值代表了不同的频率模式，可以认为系统频率响应只和主导模式相关，将主导模式定义为由靠近原点的一对共轭特征值决定的频率响应轨迹，计算结果表明，该对共轭特征值的参与因子较大，附录B 图B1 给出了IEEE 10 机39 节点系统的特征值分布示意图。

下面对系统主导特征值进行求解，令JΔx=λΔx，λ为所求的特征值，展开后得到：

对式（19）进行整理后，得到：

当全网光伏控制参数一致时，有J25=-kfkppJ23，且由于Δω≠0，式（20）成立的条件如下：

式中：det 表示求矩阵的行列式。经过推导，得到主导特征值的计算公式如式（22）所示，推导过程和参数含义见附录C。

式中：下标d 表示主导模式。由一元三次方程求解式（22）知其有解析表达式，存在一对共轭主导特征值。进一步，根据式（17），可以得到频率稳态偏差理论计算值，即-Fd。

2.2.3 最大频率偏差

非主导模式对系统频率的贡献很小（具体可参见下文算例），在分析时不妨将其忽略，将主导模式频率fd对时间求导，可得：

令式（23）为0，解得：

式中：下标max 表示频率偏差最大时的值。将式（24）代入式（17），则得到最大频率偏差的解析解表达式为：

3 电网电压特性量化分析

3.1 网络方程电压可解性

电力系统机电暂态模型可以统一描述为一组微分-代数方程（differential-algebraic equation，DAE），一般情况下，代数方程的奇异面在DAE 模型的稳定边界上，这表明代数方程的可解性对电压稳定具有重要影响。下面研究网络方程的奇异性条件，推导反映电压可解性的指标。当同步机忽略凸极效应和阻尼绕组，负荷采用恒功率类型，光伏采用定有功功率和定无功功率控制方式时，将网络方程重写为：

式中：下标E、L分别表示同步机内电势节点、恒功率节点（含光伏和恒功率负荷），YEE为同步机内电势节点和同步机内电势节点之间的节点导纳子矩阵，YEL、YLE、YLL含义类似。

由于E′q与δ为状态变量，求解网络方程就是在已知E和SL的基础求VL，则存在如下非线性方程：

采用牛顿法求解上述非线性方程，式（27）等号两侧同乘以恒功率节点电压对角矩阵的共轭diagVˉL，则有：

忽略线路电阻，即YLL≈jBLL，BLL为节点导纳子矩阵YLL的虚部，即节点电纳子矩阵，并且令IE=YLEE，得到网络方程的近似结果为：

为表述简洁，下文将省略下标L、LL。将式（29）在直角坐标系下展开，并求得从元件端口看进去电力网络侧的雅可比矩阵为：

式中：Pg、Qg分别为有功功率向量、无功功率向量；下标x、y分别表示对应变量的实部、虚部；偏导矩阵Jpx、Jpy、Jqx、Jqy的表达式见附录D式（D1）。

若计及光伏、负荷等元件的动态模型，则从电力网络侧看进去，元件的端口特性能够用类似的矩阵表示，即：

采用定功率控制模式的光伏逆变器在准稳态时间尺度下，外特性接近理想恒功率源。下文在分析网络方程可解性时，可不计光伏的出力波动，忽略雅可比矩阵J的第2 部分。由于矩阵J非对称，一般情况其特征值为复数，需要求解最小奇异值，为得到形式更简便的指标，下文对矩阵J做进一步的化简。

当系统电压可解时，其网络方程雅可比矩阵非奇异，等价于雅可比矩阵的行列式不为0，即detJ≠0。又因为雅可比矩阵的行列式等于其特征值的乘积，即detJ=λ1λ2×…×λJ（J为特征值数量），所以最小特征值λ1( )J≠0 是一个等价的可解性条件。此处的λ1( )J指雅可比矩阵中模值最小的特征值，它可以视作一种量化电压可解性的指标。

雅可比矩阵J的表达式可以重新整理成如下形式：

式中：IBV、Vxy、BB的表达式见附录D式（D2）。

由式（33）可以得到雅可比矩阵J的行列式表达式为：

式中：Γ=V-1xy IBV+BB。值得指出的是，detVxy和Γ均具有特殊性质，有助于开展理论推导，具体论证过程见附录D定理D1和定理D2。

考虑到电网安全稳定运行时，不会出现恒功率节点电压降至0 的情况，因此detVxy≠0，雅可比矩阵的奇异性取决于矩阵Γ的性质。首先将矩阵Γ写成分块矩阵的形式：

式中：b=diag(S./V2)y，g=diag(S./V2)x，(·)x、(·)y分别为取相应矩阵的实部和虚部。

一般情况下Γ为负定矩阵，其特征值全负。为方便分析，对矩阵Γ做如下变换［18］：

相合，矩阵Γ 的负定条件是B -b 和(B +b)-g (B -b)-1g均负定，这意味着λ1(Γ′)可以用来衡量电压稳定裕度，即：

通常电纳矩阵B的数量级远大于b与g，所以λ1(B-b)将比λ1[(B+b)-g(B-b)-1g]更靠近原点，只有当g足够大时，舒尔补的最小特征值λ1[(B+b)-g(B-b)-1g]才会快速增大。

矩阵B为经Kron 降阶方法消去无源节点后的电纳子矩阵，它的负矩阵是带自环的Laplacian 矩阵，满足正定性［19］，即意味着Kron 降阶电纳矩阵B满足负定性。由于b数量级远小于B，所以B-b特征值符号主要由B决定，一般为负号。

根据以上论述B-b满足负定性，其特征值符号明确，易反映出参数变化对特征值的影响，因此λ1(B-b)是一种合适的负定性指标，具有清晰的理论表达式，参数更加直观，可以用来量化机电暂态网络方程的电压可解性。

3.2 光伏渗透率对电压可解性的影响

假设提高电网光伏渗透率的方式是保持系统潮流不变，逐步将电网中同步机替换为光伏机组。本节将借助图论等数学工具，分析光伏渗透率变化对电网电压可解性的影响。

定义同步机节点集为α，则光伏与负荷节点的集合为同步机节点的补集β。随着光伏渗透率的增加，同步机逐台被替换为光伏机组，同步机节点数减少，而光伏节点数增加。不妨假设k个非空节点集有如下关系：

柯西交错定理［18］表明对实对称矩阵做加边处理后，特征值单调不减，所以B-b的最小模值特征值随光伏渗透率的增加单调不减：

式中：Bβ1为对应光伏节点β1和光伏节点之间的节点电纳子矩阵，其他类似。

根据上述分析可知，光伏渗透率的增加将恶化网络方程电压可解性，影响电压稳定性，这一结果将在算例分析中进行计算说明。

4 算例分析

4.1 光伏并网对电网频率特性的影响

考虑4 机11 节点经典算例，系统接线图见附录E 图E1，将3 号同步机替换为光伏机组后，系统稳态潮流数据和发电机数据见附录E 表E1 — E3。取光伏附加频率控制的增益常数kf= 50，系统特征值如表1 所示，其中第1 个特征值为零特征值（可以令一台同步机为转子角参考机而去除），第2 个和第3 个特征值为一对共轭主导特征值，代表频率主导模式。

表1 系统特征值Table 1 System eigenvalues

设置4 号同步机损失10 % 的有功出力，系统频率将跌落。观察COI 时域仿真实际频率、线性化模型解析频率和主导模式频率，绘制响应曲线如图4所示。由图可知，线性模型解析解和实际频率响应结果较为接近，而主导模式曲线几乎和线性模型频率曲线完全重合。故可以单独研究频率主导模式，对系统频率响应曲线做近似刻画。

图4 线性模型解析解和时域仿真实际频率响应曲线Fig.4 Analytical solution of linear model and actual frequency response curve of time-domain simulation

为验证光伏附加频率控制增益常数kf对系统频率响应的影响，分别取kf为0（相当于无附加频率控制）、5、50、500 这4 种情况，得到COI 频率及光伏输出有功功率变化曲线如图5 所示。图中，光伏输出有功功率为标幺值。

图5 不同kf取值下系统响应曲线Fig.5 System response curves under different values of kf

图5 上图表明光伏附加频率控制能够有效抑制频率持续下降，根据2.2 节的分析，系统惯量保持不变，频率在零时刻以相同的斜率下降，而最大频率偏差和稳态频率偏差受增益常数kf的影响，一般情况下kf越大，系统等效阻尼越大，频率偏差越小。在损失电源出力场景下，光伏为了保持系统频率稳定，将提高有功出力。

图5下图表明增益常数kf越大，光伏有功出力的调整速度越快，调整量越大，故频率控制效果越好。值得注意的是，增益常数kf取值有上限，文献［20］指出kf取值过大，可能导致系统不稳定。表2 列出了不同增益常数取值下COI稳态频率和最低频率的计算值，可以为频率稳定控制措施的制定提供依据。

表2 不同kf 取值下线性化系统频率计算值Table 2 Frequency calculation values of linearized system under different values of kf

4.2 光伏并网对电网电压特性的影响

4.2.1 功角差平面电压可解性

仍采用4.1节中的4机11节点算例，若忽略光伏输出功率变化，将其视为恒功率电源处理，则由于电力网络只有1 台光伏机组，求解网络方程相当于求解一个一元二次方程，电压具有解析表达式，根据方程判别式和0 的大小关系可以分析其可解性，具体求解过程见文献［12］。将1 号同步机作为转子角参考机，绘制得到网络方程雅可比矩阵最小奇异值在功角差平面的三维变化图，如图6 所示，同网络方程判别式进行比较。图中：δ21为2 号同步机和1 号同步机的功角差；δ41为4 号同步机和1 号同步机的功角差。

图6 雅可比矩阵最小奇异值和网络方程判别式在功角差平面的三维图Fig.6 Three-dimensional plot of minimum singular value of Jacobian matrix and network equation discriminant in power angle difference plane

观察图6 可知，雅可比矩阵最小奇异值和网络方程判别式在功角差平面具有一致的变化趋势，这表明根据雅可比矩阵奇异性等值条件推导得到的负定性指标可以反映网络方程的实际可解性。

4.2.2 光伏渗透率对电压可解性的影响

考虑华东电网某年的运行方式，算例系统共有8 117 个节点，其中发电机节点个数为575，负荷节点个数为3 243。按照“小电站优先替换”以及“电源替换前后潮流维持不变”的原则，对华东电网电源结构进行调整，测试不同光伏渗透率对电压可解性的影响。算例中同步机采用暂态电势恒定模型，光伏外环采用定有功功率和定无功功率控制，负荷为恒功率类型，忽略线路电阻。负定性指标λ1(B-b)与光伏渗透率之间的变化关系如图7所示。

图7 电压可解性随光伏渗透率的变化关系（华东电网）Fig.7 Relationship between voltage solvability and photovoltaic penetration（East China Power Grid）

由图7 可知，在华东电网算例中，负定性指标λ1(B-b)随着光伏渗透率的增加呈单调递增趋势。这说明在大电网中同样存在随着光伏渗透率增加，电压可解性恶化的问题。电力系统实际可接纳的光伏极限渗透率和新能源并网位置、网络拓扑、电源及负荷动态特性等多种因素有关，本文所构造的负定性指标可以对光伏模型进行拓展，开展多场景下的电压稳定性分析。

为了观察光伏渗透率对电网电压特性的影响，对于原4机11节点算例，增大负荷水平，并进一步将系统中1号同步机替换为光伏机组，设置4号同步机发生损失有功出力故障，对比暂态过程系统中不同光伏机组数目下母线3 所连光伏节点电压幅值（标幺值）与负定性指标λ1(B-b)的变化曲线，结果如图8所示。