配电网出现故障后保护通常要计算故障后的稳态电流或电压值，对于一些应用情况，如配变支路保护、用户侧保护等，希望能尽快计算故障后的稳态电流或电压值，以保护快速出口断开故障，为主线路保护配合、备用电源投入等创造良好条件。由于配电网故障后的电流电压信号含有暂态分量，为了能在短时间内(如10 ms)计算故障后的稳态电流和电压值，就必然需要在短时间内处理电流电压暂态信号。

目前，通常采用加窗傅里叶等滤波算法^[1-3]来提取稳态分量，暂态分量的存在影响了滤波算法的提取精度，增加窗的宽度可提高滤波算法的提取精度，但增加了信号处理时间。为了既不增加信号处理时间又能提高信号处理精度，一些方法将稳态分量和暂态分量组合建模，采用拟合逼近方法求取故障后的稳态分量和暂态分量，多项式拟合方法和扩展Prony算法是具有代表性的两种方法。多项式拟合方法利用一定阶数的多项式来拟合暂态分量，多项式的阶数影响暂态分量的拟合效果，增加多项式阶数原则上可促进拟合效果，但增加了待求参数个数^[4-6]。扩展Prony算法利用带指数衰减的正弦函数组合来拟合暂态分量，其待求参数较多。随着待求参数的增多，求取参数的困难性增加，如最优解问题、参数约束问题等，且增加了采样数据长度^[7-8]。

为了能在短时间内有效计算故障后的稳态电流和电压值，文中提出一种拟合逼近处理方法，可以在短时间内求解配电网故障暂态信号。方法中，电流电压暂态信号以故障后稳态分量、直流衰减分量和其他高频衰减分量的组合进行表示，且视高频衰减分量成分很小，运用最小二乘法求解稳态分量和直流衰减分量的参数，以计算故障后的稳态电流或电压值。文中就提出的配电网故障暂态信号拟合逼近处理方法进行了仿真计算分析，分析结果表明，所提方法能在短时间内有效计算故障后的稳态电流或电压值。

1 拟合逼近处理方法原理

配电网故障后电流电压暂态信号可表示为稳态分量、直流衰减分量及其他高频衰减分量的组合，即

$ \begin{array}{l} i\left( t \right) = {I_m}{\rm{sin}}(\omega t + {\varphi _i}) + {C_i}{{\rm{e}}^{-\frac{t}{{{\tau _i}}}}} + {G_i}\left( t \right), \\ u\left( t \right) = {U_m}{\rm{sin}}(\omega t + {\varphi _u}) + {C_u}{{\rm{e}}^{-\frac{t}{{{\tau _u}}}}} + {G_u}\left( t \right), \end{array} $

(1)

式中：I_m、φ_i、C_i、τ_i、G_i(t)和U_m、φ_u、C_u、τ_u、G_u(t)分别为电流电压故障后的稳态分量幅值、初相位、直流衰减分量的幅值、衰减时间常数、其他高频衰减分量。式(1)可表示为

$ \begin{array}{l} i\left( t \right) = {A_i}{\rm{cos}}\left( {\omega t} \right) + {B_i}{\rm{sin}}\left( {\omega t} \right) + {C_i}{{\rm{e}}^{-\frac{t}{{{\tau _i}}}}} + {G_i}\left( t \right), \\ u\left( t \right) = {A_u}{\rm{cos}}\left( {\omega t} \right) + {B_u}{\rm{sin}}\left( {\omega t} \right) + {C_u}{e^{-\frac{t}{{{\tau _u}}}}} + {G_u}\left( t \right), \end{array} $

(2)

式中，A_i=I_msinφ_i，B_i=I_mcosφ_i，A_u=U_msinφ_u，B_u=U_mcosφ_u。

考虑到配电网暂态高频衰减分量通常衰减很快，且远小于暂态直流衰减分量的特点^[9]，可认为直流衰减分量是故障暂态分量的主要成分，由此可将式(2)中的A、B、C、τ 4个参数作为待辨识参数，其参数的辨识可转换为求解一个带约束的非线性最小二乘问题^[10-11]

$ \left\{ \begin{array}{l} {\rm{min}}s\left( x \right) = {\left\| {\mathit{\boldsymbol{f}}\left( x \right)} \right\|^2} \in {R^4}, \\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;{x_4} > 0, \end{array} \right. $

(3)

$ f\left( x \right) = {({f_1}\left( x \right), {f_2}\left( x \right), \ldots \ldots {f_m}\left( x \right))^{\rm{T}}}, $

(4)

$ {f_i}\left( x \right) = {x_1}{\rm{cos}}(w{t_i}) + {x_2}{\rm{sin}}(w{t_i}) + {x_3}{{\rm{e}}^{-\frac{{{t_i}}}{{{x_4}}}}}-y({t_i})。$

(5)

式(3)为目标函数和约束条件，式(4)中m的取值取决于采样点数，式(5)中y(t_i)为式(2)中的电流电压采样值i(t_i)或u(t_i)，f_i(x)对应时间t_i，不同的时刻分别对应一个约束方程，f_i(x)是关于x的非线性函数，且具有一阶连续偏导数。

2 算法实现 2.1 初值选取

针对非线性最小二乘问题，传统的解决方法是利用高斯牛顿法。但高斯牛顿法存在对初值比较敏感问题，初值选取不好往往会造成收敛速度慢甚至不收敛。随着智能优化技术的发展，粒子群优化算法(PSO, particle searm optimization)具有原理简单、易于实现、无需梯度信息、参数少等优点，在连续优化和离散优化问题中都表现出良好的效果，但缺乏局部区域精细搜索能力^[10-12]。因此，文中算法选用粒子群算法求取高斯牛顿算法的初始值。

2.2 带约束粒子群算法

粒子群算法具有较好的全局收敛特点，而在此算法中，根据实际物理意义，直流衰减时间常数τ＞0，需要在粒子群算法中引入约束。

PSO算法的约束主要有罚函数法和搜索限制法，罚函数的目的是将约束优化问题转化成无约束优化问题，但罚函数的设计通常比较复杂。搜索限制法将粒子群的搜索范围限制在可行域内，进而保证解的可行性，但初始可行解空间往往很难确定^[13-14]。算法结合模拟退火算法的思想^[15]，以一定的概率接受可行域外的粒子。概率公式如下：

$ P = {{\rm{e}}^{-\frac{D}{{{T_k}}}}}, $

(6)

$ D =-{x_4}, $

(7)

式中：D表示粒子超出可行域的距离，D越大，被接受的概率就越小；T_k与模拟退火算法中的温度类似，采用比例下降方式，随着迭代次数进行，T_k减小，被接受概率越小。当粒子在搜索过程中第4个分量x₄＜0时，计算概率^P。如果

$ P > \xi, (\xi = {\rm{rand}}\left( {0, 1} \right)), $

(8)

则接受此粒子; 否则, 保留上一次粒子位置。这样既保证初始阶段粒子的搜索能力，又保证最终解在可行域内。

2.3 算法流程

步骤1:初始化粒子随机位置和随机速度，设定初始T_k和终止条件。

步骤2:计算各粒子适应值s(x_j)=‖f(x_j)‖²，j表示第j个粒子。并与自身和全局粒子历史最好位置比较，如果更好，更新其位置。

步骤3:利用${v_j}^{k + 1} = {v_j}^k + {c_1}\xi ({p_j}^k-{x_j}^k) + {c_2}\eta ({p_g}^k-{x_j}^k)$，${x_j}^{k + 1} = {x_j}^k + {v_j}^{k + 1}$对粒子群位置进行更新。

步骤4:若x_j^k+1在可行域，转步骤5；否则, 产生ξ=rand(0, 1)，计算D，若$P = {{\rm{e}}^{-\frac{D}{{{T_k}}}}} \ge \xi $，转步骤5；否则令x_j^k+1=x_j^k。

步骤5:利用T_k+1=T_k·r, r∈(0.95, 0.99)降低T_k，若达到终止条件，转步骤6；否则转步骤2。

步骤6:将粒子群算法所得最优解作高斯牛顿法初始值x₀，给定允许误差ε₁＞0，ε₂＞0, 令迭代次数k=0。

步骤7：检查是否满足收敛准则。计算f(x_k)及Jacobi矩阵▽f(x_k)，若‖▽f(x_k)^Tf(x_k)‖＜ε₁或‖f(x_k)‖²＜ε₂，迭代终止，x_k为最优解；否则步骤8。

步骤8：构造高斯-牛顿方向。取d_k=[▽f(x_k)^T▽f(x_k)]^-1▽f(x_k)^Tf(x_k)。

步骤9：求x_k+1，设λ_k初值和m₀，x_k+1=x_k+λ_kd_k。若$\frac{{\left\| {\mathit{\boldsymbol{f}}({x_k})} \right\|}}{{\left\| {f({x_{k + 1}})} \right\|}} > \delta $，存储x_k+1，转步骤10；否则λ_k=λ_k/2，m=m+1，若m＞m₀，转步骤10；否则转步骤9。

步骤10：约束条件检查。若满足k=k+1，转步骤7，否则，调整x₄初值，转步骤7。

3 仿真分析

为了验证文中所给方法的有效性，在数据窗取10 ms的情况下，就暂态高频衰减分量大小以电流信号对方法的影响为例进行仿真分析。

设电流输入信号：${i_1}\left(t \right) = 100{\rm{cos}}({\omega _0}t + {\rm{ \mathsf{ π} }}/6) + 50{{\rm{e}}^{-t/0.05}} + {G_m}{{\rm{e}}^{-t/{T_2}}}{\rm{cos}}(n{\omega _0}t + {\rm{ \mathsf{ π} }}/4)$。

式中取ω₀=2πf₀，f₀=50 Hz。分析计算结果如表 1~表 3所示。由于配网中高频分量的主要成分在500 Hz以内^[9]，因此，n在10以内取值即可。表 1是高频衰减分量n取2的计算结果，表 2是高频衰减分量n取5的计算结果，表 3是高频衰减分量n取10的计算结果。

由分析计算结果可看出：若暂态高频衰减分量为零，在数据窗取10 ms的情况下，文中方法可精确计算工频稳态分量的幅值和初相角；随着暂态高频衰减分量幅值的增大，工频稳态分量幅值和初相角的计算精度逐渐变差，甚至变得无法收敛；工频稳态分量幅值和初相角的计算精度与暂态高频衰减分量的衰减快慢相关，衰减愈快，计算精度愈高；暂态高频衰减分量的频率高低也对工频稳态分量幅值和初相角的计算精度产生影响，频率愈高，计算精度愈差。

由于通常情况下，配电网暂态高频衰减分量很小且衰减很快，因此，文中方法是一种针对实际配网运行条件的实用方法。文中方法已在配网用户侧保护得到应用^[16]。其应用通过求取电网侧电压工频分量的变化，快速判断电网侧的故障。

4 结论

考虑到配电网暂态高频衰减分量小、衰减快的特点，文中方法采用配电网故障后的稳态工频分量和直流衰减分量组合建模，减少了待求参数个数、缩短采样数据长度。结合改进的带约束的粒子群算法寻找初值，提高了算法的整体收敛性，防止算法陷入局部最优。通过模拟退火引入约束，使所求参数符合工程实际，具有较好的鲁棒性。