diff --git a/Data/fs4voltagemodel.py b/Data/fs4voltagemodel.py new file mode 100644 index 0000000..c1164cd --- /dev/null +++ b/Data/fs4voltagemodel.py @@ -0,0 +1,136 @@ +import numpy as np +import matplotlib.pyplot as plt +from dataclasses import dataclass + +# ====================================================== +# PARAMETERS +# ====================================================== +Q_rated = 2.5 * 3600 # Coulombs (2.5 Ah cell) +dt = 1.0 # timestep [s] +N = 300 # number of steps + +# ====================================================== +# VEHICLE / BATTERY STATE +# ====================================================== +@dataclass +class VehicleState: + soc: float # State of Charge [0..1] + voltage: float # Terminal voltage [V] + + # Reserved for future revised model + temp_c: float = 25.0 + hyst: float = 0.0 + v_rc: float = 0.0 + + +# ====================================================== +# VOLTAGE UPDATE TEMPLATE ← THIS IS THE DELIVERABLE +# ====================================================== +def update_voltage_template(prev_current: float, state: VehicleState) -> float: + """ + Template for revised voltage updating (ECM baseline). + + Model: + V = OCV(SOC) - I*R_internal + + Inputs: + prev_current: current from previous timestep [A] (+ discharge, - charge) + state.soc: SOC in [0,1] + + Output: + new terminal voltage [V] + + TODO later: + - Replace OCV curve with real discharge curve fit + - Make R_internal depend on temperature/SOC + - Add hysteresis / RC recovery (sag + slow rebound) + """ + soc = float(state.soc) + + # 1) OCV curve (placeholder, but reasonable): + # At SOC=1.0 -> ~4.2V, at SOC=0.0 -> ~3.0V + ocv = 3.0 + 1.2 * soc + + # 2) Internal resistance (placeholder constant) + r_internal = 0.015 # Ohms (cell-level example) + + # 3) Terminal voltage + v_new = ocv - prev_current * r_internal + + # 4) Clamp to realistic bounds to avoid weird plots + v_new = float(np.clip(v_new, 2.5, 4.25)) + + return v_new + + +# ====================================================== +# INPUT CURRENT PROFILE (synthetic) +# ====================================================== +time = np.arange(N) * dt +I = np.zeros(N) + +I[20:80] = 5.0 +I[100:150] = 2.5 +I[170:220] = -3.0 +I[250:280] = 6.0 + +# ====================================================== +# LOG ARRAYS +# ====================================================== +SOC_log = np.zeros(N) +V_log = np.zeros(N) + +# ====================================================== +# INITIAL STATE +# ====================================================== +state = VehicleState( + soc=1.0, + voltage=4.2, # reasonable initial guess +) + +SOC_log[0] = state.soc +V_log[0] = state.voltage + +# ====================================================== +# SIMULATION LOOP +# ====================================================== +for t in range(1, N): + prev_I = I[t - 1] + + # --- SOC UPDATE (Coulomb counting) --- + state.soc -= (prev_I * dt / Q_rated) + state.soc = float(np.clip(state.soc, 0.0, 1.0)) + + # --- VOLTAGE UPDATE (TEMPLATE CALL) --- + state.voltage = update_voltage_template(prev_I, state) + + # --- LOG --- + SOC_log[t] = state.soc + V_log[t] = state.voltage + +# ====================================================== +# PLOTS +# ====================================================== +plt.figure(figsize=(10, 6)) + +plt.subplot(3, 1, 1) +plt.plot(time, I, label="Current [A]") +plt.ylabel("Current (A)") +plt.grid(True) +plt.legend() + +plt.subplot(3, 1, 2) +plt.plot(time, V_log, label="Voltage [V]", color="tab:red") +plt.ylabel("Voltage (V)") +plt.grid(True) +plt.legend() + +plt.subplot(3, 1, 3) +plt.plot(time, SOC_log * 100, label="SOC [%]", color="tab:green") +plt.xlabel("Time (s)") +plt.ylabel("SOC (%)") +plt.grid(True) +plt.legend() + +plt.tight_layout() +plt.show() diff --git a/Docs/Modeling/ChemistryBatteryModel.md b/Docs/Modeling/ChemistryBatteryModel.md index 73b76af..f6f0f7d 100644 --- a/Docs/Modeling/ChemistryBatteryModel.md +++ b/Docs/Modeling/ChemistryBatteryModel.md @@ -34,4 +34,18 @@ $$V = V_0*C_4 + C_2\frac{RT}{F}*ln(\frac{C_1(SOC-0.1^3) + C_3}{1-(SOC-0.1^3)})$$ 1. $C_1$ accounts for my first assumption. It is likely that the anode and cathode volume are not the same and this allows that to be a parameter. 1. $C_2$ accounts for any error in temperature dependence, to best fit the discharge curve slope, and for any errors in by assumption abotu $n$. 1. $C_3$ is part of the $0.001M$ correction and helps shift it in the correct direction. -1. $C_4$ accounts for errors in my assumption that the nominal voltage is equal to $V_0$ and generally just allows the model to fit that value rather than take it as an input. \ No newline at end of file +1. $C_4$ accounts for errors in my assumption that the nominal voltage is equal to $V_0$ and generally just allows the model to fit that value rather than take it as an input. + +### Fix 1: + +So turns out I forgot ln rules and the $C_4$ and $C_1$ are redundant so removing $C_4$... + +### Fix 2: + +Adding $C_4$ to the bottom to do something similar to $C_3$. Also adjusting it so $C_1$ multiplies $C_3$ which makes it less separable but less confusing. + +$$V = V_0 + C_2\frac{RT}{F}*ln(\frac{C_1(SOC-0.1^3 + C_3)}{1-(SOC-0.1^3) + C_4})$$ + +This is the fit to 2.5A data from murata VTC5A cells + +![ChemBatteryModelImage](chemModel2.png) \ No newline at end of file diff --git a/Docs/Modeling/chemModel2.png b/Docs/Modeling/chemModel2.png new file mode 100644 index 0000000..c7c6f4e Binary files /dev/null and b/Docs/Modeling/chemModel2.png differ diff --git a/FullVehicleSim/Electrical/CarData.py b/FullVehicleSim/Electrical/CarData.py deleted file mode 100644 index b9b034e..0000000 --- a/FullVehicleSim/Electrical/CarData.py +++ /dev/null @@ -1,106 +0,0 @@ -import math - -#Track -Gravity:float = 9.8 # m/s^2 -AirDensity:float = 1.293 # kg/m^3 - -#Mass -VehicleMass:int = 0 -DriverMass:int = 0 -ACCMass:int = 0 -TotalMass:float = 300 #kgs should be VehicleMass + DriverMass + ACCMass. We measured vehicle plus acc to be 231 kg and assuming 70 kg driver gives ~300 kg -# FrontRearBalance:float = 0.5 # 0 indicates all weight at rear, 1 indicates all weight at front -LongitudinalDistanceFrontAxleCoM:float = 0.7477 # meters longitudinal distance from front axle to center of gravity -CenterOfMassHeight:float = 0.247718 # meters center of mass height -Wheelbase:float = 1.6547084 # meters wheelbase - -#CenterOfGravity = [] - - - -WheelDiameter:float = 0.46 #From onshape, meters -#Drivetrain -MotorV:int -MotorA:int -DrivingGearTeeth:int = 11 #From onshape -DrivenGearTeeth:int = 40 #From onshape -GearingRatio:float = DrivenGearTeeth/DrivenGearTeeth -#reflected mass inertia?? https://www.motioncontroltips.com/how-do-gearmotors-impact-reflected-mass-inertia-from-the-load/ - - -#Suspension -RollCenter=0 - - -#Aerodynamics -FrontArea:float = 0.693 # m^2, TBD from OptimumLap -DragCoef:float = 0.63 # TBD from OptimumLap - - -#Accumulator -CellsSeries:int = 28 # number of cells in series -CellsParallel:int = 20 # number of cells in parallel - -#Accumulator Cell -InternalImpedance:float = 0.0013 # internal impedance in ohms (1.3mO) - -CurrentVoltageChargeTable = ( - (0 , [(0, 4.2 ), (0.02, 4.16), (0.05, 4.14), (0.1, 4.12), (0.25, 4.07), (0.5, 3.96), (0.75, 3.87), (1, 3.78), (1.25, 3.68), (1.5, 3.61), (1.75, 3.56), (2, 3.49), (2.25, 3.37), (2.35, 3.33), (2.45, 3.27), (2.55, 2.91), (2.6, 2.5 )]), - (0.2, [(0, 4.17), (0.02, 4.15), (0.05, 4.13), (0.1, 4.11), (0.25, 4.06), (0.5, 3.95), (0.75, 3.86), (1, 3.77), (1.25, 3.67), (1.5, 3.6 ), (1.75, 3.55), (2, 3.48), (2.25, 3.36), (2.35, 3.32), (2.45, 3.26), (2.55, 2.9 ), (2.6, 2.5 )]), - (1 , [(0, 4.15), (0.02, 4.12), (0.05, 4.1 ), (0.1, 4.08), (0.25, 4.03), (0.5, 3.91), (0.75, 3.82), (1, 3.73), (1.25, 3.64), (1.5, 3.57), (1.75, 3.52), (2, 3.45), (2.25, 3.32), (2.35, 3.27), (2.45, 3.16), (2.55, 2.7 )]), - (2 , [(0, 4.14), (0.02, 4.1 ), (0.05, 4.07), (0.1, 4.05), (0.25, 3.99), (0.5, 3.87), (0.75, 3.78), (1, 3.7 ), (1.25, 3.6 ), (1.5, 3.53), (1.75, 3.47), (2, 3.4 ), (2.25, 3.28), (2.35, 3.21), (2.45, 3.01), (2.55, 2.51)]), - (3 , [(0, 4.12), (0.02, 4.07), (0.05, 4.05), (0.1, 4.02), (0.25, 3.96), (0.5, 3.84), (0.75, 3.75), (1, 3.67), (1.25, 3.57), (1.5, 3.5 ), (1.75, 3.43), (2, 3.36), (2.25, 3.24), (2.35, 3.15), (2.45, 2.91), (2.55, 2.5 )]), - (5 , [(0, 4.09), (0.02, 4.03), (0.05, 4.0 ), (0.1, 3.97), (0.25, 3.9 ), (0.5, 3.79), (0.75, 3.7 ), (1, 3.61), (1.25, 3.51), (1.5, 3.44), (1.75, 3.38), (2, 3.3 ), (2.25, 3.18), (2.35, 3.09), (2.45, 2.87), (2.55, 2.5 )]), - (10 , [(0, 4.0 ), (0.02, 3.93), (0.05, 3.89), (0.1, 3.85), (0.25, 3.79), (0.5, 3.69), (0.75, 3.59), (1, 3.5 ), (1.25, 3.41), (1.5, 3.34), (1.75, 3.27), (2, 3.19), (2.25, 3.09), (2.35, 3.01), (2.45, 2.87), (2.55, 2.5 )]), - (20 , [(0, 3.82), (0.02, 3.74), (0.05, 3.7 ), (0.1, 3.66), (0.25, 3.58), (0.5, 3.49), (0.75, 3.4 ), (1, 3.32), (1.25, 3.24), (1.5, 3.18), (1.75, 3.12), (2, 3.05), (2.25, 2.95), (2.35, 2.89), (2.45, 2.76)]), - (30.01 , [(0, 3.62), (0.02, 3.55), (0.05, 3.49), (0.1, 3.45), (0.25, 3.36), (0.5, 3.26), (0.75, 3.18), (1, 3.13), (1.25, 3.08), (1.5, 3.03), (1.75, 2.98), (2, 2.92), (2.25, 2.82), (2.35, 2.7 )]) -) # ((current, [(discharge, voltage), ...]), ...) -# 30.01 accounts for floating point errors - -def lookupLerpVoltage(current:float, discharge:float): # current in amps, discharge is Ahs discharged, this is a bilinear interpolation - for i in range(1, len(CurrentVoltageChargeTable)): - highCurrent = CurrentVoltageChargeTable[i] - - if highCurrent[0] > current: - - lowCurrent = CurrentVoltageChargeTable[i-1] - lcVoltage:float = 0 - for j in range(1, len(lowCurrent[1])): - highDischarge = lowCurrent[1][j] - if highDischarge[0] > discharge: - lowDischarge = lowCurrent[1][j-1] - lcVoltage = lowDischarge[1] + (highDischarge[1]-lowDischarge[1])*((discharge-lowDischarge[0])/(highDischarge[0]-lowDischarge[0])) - break - - hcVoltage:float = 0 - for k in range(1, len(highCurrent[1])): - highDischarge = highCurrent[1][k] - if highDischarge[0] > discharge: - lowDischarge = highCurrent[1][k-1] - hcVoltage = lowDischarge[1] + (highDischarge[1]-lowDischarge[1])*((discharge-lowDischarge[0])/(highDischarge[0]-lowDischarge[0])) - break - - return lcVoltage + (hcVoltage-lcVoltage)*((current-lowCurrent[0])/(highCurrent[0]-lowCurrent[0])) - return MinVoltage - - -MinVoltage:float = 2.5 # volts -MaxDischargeCurrent:float = 30 # amps -Capacity:float = 2.5 # Ah - -#motor -RPMtorqueTable = [(0, 178.4), (4120, 178.4), (7500, 82)] # data is [(RPM, Nm), ...] -def lookupLerpTorque(rpm:float): - if (rpm < 0): - return 0 - for i in range(1, len(RPMtorqueTable)): - tup = RPMtorqueTable[i] - if tup[0] == rpm: - return tup[1] - elif tup[0] > rpm: - prevTup = RPMtorqueTable[i-1] - return (prevTup[1] + ((tup[1]-prevTup[1]) * ((rpm-prevTup[0])/(tup[0]-prevTup[0])))) - return 0 - -#tires -LongitudinalFrictionCoef = 1.5 # coef of static friction \ No newline at end of file diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/F1_project1.py b/FullVehicleSim/Electrical/HisteresisCellModel/F1_project1.py new file mode 100644 index 0000000..f996b13 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/F1_project1.py @@ -0,0 +1 @@ +prin \ No newline at end of file diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining.py b/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining.py new file mode 100644 index 0000000..a2b5f1e --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining.py @@ -0,0 +1,77 @@ +## When we train from data, use this + +from scipy.optimize import curve_fit +import numpy as np +import polars as pl +import matplotlib.pyplot as plt +from Data.FSLib.IntegralsAndDerivatives import integrate_with_Scipy_tCol +from Data.FSLib.AnalysisFunctions import simpleTimeCol + +# df = pl.read_csv("C:/Projects/FormulaSlug/fs-data/FS-3/voltageTableVTC5A.csv") + +temps = [] + +for i in range(5): + for j in range(6): + temps.append(f"ACC_SEG{i}_TEMPS_CELL{j}") + +df = pl.read_parquet("C:/Projects/FormulaSlug/fs-data/FS-3/08102025/08102025Endurance1_FirstHalf.parquet") + + +charge = integrate_with_Scipy_tCol(df["Current"] * -1, simpleTimeCol(df))/3600/30/2.6 # Coulombs --> Ah, 30 cells --> 1 cell, 2.6Ah per cell + +df = df.with_columns( + pl.col("ACC_POWER_PACK_VOLTAGE").alias("Voltage"), + df.select(temps).mean_horizontal().alias("Temperature"), + pl.col("SME_TEMP_BusCurrent").alias("Current") +) + + +# plt.scatter(df["Charge"], df["Voltage"],c=df["Current"], label="Current") +# plt.xlabel("Charge (Ah)") +# plt.legend() +# plt.show() + +dt = 0.01 +kernel_duration = 10.0 +kernel_size = int(kernel_duration / dt) +t = np.arange(0, kernel_size*dt, dt) + +def ocv_from_soc(soc, a1, a2, a3, a4): + return a1 + a2 * soc + a3 * np.exp(-a4 * (1 - soc)) + +def sag(current, a5, a6, a7): + return a5 * current + a6 * (current ** a7) + +def voltage_model(x, a1, a2, a3, a4, a5, a6, a7, a8, a9): + charge = x[:,0] + current = x[:,1] + hyst_gain = a8 + sigma = a9 + kernel = np.exp(-(t**2) / (2 * sigma**2)) + kernel /= np.sum(kernel) + prev_curr = np.zeros((charge.shape[0], kernel_size)) + for i in range(charge.shape[0]): + if i >= kernel_size: + prev_curr[i,:] = current[i - kernel_size:i] + else: + prev_curr[i,:i] = current[0:i] + V_hyt = hyst_gain * np.sum(prev_curr * t, axis=1) + V_ocv = ocv_from_soc(charge / 2.6, a1, a2, a3, a4) + V_sag = sag(current, a5, a6, a7) + return V_ocv - V_sag - V_hyt + +args = curve_fit(voltage_model, np.column_stack((df["Charge"], df["Current"])), df["Voltage"], p0=[3.0, 0.9, 0.25, 12.0, 0.02, 0.004, 1.3, 0.015, 3.0], maxfev=10000) +args[0] + +a1, a2, a3, a4, a5, a6, a7, a8, a9 = args[0] +plt.figure(figsize=(10,6)) +plt.scatter(df["Charge"], df["Voltage"], c='blue', label="Measured Voltage", alpha=0.5) +predicted_voltage = voltage_model(np.column_stack((df["Charge"], df["Current"])), a1, a2, a3, a4, a5, a6, a7, a8, a9) +plt.scatter(df["Charge"], predicted_voltage, c='red', label="Fitted Voltage", alpha=0.5) +plt.xlabel("Charge (Ah)") +plt.ylabel("Voltage (V)") +plt.legend() +plt.show() + + \ No newline at end of file diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py b/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py new file mode 100644 index 0000000..8b5fa51 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py @@ -0,0 +1,112 @@ +## When we train from data, use this + +from scipy.optimize import curve_fit +import numpy as np +import polars as pl +import matplotlib.pyplot as plt +# from Data.FSLib.IntegralsAndDerivatives import integrate_with_Scipy_tCol +# from Data.FSLib.AnalysisFunctions import simpleTimeCol + +df = pl.read_csv("C:/Projects/FormulaSlug/fs-data/FS-3/voltageTableVTC5A.csv") +dfLowCurr = df.filter(pl.col("Current") < 1).filter(pl.col("Voltage") > 2.5) + +df.head + +R = 8.314 +T = 293 +F = 96485.3329 +n = 1 +a = 3 +cmc = 2.6 # Cell Max Capacity (Ah) +V0 = 3.7 + +# temps = [] + +# for i in range(5): +# for j in range(6): +# temps.append(f"ACC_SEG{i}_TEMPS_CELL{j}") + +# df = pl.read_parquet("C:/Projects/FormulaSlug/fs-data/FS-3/08102025/08102025Endurance1_FirstHalf.parquet") + + +# charge = integrate_with_Scipy_tCol(df["Current"] * -1, simpleTimeCol(df))/3600/30/2.6 # Coulombs --> Ah, 30 cells --> 1 cell, 2.6Ah per cell + +# df = df.with_columns( +# pl.col("ACC_POWER_PACK_VOLTAGE").alias("Voltage"), +# df.select(temps).mean_horizontal().alias("Temperature"), +# pl.col("SME_TEMP_BusCurrent").alias("Current") +# ) + + +# plt.scatter(df["Charge"], df["Voltage"],c=df["Current"], label="Current") +# plt.xlabel("Charge (Ah)") +# plt.legend() +# plt.show() + +dt = 0.01 +kernel_duration = 10.0 +kernel_size = int(kernel_duration / dt) +t = np.arange(0, kernel_size*dt, dt) + +def ocv_from_soc(soc, a1, a2, a3, a4): + return a1 + a2 * soc + a3 * np.exp(-a4 * (1 - soc)) + +def sag(current, a5, a6, a7): + return a5 * current + a6 * (current ** a7) + +def voltage_SOC_model(x, c1, c2, c3, c4): + return V0 - c2*R*T/(n*F)*np.log((c1*((x/cmc) - (0.1**a) + c3))/(1 - (x/cmc) + (0.1**a) + c4)) + +def voltage_model(x, a1, a2, a3, a4, a5, a6, a7, a8, a9): + charge = x[:,0] + current = x[:,1] + hyst_gain = a8 + sigma = a9 + kernel = np.exp(-(t**2) / (2 * sigma**2)) + kernel /= np.sum(kernel) + prev_curr = np.zeros((charge.shape[0], kernel_size)) + for i in range(charge.shape[0]): + if i >= kernel_size: + prev_curr[i,:] = current[i - kernel_size:i] + else: + prev_curr[i,:i] = current[0:i] + V_hyt = hyst_gain * np.sum(prev_curr * t, axis=1) + V_ocv = ocv_from_soc(charge / cmc, a1, a2, a3, a4) + V_sag = sag(current, a5, a6, a7) + return V_ocv - V_sag - V_hyt + +args = curve_fit(voltage_model, np.column_stack((df["Charge"], df["Current"])), df["Voltage"], p0=[3.0, 0.9, 0.25, 12.0, 0.02, 0.004, 1.3, 0.015, 3.0], maxfev=10000) +args[0] + +a1, a2, a3, a4, a5, a6, a7, a8, a9 = args[0] +plt.figure(figsize=(10,6)) +plt.scatter(df["Charge"], df["Voltage"], c='blue', label="Measured Voltage", alpha=0.5) +predicted_voltage = voltage_model(np.column_stack((df["Charge"], df["Current"])), a1, a2, a3, a4, a5, a6, a7, a8, a9) +plt.scatter(df["Charge"], predicted_voltage, c='red', label="Fitted Voltage", alpha=0.5) +plt.xlabel("Charge (Ah)") +plt.ylabel("Voltage (V)") +plt.legend() +plt.show() + +dfLowCurr1 = df.filter(pl.col("Current") < 1).filter(pl.col("Voltage") > 2.5) +# dfLowCurr2 = df.filter(pl.col("Current") < 3) + + +args1 = curve_fit(voltage_SOC_model, dfLowCurr1["Charge"], dfLowCurr1["Voltage"], p0=[1.5, 5.36, 0.0019, 3.7]) +# args2 = curve_fit(voltage_SOC_model, dfLowCurr2["Charge"], dfLowCurr2["Voltage"], p0=[1.5, 5.36, 0.0019, 3.7]) + +c1, c2, c3, c4 = args1[0] +# c5, c6, c7, c8 = args2[0] +plt.figure(figsize=(10,6)) +plt.scatter(dfLowCurr1["Charge"], dfLowCurr1["Voltage"], c='blue', label="Measured Voltage", alpha=0.5) +predicted_voltage_soc1 = voltage_SOC_model(np.arange(0, 2.6, 0.01), c1, c2, c3, c4) +# predicted_voltage_soc2 = voltage_SOC_model(np.arange(0, 2.6, 0.01), c5, c6, c7, c8) +plt.scatter(np.arange(0, 2.6, 0.01), predicted_voltage_soc1, c='red', label="Fitted Voltage 1", alpha=0.5) +# plt.scatter(np.arange(0, 2.6, 0.01), predicted_voltage_soc2, c='green', label="Fitted Voltage 2", alpha=0.5) +plt.xlabel("Charge (Ah)") +plt.ylabel("Voltage (V)") +plt.legend() +plt.show() + +args1[0] +# args2[0] \ No newline at end of file diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTrainingDischargeCurvepy b/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTrainingDischargeCurvepy new file mode 100644 index 0000000..3fa2f66 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTrainingDischargeCurvepy @@ -0,0 +1,54 @@ +## When we train from discharge curve, use this + +from scipy.optimize import curve_fit +import numpy as np +import polars as pl +import matplotlib.pyplot as plt + +df = pl.read_csv("C:/Projects/FormulaSlug/fs-data/FS-3/voltageTableVTC5A.csv") + +df = df.filter(pl.col("Voltage") > 2.5) + +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.array([charge for _ in range(kernel_size)]).T + V_hyt = hyst_gain * np.sum(prev_curr * t, axis=1) + V_ocv = ocv_from_soc((2.6 - 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/electrical_model_1.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_1.py new file mode 100644 index 0000000..529a211 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_1.py @@ -0,0 +1,37 @@ +import time +import random +import matplotlib.pyplot as plt +from collections import deque + +# Store (time, current) values +time_window = 10 # seconds +times = deque() +currents = deque() + +plt.ion() # interactive mode +fig, ax = plt.subplots() + +start_time = time.time() + +while True: + now = time.time() - start_time + current_value = random.uniform(0, 10) # replace with real sensor current + + # Add new data + times.append(now) + currents.append(current_value) + + # Remove old data (>10 seconds) + while times and times[0] < now - time_window: + times.popleft() + currents.popleft() + + # Update plot + ax.clear() + ax.plot(times, currents, color='blue') + ax.set_xlim(max(0, now - time_window), now) + ax.set_xlabel("Time (s)") + ax.set_ylabel("Current (A)") + ax.set_title("Current over last 10 seconds") + + plt.pause(0.1) diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_2.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_2.py new file mode 100644 index 0000000..eb5a25a --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_2.py @@ -0,0 +1,138 @@ +import numpy as np +import matplotlib.pyplot as plt + +# ===================================================== +# Accumulator Voltage Model (Stateful + Hysteresis) +# ===================================================== + +class AccumulatorVoltageModel: + def __init__(self, dt=1.0): + self.dt = dt # time step [s] + self.capacity_Ah = 2.6 # cell capacity [Ah] + self.SOC = 1.0 # initial SOC + + # 10-second current history buffer + self.I_hist = np.zeros(10) + + # Gaussian kernel for hysteresis + t = np.arange(10) + sigma = 3.0 + self.kernel = np.exp(-(t**2) / (2 * sigma**2)) + self.kernel /= np.sum(self.kernel) + + self.hyst_gain = 0.015 # hysteresis strength + + # ------------------------------------------------- + # Open Circuit Voltage from SOC (datasheet-like) + # ------------------------------------------------- + def ocv_from_soc(self, soc): + return 3.0 + 0.9*soc + 0.25*np.exp(-12*(1 - soc)) + + + 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) + + # Slide current history window + self.I_hist[:-1] = self.I_hist[1:] + self.I_hist[-1] = current + + # Hysteresis voltage term + V_hyst = self.hyst_gain * np.sum(self.I_hist * self.kernel) + + # Voltage calculation + V = ( + self.ocv_from_soc(self.SOC) + - self.sag(current) * (1 - self.SOC) + - V_hyst + ) + + return V + + +model = AccumulatorVoltageModel(dt=1.0) + +time = [] +voltage = [] +soc = [] +current_log = [] + +current_profile = ( + [5]*20 + # cruise + [20]*10 + # acceleration + [10]*20 + # steady drive + [0]*10 # regen / idle +) + +for t, I in enumerate(current_profile): + V = model.step(I) + time.append(t) + voltage.append(V) + soc.append(model.SOC) + current_log.append(I) + +# ===================================================== +# Plots +# ===================================================== + +plt.figure(figsize=(14,10)) + +# Voltage vs Time +plt.subplot(2,2,1) +plt.plot(time, voltage) +plt.xlabel("Time [s]") +plt.ylabel("Voltage [V]") +plt.title("Accumulator Voltage vs Time") +plt.grid(True) + +# SOC vs Time +plt.subplot(2,2,2) +plt.plot(time, soc) +plt.xlabel("Time [s]") +plt.ylabel("SOC") +plt.title("State of Charge vs Time") +plt.grid(True) + +# Current Profile +plt.subplot(2,2,3) +plt.plot(time, current_log) +plt.xlabel("Time [s]") +plt.ylabel("Current [A]") +plt.title("Vehicle Current Profile") +plt.grid(True) + +# Voltage vs SOC +plt.subplot(2,2,4) +plt.plot(soc, voltage) +plt.xlabel("SOC") +plt.ylabel("Voltage [V]") +plt.title("Voltage vs SOC") +plt.grid(True) + +plt.tight_layout() +plt.show() + +# ===================================================== +# Hysteresis Demonstration +# ===================================================== + +model = AccumulatorVoltageModel() + +# High-current history +for _ in range(10): + model.step(20) +V_high_hist = model.step(5) + +# Low-current history +model = AccumulatorVoltageModel() +for _ in range(10): + model.step(0) +V_low_hist = model.step(5) + +print("Voltage with high-current history:", round(V_high_hist, 3), "V") +print("Voltage with low-current history :", round(V_low_hist, 3), "V") diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_3.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_3.py new file mode 100644 index 0000000..76570c3 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_3.py @@ -0,0 +1,137 @@ +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 + +class AccumulatorVoltageModel: + def __init__(self, dt=1): + self.dt = dt + self.capacity_Ah = 2.8 + self.SOC = 1.0 + + self.I_hist = np.zeros(10) + + t = np.arange(10) + sigma = 0.2 # a9 + self.kernel = np.exp(-(t**2) / (2 * sigma**2)) + self.kernel /= np.sum(self.kernel) + + self.hyst_gain = 0.0 # a8 + + def ocv_from_soc(self, 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.0 * current + + 0.0 * (abs(current) ** 1.0) + ) + + + def step(self, current): + self.SOC -= (current / 30 * self.dt) / (3600 * self.capacity_Ah) + self.SOC = np.clip(self.SOC, 0.0, 1.0) + + self.I_hist[:-1] = self.I_hist[1:] + self.I_hist[-1] = current + + V_hyst = self.hyst_gain * np.dot(self.I_hist, self.kernel) + + pack_voltage = ( + self.ocv_from_soc(self.SOC) + - self.sag(current) * (1 - self.SOC) + - V_hyst + ) + + cell_voltage = pack_voltage / 30 + + return cell_voltage + +model = AccumulatorVoltageModel() + +voltage_log = [] +soc_log = [] +I_hist_log = [] + +for I in current_profile: + 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) + + +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) + +plt.title("State of Charge") +plt.xlabel("Time step") +plt.ylabel("SOC") +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() diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_5.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_5.py new file mode 100644 index 0000000..933550c --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_5.py @@ -0,0 +1,215 @@ +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 = 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") + + +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 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/sims-battery_voltage.py b/FullVehicleSim/Electrical/HisteresisCellModel/sims-battery_voltage.py new file mode 100644 index 0000000..151bd63 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/sims-battery_voltage.py @@ -0,0 +1,132 @@ +import numpy as np +import matplotlib.pyplot as plt + +# ------------------------------------------------------- +# PARAMETERS (tuned for Murata VTC5A-like behavior) +# ------------------------------------------------------- +Q_rated = 2.5 * 3600 # cell capacity [Coulombs] (2.5Ah) +R0 = 0.015 # Ohmic resistance [Ohms] +R1 = 0.02 # RC branch resistance [Ohms] +C1 = 200.0 # RC branch capacitance [Farads] +dt = 1.0 # time step [seconds] +N = 300 # number of time steps + +lambda_hyst = 0.98 # hysteresis memory decay +gamma_hyst = 0.03 # hysteresis voltage scaling [V] +T_ref = 25.0 # reference temp [°C] +alpha_R = 0.004 # resistance temp coefficient [/°C] +beta_V = 0.0008 # OCV temperature coefficient [V/°C] +m, c, hA = 0.045, 1000.0, 0.6 # mass [kg], specific heat [J/kgK], cooling term [W/K] + +# ------------------------------------------------------- +# GAUSSIAN KERNEL FOR HYSTERESIS +# ------------------------------------------------------- +kernel_size = 9 # must be odd +sigma = 2.0 + +x = np.linspace(-3, 3, kernel_size) +gaussian_kernel = np.exp(-(x**2) / (2 * sigma**2)) +gaussian_kernel /= gaussian_kernel.sum() + +# ------------------------------------------------------- +# HELPER FUNCTIONS +# ------------------------------------------------------- +def f_OCV(SOC): + """Open-circuit voltage vs SOC curve (simple polynomial fit)""" + return 3.0 + 1.2*SOC - 0.3*(SOC**2) + 0.1*(SOC**3) + +def R_int(T): + """Internal resistance vs temperature""" + return R0 * (1 + alpha_R * (T - T_ref)) + +# ------------------------------------------------------- +# TIME-SERIES INPUTS (synthetic) +# ------------------------------------------------------- +time = np.arange(N) * dt +I = np.zeros(N) +I[20:80] = 5.0 +I[100:150] = 2.5 +I[170:220] = -3.0 +I[250:280] = 6.0 + +# ------------------------------------------------------- +# INITIAL CONDITIONS +# ------------------------------------------------------- +SOC = np.zeros(N) +V_RC = np.zeros(N) +H = np.zeros(N) +T = np.zeros(N) +V_model = np.zeros(N) + +SOC[0] = 1.0 +T[0] = 25.0 +V_RC[0] = 0.0 +H[0] = 0.0 + +# ------------------------------------------------------- +# SIMULATION LOOP +# ------------------------------------------------------- +for t in range(1, N): + + # SOC update + SOC[t] = SOC[t-1] - (I[t-1] * dt / Q_rated) + SOC[t] = np.clip(SOC[t], 0, 1) + + # OCV lookup + V_OCV = f_OCV(SOC[t]) + + # RC polarization voltage + exp_factor = np.exp(-dt / (R1 * C1)) + V_RC[t] = V_RC[t-1] * exp_factor + R1 * (1 - exp_factor) * I[t-1] + + # ---------------------------- + # HYSTERESIS (Gaussian + memory) + # ---------------------------- + I_window = I[max(0, t - kernel_size + 1):t + 1] + k_window = gaussian_kernel[-len(I_window):] # trim kernel for start of time series + + I_gauss = np.sum(I_window * k_window) + + H[t] = ( + lambda_hyst * H[t-1] + + (1 - lambda_hyst) * I[t-1] + + 0.2 * I_gauss # <-- you can tune this weight + ) + + # Thermal update + T[t] = T[t-1] + (dt / (m * c)) * (I[t-1]**2 * R_int(T[t-1]) - hA * (T[t-1] - T_ref)) + + # Temperature corrections + R0_eff = R0 * (1 + alpha_R * (T[t] - T_ref)) + V_OCV_T = V_OCV + beta_V * (T[t] - T_ref) + + # Terminal voltage + V_model[t] = V_OCV_T - I[t] * R0_eff - V_RC[t] + gamma_hyst * H[t] + +# ------------------------------------------------------- +# PLOTS +# ------------------------------------------------------- +plt.figure(figsize=(10, 6)) + +plt.subplot(3,1,1) +plt.plot(time, I, label='Current [A]') +plt.ylabel('Current (A)') +plt.grid(True) +plt.legend() + +plt.subplot(3,1,2) +plt.plot(time, V_model, label='Modeled Voltage [V]', color='tab:red') +plt.ylabel('Voltage (V)') +plt.grid(True) +plt.legend() + +plt.subplot(3,1,3) +plt.plot(time, SOC*100, label='State of Charge [%]', color='tab:green') +plt.xlabel('Time (s)') +plt.ylabel('SOC (%)') +plt.grid(True) +plt.legend() + +plt.tight_layout() +plt.show() + 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/granola-archive.py b/FullVehicleSim/Electrical/granola-archive.py deleted file mode 100644 index 839740e..0000000 --- a/FullVehicleSim/Electrical/granola-archive.py +++ /dev/null @@ -1,99 +0,0 @@ -import math -from Powertrain import CarData -import copy - -class MBS: - time:float # seconds since the sim started, at start of step - deltaTime:float # seconds, time until next step - velocity:float # m/s, at start of step - batteryVoltage:float # volts, at start of step - motorSpeed:float # rpm, at start of step - throttleDemand:float # 0 to 1, at start of step - motorTorque:float # Nm, at start of step - motorPower:float # kW, at start of step - currentDraw:float # Amperes, at start of step - tractiveForce:float # Newtons, positive is forward, at start of step - netForce:float # Newtons, positive is forward, at start of step - totalDischarged:float # Ah, at start of step - deltaDischarged:float # Ah, discharged until next step - - - def __init__(self, throttle, initDischarged, maxPower): - self.throttle = throttle # percent of max torque to demand - self.initDischarged = initDischarged # how many Ahs have already been expended on the car, 2.8 yeilds a voltage of about 115 - self.maxPower = maxPower # 80kW max - #return self - - def calculateStep(self, prevStep, timeDelta, tractionForce): - self.frictionForce = tractionForce - self.deltaTime = timeDelta - self.time = prevStep.time + self.deltaTime - self.velocity = prevStep.velocity + self.deltaTime * prevStep.netForce / CarData.TotalMass - self.totalDischarged = prevStep.totalDischarged + prevStep.deltaDischarged - self.batteryVoltage = CarData.lookupLerpVoltage(prevStep.currentDraw / CarData.CellsParallel, self.totalDischarged / CarData.CellsParallel) * CarData.CellsSeries - self.motorSpeed = self.calculateMotorRPMfromVelocity(self.velocity) - self.throttleDemand = self.throttle - self.torqueDemand = self.throttleDemand * CarData.lookupLerpTorque(self.motorSpeed) - self.motorTorque = self.calculateRealTorque(self.torqueDemand, self.batteryVoltage, self.motorSpeed) - self.motorPower = self.calculateMotorPower(self.motorTorque, self.motorSpeed) - self.currentDraw = self.calculateCurrentDraw(self.motorPower, self.batteryVoltage) - self.deltaDischarged = self.calculateDischargeDelta(self.currentDraw, self.deltaTime) - - self.tractiveForce = self.calculateTractiveForce(self.motorTorque) # assuming two driving wheels - return copy.deepcopy(self) - - def firstStep(self, initialDischarged, timeDelta, slipRatio): - self.frictionForce = 0 - self.slipRatio = slipRatio - self.time = 0 - self.deltaTime = timeDelta - self.velocity = 0 - self.totalDischarged = initialDischarged - self.batteryVoltage = CarData.lookupLerpVoltage(0, self.totalDischarged / CarData.CellsParallel) * CarData.CellsSeries - self.motorSpeed = 0 - self.throttleDemand = self.throttle - self.torqueDemand = self.throttleDemand * CarData.lookupLerpTorque(self.motorSpeed) - self.motorTorque = self.calculateRealTorque(self.torqueDemand, self.batteryVoltage, self.motorSpeed) - self.motorPower = self.calculateMotorPower(self.motorTorque, self.motorSpeed) - self.currentDraw = self.calculateCurrentDraw(self.motorPower, self.batteryVoltage) - self.deltaDischarged = self.calculateDischargeDelta(self.currentDraw, self.deltaTime) - - self.tractiveForce = self.calculateTractiveForce(self.motorTorque) - - print("GOT TO CREATION OF FIRST STEP") - return self - - - def calculateMaxTireFriction(self): - return ((self.frictionForce * CarData.LongitudinalDistanceFrontAxleCoM / CarData.Wheelbase) / (1 - (CarData.CenterOfMassHeight / CarData.Wheelbase)* (self.frictionForce/9.8/CarData.TotalMass))) - - def calculateTractiveForce(self, motorTorque:float): - axleTorque = motorTorque * CarData.DrivenGearTeeth / CarData.DrivingGearTeeth # torque for both rear wheels - wheelForce = axleTorque / (CarData.WheelDiameter/2) # T = Fr so F = T/r - - return wheelForce - - def calculateVelocityfromMotorRPM(self, motorRPM:float): - return (motorRPM * CarData.DrivingGearTeeth * math.pi * CarData.WheelDiameter)/(60*CarData.DrivenGearTeeth) # turns RPM into m/s - - def calculateMotorRPMfromVelocity(self, velocity:float): - return (velocity*60*CarData.DrivenGearTeeth)/(CarData.DrivingGearTeeth * math.pi * CarData.WheelDiameter) # turns m/s into RPM - - def calculateMotorPower(self, motorTorque, motorRPM): - return motorTorque*motorRPM/9550 - - def calculateCurrentDraw(self, motorPower, motorVoltage): - return 1000.0*motorPower/motorVoltage # times 1000 to account for power being kW - - def calculateRealTorque(self, torqueDemand, batVoltage, motorRPM): - currentTorqueLim = torqueDemand - if (motorRPM != 0): - currentLim = min(CarData.MaxDischargeCurrent*CarData.CellsParallel, self.maxPower/batVoltage*1000) - currentTorqueLim = currentLim * batVoltage * 9550 / (motorRPM * 1000) - - tractionTorqueLim = self.calculateMaxTireFriction() * (CarData.WheelDiameter/2) * CarData.DrivingGearTeeth / CarData.DrivenGearTeeth - - return min(torqueDemand, currentTorqueLim, tractionTorqueLim) - - def calculateDischargeDelta(self, currentDraw, timeDelta): - return currentDraw * timeDelta / 3600 # divide by 3600 to convert from amp seconds to amp hours \ No newline at end of file diff --git a/FullVehicleSim/Electrical/granola-archives.py b/FullVehicleSim/Electrical/granola-archives.py deleted file mode 100644 index 90a6c2c..0000000 --- a/FullVehicleSim/Electrical/granola-archives.py +++ /dev/null @@ -1,136 +0,0 @@ -import sys -import os -sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..', 'TireModel'))) - -import math -# import dumpling as tire -tire = 0 -import copy - -class MBS: - def __init__(self, params=None): - if not params: - return False - - self.timeStepDelta = params["timeStepDelta"] - self.totalDischarge = params["initDischarge"] - self.velocity = params["initVelocity"] - - self.motorSpeed = 0 - self.currStep = 0 - self.throttleDemand = 1 - self.motorTorque = 0 - self.motorPower = 0 - self.currentDraw = 0 - self.deltaDischarged = 0 - self.tractiveForce = 0 - - self.cellsInParallel = 20 - self.cellsInSeries = 28 - self.maxPower = 80 - self.minVoltage = 2.5 - self.maxDischargeCurrent = 30 - self.drivenGearTeeth = 40 - self.drivingGearTeeth = 11 - self.wheelDiameter = 0.46 #meters - - - self.rearLeftTire = tire.Tire(150 * 9.8, 0.15, 0, self.velocity, 85, 45, params["Constants"], params["Mechanical-Parameters"], params["Magic-Parameters"]) - self.rearTireTraction = self.rearLeftTire.getLongForce() #1.5 * 9.8 * 150 #Temporary guess - - def calculateStep(self, throttle, currVelocity): - self.currStep += 1 - self.velocity = currVelocity - self.totalDischarge += self.deltaDischarged - self.batteryVoltage = self.lookupLerpVoltage(self.currentDraw / self.cellsInParallel, self.totalDischarge / self.cellsInParallel) * self.cellsInSeries - self.motorRPM = self.calculateMotorRPMfromVelocity() - self.throttleDemand = throttle - self.torqueDemand = self.throttleDemand * self.lookupLerpTorque(self.motorRPM) - self.motorTorque = self.calculateRealTorque() - self.motorPower = self.calculateMotorPower() - self.currentDraw = self.calculateCurrentDraw() - self.deltaDischarged = self.calculateDischargeDelta() - - - self.tractiveForce = self.calculateTractiveForce() - - return copy.deepcopy(self) - - def lookupLerpVoltage(self, current, discharge): # Stolen straight from Cole's acceleration model - - CurrentVoltageChargeTable = ( - (0 , [(0, 4.2 ), (0.02, 4.16), (0.05, 4.14), (0.1, 4.12), (0.25, 4.07), (0.5, 3.96), (0.75, 3.87), (1, 3.78), (1.25, 3.68), (1.5, 3.61), (1.75, 3.56), (2, 3.49), (2.25, 3.37), (2.35, 3.33), (2.45, 3.27), (2.55, 2.91), (2.6, 2.5 )]), - (0.2, [(0, 4.17), (0.02, 4.15), (0.05, 4.13), (0.1, 4.11), (0.25, 4.06), (0.5, 3.95), (0.75, 3.86), (1, 3.77), (1.25, 3.67), (1.5, 3.6 ), (1.75, 3.55), (2, 3.48), (2.25, 3.36), (2.35, 3.32), (2.45, 3.26), (2.55, 2.9 ), (2.6, 2.5 )]), - (1 , [(0, 4.15), (0.02, 4.12), (0.05, 4.1 ), (0.1, 4.08), (0.25, 4.03), (0.5, 3.91), (0.75, 3.82), (1, 3.73), (1.25, 3.64), (1.5, 3.57), (1.75, 3.52), (2, 3.45), (2.25, 3.32), (2.35, 3.27), (2.45, 3.16), (2.55, 2.7 )]), - (2 , [(0, 4.14), (0.02, 4.1 ), (0.05, 4.07), (0.1, 4.05), (0.25, 3.99), (0.5, 3.87), (0.75, 3.78), (1, 3.7 ), (1.25, 3.6 ), (1.5, 3.53), (1.75, 3.47), (2, 3.4 ), (2.25, 3.28), (2.35, 3.21), (2.45, 3.01), (2.55, 2.51)]), - (3 , [(0, 4.12), (0.02, 4.07), (0.05, 4.05), (0.1, 4.02), (0.25, 3.96), (0.5, 3.84), (0.75, 3.75), (1, 3.67), (1.25, 3.57), (1.5, 3.5 ), (1.75, 3.43), (2, 3.36), (2.25, 3.24), (2.35, 3.15), (2.45, 2.91), (2.55, 2.5 )]), - (5 , [(0, 4.09), (0.02, 4.03), (0.05, 4.0 ), (0.1, 3.97), (0.25, 3.9 ), (0.5, 3.79), (0.75, 3.7 ), (1, 3.61), (1.25, 3.51), (1.5, 3.44), (1.75, 3.38), (2, 3.3 ), (2.25, 3.18), (2.35, 3.09), (2.45, 2.87), (2.55, 2.5 )]), - (10 , [(0, 4.0 ), (0.02, 3.93), (0.05, 3.89), (0.1, 3.85), (0.25, 3.79), (0.5, 3.69), (0.75, 3.59), (1, 3.5 ), (1.25, 3.41), (1.5, 3.34), (1.75, 3.27), (2, 3.19), (2.25, 3.09), (2.35, 3.01), (2.45, 2.87), (2.55, 2.5 )]), - (20 , [(0, 3.82), (0.02, 3.74), (0.05, 3.7 ), (0.1, 3.66), (0.25, 3.58), (0.5, 3.49), (0.75, 3.4 ), (1, 3.32), (1.25, 3.24), (1.5, 3.18), (1.75, 3.12), (2, 3.05), (2.25, 2.95), (2.35, 2.89), (2.45, 2.76)]), - (30.01 , [(0, 3.62), (0.02, 3.55), (0.05, 3.49), (0.1, 3.45), (0.25, 3.36), (0.5, 3.26), (0.75, 3.18), (1, 3.13), (1.25, 3.08), (1.5, 3.03), (1.75, 2.98), (2, 2.92), (2.25, 2.82), (2.35, 2.7 )]) - ) - - - for i in range(1, len(CurrentVoltageChargeTable)): - highCurrent = CurrentVoltageChargeTable[i] - - if highCurrent[0] > current: - - lowCurrent = CurrentVoltageChargeTable[i-1] - lcVoltage:float = 0 - for j in range(1, len(lowCurrent[1])): - highDischarge = lowCurrent[1][j] - if highDischarge[0] > discharge: - lowDischarge = lowCurrent[1][j-1] - lcVoltage = lowDischarge[1] + (highDischarge[1]-lowDischarge[1])*((discharge-lowDischarge[0])/(highDischarge[0]-lowDischarge[0])) - break - - hcVoltage:float = 0 - for k in range(1, len(highCurrent[1])): - highDischarge = highCurrent[1][k] - if highDischarge[0] > discharge: - lowDischarge = highCurrent[1][k-1] - hcVoltage = lowDischarge[1] + (highDischarge[1]-lowDischarge[1])*((discharge-lowDischarge[0])/(highDischarge[0]-lowDischarge[0])) - break - - return lcVoltage + (hcVoltage-lcVoltage)*((current-lowCurrent[0])/(highCurrent[0]-lowCurrent[0])) - return 2.5 # Min voltage that cole came up with - - def lookupLerpTorque(self, rpm): # Another piece stolen straight from the one and only Cole Lewis' stuff - RPMtorqueTable = [(0, 178.4), (4120, 178.4), (7500, 82)] # data is [(RPM, Nm), ...] - if (rpm < 0): - return 0 - for i in range(1, len(RPMtorqueTable)): - tup = RPMtorqueTable[i] - if tup[0] == rpm: - return tup[1] - elif tup[0] > rpm: - prevTup = RPMtorqueTable[i-1] - return (prevTup[1] + ((tup[1]-prevTup[1]) * ((rpm-prevTup[0])/(tup[0]-prevTup[0])))) - return 0 - - def calculateMotorRPMfromVelocity(self): - return (self.velocity*60*self.drivenGearTeeth)/(self.drivingGearTeeth * math.pi * self.wheelDiameter) - - def calculateRealTorque(self): - currentTorqueLim = self.torqueDemand - if (self.motorSpeed != 0): - currentLim = min(self.maxDischargeCurrent * self.cellsInParallel, self.maxPower/self.batteryVoltage*1000) - currentTorqueLim = currentLim * self.batteryVoltage * 9550 / (self.motorRPM * 1000) - tractionTorqueLim = self.rearTireTraction * (self.wheelDiameter/2) * self.drivingGearTeeth/self.drivenGearTeeth - - return min(self.torqueDemand, currentTorqueLim, tractionTorqueLim) - - def calculateMotorPower(self): - return self.motorTorque * self.motorRPM / 9550 - - def calculateCurrentDraw(self): - return 1000 * self.motorPower / self.batteryVoltage - - def calculateDischargeDelta(self): - return self.currentDraw * self.timeStepDelta / 3600 - - def calculateTractiveForce(self): - axleTorque = self.motorTorque * self.drivenGearTeeth / self.drivingGearTeeth - wheelForce = axleTorque / (self.wheelDiameter/2) - return wheelForce \ No newline at end of file diff --git a/FullVehicleSim/Electrical/granola2.py b/FullVehicleSim/Electrical/granola2.py deleted file mode 100644 index 79b5b42..0000000 --- a/FullVehicleSim/Electrical/granola2.py +++ /dev/null @@ -1,36 +0,0 @@ -import quesadilla as VoltageTools -import numpy as np -def stepElectrical(worldPrev, worldNext, params, inputs): - - worldNext.wheelRPM = worldPrev.speed / params["mechanical"]["wheelCircumferance"] * 60.0 - - worldNext.wheelRotationsHz = worldPrev.speed / params["mechanical"]["wheelCircumferance"] * 2.0 * np.pi - - worldNext.rpm = worldNext.wheelRPM * params["mechanical"]["gearRatio"] - - worldNext.motorRotationHz = worldNext.wheelRotationsHz * params["mechanical"]["gearRatio"] - - worldNext.maxPower = params["electrical"]["tractiveIMax"] * worldPrev.voltage - - if worldNext.rpm > 7500: - worldNext.torque = worldPrev.drag * params["mechanical"]["wheelRadius"] - else: - if worldNext.maxPower / params["mechanical"]["maxTorque"] < worldNext.motorRotationHz: - perfectTractionTorque = params["mechanical"]["maxTorque"] * params["mechanical"]["gearRatio"] - else: - perfectTractionTorque = params["mechanical"]["maxTorque"] * params["mechanical"]["gearRatio"] - worldNext.torque = min(perfectTractionTorque, worldPrev.maxTractionTorqueAtWheel) - - worldNext.motorTorque = worldNext.torque / params["mechanical"]["gearRatio"] - worldNext.voltage = 28.0 * VoltageTools.lookup(worldPrev.charge, worldPrev.current) - - worldNext.power = worldNext.motorTorque * worldNext.motorRotationHz - - if worldNext.power / worldNext.voltage > params["electrical"]["tractiveIMax"]: - worldNext.current = params["electrical"]["tractiveIMax"] - else: - worldNext.current = worldNext.power / worldNext.voltage - - worldNext.maxTractionTorqueAtWheel = (worldPrev.lbTireTraction.getLongForcePureSlip() + worldPrev.rbTireTraction.getLongForcePureSlip()) * params["mechanical"]["wheelRadius"] - - worldNext.motorForce = worldNext.torque / params["mechanical"]["wheelRadius"] diff --git a/FullVehicleSim/Electrical/lionCellModel.py b/FullVehicleSim/Electrical/lionCellModel.py index e5fa38a..11bab38 100644 --- a/FullVehicleSim/Electrical/lionCellModel.py +++ b/FullVehicleSim/Electrical/lionCellModel.py @@ -325,4 +325,4 @@ def voltage_match(current): upper_voltage = voltage_match(upper) lower_percent = (current_per_cell - lower) / 5.0 voltage = lower_percent * upper_voltage + (1.0 - lower_percent) * lower_voltage - return voltage + return voltage \ No newline at end of file diff --git a/FullVehicleSim/Electrical/powertrain.py b/FullVehicleSim/Electrical/powertrain.py index f583df7..97bbd11 100644 --- a/FullVehicleSim/Electrical/powertrain.py +++ b/FullVehicleSim/Electrical/powertrain.py @@ -34,6 +34,59 @@ 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: + 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 diff --git a/FullVehicleSim/Electrical/salvoltagemodel.md b/FullVehicleSim/Electrical/salvoltagemodel.md new file mode 100644 index 0000000..0736742 --- /dev/null +++ b/FullVehicleSim/Electrical/salvoltagemodel.md @@ -0,0 +1,8 @@ +Nathaniel wanted me to focus on granola2.py and improve what already exists. Since granola2.py is part of the main simulation, I’ve been thinking of it like a game loop: every timestep it takes the previous world state (worldPrev) and computes the next world state (worldNext). worldPrev holds what the car was like one step ago, and worldNext is where we write updated values for the next timestep inside stepElectrical. That includes things like speed-related values, torque, power, voltage, current, and traction limits. +Inside stepElectrical, the sim uses the previous speed plus mechanical parameters (wheel circumference, gear ratio, wheel radius) to compute wheel RPM and motor RPM / motor rotation rate. Then it determines how much torque the car can apply (based on traction and max torque), and computes power from that using the relationship power = torque × rotation rate. That power requirement is what drives how much electrical current the car needs from the battery. +The specific line Nathaniel highlighted is where voltage gets computed, and that voltage ends up affecting basically everything that comes after it. Before my changes, granola2.py computed voltage directly with: + worldNext.voltage = 28.0 * VoltageTools.lookup(worldPrev.charge, worldPrev.current). + That hard-codes the sim to a specific voltage model: it uses a lookup-table cell model for voltage and then multiplies by 28 to convert cell voltage into pack voltage (28 cells in series). After voltage is computed, the sim uses it to compute the next current draw using current = power / voltage, and then clamps it to the max allowed current (tractiveIMax). So voltage directly affects current, which affects charge usage, which affects power limits and the overall behavior of the sim. +lionCellModel.py contains the battery lookup model itself. It’s basically datasheet-derived voltage tables for different discharge currents per cell, and the lookup(charge, current_draw) function converts pack-level values into per-cell values (it divides by 20 as a parallel-cell assumption), converts charge into an “Ah discharged” representation, and then interpolates between tables to estimate the cell voltage. That’s why the main sim multiplies by 28 afterward to get full pack voltage. +The main improvement I made based on Nathaniel’s request was to stop calculating voltage directly inside the main sim loop and instead route it through a dedicated function that takes the previous current and the vehicle state and returns voltage. This decouples the battery voltage model from granola2.py, keeps the sim loop cleaner, and makes it easier to swap in a revised voltage model later (ECM/RC/hysteresis/temp) without having to rewrite the main simulation logic again. +For now the template can wrap the existing lookup behavior so results don’t change, but it creates the interface needed for a revised model later once Nathaniel believes its needed. diff --git a/FullVehicleSim/Electrical/tractiveBattery.py b/FullVehicleSim/Electrical/tractiveBattery.py index d069f2f..b4d5340 100644 --- a/FullVehicleSim/Electrical/tractiveBattery.py +++ b/FullVehicleSim/Electrical/tractiveBattery.py @@ -1,2 +1,47 @@ -from FullVehicleSim.paramLoader import * +from paramLoader import Parameters, Magic +from state import VehicleState +import numpy as np +def calcVoltage(prevWorld: VehicleState) -> float: + delta = 1 / Parameters["stepsPerSecond"] + capacity_Ah = Parameters["cellCapacity_Ah"] + soc = prevWorld.charge + + sigma = Parameters["cellModelSigma"] + hystGain = Parameters["hysteresisGain"] + + # Sliding window: last 10 seconds of current + I_hist = prevWorld.current_history + + # Hysteresis kernel + t = np.arange(prevWorld.current_history.shape[0]) + 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)) + + 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 \ No newline at end of file diff --git a/FullVehicleSim/engine.py b/FullVehicleSim/engine.py index 45abde5..40a915d 100644 --- a/FullVehicleSim/engine.py +++ b/FullVehicleSim/engine.py @@ -22,7 +22,6 @@ def calculateHeading(worldArray:np.ndarray, step:int) -> np.ndarray: ]) new_heading = rotation_matrix @ initial_heading - new_heading = new_heading / np.linalg.norm(new_heading) return np.append(new_heading, 0) @@ -52,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/params.json5 b/FullVehicleSim/params.json5 index 99123d2..461bef2 100644 --- a/FullVehicleSim/params.json5 +++ b/FullVehicleSim/params.json5 @@ -50,8 +50,16 @@ "brakepadThickness": 0.007874, "brakeMass": 0.408, "maxBrakeForce": 1500, + + "seriesCells": 30, + "cellCapacity_Ah": 2.6, + "histeresisKernelLength": 10.0, // seconds }, "Magic": { + + "cellModelSigma": 3.0, + "hysteresisGain": 0.015, + "shape-factor": 0.696268618106842, "curvature-scaling-factor": 1.7177082300186157, "horizontal-shift-factor": 1.0, diff --git a/README.md b/README.md index 4f7f4b2..cafad44 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,8 @@ Vehicle Dynamics Simulations and Data Analysis for Formula Slug written from the ground up in python. +You can install dependencies either from requirements.txt or environment.yml ## Docs - [Full Vehicle Simulation Documentation](FullVehicleSim/README.md) -- [General Guide to fs-data (Accompanying data repository)]() \ No newline at end of file +- [General Guide to fs-data (Accompanying data repository)]() diff --git a/VDCalculators/JeremeyEvanAckermann b/VDCalculators/JeremeyEvanAckermann new file mode 100644 index 0000000..c3cfa0e --- /dev/null +++ b/VDCalculators/JeremeyEvanAckermann @@ -0,0 +1,147 @@ +import os, sys, json +import numpy as np +try: + _HERE = os.path.dirname(__file__) +except NameError: + _HERE = os.getcwd() +def _find_fullvehiclesim(start_dir): #CC: function to find files??? + d = os.path.abspath(start_dir) + for _ in range(10): + cand = os.path.join(d, "FullVehicleSim") + if os.path.isdir(cand): + return cand + d = os.path.dirname(d) + raise FileNotFoundError("Could not find 'FullVehicleSim' directory above current location.") +FULLSIM_PATH = _find_fullvehiclesim(_HERE) +MECH_PATH = os.path.join(FULLSIM_PATH, "Mech") +if FULLSIM_PATH not in sys.path: + sys.path.append(FULLSIM_PATH) +if MECH_PATH not in sys.path: + sys.path.append(MECH_PATH) +import traction # FullVehicleSim/Mech/traction.py +# ---------- Load Parameters + Magic ---------- +PARAMS_PATH = os.path.join(FULLSIM_PATH, "params.json") +with open(PARAMS_PATH, "r") as f: + params = json.load(f) +Parameters = params["Parameters"] +Magic = params["Magic"] +#L, Tf, Tr, h_cg (m) +#m (kg) +#g (m/s^2) +L = 1.589989 +Tf = 1.234008 +Tr = 1.186 +m = 210.92 +g = 9.81 +h_cg = 0.3048 +front_frac = 0.4632 +rear_frac = 1.0 - front_frac +b = front_frac * L +a = L - b +FzF0 = (m * g * front_frac) / 2.0 +FzR0 = (m * g * rear_frac) / 2.0 +def _to_tire_units(alpha_rad: float) -> float: + return float(np.rad2deg(alpha_rad)) +def tire_Fy(Fz, slip_ratio, alpha_rad, U, surfaceTemperature, tirePressure): + t = traction.tire.Tire( + float(Fz), + float(slip_ratio), + _to_tire_units(alpha_rad), + float(U), + float(surfaceTemperature), + float(tirePressure), + Parameters, + Magic + ) + return float(t.getLateralForce()) +def toe_out_ideal(delta_avg_rad: float) -> float: + if abs(delta_avg_rad) < 1e-12: + return 0.0 + R = L / np.tan(delta_avg_rad) + Rin = abs(R) - Tf/2.0 + Rout = abs(R) + Tf/2.0 + sgn = np.sign(delta_avg_rad) + delta_in = sgn * np.arctan(L / Rin) + delta_out = sgn * np.arctan(L / Rout) + return float(delta_in - delta_out) +def inner_outer_angles(delta_avg_rad: float, ack_percent: float): + toe = (ack_percent / 100.0) * toe_out_ideal(delta_avg_rad) + return float(delta_avg_rad + toe/2.0), float(delta_avg_rad - toe/2.0) +# ---------- 4-wheel steady-state solver (RETURNS beta, r ONLY) ---------- +def solve_yaw_rate_4wheel( + U: float, + delta_avg_rad: float, + ack_percent: float = 50.0, + use_load_transfer: bool = True, + slipRatio: float = 0.15, + surfaceTemperature: float = 80.0, + tirePressure: float = 40.0, + max_iter: int = 50, + tol: float = 1e-8 +): + """ + Solve steady-state sideslip beta and yaw rate r using: + (1) Sum(Fy) = m * U * r + (2) a*Sum(Fy_front) = b*Sum(Fy_rear) + Returns exactly: (beta, r) + """ + U = float(U) + if U < 0.2: + return 0.0, 0.0 + beta = 0.0 + r = 0.0 + def residuals(beta, r): + ay = U * r + if use_load_transfer: + phi_f = a / L + phi_r = b / L + dF_front = phi_f * m * ay * h_cg / Tf + dF_rear = phi_r * m * ay * h_cg / Tr + else: + dF_front = 0.0 + dF_rear = 0.0 + if delta_avg_rad >= 0: + Fz_fl, Fz_fr = FzF0 - dF_front, FzF0 + dF_front + Fz_rl, Fz_rr = FzR0 - dF_rear, FzR0 + dF_rear + else: + Fz_fl, Fz_fr = FzF0 + dF_front, FzF0 - dF_front + Fz_rl, Fz_rr = FzR0 + dF_rear, FzR0 - dF_rear + eps = 1.0 + Fz_fl, Fz_fr = max(eps, Fz_fl), max(eps, Fz_fr) + Fz_rl, Fz_rr = max(eps, Fz_rl), max(eps, Fz_rr) + d_in, d_out = inner_outer_angles(delta_avg_rad, ack_percent) + if delta_avg_rad >= 0: + d_fl, d_fr = d_in, d_out + else: + d_fl, d_fr = d_out, d_in + alpha_fl = d_fl - (beta + (a * r / U) + ((Tf/2.0) * r / U)) + alpha_fr = d_fr - (beta + (a * r / U) - ((Tf/2.0) * r / U)) + alpha_rl = 0.0 - (beta + (-b * r / U) + ((Tr/2.0) * r / U)) + alpha_rr = 0.0 - (beta + (-b * r / U) - ((Tr/2.0) * r / U)) + Fy_fl = tire_Fy(Fz_fl, slipRatio, alpha_fl, U, surfaceTemperature, tirePressure) + Fy_fr = tire_Fy(Fz_fr, slipRatio, alpha_fr, U, surfaceTemperature, tirePressure) + Fy_rl = tire_Fy(Fz_rl, slipRatio, alpha_rl, U, surfaceTemperature, tirePressure) + Fy_rr = tire_Fy(Fz_rr, slipRatio, alpha_rr, U, surfaceTemperature, tirePressure) + Fy_total = Fy_fl + Fy_fr + Fy_rl + Fy_rr + Fy_front = Fy_fl + Fy_fr + Fy_rear = Fy_rl + Fy_rr + f1 = Fy_total - m * U * r + f2 = a * Fy_front - b * Fy_rear + return np.array([f1, f2], dtype=float) + for _ in range(max_iter): + F = residuals(beta, r) + if np.linalg.norm(F, 2) < tol: + return float(beta), float(r) + h = 1e-6 + Fb = residuals(beta + h, r) + Fr = residuals(beta, r + h) + J = np.column_stack([(Fb - F)/h, (Fr - F)/h]) + try: + step = np.linalg.solve(J, -F) + except np.linalg.LinAlgError: + step = np.linalg.lstsq(J, -F, rcond=None)[0] + beta += float(step[0]) + r += float(step[1]) + if np.linalg.norm(step, 2) < tol: + return float(beta), float(r) + return float(beta), float(r) \ No newline at end of file diff --git a/VDCalculators/ackermannOptimizer.py b/VDCalculators/ackermannOptimizer.py new file mode 100644 index 0000000..8d051f1 --- /dev/null +++ b/VDCalculators/ackermannOptimizer.py @@ -0,0 +1,166 @@ +import os +import sys +import math +import json +import numpy +sys.path.append("../FullVehicleSim/Mech") + +from steering import * +from traction import * + +import json +magic:dict +parameters:dict +with open('../FullVehicleSim/params.json', 'r') as file: + params = json.load(file) + Magic = params["Magic"] + Parameters = params["Parameters"] + del params + + +def calculateAckermannOther(steerAngle): + #steering wheel angle --> steering rack psition --> wheel steer angle (how static ackermann affects wheel angle function) + # equations taken from "rack and pinion" section of: https://www.mathworks.com/help/vdynblks/ref/kinematicsteering.html + import scipy + import numpy + import matplotlib.pyplot as plt + #global variables + #----------------------------- + # THIS SCRIPT USES MOSTLY FS-3 VALUES. FS-3 values denoted by [3], any theoretical or FS-4 values denoted by [4] + #----------------------------- + tw = 1083.3862 #mm (simplified track width from steering axis to steering axis [steering axis is also simplified to be A-arm knuckle to A-arm knuckle]) + rackRatio = 82.55/248 #[4] mm rack displacement/deg pinion rotation + wheelInput = steerAngle * 180/3.141592 #in degrees of steering wheel movement (CW + CCW -) + rackShift = 0.0 # mm of movement of the rack from left to right (left is - right is +) + l_rack = 292.1 #[4] mm (width of steering rack casing) + l_rod = 378.9426 #[3] mm (length of tie rod as left in FS-3 master CAD) + d = 109.7788 #[3] mm (plan view distance between front axis and rack. negative because we have a front steer setup) + l_arm = 71.628 #[3] mm (length of "steer arm", which is the distance from the center of the upright toe rod pickup to the KPA) + + def rackMovement(): #returns the amount of L-R displacement (in mm) of the steering rack, with the right direction as "positive" + rackShift: float = rackRatio*wheelInput + return rackShift + def calculateAckermann(): #calculates the steer angles of both wheels + l1Left = (0.5*(tw-l_rack)) - rackMovement() #l1 is the instantaneous parallel distance from the rack knuckle to steering axis (KPA). + l1Right = (0.5*(tw-l_rack)) + rackMovement() + l_nought = (0.5*(tw-l_rack)) + beta_nought = betaTrigSolver(l_nought) #used to find the initial "beta" geometry to determine the real steer angle at the wheels + + beta_L = betaTrigSolver(l1Left) - beta_nought #additionally, because there is a static "beta" (simply just arm geometry), we must find the difference to find the actual wheel angles + beta_R = betaTrigSolver(l1Right) - beta_nought + + return beta_L, beta_R + #return beta_nought, betaTrigSolver(l1Left), betaTrigSolver(l1Right) + + def betaTrigSolver(l1): #a separate function to solve the big bad trig equation + l2 = numpy.sqrt((l1**2) + (d**2)) #l2 is the instantaneous direct distance from rack knuckle to steering axis (KPA) + atan = numpy.arctan(d/l1) #first term of the "beta" equation + + num = (l_arm**2) + (l2**2) - (l_rod**2) #just simplifying the calculation of the second term + denom = 2*l_arm*l2 + frac = num/denom + acos = numpy.arccos(frac) + beta = (numpy.pi/2) - atan - acos + return beta + #return frac + inTire, outTire = calculateAckermann() + return (inTire, outTire) + +def calculateSlipAngle(steerAngle, velocity): + Vx = velocity / math.cos(steerAngle) + Vy = velocity * math.tan(steerAngle) + return math.atan(Vy/Vx) + +def calculateYawRate(steerAngle, frontCorneringStiffnessDeg, rearCorneringStiffnessDeg, speed, stepSteerInput): + CF = frontCorneringStiffnessDeg * 180 / np.pi + CR = rearCorneringStiffnessDeg * 180 / np.pi + a = 0.853506 + b = 1.589 - a + m = 277 + I = 658.088580080000 + Y_beta = CF + CR + Y_delta = -CF + N_beta = a * CF - b * CR + N_delta = -1 * a * CF + NR_v = a**2 * CF + b**2 * CR + YR_v = a * CF - b * CR + c = -(NR_v / speed + (I * Y_beta) / (m * speed)) + k = N_beta + (Y_beta * NR_v - N_beta * YR_v) / (m * speed**2) + C2 = (Y_delta * N_beta - Y_beta * N_delta) / (m * speed) + r_inf = (C2 * stepSteerInput) / k + return r_inf + +def solver(minSteer, maxSteer, minVelocity, maxVelocity): + carMass = [277 / 4, 277 / 4 , 277 / 4 , 277 / 4 ] + res = [] + for steer in np.arange(minSteer, maxSteer, 0.1): + currLine = [] + for velocity in np.arange(minVelocity, maxVelocity, 1): + inTire, outTire = calculateAckermannOther(steer) + inCorneringAngleSlip = calculateSlipAngle(inTire, velocity) + outCorneringAngleSlip = calculateSlipAngle(outTire, velocity) + inCorneringStiff = getCorneringStiffness(carMass, (inCorneringAngleSlip, 0), 0.15, velocity, 80, 40, Parameters, Magic)[0] + outCorneringStiff = getCorneringStiffness(carMass, (outCorneringAngleSlip, 0), 0.15, velocity, 80, 40, Parameters, Magic)[0] + netFrontCorneringStiffness = (inCorneringStiff + outCorneringStiff)/2 + netRearCornerningStiffness = -70 # lol who knows bruh + + yawRate = calculateYawRate(steer, netFrontCorneringStiffness, netRearCornerningStiffness, velocity, steer) * -1 + currLine.append(yawRate) + res.append(currLine) + return res + + +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + +# --- solver inputs --- +minSteer = 0 +maxSteer = 1.75 +minVelocity = 10 +maxVelocity = 30 + +# --- run solver --- +yaw_rates = solver(minSteer, maxSteer, minVelocity, maxVelocity) + +# --- build coordinate lists for trisurf --- +steer_vals = np.arange(minSteer, maxSteer, 0.1) +velocity_vals = np.arange(minVelocity, maxVelocity, 1) + +S = [] +V = [] +Z = [] + +for i, steer in enumerate(steer_vals): + for j, velocity in enumerate(velocity_vals): + S.append(steer) + V.append(velocity) + Z.append(yaw_rates[i][j]) + +S = np.array(S) +V = np.array(V) +Z = np.array(Z) + +# --- plot --- +fig = plt.figure(figsize=(12, 8)) +ax = fig.add_subplot(111, projection='3d') + +surf = ax.plot_trisurf( + V, # x-axis + S, # y-axis + Z, # z-axis + cmap='viridis', + linewidth=0.2, + antialiased=True +) + +ax.set_xlabel("Velocity") +ax.set_ylabel("Steering Angle") +ax.set_zlabel("Yaw Rate") +ax.set_title("Yaw Rate vs Steering Angle and Velocity") + +fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, label="Yaw Rate") + +plt.tight_layout() +plt.show() + diff --git a/VDCalculators/ackermannReal.ipynb b/VDCalculators/ackermannReal.ipynb new file mode 100644 index 0000000..335ce30 --- /dev/null +++ b/VDCalculators/ackermannReal.ipynb @@ -0,0 +1,136 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "123e5913", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "defe32a4998c45f294bf70366327b645", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(FloatSlider(value=0.0, description='Deg Wheel Input', max=90.0, min=-90.0, step=1.0), Ou…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#steering wheel angle --> steering rack psition --> wheel steer angle (how static ackermann affects wheel angle function)\n", + "# equations taken from \"rack and pinion\" section of: https://www.mathworks.com/help/vdynblks/ref/kinematicsteering.html\n", + "import scipy\n", + "import numpy\n", + "import matplotlib.pyplot as plt\n", + "import ipywidgets as widgets\n", + "from ipywidgets import interact, interactive\n", + "\n", + "#global variables\n", + "#-----------------------------\n", + "# THIS SCRIPT USES MOSTLY FS-3 VALUES. FS-3 values denoted by [3], any theoretical or FS-4 values denoted by [4]\n", + "#-----------------------------\n", + "tw = 1083.3862 #mm (simplified track width from steering axis to steering axis [steering axis is also simplified to be A-arm knuckle to A-arm knuckle])\n", + "rackRatio = 82.55/248 #[4] mm rack displacement/deg pinion rotation\n", + "wheelInput = 0.0 #in degrees of steering wheel movement (CW + CCW -)\n", + "rackShift = 0.0 # mm of movement of the rack from left to right (left is - right is +)\n", + "l_rack = 292.1 #[4] mm (width of steering rack casing)\n", + "l_rod = 378.9426 #[3] mm (length of tie rod as left in FS-3 master CAD)\n", + "d = 109.7788 #[3] mm (plan view distance between front axis and rack. negative because we have a front steer setup)\n", + "l_arm = 71.628 #[3] mm (length of \"steer arm\", which is the distance from the center of the upright toe rod pickup to the KPA)\n", + "\n", + "def rackMovement(): #returns the amount of L-R displacement (in mm) of the steering rack, with the right direction as \"positive\"\n", + " rackShift: float = rackRatio*wheelInput\n", + " return rackShift\n", + "def calculateAckermann(): #calculates the steer angles of both wheels\n", + " l1Left = (0.5*(tw-l_rack)) - rackMovement() #l1 is the instantaneous parallel distance from the rack knuckle to steering axis (KPA). \n", + " l1Right = (0.5*(tw-l_rack)) + rackMovement()\n", + " l_nought = (0.5*(tw-l_rack))\n", + " beta_nought = betaTrigSolver(l_nought) #used to find the initial \"beta\" geometry to determine the real steer angle at the wheels\n", + "\n", + " beta_L = betaTrigSolver(l1Left) - beta_nought #additionally, because there is a static \"beta\" (simply just arm geometry), we must find the difference to find the actual wheel angles\n", + " beta_R = betaTrigSolver(l1Right) - beta_nought\n", + " \n", + " return beta_L, beta_R\n", + " #return beta_nought, betaTrigSolver(l1Left), betaTrigSolver(l1Right)\n", + " \n", + "def betaTrigSolver(l1): #a separate function to solve the big bad trig equation\n", + " l2 = numpy.sqrt((l1**2) + (d**2)) #l2 is the instantaneous direct distance from rack knuckle to steering axis (KPA)\n", + " atan = numpy.arctan(d/l1) #first term of the \"beta\" equation\n", + "\n", + " num = (l_arm**2) + (l2**2) - (l_rod**2) #just simplifying the calculation of the second term \n", + " denom = 2*l_arm*l2\n", + " frac = num/denom\n", + " acos = numpy.arccos(frac)\n", + " beta = (numpy.pi/2) - atan - acos\n", + " return beta\n", + " #return frac\n", + "\n", + "def update(val): #update variables based off of interact()\n", + " global wheelInput\n", + " wheelInput = val\n", + " left_angle, right_angle = calculateAckermann()\n", + " #stat, bL, bR = calculateAckermann()\n", + " print(f\"wheelInput = {wheelInput}\")\n", + " print(f\"rack movement = {rackMovement()}\")\n", + " print(\"----------------------------\")\n", + " print(f\"left wheel radians = {left_angle}\")\n", + " print(f\"right wheel radians = {right_angle}\")\n", + " print(\"----------------------------\")\n", + " print(f\"left wheel degrees = {numpy.rad2deg(left_angle)}\")\n", + " print(f\"right wheel degrees = {numpy.rad2deg(right_angle)}\")\n", + " #print(f\"Static value must be within [-1,1] = {stat}\")\n", + " #print(f\"Left value must be within [-1,1] = {bL}\")\n", + " #print(f\"Right value must be within [-1,1] = {bR}\")\n", + "\n", + "interact( #ui\n", + " update,\n", + " val=widgets.FloatSlider(value=0.0,min=-90,max=90,step=1,description=\"Deg Wheel Input\")\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20ecdfc9", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/VDCalculators/ackman.ipynb b/VDCalculators/ackman.ipynb new file mode 100644 index 0000000..3d15c8a --- /dev/null +++ b/VDCalculators/ackman.ipynb @@ -0,0 +1,235 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "3a0504ea", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "01db11029d18466a8d2e18a9ef790a0d", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(FloatSlider(value=50.0, description='Ackermann %', min=-100.0, step=5.0), IntSlider(valu…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import os, sys, json\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import ipywidgets as widgets\n", + "from ipywidgets import interact, interactive\n", + "try:\n", + " _HERE = os.path.dirname(__file__)\n", + "except NameError:\n", + " _HERE = os.getcwd()\n", + "def _find_fullvehiclesim(start_dir): #CC: function to find files??? \n", + " d = os.path.abspath(start_dir)\n", + " for _ in range(10):\n", + " cand = os.path.join(d, \"FullVehicleSim\")\n", + " if os.path.isdir(cand):\n", + " return cand\n", + " d = os.path.dirname(d)\n", + " raise FileNotFoundError(\"Could not find 'FullVehicleSim' directory above current location.\")\n", + "FULLSIM_PATH = _find_fullvehiclesim(_HERE)\n", + "MECH_PATH = os.path.join(FULLSIM_PATH, \"Mech\")\n", + "if FULLSIM_PATH not in sys.path:\n", + " sys.path.append(FULLSIM_PATH)\n", + "if MECH_PATH not in sys.path:\n", + " sys.path.append(MECH_PATH)\n", + "import traction # FullVehicleSim/Mech/traction.py\n", + "# ---------- Load Parameters + Magic ----------\n", + "PARAMS_PATH = os.path.join(FULLSIM_PATH, \"params.json\")\n", + "with open(PARAMS_PATH, \"r\") as f:\n", + " params = json.load(f)\n", + "Parameters = params[\"Parameters\"]\n", + "Magic = params[\"Magic\"]\n", + "#L, Tf, Tr, h_cg (m)\n", + "#m (kg)\n", + "#g (m/s^2)\n", + "L = 1.589989 \n", + "Tf = 1.234008\n", + "Tr = 1.186\n", + "m = 210.92\n", + "g = 9.81\n", + "h_cg = 0.3048\n", + "front_frac = 0.4632\n", + "rear_frac = 1.0 - front_frac\n", + "b = front_frac * L\n", + "a = L - b\n", + "FzF0 = (m * g * front_frac) / 2.0 #static weight on one wheel, front and rear\n", + "FzR0 = (m * g * rear_frac) / 2.0\n", + "def _to_tire_units(alpha_rad: float) -> float:\n", + " return float(np.rad2deg(alpha_rad)) #turn radians into degrees (what is alpha_rad?)\n", + "def tire_Fy(Fz, slip_ratio, alpha_rad, U, surfaceTemperature, tirePressure):\n", + " t = traction.tire.Tire(\n", + " float(Fz),\n", + " float(slip_ratio),\n", + " _to_tire_units(alpha_rad),\n", + " float(U),\n", + " float(surfaceTemperature),\n", + " float(tirePressure),\n", + " Parameters,\n", + " Magic\n", + " )\n", + " return float(t.getLateralForce())\n", + "def toe_out_ideal(delta_avg_rad: float) -> float:\n", + " if abs(delta_avg_rad) < 1e-12:\n", + " return 0.0\n", + " R = L / np.tan(delta_avg_rad)\n", + " Rin = abs(R) - Tf/2.0\n", + " Rout = abs(R) + Tf/2.0\n", + " sgn = np.sign(delta_avg_rad)\n", + " delta_in = sgn * np.arctan(L / Rin)\n", + " delta_out = sgn * np.arctan(L / Rout)\n", + " return float(delta_in - delta_out)\n", + "def inner_outer_angles(delta_avg_rad: float, ack_percent: float):\n", + " toe = (ack_percent / 100.0) * toe_out_ideal(delta_avg_rad)\n", + " return float(delta_avg_rad + toe/2.0), float(delta_avg_rad - toe/2.0)\n", + "# ---------- 4-wheel steady-state solver (RETURNS beta, r ONLY) ----------\n", + "def solve_yaw_rate_4wheel(\n", + " U: float,\n", + " delta_avg_rad: float,\n", + " ack_percent: float = 50.0,\n", + " use_load_transfer: bool = True,\n", + " slipRatio: float = 0.15,\n", + " surfaceTemperature: float = 80.0,\n", + " tirePressure: float = 40.0,\n", + " max_iter: int = 50,\n", + " tol: float = 1e-8\n", + "):\n", + " \"\"\"\n", + " Solve steady-state sideslip beta and yaw rate r using:\n", + " (1) Sum(Fy) = m * U * r\n", + " (2) a*Sum(Fy_front) = b*Sum(Fy_rear)\n", + " Returns exactly: (beta, r)\n", + " \"\"\"\n", + " U = float(U)\n", + " if U < 0.2:\n", + " return 0.0, 0.0\n", + " beta = 0.0\n", + " r = 0.0\n", + " def residuals(beta, r):\n", + " ay = U * r\n", + " if use_load_transfer:\n", + " phi_f = a / L\n", + " phi_r = b / L\n", + " dF_front = phi_f * m * ay * h_cg / Tf\n", + " dF_rear = phi_r * m * ay * h_cg / Tr\n", + " else:\n", + " dF_front = 0.0\n", + " dF_rear = 0.0\n", + " if delta_avg_rad >= 0:\n", + " Fz_fl, Fz_fr = FzF0 - dF_front, FzF0 + dF_front\n", + " Fz_rl, Fz_rr = FzR0 - dF_rear, FzR0 + dF_rear\n", + " else:\n", + " Fz_fl, Fz_fr = FzF0 + dF_front, FzF0 - dF_front\n", + " Fz_rl, Fz_rr = FzR0 + dF_rear, FzR0 - dF_rear\n", + " eps = 1.0 #magic number\n", + " Fz_fl, Fz_fr = max(eps, Fz_fl), max(eps, Fz_fr)\n", + " Fz_rl, Fz_rr = max(eps, Fz_rl), max(eps, Fz_rr)\n", + " d_in, d_out = inner_outer_angles(delta_avg_rad, ack_percent)\n", + " if delta_avg_rad >= 0:\n", + " d_fl, d_fr = d_in, d_out\n", + " else:\n", + " d_fl, d_fr = d_out, d_in\n", + " alpha_fl = d_fl - (beta + (a * r / U) + ((Tf/2.0) * r / U))\n", + " alpha_fr = d_fr - (beta + (a * r / U) - ((Tf/2.0) * r / U))\n", + " alpha_rl = 0.0 - (beta + (-b * r / U) + ((Tr/2.0) * r / U))\n", + " alpha_rr = 0.0 - (beta + (-b * r / U) - ((Tr/2.0) * r / U))\n", + " Fy_fl = tire_Fy(Fz_fl, slipRatio, alpha_fl, U, surfaceTemperature, tirePressure)\n", + " Fy_fr = tire_Fy(Fz_fr, slipRatio, alpha_fr, U, surfaceTemperature, tirePressure)\n", + " Fy_rl = tire_Fy(Fz_rl, slipRatio, alpha_rl, U, surfaceTemperature, tirePressure)\n", + " Fy_rr = tire_Fy(Fz_rr, slipRatio, alpha_rr, U, surfaceTemperature, tirePressure)\n", + " Fy_total = Fy_fl + Fy_fr + Fy_rl + Fy_rr\n", + " Fy_front = Fy_fl + Fy_fr\n", + " Fy_rear = Fy_rl + Fy_rr\n", + " f1 = Fy_total - m * U * r\n", + " f2 = a * Fy_front - b * Fy_rear\n", + " return np.array([f1, f2], dtype=float)\n", + " for _ in range(max_iter):\n", + " F = residuals(beta, r)\n", + " if np.linalg.norm(F, 2) < tol:\n", + " return float(beta), float(r)\n", + " h = 1e-6 #magic number\n", + " Fb = residuals(beta + h, r)\n", + " Fr = residuals(beta, r + h)\n", + " J = np.column_stack([(Fb - F)/h, (Fr - F)/h])\n", + " try:\n", + " step = np.linalg.solve(J, -F)\n", + " except np.linalg.LinAlgError:\n", + " step = np.linalg.lstsq(J, -F, rcond=None)[0]\n", + " beta += float(step[0])\n", + " r += float(step[1])\n", + " if np.linalg.norm(step, 2) < tol:\n", + " return float(beta), float(r)\n", + " return float(beta), float(r)\n", + "\n", + "SPEEDS = [10, 15, 20, 30, 40, 50]\n", + "#SPEEDS = [2, 4, 6, 8, 10]\n", + "def plot_yaw_vs_steer_multi_speed(ack_percent, max_steer_deg, steer_step_deg):\n", + " steer_deg = np.arange(-max_steer_deg, max_steer_deg + steer_step_deg, steer_step_deg)\n", + " steer_rad = np.deg2rad(steer_deg)\n", + " plt.figure(figsize=(9, 6))\n", + " for U in SPEEDS:\n", + " r_vals = []\n", + " for d in steer_rad:\n", + " beta, r = solve_yaw_rate_4wheel(U, d, ack_percent) # <-- FIX HERE\n", + " r_vals.append(r)\n", + " plt.plot(steer_deg, r_vals, label=f\"{U} m/s\")\n", + " plt.axhline(0, linewidth=1)\n", + " plt.axvline(0, linewidth=1)\n", + " plt.grid(True)\n", + " plt.xlabel(\"Steering Angle (deg)\")\n", + " plt.ylabel(\"Yaw Rate (rad/s)\")\n", + " plt.title(\"Yaw rate vs steering angle (multi-speed) — Ackermann slider\")\n", + " plt.legend()\n", + " plt.show()\n", + "interact( #where is interact() defined???\n", + " plot_yaw_vs_steer_multi_speed, ack_percent=widgets.FloatSlider(min=-100, max=100, step=5, value=50, description=\"Ackermann %\"),\n", + " max_steer_deg=widgets.IntSlider(min=5, max=30, step=1, value=20, description=\"Max steer (deg)\"),\n", + " steer_step_deg=widgets.FloatSlider(min=0.1, max=5.0, step=0.1, value=1.0,description=\"Steer step (deg)\",readout_format=\".1f\",continuous_update=False)\n", + ")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/VDCalculators/traction.py b/VDCalculators/traction.py new file mode 100644 index 0000000..690e5d2 --- /dev/null +++ b/VDCalculators/traction.py @@ -0,0 +1,21 @@ +from Mech import tireState as tire + +def getTraction(tireLoad, slipAngle, slipRatio, speed, surfaceTemperature, tirePressure, Parameters, Magic): + frontLeft = tire.Tire(tireLoad[0] , 0.15, slipAngle[0], speed, 80, 40, Parameters, Magic) + frontRight = tire.Tire(tireLoad[1] , 0.15, slipAngle[0], speed, 80, 40, Parameters, Magic) + backLeft = tire.Tire(tireLoad[2] , 0.15, slipAngle[1], speed, 80, 40, Parameters, Magic) + backRight = tire.Tire(tireLoad[3] , 0.15, slipAngle[1], speed, 80, 40, Parameters, Magic) + return [(frontLeft.getLongForce(), frontLeft.getLateralForce() * 0.6), + (frontRight.getLongForce() * 0.6, frontRight.getLateralForce() * 0.6), + (backLeft.getLongForce() * 0.6, backLeft.getLateralForce() * 0.6), + (backRight.getLongForce() * 0.6, backRight.getLateralForce() * 0.6)] + +def getCorneringStiffness(tireLoad, slipAngle, slipRatio, speed, surfaceTemperature, tirePressure, Parameters, Magic): + delta = 0.1 + less = getTraction(tireLoad, tuple(x - delta for x in slipAngle), slipRatio, speed, surfaceTemperature, tirePressure, Parameters, Magic) + more = getTraction(tireLoad, tuple(x + delta for x in slipAngle), slipRatio, speed, surfaceTemperature, tirePressure, Parameters, Magic) + + front = ((more[0][1] + more[1][1]) - (less[0][1] + less[1][1])) / (2 * delta) + rear = ((more[2][1] + more[3][1]) - (less[2][1] + less[3][1])) / (2 * delta) + + return (front, rear) diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..df115fb --- /dev/null +++ b/environment.yml @@ -0,0 +1,23 @@ +name: fsae +channels: + - conda-forge + - defaults +dependencies: + - python + - polars + - matplotlib + - pandas + - pytorch=2.4.0 + - scikit-learn + - pygame + - libmagic + - tensorboardx + - pycocotools + - terminaltables + - cython + - torchaudio + - torchvision + - tensorboard + - termcolor + - pip +prefix: /Users/daniel/anaconda3/envs/fsae