{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8e6a8a74",
   "metadata": {},
   "source": [
    "# Solving a Convex Multiobjective Optimization Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2355add5",
   "metadata": {},
   "source": [
    "Given the vector-valued objective function $f:\\mathbb{R}^2\\to \\mathbb{R}^3$, we consider the following multiobjective optimization problem (MOP)\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{gathered}%\\label{eq:P}\\tag{P}\n",
    "\\min\\limits_{x\\in \\mathbb{R}^2}\n",
    "f(x) = \\left( \n",
    "\\begin{array}\n",
    "\\vert \\| x-v_1\\|^2_2\\\\\n",
    "\\|x-v_2\\|^2_2\\\\\n",
    "\\|x-v_3\\|^2_2\n",
    "\\end{array}\n",
    " \\right)\\\\\n",
    "\\begin{alignedat}{3}\n",
    "\\mbox{s.t.  }  -2 \\leq x_i \\leq 2 \\quad \\forall i\\in\\{1,2\\}.\n",
    "\\end{alignedat}\n",
    "\\end{gathered}\n",
    "\\end{equation*}\n",
    "\n",
    "Further, we set the vertices of the three parabolas as $v_1=[1,1]^T,\\ v_2=[-0.5,-1]^T$ and $v_3=[1.5,-1.5]^T$ and $\\|\\cdot\\|_2$ denotes the euclidian distance in $\\mathbb{R}^2$. \n",
    "\n",
    "Since all components of the objective function $f$ and the costraints are convex, the MOP is convex, which implies that the Pareto set consists of one connected component. Further, it will turn out that the Pareto set consists of the convex hull of the points $v_1,\\ v_2,\\ v_3$. For convex MOPs the weighted-sum method (WSM) and euclidian reference point are well suited and their application on the MOP will be demonstrated below.\n",
    "\n",
    "At first, we load the modules, define the objective function and its derivatives and store them in a list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "275e0b29",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from oppy.multOpt.mop import solve_MultiobjectiveOptimizationProblem\n",
    "import matplotlib.pyplot as plt\n",
    "from oppy.options import Options\n",
    "from oppy.conOpt import projected_armijo, projected_gradient\n",
    "\n",
    "# Specify the vertices of the parabolas\n",
    "v1 = np.array([1, 1])\n",
    "v2 = np.array([-0.5, -1])\n",
    "v3 = np.array([1.5, -1.5])\n",
    "\n",
    "# Define the cost functions\n",
    "f1 = (lambda x: (x[0] - v1[0])**2 + (x[1] - v1[1])**2)\n",
    "f2 = (lambda x: (x[0] - v2[0])**2 + (x[1] - v2[1])**2)\n",
    "f3 = (lambda x: (x[0] - v3[0])**2 + (x[1] - v3[1])**2)\n",
    "\n",
    "# Define the gradients\n",
    "df1 = (lambda x: np.array([2*(x[0] - v1[0]), 2*(x[1] - v1[1])]))\n",
    "df2 = (lambda x: np.array([2*(x[0] - v2[0]), 2*(x[1] - v2[1])]))\n",
    "df3 = (lambda x: np.array([2*(x[0] - v3[0]), 2*(x[1] - v3[1])]))\n",
    "\n",
    "# Define the Hessians\n",
    "ddf1 = (lambda x: np.array([[2, 0], [0, 2]]))\n",
    "ddf2 = (lambda x: np.array([[2, 0], [0, 2]]))\n",
    "ddf3 = (lambda x: np.array([[2, 0], [0, 2]]))\n",
    "\n",
    "# Store the cost functions, their gradients and Hessians in a list\n",
    "f = [f1, f2, f3] \n",
    "df = [df1, df2, df3] \n",
    "ddf = [ddf1, ddf2, ddf3]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf4d84cc",
   "metadata": {},
   "source": [
    "Next, we define the lower and upper bounds for the box constraints and the initial value for the minimization of the first scalar subproblem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d0c90fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the box constraints\n",
    "low = np.array([-2, -2])\n",
    "up = np.array([2, 2])\n",
    "\n",
    "# Set the initial value\n",
    "x = np.array([0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c35399e9",
   "metadata": {},
   "source": [
    "## Solution with the Weighted-Sum Method (WSM)\n",
    "\n",
    "First, we want to solve the convex MOP with the WSM. In order to do that we define the following options for the optimization. We choose the projected Gradient method with an Armijo line search as an inner solver for the scalar subproblems. Note that the only possibility to define a line search for the inner solver is the way, which is depicted below. Further, we set number of weights per dimension to 30 and mark the problem type as convex to improve the efficiency of the solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcc0f31f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the solver function for the scalar subproblems\n",
    "def optimeth_proj_grad(fhandle, dfhandle, x, low, up, ddfhandle=None):\n",
    "    \n",
    "    # Define line search\n",
    "    def linesearch_inner_proj_armijo(x, d, low, up):  \n",
    "        t, flag = projected_armijo(fhandle, x, d, low, up, options=None)\n",
    "        return t, flag\n",
    "    \n",
    "    # Set the options\n",
    "    options_proj_grad = Options(nmax=1000,\n",
    "                                disp=False, stop=False, statusdisp=True,\n",
    "                                linesearch=linesearch_inner_proj_armijo)\n",
    "    \n",
    "    # Apply the optimization algorithm\n",
    "    result_opt = projected_gradient(fhandle, dfhandle, x, low, up,\n",
    "                                    options_proj_grad)\n",
    "    return result_opt\n",
    "\n",
    "# Define the options for the WSM\n",
    "options_mop_WSM = Options(mop_method=\"WSM\", num_weights_per_dimension=30,\n",
    "                          WSM_optimeth=optimeth_proj_grad,\n",
    "                          problem_type=\"convex\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da7eaf12",
   "metadata": {},
   "source": [
    "Now we are ready to solve the MOP with the WSM and save the Pareto set and the functional values of the objectives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07ab2a71",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve the Multiobjective Optimization Problem with the WSM\n",
    "MOP, data = solve_MultiobjectiveOptimizationProblem(f, df, x, low, up,\n",
    "                                                    options_mop_WSM)\n",
    "\n",
    "# Get the Pareto set\n",
    "x_opt = MOP[\"1,2,3\"].x \n",
    "\n",
    "# Get the Pareto front\n",
    "f_values = MOP[\"1,2,3\"].func_values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "426fed4c",
   "metadata": {},
   "source": [
    "Now we can plot the resulting Pareto set and the Pareto front. Note that the Pareto set is given as the convex hull of $v_1,v_2,v_3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34396db0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the Pareto optimal parameters and the points v_i\n",
    "fig = plt.figure()\n",
    "plt.title(\"Pareto set\")\n",
    "plt.plot(x_opt[:, 0], x_opt[:, 1], 'o')\n",
    "plt.plot(v1[0], v1[1], 'rx')\n",
    "plt.plot(v2[0], v2[1], 'rx')\n",
    "plt.plot(v3[0], v3[1], 'rx')\n",
    "plt.show()\n",
    "\n",
    "# Plot the Pareto front\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "plt.title(\"Pareto Front\")\n",
    "ax.scatter(f_values[:, 0], f_values[:, 1], f_values[:, 2], 'o',label='Pareto front')\n",
    "ax.legend(loc='best')\n",
    "ax.set_xlabel('X Label')\n",
    "ax.set_ylabel('Y Label')\n",
    "ax.set_zlabel('Z Label')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfcbdc3b",
   "metadata": {},
   "source": [
    "## Solution with Euclidian Reference Point Method (ERPM)\n",
    "\n",
    "Next, we solve the problem using the ERPM. We specify the following options, where we set $h_p=0.2$ and $h_x=0.8$. Again we choose the Projected Gradient method as an inner solver for the scalar subproblems and set the problem type to convex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "669cab25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the options for the ERPM\n",
    "options_mop_ERPM = Options(mop_method=\"ERPM\", hp=0.2, hx=0.8,                          \n",
    "                           ERPM_optimeth=optimeth_proj_grad,\n",
    "                           problem_type=\"convex\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "787d0894",
   "metadata": {},
   "source": [
    "Now we are ready to solve again and save the computed Pareto set, the Pareto front and all computed reference points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c3308e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve the Multiobjective Optimization Problem with the ERPM\n",
    "MOP, data = solve_MultiobjectiveOptimizationProblem(f, df, x, low, up,\n",
    "                                                    options_mop_ERPM)\n",
    "\n",
    "# Get the Pareto set\n",
    "x_opt = MOP[\"1,2,3\"].x\n",
    "\n",
    "# Get the Pareto front\n",
    "f_values = MOP[\"1,2,3\"].func_values\n",
    "\n",
    "# Get the reference points\n",
    "z_ref = MOP[\"1,2,3\"].z_ref_all"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b5c1914",
   "metadata": {},
   "source": [
    "Finally we plot the obtained results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08a2eba6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the Pareto optimal parameters\n",
    "fig = plt.figure()\n",
    "plt.title(\"Pareto set\")\n",
    "plt.plot(x_opt[:, 0], x_opt[:, 1], 'o')\n",
    "plt.plot(v1[0], v1[1], 'rx')\n",
    "plt.plot(v2[0], v2[1], 'rx')\n",
    "plt.plot(v3[0], v3[1], 'rx')\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "plt.title(\"Pareto Front\")\n",
    "ax.scatter(f_values[:, 0], f_values[:, 1], f_values[:, 2], 'o',label='Pareto front')\n",
    "ax.scatter(z_ref[:, 0], z_ref[:, 1], z_ref[:, 2], 'o', label='Reference points')\n",
    "ax.legend(loc='best')\n",
    "ax.set_xlabel('f1')\n",
    "ax.set_ylabel('f2')\n",
    "ax.set_zlabel('f3')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}