ActiveDebrisRemoval.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "ActiveDebrisRemoval.ipynb",
      "provenance": [],
      "collapsed_sections": [
        "dJcjHJNmtmfi",
        "4aZOJqGwtrTP",
        "QW4852RdtwFn"
      ]
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "Mjyn2PjQZ7Iq"
      },
      "source": [
        "!pip install pykep\n",
        "!pip install -U TLE-tools\n",
        "!pip install astropy"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kRlT-QOgMs9K"
      },
      "source": [
        "import random\n",
        "import bisect\n",
        "import numpy\n",
        "import scipy\n",
        "import copy\n",
        "from datetime import datetime\n",
        "from datetime import timedelta\n",
        "# -- for debris\n",
        "import math\n",
        "import csv\n",
        "from google.colab import files\n",
        "from google.colab import drive\n",
        "from tletools import TLE\n",
        "from astropy import constants\n",
        "import pykep\n",
        "# -- for ploting\n",
        "import datetime\n",
        "import json\n",
        "import time\n",
        "import networkx as nx\n",
        "from matplotlib import pyplot as plt"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dJcjHJNmtmfi"
      },
      "source": [
        "# Genetic algorithm\n",
        "The implementation uses the inver-over genetic operator to optimize the static sequence of debris based on the transference cost of the arcs.\n",
        "\n",
        "Also, the implementation uses **index_frozen** to model the already deorbited debris."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kghx_D7uJ_9X"
      },
      "source": [
        "class GA:\n",
        "  def __init__(self, population, fn_fitness, subpath_fn_fitness=None):\n",
        "    self.population = population\n",
        "    self.index_frozen = -1\n",
        "    self.fitnesses = [] # fitness for each individual in population\n",
        "    self.fn_fitness = fn_fitness # fitness function for the whole path\n",
        "    self.subpath_fn_fitness = subpath_fn_fitness # fitness function for a subpath\n",
        "\n",
        "  # freezes a debris in all individuals\n",
        "  def freeze_first(self, frozen):\n",
        "    self.index_frozen += 1\n",
        "    for i in range(len(self.population)):\n",
        "      del self.population[i][self.population[i].index(frozen)]\n",
        "      self.population[i].insert(self.index_frozen, frozen)\n",
        "\n",
        "  # decay a debris in all individuals\n",
        "  def decay(self, decayed_debris):\n",
        "    for i in range(len(self.population)):\n",
        "      for x in decayed_debris:\n",
        "        if x in self.population[i]:\n",
        "          del self.population[i][self.population[i].index(x)]\n",
        "\n",
        "  # force a first debris for all individuals\n",
        "  def startBy(self, debris):\n",
        "    for i in range(len(self.population)):\n",
        "      pos = self.population[i].index(debris)\n",
        "      self.population[i] = self.population[i][pos:] + self.population[i][:pos]\n",
        "\n",
        "  # returns the best individual\n",
        "  def getBest(self):\n",
        "    self.fit_population()\n",
        "    best = min(self.fitnesses)\n",
        "    return self.population[self.fitnesses.index(best)]\n",
        "\n",
        "  # run the inverover to optimize the static case\n",
        "  \"\"\"\n",
        "  tinv : int : number of iterations\n",
        "  feach : int : milestone to run kopt on the population\n",
        "  runkopt : int : iterations of kopt\n",
        "  forn : int : how many of the best individuals goes to kopt\n",
        "  \"\"\"\n",
        "  def run_inverover(self, tinv=1000, feach=1000, runkopt=100, forn=None):\n",
        "    self.fit_population()\n",
        "    self.inver_over(tinv, feach, runkopt, forn)\n",
        "    self.fit_population()\n",
        "    best = min(self.fitnesses)\n",
        "    return self.population[self.fitnesses.index(best)]\n",
        "\n",
        "  # select a random element of the population\n",
        "  def selectUniform(self):\n",
        "    return self.population[random.randrange(0, len(self.population))]\n",
        "\n",
        "  # calculate the fitness for all individuals\n",
        "  def fit_population(self):\n",
        "    if self.index_frozen >= 0:\n",
        "      self.fitnesses = list(map(lambda x: self.subpath_fn_fitness(x[self.index_frozen:]), self.population))\n",
        "    else:\n",
        "      self.fitnesses = list(map(lambda x: self.fn_fitness(x), self.population))\n",
        "\n",
        "  # run the stochastic kopt for the population\n",
        "  \"\"\"\n",
        "  permuts : int : number of iterations\n",
        "  elite : int : how many of the best shoud be processed\n",
        "  \"\"\"\n",
        "  def koptStochastic(self, permuts=100, elite=None):\n",
        "    indexes = range(len(self.population))\n",
        "    if elite is not None:\n",
        "      indexes = numpy.array(self.fitnesses).argsort()[:elite]\n",
        "\n",
        "    for x in indexes:\n",
        "      indv = self.population[x]\n",
        "      useds = {}\n",
        "      changed = False\n",
        "\n",
        "      for _ in range(0, permuts):\n",
        "        valid = False\n",
        "        while not valid:\n",
        "          i = random.randrange(self.index_frozen+1, len(indv))\n",
        "          j = i\n",
        "\n",
        "          while j == i: j = random.randrange(self.index_frozen+1, len(indv))\n",
        "\n",
        "          if (i, j) not in useds:\n",
        "            valid = True\n",
        "\n",
        "        useds[(i, j)] = True\n",
        "\n",
        "        if j < i:\n",
        "          temp = i\n",
        "          i = j\n",
        "          j = temp\n",
        "\n",
        "        if self.subpath_fn_fitness(list(reversed(indv[i:j+1]))) < self.subpath_fn_fitness(indv[i:j+1]):\n",
        "          changed = True\n",
        "          indv = indv[0:i] + list(reversed(indv[i:j+1])) + indv[j+1:]\n",
        "\n",
        "      if changed:\n",
        "        self.population[x] = indv\n",
        "        self.fitnesses[x] = self.subpath_fn_fitness(indv[self.index_frozen+1:])\n",
        "\n",
        "  # run the ranged kopt for one individual\n",
        "  \"\"\"\n",
        "  indv : array : the individual\n",
        "  maxrange : int : the range of analysis around the individual\n",
        "  \"\"\"\n",
        "  def ranged2opt(self, indv, maxrange=10):\n",
        "    ranger = indv[len(indv)-maxrange:] + indv[self.index_frozen+1: self.index_frozen+maxrange+2]\n",
        "    if len(set(ranger)) != len(ranger):\n",
        "      return indv\n",
        "\n",
        "    fit = self.subpath_fn_fitness(ranger)\n",
        "    changed = True\n",
        "\n",
        "    while changed:\n",
        "      changed = False\n",
        "      for i in range(len(ranger)):\n",
        "        for j in range(len(ranger)):\n",
        "          new_ranger = ranger[0:i] + list(reversed(ranger[i:j+1])) + ranger[j+1:]\n",
        "          new_fit = self.subpath_fn_fitness(new_ranger)\n",
        "          if new_fit < fit:\n",
        "            fit = new_fit\n",
        "            ranger = new_ranger\n",
        "            changed = True\n",
        "            break\n",
        "        if changed:\n",
        "          break\n",
        "\n",
        "    indv[len(indv)-maxrange:] = ranger[:maxrange]\n",
        "    indv[self.index_frozen+1: self.index_frozen+maxrange+2] = ranger[maxrange:]\n",
        "\n",
        "    return indv\n",
        "\n",
        "  # run the inverover for the population\n",
        "  \"\"\"\n",
        "  tinv : int : number of iterations\n",
        "  feach : int : milestone to run kopt on the population\n",
        "  runkopt : int : iterations of kopt\n",
        "  forn : int : how many of the best individuals goes to kopt\n",
        "  \"\"\"\n",
        "  def inver_over(self, tinv, feach, runkopt, forn):\n",
        "    for w in range(tinv):\n",
        "\n",
        "      if w % feach == 0:\n",
        "        self.koptStochastic(runkopt, forn)\n",
        "\n",
        "      for i in range(len(self.population)):\n",
        "        tmp = self.population[i]\n",
        "        c1 = tmp[random.randrange(0, len(tmp))]\n",
        "\n",
        "        changed = False\n",
        "\n",
        "        while True:\n",
        "          sl = self.population[i]\n",
        "          c2 = c1\n",
        "\n",
        "          while sl == self.population[i]: sl = self.selectUniform()\n",
        "          c2 = sl[(sl.index(c1) + 1) % len(sl)]\n",
        "\n",
        "          pos_c1 = tmp.index(c1)\n",
        "          pos_c2 = tmp.index(c2)\n",
        "\n",
        "          # if the genes are adjacent\n",
        "          if c2 in [ tmp[pos_c1-1], tmp[(pos_c1 + 1) % len(tmp)] ]:\n",
        "            break\n",
        "          # elif and else reverse a subset of chromosome\n",
        "          elif pos_c2 > pos_c1:\n",
        "            changed = True\n",
        "            c1 = tmp[(pos_c2 + 1) % len(tmp)]\n",
        "            tmp = tmp[:pos_c1+1] + list(reversed(tmp[pos_c1+1:pos_c2+1])) + tmp[pos_c2+1:]\n",
        "          else:\n",
        "            changed = True\n",
        "            c1 = tmp[pos_c2-1]\n",
        "            inverted = list(reversed(tmp[pos_c1:] + tmp[:pos_c2]))\n",
        "            div_pos = len(tmp)-pos_c1\n",
        "            tmp = inverted[div_pos:] + tmp[pos_c2:pos_c1] + inverted[:div_pos]\n",
        "\n",
        "        if changed:\n",
        "          fit_tmp = self.fn_fitness(tmp)\n",
        "          if fit_tmp < self.fitnesses[i]:\n",
        "            self.population[i] = tmp\n",
        "            self.fitnesses[i] = fit_tmp"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4aZOJqGwtrTP"
      },
      "source": [
        "# Problem instance\n",
        "The active debris removal problem is going to be modeled as a complex variant of Traveling Salesman Problem (TSP), the time-dependent TSP (TDTSP).\n",
        "\n",
        "The debris are the nodes and the dynamic transference trajectories are the edges.\n",
        "\n",
        "Also, the Max Open Walk is used to find for the optimized subpath."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-KWJjRBXKDXT"
      },
      "source": [
        "class StaticDebrisTSP:\n",
        "  mu = 398600800000000 # gravitational parameter of earth\n",
        "  re = 6378100 # radius of earth\n",
        "\n",
        "  def __init__(self, debris=[], weight_matrix=[], reward_matrix=[], path_size=0, population_size=100, epoch=None, hohmanncost=False):\n",
        "    self.index_frozen = -1\n",
        "    self.debris = debris # the debris cloud\n",
        "    self.reward_matrix = reward_matrix # the removal reward per debris\n",
        "    self.kepler_elements = [] # kepler elements of the debris\n",
        "    self.decayed_debris = [] # decayed debris\n",
        "    self.hohmanncost=hohmanncost # if the cost is calculated with hohmann\n",
        "\n",
        "    if epoch is not None:\n",
        "      self.epoch = epoch\n",
        "    else:\n",
        "      epoch = pykep.epoch_from_string(\"2021-06-11 00:06:09\")\n",
        "\n",
        "    is_matrix = len(weight_matrix) != 0\n",
        "\n",
        "    # size of a indivual\n",
        "    self.size = path_size if path_size != 0 else (len(weight_matrix) if is_matrix else len(debris))\n",
        "\n",
        "    # random population that will be used just as an input for the GA\n",
        "    self.population = []\n",
        "    for i in range(0, population_size):\n",
        "      self.population.append(random.sample(range(0, self.size), self.size))\n",
        "\n",
        "    # eighter receive the weight matrix or calculate it\n",
        "    if is_matrix:\n",
        "      self.fitness_matrix = weight_matrix\n",
        "\n",
        "    else:\n",
        "      # remove decayed debris\n",
        "      i = 0\n",
        "      count = 0\n",
        "      qtd_decayed = 0\n",
        "      while count < self.size:\n",
        "        if i >= len(debris):\n",
        "          break\n",
        "\n",
        "        try:\n",
        "          self.kepler_elements.append(debris[i].osculating_elements(self.epoch))\n",
        "          count += 1\n",
        "        except:\n",
        "          self.decayed_debris.append(i)\n",
        "          qtd_decayed += 1\n",
        "        i += 1\n",
        "\n",
        "      print('Decayed debris ', qtd_decayed, 'Total ', len(self.kepler_elements))\n",
        "      if len(self.kepler_elements) < self.size:\n",
        "        raise BaseException('Insuficient size')\n",
        "\n",
        "      # fitness matrix\n",
        "      self.fitness_matrix = numpy.zeros((self.size, self.size))\n",
        "      for i in range(0, self.size):\n",
        "        for j in range(0, self.size):\n",
        "          if self.hohmanncost:\n",
        "            self.fitness_matrix[i][j] = StaticDebrisTSP.MYhohmann_impulse_aprox(self.kepler_elements[i], self.kepler_elements[j], self.epoch)\n",
        "          else:\n",
        "            try:\n",
        "              self.fitness_matrix[i][j] = pykep.phasing.three_impulses_approx(debris[i], debris[j], self.epoch, self.epoch)\n",
        "            except:\n",
        "              d1 = self.kepler_elements[i]\n",
        "              d2 = self.kepler_elements[j]\n",
        "              self.fitness_matrix[i][j] = StaticDebrisTSP.MYthree_impulse_aprox(d1[0],d1[1],d1[2],d1[3],d2[0],d2[1],d2[2],d2[3],StaticDebrisTSP.mu)\n",
        "          \n",
        "\n",
        "  # freezes the first element\n",
        "  def freeze_first(self):\n",
        "    self.index_frozen += 1\n",
        "\n",
        "  # returns if all debris were removed\n",
        "  def all_frozen(self):\n",
        "    return self.index_frozen >= (self.size-1-len(self.decayed_debris))\n",
        "\n",
        "  # transform the debris kepler elements to certain epoch\n",
        "  \"\"\"\n",
        "  dt_epoch : datetime : the target epoch\n",
        "  indexes : array : the debris that should be transformed\n",
        "  \"\"\"\n",
        "  def to_epoch(self, dt_epoch, indexes):\n",
        "    new_epoch = pykep.epoch_from_string(dt_epoch.strftime(FMT))\n",
        "    ranger = [x for x in range(0, self.size) if x in indexes]\n",
        "    \n",
        "    self.kepler_elements = list(numpy.zeros(self.size))\n",
        "    for j in ranger:\n",
        "      try:\n",
        "        self.kepler_elements[j] = debris[j].osculating_elements(new_epoch)\n",
        "      except:\n",
        "        self.decayed_debris.append(j)\n",
        "\n",
        "    for x in self.decayed_debris:\n",
        "      if x in ranger:\n",
        "        del ranger[ranger.index(x)]\n",
        "\n",
        "    for i in ranger:\n",
        "      for j in ranger:\n",
        "        if self.hohmanncost:\n",
        "          self.fitness_matrix[i][j] = StaticDebrisTSP.MYhohmann_impulse_aprox(self.kepler_elements[i], self.kepler_elements[j], new_epoch)\n",
        "        else:\n",
        "          try:\n",
        "            self.fitness_matrix[i][j] = pykep.phasing.three_impulses_approx(debris[i], debris[j], new_epoch, new_epoch)\n",
        "          except:\n",
        "            d1 = self.kepler_elements[i]\n",
        "            d2 = self.kepler_elements[j]\n",
        "            self.fitness_matrix[i][j] = StaticDebrisTSP.MYthree_impulse_aprox(d1[0],d1[1],d1[2],d1[3],d2[0],d2[1],d2[2],d2[3],StaticDebrisTSP.mu)\n",
        "\n",
        "    for x in self.decayed_debris:\n",
        "      if x in indexes:\n",
        "        del indexes[indexes.index(x)]\n",
        "\n",
        "    return indexes\n",
        "\n",
        "  # fitness is the sum cost to travel between each I and I+1 plus the last to initial\n",
        "  def fitness(self, solution):\n",
        "    fit = 0\n",
        "    for i in range(0, self.size-1):\n",
        "      fit += self.fitness_matrix[solution[i]][solution[i+1]]\n",
        "\n",
        "    fit += self.fitness_matrix[solution[self.size-1]][solution[0]]\n",
        "    return fit\n",
        "\n",
        "  # partial fitness is the sum cost to travel between each I and I+1\n",
        "  def partialFitness(self, part):\n",
        "    fit = 0\n",
        "    for i in range(0, len(part)-1):\n",
        "      fit += self.fitness_matrix[part[i]][part[i+1]]\n",
        "\n",
        "    return fit\n",
        "\n",
        "  # reward is the sum reward of the debris in the solution\n",
        "  def reward(self, solution):\n",
        "    reward = 0\n",
        "    for i in range(0, len(solution)):\n",
        "      reward += self.reward_matrix[solution[i]]\n",
        "\n",
        "    return reward\n",
        "\n",
        "  # estimate the duration of a solution\n",
        "  def duration(self, solution):\n",
        "    duration = 0\n",
        "    for i in range(0, len(solution)-1):\n",
        "      duration += self.transferDuration(solution[i], solution[i+1], StaticDebrisTSP.mu)\n",
        "\n",
        "    return duration\n",
        "\n",
        "  # fitness TD is the fitness function for a timedependent solution\n",
        "  def fitnessTD(self, solution):\n",
        "    if len(solution) < 2:\n",
        "      return 0\n",
        "      \n",
        "    fit = 0\n",
        "    for i in range(0, len(solution)-1):\n",
        "      epoch = pykep.epoch_from_string((solution[i+1][0]).strftime(FMT))\n",
        "      if self.hohmanncost:\n",
        "        d1 = debris[solution[i][1]].osculating_elements(epoch)\n",
        "        d2 = debris[solution[i+1][1]].osculating_elements(epoch)\n",
        "        fit += StaticDebrisTSP.MYhohmann_impulse_aprox(d1, d2, epoch)\n",
        "      else:\n",
        "        try:\n",
        "          fit += pykep.phasing.three_impulses_approx(debris[solution[i][1]], debris[solution[i+1][1]], epoch, epoch)\n",
        "        except:\n",
        "          d1 = debris[solution[i][1]].osculating_elements(epoch)\n",
        "          d2 = debris[solution[i+1][1]].osculating_elements(epoch)\n",
        "          fit += StaticDebrisTSP.MYthree_impulse_aprox(d1[0],d1[1],d1[2],d1[3],d2[0],d2[1],d2[2],d2[3],StaticDebrisTSP.mu)\n",
        "\n",
        "    return fit\n",
        "\n",
        "  # duration TD is the duration estimate for a timedependent solution\n",
        "  def durationTD(self, solution):\n",
        "    duration = 0\n",
        "    for i in range(0, len(solution)-1):\n",
        "      duration += (solution[i+1][0] - solution[i][0]).total_seconds() # seconds waiting for right epoch\n",
        "\n",
        "      epoch = pykep.epoch_from_string(solution[i+1][0].strftime(FMT))\n",
        "\n",
        "      duration += self.transferDurationTD(solution[i][1], solution[i+1][1], epoch, epoch, StaticDebrisTSP.mu)\n",
        "\n",
        "    return duration\n",
        "\n",
        "  # reward TD is the reward function for a timedependent solution\n",
        "  def rewardTD(self, solution):\n",
        "    reward = 0\n",
        "    for i in range(0, len(solution)):\n",
        "      reward += self.reward_matrix[solution[i][1]]\n",
        "\n",
        "    return reward\n",
        "  \n",
        "  # estimate the duration of a transfer (Hohmann) in seconds\n",
        "  def transferDuration(self, d1, d2, u):\n",
        "    d1_semi_major_axis = self.kepler_elements[d1][0]\n",
        "    d2_semi_major_axis = self.kepler_elements[d2][0]\n",
        "    transfer_semi_major_axis = (d1_semi_major_axis + d2_semi_major_axis) / 2\n",
        "    time_of_transfer = math.pi * math.sqrt((transfer_semi_major_axis**3) / u)\n",
        "    return time_of_transfer\n",
        "\n",
        "  # estimate the duration of a transfer (Hohmann) in seconds in a certain epoch\n",
        "  def transferDurationTD(self, d1, d2, epoch1, epoch2, u):\n",
        "    kepler1 = debris[d1].osculating_elements(epoch1)\n",
        "    kepler2 = debris[d2].osculating_elements(epoch2)\n",
        "    d1_semi_major_axis = kepler1[0]\n",
        "    d2_semi_major_axis = kepler2[0]\n",
        "    transfer_semi_major_axis = (d1_semi_major_axis + d2_semi_major_axis) / 2\n",
        "    time_of_transfer = math.pi * math.sqrt((transfer_semi_major_axis**3) / u)\n",
        "    return time_of_transfer\n",
        "\n",
        "  # find the constrained embedded maximal rewardable path in a solution\n",
        "  def maxOpenWalk(self, solution, cost_limit=1000, time_limit=31536000):\n",
        "    # calculate transferences\n",
        "    transfers = []\n",
        "    durations = []\n",
        "    for i in range(0, len(solution)-1):\n",
        "      sol_i = solution[i]\n",
        "      sol_j = solution[i+1]\n",
        "      transfers.append(self.fitness_matrix[sol_i][sol_j])\n",
        "      durations.append(self.transferDuration(sol_i, sol_j, StaticDebrisTSP.mu))\n",
        "\n",
        "    # calculate the maximal open walks starting at each arc\n",
        "    maxWalks = []\n",
        "    for i in range(0, len(transfers)):\n",
        "      cost = transfers[i]\n",
        "      duration = durations[i]\n",
        "      walk = [i]\n",
        "\n",
        "      for j in range(i+1, len(transfers)):\n",
        "        if (cost + transfers[j]) > cost_limit or (duration + durations[j]) > time_limit:\n",
        "          break;\n",
        "        else:\n",
        "          cost += transfers[j]\n",
        "          duration += durations[j]\n",
        "          walk.append(j)\n",
        "\n",
        "      nodes = []\n",
        "      reward = 0\n",
        "      for a in range(0, len(walk)):\n",
        "        arc = walk[a]\n",
        "        if solution[arc] not in nodes:\n",
        "          nodes.append(solution[arc])\n",
        "          reward += self.reward_matrix[solution[arc]]\n",
        "        nodes.append(solution[arc+1])\n",
        "        reward += self.reward_matrix[solution[arc+1]]\n",
        "\n",
        "      maxWalks.append({'walk': nodes, 'cost': cost, 'duration': duration, 'reward': reward})\n",
        "\n",
        "    # find the biggest open walk\n",
        "    w = 0\n",
        "    for i in range(1, len(maxWalks)):\n",
        "      if maxWalks[i]['reward'] > maxWalks[w]['reward']:\n",
        "        w = i\n",
        "\n",
        "    return maxWalks[w]\n",
        "\n",
        "  # find the constrained embedded maximal rewardable path in a timedependent solution\n",
        "  def maxOpenWalkTD(self, solution, cost_limit=1000, time_limit=31536000):\n",
        "    # calculate transferences\n",
        "    transfers = []\n",
        "    durations = []\n",
        "    for i in range(0, len(solution)-1):\n",
        "      epoch = pykep.epoch_from_string((solution[i+1][0]).strftime(FMT))\n",
        "      sol_i = solution[i][1]\n",
        "      sol_j = solution[i+1][1]\n",
        "\n",
        "      duration = (solution[i+1][0] - solution[i][0]).total_seconds() # seconds waiting for right epoch\n",
        "      duration += self.transferDurationTD(sol_i, sol_j, epoch, epoch, StaticDebrisTSP.mu)\n",
        "      durations.append(duration)\n",
        "\n",
        "      if self.hohmanncost:\n",
        "        d1 = debris[sol_i].osculating_elements(epoch)\n",
        "        d2 = debris[sol_j].osculating_elements(epoch)\n",
        "        transfers.append(StaticDebrisTSP.MYhohmann_impulse_aprox(d1, d2, epoch))\n",
        "      else:\n",
        "        try:\n",
        "          transfers.append(pykep.phasing.three_impulses_approx(debris[sol_i], debris[sol_j], epoch, epoch))\n",
        "        except:\n",
        "          d1 = debris[sol_i].osculating_elements(epoch)\n",
        "          d2 = debris[sol_j].osculating_elements(epoch)\n",
        "          transfers.append(StaticDebrisTSP.MYthree_impulse_aprox(d1[0],d1[1],d1[2],d1[3],d2[0],d2[1],d2[2],d2[3],StaticDebrisTSP.mu))\n",
        "\n",
        "    # calculate the maximal open walks starting at each arc\n",
        "    maxWalks = []\n",
        "    for i in range(0, len(transfers)):\n",
        "      cost = transfers[i]\n",
        "      duration = durations[i]\n",
        "      walk = [i]\n",
        "\n",
        "      for j in range(i+1, len(transfers)):\n",
        "        if (cost + transfers[j]) > cost_limit or (duration + durations[j]) > time_limit:\n",
        "          break;\n",
        "        else:\n",
        "          cost += transfers[j]\n",
        "          duration += durations[j]\n",
        "          walk.append(j)\n",
        "\n",
        "      nodes = []\n",
        "      reward = 0\n",
        "      for a in range(0, len(walk)):\n",
        "        arc = walk[a]\n",
        "        if solution[arc] not in nodes:\n",
        "          nodes.append(solution[arc])\n",
        "          reward += self.reward_matrix[solution[arc][1]]\n",
        "        nodes.append(solution[arc+1])\n",
        "        reward += self.reward_matrix[solution[arc+1][1]]\n",
        "\n",
        "      maxWalks.append({'walk': nodes, 'cost': cost, 'duration': duration, 'reward': reward})\n",
        "\n",
        "    # find the biggest open walk\n",
        "    w = 0\n",
        "    for i in range(1, len(maxWalks)):\n",
        "      if maxWalks[i]['reward'] > maxWalks[w]['reward']:\n",
        "        w = i\n",
        "\n",
        "    return maxWalks[w]\n",
        "\n",
        "  # estimate the hohmann cost for a transfer between two debris\n",
        "  # kepler elements order: a,e,i,W,w,M\n",
        "  def MYhohmann_impulse_aprox(kepler1, kepler2):\n",
        "    if kepler1 == kepler2:\n",
        "      return 0\n",
        "    \n",
        "    d1 = math.sqrt(StaticDebrisTSP.mu/kepler1[0]) * (math.sqrt((2*kepler2[0]) / (kepler1[0]+kepler2[0])) - 1)\n",
        "    d2 = math.sqrt(StaticDebrisTSP.mu/kepler2[0]) * (- math.sqrt((2*kepler1[0]) / (kepler1[0]+kepler2[0])) + 1)\n",
        "    dv = abs(d1 + d2)\n",
        "\n",
        "    re = - StaticDebrisTSP.mu / (2 * (StaticDebrisTSP.re + kepler2[0]))\n",
        "    rvi = math.sqrt(2 * ( (StaticDebrisTSP.mu / (StaticDebrisTSP.re + kepler2[0])) + re))\n",
        "    romega = abs(math.degrees(kepler2[2]) - math.degrees(kepler1[2]))\n",
        "    rdv = 2 * rvi * math.sin(romega/2)\n",
        "    \n",
        "    return abs(dv + rdv)\n",
        "\n",
        "  # estimate the edelbaum cost for a transfer between two debris\n",
        "  # this implementation replaces the pykep implementation, since pykep throws an exception for decayed debris\n",
        "  def MYthree_impulse_aprox(a1, e1, i1, W1, a2, e2, i2, W2, mu):\n",
        "    # radius of apocenter/pericenter starting orbit (ms)\n",
        "    ra1 = a1 * (1 + e1);\n",
        "    ra2 = a2 * (1 + e2);\n",
        "    rp1 = a1 * (1 - e2);\n",
        "    rp2 = a2 * (1 - e2);\n",
        "\n",
        "    # relative inclination between orbits\n",
        "    cosiREL = math.cos(i1) * math.cos(i2) + math.sin(i1) * math.sin(i2) * math.cos(W1 - W2)\n",
        "\n",
        "    # Strategy is Apocenter-Pericenter\n",
        "    if ra1 > ra2:\n",
        "      Vi = math.sqrt(mu * (2.0 / ra1 - 1.0 / a1));\n",
        "      Vf = math.sqrt(mu * (2.0 / ra1 - 2.0 / (rp2 + ra1)));\n",
        "\n",
        "      # Change Inclination + pericenter change\n",
        "      DV1 = math.sqrt(Vi * Vi + Vf * Vf - 2.0 * Vi * Vf * cosiREL);\n",
        "      # Apocenter Change\n",
        "      DV2 = math.sqrt(mu) * abs(math.sqrt(2.0 / rp2 - 2.0 / (rp2 + ra1)) - math.sqrt(2.0 / rp2 - 1.0 / a2));\n",
        "\n",
        "      return (DV1 + DV2)\n",
        "    \n",
        "    # Strategy is Pericenter-Apocenter\n",
        "    else:\n",
        "      Vi = math.sqrt(mu * ((2 / ra2) - (2 / (rp1 + ra2))));\n",
        "      Vf = math.sqrt(mu * ((2 / ra2) - (1 / a2)));\n",
        "\n",
        "      # Apocenter Raise\n",
        "      DV1 = math.sqrt(mu) * abs(math.sqrt((2 / rp1) - (2 / (rp1 + ra1))) - math.sqrt((2 / rp1) - (2 / (rp1 + ra2))));\n",
        "      # Change Inclination + apocenter change\n",
        "      DV2 = math.sqrt(abs((Vi * Vi) + (Vf * Vf) - (2 * Vi * Vf * cosiREL)));\n",
        "\n",
        "      return (DV1 + DV2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QW4852RdtwFn"
      },
      "source": [
        "# Instance loading\n",
        "The instances can be downloaded at SATCAT site.\n",
        "\n",
        "It is necessary to use a TXT file (TLE file) to get the debris names, codes and kepler elements, and a CSV file for the debris RCS (reward)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "awdRHHw729eV"
      },
      "source": [
        "deb_file = 'fengyun-1c-debris'\n",
        "debris = pykep.util.read_tle(tle_file=deb_file+'.txt', with_name=True)\n",
        "with open(deb_file+'.txt') as f:\n",
        "    tle_string = ''.join(f.readlines())\n",
        "\n",
        "tle_lines = tle_string.strip().splitlines()\n",
        "tle_elements = [tle_lines[i:i + 3] for i in range(0, len(tle_lines), 3)] #split in array of debris\n",
        "debris_tle = [TLE.from_lines(*tle_elements[i]) for i in range(0, len(tle_elements))]\n",
        "\n",
        "with open(deb_file+'.csv', newline='') as csvfile:\n",
        "    satcat = list(csv.reader(csvfile))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6grTld80upoA"
      },
      "source": [
        "# extract the reward for each debris\n",
        "areaDebris = []\n",
        "norad_index = satcat[0].index('NORAD_CAT_ID')\n",
        "rcs_index = satcat[0].index('RCS')\n",
        "for i in range(0, len(debris)):\n",
        "  rcs = 0\n",
        "  for j in range(1, len(satcat)):\n",
        "    if (debris_tle[i].norad == satcat[j][norad_index]):\n",
        "      if (satcat[j][rcs_index]):\n",
        "        rcs = float(satcat[j][rcs_index])\n",
        "      break\n",
        "  areaDebris.append(rcs)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZJHZk-KRmRWV"
      },
      "source": [
        "# Solution\n",
        "Here the actual solution is generated.\n",
        "\n",
        "An interpolated tree search is performed to enhance the static to a time dependent solution."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KgQqfRK5UfXr"
      },
      "source": [
        "start_epoch = \"2021-06-11 00:06:09\"\n",
        "FMT = '%Y-%m-%d %H:%M:%S'\n",
        "steps = int((24 * 60) / 10) * 7 # in days\n",
        "step_size = timedelta(minutes=10)\n",
        "removal_time = timedelta(days=1) # time taken to deorbit a debris\n",
        "winsize = 10 # range for the kopt\n",
        "\n",
        "for _ in range(10):\n",
        "  t0 = datetime.datetime.now() # to track time elapsed\n",
        "  epoch = datetime.datetime.strptime(start_epoch, FMT)\n",
        "\n",
        "  # generate the ga and problem instance\n",
        "  problem = StaticDebrisTSP(epoch=pykep.epoch_from_string(start_epoch), hohmanncost=False, debris=debris, reward_matrix=areaDebris, path_size=size, population_size=100)\n",
        "  ga = GA(population=problem.population, fn_fitness=problem.fitness, subpath_fn_fitness=problem.partialFitness)\n",
        "\n",
        "  # generate the static solution\n",
        "  curr_solution = ga.run_inverover(tinv=20000, feach=1000, runkopt=100, forn=5)\n",
        "  curr_fit = problem.partialFitness(curr_solution)\n",
        "  print('initial fit: '+str(curr_fit))\n",
        "\n",
        "  # find the static max open walk\n",
        "  path = problem.maxOpenWalk(curr_solution, 1000, 60*60*24*365) # 1km/s and 1 year\n",
        "  \n",
        "  # make the population start by best starting debris, and get the best then\n",
        "  ga.startBy(path['walk'][0])\n",
        "  curr_solution = ga.getBest()\n",
        "  curr_fit = problem.partialFitness(curr_solution)\n",
        "  print('secondal fit: '+str(curr_fit))\n",
        "\n",
        "  # use the first debris for the time dependent solution\n",
        "  solution = [(epoch, curr_solution[0])]\n",
        "  problem.freeze_first()\n",
        "  ga.freeze_first(curr_solution[0])\n",
        "\n",
        "  while not problem.all_frozen():\n",
        "    i = problem.index_frozen\n",
        "    \n",
        "    # run ranged kopt to optimize the current part of solution\n",
        "    if i > 0 and (i < len(curr_solution)-1):\n",
        "      curr_solution[i:i+winsize+1] = problem.to_epoch(epoch, curr_solution[i:i+winsize+1])\n",
        "      curr_solution[-(winsize+1):] = problem.to_epoch(epoch, curr_solution[-(winsize+1):])\n",
        "      ga.decay(problem.decayed_debris)\n",
        "      curr_solution = ga.ranged2opt(curr_solution, winsize)\n",
        "    \n",
        "    # get the next transference to be performed\n",
        "    transition = curr_solution[i:i+2]\n",
        "    \n",
        "    # validates if the debris in this transference are going to decay in during the interpolation\n",
        "    transition = problem.to_epoch(epoch + (step_size * steps), transition)\n",
        "    if len(transition) < 2:\n",
        "      curr_solution[i:i+2] = transition\n",
        "      ga.decay(problem.decayed_debris)\n",
        "      continue\n",
        "\n",
        "    # calculate the costs of the transference for the interpolation range\n",
        "    epoch_aux = epoch\n",
        "    x = []\n",
        "    y = []\n",
        "    for j in range(0, steps):\n",
        "      problem.to_epoch(epoch, transition)\n",
        "      x.append(j)\n",
        "      y.append(problem.partialFitness(transition))\n",
        "      epoch += step_size\n",
        "\n",
        "    # get the minimal cost point in the interpolated function\n",
        "    interpolator = scipy.interpolate.interp1d(x, y, kind='cubic')\n",
        "    xnew = numpy.linspace(0, steps-1, num=steps*3, endpoint=True) # num = precision\n",
        "    least = numpy.argmin(interpolator(xnew))\n",
        "\n",
        "    # get the epoch of the minimal cost transference\n",
        "    epoch = epoch_aux + (step_size * xnew[least])\n",
        "\n",
        "    # append to the time dependent solution\n",
        "    solution.append((epoch, curr_solution[i+1]))\n",
        "\n",
        "    # pushes the current epoch to after the deorbitation process\n",
        "    pykep_epoch = pykep.epoch_from_string(epoch.strftime(FMT))\n",
        "    transfer_duration = timedelta(seconds=problem.transferDurationTD(curr_solution[i], curr_solution[i+1], pykep_epoch, pykep_epoch, StaticDebrisTSP.mu))\n",
        "    epoch += removal_time + transfer_duration\n",
        "\n",
        "    # freezes the deorbited debris\n",
        "    problem.freeze_first()\n",
        "    ga.freeze_first(curr_solution[i+1])\n",
        "\n",
        "  t1 = datetime.datetime.now()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "awmaTZY23iR4"
      },
      "source": [
        "# instance results\n",
        "print(solution)\n",
        "print('elapsed time: '+ str(t1 - t0))\n",
        "print('fit: ' + str(problem.fitnessTD(solution)))\n",
        "print('dur: ' + str(problem.durationTD(solution)/60/60/24) + ' days')\n",
        "print('rew: ' + str(problem.rewardTD(solution)))\n",
        "\n",
        "# constrained (best mission) results\n",
        "path = problem.maxOpenWalkTD(solution, 1000, 60*60*24*365) # 1km/s and 1 year\n",
        "print(path)\n",
        "print('walk ' + str(len(path['walk'])))\n",
        "print('w_cost ' + str(path['cost']))\n",
        "print('w_rew ' + str(path['reward']))\n",
        "print('w_dur ' + str(path['duration']/60/60/24) + ' days')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GmwO7F7EhbmV"
      },
      "source": [
        "# Bibliography"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BH_K5qBMooYL"
      },
      "source": [
        "**Instances**\n",
        "\n",
        "TLE Derbis: https://celestrak.com/NORAD/elements/\n",
        "\n",
        "RCS: https://celestrak.com/satcat/search.php - LEGACY text\n",
        "\n",
        "Format: https://celestrak.com/satcat/satcat-format.php\n",
        "\n",
        "**Used Libs**\n",
        "\n",
        "https://esa.github.io/pykep/documentation/phasing.html#pykep.phasing.three_impulses_approx\n",
        "\n",
        "**Reference source codes**\n",
        "\n",
        "https://github.com/esa/pagmo/blob/master/src/problem/base_tsp.cpp\n",
        "https://github.com/esa/pagmo/blob/master/src/algorithm/inverover.cpp\n",
        "\n",
        "https://stackoverflow.com/questions/47982604/hamiltonian-path-using-python/47985349\n",
        "https://github.com/esa/pagmo/blob/80281d549c8f1b470e1489a5d37c8f06b2e429c0/src/util/neighbourhood.cpp\n",
        "\n",
        "https://github.com/esa/pagmo/blob/80281d549c8f1b470e1489a5d37c8f06b2e429c0/PyGMO/problem/_tsp.py\n",
        "https://github.com/esa/pagmo/blob/80281d549c8f1b470e1489a5d37c8f06b2e429c0/src/problem/base_tsp.cpp\n",
        "https://github.com/esa/pagmo/blob/80281d549c8f1b470e1489a5d37c8f06b2e429c0/src/problem/tsp.cpp\n",
        "https://github.com/esa/pagmo/blob/80281d549c8f1b470e1489a5d37c8f06b2e429c0/src/problem/tsp_cs.cpp\n",
        "https://github.com/esa/pagmo/blob/80281d549c8f1b470e1489a5d37c8f06b2e429c0/src/problem/tsp_ds.cpp\n",
        "https://github.com/esa/pykep/blob/2e1c97bea138d2c125d6695e7662991e6da30203/include/keplerian_toolbox/core_functions/three_impulses_approximation.hpp\n",
        "\n",
        "**Reference physics**\n",
        "\n",
        "https://en.wikipedia.org/wiki/Hohmann_transfer_orbit\n",
        "\n",
        "https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion#Third_law\n",
        "\n",
        "https://en.wikipedia.org/wiki/Orbital_period\n",
        "\n",
        "https://space.stackexchange.com/questions/35166/how-to-find-t₀-and-other-parameters-from-a-tle-to-calculate-an-approximate-mean/35190#35190\n",
        "\n",
        "https://space.stackexchange.com/questions/18289/how-to-get-semi-major-axis-from-tle\n",
        "\n",
        "https://ai-solutions.com/_freeflyeruniversityguide/hohmann_transfer.htm"
      ]
    }
  ]
}