粒子滤波采样

粒子滤波算法的完整建立在Gordon，Salmond和Smith提出的重采样技术上，并且一种新的采样算法（采样-重要性重采样）的发现和不断改良也对粒子滤波算法进行了丰富。

粒子滤波重采样

常用的重采样方法可以分为四类：最临近值重采样法，双线性重采样法，双立方重采样法，插值重采样法。
1）最邻近值重采样法:指的是比较目标图像与原图像的宽或者高，并且以此将原图像相对位置的像素点作为目标图像像素点的值。
2）双线性重采样法:指的是参考原图像对应像素周围四个点的值，并且由相对位置为每个点取权值来获得目标图像。
3）双立方重采样法:参考原像素点周围4*4个像素的值，并且根据这来获得目标图像。
4）插值重采样法:相比上述算法，这种方法参考更多原图像数据信息，可被用于求解对称矩阵方程组的求解以及特征值问题，重采样效果一般更好。
重采样算法是按概率进行复制和淘汰步骤，权重高的可能会被多次复制，这就保证了整体粒子数基本不变，然后进行归一化，将所有的粒子的权重都变为 1/n，继续进行下一步的预测更新步骤。虽然重采样的应用不会彻底根除粒子退化问题，但也会大大改善了粒子退化问题。
虽然重采样方法能够一定程度上减弱粒子退化问题，但是也必然会导致例子多样性的丧失，即N=1/(∑▒W_K^i )，N越小粒子退化问题也就越严重，就更加需要重采样，这样多次的重采样当然也会减慢粒子滤波的速度。

电池容量衰减模型

根据电池容量计算公式：Q=∫_tc^td▒〖I(t)dt〗,其中tc表示的是充电时间，td表示的是电池放电时间，I(t)表示的是电流，Q为锂离子电池的实时容量，而经过时间t则是记录了电池从循环开始到实验结束的时间。经过Matlab曲线拟合工具箱对四种电池容量数据进行指数拟合，能够得到比较准确的指数衰退模型。

根据电池衰减模型，我们选用双指数经验模型。

Q表示的是k时刻的电池容量，其中，k为循环次数，a,b,c,d均与锂离子电池本身有关，所以当a,b,c,d估算越准确，那么对于电池本身的模拟也就越真实，也就越能准确预测电池寿命。

算法流程

总结粒子滤波算法方法，有流程图：

使用Python3.8、Numpy1.23、matplotlib3.6软件环境对上述算法进行设计。运行Python程序可以得到CS2电池的容量预测图，在实验中以CS2_36为样本集，对样本的前600个数据进行10:1随机抽样后使用Matlab拟合工具箱对双指数模型进行拟合。

曲线拟合

有了拟合的参数和算法方程，那就开始滤波了：
注意这里python加载的数据先前就已经计算好保存下来的数据

这是我工程的文件结构。你需要在上一篇博客里下载好马里兰大学的数据，处理好数据后保存为npy格式的文件，然后在这里应用。

先上粒子滤波的效果图：

蓝色是观测值，黄色是滤波计算值，从红色线开始，绿色是没有滤波的自然观测计算值。

# 基于粒子滤波的锂离子电池RUL预测
# From: SWUST IPC14 Dai
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import sqrtm
# Python3.8、Numpy1.23、matplotlib3.6
# 除去数据异常的较大值和较小值点
def drop_outlier(array, count, bins):index = []range_ = np.arange(1, count, bins)for i in range_[:-1]:array_lim = array[i:i + bins]sigma = np.std(array_lim)mean = np.mean(array_lim)th_max, th_min = mean + sigma * 2, mean - sigma * 2idx = np.where((array_lim < th_max) & (array_lim > th_min))idx = idx[0] + iindex.extend(list(idx))return np.array(index)# 剩余容量双指数方程，状态方程
def hfun(X, k):Q = X[0] * np.exp(X[1] * k) + X[2] * np.exp(X[3] * k)return Q# 重采样步骤：
# 预测：抽取新粒子
# 更新：更新粒子权值
# 状态估计
# 多项式重采样# 重采样
def randomR(inIndex, q):outIndex = np.zeros(np.shape(inIndex))num = np.shape(q)u = np.random.rand(num[0], 1)u = np.sort(u, axis = 0)u = np.array(u)l = np.cumsum(q, axis=0)i = 0for j in np.arange(0, num[0]):while (i <= (num[0]-1)) and (u[i] <= l[j]):outIndex[i] = ji = i + 1return outIndex# capacity：阶段放电容量
# resistance：放电阶段电池的的内阻平均值
# CCCT：恒定电流充电时间
# CVCT: 恒定电压充电时间
if __name__ == "__main__":# ========================================#                加载数据# ========================================Battery_name = 'CS2_38'Battery = np.load('dataset/' + Battery_name + '.npy', allow_pickle=True)Battery = Battery.item()battery = Battery[Battery_name]names = ['capacity', 'resistance', 'CCCT', 'CVCT']# battery[1] = battery[1] /(battery[1] + 1.1)# 粒子滤波步骤# 初始化# 更新粒子状态# 权值计算和归一化# 重采样# 判断程序是否结束，迭代# ========================================#           预估数据为电池容量# ========================================N = len(battery['cycle'])pf_Number = 200                # 粒子数Prediction_Interval = 300      # 未来趋势值,预测区间的大小cita = 1e-4wa = 0.000001                  # 设定状态初始值wb = 0.01wc = 0.1wd = 0.0001Q = cita * np.diag([wa, wb, wc, wd])  # Q为过程噪声方差矩阵，diag()创建指定对角矩阵F = np.eye(4)   # F为驱动计算矩阵，eye()单位对角矩阵R = 0.001       # 观测噪声的协方差# ========= 状态方程赋初值 ==============a = -0.0000083499b = 0.055237c = 0.90097d = -0.00088543X0 = np.mat([a, b, c, d]).T  # 矩阵转置# ========= 滤波器状态初始化 ==============Xpf = np.zeros([4, N])Xpf[:, 0:1] = X0  # 对齐数组# ========= 粒子集初始化 ==============Xm = np.zeros([4, pf_Number, N])for i in np.arange(0, pf_Number - 1):sqr1 = np.array(sqrtm(Q))sqr2 = sqr1.dot(np.random.randn(4, 1))  # 矩阵乘法，直接用*是矩阵点乘Xm[:, i:i + 1, 0] = X0 + sqr2  # 对齐数组,需要将矩阵对齐后才能相加# ========= 从数据集读取观测量 =============capacity = battery[names[0]]Z = np.array(capacity)# ========= Zm为滤波器预测观测值，Zm与Xm对应 =============Zm = np.zeros([1, pf_Number, N])# ========= Zpf与Xpf对应 =============Zpf = np.zeros([1, N])         # 计算中得到的Zpf为滤波器更新得到的容量值# ========= 权值初始化 =============Weight = np.zeros([N, pf_Number])    # 计算中得到的W为更新的粒子权重# 粒子滤波算法for k in np.arange(1, N - 1):# 重要性采样for i in np.arange(0, pf_Number - 1):sqr1 = np.array(sqrtm(Q))       # 观测噪声sqr2 = sqr1.dot(np.random.randn(4, 1))  # 矩阵乘法，直接用*是矩阵点乘Xm[:, i:i + 1, k] = F.dot(Xm[:, i:i + 1, k - 1]) + sqr2# 权值重要性计算for i in np.arange(0, pf_Number - 1):Zm[0, i:i + 1, k] = hfun(Xm[:, i:i + 1, k], k)      # 观测预测Weight[k, i] = np.exp(-(Z[k] - Zm[0, i:i + 1, k:k + 1]) ** 2 / 2 / R) + 1e-99  # 重要性权值计算，乘方用 **Weight[k, :] = Weight[k, :] / sum(Weight[k, :])    # 权值归一化# 重采样# 这里的重采样以权值为传入值，返回值为采样后的索引outlndex = randomR(np.arange(0, pf_Number), Weight[k, :])# 得到新的样本for i in np.arange(0, len(outlndex)):Xm[:, i, k] = Xm[:, int(outlndex[i]), k]# 滤波后的状态更新，更新参数[a,b,c,d]Xpf[:, k] = [np.mean(Xm[0, :, k]),np.mean(Xm[1, :, k]),np.mean(Xm[2, :, k]),np.mean(Xm[3, :, k])]# 更新后的状态计算预测的容量值Zpf[0, k] = hfun(Xpf[:, k], k)# ========================================#         计算自然条件下的预测值# ========================================start = N - Prediction_Interval    # 预测的区间Zf = np.zeros(Prediction_Interval)  # 自然预测值Xf = np.zeros(Prediction_Interval)for k in np.arange(start-1, N-1):Zf[k-start] = hfun(Xpf[:, start], k)Xf[k-start] = k# 画线nax = [start, start]nay = [0, 1]plt.figure(figsize=(12, 9))plt.title('Particle filter  '+ Battery_name)  # 折线图标题plt.xlabel('Number of Cycles', fontsize=14)plt.ylim((0, 1.1))plt.ylabel(names[0], fontsize=14)plt.plot(battery['cycle'],  Z,          markersize=3)plt.plot(battery['cycle'],  Zpf[0, :],  markersize=3)plt.plot(Xf,                Zf,         markersize=3)plt.plot(nax,               nay,        linewidth=4)plt.legend(['Measured value', 'pf Predictive value', 'Natural predicted value'])plt.show()

Matlab的粒子滤波

main.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 版权声明：
%     本程序的详细中文注释请参考
%     黄小平，王岩，缪鹏程.粒子滤波原理及应用[M].电子工业出版社，2017.4
%     书中有原理介绍+例子+程序+中文注释
%     如果此程序有错误，请对提示修改
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  函数功能：粒子滤波用于电源寿命预测
% function mainload Battery_Capacity
%%load Battery_Capacity
N=length(A5Cycle);
% error('下面的参数M请参考书中的值设置，然后删除本行代码')
M=200;   %%粒子数
Future_Cycle=100;  % 未来趋势
if N>260N=260;   % 过滤大于260的数字
end%过程噪声协方差Q
cita=1e-4;
wa=0.000001;wb=0.01;wc=0.1;wd=0.0001;
Q=cita*diag([wa,wb,wc,wd]);%驱动矩阵
F=eye(4);%观测噪声协方差
R=0.001;a=-0.0000083499;b=0.055237;c=0.90097;d=-0.00088543;
X0=[a,b,c,d]';%滤波器状态初始化
Xpf=zeros(4,N);
Xpf(:,1)=X0;% 粒子集初始化
Xm=zeros(4,M,N);
for i=1:MXm(:,i,1)=X0+sqrtm(Q)*randn(4,1);
end% 观测量
Z(1,1:N)=A12Capacity(1:N,:)';Zm=zeros(1,M,N);Zpf=zeros(1,N);W=zeros(N,M);%粒子滤波算法
for k=2:N%  采样for i=1:MXm(:,i,k)=F*Xm(:,i,k-1)+sqrtm(Q)*randn(4,1);end% 权值重要性计算for i=1:MZm(1,i,k)=feval('hfun',Xm(:,i,k),k);W(k,i)=exp(-(Z(1,k)-Zm(1,i,k))^2/2/R)+1e-99;end,W(k,:)=W(k,:)./sum(W(k,:));% 重采样outIndex = randomR(1:M,W(k,:)');      % 调用外部函数% 得到新样本Xm( :,  :,  k)=Xm(  :,  outIndex,  k);% 滤波后的状态更新Xpf(:,k)=[mean(Xm(1,:,k));mean(Xm(2,:,k));mean(Xm(3,:,k));mean(Xm(4,:,k))];% 更新后的状态计算Zpf(1,k)=feval('hfun',Xpf(:,k),k);
end%预测未来电容的趋势
start=N-Future_Cycle;
for k=start:NZf(1,k-start+1)=feval('hfun',Xpf(:,start),k);Xf(1,k-start+1)=k;
endXreal=[a*ones(1,M);b*ones(1,M);c*ones(1,M);d*ones(1,M)];
figure
subplot(2,2,1);
hold on;box on;
plot(Xpf(1,:),'-r.');plot(Xreal(1,:),'-b.')
legend('粒子滤波后的a','平均值a')
subplot(2,2,2);
hold on;box on;
plot(Xpf(2,:),'-r.');plot(Xreal(2,:),'-b.')
legend('粒子滤波后的b','平均值b')
subplot(2,2,3);
hold on;box on;
plot(Xpf(3,:),'-r.');plot(Xreal(3,:),'-b.')
legend('粒子滤波后的c','平均值c')
subplot(2,2,4);
hold on;box on;
plot(Xpf(4,:),'-r.');plot(Xreal(4,:),'-b.')
legend('粒子滤波后的d','平均值d')figure
hold on;box on;
plot(Z,'-b.')
plot(Zpf,'-r.')
plot(Xf,Zf,'-g.')
bar(start,1,'y')
legend('实验测量数据','滤波估计数据','自然预测数据')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hfun.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 版权声明：
%     本程序的详细中文注释请参考
%     黄小平，王岩，缪鹏程.粒子滤波原理及应用[M].电子工业出版社，2017.4
%     书中有原理介绍+例子+程序+中文注释
%     如果此程序有错误，请对提示修改
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 函数名称：电容的观测函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=hfun(X,k)
Q=X(1)*exp(X(2)*k)+X(3)*exp(X(4)*k);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

randomR.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 版权声明：
%     本程序的详细中文注释请参考
%     黄小平，王岩，缪鹏程.粒子滤波原理及应用[M].电子工业出版社，2017.4
%     书中有原理介绍+例子+程序+中文注释
%     如果此程序有错误，请对提示修改
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function outIndex = randomR(inIndex,q);
if nargin < 2error('Not enough input arguments.');
end
outIndex=zeros(size(inIndex));
[num,col]=size(q);
u=rand(num,1);
u=sort(u);
l=cumsum(q);
i=1;
for j=1:numwhile (i<=num)&(u(i)<=l(j))outIndex(i)=j;i=i+1;end
end

锂离子电池健康状态估计（二）基于粒子滤波算法的锂电池剩余使用寿命预测，python+Matlab相关推荐

【ELM预测】基于粒子群算法PSO优化极限学习机预测含Matlab源码
1 模型为了提高空气质量预测精度,提出一种基于粒子群算法优化极限学习机的空气质量预测模型.运用粒子群算法优化极限学习机的初始权值和偏置,在保证预测误差最小的情况下实现空气质量最优预测.选择平均绝对百 ...
【回归预测-ELM预测】基于粒子群算法PSO优化极限学习机预测附matlab代码
1 内容介绍风电功率预测为电网规划提供重要的依据,研究风电功率预测方法对确保电网在安全稳定运行下接纳更多的风电具有重要的意义.针对极限学习机(ELM)回归模型预测结果受输入参数影响的问题,现将粒子群 ...
【ELAMN预测】基于粒子群算法优化ELMAN神经网络实现数据回归预测 matlab代码
1 简介风能,作为一种重要,有潜力,无污染,可再生.可持续的能源,已经成为全球发电最为迅速的能源之一,越来越受到世界各国的青睐.近年来,为缓解能源短缺问题,改善环境,实现经济乃至人类的可持续发展,世 ...
【微电网优化】基于粒子群算法求解智能微电网调度问题附matlab代码
1 简介搭建光伏,风力发电机和储能电池的数学模型.充分考虑对蓄电池的充放电保护,制定优化调度策略.应用粒子群算法(PSO)对其优化调度模型进行求解,在算法中增加了蓄电池满充满放的限制条件,同时使系统 ...
【SVM分类】基于粒子群算法优化支持向量机实现葡萄酒数据分类附matlab代码
1 简介在机器学习领域,要处理的数据的规模越来越大,而学习算法在数据的特征过多的时候,往往会产生性能上的下降.作为解决这个问题的有效手段,特征选择一直以来都得到了广泛的关注.粒子群优化算法作为一种优 ...
【优化选址】基于粒子群算法实现无线传感器网络覆盖优化附matlab代码
1 简介无线传感器网络是将大量的传感器感知节点散布在监测区域中,通过节点之间的无线信息传输形成的自组网.由于无线传感器网络工作环境复杂,传感器节点更换电源不便,网络的覆盖控制问题成为研究核心,它决定 ...
蒙特卡洛粒子滤波定位算法_基于粒子滤波的TBD算法仿真—MATLAB仿真
目标跟踪的最终目的是在最小的误差下确定目标的位置,而在无线传感器网络中要实现这个目的需要很多相关技术的支持,如定位技术.目标检测技术.估计技术.节能技术等.目标跟踪问题的求解有很多方法, 从算法的考虑 ...
蒙特卡洛粒子滤波定位算法_蒙特卡罗定位算法(基于粒子滤波的定位算法) ——原理、理解与仿真...
1 算法原理 1.1 机器人定位问题关于机器人定位,有三大问题,它们分别是: (1)"全局定位":指初始位置未知,机器人靠自身运动确定自己在地图中的位姿. (2)"位姿 ...
[转]基于粒子滤波的TBD算法仿真----MATLAB仿真
原文链接: https://blog.xxcxw.cn/2019/08/10/%e5%9f%ba%e4%ba%8e%e7%b2%92%e5%ad%90%e6%bb%a4%e6%b3%a2%e7%9a% ...
【回归预测】基于粒子滤波实现锂离子电池寿命预测附matlab代码
1 内容介绍随着现代生产生活对系统设备可靠性.安全性要求的提高,从成本.可靠性的角度考虑,电子系统正逐步由原来的定期维修变成视情维修(CBM, Condition Based Maintenanc ...

锂离子电池健康状态估计（二）基于粒子滤波算法的锂电池剩余使用寿命预测，python+Matlab