Reactive Scheduling by a Multiparametric Programming Rolling

Article pubs.acs.org/IECR

Reactive Scheduling by a Multiparametric Programming Rolling Horizon Framework: A Case of a Network of Combined Heat and Power Units Georgios M. Kopanos and Efstratios N. Pistikopoulos* Imperial College London, Department of Chemical Engineering, Centre for Process Systems Engineering, SW7 2AZ London, U.K. ABSTRACT: In this work, we introduce an approach for the reactive scheduling of production systems with bounded uncertain parameters. The proposed method follows a state-space representation for the scheduling problem, and relies on the use of a rolling horizon framework and multiparametric programming techniques. We show that by considering as uncertain parameters the set of variables that describe the state of the system at the beginning of the prediction horizon, we can effectively formulate a set of state-space multiparametric programming problems that are solved just once and off-line. In contrast to existing approaches, the repetitive solution of a new multiparametric problem after each disruptive event is avoided. The results of the parametric optimization are used in a rolling horizon basis without the need for online optimization. The proposed multiparametric programming rolling horizon (mp-RH) approach is applied in the scheduling problem of a network of combined heat and power units (i.e., a unit commitment problem type). Several case studies are solved, potential extensions of the proposed method are provided, and challenging areas wherein research is necessary are discussed. events and uncertain input data (e.g., demand fluctuations). Some excellent reviews on scheduling under uncertainty can be found in Aytug et al.,4 Li and Ierapetritou,3 and Verderame et al.5 Additionally, Sahinidis6 provided a neatly written review on theory and methodology of optimization problems under uncertainty. Recently, Subramanian et al.7 discussed that the design of process scheduling algorithms for iterative closed-loop scheduling has received limited attention so far. A major contribution of their work is that they transformed a general Mixed Integer Linear Programming (MILP) process scheduling model into a state-space model, and they showed how common scheduling disruptions can be modeled as disturbances in that state-space model. Also, Zhuge and Ierapetritou8 proposed a closed-loop strategy for the simultaneous scheduling and control of chemical processes subject to disturbances. In this work, we introduce a methodology for the reactive scheduling problem of processes with bounded-form uncertain parameters. In the same line with Subramanian et al.,7 a statespace representation for the scheduling problem is adopted. The proposed reactive scheduling approach basically relies on the use of (i) a state-space representation for the scheduling problem, (ii) parametric optimization tools, so as to achieve low response times as well as to get more insight on the scheduling problem itself, and (iii) a rolling horizon framework suitable for dynamic scheduling problems. We apply our approach in the scheduling problem of a network of Combined Heat and Power (CHP) units (i.e., a unit commitment problem type), and we show how we can formulate a multiparametric programming problem (based on a linear state-space MILP) that is solved just once and off-line. Afterward, we show how the solution of the resulting

1. INTRODUCTION The majority of the scheduling literature deals with the task of schedule generation; however, in practice, schedule generation is only one aspect of the overall scheduling process. Especially in highly dynamic production environments, reactive control is equally important for the successful implementation of scheduling systems. Scheduling approaches should not only be able to generate high-quality schedules but also to react quickly to unexpected changes and to revise schedules in a cost-effective manner.1 Since in most production environments scheduling is a dynamic process, it should be solved reactively with low response times taking into account the current state of the system and updated information so as to capture efficiently the real needs of the scheduling process. For this reason, this work addresses the scheduling problem as a reactive iterative optimization problem, which indeed shares many features with typical control problems. Uncertainty in process operations can originate from various sources. On the basis of the nature of the uncertainty in a process, uncertainty could be classified into (i) model-inherent (e.g., mass/heat transfer coefficients), (ii) process-inherent (e.g., temperature variations, processing time), (iii) external (e.g., product demands, prices), and (iv) discrete (e.g., unit availability).2 Depending on the amount of the available information, the uncertain parameters could be described by (i) probability distributions, if information about the behavior of uncertainty is available, (ii) bounded forms, if there is not enough information and just error bounds can be obtained, and (iii) fuzzy sets, if historical data are not readily available.3 In this work, demand uncertainty is considered. For this type of uncertainty, forecasting techniques based on historical data are usually used to obtain approximate ranges of uncertainty realizations (i.e., bounds), or a probability distributional form for the case in which more information is available. Uncertainty is a major concern in real-world scheduling problems, since scheduling is highly susceptible to unexpected © 2014 American Chemical Society

Received: Revised: Accepted: Published: 4366

July 25, 2013 November 13, 2013 February 17, 2014 February 17, 2014 dx.doi.org/10.1021/ie402393s | Ind. Eng. Chem. Res. 2014, 53, 4366−4386

Industrial & Engineering Chemistry Research

Article

mathematical programming, for multipurpose batch plants that uses a rolling horizon to obtain short response times and deal with the perturbations expected in the future. Fang and Xi10 proposed an event-driven rolling horizon rescheduling strategy for job shop dynamic environments, where jobs arrive continuously, machines breakdown and are repaired, and job due dates may change. Their scheduling algorithm is a hybrid of genetic algorithms and dispatching rules. Sand and Engell11 developed a rolling horizon-type two-stage stochastic programming approach for the scheduling of a polystyrene batch plant considering processing and demand uncertainty. FloresTlacuahuac and Grossmann12,13 studied the simultaneous scheduling and control of multiproduct parallel continuous reactors as a mixed-integer dynamic optimization problem. Recently, Novas and Henning14 presented a constraint programming approach for the rolling horizon-based scheduling of automated wet-etch stations, and emphasized that this scheduling problem must be addressed as an online activity, always considering the current state of the system when new orders arrive. In addition, the rolling horizon concept has found application in a number of studies that cope with the integration of production planning and scheduling under uncertainty. Sand et al.15 developed an online scheduling algorithm, based on a rolling horizon approach, for the solution of a two-level hierarchical planning and scheduling problem where uncertainty in the planning level has been modeled explicitly. van den Heever and Grossmann16 studied the production planning and scheduling of a hydrogen supply network. They proposed a solution strategy in which the planning horizon is divided into planning periods, and for each planning period the scheduling problem is solved in a rolling horizon basis. Wu and Ierapetritou17 proposed a hierarchical solution approach according to which the production for the current stage is provided to the scheduling problem. An iterative process was developed to ensure the consistency between planning and scheduling solutions, and a rolling horizon strategy was adopted to determine the optimal schedule. Sung and Maravelias18 used a computational geometry method to derive the feasible production regions for the scheduling problem, and then incorporated it into a rolling horizon planning model to obtain detailed schedules. Verderame and Floudas19 presented an approach that combines a planning with production disaggregation model with a medium-term scheduling model through a forward rolling horizon. Li and Ierapetritou 20 proposed a rolling horizon method that incorporates valid production capacity information, derived through a parametric programming method, into the planning model so as to improve the quality of the overall solution. Finally, there are some works that used a rolling horizon approach within a MPC framework to address supply chain optimization problems under uncertainty.21−25 2.2. Large-Scale Scheduling. Large-scale scheduling problems result in intractable optimization problems that are hard or even impossible to solve. In an attempt to overcome this burden, decomposition techniques based on rolling horizon (i.e., time-based decomposition) emerged and have been used in many industrial scheduling problems. In general, there is no theoretical guarantee of the optimality of the solutions produced by such methods, however they do provide a bound on the optimal solution of the complete problem. And, if the optimal solution is desired, this bound could be used to initialize the branch-and-bound solution procedure; possibly decreasing its computational effort.

multiparametric programming problem can be used in a rolling horizon fashion to cope with the reactive scheduling problem, and we discuss areas wherein further research is necessary. Of great importance is the fact that the basic methodological concept can be used to tackle large-scale optimization problems. Moreover, the main ideas of our approach could be applied to other types of scheduling problems. The remainder of the paper is organized as follows. In the next section, the use of rolling horizon in scheduling problems under uncertainty and large-scale scheduling problems is described. Then, section 3 provides some theoretical background and an overview on multiparametric programming techniques for scheduling problems. In section 4, the proposed reactive scheduling approach is described in detail, and some remarks on the proposed method are provided. Then, in section 5, the scheduling problem in a network of CHP units is formally introduced and a general state-space mathematical model for this scheduling problem is presented. In section 6, two illustrative problems and two case studies are solved and a discussion on the results is provided. A representative example of the use of the suggested approach in large-scale problems is also included in section 6. Finally, some concluding remarks are drawn and promising research directions are revealed in section 7.

2. THE ROLLING HORIZON CONCEPT The rolling (or receding) horizon principle has been broadly used in control methods, such as in Model Predictive Control (MPC) as a means to deal with the control problem when a cost function has to be minimized in a given time horizon. These approaches are based on the use of (i) receding horizons, in the form of a receding prediction horizon and a receding control horizon, and (ii) a process model that allows output predictions to be obtained. As Figure 1 illustrates, at each step these time

Figure 1. A typical rolling horizon approach: Model Predictive Control (MPC).

horizons are rolled by one discretization time interval, and the control variable values inside the control horizon are determined through an optimization procedure that minimizes a cost function that considers the errors in the prediction horizon between predicted outputs and its reference values. Then, just the control variables values obtained for the next discrete time instant are used, and the procedure is repeated successively. 2.1. Scheduling under Uncertainty. In the literature, there is a number of works that adopted a rolling horizon approach to address reactive scheduling problems. For instance, Rodrigues et al.9 presented a reactive scheduling technique, based on 4367

dx.doi.org/10.1021/ie402393s | Ind. Eng. Chem. Res. 2014, 53, 4366−4386


Article

⎧ z(θ ) := min((c + Hθ )T x + (d + Lθ )T y) ⎪ x ,y ⎪ ⎪ subject to ⎪ ⎪ mpMILP⎨ A(θ )x + E(θ )y ≤ b + Fθ ⎪ n p ⎪ x ∈  , y ∈ {0, 1} ⎪ q min max ⎪ θ ∈ Θ := {θ ∈  |θl ≤ θl ≤ θl , ⎪ l = 1, ..., q}, ⎩

Ovacik and Uzsoy26 and Chand et al.27 presented a number of rolling horizon heuristics for the dynamic parallel machine scheduling with changeover times and for the single machine scheduling problem with release dates, respectively. For largescale batch scheduling problems, Bassett et al.28 proposed two rolling horizon approaches that use a hybrid planning and scheduling model in which just a small part of the overall horizon is determined in detail at each iteration and then, a sequence of such problems is solved in forward or reverse order. In the same line, for medium-term scheduling problems in multipurpose batch plants, Dimitriadis et al.29 proposed forward and backward rolling horizon algorithms that use an aggregate mathematical model to provide a tight estimate of the production capacity over the part of the scheduling horizon not being considered in detail at each iteration. Janak et al.30 presented a mathematical model for the scheduling of a large-scale industrial batch plant, and they decomposed the complete problem by dividing the overall scheduling horizon into smaller subhorizons which were scheduled then in a rolling horizon fashion. Erdirik-Dogan and Grossmann31 developed two planning models for parallel multiproduct batch reactors with changeovers, and they applied a decomposition technique based on a rolling horizon scheme and a relaxation of the detailed planning model. Shaik et al.32 proposed a two-level decomposition scheme integrated into a rolling horizon algorithm which was used to solve an industrial case study. Lima et al.33 developed three rolling horizon algorithms for the long-term scheduling in a multiproduct single-stage continuous manufacturing glass process. Li et al.34 used a rolling horizon method for scheduling an industrial steelmaking continuous-casting process.

where x and y represent the continuous and binary decision variables, while θ denotes the vector of uncertain parameters. Additionally, c ∈ n, d ∈ p, b ∈ m, H ∈ n×q, L ∈ p×q, F ∈ m×q, A(θ) ∈ m×n, and E(θ) ∈ m×p. Moreover, A(θ) and A(θ) are affine mappings with respect to θ. In this work, the problem under consideration (described in section 5) results in a mpMILP problem with right hand side (RHS) uncertainty, as given by ⎧ z(θ ) := min(c Tx + dTy) ⎪ x ,y ⎪ ⎪ subject to ⎪ ⎪ mpMILP‐RHS⎨ Ax + Ey ≤ b + Fθ ⎪ x ∈ n , y ∈ {0, 1} p ⎪ ⎪ q min max ⎪ θ ∈ Θ := {θ ∈  |θl ≤ θl ≤ θl , ⎪ l = 1, ..., q}, ⎩

Although in this study we merely deal with RHS uncertainty, the proposed reactive scheduling methodology described in section 4 can be used to address uncertainty in the objective function coefficients and/or the LHS; assuming that there are available multiparametric programming techniques to cope with the resulting multiparametric problem and that the types of uncertainty are of bounded form. 3.2. Multiparametric Programming Algorithms Overview. Some of the earliest studies on the multiparametric Linear Programming (mpLP) problem include the works of Gal and Nedoma38 and Gal39 who developed algorithms, based on the Simplex algorithm, to solve problems considering uncertainty in the RHS and the objective function coefficients. These works proposed to cover the uncertain parameter space by critical regions, wherein a critical region is defined to be the set of parameters such that a given basis is optimal for the parametric program. There are also significant developments within the research field of process systems engineering. Pertsinidis40,41 presented two algorithms based on a master/slave iteration procedure for the solution of parametric MILP problems for the solution of single uncertain parameter problems and for problems with several uncertain parameters varying in a single direction at the same time. Acevedo and Pistikopoulos42 developed an algorithm to address mpMILP problems involving uncertainty in the RHS constraint vector. The algorithm was based on a branch-andbound method and the solution of mpLP subproblems in the branch-and-bound leaf nodes (by fixing the integer variables to their assigned values). This algorithm was proven to be computationally expensive. Dua and Pistikopoulos43 extended this work and proposed an alternative algorithm to address the general multiparametric case. In their approach, the mpMILP problem is decomposed into a series of mpLP and an MILP

3. MULTIPARAMETRIC PROGRAMMING This section provides some theoretical background and an overview on parametric optimization techniques for scheduling problems. Multiparametric programming has received considerable attention in the open literature, especially due to its important application to MPC.35,36 Nevertheless, the use of multiparametric programming in scheduling problems is rather limited so far. 3.1. Theoretical Background. Parametric programming serves as an analytic tool by mapping the uncertainty in the optimization problem to optimal alternatives, and from this point of view provides the exact mathematical solution of the optimization problem under uncertainty.3 The output of the parametric optimization provides a complete set of profiles of all the optimal inputs to the system as a function of uncertain parameters, and the regions in the space of uncertain parameters where these functions remain optimal. This fact gives further insight in the optimization problem solved, and moreover the offline output could be used to cope with online optimization problems with minimal computational effort; reducing significantly the response time when the uncertain parameters unexpectedly vary. If at the time of decision making all data are known, multiparametric programming is a powerful tool in accounting for the presence of uncertainty in a mathematical model during the optimization stage.37 Generally speaking, scheduling problems are formulated as MILP models. The general multiparametric mixed integer linear programming (mpMILP) problem can be written in the following standard form: 4368



Article

subproblem. The algorithm consists of an initialization step, and an iterative procedure involving the solution of the mpLP subproblem and the MILP subproblem. In addition, Li and Ierapetritou44 presented a branch-and-bound algorithm for the ́ et al.45 presented a new general mpMILP problem, and Faisca algorithm to solve mpMILP problems that involve varying parameters in both the RHS and the objective function. Recently, Wittmann-Holbein and Pistikopoulos46 proposed an approximate solution algorithm for mpMILP problems with uncertainty in the objective function coefficients as well as in the entries of the constraint matrices and vectors. An overview of theoretical and algorithmic advances, and applications in the areas of multiparametric programming and explicit/multiparametric MPC can be found in Pistikopoulos et al.47 3.3. Multiparametric Programming in Scheduling Problems. In the past few years, there have been some tentative attempts to use multiparametric programming in scheduling problems; however, the use of parametric optimization in scheduling problems is still limited. Ryu and Pistikopoulos48 applied parametric optimization in a zero-wait process scheduling problem, and they basically discussed on the practicality of the parametric solution. Ryu and Pistikopoulos49 presented a proactive scheduling approach based on parametric optimization. This approach gives the complete map of all optimal schedules against potential occurrences of uncertainty before the start of the process. Li and Ierapetritou50 used a multiparametric programming framework to study the effect of uncertain product demand, price, and processing time on a scheduling problem. Later, Li and Ierapetritou51 presented a reactive scheduling approach based on the previous work. Recently, Wittmann-Holbein and Pistikopoulos52 developed a proactive strategy for process scheduling problems by considering uncertainty in the RHS and the objective function as well as in the constraint matrix. Existing scheduling approaches based on multiparametric programming are designed to address single disruptive events. For this reason, successive multiple disruptions (occurring in different time instances) are addressed via the repetitive application of the method used (i.e., need to solve a new multiparametric problem after each disruptive event). That means that the parametric optimization key advantage of off-line optimization does not hold when more than a single unexpected event take place. In this work, we show that by considering the current state of the system as an uncertain parameter, we can effectively overcome this apparent drawback of previous methods.

Figure 2. Framework for reactive scheduling via multiparametric programming rolling horizon.

et al.7 The state-space representation allows us to model elegantly the reactive scheduling problem and transform it to a control problem form, similar to typical MPC problems. In simple words, the task of this initial step is to reformulate the original scheduling mathematical model to the corresponding state-space scheduling model. To achieve this, some control theory concepts should be introduced. In control theory, the decision variables that can be manipulated so as to actuate the internal system are called inputs (u), and the internal system variables that belong to the set of variables that fully describes the state of the system in any time instant are referred as states (x). In addition, uncertain parameters of the system are called disturbances (d), and output measurable variables are referred as ouputs (y). In this study, we use a state-space representation for a linear discete-time system that has the following general form:

4. PROPOSED REACTIVE SCHEDULING APPROACH [mp-RH] For reactive scheduling problems that involve bounded uncertain parameters, we introduce a general iterative approach based on multiparametric programming and a rolling horizon framework. Henceforth, the proposed reactive scheduling method will be referred as mp-RH. Figure 2 illustrates the proposed framework for reactive scheduling via multiparametric programming rolling horizon. First, the basic methodological concept of our approach is presented, and consequently some discussion on the proposed approach is provided. 4.1. Basic Methodological Concept. Step 0. State-Space Representation. The initial step of the proposed approach relies on the use of a state-space representation for the scheduling problem, in the same line with the excellent work of Subramanian

x(t + 1) = A(t )x(t ) + B(t )u(t ) + D(t )d(t ) x(t = 0) = x0 y(t ) = C(t )x(t )

where x(t) represents the states vector, u(t) is the inputs vector, d(t) corresponds to the disturbances vector, and y(t) stands for the outputs vector. In addition, for a certain initial time t = 0, the initial states vector of the system x(t = 0) is a given constant vector x0. Figure 3 shows the block diagram for the state-space representation of linear discete-time systems. To derive the state-space scheduling model the complete set of states, inputs, and disturbances should be identified and properly modeled in accordance with the state-space representation. As 4369



Article

the resolution of the new multiparametric programming problem (i.e., during the execution of the new multiparametric programming problem). In order to overcome the aforementioned drawback of existing scheduling approaches based on parametric optimization, we show here how we can derive a truly off-line scheduling method. The key concept is to consider as an additional uncertain parameter (disruption if you wish) the state of the system (i.e., the set of state variables) at the beginning of each prediction horizon x0 (i.e., the initial state of the system). For the continuous state variables appropriate value ranges (i.e., including all possible values) should be defined. At this point, it is emphasized that we will need to formulate and therefore solve (off-line and once) a number of multiparametric programming problems, if there exist binary and/or semicontinuous variables in the set of variables that describe the initial state of the system. This is the fact for the scheduling problem considered in this work. By doing all the above, the proposed method can effectively deal with multiple disruptive events without the need for solving online any new multiparametric programming problems after the occurrence of any disruption. To sum up, in the mp-RH method both disruptions and initial state variables are modeled as uncertain parameters. Step 2. Off-Line Optimization. A salient feature of the proposed approach is that the optimization is performed off-line and just once; assuming that the known parameters will not vary and the uncertain parameters will always fall within the range considered in the multiparametric programming problems. That means, that during the execution of the mp-RH method no online optimization takes place, and the execution time of the proposed method corresponds to the time spent to locate the critical region and realize the functions evaluation. Here, a

Figure 3. Block diagram for the state-space representation of linear discete-time systems.

already demonstrated by Subramanian et al.,7 introduction of new state variables (i.e., lifting variables) and additional constraints may be necessary so as to capture correctly the state of the system. Step 1. Formulating the Multiparametric Programming Problem. Once a state-space scheduling model has been derived, the next step is to reformulate the state-space scheduling model to a state-space multiparametric programming scheduling formulation. To achieve this, disruptions (e.g., demands) are commonly considered as uncertain parameters and their ranges are defined, resulting in a mpMILP model, as the one shown in subsection 3.1. And indeed, existing scheduling approaches using multiparametric programming follow this typical approach, and they suffer from the critical drawback that a new multiparametric programming problem should be solved after the occurrence of any disruptive event. It is common sense that these methods are not reliable when another unexpected event takes place before

Figure 4. Reactive scheduling via a rolling horizon framework. 4370



Article

prediction horizon must be defined for the multiparametric programming problems. Finally, the outputs of the resulting multiparametric programming problems (i.e., critical regions with the corresponding variable functions in terms of the uncertain parameters) are stored as lookup tables in order to be used in the following step. Step 3. Rolling Horizon Framework. The outputs of the previous step can be used in an iterative fashion via a rolling horizon framework, such as the one shown in Figure 4, for addressing reactive scheduling problems. This rolling horizon scheme basically involves the definition of (i) a time horizon (prediction horizon) in which uncertain parameters can be considered to be known with some certainty, and (ii) a time horizon (control horizon) for which the decisions of the optimization (of the prediction horizon) are applied (fixed). The prediction horizon length depends on the problem under consideration while the control horizon istypically equal to one time intervala subset of the prediction horizon. For instance, in Figure 4 the prediction horizon is three time intervals, and the control horizon is one time interval. It must be clear that the initial state of the system in a given prediction horizon h is equal to the final state of the system in the previous control horizon h − 1 (see Figure 4). The system receives feedback (e.g., actual demand, updated state of system) at every discrete time instant. Figure 5 shows a representative algorithm for the rolling horizon by using the output of multiparametric programming. 4.2. Future Horizon Impact. Despite the fact that the solution obtained for each one of the prediction horizons is optimal for the corresponding prediction horizon, there is no theoretical guarantee that the solution of the overall problem will be either optimal or feasible. This is because no information for the future horizon (i.e., outside the prediction horizon) is considered in the optimization of the prediction horizons. While long-term optimality may not always be necessary or meaningful, feasibility is of great importance. Therefore, longer-term information (outside of the optimized prediction horizon) should be taken into account in order to ensure the feasibility of the solution in the total scheduling horizon. In process control, this is typically achieved by adding problem-specific terminal constraints which are chosen so as to ensure that the solution of the finite horizon problem is a feasible and (near-) optimal solution to the infinite horizon problem.53 Once proper terminal constraints are defined, the mp-RH could be embedded into an integrated MPC scheme involving problem-specific forecasting techniques (for the uncertain parameters) for improving the performance of the method. This is one of our current research tasks. In the mp-RH, the values of the terminal constraints should be considered as uncertain parameters, if they may vary along the scheduling horizon. 4.3. Prediction Horizon of Varied Duration. In this work, we consider that the duration of the prediction and control horizons remain fixed throughout the overall scheduling horizon. Nevertheless, during the practical implementation of the mp-RH method varied-duration prediction and control horizons could be used. That way, the duration of each prediction (or even control) horizon will be an additional decision that the scheduler should make. Generally speaking, this decision mainly will depend on the production characteristics of the system, the current state of the system, and the forecasting of the future horizon. For instance, let us consider a nonflexible (not very responsive) production system in which scheduling decisions are desired to

Figure 5. A representative algorithm for rolling horizon via multiparametric programming.

be made in advance of a given number of time intervals. In this low-responsive system, if the uncertain parameters can be estimated within a high degree of certainty for a long horizon, then long prediction and control horizons would be more preferable. In this study, the energy system considered is characterized by frequent energy demand variations. Therefore, short control horizons and larger (if possible) prediction horizons are needed. Obviously, in the case of varied-duration prediction horizons a larger set of multiparametric programming problems should be solved beforehand. Since these computations are performed off-line, this will not affect the responsiveness of the mp-RH method in its online implementation. Finally, if control horizons greater than one time interval are used, the corresponding costs for rescheduling actions within the (already fixed) control horizon will be of particular interest too.54 4.4. Large-Scale Scheduling Problems. At this point, it is emphasized that the mp-RH approach can be considered as a time decomposition method, and thus could be used to address large-scale scheduling problems. More specifically, the total scheduling horizon is decomposed into a number of subperiods, which actually correspond to the prediction horizons of the original concept. Then, a state-space multiparametric programming problem for a subperiod is solved and its output (i.e., critical regions, variable values as functions of the uncertain parameters) 4371



Article

Figure 6. Representative structure of the energy network.

et al.57 presented a linear programming model for the operational planning in a single residential fuel cell system considering the energy market in the UK. Wakui et al.58 developed a mathematical model for the operational planning in an isolated housing complex of 20 households; wherein a residential SO fuel cell system was installed on each building. Recently, Kopanos et al.59 presented a general mathematical programming framework for the operational planning of small-scale CHP-based energy networks for the minimization of total costs (i.e., startup and operating costs, electricity production revenue, sales, and purchases). Finally, some excellent works on the design and planning problem of energy systems based on cogeneration can also be found in the literature.60−62 5.2. Problem Statement. This work addresses the reactive scheduling problem in an energy network consisting of a finite number of CHP units and a heat storage tank, as shown in Figure 6. The energy network can trade electricity with the main electrical grid, or even receive heat from external sources. The resulting scheduling problem under consideration is formally defined in terms of the following items: (i) A given scheduling prediction horizon which is divided into a set of uniform time intervals t ∈ T := {1,...,τ}, defining T̅ := {0,...,τ} time points. Time interval t starts (ends) at time point t − 1 (t). (ii) A set of CHP units i ∈ I having a ramp up limit ρ+i (>0), a ramp down limit ρi−( 1, n = 0

F(̃ i , t , n) = F(i , t )

∀ i ∈ I , t ∈ T , n ∈ Nidn: δidn > 1, n = 0 (6)

Figure 8 displays an illustrative example of how lifting-state ̃ and F̃ (i,t,n) capture the minimum running and variables S(i,t,n) shutdown history for each CHP unit i ∈ I, in accordance with constraints 4 to 6. The CHP unit index is omitted in the example, since just one CHP unit is considered.

For sake of model simplicity, we assume that the fuel consumption cost is of a linear form. Difference Equations for States. In this work, the dynamic evolution of the state of the system is given by the following set of discrete-time difference equations:

Figure 8. Example for lifting-state variables S̃(i,t,n) and F̃ (i,t,n). 4373



Article

Modeling of Startup and Shutdown Constraints. To model minimum runtime and shutdown times, we use the following sets of binary variables: (i) processing binary variables X(i,t) that indicate if the CHP unit is operating during time interval t ∈ T̅ , (ii) startup binary variables S(i,t) that describe if the CHP unit is switched on at the beginning of time interval t ∈ T, and (iii) shutdown binary variables F(i,t) that indicate if the CHP unit is switched off at the beginning of time interval t ∈ T. S(i , t ) − F(i , t ) = X(i , t ) − X(i , t − 1)

X(i , t ) ≥ S(i , t ) +

(7) up

S(̃ i , t , n)

∑

∀ i ∈ I, t ∈ T

proper connection of past decisions to the actual scheduling horizon. To achieve this, the state variables of the current prediction horizon (t ∈ T̅ := {0,...,τ}) and those of the previous control horizon (t* ∈ T̅ * := {0,...,κ}) must be linked. More specifically, the state variable values of the last time point (t* = κ) of the previous control horizon are linked to the first time point (t = 0) of the actual prediction horizon, as given by the following equations:

∀ i ∈ I , t ∈ T : δi > 1

n>0 n ∈ Niup

∑

∀ i ∈ I , t ∈ T : δidn > 1

F(̃ i , t , n)

∀ t ∈ T̅ , t * ∈ T̅ * : t = 0, t * = κ

E(i , t ) = E(i , t *)

∀ i ∈ I , t ∈ T̅ , t * ∈ T̅ * : t = 0, t * = κ

S(̃ i , t , n) = S(̃ i , t *, n − 1) ∀ i ∈ I , t ∈ T̅ , t * ∈ T̅ * , n ∈ Niup: n > 0, t = 0, t * = κ

(8)

1 − X(i , t ) ≥ F(i , t ) +

Bt + 1 = Bt *+1

F(̃ i , t , n) = F(̃ i , t *, n − 1) ∀ i ∈ I , t ∈ T̅ , t * ∈ T̅ * , n ∈ Nidn: n > 0, t = 0, t * = κ

n>0 n ∈ Nidn

(9)

(18)

Constraints 7 define S(i,t) and F(i,t) variables through processing variables X(i,t). Minimum runtime and shutdown times are modeled by constraints 8 and 9, respectively. Notice that state ̃ variables S(i,t,n) and F̃(i,t,n) are included in constraints 8 and 9 so as to properly carryover past switch on and off decisions. Electricity and Heat Constraints. The remaining constraints of the scheduling problem are

A control horizon of one time interval is usually used (i.e., κ = 1); however, longer control horizons could be defined depending on the scheduler. Notice that the RHS terms of constraints 18 are known parameters in the optimization problem of the current prediction horizon. And, importantly they fully describe the state of the system at the beginning of the actual prediction horizon. Domain of Optimization Variables. Finally, the optimization variables domain is defined by

ρi− ≤ R (i , t ) ≤ ρi+

∀ i ∈ I , t ∈ T̅ : t < τ

εiminX(i , t ) ≤ E(i , t ) ≤ εimax X(i , t )

Pt +

∑ E(i ,t) = Wt +

ζte

∀ i ∈ I , t ∈ T̅

Q t ≤ ζthYt

Pt ≤ ζteZt

∀ t ∈ T̅

∀t∈T

Dt ≤ (1 − Yt ) ∑ μi εi

E(i , t ) ≥ 0

(11)

∀ i ∈ I , t ∈ T̅ : t < τ

(12)

Bt + 1 ≥ 0

(13)

Pt , Wt , Q t , Dt ≥ 0

(14) max

∀ i ∈ I , t ∈ T̅

ρi− ≤ R (i , t ) ≤ ρi+

∀t∈T

i∈I

λ min ≤ Bt + 1 ≤ λ max

(10)

∀t∈T

∀ i ∈ I , t ∈ T̅ ∀t∈T

Zt , Yt ∈ {0, 1}

∀t∈T

X(i , t ) ∈ {0, 1}

∀ i ∈ I , t ∈ T̅

i∈I

(15)

S(i , t ) , F(i , t ) ∈ {0, 1}

∀t∈T

(16)

S(̃ i , t , n) ∈ {0, 1}

∀ i ∈ I , t ∈ T̅ , n ∈ Niup: δiup > 1

F(̃ i , t , n) ∈ {0, 1}

∀ i ∈ I , t ∈ T̅ , n ∈ Nidn: δidn > 1

Wt ≤ (1 − Zt ) ∑ εimax

∀t∈T

∀ i ∈ I, t ∈ T

(17)

(19)

Constraints 10 represent the ramp change limits for every CHP unit i ∈ I and time point t ∈ T̅ (t < τ). Electricity production level bounds for every CHP unit i ∈ I are given by constraints 11. Observe that variable X(i,t) = 1 iff εmin ≤ E(i,t) ≤ εmax i i , otherwise X(i,t) = 0. That means, that if E(i,t) is known, X(i,t) is given by 11. Especially, for time point t = 0, since E(i,0) is known X(i,0) can be directly calculated. Constraints 12 ensure that the electricity demand of the system ζet is satisfied in every time interval t ∈ T. Electricity could be acquired from (sold to) the main electrical grid. Constraints 13 give the minimum and maximum heat storage capacity of the system before the beginning of every time interval t = {1,...,τ + 1}. Constraints 14−17 give the upper bounds for all variables that describe energy interchange between the system under study and its external environment. Binary variables Zt = 0 (Yt = 0), if there are no electricity (heat) purchases in time interval t ∈ T. Linking Constraints (Initial Conditions). An important aspect of the proposed iterative reactive scheduling method is the

5.4. Incorporation of Uncertainty. In this particular scheduling problem, for a given prediction horizon, we treat as uncertain parameters of bounded-form the following parameters: the initial inventory level B1 of the heat storage tank, the initial electricity production level E(i,0) of CHP unit i ∈ I, and the demand for heat ζht and electricity ζet in each time interval t ∈ T. Join Initial Heat Storage Level with Heat Demand in First Time Interval. The initial heat storage level and the heat demand in the first time interval of the prediction horizon are two uncertain parameters of the system that could be represented by just one parameter. That way, the number of uncertain parameters considered in the mpMILP problems is reduced by one. In accordance with constraints 3, the heat balance in the heat storage tank before the beginning of the second time interval is given by: B2 = (1 −η)B1 + (∑i ∈ Iμi E(i,1)) + Q1 − D1 − ζh1, where B1 and ζh1 are two of the uncertain parameters considered. If we substitute ζh1̃ = ζh1 −(1 − η)B1, then the following equation is derived:

i∈I

4374



Article

Table 1. Data for the CHP Units i ∈ I i

μi

εmin i

εmax i

ρ+i

ρ−i

δup i

δdn i

αi

βi

ϕi

πi

1 2 3 4 5

2.50 3.50 1.50 2.00 3.00

2.00 1.00 1.25 1.50 0.75

8.00 5.50 1.75 2.50 1.25

2.880 1.650 1.300 1.875 0.785

−3.600 −2.145 −1.500 −2.310 −1.250

2 1 1 1 1

2 1 1 1 1

37.0 66.0 48.0 66.5 45.0

2.226 2.592 2.650 2.800 2.500

34 6 6 5 2

17 3 2 3 1

Table 2. Heat and Electricity Demands for All Problems (MW) example 1

example 2

case study 1

case study 2

t

ζht

ζet

ζht

ζet

ζht

ζet

ζht

ζet

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

5.92 17.34 12.36 8.12 7.97 6.81 8.95 17.56 4.14 11.50 19.97 0.00 0.00 0.00 5.22 13.88 5.71 7.25 14.37 10.40

3.56 3.49 1.62 1.78 5.51 7.56 2.46 6.16 7.09 3.08 1.44 0.00 0.00 0.00 2.75 2.93 5.55 6.64 5.84 4.44

19.29 36.08 28.76 22.53 22.31 20.60 23.75 36.41 16.68 27.51 39.95 7.07 5.59 6.57 18.27 30.99 18.99 21.25 31.72 25.88

5.68 5.57 2.71 2.95 8.66 11.80 4.00 9.65 11.09 4.95 2.44 0.05 0.34 0.18 4.45 9.72 8.72 10.4 9.17 7.03

8.43 11.64 12.32 0.85 6.34 5.14 8.03 14.67 6.20 10.00 19.96 11.57 19.82 15.25 2.61 12.79 3.19 5.00 13.38 8.71

3.89 1.82 3.20 1.33 3.37 2.27 5.52 4.93 6.39 3.08 5.63 6.29 5.39 2.99 1.60 1.72 5.49 4.82 1.22 6.55

7.69 10.37 9.20 8.20 8.17 7.90 8.40 10.43 7.27 9.00 10.99 9.31 10.96 10.05 7.52 9.56 7.64 8.00 9.68 8.74

3.91 3.62 3.81 3.55 3.84 3.68 4.15 4.06 4.27 3.80 4.16 4.26 4.13 3.78 3.59 3.60 4.14 4.05 3.53 4.29

h B2 = (∑ μi E(i ,1)) + Q 1 − D1 − (ζ1̃ − θ1)

Pt +

θ1min ≤ θ1 ≤ θ1max

min max θ(3, t ) ≤ θ(3, t ) ≤ θ(3, t )

(20)

where θ1 represents the deviation of parameter ζh1̃ from its nominal value. Lower and upper bounds on parameter ζ̃h1 are ̃ ̃ = ζh,min − (1 −η)λmax, and ζh,max = ζh,max − (1 calculated by: ζh,min 1 1 1 1 −η)λmin, respectively. Initial Electricity Production Level for Each CHP Unit. The electricity production level for every CHP unit i ∈ I at the beginning of the prediction horizon are also (state related) uncertain parameters. For t = 0, constraints 2 and 11 become

(22)

i∈I

− (ζth + θ(4, t )) min max θ(4, t ) ≤ θ(4, t ) ≤ θ(4, t )

∀ t ∈ T: t > 1 ∀ t ∈ T: t > 1

(23)

where θ(3,t) (θ(4,t)) represent the fluctuation of electricity (heat) demand from its nominal value. 5.5. Remarks on Multiparametric Programming Problems. To capture properly the potential history of minimum running and shutdown times, we additionally considered as ̃ uncertain parameters the binary lifting-state variables: S(i,0,n−1) ∀i dn ̃ ∈ I, t ∈ T̅ , n ∈ Nup : n > 0, and F ∀ ∈ I,t ∈ T , n ∈ N ̅ i (i,0,n−1) i i : n > 0. As a result, a number of mpMILP problems should be formulated and solved (off-line and once) for all possible combinations of these binary decisions made in the previous control horizon. Also, notice that these combinations should take into account all possible cases of the semicontinuous state variables, such as E(i,0) in this study (E(i,0) = 0 or E(i,0) ≥ εmin i ). The resulting mpMILP problems optimize objective 1 subject to constraints 2 and 11 for t > 0, 4−10, and 13−23. Finally, notice that the total number of uncertain parameters (no. of θs) for a given problem is given by the following expression: (no. of θs) = ion + 2PH, where ion is the total number

∀i∈I

∀i∈I

∀t∈T

Bt + 1 = (1 − η)Bt + (∑ μi E(i , t )) + Q t − Dt

εiminX(i ,0) ≤ (E(i ,0) + θ(2, i)) ≤ εimax X(i ,0) ∀ i ∈ I min max θ(2, i) = −E(i ,0) or θ(2, i) ≤ θ(2, i) ≤ θ(2, i)

∀t∈T

i∈I

i∈I

E(i ,1) = (E(i ,0) + θ(2, i)) + R (i ,0)

∑ E(i ,t) = Wt + (ζte + θ(3,t))

(21)

where θ(2,i) corresponds to the variation of the nominal parameter value of E(i,0). If CHP unit i ∈ I is operating (i.e., X(i,0) = 1), the upper and lower bounds on E(i,0) parameter are max min min given by: Emax (i,0) = εi , and E(i,0) = εi , respectively. Obviously, if X(i,0) = 0, E(i,0) = 0 and vice versa. Electricity and Heat Demand Uncertainty. In this study, energy demands are treated as uncertain parameters of bounded form, and are described by the modified version of constraints 12 and 3, as follows: 4375



Article

Table 3. Uncertain Parameters Data for All Problems (where t′ > 1) problem example 1

example 2

case study 1

case study 2

uncertain parameter ζ̃h1 − θ1

nominal value

range

−4.25 5.00 4.50 −3.75 5.00 3.25 7.00 −4.25 5.00 4.50 10.00 2.70 1.50 2.00 1.00 4.00 9.00

−24.25 ≤ θ1 ≤ 24.25 θ2 = −5.00 or −3.00 ≤ θ2 ≤ 3.00 −4.50 ≤ θ3 ≤ 4.50 −43.75 ≤ θ1 ≤ 43.75 θ(2,1) = −5.00 or −3.00 ≤ θ(2,1) ≤ 3.00 θ(2,2) = −3.25 or −2.25 ≤ θ(2,2) ≤ 2.25 −7.00 ≤ θ3 ≤ 7.00 −24.25 ≤ θ1 ≤ 24.25 θ2 = −5.00 or −3.00 ≤ θ2 ≤ 3.00 −3.50 ≤ θ(3,t) ≤ 3.50 −10.00 ≤ θ(4,t′) ≤ 10.00 −8.30 ≤ θ1 ≤ 8.30 θ(2,3) = −1.50 or −0.25 ≤ θ(2,3) ≤ 0.25 θ(2,4) = −2.00 or −0.50 ≤ θ(2,4) ≤ 0.50 θ(2,5) = −1.00 or −0.25 ≤ θ(2,5) ≤ 0.25 −0.50 ≤ θ(3,t) ≤ 0.50 −2.00 ≤ θ(4,t′) ≤ 2.00

E0 + θ2 ζe1 + θ3 ζh1̃ − θ1 E(1,0) + θ(2,1) E(2,0) + θ(2,2) ζe + θ3 ζh1 − θ1 E0 + θ2 ζet + θ(3,t) ζht′ + θ(4,t′) ζh1 − θ1 E(3,0) + θ(2,3) E(4,0) + θ(2,4) E(5,0) + θ(2,5) ζet + θ(3,t) ζht′ + θ(4,t′)

Observe that for the CHP unit 1, parameters δup = δdn = 2, and thus Nup := {0,1} and Ndn := {0,1}. For this reason, lifting-state binary variables S̃(0,0) and F̃ (0,0) are considered as additional uncertain parameters (see subsection 5.5). As a result, four mpMILP problems must be solved for all possible combinations of the initial state binary decisions. The computational results for the four mpMILP problems solved are included in Table 4. Each

of CHP unit that are on at the beginning of the prediction horizon (i.e., E(i,0) > 0), and PH is the number of time intervals considered in the prediction horizon.

6. CASE STUDIES In this section, first we present two academic examples considering a single-period prediction horizon for a single CHP unit and two CHP units. Then, a case study for the singleCHP case under different prediction horizon lengths is presented. Finally, a case study that involves three CHP units and a prediction horizon equal to three time periods is addressed. A total scheduling horizon of 20 time intervals, and a control horizon equal to 1 time interval has been used in all problems. Also, we assume a perfect oracle for energy demands (i.e., the energy demand forecast is accurate for the prediction horizon considered). In the problems under consideration, observe that there exist ̃ and F̃(i,0,n)) that are part of the set of binary variables (i.e., S(i,0,n) variables that describe the initial state of the system. For this reason, a number of mpMILP problems must be solved for all possible combinations of the initial state binary decisions. The mpMILP problems have been solved in the Parametric Optimization (POP) MATLAB toolbox, with an interface to CPLEX 9, running on a Linux workstation (Dual 4 Core Intel Xeon processor, 1.6 GHz, 4 GB RAM) with MATLAB being single threaded. The data for the CHP units can be found in Table 1. Notice that the execution time of the mp-RH corresponds to the time spent to locate the critical region and realize the functions evaluations; no optimization solver is called during the execution of the proposed method. 6.1. Illustrative Examples. For simplicity and presentation purposes a single-period prediction horizon is considered. Example 1: A Single CHP Unit and Single-Period Prediction Horizon. In this problem, we consider one CHP unit (CHP unit 1). Table 2 contains the heat and electricity demands per time interval for this example problem. Notice that there is no energy demand from time interval 12 to 14. The nominal values as well as the range of the continuous uncertain parameters are shown in Table3; index i is omitted. At the beginning of the scheduling horizon, E0 = 3.20, B1 = 3.00, and S̃(0,0) = 1. For the heat storage tank, λmin = 0, λmax = 30, and η = 0.05.

Table 4. Example 1: Computational Results for the mpMILP Problems cases

E0

S̃(0,0)

F̃ (0,0)

no. of θs

CRs

CPU s

figure

A B C D

≥ε ≥εmin =0 =0

0 1 0 0

0 0 0 1

3 3 2 2

22 14 10 2

3.41 2.45 0.91 0.79

Figure 9a Figure 9b Figure 9c Figure 9d

min

mpMILP problem involves 8 continuous and 12 binary variables. Notice that in cases C and D the CHP unit is off (i.e., X0 = 0) at the beginning of the prediction horizon, which means that θ2 = −5 (i.e., E0 = 0). A comparison of the results of the four cases shows that the number of critical regions (CRs) is reduced as the degrees of freedom decrease. Figure 9 illustrates the 48 critical regions found. Notice that for every critical region, the values of the binary variables as well as those of the continuous variables as a function of the uncertain parameters are stored. Table 5 shows an illustrative example of the stored information for the two critical regions of case D, shown in Figure 9d. Afterward, the reactive scheduling problem for the overall time horizon has been solved using the mp-RH method (objective value equals to 1849.8 m.u.). The total execution time of the method was 0.005 s; i.e., 0.0003 s per iteration. The same problem has been solved in 1.11 s by using online optimization in the same rolling horizon scheme. Table 6 shows the detailed schedule for this illustrative example; including also the binary variable values. The heat storage tank level is represented by variables Bt+1 (i.e., in time interval 1 is shown B2). Observe that the CHP unit is not operating in time intervals 14 and 15. Significant electricity purchases are reported for time interval 15 wherein the CHP unit is off but energy demand does exist; that is an apparent result of the short prediction horizon used. No heat is acquired from external sources during the scheduling horizon. 4376



Article

Table 5. Example 1: Information Stored for the Critical Regions of Case D critical region (left cyan-box)

critical region (right red-box)

−24.25 ≤ θ1 ≤ −4.25

−4.25 ≤ θ1 ≤ 24.25

−4.5 ≤ θ2 ≤ 4.5

−4.5 ≤ θ2 ≤ 4.5

continuous variables

binary variables

continuous variables

binary variables

E0 = 0 E1 = 0 B2 = 0 Q1 = −θ1 − 4.25 D1 = 0 P1 = θ2 + 4.5 W1 = 0 R0 = 0

X0 = 0 X1 = 0 S1 = 0 F1 = 0 S̃(0,0) = 0 S̃(1,0) = 0 S̃(1,1) = 0 F̃ (0,0) = 1 F̃ (1,0) = 0 F̃ (1,1) = 1 Y1 = 1 Z1 = 1

E0 = 0 E1 = 0 B2 = θ1 + 4.25 Q1 = 0 D1 = 0 P1 = θ2 + 4.5 W1 = 0 R0 = 0

X0 = 0 X1 = 0 S1 = 0 F1 = 0 S̃(0,0) = 0 S̃(1,0) = 0 S̃(1,1) = 0 F̃(0,0) = 1 F̃(1,0) = 0 F̃(1,1) = 1 Y1 = 0 Z1 = 1

Significant heat disposal is reported in time intervals that are characterized by high electricity demand and low heat demand, such as time intervals 6 and 9. An increased cost is observed in time intervals wherein electricity purchases and heat disposal take place (see time intervals 6, 9, and 15). Example 2: Two CHP Units and Single-Period Prediction Horizon. In this problem, we consider two CHP units (CHP unit 1 and 2). Table 2 gives the heat and electricity demands per time interval for this example. The nominal values and the range of the continuous uncertain parameters are included in Table 3; index t is omitted. At the beginning of the scheduling horizon, E(1,0) = 4.00, E(2,0) = 2.70, B1 = 5.00, S̃(1,0,0) = 0 and F̃(1,0,0) = 0. For the heat storage tank, λmin = 0, λmax = 50, and η = 0.05. dn Once again, notice that for CHP unit 1, parameters δup 1 = δ1 = dn N := {0,1} and := {0,1}. For this reason, lifting2, and thus Nup 1 1 state binary variables S̃(1,0,0) and F̃ (1,0,0) are considered as additional uncertain parameters (see subsection 5.5). In this case study, eight mpMILP problems should be solved (off-line and once) for all possible combinations of the initial state binary variables. The computational results for the mpMILP problems solved are presented in Table 7. Each mpMILP problem involves 11 continuous and 16 binary variables. A total number of 802 critical regions has been found. Figures 10 and 11 illustrate the critical regions for mpMILP problem cases C to H. The reactive scheduling problem for the total time horizon has been solved using the mp-RH method (objective value equals to 2536.8 m.u.). The total execution time was 0.0496 s; i.e., 0.0025 s per iteration. The same problem has been solved in 1.34 s by using online optimization in the same rolling horizon scheme. The detailed energy schedule for this problem can be found in Table 8. Electricity purchases take place only in time interval 16 (0.07 MW) while no heat is acquired from external sources throughout the total scheduling horizon. CHP unit 2 is not operating in time intervals 14 and 15. Higher cost is observed in time periods wherein heat disposal takes place (e.g., time intervals 6, 9, and 18). Electricity sales are reported in most time intervals, due to the general high heat demand compared to moderate electricity needs. Figure 12 presents the execution time (per iteration) of the proposed mp-RH method for both illustrative examples. Observe that the execution time in both examples is lower in iterations 15 and 16, because the solution is sought in multiparametric profile

Figure 9. Example 1: Critical regions for all possible cases.

Also, notice that important electricity sales are observed in time intervals that are characterized by high heat demand and relatively low electricity demand, such as time intervals 3 and 11. 4377



Article

Table 6. Example 1: Detailed Energy Schedule for All Time Intervals E0 E1 B2 Q1 D1 P1 W1 R0 X0 X1 S1 F1 ̃ S(0,0) ̃ S(1,0) ̃ S(1,1) F̃(0,0) F̃(1,0) F̃(1,1) Y1 Z1 E0 E1 B2 Q1 D1 P1 W1 R0 X0 X1 S1 F1 S̃(0,0) S̃(1,0) S̃(1,1) F̃ (0,0) F̃ (1,0) F̃ (1,1) Y1 Z1

1

2

3

4

5

6

7

8

9

3.20 6.08 12.13 0.00 0.00 0.00 2.52 2.88 1 1 0 0 1 0 1 0 0 0 0 0 11

6.08 8.00 14.18 0.00 0.00 0.00 4.51 1.92 1 1 0 0 0 0 0 0 0 0 0 0 12

8.00 8.00 21.11 0.00 0.00 0.00 6.38 0.00 1 1 0 0 0 0 0 0 0 0 0 0 13

8.00 7.22 30.00 0.00 0.00 0.00 5.44 −0.78 1 1 0 0 0 0 0 0 0 0 0 0 14

7.22 5.51 30.00 0.00 4.31 0.00 0.00 −1.71 1 1 0 0 0 0 0 0 0 0 0 0 15

5.51 7.56 30.00 0.00 10.59 0.00 0.00 2.05 1 1 0 0 0 0 0 0 0 0 0 0 16

7.56 4.18 30.00 0.00 0.00 0.00 1.72 −3.38 1 1 0 0 0 0 0 0 0 0 0 0 17

4.18 7.06 28.59 0.00 0.00 0.00 0.90 2.88 1 1 0 0 0 0 0 0 0 0 0 0 18

7.06 7.09 30.00 0.00 10.75 0.00 0.00 0.03 1 1 0 0 0 0 0 0 0 0 0 0 19

10 7.09 5.20 30.00 0.00 0.00 0.00 2.12 −1.89 1 1 0 0 0 0 0 0 0 0 0 0 20

5.20 8.00 28.53 0.00 0.00 0.00 6.56 2.80 1 1 0 0 0 0 0 0 0 0 0 0

8.00 4.40 30.00 0.00 8.10 0.00 4.40 −3.60 1 1 0 0 0 0 0 0 0 0 0 0

4.40 2.00 30.00 0.00 3.50 0.00 2.00 −2.40 1 1 0 0 0 0 0 0 0 0 0 0

2.00 0.00 28.50 0.00 0.00 0.00 0.00 −2.00 1 0 0 1 0 0 0 0 1 0 0 0

0.00 0.00 21.86 0.00 0.00 2.75 0.00 0.00 0 0 0 0 0 0 0 1 0 1 0 1

0.00 2.88 14.08 0.00 0.00 0.05 0.00 2.88 0 1 1 0 0 1 0 0 0 0 0 1

2.88 5.76 22.07 0.00 0.00 0.00 0.21 2.88 1 1 0 0 1 0 1 0 0 0 0 0

5.76 6.64 30.00 0.00 0.31 0.00 0.00 0.88 1 1 0 0 0 0 0 0 0 0 0 0

6.64 6.35 30.00 0.00 0.00 0.00 0.51 −0.29 1 1 0 0 0 0 0 0 0 0 0 0

6.35 4.76 30.00 0.00 0.00 0.00 0.32 −1.59 1 1 0 0 0 0 0 0 0 0 0 0

Table 7. Example 2: Computational Results for the mpMILP Problems cases

E(1,0)

E(2,0)

S̃(1,0,0)

F̃(0,0)

no. of θs

CRs

CPU s

figure

A B C D E F G H

≥εmin 1 ≥εmin 1 ≥εmin 1 ≥εmin 1

≥εmin 2 ≥εmin 2

0 1 0 1 0 0 0 0

0 0 0 0 0 1 0 1

4 4 3 3 3 3 2 2

285 195 65 61 121 35 28 12

56.33 35.40 10.99 10.14 21.90 4.69 4.82 2.49

Figure 10a Figure 10b Figure 10c Figure 10d Figure 11a Figure 11b

=0 =0 =0 =0

0 0 ≥εmin 2 ≥εmin 2 0 0

Figure 9c). In the second example (i.e., EX2), note that in iterations 15 and 16 the solution is sought in case C (see Table 7 or Figure 10a). 6.2. Case Study 1. In this part, we study a single-CHP case study (CHP unit 1) under different prediction horizon lengths. The main scope of this case study is to show (i) how the

cases with a reduced total number of critical regions. More specifically, in iteration 15 of the first example (i.e., EX1), the solution is sought in case D that includes just two critical regions (see Table 4 or Figure 9d), because E0 = 0 and F̃(0,0) = 1. For the same example, in iteration 16 the solution is sought in case C (E0 = 0 and F̃ (0,0) = 0) that contains 10 critical regions (see Table 4 or 4378



Article

Figure 11. Example 2: Critical regions for cases G and H.

regions were obtained. For a prediction horizon equal to two and three time intervals, 245 and 2072 critical regions have been found, respectively. The heat and electricity demands per time interval for this case study can be found in Table 2. At the beginning of the scheduling ̃ horizon, E0 = 3.20, B1 = 3.00, and S(0,0) = 1. For the overall time horizon, the reactive scheduling problem has been solved using the mp-RH method under the three different prediction horizon lengths considered. Figure 13 shows the aggregated objective function value for PH = {1,2,3}. As expected, the objective function value is improved by increasing the length of the prediction horizon. More specifically, a PH = 1 results in a total objective function value of 1700.7 m.u., a PH = 2 gives a total objective value equal to 1073.8 m.u, and a total objective value of 896.9 m.u. is obtained for PH = 3. Having as a reference the total objective value for PH = 1, the total objective value for PH = 2 and PH = 3 has been reduced by 37% and 47%, respectively. Additionally, the total objective value for PH = 3 is 16% lower than that for PH = 2. It is noted here, that if the energy demands are known with certainty for all time intervals at the beginning of the time horizon of interest (i.e., so-called perfect information case), one could solve the scheduling problem for the whole scheduling horizon at once. In this example, the objective function value for the perfect information case is equal to 831.1 m.u., which is around 7% better than that for PH = 3. With this simple example, it has been demonstrated that longer prediction horizons favor a better total objective value (under the assumptions made in this study). To continue, Figure 14 illustrates the detailed schedules for this case study for PH = {1,2,3}. Heat and electricity purchases are zero in all cases. Note that the heat storage tank level is represented by variables Bt+1 (i.e., in time interval 1 is shown B2). Observe that longer prediction horizons result in reduced heat storage levels and heat disposals. Also, notice that shorter prediction horizons result in higher total electricity production

Figure 10. Example 2: Critical regions for cases C to F.

complexity of the resulting mpMILP problems grows by increasing the prediction horizon length, and (ii) the apparent advantages of considering longer prediction horizons (under the assumption that such information is available). The nominal values and the range of the continuous uncertain parameters are given in Table 3; note that index i is omitted. For the heat storage tank, λmin = 0, λmax = 30, and η = 0.05. Table 9 presents the computational results for the mpMILP problems for the different prediction horizon lengths considered (i.e., PH = {1,2,3}). As expected, the total number of critical regions (and the corresponding computational time) increases by considering longer prediction horizons. Considering the same initial binary state cases, an increment of one time interval in the prediction horizon, on average increases the number of critical regions by a factor of 8. For a single-period prediction horizon, 31 critical 4379



Article

Table 8. Example 2: Detailed Energy Schedule for All Time Intervals t

E(1,t)

E(2,t)

Bt+1

Qt

Dt

Pt

Wt

R(1,t−1)

R(2,t−1)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

6.88 8.00 8.00 8.00 7.17 8.00 8.00 8.00 8.00 8.00 8.00 4.40 2.00 3.63 6.51 8.00 8.00 8.00 8.00 8.00

4.35 5.50 5.50 5.50 3.36 3.80 1.79 3.44 3.09 2.86 4.51 2.37 1.00 0.00 0.00 1.65 3.30 2.40 4.05 2.41

17.89 20.16 29.64 44.88 50.00 50.00 50.00 43.12 50.00 50.00 43.34 50.00 50.00 50.00 45.50 38.01 48.67 50.00 49.96 50.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00 10.20 0.00 0.00 5.09 0.00 0.00 3.38 0.41 0.00 0.00 0.00 0.00 3.39 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.07 0.00 0.00 0.00 0.00

5.55 7.93 10.79 10.55 1.87 0.00 5.79 1.79 0.00 5.91 10.07 6.72 2.66 3.45 2.06 0.00 2.58 0.00 2.88 3.38

2.88 1.12 0.00 0.00 −0.83 0.83 0.00 0.00 0.00 0.00 0.00 −3.60 −2.40 1.63 2.88 1.49 0.00 0.00 0.00 0.00

1.65 1.15 0.00 0.00 −2.15 0.45 −2.01 1.65 −0.35 −0.23 1.65 −2.15 −1.37 −1.00 0.00 1.65 1.65 −0.90 1.65 −1.64

Figure 12. Execution time of the mp-RH method per iteration for the illustrative examples.

(PH = 1 → 116.6 MW, PH = 2 → 101.9 MW, and PH = 3 → 95.5 MW) and total electricity sales (PH = 1 → 39.1 MW, PH = 2 → 24.4 MW, and PH = 3 → 18.0 MW). That means, that shorter prediction horizons generate solutions that are characterized by (unnecessary) excessive energy production. As a result, the excessive heat generated should be stored or disposed, thus resulting in increased costs. Figure 15 displays the execution time (per iteration) of the proposed mp-RH method for the different prediction horizon lengths considered (PH = {1,2,3}). The total execution time of the mp-RH method was 0.005 s for PH = 1, 0.015 for PH = 2, and 0.102 s for PH = 3. Not surprisingly, the execution time is higher

for longer prediction horizons due to the increased number of total critical regions (i.e., more search options to locate the critical region of interest). Therefore, there is a clear trade-off between the response time and the quality of the solution obtained. 6.3. Case Study 2. In this case study, we consider a network of three CHP units (CHP unit 3, 4, and 5) and a prediction horizon equal to three time intervals. For the heat storage tank, λmin = 0, λmax = 14, and η = 0.10. Table 3 contains the nominal values and the range of the continuous uncertain parameters, and Table 10 presents the computational results for the resulting eight mpMILP problems. Each mpMILP problem consists of 36 4380



Article

function value for the perfect information case is equal to 2822.1 m.u. (3% better than the one reported for PH = 3), and it was solved in 0.705 s. 6.4. Remarks on Large-Scale Optimization. It should be clear to the reader that due to our assumption of a perfect oracle for energy demands (i.e., the energy demand forecast is accurate for the prediction horizon considered), all the problems considered could also be viewed as deterministic problems that could be solved through forward time decomposition by the proposed mp-RH method. In Case Study 1, if the complete problem for the 20 time intervals is decomposed in 3-period subproblems (i.e., PH = 3 case) a solution of 896.9 m.u. is obtained in 0.102 s by the suggested mp-RH method, while the typical online optimization rolling horizon method reported the same solution in 1.386 s. The deterministic complete problem (i.e., equivalent to the perfect information case considered previously) was solved in 0.135 s reporting an optimal solution of 831.1 m.u., which is 7% better that that obtained through the time-decomposition methods. In Case Study 2, if the complete problem for the 20 time intervals is decomposed in 3-period subproblems a solution of 2910.3 m.u. is found in 0.3959 s by the mp-RH method, while the typical online optimization rolling horizon method gave the same solution in 1.302 s. The deterministic complete problem (i.e., equivalent to the perfect information case considered previously) was solved in 0.705 s reporting an optimal objective value of 2822.1 m.u., which is 3% better that that reported by the time-decomposition methods. Overall, the proposed mp-RH method is considerably faster than the typical online time-decomposition approach. As expected, the mp-RH method cannot ensure the optimality of the complete optimization problem; however, it is faster than solving the complete problem at once. As the length of the overall scheduling horizon increases (i) the computational performance of the mp-RH is expected to increase in comparison with the

Table 9. Case Study 1: Computational Results for the mpMILP Problems

PH

no. of cont vars

no. of bin vars

3

22

30

2

15

21

1

8

12

cases

E0

S̃(0,0)

A B C D A B C D A B C D

≥εmin ≥εmin =0 =0 ≥εmin ≥εmin =0 =0 ≥εmin ≥εmin =0 =0

0 1 0 0 0 1 0 0 0 1 0 0

F̃(0,0)

no. of θs

CRs

CPU s

0 0 0 1 0 0 0 1 0 0 0 1

7 7 6 6 5 5 4 4 3 3 2 2

605 652 715 100 84 84 63 14 11 10 8 2

243.82 338.72 217.96 34.75 10.38 10.70 6.55 2.26 1.69 1.06 0.98 0.78

PH: length of prediction horizon.

continuous and 36 binary variables. A total number of 16703 critical regions have been found. The heat and electricity demands per time interval for this case study are given in Table 2. At the beginning of the scheduling horizon, E(3,0) = 1.65, E(4,0) = 0, E(5,0) = 1.05, and B1 = 2.35. For the overall time horizon, the reactive scheduling problem has been solved using the proposed mp-RH method (objective value equal to 2,910.3 m.u.). The total execution time of the mp-RH method was 0.3959 s; i.e., 0.0220 s per iteration. The same problem has been solved in 1.302 s by using online optimization in the same rolling horizon scheme. Figure 16 displays the detailed schedule obtained. Note that the heat storage tank level is represented by variables Bt+1 (i.e., in time interval 1 is shown B2). Heat and electricity purchases are zero in all time intervals. No heat disposals are observed. In this problem, the objective

Figure 13. Case Study 1: Aggregated objective value profiles for PH = {1,2,3}. 4381



Article

Figure 14. Case study 1: Detailed schedules for different prediction horizon lengths (PH = {1,2,3}).

Figure 15. Case Study 1: Execution time of mp-RH method per iteration for PH={1,2,3}.

4382



Article

Table 11 presents (i) the best solution found within the predefined time limit for the complete optimization (denoted as CO), (ii) the solution obtained from the proposed mp-RH heuristic, and (iii) the solution of the complete optimization at the same time point that mp-RH terminates (denoted as CO*). No complete optimization problem has been solved to optimality within the time limit imposed (relative gaps: |T| = 100 → 0.8%, |T| = 200 → 1.2%, and |T| = 300 → 1.7%). Observe that the mp-RH heuristic was able to generate good solutions with a relative gap from the CO solution of 1.5% in minor computational times for all problems considered. In addition, as expected, mp-RH is always faster than the CO case. In |T| = 200 problem, CO* performed better than mp-RH. However, in the larger problem instance (|T| = 300), the mp-RH is considerably better than both the CO and the CO*. At this point, it should be emphasized that the salient feature of the mp-RH method is that during its implementation online optimization is avoided, a fact that ensures a more stable computational performance than any other approach that relies on online optimization. In other words, the mp-RH can be considered as a heuristic method that transforms the original optimization problem into a critical region searching problem. Obviously, the total number of critical regions affect the computational performance of the proposed method during its implementation, since the larger is the total number of critical regions, the higher is the search time to locate the critical region of interest.

Table 10. Case Study 2: Computational Results for the mpMILP Problems cases

E(3,0)

E(4,0)

E(5,0)

no. of θs

CRs

CPU s

A B C D E F G H

≥εmin 3 ≥εmin 3 ≥εmin 3 =0 ≥εmin 3 =0 =0 =0

≥εmin 4 ≥εmin 4 =0 ≥εmin 4 =0 ≥εmin 4 =0 =0

≥εmin 5 =0 ≥εmin 5 ≥εmin 5 =0 =0 ≥εmin 5 =0

9 8 8 8 7 7 7 6

2933 2020 1671 2367 2808 2520 1482 902

32666 9596 6388 9887 6122 7954 4004 1729

complete optimization problem and (ii) the gap from the optimal solution will probably grow. 6.4.1. An Example: Large-Scale Optimization via the mpRH. In this part, we consider as a reference the deterministic case of Case Study 2 with considerably longer total scheduling horizons. More specifically, three problems are addressed considering a scheduling horizon of 100, 200, and 300 time intervals, respectively. These deterministic large-scale problems have been solved using the proposed mp-RH method (forward rolling horizon with subproblems of three time intervals each) as well as the typical complete optimization approach (full optimization). A time limit of 300 and 1800 s has been set for the first and the last two complete optimization problems, respectively.

Figure 16. Case Study 2: Detailed schedule per time interval (PH = 3). 4383



Article

Table 11. Large-Scale Problems: Complete Optimization vs Proposed mp-RH method |T| = 100

|T| = 200

|T| = 300

method

objective value

CPU s

gap from CO

objective value

CPU s

gap from CO

objective value

CPU s

gap from CO

CO mp-RH CO*

14 302 14 510 14 524

300 2 2

1.5% 1.6%

28 425 28 863 28 562

1800 4 4

1.5% 0.5%

43 250 43 896 46 020

1800 6 6

1.5% 6.4%

■

7. CONCLUSIONS We introduce the reactive scheduling mp-RH method that relies on the use of a state-space model representation for the scheduling problem, a rolling horizon framework, and multiparametric programming techniques. To avoid the repetitive solution of new multiparametric problems after the occurrence of each disruptive event, we consider as uncertain parameter the complete set of variables that fully describe the initial state of the system. The mp-RH has been successfully applied in the scheduling problem of a network of CHP units. The proposed methodology could be used in other types of scheduling problems. To do this, first the scheduling problem of interest should be reformulated in a state-space form, and then the initial state variables of the system must be identified and considered as uncertain parameters of bounded form (the ones that are of continuous nature). Of great significance is also the fact that the mp-RH can be considered as a time-decomposition approach for the solution of large-scale optimization problems. In this case, different rolling horizon schemes (forward or backward) could be used. The larger is the prediction horizon considered, the more likely is it to obtain better solutions (or bounds) for the complete problem. The main limitation of the proposed method to date is the increased computational burden for solving big multiparametric programming problems (e.g., problems with many state variables and/or long prediction horizons). In this study, we solved problems with nine uncertain parameters. Notice that the proposed state-space mathematical formulation could also be used in an online rolling horizon reactive scheduling approach. Further research is necessary toward various directions. For instance, the study and the definition of proper terminal constraints and the development of forecasting techniques to take into account information of future horizons could strongly enhance the performance of the proposed method in practice. The study of different prediction horizon durations would be of particular interest too. Moreover, the mp-RH could be used in tandem with an online rolling horizon approach resulting in a new hybrid reactive scheduling approach. Finally, the study of coordinated decentralized MPC methods for solving problems involving a bigger set of units would be of great interest too.

■

NOMENCLATURE

Indices

i = CHP unit t = time interval/point n = time interval (for lifting-state variables) Sets

I = set of CHP units T = set of t = {1,...,τ} T̅ = set of t = {0,1,...,τ} dn Ndn i = set of n = {0,...,δi − 1} for CHP unit i ∈ I up Nup = {0,...,δ = set of n i i − 1} for CHP unit i ∈ I Parameters

αi = fuel consumption cost coeficcient for CHP unit i ∈ I βi = fuel consumption cost coeficcient for CHP unit i ∈ I γt = inventory cost δdn i = minimum shutdown time for CHP unit i ∈ I δup i = minimum running time for CHP unit i ∈ I εmax = maximum electricity production capacity of CHP unit i i ∈I = minimum electricity production capacity of CHP unit i εmin i ∈I ζet = total electricity demand in time interval t ∈ T (nominal value) ζht = total heat demand in time interval t ∈ T (nominal value) η = heat loss rate for the heat storage tank λmax = maximum heat storage tank capacity λmin = minimum heat storage tank capacity μi = heat to electricity production ratio for CHP unit i ∈ I νt = electricity selling price to external source in time interval t ∈T ξt = heat purchase price from external source in time interval t ∈T πi = shutdown cost for CHP unit i ∈ I ρ−i = ramp down limit for CHP unit i ∈ I (ρ−i < 0) ρ+i = ramp up limit for CHP unit i ∈ I (ρ+i > 0) ϕi = startup cost for CHP unit i ∈ I ψt = electricity purchase price from external source in time interval t ∈ T Continuous Variables

AUTHOR INFORMATION

Bt+1 = heat storage level before the beginning of time interval t ∈ T̅ Dt = heat disposed to external source in time interval t ∈ T E(i,t) = electricity production by CHP unit i ∈ I in t ∈ T̅ Pt = electricity acquired from external source in time interval t ∈T Qt = heat acquired from external source in time interval t ∈ T R(i,t) = electricity production change in CHP unit i ∈ I in t ∈ T̅ (t < τ) Wt = electricity sold to external source in time interval t ∈ T

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge the financial support from the Engineering and Physical Sciences Research Council (EPSRC) under the Research Project EP/G059071/1. The authors are also grateful for the financial support from the European Research Council (MOBILE, ERC Advanced Grant, No: 226462).

Binary Variables

S(i,t) = 1, if CHP unit i ∈ I startups at the beginning of time interval t ∈ T; i.e., X(i,t−1) = 0 and X(i,t) = 1 4384



Article

F(i,t) = 1, if CHP unit i ∈ I shutdowns at the beginning of time interval t ∈ T; i.e., X(i,t−1) = 1 and X(i,t) = 0 X(i,t) = 1, if CHP unit i ∈ I is operating in time point t ∈ T̅ Yt = 0, if there are no heat purchases in time interval t ∈ T Zt = 0, if there are no electricity purchases in time interval t ∈ T S̃(i,t,n) = lifting-state variable for minimum running time history F̃(i,t,n) = lifting-state variable for minimum shutdown time history

■

(21) Perea, E.; Grossmann, I.; Ydstie, E.; Tahmassebi, T. Dynamic modeling and classical control theory for supply chain management. Comput. Chem. Eng. 2000, 24, 1143−1149. (22) Perea-López, E.; Ydstie, B. E.; Grossmann, I. E. A model predictive control strategy for supply chain optimization. Comput. Chem. Eng. 2003, 27, 1201−1218. (23) You, F. Q.; Wassick, J. M.; Grossmann, I. E. Risk management for a global supply chain planning under uncertainty: Models and algorithms. AIChE J. 2009, 55, 931−946. (24) Liu, S.; Shah, N.; Papageorgiou, L. G. Multiechelon supply chain planning with sequence-dependent changeovers and price elasticity of demand under uncertainty. AIChE J. 2012, 58, 3390−3403. (25) Subramanian, K.; Rawlings, J. B.; Maravelias, C. T.; FloresCerrillo, J.; Megan, L. Integration of control theory and scheduling methods for supply chain management. Comput. Chem. Eng. 2013, 51, 4−20. (26) Ovacik, I. M.; Uzsoy, R. Rolling horizon procedures for dynamic parallel machine scheduling with sequence-dependent setup times. Int. J. Prod. Res. 1995, 33, 3173−3192. (27) Chand, S.; Traub, R.; Uzsoy, R. Rolling horizon procedures for the single machine deterministic total completion time scheduling problem with release dates. Ann. Oper. Res. 1997, 70, 115−125. (28) Bassett, M. H.; Pekny, J. F.; Reklaitis, G. V. Decomposition techniques for the solution of large-scale scheduling problems. AIChE J. 1996, 42, 3373−3387. (29) Dimitriadis, A. D.; Shah, N.; Pantelides, C. C. RTN-based rolling horizon algorithms for medium term scheduling of multipurpose plants. Comput. Chem. Eng. 1997, 21, S1061−S1066. (30) Janak, S. L.; Floudas, C. A.; Kallrath, J.; Vormbrock, N. Production scheduling of a large-Scale industrial batch plant. I. short-term and medium-term scheduling. Ind. Eng. Chem. Res. 2006, 45, 8234−8252. (31) Erdirik-Dogan, M.; Grossmann, I. E. Planning models for parallel batch reactors with sequence-dependent changeovers. AIChE J. 2007, 53, 2284−2300. (32) Shaik, M. A.; Floudas, C. A.; Kallrath, J.; Pitz, H.-J. Production scheduling of a large-scale industrial continuous plant: Short-term and medium-term scheduling. Comput. Chem. Eng. 2009, 33, 670−686. (33) Lima, R. M.; Grossmann, I. E.; Jiao, Y. Long-term scheduling of a single-unit multi-product continuous process to manufacture high performance glass. Comput. Chem. Eng. 2011, 35, 554−574. (34) Li, J.; Xiao, X.; Tang, Q.; Floudas, C. A. Production scheduling of a large-scale steelmaking continuous casting process via unit-specific event-based continuous-time models: short-term and medium-term scheduling. Ind. Eng. Chem. Res. 2012, 51, 7300−7319. (35) Pistikopoulos, E. N.; Georgiadis, M. C.; Dua, V. Multi-Parametric Programming: Theory, Algorithms, and Applications; Process Systems Engineering; Wiley-VCH: Weinheim, Germany, 2007; Vol. 1. (36) Pistikopoulos, E. N.; Georgiadis, M. C.; Dua, V. Multi-Parametric Model-Based Control: Theory and Applications; Process Systems Engineering; Wiley-VCH: Weinheim, Germany, 2007; Vol. 2. (37) Wallace, S. W. Decision making under uncertainty: Is sensitivity analysis of any use? Oper. Res. 2000, 48, 20−25. (38) Gal, T.; Nedoma, J. Multiparametric linear programming. Manage. Sci. 1972, 18, 406−422. (39) Gal, T. RIM multiparametric linear programming. Manage. Sci. 1975, 21, 567−575. (40) Pertsinidis, A. On the parametric optimization of mathematical programs with binary variables and its application in the chemical engineering process synthesis. Ph.D. Thesis, Carnegie Mellon University: Pittsburgh, PA, 1992. (41) Pertsinidis, A.; Grossmann, I. E.; McRae, G. J. Parametric optimization of MILP programs and a framework for the parametric optimization of MINLPs. Comput. Chem. Eng. 1998, 22, S205−S212. (42) Acevedo, J.; Pistikopoulos, E. N. A multiparametric programming approach for linear process engineering problems under uncertainty. Ind. Eng. Chem. Res. 1997, 36, 717−728. (43) Dua, V.; Pistikopoulos, E. N. An algorithm for the solution of multiparametric mixed integer linear programming problems. Ann. Oper. Res. 2000, 99, 123−139.

REFERENCES

(1) Sabuncuoglu, I.; Karabuk, S. Rescheduling frequency in an FMS with uncertain processing times and unreliable machines. J. Manuf. Syst. 1999, 18, 268−283. (2) Pistikopoulos, E. N. Uncertainty in process design and operations. Comput. Chem. Eng. 1995, 19, 553−563. (3) Li, Z.; Ierapetritou, M. G. Process scheduling under uncertainty: Review and challenges. Comput. Chem. Eng. 2008, 32, 715−727. (4) Aytug, H.; Lawley, M. A.; McKay, K.; Mohan, S.; Uzsoy, R. Executing production schedules in the face of uncertainties: a review and some future directions. Eur. J. Oper. Res. 2005, 161, 86−110. (5) Verderame, P. M.; Elia, J. A.; Li, J.; Floudas, C. A. Planning and scheduling under uncertainty: a review across multiple sectors. Ind. Eng. Chem. Res. 2010, 49, 3993−4017. (6) Sahinidis, N. V. Optimization under uncertainty: State-of-the-art and opportunities. Comput. Chem. Eng. 2004, 28, 971−983. (7) Subramanian, K.; Maravelias, C. T.; Rawlings, J. B. A state-space model for chemical production scheduling. Comput. Chem. Eng. 2012, 47, 97−110. (8) Zhuge, J.; Ierapetritou, M. G. Integration of scheduling and control with closed loop implementation. Ind. Eng. Chem. Res. 2012, 51, 8550− 8565. (9) Rodrigues, M. T. M.; Gimeno, L.; Passos, C. A. S.; Campos, M. D. Reactive scheduling approach for multipurpose chemical batch plants. Comput. Chem. Eng. 1996, 20, S1215−S1220. (10) Fang, J.; Xi, Y. G. A rolling horizon job shop rescheduling strategy in the dynamic environment. Int. J. Adv. Manuf. Technol. 1997, 13, 227− 232. (11) Sand, G.; Engell, S. Modeling and solving real-time scheduling problems by stochastic integer programming. Comput. Chem. Eng. 2004, 28, 1087−1103. (12) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous scheduling and control of multiproduct continuous parallel lines. Ind. Eng. Chem. Res. 2010, 49, 7909−7921. (13) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and control of tubular reactors: Parallel production lines. Ind. Eng. Chem. Res. 2011, 50, 8086−8096. (14) Novas, J. M.; Henning, G. P. A comprehensive constraint programming approach for the rolling horizon-based scheduling of automated wet-etch stations. Comput. Chem. Eng. 2012, 42, 189−205. (15) Sand, G.; Engell, S.; Märkert, A.; Schultz, R.; Schulz, C. Approximation of an ideal online scheduler for a multiproduct batch plant. Comput. Chem. Eng. 2000, 24, 361−367. (16) van den Heever, S. A.; Grossmann, I. E. A strategy for the integration of production planning and reactive scheduling in the optimization of a hydrogen supply network. Comput. Chem. Eng. 2003, 27, 1813−1839. (17) Wu, D.; Ierapetritou, M. G. Hierarchical approach for production planning and scheduling under uncertainty. Chem. Eng. Process. 2007, 46, 1129−1140. (18) Sung, C.; Maravelias, C. T. An attainable region approach for production planning of multiproduct processes. AIChE J. 2007, 53, 1298−1315. (19) Verderame, P. A.; Floudas, C. A. Integrated operational planning and medium-term scheduling for large-scale industrial batch plants. Ind. Eng. Chem. Res. 2008, 47, 4845−4860. (20) Li, Z.; Ierapetritou, M. G. Rolling horizon based planning and scheduling integration with production capacity consideration. Chem. Eng. Sci. 2010, 65, 5887−5900. 4385



Article

(44) Li, Z.; Ierapetritou, M. G. A new methodology for the general multiparametric mixed-integer linear programming (MILP) problems. Ind. Eng. Chem. Res. 2007, 46, 5141−5151. (45) Faísca, N. P.; Kosmidis, V. D.; Rustem, B.; Pistikopoulos, E. N. Global optimization of multi-parametric MILP problems. J. Global Optim. 2009, 45, 131−151. (46) Wittmann-Holbein, M.; Pistikopoulos, E. N. A two-stage method for the approximate solution of general multiparametric mixed-integer linear programming problems. Ind. Eng. Chem. Res. 2012, 51, 8095− 8107. (47) Pistikopoulos, E. N.; Dominguez, L.; Panos, C.; Kouramas, K.; Chinchuluun, A. Theoretical and algorithmic advances in multiparametric programming and control. Comput. Manage. Sci. 2012, 9, 183−203. (48) Ryu, J.-H.; Pistikopoulos, E. N. A novel approach to scheduling of zero-wait batch processes under processing time variations. Comput. Chem. Eng. 2007, 31, 101−106. (49) Ryu, J.-H.; Dua, V.; Pistikopoulos, E. N. Proactive scheduling under uncertainty: a parametric optimization approach. Ind. Eng. Chem. Res. 2007, 46, 8044−8049. (50) Li, Z.; Ierapetritou, M. G. Process scheduling under uncertainty using multiparametric programming. AIChE J. 2008, 53, 3183−3203. (51) Li, Z.; Ierapetritou, M. G. Reactive scheduling using parametric programming. AIChE J. 2008, 54, 2610−2623. (52) Wittmann-Holbein, M.; Pistikopoulos, E. N. Proactive scheduling of batch processes by a combined robust optimization and multiparametric programming approach. AIChE J. 2013, 59, 4184−4211. (53) Mayne, D. Q.; Rawlings, J. B.; Rao, C. V.; Scokaert, P. O. M. Constrained model predictive control: Stability and optimality. Automatica 2000, 36, 789−814. (54) Kopanos, G. M.; Capón-García, E.; Espuña, A.; Puigjaner, L. Costs for rescheduling actions: A critical issue for reducing the gap between scheduling theory and practice. Ind. Eng. Chem. Res. 2008, 47, 8785−8795. (55) Collazos, A.; Maréchal, F.; Gähler, C. Predictive optimal management method for the control of polygeneration systems. Comput. Chem. Eng. 2009, 33, 1584−1592. (56) El-Sharkh, M. Y.; Rahman, A.; Alam, M. S. Short term scheduling of multiple grid-parallel PEM fuel cells for microgrid applications. Int. J. Hydrogen Energy 2010, 35, 11099−11106. (57) Shaneb, O. A.; Taylor, P. C.; Coates, G. Optimal online operation ̀ of residential iCHP systems using linear programming. Energy Buildings 2012, 44, 17−25. (58) Wakui, T.; Yokoyama, R.; Shimizu, K. Suitable operational strategy for power interchange operation using multiple residential SOFC (solid oxide fuel cell) cogeneration systems. Energy 2010, 35, 740−750. (59) Kopanos, G. M.; Georgiadis, M. C.; Pistikopoulos, E. N. Energy production planning of a network of micro combined heat and power generators. Appl. Energy 2013, 102, 1522−1534. (60) Liu, P.; Pistikopoulos, E. N.; Li, Z. An energy systems engineering approach to the optimal design of energy systems in commercial buildings. Energy Policy 2010, 38, 4224−4231. (61) Weber, C.; Shah, N. Optimisation based design of a district energy system for an eco-town in the United Kingdom. Energy 2011, 36, 1292− 1308. (62) Mehleri, E. D.; Sarimveis, H.; Markatos, N. C.; Papageorgiou, L. G. Optimal design and operation of distributed energy systems: Application to Greek residential sector. Renew. Energy 2013, 51, 331− 342.

4386


Reactive Scheduling by a Multiparametric Programming Rolling

Recommend Documents