import numpy as np
from solcore import material, si
import matplotlib.pyplot as plt
from solcore.optics.tmm import OptiStack
from solcore.optics.tmm import calculate_rat
Import the DE implementation
from solcore.optimization import PDE, DE
from solcore.light_source import LightSource
from solcore.solar_cell import SolarCell
from solcore.structure import Junction, Layer
from solcore.solar_cell_solver import solar_cell_solver
from solcore.constants import q, kb
from solcore.material_system import create_new_material
from solcore.absorption_calculator import search_db
The "if name == "main" construction is used to avoid issues with parallel processing on Windows.
The issue arises because the multiprocessing module uses a different process on Windows than on UNIX
systems which will throw errors if this construction is not used.
Optimizing a four-junction cell using SiGeSn as the second-lowest bandgap material.
First, using a purely optical TMM simulation to calculate the photogenerated current in each sub-cell. The thing to
optimize is then the current of the current-limiting cell in the structure; in other words we want to maximize the
lowest sub-cell current, to achieve current-matching with the highest possible current.
Since differential evolution does a minimization, we are actually minimizing the negative of
this value.
Once we have good initial values for the total layer thicknesses, use full electrical simulations to determine the
n and p type layer thicknesses to calculate a maximum possible efficiency for the 4J device.
To use yabox/the optimization module for the DE, we need to define a class which sets up the problem and has an 'evaluate' function, which
will actually calculate the value we are trying to minimize for a set of parameters.
add SiGeSn optical constants to the database, and search the refractiveindex.info for the Ta2O5 data we want to use and save the pageid to use later:
create_new_material('SiGeSn', 'SiGeSn_n.txt', 'SiGeSn_k.txt', 'SiGeSn_params.txt') # Note: comment out this line after the material
# has been added to avoid being asked if you want to overwrite it.
Ta2O5_pageid = str(search_db("Ta2O5/Rodriguez-de Marcos")[0][0])
define class for the optimization:
def evaluate(self, x):
# x[0] = MgF2 thickness (anti-reflection coating)
# x[1] = Ta2O5 thickness (anti-reflection coating)
# x[2] = InGaP (top junction) thickness
# x[3] = GaAs (second junction) thickness
# x[4] = SiGeSn (third junction) thickness
# keep the thickness of the bottom cell constant; from an optical point of view, this should be infinitely thick
SC = [[x[0]] + self.MgF2, [x[1]] + self.Ta2O5, [x[2]] + self.InGaP, [x[3]] + self.GaAs, [x[4]] + self.SiGeSn, [300e3]+self.Ge]#, [x[5]] + self.Ge]
# create the OptiStack
full_stack = OptiStack(SC, no_back_reflection=False)
# calculate reflection, transmission, and absorption in each layer. We are specifying that the last layer,
# a very thick Ge substrate, should be treated incoherently, otherwise we would see very narrow, unphysical oscillations
# in the R/A/T spectra.
RAT = calculate_rat(full_stack, self.wl, no_back_reflection=False, coherent=False,
coherency_list=['c', 'c', 'c', 'c', 'c', 'i'])
# extract absorption per layer
A_InGaP = RAT['A_per_layer'][3]
A_GaAs = RAT['A_per_layer'][4]
A_SiGeSn = RAT['A_per_layer'][5]
A_Ge = RAT['A_per_layer'][6]
## calculate photo-generated currents using the AM1.5 G spectrum for each layer
Jsc_InGaP = 0.1 * q * np.trapz(A_InGaP* self.spectr, self.wl)
Jsc_GaAs = 0.1 * q * np.trapz(A_GaAs * self.spectr, self.wl)
Jsc_SiGeSn = 0.1 * q * np.trapz(A_SiGeSn* self.spectr, self.wl)
Jsc_Ge = 0.1 * q* np.trapz(A_Ge * self.spectr, self.wl)
# find the limiting current by checking which junction has the lowest current. Then take the negative since
# we need to minimize (not maximize)
limiting_Jsc = -min([Jsc_InGaP, Jsc_GaAs, Jsc_SiGeSn, Jsc_Ge])
return limiting_Jsc
def plot(self, x):
# this does basically what evaluate() does, but plots the results
SC = [[x[0]] + self.MgF2, [x[1]] + self.Ta2O5, [x[2]] + self.InGaP,
[x[3]] + self.GaAs, [x[4]] + self.SiGeSn, [300e3] + self.Ge]
full_stack = OptiStack(SC, no_back_reflection=False)
RAT = calculate_rat(full_stack, self.wl, no_back_reflection=False,
coherent=False, coherency_list=['c', 'c', 'c', 'c', 'c', 'i'])
A_InGaP = RAT['A_per_layer'][3]
A_GaAs = RAT['A_per_layer'][4]
A_SiGeSn = RAT['A_per_layer'][5]
A_Ge = RAT['A_per_layer'][6]
plt.figure()
plt.plot(self.wl, A_InGaP, label='InGaP')
plt.plot(self.wl, A_GaAs, label='A_GaAs')
plt.plot(self.wl, A_SiGeSn, label='SiGeSn')
plt.plot(self.wl, A_Ge, label = 'Ge')
plt.plot(self.wl, RAT['T'], label='T')
plt.plot(self.wl, RAT['R'], label='R')
plt.legend()
plt.xlabel('Wavelength (nm)')
plt.ylabel('R/A/T')
plt.show()
class calc_min_Jsc_DA():
def __init__(self, ARC_thickness):
self.ARC = ARC_thickness
def make_cell(self, x):
#x[0]: total InGaP thickness
#x[1]: total InGaAs thickness
#x[2]: total SiGeSn thickness
#x[3]: total Ge thickness
#x[4]: InGaP n thickness
#x[5]: InGaAs n thickness
#x[6]: SiGeSn n thickness
#x[7]: Ge n thickness
e_charge = si('1eV')
# materials
SiGeSn = material('SiGeSn')
GaAs = material('GaAs')
InGaP = material('GaInP')
Ge = material('Ge')
MgF2 = material('MgF2')()
Ta2O5 = material(Ta2O5_pageid, nk_db=True)()
AlInP = material("AlInP")
window_material = AlInP(Al=0.52)
GaInP_mobility_h = 0.03 #
GaInP_lifetime_h = 1e-8
GaInP_D_h = GaInP_mobility_h * kb * 300 / e_charge
GaInP_L_h = np.sqrt(GaInP_D_h * GaInP_lifetime_h)
GaInP_mobility_e = 0.015
GaInP_lifetime_e = 1e-8
GaInP_D_e = GaInP_mobility_e * kb * 300 / e_charge
GaInP_L_e = np.sqrt(GaInP_D_e * GaInP_lifetime_e)
top_cell_n_material = InGaP(In=0.5, Nd=si("2e18cm-3"), hole_diffusion_length=GaInP_L_h, hole_mobility=GaInP_mobility_h)
top_cell_p_material = InGaP(In=0.5, Na=si("2e17cm-3"), electron_diffusion_length=GaInP_L_e, electron_mobility=GaInP_mobility_e)
# MID CELL - GaAs
GaAs_mobility_h = 0.85 #
GaAs_lifetime_h = 1e-8
GaAs_D_h = GaAs_mobility_h * kb * 300 / e_charge
GaAs_L_h = np.sqrt(GaAs_D_h * GaAs_lifetime_h)
GaAs_mobility_e = 0.08
GaAs_lifetime_e = 1e-8
GaAs_D_e = GaAs_mobility_e * kb * 300 / e_charge
GaAs_L_e = np.sqrt(GaAs_D_e * GaAs_lifetime_e)
mid_cell_n_material = GaAs(Nd=si("1e18cm-3"), hole_diffusion_length=GaAs_L_h, hole_mobility=GaAs_mobility_h)
mid_cell_p_material = GaAs(Na=si("1e17cm-3"), electron_diffusion_length=GaAs_L_e, electron_mobility=GaAs_mobility_e)
SiGeSn.band_gap = si('0.77eV') # from PL measurement
SiGeSn_L_h = si('0.35um')
SiGeSn_L_e = si('5um')
SiGeSn_lifetime_e = 1e-6
SiGeSn_lifetime_h = 1e-6
SiGeSn_mobility_h = SiGeSn_L_h ** 2 * e_charge / (SiGeSn_lifetime_h * kb * 300)
SiGeSn_mobility_e = SiGeSn_L_e ** 2 * e_charge / (SiGeSn_lifetime_e * kb * 300)
pen_cell_n_material = SiGeSn(Nd=si("1e18cm-3"), hole_diffusion_length=SiGeSn_L_h,
relative_permittivity=16, hole_mobility=SiGeSn_mobility_h)
pen_cell_p_material = SiGeSn(Na=si("1e17cm-3"), electron_diffusion_length=SiGeSn_L_e,
relative_permittivity=16, electron_mobility=SiGeSn_mobility_e)
# Ge_mobility_h = 0.38 #
Ge_lifetime_h = 1e-6
# Ge_D_h = Ge_mobility_h*kb*300/e_charge
# Ge_L_h = np.sqrt(Ge_D_h*Ge_lifetime_h)
Ge_L_h = si('500nm')
Ge_mobility_h = Ge_L_h ** 2 * e_charge / (Ge_lifetime_h * kb * 300)
Ge_mobility_e = 0.18
Ge_lifetime_e = 1e-6
Ge_D_e = Ge_mobility_e * kb * 300 / e_charge
Ge_L_e = np.sqrt(Ge_D_e * Ge_lifetime_e)
class calc_min_Jsc():
def __init__(self):
# initialize an instance of the class; set some information which will be used in each iteration of the calculation:
# materials, wavelengths, the light source
T = 298
wl = np.linspace(300, 1900, 800)
# materials
SiGeSn = material('SiGeSn')(T=T)
GaAs = material('GaAs')(T=T)
InGaP = material('GaInP')(In=0.5, T=T)
Ge = material('Ge')(T=T)
# make these attributes of 'self' so they can be accessed by the class object
# here I am also creating lists of wavelengths and corresponding n and k data from
# the Solcore materials - the reason for this is that there is currently an issue with using the Solcore
# material class in parallel computations. Thus the information for the n and k data is saved here as a list
# rather than a material object (see the documentation of OptiStack for the different acceptable formats
# to pass optical constants for an OptiStack
self.wl = wl
self.SiGeSn = [self.wl, SiGeSn.n(self.wl*1e-9), SiGeSn.k(self.wl*1e-9)]
self.Ge = [self.wl, Ge.n(self.wl*1e-9), Ge.k(self.wl*1e-9)]
self.InGaP = [self.wl, InGaP.n(self.wl*1e-9), InGaP.k(self.wl*1e-9)]
self.GaAs = [self.wl, GaAs.n(self.wl*1e-9), GaAs.k(self.wl*1e-9)]
self.MgF2 = [self.wl, material('MgF2')().n(self.wl*1e-9), material('MgF2')().k(self.wl*1e-9)]
self.Ta2O5 = [self.wl, material(Ta2O5_pageid,
nk_db=True)().n(self.wl*1e-9), material(Ta2O5_pageid,
nk_db=True)().k(self.wl*1e-9)]
# assuming an AM1.5G spectrum
self.spectr = LightSource(source_type='standard', version='AM1.5g', x=self.wl,
output_units='photon_flux_per_nm', concentration=1).spectrum(self.wl)[1]
bot_cell_n_material = Ge(Nd=si("2e18cm-3"), hole_diffusion_length=Ge_L_h, hole_mobility=Ge_mobility_h)
bot_cell_p_material = Ge(Na=si("1e17cm-3"), electron_diffusion_length=Ge_L_e, electron_mobility=Ge_mobility_e)
# And, finally, we put everything together, adding also the surface recombination velocities. We also add some shading
# due to the metallisation of the cell = 8%, and indicate it has an area of 0.7x0.7 mm2 (converted to m2)
solar_cell = SolarCell([
# Layer(si('110nm'), material = MgF2), Layer(si('55nm'), material = ZnS),
Layer(si(self.ARC[0], 'nm'), material=MgF2), Layer(si(self.ARC[1], 'nm'), material=Ta2O5),
Junction([Layer(si(25, 'nm'), material=window_material, role='window'),
Layer(si(x[4], 'nm'), material=top_cell_n_material, role='emitter'),
Layer(si(x[0]-x[4], 'nm'), material=top_cell_p_material, role='base'),
], sn=1, sp=1, kind='DA'),
Junction([Layer(si(x[5], 'nm'), material=mid_cell_n_material, role='emitter'),
Layer(si(x[1]-x[5], 'nm'), material=mid_cell_p_material, role='base'),
], sn=1, sp=1, kind='DA'),
Junction([Layer(si(x[6], 'nm'), material=pen_cell_n_material, role='emitter'),
Layer(si(x[2]-x[6], 'nm'), material=pen_cell_p_material, role='base'),
], sn=1, sp=1, kind='DA'),
Junction([Layer(si(x[7], 'nm'), material=bot_cell_n_material, role='emitter'),
Layer(si(x[3]-x[7], 'nm'), material=bot_cell_p_material, role='base'),
], sn=1, sp=1, kind='DA'),
], shading=0.0, substrate=bot_cell_n_material)
return solar_cell
def evaluate(self, x):
light_source = LightSource(source_type='standard', version='AM1.5g')
wl = np.linspace(300, 1850, 500) * 1e-9
solar_cell = self.make_cell(x)
position = [1e-10] * 10 + [5e-8]
V = np.linspace(0, 3.5, 300)
solar_cell_solver(solar_cell, 'iv',
user_options={'voltages': V, 'light_iv': True, 'wavelength': wl, 'mpp': True,
'light_source': light_source,
'optics_method': 'TMM', 'BL_correction': True, 'position': position})
efficiency = solar_cell.iv["Eta"]
# print('Efficiency =', efficiency)
return -efficiency
def plot(self, x):
light_source = LightSource(source_type='standard', version='AM1.5g')
wl = np.linspace(300, 1850, 500) * 1e-9
solar_cell = self.make_cell(x)
position = [1e-10] * 10 + [5e-8]
V = np.linspace(0, 3.5, 300)
solar_cell_solver(solar_cell, 'iv',
user_options={'voltages': V, 'light_iv': True, 'wavelength': wl, 'mpp': True,
'light_source': light_source,
'optics_method': 'TMM', 'BL_correction': True, 'position': position})
efficiency = solar_cell.iv["Eta"]
pmax = solar_cell.iv["Pmpp"]
ff = solar_cell.iv["FF"]
voc = solar_cell.iv["Voc"]
isc = solar_cell.iv["Isc"]
plt.figure()
plt.plot(V, solar_cell.iv['IV'][1] / 10, 'k', linewidth=3, label='Total')
plt.plot(V, -solar_cell[2].iv(V) / 10, 'b', label='GaInP')
plt.plot(V, -solar_cell[3].iv(V) / 10, 'g', label='GaAs')
plt.plot(V, -solar_cell[4].iv(V) / 10, 'r', label='SiGeSn')
plt.plot(V, -solar_cell[5].iv(V) / 10, 'y', label='Ge')
plt.text(2, 10, '$\eta = $' + str(round(efficiency * 100, 1)) + '%')
plt.text(2, 8,'Pmax='+str(round(pmax,1))+'W/m$^2$')
plt.text(2, 9, 'FF = ' + str(round(ff * 100, 1)) + '%')
plt.text(2,7,'Voc='+str(round(voc,1))+'V')
plt.text(2,6, 'Jsc='+str(round(0.1*isc,1))+'mA/cm$^2$')
plt.legend()
plt.ylim(0, 18)
plt.xlim(0, 3.5)
plt.ylabel('Current (mA/cm$^2$)')
plt.xlabel('Voltage (V)')
plt.show()
# print('Efficiency =', efficiency)
solar_cell_solver(solar_cell, 'qe',
user_options={'wavelength': wl, 'optics_method': 'TMM', 'BL_correction': True, 'position': position})
plt.figure()
plt.plot(wl * 1e9, solar_cell[2].eqe(wl) * 100, 'b', label='InGaP')
plt.plot(wl * 1e9, solar_cell[3].eqe(wl) * 100, 'g', label='InGaAs')
plt.plot(wl * 1e9, solar_cell[4].eqe(wl) * 100, 'r', label='SiGeSn')
plt.plot(wl * 1e9, solar_cell[5].eqe(wl) * 100, 'y', label='Ge')
plt.plot(wl * 1e9, solar_cell.absorbed * 100, 'k--', label='Absorption')
# plt.plot(wl * 1e9, solar_cell[5].eqe(wl)*100, 'y', label='Ge')
plt.legend(loc='upper right')
plt.xlim(290, 1850)
plt.ylim(0, 100)
plt.ylabel('EQE (%)')
plt.xlabel('Wavelength (nm)')
plt.show()
def main():
# number of iterations for Differential Evolution optimization of the optical stack
maxiters=300
# make an instance of the class the DE algorithm is going to use, as defined above
DE_class = calc_min_Jsc()
# pass the function which will be minimized to the PDE (parallel differential evolution) solver. PDE calculates the
# results for each population in parallel to speed up the overall process. The bounds argument sets upper and lower bounds
# for each parameter. PDE_obj contains all the information to run the DE but does not actually invoke the calculation....
PDE_obj = PDE(DE_class.evaluate, bounds=[[10,150], [10,105], [200, 1000], [500, 10000], [500, 10000]], maxiters=maxiters)
# this will run the calculation in parallel, with all the cores available. If you don't want this, use 'DE' instead of 'PDE'
# to run the DE, use the .solve() function of the PDE object class
res = PDE_obj.solve()
# PDE_obj.solve() returns 5 things:
# res[0] is a list of the parameters which gave the minimized value
# res[1] is that minimized value
# res[2] is the evolution of the best population (the best population from each iteration
# res[3] is the evolution of the minimized value, i.e. the fitness over each iteration
# res[4] is the evolution of the mean fitness over the iterations
# best population:
best_pop = res[0]
print('parameters for best result:', best_pop, '\n', 'optimized Jsc value (mA/cm2):', -res[1])
# plot the result at these best parameters
DE_class.plot(best_pop)
best_pop_evo = res[2]
best_fitn_evo = res[3]
mean_fitn_evo = res[4]
final_fitness = res[1]
# plot evolution of the fitness of the best population per iteration
plt.figure()
plt.plot(-best_fitn_evo, '-k')
plt.xlabel('iteration')
plt.ylabel('fitness')
plt.title('Best fitness')
plt.show()
# plot evolution of the mean fitness of the population per iteration
plt.figure()
plt.plot(-mean_fitn_evo, '-k')
plt.xlabel('iteration')
plt.ylabel('fitness')
plt.title('Mean fitness')
plt.show()
# these plots show that the fitness of the best population 'jumps' every few iterations as a new best population is found,
# while the mean fitness converges slowly as the whole population gradually improves
## Now that the layer thicknesses have been optimized from an optical point of view, we want to design the device (or
# at least a simplified version, by calculating a more realistic EQE. Obviously additional parameters like the doping of the
# layers could be varied too.
# x[0]: total InGaP thickness
# x[1]: total InGaAs thickness
# x[2]: total SiGeSn thickness
# x[3]: total Ge thickness
# x[4]: InGaP n thickness
# x[5]: InGaAs n thickness
# x[6]: SiGeSn n thickness
# x[7]: Ge n thickness
# keep the ARC thicknesses fixed at the values obtained in the optical simulation
# generate upper and lower bounds: total layer thickness between 75% and 125% of values fitted in TMM calculation. Ge
# starting value is 200 um
starting_params = np.append(best_pop[2:], [200000])
lower = 0.75*starting_params
upper = 1.25*starting_params
# upper and lower bounds for the n-type (highly doped) layers
lower_ntype = [20, 20, 20, 20]
upper_ntype = [200, 300, 300, 500]
all_lower = np.append(lower, lower_ntype)
all_upper = np.append(upper, upper_ntype)
# full list of bounds
all_bounds = np.stack((all_lower, all_upper)).T
# DE calculation for the electrical simulation
maxiters_DA = 10
DE_class_DA = calc_min_Jsc_DA(best_pop[0:2])
# default population size = 5*number of params
PDE_obj_DA = PDE(DE_class_DA.evaluate, bounds=all_bounds, maxiters=maxiters_DA)
# solve, i.e. minimize the problem
res_DA = PDE_obj_DA.solve()
best_pop_DA = res_DA[0]
print('parameters for best result:', best_pop_DA, 'optimized efficiency (%)', res_DA[1]*100)
# plot the result at these best parameters
DE_class_DA.plot(best_pop_DA)
best_pop_evo = res_DA[2]
best_fitn_evo = res_DA[3]
mean_fitn_evo = res_DA[4]
final_fitness = res_DA[1]
# plot evolution of the fitness of the best population per iteration, and the mean fitness
plt.figure()
plt.plot(-best_fitn_evo, '-k')
plt.xlabel('iteration')
plt.ylabel('fitness')
plt.title('Best fitness')
plt.show()
plt.figure()
plt.plot(-mean_fitn_evo, '-k')
plt.xlabel('iteration')
plt.ylabel('fitness')
plt.title('Mean fitness')
plt.show()
if __name__ == '__main__':
main()