Problem Multi-objective optimization problem

01/2021: updated to pyGDM v1.1+

In this example, we will search for a multi-resonant gold nano-structure. The structure we seek has two specific plasmon resonances, \(\lambda_1 = 700\,\)nm for polarization along X and \(\lambda_2 = 1300\,\)nm for an incident polarization along Y. We will determine if a plasmon resonance exists via the scattering efficiency of the structure.

Load the modules

from pyGDM2 import structures
from pyGDM2 import materials
from pyGDM2 import fields
from pyGDM2 import linear
from pyGDM2 import tools
from pyGDM2 import propagators
from pyGDM2 import visu
from pyGDM2 import core

from pyGDM2 import EO
from pyGDM2.EO.models import RectangularAntenna, CrossAntenna
from pyGDM2.EO.core import run_eo

Define a multi-objective problem

Now we need to define our multi-objective problem. We do this by inheriting from EO.problems.BaseProblem and overriding objective_function:

from pyGDM2.EO.problems import BaseProblem

class ownMultiObjectiveProblem(BaseProblem):
    """maximize scattering efficiency (Q_sc) of nano-structure at two wavelengths"""

    def __init__(self, model):
        """
        In the Problem constructor, the constructor of the parent
        class "BaseProblem" must be called, passing the `model` object.:
        """
        super(self.__class__, self).__init__(model)


    def objective_function(self, params):
        """calculate the scat. efficiency for 2 incident field configs

        This method implements the multi-objective problem:
        The objective function returns an array with several objective
        function values (the "objective vector").
        Just as in a single-objective scenario, `objective_function`
        is function of the free parameters `params`.

        Use the `field_index` parameter to select the different incident
        field configurations. In this example problem we use the first
        and last available `field_index`.

        """
        ## --- generate the structure geometry as function of `params`
        self.model.generate_structure(params)

        ## --- GDM simulation and cross-section calculation
        core.scatter(self.model.sim, verbose=0)

        ## --- calculate scattering efficiency
        ##     take first and second last field config
        ##     (e.g. two crossed polarizations at two wavelengths)
        i_field1 = 0
        i_field2 = len(self.model.sim.E) - 1

        ext_cs1, sca_cs1, abs_cs1 = linear.extinct(self.model.sim, field_index=i_field1)
        ext_cs2, sca_cs2, abs_cs2 = linear.extinct(self.model.sim, field_index=i_field2)
        geom_cs = tools.get_geometric_cross_section(self.model.sim)

        ## --- scattering efficieny (scat eff. = scat CS / geometric CS)
        return [ float(sca_cs1/geom_cs), float(sca_cs2/geom_cs) ]

Optimization setup

The setup of the optimization is straightforward. There are only two small differences to the single-objective optimizations:

1. we now configure the pyGDM simulation for two wavelengths and two polarizations

2. we use an algorithm which can treat multi-objective problems. In this example we will use the pyGMO-algorithm moead

We will also use a slightly more complex geometry model: A cross-like structure, which consists of two “crossed” rectangles, hence possesses four free parameters: two widths and two lengths (we fix the height).

#==============================================================================
# Setup pyGDM part
#==============================================================================
## ---------- Setup structure
mesh = 'cube'
step = 15
material = materials.gold()  # structure material

## --- Empty dummy-geometry, will be replaced on run-time by EO trial geometries
geometry = []
struct = structures.struct(step, geometry, material)


## ---------- Setup incident field
field_generator = fields.planewave        # planwave excitation
kwargs = dict(theta = [0, 90])            # 2 polarizations
wavelengths = [700, 1300]                 # 2 wavelengths
efield = fields.efield(field_generator, wavelengths=wavelengths, kwargs=kwargs)


## ---------- environment
n1, n2 = 1.0, 1.0            # in vacuum
dyads = propagators.DyadsQuasistatic123(n1=n1, n2=n2)


## ---------- Simulation initialization
sim = core.simulation(struct, efield, dyads)


#==============================================================================
# setup evolutionary optimization
#==============================================================================
## --- output folder and file-names
results_folder = 'eo_out'
results_suffix = 'test_MO_problem'

## --- structure model and optimization problem
limits_W   = [2, 20]
limits_L   = [2, 20]
limits_pos = [-1, 1]
height = 1
model = CrossAntenna(sim, limits_W, limits_L, limits_pos, height)

problem = ownMultiObjectiveProblem(model)


## --- filename to save results
results_filename = 'emo_test.eo'

## --- size of population
population = 100         # Nr of individuals

## --- stop criteria
max_time = 1800          # seconds
max_iter = 50            # max. iterations
max_nonsuccess = 10      # max. consecutive iterations without improvement

## --- other config
generations = 1          # generations to evolve between status reports
plot_interval = 1
save_all_generations = True

##  Use multi-objective algorithm "nsga2"
import pygmo as pg
algorithm = pg.moead
algorithm_kwargs = dict(weight_generation='low discrepancy',
                        decomposition='weighted')  # optional algo. kwargs
algorithm_kwargs = dict()  # optional algo. kwargs

structure initialization - automatic mesh detection: cube
structure initialization - consistency check: 0/0 dipoles valid

/home/hans/.local/lib/python3.8/site-packages/pyGDM2/tools.py:817: UserWarning: Empty structure. Setting mesh to 'cubic'.
  warnings.warn("Empty structure. Setting mesh to 'cubic'.")
/home/hans/.local/lib/python3.8/site-packages/pyGDM2/structures.py:183: UserWarning: Emtpy structure geometry.
  warnings.warn("Emtpy structure geometry.")

Run the multi-objective optimization

eo_dict = run_eo(problem,
                 population=population,
                 algorithm=algorithm,
                 plot_interval=plot_interval,
                 generations=generations,
                 max_time=max_time, max_iter=max_iter, max_nonsuccess=max_nonsuccess,
                 filename=results_filename)

----------------------------------------------
 Starting new optimization
----------------------------------------------

/home/hans/.local/lib/python3.8/site-packages/pyGDM2/tools.py:835: UserWarning: Mesh not detected, falling back to 'cubic'.
  warnings.warn("Mesh not detected, falling back to 'cubic'.")
/home/hans/.local/lib/python3.8/site-packages/numba/core/dispatcher.py:237: UserWarning: Numba extension module 'numba_scipy' failed to load due to 'ValueError(No function '__pyx_fuse_0pdtr' found in __pyx_capi__ of 'scipy.special.cython_special')'.
  entrypoints.init_all()
/home/hans/.local/lib/python3.8/site-packages/pyGDM2/tools.py:835: UserWarning: Mesh not detected, falling back to 'cubic'.
  warnings.warn("Mesh not detected, falling back to 'cubic'.")

iter #  1, time:   56.0s, progress #  1, f_evals: 200
        - ideal: [-7.3837, -11.642]  (= minimum objective values)
        - nadir: [-5.7113, -0.093708]  (= maximum objective values)
        - 8 fronts, length best front: 5

iter #  2, time:  132.6s, progress #  2, f_evals: 300

/home/hans/.local/lib/python3.8/site-packages/pygmo/plotting/__init__.py:70: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  axes = plt.axes()

        - ideal: [-7.3837, -11.754]  (= minimum objective values)
        - nadir: [-5.6109, -0.093708]  (= maximum objective values)
        - 7 fronts, length best front: 11

iter #  3, time:  179.8s, progress #  3, f_evals: 400
        - ideal: [-7.3837, -11.754]  (= minimum objective values)
        - nadir: [-5.6109, -0.093708]  (= maximum objective values)
        - 6 fronts, length best front: 14

iter #  4, time:  226.3s, progress #  4, f_evals: 500
        - ideal: [-7.3837, -11.754]  (= minimum objective values)
        - nadir: [-5.6109, -0.093708]  (= maximum objective values)
        - 4 fronts, length best front: 23

iter #  5, time:  273.6s, progress #  5, f_evals: 600
        - ideal: [-7.3837, -11.754]  (= minimum objective values)
        - nadir: [-5.6109, -0.093708]  (= maximum objective values)
        - 6 fronts, length best front: 28

iter #  6, time:  320.9s, progress #  6, f_evals: 700
        - ideal: [-7.3837, -11.754]  (= minimum objective values)
        - nadir: [-5.6109, -0.093708]  (= maximum objective values)
        - 6 fronts, length best front: 38

iter #  7, time:  362.8s, progress #  7, f_evals: 800
        - ideal: [-7.3837, -11.754]  (= minimum objective values)
        - nadir: [-5.6109, -0.093708]  (= maximum objective values)
        - 4 fronts, length best front: 57

iter #  8, time:  411.1s, progress #  8, f_evals: 900
        - ideal: [-7.3837, -11.754]  (= minimum objective values)
        - nadir: [-5.6109, -0.093708]  (= maximum objective values)
        - 5 fronts, length best front: 63

iter #  9, time:  455.0s, progress #  9, f_evals: 1000
        - ideal: [-7.3837, -11.907]  (= minimum objective values)
        - nadir: [-1.1464, -0.093708]  (= maximum objective values)
        - 4 fronts, length best front: 78

iter # 10, time:  505.2s, progress # 10, f_evals: 1100
        - ideal: [-7.3837, -11.907]  (= minimum objective values)
        - nadir: [-1.1464, -0.093708]  (= maximum objective values)
        - 3 fronts, length best front: 87

iter # 11, time:  569.0s, progress # 11, f_evals: 1200
        - ideal: [-7.3837, -11.754]  (= minimum objective values)
        - nadir: [-5.6109, -0.093708]  (= maximum objective values)
        - 1 fronts, length best front: 100

iter # 12, time:  612.9s, progress # 11, f_evals: 1300(non-success: 1)
iter # 13, time:  663.8s, progress # 11, f_evals: 1400(non-success: 2)
iter # 14, time:  707.4s, progress # 11, f_evals: 1500(non-success: 3)
iter # 15, time:  757.0s, progress # 11, f_evals: 1600(non-success: 4)
iter # 16, time:  799.4s, progress # 11, f_evals: 1700(non-success: 5)
iter # 17, time:  852.2s, progress # 11, f_evals: 1800(non-success: 6)
iter # 18, time:  907.2s, progress # 11, f_evals: 1900(non-success: 7)
iter # 19, time:  970.3s, progress # 11, f_evals: 2000(non-success: 8)
iter # 20, time: 1013.7s, progress # 11, f_evals: 2100(non-success: 9)
iter # 21, time: 1063.8s, progress # 11, f_evals: 2200(non-success: 10)

 -------- maximum non-successful iterations reached

Reload and analyze multi-objective results

Now we will reload the results. What we actually did is a Pareto-multi-objective optimization, which calculates the so-called Pareto-front. The Pareto-front consists of all solutions that cannot be improved in any of the target functions without worsening at least one of the other objective values.

## --- load additional modules
import numpy as np
import matplotlib.pyplot as plt
import copy


#==============================================================================
# get Pareto-front (set of "Pareto-optima") from optimization
#==============================================================================
filename = "emo_test.eo"

all_paretos, all_sims, all_x = EO.tools.get_pareto_fronts(filename, iteration=-1)

## --- choose "best" Pareto-front (in case several exist):
##     index 0 is the most advanced pareto front --> the optimum solutions
pareto_f = all_paretos[0]
sims_f   = all_sims[0]

## --- plot the pareto front
EO.tools.plot_pareto_2d(pareto_f)

[<matplotlib.lines.Line2D at 0x7fcab1750400>]

Choose the solutions

Pareto-optimizations return a whole set of Pareto-optimum solutions, therefore one has to chose a specific candidate from this set of fittest candidates. We will select 3 geometries to further analyze: The two solutions on the borders of the Pareto-front, as well as one solution from its center.

## --- pick 3 Pareto-optimum solutions
sim0 = sims_f[0]
sim1 = sims_f[13]
sim2 = sims_f[-1]

visu.structure(sim0)
visu.structure(sim1)
visu.structure(sim2)

<matplotlib.collections.PathCollection at 0x7fcab15e2a30>

Analyze the Pareto-optimum solutions

## --- structures
struct0 = copy.deepcopy(sim0.struct)
struct1 = copy.deepcopy(sim1.struct)
struct2 = copy.deepcopy(sim2.struct)

## --- incident field
field_generator = fields.planewave        # planwave excitation
wavelengths = np.arange(400, 1610, 20)  # spectrum
kwargs = dict(theta = [0.0, 90.0])          # 0/90 deg polarizations
efield = fields.efield(field_generator, wavelengths=wavelengths, kwargs=kwargs)

## ---------- environment
dyads = propagators.DyadsQuasistatic123(n1=1, n2=1)


## --- spectrum simulations
sim_spec0 = core.simulation(struct0, efield, dyads)
sim_spec1 = core.simulation(struct1, efield, dyads)
sim_spec2 = core.simulation(struct2, efield, dyads)

## --- run the simulations
sim_spec0.scatter(verbose=False)
sim_spec1.scatter(verbose=False)
sim_spec2.scatter(verbose=False)

1

Calculate scattering spectra

We now calculate the scattering efficiency spectra for the three selected structures:

## -- structure #0
wl, _ext0 = tools.calculate_spectrum(sim_spec0, 0, linear.extinct)
wl, _ext90 = tools.calculate_spectrum(sim_spec0, 1, linear.extinct)
geom_cs = tools.get_geometric_cross_section(sim_spec0)
Qs0_X = _ext0.T[1] / geom_cs
Qs0_Y = _ext90.T[1] / geom_cs

## -- structure #1
wl, _ext0 = tools.calculate_spectrum(sim_spec1, 0, linear.extinct)
wl, _ext90 = tools.calculate_spectrum(sim_spec1, 1, linear.extinct)
geom_cs = tools.get_geometric_cross_section(sim_spec1)
Qs1_X = _ext0.T[1] / geom_cs
Qs1_Y = _ext90.T[1] / geom_cs

## -- structure #2
wl, _ext0 = tools.calculate_spectrum(sim_spec2, 0, linear.extinct)
wl, _ext90 = tools.calculate_spectrum(sim_spec2, 1, linear.extinct)
geom_cs = tools.get_geometric_cross_section(sim_spec2)
Qs2_X = _ext0.T[1] / geom_cs
Qs2_Y = _ext90.T[1] / geom_cs

Plot the spectra

plt.figure(figsize=(10,15))

## --- structure 0
## spectra
plt.subplot2grid((3,5), (0,0), colspan=3)
plt.title("structure '0'")
plt.plot(wl, Qs0_X, label="pol. X")
plt.plot(wl, Qs0_Y, label="pol. Y")
plt.ylabel("Q_scat")
plt.legend(loc='best', fontsize=10)

## geometry
plt.subplot2grid((3,5), (0,3), colspan=2, aspect="equal")
visu.structure(sim_spec0, show=False)


## --- structure 1
## spectra
plt.subplot2grid((3,5), (1,0), colspan=3)
plt.title("structure '1'")
plt.plot(wl, Qs1_X, label="pol. X")
plt.plot(wl, Qs1_Y, label="pol. Y")
plt.ylabel("Q_scat")

## geometry
plt.subplot2grid((3,5), (1,3), colspan=2, aspect="equal")
visu.structure(sim_spec1, show=False)


## --- structure 2
## spectra
plt.subplot2grid((3,5), (2,0), colspan=3)
plt.title("structure '2'")
plt.plot(wl, Qs2_X, label="pol. X")
plt.plot(wl, Qs2_Y, label="pol. Y")
plt.xlabel("wavelength (nm)")
plt.ylabel("Q_scat")

## geometry
plt.subplot2grid((3,5), (2,3), colspan=2, aspect="equal")
visu.structure(sim_spec2, show=False)


plt.show()

The structures on the “border” of the Pareto-front are equivalent to single-objective optimizations for either of the two target functions:

• Structure “0” corresponds to a maximization of \(Q_{\text{scat}}\) for 1300nm and Y polarization

• Structure “2” corresponds to a maximization of \(Q_{\text{scat}}\) for 700nm and X polarization

Structure “1” is the most interesting result: We take the structure on the Pareto-Front with most similar \(Q_{\text{scat}}\) for the two target configurations, and indeed, the optimization found a double-resonant structure, with addressable resonances via a rotation of the incident light polarization.

Note: The size-limit of \(20\) meshpoints in each dimension seems to be slightly too small to obtain a resonance at \(\lambda=1300\,\)nm in the double-resonant case.