modifying cost function

[Tine11L]

Hi everyone,

I'm working on a Lithium-ion battery model using UQpy for uncertainty quantification. I'm interested in implementing the Weighted Total Least Squares (WTLS) cost function for parameter estimation during the UQ process.

I have theese two functions, which I wish to improve the current loss/cost function, current one only looks at the voltage peaks, but with adding time you will get me better results. But for that I would need to somehow use the current real data (data_3) with the simulated data. Do you have any recomendations on how to do that?

-- I hope the questions was stated well, if not please don't hesitate for more details.
`
def run_uq(config_file_path, fig_save, C_file, sample_file, param_nb_down, param_nb_up):

"""
    Performs the uncertainty quantification (UQ) process based on the provided configuration file.

Args:
    config_file_path (str):     The path to the YAML configuration file.
    fig_save (str, optional):   Path to save a figure (not currently implemented). Defaults to None.
    sample_file (str):          Path to save the simulation results as a CSV file.
    param_nb_down (int):        Index of the first parameter to use (inclusive).
    param_nb_up (int):          Index of the last parameter to use (exclusive).

Returns:
    pd.DataFrame: A pandas DataFrame containing the simulation results for the specified parameters.
"""
global my_parameters_to_estimate 

# Load config
with open(config_file_path, 'r') as stream:
    config = yaml.safe_load(stream)

my_parameters_to_estimate = list(config['parameters'].items())[param_nb_down:param_nb_up]

# Generate data
var_names = [value['name'] for key, value in my_parameters_to_estimate if 'name' in value]

model = PythonModel(model_script = config['config']['model_script'],
                    model_object_name = config['config']['model_object_name'],
                    var_names = var_names)

param_true = np.array([value['reality'] for key, value in config['parameters'].items() if 'reality' in value]).reshape((1, -1))

h_func = RunModel(model = model)
h_func.run(samples = param_true)
data_clean = np.array(h_func.qoi_list[0])

# Add noise, use a RandomState for reproducible results
error_covariance = 0.0001
noise = Normal(loc = 0., scale = np.sqrt(error_covariance)).rvs(nsamples=len(data_clean), random_state=123).reshape((len(data_clean),))    

data_3 = data_clean + noise
print("data_3_1")
print(data_3)
print('Shape of data: {}'.format(data_3.shape))

# plt.plot(data_clean)
# plt.plot(data_3)
# plt.savefig(fig_save)

# Run the optimization
estimated_param_values = list(config['parameters'].values())[param_nb_down:param_nb_up]

p = []
for value in estimated_param_values:
    px = Uniform(loc=value['loc'], scale=value['scale'])
    #px = Lognormal(s=value['s'],loc=value['loc'], scale=value['scale'])
    p.append(px)

m = []
for value in estimated_param_values:
    mx = Normal(scale=value['marginal'])
    #mx = Lognormal(s=value['s'],loc=value['loc'], scale=value['scale'])
    m.append(mx)

seed = [value.get('init') for value in estimated_param_values if 'init' in value]

prior = JointIndependent(marginals = p)

# shouldn't the implemenentation of the new function be here?
inference_model = ComputationalModel(n_parameters = len(estimated_param_values),
                                     runmodel_object = h_func,
                                     error_covariance = error_covariance,
                                     prior = prior)

proposal = JointIndependent(m)

mh1 = MetropolisHastings(jump = config['config']['jump'],
                         burn_length = config['config']['burn_length'],
                         proposal = proposal,
                         seed = seed,
                         random_state = config['config']['random_state'])

# implement the new cost function, delete the "data = data_3" ??? 
bayes_estimator = BayesParameterEstimation(inference_model = inference_model,
                                       data = data_3,
                                       sampling_class = mh1,
                                       nsamples = config['config']['samples'])

print("data_3_2")
print(data_3)
s = bayes_estimator.sampler.samples

# Save the simulations
#column_names = [value['description'] for key, value in parameters_to_estimate if 'description' in value]
#print(column_names)
column_names = [value['description'] for key, value in list(config['parameters'].items())[param_nb_down:param_nb_up]if 'description' in value]

df = pd.DataFrame(s, columns = column_names)
print(df)
desired_dir = r"Sample_test"
sample_file_basename = os.path.join(desired_dir, C_file, sample_file)
df.to_csv(sample_file_basename, index = False)
print(sample_file_basename)
return s`

This is the second function
`
def Obtaining_parameters(theta, data_3, experiment = 1):

inpt = np.array(theta).reshape((-1,))   # set the input (theta) into an useable array
i = 0

# model = pybamm.lithium_ion.DFN()      # Doyle-Fuller-Newman model
# model = pybamm.lithium_ion.SPMe()     # Simple particle model with electrolyte dynamics
model = pybamm.lithium_ion.SPM()        # Simple particle model

parameter_values = pybamm.ParameterValues("Chen2020") # all the in-built parameter/constant sets

theta_0_value = get_parameter_value('theta_0', my_parameters_to_estimate)
if theta_0_value != 0:
   parameter_values["Negative electrode thickness [m]"] = inpt[i] * 1e-6
   print(f'neg elec thickness {inpt[i]}')
   i = i + 1

theta_1_value = get_parameter_value('theta_1', my_parameters_to_estimate)
if theta_1_value != 0:
    parameter_values["Positive electrode thickness [m]"] = inpt[i] * 1e-6
    print(f'pos elec thickness {inpt[i]}')
    i = i + 1

theta_2_value = get_parameter_value('theta_2', my_parameters_to_estimate)
if theta_2_value != 0:
    parameter_values["Negative particle radius [m]"] = inpt[i] * 1e-6
    print(f'neg particle radius {inpt[i]}')
    i = i + 1 

theta_3_value = get_parameter_value('theta_3', my_parameters_to_estimate)
if theta_3_value != 0:
    parameter_values["Positive particle radius [m]"] = inpt[i] * 1e-6
    print(f'pos particle radius {inpt[i]}')
    i = i + 1   

theta_4_value = get_parameter_value('theta_4', my_parameters_to_estimate)
if theta_4_value != 0:
    parameter_values["Separator thickness [m]"] = inpt[i] * 1e-6
    print(f'separator thickness {inpt[i]}')
    i = i + 1  

theta_5_value = get_parameter_value('theta_5', my_parameters_to_estimate)
if theta_5_value != 0:
    parameter_values["Negative electrode diffusivity [m2.s-1]"] = inpt[i] * 1e-14
    print(f'neg elec diffusivity {inpt[i]}')
    i = i + 1  

theta_6_value = get_parameter_value('theta_6', my_parameters_to_estimate)
if theta_6_value != 0:
    parameter_values["Positive electrode diffusivity [m2.s-1]"] = inpt[i] * 1e-15
    print(f'pos elec diffusivity {inpt[i]}')
    i = i + 1 

# theta_1_description = get_parameter_value('theta_1', my_parameters_to_estimate, key='description')

parameter_values.set_initial_stoichiometries(1) # depending on the current SOC
C = "1"

if experiment == 1:
    C = "1"

elif experiment == 0.5:
    C = "2"

elif experiment == 0.2:
    C = "5"

elif experiment == 0.1:
    C = "10"

elif experiment == 0.05:
    C = "20"

discharge_instructions = [
("Rest for 5 minutes", f"Discharge at C/{C} until 2.5 V"),  
]

experiment = pybamm.Experiment(
    [
        (
            discharge_instructions[0][1]
            
        )
    ]
)

sim = pybamm.Simulation(model, parameter_values = parameter_values, experiment = experiment)
sim.solve([0, 3600]) # time to solve the model?

solution = sim.solution # from the "array" of variables, you can access these
t = solution["Time [s]"].entries # entries - corresponded data
V = solution["Voltage [V]"].entries

return V`