diff --git a/Data/suspensionDamingView.py b/Data/suspensionDamingView.py new file mode 100644 index 0000000..b1a7f1f --- /dev/null +++ b/Data/suspensionDamingView.py @@ -0,0 +1,122 @@ +import polars as pl +import matplotlib.pyplot as plt +import numpy as np +import scipy.fft as fft +from statsmodels.nonparametric.smoothers_lowess import lowess + +df = pl.read_parquet("../fs-data/FS-3/03162026/2_steeper_regen_curve.parquet").fill_null(strategy="forward").fill_null(strategy="backward") + + +[x for x in df.columns if "VDM" in x] + + +FR_base_value = df["TPERIPH_FR_DATA_SUSTRAVEL"][:100].mean() +FL_base_value = df["TPERIPH_FL_DATA_SUSTRAVEL"][:100].mean() +BR_base_value = df["TPERIPH_BR_DATA_SUSTRAVEL"][:100].mean() +BL_base_value = df["TPERIPH_BL_DATA_SUSTRAVEL"][:100].mean() + +t = df["Time_ms"]/1000.0 +fig = plt.figure() +ax = fig.add_subplot(111) +ax.plot(t, df["TPERIPH_FR_DATA_SUSTRAVEL"] - FR_base_value, label="TPERIPH_FR_DATA_SUSTRAVEL") +ax.plot(t, df["TPERIPH_FL_DATA_SUSTRAVEL"] - FL_base_value, label="TPERIPH_FL_DATA_SUSTRAVEL") +ax.plot(t, df["TPERIPH_BR_DATA_SUSTRAVEL"] - BR_base_value, label="TPERIPH_BR_DATA_SUSTRAVEL") +ax.plot(t, df["TPERIPH_BL_DATA_SUSTRAVEL"] - BL_base_value, label="TPERIPH_BL_DATA_SUSTRAVEL") +ax.plot(t, df["VDM_Z_AXIS_YAW_RATE"]/2, label="VDM_Z_AXIS_YAW_RATE") +ax.plot(t, df["SME_TRQSPD_Speed"]/1000, label="SME_TRQSPD_Speed") +ax.legend() +ax.set_xlabel("Time [s]") +ax.set_ylabel("TPERIPH Data") +ax.set_title("Suspension Damping View") +ax.grid(True) +plt.show() + +fig = plt.figure() +ax = fig.add_subplot(111) +ax.plot(t, -1*df["TMAIN_DATA_STEERING"], label="TMAIN_DATA_STEERING") +ax.plot(t, df["VDM_Z_AXIS_YAW_RATE"], label="VDM_Z_AXIS_YAW_RATE") +ax.legend() +ax.set_xlabel("Time [s]") +ax.set_ylabel("Steering and Yaw Rate") +ax.set_title("Steering and Yaw Rate") +ax.grid(True) +plt.show() + +df = df.filter(pl.col("SME_TRQSPD_Speed") > 5000) + +suspension_FR_fft = fft.fft((df["TPERIPH_FR_DATA_SUSTRAVEL"] - FR_base_value).to_numpy()) +suspension_FL_fft = fft.fft((df["TPERIPH_FL_DATA_SUSTRAVEL"] - FL_base_value).to_numpy()) +suspension_BR_fft = fft.fft((df["TPERIPH_BR_DATA_SUSTRAVEL"] - BR_base_value).to_numpy()) +suspension_BL_fft = fft.fft((df["TPERIPH_BL_DATA_SUSTRAVEL"] - BL_base_value).to_numpy()) +suspension_FR_freq = fft.fftfreq(len(suspension_FR_fft), d=0.01) +suspension_FL_freq = fft.fftfreq(len(suspension_FL_fft), d=0.01) +suspension_BR_freq = fft.fftfreq(len(suspension_BR_fft), d=0.01) +suspension_BL_freq = fft.fftfreq(len(suspension_BL_fft), d=0.01) + +fig = plt.figure() +ax1 = fig.add_subplot(221) +ax1.scatter(suspension_FR_freq, np.abs(suspension_FR_fft), label="TPERIPH_FR_DATA_SUSTRAVEL", s=0.5) +ax1.set_title("FR Suspension FFT") +ax1.set_xlabel("Frequency [Hz]") +ax1.set_ylabel("Magnitude") +ax1.set_yscale("log") +ax1.grid(True) +ax2 = fig.add_subplot(222) +ax2.scatter(suspension_FL_freq, np.abs(suspension_FL_fft), label="TPERIPH_FL_DATA_SUSTRAVEL", s=0.5) +ax2.set_title("FL Suspension FFT") +ax2.set_xlabel("Frequency [Hz]") +ax2.set_ylabel("Magnitude") +ax2.set_yscale("log") +ax2.grid(True) +ax3 = fig.add_subplot(223) +ax3.scatter(suspension_BR_freq, np.abs(suspension_BR_fft), label="TPERIPH_BR_DATA_SUSTRAVEL", s=0.5) +ax3.set_title("BR Suspension FFT") +ax3.set_xlabel("Frequency [Hz]") +ax3.set_ylabel("Magnitude") +ax3.set_yscale("log") +ax3.grid(True) +ax4 = fig.add_subplot(224) +ax4.scatter(suspension_BL_freq, np.abs(suspension_BL_fft), label="TPERIPH_BL_DATA_SUSTRAVEL", s=0.5) +ax4.set_title("BL Suspension FFT") +ax4.set_xlabel("Frequency [Hz]") +ax4.set_ylabel("Magnitude") +ax4.set_yscale("log") +ax4.grid(True) +plt.tight_layout() +plt.show() + +yawRate_fft = fft.fft(df["VDM_Z_AXIS_YAW_RATE"].to_numpy()) +yawRate_freq = fft.fftfreq(len(yawRate_fft), d=0.01) +pitchRate_fft = fft.fft(df["VDM_Y_AXIS_YAW_RATE"].to_numpy()) +pitchRate_freq = fft.fftfreq(len(pitchRate_fft), d=0.01) +rollRate_fft = fft.fft(df["VDM_X_AXIS_YAW_RATE"].to_numpy()) +rollRate_freq = fft.fftfreq(len(rollRate_fft), d=0.01) + + +meanSegments = 101 +fig = plt.figure() +ax = fig.add_subplot(111) +ax.scatter(yawRate_freq, np.convolve(np.ones(meanSegments)*1/meanSegments, np.abs(yawRate_fft), mode='same'), label="Yaw", s=0.5) +ax.scatter(pitchRate_freq, np.convolve(np.ones(meanSegments)*1/meanSegments, np.abs(pitchRate_fft), mode='same'), label="Pitch", s=0.5) +ax.scatter(rollRate_freq, np.convolve(np.ones(meanSegments)*1/meanSegments, np.abs(rollRate_fft), mode='same'), label="Roll", s=0.5) +ax.legend() +ax.set_title("Yaw Rate FFT") +ax.set_xlabel("Frequency [Hz]") +ax.set_ylabel("Magnitude") +ax.set_yscale("log") +ax.grid(True) +plt.show() + +## heave based on z acceleration +heave_fft = fft.fft(df["VDM_Z_AXIS_ACCELERATION"].to_numpy()) +heave_freq = fft.fftfreq(len(heave_fft), d=0.01) +fig = plt.figure() +ax = fig.add_subplot(111) +ax.scatter(heave_freq, np.convolve(np.ones(meanSegments)*1/meanSegments, np.abs(heave_fft), mode='same'), label="Heave", s=0.5) +ax.legend() +ax.set_title("Heave FFT") +ax.set_xlabel("Frequency [Hz]") +ax.set_ylabel("Magnitude") +ax.set_yscale("log") +ax.grid(True) +plt.show() \ No newline at end of file diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py b/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py index 3c4713e..8b5fa51 100644 --- a/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py +++ b/FullVehicleSim/Electrical/HisteresisCellModel/batteryModelTraining2.py @@ -4,11 +4,11 @@ import numpy as np import polars as pl import matplotlib.pyplot as plt -from Data.FSLib.IntegralsAndDerivatives import integrate_with_Scipy_tCol -from Data.FSLib.AnalysisFunctions import simpleTimeCol +# from Data.FSLib.IntegralsAndDerivatives import integrate_with_Scipy_tCol +# from Data.FSLib.AnalysisFunctions import simpleTimeCol df = pl.read_csv("C:/Projects/FormulaSlug/fs-data/FS-3/voltageTableVTC5A.csv") -dfLowCurr = df.filter(pl.col("Current") < 3).filter(pl.col("Voltage") > 2.5) +dfLowCurr = df.filter(pl.col("Current") < 1).filter(pl.col("Voltage") > 2.5) df.head @@ -88,7 +88,7 @@ def voltage_model(x, a1, a2, a3, a4, a5, a6, a7, a8, a9): plt.legend() plt.show() -dfLowCurr1 = df.filter(pl.col("Current") < 3).filter(pl.col("Voltage") > 2.5) +dfLowCurr1 = df.filter(pl.col("Current") < 1).filter(pl.col("Voltage") > 2.5) # dfLowCurr2 = df.filter(pl.col("Current") < 3) diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model10.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model10.py new file mode 100644 index 0000000..f97fca2 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model10.py @@ -0,0 +1,133 @@ +import polars as pl +import numpy as np +import matplotlib.pyplot as plt +import scipy.io +from scipy.integrate import cumulative_simpson +from scipy.signal import convolve +from scipy.interpolate import CubicSpline as cs +from scipy.interpolate import LinearNDInterpolator as Lnd_interp +from scipy.optimize import curve_fit +from numba import njit + +# 1. Load Data +matData = scipy.io.loadmat("../fs-data/FS-4/P30B Cell Data/Molicel_INR18650P30B_measurement.mat", squeeze_me=True, struct_as_record=False) +measurements = matData["measurement"] + +DELTA = 0.1 # Stick to 0.1s to capture high-rate dynamic transients accurately +KERNEL_DURATION = 30 +KERNEL_SIZE = int(KERNEL_DURATION / DELTA) +Q_NOMINAL = 3.0 + +def func_SOC(I, t, starting_SOC): + integrated_ah = cumulative_simpson(y=I, x=t, initial=0) / 3600.0 + soc = starting_SOC - (integrated_ah / Q_NOMINAL) # Discharge decreases SOC + return np.clip(soc, 0.0, 1.0) + +lowCurrMeasurements = [m for m in measurements.fu.DCC if "-0.1" in m.name] +lowVs = np.concatenate([sim.V for sim in lowCurrMeasurements]) +lowTs = np.concatenate([sim.T_surf[0, :] for sim in lowCurrMeasurements]) +lowSOC = np.concatenate([func_SOC(sim.I, sim.t, 1.0) for sim in lowCurrMeasurements]) + +voltage_curve = Lnd_interp(np.array([lowVs, lowTs]).T, lowSOC, rescale=True) + +def SOC_lookup(V, T): + SOC = voltage_curve(V, T) + if np.isnan(SOC): + if V > 4.1: return 1.0 + if T < -15.0: return voltage_curve(V, -15.0) + if T > 39.0: return voltage_curve(V, 39.0) + return np.clip(SOC, 0.0, 1.0) + +def buildDF(measurement_set): + frames = [] + for sim in measurement_set: + cubic = cs(sim.t, np.array([sim.I, sim.V, sim.T_surf[0, :], sim.T_surf[1, :]]), axis=1) + t = np.arange(sim.t[0], sim.t[-1], DELTA) + I, V, T_surf1, T_surf2 = cubic(t) + starting_SOC = SOC_lookup(V[0], T_surf1[0]) + SOC = func_SOC(I, t, starting_SOC) + frames.append(pl.DataFrame({ + "I": I, "V": V, "t": t, + "T_surf1": T_surf1, "T_surf2": T_surf2, "SOC": SOC, + })) + return pl.concat(frames) + +all_sims = np.concatenate([measurements.fu.DCC, measurements.fu.CHC, measurements.fu.DCP, measurements.fu.CHP]) +fullData = buildDF(all_sims) + +# OCV Arguments +ocv_args = np.array([2.33586300e+00, 8.70615991e+00, -3.26789721e+01, 7.63836093e+01, -9.98902031e+01, 6.83036001e+01, -1.90310137e+01, 5.90707089e-04]) + +def OCV_terminal_voltage(SOC, T, a0, a1, a2, a3, a4, a5, a6, K_T): + return (a0 + a1 * SOC + a2 * SOC**2 + a3 * SOC**3 + a4 * SOC**4 + a5 * SOC**5 + a6 * SOC**6) * (1 + K_T * (T - 25)) + +@njit +def func_V_X_fast_dynamic(tauX, RX_arr, I, delta): + VX_arr = np.zeros_like(I) + if tauX <= 1e-5: return VX_arr + alpha = np.exp(-delta / tauX) + for i in range(1, len(I)): + beta = RX_arr[i] * (1.0 - alpha) + VX_arr[i] = VX_arr[i - 1] * alpha + I[i] * beta + return VX_arr + +def func_V_h(SOC, M, I, kernel): + # Fix: Hysteresis tracks current DIRECTION (sign), not raw multi-ampere values + I_sign = np.sign(I) + padded_I = np.pad(I_sign, (KERNEL_SIZE - 1, 0), mode='constant') + convolution = convolve(padded_I, kernel[::-1], mode='valid') + return M * SOC * convolution + +def terminal_voltage_dynamic(SOC, T, I, + r0_d0, r0_d1, r0_d2, + r0_c0, r0_c1, r0_c2, + r1_0, r1_1, + k_r, C1, R2, C2, tau_H, M, ocv_params): + + temp_factor = (1.0 + k_r * (T - 25.0)) + + R0_d_stream = (r0_d0 + r0_d1 * SOC + r0_d2 / (SOC + 0.01)) * temp_factor + R0_c_stream = (r0_c0 + r0_c1 * SOC + r0_c2 / (SOC + 0.01)) * temp_factor + R1_stream = (r1_0 + r1_1 * (1.0 - SOC)) * temp_factor + + tau1 = np.mean(R1_stream) * C1 + tau2 = R2 * C2 + + exponent = -np.arange(0, KERNEL_SIZE) * DELTA / tau_H + kernel = np.exp(exponent) / np.sum(np.exp(exponent)) + + V_h = func_V_h(SOC, M, I, kernel) + V_X1 = func_V_X_fast_dynamic(tau1, R1_stream, I, DELTA) + V_X2 = func_V_X_fast_dynamic(tau2, np.full_like(I, R2), I, DELTA) + + V_OCV = OCV_terminal_voltage(SOC, T, *ocv_params) + V_R0 = np.where(I >= 0, R0_d_stream * I, R0_c_stream * I) + + # FIX: Subtract transient polarization voltages (V_X1, V_X2) from V_OCV + V_est = V_OCV - V_R0 - V_X1 - V_X2 + V_h + + print(f"Current Loop Step MAE: {np.mean(np.abs(V_est - fullData['V'].to_numpy())) * 1000.0:.2f} mV") + return V_est + +def wrapper_dynamic(X, r0_d0, r0_d1, r0_d2, r0_c0, r0_c1, r0_c2, r1_0, r1_1, k_r, C1, R2, C2, tau_H, M): + return terminal_voltage_dynamic(X[0], X[1], X[2], r0_d0, r0_d1, r0_d2, r0_c0, r0_c1, r0_c2, r1_0, r1_1, k_r, C1, R2, C2, tau_H, M, ocv_args) + +# Strictly enforce physical boundaries to prevent optimization divergence +# Layout: [r0_d0, r0_d1, r0_d2, r0_c0, r0_c1, r0_c2, r1_0, r1_1, k_r, C1, R2, C2, tau_H, M] +low_b = [1e-4, -0.05, 1e-5, 1e-4, -0.05, 1e-5, 1e-4, 1e-4, -0.05, 10.0, 1e-4, 100.0, 0.5, 1e-4] +upr_b = [0.1, 0.05, 0.01, 0.1, 0.05, 0.01, 0.05, 0.05, 0.0, 5000, 0.1, 50000, 70.0, 0.1] + +p0_dynamic = [0.012, 0.0, 0.0005, 0.010, 0.0, 0.0005, 0.004, 0.002, -0.01, 200, 0.008, 4000, 25, 0.015] + +print("Starting bounded parameter estimation...") +args_dynamic, _ = curve_fit( + wrapper_dynamic, + (fullData["SOC"].to_numpy(), fullData["T_surf1"].to_numpy(), fullData["I"].to_numpy()), + fullData["V"].to_numpy(), + p0=p0_dynamic, + bounds=(low_b, upr_b), # Crucial bounds restored + maxfev=5000 +) + +print("\nFinal Optimized Parameters Array:") +print(args_dynamic) \ No newline at end of file diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model11.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model11.py new file mode 100644 index 0000000..3be52ef --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model11.py @@ -0,0 +1,156 @@ +import polars as pl +import numpy as np +import matplotlib.pyplot as plt +from scipy.optimize import curve_fit, differential_evolution +from scipy.integrate import cumulative_simpson +from scipy.signal import convolve +import scipy.io +from numba import njit +from scipy.interpolate import LinearNDInterpolator as Lnd_interp +from scipy.interpolate import CubicSpline as cs + +# ========================================== +# 1. GLOBAL BASELINE & DATA LOADING +# ========================================== +matData = scipy.io.loadmat("../fs-data/FS-4/P30B Cell Data/Molicel_INR18650P30B_measurement.mat", squeeze_me=True, struct_as_record=False) +measurements = matData["measurement"] + +delta = 0.1 +kernel_duration = 30 +kernel_size = int(kernel_duration / delta) + +lowCurrMeasurements = [m for m in measurements.fu.DCC if "-0.1" in m.name] + +def func_SOC(I, t, starting_SOC): + Q_nominal = 3.0 + integrated_ah = cumulative_simpson(y=I, x=t, initial=0) / 3600.0 + soc = starting_SOC + (integrated_ah / Q_nominal) + return np.clip(soc, 0.0, 1.0) + +lowVs = np.concatenate([sim.V for sim in lowCurrMeasurements]) +lowTs = np.concatenate([sim.T_surf[0, :] for sim in lowCurrMeasurements]) +lowSOC = np.concatenate([func_SOC(sim.I, sim.t, 1.0) for sim in lowCurrMeasurements]) + +voltage_curve = Lnd_interp(np.array([lowVs, lowTs]).T, lowSOC, rescale=True) + +def SOC_lookup(V, T): + SOC = voltage_curve(V, T) + if np.isnan(SOC): + if V > 4.1: return 1.0 + if T < -15.0: return voltage_curve(V, -15.0) + if T > 39.0: return voltage_curve(V, 39.0) + return np.clip(SOC, 0.0, 1.0) + +def buildDF(measurement_set): + frames = [] + for i, sim in enumerate(measurement_set): + cubic = cs(sim.t, np.array([sim.I, sim.V, sim.T_surf[0, :], sim.T_surf[1, :]]), axis=1) + t = np.arange(sim.t[0], sim.t[-1], delta) + I, V, T_surf1, T_surf2 = cubic(t) + starting_SOC = SOC_lookup(V[0], T_surf1[0]) + SOC = func_SOC(I, t, starting_SOC) + frames.append(pl.DataFrame({ + "I": I, "V": V, "t": t, + "T_surf1": T_surf1, "T_surf2": T_surf2, "SOC": SOC, + "run": np.ones(len(I)) * i + })) + return pl.concat(frames) + +all_sims = np.concatenate([measurements.fu.DCC, measurements.fu.CHC, measurements.fu.DCP, measurements.fu.CHP]) +aboveFreezingSims = [m for m in all_sims if m.T_surf[0, 0] > 5] +notHighPulse = [m for m in aboveFreezingSims if not ("-40.000C" in m.name or "-30.000C" in m.name)] +fullData = buildDF(notHighPulse) + +# ========================================== +# 2. CORE PHYSICS & ECM FUNCTIONS +# ========================================== +ocv_args = np.array([2.33586300e+00, 8.70615991e+00, -3.26789721e+01, 7.63836093e+01, -9.98902031e+01, 6.83036001e+01, -1.90310137e+01, 5.90707089e-04]) + +def OCV_terminal_voltage(SOC, T, a0, a1, a2, a3, a4, a5, a6, K_T): + return (a0 + a1 * SOC + a2 * SOC**2 + a3 * SOC**3 + a4 * SOC**4 + a5 * SOC**5 + a6 * SOC**6) * (1 + K_T * (T - 25)) + +def func_V_h(SOC, M, I, kernel): + padded_I = np.pad(I, (kernel_size - 1, 0), mode='constant') + convolution = convolve(padded_I, kernel, mode='valid') + return M * SOC * convolution + +@njit +def func_V_X_fast(tauX, RX, I, delta): + VX_arr = np.zeros_like(I) + alpha = np.exp(-delta / tauX) + beta = RX * (1.0 - alpha) + for i in range(1, len(I)): + VX_arr[i] = VX_arr[i - 1] * alpha + I[i] * beta + return VX_arr + +def terminal_voltage(SOC, T, I, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, ocv_params): + tau1 = R1 * C1 + tau2 = R2 * C2 + + exponent = -np.arange(0, kernel_size) * delta / tau_H + kernel = np.exp(exponent) / np.sum(np.exp(exponent)) + + V_h = func_V_h(SOC, M, I, kernel) + V_X1 = func_V_X_fast(tau1, R1, I, delta) + V_X2 = func_V_X_fast(tau2, R2, I, delta) + + V_OCV = OCV_terminal_voltage(SOC, T, *ocv_params) + V_R0 = np.where(I <= 0, R0_discharge * I, R0_charge * I) + + return V_OCV - V_R0 + V_X1 + V_X2 + V_h + +# ========================================== +# 3. CLEAN WRAPPERS FOR OPTIMIZATION +# ========================================== +def wrapper_post_OCV(X, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M): + modelVoltages = [] + for run_id in np.unique(X[3]): + mask = X[3] == run_id + modelVoltage = terminal_voltage(X[0, mask], X[1, mask], X[2, mask], R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, ocv_args) + modelVoltages.append(modelVoltage) + return np.concatenate(modelVoltages) + +def evolution_wrapper(params, X): + R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M = params + v = wrapper_post_OCV(X, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M) + return np.mean((v - X[4]) ** 2) + +# ========================================== +# 4. WINDOWS-SAFE EXECUTION GUARANTEE +# ========================================== +if __name__ == '__main__': + # 1. Filter out data chunks for optimization + lessData = fullData.filter(pl.col("t") < 1000) + + # Pack array structure cleanly + data_payload = np.array(( + lessData["SOC"].to_numpy(), + lessData["T_surf1"].to_numpy(), + lessData["I"].to_numpy(), + lessData["run"].to_numpy(), + lessData["V"].to_numpy() + )) + + # 2. Set up parameter boundary configurations + bounds = ( + [1e-7, 1e-7, 1e-7, 1e-3, 1e-7, 10.0, 1e-5, -1], # Lower + [0.2, 0.2, 0.2, 5000.0, 0.2, 50000.0, 1e4, 1e3] # Upper + ) + boundsTupled = tuple(zip(*bounds)) + + print("Launching Multi-Threaded Differential Evolution Loop...") + + # 3. Run Global Optimization Engine Safely + argsGlobal = differential_evolution( + evolution_wrapper, + boundsTupled, + args=(data_payload,), + maxiter=35, # Balance iteration density on dynamic pools + popsize=10, # Initial generation seed size + workers=-1, # Utilizes all available logical processor cores + updating='deferred' + ) + + print("\nOptimization Complete!") + print("Optimal Fit Found Parameters:") + print(argsGlobal.x) \ No newline at end of file diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model7.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model7.py new file mode 100644 index 0000000..e8f0600 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model7.py @@ -0,0 +1,816 @@ +"""Electrical Model 7 + +This file is a medium-complexity lithium-ion cell model built for interactive +work rather than CLI execution. + +Design goals: +- Use all available experiment groups that matter here: CHC, DCP, CHP, and + DCC when present. +- Keep the model functional and easy to inspect from an editor or notebook. +- Capture the main battery effects without jumping straight to a full single + particle model. +- Explain each part of the model directly in code. + +Why this approach instead of a full SPM: +Single particle models are powerful, but they are substantially more complex +to identify, tune, and maintain. They usually need more electrochemical state +structure than this dataset clearly provides. For the current use case, a +structured equivalent-circuit model with SOC-dependent OCV, ohmic loss, two +polarization branches, and a hysteresis state is a better middle ground. + +The model below is built in layers: +1. Stack all experiments into a single training set. +2. Estimate a smooth voltage curve from SOC and temperature. +3. Add current-dependent ohmic loss. +4. Add fast and slow polarization states. +5. Add a hysteresis state that depends on current direction and SOC. +6. Fit the whole structure with nonlinear least squares. +""" + +from __future__ import annotations + +import time +from pathlib import Path +from typing import Dict, List, Sequence, Tuple + +import numpy as np +import scipy.io +from scipy.integrate import cumulative_simpson +from scipy.optimize import least_squares + + +# ============================================================ +# PATH RESOLUTION AND DATA LOADING +# ============================================================ + +print("[electrical_model7] attempting to load .mat data (top-level) ...") +matData = scipy.io.loadmat("../fs-data/FS-4/P30B Cell Data/Molicel_INR18650P30B_measurement.mat", squeeze_me=True, struct_as_record=False) +print("[electrical_model7] loaded .mat data (top-level)") + +def LoadMeasurementData(matPath: Path | None = None): + """Load the MAT file and return the `measurement` structure.""" + print("[electrical_model7] LoadMeasurementData: returning measurement struct") + matData = scipy.io.loadmat("../fs-data/FS-4/P30B Cell Data/Molicel_INR18650P30B_measurement.mat", squeeze_me=True, struct_as_record=False) + print("[electrical_model7] LoadMeasurementData: loaded file") + return matData["measurement"] + + +# ============================================================ +# DATA STACKING +# ============================================================ + + +def NormalizeExperimentCollection(experimentCollection) -> List[object]: + """Return a Python list even if MATLAB stored the group as a scalar object.""" + + if isinstance(experimentCollection, np.ndarray): + items = [item for item in experimentCollection.reshape(-1)] + print(f"[electrical_model7] NormalizeExperimentCollection: normalized array -> {len(items)} items") + return items + return [experimentCollection] + + +def SafeSurfaceTemperature(experiment) -> Tuple[np.ndarray, np.ndarray]: + """Extract the two surface temperature channels in a shape-safe way.""" + + surfaceTemperature = np.asarray(experiment.T_surf, dtype=float) + if surfaceTemperature.ndim == 1: + surfaceOne = surfaceTemperature.reshape(-1) + surfaceTwo = surfaceTemperature.reshape(-1) + return surfaceOne, surfaceTwo + + if surfaceTemperature.shape[0] == 1: + surfaceOne = surfaceTemperature[0, :] + surfaceTwo = surfaceTemperature[0, :] + return surfaceOne, surfaceTwo + + surfaceOne = surfaceTemperature[0, :] + surfaceTwo = surfaceTemperature[1, :] + return surfaceOne, surfaceTwo + + +def GuessInitialSoc(groupName: str) -> float: + """Give each experiment family a conservative starting SOC guess.""" + return 1.0 + + +def BuildExperimentRecord(experiment, groupName: str, experimentIndex: int) -> Dict[str, np.ndarray]: + """Convert one raw experiment into a compact record used by the fitter.""" + + currentSeries = np.asarray(experiment.I, dtype=float).reshape(-1) + voltageSeries = np.asarray(experiment.V, dtype=float).reshape(-1) + timeSeries = np.asarray(experiment.t, dtype=float).reshape(-1) + surfaceOne, surfaceTwo = SafeSurfaceTemperature(experiment) + temperatureSeries = 0.5 * (surfaceOne + surfaceTwo) + deltaTimeSeries = np.concatenate([[0.0], np.diff(timeSeries)]) + throughputAh = np.concatenate([[0.0], cumulative_simpson(y=currentSeries, x=timeSeries) / 3600.0]) + + record = { + "groupName": np.array(groupName), + "experimentIndex": np.array(experimentIndex, dtype=int), + "initialSocGuess": np.array(GuessInitialSoc(groupName), dtype=float), + "current": currentSeries, + "voltage": voltageSeries, + "time": timeSeries, + "temperature": temperatureSeries, + "surfaceOne": surfaceOne, + "surfaceTwo": surfaceTwo, + "deltaTime": deltaTimeSeries, + "throughputAh": throughputAh, + } + print(f"[electrical_model7] BuildExperimentRecord: {groupName} idx={experimentIndex} samples={len(currentSeries)}") + return record + + +def StackExperimentGroups(measurementData, groupNames: Sequence[str] = ("CHC", "DCC", "DCP", "CHP")) -> List[Dict[str, np.ndarray]]: + """Stack all requested experiment groups into a flat Python list.""" + + stackedRecords: List[Dict[str, np.ndarray]] = [] + + total = 0 + for groupName in groupNames: + print(f"[electrical_model7] StackExperimentGroups: checking group {groupName}") + if not hasattr(measurementData.fu, groupName): + print(f"[electrical_model7] StackExperimentGroups: group {groupName} not found, skipping") + continue + + experimentCollection = NormalizeExperimentCollection(getattr(measurementData.fu, groupName)) + print(f"[electrical_model7] StackExperimentGroups: found {len(experimentCollection)} experiments in {groupName}") + for experimentIndex, experiment in enumerate(experimentCollection): + stackedRecords.append(BuildExperimentRecord(experiment, groupName, experimentIndex)) + total += 1 + + print(f"[electrical_model7] StackExperimentGroups: stacked total {total} experiments") + return stackedRecords + + +# ============================================================ +# MODEL EQUATIONS +# ============================================================ + + +def Sigmoid(value): + """Smoothly map an unconstrained raw value to [0, 1].""" + + return 1.0 / (1.0 + np.exp(-value)) + + +def SoftPlus(value): + """Smoothly map an unconstrained raw value to a positive number.""" + + return np.log1p(np.exp(-np.abs(value))) + np.maximum(value, 0.0) + + +def Logit(value): + """Inverse of Sigmoid, with clipping for numerical safety.""" + + clippedValue = np.clip(value, 1e-6, 1.0 - 1e-6) + return np.log(clippedValue / (1.0 - clippedValue)) + + +class FitInterrupted(RuntimeError): + """Raised when a fit should stop early but return the best result so far.""" + + +class ResidualAttemptTracker: + """Track residual evaluations, best-so-far state, and stop conditions.""" + + def __init__( + self, + experimentRecords: Sequence[Dict[str, np.ndarray]], + *, + statusEvery: int = 1, + checkpointPath: Path | str | None = None, + stopSignalPath: Path | str | None = None, + wallClockLimitSeconds: float | None = None, + verbose: bool = True, + ) -> None: + self.experimentRecords = experimentRecords + self.statusEvery = max(1, int(statusEvery)) + self.checkpointPath = Path(checkpointPath) if checkpointPath is not None else None + self.stopSignalPath = Path(stopSignalPath) if stopSignalPath is not None else None + self.wallClockLimitSeconds = wallClockLimitSeconds + self.verbose = verbose + self.attemptCount = 0 + self.bestAttempt = 0 + self.bestCost = np.inf + self.bestParameterVector: np.ndarray | None = None + self.bestResidualVector: np.ndarray | None = None + self.startTime = time.perf_counter() + self.parameterNames: List[str] = [] + + def _should_stop(self) -> str | None: + if self.stopSignalPath is not None and self.stopSignalPath.exists(): + return f"stop signal detected at {self.stopSignalPath}" + + if self.wallClockLimitSeconds is not None: + elapsedSeconds = time.perf_counter() - self.startTime + if elapsedSeconds >= self.wallClockLimitSeconds: + return f"wall clock limit reached after {elapsedSeconds:.1f}s" + + return None + + def _writeCheckpoint(self, reason: str) -> None: + if self.checkpointPath is None or self.bestParameterVector is None: + return + + self.checkpointPath.parent.mkdir(parents=True, exist_ok=True) + np.savez( + self.checkpointPath, + reason=np.array(reason), + attemptCount=np.array(self.attemptCount, dtype=int), + bestAttempt=np.array(self.bestAttempt, dtype=int), + bestCost=np.array(self.bestCost, dtype=float), + elapsedSeconds=np.array(time.perf_counter() - self.startTime, dtype=float), + parameterNames=np.array(self.parameterNames, dtype=object), + bestParameterVector=self.bestParameterVector, + bestResidualVector=self.bestResidualVector, + ) + + def evaluate(self, parameterVector: np.ndarray) -> np.ndarray: + stopReason = self._should_stop() + if stopReason is not None: + raise FitInterrupted(stopReason) + + self.attemptCount += 1 + if self.verbose and (self.attemptCount == 1 or self.attemptCount % self.statusEvery == 0): + if np.isfinite(self.bestCost): + print( + f"[electrical_model7] FitModel7: attempt {self.attemptCount}, best cost so far {self.bestCost:.6f}" + ) + else: + print(f"[electrical_model7] FitModel7: attempt {self.attemptCount}, no best yet") + + residualVector = BuildResidualVector( + parameterVector, + self.experimentRecords, + verbose=False, + attemptNumber=self.attemptCount, + ) + cost = 0.5 * float(np.dot(residualVector, residualVector)) + + if cost < self.bestCost: + self.bestCost = cost + self.bestAttempt = self.attemptCount + self.bestParameterVector = parameterVector.copy() + self.bestResidualVector = residualVector.copy() + if self.verbose: + print( + f"[electrical_model7] FitModel7: new best at attempt {self.bestAttempt} with cost {self.bestCost:.6f}" + ) + self._writeCheckpoint("improved") + + return residualVector + + +def BuildParameterVector(experimentRecords: Sequence[Dict[str, np.ndarray]]) -> Tuple[np.ndarray, List[str]]: + """Create an initial parameter vector and the matching parameter names.""" + print(f"[electrical_model7] BuildParameterVector: creating initial vector for {len(experimentRecords)} experiments") + parameterNames = [ + "capacityRaw", + "ocv0", + "ocv1", + "ocv2", + "ocv3", + "ocv4", + "ocv5", + "ocvTemp", + "hyst0", + "hyst1", + "hyst2", + "hystTemp", + "r0Base", + "r0Soc1", + "r0Soc2", + "r0Temp", + "rFastRaw", + "tauFastRaw", + "rSlowRaw", + "tauSlowRaw", + "hystGainRaw", + "hystTauRaw", + "hystShapeRaw", + "hystCurrentScaleRaw", + "tempVoltageCoeff", + ] + + startSocNames = [f"startSocRaw_{index}" for index in range(len(experimentRecords))] + parameterNames.extend(startSocNames) + + parameterVector = np.array( + [ + Logit(3.0 / 4.0), + 3.05, + 1.25, + -0.95, + 0.65, + -0.25, + 0.12, + 0.001, + 0.03, + -0.02, + 0.01, + 0.0005, + np.log(np.exp(0.018) - 1.0), + 0.004, + -0.003, + 0.0006, + np.log(np.exp(0.006) - 1.0), + np.log(np.exp(1.5) - 1.0), + np.log(np.exp(0.0025) - 1.0), + np.log(np.exp(30.0) - 1.0), + np.log(np.exp(0.020) - 1.0), + np.log(np.exp(80.0) - 1.0), + np.log(np.exp(1.2) - 1.0), + np.log(np.exp(5.0) - 1.0), + 0.001, + ] + + [Logit(record["initialSocGuess"]) for record in experimentRecords], + dtype=float, + ) + + print(f"[electrical_model7] BuildParameterVector: parameter vector length {parameterVector.size}") + return parameterVector, parameterNames + + +def UnpackParameters(parameterVector: np.ndarray, experimentCount: int) -> Dict[str, np.ndarray]: + """Split a raw parameter vector into transformed model parameters.""" + print("[electrical_model7] UnpackParameters: unpacking parameter vector") + cursor = 0 + + capacityRaw = parameterVector[cursor] + cursor += 1 + + ocvCoefficients = parameterVector[cursor:cursor + 6] + cursor += 6 + + ocvTemperatureCoeff = parameterVector[cursor] + cursor += 1 + + hystCoefficients = parameterVector[cursor:cursor + 4] + cursor += 4 + + resistanceCoefficients = parameterVector[cursor:cursor + 4] + cursor += 4 + + rFastRaw = parameterVector[cursor] + cursor += 1 + tauFastRaw = parameterVector[cursor] + cursor += 1 + rSlowRaw = parameterVector[cursor] + cursor += 1 + tauSlowRaw = parameterVector[cursor] + cursor += 1 + hystGainRaw = parameterVector[cursor] + cursor += 1 + hystTauRaw = parameterVector[cursor] + cursor += 1 + hystShapeRaw = parameterVector[cursor] + cursor += 1 + hystCurrentScaleRaw = parameterVector[cursor] + cursor += 1 + tempVoltageCoeff = parameterVector[cursor] + cursor += 1 + + startSocRaw = parameterVector[cursor:cursor + experimentCount] + + return { + "capacityAh": SoftPlus(capacityRaw) + 0.5, + "ocvCoefficients": ocvCoefficients, + "ocvTemperatureCoeff": ocvTemperatureCoeff, + "hystCoefficients": hystCoefficients, + "resistanceCoefficients": resistanceCoefficients, + "rFast": SoftPlus(rFastRaw) + 1e-4, + "tauFast": SoftPlus(tauFastRaw) + 1e-3, + "rSlow": SoftPlus(rSlowRaw) + 1e-4, + "tauSlow": SoftPlus(tauSlowRaw) + 1e-3, + "hystGain": SoftPlus(hystGainRaw) + 1e-4, + "hystTau": SoftPlus(hystTauRaw) + 1e-3, + "hystShape": 1.0 + SoftPlus(hystShapeRaw), + "hystCurrentScale": SoftPlus(hystCurrentScaleRaw) + 0.05, + "tempVoltageCoeff": tempVoltageCoeff, + "startSoc": Sigmoid(startSocRaw), + } + + +def EvaluateOcvSurface(socValue: float, temperatureDelta: float, currentDirection: float, parameters: Dict[str, np.ndarray]) -> float: + """Open-circuit voltage surface. + + The base curve captures the SOC shape, the temperature term nudges the + entire curve with temperature, and the direction term allows a small, + smooth charge/discharge split. + """ + + eps = 1e-9 + clippedSoc = np.clip(socValue, eps, 1.0 - eps) + socBasis = np.array([ + 1.0, + clippedSoc, + clippedSoc ** 2, + clippedSoc ** 3, + np.log(clippedSoc), + np.log(1.0 - clippedSoc), + ]) + + baseCurve = float(np.dot(parameters["ocvCoefficients"], socBasis)) + temperatureShift = parameters["ocvTemperatureCoeff"] * temperatureDelta + + directionBasis = np.array([1.0, clippedSoc, clippedSoc ** 2, temperatureDelta]) + directionShift = float(np.dot(parameters["hystCoefficients"], directionBasis)) + + return baseCurve + temperatureShift + currentDirection * directionShift + + +def EvaluateResistanceSurface(socValue: float, temperatureDelta: float, parameters: Dict[str, np.ndarray]) -> float: + """SOC and temperature dependent ohmic resistance.""" + + resistanceBasis = np.array([1.0, socValue, socValue ** 2, temperatureDelta]) + rawResistance = float(np.dot(parameters["resistanceCoefficients"], resistanceBasis)) + return SoftPlus(rawResistance) + 1e-5 + + +def EvaluateHysteresisDrive(socValue: float, temperatureDelta: float, currentDirection: float, parameters: Dict[str, np.ndarray]) -> float: + """Drive term for the hysteresis state. + + This term is not the hysteresis voltage itself. It is the direction-aware + quantity that the state relaxes toward. + """ + + socGate = 1.0 / (1.0 + np.exp(-(socValue - 0.5) / 0.08)) + socBasis = np.array([1.0, socValue, socValue ** 2, temperatureDelta]) + driveShape = float(np.dot(parameters["hystCoefficients"], socBasis)) + return currentDirection * socGate * driveShape + + +def SimulateExperiment( + experimentRecord: Dict[str, np.ndarray], + parameters: Dict[str, np.ndarray], + experimentIndex: int, + verbose: bool = True, +) -> Dict[str, np.ndarray]: + """Run the battery model forward for one experiment. + + The model is causal. Each step updates SOC, the two polarization states, + and the hysteresis state from the previous sample. + """ + + currentSeries = experimentRecord["current"] + voltageSeries = experimentRecord["voltage"] + timeSeries = experimentRecord["time"] + temperatureSeries = experimentRecord["temperature"] + sampleCount = currentSeries.shape[0] + + if verbose: + print(f"[electrical_model7] SimulateExperiment: start exp={experimentIndex} samples={sampleCount}") + progressInterval = max(1, sampleCount // 5) + + predictedVoltage = np.empty(sampleCount, dtype=float) + socSeries = np.empty(sampleCount, dtype=float) + fastPolarization = np.empty(sampleCount, dtype=float) + slowPolarization = np.empty(sampleCount, dtype=float) + hysteresisState = np.empty(sampleCount, dtype=float) + + socSeries[0] = float(parameters["startSoc"][experimentIndex]) + fastPolarization[0] = 0.0 + slowPolarization[0] = 0.0 + hysteresisState[0] = 0.0 + + for sampleIndex in range(sampleCount): + if verbose and sampleIndex % progressInterval == 0: + print(f"[electrical_model7] SimulateExperiment: exp={experimentIndex} progress {sampleIndex}/{sampleCount}") + if sampleIndex > 0: + deltaTime = max(float(timeSeries[sampleIndex] - timeSeries[sampleIndex - 1]), 1e-6) + previousCurrent = float(currentSeries[sampleIndex - 1]) + + socSeries[sampleIndex] = np.clip( + socSeries[sampleIndex - 1] - (previousCurrent * deltaTime) / (3600.0 * parameters["capacityAh"]), + 0.0, + 1.0, + ) + + fastAlpha = np.exp(-deltaTime / parameters["tauFast"]) + slowAlpha = np.exp(-deltaTime / parameters["tauSlow"]) + hystAlpha = np.exp(-deltaTime / parameters["hystTau"]) + + fastPolarization[sampleIndex] = ( + fastAlpha * fastPolarization[sampleIndex - 1] + + parameters["rFast"] * (1.0 - fastAlpha) * previousCurrent + ) + slowPolarization[sampleIndex] = ( + slowAlpha * slowPolarization[sampleIndex - 1] + + parameters["rSlow"] * (1.0 - slowAlpha) * previousCurrent + ) + + currentDirection = np.tanh(previousCurrent / parameters["hystCurrentScale"]) + temperatureDelta = float(temperatureSeries[sampleIndex - 1] - 298.15) + hysteresisDrive = EvaluateHysteresisDrive( + socValue=socSeries[sampleIndex], + temperatureDelta=temperatureDelta, + currentDirection=currentDirection, + parameters=parameters, + ) + hysteresisState[sampleIndex] = ( + hystAlpha * hysteresisState[sampleIndex - 1] + + (1.0 - hystAlpha) * hysteresisDrive + ) + else: + currentDirection = np.tanh(float(currentSeries[sampleIndex]) / parameters["hystCurrentScale"]) + socSeries[sampleIndex] = socSeries[0] + fastPolarization[sampleIndex] = 0.0 + slowPolarization[sampleIndex] = 0.0 + temperatureDelta = float(temperatureSeries[sampleIndex] - 298.15) + hysteresisState[sampleIndex] = EvaluateHysteresisDrive( + socValue=socSeries[sampleIndex], + temperatureDelta=temperatureDelta, + currentDirection=currentDirection, + parameters=parameters, + ) + + currentDirection = np.tanh(float(currentSeries[sampleIndex]) / parameters["hystCurrentScale"]) + temperatureDelta = float(temperatureSeries[sampleIndex] - 298.15) + ocvValue = EvaluateOcvSurface( + socValue=float(socSeries[sampleIndex]), + temperatureDelta=temperatureDelta, + currentDirection=currentDirection, + parameters=parameters, + ) + resistanceValue = EvaluateResistanceSurface( + socValue=float(socSeries[sampleIndex]), + temperatureDelta=temperatureDelta, + parameters=parameters, + ) + + predictedVoltage[sampleIndex] = ( + ocvValue + - float(currentSeries[sampleIndex]) * resistanceValue + - fastPolarization[sampleIndex] + - slowPolarization[sampleIndex] + + parameters["hystGain"] * hysteresisState[sampleIndex] + + parameters["tempVoltageCoeff"] * temperatureDelta + ) + + return { + "predictedVoltage": predictedVoltage, + "soc": socSeries, + "fastPolarization": fastPolarization, + "slowPolarization": slowPolarization, + "hysteresisState": hysteresisState, + "measuredVoltage": voltageSeries, + } + + +def BuildResidualVector( + parameterVector: np.ndarray, + experimentRecords: Sequence[Dict[str, np.ndarray]], + verbose: bool = True, + attemptNumber: int | None = None, +) -> np.ndarray: + """Convert the full fit problem into one stacked residual vector.""" + if verbose: + if attemptNumber is None: + print(f"[electrical_model7] BuildResidualVector: building residuals for {len(experimentRecords)} experiments") + else: + print( + f"[electrical_model7] BuildResidualVector: attempt {attemptNumber}, building residuals for {len(experimentRecords)} experiments" + ) + parameters = UnpackParameters(parameterVector, experimentCount=len(experimentRecords)) + residualParts: List[np.ndarray] = [] + + for experimentIndex, experimentRecord in enumerate(experimentRecords): + simulationResult = SimulateExperiment(experimentRecord, parameters, experimentIndex, verbose=verbose) + residualParts.append(simulationResult["predictedVoltage"] - simulationResult["measuredVoltage"]) + + return np.concatenate(residualParts) + + +def FitModel7( + experimentRecords: Sequence[Dict[str, np.ndarray]], + verbose: bool = True, + statusEvery: int = 1, + wallClockLimitSeconds: float | None = None, + checkpointPath: Path | str | None = None, + stopSignalPath: Path | str | None = None, + maxNfev: int | None = 20, +): + """Fit the medium-complexity lithium-ion model to all stacked experiments.""" + print(f"[electrical_model7] FitModel7: starting fit for {len(experimentRecords)} experiments") + initialParameterVector, parameterNames = BuildParameterVector(experimentRecords) + print(f"[electrical_model7] FitModel7: initial parameter vector length {initialParameterVector.size}") + + tracker = ResidualAttemptTracker( + experimentRecords, + statusEvery=statusEvery, + checkpointPath=checkpointPath, + stopSignalPath=stopSignalPath, + wallClockLimitSeconds=wallClockLimitSeconds, + verbose=verbose, + ) + tracker.parameterNames = parameterNames + + fitResult = None + stoppedEarly = False + stopReason = None + + try: + fitResult = least_squares( + fun=tracker.evaluate, + x0=initialParameterVector, + method="trf", + loss="soft_l1", + f_scale=0.01, + max_nfev=maxNfev, + ) + print("[electrical_model7] FitModel7: least_squares finished") + except (KeyboardInterrupt, FitInterrupted) as interruption: + stoppedEarly = True + stopReason = str(interruption) + print(f"[electrical_model7] FitModel7: stopped early -> {stopReason}") + + if tracker.bestParameterVector is None: + tracker.bestParameterVector = initialParameterVector.copy() + tracker.bestResidualVector = BuildResidualVector(initialParameterVector, experimentRecords, verbose=False) + tracker.bestCost = 0.5 * float(np.dot(tracker.bestResidualVector, tracker.bestResidualVector)) + tracker.bestAttempt = 0 + + bestParameterVector = tracker.bestParameterVector + bestResidualVector = tracker.bestResidualVector + if checkpointPath is not None: + tracker._writeCheckpoint(stopReason or "completed") + + fittedParameters = UnpackParameters(bestParameterVector, experimentCount=len(experimentRecords)) + fittedResults = [] + for experimentIndex, experimentRecord in enumerate(experimentRecords): + fittedResults.append(SimulateExperiment(experimentRecord, fittedParameters, experimentIndex, verbose=False)) + + if verbose: + residualVector = bestResidualVector + rmse = float(np.sqrt(np.mean(residualVector ** 2))) + mae = float(np.mean(np.abs(residualVector))) + if fitResult is not None: + print("Fit success:", fitResult.success) + print("Message:", fitResult.message) + else: + print("Fit stopped early before least_squares returned a result object") + print("Stopped early:", stoppedEarly) + print("Stop reason:", stopReason) + print("Attempts:", tracker.attemptCount) + print("Best attempt:", tracker.bestAttempt) + print(f"RMSE: {rmse:.6f} V") + print(f"MAE: {mae:.6f} V") + + return { + "fitResult": fitResult, + "parameterNames": parameterNames, + "parameterVector": bestParameterVector, + "fittedParameters": fittedParameters, + "experimentRecords": list(experimentRecords), + "fittedResults": fittedResults, + "attemptCount": tracker.attemptCount, + "bestAttempt": tracker.bestAttempt, + "bestCost": tracker.bestCost, + "stoppedEarly": stoppedEarly, + "stopReason": stopReason, + } + + +# ============================================================ +# INTERACTIVE HELPERS +# ============================================================ + + +def BuildDefaultTrainingSet(groupNames: Sequence[str] = ("CHC", "DCC", "DCP", "CHP")) -> List[Dict[str, np.ndarray]]: + """Load the MAT file and stack all available experiment groups.""" + print(f"[electrical_model7] BuildDefaultTrainingSet: groups={groupNames}") + measurementData = LoadMeasurementData() + records = StackExperimentGroups(measurementData, groupNames=groupNames) + print(f"[electrical_model7] BuildDefaultTrainingSet: total records={len(records)}") + return records + + +def PlotExperimentFit(experimentRecord: Dict[str, np.ndarray], simulationResult: Dict[str, np.ndarray], titlePrefix: str = "Model 7"): + """Plot measured and modeled voltage for one experiment.""" + + print(f"[electrical_model7] PlotExperimentFit: plotting experiment {experimentRecord.get('groupName')} idx={experimentRecord.get('experimentIndex')}") + + import matplotlib.pyplot as plt + + timeSeries = experimentRecord["time"] + measuredVoltage = simulationResult["measuredVoltage"] + predictedVoltage = simulationResult["predictedVoltage"] + currentSeries = experimentRecord["current"] + + figure, axes = plt.subplots(3, 1, figsize=(12, 10), sharex=True) + + axes[0].plot(timeSeries, currentSeries, color="tab:blue") + axes[0].set_ylabel("Current [A]") + axes[0].grid(True) + + axes[1].plot(timeSeries, measuredVoltage, label="Measured", color="tab:orange") + axes[1].plot(timeSeries, predictedVoltage, label="Modeled", color="tab:green", alpha=0.85) + axes[1].set_ylabel("Voltage [V]") + axes[1].legend() + axes[1].grid(True) + + axes[2].plot(timeSeries, simulationResult["soc"], color="tab:purple") + axes[2].set_ylabel("SOC") + axes[2].set_xlabel("Time [s]") + axes[2].grid(True) + + groupName = str(experimentRecord["groupName"]) + experimentIndex = int(experimentRecord["experimentIndex"]) + figure.suptitle(f"{titlePrefix}: {groupName}_{experimentIndex}") + figure.tight_layout() + return figure, axes + + +def SummarizeModelFit(fitPackage: Dict[str, object]) -> Dict[str, float]: + """Return high-level error metrics for the full stacked fit.""" + + print("[electrical_model7] SummarizeModelFit: computing summary metrics") + experimentRecords = fitPackage["experimentRecords"] + fittedParameters = fitPackage["fittedParameters"] + + residualParts = [] + for experimentIndex, experimentRecord in enumerate(experimentRecords): + simulationResult = SimulateExperiment(experimentRecord, fittedParameters, experimentIndex) + residualParts.append(simulationResult["predictedVoltage"] - simulationResult["measuredVoltage"]) + + residualVector = np.concatenate(residualParts) + return { + "rmse": float(np.sqrt(np.mean(residualVector ** 2))), + "mae": float(np.mean(np.abs(residualVector))), + "maxAbsError": float(np.max(np.abs(residualVector))), + } + + +# ============================================================ +# USAGE EXAMPLE +# ============================================================ + +# The file intentionally does not auto-run a fit on import. +if __name__ == "__main__": + print("[electrical_model7] Example run: building default training set...") + experimentRecords = BuildDefaultTrainingSet() + print("[electrical_model7] Example run: fitting model to stacked experiments...") + fitPackage = FitModel7(experimentRecords) + print("[electrical_model7] Example run: summarizing fit metrics...") + metrics = SummarizeModelFit(fitPackage) + print(f"[electrical_model7] Example run: metrics={metrics}") + firstExperimentResult = fitPackage["fittedResults"][0] + print("[electrical_model7] Example run: plotting first experiment fit...") + PlotExperimentFit(experimentRecords[0], firstExperimentResult) + + +# ============================================================ +# COMPACT MODEL EQUATIONS AND PARAMETERS +# ============================================================ +# Use this block as the single reference form for testing the model. +# +# State updates for sample k -> k+1: +# SOC[k+1] = clip(SOC[k] - I[k] * dt / (3600 * C), 0, 1) +# a_f = exp(-dt / tau_f) +# a_s = exp(-dt / tau_s) +# a_h = exp(-dt / tau_h) +# x_f[k+1] = a_f * x_f[k] + R_f * (1 - a_f) * I[k] +# x_s[k+1] = a_s * x_s[k] + R_s * (1 - a_s) * I[k] +# dir[k] = tanh(I[k] / I_scale) +# T_delta = T[k] - 298.15 +# g_soc = 1 / (1 + exp(-(SOC[k] - 0.5) / 0.08)) +# h_drive = dir[k] * g_soc * (h0 + h1*SOC[k] + h2*SOC[k]^2 + hT*T_delta) +# h[k+1] = a_h * h[k] + (1 - a_h) * h_drive +# +# Voltage model: +# OCV(SOC, T, dir) = +# ocv0 + ocv1*SOC + ocv2*SOC^2 + ocv3*SOC^3 + ocv4*ln(SOC) + ocv5*ln(1-SOC) +# + ocvT * T_delta +# + dir * (d0 + d1*SOC + d2*SOC^2 + dT*T_delta) +# +# R0(SOC, T) = softplus(r0 + r1*SOC + r2*SOC^2 + rT*T_delta) + 1e-5 +# +# V_hat[k] = OCV(SOC[k], T_delta[k], dir[k]) +# - I[k] * R0(SOC[k], T_delta[k]) +# - x_f[k] - x_s[k] +# + h_g * h[k] +# + vT * T_delta[k] +# +# Parameters (raw vector order before transforms): +# capacityRaw +# ocv0, ocv1, ocv2, ocv3, ocv4, ocv5, ocvTemp +# hyst0, hyst1, hyst2, hystTemp +# r0Base, r0Soc1, r0Soc2, r0Temp +# rFastRaw, tauFastRaw, rSlowRaw, tauSlowRaw +# hystGainRaw, hystTauRaw, hystShapeRaw, hystCurrentScaleRaw +# tempVoltageCoeff +# startSocRaw_0 ... startSocRaw_(N-1) +# +# Transforms: +# C = softplus(capacityRaw) + 0.5 +# R_f = softplus(rFastRaw) + 1e-4 +# tau_f = softplus(tauFastRaw) + 1e-3 +# R_s = softplus(rSlowRaw) + 1e-4 +# tau_s = softplus(tauSlowRaw) + 1e-3 +# h_g = softplus(hystGainRaw) + 1e-4 +# tau_h = softplus(hystTauRaw) + 1e-3 +# h_shape = 1 + softplus(hystShapeRaw) +# I_scale = softplus(hystCurrentScaleRaw) + 0.05 +# SOC0_i = sigmoid(startSocRaw_i) diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_modelG13.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_modelG13.py new file mode 100644 index 0000000..9a4f416 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_modelG13.py @@ -0,0 +1,194 @@ +import polars as pl +import numpy as np +import matplotlib.pyplot as plt +from scipy.optimize import curve_fit +from scipy.integrate import cumulative_simpson +import scipy.io +from numba import njit +from scipy.interpolate import LinearNDInterpolator as Lnd_interp +from scipy.interpolate import CubicSpline as cs + +# 1. Load Data +matData = scipy.io.loadmat("../fs-data/FS-4/P30B Cell Data/Molicel_INR18650P30B_measurement.mat", squeeze_me=True, struct_as_record=False) +measurements = matData["measurement"] + +delta = 0.1 + +lowCurrMeasurements = [m for m in measurements.fu.DCC if "-0.1" in m.name] + +def func_SOC(I, t, starting_SOC): + Q_nominal = 3.0 # Ah capacity of Molicel P30B + integrated_ah = cumulative_simpson(y=I, x=t, initial=0) / 3600.0 + soc = starting_SOC + (integrated_ah / Q_nominal) + return np.clip(soc, 0.0, 1.0) + +lowVs = np.concatenate([sim.V for sim in lowCurrMeasurements]) +lowTs = np.concatenate([sim.T_surf[0, :] for sim in lowCurrMeasurements]) +lowSOC = np.concatenate([func_SOC(sim.I, sim.t, 1.0) for sim in lowCurrMeasurements]) + +voltage_curve = Lnd_interp(np.array([lowVs, lowTs]).T, lowSOC, rescale=True) + +def SOC_lookup(V, T): + # Proactively clip boundaries to prevent out-of-bounds NaNs in 2D interpolation + V_clipped = np.clip(V, lowVs.min(), lowVs.max()) + T_clipped = np.clip(T, lowTs.min(), lowTs.max()) + SOC = voltage_curve(V_clipped, T_clipped) + return np.clip(SOC, 0.0, 1.0) + +def buildDF(measurement_set): + frames = [] + for i, sim in enumerate(measurement_set): + cubic = cs(sim.t, np.array([sim.I, sim.V, sim.T_surf[0, :], sim.T_surf[1, :]]), axis=1) + t = np.arange(sim.t[0], sim.t[-1], delta) + I, V, T_surf1, T_surf2 = cubic(t) + starting_SOC = SOC_lookup(V[0], T_surf1[0]) + SOC = func_SOC(I, t, starting_SOC) + frames.append(pl.DataFrame({ + "I": I, "V": V, "t": t, + "T_surf1": T_surf1, "T_surf2": T_surf2, "SOC": SOC, + "run": np.ones(len(I)) * i + })) + return pl.concat(frames) + +all_sims = np.concatenate([measurements.fu.DCC, measurements.fu.CHC, measurements.fu.DCP, measurements.fu.CHP]) +aboveFreezingSims = [m for m in all_sims if m.T_surf[0, 0] > -2.0] +notHighPulse = [m for m in aboveFreezingSims if not ("-40.000C" in m.name or "-30.000C" in m.name)] +fullData = buildDF(notHighPulse) + +# 2. Extract and Freeze OCV +def OCV_terminal_voltage(SOC, T, a0, a1, a2, a3, a4, a5, a6, K_T): + return (a0 + a1 * SOC + a2 * SOC**2 + a3 * SOC**3 + a4 * SOC**4 + a5 * SOC**5 + a6 * SOC**6) * (1 + K_T * (T - 25.0)) + +ocv_args = np.array([2.33586300e+00, 8.70615991e+00, -3.26789721e+01, 7.63836093e+01, + -9.98902031e+01, 6.83036001e+01, -1.90310137e+01, 5.90707089e-04]) + +# 3. Optimized Discrete State-Space Tracking +@njit +def func_V_h_fast(tau_H, M, I, delta): + Vh_arr = np.zeros_like(I) + alpha = np.exp(-delta / tau_H) + beta = 1.0 - alpha + current_direction = np.sign(I) + for i in range(1, len(I)): + Vh_arr[i] = Vh_arr[i - 1] * alpha + current_direction[i] * beta + return M * Vh_arr + +@njit +def func_V_X_fast(tauX, RX, I, delta): + VX_arr = np.zeros_like(I) + # Safely handles array inputs generated by dynamic temperature/current conditions + alpha = np.exp(-delta / tauX) + beta = RX * (1.0 - alpha) + for i in range(1, len(I)): + VX_arr[i] = VX_arr[i - 1] * alpha[i] + I[i] * beta[i] + return VX_arr + +def terminal_voltage(SOC, T, I, r0_d, r0_c, R1_ref, C1, R2_ref, C2, tau_H, M, E_a, k_I): + T_kelvin = T + 273.15 + arrhenius_mult = np.exp((E_a * 1000.0) * (1.0 / T_kelvin - 1.0 / 298.15)) + + # Vectorized evaluation using scaled vector definitions + R0_d_base = (r0_d[0] + r0_d[1] * SOC + r0_d[2] * SOC**2) * arrhenius_mult + R0_c_base = (r0_c[0] + r0_c[1] * SOC + r0_c[2] * SOC**2) * arrhenius_mult + + R1 = R1_ref * arrhenius_mult * (1.0 + k_I * np.exp(-np.abs(I))) + R2 = R2_ref * arrhenius_mult + + tau1 = R1 * C1 + tau2 = R2 * C2 + + V_h = func_V_h_fast(tau_H, M, I, delta) + V_X1 = func_V_X_fast(tau1, R1, I, delta) + V_X2 = func_V_X_fast(tau2, R2, I, delta) + + V_OCV = OCV_terminal_voltage(SOC, T, *ocv_args) + V_R0 = np.where(I <= 0, R0_d_base * I, R0_c_base * I) + + return V_OCV + V_R0 + V_X1 + V_X2 + V_h + +# 4. Nonlinear Multi-Variable Optimization Wrapper +def wrapper_post_OCV(X, r0_d0, r0_d1, r0_d2, r0_c0, r0_c1, r0_c2, R1_scaled, C1_scaled, R2_scaled, C2_scaled, tau_H, M_scaled, E_a_scaled, k_I): + r0_d = np.array([r0_d0, r0_d1, r0_d2]) * 1e-3 + r0_c = np.array([r0_c0, r0_c1, r0_c2]) * 1e-3 + R1_ref = R1_scaled * 1e-3 + C1 = C1_scaled * 1e2 + R2_ref = R2_scaled * 1e-3 + C2 = C2_scaled * 1e4 + M = M_scaled * 1e-2 + E_a = E_a_scaled * 1.0 + + modelVoltages = np.zeros(X.shape[1]) + run_ids = X[3, :] + + for id in np.unique(run_ids): + mask = run_ids == id + + # Pass properly scaled R0 coefficient vectors directly + modelVoltages[mask] = terminal_voltage( + X[0, mask], X[1, mask], X[2, mask], + r0_d, r0_c, + R1_ref, C1, R2_ref, C2, tau_H, M, E_a, k_I + ) + + actual_voltages = X[4, :] + current_mae = np.mean(np.abs(modelVoltages - actual_voltages)) + print(f"Iteration MAE: {current_mae:.5f} V") + + return modelVoltages + +# 5. Boundaries and Initial Guess Configuration +scaled_bounds = ( + [2.0, -15.0, 0.0, 2.0, -15.0, 0.0, 0.5, 1.0, 1.0, 0.5, 0.5, 0.01, 1.0, 0.01], + [30.0, 5.0, 40.0, 30.0, 5.0, 40.0, 20.0, 10.0, 30.0, 5.0, 120.0, 5.00, 5.0, 2.00] +) + +p0_scaled = [ + 14.0, -5.0, 10.0, + 16.0, -5.0, 10.0, + 6.0, 4.0, + 12.0, 1.5, + 15.0, 0.8, + 2.5, 0.2 +] + +lessData = fullData.gather_every(10) + +X_data = np.array(( + lessData["SOC"].to_numpy(), + lessData["T_surf1"].to_numpy(), + lessData["I"].to_numpy(), + lessData["run"].to_numpy(), + lessData["V"].to_numpy() +)) + +print(f"Starting advanced non-linear optimization with {X_data.shape} points...") +optimized_scaled_args, _ = curve_fit( + wrapper_post_OCV, + X_data, + lessData["V"].to_numpy(), + p0=p0_scaled, + bounds=scaled_bounds, + maxfev=3500, + method='trf', + x_scale='jac' +) + +final_params = { + "r0_discharge_poly": [optimized_scaled_args[0]*1e-3, optimized_scaled_args[1]*1e-3, optimized_scaled_args[2]*1e-3], + "r0_charge_poly": [optimized_scaled_args[3]*1e-3, optimized_scaled_args[4]*1e-3, optimized_scaled_args[5]*1e-3], + "R1_ref": optimized_scaled_args[6] * 1e-3, + "C1": optimized_scaled_args[7] * 1e2, + "R2_ref": optimized_scaled_args[8] * 1e-3, + "C2": optimized_scaled_args[9] * 1e4, + "tau_H": optimized_scaled_args[10], + "M": optimized_scaled_args[11] * 1e-2, + "E_a_over_Rg": optimized_scaled_args[12] * 1000.0, + "k_I_current_scale": optimized_scaled_args[13] +} + +print("\nOptimized Multidimensional Parameters:") +for k, v in final_params.items(): + if isinstance(v, (list, tuple, np.ndarray)): + print(f"{k}: {[f'{float(x):.6e}' for x in v]}") + else: + print(f"{k}: {float(v):.6e}") \ No newline at end of file diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_3.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_3.py index 062a608..76570c3 100644 --- a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_3.py +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_3.py @@ -1,70 +1,70 @@ +import matplotlib +matplotlib.use("MacOSX") + import numpy as np import matplotlib.pyplot as plt +import pandas as pd + + +PARQUET_PATH = "/Users/evajain/Downloads/08102025Endurance1_FirstHalf (1).parquet" +df = pd.read_parquet(PARQUET_PATH) + +current_profile = df["SME_TEMP_BusCurrent"].to_numpy(dtype=float) + +print("Loaded current samples:", len(current_profile)) + +if np.max(np.abs(current_profile)) > 10000: + print("Current appears to be in mA → converting to A") + current_profile *= 1e-3 -# ===================================================== -# Accumulator Voltage Model -# ===================================================== class AccumulatorVoltageModel: - def __init__(self, dt=1.0): + def __init__(self, dt=1): self.dt = dt - self.capacity_Ah = 2.6 + self.capacity_Ah = 2.8 self.SOC = 1.0 - # Sliding window: last 10 seconds of current self.I_hist = np.zeros(10) - # Hysteresis kernel t = np.arange(10) - sigma = 3.0 + sigma = 0.2 # a9 self.kernel = np.exp(-(t**2) / (2 * sigma**2)) self.kernel /= np.sum(self.kernel) - self.hyst_gain = 0.015 + self.hyst_gain = 0.0 # a8 def ocv_from_soc(self, soc): - return 3.0 + 0.9 * soc + 0.25 * np.exp(-12 * (1 - soc)) + return ( + 4.5 # a1 + + 2.0 * soc # a2 + + 1.0 * np.exp(-1.0 * (1 - soc)) # a3, a4 + ) def sag(self, current): - return 0.02 * current + 0.004 * (current ** 1.3) + return ( + 0.0 * current + + 0.0 * (abs(current) ** 1.0) + ) def step(self, current): - - # Update SOC - self.SOC -= (current * self.dt) / (3600 * self.capacity_Ah) + self.SOC -= (current / 30 * self.dt) / (3600 * self.capacity_Ah) self.SOC = np.clip(self.SOC, 0.0, 1.0) - # -------- Sliding array logic -------- - self.I_hist[:-1] = self.I_hist[1:] # shift old values - self.I_hist[-1] = current # add new current + self.I_hist[:-1] = self.I_hist[1:] + self.I_hist[-1] = current - # Hysteresis voltage - V_hyst = self.hyst_gain * np.sum(self.I_hist * self.kernel) + V_hyst = self.hyst_gain * np.dot(self.I_hist, self.kernel) - # Terminal voltage - voltage = ( + pack_voltage = ( self.ocv_from_soc(self.SOC) - self.sag(current) * (1 - self.SOC) - V_hyst ) - return voltage + cell_voltage = pack_voltage / 30 + return cell_voltage -# ===================================================== -# Vehicle Current Template -# (Can be replaced with vehicle state logic) -# ===================================================== -current_profile = ( - [5]*10 + # cruise - [20]*10 + # acceleration - [10]*10 + # steady - [0]*10 # idle / regen -) - -# ===================================================== -# Simulation -# ===================================================== model = AccumulatorVoltageModel() voltage_log = [] @@ -72,46 +72,63 @@ def step(self, current): I_hist_log = [] for I in current_profile: - V = model.step(I) - voltage_log.append(V) + voltage_log.append(model.step(I)) soc_log.append(model.SOC) I_hist_log.append(model.I_hist.copy()) I_hist_log = np.array(I_hist_log) -# ===================================================== -# Plots -# ===================================================== -plt.figure(figsize=(14,10)) -# Voltage -plt.subplot(2,2,1) -plt.plot(voltage_log) -plt.title("Accumulator Voltage") +plt.figure(figsize=(14, 10)) + +plt.subplot(2, 2, 1) +plt.plot(voltage_log, label="Model (cell)") +plt.plot(df["ACC_POWER_PACK_VOLTAGE"] / 30, label="Measured (cell)") +plt.title("Cell Voltage") plt.xlabel("Time step") plt.ylabel("Voltage [V]") +plt.legend() +plt.grid(True) + +plt.subplot(2, 2, 2) + +soc_plot = np.array(soc_log, dtype=float) + + +start_candidates = np.where((soc_plot <= 0.905) & (soc_plot >= 0.85))[0] +# Find an end index where SOC reaches ~0.60 (or just below) +end_candidates = np.where(soc_plot <= 0.60)[0] + +if len(start_candidates) > 0 and len(end_candidates) > 0: + i_start = start_candidates[0] + i_end = end_candidates[0] + + if i_end > i_start: + soc_plot[i_start:i_end+1] = np.linspace( + soc_plot[i_start], soc_plot[i_end], i_end - i_start + 1 + ) + +plt.plot(soc_plot) +plt.title("State of Charge") +plt.xlabel("Time step") +plt.ylabel("SOC") plt.grid(True) -# SOC -plt.subplot(2,2,2) -plt.plot(soc_log) plt.title("State of Charge") plt.xlabel("Time step") plt.ylabel("SOC") plt.grid(True) -# Sliding current window -plt.subplot(2,2,3) -plt.imshow(I_hist_log.T, aspect='auto') -plt.title("Sliding 10-Second Current Window") +plt.subplot(2, 2, 3) +plt.imshow(I_hist_log.T, aspect="auto") +plt.title("Sliding Current Window") plt.xlabel("Time step") -plt.ylabel("History index (old → new)") +plt.ylabel("History index") plt.colorbar(label="Current [A]") -# Input current -plt.subplot(2,2,4) +plt.subplot(2, 2, 4) plt.plot(current_profile) -plt.title("Vehicle Current Input") +plt.title("Input Current") plt.xlabel("Time step") plt.ylabel("Current [A]") plt.grid(True) diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_5.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_5.py new file mode 100644 index 0000000..b502bcf --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_5.py @@ -0,0 +1,236 @@ +import matplotlib +matplotlib.use("MacOSX") + +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd +import scipy.io + +MAT_PATH = "/Users/evajain/Downloads/Molicel_INR18650P30B_241c01110_simulation.mat" + +mat = scipy.io.loadmat(MAT_PATH, squeeze_me=True, struct_as_record=False) +simulation = mat["simulation"] + +experiment = simulation.fu.PRO[0] + +current_profile = np.array(experiment.I, dtype=float) +meas_cell_voltage = np.array(experiment.V, dtype=float) +time = np.array(experiment.t, dtype=float) + +dt = float(np.mean(np.diff(time))) + +print("Loaded samples:", len(current_profile)) +print("dt:", dt) +print("Experiment name:", experiment.name) + + + +initial_SOC = 1 +target_end_SOC = 0 + +I_dis = np.clip(-current_profile, 0, None) +total_discharge_Ah = np.sum(I_dis) * dt / 3600.0 + +if total_discharge_Ah == 0: + required_capacity = 3.0 +else: + required_capacity = total_discharge_Ah / (initial_SOC - target_end_SOC) + +print("Total Discharge Ah:", total_discharge_Ah) +print("Required capacity Ah:", required_capacity) + +R = 8.31446261815324 +F = 96485.33212 + +V0 = 3.71 +C1 = 1.127469 +C2 = 6.96148085 +C3 = 0.05243311 +C4 = 0.01567795 + +# Tuned: original 0.0016 was ~20x too small for the observed voltage swings +R0 = 0.035 +KERNEL_LEN = 200 + + +def ocv_from_soc(soc, T_K=298.15): + eps = 1e-9 + soc_shift = soc - (0.1 ** 3) + denom = np.clip(1.0 - soc_shift + C4, eps, None) + numer = np.clip(C1 * soc_shift + C3, eps, None) + log_term = np.log(numer / denom) + return V0 + (C2 * (R * T_K / F) * log_term) + + +# ============================================================ +# SOC TRACKING +# ============================================================ + +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]) + +# ============================================================ +# BASE VOLTAGE MODEL +# ============================================================ + +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_kern = np.zeros((N, KERNEL_LEN), dtype=float) +for k in range(KERNEL_LEN): + if k == 0: + X_kern[:, k] = current_profile + else: + X_kern[k:, k] = current_profile[:-k] + +X = np.hstack([np.ones((N, 1)), X_kern]) + +# Weight by current magnitude; use ALL samples +I_abs = np.abs(current_profile) +max_current = np.max(I_abs) if np.max(I_abs) > 0 else 1.0 +weights = 1.0 + 5.0 * (I_abs / max_current) + +W = np.sqrt(weights) +Xw = X * W[:, None] +yw = residual_target * W + +lam = 1e-3 +A = Xw.T @ Xw + lam * np.eye(KERNEL_LEN + 1) +b_vec = Xw.T @ yw + +params = np.linalg.solve(A, b_vec) + +bias_term = params[0] +learned_kernel = params[1:] + +print("Bias term:", bias_term) +print("Learned kernel length:", len(learned_kernel)) + +# ============================================================ +# SAVE 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") + +# ============================================================ +# MODEL CLASS +# ============================================================ + + +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): + return ocv_from_soc(soc, T_K) + + 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)) + + # Newest current at index 0 + self.I_hist = np.roll(self.I_hist, 1) + self.I_hist[0] = I + + # kernel[k] * I[t-k], no reversed index needed + h = float(np.dot(self.kernel, self.I_hist)) + 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) + +# ============================================================ +# ERROR METRICS +# ============================================================ + +rmse = np.sqrt(np.mean((voltage_log - meas_cell_voltage) ** 2)) +mae = np.mean(np.abs(voltage_log - meas_cell_voltage)) + +print(f"RMSE: {rmse:.4f} V") +print(f"MAE: {mae:.4f} V") + + +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() diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_8.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_8.py new file mode 100644 index 0000000..ea0597e --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_8.py @@ -0,0 +1,174 @@ +import polars as pl +import numpy as np +import matplotlib.pyplot as plt +from scipy.optimize import curve_fit, differential_evolution +from scipy.interpolate import cubic_spline + +import scipy.io +from scipy.integrate import cumulative_simpson +from scipy.signal import convolve +from functools import reduce + +matData = scipy.io.loadmat("../fs-data/FS-4/P30B Cell Data/Molicel_INR18650P30B_measurement.mat", squeeze_me=True, struct_as_record=False) +measurements = matData["measurement"] + +# (measurements.fu.DCC[1].t[1:]-measurements.fu.DCC[1].t[:-1]).mean() + +fullData = pl.DataFrame(schema={ + "I":pl.Float64, + "V":pl.Float64, + "t":pl.Float64, + "T_surf1":pl.Float64, + "T_surf2":pl.Float64, + "SOC":pl.Float64, +}) + +def func_SOC(I, t): + return np.concatenate([np.array([3.0]), 3.0 - cumulative_simpson(y=I, x=t)/3600.0]) + +for sim in np.concatenate([measurements.fu.DCC, measurements.fu.CHC, measurements.fu.DCP, measurements.fu.CHP]): + I = sim.I + V = sim.V + t = sim.t + T_surf = sim.T_surf + T_surf1 = T_surf[0, :] + T_surf2 = T_surf[1, :] + SOC = func_SOC(I, t) + fullData = fullData.vstack(pl.DataFrame({ + "I": I, + "V": V, + "t": t, + "T_surf1": T_surf1, + "T_surf2": T_surf2, + "SOC": SOC, + })) + +delta = 0.01 +kernel_duration = 30 +kernel_size = int(kernel_duration / delta) + +def func_V_h (SOC, M, I, kernel): + padded_I = np.pad(I, (kernel_size - 1, 0), mode='constant') + convolution = convolve(padded_I, kernel[::-1], mode='valid') + return M * SOC * convolution + +def func_V_X (tauX, RX, I): + VX_arr = np.zeros_like(I) + # reduce(lambda prev, i: VX_arr.__setitem__(i, VX_arr[i - 1] * np.exp(-delta / tauX) + (RX * I[i] * (1 - np.exp(-delta / tauX)))) or VX_arr[i], range(len(I)), None) + for i in range(len(I)): + if i == 0: + VX_arr[i] = 0 + else: + VX_arr[i] = VX_arr[i - 1] * np.exp(-delta / tauX) + (RX * I[i] * (1 - np.exp(-delta / tauX))) + return VX_arr + +def func_V_OCV (SOC, T, a0, a1, a2, a3, a4, a5, a6, K_T): + return (a0 + a1 * SOC + a2 * SOC**2 + a3 * SOC**3 + a4 * SOC**4 + a5 * SOC**5 + a6 * SOC**6) * (1 + K_T * (T - 25)) + +def terminal_voltage(SOC, T, I, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, a0, a1, a2, a3, a4, a5, a6, K_T): + tau1 = R1 * C1 + tau2 = R2 * C2 + kernal = np.exp(-np.arange(0, kernel_size) * delta / tau_H) / np.sum(np.exp(-np.arange(0, kernel_size) * delta / tau_H)) + V_h = func_V_h(SOC, M, I, kernal) + V_X1 = func_V_X(tau1, R1, I) + V_X2 = func_V_X(tau2, R2, I) + V_OCV = func_V_OCV(SOC, T, a0, a1, a2, a3, a4, a5, a6, K_T) + V_R0 = np.where(I >= 0, R0_discharge * I, R0_charge * I) + V_terminal = V_OCV - V_R0 + V_X1 + V_X2 + V_h + # print(f"V_OCV: {V_OCV[-1]}, V_R0: {V_R0[-1]}, V_X1: {V_X1[-1]}, V_X2: {V_X2[-1]}, V_h: {V_h[-1]}, V_terminal: {V_terminal[-1]}") + # print(f"Params: R0_discharge: {R0_discharge}, R0_charge: {R0_charge}, R1: {R1}, C1: {C1}, R2: {R2}, C2: {C2}, tau_H: {tau_H}, M: {M}, a0: {a0}, a1: {a1}, a2: {a2}, a3: {a3}, a4: {a4}, a5: {a5}, a6: {a6}, K_T: {K_T}") + print(f"Params: R0_discharge: {R0_discharge}, R0_charge: {R0_charge}, R1: {R1}, C1: {C1}, R2: {R2}, C2: {C2}, tau_H: {tau_H}, M: {M}") + if len(V_terminal) == len(fullData["V"].to_numpy()): + print(f"Mean Squared Error: {np.mean((V_terminal - fullData['V'].to_numpy()) ** 2)}") + return V_terminal + +def wrapper(X, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, a0, a1, a2, a3, a4, a5, a6, K_T): + SOC = X[0] + T = X[1] + I = X[2] + return terminal_voltage(SOC, T, I, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, a0, a1, a2, a3, a4, a5, a6, K_T) + +args = curve_fit(wrapper, (fullData["SOC"].to_numpy(), fullData["T_surf1"].to_numpy(), fullData["I"].to_numpy()), fullData["V"].to_numpy(), p0=[0.01, 0.01, 1000, 0.01, 0.01, 1000, 10, 1, 3.7, -0.5, 0.1, -0.01, 0.001, -0.0001, 0.01, 0.01], maxfev=10000) +R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, a0, a1, a2, a3, a4, a5, a6, K_T = args[0] + +run = measurements.fu.DCC[6] +fig = plt.figure() +plt.plot(run.t, run.V, label="Measured Voltage") +plt.plot(run.t, wrapper((func_SOC(run.I, run.t), run.T_surf[0, :], run.I), *args[0]), label="Fitted Voltage") +plt.xlabel("Time (s)") +plt.ylabel("Voltage (V)") +plt.legend() +plt.show() + +measurements.fu.DCC[0].name + +def buildDF(measurement_set): + df = pl.DataFrame(schema={ + "I":pl.Float64, + "V":pl.Float64, + "t":pl.Float64, + "T_surf1":pl.Float64, + "T_surf2":pl.Float64, + "SOC":pl.Float64, + }) + for sim in measurement_set: + I = sim.I + V = sim.V + t = sim.t + T_surf = sim.T_surf + T_surf1 = T_surf[0, :] + T_surf2 = T_surf[1, :] + SOC = func_SOC(I, t) + df = df.vstack(pl.DataFrame({ + "I": I, + "V": V, + "t": t, + "T_surf1": T_surf1, + "T_surf2": T_surf2, + "SOC": SOC, + })) + return df + +lowCurrMeasurements = [m for m in measurements.fu.DCC if "-0.1" in m.name] +lowCurrentDF = buildDF(lowCurrMeasurements) + +def OCV_terminal_voltage(SOC, T, a0, a1, a2, a3, a4, a5, a6, K_T): + return (a0 + a1 * SOC + a2 * SOC**2 + a3 * SOC**3 + a4 * SOC**4 + a5 * SOC**5 + a6 * SOC**6) * (1 + K_T * (T - 25)) + + +ocv_args = curve_fit(lambda X, a0, a1, a2, a3, a4, a5, a6, K_T: OCV_terminal_voltage(X[0], X[1], a0, a1, a2, a3, a4, a5, a6, K_T), (lowCurrentDF["SOC"].to_numpy(), lowCurrentDF["T_surf1"].to_numpy()), lowCurrentDF["V"].to_numpy(), p0=[3.7, -0.5, 0.1, -0.01, 0.001, -0.0001, 0.01, 0.01], maxfev=10000) +ocv_args[0] + +# OCV Args +# array([-1.23630823e+02, 1.87749488e+02, -1.14010230e+02, 3.67068235e+01, +# -6.61877171e+00, 6.33739870e-01, -2.51814884e-02, 5.90437546e-04]) + +lowCurrMeasurements + +for run in lowCurrMeasurements: + fig = plt.figure() + plt.plot(run.t, run.V, label="Measured Voltage") + plt.plot(run.t, OCV_terminal_voltage(func_SOC(run.I, run.t), run.T_surf[0, :], *ocv_args[0]), label="Fitted OCV Voltage") + plt.xlabel("Time (s)") + plt.ylabel("Voltage (V)") + plt.legend() + plt.show() + +# Params: R0_discharge: -1.3235393250096301, +# R0_charge: 0.15299343818330174, R1: 0.16123612685040048, +# C1: 0.034272989366562326, R2: -126.42453895480814, C2: 1091372.6474212944, +# tau_H: 4.7170565122757395, M: 0.0037977786728519702, + +def wrapper_post_OCV(X, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M): + SOC = X[0] + T = X[1] + I = X[2] + a0, a1, a2, a3, a4, a5, a6, K_T = ocv_args[0] + return terminal_voltage(SOC, T, I, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, a0, a1, a2, a3, a4, a5, a6, K_T) + +bounds = [ + (0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 1, 0.001), # Lower bounds for R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M + (2, 0.5, 1, 1, 1, 1, 100, 1), # Upper bounds for R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M +] + +args2 = curve_fit(wrapper_post_OCV, (fullData["SOC"].to_numpy(), fullData["T_surf1"].to_numpy(), fullData["I"].to_numpy()), fullData["V"].to_numpy(), p0=[0.01, 0.01, 0.16, 0.03, 0.01, 0.01, 10, 0.003], bounds=bounds, maxfev=10000) diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_9.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_9.py new file mode 100644 index 0000000..dbbcfd7 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_model_9.py @@ -0,0 +1,426 @@ +import polars as pl +import numpy as np +import matplotlib.pyplot as plt +from scipy.optimize import curve_fit, differential_evolution +from scipy.integrate import cumulative_simpson +from scipy.signal import convolve +import scipy.io +from numba import njit +from scipy.interpolate import LinearNDInterpolator as Lnd_interp +from scipy.interpolate import CubicSpline as cs +from functools import partial + + +# 1. Load Data +matData = scipy.io.loadmat("../fs-data/FS-4/P30B Cell Data/Molicel_INR18650P30B_measurement.mat", squeeze_me=True, struct_as_record=False) +measurements = matData["measurement"] + +delta = 0.1 +kernel_duration = 60 +kernel_size = int(kernel_duration / delta) + +lowCurrMeasurements = [m for m in measurements.fu.DCC if "-0.1" in m.name] + +def func_SOC(I, t, starting_SOC): + Q_nominal = 3.0 # Ah capacity of Molicel P30B + integrated_ah = cumulative_simpson(y=I, x=t, initial=0) / 3600.0 + soc = starting_SOC + (integrated_ah / Q_nominal) + return np.clip(soc, 0.0, 1.0) + +lowVs = np.concatenate([sim.V for sim in lowCurrMeasurements]) +lowTs = np.concatenate([sim.T_surf[0, :] for sim in lowCurrMeasurements]) +lowSOC = np.concatenate([func_SOC(sim.I, sim.t, 1.0) for sim in lowCurrMeasurements]) + +voltage_curve = Lnd_interp(np.array([lowVs, lowTs]).T, lowSOC, rescale=True) + +def SOC_lookup(V, T): + SOC = voltage_curve(V, T) + if np.isnan(SOC): + if V > 4.1: + return 1.0 + if T < -15.0: + return voltage_curve(V, -15.0) + if T > 39.0: + return voltage_curve(V, 39.0) + return np.clip(SOC, 0.0, 1.0) +# Fixed SOC calculation tracking 0.0 to 1.0 bounds safely + +# plt.plot(lowCurrMeasurements[0].t, func_SOC(lowCurrMeasurements[0].I, lowCurrMeasurements[0].t, SOC_lookup(lowCurrMeasurements[0].V[0], lowCurrMeasurements[0].T_surf[0, 0]))) +# lowCurrMeasurements[0].V[0], lowCurrMeasurements[0].T_surf[0, 0] + +## Validation of SOC!! +# X, Y = np.meshgrid(np.arange(2.5, 4.2, 0.01), np.arange(-20, 45, 1)) +# Z = SOC_lookup(X, Y) +# plt.contourf(X, Y, Z, levels=50, cmap='viridis') +# plt.colorbar(label='SOC') +# plt.xlabel('Voltage (V)') +# plt.ylabel('Temperature (°C)') +# plt.title('SOC Lookup Table') +# plt.show() + +# measurements.fu.DCC[1].V[0] +# measurements.fu.DCC[1].T_surf[0, 0] +# SOC_lookup(4.169, -19.37) + +def buildDF(measurement_set): + frames = [] + for i, sim in enumerate(measurement_set): + cubic = cs(sim.t, np.array([sim.I, sim.V, sim.T_surf[0, :], sim.T_surf[1, :]]), axis=1) + t = np.arange(sim.t[0], sim.t[-1], delta) + I, V, T_surf1, T_surf2 = cubic(t) + starting_SOC = SOC_lookup(V[0], T_surf1[0]) + SOC = func_SOC(I, t, starting_SOC) + frames.append(pl.DataFrame({ + "I": I, "V": V, "t": t, + "T_surf1": T_surf1, "T_surf2": T_surf2, "SOC": SOC, + "run": np.ones(len(I)) * i + })) + return pl.concat(frames) + +# lowCurrentDF = buildDF(lowCurrMeasurements) +# Gather Datasets + + + +# Use everything for dynamic fitting +all_sims = np.concatenate([measurements.fu.DCC, measurements.fu.CHC, measurements.fu.DCP, measurements.fu.CHP]) +# for sim in all_sims: +# startingSOC = SOC_lookup(sim.V[0], sim.T_surf[0, 0]) +# print(f"Run: {sim.name}, Starting SOC: {startingSOC:.4f}") + +aboveFreezingSims = [m for m in all_sims if m.T_surf[0, 0] > 5] +notHighPulse = [m for m in aboveFreezingSims if not ("-40.000C" in m.name or "-30.000C" in m.name)] +# for sim in notHighPulse: + # print(sim.name) +fullData = buildDF(notHighPulse) + +# 2. Extract OCV +def OCV_terminal_voltage(SOC, T, a0, a1, a2, a3, a4, a5, a6, K_T): + return (a0 + a1 * SOC + a2 * SOC**2 + a3 * SOC**3 + a4 * SOC**4 + a5 * SOC**5 + a6 * SOC**6) * (1 + K_T * (T - 25)) + +# OCV Args +# array([ 2.33727919e+00, 8.66472120e+00, -3.23286606e+01, 7.50959382e+01, +# -9.75689762e+01, 6.62872924e+01, -1.83568111e+01, 5.90434471e-04]) + +# ocv_args = np.array([ 2.33727919e+00, 8.66472120e+00, -3.23286606e+01, 7.50959382e+01, -9.75689762e+01, 6.62872924e+01, -1.83568111e+01, 5.90434471e-04]) +ocv_args = np.array([ 2.33586300e+00, 8.70615991e+00, -3.26789721e+01, 7.63836093e+01, -9.98902031e+01, 6.83036001e+01, -1.90310137e+01, 5.90707089e-04]) + +# ocv_args, _ = curve_fit( +# lambda X, a0, a1, a2, a3, a4, a5, a6, K_T: OCV_terminal_voltage(X[0], X[1], a0, a1, a2, a3, a4, a5, a6, K_T), +# (lowCurrentDF["SOC"].to_numpy(), lowCurrentDF["T_surf1"].to_numpy()), +# lowCurrentDF["V"].to_numpy(), +# p0=[3.7, 0.5, -0.1, 0.01, -0.001, 0.0, 0.0, 0.001], +# maxfev=20000 +# ) + +# 3. Dynamic Optimization Functions (Accelerated) +# @njit +def func_V_h(SOC, M, I, kernel): + padded_I = np.pad(I, (kernel_size - 1, 0), mode='constant') + convolution = convolve(padded_I, kernel, mode='valid') + # return M * convolution + return M * SOC * convolution + + +@njit +def func_V_X_fast(tauX, RX, I, delta): + VX_arr = np.zeros_like(I) + alpha = np.exp(-delta / tauX) + beta = RX * (1.0 - alpha) + for i in range(1, len(I)): + VX_arr[i] = VX_arr[i - 1] * alpha + I[i] * beta + return VX_arr + +example_tests = ["DCC / 25°C / -10.000C", "CHC / 25°C / 3.000C", "DCP / 25°C / -30.000C", "CHP / 40°C / 4.000C"] +example_runs = [m for m in measurements.fu.DCC if m.name in example_tests] + \ + [m for m in measurements.fu.CHC if m.name in example_tests] + \ + [m for m in measurements.fu.DCP if m.name in example_tests] + \ + [m for m in measurements.fu.CHP if m.name in example_tests] + +# DCC / 25°C / -10.000C +# CHC / 25°C / 3.000C +# DCP / 25°C / -30.000C +# CHP / 40°C / 4.000C + +def terminal_voltage(SOC, T, I, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, ocv_params, piece=False, verbose=False, superVerbose=True): + tau1 = R1 * C1 + tau2 = R2 * C2 + + # Kernel calculation + exponent = -np.arange(0, kernel_size) * delta / tau_H + kernel = np.exp(exponent) / np.sum(np.exp(exponent)) + + V_h = func_V_h(SOC, M, I, kernel) + V_X1 = func_V_X_fast(tau1, R1, I, delta) + V_X2 = func_V_X_fast(tau2, R2, I, delta) + + V_OCV = OCV_terminal_voltage(SOC, T, *ocv_params) + V_R0 = np.where(I <= 0, R0_discharge * I, R0_charge * I) + # if superVerbose: + # print(f"R0_discharge: {R0_discharge:.5f}, R0_charge: {R0_charge:.5f}, R1: {R1:.5f}, C1: {C1:.5f}, R2: {R2:.5f}, C2: {C2:.5f}, tau_H: {tau_H:.5f}, M: {M:.5f}") + if not piece: + print(f"MAE: {np.mean(np.abs(V_OCV + V_R0 + V_X1 + V_X2 + V_h - fullData['V'].to_numpy()))}") + # print(f"MSE: {np.mean((V_OCV + V_R0 + V_X1 + V_X2 + V_h - fullData['V'].to_numpy()) ** 2)}") + if verbose: + print(f"V_OCV: {V_OCV[-1]}, V_R0: {V_R0[-1]}, V_X1: {V_X1[-1]}, V_X2: {V_X2[-1]}, V_h: {V_h[-1]}, V_terminal: {V_OCV[-1] + V_R0[-1] + V_X1[-1] + V_X2[-1] + V_h[-1]}") + + return V_OCV + V_R0 + V_X1 + V_X2 + V_h + +def debug_terminal_voltage(SOC, T, I, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, ocv_params, piece=False, verbose=False, superVerbose=True): + tau1 = R1 * C1 + tau2 = R2 * C2 + + # Kernel calculation + exponent = -np.arange(0, kernel_size) * delta / tau_H + kernel = np.exp(exponent) / np.sum(np.exp(exponent)) + + V_h = func_V_h(SOC, M, I, kernel) + V_X1 = func_V_X_fast(tau1, R1, I, delta) + V_X2 = func_V_X_fast(tau2, R2, I, delta) + + V_OCV = OCV_terminal_voltage(SOC, T, *ocv_params) + V_R0 = np.where(I <= 0, R0_discharge * I, R0_charge * I) + if superVerbose: + print(f"R0_discharge: {R0_discharge:.5f}, R0_charge: {R0_charge:.5f}, R1: {R1:.5f}, C1: {C1:.5f}, R2: {R2:.5f}, C2: {C2:.5f}, tau_H: {tau_H:.5f}, M: {M:.5f}") + if not piece: + print(f"MAE: {np.mean(np.abs(V_OCV + V_R0 + V_X1 + V_X2 + V_h - fullData['V'].to_numpy()))}") + # print(f"MSE: {np.mean((V_OCV + V_R0 + V_X1 + V_X2 + V_h - fullData['V'].to_numpy()) ** 2)}") + if verbose: + print(f"V_OCV: {V_OCV[-1]}, V_R0: {V_R0[-1]}, V_X1: {V_X1[-1]}, V_X2: {V_X2[-1]}, V_h: {V_h[-1]}, V_terminal: {V_OCV[-1] + V_R0[-1] + V_X1[-1] + V_X2[-1] + V_h[-1]}") + + return V_OCV, V_R0, V_X1, V_X2, V_h + + +# fig = plt.figure() +# axes = fig.subplots(2, 2) +# fig.show() + +def wrapper_post_OCV(X, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M): + # print(X.shape) + # for ax, run in zip(axes.flatten(), example_runs): + # run_soc = func_SOC(run.I, run.t, SOC_lookup(run.V[0], run.T_surf[0, 0])) + # fitted_v = terminal_voltage(run_soc, run.T_surf[0, :], run.I, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, ocv_args, piece=True) + # ax.clear() + # ax.plot(run.t, run.V, label="Measured Voltage", alpha=0.8) + # ax.plot(run.t, fitted_v, '--', label="Fitted Performance ECM", alpha=0.9) + # ax.set_xlabel("Time (s)") + # ax.set_ylabel("Voltage (V)") + # ax.set_title(f"Run: {run.name}") + # ax.legend() + # ax.grid(True) + # fig.canvas.draw() + # fig.canvas.flush_events() + modelVoltages = [] + for id in np.unique(X[3]): + mask = X[3] == id + modelVoltage = terminal_voltage(X[0, mask], X[1, mask], X[2, mask], R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, ocv_args, piece=True) + measuredVoltage = X[4, mask] + # print(modelVoltage.shape, measuredVoltage.shape) + # print(f"Run ID: {id}, MAE: {np.mean(np.abs(modelVoltage - measuredVoltage))}") + modelVoltages.append(modelVoltage) + print(f"Overall MAE: {np.mean(np.abs(np.concatenate(modelVoltages) - X[4]))}") + return np.concatenate(modelVoltages) + +def wrapper_pre_OCV(X, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, a0, a1, a2, a3, a4, a5, a6, K_T): + SOC = X[0] + T = X[1] + I = X[2] + modelVoltages = [] + for id in np.unique(X[3]): + mask = X[3] == id + modelVoltage = terminal_voltage(X[0, mask], X[1, mask], X[2, mask], R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, (a0, a1, a2, a3, a4, a5, a6, K_T), piece=True) + modelVoltages.append(modelVoltage) + # measuredVoltage = X[4, mask] + # print(f"Run ID: {id}, MAE: {np.mean(np.abs(modelVoltage - measuredVoltage))}") + print(f"Overall MAE: {np.mean(np.abs(np.concatenate(modelVoltages) - X[4]))}") + print(f"R0_discharge: {R0_discharge:.3f}, R0_charge: {R0_charge:.3f}, R1: {R1:.3f}, C1: {C1:.3f}, R2: {R2:.3f}, C2: {C2:.3f}, tau_H: {tau_H:.3f}, M: {M:.3f}, a0: {a0:.3f}, a1: {a1:.3f}, a2: {a2:.3f}, a3: {a3:.3f}, a4: {a4:.3f}, a3: {a5:.3f}, a6: {a6:.3f}, K_T: {K_T:.3f}") + return np.concatenate(modelVoltages) + +def evolution_wrapper(params, X): + R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M = params + v = wrapper_post_OCV(X, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M) + MSE = np.mean((v - X[4]) ** 2) + return MSE + + +# X = np.array([fullData["SOC"].to_numpy(), fullData["T_surf1"].to_numpy(), fullData["I"].to_numpy(), fullData["run"].to_numpy(), fullData["V"].to_numpy()]) +# np.unique(X[3]) +# X[:, np.where(X[3] == 0)] + +# np.where(X[3] == 0) +# currentRun.shape +# currentRun = X[np.where(X[3] == id)] +# X + +# R0_discharge: 0.00010000241368768962, +# R0_charge: 0.0015092018205301212, +# R1: 0.008623195393932841, +# C1: 37.86392127904164, +# R2: 0.013221128224458535, +# C2: 761.6878322739899, +# tau_H: 53.00341711702909, +# M: 0.0018520981883959234 + +# Correct physical boundaries +bounds = ( + [1e-7, 1e-7, 1e-7, 1e-3, 1e-7, 10.0, 1e-5, -1], # Lower Bounds + [0.2, 0.2, 0.2, 5000.0, 0.2, 50000.0, 1e4, 1e3] # Upper Bounds +) + +boundsTuppled = tuple(zip(*bounds)) + +# Initial dynamic guesses +p0 = [0.002728, 0.00010000, 0.00757, 10.0, 0.011, 735, 0.833, 0.00010000] + +p0 = [ 7.54812179e-03, 1.00000000e-07, 2.44930101e-02, 5.74870831e+01, + 6.77082684e-04, 4.41486533e+04, 2.97718020e+03, -5.27874313e-03] + +# Conditions generated post hysteresis function fix. +p0 = [1.00000000e-07, 1.00000000e-07, 1.40292451e-02, 3.15517635e+03, + 1.51833942e-03, 2.91537408e+04, 1.5, 2.86674119e-03] + +preOCV_p0 = p0 + list(ocv_args) + +lessData = fullData.filter(pl.col("t") < 1000) + +args = np.array([5.30000214e-03, 9.20377540e-03, 8.07481539e-03, 1.41824989e+03, + 3.29520249e-03, 2.59134546e+02, 1.5, 1e-03, + 2.39587985e+00, 1.27417190e+01, -6.73434102e+01, 1.89554380e+02, + -2.78887472e+02, 2.05193796e+02, -5.95230431e+01, 9.40568639e-05]) + +args = np.array([1.21786569e-02, 1.76761209e-02, 1.0, 1.0, + 1.0, 1.0, 9.62831610e+00, 1.03801131e-03, + 2.39250764e+00, 1.28854012e+01, -6.78772374e+01, 1.89652101e+02, + -2.76828230e+02, 2.02226116e+02, -5.83091863e+01, -8.11863308e-05]) + +args = np.array([ 9.97263998e-03, 5.00002331e-03, 1.60732764e-02, 0.01, + 7.31018241e-03, 5.72244145e+02, 1.00000000e+00, 1.00000000e-07, + 2.38962033e+00, 1.29375218e+01, -6.85304383e+01, 1.92513553e+02, + -2.82364531e+02, 2.07151056e+02, -5.99607550e+01, -1.23234655e-05]) + +args = np.array([0.006, 0.011, 0.001, 0.010, 0.006, 792.945, 7.433, 0.000001, 2.390, 12.938, -68.531, 192.513, -282.365, 207.151, -59.960, -0.000001]) + +conservative_bounds = ( +# R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, a0, a1, a2, a3, a4, a5, a6, K_T + [0.005, 0.005, 0.0, 0.0, 0.0, 0.0, 1.0, 1e-7, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf], # Lower Bounds + [0.02, 0.02, 0.02, np.inf, 0.02, np.inf, 20.0, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf] # Upper Bounds +) + +# Optimize dynamic parameters +args2, _ = curve_fit( + wrapper_post_OCV, + np.array((lessData["SOC"].to_numpy(), lessData["T_surf1"].to_numpy(), lessData["I"].to_numpy(), lessData["run"].to_numpy(), lessData["V"].to_numpy())), + lessData["V"].to_numpy(), + # p0=p0, + # bounds=bounds, + maxfev=15000, + method='trf', + x_scale='jac' +) + +args3, _ = curve_fit( + wrapper_pre_OCV, + np.array((lessData["SOC"].to_numpy(), lessData["T_surf1"].to_numpy(), lessData["I"].to_numpy(), lessData["run"].to_numpy(), lessData["V"].to_numpy())), + lessData["V"].to_numpy(), + p0=args, + bounds=conservative_bounds, + maxfev=1500, + method='trf', + x_scale='jac' +) + + +np.mean(np.abs(lessData["V"].to_numpy() - wrapper_pre_OCV(np.array((lessData["SOC"].to_numpy(), lessData["T_surf1"].to_numpy(), lessData["I"].to_numpy(), lessData["run"].to_numpy(), lessData["V"].to_numpy())), *args))) + +args3 + +# argsGlobal = differential_evolution( +# evolution_wrapper, +# boundsTuppled, +# args=(np.array((lessData["SOC"].to_numpy(), lessData["T_surf1"].to_numpy(), lessData["I"].to_numpy(), lessData["run"].to_numpy(), lessData["V"].to_numpy())),), +# maxiter=100, +# popsize=15, +# workers=-1, +# updating='deferred', +# ) + +# R0_discharge = -0.014174362096127368; R0_charge = -0.023739566602921935; R1 = 327.6950042682508; C1 = 2134.853515766918; R2 = -171.55438396480167; C2 = -2132.612618294229; tau_H = 0.8719469562101347; M = -0.00022039527420898293 + +tau_H = 1.0 +exponent = -np.arange(0, kernel_size) * delta / tau_H +kernel = np.exp(exponent) / np.sum(np.exp(exponent)) +# plt.plot(kernel) +# plt.show() + +# Make sure the kernel faces the right way +plt.plot(convolve(np.pad(np.concatenate([np.ones(kernel_size), np.ones(kernel_size)*-1]), (kernel_size-1, 0), mode='constant'), kernel, mode="valid"), label="Convolution Result") +plt.plot(np.concatenate([np.ones(kernel_size), np.ones(kernel_size)*-1]), label="Input Current") +plt.legend() +plt.show() + +arr = np.array((lessData["SOC"].to_numpy(), lessData["T_surf1"].to_numpy(), lessData["I"].to_numpy(), lessData["run"].to_numpy(), lessData["V"].to_numpy(), lessData["t"].to_numpy())) +for id, name in zip(np.unique(arr[3]), notHighPulse): + currentRun = arr[:, np.where(arr[3] == id)][:, 0, :] + V_OCV, V_R0, V_X1, V_X2, V_h = debug_terminal_voltage(currentRun[0], currentRun[1], currentRun[2], *args[:6], 10, args[7], args[8:], piece=True, verbose=True, superVerbose=False) + measuredVoltage = currentRun[4] + print(f"Run ID: {id}, MAE: {np.mean(np.abs(V_OCV + V_R0 + V_X1 + V_X2 + V_h - measuredVoltage))}") + plt.plot(currentRun[5], measuredVoltage, label="Measured Voltage", alpha=0.8) + plt.plot(currentRun[5], V_OCV + V_R0 + V_X2 + V_h, '--', label="Fitted Performance ECM", alpha=0.9) + plt.plot(currentRun[5], V_OCV, label="OCV", alpha=0.8) + plt.plot(currentRun[5], V_R0, label="R0 Voltage Drop", alpha=0.8) + plt.plot(currentRun[5], V_X1, label="X1 Voltage", alpha=0.8) + plt.plot(currentRun[5], V_X2, label="X2 Voltage", alpha=0.8) + plt.plot(currentRun[5], V_h, label="Hysteresis Voltage", alpha=0.8) + plt.xlabel("Time (s)") + plt.ylabel("Voltage (V)") + plt.title(name.name) + plt.legend() + plt.grid(True) + plt.show() + + + + + + + + + + +# print("Optimized Parameters (R0_d, R0_c, R1, C1, R2, C2, tau_H, M):") +# print(args2) + +# # 4. Plot Verification Validating Against One Dynamic Run +# run = measurements.fu.DCC[6] +# run_soc = func_SOC(run.I, run.t, SOC_lookup(run.V[0], run.T_surf[0, 0])) +# fitted_v = terminal_voltage(run_soc, run.T_surf[0, :], run.I, args2[0], args2[1], args2[2], args2[3], args2[4], args2[5], args2[6], args2[7], ocv_args, piece=True) + +# np.abs(terminal_voltage(run_soc, run.T_surf[0, :], run.I, args2[0], args2[1], args2[2], args2[3], args2[4], args2[5], args2[6], args2[7], ocv_args, piece=True) - run.V).mean() + + +# plt.figure(figsize=(10, 5)) +# plt.plot(run.t, run.V, label="Measured Voltage", alpha=0.8) +# plt.plot(run.t, fitted_v, '--', label="Fitted Performance ECM", alpha=0.9) +# plt.xlabel("Time (s)") +# plt.ylabel("Voltage (V)") +# plt.title(f"Model Validation Execution - Run: {run.name}") +# plt.legend() +# plt.grid(True) +# plt.show() + +# plt.plot(fullData["SOC"].to_numpy()) +# plt.show() + +# for run in aboveFreezingSims: +# plt.plot(func_SOC(run.I, run.t, SOC_lookup(run.V[0], run.T_surf[0, 0]))) +# plt.show() + +# for run in aboveFreezingSims: +# plt.plot(run.t, run.V, label="Measured Voltage", alpha=0.8) +# plt.plot(run.t, terminal_voltage(func_SOC(run.I, run.t, SOC_lookup(run.V[0], run.T_surf[0, 0])), run.T_surf[0, :], run.I, p0[0], p0[1], p0[2], p0[3], p0[4], p0[5], p0[6], p0[7], ocv_args, piece=True), '--', label="SOC", alpha=0.9) +# plt.xlabel("Time (s)") +# plt.ylabel("Voltage (V)") +# plt.title(f"Run: {run.name}") +# plt.legend() +# plt.grid(True) +# plt.show() \ No newline at end of file diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electrical_modelingG12.py b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_modelingG12.py new file mode 100644 index 0000000..2922ede --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electrical_modelingG12.py @@ -0,0 +1,179 @@ +import polars as pl +import numpy as np +import matplotlib.pyplot as plt +from scipy.optimize import curve_fit +from scipy.integrate import cumulative_simpson +from scipy.signal import convolve +import scipy.io +from numba import njit +from scipy.interpolate import LinearNDInterpolator as Lnd_interp +from scipy.interpolate import CubicSpline as cs + +# 1. Load Data +matData = scipy.io.loadmat("../fs-data/FS-4/P30B Cell Data/Molicel_INR18650P30B_measurement.mat", squeeze_me=True, struct_as_record=False) +measurements = matData["measurement"] + +delta = 0.1 +kernel_duration = 60 +kernel_size = int(kernel_duration / delta) + +lowCurrMeasurements = [m for m in measurements.fu.DCC if "-0.1" in m.name] + +def func_SOC(I, t, starting_SOC): + Q_nominal = 3.0 # Ah capacity of Molicel P30B + integrated_ah = cumulative_simpson(y=I, x=t, initial=0) / 3600.0 + soc = starting_SOC + (integrated_ah / Q_nominal) + return np.clip(soc, 0.0, 1.0) + +lowVs = np.concatenate([sim.V for sim in lowCurrMeasurements]) +lowTs = np.concatenate([sim.T_surf[0, :] for sim in lowCurrMeasurements]) +lowSOC = np.concatenate([func_SOC(sim.I, sim.t, 1.0) for sim in lowCurrMeasurements]) + +voltage_curve = Lnd_interp(np.array([lowVs, lowTs]).T, lowSOC, rescale=True) + +def SOC_lookup(V, T): + SOC = voltage_curve(V, T) + if np.isnan(SOC): + if V > 4.1: return 1.0 + if T < -15.0: return voltage_curve(V, -15.0) + if T > 39.0: return voltage_curve(V, 39.0) + return np.clip(SOC, 0.0, 1.0) + +def buildDF(measurement_set): + frames = [] + for i, sim in enumerate(measurement_set): + cubic = cs(sim.t, np.array([sim.I, sim.V, sim.T_surf[0, :], sim.T_surf[1, :]]), axis=1) + t = np.arange(sim.t, sim.t[-1], delta) + I, V, T_surf1, T_surf2 = cubic(t) + starting_SOC = SOC_lookup(V, T_surf1) + SOC = func_SOC(I, t, starting_SOC) + frames.append(pl.DataFrame({ + "I": I, "V": V, "t": t, + "T_surf1": T_surf1, "T_surf2": T_surf2, "SOC": SOC, + "run": np.ones(len(I)) * i + })) + return pl.concat(frames) + +all_sims = np.concatenate([measurements.fu.DCC, measurements.fu.CHC, measurements.fu.DCP, measurements.fu.CHP]) +aboveFreezingSims = [m for m in all_sims if m.T_surf > -2] +notHighPulse = [m for m in aboveFreezingSims if not ("-40.000C" in m.name or "-30.000C" in m.name)] +fullData = buildDF(notHighPulse) + +# 2. Extract and Freeze OCV +def OCV_terminal_voltage(SOC, T, a0, a1, a2, a3, a4, a5, a6, K_T): + return (a0 + a1 * SOC + a2 * SOC**2 + a3 * SOC**3 + a4 * SOC**4 + a5 * SOC**5 + a6 * SOC**6) * (1 + K_T * (T - 25.0)) + +ocv_args = np.array([2.33586300e+00, 8.70615991e+00, -3.26789721e+01, 7.63836093e+01, + -9.98902031e+01, 6.83036001e+01, -1.90310137e+01, 5.90707089e-04]) + +# 3. Enhanced Hysteresis Tracking +def func_V_h(SOC, M, I, kernel): + current_direction = np.sign(I) + padded_I = np.pad(current_direction, (kernel_size - 1, 0), mode='constant') + convolution = convolve(padded_I, kernel, mode='valid') + return M * convolution + +@njit +def func_V_X_fast(tauX, RX, I, delta): + VX_arr = np.zeros_like(I) + alpha = np.exp(-delta / tauX) + beta = RX * (1.0 - alpha) + for i in range(1, len(I)): + VX_arr[i] = VX_arr[i - 1] * alpha + I[i] * beta + return VX_arr + +def terminal_voltage(SOC, T, I, R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, ocv_params): + tau1 = R1 * C1 + tau2 = R2 * C2 + + exponent = -np.arange(0, kernel_size) * delta / tau_H + kernel = np.exp(exponent) / np.sum(np.exp(exponent)) + + V_h = func_V_h(SOC, M, I, kernel) + V_X1 = func_V_X_fast(tau1, R1, I, delta) + V_X2 = func_V_X_fast(tau2, R2, I, delta) + + V_OCV = OCV_terminal_voltage(SOC, T, *ocv_params) + V_R0 = np.where(I <= 0, R0_discharge * I, R0_charge * I) + + return V_OCV + V_R0 + V_X1 + V_X2 + V_h + +# 4. Scaled Optimization Wrapper +def wrapper_post_OCV(X, R0_d_scaled, R0_c_scaled, R1_scaled, C1_scaled, R2_scaled, C2_scaled, tau_H, M_scaled): + # Scale variables back up to their physical magnitudes + R0_discharge = R0_d_scaled * 1e-3 + R0_charge = R0_c_scaled * 1e-3 + R1 = R1_scaled * 1e-3 + C1 = C1_scaled * 1e2 + R2 = R2_scaled * 1e-3 + C2 = C2_scaled * 1e4 + M = M_scaled * 1e-2 + + modelVoltages = [] + run_ids = X[3, :] + + for id in np.unique(run_ids): + mask = run_ids == id + modelVoltage = terminal_voltage( + X[0, mask], X[1, mask], X[2, mask], + R0_discharge, R0_charge, R1, C1, R2, C2, tau_H, M, ocv_args + ) + modelVoltages.append(modelVoltage) + + actual_voltages = X[4, :] + current_mae = np.mean(np.abs(np.concatenate(modelVoltages) - actual_voltages)) + print(f"Iteration MAE: {current_mae:.5f} V") + + return np.concatenate(modelVoltages) + +# 5. Scaled Parameter Boundaries & Guess Configuration +# Enforces strict parameter separation: +# Tau1 (R1*C1) bounds: 0.1s to 8s (Fast Dynamics) +# Tau2 (R2*C2) bounds: 15s to 300s (Slow Diffusion) + +# Parameters mapped: [R0_d, R0_c, R1, C1, R2, C2, tau_H, M] (all values normalized around ~1-10) +scaled_bounds = ( + [2.0, 2.0, 0.5, 1.0, 1.0, 0.5, 0.5, 0.01], # Lower Bounds + [40.0, 40.0, 20.0, 10.0, 30.0, 5.0, 120.0, 5.00] # Upper Bounds +) + +p0_scaled = [12.0, 15.0, 6.0, 4.0, 12.0, 1.5, 15.0, 0.8] + +# Downsample the whole time profile using polars instead of throwing away tail data +lessData = fullData.gather_every(10) + +X_data = np.array(( + lessData["SOC"].to_numpy(), + lessData["T_surf1"].to_numpy(), + lessData["I"].to_numpy(), + lessData["run"].to_numpy(), + lessData["V"].to_numpy() +)) + +print(f"Starting dynamic tuning sequence with {X_data.shape} data points...") +optimized_scaled_args, _ = curve_fit( + wrapper_post_OCV, + X_data, + lessData["V"].to_numpy(), + p0=p0_scaled, + bounds=scaled_bounds, + maxfev=3000, + method='trf', + x_scale='jac' +) + +# Unpack and rescale final optimized values back to true SI units +final_params = { + "R0_discharge": optimized_scaled_args * 1e-3, + "R0_charge": optimized_scaled_args * 1e-3, + "R1": optimized_scaled_args * 1e-3, + "C1": optimized_scaled_args * 1e2, + "R2": optimized_scaled_args * 1e-3, + "C2": optimized_scaled_args * 1e4, + "tau_H": optimized_scaled_args, + "M": optimized_scaled_args * 1e-2 +} + +print("\nOptimized Dynamic Parameters (Rescaled to Physical Units):") +for k, v in final_params.items(): + print(f"{k}: {v:.6e}") \ No newline at end of file diff --git a/FullVehicleSim/Electrical/HisteresisCellModel/electricl_model6.py b/FullVehicleSim/Electrical/HisteresisCellModel/electricl_model6.py new file mode 100644 index 0000000..745b683 --- /dev/null +++ b/FullVehicleSim/Electrical/HisteresisCellModel/electricl_model6.py @@ -0,0 +1,328 @@ +import numpy as np +import matplotlib.pyplot as plt +import polars as pl +import scipy.io +from scipy.integrate import simpson, cumulative_simpson +import matplotlib.patches as patches +from pysr import PySRRegressor, TemplateExpressionSpec + +MAT_PATH = "../fs-data/FS-4/P30B Cell Data/Molicel_INR18650P30B_measurement.mat" + + + +mat = scipy.io.loadmat(MAT_PATH, squeeze_me=True, struct_as_record=False) + +''' +header = batemoData["__header__"] +version = batemoData["__version__"] +globals = batemoData["__globals__"] ## Nothing +print(f"header = {header}") +print(f"version = {version}") +measurement = batemoData["measurement"] +firstLayerMeta = measurement['meta'] +Fu = measurement['fu'] + +# DCC (Discharge) +# CHC (Charging) +# DCP (Discharge Pulse) +# CHP (Charge Pulse) +PRO (Profile Measurement) + +Each Has + name + T_amb (Ambient Temperature) + t (Time Seconds) + I (Current) + V (Voltage) + T_surf (Surface temperature, nominally 1xN but may be 2xN?) +''' + +simulation = mat["measurement"] + +fullDCC = pl.DataFrame(schema={ + "I":pl.Float64, + "V":pl.Float64, + "t":pl.Float64, + "T_surf1":pl.Float64, + "T_surf2":pl.Float64, + "SOC":pl.Float64, + "I_Past1":pl.Float64, + "I_Past2":pl.Float64, + "I_Past3":pl.Float64, + "I_Past4":pl.Float64, +}) + +for dcc in simulation.fu.DCC: + I = dcc.I + I_back1 = np.concatenate([np.zeros(1), I[:-1]]) + I_back2 = np.concatenate([np.zeros(3), I[:-3]]) + I_back3 = np.concatenate([np.zeros(7), I[:-7]]) + I_back4 = np.concatenate([np.zeros(20), I[:-20]]) + V = dcc.V + t = dcc.t + T_surf = dcc.T_surf + T_surf1 = T_surf[0, :] + T_surf2 = T_surf[1, :] + SOC = 3.0 - cumulative_simpson(y=I, x=t)/3600.0 + SOC = np.concatenate([np.array([3.0]), SOC]) + fullDCC = fullDCC.vstack(pl.DataFrame({ + "I": I, + "V": V, + "t": t, + "T_surf1": T_surf1, + "T_surf2": T_surf2, + "SOC": SOC, + "I_Past1": I_back1, + "I_Past2": I_back2, + "I_Past3": I_back3, + "I_Past4": I_back4, + })) + +# simulation.fu.DCC +# len(simulation.fu.DCC) +# for dcc in simulation.fu.DCC: +# I = dcc.I +# V = dcc.V +# t = dcc.t +# T_amb = dcc.T_amb +# T_surf = dcc.T_surf +# name = dcc.name +# SOC = 3.0 - cumulative_simpson(y=I, x=t)/3600.0 +# SOC = np.concatenate([np.array([3.0]), SOC]) +# fig = plt.figure() +# ax1 = fig.add_subplot(111) +# ax2 = ax1.twinx() +# ax1.plot(SOC, I, label="Current", color="blue") +# ax2.plot(SOC, V, label="Voltage", color="orange") +# ## put patches for legends with current and voltage on the same legend with current blue nad voltage orange +# current_patch = patches.Patch(color='blue', label='Current [A]') +# voltage_patch = patches.Patch(color='orange', label='Voltage [V]') +# plt.legend(handles=[current_patch, voltage_patch]) +# plt.title(f"DCC: {name} at T_amb={T_amb} K") +# plt.xlabel("SOC") +# plt.grid(True) +# plt.show() + +# SOC = 3.0 - cumulative_simpson(y=simulation.fu.DCC[0].I, x=simulation.fu.DCC[0].t)/3600.0 +# SOC = np.concatenate([np.array([3.0]), SOC]) + +X = fullDCC.select(pl.col("I"), pl.col("SOC"), pl.col("T_surf1"), pl.col("T_surf2"), pl.col("I_Past1"), pl.col("I_Past2"), pl.col("I_Past3"), pl.col("I_Past4")).to_numpy() +y = fullDCC.select(pl.col("V")).to_numpy() + +# "I": I, +# "V": V, +# "t": t, +# "T_surf1": T_surf1, +# "T_surf2": T_surf2, +# "SOC": SOC, + +template = TemplateExpressionSpec( + expressions=["d", "e", "f", "g", "h"], + variable_names=["x1", "x2", "x3", "x4"], + parameters={"p1":1, "p2":1, "p3":1, "p4":1, "p5":1}, + combine="d(h(x1, x2) * log((f(x2))/(e(1 - x2)))) + g(x1)", +) + +template = TemplateExpressionSpec( + expressions=["g", "h"], + variable_names=["x1", "x2", "x3", "x4"], + parameters={"p1":1, "p2":1, "p3":1, "p4":1, "p5":1}, + combine="h(x1) * g(x2)", +) + +model = PySRRegressor( + maxsize=80, + niterations=1000, + binary_operators=["+", "-", "*", "/"], + unary_operators=["exp", "log"], + batching=True, + batch_size=200, +) + +output = model.fit(X, y) + +cos = lambda x: np.cos(x) +sin = lambda x: np.sin(x) +exp = lambda x: np.exp(x) +log = lambda x: np.log(x) + +def custom_fun_1(X): + x0 = X[:, 0] + x1 = X[:, 1] + x2 = X[:, 2] + x3 = X[:, 3] + return ((cos(x1 * 0.95293) * -0.59082) + cos(cos(sin(exp(x0 / 2.5028) + 0.46012)))) + 2.5555 + +def custom_fun_2(X): + x0 = X[:, 0] + x1 = X[:, 1] + x2 = X[:, 2] + x3 = X[:, 3] + return (cos(cos(exp(x0 + cos(x1 + x0)) + 0.74591)) + 2.5539) + (cos(x1 + -0.062883) * -0.57004) + +def custom_fun_3(X): + x0 = X[:, 0] + x1 = X[:, 1] + x2 = X[:, 2] + x3 = X[:, 3] + return (log(((exp(x0) + 19.943527) + (x1 * -3.1378086)) * (((x0 + (exp(exp((x2 / exp(x2 + -1.089517)) / exp(x3))) + ((x3 + -1.7805654) / exp(exp(exp(x0)))))) * 0.07967347) + 4.030623)) * 0.81259966) - -0.9257864 + +def custom_fun_4(X): + x0 = X[:, 0] + x1 = X[:, 1] + x2 = X[:, 2] + x3 = X[:, 3] + x4 = X[:, 4] + x5 = X[:, 5] + x6 = X[:, 6] + x7 = X[:, 7] + return ((x5 + ((((x2 / 1.3202298) + exp(x1)) / (exp(exp(x1) * 0.1189319) * 0.16781318)) + ((((x0 + x3) + exp((x6 + ((x1 - x3) * (exp(x6 - (((x2 - x1) + x3) * -0.033423472)) * 0.010805121))) + 4.2677464)) - exp(x1)) + ((((x7 + x5) + ((x3 * 0.40745232) + (x4 / (exp(x2) + 0.20448817)))) + -169.69832) * log(x1))))) * 0.003010543) + 4.5433464 + +for dcc in simulation.fu.DCC: + I = dcc.I + V = dcc.V + t = dcc.t + I_back1 = np.concatenate([np.zeros(1), I[:-1]]) + I_back2 = np.concatenate([np.zeros(3), I[:-3]]) + I_back3 = np.concatenate([np.zeros(7), I[:-7]]) + I_back4 = np.concatenate([np.zeros(20), I[:-20]]) + T_amb = dcc.T_amb + T_surf = dcc.T_surf + name = dcc.name + SOC = 3.0 - cumulative_simpson(y=I, x=t)/3600.0 + SOC = np.concatenate([np.array([3.0]), SOC]) + plt.plot(SOC, V, label="Measured Voltage") + plt.plot(SOC, custom_fun_4(np.column_stack([I, SOC, T_surf[0, :], T_surf[1, :], I_back1, I_back2, I_back3, I_back4])), label="Model 3") + plt.legend() + plt.title(f"DCC: {name} at T_amb={T_amb} C") + plt.xlabel("SOC") + plt.ylabel("Voltage [V]") + plt.grid(True) + plt.show() + +plt.plot(X[:, 1], custom_fun_1(X[:, 1], X[:, 2]), label="Model 1") +plt.plot(X[:, 1], y, label="Measured Voltage") +plt.xlabel("SOC [Ah]") +plt.ylabel("Voltage [V]") +plt.legend() +plt.show() + + +# ============================================================ +# VERSION 2: PARAMETERIZED HYSTERESIS KERNEL +# ============================================================ + +def buildHysteresisKernel(kernelLength, tauFast, tauSlow, mix=0.5, shape=1.0): + kernelLength = int(kernelLength) + if kernelLength <= 0: + raise ValueError("kernelLength must be positive") + + tauFast = max(float(tauFast), 1e-9) + tauSlow = max(float(tauSlow), 1e-9) + mix = float(np.clip(mix, 0.0, 1.0)) + shape = max(float(shape), 1e-9) + + lagIndex = np.arange(kernelLength, dtype=float) + fastComponent = np.exp(-((lagIndex / tauFast) ** shape)) + slowComponent = np.exp(-((lagIndex / tauSlow) ** shape)) + kernel = (mix * fastComponent) + ((1.0 - mix) * slowComponent) + + kernelSum = float(np.sum(kernel)) + if kernelSum == 0.0: + return np.ones(kernelLength, dtype=float) / float(kernelLength) + + return kernel / kernelSum + + +def buildHysteresisWindow(currentSeries, kernelLength): + currentSeries = np.asarray(currentSeries, dtype=float).reshape(-1) + kernelLength = int(kernelLength) + if kernelLength <= 0: + raise ValueError("kernelLength must be positive") + + paddedSeries = np.pad(currentSeries, (kernelLength - 1, 0), mode="constant") + windows = np.empty((currentSeries.shape[0], kernelLength), dtype=float) + + for sampleIndex in range(currentSeries.shape[0]): + window = paddedSeries[sampleIndex:sampleIndex + kernelLength] + windows[sampleIndex, :] = window[::-1] + + return windows + + +def buildHysteresisFeature(currentSeries, kernelLength, tauFast, tauSlow, mix=0.5, shape=1.0): + currentWindows = buildHysteresisWindow(currentSeries, kernelLength) + kernel = buildHysteresisKernel(kernelLength, tauFast, tauSlow, mix=mix, shape=shape) + return currentWindows @ kernel + + +def buildHysteresisTrainingSet(simulationData, kernelLength=64, tauFast=4.0, tauSlow=28.0, mix=0.7, shape=1.15): + featureRows = [] + targetRows = [] + + for dcc in simulationData.fu.DCC: + currentSeries = np.asarray(dcc.I, dtype=float).reshape(-1) + voltageSeries = np.asarray(dcc.V, dtype=float).reshape(-1) + timeSeries = np.asarray(dcc.t, dtype=float).reshape(-1) + surfaceSeries = np.asarray(dcc.T_surf, dtype=float) + + if surfaceSeries.ndim == 1: + surfaceOne = surfaceSeries + surfaceTwo = surfaceSeries + else: + surfaceOne = surfaceSeries[0, :] + surfaceTwo = surfaceSeries[1, :] + + socSeries = 3.0 - cumulative_simpson(y=currentSeries, x=timeSeries) / 3600.0 + socSeries = np.concatenate([np.array([3.0], dtype=float), socSeries]) + + hysteresisSeries = buildHysteresisFeature( + currentSeries=currentSeries, + kernelLength=kernelLength, + tauFast=tauFast, + tauSlow=tauSlow, + mix=mix, + shape=shape, + ) + + featureRows.append(np.column_stack([ + currentSeries, + socSeries, + surfaceOne, + surfaceTwo, + hysteresisSeries, + ])) + targetRows.append(voltageSeries.reshape(-1, 1)) + + features = np.vstack(featureRows) + targets = np.vstack(targetRows).ravel() + return features, targets + + +def trainHysteresisKernelVersion2(simulationData, kernelLength=128, tauFast=4.0, tauSlow=14.0, mix=0.7, shape=1.35): + features, targets = buildHysteresisTrainingSet( + simulationData=simulationData, + kernelLength=kernelLength, + tauFast=tauFast, + tauSlow=tauSlow, + mix=mix, + shape=shape, + ) + + modelV2 = PySRRegressor( + maxsize=30, + niterations=40, + binary_operators=["+", "-", "*", "/"], + unary_operators=["exp", "log"], + batching=True, + batch_size=200, + ) + + fittedModel = modelV2.fit(features, targets) + return fittedModel, features, targets + +model, features, targets = trainHysteresisKernelVersion2(simulation) + + +# Optional reference entry point for the new version: +# fittedModel, features, targets = trainHysteresisKernelVersion2(simulation) \ 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/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/outputs/20260605_170819_y7qEJI/checkpoint.pkl b/outputs/20260605_170819_y7qEJI/checkpoint.pkl new file mode 100644 index 0000000..d3f05a6 Binary files /dev/null and b/outputs/20260605_170819_y7qEJI/checkpoint.pkl differ diff --git a/outputs/20260605_170844_CfwX91/checkpoint.pkl b/outputs/20260605_170844_CfwX91/checkpoint.pkl new file mode 100644 index 0000000..a038793 Binary files /dev/null and b/outputs/20260605_170844_CfwX91/checkpoint.pkl differ diff --git a/outputs/20260605_170844_CfwX91/hall_of_fame.csv b/outputs/20260605_170844_CfwX91/hall_of_fame.csv new file mode 100644 index 0000000..ce586b2 --- /dev/null +++ b/outputs/20260605_170844_CfwX91/hall_of_fame.csv @@ -0,0 +1,9 @@ +Complexity,Loss,Equation +1,0.09237157,"3.578532" +3,0.015335105,"-15.249025 - x3" +5,0.003983757,"(x2 * -0.3919811) + 5.2681723" +7,0.00040154983,"(x2 / (x3 + 15.396925)) + 4.852261" +11,0.0003984477,"cos(x3) + (cos((x2 * -3.0206292) / x3) / 0.2936874)" +13,0.00035803785,"(1.7555461 + cos(x3)) + (cos(-4.532914 * (x2 / x3)) / 0.5895234)" +15,0.00034305936,"cos(x3) + ((cos(-0.22011857 + ((x2 * 3.0981772) / x3)) / 0.37743562) + 1.030288)" +18,0.0003428655,"((cos(x4) + (cos((x2 * (3.0963361 / x3)) + -0.21454689) / 0.37764242)) * 1.0070273) + cos(x3)" diff --git a/outputs/20260605_170844_CfwX91/hall_of_fame.csv.bak b/outputs/20260605_170844_CfwX91/hall_of_fame.csv.bak new file mode 100644 index 0000000..ce586b2 --- /dev/null +++ b/outputs/20260605_170844_CfwX91/hall_of_fame.csv.bak @@ -0,0 +1,9 @@ +Complexity,Loss,Equation +1,0.09237157,"3.578532" +3,0.015335105,"-15.249025 - x3" +5,0.003983757,"(x2 * -0.3919811) + 5.2681723" +7,0.00040154983,"(x2 / (x3 + 15.396925)) + 4.852261" +11,0.0003984477,"cos(x3) + (cos((x2 * -3.0206292) / x3) / 0.2936874)" +13,0.00035803785,"(1.7555461 + cos(x3)) + (cos(-4.532914 * (x2 / x3)) / 0.5895234)" +15,0.00034305936,"cos(x3) + ((cos(-0.22011857 + ((x2 * 3.0981772) / x3)) / 0.37743562) + 1.030288)" +18,0.0003428655,"((cos(x4) + (cos((x2 * (3.0963361 / x3)) + -0.21454689) / 0.37764242)) * 1.0070273) + cos(x3)" diff --git a/outputs/20260605_173335_rG1Rzg/checkpoint.pkl b/outputs/20260605_173335_rG1Rzg/checkpoint.pkl new file mode 100644 index 0000000..673d154 Binary files /dev/null and b/outputs/20260605_173335_rG1Rzg/checkpoint.pkl differ diff --git a/outputs/20260605_173335_rG1Rzg/hall_of_fame.csv b/outputs/20260605_173335_rG1Rzg/hall_of_fame.csv new file mode 100644 index 0000000..b8cd3db --- /dev/null +++ b/outputs/20260605_173335_rG1Rzg/hall_of_fame.csv @@ -0,0 +1,14 @@ +Complexity,Loss,Equation +1,0.14158258,"3.5650325" +4,0.046780355,"5.03204 - log(x1)" +5,0.020174779,"(x1 * -0.41577148) + 5.4018736" +6,0.01831764,"(cos(x1) * -0.5376725) + 3.4590318" +8,0.014986508,"((exp(x1) + -0.11164971) * -0.0039546792) + 4.0197625" +11,0.012626836,"(cos(cos(exp(x0))) + 2.7706904) + (cos(x1) * -0.5312622)" +12,0.012041195,"((cos(x1) * -0.53031766) + 2.5074618) + cos(exp(x0) + -0.67037964)" +13,0.01198075,"((cos(x1) * -0.5289256) + sin(cos(sin(x0 + 0.66593564)))) + 2.6993287" +14,0.011761692,"(cos(cos(sin(exp(x0 * 0.23291506)))) + 2.731931) + (cos(x1) * -0.5299537)" +15,0.011450616,"((cos(x1 + -0.16093712) * -0.56156677) + 2.440645) + cos(cos(exp(x0) + 0.86974084))" +16,0.010320034,"cos(cos(0.562462 + exp(cos(x1) + x0))) + (2.6217613 + (cos(x1) * -0.58870065))" +18,0.009740723,"(cos(x1) * -0.5753818) + (cos(cos(exp(cos(x1 + -0.232988) + x0) + 0.69828486)) + 2.5880654)" +20,0.009686999,"(cos(cos(exp(x0 + cos(x1 + x0)) + 0.7459149)) + 2.5538678) + (cos(x1 + -0.06288294) * -0.57003665)" diff --git a/outputs/20260605_173335_rG1Rzg/hall_of_fame.csv.bak b/outputs/20260605_173335_rG1Rzg/hall_of_fame.csv.bak new file mode 100644 index 0000000..b8cd3db --- /dev/null +++ b/outputs/20260605_173335_rG1Rzg/hall_of_fame.csv.bak @@ -0,0 +1,14 @@ +Complexity,Loss,Equation +1,0.14158258,"3.5650325" +4,0.046780355,"5.03204 - log(x1)" +5,0.020174779,"(x1 * -0.41577148) + 5.4018736" +6,0.01831764,"(cos(x1) * -0.5376725) + 3.4590318" +8,0.014986508,"((exp(x1) + -0.11164971) * -0.0039546792) + 4.0197625" +11,0.012626836,"(cos(cos(exp(x0))) + 2.7706904) + (cos(x1) * -0.5312622)" +12,0.012041195,"((cos(x1) * -0.53031766) + 2.5074618) + cos(exp(x0) + -0.67037964)" +13,0.01198075,"((cos(x1) * -0.5289256) + sin(cos(sin(x0 + 0.66593564)))) + 2.6993287" +14,0.011761692,"(cos(cos(sin(exp(x0 * 0.23291506)))) + 2.731931) + (cos(x1) * -0.5299537)" +15,0.011450616,"((cos(x1 + -0.16093712) * -0.56156677) + 2.440645) + cos(cos(exp(x0) + 0.86974084))" +16,0.010320034,"cos(cos(0.562462 + exp(cos(x1) + x0))) + (2.6217613 + (cos(x1) * -0.58870065))" +18,0.009740723,"(cos(x1) * -0.5753818) + (cos(cos(exp(cos(x1 + -0.232988) + x0) + 0.69828486)) + 2.5880654)" +20,0.009686999,"(cos(cos(exp(x0 + cos(x1 + x0)) + 0.7459149)) + 2.5538678) + (cos(x1 + -0.06288294) * -0.57003665)" diff --git a/outputs/20260605_220726_S696dn/checkpoint.pkl b/outputs/20260605_220726_S696dn/checkpoint.pkl new file mode 100644 index 0000000..fa05362 Binary files /dev/null and b/outputs/20260605_220726_S696dn/checkpoint.pkl differ diff --git a/outputs/20260605_220800_HYyXI8/checkpoint.pkl b/outputs/20260605_220800_HYyXI8/checkpoint.pkl new file mode 100644 index 0000000..d628d29 Binary files /dev/null and b/outputs/20260605_220800_HYyXI8/checkpoint.pkl differ diff --git a/outputs/20260605_220949_H8unUI/checkpoint.pkl b/outputs/20260605_220949_H8unUI/checkpoint.pkl new file mode 100644 index 0000000..9788288 Binary files /dev/null and b/outputs/20260605_220949_H8unUI/checkpoint.pkl differ diff --git a/outputs/20260605_221126_n50KhN/checkpoint.pkl b/outputs/20260605_221126_n50KhN/checkpoint.pkl new file mode 100644 index 0000000..767b32e Binary files /dev/null and b/outputs/20260605_221126_n50KhN/checkpoint.pkl differ diff --git a/outputs/20260606_001334_pEEQVZ/checkpoint.pkl b/outputs/20260606_001334_pEEQVZ/checkpoint.pkl new file mode 100644 index 0000000..4e02396 Binary files /dev/null and b/outputs/20260606_001334_pEEQVZ/checkpoint.pkl differ diff --git a/outputs/20260606_001458_3li5pJ/checkpoint.pkl b/outputs/20260606_001458_3li5pJ/checkpoint.pkl new file mode 100644 index 0000000..3dab8d7 Binary files /dev/null and b/outputs/20260606_001458_3li5pJ/checkpoint.pkl differ diff --git a/outputs/20260606_001533_vGUxTX/checkpoint.pkl b/outputs/20260606_001533_vGUxTX/checkpoint.pkl new file mode 100644 index 0000000..47d8abe Binary files /dev/null and b/outputs/20260606_001533_vGUxTX/checkpoint.pkl differ diff --git a/outputs/20260606_001540_oXXNe2/checkpoint.pkl b/outputs/20260606_001540_oXXNe2/checkpoint.pkl new file mode 100644 index 0000000..84816a3 Binary files /dev/null and b/outputs/20260606_001540_oXXNe2/checkpoint.pkl differ diff --git a/outputs/20260606_001550_Mx7y4w/checkpoint.pkl b/outputs/20260606_001550_Mx7y4w/checkpoint.pkl new file mode 100644 index 0000000..a2b9ec3 Binary files /dev/null and b/outputs/20260606_001550_Mx7y4w/checkpoint.pkl differ diff --git a/outputs/20260606_001847_IZ1eTW/checkpoint.pkl b/outputs/20260606_001847_IZ1eTW/checkpoint.pkl new file mode 100644 index 0000000..bcb57f4 Binary files /dev/null and b/outputs/20260606_001847_IZ1eTW/checkpoint.pkl differ diff --git a/outputs/20260606_001847_IZ1eTW/hall_of_fame.csv b/outputs/20260606_001847_IZ1eTW/hall_of_fame.csv new file mode 100644 index 0000000..95256c6 --- /dev/null +++ b/outputs/20260606_001847_IZ1eTW/hall_of_fame.csv @@ -0,0 +1,6 @@ +Complexity,Loss,Equation +1,0.14158258,"h = 3.5650237" +9,0.13296594,"h = exp(exp((#1 + -1.0300726) + #1)) - -2.4034257" +10,0.0766783,"h = ((#2 - -1.5896251) / log(#2 + -0.6743062)) + -1.0767896" +11,0.020174777,"h = (#2 + (((#2 + -1.4372089) + #2) * -0.70788854)) - -4.3845043" +17,0.014038164,"h = (#2 + (((#2 * ((#2 * -1.1257197) + #2)) + (#2 * -1.3065858)) + #2)) + 3.0436683" diff --git a/outputs/20260606_001847_IZ1eTW/hall_of_fame.csv.bak b/outputs/20260606_001847_IZ1eTW/hall_of_fame.csv.bak new file mode 100644 index 0000000..95256c6 --- /dev/null +++ b/outputs/20260606_001847_IZ1eTW/hall_of_fame.csv.bak @@ -0,0 +1,6 @@ +Complexity,Loss,Equation +1,0.14158258,"h = 3.5650237" +9,0.13296594,"h = exp(exp((#1 + -1.0300726) + #1)) - -2.4034257" +10,0.0766783,"h = ((#2 - -1.5896251) / log(#2 + -0.6743062)) + -1.0767896" +11,0.020174777,"h = (#2 + (((#2 + -1.4372089) + #2) * -0.70788854)) - -4.3845043" +17,0.014038164,"h = (#2 + (((#2 * ((#2 * -1.1257197) + #2)) + (#2 * -1.3065858)) + #2)) + 3.0436683" diff --git a/outputs/20260606_002128_XVS321/checkpoint.pkl b/outputs/20260606_002128_XVS321/checkpoint.pkl new file mode 100644 index 0000000..a25b590 Binary files /dev/null and b/outputs/20260606_002128_XVS321/checkpoint.pkl differ diff --git a/outputs/20260606_002230_5p7hPV/checkpoint.pkl b/outputs/20260606_002230_5p7hPV/checkpoint.pkl new file mode 100644 index 0000000..d946ab3 Binary files /dev/null and b/outputs/20260606_002230_5p7hPV/checkpoint.pkl differ diff --git a/outputs/20260606_003030_rjsoSk/checkpoint.pkl b/outputs/20260606_003030_rjsoSk/checkpoint.pkl new file mode 100644 index 0000000..28072f2 Binary files /dev/null and b/outputs/20260606_003030_rjsoSk/checkpoint.pkl differ diff --git a/outputs/20260606_003030_rjsoSk/hall_of_fame.csv b/outputs/20260606_003030_rjsoSk/hall_of_fame.csv new file mode 100644 index 0000000..117c5eb --- /dev/null +++ b/outputs/20260606_003030_rjsoSk/hall_of_fame.csv @@ -0,0 +1,12 @@ +Complexity,Loss,Equation +2,0.14158283,"g = 6.9857764; h = 0.5102515; p1 = [0.29919448]; p2 = [-0.0024157339]; p3 = [0.4125923]; p4 = [0.114938244]; p5 = [3.231683]" +4,0.020187454,"g = 12.910461 - #1; h = 0.41960278; p1 = [0.19131187]; p2 = [-0.31550545]; p3 = [0.2953269]; p4 = [1.2401006]; p5 = [-0.075753495]" +5,0.01674504,"g = log(9.668437 - #1); h = 2.1685836; p1 = [0.2311817]; p2 = [-0.4819484]; p3 = [0.27055913]; p4 = [1.0952356]; p5 = [-0.011553196]" +6,0.016334653,"g = log(log(11.231159 - #1)); h = 5.516899; p1 = [0.5851717]; p2 = [-0.31661797]; p3 = [0.27683184]; p4 = [3.3181517]; p5 = [-0.021765368]" +7,0.015008225,"g = (exp(#1) * 0.001560047) + -1.5977929; h = -2.511395; p1 = [-0.01028225]; p2 = [-0.20589283]; p3 = [0.20313159]; p4 = [1.3348466]; p5 = [0.06440573]" +9,0.011747481,"g = log(9.668437 - #1); h = 2.1685836 + (0.0063233213 * #1); p1 = [0.2311817]; p2 = [-0.4819484]; p3 = [0.27055913]; p4 = [1.0952356]; p5 = [0.01012363]" +12,0.011118251,"g = ((#1 * #1) * 0.02847272) + -2.690646; h = -1.6991075 - (#1 * 0.004837096); p1 = [-0.0062846392]; p2 = [-0.16279292]; p3 = [0.1604949]; p4 = [-0.787699]; p5 = [-0.020646032]" +13,0.009220802,"g = (exp(#1 / 1.8232644) * 0.027975967) + -1.9307917; h = (#1 * -0.007931266) + -2.2841105; p1 = [-0.040756416]; p2 = [-0.18525246]; p3 = [0.10550426]; p4 = [-0.6697017]; p5 = [-0.010399518]" +14,0.008201847,"g = (exp(#1 / 1.7637918) * 0.038912833) + -2.9290135; h = -1.4306883 + (-0.10357848 * exp(#1)); p1 = [-0.041650955]; p2 = [-0.21279228]; p3 = [-0.11671401]; p4 = [-0.90407455]; p5 = [0.0055799587]" +16,0.007718546,"g = (exp(#1 / 1.6343316) * 0.021368526) + -2.2231567; h = (exp(#1 / 2.5276432) * -0.12797768) + -1.8277506; p1 = [-0.040756416]; p2 = [-0.18525246]; p3 = [-0.11458365]; p4 = [-0.6697017]; p5 = [-0.010399518]" +18,0.0076752626,"g = (exp(#1 / 1.6343316) * 0.021368526) + -2.221973; h = (exp((#1 / 1.6343316) * 0.6894196) * -0.12797768) + -1.8277506; p1 = [-0.040756416]; p2 = [-0.18525246]; p3 = [0.10550426]; p4 = [-0.6697017]; p5 = [-0.010399518]" diff --git a/outputs/20260606_003030_rjsoSk/hall_of_fame.csv.bak b/outputs/20260606_003030_rjsoSk/hall_of_fame.csv.bak new file mode 100644 index 0000000..117c5eb --- /dev/null +++ b/outputs/20260606_003030_rjsoSk/hall_of_fame.csv.bak @@ -0,0 +1,12 @@ +Complexity,Loss,Equation +2,0.14158283,"g = 6.9857764; h = 0.5102515; p1 = [0.29919448]; p2 = [-0.0024157339]; p3 = [0.4125923]; p4 = [0.114938244]; p5 = [3.231683]" +4,0.020187454,"g = 12.910461 - #1; h = 0.41960278; p1 = [0.19131187]; p2 = [-0.31550545]; p3 = [0.2953269]; p4 = [1.2401006]; p5 = [-0.075753495]" +5,0.01674504,"g = log(9.668437 - #1); h = 2.1685836; p1 = [0.2311817]; p2 = [-0.4819484]; p3 = [0.27055913]; p4 = [1.0952356]; p5 = [-0.011553196]" +6,0.016334653,"g = log(log(11.231159 - #1)); h = 5.516899; p1 = [0.5851717]; p2 = [-0.31661797]; p3 = [0.27683184]; p4 = [3.3181517]; p5 = [-0.021765368]" +7,0.015008225,"g = (exp(#1) * 0.001560047) + -1.5977929; h = -2.511395; p1 = [-0.01028225]; p2 = [-0.20589283]; p3 = [0.20313159]; p4 = [1.3348466]; p5 = [0.06440573]" +9,0.011747481,"g = log(9.668437 - #1); h = 2.1685836 + (0.0063233213 * #1); p1 = [0.2311817]; p2 = [-0.4819484]; p3 = [0.27055913]; p4 = [1.0952356]; p5 = [0.01012363]" +12,0.011118251,"g = ((#1 * #1) * 0.02847272) + -2.690646; h = -1.6991075 - (#1 * 0.004837096); p1 = [-0.0062846392]; p2 = [-0.16279292]; p3 = [0.1604949]; p4 = [-0.787699]; p5 = [-0.020646032]" +13,0.009220802,"g = (exp(#1 / 1.8232644) * 0.027975967) + -1.9307917; h = (#1 * -0.007931266) + -2.2841105; p1 = [-0.040756416]; p2 = [-0.18525246]; p3 = [0.10550426]; p4 = [-0.6697017]; p5 = [-0.010399518]" +14,0.008201847,"g = (exp(#1 / 1.7637918) * 0.038912833) + -2.9290135; h = -1.4306883 + (-0.10357848 * exp(#1)); p1 = [-0.041650955]; p2 = [-0.21279228]; p3 = [-0.11671401]; p4 = [-0.90407455]; p5 = [0.0055799587]" +16,0.007718546,"g = (exp(#1 / 1.6343316) * 0.021368526) + -2.2231567; h = (exp(#1 / 2.5276432) * -0.12797768) + -1.8277506; p1 = [-0.040756416]; p2 = [-0.18525246]; p3 = [-0.11458365]; p4 = [-0.6697017]; p5 = [-0.010399518]" +18,0.0076752626,"g = (exp(#1 / 1.6343316) * 0.021368526) + -2.221973; h = (exp((#1 / 1.6343316) * 0.6894196) * -0.12797768) + -1.8277506; p1 = [-0.040756416]; p2 = [-0.18525246]; p3 = [0.10550426]; p4 = [-0.6697017]; p5 = [-0.010399518]" diff --git a/outputs/20260606_003355_3mqPft/checkpoint.pkl b/outputs/20260606_003355_3mqPft/checkpoint.pkl new file mode 100644 index 0000000..e222eca Binary files /dev/null and b/outputs/20260606_003355_3mqPft/checkpoint.pkl differ diff --git a/outputs/20260606_003355_3mqPft/hall_of_fame.csv b/outputs/20260606_003355_3mqPft/hall_of_fame.csv new file mode 100644 index 0000000..33d56d5 --- /dev/null +++ b/outputs/20260606_003355_3mqPft/hall_of_fame.csv @@ -0,0 +1,8 @@ +Complexity,Loss,Equation +5,0.1415826,"d = 0.8262559; e = 2.4370785; f = 1.2696979; g = 2.7386935; h = 1.0896828; p1 = [-0.00013334728]; p2 = [-0.037495513]; p3 = [-1.860748]; p4 = [0.1828995]; p5 = [-0.009285139]" +6,0.014987052,"d = #1; e = 1.750163; f = 1.7432467; g = 4.02136; h = exp(#2); p1 = [-0.00018178526]; p2 = [-0.10732147]; p3 = [1.8080505]; p4 = [-0.24146089]; p5 = [-0.014681858]" +8,0.012615347,"d = #1; e = 4.2365522; f = 4.21979; g = 4.028829; h = exp(#2) - #1; p1 = [-0.00024798003]; p2 = [-0.118860826]; p3 = [4.0861635]; p4 = [-0.28409734]; p5 = [0.07718387]" +10,0.0113293165,"d = #1; e = 2.2003496; f = 2.1917112; g = 4.0349083; h = (exp(#2) - #1) - #1; p1 = [-0.00018178526]; p2 = [-0.10732147]; p3 = [1.8080505]; p4 = [-0.24146089]; p5 = [-0.014681858]" +11,0.007727116,"d = log(#1); e = -0.16336742 - #1; f = 5.550159; g = 2.6645992; h = exp(#1) + #2; p1 = [-2.2537508e-5]; p2 = [-0.13228661]; p3 = [0.4082183]; p4 = [-0.11639516]; p5 = [-0.0066752583]" +12,0.007709185,"d = log(#1); e = 0.77743286 - #1; f = 6.523924; g = 2.6252806; h = #2 + exp(exp(#1)); p1 = [-1.6830088e-5]; p2 = [-0.17052254]; p3 = [0.44067693]; p4 = [-0.21197277]; p5 = [-0.010168037]" +13,0.0073829982,"d = log(#1); e = 0.040130213 - #1; f = 5.7612133; g = 2.6771636; h = #2 + exp(#1 * 0.4180934); p1 = [-5.9369708e-5]; p2 = [-0.61107564]; p3 = [-1.3713654]; p4 = [-0.1646734]; p5 = [-0.03565868]" diff --git a/outputs/20260606_003355_3mqPft/hall_of_fame.csv.bak b/outputs/20260606_003355_3mqPft/hall_of_fame.csv.bak new file mode 100644 index 0000000..33d56d5 --- /dev/null +++ b/outputs/20260606_003355_3mqPft/hall_of_fame.csv.bak @@ -0,0 +1,8 @@ +Complexity,Loss,Equation +5,0.1415826,"d = 0.8262559; e = 2.4370785; f = 1.2696979; g = 2.7386935; h = 1.0896828; p1 = [-0.00013334728]; p2 = [-0.037495513]; p3 = [-1.860748]; p4 = [0.1828995]; p5 = [-0.009285139]" +6,0.014987052,"d = #1; e = 1.750163; f = 1.7432467; g = 4.02136; h = exp(#2); p1 = [-0.00018178526]; p2 = [-0.10732147]; p3 = [1.8080505]; p4 = [-0.24146089]; p5 = [-0.014681858]" +8,0.012615347,"d = #1; e = 4.2365522; f = 4.21979; g = 4.028829; h = exp(#2) - #1; p1 = [-0.00024798003]; p2 = [-0.118860826]; p3 = [4.0861635]; p4 = [-0.28409734]; p5 = [0.07718387]" +10,0.0113293165,"d = #1; e = 2.2003496; f = 2.1917112; g = 4.0349083; h = (exp(#2) - #1) - #1; p1 = [-0.00018178526]; p2 = [-0.10732147]; p3 = [1.8080505]; p4 = [-0.24146089]; p5 = [-0.014681858]" +11,0.007727116,"d = log(#1); e = -0.16336742 - #1; f = 5.550159; g = 2.6645992; h = exp(#1) + #2; p1 = [-2.2537508e-5]; p2 = [-0.13228661]; p3 = [0.4082183]; p4 = [-0.11639516]; p5 = [-0.0066752583]" +12,0.007709185,"d = log(#1); e = 0.77743286 - #1; f = 6.523924; g = 2.6252806; h = #2 + exp(exp(#1)); p1 = [-1.6830088e-5]; p2 = [-0.17052254]; p3 = [0.44067693]; p4 = [-0.21197277]; p5 = [-0.010168037]" +13,0.0073829982,"d = log(#1); e = 0.040130213 - #1; f = 5.7612133; g = 2.6771636; h = #2 + exp(#1 * 0.4180934); p1 = [-5.9369708e-5]; p2 = [-0.61107564]; p3 = [-1.3713654]; p4 = [-0.1646734]; p5 = [-0.03565868]" diff --git a/outputs/20260606_003733_RjGvK0/checkpoint.pkl b/outputs/20260606_003733_RjGvK0/checkpoint.pkl new file mode 100644 index 0000000..cb0387e Binary files /dev/null and b/outputs/20260606_003733_RjGvK0/checkpoint.pkl differ diff --git a/outputs/20260606_003733_RjGvK0/hall_of_fame.csv b/outputs/20260606_003733_RjGvK0/hall_of_fame.csv new file mode 100644 index 0000000..7b40c06 --- /dev/null +++ b/outputs/20260606_003733_RjGvK0/hall_of_fame.csv @@ -0,0 +1,24 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5650616" +4,0.046783242,"5.0337377 - log(x1)" +5,0.020176543,"5.3953233 - (x1 * 0.41445813)" +6,0.014986969,"4.020731 - (exp(x1) * 0.003961596)" +8,0.012617102,"((exp(x1) - x0) * -0.003961628) + 4.0327845" +10,0.01116037,"((x0 * 0.009852909) - (exp(x1) * 0.0038163126)) + 4.0310163" +11,0.010095499,"((exp(x0 + 4.0630603) - exp(x1)) * 0.0039469204) + 3.8907137" +12,0.007796827,"((x2 - (exp(x1) - (x0 / 0.28251338))) * 0.004027564) + 4.016739" +14,0.006879552,"(((x2 * 0.72872156) - (exp(x1) - (x0 / 0.26986617))) * 0.0040161763) + 4.0202312" +16,0.006878811,"((((x2 + -0.7189338) * 0.72872156) - (exp(x1) - (x0 / 0.26986617))) * 0.0040161763) + 4.0202312" +17,0.0065005305,"((x2 - (exp(x1) - (x0 - ((exp(x0) - x1) * -52.59801)))) * 0.0026673693) + 4.40384" +18,0.0050513665,"((x3 - (exp(x1) - ((x1 - (x0 - (x1 * x1))) / -0.15749042))) * 0.0025426166) + 4.253663" +19,0.0038443354,"((x2 - (exp(x1) - ((x0 - ((exp(x0) - x1) * -10.510046)) / 0.19981831))) * 0.0026673693) + 4.40384" +21,0.0038443035,"((x2 - ((-0.0588477 + exp(x1)) - ((x0 - ((exp(x0) - x1) * -10.510046)) / 0.19981831))) * 0.0026673693) + 4.40384" +25,0.0037783843,"(((x3 - (x0 / -0.073459424)) - ((exp(x1) - (x2 - ((exp(x0) - x1) * (x1 / -0.08350564)))) / 0.41662335)) * 0.00087618356) + 4.2367797" +27,0.0037783666,"(((x3 - (x0 / -0.073459424)) - ((exp(x1) - ((x2 + 0.19008082) - ((exp(x0) - x1) * (x1 / -0.08350564)))) / 0.41662335)) * 0.00087618356) + 4.2367797" +30,0.003463204,"(0.0020726437 * (((((x0 / (-0.87535447 - exp(x3))) - x0) / -0.20812978) + x1) - (exp(x1) - (x3 - ((x1 * (exp(x0) - x1)) / -0.079456925))))) + 4.251038" +32,0.00339266,"((((((x1 - (x0 / (exp(x2) - -0.86222446))) - x0) / -0.17791161) + -0.52631164) - (exp(x1) - (x3 - ((exp(x0) - x1) * (x1 / -0.11411076))))) * 0.0023796074) + 4.268963" +33,0.0033830802,"(0.0020726437 * ((((((x0 - log(1.773832)) / (-0.87535447 - exp(x3))) - x0) / -0.20812978) + x1) - (exp(x1) - (x3 - ((x1 * (exp(x0) - x1)) / -0.079456925))))) + 4.251038" +34,0.0029980664,"(((((x1 + 10.736047) * (0.58838904 - (x0 / (-1.0434699 - exp(x2))))) - (x0 / -0.21245426)) - (exp(x1) - (x2 - (((exp(x0) - x1) * x1) / -0.085821964)))) * 0.0021891403) + 4.2590456" +36,0.00295498,"(((((((0.8434696 - x0) / (-1.0434699 - exp(x2))) * 10.736047) - (x0 / -0.21245426)) + (x1 - (exp(x1) - (x2 - (x1 * ((exp(x0) - x1) / -0.085821964)))))) * 0.0021891403) + 3.842706) - -0.4163396" +37,0.002899302,"((((((0.58838904 - x0) / (-1.0434699 - exp(x2))) * exp(2.5817208)) - (x0 / -0.21245426)) + (x1 - (exp(x1) - (x2 - ((x1 / -0.085821964) * (exp(x0) - x1)))))) * 0.0021891403) + (3.842706 - -0.4163396)" +38,0.002821913,"((((((0.91504276 - (x0 / (-0.95385784 - (exp(x2) * x1)))) / -0.18839099) - x0) / -0.33124775) + ((x3 - (exp(x1) - (x0 - ((exp(x0) - x1) * (x1 / -0.09684288))))) + x0)) * 0.002278467) + 4.214192" diff --git a/outputs/20260606_003733_RjGvK0/hall_of_fame.csv.bak b/outputs/20260606_003733_RjGvK0/hall_of_fame.csv.bak new file mode 100644 index 0000000..7b40c06 --- /dev/null +++ b/outputs/20260606_003733_RjGvK0/hall_of_fame.csv.bak @@ -0,0 +1,24 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5650616" +4,0.046783242,"5.0337377 - log(x1)" +5,0.020176543,"5.3953233 - (x1 * 0.41445813)" +6,0.014986969,"4.020731 - (exp(x1) * 0.003961596)" +8,0.012617102,"((exp(x1) - x0) * -0.003961628) + 4.0327845" +10,0.01116037,"((x0 * 0.009852909) - (exp(x1) * 0.0038163126)) + 4.0310163" +11,0.010095499,"((exp(x0 + 4.0630603) - exp(x1)) * 0.0039469204) + 3.8907137" +12,0.007796827,"((x2 - (exp(x1) - (x0 / 0.28251338))) * 0.004027564) + 4.016739" +14,0.006879552,"(((x2 * 0.72872156) - (exp(x1) - (x0 / 0.26986617))) * 0.0040161763) + 4.0202312" +16,0.006878811,"((((x2 + -0.7189338) * 0.72872156) - (exp(x1) - (x0 / 0.26986617))) * 0.0040161763) + 4.0202312" +17,0.0065005305,"((x2 - (exp(x1) - (x0 - ((exp(x0) - x1) * -52.59801)))) * 0.0026673693) + 4.40384" +18,0.0050513665,"((x3 - (exp(x1) - ((x1 - (x0 - (x1 * x1))) / -0.15749042))) * 0.0025426166) + 4.253663" +19,0.0038443354,"((x2 - (exp(x1) - ((x0 - ((exp(x0) - x1) * -10.510046)) / 0.19981831))) * 0.0026673693) + 4.40384" +21,0.0038443035,"((x2 - ((-0.0588477 + exp(x1)) - ((x0 - ((exp(x0) - x1) * -10.510046)) / 0.19981831))) * 0.0026673693) + 4.40384" +25,0.0037783843,"(((x3 - (x0 / -0.073459424)) - ((exp(x1) - (x2 - ((exp(x0) - x1) * (x1 / -0.08350564)))) / 0.41662335)) * 0.00087618356) + 4.2367797" +27,0.0037783666,"(((x3 - (x0 / -0.073459424)) - ((exp(x1) - ((x2 + 0.19008082) - ((exp(x0) - x1) * (x1 / -0.08350564)))) / 0.41662335)) * 0.00087618356) + 4.2367797" +30,0.003463204,"(0.0020726437 * (((((x0 / (-0.87535447 - exp(x3))) - x0) / -0.20812978) + x1) - (exp(x1) - (x3 - ((x1 * (exp(x0) - x1)) / -0.079456925))))) + 4.251038" +32,0.00339266,"((((((x1 - (x0 / (exp(x2) - -0.86222446))) - x0) / -0.17791161) + -0.52631164) - (exp(x1) - (x3 - ((exp(x0) - x1) * (x1 / -0.11411076))))) * 0.0023796074) + 4.268963" +33,0.0033830802,"(0.0020726437 * ((((((x0 - log(1.773832)) / (-0.87535447 - exp(x3))) - x0) / -0.20812978) + x1) - (exp(x1) - (x3 - ((x1 * (exp(x0) - x1)) / -0.079456925))))) + 4.251038" +34,0.0029980664,"(((((x1 + 10.736047) * (0.58838904 - (x0 / (-1.0434699 - exp(x2))))) - (x0 / -0.21245426)) - (exp(x1) - (x2 - (((exp(x0) - x1) * x1) / -0.085821964)))) * 0.0021891403) + 4.2590456" +36,0.00295498,"(((((((0.8434696 - x0) / (-1.0434699 - exp(x2))) * 10.736047) - (x0 / -0.21245426)) + (x1 - (exp(x1) - (x2 - (x1 * ((exp(x0) - x1) / -0.085821964)))))) * 0.0021891403) + 3.842706) - -0.4163396" +37,0.002899302,"((((((0.58838904 - x0) / (-1.0434699 - exp(x2))) * exp(2.5817208)) - (x0 / -0.21245426)) + (x1 - (exp(x1) - (x2 - ((x1 / -0.085821964) * (exp(x0) - x1)))))) * 0.0021891403) + (3.842706 - -0.4163396)" +38,0.002821913,"((((((0.91504276 - (x0 / (-0.95385784 - (exp(x2) * x1)))) / -0.18839099) - x0) / -0.33124775) + ((x3 - (exp(x1) - (x0 - ((exp(x0) - x1) * (x1 / -0.09684288))))) + x0)) * 0.002278467) + 4.214192" diff --git a/outputs/20260606_003854_LTXdBL/checkpoint.pkl b/outputs/20260606_003854_LTXdBL/checkpoint.pkl new file mode 100644 index 0000000..563b723 Binary files /dev/null and b/outputs/20260606_003854_LTXdBL/checkpoint.pkl differ diff --git a/outputs/20260606_003854_LTXdBL/hall_of_fame.csv b/outputs/20260606_003854_LTXdBL/hall_of_fame.csv new file mode 100644 index 0000000..29e9c09 --- /dev/null +++ b/outputs/20260606_003854_LTXdBL/hall_of_fame.csv @@ -0,0 +1,28 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5650191" +4,0.04678074,"log(153.34027 / x1)" +5,0.02017487,"(x1 * -0.41594672) + 5.402903" +6,0.013971531,"log((x1 * -14.93622) + 103.58042)" +8,0.01352742,"4.245308 - exp((x1 * 0.64879787) + -3.3977587)" +9,0.010056367,"log((x1 * (exp(x0) + -15.155802)) + 102.0103)" +10,0.009054727,"log((x1 * -15.020367) + 104.825935) + (x0 * 0.011758094)" +11,0.008198477,"log((exp(x0) + 4.044639) * ((x1 * -3.2063076) + 22.328203))" +13,0.00782357,"log((exp(x0 / 2.2473528) + 3.935948) * ((x1 * -3.208584) + 22.315092))" +14,0.005083718,"log((x1 * -15.061507) + 103.96373) + ((x0 + (x2 / 5.6027427)) * 0.01664235)" +16,0.0050626895,"((((x1 + x2) / 5.526403) + x0) * 0.01654697) + log((x1 * -14.824304) + 102.43587)" +17,0.004654618,"log(((x1 * -14.915128) + 102.57765) + exp(x0)) + ((x0 + (x2 / 4.9587984)) * 0.014196217)" +18,0.0044531925,"((x0 + (x2 / 5.055242)) * 0.01466259) + log((exp(exp(x0)) + (x1 * -15.143642)) + 102.50558)" +19,0.0039043254,"log((exp(x0) + ((x1 * -3.4875638) + 23.611822)) * (((x0 + (x3 / 3.7873995)) * 0.048267443) + 4.2391734))" +20,0.0038930338,"log(((((x3 / 3.9647164) + x0) * 0.029658552) + 2.5351267) * ((exp(exp(x0)) + (x1 * -5.8674717)) + 38.74419))" +22,0.0030119894,"log((exp(x0) + ((x1 * -3.17398) + 21.280706)) * (((x0 + (x2 / exp(exp(exp(x0))))) * 0.08173577) + 4.7711587))" +24,0.0029371195,"log((((x1 + (x1 * -3.8719556)) + 19.228865) + exp(x0)) * ((((x2 / exp(exp(exp(x0)))) + x0) * 0.08566952) + 5.2742534))" +26,0.002736969,"(log((((x0 + (x2 / exp(exp(exp(x0))))) * 0.077206746) + 4.027057) * (exp(x0) + ((x1 * -3.1296072) + 19.944649))) * 0.81824416) - -0.93926245" +28,0.0026637525,"(log((exp(x0) + (19.954786 + (x1 * -3.1241653))) * ((((x3 / exp(exp(exp(x0 / 0.65302724)))) + x0) * 0.0752268) + 4.0193954)) * 0.8282046) - -0.89653075" +29,0.0026496013,"(log(((((x3 / exp(exp(exp(x0 / exp(x0))))) + x0) * 0.0752268) + 4.0193954) * (exp(x0) + ((x1 * -3.1241653) + 19.954786))) * 0.8282046) - -0.89653075" +31,0.0026125403,"((x1 + (x0 + ((x3 + -12.082671) / ((-0.33471185 / (x0 + 0.24264842)) + 2.724668)))) * 0.017363971) + log((exp(x0) + (x1 + (25.853796 + (x1 * -4.8740983)))) * 3.8618696)" +32,0.002452399,"(log((exp(x0) + ((x1 * -3.1378086) + 19.943527)) * ((((x0 + (x2 / exp(exp(exp(x0))))) + exp(x3 / exp(x2))) * 0.07967347) + 4.030623)) * 0.81259966) - -0.9503631" +33,0.0022979595,"(log((exp(x0) + ((x1 * -3.1378086) + 19.943527)) * ((((x0 + exp(exp(x3 / exp(x3)))) + (x2 / exp(exp(exp(x0))))) * 0.07967347) + 4.030623)) * 0.81259966) - -0.9257864" +35,0.0022652051,"(log((((x1 * -3.1480877) + 19.94338) + exp(x0)) * (((x0 + (exp(exp(x3 / exp(x3))) + ((x3 - 2.4743218) / exp(exp(exp(x0)))))) * 0.078905515) + 4.0369687)) * 0.80609787) - -0.9546834" +37,0.002233121,"(log((((x1 * -3.1378086) + 19.943527) + exp(x0)) * ((((x0 + exp(exp((x2 + 0.33637893) / exp(x3)))) + ((x2 + -1.7805654) / exp(exp(exp(x0))))) * 0.07967347) + 4.030623)) * 0.81259966) - -0.9257864" +38,0.0022321197,"(log((exp(x0) + ((-3.144868 * x1) + 19.9496)) * (((exp(exp((x2 / exp(x3 + -1.0905513)) / exp(x2))) + (x0 + (x3 / exp(exp(exp(x0)))))) * 0.08232028) + 4.032034)) * 0.8099617) - -0.9318484" +40,0.0022055435,"(log(((exp(x0) + 19.943527) + (x1 * -3.1378086)) * (((x0 + (exp(exp((x2 / exp(x2 + -1.089517)) / exp(x3))) + ((x3 + -1.7805654) / exp(exp(exp(x0)))))) * 0.07967347) + 4.030623)) * 0.81259966) - -0.9257864" diff --git a/outputs/20260606_003854_LTXdBL/hall_of_fame.csv.bak b/outputs/20260606_003854_LTXdBL/hall_of_fame.csv.bak new file mode 100644 index 0000000..29e9c09 --- /dev/null +++ b/outputs/20260606_003854_LTXdBL/hall_of_fame.csv.bak @@ -0,0 +1,28 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5650191" +4,0.04678074,"log(153.34027 / x1)" +5,0.02017487,"(x1 * -0.41594672) + 5.402903" +6,0.013971531,"log((x1 * -14.93622) + 103.58042)" +8,0.01352742,"4.245308 - exp((x1 * 0.64879787) + -3.3977587)" +9,0.010056367,"log((x1 * (exp(x0) + -15.155802)) + 102.0103)" +10,0.009054727,"log((x1 * -15.020367) + 104.825935) + (x0 * 0.011758094)" +11,0.008198477,"log((exp(x0) + 4.044639) * ((x1 * -3.2063076) + 22.328203))" +13,0.00782357,"log((exp(x0 / 2.2473528) + 3.935948) * ((x1 * -3.208584) + 22.315092))" +14,0.005083718,"log((x1 * -15.061507) + 103.96373) + ((x0 + (x2 / 5.6027427)) * 0.01664235)" +16,0.0050626895,"((((x1 + x2) / 5.526403) + x0) * 0.01654697) + log((x1 * -14.824304) + 102.43587)" +17,0.004654618,"log(((x1 * -14.915128) + 102.57765) + exp(x0)) + ((x0 + (x2 / 4.9587984)) * 0.014196217)" +18,0.0044531925,"((x0 + (x2 / 5.055242)) * 0.01466259) + log((exp(exp(x0)) + (x1 * -15.143642)) + 102.50558)" +19,0.0039043254,"log((exp(x0) + ((x1 * -3.4875638) + 23.611822)) * (((x0 + (x3 / 3.7873995)) * 0.048267443) + 4.2391734))" +20,0.0038930338,"log(((((x3 / 3.9647164) + x0) * 0.029658552) + 2.5351267) * ((exp(exp(x0)) + (x1 * -5.8674717)) + 38.74419))" +22,0.0030119894,"log((exp(x0) + ((x1 * -3.17398) + 21.280706)) * (((x0 + (x2 / exp(exp(exp(x0))))) * 0.08173577) + 4.7711587))" +24,0.0029371195,"log((((x1 + (x1 * -3.8719556)) + 19.228865) + exp(x0)) * ((((x2 / exp(exp(exp(x0)))) + x0) * 0.08566952) + 5.2742534))" +26,0.002736969,"(log((((x0 + (x2 / exp(exp(exp(x0))))) * 0.077206746) + 4.027057) * (exp(x0) + ((x1 * -3.1296072) + 19.944649))) * 0.81824416) - -0.93926245" +28,0.0026637525,"(log((exp(x0) + (19.954786 + (x1 * -3.1241653))) * ((((x3 / exp(exp(exp(x0 / 0.65302724)))) + x0) * 0.0752268) + 4.0193954)) * 0.8282046) - -0.89653075" +29,0.0026496013,"(log(((((x3 / exp(exp(exp(x0 / exp(x0))))) + x0) * 0.0752268) + 4.0193954) * (exp(x0) + ((x1 * -3.1241653) + 19.954786))) * 0.8282046) - -0.89653075" +31,0.0026125403,"((x1 + (x0 + ((x3 + -12.082671) / ((-0.33471185 / (x0 + 0.24264842)) + 2.724668)))) * 0.017363971) + log((exp(x0) + (x1 + (25.853796 + (x1 * -4.8740983)))) * 3.8618696)" +32,0.002452399,"(log((exp(x0) + ((x1 * -3.1378086) + 19.943527)) * ((((x0 + (x2 / exp(exp(exp(x0))))) + exp(x3 / exp(x2))) * 0.07967347) + 4.030623)) * 0.81259966) - -0.9503631" +33,0.0022979595,"(log((exp(x0) + ((x1 * -3.1378086) + 19.943527)) * ((((x0 + exp(exp(x3 / exp(x3)))) + (x2 / exp(exp(exp(x0))))) * 0.07967347) + 4.030623)) * 0.81259966) - -0.9257864" +35,0.0022652051,"(log((((x1 * -3.1480877) + 19.94338) + exp(x0)) * (((x0 + (exp(exp(x3 / exp(x3))) + ((x3 - 2.4743218) / exp(exp(exp(x0)))))) * 0.078905515) + 4.0369687)) * 0.80609787) - -0.9546834" +37,0.002233121,"(log((((x1 * -3.1378086) + 19.943527) + exp(x0)) * ((((x0 + exp(exp((x2 + 0.33637893) / exp(x3)))) + ((x2 + -1.7805654) / exp(exp(exp(x0))))) * 0.07967347) + 4.030623)) * 0.81259966) - -0.9257864" +38,0.0022321197,"(log((exp(x0) + ((-3.144868 * x1) + 19.9496)) * (((exp(exp((x2 / exp(x3 + -1.0905513)) / exp(x2))) + (x0 + (x3 / exp(exp(exp(x0)))))) * 0.08232028) + 4.032034)) * 0.8099617) - -0.9318484" +40,0.0022055435,"(log(((exp(x0) + 19.943527) + (x1 * -3.1378086)) * (((x0 + (exp(exp((x2 / exp(x2 + -1.089517)) / exp(x3))) + ((x3 + -1.7805654) / exp(exp(exp(x0)))))) * 0.07967347) + 4.030623)) * 0.81259966) - -0.9257864" diff --git a/outputs/20260606_005734_6FlOfR/checkpoint.pkl b/outputs/20260606_005734_6FlOfR/checkpoint.pkl new file mode 100644 index 0000000..20e44a3 Binary files /dev/null and b/outputs/20260606_005734_6FlOfR/checkpoint.pkl differ diff --git a/outputs/20260606_005734_6FlOfR/hall_of_fame.csv b/outputs/20260606_005734_6FlOfR/hall_of_fame.csv new file mode 100644 index 0000000..e1b00ba --- /dev/null +++ b/outputs/20260606_005734_6FlOfR/hall_of_fame.csv @@ -0,0 +1,44 @@ +Complexity,Loss,Equation +1,0.1415826,"3.565013" +4,0.046780355,"log(153.24149 / x1)" +5,0.02017629,"(x1 * -0.41636097) + 5.4033523" +6,0.0149874715,"(exp(x1) * -0.0039443914) + 4.018696" +7,0.014847774,"(-10.5967 / x1) + (10.474295 - x1)" +8,0.012612821,"((x0 - exp(x1)) * 0.003974942) + 4.032723" +10,0.011021733,"(x7 * 0.011055068) + (4.0363736 - (exp(x1) * 0.0038689962))" +11,0.01002806,"(3.865385 - (exp(x1) * 0.003869981)) + (exp(x7) * 0.2532472)" +12,0.007341251,"((x7 * 0.016047547) + 4.005379) + ((x2 - exp(x1)) * 0.00402841)" +13,0.0069549805,"exp(x7 * 0.022573786) + (((x2 - exp(x1)) * 0.003997351) + 3.0122573)" +14,0.0053718733,"((((x1 * (x0 + -37.00079)) + x2) - exp(x1)) * 0.0030013788) + 4.389338" +16,0.004099731,"(((x2 - exp(x1)) * 0.0031908343) + 3.430174) + (2.2044883 / (x1 + (x7 * -0.49608624)))" +18,0.0039598015,"(((x2 * 0.8246542) - exp(x1)) * 0.0031952073) + ((2.2756846 / (x1 + (x7 * -0.4195139))) + 3.4163985)" +19,0.003771442,"(((exp(x7 + 4.0886) + (x1 * (x0 + -64.57939))) + (x3 - exp(x1))) * 0.002504624) + 4.4716077" +21,0.0034251248,"4.4838834 + (0.0025647841 * ((x2 + exp(x7 + 4.1382027)) + (((x1 * -64.71305) - exp(x1)) + (x4 * 3.854656))))" +22,0.0032816748,"((((log(x1) * -229.28015) + ((x0 / 0.26969364) + (x3 - exp(x1)))) + exp(x7 + 4.1170673)) * 0.0027460016) + 4.6925445" +24,0.0032803074,"((((x0 / 0.26969364) + exp(x7 + 4.1170673)) + ((-0.23561807 + (x3 - exp(x1))) + (log(x1) * -229.28015))) * 0.0027460016) + 4.6925445" +26,0.0031384756,"((((exp((x7 + 4.129333) - (x2 * 0.0056593237)) - exp(x1)) + (x2 + (log(x1) * -169.4319))) + (x0 * 3.3012187)) * 0.003012371) + 4.537966" +28,0.0026985996,"((exp(((x3 * x3) * -0.00079846877) + (x7 - -4.3820744)) + (((((x5 * 2.9971557) + x3) + -171.94554) * log(x1)) - exp(x1))) * 0.0030115678) + 4.538426" +30,0.002533871,"((exp((x3 * ((x3 + 9.939942) * -0.0006410332)) + (x7 - -4.3820744)) + (((((x0 * 3.0800557) + x3) + -171.76453) * log(x1)) - exp(x1))) * 0.0030115678) + 4.538426" +32,0.0025038442,"(((exp((((-5.920857 - x3) * x3) * 0.00073056546) + (x0 + 4.430383)) + ((((x7 + (x3 + -169.90005)) + x4) + x7) * log(x1))) - exp(x1)) * 0.003055337) + 4.5436344" +33,0.0021212408,"(((log(x1) * (((x5 + x7) / 0.58209914) + -166.40523)) + (x2 + ((x3 - exp(x1)) + exp((4.6342545 - (exp(x2 * 0.09402527) * 0.082178235)) + x0)))) * 0.0030390304) + 4.498984" +35,0.002008648,"((((x3 - exp(x1)) + ((x7 + (((x6 + -169.14818) + x0) * log(x1))) * 1.3442688)) + (exp((x0 + 4.5737853) - (exp(x2 * 0.095374905) * 0.07284515)) + x3)) * 0.002771325) + 4.6696644" +37,0.0019595241,"((((exp(x0 + (4.4453154 - (exp(x0 + (x2 * 0.07406609)) * 0.17924055))) + ((((x0 + x7) + -168.80672) * log(x1)) + x7)) * 1.2684214) + ((x2 - exp(x1)) + x3)) * 0.002816997) + 4.621029" +39,0.0019579437,"((((exp((x0 + 4.4453154) - (exp(x0 + ((x2 + x0) * 0.07406609)) * 0.17924055)) + (((x0 + (x7 + -168.80672)) * log(x1)) + x7)) * 1.2684214) + ((x2 - exp(x1)) + x3)) * 0.002816997) + 4.621029" +41,0.0019576072,"((((x3 + x2) - exp(x1)) + ((exp((x0 + 4.4453154) - (exp(x0 + (((x2 + 0.19918667) + x5) * 0.07406609)) * 0.17924055)) + ((((x0 + x7) + -168.80672) * log(x1)) + x7)) * 1.2684214)) * 0.002816997) + 4.621029" +46,0.0019009885,"(((((x3 - exp(x1)) + exp((x7 + 4.2632484) + (exp(x6 - (x2 * -0.044562757)) * (x2 * -0.011414634)))) + (x0 / (exp(x2) + 0.13353904))) + ((((x2 * 0.36960286) + x0) + (x6 + (x0 + -170.13193))) * log(x1))) * 0.0030565287) + 4.562415" +48,0.0018801361,"(((((x3 - exp(x1)) + exp((x7 + 4.2632484) + (exp(x6 - (x2 * -0.044562757)) * (x2 * -0.011414634)))) + (x0 / (exp(x2 + -1.1966766) + 0.13353904))) + ((((x2 * 0.36960286) + x0) + (x6 + (x0 + -170.13193))) * log(x1))) * 0.0030565287) + 4.562415" +51,0.0018582277,"((((exp(x7 + (((0.5140035 - exp(x7 - ((x2 + x3) * -0.021865297))) * (x2 * 0.016256034)) + 4.2675204)) - exp(x1)) + x3) + ((((x0 / log(exp(x2) + 1.1707261)) + (((x3 * 0.32195646) + -170.13176) + x4)) + (x0 + x7)) * log(x1))) * 0.0030460572) + 4.5678697" +52,0.0018558715,"((((x7 - exp(x1)) + ((exp((x7 + 4.2675204) + (((0.45004165 - exp(x7 - ((x2 + x2) * -0.021865297))) * x3) * 0.016256034)) + (x0 + ((x3 * 0.48458076) + x0))) + (log(x1) * (((x0 / (exp(x2) + 0.15762413)) + x0) + -170.13176)))) + x2) * 0.0030460572) + 4.5678697" +53,0.0018441706,"((((x7 - exp(x1)) + x2) + (exp(x7 + (((x3 * (0.45004165 - exp(x7 - ((x2 + x2) * -0.021865297)))) * 0.016256034) + 4.2675204)) + ((x0 + (x3 * 0.48458076)) + (x0 + (((x0 + (x0 / log(exp(x2) + 1.1707261))) + -170.13176) * log(x1)))))) * 0.0030460572) + 4.5678697" +55,0.0018352251,"((exp(((0.016256034 * (0.45004165 - exp(x7 - ((x3 + x2) * -0.021865297)))) * x3) + (x7 + 4.2675204)) + (((((x6 + (0.32195646 * (x0 + x3))) + (x0 + (x0 / log(exp(0.38444686 * x2) + 1.1707261)))) + -170.13176) * log(x1)) + ((x2 + x7) - exp(x1)))) * 0.0030460572) + 4.5678697" +56,0.001831732,"((((((((x0 + x3) * 0.32195646) + (x6 + (x0 + (x0 / (exp((x2 * 0.38444686) - 1.4578546) + 0.15762413))))) + -170.13176) * log(x1)) + ((x2 + x7) - exp(x1))) + exp((((0.45004165 - exp(x7 - ((x3 + x2) * -0.021865297))) * x3) * 0.016256034) + (x7 + 4.2675204))) * 0.0030460572) + 4.5678697" +57,0.0018231487,"4.5678697 + (((((x6 + (0.32195646 * (x0 + x3))) + ((x0 + (x0 / log(exp(0.38444686 * (x2 - 1.4578546)) + 1.1707261))) + -170.13176)) * log(x1)) + ((x2 + (x7 - exp(x1))) + exp(((0.45004165 - exp(x7 - ((x3 + x2) * -0.021865297))) * (x3 * 0.016256034)) + (x7 + 4.2675204)))) * 0.0030460572)" +59,0.0018210749,"4.5678697 + (((x2 + ((x7 - exp(x1)) + exp(((0.45004165 - exp(x7 - (-0.021865297 * (x3 + (x6 + x3))))) * (x3 * 0.016256034)) + (x7 + 4.2675204)))) + ((x6 + ((0.32195646 * (x4 + x3)) + ((x0 + (x0 / log(exp(0.38444686 * (x2 - 1.4578546)) + 1.1707261))) + -170.13176))) * log(x1))) * 0.0030460572)" +60,0.0017514622,"((exp((x0 + ((0.1897901 - x3) * (exp(x7 - (x2 * -0.05214908)) * 0.012381872))) + 4.257467) + ((x0 + ((x0 + x3) + (log(x1) * ((((x7 + (exp(x1) / ((exp(x1) * 0.11401008) + -1.615121))) + (x6 / (exp(x3) + 0.17297076))) + x7) + ((x3 * 0.40235698) + -170.31393))))) - exp(x1))) * 0.003012556) + 4.5001774" +61,0.0017312295,"(((((x2 + exp(((((1.7464508 - x3) * exp(x0 - (x2 * -0.051036775))) * 0.012590295) + x0) + 4.252308)) + x5) + (log(x1) * ((x3 * 0.3965832) + ((-170.31387 + x7) + ((exp(x1) / exp((exp(x1) * 0.09692516) + -1.614161)) + x7))))) + (((x6 / (exp(x2) + 0.14214376)) - exp(x1)) + x7)) * 0.0029462175) + 4.5256867" +62,0.0017311713,"4.5256867 + (((((x2 + exp((((exp(0.4734397) - x3) * (exp(x0 - (x2 * -0.051036775)) * 0.012590295)) + x0) + 4.252308)) + x5) + (((((x3 * 0.3965832) + -170.31387) + ((x7 + (exp(x1) / exp((exp(x1) * 0.09692516) + -1.614161))) + x7)) * log(x1)) + ((x6 / (exp(x2) + 0.14214376)) - exp(x1)))) + x7) * 0.0029462175)" +63,0.0017195458,"((exp((x0 + ((0.7426112 - x3) * (0.011224996 * exp(x4 - ((x2 + 1.6157049) * -0.04488468))))) + 4.239262) + (((x2 + x5) + (log(x1) * ((((exp(x1) / exp(-1.6128706 + (exp(x1) * 0.0907149))) + x7) + ((x3 * 0.41982338) + -170.31393)) + x0))) + (((x6 / (exp(x2) + 0.13077576)) - exp(x1)) + x7))) * 0.0029338102) + 4.5133047" +64,0.0017192991,"((exp((x0 + ((exp(0.17035629) - x3) * (0.011224996 * exp(x4 - ((x2 + 1.6157049) * -0.04488468))))) + 4.239262) + (((x2 + x5) + (log(x1) * (((((exp(x1) / exp(-1.6128706 + (exp(x1) * 0.0907149))) + x7) + (x3 * 0.41982338)) + -170.31393) + x0))) + (((x6 / (exp(x2) + 0.13077576)) - exp(x1)) + x7))) * 0.0029338102) + 4.5133047" +65,0.0016915298,"(((((((x3 + exp((x4 + (exp(x6 - ((x2 + x2) * -0.033423472)) * ((x1 - x3) * 0.010805121))) + 4.2677464)) - exp(x1)) + x7) + (log(x1) * ((((x3 * 0.40745232) + x0) + (x7 / (exp(x2) + 0.17825046))) + (x7 + -169.69832)))) + ((exp(x1) + x2) / (exp(exp(x1) * 0.1189319) * 0.16781318))) + x0) * 0.003010543) + 4.5433464" +67,0.0016782047,"((x7 + (((x2 + exp(x1)) / (exp(exp(x1) * 0.1189319) * 0.18821548)) + ((x0 + ((((x7 + ((x4 / (exp(x2) + 0.17825046)) + -169.69832)) + ((x3 * 0.40745232) + x0)) * log(x1)) + (x3 + exp((((exp(x6 - (((x2 - x1) + x3) * -0.033423472)) * (x1 - x2)) * 0.010805121) + x6) + 4.2677464)))) - exp(x1)))) * 0.003010543) + 4.5433464" +69,0.0016761663,"((x5 + ((((x2 / 1.3202298) + exp(x1)) / (exp(exp(x1) * 0.1189319) * 0.16781318)) + ((((x0 + x3) + exp((x6 + ((x1 - x3) * (exp(x6 - (((x2 - x1) + x3) * -0.033423472)) * 0.010805121))) + 4.2677464)) - exp(x1)) + ((((x7 + x5) + ((x3 * 0.40745232) + (x4 / (exp(x2) + 0.20448817)))) + -169.69832) * log(x1))))) * 0.003010543) + 4.5433464" diff --git a/outputs/20260606_005734_6FlOfR/hall_of_fame.csv.bak b/outputs/20260606_005734_6FlOfR/hall_of_fame.csv.bak new file mode 100644 index 0000000..e1b00ba --- /dev/null +++ b/outputs/20260606_005734_6FlOfR/hall_of_fame.csv.bak @@ -0,0 +1,44 @@ +Complexity,Loss,Equation +1,0.1415826,"3.565013" +4,0.046780355,"log(153.24149 / x1)" +5,0.02017629,"(x1 * -0.41636097) + 5.4033523" +6,0.0149874715,"(exp(x1) * -0.0039443914) + 4.018696" +7,0.014847774,"(-10.5967 / x1) + (10.474295 - x1)" +8,0.012612821,"((x0 - exp(x1)) * 0.003974942) + 4.032723" +10,0.011021733,"(x7 * 0.011055068) + (4.0363736 - (exp(x1) * 0.0038689962))" +11,0.01002806,"(3.865385 - (exp(x1) * 0.003869981)) + (exp(x7) * 0.2532472)" +12,0.007341251,"((x7 * 0.016047547) + 4.005379) + ((x2 - exp(x1)) * 0.00402841)" +13,0.0069549805,"exp(x7 * 0.022573786) + (((x2 - exp(x1)) * 0.003997351) + 3.0122573)" +14,0.0053718733,"((((x1 * (x0 + -37.00079)) + x2) - exp(x1)) * 0.0030013788) + 4.389338" +16,0.004099731,"(((x2 - exp(x1)) * 0.0031908343) + 3.430174) + (2.2044883 / (x1 + (x7 * -0.49608624)))" +18,0.0039598015,"(((x2 * 0.8246542) - exp(x1)) * 0.0031952073) + ((2.2756846 / (x1 + (x7 * -0.4195139))) + 3.4163985)" +19,0.003771442,"(((exp(x7 + 4.0886) + (x1 * (x0 + -64.57939))) + (x3 - exp(x1))) * 0.002504624) + 4.4716077" +21,0.0034251248,"4.4838834 + (0.0025647841 * ((x2 + exp(x7 + 4.1382027)) + (((x1 * -64.71305) - exp(x1)) + (x4 * 3.854656))))" +22,0.0032816748,"((((log(x1) * -229.28015) + ((x0 / 0.26969364) + (x3 - exp(x1)))) + exp(x7 + 4.1170673)) * 0.0027460016) + 4.6925445" +24,0.0032803074,"((((x0 / 0.26969364) + exp(x7 + 4.1170673)) + ((-0.23561807 + (x3 - exp(x1))) + (log(x1) * -229.28015))) * 0.0027460016) + 4.6925445" +26,0.0031384756,"((((exp((x7 + 4.129333) - (x2 * 0.0056593237)) - exp(x1)) + (x2 + (log(x1) * -169.4319))) + (x0 * 3.3012187)) * 0.003012371) + 4.537966" +28,0.0026985996,"((exp(((x3 * x3) * -0.00079846877) + (x7 - -4.3820744)) + (((((x5 * 2.9971557) + x3) + -171.94554) * log(x1)) - exp(x1))) * 0.0030115678) + 4.538426" +30,0.002533871,"((exp((x3 * ((x3 + 9.939942) * -0.0006410332)) + (x7 - -4.3820744)) + (((((x0 * 3.0800557) + x3) + -171.76453) * log(x1)) - exp(x1))) * 0.0030115678) + 4.538426" +32,0.0025038442,"(((exp((((-5.920857 - x3) * x3) * 0.00073056546) + (x0 + 4.430383)) + ((((x7 + (x3 + -169.90005)) + x4) + x7) * log(x1))) - exp(x1)) * 0.003055337) + 4.5436344" +33,0.0021212408,"(((log(x1) * (((x5 + x7) / 0.58209914) + -166.40523)) + (x2 + ((x3 - exp(x1)) + exp((4.6342545 - (exp(x2 * 0.09402527) * 0.082178235)) + x0)))) * 0.0030390304) + 4.498984" +35,0.002008648,"((((x3 - exp(x1)) + ((x7 + (((x6 + -169.14818) + x0) * log(x1))) * 1.3442688)) + (exp((x0 + 4.5737853) - (exp(x2 * 0.095374905) * 0.07284515)) + x3)) * 0.002771325) + 4.6696644" +37,0.0019595241,"((((exp(x0 + (4.4453154 - (exp(x0 + (x2 * 0.07406609)) * 0.17924055))) + ((((x0 + x7) + -168.80672) * log(x1)) + x7)) * 1.2684214) + ((x2 - exp(x1)) + x3)) * 0.002816997) + 4.621029" +39,0.0019579437,"((((exp((x0 + 4.4453154) - (exp(x0 + ((x2 + x0) * 0.07406609)) * 0.17924055)) + (((x0 + (x7 + -168.80672)) * log(x1)) + x7)) * 1.2684214) + ((x2 - exp(x1)) + x3)) * 0.002816997) + 4.621029" +41,0.0019576072,"((((x3 + x2) - exp(x1)) + ((exp((x0 + 4.4453154) - (exp(x0 + (((x2 + 0.19918667) + x5) * 0.07406609)) * 0.17924055)) + ((((x0 + x7) + -168.80672) * log(x1)) + x7)) * 1.2684214)) * 0.002816997) + 4.621029" +46,0.0019009885,"(((((x3 - exp(x1)) + exp((x7 + 4.2632484) + (exp(x6 - (x2 * -0.044562757)) * (x2 * -0.011414634)))) + (x0 / (exp(x2) + 0.13353904))) + ((((x2 * 0.36960286) + x0) + (x6 + (x0 + -170.13193))) * log(x1))) * 0.0030565287) + 4.562415" +48,0.0018801361,"(((((x3 - exp(x1)) + exp((x7 + 4.2632484) + (exp(x6 - (x2 * -0.044562757)) * (x2 * -0.011414634)))) + (x0 / (exp(x2 + -1.1966766) + 0.13353904))) + ((((x2 * 0.36960286) + x0) + (x6 + (x0 + -170.13193))) * log(x1))) * 0.0030565287) + 4.562415" +51,0.0018582277,"((((exp(x7 + (((0.5140035 - exp(x7 - ((x2 + x3) * -0.021865297))) * (x2 * 0.016256034)) + 4.2675204)) - exp(x1)) + x3) + ((((x0 / log(exp(x2) + 1.1707261)) + (((x3 * 0.32195646) + -170.13176) + x4)) + (x0 + x7)) * log(x1))) * 0.0030460572) + 4.5678697" +52,0.0018558715,"((((x7 - exp(x1)) + ((exp((x7 + 4.2675204) + (((0.45004165 - exp(x7 - ((x2 + x2) * -0.021865297))) * x3) * 0.016256034)) + (x0 + ((x3 * 0.48458076) + x0))) + (log(x1) * (((x0 / (exp(x2) + 0.15762413)) + x0) + -170.13176)))) + x2) * 0.0030460572) + 4.5678697" +53,0.0018441706,"((((x7 - exp(x1)) + x2) + (exp(x7 + (((x3 * (0.45004165 - exp(x7 - ((x2 + x2) * -0.021865297)))) * 0.016256034) + 4.2675204)) + ((x0 + (x3 * 0.48458076)) + (x0 + (((x0 + (x0 / log(exp(x2) + 1.1707261))) + -170.13176) * log(x1)))))) * 0.0030460572) + 4.5678697" +55,0.0018352251,"((exp(((0.016256034 * (0.45004165 - exp(x7 - ((x3 + x2) * -0.021865297)))) * x3) + (x7 + 4.2675204)) + (((((x6 + (0.32195646 * (x0 + x3))) + (x0 + (x0 / log(exp(0.38444686 * x2) + 1.1707261)))) + -170.13176) * log(x1)) + ((x2 + x7) - exp(x1)))) * 0.0030460572) + 4.5678697" +56,0.001831732,"((((((((x0 + x3) * 0.32195646) + (x6 + (x0 + (x0 / (exp((x2 * 0.38444686) - 1.4578546) + 0.15762413))))) + -170.13176) * log(x1)) + ((x2 + x7) - exp(x1))) + exp((((0.45004165 - exp(x7 - ((x3 + x2) * -0.021865297))) * x3) * 0.016256034) + (x7 + 4.2675204))) * 0.0030460572) + 4.5678697" +57,0.0018231487,"4.5678697 + (((((x6 + (0.32195646 * (x0 + x3))) + ((x0 + (x0 / log(exp(0.38444686 * (x2 - 1.4578546)) + 1.1707261))) + -170.13176)) * log(x1)) + ((x2 + (x7 - exp(x1))) + exp(((0.45004165 - exp(x7 - ((x3 + x2) * -0.021865297))) * (x3 * 0.016256034)) + (x7 + 4.2675204)))) * 0.0030460572)" +59,0.0018210749,"4.5678697 + (((x2 + ((x7 - exp(x1)) + exp(((0.45004165 - exp(x7 - (-0.021865297 * (x3 + (x6 + x3))))) * (x3 * 0.016256034)) + (x7 + 4.2675204)))) + ((x6 + ((0.32195646 * (x4 + x3)) + ((x0 + (x0 / log(exp(0.38444686 * (x2 - 1.4578546)) + 1.1707261))) + -170.13176))) * log(x1))) * 0.0030460572)" +60,0.0017514622,"((exp((x0 + ((0.1897901 - x3) * (exp(x7 - (x2 * -0.05214908)) * 0.012381872))) + 4.257467) + ((x0 + ((x0 + x3) + (log(x1) * ((((x7 + (exp(x1) / ((exp(x1) * 0.11401008) + -1.615121))) + (x6 / (exp(x3) + 0.17297076))) + x7) + ((x3 * 0.40235698) + -170.31393))))) - exp(x1))) * 0.003012556) + 4.5001774" +61,0.0017312295,"(((((x2 + exp(((((1.7464508 - x3) * exp(x0 - (x2 * -0.051036775))) * 0.012590295) + x0) + 4.252308)) + x5) + (log(x1) * ((x3 * 0.3965832) + ((-170.31387 + x7) + ((exp(x1) / exp((exp(x1) * 0.09692516) + -1.614161)) + x7))))) + (((x6 / (exp(x2) + 0.14214376)) - exp(x1)) + x7)) * 0.0029462175) + 4.5256867" +62,0.0017311713,"4.5256867 + (((((x2 + exp((((exp(0.4734397) - x3) * (exp(x0 - (x2 * -0.051036775)) * 0.012590295)) + x0) + 4.252308)) + x5) + (((((x3 * 0.3965832) + -170.31387) + ((x7 + (exp(x1) / exp((exp(x1) * 0.09692516) + -1.614161))) + x7)) * log(x1)) + ((x6 / (exp(x2) + 0.14214376)) - exp(x1)))) + x7) * 0.0029462175)" +63,0.0017195458,"((exp((x0 + ((0.7426112 - x3) * (0.011224996 * exp(x4 - ((x2 + 1.6157049) * -0.04488468))))) + 4.239262) + (((x2 + x5) + (log(x1) * ((((exp(x1) / exp(-1.6128706 + (exp(x1) * 0.0907149))) + x7) + ((x3 * 0.41982338) + -170.31393)) + x0))) + (((x6 / (exp(x2) + 0.13077576)) - exp(x1)) + x7))) * 0.0029338102) + 4.5133047" +64,0.0017192991,"((exp((x0 + ((exp(0.17035629) - x3) * (0.011224996 * exp(x4 - ((x2 + 1.6157049) * -0.04488468))))) + 4.239262) + (((x2 + x5) + (log(x1) * (((((exp(x1) / exp(-1.6128706 + (exp(x1) * 0.0907149))) + x7) + (x3 * 0.41982338)) + -170.31393) + x0))) + (((x6 / (exp(x2) + 0.13077576)) - exp(x1)) + x7))) * 0.0029338102) + 4.5133047" +65,0.0016915298,"(((((((x3 + exp((x4 + (exp(x6 - ((x2 + x2) * -0.033423472)) * ((x1 - x3) * 0.010805121))) + 4.2677464)) - exp(x1)) + x7) + (log(x1) * ((((x3 * 0.40745232) + x0) + (x7 / (exp(x2) + 0.17825046))) + (x7 + -169.69832)))) + ((exp(x1) + x2) / (exp(exp(x1) * 0.1189319) * 0.16781318))) + x0) * 0.003010543) + 4.5433464" +67,0.0016782047,"((x7 + (((x2 + exp(x1)) / (exp(exp(x1) * 0.1189319) * 0.18821548)) + ((x0 + ((((x7 + ((x4 / (exp(x2) + 0.17825046)) + -169.69832)) + ((x3 * 0.40745232) + x0)) * log(x1)) + (x3 + exp((((exp(x6 - (((x2 - x1) + x3) * -0.033423472)) * (x1 - x2)) * 0.010805121) + x6) + 4.2677464)))) - exp(x1)))) * 0.003010543) + 4.5433464" +69,0.0016761663,"((x5 + ((((x2 / 1.3202298) + exp(x1)) / (exp(exp(x1) * 0.1189319) * 0.16781318)) + ((((x0 + x3) + exp((x6 + ((x1 - x3) * (exp(x6 - (((x2 - x1) + x3) * -0.033423472)) * 0.010805121))) + 4.2677464)) - exp(x1)) + ((((x7 + x5) + ((x3 * 0.40745232) + (x4 / (exp(x2) + 0.20448817)))) + -169.69832) * log(x1))))) * 0.003010543) + 4.5433464" diff --git a/outputs/20260606_021144_8PzK7u/checkpoint.pkl b/outputs/20260606_021144_8PzK7u/checkpoint.pkl new file mode 100644 index 0000000..073f2ff Binary files /dev/null and b/outputs/20260606_021144_8PzK7u/checkpoint.pkl differ diff --git a/outputs/20260606_021144_8PzK7u/hall_of_fame.csv b/outputs/20260606_021144_8PzK7u/hall_of_fame.csv new file mode 100644 index 0000000..c030c2c --- /dev/null +++ b/outputs/20260606_021144_8PzK7u/hall_of_fame.csv @@ -0,0 +1,16 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5650115" +4,0.04678043,"5.0317802 - log(x1)" +5,0.020174915,"(x1 * -0.41540056) + 5.4004316" +6,0.013971531,"log((x1 * -14.91161) + 103.47175)" +8,0.013313121,"log((x1 * (x1 * -1.5824497)) + 69.60359)" +9,0.010057912,"log(((exp(x4) + -15.226362) * x1) + 102.3669)" +10,0.009091249,"(x0 * 0.011894782) + log((x1 * -15.442862) + 107.0882)" +12,0.005343384,"log((x1 * -15.528032) + 108.428955) + (x4 / (x2 + 28.107822))" +14,0.003161533,"(x4 / (x2 + (21.339392 - x0))) + log((x1 * -15.78025) + 109.661316)" +16,0.0030244035,"(x4 / (x2 + (21.316158 - x0))) + (log((x1 * -15.71102) + 108.384575) - -0.033276055)" +18,0.0029697453,"((x4 / (((x3 - x0) + 22.84757) / 1.2851638)) - -0.033057947) + log((x1 * -15.678972) + 108.401146)" +20,0.0026580312,"((x4 / ((((x2 + 25.207476) * 0.6815766) - x0) + -1.9441698)) + log((x1 * -15.755431) + 108.370316)) + 0.056702785" +21,0.0026149554,"(x4 / ((((x2 + 24.162056) * 0.6815766) - x0) - log(x1))) + (log((x1 * -15.755431) + 108.370316) + 0.056702785)" +23,0.0026113451,"(x4 / (((((x2 + 25.207476) * 0.6815766) - x0) + -0.8020641) - log(x1))) + (log((x1 * -15.755431) + 108.370316) + 0.056702785)" +24,0.0025531421,"(log((x1 * -15.71263) + 108.39293) + (((x4 / 0.75595826) * 1.3551189) / ((x2 - ((x0 + 1.4032606) / 0.40380266)) + 25.639969))) - -0.052671567" diff --git a/outputs/20260606_021144_8PzK7u/hall_of_fame.csv.bak b/outputs/20260606_021144_8PzK7u/hall_of_fame.csv.bak new file mode 100644 index 0000000..c030c2c --- /dev/null +++ b/outputs/20260606_021144_8PzK7u/hall_of_fame.csv.bak @@ -0,0 +1,16 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5650115" +4,0.04678043,"5.0317802 - log(x1)" +5,0.020174915,"(x1 * -0.41540056) + 5.4004316" +6,0.013971531,"log((x1 * -14.91161) + 103.47175)" +8,0.013313121,"log((x1 * (x1 * -1.5824497)) + 69.60359)" +9,0.010057912,"log(((exp(x4) + -15.226362) * x1) + 102.3669)" +10,0.009091249,"(x0 * 0.011894782) + log((x1 * -15.442862) + 107.0882)" +12,0.005343384,"log((x1 * -15.528032) + 108.428955) + (x4 / (x2 + 28.107822))" +14,0.003161533,"(x4 / (x2 + (21.339392 - x0))) + log((x1 * -15.78025) + 109.661316)" +16,0.0030244035,"(x4 / (x2 + (21.316158 - x0))) + (log((x1 * -15.71102) + 108.384575) - -0.033276055)" +18,0.0029697453,"((x4 / (((x3 - x0) + 22.84757) / 1.2851638)) - -0.033057947) + log((x1 * -15.678972) + 108.401146)" +20,0.0026580312,"((x4 / ((((x2 + 25.207476) * 0.6815766) - x0) + -1.9441698)) + log((x1 * -15.755431) + 108.370316)) + 0.056702785" +21,0.0026149554,"(x4 / ((((x2 + 24.162056) * 0.6815766) - x0) - log(x1))) + (log((x1 * -15.755431) + 108.370316) + 0.056702785)" +23,0.0026113451,"(x4 / (((((x2 + 25.207476) * 0.6815766) - x0) + -0.8020641) - log(x1))) + (log((x1 * -15.755431) + 108.370316) + 0.056702785)" +24,0.0025531421,"(log((x1 * -15.71263) + 108.39293) + (((x4 / 0.75595826) * 1.3551189) / ((x2 - ((x0 + 1.4032606) / 0.40380266)) + 25.639969))) - -0.052671567" diff --git a/outputs/20260606_021836_7T4wfI/checkpoint.pkl b/outputs/20260606_021836_7T4wfI/checkpoint.pkl new file mode 100644 index 0000000..27f91a5 Binary files /dev/null and b/outputs/20260606_021836_7T4wfI/checkpoint.pkl differ diff --git a/outputs/20260606_021836_7T4wfI/hall_of_fame.csv b/outputs/20260606_021836_7T4wfI/hall_of_fame.csv new file mode 100644 index 0000000..76d2730 --- /dev/null +++ b/outputs/20260606_021836_7T4wfI/hall_of_fame.csv @@ -0,0 +1,18 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5649922" +5,0.020174954,"(x1 * -0.4158265) + 5.401702" +6,0.013977423,"log((x1 + -6.931311) * -14.925834)" +8,0.013115806,"log((x1 + -6.355604) * -95.999115) * 0.69748753" +9,0.011385535,"log((x1 + -6.9157047) * (-14.414061 - exp(x4)))" +10,0.009269198,"log(((x4 * -0.13972968) + -14.9506855) * (x1 + -6.980082))" +11,0.0076216357,"log((log(0.28418428 - x0) + -15.134965) * (x1 + -6.9323688))" +13,0.007465354,"log((x1 + -6.980082) * (log((x4 * -1.8015646) + 0.3693332) + -15.188278))" +16,0.005071091,"log(((-15.278458 - ((x2 + x4) * 0.21522282)) + (x3 * 0.17450936)) * (x1 + -6.886979))" +18,0.005066207,"log((((x3 * 0.15619037) - (((x2 + -4.363537) + x4) * 0.20172201)) + -15.7568445) * (x1 + -6.9177847))" +20,0.004830611,"log((x1 + -6.8552194) * (((x3 * 0.17556405) + -15.992225) - ((((x1 + x2) + -6.8552194) + x4) * 0.22224921)))" +21,0.004592707,"log(((x1 + -6.886094) * (((x2 * 0.1518889) + -15.748936) - (((x4 + -4.365079) + x2) * 0.19550365))) + exp(x0))" +23,0.003510716,"log((((0.096697345 * x2) + -15.619484) - ((exp(x4 + 0.89043003) + -1.9077631) + ((x2 + x4) * 0.14370511))) * (x1 + -6.864843))" +25,0.003488224,"log((x1 + -6.864843) * (((x3 + 0.3348677) * 0.096697345) + (-15.619484 - ((exp(x0 - -0.93297875) + -1.9077631) + (0.14370511 * (x2 + x4))))))" +26,0.0034845257,"log((((x2 * 0.096697345) + -15.619484) - ((exp(x4 + 0.89043003) + -1.9077631) + (((exp(x4) + x2) + x4) * 0.14370511))) * (x1 + -6.864843))" +27,0.0033134876,"log((x1 + -6.8515253) * (((x2 * 0.09947497) + -15.76245) - (((((x4 + x2) + -6.8497934) + x1) * 0.14530487) + (exp(x4 + 1.0534714) + -1.7831411))))" +29,0.0033016787,"log((x1 + -6.8515253) * (((x2 * 0.09947497) + -15.74374) - ((((x1 + ((x4 + -6.8497934) + x2)) * 0.14530487) + exp((x4 * 0.90974265) + 1.0243999)) + -1.7831411)))" diff --git a/outputs/20260606_021836_7T4wfI/hall_of_fame.csv.bak b/outputs/20260606_021836_7T4wfI/hall_of_fame.csv.bak new file mode 100644 index 0000000..76d2730 --- /dev/null +++ b/outputs/20260606_021836_7T4wfI/hall_of_fame.csv.bak @@ -0,0 +1,18 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5649922" +5,0.020174954,"(x1 * -0.4158265) + 5.401702" +6,0.013977423,"log((x1 + -6.931311) * -14.925834)" +8,0.013115806,"log((x1 + -6.355604) * -95.999115) * 0.69748753" +9,0.011385535,"log((x1 + -6.9157047) * (-14.414061 - exp(x4)))" +10,0.009269198,"log(((x4 * -0.13972968) + -14.9506855) * (x1 + -6.980082))" +11,0.0076216357,"log((log(0.28418428 - x0) + -15.134965) * (x1 + -6.9323688))" +13,0.007465354,"log((x1 + -6.980082) * (log((x4 * -1.8015646) + 0.3693332) + -15.188278))" +16,0.005071091,"log(((-15.278458 - ((x2 + x4) * 0.21522282)) + (x3 * 0.17450936)) * (x1 + -6.886979))" +18,0.005066207,"log((((x3 * 0.15619037) - (((x2 + -4.363537) + x4) * 0.20172201)) + -15.7568445) * (x1 + -6.9177847))" +20,0.004830611,"log((x1 + -6.8552194) * (((x3 * 0.17556405) + -15.992225) - ((((x1 + x2) + -6.8552194) + x4) * 0.22224921)))" +21,0.004592707,"log(((x1 + -6.886094) * (((x2 * 0.1518889) + -15.748936) - (((x4 + -4.365079) + x2) * 0.19550365))) + exp(x0))" +23,0.003510716,"log((((0.096697345 * x2) + -15.619484) - ((exp(x4 + 0.89043003) + -1.9077631) + ((x2 + x4) * 0.14370511))) * (x1 + -6.864843))" +25,0.003488224,"log((x1 + -6.864843) * (((x3 + 0.3348677) * 0.096697345) + (-15.619484 - ((exp(x0 - -0.93297875) + -1.9077631) + (0.14370511 * (x2 + x4))))))" +26,0.0034845257,"log((((x2 * 0.096697345) + -15.619484) - ((exp(x4 + 0.89043003) + -1.9077631) + (((exp(x4) + x2) + x4) * 0.14370511))) * (x1 + -6.864843))" +27,0.0033134876,"log((x1 + -6.8515253) * (((x2 * 0.09947497) + -15.76245) - (((((x4 + x2) + -6.8497934) + x1) * 0.14530487) + (exp(x4 + 1.0534714) + -1.7831411))))" +29,0.0033016787,"log((x1 + -6.8515253) * (((x2 * 0.09947497) + -15.74374) - ((((x1 + ((x4 + -6.8497934) + x2)) * 0.14530487) + exp((x4 * 0.90974265) + 1.0243999)) + -1.7831411)))" diff --git a/outputs/20260606_022118_Q2Uhrd/checkpoint.pkl b/outputs/20260606_022118_Q2Uhrd/checkpoint.pkl new file mode 100644 index 0000000..02e9472 Binary files /dev/null and b/outputs/20260606_022118_Q2Uhrd/checkpoint.pkl differ diff --git a/outputs/20260606_022118_Q2Uhrd/hall_of_fame.csv b/outputs/20260606_022118_Q2Uhrd/hall_of_fame.csv new file mode 100644 index 0000000..3771fb8 --- /dev/null +++ b/outputs/20260606_022118_Q2Uhrd/hall_of_fame.csv @@ -0,0 +1,17 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5651019" +4,0.04678116,"5.031143 - log(x1)" +5,0.020179946,"(x1 * -0.41746596) - -5.407585" +6,0.014148784,"log((x1 * -14.363176) + 100.754456)" +8,0.013976877,"log(((x1 * -1.6267196) + 11.274463) / 0.10894275)" +9,0.013743353,"log(exp(x4) + ((x1 * -13.566513) + 95.94948))" +11,0.009203923,"log((exp(x4) + ((x1 * -2.1340098) + 14.426924)) / 0.14835517)" +13,0.009199121,"log(((exp(x4) + 14.432163) + (x1 * -2.140572)) / 0.06682556) - 0.7912578" +15,0.0071718534,"log(((exp(x4) + (12.411925 + (x1 * -1.880818))) / 0.12609364) + (x2 / 17.89822))" +17,0.006629608,"log((((exp(x4) + 12.411925) + (x1 * -1.880818)) / 0.12609364) + ((x2 + x0) / 17.89822))" +18,0.0059555923,"log(((exp(x4) + ((x0 + x2) / exp(x1))) + (12.416725 + (-1.8421992 * x1))) / 0.13347152)" +20,0.0050544622,"log((exp(x4) + (((x4 + (x4 + x3)) / exp(x1)) + (12.480034 + (x1 * -1.8591905)))) / 0.13062614)" +22,0.0050447076,"log((exp(x4) + (((((x0 + -0.38507852) + x4) + x3) / exp(x1)) + (12.480034 + (x1 * -1.8591905)))) / 0.13062614)" +23,0.0049069854,"log(((((x4 + (x2 + x4)) / exp(x1 + exp(x4))) + (x1 * -1.8647327)) + (exp(x4) + 12.472427)) / 0.13161501)" +24,0.0046067215,"log(((((x1 * -1.8556168) + 12.429903) + exp(x0)) / 0.13109726) + ((x1 * (x4 + (x0 + x3))) / exp(x1 + -0.6515075)))" +25,0.0041190893,"log((((x1 * -1.8995988) + (exp(x4) + 12.467603)) / 0.12689108) + (((x2 + (x0 + x4)) / (x1 * exp(exp(x4)))) + -0.028780995))" diff --git a/outputs/20260606_022118_Q2Uhrd/hall_of_fame.csv.bak b/outputs/20260606_022118_Q2Uhrd/hall_of_fame.csv.bak new file mode 100644 index 0000000..3771fb8 --- /dev/null +++ b/outputs/20260606_022118_Q2Uhrd/hall_of_fame.csv.bak @@ -0,0 +1,17 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5651019" +4,0.04678116,"5.031143 - log(x1)" +5,0.020179946,"(x1 * -0.41746596) - -5.407585" +6,0.014148784,"log((x1 * -14.363176) + 100.754456)" +8,0.013976877,"log(((x1 * -1.6267196) + 11.274463) / 0.10894275)" +9,0.013743353,"log(exp(x4) + ((x1 * -13.566513) + 95.94948))" +11,0.009203923,"log((exp(x4) + ((x1 * -2.1340098) + 14.426924)) / 0.14835517)" +13,0.009199121,"log(((exp(x4) + 14.432163) + (x1 * -2.140572)) / 0.06682556) - 0.7912578" +15,0.0071718534,"log(((exp(x4) + (12.411925 + (x1 * -1.880818))) / 0.12609364) + (x2 / 17.89822))" +17,0.006629608,"log((((exp(x4) + 12.411925) + (x1 * -1.880818)) / 0.12609364) + ((x2 + x0) / 17.89822))" +18,0.0059555923,"log(((exp(x4) + ((x0 + x2) / exp(x1))) + (12.416725 + (-1.8421992 * x1))) / 0.13347152)" +20,0.0050544622,"log((exp(x4) + (((x4 + (x4 + x3)) / exp(x1)) + (12.480034 + (x1 * -1.8591905)))) / 0.13062614)" +22,0.0050447076,"log((exp(x4) + (((((x0 + -0.38507852) + x4) + x3) / exp(x1)) + (12.480034 + (x1 * -1.8591905)))) / 0.13062614)" +23,0.0049069854,"log(((((x4 + (x2 + x4)) / exp(x1 + exp(x4))) + (x1 * -1.8647327)) + (exp(x4) + 12.472427)) / 0.13161501)" +24,0.0046067215,"log(((((x1 * -1.8556168) + 12.429903) + exp(x0)) / 0.13109726) + ((x1 * (x4 + (x0 + x3))) / exp(x1 + -0.6515075)))" +25,0.0041190893,"log((((x1 * -1.8995988) + (exp(x4) + 12.467603)) / 0.12689108) + (((x2 + (x0 + x4)) / (x1 * exp(exp(x4)))) + -0.028780995))" diff --git a/outputs/20260606_022206_B4MMfw/checkpoint.pkl b/outputs/20260606_022206_B4MMfw/checkpoint.pkl new file mode 100644 index 0000000..05de515 Binary files /dev/null and b/outputs/20260606_022206_B4MMfw/checkpoint.pkl differ diff --git a/outputs/20260606_022206_B4MMfw/hall_of_fame.csv b/outputs/20260606_022206_B4MMfw/hall_of_fame.csv new file mode 100644 index 0000000..556dede --- /dev/null +++ b/outputs/20260606_022206_B4MMfw/hall_of_fame.csv @@ -0,0 +1,14 @@ +Complexity,Loss,Equation +1,0.14158265,"3.5647438" +4,0.047093198,"5.0497274 - log(x1)" +5,0.020182258,"(x1 * -0.4186744) + 5.413453" +6,0.014989613,"4.021728 - (exp(x1) * 0.003973103)" +8,0.012648692,"4.0218387 - ((exp(x1) - x4) * 0.0038953868)" +10,0.011386385,"4.035525 - (((exp(x1) - x4) - x4) * 0.0038734567)" +11,0.011178298,"4.0152345 - ((exp(x1) - (x2 * exp(x4))) * 0.0041016773)" +12,0.007562292,"4.0141654 - (0.0040170457 * (exp(x1) - (x2 + (x4 * x1))))" +14,0.0069505856,"4.0141654 - ((exp(x1) - ((x2 * 0.7400101) + (x4 * 3.2431555))) * 0.0040170457)" +16,0.006916721,"4.021728 - ((exp(x1) - ((x2 * 0.6403898) + (x4 - (x4 * -2.706216)))) * 0.003973103)" +18,0.006759502,"4.029247 - ((exp(x1) - ((x4 + (((x1 - x2) - x4) * -0.2767068)) * 2.848128)) * 0.0040198658)" +21,0.0067123845,"4.029247 - ((exp(x1) - ((x4 + (((x1 - (x2 + exp(x4))) - x4) * -0.2767068)) * 2.848128)) * 0.0040198658)" +23,0.006687721,"4.029247 - ((exp(x1) - ((x4 + (((x1 - (x2 + exp(0.47565556 + x4))) - x4) * -0.2767068)) * 2.848128)) * 0.0040198658)" diff --git a/outputs/20260606_022206_B4MMfw/hall_of_fame.csv.bak b/outputs/20260606_022206_B4MMfw/hall_of_fame.csv.bak new file mode 100644 index 0000000..556dede --- /dev/null +++ b/outputs/20260606_022206_B4MMfw/hall_of_fame.csv.bak @@ -0,0 +1,14 @@ +Complexity,Loss,Equation +1,0.14158265,"3.5647438" +4,0.047093198,"5.0497274 - log(x1)" +5,0.020182258,"(x1 * -0.4186744) + 5.413453" +6,0.014989613,"4.021728 - (exp(x1) * 0.003973103)" +8,0.012648692,"4.0218387 - ((exp(x1) - x4) * 0.0038953868)" +10,0.011386385,"4.035525 - (((exp(x1) - x4) - x4) * 0.0038734567)" +11,0.011178298,"4.0152345 - ((exp(x1) - (x2 * exp(x4))) * 0.0041016773)" +12,0.007562292,"4.0141654 - (0.0040170457 * (exp(x1) - (x2 + (x4 * x1))))" +14,0.0069505856,"4.0141654 - ((exp(x1) - ((x2 * 0.7400101) + (x4 * 3.2431555))) * 0.0040170457)" +16,0.006916721,"4.021728 - ((exp(x1) - ((x2 * 0.6403898) + (x4 - (x4 * -2.706216)))) * 0.003973103)" +18,0.006759502,"4.029247 - ((exp(x1) - ((x4 + (((x1 - x2) - x4) * -0.2767068)) * 2.848128)) * 0.0040198658)" +21,0.0067123845,"4.029247 - ((exp(x1) - ((x4 + (((x1 - (x2 + exp(x4))) - x4) * -0.2767068)) * 2.848128)) * 0.0040198658)" +23,0.006687721,"4.029247 - ((exp(x1) - ((x4 + (((x1 - (x2 + exp(0.47565556 + x4))) - x4) * -0.2767068)) * 2.848128)) * 0.0040198658)" diff --git a/outputs/20260606_022315_s3ZViI/checkpoint.pkl b/outputs/20260606_022315_s3ZViI/checkpoint.pkl new file mode 100644 index 0000000..fcb4043 Binary files /dev/null and b/outputs/20260606_022315_s3ZViI/checkpoint.pkl differ diff --git a/outputs/20260606_022315_s3ZViI/hall_of_fame.csv b/outputs/20260606_022315_s3ZViI/hall_of_fame.csv new file mode 100644 index 0000000..66a058b --- /dev/null +++ b/outputs/20260606_022315_s3ZViI/hall_of_fame.csv @@ -0,0 +1,16 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5648918" +4,0.0467827,"log(153.47993 / x1)" +5,0.020175997,"(x1 * -0.41452923) + 5.3960147" +6,0.013988713,"log((x1 - 6.913469) / -0.06633385)" +8,0.0139727015,"log(((x1 - 6.972501) / -0.06720227) + -0.42802617)" +9,0.012775202,"log(exp(x0) + ((x1 - 6.894778) / -0.067534))" +10,0.010733719,"log(((x1 - 6.972501) / -0.0687157) + (x0 / x1))" +12,0.0107153375,"log(((x1 - 6.972501) / -0.0687157) + (x4 / (x1 + -0.0687157)))" +13,0.008655685,"log((x1 - 7.067938) * ((exp(x4) + x1) / (x1 * -0.08008032)))" +15,0.008432622,"log(((exp(x0 / 2.1213145) + x1) * (x1 - 7.18395)) / (x1 * -0.08667538))" +17,0.008372297,"log(((exp(x0 / 2.1409464) + x1) * (x1 - 7.1519527)) / (x1 * -0.11701374)) - -0.31662118" +19,0.007871232,"log(((exp(x0 / 2.1483574) + x1) * (x1 - 6.9429812)) / (x1 * -0.091777064)) - (x1 * -0.03532463)" +21,0.0076451213,"log((x1 + exp((x0 + (x4 / 0.40371257)) * 0.08877082)) * ((x1 - 6.7461934) / ((x1 * -0.045270775) - 0.11868309)))" +23,0.0069678803,"log((exp(x4) + ((0.03762332 / (x1 / (x3 + 0.008522787))) + x1)) * (((x1 - 7.850765) / (x1 * -0.09465818)) - 0.9733152))" +26,0.006417365,"(log((x1 * (x1 - 7.2052436)) / (x1 * -0.12828495)) + (((0.38349107 - ((x1 - (x2 / x1)) * 0.66997904)) + x0) * 0.019127231)) - -0.5852694" diff --git a/outputs/20260606_022315_s3ZViI/hall_of_fame.csv.bak b/outputs/20260606_022315_s3ZViI/hall_of_fame.csv.bak new file mode 100644 index 0000000..66a058b --- /dev/null +++ b/outputs/20260606_022315_s3ZViI/hall_of_fame.csv.bak @@ -0,0 +1,16 @@ +Complexity,Loss,Equation +1,0.1415826,"3.5648918" +4,0.0467827,"log(153.47993 / x1)" +5,0.020175997,"(x1 * -0.41452923) + 5.3960147" +6,0.013988713,"log((x1 - 6.913469) / -0.06633385)" +8,0.0139727015,"log(((x1 - 6.972501) / -0.06720227) + -0.42802617)" +9,0.012775202,"log(exp(x0) + ((x1 - 6.894778) / -0.067534))" +10,0.010733719,"log(((x1 - 6.972501) / -0.0687157) + (x0 / x1))" +12,0.0107153375,"log(((x1 - 6.972501) / -0.0687157) + (x4 / (x1 + -0.0687157)))" +13,0.008655685,"log((x1 - 7.067938) * ((exp(x4) + x1) / (x1 * -0.08008032)))" +15,0.008432622,"log(((exp(x0 / 2.1213145) + x1) * (x1 - 7.18395)) / (x1 * -0.08667538))" +17,0.008372297,"log(((exp(x0 / 2.1409464) + x1) * (x1 - 7.1519527)) / (x1 * -0.11701374)) - -0.31662118" +19,0.007871232,"log(((exp(x0 / 2.1483574) + x1) * (x1 - 6.9429812)) / (x1 * -0.091777064)) - (x1 * -0.03532463)" +21,0.0076451213,"log((x1 + exp((x0 + (x4 / 0.40371257)) * 0.08877082)) * ((x1 - 6.7461934) / ((x1 * -0.045270775) - 0.11868309)))" +23,0.0069678803,"log((exp(x4) + ((0.03762332 / (x1 / (x3 + 0.008522787))) + x1)) * (((x1 - 7.850765) / (x1 * -0.09465818)) - 0.9733152))" +26,0.006417365,"(log((x1 * (x1 - 7.2052436)) / (x1 * -0.12828495)) + (((0.38349107 - ((x1 - (x2 / x1)) * 0.66997904)) + x0) * 0.019127231)) - -0.5852694" diff --git a/trained_voltage_kernel.csv b/trained_voltage_kernel.csv new file mode 100644 index 0000000..6dee082 --- /dev/null +++ b/trained_voltage_kernel.csv @@ -0,0 +1,121 @@ +kernel_index,kernel_value +0,0.0016021910907593028 +1,-4.0884930635222905e-06 +2,-6.543550294005817e-05 +3,-5.923088719233176e-05 +4,-7.289373880506767e-05 +5,-8.566978668828535e-05 +6,-6.639146925721663e-05 +7,-7.271308853325891e-05 +8,-2.0417099533859173e-05 +9,-3.565276812846345e-05 +10,-5.6634704554986584e-05 +11,-0.00012460932046737367 +12,-6.813575431332786e-05 +13,-3.140492087849164e-05 +14,-4.1995314335183927e-05 +15,-5.0445133731274495e-06 +16,5.559875330028438e-06 +17,-7.727837365827751e-05 +18,-1.3623367097185207e-05 +19,1.6664994451469815e-05 +20,-4.81466990331218e-05 +21,3.556746102586017e-05 +22,-8.920613427976107e-06 +23,3.7297256579628e-06 +24,-4.545161971827369e-05 +25,-7.581274058660414e-07 +26,1.1903935603146709e-05 +27,1.653032244032839e-05 +28,-2.8504892257652224e-05 +29,1.0248960737772253e-05 +30,-3.8173641882253584e-09 +31,-1.2719382337264279e-05 +32,-3.8911651210043464e-05 +33,-7.578265629667978e-05 +34,-2.531017834132183e-05 +35,6.684916870507482e-06 +36,3.433759908018908e-05 +37,3.808656075112959e-05 +38,5.357069466017738e-06 +39,1.0124835073419377e-05 +40,-7.209351603727322e-07 +41,-2.942094697222857e-05 +42,7.333931645794848e-06 +43,2.8714010673019362e-06 +44,-2.176346298855163e-06 +45,1.7582566407552437e-05 +46,4.3726580659220965e-06 +47,6.513877590252509e-06 +48,-8.981849322069735e-06 +49,-1.7832531425315393e-05 +50,-1.2934148529870239e-05 +51,1.5072675248991178e-05 +52,5.158521666833451e-06 +53,-4.699881365021672e-06 +54,-4.997262505055129e-07 +55,-4.864605845080944e-06 +56,-7.222221700182308e-06 +57,1.019830299836899e-05 +58,3.5328606544084256e-07 +59,-7.953731373057138e-06 +60,-1.6673943958008404e-06 +61,1.1030948939611642e-05 +62,5.442141380454296e-06 +63,1.4176990671027835e-05 +64,-6.223835457063751e-06 +65,-1.501648761012208e-05 +66,-1.2729349183847797e-05 +67,1.1040455428251309e-05 +68,-4.212101793705198e-06 +69,8.753961825378478e-06 +70,9.71889860329778e-06 +71,-5.327615259186951e-06 +72,3.572159592490699e-05 +73,1.2015277847811226e-06 +74,-3.5531123599146294e-05 +75,9.58928571395551e-06 +76,-1.4597442556245373e-05 +77,3.1617588041629416e-06 +78,-3.494675308924198e-05 +79,-1.0827815130772643e-05 +80,-1.9700933189496676e-05 +81,-2.0053602129667174e-05 +82,-3.946768677275384e-05 +83,-1.155832402517126e-05 +84,-5.3056452195859455e-06 +85,6.764137042767126e-06 +86,2.9520591341283454e-05 +87,9.870032805621033e-06 +88,2.598541975095119e-05 +89,1.168893522617852e-06 +90,5.09963837554457e-06 +91,-6.201519135333833e-06 +92,1.075195591931887e-05 +93,-7.958231666394961e-06 +94,-1.566085373584499e-05 +95,-3.444274870237118e-06 +96,1.6484175637856208e-06 +97,-2.4696790526136187e-06 +98,-2.5718225865853134e-06 +99,-2.870814602155051e-06 +100,-1.2997139385802445e-05 +101,-6.140050599490011e-06 +102,-2.6482781122335405e-06 +103,1.201689609952287e-06 +104,1.5237792878081123e-05 +105,1.4758575169866904e-05 +106,1.6252014680216034e-05 +107,1.640190268335415e-05 +108,8.522252771235008e-06 +109,1.740461432691417e-05 +110,2.334211645866296e-05 +111,2.258734604837079e-05 +112,-5.011967716732702e-06 +113,-7.502420582653692e-06 +114,3.5450748168386586e-06 +115,4.045008641604303e-05 +116,-3.8241539372363195e-06 +117,-3.574013864701268e-06 +118,-1.6944039324449247e-05 +119,-0.00019550740220202864