Mixed problems

The pymoo interface supports the solution of mixed integer/float problems via the genetic algorithm MixedVariableGA. Here we demonstrate this by maximizing the number if points that fit inside a circle and lie on a regular grid.

These are the imports:

import numpy as np
import matplotlib.pyplot as plt

from iwopy import Problem, Objective
from iwopy.interfaces.pymoo import Optimizer_pymoo

The problem describes points on a regular grid of size nx * ny. Points that lie outside of the given radius will be marked by a False value in a boolean array called valid:

class GridProblem(Problem):
    def __init__(self, n_row_max, radius, min_dist, ctol=1e-3):
        super().__init__(name="grid_problem")

        self.n_row_max = n_row_max
        self.radius = float(radius)
        self.min_dist = float(min_dist)
        self.max_dist = 2 * radius

        self.xy = None
        self.valid = None

    def initialize(self, verbosity=1):
        super().initialize(verbosity)
        self.apply_individual(self.initial_values_int(), self.initial_values_float())

    def var_names_int(self):
        return ["nx", "ny"]

    def initial_values_int(self):
        return [2, 2]

    def min_values_int(self):
        return [1, 1]

    def max_values_int(self):
        return [self.n_row_max, self.n_row_max]

    def var_names_float(self):
        return ["x0", "y0", "dx", "dy", "alpha"]

    def initial_values_float(self):
        return [0.0, 0.0, self.min_dist, self.min_dist, 0.0]

    def min_values_float(self):
        return [-2 * self.radius, -2 * self.radius, self.min_dist, self.min_dist, 0.0]

    def max_values_float(self):
        return [self.radius, self.radius, self.max_dist, self.max_dist, 90.0]

    def apply_individual(self, vars_int, vars_float):
        """Create points on a regular grid, and evaluate their validity"""

        nx, ny = vars_int
        x0, y0, dx, dy, alpha = vars_float

        a = np.deg2rad(alpha)
        nax = np.array([np.cos(a), np.sin(a), 0.0])
        naz = np.array([0.0, 0.0, 1.0])
        nay = np.cross(naz, nax)

        self.xy = np.zeros((nx, ny, 2))
        self.xy[:] = np.array([x0, y0])[None, None, :]
        self.xy[:] += np.arange(nx)[:, None, None] * dx * nax[None, None, :2]
        self.xy[:] += np.arange(ny)[None, :, None] * dy * nay[None, None, :2]

        self.valid = np.linalg.norm(self.xy, axis=-1) <= self.radius

        return self.xy, self.valid

    def apply_population(self, vars_int, vars_float):
        """Create points on a regular grid, and evaluate their validity"""

        n_pop = vars_int.shape[0]
        nx = vars_int[:, 0]
        ny = vars_int[:, 1]
        x0 = vars_float[:, 0]
        y0 = vars_float[:, 1]
        dx = vars_float[:, 2]
        dy = vars_float[:, 3]
        alpha = vars_float[:, 4]

        a = np.deg2rad(alpha)
        nax = np.stack([np.cos(a), np.sin(a), np.zeros(a.shape)], axis=-1)
        naz = np.zeros_like(nax)
        naz[:, 2] = 1
        nay = np.cross(naz, nax)

        mx = np.max(nx)
        my = np.max(ny)
        self.xy = np.full((n_pop, mx, my, 2), -2 * self.radius)
        for i in range(n_pop):
            self.xy[i, : nx[i], : ny[i]] = np.array([x0[i], y0[i]])[None, None, :]
            self.xy[i, : nx[i], : ny[i]] += (
                np.arange(nx[i])[:, None, None]
                * dx[i, None, None, None]
                * nax[i, None, None, :2]
            )
            self.xy[i, : nx[i], : ny[i]] += (
                np.arange(ny[i])[None, :, None]
                * dy[i, None, None, None]
                * nay[i, None, None, :2]
            )

        self.valid = np.linalg.norm(self.xy, axis=-1) <= self.radius

        return self.xy, self.valid

    def get_fig(self, xy=None, valid=None):
        if xy is None:
            xy = self.xy
        if valid is None:
            valid = self.valid

        nx, ny = xy.shape[:2]
        xy = xy.reshape(nx * ny, 2)[valid.reshape(nx * ny)]

        fig, ax = plt.subplots()
        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 = {len(xy)}, min_dist = {self.min_dist}")

        return fig

Notice that the calculate_individual and calculate_population functions return the current grid points and the validity array. This problem has two integer variables nx and ny, determining the number of grid points in x and y directions, respectively, and 5 float type variables.

Our objective is the maximization of the total number of points within the circle:

class MaxN(Objective):
    def __init__(self, problem):
        super().__init__(
            problem,
            "maxN",
            vnames_int=problem.var_names_int(),
            vnames_float=problem.var_names_float(),
        )

    def n_components(self):
        return 1

    def maximize(self):
        return [True]

    def calc_individual(self, vars_int, vars_float, problem_results, cmpnts=None):
        xy, valid = problem_results
        return np.sum(valid)

    def calc_population(self, vars_int, vars_float, problem_results, cmpnts=None):
        xy, valid = problem_results
        return np.sum(valid, axis=(1, 2))[:, None]

This objective makes use of the returned validity array. We can now solve this with the appropriate algorithm from pymoo:

problem = GridProblem(n_row_max=1000, radius=5, min_dist=0.5)
problem.add_objective(MaxN(problem))
problem.initialize()
Problem 'grid_problem' (GridProblem): Initializing
--------------------------------------------------
  n_vars_int  : 2
  n_vars_float: 5
--------------------------------------------------
  n_objectives: 1
  n_obj_cmptns: 1
--------------------------------------------------
  n_constraints: 0
  n_con_cmptns: 0
--------------------------------------------------
solver = Optimizer_pymoo(
    problem,
    problem_pars=dict(
        vectorize=True,
    ),
    algo_pars=dict(
        type="MixedVariableGA",
        pop_size=40,
        seed=42,
    ),
    setup_pars=dict(),
    term_pars=("n_eval", 1000),
)
solver.initialize()
solver.print_info()
Loading pymoo
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
File ~/work/iwopy/iwopy/iwopy/utils/load.py:27, in import_module(name, package, hint)
     26 try:
---> 27     return importlib.import_module(name, package)
     28 except ModuleNotFoundError:

File ~/work/_temp/uv-python-dir/cpython-3.10.19-linux-x86_64-gnu/lib/python3.10/importlib/__init__.py:126, in import_module(name, package)
    125         level += 1
--> 126 return _bootstrap._gcd_import(name[level:], package, level)

File <frozen importlib._bootstrap>:1050, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:1027, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:992, in _find_and_load_unlocked(name, import_)

File <frozen importlib._bootstrap>:241, in _call_with_frames_removed(f, *args, **kwds)

File <frozen importlib._bootstrap>:1050, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:1027, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:992, in _find_and_load_unlocked(name, import_)

File <frozen importlib._bootstrap>:241, in _call_with_frames_removed(f, *args, **kwds)

File <frozen importlib._bootstrap>:1050, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:1027, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:1004, in _find_and_load_unlocked(name, import_)

ModuleNotFoundError: No module named 'pymoo'

During handling of the above exception, another exception occurred:

ModuleNotFoundError                       Traceback (most recent call last)
Cell In[5], line 14
      1 solver = Optimizer_pymoo(
      2     problem,
      3     problem_pars=dict(
   (...)
     12     term_pars=("n_eval", 1000),
     13 )
---> 14 solver.initialize()
     15 solver.print_info()

File ~/work/iwopy/iwopy/iwopy/interfaces/pymoo/optimizer.py:159, in Optimizer_pymoo.initialize(self, verbosity)
    149 """
    150 Initialize the object.
    151 
   (...)
    156 
    157 """
    158 if self.problem.n_objectives <= 1:
--> 159     self.pymoo_problem = SingleObjProblemTemplate.get_class()(
    160         self.problem, **self.problem_pars
    161     )
    162 else:
    163     self.pymoo_problem = MultiObjProblemTemplate.get_class()(
    164         self.problem, **self.problem_pars
    165     )

File ~/work/iwopy/iwopy/iwopy/interfaces/pymoo/problem.py:278, in SingleObjProblemTemplate.get_class(cls)
    273 @classmethod
    274 def get_class(cls):
    275     """
    276     Creates the class, dynamically derived from pymoo.Problem
    277     """
--> 278     imports.load()
    279     attrb = {v: d for v, d in cls.__dict__.items()}
    280     init0 = cls.__init__

File ~/work/iwopy/iwopy/iwopy/interfaces/pymoo/imports.py:51, in load(verbosity)
     48 if verbosity:
     49     print("Loading pymoo")
---> 51 Callback = import_module(
     52     "pymoo.core.callback", hint="pip install pymoo"
     53 ).Callback
     54 Problem = import_module("pymoo.core.problem", hint="pip install pymoo").Problem
     55 Real = import_module("pymoo.core.variable", hint="pip install pymoo").Real

File ~/work/iwopy/iwopy/iwopy/utils/load.py:31, in import_module(name, package, hint)
     29 mdl = name if package is None else f"{package}.{name}"
     30 hint = hint if hint is not None else f"pip install {name}"
---> 31 raise ModuleNotFoundError(f"Module '{mdl}' not found, maybe try '{hint}'")

ModuleNotFoundError: Module 'pymoo.core.callback' not found, maybe try 'pip install pymoo'

Here tol is an IPOPT parameter that defines the convergence tolerance. Now we are finally ready - let’s solve the problem!

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

These are the results:

print(results)

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

Clearly the circle is fully filled with points on a regular grid.