CN104009472B - 一种基于协同粒子群的配电网状态估计方法 - Google Patents

Abstract

本发明提供一种基于协同粒子群的配电网状态估计方法，包括以下步骤：步骤l：载入配网系统及其量测配置；步骤2：初始化主种群M和协同种群集合C；步骤3：判断是否满足协同条件，若满足执行步骤4，否则执行步骤6；步骤4：对主种群M中粒子进行求解；步骤5：判断是否满足收敛条件，若满足则状态估计结束，否则转步骤3；步骤6：对协同种群集合C进行求解。本发明以节点电压的实部和虚部为状态变量，计算支路电流、支路首端复功率和节点注入复功率，运用最小二乘方法计算目标函数以衡量计算量和测量量之间的接近程度，基于粒子问的交互更新状态变量，基于种群问的协同提升算法寻优能力；实验结果表明该方法具有良好的适应性和收敛性。

Description

一种基于协同粒子群的配电网状态估计方法

技术领域

本发明涉及一种估计方法，具体涉及一种基于协同粒子群的配电网状态估计方法。

背景技术

配电网状态估计(DistributionStateEstimation，DSE)是一种利用测量数据的相关性和冗余度，应用计算机技术，采用数学处理的方法来对运行参数进行预测、拟合、纠错处理，以提高数据的可靠性与完整性，从而有效地获得配电网实时状态信息。DSE是配电网各种分析仿真高级应用的数据基础。

输电网状态估计有较为丰富的理论基础，但配电网具有不同于输电网的显著特点。配电网实时量测配置少，数据冗余度不足；网络负荷较重、伪量测精度较低、使用了大量的控制设备，使得算法收敛困难；因此需要对配电网状态估计进行深入研究。

配电网状态估计的研究主要集中在估计精度和收敛性能领域。如基于等效功率变换提升数据吸纳能力，运用模糊匹配方法提升计算速度和鲁棒性，基于测量不确定度方法提升处理不良数据的能力；它们从不同角度提高了配电网状态估计性能，但未提及量测配置和收敛性能。在考虑量测配置的研究中，有的分析了导纳矩阵与收敛性之间的关系，但没有说明量测误差对收敛性的影响；有的提出基于网络分解的区域解耦配电网状态估计算法，详细分析了量测配置的特性，但实验对象比较简单；有的提出一种适用于支路电流幅值量测为主的配网状态估计方法，对网络可观测性、算法性能进行了探讨，但没有对支路电流或功率量测覆盖度指标进行实验分析。由此可见，开展配电网状态估计研究需要对大量伪量测环境下算法的收敛性和准确性进行深入研究。

粒子群优化算法(ParticleSwarmOptimization，PSO)，以其突出的准确性、可靠性、鲁棒性、效率、多样性和相关性，是实现配电网状态估计比较好的求解方法。但配电网结构复杂、设备多样、数据量大，容易导致PSO陷入“维数灾”的困境。当前的粒子群状态估计算法，有的未能详细说明状态估计与粒子群优化的混合算法模型，有的严格控制了状态变量的维数且实验对象网络规模较小，有的没有提出提升PSO全局寻优能力的办法。因此，对于配电网状态估计粒子群算法，仍需对其收敛性、寻优能力进行深入的研究。

发明内容

为了克服上述现有技术的不足，本发明提供一种基于协同粒子群的配电网状态估计方法，以节点电压的实部和虚部为状态变量，计算支路电流、支路首端复功率和节点注入复功率，运用最小二乘方法计算目标函数以衡量计算量和测量量之间的接近程度，基于粒子间的交互来更新状态变量，基于种群间的协同提升算法寻优能力；实验结果表明该方法法具有良好的精度、适应性和收敛性。

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

提供一种基于协同粒子群的配电网状态估计方法，所述方法包括以下步骤：

步骤1：载入配网系统及其量测配置；

步骤2：初始化主种群M和协同种群集合C；

步骤3：判断是否满足协同条件，若满足执行步骤4，否则执行步骤6；

步骤4：对主种群M中粒子进行求解；

步骤5：判断是否满足收敛条件，若满足则状态估计结束，否则转步骤3；

步骤6：对协同种群集合C进行求解。

所述步骤1中，进行拓扑分析，确定搜索方向，进行节点编号，生成节点导纳矩阵Y与节点-支路关联矩阵A。

所述步骤2中，基于量测的节点注入复功率与平衡节点电压进行潮流分析算出节点电压向量V，以之作为初始化参考值；初始化主种群M，包括粒子当前位置、个体最优位置和全局最优位置，并计算其对应适应度，初始速度V(0)＝0；

适应度函数F(x)表示为：

$F\left(x\right)=f\left(x\right)+k_f\underset{k=1}{\overset{N_{\mathrm{ue}}}{\mathrm{\Σ}}}\left(\max \right[0,g_k^j\left(x\right)\left]\right)---\left(1\right)$

其中，f(x)为配电网状态估计函数；k_f为惩罚因子；N_ue为不等式约束数；为惩罚函数，分为以下几种情况：

1)负荷有功功率对应的惩罚函数和表达式为：

$\left\{\begin{array}{c}{g}_k^1\left(x\right)=P_{\mathrm{LD},\min }^l-P_{\mathrm{LD}}^l\\ g_k^2\left(x\right)=P_{\mathrm{LD}}^l-P_{\mathrm{LD},\max }^l\end{array}\right.---\left(2\right)$

其中，为第l个负荷的有功功率，和分别为第l个负荷的有功功率的下限和上限；

2)负荷无功功率对应的惩罚函数和表达式为：

$\left\{\begin{array}{c}{g}_k^3\left(x\right)=Q_{\mathrm{LD},\min }^l-Q_{\mathrm{LD}}^l\\ g_k^4\left(x\right)=Q_{\mathrm{LD}}^l-Q_{\mathrm{LD},\max }^l\end{array}\right.---\left(3\right)$

其中，为第l个负荷的无功功率，分别为第l个负荷的无功功率的下限和上限；

3)支路电流幅值对应的惩罚函数表达式为：

$g_k^5\left(x\right)=I_{\mathrm{LN}}^b-I_{\mathrm{LN},\max }^b---\left(6\right)$

其中，为第b条支路的支路电流幅值，为第b条支路的支路电流幅值的最大值；

4)节点电压幅值对应的惩罚函数和表达式为：

$\left\{\begin{array}{c}{g}_k^6\left(x\right)=V_{\mathrm{ND},\min }^n-V_{\mathrm{ND}}^n\\ g_k^7\left(x\right)=V_{\mathrm{ND}}^n-V_{\mathrm{ND},\max }^n\end{array}\right.---\left(7\right)$

其中，为第n个节点的节点电压幅值，分别为第n个节点的节点电压幅值的下限和上限；

第l个负荷的有功功率和无功功率第b条支路电流幅值第n个节点电压幅值均为状态变量x的函数，通过下式计算得到；

$\left\{\begin{array}{c}{P}_{\mathrm{LD}}+{\mathrm{iQ}}_{\mathrm{LD}}=S_n\\ I_{\mathrm{LN}}=\left|I\right|\\ V_{\mathrm{ND}}=\left|V\right|\end{array}\right.---\left(8\right)$

其中，负荷有功功率向量其中n_n为节点总数；负荷无功功率向量 $Q_{\mathrm{LD}}={[Q_{\mathrm{LD},}^1{Q}_{\mathrm{LD}}^2,\·\·\·,Q_{\mathrm{LD}}^l,\·\·\·,Q_{\mathrm{LD}}^{n_n}]}^T;$ 支路电流幅值向量 $I_{\mathrm{LN}}={[I_{\mathrm{LN}}^1,I_{\mathrm{LN}}^2,\·\·\·,I_{\mathrm{LN}}^b,\·\·\·,I_{\mathrm{LN}}^{b_n}]}^T,$ 其中b_n为支路总数；I为支路电流向量；节点电压幅值向量V为节点电压向量；

因状态变量之间无耦合关系，协同种群集合C中种群数量n_x′等于搜索空间的维度数n_x，即种群C_j仅对一维状态变量进行寻优，j＝1,2,...,n′_x；采用与主种群M初始化相同的方法对协同种群集合C进行初始化。

所述步骤3中，协同条件为：mod(it,ci)＝＝ci，其中it为当前迭代次数，ci为协同间隔次数，mod为求取余数函数，符号＝＝表示等于。

所述步骤4中，对主种群M进行求解的过程包括更新粒子的个体最优位置y_i(t+1)和全局最优位置计算粒子t+1时刻的速度v_i(t+1)以及粒子t+1时刻的位置x_i(t+1)；

1)粒子的个体最优位置是第i个粒子从开始到现在达到过的最佳位置，t时刻粒子的个体最优位置y_i(t)表示为：

$y_i\left(t\right)={[y_{i1}\left(t\right),y_{i2}\left(t\right),\·\·\·,y_{\mathrm{ij}}\left(t\right),\·\·\·,{y_{\mathrm{in}}}_x]}^T---\left(9\right)$

其中，y_ij(t)为第i个粒子t时刻在第j维度上的个体最优位置；

于是t+1时刻粒子的个体最优位置y_i(t+1)由下式计算：

$y_i(t+1)=\left\{\begin{array}{ccc}{y}_i\left(t\right),& \mathrm{if}& F\left(x_i(t+1)\right)\&GreaterEqual;F\left(y_i\left(t\right)\right)\\ x_i(t+1),& \mathrm{if}& F\left(x_i(t+1)\right)<F\left(y_i\left(t\right)\right)\end{array}\right.---\left(10\right)$

其中，x_i(t+1)为第i个粒子t+1时刻在搜索空间的位置，F(x_i(t+1))为第i个粒子t时刻在搜索空间的位置对应的适应度函数，F(y_i(t))为t时刻粒子的个体最优位置对应的适应度函数，其取值为：

$F\left(y_i\left(t\right)\right)=\min \{F\left(y_1\left(t\right)\right),F\left(y_2\left(t\right)\right),...,F\left({y_n}_s\left(t\right)\right)\}---\left(11\right)$

其中，n_s为单个种群中粒子的数量；

2)t+1时刻的全局最优位置表示为：

$\hat{y}(t+1)=y_i\left(t\right)---\left(12\right)$

在t时刻，对于任意第i个粒子t时刻的全局最优位置有：

${\hat{y}}_i\left(t\right)=[{\hat{y}}_{i1}\left(t\right),{\hat{y}}_{i2}\left(t\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\hat{y}}_{\mathrm{ij}}\left(t\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\widehat{y_{\mathrm{in}}}}_x{]}^T---\left(13\right)$

其中，为第i个粒子在第j维度上的全局最优位置；

3)计算粒子t+1时刻的速度v_i(t+1)；

在t时刻，粒子t时刻的速度其中v_ij(t)为第i个粒子第j维度上的速度，于是t+1时刻第i个粒子在第j维度上的速度v_ij(t+1)表示为：

$v_{\mathrm{ij}}(t+1)={\mathrm{wv}}_{\mathrm{ij}}\left(t\right)+c_1{r}_{1j}\left(t\right)[y_{\mathrm{ij}}\left(t\right)-x_{\mathrm{ij}}\left(t\right)]+c_2{r}_{2j}\left(t\right)[{\hat{y}}_{\mathrm{ij}}\left(t\right)-x_{\mathrm{ij}}\left(t\right)]---\left(14\right)$

其中，w为惯性因子，衡量t时刻的速度对t+1时刻移动的影响；c₁和c₂分别为认知因子和社会因子，是正数的加速度常量，分别用来度量认知成分和社会成分对于速度更新的贡献；r_1j(t)、r_2j(t)～U(0，1)都是在区间[0，1]中均匀抽取的随机数，代表不确定性因素；x_ij(t)为第i个粒子t时刻在第j维度上在搜索空间的位置，j＝1,2,…,n_x，n_x为搜索空间的维度数；

4)粒子t+1时刻的位置x_i(t+1)表示为：

x_i(t+1)＝x_i(t)+v_i(t+1)(15)

其中，x_i(t)第i个粒子t时刻在搜索空间的位置，i＝1,2,…,n_s，其中n_s为单个种群中粒子的数量；x_i(t)表示为：

$x_i\left(t\right)={[x_{i1}\left(t\right),x_{i2}\left(t\right),\&CenterDot;\&CenterDot;\&CenterDot;,x_{\mathrm{ij}}\left(t\right),\&CenterDot;\&CenterDot;\&CenterDot;,x_{{\mathrm{in}}_x}]}^T---\left(16\right)$

其中，t代表的是离散时间点，也代表迭代次数；初始位置x_i(0)～U(x_min,x_max)，x_min和x_max分别为粒子位置的最小值向量和最大值向量。

所述步骤5中，收敛条件为迭代次数超过最大次数限制，或目标函数值小于给定的阈值，即：

$\mathrm{it}>\mathrm{it}\_\lim \left|\right|F(M.\hat{y})<\mathrm{\&epsiv;}---\left(17\right)$

其中：it为当前迭代次数，it_lim为迭代次数限制，ε为收敛标准，符号||表示逻辑关系或，为主种群M中全局最优位置对应的适应度值。

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

步骤1：更新协同种群集合C；

对种群C_j，某个粒子i，式中i～U(1，n_s/2)且i∈N，其中，C_j.x_i为种群C_j中第i粒子在搜索空间的位置，为主种群M全局最优位置中第j维度上的数值，n_s为单个种群中粒子的数量，N为自然数；

步骤2：定义向量 为主种群M中全局最优位置，b_j为向量b中第j维度上的数值，种群C_j状态变量定义为：其中n_x′为协同种群集合C中种群数量；并按照步骤4的方法对每个种群C_j进行求解，计算完成后为种群C_j的全局最优位置；

步骤3：更新主种群M：M.x_i＝b，其中M.x_i为主种群M中第i个粒子在搜索空间的位置，i～U(1,n_s/2)且k∈N。

配电网状态估计以f(x)为最小为目标函数，具体有：

$\min f\left(x\right)=\underset{m=1}{\overset{n_m}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_m{[z_m-h_m\left(x\right)]}^2---\left(18\right)$

其中，n_m为量测量的数量，ω_m第m个量测量的测量权重，z_m为第m个量测量，h_m为第m个量测量的状态函数；

1)所述目标函数的不等式约束包括负荷有功功率约束、负荷无功功率约束、支路电流约束和节点电压约束，分别表示为：

$\left\{\begin{array}{c}{P}_{\mathrm{LD},\min }^l\&le;P_{\mathrm{LD}}^l\&le;P_{\mathrm{LD},\max }^l\\ Q_{\mathrm{LD},\min }^l\&le;Q_{\mathrm{LD}}^l\&le;Q_{\mathrm{LD},\max }^l\\ I_{\mathrm{LN}}^b\&le;I_{\mathrm{LN},\max }^b\\ V_{\mathrm{ND},\min }^n\&le;V_{\mathrm{ND}}^n\&le;V_{\mathrm{ND},\max }^n\end{array}\right.---\left(19\right)$

2)所述目标函数的等式约束为潮流约束，基于状态变量运用潮流相关算法求取测量量对应的电气量，其值满足潮流约束；选择节点电压为状态变量，状态方程为：

$\left\{\begin{array}{c}V=V\\ I=Y_b{A}^{T}V\\ S_b=V\left(N_f\right).*\stackrel{\&OverBar;}{I}\\ S_n=V.*\stackrel{\&OverBar;}{(Y*V)}\end{array}\right.---\left(20\right)$

其中，V为n_n×1维节点电压向量，n_n为节点总数，搜索空间的维度数n_x＝2*n_n；I为b_n×1维支路电流向量，b_n为支路总数；Y_b为支路导纳矩阵；A为节点-支路关联矩阵；S_b为b_n×1维支路首端复功率向量，N_f为支路首端节点编号向量，V(N_f)为b_n×1维支路首端电压向量，符号.*表示矩阵对应元素相乘；S_n为n_n×1维节点注入复功率，Y为节点导纳矩阵。

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

1.该方法简单、有效，精度能够满足应用需求，易于实现；在网络非线性程度较强的情况下，计算速度超过经典数学方法；

2.实用性好，可以在缺少实时量测的情况下，充分利用精度较低的伪量测数据，得到较为理想的计算结果；

3.收敛性好，能够在可接受的时间内处理大规模复杂配电网的状态估计问题。

附图说明

图1是基于协同粒子群的配电网状态估计方法流程图；

图2是精确分析相对误差概率分布曲线图；

图3是负荷倍率与迭代次数之间的曲线图。

具体实施方式

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

本发明提供一种基于协同粒子群的配电网状态估计方法，协同粒子群算法(CooperativeParticleSwarmOptimization,CPSO)最初由FransvandenBergh提出，协同粒子群在搜索质量、鲁棒性等性能方面有了重要的性能提升。提出的基于协同粒子群的配电网状态估计方法步骤如下：

步骤1：载入配网系统及其量测配置；

步骤2：初始化主种群M和协同种群集合C；

步骤3：判断是否满足协同条件，若满足执行步骤4，否则执行步骤6；

步骤4：对主种群M中粒子进行求解；

步骤5：判断是否满足收敛条件，若满足则状态估计结束，否则转步骤3；

步骤6：对协同种群集合C进行求解。

所述步骤1中，进行拓扑分析，确定搜索方向，进行节点编号，生成节点导纳矩阵Y与节点-支路关联矩阵A。

所述步骤2中，基于量测的节点注入复功率与平衡节点电压进行潮流分析算出节点电压向量V，以之作为初始化参考值；初始化主种群M，包括粒子当前位置、个体最优位置和全局最优位置，并计算其对应适应度，初始速度V(0)＝0；

适应度函数F(x)表示为：

$F\left(x\right)=f\left(x\right)+k_f\underset{k=1}{\overset{N_{\mathrm{ue}}}{\mathrm{\&Sigma;}}}\left(\max \right[0,g_k^j\left(x\right)\left]\right)---\left(1\right)$

其中，f(x)为配电网状态估计函数；k_f为惩罚因子；N_ue为不等式约束数；为惩罚函数，分为以下几种情况：

1)负荷有功功率对应的惩罚函数和表达式为：

$\left\{\begin{array}{c}{g}_k^1\left(x\right)=P_{\mathrm{LD},\min }^l-P_{\mathrm{LD}}^l\\ g_k^2\left(x\right)=P_{\mathrm{LD}}^l-P_{\mathrm{LD},\max }^l\end{array}\right.---\left(2\right)$

其中，为第l个负荷的有功功率，和分别为第l个负荷的有功功率的下限和上限；

2)负荷无功功率对应的惩罚函数和表达式为：

$\left\{\begin{array}{c}{g}_k^3\left(x\right)=Q_{\mathrm{LD},\min }^l-Q_{\mathrm{LD}}^l\\ g_k^4\left(x\right)=Q_{\mathrm{LD}}^l-Q_{\mathrm{LD},\max }^l\end{array}\right.---\left(3\right)$

其中，为第l个负荷的无功功率，分别为第l个负荷的无功功率的下限和上限；

3)支路电流幅值对应的惩罚函数表达式为：

$g_k^5\left(x\right)=I_{\mathrm{LN}}^b-I_{\mathrm{LN},\max }^b---\left(6\right)$

其中，为第b条支路的支路电流幅值，为第b条支路的支路电流幅值的最大值；

4)节点电压幅值对应的惩罚函数和表达式为：

$\left\{\begin{array}{c}{g}_k^6\left(x\right)=V_{\mathrm{ND},\min }^n-V_{\mathrm{ND}}^n\\ g_k^7\left(x\right)=V_{\mathrm{ND}}^n-V_{\mathrm{ND},\max }^n\end{array}\right.---\left(7\right)$

其中，为第n个节点的节点电压幅值，分别为第n个节点的节点电压幅值的下限和上限；

第l个负荷的有功功率和无功功率第b条支路电流幅值第n个节点电压幅值均为状态变量x的函数，通过下式计算得到；

$\left\{\begin{array}{c}{P}_{\mathrm{LD}}+{\mathrm{iQ}}_{\mathrm{LD}}=S_n\\ I_{\mathrm{LN}}=\left|I\right|\\ V_{\mathrm{ND}}=\left|V\right|\end{array}\right.---\left(8\right)$

其中，负荷有功功率向量其中n_n为节点总数；负荷无功功率向量 $Q_{\mathrm{LD}}={[Q_{\mathrm{LD},}^1{Q}_{\mathrm{LD}}^2,\&CenterDot;\&CenterDot;\&CenterDot;,Q_{\mathrm{LD}}^l,\&CenterDot;\&CenterDot;\&CenterDot;,Q_{\mathrm{LD}}^{n_n}]}^T;$ 支路电流幅值向量 $I_{\mathrm{LN}}={[I_{\mathrm{LN}}^1,I_{\mathrm{LN}}^2,\&CenterDot;\&CenterDot;\&CenterDot;,I_{\mathrm{LN}}^b,\&CenterDot;\&CenterDot;\&CenterDot;,I_{\mathrm{LN}}^{b_n}]}^T,$ 其中b_n为支路总数；I为支路电流向量；节点电压幅值向量V为节点电压向量；

因状态变量之间无耦合关系，协同种群集合C中种群数量n_x′等于搜索空间的维度数n_x，即种群C_j仅对一维状态变量进行寻优，j＝1,2,...,n′_x；采用与主种群M初始化相同的方法对协同种群集合C进行初始化。

所述步骤3中，协同条件为：mod(it,ci)＝＝ci，其中it为当前迭代次数，ci为协同间隔次数，mod为求取余数函数，符号＝＝表示等于。

粒子群优化，是由Kennedy和Eberhart基于一种社会心理学模型中的社会影响和社会学习提出的一种智能算法。粒子群中的每一个个体遵循简单的行为，即效仿相邻个体的成功经验而展开的累积行为。粒子群优化模型中，个体最优点代表从仿真开始被这个个体经历过的最好位置，邻域最优点是被这个个体的所有邻居经历过的最好位置，这两个最优点被作为吸引子；个体具有个体最优点和邻域最优点的记忆，它根据一些简单的规则按一定的比例利用最优点与当前位置的距离来调整粒子的位置，使得群体在一定的迭代次数内聚集到目标附近。粒子群中个体被称为“粒子”是因为既需要将个体描述为没有质量、没有体积的，同时也需要描述它的速度和加速状态。

一个粒子群优化算法维持着一个一定数量粒子的种群，其中每个粒子都代表了问题的一个潜在解。粒子们在多维空间中飞行，粒子的位置变化由加入的速度引起，速度由自身和周围邻居的信息决定。

所述步骤4中，对主种群M进行求解的过程包括更新粒子的个体最优位置y_i(t+1)和全局最优位置计算粒子t+1时刻的速度v_i(t+1)以及粒子t+1时刻的位置x_i(t+1)；

1)粒子的个体最优位置是第i个粒子从开始到现在达到过的最佳位置，t时刻粒子的个体最优位置y_i(t)表示为：

$y_i\left(t\right)={[y_{i1}\left(t\right),y_{i2}\left(t\right),\&CenterDot;\&CenterDot;\&CenterDot;,y_{\mathrm{ij}}\left(t\right),\&CenterDot;\&CenterDot;\&CenterDot;,{y_{\mathrm{in}}}_x]}^T---\left(9\right)$

其中，y_ij(t)为第i个粒子t时刻在第j维度上的个体最优位置；

于是t+1时刻粒子的个体最优位置y_i(t+1)由下式计算：

$y_i(t+1)=\left\{\begin{array}{ccc}{y}_i\left(t\right),& \mathrm{if}& F\left(x_i(t+1)\right)\&GreaterEqual;F\left(y_i\left(t\right)\right)\\ x_i(t+1),& \mathrm{if}& F\left(x_i(t+1)\right)<F\left(y_i\left(t\right)\right)\end{array}\right.---\left(10\right)$

其中，x_i(t+1)为第i个粒子t+1时刻在搜索空间的位置，F(x_i(t+1))为第i个粒子t时刻在搜索空间的位置对应的适应度函数，F(y_i(t))为t时刻粒子的个体最优位置对应的适应度函数，其取值为：

$F\left(y_i\left(t\right)\right)=\min \{F\left(y_1\left(t\right)\right),F\left(y_2\left(t\right)\right),...,F\left({y_n}_s\left(t\right)\right)\}---\left(11\right)$

其中，n_s为单个种群中粒子的数量；

2)t+1时刻的全局最优位置表示为：

$\hat{y}(t+1)=y_i\left(t\right)---\left(12\right)$

在t时刻，对于任意第i个粒子t时刻的全局最优位置有：

${\hat{y}}_i\left(t\right)=[{\hat{y}}_{i1}\left(t\right),{\hat{y}}_{i2}\left(t\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\hat{y}}_{\mathrm{ij}}\left(t\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\widehat{y_{\mathrm{in}}}}_x{]}^T---\left(13\right)$

其中，为第i个粒子在第j维度上的全局最优位置；

3)计算粒子t+1时刻的速度v_i(t+1)；

在t时刻，粒子t时刻的速度其中v_ij(t)为第i个粒子第j维度上的速度，于是t+1时刻第i个粒子在第j维度上的速度v_ij(t+1)表示为：

$v_{\mathrm{ij}}(t+1)={\mathrm{wv}}_{\mathrm{ij}}\left(t\right)+c_1{r}_{1j}\left(t\right)[y_{\mathrm{ij}}\left(t\right)-x_{\mathrm{ij}}\left(t\right)]+c_2{r}_{2j}\left(t\right)[{\hat{y}}_{\mathrm{ij}}\left(t\right)-x_{\mathrm{ij}}\left(t\right)]---\left(14\right)$

其中，w为惯性因子，衡量t时刻的速度对t+1时刻移动的影响；c₁和c₂分别为认知因子和社会因子，是正数的加速度常量，分别用来度量认知成分和社会成分对于速度更新的贡献；r_1j(t)、r_2j(t)～U(0，1)都是在区间[0，1]中均匀抽取的随机数，代表不确定性因素；x_ij(t)为第i个粒子t时刻在第j维度上在搜索空间的位置，j＝1,2,…,n_x，n_x为搜索空间的维度数；

4)粒子t+1时刻的位置x_i(t+1)表示为：

x_i(t+1)＝x_i(t)+v_i(t+1)(15)

其中，x_i(t)第i个粒子t时刻在搜索空间的位置，i＝1,2,…,n_s，其中n_s为单个种群中粒子的数量；x_i(t)表示为：

$x_i\left(t\right)={[x_{i1}\left(t\right),x_{i2}\left(t\right),\&CenterDot;\&CenterDot;\&CenterDot;,x_{\mathrm{ij}}\left(t\right),\&CenterDot;\&CenterDot;\&CenterDot;,x_{{\mathrm{in}}_x}]}^T---\left(16\right)$

其中，t代表的是离散时间点，也代表迭代次数；初始位置x_i(0)～U(x_min,x_max)，x_min和x_max分别为粒子位置的最小值向量和最大值向量。

所述步骤5中，收敛条件为迭代次数超过最大次数限制，或目标函数值小于给定的阈值，即：

$\mathrm{it}>\mathrm{it}\_\lim \left|\right|F(M.\hat{y})<\mathrm{\&epsiv;}---\left(17\right)$

其中：it为当前迭代次数，it_lim为迭代次数限制，ε为收敛标准，符号||表示逻辑关系或，为主种群M中全局最优位置对应的适应度值。

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

步骤1：更新协同种群集合C；

对种群C_j，某个粒子i，式中i～U(1，n_s/2)且i∈N，其中，C_j.x_i为种群C_j中第i粒子在搜索空间的位置，为主种群M全局最优位置中第j维度上的数值，n_s为单个种群中粒子的数量，N为自然数；

步骤2：定义向量 为主种群M中全局最优位置，b_j为向量b中第j维度上的数值，种群C_j状态变量定义为：其中n_x′为协同种群集合C中种群数量；并按照步骤4的方法对每个种群C_j进行求解，计算完成后为种群C_j的全局最优位置；

步骤3：更新主种群M：M.x_i＝b，其中M.x_i为主种群M中第i个粒子在搜索空间的位置，i～U(1,n_s/2)且k∈N。

状态估计中使用的量测可分为三大类：实时量测、伪量测、虚拟量测。实时量测、虚拟量测精度较高，伪量测精度较低。但实时量测装置成本较高，PMU仅安装在变电站内，RTU设置在出线开关、分段开关和联络开关上，我国发达地区一条馈线上一般仅配置4～5个实时量测点；受通信条件的限制，仅变电站、开关站、环网柜中的实时数据能够上传。配电网负荷数量多、密度大，浮游节点的比例约为10％。配电变压器、负荷数据来自营销自动化系统，上传周期一般为5～30分钟，不能作为实时数据。综上所述，实时量测、虚拟量测数量太少不足以进行状态估计，因此必须大量使用精度较低的伪量测数据。

配电网状态估计以f(x)为最小为目标函数，具体有：

$\min f\left(x\right)=\underset{m=1}{\overset{n_m}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_m{[z_m-h_m\left(x\right)]}^2---\left(18\right)$

其中，n_m为量测量的数量，ω_m第m个量测量的测量权重，z_m为第m个量测量，h_m为第m个量测量的状态函数；

1)所述目标函数的不等式约束包括负荷有功功率约束、负荷无功功率约束、支路电流约束和节点电压约束，分别表示为：

$\left\{\begin{array}{c}{P}_{\mathrm{LD},\min }^l\&le;P_{\mathrm{LD}}^l\&le;P_{\mathrm{LD},\max }^l\\ Q_{\mathrm{LD},\min }^l\&le;Q_{\mathrm{LD}}^l\&le;Q_{\mathrm{LD},\max }^l\\ I_{\mathrm{LN}}^b\&le;I_{\mathrm{LN},\max }^b\\ V_{\mathrm{ND},\min }^n\&le;V_{\mathrm{ND}}^n\&le;V_{\mathrm{ND},\max }^n\end{array}\right.---\left(19\right)$

2)所述目标函数的等式约束为潮流约束，基于状态变量运用潮流相关算法求取测量量对应的电气量，其值满足潮流约束；选择节点电压为状态变量，状态方程为：

$\left\{\begin{array}{c}V=V\\ I=Y_b{A}^{T}V\\ S_b=V\left(N_f\right).*\stackrel{\&OverBar;}{I}\\ S_n=V.*\stackrel{\&OverBar;}{(Y*V)}\end{array}\right.---\left(20\right)$

其中，V为n_n×1维节点电压向量，n_n为节点总数，搜索空间的维度数n_x＝2*n_n；I为b_n×1维支路电流向量，b_n为支路总数；Y_b为支路导纳矩阵；A为节点-支路关联矩阵；S_b为b_n×1维支路首端复功率向量，N_f为支路首端节点编号向量，V(N_f)为b_n×1维支路首端电压向量，符号.*表示矩阵对应元素相乘；S_n为n_n×1维节点注入复功率，Y为节点导纳矩阵。

为了验证提出的协同PSO算法在配电网状态估计中的有效性，采用如下配置条件开展实验。量测配置如表1所示，V_b为平衡节点的复电压，P_LD、Q_LD为负荷有功、无功，P_b、Q_b为支路首端有功、无功，特征形式为N(μ,σ²)，μ为期望，σ为标准差；表中数据的物理含义为测量值偏离真值的百分比。参数配置如表2所示，C_j表示第j个子种群，σ²为正态分布的方差，n_x为搜索空间的维度数，c_k为常数值，ci为协同间隔次数。

表1

1.精确性分析

以IEEE33-bus配电系统为实验对象，量测配置设置为N(0,0.03²)，基于节点电压的加权最小二乘状态估计算法(WeightedLeastSquares,WLS)与CPSO计算结果的相对误差分布情况如附图2所示，详细的概率统计数值见表3。CPSO与WLS标准差均小于0.03，均发挥了滤波作用；二者在均值、标准差、6σ、3σ、2σ、1σ上的差值都很小；表明CPSO与WLS计算精度非常接近。

表3

2.适用性分析

以IEEE33节点配电系统为基础实验对象，所有负荷增大给定倍数，量测配置设置为N(0,0.03²)。负荷倍率与迭代次数之间的关系如附图3所示；计算时间与结果精度也发生相应变化，详细数据见表4，表中NaN表示无法完成数值计算。

表4

由表可知，当负荷倍率小于等于3.6时，CPSO与WLS都能正常收敛。在正常收敛阶段，CPSO与WLS均能计算出正确的结果；随着倍率增大，CPSO迭代次数在一定范围内较稳定地波动，WLS迭代次数则不断增加，使得计算时间在倍率接近3时超过CPSO。在不正常收敛阶段，CPSO表现为达到最大迭代次数仍然搜索不到满足需求的结果；WLS则以较少次数和计算时间停止迭代，但计算结果不是正确的解，或者信息矩阵成为奇异矩阵，无法完成数值计算。由此可见，CPSO、WLS均需满足一定的限制条件，CPSO适用于负荷较重、控制装置较多的网络，WLS适用于负荷较轻、控制装置较少的网络。

3.收敛性分析

以IEEE33、PG&E69、实际292节点、588节点、1180节点配电网络为实验对象，选择基本粒子群算法(PSO)、多启动粒子群算法(MPSO)、基于遗传思想的混合粒子群算法(GA-PSO)、基于模拟退火的混合粒子群算法(SA-PSO)、基于混沌思想的混合粒子群算法(Chaos-PSO)作为比较对象，量测配置设置为N(0,0.03²)。

各算法性能指标如表5所示：

表5

实验次数为10，为平均计算时间，c为平均迭代次数，r为目标函数满足收敛标准的实验百分数。

$\left\{\begin{array}{c}\stackrel{\&OverBar;}{t}=\frac{1}{T}\underset{h=1}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{Time}}_h\\ c=\frac{1}{T}\underset{\stackrel{\&OverBar;}{t}=h}{\overset{T}{\mathrm{\&Sigma;}}}{\mathrm{it}}_h\\ r=\frac{T_{\mathrm{qia}}}{T}*100\%\end{array}\right.---\left(21\right)$

式中：Time_h为第h次试验计算消耗的时间(秒)，it_h为第h次试验的迭代次数，T_qua为适应度函数值小于收敛标准ε的实验次数。

在所有实验对象特别是大规模配电系统中，CPSO均以较少的迭代次数达到了目标函数收敛标准，计算时间也相对较短，从而说明提出算法的收敛能力非常突出。

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

Claims (8)

1.一种基于协同粒子群的配电网状态估计方法，其特征在于：所述方法包括以下步骤：

步骤1：载入配网系统及其量测配置；

步骤2：初始化主种群M和协同种群集合C；

步骤3：判断是否满足协同条件，若满足执行步骤4，否则执行步骤6；

步骤4：对主种群M中粒子进行求解；

步骤5：判断是否满足收敛条件，若满足则状态估计结束，否则转步骤3；

步骤6：对协同种群集合C进行求解。

2.根据权利要求1所述的基于协同粒子群的配电网状态估计方法，其特征在于：所述步骤1中，进行拓扑分析，确定搜索方向，进行节点编号，生成节点导纳矩阵Y与节点-支路关联矩阵A。

3.根据权利要求1所述的基于协同粒子群的配电网状态估计方法，其特征在于：所述步骤2中，基于量测的节点注入复功率与平衡节点电压进行潮流分析算出节点电压向量V，以之作为初始化参考值；初始化主种群M，包括粒子当前位置、个体最优位置和全局最优位置，并计算其对应适应度，初始速度V(0)＝0；

适应度函数F(x)表示为：

$F\left(x\right)=f\left(x\right)+k_f\underset{k=1}{\overset{N_{ue}}{\&Sigma;}}(max\&lsqb;0,g_k^j(x)\&rsqb;)---\left(1\right)$

其中，f(x)为配电网状态估计函数；k_f为惩罚因子；N_ue为不等式约束数；为惩罚函数，分为以下几种情况：

1)负荷有功功率对应的惩罚函数和表达式为：

$\left\{\begin{array}{c}{g}_k^1\left(x\right)=P_{LD,min}^l-P_{LD}^l\\ g_k^2\left(x\right)=P_{LD}^l-P_{LD,\max }^l\end{array}\right.---\left(2\right)$

其中，为第l个负荷的有功功率，和分别为第l个负荷的有功功率的下限和上限；

2)负荷无功功率对应的惩罚函数和表达式为：

$\left\{\begin{array}{c}{g}_k^3\left(x\right)=Q_{LD,\min }^l-Q_{LD}^l\\ g_k^4\left(x\right)=Q_{LD}^l-Q_{LD,\max }^l\end{array}\right.---\left(3\right)$

其中，为第l个负荷的无功功率，分别为第l个负荷的无功功率的下限和上限；

3)支路电流幅值对应的惩罚函数表达式为：

$g_k^5\left(x\right)=I_{LN}^b-I_{LN,max}^b---\left(6\right)$

其中，为第b条支路的支路电流幅值，为第b条支路的支路电流幅值的最大值；

4)节点电压幅值对应的惩罚函数和表达式为：

$\left\{\begin{array}{c}{g}_k^6\left(x\right)=V_{ND,min}^n-V_{ND}^n\\ g_k^7\left(x\right)=V_{ND}^n-V_{ND,\max }^n\end{array}\right.---\left(7\right)$

其中，为第n个节点的节点电压幅值，分别为第n个节点的节点电压幅值的下限和上限；

第l个负荷的有功功率和无功功率第b条支路电流幅值第n个节点电压幅值均为状态变量x的函数，通过下式计算得到；

$\left\{\begin{array}{c}{P}_{LD}+{\mathrm{iQ}}_{LD}=S_n\\ I_{LN}=\left|I\right|\\ V_{ND}=\left|V\right|\end{array}\right.---\left(8\right)$

其中，负荷有功功率向量其中n_n为节点总数；负荷无功功率向量 $Q_{LD}={\&lsqb;Q_{LD}^1,Q_{LD}^2,...,Q_{LD}^l,...,Q_{LD}^{n_n}\&rsqb;}^T;$ 支路电流幅值向量 $I_{LN}={\&lsqb;I_{LN}^1,I_{LN}^2,...,I_{LN}^b,...,I_{LN}^{b_n}\&rsqb;}^T,$ 其中b_n为支路总数；I为支路电流向量；节点电压幅值向量V为节点电压向量；

因状态变量之间无耦合关系，协同种群集合C中种群数量n′_x等于搜索空间的维度数n_x，即种群C_j仅对一维状态变量进行寻优，j＝1,2,...,n′_x；采用与主种群M初始化相同的方法对协同种群集合C进行初始化。

4.根据权利要求1所述的基于协同粒子群的配电网状态估计方法，其特征在于：所述步骤3中，协同条件为：mod(it,ci)＝＝0，其中it为当前迭代次数，ci为协同间隔次数，mod为求取余数函数，符号＝＝表示等于。

5.根据权利要求1所述的基于协同粒子群的配电网状态估计方法，其特征在于：所述步骤4中，对主种群M进行求解的过程包括更新粒子的个体最优位置y_i(t+1)和全局最优位置计算粒子t+1时刻的速度v_i(t+1)以及粒子t+1时刻的位置x_i(t+1)；

1)粒子的个体最优位置是第i个粒子从开始到现在达到过的最佳位置，t时刻粒子的个体最优位置y_i(t)表示为：

$y_i\left(t\right)={\&lsqb;y_{i1}\left(t\right),y_{i2}\left(t\right),...,y_{ij}\left(t\right),...,y_{{\mathrm{in}}_x}\&rsqb;}^T---\left(9\right)$

其中，y_ij(t)为第i个粒子t时刻在第j维度上的个体最优位置；

于是t+1时刻粒子的个体最优位置y_i(t+1)由下式计算：

$y_i(t+1)=\left\{\begin{array}{ccc}{y}_i\left(t\right),& if& F\left(x_i\right(t+1\left)\right)\&GreaterEqual;F\left(y_i\right(t\left)\right)\\ x_i(t+1),& if& F\left(x_i\right(t+1\left)\right)<F\left(y_i\right(t\left)\right)\end{array}\right.---\left(10\right)$

其中，x_i(t+1)为第i个粒子t+1时刻在搜索空间的位置，F(x_i(t+1))为第i个粒子t时刻在搜索空间的位置对应的适应度函数，F(y_i(t))为t时刻粒子的个体最优位置对应的适应度函数，其取值为：

$F\left(y_i\right(t\left)\right)=min\{F\left(y_1\right(t\left)\right),F\left(y_2\right(t\left)\right),...,F\left(y_{n_s}\right(t\left)\right)\}---\left(11\right)$

其中，n_s为单个种群中粒子的数量；

2)t+1时刻的全局最优位置表示为：

$\hat{y}(t+1)=y_i\left(t\right)---\left(12\right)$

在t时刻，对于任意第i个粒子t时刻的全局最优位置有：

${\hat{y}}_i\left(t\right)={\&lsqb;{\hat{y}}_{i1}\left(t\right),{\hat{y}}_{i2}\left(t\right),...,{\hat{y}}_{ij}\left(t\right),...,{\hat{y}}_{{\mathrm{in}}_x}\&rsqb;}^T---\left(13\right)$

其中， 为第i个粒子在第j维度上的全局最优位置；

3)计算粒子t+1时刻的速度v_i(t+1)；

在t时刻，粒子t时刻的速度其中v_ij(t)为第i个粒子第j维度上的速度，于是t+1时刻第i个粒子在第j维度上的速度v_ij(t+1)表示为：

$v_{ij}(t+1)={\mathrm{wv}}_{ij}\left(t\right)+c_1{r}_{1j}\left(t\right)\&lsqb;y_{ij}\left(t\right)-x_{ij}\left(t\right)\&rsqb;+c_2{r}_{2j}\left(t\right)\&lsqb;{\hat{y}}_{ij}\left(t\right)-x_{ij}\left(t\right)\&rsqb;---\left(14\right)$

其中，w为惯性因子，衡量t时刻的速度对t+1时刻移动的影响；c₁和c₂分别为认知因子和社会因子，是正数的加速度常量，分别用来度量认知成分和社会成分对于速度更新的贡献；r_1j(t)、r_2j(t)～U(0，1)都是在区间[0，1]中均匀抽取的随机数，代表不确定性因素；x_ij(t)为第i个粒子t时刻在第j维度上在搜索空间的位置，j＝1,2,…,n_x，n_x为搜索空间的维度数；

4)粒子t+1时刻的位置x_i(t+1)表示为：

x_i(t+1)＝x_i(t)+v_i(t+1)(15)

其中，x_i(t)第i个粒子t时刻在搜索空间的位置，i＝1,2,…,n_s，其中n_s为单个种群中粒子的数量；x_i(t)表示为：

$x_i\left(t\right)={\&lsqb;x_{i1}\left(t\right),x_{i2}\left(t\right),...,x_{ij}\left(t\right),...,x_{{\mathrm{in}}_x}\&rsqb;}^T---\left(16\right)$

其中，t代表的是离散时间点，也代表迭代次数；初始位置x_i(0)～U(x_min,x_max)，x_min和x_max分别为粒子位置的最小值向量和最大值向量。

6.根据权利要求1所述的基于协同粒子群的配电网状态估计方法，其特征在于：所述步骤5中，收敛条件为迭代次数超过最大次数限制，或目标函数值小于给定的阈值，即：

$it>it\_\lim \left|\right|F(M.\hat{y})<\mathrm{\&epsiv;}---\left(17\right)$

其中：it为当前迭代次数，it_lim为迭代次数限制，ε为收敛标准，符号||表示逻辑关系或，为主种群M中全局最优位置对应的适应度值。

7.根据权利要求1所述的基于协同粒子群的配电网状态估计方法，其特征在于：所述步骤6包括以下步骤：

步骤(1)：更新协同种群集合C；

对种群C_j，某个粒子i，式中i～U(1，n_s/2)且i∈N，其中，C_j.x_i为种群C_j中第i粒子在搜索空间的位置，为主种群M全局最优位置中第j维度上的数值，n_s为单个种群中粒子的数量，N为自然数；

步骤(2)：定义向量 为主种群M中全局最优位置，b_j为向量b中第j维度上的数值，种群C_j状态变量定义为：其中n′_x为协同种群集合C中种群数量；并按照步骤4的方法对每个种群C_j进行求解，计算完成后 为种群C_j的全局最优位置；

步骤(3)：更新主种群M：M.x_i＝b，其中M.x_i为主种群M中第i个粒子在搜索空间的位置，i～U(1,n_s/2)且k∈N。

8.根据权利要求3所述的基于协同粒子群的配电网状态估计方法，其特征在于：配电网状态估计以f(x)为最小为目标函数，具体有：

$\min f\left(x\right)=\underset{m=1}{\overset{n_m}{\&Sigma;}}{\mathrm{\&omega;}}_m{\&lsqb;z_m-h_m\left(x\right)\&rsqb;}^2---\left(18\right)$

其中，n_m为量测量的数量，ω_m第m个量测量的测量权重，z_m为第m个量测量，h_m为第m个量测量的状态函数；

1)所述目标函数的不等式约束包括负荷有功功率约束、负荷无功功率约束、支路电流约束和节点电压约束，分别表示为：

$\left\{\begin{array}{c}{P}_{LD,min}^l\&le;P_{LD}^l\&le;P_{LD,max}^l\\ Q_{LD,min}^l\&le;Q_{LD}^l\&le;Q_{LD,max}^l\\ I_{LN}^b\&le;I_{LN,\max }^b\\ V_{ND,\min }^n\&le;V_{ND}^n\&le;V_{ND,max}^n\end{array}\right.---\left(19\right)$

2)所述目标函数的等式约束为潮流约束，基于状态变量运用潮流相关算法求取测量量对应的电气量，其值满足潮流约束；选择节点电压为状态变量，状态方程为：

$\left\{\begin{array}{c}V=V\\ I=Y_b{A}^{T}V\\ S_b=V\left(N_f\right).*\stackrel{\&OverBar;}{I}\\ S_n=V.*\stackrel{\&OverBar;}{(Y*V)}\end{array}\right.---\left(20\right)$

其中，V为n_n×1维节点电压向量，n_n为节点总数，搜索空间的维度数n_x＝2*n_n；I为b_n×1维支路电流向量，b_n为支路总数；Y_b为支路导纳矩阵；A为节点-支路关联矩阵；S_b为b_n×1维支路首端复功率向量，N_f为支路首端节点编号向量，V(N_f)为b_n×1维支路首端电压向量，符号.*表示矩阵对应元素相乘；S_n为n_n×1维节点注入复功率，Y为节点导纳矩阵。