连续潮流的理论与编程

最新推荐文章于 2024-07-02 14:10:02 发布

圆周率的距离

最新推荐文章于 2024-07-02 14:10:02 发布

阅读量7.7k

点赞数 29

分类专栏：电力系统基本理论

本文链接：https://blog.csdn.net/qq_38642223/article/details/107959138

版权

电力系统基本理论专栏收录该内容

3 篇文章 31 订阅

订阅专栏

Index

连续潮流的理论简介及其编程
6. 程序链接

连续潮流的理论简介及其编程

1. 问题背景与连续潮流的基础

1.1 问题背景

根据我国《电力系统安全稳定导则》，电力系统的电压稳定性，是指电力系统在遭受小的或大的干扰后，系统电压能够保持或恢复到允许的范围内，不发生电压崩溃的能力。类似于功角稳定性的分析，通常将电压稳定性划分为小干扰电压稳定性和大干扰电压稳定性。而电力系统的电压稳定性问题，从本质上仍然归属于动态稳定问题，其数学模型与研究电力系统暂态稳定的数学模型基本相同，可以用如下的代数方程描述： $\Bigg\{ \begin{alignedat}{2} \frac {d\bm{x}} {dt} = \bm{f}(\bm{x}, \bm{V}) \\ \bm{I}(\bm{x}, \bm{V}) = \bm{Y} \bm{V} \end{alignedat} \tag {1}$

其中， $\bm{x}$ 是系统的状态变量向量， $\bm{V}$ 是系统的节点电压向量， $\bm{I}$ 是节点的注入电流向量， $\bm{Y}$ 则是系统的导纳矩阵。
与对暂态过程的研究类似，可以利用数值积分方法对式(1)进行求解，就可以对实际的电力系统进行电压稳定性分析，虽然这种方法的准确性较高，然而其计算时间较长，实用性较差。与之相对的，通过求取PV曲线（或QV曲线）来判断电力系统电压稳定性的静态分析方法则较为成熟。其中，连续潮流(Continuation Power Flow)，有文献译作延拓潮流，由于可以作出完整的PV曲线，因而在实际系统中得到较多应用。

1.2 连续潮流的基本原理与分类

在连续潮流法出现之前，电力行业通常通过大量的潮流计算（采用常规潮流的程序，如牛顿-拉夫逊法等）得出指定负荷节点的PV或QV曲线的方式，对电力系统的电压稳定性进行静态分析。这样的方法存在着一些不足：首先，费时；其次，需要谨慎地选择所观测的母线；最后，由于在临界点（PV曲线的“鼻点”）及其附近，常规潮流的Jacobi矩阵趋于奇异，导致迭代难以收敛，因此借助常规潮流是难以得到完整的PV曲线的。
而基于连续方法（Continuation Method，又称延拓法），是跟踪非线性代数方程组随参数变化的解曲线的基本方法。连续潮流使用连续方法跟踪负荷及发电机功率变化下的电力系统稳态行为，这就需要将常规潮流的Jacobi矩阵增加一阶，得到的增广Jacobi矩阵在临界点（PV曲线的“鼻点”）及其附近不再是奇异的，可以使迭代过程顺利穿过临界点，进而可以作出完整的PV曲线。目前而言，连续潮流的模型分类大致有四种：负荷型、支路型、故障型、控制型。其中又以负荷型的研究较为广泛，也是本文主要讨论的重点，下文中的“连续潮流”统一指代负荷型连续潮流。

2. 连续潮流的数学模型：引入参数的潮流方程

如前所述，PV曲线或QV曲线实际上反映的是负荷及发电机功率变化下的电力系统稳态行为，我们往往使用参数 $\lambda$ 的变化来表征这些变化。那么，对于一个含有 $n$ 个节点的网络（设其包含 $n_1$ 个PQ节点， $n_2$ 个PV节点）而言，系统参数化后的潮流方程可以简写为：
$\bm{f}(\bm{x}, {\lambda}) = 0 \tag {2}$

$\leq {\lambda} \leq {\lambda}_{critical} \tag {3}$

$\Bigg\{\begin{alignedat}{2} P_{Li}({\lambda}) = P_{Li0} + {\lambda}(k_{Li}S_{{\Delta}base}cos{\varphi}_{i}) \\ Q_{Li}({\lambda}) = Q_{Li0} + {\lambda}(k_{Li}S_{{\Delta}base}sin{\varphi}_{i}) \end{alignedat} \tag {4}$

$P_{Gi} = P_{Gi0}(1 + {\lambda}k_{Gi}) \tag {5}$

为了简化分析，假设所有负荷节点都保持初始状态下各自的功率因数，然后按照负荷各自的功率因数同步增加所有负荷，则式(4)可以简化为：
$\Bigg\{\begin{alignedat}{2} P_{Li}({\lambda}) = P_{Li0}(1 + {\lambda}k_{Li}) \\ Q_{Li}({\lambda}) = Q_{Li0}(1 + {\lambda}k_{Li}) \end{alignedat} \tag {6}$

常规潮流方程 ${\bm f}({\bm x}) = {\bm 0}$ 在极坐标下可以表示为：
$\begin{alignedat}{2} P_{Gi} - P_{Li} - U_{i} \small{\sum_{\mathclap{ i \isin j}}} [U_{j}(G_{ij}cos\theta_{ij} + B_{ij}sin\theta_{ij})] = 0 \\ Q_{Gi} - Q_{Li} - U_{i} \small{\sum_{\mathclap{ i \isin j}}} [U_{j}(G_{ij}sin\theta_{ij} - B_{ij}cos\theta_{ij})] = 0 \end{alignedat} \tag {7}$

将式(5)、(6)带入(7)中，可以得到负荷型连续潮流的数学模型 ${\bm f}({\bm x}, {\lambda}) = {\bm 0}$ 在极坐标下的表达式：
$\begin{alignedat}{2} P_{Gi0}(1 + {\lambda}k_{Gi}) - P_{Li0}(1 + {\lambda}k_{Li}) - U_{i} \small{\sum_{\mathclap{ i \isin j}}} [U_{j}(G_{ij}cos\theta_{ij} + B_{ij}sin\theta_{ij})] = 0 \\ Q_{Gi0} - Q_{Li0}(1 + {\lambda}k_{Li}) - U_{i} \small{\sum_{\mathclap{ i \isin j}}} [U_{j}(G_{ij}sin\theta_{ij} - B_{ij}cos\theta_{ij})] = 0 \end{alignedat} \tag {8}$

3. 连续潮流的主要步骤及程序流程图

假设我们研究一个含有 $n$ 个节点的网络，其中包含 $n_1$ 个PQ节点， $n_2$ 个PV节点。文末则提供了相关代码的下载链接。

3.1 预估(Predictor)

令 ${\bm x} =[{\bm \theta}, {\bm V}]^T$ ，则连续潮流的数学模型 ${\bm f}({\bm x}, {\lambda}) = {\bm 0}$ 可以写作 ${\bm f}({\bm \theta}, {\bm V}, {\lambda}) = {\bm 0}$ ，在取 ${\lambda}^{(0)} = 0$ （相当于常规潮流计算）后，首先可以得到初始解向量 $[{\bm \theta}^{(0)}, {\bm V}^{(0)}, {\lambda}^{(0)}]^T$ ，其中 ${\bm V}^{(0)}$ 为所有节点的电压幅值向量， ${\bm \theta}^{(0)}$ 为所有节点的电压相角向量。需要注意的是， ${\bm f}$ 与 ${\bm x}$ 的维数均为 $2n_1+n_2)$ ，然而连续潮流的未知量个数却有 $2n_1+n_2+1)$ 个（在 ${\bm x}$ 的基础上增加了一个未知量 ${\lambda}$ ），即未知量的个数多于方程的个数，这也是为什么要采用连续法求解这一类问题。
预估的目的是找出下一个解的近似值，并将该近似值作为校正步骤的初值。假设在连续过程中的第 $i$ 步，已知 ${\bm f}({\bm x}, {\lambda}) = {\bm 0}$ 的解 $[{\bm \theta}^{(i)}, {\bm V}^{(i)}, {\lambda}^{(i)}]^T$ ，那么预估步骤就是寻找第 $(i + 1)$ 步的解 $[{\bm \theta}^{(i+1)}, {\bm V}^{(i+1)}, {\lambda}^{(i+1)}]^T$ 的近似解（校正步骤初值） $[{\tilde {\bm \theta}}^{(i+1)}, {\tilde {\bm V}}^{(i+1)}, {\tilde \lambda}^{(i+1)}]^T$ 。采用的预估方法为：切线法。由于本文采用的参数化方法为局部参数化，因此接下来的切线法仅是针对局部参数化而言的。
对 ${\bm f}({\bm x}, {\lambda}) = 0$ ，即对式(8)的两边，对 ${\bm \theta}$ 、 ${\bm V}$ 、 ${\lambda}$ 取微分可得：
${\bm f}} = {\frac {\partial {\bm f}} {\partial {\bm \theta}}} {d {\bm \theta}} + {\frac {\partial {\bm f}} {\partial {\bm V}}} {d {\bm V}} + {\frac {\partial {\bm f}} {\partial {\lambda}}} {d {\lambda}} = {\bm 0} \tag {9a}$

写成矩阵形式为：
$\begin{bmatrix} \Large {\frac {\partial {\bm f}} {\partial {\bm \theta}}} & \Large {\frac {\partial {\bm f}} {\partial {\bm V}}} & \Large {\frac {\partial {\bm f}} {\partial {\lambda}}} \end{bmatrix} \begin{bmatrix} {d {\bm \theta}} \\ {d {\bm V}} \\{d {\lambda}} \end{bmatrix} = \bm{0} \tag {9b}$

式(9b)中，矩阵 $[{\Large {\frac {\partial {\bm f}} {\partial {\bm \theta}}}}, {\Large {\frac {\partial {\bm f}}{\partial {\bm V}}}}]$ 就是常规潮流计算的Jacobi矩阵，维数为 $(2n_1+n_2) {\times} (2n_1+n_2)$ ；而列向量 $[{\Large {\frac {\partial {\bm f}} {\partial {\lambda}}}}]$ 则是式(8)对 ${\lambda}$ 求偏导，其维数显然为 $(2n_1+n_2) {\times} 1$ 。因此， $[{\Large {\frac {\partial {\bm f}} {\partial {\bm \theta}}}}, {\Large {\frac {\partial {\bm f}}{\partial {\bm V}}}}, {\Large {\frac {\partial {\bm f}} {\partial {\lambda}}}}]$ 的维数为 $(2n_1+n_2) {\times} (2n_1+n_2+1)$ ，切向量 ${\bm \theta}}, {d {\bm V}}, {d {\lambda}}]^T$ 的维数为 $(2n_1+n_2+1) {\times} 1$ 。
显然，方程的个数比未知量的个数少1，这就需要添加一个方程。具体的解决方法是：指定切向量中的某一个分量为 $+ 1$ 还是 $- 1$ ，并称指定的这个分量为连续参数。在计算的最开始，选择 $\lambda$ 作为连续参数，因此切向量中的相应分量为 $+ 1$ ，而在后续的计算中，则需要根据3.3节所述的方式选择连续参数，其斜率的符号（连续参数增大/减小的趋势）决定切向量中相应分量的符号。通过引入连续参数，式(9b)变为：
$\begin{bmatrix} \Large {\frac {\partial {\bm f}} {\partial {\bm \theta}}} & \Large {\frac {\partial {\bm f}} {\partial {\bm V}}} & \Large {\frac {\partial {\bm f}} {\partial {\lambda}}} \\ {\bm{e}_k} \end{bmatrix} \begin{bmatrix} {d {\bm \theta}} \\ {d {\bm V}} \\ {d {\lambda}} \end{bmatrix} = \begin{bmatrix} \bm{0} \\ {\pm1} \end{bmatrix} \tag {10}$

式(10)中， ${\bm {e}_k}$ 为一行向量，维数为 $2n_1+n_2+1)$ ，从而使方程组维数与未知量个数匹配。其中，除了连续参数所对应的元素（设为第 $k$ 个元素，其索引为 $k$ ）为 $1$ 外，其余所有元素均为 $0$ 。
根据式(10)，容易解出切向量 ${\bm \theta}}, {d {\bm V}}, {d {\lambda}}]^T$ ，当采用局部参数化方式时，可以根据切向量及式(11)进行解的预估：
$\begin{bmatrix} {{\tilde {\bm \theta}}^{(i+1)}} \\ {{\tilde {\bm V}}^{(i+1)}} \\ {{\tilde \lambda}}^{(i+1)} \\ \end{bmatrix} = \begin{bmatrix} {{{\bm \theta}}^{(i)}} \\ {{{\bm V}}^{(i)}} \\ {{\lambda}}^{(i)} \\ \end{bmatrix} + {\sigma} \begin{bmatrix} {d {\bm \theta}} \\ {d {\bm V}} \\ {d {\lambda}} \\ \end{bmatrix} \tag {11}$

式(11)中， $[{\tilde {\bm \theta}}^{(i+1)}, {\tilde {\bm V}}^{(i+1)}, {\tilde \lambda}^{(i+1)}]^T$ 为PV曲线上第 $(i + 1)$ 个潮流解的预估值； $[{\bm \theta}^{(i)}, {\bm V}^{(i)}, {\lambda}^{(i)}]^T$ 为第 $i$ 个潮流解； ${\bm \theta}}, {d {\bm V}}, {d {\lambda}}]^T$ 为式(10)解出的切向量； ${\sigma}$ 为步长，合适的步长可以使预估值落在校正环节的收敛域内（离真实解足够近）。例如，假设在获取初始解 $[{\bm \theta}^{(0)}, {\bm V}^{(0)}, {\lambda}^{(0)}]^T$ 之后，结合由式(10)解出的切向量，选择合适的步长 ${\sigma}$ ，就可以沿着切线方向，得到PV曲线上的下一个潮流解 $[{\bm \theta}^{(1)}, {\bm V}^{(1)}, {\lambda}^{(1)}]^T$ 的预估值 $[{\tilde {\bm \theta}}^{(1)}, {\tilde {\bm V}}^{(1)}, {\tilde \lambda}^{(1)}]^T$ 。

3.2 校正(Corrector)

在根据式(11)获取真实解的近似解之后，就可以将其作为式(8)的初始解，进行解的校正，通常采用的方式是牛顿法和拟牛顿法等。然而式(8)对应的方程组 ${\bm f} ({\bm \theta}, {\bm V}, {\lambda}) = {\bm 0}$ 的个数仍然比未知量的个数少1，因此需要增加一个方程。
对于局部参数化，可以通过固定连续参数的方式获得附加方程：
$x_k - {\tilde x}_k = 0 \tag {12}$

式(12)中， $x_k$ 为连续参数， $k$ 为连续参数对应的索引， ${\tilde x}_k$ 为式(11)计算出的连续参数的预估值。
将式(12)与式(8)联立，可以得到增广潮流方程：
$\begin{bmatrix} {{\bm f}({\bm \theta}, {\bm V}, {\lambda})} \\ {x_k - {\tilde x_k}} \end{bmatrix} = {\bm 0} \tag {13}$

此时由于连续参数已经确定，即 $k$ 与$$已知，因此对于式(13)而言，只需对基于牛顿-拉夫逊法的常规潮流计算程序作些许修正即可用于对式(13)的求解，具体而言，仅仅只需将常规潮流计算程序的Jacobi矩阵及不平衡向量各自增加一阶即可。对增广Jacobi矩阵而言，在常规潮流的Jacobi矩阵基础上，增加一个列向量 $[\Large {\frac {\partial {\bm f}} {\partial {\lambda}}}]$ ，维数为 $(2n_1+n_2) {\times} 1$ ；增加一个行向量 ${\bm e}_k$ ，维数为 $2n_1+n_2+1)$ ，显然，除了第 $k$ 个元素为1外， ${\bm e}_k$ 的其他元素均为0，故而增广Jacobi矩阵的维数为 $(2n_1+n_2+1) {\times} (2n_1+n_2+1)$ 。在临界点附近，常规潮流的Jacobi矩阵接近奇异，但是附加方程的存在，保证在临界点附近的增广Jacobi矩阵非奇异，需要指出的是，在临界点以外，仍然可以使用连续法，从而越过临界点，作出PV曲线的下半段。

3.3 选择连续参数(Select the Continuation Parameter)

在预估步骤中，求解出切向量 ${\bm t}=[{d {\bm \theta}}, {d {\bm V}}, {d {\lambda}}]^T$ 之后，可以按照式(14)进行连续参数的选择：
$y_k=max{\{|t_1|, |t_2|, ..., |t_m|\}}, m=2n_1+n_2+1 \tag {14}$

式(14)中， $\large {y_k}$ 的下标 $k$ 就是在后续的预估、校正环节中所对应的索引 $k$ 。连续潮流计算的起始阶段，PV曲线较为平坦，换言之， $\lambda$ 的变化率较大，而电压幅值和相角在负载变化下则保持相当恒定，故而此时选择 $\lambda$ 作为连续参数（显然此时 $k = m$ ）。而在PV曲线到达临界点附近时，则电压幅值和相角的变化率较大，因此选择 $\lambda$ 作为连续参数则不是较好的选择，此时按照式(14)选择连续参数及索引 $k$ 。
一旦连续参数被确定下来，式(10)中，等号右侧向量的最后一项取 $+ 1$ 还是 $- 1$ 也可以随之确定：若连续参数在连续潮流计算过程中是增加的，则取 $+ 1$ ；若连续参数在连续潮流计算过程中是减小的，则取 $- 1$ 。一般而言，选择 $\lambda$ 作为连续参数时，式(10)等号右侧向量的最后一项取 $+ 1$ ，而在选择电压幅值作为连续参数时，式(10)等号右侧向量的最后一项取 $- 1$ （在PV曲线中，电压幅值一直是降低的，但是，电压相角不一定具有这样的性质）。
最后需要指出两点：(1) 每次进行“预估-校正”时，都要重新进行连续参数的选择；(2) 在同一“预估-校正”流程中，校正环节中使用的索引 $k$ 与预估环节中使用的索引 $k$ 是相同的。

3.4 关于临界点(Critical Point)与步长

在连续潮流计算中，在预估、校正环节中，需要额外注意临界点的位置及步长选择等问题。
要判断连续潮流程序是否已经计算到临界点，可以借用切向量 ${\bm \theta}}, {d {\bm V}}, {d {\lambda}}]^T$ 中，关于 ${\lambda}$ 的分量 ${\lambda}}$ 所具备的性质：在PV曲线的上半段为正值，在PV曲线的临界点处为０，越过临界点后（在PV曲线的下半段）为负值。因此在根据式(10)解出切向量 ${\bm \theta}}, {d {\bm V}}, {d {\lambda}}]^T$ 后，可以根据 ${\lambda}}$ 的符号立刻判断出连续潮流计算程序是否计算至临界点，进而可以确定程序对于PV曲线的跟踪状态。
另一方面，步长 ${\sigma}$ 的选择决定着连续潮流的效率甚至是成功与否：若步长太小，虽然更加安全，但是需要绘制大量不必要的点，造成计算时间过长；而若步长过大，则会使预估值偏离真实值，甚至可能造成校正步骤的迭代不收敛（牛顿－拉夫逊法对初值的要求较高）。然而我们事先并不知道PV曲线的走势， ${\sigma}$ 很大程度上要根据具体情况选择，因此步长控制较为困难。

3.5 连续潮流的程序流程图

综上所述，连续潮流的程序流程图如下：

4. 程序运行结果

这里选择的算例为IEEE 14节点，其数据来源为matpower 6.0。为简化分析，令所有负荷的增长系数 $k_{Li} = 1$ 、所有发电机有功出力的增长系数 $k_{Gi} = 1$ ，同时不考虑发电机因为功率越限问题而发生的PV->PQ节点等问题。选择Bus 5作为电压观测点，其PV曲线如下图所示。
在这里插入图片描述
PV曲线局部：起始阶段。此时PV曲线较为平坦，选择 ${\lambda}$ 作为连续参数，故而在上图中，体现为垂直校正：预估点与校正点的横坐标相同（因为此时连续参数 ${\lambda}$ 是固定的）。

PV曲线局部：临界点附近。此时PV曲线较为陡峭，在这一阶段，程序经过判断后选择Bus 5的电压幅值作为连续参数，故而在上图中，体现为水平校正：预估点与校正点的纵坐标相同（因为此时连续参数Bus 5的电压幅值是固定的）。
在这里插入图片描述
对于Bus 5的PV曲线则如上图所示。

5. 参考文献

[1] 赵晋泉, 张伯明. 连续潮流及其在电力系统静态稳定分析中的应用[J]. 电力系统自动化, 2002, 29(11): 91-97.
[2] AJJARAPU V, CHRISTY C. The Continuation Power Flow: A Tool for Steady State Voltage Stability Analysis[J]. IEEE Trans on Power System, 1992, 7(1): 416-423.
[3] 王锡凡. 现代电力系统分析[M]. 北京: 科学出版社, 2003.
[4] 古婷婷. 基于改进连续潮流法的静态电压稳定分析[D]. 西安: 西安理工大学.
[5] 安俊俊. 基于改进的连续潮流的静态电压稳定分析的研究[D]. 郑州: 郑州大学.
[6] AJJARAPU V. Computational Techniques for Voltage Stability Assessment and Control[M]. New York: Springer, 2006.