Electrostatic potential minimization

The main purpose of iwopy is to be helpful when attacking more complicated optimization tasks than the minimization of simple analytical functions. As an example, we consider the minimization of an inverse distance type potential for a fully coupled system of \(N\) particles in two dimensions. We can imagine such a system to be composed of \(N\) individual point charges, each of them carrying the same electric unit charge. The potential that we want to minimize is then represented by $\( V = \sum_{i\neq j} \frac 1 {|\mathrm r_i - \mathrm r_j|} \ , \quad \text{where} \quad \mathrm r_i = \left( \begin{array}{c} x_i \\ y_i \end{array} \right) \)\( denotes the two dimensional position vector of charge \)i\( and the sum is over all unequal index pairs. This potential favours large distances between charges, hence we need to constrain them to a certain area for a meaningful solution. We confine them to a circle of radius \)R\( by imposing \)N\( constraints, \)\( |\mathrm r_i| \leq R \ , \quad \text{for all} \quad i = 0 \ldots N-1 \ . \)$ These are the required imports for this example:

import matplotlib.pyplot as plt
import numpy as np

from iwopy import Constraint, Objective, Problem, create_optimizer

We start by creating a specific class that describes the variables of our problem:

class ChargesProblem(Problem):
    def __init__(self, xy_init, radius):
        super().__init__(name="charges_problem")
        self.xy_init = xy_init
        self.n_charges = len(xy_init)
        self.radius = radius

    def var_names_float(self):
        """Defines the variable names"""
        vnames = []
        for i in range(self.n_charges):
            vnames += [f"x{i}", f"y{i}"]
        return vnames

    def initial_values_float(self):
        """Returns inital values, as given to constructor"""
        return self.xy_init.reshape(2 * self.n_charges)

    def min_values_float(self):
        """Minimal values for a square of size 2*radius"""
        return np.full(2 * self.n_charges, -self.radius)

    def max_values_float(self):
        """Maximal values for a square of size 2*radius"""
        return np.full(2 * self.n_charges, self.radius)

    def apply_individual(self, vars_int, vars_float):
        """Returns (x, y) variables for each charge"""
        return vars_float.reshape(self.n_charges, 2)

    def apply_population(self, vars_int, vars_float):
        """Returns (x, y) variables for each charge per individual"""
        n_pop = len(vars_float)
        return vars_float.reshape(n_pop, self.n_charges, 2)

    def get_fig(self, xy):
        """Get a figure showing the charges locations"""
        fig, ax = plt.subplots(figsize=(6, 6))
        ax.scatter(xy[:, 0], xy[:, 1], color="orange")
        ax.add_patch(plt.Circle((0, 0), self.radius, color="darkred", fill=False))
        ax.set_aspect("equal", adjustable="box")
        ax.set_xlabel("x")
        ax.set_ylabel("y")
        ax.set_title(f"N = {self.n_charges}")
        return fig

There are \(2 N\) variables for this problem, and we order them as \((x_0, y_0, \ldots, x_{N-1}, y_{N-1})\). This is convenient, since the numpy reshaping operation vars_float.reshape(n_charges, 2) then represents an array of \((x_i, y_i)\) type arrays, which is handy for calculations.

Notice the two functions apply_individual and apply_population. They are being called during optimization whenever new variables are being set by the optimizer, and their purpose is to update the data in the problem class accordingly (excluding objectives and constraints, which we will deal with shortly).

The difference between the two functions is that apply_individual evaluates a single vector of new problem variables, whereas apply_population handles a full population of such vectors. The implementation of the latter function is not strictly required (alternatively a loop over individuals of the population will be evaluated), but it provides a chance for a vast speed-up by vectorized evaluation. This is particularly beneficial for genetic algorithms or other easily vectorizable approaches (in fact also the gradient calculation via iwopy discretizations makes use of this vectorization). Of course the chosen optimizer must be able to support vectorized evaluation for this to work. Both functions return any kind of evaluation data which will be forwarded to the objective and constraint evaluation.

So far we have only defined the problem variables. Objectives and constraints can be added to the problem via the add_objective and add_constraint functions, respectively. First, we implement the inverse distance potential as our objective:

class MinPotential(Objective):
    def __init__(self, problem):
        """Define the same variable names and ordering as in the problem"""
        super().__init__(problem, "potential", vnames_float=problem.var_names_float())
        self.n_charges = problem.n_charges

    def n_components(self):
        """The potential is a scalar function, hence one component"""
        return 1

    def maximize(self):
        """Indicates that the single component is to be minimized"""
        return [False]

    def calc_individual(self, vars_int, vars_float, problem_results, components=None):
        """This evaluates the potential. See problem.apply_individual
        for problem results"""
        xy = problem_results
        value = 0.0
        for i in range(1, self.n_charges):
            dist = np.maximum(np.linalg.norm(xy[i - 1, None] - xy[i:], axis=-1), 1e-10)
            value += 2 * np.sum(1 / dist)
        return value

    def calc_population(self, vars_int, vars_float, problem_results, components=None):
        """This evaluates the potential in vectorized manner.
        See problem.apply_population for problem results"""
        xy = problem_results
        n_pop = len(xy)
        value = np.zeros((n_pop, 1))
        for i in range(1, self.n_charges):
            dist = np.maximum(
                np.linalg.norm(xy[:, i - 1, None] - xy[:, i:], axis=-1), 1e-10
            )
            value[:, 0] += 2 * np.sum(1 / dist, axis=1)
        return value

Notice that the functions calc_individual and its (again optional) vectorized correspondent calc_population do not directly make use of the variables vector vars_float, which they could (and it would be perfectly fine, even intended), but instead of the problem_results provided by the problem functions apply_individual and apply_population.

Next, we implement the \(N\) radial constraints:

class MaxRadius(Constraint):
    def __init__(self, problem, tol=1e-2):
        """Define the same variable names and ordering as in the problem"""
        super().__init__(
            problem, "radius", vnames_float=problem.var_names_float(), tol=tol
        )
        self.n_charges = problem.n_charges
        self.radius = problem.radius

    def n_components(self):
        """One constraint per charge"""
        return self.n_charges

    def vardeps_float(self):
        """Boolean array that defines which component depends
        on which variable (optional, default is on all)"""
        deps = np.zeros((self.n_components(), self.n_charges, 2), dtype=bool)
        np.fill_diagonal(deps[..., 0], True)
        np.fill_diagonal(deps[..., 1], True)
        return deps.reshape(self.n_components(), 2 * self.n_charges)

    def calc_individual(self, vars_int, vars_float, problem_results, components=None):
        """This evaluates the constraints, negative values are valid.
        See problem.apply_individual for problem results"""
        components = np.s_[:] if components is None else components
        xy = problem_results[components]
        r = np.linalg.norm(xy, axis=-1)
        return r - self.radius

    def calc_population(self, vars_int, vars_float, problem_results, components=None):
        """This evaluates the constraints in vectorized manner,
        negative values are valid. See problem.apply_population for
        problem results"""
        components = np.s_[:] if components is None else components
        xy = problem_results[:, components]
        r = np.linalg.norm(xy, axis=-1)
        return r - self.radius

Note that by default, negative constraint values represent validity. This can be changed by overloading the function get_bounds of the base class, see API, but is not recommended.

The vardeps_float mask defines structural zeros in the component-by-variable Jacobian. iwopy uses it during analytical and numerical derivative assembly for all gradient-based optimizers. PyGMO/IPOPT and SciPy trust-constr also use sparse constraint Jacobians at the native boundary. The genetic algorithm used here does not request derivatives, so the mask is included for completeness.

Now let’s add the objective and the constraints and initialize the problem:

N = 20
R = 5

# generate random initial coordinates of the N charges:
xy = np.random.uniform(-R / 10.0, R / 10.0, (N, 2))

problem = ChargesProblem(xy, R)
problem.add_objective(MinPotential(problem))
problem.add_constraint(MaxRadius(problem))
problem.initialize()

fig = problem.get_fig(xy)
plt.show()

Problem 'charges_problem' (ChargesProblem): Initializing
--------------------------------------------------------
../_images/825b2fd109173cd2d83398b312b3cb8bc97a76652ac03518e4f967f8aa25b57c.png

Finally, we can now setup the genetic algorithm GA from pymoo and solve the problem, here in vectorized form (flag vectorize=True):

solver = create_optimizer(
    family="pymoo",  # Backend family.
    optimizer="GA",  # Genetic algorithm.
    problem=problem,
    vectorize=True,  # Evaluate each population together.
    pop_size=100,  # Population size; a production starting point is 100.
    seed=1,  # Reproducible search.
    termination="default",  # Native default convergence criteria.
    max_iterations=300,  # Notebook limit; start around 500 in production.
    ftol=1e-6,  # Objective convergence tolerance.
    xtol=1e-6,  # Variable convergence tolerance.
)
solver.initialize()
solver.print_info()

results = solver.solve(verbosity=0)
solver.finalize(results)

Loading pymoo

pymoo successfully loaded
Initializing Optimizer_pymoo
Selecting sampling: float_random (FloatRandomSampling)
Selecting algorithm: GA (GA)
Selecting termination: default (DefaultSingleObjectiveTermination)

Problem:
--------
  name         : charges_problem
  n_vars_int   : 0
  n_vars_float : 40
  n_objectives : 1
  n_obj_cmptns : 1
  n_constraints: 1
  n_con_cmptns : 20

Algorithm:
----------
  pop_size: 100

Termination:
------------
  type: default
  n_max_gen: 300
  ftol: 1e-06
  xtol: 1e-06

Optimizer_pymoo: Optimization run finished
  Success: True
  Best potential = 78.63202076960529

print(results)

xy = results.problem_results
fig = problem.get_fig(xy)
plt.show()

Results problem 'charges_problem':
-----------------------------------
  Float variables:
    0: x0  = 1.143531e+00
    1: y0  = -4.867063e+00
    2: x1  = 2.299241e+00
    3: y1  = 4.439655e+00
    4: x2  = 4.974626e+00
    5: y2  = -5.019755e-01
    6: x3  = -4.981757e+00
    7: y3  = -4.226352e-01
    8: x4  = 4.216121e-01
    9: y4  = 4.981225e+00
    10: x5  = 4.560218e+00
    11: y5  = -2.047932e+00
    12: x6  = -3.362242e+00
    13: y6  = 3.691631e+00
    14: x7  = -2.509514e+00
    15: y7  = 1.169439e+00
    16: x8  = -1.772239e+00
    17: y8  = 4.674721e+00
    18: x9  = 1.082625e+00
    19: y9  = -2.407658e+00
    20: x10 = -4.743970e+00
    21: y10 = 1.578952e+00
    22: x11 = -2.661491e+00
    23: y11 = -4.232701e+00
    24: x12 = 2.041850e+00
    25: y12 = 4.783230e-01
    26: x13 = 3.166017e+00
    27: y13 = -3.869272e+00
    28: x14 = -8.736933e-01
    29: y14 = -4.921881e+00
    30: x15 = -1.998249e+00
    31: y15 = -1.832655e+00
    32: x16 = 1.243952e-01
    33: y16 = 2.600018e+00
    34: x17 = 4.131653e+00
    35: y17 = 2.812240e+00
    36: x18 = 4.865998e+00
    37: y18 = 1.143517e+00
    38: x19 = -4.348955e+00
    39: y19 = -2.459572e+00
-----------------------------------
  Objectives:
    0: potential = 7.863202e+01
-----------------------------------
  Constraints:
    0: radius_0  = -4.037685e-04
    1: radius_1  = -2.957720e-04
    2: radius_2  = -1.115928e-04
    3: radius_3  = -3.473239e-04
    4: radius_4  = -9.644986e-04
    5: radius_5  = -1.038249e-03
    6: radius_6  = -6.723411e-03
    7: radius_7  = -2.231381e+00
    8: radius_8  = -6.151979e-04
    9: radius_9  = -2.360134e+00
    10: radius_10 = -1.657219e-04
    11: radius_11 = -7.120214e-05
    12: radius_12 = -2.902872e+00
    13: radius_13 = -5.072941e-04
    14: radius_14 = -1.174803e-03
    15: radius_15 = -2.288612e+00
    16: radius_16 = -2.397008e+00
    17: radius_17 = -2.075239e-03
    18: radius_18 = -1.443679e-03
    19: radius_19 = -3.710411e-03
-----------------------------------
  Success: True
-----------------------------------
../_images/1d34f7c628e7f171b6b3079699afa9b079bff2a5f024017ffef68ed5df0a85c1.png

Note that this problem has many equivalent solutions, and the displayed result is the best solution that the genetic algorithm was able to find within its convergence criteria. With heuristic methods there is generally no guarantee or proof that a global optimum has been found. However, the objective function has clearly been minimized, and therefore the “design” was improved which is often the goal in engineering tasks.

For comparison, let’s also solve the problem using iwopy’s native SLSQP optimizer, now for 50 charges. SLSQP is a gradient-based constrained optimizer backed by SciPy. The LocalFD wrapper supplies finite-difference derivatives and evaluates the perturbations through iwopy’s vectorized population interface. The solve is capped at 100 iterations to keep the example runtime bounded.

from iwopy import LocalFD

N = 50
R = 5

np.random.seed(42)
xy = np.random.uniform(-R / 10.0, R / 10.0, (N, 2))

problem = ChargesProblem(xy, R)
problem.add_objective(MinPotential(problem))
problem.add_constraint(MaxRadius(problem))

gproblem = LocalFD(problem, deltas=1e-2, fd_order=1)
gproblem.initialize()

solver = create_optimizer(
    family="iwopy",  # iwopy optimizer implementation.
    optimizer="SLSQP",  # Scaled local gradient search.
    problem=gproblem,
    tol=1e-6,  # Native solver tolerance.
    max_iterations=100,  # Notebook limit; start around 500 in production.
    vectorize=True,  # Batch finite-difference gradient evaluations.
)
solver.initialize()

results = solver.solve(verbosity=0)
solver.finalize(results)
print(results)

Problem 'charges_problem' (ChargesProblem): Initializing
--------------------------------------------------------
Problem 'charges_problem_fd' (LocalFD): Initializing
----------------------------------------------------

SLSQP: Optimization run finished
  Success: False
  Best potential = 594.0014353404423
Results problem 'charges_problem_fd':
--------------------------------------
  Float variables:
    0: x0  = 1.364443e+00
    1: y0  = 4.801353e+00
    2: x1  = 1.880671e+00
    3: y1  = 5.859686e-01
    4: x2  = 5.624808e-01
    5: y2  = -4.941567e+00
    6: x3  = 1.284452e-01
    7: y3  = -3.613355e+00
    8: x4  = -4.289386e-02
    9: y4  = 2.276008e+00
    10: x5  = -2.921010e+00
    11: y5  = 4.070232e+00
    12: x6  = -7.347943e-01
    13: y6  = -4.925829e+00
    14: x7  = 2.383899e+00
    15: y7  = 2.834112e+00
    16: x8  = 4.850614e+00
    17: y8  = -1.068969e+00
    18: x9  = -3.615992e+00
    19: y9  = 3.464299e+00
    20: x10 = -3.365053e+00
    21: y10 = -3.602091e+00
    22: x11 = -1.681125e+00
    23: y11 = -4.523549e+00
    24: x12 = 7.071707e-02
    25: y12 = 4.916371e+00
    26: x13 = 3.752410e+00
    27: y13 = 2.182559e-01
    28: x14 = 2.972416e+00
    29: y14 = -2.623446e+00
    30: x15 = -2.234937e+00
    31: y15 = 2.792316e+00
    32: x16 = -2.797864e+00
    33: y16 = -8.248501e-01
    34: x17 = 4.872613e+00
    35: y17 = 1.072300e+00
    36: x18 = -3.988664e+00
    37: y18 = 5.714319e-01
    38: x19 = -1.614136e+00
    39: y19 = -3.268282e+00
    40: x20 = -4.910257e+00
    41: y20 = 7.186934e-01
    42: x21 = -4.502307e+00
    43: y21 = -2.142219e+00
    44: x22 = 3.137615e+00
    45: y22 = 1.570241e+00
    46: x23 = 4.445302e+00
    47: y23 = 2.040674e+00
    48: x24 = -4.509826e+00
    49: y24 = -1.266617e+00
    50: x25 = 4.584564e+00
    51: y25 = -1.996534e+00
    52: x26 = 3.230365e+00
    53: y26 = 3.779994e+00
    54: x27 = 3.971590e+00
    55: y27 = 3.038240e+00
    56: x28 = -2.589957e+00
    57: y28 = -4.283697e+00
    58: x29 = -4.644107e+00
    59: y29 = 1.802595e+00
    60: x30 = -1.005920e+00
    61: y30 = 3.762524e+00
    62: x31 = 3.440034e+00
    63: y31 = -3.638265e+00
    64: x32 = -2.069594e+00
    65: y32 = 4.530066e+00
    66: x33 = -1.028289e+00
    67: y33 = 4.881546e+00
    68: x34 = 2.321657e+00
    69: y34 = 4.432296e+00
    70: x35 = -5.729551e-01
    71: y35 = -1.390764e+00
    72: x36 = -4.923851e+00
    73: y36 = -4.131104e-01
    74: x37 = -1.906142e+00
    75: y37 = 6.482256e-01
    76: x38 = 2.661702e+00
    77: y38 = -4.221638e+00
    78: x39 = -2.800183e+00
    79: y39 = -2.307495e+00
    80: x40 = 4.105479e+00
    81: y40 = -2.779350e+00
    82: x41 = 1.654198e+00
    83: y41 = -4.716092e+00
    84: x42 = 1.531609e+00
    85: y42 = -3.638168e+00
    86: x43 = 4.961507e+00
    87: y43 = -8.485057e-02
    88: x44 = 2.781662e+00
    89: y44 = -1.249452e+00
    90: x45 = -3.704838e+00
    91: y45 = 2.243721e+00
    92: x46 = 9.009654e-01
    93: y46 = 3.418176e+00
    94: x47 = 7.387893e-01
    95: y47 = -1.893857e+00
    96: x48 = -1.009258e-01
    97: y48 = 7.186303e-01
    98: x49 = -4.066497e+00
    99: y49 = -2.915331e+00
--------------------------------------
  Objectives:
    0: potential = 5.940014e+02
--------------------------------------
  Constraints:
    0: radius_0  = -8.537778e-03
    1: radius_1  = -3.030156e+00
    2: radius_2  = -2.652345e-02
    3: radius_3  = -1.384363e+00
    4: radius_4  = -2.723587e+00
    5: radius_5  = 9.899430e-03
    6: radius_6  = -1.966681e-02
    7: radius_7  = -1.296601e+00
    8: radius_8  = -3.299374e-02
    9: radius_9  = 7.671014e-03
    10: radius_10 = -7.063502e-02
    11: radius_11 = -1.741659e-01
    12: radius_12 = -8.312001e-02
    13: radius_13 = -1.241248e+00
    14: radius_14 = -1.035441e+00
    15: radius_15 = -1.423414e+00
    16: radius_16 = -2.083080e+00
    17: radius_17 = -1.079290e-02
    18: radius_18 = -9.706108e-01
    19: radius_19 = -1.354852e+00
    20: radius_20 = -3.742613e-02
    21: radius_21 = -1.403286e-02
    22: radius_22 = -1.491399e+00
    23: radius_23 = -1.086748e-01
    24: radius_24 = -3.156802e-01
    25: radius_25 = 4.378355e-04
    26: radius_26 = -2.771618e-02
    27: radius_27 = 4.428824e-04
    28: radius_28 = 5.789933e-03
    29: radius_29 = -1.832554e-02
    30: radius_30 = -1.105329e+00
    31: radius_31 = 7.076165e-03
    32: radius_32 = -1.956599e-02
    33: radius_33 = -1.132547e-02
    34: radius_34 = 3.532839e-03
    35: radius_35 = -3.495839e+00
    36: radius_36 = -5.884897e-02
    37: radius_37 = -2.986651e+00
    38: radius_38 = -9.320261e-03
    39: radius_39 = -1.371562e+00
    40: radius_40 = -4.220400e-02
    41: radius_41 = -2.210597e-03
    42: radius_42 = -1.052584e+00
    43: radius_43 = -3.776727e-02
    44: radius_44 = -1.950611e+00
    45: radius_45 = -6.687053e-01
    46: radius_46 = -1.465079e+00
    47: radius_47 = -2.967144e+00
    48: radius_48 = -4.274317e+00
    49: radius_49 = 3.554396e-03
--------------------------------------
  Success: False
--------------------------------------

xy = results.problem_results
fig = problem.get_fig(xy)
plt.show()
../_images/2e40006da9ad436f5984b4e773d43d9d3af4c6e685abed3191fe014ab9a4b8f1.png