examples/advanced_examples/differential_evolution_4Jcell.ipynb

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