Enable speed optimisation for a fixed route

WeatherRoutingTool/algorithms/genetic/__init__.py

@@ -162,18 +162,30 @@ def optimize(
 
     # FIXME temporary consistency check
     def consistency_check(self, res, problem):
+        """
+        Temporary consistency check to uncover memory issues.
+        """
         X = res.X
+        res_objs = res.F
         i_route = 0
+
+        # solve shape issue in case there is only one objective
+        if self.n_objs == 1:
+            res_objs = np.array([res.F])
+            X = [X]
+
         for route in X:
             fuel_dict = problem.get_power(route[0])
 
-            i_obj = 0
+            # ordering of objective values in res.F is defined by RoutingProblem.get_objectives()
             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)
+                    i_obj = 1
+                    if self.n_objs == 1:
+                        i_obj = 0
+                    np.testing.assert_equal(fuel_dict["fuel_sum"].value, res_objs[i_route, i_obj], 5)
                 else:
-                    np.testing.assert_equal(fuel_dict["time_obj"], res.F[i_route, i_obj], 5)
-                i_obj += 1
+                    np.testing.assert_equal(fuel_dict["time_obj"], res_objs[i_route, 0], 5)
             i_route += 1
 
     def terminate(self, res: Result, problem: RoutingProblem):
@@ -196,6 +208,7 @@ def terminate(self, res: Result, problem: RoutingProblem):
 
             self.plot_running_metric(res)
             self.plot_population_per_generation(res, best_route)
+            self.plot_speed_per_generation(res, best_route)
             self.plot_convergence(res)
             self.plot_coverage(res, best_route)
             self.plot_objective_space(res, best_index)
@@ -343,8 +356,62 @@ def plot_running_metric(self, res):
         plt.cla()
         plt.close()
 
+    def plot_speed_per_generation(self, res, best_route) -> None:
+        """Plot line diagrams of speed vs. travel distance for each individual in one generation.
+
+        :param res: Result of GA minimization
+        :type res: pymoo.core.result.Result
+        :param best_route: Optimum route
+        :type best_route: np.ndarray
+        """
+        history = res.history
+
+        for igen in range(len(history)):
+            plt.clf()
+            plt.close('all')
+
+            fig, ax = plt.subplots(figsize=graphics.get_standard('fig_size'))
+            plt.rcParams['font.size'] = graphics.get_standard('font_size')
+
+            last_pop = history[igen].pop.get('X')
+            objs = []
+            for iroute in range(0, last_pop.shape[0]):
+                hist_values = utils.get_hist_values_from_route(last_pop[iroute, 0], self.departure_time)
+
+                new_line = ax.plot(
+                    hist_values["bin_centres"].to(u.km).value,
+                    hist_values["bin_contents"].to(u.m / u.second).value,
+                    color="blue",
+                    alpha=0.3,
+                    linestyle='-',
+                    zorder=2
+                )
+                objs.append(new_line)
+
+            if igen == (self.n_generations - 1):
+                hist_values_best_route = utils.get_hist_values_from_route(best_route, self.departure_time)
+                ax.plot(
+                    hist_values_best_route["bin_centres"].to(u.km).value,
+                    hist_values_best_route["bin_contents"].to(u.m / u.second).value,
+                    color="firebrick",
+                    linewidth=3
+                )
+            left, right = plt.xlim()
+            ax.set_xlim(-100, right)
+            ax.set_ylim(0, 10)
+
+            plt.ylabel("speed (m/s)")
+            plt.xlabel('travel distance (km)')
+            plt.xticks()
+            plt.tight_layout()
+            ax.legend()
+
+            figname = f"genetic_algorithm_speed {igen:02}.png"
+            plt.savefig(os.path.join(self.figure_path, figname))
+            plt.close(fig)
+
     def plot_population_per_generation(self, res, best_route):
-        """Plot figures and save them in WRT_FIGURE_PATH
+        """Plot routes for each individual in one generation on a map.
 
         :param res: Result of GA minimization
         :type res: pymoo.core.result.Result
@@ -413,7 +480,6 @@ def plot_population_per_generation(self, res, best_route):
             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()
 
             figname = f"genetic_algorithm_generation {igen:02}.png"

WeatherRoutingTool/algorithms/genetic/crossover.py

@@ -281,6 +281,55 @@ def crossover(
         return o1, o2
 
 
+class TwoPointCrossoverSpeed(OffspringRejectionCrossover):
+    """
+    Class for two-point crossover of ship speed.
+
+    The ship speed of a random sequence of one chromosome is replaced by the average ship speed of a random sequence
+    of another individual.
+    """
+
+    def __init__(self, **kw):
+        super().__init__(**kw)
+
+    def crossover(self, p1, p2) -> tuple[np.ndarray, np.ndarray]:
+        """
+        Crossover implementation for two-point crossover of ship speed.
+
+        :param p1: the first chromosome
+        :param p2: the second chromosome
+        :return: the two offspring chromosomes
+        :rtype: tuple[np.ndarray, np.ndarray]
+        """
+        r1 = deepcopy(p1)
+        r2 = deepcopy(p2)
+
+        p1x1 = np.random.randint(1, p1.shape[0] - 4)
+        p1x2 = p1x1 + np.random.randint(3, p1.shape[0] - p1x1 - 1)
+
+        p2x1 = np.random.randint(1, p2.shape[0] - 4)
+        p2x2 = p2x1 + np.random.randint(3, p2.shape[0] - p2x1 - 1)
+
+        speed1 = p1[:, -1]
+        speed2 = p2[:, -1]
+        av_speed_seg1 = np.average(speed1[p1x1:p1x2 + 1])
+        av_speed_seg2 = np.average(speed2[p2x1:p2x2 + 1])
+
+        new_speed1 = np.concatenate([
+            speed1[:p1x1],
+            np.full(p1x2 - p1x1, av_speed_seg2),
+            speed1[p1x2:], ])
+        new_speed2 = np.concatenate([
+            speed2[:p2x1],
+            np.full(p2x2 - p2x1, av_speed_seg1),
+            speed2[p2x2:], ])
+
+        r1[:, -1] = new_speed1
+        r2[:, -1] = new_speed2
+
+        return r1, r2
+
+
 # factory
 # ----------
 class CrossoverFactory:
@@ -292,12 +341,12 @@ def get_crossover(config: Config, constraints_list: ConstraintsList):
 
         if config.GENETIC_CROSSOVER_TYPE == "speed":
             logger.debug('Setting crossover type of genetic algorithm to "speed".')
-            return SpeedCrossover(
+            return TwoPointCrossoverSpeed(
                 config=config,
                 departure_time=departure_time,
                 constraints_list=constraints_list,
                 prob=.5,
-                crossover_type="Speed crossover")
+                crossover_type="TP Crossover speed")
 
         if config.GENETIC_CROSSOVER_TYPE == "waypoints":
             logger.debug('Setting crossover type of genetic algorithm to "random".')

WeatherRoutingTool/algorithms/genetic/mutation.py

@@ -147,6 +147,9 @@ class NoMutation(MutationBase):
     def _do(self, problem, X, **kw):
         return X
 
+    def print_mutation_statistics(self):
+        print('No mutation.')
+
 
 class RandomPlateauMutation(MutationConstraintRejection):
     """
@@ -574,7 +577,7 @@ class GaussianSpeedMutation(MutationConstraintRejection):
     n_updates: int
     config: Config
 
-    def __init__(self, n_updates: int = 10, **kw):
+    def __init__(self, n_updates: int = 5, **kw):
         super().__init__(
             mutation_type="GaussianSpeedMutation",
             **kw
@@ -583,7 +586,7 @@ def __init__(self, n_updates: int = 10, **kw):
         # FIXME: these numbers should be carefully evaluated
         # ~99.7 % in interval (0, BOAT_SPEED_MAX)
         self.mu = 0.5 * self.config.BOAT_SPEED_BOUNDARIES[1]
-        self.sigma = self.config.BOAT_SPEED_BOUNDARIES[1] / 6
+        self.sigma = 1.
 
     def mutate(self, problem, rt, **kw):
         rt_new = copy.deepcopy(rt)
@@ -600,6 +603,7 @@ def mutate(self, problem, rt, **kw):
                 new = old_speed
             rt_new[i][2] = new
 
+        rt_new[:, 2] = utils.smoothen_speed(rt_new[:, 2], 1)
         return rt_new
 
 

WeatherRoutingTool/algorithms/genetic/population.py

@@ -190,11 +190,53 @@ def __init__(self, config: Config, routes_dir: Path, default_route, constraints_
             constraints_list=constraints_list,
             pop_size=pop_size
         )
+        self.sole_speed_mutation = config.GENETIC_MUTATION_TYPE == "speed" and config.GENETIC_CROSSOVER_TYPE == "speed"
+        self.min_boat_speed = config.BOAT_SPEED_BOUNDARIES[0]
+        self.max_boat_speed = config.BOAT_SPEED_BOUNDARIES[1]
 
         if not routes_dir.exists() or not routes_dir.is_dir():
             raise FileNotFoundError("Routes directory not found")
         self.routes_dir = routes_dir
 
+    @staticmethod
+    def spread_velocity(min_boat_speed: float, max_boat_speed: float, boat_speed: float, pop_size: float) -> np.ndarray:
+        """
+        Calculate velocity spread for individuals in case of pure speed optimisation.
+
+        The following steps are performed to obtain values for the boat speed that are accumulated around the original
+        boat speed:
+            - sample a numpy array from a gaussian distribution (mean = original boat speed,
+              sigma = 1/4 possible value range of boat speed)
+            - cut values below the minimum velocity and above the maximum velocity
+            - determine `pop_size` quantiles with equally spaced probabilities
+
+        :param min_boat_speed: Minimum boat speed
+        :type min_boat_speed: float
+        :param max_boat_speed: Maximum boat speed
+        :type max_boat_speed: float
+        :param boat_speed: Boat speed
+        :type boat_speed: float
+        :param pop_size: Population size
+        :type pop_size: int
+        :return: array of quantiles
+        :rtype: np.ndarray
+        """
+        std_dev = (max_boat_speed - min_boat_speed) / 4
+        gaussian_sample = np.random.normal(boat_speed, std_dev, 1000)
+        gaussian_sample[gaussian_sample < min_boat_speed] = np.nan
+        gaussian_sample[gaussian_sample > max_boat_speed] = np.nan
+        gaussian_sample = gaussian_sample[~np.isnan(gaussian_sample)]
+
+        quant_size = 100. / pop_size * 0.01
+
+        quantiles = np.full(pop_size, np.nan)
+        quant_sum = 0
+        for q in range(pop_size):
+            quantiles[q] = np.quantile(gaussian_sample, q=quant_sum)
+            quant_sum += quant_size
+
+        return quantiles
+
     def generate(self, problem, n_samples, **kw):
         logger.debug(f"Population from geojson routes: {self.routes_dir}")
 
@@ -205,15 +247,22 @@ def generate(self, problem, n_samples, **kw):
         # FIXME: add test in config.py and raise exception depending on configuration (not only speed optimization...)
 
         X = np.full((n_samples, 1), None, dtype=object)
+        quantiles = None
+
+        # determine quantiles for pure speed optimisation
+        if self.sole_speed_mutation:
+            quantiles = self.spread_velocity(self.min_boat_speed, self.max_boat_speed, self.boat_speed.value,
+                                             self.pop_size)
 
+        # obtain list of filenames
         files = []
         for file in os.listdir(self.routes_dir):
             if match(r"route_[0-9]+\.(json|geojson)$", file.lower()):
                 files.append(file)
-
         if len(files) == 0:
             raise ValueError(f"Couldn't find any route in {self.routes_dir} for the initial population.")
 
+        # read route(s) from file
         for i, file in enumerate(files):
             path = os.path.join(self.routes_dir, file)
             if not os.path.exists(path):
@@ -227,6 +276,12 @@ def generate(self, problem, n_samples, **kw):
             X[added_routes, 0] = np.copy(X[0, 0])
             added_routes += 1
 
+        # mutate velocity in case of pure speed optimisation
+        if self.sole_speed_mutation:
+            for i, (rt,) in enumerate(X):
+                rt[:, -1] = quantiles[i]
+                rt[-1, -1] = -99.
+
         return X