基于改进匹配追踪的智能配电网潮流雅可比矩阵估计方法与流程

本发明涉及一种雅可比矩阵估计方法。特别是涉及一种基于改进匹配追踪的智能配电网潮流雅可比矩阵估计方法。

背景技术：

大量分布式电源、电动汽车、电力电子装置的接入，给配电网的运行带来更大灵活性的同时，也对配电网的可观性和可控性提出了更高的要求。同步相量量测技术在配电网层面的发展，能够有效提高智能配电网的可观性，也为解决智能配电网的运行控制问题提供了新的思路。通过潮流雅可比矩阵求逆可以得到电压-功率灵敏度关系，而电压-功率灵敏度能够准确反映系统中电压与功率的变化关系，在配电系统的运行控制中发挥着重要的作用。采用离线潮流计算的电力系统灵敏度参数，由于存在元件参数不准确、信息更新不及时、难以追踪系统运行点及相关拓扑变化等缺陷，往往会造成较大的计算误差，极大地影响了电压控制的效果。同步相量量测技术的发展使得通过量测数据实现系统灵敏度参数的在线估计成为可能。

同步相量量测单元能够实现有功功率、无功功率、电压相角、电压幅值、系统频率等电气量的高精度同步量测，已经广泛应用于电力系统的参数辨识与运行控制等方面。通过同步相量量测装置采集的多时间段面的量测信息，采用最小二乘估计方法，能够实现潮流雅可比矩阵的精确估计。同时根据潮流方程的特点，雅可比矩阵本身具有较大的稀疏性，在参数估计问题中考虑稀疏性能够在减少量测组数需求的同时提高估计的精度。利用稀疏性估计系统参数，需要将原最小二乘估计问题转换为稀疏恢复问题，正交匹配追踪类的贪婪算法是解决稀疏恢复问题的有效手段。但是现有的正交匹配追踪类算法在求解稀疏恢复问题时，要求传感矩阵满足有限等距条件，恢复的成功率依赖于传感矩阵中各列间具有较小的相关性，但是雅可比矩阵估计问题的传感矩阵中各个节点电压变化向量间具有较强的相关性，从而造成采用一般稀疏恢复方法在求解恢复问题时，在量测组数较少的情况下无法保证恢复的成功率。因此结合雅可比矩阵的特殊性，需要提出更有针对性的改进方法，以进一步提高求解效率，降低问题对量测组数的依赖性。

技术实现要素：

本发明所要解决的技术问题是，提供一种在量测数据较少的情况下能够实现潮流雅可比矩阵精确估计的基于改进匹配追踪的智能配电网潮流雅可比矩阵估计方法。

本发明所采用的技术方案是：一种基于改进匹配追踪的智能配电网潮流雅可比矩阵估计方法，包括如下步骤：

1)获取配电网的节点数，将源节点编号为0，其他节点依次编号为1,…,i,…,n，输入网络最大度的保守估计值dmax，设定残差阈值ε及最大迭代次数m；

2)获取系统各节点同步相量量测装置的有功功率、无功功率、电压幅值和电压相角当前时刻的量测数据和c组历史量测数据，其中，c是大于1的整数；

3)将c组历史量测数据分别与当前量测值做差，得到c组有功功率、无功功率、电压幅值和电压相角的变化向量，利用节点1～节点n电压量测的变化向量构造传感矩阵，初始化潮流雅可比矩阵的行号m＝1；

4)若m不大于n，选取节点i的有功功率量测变化向量初始化残差向量，分别计算传感矩阵的第q列和第q+n列与传感矩阵所有列的相关性，得到电压相角的相关系数向量uθ和电压幅值的相关系数向量uu，此时i＝m，q＝m；若m大于n则用节点i的无功功率量测变化向量初始化残差向量，分别计算传感矩阵的第q-n列和第q列与传感矩阵所有列的相关性，得到电压相角的相关系数向量uθ和电压幅值的相关系数向量uu，此时i＝m-n，q＝m；初始化传感矩阵列号索引集合为空集，初始化迭代次数n＝1；

5)分别选取电压相角的相关系数向量uθ和电压幅值的相关系数向量uu中最大的z个数值，其中z＝dmax+1，将z个数值对应的传感矩阵中的列号索引构成集合ωn，更新传感矩阵列号索引集合；

6)利用传感矩阵列号索引集合求最小二乘解，更新残差向量；

7)若更新后的残差向量的2范数r小于残差阈值ε，则进入步骤11)；否则：如果迭代次数n＝1，分别选取电压相角的相关系数向量uθ和电压幅值的相关系数向量uu中最大的2个值，将2个值对应的传感矩阵中的索引构成集合ωn，更新列号索引集合，迭代次数n＝n+1，返回步骤6)；如果迭代次数n＝2，则迭代次数n＝n+1，进入步骤8)；

8)分别计算更新后的残差向量与传感矩阵各列的相关性，得到残差相关系数向量u，选取残差相关系数向量u中最大的2列，将2列对应的传感矩阵中的索引构成集合ωn，更新传感矩阵列号索引集合；

9)利用传感矩阵列号索引集合求最小二乘解，选取最小二乘解中绝对值最大的4项，将4项对应的传感矩阵中的索引更新集合ωn，重新构造列号索引集合，再次更新残差向量；

10)若再次更新的残差向量的2范数r小于残差阈值ε或者迭代次数超过设定的最大迭代次数m，则进入步骤11)；否则n＝n+1，返回步骤8)；

11)输出最小二乘估计结果，根据传感矩阵列号索引集合恢复出2n维向量作为雅可比矩阵第m行的估计结果，m＝m+1；若m大于2n，停止迭代，输出雅可比矩阵估计结果，否则，返回步骤4)。

步骤3)中

(1)所述的c组有功功率、无功功率、电压幅值和电压相角的变化向量表示为：

δpi[k]＝pi(k)-pi(0)、δqi[k]＝qi(k)-qi(0)、δvi[k]＝vi(k)-vi(0)和δθi[k]＝θi(k)-θi(0)，k＝1,2,…,c，pi(0)、qi(0)、θi(0)、vi(0)分别表示节点i当前时刻的有功功率、无功功率、电压相角和电压幅值的量测值；pi(k)、qi(k)、θi(k)、vi(k)分别表示节点i第k组的历史量测值；

(2)所述构造传感矩阵a如下：

式中，表示由电压相角和电压幅值量测变化向量构成的矩阵，δθi＝[δθi[1],…,δθi[c]]^t表示节点i的c组电压相角量测变化量组成的列向量，δvi＝[δvi[1],…,δvi[c]]^t表示节点i的c组电压幅值量测变化量组成的列向量；表示中的元素，ap,q表示传感矩阵a中第p行第q列的元素。

步骤4)中

(1)所述的选取节点i的有功功率量测变化向量初始化残差向量表示为：

r0＝δpi

式中，r0表示初始的残差向量，δpi＝[δpi[1],…,δpi[c]]^t表示节点i的c组有功功率量测变化量组成的列向量；

(2)用节点i的无功功率量测变化向量初始化残差向量表示为：

r0＝δqi

式中，δqi＝[δqi[1],…,δqi[c]]^t表示节点i的c组无功功率量测变化量组成的列向量。

(3)电压相角的相关系数向量uθ的计算方法为：

若1≤m≤n：

uθ＝abs[a^taq]

若n＜m≤2n：

uθ＝abs[a^taq-n]

式中，abs[·]表示取绝对值运算，a为传感矩阵，aq和aq-n分别表示矩阵a的第q和q-n列；

(4)电压幅值的相关系数向量uu的计算方法为：

若1≤m≤n：

uu＝abs[a^taq+n]

若n＜m≤2n：

uu＝abs[a^taq]

步骤5)所述的更新传感矩阵列号索引集合表示为：

λn＝λn-1∪ωn

式中，λn为第n次迭代的列号索引集合，当n＝1时，λn-1表示初始的列号索引集合，ωn表示新选入的索引构成的集合。

步骤6)包括：

(1)最小二乘解表示为：

式中，表示第n次迭代时的最小二乘解，r0为初始的残差向量，表示第n次迭代的列号索引集合λn中元素对应的传感矩阵各列构成的矩阵；

(2)更新后的残差向量表示为：

式中，rn为第n次迭代的残差向量。

步骤8)所述的残差相关系数向量u计算方式为：

u＝abs[a^trn-1]

式中，u为相关系数向量，abs[·]表示取绝对值运算，rn-1表示第n-1次迭代时的残差向量，当n＝1时，rn-1表示初始残差向量，a为传感矩阵。

步骤9)中

(1)所述的重新构造列号索引集合为

λn＝ωn

(2)所述的再次更新残差向量表示为：

式中，λn为第n次迭代的列号索引集合，ωn表示新选入的索引构成的集合，表示集合λn中元素对应的传感矩阵各列构成的矩阵，rn为第n次迭代的残差向量，r0为初始的残差向量。

步骤11)所述的根据传感矩阵列号索引集合恢复出2n维向量作为雅可比矩阵第m行的估计结果表示为：

式中，表示第n次迭代时的最小二乘解，表示2n维恢复向量，表示的第g个元素，表示列号索引集合λn中各元素所对应的向量的元素组成的向量，对于其他不在列号索引集合中的元素取值为0；矩阵第q列的模值，表示雅可比矩阵y的第m行第q个元素的估计解，其中，g＝q。

本发明的基于改进匹配追踪的智能配电网潮流雅可比矩阵估计方法，利用了矩阵的稀疏性质，采用压缩感知技术将估计问题转化成稀疏恢复问题，通过求解稀疏恢复问题，以少量的量测组数实现了雅可比矩阵的精确估计。并且针对现有稀疏恢复算法在求解雅可比矩阵估计问题的局限性，考虑电压幅值和相角量测变化向量之间的相关关系对算法进行了改进，进一步减少了对量测组数的要求，提高了估计的成功率和计算效率。

附图说明

图1是本发明基于改进匹配追踪的智能配电网潮流雅可比矩阵估计方法的流程图；

图2是ieee33节点算例图。

具体实施方式

下面结合实施例和附图对本发明的基于改进匹配追踪的智能配电网潮流雅可比矩阵估计方法做出详细说明。

如图1所示,本发明的基于改进匹配追踪的智能配电网潮流雅可比矩阵估计方法，包括如下步骤：

1)，获取配电网的节点数，将源节点编号为0，其他节点依次编号为1,…,i,…,n，输入网络最大度的保守估计值dmax，设定残差阈值ε及最大迭代次数m；

2)获取系统各节点同步相量量测装置的有功功率、无功功率、电压幅值和电压相角当前时刻的量测数据和c组历史量测数据，其中，c是大于1的整数；

所述的c组历史量测数据生成方法如下：

(1)采用如下的公式生成节点i的第k组有功功率量测数据,

式中，pi(k)表示节点i的第k组有功功率，pi(0)表示当前节点i的有功功率量测，是服从均值为0正态分布的随机数，分别用来模拟不同量测时刻相对于当前时刻的功率变化和量测误差。

(2)采用如下公式生成节点i的第k组无功量测数据,

式中，qi(k)表示节点i的第k组无功功率，qi(0)表示当前节点i的无功功率量测，是服从均值为0正态分布的随机数。

(3)在得到节点i的第k组有功功率和无功功率量测数据后，通过潮流计算求得对应的电压相角θi(k)和幅值vi(k)作为节点i的第k组电压相角和幅值量测值。

3)将c组历史量测数据分别与当前量测值做差，得到c组有功功率、无功功率、电压幅值和电压相角的变化向量，利用节点1～节点n电压量测的变化向量构造传感矩阵，初始化潮流雅可比矩阵的行号m＝1；其中,

所述的c组有功功率、无功功率、电压幅值和电压相角的变化向量表示为：

δpi[k]＝pi(k)-pi(0)、δqi[k]＝qi(k)-qi(0)、δvi[k]＝vi(k)-vi(0)和δθi[k]＝θi(k)-θi(0)，k＝1,2,…,c，pi(0)、qi(0)、θi(0)、vi(0)分别表示节点i当前时刻的有功功率、无功功率、电压相角和电压幅值的量测值；pi(k)、qi(k)、θi(k)、vi(k)分别表示节点i第k组的历史量测值；

(2)所述构造传感矩阵a如下：

式中，表示由电压相角和电压幅值量测变化向量构成的矩阵，δθi＝[δθi[1],…,δθi[c]]^t表示节点i的c组电压相角量测变化量组成的列向量，δvi＝[δvi[1],…,δvi[c]]^t表示节点i的c组电压幅值量测变化量组成的列向量；表示中的元素，ap,q表示传感矩阵a中第p行第q列的元素。

4)若m不大于n，选取节点i的有功功率量测变化向量初始化残差向量，分别计算传感矩阵的第q列和第q+n列与传感矩阵所有列的相关性，得到电压相角的相关系数向量uθ和电压幅值的相关系数向量uu，此时i＝m，q＝m；若m大于n则用节点i的无功功率量测变化向量初始化残差向量，分别计算传感矩阵的第q-n列和第q列与传感矩阵所有列的相关性，得到电压相角的相关系数向量uθ和电压幅值的相关系数向量uu，此时i＝m-n，q＝m；初始化传感矩阵列号索引集合为空集，初始化迭代次数n＝1；其中，

(1)所述的选取节点i的有功功率量测变化向量初始化残差向量表示为：

r0＝δpi

式中，r0表示初始的残差向量，δpi＝[δpi[1],…,δpi[c]]^t表示节点i的c组有功功率量测变化量组成的列向量；

(2)用节点i的无功功率量测变化向量初始化残差向量表示为：

r0＝δqi

式中，δqi＝[δqi[1],…,δqi[c]]^t表示节点i的c组无功功率量测变化量组成的列向量。

(3)电压相角的相关系数向量uθ的计算方法为：

若1≤m≤n：

uθ＝abs[a^taq]

若n＜m≤2n：

uθ＝abs[a^taq-n]

式中，abs[·]表示取绝对值运算，a为传感矩阵，aq和aq-n分别表示矩阵a的第q和q-n列；

(4)电压幅值的相关系数向量uu的计算方法为：

若1≤m≤n：

uu＝abs[a^taq+n]

若n＜m≤2n：

uu＝abs[a^taq]

5)分别选取电压相角的相关系数向量uθ和电压幅值的相关系数向量uu中最大的z个数值，其中z＝dmax+1，将z个数值对应的传感矩阵中的列号索引构成集合ωn，更新传感矩阵列号索引集合；

所述的更新传感矩阵列号索引集合表示为：

λn＝λn-1∪ωn

式中，λn为第n次迭代的列号索引集合，当n＝1时，λn-1表示初始的列号索引集合，ωn表示新选入的索引构成的集合。

6)利用传感矩阵列号索引集合求最小二乘解，更新残差向量；包括：

(1)最小二乘解表示为：

式中，表示第n次迭代时的最小二乘解，r0为初始的残差向量，表示第n次迭代的列号索引集合λn中元素对应的传感矩阵各列构成的矩阵；

(2)更新后的残差向量表示为：

式中，rn为第n次迭代的残差向量。

7)若更新后的残差向量的2范数r小于残差阈值ε，则进入步骤11)；否则：如果迭代次数n＝1，分别选取电压相角的相关系数向量uθ和电压幅值的相关系数向量uu中最大的2个值，将2个值对应的传感矩阵中的索引构成集合ωn，更新列号索引集合，迭代次数n＝n+1，返回步骤6)；如果迭代次数n＝2，则迭代次数n＝n+1，进入步骤8)；

8)分别计算更新后的残差向量与传感矩阵各列的相关性，得到残差相关系数向量u，选取残差相关系数向量u中最大的2列，将2列对应的传感矩阵中的索引构成集合ωn，更新传感矩阵列号索引集合；

所述的残差相关系数向量u计算方式为：

u＝abs[a^trn-1]

式中，u为相关系数向量，abs[·]表示取绝对值运算，rn-1表示第n-1次迭代时的残差向量，当n＝1时，rn-1表示初始残差向量，a为传感矩阵。

9)利用传感矩阵列号索引集合求最小二乘解，选取最小二乘解中绝对值最大的4项，将4项对应的传感矩阵中的索引更新集合ωn，重新构造列号索引集合，再次更新残差向量；其中,

所述的重新构造列号索引集合为

λn＝ωn

所述的再次更新残差向量表示为：

式中，λn为第n次迭代的列号索引集合，ωn表示新选入的索引构成的集合，表示集合λn中元素对应的传感矩阵各列构成的矩阵，rn为第n次迭代的残差向量，r0为初始的残差向量。

10)若再次更新的残差向量的2范数r小于残差阈值ε或者迭代次数超过设定的最大迭代次数m，则进入步骤11)；否则n＝n+1，返回步骤8)；

11)输出最小二乘估计结果，根据传感矩阵列号索引集合恢复出2n维向量作为雅可比矩阵第m行的估计结果，m＝m+1；若m大于2n，停止迭代，输出雅可比矩阵估计结果，否则，返回步骤4)。

所述的根据传感矩阵列号索引集合恢复出2n维向量作为雅可比矩阵第m行的估计结果表示为：

式中，表示第n次迭代时的最小二乘解，表示2n维恢复向量，表示的第g个元素，表示列号索引集合λn中各元素所对应的向量的元素组成的向量，对于其他不在列号索引集合中的元素取值为0；矩阵第q列的模值，表示雅可比矩阵y的第m行第q个元素的估计解，其中，g＝q。

下面给出具体实例：

首先输入ieee33节点算例网络拓扑连接关系如图2所示，其中节点0为平衡节点，其他节点1～32为pq节点，系统的基准容量为1mva，基准电压为12.66kv，各个pq节点的当前功率量测如表1所示。输入网络的最大度的估计值为4，模拟量测功率变化和误差随机数的标准差分别设为0.01和0.025％。采用下式计算雅可比矩阵的误差。

式中，分别表示雅可比矩阵参数的估计值，ji,j为采用精确线路参数的计算值。

为验证本发明方法的先进性，采取如下两种场景进行分析：

场景1，设置量测组数为30、35、40、45、50、55、60，分别采用最小二乘方法、正交匹配追踪算法和本专利所提出算法进行雅可比矩阵估计；

场景2，设置量测组数为40、45、50、55、60，采用本发明的方法对雅可比矩阵进行估计；设置量测组数为100、500、1000、1500、2000，直接采用最小二乘方法进行估计。

场景1估计结果如表2所示，场景2估计结果如表3和表4所示。

执行优化计算的计算机硬件环境为intel(r)xeon(r)cpue5-1620，主频为3.70ghz，内存为32gb；软件环境为windows7操作系统，采用matlab的matpower工具包计算潮流。

从表2中可以看出，当量测组数少于60组时最小二乘计算方法无法实现雅可比矩阵任意一行的估计，正交匹配追踪算法能够利用雅可比矩阵的稀疏性，实现雅可比矩阵中一部分行的估计，而且随着量测组数的增加，能够成功估计的行数增加，但是即使量测组数增加为60组时，仍然不能实现所有行的估计；本发明的方法在量测组数增加到40时，即可实现雅可比矩阵所有行的估计，减少了估计问题对量测组数的依赖，提高了估计的成功率。通过表3和表4的对比可以看出，本发明所提出的估计方法在减少对量测组数依赖的同时也能够保证估计的精确性，当量测组数为60组时的估计精度与采用最小二乘方法2000组量测的估计精度基本相同。

表1ieee33节点算例pq节点当前功率量测值

表2场景1的估计结果

表3场景2的估计结果

表4场景2的估计结果