基于线性规划的直角坐标形式的区间潮流计算方法与流程

本发明涉及求解电力系统区间潮流计算的技术，尤其涉及一种基于线性规划的直角坐标形式的区间潮流计算方法，该方法在涉及区间潮流方程求解时，将电压相量表示成直角坐标形式，并采用仿射算术将直角坐标下的区间潮流方程问题转化为线性规划问题的求解，通过求解线性规划来获取区间潮流解的上下边界的领域。

背景技术：

新能源机组并入电网后，其出力的不确定性使电压频繁地波动，给电网安全运行带来威胁。为分析新能源出力对电网潮流(主要是电网电压)的影响，传统的做法是采用蒙特卡洛(montecarlo)方法对新能源机组出力进行模拟，在每一个随机生成的场景下进行一次潮流计算，得到电网各节点的电压，最后将所有场景下潮流计算得到的电压进行统计，得到每个节点电压的最大和最小值，即电压的分布区间，若所有电压的分布区间都在电压上下限内，则在当前的新能源机组出力下电网是安全的，否则，应采取相应的无功电压控制策略来控制越限的节点电压，使其处于正常运行范围。随着现代电力系统的不断发展，电网规模的不断扩大，以及越来越多新能源的接入，传统的蒙特卡洛方法在抽样数目上会大大增加，潮流计算的规模也增大，计算时间急剧增长，不再适合用来分析新能源出力对电网电压的影响。因此，需要探索新的分析方法，区间潮流分析为解决这一问题提供了新的思路，但其估计的区间范围过于保守且收敛性得不保证，基于线性规划和仿射算术的区间潮流方法很好地解决这一问题，但极坐标形式下的潮流方法存在三角函数等非仿射形式的区间计算，对区间的估计精度造成了影响。因此，本发明提出了基于线性规划的直角坐标形式的区间潮流方法，避免了三角函数的仿射逼近运算，缩小了区间潮流解的范围。

区间分析计算最早由moore于1966年提出，后来在工程领域得到了广泛应用。近年来，区间分析逐渐被用于区间潮流计算，例如，文[1]将krawkzyk-moore算子用于配电网三相潮流计算，该方法既可以计算负荷确定的潮流，也可以计算含不确定性负荷的区间潮流。文[2]将krawkzyk算法用于直流区间潮流的计算，用于处理高压输电系统潮流计算中的不确定性问题，包括网络参数、机组出力和负荷的不确定性。文[3]将基于仿射算术的krawkzyk-moore算法用于电力系统潮流计算,有效缩小了区间潮流的解的范围。文[4]提出了基于线性规划的极角坐标形式的区间潮流计算方法，此方法在计算时间、收敛性能和精度都取得了突破。

参考文献：

文[1]王成山,王守相.基于区间算法的配电网三相潮流计算及算例分析[j].中国电机工程学报,2002,22(3):58-62.

文[2]郑志杰,王守相,赵龙,等.基于krawczyk算法的直流潮流区间算法[j].电力系统自动化,2012,36(20):50-53.

文[3]丁涛,崔翰韬,顾伟,等.基于区间和仿射运算的不确定潮流算法[j].电力系统自动化,2012,36(13):51-55.

文[4]vaccaroa,canizaresc,villaccid.anaffinearithmetic-basedmethodologyforreliablepowerflowanalysisinthepresenceofdatauncertainty[j].powersystems,ieeetransactionson,2010,25(2):624-632.

技术实现要素：

本发明的目的在于克服现有技术的缺点与不足，提供一种基于线性规划的直角坐标形式的区间潮流计算方法，该方法在求解区间潮流方程时，先将电压相量表示成直角坐标形式，进而将潮流方程改写成直角坐标形式。然后，将电压的实部和虚部表示成仿射形式，此仿射形式的噪声包含了各新能机组出力对电压实部或虚部产生的变化量(采用偏导数乘以新能源机组的变化量表示)。为了让初始的仿射形式下的电压实部和虚部包含实际潮流所在的区间，在电压实部和虚部的噪声元素系数上乘以一个放大系数，再将电压的实部和虚部代入直角坐标形式下的潮流方程中，构造线性规划，分别最小化和最大化电压的实部和虚部的范围，得到压缩后区间，作为最终区间算法估计的区间。由于潮流方程是直角坐标形式，避免了极坐标形式中三角函数的仿射计算，简化了仿射计算，提高了计算的精度。本发明的目的通过下述技术方案实现：

本发明的目的通过以下技术方案实现：一种基于线性规划的直角坐标形式的区间潮流计算方法，具体可以包括以下步骤：

步骤1、建立直角坐标形式的潮流方程，计算区间中点潮流解；假设风电机组的有功和无功出力分别可以表示成区间和(i∈nw)，nw代表所有风电机组组成的集合，pi为有功出力下限，为有功出力上限，qi为无功出力下限，为无功出力上限，方便本方法的阐述，不考虑风电机组运行的控制模式；建立直角坐标形式下的潮流方程，对区间中点值进行潮流计算，得到潮流解x⁽⁰⁾＝[e⁽⁰⁾,f⁽⁰⁾]，其中，e⁽⁰⁾为节点电压的实部，f⁽⁰⁾为节点电压的虚部；

步骤2、估计潮流所在的波动区间；根据潮流方程的雅可比矩阵的逆矩阵可以计算电压实部e和虚部f在x⁽⁰⁾＝[e⁽⁰⁾,f⁽⁰⁾]点对各个风电机组出力的灵敏度矩阵；即可以将电压实部e和虚部f的变化量表示成：其中x＝[e,f]不包含平衡节点的电压实部和虚部，w为风电机组的有功、无功出力以及节点电压平方构成的向量，δw为节点功率和电压平方的变化量，为潮流解对节点功率和电压平方的偏导数，w⁽⁰⁾为有功出力、无功出力以及节点电压平方的区间中点值，同时为了使估计的电压区间包含实际潮流的电压波动的范围，需要对变化量乘以一个放大系数α(大约1-10较合适)，这样，可以得到初始估计的电压波动区间为：x⁽⁰⁾为初始潮流解，δx为潮流解的变化量；

步骤3、建立电压的仿射坐标形式；由于电压主要变化主要受到风电机组的影响，因此，我们可以将风电出力看成噪声的来源，进而将电压的实部和虚部可以表示成含噪声的仿射坐标形式；

步骤4、计算节点功率的仿射坐标形式；将步骤3中得到的电压仿射坐标形式代入直角坐标潮流方程中，得到了节点功率的仿射坐标形式；

步骤5、采用线性规划压缩潮流解区间；将得到的节点功率的仿射形式转化成区间形式，构造一个线性规划，分别最小化和最大化噪声元，进而得到区间长度最小的电压实部和虚部。

在上述步骤1中，建立直角坐标形式的潮流方程，计算区间中点潮流解的步骤具体为：

1)将节点电压写成直角坐标形式，即n为系统节点个数，表示第i个节点的节点电压，ei和fi分别为实部和虚部；

对于pq节点，列写潮流方程如下：

式中，pq为pq节点集合，δpi为节点i的有功不平衡量，pli为节点的有功负荷，pi为节点j的有功出力，gij为导纳矩阵的第i行第j列元素的虚部，bij为导纳矩阵的第i行第j列元素的虚部，δqi节点i的无功不平衡量，qi为节点i的无功出力，qli为节点i的无功负荷，ei和ej为节点i和节点j的电压的实部，fi和fj分别节点i和节点j的电压的虚部；

对于pv节点，列写潮流方程如下：

式中，pv表示pv节点集合，δpi为节点i的有功不平衡量，pli为节点的有功负荷，pi为节点j的有功出力，gij为导纳矩阵的第i行第j列元素的虚部，bij为导纳矩阵的第i行第j列元素的虚部，ei和ej为节点i和节点j的电压的实部，fi和fj分别节点i和节点j的电压的虚部；

2)将风电机组的出力用区间的中点值代替，即和i∈nw，nw代表所有风电机组组成的集合，pi为有功出力下限，为有功出力上限，qi为无功出力下限，为无功出力上限；采用牛顿拉夫逊法便可以得到区间中点的潮流解x⁽⁰⁾＝[e⁽⁰⁾,f⁽⁰⁾]，e⁽⁰⁾为节点电压的实部，f⁽⁰⁾为节点电压的虚部。

所述步骤2中的估计潮流所在的波动区间，需要计算区间中点潮流解处的雅可比矩阵的逆矩阵；具有步骤如下：

1)已知区间中点潮流解为x⁽⁰⁾＝[e⁽⁰⁾,f⁽⁰⁾]，可以得到功率变化量的表达式：

式中，δs＝[δp,δq,δv²]^t为节点注入功率和节点电压平方的变化量，δp为有功功率的变化量，δq为无功功率的变化量，δv²为节点电压幅值平方的变化量，为雅可比矩阵，为节点有功注入功率对电压实部的偏导数，为有功注入功率对电压实部的偏导数，为节点电压的平方对电压实部的偏导数，为节点有功注入功率对电压虚部的偏导数，为有功注入功率对电压虚部的偏导数，为节点电压的平方对电压虚部的偏导数，δx＝[δe,δf]^t为潮流变化量，δe为节点电压实部的变化量，δf为节点电压虚部的变化量；

2)根据1)中的式子，可以推出相应的电压实部和虚部变化量的表达式如下：

式中，h＝j^-1为雅可比矩阵的逆矩阵，δs＝[δp,δq,δv²]^t为节点注入功率和节点电压平方的变化量，δp为有功功率的变化量，δq为无功功率的变化量，δv²为节点电压幅值平方的变化量，δx＝[δe,δf]^t为潮流变化量，δe为节点电压实部的变化量，δf为节点电压虚部的变化量，s⁽⁰⁾为节点注入功率和节点电压平方所在区间的中点值，同时我们可以将h写成：

式中，为电压实部对节点注入有功功率的偏导数，为电压虚部对节点有功注入功率的偏导数，为电压实部对节点注入无功功率的偏导数，为电压虚部对节点注入无功功率的偏导数，为电压实部对节点电压的平方的偏导数，为电压虚部对节点电压的平方的偏导数；若只考虑风电机组的出力发生波动，即忽略节点电压的变化，则可以将电压实部和虚部的变化量写成w＝[p,q]，p和q分别为风电机组的有功和无功出力为风电机组的有功和无功出力，δw为节点功率的变化量，为潮流解对节点功率的偏导数，w⁽⁰⁾为有功和无功出力区间中点值；

3)为了使估计的潮流区间包含实际潮流解的波动范围，在变化量上加放大系数α，放大系数的选取比较宽松，但α太大会导致估计的区间太大，导致最终获取的潮流区间的放大，α太小会导致估计的区间不能完全包含实际潮流解，一般需要根据风电波动区间的幅度选取(一般可选1-10之间)，进而可以得到估计的潮流区间x⁽⁰⁾为初始潮流解，δx为潮流解的变化量。

所述步骤3中建立电压的仿射坐标形式的步骤具体如下：

1)确定噪声元的个数；由于只有风电机组的出力具有波动性，它们是引起电压波动的根源，因此噪声元的个数包含了2k个，k为风电机组的个数；

2)形成仿射坐标。步骤2中已经计算出了电压实部与虚部变化量的表达式，因此，可以将电压实部和虚部表示成如下的仿射形式：

式中，δpj为节点j的注入有功功率变化量，δqj为节点j的注入无功功率变化量，为节点j注入有功功率区间的中点值，为节点j注入无功功率区间的中点值，为节点i的电压实部对节点j的注入有功功率的偏导数，为节点i的电压实部对节点j的注入无功功率的偏导数，为节点i的电压虚部对节点j的注入有功功率的偏导数，为节点i的电压虚部对节点j的注入无功功率的偏导数；为节点的节点注入有功功率对应的噪声元，为节点的节点注入有功功率对应的噪声元。

所述步骤4中计算节点功率的仿射坐标形式的步骤具体如下：

1)将步骤3中电压的仿射坐标形式代入到潮流方程中；根据仿射算术的则运算法则，即：

xy＝x0y0+(x1y0+x0y1)ε1+(x2y0+x0y2)ε2+…+(xny0+x0yn)εn+zkεk，

式中，εi为第i个噪声元，xi和yi为区间和对应的第i个位置的噪声系数，zk为新增加的噪声元εk对应的系数，x0和y0分别为区间和对应的区间中点值；

便可以计算出潮流方程的仿射形式：

和

式中，分别为有功功率、无功功率和电压平方的仿射形式，均为仿射计算中有功功率、无功功率和电压平方表达式出现的常数项，为仿射计算后噪声元的系数，分别为新产生噪声元的对节点注入有功功率、无功功率和电压平方对应的系数，nn为新产生噪声元的集合，nw代表所有风电机组组成的集合；

2)将常数项替换；将功率表示成仿射形式是为了计算出现功率的波动范围，而产生的常数项实际上应该是原来波动区间的中点值，因此，我们需要将相应的功率和电压仿射式中产生的常数进行替换，公式为：

式中，分别为有功功率、注入无功功率和电压区间中点值，若为常数，则区间中点为本身；分别为节点的注入有功功率、注入无功功率和电压平方。pli和qli分别为节点i的有功负荷和无功负荷。

所述步骤5中采用线性规划压缩潮流解区间的具体步骤如下：

1)将新增加的噪声元部分，按照仿射运算的逆运算，将其转化成区间形式；例如：其中，即为所有系数的绝对值的乘积；

2)将潮流方程进行排序，依次为功潮流方程、无功潮流方程和电压潮流方程，将相应的表达式表示成矩阵向量形式：

s＝ay+b，

式中，a矩阵中分别代表由组成元素形成的矩阵，为仿射计算后噪声元的系数；b向量为新增噪声元合成的区间向量，[-b^p,b^p]为有功功率合成区间向量，[-b^q,b^q]为无功功率合成区间向量，[-b^v,b^v]为节点电压合成区间向量；为有功功率对应的噪声元，为无功功率对应的噪声元；分别表示有功功率，无功功率和节点电压平方对应的区间；

3)压缩噪声元向量；对噪声元向量初始时我们将它们的限定在区间[-1,1]中，并乘以了相应的放大系数，目的是包含实际的潮流区间其中，pmin、qmin和分别为注入有功功率、无功功率以及电压平方的下限，pmax、qmax和分别为注入有功功率、无功功率以及电压平方的上限；但实际上y可以在包含潮流区间的基础上进一步被压缩，直到使ay+b＝f^sp，因此我们建立如下线性规划：

以及

式中，c＝f^sp-b，inf(·)表示求下限，sup(·)表示求上限，nw代表所有风电机组组成的集合，n为系统的节点总数，为节点的节点注入有功功率对应的噪声元，为节点的节点注入有功功率对应的噪声元，为实数，取值范围为-1到1之间，aij和aik为系数矩阵a的元素；

一共需要求解4k个单目标线性规划，才能得到所有压缩的噪声元，即：进一步，我们将其代入了电压实部和虚部的仿射形式中，得到：

式中，和分别为电压实部和虚部的仿射形式，和分别为电压实部和虚部的区间中点值，和分别为电压实部和虚部噪声元对节点注入有功功率的对应的系数，和分别为电压实部和虚部噪声元对节点注入无功功率的对应的系数，和分别为经过仿射后的新的噪声元，nw代表所有风电机组组成的集合；

采用仿射逆运算，便可以得到电压实部和虚部所在区间，即和emin为实部的下限，emax为实部的上限，fmin为虚部的下限，fmax为虚部的上限；

4)求电压幅值和相角范围；对pv节点，我们不需要求电压幅值的范围，只需要求其相角的范围。对于pq节点，电压幅值和相角的范围都需要重新确定；对电压幅值，有对于电压相角，arctan为正切三角函数的反函数，为电压幅值的区间，和分别为电压实部和虚部所在区间。

本发明的目的也可以通过以下技术方案实现：一种基于线性规划的直角坐标形式的区间潮流计算方法，考虑风电出力都在相应的区间内变化，可以包括以下步骤：

步骤1、建立直角坐标形式的潮流方程，计算区间中点潮流解。假设风电机组的有功和无功出力分别可以表示成区间和(i∈nw)，nw代表所有风电机组组成的集合，pi为有功出力下限，为有功出力上限，qi为无功出力下限，为无功出力上限，方便本方法的阐述，不考虑风电机组运行的控制模式。建立直角坐标形式下的潮流方程，对区间中点值进行潮流计算，得到潮流解x⁽⁰⁾＝[e⁽⁰⁾,f⁽⁰⁾]，其中，e⁽⁰⁾为节点电压的实部，f⁽⁰⁾为节点电压的虚部。

步骤2、估计潮流所在的波动区间。根据潮流方程的雅可比矩阵的逆矩阵可以计算电压实部e和虚部f在x⁽⁰⁾＝[e⁽⁰⁾,f⁽⁰⁾]点对各个风电机组出力的灵敏度矩阵。即可以将电压实部e和虚部f的变化量表示成：其中x＝[e,f]不包含平衡节点的电压实部和虚部，w为风电机组的有功、无功出力以及节点电压平方构成的向量，δw为节点功率和电压平方的变化量，为潮流解对节点功率和电压平方的偏导数，w⁽⁰⁾为有功出力、无功出力以及节点电压平方的区间中点值，同时为了使估计的电压区间包含实际潮流的电压波动的范围，需要对变化量乘以一个放大系数α(大约1-10较合适)，这样，可以得到初始估计的电压波动区间为：x⁽⁰⁾为初始潮流解，δx为潮流解的变化量。

步骤3、建立电压的仿射坐标形式。由于电压主要受到风电机组的影响，因此，我们可以将风电出力看成噪声的来源，进而将电压的实部和虚部可以表示成含噪声的仿射坐标形式。

步骤4、计算节点功率的仿射坐标形式。将步骤3中得到的电压仿射坐标形式代入直角坐标潮流方程中，得到了节点功率的仿射坐标形式。

步骤5、采用线性规划压缩潮流解区间。将得到的节点功率的仿射形式转化成区间形式，构造一个线性规划，分别最小化和最大化噪声元，进而得到区间长度最小的电压实部和虚部。

在上述步骤1中，建立直角坐标形式的潮流方程，计算区间中点潮流解的步骤具体为：

1)将节点电压写成直角坐标形式，即n为系统节点个数，表示第i个节点的节点电压，ei和fi分别为实部和虚部。对于pq节点，列写潮流方程如下：

式中，pq为pq节点集合，δpi为节点i的有功不平衡量，pli为节点的有功负荷，pi为节点j的有功出力，gij为导纳矩阵的第i行第j列元素的虚部，bij为导纳矩阵的第i行第j列元素的虚部，δqi节点i的无功不平衡量，qi为节点i的无功出力，qli为节点i的无功负荷，ei和ej为节点i和节点j的电压的实部，fi和fj分别节点i和节点j的电压的虚部。

对于pv节点，列写潮流方程如下：

式中，pv表示pv节点集合，δpi为节点i的有功不平衡量，pli为节点的有功负荷，pi为节点j的有功出力，gij为导纳矩阵的第i行第j列元素的虚部，bij为导纳矩阵的第i行第j列元素的虚部，ei和ej为节点i和节点j的电压的实部，fi和fj分别节点i和节点j的电压的虚部。

2)将风电机组的出力用区间的中点值代替，即和i∈nw，nw代表所有风电机组组成的集合，pi为有功出力下限，为有功出力上限，qi为无功出力下限，为无功出力上限。采用牛顿拉夫逊法便可以得到区间中点的潮流解x⁽⁰⁾＝[e⁽⁰⁾,f⁽⁰⁾]，e⁽⁰⁾为节点电压的实部，f⁽⁰⁾为节点电压的虚部。

所述步骤2中的估计潮流所在的波动区间，需要计算区间中点潮流解处的雅可比矩阵的逆矩阵。具有步骤如下：

1)已知区间中点潮流解为x⁽⁰⁾＝[e⁽⁰⁾,f⁽⁰⁾]，可以得到功率变化量的表达式：

式中，δs＝[δp,δq,δv²]^t为节点注入功率和节点电压平方的变化量，δp为有功功率的变化量，δq为无功功率的变化量，δv²为节点电压幅值平方的变化量，为雅可比矩阵，为节点有功注入功率对电压实部的偏导数，为有功注入功率对电压实部的偏导数，为节点电压的平方对电压实部的偏导数，为节点有功注入功率对电压虚部的偏导数，为有功注入功率对电压虚部的偏导数，为节点电压的平方对电压虚部的偏导数，δx＝[δe,δf]^t为潮流变化量，δe为节点电压实部的变化量，δf为节点电压虚部的变化量。

2)根据1)中的式子，可以推出相应的电压实部和虚部变化量的表达式如下：

式中，h＝j^-1为雅可比矩阵的逆矩阵，δs＝[δp,δq,δv²]^t为节点注入功率和节点电压平方的变化量，δp为有功功率的变化量，δq为无功功率的变化量，δv²为节点电压幅值平方的变化量，δx＝[δe,δf]^t为潮流变化量，δe为节点电压实部的变化量，δf为节点电压虚部的变化量，s⁽⁰⁾为节点注入功率和节点电压平方所在区间的中点值，同时我们可以将h写成：

式中，为电压实部对节点注入有功功率的偏导数，为电压虚部对节点有功注入功率的偏导数，为电压实部对节点注入无功功率的偏导数，为电压虚部对节点注入无功功率的偏导数，为电压实部对节点电压的平方的偏导数，为电压虚部对节点电压的平方的偏导数。若只考虑风电机组的出力发生波动，即忽略节点电压的变化，则可以将电压实部和虚部的变化量写成w＝[p,q]，p和q分别为风电机组的有功和无功出力，δw为节点功率的变化量，为潮流解对节点功率的偏导数，w⁽⁰⁾为有功和无功出力区间中点值。

3)为了使估计的潮流区间包含实际潮流解的波动范围，在变化量上加放大系数α，放大系数的选取比较宽松，但α太大会导致估计的区间太大，导致最终获取的潮流区间的放大，α太小会导致估计的区间不能完全包含实际潮流解，一般需要根据风电波动区间的幅度选取(一般可选1-10之间)，进而可以得到估计的潮流区间x⁽⁰⁾为初始潮流解，δx为潮流解的变化量。

所述步骤3中建立电压的仿射坐标形式的步骤具体如下：

1)确定噪声元的个数。由于只有风电机组的出力具有波动性，它们是引起电压波动的根源，因此噪声元的个数包含了2k个，k为风电机组的个数。

2)形成仿射坐标。步骤2中已经计算出了电压实部与虚部变化量的表达式，因此，可以将电压实部和虚部表示成如下的仿射形式：

式中，δpj为节点j的注入有功功率变化量，δqj为节点j的注入无功功率变化量，为节点j注入有功功率区间的中点值，为节点j注入无功功率区间的中点值，为节点i的电压实部对节点j的注入有功功率的偏导数，为节点i的电压实部对节点j的注入无功功率的偏导数，为节点i的电压虚部对节点j的注入有功功率的偏导数，为节点i的电压虚部对节点j的注入无功功率的偏导数。为节点的节点注入有功功率对应的噪声元，为节点的节点注入有功功率对应的噪声元。

所述步骤4中计算节点功率的仿射坐标形式的步骤具体如下：

1)将步骤3中电压的仿射坐标形式代入到潮流方程中。根据仿射算术的则运算法则，即：

式中，εi为第i个噪声元，xi和yi为区间和对应的第i个位置的噪声系数，zk为新增加的噪声元εk对应的系数，x0和y0分别为区间和对应的区间中点值。

便可以计算出潮流方程的仿射形式：

和

式中，分别为有功功率、无功功率和电压平方的仿射形式，均为仿射计算中有功功率、无功功率和电压平方表达式出现的常数项，为仿射计算后噪声元的系数，分别为新产生噪声元的对节点注入有功功率、无功功率和电压平方对应的系数，nn为新产生噪声元的集合，nw代表所有风电机组组成的集合。

2)将常数项替换。将功率表示成仿射形式是为了计算出现功率的波动范围，而产生的常数项实际上应该是原来波动区间的中点值，因此，我们需要将相应的功率和电压仿射式中产生的常数进行替换，公式为：

式中，分别为有功功率、注入无功功率和电压区间中点值，若为常数，则区间中点为本身；分别为节点的注入有功功率、注入无功功率和电压平方。pli和qli分别为节点i的有功负荷和无功负荷。

所述步骤5中采用线性规划压缩潮流解区间的具体步骤如下：

1)将新增加的噪声元部分，按照仿射运算的逆运算，将其转化成区间形式。例如：其中，即为所有系数的绝对值的乘积。

2)将潮流方程进行排序，依次为功潮流方程、无功潮流方程和电压潮流方程，将相应的表达式表示成矩阵向量形式：

s＝ay+b，

式中，a矩阵中分别代表由组成元素形成的矩阵，为仿射计算后噪声元的系数；b向量为新增噪声元合成的区间向量，[-b^p,b^p]为有功功率合成区间向量，[-b^q,b^q]为无功功率合成区间向量，[-b^v,b^v]为节点电压合成区间向量；为有功功率对应的噪声元，为无功功率对应的噪声元；分别表示有功功率，无功功率和节点电压平方对应的区间；

3)压缩噪声元向量。对噪声元向量初始时我们将它们的限定在区间[-1,1]中，并乘以了相应的放大系数，目的是包含实际的潮流区间其中，pmin、qmin和分别为注入有功功率、无功功率以及电压平方的下限，pmax、qmax和分别为注入有功功率、无功功率以及电压平方的上限。但实际上y可以在包含潮流区间的基础上进一步被压缩，直到使ay+b＝f^sp，因此我们建立如下线性规划：

以及

式中，c＝f^sp-b，inf(·)表示求下限，sup(·)表示求上限，nw代表所有风电机组组成的集合，n为系统的节点总数，为节点的节点注入有功功率对应的噪声元，为节点的节点注入有功功率对应的噪声元，为实数，取值范围为-1到1之间，aij和aik为系数矩阵a的元素。

一共需要求解4k个单目标线性规划，才能得到所有压缩的噪声元，即：进一步，我们将其代入了电压实部和虚部的仿射形式中，得到：

式中，和分别为电压实部和虚部的仿射形式，和分别为电压实部和虚部的区间中点値，和分别为电压实部和虚部噪声元对节点注入有功功率的对应的系数，和分别为电压实部和虚部噪声元对节点注入无功功率的对应的系数，和分别为经过仿射后的新的噪声元，nw代表所有风电机组组成的集合。

采用仿射逆运算，便可以得到电压实部和虚部所在区间，即和emin为实部的下限，emax为实部的上限，fmin为虚部的下限，fmax为虚部的上限。

4)求电压幅值和相角范围。对pv节点，我们不需要求电压幅值的范围，只需要求其相角的范围。对于pq节点，电压幅值和相角的范围都需要重新确定。对电压幅值，有对于电压相角，arctan为正切三角函数的反函数，为电压幅值的区间，和分别为电压实部和虚部所在区间。

本发明相对于现有技术具有如下的优点及效果：

(1)本发明可用于分析风电、光伏等新能源机组出力和负荷为不确定区间时的潮流，确定电压和平衡节点出力的波动区间，为调度运行工作人员提供电网运行的安全信息。

(2)本发明采用的仿射算术考虑了区间运算过程中区间相关性的处理，可以大幅缩减区间宽度。

(3)本发明在区间潮流的仿射计算中采用了直角坐标形式，避免了三角函数的计算，无需采用切比雪夫近似，简化了计算过程，可以进一步提高区间潮流的精度。

(4)本发明利用线性规划来获取区间潮流的上下边界，简化了区间计算。

附图说明

图1是修改后ieee30节点系统的接线图，图中增加了3台风电机组。表示发电机节点，表示负荷节点，表示平衡机节点，表示风电机组，表示变压器，表示电容器。

图2是直角坐标区间算法和montecarlo模拟方法得到的各节点电压实部的区间情况，纵坐标单位为标幺值(p.u.)，其中e1l和e1u为区间算法得到的电压实部的下限和上限，e2l和e2u为montecarlo得到的电压实部的下限和上限。可见，区间算法所得的电压实部的区间包含了montecarlo方法得到的电压实部的区间。

图3是直角坐标区间算法和montecarlo模拟方法得到的各节点电压虚部的区间情况，纵坐标单位为标幺值(p.u.)，其中f1l和f1u为区间算法得到的电压虚部的下限和上限，f2l和f2u为montecarlo得到的电压虚部的下限和上限。可见，区间算法所得的电压虚部的区间包含了montecarlo方法得到的电压虚部的区间。

图4是直角坐标区间算法和montecarlo模拟方法得到的各节点电压幅值的区间情况，纵坐标单位为标幺值(p.u.)，其中v1l和v1u为区间算法得到的电压幅值的下限和上限，v2l和v2u为montecarlo得到的电压幅值的下限和上限。可见，区间算法所得的电压幅值的区间包含了montecarlo方法得到的电压幅值的区间。

图5是直角坐标区间算法和montecarlo模拟方法得到的各节点电压相角的区间情况，纵坐标单位为度(deg.)，其中t1l和t1u为区间算法得到的电压相角的下限和上限，t2l和t2u为montecarlo得到的电压相角的下限和上限。可见，区间算法所得的电压相角的区间包含了montecarlo方法得到的电压相角的区间。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例

为便于理解本发明，下面结合附图进行阐述。

采用修改后的ieee30节点系统进行测试，系统的接线图见附图中的图1，该系统有24条传输线路，4台变压器，2个无功补偿点，9台发电机组，参数见表2，其中3台为风电机组，1台为平衡机组，24个负荷节点。三台风电机组的出力的区间如表1(表1为风电机组出力参数表)：

表1

所有参数的计算都采用标幺制，基准功率取100mv·a。

下面具体说明直角坐标区间潮流计算的算法步骤：

第一步，读取ieee30节点数据，包括了发电机、负荷、线路、变压器和接地电容参数。并设置风电有功出力和无功出力波动的区间范围。

第二步，利用支路追加法形成导纳矩阵。

第三步，取风电有功功率出力和无功功率出力的中点值进行潮流计算，采用牛顿拉夫逊法计算，得到中点处的潮流解为x⁽⁰⁾＝[e⁽⁰⁾,f⁽⁰⁾]。

第四步，计算雅可比矩阵的逆矩阵，并采用步骤3中的电压虚部和实部的变化量计算公式，计算得到电压实部和虚部的估计区间和仿射形式的系数并取放大系数α＝1。

第五步，采用区间的仿射算术将潮流方程转化为仿射形式，得到相关的仿射系数即可得到步骤5中的矩阵a和向量b等线性规划的相关参数，从而构造求解区间潮流的两类线性规划。

第六步，采用matlab的linprog函数求解线性规划，得到压缩后的噪声元，其代入了电压实部和虚部的仿射形式中，便可得到它们的波动区间。

第七步，根据电压的实部和虚部的变化区间，计算出电压幅值和相角的波动的区间。

为了进一步验证算法的有效性，我们采用了montecarlo法(mc)对风电出力区间进行模拟潮流计算，假设其出力在区间内服从均匀分布，并且各风电机组出力的样本相互独立，抽取了5000个样本，统计电压实部、电压虚部、电压幅值和电压相角的最大和最小值。为方便作图，对节点编号进行重排，1-5号为常规发电机节点，30号为平衡节点，其余节点按原来编号从小到大排列，并将得到的结果和直角坐标的区间潮流计算结果进行比较，得到各节点电压实部和虚部的区间分布情况如图2和图3所示，可见区间算法所得的电压实部和虚部所在的区间包含了montecarlo方法得到的电压实部和虚部所在的区间。电压幅值和相角所在的区间分布情况如图4和图5所示，区间算法所得的电压幅值和电压相角的区间包含了montecarlo方法得到的电压幅值和相角的区间。从以上分析可知，区间潮流算法所得到的潮流区间包含了montecarlo模拟方法得到的潮流区间，两者相差很小，符合区间算法的计算要求和有效性。但计算时间方面，直角坐标区间潮流算法只需要3s,而montecarlo模拟需要大约2min，在计算效率上远不及直角坐标区间潮流算法，进一步验证了它的有效性。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。