Hybrid Bilevel-Lagrangean Decomposition Scheme for the Integration of Planning and Scheduling of a Network of Batch Plants

Motivated by a real-world industrial problem, this work deals with the integration of planning and scheduling in the operation of a network of batch plants. The network consists of single-stage, multiproduct batch plants located in diﬀerent sites, which can exchange intermediate products in order to blend them to obtain ﬁnished products. The time horizon is given and divided into multiple time periods, at the end of which the customer demands have to be exactly satisﬁed. The planning model is a simpliﬁed and aggregate formulation derived from the detailed precedence-based scheduling formulation. Traveling Salesman Problem (TSP) constraints are incorporated at the planning level in order to predict the sequence-dependent changeovers between groups of products, within and across time periods, without requiring the detailed timing of operations, which is performed at the scheduling level. In an eﬀort to avoid solving the full-space, rigorous scheduling model, especially for large problem sizes, two decomposition strategies are investigated: Bilevel and Temporal Lagrangean. We demonstrate that Bilevel Decomposition is eﬃcient for small to medium problem instances and that further decomposition of the planning problem, yielding a hybrid decomposition scheme, is advantageous for tackling a large-scale industrial test case.


Introduction
The Process Systems Engineering (PSE) and Operations Research (OR) communities have contributed to the development of models and solution strategies to tackle the decisionmaking on three different time scales in Enterprise-Wide Optimization (EWO) problems. The long-term decisions (years) belong to the strategic level and are usually related to investments, plant capacity expansion and retrofit; the medium-term decisions (months, years) pertain to planning and commonly refer to production planning and inventory control; the short-term decisions (days, weeks) represent the scheduling of operations, i.e. determining the detailed timing and sequencing of operations. Therefore, the integration of different time scales within a mathematical optimization framework becomes valuable for the decision-making process in an enterprise. In order to make this integration efficient for solving real-world problems, special-purpose models and decomposition algorithms need to be investigated. Reviews by [1] and [2] provide more details on the relevance of EWO activities in academic research and industrial practice. The focus of this work is on integrated planning and scheduling for batch processes. Excellent reviews on each of these levels are available in the literature, for instance [3], [4], and [5]. Planning models are usually posed as Linear Programming (LP) or Mixed-Integer Linear Programming (MILP) problems. For instance, [6] present several MIP formulations of the lot-sizing production planning problem. The main characteristics of such models are captured in constraints related to production targets, production and inventory costs, and material/inventory balances. Usually, those considerations are satisfactory at a planning level, which may span a time horizon of months to years.
However, when applied to multi-product plants such planning models lack a rigorous treatment of the production sequencing, which if ignored may underestimate the total costs and yield a plan that is not feasible. Without adding considerable complexity to the planning model, [7] proposed a Detailed Planning (DP) model that utilizes Traveling Salesman Problem (TSP) constraints to generate a cyclic schedule, which is broken in one link to yield an optimal sequence with minimum changeover times. The formulation allows sequencedependent changeovers across time periods. The planning model used in this work is based on the DP model. [8] have also investigated the estimation of the sequencing using TSP constraints.
At the scheduling level, a number of models that use either discrete-or continuous-time representation of events have been proposed. A review of both representations including their mathematical formulations can be found in [9]. We will focus on continuous-time representation models in this work. For sequential batch/continuous processes, two event representations have received a great deal of attention: time slots and precedence-based.
The main idea in the use of time slots is that each product can be assigned to a specific slot that has variable length. In some cases, especially in continuous processes, the number of time slots to be allocated for each reactor is known a priori. However, in batch processes that may not be true and additional time slots have to be allocated, thus increasing the size of the problem [10]. Likewise in the aforementioned DP formulation, the works by [11] and [12] has the desirable feature of sequence-dependent changeovers across adjacent time periods, which even though adds more complexity to the model due to the larger number of binary and continuous variables, renders it more realistic by not imposing a "hard" barrier for changeovers across time periods.
Unlike time slot models, precedence-based scheduling models effectively model changeovers by means of Big-M constraints. The two types of precedence relationships are local (immediate) and global (general) that are used to track the relative position of products in the production sequence. Four types of precedence-based models have been proposed: Unit-Specific Immediate Precedence (USIP) [13], Immediate Precedence (IP) [14], General Precedence (GP) [15], and Unit-Specific General Precedence (USGP) [16]. Briefly, IP and USIP differ in that the latter only takes into account the immediate precedence of products that are assigned to the same processing unit, whereas the former does not. The GP model generalizes the concept of precedence by accounting for all precedence relations (immediate or not) and requires fewer binary variables than the immediate precedence models, but cannot account for changeovers costs. To overcome that limitation, the USGP model was proposed. This work presents two extensions to the USGP model, namely: variable number of batches and sequence-dependent changeovers across time periods.
When attempting to solve large-scale industrial problems, some authors have identified the need to decompose the problem into subproblems. For instance, [17] proposed a Bilevel Decomposition algorithm, which involves solving an aggregate formulation in the upper level and a detailed formulation in the lower level with some decisions fixed by the upper level. Lagrangean Decomposition, a particular case of Lagrangean Relaxation [18], has also been investigated by [19] who decomposed their planning problem temporally and spatially in order to solve large instances. A review of methods and decomposition approaches for the integration of planning and scheduling can be found in the work by [20].
In this work, we deal with the integration of planning and scheduling of a network of batch plants motivated by a real-world industrial problem. The network consists of singlestage, multiproduct batch plants located in different sites, which can exchange intermediate products in order to blend them to obtain final products. Sequence-dependent changeover data are given in terms of groups or families of products, which are applied at the planning level. At the scheduling level a detailed model to obtain the timing and sequencing of batches of individual products is applied. It should be noted that at the scheduling level there can also be sequence-dependent changeovers between products belonging to the same group that are accounted for but these are of smaller magnitude than changeovers among different products. Hence, changeovers among products belonging to the same group are not considered at the planning level. The time horizon is given and divided into multiple time periods, at the end of which the customer demands have to be satisfied. We show that Bilevel Decomposition, Lagrangean Decomposition and a hybrid of those methods can be efficient when solving large-scale problems.
The remainder of the paper is organized as follows. In section 2, we define the manufacturing problem. In section 3, we present the complete planning and scheduling formulations and outline the major assumptions in our models. In Section 4, we discuss how the Bilevel Decomposition, Lagrangean Decomposition and a hybrid of those schemes are integrated. Section 5 contains computational results of three problem instances of increasing size. Section 6 summarizes the main conclusions of the paper.

Problem Definition
The problem addressed in this work concerns the planning and scheduling of a network of multi-product, single-stage batch plants as shown in Figure 1. The problem statement is as follows: a set of raw materials is purchased and transported to each plant, which transforms them into products to fulfill specified customer demands. Shipments of intermediate products between plants are allowed, hence the network structure depicted in Figure 1. Products are classified into groups or families, and we identify a product i belonging to group g through the subset I g . Sequence-dependent changeover data are given between groups of products, which are considered at the planning level. For the scheduling level, changeover times and costs associated with the transition from one product to another of the same group are also given. As indicated above, the sequence-dependent changeovers between different groups of products are significantly higher than the ones from one product to another in the same group. Additionally, the assignment of a product to a plant that does not normally produce it incurs a scale-up cost, which for example can be associated with the testing of new batches of the product in order to meet the quality standards set by the customer. The goal is to minimize total costs over a time horizon given by time periods of months in which the demands are specified at the end of each period.   Given batch sizes and fixed processing times, sequence-dependent changeover times and costs, transportation and inventory costs, and demand forecasts over a time horizon consisting of several time periods, the objectives are to determine: (c) The detailed timing of operations and sequencing of products in each unit, taking into consideration sequence-dependent changeover times and costs.
The goal is to minimize total cost -including operating, transportation, inventory, changeover costs -in order to meet customer demands at the end of each time period.
In the next section, we present the planning and scheduling models. The planning formulation not only contains material and inventory balances for products (production planning), but also estimates the sequencing of different groups of products through the use of TSP constraints. The scheduling formulation contains equivalent material and inventory balances for products and obtains the detailed production schedule of products, including the timing of operations and accounting for all possible changeovers. The two formulations are integrated through decomposition schemes (Bilevel and Lagrangean) in Section 4.

Problem Formulation
The plants' topology is sequential, i.e. single-staged plant with parallel units. The following assumptions are made: i. Inventories are stored in dedicated vessels with finite capacity; ii. Times for material transfer throughout the plant are considered to be negligible compared to processing and changeover times; iii. Batch sizes are fixed and given; iv. Batch processing times are given and unit-dependent; v. Demands that must be exactly satisfied are specified at the end of each period; vi. Sequence-dependent changeovers are group-and unit-dependent; vii. All data are deterministic.
The kernel of the planning formulation is based on the Detailed Planning model proposed by [7], which involves the classical Traveling Salesman Problem (TSP) sequencing constraints. The scheduling model is based on a modification of the Unit-Specific General Precedence (USGP) model proposed by [16]. Also, sequence-dependent changeovers between products families or groups have been included at the planning level [12].

Assignment and Sequencing within Time Periods
The assignments of group g to unit u in plant l in time period t are given by binary variables yg gult defined as: Likewise for products, the assignments of product i to unit u in plant l in time period t are denoted by binary variables yp iult defined as: Therefore, constraints (1) state that if product i is assigned to unit u at time period t, then the group to which it belongs is also assigned to the same unit and time period for plant l. In addition, according to constraints (2), yg gult is forced to zero if no product i of group g is assigned to any unit u at time period t in plant l.
Constraints (3) activate the binary variables responsible for identifying the products that have been assigned to unit u of plant l that does not normally produce it. When this happens, it will incur scale-up costs that for example can be related to the testing phase of new batches of the product in order ensure its quality.
The main idea for the sequencing of products groups is to generate a cyclic schedule for each time period that minimizes changeover times among the assigned groups, and then to determine the optimal sequence of groups by breaking one of the links in the cycle as illustrated in Figure 3. The 0-1 variables zg gg ult represent the changeovers between groups g and g in unit u in plant l in time period t. The cyclic schedules are generated with constraints (4) and (5), the assignment constraints of the TSP constraints, which state that for each plant l, group g is assigned to unit u during period t if and only if there is exactly one transition from group g to product g in unit u in time period t, and group g is assigned to unit u in period t if and only if there is exactly one transition from any group g to group g in unit u in time period t, respectively.
Only one link within a cyclic schedule is allowed to be broken (zzg gg ult = 1) and that is represented by constraints in (6). However, a link can only be broken if the corresponding pair is selected in the cycle as according to constraints (7).
In equations (4) and (5) the groups g and g can refer to the same group. In order to properly allow self changeovers constraints (8) -(10) have to be considered. These constraints state that for each plant l, if group g is assigned to unit u in time period t, and none of the groups g different from g are assigned to the same unit in the same time period, then batches belonging to group g can be followed by another set of batches also belonging to the same group. Also, if batches of group g are followed by other batches of the same group in unit u in period t, then only group g is assigned to unit u for period t and none of the groups g other than g are assigned to the same unit in the same time period.

Sequencing of Groups across Time Periods
In order to model sequence-dependent changeovers across time periods, it is necessary to identify the first and last products groups in each sequence. Therefore, two binary variables are introduced: ygf gult and ygl gult (see Figure 4).
Head Tail Figure 4: Breaking the cyclic schedule to determine the first and last groups of the sequence.
If at least one of the links between groups g and g is broken, then group g becomes the first group in the optimal sequence for unit u during time period t as stated in constraints (11). Similarly, constraints (12) imply that if at least one of the links between group g to any group g is broken, then group g becomes the last group in the optimal sequence for unit u during time period t.
Moreover, exactly one group can be the first and the last group to be processed in each unit as represented in equations (13) and (14).
Sequence-dependent changeovers across time periods are modeled as follows. For each plant l, equations (15) state that exactly one changeover occurs from group g to group g at the end of time period t in unit u if and only if group g is produced last in time period t. Similarly, according to equations (16), exactly one changeover from group g to group g occurs at the beginning of time period t + 1 in unit u if and only if group g is produced first in unit u in time period t + 1.
wheret is the last time period.

Time Balances
The time balance constraints (17) enforce that the total allocation of production times plus the total changeover times do not exceed the length of each time period. The usage time of each unit u in plant l in time period t is then calculated by adding the batch time to the changeover times, which in turn comprise the changeovers within and across time periods minus the changeover time of the broken link in the cyclic schedule.

Material and Inventory Balances
Upper and lower bounds on the production amounts for each product are enforced in constraints (18) and (19).   (20) - (24). If an intermediate product i is not to be blended with other intermediate products to form blend i , then its blending ratio is set to zero, i.e. BR ii = 0, which in the implementation is represented with a subset with which the corresponding constraint is eliminated from the model.
The material balances on the three nodes in Figure 5 for finished and intermediate products are expressed in constraints (25) - (27).
The inventory balances for all products are performed and their capacities are enforced in constraints (28) -(31).

Capacities
Each plant has capacity limitations for individual products for the entire time horizon as stated in constraints (32).
Upper bounds on the number of scale-up decisions per each plant and for all plants are imposed by constraints (33) and (34).
The number of assignments of products to plants is limited based on operational constraints through constraints (35).

Demand Satisfaction
The demands of each customer l are met according to constraints (36) -(38), which enforce that exactly one plant can sell a certain product to a customer.

Objective Function
The objective is to minimize the total cost, which includes the following terms: i.

Modified USGP Scheduling Model
We present two main enhancements to the Unit-Specific General Precedence (USGP) model proposed by [16]: the treatment of the number of batches as variable and the introduction of binary variables that model sequence-dependent changeovers across time periods. As opposed to the planning model, the scheduling model provides more detailed information concerning the timing and sequencing of individual products batches by accounting for all possible changeovers.

Assignment Constraints
Constraints (40) imply that every unit in a given plant and time period must process at least one product.
i∈I ul

Sequence-Dependent Changeovers Constraints
The planning model is an aggregate formulation based on the rigorous scheduling model and is only concerned with sequence-dependent changeovers between groups of products within and across adjacent periods. At the scheduling level, we compute the sequence of products rigorously. Therefore, the product-by-product changeover times, CTp ii ul , and costs, CCp ii ul , are computed based on the corresponding group-by-group parameters, CT gg ul and CC gg ul . The original USGP model defines two sets of 0-1 variables for sequence-dependent changeovers within time periods: zp ii ult (global or general precedence) and xp ii ult (local or immediate precedence). In this work, we introduce the 0-1 variables for sequence-dependent changeovers across time periods, wp ii ult .

Sequence-Dependent Changeovers within Time Periods
The big-M constraints (41) ensure that the start time of a product is at least the end time of another product plus the changeover time between them.

Sequence-Dependent Changeovers across Time Periods
In order to model changeovers across time periods, we introduce two sets of 0-1 variables that identify the first, ypf iult , and last, ypl iult , products assigned to a given unit in a plant. Constraints (42) and (43) guarantee that a product can be the first or last of a unit if it is assigned to the same unit. Moreover, there can be exactly one first product and exactly one last product assigned to a given unit of a plant as represented in equations (44) and (45), respectively.
With respect to timing of operations, the main idea is to obtain the start and end times of each unit in each plant, T sU ult and T eU ult respectively, and ensure that the changeover time from the last product in a time period and the first product in the subsequent time period is taken into consideration (see Figure 6). Constraints (46) ensure that the start time of product i in time period t + 1 is at least the end time in the previous time period t plus a changeover time in unit u. Similarly, constraints (47) guarantee that the start time in unit u in time period t + 1 is at least the end time of product i in time period t plus a changeover time.

Sequencing-Allocation Constraints
The global sequencing variables, zp ii ult , are activated by the following constraints. The logic proposition (48) states that product i precedes product i in unit u of plant l and time period t or product i precedes product i in the same unit, plant, and period if and only if both products are assigned to the same unit.
The constraints correspondent to the above proposition are as follows: The local sequencing variables, xp ii ult , are activated through the Big-M constraints (54) and (55), which state that two products i and i are consecutive only in the case that the binary variable zp ii ult = 1 and when there is no other product i between products i and i , and vice versa. The idea is to count how many products i , where i = i = i , are in between products i and i and, thus, track the relative position of products that are assigned to a unit.
where pos ii ult is an nonnegative continuous auxiliary "position" variable and M is a Big-M value, which could be the total number of products assigned to the respective unit, plant, and time period. Therefore, if pos ii ult = 0, that is there are no products between products i and i , then we force xp ii ult = 1 to ensure that product i immediately precedes product i .
Lastly, the assignment variables relative to the first and last product for each time period are related to the sequencing variables for the changeovers across time periods as follows:

Time Balance Constraints
Equations (58) state that the end time of product i in unit u of plant l in time period t is the sum of its start time and processing time.
where P T iult can be modeled as the following set of disjunctions: Given that there is a one-to-one correspondence between Y P iult (logical variables) and yp iult (binary variables), the convex hull [21] reformulation yields the following linear constraints: where N B1 iult and N B2 iult are disaggregated variables and NBULP iult represents the number of batches of products obtained in the Upper Level Planning problem. In addition, constraints (64) -(67) force the start and end times of a given unit in a time period to coincide with the start and end times of the first and last products, respectively.
Finally, constraints (68) and (69) ensure that the start and end times lie within the current time period length.

Changeover Costs within Time Periods
Given the changeover costs for each unit, the costs associated with the changeovers within time periods are given by:

Changeover Costs across Adjacent Time Periods
Similarly, the changeover costs across adjacent time periods are calculated as follows:

Material Balances
The constraints for material and inventory balances are equivalent for the planning model, equations (20) and (25)

Objective Function
The objective function considers the same terms as in the ULP model except for the changeover costs and it can be written as follows (see section 3.

Bilevel Decomposition
In order to solve both the planning and scheduling problems simultaneously for medium-tolarge problem instances in practical time, a decomposition scheme can be employed. One approach for decomposing this problem is the Bilevel Decomposition (BD) consisting of an Upper Level Planning (ULP) and a Lower Level Scheduling (LLS) problems that yield lower and upper bounds on the total cost, respectively [10]. The problems are solved iteratively until the relative difference between the lower and upper bounds is less than a pre-specified tolerance. Integer inequalities or cuts are added to the upper level problem to ensure the generation of new schedules, and/or to avoid infeasible ones at the lower level problem. A schematic of this strategy is shown in Figure 7.

Upper Level Planning (ULP) Determine lower bound (LB) on cost.
Solve for entire network. The main motivation behind the BD is that the ULP problem is less complex and, thus, easier to solve than the LLS problem. Therefore, by solving the aggregate and simplified planning problem and iteratively attempting to match its solution with the one obtained with a rigorous and detailed scheduling model, we strive to arrive at a feasible production plan without the burden of solving the full-space scheduling model.
In addition to computing a lower bound on the total cost, the ULP problem determines which products may be processed in a given unit and also the number of batches of each product on every unit and time period, which will be inputs to the LLS problem. Hence, the scheduling problem at the lower level is solved only for the product assignments predicted by the upper level, which may significantly reduce the number of constraints and variables in the LLS problem.
It may be advantageous to add different integer cuts to the ULP problem after each iteration, such as superset and subset cuts [22] and symmetry-breaking cuts [10]. However, from our experience gained from solving the present problem, at most two iterations were required to converge the ULP and LLS problems to less then 1% relative error even for large problems. Hence, the integer cuts to exclude the assignments at the previous BD iteration take the form: where k is the index of BD iterations that are contained in set BD k (maximum of 10 iterations were used in all computational experiments), B k iult = {(i, u, l, t) ∈ I iul ×U l ×LP l ×T : yp iult = 1 at iteration k−1} and N k iult = {(i, u, l, t) ∈ I iul ×U l ×LP l ×T : yp iult = 0 at iteration k−1}.

Lagrangean Decomposition
Problems with multiperiod formulations or with variables associated with spatially distributed entities, such as plants and customers, are candidates that can undergo Lagrangean Decomposition (LD), which is a special case of Lagrangean Relaxation (LR) where the original problem can be decomposed into subproblems with common variables by splitting them first, and then dualizing the copy constraints. In this work, since the ULP can become the bottleneck as it has to be solved for the entire network, unlike the LLS that can be solved for each plant, we decompose the ULP problem using LD inside the BD loop discussed in the previous section. This proved to be computationally advantageous when solving the large-scale industrial problem. Two types of LD have been explored in the planning and scheduling literature: Temporal Lagrangean Decomposition (TLD) and Spatial Lagrangean Decomposition (SLD). As the names imply, the former is characterized by decomposing the original problem into subproblems corresponding to each time period, whereas the latter is decomposed by the spatially distributed sites. Due to the presence of binary variables in MILP problems, i.e. a source of nonconvexity, there is a duality gap between the solution of the Lagrangean dual problem and the primal problem. Moreover, it has been shown by [19] that the temporal dual bound is at least as tight as the spatial dual bound. Therefore, in this work we focus on the TLD only. Figure 8 shows a schematic of TLD. Each time period corresponds to an independent subproblem that can be solved in parallel with the other subproblems. The variables z t in this figure represent inventory variables, IN V ilt and BIN V ilt , and the assignment variables ygf gult in the ULP formulation because they appear in terms in time period t and t + 1, that is they link consecutive time periods. Since we need to separate time periods into unique subproblems, we introduce the copy variables IN V ilt = IN V ilt , BIN V ilt = BIN V ilt , and ygf gult = ygf gult . The next step is to dualize these equations by multiplying each of them by the respective Lagrange multipliers λIN V ilt , λBIN V ilt and λygf gult , respectively, and adding the terms to the objective function.
For fixed multipliers, the planning model is separable into time periods and can be written as follows: and subject to constraints (1) - (15), (17) - (27), (29), (31) -(35), and the following constraints, which are rewritten by including the copy variables: At the end of a TLD iteration, the multipliers must be updated for the next iteration. Different strategies have been proposed, such as cutting plane [23], subgradient [24], boxstep [25], bundle [26], analytic center cutting plane methods [27], and volume algorithm [28]. In this work, we use the subgradient method in which the update formula for a multiplier λX p at iteration p ∈ T LD p of the TLD and associated with the equality constraint of generic variable X is given by: where p is the step size that can be modified at each iteration according to some criterion and usually lies in the interval (0, 2], U B TLD is the best upper bound in the TLD scheme up to iteration p obtained from solving the original planning problem in the reduced space (i.e. by fixing variables from the solution of the Lagrangean subproblems), LB TLD is the best lower bound in the TLD scheme up to iteration p obtained from summing all the optimal objective function values of the Lagrangean subproblems, and den p is the sum of squared of the differences between the original variable X p and its copy X p and for the planning formulation is defined as follows: Lastly, the multipliers may be initialized prior to the first TLD iteration to zero or to other values, for example the marginal values of the respective equations from the relaxed TLD model (integrality conditions are relaxed). From our experience, the best choice of initial values may vary among different problems. The next section contains the full description of the algorithm implemented when TLD is applied to the ULP problem within the BD loop.

Hybrid BD-TLD Scheme
The two decomposition methods, BD and TLD, are combined in an effort to solve large-scale industrial problems. Figure 9 shows the schematic of the hybrid approach.
The main steps of the hybrid decomposition scheme is given in Figure 10. We underscore the advantage of creating smaller and independent subproblems through TLD, because they can be solved in shorter time as compared to the full problem and also in parallel. Particularly in this work, since each subproblem corresponds to a time period, it means that it is possible to solve all subproblems at the same time in a multi-core computer instead of solving each subproblem one at a time. To our knowledge, parallel "solve statements" have become supported by two of the major modeling platforms, GAMS [29] and AIMMS [30]. We used GAMS to implement our models and algorithms. For more information on how to solve problems in parallel in GAMS, see the website: http://interfaces.gams-software.com/ doku.php?id=the_gams_grid_computing_facility.

Computational Results
Three example problems of increasing size and complexity were solved to demonstrate the efficiency of the decomposition approaches discussed in the previous section. Example 1 is a small-scale problem for which we provide all data as well as present the optimal Gantt chart obtained. Example 2 is a medium-scale problem and Example 3 is a large-scale industrial problem. All examples were solved using BD and the full-space scheduling model. We used the hybrid BD-TLD algorithm only to solve Example 3.
All models were implemented in GAMS 23.8.2 and solved with Gurobi 4.6.1. All computational experiments were performed in a Dell PowerEdge T410 server with 6 Intel R Xeon R 2.67 GHz CPUs (total 12 threads), 16 GB of RAM, and running Ubuntu Server 12.04 LTS. The Gurobi's option threads 0 was enabled for all computational runs, which means that all threads were used for parallel processing. The maximum allowed wall time for all problems was 24 hours.

Example 1
Example 1 considers a small-scale problem, and it allows us to analyze not only the computational benefit of Bilevel Decomposition (BD), but also to gain insight on the optimal schedule represented by a Gantt chart. All the data are given in Appendix A and the characteristics of Example 1 are as follows: • All models were solved to optimality Figure 11 shows the optimal schedule represented by a Gantt chart for the first time period (week) only, where the letters inside the colored blocks represent the products and the numbers within parentheses denote their number of batches. The same result was obtained after solving both the LLS and the full-space scheduling problems. In both units U12 and U21, the last sequence-dependent changeovers occur across the first time period. We emphasize that the number of batches (integer variables) were optimized simultaneously with the planning and scheduling decisions in each model.  Table 1 shows the objective function breakdown for all models. Notice that the ULP problem predicted lower changeover costs than the scheduling models due to the aggregation of products into groups of products in the former model so that changeovers between products belonging to the same group are neglected. Moreover, the solution obtained for the LLS model was numerically the same as the one obtained for the full-space scheduling model denoted as "FS", which may be expected for small and "well-posed" problems. The BD gap between the LLS and the ULP solutions was 0.071% and only one BD iteration was needed to achieve convergence. The problems sizes as well as computational statistics are shown in Table 2. The row "Nodes" indicates how many tree nodes in the branch-and-cut method were explored by the solver. Using BD, the total wall time was 1.17 seconds as opposed to solving the full-space (FS) scheduling model, which took 44.98 seconds. Notice that the LLS problem contains fewer variables and constraints than the FS problem, which facilitated its solution yielding only 29 nodes explored in the reduced space as opposed to nearly 95,000 nodes in the full space.

Example 2
Example 2 is a medium-scale problem. We compared the effectiveness of applying only the BD with solving the FS problem. The computational benefit of using BD becomes more evident as the problem grows in size and complexity as will be shown in the results below. The characteristics of Example 2 are as follows: • All models were solved to 0.5% optimality gap Table 3 shows the objective function breakdown for all models. Notice that the final solution of the BD, i.e. the solution to the LLS problem yields a lower total cost than the one obtained solving the FS problem. Also, no shipments between plants were observed in any of the problems. The BD gap between the LLS and the ULP solutions was 0.2% and only one BD iteration was needed to achieve such convergence. The problems sizes as well as computational statistics are shown in Table 4. Using BD, the total wall time was 3.96 seconds for a gap of 0.2% as opposed to solving the full-space (FS) scheduling model, which took approximately 3 hours and 25 minutes for 0.5% optimality gap. The difference in terms of number of variables and constraints between performing the BD and not decomposing the problem becomes considerably more expensive as the problem instance increases. That is directly translated into 0 nodes necessary to solve the LLS problem compared to nearly 60,000 nodes explored in the solution of the FS problem.

Example 3
Example 3 is a large-scale industrial problem. We identified the need to decompose the ULP problem and use the hybrid BD-TLD algorithm in order to get a solution in practical time. The following characteristics provide an idea of the problem size of Example 3. Exact characteristics are not provided due to confidentiality reasons: • More than 5 plants • All models were solved to 2% optimality gap Table 5 shows the objective function breakdown for all models except the FS problem, which could not be solved due to the excessive RAM required. A maximum of 30 TLD iterations were enforced. The BD gap between the LLS and the ULP solutions was 0.13% and only one BD iteration was needed to achieve such convergence. The problems sizes as well as computational statistics are shown in Table 6. Using hybrid BD-TLD, the total wall time was around 1 hour and 15 minutes. Using only the BD, the wall times were 6 hours and 12 minutes to solve the ULP problem and 1 hour and 26 minutes to solve the LLS problem (7 hours and 38 minutes total). It can be noted that the decomposition in the ULP problem allowed us to obtain a solution with significantly less computation time. We show the size of the FS problem even though we could not solve it. It is clear that decomposition approaches enable tackling real-world problems.  Figure 12 shows the evolution of the upper and lower bounds on the objective function value obtained when applying TLD to the ULP problem. The best lower and upper bounds after 30 iterations were $9, 363, 323.00 and $9, 440, 580.00, which represent a final TLD gap of 0.82%.

Conclusions
In this work, we developed a model for the integration of planning and scheduling in the operation of a network of multiproduct batch plants. Each plant contains single-stage processing units in parallel and the products are classified into groups or families. The planning model not only incorporates production and capacity constraints, but also estimates the sequencing of groups of products through Traveling Salesman Problem (TSP) constraints. The scheduling formulation employs continuous time representation in a precedence-based framework. It solves the batching and the scheduling problem simultaneously, i.e. the number of batches is a decision variable in the model, and allows sequence-dependent changeovers across time periods. Both models are multi-period Mixed-Integer Linear Programs (MILPs).
The integration of planning and scheduling was performed via two decomposition methods: Bilevel and Lagrangean. For small to medium problems, it was observed that Bilevel Decomposition (BD) represents an attractive and rigorous decomposition strategy to solve them in practical computational time. However, when attempting to solve a real-world problem, the planning level became computationally expensive and, thus, became amenable to further decomposition. We performed Temporal Lagrangean Decomposition (TLD) in the planning level in an inner loop of a BD framework to obtain a solution in a reasonable time.
The computational results of the examples have shown that detailed model formulations applied to large problems can be tackled through efficient modeling decomposition strategies and advances in computing, such as in optimization solvers capabilities, multi-core computer architecture and parallel computing.