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:
In [1]:
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:
In [2]:
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., self.min_dist, self.min_dist, 0.]
def min_values_float(self):
return [-2*self.radius, -2*self.radius, self.min_dist, self.min_dist, 0.]
def max_values_float(self):
return [self.radius, self.radius, self.max_dist, self.max_dist, 90.]
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.])
naz = np.array([0., 0., 1.])
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:
In [3]:
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:
In [4]:
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
--------------------------------------------------
In [5]:
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()
Initializing Optimizer_pymoo
Selecting algorithm: MixedVariableGA (MixedVariableGA)
Problem:
--------
vectorize: True
Algorithm:
----------
type: MixedVariableGA
pop_size: 40
seed: 42
Termination:
------------
n_eval: 1000
Here tol is an IPOPT parameter that defines the convergence tolerance. Now we are finally ready - let’s solve the problem!
In [6]:
results = solver.solve()
solver.finalize(results)
=================================================
n_gen | n_eval | f_avg | f_min
=================================================
1 | 40 | -3.100000E+00 | -1.800000E+01
2 | 80 | -5.425000E+00 | -2.400000E+01
3 | 120 | -1.132500E+01 | -5.100000E+01
4 | 160 | -2.195000E+01 | -6.800000E+01
5 | 200 | -3.755000E+01 | -1.120000E+02
6 | 240 | -6.500000E+01 | -1.530000E+02
7 | 280 | -9.782500E+01 | -1.710000E+02
8 | 320 | -1.257500E+02 | -2.590000E+02
9 | 360 | -1.531750E+02 | -2.590000E+02
10 | 400 | -1.788000E+02 | -2.680000E+02
11 | 440 | -2.088500E+02 | -2.740000E+02
12 | 480 | -2.367500E+02 | -2.750000E+02
13 | 520 | -2.557750E+02 | -2.750000E+02
14 | 560 | -2.630750E+02 | -2.800000E+02
15 | 600 | -2.696000E+02 | -3.120000E+02
16 | 640 | -2.793250E+02 | -3.120000E+02
17 | 680 | -2.868250E+02 | -3.120000E+02
18 | 720 | -2.933250E+02 | -3.120000E+02
19 | 760 | -2.991250E+02 | -3.120000E+02
20 | 800 | -3.042500E+02 | -3.150000E+02
21 | 840 | -3.078500E+02 | -3.150000E+02
22 | 880 | -3.110000E+02 | -3.150000E+02
23 | 920 | -3.125500E+02 | -3.150000E+02
24 | 960 | -3.135250E+02 | -3.160000E+02
25 | 1000 | -3.140250E+02 | -3.160000E+02
Optimizer_pymoo: Optimization run finished
Success: True
Best maxN = 316.0
These are the results:
In [7]:
print(results)
fig = problem.get_fig()
plt.show()
plt.close(fig)
Results problem 'grid_problem':
--------------------------------
Integer variables:
0: nx = 721
1: ny = 337
--------------------------------
Float variables:
0: x0 = -9.852999e+00
1: y0 = -8.461444e+00
2: dx = 5.024268e-01
3: dy = 5.000655e-01
4: alpha = 6.319982e+00
--------------------------------
Objectives:
0: maxN = 3.160000e+02
--------------------------------
Success: True
--------------------------------
Clearly the circle is fully filled with points on a regular grid.