CN103366097B - 一种基于类扩展变量内点法的最优潮流计算方法 - Google Patents

Abstract

本发明提供一种基于类扩展变量内点法的最优潮流计算方法，包括以下步骤：引入非凸非线性规划中类扩展变量，并确定引入类扩展变量后的拉格朗日函数；构建引入类扩展变量后的拉格朗日函数的KKT条件方程；求解KKT条件方程；根据预测-校正法进行最优潮流计算。本发明在原优化问题因约束条件限制无解时给出约束需放开的幅度以及对应的最优解，且在某些用传统预测-校正法不能收敛的优化问题中也能有解，扩大了优化问题的解空间，提高了收敛性。本发明对新方法在大规模方程中的数值问题给出了实用化的处理，其有效性通过电力系统潮流优化问题的算例得到了验证。

Description

一种基于类扩展变量内点法的最优潮流计算方法

技术领域

本发明属于电力系统技术领域，具体涉及一种基于类扩展变量内点法的最优潮流计算方法。

背景技术

内点算法是指每个迭代点都是可行域的内点的算法，Karmarkar于1984年首次提出了求解线性规划的内点算法。1989年后内点法求解非线性规划逐渐成为研究的热点，并在电力系统优化问题中得到大量应用。

非线性内点法具有较好的收敛性和较快的计算速度，已经成为求解电力系统优化问题的有力工具。但是，不论原始对偶内点法，还是预测-校正内点法、多中心校正内点法或者其它内点法的改进算法，其解空间并不理想，具体表现在以下两个方面：

（1）不等式约束的解空间和等式约束的解空间无重合点，使计算溢出。在优化计算前验证这点是困难的，特别是对大规模的优化问题。更严格的约束条件往往意味着更好的物理性能，因此内点法的使用者不得不在保证收敛和提高性能间进行小心的平衡。

（2）不等式约束的解空间和等式约束的解空间虽有重合点，但因数值稳定性低使优化问题无解。

为了缓解上述两方面的困难，在计算过程中扩大内点法的解空间成为重要的手段。有些学者利用已知的可能越限或允许越限的不等式约束信息，进行部分约束松弛，取得一定效果。然而，导致内点法不收敛的关键约束经常不能在计算前获得，使得预设的约束松弛难以从根本上改善内点法优化算法的收敛性。另外一类方法是在求解最优潮流问题时采用了一种参数化电力需求或运行限制的方法。虽能扩大等式和不等式约束的解空间，但是由于限定新增的参数不小于0，对迭代的收敛性有一定的影响；并且也未考虑不同类别的约束对计算过程以及物理实质影响的差异，采用了单一的参数化因子。

发明内容

为了克服上述现有技术的不足，本发明提供一种基于类扩展变量内点法的最优潮流计算方法，在原优化问题因约束条件限制无解时给出约束需放开的幅度以及对应的最优解，且在某些用传统预测-校正法不能收敛的优化问题中也能有解，扩大了优化问题的解空间，提高了收敛性。本发明对新方法在大规模方程中的数值问题给出了实用化的处理，其有效性通过电力系统潮流优化问题的算例得到了验证。

为了实现上述发明目的，本发明采取如下技术方案：

提供一种基于类扩展变量内点法的最优潮流计算方法，所述方法包括以下步骤：

步骤1：引入非凸非线性规划中类扩展变量，并确定引入类扩展变量后的拉格朗日函数；

步骤2：构建引入类扩展变量后的拉格朗日函数的KKT条件方程；

步骤3：求解KKT条件方程；

步骤4：根据预测-校正法进行最优潮流计算。

所述非凸非线性规划表示为：

$\left\{\begin{array}{c}\underset{x}{\min f}\left(x\right)\\ s.t.h\left(x\right)=0\\ \underset{\‾}{g}\≤g\left(x\right)\≤\stackrel{\‾}{g}\end{array}\right.---\left(1\right)$

其中，x为自变量构成的列向量，f(x)为优化目标函数，h(x)为等式约束矩阵，g(x)为不等式约束矩阵，和g分别为不等式约束的上限和下限列向量。

将不等式约束矩阵按类别增加类扩展变量，并在目标函数中增加带罚因子的扩展变量二次项，则式（1）变为：

$\left\{\begin{array}{c}\underset{x}{\min }(f\left(x\right)+\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ξ}}_k({\mathrm{\α}}_k^2+{\mathrm{\β}}_k^2))\\ s.t.h\left(x\right)=0\\ \underset{\‾}{g}-\mathrm{T\α}\≤g\left(x\right)\≤\stackrel{\‾}{g}+\mathrm{T\β}\end{array}---\left(2\right)\right.$

n为设置扩展变量的不等式约束的种类，α和β为类扩展变量向量，分别代表下扩展变量和上扩展变量，维度与设置扩展变量的不等式约束的种类相等，即为n；α_k、β_k分别表示α、β的第k个元素；ξ_k为第k个类扩展变量二次项的罚因子，T为不等式约束和类扩展变量的联系矩阵，其为r×n维矩阵，形式如下：

$T=\left[\begin{array}{cccc}1& 0& \·\·\·& 0\\ 1& 0& \·\·\·& 0\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ 0& 0& \·\·\·& 1\\ 0& 0& \·\·\·& 0\end{array}\right]---\left(3\right)$

若式（3）中第i个不等式约束中增加了第j个类扩展变量，则T_i,j为1，T第i行的其它元素均为0，且T中非0元均设为1；

根据类扩展变量内点法，用对数障碍函数化不等式约束为等式约束，用拉格朗日乘子法处理等式约束，与式（3）对应引入类扩展变量后的拉格朗日函数为：

$\begin{array}{c}{L}_g=f\left(x\right)+\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ξ}}_k({\mathrm{\α}}_k^2+{\mathrm{\β}}_k^2)-y^{T}h\left(x\right)-\\ z^T(g\left(x\right)-1-\underset{\‾}{g}+\mathrm{T\α})-\mathrm{\μ}\underset{j=1}{\overset{r}{\mathrm{\Σ}}}\ln l_j-\\ w^T(g\left(x\right)+u-\stackrel{\‾}{g}-\mathrm{T\β})-\mathrm{\μ}\underset{j=1}{\overset{r}{\mathrm{\Σ}}}\ln u_j\end{array}---\left(4\right)$

其中，L_g为拉格朗日函数；μ为障碍参数；r为不等式约束的数目；l和u均为松弛因子列向量；y、z和w分别为对应等式约束和不等式约束的拉格朗日乘子列向量；l_j和u_j分别表示向量l和u的第j个元素；其中列向量l、u、z和w的维度等于不等式约束的数目，即为r，y的维度为等式约束的数目m。

所述步骤2中，构建引入类扩展变量后的拉格朗日函数的KKT条件方程，有：

L_x＝▽_xf(x)-▽_xh(x)y-▽_xg(x)(z+w)＝0（5）

L_y＝h(x)＝0（6）

L_z＝g(x)-l-g+Tα＝0（7）

$L_w=g\left(x\right)+u-\stackrel{\‾}{g}-\mathrm{T\β}=0---\left(8\right)$

L_l＝LZE-μE（9）

L_u＝UWE+μE（10）

L_α＝2ξα-T^Tz（11）

L_β＝2ξβ+T^Tw（12）

式（5）-（12）中，L_x、L_y、L_z、L_w、L_l、L_u、L_α和L_β表示拉格朗日函数对变量x、y、z、w、l、u、α和β分别取偏导数，形成的矩阵称为雅可比矩阵；其中，L＝diag(l₁,l₂,…,l_r)，Z＝diag(z₁,z₂,…,z_r)，U＝diag(u₁,u₂,…,u_r)，W＝diag(w₁,w₂,…,w_r)，ξ＝diag(ξ₁,ξ₂,…,ξ_n)，E＝[1,1,…,1]^T，E为r维向量，▽_x为对自变量x的一阶导数。

所述步骤3中，用牛顿法求解式（5）-（12），先进行线性化处理，有：

$-[{\▿}_x^{2}f\left(x\right)-{\▿}_x^{2}h\left(x\right)y-{\▿}_x^{2}g\left(x\right)(z+w)]\mathrm{\Δx}+{\▿}_{x}h\left(x\right)\mathrm{\Δy}+{\▿}_{x}g\left(x\right)(\mathrm{\Δz}+\mathrm{\Δw})=L_x---\left(13\right)$

${\▿}_x^{T}h\left(x\right)\mathrm{\Δx}=-L_y---\left(14\right)$

${\▿}_x^{T}g\left(x\right)\mathrm{\Δx}-\mathrm{\Δ}1+\mathrm{T\Δ\α}=-L_z---\left(15\right)$

${\▿}_x^{T}g\left(x\right)\mathrm{\Δx}+\mathrm{\Δu}-\mathrm{T\Δ\β}=-L_w---\left(16\right)$

LΔz+ZΔl＝-L_l（17）

uΔw+wΔu＝-L_u（18）

2ξΔα-T^TΔz＝-L_α（19）

2ξΔβ+T^TΔw＝-L_β（20）

式（13）-（20）中，Δx、Δl、Δu、Δy、Δz、Δw、Δα和Δβ分别为x、l、u、y、z、w、α和β的微分列向量，表示x的二阶偏导数形成海森矩阵；经过推导，得到如下解耦的线性方程组：

$\left[\begin{array}{cccc}{H}^{\′}& {\▿}_{x}h\left(x\right)& F1& F2\\ {\mathrm{\Δ}}_x^{T}h\left(x\right)& 0& 0& 0\\ -{F1}^T& 0& F3& 0\\ -F{2}^T& 0& 0& F4\end{array}\right]\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δ\α}\\ \mathrm{\Δ\β}\end{array}\right]=\left[\begin{array}{c}{L}_x^{\′}\\ -L_y\\ -L_{\mathrm{\α}}^{\′}\\ -L_{\mathrm{\β}}^{\′}\end{array}\right]---\left(21\right)$

$\left[\begin{array}{cc}L& Z\\ I& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\Δz}\\ \mathrm{\Δl}\end{array}\right]=\left[\begin{array}{c}-L_1\\ -L_z^{\′}-Z{L}^{-1}({\▿}_x^{T}g\left(x\right)\mathrm{\Δx}+\mathrm{T\Δ\α})\end{array}\right]---\left(22\right)$

方程组（21）-（23）中，有：

$H^{\′}=-[{\▿}_x^{2}f\left(x\right)-{\▿}_x^{2}h\left(x\right)y-{\▿}_x^{2}g\left(x\right)(z+w)]+{\▿}_{x}g\left(x\right)(W{U}^{-1}-{\mathrm{ZL}}^{-1}){\▿}_x^{T}g\left(x\right)---\left(24\right)$

F1＝-▽_xg(x)ZL^-1T（25）

F2＝-▽_xg(x)WU^-1T（26）

F3＝2ξ+T^TZL^-1T（27）

F4＝2ξ-T^TWU^-1T（28）

L′_x＝L_x+▽_xg(x)(L′_Z-L′_W)（29）

L′_α＝L_α+T^TL′_z（30）

L′_β＝L_β+T^TL′_w（31）

L′_z＝L^-1(ZL_z+L_l)（32）

L′_w＝U^-1(WL_w-L_u)（33）。

所述步骤4包括以下步骤：

步骤4-1：初始化：置迭代次数i为0，给定自变量x⁰、松弛因子l⁰和u⁰、拉格朗日乘子y⁰、z⁰和w⁰、以及类扩展变量的初值α⁰和β⁰；

步骤4-2：计算目标函数和约束的雅克比矩阵及海森矩阵，形成方程（21）-（23）形成的线性系统；

步骤4-3：对式（21）进行LU三角分解；

步骤4-4：令μ＝0，求解式（21），分别得到微分列向量Δx、Δy、Δα和Δβ各自的预测值向量Δx_af、Δy_af、Δα_af和Δβ_af，通过回代求解式（22）和（23），得到微分列向量Δl、Δu、Δz和Δw各自的预测值向量Δl_af、Δu_af、Δz_af和Δw_af；

步骤4-5：根据式（31）计算原变量仿射步长λ_afp和对偶变量仿射步长λ_afd：

$\left\{\begin{array}{c}\\ \end{array}\begin{array}{cc}{\mathrm{\&lambda;}}_{\mathrm{afp}}=\min \left(0.9995\underset{i}{\min }\right)(\frac{-l_i}{\mathrm{\&Delta;}{l}_{\mathrm{afi}}},\frac{-u_i}{\mathrm{\&Delta;}{u}_{\mathrm{afi}}}),1)& \mathrm{\&Delta;}{l}_{\mathrm{afi}}<0,\mathrm{\&Delta;}{u}_{\mathrm{afi}}<0\\ {\mathrm{\&lambda;}}_{\mathrm{afd}}=\min (0.9995\underset{i}{\min }(\frac{-z_i}{\mathrm{\&Delta;}{z}_{\mathrm{afi}}},\frac{-w_i}{\mathrm{\&Delta;}{\mathrm{\&omega;}}_{\mathrm{afi}}}),1)& \mathrm{\&Delta;}{z}_{\mathrm{afi}}<0,\mathrm{\&Delta;}{w}_{\mathrm{afi}}<0\end{array}\right.---\left(34\right)$

其中，l_i、u_i、z_i和w_i分别为列向量l、u、z和w的第i个元素，Δl_afi、Δu_afi、Δz_afi、Δw_afi分别为向量Δl_af、Δu_af、Δz_af和Δw_af的第i个元素；

更新中心参数σ和障碍参数的预测值μ_af：

$\left\{\begin{array}{c}\mathrm{\&sigma;}=\min ({\left(\frac{{\mathrm{\&rho;}}_{\mathrm{af}}}{\mathrm{\&rho;}}\right)}^3,0.1)\\ {\mathrm{\&mu;}}_{\mathrm{af}}=\mathrm{\&sigma;}\frac{\mathrm{\&rho;}}{2r}\end{array}\right.---\left(34\right)$

其中，ρ为对偶间隙，ρ_af为仿射对偶间隙，分别表示为：

ρ＝l^Tz-u^Tw（36）

ρ_af＝(l+λ_afpΔl_af)^T(z+λ_afdΔz_af)-(u+λ_afpΔu_af)^T(w+λ_afdΔw_af)（37）

步骤4-6：考虑到二次项的影响，将式（17）和（18）分别改为：

LΔz+ZΔl+ΔlΔz＝-L_l（38）

uΔw+wΔu+ΔwΔu＝-L_u（39）

为维持三角分解的因子矩阵不变，将式（38）和（39）新增的二阶项移至方程右边，则式（9）、（10）分别改为：

L_l＝LZE-μE+ΔlΔz（40）

L_u＝UWE+μE+ΔwΔu（41）；

步骤4-7：令μ＝μ_af，diag(Δl)Δz＝diag(Δl_af)Δz_af，diag(Δu)Δw＝diag(Δu_af)Δw_af，采用步骤4-3所得分解后的LU矩阵求解式（21），得到微分列向量Δx、Δy、Δα和Δβ，通过回代求解式（22）和（23），得到微分列向量Δl、Δu、Δz和Δw；按式（34）更新校正原变量仿射步长λ_afp和对偶变量仿射步长λ_afd，更新后的原变量仿射步长和对偶变量仿射步长分别用λ_p和λ_d表示；

步骤4-8：x、l、u、α、β、y、z和w分别更新为：

$\left\{\begin{array}{c}x=x+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;x}\\ l=l+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;l}\\ u=u+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;u}\\ \mathrm{\&alpha;}=\mathrm{\&alpha;}+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;\&alpha;}\\ \mathrm{\&beta;}=\mathrm{\&beta;}+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;\&beta;}\\ y=y+{\mathrm{\&lambda;}}_d\mathrm{\&Delta;y}\\ z=z+{\mathrm{\&lambda;}}_d\mathrm{\&Delta;z}\\ w=w+{\mathrm{\&lambda;}}_d\mathrm{\&Delta;w}\end{array}\right.---\left(42\right)$

步骤4-9：根据式（36）计算对偶间隙ρ，判断是否满足对偶间隙ρ小于对偶间隙给定值，且满足最大潮流偏差小于最大潮流偏差给定值，则最优潮流收敛；若不满足，转步骤4-2继续执行。

与现有技术相比，本发明的有益效果在于：

1、某些情况下，由于数值稳定性的原因使传统的内点法不能收敛，本发明能给出优化解；

2、某些情况下，当不等式约束的解空间和等式约束的解空间虽不相交，但相距不远，本发明能扩大不等式约束的解空间，使优化问题有解；

3、不等式约束的解空间扩大后，本发明能指出造成原优化问题无解的不等式约束，为进一步分析调整提供依据；

4、在原优化问题因约束条件限制无解时给出约束需放开的幅度以及对应的最优解；在某些用传统预测-校正法不能收敛的优化问题中也能有解，扩大了优化问题的解空间，提高了收敛性；

5、同时对新方法在大规模方程中的数值问题给出了实用化的处理，算法的收敛性强、解空间大、对传统内点法的改动不大，易于实现。

附图说明

图1是根据预测-校正法进行最优潮流计算流程图。

具体实施方式

下面结合附图对本发明作进一步详细说明。

图1为根据预测-校正法进行最优潮流计算流程图。

用IEEE30节点算例校验本文算法的收敛性。以网损最小为优化目标。分别采用经典的预测-校正内点法和以预测-校正法为基础的类扩展内点法计算相同算例，进行对比验证。算例中这二种方法分别称为方法1和方法2。收敛标准为ρ<1D-5和‖h(x)‖_∞<2D-4。

将电压选为扩展类，取下扩展变量α和上扩展变量β的初值为0.01，罚因子ξ为10⁶。对用同样的标幺电压约束，并不断缩小约束带，对比二种算法如表1。表中结果均采用标幺值，100MW。

两种算法在IEEE30数据计算结果比较如表1，由表1可见，当电压约束较松，预测校正法的迭代次数相同，目标函数值接近，类扩展变量接近0。当电压约束较紧，为1±0.02时，法的收敛性开始下降，迭代次数增多。此时两种算法的结果仍十分接近，类扩展变量仍接近0进一步收紧，约束带小于1±0.018后预测校正法不能收敛，类扩展变量法保持良好的收敛性带越小，类扩展变量越大，且类扩展变量与电压限值之和也越大，使目标函数减小。

在类扩展变量不接近0的算例中，突破电压限值的母线呈现明显的一致性。电压越上限2、9、29号母线；越下限最多的为26号母线，其次为23号母线。不同电压约束下越限母线固定，说明类扩展法能有效识别运行状况恶劣的变量，为采取调整和控制措施提供依据。

表1

电压约束带为1±0.01，罚因子ξ取为不同的值，类扩展法的计算结果如表2所示。由表2可见，在相当大的范围内（10⁴-10¹²），罚因子的变化对计算精度没有明显的影响。但罚因子有一较优的区域（10⁵附近），小于该区域较多，类扩展变量增大，精度下降；大于该区域较多，精度没有明显改善，迭代次数增加。

表2

电压约束带为1±0.01，罚因子ξ取10⁶，改变下、上扩展变量的初值，计算结果对比如表3。

不同扩展变量初值计算结果比较如表3；

表3

由表3可见，不同的扩展变量初值对计算结果没有明显影响，但不合适的初值将增加迭代次数。表2和表3说明，类扩展内点法计算结果对初值和罚因子不太敏感，但选择合适的参数有助于提升计算效率和精度。

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

Claims (3)

1.一种基于类扩展变量内点法的最优潮流计算方法，其特征在于：所述方法包括以下步骤：

步骤1：引入非凸非线性规划中类扩展变量，并确定引入类扩展变量后的拉格朗日函数；

步骤2：构建引入类扩展变量后的拉格朗日函数的KKT条件方程；

步骤3：求解KKT条件方程；

步骤4：根据预测-校正法进行最优潮流计算；

所述非凸非线性规划表示为：

$\left\{\begin{array}{c}\underset{x}{\min }f\left(x\right)\\ \begin{array}{cc}s.t.& h\left(x\right)=0\end{array}\\ \underset{\&OverBar;}{g}\&le;g\left(x\right)\&le;\stackrel{\&OverBar;}{g}\end{array}\right.---\left(1\right)$

其中，x为自变量构成的列向量，f(x)为优化目标函数，h(x)为等式约束矩阵，g(x)为不等式约束矩阵，和g分别为不等式约束的上限和下限列向量；

将不等式约束矩阵按类别增加类扩展变量，并在目标函数中增加带罚因子的扩展变量二次项，则式(1)变为：

$\left\{\begin{array}{c}\underset{x}{\min }\left(f\left(x\right)+\underset{k=1}{\overset{n}{\&Sigma;}}{\mathrm{\&xi;}}_k\left({\mathrm{\&alpha;}}_k^2+{\mathrm{\&beta;}}_k^2\right)\right)\\ \begin{array}{cc}s.t.& h\left(x\right)=0\end{array}\\ \underset{\&OverBar;}{g}-T\mathrm{\&alpha;}\&le;g\left(x\right)\&le;\stackrel{\&OverBar;}{g}+T\mathrm{\&beta;}\end{array}\right.---\left(2\right)$

n为设置扩展变量的不等式约束的种类，α和β为类扩展变量向量，分别代表下扩展变量和上扩展变量，维度与设置扩展变量的不等式约束的种类相等，即为n；α_k、β_k分别表示α、β的第k个元素；ξ_k为第k个类扩展变量二次项的罚因子，T为不等式约束和类扩展变量的联系矩阵，其为r×n维矩阵，形式如下：

$T=\left[\begin{array}{cccc}1& 0& ...& 0\\ 1& 0& ...& 0\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 0& 0& ...& 1\\ 0& 0& ...& 0\end{array}\right]---\left(3\right)$

若式(3)中第i个不等式约束中增加了第j个类扩展变量，则T_i,j为1，T第i行的其它元素均为0，且T中非0元均设为1；

根据类扩展变量内点法，用对数障碍函数化不等式约束为等式约束，用拉格朗日乘子法处理等式约束，与式(3)对应引入类扩展变量后的拉格朗日函数为：

$\begin{array}{c}{L}_g=f\left(x\right)+\underset{k=1}{\overset{n}{\&Sigma;}}{\mathrm{\&xi;}}_k\left({\mathrm{\&alpha;}}_k^2+{\mathrm{\&beta;}}_k^2\right)-y^{T}h\left(x\right)-\\ z^T\left(g\left(x\right)-l-\underset{\&OverBar;}{g}+T\mathrm{\&alpha;}\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\&Sigma;}}\ln l_j-\\ w^T\left(g\left(x\right)+u-\stackrel{\&OverBar;}{g}-T\mathrm{\&beta;}\right)-\mathrm{\&mu;}\underset{j=1}{\overset{r}{\&Sigma;}}\ln u_j\end{array}---\left(4\right)$

其中，L_g为拉格朗日函数；μ为障碍参数；r为不等式约束的数目；l和u均为松弛因子列向量；y、z和w分别为对应等式约束和不等式约束的拉格朗日乘子列向量；l_j和u_j分别表示向量l和u的第j个元素；其中列向量l、u、z和w的维度等于不等式约束的数目，即为r，y的维度为等式约束的数目m；

所述步骤2中，构建引入类扩展变量后的拉格朗日函数的KKT条件方程，有：

$L_x={\&dtri;}_{x}f\left(x\right)-{\&dtri;}_{x}h\left(x\right)y-{\&dtri;}_{x}g\left(x\right)(z+w)=0---\left(5\right)$

L_y＝h(x)＝0(6)

L_z＝g(x)-l-g+Tα＝0(7)

$L_w=g\left(x\right)+u-\stackrel{\&OverBar;}{g}-T\mathrm{\&beta;}=0---\left(8\right)$

L_l＝LZE-μE(9)

L_u＝UWE+μE(10)

L_α＝2ξα-T^Tz(11)

L_β＝2ξβ+T^Tw(12)

式(5)-(12)中，L_x、L_y、L_z、L_w、L_l、L_u、L_α和L_β表示拉格朗日函数对变量x、y、z、w、l、u、α和β分别取偏导数，形成的矩阵称为雅可比矩阵；其中，L＝diag(l₁,l₂,…,l_r)，Z＝diag(z₁,z₂,…,z_r)，U＝diag(u₁,u₂,…,u_r)，W＝diag(w₁,w₂,…,w_r)，ξ＝diag(ξ₁,ξ₂,…,ξ_n)，E＝[1,1,…,1]^T，E为r维向量，为对自变量x的一阶导数。

2.根据权利要求1所述的基于类扩展变量内点法的最优潮流计算方法，其特征在于：所述步骤3中，用牛顿法求解式(5)-(12)，先进行线性化处理，有：

$-\&lsqb;{\&dtri;}_x^{2}f\left(x\right)-{\&dtri;}_x^{2}h\left(x\right)y-{\&dtri;}_x^{2}g\left(x\right)(z+w)\&rsqb;\mathrm{\&Delta;}x+{\&dtri;}_{x}h\left(x\right)\mathrm{\&Delta;}y+{\&dtri;}_{x}g\left(x\right)(\mathrm{\&Delta;}z+\mathrm{\&Delta;}w)=L_x---\left(13\right)$

${\&dtri;}_x^{T}h\left(x\right)\mathrm{\&Delta;}x=-L_y---\left(14\right)$

${\&dtri;}_x^{T}g\left(x\right)\mathrm{\&Delta;}x-\mathrm{\&Delta;}l+T\mathrm{\&Delta;}\mathrm{\&alpha;}=-L_z---\left(15\right)$

${\&dtri;}_x^{T}g\left(x\right)\mathrm{\&Delta;}x+\mathrm{\&Delta;}u-T\mathrm{\&Delta;}\mathrm{\&beta;}=-L_w---\left(16\right)$

LΔz+ZΔl＝-L_l(17)

uΔw+wΔu＝-L_u(18)

2ξΔα-T^TΔz＝-L_α(19)

2ξΔβ+T^TΔw＝-L_β(20)

式(13)-(20)中，Δx、Δl、Δu、Δy、Δz、Δw、Δα和Δβ分别为x、l、u、y、z、w、α和β的微分列向量，表示x的二阶偏导数形成海森矩阵；经过推导，得到如下解耦的线性方程组：

$\left[\begin{array}{cccc}{H}^{\&prime;}& {\&dtri;}_{x}h\left(x\right)& F1& F2\\ {\&dtri;}_x^{T}h\left(x\right)& 0& 0& 0\\ -F{1}^T& 0& F3& 0\\ -F{2}^T& 0& 0& F4\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}x\\ \mathrm{\&Delta;}y\\ \mathrm{\&Delta;}\mathrm{\&alpha;}\\ \mathrm{\&Delta;}\mathrm{\&beta;}\end{array}\right]=\left[\begin{array}{c}{L}_x^{\&prime;}\\ -L_y\\ -L_{\mathrm{\&alpha;}}^{\&prime;}\\ -L_{\mathrm{\&beta;}}^{\&prime;}\end{array}\right]---\left(21\right)$

$\left[\begin{array}{cc}L& Z\\ I& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}z\\ \mathrm{\&Delta;}l\end{array}\right]=\left[\begin{array}{c}-L_1\\ -L_z^{\&prime;}-{\mathrm{ZL}}^{-1}\left({\&dtri;}_x^{T}g\right(x)\mathrm{\&Delta;}x+T\mathrm{\&Delta;}\mathrm{\&alpha;})\end{array}\right]---\left(22\right)$

$\left[\begin{array}{cc}U& W\\ I& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}w\\ \mathrm{\&Delta;}u\end{array}\right]=\left[\begin{array}{c}-L_u\\ L_w^{\&prime;}+{\mathrm{WU}}^{-1}\left({\&dtri;}_x^{T}g\right(x)\mathrm{\&Delta;}x-T\mathrm{\&Delta;}\mathrm{\&beta;})\end{array}\right]---\left(23\right)$

方程组(21)-(23)中，有：

$H^{\&prime;}=-\&lsqb;{\&dtri;}_x^{2}f\left(x\right)-{\&dtri;}_x^{2}h\left(x\right)y-{\&dtri;}_x^{2}g\left(x\right)(z+w)\&rsqb;+{\&dtri;}_{x}g\left(x\right)({\mathrm{WU}}^{-1}-{\mathrm{ZL}}^{-1}){\&dtri;}_x^{T}g\left(x\right)---\left(24\right)$

$F1=-{\&dtri;}_{x}g\left(x\right){\mathrm{ZL}}^{-1}T---\left(25\right)$

$F2=-{\&dtri;}_{x}g\left(x\right){\mathrm{WU}}^{-1}T---\left(26\right)$

F3＝2ξ+T^TZL^-1T(27)

F4＝2ξ-T^TWU^-1T(28)

$L_x^{\&prime;}=L_x+{\&dtri;}_{x}g\left(x\right)(L_Z^{\&prime;}-L_W^{\&prime;})---\left(29\right)$

L′_α＝L_α+T^TL′_z(30)

L′_β＝L_β+T^TL′_w(31)

L′_z＝L^-1(ZL_z+L_l)(32)

L′_w＝U^-1(WL_w-L_u)(33)。

3.根据权利要求2所述的基于类扩展变量内点法的最优潮流计算方法，其特征在于：所述步骤4包括以下步骤：

步骤4-1：初始化：置迭代次数i为0，给定自变量x⁰、松弛因子l⁰和u⁰、拉格朗日乘子y⁰、z⁰和w⁰、以及类扩展变量的初值α⁰和β⁰；

步骤4-2：计算目标函数和约束的雅克比矩阵及海森矩阵，形成方程(21)-(23)形成的线性系统；

步骤4-3：对式(21)进行LU三角分解；

步骤4-4：令μ＝0，求解式(21)，分别得到微分列向量Δx、Δy、Δα和Δβ各自的预测值向量Δx_af、Δy_af、Δα_af和Δβ_af，通过回代求解式(22)和(23)，得到微分列向量Δl、Δu、Δz和Δw各自的预测值向量Δl_af、Δu_af、Δz_af和Δw_af；

步骤4-5：根据式(31)计算原变量仿射步长λ_afp和对偶变量仿射步长λ_afd：

$\left\{\begin{array}{cc}{\mathrm{\&lambda;}}_{afp}=\min \left(0.9995\underset{i}{\min }\left(\frac{-l_i}{{\mathrm{\&Delta;l}}_{afi}},\frac{-u_i}{{\mathrm{\&Delta;u}}_{afi}}\right),1\right)& {\mathrm{\&Delta;l}}_{afi}<0,{\mathrm{\&Delta;u}}_{afi}<0\\ {\mathrm{\&lambda;}}_{afd}=\min \left(0.9995\underset{i}{\min }\left(\frac{-z_i}{{\mathrm{\&Delta;z}}_{afi}},\frac{-w_i}{{\mathrm{\&Delta;w}}_{afi}}\right),1\right)& {\mathrm{\&Delta;z}}_{afi}<0,{\mathrm{\&Delta;w}}_{afi}>0\end{array}\right.---\left(34\right)$

其中，l_i、u_i、z_i和w_i分别为列向量l、u、z和w的第i个元素，Δl_afi、Δu_afi、Δz_afi、Δw_afi分别为向量Δl_af、Δu_af、Δz_af和Δw_af的第i个元素；

更新中心参数σ和障碍参数的预测值μ_af：

$\left\{\begin{array}{c}\mathrm{\&sigma;}=\min ({\left(\frac{{\mathrm{\&rho;}}_{af}}{\mathrm{\&rho;}}\right)}^3,0.1)\\ {\mathrm{\&mu;}}_{af}=\mathrm{\&sigma;}\frac{\mathrm{\&rho;}}{2r}\end{array}\right.---\left(35\right)$

其中，ρ为对偶间隙，ρ_af为仿射对偶间隙，分别表示为：

ρ＝l^Tz-u^Tw(36)

ρ_af＝(l+λ_afpΔl_af)^T(z+λ_afdΔz_af)-(u+λ_afpΔu_af)^T(w+λ_afdΔw_af)(37)

步骤4-6：考虑到二次项的影响，将式(17)和(18)分别改为：

LΔz+ZΔl+ΔlΔz＝-L_l(38)

uΔw+wΔu+ΔwΔu＝-L_u(39)

为维持三角分解的因子矩阵不变，将式(38)和(39)新增的二阶项移至方程右边，则式(9)、(10)分别改为：

L_l＝LZE-μE+ΔlΔz(40)

L_u＝UWE+μE+ΔwΔu(41)；

步骤4-7：令μ＝μ_af，diag(Δl)Δz＝diag(Δl_af)Δz_af，diag(Δu)Δw＝diag(Δu_af)Δw_af，采用步骤4-3所得分解后的LU矩阵求解式(21)，得到微分列向量Δx、Δy、Δα和Δβ，通过回代求解式(22)和(23)，得到微分列向量Δl、Δu、Δz和Δw；按式(34)更新校正原变量仿射步长λ_afp和对偶变量仿射步长λ_afd，更新后的原变量仿射步长和对偶变量仿射步长分别用λ_p和λ_d表示；

步骤4-8：x、l、u、α、β、y、z和w分别更新为：

$\left\{\begin{array}{c}x=x+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;}x\\ l=l+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;}l\\ u=u+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;}u\\ \mathrm{\&alpha;}=\mathrm{\&alpha;}+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;}\mathrm{\&alpha;}\\ \mathrm{\&beta;}=\mathrm{\&beta;}+{\mathrm{\&lambda;}}_p\mathrm{\&Delta;}\mathrm{\&beta;}\\ y=y+{\mathrm{\&lambda;}}_d\mathrm{\&Delta;}y\\ z=z+{\mathrm{\&lambda;}}_d\mathrm{\&Delta;}z\\ w=w+{\mathrm{\&lambda;}}_d\mathrm{\&Delta;}w\end{array}\right.---\left(42\right)$

步骤4-9：根据式(36)计算对偶间隙ρ，判断是否同时满足对偶间隙ρ小于对偶间隙给定值，且最大潮流偏差小于最大潮流偏差给定值，则最优潮流收敛；若不满足，转步骤4-2继续执行。