1 Introduction

Location analysis is a field of operations research with a long tradition, which deals with distribution of spatial units to meet specific objectives and requirements [5, 11]. It is widely applied in many domains of engineering, for example to design various kinds of networks (distribution, telecommunications). The key element in the location problems are utilities that express an abstract measure of distance between the suppliers and clients of the considered services. If the individual clients are independent of each other, in addition to the global efficiency, the distribution of distances plays an important role [20]. Justice (equity of distribution) becomes an additional criterion for assessing the resulting solution. This approach is especially important in decisions concerning the location of public facilities, for example hospitals, crisis management centers, schools [12], where clients (citizens) have the right to a fair public access in accordance with regulations.

Demand weights are used in location problems to express the client demand for a service thus defining the location decision output as distances distributed according to measures defined by the demand weights. Note that the model of such distribution weights allows us for a clear interpretation of demand weights as the client repetitions at the same place. Splitting a client into two clients sharing the demand at the same geographical point does not cause any change of the final distribution of distances. Therefore, the distribution model of weights is important to accommodate various demand coefficients in location problems.

Numerous models for the discrete location problem were developed. Many of them differ only in the aggregation function. It is immediately apparent when we take into account effect of the siting facilities on individuals or groups [14] and consider the multicriteria model with objectives corresponding to these individual evaluations (impacts) [20]. The most commonly used aggregation is based on the weighted mean, called the median concept, where positive demand weights \(p_i\) (\(i=1,\ldots ,m\)) are allocated to several clients

$$\begin{aligned} \mu _\mathbf{p}({\mathbf z}) = \sum _{i=1}^{m}\ p_i z_i. \end{aligned}$$
(1)

The weights are typically normalized to the total 1 (\(\sum _{i=1}^{m}\ p_i=1\)). When all weights are equal we obtain simple arithmetic average. The average objective is equivalent to the total sum, which aims to global efficiency and it might discriminate isolated and low populated sites. To overcome these difficulties, especially when the equity of distribution is important, another popular approach, called the center concept [8], is used. This objective is independent of the demand weights and corresponds to the worst outcome (the situation of the client in the worst position):\(M({\mathbf z}) = \max _{i=1,\ldots ,m} z_i\). However, the center criterion might lead to substantial increase in total distance, and thus deteriorate global efficiency. Additionally considering only maximal distance limits the possibility to differentiate various feasible solutions [13].

During the last decade a new type of objective function in location theory, called ordered median (OM) function has been developed and analyzed. It originates from early models [19, 29] of compensatory extensions of the lexicographic center approach, thus representing weighted sum of the ordered outcomes (distances). The ordered median location problems (OMP) were formulated for locations on networks [17], on the plane [30] and for general discrete location problems [16]. Some special classes of the ordered solution concepts such as k-centrum and conditional median were independently developed for location problems [27, 31]. The general OM methodology was developed [18] unifying various location models. Exact and approximate solution methods were studied [3, 4, 10]. The OM objective function unifies and generalizes most common objective functions used in location theory. In fact, the ordered median function corresponds to the Ordered Weighted Averaging (OWA) aggregation, developed by Yager [33], with the non-negative preference weights. The OWA operator is a weighted average with weights allocated to the ordered distances (i.e. to the largest distance, the second largest and so on) rather than to the distances of specific clients. When applying to optimization problems with attributes modeled by variables the OWA operator is nonlinear. Yager [34] has shown that the nature of the non-linearity introduced by the ordering operations allows one to convert the OWA optimization into a mixed integer programming problem. In [24] there was shown that the OWA optimization with monotonic weights can be formed as a standard linear program of higher dimension, thus leading to efficient solution techniques for many related problems [26]. We compared different MILP formulation of the OMP for any non-negative preference weights and examined possible improvements of the computational performance by introducing various valid inequalities [22]. MILP formulations and valid inequalities for the OWA aggregation were also studied for different combinatorial optimization problems [7].

The OM approach allows to model various unweighted location problems. On the other hand, it does not allow to allocate any demand weights to specific clients and the (weighted) median solution concept (1) cannot be expressed in terms of the OMP. Typical ordered median model allows weighting of several clients only by straightforward rescaling of the distance values. However, the OM approach might be extended by the incorporation of the demand weights by rescaling accordingly clients measure within the distribution of distances, which is equivalent to the so-called Weighted OWA (WOWA) aggregation [32] using two sets of weights: the preference (OWA type) weights and the demand (distribution measure) weights. Such a Weighted Ordered Median Problem (WOMP) covers as special cases both the weighted median solution concept defined with the demand weights (in the case of equal all the preference weights), as well as the OM solution concept defined with the preference weights (in the case of equal all the demand weights).

This paper studies basic properties of the WOMP taking into account the demand weights following the WOWA aggregation rules. Linear programming formulations were introduced for optimization of the WOWA objective with monotonic preference weights thus representing the equitable preferences. We extend it to general MILP models of the WOMP for any non-negative preference weights. We examined the computational performance of WOMP and consider the possibility of improving it by introducing various additional constraints. The paper is organized as follows. In the next section we use the location problem as the multiobjective optimization problem with objectives corresponding individual clients evaluations of the location schemes to introduce the OMP and WOMP concepts. The way the demand weights are included in the problem and their interpretation is discussed. In Sect. 3 we analyze mathematical programming formulations for the WOMP and possible reinforcements of the model with valid inequalities. Section 4 describes the computational experiments and analyzes the obtained results. In Sect. 5 we conclude with main observations and propose some future research steps.

2 Problem description

We consider discrete location problem [15], which can also be defined as network location problem, where facilities are allowed to be placed only on vertices (or subset of vertices) of the underlying network. We assume no capacity limit of facilities. There is given a set of m sites (e.g. clients) and a set of potential facility locations. Without loss of generality it can be assumed that these two sets are identical. We have to place n facilities (\(n \le m\)) and assign them to clients to meet the demand. We aim at optimizing a given objective function, which is usually based on distances (costs) between the clients and the facilities. Because we consider unlimited capacities each client is assigned the closest facility. Formally the model can be expressed in the following form:

$$\begin{aligned}&\min \quad (z_1, z_2, \ldots , z_m) \end{aligned}$$
(2a)
$$\begin{aligned}&\text{ s.t. }\quad z_{i} = \sum _{j=1}^{m}{c_{ij} x_{ij}} \quad \text {for} \quad i = 1,2,\ldots ,m, \end{aligned}$$
(2b)
$$\begin{aligned}&\sum _{j=1}^{m} y_{j} = n, \end{aligned}$$
(2c)
$$\begin{aligned}&\sum _{j=1}^{m} x_{ij} = 1 \text {for} \quad i = 1,2,\ldots ,m, \end{aligned}$$
(2d)
$$\begin{aligned}&0 \le x_{ij} \le y_{j} \text {for} \quad i,j = 1,2,\ldots ,m, \end{aligned}$$
(2e)
$$\begin{aligned}&y_{j} \in \{0,1\} \text {for} \quad j = 1,2,\ldots ,m, \end{aligned}$$
(2f)

where \(c_{ij}\) denotes the cost of satisfying the total demand of client i from facility j. The main decisions are described by binary variables: \(y_{j}\) (\(j=1,2,\ldots ,m\)) is equal to 1 if a facility is placed at site j and equal to 0 otherwise. There are also binary variables that represent allocation decisions: \(x_{ij}\) (\(i,j=1,2,\ldots ,m\)) is equal to 1 if the demand of client i is satisfied by facility j and 0 otherwise. Due to lack of capacity restriction each client will be assigned to the closest facility and therefore variables \(x_{ij}\) can be relaxed to continuous variables. The auxiliary variable \(z_{i}\) (2b) expresses the cost of satisfying the demand of client i. Constraint (2c) enforces that exactly n facilities are placed. The requirement that full demand of each client is satisfied is modeled with constraint (2d). Constraint (2e) ensures that the clients are assigned to the existing facilities. Thus constraints (2c)–(2f) define the set of feasible solutions \({\mathcal F}\), which is mapped into the set of attainable outcome (cost) vectors \(\mathbf {z}\) by constraint (2b).

Further, for each client \(i=1,2,\ldots ,m\) there is also given weight \(p_i\), which determines the demand for service. We want to obtain efficient solutions of problem (2) in the sense of outcomes \(z_i = f_i(\mathbf {x})\) distribution with measures \(p_i\) for \(i=1,2,\ldots ,m\). So, intuitively, we can imagine that weight \(p_i\) scales number of clients within one location with the same value of outcome (distance) \(z_i\). It has also a direct interpretation, where different locations correspond to cities and weights \(p_i\) express number of clients in these cities. It differs substantially from standard approach, which uses weights \(p_i\) to scale the distances, thus to define the outcomes as \(z_i = p_i f_i(\mathbf {x})\) for \(i=1,2,\ldots ,m\) with a uniform distribution (with single client at each site).

Example 1

To illustrate the difference between these two approaches let us consider a simple problem with 3 locations, where we have to place one facility to minimize the average distance within one third of clients in the worst position. In other words we need to minimize the average of one third of the largest outcomes.

In Table 1 all three feasible solutions are presented for a given distance (cost) matrix \(\mathbf {c}\) and demand weights \(\mathbf {p}\). In case of outcomes distribution integer weights could be interpreted as client multiplication within one location. Thus we can consider extended vector \({\mathbf z}\), where each component corresponds to single client after multiplication. The optimal decision (with the lowest value of objective) in the sense of outcomes distribution is to place the facility in the second location, while in case of distance scaling in the first location.

Table 1 Comparison of outcomes distribution and distances scaling weighting schemes
Full size table

In practice, the distance scaling may be implemented within the individual objective functions \(f_i\). It leads to an equivalent problem without explicit weights but with transformed distance matrix (rows multiplied by weights). Therefore, it can be solved by the basic formulation of the location problem.

Direct interpretation of integer weights within optimization of outcomes distribution allows to disaggregate the problem to basic form, where demand weights for all clients are equal to \(p_i=1\). Similarly, one can proceed with any rational weights by disaggregation to clients with equal demand weights (not necessary equal to 1). Such transformation is possible, but in practice usually causes significant increase in size of the problem (number of clients) and thus made the problem impossible to solve. Our approach can directly take into account the demands weights, without the need for disaggregation.

Specific solution concepts depend on aggregations of the multiple objective outputs (2a). In particular the median solution concept is defined by the mean aggregation (1). The OM concept is based on the OWA aggregation of attributes \({\mathbf z}=(z_1,\ldots ,z_m)\). For a given preference weights \({\mathbf w}=(w_1,\ldots ,w_m) \) such that \(w_i \ge 0\) for \(i=1,\ldots ,m\) and \(\sum _{i=1}^{m}\ w_i =1\) the OWA aggregation takes the form:

$$\begin{aligned} A_\mathbf{w}({{\mathbf z}}) = \sum _{i=1}^{m}\ w_i \theta _i({{\mathbf z}}), \end{aligned}$$
(3)

where \(\Theta ({{\mathbf z}}) = (\theta _1({{\mathbf z}}),\theta _2({{\mathbf z}}),\ldots ,\theta _m({{\mathbf z}}))\) is the ordering map, i.e. \(\theta _1({{\mathbf z}}) \ge \theta _2({{\mathbf z}}) \ge \cdots \ge \theta _m({{\mathbf z}})\) and there exists a permutation \(\tau \) of set I such that \(\theta _i({{\mathbf z}}) =z_{\tau (i)}\) for \(i=1,\ldots ,m\).

In the case of decreasing weights \(w_1\ge w_2 \ge \cdots \ge w_m\), the OWA aggregation is a convex function thus, when minimized in the OMP it models the so-called equitable preferences [21]. The latter are important for many locations problems related to public facilities and thus requiring modeling the equity preferences. On the other hand, the weighted mean (1) aggregation, the standard criterion of the median location problems, cannot be expressed as an OMP.

Following the concept of demand distribution weights, the ordered averaging model enables one to introduce demand weights by rescaling accordingly its measure within the distribution of achievements as

$$\begin{aligned} A_\mathbf{w,p}({{\mathbf z}}) = \sum _{k=1}^m w_k Q_k({{\mathbf z}}), \end{aligned}$$
(4)

where \(Q_k({{\mathbf z}})\) express the conditional means within the quantile interval \([\frac{k-1}{m},\frac{k}{m}]\). The latter can be formally represented as \(Q_k({{\mathbf z}})= m\int _{(k-1)/m}^{k/m} F_{{\mathbf z}}^{(-1)}(\xi )\,d\xi \) where quantile function \(F_{{\mathbf z}}^{(-1)}(\xi )\) is the left-continuous inverse of the left-continuous right tail cumulative distribution function (cdf):

$$\begin{aligned} F_{{\mathbf z}}(d) = \sum _{i=1}^m p_i \delta _i(d) \quad \text{ where } \quad \delta _i(d) = \left\{ \begin{array}{ll} 1 &{} \text{ if } z_{i} \ge d\\ 0 &{} \text{ otherwise } \end{array}, \right. \end{aligned}$$
(5)

which for any real (outcome) value d provides the measure of outcomes greater or equal to d. That means, \(F_{{\mathbf z}}^{(-1)}(\xi ) = \sup \ \{ \eta : F_{{\mathbf z}}(\eta ) \ge \xi \}\) for \(0 < \xi \le 1\). Putting formula on \(Q_k({{\mathbf z}})\) into (4) one gets:

$$\begin{aligned} A_\mathbf{w,p}({{\mathbf z}}) = \sum _{k=1}^m w_k m\int _{(k-1)/m}^{k/m} F_{{\mathbf z}}^{(-1)}(\xi )\,d\xi . \end{aligned}$$
(6)

As shown in [25], provided the preference weights are normalized (\(\sum _{i=1}^{m}\ w_i =1\)), the aggregation (6) meets the standard WOWA definition introduced by Torra [32]. Thus one may treat formula (6) as an alternative definition of the WOWA aggregation.

Applying the WOWA aggregation (6) to multiple objective outputs (2a) of the location problem we get the Weighted Ordered Median Problem (WOMP). In the case of equal demand weights \(p_i=1/m\), formula (6) represents the standard OMP criterion (3), since \(F_{{\mathbf z}}^{(-1)}(\xi )=\theta _k({\mathbf z})\) for \((k-1)/m \le \xi < k/m\). On the other hand, for equal preference weights \(w_k=1/m\) one gets \(A_\mathbf{w,p}({{\mathbf z}})= \int _0^1 F_{{\mathbf z}}^{(-1)}(\xi ) d\xi = \mu _{{\mathbf p}}({\mathbf z})\) thus reducing WOMP to the standard median model.

Example 2

To illustrate the concept of the WOMP let us consider a location problem with 5 sites (\(m=5\)) and the normalized demand weights \(\mathbf {p} = (0.1 ,0.2 ,0.2 ,0.4 ,0.1)\). Thus the demand needs of the second and third clients are twice the demand of the first client, and the fourth client has four times bigger demand than the first one (the demand needs of the fifth and first clients are equal). Furthermore, assume the preference weights \(\mathbf {w} = (0.4 ,0.3 ,0.15 ,0.1 ,0.05)\).

Let us consider a feasible solution with the cost (distance) vector \(\mathbf {z} = (1 ,3 ,2 ,4 ,5)\). Figure 1 illustrates the corresponding quantile function \(F_{{\mathbf z}}^{(-1)}(\xi )\), which expresses the distribution of values \(z_i\) according to measures \(p_i\). Based on the quantile function, we can calculate the averages of the ordered cost vector for the consecutive equal demand portions of 1 / 5 (corresponding to integrals \(\int _{(k-1)/5}^{k/5} F_{\mathbf {z}}^{(-1)}(\xi )\,d\xi \) for \(k=1,\ldots ,5\)). Finally, according to formula (6), we get \(A_{\mathbf {w}, \mathbf {p}}(\mathbf {z}) = 5 \cdot (0.4 \cdot 0.9 + 0.3 \cdot 0.8 + 0.15 \cdot 0.7 + 0.1 \cdot 0.5 + 0.05 \cdot 0.3) = 3.85\).

Fig. 1
figure 1

Quantile function \(F_{\mathbf {y}}^{(-1)}(\xi )\) for Example 2

Full size image

3 Optimization models for WOMP

Formally, we define the Weighted Ordered Median Problem (WOMP) as

$$\begin{aligned} \min \{ A_\mathbf{w,p}({\mathbf z}) :\ {\mathbf z}=\mathbf{f}({\mathbf x}), \ {{\mathbf x}} \in {\mathcal F} \}, \end{aligned}$$
(7)

where \(A_\mathbf{w,p}({\mathbf z})\) given by (6) is applied to the location problem (2). As (6) is equivalent to the WOWA aggregation we can exploit the results of [25] to formulate optimization model for (7) with decreasing preference weights \(w_1\ge w_2 \ge \cdots \ge w_m\).

The integrals on intervals in (6) can be replaced by the left-tail integrals according to \(\int _{(k-1)/m}^{k/m} F_{{\mathbf z}}^{(-1)}(\xi )\,d\xi = L({\mathbf z},{\mathbf p}, \frac{k}{m}) - L({\mathbf z},{\mathbf p},\frac{k-1}{m})\), where

$$\begin{aligned} L({\mathbf z},{\mathbf p},0) =0 \quad \text{ and } \quad L({\mathbf z},{\mathbf p},\alpha ) = \int _0^\alpha F_{{\mathbf z}}^{(-1)}(\xi )\, d\xi \quad \text{ for } \ 0 < \alpha \le 1. \end{aligned}$$
(8)

Therefore, the WOMP criterion may be expressed as

$$\begin{aligned} A_\mathbf{w,p}({\mathbf z}) = \sum _{k=1}^m m w_k \left( L\left( {\mathbf z},{\mathbf p},\frac{k}{m}\right) - L\left( {\mathbf z},{\mathbf p},\frac{k-1}{m}\right) \right) = \sum _{k=1}^m w^\prime _k L\left( {\mathbf z},{\mathbf p},\frac{k}{m}\right) , \end{aligned}$$
(9)

where \(w_m^\prime =m w_m\), \(w_k^\prime = m(w_k - w_{k+1})\).

Graphs of functions \(L({\mathbf z},{\mathbf p},\alpha )\) (with respect to \(\alpha \)) are concave curves, the so-called (upper) absolute Lorenz curves [21]. Due to formula (8), as quantile function \(F_{{\mathbf z}}^{(-1)}\) represents the distribution of ordered outcomes, the Lorenz term \(L({\mathbf z},{\mathbf p},\alpha )\) expresses the weighted mean of \(\alpha \) portion of the largest \({\mathbf z}\) components. Thus, as noted in [25], values of function \(L({\mathbf z},{\mathbf p},\alpha )\) for any \(0 \le \alpha \le 1\) can be found by optimization:

$$\begin{aligned} L({\mathbf z},{\mathbf p},\alpha ) = \max _{u_{i}}\ \left\{ \sum _{i=1}^{m}\ z_i u_{i} : \sum _{i=1}^{m}\ u_{i}=\alpha , \quad 0 \le u_{i} \le p_i \quad \forall \ i\ \right\} . \end{aligned}$$
(10)

The above problem is an LP for a given outcome vector \({\mathbf z}\) while it becomes nonlinear for \({\mathbf z}\) being a vector of variables. This difficulty can be overcome by taking advantages of the LP dual to (10). Introducing dual variable t corresponding to the equation \(\sum _{i=1}^{m}\ u_{i}=\alpha \) and variables \(d_{i}\) corresponding to upper bounds on \(u_{i}\) one gets the following LP dual of problem (10):

$$\begin{aligned} L({\mathbf z}, {\mathbf p},\alpha )&= \min _{t,d_{i}}\ \left\{ \alpha t +\sum _{i=1}^{m}\ p_i d_{i} : t + d_{i}\ge z_i,\ d_{i} \ge 0 \quad \forall \ i \right\} \end{aligned}$$
(11a)
$$\begin{aligned}&= \min _{t}\ \left\{ \alpha t +\sum _{i=1}^{m}\ p_i \max \{ z_i - t, 0 \}\ \right\} , \end{aligned}$$
(11b)

where the optimal value \(\bar{t}\) is the \(\alpha \)-quantile of distribution of values \(z_i\) with respect to the measures \(p_i\). Equation (11a) enables the following statement.

Proposition 1

For any vector \({\mathbf z}\), value \(\varrho \) fulfills inequality \({L}({\mathbf z},{\mathbf p},\alpha ) \le \varrho \) if and only if there exist t and \(d_i\) (\(i=1,\ldots ,m\)) such that

$$\begin{aligned} \alpha t + \sum _{i=1}^{m}\ p_i d_{i} \le \varrho \quad \text{ and } \quad t + d_{i}\ge z_i,\ d_{i} \ge 0 \quad \forall \ i. \end{aligned}$$

Following (9), in the case of equitable preferences specified by decreasing weights \(w_1\ge w_2 \ge \cdots \ge w_m\), the WOMP criterion takes the form \(A_\mathbf{w,p}({{\mathbf z}}) = \sum _{k=1}^{m} w^\prime _{k} {L}({\mathbf z},{\mathbf p},\frac{k}{m})\) with positive weights \(w_k^\prime \). Therefore, the following assertion can be proven.

Proposition 2

Optimization problem (7) with decreasing weights \(w_1\ge w_2 \ge \ldots \ge w_m\) (equitable WOMP) may be expressed as the following problem with auxiliary linear inequalities:

Above model with linear WOMP criterion is further depicted as MWLP.

In general case of WOMP with non-monotonic weights \(w_i\), one may get negative coefficient \(w^\prime _k\) in formula (9). Therefore, one cannot rely on minimization of only upper bounds \(\varrho _k\) as in Proposition 2. For negative coefficients one needs to use lower bounds on the corresponding Lorenz terms. Following (11b) and taking into account that optimal value \(\bar{t}\) is the corresponding quantile, thus one of the values \(z_i\), we get that \({L}({\mathbf z},{\mathbf p},\alpha ) \ge \varrho \) if and only if

$$\begin{aligned} \varrho \le \alpha z_{i^\prime } +\sum _{i=1}^{m}\ p_i \max \{ z_i - z_{i^\prime }, 0 \} \quad \quad \text{ for } \ i^\prime =1,\ldots ,m. \end{aligned}$$

Proposition 3

For any vector \({\mathbf z}\), value \(\varrho \) fulfills inequality \({L}({\mathbf z},{\mathbf p},\alpha ) \ge \varrho \) if and only if there exist \(u_{ii^\prime }\) and \(\bar{d}_{ii^\prime }\) (\(i^\prime , i=1,\ldots ,m\)) such that

$$\begin{aligned}&\varrho \le \alpha z_{i^\prime } + \sum _{\begin{array}{c} i=1 \\ i \ne i^\prime \end{array}}^{m}\ p_i \bar{d}_{ii^\prime }&\qquad \qquad \text {for} \quad i^\prime =1,\ldots ,m, \end{aligned}$$
(12a)
$$\begin{aligned}&\bar{d}_{ii^\prime } \le z_i - z_{i^\prime } + M u_{ii^\prime }&\qquad \qquad \text {for} \quad i^\prime \ne i = 1,\ldots ,m, \end{aligned}$$
(12b)
$$\begin{aligned}&\bar{d}_{ii^\prime } \le M (1 - u_{ii^\prime } )&\qquad \qquad \text {for} \quad i^\prime \ne i = 1,\ldots ,m, \end{aligned}$$
(12c)
$$\begin{aligned}&u_{ii^\prime } \in \{0,1\}&\qquad \qquad \text {for} \quad i^\prime \ne i = 1,\ldots ,m. \end{aligned}$$
(12d)

M is a large constant. Variables \(\bar{d}_{ii^\prime }\) correspond to \(\max \{ z_i - z_{i^\prime }, 0 \}\), which is modeled by binary variables \(u_{ii^\prime }\) representing pairwise comparisons of values \(z_i\) and \(z_{i^\prime }\). Exactly, \(u_{ii^\prime }=1\) when \(z_i < z_{i^\prime }\) and \(u_{ii^\prime }=0\) otherwise. It may be further modified to reduce number of variables and constraints by taking advantages of the symmetry for variables \(\bar{d}_{ii^\prime }\) and \(\bar{d}_{i^\prime i}\). This allows us to form a model for general WOMP.

Proposition 4

Optimization problem (7) with any non-negative preference weights \(\mathbf {w}\) (general WOMP) may be expressed as the following problem with auxiliary linear inequalities and binary variables:

$$\begin{aligned} \displaystyle \min _{\varrho _k,t_k,d_{ik},z_i,\bar{d}_{ii^\prime },u_{ii^\prime }} \quad&\displaystyle \sum _{k=1}^{m} w_k^\prime \varrho _k \end{aligned}$$
(13a)
$$\begin{aligned} \text{ s.t. }\quad&\displaystyle \frac{k}{m} t_{k} + \sum _{i=1}^{m}\ p_i d_{ik} \le \varrho _k&\qquad \qquad&\text {for} \ k=1,\ldots ,m, \end{aligned}$$
(13b)
$$\begin{aligned}&t_{k} + d_{ik}\ge z_i,\ d_{ik} \ge 0&\qquad \qquad \text {for} \ i,k=1,\ldots ,m, \end{aligned}$$
(13c)
$$\begin{aligned}&\displaystyle \varrho _k \le \frac{k}{m} z_{i^\prime } + \sum _{\begin{array}{c} i=1 \\ i \ne i^\prime \end{array}}^{m}\ p_i \bar{d}_{ii^\prime }&\qquad \qquad \text {for}\ i^\prime ,k=1,\ldots ,m, \end{aligned}$$
(13d)
$$\begin{aligned} \displaystyle&\displaystyle \bar{d}_{ii^\prime } \le z_i - z_{i^\prime } + M u_{ii^\prime }&\qquad \qquad \text {for}\ i<i^\prime = 1,\ldots ,m, \end{aligned}$$
(13e)
$$\begin{aligned}&\bar{d}_{ii^\prime } \le M (1- u_{ii^\prime } )&\qquad \qquad \text {for}\ i<i^\prime = 1,\ldots ,m, \end{aligned}$$
(13f)
$$\begin{aligned}&\displaystyle \bar{d}_{ii^\prime } \le z_i - z_{i^\prime } + \bar{d}_{i^\prime i}&\qquad \qquad \text {for}~i>i^\prime = 1,\ldots ,m, \end{aligned}$$
(13g)
$$\begin{aligned}&u_{ii^\prime } \in \{0,1\}&\qquad \qquad \text {for}\ i<i^\prime = 1,\ldots ,m \end{aligned}$$
(13h)
$$\begin{aligned}&{\mathbf z}=\mathbf{f}({\mathbf x}), \ {{\mathbf x}} \in {\mathcal F}. \end{aligned}$$
(13i)

All constraints (13b)–(13i) together represent a valid model for general WOMP. However, there is no need to use both upper and lower bound constraints for all k. The corresponding upper constraints (13b)–(13c) (the linear part), the same as in model MWLP, need to be used only for \(w^\prime _k \ge 0\). In the non-linear (integer) part of the model (13d)–(13h) the corresponding lower constraint (13d) need to be applied only for \(w^\prime _k < 0\). Constraints (13e)–(13h), which are modified version of (12b)–(12d), may be skipped only in the case of all \(w^\prime _k \ge 0\) (equitable preferences). This model with MILP formulation of WOMP criterion is further denoted as MWMIP.

Model MWMIP is also consistent with minimization of \({\mathbf z}\) — lower values of \({\mathbf z}\) lead to lower value of the objective function, even though the integer part (13d)–(13h) alone is not consistent with minimization of \({\mathbf z}\). Firstly, observe that the linear part of the model (constraints (13b)–(13c) for given \({\mathbf z}\) and k correctly determines the corresponding Lorenz term by minimization of its upper bound \(\varrho _k\) (for k where \(w'_k \ge 0\)) and it is also consistent with minimization of \({\mathbf z}\). Secondly, the integer part of the model for a given \({\mathbf z}\) vector and k correctly determines the corresponding Lorenz term by maximization of its lower bound \(\varrho _k\) (for k where \(w'_k < 0\)). Vector \({\mathbf z}\) is common for both linear and integer parts of the model. Decreasing \({\mathbf z}\) leads to lower value of \(\varrho _k\), which improves (decreases) value of the objective function components for k where \(w'_k \ge 0\) and deteriorates (increases) value of the objective function components for k where \(w'_k < 0\). The objective function (13a) as a whole is equivalent to (6) (with \(w_k \ge 0\)), which is increasing with respect to \({\mathbf z}\) (but not strictly increasing). Thus minimizing the objective function leads to minimization of \(z_i\) components. It does not concern only components that correspond to preference weights \(w_k = 0\). However, similar shortcoming concerns also simpler approaches like the center criterion. Appropriate value of such \(z_i\) components can be easily determined based on identified solution (facility locations).

Some valid inequalities can be used to strengthen the MWMIP model. First, variables \(\bar{d}_{ii^\prime }\) should be non-negative, that is,

$$\begin{aligned} \bar{d}_{ii^\prime } \ge 0 \quad \text {for} \quad i, i^\prime = 1,2,\ldots ,m. \end{aligned}$$
(14)

Consider formulation (13) for a given \(\mathbf {z}\) vector. Variable \(\bar{d}_{ii^\prime }\) is included in integer part of (13) and it is used to determine \(\varrho _k\) value for k where \(w_k < 0\). As the objective function is minimized, variables \(\varrho _k\) for k where \(w_k < 0\) are maximized with upper limit defined by constraint (13d). Hence, the right hand side of the constraint is also maximized. Moreover, by constraint (13g), there is positive dependency between variable \(\bar{d}_{ii^\prime }\) and its symmetric counterpart \(\bar{d}_{i^\prime i}\) — both can be increased without violation of constraint (13g). Thus there is always optimal solution, which satisfies constraint (14).

Next,we may notice that the linear constraints on variables \(\bar{d}_{ii^\prime }\) may be additionally strengthen by adding several transitivity relations on pairwise comparisons. When \(z_i < z_{i^\prime }\) and \(z_{i^\prime } < z_{i^{\prime \prime }}\) then \(z_{i} < z_{i^{\prime \prime }}\), which is equivalent to the following constraint

$$\begin{aligned} u_{ii^{\prime \prime }} \ge u_{ii^{\prime }} + u_{i^{\prime }i^{\prime \prime }} - 1 \quad \text {for} \quad i,i^{\prime },i^{\prime \prime } = 1,2,\ldots ,m; i < i^{\prime } < i^{\prime \prime }. \end{aligned}$$
(15)

Constraint (15) can be regarded as a lower bound on binary variables arising from the transitivity relations. Similarly, one can add upper bound, which corresponds to the following relationship: if \(z_i \ge z_{i^{\prime }}\) and \(z_{i^{\prime }} \ge z_{i^{\prime \prime }}\) then \(z_{i} \ge z_{i^{\prime \prime }}\). The equivalent constraint can be stated as

$$\begin{aligned} u_{ii^{\prime \prime }} \le u_{ii^{\prime }} + u_{i^{\prime }i^{\prime \prime }} \quad \text {for} \quad i,i^{\prime },i^{\prime \prime } = 1,2,\ldots ,m; i < i^{\prime } < i^{\prime \prime }. \end{aligned}$$
(16)

It should be emphasized, however, that the transitivity relation generates huge number of constraints.

4 Computational tests

The experimental procedure has been analogous to that presented in [4]. In order to check the computational performance of the presented models and their different formulations, we have applied them to various instances of the location problem. To generate various instances we have considered some parameters characterizing the location problem and have determined their sets of possible values. Then, based on combinations of these parameters various instances of the location problems have been defined. We have considered the following parameters: the number of sites (clients) m, the number of facilities to be placed n and the type of problem defined by the vector of preference weights \(\mathbf {w}\). Besides these, we have also generated additional parameter \(\mathbf {p}\) corresponding to the demand requirements.

The size of the problem is determined by the number of sites (clients) — seven values are considered: \(m \in \{8,10,12,15,20,25,30 \}\). Due to computational complexity, general WOMP formulations are tested only on smaller sizes. The second parameter, the number of facilities n, is defined as proportional to the problem size. The following cases are examined: \(\left\lceil \frac{m}{4}\right\rceil \), \(\left\lceil \frac{m}{3}\right\rceil \), \(\left\lceil \frac{m}{2}\right\rceil \) and \(\left\lceil \frac{m}{2} + 1\right\rceil \) , where \(\left\lceil a \right\rceil \) is the smallest integer value not smaller than a. Equitable WOMP formulation is additionally evaluated on larger problems with \(m \in \{100,200 \}\) from OR-library [2].

Problem type corresponds to objective function, which is defined by the preference weighting vector \(\mathbf {w}\). This vector determines the structure and thus the complexity of the problem. We consider 6 problem types, which are described in Table 2 with respect to the number of clients m and the number of facilities n. The n-median and n-center are the most popular objective functions in multicriteria optimization. The k-centrum and \(k_1 + k_2\)-trimmed mean are less popular but also known in the field. Actually, with demand weights both T2 and T3 represent various conditional median problems [27]. As the last types T5 and T6 we consider, respectively, problems with decreasing and increasing weights.

Table 2 Problem types defined by the weighting vector \(\mathbf {w}\) with respect to the number of clients m and the number of facilities n (\(\lceil a \rceil \), \(\lfloor a \rfloor \) denote ceil and floor of a, respectively)
Full size table

The demand weights \(\mathbf {p}\) have been generated according to the Zipf distribution [35], which is typical distribution of company sizes [1] as well as population of the largest cities [6]. According to Zipf distribution the size of any object (phenomenon) is inversely proportional to its rank, when ordering the objects from the biggest to the smallest ones. Formally, it means \(p_i \sim 1/i^b\), where \(p_i\) is the size of an object in the i-th ranking position. The exponent b is very close to 1 and for the sake of simplicity it is assumed that \(b=1\). We presume that the locations are ordered by decreasing demand size, i.e. the normalized demand weights are given as \(p_i = {1}/({i \sum _{k=1}^m \frac{1}{k}})\) for \(i=1,\ldots ,m\).

For each size case we have generated 15 cost matrices with zeros on the main diagonal and the remaining entries randomly generated from a discrete uniform distribution in the interval \(\left[ 1,100\right] \). These matrices have been assigned to each combination of parameters with the corresponding problem size.

The experimental procedure has been implemented in C++ on a machine with the Intel Core2 Duo 2.53 GHz (mobile) and 3 GB of RAM. IBM ILOG CPLEX Optimization Studio (including the solver CPLEX) version 12.4 [9] has been used to solve optimization problems. A time limit of 600 seconds has been imposed as the maximum solution time for a single instance of the location problem. While presenting the average computational times for small problems, upper index in front of the time is the number of instances (of 15) that exceeded the time limit. For large problems if the instance was not solved at all, minus sign is placed.

Table 3 Average solution times [s] for MWLP (small problems)
Full size table

When the preference weights \(w_k\) (\(k=1,\ldots ,m\)) are non-increasing, then all weights \(w_k^\prime \) (\(k=1,\ldots ,m\)) are non-negative and the model reduces to MWLP. Computational results for \(m=25\) and \(m=30\) are presented in Table 3. As one can see, model MWLP copes quite well with problems up to 30 locations. The longest times concerns n-center problems and are about a few seconds. For other types of problems with non-increasing weights the solution times are shorter, reaching the shortest values for n-median problems. In Table 4 the results of MWLP are also presented for larger problems with 100 and 200 locations from OR-library [2]. The worst results relate to problems of types T2 and T3 while T5 was solved only for 100 locations [23]. It shows the computational limits of linear formulation with general MILP solver.

Table 4 Solution times [s] for MWLP (large problems)
Full size table
Table 5 Average solution times [s] for general WOMP
Full size table

When the preference weights are non-decreasing or non-monotonic the binary variables are required in the model. This leads to MILP models of WOMP criterion, whose computational complexity is significantly greater than LP models. We have tested computational performance of MILP models on problems with 8 and 10 locations. Three MILP models have been analyzed:

  • MWMIP—basic model (13a)–(13i),

  • MWMIP\(_1\)—model (13a)–(13i) with constraint (14),

  • MWMIP\(_2\)—model (13a)–(13i) with constraints (14) and (15)–(16).

The results are presented in Table 5. Model MWMIP\(_1\) achieves significantly shorter times than the basic model MWMIP. One can see that the problems with increasing preference weights (T6) are relatively easy to solve by the model MWMIP\(_1\). On the other hand, for problems T4 constraints (15)–(16), arising from the transitivity relation, allow to achieve slightly better results.

5 Conclusions

This paper has investigated the Weighted Ordered Median Problem (WOMP), which extends the Order Median Problem by taking into account the demand requirements according to the WOWA aggregation. This approach allows to obtain the optimal solution in terms of the distribution of outcomes given by the demand weights. In case of non-increasing preference weights, thus consistent with equitable WOMP, the WOWA criterion can be formulated by LP constructs. This formulation is based on the piecewise linear Lorenz function which expresses the weighted average of the largest costs within the fixed demand portion. In general, when the preference weights do not satisfy the monotonicity condition, we have proposed the extended formulation, which can be applied for any non-negative preference weights. However, this flexibility requires the binary variables and related constraints, which substantially increase the computational complexity, and thus significantly limit the maximum size of problems that can be solved.

The equitable WOMP (with LP model of WOWA criterion) have performed very well with small problems, up to 30 locations, which have been solved in a few hundreds to a few seconds. However, larger problems, about 100 locations, may cause some difficulties. The general WOMP (with MILP model of WOWA criterion), for any non-negative preference weights, might be strengthen by introducing the valid inequalities. Some of the proposed valid inequalities have allowed for several times reduction of the solution times, and thus all problems with 10 locations have been solved. Nevertheless, in the case of non-monotonic preference there is a need for the use of approximate method for problems of larger size. At present we are working on the adaptation of called Variable Neighborhood Search (VNS) metaheuristic [28], which was earlier successfully applied to Order Median Problems.

Models for WOWA optimization developed for WOMP can also be considered for various other problems not related to the location analysis. Although, the use of WOWA as an optimization criterion for other applications has not yet been widely recognized and studied. Both the solution properties and computational techniques for specific applications should be analyzed.