Rolling Operation Algorithm for Solving Complex ... - ACS Publications

3898

Ind. Eng. Chem. Res. 2009, 48, 3898–3908

Rolling Operation Algorithm for Solving Complex Scheduling Problems Barbara Czuczai,† Tivadar Farkas,‡ Endre Rev,*,† and Zoltan Lelkes† Department of Chemical and EnVironmental Process Engineering, H-1111, Budapest UniVersity of Technology and Economics, Budapest, Budafoki ut 8, Hungary, and HAS-BUTE Materials Structure and Modeling Research Group, Department of Chemical Engineering, Budapest UniVersity of Technology and Economics, H-1521 Budapest, Hungary

Large-scale scheduling problems have remained difficult to solve despite the intense research in the area and many successful algorithms. In this paper, a new decomposition algorithm is presented which finds a good feasible solution even for large-scale problems in a reasonable time. In each iteration, a maximum allowed number of operations is first assigned to available units, then the algorithm alternates between a production quantity maximization step and a production time minimization step in each iteration. The performance of the “rolling operation” algorithm is illustrated using an industrial example. 1. Introduction Process scheduling problems usually involve a diverse data structure. (1) A set of different product orders are considered with specified amount to be produced before their deadlines. (2) A set of different raw materials with given amount (e.g., mass or volume) and release time data can be applied to satisfy the product orders. (3) A set of tasks, i.e., single activities consuming raw materials or intermediates and producing intermediates or end products constitute the technology used in the production, together with (4) processing time, etc., data in the production recipes. (One production recipe contains the characteristics of one task.) (5) A process recipe contains data about the process itself, i.e., relations and connections between tasks. (6) Plant specifications contain information about the characteristics of available processing units and storage tanks (from now on, they are referred to as the common term “units”), and the allowed connections between them. (7) Detailed technological constraints can be derived from the above data (capacity, storage restrictions, etc.). The goal is to find either a solution optimizing an objective function (makespan, etc.), or only a feasible solution that satisfies the constraints. Complicated problems allow splitting the product of single operations and mixing the products of separated operations; this implies that material balances have to be taken into account in the model. Such problems are referred to as general networkrepresented processes.1 The present paper deals with this general and complicated problem type. Process scheduling problems have been intensively researched; many articles are published in the chemical engineering and operational research literature. Several MILP models have been developed during the past 15 years for general network-represented processes. The earliest model, invented by Kondili et al.,2 applies discrete time representation, monolithic material balance, and global time intervals as event representation. Here the operations can start and end only at discrete time moments. In order to improve the performance of the original model, several other researchers elaborated different MILP formulations employing continuous time representation. * To whom correspondence should be addressed. E-mail: ufo@ mail.bme.hu. † Department of Chemical and Environmental Process Engineering, Budapest University of Technology and Economics. ‡ HAS-BUTE Materials Structure and Modeling Research Group, Budapest University of Technology and Economics.

One of the most recent works of this kind is published by Maravelias.3 He showed that the discrete time representation is a special case of the continuous one. He proposed a new mixedtime representation which combines the modeling advantages of both discrete and continuous time models. In another paper, Sundaramoorthy and Karimi4 presented a new slot-based model for handling multipurpose batch plants. Their new model represents time in terms of ordered blocks of unknown variable length. Although their model is considered in the literature as different from the others, this concept does not differ significantly from those being event based: each of these models monitors the level of shared resources at particular event points. Most of the researchers agree that the best results can be achieved with models using unit-specific events introduced first by Ierapetritou and Floudas.5 This is so due to the smaller number of binary variables and the better relaxation. The unitspecific event concept allows different tasks to start at different moments in different units for the same event point. Consequently, (1) the number of necessary event points is less than that in the case of the global event-based representations; and (2) the event point of finishing an operation is known from its start. The original version of this model has been updated and improved.6 A useful survey on different MILP models can be found in Mendez et al.1 Based on a different concept, Czuczai et al.7 recently developed a new model. It introduces the notion of operations of units. The material transport between operations is represented by binary variables, and these variables are used for modeling precedence relations between the operations. All the above models are able to handle small-scale problems only. In the past decade, growing interest has shown up to solve real size scheduling problems. In order to deal with industrialscale problems, an abundant number of MILP-based heuristics8,9 and decomposition approaches have been developed to enhance the efficiency of the models. Three types of decomposition exist according to Bassett et al.:10 those that can be decomposed are (1) time, (2) units, and (3) tasks/resources. The most widely studied decomposition method is the first one, which decomposes the problem into several subproblems along the time axis and solves them by an iterative way. In the first iteration a subproblem is solved over the first part of the horizon, while a remaining portion of the horizon is modeled

10.1021/ie801037k CCC: $40.75  2009 American Chemical Society Published on Web 03/12/2009

Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3899

Figure 1. Rolling horizon method.

in an aggregated way. In the next iteration the solution belonging to the preceding part of the horizon is fixed, and the next portion is solved in details. The actual detailed portion of the horizon is rolling in time, as is shown in Figure 1. Therefore, this kind of decomposition method is called “rolling horizon”. The subproblems, being smaller than the original one, are less difficult to solve and it is possible to obtain an optimal solution for them. However, the detailed solution built from them may not be proven optimal. Several time-based decomposition approaches are presented by Bassett et al.10 They developed an iterative scheme in which a planning problem is solved using an aggregate formulation as master problem which assigns production quantities to subperiods, and then a detailed scheduling problem is solved for each subperiod using the results of the planning problem as boundary condition. Their aim was to avoid the artificial end effect: as operations must be finished before the end of the actual subperiod, the production capacity of the plant is artificially decreased by the algorithm, making the scheduling problem infeasible in some cases although it might have been solvable originally. A similar time-based decomposition is presented by Wu and Ierapetritou.11 The unit-specific event formulation5,6 developed by Ierapetritou and Floudas is used in the resulting subproblems, and different heuristic-based decomposition methods are combined with a Lagrangean-relaxation approach in an iterative scheme. Janak et al.12 developed a two-level decomposition formulation for medium-term scheduling. In this method the first level decides about the number of days, and the products that are taken into account in the second (short-term) scheduling level. A unit-decomposition-based algorithm is published by Kelly et al.13 In their method the subproblems are assigned to unit groups consisting of units with similar characteristics. The subproblems belonging to unit groups are solved separately. The “required production method”, based on the partition of tasks, was published by Wu and Ierapetritou.11 Here the necessary amount of intermediates required to satisfy the product orders is targeted. A resource-based algorithm published in the

Figure 2. Process flowsheet of motivating example.

same article is based on the decomposition of the event points. First, a few event points are used and the scheduling problem is solved. Then the task sequence is fixed, based on the solution of the problem with those few event points, and additional event points are considered, and that problem is solved. The iteration terminates when no further improvement of the objective function is achieved. The decomposition strategy by Harjunkoski and Grossmann14 can also be categorized as a task-decomposition scheme. They first presort the products into product families and then disaggregate the product families into groups and schedule each of the groups independently. All these decomposition algorithms share the common disadvantage that their application may probably lead to a suboptimal and, more importantly, infeasible solution, especially when tight deadlines are present. In the case of time decomposition algorithms, an aggregate planning model is used to set up production goals in the subproblems; this generally overestimates the plant capacity, causing infeasibility in the subproblems. In this paper, a new decomposition algorithm is presented which can further reduce the above-mentioned drawbacks. The algorithm divides the original problem into subproblems; these subproblems are easy to solve due to their small size. The new decomposition is named “rolling operation” algorithm, as it is similar to the “rolling horizon” methods. A motivating example from the beer industry is first presented, highlighting its special characteristic inspiring the new algorithm. Then the model by Czuczai et al.7 is summarized so that the new concept applied in the present article can be introduced. The algorithm itself is then explained. Finally, the results of the industrial example are shown. 2. Motivating Example The example inspiring the algorithm, and used for studying its performance, derives from the beer industry, where there are processes attributed with special characteristics. The example process consists of three stages: (1) filtration, (2) BBT task, and (3) packaging. BBT itself means “bright beer tank”. The BBT tank will be distinguished from the BBT task, although the BBT abbreviation itself means “tank”. The flowsheet is shown in Figure 2. Raw beer is stored in raw beer tanks. The raw beer can be filtered in several different filters. A raw beer tank can be connected to only one filter at the same time, and a filter can be connected to only one raw beer tank at the same time. During filtration, the beer is immediately loaded to a BBT tank. A filter can be connected to only one BBT tank at the

3900 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009

Figure 3. (a) Fragment of a full superstructure: any connection is allowed. (b) Fragment of a simplified substructure. (c) Fragment of a solution substructure.

same time, and a BBT tank can be connected to only one filter at the same time. Filtration is a continuous task; the beer flows continuously through the filter; the raw beer tank and the BBT tank must be exclusively connected to the filter during the whole filtration process. The beer can be accumulated in the BBT tank from several filtration operations during a period for which an upper bound (maximum loading time) is specified. BBT itself is a batch task. After some minimum resting time in the BBT tank, the tank is emptied, and the beer is directly loaded to packlines. The load of the BBT tank can be packed in several different packing operations. The time period used for emptying the BBT tank is constrained from above by an upper bound (maximum emptying time). A BBT tank can be connected to only one packline at the same time, and a packline can be connected to only one BBT tank at the same time. Packaging is also a continuous task; the BBT tank and the packline must be exclusively connected to each other during the whole packaging process. Since filtration and packaging are continuous tasks, the starting and finishing time moments must be synchronized with

the corresponding loading and emptying operations of the actual BBT tank. From viewpoint of unit connectivity, this problem belongs to the “operational restricted” group according to the classification of Czuczai et al.7 because each unit can be connected to only one other unit of a particular type at the same time. The beer can be stored in the BBT tank, before and after the BBT “resting” task, only for specified maximum time periods. Therefore, this problem belongs to the “finite storage time” group. Deadlines (assigned by the warranty periods) for consuming the raw beer charges stored in each raw beer tank are given. Deadlines for satisfying the product orders are also known. There are several orders for each product. The industrial example contains 20 raw beer tanks (rbt01 to rbt20), 2 filters (f1 and f2), 9 BBT tanks (b1 to b9), and 3 packlines (BotL, CanL, and KegL). Three product types exist: bottles, cans, and kegs; each packline is able to pack one product type only; their names refer to the pack type. Twenty-four product orders for 15 products, and by different deadlines, are to be satisfied. The deadlines are specified such a way that for


satisfying the product orders, the machines should work almost permanently; and finding a feasible solution is difficult because of the tightness of the production plan. Because simultaneous connections between units are forbidden, connections between units (and connections between operations of units) are to be modeled. This inspired the development of the model by Czuczai et al.,7 where a number of operations (of different types) are assigned to units; binary variables are designated for modeling the type of task executed in an operation of a unit; and different types of binary variables are designated for modeling the existing connections between operations of different units. The model uses the natural concept of “operation” instead of event points used in other works; this leads to a different model formulation and a different decomposition algorithm. 3. Mathematical Model This section provides a short overview of the model given by Czuczai et al.,7 in order to make its concept understood, needed for explaining the algorithm presented in section 4. 3.1. Superstructure. 3.1.1. Operations. No event points are considered in the model. Instead, the notion of operation of a unit is introduced; see the rectangles in the superstructure fragment in Figure 3a. The operations of units are ordered along the time axis. We always speak about outh operation of type it of unit u. The type it of a unit can be either continuous or batch. We use a language that a task is processed in the outh operation of a unit. Binary variables ytitu,ou,i are used to denote the existence it , of task i in operation ou of unit u. Continuous variables Vtu,ou,i tsitu,ou, and tf itu,ou are assigned to each operation modeling its size, starting time moment, and finishing time moment, respectively. 3.1.2. Connections. As it is detailed below, the model has a new attribute which distinguishes it from previous literature models: it directly handles the connections between operations, denoted by arrows between the rectangles in the superstructure it1,it2 fragment in Figure 3a. Binary variables ycu1,ou1,i1,u2,ou2,i2 are assigned to connections between operations. A value of 1 represents existing connection between the two operations (i.e., if there is material flow between them, one uses the product of the other one as raw material); and zero represents no connection. The precedence constraints (i.e., which operation finishes earlier and which one starts later) are handled through these binary variables. The value of the binary variable of a particular connection between operations specifies the precedence relation between those operations, i.e., the equality and inequality constraints between their starting and finishing time moments. This connectivity variable is very similar to the one proposed by Prasad et al.15 There binary variables are used for modeling the connection between two units at time slots. 3.2. Formulation. For sake of simplicity and saving space, only the algebraic equations containing binary variables are presented here in a concise form. 3.2.1. Operations. 3.2.1.1. Only One Task in a Given Operation. If task i is being processed in outh operation of type it of unit u, then other tasks cannot be processed in the same operation.

Vtitu,ou,i

e CAPu,i·ytitu,ou,i u ∈ U; it ∈ TT; ou ∈ OUu,it ; i ∈ Iu,it

(2)

3.2.2. Connections. 3.2.2.1. Only One Connection during a Given Operation. In some cases, consuming the product of an operation is prohibited by more than one operation, independently from their type or whether those are performed in the same or different units. This is described by the following algebraic equation for inlet connections:

∑ ∑

∑

∑

ycit1,it2 u1,ou1,i1,u2,ou2,i2 e 1

u2∈U it2∈TT ou2∈OUu2,it2 i2∈Iu2,it2

u1 ∈ U; it1 ∈ TT; ou1 ∈ OUu1,it1 ; i1 ∈ Iu1,it1 ; i1 ∈ IOCIN (3)

3.2.2.2. Capacity Constraints between Operations. If there is no material sent from ou1th operation of type it1 of unit u1 to ou2th operation of type it2 of unit u2, then the material flow between the two operations is zero. If there is a material transport between them, then the material flow is constrained by the capacities of the two units. it1,it2 P C Vcu1,ou1,i1,u2,ou2,i2,m e min(Fm,i1 ·CAPu1,i1, Fm,i2 ·CAPu2,i2)· it1,it2 ycu1,ou1,i1,u2,ou2,i2 u1, u2 ∈ U; it1, it2 ∈ TT; ZOCu1,it1,ou1,u2,it2,ou2 ) 1;i1 ∈ Iu1,it1 ; i2 ∈ Iu2,it2 ; m ∈ MINi2 ∩ MOUTi1 (4)

3.2.3. Precedence Constraints. Processing tasks and storage tasks, as well as processing operations and storage operations, are distinguished in our terminology;7 constraints are related to connections between some combinations. Only some of them are listed here. 3.2.3.1. Precedence between Batch Processing Operations. If the material produced in batch operation ou1 of type it1 of unit u1 is used by another batch operation ou2 of type it2 of unit u2, then the starting time moment of operation ou2 of type it2 of unit u2 must not be earlier than the finishing time moment of operation ou1 of type it1 of unit u1. it tfu1,ou1 - tsitu2,ou2 e THOR·(1 -

∑ ∑

it,it ycu1,ou1,i1,u2,ou2,i2 )

i1∈Iu1,it i2∈Iu2,it

u1, u2 ∈ U; it ∈ TT; it ) BATCH; ou1 ∈ OUu1,it ; ou2 ∈ OUu2,it (5) 3.2.3.2. Synchronization between Continuous Processing Operations. If the material produced in continuous operation ou1 of type it1 of unit u1 is used by another continuous operation ou2 of type it2 of unit u2, the starting time moment of operation ou2 of type it2 of unit u2 must be equal to the starting time moment of operation ou1 of type it1 of unit u1. tsitu1,ou1 - tsitu2,ou2 e THOR·(1 -

∑ ∑

ycit,it u1,ou1,i1,u2,ou2,i2)


u1, u2 ∈ U; it ∈ TT; it ) CONT; ou1 ∈ OUu1,it ; ou2 ∈ OUu2,it (6) -THOR·(1 -

∑ ∑

it ycit,it u1,ou1,i1,u2,ou2,i2) e tsu1,ou1 -


(1)

tsitu2,ou2 u1, u2 ∈ U; it ∈ TT; it ) CONT; ou1 ∈ OUu1,it ; ou2 ∈ OUu2,it (7)

3.2.1.2. Capacity Constraints. If task i is not processed in outh operation of type it of unit u, then the size of that operation is zero. If task i is processed in outh operation of type it of unit u, then the maximum size of that operation is determined by the capacity of the unit.

For forcing the equality of finishing time moments, equations it it with tf u,ou . of the same form are given by substituting tsu,ou Similarly to the above presented equations, precedence relation may exist between batch, storage, loading, and emptying operations, in every matching. Those all require describing

∑ yt

it u,ou,i

e1

u ∈ U; it ∈ TT; ou ∈ OUu,it

i∈Iu,it


different combinations of relationships between starting and finishing time moments of different types of operations. Besides the above detailed constraints, the most important equations not presented here are (1) the inequality constraints that exclude overlapping of different operations of a unit; (2) the equations calculating the finishing time moments of operations from their starting time moments; (3) the material balances between operations; and (4) the constraints enforcing to meet the product demands and to satisfy the deadlines. 3.2.4. Simplifying the Superstructure. So far in the model, all the connections apart from those which are technologically forbidden are able to exist in the superstructure; see Figure 3a. This possibility is modeled by the binary parameter ZOCu1,it1,ou1,u2,it2,ou2. This results in a large number of possible connections which may lead to a large number of binary variables. This may make the solution process difficult. Therefore, it is worthwhile to tighten the set of the possible connections in a way that no potential optimum is excluded. The studied process contains (mainly continuous) operations allowed to be connected to one other operation only. In this situation, and in similar problems, there is a possibility to simplify the superstructure. If only those connections are allowed which connect operations of same sequence number, as it is shown in Figure 3b, then the number of binary variables is considerably decreased. This is made by setting the value of the corresponding ZOC parameters to zero. In order to avoid excluding a potential optimum, enough number of operations should be included in the superstructure. It can also be ensured that the optimum solution may be found for general network represented processes in multipurpose plants, but the necessary number of available operations is larger in that case. When there are uneven number of units per stage or even when there are no stages in the process, the solution structure will include holes, as depicted in Figure 3c. 4. Rolling Operation As is shown by Czuczai et al.,7 the model summarized in section 3 works fine for small-scale problems. However, the problem may become unsolvable even for medium-scale processes because of the combinatorial difficulties caused by the need of modeling the connections between units. Therefore, a new decomposition algorithm has been developed, and is presented here. The algorithm also uses the notion of “operation of a unit”. The rolling operation algorithm is similar in some sense to the rolling horizon. The original problem is divided into subproblems. In each subproblem, every unit has a maximum allowed number of operations assigned to. The progress of the algorithm iteration by iteration is depicted in Figure 4. This approach is similar to Wu and Ierapetritou’s approach11 in which the set of event points necessary to satisfy the requested product orders is distributed over several subproblems, and in each iteration the new subproblem gets the fixed production data (fixed task sequence in form of fixed binary variables) from the previous subproblem. The BBT operations described in section 2 accumulate and distribute beer from and to other operations, and they may be considered as central operations of the process. This is the source of the idea to divide the scheduling problem to subproblems which are organized around BBT operations. Subproblems are formed by restricting the number of BBT operations applied to each BBT tank to a very small number, e.g., one or two only. The constraints enforcing to satisfy all the orders are relaxed, and weighted production quantity constrained by the product

orders is first maximized, then the production time necessary for producing that quantity is minimized, finally the order and production data are updated and prepared for the next cycle. The procedure ends when all the original orders are satisfied. The scheme of the new algorithm is shown in Figure 5. The first step of solving a subproblem is maximization of a weighted total satisfaction of orders. The orders are weighted with the inverse of their deadlines. In this way, the orders with earliest deadlines are satisfied earlier. The mathematical form of the objective function maximizing the production quantity is given in eq 8. objV ) WRBTO ·

Vconsrbto - VRAWrbto + TSTAND rbto - TSTART rbto∈RM Vordord - VREQord WLPD · (8) TTARGord - TSTART ord∈ORDS

∑

∑

By maximizing the production quantity, the solution is forced to reach maximum capacity utilization. The makespan is not considered at all, i.e., the operations can take any time; the only requirement is that the outh operation has to be finished before the deadline of order ord if the outh operation produces product for satisfying order ord. The subproblem includes all the constraints detailed in the paper of Czuczai et al.,7 except the two enforcing the satisfaction of the product orders and consumption of raw materials. The relaxation of the latter two constraints is necessary because there are not enough number of available operations for satisfying all the orders and consuming all the raw materials. The second step of solving the subproblem is minimization of the time used for producing the products determined in the first step. There are, however, alternatives for interpreting this production. (a) The stricter interpretation fixes all the intermediate product amounts, i.e., the amounts of intermediates and end products produced in each operation. Then minimization of the total production time gives rise to push all operations to their possible earliest starting moments. (b) A looser interpretation fixes the production amounts of the final products only, and thus leaves some degree of freedom for rearranging the production structure. These alternatives are below referred to as case a and case b, respectively. An even stricter interpretation is fixing all the binary variables in the model, including those denoting connections between operations. This latter version gives rise to an LP problem in the second step. Choosing the objective function of this problem may also have a significant effect on the performance of the algorithm. The best objective function seems to be the sum of the finishing times of all the operations compared to the starting time of the current subproblem. The mathematical form of the objective function minimizing the production time is given in eq 9. This way all of the operations are dragged to left along the time axis. objT )

∑ ∑ ∑ (tf

it u,ou

u

it

- TSTART)

(9)

ou

The production time minimizing subproblem includes the same constraints as its preceding production quantity maximization problem, the only difference is the objective function. After solving the actual subproblem, the production data are compared to the orders. If all the orders are satisfied, the algorithm stops. If there are still unsatisfied orders, the order data are updated; i.e., the resulting production of products is subtracted from the amounts of product orders (and the raw


material consumption is subtracted from the available amounts of raw materials). The weights of the orders and of the raw materials are also updated in each subproblem by subtracting the actual starting time moment of the subproblem from the deadline of the product order or raw material because otherwise the difference between orders with unity difference in deadline would monotonically decrease while proceeding further in the solution process. Since the subproblems may still be large enough to form an obstacle against finding optimum within sufficiently short time, small integrality gap is allowed between the lower and upper bound at the end of the Branch and Bound algorithm. The minimum number of necessary iterations can be estimated from the capacity of the units and the quantity of product orders, using the following equation:

∑

NMINIT )

P max(Fm,i u,it,i

· CAPu,i · NOu,it)

no. of variables

no. of binary variables

no. of constraints

lower bound for makespan (h)

25818

9116

50747

100.14

Table 2. Problem Statistics for the Subproblems with One BBT-Operation/BBT-Tank no. of variables

no. of binary variables

no. of constraints

3953

1299

7504

Table 3. Solution Data of Subproblems with One BBT-Operation, Case a max volume subproblems objective value (kg)

VREQord

ord∈ORDMm

Table 1. Problem Statistics for Example Problem, without Decomposition

(10)

However, it is not required to estimate it accurately; it is enough to give an upper bound which will be sufficient to satisfy the product orders, since the algorithm will automatically stop after complying with all the orders. It is possible, similarly to Basset et al.’s10 approach, to perceive the production quantity maximization subproblem as the production planning model, and the production time minimization subproblem as the scheduling phase. Basset et al.’s10 approach uses an aggregated model for the planning phase, which assigns the necessary product amounts to scheduling subproblems. The present decomposition algorithm differs from others in the literature by using the same model both for assigning the required production quantity to subproblems and for sequencing the tasks. While some literature algorithms10,12 use an aggregate planning model for providing a tight upper bound for the plant capacity, the new approach presented here uses the same rigorous model. Therefore, the plant capacity is not overestimated, and the planning phase does not set an unrealistic production goal which could cause infeasibility in the scheduling phase. Although the performance of the algorithm is tested in the next section using a more or less sequential problem, theoretically it would be possible to use it for scheduling a networkrepresented process in multipurpose plant, by omitting the simplification of the superstructure (i.e., allowing all the technologically reasonable connections between operations). For general network-represented processes, the size of the subproblems probably must be larger in order to ensure feasibility, and the solution would take much longer time. Therefore, it is not straightforward to use the algorithm for those processes in its present form. Future research work will investigate possibilities for extending its validity for any process. 5. Results In this section, the results obtained by applying the rolling operation algorithm, detailed in section 4, to a middle-sized industrial example are presented. The example on which the algorithm is tested is the one detailed at the end of section 2. The problem statistics are shown in Table 1 for the case when it is tried to solve in one step without decomposition. The MILP solver is unable to provide even a feasible solution in reasonable time for that problem. However, it is possible to calculate an absolute lower bound

1 2 3 4 5 6

min time subproblems

solution time (CPU s)

objective value (h)


56.88 33.14 7.70 2.52 2.97 0.24 150.49 157.91

6177.85 7451.42 5754.36 5020.54 3500.56 1043.95

0.16 0.17 0.19 0.19 0.19 0.19

-2534.24 -1474.19 -1130.10 -403.58 -207.43 0 final objective total solution time

for the makespan, using the product order quantity, and production rate data.

[

∑

TLOW ) max u,it

(

∑

i∈{i∈Iu,it|

m∈{m∈MORD| P >0}|)1} |{i∈T|Fm,i

∑

P >0, Fm,i

|{u∈U|CAPu,i>0}|)1}

VREQord

ord∈ORDMm P Fm,i

)]

· min (RATECu1,i1,u,i) u1∈U, i1∈∪Iu1,it1 it1

(11)

Its value for the present example is 100.14 h. The core of the algorithm is written in section 4. The performance of the algorithm may strongly depend on the following characteristics: (1) the size of the subproblems (allowed number of operations); (2) the allowed maximum integrality gap; and (3) the way how the production data are fixed for preparing the minimization of production time (cases a and b). Allowing the maximum integrality gap to be remain wide can speed up the solution process but may lead to worse solution. Fixing some of the binary variables (case a) shortens the solution time but the production time in the solution is longer because the minimization subproblem is less flexible. The value of the objective function may be better if only the total production amount is fixed for each order (case b) but the solution time is much longer in that case. The studied problems are formulated with AIMMS 3. 7,16 and CPLEX 10.017 is applied as MILP solver, using a 1,83 GHz AMD Sempron processor and 768 MB RAM. Table 2 presents the problem statistics for the subproblems while solving the motivating example with one central operation per subproblem. The type and size in each operation determined by the production quantity maximization subproblem are fixed in the time minimization subproblem (case a above), and the solution data are presented in Table 3. The specified maximum


Figure 4. Visualization of the advance of the rolling operation algorithm.

integrality gap is 20% in the production quantity maximization subproblem, and 0% in the production time minimization subproblem. Some binary variables are fixed in the second step, i.e., in minimizing the production time, because existence of operations with fixed amount is not a question. In contrast, all the binary variables are free in the first step, i.e., in maximizing the production quantity. Hence the first step takes much longer time than the second step. The total solution time is 157.9 CPU s; the objective value (the finishing time of the last operation) is 150.5 h. The Gantt chart of the solution is shown in Figure 6; the identifiers of the units listed at the left-hand side in the figure are the same as used in section 2. The bottleneck unit is the packline producing kegs (KegL); this unit has to work almost permanently to achieve a feasible schedule. Applying so many available BBT tanks is useless; and a plant with such an arrangement would be badly designed. Some of the raw beer batches are not used at all; this is due to the fact that their deadlines are farther then the others, and the others are consumed

Figure 6. Gantt chart of solution with one BBT operation, case a.

Figure 5. Rolling operation algorithm.

earlier. This makes the solution even more difficult, as there are raw beer batches not to be used at all. Table 4 presents the solution data for subproblems while solving the motivating example with fixing the total production quantity for each order and consumption of raw materials, but not fixing the production quantity in each operation, between


solution presented in Figure 7 is better; i.e., the makespan is shorter. These are due to the fact that the binary variables of the production time minimization subproblem are all free, and the solution obtained in the production quantity maximization subproblem can be well rearranged to reach shorter production time. The total solution time here is longer by 20% than in case a, but the objective value (115.8 h) is better by 23%.

Table 4. Solution Data of Subproblems with One BBT Operation, Case b max volume subproblems (gap 20%) objective value (kg) -2534.24 -1474.19 -782.07 -444.48 -137.30 0 final objective (h) total solution time (CPU s)

1 2 3 4 5 6

min time subproblems (gap 10%)


objective value (h)


56.70 37.53 2.67 5.08 0.89 0.34 115.77 195.32

6192.57 7209.88 7234.25 3547.02 3459.17 1824.86

2.20 3.52 3.06 12.70 16.48 0.72

Effect of the Specified Integrality Gap. The algorithm was also tested at different integrality gap (also known as relative gap) values to study if the quality of the solution could be improved for exchange of the solution time. Table 5 shows the solution data with (1) 10% gap for the maxV subproblem and 5% gap for the minT subproblem, and (2) 0% gap for both of them. In the latter case, the whole solution process requires almost 2 h, which is hardly reasonable, while the quality of the solution is not improved similarly. This is due to the fact that the used CPLEX solver can find very good integer solution at the beginning of the solution process and it is not necessary to continue searching to proven optimality.

the production quantity maximization and the production time minimization step. That is, results of considering case b are shown in Table 4. The specified maximum integrality gap is 20% in the production quantity maximization subproblems, and 10% in the production time minimization subproblems. Solution of the minimization of production time is much slower than in the previous case, whereas the quality of the

Figure 7. Gantt chart of solution with one BBT operation, case b. Table 5. Solution Data with Different Specified Gap Values, Case b (1) max volume subproblems (gap 10%) 1 2 3 4 5 6

-2534.24 -1521.57 -1034.28 -333.14 -135.36 0 final objective (h) total solution time (CPU s)

116.84 117.53 1.83 4.967 1.19 0.23 112.69 552.28

(2) min time subproblems (gap 5%) 6127.60 6749.04 5168.92 5038.58 2662.27 2580.69

10.97 59.50 62.49 65.39 55.52 0.30

max volume subproblems (gap 0%) -2534.24 225.41 -1474.19 4778.09 -998.74 23.17 -315.29 25.16 -129.86 0.83 0 0.25 final objective (h) total solution time (CPU s)

min time subproblems (gap 0%) 6102.57 6819.83 5212.97 5159.35 2959.01 2651.02 107.42 6787.11

208.09 122.48 1046.30 32.08 269.89 0.72


Figure 8. Gantt chart of solution found with local search. Table 6. Solution Data with Tightened Deadlines, with Cases a and b 1 BBT operation, case a (gapmaxV 20%, gapminT 0%) problem deadline deadline deadline deadline

9h 12 h 15 h 18 h

solution feasibility


feasible infeasible infeasible infeasible

372.36 395.70 283.41 231.71

1 BBT operation, case b (gapmaxV 20%, gapminT 10%) solution feasibility

solution feasibility


202.36 455.44 208.72 300.14

feasible feasible feasible infeasible

4922.86 3938.82 3814.45 2958.69

Table 7. Problem Statistics for the Subproblems with Two BBT-Operations/BBT-Tank no. of binary variables

no. of constraints

7537

2598

14622

2 BBT operations, case b (solution time limit 3600 s)


feasible feasible feasible infeasible

no. of variables

1 BBT operation, case b (gapmaxV 0%, gapminT 0%)

Performance of the model and the algorithm are also studied with further tightened deadlines. Four different problems are created, in which all the product order and raw material deadlines are brought forward with 9, 12, 15, and 18 h, respectively. Table 6 shows the solution feasibility and solution times with cases a and b with different integrality gaps and increased number of BBT operations in the four different problems. Case b provides results with better reliability than case a; therefore, it is worthwhile to use this approach when tight deadlines are present. Effect of the Number of Operations. The algorithm was also tested with increased number of operations, two BBToperations/BBT-tank instead of one. Table 7 shows model statistics for the subproblems. The subproblems are, not surprisingly, two times larger that in the previous cases with one BBT-operation/BBT-tank.

solution feasibility feasible feasible infeasible feasible

solution time (CPU s) 11962.57 11187.51 13097.39 13906.28

Table 8. Solution Data of the Local Search starting solution (h)

radius of the environment (-)

running time (CPU s)

found best objective function (h)

107.42

10

100001.20

104.52

The solution times are steeply increased due to the enlarged complexity, and, therefore, the zero integrality gap could not be reached even for this decomposed problem. 3600 CPU s was applied as solution time limit. The results of different deadlines were compared to the ones obtained with one BBT-operation/ BBT-tank. Now it is possible to find a feasible solution even for the tightest problem (deadline 18 h). However, as it could be expected, the reliability of the algorithm is worse than with 1 BBT-operation, due to the large remaining gap; and the solution times are hardly acceptable. Improving the Solution with Local Search. It is probable that the found solution is not globally optimal, but when tight deadlines are present then the number of existing feasible integer solutions is not large, and thus the found feasible solutions are close to the optimum. This closeness can be expressed using the distance between the binary vectors. Therefore, searching for better solution in a small environment of the current solution can possibly lead to the optimum. The distance between the


binary vectors in the present case can be restricted by the following inequality:

∑

∑

(PYj - yj) +

j,pyj)1

(yj - PYj) e MAXR

(12)

j,pyj)0

where yj are the binary variables of the model, pyj are parameter values of the binary variables of the current solution, and MAXR is the allowed radius of the search. In the present case, only binary variables denoting the execution of a task in an operation are taken into account for calculating the distance between two solutions:

∑∑ ∑ ∑

∑

(PYTu,it,ou,i - ytu,it,ou,i) +

u∈U it∈TT i∈Iu,it ou∈OUu,it pytu,it,ou,i)1

∑∑ ∑ ∑

∑

IOCIN ⊆ T ) tasks which may have inlet transport from only one operation M ) materials MINi ⊆M ) materials that are used in task i MOUTi ⊆M ) materials that are produced in task i MORD ⊆M ) materials that are ordered (i.e., products) O ) operations ORDS ) product orders ORDMm ⊆ ORDS ) product orders belonging to material m OUu,it ⊆ O ) operations of type it of unit u (u ∈ U, it ∈ TT) RM ⊆ M ) raw materials T ) tasks TT ) task types () {BATCH, CONT}) U ) units (including storage tanks) Indices

(ytu,it,ou,i - PYTu,it,ou,i) e

u∈U it∈TT i∈Iu,it ou∈OUu,it pytu,it,ou,i)0

MAXR (13) Adding the above inequality to the model and solving it, it is possible to find better solution than the current one in its arbitrary environment. Table 8 shows the results of this search, and Figure 8 depicts the Gantt chart of the found best solution. As the absolute lower bound of the makespan for the whole problem is 100.14 h, the relative gap between the lower bound and the best found solution is only 4.2%. As can also be seen in the Gantt chart, the packline packing kegs (KegL) is almost constantly working, getting close to the theoretical optimum. 6. Conclusions A new rolling operation decomposition algorithm is developed based on the MILP model of Czuczai et al.7 The new algorithm divides the original process scheduling problem to several subproblems by constraining the maximum number of allowed operations, so that their sizes are considerably smaller. In each subproblem, the sum of weighted production quantity is first maximized, so that the orders with earlier deadlines are satisfied first, and then the production quantity is fixed and the makespan is minimized. The algorithm stops when all the product orders are satisfied and/or all the raw materials are consumed. Application of the new algorithm is demonstrated on a middle-sized industrial example. The algorithm reliably provides us with good feasible solution even in the case of very tight scheduling problems. The example problem is solved using different settings: with different methods for transferring data between the subproblems, with different maximum allowed gaps, and with increased size of subproblems. Effect of tightening the deadlines is also studied. By using a local search method for improving the solution, it is possible to reach a very modest integrality gap between the lower bound and the best found solution, which shows that the optimal solution is well approached. Acknowledgment This research has been supported by OTKA K062099. Notation Sets Iu,it ⊆ T ) tasks of type it of unit u (u ⊆ U, it ⊆ TT)

i,i1,i2,i3 ) for tasks it, it1,it2,it3 ) for task types j ) index of general binary variable m ) for materials ord ) for orders ou,ou1,ou2,ou3 ) for operations rbt ) for raw beer tanks rbto ) for raw materials u,u1,u2,u3 ) for units Parameters CAPu,i ) capacity of unit u in task i (kg) NMINIT ) minimum number of iterations (-) PYj ) values of general binary variables in the current solution prior local search (-) PYTu,it,ou,i ) values of binary variables in the current solution prior local search (-) C Fm,i ) coefficients of material m in the raw material of task i (-) P Fm,i ) coefficients of material m in the product of task i (-) TLOW ) lower bound of makespan (h) TSTANDrbto ) target deadline of raw material rbto (h) TSTART ) starting time of the current subproblem (h) THOR ) length of scheduling time horizon (h) TTARGord ) target deadline of order ord (h) VRAWrbto ) amount of raw material rbto (kg) VREQord ) product demand of order ord (kg) WRBTO ) weight of the raw materials in the production quantity maximization subproblem (-) WLPD ) weight of the product orders in the production quantity maximization subproblem (-) ZIIi1,i2 ) binary parameter denoting the connections between tasks in the production recipe (-) ZOCu1,it1,ou1,u2,it2,ou2 ) binary parameter denoting the allowed connections between operations (-) Variables objT ) objective function of the production time minimization subproblem (h) objV ) objective function of the production quantity maximization subproblem (kg) it tfu,ou |(u ∈U,it ∈TT,ou ∈OUu,it) ) finishing time of outh operation of type it of unit u (h) it tsu,ou |(u ∈U,it ∈TT,ou ∈OUu,it) ) starting time of outh operation of type it of unit u (h) it1,it2 Vcu1,ou1,i1,u2,ou2,i2,m |(u1,u2 ∈U,it1,it2 ∈TT,i1 ∈Iu1,it1,i2 ∈Iu2,it2,m ∈MOUTi1∩MINi2,ZOCu1,it1,ou1,u2,it2,ou2 ) 1,ZIIi1,i2 ) 1) ) material flow of material m from ou1th operation of type it1 of unit u1 to ou2th operation of type it2 of unit u2 (kg)

3908 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 Vconsrbto|(rbto∈RM) ) amount of material consumed of raw material rbto in the current subproblem (kg) Vordord|(ord∈ORDS) ) amount of material produced for satisfying order ord in the current subproblem (kg) it Vtu,ou,i |(u∈U,it∈TT,ou∈OUu,it,i∈Iu,it) ) size of task i processed in outh operation of type it of unit u (kg) Binary Variables it1,it2 ycu1,ou1,i1,u2,ou2,i2 |(u1,u2∈U,it1,it2∈TT,ou1∈OUu1,it1,ou2∈OUu2,it2, i1∈Iu1,it1,i2∈Iu2,it2,ZOCu1,it1,ou1,u2,it2,ou2)1,ZIIi1,i2)1) ) 1, if the material produced in task i1 in ou1th operation of type it1 of unit u1 is used by task i2 in ou2th operation of type it2 of unit u2, 0 otherwise. If unit u2 cannot receive material from unit u1, then that binary variable is not created. it ytu,ou,i |(u∈U,it∈TT,ou∈OUu,it,i∈Iu,it) ) 1, if the task i is processed in outh operation of type it of unit u, 0 otherwise. If unit u cannot process task i, then that binary variable is not created.

Literature Cited (1) Mendez, C. A.; Cerda, J.; Grossmann, I. E.; Harjunkoski, I.; Fahl, M. State-of-the-Art Review of Optimization Methods for Short-Term Scheduling of Batch Processes. Comput. Chem. Eng. 2006, 30, 913–946. (2) Kondili, E.; Pantelides, C. C.; Sargent, R. W. H. A General Algorithm for Short-Term Scheduling of Batch Operations-1. MILP Formulation. Comput. Chem. Eng. 1993, 17, 211–227. (3) Maravelias, C. T. Mixed-Time Representation for State-Task Network Models. Ind. Eng. Chem. Res. 2005, 44, 9129–9145. (4) Sundaramoorthy, A.; Karimi, I. A. A Simpler Better Slot-Based Continuous-Time Formulation for Short-Term Scheduling in Multipurpose Batch Plants. Chem. Eng. Sci. 2005, 60, 2679–2702. (5) Ierapetritou, M. G.; Floudas, C. A. Effective Continuous-Time Formulation for Short-Term Scheduling. 2. Continuous and semicontinuous processes. Ind. Eng. Chem. Res. 1998, 37, 4360–4374.

(6) Janak, S. L.; Lin, X.; Floudas, C. A. Enhanced Continuous-Time Unit-Specific Event-Based Formulation for Short-Term Scheduling of Multipurpose Batch Processes: Resource Constraints and Mixed Storage Policies. Ind. Eng. Chem. Res. 2004, 43, 2516–2533. (7) Czuczai, B.; Farkas, T.; Rev, E.; Lelkes, Z. New MILP Model for Solving Scheduling Problems with Special Characteristics. Submitted for publication in Ind. Eng. Chem. Res. 2008. (8) Mendez, C. A.; Cerda, J. An MILP-Based Approach to the ShortTerm Scheduling of Make-and-Pack Continuous Production Plants. OR Spectrum 2002, 24, 403–429. (9) Pan, M.; Qian, Y.; Li, X. Novel Precedence-Based and Heuristic Approach for Short-Term Scheduling of Multipurpose Batch Plants. Chem. Eng. Sci. 2008, 63, 4313–4332. (10) Bassett, M. H.; Pekny, J. F.; Reklaitis, G. V. Decomposition Techniques for the Solution of Large-Scale Scheduling Problems. AIChE J. 1996, 42, 3373–3387. (11) Wu, D.; Ierapetritou, M. G. Decomposition Approaches for the Efficient Solution of Short-Term Scheduling Problems. Comput. Chem. Eng. 2003, 27, 1261–1276. (12) Janak, S. L.; Floudas, C. A.; Kallrath, J.; Vormbrock, N. Production Scheduling of a Large-Scale Industrial Batch Plant I. Short-Term and Medium-Term Scheduling. Ind. Eng. Chem. Res. 2006, 45, 8234–8252. (13) Kelly, J. D.; Mann, J. L. Flowsheet Decomposition Heuristic for Scheduling: a Relax-and-Fix Method. Comput. Chem. Eng. 2004, 28, 2193– 2200. (14) Harjunkoski, I.; Grossmann, I. E. A Decomposition Approach for the Scheduling of a Steel Plant Production. Comput. Chem. Eng. 2001, 25, 1647–1660. (15) Prasad, P.; Maravelias, C. T.; Kelly, J. Optimization of Aluminum Smelter Casthouse Operations. Ind. Eng. Chem. Res. 2006, 45, 7603–7617. (16) Bisschop, J.; Roelofs, M. AIMMS-The Language Reference, 2007. (17) ILOG CPLEX 10.0 User’s Manual, 2006.

ReceiVed for reView July 4, 2008 ReVised manuscript receiVed January 15, 2009 Accepted January 19, 2009 IE801037K

Rolling Operation Algorithm for Solving Complex ... - ACS Publications

Recommend Documents