CN113810024B - 一种基于混合概率选择算子的代价参考粒子滤波方法 - Google Patents

一种基于混合概率选择算子的代价参考粒子滤波方法 Download PDF

Info

Publication number
CN113810024B
CN113810024B CN202111017989.3A CN202111017989A CN113810024B CN 113810024 B CN113810024 B CN 113810024B CN 202111017989 A CN202111017989 A CN 202111017989A CN 113810024 B CN113810024 B CN 113810024B
Authority
CN
China
Prior art keywords
particle
cost
particles
value
subset
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202111017989.3A
Other languages
English (en)
Other versions
CN113810024A (zh
Inventor
张新雨
王晓璇
雷碧玉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian University of Technology
Original Assignee
Xian University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian University of Technology filed Critical Xian University of Technology
Priority to CN202111017989.3A priority Critical patent/CN113810024B/zh
Publication of CN113810024A publication Critical patent/CN113810024A/zh
Application granted granted Critical
Publication of CN113810024B publication Critical patent/CN113810024B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H17/00—Networks using digital techniques
    • H03H17/02—Frequency selective networks
    • H03H17/0248—Filters characterised by a particular frequency response or filtering method
    • H03H17/0255—Filters based on statistics
    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H17/00—Networks using digital techniques
    • H03H17/02—Frequency selective networks
    • H03H17/0248—Filters characterised by a particular frequency response or filtering method
    • H03H17/0261—Non linear filters

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Mathematical Physics (AREA)
  • Nonlinear Science (AREA)
  • Probability & Statistics with Applications (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于混合概率选择算子的CRPF方法，首先随机产生N个实数作为系统状态的样本，所述样本称为粒子，并定义每一粒子对应的代价值均为0，由所述样本与代价值组成的集合称为粒子—代价集合；同时，分别随机产生粒子组成1号子集合与2号子集合，每一个子集合的数量为N；分别计算所有子集合中粒子的代价函数值、风险函数值及权值；对子集合进行重采样、更新、信息交互、粒子选择操作，得到最终的粒子—代价集合并对所得到的集合进行更新；最后系统滤波，最终得到系统状态的最优估计值。本发明解决了现有技术中对噪声统计特性未知系统滤波准确度较低的问题。

Description

一种基于混合概率选择算子的代价参考粒子滤波方法
技术领域
本发明属于非线性滤波技术领域，具体涉及一种基于混合概率选择算子的代价参考粒子滤波方法。
背景技术
在很多的工业控制过程中，所有通过传感器测量得到的数据均含有噪声且不能被完全消除，同时有很多内部的系统状态是无法直接测量的。因此，在系统状态分析过程中，需要对传感器所测得的含有大量噪声的数据进行滤波处理，进而尽可能得到被测量的最优值或无法测得的系统状态最优估计值。例如在非线性、非高斯的锂离子电池寿命预测中，锂电池的剩余寿命无法直接进行在线测量，只能通过在线测量的充电电压、电流等相关参量对其进行估计。因为现代滤波技术具有预测和估计的作用，因此在此领域具有较大优势。代价参考粒子滤波(Cost Reference Particle Filter,CRPF)能够处理非线性和噪声统计特性未知的滤波问题，但标准CRPF使用多项式重采样、随机重采样、系统重采样等传统重采样方法进行重采样。这些重采样方法虽然可以增加有效粒子数量，但是由于一些权值大的粒子被大量复制，使得粒子过于集中于同一区域，导致重采样结果出现粒子多样性匮乏这一情况。故标准CRPF的重采样过程存在着粒子多样性匮乏的问题，进而影响滤波的精度与准确度。因此，需要对标准CRPF做进一步地优化，改善CRPF的滤波精度与准确度，进而提高对锂离子剩余寿命估计的精度。
发明内容
本发明的目的是提供一种基于混合概率选择算子的代价参考粒子滤波方法，解决了现有技术中对噪声统计特性未知系统滤波准确度较低的问题。
本发明所采用的技术方案是，一种基于混合概率选择算子的代价参考粒子滤波方法，具体按照以下步骤实施：
步骤1、随机产生N个实数作为系统状态的样本，所述样本称为粒子，并定义每一粒子对应的代价值均为0，由所述样本与代价值组成的集合称为粒子—代价集合；同时，分别随机产生粒子组成1号子集合与2号子集合，每一个子集合的数量为N；
步骤2、分别计算所有子集合中粒子的代价函数值、风险函数值及权值；
步骤3、重采样子集合，重采样方法选取多项式重采样；
步骤4、更新子集合；
步骤5、子集合之间粒子的信息交互，若到达设定的交换步数，则执行此项操作；否则，跳过该步骤，执行步骤6；
步骤6、生成最终用于系统滤波的粒子—代价集合；
步骤7、更新最终的粒子—代价集合；
步骤8、系统滤波，最终得到系统状态的最优估计值。
本发明的特点还在于，
步骤1生成初始规模为N的粒子—代价集合E0，表示为
Figure GDA0004266232870000021
其中，/>
Figure GDA0004266232870000022
表示系统状态随机估计值，/>
Figure GDA0004266232870000023
即/>
Figure GDA0004266232870000024
服从均匀分布U(I0)，初始代价函数值/>
Figure GDA0004266232870000025
i＝1,2,…,N，N为系统状态随机估计值的采样样本数；
同时，生成两个大小为N的子集合
Figure GDA0004266232870000031
其中，/>
Figure GDA0004266232870000032
服从均匀分布U(I0)，j为子集合个数，j＝1,2，当j＝1时，/>
Figure GDA0004266232870000033
表示1号子集合，当j＝2时，/>
Figure GDA0004266232870000034
表示2号子集合。
步骤2具体如下：
步骤2.1、计算t时刻子集合中粒子
Figure GDA0004266232870000035
的代价函数值/>
Figure GDA0004266232870000036
j＝1,2，i＝1,2,…,N，t＝1,2,…,T，T为时间序列长度，代价函数的计算方式如下：
Figure GDA0004266232870000037
Figure GDA0004266232870000038
在公式(1)、(2)中，
Figure GDA0004266232870000039
表示t时刻j号子集合中第i个粒子/>
Figure GDA00042662328700000310
的代价函数值；λ表示遗忘因子，0＜λ＜1；/>
Figure GDA00042662328700000311
表示t-1时刻j号子集合中第i个粒子/>
Figure GDA00042662328700000312
的代价值函数值；/>
Figure GDA00042662328700000313
表示粒子/>
Figure GDA00042662328700000314
的代价增量，/>
Figure GDA00042662328700000315
表示粒子/>
Figure GDA00042662328700000316
的代价增量函数；yt表示t时刻传感器的测量值，h(·)表示传感器估计值的观测函数，预设参数q满足q≥1；
步骤2.2、由子集合粒子的代价函数值
Figure GDA00042662328700000317
计算得到风险函数值/>
Figure GDA00042662328700000318
风险函数的计算公式为：
Figure GDA00042662328700000319
对于公式(3)，
Figure GDA00042662328700000320
表示t时刻j号子集合中第i个粒子/>
Figure GDA00042662328700000321
的风险函数值；0＜λ＜1，i＝1,2,…,N，q≥1；f(·)表示系统状态的状态转移函数；
步骤2.3、由子集合粒子的代价函数值
Figure GDA00042662328700000322
或风险函数值/>
Figure GDA00042662328700000323
计算粒子权值/>
Figure GDA00042662328700000324
计算方法如下：
Figure GDA00042662328700000325
亦可通过如下方法计算：
Figure GDA0004266232870000041
在公式(4)、(5)中，
Figure GDA0004266232870000042
即为所求的t时刻j号子集合中第i个粒子/>
Figure GDA0004266232870000043
的粒子权值，粒子权值调节参数β＞1；公式(4)中/>
Figure GDA0004266232870000044
为粒子/>
Figure GDA0004266232870000045
的代价值函数值；公式(5)中/>
Figure GDA0004266232870000046
表示粒子/>
Figure GDA0004266232870000047
的风险函数值，/>
Figure GDA0004266232870000048
表示t时刻j号子集合N个随机估计样本中风险函数的最小值，δ为调节参数，0＜δ＜1。
步骤3具体如下：
步骤3.1、在区间[0,1]上，随机生成服从均匀分布的随机数集合{ui}i＝1:N，且集合中的各元素满足独立同分布；
步骤3.2、权值的累积值Ii＝cdf{ui}，法则cdf表示粒子权重累积分布函数，同时对于随机数ui，存在
Figure GDA0004266232870000049
步骤3.3、令粒子权值wk＝1/N，在进行多项式重采样操作后，得到复制后粒子的数目集合{vi}i＝1:N，0≤vi≤m，vi为第i个粒子经过重采样后被复制的次数；
步骤3.4、由公式(4)知，粒子的代价函数值越小，权值越大，因此重采样步骤中权值大的粒子被复制的次数越多，同时保存大权值粒子相应的代价函数值；经过多项式重采样操作后，最终得到t时刻粒子—代价函数值的集合
Figure GDA00042662328700000410
j＝1,2。
步骤4具体如下：
步骤4.1、1号子集合
Figure GDA00042662328700000411
的粒子以及粒子对应的代价函数值与权值的更新如下：
当t＝1时，1号子集合的粒子更新为
Figure GDA0004266232870000051
当t≥2时，粒子基于高斯分布更新，即/>
Figure GDA0004266232870000052
表示经历多项式重采样后t-1时刻粒子的状态值，/>
Figure GDA0004266232870000053
表示协方差矩阵，恒等函数/>
Figure GDA0004266232870000054
与方差σt的计算方式如下：
Figure GDA0004266232870000055
Figure GDA0004266232870000056
公式(6)、(7)中，cov(·)表示协方差的运算。
代价函数值的更新由公式(1)、(2)实现；权值按照公式(4)或公式(5)进行更新；更新后，得到t时刻新的1号粒子—代价集合
Figure GDA0004266232870000057
步骤4.2、2号子集合
Figure GDA0004266232870000058
粒子以及粒子对应的代价函数值与权值的更新如下：
当t＝1时，2号子集合的粒子更新为
Figure GDA0004266232870000059
当t≥2时，2号子集合中粒子的更新是基于柯西分布实现的，即/>
Figure GDA00042662328700000510
γ(γ＞0)表示尺度参数，l0表示位置参数，γ与l0的取值与系统有关；代价函数值和权值的更新方法与1号子集合相同，最终，得到t时刻新的2号粒子—代价集合/>
Figure GDA00042662328700000511
步骤5中假设两个子集合到达m时刻，m∈{1,2,…,T}，粒子进行信息交互操作，具体操作如下：
步骤5.1、取1号子集合的粒子—代价集合
Figure GDA00042662328700000512
以及排序后2号子集合/>
Figure GDA00042662328700000513
前M个粒子及粒子对应代价值的集合/>
Figure GDA0004266232870000061
M∈{1,2,…,N}，且集合/>
Figure GDA0004266232870000062
中粒子/>
Figure GDA0004266232870000063
对应代价值/>
Figure GDA0004266232870000064
满足条件/>
Figure GDA0004266232870000065
将上述两个集合合并为新的粒子—代价集合/>
Figure GDA0004266232870000066
步骤5.2、将集合
Figure GDA0004266232870000067
粒子的代价值按升序排序，得到集合/>
Figure GDA0004266232870000068
即
Figure GDA0004266232870000069
同时集合/>
Figure GDA00042662328700000610
中粒子/>
Figure GDA00042662328700000611
的代价值/>
Figure GDA00042662328700000612
必须满足条件
Figure GDA00042662328700000613
然后，取集合/>
Figure GDA00042662328700000614
前N个粒子及粒子相应的代价值组成集合/>
Figure GDA00042662328700000615
作为1号子集合m时刻的粒子集合；
步骤5.3、对2号子集合
Figure GDA00042662328700000616
进行步骤5.1和步骤5.2的类似操作，得到第m代2号子集合的粒子—代价集合/>
Figure GDA00042662328700000617
且集合/>
Figure GDA00042662328700000618
中每个粒子/>
Figure GDA00042662328700000619
的代价值/>
Figure GDA00042662328700000620
满足/>
Figure GDA00042662328700000621
经过步骤5最终得到1号子集合的粒子—代价集合
Figure GDA00042662328700000622
2号子集合的粒子—代价集合/>
Figure GDA00042662328700000623
步骤6具体如下：
在粒子—代价集合Ej(j＝1,2)中，分别取t时刻第i个粒子及粒子相应的代价函数值
Figure GDA00042662328700000624
若1号子集合E1中粒子/>
Figure GDA00042662328700000625
的代价函数值/>
Figure GDA00042662328700000626
大于2号子集合E2中粒子/>
Figure GDA00042662328700000627
的代价函数值/>
Figure GDA00042662328700000628
即/>
Figure GDA00042662328700000629
则选取/>
Figure GDA00042662328700000630
进入最终的粒子—代价集合/>
Figure GDA00042662328700000631
否则选取/>
Figure GDA00042662328700000632
进入集合/>
Figure GDA00042662328700000633
经过上述操作后，得到进行系统滤波所需的粒子—代价集合/>
Figure GDA00042662328700000634
步骤7具体如下：
参照步骤4中1号子集合的更新方法，对粒子—代价集合
Figure GDA00042662328700000635
进行更新操作，得到更新后的粒子—代价集合/>
Figure GDA00042662328700000636
步骤8具体如下：
按照步骤2.3中的方法计算由步骤7所得粒子—代价集合E中粒子对应的权值，并进行加权平均处理，最终得到t时刻系统状态的最优估计如下
Figure GDA0004266232870000071
公式(8)中，
Figure GDA0004266232870000072
即为系统t时刻系统状态的最优估计值，即为整个滤波过程需要求得的最终结果，/>
Figure GDA0004266232870000073
为t时刻系统状态的随机估计值，/>
Figure GDA0004266232870000074
为随机估计值/>
Figure GDA0004266232870000075
对应的权值，t＝1,2,…,T。
本发明的有益效果是，一种基于混合概率选择算子的代价参考粒子滤波方法，通过基于混合概率选择算子的代价参考粒子滤波方法对噪声统计特性未知的系统进行处理。对两个规模相同的子集合进行标准CRPF操作，但分别基于高斯分布和柯西分布更新子集合的粒子集合，并且让两个子集合在规定的步数进行子集合之间的信息交互操作。在此之后，将两个子基集合中粒子相应的代价函数值进行比较，选择代价值小的粒子进入用于最终进行系统滤波的粒子—代价集合，以此取代标准CRPF方法的重采样环节。然后，对经过选择得到的粒子—代价集合进行更新操作，最后对系统进行滤波处理。本发明通过子集合的信息交互以及粒子选择操作，改善粒子多样性，提高滤波的精度与准确度，使得系统滤波的效果更好。
附图说明
图1是本发明的信息交互操作示意图；
图2是本发明的粒子选择操作示意图；
图3是本发明实施例中电池容量衰退数据分布图；
图4是本发明与标准CRPF的滤波结果比较图；
图5是本发明与标准CRPF的绝对偏差对比图；
图6是本发明与标准CRPF的MAPE对比图。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
基于混合概率选择算子的代价参考粒子滤波方法在处理噪声统计特性未知的系统时，在混合概率选择算子中，服从高斯分布的选择算子可提高局部搜索的能力，服从柯西分布的算子可提高全局搜索的能力，基于上述两种分布的混合概率选择算子可提高收敛速度和滤波精度。首先，按照标准CRPF的操作初始化一个规模为N用于最终系统滤波的粒子—代价集合(粒子即随机产生的系统状态值)与两个规模为N的子集合，两个子集合分别为1号子集合、2号子集合，其中，1号子集合在粒子集合的更新阶段基于高斯分布进行更新，2号子集合基于柯西分布进行类似操作。然后，分别输出两个子集合经过更新操作得到的粒子—代价子集合，除此之外子集合内部的操作与标准CRPF相同。当到达指定步数时，对两个子集合进行信息交互操作。以1号子集合为例，将1号子集合的全体与代价值(即代价函数值)按升序排序后的2号子集合对应的前M个粒子合成一个新的粒子—代价集合。再将该新集合的代价值按升序排序，取该新集合的代价值大小为前N个对应的粒子组成集合作为1号子集合的新一代，然后对2号子集合进行类似操作。接下来，对两个子集合进行粒子选择操作，即比较每一步两个子集合中粒子对应的代价函数值，为实现最终系统滤波的粒子—代价集合选取最优粒子及最优粒子对应的代价值。具体选取标准为：在每一次粒子选择的过程中，若1号子集合中粒子的代价值大于2号子集合的代价值，则将2号子集合中的粒子及粒子相应的代价值选入最终用于系统滤波的粒子—代价集合；否则选取1号子集合中的粒子及粒子相应的代价值。所述最终用于滤波的粒子—代价集合是通过上述信息交互与粒子选择的过程生成的，即按照代价函数值大小筛选分别基于高斯分布与柯西分布更新所生成两个子集合中的粒子，构成的粒子集合中既含有服从高斯分布的粒子又含有服从柯西分布的粒子，实现混合概率选择算子的运用。完成粒子选择的操作后，对所得粒子—代价集合进行更新和权值的加权平均处理，最终得到每一时刻系统状态的最优估计值，完成系统的滤波过程。
本发明一种基于混合概率选择算子的代价参考粒子滤波方法，具体按照以下步骤实施：
步骤1、如图1所示，随机产生N个实数作为系统状态的样本，所述样本称为粒子，并定义每一粒子对应的代价值均为0，由所述样本与代价值组成的集合称为粒子—代价集合；同时，分别随机产生粒子组成1号子集合与2号子集合，每一个子集合的数量为N；
步骤1生成初始规模为N的粒子—代价集合E0，表示为
Figure GDA0004266232870000091
其中，/>
Figure GDA0004266232870000092
表示系统状态随机估计值，/>
Figure GDA0004266232870000093
(I0的范围与系统有关)，即/>
Figure GDA0004266232870000094
服从均匀分布U(I0)，初始代价函数值/>
Figure GDA0004266232870000095
i＝1,2,…,N，N为系统状态随机估计值的采样样本数；
同时，生成两个大小为N的子集合
Figure GDA0004266232870000096
其中，/>
Figure GDA0004266232870000097
服从均匀分布U(I0)，j为子集合个数，j＝1,2，当j＝1时，/>
Figure GDA0004266232870000098
表示1号子集合，当j＝2时，/>
Figure GDA0004266232870000099
表示2号子集合，其余参数含义均相同。
步骤2、分别计算所有子集合中粒子的代价函数值、风险函数值及权值；
步骤2具体如下：
步骤2.1、计算t时刻子集合中粒子
Figure GDA00042662328700000910
的代价函数值/>
Figure GDA00042662328700000911
j＝1,2，i＝1,2,…,N，t＝1,2,…,T，T为时间序列长度，代价函数的计算方式如下：
Figure GDA0004266232870000101
Figure GDA0004266232870000102
在公式(1)、(2)中，
Figure GDA0004266232870000103
表示t时刻j号子集合中第i个粒子/>
Figure GDA0004266232870000104
的代价函数值；λ表示遗忘因子，0＜λ＜1；/>
Figure GDA0004266232870000105
表示t-1时刻j号子集合中第i个粒子/>
Figure GDA0004266232870000106
的代价值函数值；/>
Figure GDA0004266232870000107
表示粒子/>
Figure GDA0004266232870000108
的代价增量，/>
Figure GDA0004266232870000109
表示粒子/>
Figure GDA00042662328700001010
的代价增量函数；yt表示t时刻传感器的测量值，h(·)表示传感器估计值的观测函数，预设参数q满足q≥1；
步骤2.2、由子集合粒子的代价函数值
Figure GDA00042662328700001011
计算得到风险函数值/>
Figure GDA00042662328700001012
风险函数的计算公式为：
Figure GDA00042662328700001013
对于公式(3)，
Figure GDA00042662328700001014
表示t时刻j号子集合中第i个粒子/>
Figure GDA00042662328700001015
的风险函数值；0＜λ＜1，i＝1,2,…,N，q≥1；f(·)表示系统状态的状态转移函数；
步骤2.3、由子集合粒子的代价函数值
Figure GDA00042662328700001016
或风险函数值/>
Figure GDA00042662328700001017
计算粒子权值/>
Figure GDA00042662328700001018
计算方法如下：
Figure GDA00042662328700001019
亦可通过如下方法计算：
Figure GDA00042662328700001020
在公式(4)、(5)中，
Figure GDA00042662328700001021
即为所求的t时刻j号子集合中第i个粒子/>
Figure GDA00042662328700001022
的粒子权值，粒子权值调节参数β＞1；公式(4)中/>
Figure GDA00042662328700001023
为粒子/>
Figure GDA00042662328700001024
的代价值函数值；公式(5)中/>
Figure GDA00042662328700001025
表示粒子/>
Figure GDA00042662328700001026
的风险函数值，/>
Figure GDA00042662328700001027
表示t时刻j号子集合N个随机估计样本中风险函数的最小值，δ为调节参数，0＜δ＜1。
步骤3、重采样子集合，重采样方法选取多项式重采样；
步骤3具体如下：
步骤3.1、在区间[0,1]上，随机生成服从均匀分布的随机数集合{ui}i＝1:N，且集合中的各元素满足独立同分布；
步骤3.2、权值的累积值Ii＝cdf{ui}，法则cdf表示粒子权重累积分布函数，同时对于随机数ui，存在
Figure GDA0004266232870000111
步骤3.3、令粒子权值wk＝1/N，在进行多项式重采样操作后，得到复制后粒子的数目集合{vi}i＝1:N，0≤vi≤m，vi为第i个粒子经过重采样后被复制的次数；
步骤3.4、由公式(4)知，粒子的代价函数值越小，权值越大，因此重采样步骤中权值大的粒子被复制的次数越多，同时保存大权值粒子相应的代价函数值；经过多项式重采样操作后，最终得到t时刻粒子—代价函数值的集合
Figure GDA0004266232870000112
j＝1,2。
步骤4、更新子集合；
步骤4具体如下：
步骤4.1、1号子集合
Figure GDA0004266232870000113
的粒子以及粒子对应的代价函数值与权值的更新如下：
当t＝1时，1号子集合的粒子更新为
Figure GDA0004266232870000114
当t≥2时，粒子基于高斯分布更新，即/>
Figure GDA0004266232870000115
表示经历多项式重采样后t-1时刻粒子的状态值，/>
Figure GDA0004266232870000116
表示协方差矩阵，恒等函数/>
Figure GDA0004266232870000117
与方差σt的计算方式如下：
Figure GDA0004266232870000121
Figure GDA0004266232870000122
公式(6)、(7)中，cov(·)表示协方差的运算。
代价函数值的更新由公式(1)、(2)实现；权值按照公式(4)或公式(5)进行更新；更新后，得到t时刻新的1号粒子—代价集合
Figure GDA0004266232870000123
步骤4.2、2号子集合
Figure GDA0004266232870000124
粒子以及粒子对应的代价函数值与权值的更新如下：
当t＝1时，2号子集合的粒子更新为
Figure GDA0004266232870000125
当t≥2时，2号子集合中粒子的更新是基于柯西分布实现的，即/>
Figure GDA0004266232870000126
γ(γ＞0)表示尺度参数，l0表示位置参数，γ与l0的取值与系统有关；代价函数值和权值的更新方法与1号子集合相同，最终，得到t时刻新的2号粒子—代价集合/>
Figure GDA0004266232870000127
步骤5、子集合之间粒子的信息交互，若到达设定的交换步数，则执行此项操作；否则，跳过该步骤，执行步骤6；
步骤5中假设两个子集合到达m时刻，m∈{1,2,…,T}，粒子进行信息交互操作，由图1可知具体操作如下：
步骤5.1、取1号子集合的粒子—代价集合
Figure GDA0004266232870000128
以及排序后2号子集合/>
Figure GDA0004266232870000129
前M个粒子及粒子对应代价值的集合/>
Figure GDA00042662328700001210
M∈{1,2,…,N}，且集合/>
Figure GDA00042662328700001211
中粒子/>
Figure GDA00042662328700001212
对应代价值/>
Figure GDA00042662328700001213
满足条件/>
Figure GDA00042662328700001214
将上述两个集合合并为新的粒子—代价集合/>
Figure GDA0004266232870000131
步骤5.2、将集合
Figure GDA0004266232870000132
粒子的代价值按升序排序，得到集合/>
Figure GDA0004266232870000133
即
Figure GDA0004266232870000134
同时集合/>
Figure GDA0004266232870000135
中粒子/>
Figure GDA0004266232870000136
的代价值/>
Figure GDA0004266232870000137
必须满足条件
Figure GDA0004266232870000138
然后，取集合/>
Figure GDA0004266232870000139
前N个粒子及粒子相应的代价值组成集合/>
Figure GDA00042662328700001310
作为1号子集合m时刻的粒子集合；
步骤5.3、对2号子集合
Figure GDA00042662328700001311
进行步骤5.1和步骤5.2的类似操作，得到第m代2号子集合的粒子—代价集合/>
Figure GDA00042662328700001312
且集合/>
Figure GDA00042662328700001313
中每个粒子/>
Figure GDA00042662328700001314
的代价值/>
Figure GDA00042662328700001315
满足/>
Figure GDA00042662328700001316
经过步骤5最终得到1号子集合的粒子—代价集合
Figure GDA00042662328700001317
2号子集合的粒子—代价集合/>
Figure GDA00042662328700001318
步骤6、生成最终用于系统滤波的粒子—代价集合；
步骤6具体如下：
如图2所示，在粒子—代价集合Ej(j＝1,2)中，分别取t时刻第i个粒子及粒子相应的代价函数值
Figure GDA00042662328700001319
若1号子集合E1中粒子/>
Figure GDA00042662328700001320
的代价函数值/>
Figure GDA00042662328700001321
大于2号子集合E2中粒子/>
Figure GDA00042662328700001322
的代价函数值/>
Figure GDA00042662328700001323
即/>
Figure GDA00042662328700001324
则选取/>
Figure GDA00042662328700001325
进入最终的粒子—代价集合
Figure GDA00042662328700001326
否则选取/>
Figure GDA00042662328700001327
进入集合/>
Figure GDA00042662328700001328
经过上述操作后，得到进行系统滤波所需的粒子—代价集合/>
Figure GDA00042662328700001329
步骤7、更新最终的粒子—代价集合；
步骤7具体如下：
参照步骤4中1号子集合的更新方法，对粒子—代价集合
Figure GDA00042662328700001330
进行更新操作，得到更新后的粒子—代价集合/>
Figure GDA0004266232870000141
步骤8、系统滤波，最终得到系统状态的最优估计值。
步骤8具体如下：
按照步骤2.3中的方法计算由步骤7所得粒子—代价集合E中粒子对应的权值，并进行加权平均处理，最终得到t时刻系统状态的最优估计如下
Figure GDA0004266232870000142
公式(8)中，
Figure GDA0004266232870000143
即为系统t时刻系统状态的最优估计值，即为整个滤波过程需要求得的最终结果，/>
Figure GDA0004266232870000144
为t时刻系统状态的随机估计值，/>
Figure GDA0004266232870000145
为随机估计值/>
Figure GDA0004266232870000146
对应的权值，t＝1,2,…,T。/>
Figure GDA0004266232870000147
作为系统的滤波结果可以为系统后续的分析与控制提供必要的数据支撑与基本保障。
本发明公开了一种基于混合概率选择算子的代价参考粒子滤波方法，对噪声统计特性未知的系统进行处理，采用分别基于高斯分布与柯西分布更新的子集合，对子集合进行信息交互与粒子选择操作，增加粒子的多样性，得到进行系统滤波所需的粒子—代价集合后更新该集合，最后，计算更新后的粒子集合中粒子对应的权值，并进行加权平均处理，最终得到系统状态的最优估计值，即为最终的系统滤波将结果，实现对传感器所采集的观测数据的滤波操作。本发明无需噪声先验分布函数即可在过程噪声和测量噪声未知的情况下进行滤波，提高了滤波结果的精度与准确性，减小了系统滤波偏差，改善了滤波效果，有利于对系统进行进一步的分析与控制。本发明解决了现有技术中对噪声统计特性未知系统滤波准确度较低的问题。
实施例
本实施例对型号Li-ion 18650、额定容量2Ah的锂电池进行电池容量预测，电池容量衰退数据如图3所示。
为更好地预测电池的容量衰退过程，需确定系统状态空间模型，状态方程如公式(9)所示，观测方程如公式(10)所示。模型参数的初始值设置为a1＝1.942，b1＝-2.052×10-3，c1＝1.57×10-7，d1＝0.07406。
Figure GDA0004266232870000151
Figure GDA0004266232870000152
在公式(9)中，t＝2,3,...,168，at、ct为电池内阻抗相关参数在t时刻的估计值，bt、dt为电池老化速率相关参数在t时刻的估计值，wa、wb、wc、wd分别代表a、b、c、d的过程噪声；公式(10)中Qt为t时刻电池容量的估计值，测量噪声vt～N(0,1)，N(0,1)即标准正态分布，即均值为0，方差为1的高斯分布。
设置仿真步数T为168，将标准CRPF与基于混合概率选择算子的代价参考粒子滤波方法(本发明)粒子集合的粒子数均设置为N＝150，且重采样步骤均使用多项式重采样的方法。设本发明中1号子集合与2号子集合每G＝10代进行子集合间的信息交互，且在交流互换时取对方子集合的优秀个体数量为M＝70。特别地，本发明的2号子集合在更新粒子集合时，基于公式(11)服从柯西分布进行更新。
Figure GDA0004266232870000153
在公式(11)中，
Figure GDA0004266232870000161
分别表示参数a、b、c、d经过更新得到的2号子集合t时刻第i个粒子。
设置本发明步骤中涉及到的参数q＝2，遗忘因子λ＝1×10-6，参数β＝1.5，调节参数δ＝0.001。
完成图1所示的信息交互操作，以及图2所示的粒子选择操作；对图3中B0007电池容量衰退数据进行处理；得到系统的滤波结果，如图4所示。
为了更加准确分析本发明的有效性，选取绝对偏差、RMSE(即均方根误差)以及MAPE(即平均绝对百分比误差)作为系统指标进行比较分析。图5为标准CRPF与本发明的偏差对比图；图6为标准CRPF与本发明每一周期的MAPE对比图；由表1可知：改进后，本发明完整滤波过程的RMSE减少至标准CRPF方法RMSE的51.482％，本发明整体的MAPE与标准CRPF方法的MAPE相比减少了0.50775％.
表1
Figure GDA0004266232870000162
由上述实验结果可知，与标准CRPF相比，本发明(即基于混合概率选择算子的代价参考粒子滤波方法)的偏差更小，滤波效果更好，滤波的准确性与精度更高。

Claims (4)

1.一种基于混合概率选择算子的代价参考粒子滤波方法，其特征在于，具体按照以下步骤实施：
步骤1、随机产生N个实数作为系统状态的样本，所述样本称为粒子，并定义每一粒子对应的代价值均为0，由所述样本与代价值组成的集合称为粒子—代价集合；同时，分别随机产生粒子组成1号子集合与2号子集合，每一个子集合的数量为N；
所述步骤1生成初始规模为N的粒子—代价集合E0，表示为
Figure FDA0004266232850000011
其中，/>
Figure FDA0004266232850000012
表示系统状态随机估计值，/>
Figure FDA0004266232850000013
即/>
Figure FDA0004266232850000014
服从均匀分布U(I0)，初始代价函数值
Figure FDA0004266232850000015
N为系统状态随机估计值的采样样本数；
同时，生成两个大小为N的子集合
Figure FDA0004266232850000016
其中，/>
Figure FDA0004266232850000017
服从均匀分布U(I0)，j为子集合个数，j＝1,2，当j＝1时，/>
Figure FDA0004266232850000018
表示1号子集合，当j＝2时，/>
Figure FDA0004266232850000019
表示2号子集合；
步骤2、分别计算所有子集合中粒子的代价函数值、风险函数值及权值；
所述步骤2具体如下：
步骤2.1、计算t时刻子集合中粒子
Figure FDA00042662328500000110
的代价函数值/>
Figure FDA00042662328500000111
Figure FDA00042662328500000112
T为时间序列长度，代价函数的计算方式如下：
Figure FDA00042662328500000113
Figure FDA00042662328500000114
在公式(1)、(2)中，
Figure FDA00042662328500000115
表示t时刻j号子集合中第i个粒子/>
Figure FDA00042662328500000116
的代价函数值；λ表示遗忘因子，0＜λ＜1；/>
Figure FDA00042662328500000117
表示t-1时刻j号子集合中第i个粒子/>
Figure FDA0004266232850000021
的代价值函数值；/>
Figure FDA0004266232850000022
表示粒子/>
Figure FDA0004266232850000023
的代价增量，/>
Figure FDA0004266232850000024
表示粒子/>
Figure FDA0004266232850000025
的代价增量函数；yt表示t时刻传感器的测量值，h(·)表示传感器估计值的观测函数，预设参数q满足q≥1；
步骤2.2、由子集合粒子的代价函数值
Figure FDA00042662328500000222
计算得到风险函数值/>
Figure FDA0004266232850000026
风险函数的计算公式为：
Figure FDA0004266232850000027
对于公式(3)，
Figure FDA0004266232850000028
表示t时刻j号子集合中第i个粒子/>
Figure FDA0004266232850000029
的风险函数值；0＜λ＜1，i＝1,2,…,N，q≥1；f(·)表示系统状态的状态转移函数；
步骤2.3、由子集合粒子的代价函数值
Figure FDA00042662328500000210
或风险函数值/>
Figure FDA00042662328500000211
计算粒子权值/>
Figure FDA00042662328500000212
计算方法如下：
Figure FDA00042662328500000213
亦可通过如下方法计算：
Figure FDA00042662328500000214
在公式(4)、(5)中，
Figure FDA00042662328500000215
即为所求的t时刻j号子集合中第i个粒子/>
Figure FDA00042662328500000216
的粒子权值，粒子权值调节参数β＞1；公式(4)中/>
Figure FDA00042662328500000217
为粒子/>
Figure FDA00042662328500000218
的代价值函数值；公式(5)中/>
Figure FDA00042662328500000219
表示粒子/>
Figure FDA00042662328500000220
的风险函数值，/>
Figure FDA00042662328500000221
表示t时刻j号子集合N个随机估计样本中风险函数的最小值，δ为调节参数，0＜δ＜1；
步骤3、重采样子集合，重采样方法选取多项式重采样；
步骤3具体如下：
步骤3.1、在区间[0,1]上，随机生成服从均匀分布的随机数集合{ui}i＝1:N，且集合中的各元素满足独立同分布；
步骤3.2、权值的累积值Ii＝cdf{ui}，法则cdf表示粒子权重累积分布函数，同时对于随机数ui，存在
Figure FDA0004266232850000031
步骤3.3、令粒子权值wk＝1/N，在进行多项式重采样操作后，得到复制后粒子的数目集合{vi}i＝1:N，0≤vi≤m，vi为第i个粒子经过重采样后被复制的次数；
步骤3.4、由公式(4)知，粒子的代价函数值越小，权值越大，因此重采样步骤中权值大的粒子被复制的次数越多，同时保存大权值粒子相应的代价函数值；经过多项式重采样操作后，最终得到t时刻粒子—代价函数值的集合
Figure FDA0004266232850000032
步骤4、更新子集合；
步骤4具体如下：
步骤4.1、1号子集合
Figure FDA0004266232850000033
的粒子以及粒子对应的代价函数值与权值的更新如下：
当t＝1时，1号子集合的粒子更新为
Figure FDA0004266232850000034
当t≥2时，粒子基于高斯分布更新，即/>
Figure FDA0004266232850000035
Figure FDA0004266232850000036
表示经历多项式重采样后t-1时刻粒子的状态值，/>
Figure FDA0004266232850000037
表示协方差矩阵，恒等函数/>
Figure FDA0004266232850000038
与方差σt的计算方式如下：
Figure FDA0004266232850000039
Figure FDA00042662328500000310
公式(6)、(7)中，cov(·)表示协方差的运算；
代价函数值的更新由公式(1)、(2)实现；权值按照公式(4)或公式(5)进行更新；更新后，得到t时刻新的1号粒子—代价集合
Figure FDA0004266232850000041
步骤4.2、2号子集合
Figure FDA0004266232850000042
粒子以及粒子对应的代价函数值与权值的更新如下：
当t＝1时，2号子集合的粒子更新为
Figure FDA0004266232850000043
当t≥2时，2号子集合中粒子的更新是基于柯西分布实现的，即/>
Figure FDA0004266232850000044
γ(γ＞0)表示尺度参数，l0表示位置参数，γ与l0的取值与系统有关；代价函数值和权值的更新方法与1号子集合相同，最终，得到t时刻新的2号粒子—代价集合/>
Figure FDA0004266232850000045
步骤5、子集合之间粒子的信息交互，若到达设定的交换步数，则执行此项操作；否则，无需执行信息交互操作；
步骤5中两个子集合每G代进行一次子集合之间的信息交互操作，若t≠m时，则无需对子集合进行操作，m为G的整数倍，m＜T，令
Figure FDA0004266232850000046
若t＝m，则进行信息交互操作，具体操作如下：
步骤5.1、取1号子集合的粒子—代价集合
Figure FDA0004266232850000047
以及排序后2号子集合/>
Figure FDA0004266232850000048
前M个粒子及粒子对应代价值的集合/>
Figure FDA0004266232850000049
且集合/>
Figure FDA00042662328500000410
中粒子/>
Figure FDA00042662328500000411
对应代价值/>
Figure FDA00042662328500000412
满足条件/>
Figure FDA00042662328500000413
将上述两个集合合并为新的粒子—代价集合/>
Figure FDA00042662328500000414
步骤5.2、将集合
Figure FDA00042662328500000415
粒子的代价值按升序排序，得到集合/>
Figure FDA00042662328500000416
即/>
Figure FDA0004266232850000051
同时集合/>
Figure FDA0004266232850000052
中粒子/>
Figure FDA0004266232850000053
的代价值/>
Figure FDA0004266232850000054
必须满足条件
Figure FDA0004266232850000055
然后，取集合/>
Figure FDA0004266232850000056
前N个粒子及粒子相应的代价值组成集合/>
Figure FDA0004266232850000057
作为1号子集合m时刻的粒子集合；
步骤5.3、对2号子集合
Figure FDA0004266232850000058
进行步骤5.1和步骤5.2的类似操作，得到第m代2号子集合的粒子—代价集合/>
Figure FDA0004266232850000059
且集合/>
Figure FDA00042662328500000510
中每个粒子/>
Figure FDA00042662328500000511
的代价值/>
Figure FDA00042662328500000512
满足/>
Figure FDA00042662328500000513
经过步骤5最终得到1号子集合的粒子—代价集合
Figure FDA00042662328500000514
2号子集合的粒子—代价集合/>
Figure FDA00042662328500000515
步骤6、生成最终用于系统滤波的粒子—代价集合；
步骤7、更新最终的粒子—代价集合；
步骤8、系统滤波，最终得到系统状态的最优估计值。
2.根据权利要求1所述的一种基于混合概率选择算子的代价参考粒子滤波方法，其特征在于，所述步骤6具体如下：
在粒子—代价集合Ej(j＝1,2)中，分别取t时刻第i个粒子及粒子相应的代价函数值
Figure FDA00042662328500000516
若1号子集合E1中粒子/>
Figure FDA00042662328500000517
的代价函数值/>
Figure FDA00042662328500000518
大于2号子集合E2中粒子/>
Figure FDA00042662328500000519
的代价函数值/>
Figure FDA00042662328500000520
即/>
Figure FDA00042662328500000521
则选取/>
Figure FDA00042662328500000522
进入最终的粒子—代价集合/>
Figure FDA00042662328500000523
否则选取/>
Figure FDA00042662328500000524
进入集合/>
Figure FDA00042662328500000525
经过上述操作后，得到进行系统滤波所需的粒子—代价集合/>
Figure FDA00042662328500000526
3.根据权利要求2所述的一种基于混合概率选择算子的代价参考粒子滤波方法，其特征在于，所述步骤7具体如下：
参照步骤4中1号子集合的更新方法，对粒子—代价集合
Figure FDA00042662328500000527
进行更新操作，得到更新后的粒子—代价集合/>
Figure FDA0004266232850000061
4.根据权利要求3所述的一种基于混合概率选择算子的代价参考粒子滤波方法，其特征在于，所述步骤8具体如下：
按照步骤2.3中的方法计算由步骤7所得粒子—代价集合E中粒子对应的权值，并进行加权平均处理，最终得到t时刻系统状态的最优估计如下
Figure FDA0004266232850000062
公式(8)中，
Figure FDA0004266232850000063
即为系统t时刻系统状态的最优估计值，即为整个滤波过程需要求得的最终结果，/>
Figure FDA0004266232850000064
为t时刻系统状态的随机估计值，/>
Figure FDA0004266232850000065
为随机估计值/>
Figure FDA0004266232850000066
对应的权值，t＝1,2,…,T。
CN202111017989.3A 2021-08-30 2021-08-30 一种基于混合概率选择算子的代价参考粒子滤波方法 Active CN113810024B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111017989.3A CN113810024B (zh) 2021-08-30 2021-08-30 一种基于混合概率选择算子的代价参考粒子滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111017989.3A CN113810024B (zh) 2021-08-30 2021-08-30 一种基于混合概率选择算子的代价参考粒子滤波方法

Publications (2)

Publication Number Publication Date
CN113810024A CN113810024A (zh) 2021-12-17
CN113810024B true CN113810024B (zh) 2023-07-14

Family

ID=78894524

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111017989.3A Active CN113810024B (zh) 2021-08-30 2021-08-30 一种基于混合概率选择算子的代价参考粒子滤波方法

Country Status (1)

Country Link
CN (1) CN113810024B (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO1999046731A1 (en) * 1998-03-13 1999-09-16 The University Of Houston System Methods for performing daf data filtering and padding
GB201010718D0 (en) * 2010-06-25 2010-08-11 Total Sa An improved process for characterising the evolution of an oil or gas reservoir over time
GB201619902D0 (en) * 2016-11-24 2017-01-11 Univ Oxford Innovation Ltd Patient status monitor and method of monitoring patient status
CN106501800A (zh) * 2016-10-28 2017-03-15 中国人民解放军信息工程大学 基于代价参考粒子滤波的mimo雷达目标检测前跟踪方法
CN110568494A (zh) * 2019-09-12 2019-12-13 电子科技大学 基于广义极值分布的叠前非高斯avo反演方法
CN113239628A (zh) * 2021-06-02 2021-08-10 哈尔滨工程大学 基于量子海鸥演化机制加权Myriad滤波器设计方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO1999046731A1 (en) * 1998-03-13 1999-09-16 The University Of Houston System Methods for performing daf data filtering and padding
GB201010718D0 (en) * 2010-06-25 2010-08-11 Total Sa An improved process for characterising the evolution of an oil or gas reservoir over time
CN106501800A (zh) * 2016-10-28 2017-03-15 中国人民解放军信息工程大学 基于代价参考粒子滤波的mimo雷达目标检测前跟踪方法
GB201619902D0 (en) * 2016-11-24 2017-01-11 Univ Oxford Innovation Ltd Patient status monitor and method of monitoring patient status
CN110568494A (zh) * 2019-09-12 2019-12-13 电子科技大学 基于广义极值分布的叠前非高斯avo反演方法
CN113239628A (zh) * 2021-06-02 2021-08-10 哈尔滨工程大学 基于量子海鸥演化机制加权Myriad滤波器设计方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Silicon Melt Liquid Level Detection Based on Parallel CRPF;Xinyu Zhang;《2018 Chinese Automation Congress (CAC)》;全文 *
基于CRPF的MIMO雷达目标检测前跟踪算法;秦文利;《四川大学学报自然科学版》;全文 *
智能优化的代价评估粒子滤波算法;王进花;曹洁;李伟;;系统工程与电子技术(12);全文 *
结合质心思想和柯西变异策略的粒子群优化算法;吕立国;季伟东;;计算机应用(05);全文 *

Also Published As

Publication number Publication date
CN113810024A (zh) 2021-12-17

Similar Documents

Publication Publication Date Title
CN110084367B (zh) 一种基于lstm深度学习模型的土壤墒情预测方法
Zhang et al. A novel approach of battery pack state of health estimation using artificial intelligence optimization algorithm
CN110135630B (zh) 基于随机森林回归和多步寻优的短期负荷需求预测方法
CN110457789B (zh) 一种锂离子电池剩余寿命预测方法
CN111563706A (zh) 一种基于lstm网络的多变量物流货运量预测方法
CN101893884B (zh) 密炼机的橡胶混炼过程中质量指标数据的软测量方法
CN111260136A (zh) 一种基于arima-lstm组合模型的楼宇短期负荷预测方法
CN111861013B (zh) 一种电力负荷预测方法及装置
CN107957555A (zh) 一种估算动力锂电池SoC的新方法
CN106600037B (zh) 一种基于主成分分析的多参量辅助负荷预测方法
CN109492748B (zh) 一种基于卷积神经网络的电力系统的中长期负荷预测模型建立方法
CN112288137A (zh) 一种计及电价和Attention机制的LSTM短期负荷预测方法及装置
CN102073922A (zh) 基于影响因素筛选的短期负荷预测方法
CN112946484A (zh) 一种基于bp神经网络的soc估计方法、系统、终端设备及可读存储介质
CN112381279B (zh) 一种基于vmd和bls组合模型的风电功率预测方法
CN111310348A (zh) 一种基于pso-lssvm的材料本构模型预测方法
CN113139605A (zh) 基于主成分分析和lstm神经网络的电力负荷预测方法
CN113406503A (zh) 基于深度神经网络的锂电池soh在线估算方法
CN112766537A (zh) 一种短期电负荷预测方法
CN110807508B (zh) 计及复杂气象影响的母线峰值负荷预测方法
CN116106761A (zh) 基于典型相关分析的锂离子电池电量实时估计方法
CN115308558A (zh) Cmos器件寿命预测方法、装置、电子设备及介质
CN115096357A (zh) 一种基于ceemdan-pca-lstm的室内环境质量预测方法
CN110826794A (zh) 基于pso优化svm的电厂耗煤基准值滚动预测方法和装置
CN116826745B (zh) 一种电力系统背景下的分层分区短期负荷预测方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant