practice-pykep-lt_direct_pl2pl.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7e3a108e",
   "metadata": {},
   "source": [
    "# Practice using pykep's `direct_pl2pl` UDP and pygmo\n",
    "\n",
    "pykep doc: https://esa.github.io/pykep/documentation/trajopt.html#pykep.trajopt.direct_pl2pl\n",
    "\n",
    "Yuri Shimane, 2022.04.21"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "206b1f98",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pygmo as pg\n",
    "import pykep as pk\n",
    "\n",
    "%matplotlib notebook\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4a0774e",
   "metadata": {},
   "source": [
    "### Launch epoch\n",
    "\n",
    "First, a word on launch epoch. Pykep's UDPs usually want this in mjd2000. So we first create `pk.epoch` objects, which we then convert to floats. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "3b1c86be",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10958.515972222222"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epoch1 = pk.epoch_from_string(\"2030-01-01 12:23:00.000\")\n",
    "epoch2 = pk.epoch_from_string(\"2032-01-01 12:23:00.000\")\n",
    "epoch1.mjd2000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "c70cd887",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we construct our problem\n",
    "prob = pk.trajopt.direct_pl2pl(\n",
    "    p0 = \"earth\",\n",
    "    pf = \"mars\",\n",
    "    mass = 2000,   # kg\n",
    "    thrust = 0.4,  # thrust in N\n",
    "    isp = 3000,    # isp in seconds\n",
    "    nseg = 30,     # number of segments in Sims-Flanagan transcription\n",
    "    t0 = [epoch1.mjd2000, epoch2.mjd2000],  # Launch epochs bounds [mjd2000]\n",
    "    tof = [150, 500],  # \n",
    "    vinf_dep = 2.0,   # Launch DV, km/s\n",
    "    vinf_arr = 0.0,   # Arrival DV, km/s\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "576e7bb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "class add_gradient:\n",
    "    \"\"\"Construct pygmo UDP with gradient\n",
    "    https://esa.github.io/pygmo2/tutorials/coding_udp_constrained.html?highlight=add_gradient\n",
    "    \"\"\"\n",
    "    def __init__(self, prob):\n",
    "        self.prob = pg.problem(prob)\n",
    "\n",
    "    def fitness(self, x):\n",
    "        return self.prob.fitness(x)\n",
    "\n",
    "    def get_bounds(self):\n",
    "        return self.prob.get_bounds()\n",
    "\n",
    "    def get_nec(self):\n",
    "        return self.prob.get_nec()\n",
    "\n",
    "    def get_nic(self):\n",
    "        return self.prob.get_nic()\n",
    "\n",
    "    def get_nobj(self):\n",
    "        return self.prob.get_nobj()\n",
    "        \n",
    "    def gradient(self, x):\n",
    "        return pg.estimate_gradient(lambda x: self.fitness(x), x) # we here use the low precision gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "aa9e053c",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob_with_gradient = add_gradient(prob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "id": "3dd0496e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create ipopt algorithm class\n",
    "uda = pg.ipopt() # pylint: disable=no-member\n",
    "uda.set_integer_option(\"print_level\", 5)\n",
    "uda.set_integer_option(\"acceptable_iter\", 4)\n",
    "uda.set_integer_option(\"max_iter\", 5000)      # probably need thousands of iterations\n",
    "\n",
    "uda.set_numeric_option(\"tol\", 1e-5)\n",
    "uda.set_numeric_option(\"dual_inf_tol\", 1e-6)\n",
    "uda.set_numeric_option(\"constr_viol_tol\", 1e-5)\n",
    "uda.set_numeric_option(\"compl_inf_tol\", 1e-6)\n",
    "\n",
    "uda.set_numeric_option(\"acceptable_tol\", 1e-3)\n",
    "uda.set_numeric_option(\"acceptable_dual_inf_tol\", 1e-2)\n",
    "uda.set_numeric_option(\"acceptable_constr_viol_tol\", 1e-5)\n",
    "uda.set_numeric_option(\"acceptable_compl_inf_tol\", 1e-5)\n",
    "\n",
    "algo = pg.algorithm(uda)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "a09d5f86",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Algorithm name: Ipopt: Interior Point Optimization [deterministic]\n",
       "\tC++ class name: class pagmo::ipopt\n",
       "\n",
       "\tThread safety: none\n",
       "\n",
       "Extra info:\n",
       "\tLast optimisation return code: Solve_Succeeded (value = 0)\n",
       "\tVerbosity: 0\n",
       "\tIndividual selection policy: best\n",
       "\tIndividual replacement policy: best\n",
       "\tInteger options: {acceptable_iter : 4,  max_iter : 5000,  print_level : 5}\n",
       "\tNumeric options: {acceptable_compl_inf_tol : 1e-05,  acceptable_constr_viol_tol : 1e-05,  acceptable_dual_inf_tol : 0.01,  acceptable_tol : 0.001,  compl_inf_tol : 1e-06,  ...}"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "algo"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3a5c7a8",
   "metadata": {},
   "source": [
    "# Create a population\n",
    "\n",
    "We now create a population. Here, we are not trying to do a thorough global search, so we initialize the population with `size=1`. This is equivalent to using the algorithm with a randomized initial guess, once. \n",
    "If we want to use an initial guess that we have a priori (say `x0`), then we should initialize the population with `size=0` then push back:\n",
    "\n",
    "```python\n",
    "# ... assume we have some initial guess vector x0 that we'd like to use\n",
    "pop = pg.population(prob_with_gradient, size=0)   # initialize with size=0\n",
    "pop.push_back(xopt)                               # now size is +1\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "927abaf6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Problem name: <class '__main__.add_gradient'>\n",
       "\tC++ class name: class pybind11::object\n",
       "\n",
       "\tGlobal dimension:\t\t\t99\n",
       "\tInteger dimension:\t\t\t0\n",
       "\tFitness dimension:\t\t\t40\n",
       "\tNumber of objectives:\t\t\t1\n",
       "\tEquality constraints dimension:\t\t7\n",
       "\tInequality constraints dimension:\t32\n",
       "\tTolerances on constraints: [0, 0, 0, 0, 0, ... ]\n",
       "\tLower bounds: [10958.5, 150, 200, -2000, -2000, ... ]\n",
       "\tUpper bounds: [11688.5, 500, 2000, 2000, 2000, ... ]\n",
       "\tHas batch fitness evaluation: false\n",
       "\n",
       "\tHas gradient: true\n",
       "\tUser implemented gradient sparsity: false\n",
       "\tExpected gradients: 3960\n",
       "\tHas hessians: false\n",
       "\tUser implemented hessians sparsity: false\n",
       "\n",
       "\tFitness evaluations: 1\n",
       "\tGradient evaluations: 0\n",
       "\n",
       "\tThread safety: none\n",
       "\n",
       "Population size: 1\n",
       "\n",
       "List of individuals: \n",
       "#0:\n",
       "\tID:\t\t\t12692838010386895863\n",
       "\tDecision vector:\t[11246.6, 169.257, 1190.45, 1782.38, 562.858, ... ]\n",
       "\tFitness vector:\t\t[-1190.45, 1.37493, 0.758595, -0.104617, -0.825239, ... ]\n",
       "\n",
       "Champion decision vector: [11246.6, 169.257, 1190.45, 1782.38, 562.858, ... ]\n",
       "Champion fitness: [-1190.45, 1.37493, 0.758595, -0.104617, -0.825239, ... ]"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pop = pg.population(prob_with_gradient, size=1)\n",
    "pop"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca9cb86f",
   "metadata": {},
   "source": [
    "**CAVEAT**: when running with Jupyter notebooks, the algorithm iteration messages are shown on the shell session running Jupyter notebook, not on the notebook `Out` cell itself..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a68ed7ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# solving with ipopt\n",
    "pop = algo.evolve(pop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be06d507",
   "metadata": {},
   "outputs": [],
   "source": [
    "algo"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40565da8",
   "metadata": {},
   "source": [
    "### Visualize result\n",
    "\n",
    "See also: https://esa.github.io/pykep/documentation/trajopt.html#pykep.trajopt.direct_pl2pl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70a2e1fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "xopt, fopt = pop.champion_x, pop.champion_f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee7ba78b",
   "metadata": {},
   "outputs": [],
   "source": [
    "pop.champion_f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c194fe06",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob.pretty(xopt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dfe75b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob.plot_traj(xopt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdc5ab98",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob.plot_control(xopt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16644998",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f76ecfd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "853e1ba5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d6f0167",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}