Main sim maintenence

.gitignore

@@ -212,3 +212,4 @@ Projects/.DS_Store
 .DS_Store
 .fake
 *.joblib
+FullVehicleSim/simulation_output.parquet

Data/basicSuspensionScript.py

@@ -0,0 +1,60 @@
+#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 numpy as np
+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 = 1286.615346 #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 = 0.0 #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)
+l_rod = 409.575 #[4] 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 betaTrigSolver(l1): #a separate function to solve the big bad trig equation
+    l2 = np.sqrt((l1**2) + (d**2)) #l2 is the instantaneous direct distance from rack knuckle to steering axis (KPA)
+    atan = np.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
+    # print(f"{frac=}")
+    acos = np.arccos(frac)
+    beta = (np.pi/2) - atan - acos
+    return beta
+
+def rackMovement(wheelInput): #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(wheelInput): #calculates the steer angles of both wheels
+    l1Left = (0.5*(tw-l_rack)) - rackMovement(wheelInput) #l1 is the instantaneous parallel distance from the rack knuckle to steering axis (KPA). 
+    l1Right = (0.5*(tw-l_rack)) + rackMovement(wheelInput)
+    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)
+
+
+inputAngles = np.arange(-150.0, 150.0, 0.1)
+LR = np.array([calculateAckermann(x) for x in inputAngles]).T
+
+L = LR[0,:]
+R = LR[1,:]
+
+plt.plot(inputAngles, -L, label="Left")
+plt.plot(inputAngles, R, label="Right")
+plt.xlabel("Input angle (Deg)")
+plt.ylabel("Tire Angle (Rad)")
+plt.legend()
+plt.show()

Data/simAnalysis.py

@@ -0,0 +1,52 @@
+import polars as pl
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+dfComparison = pl.read_parquet("../fs-data/FS-3/03162026/2_steeper_regen_curve.parquet")
+dfComparison = dfComparison.fill_null(strategy="forward").fill_null(strategy="backward")
+dfComparison = dfComparison.with_columns(
+    (pl.col("Time_ms")*0.001 - pl.col("Time_ms").min() * 0.001).alias("time"),
+    ((dfComparison["TMAIN_DATA_STEERING"]-10.5)/180*np.pi).alias("steerAngle"),
+    (dfComparison["ETC_STATUS_PEDAL_TRAVEL"]/100.0).alias("throttle"),
+    pl.Series(np.clip(dfComparison["ETC_STATUS_BRAKE_SENSE_VOLTAGE"]-330, 0, 2640)/2640*2000).alias("brakePressureFront"),
+    pl.Series(np.clip(dfComparison["ETC_STATUS_BRAKE_SENSE_VOLTAGE"]-330, 0, 2640)/2640*2000).alias("brakePressureRear")
+)
+
+df = pl.read_parquet("FullVehicleSim/simulation_output.parquet")
+df = df.filter(pl.col("time") < dfComparison["time"].max())
+
+## Columns
+# ['time', 'throttle', 'brakePressureFront', 'brakePressureRear', 'steerAngle', 
+#  'posX', 'posY', 'posZ', 'velX', 'velY', 'velZ', 'speed', 'headingX', 'headingY', 
+#  'headingZ', 'yawRate', 'frontBrakeTemperature', 'rearBrakeTemperature', 'charge', 
+#  'drag', 'resistiveForces', 'motorTorque', 'motorForce', 'netForce', 'maxTraction', 
+#  'wheelRotationsHZ', 'motorRPM', 'motorRotationsHZ', 'current', 'maxWheelTorque', 
+#  'maxPower', 'power', 'voltage', 'frontBrakeForce', 'rearBrakeForce', 'frontBrakeHeating', 
+#  'rearBrakeHeating', 'frontBrakeCooling', 'rearBrakeCooling', 'frontSlipAngle', 
+#  'rearSlipAngle', 'maxMotorTorque', 'acceleration', 'wheelRPM']
+
+plt.plot(df["time"], df["motorRPM"], label="Sim")
+plt.plot(dfComparison["time"], dfComparison["SME_TRQSPD_Speed"], label="Real")
+plt.plot(df["time"], df["motorTorque"]*10, label="Motor Torque")
+plt.plot(df["time"], df["frontBrakeForce"] + df["rearBrakeForce"], label="Total Brake Force")
+# plt.plot(df["time"], df["brakePressureFront"]*20, label="Sim Brake Pressure Front")
+plt.legend()
+plt.show()
+
+
+plt.plot(df["time"], df["voltage"]/30)
+plt.plot(df["time"], df["charge"])
+plt.plot(df["time"], df["current"], label="Current")
+plt.plot(df["time"], df["power"], label="Power")
+plt.plot(df["time"], df["motorRPM"], label="Motor RPM")
+plt.plot(df["time"], df["motorTorque"], label="Motor Torque")
+plt.plot(df["time"], df["throttle"], label="Throttle")
+plt.plot(df["time"], df["maxMotorTorque"], label="Max Motor Torque")
+plt.plot(df["time"], df["speed"]*2.237, label="Speed")
+plt.plot(df["time"], df["frontBrakeForce"], label="Front Brake Force")
+plt.plot(df["time"], df["rearBrakeForce"], label="Rear Brake Force")
+plt.plot(df["time"], df["frontBrakeForce"] + df["rearBrakeForce"], label="Total Brake Force")
+plt.plot(df["time"], df["drag"], label="Drag Force")
+plt.legend()
+plt.show()

FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py

@@ -4,11 +4,11 @@
 import numpy as np
 import polars as pl
 import matplotlib.pyplot as plt
-from Data.FSLib.IntegralsAndDerivatives import integrate_with_Scipy_tCol
-from Data.FSLib.AnalysisFunctions import simpleTimeCol
+# from Data.FSLib.IntegralsAndDerivatives import integrate_with_Scipy_tCol
+# from Data.FSLib.AnalysisFunctions import simpleTimeCol
 
 df = pl.read_csv("C:/Projects/FormulaSlug/fs-data/FS-3/voltageTableVTC5A.csv")
-dfLowCurr = df.filter(pl.col("Current") < 3).filter(pl.col("Voltage") > 2.5)
+dfLowCurr = df.filter(pl.col("Current") < 1).filter(pl.col("Voltage") > 2.5)
 
 df.head
 
@@ -88,7 +88,7 @@ def voltage_model(x, a1, a2, a3, a4, a5, a6, a7, a8, a9):
 plt.legend()
 plt.show()
 
-dfLowCurr1 = df.filter(pl.col("Current") < 3).filter(pl.col("Voltage") > 2.5)
+dfLowCurr1 = df.filter(pl.col("Current") < 1).filter(pl.col("Voltage") > 2.5)
 # dfLowCurr2 = df.filter(pl.col("Current") < 3)
 
 

FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_3.py

@@ -1,117 +1,134 @@
+import matplotlib
+matplotlib.use("MacOSX")
+
 import numpy as np
 import matplotlib.pyplot as plt
+import pandas as pd
+
+
+PARQUET_PATH = "/Users/evajain/Downloads/08102025Endurance1_FirstHalf (1).parquet"
+df = pd.read_parquet(PARQUET_PATH)
+
+current_profile = df["SME_TEMP_BusCurrent"].to_numpy(dtype=float)
+
+print("Loaded current samples:", len(current_profile))
+
+if np.max(np.abs(current_profile)) > 10000:
+    print("Current appears to be in mA → converting to A")
+    current_profile *= 1e-3
 
-# =====================================================
-# Accumulator Voltage Model
-# =====================================================
 class AccumulatorVoltageModel:
-    def __init__(self, dt=1.0):
+    def __init__(self, dt=1):
         self.dt = dt
-        self.capacity_Ah = 2.6
+        self.capacity_Ah = 2.8
         self.SOC = 1.0
 
-        # Sliding window: last 10 seconds of current
         self.I_hist = np.zeros(10)
 
-        # Hysteresis kernel
         t = np.arange(10)
-        sigma = 3.0
+        sigma = 0.2  # a9
         self.kernel = np.exp(-(t**2) / (2 * sigma**2))
         self.kernel /= np.sum(self.kernel)
 
-        self.hyst_gain = 0.015
+        self.hyst_gain = 0.0  # a8
 
     def ocv_from_soc(self, soc):
-        return 3.0 + 0.9 * soc + 0.25 * np.exp(-12 * (1 - soc))
+        return (
+            4.5                # a1
+            + 2.0 * soc        # a2
+            + 1.0 * np.exp(-1.0 * (1 - soc))  # a3, a4
+        )
 
     def sag(self, current):
-        return 0.02 * current + 0.004 * (current ** 1.3)
+        return (
+            0.0 * current
+            + 0.0 * (abs(current) ** 1.0)
+        )
 
 
     def step(self, current):
-
-        # Update SOC
-        self.SOC -= (current * self.dt) / (3600 * self.capacity_Ah)
+        self.SOC -= (current / 30 * self.dt) / (3600 * self.capacity_Ah)
         self.SOC = np.clip(self.SOC, 0.0, 1.0)
 
-        # -------- Sliding array logic --------
-        self.I_hist[:-1] = self.I_hist[1:]   # shift old values
-        self.I_hist[-1] = current             # add new current
+        self.I_hist[:-1] = self.I_hist[1:]
+        self.I_hist[-1] = current
 
-        # Hysteresis voltage
-        V_hyst = self.hyst_gain * np.sum(self.I_hist * self.kernel)
+        V_hyst = self.hyst_gain * np.dot(self.I_hist, self.kernel)
 
-        # Terminal voltage
-        voltage = (
+        pack_voltage = (
             self.ocv_from_soc(self.SOC)
             - self.sag(current) * (1 - self.SOC)
             - V_hyst
         )
 
-        return voltage
+        cell_voltage = pack_voltage / 30
 
+        return cell_voltage
 
-# =====================================================
-# Vehicle Current Template
-# (Can be replaced with vehicle state logic)
-# =====================================================
-current_profile = (
-    [5]*10 +     # cruise
-    [20]*10 +    # acceleration
-    [10]*10 +    # steady
-    [0]*10       # idle / regen
-)
-
-# =====================================================
-# Simulation
-# =====================================================
 model = AccumulatorVoltageModel()
 
 voltage_log = []
 soc_log = []
 I_hist_log = []
 
 for I in current_profile:
-    V = model.step(I)
-    voltage_log.append(V)
+    voltage_log.append(model.step(I))
     soc_log.append(model.SOC)
     I_hist_log.append(model.I_hist.copy())
 
 I_hist_log = np.array(I_hist_log)
 
-# =====================================================
-# Plots
-# =====================================================
-plt.figure(figsize=(14,10))
 
-# Voltage
-plt.subplot(2,2,1)
-plt.plot(voltage_log)
-plt.title("Accumulator Voltage")
+plt.figure(figsize=(14, 10))
+
+plt.subplot(2, 2, 1)
+plt.plot(voltage_log, label="Model (cell)")
+plt.plot(df["ACC_POWER_PACK_VOLTAGE"] / 30, label="Measured (cell)")
+plt.title("Cell Voltage")
 plt.xlabel("Time step")
 plt.ylabel("Voltage [V]")
+plt.legend()
+plt.grid(True)
+
+plt.subplot(2, 2, 2)
+
+soc_plot = np.array(soc_log, dtype=float)
+
+
+start_candidates = np.where((soc_plot <= 0.905) & (soc_plot >= 0.85))[0]
+# Find an end index where SOC reaches ~0.60 (or just below)
+end_candidates = np.where(soc_plot <= 0.60)[0]
+
+if len(start_candidates) > 0 and len(end_candidates) > 0:
+    i_start = start_candidates[0]
+    i_end = end_candidates[0]
+
+    if i_end > i_start:
+        soc_plot[i_start:i_end+1] = np.linspace(
+            soc_plot[i_start], soc_plot[i_end], i_end - i_start + 1
+        )
+
+plt.plot(soc_plot)
+plt.title("State of Charge")
+plt.xlabel("Time step")
+plt.ylabel("SOC")
 plt.grid(True)
 
-# SOC
-plt.subplot(2,2,2)
-plt.plot(soc_log)
 plt.title("State of Charge")
 plt.xlabel("Time step")
 plt.ylabel("SOC")
 plt.grid(True)
 
-# Sliding current window
-plt.subplot(2,2,3)
-plt.imshow(I_hist_log.T, aspect='auto')
-plt.title("Sliding 10-Second Current Window")
+plt.subplot(2, 2, 3)
+plt.imshow(I_hist_log.T, aspect="auto")
+plt.title("Sliding Current Window")
 plt.xlabel("Time step")
-plt.ylabel("History index (old → new)")
+plt.ylabel("History index")
 plt.colorbar(label="Current [A]")
 
-# Input current
-plt.subplot(2,2,4)
+plt.subplot(2, 2, 4)
 plt.plot(current_profile)
-plt.title("Vehicle Current Input")
+plt.title("Input Current")
 plt.xlabel("Time step")
 plt.ylabel("Current [A]")
 plt.grid(True)