OPV 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.
[1]:
# 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:11:27] ax.utils.notebook.plotting: Injecting Plotly library into cell. Do not overwrite or delete cell.
[INFO 08-13 18:11:27] 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
[2]:
# 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','JVrealOPV')))
simulation_setup_filename = 'simulation_setup_PM6_L8BO.txt'
simulation_setup = os.path.join(session_path, simulation_setup_filename)
# path to the layer files defined in the simulation_setup file
l1 = 'ZnO.txt'
l2 = 'PM6_L8BO.txt'
l3 = 'BM_HTL.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)
# 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_PM6_L8BO.dat'), sep=' ')
X = df[['Vext','Gfrac']].values
y = df[['Jext']].values.reshape(-1)
Gfracs = pd.unique(df['Gfrac'])
X_1sun = df[['Vext','Gfrac']].values[df['Gfrac'] == 1.0]
y_1sun = df[['Jext']].values[df['Gfrac'] == 1.0].reshape(-1)
# get 1sun Jsc as it will help narrow down the range for Gehp
Jsc_1sun = np.interp(0.0, X_1sun[:,0], y_1sun)
minJ = np.min(y_1sun)
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))
# get Voc
Voc_1sun = np.interp(0.0, y_1sun, X_1sun[:,0])
# Filter some data
Vmin = -0.1
Jmax = abs(Jsc_1sun)
idx = np.where((X[:,0] > Vmin) & (y < Jmax))
y = y[idx]
X = X[idx]
# plot the data for each Gfrac
plt.figure(figsize=(10,10))
viridis = plt.get_cmap('viridis', len(Gfracs))
for i, Gfrac in enumerate(Gfracs):
plt.plot(X[X[:,1] == Gfrac][:,0], y[X[:,1] == Gfrac], '-', label = f'Gfrac = {Gfrac}', color = viridis(i))
plt.grid()
plt.xlabel('Voltage [V]')
plt.ylabel('Current density [A m$^{-2}$]')
plt.legend()
plt.show()
Define the parameters for the simulation
[3]:
params = [] # list of parameters to be optimized
mun = FitParam(name = 'l2.mu_n', value = 7e-8, bounds = [1e-9,1e-6], log_scale = True, value_type = 'float', fscale = None, rescale = False, display_name=r'$\mu_n$', unit='m$^2$ V$^{-1}$s$^{-1}$', axis_type = 'log', force_log = True)
params.append(mun)
mup = FitParam(name = 'l2.mu_p', value = 5e-8, bounds = [1e-9,1e-6], log_scale = True, value_type = 'float', fscale = None, rescale = False, display_name=r'$\mu_p$', unit=r'm$^2$ V$^{-1}$s$^{-1}$', axis_type = 'log', force_log = True)
params.append(mup)
bulk_tr = FitParam(name = 'l2.N_t_bulk', value = 1e20, bounds = [1e16,8e22], log_scale = True, value_type = 'float', fscale = None, rescale = False, display_name=r'$N_{T}$', unit=r'm$^{-3}$', axis_type = 'log', force_log = True)
params.append(bulk_tr)
preLangevin = FitParam(name = 'l2.preLangevin', value = 1e-2, bounds = [1e-3,1], log_scale = True, value_type = 'float', fscale = None, rescale = False, display_name=r'$\gamma_{pre}$', unit=r'', axis_type = 'log', force_log = True)
params.append(preLangevin)
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)
G_ehp = FitParam(name = 'l2.G_ehp', value = G_ehp_calc, bounds = [G_ehp_calc*0.95,G_ehp_max*1.05], log_scale = False, value_type = 'float', fscale = None, rescale = False, display_name=r'$G_{ehp}$', unit=r'm$^{-3}$ s$^{-1}$', axis_type = 'linear', force_log = False)
params.append(G_ehp)
# save the original parameters for later
params_orig = copy.deepcopy(params)
Run the optimization
[4]:
# 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 = True, max_jobs = 3, metric = metric, loss = loss)
[5]:
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,45], 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)
[6]:
# optimizer.optimize() # run the optimization with ax
optimizer.optimize_turbo() # run the optimization with turbo
[INFO 08-13 18:11:28] optimpv.axBOtorchOptimizer: Starting optimization with 46 batches and a total of 100 trials
[INFO 08-13 18:11:28] optimpv.axBOtorchOptimizer: Starting Sobol batch 1 with 10 trials
[INFO 08-13 18:11:30] optimpv.axBOtorchOptimizer: Finished Sobol
[INFO 08-13 18:11:32] optimpv.axBOtorchOptimizer: Finished Turbo batch 2 with 2 trials with current best value: 8.12e-02, TR length: 8.00e-01
[INFO 08-13 18:11:32] optimpv.axBOtorchOptimizer: Finished Turbo batch 3 with 2 trials with current best value: 8.12e-02, TR length: 8.00e-01
[INFO 08-13 18:11:33] optimpv.axBOtorchOptimizer: Finished Turbo batch 4 with 2 trials with current best value: 8.12e-02, TR length: 8.00e-01
[INFO 08-13 18:11:34] optimpv.axBOtorchOptimizer: Finished Turbo batch 5 with 2 trials with current best value: 8.12e-02, TR length: 8.00e-01
[INFO 08-13 18:11:34] optimpv.axBOtorchOptimizer: Finished Turbo batch 6 with 2 trials with current best value: 8.12e-02, TR length: 8.00e-01
[INFO 08-13 18:11:35] optimpv.axBOtorchOptimizer: Finished Turbo batch 7 with 2 trials with current best value: 8.12e-02, TR length: 8.00e-01
[INFO 08-13 18:11:35] optimpv.axBOtorchOptimizer: Finished Turbo batch 8 with 2 trials with current best value: 8.12e-02, TR length: 8.00e-01
[INFO 08-13 18:11:36] optimpv.axBOtorchOptimizer: Finished Turbo batch 9 with 2 trials with current best value: 8.12e-02, TR length: 8.00e-01
[INFO 08-13 18:11:36] optimpv.axBOtorchOptimizer: Finished Turbo batch 10 with 2 trials with current best value: 7.10e-02, TR length: 8.00e-01
[INFO 08-13 18:11:37] optimpv.axBOtorchOptimizer: Finished Turbo batch 11 with 2 trials with current best value: 7.10e-02, TR length: 8.00e-01
[INFO 08-13 18:11:38] optimpv.axBOtorchOptimizer: Finished Turbo batch 12 with 2 trials with current best value: 7.10e-02, TR length: 8.00e-01
[INFO 08-13 18:11:38] optimpv.axBOtorchOptimizer: Finished Turbo batch 13 with 2 trials with current best value: 7.10e-02, TR length: 8.00e-01
[INFO 08-13 18:11:39] optimpv.axBOtorchOptimizer: Finished Turbo batch 14 with 2 trials with current best value: 7.10e-02, TR length: 8.00e-01
[INFO 08-13 18:11:40] optimpv.axBOtorchOptimizer: Finished Turbo batch 15 with 2 trials with current best value: 5.78e-02, TR length: 8.00e-01
[INFO 08-13 18:11:40] optimpv.axBOtorchOptimizer: Finished Turbo batch 16 with 2 trials with current best value: 5.31e-02, TR length: 8.00e-01
[INFO 08-13 18:11:41] optimpv.axBOtorchOptimizer: Finished Turbo batch 17 with 2 trials with current best value: 3.74e-02, TR length: 1.60e+00
[INFO 08-13 18:11:42] optimpv.axBOtorchOptimizer: Finished Turbo batch 18 with 2 trials with current best value: 3.74e-02, TR length: 1.60e+00
[INFO 08-13 18:11:42] optimpv.axBOtorchOptimizer: Finished Turbo batch 19 with 2 trials with current best value: 3.74e-02, TR length: 1.60e+00
[INFO 08-13 18:11:43] optimpv.axBOtorchOptimizer: Finished Turbo batch 20 with 2 trials with current best value: 3.69e-02, TR length: 1.60e+00
[INFO 08-13 18:11:44] optimpv.axBOtorchOptimizer: Finished Turbo batch 21 with 2 trials with current best value: 3.69e-02, TR length: 1.60e+00
[INFO 08-13 18:11:45] optimpv.axBOtorchOptimizer: Finished Turbo batch 22 with 2 trials with current best value: 3.69e-02, TR length: 1.60e+00
[INFO 08-13 18:11:46] optimpv.axBOtorchOptimizer: Finished Turbo batch 23 with 2 trials with current best value: 2.25e-02, TR length: 1.60e+00
[INFO 08-13 18:11:46] optimpv.axBOtorchOptimizer: Finished Turbo batch 24 with 2 trials with current best value: 2.25e-02, TR length: 1.60e+00
[INFO 08-13 18:11:47] optimpv.axBOtorchOptimizer: Finished Turbo batch 25 with 2 trials with current best value: 2.25e-02, TR length: 1.60e+00
[INFO 08-13 18:11:48] optimpv.axBOtorchOptimizer: Finished Turbo batch 26 with 2 trials with current best value: 2.25e-02, TR length: 1.60e+00
[INFO 08-13 18:11:49] optimpv.axBOtorchOptimizer: Finished Turbo batch 27 with 2 trials with current best value: 2.25e-02, TR length: 1.60e+00
[INFO 08-13 18:11:49] optimpv.axBOtorchOptimizer: Finished Turbo batch 28 with 2 trials with current best value: 1.70e-02, TR length: 1.60e+00
[INFO 08-13 18:11:50] optimpv.axBOtorchOptimizer: Finished Turbo batch 29 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:51] optimpv.axBOtorchOptimizer: Finished Turbo batch 30 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:52] optimpv.axBOtorchOptimizer: Finished Turbo batch 31 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:52] optimpv.axBOtorchOptimizer: Finished Turbo batch 32 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:53] optimpv.axBOtorchOptimizer: Finished Turbo batch 33 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:54] optimpv.axBOtorchOptimizer: Finished Turbo batch 34 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:54] optimpv.axBOtorchOptimizer: Finished Turbo batch 35 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:55] optimpv.axBOtorchOptimizer: Finished Turbo batch 36 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:56] optimpv.axBOtorchOptimizer: Finished Turbo batch 37 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:57] optimpv.axBOtorchOptimizer: Finished Turbo batch 38 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:57] optimpv.axBOtorchOptimizer: Finished Turbo batch 39 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:58] optimpv.axBOtorchOptimizer: Finished Turbo batch 40 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:11:59] optimpv.axBOtorchOptimizer: Finished Turbo batch 41 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:12:00] optimpv.axBOtorchOptimizer: Finished Turbo batch 42 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:12:01] optimpv.axBOtorchOptimizer: Finished Turbo batch 43 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:12:02] optimpv.axBOtorchOptimizer: Finished Turbo batch 44 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:12:02] optimpv.axBOtorchOptimizer: Finished Turbo batch 45 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:12:03] optimpv.axBOtorchOptimizer: Finished Turbo batch 46 with 2 trials with current best value: 1.15e-02, TR length: 1.60e+00
[INFO 08-13 18:12:04] optimpv.axBOtorchOptimizer: Finished Turbo batch 47 with 2 trials with current best value: 7.35e-03, TR length: 1.60e+00
[INFO 08-13 18:12:04] optimpv.axBOtorchOptimizer: Turbo is terminated as the max number (100) of trials is reached
[7]:
# 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: 2.7549362125794885e-08
l2.mu_p fitted value: 3.115685253947716e-08
l2.N_t_bulk fitted value: 2.4287616095425905e+19
l2.preLangevin fitted value: 0.0013405850516945467
R_series fitted value: 6.218517272918729e-05
R_shunt fitted value: 41.93809160351753
l2.G_ehp fitted value: 1.467606101647734e+28
SimSS command line:
./simss -l2.mu_n 2.7549362125794885e-08 -l2.mu_p 3.115685253947716e-08 -l2.N_t_bulk 2.4287616095425905e+19 -l2.preLangevin 0.0013405850516945467 -R_series 6.218517272918729e-05 -R_shunt 41.93809160351753 -l2.G_ehp 1.467606101647734e+28
[8]:
# 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.007349574643789233 at iteration 101
[13]:
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.
[10]:
# 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')
[11]:
# rerun the simulation with the best parameters
yfit = jv.run(parameters={}) # run the simulation with the best parameters
viridis = plt.get_cmap('viridis', len(Gfracs))
plt.figure(figsize=(10,10))
linewidth = 2
for idx, Gfrac in enumerate(Gfracs[::-1]):
plt.plot(X[X[:,1]==Gfrac,0],y[X[:,1]==Gfrac],label='Gfrac = '+str(Gfrac),color=viridis(idx),alpha=0.5,linewidth=linewidth)
plt.plot(X[X[:,1]==Gfrac,0],yfit[X[:,1]==Gfrac],label='Gfrac = '+str(Gfrac)+' fit',linestyle='--',color=viridis(idx),linewidth=linewidth)
plt.xlabel('Voltage [V]')
plt.ylabel('Current density [A m$^{-2}$]')
plt.legend()
plt.show()
[12]:
# 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('ZnO',session_path)
sim.clean_up_output('PM6_L8BO',session_path)
sim.clean_up_output('BM_HTL',session_path)
sim.clean_up_output('simulation_setup_PM6_L8BO',session_path)