☰
电池连续时间模型与一阶RC等效电路参数估计实战
2026/10/15 5:48:32 网站建设 项目流程

1. 先吃透第一问:这道题真正要交的是“连续时间模型”,不是一个查表程序

很多人拿到智能手机电池的放电数据表,第一反应是拉一条平滑曲线,然后说“看,拟合得挺好的”。这种做法在建模竞赛题里基本拿不到高分,因为题目关键词已经写得很明确了:连续时间建模。它要你做的不是把电压和时间的关系用多项式硬凑出来,而是把电池内部那段无法直接观测的物理过程,用一组连续变化的微分方程描述出来,再通过数据反解参数。

先看第一问最常见的问法:给了某型号智能手机电池在恒流放电条件下的电压-时间记录,要求建立连续时间模型,解释电压下降过程,并给出关键参数。这里的“连续时间”有很实际的工程意义——你手里只有每隔几秒或几分钟采样一次的数据,但你希望回答任意时刻 t 的电压是多少,剩余电量还剩多少,什么时候会掉到截止电压。这是离散查表和连续模型最本质的区别:离散表只能告诉你“记录过的时刻”,连续模型能回答“没记录过的时刻”。

所以第一问的交付物,我建议拆成四块:

  • 一套状态方程或者微分方程,明确写出状态变量是什么(通常至少包含荷电状态 SoC 和一个极化电压分量);
  • 观测方程,把不可直接测量的状态变量和可以直接测量的端电压 V(t) 连起来;
  • 基于实测数据的参数估计结果,比如欧姆内阻、极化电阻、时间常数、开路电压曲线系数;
  • 模型验证,至少要画出拟合曲线与原始数据的对照图,并计算均方根误差 RMSE。

很多队伍第一问只交第三块和第四块,漏掉最核心的第一块和第二块,这会让整道题的逻辑链断掉。后续的问法一旦要求预测剩余使用寿命、讨论温度影响,如果没有状态方程打底,基本就接不上去了。

1.1 连续模型和“逐点查表”到底差在哪

举一个很朴素的理解方式。假设你手里有一组每分钟记录一次的电压表,电池从 4.2V 开始放电,两个小时后到 3.4V。如果只用插值或者分段线性拟合,你在第 61 秒、第 123 秒这些非采样时刻给的电压值其实没有任何物理含义,纯粹是数学插值的产物。

连续时间模型则不同。你假设电压的变化不是无缘无故发生的,而是由电池内部的状态——例如剩余电量下降、极化电压建立、内阻分压变化——共同推动的。这些状态随时间连续演化,所以电压随时间连续演化。你把演化规则写清楚,剩下的工作就是定量估计演化规则里的系数。

从数学形式上看,离散数据经常写成时间序列 ({t_k, V_k}),连续时间模型则写成:

[ \frac{dx}{dt} = f(x, u, \theta), \quad V = g(x, u, \theta) ]

其中 (x) 是状态向量,(u) 是外部输入,比如放电电流和工作温度,(\theta) 是待定参数。第一问的难点不在 f 和 g 写不出来,而在怎么用真实数据把 (\theta) 估准。这也是本文下面要重点展开的地方。

1.2 第一问常见的模型框架

针对智能手机锂电池,实用且不过度复杂的方案是一阶 RC 等效电路模型。它把电池看成三部分组成:一个开路电压源 OCV,描述电池平衡状态下的电压;一个串联欧姆内阻 R0,描述电流瞬间变化时立刻产生的压降;一组并联的极化电阻 R1 和极化电容 C1,描述电流变化后电压逐步趋稳的过程。

相比把电压直接写成时间的高阶多项式,一阶 RC 模型的好处是每一项都有物理对应物,评审问起来你能说清楚R0为什么存在、C1的物理含义是什么。相比复杂到十几个参数的电化学模型,一阶 RC 又足够简洁,用题目给的数据就能稳定辨识。

2. 数学化路径:从等效电路到状态方程,再到可拟合的解析形式

这一节是整个第一问的灵魂。我先把完整的一阶 RC 模型写出来,再讲怎么在连续时间框架下化简和求解。

先定义荷电状态 SoC,用 (s) 表示,取值范围 0 到 1。额定容量记为 (Q_n),单位是安时,放电电流为 (I(t)),单位是安。按电流方向定义,放电时电流为正。那么 SoC 的变化率是:

[ \frac{ds}{dt} = -\frac{I(t)}{Q_n} ]

这里需要特别注意单位。如果时间用小时,容量用安时,电流用安,等式两边自然匹配;如果 t 用秒,右边就要除以 3600,否则 SoC 下降速度会快得离谱。这个坑我放在后面详细说,代码里也会处理。

极化电压 (U_1) 描述 RC 支路两端的电压,它的微分方程来自 RC 电路的充放电规律:

[ \frac{dU_1}{dt} = \frac{I(t)}{C_1} - \frac{U_1}{R_1 C_1} ]

当电流发生阶跃时,(U_1) 会按指数规律趋近于稳态值 (I R_1),时间常数为 (\tau = R_1 C_1)。这就是电池电压在电流跳变后会“慢慢到位”的原因。

观测方程把状态量和端电压连起来:

[ V(t) = \text{OCV}(s(t)) - I(t) R_0 - U_1(t) ]

开路电压 OCV 是 SoC 的非线性函数。工程上常用多项式近似,最常见的做法是:

[ \text{OCV}(s) = a_0 + a_1 s + a_2 s^2 + a_3 s^3 ]

如果你手里的数据覆盖了从满电到放空的全过程,那么这个三阶多项式通常够用。但如果只是第一问的一小段放电数据,拟合太多多项式系数会导致严重过拟合,我建议第一问先用线性简化:

[ \text{OCV}(s) = V_0 + k(s - 1) ]

它的含义是:开路电压围绕满电电压 (V_0) 随 SoC 近似线性下降,斜率 k 通常为负值。真实电池在小范围 SoC 区间内,这个近似是合理的。

2.1 恒流放电下的解析解

第一问如果给定的是恒定电流放电数据,那么状态方程化简非常舒服。此时 (I(t) = I) 是常数,SoC 的方程可以直接积分:

[ s(t) = 1 - \frac{I}{Q_n} t ]

这里的 t 单位必须与 Q_n、I 匹配。在完整模型里,(U_1) 的微分方程是一阶线性常系数方程,从零初始状态出发的解是:

[ U_1(t) = I R_1 \left( 1 - e^{-t/\tau} \right) ]

于是端电压可以写成显式解析式:

[ V(t) = V_0 + k \left( -\frac{I}{Q_n} t \right) - I R_0 - I R_1 \left( 1 - e^{-t/\tau} \right) ]

注意这里我把 (s-1) 代入后得到的就是 (-I t / Q_n)。所以最终拟合参数为:(V_0)、(k)、(R_0)、(R_1)、(\tau)。这个式子看起来不起眼,但它把连续时间模型转化成了一个可以直接进行最小二乘拟合的形式,而且每个参数都有明确物理意义。

2.2 参数可辨识性:单独一条曲线能认出哪些参数

这里有个很多队伍想不到的点。观察上面的解析式,(R_1) 和 (C_1) 从不单独出现,只以乘积 (\tau = R_1 C_1) 的形式出现。这意味着:如果你只有一条恒流放电曲线,你最多只能稳定估计出时间常数 (\tau),无法把极化电阻和极化电容分开识别。强行同时拟合 R1 和 C1,结果会非常不稳定,出现一个参数很大、另一个参数很小的“补偿效应”。

这不是数学软件的锅,是物理模型的固有限制。要让 R1 和 C1 分别可辨识,通常需要两类额外信息:一是多组不同电流的放电数据,因为 RC 支路的响应幅度和电流大小相关,而时间常数与电流无关;二是动态电流工况下电压的瞬态响应,例如脉冲放电,可以从初始压降和恢复过程分别提取 R0 和 RC 参数。

竞赛题通常不会在第一问就逼你解决这个识别问题,但写论文时把这一点讲清楚,会显得你对模型的理解比单纯调库深一层。我在自己处理第一问时,通常直接拟合五参数集 ({V_0, k, R_0, R_1, \tau}),在正文里明确说明:由于数据是单电流恒流放电,极化电容的绝对数值不可单独辨识,后续如有多工况数据再进行扩展。

3. 第一问 Python 实现:数据准备、解析拟合、ODE 验证三件套

现在上代码。下面的代码以“一张从题面提取的放电记录表”为输入,表格至少包含两列:时间 t(单位统一为秒)和电压 V(单位统一为伏)。如果题面给的是分钟或者小时,读取后先转换。

为了演示流程,我先用已知参数生成一条带噪声的模拟放电曲线,再当作“题面数据”来做参数估计。真实比赛时,把模拟数据替换成你的数据文件即可。

3.1 生成模拟数据并统一单位

import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.optimize import curve_fit # 电池额定参数,来自题面或电池规格 Qn = 3000.0 / 1000.0 # 额定容量,单位 Ah,这里假设 3000 mAh I = 1.5 # 恒流放电电流,单位 A V_cutoff = 3.0 # 放电截止电压,单位 V # 真实模型参数(用于生成模拟观测) true_params = { "V0": 4.20, "k": -0.45, "R0": 0.08, "R1": 0.03, "tau": 800.0, } # 模拟采样:从 0 到 7200 秒,每 2 秒采一个点 t = np.linspace(0, 7200, 3601) def discharge_voltage(t, V0, k, I, Qn, R0, R1, tau): # 注意 t 以秒为单位,Qn 用 Ah,I 用 A,SoC 变化需要换算成小时 soc = 1.0 - I / Qn * (t / 3600.0) u1 = I * R1 * (1.0 - np.exp(-t / tau)) return V0 + k * (soc - 1.0) - I * R0 - u1 V_true = discharge_voltage( t, true_params["V0"], true_params["k"], I, Qn, true_params["R0"], true_params["R1"], true_params["tau"], ) # 加一点观测噪声,模拟真实采样误差 rng = np.random.default_rng(42) V_obs = V_true + rng.normal(0, 0.005, size=len(t)) df = pd.DataFrame({"t_s": t, "V": V_obs}) df.to_csv("discharge_data.csv", index=False)

把数据存成discharge_data.csv后,读取代码就是标准的 pandas 操作:

df = pd.read_csv("discharge_data.csv") t_data = df["t_s"].values V_data = df["V"].values

这里有一个很多人容易忽略的细节:模拟数据的 SoC 不能出现负值。7200 秒等于 2 小时,放电电流 1.5A,容量 3Ah,恰好能把电池放空到 SoC=0。如果你的数据时间更长,就要把放电截断那一段滤掉,否则 SoC 变成负值而解析式还继续算,会严重扭曲拟合结果。

3.2 方法一:用解析式做最小二乘参数估计

连续时间模型的好处在这里体现出来——恒流放电下解析式可以直接喂给curve_fit:

# 注意这里 q_n 和 i_val 被固定为已知量,不参与拟合 def fit_func(t, V0, k, R0, R1, tau): return discharge_voltage(t, V0, k, I, Qn, R0, R1, tau) p0 = [4.1, -0.5, 0.06, 0.05, 900.0] # 初始猜测 bounds = ( [3.5, -2.0, 0.0, 0.0, 10.0], # 下界 [4.6, -0.05, 0.2, 0.5, 5000.0], # 上界 ) popt, pcov = curve_fit( fit_func, t_data, V_data, p0=p0, bounds=bounds, maxfev=20000, ) V0_hat, k_hat, R0_hat, R1_hat, tau_hat = popt print("V0 =", V0_hat) print("k =", k_hat) print("R0 =", R0_hat) print("R1 =", R1_hat) print("tau =", tau_hat)

观察输出参数会发现,只要初值不是太离谱,估计值能回到真实参数附近。这里的关键是边界约束,特别是 R0、R1 的下界必须设为 0。物理上电阻不可能为负,如果不加这个约束,最小二乘很可能为了吸收末端曲线形状而给出负电阻,这在论文里根本说不过去。

拟合完成后,评价指标按下面方式计算:

V_pred = fit_func(t_data, *popt) residual = V_data - V_pred rmse = np.sqrt(np.mean(residual ** 2)) r2 = 1 - np.sum(residual ** 2) / np.sum((V_data - np.mean(V_data)) ** 2) print("RMSE =", rmse) print("R^2 =", r2) plt.figure(figsize=(8, 4)) plt.plot(t_data / 60.0, V_data, ".", markersize=1, label="observed") plt.plot(t_data / 60.0, V_pred, "-", linewidth=2, label="fitted") plt.xlabel("time (min)") plt.ylabel("voltage (V)") plt.legend() plt.grid(True) plt.show()

RMSE 通常要做到毫伏量级才算合理。如果你的残差出现系统性 U 形或者波浪形,说明模型结构不够,后面要往更高阶 OCV 多项式或者二阶 RC 方向走。

3.3 方法二:用 solve_ivp 验证完整连续模型

解析式只是恒流的特例。为了证明你建立的确实是“连续时间模型”,也为了方便后续处理非恒流工况,我建议再加一段基于solve_ivp的 ODE 求解。这既是一道保险,也能在论文里展示你的模型可以延展到动态电流剖面。

from scipy.integrate import solve_ivp def battery_ode(t, y, I_profile, t_profile, Qn, R0, R1, C1, ocv_coeff): s, u1 = y # 如果电流随时间变化,用插值获取当前时刻的电流 I_now = np.interp(t, t_profile, I_profile) ds_dt = -I_now / Qn / 3600.0 # 转化为每秒变化量 du1_dt = I_now / C1 - u1 / (R1 * C1) return [ds_dt, du1_dt] def battery_output(y, I_now, R0, ocv_coeff): s, u1 = y ocv = np.polyval(ocv_coeff, s) return ocv - I_now * R0 - u1

调用时把 OCV 多项式系数传进去,初始状态设为满电 (s=1, u1=0):

# 使用前面估计出的 tau,并假设 C1 = 1000 F(此时 R1 = tau / C1) C1_assume = 1000.0 R1_from_tau = tau_hat / C1_assume # 这里先用线性 OCV 近似折算成三阶多项式 ocv_coeff = np.poly1d([0.0, 0.0, k_hat, V0_hat]) current_profile = np.full_like(t_data, I) sol = solve_ivp( battery_ode, [t_data[0], t_data[-1]], [1.0, 0.0], t_eval=t_data, args=(current_profile, t_data, Qn, R0_hat, R1_from_tau, C1_assume, ocv_coeff), method="RK45", rtol=1e-6, atol=1e-8, ) V_ode = np.array([battery_output(y, I, R0_hat, ocv_coeff) for y in sol.y.T])

这里的C1_assume只是一个示范值,因为前面说过单恒流放电无法独立辨识 R1 和 C1。ODE 求解代码存在的意义,是给后续多工况扩展留好接口。我把这一块写在论文里,能直接回应“连续时间建模”的要求,比单纯甩一个curve_fit有说服力得多。

3.4 参数结果如何写进第一问答案

论文里不要只贴一个popt数组。要配一张参数表,并注明每个参数的物理意义和单位:

参数含义单位估计值
V0满电开路电压V4.198
k开路电压随 SoC 的变化斜率V-0.453
R0欧姆内阻Ω0.081
R1极化电阻Ω0.031
tauRC 时间常数s812.5

然后附上一句话结论:在 1.5A 恒流放电条件下,模型预测电压与实测数据的 RMSE 约为 x mV,说明一阶 RC 模型可以较准确描述该电池从满电到截止电压的动态过程。第一问到这里就算完整闭环了。

4. 拟合电池电压曲线最容易翻车的四件事

代码能跑通只是第一步,真正拉分的是你能不能处理拟合中的异常现象。下面四个坑我几乎每次都会遇到,按出现频率排序。

4.1 时间单位混用,让时间常数面目全非

这是最容易犯的错。题面如果按分钟给数据,你读取后直接用分钟代入秒单位的公式,指数项会变成 (e^{-t/800}) 和 (e^{-t/48000}) 的差异,时间常数估计直接偏移几十倍。我的处理习惯是进入拟合之前,统一用t_s = t_min * 60.0或t_s = t_hour * 3600.0完成转换,然后在 DataFrame 里保留原始单位列和标准秒单位列,方便后面写论文时对照。

还有 SoC 的变化率单位。很多队伍的公式写成 (ds/dt=-I/Q_n),但代码里 t 是秒,容量是毫安时,电流是安,三个单位互相打架,拟合出的 k 能偏出几倍。建议一律把容量转成安时、电流转成安、时间转成秒,公式里显式除以 3600。

4.2 初值选择不当,收敛到没有物理意义的解

curve_fit本质是局部最小二乘,对初值非常敏感。电池模型常见的情况是 R0 和 k 互相补偿:一个变大,一个变小,都能拟合出相近的电压曲线,但物理含义完全错误。我常用的初值策略是:

  • V0 从数据的第一点取,通常在 4.1 到 4.2;
  • k 初始设为 -0.5,因为一块电池从满电到放空,开路电压大约下降 0.5V 左右;
  • R0 初始设 0.05 到 0.1 欧姆,这个量级对手机电池接近真实;
  • R1 初始设 0.02 到 0.05 欧姆;
  • tau 初始设 500 到 1000 秒,锂电池的极化时间常数一般在分钟量级。

如果不确定,可以用多组初值并行拟合,选 RMSE 最小且参数都在物理合理范围内的解。

4.3 只看整体 RMSE,末端电压被淹没

电池放电数据里,高压段通常变化缓慢,低压段变化剧烈。如果直接用绝对电压误差做损失,拟合算法会把主要精力放在占样本比例更高的高压段,末端接近截止电压的行为反而被忽略。而竞赛题往往最关心“什么时候降到截止电压”,也就是末端行为。

我的做法分两步:先做一次全段拟合,然后看残差。如果末端残差明显偏大,就加大末段权重,例如给每点权重加一个随 SoC 减小的因子:

[ w_i = \frac{1}{\max(V_i - V_{cutoff}, 0.1)^2} ]

这样算法会更重视低压段。更严谨的做法是直接以“预测放电终止时间误差”为目标做拟合,因为第一问最终要回答的关键问题之一就是剩余可放电时间。

4.4 残差呈明显规律性时,不要急着加噪声项

拟合完如果画残差图,发现残差不是随机上下波动,而是呈现系统性的一撇或波浪形,这通常不是随机噪声,而是模型结构缺陷。常见的有两类:一是真实电池的开路电压曲线并非线性,二阶和三阶项被忽略,残差会随 SoC 呈现整体弯曲;二是电池存在二阶极化过程,一个 RC 支路不够,残差会呈现双时间常数特征。

遇到这种情况,不要急着在模型里加正则化或者用更高阶多项式硬卷。先判断物理上缺什么,再决定是否把 OCV 多项式升到二阶,或者把一阶 RC 扩成二阶 RC。哪怕是简单的 A/B 对比,也要说明“为什么增加这个结构能修正残差”,这才是建模题想要的推理过程。

5. 第一问交完之后,模型怎么接住后面几问

一道建模竞赛题通常不只一问。电池建模的后续方向几乎可以预判:预测剩余放电时间、估计循环寿命、评估温度对容量的影响、给出不确定性区间。第一问建立的连续时间模型,恰好能为这些扩展提供统一框架。

5.1 从电压模型到放电终止时间预测

第一问拟合出的 (V(t)) 解析式,可以直接解方程 (V(t)=V_{cutoff}) 求出剩余放电时间。

from scipy.optimize import brentq def voltage_minus_cutoff(t): return fit_func(t, *popt) - V_cutoff t_end = brentq(voltage_minus_cutoff, 0, t_data[-1]) print("预测放电终止时间(秒):", t_end)

这里必须注意brentq的搜索区间内要保证电压确实能降到截止值,否则会报错。如果模型预测的终止时间超出数据长度,说明模型外推能力有限,要在论文里诚实标注外推风险,不要硬给一个精确到小数点后两位的数字。

5.2 从“单次放电模型”到“循环寿命模型”

如果后续问电池经过多次充放电循环后容量如何衰减,那就要把 SoC 层面转到循环次数 n 或总放电容量层面。常见连续形式:

[ \frac{dQ}{dn} = -\alpha Q ]

解得 (Q(n)=Q_0 e^{-\alpha n})。更贴合实际的是带线性项的模型:

[ Q(n)=Q_0\left(1-\beta n\right)e^{-\alpha n} ]

前一段代表早期快速衰减,后一段代表长期线性老化。这种模型的好处仍然是连续函数,插值和预测都比离散表方便。第一问建立的参数估计流程可以原样复用:把“电压-时间”换成“容量-循环数”,把curve_fit的目标函数换成上面的指数形式即可。

5.3 温度因素与不确定性量化

真实电池的老化速度随温度明显变化,工程上常用阿伦尼乌斯公式修正衰减系数:

[ \alpha(T) = A \exp\left(-\frac{E_a}{R T}\right) ]

其中 (E_a) 是活化能,R 是气体常数,T 是绝对温度。如果题面提供了不同温度下的容量衰减数据,就能把阿伦尼乌斯项嵌入前面的一阶衰减模型。

如果题目要求给出预测区间,我会用两种办法:简单做法是假设参数服从拟合协方差矩阵给出的多元正态分布,对参数做几千次蒙特卡洛采样,得到电压或寿命预测的置信带;严格做法是用贝叶斯 MCMC,但竞赛时间有限,蒙特卡洛已经足够支撑一篇论文的讨论。

我的习惯是第一问先把单次放电的连续模型和处理流程吃透,参数、残差、RMSE、外推风险全部理清,后面几问不管怎么变,本质都落在这套“定义状态变量、写出状态方程、数据反解参数、验证模型结构”的循环里。题目可以换,电池可以换,但这套连续建模的思路不会变。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询