{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b1e85461",
   "metadata": {},
   "source": [
    "# Multiojective Optimization of the Fonseca Fleming function\n",
    "\n",
    "For $n\\in \\mathbb{N}$, we consider the box-constrained bicriterial optimization of the Fonseca-Fleming function\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{gathered}%\\label{eq:P}\\tag{P}\n",
    "\\min\\limits_{x\\in \\mathbb{R}^n}\n",
    "f(x)=\\left( \\begin{array}{c}\n",
    "1-\\exp\\bigg(-\\sum\\limits_{i=1}^n\\Big(x_i-\\frac{1}{\\sqrt{n}}\\Big) ^2\\bigg)\\\\\n",
    "1-\\exp\\bigg(-\\sum\\limits_{i=1}^n\\Big(x_i+\\frac{1}{\\sqrt{n}}\\Big)^2\\bigg)\n",
    "\\end{array}\n",
    "\\right)\\\\\n",
    "\\begin{alignedat}{3}\n",
    "\\mbox{s.t.  }  -2 \\leq x_i \\leq 2 \\quad \\forall i\\in\\{1,\\ldots,n\\}.\n",
    "\\end{alignedat}\n",
    "\\end{gathered}\n",
    "\\end{equation*}\n",
    "\n",
    "This is a non-convex MOP since the objective function is not convex. For these types of problems the canonical choice of our optimization package is the Pascoletti-Serafini method (PSM), which will give the best results.  We will see that the WSM cannot compute the Pareto set in this case. However the ERPM can compute the Pareto front,  but not in an uniformly discretized manner.\n",
    "\n",
    "First, we import the modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52470953",
   "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 l_bfgs_b, augmented_lagrangian, projected_armijo, projected_gradient"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4d65d9b",
   "metadata": {},
   "source": [
    "Next, we set $n=2$ for simplicity, specify the cost functions and their derivatives and store them in a list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb39aaa6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set n\n",
    "n = 2\n",
    "\n",
    "# Define the cost functions\n",
    "def f1(x):\n",
    "    exponent = 0\n",
    "    for i in range(n):\n",
    "        exponent -= (x[i] - 1/np.sqrt(n))**2\n",
    "    return 1 - np.exp(exponent)\n",
    "\n",
    "\n",
    "def f2(x):\n",
    "    exponent = 0\n",
    "    for i in range(n):\n",
    "        exponent -= (x[i] + 1/np.sqrt(n))**2\n",
    "    return 1 - np.exp(exponent)\n",
    "\n",
    "# Define the derivatives\n",
    "def df1(x):\n",
    "    grad_val = np.zeros(x.shape[0])\n",
    "    exponent = 0\n",
    "    for i in range(x.shape[0]):\n",
    "        exponent -= (x[i] - 1/np.sqrt(n))**2\n",
    "    for i in range(x.shape[0]):\n",
    "        grad_val[i] = 2*(x[i] - 1/np.sqrt(n)) * np.exp(exponent)\n",
    "    return grad_val\n",
    "\n",
    "\n",
    "def df2(x):\n",
    "    grad_val = np.zeros(x.shape[0])\n",
    "    exponent = 0\n",
    "    for i in range(x.shape[0]):\n",
    "        exponent -= (x[i] + 1/np.sqrt(n))**2\n",
    "    for i in range(x.shape[0]):\n",
    "        grad_val[i] = 2*(x[i] + 1/np.sqrt(n)) * np.exp(exponent)\n",
    "    return grad_val\n",
    "\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",
    "\n",
    "# Store the cost functions, their gradients and Hessians in a list\n",
    "f = [f1, f2]  \n",
    "df = [df1, df2] \n",
    "ddf = [ddf1, ddf2] "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca876dd1",
   "metadata": {},
   "source": [
    "Further, we set the box constraints and define the initial guess for the solution of the first subproblem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95136883",
   "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 for the solution of the first subproblem\n",
    "x = np.array([0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62d72e0d",
   "metadata": {},
   "source": [
    "## Solution with the Pascoletti-Serafini method (PSM)\n",
    "\n",
    "We want to solve the problem with the PSM. As an solver for the scalar subproblems of the PSM we choose the Augmented Lagrange method.  Note that we need to extend the bounds for the Augmented Lagrangian, since it also approximates the multiplier of the box constraints and further the Augmented Lagrange method uses an inner solver, for which we choose the L-BFGS algorithm. The whole procedure is depicted below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6a6b5a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extend the box constraints for the augmented Lagrange\n",
    "low_aug = np.hstack((low, -np.Inf)) \n",
    "up_aug = np.hstack((up, np.Inf))\n",
    "\n",
    "# Inner optimization for the augmented Lagrangian algorithm\n",
    "def optimeth_inner_opt_aug_lagr(x, fhandle, dfhandle, eps):\n",
    "    # Define the line search algorithm\n",
    "    def default_linesearch(x, d, low, up):\n",
    "        t, flag = projected_armijo(fhandle, x, d, low, up, options=None)\n",
    "        return t, flag\n",
    "    # Set the options\n",
    "    options_l_bfgs = Options(tol_abs=eps, tol_rel=eps,\n",
    "                             linesearch=default_linesearch)\n",
    "    # Apply the optimization algorithm\n",
    "    result_opt = l_bfgs_b(fhandle, dfhandle, x, low_aug, up_aug,\n",
    "                          options=options_l_bfgs)\n",
    "    return result_opt\n",
    "\n",
    "\n",
    "# Outer optimization algorithm for the Pascoletti-Serafini problems (Augmented\n",
    "# Lagrangian algorithm)\n",
    "def optimeth_augm_lagr(fhandle, dfhandle, ghandle, dghandle, x,\n",
    "                       ddfhandle=None, ddghandle=None):\n",
    "    # Set the options\n",
    "    options_augm_lagr = Options(optimeth=optimeth_inner_opt_aug_lagr)\n",
    "    # Apply the optimization algorithm\n",
    "    result_opt = augmented_lagrangian(fhandle, dfhandle,\n",
    "                                      ghandle=ghandle,\n",
    "                                      dghandle=dghandle, x=x,\n",
    "                                      options=options_augm_lagr)\n",
    "    return result_opt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e293d99f",
   "metadata": {},
   "source": [
    "All in all, we specify the following options and set the problem_type to non-convex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a216e45",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the options for the PSM\n",
    "options_mop_PSM = Options(mop_method=\"PSM\", hp=0.2, hx=0.05,\n",
    "                          r_PSM=np.ones(len(f)),                         \n",
    "                          PSM_optimeth=optimeth_augm_lagr,\n",
    "                          problem_type=\"non-convex\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "734f0f69",
   "metadata": {},
   "source": [
    "Now we are ready to solve and store the computed data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06f691cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve the Multiobjective Optimization Problem with the PSM\n",
    "MOP, data = solve_MultiobjectiveOptimizationProblem(f, df, x, low, up,\n",
    "                                                    options_mop_PSM)\n",
    "# Get the Pareto optimal parameters\n",
    "x_opt = MOP[\"1,2\"].x\n",
    "\n",
    "# Get the Pareto front\n",
    "f_values = MOP[\"1,2\"].func_values\n",
    "\n",
    "# Get the reference points\n",
    "z_ref = MOP[\"1,2\"].z_ref_all"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bbed054",
   "metadata": {},
   "source": [
    "Next, we plot the computed Pareto set and the Pareto front with all reference points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13f971f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the Pareto set\n",
    "fig = plt.figure()\n",
    "plt.title(\"Pareto set\")\n",
    "plt.plot(x_opt[:, 0], x_opt[:, 1], 'o')\n",
    "plt.show()\n",
    "\n",
    "# Plot the Pareto front and the reference points\n",
    "fig = plt.figure()\n",
    "plt.title(\"Pareto Front\")\n",
    "plt.plot(f_values[:, 0], f_values[:, 1], 'o',label='Pareto front')\n",
    "plt.plot(z_ref[:, 0], z_ref[:, 1], 'o',label='Reference points')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('f1')\n",
    "plt.ylabel('f2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e69cf899",
   "metadata": {},
   "source": [
    "We observe that the PSM computes the Pareto set and the discretization of the Pareto front is uniformly.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba68b993",
   "metadata": {},
   "source": [
    "## Solution with Euclidian Reference Point Method (ERPM)\n",
    "\n",
    "Next, we try to solve the problem using the ERPM. Note that from a theoretical point of view, the ERPM fits only for convex problems. However, in this special case we can compute the Pareto front, but not homogeneously discretized. \n",
    "\n",
    "We define the Projected Gradient method as an inner solver for the ERPM and choose the following options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae4c63c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the Projected Gradient method as an inner solver for the ERPM\n",
    "def optimeth_proj_grad(fhandle, dfhandle, x, low, up, ddfhandle=None):\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",
    "    options_proj_grad = Options(nmax=1000,\n",
    "                                disp=False, stop=False, statusdisp=True,\n",
    "                                linesearch=linesearch_inner_proj_armijo)\n",
    "    result_opt = projected_gradient(fhandle, dfhandle, x, low, up,\n",
    "                                    options_proj_grad)\n",
    "    return result_opt\n",
    "\n",
    "# Options for the ERPM\n",
    "options_mop_ERPM = Options(mop_method=\"ERPM\", hp=0.01, hx=0.05,\n",
    "                           ERPM_optimeth=optimeth_proj_grad,\n",
    "                           problem_type=\"non-convex\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d661461c",
   "metadata": {},
   "source": [
    "We are ready to solve the problem and store the computed data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e75c6d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve the problem with the ERPM\n",
    "MOP, data = solve_MultiobjectiveOptimizationProblem(f, df, x, low, up,\n",
    "                                                    options_mop_ERPM)\n",
    "\n",
    "# Get the Pareto optimal parameters\n",
    "x_opt = MOP[\"1,2\"].x\n",
    "\n",
    "# Get the Pareto front\n",
    "f_values = MOP[\"1,2\"].func_values\n",
    "\n",
    "# Get the reference points\n",
    "z_ref = MOP[\"1,2\"].z_ref_all"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c675f9f7",
   "metadata": {},
   "source": [
    "And we plot the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da632d93",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the Pareto set\n",
    "fig = plt.figure()\n",
    "plt.title(\"Pareto set\")\n",
    "plt.plot(x_opt[:, 0], x_opt[:, 1], 'o')\n",
    "plt.show()\n",
    "\n",
    "# Plot the Pareto front and the reference points\n",
    "fig = plt.figure()\n",
    "plt.title(\"Pareto Front\")\n",
    "plt.plot(f_values[:, 0], f_values[:, 1], 'o',label='Pareto front')\n",
    "plt.plot(z_ref[:, 0], z_ref[:, 1], 'o',label='Reference points')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('f1')\n",
    "plt.ylabel('f2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3da8340e",
   "metadata": {},
   "source": [
    "One can see the structure of the Pareto front, but the ERPM fails in discretizing the Pareto front homogeneously.\n",
    "\n",
    "## Solution with Weighted-Sum Method (WSM)\n",
    "\n",
    "This example should indicate what can go wrong when using the WSM to solve a non-convex problem. We specify the following options and solve the problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b17c433",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify 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=\"non-convex\")\n",
    "\n",
    "# Solve the MOP with the WSM\n",
    "MOP, data = solve_MultiobjectiveOptimizationProblem(f, df, x, low, up,\n",
    "                                                    options_mop_WSM)\n",
    "\n",
    "# Get the Pareto optimal parameters\n",
    "x_opt = MOP[\"1,2\"].x\n",
    "\n",
    "# Get the Pareto front\n",
    "f_values = MOP[\"1,2\"].func_values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a94683e4",
   "metadata": {},
   "source": [
    "Finally we plot the results. We observe that the WSM fails to compute the Pareto set/front. Therefore one should pay attention when using the WSM for solving non-convex problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34f1ee63",
   "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.show()\n",
    "\n",
    "# Plot the Pareto front and the reference points\n",
    "fig = plt.figure()\n",
    "plt.title(\"Pareto Front\")\n",
    "plt.plot(f_values[:, 0], f_values[:, 1], 'o',label='Pareto front')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('f1')\n",
    "plt.ylabel('f2')\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
}