A Decomposition Scheme for Integration of Production Scheduling

One views demand response in industrial processes as the management of operating policies to accommodate variations in the availability and pricing of...
0 downloads 8 Views 2MB Size
Download PDF
Article pubs.acs.org/IECR

A Decomposition Scheme for Integration of Production Scheduling and Control: Demand Response to Varying Electricity Prices Chudong Tong,*,† Ahmet Palazoglu,‡ and Nael H. El-Farra‡ †

Faculty of Electrical Engineering & Computer Science, Ningbo University, Zhejiang, 315211, P. R. China Department of Chemical Engineering and Materials Science, University of California, Davis, One Shields Avenue, Davis, California 95616, United States

‡

ABSTRACT: One views demand response in industrial processes as the management of operating policies to accommodate variations in the availability and pricing of electricity from the grid. To find the optimal operating regimes and to facilitate the transitions in between such regimes in a costeffective manner, a mixed-integer nonlinear programming problem is often formulated. A full solution to this problem involves not only finding the optimal operating regimes but also determining the optimal control policies to ensure minimal transition costs; such solutions are often mathematically intractable. In this paper, a decomposition strategy is proposed to arrive at a feasible, albeit suboptimal, solution. The scheduling problem is first solved using a prefixed transition time to minimize the total operating cost; only the first scheduled transition of the obtained solution is implemented in the process, and optimization of controller parameters is then carried out to economically shape the transient profiles. When the first scheduled transition with the optimized controller parameters is completed, the scheduling and control problems are reformulated and resolved again in a receding-horizon optimization manner. A conceptual case study of a continuous stirred tank reactor demonstrates the key contributions and effectiveness of the strategy.

1. INTRODUCTION As the availability and pricing of electricity vary during the day due to intermittency in renewable resources or real-time electricity tariff policies, the chemical industry, especially those with power-intensive processes, need to rethink operating strategies to maximize economic viability. The implementation of demand response (DR) strategies, which focus on demandside activities, has been at the center of all pragmatic energy policy decisions and is the key to minimizing production costs and maximizing profits. The economic potential of DR for industrial processes has been acknowledged by numerous studies.1−4 In developing DR strategies, the target lies in scheduling the production rates, production times, and inventory levels at each time step to achieve specific objectives. Therefore, process-control-related activities are closely linked to the responsive demand augmentations, given that control actions involve real-time manipulation of selected process variables to hold key product quality and production rates near their desired target values.5,6 Moreover, as the plant operating schedule varies, process control facilitates the dynamic transition periods between different operating modes. Therefore, it becomes imperative to develop operating methodologies that can effectively integrate production scheduling and control functions. There have been several early studies addressing the integrated scheduling and control problems. The approaches in the existing literature can be categorized into (i) simultaneous modeling and (ii) decomposition methods.7 Mahadevan et al.8 solved the control problem under different scheduling policies for polymer processes by analyzing the © XXXX American Chemical Society

grade transition scheduling problem from a robust closed-loop control point of view. They obtained grade schedules by defining easy and hard to carry out transitions, respectively. However, the problem was not attempted as a mixed-integer dynamic optimization (MIDO) problem. Flores-Tlacuahuac and Grossmann9 proposed a simultaneous scheduling and control formulation by explicitly incorporating into the scheduling model process dynamics in the form of differential/algebraic constraints. Although the original simultaneous scheduling and control problems cast as a MIDO problem is transformed into a mixed-integer nonlinear programming (MINLP) problem, the authors did not address the problem in the closed-loop context (e.g., no control system was present). Zhuge and Ierapetritou7 extended the work of Flores-Tlacuahuac and Grossmann9 and introduced a closedloop implementation to guarantee effective responses to process disturbances. Nevertheless, optimizing the dynamic profile of the transitions between different operating modes through controller parameter tuning was not considered in these two studies. Recently, Chu and You10 proposed a novel integration method to solve the scheduling and control problem simultaneously. In their work, a set of controller candidates is generated offline for all available transitions, and the online control-related implementation is carried out by selecting the optimal controller instead of real-time optimizing Received: Revised: Accepted: Published: A

March 1, 2017 June 19, 2017 July 11, 2017 July 12, 2017 DOI: 10.1021/acs.iecr.7b00869 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research controller parameters. From this viewpoint, the problem of controller parameter tuning is indeed not implemented with the scheduling problem simultaneously. Additionally, Chu and You11 presented a moving horizon approach of integrating scheduling and control; the integrated problem is solved to determine the reference trajectories obtained beforehand for advanced controllers. Furthermore, some review papers are available in the literature discussing the integration of scheduling and control.12−14 On the other hand, in the class of decomposition methods, the integration of scheduling and control can be modeled by formulating two subproblems: a mixed-integer linear/nonlinear programming problem (scheduling) and a dynamic optimization problem (control). Nystrom et al.15 iterated between these two subproblems and updated the solution until a target gap solution is obtained. The same authors later applied the same approach to parallel polymerization lines with multiple units.16 However, this approach is only feasible if a small number of production transitions must be considered. An extensive literature review on the solution of scheduling and control problems can be found elsewhere.17,18 In our previous work, a DR problem was posed through a MINLP formulation, where optimal operating policies (modes) are found in response to varying electricity prices.19 Our key assumption in that paper was that the control policies were predetermined and no attempt was made to optimize the transition periods. By relaxing that assumption, the goal of the present work is to offer a new DR formulation for the integration of scheduling and control under the consideration of time-varying electricity prices. As will be stated in the next section, the formulation of the problem results in a rather complex decision-making problem, which is difficult to solve in a single formulation. To tackle this task, a decomposition strategy is proposed to simplify the modeling effort and generate a feasible solution. While one of the advantages of decomposition methods lies in the relative ease of solving simpler subproblems, the main disadvantage is that such a procedure generally leads to suboptimal solutions.17 It should be stressed here that the special simultaneous scheduling and control problem considered in the current work cannot be easily solved through a single optimization formulation. Furthermore, the incorporation of energy constraints related to electricity pricing places additional difficulty in solving the problem. The so-called “suboptimal” solution that results from the proposed decomposition strategy is the only feasible solution that can be obtained. This point will be articulated further in the paper, along with the detailed formulation of the proposed decomposition strategy.

Figure 1. A schematic visualization of the problem statement.

with the production rates. The problem is to determine production rates and storage levels on an hourly basis and tune the controller parameters to improve the dynamic transition performance (i.e., mainly concerning transition time and actuator dynamics), so that the objective of minimizing the total operating costs (mostly due to energy costs) is achieved. At this point, it should be emphasized that there is no previous work that simultaneously addresses the scheduling of a whole range of production rates (i.e., the manufacturing facility can be operated at any continuous production rate as long as it is bounded) and the optimization of the controller parameters to optimally shape dynamic transition profiles between different production rates. The objective is to guarantee more economical process operations through interaction between scheduling and control problems. Moreover, the formulated model also needs to satisfy the hourly demand of product with the desired quality and to minimize the total production cost (i.e., energy cost for both steady-state and transition periods as well as raw material waste during transitions). However, given that the previously studied production scheduling and control problems are dealing with different operating levels, it is not easy to integrate them effectively. A major challenge is the computational intractability of a single optimization model dealing with simultaneous scheduling and control. The scheduling problem is solved to define the hourly production rates, while the control-related problem considering the transition dynamics, on the other hand, tends to be more detailed in nature. If a single optimization model is formulated with a detailed time scale considered, the resulting model would come with an intractable size. In addition, a key aspect of controller tuning is the consideration of stability. The objective of integrating scheduling and control, however, mainly focuses on process economics. Therefore, it is difficult to embed the objective of control performance within the final economical objective function. To overcome the above difficulties, a decomposition scheme is proposed to address the integration problem of scheduling and control, and a receding-horizon optimization strategy is employed to get the solution at each scheduling step. 2.2. Solution Overview. In the current work, a conceptual case study will be later carried out on a nonisothermal continuous stirred tank reactor (CSTR) where the hourly energy load (Wh) is assumed to be roughly proportional to the inlet feed flow rate (Fh) (or production rate in this case). Although the CSTR is not necessarily considered as a powerintensive unit in actual operation, it is still found useful as a way of illustrating the proposed problem formulation conceptually. As shown in Figure 2, the relationship between energy load (Wh) and production rate (Fh) can be expressed by a linear relationship

2. FORMULATION OVERVIEW 2.1. Generic Problem Definition. We start with the manufacture of a single product on a single continuously operating unit (power-intensive unit) as shown in Figure 1. The unit is required to satisfy hourly demand for the product that is specified as a constant rate, and a storage tank provides flexibility and feasibility for the DR strategy. It is assumed that the production costs vary every hour and are related to hourly varying electricity prices, which are assumed to be available for a typical week. DR is achieved through mechanisms that discourage the energy load when the real-time electricity price is high and vice versa. Therefore, a range of continuous production rates with upper and lower limits needs to be specified a priori given that the energy load is closely related

Wh = α(Fh − β) B

(1) DOI: 10.1021/acs.iecr.7b00869 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Remark 1: It should be noted that the prespecified fixed transition time (Ts) in step 1 plays an important role in the overall decomposition framework. If Ts is chosen to be too small, it would lead to an infeasible solution given that the actual storage level and the loss of production associated with the optimized transition time in step 2 cannot manage to meet the hourly demand. On the other hand, with a large value of Ts, the resulting solutions from the scheduling step might be too aggressive and lead to an inferior suboptimal final solution. Therefore, we propose to heuristically determine Ts with the consideration of a possible maximum transition time. In addition, the existing optimized transition times are not suggested to be used in updating the prespecified transition time Ts, given that the optimized times might be rather too small or too large for future scheduled transitions. Therefore, the prespecified transition time Ts is not updated or corrected by one of the existing optimized transition times obtained from step 2. Remark 2: It is assumed that the steady-state periods yield product(s) with the desired quality specifications (on-spec). The transition period results are due to the scheduling of production rates and cannot be avoided even if the operation yields off-spec product(s). From this point of view, it is imperative to minimize the transition time with respect to the minimization of production costs. With the actual optimized transition time, the corresponding storage level and loss of production can then be updated, and a subsequent scheduling formulation based on the new updated conditions is carried out for better economic operating decisions. Hence, information is shared between the scheduling and control subproblems in the proposed approach. Remark 3: Unlike other previous methods addressing the scheduling and control problems sequentially, the proposed decomposition scheme works by using a receding-horizon strategy. First, the optimization model in scheduling step is regarded as a predictive model with a constant predicted transition time Ts, and the remaining hours of the scheduled week is considered as the prediction horizon, which shrinks as time progresses. Second, the deviations of storage level as well transition time are determined for the first scheduled transition and then new solutions are generated in the case of updated variables using the predictive model. Notably, only the first scheduled transition is realized and the storage level associated with the optimized transition time is updated. This sort of integration proposed in this work results in more economical process operations, which will be illustrated in section 4.

Figure 2. Conceptual representation of the power relationship.

where α and β are constant parameters. The formulation of eq 1 is motivated by the so-called Willan’s line, which is applicable for modeling power output of turbines.20−22 As analyzed in the previous subsection, the present problem would result in a complex and unsolvable MINLP problem if it is formulated via a single model. One feasible way to solve this problem is to decompose it into two subproblems, a production scheduling problem (MINLP) and a control problem (NLP), and the interaction between the scheduling and control subproblems should also be considered to ensure economically meaningful process operation. The proposed framework is depicted in Figure 3, and a detailed description is given here:

Figure 3. Flowchart of the proposed decomposition scheme.

3. MATHEMATICAL PROGRAMMING MODELS In this section, the detailed mathematical models for both scheduling and control subproblems will be discussed. For convenience, the starting time point is specified at the beginning of each week and a continuous stable production rate is assumed to be reached before the DR is implemented. The definition of symbols can be found in the Nomenclature section. 3.1. Production Scheduling Model. As stated before, the goal of the scheduling problem is to improve the operational and economic objectives. Therefore, the objective function is defined as follows

(1) The scheduling problem is solved using a prespecified transition time (Ts) between different production rates to minimize the total operating cost, which consist of energy expense, inventory cost, and raw material waste. (2) The first scheduled transition of the solution obtained from step 1 is implemented in the process, and optimization of controller parameters is then carried out between the set of two different production rates with the objective of minimizing transition time. (3) With the optimized transition time, the actual storage level and the loss of production for the hour of the scheduled transition are updated. (4) The hour of the first scheduled transition is set to be the new initial time instant, and the scheduling problem is reformulated with the initial condition calculated from step 3. (5) Steps 1−4 are repeated until the final solution is obtained on an hourly basis for a whole week.

Jschedule = min{Φ1 + Φ2 + Φ3}

(2)

with C

DOI: 10.1021/acs.iecr.7b00869 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research H

Φ1 =

H

∑ Whδhel =

∑ α(Fh − β)δhel

h=1

h=1

globally optimal solutions to general nonlinear problems with continuous and/or discrete variables. 3.2. Controller Parameter Optimization Model. Among numerous control methods, the proportional-integral-derivative (PID) control algorithm remains the most popular approach at the regulatory level due to its simple structure and robust performance in the vast majority of applications.24,25 In this work, a PI controller is adopted to hold a key product quality for the CSTR unit, and the parameters of proportional and integral terms are optimally tuned by taking into account the economic process operation. 3.2.1. Performance Criterion. In general, it is critical to have a good estimate for the performance criterion. There are several performance criteria for tuning PID controllers, such as the integral of absolute error (IAE), the integral of squared error (ISE), and the integral of time-weighted squared error (ITSE).26,27 However, they are not very useful for our purpose of minimizing transition time or other relevant objectives. In this paper, a new performance criterion is proposed for tuning a PI controller that can be computed numerically

(3)

H

Φ2 =

raw ∑ zhTF s hδ h=1

(4)

H

Φ3 =

∑ Shδ S h=1

(5)

representing the energy cost, cost of raw material waste, and inventory cost, respectively. It can be seen that the scheduling formulation is a discrete time model with a granularity of 1 h. In practice, the manufacturing facility will have also a constraint on the processing capacity, which is given as

Fmin ≤ Fh ≤ Fmax

(6)

To model the transitional constraints efficiently, the binary variable zh, regarded as a transition indicator, is introduced ⎧1, if |Fh − Fh − 1| > 0 zh = ⎨ ⎩ 0, otherwise ⎪

⎪

Jcontrol = ts + ϕu (7)

with

where zh will be 1 if there is a change in the production rate from time step h − 1 to h. However, from a control perspective, the rate of change of production rates might also be restricted. The following constraints can be formulated in terms of the transition indicator zh and production rate Fh zhΔFmin ≤ |Fh − Fh − 1|

(8)

|Fh − Fh − 1| ≤ zhΔFmax

(9)

ts = max flag(k) K

ϕu =

(15)

where ts is the transition time, ϕu is defined to ensure the stable action of the actuator, K is chosen sufficiently large so that Δu(k) is negligible after sampling time K, and t(k) is the sampling time at sampling step k. The flag(k) is a discontinuous function defined to represent the transition time, which is given as ⎧ k Δt , if |x − x ̅ | < ε flag(k) = ⎨ otherwise ⎩ 0,

(10)

(16)

where x and x̅ are product-quality indices under transient and steady-state conditions, respectively. The transition process is defined as the time period when no on-spec product is produced constantly. 3.2.2. Dynamic Model Discretization. To address the optimal tuning of controller parameters, the dynamic model representing the system behavior is discretized using the Runge−Kutta methods.28 Let the first principal model of the given plant be specified as follows: xi̇ = fi (t ,x1 ,x 2 ,...,xn ,u)

(11)

Given that satisfying the hourly demand is a hard constraint, the storage tank is initially charged with a certain amount of onspec product S0 at the beginning of the weekly schedule so as to make the DR problem feasible. To avoid depleting the product in the storage level at end of the week, an additional constraint is placed on the storage level at the last hour:

SH ≥ S0

∑ [|Δu(k)| t(k)] k=1

where S̃h is an auxiliary variable for Sh, which is defined to take the value assigned to the same variable but one time step backward; i.e., S̃h = Sh−1. The initial storage level at the beginning of each scheduling formulation is thus assigned to S̃h at the time index h = 1. Note that the term (1 − zhTs)Fh in eq 10 indicates that a certain amount of off-spec product is generated and sent to waste when zh = 1. Additionally, the storage level is constrained by the storage capacity:

0 ≤ Sh ≤ Smax ∀ h

(14)

k = 1,2,... K

Note that the constraint in eq 8 can prevent the whole system from being highly oscillatory, and a large change in production rates is prohibited by constraint in eq 9. A mass balance has to be enforced to describe the relationship between production rates (Fh), storage levels (Sh), and the demand (d) Sh = Sh̃ + (1 − zhTs)Fh − d

(13)

i = 1, 2, ..., n

(17)

The state of the next sampling step is calculated as xi(k+1) = xi(k) +

Δt (k1 + 2k 2 + 2k 3 + k4) 6

(18)

with

(12)

The problem of production scheduling is naturally cast as a MINLP problem. The solution is generated by the LINDOGLOBAL solver implemented in GAMS23 given that it supports nonsmooth functions like absolute value, or even discontinuous mathematical functions, and finds guaranteed D

k1 = fi (t(k),x1(k),xi(k),...,xn(k),u(k))

(19)

⎛ ⎞ Δt Δt k 2 = fi ⎜t(k)+ ,x1(k),xi(k)+ k1 ,...,xn(k),u(k)⎟ ⎝ ⎠ 2 2

(20)

⎛ ⎞ Δt Δt k 3 = fi ⎜t(k)+ ,x1(k),xi(k)+ k 2 ,...,xn(k),u(k)⎟ ⎝ ⎠ 2 2

(21)

DOI: 10.1021/acs.iecr.7b00869 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research k4 = fi (t(k)+Δt ,x1(k),xi(k)+Δtk 3 ,...,xn(k),u(k))

optimization of PI controller parameters aims to improve the dynamic profile of the transient process with respect to minimizing the transition cost. Without optimized parameters, the PI controller itself can still achieve closed-loop control performance, and this is particularly important for online controller parameter tuning if the optimized parameters cannot be reached on time. Remark 4: On the basis of the mathematical formulations given above, the resulting MINLP model has a complex structure and is inefficient for solving problems of practical interest if only a single optimization model is formulated. To overcome computational deficiencies in the original complex problem, the proposed decomposition scheme is decomposed in a sequential manner, as shown in Figure 3. The advantage of such a strategy is that a series of simpler subproblems is solved instead of one complex problem. Additionally, the proposed interaction between scheduling and control enables us to generate a feasible solution but, more importantly, enables us to solve the problem more economically. This will be illustrated in the next section.

(22)

where k1, k2, k3, k4 are intermediate increments that are used in eq 18. 3.2.3. PI Controller. The velocity form of the discrete approximation of a PI controller is given by29 Δu(k) = k p[e(k) − e(k − 1)] + kIΔte(k)

(23)

where e(k) = r(k) − y(k) is the error and kP and kI are the proportional and the integral parameters, respectively. The control input is thus calculated as u(k) = u(k − 1) + Δu(k)

(24)

Considering the dynamics of the actuator, the actual input is bounded as umin ≤ u(k) ≤ umax

(25)

3.2.4. Ensemble Strategy of Initialization. The optimization of PI controller parameters is cast as a NLP problem, which includes nonsmooth and discontinuous functions shown above. It is important to note that, however, kP and kI are the only free decisions and they directly influence other variables. Given that there are plenty of intermittent variables and equality relations in this formulation, we can define some user-defined functions to avoid this problem. The number of decision variables is thus reduced to four, that is, kP, kI, ts, and ϕu. The NLP problem developed in this work is solved using NLOPT solver implemented in MATLAB.30 However, the NLP problem is sensitive to the initial point, and in this work, we propose an ensemble initialization strategy to increase the probability of reaching the global optimum. To account for system stability, the controller parameters kP and kI are assumed to lie in the interval [kPmin,kPmax] and [kImin,kImax], respectively. The intervals for controller parameters can be determined by the stability analysis of the control system or by prior knowledge provided by process operators. Even if the determined intervals for controller parameters accidentally cover some realizations that lead to an unstable system, the optimization of controller parameters could still search for a realization that ensures the stability of the controlled system. The visualization of the ensemble strategy is depicted in Figure 4, where the search regions are divided into three subregions,

4. CONCEPTUAL CASE STUDY 4.1. Plant Description and Component Data. To motivate the discussion, consider a single product to be manufactured by the reaction A → B in a single CSTR as shown in Figure 5. The model equations are given as follows: dCA F = f1 (CA ,T ,Tj) = (CAf − CA ) − k 0e−ΔE / RT CA dt V (26)

dT = f2 (CA ,T ,Tj) dt UAc F ΔH −ΔE / RT k 0e CA − (T − Tj) = (Tf − T ) − V VρCp ρCp (27)

dTj dt

= f3 (CA ,T ,Tj) =

Fj Vj

(Tjf − Tj) +

UA c (T − Tj) VjρjCp j (28)

Figure 4. Visualization of the ensemble initialization strategy.

Figure 5. CSTR process with the PI controller.

and nine different starting points are obtained. The final solution is determined by the one that achieves the minimum value of the objective function defined in eq 13. It should be emphasized that there is no dominant rule for guiding the definition of search regions; the number of search regions is determined by a trade-off consideration between the calculation burden and the “goodness” of the final solution. Moreover, the

The inlet concentration C Af = 8 mol/m 3 and the concentration of the product is within 2% of the steady-state value C̅ A = 0.6 mol/m3. A constraint on the plant capacity is placed as 0.5 m3/h ≤ F ≤ 1.2 m3/h. The values of the steadystate conditions for different production rates (i.e., feed flow rate F) can be calculated by forcing eqs 26−28 to be 0. The sampling interval of the CSTR is Δt = 1/60 h. Note that the E

DOI: 10.1021/acs.iecr.7b00869 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Table 1. Electricity Prices for the Given Week in dollars/MWh Monday 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

74.47 61.67 54.24 50.16 51.65 60.61 73.33 78.52 89.19 106.38 126.38 142.76 156.82 170.55 182.34 199.88 209.06 199.48 166.29 146.04 151.89 130.88 97.82 84.87

Tuesday 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

73.76 60.81 54.12 49.5 50.36 59.69 75.07 78.9 88.78 103.88 120.89 140.41 153.94 169.68 184.76 199.18 208.01 197.97 159.68 137.65 141.29 124.22 93.38 83.82

Wednesday 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72

71.59 59.65 52.42 48.43 49.03 58.47 75.26 76.61 83.83 95.67 110.56 126.24 136.1 148.15 160.21 173.83 179 173.2 143.2 122.91 132.37 116.02 85.15 79.99

Thursday 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

Friday

71.19 59.56 52.94 48.94 50.05 58.71 73.37 78.01 85.49 98.56 113.3 129.73 138.12 152.34 164.73 179.33 186.55 176.83 146.87 127.28 134.18 120.52 88.85 86.48

97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120

75.49 62.24 53.82 48.42 48.97 57.88 70.76 77.85 84.63 99.05 115.4 131.01 141.74 154.11 162.69 177.7 184.5 171.74 143.34 126.97 131.61 116.28 90.41 82.57

Saturday 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144

72.49 61.22 51.37 46.39 45.08 47.26 48.56 57.43 73.34 97.7 116.08 129.9 141.53 147.1 155.32 167.39 171.66 162.44 139.72 123.53 131.48 119.18 93.07 83.19

Sunday 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168

71.68 58.8 49.84 44.37 42.73 43.03 42.95 50.6 66.56 86.12 102.67 121.61 134.27 140.35 148.37 159.65 168.15 162.58 140.37 125.69 135.32 119.88 93.74 83.11

Figure 6. (Top) Electricity prices profiles for a given week. (Middle) Production rates profiles obtained from the decomposition scheme (black solid line) and the one with a fixing transition time (red dashed line). (Bottom) inventory profile from the decomposition scheme.

4.2. Results. The objective of this case study is to maintain an economical process operation (i.e., production rates and storage level at each time step) and optimal PI tuning parameters for all scheduled transitions. The prefixed transition time is taken as Ts = 0.5 h and a comparison study of different transition times will be discussed later. The proposed decomposition scheme is implemented within a recedinghorizon framework and only the first scheduled transition of

demand assumes a constant rate for the entire week, which is specified as d = 0.8 m3/h. The CSTR is operated at a production rate of F = 0.8 m3/h before the proposed DR strategy is initiated. The electricity prices varying on an hourly basis for a typical week are tabulated in Table 1, which are based on PJM data (www.pjm.com). Our production decisions do not influence the electricity prices. F

DOI: 10.1021/acs.iecr.7b00869 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 7. Transition profiles for the scheduled transitions at times 35 and 120 h, respectively.

Figure 8. (Top) Electricity prices for a given week. (Bottom) Production rates profiles for different prespecified transition times, 0.5 h (black solid line), 0.6 h (red dashed−dotted line), and 0.7 h (blue dashed line), respectively.

Figure 9. Inventory profiles generated with Ts = 0.5 h (top), Ts = 0.6 h (middle), and Ts = 0.7 h (bottom), respectively.

the solution obtained from the scheduling phase is input to the process and the optimization of PI parameters is then initiated to shape the dynamic transition profile. The application of the decomposition scheme results in the plots of Figure 6. As one would expect, the solution is heavily influenced by the market prices of electricity, in which the process tends to be operated at high production rates during the off-peak hours and vice versa.

For comparison purposes, an additional case without considering transient process optimization is considered here (i.e., the transition time is fixed for all transitions). The compared case is just a scheduling problem, no recedinghorizon optimization is involved, and the corresponding solution is actually obtained from the first run of the scheduling problem in the proposed formulation. Given that the actual transitions implemented in the compared case are not G

DOI: 10.1021/acs.iecr.7b00869 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

10.53% compared with having no DR strategy. As Ts increases, the solution becomes “more suboptimal”. However, it is important to note that the economic gain is still impressive, even though the Ts is chosen to be a larger value, 0.7 h. Furthermore, we also studied two additional values of Ts, 0.4 and 0.8 h; the final solutions, however, became infeasible and come with no transition, respectively. With the underestimated Ts, some of the scheduled transitions cannot be realized with respect to satisfying the hourly demand hard constraint. While the process decides to stay constant at the nominal steady state (i.e., Fh = 0.8 m3/h), all the time since the transition is costly. As concluded from this comparison, the proposed decomposition method for the integration of scheduling and control problems is sensitive to the prefixed transition time, a possible way of specifying a “good” value for the transition time would be relying on process prior knowledge, and the value should neither be too large or too small.

considered in detail, the real loss of production cannot be reached and the inventory profile has no real meaning; it is thus not compared here. According to the middle plot of Figure 6, the curve is similar to the proposed decomposition case but a bit more exaggerated. This is because the fixed transition time (0.5 h) is overestimated compared to the optimized one. Given that the minimum and maximum change rates of production rates are restricted by eqs 8 and 9, the production rates cannot change from the highest level (1.2 m3/h) to the lowest level (0.5 m3/h) directly and vice versa (Figure 6). The inventory profile given at the bottom of Figure 6 clearly shows how the demand is satisfied from the storage tank during plant underutilization. To evaluate the achieved performance by the tuned PI controller, the transition profiles at time steps of 35 h (left plot) and 120 h (right plot) are depicted in Figure 7, where the production rates change from F34 = 1.065 86 m3/h to F35 = 0.565 86 m3/h and from F119 = 0.929 56 m3/h to F120 = 1.1994 m3/h, respectively. Note that the black dashed lines define the region of on-spec product. The tuned parameters of the PI controller for these two transition profiles are kP = −0.1546, kI = −2.4992 and kP = −0.2078, kI = −3.6291, respectively. 4.3. Comparisons. The main argument for the determination of the prefixed transition time Ts is that it is required to ensure the feasibility of the integration of scheduling and control with respect to meeting hourly demand. Without enough practical experience, it might be difficult to choose a meaningful value. As stated before, Ts can be neither too small nor too large. Therefore, a comparison study is presented here to illustrate the influence of different values of Ts, which are given as 0.5, 0.6, and 0.7 h. The solutions given in Figure 8 show a similar pattern with subtle differences. Within the time period from 105 to 125 h, the operation exhibits relatively more oscillations. This is mainly due to the suboptimal nature of the proposed decomposition scheme. Additionally, it is interesting to observe that the solution associated with Ts = 0.7 h tends to stay at each scheduled production rate for a relatively longer time given that it is more costly to make such a transition considered in the scheduling problem. The inventory profiles corresponding to different values of Ts are depicted in Figure 9, from which it is easy to observe how the storage tank is used to meet the hourly demand. The corresponding economic impact of the optimization strategy for different values of Ts are reported in Figure 10, in which the total costs are decreased by 13.00%, 12.28%, and

5. CONCLUSIONS A novel decomposition scheme for solving the integration of scheduling and control problems under the consideration of time-varying electricity prices has been presented. We described a problem formulation that simultaneously takes into account the production scheduling with respect to minimizing production cost, the optimization of controller parameters to realize the purpose of reducing the time and the lost product during transitions, and the integration of scheduling and control to achieve a more economical operating policy. Instead of formulating one complex and unsolvable MINLP problem, we proposed a novel decomposition scheme that generates smaller subproblems that can often be solved easily. From the point of view of scheduling, control is a necessary means to implement the planned schedules. If the transient process could better controlled, less transition time and lost product could be observed. In terms of process control in the scheduling problem, the schedule provides targets that may be more or less difficult to realize. The scheduling and control problems may be integrated to simultaneously optimize the process operation in an economic manner. The proposed strategy is implemented within a receding-horizon optimization framework, in which only the first transition is carried out with optimal tuned PI parameters and a subsequent scheduling problem is then reformulated with the new updated storage level corresponding to optimized transition time. This implementation was successfully applied to the CSTR process. The achieved economic impact in the CSTR process, i.e., more than 10% reduction in total operating cost compared to the operating strategy without DR, demonstrated the feasibility and validation of the proposed method. It must be stressed that the current work only considered the optimal tuning of controller parameters to improve transition profiles. However, real processes are subject to disturbances or faults, and it is necessary to develop a mechanism that can effectively respond to process disturbances and faults in the coming future. Furthermore, from the case study simulated in the current work, the ability to deal with larger dimensional systems featuring multiple units, advanced control strategies, and long-term resource planning needs to be tackled. In addition, another challenge is to meet all practical requirements and find appropriate applications in industry. Again, the decomposition strategy seems to be the direct solution of complex systems with the above-mentioned characteristics. Other future work will consider the effect of uncertainty in the

Figure 10. Comparison of total costs. H

DOI: 10.1021/acs.iecr.7b00869 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research problem data, different participation of electricity marketing prices, and the integration of planning, scheduling, and control problems.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: +86 18815280932. ORCID

Chudong Tong: 0000-0002-4973-8060 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was sponsored by the National Natural Science Foundation of China (61503204), the US National Science Foundation (CBET-1438456), the K.C. Wong Magna Fund in Ningbo University, the Natural Science Foundation of Zhejiang Province (Y16F030001 and LY14F030004), and the Natural Science Foundation of Ningbo City (2016A610092).

■

■

ΔE = activation energy, 3.1 × 1010 J/mol R = universal gas constant, 3780 J/mol·K T = reactor temperature, K Tf = temperature of reactor feed stream, 298 K ρ = density of reactor contents, 800 kg/m3 Cp = heat capacity of reactor contents, 2540 J/kg·K ΔH = heat of reaction, −3.1 × 107 J/mol U = heat transfer coefficient, 3 × 106 J/h·K·m3 Ac = area available for heat transfer, 20 m2 Tj = temperature of jacket outlet stream, K Fj = jacket flow rate, m3/h Vj = jacket volume, 0.2 m3 Tjf = temperature of jacket feed stream, 298 K ρj = density of jacket liquid, 1000 kg/m3 Cpj = heat capacity of jacket liquid, 4186 J/kg·K

REFERENCES

(1) Paulus, M.; Borggrefe, F. The potential of demand-side management in energy-intensive industries for electricity markets in Germany. Appl. Energy 2011, 88, 432. (2) 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. (3) Wang, Y.; Li, L. Time-of-use based electricity demand response for sustainable manufacturing systems. Energy 2013, 63, 233. (4) 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. (5) Mendoza-Serrano, D. I.; Chmielewski, D. J. Demand response for chemical manufacturing using economic MPC. In Proceedings of the 2013 American Control Conference; IEEE, 2013; p 6655. (6) Harjunkoski, I.; Nyström, R.; Horch, A. Integration of scheduling and control − theory or practice? Comput. Comput. Chem. Eng. 2009, 33, 1909. (7) Zhuge, J.; Ierapetritou, M. G. Integration of scheduling and control with closed loop implementation. Ind. Eng. Chem. Res. 2012, 51, 8550. (8) Mahadevan, R.; Doyle, F. J.; Allcock, A. C. Control-relevant scheduling of polymer grade transitions. AIChE J. 2002, 48, 1754. (9) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous cyclic scheduling and control of a multiproduct CSTR. Ind. Eng. Chem. Res. 2006, 45, 6698. (10) Chu, Y.; You, F. Integration of scheduling and control with online closed-loop implementation: Fast computational strategy and large-scale global optimization algorithm. Comput. Chem. Eng. 2012, 47, 248. (11) Chu, Y.; You, F. Moving horizon approach of integrating scheduling and control for sequential batch processes. AIChE J. 2014, 60, 1654−1671. (12) Chu, Y.; You, F. Model-based integration of control and operations: overview, chanllenges, advances, and opportuniti es. Comput. Chem. Eng. 2015, 83, 2. (13) Baldea, M.; Harjunkoski, I. Integrated production scheduling and process control: A systematic review. Comput. Chem. Eng. 2014, 71, 377. (14) Engell, S.; Harjunkoski, I. Optimal operation: Scheduling, advanced control and their integration. Comput. Chem. Eng. 2012, 47, 121−133. (15) Nystrom, R. H.; Harjunkoski, I.; Kroll, A.; Franke, R. Production campaign planning including grade transition sequencing and dynamic optimization. Comput. Chem. Eng. 2005, 29, 2163. (16) Nystrom, R. H.; Harjunkoski, I.; Kroll, A. Production optimization for continuously operated processes with optimal operation and scheduling of multiple units. Comput. Chem. Eng. 2006, 30, 392.

NOMENCLATURE

Indices

h = hourly time step indices (1, 2, ..., H) k = sampling step indices (1, 2, ..., K) Parameters

α = constant parameter in Willan’s line, 0.8 β = constant parameter in Willan’s line, 0.4 δelh = hourly electricity price, $/MW δraw = unit raw material price, 1.0 $/m3/h δS = unit storage price, 0.1 $/m3/h Ts = prefixed transition time, h Fmin = minimum production rate, 0.5 m3/h Fmax = maximum production rate, 1.2 m3/h ΔFmin = minimum change rate, 0.1 m3/h ΔFmax = maximum change rate, 0.5 m3/h d = hourly demand, 0.8 m3/h Smax = maximum storage capacity, 4m3 S0 = initial storage level at the beginning of the given week, 0.5 m3 Δt = sampling interval, 1 min umin = minimum manipulated input, 0.0001 m3/h umax = maximum manipulated input, 1 m3/h kPmin = minimum value of controller gain kPmax = maximum value of controller gain kImin = minimum value of controller integral term kImax = maximum value of controller integral term

Decision Variables

Fh = production rates on hourly basis zh = transition indicator Sh = storage level on hourly basis kP = parameter of proportional term kI = parameter of integral term ts = transition time ϕu = accumulated function of actuator action Other Symbols

F = inlet flow rate (production rate), m3/h V = reactor volume, 1 m3 CAf = concentration of species A in reactor feed stream, 8 mol/m3 CA = concentration of species A in reactor, mol/m3 k0 = pre-exponential factor, 7.08 × 1010 h−1 I

DOI: 10.1021/acs.iecr.7b00869 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research (17) Flores-Tlacuahuac, A.; Grossmann, I. E. Simultaneous scheduling and control of multiproduct continuous parallel lines. Ind. Eng. Chem. Res. 2010, 49, 7909. (18) Mitra, K.; Gudi, R.; Patwardhan, S.; Sardar, G. Resilency issues in integration of scheduling and control. Ind. Eng. Chem. Res. 2010, 49, 222. (19) Tong, C.; Palazoglu, A.; El-Farra, N. H.; Yan, X. Energy demand management for process systems through production scheduling and control. AIChE J. 2015, 61, 3756. (20) Willans, P. W. Economy trials of a non-condensing steamengine: simple, compound and triple. Minutes of the Proceedings. 1888, 93, 128. (21) Varbanov, P. S.; Doyle, S.; Smith, R. Modeling and optimization of utility systems. Chem. Eng. Res. Des. 2004, 82, 561. (22) Mitra, S.; Sun, L.; Grossmann, I. E. Optimal scheduling of industrial combined heat and power plants under time-sensitive electricity prices. Energy 2013, 54, 194. (23) Brook, A.; Kendrick, D.; Meeraus, A. GAMS: A User Guide, release 23.9.1; The Scientific Press, South San Francisco, CA, 2012. (24) Liu, G. P.; Daley, S. Optimal-tuning PID control for industrial systems. Control Eng. Pract. 2001, 9, 1185. (25) Gaing, Z. L. A particle swarm optimization approach for optimum design of PID controller in AVR system. IEEE Trans. Energy Conversion 2004, 19, 384. (26) Krohling, R. A.; Rey, J. P. Design of optimal disturbance rejection PID controllers using genetic algorithm. IEEE Trans. Evol. Comput. 2001, 5, 78. (27) Zamani, M.; Karimi-Ghartemani, M.; Sadati, N.; Parniani, M. Design of a fractional order PID controller for an AVR using particle swarm optimization. Control Eng. Pract. 2009, 17, 1380. (28) Butcher, J. C. Numerical Methods for Ordinary Differential Equations; Wiley: Chichester, 2003. (29) Romagnoli, J. A.; Palazoglu, A. Introduction to Process Control, 2nd ed.; CRC Press: Boca Raton, FL, 2012. (30) Currie, J.; Wilson, D. I. OPTI Toolbox, 2012.

J

DOI: 10.1021/acs.iecr.7b00869 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX