A Comparison of Genetic Representations and Initialisation Methods for the Multi-objective Shortest Path Problem on Multigraphs

This paper compares different solution approaches for the multi-objective shortest path problem (MSPP) on multigraphs. Multigraphs as a modelling tool are able to capture different available trade-offs between objectives for a given section of a route. For this reason, they are increasingly popular in modelling transportation problems with multiple conflicting objectives (e.g., travel time and fuel consumption), such as time-dependent vehicle routing, multi-modal transportation planning, energy-efficient driving, and airport operations. The multigraph MSPP is more complex than the NP-hard simple graph MSPP. Therefore, approximate solution methods are often needed to find a good approximation of the true Pareto front in a given time budget. Evolutionary algorithms have been successfully applied for the simple graph MSPP. However, there has been limited investigation of their applications to the multigraph MSPP. Here, we extend the most popular genetic representations to the multigraph case and compare the achieved solution qualities. Two heuristic initialisation methods are also considered to improve the convergence properties of the algorithms. The comparison is based on a diverse set of problem instances, including both bi-objective and triple objective problems. We found that the metaheuristic approach with heuristic initialisation provides good solutions in shorter running times compared to an exact algorithm. The representations were all found to be competitive. The results are encouraging for future application to the time-constrained multigraph MSPP.


Introduction
There is substantial evidence [1][2][3][4][5] that modelling transportation problems as multigraphs offer benefits with regards to time, cost, environmental impact, and flexibility in multiple practical settings. The adoption of multiple parallel edges between pairs of nodes is rooted in the multi-objective nature of the problems. When there are multiple ways to traverse a section of the route, offering different trade-offs between the objectives, they should all be modelled in the optimisation process to find a set of Pareto-optimal solutions. The parallel edges in practice might correspond to different physical routes as in multi-objective vehicle routing problems [2,6], different modes of transport [3], or the same physical route traversed with different speed profiles such as in energy-efficient driving [4,5] and in the airport ground movement problem [1].
The multi-objective shortest path problem (MSPP) on multigraphs can be decomposed to two intertwined NPhard problems, the MSPP on simple graphs [7] and the fixed sequence arc selection problem (FSASP) [6]. In the FSASP, the sequence of nodes to be traversed is fixed, but there are multiple parallel edges between any neighbouring nodes in the sequence. Consequently, finding the exact solutions for a multigraph MSPP often requires an unacceptably long time, in particular, in real-world applications.
For the simple graph MSPP, metaheuristics and, in particular, genetic algorithms have been used with success [8][9][10] to obtain a good representation of the Pareto-optimal solutions in a short time. Genetic algorithms imitate the biological process of natural selection to find near-optimal solutions to computationally difficult problems. A key aspect of this process is the chosen representation method, that determines how the candidate solutions are encoded to the so-called chromosomes for the evolutionary process. The representation method has an important role in determining the search landscape and thus affects the effectiveness of the search [11]. Moreover, the representation also determines the kind of crossover and mutation operators that can be used.
Several different representation methods and genetic operators have been proposed for the simple graph MSPP. The multigraph MSPP, on the other hand, received considerably less attention, even though it is a relevant issue for a wide range of real-world problems and a natural extension to the body of research investigating the simple graph MSPP. To our knowledge, the only previous attempts at applying metaheuristics to the multigraph MSPP were in the context of multi-modal transportation [12,13]. Both of these studies used genetic algorithms and extend the direct variable length [14] representation. There are three other representations proposed for the simple graph shortest path problems, the direct fixed length [15], random key [16], and integer-valued priority [10] representations. In this paper, we investigate and extend all the four mentioned representations.
The representations introduced for the simple graph problem cannot be used without modification for the multigraph MSPP. The extension requires the incorporation of the choice between parallel edges into the chromosomes. Figure 1 illustrates an example of a solution path in a multigraph.
Heuristic initialisation techniques have been proven to be useful in different combinatorial optimisation problems [17][18][19][20], as they can lead to quicker convergence by starting the evolutionary process with an already high-quality initial population. This is done by locating promising areas of the search space through utilising a priori information about the specific problem. The initial population is then drawn partially or exclusively from those promising areas.
A common concern about such initialisation methods is that they might lead to premature convergence caused by a lack of sufficient diversity. This explains why there are relatively few studies utilising heuristic initialisation for the MSPP. The problem of premature convergence can be mitigated by introducing enough randomness into the population, while still preserving higher quality compared to a purely random population.
We employ a heuristic initialisation method utilising single objective search. This method was introduced for the case of direct representation in [9], and here, we extend it to the case of the other representations. Moreover, we also introduce a novel heuristic initialisation method based on the idea of discouraging detours, that is applicable to all the representations used.

Contributions
The main contributions of this paper: -Four representations that were originally proposed for simple graph shortest path problems are adapted and extended to the multigraph MSPP. Three of these representations have not been previously applied to the multigraph MSPP. -A novel indirect method of incorporating the choice of parallel edges into the chromosomes is proposed. This method makes it possible to adopt the priority-based representations to the multigraph problem with inhomogeneous numbers of parallel edges without the use of any repairing mechanism. The method can also be used with the direct representations. -The direct [12] and indirect ways of encoding parallel edges are compared for the direct variable length representation. -A heuristic initialisation method based on graph structure is proposed and shown to provide additional benefits when used together with the existing heuristic initialisation method [9]. -The representation and initialisation methods are compared on a diverse problem set using dominance-compliant quality indicators. It is observed that some representations are better suited for specific graph types.
This study extends our previous work [21]. The first main addition compared to [21] is the inclusion of more variants of different representations. Variants include (1) different crossover operators, (2) direct way of encoding parallel edges for the direct variable length representations, and (3) adapting an existing heuristic initialisation method [9] to priority-based representations. Second, a larger set of problem instances is used for the comparison, and the performance of the exact algorithm on these instances is described in more detail.

Structure of Paper
In the rest of the paper, "" describes the multigraph MSPP problem. Related work on representation methods for the simple and multigraph MSPP, and initialisation techniques are summarised in "Related Work". "Proposed Representations and Initialisation Methods" describes our approach to extend the main representations to the multigraph case, the genetic operators used, and the proposed heuristic initialisation technique. " Details" contains the implementation details of the algorithms, description of test instances, and the evaluation of the approximate solution sets. In "Results", numerical experiments and their results are presented. Finally, conclusion is drawn and future research directions are discussed in "Conclusion and Future Work".

Problem Description
The multigraph MSPP is defined by a multigraph network G = (V, E) and a multi-dimensional cost-vector associated with each edge in G. The network is assumed to be undirected in this paper. V = 1, … , n represents the set of nodes, and E represents the set of edges. Given that G is a multigraph, there might be multiple edges in E connecting the same nodes. Thus, an edge in the network is denoted by e = (u, v, i) , where u, v ∈ V and i is a parallel edge index, which differentiates between the edges with the same endpoints u and v. These edges are numbered as i (u, v), where i is starting from 1 up to l(u, v), where l(u, v) is the number of parallel edges between the two nodes u and v. There are two special nodes v O , v D ∈ V , the origin and destination nodes, respectively. The costs associated with each edge according to k objectives considered are given by a k-dimensional cost-vector cost(e) = (c 1 (e), c 2 (e), … , c k (e)) . For a valid path P between v O and v D , the corresponding cost-vector can be calculated according to Eq. (1).
We are looking for the minimum cost path between v O and v D considering all k objectives. The solution is a set of valid paths with non-dominated cost-vectors, i.e., the Paretooptimal solutions. A solution path P 1 is said to be Paretooptimal if there is no solution path that is at least as good as P 1 according to all k objectives and better according to at least one objective. (1)

Related Work
In this section, we review the genetic representation methods and operators proposed for simple graph shortest path problems and the multigraph MSPP. Initialisation techniques are also discussed.

Genetic Representations for Paths in Simple Graphs
In the single objective shortest path problem (SSPP) and simple graph MSPP, candidate solutions need to encode a sequence of neighbouring nodes in a graph. Several genetic representation methods have been introduced for these problems. These methods can be classified as direct representations and priority-based representations.

Direct Variable Length Encoding
The most straightforward way of encoding a path in a graph is the direct variable length representation proposed in [22] for the SSPP. Chromosomes consist of lists of node IDs, that form a path starting with the origin node. An arbitrary list of nodes usually will not correspond to a feasible path in the graph, and this necessitates the use of problem specific genetic operators, which maintain feasibility of the path. Ahn and Ramakrishna introduced a crossover based on common nodes in the solution paths and a mutation based on a random walk for the SSPP [14], which were adapted by multiple authors for the multiobjective problem [8,9,23].
The main advantage of this representation is that it gives a one-to-one mapping, which is usually preferable over one-ton mapping, since it avoids introducing plateaus in the search space. In one-to-n mapping, several different chromosomes might encode the same solution path, and thus have the same associated fitnesses, forming a plateau. The algorithm then might rediscover the same solutions over and over again through different chromosomes and waste computational resources.
A disadvantage is that the genetic operators might lead to loop formation, and thus offspring need to be checked and repaired after mutation and crossover. Also, according to [24,25], this representation is not suitable for large networks.
An example of a chromosome with direct variable length representation that encodes the solution path in Fig. 1 without specifying the choice between parallel edges is:

Direct Fixed Length Encoding
Another node ID-based representation was proposed by Inagaki et al. [15]. The length of the chromosome equals n, the [1,6,2,4].
number of nodes in the network, and the node IDs are the numbers from 1 to n. The chromosomes incorporate a pointer to a neighbouring node for each node. The locus of a gene corresponds to a node in the network with the same ID, and the value of the gene is the ID of a neighbour that should be the next node in the path. The path is decoded by following the pointers until the destination is reached or a loop is formed.
Note that this is an one-to-n representation, because, generally, there are genes whose values do not affect the decoding process. The crossover operator applied in Ref. [15] is essentially uniform crossover. This approach is deemed ineffective and requires large population sizes [14,24].
An example of a chromosome with direct fixed length representation that encodes the solution path in Fig. 1 without specifying the choice between parallel edges is:

Integer-Valued Priority-Based Encoding
Gen et al. proposed a priority-based encoding technique [26] that uses integer priority values to encode solution paths indirectly. Chromosomes consist of some permutation of the integers from 1 to n. The priority assigned to a node with ID i is given by the value of a gene at locus i.
A path is decoded from the chromosome by starting at the origin node and step-by-step moving to the neighbouring node with the highest priority, given it is not yet in the path. If it is already in the path, the neighbour with the next highest priority is chosen instead, if there is one.
The main advantage of this representation is that a random permutation of the priorities will always be decoded to some valid path starting from the origin node. This means that more traditional crossover operators can be used and expected to produce feasible paths, unlike in the case of the direct representations.
In Ref. [26], position-based crossover (PX) was used. An extension of the one-point crossover for integer-valued priority representation, the weight mapping crossover (WMX), was proposed specifically for this problem [10]. This approach has been shown to perform well in comparison to the direct variable length representation for the bi-objective problem.
An example of a chromosome with integer-valued priority representation that encodes the solution path in Fig. 1 without specifying the choice between parallel edges is:

Random Key-Based Encoding
Gen and Lin proposed another similar encoding technique using floating-point numbers instead of integers as priorities [27], resulting in a random key representation for the SSPP. Random keys were found to be a powerful method for permutation representation in other combinatorial optimisation problems [28].
The advantage of random key encoding is that we do not have to maintain distinct priority values, as most of the time, the value of two priorities compared will not be equal. This makes it possible to use different crossover operators that are not designed specifically for the MSPP, such as arithmetical crossover, uniform crossover, or two-point crossover. In Ref. [27], arithmetical crossover is employed, where the offspring are calculated as the weighted average of the two parent chromosomes. The authors report higher search capability, enhanced rate of reaching optimal solutions, and improved computation time compared to the integer-valued prioritybased encoding and the direct variable length encoding.
An example of a chromosome with random key representation that encodes the solution path in Fig. 1 without specifying the choice between parallel edges is:

Genetic Representations for Paths in Multigraphs
To encode solution paths in multigraphs, the indices of parallel edges used in the path also need to be specified.
Abbaspour and Samadzadegan extended the direct variable length representation to the multi-modal transportation problem [12] by indicating the mode of travel for each edge in the path. The chromosome is twice as long as the number of nodes in the network. The genes at odd loci correspond to node IDs, as in Ref. [14]. The genes at even loci are used to indicate the mode of travel between the nodes.
Yu and Lu [13] proposed a slightly different approach, where they only indicated the changes in the mode of travel and do not specify it for each edge separately. Genes indicating the mode of travel for the consecutive node IDs have a negative value, to differentiate them from the node IDs. However, only a limited number of parallel edges were used. When the number of parallel edges grows, the possible advantage of a shorter chromosome is diminished by the burden of maintaining the chromosome structure with more complicated operators.

Initialisation
Most mentioned works applied random initialisation, which implies a random walk starting from the origin node, random pointers to neighbours, or random priorities. There are two notable examples when heuristic initialisation methods were used, both with the direct variable length representation. (1) Information about the spatial location of the nodes in the SN Computer Science graph is utilised [23]. The initial solutions are generated in a way that each edge moves away from the origin node and closer to the destination node measured by Euclidean distance. However, spatial information is not always available and in some cases might be misleading. (2) Single objective search is utilised [9]. Shortest paths are found according to weighted aggregations of smaller subsets of the objectives with the Dijkstra's algorithm. These are included in the initial population along with randomly generated chromosomes.
A heuristic initialisation approach was proposed for solving the SSPP with particle swarm optimisation and priority-based representation in Ref. [25]. Similar to [23], this method aims to decrease the risk of loop formation, but it uses information about the structure of the graph instead of spatial information. Priorities were randomly assigned and node IDs were assigned in a way that the higher the IDs are, the closer a node is to the destination. Then, detours in solution paths can be detected through finding decreasing node IDs. Not all detours were prohibited, only the ones where the decrease in node IDs is above some pre-specified limit.

Encoding Parallel Edges
Representations for the multigraph MSPP need to encode the parallel indices for each pair of consecutive nodes in a solution path, not just the node sequence. The previously proposed methods for encoding mode of transport in the multimodal transportation literature employ the direct variable length representation, which allows for direct inclusion of the parallel edge indices. Extending the priority-based representations to the multigraph case is less straightforward.
The multigraph might contain varied numbers of parallel edges between pairs of nodes, as in the multi-modal transportation problem. This does not cause a problem in the direct representations, because once the parallel edge indices are initialised to be feasible, and the crossover and mutation operators do not ruin feasibility. This is not the case with the priority-based representations.
Following a similar strategy for the priority-based representations to assign parallel edge indices to node IDs may lead to infeasible solutions after genetic operators. Take an example of a crossover that results in a solution path in which a node u is followed by another node v. In both of the parent solution paths, the node u might be followed by nodes other than v. Then, the parallel indices assigned to node u in the parents might be higher than what is available between u and v.
One possible way of overcoming this problem would be to apply a repairing mechanism after the solution path is decoded. Apparently, the choice between parallel edges would be mostly governed by this repairing mechanism, not taking full advantage of the exploration and exploitation capabilities rendered by the mechanism of the evolutionary process. Also, the computational burden would increase. Here, we propose an alternative method for representing the parallel indices that does not require repairing.
Instead of encoding the indices of parallel edges directly, we use a floating-point number r between 0 and 1, denoted the parallel edge indicator. Parallel edge indicators are assigned to node IDs, and they encode which parallel edge to use when leaving the given node. Given two neighbouring nodes u, v and the number of parallel edges between them l(u, v), the index of the chosen parallel edge can be calculated as ⌊r(u) * l(u, v)⌋ + 1 when moving from u to v (indexing starts at 1). Any random value of parallel edge indicators can be decoded to an index of a parallel edge that is available between two given nodes. The parallel edge indicators are employed with all four representations as a way of extending them to the multigraph problem.
One possible drawback is the increased number of alternative chromosomes that can encode the same solution path, leading to the increased ambiguity of the representations. To investigate the effect of this ambiguity, the directly encoded parallel edge indices [12] are also implemented with the direct variable length representation. This provides a way of encoding solution paths in the multigraph without any ambiguity, providing an interesting case for comparison with all other representations. The directly encoded parallel edge indices are also applicable to the direct fixed length representation also. We did not include this in our investigation, because that representation already includes some ambiguity, as noted in "Direct Fixed Length Encoding".
The directly encoded parallel edge indices are also represented as floating-point numbers assigned to node IDs. The value of a parallel edge index is divided by 100 to get a floating-point number. Therefore, up to 100 parallel edges can be represented between two nodes.

Direct Variable Length Representation for the Multigraph MSPP
We implement two versions of this representation, abbreviated DirVarL-indir and DirVarL-dir. The difference is in the encoding of the choice between parallel edges, which is done indirectly with parallel edge indicators for DirVarLindir, and for DirVarL-dir, it is done directly with parallel edge indices. Note that DirVarL-dir is similar to the solution representation used in Ref. [12], but we use different genetic operators. [12] employed an expensive mutation operator that includes using the Dijkstra's algorithm and a repairing mechanism. Instead, we use the mutation suggested in [14], as it is most often used with the direct variable length representation.
The base of the chromosome is a sequence of node IDs, that are integers, such as in [14]. As the parallel edge indicators are floating-point numbers below 1, it is convenient to add them to the node IDs. This way, the chromosome is the same length as in the simple graph problem. The integer parts indicate the node sequence of the path and the fractional parts indicate which of the parallel edges to use when leaving the nodes.
The solution path can be decoded according to Algorithm 1. The algorithm loops through the chromosome (lines 3-13) and collects the nodes and indices determining the solution path. If the choice of parallel edges is encoded directly, lines 8-9 are executed, and if they are encoded indirectly, line 11-12 is used to decode the parallel indices.
An example of a chromosome with DirVarL-indir that encodes the solution path in Fig. 1 is: where the value of the first gene encodes that the first node in the solution path is 1 and the parallel edge to be used to reach node 6 is calculated as ⌊r(1) * l(1, 6)⌋ + 1 = ⌊0.38 * 3⌋ + 1 = 2 . An example of a chromosome with DirVarL-dir that encodes the solution path in Fig. 1

Direct Fixed Length Representation for the Multigraph MSPP
The parallel edge indicators are added to the genes to extend the representation from [15]. The solution path can be decoded according to Algorithm 2. The loop in lines 3-12 is executed until the destination node is reached by the path, or a node would appear a second time in the path. The next node to add to the path is found in line 5, and if it is not yet in the node sequence, the parallel index to use is decoded in lines 10-12 and added to the parallel index sequence.
We abbreviate the direct fixed length representation for the multigraph MSPP using parallel edge indicators DirFixL.
An example of a chromosome with DirFixL that encodes the solution path in Fig. 1 is:

Integer-Valued Priority-Based Representation for the Multigraph MSPP
The parallel edge indicators are added to the genes to extend the representation from [10]. Then, the priority value of node v can be found as the integer part of the vth gene in the chromosome, and the parallel edge indicator for node v can be found as the fractional part of this gene.
The node sequence and sequence of parallel indices can be decoded according to Algorithm

SN Computer Science
3-14 is executed until the destination node is reached by the path, or any node that only has neighbours that are already in the path. The next node to add to the path is found in line 10, and the parallel index to use is decoded in lines 12-14 and added to the parallel index sequence.
An example of a chromosome from this representation that encodes the solution path in Fig. 1 is:

Random Key-Based Representation for the Multigraph MSPP
Given that the priorities in this representation [27] are floating-point numbers, we cannot use the same strategy to include the parallel edge indicators. They are kept separately, making the genes two-dimensional. The first value of the gene at locus i encodes the priority value of the node with ID i, while the second of the gene at locus i encodes the parallel edge indicator of the node with ID i. The sequence of nodes and parallel edge indices can be decoded from a chromosome according to Algorithm 3.
An example of a chromosome from this representation that encodes the solution path in Fig. 1 is:

Variation Operators for the Representations
Different representations require different genetic operators. This section presents the genetic operators employed in this paper with each of the representations. For some representations, multiple crossover operators are considered, resulting in different variants of the same representations. The variants are summarised in Table 1 with example chromosomes that encode the same solution path for each representation.
There are two variants for the direct variable length representation, the difference being the method used for encoding the parallel edges. Both variants use the same genetic operators. Unlike the simple graph case, crossovers can be conducted on any two parents that have at least one node in common apart from the origin and destination node. Even if the resulting node sequence is the same as a parent node sequence, the parallel edges might change. The mutation operator generates new partial solutions by a random walk from a randomly chosen node in the path.
There is only one variant for the direct fixed length representation. Uniform crossover is employed. In mutation, the integer part of each gene is changed with probability 0.5. The integer part at locus i is changed to a random neighbour of the ith node. There is no need to change the fractional part, i.e., the parallel edge indicators, because they indicate a decision between a different set of parallel edges.  We include two variants for the integer priority representation, with two different crossover operators, PX [26] and WMX [10], which was specifically introduced for the MSPP. In both cases, insertion mutation is used, meaning that a randomly picked gene is removed from the chromosome and inserted back at a new locus. We abbreviate the integer-valued priority-based representation for the multigraph MSPP using parallel edge indicators paired with WMX IntPri-WMX and IntPri-PX when it is paired with PX.
We include three variants for the random key representation paired with three different crossover operators. In Ref. [27], arithmetic crossover was used. We also investigate twopoint crossover and uniform crossover. Insertion mutation is used in all these cases. Both the mutation and crossover mechanisms are independent of the values of the genes, and thus, they are straightforward to apply on the two-dimensional genes described in "Random Key Based Representation for the Multigraph MSPP". We abbreviate the random key-based representation for the multigraph MSPP using parallel edge indicators RanKey-arithX, RanKey-uniX, or RanKey-2ptX when it is paired with arithmetic, uniform, or two-point crossover, respectively.

Heuristic Initialisation Based on Graph Structures
Inspired by [25], we propose a novel heuristic initialisation technique that aims to discourage detours (HeurI1). The method incorporates knowledge about the network structure and randomisation to provide a diverse initial population of high quality. The difference is that we incorporate the roles of IDs and priorities into a single priority value assigned to each node in a semi-random way. The idea is to give higher priorities to nodes closer to the destination node with higher priorities than to the nodes far from it, thereby discouraging detours.
The method is first introduced for random key representation. Then, the initialised solutions can be easily converted to the other three representations. The priority p is assigned to node v according to (2). The hopcount (the minimum number of edges in a path) between node v and node v D in G is denoted h(v, v D , G) , and max denotes the maximum value of the randomisation coefficient : The likelihood of detours appearing in the decoded paths can be controlled by the parameter max . The higher max is, the more random the priorities are, and the less prominent is the effect of the heuristic initialisation compared to a purely random one. With a value of max > 1 , small detours are possible, a path might move from a node to another one with the same hopcount, as depicted in Fig. 2. If max > 2 , moving to a node that is at higher hopcount from the destination is possible and becomes more probable with the increase of max . The value of max can be optimised with the rest of the parameters for the experiments.
The resulting values need to be transformed according to the representations before they are fed to the algorithms as initial populations. For random key encoding, they need to be normalised to fit the appropriate intervals. For integervalued priority encoding, the priorities are converted to integers by sorting them into increasing order and assign to each node the rank of its priority value. For direct encodings, the priorities are converted to node-based representations. This way, the method can be used with all four representations. The parallel edge indicators are assigned randomly.
A crucial point is that the heuristic initialisation method should be easily computable compared to the original problem being solved. The proposed method makes use of the hopcount of each node in the graph from the destination node, which can be computed in O(V + E) time. This is significantly lower than solving the multigraph MSPP, which is in general NP-hard.

Heuristic Initialisation Based on Single Objective Search
The second heuristic initialisation method (HeurI2) is adapted from [9], it returns an already good solution by using single objective search (Dijkstra's algorithm). We use this method to initialise five solutions for the bi-objective problem with weights for the objective values equally distributed in the interval (0, 1). For the triple objective problem, we initialise seven solutions with this method: one with each single objective, one with each pair of objectives with equal weights, and one considering all objectives with equal weights. The rest of the candidates are initialised using the HeurI1 method, or a purely random initialisation. Because the difference in the hopcounts of nodes 2 and 4 from node 5 is smaller than the difference of the random numbers associated with these nodes, a detour is formed.

SN Computer Science
The HeurI2 method is most conveniently used with the direct variable length representation, as in Ref. [9], because Dijkstra's algorithm returns a solution path. Here, we also use it with the other three representations. To do this, we first need to translate the solution path with the parallel indices to the respective representations. There are multiple ways to do this, because only the direct variable length representation provides one-to-one encoding.
In Algorithm 4, we describe one of the possible methods for producing an integer priority-based chromosome that encodes a given solution path. In the first stage, the priorities of the nodes that appear in the path are set (lines 5-10). These priorities are increasing from the destination node towards the origin node. The rest of the genes are filled up with the remaining priority values, and random parallel edge indicators in lines (13)(14)(15). This way, at any point in the decoding process, the neighbour that follows the current node in the path has the highest priority in the whole graph except for nodes already in the path. Parallel edge indicators are set to encode the required parallel indices between nodes in the solution path (lines [6][7][8]. The translation to the remaining two representations is done along the same lines.
that corresponds to the ith node in nodeSequence with

Implementation Details
NSGA-II [29] is employed to compare the representation and initialisation methods. The selection and elitism mechanisms are defined by NSGA-II. This algorithm scales well for two and three objective problems [30]. All numerical tests are performed on Queen Mary's Apocrita HPC facility [31]. The methods are implemented in Python, for the NSGA-II implementation, the inspyred package [32] was used. The parameters for the NSGA-II and the initialisation were tuned with the use of the irace package [33], separately for all the variants of the representations. The tuned parameters are shown in Table 2.
The fitness of a valid path is calculated according to (1). It might happen that some candidates encode paths that do not reach the destination node. In these cases, a penalty function is used, which assigns a large cost to such candidates. The penalty is larger the further away the path ends from the destination node, measured by hopcount. The fitness of an infeasible path P ′ that does not reach the destination node is calculated according to Eq. (3), where cost max is the k-dimensional vector where each component equals the maximum value of any cost component in the given instance:

Test Instances
The algorithms are evaluated using 32 test instances, 16 for the 2 objective problem, and 16 for the 3 objective problem. These instances differ in the graph type, the maximum number of parallel edges in the multigraph, and the correlation between the objectives.
We use Waxman networks [34] with 100 and 196 nodes and square grid networks with the same number of nodes (10 by 10 and 14 by 14 nodes). The origin and destination nodes are specified as two endpoints of a diameter (largest hopcount) of the network, to ensure that they are not too close, to avoid setting a trivial problem. Each edge in these simple graphs is converted to a multi-edge by assigning a cost matrix to it. The cost matrix contains a cost-vector in each of its rows, corresponding to a parallel edge. The number of rows of these cost matrices is randomly chosen between 1 and l max , where l max the maximum allowed number of parallel edges, 5 or 10 in this case. All parallel edges between the same two nodes have non-dominated cost-vectors.
We included instances with uncorrelated objectives and with negative correlation between the objectives. Negative correlation is the case where a multi-objective approach is essential for real-world applications and also these kinds of instances are the most challenging for exact solution approaches [35].
The cost assignment method described in Ref. [35] was used, to generate pairs of cost components with a specified negative correlation multiplier . The method is described by (4), where the first cost component ( C 1 ) is randomly generated from a uniform distribution within the interval specified by C min and C max . (C 2 ) * is a randomly generated value from the same interval, and C 2 is the second cost component. In our instances, the interval is specified by C min = 10 and C max = 1000: Note that for positive correlation multipliers, there is a separate formula [36]. However, we do not include instances with positive correlation, because they are generally easier to solve by exact approaches. In the case of three objectives, the third cost component is also calculated from the first cost component with the same value of .

Evaluation of Approximate Solutions
The proposed representations and their variants (listed in Table 1) are tested empirically and their performances are compared to a reference front using a set of quality indicators. The reference front is the true Pareto front, when it is available, and an approximation of it otherwise.

Reference Fronts
The true Pareto front was found by a state-of-the-art exact algorithm NAMOA* [37], a multi-objective variant of the A* algorithm, which was adapted to the multigraph problem. This algorithm uses heuristic functions to speed up the search. Here, we used a heuristic function proposed in [38] defined as, h TC (n) = (c 1 (n), c 2 (n), … , c q (n)) , where c i (n) is the optimal scalar cost of a path from node n to the destination node, considering only the ith cost component.
We allowed 1 day for the execution time of NAMOA*. When this was not enough, we approximated the Pareto (4) front using the solutions already found by NAMOA*, if any, and the approximate solutions returned by variants of NSGA-II. There are some bounds on the quality of Pareto fronts approximated this way. We know that they contain at least as many members of the true Pareto front as the number of objectives, because solutions initialised with the HeurI2 method are included. When NAMOA* is stopped prematurely, it either does not return any solutions, or it returns a subset of the Pareto-optimal solutions, these are also included in the reference front when they are available. Thus, all other members of the approximate front are nondominated by at least some members of the true Pareto front.

Quality Indicators
We use the multiplicative Epsilon indicator, the R3 indicator, and the relative hypervolume (RHV) indicator to evaluate the approximations of the Pareto fronts. A detailed description of these indicators can be found in Ref. [39]. The use of these metrics is in line with the recommendations in [40]. All three indicators signal higher qualities of the approximation front by lower values. When the approximate front is fully converged to the real Pareto front, the R3 metric and the RHV indicators have a value of 0, and the Epsilon indicator has a value of 1. Table 3 presents the experimental data regarding solving the 32 test instances using NAMOA*. In particular, it shows the high execution times of NAMOA*; for grid graphs, this time often exceeds 24 h, and overall, it is only twice below 10 s.

Finding Exact Solutions for the Multigraph MSPP Instances
We can see the size of the reference front (true Pareto front or the approximated Pareto front). According to both running time and reference front size, the bi-objective instances on Waxman graphs can be considered the easiest and the instances on grid graphs with three objectives the most difficult. The last two columns also show that if NAMOA* is stopped prematurely after 10 or 20 s, it does not return any solutions in the case of most test instances. When it does return solutions after this short time, we can see that they are rarely a good approximation of the reference front, as measured by the RHV quality metric. These results are compared to NSGA-II in " Differences Across the Instances".

SN Computer Science
In the following, we present results about using NSGA-II with the different representations and show that high-quality approximate solutions can be found in only 10 s.

Initialisation Methods
Four different initialisation procedures were used.
We combined all eight variants of NSGA-II with each initialisation procedure and compared the achieved solution qualities according to the quality indicators in Table 4. We can see that H1 + H2 and H2 + R are visibly better than the Table 3 The performance of NAMOA* on the 32 problem instances. The first five columns describe the problem instances The column "NAMOA* time" specifies the running time of the NAMOA* algorithm on the given instances in seconds, if this time is less than a day. The column "Reference front size" specifies the number of solutions in the Pareto front, or when the algorithm does not finish in a day, the number of solutions in the approximation of the Pareto front. The last two columns describe the solution quality measured by relative hypervolume reached by NAMOA* when it is stopped prematurely after 10 or 20 s other two in all cases. One main reason for this is that the H2 method finds the extreme solutions in the initialisation step, while these are often not found at all using a method without H2. Apart from the direct variable length representation, the H1 + H2 method significantly outperforms the H2 + R. On the other hand, H2 + R is significantly better than H1 + H2 only in case of the DirVarL-dir variant. The p values for the statistical tests are indicated in Table 4. The general patterns of the results are similar in case of the 20 s time budget. For details, see the Appendix Table 6.
One possible explanation for why the H1 + H2 initialisation method is not better than the H2 + R method only in case of the direct variable length representation is the following. In the case of the other three representations, the H1 method encodes some extra information about which direction to choose at each node in the graph, not just the nodes in the encoded path. This is possible because of the inherent ambiguity of these representations. If such a chromosome is modified through mutation or crossover, the extra information might make it more probable that the resulting new chromosome will be a valid path. While the direct variable length representation encodes the path without ambiguity, and thus, the extra information provided by the H1 initialisation is lost when the chromosome is converted to this representation.

Comparison of NSGA-II Variants
In the following, only the best initialisations are considered for each variant. Table 4 also shows that, on average, the RanKey-2ptX variant reaches the best results (with H1 + H2 initialisation) considering all instances. This is confirmed by the one-sided Wilcoxon signed-rank test with p values below 10 −7 for all three quality indicators. This also holds for the 20 s time budget, see Appendix Tables 6 and 7.

Performance of Variants on Average Across All Instances
These results might be surprising given that RanKey-2ptX has the most ambiguity out of the four representations. The results suggest that a different search landscape and more effective genetic operators can balance the increased ambiguity of representations. The extra information encoded in the chromosomes with some ambiguity-as discussed in "Initialisation Methods"-might also contribute.
In the following, we identify the best variant for each of the four representations, based on the data presented in Table 4 for the further analysis.
There were two variants considered for the direct variable length representation. The DirVarL-dir variant was better according to the averages of all three quality indicators than the DirVarL-indir. The one-sided Wilcoxon signed-rank test Table 4 Illustration of the role of the heuristic initialisation methods in the solution quality achieved by variants of NSGA-II using different genetic representations and crossover operators, described by the R3 metric The lowest values of the quality indicators corresponding to the best initialisation technique are highlighted in bold in each sub-column Lower values correspond to better quality. The average values are calculated from 50 runs for each of the 32 instances. The time budget for NSGA-II is specified as 10 s. The one-sided Wilcoxon signed-rank test was used to decide statistical significance between H1 + H2 and H2 + R. Results are indicated as (*)p < 0.05 , (**)p < 0.005 , and (***)p < 0.0005 confirms the statistical significance in all cases with p values below 10 −35 for both values of the time budget. This difference confirms the expectations that increasing the ambiguity of the representation of the parallel edge indices decreases the effectiveness of the search to some degree. There were two variants considered for the integer priority representation. The IntPri-WMX variant outperformed the IntPri-PX variant according to all quality indicators on average. The one-sided Wilcoxon signed-rank test confirms the statistical significance in all cases with p values below 10 −22 for both values of the time budget. This is not surprising given that WMX was specifically introduced for the shortest path problem, while PX not.
There were three variants considered for the random key representation. The RanKey-2ptX variant outperformed both the RanKey-arithX and RanKey-uniX variants among all others, as we have already seen in the overall comparison of the variants.
Taking the average of the quality indicators across all instances as the basis of comparison might cover up important differences between the representations. It can be the case that the RanKey-2ptX variant is consistently outperformed by another variant for a subgroup of the instances. This is investigated next.

Performance Differences Across the Instances
In the following, we only consider the best variant for each representation with the best initialisation method. The four representations are compared in detail on the 32 instances separately. Here, the representations are compared with a time budget of 10 s. For the sake of comparison, the solution quality achieved by NAMOA* in 10 s is also included. The RHV values are shown, because this quality indicator has the highest correlation with the other two. Note that the value of the Spearman's rank correlation coefficient is above 0.91 between any two of the three quality indicators, and this suggests that using any of the three would result in roughly the same order between the evaluated variants.
In Table 5, we can see that while RanKey was shown to be the best on average in "Performance of Variants on Average Across All Instances", in fact, three of the representations are shown to be competitive when compared separately for each of the instances. IntPri-WMX is outperformed by the other three in the majority of cases. DirVarL-dir seems to be particularly suitable for grid instances with three objectives, while generally less successful in case of Waxman networks and grids with two objectives. RanKey-2ptX shows an overall good performance. However, the grid instances with more than two objectives seem to be a weak point.
We can also see in Table 5 that NAMOA* outperformed NSGA-II with the 10 s time budget in only two cases, when it found the whole Pareto front within 10 s. In other cases, even if NAMOA* returns some solutions, the quality is magnitudes lower than with the NSGA-II variants. Thus, we can conclude that for these test instances, we should not choose NAMOA* as a solver algorithm if the time budget is small.
The general patterns of the results are similar with the other quality indicators and with the time budget of 20 s. For the detailed results, see Appendix Tables 6,7,8,9,10,11,and 12.

Conclusion and Future Work
A systematic investigation of design choices of genetic representations for the multigraph MSPP was presented in this paper. Four main representations were investigated, some with multiple variants. The resulting eight variants were compared using NSGA-II.
The proposed encoding schemes for incorporating the choice between parallel edges are through a floating-point number as parallel edge indicators, instead of including the parallel edge indices directly. This approach is necessary for the extension of priority-based representations to the multigraph MSPP when inhomogeneous numbers of parallel edges between pairs of nodes are possible.
Multiple different initialisation methods are also considered. Apart from purely random initialisation, two heuristic initialisation methods are investigated. One of these is an existing method based on single objective search that was previously only used with the direct variable length representation. We adapted it to be applicable to the other representations. We also propose a novel heuristic initialisation method that makes use of a priori knowledge about the graph structure and can be used with all representations.
The representations and their variants combined with different initialisation methods are compared according to three quality indicators. It is found that our proposed initialisation method provides additional benefits when used together with the other heuristic initialisation methods in three out of the four representations.
Regarding the average of the quality indicators across all instances, the best variant of the algorithm was the random keys representation with two-point crossover. This variant employs the novel heuristic initialisation proposed, and the parallel edge indicators to encode the choice between parallel edges, which proves the effectiveness of the proposed methods.
When breaking down the test instances further, it was revealed that for grid-based problem instances with more than two objectives, the best approach is the direct variable length representation, while the random keys representation is generally the best for the other problem instances. This might be because the crossover operator of the direct variable length representation requires a common node between the two parent paths, and this is more likely to exist in grid graphs. One future research direction is to understand this difference in depth, and design an algorithm that incorporates the strengths of both representations.
The solution qualities achieved this way gave a good representation of the reference fronts in only 10 s, while obtaining the exact Pareto front required significantly longer times for most of the investigated instances. The running times are often important in practical applications.
We saw that the ambiguity of representing the choice of parallel edges makes NSGA-II less effective. However, the inherent ambiguity of the random key representation in encoding node sequences does not seem to cause a problem. One future direction is to understand how can the representation with the most amount of ambiguity be the most successful.
Additionally, future work includes the extension to constrained problems, particularly the investigation of constraint handling methods. The extension will allow Table 5 Comparing the best variants of NSGA-II for the four representations, separately for the 32 test instances with a time budget of 10 s The best values among the representations are shown in bold in each row The values of the average RHV quality indicator of 50 runs with each representation is shown. The results of NAMOA* regarding the RHV quality indicator from Table 3  the proposed approaches to be tested on real-world problems and to be compared with other exact and approximate algorithms.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are  included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.