LAGRANGEAN RELAXATION BASED APPROXIMATE APPROACH TO THE CAPACITATED LOCATION PROBLEM LAGRANGEAN RELAXATION BASED APPROXIMATE APPROACH TO THE CAPACITATED LOCATION PROBLEM

Cost optimal design of the most of distribution and servicing systems consists in decisions on number and locations of facilities, from which customer’s demands are satisfied. When such public or private servicing system is designed, some limits are put on abilities of the particular facilities. These constraints follow from the limited ability of the facilities to satisfy demands of customers. In spite of the fact that capacity of a facility is a very vague value, the capacities of the located facilities bring serious complications concerning a solving technique, which may be used to solve a realsized facility location problem. In contrast to an uncapacitated facility location problem, which can be solved exactly for real-sized instances containing hundreds of possible locations and thousands of customers, the capacitated location problem resists to all attempts to solve it exactly in reasonable time. To avoid this complication and to obtain a good solution of the capacitated facility location problem in sensible time, we employ approximate methods based on Lagrangean relaxation of the capacity constraints, which has several advantageous properties. One of these properties is that the relaxed problem known as the uncapacitated location problem can be solved exactly by an algorithm based on Erlenkotter’s approach [2]. Real sized instances of the uncapacitated problem were broadly tested [6], [4] and it was proved that it is possible to obtain their optimal solution in a sensible time. The next useful property of the Lagrangean relaxation, which can be exploited with the previous one, is that the objective function value of optimal solution of the relaxed problem provides lower bound of the optimal solution of the original problem. In this way, the obtained approximate solution can be compared with the estimated value of the optimal one. The third property is that any optimal solution of the Lagrangean relaxation preserves so called compactness of particular customer clusters, which arise from a customer assignment to the individual facility locations. This property is not directly connected with the facility location, but it has great impact on applications following up with a distribution system design. This property enables good service of the customers by vehicle routes starting and terminating at the facility location.


Introduction
Cost optimal design of the most of distribution and servicing systems consists in decisions on number and locations of facilities, from which customer's demands are satisfied. When such public or private servicing system is designed, some limits are put on abilities of the particular facilities. These constraints follow from the limited ability of the facilities to satisfy demands of customers. In spite of the fact that capacity of a facility is a very vague value, the capacities of the located facilities bring serious complications concerning a solving technique, which may be used to solve a realsized facility location problem. In contrast to an uncapacitated facility location problem, which can be solved exactly for real-sized instances containing hundreds of possible locations and thousands of customers, the capacitated location problem resists to all attempts to solve it exactly in reasonable time. To avoid this complication and to obtain a good solution of the capacitated facility location problem in sensible time, we employ approximate methods based on Lagrangean relaxation of the capacity constraints, which has several advantageous properties. One of these properties is that the relaxed problem known as the uncapacitated location problem can be solved exactly by an algorithm based on Erlenkotter's approach [2]. Real sized instances of the uncapacitated problem were broadly tested [6], [4] and it was proved that it is possible to obtain their optimal solution in a sensible time. The next useful property of the Lagrangean relaxation, which can be exploited with the previous one, is that the objective function value of optimal solution of the relaxed problem provides lower bound of the optimal solution of the original problem. In this way, the obtained approximate solution can be compared with the estimated value of the optimal one. The third property is that any optimal solution of the Lagrangean relaxation preserves so called compactness of particular customer clusters, which arise from a customer assignment to the individual facility locations. This property is not directly connected with the facility location, but it has great impact on applications following up with a distribution system design. This property enables good service of the customers by vehicle routes starting and terminating at the facility location.
We present two methods for obtaining suitable values of Lagrangean multipliers. The classical one is based on a sub-gradient method applied on capacity constraints after their special adjustment. The second method is designed as an adaptive method with random experiments for determination of candidates for move from a current solution to the next one. Both approaches make use of the strengthening constraint, which considerably improves quality of the obtained lower bound. Both methods were tested, compared and the associated results are reported in the section "Numerical experiments" of this paper and furthermore some explanation of their various behaviors is suggested.

Mathematical model of the solved problem
Mathematical programming approach to the capacitated facility location problem originates from the assumption that goods distribution is performed from a primary source via warehouses or other goods transshipment places to the particular customers. The number and positions of the warehouses, which will be called facilities in this paper, should be determined so that the total yearly cost for both located facilities and customer demand satisfaction are minimal. TO THE CAPACITATED LOCATION PROBLEM   LAGRANGEAN RELAXATION BASED APPROXIMATE APPROACH  TO THE CAPACITATED LOCATION PROBLEM   Jaroslav Janáček -Lýdia Gábrišová * When a distribution system is to be designed, limits on terminal capability often must be taken into account. These capacity constraints in this and other facility location problems constitute severe obstacles in exact solving processes. Within this paper, we focused on study of approximate methods based on Lagrangean relaxation of the capacity constraints, which has several advantageous properties. The first of them is that the relaxed problem, known as the uncapacitated location problem, can be solved exactly even for real sized instances [6], [4]. The second useful property of the Lagrangean relaxation is that the objective function value of the optimal solution of the relaxed problem provides lower bound of the optimal solution of the original problem. We present two methods for obtaining suitable values of Lagrangean multipliers. The classical one is based on a sub-gradient method applied on capacity constraints after their special adjustment. The second method is designed as an adaptive method with random experiments for determination of candidates for move from the current solution to the next one. These two methods were tested, compared and the associated results are reported in the concluding part of this paper. The problem is described by finite set I of possible facility locations and finite set J of the customers, whose demands should be satisfied. Associated costs and charges are connected to particular elements of these sets and to their pairs. A fixed charge for location of a facility at possible location iʦI is denoted by f i . This fixed charge includes all costs connected with keeping this facility at the location for one year. This charge does not include items, size of which depends on amount of demands, which are satisfied via this location.

LAGRANGEAN RELAXATION BASED APPROXIMATE APPROACH
The cost of j-th customer yearly demand satisfaction via facility located at place i is denoted by coefficient c ij . This coefficient includes all transportation costs for goods transport from the primary source to the facility location, from this location to the customer and manipulating cost in the facility.
It is presumed that a facility may be placed only at some place of the above-introduced finite set I of possible locations. To model the decision on placing or not placing a facility at location i, variable y i ʦ {0, 1} is introduced for each location i from set I.
Let us assume that each customer jʦJ should be supplied by a yearly amount of goods b j . To be able to express that a customer is assigned to a given facility location and that he is supplied via this location, another set of zero-one variables is established. Variable z ij ʦ {0, 1} models a decision on assigning or not assigning customer j to facility location i. Each solution of the following integer programming problem will be described by zero -one variables y i and z ij . This zero -one range of the decision variables is assumed in each model within this paper and that is why this associated obligatory constraint is not resumed at the following models.
Let us consider that a i denotes the capacity of a facility located at i then the complete model of the cost minimal capacitated facility location problem can be formulated as follows: In this integer programming model, constraints (2) ensure that each customer demand must be satisfied from exactly one facility location and constraints (3) force out the placement of a facility at location i whenever a customer is assigned to this facility location. Constraints (4) ensure that the total demand satisfied via facility location i doesn't exceed given capacity a i . Having omitted or relaxed constraints (4), the problem (1)-(3) is known as the uncapacitated facility location problem and it can be effectively solved making use of implementation of the branch and bound method with Erlenkotter's lower bounding [1]. Computa-tion behaviour of the technique was broadly examined in [4], [6], [7], [8] and it was shown that this approach is able to manage large size problems of practice. The capacitated problem (1)-(4) loses integrality property of variables z ij due to the capacity constraints and so it constitutes a very hard problem for exact solving when a real-sized instance of the problem is considered.
Model (1)-(4) can be replaced by following formulation, which is equivalent from the point of set of feasible (zero -one) solutions: In this formulation, constraints (7) become surplus, but we let them in the model for a better numerical convergence of the branch and bound method solving the relaxed problem. The performed reformulation of constraints (8) has proved to be a very efficient form for subsequent Lagrangean relaxation.
The above -mentioned feasible solution equality of the models (1)-(4) and (5)-(8) follows from the following proofs: Each zero -one solution described by variables y i and z ij , which satisfies (3) and (4), must satisfy (8). For y i ϭ 1 constraint (8) takes the same form as (4) and for y i ϭ 0 is z ij ϭ 0 for each j ʦ J and it follows that left -hand -side of (8) is equal to zero as well as the righthand -side. Let us assume that (7) and (8) hold. If y i ϭ 1, then (8) takes the same form as (4). If y i ϭ 0, then z ij ϭ 0 for each j ʦ J and (4) holds also.

Lagrangean relaxation and strengthening constraint
The studied Lagrangean relaxation of problem (5)-(8) consists in removing constraints (8) from the model and in embedding a measure of their dissatisfaction into the objective function. The measure of the capacity constraint dissatisfaction is enumerated as the difference between the capacity and the demand, which is modeled by the expression on the left-hand-side of the constraint. This overload of a facility is weighed by a nonnegative weight called Lagrangean multiplier. This multiplier represents a penalty paid for breaking the constraints.
Having denoted u i the Lagrangean multiplier of the capacity constraint i ʦ I, the Lagrangean relaxation takes the form: For fixed nonnegative values of all Lagrangean multipliers, the objective function of the Lagrangean relaxation has the form: which is the objective function of an uncapacitated location problem. Due zero -one variables y i and z ij , it is obvious that an optimal solution exists for each setting of the Lagrangean multiplier values.
The preliminary experiments with this model showed that the set of feasible solutions enlarges too much in this way of relaxation. This consequence can be minimized to some extent by addition of a simpler constraint, which is redundant with respect to the relaxed ones, but which diminishes the set of feasible solutions of the relaxed problem. The strengthening constraint is derived from the capacities and is laid upon the number of located facilities. It is obvious that at least such a number of facilities must be located to cover all customer demands. If a denotes max {a i :i ʦ I} and B denotes sum of all customer demands, then minimal number r of the necessary located facilities is r ϭ B / a and the strengthening constraint can be constructed in this way: To keep the form of the uncapacitated location problem, this constraint can be relaxed as well as the other capacity constraints using nonnegative Lagrangean multiplier v. After this reformulation we get To determine value v, for which the optimal solution of (14)-(16) satisfies complementary constraint (17), a simple dichotomy algorithm can be used.
It is true that the dichotomy algorithm does not ensure satisfaction of constraint (17) in general, but we did not meet any real instance of the solved problem, in which the algorithm failed.

Sub-gradient method
The sub-gradient method used in our approach was designed with the goal to maximize a lower bound of the objective function value of an optimal solution of the original capacitated facility location problem [1]. The method improves the lower bound by particular steps, in which it adjusts the values of the Lagrangean multipliers u i . During this iterative process, the overloads of facil-ities with dissatisfied capacity constraints are penalized by higher values of the associated Lagrangean multipliers and hence forced to minimize size of dissatisfaction. This process may lead to a complete capacity constraint satisfaction or, which is a more frequent case, to such solution, which is only slightly infeasible.
The sub-gradient method maximizes function F(u), the value of which is equal to the minimal value of expression (14) subject to (15), (16,) (17) and integer nonnegative values of variables y i and z ij for the given values of u ϭ Ͻu 1 , u 2 , …Ͼ.
Value F(u) for the given u is obtained using the above mentioned dichotomy algorithm applied on the branch and bound method developed for the uncapacitated location problem.
The algorithm of the sub-gradient method starts from an initial point from nonnegative part of |I| dimensional space of the Lagrangean multipliers. This point is usually set to the origin of the coordinates system i.e. its coordinates are set to zeros.
A move from current point u k to next point u kϩ1 is made along a direction, in which the function increases in the neighborhood of point u k . The steepest increase of F(u) in the neighborhood of point uk can be found along the gradient of this function at point u k . The coordinates of the gradient can be enumerated as values of partial derivatives of (17) or (14) or (9) by the individual multipliers u i with consequent substitution of values u i k . In this way, the i-th component takes the value of the associated located facility overload.
To follow the direction of the gradient, the move should be performed in accordance with equality Unfortunately, keeping this formula, the requirement of nonnegative multipliers should be broken. That is why the following formula is used in the algorithm That is why the move is not performed exactly along the gradient, but along the direction of sub-gradient. The length of the step is given by parameter ␣, the value of which is chosen from interval Ͻ␣ min , ␣ max Ͼ. After each step, values of F(u k ) and F(u kϩ1 ) are compared, and if F(u k ) Ն F(u kϩ1 ) holds, then the return to u k is done and the move is repeated with a lower value of ␣. The process terminates if parameter ␣ reaches value of ␣ min or if the resulting improvement of the last step is less than the given value ⑀.

Adaptive method SUPRA
The method SUPRA was designed to overcome the weak side of the sub-gradient method, which emerges when the minimization The method also improves the lower bound by particular steps in which it adjusts the values of the Lagrangean multipliers u i . Similarly to the sub-gradient method, the overloads of facilities with dissatisfied capacity constraints are penalized by higher values of the associated Lagrangean multipliers and hence forced to minimize the size of dissatisfaction.
Maximization of function F(u) is performed where the value of this function is also equal to the minimal value of expression (14) subject to (17) and integer nonnegative values of variables y i and z ij for the given values of u ϭ Ͻu 1 , u 2 , … Ͼ. As in the previous approach, the dichotomy algorithm applied on the branch and bound method is employed here to obtain value F(u) for the given u.
The algorithm of method SUPRA also starts from an initial point, which can be set to the origin of the coordinates system, i.e. its coordinates are set to zeros but a move from the current point u k to the next point u kϩ1 is made in a completely different way, which does not utilize the property of the gradient computed from the facility overloads.
The move from the current point u k to point u kϩ1 is the best found move obtained as a result of two phases, which explore neighborhood of the current point uk.
The first phase consists of a random generation of changes r j for j ϭ 1, …, s , from which s trial points u kj are derived in accordance with the formula u kj ϭ P(u k ϩ r j ), where P denoted a projection of a general point to the nonnegative part of |I| dimensional space of the Lagrangean multipliers. The values F(u kj ) are compared and the best found solution is updated. Furthermore, the statistical gradient r is enumerated in accordance with the equation: The random generation of the changes is based on an adaptive principle, in which a "memory" is modeled by vector w of memory coefficients, which are connected with the particular coordinates of the space of the Lagrangean multipliers. The generation of r j is done using equation r j ϭ w ϩ x, where x is obtained as a realization of |I|-th dimensional random variable with a uniform probability distribution of their components over interval ϽϪ2A, 2A Ͼ.
The vector of the memory coefficients is initialized to zero values at the beginning of the optimization process, and whenever new point u kj is constituted and its objective function value F(u kj ) is found, the vector of memory coefficients is updated in accordance with the rule w i ϭ Bw i ϩ C.(F(u kj ) Ϫ F(u k ))(u i kj Ϫ u i k ), where B is the parameter of oblivion and C is the parameter of learning intensity. Both parameters were chosen from interval (0, 1).
The second phase utilizes the statistical gradient enumerated in the first phase and explores the neighborhood of the current point u k in the associated direction. The statistical gradient is normalized to the unit length and the first experiment is undertaken with point u ϭ u k ϩ4.A.r/ͰrͰ. Then the step length is diminished and the new point is tested. This searching process is repeated until the number of failed of attempts to improve solution reaches the prescribed maximal value.
The best-found solution in both phases represents the resulting move to u kϩ1 subject to condition that the new point has a better function value than the current one.
The whole searching process of the method SUPRA is repeated until one of the following conditions is met. The first condition tests the equality of given number N and the total number of move designs and the second one tests the equality of given number Nb and the number of move designs performed from the last improvement of the best-found solution.

Numerical experiments
To test and compare both approaches, the associated algorithms were implemented using Delphi 7 programming environment. To perform the numerical experiments, Pentium 4, 3.0 GHz, 1 GB was used. The approach to the capacitated facility location problem was tested with data originating at the Slovak road network with 2907 dwelling places, where each of them represents one customer. Seventy-one centres of the former districts formed the set of the possible facility locations. The testing problems are referred as LAS*071xy2911.txt, where * denotes D, B, A, F, x ϭ 3, …, 9 and y ϭ 2, …, 5. Then the number of the test problems is 28. The set of the test problems is partitioned into four groups in accordance with different values of y, which corresponds to the ratio cost coefficients used in enumeration of c ij .
In the performed experiments the following attributes of the resulting solutions are observed and evaluated in Tables 1 -4. The first attribute of the solution quality is resulting lower bound (14) of the objective function value of an optimal solution of the original problem (LB), but it must be taken into consideration together with some measure of capacity constraint violation, because this way of handling with capacity constraint admits some violation. Two parameters were chosen to evaluate this possible capacity constraint violation and they are considered to be the second and third quality characteristics of a solution. The first parameter (Max Over.) is ratio of the maximum value of positive differences between the workloads of facilities and their capacities and the associated capacity. The second parameter (Sum Over.) is ratio of the sum of positive differences between the workloads of facilities and their capacities and the sum of all located facility capacities. The computational time in seconds is reported in column (Time).

Conclusions
As concerns lower bound of the original problem, the abovedescribed approaches to the capacitated facility location problem Comparing the lower bounds provided by both approaches on the groups of tested instances has shown that there is no unambiguous winner.
When applying the approaches on the first two groups of instances, the SUPRA wins in each case. As concerns the third group, the SUPRA wins only in half of the instances and the SUPRA completely fails in the fourth group.
This effect can be caused by bias setting of the parameters A, B, C of the adaptive process of SUPRA. The excellent performance of the SUPRA considering facility overload minimization is worth noticing. In most of the solved problems, SUPRA achieves half of the maximal overload produced by the sub-gradient method.
The negligible computational time of the sub-gradient method in comparison with the computational time of SUPRA evokes an idea for the future research, where composed heuristic should be designed and studied. This heuristic could combine both abovementioned approaches so that it could use the sub -gradient method in the first phase and then utilize the SUPRA to find a new improving direction. The group location problems for y ϭ 2 Table 1 Problem The group location problems for y ϭ 3 The group location problems for y ϭ 4 Table 3 Problem The group location problems for y ϭ 5 Table 4