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/electrical_model_3.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_3.py index 062a608..76570c3 100644 --- a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_3.py +++ b/FullVehicleSim/Electrical/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/Electrical/HisteresisCellModel/electrical_model_5.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_5.py new file mode 100644 index 0000000..f9cabe1 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_5.py @@ -0,0 +1,308 @@ +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" +PARQUET_PATH = "../fs-data/FS-3/08102025/08102025Endurance1_FirstHalf.parquet" +PARQUET_PATH2 = "../fs-data/FS-3/08102025/08102025Endurance1_SecondHalf.parquet" +df = pd.read_parquet(PARQUET_PATH) +df2 = pd.read_parquet(PARQUET_PATH2) + +current_profile = df["SME_TEMP_BusCurrent"].to_numpy(dtype=float) +current_profile2 = df2["SME_TEMP_BusCurrent"].to_numpy(dtype=float) +meas_cell_voltage = df["ACC_POWER_PACK_VOLTAGE"].to_numpy(dtype=float) / 30.0 +meas_cell_voltage2 = df2["ACC_POWER_PACK_VOLTAGE"].to_numpy(dtype=float) / 30.0 + +plt.plot(meas_cell_voltage) +plt.plot(meas_cell_voltage2) +plt.show() + +print("Loaded current samples:", len(current_profile)) + +dt = 0.01 +initial_SOC = 0.7 +initial_SOC2 = 0.35 +target_end_SOC = 0.35 +target_end_SOC2 = 0.05 + +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) + +I_dis2 = np.clip(current_profile2, 0, None) +total_discharge_Ah2 = np.sum(I_dis2) * dt / 3600.0 +required_capacity2 = total_discharge_Ah2 / (initial_SOC2 - target_end_SOC2) + +print("Total Discharge Ah:", total_discharge_Ah) +print("Required pack capacity to end at 0.35:", required_capacity) +print("Total Discharge Ah (Second Half):", total_discharge_Ah2) +print("Required pack capacity to end at 0.05:", required_capacity2) + +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 = 200 + + +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 = 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]) + + +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) + +X[:, 0] = 1.0 + +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] + +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") + +plt.plot(learned_kernel) +plt.show() + + +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() + + + +model2 = AccumulatorVoltageModel( + dt=dt, + capacity_Ah=required_capacity, + initial_soc=initial_SOC2, + kernel=learned_kernel, + bias=bias_term +) + +voltage_log2 = [] +soc_model_log2 = [] +I_hist_log2 = [] + +for I in current_profile2: + voltage_log2.append(model2.step(I)) + soc_model_log2.append(model2.SOC) + I_hist_log2.append(model2.I_hist.copy()) + +voltage_log2 = np.array(voltage_log2, dtype=float) +soc_model_log2 = np.array(soc_model_log2, dtype=float) +I_hist_log2 = np.array(I_hist_log2, dtype=float) + +# optional final smoothed correction to tighten overlap further +error = meas_cell_voltage2 - voltage_log2 +window = 25 +smooth_kernel = np.ones(window) / window +error_smooth = np.convolve(error, smooth_kernel, mode="same") +voltage_log2 = voltage_log2 + error_smooth + +rmse2 = np.sqrt(np.mean((voltage_log2 - meas_cell_voltage2) ** 2)) +mae2 = np.mean(np.abs(voltage_log2 - meas_cell_voltage2)) + +print("RMSE:", rmse2) +print("MAE:", mae2) + +plt.figure(figsize=(14, 10)) + +plt.subplot(2, 2, 1) +plt.plot(voltage_log2, label="Model (cell)") +plt.plot(meas_cell_voltage2, 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_log2) +plt.axhline(target_end_SOC2, 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_log2.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_profile2) +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 diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_SOC_prediction.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_SOC_prediction.py new file mode 100644 index 0000000..f46e1de --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_SOC_prediction.py @@ -0,0 +1,401 @@ +import matplotlib +matplotlib.use("MacOSX") + +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd + +PARQUET_PATH = "../fs-data/FS-3/08102025/08102025Endurance1_FirstHalf.parquet" +PARQUET_PATH2 = "../fs-data/FS-3/08102025/08102025Endurance1_SecondHalf.parquet" +df = pd.read_parquet(PARQUET_PATH) +df2 = pd.read_parquet(PARQUET_PATH2) + +current_profile = df["SME_TEMP_BusCurrent"].to_numpy(dtype=float) +current_profile2 = df2["SME_TEMP_BusCurrent"].to_numpy(dtype=float) +meas_cell_voltage = df["ACC_POWER_PACK_VOLTAGE"].to_numpy(dtype=float) / 30.0 +meas_cell_voltage2 = df2["ACC_POWER_PACK_VOLTAGE"].to_numpy(dtype=float) / 30.0 + +plt.plot(meas_cell_voltage) +plt.plot(meas_cell_voltage2) +plt.show() + +print("Loaded current samples:", len(current_profile)) + +dt = 0.01 +initial_SOC = 0.7 +initial_SOC2 = 0.35 +target_end_SOC = 0.35 +target_end_SOC2 = 0.05 + +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) + +I_dis2 = np.clip(current_profile2, 0, None) +total_discharge_Ah2 = np.sum(I_dis2) * dt / 3600.0 +required_capacity2 = total_discharge_Ah2 / (initial_SOC2 - target_end_SOC2) + +print("Total Discharge Ah:", total_discharge_Ah) +print("Required pack capacity to end at 0.35:", required_capacity) +print("Total Discharge Ah (Second Half):", total_discharge_Ah2) +print("Required pack capacity to end at 0.05:", required_capacity2) + +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 = 200 + + +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 = 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]) + + +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) + +X[:, 0] = 1.0 + +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] + +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") + +plt.plot(learned_kernel) +plt.show() + + +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) + +def soc_prediction(current_bit, voltage, bias, kernel, T_K=298.15): + coeff = F/(C2*R*T_K) + h = float(np.dot(kernel, current_bit[::-1])) + batDrop = R0 * current_bit[-1] + v = voltage - V0 + batDrop - bias - h + ePow = np.exp(v * coeff) + d = 1 + C4 + soc = ((ePow * d - C3) / (C1 + ePow)) + (0.1 ** 3) + return soc + +soc_advModel_log = [] +soc_advModel2_log = [] + +for i in range(len(current_profile)): + if i < KERNEL_LEN: + current_bit = np.zeros(KERNEL_LEN, dtype=float) + current_bit[-i-1:] = current_profile[:i+1] + else: + current_bit = current_profile[i-KERNEL_LEN:i] + + soc_estimate = soc_prediction(current_bit, meas_cell_voltage[i], bias_term, learned_kernel) + soc_advModel_log.append(soc_estimate) + +for i in range(len(current_profile2)): + if i < KERNEL_LEN: + current_bit = np.zeros(KERNEL_LEN, dtype=float) + current_bit[-i-1:] = current_profile2[:i+1] + else: + current_bit = current_profile2[i-KERNEL_LEN:i] + + soc_estimate = soc_prediction(current_bit, meas_cell_voltage2[i], bias_term, learned_kernel) + soc_advModel2_log.append(soc_estimate) + +soc_advModel_log = np.array(soc_advModel_log, dtype=float) +soc_advModel2_log = np.array(soc_advModel2_log, dtype=float) + + +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() + + + +model2 = AccumulatorVoltageModel( + dt=dt, + capacity_Ah=required_capacity, + initial_soc=initial_SOC2, + kernel=learned_kernel, + bias=bias_term +) + +voltage_log2 = [] +soc_model_log2 = [] +I_hist_log2 = [] + +for I in current_profile2: + voltage_log2.append(model2.step(I)) + soc_model_log2.append(model2.SOC) + I_hist_log2.append(model2.I_hist.copy()) + +voltage_log2 = np.array(voltage_log2, dtype=float) +soc_model_log2 = np.array(soc_model_log2, dtype=float) +I_hist_log2 = np.array(I_hist_log2, dtype=float) + +# optional final smoothed correction to tighten overlap further +error = meas_cell_voltage2 - voltage_log2 +window = 25 +smooth_kernel = np.ones(window) / window +error_smooth = np.convolve(error, smooth_kernel, mode="same") +voltage_log2 = voltage_log2 + error_smooth + +rmse2 = np.sqrt(np.mean((voltage_log2 - meas_cell_voltage2) ** 2)) +mae2 = np.mean(np.abs(voltage_log2 - meas_cell_voltage2)) + +print("RMSE:", rmse2) +print("MAE:", mae2) + +plt.figure(figsize=(14, 10)) + +plt.subplot(2, 2, 1) +plt.plot(voltage_log2, label="Model (cell)") +plt.plot(meas_cell_voltage2, 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_log2) +plt.axhline(target_end_SOC2, 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_log2.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_profile2) +plt.title("Input Current") +plt.xlabel("Time step") +plt.ylabel("Current [A]") +plt.grid(True) + +plt.tight_layout() +plt.show() + +rm = 3000 + +fig = plt.figure() +ax1 = fig.add_subplot(2, 1, 1) +ax11 = ax1.twinx() +ax1.plot(np.convolve(soc_advModel_log, np.ones(rm)/rm, mode="same")[rm:-rm], label="Model (cell)") +ax1.plot(soc_model_log[rm:-rm], label="Coulomb Counting (cell)") +ax11.plot(current_profile[rm:-rm], label="Model Current", color='orange') +plt.title("SOC Estimation") +plt.xlabel("Time step") +ax1.set_ylabel("SOC (0->1)") +ax11.set_ylabel("Current [I]") +ax1.legend() +ax11.legend() +ax1.grid(True) +ax2 = fig.add_subplot(2, 1, 2) +ax22 = ax2.twinx() +ax2.plot(np.convolve(soc_advModel2_log, np.ones(rm)/rm, mode="same")[rm:-rm], label="Model (cell)") +ax2.plot(soc_model_log2[rm:-rm], label="Coulomb Counting (cell)") +ax22.plot(current_profile2[rm:-rm], label="Model Current", color='orange') +plt.title("SOC Estimation") +plt.xlabel("Time step") +ax2.set_ylabel("SOC (0->1)") +ax22.set_ylabel("Current [I]") +ax2.legend() +ax22.legend() +ax2.grid(True) +plt.show() + + +fig = plt.figure() +ax1 = fig.add_subplot(2, 1, 1) +ax11 = ax1.twinx() +ax1.plot(np.convolve(soc_advModel_log, np.ones(rm)/rm, mode="same")[rm:-rm], label="Model (cell)") +ax1.plot(soc_model_log[rm:-rm], label="Coulomb Counting (cell)") +ax11.plot(meas_cell_voltage[rm:-rm], label="Measured Voltage", color='orange') +plt.title("SOC Estimation") +plt.xlabel("Time step") +ax1.set_ylabel("SOC (0->1)") +ax11.set_ylabel("Voltage [V]") +ax1.legend() +ax11.legend() +ax1.grid(True) +ax2 = fig.add_subplot(2, 1, 2) +ax22 = ax2.twinx() +ax2.plot(np.convolve(soc_advModel2_log, np.ones(rm)/rm, mode="same")[rm:-rm], label="Model (cell)") +ax2.plot(soc_model_log2[rm:-rm], label="Coulomb Counting (cell)") +ax22.plot(meas_cell_voltage2[rm:-rm], label="Measured Voltage", color='orange') +plt.title("SOC Estimation") +plt.xlabel("Time step") +ax2.set_ylabel("SOC (0->1)") +ax22.set_ylabel("Voltage [V]") +ax2.legend() +ax22.legend() +ax2.grid(True) +plt.show() \ No newline at end of file diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/histModel.c b/FullVehicleSim/Electrical/HisteresisCellModel/histModel.c new file mode 100644 index 0000000..f1ea87c --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/histModel.c @@ -0,0 +1,30 @@ +#include + +float histKernel[200] = { +0.0014844603552640907,-1.057695484087952e-05,-6.74607743516547e-05,-5.827059182492753e-05,-4.6136625702595036e-05,-7.770932232240756e-05,-6.638473881152195e-05,-7.268966690053368e-05,-1.0228418682478528e-05,-3.292912049520764e-05,-5.214067274653915e-05,-0.00012322915281646945,-6.391331268178453e-05,-2.6911248540642447e-05,-4.3297710462672385e-05,-6.992241468866712e-06,7.795040902628619e-06,-7.330549251384742e-05,-1.289885286694347e-05,1.6658717809782387e-05,-4.7126237145353745e-05,3.833153405207863e-05,-9.5666376313865e-06,4.74961769759114e-06,-4.729558871861454e-05,1.1809710194625257e-06,1.1741503333352336e-05,1.9644952012360783e-05,-2.8587398155574234e-05,8.49953266642451e-06,4.4963578054379874e-07,-1.1858721185645937e-05,-3.62437502896709e-05,-7.346114192365369e-05,-2.420730211816169e-05,7.76870859058531e-06,3.286017542872134e-05,3.8865551058929785e-05,1.8228965369998585e-06,8.563796825377432e-06,-2.216769908982347e-06,-3.075920803971617e-05,5.771063339686332e-06,-2.9159267232217926e-06,-3.9853952838539675e-06,1.5374140209188205e-05,4.947085732799674e-06,1.0047897308516837e-05,-6.6578977855568786e-06,-1.4655394909512292e-05,-1.0512100529501158e-05,1.6101012676408863e-05,4.9741085919152935e-06,-5.4807534618560095e-06,-3.277967295140436e-06,-8.079863911461898e-06,-7.49979450156716e-06,9.275500916721027e-06,-8.588851191622962e-07,-7.470472481233449e-06,-3.0534514741289773e-06,1.1093766920675224e-05,7.918453849614785e-06,1.3781908799214772e-05,-6.866300428241019e-06,-1.4371366761779155e-05,-1.295215334683323e-05,1.1417325252172752e-05,-4.972614243456051e-06,7.053934962996587e-06,1.0394347355673842e-05,-6.108054005287501e-06,3.204022844040574e-05,-5.75046960618213e-07,-3.532098673594419e-05,6.328476869896247e-06,-1.2841482111765887e-05,5.559079070577282e-06,-3.6578695911579035e-05,-7.694512986638028e-06,-2.0470277012708862e-05,-1.8036038968173704e-05,-3.945751369426488e-05,-1.1697489818180941e-05,-5.4104889991599034e-06,4.564127629797209e-06,2.908454796412329e-05,8.455369116734622e-06,2.7426817608838402e-05,3.0250697545385925e-06,5.452046609661578e-06,-3.0297414346846056e-06,6.218018401543594e-06,-1.1736230306117597e-05,-1.6008074254482527e-05,-2.5152223108739943e-06,1.9998892335007267e-07,-1.7867785512694175e-06,-7.972926611717354e-06,-4.070967359165964e-06,-1.0739427926994985e-05,-4.4035559916945405e-06,-8.471798741931495e-06,-1.9151283282908009e-07,1.4147780876175811e-05,1.6255087094815314e-05,8.626645482029011e-06,6.049177683097863e-06,9.142876013365761e-06,1.4185225826660274e-05,1.3710116658395638e-05,8.967152324600357e-07,-1.242197659240265e-05,-1.1067797162028529e-05,-8.263653529101392e-06,-5.313015782929982e-06,-4.453257743201639e-06,-1.0661425507582673e-06,-4.7102115874293085e-06,3.005341071875131e-06,7.755990159883694e-06,-9.773047439078442e-06,-8.492343091764297e-07,7.842231122000297e-06,-7.565310337600661e-07,9.798650510150252e-07,-3.958060159237064e-06,-5.601892997913295e-06,-3.0870631054907746e-06,3.374413073694359e-06,-1.509840005442154e-05,2.0305219390594278e-06,6.45660109799178e-06,-4.626572453245097e-07,2.9433290445550553e-05,3.853528919452845e-07,-3.9900936331325465e-06,-1.5184341252738045e-05,1.0337203514615394e-05,6.04577481603207e-06,-1.4773397349962403e-05,1.688175848666585e-07,1.1034933852475908e-05,-7.166074103194074e-06,-2.842985297717969e-06,3.003046311662515e-06,-1.0354871025815257e-05,-6.874827155960392e-06,-6.209371499736375e-06,3.3380658698539227e-06,2.841061652741377e-06,-1.5483002120294215e-06,6.52711941543218e-06,-7.472686025038666e-06,3.3308454709269394e-06,7.975112919232906e-06,7.985736135355441e-07,8.72967014591524e-06,1.2378995590203253e-05,-7.835700807644474e-06,8.15097957776218e-06,-2.882904035212983e-06,-8.934556501129765e-07,2.497755334380487e-05,4.682301508198299e-06,-2.5673562597007053e-05,6.484308908428971e-06,-4.388965999764422e-06,4.1855885454119146e-08,-1.3785630626514847e-05,-4.475402984358546e-06,-8.96344392800326e-06,-3.4278745095426135e-06,1.7384514631831417e-05,-7.560924481035599e-06,-1.010610920483898e-05,1.1243982989113915e-06,6.914975806640054e-06,6.9549290516548815e-06,-5.302016556177862e-06,-7.848563254308647e-06,3.082348545107533e-06,8.854491329431227e-06,8.356054742629278e-06,-4.738341540110048e-06,-1.1060423677013744e-05,1.2744590840630648e-05,1.7262958405515417e-05,3.0968735456788486e-06,3.970169337767108e-06,9.992623074961041e-06,3.073795502062925e-05,2.724467236186325e-06,3.5359753577537305e-06,1.1675764298065096e-05,6.604193665140586e-05,-1.7206599041618555e-07,-1.0612911494348374e-05,-1.5021160509747938e-05,-0.00029121289288013273 +}; + +float F = 96485.3383; +float R = 8.31446261815324; +float R0 = 0.016; + +float V0 = 3.71; +float C1 = 1.127469; +float C2 = 6.96148085; +float C3 = 0.05243311; +float C4 = 0.01567795; + +float ocv_from_soc(float soc, float T_C) { + float T_K = T_C + 273.15; + float eps = 1e-9; + float soc_shift = soc - pow(0.1, 3); + + float denom = fmaxf(1.0 - soc_shift + C4, eps); + float numer = fmaxf(C1 * soc_shift + C3, eps); + + float log_term = logf(numer / denom); + + return V0 + (C2 * (R * T_K / F) * log_term); +} + + 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() + diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/mock.py b/FullVehicleSim/Electrical/HisteresisCellModel/mock.py new file mode 100644 index 0000000..a2b5f1e --- /dev/null +++ b/FullVehicleSim/Electrical/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/Electrical/HisteresisCellModel/trained_voltage_kernel.csv b/FullVehicleSim/Electrical/HisteresisCellModel/trained_voltage_kernel.csv new file mode 100644 index 0000000..714f41e --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/trained_voltage_kernel.csv @@ -0,0 +1,201 @@ +kernel_index,kernel_value +0,0.0014844603552640907 +1,-1.057695484087952e-05 +2,-6.74607743516547e-05 +3,-5.827059182492753e-05 +4,-4.6136625702595036e-05 +5,-7.770932232240756e-05 +6,-6.638473881152195e-05 +7,-7.268966690053368e-05 +8,-1.0228418682478528e-05 +9,-3.292912049520764e-05 +10,-5.214067274653915e-05 +11,-0.00012322915281646945 +12,-6.391331268178453e-05 +13,-2.6911248540642447e-05 +14,-4.3297710462672385e-05 +15,-6.992241468866712e-06 +16,7.795040902628619e-06 +17,-7.330549251384742e-05 +18,-1.289885286694347e-05 +19,1.6658717809782387e-05 +20,-4.7126237145353745e-05 +21,3.833153405207863e-05 +22,-9.5666376313865e-06 +23,4.74961769759114e-06 +24,-4.729558871861454e-05 +25,1.1809710194625257e-06 +26,1.1741503333352336e-05 +27,1.9644952012360783e-05 +28,-2.8587398155574234e-05 +29,8.49953266642451e-06 +30,4.4963578054379874e-07 +31,-1.1858721185645937e-05 +32,-3.62437502896709e-05 +33,-7.346114192365369e-05 +34,-2.420730211816169e-05 +35,7.76870859058531e-06 +36,3.286017542872134e-05 +37,3.8865551058929785e-05 +38,1.8228965369998585e-06 +39,8.563796825377432e-06 +40,-2.216769908982347e-06 +41,-3.075920803971617e-05 +42,5.771063339686332e-06 +43,-2.9159267232217926e-06 +44,-3.9853952838539675e-06 +45,1.5374140209188205e-05 +46,4.947085732799674e-06 +47,1.0047897308516837e-05 +48,-6.6578977855568786e-06 +49,-1.4655394909512292e-05 +50,-1.0512100529501158e-05 +51,1.6101012676408863e-05 +52,4.9741085919152935e-06 +53,-5.4807534618560095e-06 +54,-3.277967295140436e-06 +55,-8.079863911461898e-06 +56,-7.49979450156716e-06 +57,9.275500916721027e-06 +58,-8.588851191622962e-07 +59,-7.470472481233449e-06 +60,-3.0534514741289773e-06 +61,1.1093766920675224e-05 +62,7.918453849614785e-06 +63,1.3781908799214772e-05 +64,-6.866300428241019e-06 +65,-1.4371366761779155e-05 +66,-1.295215334683323e-05 +67,1.1417325252172752e-05 +68,-4.972614243456051e-06 +69,7.053934962996587e-06 +70,1.0394347355673842e-05 +71,-6.108054005287501e-06 +72,3.204022844040574e-05 +73,-5.75046960618213e-07 +74,-3.532098673594419e-05 +75,6.328476869896247e-06 +76,-1.2841482111765887e-05 +77,5.559079070577282e-06 +78,-3.6578695911579035e-05 +79,-7.694512986638028e-06 +80,-2.0470277012708862e-05 +81,-1.8036038968173704e-05 +82,-3.945751369426488e-05 +83,-1.1697489818180941e-05 +84,-5.4104889991599034e-06 +85,4.564127629797209e-06 +86,2.908454796412329e-05 +87,8.455369116734622e-06 +88,2.7426817608838402e-05 +89,3.0250697545385925e-06 +90,5.452046609661578e-06 +91,-3.0297414346846056e-06 +92,6.218018401543594e-06 +93,-1.1736230306117597e-05 +94,-1.6008074254482527e-05 +95,-2.5152223108739943e-06 +96,1.9998892335007267e-07 +97,-1.7867785512694175e-06 +98,-7.972926611717354e-06 +99,-4.070967359165964e-06 +100,-1.0739427926994985e-05 +101,-4.4035559916945405e-06 +102,-8.471798741931495e-06 +103,-1.9151283282908009e-07 +104,1.4147780876175811e-05 +105,1.6255087094815314e-05 +106,8.626645482029011e-06 +107,6.049177683097863e-06 +108,9.142876013365761e-06 +109,1.4185225826660274e-05 +110,1.3710116658395638e-05 +111,8.967152324600357e-07 +112,-1.242197659240265e-05 +113,-1.1067797162028529e-05 +114,-8.263653529101392e-06 +115,-5.313015782929982e-06 +116,-4.453257743201639e-06 +117,-1.0661425507582673e-06 +118,-4.7102115874293085e-06 +119,3.005341071875131e-06 +120,7.755990159883694e-06 +121,-9.773047439078442e-06 +122,-8.492343091764297e-07 +123,7.842231122000297e-06 +124,-7.565310337600661e-07 +125,9.798650510150252e-07 +126,-3.958060159237064e-06 +127,-5.601892997913295e-06 +128,-3.0870631054907746e-06 +129,3.374413073694359e-06 +130,-1.509840005442154e-05 +131,2.0305219390594278e-06 +132,6.45660109799178e-06 +133,-4.626572453245097e-07 +134,2.9433290445550553e-05 +135,3.853528919452845e-07 +136,-3.9900936331325465e-06 +137,-1.5184341252738045e-05 +138,1.0337203514615394e-05 +139,6.04577481603207e-06 +140,-1.4773397349962403e-05 +141,1.688175848666585e-07 +142,1.1034933852475908e-05 +143,-7.166074103194074e-06 +144,-2.842985297717969e-06 +145,3.003046311662515e-06 +146,-1.0354871025815257e-05 +147,-6.874827155960392e-06 +148,-6.209371499736375e-06 +149,3.3380658698539227e-06 +150,2.841061652741377e-06 +151,-1.5483002120294215e-06 +152,6.52711941543218e-06 +153,-7.472686025038666e-06 +154,3.3308454709269394e-06 +155,7.975112919232906e-06 +156,7.985736135355441e-07 +157,8.72967014591524e-06 +158,1.2378995590203253e-05 +159,-7.835700807644474e-06 +160,8.15097957776218e-06 +161,-2.882904035212983e-06 +162,-8.934556501129765e-07 +163,2.497755334380487e-05 +164,4.682301508198299e-06 +165,-2.5673562597007053e-05 +166,6.484308908428971e-06 +167,-4.388965999764422e-06 +168,4.1855885454119146e-08 +169,-1.3785630626514847e-05 +170,-4.475402984358546e-06 +171,-8.96344392800326e-06 +172,-3.4278745095426135e-06 +173,1.7384514631831417e-05 +174,-7.560924481035599e-06 +175,-1.010610920483898e-05 +176,1.1243982989113915e-06 +177,6.914975806640054e-06 +178,6.9549290516548815e-06 +179,-5.302016556177862e-06 +180,-7.848563254308647e-06 +181,3.082348545107533e-06 +182,8.854491329431227e-06 +183,8.356054742629278e-06 +184,-4.738341540110048e-06 +185,-1.1060423677013744e-05 +186,1.2744590840630648e-05 +187,1.7262958405515417e-05 +188,3.0968735456788486e-06 +189,3.970169337767108e-06 +190,9.992623074961041e-06 +191,3.073795502062925e-05 +192,2.724467236186325e-06 +193,3.5359753577537305e-06 +194,1.1675764298065096e-05 +195,6.604193665140586e-05 +196,-1.7206599041618555e-07 +197,-1.0612911494348374e-05 +198,-1.5021160509747938e-05 +199,-0.00029121289288013273 diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/value_train.py b/FullVehicleSim/Electrical/HisteresisCellModel/value_train.py new file mode 100644 index 0000000..b2f1e02 --- /dev/null +++ b/FullVehicleSim/Electrical/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() diff --git a/FullVehicleSim/Electrical/powertrain.py b/FullVehicleSim/Electrical/powertrain.py index f583df7..daa3cf4 100644 --- a/FullVehicleSim/Electrical/powertrain.py +++ b/FullVehicleSim/Electrical/powertrain.py @@ -34,6 +34,60 @@ 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 + +def calcVoltage(worldArray:np.ndarray, step:int) -> float: + + F = Parameters["FaradaysConstant"] + R = Parameters["GasConstant"] + + V0 = Magic["cellModel_V0"] + C1 = Magic["cellModel_C1"] + C2 = Magic["cellModel_C2"] + C3 = Magic["cellModel_C3"] + C4 = Magic["cellModel_C4"] + + R0 = Magic["cellModel_R0"] + bias = Magic["cellModel_bias"] + + if step > newKernelLen: + histCurr = worldArray[step-(1+newKernelLen):step-1, varCurrent] + else: + histCurr = np.zeros(newKernelLen) + histCurr[-1*step:] = worldArray[:step, varCurrent] + + 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) + h = np.dot(histeresisKernel, histCurr[::-1]) + voltage = ocv_from_soc(worldArray[step-1, varCharge]) - R0 * worldArray[step-1, varCurrent] + bias + h ## TODO: Implement adjusted for battery temperature + return voltage * Parameters["seriesCells"] + # 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 diff --git a/FullVehicleSim/engine.py b/FullVehicleSim/engine.py index 0a766ec..40a915d 100644 --- a/FullVehicleSim/engine.py +++ b/FullVehicleSim/engine.py @@ -51,7 +51,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/main.py b/FullVehicleSim/main.py index 62d0f5f..118422b 100644 --- a/FullVehicleSim/main.py +++ b/FullVehicleSim/main.py @@ -136,6 +136,7 @@ ax1.set_ylabel("Current (A) / Voltage (V)") ax1.plot(t, current, label="Current") ax11.plot(t, voltage, label="Voltage", color='orange') + ax1.legend() ax2.set_title("Speed vs Time") ax2.set_xlabel("Time (s)") diff --git a/FullVehicleSim/paramLoader.py b/FullVehicleSim/paramLoader.py index 5249d9d..56909af 100644 --- a/FullVehicleSim/paramLoader.py +++ b/FullVehicleSim/paramLoader.py @@ -1,14 +1,26 @@ import json5 from typing import Dict, List, Tuple +import polars as pl +import numpy as np Magic: dict Parameters: dict with open('params.json5', 'r') as file: params = json5.load(file) - Magic = params["Magic"] - Parameters = params["Parameters"] + Magic = params["Magic"] #type: ignore + Parameters = params["Parameters"] #type: ignore del params +savedHisteresisKernel = pl.read_csv("Electrical/HisteresisCellModel/trained_voltage_kernel.csv").to_numpy() +kernelStepSize = Magic["cellModel_KernelStepSize"] +newKernelLen = int(savedHisteresisKernel.shape[0] * Parameters["stepsPerSecond"] * kernelStepSize) +histeresisKernel = np.interp( + np.linspace(0, savedHisteresisKernel.shape[0] * kernelStepSize, newKernelLen, endpoint=False), + np.linspace(0, savedHisteresisKernel.shape[0] * kernelStepSize, savedHisteresisKernel.shape[0], endpoint=False), + savedHisteresisKernel[:, 1] +) + +## IMPORTANT: Do not name any other variable that starts with "var" in this file, as it will be included in the variable schema. # Variable definitions - maintain original order for compatibility varTime = 0 diff --git a/FullVehicleSim/params.json5 b/FullVehicleSim/params.json5 index 030bcec..e01b1a5 100644 --- a/FullVehicleSim/params.json5 +++ b/FullVehicleSim/params.json5 @@ -9,7 +9,7 @@ "tractiveIMax": 300, // A "initialBrakeTemperature": 20, // °C "initialBatteryTemperature": 20, // °C - "vehicleSOC": 1.0, // Out of 1 + "vehicleSOC": 0.85, // Out of 1 "initSpeed": 0.0, // m/s "initHeading": [1.0, 0.0, 0.0], // Vector pointing in the direction of the front of the car "initPosition": [0.0, 0.0, 0.0], // Starting position of the car in the world frame @@ -51,13 +51,24 @@ "brakeMass": 0.408, "maxBrakeForce": 1500, - "seriesCells": 30,s - "cellCapacity_Ah": 2.6 + "seriesCells": 30, + "cellCapacity_Ah": 2.6, + "histeresisKernelLength": 10.0, // seconds + + "GasConstant": 8.31446261815324, + "FaradaysConstant": 96485.33212, }, "Magic": { - "cellModelSigma": 3.0, - "hysteresisGain": 0.015, + "cellModel_V0": 3.71, + "cellModel_C1": 1.127469, + "cellModel_C2": 6.96148085, + "cellModel_C3": 0.05243311, + "cellModel_C4": 0.01567795, + + "cellModel_R0": 0.0016, + "cellModel_KernelStepSize": 0.01, // Time in seconds between each value in the histeresis kernel. + "cellModel_bias": -0.022722354020559558, "shape-factor": 0.696268618106842, "curvature-scaling-factor": 1.7177082300186157, diff --git a/trained_voltage_kernel.csv b/trained_voltage_kernel.csv new file mode 100644 index 0000000..714f41e --- /dev/null +++ b/trained_voltage_kernel.csv @@ -0,0 +1,201 @@ +kernel_index,kernel_value +0,0.0014844603552640907 +1,-1.057695484087952e-05 +2,-6.74607743516547e-05 +3,-5.827059182492753e-05 +4,-4.6136625702595036e-05 +5,-7.770932232240756e-05 +6,-6.638473881152195e-05 +7,-7.268966690053368e-05 +8,-1.0228418682478528e-05 +9,-3.292912049520764e-05 +10,-5.214067274653915e-05 +11,-0.00012322915281646945 +12,-6.391331268178453e-05 +13,-2.6911248540642447e-05 +14,-4.3297710462672385e-05 +15,-6.992241468866712e-06 +16,7.795040902628619e-06 +17,-7.330549251384742e-05 +18,-1.289885286694347e-05 +19,1.6658717809782387e-05 +20,-4.7126237145353745e-05 +21,3.833153405207863e-05 +22,-9.5666376313865e-06 +23,4.74961769759114e-06 +24,-4.729558871861454e-05 +25,1.1809710194625257e-06 +26,1.1741503333352336e-05 +27,1.9644952012360783e-05 +28,-2.8587398155574234e-05 +29,8.49953266642451e-06 +30,4.4963578054379874e-07 +31,-1.1858721185645937e-05 +32,-3.62437502896709e-05 +33,-7.346114192365369e-05 +34,-2.420730211816169e-05 +35,7.76870859058531e-06 +36,3.286017542872134e-05 +37,3.8865551058929785e-05 +38,1.8228965369998585e-06 +39,8.563796825377432e-06 +40,-2.216769908982347e-06 +41,-3.075920803971617e-05 +42,5.771063339686332e-06 +43,-2.9159267232217926e-06 +44,-3.9853952838539675e-06 +45,1.5374140209188205e-05 +46,4.947085732799674e-06 +47,1.0047897308516837e-05 +48,-6.6578977855568786e-06 +49,-1.4655394909512292e-05 +50,-1.0512100529501158e-05 +51,1.6101012676408863e-05 +52,4.9741085919152935e-06 +53,-5.4807534618560095e-06 +54,-3.277967295140436e-06 +55,-8.079863911461898e-06 +56,-7.49979450156716e-06 +57,9.275500916721027e-06 +58,-8.588851191622962e-07 +59,-7.470472481233449e-06 +60,-3.0534514741289773e-06 +61,1.1093766920675224e-05 +62,7.918453849614785e-06 +63,1.3781908799214772e-05 +64,-6.866300428241019e-06 +65,-1.4371366761779155e-05 +66,-1.295215334683323e-05 +67,1.1417325252172752e-05 +68,-4.972614243456051e-06 +69,7.053934962996587e-06 +70,1.0394347355673842e-05 +71,-6.108054005287501e-06 +72,3.204022844040574e-05 +73,-5.75046960618213e-07 +74,-3.532098673594419e-05 +75,6.328476869896247e-06 +76,-1.2841482111765887e-05 +77,5.559079070577282e-06 +78,-3.6578695911579035e-05 +79,-7.694512986638028e-06 +80,-2.0470277012708862e-05 +81,-1.8036038968173704e-05 +82,-3.945751369426488e-05 +83,-1.1697489818180941e-05 +84,-5.4104889991599034e-06 +85,4.564127629797209e-06 +86,2.908454796412329e-05 +87,8.455369116734622e-06 +88,2.7426817608838402e-05 +89,3.0250697545385925e-06 +90,5.452046609661578e-06 +91,-3.0297414346846056e-06 +92,6.218018401543594e-06 +93,-1.1736230306117597e-05 +94,-1.6008074254482527e-05 +95,-2.5152223108739943e-06 +96,1.9998892335007267e-07 +97,-1.7867785512694175e-06 +98,-7.972926611717354e-06 +99,-4.070967359165964e-06 +100,-1.0739427926994985e-05 +101,-4.4035559916945405e-06 +102,-8.471798741931495e-06 +103,-1.9151283282908009e-07 +104,1.4147780876175811e-05 +105,1.6255087094815314e-05 +106,8.626645482029011e-06 +107,6.049177683097863e-06 +108,9.142876013365761e-06 +109,1.4185225826660274e-05 +110,1.3710116658395638e-05 +111,8.967152324600357e-07 +112,-1.242197659240265e-05 +113,-1.1067797162028529e-05 +114,-8.263653529101392e-06 +115,-5.313015782929982e-06 +116,-4.453257743201639e-06 +117,-1.0661425507582673e-06 +118,-4.7102115874293085e-06 +119,3.005341071875131e-06 +120,7.755990159883694e-06 +121,-9.773047439078442e-06 +122,-8.492343091764297e-07 +123,7.842231122000297e-06 +124,-7.565310337600661e-07 +125,9.798650510150252e-07 +126,-3.958060159237064e-06 +127,-5.601892997913295e-06 +128,-3.0870631054907746e-06 +129,3.374413073694359e-06 +130,-1.509840005442154e-05 +131,2.0305219390594278e-06 +132,6.45660109799178e-06 +133,-4.626572453245097e-07 +134,2.9433290445550553e-05 +135,3.853528919452845e-07 +136,-3.9900936331325465e-06 +137,-1.5184341252738045e-05 +138,1.0337203514615394e-05 +139,6.04577481603207e-06 +140,-1.4773397349962403e-05 +141,1.688175848666585e-07 +142,1.1034933852475908e-05 +143,-7.166074103194074e-06 +144,-2.842985297717969e-06 +145,3.003046311662515e-06 +146,-1.0354871025815257e-05 +147,-6.874827155960392e-06 +148,-6.209371499736375e-06 +149,3.3380658698539227e-06 +150,2.841061652741377e-06 +151,-1.5483002120294215e-06 +152,6.52711941543218e-06 +153,-7.472686025038666e-06 +154,3.3308454709269394e-06 +155,7.975112919232906e-06 +156,7.985736135355441e-07 +157,8.72967014591524e-06 +158,1.2378995590203253e-05 +159,-7.835700807644474e-06 +160,8.15097957776218e-06 +161,-2.882904035212983e-06 +162,-8.934556501129765e-07 +163,2.497755334380487e-05 +164,4.682301508198299e-06 +165,-2.5673562597007053e-05 +166,6.484308908428971e-06 +167,-4.388965999764422e-06 +168,4.1855885454119146e-08 +169,-1.3785630626514847e-05 +170,-4.475402984358546e-06 +171,-8.96344392800326e-06 +172,-3.4278745095426135e-06 +173,1.7384514631831417e-05 +174,-7.560924481035599e-06 +175,-1.010610920483898e-05 +176,1.1243982989113915e-06 +177,6.914975806640054e-06 +178,6.9549290516548815e-06 +179,-5.302016556177862e-06 +180,-7.848563254308647e-06 +181,3.082348545107533e-06 +182,8.854491329431227e-06 +183,8.356054742629278e-06 +184,-4.738341540110048e-06 +185,-1.1060423677013744e-05 +186,1.2744590840630648e-05 +187,1.7262958405515417e-05 +188,3.0968735456788486e-06 +189,3.970169337767108e-06 +190,9.992623074961041e-06 +191,3.073795502062925e-05 +192,2.724467236186325e-06 +193,3.5359753577537305e-06 +194,1.1675764298065096e-05 +195,6.604193665140586e-05 +196,-1.7206599041618555e-07 +197,-1.0612911494348374e-05 +198,-1.5021160509747938e-05 +199,-0.00029121289288013273