# A Simulation Approach for Competitive Facility Location with Random Utility
Maximizing Customers

Robin Legault DIRO, Université de Montréal, Montreal, Canada,
<EMAIL_ADDRESS>Emma Frejinger DIRO, Université de Montréal,
Montreal, Canada<EMAIL_ADDRESS>
###### Abstract

The competitive facility location problem formalizes the managerial problem
faced by a firm that enters a market already occupied by competitors in which
customers are assumed to choose the facility that maximizes their utility.
Given a discrete choice model that expresses the probability for a customer to
select each location as a function of its attributes, and given the
distribution of these attributes in the population, the objective is to
identify facility locations that maximize the expected market share captured
by the firm.

We present a simulation approach for solving this problem with any discrete
choice model. Our method exploits the fact that it is possible to aggregate
customers to reduce the dimension of the optimization problem without
affecting its optimal solution. Our experiments indicate that solving the
resulting 0-1 linear program produces high quality solutions for problems
based on the mixed multinomial logit model more efficiently than existing
methods.

In many cases, our approach also leads to near-optimal solutions for large-
scale instances based on the multinomial logit model in a fraction of the
computing time required by the state-of-the-art exact method from the
literature. To interpret these results, we introduce an entropy measure to
characterize the properties that influence the performance of our method on
different types of instances. We finally propose potential uses of this novel
information-theoretic perspective in the broader context of optimization
problems based on random utility maximization models.

Keywords: Competitive facility location, Combinatorial optimization, Random
utility maximization, Simulation, Sample average approximation, Clustering

##  1 Introduction

In the last two decades, a growing body of literature has proposed efficient
exact methods for solving the competitive facility location problem when
customer demand is represented by a multinomial logit (MNL) model (Benati and
Hansen, 2002, Haase, 2009, Aros-Vera et al., 2013, Freire et al., 2016, Ljubić
and Moreno, 2018, Mai and Lodi, 2020, Zhang et al., 2012). The MNL model was
also used to study a variant of the problem that involves congestion Dan and
Marcotte (2019). However, it implies restrictive proportional substitution
patterns. Although empirical studies underscore the need to capture unobserved
taste variations and spatial correlation to provide more realistic modeling of
location-related behavior (Bhat and Guo, 2004, Miyamoto et al., 2004, Müller
et al., 2012), relatively little attention has been devoted to the study of
the competitive facility location problem under more general discrete choice
models. In this work we aim to bridge this gap by proposing an efficient
simulation approach that can be used for any type of random utility
maximization (RUM) discrete choice model.

The mixed multinomial logit (MMNL) model has been considered in some facility
location studies, e.g., Müller et al. (2009) and Haase and Müller (2013).
However, in these works, the MMNL-based problems are solved approximately
using mixed-integer linear programming (MILP) approaches which have been shown
to be dominated by the branch-and-cut method of Ljubić and Moreno (2018) and
the multicut outer-approximation (MOA) algorithm of Mai and Lodi (2020). A
simulation-based approach has been proposed by Haase et al. (2016) to
approximate the MMNL-based competitive facility location problem by a
deterministic MILP. Their method exploits the fact that, for a given sample of
realizations of the model’s random variables, binary coefficients indicating
whether each candidate location would be preferred by a customer to the most
attractive competing location can be computed outside of the optimization
problem. This makes it possible to represent the behavior of each simulated
customer by a single binary variable indicating whether they select a facility
of the decision-maker or a competing location. This type of approach offers a
direct way to integrate flexible discrete choice models in MILP and is now
receiving increasing attention in several fields of application (Paneque et
al., 2021).

In the following, we outline the three main contributions offered by this
paper. We build on this previous work to propose an efficient simulation-based
method for solving the competitive facility location problem under general RUM
models. In existing methods, a set of realizations of the random terms of the
discrete choice model are drawn to produce a sample of simulated customers.
The resulting facility location problem is then directly solved using a MILP
solver. Our first contribution is to add a clustering step to this procedure.
By aggregating the simulated customers whose preference profile is described
by the same coefficient vectors, the dimension of the resulting 0-1 program
can be significantly reduced. We prove that this simple clustering method
affects neither the solution nor the optimal value of the simulated problem
while systematically reducing the computing time of the algorithm.

Our second contribution is to present a detailed computational study in which
the state-of-the-art method by Mai and Lodi (2020) for large maximum capture
facility location problems under the MNL model is compared to our simulation-
based method. Our experiments are performed on two sets of MNL-based benchmark
instances from the literature and on new MMNL-based problems. For real-life
large-scale instances from a park-and-ride location problem in New York City
defined under the MNL model, our method can reduce the resolution time by up
to two orders of magnitude compared to the MOA algorithm while producing near-
optimal solutions (relative optimality gap of less than 0.01% on average).
More importantly, our simulation experiments show that our method requires
several hundred times less computation than the MOA algorithm to identify a
solution of the same quality under the MMNL model. For a given computational
budget, our approach therefore leads to much better solutions. The relative
performance of our method compared to the MOA algorithm can however vary
drastically depending on the parameters of the discrete choice model.

Our third contribution is motivated by this last observation. We introduce a
novel entropy measure that enables us to present an intelligible
interpretation of our computational results. This type of information-
theoretic characterization of instances could be used to quantify the
stochasticity of any optimization problem integrating a random utility-based
discrete choice model. This new perspective can therefore offer valuable
insights allowing to better understand the conditions under which a
simulation-based approach has the potential to improve upon the state-of-the-
art methods for a given application.

The paper is structured as follows. In Section 2, we formulate the maximum
capture problem under a general random utility maximization framework. The MOA
algorithm by Mai and Lodi (2020) is briefly discussed in Section 3. Our
simulation-based approach is then presented in Section 4. In Section 5, we
introduce an information-theoretic perspective on the problem. Computational
experiments are presented in Section 6, and Section 7 concludes the paper.

##  2 Maximum Capture Facility Location Problem with Random Utilities

We consider a probabilistic facility location problem in a competitive market.
Given a set $D$ of available locations on which a firm can install new
facilities and a set $E$ of locations that are already occupied by competing
facilities, the goal is to identify an admissible configuration $x\in X$ that
maximizes the expected market share captured by the new facilities. For each
$d\in D$, the binary variable $x_{d}$ equals one if and only if the company
decides to open a facility at that location. The set of admissible
configurations $X\subseteq\\{0,1\\}^{|D|}$ is specified by linear business
constraints, such as the number of facilities that can be open or the firm’s
budget.

The customers’ behavior is represented by a RUM model. The probability that a
customer chooses a given alternative thus corresponds to the probability that
its utility is the highest among all available alternatives.

The attributes of the customers among the population, such as their location
and preference profile, are modeled by a random variable $\theta$ with support
$\Theta$. A random term $\varepsilon$ also affects the utility of each
alternative. For the decision vector $x\in X$, the open alternatives are given
by the set $C(x)=D(x)\cup E$, where $D(x)=\\{d\in D:x_{d}=1\\}$. The utility
of each available alternative $c\in C(x)$ is defined as
$u_{c}(\theta,\varepsilon)=v_{c}(\theta)+\varepsilon_{c}$, where the random
variables $\theta$ and $\varepsilon$ are assumed to be independent and where
$v_{c}:\Theta\rightarrow\mathbb{R}$ is any real function for each $c\in C(x)$.
Hence, the expected market share captured by a set of locations depends on the
distribution of both $\theta$ and $\varepsilon$.

In most of the literature, it is assumed that a finite number of customers
$n\in N$, each described by a matrix of observed attributes $A_{n}$, form the
population. The relation of customer $n$ to alternative $c$ is specified by a
vector $a_{nc}$ of features (e.g., the distance from the facility, or a level
of appreciation of their service) and by a vector of coefficients $\beta$.
Using the above notation, $v_{c}(\theta)$ can then be expressed as the linear
function $v_{c}(\theta=A_{n})=a_{nc}^{T}\beta$. The distribution of $\theta$
is thus given by the empirical distribution of the matrices $A_{n}$. To use
machine learning terminology, this interpretation underlies a conditional, or
discriminative perspective on the attributes of the customers. In other words,
the observed set of customers is considered fixed instead of being seen as a
random sample drawn from a probability distribution.

For the sake of generality, the notation adopted in this paper instead
proposes a generative perspective, in which all the attributes of the
customers and, eventually, their associated coefficients, are modeled as
(possibly degenerate) random variables. The general maximum capture problem
can therefore be written as the problem of maximizing the expected probability
of capture with respect to the joint distribution of the random vector
$\theta$,

$\displaystyle\max_{\begin{subarray}{c}x\in
X\end{subarray}}\mathbb{E}_{\theta}\Big{[}\mathbb{P}_{\varepsilon}\big{[}\operatorname*{arg\,max}_{c\in
C(x)}\\{u_{c}(\theta,\varepsilon)\\}\in D(x)\big{]}\Big{]}$ (1)
$\displaystyle=$ $\displaystyle\max_{\begin{subarray}{c}x\in
X\end{subarray}}\int_{\theta}\mathbb{P}_{\varepsilon}\big{[}\operatorname*{arg\,max}_{c\in
C(x)}\\{u_{c}(\theta,\varepsilon)\\}\in D(x)\big{]}p(\theta)d\theta$ (2)
$\displaystyle=$ $\displaystyle\max_{\begin{subarray}{c}x\in
X\end{subarray}}\int_{\theta}\sum_{c\in
D(x)}\mathbb{P}_{\varepsilon}\big{[}u_{c}(\theta,\varepsilon)\geq
u_{k}(\theta,\varepsilon),\ \forall k\in C(x)\big{]}p(\theta)d\theta$ (3)
$\displaystyle=$ $\displaystyle\max_{\begin{subarray}{c}x\in
X\end{subarray}}\int_{\theta}\sum_{c\in
D(x)}\mathbb{P}_{\varepsilon}\big{[}\varepsilon_{k}-\varepsilon_{c}\leq
v_{c}(\theta)-v_{k}(\theta),\ \forall k\in C(x)\big{]}p(\theta)d\theta.$ (4)

In order for (4) to be directly embedded in an optimization problem, the
cumulative distribution function (c.d.f.) of the difference
$\varepsilon_{k}-\varepsilon_{c}$ must admit an analytical form. This has
motivated the development of increasingly flexible models belonging to the
generalized extreme value (GEV) class, which all result in closed-form
expressions for the choice probabilities. Prominent members of this family
include the nested logit model (Williams, 1977), the paired combinatorial
logit model (Chu, 1989), and the cross-nested logit model (Vovsha, 1997).

A large part of the RUM literature is based on the assumption that the random
terms $\varepsilon_{c}$ are independent and identically distributed (i.i.d.)
type I extreme value (or Gumbel) variables. In the resulting mixed multinomial
logit (MMNL) model (McFadden and Train, 2000), the conditional probabilities
can be expressed as

$\displaystyle\mathbb{P}_{\varepsilon}\big{[}u_{c}(\theta,\varepsilon)\geq
u_{k}(\theta,\varepsilon),\ \forall k\in C(x)|\
\theta]=\frac{e^{v_{c}(\theta)}}{\sum\limits_{k\in C(x)}e^{v_{k}(\theta)}}.$
(5)

McFadden and Train (2000) proved that, under mild regularity conditions, any
RUM model can be approximated with arbitrary precision by a MMNL model. This
approximation property widely motivates the application of the MMNL model,
although its estimation can be challenging and is still an active research
topic (Bansal et al., 2020).

Under the MMNL model, the maximum capture problem can be written as

$\displaystyle\max_{\begin{subarray}{c}x\in
X\end{subarray}}\int_{\theta}\frac{\sum_{c\in
D(x)}e^{v_{c}(\theta)}}{\sum_{k\in C(x)}e^{v_{k}(\theta)}}p(\theta)d\theta.$
(6)

If $\theta$ has finite support, the problem becomes

$\displaystyle\max_{\begin{subarray}{c}x\in
X\end{subarray}}\sum_{\theta}\frac{\sum_{c\in
D(x)}e^{v_{c}(\theta)}}{\sum_{k\in C(x)}e^{v_{k}(\theta)}}p(\theta).$ (7)

In the latter case, this optimization problem is mathematically equivalent to
the maximum capture problem under the MNL model, in which the utility
functions $v_{c}$ are assumed to be deterministic. It can then be solved
exactly by nonlinear programming (MINLP) solvers, for example from the BONMIN
package (Bonami et al., 2008) or by specialized algorithms, including the
Branch-and-Cut method (BC) presented in Ljubić and Moreno (2018) and the
multicut outer-approximation (MOA) algorithm of Mai and Lodi (2020).

##  3 Exact Resolution Under the MNL Model

Denoting by $N$ the set of $|\Theta|$ possible realizations of
$\theta=(A,\beta)$, whose respective probabilities $q_{n}:=p(\theta),~{}n\in
N$, correspond to the market share represented by customers associating a
utility $v_{nc}:=a_{nc}^{T}\beta$ to each alternative $c\in C(x)$, we obtain
the classical formulation of the problem under the MNL model

$\displaystyle\max_{\begin{subarray}{c}x\in X\end{subarray}}\sum_{n\in
N}q_{n}\frac{\sum_{c\in D(x)}e^{v_{nc}}}{\sum_{k\in C(x)}e^{v_{nk}}}.$ (8)

This model is usually reexpressed as

$\displaystyle\max_{\begin{subarray}{c}x\in X\end{subarray}}\sum_{n\in
N}q_{n}\frac{\sum_{c\in D(x)}e^{v_{nc}}}{\sum_{c\in D(x)}e^{v_{nc}}+W_{nc}},$
(9)

where the term $W_{nc}:=\sum_{c\in E}e^{v_{nc}}$, aggregates the deterministic
utility of all the competing facilities, which do not depend on the decision
variables $x$.

The recently proposed MOA algorithm (Mai and Lodi, 2020) uses an outer-
approximation scheme (Duran and Grossmann, 1986) to solve this problem
exactly. To do so, model (9) is first reformulated as a minimization problem,
with objective function $G(x)=-\sum_{n\in N}q_{n}\frac{\sum_{c\in
D(x)}e^{v_{nc}}}{\sum_{k\in C(x)}e^{v_{nk}}}$. The main idea of the MOA
algorithm is to partition the set $N$ of customers into $T$ subsets, so that
$G(x)$ can be expressed as a sum of $T$ convex and continuously differentiable
functions $g_{t}(x)$. At each iteration of the algorithm, the master problem
$\min_{x\in X}\\{\sum_{t=1}^{T}\phi_{t}|\phi_{t}\geq L_{t},\
\Pi_{t}x-\boldsymbol{1}\phi_{t}\leq\pi_{0t}\ \forall t\\}$ is solved, where
$\Pi_{t}x-\boldsymbol{1}\phi_{t}\leq\pi_{0t}$ is the set of subgradient cuts
corresponding to $g_{t}(x)$ and $L_{t}$ is a lower bound on $g_{t}(x)$. Its
optimal solution $x^{*}$ is then used to add up to $T$ new subgradient cuts
$\phi_{t}\geq\nabla g_{t}(x^{*})(x-x^{*})+g_{t}(x^{*})$, $t=1,\dots,T$ to the
master problem, where $\nabla g_{t}(x^{*})$ is the gradient of $g_{t}$
evaluated at $x^{*}$. This multicut approach generalizes an earlier single-cut
algorithm by Bonami et al. (2008). The single-cut version is obtained by using
only $T=1$ objective function $g_{1}(x)=G(x)$, while taking
$T\in\\{2,\dots,|N|-1\\}$ leads to what is sometimes referred to as a hybrid
approach (Birge and Louveaux, 1988) in the stochastic programming literature.

##  4 Simulation-Based Scheme

Our goal is to identify an efficient method to approximate the optimal
solution of the general model (1). If $\theta$ has finite support and the
random terms $\varepsilon_{c}$ are i.i.d. Gumbel variables, instances of
relatively large size can be efficiently solved to optimality, for example, by
using the MOA algorithm (Mai and Lodi, 2020). However, if the support of
$\theta$ is infinite, the exact methods developed for the MNL model only allow
to solve the maximum capture facility location problem associated with a
finite set of realizations $\\{\psi_{n}\\}_{n\in N}$ sampled from $\theta$.
Solving exactly the resulting problem,

$\displaystyle\max_{\begin{subarray}{c}x\in
X\end{subarray}}\frac{1}{|N|}\sum_{n\in
N}\mathbb{P}_{\varepsilon}\big{[}\operatorname*{arg\,max}_{c\in
C(x)}\\{u_{c}(\psi_{n},\varepsilon)\\}\in D(x)\big{]},$ (10)

remains the best way to approximate the optimal solution of model (1) assuming
that the realizations $\\{\psi_{n}\\}_{n\in N}$ are the only ones available.
It is generally not the case if the sampling distribution of $\theta$ is known
(or estimated with sufficient precision). In this generative perspective, it
could be preferable to solve model (10) approximately if an efficient
heuristic allowed to manage a significantly larger sample for the same
computational budget.

This idea motivates the use of a simulation-based approach in which a set of
realizations $\\{\xi_{ns}\\}_{s\in S}$ is drawn from the distribution of
$\varepsilon$ for each customer $n\in N$. The problem associated with the
empirical distribution $\\{\psi_{n}\\}_{n\in N}$, $\\{\xi_{ns}\\}_{n\in N,s\in
S}$ is given by

$\displaystyle\max_{\begin{subarray}{c}x\in
X\end{subarray}}\frac{1}{|N||S|}\sum_{n\in N}\sum_{s\in
S}\mathds{1}\big{[}\operatorname*{arg\,max}_{c\in
C(x)}\\{u_{c}(\psi_{n},\xi_{ns})\\}\in D(x)\big{]}.$ (11)

As $|S|\rightarrow\infty$, this problem approximates model (10) with an
arbitrarily high precision.

Model (11) can be rewritten as

$\displaystyle\max_{\begin{subarray}{c}x\in
X\end{subarray}}\frac{1}{|N||S|}\sum_{n\in N}\sum_{s\in S}y_{ns}$ (12) s.t.
$\displaystyle y_{ns}\leq\sum_{c\in D}a_{nsc}x_{c},\hskip 14.22636pt$
$\displaystyle\forall n\in N,\forall s\in S$ (13) $\displaystyle
y_{ns}\in\\{0,1\\},$ $\displaystyle\forall n\in N,\forall s\in S$ (14)

where the binary decision variable $y_{ns}$ equals 1 if the simulated customer
$(n,s)\in N\times S$ selects one of the new facilities of the firm. Binary
coefficients $a_{nsc}:=\mathds{1}[u_{c}(\psi_{n},\xi_{ns})\geq
u_{k}(\psi_{n},\xi_{ns})\ \forall k\in E]$ indicate whether the candidate
location $c\in D$ would be preferred by customer $n\in N$ to all the competing
facilities $k\in E$ under scenario $s\in S$. Constraint (13) requires the firm
to open at least one such facility $c\in D$ to capture the demand of simulated
customer $(n,s)$.

Since the coefficients $a_{nsc}$ are computed outside of the optimization
problem, model (12)-(14) is a 0-1 linear programming problem on a set of
$|N|\cdot|S|+|D|$ decision variables. As pointed out by Haase et al. (2016)
(co-authors applied this type of simulation-based approach in the context of
MMNL models in Haase et al., 2018), this formulation corresponds to a
classical deterministic maximum capture facility location problem, such as the
one studied in Revelle (1986).

When a relatively small set of scenarios $S$ is considered, model (12)-(14)
generally requires significantly less computation to be solved to optimality
compared to a naive formulation of model (11). Furthermore, given fixed sets
or realizations $\\{\psi_{n}\\}_{n\in N}$ and $\\{\xi_{ns}\\}_{n\in N,s\in
S}$, the dimension of the simulation-based model can be reduced by
preprocessing that proceeds in two steps:

  1. 1.

The simulated customers are grouped according to their parametric
representation $a_{ns}$. More precisely, two customers
$(n_{1},s_{1}),(n_{2},s_{2})\in N\times S$ are placed in the same cluster $p$
if and only if $a_{n_{1}s_{1}c}=a_{n_{2}s_{2}c}$ $\forall c\in D$.

  2. 2.

If the resulting partition $P$ of $N\times S$ contains a subset $p_{0}$
composed of simulated customers who will be captured by the competitor’s
facilities regardless of the firm’s decisions, then this subset $p_{0}$ is
removed from $P$. More formally, subset $p_{0}$ is composed of simulated
customers $(n,s)\in N\times S$ such that $a_{nsc}=0$ $\forall c\in D$.

This clustering method reduces the number of decision variables as well as the
number of constraints of model (12)-(14) to $|P|$. Let $q_{p}$ denote the
cardinality of each subset $p\in P$. The reduced size model can be written as

$\displaystyle\max_{\begin{subarray}{c}x\in X\end{subarray}}\frac{\sum_{p\in
P}q_{p}y_{p}}{\sum_{p\in P}q_{p}}$ (15) s.t. $\displaystyle
y_{p}\leq\sum_{c\in D}a_{pc}x_{c},\hskip 14.22636pt$ $\displaystyle\forall
p\in P$ (16) $\displaystyle y_{p}\in\\{0,1\\},$ $\displaystyle\forall p\in P$
(17)

where, for each subset $p\in P$, $a_{pc}=a_{nsc}$ $\forall(n,s)\in p$. As we
prove in the following proposition, this simple and inexpensive preprocessing
step can produce a reduction in the resolution time of the simulation-based
problem without affecting its optimal solution.

###### Proposition 4.1

The optimal value and the optimal configuration $x^{*}\in X$ of model
(12)-(14) are identical to those of its reformulation (15)-(17).

##### Proof of Proposition 4.1 :

The set $X$ of feasible configurations is the same for both models. Let $x\in
X$ be a feasible decision vector. Also, let $(n,s)\in N\times S$ be an
arbitrary simulated customer and $p\in P$ the subset to which it has been
assigned. By definition of $P$, $\sum_{c\in D}a_{nsc}x_{c}=\sum_{c\in
D}a_{pc}x_{c}$. Therefore, the decision variable $y_{ns}$ can be set to $1$ in
the first model if and only if $y_{p}$ can be set to $1$ in the second. This
being true for each $p\in P$, and for each of the $q_{p}$ simulated customer
$(n,s)\in p$, then, for any fixed solution $x\in X$, $\max_{y}\\{\sum_{n\in
N}\sum_{s\in S}y_{ns}\text{ s.t.
}(\ref{model:simul_2,con:1}),(\ref{model:simul_2:fin})\\}$ =
$\max_{y}\\{\sum_{p\in P}q_{p}y_{p}\text{ s.t.
}(\ref{model:simul_3,con:1}),(\ref{model:simul_3:fin})\\}$. $\square$

To summarize, in presence of a non-trivial distribution $\theta$ defined on an
infinite support, the principal exact methods of the literature for MNL-based
problems can be used to solve model (10) to optimality if
$\varepsilon_{c}\overset{i.i.d}{\sim}\text{Gumbel}$, but are generally not
appropriate otherwise. If this assumption on the random terms $\varepsilon$ is
respected, drawing a large enough sample of realizations $\\{\psi_{n}\\}_{n\in
N}$ from $\theta$ approximates the optimal solution to (1) with arbitrary
precision. The computational time required for the empirical distribution
$\\{\psi_{n}\\}_{n\in N}$ to adequately approximate the true distribution of
$\theta$ may, however, be prohibitive. Solving (10) to optimality in
reasonable time could therefore force one to use a small sample size which
would be likely to produce a solution that suffers from an important
optimization bias.

The motivation behind the simulation approach is to solve problem (11)
efficiently so that a significantly larger set $\\{\psi_{n}\\}_{n\in N}$ can
be handled with the same computational resources. Although it introduces an
error in the estimation of the conditional probability
$\mathbb{P}_{\varepsilon}\big{[}\operatorname*{arg\,max}_{c\in
C(x)}\\{u_{c}(\psi_{n},\varepsilon)\\}\in D(x)|x;\psi_{n}\big{]}$ associated
with each realization $\psi_{n}$, the simulation approach could improve the
approximation of the optimal solution to model (1) by mitigating the gap
between the empirical distribution of the customers obtained by simulation and
the theoretical distribution of $\theta$.

##  5 Information-theoretic Characterization of Instances

The main assumption underlying the simulation approach is that the stochastic
model describing customers’ choices can be approximated efficiently from a
finite set of realizations of their preferences. In this section, we introduce
an information-theoretic approach for measuring the extent to which this
assumption is respected by a given problem.

In the degenerate case where the behavior of a customer is fully determined by
its attributes (i.e., $u_{c}(\theta,\varepsilon)=v_{c}(\theta)$), model (11)
is exactly equivalent to model (10), since the single possible (degenerate)
scenario $\\{\xi_{n1}\\}_{n\in N}$ properly depicts the behavior of each
customer $n\in N$ almost surely. In the opposite case, where the attributes of
a customer do not provide any information about its behavior (i.e.,
$u_{c}(\theta,\varepsilon)=\varepsilon$), representing each simulated customer
in a deterministic way radically transforms the structure of models (10) and
(11).

Let us define the Multinoulli random variable $Y$, with support $D$, as

$\displaystyle Y:=\operatorname*{arg\,max}_{c\in
D}\\{u_{c}(\theta,\varepsilon)\\}.$ (18)

The two extreme cases mentioned above respectively minimize and maximize the
expectation of the conditional entropy of $Y$ given $\theta$, which we denote
by

$\displaystyle\mathcal{H}($
$\displaystyle\theta,\varepsilon):=\mathbb{E}_{\theta}[H(Y|\theta)]$ (19)
$\displaystyle=\int_{\theta}H(Y|\theta)p(\theta)d(\theta)$ (20)
$\displaystyle=-\int_{\theta}\Big{(}\sum_{c\in
D}\mathbb{P}_{\varepsilon}(Y=c|\theta)\log\mathbb{P}_{\varepsilon}(Y=c|\theta)\Big{)}p(\theta)d(\theta)$
(21)

Indeed, $Y|\theta$ is a constant in the first (degenerate) case and thus
$\mathcal{H}(\theta,\varepsilon)=H(Y|\theta)=0$. In the second case, assuming
that the random terms $\varepsilon_{c}$ are i.i.d., it follows a uniform
distribution over the finite set $D$. Consequently,
$\mathcal{H}(\theta,\varepsilon)=H(Y|\theta)=\log|D|$ reaches the theoretical
upper bound on the entropy of a discrete random variable defined on a finite
support of cardinality $|D|$.

If little uncertainty remains in the behavior of customers after observing
their attributes $\psi_{n}$, the objective function $Z_{NS}(\cdot)$ that is
maximized for a specific set of scenarios $S$ in the simulation-based approach
is likely to share the same maximizer $x^{*}_{N}$ as problem model (10).

Furthermore, considering only a small number $|D|$ of potential new locations
limits the set of possible preference profiles. This can make it possible to
approximate model (10) precisely with relatively few scenarios and to limit
the number of clusters $|P|$ that is necessary to describe exactly the
parametric representation $a_{ns}$ of each simulated customer. The expected
entropy $\mathcal{H}(\theta,\varepsilon)$, which is bounded above by
$\log|D|$, aims to capture the global level of randomness of problem (10). We
can thus expect the simulation approach to perform well when
$\mathcal{H}(\theta,\varepsilon)$ is low. However, it is likely that a
prohibitive number of scenarios will be required to obtain a satisfactory
solution in cases where $\mathcal{H}(\theta,\varepsilon)$ is very high.

In the MNL model, $\theta$ is assumed to be defined on a finite support
corresponding to an observed sample of customers $N$. Under this conditional
perspective, models (1) and (10) are regarded as being equivalent and the
expected entropy $\mathcal{H}(\theta,\varepsilon)$ can be calculated exactly
as

$\mathcal{H}(\theta,\varepsilon)=-\frac{1}{|N|}\sum_{n\in N}\sum_{c\in
D}\frac{e^{v_{nc}}}{\sum\limits_{k\in
D}e^{v_{nk}}}\log\Bigg{(}\frac{e^{v_{nc}}}{\sum\limits_{k\in
D}e^{v_{nk}}}\Bigg{)}.$ (22)

For MMNL-based problems, it can be estimated by computing the average entropy
$H(X|\psi_{n})$ over a large set $\widetilde{N}$ of randomly generated
customers

$\hat{\mathcal{H}}(\theta,\varepsilon)=-\frac{1}{|\widetilde{N}|}\sum_{n\in\widetilde{N}}\sum_{c\in
D}\frac{e^{v_{c}(\psi_{n})}}{\sum\limits_{k\in
D}e^{v_{k}(\psi_{n})}}\log\Bigg{(}\frac{e^{v_{c}(\psi_{n})}}{\sum\limits_{k\in
D}e^{v_{k}(\psi_{n})}}\Bigg{)}.$ (23)

Note that additional measures of randomness would be necessary to characterize
the properties of a model in a perfectly satisfactory way. For example, the
expected entropy of a Bernoulli variable indicating whether the potential
alternative of maximal utility is an element of $D$ could help to determine
how many scenarios should be drawn for the optimal solution of model (11) to
correspond to that of (10) with high probability. In this case the Bernoulli
variable can be expressed as

$\displaystyle\mathbb{E}_{\theta}\Big{[}{H\big{(}\mathds{1}\big{[}\operatorname*{arg\,max}_{c\in
D\cup E}\\{u_{c}(\theta,\varepsilon)\\}\in
D\big{]}\big{|}\theta\big{)}}\Big{]}.$ (24)

If the existing facilities are extremely competitive, a high number of draws
may be necessary before producing a single simulated customer that could be
captured by a new facility, and hence indicate which location the firm should
prioritize. Although the value of $\mathcal{H}(\theta,\varepsilon)$ does not
contain such information, we will see in the next section that this entropy
measure effectively summarizes important features of the problems and enables
us to present an intuitive interpretation of our empirical results.

##  6 Computational Experiments

The purpose of this computational study is twofold. First, we want to compare
the performance of our resolution technique with the state-of-the-art exact
method for MNL-based problems. To assess the relevance of our approach under
the MNL model, our first experiments aim at verifying if the use of simulation
allows to obtain good quality solutions in a computation time significantly
lower than the one required by the MOA algorithm by Mai and Lodi (2020). In
the second part of this section, the goal is to compare both methods under
more general discrete choice models. In order to use the MOA algorithm in this
context, the random terms $\varepsilon_{c}$ are assumed to be i.i.d standard
Gumbel variables, so that model (1) can be reexpressed as (6). Heuristic
solutions to the MMNL-based problem can then be obtained by solving either
problem (10) using the MOA algorithm or problem (11) through formulations
(12)-(14) or (15)-(17).

The experiments are conducted on a PC with processor Intel(R) i7-10875H CPU @
2.30GHz along with 32 GB of RAM operated with Windows 10 Pro. The simulation-
based methods are coded in Julia. We also ported the original MATLAB code of
the MOA algorithm in Julia. All the implemented methods are linked to IBM-ILOG
CPLEX 20.1.0 optimization routines under default settings.

###  6.1 MNL-based Problems

In this section, $\theta$ is assumed to be defined on a finite support
corresponding to an observed sample of customers $N$. Under this conditional
perspective, models (1) and (10) are equivalent. In the following section, we
describe the problem instances and outline the experimental setup. In Section
6.1.2, we delineate the performance measures that we use for reporting the
results (Section 6.1.3).

####  6.1.1 Problem Instances and Experimental Setup

For our experiments, we consider the following classical datasets on which the
most recent exact methods for the MNL-based competitive facility location
problem have been benchmarked (Freire et al., 2016, Ljubić and Moreno, 2018,
Mai and Lodi, 2020). A more detailed description is provided by Freire et al.
(2016).

  * •

HM14: This dataset has been proposed by Haase and Müller (2014). To produce
these problems, the location of the facilities as well as the customers have
been uniformly generated on the set $[0,30]\times[0,30]$. This dataset
includes instances with $|D|\in\\{25,50,100\\}$ available locations and
$|N|\in\\{50,100,200,400\\}$ customers.

  * •

NYC: This problem comes from a large-scale park-and-ride location problem in
New York City. It includes $|D|=59$ available locations and $|N|=82,341$
customers and is generally considered to produce the most challenging
instances of the competitive facility location literature.

For both these sets of problems, the random utility of an available
alternative for customer $n\in N$ is given by $u_{nc}(\varepsilon_{nc})=-\beta
v_{nc}+\varepsilon_{nc}$ for $c\in D$ and by
$u_{nc}(\varepsilon_{nc})=-\alpha\beta v_{nc}+\varepsilon_{nc}$ for the
competing facilities $c\in E$. The perceived utility $v_{nc}$ of candidate
locations $c\in D$ corresponds to the distance between $n$ and $c$ in HM14,
and depends on several factors (including travel time, tolls, auto costs, and
waiting time) in NYC. Thus, parameter $\beta$ controls the importance of the
deterministic term in the total utility and $\alpha$ represents the
competitiveness of the existing facilities (small values of $\alpha$ lead to
attractive competitors).

Like the previous studies, we consider a unique business constraint, namely
that a maximum of $r$ new facilities can be installed by the firm. The
feasible set is thus given by $X=\\{x\in\\{0,1\\}^{|D|}:\sum_{d\in D}\leq
r\\}$. Mai and Lodi (2020) point out that opening a new location cannot reduce
the market share captured by the firm under the MNL model. The inequality
constraint can therefore be replaced by an equality constraint to reduce the
feasible set $X$.

For each of the 12 problems
$(|N|,|D|)\in\\{25,50,100\\}\times\\{50,100,200,400\\}$ of the HM14 dataset,
an instance is solved for each configuration
$(\beta,\alpha,r)\in\\{1,5,10\\}\times\\{0.05,0.1,0.2\\}\times\\{2,3,4,5,6,7,8,9,10\\}$.
These same parameters were used in Ljubić and Moreno (2018) and Mai and Lodi
(2020), except for $\alpha$, which was taken in $\\{0.01,0.1,1\\}$ in both
studies. However, setting $\alpha=1$ produces some degenerate instances in
which more than 99.99% of the market can be captured by a single location. The
parameters we consider produce more reasonable problems for which the optimal
value ranges from 2% to 90%. Similarly, we produced an instance of the NYC
problem for each configuration
$(\beta,\alpha,r)\in\\{0.25,0.5,0.75,1,1.25,1.5,1.75,2\\}\times\\{0.75,1,1.25\\}\times\\{2,3,4,5,6,7,8,9,10\\}$.
The optimal value of the resulting problems lies between 2% and 75%. In total,
these two datasets are respectively composed of 972 and 216 instances.

The simulation-based methods have been applied with $|S|\in\\{1,10,100\\}$
scenarios on the HM14 dataset and with $|S|\in\\{1,2,5\\}$ scenarios on the
NYC dataset. The MOA algorithm has been executed with $T=|N|$ groups for the
HM14 instances since our preliminary experiments have shown that the multicut
version was more efficient than the single-cut and hybrid versions of the
algorithm on this dataset. For the NYC instances, we choose $T=20$, as
prescribed by Mai and Lodi (2020).

The instances were grouped according to the parameters having a direct
influence on the average entropy $\mathcal{H}(\theta,\varepsilon)$, namely
$\beta$ and $|D|$ for HM14 and $\beta$ for NYC. It is worth noticing that, for
both groups of problems, given a set $D$ of available locations,
$\mathcal{H}(\theta,\varepsilon)$ varies only superficially between instances,
i.e., from one set of customers $N$ to another. This follows from the fact
that $\mathcal{H}(\theta,\varepsilon)$ corresponds to the estimator of the
expected entropy $\mathcal{H}(\theta^{*},\varepsilon)$ given in (22), for
$\widetilde{N}=N$, where $\theta^{*}$ is the random variable associated with
the (potentially unknown) generative process that produced the customers of a
dataset.

####  6.1.2 Performance Measures

For each group of instances, we present the computing time (second set of
columns in Tables 1 and 2) required by the MOA algorithm and by the
simulation-based method without clustering (SB) and with clustering (SBC). In
the third set of columns, we report the relative optimality gap. For a
feasible solution $x\in X$, it is defined as

$\displaystyle\textit{RGap}(x;N)=100\cdot\frac{Z_{N}(x^{*}_{N})-Z_{N}(x)}{Z_{N}(x^{*}_{N})},$
(25)

where $Z_{N}(\cdot)$ and $x^{*}_{N}$ respectively denote the objective
function and the optimal solution of problem (10). The relative optimality gap
of the solution obtained by SB and SBC can be expressed as
$\textit{RGap}(x^{*}_{NS};N)$, where $x^{*}_{NS}$ denotes the optimal solution
to problem (11) for a set of realizations $\\{\xi_{ns}\\}_{n\in N,s\in S}$ of
the random terms $\varepsilon$, and where the set $\\{\psi_{n}\\}_{n\in N}$
corresponds to the observed sample of customers $N$ on which the MNL-based
problem is defined.

We also report the reduction in the number of decision variables $y$ that
results from applying our clustering technique to the simulation-based problem
(12)-(14) to obtain reformulation (15)-(17). This reduction, in percentage, is
given by $100\cdot(1-\frac{|P|}{|N||S|})$.

####  6.1.3 Results

The numerical results for HM14 and NYC are presented in Tables 1 and 2,
respectively. They suggest that the simulation-based approach, for an
appropriate choice of $|S|$, offers a clear computational advantage over the
MOA algorithm and only causes a slight loss of accuracy in the resolution of
problem (11) for most types of instances.

On the HM14 dataset, taking $|S|=100$ scenarios leads to an average relative
optimality gap of less than 0.4% for each group of instances while dividing
the computation time by approximately ten.

The quality of the solution is even better for the NYC problems. Even with
only $|S|=1$ scenario per customer, the relative optimality gap remains as low
as 0.01% on average. This is understandable because of the high number of
customers in these instances. Since $|N|=82,341$, many customers share very
similar sets of attributes. Taken together, the unique scenarios sampled for
each of these customers are thus likely to faithfully depict the preference
profile of the group. This causes the optimal solutions to problems (10) and
(11) to be identical for practically all the instances of this dataset.

Table 1: Average numerical results for HM14 dataset, grouped by $\beta$ and
$D$ (108 instances per row)

Parameters and entropy | Computing time (s) |  RGap (%) | Size reduction (%)
---|---|---|---
SB, $|S|$ | SBC, $|S|$ |  MOA | SB/SBC, $|S|$ | SBC, $|S|$
$\beta$ | $|D|$ | $\mathcal{H}(\theta,\varepsilon)$ | 1 | 10 | 100 | 1 | 10 | 100 | 1 | 10 | 100 | 1 | 10 | 100
10 | 25 | 0.09 | 0.01 | 0.02 | 0.11 | 0.01 | 0.01 | 0.01 | 0.05 | 3.06 | 0.17 | 0.03 | 81.21 | 97.31 | 99.67
10 | 50 | 0.11 | 0.01 | 0.02 | 0.13 | 0.01 | 0.01 | 0.03 | 0.10 | 3.78 | 0.21 | 0.02 | 66.02 | 94.71 | 99.22
10 | 100 | 0.17 | 0.01 | 0.03 | 0.16 | 0.01 | 0.02 | 0.04 | 0.13 | 3.10 | 0.48 | 0.01 | 45.19 | 88.01 | 97.95
5 | 25 | 0.17 | 0.01 | 0.02 | 0.11 | 0.01 | 0.01 | 0.01 | 0.09 | 5.39 | 0.48 | 0.06 | 80.88 | 96.64 | 99.54
5 | 50 | 0.25 | 0.01 | 0.02 | 0.13 | 0.01 | 0.01 | 0.03 | 0.16 | 8.77 | 0.51 | 0.03 | 64.98 | 92.41 | 98.76
5 | 100 | 0.36 | 0.01 | 0.03 | 0.18 | 0.01 | 0.02 | 0.06 | 0.37 | 6.44 | 0.64 | 0.06 | 44.17 | 82.25 | 95.91
1 | 25 | 0.78 | 0.01 | 0.02 | 0.13 | 0.01 | 0.01 | 0.03 | 0.53 | 8.53 | 1.60 | 0.38 | 72.33 | 90.55 | 97.55
1 | 50 | 1.22 | 0.01 | 0.04 | 0.23 | 0.01 | 0.03 | 0.13 | 1.77 | 11.79 | 3.53 | 0.13 | 49.50 | 74.34 | 90.08
1 | 100 | 1.68 | 0.01 | 0.12 | 2.46 | 0.01 | 0.11 | 2.24 | 31.14 | 12.71 | 4.11 | 0.29 | 29.23 | 45.67 | 65.69

Table 2: Average numerical results for NYC dataset, grouped by $\beta$ (27
instances per row)

Parameters and entropy | Computing time (s) |  RGap (%) | Size reduction (%)
---|---|---|---
SB, $|S|$ | SBC, $|S|$ |  MOA | SB/SBC, $|S|$ | SBC, $|S|$
$\beta$ | $\mathcal{H}(\theta,\varepsilon)$ | 1 | 2 | 5 | 1 | 2 | 5 | 1 | 2 | 5 | 1 | 2 | 5
2.00 | 0.19 | 0.71 | 1.36 | 3.86 | 0.14 | 0.31 | 0.81 | 16.94 | 0.01 | 0.01 | 0.00 | 96.77 | 97.77 | 98.67
1.75 | 0.22 | 0.70 | 1.38 | 3.87 | 0.15 | 0.34 | 0.82 | 13.64 | 0.01 | 0.00 | 0.00 | 96.62 | 97.64 | 98.55
1.50 | 0.26 | 0.71 | 1.40 | 3.93 | 0.16 | 0.36 | 0.85 | 9.37 | 0.00 | 0.00 | 0.00 | 96.40 | 97.43 | 98.39
1.25 | 0.31 | 0.70 | 1.41 | 3.95 | 0.17 | 0.36 | 0.92 | 6.41 | 0.00 | 0.00 | 0.00 | 96.10 | 97.13 | 98.16
1.00 | 0.40 | 0.72 | 1.47 | 4.07 | 0.17 | 0.44 | 0.98 | 3.41 | 0.01 | 0.01 | 0.00 | 95.56 | 96.64 | 97.76
0.75 | 0.54 | 0.83 | 1.51 | 4.16 | 0.27 | 0.48 | 1.23 | 2.14 | 0.01 | 0.00 | 0.00 | 94.53 | 95.72 | 96.98
0.50 | 0.82 | 0.87 | 1.76 | 5.06 | 0.40 | 0.69 | 2.04 | 1.24 | 0.04 | 0.00 | 0.01 | 92.25 | 93.59 | 95.12
0.25 | 1.62 | 1.71 | 3.37 | 9.19 | 1.14 | 2.16 | 6.08 | 0.64 | 0.00 | 0.00 | 0.00 | 84.67 | 86.37 | 88.42

As expected, the average entropy $\mathcal{H}(\theta,\varepsilon)$ has a
direct impact on the overall performance of the simulation-based methods. As
the level of randomness of model (10) decreases, model (11) becomes easier to
solve for any number of scenarios and SBC becomes significantly more efficient
than SB. While the clustering step reduces the computing time of the
simulation-based method by approximately 34% on the NYC instances when
$\beta=0.25$, this reduction increases to 79% for $\beta=2$. Overall, SBC
requires about seven times more computations on the highest entropy instances
compared to the lowest entropy problems.

Interestingly, the entropy has a completely different effect on the MOA
algorithm. The HM14 instances seem to become increasingly difficult for all
approaches as $\mathcal{H}(\theta,\varepsilon)$ becomes larger. However, the
computation times of MOA and SBC react to entropy in perfectly opposite ways
for the NYC dataset. This observation can possibly be explained by the fact
that the MOA uses a small number of cuts ($T=20$) in this setting. Since
customers are simply clustered into groups of the same size, the single-cut
and hybrid versions of the MOA may lead to stronger cuts when the entropy is
high. Indeed, the choice probabilities of all the customers become more
similar as $\beta$ increases, which may reduce the benefit associated with
adding cuts based on each concave component of the objective function at each
iteration. Using this aggregation method may therefore become more sound as
the sensitivity of customers about the perceived utilities, and thus the
heterogeneity of their behavior, decreases.

To summarize, the results presented in this section demonstrate that our
simulation approach can be of practical interest for solving large-scale MNL-
based instances of the competitive facility location problem in situations
where using an exact method is too computationally demanding. Furthermore, we
have shown that the entropy measure introduced in Section 5 provides a useful
characterization of the instances. In particular, it allowed us to illustrate
the complementarity of the simulation-based and MOA methods for the real-life
large-scale instances from the NYC dataset.

###  6.2 MMNL-based Problems

In this section, we adopt a generative perspective to study the resolution of
model (1) when $\theta$ is a random variable with infinite support. The
section has the same structure as the previous one on MNL-based problems. That
is, we start by describing the problem instances and outlining the
experimental setup. We then delineate the performance measures (Section 6.2.2)
followed by the results (Section 6.2.3).

####  6.2.1 Problem Instances and Experimental Setup

Our experiments are carried out on two sets of problems, which we denote by
HM14-MMNL and MMNL-3.

In the HM14-MMNL dataset, instead of adopting the conditional perspective
which leads to the HM14 MNL-based instances, we directly consider the
generative model presented by Haase and Müller (2014). A set $D$ of $50$
candidate locations as well as $|E|=10$ competing facilities are first
uniformly generated on the rectangular region $[0,30]\times[0,30]$. These
locations are then kept fixed throughout the experiments. The spatial
distribution of the customers is given by the random variable
$\theta=(\theta_{1},\theta_{2})$, where
$\theta_{1},\theta_{2}\overset{i.i.d}{\sim}\text{Unif}(0,30)$. The component
$v_{c}(\theta)$ of the utility functions is defined as the Manhattan distance
between $\theta$ and location $c$ for each open location $c\in D\cup E$. The
expected entropy $\mathcal{H}(\theta,\varepsilon)$ of the resulting problem is
approximately 0.96, which corresponds to a slightly higher level of randomness
compared to the average entropy of the classical HM14 instances.

In the new group of problems we denote by MMNL-3, three types of locations
$l\in L$ (we could think of commercial spaces on a shopping street, in a mall
or in a rest area, for example) as well as three types of customers
$k\in\mathcal{K}$ (e.g., teenagers, adults and seniors) are considered. We
suppose that a fixed set of competing facilities and available locations, each
belonging to a specific type of locations, is given. In practice, we generated
the locations randomly on the region $[-15,15]\times[-15,15]$. The
distribution of the customers’ attributes among the population is given by the
random variable $\theta=(\theta_{1},\theta_{2},K)$. The pair
$(\theta_{1},\theta_{2})$ denotes a position in the plane and the categorical
random variable $K$ is the type of the customer. The coefficients
$\\{\delta_{k}\\}_{k\in\mathcal{K}}$ and
$\\{\gamma_{kl}\\}_{k\in\mathcal{K},l\in L}$ respectively control the level of
aversion of each customer type to travel distances, and to each type of
location. A scaling parameter $\beta$ controls the weight of function $v_{c}$
in the total utility
$u_{c}(\theta,\varepsilon)=v_{c}(\theta)+\varepsilon_{c}$. To summarize, for a
facility $c$ of type $l$, $v_{c}(\theta)=-\beta(\delta_{K}M_{c}+\gamma_{Kl})$,
where $M_{c}$ denotes the Manhattan distance separating position
$(\theta_{1},\theta_{2})$ from location $c$. The following parameters are kept
unchanged for all the computational experiments:

$\begin{cases}\delta&=[3,\ 1,\ 2],\\\ \gamma&=\begin{bmatrix}20&60&30\\\
40&20&60\\\ 60&40&20\end{bmatrix}.\end{cases}$

To define the distribution of $\theta$ in this model, we suppose that the
population is divided into four neighborhoods $j\in\\{1,2,3,4\\}$. The part of
the customers residing in neighborhood $j$ is given by $\pi_{j}$, and
$\rho_{jk}$ specifies the proportion of customers of type $k$ in neighborhood
$j$. The spatial distribution of the population in neighborhood $j$ is given
by a bivariate normal with mean $\mu_{j}$ and covariance matrix $\Sigma_{j}$.
The following parameters are used in our experiments:

$\begin{cases}\pi&=[0.4,\ 0.3,\ 0.2,\ 0.1],\\\
\rho&=\begin{bmatrix}0.2&0.7&0.1\\\ 0.3&0.4&0.3\\\ 0.3&0.4&0.3\\\
0.0&0.2&0.8\\\ \end{bmatrix},\ \mu=\begin{bmatrix}\begin{pmatrix}2\\\ -2\\\
\end{pmatrix},\begin{pmatrix}-10\\\ -10\\\ \end{pmatrix},\begin{pmatrix}-4\\\
10\\\ \end{pmatrix},\begin{pmatrix}12\\\ -5\\\ \end{pmatrix}\end{bmatrix},\\\
\Sigma&=\begin{bmatrix}\begin{bmatrix}9&1\\\ 1&9\\\
\end{bmatrix},\begin{bmatrix}9&-6\\\ -6&9\\\
\end{bmatrix},\begin{bmatrix}16&1\\\ 1&4\\\
\end{bmatrix},\begin{bmatrix}2&0\\\ 0&21\\\
\end{bmatrix}\end{bmatrix}.\end{cases}$

As in the previous section, we only consider a constraint on the maximum
number of new facilities that can be installed by the firm. This bound is set
to $r=5$ for the HM14-MMNL problems and to $r=10$ for MMNL-3.

For the HM14-MMNL instances, 50 instances are generated for each sample size
$|N|\in\\{125,250,500,1000,2000\\}$. Problems (10) and (11) are defined on the
same set of simulated customers $N$ and a set of $|S|\in\\{1,10,100\\}$
scenarios are considered for the simulation approach.

Figure 1: Visualization of the distribution of $\theta$ for the MMNL-3 problem
on a sample of $|N|=10,000$ realizations $\\{\psi_{n}\\}_{n\in N}$.

We then consider problem from the MMNL-3 dataset with 10 competing facilities
and 20 open alternatives of each type $l\in L$, for a total of $|D|=60$
available locations and $|E|=30$ competitors. Figure 1 presents the location
and the type of each facility and illustrates the distribution of $\theta$ on
a sample of $10,000$ customers. The MOA algorithm is used to solve model (10)
for $|N|\in\\{50,100,200,400,800,1600,3200\\}$. As our primary objective is
not to solve problem (10) exactly for a finite random sample of customers, but
instead to efficiently find a near-optimal solution to problem (1), a single
scenario is considered for each simulated customer $n\in N$ for this dataset.
Since problem (11) is significantly less challenging than problem (10) when
$|S|=1$, larger sample sizes
$|N|\in\\{100,400,1600,6400,25600,102400,409600\\}$ are generated for the
simulation method. The values tested for the sensitivity parameter of the
customers to perceived utility are $\beta\in\\{0.25,0.5,1\\}$.

The expected entropy $\mathcal{H}(\theta,\varepsilon)$ of the resulting
problems is approximately 1.59, 0.90 and 0.46 respectively. As for the
HM14-MMNL instances, 50 samples are generated for each configuration. The
numerical results are presented in Table 4. The multi-cut version is the most
efficient variant of the MOA algorithm for both sets of MMNL-based problems.
Hence, all the reported results in this section have been obtained with
$T=|N|$.

####  6.2.2 Performance Measures

Under the generative perspective we adopt in this section, the performance of
a method is measured by the quality of the best feasible solution to the
stochastic problem (1) it can produce for a given computational budget. Let
$x\in X$ be a feasible configuration, its value $Z(x)$ can be estimated by

$\displaystyle\hat{Z}(x)=Z_{\widetilde{N}}(x),$ (26)

where $\widetilde{N}$ is a very large sample of realizations
$\\{\psi_{n}\\}_{n\in\widetilde{N}}$ drawn from the distribution of $\theta$.
In our experiments, we use $|\widetilde{N}|=1,000,000$. This same set of
simulated customers is used to compute the estimator
$\hat{\mathcal{H}}(\theta,\varepsilon)$ of the expected entropy of $Y|\theta$.

In this context, given a sample $N$ of simulated customers whose attributes
correspond to realizations $\\{\psi_{n}\\}_{n\in N}$ of the random variable
$\theta$, the MOA algorithm is used to compute the optimal solution
$x^{*}_{N}$ to model (10), with objective value $Z_{N}(x^{*}_{N})$. For the
simulation-based approach, in addition to a sample of customers $N$, a set of
realizations $\\{\xi_{ns}\\}_{s\in S}$ is drawn from the distribution of
$\varepsilon$ for each of them. The optimal solution of the resulting problem
(11) is denoted by $x^{*}_{NS}$ and its value is given by
$Z_{NS}(x^{*}_{NS})$. We will therefore be interested in comparing the
estimated values $\hat{Z}(x^{*}_{N})$ and $\hat{Z}(x^{*}_{NS})$ for different
samples sizes $|N|$ and $|S|$.

Sample average approximation estimates are known to be biased. On average,
they provide optimal values that are not attainable for the original
stochastic problem (Norkin et al., 1998, Shapiro, 2013). In our case, denoting
by $x^{*}$ an optimal solution to model (1), this means that, for any sample
size $|N|$ and any number of scenarios $|S|$,

$\displaystyle\mathbb{E}[Z_{NS}(x^{*}_{NS})]>\mathbb{E}[Z_{N}(x^{*}_{N})]>Z(x^{*}).$
(27)

We can therefore consider that solving models (10) for a sample of simulated
customers of size $|N|$ or model (11) for a configuration of sample sizes
$(|N|,|S|)$ consistently provides a quasi-optimal solution to problem (1) if
the following relative generalization gap is very close to zero on average.
For model (10), we define this measure as

$\displaystyle{\it
RGenGap}_{N}(x^{*}_{N})=100\cdot\frac{Z_{N}(x^{*}_{N})-\hat{Z}(x^{*}_{N})}{Z_{N}(x^{*}_{N})}.$
(28)

Similarly, it is given by

$\displaystyle{\it
RGenGap}_{NS}(x^{*}_{NS})=100\cdot\frac{Z_{NS}(x^{*}_{NS})-\hat{Z}(x^{*}_{NS})}{Z_{NS}(x^{*}_{NS})}.$
(29)

for model (11).

####  6.2.3 Results

The purpose of the experiments carried out on the HM14-MMNL dataset is to
study the respective impact of $|N|$ and $|S|$ on the quality of the solutions
to problem (1) provided by models (10) and (11). The main conclusion to be
drawn from the results in Table 3 is that the large computational effort
required to solve model (10) exactly only produces a marginal improvement in
the quality of the solution to problem (1). Faithfully representing the
distribution of $\theta$ through a large sample $N$ of simulated customers
appears to be more critical than approximating the behavior of these customers
with high fidelity. In particular, we notice that the MOA algorithm does not
allow to close the relative generalization gap in reasonable time, despite the
relatively small size of this problem. Indeed, when $|N|=2000$, solutions
$x^{*}_{N}$ lead to a value approximately $2.39\%$ larger for problem (10)
than for problem (1) and require more than one minute to be computed on
average. Solving model (11) with the SBC method for $|S|=100$ scenarios
produces solutions of almost identical quality nearly 100 times faster.

Table 3: Average numerical results for HM14-MMNL model (50 randomly generated
instances per row)

| Computing time (s) | $Z_{N}(x)$ | $\hat{Z}(x)$
---|---|---|---
| SB, $|S|$ | SBC, $|S|$ |  MOA |  $x=x^{*}_{NS}$, $|S|$ |  $x=x^{*}_{N}$ |  $x=x^{*}_{NS}$, $|S|$ |  $x=x^{*}_{N}$
$N$ | 1 | 10 | 100 | 1 | 10 | 100 | 1 | 10 | 100 | 1 | 10 | 100
125 | 0.02 | 0.04 | 0.20 | 0.02 | 0.03 | 0.12 | 0.59 | 49.08 | 51.08 | 51.49 | 51.53 | 43.33 | 44.13 | 44.24 | 44.30
250 | 0.02 | 0.06 | 0.39 | 0.02 | 0.05 | 0.22 | 1.79 | 48.66 | 49.60 | 49.79 | 49.81 | 44.37 | 44.98 | 44.99 | 44.91
500 | 0.03 | 0.11 | 0.91 | 0.02 | 0.09 | 0.57 | 5.70 | 47.60 | 48.38 | 48.49 | 48.50 | 45.19 | 45.47 | 45.56 | 45.58
1000 | 0.05 | 0.22 | 1.21 | 0.04 | 0.17 | 0.49 | 23.44 | 47.13 | 47.57 | 47.65 | 47.66 | 45.48 | 45.73 | 45.75 | 45.80
2000 | 0.07 | 0.37 | 2.63 | 0.06 | 0.32 | 0.76 | 66.95 | 46.78 | 47.14 | 47.19 | 47.19 | 45.82 | 46.05 | 46.05 | 46.06

Table 4: Average numerical results for MMNL-3 dataset (50 randomly generated
instances per row)

Parameter and entropy | MOA | SB/SBC, $|S|=1$
---|---|---
$|N|$ |  Time (s) |  $Z_{N}(x^{*}_{N})$ |  $\hat{Z}(x^{*}_{N})$ |  RGen $\textit{Gap}_{N}$ ($\%$) |  $|N|$ | Time (s) |  $Z_{NS}(x^{*}_{NS})$ |  $\hat{Z}(x^{*}_{NS})$ |  RGen ${\it Gap}_{NS}$ ($\%$)
$\beta$ | $\hat{\mathcal{H}}(\theta,\varepsilon)$ | SB | SBC
1 |  0.46 | 50 | 0.13 | 55.51 | 41.66 | 24.96 | 100 | 0.01 | 0.01 | 52.28 | 42.16 | 19.35
100 | 0.15 | 51.21 | 43.71 | 14.63 | 400 | 0.02 | 0.01 | 49.62 | 45.20 | 8.92
200 | 0.24 | 49.05 | 44.68 | 8.90 | 1600 | 0.02 | 0.02 | 47.55 | 46.25 | 2.74
400 | 0.43 | 47.61 | 45.59 | 4.23 | 6400 | 0.06 | 0.02 | 46.99 | 46.54 | 0.96
800 | 1.12 | 47.40 | 46.00 | 2.93 | 25600 | 0.18 | 0.04 | 46.78 | 46.65 | 0.30
1600 | 2.73 | 47.48 | 46.43 | 2.21 | 102400 | 0.75 | 0.15 | 46.79 | 46.67 | 0.25
|  | 3200 | 9.38 | 47.11 | 46.54 | 1.21 | 409600 | 3.54 | 0.51 | 46.75 | 46.68 | 0.15
0.5 |  0.90 | 50 | 0.18 | 52.49 | 40.54 | 22.77 | 100 | 0.01 | 0.01 | 51.06 | 39.87 | 21.92
100 | 0.28 | 48.68 | 41.79 | 14.15 | 400 | 0.02 | 0.02 | 47.99 | 42.87 | 10.67
200 | 0.78 | 46.85 | 42.84 | 8.55 | 1600 | 0.03 | 0.02 | 45.65 | 43.75 | 4.15
400 | 1.60 | 45.41 | 43.38 | 4.47 | 6400 | 0.09 | 0.06 | 44.59 | 44.01 | 1.29
800 | 6.26 | 45.09 | 43.69 | 3.11 | 25600 | 0.25 | 0.11 | 44.46 | 44.20 | 0.58
1600 | 19.55 | 45.03 | 44.02 | 2.24 | 102400 | 1.00 | 0.32 | 44.39 | 44.23 | 0.35
|  | 3200 | 63.79 | 44.70 | 44.09 | 1.37 | 409600 | 3.98 | 0.80 | 44.29 | 44.26 | 0.07
0.25 |  1.59 | 50 | 0.27 | 47.76 | 38.34 | 19.71 | 100 | 0.01 | 0.01 | 49.16 | 36.75 | 25.25
100 | 0.50 | 44.73 | 39.51 | 11.67 | 400 | 0.02 | 0.02 | 45.75 | 39.25 | 14.22
200 | 1.66 | 43.28 | 40.08 | 7.39 | 1600 | 0.04 | 0.03 | 42.79 | 40.57 | 5.19
400 | 4.60 | 42.10 | 40.62 | 3.52 | 6400 | 0.16 | 0.13 | 41.61 | 40.94 | 1.60
800 | 18.26 | 41.82 | 40.81 | 2.42 | 25600 | 0.96 | 0.82 | 41.47 | 41.11 | 0.89
1600 | 48.42 | 41.74 | 40.99 | 1.80 | 102400 | 2.61 | 1.93 | 41.36 | 41.15 | 0.49
|  | 3200 | 179.97 | 41.55 | 41.04 | 1.22 | 409600 | 11.9 | 7.80 | 41.19 | 41.16 | 0.09

Table 4 presents the numerical results for MMNL-3. These results clearly
demonstrate the value of an effective simulation-based method to solve the
competitive facility location problem under the MMNL model. For the same
number of simulated customers, solving model (10) provides better solutions
than SB and SBC for problem (1), since only one scenario was generated for
each customer $(|S|=1)$. However, the resulting computational cost can be
several orders of magnitude higher. When the entropy is relatively low
($\mathcal{H}(\theta,\varepsilon)<1$ in this example), SBC makes it possible
to solve model (11) for several hundred thousand simulated customers in less
than one second, which enables us to practically close the relative
generalization gap. On the contrary, problem (10) rapidly becomes challenging
when more than a few thousand simulated customers are considered. As a result,
the MOA algorithm did not achieve an average relative generalization gap below
1% for any of the above problems.

In conclusion, the computational experiments presented in this section show
that the strategy of using the MOA algorithm to solve MMNL-based problems is
dominated by our simulation approach.

##  7 Conclusion

In this paper, we have proposed a simulation approach for solving the
competitive facility location problem with random utility maximizing customers
that can be used for any type of RUM discrete choice model. It is based on a
new reformulation in which simulated customers sharing the same preference
profile are aggregated. This clustering step systematically reduces the
resolution time without impacting the quality of the solution.

We compared our approach with the state-of-the-art exact method for MNL-based
problems on a large set of benchmark instances from the literature and on new
MMNL-based problems. Our experiments showed that our resolution method
drastically outperforms the existing approaches when the customers’ behavior
is described by a MMNL model. It can also provide near-optimal solutions for
large-scale problems under the MNL model with a significantly reduced
computation budget compared to the MOA algorithm by Mai and Lodi (2020) when
observable attributes are strong predictors of a customer’s behavior.

To formalize this idea, we proposed an entropy measure to quantify the level
of randomness of the customers’ behavior. This information-theoretic approach
could be adapted to analyze the structure of other optimization problems based
on random discrete choice models. Consequently, it constitutes a new tool to
better understand the conditions under which a simulation-based method has the
potential to improve the state-of-the-art algorithms for a given application.

In the case of stochastic facility location problems, our entropy measure
offers valuable information that could be exploited in the development of
adaptive algorithms. This line of research could be worth exploring, since our
results indicate that our simulation approach and the MOA algorithm can react
oppositely to the entropy of a problem. One could take advantage of this
complementarity to design algorithms that are more robust than the existing
methods.

Among the new research questions raised by this work, it could also be of
particular interest to use our simulation approach as a heuristic step within
an exact method for the MNL model. For example, identifying a feasible
solution of good quality by simulation could allow the MOA algorithm to
rapidly generate strong subgradient cuts and hopefully, to converge after a
significantly reduced number of iterations. The development of alternative
clustering methods is another promising topic for future work.

## Acknowledgments

This research was supported by FRQNT and IVADO through scholarships to the
first author. We wish to thank Tien Mai and Andrea Lodi for making available
the source code of the MOA algorithm as well as their benchmark instances.

## References

  * Aros-Vera et al. (2013) Aros-Vera, F., Marianov, V., and Mitchell, J. E.  p-Hub approach for the optimal park-and-ride facility location problem.  _European Journal of Operational Research_ , 226(2):277–285, 2013.
  * Bansal et al. (2020) Bansal, P., Krueger, R., Bierlaire, M., Daziano, R. A., and Rashidi, T. H.  Bayesian estimation of mixed multinomial logit models: Advances and simulation-based evaluations.  _Transportation Research Part B: Methodological_ , 131:124–142, 2020.
  * Benati and Hansen (2002) Benati, S. and Hansen, P.  The maximum capture problem with random utilities: Problem formulation and algorithms.  _European Journal of Operational Research_ , 143(3):518–530, 2002.
  * Bhat and Guo (2004) Bhat, C. R. and Guo, J.  A mixed spatially correlated logit model: formulation and application to residential choice modeling.  _Transportation Research Part B: Methodological_ , 38(2):147–168, 2004.
  * Birge and Louveaux (1988) Birge, J. R. and Louveaux, F. V.  A multicut algorithm for two-stage stochastic linear programs.  _European Journal of Operational Research_ , 34(3):384–392, 1988.
  * Bonami et al. (2008) Bonami, P., Biegler, L. T., Conn, A. R., Cornuéjols, G., Grossmann, I. E., Laird, C. D., Lee, J., Lodi, A., Margot, F., Sawaya, N., and Wächter, A.  An algorithmic framework for convex mixed integer nonlinear programs.  _Discrete Optimization_ , 5(2):186–204, 2008\.
  * Chu (1989) Chu, C.  A paired combinatorial logit model for travel demand analysis.  In _Proceedings of the Fifth World Conference on Transportation Research, 1989_ , volume 4, pages 295–309, 1989.
  * Dan and Marcotte (2019) Dan, T. and Marcotte, P.  Competitive facility location with selfish users and queues.  _Operations Research_ , 67(2):479–497, 2019.
  * Duran and Grossmann (1986) Duran, M. A. and Grossmann, I. E.  An outer-approximation algorithm for a class of mixed-integer nonlinear programs.  _Mathematical Programming_ , 36(3):307–339, 1986\.
  * Freire et al. (2016) Freire, A. S., Moreno, E., and Yushimito, W. F.  A branch-and-bound algorithm for the maximum capture problem with random utilities.  _European Journal of Operational Research_ , 252(1):204–212, 2016.
  * Haase et al. (2016) Haase, K., Müller, S., Krohn, R., and Hensher, D.  The maximum capture problem with flexible substitution patterns. Hamburg University Working paper.  2016\.
  * Haase (2009) Haase, K.  Discrete location planning.  2009\.  URL http://hdl.handle.net/2123/19420.
  * Haase and Müller (2013) Haase, K. and Müller, S.  Management of school locations allowing for free school choice.  _Omega_ , 41(5):847–855, 2013.
  * Haase and Müller (2014) Haase, K. and Müller, S.  A comparison of linear reformulations for multinomial logit choice probabilities in facility location models.  _European Journal of Operational Research_ , 232(3):689–691, 2014.
  * Haase et al. (2018) Haase, K., Müller, S., Krohn, R., and Hensher, D.  School location planning with flexible substitution patterns.  doi:10.13140/RG.2.2.28100.94084, 2018.
  * Ljubić and Moreno (2018) Ljubić, I. and Moreno, E.  Outer approximation and submodular cuts for maximum capture facility location problems with random utilities.  _European Journal of Operational Research_ , 266(1):46–56, 2018.
  * Mai and Lodi (2020) Mai, T. and Lodi, A.  A multicut outer-approximation approach for competitive facility location under random utilities.  _European Journal of Operational Research_ , 284:874–881, 2020.
  * McFadden and Train (2000) McFadden, D. and Train, K.  Mixed MNL models for discrete response.  _Journal of Applied Econometrics_ , 15(5):447–470, 2000.
  * Miyamoto et al. (2004) Miyamoto, K., Vichiensan, V., Shimomura, N., and Páez, A.  Discrete choice model with structuralized spatial effects for location analysis.  _Transportation Research Record_ , 1898(1):183–190, 2004.
  * Müller et al. (2009) Müller, S., Haase, K., and Kless, S.  A multiperiod school location planning approach with free school choice.  _Environment and Planning A_ , 41(12):2929–2945, 2009.
  * Müller et al. (2012) Müller, S., Haase, K., and Seidel, F.  Exposing unobserved spatial similarity: Evidence from g erman school choice data.  _Geographical Analysis_ , 44(1):65–86, 2012.
  * Norkin et al. (1998) Norkin, V. I., Pflug, G. C., and Ruszczyński, A.  A branch and bound method for stochastic global optimization.  _Mathematical programming_ , 83(1):425–450, 1998\.
  * Paneque et al. (2021) Paneque, M. P., Bierlaire, M., Gendron, B., and Azadeh, S. S.  Integrating advanced discrete choice models in mixed integer linear optimization.  _Transportation Research Part B: Methodological_ , 146:26–49, 2021.
  * Revelle (1986) Revelle, C.  The maximum capture or “sphere of influence” location problem: Hotelling revisited on a network.  _Journal of Regional Science_ , 26:343 – 358, 1986.
  * Shapiro (2013) Shapiro, A.  Sample Average Approximation.  In Gass, S. I. and Fu, M. C., editors, _Encyclopedia of Operations Research and Management Science_ , pages 1350–1355. Springer US, Boston, MA, 2013.
  * Vovsha (1997) Vovsha, P.  _The cross-nested logit model: application to mode choice in the Tel-Aviv metropolitan area_.  Transportation Research Board, 1997.
  * Williams (1977) Williams, H. C.  On the formation of travel demand models and economic evaluation measures of user benefit.  _Environment and planning A_ , 9(3):285–344, 1977\.
  * Zhang et al. (2012) Zhang, Y., Berman, O., and Verter, V.  The impact of client choice on preventive healthcare facility network design.  _OR Spectr._ , 34(2):349–370, 2012.