CN103886193A - 一种电力系统模糊自适应抗差估计方法 - Google Patents

Abstract

本发明公开了一种电力系统模糊自适应抗差估计方法，本发明提出一种计及量测权重不确定性的电力系统模糊自适应抗差估计算法(FARE)，该算法提出了测点劣质性的模糊隶属度概念，并根据测点的隶属度在线修正测点的标准差，实现了对粗差的自适应。以现代内点法求解，收敛性好，结果不易受初值影响，适用于求解大规模优化问题。有效的地辨识了不良数据，避免了残差污染及残差淹没，对状态量及量测量有着极佳地抗差估计性能，具有很好地工程运用背景。

Description

一种电力系统模糊自适应抗差估计方法

技术领域

本发明涉及一种电力系统模糊自适应抗差估计方法，属于电力系统监测、分析和控制技术领域。 

背景技术

状态估计是电力系统分析与控制的基础，是能量管理系统核心基础模块，其主要任务是根据数据采集与监控系统(SCADA)提供的实时信息，估计出电网各母线电压(幅值、相角)和功率，并且包含不良数据检测与辨识功能。 

传统的状态估计算法由Schweppe等人于1970年创立。根据求解目标函数的不同，状态估计器包括：加权最小二乘(WLS)、非二次准则(QC、QL等)、加权最小绝对值(WLAV)、最小中位数平方(LMS)、最小截平方(LTS)等。目前应用最广泛的是WLS估计准则，优点在于模型简单、计算量小、对理想的正态分布量测量，能达到最优估计结果，但估计结果容易受量测不良数据影响，需专门的不良数据检测与辨识程序(如基于最大加权残差准则)以减小不良数据对状态估计结果的影响，当量测存在大规模不良数据时，迭代次数多，且未必能有效辨识出所有不良数据。QC、QL、WLAV、LMS、LTS等估计器均属于抗差估计范畴，其优点在于状态估计算法无需额外的不良数据检测与辨识程序，但也存在需主观设定加权因子、计算量大、某些特定情况下抗差性能薄弱等缺陷。 

目前国内外绝大多数状态估计算法均是假设量测方差σ²已知，而在实际工程中，大规模电力系统中存在众多量测仪表，对量测仪表进行重新校核必须付出很大的经济代价，因而很难系统地重新校核量测仪表；此外随着设备的老化、运行环境的时刻变化，仪表的量测精度难以持续保持稳定。也就是说，量测方差是随时间变化的，难以精确估计的，即量测权重存在不确定性。 

发明内容

发明目的：本发明提出一种电力系统模糊自适应抗差估计方法，有效地辨识了不良数据，避免了残差污染及残差淹没，有着极佳地抗差估计性能。 

技术方案：本发明采用的技术方案为一种电力系统模糊自适应抗差估计方法，包括以下步骤： 

1）定义第i个测点劣质性的模糊隶属度函数为： 

$v_i\left(\right|r_i|,{\mathrm{\σ}}_i)=1-1/(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\σ}}_i}\right)}^b)\∀i=\mathrm{1,2},\mathrm{K\; m}$

式中：r_i为测点i的残差，σ_i为测点i的量测标准差，v_i(|r_i|,σ_i)为测点i的模糊隶属度，a,b为大于0的常数； 

2）以最小化测点劣质性的加权模糊隶属度之和为优化目标，提出以下优化模型： 

$\left\{\begin{array}{c}\min J=\underset{i=1}{\overset{i=m}{\mathrm{\Σ}}}\frac{1}{{{\mathrm{\σ}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\σ}}_i}\right)}^b))\\ s.t.r=z-h\left(x\right)\end{array}\right.$

式中，

为测点i的权重； 

3）考虑量测权重不确定性，即量测标准差的不确定性，基于测点模糊隶属度修正量测标准差，以实现对量测粗差的自适应，对于第k+1次迭代，令： 

${{\mathrm{\σ}}_i}^{(k+1)}={{\mathrm{\σ}}_i}^{\left(k\right)}\·f\left(v_i\right)\∀i=\mathrm{1,2},\mathrm{K\; m}$

4）引入非负松弛因子l,u，则步骤2)中优化模型可等价为： 

$\left\{\begin{array}{c}\min J(l,u)=\underset{i=1}{\overset{i=m}{\mathrm{\Σ}}}\frac{1}{{{\mathrm{\σ}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\σ}}_i}\right)}^b))\\ s.t.z-h\left(x\right)+l-u=0(l,u)>0\end{array}\right.$

5）将步骤3）中的等式约束设为障碍函数，可得以下拉格朗日函数： 

$L(l,u,x,\mathrm{\α},\mathrm{\β},\mathrm{\λ})=\underset{i=1}{\overset{i=m}{\mathrm{\Σ}}}\frac{1}{{{\mathrm{\σ}}_i}^2}/(1-(1+{\left(\frac{|l_i+u_i|}{a{\mathrm{\σ}}_i}\right)}^b))+\mathrm{\λ}(z-h\left(x\right)+l-u)+{\mathrm{\α}}^{T}l+{\mathrm{\β}}^{T}u$

式中：λ、α、β为m维拉格朗日乘子，即对偶变量；l,u为原变量； 

6）求解步骤4）中拉格朗日函数的KKT条件得： 

$\left\{\begin{array}{c}{L}_l=\frac{\∂J(l,u)}{\∂l}+\mathrm{\λ}+\mathrm{\α}=0\\ L_u=\frac{\∂J(l,u)}{\∂l}-\mathrm{\λ}+\mathrm{\β}=0\\ L_x=\▿h\left(x\right)\mathrm{\λ}=0\\ L_{\mathrm{\α}}=l\⇒{L_{\mathrm{\α}}}^u=\mathrm{ALe}-\mathrm{\μe}=0\\ L_{\mathrm{\β}}=u\⇒{L_{\mathrm{\β}}}^u=\mathrm{BUe}-\mathrm{\μe}=0\\ L_{\mathrm{\λ}}=z-h\left(x\right)+l-u=0\end{array}\right.$

上式中，▽h(x)为h(x)的雅克比矩阵，L＝diag(l₁,l₂,…,l_m)，U＝diag(u₁,u₂…,u_m)；A＝diag(α₁,α₂,…,α_m)，B＝diag(β₁,β₂,…,β_m)，e＝[1,…1]^T，μ为扰动因子，且满足

8）将步骤5）中的KKT条件用牛顿—拉弗森法线性化后得到如下修正方程： 

$\left[\begin{array}{cccccc}\frac{{\∂}^{2}J(l,u)}{\∂l^2}& \frac{{\∂}^{2}J(l,u)}{\∂l\∂u}& I& I& 0& 0\\ \frac{{\∂}^{2}J(l,u)}{\∂u\∂l}& \frac{{\∂}^{2}J(l,u)}{\∂u^2}& -I& 0& I& 0\\ I& -I& 0& 0& 0& -{\▿}^{T}h\left(x\right)\\ A& 0& 0& L& 0& 0\\ 0& B& 0& 0& U& 0\\ 0& 0& \▿h\left(x\right)& 0& 0& {\▿}^{2}h\left(x\right)\mathrm{\λ}\end{array}\right]\·\left[\begin{array}{c}\mathrm{dl}\\ \mathrm{du}\\ \mathrm{dx}\\ \mathrm{d\α}\\ \mathrm{d\β}\\ \mathrm{d\λ}\end{array}\right]=-\left[\begin{array}{c}{L}_l\\ L_u\\ L_x\\ L_{\mathrm{\α}}\\ L_{\mathrm{\β}}\\ L_{\mathrm{\λ}}\end{array}\right]$

式中：▽²h(x)为h(x)的海森矩阵； 

8）根据原—对偶内点法的迭代步骤，更新状态量，直至收敛。 

有益效果：本发明提出一种计及量测权重不确定性的电力系统模糊自适应抗差估计算法(FARE)，该算法提出了测点劣质性的模糊隶属度概念，并根据测点的隶属度在线修正测点的标准差，实现了对粗差的自适应。以现代内点法求解，收敛性好，结果不易受初值影响，适用于求解大规模优化问题。有效的地辨识了不良数据，避免了残差污染及残差淹没，对状态量及量测量有着极佳地抗差估计性能，具有很好地工程运用背景。 

附图说明

图1为本发明的工作流程图； 

图2为以IEEE14为标准测试节点，本发明方法FARE与WLAV、QC、WLS算法 量测量抗差性能的测试效果图，量测的加权残差均值|r/σ|取值分布图； 

图3为以IEEE118为标准测试节点，本发明方法FARE与WLAV、QC、WLS算法量测量抗差性能的测试效果图，量测的加权残差均值|r/σ|取值分布图。 

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等同形式的修改均落于本申请所附权利要求所限定的范围。 

1.建立电力系统模糊自适应抗差估计模型 

电力系统状态估计量测方程为： 

z＝h(x)+e 

式中：x为状态向量，设系统节点数为n，则x包含节点电压的幅值和相角，其维数为2n-1；z为m维量测向量；h(x)为量测函数向量；e为量测误差向量。 

状态估计的残差方程为： 

$r=z-h\left(\hat{x}\right)$

式中：r为m维残差向量，

为状态量x的估计值。 

则对于第i个测点，若其估计残差r_i很小，则可认为该测点为优质测点；反之，若其估计残差r_i很大，则可认为该测点为劣质测点。由于各个测点的量测精度不尽相同，定义σ_i为测点i的量测标准差，因而以测点的加权残差r_i/σ_i的大小区分优劣测点更为合理。 

定义事件y为测点的优劣，若集合Y包含了所有事件y，那么元素y与集合Y的关系可用一个特征函数——隶属度函数v(y)表示，对于经典的数据集合理论，有： 

但实际测点的优劣只是相对的概念，即不存在绝对的优与劣，相比于经典的数据集合理论，模糊集合允许隶属度取[0,1]上的任何值，本文选取连续可微的钟形隶属度函数，对于第i个测点，其劣质性的模糊隶属度函数为： 

$v_i\left(\right|r_i|,{\mathrm{\σ}}_i)=1-1/(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\σ}}_i}\right)}^b)\∀i=\mathrm{1,2},\mathrm{K\; m}$

式中：a、b为大于0的模糊隶属度特征参数。 

以加权最小二乘(WLS)为主的电力系统状态估计算法以最小加权残差平方和为优化目标，类似地，本文基于测点劣质性的模糊隶属度，以最小化测点劣质性的加权模糊隶属度之和为优化目标，提出以下优化模型： 

$\left\{\begin{array}{c}\min J=\underset{i=1}{\overset{i=m}{\mathrm{\Σ}}}\frac{1}{{{\mathrm{\σ}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\σ}}_i}\right)}^b))\\ s.t.r=z-h\left(x\right)\end{array}\right.$

式中：

为测点i的权重。 

目前国内外绝大多数状态估计算法均是假设量测方差σ²已知，而在实际工程中，大规模电力系统中存在众多量测仪表，对量测仪表进行重新校核必须付出很大的经济代价，因而很难系统地重新校核量测仪表；此外由于使用时间的增加、时刻变化的外部环境，仪表的量测精度难以持续保持稳定。也就是说，量测方差是难以精确估计的，即量测权重存在不确定性。国内实际系统的状态估计一般主观设置各量测的权重因子，调试及维护复杂，且主观设定的权重因子未必与实际量测方差吻合。 

为此，本文基于测点模糊隶属度修正量测标准差，即对于求解优化目标过程中的第k+1次迭代，令： 

${{\mathrm{\σ}}_i}^{(k+1)}={{\mathrm{\σ}}_i}^{\left(k\right)}\·f\left(v_i\right)\∀i=\mathrm{1,2},\mathrm{K\; m}$

量测标准差修正函数f(v)应遵循2个原则：①当测点模糊隶属度v_i接近0时，此时该测点为优质测点，应保持σ_i近似不变；②当测点模糊隶属度v_i接近1时，此时该测点为劣质测点，应增大σ_i，即出现量测粗差时，应降低该测点在迭代中的权重，减小其对状态估计结果的影响。 

基于上述准则，可选取修正函数f(v)＝[1-(1/v)]^1/b。 

2.基于原对偶内点法的FARE模型的求解 

参照图1，本发明首先获取电力系统的网络参数、拓扑结构、量测参数，网 络参数包括：母线编号，节点并联电容电抗，支路首端编号和末端编号，支路电阻电抗，电抗，对地充电电容，变压器电阻、电抗、标准变比；拓扑结构主要包括连接两个电气元件间的开关状态；量测参数包括：节点注入有功、无功量测，支路首端及末端有功、无功量测，以及母线电压幅值量测。 

在获取状态估计程序需要的数据后，进行程序初始化，包括：状态量电压幅值与相角、拉格朗日乘子、原对偶松弛变量、原对偶中心参数、最大迭代次数以及收敛精度。并形成节点导纳矩阵。 

引入m维非负松弛因子l,u，则优化目标可等价为： 

$\left\{\begin{array}{c}\min J(l,u)=\underset{i=1}{\overset{i=m}{\mathrm{\Σ}}}\frac{1}{{{\mathrm{\σ}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\σ}}_i}\right)}^b))\\ s.t.z-h\left(x\right)+l-u=0(l,u)>0\end{array}\right.$

将上式中的等式约束设为障碍函数，可得以下拉格朗日函数： 

$L(l,u,x,\mathrm{\α},\mathrm{\β},\mathrm{\λ})=\underset{i=1}{\overset{i=m}{\mathrm{\Σ}}}\frac{1}{{{\mathrm{\σ}}_i}^2}/(1-(1+{\left(\frac{|l_i+u_i|}{a{\mathrm{\σ}}_i}\right)}^b))+\mathrm{\λ}(z-h\left(x\right)+l-u)+{\mathrm{\α}}^{T}l+{\mathrm{\β}}^{T}u$

式中：λ、α、β为m维拉格朗日乘子，即对偶变量；l,u为原变量。 

求解其KKT条件可得： 

$L_l=\frac{\∂J(l,u)}{\∂l}+\mathrm{\λ}+\mathrm{\α}=0$

$L_u=\frac{\∂J(l,u)}{\∂l}-\mathrm{\λ}+\mathrm{\β}=0$

L_x＝▽h(x)λ＝0 

$L_{\mathrm{\α}}=l\⇒{L_{\mathrm{\α}}}^u=\mathrm{ALe}-\mathrm{\μe}=0$

$L_{\mathrm{\β}}=u\⇒{L_{\mathrm{\β}}}^u=\mathrm{BUe}-\mathrm{\μe}=0$

L_λ＝z-h(x)+l-u＝0 

上式中，▽h(x)为h(x)的雅克比矩阵，L＝diag(l₁,l₂,…,l_m)，U＝diag(u₁,u₂…,u_m)；A＝diag(α₁,α₂,…,α_m)，B＝diag(β₁,β₂,…,β_m)， e＝[1,…1]^T，μ为扰动因子，且满足

KKT条件是非线性方程组，可用牛顿—拉弗森法求解，将其线性化后得到修正方程组： 

$\frac{{\∂}^{2}J(l,u)}{\∂l^2}\mathrm{dl}+\frac{{\∂}^{2}J(l,u)}{\∂l\∂u}\mathrm{du}+\mathrm{d\λ}+\mathrm{d\α}=-L_l$

$\frac{{\∂}^{2}J(l,u)}{\∂l\∂u}\mathrm{dl}+\frac{{\∂}^{2}J(l,u)}{\∂u^2}\mathrm{du}-\mathrm{d\λ}+\mathrm{d\β}=-L_u$

▽²h(x)λdx+▽h(x)dλ＝-L_x

-▽^Th(x)dx+dl-du＝-L_λ

$\mathrm{Ld\α}+\mathrm{Adl}=-L_{\mathrm{\α}}^u$

$\mathrm{Ud\β}+\mathrm{Bdu}=-L_{\mathrm{\β}}^u$

式中：▽²h(x)为h(x)的海森矩阵。写成矩阵形式可得： 

$\left[\begin{array}{cccccc}\frac{{\∂}^{2}J(l,u)}{\∂l^2}& \frac{{\∂}^{2}J(l,u)}{\∂l\∂u}& I& I& 0& 0\\ \frac{{\∂}^{2}J(l,u)}{\∂u\∂l}& \frac{{\∂}^{2}J(l,u)}{\∂u^2}& -I& 0& I& 0\\ I& -I& 0& 0& 0& -{\▿}^{T}h\left(x\right)\\ A& 0& 0& L& 0& 0\\ 0& B& 0& 0& U& 0\\ 0& 0& \▿h\left(x\right)& 0& 0& {\▿}^{2}h\left(x\right)\mathrm{\λ}\end{array}\right]\·\left[\begin{array}{c}\mathrm{dl}\\ \mathrm{du}\\ \mathrm{dx}\\ \mathrm{d\α}\\ \mathrm{d\β}\\ \mathrm{d\λ}\end{array}\right]=-\left[\begin{array}{c}{L}_l\\ L_u\\ L_x\\ L_{\mathrm{\α}}\\ L_{\mathrm{\β}}\\ L_{\mathrm{\λ}}\end{array}\right]$

式中： $\frac{{\∂}^{2}J(l,u)}{\∂l\∂u}=\frac{{\∂}^{2}J(l,u)}{\∂u\∂l},\frac{{\∂}^{2}J(l,u)}{\∂u^2}=\frac{{\∂}^{2}J(l,u)}{\∂l^2},$ 且由上式可计算[dl du dx]^T、[dα dβ dλ]^T。 

计算对偶间隙Gap＝α^Tl+β^Tu及扰动因子μ＝δ·Gap/2m。此外，为了保证松弛因子恒为非负，依据下式原、对偶变量的迭代步长(θ_r,θ_D): 

${\mathrm{\&theta;}}_r=0.9995\min \{\underset{i}{\min }(-\frac{l_i}{{\mathrm{dl}}_i}:{\mathrm{dl}}_i<0;-\frac{u_i}{{\mathrm{du}}_i}:{\mathrm{du}}_i<0),1\}$

${\mathrm{\&theta;}}_D=0.9995\min \{\underset{i}{\min }(-\frac{{\mathrm{\&alpha;}}_i}{{\mathrm{d\&alpha;}}_i}:{\mathrm{d\&alpha;}}_i<0;-\frac{{\mathrm{\&beta;}}_i}{{\mathrm{d\&beta;}}_i}:{\mathrm{d\&beta;}}_i<0),1\}$

修正原、对偶变量： 

$\left[\begin{array}{c}{l}^{(k+1)}\\ u^{(k+1)}\\ x^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{l}^{\left(k\right)}\\ u^{\left(k\right)}\\ x^{\left(k\right)}\end{array}\right]+{\mathrm{\&theta;}}_r\left[\begin{array}{c}\mathrm{dl}\\ \mathrm{du}\\ \mathrm{dx}\end{array}\right],\left[\begin{array}{c}{\mathrm{\&lambda;}}^{(k+1)}\\ {\mathrm{\&alpha;}}^{(k+1)}\\ {\mathrm{\&beta;}}^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&lambda;}}^{\left(k\right)}\\ {\mathrm{\&alpha;}}^{\left(k\right)}\\ {\mathrm{\&beta;}}^{\left(k\right)}\end{array}\right]+{\mathrm{\&theta;}}_r\left[\begin{array}{c}\mathrm{d\&lambda;}\\ \mathrm{d\&alpha;}\\ \mathrm{d\&beta;}\end{array}\right]$

修正量测标准差：σ^(k+1)＝σ^(k)·(1/v)^1/b。 

重复上述修正过程直至对偶间隙Gap达收敛精度或迭代次数超过设定的最大迭代次数。 

3.算例分析 

本发明的测试算例包括IEEE14、IEEE30、IEEE57、IEEE118节点，波兰2383,2746节点(以WP-2383，WP-2746表示)，量测数据由严格潮流结果添加随机量测误差得到，考虑量测标准差的不确定性，电压量测标准差服从[0.002,0.005]的均匀分布，功率量测标准差服从[0.004,0.01]的均匀分布。量测冗余度

介于3.5～4.5之间。不良数据在量测数据的基础上加减[5,20]倍的最大量测标准差，且随机添加10%不良数据。 

为验证本发明提出的FARE算法的优越性，将FARE与以往提出的的状态估计算法，包括WLS、WLAV、QC，从不良数据辨识、抗差性能、收敛性能等方面进行对比研究。WLS、WLAV、QC算法的电压量测权重取4，有功、无功量测权重取1。对于FARE算法的量测标准差初值σ⁽⁰⁾，电压量测取0.001，有功、无功量测取0.002；模糊隶属度参数取a=2.5，b=3。 

3.1量测量抗差性能 

为比较各种算法的不良数据辨识以及抵御残差污染淹没的性能，在量测真值 

量测标准差真值

已知的情况下，定义： 

${\mathrm{\&tau;}}_i=|\widehat{z_i}-{\stackrel{\&OverBar;}{z}}_i|/{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_i$

式中：

为第i个测点的估计值，指标τ衡量了量测估计值偏离量测真值的程度。 

以IEEE14、IEEE118为测试算例，4种状态估计算法下τ的百分比分布见说明书附图2、图3。 

由图2、图3可知，WLS算法估计下的τ＞3的百分比在20%以上，而量测不良数据比例仅为10%，这说明出现了很严重的残差污染淹没现象，因而基于 WLS的QC算法也很难有效辨识不良数据及抵御残差污染淹没。相比于WLAV、QC、WLS算法，FARE算法估计下τ的分布主要集中于0～1之间，且τ＞3的比例明显少于其余3种算法，因而FARE算法下的量测估计值更接近于量测真值，有着很好地不良数据辨识及抵御残差污染淹没性能。 

3.2状态量抗差性能 

为了比较4种算法的状态量抗差性能，采用以下2个指标来评价不同状态估计算法对状态量的估计性能: 

1)电压幅值的均方误差： 

$\mathrm{MSE}1=\underset{i=1}{\overset{i=n}{\mathrm{\&Sigma;}}}{(V_i-\stackrel{}{\widehat{V_i}})}^2/n$

式中：V_i为节点i电压幅值的真值，为节点i电压幅值的估计值。 

2)电压相角的均方误差： 

$\mathrm{MSE}2=\underset{i=2}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&theta;}}_i-{\widehat{\mathrm{\&theta;}}}_i)}^2/(n-1)$

式中：θ_i为节点i电压相角的真值，

为节点i相角的估计值，默认节点1为平衡节点。 

在6个系统系统均含10%不良数据的情况下，表1给出了4种状态估计算法的MSE1、MSE2的指标。 

比较表1中MSE1和MSE2的大小，可知各算法对状态量的估计性能：FARE>WLAV>QC>WLS，且FARE明显优于其余3种算法。此外，FARE算法对电压相角的估计性能要好于对电压幅值的估计性能。 

综合量测量和状态量的抗差性能可知，在选取的4种状态估计算法中，本发明提出的FARE是最优质的电力系统状态估计算法。 

表1：含10%不良数据下不同估计算法的MSE1及MSE2指标性能 

Claims (1)

1.一种电力系统模糊自适应抗差估计方法，其特征在于，包括以下步骤：

1）定义第i个测点劣质性的模糊隶属度函数为：

$v_i\left(\right|r_i|,{\mathrm{\&sigma;}}_i)=1-1/(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\&sigma;}}_i}\right)}^b)\&ForAll;i=\mathrm{1,2},\mathrm{K\; m}$

式中：r_i为测点i的残差，σ_i为测点i的量测标准差，v_i(|r_i|,σ_i)为测点i的模糊隶属度，a,b为大于0的常数；

2）以最小化测点劣质性的加权模糊隶属度之和为优化目标，提出以下优化模型：

$\left\{\begin{array}{c}\min J=\underset{i=1}{\overset{i=m}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\&sigma;}}_i}\right)}^b))\\ s.t.r=z-h\left(x\right)\end{array}\right.$

式中，

为测点i的权重；

3）考虑量测权重不确定性，即量测标准差的不确定性，基于测点模糊隶属度修正量测标准差，以实现对量测粗差的自适应，对于第k+1次迭代，令：

${{\mathrm{\&sigma;}}_i}^{(k+1)}={{\mathrm{\&sigma;}}_i}^{\left(k\right)}\&CenterDot;f\left(v_i\right)\&ForAll;i=\mathrm{1,2},\mathrm{K\; m}$

4）引入非负松弛因子l,u，则步骤2)中优化模型可等价为：

$\left\{\begin{array}{c}\min J(l,u)=\underset{i=1}{\overset{i=m}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\&sigma;}}_i}\right)}^b))\\ s.t.z-h\left(x\right)+l-u=0(l,u)>0\end{array}\right.$

5）将步骤3）中的等式约束设为障碍函数，可得以下拉格朗日函数：

$L(l,u,x,\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&lambda;})=\underset{i=1}{\overset{i=m}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_i}^2}/(1-(1+{\left(\frac{|l_i+u_i|}{a{\mathrm{\&sigma;}}_i}\right)}^b))+\mathrm{\&lambda;}(z-h\left(x\right)+l-u)+{\mathrm{\&alpha;}}^{T}l+{\mathrm{\&beta;}}^{T}u$

式中：λ、α、β为m维拉格朗日乘子，即对偶变量；l,u为原变量；

6）求解步骤4）中拉格朗日函数的KKT条件得：

$\left\{\begin{array}{c}{L}_l=\frac{\&PartialD;J(l,u)}{\&PartialD;l}+\mathrm{\&lambda;}+\mathrm{\&alpha;}=0\\ L_u=\frac{\&PartialD;J(l,u)}{\&PartialD;l}-\mathrm{\&lambda;}+\mathrm{\&beta;}=0\\ L_x=\&dtri;h\left(x\right)\mathrm{\&lambda;}=0\\ L_{\mathrm{\&alpha;}}=l\&DoubleRightArrow;{L_{\mathrm{\&alpha;}}}^u=\mathrm{ALe}-\mathrm{\&mu;e}=0\\ L_{\mathrm{\&beta;}}=u\&DoubleRightArrow;{L_{\mathrm{\&beta;}}}^u=\mathrm{BUe}-\mathrm{\&mu;e}=0\\ L_{\mathrm{\&lambda;}}=z-h\left(x\right)+l-u=0\end{array}\right.$

上式中，▽h(x)为h(x)的雅克比矩阵，L＝diag(l₁,l₂,…,l_m)，U＝diag(u₁,u₂…,u_m)；A＝diag(α₁,α₂,…,α_m)，B＝diag(β₁,β₂,…,β_m)，e＝[1,…1]^T，μ为扰动因子，且满足

7）将步骤5）中的KKT条件用牛顿—拉弗森法线性化后得到如下修正方程：

$\left[\begin{array}{cccccc}\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;l^2}& \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;l\&PartialD;u}& I& I& 0& 0\\ \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;u\&PartialD;l}& \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;u^2}& -I& 0& I& 0\\ I& -I& 0& 0& 0& -{\&dtri;}^{T}h\left(x\right)\\ A& 0& 0& L& 0& 0\\ 0& B& 0& 0& U& 0\\ 0& 0& \&dtri;h\left(x\right)& 0& 0& {\&dtri;}^{2}h\left(x\right)\mathrm{\&lambda;}\end{array}\right]\&CenterDot;\left[\begin{array}{c}\mathrm{dl}\\ \mathrm{du}\\ \mathrm{dx}\\ \mathrm{d\&alpha;}\\ \mathrm{d\&beta;}\\ \mathrm{d\&lambda;}\end{array}\right]=-\left[\begin{array}{c}{L}_l\\ L_u\\ L_x\\ L_{\mathrm{\&alpha;}}\\ L_{\mathrm{\&beta;}}\\ L_{\mathrm{\&lambda;}}\end{array}\right]$

式中：▽²h(x)为h(x)的海森矩阵；

8）根据原—对偶内点法的迭代步骤，更新状态量，直至收敛。

Images