I have a battery simulation script where the battery’s state-of-charge (SOC) is updated iteratively. The charge and discharge values (how much energy is added/taken from the battery) at each timestep depend amongst other things, on the SOC at the previous timestep.

# Parameters
E_max = 100.0              # Maximum battery capacity
P_B_max = 10.0             # Maximum power for charge/discharge
SOC_min = 20.0             # Minimum state of charge
total_timesteps = 10       # Number of time steps

surplus = np.array([5, 3, 0, 7, 8, 0, 0, 5, 4, 3])# Excess production available for charging
remaining = np.array([0, 2, 5, 1, 0, 4, 6, 0, 2, 5])# Unmet load that requires discharge

SOC = np.zeros(total_timesteps + 1)
SOC[0] = 0.5 * E_max

charge = np.zeros(total_timesteps)
discharge = np.zeros(total_timesteps)

for t in range(total_timesteps):
    current_SOC = SOC[t]
    charge[t] = min(surplus[t], P_B_max, (E_max - current_SOC))
    discharge[t] = min(remaining[t], P_B_max, (current_SOC - SOC_min))
    SOC[t + 1] = np.clip(current_SOC + charge[t] - discharge[t],SOC_min, E_max)

I want to vectorize this calculation to remove the explicit for loop over t. However, I ve found it challenging because:

Charge and discharge values depend on the previous SOC. The amount of charge is constrained by how much available capacity remains in the battery. The amount of discharge is constrained by how much energy is available above the minimum SOC.

SOC updates sequentially, meaning that the SOC at time t+1 depends on SOC at time t, which was influenced by previous charge and discharge decisions.

I’m looking for an efficient way to vectorize or accelerate this calculation, ideally avoiding an explicit loop over timesteps altogether, while maintaining the correct behavior. Any guidance would be much appreciated appreciated!

I have a battery simulation script where the battery’s state-of-charge (SOC) is updated iteratively. The charge and discharge values (how much energy is added/taken from the battery) at each timestep depend amongst other things, on the SOC at the previous timestep.

# Parameters
E_max = 100.0              # Maximum battery capacity
P_B_max = 10.0             # Maximum power for charge/discharge
SOC_min = 20.0             # Minimum state of charge
total_timesteps = 10       # Number of time steps

surplus = np.array([5, 3, 0, 7, 8, 0, 0, 5, 4, 3])# Excess production available for charging
remaining = np.array([0, 2, 5, 1, 0, 4, 6, 0, 2, 5])# Unmet load that requires discharge

SOC = np.zeros(total_timesteps + 1)
SOC[0] = 0.5 * E_max

charge = np.zeros(total_timesteps)
discharge = np.zeros(total_timesteps)

for t in range(total_timesteps):
    current_SOC = SOC[t]
    charge[t] = min(surplus[t], P_B_max, (E_max - current_SOC))
    discharge[t] = min(remaining[t], P_B_max, (current_SOC - SOC_min))
    SOC[t + 1] = np.clip(current_SOC + charge[t] - discharge[t],SOC_min, E_max)

I want to vectorize this calculation to remove the explicit for loop over t. However, I ve found it challenging because:

Charge and discharge values depend on the previous SOC. The amount of charge is constrained by how much available capacity remains in the battery. The amount of discharge is constrained by how much energy is available above the minimum SOC.

SOC updates sequentially, meaning that the SOC at time t+1 depends on SOC at time t, which was influenced by previous charge and discharge decisions.

I’m looking for an efficient way to vectorize or accelerate this calculation, ideally avoiding an explicit loop over timesteps altogether, while maintaining the correct behavior. Any guidance would be much appreciated appreciated!

• python
• numpy
• vectorization

• While there might be a convoluted way of vectorising that with Numpy, I advise you to use modules like Numba or Cython for such codes (having loop-carried dependencies). The code will be more readable (similar to your current loop) and faster (MUCH faster than your current loop). – Jérôme Richard Commented Mar 4 at 16:39
• 'Vectorize' in numpy means working with whole arrays (in your case all time values), so compiled numpy methods can be used. With few exceptions like cumsum those methods are not sequential. It's hard to write general purpose building blocks that do custom things at each 'time step'. – hpaulj Commented Mar 4 at 17:05
• Is the simulation taking into account battery physics? Except for the limits (100, 20) It's not clear why charge/discharge depends on SOC. Is the script simulating a battery charger instead? – LMC Commented Mar 4 at 17:27

1 Answer 1

This is a time evolution, where the change depends on the current state, so can't be parallelised in time. Luckily initial value problems like this are common and python libraries exist which can solve them very efficiently and more accurately than the timstepping system.

Scipy provides an initial value problem solver solve_ivp which can use a range of different methods to solve.

Here is an example making use of solve_ivp to find the SOC at a range of times and plotting the resultant SOC using matplotlib.

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
import math

# Parameters
E_max = 100.0              # Maximum battery capacity
P_B_max = 10.0             # Maximum power for charge/discharge
SOC_min = 20.0             # Minimum state of charge
total_timesteps = 10       # Number of time steps

surplus = np.array([5, 3, 0, 7, 8, 0, 0, 5, 4, 3])# Excess production available for charging
remaining = np.array([0, 2, 5, 1, 0, 4, 6, 0, 2, 5])# Unmet load that requires discharge

def change_to_SOC(t, current_SOC) :
    charge = min(surplus[math.floor(t)], P_B_max, (E_max - current_SOC))
    discharge = min(remaining[math.floor(t)], P_B_max, (current_SOC - SOC_min))
    return charge-discharge

initial_SOC= 0.5 * E_max
t_eval = np.arange(0, 10, 0.01)
soc_plot = solve_ivp(change_to_SOC, [0,9.999], [initial_SOC], method = 'RK45', t_eval=t_eval)


plt.figure(figsize = (12, 4))
plt.plot(soc_plot.t, soc_plot.y[0])
plt.xlabel('t')
plt.ylabel('SOC(t)')
plt.show()

It plots

The result isn't exactly the same as yours because this is solving continuously in between steps using an Explicit Runge-Kutta method rather than just jumping at the integer timesteps.

Edit: I had to take a liberty by taking the floor of the time and using that to deduce the surplus and remaining, with knowledge of the source of these numbers, you may know better how to transition between times