做电池管理系统、电机控制或者任何“模型里带未知参数”的工程大概率都经历过一个尴尬阶段模型方程写得明明白白但参数要么拿不准要么跑着跑着就漂了。电池内阻随温度、SOC、循环次数变化电机电感电阻随工况漂移如果只用一组固定参数去跑模型误差会越滚越大。这个时候就需要在线参数辨识。传统做法里最小二乘适合离线处理递推最小二乘能实现在线但很难处理强非线性系统扩展卡尔曼滤波EKF理论上可以上实际调试时你要手动推导雅可比矩阵线性化误差一大估计结果直接给你表演“跑偏”。无迹卡尔曼滤波UKF解决的就是这个痛点——它不需要算雅可比矩阵用一组精心挑选的sigma点去传播状态分布非线性系统的参数辨识瞬间变得简单许多。这篇文章我会把UKF做参数辨识的原理拆开讲清楚然后用锂电池一阶RC等效电路模型的参数辨识作为实战案例从建模、数据生成、代码实现到调参避坑一次性走完整个流程。适合正在做BMS算法、控制系统状态估计或者想给自己的模型加上在线自校正能力的工程师和研究生参考。1. 选型思路为什么参数辨识要用UKF而不是EKF1.1 参数辨识问题的本质就是状态估计假如我们有一组观测数据z和一个带未知参数的模型yh(x,u,θ)其中θ是待辨识的参数。所谓参数辨识就是要找一组θ让模型的预测输出尽量贴近观测值。这个问题的经典解法是最小二乘把θ当作待优化变量构建目标函数min Σ(yi-zi)^2然后求解析解或者迭代解。最小二乘在离线场景下非常好用数据攒够了一次性算出最优值。但如果你要求“在线”——也就是数据一边进来、参数一边更新——最小二乘就要升级成递推最小二乘RLS。RLS处理线性系统很高效可是现实中的系统大多是非线性的OCV随SOC的变化曲线是一条弯弯绕绕的曲线电机磁链和电流的关系也不是线性的RLS在这种场景下强行线性化效果往往不理想。另一种思路是把θ当成状态变量放进状态向量里面和系统的原始状态一起估计。这就是“状态扩维”。这样一来参数辨识就变成了一个状态估计问题而状态估计正好是卡尔曼滤波的主场模型预测一步量测更新一步参数在每一步都被“校正”一点最终收敛到真实值附近。这种思路最妙的地方在于它把“动态系统状态不可测”和“模型参数未知”这两个问题统一成了一个框架。你不需要单独设计参数更新规则只需要给参数部分加一个简单的随机游走模型θ(k1)θ(k)w让卡尔曼滤波自己去判断参数该怎么动。w的方差大说明参数变化快滤波就更愿意跟着观测去修正参数w的方差小说明参数稳定滤波就更保守。这个“信任分配”机制是整个参数辨识方法的核心。1.2 EKF、UKF、粒子滤波的取舍判断明确了状态扩维思路之后下一步就是选滤波器。最常见的在线状态估计选项有三个EKF、UKF、粒子滤波PF。EKF的思路是把非线性函数在当前估计点做一阶泰勒展开用雅可比矩阵近似线性化。它在系统非线性程度不高的时候表现不错但有两个很难受的点一是雅可比矩阵的推导非常容易出错尤其是系统复杂或者状态维数高的时候一个符号错了调试时间按小时起步二是线性化误差会直接进入协方差传播导致估计结果有偏严重时直接发散。我有一次用EKF辨识电机参数电磁方程里两个微小变量的交叉项被我漏掉了结果滤波估计值一直往一个方向偏最后排查了整整一个下午。粒子滤波的思路是用大量随机粒子去近似状态的后验分布理论上很漂亮几乎不限制系统类型但代价是计算量巨大。工程上在线跑PF动不动就要几千个粒子每个粒子都要过一遍状态方程和量测方程算力预算很容易超标而且粒子退化之后还要做重采样实现复杂度也不低。UKF正好卡在中间既不用算雅可比矩阵也不需要成千上万个粒子只需要2n1个sigma点n是状态维数就能把非线性变换后的分布均值和协方差近似到二阶精度。这个精度在绝大多数工程问题里已经够用。实践下来UKF的代码实现难度和EKF差不多但鲁棒性好很多尤其是系统非线性强、或者状态初值偏离较大的时候优势很明显。所以我的建议是系统简单、非线性弱可以继续用EKF系统非线性明显或者你根本不想推雅可比矩阵直接上UKF除非状态分布严重非高斯、或者对精度有极端要求再考虑粒子滤波。方法是否需要雅可比矩阵计算量非线性适应能力实现难度EKF需要手推易错低弱线性化误差大中UKF不需要中强二阶精度低粒子滤波不需要高最强高1.3 无迹变换到底在算什么要理解UKF必须先理解无迹变换Unscented Transform。简单说sigma点是一组精心选取的确定性采样点。你不是随机撒点而是根据当前状态均值和协方差对称地在均值周围摆上2n1个点。这些点经过非线性变换之后用加权平均的方式算出变换后分布的均值和协方差。这个过程可以理解成派几个侦察兵先穿过非线性函数然后根据侦察兵带回来的信息推断整个分布变成了什么样。这里有一个很关键的对比EKF是“用一个点均值去近似整个分布然后把函数线性化”UKF是“用一组点去近似分布然后让完整的非线性函数直接作用在这些点上”。后者没有做任何线性化近似所以它天然适用于强非线性场景。sigma点的选取是确定性的不带随机性不会像粒子滤波那样受采样噪声影响所以计算量小、结果稳定。这也就是为什么UKF在参数辨识上这么顺手——参数和状态之间往往是非线性的耦合关系无迹变换正好把这一层非线性绕过去了。2. 核心代码骨架从sigma点到预测更新2.1 状态扩维给模型“装”上待辨识参数实现层面第一步就是状态扩维。假设原系统状态是x维度n_x要辨识的参数维度n_θ那么扩维之后的状态维度是nn_xn_θ。在UKF里所有sigma点、协方差矩阵、过程噪声矩阵都要跟着扩展到n维。我以电池一阶RC模型为例这个模型的状态变量是SOC和极化电压U1待辨识参数是欧姆内阻R0、极化内阻R1、极化电容C1扩维后状态向量就是[SOC, U1, R0, R1, C1]^T。初始协方差P0里参数对应的对角元素要设得比状态部分大一些意思是对参数的初始估计不太自信让滤波器有足够的自由度去搜索真实值。这里要特别注意量纲差异。C1的量级是干法拉约1000F它的方差设置就要比R0约0.05Ω大很多否则滤波器一动就超出合理范围数值上容易出问题。如果你把不同物理量的方差都设成同一个值UKF几乎必发散原因是Cholesky分解对矩阵正定性要求较高量纲差异过大时矩阵条件数爆炸数值不稳定。实践中的做法是分别设置让每个对角元素跟对应状态量的量级匹配。import numpy as np # 扩维状态: [SOC, U1, R0, R1, C1] x0 np.array([0.85, 0.0, 0.08, 0.05, 800.0]) P0 np.diag([1e-2, 1e-4, 1e-4, 1e-4, 1e4])2.2 sigma点采样与权重计算sigma点采样的标准公式是这样的设状态维度为n缩放参数λα²(nκ)-n其中α决定sigma点离均值的距离κ是次级缩放参数β和先验分布相关。对于高斯分布典型取值为α1e-3、κ0、β2。为什么β取2因为在高斯假设下这个取值能让协方差估计达到最优精度。然后生成2n1个sigma点第一个是均值本身剩下2n个点沿协方差矩阵的Cholesky分解方向对称分布。权重分两组Wm用于计算均值Wc用于计算协方差第一个点的权重和其他点不同这是为了保证加权后的协方差无偏。alpha 1e-3 beta 2.0 kappa 0.0 n len(x0) lam alpha**2 * (n kappa) - n Wm np.zeros(2*n 1) Wc np.zeros(2*n 1) Wm[0] lam / (n lam) Wc[0] lam / (n lam) (1 - alpha**2 beta) Wm[1:] 1 / (2 * (n lam)) Wc[1:] 1 / (2 * (n lam)) def sigma_points(x, P, lam): L np.linalg.cholesky((n lam) * P) chi np.zeros((2*n 1, n)) chi[0] x for i in range(n): chi[i1] x L[i] chi[ni1] x - L[i] return chialpha的取值很有意思取1e-3意味着sigma点非常靠近均值这在大噪声、强非线性系统里反而有可能让采样点过于集中丢失尾部信息。我见过一些工程实现把alpha调到0.1甚至0.5效果反而更稳。但从理论角度alpha越小高阶项误差控制越好。我的习惯是先用默认值如果发现滤波收敛慢就把alpha往大调一档试试。kappa取0时λ可能是负的当α很小时但nλ通常会保持正数只需在实现里加个判断即可。2.3 预测与更新两步循环的实现细节UKF主循环分成两步预测时间更新和校正量测更新。预测步做的事情是把每个sigma点都扔进状态方程f得到传播后的sigma点集合然后用权重加权得到先验状态均值和先验协方差。注意这里Q矩阵一定要加进去它表示模型本身的不确定性。校正步稍微复杂一点要用先验均值和协方差重新生成一组sigma点再把这组点扔进量测方程h得到预测的量测值然后计算量测协方差S和互协方差Pxz最终算出卡尔曼增益K完成状态和协方差的更新。为什么校正步要重新生成sigma点而不是直接用传播后的那组严格来说UKF算法要求在每个阶段都用当前最新的高斯分布生成新的sigma点这样能保证量测更新的统计一致性。工程上有些实现会偷懒直接复用传播后的sigma点结果差别通常不大但严谨起见我还是按标准流程写。def f_func(chi, I, dt, Q_bat): chi_next np.zeros_like(chi) for i in range(chi.shape[0]): soc, u1, r0, r1, c1 chi[i] soc_next soc - I * dt / (3600.0 * Q_bat) u1_next np.exp(-dt/(r1*c1)) * u1 r1*(1 - np.exp(-dt/(r1*c1))) * I # 参数随机游走保持不变 chi_next[i] [soc_next, u1_next, r0, r1, c1] return chi_next def h_func(chi, I): soc, u1, r0, r1, c1 chi.T return ocv(soc) - r0 * I - u1 def ukf_predict(chi, Wm, Wc, Q): x_pred chi Wm d chi - x_pred P_pred (Wc[:, None] * d).T d Q return x_pred, P_pred def ukf_update(x_pred, P_pred, z, I, Wm, Wc, R): chi2 sigma_points(x_pred, P_pred, lam) Z h_func(chi2, I) z_pred Z Wm dZ Z - z_pred S (Wc[:, None] * dZ).T dZ R dX chi2 - x_pred Pxz (Wc[:, None] * dX).T dZ K Pxz np.linalg.inv(S) x_new x_pred K * (z - z_pred) P_new P_pred - K S K.T P_new (P_new P_new.T) / 2 # 强制对称化防止数值误差累积 return x_new, P_newP_new的对称化处理是我每次必写的一行。反复迭代之后矩阵乘法带来的舍入误差会让P阵逐渐失去对称性对称阵一旦不对称后面的Cholesky分解迟早会报错。顺手加上这一行能省掉很多排查时间。另外在更新步如果遇到S矩阵接近奇异可以在S上加一个极小的对角阵比如1e-12 * I来保证可逆这也是工程上常用的技巧。3. 实战锂电池一阶RC模型参数在线辨识3.1 Thevenin模型离散化与工况数据生成电池参数辨识最经典的场景之一是一阶RC等效电路模型Thevenin模型。这个模型包含一个开路电压源OCV(SOC)、一个欧姆内阻R0和一个RC并联网络R1、C1它能够描述电池的端电压瞬降和极化弛豫特性。模型方程如下端电压Ut OCV(SOC) - R0 * I - U1极化电压动态dU1/dt -U1/(R1*C1) I/C1离散化之后用指数积分形式可以写成SOC(k1) SOC(k) - I(k) * dt / (3600 * Q_bat)U1(k1) exp(-dt/(R1C1)) * U1(k) R1(1 - exp(-dt/(R1*C1))) * I(k)其中Q_bat是电池容量Ahdt是采样周期。R1*C1是极化时间常数这个值直接决定了电池电压的“弛豫速度”。OCV与SOC的关系一般通过实验标定我用一个三次多项式近似对示例来说足够了。生成仿真数据的时候我特意设计了一个多频叠加电流工况而不是简单的恒流或者单频正弦。为什么因为参数辨识需要“持续激励”——输入的电流信号要足够丰富让系统动态信息都暴露出来。恒流工况下U1进入稳态R1和C1的信息就全部丢失了单频正弦只能激励到某个频段多个频率叠加则能让极化过程在多个时间尺度上都被激活。这一点在后面调参时体会会更深。T 1500 dt 1.0 t np.arange(T) I (2.0*np.sin(0.005*t) 1.5*np.sin(0.02*t) 1.0*np.sin(0.08*t) np.random.normal(0, 0.2, T)) Q_bat 7.0 R0_true, R1_true, C1_true 0.05, 0.03, 1000.0 def ocv(soc): return 3.0 1.2*soc - 0.8*soc**2 0.4*soc**3 SOC_true np.zeros(T) U1_true np.zeros(T) V_true np.zeros(T) SOC_true[0] 0.9 U1_true[0] 0.0 for k in range(T-1): SOC_true[k1] SOC_true[k] - I[k] * dt / (3600.0 * Q_bat) U1_true[k1] (np.exp(-dt/(R1_true*C1_true)) * U1_true[k] R1_true*(1 - np.exp(-dt/(R1_true*C1_true))) * I[k]) for k in range(T): V_true[k] ocv(SOC_true[k]) - R0_true * I[k] - U1_true[k] V_meas V_true np.random.normal(0, 0.01, T)3.2 滤波参数初始化P0、Q、R的取值逻辑滤波器能不能收敛一半的功夫在初始化上。P0、Q、R三个矩阵分别代表初始状态不确定性、过程噪声方差、量测噪声方差。先说P0我给了[SOC, U1, R0, R1, C1] [1e-2, 1e-4, 1e-4, 1e-4, 1e4]。SOC初始给0.85而真实值是0.9偏差0.05方差1e-2意味着标准差0.1留了足够裕度U1初始0给1e-4是差不多的量级R0给0.08而真实值0.05偏差0.03方差1e-4对应标准差0.01稍微偏小但还能接受C1给800而真实值1000偏差200方差1e4对应标准差100这个裕度就合理了。Q矩阵表示模型不确定性。状态部分给1e-6左右因为SOC和U1的状态方程我比较信任主要误差来自工况随机波动参数部分每一项都不同R0和R1的过程噪声给1e-8C1给1e-2。这里最容易被忽略的是量纲匹配——C1是千法拉级如果它的过程噪声和R0一样是1e-8那滤波器几乎不会去更新C1因为它觉得C1很“稳定”结果就是C1永远停在初值800附近。把C1的Q调大到1e-2标准差0.1滤波器才愿意在合理范围内移动C1。R就是量测噪声端电压的测量噪声标准差大约0.01V所以R给1
做电池管理系统、电机控制或者任何“模型里带未知参数”的工程大概率都经历过一个尴尬阶段模型方程写得明明白白但参数要么拿不准要么跑着跑着就漂了。电池内阻随温度、SOC、循环次数变化电机电感电阻随工况漂移如果只用一组固定参数去跑模型误差会越滚越大。这个时候就需要在线参数辨识。传统做法里最小二乘适合离线处理递推最小二乘能实现在线但很难处理强非线性系统扩展卡尔曼滤波EKF理论上可以上实际调试时你要手动推导雅可比矩阵线性化误差一大估计结果直接给你表演“跑偏”。无迹卡尔曼滤波UKF解决的就是这个痛点——它不需要算雅可比矩阵用一组精心挑选的sigma点去传播状态分布非线性系统的参数辨识瞬间变得简单许多。这篇文章我会把UKF做参数辨识的原理拆开讲清楚然后用锂电池一阶RC等效电路模型的参数辨识作为实战案例从建模、数据生成、代码实现到调参避坑一次性走完整个流程。适合正在做BMS算法、控制系统状态估计或者想给自己的模型加上在线自校正能力的工程师和研究生参考。1. 选型思路为什么参数辨识要用UKF而不是EKF1.1 参数辨识问题的本质就是状态估计假如我们有一组观测数据z和一个带未知参数的模型yh(x,u,θ)其中θ是待辨识的参数。所谓参数辨识就是要找一组θ让模型的预测输出尽量贴近观测值。这个问题的经典解法是最小二乘把θ当作待优化变量构建目标函数min Σ(yi-zi)^2然后求解析解或者迭代解。最小二乘在离线场景下非常好用数据攒够了一次性算出最优值。但如果你要求“在线”——也就是数据一边进来、参数一边更新——最小二乘就要升级成递推最小二乘RLS。RLS处理线性系统很高效可是现实中的系统大多是非线性的OCV随SOC的变化曲线是一条弯弯绕绕的曲线电机磁链和电流的关系也不是线性的RLS在这种场景下强行线性化效果往往不理想。另一种思路是把θ当成状态变量放进状态向量里面和系统的原始状态一起估计。这就是“状态扩维”。这样一来参数辨识就变成了一个状态估计问题而状态估计正好是卡尔曼滤波的主场模型预测一步量测更新一步参数在每一步都被“校正”一点最终收敛到真实值附近。这种思路最妙的地方在于它把“动态系统状态不可测”和“模型参数未知”这两个问题统一成了一个框架。你不需要单独设计参数更新规则只需要给参数部分加一个简单的随机游走模型θ(k1)θ(k)w让卡尔曼滤波自己去判断参数该怎么动。w的方差大说明参数变化快滤波就更愿意跟着观测去修正参数w的方差小说明参数稳定滤波就更保守。这个“信任分配”机制是整个参数辨识方法的核心。1.2 EKF、UKF、粒子滤波的取舍判断明确了状态扩维思路之后下一步就是选滤波器。最常见的在线状态估计选项有三个EKF、UKF、粒子滤波PF。EKF的思路是把非线性函数在当前估计点做一阶泰勒展开用雅可比矩阵近似线性化。它在系统非线性程度不高的时候表现不错但有两个很难受的点一是雅可比矩阵的推导非常容易出错尤其是系统复杂或者状态维数高的时候一个符号错了调试时间按小时起步二是线性化误差会直接进入协方差传播导致估计结果有偏严重时直接发散。我有一次用EKF辨识电机参数电磁方程里两个微小变量的交叉项被我漏掉了结果滤波估计值一直往一个方向偏最后排查了整整一个下午。粒子滤波的思路是用大量随机粒子去近似状态的后验分布理论上很漂亮几乎不限制系统类型但代价是计算量巨大。工程上在线跑PF动不动就要几千个粒子每个粒子都要过一遍状态方程和量测方程算力预算很容易超标而且粒子退化之后还要做重采样实现复杂度也不低。UKF正好卡在中间既不用算雅可比矩阵也不需要成千上万个粒子只需要2n1个sigma点n是状态维数就能把非线性变换后的分布均值和协方差近似到二阶精度。这个精度在绝大多数工程问题里已经够用。实践下来UKF的代码实现难度和EKF差不多但鲁棒性好很多尤其是系统非线性强、或者状态初值偏离较大的时候优势很明显。所以我的建议是系统简单、非线性弱可以继续用EKF系统非线性明显或者你根本不想推雅可比矩阵直接上UKF除非状态分布严重非高斯、或者对精度有极端要求再考虑粒子滤波。方法是否需要雅可比矩阵计算量非线性适应能力实现难度EKF需要手推易错低弱线性化误差大中UKF不需要中强二阶精度低粒子滤波不需要高最强高1.3 无迹变换到底在算什么要理解UKF必须先理解无迹变换Unscented Transform。简单说sigma点是一组精心选取的确定性采样点。你不是随机撒点而是根据当前状态均值和协方差对称地在均值周围摆上2n1个点。这些点经过非线性变换之后用加权平均的方式算出变换后分布的均值和协方差。这个过程可以理解成派几个侦察兵先穿过非线性函数然后根据侦察兵带回来的信息推断整个分布变成了什么样。这里有一个很关键的对比EKF是“用一个点均值去近似整个分布然后把函数线性化”UKF是“用一组点去近似分布然后让完整的非线性函数直接作用在这些点上”。后者没有做任何线性化近似所以它天然适用于强非线性场景。sigma点的选取是确定性的不带随机性不会像粒子滤波那样受采样噪声影响所以计算量小、结果稳定。这也就是为什么UKF在参数辨识上这么顺手——参数和状态之间往往是非线性的耦合关系无迹变换正好把这一层非线性绕过去了。2. 核心代码骨架从sigma点到预测更新2.1 状态扩维给模型“装”上待辨识参数实现层面第一步就是状态扩维。假设原系统状态是x维度n_x要辨识的参数维度n_θ那么扩维之后的状态维度是nn_xn_θ。在UKF里所有sigma点、协方差矩阵、过程噪声矩阵都要跟着扩展到n维。我以电池一阶RC模型为例这个模型的状态变量是SOC和极化电压U1待辨识参数是欧姆内阻R0、极化内阻R1、极化电容C1扩维后状态向量就是[SOC, U1, R0, R1, C1]^T。初始协方差P0里参数对应的对角元素要设得比状态部分大一些意思是对参数的初始估计不太自信让滤波器有足够的自由度去搜索真实值。这里要特别注意量纲差异。C1的量级是干法拉约1000F它的方差设置就要比R0约0.05Ω大很多否则滤波器一动就超出合理范围数值上容易出问题。如果你把不同物理量的方差都设成同一个值UKF几乎必发散原因是Cholesky分解对矩阵正定性要求较高量纲差异过大时矩阵条件数爆炸数值不稳定。实践中的做法是分别设置让每个对角元素跟对应状态量的量级匹配。import numpy as np # 扩维状态: [SOC, U1, R0, R1, C1] x0 np.array([0.85, 0.0, 0.08, 0.05, 800.0]) P0 np.diag([1e-2, 1e-4, 1e-4, 1e-4, 1e4])2.2 sigma点采样与权重计算sigma点采样的标准公式是这样的设状态维度为n缩放参数λα²(nκ)-n其中α决定sigma点离均值的距离κ是次级缩放参数β和先验分布相关。对于高斯分布典型取值为α1e-3、κ0、β2。为什么β取2因为在高斯假设下这个取值能让协方差估计达到最优精度。然后生成2n1个sigma点第一个是均值本身剩下2n个点沿协方差矩阵的Cholesky分解方向对称分布。权重分两组Wm用于计算均值Wc用于计算协方差第一个点的权重和其他点不同这是为了保证加权后的协方差无偏。alpha 1e-3 beta 2.0 kappa 0.0 n len(x0) lam alpha**2 * (n kappa) - n Wm np.zeros(2*n 1) Wc np.zeros(2*n 1) Wm[0] lam / (n lam) Wc[0] lam / (n lam) (1 - alpha**2 beta) Wm[1:] 1 / (2 * (n lam)) Wc[1:] 1 / (2 * (n lam)) def sigma_points(x, P, lam): L np.linalg.cholesky((n lam) * P) chi np.zeros((2*n 1, n)) chi[0] x for i in range(n): chi[i1] x L[i] chi[ni1] x - L[i] return chialpha的取值很有意思取1e-3意味着sigma点非常靠近均值这在大噪声、强非线性系统里反而有可能让采样点过于集中丢失尾部信息。我见过一些工程实现把alpha调到0.1甚至0.5效果反而更稳。但从理论角度alpha越小高阶项误差控制越好。我的习惯是先用默认值如果发现滤波收敛慢就把alpha往大调一档试试。kappa取0时λ可能是负的当α很小时但nλ通常会保持正数只需在实现里加个判断即可。2.3 预测与更新两步循环的实现细节UKF主循环分成两步预测时间更新和校正量测更新。预测步做的事情是把每个sigma点都扔进状态方程f得到传播后的sigma点集合然后用权重加权得到先验状态均值和先验协方差。注意这里Q矩阵一定要加进去它表示模型本身的不确定性。校正步稍微复杂一点要用先验均值和协方差重新生成一组sigma点再把这组点扔进量测方程h得到预测的量测值然后计算量测协方差S和互协方差Pxz最终算出卡尔曼增益K完成状态和协方差的更新。为什么校正步要重新生成sigma点而不是直接用传播后的那组严格来说UKF算法要求在每个阶段都用当前最新的高斯分布生成新的sigma点这样能保证量测更新的统计一致性。工程上有些实现会偷懒直接复用传播后的sigma点结果差别通常不大但严谨起见我还是按标准流程写。def f_func(chi, I, dt, Q_bat): chi_next np.zeros_like(chi) for i in range(chi.shape[0]): soc, u1, r0, r1, c1 chi[i] soc_next soc - I * dt / (3600.0 * Q_bat) u1_next np.exp(-dt/(r1*c1)) * u1 r1*(1 - np.exp(-dt/(r1*c1))) * I # 参数随机游走保持不变 chi_next[i] [soc_next, u1_next, r0, r1, c1] return chi_next def h_func(chi, I): soc, u1, r0, r1, c1 chi.T return ocv(soc) - r0 * I - u1 def ukf_predict(chi, Wm, Wc, Q): x_pred chi Wm d chi - x_pred P_pred (Wc[:, None] * d).T d Q return x_pred, P_pred def ukf_update(x_pred, P_pred, z, I, Wm, Wc, R): chi2 sigma_points(x_pred, P_pred, lam) Z h_func(chi2, I) z_pred Z Wm dZ Z - z_pred S (Wc[:, None] * dZ).T dZ R dX chi2 - x_pred Pxz (Wc[:, None] * dX).T dZ K Pxz np.linalg.inv(S) x_new x_pred K * (z - z_pred) P_new P_pred - K S K.T P_new (P_new P_new.T) / 2 # 强制对称化防止数值误差累积 return x_new, P_newP_new的对称化处理是我每次必写的一行。反复迭代之后矩阵乘法带来的舍入误差会让P阵逐渐失去对称性对称阵一旦不对称后面的Cholesky分解迟早会报错。顺手加上这一行能省掉很多排查时间。另外在更新步如果遇到S矩阵接近奇异可以在S上加一个极小的对角阵比如1e-12 * I来保证可逆这也是工程上常用的技巧。3. 实战锂电池一阶RC模型参数在线辨识3.1 Thevenin模型离散化与工况数据生成电池参数辨识最经典的场景之一是一阶RC等效电路模型Thevenin模型。这个模型包含一个开路电压源OCV(SOC)、一个欧姆内阻R0和一个RC并联网络R1、C1它能够描述电池的端电压瞬降和极化弛豫特性。模型方程如下端电压Ut OCV(SOC) - R0 * I - U1极化电压动态dU1/dt -U1/(R1*C1) I/C1离散化之后用指数积分形式可以写成SOC(k1) SOC(k) - I(k) * dt / (3600 * Q_bat)U1(k1) exp(-dt/(R1C1)) * U1(k) R1(1 - exp(-dt/(R1*C1))) * I(k)其中Q_bat是电池容量Ahdt是采样周期。R1*C1是极化时间常数这个值直接决定了电池电压的“弛豫速度”。OCV与SOC的关系一般通过实验标定我用一个三次多项式近似对示例来说足够了。生成仿真数据的时候我特意设计了一个多频叠加电流工况而不是简单的恒流或者单频正弦。为什么因为参数辨识需要“持续激励”——输入的电流信号要足够丰富让系统动态信息都暴露出来。恒流工况下U1进入稳态R1和C1的信息就全部丢失了单频正弦只能激励到某个频段多个频率叠加则能让极化过程在多个时间尺度上都被激活。这一点在后面调参时体会会更深。T 1500 dt 1.0 t np.arange(T) I (2.0*np.sin(0.005*t) 1.5*np.sin(0.02*t) 1.0*np.sin(0.08*t) np.random.normal(0, 0.2, T)) Q_bat 7.0 R0_true, R1_true, C1_true 0.05, 0.03, 1000.0 def ocv(soc): return 3.0 1.2*soc - 0.8*soc**2 0.4*soc**3 SOC_true np.zeros(T) U1_true np.zeros(T) V_true np.zeros(T) SOC_true[0] 0.9 U1_true[0] 0.0 for k in range(T-1): SOC_true[k1] SOC_true[k] - I[k] * dt / (3600.0 * Q_bat) U1_true[k1] (np.exp(-dt/(R1_true*C1_true)) * U1_true[k] R1_true*(1 - np.exp(-dt/(R1_true*C1_true))) * I[k]) for k in range(T): V_true[k] ocv(SOC_true[k]) - R0_true * I[k] - U1_true[k] V_meas V_true np.random.normal(0, 0.01, T)3.2 滤波参数初始化P0、Q、R的取值逻辑滤波器能不能收敛一半的功夫在初始化上。P0、Q、R三个矩阵分别代表初始状态不确定性、过程噪声方差、量测噪声方差。先说P0我给了[SOC, U1, R0, R1, C1] [1e-2, 1e-4, 1e-4, 1e-4, 1e4]。SOC初始给0.85而真实值是0.9偏差0.05方差1e-2意味着标准差0.1留了足够裕度U1初始0给1e-4是差不多的量级R0给0.08而真实值0.05偏差0.03方差1e-4对应标准差0.01稍微偏小但还能接受C1给800而真实值1000偏差200方差1e4对应标准差100这个裕度就合理了。Q矩阵表示模型不确定性。状态部分给1e-6左右因为SOC和U1的状态方程我比较信任主要误差来自工况随机波动参数部分每一项都不同R0和R1的过程噪声给1e-8C1给1e-2。这里最容易被忽略的是量纲匹配——C1是千法拉级如果它的过程噪声和R0一样是1e-8那滤波器几乎不会去更新C1因为它觉得C1很“稳定”结果就是C1永远停在初值800附近。把C1的Q调大到1e-2标准差0.1滤波器才愿意在合理范围内移动C1。R就是量测噪声端电压的测量噪声标准差大约0.01V所以R给1
