Competitive facility location problem with attractiveness adjustment of the follower on the closed supply chain

Abstract

In this paper, the problem examined concerns a firm entering a new facility in the market where facilities belonging to the firm and another competitor firm exist. The market’s reverse logistics network that has attracted growing attention with stringent pressures from environmental and social requirements is considered. The firm wants to maximize its profit by finding the best location and the most attractive facility to open. The other firm (the competitor) can counteract this situation and attempt to maximize its own profit by adjusting the attractiveness of its existing facilities. The demand is assumed to be aggregated at certain points in the plane and the facilities of the firm can be located at pre-determined candidate sites. To do this, Huff’s gravity-based rule in modeling the behavior of the customers is employed where the fraction of customers at a demand point that visit a certain facility is proportional to the facility attractiveness and inversely proportional to the distance between the facility site and demand point. A mathematical bi-level mixed-integer nonlinear programming model where the firm entering the market is the leader and the competitor is the follower is delivered. Furthermore, for finding the optimal solution of this model, it is converted into a one-level mixed-integer nonlinear program so that it can be solved by global optimization methods. The computational results on some examples obtained from random generated problems show that the method is able to solve the Bi-Level Programming Problem (BLPP) efficiently in a reasonable time.

Keywords:

• competitive location problems
• attractiveness adjustment
• closed supply chain
• leader–follower problem
• convexify method
• B&B method

Public Interest Statement

Competitive facility location (CFL) problems differ from the classical facility location problems considered in the operations research area because they incorporate the competition among the facilities belonging to different firms. The new facility or facilities to be located by a firm have to compete with the facilities of the other firm(s) that are already (or will be) present in the market in order to capture market share. The nature of the competition is the primary factor determining the class of the CFL problem.

In this paper, the proposed model is applicable for sales and marketing systems which encompass at least two firms as competitors. The solution methodology and the proposed model can be used separately as a basic way for using in the same work in mathematics, engineering science, or even management branches.

1. Introduction

In many competitive location problems, a firm wants to open a set of new facilities to provide goods to customers of a given geographical area where other competing firms offering the same goods are already present. The new facility or facilities to be located by a firm have to compete with the facilities of other firm(s) that are already (or will be) present in the market in order to capture market share. Many factors influence such decisions, including the available budget and expansion strategy of the firm, the location of the firm’s existing facilities and their performance, location of competitors’ facilities, the market size, demographics (e.g. population, age, gender, purchasing power) of the demand-generating population, target customer segments, availability and logistical suitability of candidate locations, and existence of other nearby attractions for generating traffic.

Huff (1964, 1966) proposed one basic formulation for estimating the market share among the competing facilities. He assumes that the probability that a customer patronizes a certain facility is proportional to its floor area and inversely proportional to a certain power of the distance from it. Huff’s model has been improved in (Cooper & Nakanishi, 1974) by replacing the floor area by a product of attraction factors. One application of this model has been found in (Jane & Mahajan, 1979) and then (Hadgson, 1981) which suggested to replace the power of the distance by an exponential function of the distance which accelerates the distance decay. Literature reviews encompass “the probability of consumers patronizing a facility” called attractiveness. It means that a customer patronizing certain facility is proportional to its attraction toward the facility. This paradigm is expressed as a parameter called the facility quality (depends on the characteristics of the facility), divided by a non-negative nondecreasing function of the distance from it. (For more details see (Drezner, 1994), (Fernández, Pelegrín, Plastria, & Tóth, 2007)) In their models, the market share captured by a facility depends on its quality and its location.

In this paper, it is supposed that the firm as leader wants to locate a single new facility in a given area of the plane. In this chain, already m + s facilities exist which offer the same good or product. It is further supposed that the first m facilities belong to the leader and the next s facilities belong to the competitor firm as follower. This problem includes competitive facility location CFL problem. It was assumed that the demand is inelastic and concentrated at n demand points; these demand points include both demand points for new products (delivered from facilities to customers) and used products (delivered from customers to facilities).

The primary works comprising of the CFL problem are (Hakimi, 1983) and (Revelle, 1986). In these problems, the new facilities must be located on a network space to compete with a number of existing (competitor) facilities. The authors use the assumption that a customer always visits the nearest open facility. According to Huff, customers’ patronization of different facilities depends on the attractiveness of each facility and is also inversely proportional to the distance to the facility. This means that customers may not always choose the nearest facility; they may sometimes visit a distant facility because it is more attractive. This idea has been extended a step further by introducing additional attributes for GIS-Based by (Nakanishi & Cooper, 1974) (a geographic information system or geographical information system (GIS) is a system designed to capture, store, manipulate, analyze, manage, and present all types of spatial or geographical data).The optimization attitude for Competitive Multi-Facility Location-Routing Problem is calculating facility attractiveness in the form of a multiplicative interaction model. In this work, this rule is applied. Moreover, our considered chain includes forward and reverses logistics operations via three kinds of facilities; forward processing facility, backward processing facility, and hybrid processing facility. New products are delivered to a number of geographically dispersed customers from plants via forward processing facilities. Used products are taken back from the customers and shipped to the plants via backward processing facilities for the purpose of recovery or safe disposal. In this paper, instead of only handling separate forward processing and collection facilities, a type of intermediate depot, namely hybrid processing facility, is also taken into account. Both new products and used products can be transferred via hybrid processing facilities. Advantages of building such hybrid processing facilities might include cost savings and pollution reduction as a result of sharing material handling equipment and infrastructure.

We found how to manipulate the B&B method by reviewing existing B&B methods in (Dempe & Gadhi, 2007) for solving Bi-level programming problem (BLPP). However, B&B methods require linear or convex quadratic nonlinear functions in the lower level problem of the bi-level programming (BP) models. In our bi-level model, the objective function of the follower’s problem is nonlinear and non-quadratic. Thus, for applying a solution method, the mixed-integer nonlinear BP formulation is firstly converted into a one-level mixed-integer nonlinear programming (MINLP) problem by using the fact that the follower’s problem is a continuous concave maximization problem. Therefore, the Karush–Kuhn–Tucker optimality conditions make it necessary and efficient. Although the NLP relaxation of resulting equivalent one-level MINLP problem is not a concave programming problem, we transformed it by introducing new variables so that it can be solved for global optimality by applying the General Structure Mixed-Integer Nonlinear a BB (GMIN-aBB) algorithm which is based on parameters and branch and bound method.

The rest of this paper is organized as follows: Section 2 describes the model formulation. In Section 3, an insight into the solution method is provided. In Section 4, concluding remarks and related examples are summarized and finally in Section 5 some conclusions are drawn.

2. Model description

In this section, the problem faced by the firm during the course of entering a single facility into the market is at first described. Then, BP formulation is delivered in which the market entrant firm that attempts to enter (locate) a facility is the leader and the other competing firms (referred to as the competitor) existing in the market are the followers. The firm wishes to determine the optimal location and attractiveness of the new kind of facility to be opened so as to maximize its profit. In addition, there are m existing facilities belonging to the leader and sexisting facility belongs to the competitor. It is assumed that the follower reacts to the market entry of the new facility (by leader) via adjusting the attractiveness of some or all of its existing facilities for maximizing its profit. Attractiveness level of a facility is determined by a function of its attributes related to the facility. For example, when the considered facility is a shopping mall, these factors include the variety of stores, availability of food, court/restaurants, adequate parking, accessibility by public transportation, selling price of products and existence of movie theaters. (Note that decreasing the attractiveness to zero means that the facility is shut down).

Let J denote the set of customer locations and j ∊ J indicate individual customer locations which is considered as aggregate population demand centers in our model. Furthermore, let E, M and S denote the set of potential location, leader and follower facilities, respectively. In our problem setting, there are n₁ demand points of new products and there are n₂ = n – n₁ demand points of used products which are indexed by j.

The other parameters and decision variables of the model are as follows:

hj	=	Annual buying power at point j of both new and used products
cl	=	Unit attractiveness cost at candidate site l
${\overline{c}}_k$	=	Unit attractiveness cost or revenue of follower’s facility at site k
$f_l^f$	=	Annualized fixed cost of opening a forward processing facility at site l
$f_l^h$	=	Annualized fixed cost of opening a hybrid processing facility at site l
$f_l^r$	=	Annualized fixed cost of opening a backward processing facility at site l
$u_l^f(u_l^h,u_l^r)$	=	Maximum attractiveness level for a forward (hybrid, backward) processing facility to be opened at site l
dlj	=	Euclidean distance between candidate site land point j
${\overline{d}}_{kj}$	=	Euclidean distance between existing facility site k and point j
${\tilde{A}}_k$	=	Current attractiveness level of competitor’s facility at site k
${\overline{A}}_k$	=	Maximum attractiveness level of competitor’s facility at site k
Decision variables	=
$g_l^f(g_l^h,g_l^r)$	=	Attractiveness of the forward (hybrid, backward) processing facility opened at site l
xl	=	Binary variable which is equal to one if a facility is opened at site l, and zero if otherwise
Ak	=	New attractiveness level of competitor’s facility at site k

A location solution vector X consists of binary decision variables $x_l^f(x_l^h,x_l^r)$ where l ∊ E and $x_l^f=1(x_l^h=1,x_l^r=1)$ indicates facility l is open for forward (hybrid, returned) processing facility, i.e. new facility is established in candidate location j.

We assume that each population center patronizes an open facility belonging to set E, M, and S proportional to the facility’s attractiveness score and inversely proportional to the distance between the demand center and the facility. In other words, each demand center has a “utility” towards an open facility which is calculated as shown in the following formulation:$$u_{ij}=\frac{A_i}{d_{ij}^2}$$

Hence the patronage probabilities of a demand center j against all open facilities are formulated as follows:

$$p_{ij}=\frac{u_{ij}}{{\sum }_i{u}_{ij}}$$

Under this patronization behavior, customers visit facilities according to probabilities p_ij and contribute to the overall revenues that are collected by the leader firm and its competitor follower. Consequently, the proportion p_ij of customers at point j (formulation is the same for both the new demand and used product) who visit a new facility at site l is expressed as:

$$p_{lj}=\frac{\frac{g_l^f+g_l^h}{d_{lj}^2}}{\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}}$$

Now, we can formulate the following bi-level MINLP model P as follows:$$P:max\sum _{j=1}^{n_1}{h}_j\frac{{\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}}{{\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}}+\sum _{j=n_1+1}^n{h}_j\frac{{\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}}{{\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}}$$(1) $$-{\sum }_{l\in E}{f}^f{x}_l^f+f^h{x}_l^h+f^r{x}_l^r-\sum _{l\in E}(g_l^f+g_l^h+g_l^r)c_l$$(1) (2) $$S.T{g}_l^f\le {\sum }_{l\in E}{u}_l^f{x}_l^f\forall l$$(2) (3) $$g_l^h\le {\sum }_{l\in E}{u}_l^h{x}_l^h\forall l$$(3) (4) $$g_l^r\le \sum _{l\in E}{u}_l^r{x}_l^r\forall l$$(4) (5) $${\sum }_{l\in E}{x}_l^f+x_l^h+x_l^r=1$$(5) (6) $$x_l^f,x_l^h,x_l^r\in \left\{0,1\right\},g_l^f,g_l^h,g_l^r\ge 0$$(6) $$max\sum _{j=1}^{n_1}{h}_j\frac{{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{ij}^2}}{{\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}}+\sum _{j=n_1+1}^{n_2}{h}_j\frac{{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{ij}^2}}{{\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}}$$(7) $$-{\sum }_{k=1}^s{\overline{c}}_k(A_k-{\tilde{A}}_k)$$(7) (8) $$S.T{A}_k\le {\overline{A}}_k\forall k$$(8) (9) $$A_k\ge 0$$(9)

The objective function (1) of the leader firm has four components. The first one (resp. the second one) represents the revenue of new (resp. used) products which is collected by the new (resp. pervious) facility of leader firm. The third and the forth components represent the fixed cost and attractiveness cost associated with opening the new facility, respectively. Constraints (2)–(4) along with the binary location and restrictions (5) and (6) on the location variables $x_l^f(x_l^h,x_l^r)$ and non-negativity restrictions existing in (6) on attractiveness variables $g_l^f(g_l^h,g_l^r)$ ensure that if no facility is opened at site l, then the corresponding attractiveness $g_l^f(g_l^h,g_l^r)$ of the facility is zero and if a facility is opened at site i, then its attractiveness $g_l^f(g_l^h,g_l^r)$ cannot exceed the maximum level $u_l^f(u_l^h,u_l^r)$. The objective function of the competitor (7) contains the summation of three terms: the revenue collected by the competitor’s facilities and the cost or revenue associated with re-adjusting the attractiveness levels. Constraints (8) and (9) ensure that the new attractiveness A_k of an existing facility at site k is between zero and an upper limit A_k.

For the solution of the above MINLP model, the GMIN-aBB algorithm was employed which will be explained in detail in Section 3. This algorithm performs a pre-processing step where all the terms in the objective function as well as in the constraints are grouped into different classes such as linear, fractional, concave, bilinear, univariate convex, and general nonconcave. Then, for each term in all classes, with the exception of the linear and concave ones, a concave over-estimator is generated.

3. Solution methodologies

Global optimization of BLPP is studied in (Gümüş & Floudas, 2005). Firstly, they propose a convex relaxation of the inner problem followed by its equivalent representation via necessary and sufficient optimality conditions. The, they introduce a BB global optimal principles presented as branch and bound framework. Their problems involve twice differentiable continuous functions. The first rigorous global optimization approach for the calculation of the flexibility test which is bi-level nonlinear optimization model is introduced by (Floudas et al., 1999). They demonstrate the applicability of their model to a heat exchanger network problem, a pump and pip problem, a reactor-cooler system, and a prototype process flow sheet model. (Pistikopoulos, Georgiadis, & Dua, 2007) introduce methods based on parametric programming to transform the bi-level problem into a family of single level optimization problems which can be solved for global optimality. The authors present computational results on several small benchmark linear–linear, linear–quadratic, quadratic–linear, and quadratic–quadratic type problems. (Gümüş & Floudas, 2005) proposes two new methods for BLPP. The first one is applicable when the outer problem includes mixed integer nonlinear and inner problem encompasses twice differentiable continuous problems; the second method is applied to problems including the inner problem featuring functions which are mixed integer nonlinear in the outer variables, and linear, polynomial or multi-linear in the inner integer variables and linear in the inner continuous variables.

In this work, the aBB global optimal approach is proposed for the efficient solution of the proposed problem.

The approach is based on the concept of feasible domain relaxation and the basic principles of the global optimization algorithm a BB. The main feature of the presented approach is that it guarantees convergence to the global optimal solution for problems that are described by general nonconvex constraints. The basic framework of the applied method is described as followed:

(1)	If any of the constraints are nonconvex, the inner level of the feasibility problem is nonconvex. Therefore, the Karush–Kuhn–Tucker (KKT) optimality conditions are necessary but not sufficient to guarantee the global optimum of the inner level. Local or even suboptimal solutions may be obtained. Hence, the first step of the proposed framework involves the convexification of the feasible region defined by the inner constraints of the original problem. The basic underestimation principles of the aBB global optimization approach are utilized to underestimate the constraints nonconvex. For the convexified problem, the KKT optimality conditions are necessary and sufficient, assuming that the linear independence constraint qualification is satisfied. Under these conditions, the convexified inner level optimization problem is equivalent to its KKT optimality conditions.
(2)	The convexified inner level optimization problem is replaced with its equivalent set of equations that are defined by the KKT optimality conditions; this transforms the bilevel problem into a single level optimization problem. Note that this problem is in general a nonconvex optimization problem due to the complementarity conditions. As a result, further convex underestimation may be needed which takes place according to aBB principles and the solution of the convexified single stage problem provides an upper bound of the original problem P.
(3)	A lower bound to the problem P is determined through a feasible solution of the original Mixed Integer Nonlinear Programming. MINLP formulation is obtained by substituting the inner problem by the KKT optimality conditions (conditions which are only necessary) as proposed by Grossmann and Floudas.
(4)	The next step after establishing an upper and a lower bound on the global solution is to refine them. This is accomplished by successfully partitioning the initial region of the variables into smaller ones. The partitioning strategy involves the successive subdivision of a hyperectangle into two sub-rectangles by halving at the middle point of the longest side of the initial rectangle (bisection). In any iteration, the lower bound of the feasibility test and the flexibility index problem is the maximum over the minima found in every sub-rectangle composing of the initial rectangle. Consequently, on-decreasing sequence of lower bounds is generated by halving the sub-rectangle that is responsible for the infimum over the minima obtained at any iteration. A non-increasing sequence of upper bounds is derived by solving the nonconvex MINLP single optimization problem obtained after the substitution of the inner problem by the KKT optimality conditions, and selecting as an upper bound the minimum over all previously determined upper bounds. If at any iteration the solution of the convexified MINLP in any sub-rectangle is found to be greater than the upper bound or the solution is not feasible, this sub-rectangle is fathomed since the global solution cannot be found inside it.

Thus, as described above, the proposed procedure for solving the problem involves the following steps:

Step 1. Set the lower bound LB = − ∞, upper bound UP = ∞ and select a tolerance ɛ.

Step 2. Substitute the convexified constraints instead of the inner constraints: we will show that the inner problem of P has concave objective function with convex constraints, and therefore it does not need any convexification in this step.

Proposition 1. The objective function

$$max\sum _{j=1}^{n_1}{h}_j\frac{{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{ij}^2}}{{\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}}+\sum _{j=n_1+1}^{n_2}{h}_j\frac{{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{ij}^2}}{{\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}}-{\sum }_{k=1}^s{\overline{c}}_k(A_k-{\tilde{A}}_k)$$

is concave.

The constraints of inner problem are linear and so are convex.

Step 3. Substitute the inner optimization problem of the original problem by the KKT optimality conditions:$$\begin{array}{c}\sum _{j=1}^{n_1}{h}_j\frac{\frac{1}{{\overline{d}}_{kj}^2}{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}+\frac{1}{{\overline{d}}_{kj}^2}\left({\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}\right)}{{\left[{\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}\right]}^2}+\\ \sum _{j=n_1+1}^{n_2}{h}_j\frac{\frac{1}{{\overline{d}}_{kj}^2}{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}+\frac{1}{{\overline{d}}_{kj}^2}\left({\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}\right)}{{\left[{\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}\right]}^2}-{\overline{c}}_k+{\mathit{\lambda }}_{1k}-{\mathit{\lambda }}_{2k}=0\forall k\\ \end{array}$$

$$A_k-{\overline{A}}_k+s_{1k}=0\forall k$$

$$-A_k+s_{2k}=0\forall k$$

$${\mathit{\lambda }}_{1k}{s}_{1k}=0\forall k$$

$${\mathit{\lambda }}_{2k}{s}_{2k}=0\forall k$$

$${\mathit{\lambda }}_{1k},{\mathit{\lambda }}_{2k},s_{1k},s_{2k}\ge 0$$

where ${\mathit{\lambda }}_1=({\mathit{\lambda }}_{11},{\mathit{\lambda }}_{12},...,{\mathit{\lambda }}_{1s})$ and ${\mathit{\lambda }}_2=({\mathit{\lambda }}_{21},{\mathit{\lambda }}_{22},...,{\mathit{\lambda }}_{2s})$ are Lagrangean multiplier vectors and s₁ = (s₁₁, s₁₂, …, s_1s) and s₂ = (s₂₁, s₂₂, …, s_2s) are slack variables corresponding to constraint sets and , respectively. However, we can remove the nonlinearity caused by the complementarily constraints using the active set strategy suggested by Grossmann and Floudas (Grossmann & Floudas, 1987). For example, consider nonlinearity phrase, Grossmann and Floudas condition for this phrase stated as:$$\mathit{\lambda }\le Ma$$$$g(z)\le M(1-a)$$

where $a\in \left\{0,1\right\}$, and M is an upper bound on the slack variables s_1k, s_2k and y_1k, y_2k for k = 1, 2, …, s. Consequently, the problem P transformed the following one-level problem:$$P^{\prime }:max\sum _{j=1}^{n_1}{h}_j\frac{{\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}}{{\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}}+\sum _{j=n_1+1}^n{h}_j\frac{{\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}}{{\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}}$$(10) $$-{\sum }_{l\in E}{f}^f{x}_l^f+f^h{x}_l^h+f^r{x}_l^r-\sum _{l\in E}(g_l^f+g_l^h+g_l^r)c_l$$(10) (11) $$S.T{g}_l^f\le {\sum }_{l\in E}{u}_l^f{x}_l^f$$(11) (12) $$g_l^h\le \sum _{l\in E}{u}_l^h{x}_l^h$$(12) (13) $$g_l^r\le {\sum }_{l\in E}{u}_l^r{x}_l^r$$(13) (14) $$\sum _{l\in E}{x}_l^f+x_l^h+x_l^r=1$$(14) (15) $$x_l^f,x_l^h,x_l^r\in \left\{0,1\right\}{g}_l^f,g_l^h,g_l^r\ge 0$$(15) $$\sum _{j=1}^{n_1}{h}_j\frac{\frac{1}{{\overline{d}}_{kj}^2}{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}+\frac{1}{{\overline{d}}_{kj}^2}({\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2})}{{\left[{\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}\right]}^2}+$$(16) $$\sum _{j=n_1+1}^{n_2}{h}_j\frac{\frac{1}{{\overline{d}}_{kj}^2}{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}+\frac{1}{{\overline{d}}_{kj}^2}({\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2})}{{\left[{\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}\right]}^2}-{\overline{c}}_k+{\mathit{\lambda }}_{1k}-{\mathit{\lambda }}_{2k}=0$$(16) (17) $$A_k-{\overline{A}}_k+s_{1k}=0\forall k$$(17) (18) $$-A_k+s_{2k}=0\forall k$$(18) (19) $${\mathit{\lambda }}_{1k}-M{a}_{1k}\le 0\forall k$$(19) (20) $${\mathit{\lambda }}_{2k}-M{a}_{2k}\le 0\forall k$$(20)

(21) $$s_{1k}-M(1-a_{1k})\le 0\forall k$$(21)

(22) $$s_{2k}-M(1-a_{2k})\le 0\forall k$$(22)

(23) $${\mathit{\lambda }}_{1k},{\mathit{\lambda }}_{2k},s_{1k},s_{2k},A_k\ge 0,a_{1k},a_{2k}\in \{0,1\}$$(23)

Step 4. The resulting formulation P′ is an MINLP model. For its solution, the GMIN-aBB algorithm was employed. To do this, both upper bound and lower bound of the problem P were obtained. First, it was necessary to get rid of nonconvex terms. To do so, new variables v_1j,v_2j for j = 1, 2, …, n₁ and new variables v_3j, v_4j for j = n₁ + 1, 2, …, n were defined as follows:(24) $$v_{1j}=\frac{{\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}}{v_{2j}}$$(24) (25) $$And{v}_{2j}={\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+\sum _{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+\sum _{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}$$(25) (26) $$v_{3j}=\frac{{\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}}{v_{4j}}$$(26) (27) $$\text{And}{v}_{4j}={\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+\sum _{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+\sum _{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}$$(27) (28) $$v_{5j}=\frac{\frac{1}{{\overline{d}}_{kj}^2}{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}}{v_{2j}}$$(28) (29) $$\text{And}{v}_{6j}=\frac{\frac{1}{{\overline{d}}_{kj}^2}{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}}{v_{4j}}$$(29)

Therefore, by using these new variables P′ is written as below:(30) $$P^{″}:\text{max}\sum _{j=1}^{n_1}{h}_j{v}_{1j}+\sum _{j=n_1+1}^n{h}_j{v}_{3j}-\sum _{l\in E}{f}^f{x}_l^f+f^h{x}_l^h+f^r{x}_l^r$$(30) (31) $$S.T{g}_l^f\le {\sum }_{l\in E}{u}_l^f{x}_l^f\forall l$$(31) (32) $$g_l^h\le {\sum }_{l\in E}{u}_l^h{x}_l^h\forall l$$(32) (33) $$g_l^r\le {\sum }_{l\in E}{u}_l^r{x}_l^r\forall l$$(33) (34) $${\sum }_{l\in E}{x}_l^f+x_l^h+x_l^r=1$$(34) (35) $$x_l^f,x_l^h,x_l^r\in \left\{0,1\right\},g_l^f,g_l^h,g_l^r\ge 0$$(35) (36) $$\sum _{j=1}^{n_1}{h}_j\frac{1}{{\overline{d}}_{kj}^2}\left(\frac{v_{5j}}{v_{2j}}+\frac{v_{1j}}{v_{2j}}\right)+\sum _{j=n_1+1}^{n_2}{h}_j\frac{1}{{\overline{d}}_{kj}^2}\left(\frac{v_{6j}}{v_{4j}}+\frac{v_{3j}}{v_{4j}}\right)-{\overline{c}}_k+{\mathit{\lambda }}_{1k}-{\mathit{\lambda }}_{2k}=0\forall k$$(36) (37) $$A_k-{\overline{A}}_k+s_{1k}=0\forall k$$(37) (38) $$-A_k+s_{2k}=0\forall k$$(38) (39) $${\mathit{\lambda }}_{1k}-M{a}_{1k}\le 0\forall k$$(39) (40) $${\mathit{\lambda }}_{2k}-M{a}_{2k}\le 0\forall k$$(40) (41) $${\mathit{\lambda }}_{2k}-M{a}_{2k}\le 0\forall k$$(41) (42) $${\mathit{\lambda }}_{2k}-M{a}_{2k}\le 0\forall k$$(42) (43) $$s_{1k}-M(1-a_{1k})\le 0\forall k$$(43) (44) $$s_{2k}-M(1-a_{2k})\le 0\forall k$$(44) (45) $$-v_{1j}{v}_{2j}+{\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}=0\forall j$$(45) (45) $$-v_{3j}{v}_{4j}-{\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}=0\forall j$$(45) (47) $$-v_{2j}+{\sum }_{l\in E}\frac{g_l^f+g_l^h}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}=0\forall j$$(47) (48) $$v_{4j}-{\sum }_{l\in E}\frac{g_l^h+g_l^r}{d_{lj}^2}+{\sum }_{i=1}^m\frac{A_i}{{\overline{d}}_{ij}^2}+{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}=0\forall j$$(48) (49) $$-v_{5j}{v}_{2j}+\frac{1}{{\overline{d}}_{kj}^2}{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}=0\forall j$$(49) (50) $$-v_{6j}{v}_{4j}+\frac{1}{{\overline{d}}_{kj}^2}{\sum }_{k=1}^s\frac{A_k}{{\overline{d}}_{kj}^2}=0\forall j$$(50) (51) $${\mathit{\lambda }}_{1k},{\mathit{\lambda }}_{2k},s_{1k},s_{2k},v_{1j},v_{2j},v_{3j},v_{4j},v_{5j},v_{6j},A_k\ge 0,a_{1k},a_{2k}\in \{0,1\}$$(51)

In the above formulation, nonconvex terms are eliminated, but the cost of this work is introducing new bilinear and fractional terms into the formulation. It helps to avoid the computationally intensive calculations used in the convexification procedure of the terms. To obtain the upper bound, the relaxed version of P″ needed to be solved. Firstly, it was necessary to convexify the nonconvex in it except for the linear and convex ones. The terms for which an over-estimator should be constructed are the fractional and bi-linear terms. To do this, the procedure given by ( Floudas, 2000 ) was used:

Let v_1jv_2j = w_1j, v_3jv_4j = w_2j, $\frac{v_{1j}}{v_{2j}}=w_{3j}$, $\frac{v_{3j}}{v_{4j}}=w_{4j}$, $\frac{v_{5j}}{v_{2j}}=w_{5j}$, $\frac{v_{6j}}{v_{4j}}=w_{6j}$, v_5jv_2j = w_7j, and v_6jv_2j = w_8j. Then, to overestimate the fractional terms the following linear constraints are added to P″:(52) $$\frac{v_{1j}^L}{v_{2j}}+\frac{v_{1j}}{v_{2j}^U}-\frac{v_{1j}^L}{v_{2j}^U}-w_{3j}\le 0\forall j$$(52) (53) $$\frac{v_{3j}^L}{v_{4j}}+\frac{v_{3j}}{v_{4j}^U}-\frac{v_{3j}^L}{v_{4j}^U}-w_{4j}\le 0\forall j$$(53) (54) $$\frac{v_{1j}^U}{v_{2j}}+\frac{v_{1j}}{v_{2j}^L}-\frac{v_{1j}^U}{v_{2j}^L}-w_{3j}\le 0\forall j$$(54) (55) $$\frac{v_{3j}^U}{v_{4j}}+\frac{v_{3j}}{v_{4j}^L}-\frac{v_{3j}^U}{v_{4j}^L}-w_{4j}\le 0\forall j$$(55) (56) $$\frac{v_{5j}^L}{v_{2j}}+\frac{v_{5j}}{v_{2j}^U}-\frac{v_{5j}^L}{v_{2j}^U}-w_{5j}\le 0\forall j$$(56) (57) $$\frac{v_{6j}^L}{v_{4j}}+\frac{v_{6j}}{v_{4j}^U}-\frac{v_{6j}^L}{v_{4j}^U}-w_{6j}\le 0\forall j$$(57) (58) $$\frac{v_{5j}^U}{v_{2j}}+\frac{v_{5j}}{v_{2j}^L}-\frac{v_{5j}^U}{v_{2j}^L}-w_{5j}\le 0\forall j$$(58) (59) $$\frac{v_{6j}^U}{v_{4j}}+\frac{v_{6j}}{v_{4j}^L}-\frac{v_{6j}^U}{v_{4j}^L}-w_{6j}\le 0\forall j$$(59)

And the bilinear terms are overestimated using the following linear constraints:(60) $$v_{1j}^L{v}_{2j}+v_{1j}{v}_{2j}^L-v_{1j}^L{v}_{2j}^L-w_{1j}\le 0\forall j$$(60) (61) $$v_{1j}^U{v}_{2j}+v_{1j}{v}_{2j}^U-v_{1j}^U{v}_{2j}^U-w_{1j}\le 0\forall j$$(61)

(62) $$-v_{1j}^L{v}_{2j}-v_{1j}{v}_{2j}^L+v_{1j}^L{v}_{2j}^L+w_{1j}\le 0\forall j$$(62) (63) $$-v_{1j}^U{v}_{2j}-v_{1j}{v}_{2j}^U+v_{1j}^U{v}_{2j}^U+w_{1j}\le 0\forall j$$(63) (64) $$v_{3j}^L{v}_{4j}+v_{3j}{v}_{4j}^L-v_{3j}^L{v}_{4j}^L-w_{2j}\le 0\forall j$$(64) (65) $$v_{3j}^U{v}_{4j}+v_{3j}{v}_{4j}^U-v_{3j}^U{v}_{4j}^U-w_{2j}\le 0\forall j$$(65) (66) $$-v_{3j}^L{v}_{4j}-v_{3j}{v}_{4j}^L+v_{3j}^L{v}_{4j}^L+w_{2j}\le 0\forall j$$(66) (67) $$-v_{3j}^U{v}_{4j}-v_{3j}{v}_{4j}^U+v_{3j}^U{v}_{4j}^U+w_{2j}\le 0\forall j$$(67) (68) $$v_{5j}^L{v}_{2j}+v_{5j}{v}_{2j}^L-v_{5j}^L{v}_{2j}^L-w_{7j}\le 0\forall j$$(68) (69) $$v_{5j}^U{v}_{2j}+v_{5j}{v}_{2j}^U-v_{5j}^U{v}_{2j}^U-w_{7j}\le 0\forall j$$(69) (70) $$-v_{5j}^L{v}_{2j}-v_{5j}{v}_{2j}^L+v_{5j}^L{v}_{2j}^L+w_{7j}\le 0\forall j$$(70) (71) $$-v_{5j}^U{v}_{2j}-v_{5j}{v}_{2j}^U+v_{5j}^U{v}_{2j}^U+w_{7j}\le 0\forall j$$(71) (72) $$v_{6j}^L{v}_{4j}+v_{6j}{v}_{4j}^L-v_{6j}^L{v}_{4j}^L-w_{8j}\le 0\forall j$$(72) (73) $$v_{6j}^U{v}_{4j}+v_{6j}{v}_{4j}^U-v_{6j}^U{v}_{4j}^U-w_{8j}\le 0\forall j$$(73) (74) $$-v_{6j}^L{v}_{4j}-v_{6j}{v}_{4j}^L+v_{6j}^L{v}_{4j}^L+w_{8j}\le 0\forall j$$(74) (75) $$-v_{6j}^U{v}_{4j}-v_{6j}{v}_{4j}^U+v_{6j}^U{v}_{4j}^U+w_{8j}\le 0\forall j$$(75)

where $v_{tj}^L$ and $v_{tj}^U$ are the lower and upper bounds of the on $v_{tj},t=1,...6$, respectively. Therefore, the upper bound of the problem P″ will be obtained by solving the following convexified problem P‴ at a given node of the B&B tree when binary variables are relaxed.

$$P^{″\prime }:max\sum _{j=1}^{n_1}{h}_j{v}_{1j}+\sum _{j=n_1+1}^n{h}_j{v}_{3j}-\sum _{l\in E}{f}^f{x}_l^f+f^h{x}_l^h+f^r{x}_l^r-\sum _{l\in E}(g_l^f+g_l^h+g_l^r)c_l$$(76) $$S.T\sum _{j=1}^{n_1}{h}_j\frac{1}{d_{kj}^2}\left(w_{5j}+w_{3j}\right)+\sum _{j=n_1+1}^{n_2}{h}_j\frac{1}{d_{kj}^2}\left(w_{6j}+w_{4j}\right)-{\overline{c}}_k+{\mathit{\lambda }}_{1k}-{\mathit{\lambda }}_{2k}=0\forall k$$(76) (77) $${\mathit{\lambda }}_{1k},{\mathit{\lambda }}_{2k},s_{1k},s_{2k},v_{tj},w_{lj},t=1,..6,l=1,..,8,A_k\ge 0,0\le a_{1k},a_{2k}\le 1$$(77)

And Equations (31)–(35) and (37)–(50) and (52)–(75)

In addition, at the root node of the tree, a local optimal solution is found for the original NLP, which gives a lower bound. In fact a lower bound is obtained by finding a local optimal solution to the problem using, for example, a commercial solver CPLEX that can handle a nonlinear programming problem. To guarantee ɛ-convergence from the global maximum, i.e. the difference between the lower and upper bounds is within ɛ of the optimal objective value, the feasible region of the problem is divided at each node such that the rectangles defined by the lower and upper limits on the decision variables are partitioned into smaller ones. This partitioning constitutes the branching mechanism in the B&B tree only.

Algorithm

(1)	Set the lower bound LB = − ∞, Obtain the relaxed problem P′ by relaxing all integer variables. While there are active nodes in the tree, repeat the following steps:
(2)	Select an active node according to some branching rule. (Branching at a node is performed by considering the solution at that node and selecting the relaxed binary variable whose value is the closest to 0.5.)
(3)	Apply the aBB algorithm and obtain an upper bound at the active node from the solution of P′.
(4)	If a feasible solution is obtained, let LB = UB, where LB indicates a lower bound to the problem. Update ZLB by $Z_{LB}=\text{max}\left\{Z_{LB},LB\right\}$, prune that node and backtrack. Otherwise, if an infeasible solution is obtained or UP ≤ ZLB, prune that node and backtrack. Report the solution with the objective value ZLB as the optimal solution.

4. Computational analysis

The performance of the presented approach was verified experimentally. It was implemented in MATLAB software. All tests were conducted for a fixed number of iterations in random test instances for assessment. The instances were classified based on the number of potential facilities and number of customers. The following policy was applied for generating test instances (for each data-set, five different data instances were created):

•	The number of demand points belonging to the {5, 10, 20, 30} where every point encompasses both used and new products
•	The number of candidate sites s belonging to {2, 5, 7}
•	The number of existing facility of competitors belonging to {2, 3}
•	Points representing the facilities and clients are uniformly distributed in [0,100] × [0,100]
•	The annualized buying power hj of customers at point j of new product are uniformly distributed in [10, 10,000] with normal probability distributions in N(50.05, 16. 52) and the annualized buying power hj of customers at point j of used product are uniformly distributed in [5, 5,000] and have normal probability distributions in N(40, 10.52)
•	The attractiveness cost ci of candidate facility for the leader at point i are real-valued quantities and uniformly distributed in [0.5,5] with normal probability distributions in N(5, 1.52)
•	The attractiveness cost \tilde{c}_{k} of facility of follower at point k are real-valued quantities and are uniformly distributed in [0.5, 5] with normal probability distributions in N(5, 1.52)
•	The current attractiveness cost \tilde{A}_{k} of facility at point k are real-valued quantities and are uniformly distributed in [10, 1,000] with normal probability distributions in N(500, 1402)
•	The fixed costs f_{i}^{f} ,f_{i}^{h}and f_{i}^{r} are set as 800ci,1,100ci and 550ci respectively
•	Upper bounds u_{l}^{f} ,u_{l}^{h}and u_{l}^{r} are taken as 7,500ci, 8,000ci, and 7,000ci respectively, and upper bound \bar{A}_{k} is assigned a value of 7,500\tilde{c}_{k}.

As can be seen, GMIN-aBB is computationally expensive especially for very small values of the convergence parameter ɛ. To improve the running time efficiency of GMIN-aBB, instead of terminating the aBB algorithm when the difference Z_UB − Z_LB ≤ ɛ%, the iterations are stopped as soon as the improvement in the gap Z_UB − Z_LB between two non-successive iterations becomes less than a user-specified threshold value δ, namely the exploration of BB tree terminated when$$\frac{\left(Z_{UB}^{t-k}-Z_{LB}^{t-k}\right)-\left(Z_{UB}^t-Z_{LB}^t\right)}{Z_{UB}^{t-k}-Z_{LB}^{t-k}}\le \mathit{\delta },$$

where $Z_{UB}^t$and $Z_{UB}^{t-k}$are the upper bounds at iterations t and t – k, respectively, while $Z_{LB}^t$and $Z_{LB}^{t-k}$ are the best lower bounds at iterations t and t–k, respectively. k is a user-specified parameter that represents the number of iterations during which the improvement in the gap between Z_UB and Z_LB is measured. In our examples, k = 4 and δ = 0.1.

We present the results on the test instances in the following Tables. In these tables, $Z_L^{best}$ indicates the best lower bound of the leader obtained by the GMIN-aBB algorithm run with the above-mentioned approximation strategy, and $Z_F^{best}$, stands for the profit realized by the follower as a reaction to the leader.

The computational results in Table 1 show that the GMIN-aBB method solves the proposed model very efficiently and the model is extremely efficient in finding integer solutions. Most of the time, GMIN-aBB method was not invoked to obtain integer solutions, which shows that our model is integer- friendly. It is worth noting that our method for obtaining a lower bound often generated an optimal solution or a high-quality approximate solution with only one or two locations different from the optimal location. Furthermore, for showing the affectedness of attractiveness adjustment of the follower facilities, the advantages of the leader and the disadvantages of the competitor were firstly quantified when the leader takes into consideration the reaction of the competitor (fixed A_k). In the second stage, the aim was to quantify the benefit of the leader and the benefit of the competitor when the leader takes into consideration the reaction of the competitor but when the upper bound of the attractiveness is changed to half $({\overline{A}}_k/2)$. For each instance, the percentage of loss from the competitor was computed.

Table 1. computational results obtained by GMIN-aBB

Instance(n, m, r)	$Z_L^{best}$	$Z_F^{best}$	Computation time solution
(5, 2, 2)	162,183.6007	235,145.6705	12.929
	36,349.04157	157,060.0119	31.402
	66,631.85411	58,711.14141	38.603
	89,943.47257	223,194.9707	38.526
	138,772.4512	56,199.83209	11.546
(5, 2, 3)	89,739.0741	75,874.64574	9.728
	173,197.0197	202,176.3623	7.789
	127,315.7748	173,560.5099	35.743
	163,068.2135	285,158.7187	18.609
	154,481.1468	225,034.012	7.734
(5, 5, 2)	203,396.9806	144,857.8099	23.755
	105,957.5555	253,223.8754	32.685
	86,379.93404	68,352.56194	5.314
	40,739.87958	172,223.5862	5.56
	201,139.1605	214,540.4041	7.532
(5, 5, 3)	62,911.88963	90,561.81037	22.072
	113,893.9911	247,933.2591	28.99
	85,188.16331	52,256.00463	38.594
	84,427.9808	273,325.108	15.525
	175,078.0963	301,299.9646	45.921
(5, 7, 2)	138,233.9526	221,840.5653	19.901
	233,885.4458	44,592.87104	34.749
	263,980.0847	182,909.0165	53.742
	122,832.3421	172,596.3891	48.382
	241,437.2665	269,588.1328	56.482
(5, 7, 3)	69,482.35423	53,840.55126	52.555
	193,436.4351	211,785.1145	67.448
	218,915.7749	175,336.0339	28.899
	153,821.2247	214,380.0306	17.932
	134,490.9615	181,012.1815	17.941
(10, 2, 2)	70,604.33593	271,916.4984	78.562
	188,830.8969	291,307.906	72.487
	52,431.17585	264,983.3772	99.124
	223,174.397	184,677.3347	57.186
	171,770.5513	55,232.11059	78.522
(10, 2, 3)	114,472.1819	129,969.3082	60.394
	59,810.44828	179,560.6885	59.202
	96,562.32776	130,097.7135	101.514
	75,264.63391	275,796.4832	103.686
	137,495.253	209,843.7231	97.645
(10, 5, 2)	330,410.1511	222,858.3721	121.385
	281,186.4398	317,812.4356	90.415
	122,673.6612	339,240.4703	98.017
	212,524.2916	282,221.2561	124.526
	108,363.6277	111,240.4921	52.729
(10, 5, 2)	335,998.9269	141,323.166	124.512
	299,940.1545	226,142.0719	62.797
	284,357.9325	134,539.5229	95.809
	117,829.0574	141,430.8419	92.808
	106,278.132	51,260.48358	137.385
(10, 7, 2)	109,500.8225	276,289.3792	118.054
	296,213.3181	128,327.8353	87.64
	363,374.5414	125,776.3877	150.382
	375,706.9231	227,135.0875	157.824
	397,676.5038	251,334.3833	76.778
(10, 7, 3)	372,090.7241	168,049.5399	136.328
	189,520.9862	295,089.6073	103.129
	67,882.32438	300,067.2645	92.484
	186,791.8758	477,048.4731	41.554
	342,006.7619	281,694.029	102.742
(20, 2, 2)	199,654.8371	359,331.5461	48.299
	62,380.77262	236,306.7942	77.227
	95,416.56812	81,357.14524	96.299
	120,847.4246	340,503.0967	96.097
	174,115.4013	77,400.55178	34.636
(20, 2, 3)	118,277.0718	114,184.2124	165.227
	269,544.5983	93,956.01716	72.66
	186,384.8419	385,585.4354	54.061
	251,186.137	206,843.1408	80.703
	235,622.0786	93,386.83148	103.082
(20, 5, 2)	256,457.0542	271,809.5933	13.227
	521,261.3496	363,142.3941	38.208
	69,508.14009	83,214.97575	54.142
	323,328.3122	85,730.53802	43.77
	426,734.0701	105,891.3373	68.188
(20, 5, 3)	294,766.428	319,680.5018	126.679
	513,666.5635	376,284.4593	229.163
	241,483.9162	454,857.6198	261.407
	548,986.0524	266,111.7024	110.178
	587,898.447	514,811.3983	237.254
(20, 7, 2)	42,709.65389	315,852.3485	247.766
	370,151.9575	485,110.7658	284.165
	345,738.3904	437,345.0982	191.173
	575,351.498	509,527.9404	163.559
	86,137.15274	467,642.1405	142.926
(20, 7, 3)	556,747.484	546,487.4743	389.231
	448,144.1009	342,408.2299	105.065
	564,479.2747	284,345.5649	460.78
	465,056.7041	284,392.409	353.689
	199,654.8371	359,331.5461	128.299
(30, 2, 2)	197,218.2698	400,247.3211	341.568
	468,154.1387	431,408.0801	316.507
	155,571.4447	389,106.2362	426.387
	546,857.9932	260,059.5321	253.391
	429,057.5134	52,049.28754	341.404
(30, 2, 3)	513,975.5685	262,227.7069	208.015
	343,157.6509	293,709.8753	201.347
	458,007.2743	262,309.2226	492.373
	391,451.981	354,803.3854	424.872
	585,922.6657	312,934.4984	412.929
(30, 5, 2)	824,098.3337	182,133.6065	231.402
	759,456.5808	288,156.2301	238.603
	531,479.4225	130,667.8153	438.526
	479,196.0078	348,556.6453	211.546
	402,414.2866	270,371.3102	29.728
(30, 5, 3)	1018,464.901	245,159.8741	47.789
	907,927.4439	333,519.9054	135.743
	551,968.9233	353,459.8821	148.609
	753,738.76	300,400.3355	467.734
	519,834.1113	141,293.2357	203.755
(30, 7, 2)	726,401.7239	311,873.8549	405.314
	682,214.9942	178,286.8042	215.56
	209,986.6841	188,336.6445	527.532
	177,231.5601	56,838.20522	422.072
	360,750.0297	473,714.6627	728.99
(30, 7, 3)	707,342.0875	223,085.925	838.594
	832,012.6217	218,764.0853	415.525
	854,905.0668	390,453.3115	45.921
	895,686.9233	431,443.9555	19.901
	162,183.6007	235,145.6705	12.929

The conclusion drawn from the results is that the competitor almost takes disadvantage of the leader’s effort when the attractiveness of competitor’s facility is decreased or fixed. In fact, with A_k decreasing, the objective value increases while the market share captured by Firm B shrinks, and the optimal locations tend to be nodes which are closer to the warehouse/plant since the transportation cost becomes the main cost (Table 2).

Table 2. The loss of the competitor, in the case of $(fixed{\overline{A}}_k)$ and $({\overline{A}}_k/2)$

Instance (n, m, r)	Loss (%) $(\text{fixed}{A}_k)$	Loss (%) $({\overline{A}}_k/2)$	Instance (n, m, r)	Loss (%) $(\text{fixed}{A}_k)$	Loss (%) $({\overline{A}}_k/2)$
(5, 2, 2)	−71.2528	−26.7198	(20, 2, 2)	−65.0711	−24.4017
	−54.6343	−20.4879		−5.15566	0
	−43.8187	−16.432		−52.3502	−19.6313
	3.29196	1.23449		−4.4203	−1.65761
	−72.4594	−27.1723		−18.5769	−6.96633
(5, 2, 3)	−55.606	−20.8522	(20, 2, 3)	−40.0674	−15.0253
	0	0		−21.0946	−7.91047
	−23.7305	−8.89895		−78.0695	−29.2761
	−38.513	−14.4424		−49.0011	−18.3754
	−6.34794	0		−24.699	−9.26214
(5, 5, 2)	−5.82299	−2.18362	(20, 5, 2)	−53.3279	−19.998
	−18.7321	−7.02455		4.55657	3.70872
	−56.6354	−21.2383		−2.64937	−0.99351
	−2.28611	−0.85729		−67.5641	−25.3365
	−60.41	−22.6538		−15.3931	−5.77243
(5, 5, 3)	−8.75085	−3.28157	(20, 5, 3)	−27.1221	−10.1708
	−12.6961	−4.76103		−21.8111	−8.17915
	−15.6241	−5.85904		−14.8479	−5.56798
	−56.4283	−21.1606		−68.186	−25.5698
	−63.1049	−23.6643		−7.77452	−2.91545
(5, 7, 2)	−36.635	−13.7381	(20, 7, 2)	−62.7091	−23.5159
	−58.6022	−21.9758		−4.28833	0
	−25.311	−9.49163		−2.02326	−0.75872
	−28.8598	−10.8224		−7.85695	−2.94636
	2.16348	0.8113		−77.3864	−29.0199
(5, 7, 3)	−23.4194	−8.78227	(20, 7, 3)	−65.1301	−24.4238
	−9.48925	0		−50.1543	−18.8079
	−54.4435	−20.4163		−41.2046	−15.4517
	−15.1416	−5.67811		−14.8676	−5.57536
	−41.2847	−15.4818		−70.3426	−26.3785
(10, 2, 2)	−32.813	−12.3049	(30, 2, 2)	−69.0811	−25.9054
	−11.2048	−4.20178		−71.9083	−26.9656
	0	0		−2.18569	−0.81963
	−11.7362	−4.40109		−70.5922	−26.4721
	−47.7947	−17.923		−57.2363	−21.4636
(10, 2, 3)	−67.9163	−25.4686	(30, 2, 3)	−36.1381	−13.5518
	−41.6369	−15.6138		−5.20567	0
	−55.1919	−20.697		−44.5298	−16.6987
	−20.8229	−7.80859		−40.7943	−15.2979
	−72.8108	−27.304		−3.7339	−1.40021
(10, 5, 2)	−51.3014	−19.238	(30, 5, 2)	−73.7731	−27.6649
	−78.657	−29.4964		−73.1438	−27.4289
	−13.7076	−5.14034		−66.9033	−25.0888
	−37.6759	−14.1285		−64.4928	−24.1848
	−65.2035	−24.4513		−71.6976	−26.8866
(10, 5, 2)	−20.0422	−7.51584	(30, 5, 3)	−37.2121	−13.9546
	−66.7327	−25.0248		−63.3512	−23.7567
	−23.6039	−8.85146		−76.108	−28.5405
	−54.6415	−20.4906		−68.0937	−25.5351
	−21.5778	−8.09167		−48.574	−18.2152
(10, 7, 2)	−1.36702	−0.51263	(30, 7, 2)	−77.2526	−28.9697
	−38.0556	−14.2708		−61.5623	−23.0859
	−24.5988	−9.22454		−69.5192	−26.0697
	−1.97686	−0.74132		−39.0315	−14.6368
	−62.3955	−23.3983		−13.8508	−5.19405
(10, 7, 3)	−60.4663	−22.6749	(30, 7, 3)	−63.5692	−23.8384
	−27.6192	−10.3572		−32.8455	−12.3171
	−5.1332	0		−60.7151	−22.7682
	−73.6578	−27.6217		−2.09586	0
	−1.82969	−0.68613		−66.4812	−24.9304

5. Conclusion

In this paper, we considered a relatively new discrete CFL problem in a closed chain with two competitors, which at first attempts to add new facility to its facilities so as to maximize own profits. It is assumed that there is a second competitor in the market which reacts to the opening of the new facility by adjusting the new attractiveness levels of its facilities with the objective of maximizing its own profits. As the market includes forward and reverse logistics network, we considered three kinds of facilities in this market: forward, backward, and hybrid processing facilities. A programming model is proposed which is a bi-level model and solves the problem outlined in this paper by using a global optimization method called GMIN-aBB after converting the bi-level model into an equivalent one-level MINLP model. This conversion is possible since the objective function of the competitor, which is the follower in the bi-level programming formulation, is concave in terms of the attractiveness variables when the locations and the attractiveness levels of the leader firm are fixed. We applied this method on the random generated examples and the computation results show the efficiency of the proposed method. For future research, one can consider an extension of the proposed bi-level model where the competitor reacts to the market entrant firm in the case of adjusting the attractiveness levels and also opening new facilities and/or closing existing ones. Furthermore, a combination of available algorithms such as branch and bound and multi-parametric could be helpful.

Acknowledgments

The author is grateful to the referees for their constructive suggestions for improving this paper.

Additional information

Funding

Funding. This research was supported by Kurdistan University [grant number 3732194681].

Notes on contributors

Arsalan Rahmani

He is a faculty member of the Mathematics Department at the University of Kurdistan, Sanandaj, Iran. Other information is outlined below:

Field of Specialty:

Combinatorial Optimization, Mathematical Modeling, Stochastic Programming.

2011–2015

PhD in Applied Mathematics (Operations Research).

Department of Mathematics and Computer Science, Amirkabir University of Technology, Tehran, Iran.

2009–2011

Msc in Applied Mathematics (Operations Research).

Department of Mathematics and Computer Science, Amirkabir University of Technology, Tehran, Iran.

2008

BSc in Pure Mathematics.

University of Kurdistan, Sanandaj, Iran.

Member of Operations Research Society Iran.