{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Multi objective geometric chain" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we demonstrate how to solve a simple multi-objective problem via `pymoo`'s [NSGA2](https://pymoo.org/algorithms/moo/nsga2.html) algorithm within the `iwopy` framework.\n", "\n", "The problem that we want to attack here is based on a chain of geometrically touching circles in two dimensions, for example representing a chain of marbles on a string. The marbles may have different radii, and each of them is allowed to touch two neighbours (except the outermost marbles, which have only one neighbour).\n", "\n", "Our two objectives for this chain of blobs are:\n", "\n", "- Maximize the extension of the overall chain in x-direction,\n", "- Maximize the extension of the overall chain in y-direction.\n", "\n", "Obviously these objectives are contradictory, and hence we can expect a nice and clean Pareto front when looking at solutions.\n", "\n", "Here are the required imports:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from iwopy import Problem, Constraint, Objective\n", "from iwopy.interfaces.pymoo import Optimizer_pymoo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we create a class that describes our chain of blobs. Actually, since we are planning to use `iwopy`'s fast vectorization capabilities, we directly implement a _population_ of chains. This means that each object of the class represents not one but `n_pop` chains, each containing `N` blobs. Note that the _individuals_ of a _population_ are completely independent of one another.\n", "\n", "The geometry of a chain is descibed by `N-1` angles in degrees, giving the direction from one blob centre to the next. When combined with the information of the radii of the blobs and the location of the first blob in the chain, these angles fully determine the location of each blob. They will later serve as the optimization variables of the problem:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class ChainPopulation:\n", " \"\"\" A polulation of chains with N blobs \"\"\"\n", "\n", " def __init__(self, n_pop, N, radii=1., xy0=0., alpha=0.):\n", " self.N = N\n", " self.n_pop = n_pop\n", "\n", " self.radii = np.zeros(N)\n", " self.radii[:] = radii\n", "\n", " self.xy = np.zeros((n_pop, N, 2))\n", " self.xy[:, 0] = xy0\n", "\n", " self.alpha = np.zeros((n_pop, N-1))\n", " self.dists = np.zeros((n_pop, N, N))\n", " self.set_alpha(alpha)\n", " \n", " def set_alpha(self, alpha):\n", " \"\"\" Set new alpha values and update data \"\"\"\n", " self.alpha[:] = alpha\n", " arad = self.alpha*np.pi/180.\n", " uv = np.stack([np.cos(arad), np.sin(arad)], axis=-1)\n", " for i in range(1, self.N):\n", " self.xy[:, i] = self.xy[:, i-1] + uv[:, i-1] * (\n", " self.radii[i-1] + self.radii[i])\n", " \n", " for i in range(self.N):\n", " d = self.xy - self.xy[:, i, None]\n", " self.dists[:, i] = np.linalg.norm(d, axis=-1)\n", "\n", " def get_fig(self, i=0, ax=None, title=None):\n", " \"\"\" Visualize the chain for a selected individual \"\"\"\n", " if ax is None:\n", " fig, ax = plt.subplots()\n", " else:\n", " fig = ax.get_figure()\n", " xy = self.xy[i]\n", " for pi, pxy in enumerate(xy):\n", " ax.add_patch(plt.Circle(pxy, self.radii[pi], color='orange'))\n", " rmax = np.max(self.radii)\n", " xy_imin = np.argmin(xy, axis=0)\n", " xy_imax = np.argmax(xy, axis=0)\n", " xy_min = xy[xy_imin, range(2)] - rmax\n", " xy_max = xy[xy_imax, range(2)] + rmax\n", " xy_del = xy_max - xy_min\n", " ax.set_xlim((xy_min[0] - 0.1*xy_del[0], xy_max[0] + 0.1*xy_del[0]))\n", " ax.set_ylim((xy_min[1] - 0.1*xy_del[1], xy_max[1] + 0.1*xy_del[1]))\n", " ax.set_aspect(\"equal\", adjustable=\"box\")\n", " ax.set_xlabel(\"x\")\n", " ax.set_ylabel(\"y\")\n", " ax.set_title(f\"N = {self.N}\" if title is None else title)\n", " return fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note the `get_fig` function which was added for results plotting in the end.\n", "\n", "Next, we define an optimization problem whose optimization variables are the angles of the chain of blobs. Our intention is the evaluation in vectorized form, as mitigated by the `apply_population` function. Here this simply updates the `alpha` values of the chain to those given by the optimizer:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class ChainProblem(Problem):\n", "\n", " def __init__(self, chain):\n", " super().__init__(name=\"chain_problem\")\n", " self.chain = chain\n", "\n", " def var_names_float(self):\n", " \"\"\" The variable names \"\"\"\n", " return [f\"alpha_{i:04}\" for i in range(self.chain.N - 1)]\n", "\n", " def initial_values_float(self):\n", " \"\"\" The initial variable values \"\"\"\n", " return self.chain.alpha[:-1]\n", "\n", " def min_values_float(self):\n", " \"\"\" The minimal variable values: 0 degrees \"\"\"\n", " return np.full(self.chain.N - 1, 0.)\n", "\n", " def max_values_float(self):\n", " \"\"\" The maximal variable values: 360 degrees \"\"\"\n", " return np.full(self.chain.N - 1, 360.)\n", "\n", " def apply_individual(self, vars_int, vars_float):\n", " \"\"\" Apply new variables from the optimizer \"\"\"\n", " self.chain.set_alpha(vars_float[None, :])\n", "\n", " def apply_population(self, vars_int, vars_float):\n", " \"\"\" Apply new variables from the optimizer \"\"\"\n", " self.chain.set_alpha(vars_float)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to avoid that the chain is crossing, or that any of the blobs that are not neighbours crash into each other. Let's introduce a constraint that prevent such solutions:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class NoCrossing(Constraint):\n", " \"\"\" The chain blobs are not allowed to cross \"\"\"\n", "\n", " def __init__(self, problem, tol=1e-3):\n", " super().__init__(\n", " problem, \"nocross\", vnames_float=problem.var_names_float(), tol=tol\n", " )\n", " self.chain = problem.chain\n", "\n", " def n_components(self):\n", " \"\"\" Each blob touches two neighbours only \"\"\"\n", " N = self.chain.N\n", " return int((N**2 - N - 2*(N - 1))/2)\n", "\n", " def calc_individual(self, vars_int, vars_float, problem_results, cmpnts=None):\n", " \"\"\" Positive values for too nearby 3rd and higher neighbours \"\"\"\n", " rmin = np.min(self.chain.radii)\n", " values = np.zeros(self.n_components())\n", " i0 = 0\n", " for i in range(self.chain.N - 2):\n", " i1 = i0 + self.chain.N - 2 - i\n", " meet = self.chain.dists[0, i, i+2:] - self.chain.radii[i] - self.chain.radii[i+2:]\n", " values[i0:i1] = 0.1 * rmin - meet\n", " i0 = i1\n", " \n", " return values\n", "\n", " def calc_population(self, vars_int, vars_float, problem_results, cmpnts=None):\n", " \"\"\" Positive values for too nearby 3rd and higher neighbours \"\"\"\n", " rmin = np.min(self.chain.radii)\n", " values = np.zeros((self.chain.n_pop, self.n_components()))\n", " i0 = 0\n", " for i in range(self.chain.N - 2):\n", " i1 = i0 + self.chain.N - 2 - i\n", " meet = self.chain.dists[:, i, i+2:] - self.chain.radii[i] - self.chain.radii[None, i+2:]\n", " values[:, i0:i1] = 0.1 * rmin - meet\n", " i0 = i1\n", " \n", " return values" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For `N` blobs in the chain, this defines `(N**2 - N - (N - 1))/2` costraint component functions. Imagine a matrix with `N x N` entries. Since only 3rd neighbours onwards correspond to a constraint, and there is no need to repeat constraints for backward located blobs on the chain, this constraint number represents the upper-diagonal content of that matrix when ignoring the diagonal and also the secondary diagonal.\n", "\n", "The objectives are straight-forward measures of the extension of given directions:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "class MaxStretch(Objective):\n", " \"\"\" Aim for maximal stretch along a given direction \"\"\"\n", "\n", " def __init__(self, problem, direction=np.array([0., 1.]), name=\"stretch\"):\n", " super().__init__(problem, name, vnames_float=problem.var_names_float())\n", " self.chain = problem.chain\n", " self.direction = direction\n", "\n", " def n_components(self):\n", " \"\"\" There is only one component for this objective \"\"\"\n", " return 1\n", "\n", " def maximize(self):\n", " \"\"\" The stretch length is to be maximized \"\"\"\n", " return [True]\n", "\n", " def calc_individual(self, vars_int, vars_float, problem_results, cmpnts=None):\n", " \"\"\" Calculate the stretch length \"\"\"\n", " u = np.einsum('cd,d->c', self.chain.xy[0], self.direction)\n", " return np.max(u + self.chain.radii) - np.min(u + self.chain.radii)\n", "\n", " def calc_population(self, vars_int, vars_float, problem_results, cmpnts=None):\n", " \"\"\" Calculate the stretch length \"\"\"\n", " u = np.einsum('pcd,d->pc', self.chain.xy, self.direction)[:, :, None]\n", " return np.max(u + self.chain.radii[None, :, None], axis=1) - np.min(u - self.chain.radii[None, :, None], axis=1)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This completes our elements for the description of our optimization problem. So let's create the problem and add objectives and constraints, here for maximal radius `r = 5` and `n = 10` blobs per chain:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "r = 5.0\n", "n = 10\n", "n_pop = 100\n", "\n", "radii = np.random.uniform(r/2., r, n)\n", "chain = ChainPopulation(n_pop, n, radii)\n", "\n", "problem = ChainProblem(chain)\n", "problem.add_constraint(NoCrossing(problem))\n", "problem.add_objective(MaxStretch(problem, direction=np.array([1.,0.]), name=\"stretch_x\"))\n", "problem.add_objective(MaxStretch(problem, direction=np.array([0.,1.]), name=\"stretch_y\"))\n", "problem.initialize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that we selected `n_pop = 100` individuals per population. \n", "\n", "We would like to solve this problem using `pymoo`'s `NSGA2` algorithm. Here we create the corresponding solver object and initialize it:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "solver = Optimizer_pymoo(\n", " problem,\n", " problem_pars=dict(\n", " vectorize=True,\n", " ),\n", " algo_pars=dict(\n", " type=\"NSGA2\",\n", " pop_size=n_pop,\n", " seed=42,\n", " ),\n", " setup_pars=dict(),\n", " term_pars=dict(\n", " type=\"default\",\n", " n_max_gen=200,\n", " ftol=0,\n", " xtol=0,\n", " ),\n", ")\n", "solver.initialize()\n", "solver.print_info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that we selected `n_max_gen = 200` generations in this setup. We are now ready to go! Let's run the solver:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "results = solver.solve(verbosity=0)\n", "solver.finalize(results)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All solution data is now contained in the `results` object. We can ask it to create a Pareto front plot:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ax = results.plot_pareto()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One of the nice things about `NSGA2` is that it produces a well covered Pareto front in the final population. We can have a look at individuals, picking them according to weights of the objectives:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for w in [[1, 0], [0.5, 0.5], [0, 1]]:\n", " i = results.find_pareto_objmix(w, max=True)\n", " fig = chain.get_fig(i, title=f\"Weights stretch x, y: {w}\")\n", "\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3.8.10 ('iwopyi')", "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.8.10 (default, Nov 14 2022, 12:59:47) \n[GCC 9.4.0]" }, "orig_nbformat": 4, "vscode": { "interpreter": { "hash": "0e385321928b3b4f0753cf8f84cadb34ba8fa899d98e78ebdf797f13a2c801ba" } } }, "nbformat": 4, "nbformat_minor": 2 }