An Improved Energy-Awareness Formulation for General Precedence

Article pubs.acs.org/IECR

An Improved Energy-Awareness Formulation for General Precedence Continuous-Time Scheduling Models Hubert Hadera,†,‡ Rachid Labrik,†,§ Guido Sand,† Sebastian Engell,‡ and Iiro Harjunkoski*,† †

ABB Corporate Research, Wallstadter Strasse 59, 68526 Ladenburg, Germany Technical University of Dortmund, Emil-Figge-Strasse 70, 44221 Dortmund, Germany § KTH Royal Institute of Technology, Brinellvägen 8, 11428 Stockholm, Sweden ‡

ABSTRACT: Efficient frameworks for scheduling of energy-intensive industries are important enablers for the industrial implementation of the demand-side management concept (DSM). DSM enables enterprises to cut their energy cost as well as to support the grid in dealing with fluctuating availability of energy. In this work, we present an improvement of the energy-aware production scheduling strategy reported in Hadera et al.,1 which is based upon a continuous-time scheduling formulation. The improvement concerns the modeling of the coupling between the scheduling problem and the computation of the energy bill. The resulting optimization problem is solved by the heuristic bilevel decomposition from Hadera et al.2 Numerical studies showed that although the proposed formulation and heuristic decomposition does not guarantee to obtain the global minimum it does provide good quality solutions within reasonable solution times. The new energy-awareness formulation reduces the number of binary variables and constraints and results in shorter solution times.

1. INTRODUCTION Energy has become a critical resource for most industries, in particular for the process industries. Taking its cost into consideration when optimizing the schedule of a production process can lead to significant cost savings in the electricity bill by exploiting the fluctuating price of electricity and the structure of the procurement contracts. The fluctuations of the price of electricity are caused by several factors. The increasing generation from renewables is associated with a high level of uncertainty and unpredictability, especially if electricity from renewables can be fed into the grid without limitations and at fixed prices, as in Germany. New markets have been created to deal with the demand and supply of electricity with the goal to arrive at prices that reflect the market situation. This changes the business model of the utilities and the energy purchasing policies of the consumers. Apart from this, the demand of electricity is increasing in many regions, which contributes to volatility of the availability and price of electricity. Therefore, adjusting the use of electric energy, for example, to take advantage of off-peak tariffs, can lead to a better control of the cost of electricity. Consequently, production planning should be done including the knowledge about the pricing of electricity. The same technique can also be used to reduce the CO2 footprint of the production by the use of “green” electricity when it is available and reducing consumption in a period with a high share of generation from fossil sources. The implementation of strategies to reduce the consumption of electricity and the energy bill is called industrial demand-side management (iDSM). It has two components: The first one is to improve energy efficiency, which is mainly achieved by technological improvements and often related to the design or to the use of particular pieces of equipment. The second component is demand response management and can be achieved by intelligent scheduling of the production operations, © 2016 American Chemical Society

considering information on the availability and the price of electricity. This technique has the advantage of low upfront cost because it does not require significant investments into hardware. Optimizing the production planning regarding the electricity prices in this domain can bring benefits both for the plant and for the supplier as noted by many industrial and academic studies.3−5 In contrast, production scheduling under constraints on the consumption of energy or coupled with optimizing the energy bill creates big challenges for research as well as for industrial solution providers, as pointed out in the review by Harjunkoski et al.6 One major challenge with regard to optimal scheduling in the process industries is scalability (i.e., the applicability of approaches to scheduling of real-world problems). This field has gained considerable attention in the past 40 years. Proposed solutions range from dispatching rules over bottleneck heuristics to sophisticated parallel branch and bound algorithms. Rigorous solutions of production planning and scheduling problems are based upon mixed integer linear programming as presented, for example, by Pochet and Wolsey.7 Méndez et al.8 presented a classification of scheduling problems for batch processes and point out the features that characterize the corresponding optimization models. These problems involve many on−off decisions (binary variables) and many critical constraints, which often makes the problems large and complicated to solve. Large-scale problems arise especially for batch processes with multiple stages and multiple parallel machines at each stage. An example for this is the optimal scheduling of stainless-steel plants, which has been recognized Received: Revised: Accepted: Published: 1336

September 1, 2015 January 8, 2016 January 14, 2016 January 14, 2016 DOI: 10.1021/acs.iecr.5b03239 Ind. Eng. Chem. Res. 2016, 55, 1336−1346

Article

Industrial & Engineering Chemistry Research

Figure 1. Plant layout and product flow.

as a very challenging industrial scheduling problem9 and is also considered in this contribution. An important issue in energy-price-aware scheduling is how to deal with multiple complicated energy contracts. Most studies assume single fixed consumption versus price curves. Castro et al.10 developed an RTN-based discrete-time model for dealing with a single price curve and load commitment problem of a steel plant. The model discretizes the time horizon into time intervals and faces computational problems when the discretization is done with the exactness that is desired by the industry for scheduling problems, for example, using 10 min intervals. Another discrete-time approach to steel plant scheduling has been proposed by Ashok,11 where time-ofuse (TOU) pricing is considered. There are also several studies that investigate discrete-time approaches for air separation plants. For example, Karwan and Keblis12 discuss air separation plant processes with a focus on different electricity tariff designs and plant utilizations. Recently, Mitra et al.13 developed a discrete-time MILP model to address production planning in air separation plants under time-sensitive electricity prices for continuous power-intensive processes and designed a crossdecomposition scheme for solving large instances.14 Continuous-time models do not suffer from the blow-up of the number of variables with increasing precision of the timing of the operations and have been used successfully for steel plant scheduling (e.g., Harjunkoski and Grossmann9). In contrast, it is challenging to track the consumption of resources on a fixed discrete grid as it arises from electricity pricing (e.g., in 1 h intervals). Consequently, research has been done recently on how to embed energy-awareness in precedence-based continuous-time formulations of the production scheduling problem. Nolde and Morari15 developed a model with six binaries per task that captures how much electricity a given production task is consuming during prespecified price-related time intervals. Hadera and Harjunkoski16 expanded this model to parallel machines and the concept of stages on the basis of the scheduling study by Harjunkoski and Sand.17 Castro et al.18 improved the strategy how to account for electricity consumption by formulating tighter constraints that were derived from a general-disjunctive programming (GDP) formulation. Another approach for solving the stainless-steel production scheduling problem with energy cost awareness was proposed by Hait and Artigues.19 They introduced two binaries that when combined with constraint programming (CP) led to a heuristic for solving a problem with a single price curve and operator availability constraints. Other continuous-time approaches were studied for example by Castro et al.,20 where an RTN-based approach was used to respond to a single price

curve and to schedule cement production. In a follow up study,21 the computational performance limitations were addressed by a rolling-horizon approach. The scheduling models that are developed in this paper are targeting energy-price-aware scheduling of large industrial processes. As in Hadera and Harjunkoski16 and Hadera et al.,1 we tackle the stainless-steel production process that poses the challenges of large-scale scheduling formulations, embedding energy-consumption constraints into continuous-time scheduling formulations and dealing with multiple timesensitive energy contracts as well as with the option of generating and/or selling electric energy to the market. We present an improved monolithic scheduling and energy procurement optimization model where the number of binary variables is reduced. However, the resulting problems are not solvable within reasonable computation times for industrial scale problems. To overcome this limitation, the combined scheduling and energy purchase and generation optimization problem is solved using a heuristic bilevel decomposition,2 which provides satisfactory solutions in reasonable computation times. In the remaining sections of the paper, the industrial case is first described in section 2, including the energy procurement optimization problem. In section 3, the proposed improved formulation for coupling energy-awareness with continuoustime scheduling models is presented. The resulting formulation is then used in numerical experiments (section 4) both for the solution of a monolithic model and for a heuristic bilevel scheme for the evaluation of the benefits of the more compact formulation. On the basis of the numerical experiments, conclusions are drawn, and suggestions for further research are made in section 5.

2. PROBLEM DESCRIPTION In this section, an energy-intensive industrial process, the melt shop of a stainless-steel plant, is described. The process consumes a large amount of electric energy, the cost of which is included into the optimization. The energy bill consists of the cost of electricity purchase and of load deviation penalties. The bill can be reduced by generating electric power and by selling a surplus of electricity on the market in times where this is advantageous. 2.1. Steel Plant Use Case and Scheduling Model. In the melt shop, batches of material are produced that are called heats. There are four production stages as shown in Figure 1. A heat goes through all of the four process stages, and the process is considered a batch process with a semicontinuous last stage. The production starts with melting scrap metal in energy1337

DOI: 10.1021/acs.iecr.5b03239 Ind. Eng. Chem. Res. 2016, 55, 1336−1346

Article


contracts and from technical constraints (e.g., the installed generation capacity and the restrictions during ramp-up). An important aspect of energy management to consider is the load deviation problem. The plant must commit ahead of the scheduling horizon to consume certain amounts of electricity in the different contracts which may differ on an hourly basis. If the actual consumption differs from the preagreed values by more than an assumed margin, then the plant is obliged to pay financial penalties. Both the provider and the plant benefit from this. The provider knows in advance a good approximation of the load levels to be balanced with the supply, which leads to minimization of his operating cost. In return, the consumer (the steel plant) gets a considerable reduction of the price of electricity from the provider. Often the load deviation penalty is related to one single contract with a specific pricing scheme such as time-of-use. The plant is assumed to have the opportunity to sell electricity on the day-ahead spot market in order to reduce the net bill. This electricity may come from generation but also from purchase under other contracts. All the contracts as well as the selling activity have specific thresholds of the maximum amounts that can be bought or sold at each hour as further explained in section 4.

intensive electric arc furnaces (EAF). Then, chemistry adjustments and reduction of carbon content is done at the argon oxygen decarburization (AOD) stage. For a heat to be cast, it needs to be first prepared at the ladle furnace (LF) to ensure proper temperature levels. Lastly, at the continuous casting (CC) stage the heats are processed to form steel slabs. Complex operational constraints must be satisfied when producing a given number of heats of specified material by the end of a given scheduling horizon. This includes the fact that at the CC stage a specified number of heats must be processed in a specified sequence without interruptions. There are nonidentical parallel units (called machines) at each stage: Each product must be processed at each stage, and the defined order of flow must be respected. The model includes machine specific setup times, minimum transportation times that are dependent on the units that are used, and product- and stagespecific maximum hold-up times. The scheduling problem is to find an optimal production schedule for the next 24 h that satisfies numerous operational constraints and at the same time minimizes the sum of the completion times of all tasks of the overall schedule and the cost of energy with a suitable weighting factor. The plant needs to fulfill the given orders; hence, the number of products that must be delivered during the scheduling horizon is fixed. 2.2. Energy Optimization Problem. The energy demand of the production process described in section 2.1 must be satisfied at all times. The plant is assumed to have access to multiple sources of electricity (Figure 2):

3. MODEL FORMULATION To address the problem described in the precious section, an MILP model is set up that consists of several components, similar to Hadera et al.2 3.1. Monolithic Formulation. The scheduling part of the model describes the process rules and constraints using a general precedence approach with assignment and sequencing variables. The continuous-time representation is used because of its exactness of the starting times. The scheduling model is extended by the energy-awareness strategy that models the consumption of electricity in the time intervals related to energy pricing and accounts for load deviation penalties. Once the electricity consumption is computed, it is employed in the minimum-cost flow network problem for the optimization of the purchase and sale of electricity. The load deviation constraints and cost are formulated using the knowledge of the consumption of electricity. The first component, the scheduling part, is identical to the scheduling formulation provided in Hadera et al.,1 and here, only the key components are discussed for completeness sake. The model assigns products p to machines m using binary variables Xm,p and sequences them per stage st with binary variables Vst,p,p′ providing the scheduling component to the objective function (summation of the completion times of the tasks). The model ensures that all production constraints are satisfied; this concerns in particular the continuous casting of heats p within heat groups hg at the final CC stage st4. The variable waiting time wp,st that a heat spends in between processing on two stages is bounded by a minimum min transportation time tm,m′ and maximum hold-up time tmax p,st . Before the process can start, a machine specific setup time tsetup m is respected. The second part of the monolithic model, the energy cost awareness extension, is linked to the scheduling formulation via s the variables of task start times tm,p and provides the information on how much electricity the process consumes (expressed by a continuous variable qs) in a particular time interval s related to pricing cs,i,j and penalties for over- and underconsumption (cso and csu, respectively). Having the

Figure 2. Conceptual model for energy procurement optimization based on Hadera et al.2

• long-term contract (base load or base contract) with fixed price and constant amount delivered at all times • time-of-use (TOU) contract with two price levels (e.g., day and night tariff) • day-ahead (spot) market with hourly varying prices • onsite generation with constant price, with constraints on minimum and maximum runtime of the power plant, a given ramp-up time, and start-up cost In Figure 2, the flow network representation of the energy sources (purchase contracts) and sinks (sales and steel process demand) is shown. The network must always be balanced at the balance node such that the electricity drawn from all sources is matched by the electricity used by the sinks within all time intervals (s) related to energy pricing. All flows between the nodes are represented by arcs that have a certain cost associated as well as capacity bounds that result from the 1338


Article

Industrial & Engineering Chemistry Research knowledge of consumption, the minimum-cost flow network problem can be formulated in order to optimize the electricity purchasing and selling strategy. The flow network model also addresses the constraints that concern the start-up and the shut-down restrictions of the onsite power generation plant as well as the capacity limitations of the purchase and the sales according to different electricity contracts. The flows through the network give rise to the second component of the objective function, the net electricity consumption cost μ. The latter consists of the purchase cost, the earnings from sales, and the start-up cost of the onsite generation plant. Using the resulting consumption variables, the model can respond to the committed load problem and be used to minimize the deviation penalties δ in the objective function. Further details concerning the constraints of the flow network and the load deviation problem can be found in Hadera et al.1 The complete monolithic model notation is shown in Table 1. 3.2. Improvement of the Coupling between the Continuous-Time Scheduling Model and the IntervalBased Energy Cost Calculation. To account for different energy costs in a scheduling problem, a formulation is needed that links the energy optimization to the scheduling problem by computing the electricity consumption in a given fixed time interval. Hadera and Harjunkoski16 and Hadera et al.2 introduced two different energy-aware strategies for continuous-time models. In this work, we improve the latter that is based on event binaries. The binaries capture the time slot in which a task started or finished by the binaries Ysp,s,st and Yfp,s,st (eqs 1−4). t ps , st ≥ Ts − 1 − Ts − 1(1 − Y ps , st , s)

Table 1. Model Notation

P HG HG(P) L(P), F(P) M EAF, AOD, LF, CC S ST SM(ST, M) nodes I and J Pur(node) Dem(node) Gen(node) Bal(node) Sale(node) ARCi,j,s Tp,m t setup m t min m,m′ t max p,st as Ts hp,m

∀ p ∈ P , st ∈ ST , s ∈ S

cs,i,j min max f s,i,j , f s,i,j rmin, dmin cstart k

(1) t ps , st < Ts + (M − Ts)(1 − Y ps , st , s)

∀ p ∈ P , st ∈ ST , s ∈ S (2)

t pf , st > Ts − 1 − Ts − 1(1 − Ypf , st , s)

∀ p ∈ P , st ∈ ST , s ∈ S

c

(3)

t pf , st

≤ Ts + (M − Ts)(1 −

Ypf , st , s)

t sm,p, t fm,p t sp,st, t fp,st wp,st qs Xm,p

∀ p ∈ P , st ∈ ST , s ∈ S

Y ps , st , s , Ypf , st , s ∈ {0, 1} (4)

Using these variables, it is possible to capture different cases of how a task contributes to the resource consumption in a particular time slot of interest (Figure 3). Each production task may require a certain consumption of electricity that is accounted for in this time slot. How much electric power is needed in a time slot can be captured by the set of equations shown in Figure 3. For each case, a detailed description is provided in section 3.2.1. The improvement of the energy-aware formulation is due to two modifications. First, the auxiliary pseudobinaries ysaux p,m,st,s and 1 yfaux p,m,st,s that were used by Hadera et al. are removed. These two variables were designed to take the value of the event binary into account only for those machine−product combinations that were assigned to perform a given task (assignment binary Xp,m = true). The improvement embeds the elimination of the not-assigned tasks in the constraints for capturing the energy consumption as it is explained below in the explanations of the time−task-slot relationship. Second, the way how the model captures the amount of time that a given task spends in a given energy-related time slot are

Vst,p,p′ Ysp,st,s, Yfp,st,s Gs,i,j s gs,i,j faux ysaux p,m,st,s, yp,m,st,s

ap,m,st,s, bp,m,st,s, cp,m,st,s, dp,m,st,s op,m,s bs bos , bus cos , cus fs,i,j cgen s μ δ 1339

Sets heats (products) to be produced heat groups with defined sequence of casting subset of heats p mapped to corresponding heat group hg subset of heats p that are cast, respectively, last or first in casting sequence equipment (machines) subsets of equipment time intervals production stages production stage st mapped to corresponding equipment m nodes in flow network denoting sources and sinks of electricity purchase contracts node production process electricity demand node onsite generation node balancing node electricity sale sink node defined arc between nodes i and j in time slot s Parameters processing duration of heat p on equipment m setup time for machine m minimum transport time from equipment m to m′ maximum hold-up time after stage st preagreed (committed) load curve electricity consumption time slot boundary specific power consumption of processing heat p on equipment m electricity cost of flow from i to j in time slot s minimum and maximum flow between nodes i and j minimum run- and down-time of onsite generation startup cost of onsite generation coefficient of delivered power reduction due to startup of onsite generation coefficient of task start time weight in the objective function Variables starting and finishing time of heat p on equipment m starting and finishing time of heat p at stage st waiting time of heat p after stage st electricity consumed in time slot s binary variable, true when heat p is assigned for processing on equipment m binary variable, true when heat p′ is processed after heat p at stage st binary variable, true when heat p starts or finishes at stage st in slot s binary variable, true when generation is running in time slot s pseudocontinuous positive variable denoting if onsite generation start-up occurred in time slot s auxiliary continuous variable which is true when heat p is assigned for processing on machine m and started or finished processing at stage st in time slot s positive variables accounting for processing time of heat p on equipment m at stage st spent within a slot s positive variable accounting for the processing time of heat p on equipment m within slot s in the improved model buffer for the allowed deviation from committed load in time slot s upper and lower bounds for the buffer in time slot s actual over- and underconsumption in time slot s flow from node i to j in time slot s cost of onsite generation in slot s net electricity consumption cost deviation penalties cost DOI: 10.1021/acs.iecr.5b03239 Ind. Eng. Chem. Res. 2016, 55, 1336−1346

Article


Figure 3. Task−time-slot relations depending on the event binaries.

changed. For a given task that is a combination of a product p and a machine m we introduce a variable op,m,s which will denote how much processing time a given task spends in a particular time slot s. (See Figure 4 for an example.) Note that the

op , m , s ≥ τp , m(Y ps , s , st + Ypf , s , st + X p , m − 2) ∀ p ∈ P , m ∈ M , s ∈ S , st ∈ ST , {st , m} ∈ SM

3.2.1.2. A Task Starts before and Finishes within the Time Slot (Ysp,s,st = 0 and Yfp,s,st = 1). In the second case, the duration is the finishing time minus the time slot boundary, as in eqs 7 and 8. Because we know from the scheduling model that unassigned tasks have starting and finishing times (tsp,m, tfp,m) equal to zero, the inequalities remain valid because of the assignment variable in the expression in eq 5. The M constant should take a value larger than the upper bound of the last time slot, which is equal to the upper bound of the scheduling horizon. op , m , s ≥ t pf , m − τs − 1 − (M − τs − 1)(1 − Ypf , st , s + Y ps , st , s)

Figure 4. Example of the calculation of the time contribution variable.

∀ p ∈ P , m ∈ M , s ∈ S , st ∈ ST , {st , m} ∈ SM

duration of the task is d1 + d2 which corresponds to ∑s op,m,s. This single continuous variable will replace the four continuous variables ap,m,st,s, bp,m,st,s, cp,m,st,s, and dp,m,st,s that were used in Hadera et al.2 3.2.1. Task−Time-Slot Relations. We are considering uniform time intervals of 1 h that are related to energy pricing and to the committed load curve; however, for other interval durations for each time slot (τs − τs−1) the formulation can still be applied, with some exceptions that are explicitly mentioned below. First, a global constraint (valid for all cases below) is applied to enforce an upper bound on the positive variable op,m,s. The time contribution of a given task (p,m) to processing within a given slot cannot be longer than the length of the slot itself or the total processing time of that task. Therefore, the bound has the value of either the length of the time slot (τs − τs−1) or the processing time. 0 ≤ op , m , s ≤ min{τp , m , τs − τs − 1}X p , m

(6)

(7)

op , m , s ≤ t pf , m − τs − 1X p , m + (M + τs − 1)(1 − Ypf , st , s + Y ps , st , s) ∀ p ∈ P , m ∈ M , s ∈ S , st ∈ ST , {st , m} ∈ SM

(8)

3.2.1.3. A Task Starts within and Finishes after the Time Slot (Ysp,s,st = 1 and Yfp,s,st = 0). Similar to the previous case, we derive the constraints eqs 9 and 10 using the expression τsXp,m − tsp,m that captures the exact amount of the contribution of the processing time. op , m , s ≥ τsX p , m − t ps , m − (M + τs)(1 − Y ps , st , s + Ypf , st , s) ∀ p ∈ P , m ∈ M , s ∈ S , st ∈ ST , {st , m} ∈ SM

(9)

op , m , s ≤ τs − t ps , m + (M − τs)(1 − Y ps , st , s + Ypf , st , s) ∀ p ∈ P , m ∈ M , s ∈ S , st ∈ ST , {st , m} ∈ SM

(10)

3.2.1.4. A Task Spans the Whole Time Slot but Starts before Its Beginning and Ends after Its End (Ysp,s,st = 0, Yfp,s,st = 0 and Ysp,s′,st = 1, Yfp,s″,st = 1). Here, we know that the task contribution will be equal to the length of the time slot for assigned tasks (τs − τs−1)Xp,m. To ensure that this contribution is captured for all the slots that are covered by the task, we need to identify the slots in between the time slots where a tasks starts (s′) and finishes (s″) and s′ < s″. Therefore, we can formulate a constraint (eq 11a) that expresses that for all of the slots with the task spanning ∑z s=″−1 s′+1 op,m,z the processing contribution in each slot (assuming (s″ − s′ − 1) slots) is equal to the length of the slot (τs − τs−1)Xp,m; however, this is so only

∀ p ∈ P, m ∈ M, s ∈ S (5)

We identify four different cases of how the task can contribute to the consumption of electricity. 3.2.1.1. A Task Is Processed Entirely within the Time Slot (Ysp,s,st = 1 and Yfp,s,st = 1). The utilization op,m,s within the time slot then is equal to the task processing time τp,m. The consumption must be accounted for only if both event binaries are true, and if the task is assigned to machine m, as in eq 6. Including the assignment Xp,m in eqs 5 and 6 will set the variable op,m,s to zero in the case where a task is not assigned. 1340


Article

Industrial & Engineering Chemistry Research if both respective event binaries are true (Yfp,s″,st + Ysp,s′,st − 1). Because only one among all available parallel machines at each stage can be chosen to process a given product, the nonlinearity in eq 11a can be easily eliminated by summation over machines, leading to the constraint eq 11b. Therefore, the resulting constraint (eq 11) will ensure that the processing contribution is τs″−1 − τs′ in all slots strictly between s′ and s″. Note that replacing the equality constraint by the inequality constraint was possible because τs″−1 − τs′ is the upper bound of the summation of op,m,z over z between s′ + 1 and s″ − 1. In addition, if the length of the slots is uniform, τs − τs−1 = τ, then we know that the constraint should be employed for s″ ≤ s′ + ⌈(maxm∈SMst,m{τp,m})/τ⌉, which defines a window of slots between which the task might occur.

assignment as in eqs 14 and 15. The constraints express that the event binary can hold true only if the task exists (Xp,m = 1).

∑ Y ps ,s ,st = ∑ s∈S

∑ Ypf,s ,st = ∑ s∈S

∀ p ∈ P , st ∈ ST (14)

Xp,m

∀ p ∈ P , st ∈ ST (15)

m ∈ SMst , m

Another set of constraints relates the occurrence of the event binaries with the time contribution op,m,s. For the time horizon where the task has not occurred yet, the sum of start event binaries is zero as well as is the time contribution (eq 16). op , m , s ≤ ∑ Y ps , s ′ , st ∀ p ∈ P , st ∈ ST , s∈ ∑ τ − τ s s−1 m ∈ SM s ′≤ s

s ″− 1

∑

Xp,m

m ∈ SMst , m

op , m , z = (Ypf , s ″ , st + Y ps , s ′ , st − 1)(τs ″− 1 − τs ′)X p , m

st , m

(16)

z = s ′+ 1

Similarly for the finish event binaries, in all the slots after the one for which the binary holds true the time contribution is zero (eq 17). Applying the same thinking, constraints eq 18 and 19 are developed, fixing to zero some of the redundant binary values. op , m , s ≤ ∑ Ypf , s ′ st ∀ p ∈ P , st ∈ ST , s ∈ S ∑ τ − τ s s−1 m ∈ SM s ′≥ s

∀ p ∈ P , m ∈ M , s′, s″ ∈ S , st ∈ ST , {st , m} ∈ SM , s′ < s″ (11a) s ″− 1

∑

∑

op , m , z = (Ypf , s ″ , st + Y ps , s ′ , st − 1)(τs ″− 1 − τs ′)

m ∈ SMst , m z = s ′+ 1

∀ p ∈ P , s′, s″ ∈ S , st ∈ ST , {st , m} ∈ SM , s′ < s″ (11b)

st , m

(17)

(Ypf , s ″ , st + Y ps , s ′ , st − 1)(τs ″− 1 − τs ′)

∑ Y ps ,s ′ ,st ≤ 1− ∑

s ″− 1

≤

∑

∑

op , m , z + M(2 − Y ps , s ′ , st − Ypf , s ″ , st )

s ′> s

m ∈ SMst , m z = s ′+ 1

m ∈ SMst , m

op , m , s τs − τs − 1

∀ p ∈ P , st ∈ ST , s ∈ S

∀ p ∈ P , s′, s″ ∈ S , st ∈ ST , s′ < s″ , s″ , ≤ s′ + ⌈( max {τp , m})/τ ⌉

(18)

m ∈ SMst , m

(11)

∑ Ypf,s ′ ,st ≤ 1− ∑

The above constraint can be considered as a tightening constraint because it has no influence on the proper computation of the variable op,m,s; however, by creating a relatively small number of equations, it helps to decrease the computational time. In any other cases, a task does not contribute to the consumption of electricity in the time slot considered because the task either does not exist or is processed in some other time slot(s). The constraint eq 12 together with eq 11 expresses both situations.

∑ op, m,s = τp,m Xp,m

s ′< s

∑ p∈P ,m∈M

hp , mop , m , s /60

(19)

Another constraint that links the event binaries is useful for handling large scale instances. In eq 20, the fact is that if the event start of a task has not occurred then the event finish cannot have occurred. s

∑ s ′= 1

s

Ypf , s ′ , st ≤ ∑ Y ps , s ′ , st s=1

∀ p ∈ P , st ∈ ST , s ∈ S (20)

We observed that this constraint helps to find good solutions faster when solving larger instances; however for smaller easily tractable instances, it slows down the computational performance. This might be related to the fact that although the constraint reduces the search space somewhat it is creating additional equations that have to be considered by the solver. 3.3. Energy Cost Optimization. The energy purchase and sale optimization that extends the scheduling formulation is described by eqs 12−22 in Hadera et al.1 These equations provide the best purchase and sales strategy for a given electricity consumption curve that results from the scheduling solution. The minimization of the penalties δ paid due to the load commitment problem is considered by eqs 30−32 in Hadera and Harjunkoski.16 The final objective function of the problem is to minimize the net energy purchase/sale costs μ, deviation penalties δ, and the completion times of all tasks tsp,m in a weighted manner (eq 21), similar to that in Hadera et al.:2

(12)

With the above constraints (eqs 1−12) the variable op,m,s will indicate the length of time a given task is processed within a particular time slot s. To account for the consumption of electricity qs in a given time slot s one can simply multiply op,m,s by the specific electricity consumption of the task hp,m. (See eq 13.) The expression is divided by 60 in order to convert MWmin to MWh. qs =

τs − τs − 1

∀ p ∈ P , st ∈ ST , s ∈ S

∀ p ∈ P, m ∈ M

s∈S

m ∈ SMst , m

op , m , s

∀s∈S (13)

3.2.2. Tightening Constraints. This above formulation is sufficient to capture the use of electric energy. From problem knowledge, the model can be tightened by adding the following set of constraints. Following the original formulation in Hadera et al.,1 the event binaries can be related to the machine 1341


Article

Industrial & Engineering Chemistry Research min(μ + δ + c

t ps , m)

∑

considers a committed load curve that was calculated for more products than the production target. This scenario simulates a realistic situation, for example, one where because of equipment breakdown the daily production goal was reduced on short notice. The electricity contracts are as discussed in section 2.2. In Table 3, we report the model statistics and the economic assessment of optimization runs with a CPU-time limit of 600s for the different formulations. As expected, it can be seen that the total numbers of variables and equations are reduced for the improved formulation given above. In all the cases, this improves the computational performance for large-scale monolithic problems. The improvement compared to the original formulation2 is about 1−3% for the two larger cases and 6−9% for the two smaller cases. An economic assessment of the obtained solutions indicates that for most of the problem instances the original and the improved formulation led to similar purchase structures. In Table 4, results for the same problem instances are shown with a CPU-time limit of 3600s. Again, the proposed formulation performed best. Interestingly, for scenario 1 the improvement increased to 6% compared to the original formulation, whereas for scenario 2, the two models performed the same. For the smaller instances, the gap is around 9−18% better than with the six binaries model15,16 and 2−3% better than the original event binaries model.2 The corresponding economic assessment of solutions shown in Table 4 again shows similar purchasing structures, except for scenario 3. For the latter, the original formulation resulted in no use of onsite generation, whereas the improved formulation recognized the benefit of it. To assess more clearly the benefit of the reduced number of binary variables, the performance of the monolithic models is compared for instances where all scheduling binaries have been fixed (sequencing and assignment) in a feasible manner. The complexity of the scheduling part of the problem is then reduced to an LP, and integer variables only result from taking into account the energy-related optimization. In Table 5, the results of optimization runs for the same scenarios are shown. A limit of 2% optimality gap was enforced to avoid long solution times of the models due to long computation times for improving a very small gap. The results show that the proposed formulation outperforms the alternatives for the two larger problem instances. However, for the two smaller ones the six binaries model performs better. This could be related to

(21)

p∈P ,m∈M

4. NUMERICAL EXPERIMENTS The numerical case study considers the same case and input data as described in Hadera et al.2 In the previous study, the Table 2. Investigated Problem Instances scenario

horizon (h)

products

1 2 3 4

24 24 24 24

20 20 16 16

name NM HM RM BH BR

electricity sources and sinks all possible, day-ahead with high prices all possible, day-ahead with low prices all possible, day-ahead with high prices all possible, day-ahead with high prices, overcommitted load curve (as for 20 products) model type

literature monolithic six binaries model (Nolde and Morari,15 Hadera and Harjunkoski)16 original monolithic event binaries model (Hadera et al.2) improved monolithic event binaries model bilevel heuristic decomposition using the original event binaries model (Hadera et al.2) bilevel heuristic decomposition using the improved event binaries model

monolithic formulation was tested with two energy-awareness extensions (one from Nolde and Morari15 and the other the new formulation developed in this paper), and the heuristic bilevel scheme was tested using the event binaries model previously developed by the authors. 4.1. Improvement of the Monolithic Energy-Aware Formulation. The new coupling between the continuous-time scheduling model and the energy optimization was compared to the two models presented in Hadera et al.2 We make the same assumptions that there is a priori knowledge concerning the assignment on the casting machines and of the sequence of the heats that must be processed on each caster. Using the same problem instances and input data as the previous study, the improvement by the new formulation can be assessed. The different scenarios consider varying number of products to be delivered as well as different pricing scenarios taken from historical data that are available publicly. The last scenario

Table 3. Comparison of the Energy-Awareness Formulations in Monolithic Models−600s Computation Limit model statistics

scenario

binary vars

total vars

NM1 HM1 RM1 NM2 HM2 RM2 NM3 HM3 RM3 NM4 HM4 RM4

13017 4065 4065 13017 4065 4065 10181 3229 3229 10181 3229 3229

29508 29508 10308 29508 29508 10308 23428 23428 8068 23428 23428 8068

economic assessment

equations

MIP solution

relative gap

lead times (min)

net cost (€)

electricity purchase (€)

deviation penalties (€)

day-ahead market (MWh)

TOU (MWh)

onsite generation (MWh)

98095 102335 51375 98095 102335 51375 77136 80528 39760 77136 80528 39760

313128 247838 244870 223887 200038 193348 174227 155226 147701 234643 204173 183590

43.78% 29.30% 28.30% 32.30% 24.90% 21.79% 32.63% 22.81% 16.54% 31.53% 22.50% 13.16%

60784 51990 53138 61210 53946 53273 42283 33038 32251 43598 36152 37166

149832 133972 143209 119440 120989 117602 86071 96508 98162 96266 72846 99105

161098 151905 151053 98505 98592 99133 88984 128208 131569 134759 140130 133488

102512 61876 48522 43236 25103 22474 45872 25679 17286 94780 95175 47319

172.55 173.49 223.07 1608.92 1514.51 1580.79 1417.98 82.60 100.55 77.21 142.48 77.73

1471.83 1421.44 1353.58 4.00 95.78 91.95 54.27 1241.16 1270.68 1266.04 1318.80 1257.39

912 952 952 352 432 352 0 952 952 952 952 952

1342


Article

Industrial & Engineering Chemistry Research Table 4. Comparison of the Energy-Awareness Formulations in Monolithic Models−3600s Computation Limit model statistics

scenario

binary vars

total vars


13017 4065 4065 13017 4065 4065 10181 3229 3229 10181 3229 3229

29508 29508 10308 29508 29508 10308 23428 23428 8068 23428 23428 8068

economic assessment

equations

MIP solution

relative gap

lead times (min)

98095 102335 51375 98095 102335 51375 77136 80528 39760 77136 80528 39760

290708 241136 223028 222167 180023 180236 156987 146339 146906 221454 180965 175887

38.97% 26.80% 20.95% 31.20% 16.10% 16.12% 24.63% 17.93% 15.64% 27.20% 12.10% 9.19%

58763 50796 50908 60282 51598 52433 39444 32753 31879 42269 35396 36053

net cost (€)


deviation penalties (€)

day-ahead market (MWh)

TOU (MWh)


146373 142452 142288 116872 115937 115944 83531 95840 108889 100104 94723 101195

162528 14396 151919 98942 98229 100968 89620 130217 133085 127986 128534 120442

85572 47888 29831 45014 12489 11859 34012 17746 6137 79081 50846 38639

203.72 234.22 173.40 1649.35 1635.81 1661.49 1547.37 83.77 119.33 110.87 46.66 19.60

1466.61 1349.66 1425.09 16.83 56.11 74.62 0.60 1270.35 1240.03 1164.93 1256.92 1177.27

952 952 952 312 352 312 0 952 952 952 952 952

problem data that helped to steer the branch-and-bound solver faster toward the 2% gap solution. On the basis of all these tests, it can be concluded that the proposed energy-awareness strategy is more compact and performs better overall when used in the monolithic model than the strategies proposed before. 4.2. Application of the New Formulation within the Heuristic Bilevel Decomposition. Because the monolithic formulation of the problem is unable to cope with problem instances of realistic size, we apply the proposed improved model with event binaries within the heuristic bilevel decomposition that was developed by Hadera et al.2 In the previous study, the original event binaries model was tested on the same problem instances as the monolithic model. We follow the same approach here to see if there is an improvement when applying the new event binaries formulation within the decomposition scheme. Compared to the original heuristic decomposition approach, only the energyawareness part of step is replaced. Thus, the new formulation does not change the decomposition procedure nor the validity

Table 5. Comparison of the Energy-Awareness in Monolithic Models with Fixed Assignment and Sequencing Binaries− Computation Times for 2% Gap model statistics−2% gap scenario

binary vars

total vars

equations

MIP solution

CPUs


12861 3593 3593 12861 3593 3593 10069 2865 2865 10069 2865 2865

29508 29508 10308 29508 29508 10308 23428 23428 8068 23428 23428 8068

98095 102335 51375 98095 102335 51375 77136 80528 39760 77136 80528 39760

195020 194410 194197 166687 166550 166667 129260 135012 135241 174992 174922 174769

64.90 62.47 8.98 7.46 6.48 3.37 7.98 24.63 11.69 39.05 56.23 58.74

Figure 5. Heuristic bilevel decomposition algorithm based on Hadera et al.2 1343


Article

4(3) 5(3)

3(1) 5(1)

5(3) 5(4)

227293 227293

160707 160707

169197 169197

153727 153726

165285 165402

134749 134592

176244 173873

51375 102335 51375 51375 102335 51375 39760 80528 39760 39760 80528 39760

244870 193845 193667 193348 165196 165196 147701 134588 134577 183590 176006 173873

193888 193832

28.30% 9.89% 9.30% 21.79% 9.09% 9.09% 16.54% 9.87% 9.87% 13.16% 8.71% 8.61%

53138 46176 46576 53273 45472 45613 32251 30626 30885 37166 35289 31987

143209 147087 146259 117602 119443 119338 98162 103057 102839 99105 98990 99814

151053 156853 157031 99133 101531 103838 131569 133980 134108 133488 123190 126177

48522 582 832 22474 281 245 17286 904 853 47319 41727 42070

223.065 176.349 179.396 1580.79 1512.868 226.833 100.554 133.033 133.759 77.735 13.13 79.791

1353.578 1456.514 1457.204 91.948 196.161 1521.583 1270.681 1243.566 1243.567 1257.39 1228.59 1228.357

952 952 952 352 392 352 952 952 952 952 952 952

5(5) 5(3)

and feasibility of the solutions as explained by Hadera et al.2 Below, the most important characteristics of the decomposition approach are highlighted; further details can be found in the paper by Hadera et al.2 Similar to the previous study, before running the decomposition algorithm we eliminate some event binary variables by obtaining upper and lower bounds of start times of tasks. To find good quality solutions for the scheduling problem, the heuristic bilevel decomposition considers two separate problems: an upper level approximation UL1 (two distinct steps with UL2 from second iteration of the algorithm on) and a lower-level LL full-scope problem, which is solved with a reduced number of degrees of freedom (Figure 5). The approximation on the upper level considers only two production stages, the most energy-intensive one (EAF) and the most constrained one (CC because of the semicontinuous production mode). This simplification makes the problem tractable. The upper-level problem provides guidance to the lower-level problem on how the jobs in two most important stages should be sequenced and how batches are assigned to machines. The lower-level problem considers all production stages (AOD and LF in addition) with fixed product sequences and assignments to EAFs obtained from UL1. In addition, the range of the event start binary variables on the lower level is restricted as explained in the previous study.2 After each iteration, the full schedule obtained at the lower level problem LL is evaluated within UL2 to find better assignments at the EAF stage for the upper level problem. This is done by fixing all binary variables in UL2 to the values obtained in the solution of the LL problem, except for the assignments of EAFs. As explained in the previous study,2 both full-scope problems LL and UL2 always provide a valid feasible solution (upper bound). The solution of UL2 is always at least as good as the solution of LL. Previously found solutions of UL2 are cut off from the search space of UL1 using cuts as described in Hadera et al.2 Moreover, in the case when within the assigned solution time the gap reaches a certain level, an additional restriction on the event start binary is enforced.2 Because in the UL1 problem the assignments of EAFs are fixed from the second iteration on for the purpose of speeding up the solution times, its solution does not provide an increasing lower bound to the original problem;2 thus, it is a heuristic decomposition strategy. In Table 6, we compare the results obtained with the heuristic bilevel decomposition when running the optimization with a stopping criterion of 600s computation time as in the previous study.2 The values of the relative gap are obtained from the upper and lower bounds of the original problem. The gap is computed by comparing the solution with best bound (LP relaxation reported by the solver) that was obtained from optimization runs of the corresponding monolithic model with the heuristic decomposition solution being provided as the initial solution to the solver. It can be observed that even though the formulation of the energy-awareness has been improved this does not impact the solution quality of the heuristic decomposition much. For scenarios 2−3, both strategies performed similarly; for scenarios 1 and 4, the proposed formulation improved the solution slightly. This is probably due to the fact that the performance of the algorithm depends on the quality of the solution of the upper-level problem. The improved formulation does not have much impact on the quality of the approximation of the relaxed problem. The heuristic decomposition yields very similar solutions, also with respect to the energy procurement

4065 1458 1458 4065 1458 1458 3229 1276 1276 3229 1276 1276

scenario

RM1 BH1 BR1 RM2 BH2 BR2 RM3 BH3 BR3 RM4 BH4 BR4

10308 29508 10308 10308 29508 10308 8068 23428 8068 8068 23428 8068

day-ahead market (MWh) deviation penalties (€) binary vars UL2

total vars UL2

equations UL2

MIP solution UL2

MIP solution LL

MIP solution UL1

relative gap

lead times (min)

net electricity cost (€)


economic assessment model statistics

Table 6. Results of the Heuristic Bilevel Decomposition for Different Model Formulations−600s Computation Limit

TOU (MWh)


no. of iterations (best)


1344



■

structure. With fixed scheduling binaries, the improved model is considerably faster. As a result, the heuristic decomposition is able to perform larger numbers of iterations compared to the other models, as in cases 3 and 4. This increases the chances of finding better solutions compared to the original problem formulation.

REFERENCES

(1) Hadera, H.; Harjunkoski, I.; Grossmann, I. E.; Sand, G.; Engell, S. Steel production scheduling under time-sensitive electricity cost. Comput.-Aided Chem. Eng. 2014, 33, 373−378. (2) Hadera, H.; Harjunkoski, I.; Sand, G.; Grossmann, I. E.; Engell, S. Optimization of steel production scheduling with complex timesensitive electricity cost. Comput. Chem. Eng. 2015, 76, 117−136. (3) Benefits of Demand Response in Electricity Markets and Recommendations for Achieving Them: A Report to the United States Congress Pursuant to Section 1252 of the Energy Policy Act of 2005; U.S. Department of Energy: Washington, DC, 2006. http://energy.gov/oe/ downloads/benefits-demand-response-electricity-markets-and-recommendations-achieving-them-report (accessed February 02, 2014). (4) NERC, Data Collection for Demand-Side Management for Quantifying its Influence on Reliability, 2007, http://www.nerc.com/ docs/pc/drdtf/NERC_DSMTF_Report_040308.pdf, accessed 02.02.2014. (5) dena (German Energy Agency) dena Grid Study II: Integration of Renewable Energy Sources in the German Power Supply System from 2015−2020 with an Outlook to 2025; Deutsche Energie-Agentur GmbH (dena): Berlin, Germany, 2011. (6) Harjunkoski, I.; Maravelias, C.; Bongers, P.; Castro, P.; Engell, S.; Grossmann, I. E.; Hooker, J.; Méndez, C.; Sand, G.; Wassick, J. Scope for industrial applications of production scheduling models and solution methods. Comput. Chem. Eng. 2014, 62 (5), 161−193. (7) Pochet, Y.; Wolsey, L. A. Production Planning by Mixed Integer Programming; Springer: New York, 2006. (8) Méndez, 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 (6−7), 913−946. (9) Harjunkoski, I.; Grossmann, I. E. A decomposition approach for the scheduling of a steel plant production. Comput. Chem. Eng. 2001, 25 (11−12), 1647−1660. (10) Castro, P.; Sun, L.; Harjunkoski, I. Resource−Task Network Formulations for Industrial Demand Side Management of a Steel Plant. Ind. Eng. Chem. Res. 2013, 52 (36), 13046−13058. (11) Ashok, S. Peak-load management in steel plants. Appl. Energy 2006, 83 (5), 413−424. (12) Karwan, M. H.; Keblis, M. F. Operations planning with real time pricing of a primary input. Computers & Operations Research 2007, 34 (3), 848−867. (13) Mitra, S.; Grossmann, I. E.; Pinto, J. M.; Arora, N. Optimal production planning under time-sensitive electricity prices for continuous power-intensive processes. Comput. Chem. Eng. 2012, 38, 171−184. (14) Mitra, S.; Pinto, J. M.; Grossmann, I. E. Optimal Multi-scale Capacity Planning for Power-Intensive Continuous Processes under Time-sensitive Electricity Prices and Demand Uncertainty, Part II: Enhanced Hybrid Bi-level Decomposition. Comput. Chem. Eng. 2014, 65, 102. (15) Nolde, K.; Morari, M. Electrical load tracking scheduling of a steel plant. Comput. Chem. Eng. 2010, 34 (11), 1899−1903. (16) Hadera, H.; Harjunkoski, I. Continuous-time Batch Scheduling Approach for Optimizing Electricity Consumption Cost. Comput.Aided Chem. Eng. 2013, 32, 403−408. (17) Harjunkoski, I.; Sand, G. Flexible and configurable MILP models for meltshop scheduling optimization. Comput.-Aided Chem. Eng. 2008, 25, 677−682. (18) Castro, P.; Grossmann, I. E.; Veldhuizen, P.; Esplin, D. Optimal maintenance scheduling of a gas engine power plant using generalized disjunctive programming. AIChE J. 2014, 60, 2083−2097. (19) Hait, A.; Artigues, C. A hybrid CP/MILP method for scheduling with energy costs. European Journal of Industrial Engineering 2011, 5 (4), 471−489. (20) Castro, P.; Harjunkoski, I.; Grossmann, I. E. New continuoustime scheduling formulation for continuous plants under variable electricity cost. Ind. Eng. Chem. Res. 2009, 48 (14), 6701−6714.

5. DISCUSSION AND FURTHER WORK In this paper we proposed an improvement of the energyawareness formulation1,2 that is used to optimize simultaneously the production schedule (in the sense of minimizing the sum of all completion times) and the cost of energy procurement in industrially relevant problems. The proposed improvement reduces the number of binary variables and constraints. It helps to create a tighter model leading to improvements in the computational performance of a monolithic formulation of the overall problem. However, the monolithic problem formulation is still intractable for realistic problem sizes and scheduling horizons because of the complexity of the scheduling part; therefore following previous work, we also used the new formulation within the solution that employs a heuristic bilevel decomposition. There are three main directions for further improvements. One concerns the inclusion of uncertainties in the process and in the energy prices. Expanding the model according to robust optimization concepts would make the final schedule less vulnerable to changes in the equipment availability or the exact level of predicted electricity prices (especially for the spot market prices and for horizons longer than 1 day). Strategies such as stochastic programming would most likely render the model size intractable; however with proper adaptations of the heuristic decomposition, they might also be considered. In particular, several successful rigorous models and heuristics have been introduced in the recent years so that industrial size problems can be solved (for example, Wu and Ierapetritou,22 Castro et al.,23 and Xu et al.24). A second direction of improvements relates to solving larger time horizons with better quality solutions while preserving the rigorous solution approach. Third, from an implementation point of view, a functional separation of the energy-related problem from the scheduling problem is an interesting line of research suggested by Hadera et al.25 Such functional decomposition makes it possible to couple separate software solutions for scheduling and for energy procurement, leading to reduced development and maintenance cost and having the benefit of potentially exploiting existing partial solutions.

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address

H.H.: BASF SE, Carl-Bosh-Strasse 38, 67056 Ludwigshafen, Germany. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS H. H. acknowledges the Marie Curie FP7-ITN project “Energy savings from smart operation of electrical, process and mechanical equipment− ENERGY-SMARTOPS”, Contract No: PITN-GA-2010-264940 for financial support. 1345


Article

Industrial & Engineering Chemistry Research (21) Castro, P.; Harjunkoski, I.; Grossmann, I. E. Optimal scheduling of continuous plants with energy constraints. Comput. Chem. Eng. 2011, 35 (2), 372−387. (22) Wu, D.; Ierapetritou, M. G. Decomposition approaches for the efficient solution of short-term scheduling problems. Comput. Chem. Eng. 2003, 27 (8), 1261−1276. (23) Castro, P.; Harjunkoski, I.; Grossmann, I. E. Rolling-Horizon Algorithm for Scheduling under Time-Dependent Utility Pricing and Availability. Comput.-Aided Chem. Eng. 2010, 28, 1171−1176. (24) Xu, C.; Sand, G.; Harjunkoski, I.; Engell, S. A new heuristic for plant-wide schedule coordination problems: The intersection coordination heuristic. Comput. Chem. Eng. 2012, 42, 152−167. (25) Hadera, H.; Wide, P.; Harjunkoski, I.; Mäntysaari, J.; Ekström, J.; Sand, G.; Engell, S. A Mean Value Cross Decomposition Strategy for Demand-side Management of a Pulping Process. Comput.-Aided Chem. Eng. 2015, 37, 1931−1936.

1346


An Improved Energy-Awareness Formulation for General Precedence

Recommend Documents