diff --git a/WeatherRoutingTool/algorithms/genetic/__init__.py b/WeatherRoutingTool/algorithms/genetic/__init__.py index c514fb0a..6887d406 100644 --- a/WeatherRoutingTool/algorithms/genetic/__init__.py +++ b/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 diff --git a/WeatherRoutingTool/algorithms/genetic/crossover.py b/WeatherRoutingTool/algorithms/genetic/crossover.py index c525bf0f..8de38e72 100644 --- a/WeatherRoutingTool/algorithms/genetic/crossover.py +++ b/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") + ]) diff --git a/WeatherRoutingTool/algorithms/genetic/mcdm.py b/WeatherRoutingTool/algorithms/genetic/mcdm.py new file mode 100644 index 00000000..602fbe1f --- /dev/null +++ b/WeatherRoutingTool/algorithms/genetic/mcdm.py @@ -0,0 +1,229 @@ +import numpy as np +import pandas as pd + + +class MCDM: + """Base Class for Multi-Criteria Decision Making (MCDM). + + This Class implements the base functionality for selecting a single solution from the set of non-dominated + solutions found by the genetic algorithm. + """ + + def __init__(self, objectives: dict): + self.n_objs = len(objectives) + self.objectives = objectives + + def get_best_compromise(self, solutions): + pass + + +class RMethod(MCDM): + """ + Implements the R-Method for MCDM. + + The R-Method ranks alternatives by calculating a composite weight derived from + rank-based objective importance and solution performance rankings. It is + designed to identify the best compromise solution in multi-objective scenarios. + + If two objectives are supposed to have equal rank, the user needs to provide the mean between the two + possible ranks, e.g. weight=1.5 for both objectives for optimisation with two objectives. + + References: + R.V. Rao and R.J. Lakshmi, "Ranking of Pareto-optimal solutions and selecting the best solution in multi-and + many-objective optimization problems using R-method". Soft Computing Letters 3 (2021) 100015 + + :param objectives: dictionary of objective names (keys) and their weights (values). + :type objectives: dict + """ + + def __init__(self, objectives: dict[str, int]): + """ + Initialises the RMethod instance. + + :param objectives: dictionary of objective names (keys) and their weights (values). + :type objectives: dict + """ + super().__init__(objectives) + + self.objective_weights = {} + + for obj_str in self.objectives: + self.objective_weights[obj_str] = self.get_weigths_from_rankarr( + np.array([self.objectives[obj_str]]), + self.n_objs + ) + + def rank_solutions(self, obj: np.ndarray, dec: bool = False) -> np.ndarray: + """ + Ranks array content according to increasing (dec = False) or decreasing (dec = True) values. + + :param obj: Array to be ranked. + :type obj: numpy.ndarray + :param dec: If True, rank in descending order (highest value gets rank 1). + Defaults to False. + :type dec: bool + :return: array of ranks + :rtype: numpy.ndarray + """ + rank_ind = np.argsort(obj) + if dec: + rank_ind = rank_ind[::-1] + rank = np.argsort(rank_ind) + rank = rank + 1 + return rank + + def get_composite_weight(self, sol_weight_list: list[np.ndarray], obj_weight_list: list[float]) -> np.ndarray: + """ + Calculate the composite weight for all non-dominated solutions based on solution weights and objective weights. + + Note: The current implementation is limited to problems with exactly two objectives. + + :param sol_weight_list: List of weights for each solution based on their performance wrt. each objective. + :type sol_weight_list: list[numpy.ndarray] + :param obj_weight_list: List of objective weights based on user ranking. + :type obj_weight_list: list[float] + :raises NotImplementedError: If the number of objectives is greater than two. + :return: Calculated composite weights for all solutions. + :rtype: numpy.ndarray + """ + sign = [1, -1] + denominator = 0 + summands = 0 + product = 1 + + if len(sol_weight_list) > 2: + raise NotImplementedError('Calculation of the composite weight of the R-method is not implemented for' + 'more than two objectives.') + + for i in range(len(sol_weight_list)): + denominator = denominator + sign[i] * 1. / obj_weight_list[i] * sol_weight_list[i] + product = product * sol_weight_list[i] + denominator = np.abs(denominator) + 0.2 + + for i in range(len(sol_weight_list)): + summands = summands + sol_weight_list[i] / denominator * obj_weight_list[i] * obj_weight_list[i] + + composite_weight = product + summands + return composite_weight + + def get_best_compromise(self, solutions: np.ndarray) -> int: + """ + Find the index of the best compromise solution from a set of candidates. + + This method normalises the objective values for each solution, ranks the latter with respect to each objective, + and calculates composite weights based on solution and objective weights. + + :param solutions: 2D array of objective values where rows are alternative solutions and columns are + objective values. + :type solutions: numpy.ndarray + :return: The index of the optimal solution in the provided array. + :rtype: int + """ + debug = False + sol_weight_list = [] + obj_weight_list = [] + + if self.n_objs == 1: + return solutions.argmin() + + if debug: + print('solutions: ', solutions) + print('solutions shape: ', solutions.shape) + + rmethod_table = pd.DataFrame() + pd.set_option('display.max_rows', None) + pd.set_option('display.max_columns', None) + + i_obj = 0 + norm = 1. + for obj_str in self.objectives: + objective_values = solutions[:, i_obj] + max_value = np.max(objective_values) + if i_obj == 0: + norm = max_value + else: + objective_values = objective_values * norm * 1. / max_value + rmethod_table[obj_str + '_obj'] = objective_values + rmethod_table[obj_str + '_rank'] = self.rank_solutions(objective_values) + rmethod_table[obj_str + '_weight'] = self.get_weigths_from_rankarr( + rmethod_table[obj_str + '_rank'].to_numpy(), + len(solutions)) + i_obj += 1 + sol_weight_list.append(rmethod_table[obj_str + '_weight'].to_numpy()) + obj_weight_list.append(self.objective_weights[obj_str]) + + if debug: + print('rmethod table:', rmethod_table) + + rmethod_table['composite_weight'] = self.get_composite_weight( + sol_weight_list=sol_weight_list, + obj_weight_list=obj_weight_list, + ) + rmethod_table['composite_rank'] = self.rank_solutions(rmethod_table['composite_weight'], True) + best_ind = np.argmax(rmethod_table['composite_rank'].to_numpy()) + + if debug: + print('rmethod table:', rmethod_table) + print('best index: ', rmethod_table.iloc[best_ind]) + return best_ind + + def get_rank_sum(self, rank_max: int) -> float: + """ + Calculate the reciprocal of the harmonic sum for a given rank. + + :param rank_max: The rank to evaluate the sum for. + :type rank_max: int + :return: Reciprocal of the sum of (1/k) for k from 1 to rank_max. + :rtype: float + """ + rank_sum = 0 + for rk in range(1, rank_max + 1): + rank_sum += 1 / rk + return 1 / rank_sum + + def get_weight_from_rank(self, rank: int, n_parts: int) -> float: + """ + Compute a normalized weight for an individual rank of a solutions or objective. + + :param rank: The specific rank of the item. + :type rank: int + :param n_parts: Total number of solutions/objectives in the set. + :type n_parts: int + :return: Normalized weight. + :rtype: float + """ + numerator = self.get_rank_sum(rank) + denominator_sum = 0. + + for j in range(1, n_parts + 1): + temp = self.get_rank_sum(j) + denominator_sum += temp + return numerator / denominator_sum + + def get_weigths_from_rankarr(self, rank_arr: np.ndarray, n_parts: int) -> np.ndarray: + """ + Convert an array of ranks into an array of weights. + + Objectives or solutions can receive fractional ranks, if two of them are supposed to have equal ranks. These + fractional ranks are handled by averaging the weights of the floor and ceiling integer ranks. + + :param rank_arr: Array of ranks. + :type rank_arr: numpy.ndarray + :param n_parts: Total number of solutions/objectives in the set. + :type n_parts: int + :return: Array containing the derived weights. + :rtype: numpy.ndarray + """ + weight_array = np.full(rank_arr.shape, -99.) + + for irank in range(0, rank_arr.shape[0]): + if rank_arr[irank] % 1. != 0.: + smaller = int(np.floor(rank_arr[irank])) + larger = int(np.ceil(rank_arr[irank])) + smaller_weight = self.get_weight_from_rank(smaller, n_parts) + larger_weight = self.get_weight_from_rank(larger, n_parts) + weight_array[irank] = (smaller_weight + larger_weight) / 2 + else: + weight_array[irank] = self.get_weight_from_rank(int(rank_arr[irank]), n_parts) + + return weight_array diff --git a/WeatherRoutingTool/algorithms/genetic/mutation.py b/WeatherRoutingTool/algorithms/genetic/mutation.py index 6fe4ddf0..66b3a864 100644 --- a/WeatherRoutingTool/algorithms/genetic/mutation.py +++ b/WeatherRoutingTool/algorithms/genetic/mutation.py @@ -1,3 +1,4 @@ +import copy import logging import math import os @@ -435,20 +436,22 @@ def mutate(self, problem, rt, **kw): # test shape of input route assert len(rt.shape) == 2 assert rt.shape[1] == 3 - route_length = rt.shape[0] + rt_new = copy.deepcopy(rt) + + route_length = rt_new.shape[0] # only mutate routes that are long enough if route_length <= self.min_length: - return rt + return rt_new start = np.random.randint(0, route_length - self.min_length) length = np.random.randint(self.min_length, min(self.max_length, route_length - start)) end = start + length n_points = length - rt = np.concatenate([rt[:start], self.bezier_curve(rt[start:end], n_points), rt[end:]], axis=0) + rt_new = np.concatenate([rt_new[:start], self.bezier_curve(rt_new[start:end], n_points), rt_new[end:]], axis=0) - return rt + return rt_new class RandomWalkMutation(MutationConstraintRejection): diff --git a/WeatherRoutingTool/algorithms/genetic/patcher.py b/WeatherRoutingTool/algorithms/genetic/patcher.py index 48dabb21..8b689522 100644 --- a/WeatherRoutingTool/algorithms/genetic/patcher.py +++ b/WeatherRoutingTool/algorithms/genetic/patcher.py @@ -182,7 +182,8 @@ def _setup_configuration(self): "ROUTE_PATH", "BOAT_UNDER_KEEL_CLEARANCE", "BOAT_DRAUGHT_AFT", - "BOAT_DRAUGHT_FORE" + "BOAT_DRAUGHT_FORE", + "GENETIC_POPULATION_SIZE" ], ) cfg_path = Path(os.path.dirname(__file__)) / "configs" / "config.isofuel_single_route.json" @@ -208,7 +209,8 @@ def _setup_configuration(self): # set config path to patcher configuration cfg.CONFIG_PATH = cfg_path self.config = cfg - print('self.config: ', cfg) + self.config.ISOCHRONE_NUMBER_OF_ROUTES = cfg_select["GENETIC_POPULATION_SIZE"] + return def _setup_components(self) -> tuple[WeatherCond, Boat, WaterDepth, ConstraintsList]: """ diff --git a/WeatherRoutingTool/algorithms/genetic/population.py b/WeatherRoutingTool/algorithms/genetic/population.py index eb5b44a6..24ee40b0 100644 --- a/WeatherRoutingTool/algorithms/genetic/population.py +++ b/WeatherRoutingTool/algorithms/genetic/population.py @@ -269,6 +269,7 @@ def generate(self, problem, n_samples, **kw): Access i'th route as ``X[i,0]`` and the j'th coordinate pair off the i'th route as ``X[i,0][j, :]``. :rtype: np.array """ + boat_speed = self.boat_speed if self.boat_speed_from_arrival_time: boat_speed = 6 * u.meter / u.second # add dummy speed, will be recalculated @@ -297,6 +298,7 @@ def __init__(self, config: Config, default_route, constraints_list, pop_size): constraints_list=constraints_list, pop_size=pop_size ) + self.algo = GcrSliderAlgorithm(config) def generate(self, problem, n_samples, **kw): diff --git a/WeatherRoutingTool/algorithms/genetic/problem.py b/WeatherRoutingTool/algorithms/genetic/problem.py index 0fb98ddc..adcbbc5a 100644 --- a/WeatherRoutingTool/algorithms/genetic/problem.py +++ b/WeatherRoutingTool/algorithms/genetic/problem.py @@ -1,3 +1,4 @@ +import datetime import logging import numpy as np @@ -6,48 +7,122 @@ from WeatherRoutingTool.algorithms.data_utils import get_speed_from_arrival_time from WeatherRoutingTool.algorithms.genetic.utils import get_constraints +from WeatherRoutingTool.constraints.constraints import ConstraintsList from WeatherRoutingTool.routeparams import RouteParams +from WeatherRoutingTool.ship.ship import Boat +import WeatherRoutingTool.algorithms.genetic.utils as utils logger = logging.getLogger('WRT.Genetic') class RoutingProblem(ElementwiseProblem): - """GA definition of the Weather Routing Problem""" + """ + Definition of the Weather Routing Problem. - def __init__(self, departure_time, arrival_time, boat, boat_speed, constraint_list): - super().__init__(n_var=1, n_obj=1, n_constr=1) + This class defines the optimization problem for finding the best weather-dependent + route using the pymoo framework. It handles the evaluation of fuel consumption, + arrival time accuracy, and navigational constraints. + + :param departure_time: The time of departure. + :type departure_time: datetime.datetime + :param arrival_time: The desired time of arrival. + :type arrival_time: datetime.datetime + :param boat: Boat object for calculating fuel and power consumption. + :type boat: Boat + :param boat_speed: Boat speed. Only used to set self.boat_speed_from_arrival_time. + :type boat_speed: float + :param constraint_list: List of constraints to be checked. + :type constraint_list: ConstraintsList + :param objectives: dictionary of objective names and respective user weights. + :type objectives: dict + + """ + + def __init__(self, + departure_time: datetime.datetime, + arrival_time: datetime.datetime, + boat: Boat, + boat_speed: float, + constraint_list: ConstraintsList, + objectives: dict + ): + super().__init__( + n_var=1, + n_obj=len(objectives), + n_constr=1 + + ) self.boat = boat self.constraint_list = constraint_list self.departure_time = departure_time self.arrival_time = arrival_time - self.boat_speed = boat_speed self.boat_speed_from_arrival_time = False - if boat_speed.value == -99.: + self.objectives = objectives + + if boat_speed is None: self.boat_speed_from_arrival_time = True - def _evaluate(self, x, out, *args, **kwargs): - """Overridden function for population evaluation + def get_objectives(self, obj_dict: dict) -> np.array: + """ + Convert dictionary of objective values into a numpy array for pymoo. + + :param obj_dict: Dictionary containing calculated metrics like 'time_obj' and 'fuel_sum'. + :type obj_dict: dict + :return: A column-stacked array of objective values. + :rtype: np.ndarray + :raises ValueError: If no valid objectives are specified or found in the dictionary. + """ + objective_keys = list(self.objectives.keys()) + objs = None + if "arrival_time" in objective_keys: + objs = np.column_stack([obj_dict["time_obj"]]) + if "fuel_consumption" in objective_keys: + if objs is None: + objs = np.column_stack([obj_dict["fuel_sum"].value]) + else: + objs = np.column_stack((objs, [obj_dict["fuel_sum"].value])) + + if objs is None: + raise ValueError('Please specify an objective for the genetic algorithm.') + + return objs + + def _evaluate(self, x: np.ndarray, out: dict, *args, **kwargs) -> None: + """Overridden function for population evaluation. :param x: numpy matrix with shape (rows: number of solutions/individuals, columns: number of design variables) - :type x: numpy matrix + :type x: np.ndarray :param out: out['F']: function values, vector of length of number of solutions out['G']: constraints :type out: dict - :param *args: - :param **kwargs: """ # logger.debug(f"RoutingProblem._evaluate: type(x)={type(x)}, x.shape={x.shape}, x={x}") - fuel, _ = self.get_power(x[0]) - constraints = get_constraints(x[0], self.constraint_list) - # print(costs.shape) - out['F'] = np.column_stack([fuel]) + obj_dict = self.get_power(x[0]) + constraints = utils.get_constraints(x[0], self.constraint_list) + out['F'] = self.get_objectives(obj_dict) out['G'] = np.column_stack([constraints]) - def get_power(self, route): + def get_power(self, route: np.array) -> dict: + """ + Calculate objective values for fuel consumption and arrival-time accuracy for a specific route. + + This method extracts speed data, calculates weather-dependent ship parameters, and determines the deviation + from the target arrival time. + + :param route: A 2D numpy array where columns represent [latitude, longitude, speed]. + :type route: np.ndarray + :return: A dictionary containing the total fuel consumption ('fuel_sum'), further ship parameters + ('shipparams'), and the objective value for the arrival-time accuracy. + :rtype: dict + """ + debug = False + bs = route[:, 2] - bs = bs[:-1] * u.meter/u.second + bs = bs[:-1] * u.meter / u.second + fuel_obj = None + time_obj = None if self.boat_speed_from_arrival_time: bs = get_speed_from_arrival_time( @@ -56,6 +131,7 @@ def get_power(self, route): departure_time=self.departure_time, arrival_time=self.arrival_time, ) + bs = np.full(route[:, 1].shape[0] - 1, bs) * u.meter / u.second route_dict = RouteParams.get_per_waypoint_coords( route[:, 1], @@ -71,6 +147,35 @@ def get_power(self, route): speed=bs, ) - fuel = shipparams.get_fuel_rate() - fuel = fuel * route_dict['travel_times'] - return np.sum(fuel), shipparams + if "fuel_consumption" in self.objectives.keys(): + fuel = shipparams.get_fuel_rate() + fuel = fuel * route_dict['travel_times'] + fuel_spread = np.max(fuel) - np.min(fuel) + if debug: + print('max fuel: ', np.max(fuel)) + print('min fuel: ', np.min(fuel)) + print('fuel max spread: ', fuel_spread) + + print('last start_time: ', route_dict['start_times'][-1]) + print('last travel time: ', route_dict['travel_times'][-1].value) + fuel_obj = np.sum(fuel) + + if "arrival_time" in self.objectives.keys(): + real_arrival_time = route_dict['start_times'][-1] + datetime.timedelta( + seconds=route_dict['travel_times'][-1].value) + time_diff = np.abs(self.arrival_time - real_arrival_time).total_seconds() / 60 + + # set minimal time difference to 1 minute + if time_diff < 1: + time_diff = 1 + + time_obj = time_diff * time_diff * time_diff * time_diff + if debug: + print('departure time: ', self.departure_time) + print('planned arrival time:', self.arrival_time) + + print('real arrival time: ', real_arrival_time) + print('time_diff: ', time_diff) + print('time obj.: ', time_obj) + + return {"fuel_sum": fuel_obj, "shipparams": shipparams, "time_obj": time_obj} diff --git a/WeatherRoutingTool/algorithms/routingalg.py b/WeatherRoutingTool/algorithms/routingalg.py index fc726543..777b198f 100644 --- a/WeatherRoutingTool/algorithms/routingalg.py +++ b/WeatherRoutingTool/algorithms/routingalg.py @@ -52,7 +52,7 @@ def __init__(self, config): self.boat_speed = config.BOAT_SPEED if self.boat_speed is not None: - self.boat_speed = self.boat_speed * u.meter / u.second + self.boat_speed = self.boat_speed * u.meter/u.second def get_boat_speed(self): return self.boat_speed diff --git a/WeatherRoutingTool/config.py b/WeatherRoutingTool/config.py index 79b670aa..d75b9fb0 100644 --- a/WeatherRoutingTool/config.py +++ b/WeatherRoutingTool/config.py @@ -4,7 +4,7 @@ import sys from datetime import datetime, timedelta from pathlib import Path -from typing import Annotated, List, Literal, Optional, Self, Union +from typing import Annotated, List, Literal, Optional, Self, Union, Dict import pandas as pd import xarray as xr @@ -56,11 +56,12 @@ class Config(BaseModel): ALGORITHM_TYPE: Literal[ 'dijkstra', 'gcr_slider', 'genetic', 'genetic_shortest_route', 'isofuel', 'speedy_isobased' ] = 'isofuel' - ARRIVAL_TIME: datetime = None # arrival time at destination, format: 'yyyy-mm-ddThh:mmZ' + ARRIVAL_TIME: datetime | None = None # arrival time at destination, format: 'yyyy-mm-ddThh:mmZ' BOAT_TYPE: Literal['CBT', 'SAL', 'speedy_isobased', 'direct_power_method'] = 'direct_power_method' - BOAT_SPEED: float = None # boat speed [m/s] - BOAT_SPEED_MAX: float = 10 # maximum possible boat speed [m/s] + BOAT_SPEED: float | None = None # boat speed [m/s] + BOAT_SPEED_BOUNDARIES: Annotated[list[Union[float, float]], Field(min_length=2, max_length=2)] = [1., 10.] + # minimum and maximum possible boat speed [m/s] CONSTRAINTS_LIST: List[Literal[ 'land_crossing_global_land_mask', 'land_crossing_polygons', 'seamarks', 'water_depth', 'on_map', 'via_waypoints', 'status_error' @@ -106,11 +107,14 @@ class Config(BaseModel): 'waypoints_infill', 'constraint_violation', 'no_repair' ]] = ["waypoints_infill", "constraint_violation"] GENETIC_MUTATION_TYPE: Literal[ - 'random', 'rndm_walk', 'rndm_plateau', 'route_blend', 'percentage_change_speed', 'gaussian_speed', 'no_mutation' + 'random', 'speed', 'waypoints', 'rndm_walk', 'rndm_plateau', 'route_blend', + 'percentage_change_speed', 'gaussian_speed', 'no_mutation' ] = 'random' - GENETIC_CROSSOVER_TYPE: Literal['random', 'speed'] = 'random' + GENETIC_CROSSOVER_TYPE: Literal['random', 'speed', "waypoints"] = 'random' GENETIC_CROSSOVER_PATCHER: Literal['gcr', 'isofuel'] = 'isofuel' GENETIC_FIX_RANDOM_SEED: bool = False + GENETIC_OBJECTIVES: Dict[str, float] = {"arrival_time": 1.5, "fuel_consumption": 1.5} # dictionary for configuring + # the objectives and objective weights INTERMEDIATE_WAYPOINTS: Annotated[ list[Annotated[list[Union[int, float]], Field(min_length=2, max_length=2)]], @@ -493,11 +497,37 @@ def check_boat_speed(cls, v): def check_speed_determination(self) -> Self: logger.info(f'arrival time: {self.ARRIVAL_TIME}') logger.info(f'speed: {self.BOAT_SPEED}') - if self.ARRIVAL_TIME is None and self.BOAT_SPEED is None: - raise ValueError('Please specify EITHER the boat speed OR the arrival time.') - if self.ARRIVAL_TIME is not None and self.BOAT_SPEED is not None: - raise ValueError('Please specify EITHER the boat speed OR the arrival time but not both.') - if self.ARRIVAL_TIME is not None and self.ALGORITHM_TYPE not in ['genetic', 'gcr_slider']: - raise ValueError('The determination of the speed from the arrival time is only possible for the' - ' genetic and the gcr slider algorithm') + + mutate_only_waypoints = ((self.GENETIC_MUTATION_TYPE == "waypoints") + or (self.GENETIC_MUTATION_TYPE == "rndm_walk") + or (self.GENETIC_MUTATION_TYPE == "rndm_plateau") + or (self.GENETIC_MUTATION_TYPE == "route_blend")) + + crossover_only_waypoints = self.GENETIC_CROSSOVER_TYPE == "waypoints" + + if self.ALGORITHM_TYPE == "genetic": + # run mode: route optimisation with constant speed + if mutate_only_waypoints and crossover_only_waypoints: + logger.info('Algorithm run mode: only waypoint optimisation.') + if self.ARRIVAL_TIME is None and self.BOAT_SPEED is None: + raise ValueError('Please specify EITHER the boat speed OR the arrival time.') + if self.ARRIVAL_TIME is not None and self.BOAT_SPEED is not None: + raise ValueError('Please specify EITHER the boat speed OR the arrival time but not both.') + + # run modes: speed optimisation for fixed route as well as simultaneous waypoint and speed optimisation + else: + logger.info('Algorithm run mode: speed optimisation for fixed route or simultaneous ' + 'waypoint and speed optimisation.') + if self.ARRIVAL_TIME is None or self.BOAT_SPEED is None: + raise ValueError('Please provide a valid arrival time and boat speed.') + + if self.GENETIC_MUTATION_TYPE == "speed" and self.GENETIC_CROSSOVER_TYPE == "speed": + raise NotImplementedError("Pure speed optimisation of single routes is not yet implemented but planned " + "for the future.") + else: + if self.BOAT_SPEED is None: + raise ValueError('Please provide a valid boat speed.') + if self.ARRIVAL_TIME: + logger.warning('You specified an arrival time. The arrival time is only used as input for the ' + 'optimisation process for the genetic algorithm.') return self diff --git a/WeatherRoutingTool/execute_routing.py b/WeatherRoutingTool/execute_routing.py index 9958b529..ed185ab7 100644 --- a/WeatherRoutingTool/execute_routing.py +++ b/WeatherRoutingTool/execute_routing.py @@ -66,7 +66,7 @@ def execute_routing(config, ship_config): # routing min_fuel_route, error_code = alg.execute_routing(boat, wt, constraint_list) # min_fuel_route.print_route() - min_fuel_route.write_to_geojson(routepath / f"{min_fuel_route.route_type}.geojson") + min_fuel_route.write_to_geojson(str(routepath) + '/' + str(min_fuel_route.route_type) + ".json") if config.ROUTE_POSTPROCESSING: postprocessed_route = RoutePostprocessing(min_fuel_route, boat) diff --git a/WeatherRoutingTool/routeparams.py b/WeatherRoutingTool/routeparams.py index 064db28d..02cbb88a 100644 --- a/WeatherRoutingTool/routeparams.py +++ b/WeatherRoutingTool/routeparams.py @@ -417,6 +417,26 @@ def plot_power_vs_dist(self, color, label, power_type, ax, bin_center_mean=None, plt.xlabel('travel distance (km)') plt.xticks() + def plot_speed_vs_dist(self, color, label, ax): + speed = self.ship_params_per_step.get_speed() + dist = self.dists_per_step + + hist_values = graphics.get_hist_values_from_widths(dist, speed, "speed") + + # only for power: also plot bin showing weighted mean. This does not make sense for fuel. + plt.ylabel("speed (m/s)") + plt.bar( + hist_values["bin_centres"].to(u.km).value, + hist_values["bin_contents"].to(u.m/u.second).value, + hist_values["bin_widths"].to(u.km).value, + alpha=0.5, color=color, edgecolor=color, label=label, linewidth=2) + + left, right = plt.xlim() + ax.set_xlim(-100, right) + + plt.xlabel('travel distance (km)') + plt.xticks() + # TODO check whether correct: Why do we see steps and no smooth curve? def plot_acc_power_vs_dist(self, color, label, power_type): power = self.get_power_type(power_type) diff --git a/WeatherRoutingTool/utils/graphics.py b/WeatherRoutingTool/utils/graphics.py index a446b3e2..6d54c58c 100644 --- a/WeatherRoutingTool/utils/graphics.py +++ b/WeatherRoutingTool/utils/graphics.py @@ -10,6 +10,8 @@ import pandas as pd from astropy import units as u from geovectorslib import geod +from matplotlib.collections import LineCollection +from matplotlib.colors import Normalize from matplotlib.figure import Figure from PIL import Image @@ -230,6 +232,8 @@ def get_hist_values_from_widths(bin_widths, contend_unnormalised, power_type): contents = np.array([]) if power_type == 'fuel': contents = contents * u.kg / u.meter + elif power_type == 'speed': + contents = contents * u.meter / u.second else: contents = contents * u.Watt cent_temp = 0 * u.meter @@ -373,3 +377,18 @@ def write_graph_to_csv(path, x, y): for i in range(0, len(x)): graphwriter.writerow([x[i], y[i]]) + + +def get_route_lc(X): + lats = X[:, 0] + lons = X[:, 1] + speed = X[:, 2] + + points = np.array([lons, lats]).T.reshape(-1, 1, 2) + segments = np.concatenate([points[:-1], points[1:]], axis=1) + + norm = Normalize(vmin=0, vmax=20) + lc = LineCollection(segments, cmap='viridis', norm=norm, transform=ccrs.Geodetic()) + lc.set_array(speed) + lc.set_linewidth(1) + return lc diff --git a/docs/source/algorithms/genetic.rst b/docs/source/algorithms/genetic.rst index 8932adf0..b4c994dc 100644 --- a/docs/source/algorithms/genetic.rst +++ b/docs/source/algorithms/genetic.rst @@ -3,12 +3,53 @@ Genetic Algorithm ================= +The Weather Routing Tool makes use of `pymoo `__ +as the supporting library for the Genetic algorithm’s implementation. The Genetic algorithm + - considers weather conditions, + - considers constraints, + - can be used for waypoint optimisation as well as combined waypoint and boat speed optimisation (degrees of freedom: waypoints and/or speed) and + - can be used for optimisation of fuel consumption and arrival-time accuracy (objectives: fuel consumption and/or arrival-time accuracy). +Adding functionality for sole speed optimisation for a fixed route is planned for the near future. -Foundation ----------- -The Weather Routing Tool makes use of `pymoo `__ -as the supporting library for the Genetic algorithm’s implementation. +Different Run Modes +----------------------- +*Degrees of Freedom (DOF)* + +The DOF can be specified by setting the config variables ``GENETIC_MUTATION_TYPE`` and ``GENETIC_CROSSOVER_TYPE``. + +- pure waypoint optimisation: In case pure waypoint optimisation is requested, both config variables need to be set to ``"waypoints"``. + ``GENETIC_MUTATION_TYPE`` can also be ``"rndm_walk"``, ``"rndm_plateau"`` or ``"route_blend"``. The boat speed is taken from the + user input to ``BOAT_SPEED`` and is left constant. + +- pure speed optimisation (NOT YET IMPLEMENTED!): In case pure speed optimisation is requested, + both config variables need to be set to ``"speed"``. ``GENETIC_MUTATION_TYPE`` can also be ``"percentage_change_speed"`` or ``"gaussian_speed"``. + The boat speed of the initial population is read from the user input to ``BOAT_SPEED``. The waypoints of the route to be + optimised are read from a GeoJSON file. Only speed optimisation of a single route is allowed, meaning only one GeoJSON file can + be provided as initial population. + +- waypoint and speed optimisation: Any other combination of both config variables results in mixed speed and waypoint optimisation. The initial + population differs in waypoints but is generated with constant speed from the user input to ``BOAT_SPEED``. All generation methods for the initial + population are allowed. + +*Objectives* + +The objectives can be specified by setting the config variable ``GENETIC_OBJECTIVES``. Currently only the optimisation of the total fuel +consumption (``"fuel_consumption"``) and/or arrival-time accuracy (``"arrival_time"``) is possible. In case fuel consumption shall be optimised, the algorithm minimises the total amount of fuel that is consumed for +a route. In case the arrival-time accuracy shall be optimised, the algorithm minimises the following function of the real arrival time (t_real) +and the planned arrival time (t_planned): + +:math:`(t_{planned} - t_{real})^4.` + +Along with the objective keys, integer weights are to be specified that rank the objectives according to their importance. +E.g. ``GENETIC_OBJECTIVES={"fuel_consumption": 2, "arrival_time": 1}`` refers to optimisation of fuel consumption and arrival-time +accuracy with an emphasis on fuel-consumption optimisation. In case both objectives +are to be considered of equal importance, the mean values of the maximum possible rank shall be provided e.g. +``GENETIC_OBJECTIVES={"fuel_consumption": 1.5, "arrival_time": 1.5}``. + +The selection of the final solution is done using methods for Multi-Criteria Decision-Making (MCDM). Further details are given in the respective +section below. + General Concept --------------- @@ -320,43 +361,15 @@ Route Patching process. -**Implementation Notes:** - -The intuition behind having Route Patching implementations setup as -classes follows the following: - a. Route patching can be quite expensive during both the preparation - (defining map, loading configs, etc.) and the execution stage (patching - between point A and point B). An Object Oriented implementation of the same - helps separate the two processes, avoids redundancy and can contribute to the - overall speed in the longer run. - - b. Implementation consistency makes it easier to swap between different - Patching implementations and maintains clean code - -Config Parameters ------------------ - -1. ``GENETIC_NUMBER_GENERATIONS`` — Max number of generations. - -2. ``GENETIC_NUMBER_OFFSPRINGS`` — Number of offsprings. - -3. ``GENETIC_POPULATION_SIZE`` — Population size of the genetic algorithm. - -4. ``GENETIC_POPULATION_TYPE`` — Population generation method for the. - genetic algorithm - - a. ``GENETIC_POPULATION_PATH`` — Path to population directory when. - ``GENETIC_POPULATION_TYPE`` is “\ *from_geojson*\ ” - -5. ``GENETIC_REPAIR_TYPE`` - Repair strategy for genetic algorithm. - -6. ``GENETIC_MUTATION_TYPE`` - Mutation strategy. - -7. ``GENETIC_CROSSOVER_PATCHER`` - Patching strategy for crossover. +Multi-Criteria Decision-Making +------------------------------ +Currently, the *R-Method* by R.V. Rao and R.J. Lakshmi is implemented for MCDM. This method starts from the ranking of the objectives by the user via the config +variable ``GENETIC_OBJECTIVES``. These ranks are converted into weights :math:`w^\text{o}_{obj}` with *obj* = time, fuel. After the optimisation with the Genetic algorithm, the alternative solutions of the set of non-dominated solutions are ranked separately according to their performance +with respect to each objective. As for the objective ranks, the solution ranks are converted into weights :math:`w^s_{obj}` for solution *s* with *obj* = time, fuel. From the objective weights and the solution weights, a composite weight is determined for each solution. +In contrast to the original method, the WRT uses the following function to determine the composite weight as it was found to perform better for MCDM with two objectives -8. ``GENETIC_FIX_RANDOM_SEED`` - Handling of random seed. +:math:`c_s = \frac{w^s_\text{fuel} w^\text{o}_\text{fuel} + w^s_\text{time} w^\text{o}_\text{time}}{1/w^\text{o}_\text{fuel} w^s_\text{fuel} - 1/w^\text{o}_\text{time} w^s_\text{time} + 0.2}.` -For more details on the configuration variables, please check the general section on Configuration. Useful References ----------------- diff --git a/docs/source/algorithms/overview.rst b/docs/source/algorithms/overview.rst index b429bd12..3bf10b1e 100644 --- a/docs/source/algorithms/overview.rst +++ b/docs/source/algorithms/overview.rst @@ -3,14 +3,20 @@ Algorithm overview ================== -+------------------------+-----------------------+-------------------------+---------+-------------+ -| Algorithm / Feature | Waypoint optimization | Ship speed optimization | Weather | Constraints | -+========================+=======================+=========================+=========+=============+ -| Genetic | Yes | Yes | Yes | Yes | -+------------------------+-----------------------+-------------------------+---------+-------------+ -| Isofuel | Yes | No | Yes | Yes | -+------------------------+-----------------------+-------------------------+---------+-------------+ -| GCR Slider | Yes | No | No | Partially | -+------------------------+-----------------------+-------------------------+---------+-------------+ -| Dijkstra | Yes | No | No | Partially | -+------------------------+-----------------------+-------------------------+---------+-------------+ \ No newline at end of file +The implementations of the individual algorithms meet different levels of sophistication. In particular, the algorithms +differ in the possible *degrees of freedom (DOF)* that can be manipulated (e.g. boat speed, waypoints) and the *objectives* that +can be optimised (fuel, arrival-time accuracy, distance). The table below summarises the possible run modes for all available algorithms. +Information on the configurations for the different run modes can be found in the sections that describe the functionality of each algorithm in detail. + ++------------------------+-------------------------------------------------+-----------------------------------------+---------+-------------+ +| Algorithm / Feature | Degrees of Freedom | Objectives | Weather | Constraints | ++========================+=======================+=========================+=========================================+=========+=============+ +| Genetic | speed, waypoints, speed & waypoints | fuel consumption, arrival-time accuracy | Yes | Yes | ++------------------------+-------------------------------------------------+-----------------------------------------+---------+-------------+ +| Isofuel | waypoints | fuel consumption | Yes | Yes | ++------------------------+-------------------------------------------------+-----------------------------------------+---------+-------------+ +| GCR Slider | waypoints | distance | No | Partially | ++------------------------+-------------------------------------------------+-----------------------------------------+---------+-------------+ +| Dijkstra | waypoints | distance | No | Partially | ++------------------------+-------------------------------------------------+-----------------------------------------+---------+-------------+ + diff --git a/docs/source/api-documentation/WeatherRoutingTool.algorithms.rst b/docs/source/api-documentation/WeatherRoutingTool.algorithms.rst index 99f918d0..b7e387d7 100644 --- a/docs/source/api-documentation/WeatherRoutingTool.algorithms.rst +++ b/docs/source/api-documentation/WeatherRoutingTool.algorithms.rst @@ -60,6 +60,14 @@ WeatherRoutingTool.algorithms.genetic.repair module :undoc-members: :show-inheritance: +WeatherRoutingTool.algorithms.genetic.mcdm module +--------------------------------------------------- + +.. automodule:: WeatherRoutingTool.algorithms.genetic.mcdm + :members: + :undoc-members: + :show-inheritance: + WeatherRoutingTool.algorithms.genetic.utils module --------------------------------------------------- diff --git a/docs/source/configuration.rst b/docs/source/configuration.rst index 5132156d..21dfa89b 100644 --- a/docs/source/configuration.rst +++ b/docs/source/configuration.rst @@ -175,8 +175,6 @@ The following lists contain information on each variable which can be set. The c +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ | FACTOR_WIND_FORCES | multiplication factor for the added resistance in wind model | +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| GENETIC_MUTATION_TYPE | type for mutation (options: 'grid_based') | -+----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ | GENETIC_NUMBER_GENERATIONS | number of generations for genetic algorithm (default: 20) | +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ | GENETIC_NUMBER_OFFSPRINGS | number of offsprings for genetic algorithm (default: 2) | @@ -185,17 +183,19 @@ The following lists contain information on each variable which can be set. The c +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ | GENETIC_POPULATION_TYPE | type for initial population (options: 'grid_based', 'from_geojson', 'isofuel', 'gcrslider'; default: 'isofuel') | +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| GENETIC_POPULATION_PATH | path to initial population | +| GENETIC_POPULATION_PATH | path to initial population for input via geojson | +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ | GENETIC_REPAIR_TYPE | repair strategy for genetic algorithm (options: 'waypoints_infill', 'constraint_violation', 'no_repair', default: 'waypoints_infill' and 'constraint_violation') | +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| GENETIC_MUTATION_TYPE | options: 'random', 'rndm_walk', 'rndm_plateau', 'route_blend', 'no_mutation' (default: 'random') | +| GENETIC_MUTATION_TYPE | options: 'random', 'speed', 'waypoints', 'rndm_walk', 'rndm_plateau', 'route_blend', 'percentage_change_speed', 'gaussian_speed', 'no_mutation' (default: 'random') | +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| GENETIC_CROSSOVER_TYPE | options: 'random', 'speed' (default: 'random') | +| GENETIC_CROSSOVER_TYPE | options: 'random', 'speed', 'waypoints', (default: 'random') | +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ | GENETIC_CROSSOVER_PATCHER | patching strategy for crossover (options: 'gcr', 'isofuel', default: 'isofuel') | +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ -| GENETIC_FIX_RANDOM_SEED | options: True, False (default: False) | +| GENETIC_OBJECTIVES | dictionary of the objectives and objective weights; possible keys: "arrival_time", "fuel_consumption" (default: {"arrival_time": 1.5, "fuel_consumption": 1.5}) | ++----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ +| GENETIC_FIX_RANDOM_SEED | configuration of numpy random seed (options: True, False, default: False) | +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ | INTERMEDIATE_WAYPOINTS | coordinates for intermediate waypoints [[lat_one,lon_one], [lat_two,lon_two] ... ] (default: []) | +----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/scripts/compare_objectives.py b/scripts/compare_objectives.py new file mode 100644 index 00000000..461b01ed --- /dev/null +++ b/scripts/compare_objectives.py @@ -0,0 +1,121 @@ +from datetime import datetime, timedelta + +import argparse +import matplotlib.pyplot as plt +import numpy as np +import matplotlib.ticker as ticker + + +def time_formatter(x, pos): + """Converts total minutes to HH:MM format.""" + hours = int(x) // 60 + minutes = int(x) % 60 + return f'{hours:02d}:{minutes:02d}' + + +def time_str_convertert(time_str, days): + t = datetime.strptime(time_str, '%H:%M:%S') + delta = timedelta(days=days, hours=t.hour, minutes=t.minute, seconds=t.second) + t_min = delta.total_seconds() * 1. / 60 + print('t: ', delta) + print('t_min: ', t_min) + return t_min + + +def compare_time_obj(figure_dir: str): + # Sample data + labels = ['fuel opt.', 'time opt.', 'fuel:time opt. 1:1'] + values = [ + time_str_convertert(time_str="8:36:36", days=1), + time_str_convertert(time_str="10:12:00", days=1), + time_str_convertert(time_str="8:58:25", days=1), + ] # Time in minutes + + departure_time = datetime.strptime("2025-12-05T13:48Z", '%Y-%m-%dT%H:%MZ') + arrival_time = datetime.strptime("2025-12-07T00:00Z", '%Y-%m-%dT%H:%MZ') + optimal_travel_time = (arrival_time - departure_time).total_seconds() * 1. / 60 + + arrival_time = 140 # Arrival time in minutes + + # Create the plot + fig, ax = plt.subplots(figsize=(10, 6)) + + # Generate the bar plot + x_pos = np.arange(len(labels)) + ax.bar(x_pos, values, color='skyblue', edgecolor='navy', alpha=0.7, label='Recorded Time') + + # Set manual labels for x-axis + ax.set_xticks(x_pos) + ax.set_xticklabels(labels) + + # Add the constant 'arrival time' line + ax.axhline(y=optimal_travel_time, color='red', linestyle='--', linewidth=2, label='Arrival Time') + + # Add the shaded area (+/- 30 minutes) + ax.fill_between([-0.5, len(labels) - 0.5], + optimal_travel_time - 30, + optimal_travel_time + 30, + color='gray', alpha=0.2, label=r'Arrival Window ($\pm 30$ min)') + + # Apply the custom formatter + ax.yaxis.set_major_formatter(ticker.FuncFormatter(time_formatter)) + + # Set tick frequency to every 30 minutes for clarity + ax.yaxis.set_major_locator(ticker.MultipleLocator(30)) + ax.set_ylim(25 * 60, 35 * 60) + + # General formatting + ax.set_ylabel('Time (HH:MM)') + ax.set_xlabel('Run ID') + ax.legend(loc='upper left') + ax.grid(axis='y', linestyle='--', alpha=0.5) + + figure_path = f'{figure_dir}/arrival_time.png' + plt.tight_layout() + plt.savefig(figure_path) + + +def compare_fuel_obj(figure_dir: str): + # Sample data + labels = ['fuel opt.', 'time opt.', 'fuel:time opt. 1:1'] + values = [ + 26567.46247030554, + 27755.58059942486, + 26742.834078860997, + ] # Time in minutes + + # Create the plot + fig, ax = plt.subplots(figsize=(10, 6)) + + # Generate the bar plot + x_pos = np.arange(len(labels)) + ax.bar(x_pos, values, color='skyblue', edgecolor='navy', alpha=0.7, label='fuel consumption') + + # Set manual labels for x-axis + ax.set_xticks(x_pos) + ax.set_xticklabels(labels) + ax.set_ylim(20000, 30000) + + # General formatting + ax.set_ylabel('fuel consumption (kg)') + ax.set_xlabel('Run ID') + ax.legend(loc='upper left') + ax.grid(axis='y', linestyle='--', alpha=0.5) + + figure_path = f'{figure_dir}/fuel_consumption.png' + plt.tight_layout() + plt.savefig(figure_path) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description='Weather Routing Tool') + required_args = parser.add_argument_group('required arguments') + required_args.add_argument('--figure-dir', help="Figure directory (absolute path).", + required=True, type=str) + + args = parser.parse_args() + figure_dir = args.figure_dir + + # Compare variations of resistances for specific routes + compare_time_obj(figure_dir) + compare_fuel_obj(figure_dir) diff --git a/scripts/compare_routes.py b/scripts/compare_routes.py index 47120c87..9cb38e7c 100644 --- a/scripts/compare_routes.py +++ b/scripts/compare_routes.py @@ -27,6 +27,20 @@ def plot_power_vs_dist(rp_list, rp_str_list, scenario_str, power_type='fuel'): plt.savefig(figurefile + '/' + power_type + '_vs_dist.png') +def plot_speed_vs_dist(rp_list, rp_str_list, scenario_str): + fig, ax = plt.subplots(figsize=(12, 8), dpi=96) + # ax.set_ylim(4000, 5500) + for irp in range(0, len(rp_list)): + rp_list[irp].plot_speed_vs_dist(graphics.get_colour(irp), rp_str_list[irp], ax) + + ax.legend(loc='upper left', frameon=False) + ax.tick_params(top=True, right=True) + # ax.tick_params(labelleft=False, left=False, top=True) # hide y labels + ax.text(0.95, 0.96, scenario_str, verticalalignment='top', horizontalalignment='right', + transform=ax.transAxes) + plt.savefig(figurefile + '/speed_vs_dist.png') + + def plot_acc_power_vs_dist(rp_list, rp_str_list, power_type='fuel'): fig, ax = plt.subplots(figsize=(12, 8), dpi=96) for irp in range(0, len(rp_list)): @@ -89,7 +103,8 @@ def plot_power_vs_dist_ratios(rp_list, rp_str_list, scenario_str, power_type='fu 'fuel_vs_lat': False, 'power_vs_dist_showing_weather': False, 'power_vs_dist_ratios': False, - 'fuel_vs_dist_ratios': False + 'fuel_vs_dist_ratios': False, + 'speed_vs_dist': False } parser = argparse.ArgumentParser(description='Weather Routing Tool') @@ -188,7 +203,7 @@ def plot_power_vs_dist_ratios(rp_list, rp_str_list, scenario_str, power_type='fu # plotting routes in depth profile if hist_dict['route']: fig, ax = plt.subplots(figsize=graphics.get_standard('fig_size')) - depth = xr.open_dataset(depth_path) + # depth = xr.open_dataset(depth_path) ax.axis('off') ax.xaxis.set_tick_params(labelsize='large') fig, ax = graphics.generate_basemap( @@ -211,6 +226,9 @@ def plot_power_vs_dist_ratios(rp_list, rp_str_list, scenario_str, power_type='fu if hist_dict['power_vs_dist']: plot_power_vs_dist(rp_list, rp_str_list, scenario_str, 'power') + if hist_dict['speed_vs_dist']: + plot_speed_vs_dist(rp_list, rp_str_list, scenario_str) + if hist_dict['fuel_vs_dist']: plot_power_vs_dist(rp_list, rp_str_list, scenario_str, 'fuel') @@ -251,7 +269,7 @@ def plot_power_vs_dist_ratios(rp_list, rp_str_list, scenario_str, power_type='fu plot_power_vs_dist_ratios(rp_list, rp_str_list, scenario_str, 'fuel') ## - # write fuel consumption, travel distance and time + # write fuel consumption, trafvel distance and time print('Full fuel consumption:') for irp in range(0, len(rp_list)): print(rp_str_list[irp] + ': ' + str(rp_list[irp].get_full_fuel())) diff --git a/tests/test_config.py b/tests/test_config.py index 10d2541d..d680ef6a 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -131,3 +131,76 @@ def test_weather_start_time_compatibility(): Config.assign_config(init_mode="from_dict", config_dict=config_data) assert "Weather data does not cover the full routing time range." in str(excinfo.value) + + +@pytest.mark.parametrize("boat_speed,arrival_time,mut_type,cross_type,ierr", [ + (7, "2025-12-07T00:00Z", "waypoints", "waypoints", 0), + (None, None, "waypoints", "waypoints", 1), +]) +def test_boat_speed_arrival_time_waypoint_optimisation_failure(boat_speed, arrival_time, mut_type, cross_type, ierr): + config_data, _ = load_example_config() + config_data["BOAT_SPEED"] = boat_speed + config_data["ARRIVAL_TIME"] = arrival_time + config_data["GENETIC_MUTATION_TYPE"] = mut_type + config_data["GENETIC_CROSSOVER_TYPE"] = cross_type + config_data["ALGORITHM_TYPE"] = "genetic" + error_str_list = [ + "Please specify EITHER the boat speed OR the arrival time but not both.", + "Please specify EITHER the boat speed OR the arrival time.", + ] + + with pytest.raises(ValueError) as excinfo: + Config.assign_config(init_mode="from_dict", config_dict=config_data) + + assert error_str_list[ierr] in str(excinfo.value) + + +@pytest.mark.parametrize("boat_speed,arrival_time,mut_type,cross_type", [ + (7, None, "waypoints", "waypoints"), + (None, "2025-12-07T00:00Z", "waypoints", "waypoints"), +]) +def test_boat_speed_arrival_time_waypoint_optimisation_success(boat_speed, arrival_time, mut_type, cross_type, ): + config_data, _ = load_example_config() + config_data["BOAT_SPEED"] = boat_speed + config_data["ARRIVAL_TIME"] = arrival_time + config_data["GENETIC_MUTATION_TYPE"] = mut_type + config_data["GENETIC_CROSSOVER_TYPE"] = cross_type + config_data["ALGORITHM_TYPE"] = "genetic" + + Config.assign_config(init_mode="from_dict", config_dict=config_data) + + +@pytest.mark.parametrize("boat_speed,arrival_time,mut_type,cross_type", [ + (7, None, "random", "random"), + (None, "2025-12-07T00:00Z", "random", "random"), + (None, "2025-12-07T00:00Z", "waypoints", "random"), + (None, "2025-12-07T00:00Z", "random", "waypoints"), +]) +def test_boat_speed_arrival_time_speed_optimisation_failure(boat_speed, arrival_time, mut_type, cross_type, ): + config_data, _ = load_example_config() + config_data["BOAT_SPEED"] = boat_speed + config_data["ARRIVAL_TIME"] = arrival_time + config_data["GENETIC_MUTATION_TYPE"] = mut_type + config_data["GENETIC_CROSSOVER_TYPE"] = cross_type + config_data["ALGORITHM_TYPE"] = "genetic" + + with pytest.raises(ValueError) as excinfo: + Config.assign_config(init_mode="from_dict", config_dict=config_data) + + assert "Please provide a valid arrival time and boat speed." in str(excinfo.value) + + +@pytest.mark.parametrize("boat_speed,arrival_time,mut_type,cross_type", [ + (7, "2025-12-07T00:00Z", "random", "random"), + (7, "2025-12-07T00:00Z", "random", "waypoints"), + (7, "2025-12-07T00:00Z", "waypoints", "random"), +]) +def test_boat_speed_arrival_time_speed_optimisation_success(boat_speed, arrival_time, mut_type, cross_type, ): + config_data, _ = load_example_config() + config_data["BOAT_SPEED"] = boat_speed + config_data["ARRIVAL_TIME"] = arrival_time + config_data["GENETIC_MUTATION_TYPE"] = mut_type + config_data["GENETIC_CROSSOVER_TYPE"] = cross_type + config_data["ALGORITHM_TYPE"] = "genetic" + + Config.assign_config(init_mode="from_dict", config_dict=config_data) diff --git a/tests/test_genetic.py b/tests/test_genetic.py index 8a6ac01c..7042c075 100644 --- a/tests/test_genetic.py +++ b/tests/test_genetic.py @@ -8,12 +8,11 @@ import matplotlib.pyplot as pyplot import pytest from astropy import units as u -from matplotlib.collections import LineCollection -from matplotlib.colors import Normalize import tests.basic_test_func as basic_test_func import WeatherRoutingTool.utils.graphics as graphics -from WeatherRoutingTool.algorithms.genetic.crossover import SinglePointCrossover +from WeatherRoutingTool.algorithms.genetic import Genetic +from WeatherRoutingTool.algorithms.genetic.crossover import SinglePointCrossover, SpeedCrossover from WeatherRoutingTool.algorithms.genetic.patcher import PatcherBase, GreatCircleRoutePatcher, IsofuelPatcher, \ GreatCircleRoutePatcherSingleton, IsofuelPatcherSingleton, PatchFactory from WeatherRoutingTool.algorithms.genetic.population import IsoFuelPopulation @@ -99,21 +98,6 @@ def get_dummy_route_input(length='long'): ''' -def get_route_lc(X): - lats = X[:, 0] - lons = X[:, 1] - speed = X[:, 2] - - points = np.array([lons, lats]).T.reshape(-1, 1, 2) - segments = np.concatenate([points[:-1], points[1:]], axis=1) - - norm = Normalize(vmin=10, vmax=20) - lc = LineCollection(segments, cmap='viridis', norm=norm, transform=ccrs.Geodetic()) - lc.set_array(speed) - lc.set_linewidth(3) - return lc - - @pytest.mark.manual def test_random_plateau_mutation(plt): dirname = os.path.dirname(__file__) @@ -139,10 +123,10 @@ def test_random_plateau_mutation(plt): show_depth=False, show_gcr=False ) - old_route_one_lc = get_route_lc(old_route[0, 0]) - old_route_two_lc = get_route_lc(old_route[1, 0]) - new_route_one_lc = get_route_lc(new_route[0, 0]) - new_route_two_lc = get_route_lc(new_route[1, 0]) + old_route_one_lc = graphics.get_route_lc(old_route[0, 0]) + old_route_two_lc = graphics.get_route_lc(old_route[1, 0]) + new_route_one_lc = graphics.get_route_lc(new_route[0, 0]) + new_route_two_lc = graphics.get_route_lc(new_route[1, 0]) ax.add_collection(old_route_one_lc) ax.add_collection(old_route_two_lc) ax.add_collection(new_route_one_lc) @@ -214,10 +198,10 @@ def test_bezier_curve_mutation(plt): show_gcr=False ) - old_route_one_lc = get_route_lc(old_route[0, 0]) - old_route_two_lc = get_route_lc(old_route[1, 0]) - new_route_one_lc = get_route_lc(new_route[0, 0]) - new_route_two_lc = get_route_lc(new_route[1, 0]) + old_route_one_lc = graphics.get_route_lc(old_route[0, 0]) + old_route_two_lc = graphics.get_route_lc(old_route[1, 0]) + new_route_one_lc = graphics.get_route_lc(new_route[0, 0]) + new_route_two_lc = graphics.get_route_lc(new_route[1, 0]) ax.add_collection(old_route_one_lc) ax.add_collection(old_route_two_lc) ax.add_collection(new_route_one_lc) @@ -308,8 +292,8 @@ def test_constraint_violation_repair(plt): show_depth=False, show_gcr=False ) - old_route_lc = get_route_lc(old_route[0, 0]) - new_route_lc = get_route_lc(new_route) + old_route_lc = graphics.get_route_lc(old_route[0, 0]) + new_route_lc = graphics.get_route_lc(new_route) ax.add_collection(old_route_lc) ax.add_collection(new_route_lc) @@ -389,3 +373,44 @@ def test_single_point_crossover(plt): ax.plot(old_route[1, 0][:, 1], old_route[0, 0][:, 0], color="orange", transform=input_crs, marker='o') plt.saveas = "test_single_point_crossoverr.png" + + +def test_speed_crossover(plt): + dirname = os.path.dirname(__file__) + configpath = os.path.join(dirname, 'config.isofuel_single_route.json') + config = Config.assign_config(Path(configpath)) + default_map = Map(32., 15, 36, 29) + constraint_list = basic_test_func.generate_dummy_constraint_list() + departure_time = datetime(2025, 4, 1, 11, 11) + + X = get_dummy_route_input() + + sp = SpeedCrossover(config=config, departure_time=departure_time, constraints_list=constraint_list) + o1, o2 = sp.crossover(X[0, 0], X[1, 0]) + + # plot figure with original and mutated routes + fig, ax = graphics.generate_basemap( + map=default_map.get_var_tuple(), + depth=None, + start=(35.199, 15.490), + finish=(32.737, 28.859), + title='', + show_depth=False, + show_gcr=False + ) + old_X1_lc = graphics.get_route_lc(X[0, 0]) + old_X2_lc = graphics.get_route_lc(X[1, 0]) + + new_X1_lc = graphics.get_route_lc(o1) + new_X2_lc = graphics.get_route_lc(o2) + + ax.add_collection(old_X1_lc) + ax.add_collection(old_X2_lc) + ax.add_collection(new_X1_lc) + ax.add_collection(new_X2_lc) + + cbar = fig.colorbar(old_X2_lc, ax=ax, orientation='vertical', pad=0.15, shrink=0.7) + cbar.set_label('Geschwindigkeit ($m/s$)') + + pyplot.tight_layout() + plt.saveas = "test_speed_crossover.png" diff --git a/tests/test_genetic_mcdm.py b/tests/test_genetic_mcdm.py new file mode 100644 index 00000000..ed7c0f5d --- /dev/null +++ b/tests/test_genetic_mcdm.py @@ -0,0 +1,68 @@ +import numpy as np +import pytest + +from WeatherRoutingTool.algorithms.genetic.mcdm import RMethod + + +@pytest.mark.parametrize("obj_fuel,obj_time", [(1, 1), (1, 2), (2, 1)]) +def test_weight_determination_for_solution_selection(plt, obj_fuel, obj_time): + fuel_weight = np.random.rand(1, 10000) * 0.1 + time_weight = np.random.rand(1, 10000) * 0.1 + + objective_dict = { + "fuel_consumption": obj_fuel, + "arrival_time": obj_time, + } + mcdm = RMethod(objective_dict) + composite_weight = mcdm.get_composite_weight( + sol_weight_list=[time_weight, fuel_weight], + obj_weight_list=[obj_time, obj_fuel] + ) + + fig = plt.figure(figsize=(10, 7)) + ax = fig.add_subplot(111, projection='3d') + + # Plot the scatter points + ax.set_xlim(fuel_weight.max(), fuel_weight.min()) + ax.scatter(fuel_weight, time_weight, composite_weight, c=composite_weight, cmap='viridis', marker='o', s=40, + alpha=0.6, edgecolors='w') + + # Set labels and title + ax.set_title('3D Scatter Plot of Point Selections') + ax.set_xlabel('fuel weight') + ax.set_ylabel('time weight') + ax.set_zlabel('composite weight') + + plt.saveas = f"test_composite_weight_fuel{obj_fuel}_time{obj_time}.png" + + +@pytest.mark.parametrize("rank,out", [ + (1, 1.), + (2, 0.666666), + (3, 0.545454), + (4, 0.48), +]) +def test_get_rank_sum(rank, out): + objective_dict = { + "fuel_consumption": 1, + "arrival_time": 1, + } + mcdm = RMethod(objective_dict) + res = mcdm.get_rank_sum(rank) + assert np.isclose(res, out) + + +@pytest.mark.parametrize("rank,n_parts,out", [ + (4, 4, 0.48 / 2.69212), + (50, 50, 0.2222614 / 15.287014), + (1.5, 4, 0.309545), +]) +def test_get_weigth_from_rank(rank, out, n_parts): + objective_dict = { + "fuel_consumption": 1, + "arrival_time": 1, + } + rank_arr = np.array([rank]) + mcdm = RMethod(objective_dict) + res = mcdm.get_weigths_from_rankarr(rank_arr, n_parts) + assert np.isclose(res, out) diff --git a/tests/test_utils.py b/tests/test_utils.py index 236c8b26..dacf5e54 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -5,6 +5,7 @@ from astropy import units as u import WeatherRoutingTool.utils.unit_conversion as unit +import WeatherRoutingTool.algorithms.genetic.utils as gen_utils from WeatherRoutingTool.utils.maps import Map