28 fs 3 histeresis electrical model

Data/suspensionDamingView.py

@@ -0,0 +1,122 @@
+import polars as pl
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.fft as fft
+from statsmodels.nonparametric.smoothers_lowess import lowess
+
+df = pl.read_parquet("../fs-data/FS-3/03162026/2_steeper_regen_curve.parquet").fill_null(strategy="forward").fill_null(strategy="backward")
+
+
+[x for x in df.columns if "VDM" in x]
+
+
+FR_base_value = df["TPERIPH_FR_DATA_SUSTRAVEL"][:100].mean()
+FL_base_value = df["TPERIPH_FL_DATA_SUSTRAVEL"][:100].mean()
+BR_base_value = df["TPERIPH_BR_DATA_SUSTRAVEL"][:100].mean()
+BL_base_value = df["TPERIPH_BL_DATA_SUSTRAVEL"][:100].mean()
+
+t = df["Time_ms"]/1000.0
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.plot(t, df["TPERIPH_FR_DATA_SUSTRAVEL"] - FR_base_value, label="TPERIPH_FR_DATA_SUSTRAVEL")
+ax.plot(t, df["TPERIPH_FL_DATA_SUSTRAVEL"] - FL_base_value, label="TPERIPH_FL_DATA_SUSTRAVEL")
+ax.plot(t, df["TPERIPH_BR_DATA_SUSTRAVEL"] - BR_base_value, label="TPERIPH_BR_DATA_SUSTRAVEL")
+ax.plot(t, df["TPERIPH_BL_DATA_SUSTRAVEL"] - BL_base_value, label="TPERIPH_BL_DATA_SUSTRAVEL")
+ax.plot(t, df["VDM_Z_AXIS_YAW_RATE"]/2, label="VDM_Z_AXIS_YAW_RATE")
+ax.plot(t, df["SME_TRQSPD_Speed"]/1000, label="SME_TRQSPD_Speed")
+ax.legend()
+ax.set_xlabel("Time [s]")
+ax.set_ylabel("TPERIPH Data")
+ax.set_title("Suspension Damping View")
+ax.grid(True)
+plt.show()
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.plot(t, -1*df["TMAIN_DATA_STEERING"], label="TMAIN_DATA_STEERING")
+ax.plot(t, df["VDM_Z_AXIS_YAW_RATE"], label="VDM_Z_AXIS_YAW_RATE")
+ax.legend()
+ax.set_xlabel("Time [s]")
+ax.set_ylabel("Steering and Yaw Rate")
+ax.set_title("Steering and Yaw Rate")
+ax.grid(True)
+plt.show()
+
+df = df.filter(pl.col("SME_TRQSPD_Speed") > 5000)
+
+suspension_FR_fft = fft.fft((df["TPERIPH_FR_DATA_SUSTRAVEL"] - FR_base_value).to_numpy())
+suspension_FL_fft = fft.fft((df["TPERIPH_FL_DATA_SUSTRAVEL"] - FL_base_value).to_numpy())
+suspension_BR_fft = fft.fft((df["TPERIPH_BR_DATA_SUSTRAVEL"] - BR_base_value).to_numpy())
+suspension_BL_fft = fft.fft((df["TPERIPH_BL_DATA_SUSTRAVEL"] - BL_base_value).to_numpy())
+suspension_FR_freq = fft.fftfreq(len(suspension_FR_fft), d=0.01)
+suspension_FL_freq = fft.fftfreq(len(suspension_FL_fft), d=0.01)
+suspension_BR_freq = fft.fftfreq(len(suspension_BR_fft), d=0.01)
+suspension_BL_freq = fft.fftfreq(len(suspension_BL_fft), d=0.01)
+
+fig = plt.figure()
+ax1 = fig.add_subplot(221)
+ax1.scatter(suspension_FR_freq, np.abs(suspension_FR_fft), label="TPERIPH_FR_DATA_SUSTRAVEL", s=0.5)
+ax1.set_title("FR Suspension FFT")
+ax1.set_xlabel("Frequency [Hz]")
+ax1.set_ylabel("Magnitude")
+ax1.set_yscale("log")
+ax1.grid(True)
+ax2 = fig.add_subplot(222)
+ax2.scatter(suspension_FL_freq, np.abs(suspension_FL_fft), label="TPERIPH_FL_DATA_SUSTRAVEL", s=0.5)
+ax2.set_title("FL Suspension FFT")
+ax2.set_xlabel("Frequency [Hz]")
+ax2.set_ylabel("Magnitude")
+ax2.set_yscale("log")
+ax2.grid(True)
+ax3 = fig.add_subplot(223)
+ax3.scatter(suspension_BR_freq, np.abs(suspension_BR_fft), label="TPERIPH_BR_DATA_SUSTRAVEL", s=0.5)
+ax3.set_title("BR Suspension FFT")
+ax3.set_xlabel("Frequency [Hz]")
+ax3.set_ylabel("Magnitude")
+ax3.set_yscale("log")
+ax3.grid(True)
+ax4 = fig.add_subplot(224)
+ax4.scatter(suspension_BL_freq, np.abs(suspension_BL_fft), label="TPERIPH_BL_DATA_SUSTRAVEL", s=0.5)
+ax4.set_title("BL Suspension FFT")
+ax4.set_xlabel("Frequency [Hz]")
+ax4.set_ylabel("Magnitude")
+ax4.set_yscale("log")
+ax4.grid(True)
+plt.tight_layout()
+plt.show()
+
+yawRate_fft = fft.fft(df["VDM_Z_AXIS_YAW_RATE"].to_numpy())
+yawRate_freq = fft.fftfreq(len(yawRate_fft), d=0.01)
+pitchRate_fft = fft.fft(df["VDM_Y_AXIS_YAW_RATE"].to_numpy())
+pitchRate_freq = fft.fftfreq(len(pitchRate_fft), d=0.01)
+rollRate_fft = fft.fft(df["VDM_X_AXIS_YAW_RATE"].to_numpy())
+rollRate_freq = fft.fftfreq(len(rollRate_fft), d=0.01)
+
+
+meanSegments = 101
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.scatter(yawRate_freq, np.convolve(np.ones(meanSegments)*1/meanSegments, np.abs(yawRate_fft), mode='same'), label="Yaw", s=0.5)
+ax.scatter(pitchRate_freq, np.convolve(np.ones(meanSegments)*1/meanSegments, np.abs(pitchRate_fft), mode='same'), label="Pitch", s=0.5)
+ax.scatter(rollRate_freq, np.convolve(np.ones(meanSegments)*1/meanSegments, np.abs(rollRate_fft), mode='same'), label="Roll", s=0.5)
+ax.legend()
+ax.set_title("Yaw Rate FFT")
+ax.set_xlabel("Frequency [Hz]")
+ax.set_ylabel("Magnitude")
+ax.set_yscale("log")
+ax.grid(True)
+plt.show()
+
+## heave based on z acceleration
+heave_fft = fft.fft(df["VDM_Z_AXIS_ACCELERATION"].to_numpy())
+heave_freq = fft.fftfreq(len(heave_fft), d=0.01)
+fig = plt.figure()
+ax = fig.add_subplot(111)
+ax.scatter(heave_freq, np.convolve(np.ones(meanSegments)*1/meanSegments, np.abs(heave_fft), mode='same'), label="Heave", s=0.5)
+ax.legend()
+ax.set_title("Heave FFT")
+ax.set_xlabel("Frequency [Hz]")
+ax.set_ylabel("Magnitude")
+ax.set_yscale("log")
+ax.grid(True)
+plt.show()

FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py

@@ -4,11 +4,11 @@
 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
+# 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)
+dfLowCurr = df.filter(pl.col("Current") < 1).filter(pl.col("Voltage") > 2.5)
 
 df.head
 
@@ -88,7 +88,7 @@ def voltage_model(x, a1, a2, a3, a4, a5, a6, a7, a8, a9):
 plt.legend()
 plt.show()
 
-dfLowCurr1 = df.filter(pl.col("Current") < 3).filter(pl.col("Voltage") > 2.5)
+dfLowCurr1 = df.filter(pl.col("Current") < 1).filter(pl.col("Voltage") > 2.5)
 # dfLowCurr2 = df.filter(pl.col("Current") < 3)
 
 

FullVehicleSim/Electrical/HisteresisCellModel/electrical_model10.py

@@ -0,0 +1,133 @@
+import polars as pl
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.io
+from scipy.integrate import cumulative_simpson
+from scipy.signal import convolve
+from scipy.interpolate import CubicSpline as cs
+from scipy.interpolate import LinearNDInterpolator as Lnd_interp
+from scipy.optimize import curve_fit
+from numba import njit
+
+# 1. Load Data
+matData = scipy.io.loadmat("../fs-data/FS-4/P30B Cell Data/Molicel_INR18650P30B_measurement.mat", squeeze_me=True, struct_as_record=False)
+measurements = matData["measurement"]
+
+DELTA = 0.1  # Stick to 0.1s to capture high-rate dynamic transients accurately
+KERNEL_DURATION = 30
+KERNEL_SIZE = int(KERNEL_DURATION / DELTA)
+Q_NOMINAL = 3.0  
+
+def func_SOC(I, t, starting_SOC):
+    integrated_ah = cumulative_simpson(y=I, x=t, initial=0) / 3600.0
+    soc = starting_SOC - (integrated_ah / Q_NOMINAL)  # Discharge decreases SOC
+    return np.clip(soc, 0.0, 1.0)
+
+lowCurrMeasurements = [m for m in measurements.fu.DCC if "-0.1" in m.name]
+lowVs = np.concatenate([sim.V for sim in lowCurrMeasurements])
+lowTs = np.concatenate([sim.T_surf[0, :] for sim in lowCurrMeasurements])
+lowSOC = np.concatenate([func_SOC(sim.I, sim.t, 1.0) for sim in lowCurrMeasurements])
+
+voltage_curve = Lnd_interp(np.array([lowVs, lowTs]).T, lowSOC, rescale=True)
+
+def SOC_lookup(V, T):
+    SOC = voltage_curve(V, T)
+    if np.isnan(SOC):
+        if V > 4.1: return 1.0
+        if T < -15.0: return voltage_curve(V, -15.0)
+        if T > 39.0: return voltage_curve(V, 39.0)
+    return np.clip(SOC, 0.0, 1.0)
+
+def buildDF(measurement_set):
+    frames = []
+    for sim in measurement_set:
+        cubic = cs(sim.t, np.array([sim.I, sim.V, sim.T_surf[0, :], sim.T_surf[1, :]]), axis=1)
+        t = np.arange(sim.t[0], sim.t[-1], DELTA)
+        I, V, T_surf1, T_surf2 = cubic(t)
+        starting_SOC = SOC_lookup(V[0], T_surf1[0])
+        SOC = func_SOC(I, t, starting_SOC)
+        frames.append(pl.DataFrame({
+            "I": I, "V": V, "t": t, 
+            "T_surf1": T_surf1, "T_surf2": T_surf2, "SOC": SOC,
+        }))
+    return pl.concat(frames)
+
+all_sims = np.concatenate([measurements.fu.DCC, measurements.fu.CHC, measurements.fu.DCP, measurements.fu.CHP])
+fullData = buildDF(all_sims)
+
+# OCV Arguments
+ocv_args = np.array([2.33586300e+00, 8.70615991e+00, -3.26789721e+01, 7.63836093e+01, -9.98902031e+01, 6.83036001e+01, -1.90310137e+01, 5.90707089e-04])
+
+def OCV_terminal_voltage(SOC, T, a0, a1, a2, a3, a4, a5, a6, K_T):
+    return (a0 + a1 * SOC + a2 * SOC**2 + a3 * SOC**3 + a4 * SOC**4 + a5 * SOC**5 + a6 * SOC**6) * (1 + K_T * (T - 25))
+
+@njit
+def func_V_X_fast_dynamic(tauX, RX_arr, I, delta):
+    VX_arr = np.zeros_like(I)
+    if tauX <= 1e-5: return VX_arr
+    alpha = np.exp(-delta / tauX)
+    for i in range(1, len(I)):
+        beta = RX_arr[i] * (1.0 - alpha)
+        VX_arr[i] = VX_arr[i - 1] * alpha + I[i] * beta
+    return VX_arr
+
+def func_V_h(SOC, M, I, kernel):
+    # Fix: Hysteresis tracks current DIRECTION (sign), not raw multi-ampere values
+    I_sign = np.sign(I)
+    padded_I = np.pad(I_sign, (KERNEL_SIZE - 1, 0), mode='constant')
+    convolution = convolve(padded_I, kernel[::-1], mode='valid')
+    return M * SOC * convolution
+
+def terminal_voltage_dynamic(SOC, T, I, 
+                             r0_d0, r0_d1, r0_d2,  
+                             r0_c0, r0_c1, r0_c2,  
+                             r1_0, r1_1,           
+                             k_r, C1, R2, C2, tau_H, M, ocv_params):
+
+    temp_factor = (1.0 + k_r * (T - 25.0))
+
+    R0_d_stream = (r0_d0 + r0_d1 * SOC + r0_d2 / (SOC + 0.01)) * temp_factor
+    R0_c_stream = (r0_c0 + r0_c1 * SOC + r0_c2 / (SOC + 0.01)) * temp_factor
+    R1_stream = (r1_0 + r1_1 * (1.0 - SOC)) * temp_factor
+
+    tau1 = np.mean(R1_stream) * C1  
+    tau2 = R2 * C2
+
+    exponent = -np.arange(0, KERNEL_SIZE) * DELTA / tau_H
+    kernel = np.exp(exponent) / np.sum(np.exp(exponent))
+
+    V_h = func_V_h(SOC, M, I, kernel)
+    V_X1 = func_V_X_fast_dynamic(tau1, R1_stream, I, DELTA)
+    V_X2 = func_V_X_fast_dynamic(tau2, np.full_like(I, R2), I, DELTA)
+
+    V_OCV = OCV_terminal_voltage(SOC, T, *ocv_params)
+    V_R0 = np.where(I >= 0, R0_d_stream * I, R0_c_stream * I)
+
+    # FIX: Subtract transient polarization voltages (V_X1, V_X2) from V_OCV
+    V_est = V_OCV - V_R0 - V_X1 - V_X2 + V_h
+
+    print(f"Current Loop Step MAE: {np.mean(np.abs(V_est - fullData['V'].to_numpy())) * 1000.0:.2f} mV")
+    return V_est
+
+def wrapper_dynamic(X, r0_d0, r0_d1, r0_d2, r0_c0, r0_c1, r0_c2, r1_0, r1_1, k_r, C1, R2, C2, tau_H, M):
+    return terminal_voltage_dynamic(X[0], X[1], X[2], r0_d0, r0_d1, r0_d2, r0_c0, r0_c1, r0_c2, r1_0, r1_1, k_r, C1, R2, C2, tau_H, M, ocv_args)
+
+# Strictly enforce physical boundaries to prevent optimization divergence
+# Layout: [r0_d0, r0_d1, r0_d2, r0_c0, r0_c1, r0_c2, r1_0, r1_1, k_r, C1, R2, C2, tau_H, M]
+low_b = [1e-4, -0.05, 1e-5, 1e-4, -0.05, 1e-5, 1e-4, 1e-4, -0.05, 10.0, 1e-4, 100.0, 0.5, 1e-4]
+upr_b = [0.1,   0.05,  0.01, 0.1,   0.05,  0.01, 0.05, 0.05,  0.0,   5000, 0.1,  50000, 70.0, 0.1]
+
+p0_dynamic = [0.012, 0.0, 0.0005, 0.010, 0.0, 0.0005, 0.004, 0.002, -0.01, 200, 0.008, 4000, 25, 0.015]
+
+print("Starting bounded parameter estimation...")
+args_dynamic, _ = curve_fit(
+    wrapper_dynamic,
+    (fullData["SOC"].to_numpy(), fullData["T_surf1"].to_numpy(), fullData["I"].to_numpy()),
+    fullData["V"].to_numpy(),
+    p0=p0_dynamic,
+    bounds=(low_b, upr_b),  # Crucial bounds restored
+    maxfev=5000
+)
+
+print("\nFinal Optimized Parameters Array:")
+print(args_dynamic)