Feature/15 add time objective

WeatherRoutingTool/algorithms/genetic/__init__.py

@@ -6,14 +6,18 @@
 import cartopy.crs as ccrs
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 from astropy import units as u
+from matplotlib.ticker import ScalarFormatter
 from pymoo.algorithms.moo.nsga2 import NSGA2
 from pymoo.core.result import Result
+from pymoo.optimize import minimize
 from pymoo.termination import get_termination
 from pymoo.util.running_metric import RunningMetric
 
 import WeatherRoutingTool.utils.formatting as formatting
 import WeatherRoutingTool.utils.graphics as graphics
+import WeatherRoutingTool.algorithms.genetic.mcdm as MCDM
 from WeatherRoutingTool.algorithms.genetic.population import PopulationFactory
 from WeatherRoutingTool.algorithms.genetic.crossover import CrossoverFactory
 from WeatherRoutingTool.algorithms.genetic.mutation import MutationFactory
@@ -48,6 +52,8 @@ def __init__(self, config: Config):
 
         self.n_generations = config.GENETIC_NUMBER_GENERATIONS
         self.n_offsprings = config.GENETIC_NUMBER_OFFSPRINGS
+        self.objectives = config.GENETIC_OBJECTIVES
+        self.n_objs = len(config.GENETIC_OBJECTIVES)
 
         # population
         self.pop_type = config.GENETIC_POPULATION_TYPE
@@ -83,7 +89,9 @@ def execute_routing(
             arrival_time=self.arrival_time,
             boat_speed=self.boat_speed,
             boat=boat,
-            constraint_list=constraints_list, )
+            constraint_list=constraints_list,
+            objectives=self.objectives
+        )
 
         initial_population = PopulationFactory.get_population(
             self.config, boat, constraints_list, wt, )
@@ -152,16 +160,35 @@ def optimize(
 
         return res
 
+    # FIXME temporary consistency check
+    def consistency_check(self, res, problem):
+        X = res.X
+        i_route = 0
+        for route in X:
+            fuel_dict = problem.get_power(route[0])
+
+            i_obj = 0
+            for obj_str in self.objectives:
+                if obj_str == "fuel_consumption":
+                    np.testing.assert_equal(fuel_dict["fuel_sum"].value, res.F[i_route, i_obj], 5)
+                else:
+                    np.testing.assert_equal(fuel_dict["time_obj"], res.F[i_route, i_obj], 5)
+                i_obj += 1
+            i_route += 1
+
     def terminate(self, res: Result, problem: RoutingProblem):
         """Genetic Algorithm termination procedures"""
 
         super().terminate()
+        self.consistency_check(res, problem)
 
-        best_index = res.F.argmin()
-        # ensure res.X is of shape (n_sol, n_var)
+        mcdm = MCDM.RMethod(self.objectives)
+        best_index = mcdm.get_best_compromise(res.F)
         best_route = np.atleast_2d(res.X)[best_index, 0]
 
-        fuel, ship_params = problem.get_power(best_route)
+        fuel_dict = problem.get_power(best_route)
+        fuel = fuel_dict["fuel_sum"]
+        ship_params = fuel_dict["shipparams"]
         logger.info(f"Best fuel: {fuel}")
 
         if self.figure_path is not None:
@@ -171,6 +198,7 @@ def terminate(self, res: Result, problem: RoutingProblem):
             self.plot_population_per_generation(res, best_route)
             self.plot_convergence(res)
             self.plot_coverage(res, best_route)
+            self.plot_objective_space(res, best_index)
 
         lats = best_route[:, 0]
         lons = best_route[:, 1]
@@ -212,6 +240,32 @@ def terminate(self, res: Result, problem: RoutingProblem):
         self.check_positive_power()
         return route
 
+    def plot_objective_space(self, res, best_index):
+        F = res.F
+        fig, ax = plt.subplots(figsize=(7, 5))
+
+        if self.n_objs == 2:
+            ax.scatter(F[:, 0], F[:, 1], s=30, facecolors='none', edgecolors='blue')
+        else:
+            return
+
+        ax.plot(F[best_index, 0], F[best_index, 1], color='red', marker='o')
+        ax.set_xlabel('f1', labelpad=10)
+        ax.set_ylabel('f2', labelpad=10)
+        ax.grid(True, linestyle='--', alpha=0.7)
+        plt.title("Objective Space")
+
+        formatter = ScalarFormatter(useMathText=True)
+        formatter.set_scientific(True)
+        formatter.set_powerlimits((-1, 1))  # Force scientific notation
+
+        ax.xaxis.set_major_formatter(formatter)
+        ax.yaxis.set_major_formatter(formatter)
+
+        plt.savefig(os.path.join(self.figure_path, 'genetic_objective_space.png'))
+        plt.cla()
+        plt.close()
+
     def print_init(self):
         """Log messages to print on algorithm initialization"""
 
@@ -300,6 +354,7 @@ def plot_population_per_generation(self, res, best_route):
         input_crs = ccrs.PlateCarree()
         history = res.history
         fig, ax = plt.subplots(figsize=graphics.get_standard('fig_size'))
+        route_lc = None
 
         for igen in range(len(history)):
             plt.rcParams['font.size'] = graphics.get_standard('font_size')
@@ -328,19 +383,24 @@ def plot_population_per_generation(self, res, best_route):
                     ax.plot(
                         last_pop[iroute, 0][:, 1],
                         last_pop[iroute, 0][:, 0],
-                        **(marker_kw if igen != self.n_generations - 1 else {}),
+                        # **(marker_kw if igen != self.n_generations - 1 else {}),
                         color="firebrick",
                         label=f"full population [{last_pop.shape[0]}]",
+                        linewidth=0,
                         transform=input_crs)
 
                 else:
                     ax.plot(
                         last_pop[iroute, 0][:, 1],
                         last_pop[iroute, 0][:, 0],
-                        **(marker_kw if igen != self.n_generations - 1 else {}),
+                        # **(marker_kw if igen != self.n_generations - 1 else {}),
                         color="firebrick",
+                        linewidth=0,
                         transform=input_crs)
 
+                route_lc = graphics.get_route_lc(last_pop[iroute, 0])
+                ax.add_collection(route_lc)
+
             if igen == (self.n_generations - 1):
                 ax.plot(
                     best_route[:, 1],
@@ -350,6 +410,9 @@ def plot_population_per_generation(self, res, best_route):
                     label="best route",
                     transform=input_crs
                 )
+            cbar = fig.colorbar(route_lc, ax=ax, orientation='vertical', pad=0.15, shrink=0.7)
+            cbar.set_label('Geschwindigkeit ($m/s$)')
+            plt.tight_layout()
 
             ax.legend()
 
@@ -387,27 +450,33 @@ def plot_convergence(self, res):
         """Plot the convergence curve (best objective value per generation)."""
 
         best_f = []
+        is_initialised = False
 
         for algorithm in res.history:
-            # For single-objective, take min of F; for multi-objective, take min of first objective
             F = algorithm.pop.get('F')
-            if F.ndim == 2:
-                best_f.append(np.min(F[:, 0]))
-            else:
-                best_f.append(np.min(F))
 
-        n_gen = np.arange(1, len(best_f) + 1)
+            for iobj in range(self.n_objs):
+                if not is_initialised:
+                    best_f.append([])
+                best_f[iobj].append(np.min(F[:, iobj]))
+            is_initialised = True
 
-        # plot png
-        plt.figure(figsize=graphics.get_standard('fig_size'))
-        plt.plot(n_gen, best_f, marker='o')
-        plt.xlabel('Generation')
-        plt.ylabel('Best Objective Value')
-        plt.title('Convergence Plot')
-        plt.grid(True)
-        plt.savefig(os.path.join(self.figure_path, 'genetic_algorithm_convergence.png'))
-        plt.cla()
-        plt.close()
+        n_gen = np.arange(1, len(best_f[0]) + 1)
 
-        # write to csv
-        graphics.write_graph_to_csv(os.path.join(self.figure_path, 'genetic_algorithm_convergence.csv'), n_gen, best_f)
+        # plot png
+        i_obj = 0
+        for obj_str in self.objectives:
+            fig_path_name = 'genetic_algorithm_convergence' + obj_str
+            plt.figure(figsize=graphics.get_standard('fig_size'))
+            plt.plot(n_gen, best_f[i_obj], marker='o')
+            plt.xlabel('Generation')
+            plt.ylabel('Best Objective Value ' + obj_str)
+            plt.title('Convergence Plot')
+            plt.grid(True)
+            plt.savefig(os.path.join(self.figure_path, fig_path_name + '.png'))
+            plt.cla()
+            plt.close()
+
+            # write to csv
+            graphics.write_graph_to_csv(os.path.join(self.figure_path, fig_path_name + '.csv'), n_gen, best_f[i_obj])
+            i_obj += 1

WeatherRoutingTool/algorithms/genetic/crossover.py

@@ -249,10 +249,6 @@ class SpeedCrossover(OffspringRejectionCrossover):
     """
 
     def __init__(self, **kw):
-        # for now, we don't want to allow repairing routes for speed crossover
-        config = deepcopy(kw['config'])
-        config.GENETIC_REPAIR_TYPE = ["no_repair"]
-        kw['config'] = config
         super().__init__(**kw)
         self.threshold = 50000  # in m
         self.percentage = 0.5
@@ -262,24 +258,27 @@ def crossover(
             p1: np.ndarray,
             p2: np.ndarray
     ) -> tuple[np.ndarray, np.ndarray]:
+        o1 = deepcopy(p1)
+        o2 = deepcopy(p2)
+
         # Find points between parents with a distance below the specified threshold.
         # There should always be one candidate (source). The destination has to be ignored.
         crossover_candidates = []
-        for m in range(0, len(p1)-1):
-            coord1 = p1[m, 0:2]
-            for n in range(0, len(p2)-1):
-                coord2 = p2[n, 0:2]
+        for m in range(0, len(o1) - 1):
+            coord1 = o1[m, 0:2]
+            for n in range(0, len(o2) - 1):
+                coord2 = o2[n, 0:2]
                 d = geod.Inverse(coord1[0], coord1[1], coord2[0], coord2[1])["s12"]
                 if d < self.threshold:
                     crossover_candidates.append((m, n))
         # Swap speed values for a subset of candidate points
-        indices = random.sample(range(0, len(crossover_candidates)), ceil(self.percentage*len(crossover_candidates)))
+        indices = random.sample(range(0, len(crossover_candidates)), ceil(self.percentage * len(crossover_candidates)))
         for idx in indices:
             c = crossover_candidates[idx]
-            speed1 = p1[c[0], -1]
-            p1[c[0], -1] = p2[c[1], -1]
-            p2[c[1], -1] = speed1
-        return p1, p2
+            speed1 = o1[c[0], -1]
+            o1[c[0], -1] = o2[c[1], -1]
+            o2[c[1], -1] = speed1
+        return o1, o2
 
 
 # factory
@@ -300,7 +299,7 @@ def get_crossover(config: Config, constraints_list: ConstraintsList):
                 prob=.5,
                 crossover_type="Speed crossover")
 
-        if config.GENETIC_CROSSOVER_TYPE == "random":
+        if config.GENETIC_CROSSOVER_TYPE == "waypoints":
             logger.debug('Setting crossover type of genetic algorithm to "random".')
             return RandomizedCrossoversOrchestrator(
                 opts=[
@@ -319,3 +318,29 @@ def get_crossover(config: Config, constraints_list: ConstraintsList):
                         prob=.5,
                         crossover_type="SP crossover")
                 ])
+
+        if config.GENETIC_CROSSOVER_TYPE == "random":
+            logger.debug('Setting crossover type of genetic algorithm to "random".')
+            return RandomizedCrossoversOrchestrator(
+                opts=[
+                    TwoPointCrossover(
+                        config=config,
+                        patch_type=config.GENETIC_CROSSOVER_PATCHER + "_singleton",
+                        departure_time=departure_time,
+                        constraints_list=constraints_list,
+                        prob=.5,
+                        crossover_type="TP crossover"),
+                    SinglePointCrossover(
+                        config=config,
+                        patch_type=config.GENETIC_CROSSOVER_PATCHER + "_singleton",
+                        departure_time=departure_time,
+                        constraints_list=constraints_list,
+                        prob=.5,
+                        crossover_type="SP crossover"),
+                    SpeedCrossover(
+                        config=config,
+                        departure_time=departure_time,
+                        constraints_list=constraints_list,
+                        prob=.5,
+                        crossover_type="Speed crossover")
+                ])