Nash equilibrium sorting genetic algorithm for simultaneous competitive maximal covering location with multiple players

In this article, a new genetic algorithm (GA), called the Nash equilibrium sorting genetic algorithm (NESGA,) is introduced to identify Nash equilibria for the competitive maximal covering location problem with two and three competitors, which is a combinatorial game theory problem where it is computationally intractable to enumerate all decision options of the competitors. Although GAs are widely used in combinatorial optimization, their applications to non-cooperative, simultaneous games have been limited owing to challenges in guiding the evolutionary search to Nash equilibria, which are not necessarily on the Pareto front of the search space and cannot be found by current multi-objective GAs. A novel fitness assignment strategy is proposed to help the NESGA to converge to multiple Nash equilibria. Computational experiments show that the NESGA can discover multiple Nash equilibria in a single run and outperform other game-theoretic GAs.


Introduction
Facility location problems are among the classical operations research problems with the broadest range of application areas. One of the most widely studied facility location problems is the maximal covering location problem (MCLP) (Church and ReVelle 1974), where a decision maker seeks to cover as much customer demand as possible by opening multiple facilities, selected from a set of candidate points. Customer demands are concentrated at a set of discrete points, and a customer point is covered if there exists at least one facility within a critical distance from the customer point. The MCLP has extensively been studied in the literature, and several review articles (Drezner and Hamacher 2001;Farahani et al. 2012;Melo, Nickel, and Saldanha-da-Gama 2009;Snyder 2006) summarize the existing work. The competitive version of the problem, where two or more competitors try to maximize their demand coverage in the same target market, is applicable when a firm determines facility locations in a competitive market. Regarding the nature of competition, the competitive MCLP can be modelled in two ways: simultaneous or sequential (Kress and Pesch 2012). In simultaneous models, the competitors make decisions without being aware of one another's decisions. It is assumed that the competitors have the knowledge of all possible decisions that their opponents may make, and they maximize their expected payoff by choosing the decision that is the best response to their opponents' possible decisions. Simultaneous games are useful for studying competitive behaviours as well as analysing strategies that would emerge as a response to market conditions, rules and principles.
CONTACT Sadan Kulturel  Supplemental data for this article can be accessed at https://doi.org/10.1080/ 0305215X.2021.1957861 Sequential models, which are more prevalent in the literature, are based on a Stackelberg game in which a competitor (i.e. leader) first locates its facilities, and then the second one (i.e. follower) locates its facilities optimally based on the leader's decision. This NP-hard problem is usually formulated as an integer bi-level optimization problem. Sequential models are strategic decision-making tools at the organizational level compared to the macro-level use of simultaneous models. Several heuristics (Alekseeva and Kochetov 2013;Davydov, Kochetov, and Carrizosa 2014;Seyhan, Snyder, and Zhang 2018) and exact approaches (Alekseeva, Kochetov, and Plyasunov 2015;Roboredo and Pessoa 2013) have been proposed to solve the sequential competitive MCLP.
In this article, the simultaneous version of the competitive MCLP is studied with two different demand capture models and two/three competitors, which has not been studied in the literature. A new approach, called the Nash equilibrium sorting genetic algorithm (NESGA), is proposed to find multiple Nash equilibria (NE 1 ) in the non-cooperative, simultaneous competitive MCLP and extend the problem to three competitors for the first time. In a multi-player, non-cooperative, simultaneous game, a decision profile is called NE if no player can improve its payoff by changing its decision option independently while the other players keep theirs unchanged. Simultaneous game theory problems are challenging for evolutionary algorithms because a problem may have no NE or multiple ones that are unlikely to be located on the Pareto-optimal front of the objective space. This aspect of non-cooperative simultaneous games makes it difficult for the existing single-/multi-objective evolutionary algorithms to converge to an NE. The current game-theoretic genetic algorithms (GAs) have not been applied to non-cooperative simultaneous games with three competitors. The NESGA utilizes a sorting-based fitness assignment strategy and multiple populations to discover multiple NE in a single run. The proposed sorting-based fitness assignment strategy can handle problems with multiple competitors.
The remaining sections of the article are organized as follows. Section 2 presents a review of the competitive MCLP and current game-theoretic evolutionary algorithms. Section 3 presents the mixed integer programming (MIP) formulations of the competitive MCLP with the closest-take-all and Huff gravity-based demand capture models. Section 4 presents the NESGA and its fitness evaluation method. In Section 5, it is shown that the NESGA can converge to multiple NE in a single run.

Background and previous work
Simultaneous games, where the players make their decisions simultaneously without knowledge of the other players' decisions, are challenging to solve by evolutionary algorithms owing to the limitations imposed by finite population size (Adami, Schossau, and Hintze 2016). One of the challenges is identifying and maintaining NE during the search since a problem may have no NE, a single one or multiple ones. Sefrioui and Perlaux (2000) proposed an algorithm called NashGA for non-cooperative simultaneous games. In the NashGA, each player has its population representing a set of decisions of the player, and the best decision of a population at the current iteration is used to evaluate the decisions of the other populations in the next generation. Pace (2009) also proposed a multi-population GA to discover dominated strategies for the Traveller's dilemma game, where the fitness of each solution in a population is evaluated against several randomly selected solutions from the other populations. Chen et al. (2002) used a negotiation policy to find the equilibrium in the payoff matrix of two players. Konak, Kulturel-Konak, and Snyder (2015) proposed the game-theoretic GA, where each individual of the leader's population is evaluated against every individual of the follower's population.
Several articles in the literature present a review of the competitive MCLP. Plastria (2001) introduced a taxonomy of the competitive MCLP in terms of the competition type, the market structure and the decision space for static problems. Kress and Pesch (2012) provided an extensive review of the sequential competitive MCLP. Alekseeva and Kochetov (2013) also presented a review of Stackelberg models and solution methods for the competitive MCLP. Eiselt, Marianov, and Drezner (2019) reviewed the latest trends in the competitive MCLP. The analysis of these review articles can be summarized in two main points: (1) the literature mainly focuses on sequential games, in particular, Stackelberg-type games with two competitors; and (2) the existing approaches to the problem utilize metaheuristics or combinations of metaheuristics and mathematical programming methods since the problem is p 2 complete (Noltemeier, Spoerhase, and Wirth 2007;Plastria and Vanhaverbeke 2008). The competitive MCLP was first defined by ReVelle (1986) for determining the location of a single facility considering the existence of other competitor facilities (although similar ideas date back to Harold 1929). This formulation represents a static view of the problem. Serra and ReVelle (1994) extended the static problem defined by ReVelle (1986) as a Stackelberg game and described an iterative heuristic procedure such that the leader moves the location of a single facility to an unassigned location, and the follower solves the MCLP problem for the leader's new solution. Plastria and Vanhaverbeke (2008) considered MaxMin, MinRegret and Stackelberg-type objectives and presented integer programming formulations that could find exact Stackelberg equilibria when the follower had only a single facility to open. Later, Seyhan, Snyder, and Zhang (2018) extended this approach to multiple follower facilities, where the follower's problem is assumed to be solved by a greedy algorithm (Kuehn and Hamburger 1963).  defined the competitive MCLP on the Euclidian plane, where the follower problem is solved with an exact approach and the leader decision options are chosen by a variable neighbourhood search heuristic. Küçükaydin, Aras, and Kuban Altınel (2011) introduced a bi-level mixed-integer nonlinear formulation to determine the leader and follower decisions considering the Huff gravity-based model. Barmak et al. (2016) and Nasiri et al. (2018) used a hybrid of heuristics and exact methods to solve a capacitated competitive MCLP.
In the context of the MCLP, a simultaneous game was first considered by Labbé and Louis Hakimi (1991), who studied the conditions for the existence of an NE for two competitors that aim to place a single facility on a network. Rhim, Ho, and Karmarkar (2003) analysed pure-strategy NE for an integrated facility location problem where competitors decide not only the locations to open facilities but also the capacity and pricing options of each facility. Godinho and Dias (2010) defined an iterative algorithm to solve a two-competitor simultaneous MCLP. Their algorithm starts from a feasible decision and sequentially computes the best response of each competitor to the last decision of the other until a pure NE is found, or a cycle of decisions is observed. Guo and Lai (2017) presented a simultaneous price and location game in which bricks-and-mortar stores compete with an online retailer. Cardinal and Hoefer (2010) extended the set covering problem on networks to noncooperative games with cost-sharing among multiple players and identified NE in facility location games where players offer payments for connecting and opening existing facilities. Another related problem is to identify a saddle point in a continuous plane for the location of a new facility under the gravity-based Huff model (Drezner and Drezner 2004;Drezner 2014;Díaz-Báñez et al. 2011).
The competitive MCLP was also formulated as a multi-objective problem (Lančinskas and Żilinskas 2013;Redondo et al. 2015;Fernández and Tóth 2009;Konak, Kulturel-Konak, and Snyder 2017;Rohaninejad et al. 2017) or considering cooperation among competitors (Mallozzi 2011). In these approaches, the payoffs of two competitors are simultaneously optimized using a multi-objective GA to identify Pareto-optimal decision profiles. Note that the objective of the NESGA is to discover NE in non-cooperative simultaneous games where an NE does not necessarily reside on the Pareto-optimal front (such examples can be found in the decision profiles given in the Supplementary Material). A Pareto-optimal decision profile is also not necessarily an NE. In this sense, finding an NE for non-cooperative games is more challenging than discovering NE on the Pareto-optimal front, which provides a target to evaluate alternative decision profiles.

Competitive maximal covering location problem
There is a set of competitors, denoted by set K, targeting the demands of n customer points distributed over a region. The candidate facility locations and customer points are represented by sets I and J, respectively. Facilities are opened at the customer points in this article (i.e. I = J). Distance d ij represents the Euclidian distance between two points i and j, and customers can only access facilities located within a threshold distance of D. Each customer point j has a deterministic demand of w j . Each competitor k aims to maximize its total demand capture π k by opening α k number of facilities. Customers are indifferent among the competitors and locations. The only factor that determines the decision of a customer is the distance to the facilities.
Let X k = [x k1 , x k2 , . . . , x kn ] represent a decision option of competitor k such that x ki = 1 if competitor k locates a facility at point i, and x ki = 0 otherwise, and k is the set of all decision options for competitor k. In the model, the competitors can co-locate their facilities at the same point. Let z kij be a binary variable such that z kij = 1 if point j is served by competitor k's facility located at point i, and z kij = 0 otherwise. The values of variable z kij depend on where the competitors locate their facilities as well as the assumptions about customer behaviour. Two demand capture models, namely 'closest-take-all' and 'Huff gravity-based' models, are considered.
The objective is to determine a decision for each competitor such that no competitor will benefit from changing its decision unilaterally while the other competitors keep their decisions unchanged. In other words, a set of NE decision profiles [X * 1 , X * 2 ] and [X * 1 , X * 2 , X * 3 ] is aimed to be identified for problems with |K| = 2 and |K| = 3, respectively, such that X * k is the best decision option of competitor k for the given best decision options of the other competitor(s). Let s represent a decision profile (e.g. s = [X 1 , X 2 , X 3 ]) and s −k denote a partial decision profile excluding competitor k (e.g. s −1 = [X 2 , X 3 ]), and π k (X k , s −k ) be the payoff of decision option X k for a given partial decision profile s −k .

Closest-take-all model
In the closest-take-all demand capture model, a customer always prefers the closest accessible facility, i.e. all demand of a customer point is captured by the closest accessible facility. If there exist multiple accessible facilities that are located at the same closest distance to a customer point, then the demand of such a point is shared equally by those facilities. Below, an MIP formulation is introduced for threecompetitor games by extending an earlier model (Konak, Kulturel-Konak, and Snyder 2017) to verify whether a candidate NE returned by the NESGA is a true NE or not. Additional notation for the formulation is given as follows: A(j) set of points that can cover customer point j, set of points that are closer to point j than point i, Problem MCLP (closest-take-all): subject to: Although the objective functions of the above model are nonlinear owing to the interactions among the competitors, the model becomes linear when it is solved only for a competitor k with a given s −k . Constraint (2) makes sure that point j can be covered by at most one facility of competitor k. Constraint (3) requires that competitor k locates a facility at point i if this point covers any customers. Constraint (4) sets the number of facilities that competitor k opens. Constraint (5) states that competitor k cannot cover customer point j through point i if there exists a facility closer to point j than point i.

Huff gravity-based model
In this model (Huff 1966), the probability that a customer prefers a specific facility is proportional to an attractiveness measure of the facility and inversely proportional to the distance. Let w kij represent the amount of demand w j captured by competitor k's facility located at point i. Assuming that all competitors and points have the same level of attractiveness, the ratio holds between two facilities located at points i and l for i = j, l = j, d ij ≤ D and d lj ≤ D, where λ ≥ 1 is called the distance decay parameter. Note that the above ratio is also valid for the facilities of the same competitor (i.e. k = r). The Huff gravity-based model implicitly assumes that w j = k∈K i∈I w kij z kij if point j is covered by at least one facility. Since the customers have no preferences, the facilities located at point j share w j equally. Below, an MIP formulation is introduced to verify whether a candidate NE is a true NE or not for the Huff gravity-based model.
Problem MCLP (Huff gravity-based model): subject to: Constraint (8) sets the number of facilities that competitor k opens. Constraint (9) states that competitor k must serve node j through its facility located at node i if no facility is located at node j. Constraint (10) makes sure that competitor k can serve node j through node i only if the competitor has a facility located at node i. Constraint (11) forces competitor k to serve node j if the competitor has a facility located at node j. Constraint (12) states that competitor k cannot serve node j through any other node i if node j has a facility located there. Constraints (13) and (14) are, respectively, upper bounds on the amount of demand w j captured by competitor k's facility located at point i and the total demand captured from node j. Constraint (15) ensures that the total demand of node j is captured if node j is served by any facilities. Constraints (16) and (17) model the demand capture-distance relations given in Equation (6). Constraint (18) calculates the demand capture from node j through a facility located at node j. The model becomes linear when it is solved for one competitor at a time to verify a candidate NE.

Description of the NESGA
In the NESGA, each competitor k ∈ K has its own population P k = [X k1 , . . . ,X kμ ] of μ decision options, where X kq = [x k1q , . . . ,x knq ] denotes the qth decision option of competitor k. The NESGA aims to converge the whole P k toward a subdecision space that includes NE for competitor k. To describe the NESGA, the game theory notation defined before by replacing k with P k is used, as follows: X k a decision option in P k X kq the qth decision option in P k X * k a decision option that is a part of a candidate NE in P k S the set of all decision profiles for all competitors at an iteration (e.g. S = P 1 × P 2 × P 3 ) S −k the set of all partial decision profiles that exclude P k at the current iteration S * −k the set of all partial candidate NE decision profiles that exclude P k s a decision profile in S (e.g. the subset of S * −k such that X * k is a part of a candidate NE. π k (X kq , s −k ) payoff of decision option X kq with respect to a partial decision profile s −k .
The following definition is used in the calculation of fitness values.
Definition 2: A decision profile s * = [X * 1 , X * 2 , X * 3 ] is called a candidate NE during an iteration of the NESGA if π k (X * k , s * −k ) ≥ π k (X kq , s * −k ) ∀X kq ∈ P k holds for each competitor k. The set of all candidate NE decision profiles at the current iteration is denoted as N.
In a generation, the NESGA compares each decision option in a population to every decision option in the other populations to construct a partial payoff matrix. Because multiple copies of a decision option are irrelevant in the payoff matrix, the NESGA's populations include only unique individuals. The overall procedure of the NESGA is given below, where g max is the maximum number of generations.

Fitness assignment
The fitness assignment of the NESGA is designed to maintain candidate NE in the current populations. Two concerns need to be addressed in the assignment of fitness values: (1) Assigning fitness values when the current populations do not include any candidate NE so that the populations retain decision options that are likely to produce promising offspring in future iterations.
(2) If there are multiple candidate NE at the current iteration, how to retain them for future generations and assign fitness values to other decision options.
To address these concerns, a hierarchical approach is proposed to rank decision options in each population. The first component of the fitness function depends on whether a decision option is part of a candidate NE. If a decision option X kq is a part of the candidate NE in P k , then its maximum payoff across the decision profiles in which it is a part of NE is used as its first fitness value: The second component of the fitness function depends on the existence of a candidate NE in the population, as follows: If there is no current candidate NE in the populations (i.e. N = Ø), then the maximum possible payoff of a decision option is used as its second fitness value. In other words, the maximax rule, which represents a 'risk-taking' behaviour, is used. If the populations have candidate NE (i.e. N = Ø), then the fitness values of the other decision options are evaluated considering only the set of the candidate NE. In this case, the maximax rule is also used but considering the payoff of the decision options only with respect to the candidate NE, not the entire payoff matrix. This approach intensifies the search around the decision options that are part of a candidate NE in each population and tests their fitness. After evaluating both fitness components, the populations are sorted in the descending order of Equation (19), breaking ties by Equation (20). Thereby, the NESGA maintains candidate NE over the generations and returns set N at the termination of the search.

Solution encoding and crossover/mutation
The NESGA can be used with any encoding schema and evolutionary operators. In this article, basic GA operators were utilized in the NESGA so that its performance, especially its fitness assignment strategy, could be studied without the effect of advanced encoding schemas, solution improvement methods, constraint handling techniques and parameter optimization. To compare the performance of the NESGA to the NashGA under the same conditions, both algorithms used a random key encoding (Bean 1994) and the same crossover/mutation operators. Let R kq = {r k1q , . . . ,r knq } represent the chromosome of the qth decision option in P k such that r kiq ∈ [0,1] is the random key corresponding to decision variable x ki . Given R kq , the decoding algorithm sets x ki = 1 for the largest α k random keys and x ki = 0 for all others to construct a decision option X kq . Thereby, the decoding algorithm ensures that a decision option X kq always has α k facilities opened.
The NESGA used a uniform crossover to produce an offspring from two randomly selected parents such that the offspring randomly inherits each gene from either of the parents with a probability of 0.5. The mutation involves assigning a random number to a gene with a probability of 1/n. Identical solutions are not allowed in the population. Hence, if a generated offspring already exists in the population, it is discarded, and a new one is generated. The procedure of the crossover operator is given as follows: Randomly uniformly select two decision options, a and b, from P k

Computational experiments
In the computational experiments, it is investigated whether the NESGA can converge to true NE, and its performance is compared to the NashGA (Sefrioui and Perlaux 2000). To demonstrate that the NESGA can discover multiple NE in a single run, competitive MCLP instances with known NE are needed. Therefore, new test problem instances with n = {30,60,120} were created, and their true NE were identified through enumeration for small-size instances. The data of the test problems and their attributes are given in the Supplementary Material. The problem instances have different values of parameters α k , D and λ, as given in the following tables.
If enumeration was not feasible, the MIPs given above were used to verify whether a candidate NE s * with the payoff value π k (X * k , s * −k ) was a true NE or not, as follows. First, the optimal payoff value π * k of s * for each competitor k was found using the MIP model after fixing the decision variables of the other competitors according to s * −k . Then, the following relative error (RE) metric was calculated.

RE(s
By the definition of the NE, if RE(s * ) = 0, then the candidate s * is a true NE. The mean relative error (MRE) for all candidate NE returned in a run can be calculated as follows: The NESGA and NashGA were coded in MATLAB ® R2018a. The enumeration and optimization algorithms to test equilibrium were coded in AMPL, and the resulting MIPs were solved by CPLEX 12.0 on the Pennsylvania State University's Institute for Computational and Data Sciences' Roar supercomputer. All MATLAB codes were run on a personal computer with a 4-core Intel i5 @ 3.5 GHz processor.

Computational experiments with the closest-take-all model
First, small problems with different α k values were used, as given in Table 1, to investigate whether the NESGA and the NashGA could discover NE, especially in problem instances with multiple NE. An exhaustive enumeration was used to determine NE for n = 30 and |K| = 3. The MIP model was also used to find the set of NE decision profiles for other problem instances with |K| = 2 where exhaustive enumeration was impossible, as follows. First, 1 was enumerated, and for each X 1 ∈ 1 , the optimal decision X * 2 of competitor 2 was found using the MIP model. Then, the optimal response X * 1 to X * 2 was found, and if X * 1 = X 1 , the decision profile (X * 1 , X * 2 ) is an NE. Table 1 summarizes the performance of the NESGA and NashGA in terms of discovering NE in 30 random replications. The parameters of the NESGA and NashGA used in the experiments are also given in Table 1. Note that the NashGA was run for a larger g max because the NESGA and NashGA require (2μ) |K| and 2μ|K| objective function evaluations in a replication, respectively. These instances have different numbers of NE, ranging from one to eight (see Tables S1 and S2 in the Supplementary Material for details). Table 1 provides a summary of results as follows: the number of NE (nNE) discovered using enumeration, the number of unique candidate NE (nUCNE) found by the algorithms in 30 replications, the total number of true NE (nTNE) returned in 30 replications, the number of replications in which at least one true NE was discovered (nRTNE) and average central processing unit (CPU) time in seconds. A careful analysis of the results in Table 1 can reveal the pros and cons of the NESGA and NashGA. The results showed that the NESGA could discover multiple NE in a single run. The NESGA was able to discover all NE in the majority of instances with multiple NE. The NashGA also performed well for the instances with |K| = 2, especially if the instances had a limited number of NE, since the NashGA concentrates the search on a single solution. For the instances with |K| = 3, NE was found only for small values of α k with an enumeration. However, such instances provided strong evidence that the NESGA could also discover multiple NE in a single run for |K| = 3. In instance (120,1,1,1,0.25), for example, the NESGA discovered all six symmetrical NE (see Table S2 Table 1. Comparison of the Nash equilibrium sorting genetic algorithm (NESGA) and Nash genetic algorithm (NashGA) on problem instances with known Nash equilibria (NE). NESGA (|K| = 2: µ = 40, g max = 100; |K| = 3: µ = 10, g max = 200) NashGA (|K| = 2: µ = 40, g max = 4000; |K| = 3: µ = 10, g max = 28,000) Problem instance (n,α 1 ,α 2 ,α 3, D) Note: nNE = number of NE discovered using enumeration; nUCNE = number of unique candidate NE; nTNE = total number of true NE returned in 30 replications; nRTNE = number of replications in which at least one true NE was discovered; avg. CPU = average central processing unit time. Table 2. Comparison of the Nash equilibrium sorting genetic algorithm (NESGA) and Nash genetic algorithm (NashGA) on the closest-take-all model problem instances with unknown Nash equilibria (NE).
NESGA (|K| = 2: µ = 40, g max = 200; |K| = 3: µ = 10, g max = 200) NashGA (|K| = 2: µ = 40, g max = 8000; |K| = 3: µ = 10, g max = 28,000) Problem instance (n,α 1 ,α 2, α 3 ,D) in the Supplementary Material) and returned 5.26 true NE per replication (158/30) on average, while the NashGA was not able to find any of these six symmetrical NE. Table 2 presents the results for instances with a larger number of facilities, for which the MIP model was used to verify whether a candidate NE returned by the algorithms was a true NE or not, as well as to calculate the MRE metric. In Table 2, the nVNE column indicates the number of candidate NE that were verified as a true NE, and the other columns have the same meaning as given in Table 1. In addition, the average MRE over 30 replications is provided (the 'Avg. MRE' column). Both algorithms were able to find NE in multiple replications for |K| = 2. However, the NashGA failed to discover an NE for the instances with |K| = 3. Similarly to the results in Table 1, the NESGA returned multiple NE in a single run. Overall, the NESGA identified a slightly larger number of NE in a larger number of replications compared to the NashGA.
In summary, the NESGA was able to discover multiple NE in a single run while the NashGA concentrates the search on a single candidate NE solution. Therefore, the NESGA outperformed the NashGA in instances with |K| = 3 or when the problem instance had multiple NE.

Computational experiments with the Huff gravity-based model
The computational experiments in this section focused on evaluating the performance of the NESGA using the Huff gravity-based model. Table 3 summarizes the results of the problem instances with |K| = 2 for which NE could be determined through enumeration (see Table S3 in the Supplementary Material for the details of the results). The fitness assignment mechanism of the NESGA encourages the population to retain multiple candidate NE, whereas the NashGA intensifies the search around one NE. This advantage of the NashGA became apparent for problem instances with a single NE. For example, in instance (30,2,4,0,0.5) with a single NE, the NashGA returned the NE in all replications, whereas the NESGA did so in 24 replications. On the other hand, the NESGA became advantageous when a problem instance had multiple NE. For example, in instance (30,2,4,0,0.17) with four NE, the NashGA could not find one of them while the NESGA returned this NE in 21 replications [see (30,2,4,0,0.17)-III in Table S2 in the Supplementary Material]. This pattern was observed for the instances of n = 60. In instance (60,2,5,0.18), the NESGA was able to find all nine NE in multiple replications but the NashGA missed some of them. On the other hand, the NashGA performed well for instance (60,2,7,0,0.18), where there is a single NE. In summary, the NESGA was shown to be a more robust algorithm in terms of discovering the whole spectrum of NE. Table 4 presents the results of two small-size problems with |K| = 3. The instance with D = 0.3 has 11 NE, which were identified by exhaustive enumeration. The NESGA found all of these 11 NE in multiple replications. In this instance, the NESGA returned a total of 286 unique candidate NE in 30 replications, and 257 of them were true NE. As seen in Table 4, the NESGA discovered each of the NE in multiple random replications, ranging from 17 to 28. The instance with D = 0.5 has three NE, and the NESGA was also able to find these NE in multiple replications. For this instance, the NESGA returned a total of 84 unique candidate NE in 30 replications, and 57 of them were true NE, whereas the NashGA could not find any of the listed NE. Table 5 summarizes the results of larger problem instances for which at least one true NE was identified. The NESGA was able to identify a true NE in all problem instances, while the NashGA failed to do so in the problems with |K| = 3.

Convergence of the NESGA experiments with the Huff gravity-based model
Since NE may be inferior points in the solution space, demonstrating the convergence of the NESGA requires a different approach than a typical heuristic algorithm does. The MRE metric was used to demonstrate the convergence of the NESGA using instances (60,4,3,2,0.30) and Table 4. Performance of the Nash equilibrium sorting genetic algorithm (NESGA) in simultaneous problems with three competitors (λ = 1, µ = 10, g max = 100 for all).
Problem instance (n,α 1 ,α 2 ,α 3, D) N E π 1 π 2 π 3 Decision profile Note: NE = Nash equilibrium. (60,4,6,0,0.18) from Table 5. Figure 1 illustrates the average and 95% confidence interval (CI) of the MRE observed in 30 random replications of the algorithms for the different levels of the stopping criteria g max = {10,20, . . . ,100}. In both instances, the average MRE was reduced, and the 95% CI became narrower gradually from g max = 10 to g max = 100, indicating that the NESGA converged. The NashGA failed to converge in the case of instance (60,4,3,2,0.30). The failure of the NashGA's convergence in this instance with |K| = 3 is parallel to the results observed in the above tables and can be intuitively explained as follows. To evaluate the fitness of the current populations, the NashGA utilizes the best decision options of the previous populations. In a simultaneous game with |K| = 2, the best response to a given decision of a competitor can always be identified. However, this claim cannot be made in a game with |K| = 3 because there may exist many non-dominated responses to a decision of a competitor by the other two competitors. Hence, the fitness assignment approach of the NashGA does not effectively extend to problems with |K| = 3.

Conclusion
In this article, a new GA called the NESGA was proposed to discover NE in the case of the competitive MCLP with two and three competitors, for the first time. Identifying NE is important to study competitor behaviours that would emerge as a response to market competition, conditions, rules and Figure 1. Illustration of the convergence of the Nash equilibrium sorting genetic algorithm (NESGA) and Nash genetic algorithm (NashGA) using the average (middle dark line) and 95% confidence interval (CI) (grey shaded area) of the mean relative error (MRE) metric.
principles in the context of the competitive MCLP. The NESGA is a multi-population GA with a hierarchical fitness assignment method that ranks population members with respect to all candidate NE in the population. The computational experiments showed that the NESGA could discover multiple NE in a single run, and the proposed fitness evaluation method was effective in problems with three competitors. The NESGA performed robustly over a wide range of problem instances, including ones with three competitors for which the NashGA failed to converge. In addition, a novel way was proposed to demonstrate the convergence of evolutionary algorithms in the case of simultaneous game theory problems. Certainly, the NESGA can be applied to any type of combinatorial game theory problem. Studying the performance of the proposed fitness evaluation method in other combinatorial game theory problems would be an interesting avenue for further research.

Note
1. Throughout the article, the same form of the abbreviation NE is used for both singular 'Nash equilibrium' and plural 'Nash equilibria'.

Disclosure statement
No potential conflict of interest was reported by the authors.

Data availability statement
Data are included in the Supplementary Material.