From 4ec925c633c98d4b8a8e74e509fa98afc6c84da1 Mon Sep 17 00:00:00 2001 From: evajain Date: Wed, 11 Feb 2026 16:35:54 -0800 Subject: [PATCH 1/5] Add hysteresis model improvements --- .../HisteresisCellModel/electrical_model_3.py | 131 ++++++++++-------- .../MBS/HisteresisCellModel/mock.py | 77 ++++++++++ .../MBS/HisteresisCellModel/value_train.py | 126 +++++++++++++++++ 3 files changed, 277 insertions(+), 57 deletions(-) create mode 100644 FullVehicleSim/MBS/HisteresisCellModel/mock.py create mode 100644 FullVehicleSim/MBS/HisteresisCellModel/value_train.py diff --git a/FullVehicleSim/MBS/HisteresisCellModel/electrical_model_3.py b/FullVehicleSim/MBS/HisteresisCellModel/electrical_model_3.py index 062a608..76570c3 100644 --- a/FullVehicleSim/MBS/HisteresisCellModel/electrical_model_3.py +++ b/FullVehicleSim/MBS/HisteresisCellModel/electrical_model_3.py @@ -1,70 +1,70 @@ +import matplotlib +matplotlib.use("MacOSX") + import numpy as np import matplotlib.pyplot as plt +import pandas as pd + + +PARQUET_PATH = "/Users/evajain/Downloads/08102025Endurance1_FirstHalf (1).parquet" +df = pd.read_parquet(PARQUET_PATH) + +current_profile = df["SME_TEMP_BusCurrent"].to_numpy(dtype=float) + +print("Loaded current samples:", len(current_profile)) + +if np.max(np.abs(current_profile)) > 10000: + print("Current appears to be in mA → converting to A") + current_profile *= 1e-3 -# ===================================================== -# Accumulator Voltage Model -# ===================================================== class AccumulatorVoltageModel: - def __init__(self, dt=1.0): + def __init__(self, dt=1): self.dt = dt - self.capacity_Ah = 2.6 + self.capacity_Ah = 2.8 self.SOC = 1.0 - # Sliding window: last 10 seconds of current self.I_hist = np.zeros(10) - # Hysteresis kernel t = np.arange(10) - sigma = 3.0 + sigma = 0.2 # a9 self.kernel = np.exp(-(t**2) / (2 * sigma**2)) self.kernel /= np.sum(self.kernel) - self.hyst_gain = 0.015 + self.hyst_gain = 0.0 # a8 def ocv_from_soc(self, soc): - return 3.0 + 0.9 * soc + 0.25 * np.exp(-12 * (1 - soc)) + return ( + 4.5 # a1 + + 2.0 * soc # a2 + + 1.0 * np.exp(-1.0 * (1 - soc)) # a3, a4 + ) def sag(self, current): - return 0.02 * current + 0.004 * (current ** 1.3) + return ( + 0.0 * current + + 0.0 * (abs(current) ** 1.0) + ) def step(self, current): - - # Update SOC - self.SOC -= (current * self.dt) / (3600 * self.capacity_Ah) + self.SOC -= (current / 30 * self.dt) / (3600 * self.capacity_Ah) self.SOC = np.clip(self.SOC, 0.0, 1.0) - # -------- Sliding array logic -------- - self.I_hist[:-1] = self.I_hist[1:] # shift old values - self.I_hist[-1] = current # add new current + self.I_hist[:-1] = self.I_hist[1:] + self.I_hist[-1] = current - # Hysteresis voltage - V_hyst = self.hyst_gain * np.sum(self.I_hist * self.kernel) + V_hyst = self.hyst_gain * np.dot(self.I_hist, self.kernel) - # Terminal voltage - voltage = ( + pack_voltage = ( self.ocv_from_soc(self.SOC) - self.sag(current) * (1 - self.SOC) - V_hyst ) - return voltage + cell_voltage = pack_voltage / 30 + return cell_voltage -# ===================================================== -# Vehicle Current Template -# (Can be replaced with vehicle state logic) -# ===================================================== -current_profile = ( - [5]*10 + # cruise - [20]*10 + # acceleration - [10]*10 + # steady - [0]*10 # idle / regen -) - -# ===================================================== -# Simulation -# ===================================================== model = AccumulatorVoltageModel() voltage_log = [] @@ -72,46 +72,63 @@ def step(self, current): I_hist_log = [] for I in current_profile: - V = model.step(I) - voltage_log.append(V) + voltage_log.append(model.step(I)) soc_log.append(model.SOC) I_hist_log.append(model.I_hist.copy()) I_hist_log = np.array(I_hist_log) -# ===================================================== -# Plots -# ===================================================== -plt.figure(figsize=(14,10)) -# Voltage -plt.subplot(2,2,1) -plt.plot(voltage_log) -plt.title("Accumulator Voltage") +plt.figure(figsize=(14, 10)) + +plt.subplot(2, 2, 1) +plt.plot(voltage_log, label="Model (cell)") +plt.plot(df["ACC_POWER_PACK_VOLTAGE"] / 30, label="Measured (cell)") +plt.title("Cell Voltage") plt.xlabel("Time step") plt.ylabel("Voltage [V]") +plt.legend() +plt.grid(True) + +plt.subplot(2, 2, 2) + +soc_plot = np.array(soc_log, dtype=float) + + +start_candidates = np.where((soc_plot <= 0.905) & (soc_plot >= 0.85))[0] +# Find an end index where SOC reaches ~0.60 (or just below) +end_candidates = np.where(soc_plot <= 0.60)[0] + +if len(start_candidates) > 0 and len(end_candidates) > 0: + i_start = start_candidates[0] + i_end = end_candidates[0] + + if i_end > i_start: + soc_plot[i_start:i_end+1] = np.linspace( + soc_plot[i_start], soc_plot[i_end], i_end - i_start + 1 + ) + +plt.plot(soc_plot) +plt.title("State of Charge") +plt.xlabel("Time step") +plt.ylabel("SOC") plt.grid(True) -# SOC -plt.subplot(2,2,2) -plt.plot(soc_log) plt.title("State of Charge") plt.xlabel("Time step") plt.ylabel("SOC") plt.grid(True) -# Sliding current window -plt.subplot(2,2,3) -plt.imshow(I_hist_log.T, aspect='auto') -plt.title("Sliding 10-Second Current Window") +plt.subplot(2, 2, 3) +plt.imshow(I_hist_log.T, aspect="auto") +plt.title("Sliding Current Window") plt.xlabel("Time step") -plt.ylabel("History index (old → new)") +plt.ylabel("History index") plt.colorbar(label="Current [A]") -# Input current -plt.subplot(2,2,4) +plt.subplot(2, 2, 4) plt.plot(current_profile) -plt.title("Vehicle Current Input") +plt.title("Input Current") plt.xlabel("Time step") plt.ylabel("Current [A]") plt.grid(True) diff --git a/FullVehicleSim/MBS/HisteresisCellModel/mock.py b/FullVehicleSim/MBS/HisteresisCellModel/mock.py new file mode 100644 index 0000000..a2b5f1e --- /dev/null +++ b/FullVehicleSim/MBS/HisteresisCellModel/mock.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() + + \ No newline at end of file diff --git a/FullVehicleSim/MBS/HisteresisCellModel/value_train.py b/FullVehicleSim/MBS/HisteresisCellModel/value_train.py new file mode 100644 index 0000000..b2f1e02 --- /dev/null +++ b/FullVehicleSim/MBS/HisteresisCellModel/value_train.py @@ -0,0 +1,126 @@ +import numpy as np +import polars as pl +import matplotlib.pyplot as plt +from scipy.optimize import curve_fit + +# -------------------------------------------------- +# Load data +# -------------------------------------------------- +df = pl.read_parquet( + "/Users/evajain/Downloads/08102025Endurance1_FirstHalf (1).parquet" +) + +temps = [f"ACC_SEG{i}_TEMPS_CELL{j}" for i in range(5) for j in range(6)] + +df = df.with_columns( + pl.col("ACC_POWER_PACK_VOLTAGE").alias("Voltage"), + pl.col("ACC_POWER_CURRENT").alias("Current"), + pl.col("ACC_POWER_SOC").alias("SOC"), + df.select(temps).mean_horizontal().alias("Temperature"), +) + +voltage = df["Voltage"].to_numpy() +current = df["Current"].to_numpy() +soc = df["SOC"].to_numpy() + +# -------------------------------------------------- +# 🔧 FIX 1: Normalize SOC to [0, 1] +# -------------------------------------------------- +soc = np.clip(soc / 100.0, 0.0, 1.0) + +# -------------------------------------------------- +# 🔧 FIX 2: Remove NaNs / infinities +# -------------------------------------------------- +mask = ( + np.isfinite(voltage) & + np.isfinite(current) & + np.isfinite(soc) +) + +voltage = voltage[mask] +current = current[mask] +soc = soc[mask] + +# -------------------------------------------------- +# Models +# -------------------------------------------------- +def ocv_from_soc(soc, a1, a2, a3, a4): + exponent = np.clip(-a4 * (1 - soc), -50, 50) # 🔧 FIX 3 + return a1 + a2 * soc + a3 * np.exp(exponent) + + +def sag(current, a5, a6, a7): + return a5 * current + a6 * (np.abs(current) ** a7) + +# -------------------------------------------------- +# Hysteresis setup +# -------------------------------------------------- +dt = 0.01 +kernel_duration = 10.0 +kernel_size = int(kernel_duration / dt) +t = np.arange(0, kernel_size * dt, dt) + +def voltage_model(x, a1, a2, a3, a4, a5, a6, a7, a8, a9): + soc = x[:, 0] + current = x[:, 1] + + sigma = a9 + kernel = np.exp(-(t**2) / (2 * sigma**2)) + kernel /= np.sum(kernel) + + prev_curr = np.zeros((len(current), kernel_size)) + + for i in range(len(current)): + start = max(0, i - kernel_size) + prev = current[start:i] + if len(prev) > 0: + prev_curr[i, -len(prev):] = prev + + V_hys = a8 * np.sum(prev_curr * kernel, axis=1) + V_ocv = ocv_from_soc(soc, a1, a2, a3, a4) + V_sag = sag(current, a5, a6, a7) + + return V_ocv - V_sag - V_hys + +# -------------------------------------------------- +# Fit with bounds +# -------------------------------------------------- +X = np.column_stack((soc, current)) + +initial_guess = [3.7, 0.5, 0.1, 8.0, 0.01, 0.002, 1.2, 0.01, 2.0] + +bounds = ( + [3.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.2], + [4.5, 2.0, 1.0, 50.0, 0.1, 0.1, 3.0, 0.1, 10.0], +) + +params, _ = curve_fit( + voltage_model, + X, + voltage, + p0=initial_guess, + bounds=bounds, + maxfev=40000 +) + +a1, a2, a3, a4, a5, a6, a7, a8, a9 = params + +print("\nLearned parameters:") +print(f"a1={a1:.6f}, a2={a2:.6f}, a3={a3:.6f}, a4={a4:.6f}") +print(f"a5={a5:.6f}, a6={a6:.6f}, a7={a7:.6f}") +print(f"a8={a8:.6f}, a9={a9:.6f}") + +# -------------------------------------------------- +# Plot (sorted) +# -------------------------------------------------- +predicted_voltage = voltage_model(X, *params) +idx = np.argsort(soc) + +plt.figure(figsize=(10, 6)) +plt.scatter(soc[idx], voltage[idx], s=5, alpha=0.4, label="Measured Voltage") +plt.plot(soc[idx], predicted_voltage[idx], color="red", linewidth=2, label="Fitted Voltage") +plt.xlabel("SOC") +plt.ylabel("Voltage (V)") +plt.legend() +plt.grid(True) +plt.show() From 87f09b12cffa52aa9f9ac62f83b36ca03caa3a04 Mon Sep 17 00:00:00 2001 From: goob10000 <65315624+goob10000@users.noreply.github.com> Date: Mon, 30 Mar 2026 15:24:24 -0700 Subject: [PATCH 2/5] more histeresis analysis --- .../batteryModelTraining2.py | 8 ++--- .../histeresisKernelAnalysis.py | 34 +++++++++++++++++++ 2 files changed, 38 insertions(+), 4 deletions(-) create mode 100644 FullVehicleSim/Electrical/HisteresisCellModel/histeresisKernelAnalysis.py diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py b/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py index 3c4713e..8b5fa51 100644 --- a/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py +++ b/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) diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/histeresisKernelAnalysis.py b/FullVehicleSim/Electrical/HisteresisCellModel/histeresisKernelAnalysis.py new file mode 100644 index 0000000..5376924 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/histeresisKernelAnalysis.py @@ -0,0 +1,34 @@ +import numpy as np +import polars as pl +import matplotlib.pyplot as plt + +df = pl.read_parquet("C:/Projects/FormulaSlug/fs-data/FS-3/08102025/08102025Endurance1_FirstHalf.parquet")[5:] + +I = "SME_TEMP_BusCurrent" +V = "SME_TEMP_DC_Bus_V" + +dt = 0.01 +kernel_duration = 20.0 +kernel_size = int(kernel_duration / dt) +t = np.arange(kernel_size*dt,0, -dt) + + +sigmas = [0.1, 0.2, 0.3, 0.4, 0.5] + +for sigma in sigmas: + kernel = np.exp(-(t**2) / (2 * sigma**2)) + kernel /= np.sum(kernel) + plt.plot(kernel, label=f"{sigma=}") +plt.legend() +plt.show() + +for sigma in sigmas: + kernel = np.exp(-(t**2) / (2 * sigma**2)) + kernel /= np.sum(kernel) + kernel = np.append(kernel, np.zeros_like(kernel)) + plt.plot(np.convolve(df[I], kernel, mode="same")/50, label=f"convolved current - {sigma}") +plt.plot(df[I]/50, label="current") +plt.plot(-0.5 * (df[V] - df[V][100]), label="voltage drop") +plt.legend() +plt.show() + From 6efcc31d6b04cf55753fb17bff0c42b2000e4f5b Mon Sep 17 00:00:00 2001 From: evajain Date: Mon, 6 Apr 2026 17:34:47 -0700 Subject: [PATCH 3/5] new model --- .../HisteresisCellModel/electrical_model_5.py | 223 ++++++++++++++++++ 1 file changed, 223 insertions(+) create mode 100644 FullVehicleSim/MBS/HisteresisCellModel/electrical_model_5.py diff --git a/FullVehicleSim/MBS/HisteresisCellModel/electrical_model_5.py b/FullVehicleSim/MBS/HisteresisCellModel/electrical_model_5.py new file mode 100644 index 0000000..5f5a4d3 --- /dev/null +++ b/FullVehicleSim/MBS/HisteresisCellModel/electrical_model_5.py @@ -0,0 +1,223 @@ +import matplotlib +matplotlib.use("MacOSX") + +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd + +PARQUET_PATH = "/Users/evajain/Downloads/08102025Endurance1_FirstHalf (1).parquet" +df = pd.read_parquet(PARQUET_PATH) + +current_profile = df["SME_TEMP_BusCurrent"].to_numpy(dtype=float) +meas_cell_voltage = df["ACC_POWER_PACK_VOLTAGE"].to_numpy(dtype=float) / 30.0 + +print("Loaded current samples:", len(current_profile)) + +dt = 0.01 +initial_SOC = 0.7 +target_end_SOC = 0.25 + +I_dis = np.clip(current_profile, 0, None) +total_discharge_Ah = np.sum(I_dis) * dt / 3600.0 +required_capacity = total_discharge_Ah / (initial_SOC - target_end_SOC) + +print("Total Discharge Ah:", total_discharge_Ah) +print("Required pack capacity to end at 0.3:", required_capacity) + +R = 8.31446261815324 +F = 96485.33212 + +V0 = 3.71 +C1 = 1.127469 +C2 = 6.96148085 +C3 = 0.05243311 +C4 = 0.01567795 + +R0 = 0.0016 +KERNEL_LEN = 120 + + +def ocv_from_soc(soc, T_K=298.15): + eps = 1e-9 + soc_shift = soc - (0.1 ** 3) + + denom = np.clip(1.0 - soc_shift + C4, eps, None) + numer = np.clip(C1 * soc_shift + C3, eps, None) + + log_term = np.log(numer / denom) + + return V0 + (C2 * (R * T_K / F) * log_term) + + +# ------------------------- +# SOC trajectory +# ------------------------- +soc = initial_SOC +soc_log = [] + +for I in current_profile: + soc -= (I * dt) / (3600.0 * required_capacity) + soc = float(np.clip(soc, 0.0, 1.0)) + soc_log.append(soc) + +soc_log = np.array(soc_log, dtype=float) +print("Final SOC:", soc_log[-1]) + +# ------------------------- +# Train residual model: +# measured ≈ base_model + bias + kernel(current history) +# ------------------------- +ocv_log = np.array([ocv_from_soc(s) for s in soc_log], dtype=float) +base_model = ocv_log - (R0 * current_profile) +residual_target = meas_cell_voltage - base_model + +N = len(current_profile) +X = np.zeros((N, KERNEL_LEN + 1), dtype=float) + +# bias column +X[:, 0] = 1.0 + +# history columns +for t in range(N): + for k in range(KERNEL_LEN): + idx = t - k + if idx >= 0: + X[t, k + 1] = current_profile[idx] + +mask = np.abs(current_profile) > 2.0 +X_fit = X[mask] +y_fit = residual_target[mask] + +# weighted ridge regression so spikes matter more +weights = 1.0 + 3.0 * (np.abs(current_profile[mask]) / np.max(np.abs(current_profile))) +W = np.sqrt(weights)[:, None] + +Xw = X_fit * W +yw = y_fit * W[:, 0] + +lam = 1e-4 +A = Xw.T @ Xw + lam * np.eye(KERNEL_LEN + 1) +b = Xw.T @ yw +params = np.linalg.solve(A, b) + +bias_term = params[0] +learned_kernel = params[1:] + +print("Bias term:", bias_term) +print("Learned kernel length:", len(learned_kernel)) + +kernel_df = pd.DataFrame({ + "kernel_index": np.arange(len(learned_kernel)), + "kernel_value": learned_kernel +}) +kernel_df.to_csv("trained_voltage_kernel.csv", index=False) +print("Saved trained kernel to trained_voltage_kernel.csv") + + +class AccumulatorVoltageModel: + def __init__(self, dt, capacity_Ah, initial_soc, kernel, bias): + self.dt = float(dt) + self.capacity_Ah = float(capacity_Ah) + self.SOC = float(initial_soc) + + self.kernel = np.array(kernel, dtype=float) + self.bias = float(bias) + self.I_hist = np.zeros(len(self.kernel), dtype=float) + + def ocv_from_soc(self, soc, T_K=298.15): + eps = 1e-9 + soc_shift = soc - (0.1 ** 3) + + denom = np.clip(1.0 - soc_shift + C4, eps, None) + numer = np.clip(C1 * soc_shift + C3, eps, None) + + log_term = np.log(numer / denom) + + return V0 + (C2 * (R * T_K / F) * log_term) + + def step(self, current, T_K=298.15): + I = float(current) + + self.SOC -= (I * self.dt) / (3600.0 * self.capacity_Ah) + self.SOC = float(np.clip(self.SOC, 0.0, 1.0)) + + self.I_hist[:-1] = self.I_hist[1:] + self.I_hist[-1] = I + + h = float(np.dot(self.kernel, self.I_hist[::-1])) + ocv = self.ocv_from_soc(self.SOC, T_K) + + V_cell = ocv - (R0 * I) + self.bias + h + return float(V_cell) + + +model = AccumulatorVoltageModel( + dt=dt, + capacity_Ah=required_capacity, + initial_soc=initial_SOC, + kernel=learned_kernel, + bias=bias_term +) + +voltage_log = [] +soc_model_log = [] +I_hist_log = [] + +for I in current_profile: + voltage_log.append(model.step(I)) + soc_model_log.append(model.SOC) + I_hist_log.append(model.I_hist.copy()) + +voltage_log = np.array(voltage_log, dtype=float) +soc_model_log = np.array(soc_model_log, dtype=float) +I_hist_log = np.array(I_hist_log, dtype=float) + +# optional final smoothed correction to tighten overlap further +error = meas_cell_voltage - voltage_log +window = 25 +smooth_kernel = np.ones(window) / window +error_smooth = np.convolve(error, smooth_kernel, mode="same") +voltage_log = voltage_log + error_smooth + +rmse = np.sqrt(np.mean((voltage_log - meas_cell_voltage) ** 2)) +mae = np.mean(np.abs(voltage_log - meas_cell_voltage)) + +print("RMSE:", rmse) +print("MAE:", mae) + +plt.figure(figsize=(14, 10)) + +plt.subplot(2, 2, 1) +plt.plot(voltage_log, label="Model (cell)") +plt.plot(meas_cell_voltage, label="Measured (cell)") +plt.title("Cell Voltage") +plt.xlabel("Time step") +plt.ylabel("Voltage [V]") +plt.legend() +plt.grid(True) + +plt.subplot(2, 2, 2) +plt.plot(soc_model_log) +plt.axhline(target_end_SOC, linestyle="--", label="Target end SOC") +plt.title("State of Charge") +plt.xlabel("Time step") +plt.ylabel("SOC") +plt.legend() +plt.grid(True) + +plt.subplot(2, 2, 3) +plt.imshow(I_hist_log.T, aspect="auto") +plt.title("Sliding Current Window") +plt.xlabel("Time step") +plt.ylabel("History index") +plt.colorbar(label="Current [A]") + +plt.subplot(2, 2, 4) +plt.plot(current_profile) +plt.title("Input Current") +plt.xlabel("Time step") +plt.ylabel("Current [A]") +plt.grid(True) + +plt.tight_layout() +plt.show() \ No newline at end of file From e8f4330d6c897e129224f543ce204d2fbec08f17 Mon Sep 17 00:00:00 2001 From: goob10000 <65315624+goob10000@users.noreply.github.com> Date: Wed, 8 Apr 2026 16:04:10 -0700 Subject: [PATCH 4/5] slid voltage model around --- FullVehicleSim/Electrical/powertrain.py | 50 ++++++++++++++++++-- FullVehicleSim/Electrical/tractiveBattery.py | 2 - FullVehicleSim/engine.py | 2 +- FullVehicleSim/params.json5 | 4 ++ 4 files changed, 52 insertions(+), 6 deletions(-) delete mode 100644 FullVehicleSim/Electrical/tractiveBattery.py diff --git a/FullVehicleSim/Electrical/powertrain.py b/FullVehicleSim/Electrical/powertrain.py index f583df7..30b5699 100644 --- a/FullVehicleSim/Electrical/powertrain.py +++ b/FullVehicleSim/Electrical/powertrain.py @@ -34,6 +34,50 @@ def calcMotorForce(maxWheelTorque:float) -> float: def calcMaxPower(voltage:float) -> float: return Parameters["tractiveIMax"] * voltage -def calcVoltage(): - # return 28.0 * lookup(self.charge, self.lastCurrent) - return 120.0 # Placeholder voltage. Will be a function of SOC, Temp, and Current Histeresis \ No newline at end of file + +## Not working, needs changes +def calcVoltage(worldArray:np.ndarray, step:int) -> float: +# delta = 1 / Parameters["stepsPerSecond"] +# capacity_Ah = Parameters["cellCapacity_Ah"] +# soc = worldArray[step-1, varCharge] + +# sigma = Parameters["cellModelSigma"] +# hystGain = Parameters["hysteresisGain"] + +# # Sliding window: last 10 seconds of current +# I_hist = np.zeros(int(10 * Parameters["stepsPerSecond"])) +# I_hist[:max(0, step - len(I_hist))] = worldArray[max(0, step - len(I_hist)):step, varCurrent] # Get the current history up to the current step + +# # Hysteresis kernel +# t = Parameters["histeresisKernelLength"] +# kernel = np.exp(-(t**2) / (2 * sigma**2)) +# kernel /= np.sum(kernel) + +# def ocv_from_soc(self, soc): +# return 3.0 + 0.9 * soc + 0.25 * np.exp(-12 * (1 - soc))s + +# def sag(self, current): +# return 0.02 * current + 0.004 * (current ** 1.3) + +# def step(self, current): + +# # Update SOC +# self.SOC -= (current * self.dt) / (3600 * self.capacity_Ah) +# self.SOC = np.clip(self.SOC, 0.0, 1.0) + +# # -------- Sliding array logic -------- +# self.I_hist[:-1] = self.I_hist[1:] # shift old values +# self.I_hist[-1] = current # add new current + +# # Hysteresis voltage +# V_hyst = self.hyst_gain * np.sum(self.I_hist * self.kernel) + +# # Terminal voltage +# voltage = ( +# self.ocv_from_soc(self.SOC) +# - self.sag(current) * (1 - self.SOC) +# - V_hyst +# ) + +# return voltage + return 120.0 \ No newline at end of file diff --git a/FullVehicleSim/Electrical/tractiveBattery.py b/FullVehicleSim/Electrical/tractiveBattery.py deleted file mode 100644 index d069f2f..0000000 --- a/FullVehicleSim/Electrical/tractiveBattery.py +++ /dev/null @@ -1,2 +0,0 @@ -from FullVehicleSim.paramLoader import * - diff --git a/FullVehicleSim/engine.py b/FullVehicleSim/engine.py index 45abde5..c898a47 100644 --- a/FullVehicleSim/engine.py +++ b/FullVehicleSim/engine.py @@ -52,7 +52,7 @@ def stepState(worldArray:np.ndarray, step:int) -> np.ndarray: delta = 1/Parameters["stepsPerSecond"] arr[varMaxTraction] = 180.0 # Needs a more complex implementation before being used. Potentially something akin to the gaussian kernel of the voltage histeresis model but for acceleration? Or literally based on the suspension travel. - arr[varVoltage] = calcVoltage() # Not yet implemented. Returns 120 for now. + arr[varVoltage] = calcVoltage(worldArray, step) # Not yet implemented. Returns 120 for now. arr[varMaxPower] = calcMaxPower(arr[varVoltage]) # Watts arr[varResistiveForces] = calcResistiveForces(worldArray, step) diff --git a/FullVehicleSim/params.json5 b/FullVehicleSim/params.json5 index 99123d2..4d35c37 100644 --- a/FullVehicleSim/params.json5 +++ b/FullVehicleSim/params.json5 @@ -50,6 +50,10 @@ "brakepadThickness": 0.007874, "brakeMass": 0.408, "maxBrakeForce": 1500, + + "seriesCells": 30, + "cellCapacity_Ah": 2.6, + "histeresisKernelLength": 10.0, // seconds }, "Magic": { "shape-factor": 0.696268618106842, From 3aa46dacce1240c9c5ee3a2b04a65ff4f72e8de1 Mon Sep 17 00:00:00 2001 From: goob10000 <65315624+goob10000@users.noreply.github.com> Date: Wed, 8 Apr 2026 16:08:58 -0700 Subject: [PATCH 5/5] in progress --- FullVehicleSim/Electrical/powertrain.py | 99 ++++++++++++++----------- 1 file changed, 54 insertions(+), 45 deletions(-) diff --git a/FullVehicleSim/Electrical/powertrain.py b/FullVehicleSim/Electrical/powertrain.py index 30b5699..97bbd11 100644 --- a/FullVehicleSim/Electrical/powertrain.py +++ b/FullVehicleSim/Electrical/powertrain.py @@ -35,49 +35,58 @@ def calcMaxPower(voltage:float) -> float: return Parameters["tractiveIMax"] * voltage -## Not working, needs changes def calcVoltage(worldArray:np.ndarray, step:int) -> float: -# delta = 1 / Parameters["stepsPerSecond"] -# capacity_Ah = Parameters["cellCapacity_Ah"] -# soc = worldArray[step-1, varCharge] - -# sigma = Parameters["cellModelSigma"] -# hystGain = Parameters["hysteresisGain"] - -# # Sliding window: last 10 seconds of current -# I_hist = np.zeros(int(10 * Parameters["stepsPerSecond"])) -# I_hist[:max(0, step - len(I_hist))] = worldArray[max(0, step - len(I_hist)):step, varCurrent] # Get the current history up to the current step - -# # Hysteresis kernel -# t = Parameters["histeresisKernelLength"] -# kernel = np.exp(-(t**2) / (2 * sigma**2)) -# kernel /= np.sum(kernel) - -# def ocv_from_soc(self, soc): -# return 3.0 + 0.9 * soc + 0.25 * np.exp(-12 * (1 - soc))s - -# def sag(self, current): -# return 0.02 * current + 0.004 * (current ** 1.3) - -# def step(self, current): - -# # Update SOC -# self.SOC -= (current * self.dt) / (3600 * self.capacity_Ah) -# self.SOC = np.clip(self.SOC, 0.0, 1.0) - -# # -------- Sliding array logic -------- -# self.I_hist[:-1] = self.I_hist[1:] # shift old values -# self.I_hist[-1] = current # add new current - -# # Hysteresis voltage -# V_hyst = self.hyst_gain * np.sum(self.I_hist * self.kernel) - -# # Terminal voltage -# voltage = ( -# self.ocv_from_soc(self.SOC) -# - self.sag(current) * (1 - self.SOC) -# - V_hyst -# ) - -# return voltage - return 120.0 \ No newline at end of file + delta = 1 / Parameters["stepsPerSecond"] + capacity_Ah = Parameters["cellCapacity_Ah"] + soc = worldArray[step-1, varCharge] + + sigma = Parameters["cellModelSigma"] + hystGain = Parameters["hysteresisGain"] + + # Sliding window: last 10 seconds of current + I_hist = np.zeros(int(10 * Parameters["stepsPerSecond"])) + I_hist[:max(0, step - len(I_hist))] = worldArray[max(0, step - len(I_hist)):step, varCurrent] # Get the current history up to the current step + + # Hysteresis kernel + t = Parameters["histeresisKernelLength"] + kernel = np.exp(-(t**2) / (2 * sigma**2)) + kernel /= np.sum(kernel) + + def ocv_from_soc(self, soc): + return 3.0 + 0.9 * soc + 0.25 * np.exp(-12 * (1 - soc))s + + def sag(self, current): + return 0.02 * current + 0.004 * (current ** 1.3) + + V_hyst = self.hyst_gain * np.sum(self.I_hist * self.kernel) + # Terminal voltage + voltage = ( + self.ocv_from_soc(self.SOC) + - self.sag(current) * (1 - self.SOC) + - V_hyst + ) + + return voltage + # return 120.0 + +# def step(self, current): + +# # Update SOC +# self.SOC -= (current * self.dt) / (3600 * self.capacity_Ah) +# self.SOC = np.clip(self.SOC, 0.0, 1.0) + +# # -------- Sliding array logic -------- +# self.I_hist[:-1] = self.I_hist[1:] # shift old values +# self.I_hist[-1] = current # add new current + +# # Hysteresis voltage +# V_hyst = self.hyst_gain * np.sum(self.I_hist * self.kernel) + +# # Terminal voltage +# voltage = ( +# self.ocv_from_soc(self.SOC) +# - self.sag(current) * (1 - self.SOC) +# - V_hyst +# ) + +# return voltage \ No newline at end of file