CN103294890A - 基于原对偶内点法的非线性规划潮流计算方法 - Google Patents

Abstract

本发明提供了一种基于原对偶内点法的非线性规划潮流计算方法，包括如下步骤：(1)构造并求解无调整量的非线性规划潮流模型，得到目标函数值；(2)判断目标函数值是否为0，若为0则直接输出潮流结果；否则执行步骤3；(3)对可调整的调整量进行调整；(4)构造并求解有调整量的非线性规划潮流模型，得到目标函数值；(5)判断目标函数值是否为0，若为0，则输出调整量的调整信息和潮流计算结果，程序结束；否则跳转到步骤3。本发明提供的基于原对偶内点法的非线性规划潮流计算方法，在传统算法不收敛时，判断潮流解是否存在，准确得到潮流有解还是无解的信息；而且在潮流解不存在的情况下，给出潮流调整措施，并能自动进行潮流数据调整，得到正常的潮流解。

Description

基于原对偶内点法的非线性规划潮流计算方法

技术领域

本发明属于电力系统仿真分析领域，具体涉及一种基于原对偶内点法的非线性规划潮流计算方法。

背景技术

潮流计算是电力系统最基本的分析计算。牛顿-拉夫逊法(简称牛顿法)是求解非线性方程组的经典方法，特别是稀疏存储技术和节点优化编号技术的引入后，因其收敛性好、计算速度快、内存需求少的优点，广泛应用于潮流计算领域，至今仍然是求解潮流计算问题的主流算法。

但牛顿法难以处理潮流计算不收敛情况。潮流计算不收敛的原因从原理上可分为两种情况：一种是潮流解(实数解)不存在，这种情况不收敛是正常的，需要计算人员调整潮流数据本身；一种是潮流解存在，但是由于计算方法本身的问题没有计算出来，这种情况是不正常的，需要更换或改进潮流算法以得到潮流解。牛顿法无法进行不收敛原因分析，如果出现不收敛情况无法判断是潮流解不存在还是算法本身的问题。

根据非线性规划的原理提出了最优乘子法，将最优乘子与常规的牛顿法相结合，在牛顿法状态变量的方向上乘以一个最优乘子，通过最优乘子信息可判断潮流解是否存在。如果潮流迭代计算过程中，最优乘子趋近于0，则说明是潮流解不存在；如果最优乘子趋近于1，则说明是算法本身问题。

现有经典的潮流计算方法--牛顿法最大的缺点是无法进行不收敛原因分析，如果出现不收敛情况无法判断是否是算法本身的问题。尤其是近年来随着我国电网规模的不断扩大，在实际电网分析中常常出现采用牛顿法及其改进方法都不能很好收敛的情况，这种情况下牛顿法本身无法做出任何处理。

最优乘子法理论上能给出潮流解存在与否的信息，但仅仅得到潮流解不存在的信息是不够的，潮流解不存在并不是潮流计算人员希望的最终结果，潮流计算人员往往只能自行调整潮流数据以期望得到正常的潮流解。在潮流解不存在的情况下，最优乘子法具有不能给出进一步的信息来指导用户进行潮流数据调整以得到潮流解的问题。

发明内容

为克服上述缺陷，本发明提供了一种基于原对偶内点法的非线性规划潮流计算方法，在传统算法不收敛时，判断潮流解是否存在，准确得到潮流有解还是无解的信息；而且在潮流解不存在的情况下，给出潮流调整措施，并能自动进行潮流数据调整，得到正常的潮流解。

为实现上述目的，本发明提供一种基于原对偶内点法的非线性规划潮流计算方法，其改进之处在于，所述方法包括如下步骤：

(1).构造并求解无调整量的非线性规划潮流模型，得到目标函数值；

(2).判断目标函数值是否为0，若为0则直接输出潮流结果；否则执行步骤3；

(3).对可调整的调整量进行调整；

(4).构造并求解有调整量的非线性规划潮流模型，得到目标函数值；

(5).判断目标函数值是否为0，若为0，则输出调整量的调整信息和潮流计算结果，程序结束；否则跳转到步骤3。

本发明提供的优选技术方案中，在所述步骤1中，无调整量的情况下只求解潮流方程；电力系统的潮流方程可表述为：

$P_i-V_i\underset{j=1}{\overset{N_B}{\mathrm{\Σ}}}{V}_j[G_{\mathrm{ij}}\cos ({\mathrm{\θ}}_i-{\mathrm{\θ}}_j)+B_{\mathrm{ij}}\sin ({\mathrm{\θ}}_i-{\mathrm{\θ}}_j)]=0(i=1,..,N_B)---\left(1\right)$

$Q_i-V_i\underset{j=1}{\overset{N_B}{\mathrm{\Σ}}}{V}_j[G_{\mathrm{ij}}\sin ({\mathrm{\θ}}_i-{\mathrm{\θ}}_j)+B_{\mathrm{ij}}\cos ({\mathrm{\θ}}_i-{\mathrm{\θ}}_j)]=0(i=1,..,N_B)---\left(2\right)$

其中，P_i、Q_i分别为第i个节点的注入有功功率、注入无功功率，G_ij、B_ij为第i个节点和第j个节点之间的节点导纳的电导分量和电纳分量，为已知量；V_i、θ_i分别为第i个节点的电压幅值和电压相角，是待求量；N_B为节点总数；

将待求量统一用X表示，潮流方程即可统一用g(X)＝0表示；

引入潮流方程误差量，构造方程误差量平方和的目标函数，可构造非线性规划潮流模型：

$\left\{\begin{array}{cc}\mathrm{Min}& \mathrm{\Σ}{\mathrm{\ϵ}}^2\\ s.t.& g\left(X\right)-\mathrm{\ϵ}=0\end{array}\right.---\left(3\right)$

其中，ε为每个潮流方程的误差量，当潮流方程有解时，每个潮流方程的ε均为0，即目标函数为0，求解的待求量X为正常的潮流结果；若目标函数不为0，则说明原始潮流方程无解，求解的待求量X为潮流方程误差最小的最小二乘解，不是正常的潮流结果，潮流解不存在。

本发明提供的第二优选技术方案中，在所述步骤1中，采用原对偶内点法求解潮流模型。

本发明提供的第三优选技术方案中，原对偶内点法是在非线性规划问题中的不等式约束加入松弛变量，如公式4所示：

$\left\{\begin{array}{cc}\mathrm{Min}& f\left(X\right)\\ s.t.& g\left(X\right)=0\\ & h\left(X\right)-S_l-h_l=0\\ & h\left(X\right)+S_u-h_u=0\end{array}\right.---\left(4\right)$

由式(4)可得拉格朗日函数：

$f\left(X\right)+{\mathrm{\λ}}^{T}g\left(X\right)+{\mathrm{\π}}_i^T(h\left(X\right)-S_l-h_l)+$

$\left(5\right)$

${\mathrm{\π}}_u^T(h\left(X\right)+S_u-h_u)-\mathrm{\μ}(\underset{i}{\mathrm{\Σ}}{\log }_e\left(S_{\mathrm{li}}\right)+\underset{i}{\mathrm{\Σ}}{\log }_e\left(S_{\mathrm{ui}}\right))$

其中，λ为拉格朗日乘子，S_l、S_u为上、下限松弛变量，π_l、π_u为对偶变量，μ为障碍因子。

根据K-T条件求解式(5)可得：

$\left\{\begin{array}{c}\mathrm{\Δ}{S}_l={\▿}_{X}h\left(X\right)\mathrm{\ΔX}\\ \mathrm{\Δ}{S}_u=-{\▿}_{X}h\left(X\right)\mathrm{\ΔX}\\ \mathrm{\Δ}{\mathrm{\π}}_l=-{{[S}_l]}^{-1}\left(\right[S_l]{\mathrm{\π}}_l+\mathrm{\μ}.e)+\left[{\mathrm{\π}}_l\right]{\▿}_{X}h\left(X\right)\mathrm{\ΔX})\\ \mathrm{\Δ}{\mathrm{\π}}_u=-{{[S}_u]}^{-1}\left(\right[S_u]{\mathrm{\π}}_u-\mathrm{\μ}.e)-\left[{\mathrm{\π}}_u\right]{\▿}_{X}h\left(X\right)\mathrm{\ΔX})\\ \left[\begin{array}{c}\mathrm{\ΔX}\\ \mathrm{\Δ\λ}\end{array}\right]=-{\left[\begin{array}{cc}H& J^T\\ J& 0\end{array}\right]}^{-1}\left[\begin{array}{c}t\\ g\left(X\right)\end{array}\right]\end{array}\right.---\left(6\right)$

其中，H阵被称为海森阵：

$H=\frac{\∂f{\left(X\right)}^2}{\∂X^2}+\mathrm{\λ}\frac{{\left[\frac{\∂g\left(X\right)}{\∂X}\right]}^T}{\∂X}+({\mathrm{\π}}_l^T+{\mathrm{\π}}_u^T){{\▿}^2}_{X}h\left(X\right)+{{\▿}^T}_{x}h\left(X\right)---\left(7\right)$

$(-{\left[S_l\right]}^{-1}{\mathrm{\π}}_l-{\left[S_u\right]}^{-1}{\mathrm{\π}}_u){\▿}_{X}h\left(X\right)$

$t=\frac{\∂f\left(X\right)}{\∂X}+J^T\mathrm{\λ}+{{\▿}^T}_{x}h\left(X\right)({\mathrm{\π}}_l+{\mathrm{\π}}_u----(8)$

$[{S_l]}^{-1}\left(\right[S_l]{\mathrm{\π}}_l+\mathrm{\μe})-{\left[S_u\right]}^{-1}(S_u{\mathrm{\π}}_u-\mathrm{\μe}))$

求解以上线性方程组，更新原、对偶变量：

$\left\{\begin{array}{c}{X}^{(K+1)}=X^{\left(K\right)}+\mathrm{\σ}{{\mathrm{\α}}_P}^{\left(K\right)}\mathrm{\Δ}{X}^{\left(K\right)}\\ {S_l}^{(K+1)}={S_l}^{\left(K\right)}+\mathrm{\σ}{{\mathrm{\α}}_P}^{\left(K\right)}\mathrm{\Δ}{S_l}^{\left(K\right)}\\ {S_u}^{(K+1)}={S_u}^{\left(K\right)}+\mathrm{\σ}{{\mathrm{\α}}_P}^{\left(K\right)}\mathrm{\Δ}{S_u}^{\left(K\right)}\\ {{\mathrm{\π}}_l}^{(K+1)}={{\mathrm{\π}}_l}^{\left(K\right)}+\mathrm{\σ}{{\mathrm{\α}}_D}^{\left(K\right)}\mathrm{\Δ}{{\mathrm{\π}}_V}^{\left(K\right)}\\ {{\mathrm{\π}}_u}^{(K+1)}={{\mathrm{\π}}_u}^{\left(K\right)}+\mathrm{\σ}{{\mathrm{\α}}_D}^{\left(K\right)}\mathrm{\Δ}{{\mathrm{\π}}_u}^{\left(K\right)}\\ {\mathrm{\λ}}^{(K+1)}={\mathrm{\λ}}^{\left(K\right)}+\mathrm{\σ}{{\mathrm{\α}}_D}^{\left(K\right)}\mathrm{\Δ}{\mathrm{\λ}}^{\left(K\right)}\end{array}\right.---\left(9\right)$

其中，σ为常数，一般0.99≤σ≤0.9995；α_P，α_D为原、对偶变量的步长因子，其设置是为了保证所有变量可行性，一般由下式确定：

${{\mathrm{\&alpha;}}_P}^{\left(K\right)}=\min \{\underset{{{\mathrm{\&Delta;S}}_{\mathrm{li}}}^{\left(K\right)}<0}{\min }(-\frac{{S_{\mathrm{li}}}^{\left(K\right)}}{{{\mathrm{\&Delta;S}}_{\mathrm{li}}}^{\left(K\right)}}),\underset{{{\mathrm{\&Delta;S}}_{\mathrm{ui}}}^{\left(K\right)}<0}{\min }(-\frac{{S_{\mathrm{ui}}}^{\left(K\right)}}{{{\mathrm{\&Delta;S}}_{\mathrm{ui}}}^{\left(K\right)}}),1\}---\left(10\right)$

${{\mathrm{\&alpha;}}_D}^{\left(K\right)}=\min \{\underset{{{\mathrm{\&Delta;\&pi;}}_{\mathrm{li}}}^{\left(K\right)}>0}{\min }(-\frac{{{\mathrm{\&pi;}}_{\mathrm{li}}}^{\left(K\right)}}{{{\mathrm{\&Delta;\&pi;}}_{\mathrm{li}}}^{\left(K\right)}}),\underset{{{\mathrm{\&Delta;\&pi;}}_{\mathrm{ui}}}^{\left(K\right)}<0}{\min }(-\frac{{{\mathrm{\&pi;}}_{\mathrm{vi}}}^{\left(K\right)}}{{{\mathrm{\&Delta;\&pi;}}_i}^{\left(K\right)}}),1\}---\left(11\right)$

障碍因子μ的选取直接影响算法的收敛性，无论原变量还是对偶变量，其转移方向都受障碍因子影响；障碍因子采用互补间隙法确定：

${\mathrm{\&rho;}}^{(K+1)}=-{S_l}^{(K+1)T}{{{\mathrm{\&pi;}}_l}^{(K+1)}}_{}+{S_u}^{(K+1)T}{{\mathrm{\&pi;}}_u}^{(K+1)}$

μ^(K+1)＝ρ^(K+1)/2nβ  (13)

其中，ρ为互补间隙，β为加速因子，n为同一类型变量的数目。

本发明提供的第四优选技术方案中，在所述步骤3中，调整量包括：变压器变比、发电机有功无功功率，负荷有功无功功率和无功补偿容量。

本发明提供的第五优选技术方案中，对变压器的变比的调整包括：取变比的倒数(1/T_k)为调整量；目标函数、节点平衡方程约束与变比调整量的关系为一次函数关系，避免变比为分母的强非线性问题；

设(i，k)支路为变压器，变比为1∶T_k，令T_rec＝1/T_k，节点i注入的有功功率为：

$P_i(V,\mathrm{\&theta;})=V_i\underset{\underset{j\&NotEqual;k}{j=1}}{\overset{N_B}{\mathrm{\&Sigma;}}}{V}_j(G_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}})$

$+\frac{V_i{V}_k}{T_k}(G_{\mathrm{ik}}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}+B_{\mathrm{ik}}\sin {\mathrm{\&theta;}}_{\mathrm{ik}})---\left(14\right)$

则潮流等式方程对T_k、T_rec的二阶导数为：

$\frac{{\&PartialD;}^2{P}_i}{\&PartialD;{T_k}^2}=\frac{2V_i{V}_k}{{T_k}^3}(G_{\mathrm{ik}}\cos ({\mathrm{\&theta;}}_{\mathrm{ik}}+\mathrm{\&delta;})+B_{\mathrm{ik}}\sin ({\mathrm{\&theta;}}_{\mathrm{ik}}+\mathrm{\&delta;}))---\left(15\right)$

${\&PartialD;}^2{P}_i/\&PartialD;{T_{\mathrm{rec}}}^2=0---\left(16\right)$

可见，采用1/T_k为调整量避免了T_k为调整量的强非线性。

本发明提供的第六优选技术方案中，发电机有功功率Pg和无功功率Qg作为调整量，根据发电机功率的上下限约束条件对发电机有功无功功率进行调整：

Pg_min≤Pg≤Pg_max  (17)

Qg_min≤Qg≤Qg_max  (18)

本发明提供的第六优选技术方案中，负荷有功功率Pl和无功功率Ql作为调整量，根据负荷功率的上下限约束条件对负荷有功、无功功率进行调整：

Pl_min≤Pl≤Pl_max  (19)

Ql_min≤Ql≤Ql_max  (20)

本发明提供的第七优选技术方案中，无功补偿包括：并联电容电抗器的投切和未知的负荷侧无功补偿。

本发明提供的第八优选技术方案中，对并联电容电抗器的投切的调整通过控制支路的导纳值实现；设i支路为可控支路，阻抗为R+jX，导纳值G+jB，则节点i注入的功率为：

$P_i(V,\mathrm{\&theta;})=V_i\underset{\underset{j\&NotEqual;i}{j=1}}{\overset{N_B}{\mathrm{\&Sigma;}}}{V}_j(G_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}})+{V_i}^2(-\underset{\underset{j\&NotEqual;i}{j=1}}{\overset{N_B}{\mathrm{\&Sigma;}}}{G}_{\mathrm{ij}}+G_i)---\left(21\right)$

其中， $G=\frac{R}{R^2+X^2}$

则潮流有功等式方程对X、B的一阶导数为：

$\frac{\&PartialD;P_i}{\&PartialD;X}=-\frac{2\mathrm{RX}}{{(R^2+X^2)}^2}{V_i}^2---\left(22\right)$

$\&PartialD;P_i/\&PartialD;B=0---\left(23\right)$

取导纳B为调整量避免了X为调整量的强非线性。

本发明提供的第九优选技术方案中，利用负荷的无功功率Ql，根据补偿设备总的容量限制，调整未知的负荷侧的无功补偿。

本发明提供的第十优选技术方案中，在所述步骤4中，通过优化各个调整量使潮流调整到有解的状态，其数学模型如下：

$\left\{\begin{array}{cc}\mathrm{Min}& \mathrm{\&Sigma;}{\mathrm{\&epsiv;}}^2\\ s.t.& g(X,u)-\mathrm{\&epsiv;}=0\\ & h_l\&le;h\left(u\right)\&le;h_u\end{array}\right.---\left(24\right)$

其中，X为系统待求量，u为调整量；g(X，u)为潮流方程；h(u)为不等式约束方程，即调整量的调节范围约束，h_l、h_u为调整量的下限和上限。

与现有技术比，本发明提供的一种基于原对偶内点法的非线性规划潮流计算方法，解决判断潮流解是否存在的问题以及潮流解不存在的潮流调整问题。与传统的牛顿法相比，可以准确得到潮流有解还是无解的信息；与最优乘子法相比，可以在潮流解不存在的情况下，可以给出潮流调整措施，并能自动进行潮流数据调整，得到正常的潮流解。

附图说明

图1为一般矩阵排列方式。

图2为改进后矩阵排列方式。

图3为基于原对偶内点法的非线性规划潮流计算方法的流程图。

具体实施方式

如图3所示，一种基于原对偶内点法的非线性规划潮流计算方法，包括如下步骤：

(1).构造并求解无调整量的非线性规划潮流模型，得到目标函数值；

(2).判断目标函数值是否为0，若为0则直接输出潮流结果；否则执行步骤3；

(3).对可调整的调整量进行调整；

(4).构造并求解有调整量的非线性规划潮流模型，得到目标函数值；

(5).判断目标函数值是否为0，若为0，则输出调整量的调整信息和潮流计算结果，程序结束；否则跳转到步骤3。

在所述步骤1中，无调整量的情况下只求解潮流方程；电力系统的潮流方程可表述为：

$P_i-V_i\underset{j=1}{\overset{N_B}{\mathrm{\&Sigma;}}}{V}_j[G_{\mathrm{ij}}\cos ({\mathrm{\&theta;}}_i-{\mathrm{\&theta;}}_j)+B_{\mathrm{ij}}\sin ({\mathrm{\&theta;}}_i-{\mathrm{\&theta;}}_j)]=0(i=1,..,N_B)---\left(1\right)$

$P_i-V_i\underset{j=1}{\overset{N_B}{\mathrm{\&Sigma;}}}{V}_j[G_{\mathrm{ij}}\cos ({\mathrm{\&theta;}}_i-{\mathrm{\&theta;}}_j)-B_{\mathrm{ij}}\sin ({\mathrm{\&theta;}}_i-{\mathrm{\&theta;}}_j)]=0(i=1,..,N_B)---\left(2\right)$

其中，P_i、Q_i分别为第i个节点的注入有功功率、注入无功功率，G_ij、B_ij为第i个节点和第j个节点之间的节点导纳的电导分量和电纳分量，为已知量；V_i、θ_i分别为第i个节点的电压幅值和电压相角，是待求量；N_B为节点总数；

将待求量统一用X表示，潮流方程即可统一用g(X)＝0表示；

引入潮流方程误差量，构造方程误差量平方和的目标函数，可构造非线性规划潮流模型：

$\left\{\begin{array}{cc}\mathrm{Min}& \mathrm{\&Sigma;}{\mathrm{\&epsiv;}}^2\\ s.t.& g\left(X\right)-\mathrm{\&epsiv;}=0\end{array}\right.---\left(3\right)$

其中，ε为每个潮流方程的误差量，当潮流方程有解时，每个潮流方程的ε均为0，即目标函数为0，求解的待求量X为正常的潮流结果；若目标函数不为0，则说明原始潮流方程无解，求解的待求量X为潮流方程误差最小的最小二乘解，不是正常的潮流结果，潮流解不存在。

在所述步骤1中，采用原对偶内点法求解潮流模型。

原对偶内点法是在非线性规划问题中的不等式约束加入松弛变量，如公式4所示：

$\left\{\begin{array}{cc}\mathrm{Min}& f\left(X\right)\\ s.t.& g\left(X\right)=0\\ & h\left(X\right)-S_l-h_l=0\\ & h\left(X\right)+S_u-h_u=0\end{array}\right.---\left(4\right)$

由式(4)可得拉格朗日函数：

$f\left(X\right)+{\mathrm{\&lambda;}}^{T}g\left(X\right)+{\mathrm{\&pi;}}_i^T(h\left(X\right)-S_l-h_l)+$

$\left(5\right)$

${\mathrm{\&pi;}}_u^T(h\left(X\right)+S_u-h_u)-\mathrm{\&mu;}(\underset{i}{\mathrm{\&Sigma;}}{\log }_e\left(S_{\mathrm{li}}\right)+\underset{i}{\mathrm{\&Sigma;}}{\log }_e\left(S_{\mathrm{ui}}\right))$

其中，λ为拉格朗日乘子，S_l、S_u为上、下限松弛变量，π_l、π_u为对偶变量，μ为障碍因子。

根据K-T条件求解式(5)可得：

$\left\{\begin{array}{c}\mathrm{\&Delta;}{S}_l={\&dtri;}_{X}h\left(X\right)\mathrm{\&Delta;X}\\ \mathrm{\&Delta;}{S}_u=-{\&dtri;}_{X}h\left(X\right)\mathrm{\&Delta;X}\\ \mathrm{\&Delta;}{\mathrm{\&pi;}}_l=-{{[S}_l]}^{-1}\left(\right[S_l]{\mathrm{\&pi;}}_l+\mathrm{\&mu;}.e)+\left[{\mathrm{\&pi;}}_l\right]{\&dtri;}_{X}h\left(X\right)\mathrm{\&Delta;X})\\ \mathrm{\&Delta;}{\mathrm{\&pi;}}_u=-{{[S}_u]}^{-1}\left(\right[S_u]{\mathrm{\&pi;}}_u-\mathrm{\&mu;}.e)-\left[{\mathrm{\&pi;}}_u\right]{\&dtri;}_{X}h\left(X\right)\mathrm{\&Delta;X})\\ \left[\begin{array}{c}\mathrm{\&Delta;X}\\ \mathrm{\&Delta;\&lambda;}\end{array}\right]=-{\left[\begin{array}{cc}H& J^T\\ J& 0\end{array}\right]}^{-1}\left[\begin{array}{c}t\\ g\left(X\right)\end{array}\right]\end{array}\right.---\left(6\right)$

其中，H阵被称为海森阵：

$H=\frac{\&PartialD;f{\left(X\right)}^2}{\&PartialD;X^2}+\mathrm{\&lambda;}\frac{{\left[\frac{\&PartialD;g\left(X\right)}{\&PartialD;X}\right]}^T}{\&PartialD;X}+({\mathrm{\&pi;}}_l^T+{\mathrm{\&pi;}}_u^T){{\&dtri;}^2}_{X}h\left(X\right)+{{\&dtri;}^T}_{x}h\left(X\right)---\left(7\right)$

$(-{\left[S_l\right]}^{-1}{\mathrm{\&pi;}}_l-{\left[S_u\right]}^{-1}{\mathrm{\&pi;}}_u){\&dtri;}_{X}h\left(X\right)$

$t=\frac{\&PartialD;f\left(X\right)}{\&PartialD;X}+J^T\mathrm{\&lambda;}+{{\&dtri;}^T}_{x}h\left(X\right)({\mathrm{\&pi;}}_l+{\mathrm{\&pi;}}_u----(8)$

${{[S}_l]}^{-1}\left(\right[S_l]{\mathrm{\&pi;}}_l+\mathrm{\&mu;e})-{\left[S_u\right]}^{-1}(S_u{\mathrm{\&pi;}}_u-\mathrm{\&mu;e}))$

求解以上线性方程组，更新原、对偶变量：

$\left\{\begin{array}{c}{X}^{(K+1)}=X^{\left(K\right)}+\mathrm{\&sigma;}{{\mathrm{\&alpha;}}_P}^{\left(K\right)}\mathrm{\&Delta;}{X}^{\left(K\right)}\\ {S_l}^{(K+1)}={S_l}^{\left(K\right)}+\mathrm{\&sigma;}{{\mathrm{\&alpha;}}_P}^{\left(K\right)}\mathrm{\&Delta;}{S_l}^{\left(K\right)}\\ {S_u}^{(K+1)}={S_u}^{\left(K\right)}+\mathrm{\&sigma;}{{\mathrm{\&alpha;}}_P}^{\left(K\right)}\mathrm{\&Delta;}{S_u}^{\left(K\right)}\\ {{\mathrm{\&pi;}}_l}^{(K+1)}={{\mathrm{\&pi;}}_l}^{\left(K\right)}+\mathrm{\&sigma;}{{\mathrm{\&alpha;}}_D}^{\left(K\right)}\mathrm{\&Delta;}{{\mathrm{\&pi;}}_V}^{\left(K\right)}\\ {{\mathrm{\&pi;}}_u}^{(K+1)}={{\mathrm{\&pi;}}_u}^{\left(K\right)}+\mathrm{\&sigma;}{{\mathrm{\&alpha;}}_D}^{\left(K\right)}\mathrm{\&Delta;}{{\mathrm{\&pi;}}_u}^{\left(K\right)}\\ {\mathrm{\&lambda;}}^{(K+1)}={\mathrm{\&lambda;}}^{\left(K\right)}+\mathrm{\&sigma;}{{\mathrm{\&alpha;}}_D}^{\left(K\right)}\mathrm{\&Delta;}{\mathrm{\&lambda;}}^{\left(K\right)}\end{array}\right.---\left(9\right)$

其中，σ为常数，一般0.99≤σ≤0.9995；α_P，α_D为原、对偶变量的步长因子，其设置是为了保证所有变量可行性，一般由下式确定：

${{\mathrm{\&alpha;}}_P}^{\left(K\right)}=\min \{\underset{{{\mathrm{\&Delta;S}}_{\mathrm{li}}}^{\left(K\right)}<0}{\min }(-\frac{{S_{\mathrm{li}}}^{\left(K\right)}}{{{\mathrm{\&Delta;S}}_{\mathrm{li}}}^{\left(K\right)}}),\underset{{{\mathrm{\&Delta;S}}_{\mathrm{ui}}}^{\left(K\right)}<0}{\min }(-\frac{{S_{\mathrm{ui}}}^{\left(K\right)}}{{{\mathrm{\&Delta;S}}_{\mathrm{ui}}}^{\left(K\right)}}),1\}---\left(10\right)$

${{\mathrm{\&alpha;}}_D}^{\left(K\right)}=\min \{\underset{{{\mathrm{\&Delta;\&pi;}}_{\mathrm{li}}}^{\left(K\right)}>0}{\min }(-\frac{{{\mathrm{\&pi;}}_{\mathrm{li}}}^{\left(K\right)}}{{{\mathrm{\&Delta;\&pi;}}_{\mathrm{li}}}^{\left(K\right)}}),\underset{{{\mathrm{\&Delta;\&pi;}}_{\mathrm{ui}}}^{\left(K\right)}<0}{\min }(-\frac{{{\mathrm{\&pi;}}_{\mathrm{vi}}}^{\left(K\right)}}{{{\mathrm{\&Delta;\&pi;}}_i}^{\left(K\right)}}),1\}---\left(11\right)$

障碍因子μ的选取直接影响算法的收敛性，无论原变量还是对偶变量，其转移方向都受障碍因子影响；障碍因子采用互补间隙法确定：

${\mathrm{\&rho;}}^{(K+1)}=-{S_l}^{(K+1)T}{{{\mathrm{\&pi;}}_l}^{(K+1)}}_{}+{S_u}^{(K+1)T}{{\mathrm{\&pi;}}_u}^{(K+1)}---\left(12\right)$

μ^(K+1)＝ρ^(K+1)/2nβ  (13)

其中，ρ为互补间隙，β为加速因子，n为同一类型变量的数目。

在所述步骤3中，调整量包括：变压器变比、发电机有功无功功率，负荷有功无功功率和无功补偿容量。

对变压器的变比的调整包括：取变比的倒数(1/T_k)为调整量；目标函数、节点平衡方程约束与变比调整量的关系为一次函数关系，避免变比为分母的强非线性问题；

设(i，k)支路为变压器，变比为1∶T_k，令T_rec＝1/T_k，节点i注入的有功功率为：

$P_i(V,\mathrm{\&theta;})=V_i\underset{\underset{j\&NotEqual;k}{j=1}}{\overset{N_B}{\mathrm{\&Sigma;}}}{V}_j(G_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}})$

$+\frac{V_i{V}_k}{T_k}(G_{\mathrm{ik}}\cos {\mathrm{\&theta;}}_{\mathrm{ik}}+B_{\mathrm{ik}}\sin {\mathrm{\&theta;}}_{\mathrm{ik}})---\left(14\right)$

则潮流等式方程对T_k、T_rec的二阶导数为：

$\frac{{\&PartialD;}^2{P}_i}{\&PartialD;{T_k}^2}=\frac{2V_i{V}_k}{{T_k}^3}(G_{\mathrm{ik}}\cos ({\mathrm{\&theta;}}_{\mathrm{ik}}+\mathrm{\&delta;})+B_{\mathrm{ik}}\sin ({\mathrm{\&theta;}}_{\mathrm{ik}}+\mathrm{\&delta;}))---\left(15\right)$

${\&PartialD;}^2{P}_i/\&PartialD;{T_{\mathrm{rec}}}^2=0---\left(16\right)$

可见，采用1/T_k为调整量避免了T_k为调整量的强非线性。

本发明提供的第六优选技术方案中，发电机有功功率Pg和无功功率Qg作为调整量，根据发电机功率的上下限约束条件对发电机有功无功功率进行调整：

Pg_min≤Pg≤Pg_max  (17)

Qg_min≤Qg≤Qg_max  (18)

本发明提供的第六优选技术方案中，负荷有功功率Pl和无功功率Ql作为调整量，根据负荷功率的上下限约束条件对负荷有功、无功功率进行调整：

Pl_min≤Pl≤Pl_max  (19)

Ql_min≤Ql≤Ql_max  (20)

本发明提供的第七优选技术方案中，无功补偿包括：并联电容电抗器的投切和未知的负荷侧无功补偿。

对并联电容电抗器的投切的调整通过控制支路的导纳值实现；设i支路为可控支路，阻抗为R+jX，导纳值G+jB，则节点i注入的功率为：

$P_i(V,\mathrm{\&theta;})=V_i\underset{\underset{j\&NotEqual;i}{j=1}}{\overset{N_B}{\mathrm{\&Sigma;}}}{V}_j(G_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}})+{V_i}^2(-\underset{\underset{j\&NotEqual;i}{j=1}}{\overset{N_B}{\mathrm{\&Sigma;}}}{G}_{\mathrm{ij}}+G_i)---\left(21\right)$

其中， $G=\frac{R}{R^2+X^2}$

则潮流有功等式方程对X、B的一阶导数为：

$\frac{\&PartialD;P_i}{\&PartialD;X}=-\frac{2\mathrm{RX}}{{(R^2+X^2)}^2}{V_i}^2---\left(22\right)$

$\&PartialD;P_i/\&PartialD;B=0---\left(23\right)$

取导纳B为调整量避免了X为调整量的强非线性。

利用负荷的无功功率Ql，根据补偿设备总的容量限制，调整未知的负荷侧的无功补偿。

在所述步骤4中，通过优化各个调整量使潮流调整到有解的状态，其数学模型如下：

$\left\{\begin{array}{cc}\mathrm{Min}& \mathrm{\&Sigma;}{\mathrm{\&epsiv;}}^2\\ s.t.& g(X,u)-\mathrm{\&epsiv;}=0\\ & h_l\&le;h\left(u\right)\&le;h_u\end{array}\right.---\left(24\right)$

其中，X为系统待求量，u为调整量；g(X，u)为潮流方程；h(u)为不等式约束方程，即调整量的调节范围约束，h_l、h_u为调整量的下限和上限。

有调整量时的原对偶内点法改进

调整量增加后，迭代计算中LU分解时的非零注入元会急剧增加，计算时间缓慢，效率会大大降低，这主要是由于海森阵元素位置未优化。

如图1所示，一般海森阵H形成时，

ΔV、Δλp(有功潮流方程的λ)、Δλq(无功潮流方程的λ)所在行是作为一个分块4*4矩阵整体排列，但ΔTk、ΔB等调整量却放在矩阵最后，对整个矩阵的稀疏性造成很大影响。

如图2所示，实际上，ΔTk只与变压器两侧的母线有关(可取任意一侧)，ΔB只与所联母线有关，两者可移到对应母线所在ΔV、Δλp、Δλq的4*4分块矩阵上，形成5*5、6*6或更大的分块阵。

需要声明的是，本发明内容及具体实施方式意在证明本发明所提供技术方案的实际应用，不应解释为对本发明保护范围的限定。本领域技术人员在本发明的精神和原理启发下，可作各种修改、等同替换、或改进。但这些变更或修改均在申请待批的保护范围内。

Claims (12)

1.一种基于原对偶内点法的非线性规划潮流计算方法，其特征在于，所述方法包括如下步骤： 

(1).构造并求解无调整量的非线性规划潮流模型，得到目标函数值； 

(2).判断目标函数值是否为0，若为0则直接输出潮流结果；否则执行步骤3； 

(3).对可调整的调整量进行调整； 

(4).构造并求解有调整量的非线性规划潮流模型，得到目标函数值； 

(5).判断目标函数值是否为0，若为0，则输出调整量的调整信息和潮流计算结果，程序结束；否则跳转到步骤3。 

2.根据权利要求1所述的方法，其特征在于，在所述步骤1中，无调整量的情况下只求解潮流方程；电力系统的潮流方程可表述为： 

其中，P_i、Q_i分别为第i个节点的注入有功功率、注入无功功率，G_ij、B_ij为第i个节点和第j个节点之间的节点导纳的电导分量和电纳分量，为已知量；V_i、θ_i分别为第i个节点的电压幅值和电压相角，是待求量；N_B为节点总数； 

将待求量统一用X表示，潮流方程即可统一用g(X)＝0表示； 

引入潮流方程误差量，构造方程误差量平方和的目标函数，可构造非线性规划潮流模型： 

其中，ε为每个潮流方程的误差量，当潮流方程有解时，每个潮流方程的ε均为0，即目标函数为0，求解的待求量X为正常的潮流结果；若目标函数不为0， 则说明原始潮流方程无解，求解的待求量X为潮流方程误差最小的最小二乘解，不是正常的潮流结果，潮流解不存在。 

3.根据权利要求1所述的方法，其特征在于，在所述步骤1中，采用原对偶内点法求解潮流模型。 

4.根据权利要求3所述的方法，其特征在于，原对偶内点法是在非线性规划问题中的不等式约束加入松弛变量，如公式4所示： 

由式(4)可得拉格朗日函数： 

其中，λ为拉格朗日乘子，S_l、S_u为上、下限松弛变量，π_l、π_u为对偶变量，μ为障碍因子。 

根据K-T条件求解式(5)可得： 

其中，H阵被称为海森阵： 

求解以上线性方程组，更新原、对偶变量： 

其中，σ为常数，一般0.99≤σ≤0.9995；α_P，α_D为原、对偶变量的步长因子，其设置是为了保证所有变量可行性，一般由下式确定： 

障碍因子μ的选取直接影响算法的收敛性，无论原变量还是对偶变量，其转移方向都受障碍因子影响；障碍因子采用互补间隙法确定： 

μ^(K+1)＝ρ^(K+1)/2nβ  (13) 

其中，ρ为互补间隙，β为加速因子，n为同一类型变量的数目。 

5.根据权利要求1所述的方法，其特征在于，在所述步骤3中，调整量包括：变压器变比、发电机有功无功功率，负荷有功无功功率和无功补偿容量。 

6.根据权利要求5所述的方法，其特征在于，对变压器的变比的调整包括：取变比的倒数(1/T_k)为调整量；目标函数、节点平衡方程约束与变比调整量的关系为一次函数关系，避免变比为分母的强非线性问题； 

设(i，k)支路为变压器，变比为1∶T_k，令T_rec＝1/T_k，节点i注入的有功功率为： 

则潮流等式方程对T_k、T_rec的二阶导数为： 

可见，采用1/T_k为调整量避免了T_k为调整量的强非线性。 

7.根据权利要求5所述的方法，其特征在于，发电机有功功率Pg和无功功率Qg作为调整量，根据发电机功率的上下限约束条件对发电机有功无功功率进行调整： 

Pg_min≤Pg≤Pg_max   (17) 

Qg_min≤Qg≤Qg_max   (18) 。

8.根据权利要求5所述的方法，其特征在于，负荷有功功率Pl和无功功率Ql作为调整量，根据负荷功率的上下限约束条件对负荷有功、无功功率进行调整： 

Pl_min≤Pl≤Pl_max   (19) 

Ql_min≤Ql≤Ql_max   (20) 。

9.根据权利要求5所述的方法，其特征在于，无功补偿包括：并联电容电抗器的投切和未知的负荷侧无功补偿。 

10.根据权利要求9所述的方法，其特征在于，对并联电容电抗器的投切的调整通过控制支路的导纳值实现；设i支路为可控支路，阻抗为R+jX，导纳值G+jB，则节点i注入的功率为： 

   (21) 

其中，

则潮流有功等式方程对X、B的一阶导数为： 

   (22) 

   (23) 

取导纳B为调整量避免了X为调整量的强非线性。 

11.根据权利要求9所述的方法，其特征在于，利用负荷的无功功率Ql，根据补偿设备总的容量限制，调整未知的负荷侧的无功补偿。 

12.根据权利要求1所述的方法，其特征在于，在所述步骤4中，通过优化各个调整量使潮流调整到有解的状态，其数学模型如下： 

   (24) 

其中，X为系统待求量，u为调整量；g(X,u)为潮流方程；h(u)为不等式约束方程，即调整量的调节范围约束，h_l、h_u为调整量的下限和上限。 

Images