Physics-based battery models are essential tools in understanding and predicting the performance of batteries under various operating conditions. The Single Particle Model (SPM) is one of the simplest forms, representing each electrode as a single spherical particle. This model captures the basic electrochemical processes and is particularly useful for its computational efficiency. However, it oversimplifies the complex internal structures of the electrodes. The Doyle-Fuller-Newman (DFN) model, on the other hand, offers a more comprehensive approach. It includes the detailed physics of charge and mass transport within the electrodes and the electrolyte, making it a more accurate but computationally intensive model. Lastly, Equivalent Circuit Models (ECMs) represent the battery using electrical circuit elements like resistors, capacitors, and voltage sources. ECMs are widely used for battery management systems due to their simplicity and ability to fit experimental data, but they lack the detailed physics captured by models like the DFN. Each of these models serves different purposes, ranging from detailed analysis to real-time battery management applications.
In this example, we show how to model ECMs in Collimator. For SPM and DFN models, we recommend the [.code]pybamm[.code] Python library.
Equivalent Circuit Model (ECM)
The ECM considered is a Thevenin model with one RC branch. Its schematic is shown in the Figure below:
The variables and parameters of the model have the following interpretation:
Q: Total capacity of the battery (typically specified in ampere-hours)
s: Battery state-of-charge (SoC)
v_0: Open circuit voltage of the battery
v_t: Terminal voltage of the battery
vC1: Voltage across the capacitor
i: Discharge current (current to the load)
Rs: Internal resistance of the battery
R1: R1-C1 branch models slow diffusion processes in the battery (Ω)
C1: R1-C1 branch models slow diffusion processes in the battery [F]
State-of-charge (SoC) is a variable used to represent the current level of energy stored in a battery relative to its total capacity. Expressed as a percentage, SoC is an indicator of how much charge is remaining in the battery; an SoC of 100% (s=1) means the battery is fully charged, and an SoC of 0% (s=0) indicates it is completely discharged. The total capacity Q of a battery, measured typically in ampere-hours (Ah), is the maximum amount of charge it can hold when fully charged. The SoC can be calculated by dividing the current charge q (the amount of charge present at any given time) by the total capacity Q.
Note that v0, Rs, R1, and C1 typically depend on the state-of-charge s, cell temperature (T), and C-rate (rate of battery discharge relative to nominal capacity of the battery). However, in this demo we ignore dependencies on T and C-rate, and only consider dependencies of these parameters on s.
ECM equations
Based on the above Figure, the system of ordinary differential equations governing battery dynamics are:
where the factor $3600$ in the denominator converts $Q$ from Ampere-hours to Ampere-seconds. The terminal voltage $v_t$ is given by an algebraic relation:
$$ v_t = v_0 - i R_s - v_{C_1}. $$
Model implementation in Collimator using Foundational Blocks
We can implement this battery model using foundational blocks available in the Collimator block library. A screenshot of this model is shown here:
We can also build this model in code using the PyCollimator libraries:
from collimator.library import *
def make_battery_from_primitives(
soc_0=1.0, # initial value of SoC
vc1_0=0.0, # Iniital value of voltage across capacitor
Q=100.0, # Total capacity of the battery
soc_points=jnp.linspace(
0, 1, 11
), # discrete SoC points at which v_0, Rs, R1, and C1 are specified
v0_points=5.0 * jnp.ones(11), # discrete values of v0 at soc_points
Rs_points=15e-03 * jnp.ones(11), # discrete values of Rs at soc_points
R1_points=25e-03 * jnp.ones(11), # discrete values of R1 at soc_points
C1_points=3e03 * jnp.ones(11), # discrete values of C1 at soc_points
name="battery",
):
builder = collimator.DiagramBuilder()
# Creat and IO port for current input
current = builder.add(IOPort(name="discharge_current"))
# Keep track of soc through a current integrator (solve ODE for soc)
dot_soc = builder.add(Gain(-1.0 / 3600.0 / Q))
soc = builder.add(Integrator(soc_0, name="soc"))
builder.connect(current.output_ports[0], dot_soc.input_ports[0])
builder.connect(dot_soc.output_ports[0], soc.input_ports[0])
# Lookup blocks to capture the dependency of v0, Rs, R1, and C1 on soc
v0 = builder.add(LookupTable1d(soc_points, v0_points, "linear", name="v0"))
Rs = builder.add(LookupTable1d(soc_points, Rs_points, "linear", name="Rs"))
R1 = builder.add(LookupTable1d(soc_points, R1_points, "linear", name="R1"))
C1 = builder.add(LookupTable1d(soc_points, C1_points, "linear", name="C1"))
# Connect the soc to lookup blocks
builder.connect(soc.output_ports[0], v0.input_ports[0])
builder.connect(soc.output_ports[0], Rs.input_ports[0])
builder.connect(soc.output_ports[0], R1.input_ports[0])
builder.connect(soc.output_ports[0], C1.input_ports[0])
# ODE for vc1
i_over_C1 = builder.add(Product(2, operators="*/", name="i_over_C1"))
vc1_over_tau1 = builder.add(Product(3, operators="*//", name="vc1_over_tau1"))
dot_vc1 = builder.add(Adder(2, operators="+-", name="dot_vc1"))
vc1 = builder.add(Integrator(vc1_0, name="vc1"))
builder.connect(current.output_ports[0], i_over_C1.input_ports[0])
builder.connect(C1.output_ports[0], i_over_C1.input_ports[1])
builder.connect(vc1.output_ports[0], vc1_over_tau1.input_ports[0]) # feedback
builder.connect(R1.output_ports[0], vc1_over_tau1.input_ports[1])
builder.connect(C1.output_ports[0], vc1_over_tau1.input_ports[2])
builder.connect(i_over_C1.output_ports[0], dot_vc1.input_ports[0])
builder.connect(vc1_over_tau1.output_ports[0], dot_vc1.input_ports[1])
builder.connect(dot_vc1.output_ports[0], vc1.input_ports[0])
# Algebraic equation to get v (terminal voltage)
i_Rs = builder.add(Product(2, operators="**", name="i_Rs"))
vt = builder.add(Adder(3, operators="+--", name="vt"))
builder.connect(current.output_ports[0], i_Rs.input_ports[0])
builder.connect(Rs.output_ports[0], i_Rs.input_ports[1])
builder.connect(v0.output_ports[0], vt.input_ports[0])
builder.connect(i_Rs.output_ports[0], vt.input_ports[1])
builder.connect(vc1.output_ports[0], vt.input_ports[2])
# Diagram-level inport/outport
builder.export_input(current.input_ports[0]) # current
builder.export_output(soc.output_ports[0]) # soc
builder.export_output(vt.output_ports[0]) # terminal voltage
return builder.build(name=name)
...and then we can simulate our model. Let's apply a pulse discharge current to approximate non-continuous usage. Note that our dependence of $v_0$, $R_s$, $R_1$, and $C_1$ on $s$ is flat (default battery parameters) and battery capacity is 100 Ah:
builder = collimator.DiagramBuilder()
battery = make_battery_from_primitives(
name="battery"
) # Create a battery with the default parameters
builder.add(battery) # add the battery to the system
discharge_current = builder.add(
Pulse(100.0, 6 / 16.0, 16 * 60.0, name="discharge_current")
) # Add a pulse block
# We have two blocks/LeafSystems: discharge_current and battery. We need to connect the output
# of the discharge_current to the input of the battery.
builder.connect(discharge_current.output_ports[0], battery.input_ports[0])
diagram = builder.build()
context = diagram.create_context() # Create default context
# Specify which signals to record in the simulation output
recorded_signals = {
"discharge_current": discharge_current.output_ports[0],
"soc": battery.output_ports[0],
"vt": battery.output_ports[1],
}
t_span = (0.0, 9600.0) # span over which to simulate the system
# simulate the combination of diagram and context over t_span
sol = collimator.simulate(diagram, context, t_span, recorded_signals=recorded_signals)
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 6))
ax1.plot(sol.time, sol.outputs["discharge_current"], "-r", label=r"$i$")
ax2.plot(sol.time, sol.outputs["soc"], "-b", label=r"$s$")
ax3.plot(sol.time, sol.outputs["vt"], "-g", label=r"$v_t$")
ax1.set_ylabel("i [A]")
ax2.set_ylabel("soc [-]")
ax3.set_ylabel("$v_t$ [V]")
for ax in (ax1, ax2, ax3):
ax.set_xlabel("time [s]")
ax.legend(loc="best")
fig.tight_layout()
plt.show()
Running this simulation results in the expected behavior where $SoC$ decreases linearly with current:
Check out the next installment where we will estimate optimal parameters based on synthetic data.