基于改进约束ekf算法的动态振荡信号参数辨识方法
【技术领域】
[0001] 本发明涉及一种基于改进约束EKF算法的动态振荡信号参数辨识方法，属于信号 分析与参数辨识技术领域。
【背景技术】
[0002] 近年来，随着现代电网规模不断扩大，电网互联程度的日益提高，系统受到大、小 扰动后产生的动态振荡已经成为制约电网安全稳定运行的最主要的因素之一。由于这些动 态振荡信号可以提供关于电力系统运行模式的重要信息，所以发现并准确掌握这些振荡信 号特征对于电力系统安全稳定运行具有重要意义。
[0003] 鉴于动态振荡信号辨识的重要性，研究人员提出许多辨识方法，如矩阵束法、最大 似然法、普龙尼法等。然而，这些方法大都具有窗口属性，无法实现动态振荡信号的实时在 线辨识。最近，研究人员提出了一种基于EKF算法的参数辨识方法，克服了之前许多算法的 缺点，有效的实现了振荡信号的实时在线辨识。但是，该方法并未考虑待辨识参数所受到的 实际约束条件，所以该方法的收敛性较差，且未考虑工程实用背景。因此，有必要对实际约 束条件下的动态振荡信号参数辨识问题进行研究。
【发明内容】
[0004] 发明目的：针对现有技术中存在的问题，为了有效的解决动态振荡信号参数实际 约束条件下的实时辨识，克服基于传统EKF算法辨识的缺点，本发明提供一种基于改进约束 EKF算法的动态振荡信号参数辨识方法，有效的实现了动态振荡信号参数约束条件下的辨 识。
[0005] 技术方案:一种基于改进约束EKF算法的动态振荡信号参数辨识方法，该方法在计 算机中是依次按照如下步骤实现的：
[0006] (1)、获取状态分量中包含动态振荡信号待辨识参数的状态空间表达式；
[0007] (2)、初始化，包括:设定参数辨识的初值％，初始参数辨识误差协方差&，以及过程 噪声和量测噪声所满足的协方差矩阵Q和R，整体算法迭代次数最大值S和粒子群寻优最大 迭代次数Μ;
[0008] (3)、由已知的k-1时刻的状态估计值和状态估计误差协方差，利用传统扩展卡尔 曼滤波的预测步，得到k时刻的状态预测值和状态预测误差协方差，计算公式为：
[0009]
[0010]
[0011]式中/?表示k时刻的状态预测值，f( ·)对应具体问题状态方程中的非线性函数， 表示k-1时刻的状态估计向量，Uk-ι表示k-1时刻的控制输入。為表示k时刻的状态预测误 差协方差，
表示非线性函数/〇在^^处的雅克比矩阵，&4表示k_l时刻的 状态估计误差协方差，上标T表示转置，Qh是系统噪声k-Ι时刻所满足的协方差矩阵。
[0012] (4)、在上一步基础上，利用扩展卡尔曼滤波的滤波步，得到k时刻的状态估计值， 计算步骤为：
[0013]
[0014]
[0015]
[0016] 式中，Kk表示k时刻的卡尔曼滤波增益，為表示k时刻的状态预测误差协方差，上标 T表示转置，
表示非线性函数h( ·)在為处的雅克比矩阵，其中h( ·)对应具体 问题输出方程中的非线性函数。Rk是量测噪声k时刻所满足的协方差矩阵，I是和状态向量 维度相同的单位矩阵，毛表示k时刻的状态估计向量，yk是k时刻输出方程的输出量。
[0017] (5)、判定k时刻的参数辨识结果是否满足相应的实际约束条件。若满足，则直接运 用EKF再次迭代辨识。
[0018] (6)、若不满足，则运用改进的约束算法对该时刻参数辨识结果进行约束调整，得 到约束优化目标函数为：
[0019]
[0020] 式中，意表示k时刻待求寻优状态估计值，yk是k时刻输出量，WjPW2是正定的权矩 阵，为方便起见，本发明实施例中选取为单位矩阵，式中xk满足约束条件如下式：
[0021] Dx<d
[0022] D是已知的sXn常数行满秩矩阵，s是受约束条件参数的个数，η是状态向量的维 数，显然，s < η，d是已知的常数约束条件。
[0023] (7)、在步骤(6)约束优化目标函数中，对k时刻的估计输出量Λ进行求解时，为了 充分利用已知的量测数据信息，本发明采用卡尔曼平滑方法，计算公式如下：
[0024]
[0025]
[0026]
[0027] 式中，表示k-Ι时刻的平滑状态估计向量，L是卡尔曼平滑增益，九是k时刻的 估计输出量。
[0028] (8)、在上一步的基础上，借助罚函数方法，通过对原目标函数加上一个惩罚项，把 约束的优化问题转化为一个无约束的优化问题，其转换的计算公式如下：
[0029]
[0030] 式中，f(x)是原目标函数，h(k)是动态更新的惩罚值，一般情况下或者
k是当前的迭代次数。Η(χ)是惩罚因子，计算公式如下：
[0031]
[0032]式中，0(qi(x))是多级分配函数，qi(x)是与违反约束条件有关的函数，qi(x)=max {0，gi(x)}，i = l，···，s，其中gi(x)=Dx_d。γ (qi(x))表示罚函数的效力。函数取值遵循的规 则为：
[0033]
[0034]
[0035] (9)、在上一步的基础上，利用改进的粒子群算法来进行多次迭代寻优。其中，改进 的粒子群速度和位置更新所遵循的规则如下：
[0036]
[0037]
[0038] 设搜寻空间为D维的，粒子群包含N个粒子。式中
和
分别代表第i个粒子在时刻iter的速度和位置矢量。α是用于 控制和限制速度的收缩因子，(^和(：2分别为认知和社会系数(一般均可以设置为2) 41;1和 r2>1是在[0，1]范围内取值的两个独立的随机数。抑纽广代表第i个粒子截至iter时刻到历 史位置最优值。Gb estlte^代表所有粒子中截至到iter时刻的历史位置最优值。w(L)是动态 惯性权重，用于调节搜索空间范围，以期获得全局最优，其计算公式为：
[0039]
[0040] 式中，Wstart是初始惯性权重，Wend是迭代至最大次数时的惯性权重，L为粒子群当前 寻优的迭代次数，T max为粒子群寻优的迭代次数上限，一般来说惯性权值取wstart = 0.9，We3nd = 0.4时粒子群寻优效果较好。
[0041 ] (10)、若迭代次数L>M，则粒子群寻优迭代结束，此时把对k时刻状态估计值的寻 优结果Gbest1作为k时刻的状态估计值，然后对下一时刻状态进行估计；
[0042] (11)、若k = k+l < S，则迭代继续，若k = k+l>S，则迭代结束，输出辨识结果。
【附图说明】
[0043] 图1为本发明实施例的方法流程图；
[0044] 图2实施例的动态振荡信号；
[0045] 图3为实施例采用本发明所提方法的信号频率w辨识结果；
[0046]图4为实施例采用本发明所提方法的信号阻尼因子δ辨识结果。
【具体实施方式】
[0047]下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明 而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价 形式的修改均落于本申请所附权利要求所限定的范围。
[0048]如图1所示，动态振荡信号参数辨识方法，其包含如下步骤：
[0049] (1)、获取状态分量中包含动态振荡信号待辨识参数的状态空间表达式。
[0050] (2)、初始化，包括:设定参数辨识的初值4，初始参数辨识误差协方差^，以及过程 噪声和量测噪声所满足的协方差矩阵Q和R，整体算法迭代次数最大值S和粒子群寻优最大 迭代次数Μ。
[0051] (3)、由已经得到的k-Ι时刻的状态估计值和状态估计误差协方差，利用扩展卡尔 曼滤波的预测步，得到k时刻的状态预测值和状态预测误差协方差，计算公式为：
[0052]
[0053]
[0054] 式中，毛表示k时刻的状态预测值，f( ·)对应具体问题状态方程中的非线性函数， 先-1表示k-Ι时刻的状态估计向量，Uk-i表示k-1时刻的控制输入。戶i表示k时刻的状态预测误 差协方差，
表示非线性函数f( ·)在馬4处的雅克比矩阵，表示k-Ι时刻 的状态估计误差协方差，上标T表示转置，Qh是系统噪声k-Ι时刻所满足的协方差矩阵。
[0055] (4)、在上一步基础上，利用扩展卡尔曼滤波的滤波步，得到k时刻的状态估计值。
[0056] a?
[0057]
[0058]
[0059]
[0060] 式中，Kk表示k时刻的卡尔曼滤波增益，&表示k时刻的状态预测误差协方差，上标 τ表示转置
表示非线性函数h( ·)在％处的雅克比矩阵，其中h( ·)对应具体 问题输出方程中的非线性函数。Rk是量测噪声k时刻所满足的协方差矩阵，&表示k时刻的 状态估计误差协方差，I是和状态向量维度相同的单位矩阵，氧表示k时刻的状态估计向量， yk是k时刻输出方程的输出量。
[0061] (5)、判定k时刻的参数辨识结果是否满足相应的实际约束条件。若满足，则直接运 用EKF再次迭代辨识。
[0062] (6)、若不满足，则运用改进的约束算法对该时刻参数辨识结果进行约束，把约束 的状态估计问题等价转化为约束的优化问题。
[0063] 约束优化目标函数为：
[0064]
[0065] 式中，f4:表示k时刻待求寻优状态估计值，yk是k时刻输出量，Wi和W2是正定的权矩 阵，为方便起见，本发明实施例中选取为单位矩阵，式中xk满足约束条件如下式：
[0066] Dx<d
[0067] D是已知的sXn常数行满秩矩阵，s是受约束条件参数的个数，η是状态向量的维 数，显然，s < η，d是已知的常数约束条件。
[0068] 在约束优化目标函数中，对k时刻的估计输出量λ进行求解时，为了充分利用已知 的量测数据信息，采用卡尔曼平滑方法，计算公式如下：
[0069]
[0070]
[0071]
[0072] 式中Λ,η表示k-Ι时刻的平滑状态估计向量，L是卡尔曼平滑增益，Α是k时刻的 估计输出量。
[0073] (7)、在上一步的基础上，借助罚函数方法，通过对原目标函数加上一个惩罚项，把 约束优化问题转化为一个无约束优化问题。
[0074]转换的计算公式如下：
[0075] F(x)=f(x)+h(k)H(x),xeRn
[0076] 式中，f(x)是原目标函数，h(k)是动态更新的惩罚值，一般情况下或者
k是当前的迭代次数。H(x)是惩罚因子，计算公式如下：
[0077]
[0078]式中，0(qi(x))是多级分配函数，qi(x)是与违反约束条件有关的函数， qi(x)=max {0，gi(x)}，i = l，···，s，其中gi(x)=Dx_d。γ (qi(x))表示罚函数的效力。函数取值遵循的规 则为：
[0079]