import os
import sys
import time
import warnings
from typing import List, Dict, Any, Union
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import interpolate
import mph
from sko.GA import GA
from sko.PSO import PSO
from sko.DE import DE
warnings.filterwarnings('ignore') # Suppress unnecessary warnings
plt.ion() # Enable interactive mode for matplotlib
# model.study("std1").feature("time").set("tlist", "range(0,10,15000)");
#==================== Parameter Definitions ===================================
OP='DE' # Optimization method: GA, PSO, DE
C_Factor=3.6 # Coulomp to mAh transform Factor
U_column=2 # Column index of voltage data in xlsx file,array index starts from 0
Eneg_column=1 # Column index of Eneg data in xlsx file,array index starts from 0
C_rate=[0.33,0.5,1.0,2.0,3.0,3.5] # C-rates for test data
max_iter=20 # Optimization maximum iterations
size_pop=40 # Optimization population size or particle numb
WEIGHT_FACTOR_Eneg = 0.5 # Weight factor for voltage error calculation,1 for only voltage error
WEIGHT_FACTOR_Cap = -0.3 # Weight factor for voltage error calculation,1 for only voltage error
MODEL_NAME ='opt_test1_noT.mph' # COMSOL model file name
Sim_Str=['Caps', 'Ecell', 'Eneg'] # Simulation variables: time, cell voltage, negative electrode potential
Time_STR="range(0,10,15000)" # Time list for simulation
Eneg_Bool=True # Boolean flag to include negative electrode potential
Temp_func_Bool=False # Boolean flag to use test Temp data
Current_func_Bool=False # Boolean flag to use test current data
Temp_func="int3" # Temperature function name in COMSOL model
Current_func="int4" # Current function name in COMSOL model
INPUT_DIR = './input' # Input data directory
OUTPUT_DIR = './output' # Output data directory
OCV_NAME ='CELL_OCV.xlsx' # OCV data header
OCV_FILE = os.path.join(INPUT_DIR, OCV_NAME) # OCV data file
MODEL_PATH = os.path.join(os.getcwd(), MODEL_NAME) # COMSOL model file path
# Parameter bounds for genetic algorithm optimization
# Format: 'parameter_name': (lower_bound, upper_bound)
PARAM_BOUNDS = {
'brugl_pos': (1.5, 2.5), 'brugl_neg': (1.5, 2.5),
'k_pos1': (1e-11, 1e-9),
'k_neg1': (1e-10, 1e-8),
'Ds_neg1': (1e-14, 5e-12),
'Ds_pos1': (1e-16, 1e-13),
'Dl_1': (1e-11, 5e-9),
# 'Ea_Ds_pos': (10000, 100000), 'Ea_Ds_neg': (10000, 100000), 'Ea_Dl': (10000, 100000), 'Ea_k_pos': (10000, 100000), 'Ea_k_neg': (10000, 100000),
}
PARAM_UNITS = { # Define parameter units for proper COMSOL formatting
'brugl_pos': '', 'brugl_neg': '',
'k_pos1': '[m/s]',
'k_neg1': '[m/s]',
'Ds_neg1': '[m^2/s]',
'Ds_pos1': '[m^2/s]',
'Dl_1': '[m^2/s]',
# 'Ea_Ds_pos': '[J/mol]', 'Ea_Ds_neg': '[J/mol]', 'Ea_Dl': '[J/mol]', 'Ea_k_pos': '[J/mol]', 'Ea_k_neg': '[J/mol]',
}
# Common bounds for all algorithms
lb = [bounds[0] for bounds in PARAM_BOUNDS.values()]
ub = [bounds[1] for bounds in PARAM_BOUNDS.values()]
n_dim = len(PARAM_BOUNDS)
#========================== Algorithm parameter configurations ================
GA_PARAMS = {
'func': None, # Will be set dynamically
'n_dim': n_dim,
'size_pop': size_pop,
'max_iter': max_iter,
'prob_mut': 0.01,
'lb': lb,
'ub': ub
}
PSO_PARAMS = {
'func': None, # Will be set dynamically
'n_dim': n_dim,
'pop': size_pop, # Size of population/particles (consistent with GA parameter naming)
'max_iter': max_iter,
'w': 0.8, # Inertia weight
'c1': 0.5, # Cognitive parameter
'c2': 0.5, # Social parameter
'lb': lb,
'ub': ub
}
DE_PARAMS = {
'func': None, # Will be set dynamically
'n_dim': n_dim,
'size_pop': size_pop,
'max_iter': max_iter,
'lb': lb,
'ub': ub,
'F': 0.5, # Mutation factor
'prob_mut': 0.3 # Mutation probability (DE uses 'prob_mut' instead of 'CR')
}
#========================== Global Variables ==================================
model = None # COMSOL model instance
client = None # COMSOL client instance
pymodel = None # PyCOMSOL model interface
data_input = [] # List to store input datasets
SOC_rate = [] # List to store SOC values for each dataset
in_optimization_loop = True # Flag to track if we're in the main optimization loop
#========================== Helper Functions ==================================
def load_input_data() -> List[Dict[str, Any]]:
"""Load input data files from the INPUT_DIR.
Returns:
List[Dict[str, Any]]: List of dictionaries containing filename and data for each Excel file
"""
global data_input, SOC_rate
# Get all Excel files in the input directory, excluding OCV files
excel_files = [f for f in os.listdir(INPUT_DIR) if f.endswith('.xlsx') and 'OCV' not in f.upper()]
data_input = [] # Reset data_input list
for filename in excel_files: # Load each Excel file
file_path = os.path.join(INPUT_DIR, filename)
try:
df = pd.read_excel(file_path) # Read Excel file
data_array = df.values # Convert to numpy array
data_input.append({ # Add to data_input list with metadata
'filename': filename,
'data': data_array
})
print(f"Successfully loaded file: {filename}, data shape: {data_array.shape}")
except Exception as e:
print(f"Failed to load file {filename}: {str(e)}")
if data_input: # Calculate initial SOC values
calculate_initial_soc()
return data_input
def calculate_initial_soc() -> None:
"""Calculate initial State of Charge (SOC) for each dataset using OCV lookup."""
global SOC_rate
SOC_rate = [] # Reset SOC_rate list
try: # Read OCV data from Excel file
ocv_df = pd.read_excel(OCV_FILE)
ocv_data = ocv_df.values # Convert to numpy array
unique_ocv = {} # Store unique OCV data points to avoid duplicates
for row in ocv_data:
if row[1] not in unique_ocv:
unique_ocv[row[1]] = row[0]
sorted_ocv = np.array([[v, k] for k, v in unique_ocv.items()]) # Create and sort OCV data for interpolation
sorted_ocv = sorted_ocv[np.argsort(sorted_ocv[:, 1])] # Sort by voltage
for data_item in data_input:
if data_item['data'].shape[1] >= 3: # Ensure data has enough columns
initial_voltage = data_item['data'][0, int(U_column)] # Column index int(U_column) is Ecell (Cell Voltage)
f = interpolate.interp1d(sorted_ocv[:, 1], sorted_ocv[:, 0], # Create interpolation function for OCV-SOC relationship
bounds_error=False,fill_value="extrapolate")
soc = float(f(initial_voltage)) # Convert to scalar
SOC_rate.append(soc)
print(f"Initial SOC for file {data_item['filename']}: {soc}")
except Exception as e:
print(f"Failed to calculate initial SOC: {str(e)}")
SOC_rate = [0.01] * len(data_input) # Set default SOC value for each dataset
print(f"Using default SOC value of 0.01 for all datasets")
def initialize_comsol_model() -> bool:
"""Initialize COMSOL model by starting client and loading the model file.
Returns:
bool: True if initialization successful, False otherwise
"""
global client, pymodel, model
try:
print("Starting COMSOL client...")
client = mph.start()
print("Loading COMSOL model...")
# Check if model file exists before attempting to load
if os.path.exists(MODEL_PATH):
pymodel = client.load(MODEL_PATH)
print("Creating Java Object...")
model = pymodel.java
print("COMSOL model loaded successfully")
return True
else:
print(f"Model file not found: {MODEL_PATH}")
# Clean up client if model loading fails
client = None
return False
except Exception as e:
print(f"Failed to initialize COMSOL model: {str(e)}")
client = None # Ensure resources are released on failure
pymodel = None
model = None
return False
def set_model_parameters(params: Union[List[float], np.ndarray]) -> bool:
"""Set parameters in the COMSOL model with proper unit handling.
Args:
params: List or array of parameter values in the same order as PARAM_BOUNDS keys
Returns:
bool: True if parameters were set successfully, False otherwise
"""
global model
try:
if model is None: # Ensure model is initialized
print("Cannot set parameters: COMSOL model not initialized")
return False
param_names = list(PARAM_BOUNDS.keys()) # Get parameter names from PARAM_BOUNDS
if len(params) != len(param_names): # Verify parameter count matches expected number
print(f"Parameter count mismatch: expected {len(param_names)}, got {len(params)}")
return False
for i, param_name in enumerate(param_names): # Set each parameter in the model with appropriate unit
param_value = float(params[i]) # Ensure numeric value
unit = PARAM_UNITS[param_name] # Get unit from dictionary
if unit: # Set parameter with or without unit as needed
model.param().set(param_name, f'{param_value}{unit}')
else:
model.param().set(param_name, f'{param_value}')
return True
except Exception as e:
print(f"Failed to set model parameters: {str(e)}")
return False
def calculate_rmse(params: Union[List[float], np.ndarray]) -> float:
"""Calculate Root Mean Square Error (RMSE) between experimental and simulation data.
Args:
params: List or array of model parameters to evaluate
Returns:
float: Average RMSE value across all datasets in mV. Returns 666 (large value)
in case of errors to indicate optimization failure.
"""
global model, data_input, SOC_rate, in_optimization_loop
try:
if not set_model_parameters(params): # Set model parameters and check if successful
print("Failed to set model parameters in calculate_rmse") # When Set Success , Return Ture,Not True=Faluse, Not In Loop
return 666 # Return large value to indicate failure
RMSE_values = [] # Initialize storage for RMSE values
n = len(data_input)
for i in range(n): # Run Process each dataset
try:
# model.param().set('SOC', f'{SOC_rate[i]}') # Set initial SOC value for current dataset
data_array = data_input[i]['data'] # Get input data for current dataset
if data_array.shape[1] >= 3: # Ensure input data has sufficient columns
base_name = os.path.splitext(data_input[i]['filename'])[0] # Get base name directly from data_input
if Temp_func_Bool:
temp_file_path = os.path.join(INPUT_DIR, f"{base_name}_Temp.txt") # Dynamic temp file path
model.component("comp1").func(Temp_func).set("filename", temp_file_path)
if Current_func_Bool:
current_file_path = os.path.join(INPUT_DIR, f"{base_name}_Current.txt") # Dynamic current file path
model.component("comp1").func(Current_func).set("filename", current_file_path)
else: # else set C parameter
model.param().set("C", str(C_rate[i])) # c_value = base_name.replace("C", "")
print(f"Set C parameter to {C_rate[i]} for dataset{data_input[i]['filename']}")
model.study("std1").feature("time").set("tlist", Time_STR); # Set time list for current dataset
model.sol("sol1").runAll() # Run COMSOL simulation
# try: # Extract simulation data from COMSOL
sim_data = pymodel.evaluate(Sim_Str) # Evaluate simulation variables
sim_Time = sim_data[:, 0]/C_Factor # Time in seconds or Capacity in mAh
sim_Ecell = sim_data[:, 1] # Cell voltage data
sim_Eneg_10 = sim_data[:, 2] # Negative electrode potential data
# except Exception as inner_e:
# print(f"Failed to use evaluate method for dataset {data_input[i]['filename']}: {str(inner_e)}")
# continue
else:
print(f"Dataset {data_input[i]['filename']} does not have sufficient columns for processing")
continue
exp_Time = data_array[:, 0] # Experimental time data
exp_Ecell = data_array[:, int(U_column)] # Experimental cell voltage
if Eneg_Bool:
exp_Eneg_10 = data_array[:, int(Eneg_column)] # Experimental negative electrode potential
t_min = min(max(exp_Time), max(sim_Time)) # Determine interpolation time range
t_min_array = np.linspace(0, t_min, 300) # Create uniform time points for interpolation
# Create interpolation functions for simulation and experimental data
f_sim_Ecell = interpolate.interp1d(sim_Time, sim_Ecell, bounds_error=False, fill_value="extrapolate")
f_exp_Ecell = interpolate.interp1d(exp_Time, exp_Ecell, bounds_error=False, fill_value="extrapolate")
f_sim_Eneg_10 = interpolate.interp1d(sim_Time, sim_Eneg_10, bounds_error=False, fill_value="extrapolate")
# Interpolate data to match time points
sim_interp_Ecell = f_sim_Ecell(t_min_array)
exp_interp_Ecell = f_exp_Ecell(t_min_array)
sim_interp_Eneg_10 = f_sim_Eneg_10(t_min_array)
# Calculate RMSE values in millivolts
diff_Ecell = (exp_interp_Ecell - sim_interp_Ecell) * 1000 # Convert to mV
RMSE_Ecell = np.sqrt(np.mean(diff_Ecell**2)) # Voltage RMSE
if Eneg_Bool:
f_exp_Eneg_10 = interpolate.interp1d(exp_Time, exp_Eneg_10, bounds_error=False, fill_value="extrapolate")
exp_interp_Eneg_10 = f_exp_Eneg_10(t_min_array)
diff_Eneg_10 = (exp_interp_Eneg_10 - sim_interp_Eneg_10) * 1000 # Convert to mV
RMSE_Eneg_10 = np.sqrt(np.mean(diff_Eneg_10**2)) # Negative potential RMSE
else:
RMSE_Eneg_10 = 0
# Calculate weighted average of RMSE values
if WEIGHT_FACTOR_Cap > 0:
RMSE_Cap=np.sqrt((exp_Time[-1]-sim_Time[-1])**2) # Capacity RMSE
RMSE_val = RMSE_Ecell * (1-WEIGHT_FACTOR_Eneg-WEIGHT_FACTOR_Cap) + RMSE_Eneg_10 *WEIGHT_FACTOR_Eneg+RMSE_Cap*WEIGHT_FACTOR_Cap
RMSE_values.append(RMSE_val)
else:
RMSE_val = RMSE_Ecell * (1-WEIGHT_FACTOR_Eneg) + RMSE_Eneg_10 *WEIGHT_FACTOR_Eneg
RMSE_values.append(RMSE_val)
rmse_val_scalar = float(RMSE_val[0]) if isinstance(RMSE_val, np.ndarray) else float(RMSE_val)
print(f"RMSE for dataset {data_input[i]['filename']} : {rmse_val_scalar:.4f} mV") # Convert to scalar if necessary for proper printing
except Exception as e:
print(f"Error calculating RMSE for dataset {data_input[i]['filename']}: {str(e)}")
# RMSE_values.append(666) # Append large value to indicate dataset failure
print(f"Error calculating RMSE for dataset {data_input[i]['filename']}:Re-Start Calculattion")
return 666
if RMSE_values: # Calculate and return average RMSE
avg_rmse = np.mean(RMSE_values)
avg_rmse_scalar = float(avg_rmse[0]) if isinstance(avg_rmse, np.ndarray) else float(avg_rmse)
print(f"Average RMSE for dataset: {avg_rmse_scalar:.4f} mV")
return avg_rmse
else:
print("No valid RMSE values calculated")
return 666
except Exception as e:
print(f"Error occurred while calculating RMSE: {str(e)}")
return 666
def run_optimization(func, op_type='GA') -> tuple:
"""Run optimization using the specified algorithm.
Args:
func: The objective function to minimize
op_type: Algorithm type ('GA', 'PSO', 'DE')
Returns:
tuple: (best_params, best_rmse, algorithm_instance)
"""
global in_optimization_loop
# Set the objective function in the appropriate parameter dictionary
if op_type == 'GA':
params = GA_PARAMS.copy()
params['func'] = func
optimizer = GA(**params)
print(f"Running Genetic Algorithm (GA) optimization...")
print(f"Population size: {params['size_pop']}")
print(f"Maximum iterations: {params['max_iter']}")
print(f"Mutation probability: {params['prob_mut']}")
elif op_type == 'PSO':
params = PSO_PARAMS.copy()
params['func'] = func
optimizer = PSO(**params)
print(f"Running Particle Swarm Optimization (PSO)...")
print(f"Swarm size: {params['pop']}")
print(f"Maximum iterations: {params['max_iter']}")
print(f"Inertia weight (w): {params['w']}")
print(f"Cognitive parameter (c1): {params['c1']}")
print(f"Social parameter (c2): {params['c2']}")
elif op_type == 'DE':
params = DE_PARAMS.copy()
params['func'] = func
optimizer = DE(**params)
print(f"Running Differential Evolution (DE)...")
print(f"Population size: {params['size_pop']}")
print(f"Maximum iterations: {params['max_iter']}")
print(f"Mutation factor (F): {params['F']}")
print(f"Mutation probability: {params['prob_mut']}")
else:
raise ValueError(f"Unsupported optimization algorithm: {op_type}")
print(f"Number of parameters to optimize: {params['n_dim']}")
print("="*60)
try:
best_params, best_rmse = optimizer.run()
return best_params, best_rmse, optimizer
except KeyboardInterrupt:
print("\nWARNING: Optimization interrupted by user")
return None, None, optimizer
except Exception as e:
print(f"ERROR: Optimization failed: {str(e)}")
return None, None, None
def plot_optimization_results(op_instance: Any, op_type='GA') -> None:
"""Plot and save optimization convergence curve.
Args:
op_instance: The optimization instance containing generation data
op_type: Algorithm type ('GA', 'PSO', 'DE')
"""
plt.figure(figsize=(10, 6))
# Get convergence data based on algorithm type
try:
if op_type == 'GA':
convergence_data = op_instance.generation_best_Y
elif op_type == 'PSO':
convergence_data = op_instance.gbest_y_hist
elif op_type == 'DE':
convergence_data = op_instance.generation_best_Y
else:
# Default fallback for any other algorithm
convergence_data = getattr(op_instance, 'generation_best_Y', [])
if not convergence_data:
convergence_data = getattr(op_instance, 'gbest_y_history', [])
if convergence_data:
plt.plot(convergence_data)
plt.xlabel('Number of Generations' if 'generation' in str(type(convergence_data)) else 'Iterations')
plt.ylabel('RMSE Value (mV)')
plt.title(f'{op_type} Optimization Convergence Curve')
plt.grid(True)
# Save with algorithm-specific filename
filename = f'{op_type.lower()}_convergence_curve.png'
plt.savefig(os.path.join(OUTPUT_DIR, filename), dpi=300)
plt.show(block=False)
print(f"Convergence plot saved as: {filename}")
else:
print(f"Warning: No convergence data found for {op_type}")
except Exception as e:
print(f"Error plotting convergence curve for {op_type}: {str(e)}")
def save_results(op_instance, op_type='GA', best_params=None, best_rmse=None) -> None:
"""Save optimization results to output files.
Args:
op_instance: The optimization instance containing results
op_type: Algorithm type ('GA', 'PSO', 'DE')
best_params: Pre-computed best parameters (optional)
best_rmse: Pre-computed best RMSE (optional)
"""
try:
# Try to get best parameters and RMSE
if best_params is None or best_rmse is None:
if op_type == 'GA':
if best_params is None and hasattr(op_instance, 'generation_best_X'):
best_params = op_instance.generation_best_X[-1]
if best_rmse is None and hasattr(op_instance, 'generation_best_Y'):
best_rmse = op_instance.generation_best_Y[-1]
elif op_type == 'PSO':
if best_params is None and hasattr(op_instance, 'gbest_x'):
best_params = op_instance.gbest_x
if best_rmse is None and hasattr(op_instance, 'gbest_y'):
best_rmse = op_instance.gbest_y
elif op_type == 'DE':
if best_params is None and hasattr(op_instance, 'generation_best_X'):
best_params = op_instance.generation_best_X[-1]
if best_rmse is None and hasattr(op_instance, 'generation_best_Y'):
best_rmse = op_instance.generation_best_Y[-1]
# Ensure we have valid results
if best_params is None or best_rmse is None:
raise ValueError(f"Could not retrieve optimization results from {op_type} instance")
# Save results with algorithm-specific filename
output_file_path = os.path.join(OUTPUT_DIR, f'python_{op_type.lower()}_results.txt')
with open(output_file_path, 'w') as f:
f.write(f"Optimization Algorithm: {op_type}\n")
f.write(f"Best RMSE value: {float(best_rmse):.6f} mV\n\n")
f.write("Best parameter values:\n")
for i, param_name in enumerate(list(PARAM_BOUNDS.keys())):
param_val = float(best_params[i])
f.write(f"{param_name}: {param_val:.6e}\n")
print(f"Optimization results saved to: {output_file_path}")
except Exception as e:
print(f"Error saving optimization results for {op_type}: {str(e)}")
def generate_comparison_plots_with_best_params(best_params_op):
final_rmse = None
try:
final_rmse = calculate_rmse(best_params_op)
except Exception as rse:
print(f"WARNING: Failed to calculate final RMSE: {str(rse)}")
RMSE_values = []
for i in range(len(data_input)):
try:
# model.param().set('SOC', f'{SOC_rate[i]}') # Set initial SOC value for current dataset
data_array = data_input[i]['data'] # Get input data for current dataset
if data_array.shape[1] >= 3: # Ensure input data has sufficient columns
base_name = os.path.splitext(data_input[i]['filename'])[0] # Get base name directly from data_input
if Temp_func_Bool:
temp_file_path = os.path.join(INPUT_DIR, f"{base_name}_Temp.txt") # Dynamic temp file path
model.component("comp1").func(Temp_func).set("filename", temp_file_path)
if Current_func_Bool:
current_file_path = os.path.join(INPUT_DIR, f"{base_name}_Current.txt") # Dynamic current file path
model.component("comp1").func(Current_func).set("filename", current_file_path)
else: # else set C parameter
# c_value = base_name.replace("C", "")
model.param().set("C", str(C_rate[i]))
print(f"Set C parameter to {C_rate[i]} for dataset{data_input[i]['filename']}")
model.study("std1").feature("time").set("tlist", Time_STR); # Set time list for current dataset
model.sol("sol1").runAll() # Run COMSOL simulation
# try: # Extract simulation data from COMSOL
sim_data = pymodel.evaluate(Sim_Str) # Evaluate simulation variables
sim_Time = sim_data[:, 0]/C_Factor # Time in seconds
sim_Ecell = sim_data[:, 1] # Cell voltage data
sim_Eneg_10 = sim_data[:, 2] # Negative electrode potential data
# except Exception as inner_e:
# print(f"Failed to use evaluate method for dataset {C_rate[i]}.xlsx: {str(inner_e)}")
# continue
else:
print(f"Dataset {C_rate[i]}.xlsx does not have sufficient columns for processing")
continue
exp_Time = data_array[:, 0] # Experimental time data
exp_Ecell = data_array[:, int(U_column)] # Experimental cell voltage
if Eneg_Bool:
exp_Eneg_10 = data_array[:, int(Eneg_column)] # Experimental negative electrode potential
t_min = min(max(exp_Time), max(sim_Time)) # Determine interpolation time range
t_min_array = np.linspace(0, t_min, 300) # Create uniform time points for interpolation
# Create interpolation functions for simulation and experimental data
f_sim_Ecell = interpolate.interp1d(sim_Time, sim_Ecell, bounds_error=False, fill_value="extrapolate")
f_exp_Ecell = interpolate.interp1d(exp_Time, exp_Ecell, bounds_error=False, fill_value="extrapolate")
f_sim_Eneg_10 = interpolate.interp1d(sim_Time, sim_Eneg_10, bounds_error=False, fill_value="extrapolate")
# Interpolate data to match time points
sim_interp_Ecell = f_sim_Ecell(t_min_array)
exp_interp_Ecell = f_exp_Ecell(t_min_array)
sim_interp_Eneg_10 = f_sim_Eneg_10(t_min_array)
# Calculate RMSE values in millivolts
diff_Ecell = (exp_interp_Ecell - sim_interp_Ecell) * 1000 # Convert to mV
RMSE_Ecell = np.sqrt(np.mean(diff_Ecell**2)) # Voltage RMSE
if Eneg_Bool:
f_exp_Eneg_10 = interpolate.interp1d(exp_Time, exp_Eneg_10, bounds_error=False, fill_value="extrapolate")
exp_interp_Eneg_10 = f_exp_Eneg_10(t_min_array)
diff_Eneg_10 = (exp_interp_Eneg_10 - sim_interp_Eneg_10) * 1000 # Convert to mV
RMSE_Eneg_10 = np.sqrt(np.mean(diff_Eneg_10**2)) # Negative potential RMSE
else:
RMSE_Eneg_10 = 0
# Calculate weighted average of RMSE values
if WEIGHT_FACTOR_Cap > 0:
RMSE_Cap=np.sqrt((exp_Time[-1]-sim_Time[-1])**2) # Capacity RMSE
RMSE_val = RMSE_Ecell * (1-WEIGHT_FACTOR_Eneg-WEIGHT_FACTOR_Cap) + RMSE_Eneg_10 *WEIGHT_FACTOR_Eneg+RMSE_Cap*WEIGHT_FACTOR_Cap
RMSE_values.append(RMSE_val)
else:
RMSE_val = RMSE_Ecell * (1-WEIGHT_FACTOR_Eneg) + RMSE_Eneg_10 *WEIGHT_FACTOR_Eneg
RMSE_values.append(RMSE_val)
rmse_val_scalar = float(RMSE_val[0]) if isinstance(RMSE_val, np.ndarray) else float(RMSE_val)
print(f"RMSE for dataset {data_input[i]['filename']}: {rmse_val_scalar:.4f} mV") # Convert to scalar if necessary for proper printing
try: # Create and save plot
plt.figure(figsize=(18, 5))
# Cell voltage comparison plot
plt.subplot(1, 3, 1)
plt.plot(t_min_array, sim_interp_Ecell, 'r-', label='Simulation')
plt.plot(t_min_array, exp_interp_Ecell, 'bo', label='Experiment')
plt.xlabel('Time (s)')
plt.ylabel('Cell Voltage (V)')
plt.title('Cell Voltage: Simulation vs Experiment')
plt.legend()
plt.grid(True)
# Negative potential comparison plot
plt.subplot(1, 3, 2)
plt.plot(t_min_array, sim_interp_Eneg_10*1000, 'r-', label='Simulation')
if Eneg_Bool:
plt.plot(t_min_array, exp_interp_Eneg_10*1000, 'bo', label='Experiment')
plt.xlabel('Time (s)')
plt.ylabel('Negative Potential (mV)')
plt.title('Negative Potential: Simulation vs Experiment')
plt.legend()
plt.grid(True)
# Error comparison plot
plt.subplot(1, 3, 3)
plt.plot(t_min_array, diff_Ecell, 'g-', label='Cell Voltage Error')
if Eneg_Bool:
plt.plot(t_min_array, diff_Eneg_10, 'm-', label='Negative Potential Error')
plt.xlabel('Time (s)')
plt.ylabel('Error (mV)')
# Add RMSE to title if available
rmse_text = "N/A" if rmse_val_scalar is None else f"{float(rmse_val_scalar):.2f}"
plt.title(f'Errors (RMSE: {rmse_text} mV)')
plt.legend()
plt.grid(True)
plt.tight_layout()
# Save plot with standardized path
plot_filename = f'{OP}_op_result_comp_dataset_{data_input[i]["filename"]}.png'
plot_path = os.path.join(OUTPUT_DIR, plot_filename)
plt.savefig(plot_path, dpi=300)
plt.show(block=False)
print(f"Comparison plot saved as: {plot_path}")
except Exception as plot_e:
print(f"ERROR: Failed to generate plots for dataset {data_input[i]['filename']}: {str(plot_e)}")
except Exception as ei:
print(f"ERROR: Failed to process dataset {data_input[i]['filename']}: {str(ei)}")
if RMSE_values: # Calculate Optimal Average RMSE
avg_rmse = np.mean(RMSE_values)
avg_rmse_scalar = float(avg_rmse[0]) if isinstance(avg_rmse, np.ndarray) else float(avg_rmse)
print(f"Optimal Average RMSE: {avg_rmse_scalar:.4f} mV")
def main() -> int:
"""Main function to execute parameter optimization with selected algorithm.
Returns:
int: Exit code (0 for success, non-zero for errors)
"""
global client, pymodel, model, OP
start_time = None
exit_code = 0 # Default to success
# Validate optimization algorithm type
supported_algorithms = ['DE', 'GA', 'PSO']
if OP not in supported_algorithms:
print(f"ERROR: Unsupported optimization algorithm: {OP}")
print(f"Supported algorithms: {', '.join(supported_algorithms)}")
return 1
try:
start_time = time.time() # Record start time
print(f"===== Python {OP} Parameter Optimization Started =====")
# Step I: Load input data with error handling
print("\nI. Loading input data...")
try:
data_result = load_input_data()
if not data_result and not data_input: # Check both function return and global variable
print("ERROR: No valid input data found")
return 1
except Exception as e:
print(f"ERROR: Failed to load input data: {str(e)}")
return 1
# Step II: Initialize COMSOL model with error handling
print("\nII. Initializing COMSOL model...")
if not initialize_comsol_model():
print("ERROR: COMSOL model initialization failed")
return 2
# Step III & IV: Configure and run selected optimization algorithm
print(f"\nIII. Starting {OP} optimization...")
best_params, best_rmse, optimizer = run_optimization(calculate_rmse, OP)
if optimizer is None:
print(f"ERROR: Failed to initialize {OP} optimizer")
exit_code = 3
# Step V: Display optimization results with error handling
if best_params is not None and best_rmse is not None:
print("\nV. Optimization results:")
try:
best_rmse_val = float(best_rmse[0]) if isinstance(best_rmse, np.ndarray) else float(best_rmse)
print(f"Best RMSE value: {best_rmse_val:.6f} mV")
print("\nBest parameter values:")
for i, param_name in enumerate(list(PARAM_BOUNDS.keys())):
try:
param_val = float(best_params[i])
print(f"{param_name}: {param_val:.6e}")
except IndexError:
print(f"WARNING: Missing value for parameter {param_name}")
except Exception as e:
print(f"WARNING: Failed to process parameter {param_name}: {str(e)}")
except Exception as e:
print(f"ERROR: Failed to display optimization results: {str(e)}")
# Step VI: Plot and save results with error handling
if optimizer is not None:
print("\nVI. Plotting and saving results...")
try:
plot_optimization_results(optimizer, OP)
except Exception as e:
print(f"ERROR: Failed to plot optimization results: {str(e)}")
try:
save_results(optimizer, OP, best_params, best_rmse)
except Exception as e:
print(f"ERROR: Failed to save optimization results: {str(e)}")
# Step VII: Set optimal parameters and save model with error handling
if best_params is not None and pymodel is not None:
print("\nVII. Setting optimal parameters and saving model...")
try:
if set_model_parameters(best_params):
original_model_name = MODEL_PATH.split('.')[0] if '.' in MODEL_PATH else MODEL_PATH
# Add algorithm type to the optimized model filename
optimized_model_path = os.path.join(os.path.dirname(MODEL_PATH), f"{original_model_name}_Op_{OP}.mph")
pymodel.save(optimized_model_path)
print(f"Optimized model saved as: {optimized_model_path}")
else:
print("WARNING: Failed to set optimal parameters in model")
except Exception as e:
print(f"ERROR: Failed to save optimized model: {str(e)}")
# Step VIII: Generate comparison plots with optimal parameters
if best_params is not None and data_input:
print("\nVIII. Running simulation with optimal parameters...")
generate_comparison_plots_with_best_params(best_params)
# Calculate and display total execution time
if start_time:
end_time = time.time()
print(f"\nTotal optimization time: {(end_time - start_time):.2f} seconds")
if exit_code == 0:
print(f"\n===== Python {OP} Parameter Optimization Completed Successfully =====")
else:
print(f"\n===== Python {OP} Parameter Optimization Completed with Errors =====")
return exit_code
except KeyboardInterrupt:
print("\nERROR: Program interrupted by user")
return 130 # Standard exit code for keyboard interrupt
except Exception as e:
print(f"ERROR: Unexpected error during program execution: {str(e)}")
import traceback
print("Detailed error information:")
traceback.print_exc()
return 5
# Update the main execution line
if __name__ == "__main__":
sys.exit(main())
解析以上代码,给出每个函数的计算逻辑以及变量返回值的传递,重点拆接触如何对comsol中模型设置参数并且运行模型,获取comsol中计算结果与实测值计算RMS误差,调用算法寻优迭代等主要函数。给出具体解析说明,变量传递。
最新发布