CN103779863A - 一种在极坐标系下减小区间潮流保守性的方法 - Google Patents

Abstract

本发明涉及一种在极坐标系下减小区间潮流保守性的方法，包括以下步骤：(1)获取网络元件、线路参量的区间值并用仿射区间量表示出来，获取网络拓扑结构；(2)在网络拓扑结构的基础上建立导纳矩阵；(3)形成雅可比矩阵，并用三角函数法计算仿射区间量；(4)通过迭代得到减小了保守性的潮流结果。与现有技术相比，本发明具有原理简单、易于实现，通过一次计算可以得到在某种不确定输入下的所有可能的潮流结果，并且在减小结果的保守性方面可以取得较好的效果。

Description

一种在极坐标系下减小区间潮流保守性的方法

技术领域

本发明涉及一种潮流计算技术，尤其是涉及一种在极坐标系下减小区间潮流保守性的方法。

背景技术

电力系统中存在很多不确定性，比如：系统中元件的模型参数误差；负荷需求的预测误差；风力发电的随机性。这些不确定性因素对系统的影响也是多方面的，包括：暂态稳定性、电能质量、有功和无功潮流、电压稳定、频率稳定、系统备用等。

所以，在系统参数、负荷和发电存在不确定性的情况下，寻找所有可能出现的系统运行状态的范围，比求解系统在某个瞬间的状态更加合理。从而便于在此基础上有效地对系统进行分析和评估。

区间算术是Moore在他的论文Interval Arithmetic and Automatic Error analysisin digital computing[D].Stanford：Smnford University，1962中提出的处理电子计算机浮点运算的误差问题的一种方法，随着其理论的不断成熟，开始运用到各行各业。

Moore在他的论文中定义区间数的和、差、积、商，设两个区间数

和 $Y=[\underset{\‾}{y},\stackrel{\‾}{y}],$ 则：

$X+Y=[\underset{\‾}{x}+\underset{\‾}{y},\stackrel{\‾}{x}+\stackrel{\‾}{y}]$

$X-Y=[\underset{\‾}{x}-\stackrel{\‾}{y},\stackrel{\‾}{x}-\underset{\‾}{y}]$

$X\×Y=[\min (\underset{\‾}{x}\×\underset{\‾}{y},\underset{\‾}{x}\×\stackrel{\‾}{y},\stackrel{\‾}{x}\×\underset{\‾}{y},\stackrel{\‾}{x}\×\stackrel{\‾}{y}),\max (\underset{\‾}{x}\×\underset{\‾}{y},\underset{\‾}{x}\×\stackrel{\‾}{y},\stackrel{\‾}{x}\×\underset{\‾}{y},\stackrel{\‾}{x}\×\stackrel{\‾}{y})]$

X÷Y＝X×Y^-1

此处， $Y^{-1}=[1/\stackrel{\‾}{y},1/\underset{\‾}{y}]\mathrm{if}0\∉[\underset{\‾}{y},\stackrel{\‾}{y}]$

通过上面对区间数的加减乘除的定义可以看出，区间数得加减是封闭的，而乘除是不一定封闭的。

同时定义反映区间几何特征的几个基本量，当

，区间数的中点、宽度、绝对值以及两个区间数之间的距离如下：

$\left(X\right)=\mathrm{mid}\left(X\right)=(\underset{\‾}{x}+\stackrel{\‾}{x})$

$w\left(X\right)=\stackrel{\‾}{x}-\underset{\‾}{x}$

$\left|X\right|=\max \left(\right|\underset{\‾}{x}|,|\stackrel{\‾}{x}\left|\right)$

$d(X,Y)=\max \left(\right|\underset{\‾}{x}|-\underset{\‾}{y},|\stackrel{\‾}{x}-\stackrel{\‾}{y}\left|\right)$

区间运算自身的特点使一般区间计算得到的结果都具有保守性，即存在超区间，而且这种现象在多次区间计算链中更加突出，会导致“误差爆炸”现象，甚至结果会达到毫无用处的地步。比如，当X=[-2，2]时X-X=[-4，4]而不是0，f(x)=(10-x)(10+x)，x∈X=[-2，2]用标准区间算法得到的结果是[64，144]，而其实际区间是[96，100]。

在计算不确定潮流方面，区间潮流表现出了很大的优点：

1)和随机潮流相比区间潮流不需要通过大量的数据收集来确定概率函数；

2)和模糊潮流相比不需要建立复杂的模糊隶属函数；

3)和蒙特卡罗仿真法相比它不需要进行大量的重复计算。

但是Moore R E和Kearfott R B，在Introduction to interval analysis[M].Siam，2009.中提出区间运算自身具有保守性，就是经过一系列区间运算得到的结果区间比真实的区间要大。区间计算自身的这一特性使区间潮流计算的结果也具有保守性。

王守相教授在复区间潮流保守性问题的解决方案[J].电力系统自动化，2005，29(19)：25-32中通过改变电力系统中的元件模型，和用不同的区间数学表达式来减小区间保守性，包括：方形表达式、扇形表达式、圆盘表达式。但是此方法不易于实现，且运算复杂度比较高。

Comba和Stolfi在论文Affne arithmetic and its applications to computergraphics[C].Proceedings of VI SIBGRAPI(Brazilian Symposium on ComputerGraphics and Image Processing).1993：9-18中提出的一种考虑数据之间的相关性的仿射算术，考虑到区间运算的保守性源自于忽略了数据之间的相关性，所以与仿射算术的结合在一定程度上减小了区间运算的保守性。

Comba和Stolfi提出仿射运算数

(也可以是区间数)的表达形式是：

其中x₀是仿射形式的中值；ξ_i是第i个噪声元，只能确定它在[-1，1]之间，并且噪声元之间相互独立：对应的x₁是实系数，决定了噪声元ξ_i的符号和大小。仿射算术的一个重要特征是，仿射形式的数可以共享一个噪声元，表示两个数之间依赖关系的强弱。

区间形式和仿射形式的互化：

一个区间数X，可将其表示成仿射形式为：

$X=[\underset{\‾}{x},\stackrel{\‾}{x}]=(\underset{\‾}{x}+\stackrel{\‾}{x})/2+(\underset{\‾}{x}-\stackrel{\‾}{x})/2*\mathrm{\ξ}$ 其中ξ=[-1,1]

一个给定的仿射形式：

将其转化成区间数为：

$X=[x_0-\mathrm{\ξ},x_0+\mathrm{\ξ}],\mathrm{\ξ}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left|x_i\right|$

数的仿射运算，给定两个仿射形式的数：

设三个实数α，β，γ，则

αX+βY+γ=αx₀+βy₀+γ+(x₁+y₁)ξ₁+…+(x_n+y_n)ξ_n，

这样就会产生新的噪声元，使得乘积依然是仿射形式。用仿射运算求解上面的问题可以得到X-X=(0+2ξ)-(0+2ξ)=0，对x∈X=[-2，2]，f(x)=(10-x)(10+x)的结果是[96，104]，和标准的区间算法相比大大降低了保守性。这样就可以将仿射运算的计算方法运用到区间潮流中。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种原理简单、易于实现的在极坐标系下减小区间潮流保守性的方法，通过一次计算可以得到在某种不确定输入下的所有可能的潮流结果，并且在减小结果的保守性方面可以取得较好的效果。

本发明的目的可以通过以下技术方案来实现：

一种在极坐标系下减小区间潮流保守性的方法，其特征在于，包括以下步骤：

(1)获取网络元件、线路参量的区间值并用仿射区间量表示出来，获取网络拓扑结构：

(2)在网络拓扑结构的基础上建立导纳矩阵；

(3)形成雅可比矩阵，并用三角函数法计算仿射区间量：

(4)通过迭代得到减小了保守性的潮流结果。

区间算术保守性的产生是由于区间运算的过程中认为各个点的数值是相互独立的，而没有考虑数据之间的相关性，造成最后得到的结果区间总是会比实际的区间要大一些。本专利引入仿射运算，但是因为仿射运算的计算复杂度和引入噪声元的数量具有指数相关性，所以为了均衡计算精度和计算速度，此专利中只考虑数据间的一次相关性。这样将电网中的各区间参量均转化为仿射区间变量，用仿射算术的运算法则进行计算。这从一定程度上减小了区间运算的保守性。

该方法具体过程如下：

步骤一：根据电网拓扑结构和元件参数求出导纳矩阵Y：

$Y=[\underset{\‾}{G},\stackrel{\‾}{G}]+i[\underset{\‾}{B},\stackrel{\‾}{B}]$

其中G是电阻，B是电纳，和*分别表示区间数的上限和下限，将导纳矩阵Y转化为放射区间的形式为：

$\mathrm{IF}\left(Y\right)=(\underset{\‾}{G}+\stackrel{\‾}{G})/2+(\stackrel{\‾}{G}-\underset{\‾}{G})/2+\mathrm{\ξ}[(\underset{\‾}{B}+\stackrel{\‾}{B})/2+(\stackrel{\‾}{B}-\underset{\‾}{B})/2]$

其中ξ是噪声元

步骤二：设置自变量区间初始值仿射形式为，：

$\mathrm{IF}\left(X^{\left(k\right)}\right)=\left[\begin{array}{c}({\underset{\‾}{U}}^{\left(K\right)}+{\stackrel{\‾}{U}}^{\left(k\right)})/2\\ ({\underset{\‾}{\mathrm{\θ}}}^{\left(k\right)}+{\stackrel{\‾}{\mathrm{\θ}}}^{\left(k\right)})/2\end{array}\right]+\left[\begin{array}{c}({\underset{\‾}{U}}^{\left(k\right)}-{\stackrel{\‾}{U}}^{\left(k\right)})/2\\ ({\underset{\‾}{\mathrm{\θ}}}^{\left(k\right)}+{\stackrel{\‾}{\mathrm{\θ}}}^{\left(k\right)})/2\end{array}\right]\mathrm{\ξ},$

其中k=0，U^(k)和θ^(k)是第K次迭代的电压和相角：

步骤三：求区间雅可比矩阵

$F^{\′}\left(J^{\left(k\right)}\right)=\left[\begin{array}{cc}[{\underset{\‾}{H}}^{\left(k\right)},{\stackrel{\‾}{H}}^{\left(k\right)}]& [{\underset{\‾}{N}}^{\left(k\right)},{\stackrel{\‾}{N}}^{\left(k\right)}]\\ [{\underset{\‾}{K}}^{\left(k\right)},{\stackrel{\‾}{K}}^{\left(k\right)}]& [{\underset{\‾}{L}}^{\left(k\right)},{\stackrel{\‾}{L}}^{\left(k\right)}]\end{array}\right]$

其中H^(k)，N^(k)，K^(k)，L^(k)分别是第k次迭代雅可比矩阵的元素矩阵，由于矩阵元素之间的相关性，此处结合三角函数的特性，减小由忽略元素之间的相关性引起的保守性，以矩阵H为例：

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂\mathrm{\θ}}_j}=-U_i{U}_j(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}})$

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{\∂\mathrm{\θ}}=U_i\underset{j=1,i\≠j}{\overset{n}{\mathrm{\Σ}}}{U}_j(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}})$

其中，H_ij是矩阵H中的第i行第j列元素，U是电压幅值，θ是电压相角，采用三角函数分析方法，

将故雅可比矩阵中元素H可按照以下方法进行计算：

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂\mathrm{\θ}}_j}=-U_i{U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\sin ({\mathrm{\θ}}_{\mathrm{ij}}+{\mathrm{\δ}}_{\mathrm{ij}})$

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{\∂\mathrm{\θ}}=-U_i{\underset{j=1,i\≠j}{\overset{n}{\mathrm{\Σ}}}U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\sin ({\mathrm{\θ}}_{\mathrm{ij}}+{\mathrm{\δ}}_{\mathrm{ij}})$

相应的N，K，L也可以用这种方法计算得到，分别如下：

$N_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂\mathrm{U\θ}}_j}{U}_j=-U_i{U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\cos ({\mathrm{\θ}}_{\mathrm{ij}}-{\mathrm{\δ}}_{\mathrm{ij}})$

$N_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂U}_j}{U}_i=-{U_1}^2{G}_{\mathrm{ij}}-P_i$

$K_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔQ}}_i}{{\∂\mathrm{\θ}}_j}=U_i{U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\cos ({\mathrm{\θ}}_{\mathrm{ij}}-{\mathrm{\δ}}_{\mathrm{ij}})$

$K_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔQ}}_i}{{\∂\mathrm{\θ}}_j}={U_i}^2{G}_{\mathrm{ii}}-P_i$

$L_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔQ}}_i}{{\∂U}_j}{U}_j=-U_i{U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\sin ({\mathrm{\θ}}_{\mathrm{ij}}-{\mathrm{\δ}}_{\mathrm{ij}})$

$L_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂U}_j}{U}_i={U_i}^2{B}_{\mathrm{ij}}-Q_i$

将区间值J^(k)的区间中值表示为J^(k)′：

$J^{\left(k\right)\′}=m{\left(J^{\left(k\right)}\right)}^{-1}$

步骤四：用仿射运算计算得到AF(J^(k)F′(J^(k)))；

步骤五：得到仿射形式AF(I-J^(k)F′(J^(k)))和AF(X^(k)-m(X^(k)))＝(rξ)^(k)，求得AF((I-Y^(k)F′(X^(k)))(rξ)^(k))，并转化成区间形式，其中I是区间数的单位，r是噪声元的系数，代表噪声元的大小，也表示不同噪声元之间的联系的强弱，X是要求解的未知区间数；

步骤六：按迭代算子：

K(x^(k))＝m(X^(k))-Y^(k)f(m(X^(k)))+(I-Y^(k)F(X^(k)))(X^(x)-m(X^(k)))，计算出K(X^(k))得到区间形式的X^(k+1)＝X^(k)∩K(X^(k))；

步骤七：判断是否满足终止条件，即

步骤八：输出结果。

与现有技术相比，本发明具有以下优点：

在仿射运算的基础上，又结合潮流计算在极坐标下固有的特性，将仿射运算和三角函数变换应用到区间潮流的求解中，这样进一步考虑了数据之间的相关性；原理简单，易于操作和实现，经验证在减小区间潮流保守性的问题上取得了很好的效果，使得潮流结果更具有参考意义。

附图说明

图1为本发明的流程图；

图2为九节点算例网络结构图；

图3为三种方法得到的节点电压幅值区间图；

图4为三种方法得到的节点电压相角区间图；

图5为混合算法的区间潮流和标准区间潮流对初始区间的敏感性曲线图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

实施例

如图1所示，一种在极坐标系下减小区间潮流保守性的方法，包括以下步骤：

步骤一：根据电网拓扑结构和元件参数求出导纳矩阵：

$Y=[\underset{\‾}{G},\stackrel{\‾}{G}]+i[\underset{\‾}{B},\stackrel{\‾}{B}]$

其中G是电阻，B是电纳，转化为放射区间的形式为：

$\mathrm{IF}\left(Y\right)=(\underset{\‾}{G}+\stackrel{\‾}{G})/2+(\stackrel{\‾}{G}-\underset{\‾}{G})/2+\mathrm{\ξ}[(\underset{\‾}{B}+\stackrel{\‾}{B})/2+(\stackrel{\‾}{B}-\underset{\‾}{B})/2]$

步骤二：设置自变量区间初始值仿射形式为(取k=0)：

$\mathrm{IF}\left(X^{\left(k\right)}\right)=\left[\begin{array}{c}({\underset{\‾}{U}}^{\left(K\right)}+{\stackrel{\‾}{U}}^{\left(k\right)})/2\\ ({\underset{\‾}{\mathrm{\θ}}}^{\left(k\right)}+{\stackrel{\‾}{\mathrm{\θ}}}^{\left(k\right)})/2\end{array}\right]+\left[\begin{array}{c}({\underset{\‾}{U}}^{\left(k\right)}-{\stackrel{\‾}{U}}^{\left(k\right)})/2\\ ({\underset{\‾}{\mathrm{\θ}}}^{\left(k\right)}+{\stackrel{\‾}{\mathrm{\θ}}}^{\left(k\right)})/2\end{array}\right]\mathrm{\ξ},$

U^(k)和θ^(k)是第K次迭代的电压和相角。

步骤三：求区间雅可比矩阵

$F^{\′}\left(J^{\left(k\right)}\right)=\left[\begin{array}{cc}[{\underset{\‾}{H}}^{\left(k\right)},{\stackrel{\‾}{H}}^{\left(k\right)}]& [{\underset{\‾}{N}}^{\left(k\right)},{\stackrel{\‾}{N}}^{\left(k\right)}]\\ [{\underset{\‾}{K}}^{\left(k\right)},{\stackrel{\‾}{K}}^{\left(k\right)}]& [{\underset{\‾}{L}}^{\left(k\right)},{\stackrel{\‾}{L}}^{\left(k\right)}]\end{array}\right]$

由于矩阵元素之间的相关性，此处结合三角函数得特性，尽量减小由忽略元素之间的相关性引起的保守性，以矩阵H为例：

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂\mathrm{\θ}}_j}=-U_i{U}_j(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}})$

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{\∂\mathrm{\θ}}=U_i\underset{j=1,i\≠j}{\overset{n}{\mathrm{\Σ}}}{U}_j(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}})$

中G_ijsinθ_ij-Bcosθ_ij部分sinθ_ij和cosθ_ij部分具有相关性，若通过标准区间加减法计算，可能存在保守性。例如：假设G=1，B=-1.2，θ=[0.2，0.4]，则通过区间加减法可以得到：

Gsinθ-Bcosθ=1-[0.1987,0.3894]+1.2*[0.9211，0.9804]=[1.3039，1.5655]而用三角函数分析方法：

$G\sin \mathrm{\θ}-B\cos \mathrm{\θ}=\sqrt{G^2+B^2}\sin (\mathrm{\θ}+\mathrm{ac}\tan (-\frac{B}{G}))$

$\mathrm{\θ}+\mathrm{ac}\tan (-\frac{B}{G})=\left[\mathrm{1.076,1.276}\right],\sin (\mathrm{\θ}+\mathrm{ac}\tan (-\frac{B}{G}))$ 在此区间内单调上升：

$G\sin \mathrm{\θ}-B\cos \mathrm{\θ}=\sqrt{1^2+{1.2}^2}*\sin \left(\right[\mathrm{1.076,1.276}\left]\right)=\left[\mathrm{1.3747,1.4947}\right],$ 故直接的区间加减存在保守性。故雅可比矩阵中元素H可按照以下方法进行计算：

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂\mathrm{\θ}}_j}=-U_i{U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\sin ({\mathrm{\θ}}_{\mathrm{ij}}+{\mathrm{\δ}}_{\mathrm{ij}})$

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{\∂\mathrm{\θ}}=-U_i{\underset{j=1,i\≠j}{\overset{n}{\mathrm{\Σ}}}U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\sin ({\mathrm{\θ}}_{\mathrm{ij}}+{\mathrm{\δ}}_{\mathrm{ij}})$

其中

得到的区间比较精确，但是应注意G_ij和B_ij符号变化时arctan(θ)的取值范围，和此区间内三角函数的单调性。相应的N，K，L也可以用这种方法计算得到。并得到J^(k)′：

$J^{\left(k\right)\′}=m{\left(J^{\left(k\right)}\right)}^{-1}$

步骤四：用仿射运算计算得到AF(J^(k)F′(J^(k)))；

步骤五：得到仿射形式AF(I-J^(k)F′(J^(k))和AF(X^(k)-m(X^(k)))＝(rξ)^(k)，求得AF((I-Y^(k)F′(X^(k)))(rξ)^(k))，并转化成区间形式；

步骤六：按迭代算子：

K(X^( ))=m(X^(k))-Y^(k)f(m(X^(k)))+(I-Y^(k)F′(X^(k)))(X^(k)-m(X^(k)))计算出K(X^(k))得到区间形式的X^(k+1)＝X^(k)∩K(X^(k))；

步骤七：判断是否满足终止条件，即

步骤八：输出结果。

本发明是在IEEE-9节点算例的基础上进行分析。网络拓扑结构如图2所示，系统中的具体参数如表1和表2所示，其中参数均化成以100MW为基准的标幺值，表1为IEEE-9节点系统节点参数表，表2为IEEE-9节点系统线路参数；

表1

表2

假设将风力发电机(异步发电机)通过2、3节点接入系统，平衡节点仍为常规机组，考虑风电场风速的不确定性，即考虑风电节点注入功率的不确定性。在对风电场风速-功率的分析的基础上，假设超短期风电功率预测的误差是10％，即风电场节点的功率注入在±10％内波动。不考虑线路模型的不确定性和负荷节点功率的不确定性。

用基于区间和仿射的混合运算求解非线性方程组的方法，取以标准电压为中心，区间宽度为0.1的初始电压区间，以标准电压相角为中心，区间宽度为8度的初始相角区间。求解区间交流潮流得到的结果和标准区间潮流的结果对比见表3：

表3

由表3可以看出，若考虑数据之间的相关性时，基于混合算法的潮流计算得到的结果与标准潮流计算得到的结果相比，在包含了所有蒙特卡洛仿真结果的前提下，又大大的减小了超区间的宽度。从图3和图4中可以直观的看到减小保守性的效果。

从图5两种方法对初始区间的敏感度可以看出，基于混合算法的区间潮流不仅减小了结果区间的保守性，而且也放宽了对初始区间的要求，这样在初始区间的选取方面就有了很大的自由度，减小了由于初始区间选的太大而使程序无法运行的可能性。可见此算法同时提高了程序的可行性。

Claims (2)

1.一种在极坐标系下减小区间潮流保守性的方法，其特征在于，包括以下步骤：

(1)获取网络元件、线路参量的区间值并用仿射区间量表示出来，获取网络拓扑结构；

(2)在网络拓扑结构的基础上建立导纳矩阵；

(3)形成雅可比矩阵，并用三角函数法计算仿射区间量；

(4)通过迭代得到减小了保守性的潮流结果。

2.根据权利要求1所述的一种在极坐标系下减小区间潮流保守性的方法，其特征在于，该方法具体过程如下：

步骤一：根据电网拓扑结构和元件参数求出导纳矩阵Y：

$Y=[\underset{\‾}{G},\stackrel{\‾}{G}]+i[\underset{\‾}{B},\stackrel{\‾}{B}]$

其中G是电阻，B是电纳，和*分别表示区间数的上限和下限，将导纳矩阵Y转化为放射区间的形式为：

$\mathrm{IF}\left(Y\right)=(\underset{\‾}{G}+\stackrel{\‾}{G})/2+(\stackrel{\‾}{G}-\underset{\‾}{G})/2+\mathrm{\ξ}[(\underset{\‾}{B}+\stackrel{\‾}{B})/2+(\stackrel{\‾}{B}-\underset{\‾}{B})/2]$

其中ξ是噪声元

步骤二：设置自变量区间初始值仿射形式为，：

$\mathrm{IF}\left(X^{\left(k\right)}\right)=\left[\begin{array}{c}({\underset{\‾}{U}}^{\left(K\right)}+{\stackrel{\‾}{U}}^{\left(k\right)})/2\\ ({\underset{\‾}{\mathrm{\θ}}}^{\left(k\right)}+{\stackrel{\‾}{\mathrm{\θ}}}^{\left(k\right)})/2\end{array}\right]+\left[\begin{array}{c}({\underset{\‾}{U}}^{\left(k\right)}-{\stackrel{\‾}{U}}^{\left(k\right)})/2\\ ({\underset{\‾}{\mathrm{\θ}}}^{\left(k\right)}+{\stackrel{\‾}{\mathrm{\θ}}}^{\left(k\right)})/2\end{array}\right]\mathrm{\ξ},$

其中k=0，U^(k)和θ^(k)是第K次迭代的电压和相角；

步骤三：求区间雅可比矩阵

$F^{\′}\left(J^{\left(k\right)}\right)=\left[\begin{array}{cc}[{\underset{\‾}{H}}^{\left(k\right)},{\stackrel{\‾}{H}}^{\left(k\right)}]& [{\underset{\‾}{N}}^{\left(k\right)},{\stackrel{\‾}{N}}^{\left(k\right)}]\\ [{\underset{\‾}{K}}^{\left(k\right)},{\stackrel{\‾}{K}}^{\left(k\right)}]& [{\underset{\‾}{L}}^{\left(k\right)},{\stackrel{\‾}{L}}^{\left(k\right)}]\end{array}\right]$

其中H^(k)，N^(k)，K^(k)，L^(k)分别是第k次迭代雅可比矩阵的元素矩阵，由于矩阵元素之间的相关性，此处结合三角函数的特性，减小由忽略元素之间的相关性引起的保守性，以矩阵H为例：

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂\mathrm{\θ}}_j}=-U_i{U}_j(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}})$

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{\∂\mathrm{\θ}}=U_i\underset{j=1,i\≠j}{\overset{n}{\mathrm{\Σ}}}{U}_j(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}})$

其中，H_ij是矩阵H中的第i行第j列元素，U是电压幅值，θ是电压相角，采用三角函数分析方法，

将故雅可比矩阵中元素H可按照以下方法进行计算：

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂\mathrm{\θ}}_j}=-U_i{U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\sin ({\mathrm{\θ}}_{\mathrm{ij}}+{\mathrm{\δ}}_{\mathrm{ij}})$

$H_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{\∂\mathrm{\θ}}=-U_i{\underset{j=1,i\≠j}{\overset{n}{\mathrm{\Σ}}}U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\sin ({\mathrm{\θ}}_{\mathrm{ij}}+{\mathrm{\δ}}_{\mathrm{ij}})$

相应的N，K，L也可以用这种方法计算得到，分别如下：

$N_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂\mathrm{U\θ}}_j}{U}_j=-U_i{U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\cos ({\mathrm{\θ}}_{\mathrm{ij}}-{\mathrm{\δ}}_{\mathrm{ij}})$

$N_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂U}_j}{U}_i=-{U_1}^2{G}_{\mathrm{ij}}-P_i$

$K_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔQ}}_i}{{\∂\mathrm{\θ}}_j}=U_i{U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\cos ({\mathrm{\θ}}_{\mathrm{ij}}-{\mathrm{\δ}}_{\mathrm{ij}})$

$K_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔQ}}_i}{{\∂\mathrm{\θ}}_j}={U_i}^2{G}_{\mathrm{ii}}-P_i$

$L_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔQ}}_i}{{\∂U}_j}{U}_j=-U_i{U}_j\sqrt{G_{\mathrm{ij}}^2+B_{\mathrm{ij}}^2}\sin ({\mathrm{\θ}}_{\mathrm{ij}}-{\mathrm{\δ}}_{\mathrm{ij}})$

$L_{\mathrm{ij}}=\frac{{\∂\mathrm{\ΔP}}_i}{{\∂U}_j}{U}_i={U_i}^2{B}_{\mathrm{ij}}-Q_i$

将区间值J^(k)的区间中值表示为J^(k)′：

$J^{\left(k\right)\′}=m{\left(J^{\left(k\right)}\right)}^{-1}$

步骤四：用仿射运算计算得到AF(J^(k)′F′(J^(k)))；

步骤五：得到仿射形式AF(I-J^(k)′F′(J^(k)))和AF(X^(k)-m(X^(k)))＝(rξ)^(k)，求得AF((I-Y^(k)F′(X^(k)))(rξ)^(k))，并转化成区间形式，其中I是区间数的单位，ξ是引入仿射区间后噪声元单位，r是噪声元的系数，代表噪声元的大小，也表示不同噪声元之间的联系的强弱，X是要求解的未知区间数。

步骤六：按迭代算子：

K(X^(k))=m(X^(k))-Y^(k)f(m(X^(k)))+(I-Y^(k)F′(X^(k)))(X^(k)-m(X^(k)))，计算出K(X^(k))得到区间形式的X^(k+1)=X^(k)∩K(X^(k))；

步骤七：判断是否满足终止条件，即

步骤八：输出结果。

Images