fs 3 histeresis electrical model

Data/FSLib/pyproject.toml

@@ -0,0 +1,2 @@
+[project]
+name = "FSLib"

Data/fs4voltagemodel.py

@@ -0,0 +1,136 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from dataclasses import dataclass
+
+# ======================================================
+# PARAMETERS
+# ======================================================
+Q_rated = 2.5 * 3600       # Coulombs (2.5 Ah cell)
+dt = 1.0                  # timestep [s]
+N = 300                   # number of steps
+
+# ======================================================
+# VEHICLE / BATTERY STATE
+# ======================================================
+@dataclass
+class VehicleState:
+    soc: float             # State of Charge [0..1]
+    voltage: float         # Terminal voltage [V]
+
+    # Reserved for future revised model
+    temp_c: float = 25.0
+    hyst: float = 0.0
+    v_rc: float = 0.0
+
+
+# ======================================================
+# VOLTAGE UPDATE TEMPLATE  ← THIS IS THE DELIVERABLE
+# ======================================================
+def update_voltage_template(prev_current: float, state: VehicleState) -> float:
+    """
+    Template for revised voltage updating (ECM baseline).
+
+    Model:
+      V = OCV(SOC) - I*R_internal
+
+    Inputs:
+      prev_current: current from previous timestep [A] (+ discharge, - charge)
+      state.soc: SOC in [0,1]
+
+    Output:
+      new terminal voltage [V]
+
+    TODO later:
+      - Replace OCV curve with real discharge curve fit
+      - Make R_internal depend on temperature/SOC
+      - Add hysteresis / RC recovery (sag + slow rebound)
+    """
+    soc = float(state.soc)
+
+    # 1) OCV curve (placeholder, but reasonable):
+    #    At SOC=1.0 -> ~4.2V, at SOC=0.0 -> ~3.0V
+    ocv = 3.0 + 1.2 * soc
+
+    # 2) Internal resistance (placeholder constant)
+    r_internal = 0.015  # Ohms (cell-level example)
+
+    # 3) Terminal voltage
+    v_new = ocv - prev_current * r_internal
+
+    # 4) Clamp to realistic bounds to avoid weird plots
+    v_new = float(np.clip(v_new, 2.5, 4.25))
+
+    return v_new
+
+
+# ======================================================
+# INPUT CURRENT PROFILE (synthetic)
+# ======================================================
+time = np.arange(N) * dt
+I = np.zeros(N)
+
+I[20:80]   = 5.0
+I[100:150] = 2.5
+I[170:220] = -3.0
+I[250:280] = 6.0
+
+# ======================================================
+# LOG ARRAYS
+# ======================================================
+SOC_log = np.zeros(N)
+V_log = np.zeros(N)
+
+# ======================================================
+# INITIAL STATE
+# ======================================================
+state = VehicleState(
+    soc=1.0,
+    voltage=4.2,   # reasonable initial guess
+)
+
+SOC_log[0] = state.soc
+V_log[0] = state.voltage
+
+# ======================================================
+# SIMULATION LOOP
+# ======================================================
+for t in range(1, N):
+    prev_I = I[t - 1]
+
+    # --- SOC UPDATE (Coulomb counting) ---
+    state.soc -= (prev_I * dt / Q_rated)
+    state.soc = float(np.clip(state.soc, 0.0, 1.0))
+
+    # --- VOLTAGE UPDATE (TEMPLATE CALL) ---
+    state.voltage = update_voltage_template(prev_I, state)
+
+    # --- LOG ---
+    SOC_log[t] = state.soc
+    V_log[t] = state.voltage
+
+# ======================================================
+# PLOTS
+# ======================================================
+plt.figure(figsize=(10, 6))
+
+plt.subplot(3, 1, 1)
+plt.plot(time, I, label="Current [A]")
+plt.ylabel("Current (A)")
+plt.grid(True)
+plt.legend()
+
+plt.subplot(3, 1, 2)
+plt.plot(time, V_log, label="Voltage [V]", color="tab:red")
+plt.ylabel("Voltage (V)")
+plt.grid(True)
+plt.legend()
+
+plt.subplot(3, 1, 3)
+plt.plot(time, SOC_log * 100, label="SOC [%]", color="tab:green")
+plt.xlabel("Time (s)")
+plt.ylabel("SOC (%)")
+plt.grid(True)
+plt.legend()
+
+plt.tight_layout()
+plt.show()

Data/pyproject.toml

@@ -0,0 +1,2 @@
+[project]
+name = "Data"

Docs/Modeling/ChemistryBatteryModel.md

@@ -0,0 +1,47 @@
+# Deriving a Battery model from first principles/chemistry
+by Nathaniel Platt
+
+In gen chem I learned about galvanic cells AKA your typical battery. They are made up of 2 half reactions: one for the anode and one for the cathode. You take the sum of their half reaction voltages to get $V_0$ or the voltage under equilibrium (same concentration of charge in the anode and cathode.) My guess is this is the same as the nominal voltage but I'm not sure. The whole equation is written out as $$V = V_0 + \frac{RT}{nF}*ln(\frac{[anode]}{[cathode]})$$ where $[anode]$ and $[cathode]$ are the concentration of ions in the anode/cathode respectively. 
+1. $R$ is the gas constant equal to ```8.31```. 
+1. $T$ is temperature in Kelvin ```293``` for 20 deg C.
+1. $n$ is the number of electrons processed in each reaction. I don't remember what this is at the time of creating this model so I assumed ```1``` which is a fine guess and will be fixed either way by an ML constant later. It will not show up in any more equations because it is assumed to be 1.
+1. $F$ is the faraday constant or ```96,485``` 
+
+If you make some basic assumptions you can expand this model further and allow a bit of machine learning.
+
+### Assumption 1: The volumes of each electrolyte are equivalent
+
+The concentration of ions in the anode and cathode are each just $\frac{Charge}{Volume}$ and if the volumes are the same then they cancel out:
+$$\frac{\frac{Q\substack{anode}}{\cancel{V\substack{anode}}}}{\frac{Q\substack{cathode}}{\cancel{V\substack{cathode}}}}$$
+- $Q\substack{anode}$ is the charge in the anode
+- $V\substack{anode}$ is the volume of the anode
+
+so we just get $\frac{Q_a}{Q_c}$ and since $Q_a + Q_c = Q\substack{total}$ or $6Ah$ we can effectively say it is $\frac{SOC}{1-SOC}$ so with that we have 
+$$V = V_0 + \frac{RT}{F}*ln(\frac{SOC}{1-SOC})$$
+
+### Assumption 2: The cell never goes below about 0.001M in either direction
+
+This assumption is based on not going below ```2.5V``` or above ```4.2V``` which correspond to SOC close to 0 or 1. When SOC approaches 0 or 1, the value out of the natural log skyrockets in the positive or negative direction (as seen in discharge curves when they go outside of those voltage bounds). In practice batteries tend to start growing dendrites in their anode or cathode and eventally short circuit and catch on fire. To accomplish this I essentially scaled down SOC a little bit by subtracting $0.1^a$ where an $a$ of $3$ worked well. This then looks like 
+
+$$V = V_0 + \frac{RT}{F}*ln(\frac{SOC-0.1^3}{1-(SOC-0.1^3)})$$
+
+### Assumption 3: Some of the values I came up with are wrong
+
+This is where the machine learning comes in. I threw in a few constants in areas where I knew I would be wrong that the model could correct for. It comes out to 
+
+$$V = V_0*C_4 + C_2\frac{RT}{F}*ln(\frac{C_1(SOC-0.1^3) + C_3}{1-(SOC-0.1^3)})$$
+
+1. $C_1$ accounts for my first assumption. It is likely that the anode and cathode volume are not the same and this allows that to be a parameter.
+1. $C_2$ accounts for any error in temperature dependence, to best fit the discharge curve slope, and for any errors in by assumption abotu $n$.
+1. $C_3$ is part of the $0.001M$ correction and helps shift it in the correct direction.
+1. $C_4$ accounts for errors in my assumption that the nominal voltage is equal to $V_0$ and generally just allows the model to fit that value rather than take it as an input.
+
+### Fix 1:
+
+So turns out I forgot ln rules and the $C_4$ and $C_1$ are redundant so removing $C_4$...
+
+### Fix 2:
+
+Adding $C_4$ to the bottom to do something similar to $C_3$. Also adjusting it so $C_1$ multiplies $C_3$ which makes it less separable but less confusing.
+
+$$V = V_0 + C_2\frac{RT}{F}*ln(\frac{C_1(SOC-0.1^3 + C_3)}{1-(SOC-0.1^3) + C_4})$$

FullVehicleSim/.vscode/settings.json

@@ -0,0 +1,4 @@
+{
+    "python-envs.defaultEnvManager": "ms-python.python:system",
+    "python-envs.pythonProjects": []
+}

FullVehicleSim/Electrical/HisteresisCellModel/F1_project1.py

@@ -0,0 +1 @@
+prin

FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining.py

@@ -0,0 +1,77 @@
+## When we train from data, use this
+
+from scipy.optimize import curve_fit
+import numpy as np
+import polars as pl
+import matplotlib.pyplot as plt
+from Data.FSLib.IntegralsAndDerivatives import integrate_with_Scipy_tCol
+from Data.FSLib.AnalysisFunctions import simpleTimeCol
+
+# df = pl.read_csv("C:/Projects/FormulaSlug/fs-data/FS-3/voltageTableVTC5A.csv")
+
+temps = []
+
+for i in range(5):
+    for j in range(6):
+        temps.append(f"ACC_SEG{i}_TEMPS_CELL{j}")
+
+df = pl.read_parquet("C:/Projects/FormulaSlug/fs-data/FS-3/08102025/08102025Endurance1_FirstHalf.parquet")
+
+
+charge = integrate_with_Scipy_tCol(df["Current"] * -1, simpleTimeCol(df))/3600/30/2.6 # Coulombs --> Ah, 30 cells --> 1 cell, 2.6Ah per cell
+
+df = df.with_columns(
+    pl.col("ACC_POWER_PACK_VOLTAGE").alias("Voltage"),
+    df.select(temps).mean_horizontal().alias("Temperature"),
+    pl.col("SME_TEMP_BusCurrent").alias("Current")
+)
+
+
+# plt.scatter(df["Charge"], df["Voltage"],c=df["Current"], label="Current")
+# plt.xlabel("Charge (Ah)")
+# plt.legend()
+# plt.show()
+
+dt = 0.01
+kernel_duration = 10.0
+kernel_size = int(kernel_duration / dt)
+t = np.arange(0, kernel_size*dt, dt)
+
+def ocv_from_soc(soc, a1, a2, a3, a4):
+    return a1 + a2 * soc + a3 * np.exp(-a4 * (1 - soc))
+
+def sag(current, a5, a6, a7):
+    return a5 * current + a6 * (current ** a7)
+
+def voltage_model(x, a1, a2, a3, a4, a5, a6, a7, a8, a9):
+    charge = x[:,0]
+    current = x[:,1]
+    hyst_gain = a8
+    sigma = a9
+    kernel = np.exp(-(t**2) / (2 * sigma**2))
+    kernel /= np.sum(kernel)
+    prev_curr = np.zeros((charge.shape[0], kernel_size))
+    for i in range(charge.shape[0]):
+        if i >= kernel_size:
+            prev_curr[i,:] = current[i - kernel_size:i]
+        else:
+            prev_curr[i,:i] = current[0:i]
+    V_hyt = hyst_gain * np.sum(prev_curr * t, axis=1)
+    V_ocv = ocv_from_soc(charge / 2.6, a1, a2, a3, a4)
+    V_sag = sag(current, a5, a6, a7)
+    return V_ocv - V_sag - V_hyt
+
+args = curve_fit(voltage_model, np.column_stack((df["Charge"], df["Current"])), df["Voltage"], p0=[3.0, 0.9, 0.25, 12.0, 0.02, 0.004, 1.3, 0.015, 3.0], maxfev=10000)
+args[0]
+
+a1, a2, a3, a4, a5, a6, a7, a8, a9 = args[0]
+plt.figure(figsize=(10,6))
+plt.scatter(df["Charge"], df["Voltage"], c='blue', label="Measured Voltage", alpha=0.5)
+predicted_voltage = voltage_model(np.column_stack((df["Charge"], df["Current"])), a1, a2, a3, a4, a5, a6, a7, a8, a9)
+plt.scatter(df["Charge"], predicted_voltage, c='red', label="Fitted Voltage", alpha=0.5)
+plt.xlabel("Charge (Ah)")
+plt.ylabel("Voltage (V)")
+plt.legend()
+plt.show()
+
+

FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py

@@ -0,0 +1,112 @@
+## When we train from data, use this
+
+from scipy.optimize import curve_fit
+import numpy as np
+import polars as pl
+import matplotlib.pyplot as plt
+from Data.FSLib.IntegralsAndDerivatives import integrate_with_Scipy_tCol
+from Data.FSLib.AnalysisFunctions import simpleTimeCol
+
+df = pl.read_csv("C:/Projects/FormulaSlug/fs-data/FS-3/voltageTableVTC5A.csv")
+dfLowCurr = df.filter(pl.col("Current") < 3).filter(pl.col("Voltage") > 2.5)
+
+df.head
+
+R = 8.314
+T = 293
+F = 96485.3329
+n = 1
+a = 3
+cmc = 2.6 # Cell Max Capacity (Ah)
+V0 = 3.7
+
+# temps = []
+
+# for i in range(5):
+#     for j in range(6):
+#         temps.append(f"ACC_SEG{i}_TEMPS_CELL{j}")
+
+# df = pl.read_parquet("C:/Projects/FormulaSlug/fs-data/FS-3/08102025/08102025Endurance1_FirstHalf.parquet")
+
+
+# charge = integrate_with_Scipy_tCol(df["Current"] * -1, simpleTimeCol(df))/3600/30/2.6 # Coulombs --> Ah, 30 cells --> 1 cell, 2.6Ah per cell
+
+# df = df.with_columns(
+#     pl.col("ACC_POWER_PACK_VOLTAGE").alias("Voltage"),
+#     df.select(temps).mean_horizontal().alias("Temperature"),
+#     pl.col("SME_TEMP_BusCurrent").alias("Current")
+# )
+
+
+# plt.scatter(df["Charge"], df["Voltage"],c=df["Current"], label="Current")
+# plt.xlabel("Charge (Ah)")
+# plt.legend()
+# plt.show()
+
+dt = 0.01
+kernel_duration = 10.0
+kernel_size = int(kernel_duration / dt)
+t = np.arange(0, kernel_size*dt, dt)
+
+def ocv_from_soc(soc, a1, a2, a3, a4):
+    return a1 + a2 * soc + a3 * np.exp(-a4 * (1 - soc))
+
+def sag(current, a5, a6, a7):
+    return a5 * current + a6 * (current ** a7)
+
+def voltage_SOC_model(x, c1, c2, c3, c4):
+    return V0 - c2*R*T/(n*F)*np.log((c1*((x/cmc) - (0.1**a) + c3))/(1 - (x/cmc) + (0.1**a) + c4))
+
+def voltage_model(x, a1, a2, a3, a4, a5, a6, a7, a8, a9):
+    charge = x[:,0]
+    current = x[:,1]
+    hyst_gain = a8
+    sigma = a9
+    kernel = np.exp(-(t**2) / (2 * sigma**2))
+    kernel /= np.sum(kernel)
+    prev_curr = np.zeros((charge.shape[0], kernel_size))
+    for i in range(charge.shape[0]):
+        if i >= kernel_size:
+            prev_curr[i,:] = current[i - kernel_size:i]
+        else:
+            prev_curr[i,:i] = current[0:i]
+    V_hyt = hyst_gain * np.sum(prev_curr * t, axis=1)
+    V_ocv = ocv_from_soc(charge / cmc, a1, a2, a3, a4)
+    V_sag = sag(current, a5, a6, a7)
+    return V_ocv - V_sag - V_hyt
+
+args = curve_fit(voltage_model, np.column_stack((df["Charge"], df["Current"])), df["Voltage"], p0=[3.0, 0.9, 0.25, 12.0, 0.02, 0.004, 1.3, 0.015, 3.0], maxfev=10000)
+args[0]
+
+a1, a2, a3, a4, a5, a6, a7, a8, a9 = args[0]
+plt.figure(figsize=(10,6))
+plt.scatter(df["Charge"], df["Voltage"], c='blue', label="Measured Voltage", alpha=0.5)
+predicted_voltage = voltage_model(np.column_stack((df["Charge"], df["Current"])), a1, a2, a3, a4, a5, a6, a7, a8, a9)
+plt.scatter(df["Charge"], predicted_voltage, c='red', label="Fitted Voltage", alpha=0.5)
+plt.xlabel("Charge (Ah)")
+plt.ylabel("Voltage (V)")
+plt.legend()
+plt.show()
+
+dfLowCurr1 = df.filter(pl.col("Current") < 3).filter(pl.col("Voltage") > 2.5)
+# dfLowCurr2 = df.filter(pl.col("Current") < 3)
+
+
+args1 = curve_fit(voltage_SOC_model, dfLowCurr1["Charge"], dfLowCurr1["Voltage"], p0=[1.5, 5.36, 0.0019, 3.7])
+# args2 = curve_fit(voltage_SOC_model, dfLowCurr2["Charge"], dfLowCurr2["Voltage"], p0=[1.5, 5.36, 0.0019, 3.7])
+
+c1, c2, c3, c4 = args1[0]
+# c5, c6, c7, c8 = args2[0]
+plt.figure(figsize=(10,6))
+plt.scatter(dfLowCurr1["Charge"], dfLowCurr1["Voltage"], c='blue', label="Measured Voltage", alpha=0.5)
+predicted_voltage_soc1 = voltage_SOC_model(np.arange(0, 2.6, 0.01), c1, c2, c3, c4)
+# predicted_voltage_soc2 = voltage_SOC_model(np.arange(0, 2.6, 0.01), c5, c6, c7, c8)
+plt.scatter(np.arange(0, 2.6, 0.01), predicted_voltage_soc1, c='red', label="Fitted Voltage 1", alpha=0.5)
+# plt.scatter(np.arange(0, 2.6, 0.01), predicted_voltage_soc2, c='green', label="Fitted Voltage 2", alpha=0.5)
+plt.xlabel("Charge (Ah)")
+plt.ylabel("Voltage (V)")
+plt.legend()
+plt.show()
+
+args1[0]
+# args2[0]