Integrated Scheduling and Dynamic Optimization by Stackelberg Game

Article pubs.acs.org/IECR

Integrated Scheduling and Dynamic Optimization by Stackelberg Game: Bilevel Model Formulation and Efficient Solution Algorithm Yunfei Chu and Fengqi You* Department of Chemical and Biological Engineering, Northwestern University, Evanston, Illinois 60208 United States ABSTRACT: We propose a novel method to solve the integrated scheduling and dynamic optimization problem for sequential batch processes. The scheduling problem and the dynamic optimization problems are collaborated by a Stackelberg game (leader−followers game). Mathematically, the integrated problem is formulated into a bilevel program. The scheduling problem in the upper level acts as the leader, while the dynamic optimization problems in the lower level are the followers. The follower problems have their own objectives, but the leader problem can coordinate the follower problems to pursue its objective. To efficiently solve the bilevel program, we develop a decomposition algorithm. It first solves the lower-level problems to determine the response functions. The response functions are then represented by piecewise linear functions to solve the upper-level problem. The integrated method is consistent with the ISA 95 standard and can be easily implemented in an IT infrastructure following the standard. The performance of the proposed method is demonstrated in a complicated batch plant in comparison with the traditional sequential method which solves the scheduling problem and the dynamic optimization problems separately, and the simultaneous method which formulates the integrated problem into a monolithic model. dynamic optimization (MIDO)21,22 or mixed-logic dynamic optimization (MLDO)23,24 problem. Solving the integrated problem is much more computationally demanding than solving the subproblems separately. Even though some efficient algorithms based on Lagrangean decomposition25 or generalized Benders decomposition26,27 have been developed to reduce the computational burdens, the integrated methods still have difficulties in practical implementation. The monolithic model results in a centralized method which solves a plant-wide optimization problem including different products, tasks, and processing units. Implementation of the centralized method needs to dismantle the current hierarchical organization in most manufacturing enterprises. Besides the inconsistency with the current hierarchical manufacturing system, the monolithic integrated methods also lack flexibility to deal with uncertainties. The entire integrated problem needs to be resolved even if only local uncertainties occur in a single unit. Nevertheless, reoptimizing the single unit may be adequate to deal with the minor uncertainties instead of optimizing the whole plant, which is much more computationally challenging. Owing to these difficulties, application of a monolithic solution in industry is still limited.28 The objective of this work is to propose a novel framework for solving the integrated scheduling and dynamic optimization problem that addresses the aforementioned challenges. To facilitate its implementation in the current manufacturing hierarchy, we formulate the integrated problem into a distributed optimization problem instead of the centralized optimization problem. In the distributed problem, the decisionmaking is distributed among different members which optimize

1. INTRODUCTION With the impact of enterprise-wide optimization,1−3 an increasing number of methods have been presented in the literature to integrate different decision levels in the manufacturing hierarchy.4−6 The integrated methods are proven to substantially increase the production efficiency and reduce the production cost by collaboratively optimizing different decision-makings in a manufacturing process. The achieved improvement in productivity is critical for a manufacturing enterprise to stay ahead of escalating global competitions that transform once high-value specialty products into low-value commodities with lower profit margins. Because the margins can be low, every last bit of performance has to be squeezed out by a collaborative optimization approach. In the framework of enterprise-wide optimization, integrated scheduling and dynamic optimization has attracted a wide spectrum of attentions in the recent decade. The integrated methods can significantly improve the overall performance of the entire process compared to the traditional method which solves the dynamic optimization problem and the scheduling problem in a sequential way.7−9 The sequential method first determines the dynamic process for executing a task in a unit. The determined dynamic trajectories of inputs, outputs, and states return the task operational recipes that are imported to the scheduling problem as fixed parameters. Then the scheduling problem is solved with the fixed task recipes. By contrast, the integrated methods allow flexible task operations. The task operational recipes can be optimized simultaneously with scheduling decisions. The advantages of the integrated optimization have been demonstrated in various processes, including continuous processes10−14 and batch processes.15−20 Though the integrated methods outperform the sequential method in terms of the solution quality, the integrated methods encounter severe computational complexity. Most existing methods formulate the integrated problem into a monolithic model by appending the dynamic models as constraints to the scheduling model, which leads to a complicated mixed-integer © 2014 American Chemical Society

Received: Revised: Accepted: Published: 5564

December 17, 2013 February 25, 2014 March 5, 2014 March 5, 2014 dx.doi.org/10.1021/ie404272t | Ind. Eng. Chem. Res. 2014, 53, 5564−5581

Industrial & Engineering Chemistry Research

Article

Figure 1. Structures of three methods for solving scheduling and dynamic optimization problem.

flexible so that the dynamic optimization problems can be solved collaboratively with the scheduling problem. In this work we consider the sequential batch process where the scheduling problem is linked to dynamic optimization problems via task processing times and processing costs. These variables comprise the task operational recipes and also the response functions from the perspective of the Stackelberg game. The scheduling problem (leader) influences the dynamic optimization problems (followers) by setting the task processing times (inputs). The dynamic optimization problems (followers) are then solved to return the processing costs (responses) to the scheduling problem (leader). By using a piecewise linear function to represent the response functions of the processing costs dependent on the processing times, we develop an efficient decomposition algorithm to solve the bilevel integrated problem. The bilevel problem formulation and the efficient solution algorithm lead to an easy implementation strategy, which is consistent with the ISA 95 standard, a popular international standard for process integration in a manufacturing enterprise. The ISA 95 standard defines a multilevel manufacturing system which the bilevel integrated method can be readily embedded in. Specifically, the scheduling problem is solved in level 3 of manufacturing operations management and the dynamic optimization problems are solved in level 2 of automatic control. The efficient solution algorithm can be implemented with an IT infrastructure based on the ISA 95 standard where each subproblem is solved in the distributed computer system and the collaboration is achieved via communications through the network system. Actually, the bilevel method provides a general framework that can unify the sequential method and the simultaneous method as special cases. The sequential method represents a special case where all task recipes are fixed at the response functions. The simultaneous method represents another special case where the objective functions of the followers are the processing costs. The performance of the bilevel method is demonstrated in case studies on a batch plant. The bilevel method returns higher profits than the sequential method in all cases. In the special case where the objective functions of the dynamic optimization problems are the processing costs, the

their own objective functions, while all decisions are made by a single central controller in the centralized problem.29 The formulated problem has a bilevel structure consistent with the manufacturing hierarchy. Specifically, the scheduling problem is solved in the upper level while each dynamic optimization problem for a task execution is solved in the lower level. Every subproblem has its own objective. The scheduling problem and the dynamic optimization problems are collaborated by a leader−followers game (Stackelberg game).30 The scheduling problem acting as the leader is able to coordinate the dynamic optimization problems, which act as the followers, to optimize its objective. Figure 1 visualizes the distinction of the proposed bilevel method from the traditional sequential method and the existing integrated method, also known as the simultaneous method, which solves the monolithic problem. Similar to the sequential method, the bilevel method has a hierarchical structure where the dynamic optimization problems are autonomous, having their own objectives. The autonomy is, however, absent in the simultaneous method where the dynamic models are simply appended to the scheduling problem as constraints. Without their own objectives, the dynamic models are optimized to pursue the objective of the scheduling problem. Yet, the scheduling problem might have some objective other than the economic measure;31 the dynamic optimization problem also needs to achieve some operational objectives regarding safety, quality, operability, and so on. For the monolithic model with a single objective function, it is difficult to consider all different objectives presenting in a complex process. If only the economic objective is optimized, a solution with higher profit may be also a solution difficult for execution. The execution difficulties can result in frequent delays, interruptions, or failures in the task operation process, which greatly offsets the economic advantage.32 These economic losses can be more severe than the ones resulting from a noneconomically optimal, yet safe, production operation.33 Different from the simultaneous method, the bilevel method is a distributed optimization method. The distributed nature is the same as the sequential method. However, the bilevel method enables the collaboration between the scheduling problem and the dynamic optimization problems. Unlike the sequential method, the task recipes are 5565

dx.doi.org/10.1021/ie404272t | Ind. Eng. Chem. Res. 2014, 53, 5564−5581


Article

bilevel method returns the same profit with the simultaneous method but the computational time is much shorter. More importantly, the bilevel method is flexible enough to consider different objective functions of the dynamic optimization problems, which are ignored in the simultaneous method based on the monolithic model. The novelties of this work are summarized as • a new formulation of the integrated problem into a bilevel program • an efficient solution algorithm by piecewise linear representation of the response functions which link the scheduling problem to the dynamic optimization problems • an easy implementation strategy consistent with the ISA 95 enterprise integration standard The remainder of the paper is organized as follows. The problem statement is given in section 2. The integrated method is proposed in section 3, including the model formulation, the solution algorithm, and the implementation strategy. Compared to the sequential method and the simultaneous method, the proposed method is demonstrated in case studies in section 4. The conclusion is given in section 5.

Table 1. Integrated Problem production model sequential batch process where the batch sizes are fixed given (scheduling problem) order demands and due dates order value functions job execution stages (tasks) capable units for executing a task (dynamic optimization problems) objective for operating a task dynamic models initial conditions final values constraints on initial conditions, final values, and trajectories determine task assignment and sequencing task starting times and ending times task processing times and processing costs job completion times trajectories of dynamic systems objective in upper level to maximize production profit

2. PROBLEM STATEMENT This work focuses on the sequential batch process,34−37 where the batch size of executing a task is a fixed parameter. The scheduling problem and the dynamic optimization problems are linked by task processing times and processing costs, which are variables in the integrated problem. Besides these variables, the scheduling problem mainly determines the task sequencing and assignments along with task starting times. A dynamic optimization problem determines the time-dependent trajectories of inputs, states, and outputs for the dynamic system executing a task in a processing unit. The dynamic trajectories are also related to the processing time and the processing cost of the task. Because the batch sizes are fixed in the sequential batch process, the order demand of a product can be quantified by the number of batches. After a batch of materials is executed in a processing unit, the whole batch is transformed to the following unit without material splitting and mixing. The entire procedure for processing a batch of materials through different units is called a “job”.38 Each job has multiple operational stages, called “tasks”.38 When a job is completed, a batch of products is manufactured. Thus, the number of batches required to fulfill the order demand of a product is equal to the number of jobs corresponding to the product. Without loss of any generality, the order demand can be quantified by the number of jobs. Specifically, the integrated problem is stated as shown in Table 1.

and the solution algorithm is developed in section 3.2. The implementation strategy according to the ISA 95 standard is presented in section 3.3. The relationship of the proposed method with the sequential method and the simultaneous method is revealed in section 3.4. 3.1. Model Formulation. We consider the sequential batch process illustrated in Figure 2. A job indexed by j specifies the operational stages through which a final product is manufactured from raw materials. The operational stages of a job are indexed by k. The combination of (j, k) indicates a task that can be executed in a unit indexed by i. The scheduling problem assigns a task to a capable unit, sequences the assigned tasks in every unit, and determines each task starting time. The scheduling results are visualized by a Gantt chart.38 The execution procedure of a task can be represented by a dynamic system. Manipulating the process input trajectories can change the dynamic behavior of the task execution system. The task processing time PTijk and the processing cost PCijk become variables that link the dynamic model for executing task (j, k) in unit i to the scheduling problem. The two variables depend on the dynamic trajectories, which are determined by solving the dynamic optimization problem. The integrated problem is built by combining the dynamic optimization problems with the scheduling problem. Because the task processing times are variables, a continuous-time scheduling model is more convenient than a discrete-time scheduling model.39 In this work, we adopt the general precedence model40−42 to represent the scheduling problem. The scheduling problem is formulated as a mixed-integer linear program (MILP). The dynamic model for executing a task is often represented by a set of differential and algebraic equations. The differential equations can be discretized by a collocation method.43 After the discretization, the dynamic model only consists of algebraic equations and the dynamic optimization problem becomes a nonlinear program (NLP). The detailed formulation of the scheduling model and the discretized dynamic models can be seen in section A.1 of the Appendix. In this work, we focus on the bilevel structure of the integrated problem instead of the specific expressions for the

3. BILEVEL INTEGRATED METHOD In this section, we propose a novel integrated method which formulates the integrated problem into a bilevel program. The bilevel program retains the decentralized nature with the sequential method; meanwhile, it introduces a collaboration procedure which can optimize the entire batch process like the simultaneous method. It will be seen that the proposed method is general enough to absorb the sequential method and the simultaneous method as special cases. This section is organized as follows. The bilevel model formulation is given in section 3.1, 5566



Article

processing time for task (j, k) executed in unit i. Except for the processing times {PTijk} and the processing costs {PCijk}, all scheduling variables are stacked into the vector XS. Similarly, all variables of dynamic model ijk are stacked into the vector XDijk. In the bilevel formulation, the leader problem (Leader_P) is the scheduling problem with variable processing times and processing costs. The objective function given in eq 1 is linear as we consider the linear scheduling problem. Because the objective function is the production profit, it is decreasing as the processing costs increase; i.e., bijk > 0, for all (j, k) ∈ JK, i ∈ Ijk. The constraints of the scheduling model are indexed by m and expressed by inequalities 2. In a common scheduling problem, the constraints are often posed on task assignments, sequences, starting times, and processing times.38 Thus, the task processing costs are absent from constraints 2. The bilevel program has only one leader problem (Leader_P) but multiple follower problems (Follower_Pijk). A follower problem is a dynamic optimization problem on the task execution system. The follower problem is indexed by ijk, which is the combination of the task index (j, k) and the unit index i. Each follower problem has its own objective function expressed in eq 3. The constraints of dynamic model ijk are expressed by inequalities 4 indexed by n. In the scheduling problem (Leader_P), we need to determine {PTijk} and {PCijk}. However, each pair of PTijk and PCijk is related and their relationship is determined by the solution of the dynamic optimization problem (Follower_Pijk). Therefore, the scheduling problem needs to consider the results of the dynamic optimization problems. As an autonomous decision maker, each follower can determine its objective function 3 from the local perspective. The objective function evaluates the performance of the dynamic system from the perspective of safety, operability, profitability, etc. It is common in practice that the dynamic systems of different tasks are optimized by different objectives. Even for the same task the objective function can be frequently changed according to the shifted requirement. The followers have no obligation to cater for the leader’s objective, which is an economic criterion. Of course, the followers can also choose their objectives in order to merely pursue the leader’s goal. 3.2. Solution Algorithm. A general bilevel program is challenging to solve, especially for a multifollower problem.30 Fortunately, in the integrated problem, a follower problem is linked to the leader problem only via the processing time and the processing cost. The interface leads to a straightforward collaboration mechanism between the leader problem and the follower problems. The leader of the scheduling problem assigns a processing time to a follower of the dynamic optimization problem. Then the dynamic optimization problem is solved to return the processing cost to the scheduling problem. When the scheduling problem knows how each dynamic optimization problem responds to a given processing time, it can adjust the processing times to coordinate the dynamic optimization problems so that its objective of the production profit is maximized. No matter what the objective functions of the followers are, the leader can influence the followers through the linking variables and it has complete knowledge about the responses of the followers to the inputs it imposes on them. Therefore, by manipulating the linking variables, the leader can guide the followers to optimize its objective function. This idea motivates the following solution algorithm.

Figure 2. Diagram of the sequential batch process for integrated scheduling and dynamic optimization.

subproblems, so that we represent the integrated model in a general form as (Leader_P) max f (X S) +

∑ aijk PTijk − ∑ bijk PCijk ijk

ijk

(1)

∀m

(2)

s.t. gm(X S , {PTijk }) ≤ 0,

(Follower_Pijk) D min φijk (Xijk , PTijk , PCijk ),

∀ (j , k) ∈ JK, i ∈ Ijk (3)

s.t. D ηijkn(Xijk , PTijk , PCijk ) ≤ 0,

∀n

(4)

For clarity, we use curly braces to denote the collection of variables over their indices. For example, {PTijk} is the collection of all processing times for all (j, k) ∈ JK, i ∈ Ijk. Set JK contains all tasks and set Ijk denotes the units capable for executing task (j, k). A notation without the curly braces only represents a single variable. For example, PTijk indicates a single 5567



Article

we represent the response function by a piecewise linear function. Specifically, for dynamic model ijk, we first obtain the lower bound of the processing time by solving the problem (MinPT_Pijk). Then we generate some grid points starting from the minimum processing time as

The response function of a follower can be implicitly defined by solving the dynamic optimization (Response_Pijk) D γijk(ptvijk) = arg min φijk (Xijk , PTijk , PCijk ), PCijk

min ptsijkq = ptijk + q·dt ijk ,

∀ (j , k) ∈ JK, i ∈ Ijk

(5)

∀n

PTijk = ptvijk

(6) (7)

In eq 5, the response function is denoted by γijk(ptvijk), where ptvijk is the input processing time given by the leader. The value of ptvijk affects the results of the dynamic optimization problem via equality 7. The value of γijk(ptvijk) is equal to the processing cost in the optimal solution corresponding to ptvijk. After the response functions for all followers are calculated, they can be appended to the scheduling problem to replace the follower problems in the lower level. Then the bilevel integrated problem becomes a single level problem given as follows: (IntSingle_P) max f (X S) +

ijk

(8)

∀m

(9)

PCijk = γijk(PTijk ), min PTijk ≥ ptijk ,

∀ (j , k) ∈ JK, i ∈ Ijk ∀ (j , k) ∈ JK, i ∈ Ijk

∀ q ≥ 1, (j , k) ∈ JK, i ∈ Ijk

+ pcsijk(q − 1),

(10)

(16)

Each linear function on the right-hand side is a straight line passing through two adjacent points (ptsijk(q−1), pcsijk(q−1)) and (ptsijkq, pcsijkq). The cutting plane constraint 16 relies on the convexity of the response function PCjk = γjk(PTjk). If the response function is not convex, it can be expressed by a piecewise linear function as

(11)

The follower problems in the lower level are replaced by their response functions in equalities 10, which are defined implicitly by solving the response problems (Response_Pijk). Because a response problem is a dynamic optimization problem, the processing time is equal to the duration for the dynamic system to transit from one state to another state. For example, when the dynamic system represents a reaction task, the state variables include the species concentrations. The processing time is the reaction time for the concentrations, starting from the initial values, to reach the specified values. Of course, the reaction time in this particular example cannot be arbitrarily small. In general, the processing time is bounded from below for a practical system. The minimum processing time is min denoted by ptijk , which enters the optimization problem (IntSingle_P) via inequalities 11. The values of ptmin ijk for all (j, k) ∈ JK, i ∈ Ijk, are determined by solving the following problems:

NT

PTjk =

∑ DTjkq,

∀ (j , k) ∈ JK (17)

q=1 NT

PCjk =

∑

pcsjkq − pcsjk(q − 1) ptsjkq − ptsjk(q − 1)

q=1

NT

DTjkq +

∑ pcsjk(q− 1)λjkq , q=1

∀ (j , k) ∈ JK

(18)

DTjkq ≥ λjkq ptsjk(q − 1), DTjkq ≤ λjkq ptsjkq ,

(MinPT_Pijk) min ptijk = min PTijk ,

(15)

The relaxation does not change the optimal solution. The equalities are automatically enforced because maximizing the objective function 8 pushes the processing costs to their lower bounds given that the objective function is increasing as the processing costs decrease. After the relaxation, inequalities 15 can be represented by a set of cutting plane constraints as pcsijkq − pcsijk(q − 1) PCijk ≥ (PTijk − ptsijk(q − 1)) ptsijkq − ptsijk(q − 1)

s.t. gm(X S , {PTijk }) ≤ 0,

∀ (j , k) ∈ JK, i ∈ Ijk

PCijk ≥ γijk(PTijk ),


(14)

where ptsijkq is a grid point. The time step is denoted by dtijk, and the number of grid points is nT + 1. At each grid point, the response problem (Response_Pijk) is solved with ptvijk = ptsijkq. Solving the optimization problem returns the associated processing costs, denoted by pcsijkq; i.e., pcsijkq = γijk(ptsijkq). In this way, we obtain a set of discrete points {ptsijkq, pcsijkq}, q = 0, ..., nT, from the response function PCijk = γijk(PTijk). Then the response function can be represented by a piecewise linear function linking the grid points {ptsijkq, pcsijkq}, q = 0, ..., NT. However, instead of directly using piecewise linear functions to represent the response functions, we can further simplify the equalities 10 and relax them by inequalities:


q = 0, ..., n T

∀ q ≥ 1, (j , k) ∈ JK ∀ q ≥ 1, (j , k) ∈ JK

(19) (20)

NT

∀ (j , k) ∈ JK, i ∈ Ijk

∑ λjkq = 1,

(12)

∀ (j , k) ∈ JK (21)

q=1


∀n

λjkq ∈ {0, 1}, (13)

λjk(q = 0) = 0,

Though the integrated problem (IntSingle_P) is a common single-level problem, it cannot be solved directly because the response functions in equalities 10 have no closed-form expressions. To address the difficulty of the black-box functions,

∀ q , (j , k) ∈ JK ∀ (j , k) ∈ JK

(22) (23)

The piecewise linear function is formulated by the multiple choice model.44 Other models can also be used. When the response function is convex, the piecewise linear function provides an 5568



Article

overestimator which is the boundary of the feasible range defined by the cutting plane constraint 16. When the objective function is maximized, the cost components will be pushed to the boundary automatically. Therefore, the cutting plane constraint is equivalent to the piecewise linear function formulation. In this case, the cutting plane constraint is much simpler than the piecewise linear function without introducing extra binary variables. The convexity of the response function can be verified by the scatter plot of the grid points {ptsijkq, pcsijkq}, q = 0, ..., NT. Admittedly, the cutting plane constraint 16 is an approximation of the response function in inequality 15. The approximation aims to represent the response function in a closed-form expression so that the bilevel program is transformed into a singlelevel program that is much easier to solve. Considering even a linear bilevel program with only continuous variables and a single follower problem is intrinsically difficult (NP-hard),45 the bilevel integrated problem is mathematically intractable for an exact method due to the presence of a great number of follower problems, binary variables, and nonlinear constrains. Some heuristics are necessary to tackle such a complicated problem. Thus, we approximate the response functions in this work. As the response functions are one-dimensional functions, each of which is the processing cost dependent on a single variable of the processing time. The dimension does not depend on the number of follower problems, the size of the leader problem, or the size of the follower problems. Thus, the cutting plane constraint 16 can be applied to a large-scale integrated problem. In summary, by replacing the response function constraints 10 using the cutting plane constraints 16, the single-level integrated problem becomes

Each dynamic optimization problem for generating a recipe data point ptsijkq and pcsijkq is a nonlinear program, which is optimized independently with others. The minimum processing times ptmin ijk are obtained by solving the problems (MinPT_Pijk). The flowchart of the solution algorithm is visualized in Figure 3.

Figure 3. Flowchart of the solution algorithm.

(IntMILP_P) max f (X S) +

Solving the problem (IntMILP_P) determines the optimal task processing times, denoted by {PT*ijk}. They can be different from the grid values {ptsijkg}. To determine the dynamic trajectories according to PTijk * , the response problem (Response_Pijk) is solved with ptvijk = PT*ijk. The determined trajectories are sent to the control systems as references. The detailed implementation strategy will be presented in section 3.3. 3.3. Implementation Strategy. Nowadays, the manufacturing process of an industrial company is often organized as a hierarchy. Many international standards have been issued to model and implement the manufacturing hierarchy. A widely used standard is ANSI/ISA 95, or ISA 95 as it is more commonly referred to, which has been developed for global manufacturers.9 It is applied in all industries, and in all sorts of processes, including batch processes. The bilevel formulation and the decomposition algorithm facilitate the implementation of the integrated method in an ISA 95 based system. The ISA 95 standard provides consistent terminologies, models, and interfaces to integrate different levels of the manufacturing hierarchy. Many ISA 95 based systems have been implemented in real-world processes, including IT organizations and infrastructure layouts. In the hieratical system based on the ISA 95 standard, the integrated scheduling and dynamic optimization method is implemented in level 2 of the automatic control and level 3 of the manufacturing execution. The ISA 95 standard has been introduced in two review papers7,9 for integrated scheduling and dynamic optimization. Unfortunately, the implementation issue according to the standard is neglected in previous research, resulting in a gap between the theory and the practice. The main reason is that the previous methods presented in the literature are applied to the monolithic model which is inconsistent with the


ijk

s.t. gm(X S , {PTijk }) ≤ 0, PCijk ≥

∀m

pcsijkq − pcsijk(q − 1) ptsijkq − ptsijk(q − 1)

+ pcsijk(q − 1), min PTijk ≥ ptijk ,

(PTijk − ptsijk(q − 1))

∀ q ≥ 1, (j , k) ∈ JK, i ∈ Ijk

∀ (j , k) ∈ JK

This is an MILP problem and can be solved efficiently by a state-of-the-art MILP solver such as CPLEX. The recipe data {ptsijkq} and {pcsijkq} are generated by solving a set of separable optimization problems: (Recipe_Pijkq) D pcsijkq = arg min φijk (Xijk , PTijk , PCijk ), PCijk

∀ q , (j , k) ∈ JK, i ∈ Ijk


∀n

PTijk = ptsijkq min ptsijkq = ptijk + q·dt ijk

5569



Article

Figure 4. Implementation of the integrated method in an IT infrastructure according to the ISA 95 standard.

Table 2. Implementation Procedure of Figure 4 step 1 2 3 4 5 6 7 8

The lower-level problems (Recipe_Pijkq) are solved in level 2. The recipe data ptsijkq and pcsijkq for executing task (j, k) in unit i are generated by a local computer. The recipe data of different tasks can be generated by different computers using different optimization software and algorithms. The distributed optimization can be easily implemented in the distributed control systems (DCS) that are widely used in the automatic control level. After the recipe data are generated by different computers in level 2, they are transferred to level 3 through the operation information network. All recipe data for every possible combination of the task and the unit are stored in a database server in level 3. The upper-level problem (IntMILP_P) is solved by a computer in level 3. The solution is triggered by order information sent from level 4. The demanded products determine the jobs and the tasks in the scheduling problem. According to the order information, the problem (IntMILP_P) is formulated. The computer in level 3 retrieves the recipe data for the tasks corresponding to the products demanded by the order. Based on the recipe data, the upper-level problem (IntMILP_P) is solved. The solution to the upper-level problem determines the optimal schedule and the optimal task execution recipes. The execution instruction is returned to the computers in level 2. When the computers in level 2 receive the instruction, they obtain the task assignments, starting times, and processing times. Then they execute the tasks according to the instruction. According to the optimal processing times {PT*ijk} determined in level 3, the computers in level 2 solve the response problems (Response_Pijk) with ptvijk = PT*ijk to determine the dynamic trajectories. These trajectories are then sent to automatic control systems as the references. The control systems manipulate the task execution processes tracking the reference trajectories. Once the task execution processes are completed, they return the execution results to the computers in level 2.

hierarchical model stipulated by the standard. By contrast, the bilevel model formulation has a natural hierarchical structure facilitating its implementation in practice. Figure 4 shows the implementation strategy of the integrated method in an IT infrastructure based on the ISA 95 standard. The implementation procedure in Figure 4 includes the steps presented in Table 2. The distributed optimization not only facilitates implementation of the integrated method in the IT infrastructure based on the ISA 95 standard, but also ensures an easy adaption

to uncertainties. When uncertain outcomes are observed, the integrated problem needs to be resolved according to the outcomes. The monolithic integrated method needs to solve the entire problem for a minor local disruptions. However, the bilevel method only needs to solve a part of the integrated problem that is affected by the uncertainties. For example, when a unit breaks down, it is quickly replaced by another unit. The task recipe data for the replacement unit has been stored in the database server in level 3, so that the upper-level problem (IntMILP_P) can be solved immediately with the updated task 5570



Article

recipes. The computers related to other units do not need to solve the lower-level problems (Response_Pijk) again to update the task recipe data. The advantage of solving the integrated problem by the bilevel method will be demonstrated in the case studies. 3.4. Relationship with Sequential Method and Simultaneous Method. The bilevel method actually provides a general framework to solve the integrated problem. We will show that the bilevel method is general enough to incorporate the sequential method and the simultaneous method as special cases. The sequential method first solves the dynamic optimization problems to determine the task recipes. Then the scheduling problem is solved with the fixed task recipes. Specifically, the dynamic optimization problems solved by the sequential method are the response problems (Response_Pijk) with fixed processing times ptvijk. The returned processing costs are pcvijk = γijk(ptvijk). The fixed ptvijk and pcvijk are imported to the scheduling problem (Scheduling_P1) max f (X S) +

The objectives of the dynamic optimization problems are discarded and the dynamic models are appended to the scheduling problem as constraints. The integrated problem is solved to solely maximize the objective function of the scheduling problem. The monolithic problem can be solved in two stages as max

X S ,{PTijk },{PCijk }

= Smax [f (X S) + X ,{PTijk }

= Smax [f (X ) + X ,{PTijk }

ijk

s.t. ∀m

ijk

ijk

(29)


s.t.

∀ (j , k) ∈ JK, i ∈ Ijk

gm(X S , {PTijk }) ≤ 0, PCijk = hijk (PTijk ),

∀m ∀ n , (j , k) ∈ JK, i ∈ Ijk

hijk (PTijk) = min PCijk ,

∀ n , (j , k) ∈ JK, i ∈ Ijk


∀n

The problem (Simultaneous_P2) is equivalent to the problem (IntSingle_P) when the response functions are identical, which implies the objective functions of the follower problems are chosen as

(Simultaneous_P1)


ijk

PCijk ] ∑ aijk PTijk − ∑ bijk min PC

ijk

∀ (j , k) ∈ JK, i ∈ Ijk

ijk

ijk

ijk

max f (X S) +

Clearly, the sequential method is a special case of the bilevel method where the recipe data are fixed and not changed by the integrated problem (IntSingle_P). Due to the fixed recipes, there is a lack of collaboration between the scheduling problem and the dynamic optimization problems, which is the main drawback of the sequential method. The simultaneous method can also be regarded as a special case of the bilevel method. The simultaneous method solves the monolithic model as max f (X S) +

min (∑ bijk PCijk )] ∑ aijk PTijk − {PC }

In the bilevel optimization problem 27, X and {PTijk} are the variables in the outer problem and the processing costs, {PCijk}, are the variables in the inner problem. When the inner problem is solved, the variables in the outer problems are regarded as parameters. Because the first two terms in the inner objective function in problem 27 only depend on the variables of the outer problem, they can be moved out of the inner problem, leading to problem 28. From constraints 25 and 26, we can see that the processing cost PCijk is only related to the processing time PTijk. Thus, minimizing the weighted sum of the processing times in problem 28 is equal to the weighted sum of the minimized processing times in problem 29 given that the inner problem is solved with fixed processing times. Based on the bilevel program, the monolithic problem (Simultaneous_P1) is reformulated into the following problem: (Simultaneous_P2)

∑ aijk PTijk − ∑ bijk PCijk

PCijk = γijk(ptvijk),

ijk

S

(Scheduling_P2)

PTijk = ptvijk ,

ijk

(28)

To manifest the relation with the integrated problem (IntSingle_P), the scheduling problem is reformulated into the following problem:

gm(X S , {PTijk }) ≤ 0,

∑ aijk PTijk − ∑ bijk PCijk)]

ijk

S

∀m

ijk

ijk

(27)

s.t.

max f (X S) +

ijk

X ,{PTijk } {PCijk }

ijk

gm(X S , {ptvijk}) ≤ 0,

∑ aijk PTijk − ∑ bijk PCijk

= Smax [ max (f (X S) +

∑ aijk ptvijk − ∑ bijk pcvijk ijk

f (X S) +

D φijk (Xijk , PTijk , PCijk ) = PCijk ,

(24)

s.t.

∀ (j , k) ∈ JK, i ∈ Ijk (30)

gm(X S , {PTijk }) ≤ 0, D ηijkn(Xijk , PTijk , PCijk ) ≤ 0,

∀m

In this case with the special objective functions of the follower problems stipulated by eq 30, the followers only pursue the economic objective of the leader. The bilevel formulation (Simultaneous_P2) is equivalent to the single-level formulation (Simultaneous_P1) without any simplification.

(25)

∀ n , (j , k) ∈ JK, i ∈ Ijk (26) 5571



Article

be changed by manipulating process input trajectories. The operational recipes of other tasks are fixed parameters. The processing cost is the utility cost calculated from the consumed heating fluid and cooling fluid. All jobs processed by the batch plant are cast into two families. The jobs belonging to the same family have the identical dynamic models describing the reaction tasks and the identical fixed recipe parameters for other tasks. However, the process inputs can be different and the jobs belonging to the same family can have different operational recipes for the reaction tasks, despite the identical dynamic models. The different operational recipes are a distinct feature of an integrated method, as a result of collaboratively optimizing the scheduling process and the dynamic systems. By contrast, the traditional sequential method results in the same operational recipes as the dynamic models are the same. The parameters of the scheduling model and the dynamic models are listed in section A.2 in the Appendix. All optimization problems are modeled by GAMS 24 in a computer with an Intel CPU (3.10 GHz) with four cores and 8 GB of memory. The operating system is Windows 7 Professional (64 bit). The procedure of generating the recipe data by solving the dynamic optimization problems is presented in section 4.1. The performance of the proposed method is compared with the sequential method and the simultaneous method in section 4.2 when the dynamic optimization problems are solved to minimize the processing costs. In the bilevel formulation, the objective functions can be general other than the processing costs. In section 4.3, the performance of the bilevel method is demonstrated for another type of objective functions that trade off between profitability and operability. The efficiency of the proposed method enables an immediate response to the changing production environment, which is demonstrated in section 4.4. The proposed method can be applied to a largescale problem in section 4.5. 4.1. Dynamic Optimization in Lower Level. We study two objective functions for the dynamic optimization problems. The first one is the processing cost in eq 30. In this case, the bilevel method is equivalent to the simultaneous method where the dynamic optimization problems cooperate with the scheduling problem to maximize the production profit. We call it the cooperative case. The other objective function is chosen as

However, the bilevel model in eqs 1−4 is more general than the single-level model provided that the objective functions of the follower problems can be chosen as general ones rather than the economic criteria of the processing costs. In this sense, the single-level problem (Simultaneous_P1) is just a special case (Simultaneous_P2) of the bilevel problem in eqs 1−4. Because the simultaneous method solves the single-level problem (Simultaneous_P1), it can be also regarded as a special method able to solve the particular bilevel problem (Simultaneous_P2). Nevertheless, when the followers’ objective functions are not the processing times, the bilevel model is not reduced to the singlelevel model and the simultaneous method is not applicable. We should note that the bilevel method presented in this work is distinct from other decomposition methods previously developed for solving integrated problems, such as the Lagrange decomposition,25 and the generalized Benders decomposition.26,27 The major difference between these decomposition methods and the proposed one relies on the different models which they are applied to. The existing decomposition methods are applied to the monolithic integrated problem like the simultaneous method, and they are just more efficient alternatives for solving the resultant MINLP problem than the simultaneous method. By contrast, the proposed method is applied to the bilevel model. As a bilevel model is intrinsically different from a single-level model, it is not straightforward to extend the existing decomposition methods to the bilevel model. A particular solution method should be devised as we did in section 3.2. In summary, the goal of the presented framework is not to solve previously formulated problems in a more efficient way. Instead, it presents a novel formulation of the integrated problem by the bilevel program. Compared to the common monolithic model formulation, the bilevel model is consistent with the current production hierarchy in a manufacturing company and is easy to implement in an existing IT infrastructure following the international ISA 95 standard. The formulated bilevel problem can be solved efficiently by the developed decomposition method, which other decomposition methods are difficult to apply to.

4. CASE STUDIES To demonstrate the integrated method, we apply it to a batch plant in this section. The diagram of the batch plant is shown in Figure 5. There are eight processing units in the batch plant. A job consists of five operational stages: mixing, reaction, filtration, reaction, and separation. The task of the first stage is executed in a mixing tank TI. Then the reaction task of the second stage can be executed in either reactor RI or reactor RII. The outlet materials of the reaction go through the filtration task executed in a filter FI in the third stage. The fourth stage is another reaction task, which can be executed in either reactor RIII or reactor RIV. The separation task of the fifth stage can be executed in either separator SI or separator SII. All reaction tasks executed in the four reactors are described by dynamic systems. The process inputs are the flow rate of heating fluid (FH) and the flow rate of cooling fluid (FC) of a reactor jacket. Manipulating the flow rates changes the jacket temperature (TJ) and then the reactor temperature (TR) by heat exchange. The reactor temperature determines the reaction rates, which affect trajectories of species concentrations (CA, CB, CC, and CD). The dynamic models of the reaction tasks are also shown in Figure 5. The operational recipes, i.e., the processing times and the processing costs, of reaction tasks can

D φijk (Xijk , PTijk , PCijk ) = PCijk +

∀ (j , k) ∈ JK, i ∈ Ijk

∑ ΔUijk(PT)2 , (31)

where ΔUijk(T) = Uijk(T + ΔT) − Uijk(T) denotes the difference between the input values at two adjacent time points. Given that the dynamic model is discretized by the collocation method, ΔT is set as the length of the finite element. Besides the processing cost, the objective function includes a penalty term on variations of the input trajectories. In this case, the goal of solving a dynamic optimization problem is not only profitability but also operability. Adding the penalty term aims to smooth the input trajectories, making them easy to be tracked by automatic control systems.46 In practice, an abrupt change in the process input is also not desirable as the dramatic change can impose a huge disruption on the task operation process.47 The penalty term is just an example that dynamic optimization on a unit operation often needs to consider different objectives besides the economic criterion.48 According 5572



Article

Figure 5. The batch plant for case studies.

slightly larger than those in the cooperative case, as shown in Figure 6. However, the input trajectories are much smoother. They gradually increase or decrease, avoiding the abrupt changes present in the cooperative case. We can investigate the quantitative results by comparing point A for the cooperative case with point B for the noncooperative case. Point A associates with the processing cost of 75.7 m.u. (monetary unit), and the sum of squared variations is equal to 481.8 (m3/h)2. At point B, the processing cost slightly increases to 90.4 m.u. However, the sum of squared variations substantially decreases to 42.456 (m3/h)2. Though the economic performance is slightly sacrificed, the difficulty for control systems to track the reference trajectories is significantly reduced. When all recipe functions are generated in the same way as illustrated in Figure 6, they are imported into the upper level problem (IntMILP_P). Then the integrated problem is solved. The recipe functions also provide the processing times and processing costs for the sequential method, which solves the scheduling problem (Scheduling_P2). Unlike the integrated

to the objective functions in eq 31, the dynamic optimization problems do not cater for the objective of the scheduling problem. We call it the noncooperative case. The first step of the bilevel integrated method is to solve the dynamic optimization problems (Recipe_Pijkq) in the lower level. They are solved to determine the recipe data for reaction tasks in the case studies. The minimum processing time of a reaction task is first determined by solving the problem (MinPT_Pijk). From the minimum processing time, 11 time points are generated as ptsijkq, q = 0, ..., 10. The time step is 0.1 h. For each ptsijkq, the corresponding processing cost pcsijkq is determined by solving the problem (Recipe_Pijkq). Linking the 11 points of (ptsijkq, pcsijkq) forms a piecewise linear function that represents the operational recipes of the task. As an example, Figure 6 shows the recipe functions for the reaction task of the first job executed in reactor RI in both cooperative and noncooperative cases. In the noncooperative case, the penalty terms are added to the objective functions 31, so the resulting processing costs are 5573



Article

II. The parameter values can be found in Table 8 in the Appendix. For the three methods, the dynamic models are discretized by a collocation method. The number of finite elements for the discretization is 30, and each finite element contains three collocation points. The dynamic optimization problems are solved using CONOPT on the discretized models. The scheduling problem for the sequential method and the upper-level problem for the bilevel method are solved by CPLEX 12.5 to the zero gap. The MINLP problem for the simultaneous method is solved by SBB to the 1% gap. We also tried DICOPT for solving the MINLP, but it did not return a feasible solution. The MINLP is initialized by the results returned by the sequential method. The model and solution statistics for the three methods are summarized in Table 3. The bilevel method spends 8.1 s in Table 3. Solution and Model Statistics of the Three Methods in the Cooperative Case method bilevel sequential simultaneous solution

dynamic models scheduling or integrated model

CPU (s) gap (%) profit (m.u.) sales (m.u.) variable cost (m.u.) fixed cost (m.u.) equations continuous var type

8.1 0 535.5 1880.0 544.5 800.0 2716 2830 MILP

1.5 0 428.0 1915.0 687.0 800.0 2716 2830 MILP

49 043.7 1 539.9 1881.1 541.2 800.0 − − MINLP

equations continuous var binary var

863 251 132

575 187 132

22 371 22 831 132

solving the integrated problem. The computational time is slightly larger than 1.5 s required by the sequential method. However, the bilevel method returns a much higher production profit than the sequential method by 25.1%, demonstrating the advantage of the integrated method. The profit returned by the bilevel method is nearly the same as that of the simultaneous method. The difference in the profits, which is caused by the piecewise linear approximation of the task recipe functions, is negligible. However, the computational time for the simultaneous method is by 3 orders of magnitude longer, because the simultaneous method needs to solve a complicated MINLP while the bilevel method only solves an MILP upper-level problem and a set of separable dynamic optimization problems in the lower level. The Gantt charts returned by the three methods are shown in Figure 7. The three methods return different schedules and reaction task recipes with variable time lengths, indicated by the lengths of the task bars. The sequential method results in the smallest makespan, which is associated with the largest sales shown in Table 3. However, the variable cost of the sequential method is the largest, causing the lowest profit. By contrast, the two integrated methods significantly reduce the variable cost by slightly increasing the sales. A better trade-off is obtained along with a higher profit than the sequential method. 4.3. Comparison in Noncooperative Case. Besides the economic criterion, a local dynamic optimization problem may have its own objective, such as the smoothness of the input trajectories shown in Figure 6. In the simultaneous method, the local dynamic optimization problem loses its autonomy.

Figure 6. Operational recipes of task (1, 1) (the first stage of job 1) executed in reactor RI.

method, the sequential method does not change the processing times and the processing costs. In the following comparison, we use the second data point in each function (for example, point A or B in Figure 6) to assign the parameters of processing times and processing costs for the sequential method. 4.2. Comparison in Cooperative Case. In the cooperative case, a local dynamic optimization problem minimizes the processing cost on behalf of the scheduling problem to maximize the production profit. It is revealed in section 3.4 that in this special case the simultaneous method is equivalent to the bilevel integrated method. Though they can return the same optimal value, they solve the integrated problem in different ways. The simultaneous method solves the problem (Simultaneous_P1) directly by using an MINLP solver while the bilevel method solves it with the two-stage approach presented in section 3.2. Because an integrated problem often leads to a large-scale MINLP, the simultaneous method could be time-consuming.15,16 Therefore, we use a medium-size problem for the comparison to alleviate the computational burden for the simultaneous method. Nevertheless, we should note that the proposed method can be applied to a large-scale problem which will be illustrated in section 4.5. The problem we solve in this section has four jobs. Job 1 and job 2 belong to family I, while job 3 and job 4 belong to family 5574



Article

problems are not exactly consistent with that of the scheduling problem, the bilevel method can still guide the dynamic optimization problems to pursue the objective function value for the scheduling problem. However, this guided search is absent in the sequential method, where the scheduling problem cannot influence the results of the dynamic optimization problems once they are made. The outcome is that the bilevel method returns a much larger profit than the sequential method by 46%. The Gantt charts of the two methods in the noncooperative case are shown in Figure 8. The task recipes in the noncooperative

Figure 8. Scheduling results returned by the bilevel method and the sequential method in the noncooperative case. A task is represented by (job, stage) shown on the task bar.

Figure 7. Scheduling results by the three methods in the cooperative case. A task is represented by (job, stage) shown on the task bar.

case are different from those in the cooperative case because of different objectives of the dynamic optimization problems. The different task recipes can result in different schedules from those in the cooperative case. 4.4. Resolving Integrated Problem under Uncertainties. Unlike the simultaneous method which solves the monolithic problem, the bilevel method adopts a distributed optimization approach. Under uncertainties, the bilevel method can only optimize some components of the integrated problem which are affected by the uncertain outcomes while other unaffected components can remain unchanged. The distributed optimization mechanism enables the bilevel method to have a fast response to the uncertainties, facilitating its application in practice. We solve the integrated problem in the noncooperative case under two common uncertainties: processing unit breakdown and order demand change. They are two types of uncertainties affecting different components of the integrated problem. The unit breakdown represents a disruption in the dynamic optimization problem on a single unit, while the demand change represents a disruption in the scheduling problem. Before uncertainties occur, the production process is operated according to the normal schedule shown in Figure 8 (top). When

However, its autonomy can be preserved in the sequential method and the bilevel method. We solve the integrated problem by changing the objective function of a dynamic optimization to the sum of the processing cost and the penalty term on the input variations in eq 31. Except for the objective functions for the lower-level problems, the integrated problem is the same as that in section 4.2. The solutions returned by the bilevel method and the sequential method are shown in Table 4. In the noncooperative case, even though the objectives of the dynamic optimization Table 4. Solution Statistics of the Bilevel Method and the Sequential Method in the Noncooperative Case method CPU (s) profit (m.u.) sales (m.u.) variable cost (m.u.) fixed cost (m.u.)

bilevel

sequential

8.6 434.5 1870.0 635.5 800.0

1.7 297.6 1915.0 817.4 800.0 5575



Article

entire problem. The short computational times enable an immediate response to the uncertainties. 4.5. Large-Scale Case. The efficiency of the bilevel method ensures its application to large-scale problems. In this subsection, we extend the four-job problem in section 4.2 to a 10-job problem. In this problem, jobs 1−5 belong to family I and jobs 6−8 belong to family II. The parameter values can be seen in Table 7 in the Appendix. We compare the bilevel method with the sequential method and the simultaneous method. For the sequential method and the bilevel method, the expansion of the scheduling model does not incur the extra computational efforts for the dynamic optimization problems given that the new jobs are the same as the previous ones. The scheduling problem of the sequential method has 3953 equations, 943 variables, and 810 binaries. The upper-level problem of the bilevel method has 4673 equations, 293 continuous variables, and 810 binary variables. Both problems are MILPs and solved by CPLEX 12.5. For the simultaneous method, the expansion of the scheduling model also increases the number of dynamic models contained in the

reactor RI breaks down at 1.5 h, it is replaced by another reactor, RV, which is identical to reactor RI. The schedule under uncertainties is shown in Figure 9 (top). The second task of

Table 5. Solution Statistics of the Bilevel Method and the Sequential Method in the Large-Scale Casea method CPU (s) profit (m.u.) sales (m.u.) variable cost (m.u.) fixed cost (m.u.)

bilevel

sequential

1079.1 942.3 4525.4 1583.1 2000.0

85.3 744.0 4787.4 2043.4 2000.0

a

The simultaneous method does not return a feasible solution in 50 000 s.

Figure 9. Resolving the integrated problem by the bilevel method under uncertainties. The normal case is the one shown in Figure 8 (top).

job 1 is interrupted and the materials are discarded. Job 1 needs to be reproduced from the first stage. Consequently, before the restarting task (1, 1) finishes, the new task (1, 2) cannot start. The reschedule is obtained by solving the integrated problem. Not only the scheduling decisions but also the task operational recipes are changed. The replacement of reactor RI by reactor RV does not affect the dynamic optimization problems on other processing units. These units can use the same recipe data. We only need to replace the recipe data of reactor RI by those of RV. As the recipe data are all determined offline and stored in the database, replacing the recipe data is fast. When the upper-level problem is resolved online with the new recipe data, the solution only requires about 1 s. For the demand uncertainty, the production process is initially operated according to the normal schedule shown in Figure 8 (top). At 1 h, a new order demand is received for producing job 5 which belongs to the same family as job 1 and job 2. To deal with the demand change, we do not need to reoptimize the dynamic models. The recipe data of jobs 1−4 are not changed. We only need to import the recipe data of job 5 and extend the upper-level problem to include the new job. Then the extended problem is resolved online, which requires 12.9 s. The schedule is shown in Figure 9 (bottom). The two problems under uncertainties demonstrate the advantage of the bilevel method. As a distributed optimization approach, it only needs to resolve the integrated problem partially, preventing the time-consuming reoptimization of the

Figure 10. Scheduling results of the 10-job problem by the bilevel method and the sequential method. The numbers on the task bar indicate the job numbers. 5576



Article

indicates the task (j, k) precedes task (j′, k′) in unit i. The constraints in the scheduling model are listed as follows. Each task is executed exactly once:

integrated problem. The monolithic model is a large-scale MINLP that has 58 437 equations, 56 737 continuous variables, and 810 binary variables. The MINLP problem is solved by SBB. In the time limit of 50 000 s (about 13.9 h), the solver did not return a feasible solution. The solution statistics are listed in Table 5. The results of the simultaneous method are not listed as no feasible solution is found. The computational time of the sequential method is 85.3 s, and that of the bilevel method is 1079.1 s. Even for the large-scale integrated problem that the simultaneous method fails to solve, the bilevel method is able to return the optimal solution in a reasonably short period. Compared to the sequential method, the bilevel method returns a higher production profit by 26.7%. The Gantt charts returned by the bilevel method and the sequential method are shown in Figure 10.

∑ ξijk = 1,

(A1)

The relation between the assignment variables and the sequencing variables is ξijk + ξij ′ k ′ − βijkj ′ k ′ − βij ′ k ′ jk ≤ 1, ∀ (j , k) ∈ JK, (j′, k′) ∈ JK, (j , k) ≠ (j′, k′), i ∈ Ijk ∩ Ij ′ k ′

(A2)

βijkj ′ k ′ ≤ ξijk ,

5. CONCLUSION We proposed a novel integrated method based on the leader− followers game. The scheduling problem is the leader, while the dynamic optimization problems are the followers. All players have their own objective functions. The objective function of the scheduling problem is the production profit, while those of the dynamic optimization problems concern operability and profitability. Though the followers are autonomous, they can be coordinated by the leader to pursue its objective. The leader− followers game leads to a bilevel formulation of the integrated problem where the scheduling problem is in the upper level and the dynamic optimization problems are in the lower level. The bilevel program is solved efficiently by a decomposition algorithm. Each dynamic optimization problem is solved separately to determine the response function which specifies the task processing cost as a function of the processing time. Then the dynamic optimization problems are replaced by their response functions represented by piecewise linear functions. The bilevel program is transformed into a single-level program which is the scheduling problem with flexible task recipes characterized by the response functions. The model formulation and the solution algorithm lead to a distributed optimization approach, which can be easily implemented in an IT infrastructure based on the ISA 95 standard. The advantages of the bilevel method were demonstrated in the case studies. In the cooperative case, the bilevel method returned the same production profit as the simultaneous method but the computational time was reduced by more than 3 orders of magnitudes. The computational efficiency enabled its application to a large-scale problem where the simultaneous method failed to find a feasible solution in 50 000 s. Compared to the sequential method, the bilevel method returned larger profits in the cooperative case, noncooperative case, and large-scale case by 25.1, 46.0, and 26.7%, respectively. The bilevel method is a distributed optimization approach that is adaptable to the changing production environment, such as the unit breakdown or the rush change in the order demand. Under uncertainties, only the components affected by the uncertainty outcomes were reoptimized instead of the entire problem by the simultaneous method.

■

∀ (j , k) ∈ JK

i ∈ Ijk

∀ (j , k) ∈ JK, (j′, k′) ∈ JK, (j , k) ≠ (j′, k′), i ∈ Ijk ∩ Ij ′ k ′

(A3)

βijkj ′ k ′ ≤ ξij ′ k ′ , ∀ (j , k) ∈ JK, (j′, k′) ∈ JK, (j , k) ≠ (j′, k′), i ∈ Ijk ∩ Ij ′ k ′

(A4)

The relation between the starting time TSjk and the ending time TEjk is TEjk = TSjk +

∑ XPTijk ,

∀ (j , k) ∈ JK

i ∈ Ijk

(A5)

The expression of bilinear terms (XPTijk = ξijkPTijk) is ∀ (j , k) ∈ JK, i ∈ Ijk

0 ≤ XPTijk ≤ PTijk , max XPTijk ≤ ξijk ptijk ,

∀ (j , k) ∈ JK, i ∈ Ijk

(A6) (A7)

max XPTijk ≥ PTijk − (1 − ξijk)ptijk ,

∀ (j , k) ∈ JK, i ∈ Ijk

(A8)

The execution sequence of tasks belonging to the same job is TSjk ≥ TEj(k − 1),

∀ (j , k) ∈ JK, k ≥ 2

(A9)

The task sequence in a unit is TSjk + PTijk + st i − bm(1 − βijkj ′ k ′) ≤ TSj ′ k ′, ∀ (j , k) ∈ JK, (j′, k′) ∈ JK, (j , k) ≠ (j′, k′), i ∈ Ijk ∩ Ij ′ k ′

(A10)

where sti is the changeover time in unit i and bm is a sufficiently large value. Objective Function of Scheduling Problem. The objective of the scheduling problem is to maximize the production profit defined by

APPENDIX

PR = SA − CO − cf

A.1. Formulation of Integrated Scheduling and Dynamic Optimization Problem

(A11)

where PR is the profit, SA is the sales, CO is the variable cost, and cf is the fixed cost. The sales term is expressed by the sum of the order values OVj:

Scheduling Model. The general precedence model40,41 is adopted to formulate the scheduling problem. There are two major types of binary variables. The assignment variable ξijk = 1 indicates task (j, k) is assigned to unit i. The sequencing variable βijkj′k′ = 1

SA =

∑ OVj j

5577

(A12) dx.doi.org/10.1021/ie404272t | Ind. Eng. Chem. Res. 2014, 53, 5564−5581


Article

The order value OVj is a piecewise linear function49 of the job completion time CTj, which is expressed by ⎞ ⎛ 1 OVj = pjR ⎜⎜1 − max DTIIj ⎟⎟ , c j − dd j ⎠ ⎝

∀ (j , k) ∈ JK, i ∈ Ijk , r > 0, s

(A14)

κjdd j ≤ DTIjs ≤ dd j ,

∀j

(A15)

CTj ≤ c jmax ,

∀j

∀j

⎡ 88 − 7 6 ⎢ 360 ⎢ ⎢ 296 + 169 6 [css ′] = ⎢ ⎢ 1800 ⎢ ⎢ 16 − 6 ⎢⎣ 36

(A16) (A17)

∀j

(A18)

The piecewise linear function consists of a constant segment and a decreasing segment. The two segments are joined at ddj. The binary variable κj = 1 indicates CTj ≥ ddj. The variable cost, CO, is the sum of processing costs, PCijk. CO =

⎡ 16 − 6 [ds] = ⎢ ⎣ 36

∑ XPCijk

max XPCijk ≤ ξijk pcijk ,

∀ (j , k) ∈ JK, i ∈ Ijk ∀ (j , k) ∈ JK, i ∈ Ijk

dTijk

1⎤ ⎥ 9⎦

∀ (j , k) ∈ JK, i ∈ Ijk , r , s

initial conditions:

(A20)

I (r = 0) Xijk = xijk ,

(A21)

∀ (j , k) ∈ JK, i ∈ Ijk

(A26)

final conditions: (r = nijk R )

Xijk (A22)

F = xijk ,

∀ (j , k) ∈ JK, i ∈ Ijk

(A27)

processing times:

Dynamic Models. The dynamic model for executing a task can be expressed by a set of ordinary differential and algebraic equations. dXijk(Tijk)

16 + 6 36

−2 + 3 6 ⎤ ⎥ 225 ⎥ ⎥ −2 − 3 6 ⎥ ⎥ 225 ⎥ 1 ⎥ ⎥⎦ 9

(A25)

max XPCijk ≥ PCijk − (1 − ξijk)pcijk ,

∀ (j , k) ∈ JK, i ∈ Ijk

88 + 7 6 360

16 + 6 36

rs rs ρijk (Xijk , Uijk ) ≤ 0,

The bilinear terms XPCijk = ξijkPCijk can be linearized by the following constraints: 0 ≤ XPCijk ≤ PCijk ,

296 − 169 6 1800

discretized algebraic equations:

(A19)

ijk

(A24)

The length of a finite element is represented by lijk. The coefficients css′ and ds are given by a specific collocation method. In this work, we use the Radau IIA method, which has a robust performance for stiff differential equations. 5 0 The corresponding coefficients are

(A13)

∀j

CTj = TEj(k = nK ),

s

∀j

CTj = DTIj + DTIIj ,

0 ≤ DTIIj ≤ κj(c jmax − dd j),

r r−1 rs rs Xijk = Xijk + lijk ∑ dsfijk (Xijk , Uijk ),

R PTijk = lijknijk ,

∀ (j , k) ∈ JK, i ∈ Ijk

(A28)

Table 6. Parameters of Reactors (Figure 5) = πijk(Xijk(Tijk), Uijk(Tijk)), VR ρR cR VJ ρJ cJ UAJ C0A C0B C0C C0D CFC CFD T0R T0J TFR TH TC FMax H FMax C PRH PRC

∀ (j , k) ∈ JK, i ∈ Ijk 0 ≥ ρijk (Xijk(Tijk), Uijk(Tijk)),

∀ (j , k) ∈ JK, i ∈ Ijk

Each dynamic model is indexed by ijk. Correspondingly, all variables in the differential equations are indexed by ijk, including the time variable Tijk. For a compact expression, the states Xijk(Tijk) and the inputs Uijk(Tijk) are expressed by the vector forms. The vector forms are also used to express other variables and constraints in the dynamic models. To facilitate the numerical solution of the dynamic optimization problem, the differential equations are discretized by the collocation method.43 After the discretization procedure, a continuous-time trajectory is represented by a set of grid points, indexed by finite element r and collocation point s. discretized differential equations: rs r−1 rs ′ rs ′ Xijk = Xijk + lijk ∑ css ′πijk(Xijk , Uijk ), s′

∀ (j , k) ∈ JK, i ∈ Ijk , r > 0, s

(A23) 5578

RI/RII

RIII/RIV

units

5 1 × 103 2.5 1.5 1 × 103 4.186 8 × 104 1 1 0 0 0.9 − 300 300 310 350 300 15 15 5 1

5 9 × 102 3.0 1.5 1 × 103 4.186 7 × 104 − − 0.9 0 − 0.8 300 300 310 350 300 15 15 5 1

m3 kg/m3 kJ/kg·K m3 kg/m3 kJ/kg·K kJ/h·K kmol/m3 kmol/m3 kmol/m3 kmol/m3 kmol/m3 kmol/m3 K K K K K m3/h m3/h m.u./m3 m.u./m3



Article

m = constraint in the scheduling problem n = constraint in the dynamic optimization problems q = grid point for processing times and processing costs r = finite element s = collocation point

processing costs: R

(r = nijk )

PCijk = ψijk(Xijk

),

∀ (j , k) ∈ JK, i ∈ Ijk

(A29)

Integrated Model. The integrated model is summarized as

Sets

max PR (A11 − A22)

JK = set of tasks Ijk = set of units capable for executing task (j, k)

s.t. scheduling model (A1 − A10)

Parameters

aijk = coefficient before PTijk in the objective function of the scheduling problem bijk = coefficient before PCijk in the objective function of the scheduling problem bm = sufficiently large number Cmax = upper bound on the completion time of job j j css′ = collocation matrix cf = fixed cost ds = collocation vector ddj = parameter in the order value function for job j dtijk = time step for discretizing PTijk lijk = length of the finite element for dynamic optimization problem ijk nT = number of grid points nRijk = number of finite elements for dynamic optimization problem ijk pRj = constant in the order value function for job j pcsijkq = qth discrete point of PCijk prH = unit cost of heating fluid prC = unit cost of cooling fluid pcmax ijk = upper bound of PCijk ptmax ijk = upper bound of PTijk ptmin ijk = lower bound of PTijk ptsijkq = qth discrete point of PTijk ptvijk = fixed processing time for dynamic optimization problem ijk sti = changeover time in unit i xFijk = final value of Xijk xIijk = initial value of Xijk

dynamic models (A23 − A29) A.2. Data of Case Studies

Tables 6, 7, and 8 list the parameters of the reactors and reaction tasks as well as scheduling parameters for the case studies. Table 7. Parameters of Reaction Tasks task family I ER1 k1 ΔH1 ER2 k2 ΔH2

II

1 × 10 5 × 103 −1 × 104 1 × 1010 8 × 103 −5 × 104 7

units

1 × 10 6 × 103 −2 × 104 5 × 108 7 × 103 −4 × 104 8

K m3/kmol·h kJ/kmol K 1/h kJ/kmol

Table 8. Parameters of Scheduling Problem

■

parameter

value

units

selling price of job j fixed cost of job j due date of job j horizon for completing job j fixed processing time of mixing tasks fixed processing time of filtration tasks fixed processing time of separation tasks changeover time in mixing tank changeover time in reactors changeover time in filter changeover time in separators

500 200 15 20 1.0 1.0 2.5 0.2 0.5 0.3 0.4

m.u. m.u. h h h h h h h h h

Binary Variables

ξijk = equal to 1 when task (j, k) is assigned to unit i βijkj′k′ = equal to 1 when task (j,k) precedes task (j′,k′) in unit i κj = equal to 1 if CTj ≥ ddj

AUTHOR INFORMATION

Continuous Variables

Corresponding Author

CO = variable cost CTj = completion time of job j DTIj = first component of CTj DTIIj = second component of CTj OVj = order value for completing job j PCijk = processing cost of task (j, k) executed in unit i PR = production profit PTijk = processing time of task (j, k) executed in unit i SA = sales Tijk = time variable of dynamic model ijk TEjk = ending time of task (j, k) TSjk = starting time of task (j, k) Uijk = input vector of dynamic model ijk Ursijk = value of Uijk at collocation point s of finite element r ΔUijk = change in Uijk Xijk = state vector of dynamic model ijk Xrijk = value of Xijk at the beginning of finite element r Xrsijk = value of Xijk at collocation point s of finite element r

*E-mail: [email protected]. Tel.: (847) 467-2943. Fax: (847) 491-3728. Notes

The authors declare no competing financial interest.

■

NOMENCLATURE

Abbreviations

MILP = mixed-integer linear program MINLP = mixed-integer nonlinear program MIDO = mixed-integer dynamic optimization MLDO = mixed-logic dynamic optimization NLP = nonlinear program Indices

i, i′ = processing unit j, j′ = job k, k′ = operational stage (j, k) = task 5579



Article

XS = vector containing all variables in the scheduling problem except the processing times and the processing costs XDijk = vector containing all variables in the dynamic model for task (j, k) executed in unit i XPCijk = product of ξijk and PCijk XPTijk = product of ξijk and PTijk

(15) Mishra, B. V.; Mayer, E.; Raisch, J.; Kienle, A. Short-term scheduling of batch processes. A comparative study of different approaches. Ind. Eng. Chem. Res. 2005, 44 (11), 4022−4034. (16) Nie, Y. S.; Biegler, L. T.; Wassick, J. M. Integrated scheduling and dynamic optimization of batch processes using state equipment networks. AIChE J. 2012, 58 (11), 3416−3432. (17) Chu, Y.; You, F. Integrated Scheduling and Dynamic Optimization of Sequential Batch Processes with Online Implementation. AIChE J. 2013, 59 (7), 2379−2406. (18) Capon-Garcia, E.; Guillen-Gosalbez, G.; Espuna, A. Integrating process dynamics within batch process scheduling via mixed-integer dynamic optimization. Chem. Eng. Sci. 2013, 102, 139−150. (19) Chu, Y.; You, F. Integration of Scheduling and Dynamic Optimization of Batch Processes under Uncertainty: Two-Stage Stochastic Programming Approach and Enhanced Generalized Benders Decomposition Algorithm. Ind. Eng. Chem. Res. 2013, 52 (47), 16851−16869. (20) Chu, Y.; You, F. Moving horizon approach of integrating scheduling and control for sequential batch processes. AIChE J. 2014, DOI: 10.1002/aic.14359. (21) Barton, P. I.; Allgor, R. J.; Feehery, W. F.; Galan, S. Dynamic optimization in a discontinuous world. Ind. Eng. Chem. Res. 1998, 37 (3), 966−981. (22) Bansal, V.; Sakizlis, V.; Ross, R.; Perkins, J. D.; Pistikopoulos, E. N. New algorithms for mixed-integer dynamic optimization. Comput. Chem. Eng. 2003, 27 (5), 647−668. (23) Busch, J.; Oldenburg, J.; Santos, M.; Cruse, A.; Marquardt, W. Dynamic predictive scheduling of operational strategies for continuous processes using mixed-logic dynamic optimization. Comput. Chem. Eng. 2007, 31 (5−6), 574−587. (24) Oldenburg, J.; Marquardt, W.; Heinz, D.; Leineweber, D. B. Mixed-logic dynamic optimization applied to batch distillation process design. AIChE J. 2003, 49 (11), 2900−2917. (25) Terrazas-Moreno, S.; Flores-Tlacuahuac, A.; Grossmann, I. E. Lagrangean heuristic for the scheduling and control of polymerization reactors. AIChE J. 2008, 54 (1), 163−182. (26) Chu, Y.; You, F. Integration of production scheduling and dynamic optimization for multi-product CSTRs: Generalized Benders decomposition coupled with global mixed-integer fractional programming. Comput. Chem. Eng. 2013, 58 (0), 315−333. (27) Chu, Y. F.; You, F. Q. Integrated Scheduling and Dynamic Optimization of Complex Batch Processes with General Network Structure Using a Generalized Benders Decomposition Approach. Ind. Eng. Chem. Res. 2013, 52 (23), 7867−7885. (28) Kadam, J. V.; Marquardt, W. Integration of economical optimization and control for intentionally transient process operation. In Assessment and Future Directions of Nonlinear Model Predictive Control; Findeisen, R., Allgower, F., Biegler, L. T., Eds.; SpringerVerlag: Berlin, 2007; Vol. 358, pp 419−434. (29) Li, X. H.; Wang, Q. N. Coordination mechanisms of supply chain systems. Eur. J. Oper. Res. 2007, 179 (1), 1−16. (30) Colson, B.; Marcotte, P.; Savard, G. An overview of bilevel optimization. Ann. Oper. Res. 2007, 153 (1), 235−256. (31) Yue, D. J.; You, F. Q. Sustainable scheduling of batch processes under economic and environmental criteria with MINLP models and algorithms. Comput. Chem. Eng. 2013, 54, 44−59. (32) Shobrys, D. E.; White, D. C. Planning, scheduling and control systems: why cannot they work together. Comput. Chem. Eng. 2002, 26 (2), 149−160. (33) Tatjewski, P. Advanced Control of Industrial Processes: Structures and Algorithms. In Advanced Control of Industrial Processes: Structures and Algorithms; Springer: New York, 2007; pp 1−332. (34) Maravelias, C. T. General framework and modeling approach classification for chemical production scheduling. AIChE J. 2012, 58 (6), 1812−1828. (35) Mendez, C. A.; Cerda, J.; Grossmann, I. E.; Harjunkoski, I.; Fahl, M. State-of-the-art review of optimization methods for shortterm scheduling of batch processes. Comput. Chem. Eng. 2006, 30 (6− 7), 913−946.

Functions

f = objective function of the scheduling problem gm = function in constraint m of the scheduling problem ηijkn = function in constraint n of dynamic optimization problem ijk γijk = response function of dynamic optimization problem ijk φijk = objective function of dynamic optimization problem ijk ψijk = cost function in dynamic model ijk πijk = function in the state equation of dynamic model ijk ρijk = function in the path constraint of dynamic model ijk

■

REFERENCES

(1) Grossmann, I. Enterprise-wide optimization: A new frontier in process systems engineering. AIChE J. 2005, 51 (7), 1846−1857. (2) Varma, V. A.; Reklaitis, G. V.; Blau, G. E.; Pekny, J. F. Enterprisewide modeling & optimizationAn overview of emerging research challenges and opportunities. Comput. Chem. Eng. 2007, 31 (5−6), 692−711. (3) Wassick, J. M.; Agarwal, A.; Akiya, N.; Ferrio, J.; Bury, S.; You, F. Q. Addressing the operational challenges in the development, manufacture, and supply of advanced materials and performance products. Comput. Chem. Eng. 2012, 47, 157−169. (4) Harjunkoski, I.; Nystroem, R.; Horch, A. Integration of scheduling and controlTheory or practice? Comput. Chem. Eng. 2009, 33 (12), 1909−1918. (5) Chu, Y.; You, F.; Wassick, J. M.; Agarwal, A. Integrated planning and scheduling under production uncertainties: Bi-level model formulation and hybrid solution method. Comput. Chem. Eng. 2014, DOI: 10.1016/j.compchemeng.2014.02.023. (6) Yue, D. J.; You, F. Q. Planning and Scheduling of Flexible Process Networks Under Uncertainty with Stochastic Inventory: MINLP Models and Algorithm. AIChE J. 2013, 59 (5), 1511−1532. (7) Engell, S.; Harjunkoski, I. Optimal operation: Scheduling, advanced control and their integration. Comput. Chem. Eng. 2012, 47, 121−133. (8) Munoz, E.; Capon-Garcia, E.; Moreno-Benito, M.; Espuna, A.; Puigjaner, L. Scheduling and control decision-making under an integrated information environment. Comput. Chem. Eng. 2011, 35 (5), 774−786. (9) Prata, A.; Oldenburg, J.; Kroll, A.; Marquardt, W. Integrated scheduling and dynamic optimization of grade transitions for a continuous polymerization reactor. Comput. Chem. Eng. 2008, 32 (3), 463−476. (10) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and control of a multiproduct CSTR. Ind. Eng. Chem. Res. 2006, 45 (20), 6698−6712. (11) Nystrom, R. H.; Franke, R.; Harjunkoski, I.; Kroll, A. Production campaign planning including grade transition sequencing and dynamic optimization. Comput. Chem. Eng. 2005, 29 (10), 2163−2179. (12) Zhuge, J.; Ierapetritou, M. G. Integration of scheduling and control with closed loop implementation. Ind. Eng. Chem. Res. 2012, 51 (25), 8550−8565. (13) Chu, Y.; You, F. Integration of scheduling and control with online closed-loop implementation: Fast computational strategy and large-scale global optimization algorithm. Comput. Chem. Eng. 2012, 47, 248−268. (14) Terrazas-Moreno, S.; Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and optimal control of polymerization reactors. AIChE J. 2007, 53 (9), 2301−2315. 5580



Article

(36) Chu, Y. F.; Wassick, J. M.; You, F. Q. Efficient scheduling method of complex batch processes with general network structure via agent-based modeling. AIChE J. 2013, 59 (8), 2884−2906. (37) Chu, Y. F.; You, F. Q.; Wassick, J. M. Hybrid method integrating agent-based modeling and heuristic tree search for scheduling of complex batch processes. Comput. Chem. Eng. 2014, 60, 277−296. (38) Pinedo, M. L. Scheduling: Theory, Algorithms, and Systems; Springer: New York, 2012. (39) Floudas, C. A.; Lin, X. X. Continuous-time versus discrete-time approaches for scheduling of chemical processes: a review. Comput. Chem. Eng. 2004, 28 (11), 2109−2129. (40) Mendez, C. A.; Cerda, J. Dynamic scheduling in multiproduct batch plants. Comput. Chem. Eng. 2003, 27 (8−9), 1247−1259. (41) Pan, C. H. A study of integer programming formulations for scheduling problems. Int. J. Syst. Sci. 1997, 28 (1), 33−41. (42) Demir, Y.; Isleyen, S. K. Evaluation of mathematical models for flexible job-shop scheduling problems. Appl. Math. Modell. 2013, 37 (3), 977−988. (43) Cuthrell, J. E.; Biegler, L. T. On the optimization of differentialalgebraic process systems. AIChE J. 1987, 33 (8), 1257−1270. (44) Croxton, K. L.; Gendron, B.; Magnanti, T. L. A comparison of mixed-integer programming models for nonconvex piecewise linear cost minimization problems. Manage. Sci. 2003, 49 (9), 1268−1273. (45) Hansen, P.; Jaumard, B.; Savard, G. New branch-and-bound rules for linear bilevel programming. SIAM J. Sci. Stat. Comput. 1992, 13 (5), 1194−1217. (46) Srinivasan, B.; Palanki, S.; Bonvin, D. Dynamic optimization of batch processesI. Characterization of the nominal solution. Comput. Chem. Eng. 2003, 27 (1), 1−26. (47) Aziz, N.; Mujtaba, I. M. Optimal operation policies in batch reactors. Chem. Eng. J. 2002, 85 (2−3), 313−325. (48) Bonvin, D. Optimal operation of batch reactorsa personal view. J. Process Control 1998, 8 (5−6), 355−368. (49) Wassick, J. M.; Ferrio, J. Extending the resource task network for industrial applications. Comput. Chem. Eng. 2011, 35 (10), 2124−2140. (50) Hairer, E.; Wanner, G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems; Springer: New York, 2004; Vol. 2.

5581


Integrated Scheduling and Dynamic Optimization by Stackelberg Game

Recommend Documents