An evolutionary algorithm for joint bi-criteria location-scheduling problem

A new case of joint location and scheduling (ScheLoc) problem is considered. It deals with selecting a non-fixed number of locations for identical parallel executors (machines) from a given set of available sites. Simultaneously, a schedule for a set of tasks is sought. For every task, it comprises an executor carrying-out the task and the moment of time when the performance of the task is started. The locations for executors and the schedule are evaluated by two criteria: the sum of task completion times and investment costs incurred when locations for executors are selected and launched. It is justified that the joint optimization problem is strongly NP-hard. In consequence, a heuristic algorithm Alg_BC is proposed, which uses the general scheme of NSGA II provided for the multi-criteria optimization. The performance of Alg_BC is evaluated for small instances by exact solutions determined by the Matlab solver. The sensitivity analysis for bigger instances is also provided, which among others, allows examining the influence of both component criteria on results generated by the evaluated algorithm. A case study dealing with the evacuation of citizen groups from danger zones is provided as an example of the investigated bi-criteria ScheLoc problem. The usefulness of Alg_BC is confirmed as well.


Introduction
Investigation of topics of complex optimization or decision making problems can be observed in operations research and discrete optimization. In such problems, some interconnected sub-problems can be distinguished. Neglecting the interconnections among sub-problems is a simple paradigm constituting solution algorithms of the mentioned problems frequently met in the initial phase of their investigation. Advances in the methodology and computing potential of solution algorithms facilitate dealing with such problems as a whole, and, as a consequence, it is possible to obtain improved solutions. Let us mention project scheduling-location, replenishment-location, location-routing, production-inventory-location, routingscheduling as examples of complex problems, where sub-problems are solved jointly, (e.g., Rostami & Bagherpour, 2017;Pasandideh et al., 2018;Prodhon & Prins, 2014;Shahabi et al., 2018;Paraskevopoulos et al., 2017). This paper's investigations are addressed at another complex problem called location-scheduling. Scheduling of a given set of tasks and location planning of executors (machines) as performers of the tasks are interconnected parts (sub-problems) of a complex problem. It is assumed that a finite number of tasks are deployed at given original locations in a planar area. Every task has to be moved from its original location to the position of the corresponding executor site, which is not known in advance. Then, all the tasks moved to the same executor are scheduled.
Consequently, the purpose of the complex problem is to jointly deploy a set of executors in a planar area and schedule a set of tasks on previously located executors such that relevant criteria evaluating both activities are optimized. There are many applications described in the literature where both location and scheduling are required. The unloading of ships by cranes is the first example (Scholz, 2012). Unloading ships waiting at sea can be interpreted as tasks. Locations for cranes as executors at an embankment need the determination together with the order of ships' service after their reaching the embankment. Kalsch and Drezner (2010) invoke a similar application where one has to find locations for ships on a berth and schedule the loading process of containers. The operation of movable crushers in the mining industry described in Kalsch (2009) is another example. The evacuation of groups of people from danger zones can be indicated as the next example (Heßler, 2017). Each group can be treated as a task performed by a salvage crew being an executor. Locations of the so-called concentration points where groups of people have to gather are the subject of a decision. Simultaneously a schedule of the evacuation is needed. Another example concerns modern computer distributed systems where a given number of customer tasks wait for a computation service provided by selected servers. Such servers have to be previously selected from the given set of available servers, and, then, they are assigned in time to perform the customer tasks (Saldana and Suznjevic, 2017).

Related work
Such a introduced complex problem referred to as ScheLoc has been considered for the first time in Hennes and Hamacher (2002). Then, it has been discussed and developed in many works. Many particular cases have been taking into consideration depending on the kind of location area for executors, a type and number of executors, the criterion for evaluating a schedule of tasks as well as used solution algorithms. A summary of the considered cases is provided in Table 1, see also Liu et al. (2019). A majority of works deal with a discrete area for the deployment of executors where a finite set of available positions for executors is known and given a priori (e.g., Heßler and Deghdak, 2017;Liu et al., 2019;Ławrynowicz and Józefczyk, 2019). A non-empty subset of this set is to be selected. In many works, a one-element subset is assumed. As a consequence, only the location of a single executor is sought. Two task scheduling criteria, i.e., the makespan max C and the sum of completion times j C  (Pinedo, 2012), have been noticed in the literature as the evaluation of a task schedule. The criteria serve in the analyzed works as the assessment of ScheLoc as a whole. Other possibly used criteria assume in some works the supporting role for the deployment of executors, (e.g., Kalsch and Drezner, 2010;Ławrynowicz and Józefczyk, 2019). Authors proposed different heuristic algorithms due to a severe time complexity not only of the joint problem but also of its component sub-problems. Let us remember that the task scheduling sub-problem with different release dates is strongly NP-hard for a single executor and criterion j C  (Lenstra et al., 1977). Corresponding sub-problems for at least two executors and criterion max C are NP-hard (Pinedo, 2012). The other sub-problem is also NP-hard when the deployment of executors is treated as a particular case of the uncapacitated facility location problem (UFLP) or p-median problem (e.g., Krarup and Pruzan, 1983).
As has been mentioned, ScheLoc has been presented for the first time in Hennes and Hamacher (2002), where a straightforward case is discussed. Available positions for a single executor are limited to a graph. An integrated model is provided in which executor locations define release dates for tasks. The polynomial solution algorithm is presented, comprising locational and scheduling parts. Release dates are calculated in the former part based on the distance among positions of tasks and the position of the executor. Then, the scheduling part is solved by the ERD rule.
Two similar works of Elvikis et al. (2007) and Elvikis et al. (2009) focus on the same task scheduling problem max 1| ( ) | j r x C where ( ) j r x is the release date dependent on the executor's position x. The locational analysis is substantially extended in comparison with previous works. In particular, the continuous case is considered, and various distance functions are studied. The ERD rule is used as well, which uses the values of release dates calculated based on regions in the available area where the dates do not change. The authors present three polynomial algorithms. One constructive algorithm uses geometric properties. Two other propositions contain linear programming models that are solved by the Matlab solver.
The idea of calculating release dates for task scheduling is continued in Kalsch and Drezner (2010). The sum of completion times of tasks j C  is considered apart from max C but again for a single executor. The solution algorithm of the scheduling part is based on the ERD rule, which requires equal task processing times when the criterion j C  is minimized. The core investigations are concentrated on the locational analysis for nonconvex cases with the use of the big triangle small triangle approach (Drezner and Suzuki, 2004). The version with many executors has been introduced in Rajabzadeh et al. (2016), where both discrete and continuous planar areas are considered. For the criterion max C , the mathematical models are presented, and the solver GAMS is employed as a solution tool. The authors provide examples of results limited to ten tasks and three executors.
The task scheduling on many parallel identical executors, the makespan max C as the criterion, and discrete locations for executors are investigated in Heßler andDeghdak (2017), andHeßler (2017). The constructive algorithm has been proposed that is based on the iterations between location and assignment sub-problems. In the first step, the algorithm starts with a single cluster containing all jobs, and, for this cluster, a single optimal machine location is computed. In the second step, a set of jobs is removed from the existing cluster, and a new optimal location is set. The procedure continues until it reaches the number of demanded locations.
The Master Thesis of Piasecki (2018) and the following conference paper of Piasecki and Józefczyk (2018) extend the investigations on the minimization of j C  with the development of many parallel identical executors in the continuous area. The evolutionary algorithm solves the resultant strongly NP-hard problem. The results for small instances are compared with optimal outcomes generated by the Matlab solver.
A more advanced complex problem is discussed in Ławrynowicz and Józefczyk (2019). Namely, many unrelated parallel executors should be located in a discrete area, and apart from max C , the additional criterion has been introduced in the form of the sum of release dates for tasks. The joint problem is considered as the composition of two sub-problems, i.e., the location of executors (LE) and the scheduling of tasks (ST). Two approaches are proposed, solved and compared. The first one, referred to as the sequential approach, comprises two steps: solving the LE sub-problem to minimize the sum of release dates and ST to minimize max C with the release dates calculated based on the previously solved LE sub-problem. The known algorithms have been employed to solve both sub-problems. In the second approach called the systemic approach, both sub-problems are solved jointly to minimize max C . Executor locations are determined to minimize max C not the sum of release dates, like in the sequential approach. The memetic algorithm based on Tabu Search metaheuristics has been proposed as the solution tool. The conducted experimental comparison confirmed the substantial advantage of the systemic approach in comparison with the sequential one.
It is also worth pointing out works expanding the topic directly connected with ScheLoc. Liu et al. (2019) consider the stochastic task execution times and the expected value of j C  together with the total location costs of executors as two criteria. However, the authors apply a simple scalarization, and, in consequence, a single criterion optimization problem is solved. A genetic algorithm and the sample average approximation algorithm were elaborated for the large and small instances, respectively. In turn, four criteria are taken into consideration in Wesolkowski et al. (2014) while assessing locational and assignment decisions for military applications. The solution algorithm is based on the NSGA II approach for solving the resulting multi-objective optimization problem.

Contribution of Paper
In comparison with previous works, the contribution of this paper can be explained as follows.
1. First of all, a new optimization problem is solved. A single criterion has evaluated all variables (decisions) in previous works. The distances and release dates, which have been used in some works for the determination of executor locations, were auxiliary criteria. The makespan max C , or the sum of completion times j C  , was the main criterion, and it assessed the schedule of tasks. The quality of the deployment of executors was only indirectly evaluated by max C or j C  as they were the basis for the calculation of release dates. We propose to consider the criterion which assesses the cost of executor locations apart from the criterion j C  evaluating the schedule of tasks. In consequence, a new bi-criteria joint locationscheduling problem is directly solved. Consequently, we search for a Pareto front of solutions, not for a single solution like, e.g., in Liu et al. (2019), where the simplification via scalarization is applied, and a unique solution is determined.

2.
In the face of strong NP-hardness of the optimization problem, the time-efficient heuristic algorithm is proposed. It is based on the NSGA II scheme (Deb, 2001).
3. The usefulness of the algorithm has been experimentally confirmed via comparison with the Matlab solver outcomes or sensitivity analysis for small or big instances, respectively. We propose bespoke performance indices for the evaluation, being the adaptation of known indices, which allows for a better assessment of resulted Pareto fronts.

4.
A case study concerning the evaluation has been provided, which illustrates the application potential of the problem.
The remainder of the paper is organized as follows. The mathematical model is provided in Section 2, which is followed by the presentation in Section 3 of the evolutionary algorithm. Section 4 is devoted to the computational experiments, which, first of all, allowed us to evaluate the solutions determined by the proposed algorithm by the exact solutions given by the Matlab solver. A case study presenting the application of the algorithm is provided in Section 5. Final remarks complete the paper.

[ , ]
j j j a a a  is the location of the task j. The current task j is characterized by execution time j p , ready time j  , as well as transportation speed j v , and it needs to be performed by a single executor (machine) selected from a set of identical machines referred hereafter to as executors. Executors' prospective locations should be selected from a set We assume that the number m, 1 m    of employed executors, which is equal to the number of selected locations from set B, is not known a priori.
Let us introduce a vector x is equal to 1(0) if the jth task is scheduled on the ith executor as the kth (otherwise). Performance of task j by executor i follows its transportation at a distance So, the task j is composed of two parts: transportation from the initial location j a and the essential performance at the selected executor's location bi. The first part needs the time interpreted as a release date for the task j , which has to be performed once by one executor. The time j p is allotted for the execution of the second part. We employ two conflict criteria for the evaluation of decisions x and y . The first one as the sum of completion times of all the tasks evaluates the schedule x , but it certainly depends on the locational decisions y : where the auxiliary variable ( , ) ik C x y stands for the completion time of a job performed by the ith executor as the kth and depends on decision variables via constraints (6) and (7) (Kooli and Serairi, 2014). The total location (investment) cost of all locations used by executors serves as the second criterion: where i c is the location cost of i b . Obviously, launching fewer locations m decreases the investment cost (2) ( ) q y , but, simultaneously, it increases the sum of completion times (1) ( , ) q x y , and vice versa. Therefore, the trade-off must be sought, which leads to the bi-criteria problem. The following constraints imposed on decision variables x and y enable having feasible solutions: , , 1 1 ( ) 0, 1, 2,..., , 1, 2,..., 1, , {0,1}, , 1, 2,..., , 1, 2,..., .
i jik Constraints (4) and (5) ensure that each task is performed on one position of a single executor, and no more than one task can be performed on the launched location of each executor, respectively. Moreover, in a feasible solution, the kth task in order on the ith executor cannot start before its release date or the completion time of a task which is scheduled on the (k-1)th position of the same executor. Constraints (6) and (7) fulfil this condition. The next constraint guarantees that the number of launched locations equal to the number of used executors is between 1 and  . According to (9), all tasks assigned to an individual executor have to be represented by a single sequence of '1's. It means that any sub-sequence of zeros within the whole sequence can occur. In consequence, the number of equivalent solutions represented by matrix x is limited. Finally, constraint (10) defines the domains for decision variables. In consequence, the following bi-criteria optimization problem, find the schedule of tasks x and the locations of executors y minimizing a vector of criteria (4) Proof Let us consider a special case of BC_ScheLoc when 1 m    , and values of y is known and fixed. Then (2) ( ) q y is constant, and (1) ( , ) q x y along with constraints (4)-(10) can be simplified to {0,1}, , 1, 2,..., and the only optimization variable   1 , 1, , 1, is the sub-matrix of x . Please note that the task release dates j r do not depend on x for 1 with constraints (12)- (17) is the well-known task scheduling problem 1| | j j r C  , which is strongly NP-hard (Lenstra et al., 1977).  Such a high computational complexity of the bi-criteria optimization problem justifies the search for heuristic algorithms that allow obtaining solutions in an acceptable time. Therefore, an evolutionary approach has been adopted. The general scheme of the NSGA II (Deb et al., 2002) with some improvements have been used and presented in the next section.

Solution algorithm
As the problem includes two contrary criteria, a multi-objective optimization approach is applied.

Encoding
The adopted method searches for solutions iteratively and independently for each subproblem that takes into consideration the fixed number of executors. It starts with the generation of an initial population i The solutions are created by uniformly setting the values within the genome. The higher than ϑ Hamming distance between each pair of chromosomes is assumed. This diversification strategy helps to avoid premature convergence during the initial iterations.

Ensure: Set of non-dominated solutions PF
S and the Pareto front PF.

Selection, Crossover, and Mutation
The presented method preserves the NSGA II structure. A selection process is based on the crowded comparison operator described in Deb (2001). The chromosomes undergo a parallel crossover with the use of two different operators with the probability 0.95 c   (Grefenstette, 1986). For the binary parts, the Count-Preserving Crossover (CPC-2) operator has been used (Hou and Chang, 2004). This approach allows preserving  (Syswerda, 1991). The procedure is repeated until it achieves   solutions, which are passed on for further processing.
After the crossover, the mutation is applied to each offspring individually with a small probability 0.01 m   (Hesser and Männer, 1990). Due to the importance of the order of tasks performed on a single executor, the Simple Inversion Mutation operator (SIM) has been applied for the integer part of E (Holland 1975). The SIM operator selects randomly two cut points in the sequence, and it reverses the subsequences between them. Additionally, the binary parts in E are randomly modified by the Swap Mutation (SM) operator (Oliver et al., 1987). The randomly selected single gene is swapped with another one that contains the different value to preserve the constant number of bits equal to one.
The crossover and mutation operators produce  new offsprings that are stored in       1, an index of iteration. The used chromosome's form does not need to perform any repair strategies after the modifications made within E . The achieved results are aggregated in the long-term memory where   L iz is the size of Pareto front in the izth iteration if and only if they are nondominated in the Pareto's sense. For all random-based operations, the Well Equidistributed Long-period Linear pseudorandom number generator named WELL44497a was applied (Panneton et al., 2006).

Stop condition
Finally, the stop condition is fulfilled if there is no improvement in solutions representing the Pareto front through  consecutive iterations. Moreover, the parameter max  has been introduced to restrict the maximal number of iterations for each subproblem. Consequently, the NSGA II-based heuristic solution algorithm referred to as Alg_BC is proposed for solving BC_ScheLoc.

Tunning
The tuning is a two-step procedure, and the values of different parameters are fixed before the execution of the heuristics. In the first step, the algorithm's parameters are divided into two separate groups: chosen for tuning and non-tuned. During the second step, the grid search method is applied for the selection of the best set of parameters.

Computational experiments
The evaluation of Alg_BC is the most significant purpose of the conducted computational experiments. First of all, it consists in the comparison of Pareto fronts generated by Alg_BC with similar outcomes created by the Matlab solver for small instances. The results are provided in Subsection 4.2, which follows the foundations for all experiments specifying the used dataset and quality indices adequate for evaluation and comparison. The sensitivity analysis completes the experiments. It allows checking the influence of the critical parameters on the results delivered by Alg_BC. All the computations have been made using a PC with the Intel Core Xeon W CPU 14-core processor of 2.5GHz equipped with 128GB of RAM. The algorithm Alg_BC has been implemented in the R language.

Foundations of computational experiments
between (0,0) and q  , we proposed a hypervolume-based quality index AH I (Fig. 2), which is the adjustment of a wellknown hypervolume indicator H I (Zitzler et al., 2007;Laszczyk and Myszkowski, 2019 where (  For the majority of analyzed instances, Alg_BC was able to find during its 25 runs the same point q  as the solver, Table 2. The average values of (1) q  , (2) q  are worse than (1) T, T, q  no more than 17%, and 54%, respectively. The worst outcomes do not exceed 26% and 62% for (1) T, q  and (2) T, q  , respectively. The worst results for the instances (15, 8), (25,4) show the weaker performance of Alg_BC for a comparatively big number of possible locations  or when  is small compared with the number of tasks n . Undoubtedly, the small values of  support the determination of cost criterion values (2) (2) T, For such instances, executors have been deployed at all available locations, and only such a decision has allowed reaching the best values of (1) q  .

Table 2
Dependence of q  and T, q  on n and  (bold type points the results of Alg_BC equivalent to those obtained by the solver) The values of the akin single-valued index Dist are presented in Table 3. The conclusions are similar. Again, Alg_BC was able to find during its 25 runs the results of a similar quality as the solver. The difference is worse by 9% on average but no more than 13%. However, the performance of Alg_BC deteriorates as n and  grow. The hypervolume-based index AH I , which values are presented in  Table 5 are incomparable as the former ones reach dozens of hours for launched small instances, unlike the latter ones, which do not exceed several seconds. Due to the exponential time complexity, the usage of the solver is impossible for higher instances, even if the generated results are non-operational but tactical or strategic decisions where the computational time is not the most essential. n  as in Tables 2-5. The area between both fronts is the zone where we can expect to have PFs generated by Alg_BC. The size of this zone is rather small, which indicates the decent performance of Alg_BC. This size, however, is more significant for the higher number of tasks n . Moreover, it is worth noting that the size of the fronts grows when increasing the number of available executors  . We have also examined the performance of Alg_BC for bigger instances from the used dataset for various values of parameters  and  . The results are provided in Tables 6, 7 linearly when increasing  . The comparison of these times with adequate results from

Case study
The combined location and task scheduling have some real-world applications mentioned in Section 1. The evacuation of people from danger zones is one of them (Heßler, 2017). However, it is the subject of broader interest and has been investigated in many other works, e.g., in Sbayti and Mahmassani (2006), see also Osman and Ram (2017), Amideo and Scaparra (2017), Abounacer et al. (2014) for other cases and methods. Referring to these and similar works, we demonstrate in this section the advantage of modelling a specific planning problem of evacuation in the framework of the considered problem and the usefulness of applying Alg_BC for the determination of evacuation plans. The case study deals with the planning of evacuation of citizens living in the Polish medium-sized city Trzebnica located in the south-west part of Poland. Some local requirements for evacuation planning specify the usage of helicopters as transportation facilities. Therefore, it is necessary to propose locations for adequate landing zones together with the assignment of home estates to the zones along with the order of evacuation. The possible locations of landing zones are shown in Fig. 6. The purpose of planning is not only to provide the solution that ensures the prompt leaving of the danger area but also to limit investment costs. Both conditions are fulfilled, and all decisions can be made when modeling the problem as joint location and scheduling. The city has been divided into thirty areas ( 30 n  ). Moreover, eight prospective sites for the landing of helicopters have been distinguished ( 8   ). When establishing the data, we assumed that the execution times j p are proportional to the size of the population from respective areas, and each person needs 1 min to be served. Then the ready times have been calculated as 3ln j j p   . The values of other problem parameters adequate to the case study are presented in Table 8.  100, 200, 200, 300, 300, 600, 400, 800, 300, 150, 600, 200, 200, 300, 300, 100, 400, 600, 100, 100, 200, 500, 150, 300, 200, 50, 400, 300, 500, 300. ( , ) 13. 8, 15.9, 15.9, 17.1, 17.1, 19.2, 18.0, 20.1, 17.1, 15.0, 19.2, 15.9, 15.9, 17.1, 17.1, 13.8, 18.0, 19.2, 13.8, 13.8, 15.9, 18.6, 15.0, 17.1, 15.9, 11.7, 18.0, 17.1, 18.6, 17 Table 10 and Fig. 7 demonstrate the coordinates of sixteen PF points representing sixteen equivalent solutions, which details are provided in Table 9. For the readability, the sequences of tasks for individual executors (launched locations) substitute matrix x . The last row of Table 9 informs us how many times the individual locations i b occur in PF.
The variety of Alg_BC's outcomes facilitates a planner/a user to choose an appropriate solution.
In the case of an emergency, a planner is forced to minimize the evacuation time (1) q at the cost of (2) q . In consequence, he/she would be ready to accept 7 out of 8 locations for landing zones ( 7 m  ).
However, a planner is usually unable to accept that many landing zones due to a limited budget. The analysis of PF provides information on which landing zones are suggested the most in the solutions. Namely, Zone #1 ( 1 b ) is most often proposed.
Zone #4 ( 4 b ) appears in the results in the second place.
Zones #2, #7, #8 ( 2 , b 7 , b 8 b ) have the next similar numbers of occurrences in PF. However, taking into consideration the geographical distribution of landing zones, a planner should choose Zone #8 ( 8 b ) as the third landing zone. It is doubtful that more landing zones would be built for such a small city. It is worth noting that the landing zones could be launched by turns in a more extended period starting from Zone #1. The versions with one, two, or three zones are illustrated in Table 11, which like Table 9 presents the sequences of tasks.

Final remarks
We investigate a bi-criteria integrated approach for solving the location and scheduling (ScheLoc). The case is considered with a non-fixed number of identical parallel executors (machines) as well as non-preemptive independent tasks.
Two criteria evaluating both a schedule of tasks and location costs have been used. To our best knowledge, this is the first work elaborated for the bi-criteria ScheLoc problem with the sum of task completion times as the scheduling criterion. We provide a heuristic evolutionary algorithm due to the strong NP-hardness of the resulting combinatorial optimization problem. The conducted computational experiments showed the usefulness of the proposed algorithm. In many cases, the results generated by the algorithm are close and consistent with the exact outcomes provided by the Matlab solver. The average deterioration of the results does not exceed 19% when using the most comprehensive hypervolume-based quality index.