Pymoo: Layer stack optimization with Transfer Matrix Method (TMM)

This notebook demonstrate the use of multi-objective Bayesian optimization in combination with transfer matrix modeling (TMM) to optimize the thickness of the layers in a multilayer stack and the choice of materials to maximize the average visible transmittance (AVT) and the current density (Jsc) of a solar cell.

To perform the transfer matrix modeling we use a modified version of the open-source program devoloped by McGehee’s group (Stanford University) and adapted to python by Kamil Mielczarek (University of Texas).

For more information about the transfer matrix modeling, please refer to the original paper.

# Import necessary libraries
import warnings, os, sys, shutil
# remove warnings from the output
os.environ["PYTHONWARNINGS"] = "ignore"
warnings.filterwarnings('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, matplotlib, ax, logging
from ax.utils.notebook.plotting import init_notebook_plotting, render
init_notebook_plotting() # for Jupyter notebooks

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

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

Define the parameters for the simulation

params = []
d_3 = FitParam(name = 'd_3', value = 80e-9, bounds = [40e-9, 200e-9], log_scale = False, rescale = True, value_type = 'float', type='range', display_name='d_3',unit='m')
params.append(d_3)

d_6 = FitParam(name = 'd_6', value =  10e-9, bounds = [5e-9, 20e-9], log_scale = False, rescale = True, value_type = 'float', type='range', display_name='d_6',unit='m')
params.append(d_6)

d_7 = FitParam(name = 'd_7', value =  100e-9, bounds = [50e-9, 200e-9], log_scale = False, rescale = True, value_type = 'float', type='range', display_name='d_7',unit='m')
params.append(d_7)

d_8 = FitParam(name = 'd_8', value =  10e-9, bounds = [5e-9, 20e-9], log_scale = False, rescale = True, value_type = 'float', type='range', display_name='d_8',unit='m')
params.append(d_8)

d_9 = FitParam(name = 'd_9', value =  100e-9, bounds = [50e-9, 200e-9], log_scale = False, rescale = True, value_type = 'float', type='range', display_name='d_9',unit='m')
params.append(d_9)

Run the optimization

# Initialize the agent and default device stack
layers = ['SiOx' , 'ITO' , 'ZnO' , 'PCE10_FOIC_1to1' , 'MoOx' , 'Ag', 'MoOx', 'LiF','MoOx', 'LiF','Air'] # list of layers (need to be the same than the name nk_*.csv file in the matdata folder)
thicknesses =  [0 , 100e-9 , 30e-9  , 100e-9 , 9e-9 , 8e-9, 100e-9, 100e-9, 100e-9, 100e-9, 100e-9]# list of thicknesses in nm
mat_dir = os.path.join(os.path.abspath('../'),'Data','matdata') # path to the folder containing the nk_*.csv files
lambda_min = 350e-9 # start of the wavelength range
lambda_max = 800e-9 # end of the wavelength range
lambda_step = 1e-9 # wavelength step
x_step = 1e-9 # x step
activeLayer = 3 # active layer index
spectrum = os.path.join(mat_dir ,'AM15G.txt') # path to the AM15G spectrum file
photopic_file = os.path.join(mat_dir ,'photopic_curve.txt') # path to the photopic spectrum file

from optimpv.TransferMatrix.TransferMatrixAgent import TransferMatrixAgent

TMAgent = TransferMatrixAgent(params, [None,None], layers=layers, thicknesses=thicknesses, lambda_min=lambda_min, lambda_max=lambda_max, lambda_step=lambda_step, x_step=x_step, activeLayer=activeLayer, spectrum=spectrum, mat_dir=mat_dir, photopic_file=photopic_file, exp_format=['Jsc', 'AVT'],metric=[None,None],loss=[None,None],threshold=[4,0.1],minimize=[False,False])

# Define the optimizer
optimizer = PymooOptimizer(params=params, agents=TMAgent, algorithm='NSGA2', pop_size=20, n_gen=10, name='pymoo_single_obj', verbose_logging=True,max_parallelism=100, )

res = optimizer.optimize()

[INFO 08-14 10:06:33] optimpv.pymooOptimizer: Starting optimization using NSGA2 algorithm
[INFO 08-14 10:06:33] optimpv.pymooOptimizer: Population size: 20, Generations: 10
[INFO 08-14 10:06:56] optimpv.pymooOptimizer: Generation 1: Best objective = -190.590779

==========================================================
n_gen  |  n_eval  | n_nds  |      eps      |   indicator
==========================================================
     1 |       20 |     13 |             - |             -

[INFO 08-14 10:07:21] optimpv.pymooOptimizer: Generation 2: Best objective = -190.590779

     2 |       40 |     20 |  0.0162901004 |             f

[INFO 08-14 10:07:47] optimpv.pymooOptimizer: Generation 3: Best objective = -192.145989

     3 |       60 |     20 |  0.0168206297 |         ideal

[INFO 08-14 10:08:08] optimpv.pymooOptimizer: Generation 4: Best objective = -192.594881

     4 |       80 |     20 |  0.0178237386 |         ideal

[INFO 08-14 10:08:28] optimpv.pymooOptimizer: Generation 5: Best objective = -192.594881

     5 |      100 |     20 |  0.0149210086 |             f

[INFO 08-14 10:08:48] optimpv.pymooOptimizer: Generation 6: Best objective = -193.791137

     6 |      120 |     20 |  0.0124266289 |         ideal

[INFO 08-14 10:09:08] optimpv.pymooOptimizer: Generation 7: Best objective = -193.905805

     7 |      140 |     20 |  0.0092586074 |             f

[INFO 08-14 10:09:28] optimpv.pymooOptimizer: Generation 8: Best objective = -193.962205

     8 |      160 |     20 |  0.0025586666 |         ideal

[INFO 08-14 10:09:48] optimpv.pymooOptimizer: Generation 9: Best objective = -193.962205

     9 |      180 |     20 |  0.0262990919 |         ideal

[INFO 08-14 10:10:08] optimpv.pymooOptimizer: Generation 10: Best objective = -193.962205
[INFO 08-14 10:10:08] optimpv.pymooOptimizer: Optimization completed after 11 generations
[INFO 08-14 10:10:08] optimpv.pymooOptimizer: Number of function evaluations: 200
[INFO 08-14 10:10:08] optimpv.pymooOptimizer: Final population size: 20
[INFO 08-14 10:10:08] optimpv.pymooOptimizer: Pareto front size: 20

    10 |      200 |     20 |  0.0082540085 |             f

# get the best parameters and update the params list in the optimizer and the agent
optimizer.update_params_with_best_balance() # update the params list in the optimizer with the best parameters
TMAgent.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)

Best parameters:
d_3 fitted value: 1.3807462660125542e-07
d_6 fitted value: 5.883338989931872e-09
d_7 fitted value: 5.6055640742220376e-08
d_8 fitted value: 7.14078140635034e-09
d_9 fitted value: 7.657106273782278e-08

# 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 *
    fig_dens, ax_dens = plot_density_exploration(params, optimizer = optimizer, optimizer_type = 'pymoo')

from pymoo.visualization.scatter import Scatter

Scatter().add(res.F, label="Pareto Front").show()

<pymoo.visualization.scatter.Scatter at 0x7d52404d63c0>

import matplotlib
# import itertools
from itertools import combinations
comb = list(combinations(optimizer.all_metrics, 2))
threshold_list = []
for i in range(len(optimizer.agents)):
    for j in range(len(optimizer.agents[i].threshold)):
        threshold_list.append(optimizer.agents[i].threshold[j])
threshold_comb = list(combinations(threshold_list, 2))
pareto = np.asarray(res.F)

cm = matplotlib.colormaps.get_cmap('viridis')
df = optimizer.get_df_from_pymoo() # get the dataframe from the optimizer
for metric, mini in zip(optimizer.all_metrics,optimizer.all_minimize):
    if not mini:
        df[metric] = -df[metric]

dum_dic = {}
for i , metr in enumerate(optimizer.all_metrics):
    if i not in df.keys():
        if not optimizer.all_minimize[i]:
            dum_dic[metr] = -pareto[:, i]

df_pareto = pd.DataFrame(dum_dic)
dum_dic = TMAgent.run_Ax(parameters={})
best_balanced = [dum_dic[metr] for metr in optimizer.all_metrics]

for c,t_c in zip(comb,threshold_comb):
    plt.figure(figsize=(10, 10))
    plt.scatter(df[c[0]],df[c[1]],c=df.index, cmap=cm, marker='o', s=100) # plot the points with color according to the iteration
    cbar = plt.colorbar()
    cbar.set_label('Iteration')
    sorted_df = df_pareto.sort_values(by=c[0])
    plt.plot(sorted_df[c[0]],sorted_df[c[1]],'r')

    plt.xlabel(c[0])
    plt.ylabel(c[1])
    plt.xscale('log')
    plt.yscale('log')


    plt.show()