BESS Subsystem Simulation Model
This document specifies the equivalent mathematical models for the Battery Energy Storage System (BESS) subsystem in a human-readable format, aimed at beginners who want to understand the system and implement it in code.
(Note: For automated AI code generation targets like Raspberry Pi or IDEs, please refer to the 03_bess_simulation_model_ai.md version which contains formal LaTeX math blocks and a Python pseudocode simulation loop.)
Terms and Abbreviations
| Term / Abbreviation | Definition |
|---|---|
| BESS | Battery Energy Storage System |
| BMS | Battery Management System |
| PCS | Power Conversion System (Inverter/Rectifier) |
| PCC | Point of Common Coupling (Grid injection point) |
| ECM | Equivalent Circuit Model |
| LMFP | Lithium Manganese Iron Phosphate |
| SoC | State of Charge (ratio of remaining capacity, 0.0 to 1.0) |
| SoH | State of Health (ratio of remaining lifespan maximum capacity, 0.0 to 1.0) |
| SoR | State of Resistance (aging factor: internal resistance increase multiplier) |
| OCV | Open Circuit Voltage (Battery voltage at rest, no current flowing) |
| HPPC | Hybrid Pulse Power Characterization (test for measuring ECM parameters) |
| RC pair | Resistor-Capacitor circuit pair (models delayed voltage response) |
| 1C rate | Current equal to the capacity in Amp-hours (e.g., 100 Ah battery at 100 A) |
System Overview
The BESS simulation is composed of five interconnected blocks. Understanding how they connect is essential before implementing any equations.
Data Flow (left to right = power flow during discharge):
[EMS Setpoint (Preq, Qreq)]
|
v
[4. PCS Model] <--(Vdc_pcs from DC Line)
| --> Pout, Qout, Ib
v
[2. BESS DC Line] --> Vdc_pcs
| --> Ib
v
[1. Battery Array (ECM + Thermal)]
| --> Vb, SoC, Tb
v
[5. BMS]
| --> P_limit_charge, P_limit_discharge, Alarms
v
[back to PCS as Preq_limited]
In plain language:
1. The EMS sends a power request (Preq) to the PCS. Positive = charge, negative = discharge.
2. The PCS converts the AC grid power to DC current (Ib) and injects it into the battery through the DC line.
3. The DC Line is a simple cable resistance that causes a small voltage drop.
4. The Battery Array receives Ib and updates its terminal voltage (Vb), State of Charge (SoC), and temperature (Tb).
5. The BMS monitors the array and computes safe operating power limits, which are fed back to cap the EMS request.
Sign Convention (applies to all blocks):
- Ib > 0 → Charging (current flows into the battery)
- Ib < 0 → Discharging (current flows out of the battery)
- Preq > 0 → Charging power request
- Preq < 0 → Discharging power request
1. Battery Array Model
1.1 Second-Order LMFP Equivalent Circuit Model
The battery is modeled using a second-order Thevenin Equivalent Circuit Model (ECM), which is well-suited for Lithium Manganese Iron Phosphate (LMFP) cells. The model has three internal elements: - R0: An instantaneous Ohmic resistance (voltage drop appears immediately when current flows). - R1/C1 and R2/C2: Two RC pairs that model the delayed electrochemical reactions inside the cell.
Reference Data: Experimental HPPC data for LMFP/LFP cells has been included in the repository. - Dataset File: lmfp_2nd_order_ecm_parameters.csv - Reference: Derived from open-source LFP HPPC parameterization sets. It contains OCV and RC parameters (R0, R1, C1, R2, C2) at 11 SoC breakpoints (0%, 10%, ... 100%).
How to use the CSV (Linear Interpolation):
At each simulation step, use the current SoC value to interpolate all parameters from the CSV table. Find the two rows bounding the current SoC and apply:
Value = V_low + (SoC - SoC_low) / (SoC_high - SoC_low) * (V_high - V_low)
This must be done for: OCV, R0, R1, C1, R2, C2 — at every time step.
Variables & States (values that change during simulation): - SoC: State of Charge (State, 0.0 to 1.0). Initial value: 0.5 - Vc1: Voltage across the first RC pair (State, Volts). Initial value: 0.0 V - Vc2: Voltage across the second RC pair (State, Volts). Initial value: 0.0 V - Ib: Battery array terminal current (Input, Amps). Positive = charging, Negative = discharging. - Vb: Battery array terminal voltage (Output, Volts)
Parameters (lookup table values — update at each step from CSV): - Q_cap: Total absolute array capacity (Amp-seconds). Default: 360,000 As (100 Ah) - OCV(SoC): Array Open Circuit Voltage (Volts). Lookup from CSV using current SoC. - R0: Series Ohmic resistance (Ohms). Lookup from CSV. - R1, R2: Polarization resistances (Ohms). Lookup from CSV. - C1, C2: Polarization capacitances (Farads). Lookup from CSV. - I_max_charge: Maximum continuous charge current limit (Amps). Default: 100 A (1C rate) - I_max_discharge: Maximum continuous discharge current limit (Amps). Default: 100 A (1C rate)
State of Charge (SoC) Estimation — Coulomb Counting
The State of Charge represents the remaining usable energy capacity. It is calculated by accumulating (integrating) the current flowing in or out at each simulation step (Coulomb Counting).
Sign rule: Positive current (charging) increases SoC. Negative current (discharging) decreases SoC.
Discrete update formula — implement this at every time step k:
SoC[k] = SoC[k-1] + ((Ib[k] * delta_t) / Q_cap)
Where:
- delta_t = simulation time step duration in seconds (e.g., 0.01 s for 10 ms steps)
- Q_cap = battery capacity in Amp-seconds (e.g., 360,000 As)
After computing SoC[k], clamp it: SoC[k] = clamp(SoC[k], 0.0, 1.0) to prevent going outside physical bounds.
Voltage Polarization Equations — RC Pair Updates
The RC circuits capture the delayed electrochemical voltage response. Think of them as "slow" voltages that lag behind the immediate Ohmic drop. Both RC states must be updated at each step using the discrete Forward-Euler method.
Discrete update formulas (implement both at every time step k):
Vc1[k] = Vc1[k-1] + delta_t * (-(Vc1[k-1] / (R1 * C1)) + (Ib[k] / C1))
Vc2[k] = Vc2[k-1] + delta_t * (-(Vc2[k-1] / (R2 * C2)) + (Ib[k] / C2))
Note: Use R1, C1, R2, C2 values interpolated from the CSV at the current SoC.
Output Terminal Voltage
The terminal voltage is the sum of the resting voltage (OCV) plus the voltage contributions from all internal impedances. Because the sign of Ib is positive for charging, all voltage contributions add to the OCV.
Formula:
Vb = OCV(SoC) + Vc1 + Vc2 + (Ib * R0)
Physical check: During charging (Ib > 0), Vb > OCV (the charger must push harder than the resting voltage). During discharging (Ib < 0), Vb < OCV (the voltage drops under load). This is correct.
1.2 Simplified Thermal Model
A lumped-capacitance thermal model tracks the average temperature of the battery array. The battery heats up from Joule losses and cools down by transferring heat to the ambient environment.
Variables & States: - Tb: Battery average temperature (State, degrees Celsius). Initial value: 25.0 C - Tamb: Ambient environment temperature (Input). Default: 25.0 C
Parameters: - m * cp: Total heat capacity of the array (Joules / degree Celsius). Default: 10,000,000 J/C - Rth: Convective thermal resistance to ambient (degrees Celsius / Watt). Default: 0.1 C/W
Heat Generation
Heat is generated by current flowing through all three internal resistances (Joule's Law):
Formula:
P_heat = (Ib^2 * R0) + (Vc1^2 / R1) + (Vc2^2 / R2)
Temperature Change — Discrete Update
Discrete update formula (implement at every time step k):
Tb[k] = Tb[k-1] + delta_t * (1 / (m * cp)) * (P_heat - ((Tb[k-1] - Tamb) / Rth))
Explanation: The temperature rises when heat generation (P_heat) exceeds the cooling rate ((Tb - Tamb) / Rth). The cooling rate increases as temperature rises above ambient — the hotter the battery, the faster it cools.
2. BESS DC Line
Model of the cable connecting the Battery Array to the PCS. Assumes minimum characteristics — purely resistive.
Parameters: - Rdc_line: Lumped resistance of the DC cabling. Default: 0.005 Ohms
Formula:
Vdc_pcs = Vb + (Ib * Rdc_line)
Direction note: During charging (Ib > 0), Vdc_pcs > Vb — the PCS pushes a higher voltage than the battery to force current in. During discharging (Ib < 0), Vdc_pcs < Vb — voltage drops across the line as current flows out.
3. BESS AC Line
Model of the cable between the PCS AC terminals and the grid Point of Common Coupling (PCC). Uses a simplified phasor/power-flow model, avoiding complex time-domain AC simulation.
Parameters:
- Vac_bus: Grid/Bus AC Voltage RMS magnitude. Default: 400 V (line-to-line)
- Rac_line, Lac_line: Line resistance and inductance. Defaults: Rac = 0.01 Ohms, Lac = 0.0001 H
- Xac_line: Line reactance = 2 * pi * 50 Hz * Lac_line. Default ≈ 0.03142 Ohms
- delta: Phase angle (radians) synthesized by the PCS relative to the grid bus. Controlled by PCS, typically small.
- Vac_pcs: Inverter synthesized AC Voltage RMS magnitude. Set by PCS control, close to Vac_bus.
Formulas (active and reactive power flow through the AC line):
Active Power (Pac_line) = (Vac_pcs * Vac_bus * sin(delta)) / Xac_line
Reactive Power (Qac_line) = (Vac_pcs^2 - (Vac_pcs * Vac_bus * cos(delta))) / Xac_line
Note: In this simplified model, delta is the key control variable. The PCS adjusts delta to achieve the desired Pac_line.
4. BESS PCS (Power Conversion System)
The PCS is modeled as a first-order lag (low-pass filter) on the power setpoint. This avoids modeling high-frequency inverter switching and only captures the slow control loop response.
Variables & States: - Preq, Qreq: Power setpoints from EMS (Watts, VAr). Positive = Charge, Negative = Discharge. - Pout, Qout: Actual instantaneous power the PCS delivers (States). Initial value: 0.0
Parameters: - tau_p, tau_q: Time constants (seconds). Control loop response speed. Default: 0.05 s (50 ms) - efficiency_pcs: Round-trip power conversion efficiency. Default: 0.98 - S_max: Maximum apparent power limit (VA). Default: 1,000,000 VA (1 MVA)
Step 1: Limit the Input Request
Before applying the setpoint, clamp it to the safe BMS limits (from Section 5):
Preq_limited = clamp(Preq, -P_limit_discharge, +P_limit_charge)
Preq_limited = clamp(Preq_limited, -S_max, +S_max)
This is the "limiter" block. It ensures the EMS cannot request more power than the battery or PCS hardware can provide.
Step 2: Output Power Tracking — Discrete Update
The actual output power ramps toward the limited request with a time constant tau_p:
Discrete update formula (implement at every time step k):
Pout[k] = Pout[k-1] + delta_t * ((Preq_limited - Pout[k-1]) / tau_p)
Qout[k] = Qout[k-1] + delta_t * ((Qreq_limited - Qout[k-1]) / tau_q)
Step 3: DC/AC Power Balance
Convert the AC output power to the required DC current, accounting for efficiency losses.
Formulas:
If charging (positive Pout): Pdc_pcs = Pout * efficiency_pcs
If discharging (negative Pout): Pdc_pcs = Pout / efficiency_pcs
Ib = Pdc_pcs / Vdc_pcs
(Implementation Note: Vdc_pcs depends on Ib which creates an algebraic loop. Use the Vdc_pcs value from the previous time step k-1 to break the loop.)
5. BMS (Battery Management System) Signals
The BMS is a monitoring and protection block. It reads battery states and computes safe operating limits and alarms. It does not directly control the battery — it informs the PCS and EMS.
5.1 Power Limits
These limits are calculated every time step and fed back to the PCS limiter (Preq_limited).
Step 1 — Current limits (structural + electrochemical):
I_limit_charge = min(I_max_charge, (V_cell_max - OCV(SoC)) / R0)
I_limit_discharge = min(I_max_discharge, (OCV(SoC) - V_cell_min) / R0)
Explanation: The first term is the hardware wire rating. The second term is the Ohmic voltage margin — the maximum current before the terminal voltage would hit the cell voltage limit.
Step 2 — Convert to power limits:
P_limit_charge = I_limit_charge * Vb
P_limit_discharge = I_limit_discharge * Vb
Note: Both values are positive magnitudes. The PCS clamps using -P_limit_discharge and +P_limit_charge.
5.2 State of Health (SoH) and State of Resistance (SoR) Update
SoH decreases slowly over the battery's life as cycles accumulate. SoR tracks the corresponding resistance increase.
Cycle counting — update once per full charge/discharge cycle:
SoH = SoH - SoH_degradation_per_cycle
SoR = 1.0 + (1.0 - SoH) * resistance_aging_factor
Parameters: - SoH_degradation_per_cycle: Capacity fade per cycle. Typical LMFP: 0.0002 (0.02% per cycle) - resistance_aging_factor: How much R0 scales with SoH loss. Default: 0.5
Application: Multiply the CSV lookup R0 by SoR at each step: R0_effective = R0_csv * SoR
5.3 Alarms (Evaluated cyclically)
- IF
SoC<SoC_min_thresholdTHENAlarm_Min_SoC= TRUE (Default SoC_min: 0.1) - IF
SoC>SoC_max_thresholdTHENAlarm_Max_SoC= TRUE (Default SoC_max: 0.9) - IF (
Vb/N_series) <V_cell_minTHENAlarm_Min_CellVolt= TRUE (Default V_cell_min: 2.8 V) - IF (
Vb/N_series) >V_cell_maxTHENAlarm_Max_CellVolt= TRUE (Default V_cell_max: 4.0 V) - IF
Tb>T_maxTHENAlarm_Max_Temp= TRUE (Default T_max: 60 C)
6. Model Calibratable Parameters
6.1 Battery Scaling Matrix (Auto-Adjustment Logic)
The simulation must scale cell parameters to the full array level. Specify any two of the three parameters below — the third is auto-calculated:
- BESS DC Nominal Voltage (V_BESS_nom): Default: 1000 V
- Cell DC Nominal Voltage (V_cell_nom): Typical LMFP is ~3.8V. Default: 3.8 V
- Number of cells in series (N_series): Default: 263 (= 1000 / 3.8, rounded)
Initialization Auto-Adjustment Rules:
- Rule 1: If V_BESS_nom and N_series are given → V_cell_nom = V_BESS_nom / N_series
- Rule 2: If V_cell_nom and N_series are given → V_BESS_nom = V_cell_nom * N_series
- Rule 3: If V_cell_nom and V_BESS_nom are given → N_series = round(V_BESS_nom / V_cell_nom)
The CSV parameters are per-cell. Scale OCV and all voltages by N_series. Scale R0, R1, R2 by N_series. Divide C1, C2 by N_series. Scale Q_cap to array total.
6.2 Initial Battery States
| Parameter | Description | Default |
|---|---|---|
| Initial SoC | Starting charge level (0.0 to 1.0) | 0.5 |
| Initial SoH | Starting health (1.0 = new, 0.0 = end-of-life) | 1.0 |
| Initial SoR | Starting resistance multiplier (1.0 = nominal) | 1.0 |
| Initial Tb | Starting battery temperature (°C) | 25.0 |
| Initial Vc1, Vc2 | Starting RC pair voltages (V) | 0.0 |
6.3 Runtime Fault Insertion
Simulation flags for testing EMS behavior under degraded conditions:
- Fault_PCS_Trip: Forces Pout = 0, Qout = 0 (contactor open).
- Fault_BMS_SensorLoss: Freezes SoC or Tb to a static incorrect value.
- Fault_ShortCircuit: Drops Vac_bus to zero to simulate grid collapse (FRT testing).
7. PCS Calibratable Parameters
| Parameter | Description | Default |
|---|---|---|
| S_max | PCS apparent power limit (VA) | 1,000,000 VA |
| efficiency_pcs | Power conversion efficiency | 0.98 |
| tau_p | Active power response time constant (s) | 0.05 s |
| tau_q | Reactive power response time constant (s) | 0.05 s |
| Vdc_min | Min DC voltage trip threshold (V) | 850 V |
| Vdc_max | Max DC voltage trip threshold (V) | 1200 V |