CN108123434B - 一种计算pv曲线斜率以求取pv曲线运行点的方法 - Google Patents

Abstract

本发明公开了电力系统仿真计算领域中的一种通过计算PV曲线斜率以求取PV曲线对应运行点的方法。包括：输入电网参数，形成节点导纳矩阵；对节点进行分类，设定状态变量初值，并得到原始的不平衡方程；确定负荷增长方式，推导节点PV曲线斜率的表达式，并用该式对原有方程进行替换，为节点斜率设定初值，得到新方程；运用数学方法求解该方程，其中具体数学方法包括牛顿法，高斯法等。本发明提供的一种通过计算PV曲线斜率以求取PV曲线对应运行点的方法，能够对PV曲线任意点的位置精确求解，并且具有收敛性好，精度高，运算时间短，契合电力系统在线计算需求的需求。

Description

一种计算PV曲线斜率以求取PV曲线运行点的方法

技术领域

本发明属于电力系统静态电压稳定性分析领域，尤其涉及一种计算PV曲线斜率以求取PV曲线运行点的方法。

背景技术

研究电力系统的静态电压稳定性的主要内容是通过计算得到节点PV曲线，进而计算节点电压稳定裕度，得到电压崩溃点的数据。调度员以此为根据执行相应的调度指令，从而确保电力系统稳定运行。

然而，在节点电压达到其静态电压稳定极限(电压崩溃点)时处，普通牛顿-拉夫逊法潮流计算中用到的雅克比矩阵奇异，潮流无收敛解，系统运行点难以求取。

目前，针对这一问题，主要有两种解决方法，分别是连续潮流算法以及直接法。连续潮流算法的特点在于通过引入负荷参数，使得扩展的雅克比矩阵在电压崩溃点处非奇异，存在收敛的潮流解，进而得以画出完整的PV曲线。然而，这种方法需要反复迭代求解多次潮流，且迭代速度受到步长选择的影响，计算量大且耗时长，并不契合电力系统。直接法的特点是通过扩展非线性系统，使得扩展后的雅克比矩阵同样在电压崩溃点处非奇异，存在收敛潮流解，从而可以快速计算电压崩溃点。但是，这种方法只能得到电压崩溃点一个运行点的数据，而不能得到整条PV曲线的信息。另一方面，它的收敛性也受到初值的影响，不合适的初值往往会导致算法无法收敛。现在电力系统多采用上述两种方法分析系统的静态电压稳定性，但这两种方法存在的缺陷应该得到重视。

发明内容

针对现有的连续潮流法和直接法存在的缺陷，提出了一种新的针对电力系统静态电压稳定性分析的计算方法来克服已有方法的缺陷。

本发明通过计算PV曲线斜率以求取PV曲线运行点，不仅能够得到PV曲线任意点的具体信息，而且计算速度快，收敛性好，契合电力系统在线计算的需要。

一种计算PV曲线斜率以求取PV曲线运行点的方法，包括以下步骤：

步骤A：根据系统参数建模并建立目标方程，包括形成节点导纳矩阵，将系统节点类型分为PQ节点，PV节点以及平衡节点三类，分别形成不平衡量方程，并对各节点变量赋初值。其中，不平衡量方程为

式中P_i Q_i分别是发电机向i节点注入的有功功率和无功功率，U_i为节点i的电压模值G_ij B_ij为节点导纳矩阵相应位置元素的实部和虚部，θ_ij为节点i与节点j的相角差。

步骤B：设置负荷增长方式，找到负荷增加的PQ节点i，并定义该节点的斜率不平衡方程dP_i/dU_i＝C,式中C取值范围为-20-0，改变C的值即可得到PV曲线的对应点。其中，dP/dU的具体解析表达式求解过程如下，涉及变量均采用极坐标形式表达。这里假设求取第i节点的PV曲线。

步骤B1：将不平衡量方程方程式(4)展开，并在展开式两边同乘U_i并同除dU_i，得到式(5)

式中k＝2...i-1,i+1...n,l＝n-m-1...n；

m为PV节点个数，符号

代表求偏导数,d代表求全导数。不平衡方程组按此顺序排列：PQ节点有功方程，PV节点有功方程，PQ节点无功方程，节点i无功方程。将式(5)写成矩阵形式

式中J'矩阵是由J_kl'子矩阵按照前述排列方式所对应形成的方阵，且

步骤B2：针对式(5)，式子左端功率均为给定值，视作常数，不随任意状态变量改变而改变，根据求导定义，式子左端为零。此时联立求解式(5)，可得到斜率不平衡量表达式(6)

其中

利用式(6)替换掉该节点的有功功率不平衡方程，并与式(4)联立，得到扩展后的非线性方程，记为

f(U,θ)＝0 (7)

步骤C：利用牛顿法或高斯迭代法求解方程(7)，求出修正量[ΔU,Δθ]，并根据修正量最大值判断是否满足max(|[ΔU,Δθ]|)＜ε，式中ε＝0.000001，满足则输出结果,不满足则返回步骤C继续迭代计算

其中步骤C中所述利用牛顿法求解方程，包括以下步骤：

步骤C1：采用极坐标系计算，根据节点类型，为各节点变量赋初值。其中对针PQ节点电压模值与相角初值分别赋为1与0，针对PV节点，只对其相角赋初值为0，针对所求节点i的PV曲线斜率，赋初值为0，并根据公式计算不平衡量初值

步骤C2：根据已有的方程(3)表达式，形成雅克比矩阵J,代入数据求解，得到修正量[ΔU,Δθ]。

步骤C3：根据修正量最大值判断是否满足max(|[ΔU,Δθ]|)＜ε，式中ε＝0.000001，满足则输出结果,得到系统在该负荷增长条件下的状态，不满足则执行步骤C4

步骤C4：根据修正量迭代，得到各节点变量新的计算值，并根据公式代入数据计算不平衡量，执行步骤C2

步骤C中所述利用高斯迭代法求解方程，包括以下步骤：

步骤C1：根据已有非线性方程组(1)，将其整理为x＝g(x)的形式，式中x代表各节点变量，包括电压模值，相角以及第i节点的PV曲线斜率

步骤C2：采用极坐标系计算，根据节点类型，为各节点变量赋初值。其中对针PQ节点电压模值与相角初值分别赋为1与0，针对PV节点，只对其相角赋初值为0，针对所求节点i的PV曲线斜率，赋初值为0，将初值向量记为x₀

步骤C3：按照方程x_k+1＝g(x_k)迭代计算，式中k代表第k次迭代，并计算每次迭代后修正量的最大值，其中修正量计算公式为Δx＝x_k+1-x_k

步骤C4：根据修正量最大值判断是否满足max(|Δx|)＜ε，式中ε＝0.000001，满足则输出结果,得到系统在第i节点斜率为下的状态，不满足则执行步骤C3继续迭代计算

本发明具有以下优点：

(1)：与连续潮流法相比，计算电压崩溃点时达到的精度更高，且由于不存在步长选择这一环节，计算速度变快，收敛性更好。

(2)：与直接法相比，不存在初值选择这一环节，收敛性明显强于直接法。另一方面，本方法可以通过简单地更改参数来精确求取PV曲线上任意点的位置，而直接法只能计算电压崩溃点的位置。

(3)：本方法具有高效、鲁棒性强的特点，适合应用于大型电网在线计算，且提供了新的判别静态电压稳定性的指标，扩展性强。

附图说明

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

图1是本发明提供的一种通过计算PV曲线斜率以求取PV曲线对应运行点的流程图；

图2为IEEE118节点系统图。

图3为118节点算例结果图。

图4为算法效率对比图。

具体实施方式

本算例采用IEEE118节点系统进行计算，下面结合该系统对本方法的实施进行进一步说明。下面结合附图和实施例对本发明做进一步详细说明。

步骤A：输入电网参数，包括线路，发电机，变压器以及相关负荷参数，并对变压器支路做π形等值变换，形成节点导纳矩阵。对118节点系统的节点类型进行梳理：已知注入有功功率和电压的节点为PV节点，已知注入有功功率和无功功率的节点为PQ节点，已知电压幅值和相角的节点为平衡节点。考虑含118节点电力系统，包含1个平衡节点，m个PV节点。将电网负荷参数化，定义不平衡量方程

式中P_i Q_i分别是发电机向i节点注入的有功功率和无功功率，U_i为节点i的电压模值G_ij B_ij为节点导纳矩阵相应位置元素的实部和虚部，θ_ij为节点i与节点j的相角差。

步骤B：设置负荷增长方式，找到负荷增加的PQ节点i，并定义该节点的斜率不平衡方程dP_i/dU_i＝C,式中C取值范围为-20-0，改变C的值即可得到PV曲线的对应点。这里为求取电压崩溃点，设C＝0。本算例分别研究第17，20，29,48节点的静态电压稳定性，即i＝17，20,29,48.负荷增长方式设为分别只增加相应节点的有功负荷，不改变其无功负荷。收敛精度一律设为0.000001。

步骤B1：将不平衡量方程方程式(8)展开，并在展开式两边V_i同乘dV_i并同除，得到式(9)

式中k＝2...i-1,i+1...n,l＝n-m-1...n；

m为PV节点个数，符号

代表求偏导数,d代表求全导数。不平衡方程组按此顺序排列：PQ节点有功方程，PV节点有功方程，PQ节点无功方程，节点i无功方程。

将式(9)写成矩阵形式

式中J'矩阵是由J_kl'子矩阵按照前述排列方式所对应形成的方阵，且

步骤B2：针对式(9)，式子左端对应的节点功率均为给定值，视作常数，不随任意状态变量改变而改变，根据求导定义，式子左端为零。此时联立求解式(9)，可得到斜率不平衡量表达式(10)

其中

利用式(10)替换掉i节点的有功功率不平衡方程，并与式(8)联立，得到扩展后的非线性方程，记为

f(U,θ)＝0 (11)

步骤C：利用牛顿-拉夫逊法迭代求解方程(2)，求出修正量[ΔU,Δθ]，并根据修正量最大值判断是否满足max(|[ΔU,Δθ]|)＜ε，式中ε＝0.000001，满足则输出结果,不满足则返回步骤C继续迭代计算

以下是本方法的一个具体的实施方式，本方法通过Matlab程序实现，针对118节点系统，包括以下步骤：

1)输入相关电网参数，根据已有公式形成节点导纳矩阵

2)根据118节点系统节点类型，按照公式列写功率不平衡量方程

3)设置功率增长方式为只改变对应节点的有功功率，设置收敛精度为0.000001，通过数学推导得到特定PQ节点dP/dU的解析式，并将该式替换掉原来该节点的有功公共率不平衡方程，得到新的不平衡方程f(U,θ)＝0。其中斜率设定值C设为0，即对应着PV曲线的电压崩溃点。

4)运用数学方法求解新形成的不平衡方程，本算例采用牛顿-拉夫逊法求解，求出修正量[ΔU,Δθ]，并根据修正量最大值判断是否满足max(|[ΔU,Δθ]|)＜ε，式中ε＝0.000001，满足则输出结果,不满足则返回步骤4继续迭代计算。

现对算例结果进行分析。根据图2所示，本方法计算的第17,20,29,48节点的电压崩溃点电压模值分别为0.9387,0.6156,0.8222，0.6275；连续潮流法计算的电压模值则分别为0.9389,0.6156,0.8224,0.6242.以连续潮流结果作为参照，发现二者结果几乎一致，本算法准确性得到验证。另一方面，根据图4所示，与连续潮流法相比，在相同的前提条件下，本算法迭代次数更少，运行时间更短，更加符合电力系统在线计算的需求。

最后，应当指出：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员依然可以对本发明的具体实施方式进行修改或者等同替换，这些未脱离本发明精神和范围的任何修改或者等同替换，均在申请待批的权利要求保护范围之内。

Claims (1)

1.一种计算PV曲线斜率以求取PV曲线运行点的方法，其特征在于,包括以下步骤:

步骤A：根据系统参数建模并建立目标方程，包括形成节点导纳矩阵，将系统节点类型分为PQ节点，PV节点以及平衡节点三类，分别形成有功功率，无功功率不平衡量方程，并对各节点状态变量赋初值；

步骤B：根据系统负荷参数，确定系统负荷变化的基本方式，这里以全网负荷节点以同功率因数增长为例，系统潮流解将沿PV曲线变化；找到待研究的PQ节点i，并定义该节点的斜率不平衡方程为dP_i/dU_i＝C，式中C的物理意义为PV曲线对应运行点的斜率，取值范围为-20-20，改变C的值即可得到PV曲线的对应点，利用该方程替换掉该节点的有功功率不平衡方程，并与已有方程联立，得到扩展后的非线性方程组，记为

f(U,θ)＝0 (1)

所述步骤B中利用斜率方程替换有功功率不平衡量方程的过程包括：

步骤B1:确定负荷变化节点i后，将节点i对应的不平衡方程放置方程组最底端；不平衡量方程按公式展开，并在展开式两边同乘U_i并同除dU_i，得到式(2)

式中k＝2...i-1,i+1...n,l＝n-m-1...n；

P_i Q_i分别是发电机向i节点注入的有功功率和无功功率，U_i为节点i的电压模值，m为PV节点个数，n为网络节点总数；符号

代表求偏导数，d代表求全导数；不平衡方程组按此顺序排列：PQ节点有功方程，PV节点有功方程，PQ节点无功方程，节点i无功方程；将式(2)写成矩阵形式

式中J'矩阵是由J_kl'子矩阵按照前述排列方式所对应形成的方阵，且

步骤B2：以式(2)为基础，根据PV曲线斜率定义，可得出PV曲线斜率表达式：

针对式(2)，式子左端对应的节点功率均为给定值，视作常数，不随任意状态变量改变而改变，根据求导定义，式子左端为零；此时联立求解式(2)与斜率解析式，可得到斜率不平衡量表达式(3)

其中

用式(3)对该节点的有功功率方程进行替换；

步骤C：利用牛顿法或高斯迭代法求解方程(1)，求出修正量[ΔU,Δθ]，并根据修正量最大值判断是否满足max(|[ΔU,Δθ]|)＜ε，式中ε＝0.000001，满足则输出结果,不满足则返回步骤C继续迭代计算；

所述步骤C中利用牛顿法求解方程(1)，其具体步骤包括：

步骤C1：采用极坐标系计算，根据节点类型，为各节点变量赋初值；其中对针PQ节点电压模值与相角初值分别赋为1与0，针对PV节点，只对其相角赋初值为0，针对所求节点i的PV曲线斜率，赋初值为0，并根据公式计算不平衡量初值；

步骤C2：根据已有的方程(3)表达式，形成雅克比矩阵J,代入数据求解，得到修正量[ΔU,Δθ]；

步骤C3：根据修正量最大值判断是否满足max([|ΔU,Δθ]|)＜ε，式中ε＝0.000001，满足则输出结果,得到系统在该负荷增长条件下的状态，不满足则执行步骤C4；

步骤C4：根据修正量迭代，得到各节点变量新的计算值，并根据公式代入数据计算不平衡量，执行步骤C2；

所述步骤C中利用高斯迭代法求解方程(1)，其具体步骤包括：

步骤C1：根据已有非线性方程组(1)，将其整理为x＝g(x)的形式，式中x代表各节点变量，包括电压模值，相角以及第i节点的PV曲线斜率；

步骤C2：采用极坐标系计算，根据节点类型，为各节点变量赋初值，其中对针PQ节点电压模值与相角初值分别赋为1与0，针对PV节点，只对其相角赋初值为0，针对所求节点i的PV曲线斜率，赋初值为0，将初值向量记为x₀；

步骤C3：按照方程x_k+1＝g(x_k)迭代计算，式中k代表第k次迭代，并计算每次迭代后修正量的最大值，其中修正量计算公式为Δx＝x_k+1-x_k；

步骤C4：根据修正量最大值判断是否满足max(|Δx|)＜ε,式中ε＝0.000001，满足则输出结果,得到系统在第i节点斜率为C下的状态，不满足则执行步骤C3继续迭代计算。

Images