Capacitated Lot-Sizing Problem with Sequence-Dependent Setup, Setup Carryover and Setup Crossover

: Since setup operations have signiﬁcant impacts on production environments, the capacitated lot-sizing problem considering arbitrary length of setup times helps to develop ﬂexible and e ﬃ cient production plans. This study discusses a capacitated lot-sizing problem with sequence-dependent setup, setup carryover and setup crossover. A new mixed integer programming formulation is proposed. The formulation is based on three building blocks: the facility location extended formulation; the setup variables with indices for the starting and the completion time periods; and exponential number of generalized subtour elimination constraints (GSECs). A separation routine is adopted to generate the violated GSECs. Computational experiments show that the proposed formulation outperforms models from the literature.


Background
Master production plans of a main production process of a supply chain provide all the participants with the visibility regarding the operation. Production planning is referred to as "lot-sizing" in a discrete production environment [1]. The capacitated lot-sizing problem (CLSP) has been considered in various industries to generate production plans reflecting production restriction [1,2].
Generally, lot-sizing problem assumes that planning horizon is finite and can be divided into multiple time periods, which are terms as "time buckets". The time buckets are classified into smalland big-time buckets, when the comparison of the relative length of the time period is made with respect to the relative length of the production lot [3]. Models considering small-time buckets allow at most one or two products to be produced in a period and models considering big-time buckets allow more than two products to be produced in a period [1].
Setup operations prepare processing units (e.g., materials and machines) to manufacture production lots. In process industries, the setup time required for replacing production models accounts for a large proportion of capacity utilization. Moreover, variations in the setup time may be highly depending on the product models and their operating sequence. The setup time between similar products is relatively short, allowing engineers to setup multiple times in a day or a shift. However, setup operations requiring more than one day or several days may be necessary between products when the processing material or equipment needs to be changed. For solving the big-time bucket CLSP model, several modeling techniques have been developed to reflect these characteristics of setup time. Sequence dependent setup time model is used to reflect setup time that depends on product order.
Several studies [4,5] have introduced CLSP models, in which the setup state can be preserved across several periods and have used various terms such as "linked lot sizes" and "setup carryover". In this paper, the concept is called as setup carryover. The setup carryover allows a setup state to be maintained from one period to the following one-namely, a setup carryover model enables a product to be produced in the following buckets without additional setup time. In this concept, production planers in process industries often prefer to produce only a single product without model change during several days or shifts.
In some manufacturing environments, the setup started in period t can be split between the periods t and t + 1 [6,7]. In this paper, the concept is called as "setup-splitting".
In most continuous manufacturing industries such as process industries, due to the presence of long setup times, manufacturers allow the setup to be carried across more than one period [1], the concept is called as "setup crossover". Setup crossover enables the incomplete setup operation to cross over the boundaries of time periods. In the case that setup times are longer than a period length, the setup operation may be performed in more than two periods. However, if all the setup times are less than a period, the setup crossover is exactly same as the setup-splitting.

Literature Review
There are several studies published on big-bucket CLSP with setup times. However, studies on CLSPs with setup carryover are limited. Sox and Gao [4] presented two formulations for the CLSP with setup carryover. Suerie and Stadtler [5] presented a formulation for the CLSP with setup carryover and their computational tests have shown that their formulation is better than the formulation provided by [4]. Further, Menezes et al. [6] proposed CLSPs with setup-splitting considering sequence-dependent and non-triangular setups. In the formulation, the classical subtour elimination constraints, known as the Miller-Tucker-Zemlin (MTZ) constraints, introduced by Dantzig et al. [8] are used. Kopanos et al. [9] introduced a CLSPs with setup carryover and setup-splitting problem with items classified into product families. They considered sequence-dependent setups between products in different families and sequence-independent setups between products in the same family. Mohan et al. [10] extended the formulation of Suerie [5] to address setup-splitting. Camargo et al. [11] presented three formulations for the two-stage lot-sizing and scheduling problem. With computational tests, they showed that the formulation with setup-splitting is the most flexible to incorporate setup-related features. Ramya et al. [7] proposed a formulation for the CLSP with setup carryover and setup-splitting and no backorders or lost sales. Fiorotto et al. [12] compared two formulations proposed by Mohan [10] and by Menezes [6] and proposed symmetry braking constraints. With computational experiments, they showed that the proposed constraints are effective for the formulation proposed by Mohan.
Sung and Maravelias [13], Belo-Filho et al. [14] considered the CLSPs with setup carryover and setup crossover with the setup operations can be split into more than two periods. Sung and Maravelias [13] proposed a big-bucket formulation. Belo-Filho et al. [14] proposed two formulations for the case with backlog. One of the formulations involves a time index disaggregation, defining the setup variables with indices for the starting and the completion time periods of the setup operation.
For the case of sequence dependent setup, two groups of studies have been published. First group of studies consider the problems where at most one lot could be generated per a product for each time bucket. The other group of studies allow multiple lots per a product for each time bucket. The former problems usually consider triangular setup times and costs, while the latter problems usually consider non-triangular setup times or minimum lot sizes. Almada-Lobo et al. [15] presented two models for the CLSP with sequence-dependent and triangular setup times and costs using the MTZ subtour prohibition constraints. Clark et al. [16] formulated a sequencing and lot-sizing model with non-triangular setup times based on the asymmetric traveling salesman problem. Ramya et al. [1] presented two models based on MTZ subtour prohibition constraints for the CLSP with sequence-dependent setup times and setup cost and with setup carryover and setup crossover. In the model, no idle times are allowed during setup operation and production of a product. In the models of [1,15,16] at most one lot per product can be produced in a time period.
Menezes et al. [6] allowed production of multiple lots per period. Clark et al. [17] presented a stronger formulation than Menezes et al. [6] for modeling the production of multiple many a product per period by using a polynomial number of multicommodity-flow-type constraints [18]. Mahdieh et al. [19] presented formulations based on Clark's formulation for lot-sizing and scheduling with minimum lot-size, non-triangular sequence-dependent setup times and costs with setup carryover and setup-splitting.

Contribution and Organization
In this paper, a CLSP with sequence-dependent and triangular setup times and with setup carryover and setup crossover is considered. For the problem, a new mathematical formulation is proposed. The formulation is based on the facility location extended formulation [2]. Moreover, it uses the setup variables with indices for the starting and the completion time periods as Belo-Filho's formulation [14] and exponential number of generalized subtour elimination constraints (GSECs) to prevent subtours in a time bucket.
The mixed integer programming (MIP) formulations are usually compared in terms tightness or compactness [20]. The tightness of an MIP formulation is defined as the difference between the optimal values for the continuous relaxation and the original MIP problem. Tightening an MIP formulation, usually by adding cutting planes or valid inequalities, can reduce the search space that the solver needs to explore [21]. Meanwhile, the compactness refers to the quantity of data that must be processed when solving the problem [20]. Because there are many relaxations that must be solved when solving an MIP problem, the compactness of such a problem always refers to the size of these relaxations. A more compact formulation can speed up the search for the optimal solutions. The tightness and the compactness of the formulation are shown by comparison with one of Mahdieh's formulation (deferred to as MLOV-SM) [19] and one of Ramya's formulation (deferred to as MM2) [1].
The remainder of the paper is organized as follows. In Section 2, a formal description of the problem is given. In Section 3, the formulation for the problem is presented. The formulation has exponential numbers of the GSECs. The separation algorithms for the GSECs are explained in Section 4. In addition to comparison with other models, computational experiments and their results are reported in Section 5. Finally, Section 6 concludes this study and suggests some directions for further research.
During the computational experiments, the author found that the setup-splitting is not implemented properly in MLOV-SM. The modified and fixed version of MLOV-SM is reported in Appendix A.

Problem Statement
In the following, a formal description of the problem is presented with assumptions for clarifying and simplifying the mathematical model.
A planning horizon consists of finite number of time buckets. For each bucket, an index ranging from 1 to n is assigned; that is, the number of buckets is denoted by n. For each bucket t = 1, · · · , n, a positive amount of capacity denoted by C t is assigned. In the problem, the production lot sizes and schedules for a set of items, denoted by I, are planned. For each item i in I and bucket t, a non-negative demand quantity denoted by d i t is given. The demand could be satisfied through production before the due date with inventory holding cost h i t for each bucket or after the due date with backlogging cost b i t . Production of one unit of items i in bucket t consumes a positive amount of capacity denoted by p i t , namely the process time per unit. All the demands must be met until the end of the planning horizon. Let T i denote the set of buckets where a positive amount of demands is defined for item i. Considering sequence-dependent setup time and cost, production sequences were decided for each bucket. Production change from item i to j requires setup time τ ij and causes cost c ij . If capacity allows, multiple setup could occur in a bucket. With long setup time or inadequate capacities, setup operation could be performed during consecutive multiple buckets. Production of the last lot could resume in the following bucket without any burden of setup operations. The initial and final setup states are not given. The objective is to minimize the overall cost, which includes backlogging, inventory holding and setup costs.
Under the above descriptions, we can introduce the following two assumptions without loss of generality.

Assumption 1.
No idle time is allowed during setup operations.
Bucket-independent setup times and costs were considered, therefore, any idle during setup operations could be moved to the buckets before setup started or after setup finished without increasing the cost. Therefore, the Assumption 1 is acceptable.
In a bucket, at most one lot can be generated for an item.
In [14], the setup variables with indices for the starting and completion time periods of the setup operation were introduced to solve CLSP with setup carryover and setup crossover. According to Assumption 1, for each type of setup (for each item in [14]) with setup time τ, a setup variable can be defined for each ordered bucket pair (s, t) satisfying the following constraint: Proposition 1. The number of ordered bucket pairs satisfying constraints (1) are limited by 2n − 1. (Proof: refer to [14]).
For a sequence-dependent setup operation from item i to j, we can consider the setup variable indexed by bucket pairs in the following set:

Formulation
The problem formulation is based on three building blocks: the facility location extended formulation [2]; the setup variables with indices for the starting and the completion time buckets [14]; and exponential number of generalized subtour elimination constraints (GSECs) [18].
The following two variables come from the facility location extended formulation: setup state variable, y i t , for item i and bucket t; and production portion variable, x i tu , for item i and bucket t and bucket u ∈ T i (having positive demand d i u ). y i t has value one if setup states for item i in bucket t is on and item i can be produced and 0 otherwise. Moreover, x i tu denotes the proportion of the demand d i u produced in bucket t. For example, in the case d i .5) means that total demand of d i 2 and half of demand d i 3 are produced in bucket 1 and the remaining half of demand d i 3 is produced in bucket 2. In the formulation, model change or setup operation is indicated by the following setup-indicating variable: e ij st , for a pair of items i j ∈ I and an ordered pair of buckets (s, t) ∈ T ij . The variable has a value of one if setup operation from item i to j starts in bucket s and finishes in bucket t and the value is zero otherwise. If s < t, e ij st indicates setup crossover. When s + 1 < t, e ij st indicates setup crossover longer than two buckets. By Proposition 1, the number of the variables (e ij st ) is bounded by (2n − 1) × |I| × (|I| − 1). R i t denotes the remaining setup times for item i at the beginning of bucket t in the case of setup crossover. In the formulation, a setup carryover indicating variable, z i t , is used for each item i and bucket t = 1, · · · , n + 1. z i t is has value of one if the setup state for item i is maintained at the boundary between bucket t − 1 and t and zero otherwise. As mentioned in the problem statement, the initial and finial setup states are not given. Hence, z i 1 's decide initial setup state and z i n+1 's decide final setup state.
Lot sequencing with GSECs needs exponential number of constraints, but do not need any additional variables. The process of deriving the constrains is introduced at the end of this section. Table 1 presents the parameters and the variables used in the formulation. Table 1. Notation.

Types Notation Meaning
Set The formulation M1 is as follows: u∈T i i,j∈I|i j (t,u)∈T ij Processes 2020, 8, 785 Here, the objective function (3) minimizes the total cost including backlogging, inventory holding and setup cost. Constraints (4) ensure that each demand is met within the planning horizon. Constraints (5) and (6) ensure that production occurs only in the bucket where a lot is generated. Constraints (7) ensure that capacity consumption during production and setup change plus the idle time equal to the capacity of each bucket, considering the remaining setup times during setup crossover. Constraints (8) ensure that setup state is dedicated to at most one item at the beginning of a bucket. Constraints (9) ensure that a lot comprising one item is generated in a bucket when the setup of the item is carried over or a setup operation to the item ends in the bucket. Constraints (10) ensure that one lot is followed by a setup carryover operation to the next bucket or a setup change operation to another item. Constraints (11) ensure that remaining setup times are allowed only during setup crossover. From Assumption 1, idle time is not allowed during setup operation, this is reflected in constraints (12). First two terms of the left-hand side are added to lift the constraint. Constraints (13) prevent the subtours in each bucket. Constraints (14)- (19) describe the domains of the variables. From Assumption 2, the number of lots is limited to one for each item and bucket and this is reflected in constraints (17). In the formulation, (2n − 1) × |I| 2 + 2 × |I| binary integer variables, and n × ( i∈I |T i | + |I| + 1) fractional variables are defined. Moreover, total number of constraints except constraints (13) In the remainder of this section, the process of deriving the constrains (13) is introduced. Constraints (9), (10) and the variables in the constraints can be converted into a network as Figure 1; here, the constraints and variables are converted to nodes and arcs, respectively. Without constraints (13), a sequence-dependent setup usually causes subtours, indicated by the circle with arcs y i t , e   There are many studies on subtour elimination [18]. To obtain the required results, let us consider a converted network ( , ) for each bucket , as depicted in Figure 2. We generated the network using the following procedure: 1. contracting of the arcs relating to 's and merging the end nodes into one, denoted by ; 2. merging the tail nodes of arcs crossing the border between bucket − 1 and into one node, denoted by "source"; 3. and merging the head nodes of arcs crossing the border between bucket and + 1 into one node, denoted by "sink".
contracting of the arcs relating to y's and merging the end nodes into one, denoted by y i t ; 2.
merging the tail nodes of arcs crossing the border between bucket t − 1 and t into one node, denoted by "source"; 3.
and merging the head nodes of arcs crossing the border between bucket t and t + 1 into one node, denoted by "sink". Balas [22] introduced the prize-collecting traveling salesman problem to derive a model for scheduling the daily operation of a steel rolling mill. For the problem, Balas proposed the following generalized subtour elimination constraint: ≥ for all ∈ , ⊂ \{ , }, Therefore, we obtained the followings: Balas [22] introduced the prize-collecting traveling salesman problem to derive a model for scheduling the daily operation of a steel rolling mill. For the problem, Balas proposed the following generalized subtour elimination constraint: where e ij = 1 if setup changes from item i to j, 0 otherwise; δ + (S) is the set of outgoing arcs from S. Constraints (22) ensure that the number of selected outgoing arcs from S is at least equal to the number of selected arcs going out from any node k in S. Consequently, subtours can be prevented. Tacaari [18] rewrote constraints (22) as follows: where A(S) is the set of arcs between two nodes in S. In G t (V t , A t ), S ⊂ y i t , i ∈ I and A(S) ⊂ e ij tt , i j ∈ I (t, t) ∈ T ij ; therefore, by replacing the (i,j)∈δ + (i) e ij with y i t (this relation comes from constraints (10)), we can rewrite constraint (23) to fit into our formulation as follows: Through generating constraints (24) for each bucket t, we can derive constraints (13), the subtour elimination constraint.

Branch-and-Cut Algorithm
There are exponential numbers of subsets S ⊂ I; therefore, there are exponential numbers of constraints (13) in the formulation. To solve a formulation with exponential number of constraints, we usually use the branch-and-cut algorithm. In the algorithm, most of the constraints are outside of the running formulation and constraints violated by the current solution of linear programming (LP) relaxation are identified dynamically and added into the running formulation. Therefore, we can control the size of the running formulation. The procedure of identifying the violated constraints is referred to as separation (details regarding separation and branch-and-cut algorithm are available in [21]).

Separation
After obtaining an optimal solution to the LP relaxation (let e ij st , y i t denote the values of the variables e ij st , y i t , respectively, in the optimal solution) of the current running formulation, for each bucket t = 1, · · · , n, we use two separation routines to find the violated subtour elimination constraints. One routine is based on the minimum cut, and the other is based on the strongly connected components. In the routine based on the minimum cut, we generate a network G t V t , A t , illustrated in Figure 3, Here, A 1 t , A 2 t are defined as follows: For each arc e ij tt ∈ A 1 t , we assign capacity e ij tt and for each arc y i t , sink ∈ A 2 t , we assign capacity Processes 2020, 8, x 9 of 19 After obtaining an optimal solution to the LP relaxation (let ̅ , denote the values of the variables , , respectively, in the optimal solution) of the current running formulation, for each bucket = 1, ⋯ , , we use two separation routines to find the violated subtour elimination constraints. One routine is based on the minimum cut, and the other is based on the strongly connected components. In the routine based on the minimum cut, we generate a network ̅ ( , ̅ ), illustrated in Figure 3, where ≡ , ∈ > 0 ⋃{ } and ̅ ≡ ⋃ . Here, , are defined as follows: ≡ , ∈ , > 0 .
For each arc ∈ , we assign capacity ̅ and for each arc , ∈ , we assign capacity − ∑ ̅ ∀ | ∈ . For each vertex ∈ \{ }, the minimum cut between and yields a node partition ( , \ ), where ∈ and ∈ \ . If the total capacity of the arcs from to \ is less than , the constraint (24) with = is violated by the current solution. Thus, we add the constraint into the running formulation. For each time bucket = 1, ⋯ , , we determine the minimum cut | | − 1 times. The time complexity of finding a minimum cut (or maxflow) is | | | ̅ | by using Goldberg and Tarjan algorithm [23]. Therefore, the overall time complexity is | | | ̅ | .
This routine needs a considerably long computational time but always succeeds in finding the violated constraint.
In the separation routine based on the strongly connected components, we generate a network ( , ), where ≡ , ∈ > 0 and ≡ , ≠ ∈ ( , ) ∈ , ̅ > 0 . In the network, we could find some strongly connected components with more than two vertices. In each component , we select the vertex ∈ with the largest value ; then, for and = , we check if constraint (22) is violated by the current solution. If the constraint is variolated, we add it into current running formulation. The strongly connected components are found in a + operations by Tarjan's algorithm [24], and + operations are needed to check the violation. Therefore, the overall time complexity is . This routine is significantly fast. However, it successfully determines the violated constraint only in the case of integer solution. More information about the separation routine based on the strongly connected components is available in [25]. For each vertex y i t ∈ V t \{sink}, the minimum cut between y i t and sink yields a node partition S, V t \S , where y i t ∈ S and sink ∈ V t \S. If the total capacity of the arcs from S to V t \S is less than y i t , the constraint (24) with k = y i t is violated by the current solution. Thus, we add the constraint into the running formulation. For each time bucket t = 1, · · · , n, we determine the minimum cut V t − 1 times.

Implementation
The time complexity of finding a minimum cut (or maxflow) is O V t 2 A t by using Goldberg and Tarjan algorithm [23]. Therefore, the overall time complexity is O V t 3 A t . This routine needs a considerably long computational time but always succeeds in finding the violated constraint.
In the separation routine based on the strongly connected components, we generate a network In the network, we could find some strongly connected components with more than two vertices. In each component S, we select the vertex y i t ∈ S with the largest value y i t ; then, for S and k = y i t , we check if constraint (22) is violated by the current solution. If the constraint is variolated, we add it into current running formulation. The strongly connected components are found in a O V t + Â t operations by Tarjan's algorithm [24], and O V t + Â t operations are needed to check the violation. Therefore, the overall time complexity is O Â t . This routine is significantly fast. However, it successfully determines the violated constraint only in the case of integer solution. More information about the separation routine based on the strongly connected components is available in [25].

Implementation
The formulation and the branch-and-cut procedure are implemented using the Gurobi LP/MIP solver 9.0.2 with Python 3.7. The separation routines are executed by the embedded callback function. To find the minimum cut (or maxflow) and the strongly connected components, the efficient implementations in the Python library, Python-igraph 0.8.2, were used.
The following separation scenario was used in the computational experiments: • For integral solutions, the separation routine based on the strongly connected components generates the violated inequalities for each time bucket; • For fractional solution founded in i-th branch-and-bound nodes; I If i ≤ 1000, the separation routine based on the strongly connected components generates the violated inequalities for each time bucket. In the absence of a violated inequality, the separation routine based on the minimum cut is performed for each time bucket.

I
If 1000 < i ≤ 10, 000, the separation routine based on the strongly connected components generates the violated inequalities for each time bucket.

Computational Experiments
This section describes three sets of computational experiments. The first is for comparison with MLOV-SM in [19], the CLSP with sequence dependent setup, setup carryover and "setup-splitting". The second is for the analysis of performance of M1 with test instances having long setup times. The third is for comparison with MM2 in [1], the CLSP with sequence dependent setup, setup carryover and "setup crossover".
Experiments were conducted on a Windows 2010 PC with 16 GB main memory, which was equipped with an Intel(R) Core(TM) i5-9500 @ 3.00 GHz CPU. Default parameter settings for Gurobi LP/MIP optimizer are used, except 600 s of time limit.

Test Instances
For the first and second tests, three groups of artificial test instances ("T", "L", "LA") were generated. The instances in "T" have short setup times and were designed to be used in the first test. The instances in "L", "LA" have long setup times and were designed to be used in the second test. They have the following features:
Large backorder cost was given to bucket n − 1 to minimize production in the last bucket.
I c ij = τ ij 10 , for i, j = 1, · · · , m. • Demand: for problem dimension (m, n) For each group and problem dimension, initial demand set, named as (1), was generated manually according to the following criteria: All production can be finished before the last bucket.
To increase the reliability of the experiments, for each initial demand set, the due dates were randomly perturbed to generate four additional demand sets, named as (2), (3), (4), (5). The demands of an item having the same due date are merged into one.
For the last test, a test instance introduced in [1] was used. All the above test instances are available on [26].

First Test: Comparison with the Setup-Splitting Model (MLOV-SM)
As mentioned earlier, MLOV-SM is the model for the CLSP with sequence dependent setup, setup carryover and "setup-splitting". In addition, it considers situations in which multiple lots can be generated for an item within a time bucket. For direct comparison with our model (M1), MLOV-SM was modified and represented in Appendix A. The modified features are as follows: • Setup operation from i to i is not allowed; • At most one lot can be generated for an item within a time bucket; • The initial setup state is not given but determined by solving the model. M1 has (2n − 1) × m 2 + 2m binary integer variables and n × i∈I |T i | + m + 1 fractional variables, whereas MLOV-SM has nm 3 + 2nm 2 + 2nm binary variables and 3nm + 2n real variables, much more than M1. Because the setup time of the data instance in the group "T" is less than the capacities of the buckets, we can directly compare the lower bounds of the two formulation. We compare the two formulations with respect to the size, tightness of the LP relaxation bound and solution quality. Table 2 shows the tightness of the LP relaxation bound of the formulations. The column headers of the table are given as follows: • n: number of buckets; • m: number of items; • LB0: the objective value of LP relaxation without any additional cuts except GSECs; • LB: the objective value of LP relaxation with cuts generated by Gurobi; • BEST OBJ: the best integer solution value found within the time limit. Bold means optimal; • Gap0: relative gap between LB0 and the best integer solution value (BEST OBJ) at the end.
GAP0 improved: the difference in Gap0 between MLOV-SM and M1. In  Tables 3 and 4 show that M1 has much compact sizes of LP relaxation and gives much better solutions then MLOV-SM.  The second test is to see if M1 is stable in solving problems with long setup time. Tables 5 and 6 show test results for L and LA, respectively. The gap value '-' in Table 5 indicates that no feasible solution found in the time limit.  Tables 4-6 show that M1 reliably optimizes the problem instances with ten buckets and give relatively good feasible solution for the problem instances with twenty buckets.

Third Test: Comparison with the Setup Crossover Model (MM2)
Ramya, et al. [1] introduced the only research result on CLSP with sequence-dependent setup, setup carryover, and setup crossover. In the model, no idle times are allowed during setup operation and setup carryover. MM2 is a formulation introduced in [1] for the problem. But MM2 is too big to solve even problems of moderate size. MM2 has 8mn 3 + 4m 2 n 2 + 6mn 2 + 4m 2 n + 6mn + n binary integer variables and more than 150 types of constraints. Ramya, et al. [1] reported that it took more than 26,000 s to solve a test instance with 15 items and 10 buckets and that the solver failed with a file error after 8.5 h due to insufficient space in the hard disk while solving a test instance with 15 items and 15 buckets.
For the comparison with MM2, a modified version of M1 was developed and used for the test. The modified formulation, referred to as M2, needs two types of new variables described in Table 7.

Conclusions
In this paper, the CLSP with sequence dependent setup, setup carryover, and setup crossover was considered. A new mixed integer programming formulation was introduced. The formulation is based on three building blocks: the facility location extended formulation [2]; the setup variables with indices for the starting and the completion time periods [14]; and exponential number of generalized subtour elimination constraints (GSECs) [18]. A modified version of the separation routine of [25] was adopted to generate the violated GSECs.
Three groups of artificial test instances and an instance from [1] were used for computational experiments. The computational results showed that the proposed formulation outperformed the models from the literature, but the formulation has still large number of binary setup variables, O n × |I| 2 , to solve the problems with more than twenty buckets in 10 min.
To overcome this drawback, studies to reduce the number of setup variables are needed. These studies could involve using the column generation technique or heuristics based on variable fixing.
Funding: This study was supported by research fund from Chosun University, 2017.

Conflicts of Interest:
The author declares no conflict of interest.