Perovskite light-intensity dependant JV fits with SIMsalabim (real data)

This notebook is a demonstration of how to fit light-intensity dependent JV curves with drift-diffusion models using the SIMsalabim package.

# Import necessary libraries
import warnings, os, sys, shutil
# remove warnings from the output
os.environ["PYTHONWARNINGS"] = "ignore"
warnings.filterwarnings(action='ignore', category=FutureWarning)
warnings.filterwarnings(action='ignore', category=UserWarning)
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import default_rng
from scipy import constants
import torch, copy, uuid
import pySIMsalabim as sim
from pySIMsalabim.experiments.JV_steady_state import *
import ax, logging
from ax.utils.notebook.plotting import init_notebook_plotting, render
init_notebook_plotting() # for Jupyter notebooks

try:
    from optimpv import *
    from optimpv.axBOtorch.axUtils import *
except Exception as e:
    sys.path.append('../') # add the path to the optimpv module
    from optimpv import *
    from optimpv.axBOtorch.axUtils import *

[INFO 08-13 18:14:07] ax.utils.notebook.plotting: Injecting Plotly library into cell. Do not overwrite or delete cell.
[INFO 08-13 18:14:07] ax.utils.notebook.plotting: Please see
    (https://ax.dev/tutorials/visualizations.html#Fix-for-plots-that-are-not-rendering)
    if visualizations are not rendering.

Get the experimental data

# Set the session path for the simulation and the input files
session_path = os.path.join(os.path.join(os.path.abspath('../'),'SIMsalabim','SimSS'))
input_path = os.path.join(os.path.join(os.path.join(os.path.abspath('../'),'Data','simsalabim_test_inputs','JVrealPerovskite')))
simulation_setup_filename = 'simulation_setup_realPerovskite.txt'
simulation_setup = os.path.join(session_path, simulation_setup_filename)
optical_files = ['nk_glass.txt','nk_ITO.txt','nk_PTAA.txt','nk_FACsPbIBr.txt','nk_C60_1.txt','nk_Au.txt']
# path to the layer files defined in the simulation_setup file
l1 = 'PTAA.txt'
l2 = 'perovskite.txt'
l3 = 'C60.txt'
l1 = os.path.join(input_path, l1)
l2 = os.path.join(input_path, l2)
l3 = os.path.join(input_path, l3)

# copy this files to session_path
force_copy = True
if not os.path.exists(session_path):
    os.makedirs(session_path)
for file in [l1,l2,l3,simulation_setup_filename]:
    file = os.path.join(input_path, os.path.basename(file))
    if force_copy or not os.path.exists(os.path.join(session_path, os.path.basename(file))):
        shutil.copyfile(file, os.path.join(session_path, os.path.basename(file)))
    else:
        print('File already exists: ',file)
# copy the optical files to the session path
for file in optical_files:
    file = os.path.join(input_path, os.path.basename(file))
    if force_copy or not os.path.exists(os.path.join(session_path, os.path.basename(file))):
        shutil.copyfile(file, os.path.join(session_path, os.path.basename(file)))


# Show the device structure
fig = sim.plot_band_diagram(simulation_setup, session_path)
dev_par, layers = sim.load_device_parameters(session_path, simulation_setup, run_mode = False)
SIMsalabim_params  = {}
for layer in layers:
    SIMsalabim_params[layer[1]] = sim.ReadParameterFile(os.path.join(session_path,layer[2]))
L_active_Layer = float(SIMsalabim_params['l2']['L'])

# Load the JV data
df = pd.read_csv(os.path.join(input_path, 'JV_realPerovskite.dat'), sep=' ')
X = df['Vext'].values
y = df['Jext'].values
Gfracs = None

# get 1sun Jsc as it will help narrow down the range for Gehp
Jsc_1sun = np.interp(0.0, X, y)
minJ = np.max(y) # max because we have positive current density in the file
q = constants.value(u'elementary charge')
G_ehp_calc = abs(Jsc_1sun/(q*L_active_Layer))
G_ehp_max = abs(minJ/(q*L_active_Layer))

# plot the data for each Gfrac
plt.figure(figsize=(10,10))
plt.plot(X, y, '-')
plt.grid()
plt.xlabel('Voltage [V]')
plt.ylabel('Current density [A m$^{-2}$]')
plt.show()

Define the parameters for the simulation

params = [] # list of parameters to be optimized

mun = FitParam(name = 'l2.mu_n', value = 8e-5, bounds = [1e-5,1e-2], values = None, start_value = None, log_scale = True, value_type = 'float', fscale = None, rescale = False, stepsize = None, display_name=r'$\mu_n$', unit='m$^2$ V$^{-1}$s$^{-1}$', axis_type = 'log', std = 0,encoding = None,force_log = False)
params.append(mun)

mup = FitParam(name = 'l2.mu_p', value = 8e-5, bounds = [1e-5,1e-2], values = None, start_value = None, log_scale = True, value_type = 'float', fscale = None, rescale = False, stepsize = None, display_name=r'$\mu_p$', unit='m$^2$ V$^{-1}$s$^{-1}$', axis_type = 'log', std = 0,encoding = None,force_log = False)
params.append(mup)

bulk_tr = FitParam(name = 'l2.N_t_bulk', value = 5.1e20, bounds = [1e19,1e22], values = None, start_value = None, log_scale = True, value_type = 'float', fscale = None, rescale = False, stepsize = None, display_name=r'$N_{T}$', unit='s', axis_type = 'log', std = 0,encoding = None,force_log = False)
params.append(bulk_tr)

Nions = FitParam(name = 'l2.N_ions', value = 4e22, bounds = [1e22,5e22], type='range', values = None, start_value = None, log_scale = True, value_type = 'float', fscale = None, rescale = False, stepsize = None, display_name=r'$N_{ions}$', unit='m$^{-3}$', axis_type = 'log', std = 0,encoding = None,force_log = False)
params.append(Nions)

N_t_int_l1 = FitParam(name = 'l1.N_t_int', value = 5e12, bounds = [5e11,5e13], values = None, start_value = None, log_scale = True, value_type = 'float', fscale = None, rescale = False, stepsize = None, display_name=r'$N_{T,int}^{HTL}$', unit='m$^{-2}$', axis_type = 'log', std = 0,encoding = None,force_log = False)
params.append(N_t_int_l1)

N_t_int_l2 = FitParam(name = 'l2.N_t_int', value = 4e12, bounds = [5e11,5e13], values = None, start_value = None, log_scale = True, value_type = 'float', fscale = None, rescale = False, stepsize = None, display_name=r'$N_{T,int}^{ETL}$', unit='m$^{-2}$', axis_type = 'log', std = 0,encoding = None,force_log = False)
params.append(N_t_int_l2)

R_series = FitParam(name = 'R_series', value = 1e-4, bounds = [1e-6,1e-2], log_scale = True, value_type = 'float', fscale = None, rescale = False,  display_name=r'$R_{series}$', unit=r'$\Omega$ m$^2$', axis_type = 'log', force_log = False)
params.append(R_series)

R_shunt = FitParam(name = 'R_shunt', value = 1e1, bounds = [1e-2,1e2], log_scale = True, value_type = 'float', fscale = None, rescale = False,  display_name=r'$R_{shunt}$', unit=r'$\Omega$ m$^2$', axis_type = 'log', force_log = False)
params.append(R_shunt)

# save the original parameters for later
params_orig = copy.deepcopy(params)

Run the optimization

# Define the Agent and the target metric/loss function
from optimpv.DDfits.JVAgent import JVAgent
metric = 'nrmse' # can be 'nrmse', 'mse', 'mae'
loss = 'linear' # can be 'linear', 'huber', 'soft_l1'

jv = JVAgent(params, X, y, session_path, simulation_setup, parallel = False, max_jobs = 1, metric = metric, loss = loss)

from optimpv.axBOtorch.axBOtorchOptimizer import axBOtorchOptimizer
from botorch.acquisition.logei import qLogNoisyExpectedImprovement
from ax.adapter.transforms.standardize_y import StandardizeY
from ax.adapter.transforms.unit_x import UnitX
from ax.adapter.transforms.remove_fixed import RemoveFixed
from ax.adapter.transforms.log import Log
from ax.generators.torch.botorch_modular.utils import ModelConfig
from ax.generators.torch.botorch_modular.surrogate import SurrogateSpec
from gpytorch.kernels import MaternKernel
from gpytorch.kernels import ScaleKernel
from botorch.models import SingleTaskGP

model_gen_kwargs_list = None
parameter_constraints = None

model_kwargs_list = [{},{"torch_device":torch.device("cuda" if torch.cuda.is_available() else "cpu"),'botorch_acqf_class':qLogNoisyExpectedImprovement,'transforms':[RemoveFixed, Log,UnitX, StandardizeY],'surrogate_spec':SurrogateSpec(model_configs=[ModelConfig(botorch_model_class=SingleTaskGP,covar_module_class=ScaleKernel, covar_module_options={'base_kernel':MaternKernel(nu=2.5, ard_num_dims=len(params))})])}]

# Define the optimizer
optimizer = axBOtorchOptimizer(params = params, agents = jv, models = ['SOBOL','BOTORCH_MODULAR'],n_batches = [1,25], batch_size = [10,2], ax_client = None,  max_parallelism = -1, model_kwargs_list = model_kwargs_list, model_gen_kwargs_list = None, name = 'ax_opti',parallel=True)

# optimizer.optimize() # run the optimization with ax
optimizer.optimize_turbo() # run the optimization with turbo

[INFO 08-13 18:14:08] optimpv.axBOtorchOptimizer: Starting optimization with 26 batches and a total of 60 trials
[INFO 08-13 18:14:08] optimpv.axBOtorchOptimizer: Starting Sobol batch 1 with 10 trials
[INFO 08-13 18:14:10] optimpv.axBOtorchOptimizer: Finished Sobol
[INFO 08-13 18:14:13] optimpv.axBOtorchOptimizer: Finished Turbo batch 2 with 2 trials with current best value: 6.34e-02, TR length: 8.00e-01
[INFO 08-13 18:14:15] optimpv.axBOtorchOptimizer: Finished Turbo batch 3 with 2 trials with current best value: 4.00e-02, TR length: 8.00e-01
[INFO 08-13 18:14:16] optimpv.axBOtorchOptimizer: Finished Turbo batch 4 with 2 trials with current best value: 4.00e-02, TR length: 8.00e-01
[INFO 08-13 18:14:18] optimpv.axBOtorchOptimizer: Finished Turbo batch 5 with 2 trials with current best value: 4.00e-02, TR length: 8.00e-01
[INFO 08-13 18:14:19] optimpv.axBOtorchOptimizer: Finished Turbo batch 6 with 2 trials with current best value: 3.29e-02, TR length: 8.00e-01
[INFO 08-13 18:14:21] optimpv.axBOtorchOptimizer: Finished Turbo batch 7 with 2 trials with current best value: 3.29e-02, TR length: 8.00e-01
[INFO 08-13 18:14:22] optimpv.axBOtorchOptimizer: Finished Turbo batch 8 with 2 trials with current best value: 3.29e-02, TR length: 8.00e-01
[INFO 08-13 18:14:24] optimpv.axBOtorchOptimizer: Finished Turbo batch 9 with 2 trials with current best value: 3.18e-02, TR length: 8.00e-01
[INFO 08-13 18:14:25] optimpv.axBOtorchOptimizer: Finished Turbo batch 10 with 2 trials with current best value: 2.46e-02, TR length: 8.00e-01
[INFO 08-13 18:14:27] optimpv.axBOtorchOptimizer: Finished Turbo batch 11 with 2 trials with current best value: 2.46e-02, TR length: 8.00e-01
[INFO 08-13 18:14:28] optimpv.axBOtorchOptimizer: Finished Turbo batch 12 with 2 trials with current best value: 2.46e-02, TR length: 8.00e-01
[INFO 08-13 18:14:30] optimpv.axBOtorchOptimizer: Finished Turbo batch 13 with 2 trials with current best value: 2.46e-02, TR length: 8.00e-01
[INFO 08-13 18:14:31] optimpv.axBOtorchOptimizer: Finished Turbo batch 14 with 2 trials with current best value: 2.46e-02, TR length: 8.00e-01
[INFO 08-13 18:14:33] optimpv.axBOtorchOptimizer: Finished Turbo batch 15 with 2 trials with current best value: 2.46e-02, TR length: 8.00e-01
[INFO 08-13 18:14:34] optimpv.axBOtorchOptimizer: Finished Turbo batch 16 with 2 trials with current best value: 2.46e-02, TR length: 8.00e-01
[INFO 08-13 18:14:36] optimpv.axBOtorchOptimizer: Finished Turbo batch 17 with 2 trials with current best value: 1.33e-02, TR length: 8.00e-01
[INFO 08-13 18:14:37] optimpv.axBOtorchOptimizer: Finished Turbo batch 18 with 2 trials with current best value: 9.64e-03, TR length: 8.00e-01
[INFO 08-13 18:14:38] optimpv.axBOtorchOptimizer: Finished Turbo batch 19 with 2 trials with current best value: 9.64e-03, TR length: 8.00e-01
[INFO 08-13 18:14:40] optimpv.axBOtorchOptimizer: Finished Turbo batch 20 with 2 trials with current best value: 9.64e-03, TR length: 8.00e-01
[INFO 08-13 18:14:41] optimpv.axBOtorchOptimizer: Finished Turbo batch 21 with 2 trials with current best value: 9.64e-03, TR length: 8.00e-01
[INFO 08-13 18:14:42] optimpv.axBOtorchOptimizer: Finished Turbo batch 22 with 2 trials with current best value: 9.64e-03, TR length: 8.00e-01
[INFO 08-13 18:14:44] optimpv.axBOtorchOptimizer: Finished Turbo batch 23 with 2 trials with current best value: 9.64e-03, TR length: 8.00e-01
[INFO 08-13 18:14:46] optimpv.axBOtorchOptimizer: Finished Turbo batch 24 with 2 trials with current best value: 9.64e-03, TR length: 8.00e-01
[INFO 08-13 18:14:47] optimpv.axBOtorchOptimizer: Finished Turbo batch 25 with 2 trials with current best value: 9.64e-03, TR length: 8.00e-01
[INFO 08-13 18:14:49] optimpv.axBOtorchOptimizer: Finished Turbo batch 26 with 2 trials with current best value: 9.64e-03, TR length: 8.00e-01
[INFO 08-13 18:14:50] optimpv.axBOtorchOptimizer: Finished Turbo batch 27 with 2 trials with current best value: 9.64e-03, TR length: 8.00e-01
[INFO 08-13 18:14:50] optimpv.axBOtorchOptimizer: Turbo is terminated as the max number (60) of trials is reached

# get the best parameters and update the params list in the optimizer and the agent
ax_client = optimizer.ax_client # get the ax client
optimizer.update_params_with_best_balance() # update the params list in the optimizer with the best parameters
jv.params = optimizer.params # update the params list in the agent with the best parameters

# print the best parameters
print('Best parameters:')
for p in optimizer.params:
    print(p.name, 'fitted value:', p.value)

print('\nSimSS command line:')
print(jv.get_SIMsalabim_clean_cmd(jv.params)) # print the simss command line with the best parameters

Best parameters:
l2.mu_n fitted value: 0.007196673188405359
l2.mu_p fitted value: 0.0019405736081302167
l2.N_t_bulk fitted value: 6.354965308308602e+20
l2.N_ions fitted value: 4.874368835248908e+22
l1.N_t_int fitted value: 15796096401498.818
l2.N_t_int fitted value: 12437852501869.201
R_series fitted value: 0.0002189728346564998
R_shunt fitted value: 94.32845691557006

SimSS command line:
./simss -l2.mu_n 0.007196673188405359 -l2.mu_p 0.0019405736081302167 -l2.N_t_bulk 6.354965308308602e+20 -l2.N_anion 4.874368835248908e+22 -l2.N_cation 4.874368835248908e+22 -l1.N_t_int 15796096401498.818 -l2.N_t_int 12437852501869.201 -R_series 0.0002189728346564998 -R_shunt 94.32845691557006

# Plot optimization results
data = ax_client.summarize()
all_metrics = optimizer.all_metrics
plt.figure()
plt.plot(np.minimum.accumulate(data[all_metrics]), label="Best value seen so far")
plt.yscale("log")
plt.xlabel("Iteration")
plt.ylabel("Target metric: " + metric + " with " + loss + " loss")
plt.legend()
plt.title("Best value seen so far")

print("Best value seen so far is ", min(data[all_metrics[0]]), "at iteration ", int(data[all_metrics[0]].idxmin()))
plt.show()

Best value seen so far is  0.009638409702757542 at iteration  42

from ax.analysis.plotly.parallel_coordinates import ParallelCoordinatesPlot

analysis = ParallelCoordinatesPlot()

analysis.compute(
    experiment=ax_client._experiment,
    generation_strategy=ax_client._generation_strategy,
    # compute can optionally take in an Adapter directly instead of a GenerationStrategy
    adapter=None,
)

Parallel Coordinates for JV_JV_nrmse_linear

The parallel coordinates plot displays multi-dimensional data by representing each parameter as a parallel axis. This plot helps in assessing how thoroughly the search space has been explored and in identifying patterns or clusterings associated with high-performing (good) or low-performing (bad) arms. By tracing lines across the axes, one can observe correlations and interactions between parameters, gaining insights into the relationships that contribute to the success or failure of different configurations within the experiment.

# Plot the density of the exploration of the parameters
# this gives a nice visualization of where the optimizer focused its exploration and may show some correlation between the parameters
plot_dens = True
if plot_dens:
    from optimpv.posterior.posterior import *
    best_parameters = {}
    for p in optimizer.params:
        best_parameters[p.name] = p.value

    fig_dens, ax_dens = plot_density_exploration(params, optimizer = optimizer, best_parameters = best_parameters, optimizer_type = 'ax')

# rerun the simulation with the best parameters
yfit = jv.run(parameters={}) # run the simulation with the best parameters

plt.figure(figsize=(10,10))
linewidth = 2
plt.plot(X, y, 'C0-',label='Data',linewidth=linewidth, alpha=0.5)
plt.plot(X, yfit, 'C0--', label='Fitted',linewidth=linewidth)
plt.xlabel('Voltage [V]')
plt.ylabel('Current density [A m$^{-2}$]')
plt.legend()
plt.show()

# Clean up the output files (comment out if you want to keep the output files)
sim.clean_all_output(session_path)
sim.delete_folders('tmp',session_path)
# uncomment the following lines to delete specific files
sim.clean_up_output('PTAA',session_path)
sim.clean_up_output('perovskite',session_path)
sim.clean_up_output('C60',session_path)
sim.clean_up_output('simulation_setup_realPerovskite.txt',session_path)