Resiliency Issues in Integration of Scheduling and Control

222

Ind. Eng. Chem. Res. 2010, 49, 222–235

Resiliency Issues in Integration of Scheduling and Control Kishalay Mitra,† Ravindra D. Gudi,*,‡ Sachin C. Patwardhan,‡ and Gautam Sardar† Tata Consultancy SerVices, 1 Mangaldas Road, Pune 411001, India, and Indian Institute of Technology, Bombay, Powai, Mumbai 400076, India

Integration of scheduling and control in process manufacturing systems has traditionally resulted in greater coordination between the different layers of hierarchy in optimization and control. However, the effectiveness of this coordination as well as the overall responsiveness of the manufacturing system could be greatly improved if such integration explicitly accounts for parametric variations at each of the layers. This paper proposes an approach to the problem of robust integration of the scheduling and control layers in an uncertainty analysis framework so as to yield robust manufacturing systems that can perform well in the presence of the parametric variations. The uncertainty analysis is performed here using different uncertainty approaches such as the chance constrained programming, fuzzy mathematical programming, and robust optimization. In this paper, we propose to analyze the impact of the uncertainty on the different manufacturing objectives in a multiobjective Pareto sense. The deterministic integrated scheduling and control model taken from published work [FloresTlacuahuac, A.; Grossmann, I. E. Ind. Eng. Chem. Res. 2006, 45, 6698] forms the basis of our work toward addressing the uncertainty issues. We demonstrate the development of improved and robust integration of the scheduling and control tasks using the proposed frameworks. Introduction Integration of scheduling and control in process manufacturing systems has traditionally resulted in greater coordination between the different layers of hierarchy in optimization and control, though the solution of this monolith configuration is not trivial. To start with, researchers dealt with this complicated integrated scheduling and control (ISC) problem by decomposing it into individual scheduling and control tasks, even though their real life realizations might be simultaneous. In approaching the optimal control problems alone, the purpose was to attain optimal transition trajectories between different sets of product sequences while minimizing the transition time, assuming that the associated scheduling decisions are already known and fixed.1,2 On the other hand, when the scheduling problem is solved alone, the target was to obtain the optimal settings for different production parameters that optimize an objective such as profit or makespan, for the already determined optimal control assignments.3 The joint consideration of scheduling and control enables the fulfillment of opportunities where incorporation of better control could result in improved schedules of manufacturing. Likewise, a proper incorporation and solution of scheduling problems enables realistic trajectory targets for the control problem. As a result of not considering the scheduling and control together, suboptimal solutions as compared to the integrated solution could emerge.4-7 Recognizing this fact, the simultaneous solutions of design and control,8-10 scheduling and design,4,11,12and scheduling and control7,13-17 problems have been the subject of various research work in the literature, where the research focus has adopted an integrated approach. Though the work on integration between the different layers of hierarchy in optimization and control has reached significant maturity, the effectiveness of this coordination as well as the overall responsiveness of the manufacturing system could be further enhanced if such integration explicitly accounts for variations in parameters at each of these layers, since in practical * To whom correspondence should be addressed. Phone: 91-2225767204. Fax: 91-22-2572-6895. E-mail: [email protected]. † Tata Consultancy Services. ‡ Indian Institute of Technology.

applications such parametric variations are quite common. Apart from the other sources of uncertainties, the significant sources could arise from model plant mismatch and also due to inherent variations in the manufacturing and environmental data. There has been significant progress in the broad area of optimization under uncertainty or robust optimization following the seminal work of Bellman and Zadeh,18 Charnes and Cooper,19 Dantzig,20 and Tintner,21 that aim at developing approaches and methodologies to create reliable solutions, which can remain feasible in the presence of parameter uncertainty.22,23 Some of the preventive uncertainty handling approaches in the literature that have been developed with strong theoretical backgrounds are two stage stochastic programming (TSSP), chance constrained programming (CCP), fuzzy mathematical programming (FMP), and robust optimization (RO).22 In the standard two stage approach,24-32 the decision variables are intelligently partitioned into two sets. The first stage variables are generally independent of the uncertain variables and to be decided before the realization of uncertain parameters (“here and now” decisions). The second stage variables (“wait and see” decisions) are associated with the uncertain parameters and can assume different sets of values for corresponding different sets of realization of the uncertain parameters, for a single set of first stage variables. The objective here is to minimize the sum of the first stage costs and the expected value of the second stage costs by determining the first stage decision variables. Uncertainty or parametric variation in the two stage approach has been represented in the literature in two different ways. The scenario-based approach29,30,33 represents an uncertain parameter by forecasting all its future outcomes. Different probability measures are associated with these outcomes, which are used for calculating the expected second stage costs. This method has a major drawback that it leads to a combinatorial explosion in problem size with an increase in the number of uncertain parameters. On the other hand, the distribution-based approach31,32,34,35 overcomes this difficulty to some extent by assuming continuous probability distribution for random parameters. In this case, the problem has been shown to become nonlinear and would involve integration of the multivariate

10.1021/ie900380s CCC: $40.75  2010 American Chemical Society Published on Web 11/19/2009

Ind. Eng. Chem. Res., Vol. 49, No. 1, 2010

probability integrals, which is primarily resolved through the discretization of the probability space using techniques such as Monte Carlo, Gaussian quadrature, Latin hypercube, shifted Hammerslay sequence sampling,36 etc. The disadvantage of these methods is the sharp increase in computational requirements with the increase in number of uncertain parameters. Moreover, proper cost assignment to recourse activities in the aforementioned approach often either leads to actions that are difficult to carry out practically or are very hard to specify. An alternate way to deal with this problem is to measure the resiliency of the system to parametric uncertainties by proposition of reliability metric that is expressed in terms of a decision to be feasible within a user specified probability. In CCP, this is performed by associating a probability measure with constraints having uncertain parameters, which, after mathematical treatment, leads to a relatively smaller, equivalent deterministic formulation and can be solved relatively easily.37-41 Likewise other alternative approaches to represent uncertainty are FMP and RO. In FMP,42,43 the constraints having uncertain parameters are assumed to have fuzzy boundaries that can be expressed in terms of membership functions and relaxed to the extent of a decision maker’s expectations. An optimization model that maximizes the smallest membership function (max-min formulation) of all constraints provides the solution. Application of FMP is widely spread across different applications such as capacity planning,44 supply chain planning,45 production scheduling,46 bioenergy production,47 etc. to name a few. On the other hand, in RO based techniques, the idea is to find a solution which is either solution robust, if it remains close to the optimal for all scenarios, or model robust, if it remains feasible for most scenarios.48 Followed by some pioneering work49-51 and subsequent algorithmic developments,52,53 robust optimization has been applied to applications such as machine scheduling,54 logistics,55,56 and scheduling.48,52,53,57 The progress of theoretical research on convexity properties for the deterministic equivalent in each of these cases is primarily limited to various facets of linear systems and only some specific results are available for complicated nonlinear systems.39 Applications of CCP, FMP, and RO in the process systems engineering (PSE) literature have been relatively fewer. In this paper, we propose different approaches to the problem of robust integration of the scheduling and control layers in an uncertainty analysis framework so as to yield robust manufacturing systems that can perform well in the presence of the parametric variations. Considering the problem of computational intractability of TSSP with the increase in the number of uncertain parameters, the uncertainty analysis is performed here using other preventive uncertainty handling approaches such as the CCP, FMP, and RO. We propose the use of these uncertainty handling approaches for the task of integration of scheduling and control and analyze the impact of the uncertainty on the different manufacturing objectives. The objective of profit maximization has a trade-off with reliability of the solution obtained by the aforementioned approaches. This frames an ideal platform to analyze the problem in a multiobjective Pareto sense that unfolds the mentioned trade-off and helps a user to decide the single point of operation considering the benefits/sacrifices associated with each of these Pareto optimal solutions. The deterministic model on integrated scheduling and control, taken from the work of Flores and Grossmann,15 forms the basis of our work toward addressing the uncertainty issues. We demonstrate the development of improved and robust integration of the scheduling and control tasks using the proposed frameworks.

223

The rest of the paper is organized as follows. In the next section, a brief overview on the deterministic ISC model of Flores-Tlacuahuac and Grossmann15 is presented. Adaptation of the deterministic ISC problem to various preventive uncertainty handling paradigms is presented next. A Pareto evaluation of the objectives of the ISC problem under demand uncertainty are presented in the results and discussion section for the first case study of the work of Flores-Tlacuahuac and Grossmann,15 followed by summary and concluding remarks. Formulation Deterministic Integrated Scheduling and Control Problem. While considering the case of integrated scheduling and control, we intend to consider both the sequencing and optimal control aspects simultaneously. For this, a problem is considered, where given the number of products to be manufactured in a continuously stirred tank reactor (CSTR), lower bounds on demand rates for each of them, steady state operating conditions, and price of each of the products, the purpose of the deterministic model is to find the sequence in which the products should be manufactured, total cycle time, transition times, production rates, length of processing times, and amounts manufactured for each product and manipulated variables for the transition so that the profit is maximized. The generic deterministic mixed integer dynamic optimization (DMIDO) model of Flores and Grossmann15 is considered here and presented briefly as follows: DMIDO. The objective function of the deterministic model defines the profit function as the difference between the revenue and different cost components (inventory and transition costs). max profit ) φ1 - φ2 - φ3 Np Np Cpi Wi Csi (Gi - Wi/Tc) , φ2 ) , where φ1 ) Tc 2Θi i)1 i)1 Ns Nfe Ncp Crt fe,co,kΩc,Ncp φ3 ) hf,k × Tc k)1 fe)1 co)1

∑

(

∑∑

∑

∑

1 n (xfe,co,k - jx1k )2 + ...(xfe,co,k - jxnk )2+

1 m - uj1k )2 + ...(ufe,co,k - ujmk )2 (ufe,co,k

(1)

)

The aspect that (a) any product can only be manufactured once within each production wheel and (b) only one product is manufactured at each time slot are expressed in the form of product assignment by eqs 2 and 3, respectively. Equation 4 defines a backward binary variable (y′ik) denoting that such a variable for product i in slot k takes the value assigned to the same binary variable but one slot backward, k - 1. Naturally, at the first slot, this variable is the value of the same variable at the last slot as defined by eq 5 (cyclic production wheel assumption). Ns

∑y

i,k

) 1 ∀i

(2)

i,k

) 1 ∀k

(3)

k)1 Np

∑y k)1

yi,k ′ ) yi,k-1 ∀i,k * 1

(4)

yi,1 ′ ) yi,N ′ s ∀i

(5)

Equation 6 states that the total amount manufactured of each product i must be equal to or greater than the duration of the

224


production wheel times desired demand rate. The amount manufactured of product i is computed as the product of the time used (Θi) for manufacturing such product and the production rate (Gi) (eq 7). The production rate is computed from the feed streamflow rate (F0) and the conversion (Xi) (eq 8). Wi g DiTc ∀i

(6)

Wi ) GiΘi ∀i

(7)

nonlinear equations that must be satisfied. The constraint given by eq 19 is used to calculate the value of the system states at each of the discretized points (xnf,c,k) by using the monomial basis representation. The length of all finite elements is the same and computed as hf,k ) 1/Nfe. Equation 20 is used to force continuous state profiles on all the elements at the beginning of each element n ). The value of the first-order derivatives of the systems at (x0,f,k every discretized point is computed by eq 21. Ncp

Gi ) F (1 - Xi) ∀i 0

(8)

n xfe,co,k

)

n x0,fe,k

+

θtkhfe,k

∑Ω

∀n,fe,co,k;hf,k ) 1/Nfe

n ˙ fe,l,k l,cox

i)1

The processing time constraints are given by eqs 9 (an upper bound on the time used for manufacturing product i at slot k), 10 (time used for manufacturing product i), and 11 (duration time at slot k) whereas eq 12 is used for defining the binary production transition variable zi,p,k. If such a variable is equal to 1, then a dynamic transition occurs from product i to product p within slot k; otherwise zi,p,k is zero. θi,k e qmaxyik ∀i, k

(9)

Na

Θi )

∑θ

i,k

∀i

(10)

i,k

∀k

(11)

k)1 Np

∑θ

pk )

i)1

zi,p,k g yp,k ′ + yi,k - 1 ∀i,p,k

(12)

Equations 13-18 are used to express the timing relationships. Equation 13 defines the transition time from product i to product t p at slot k. Here the term tp,i needs to be computed iteratively. Equation 14 sets time to zero at the beginning of the production wheel cycle for the first slot. Equation 15 calculates the slot end time as the sum of the slot start time plus the processing time and the transition time. Equation 16 states that the start time at all the slots is just the end time of the previous slot except the first one. Equation 17 forces the end time of each of the slots to be less than the production wheel cyclic time. Equation 18 is used to obtain the time value inside each finite element and for each internal collocation point. Np

θtk

)

∑∑

∀k

Ncp

n n ) x0,fe-1,k + θtkhfe-1,k x0,fe,k

)0

) tek

e tk-1

n 1 n 1 m x˙fe,co,k ) fn(xfe,co,k ,...,xfe,co,k , ufe,co,k , ...,ufe,co,k ) ∀n,fe,co,k (21)

Fixing the initial and final values for the dynamic model equations is done through eqs 22-30. The desired value of each n ) is derived by eqs 22 and state at the beginning of slot k (xin,k 23. Equation 24 defines the values of the state variables at the end of each slot k (xjnk ). Similarly, the values of the manipulated m and ujmk ) variables at the beginning and end of each slot k (uin,k are defined by eqs 25 and 26, respectively. Equation 27 enforces the system states to take the desired state values at each slot k. Equation 28 fixes the values at the first finite element and first m ). Equation 29 determines collocation point of each slot k (u1,1,k the values of the manipulated variables at the last finite element and last collocation point of slot k (uNmfe,Ncp,k), and eq 30 determines the values of the system states at the beginning of n ). each slot (x0,1,k Np

∑x

xnin,1 )

∀n

(22)

∀n,k * 1

(23)

∀n,k

(24)

n ss,iyi,Ns

i)1

Np

xnin,k )

∑x

n ss,iyi,k-1

i)1

Np

∑x

n ss,iyi,k

i)1

(13)

∀k * 1

(14) (15) (16)

umin,1 )

∑u

∀m

(25)

∀m,k * 1

(26)

m ss,iyi,Ns

i)1

Np

xmin,k )

∑x

m ss,iyi,k-1

i)1

Np

e Tc ∀k

(17)

θtk θtk + γ ∀fe,co,k Nfe Nfe c

(18)

tfe,co,k ) (f - 1)

∀n,fe g 2,k

Np

tek ) tsk + pk + θtk ∀k tsk

n ˙ fe-1,l,k l,Ncpx

(20)

i)1 p)1

ts1

∑Ω i)1

jxnk )

Np

ttp,izi,p,k

(19)

ujmk )

∑u

m ss,iyi,k

∀m,k

(27)

i)1

The dynamic model is discretized using the method of orthogonal collocation on finite elements. In this procedure, a given slot k is divided into a number of finite elements and an adequate number of internal collocation points is selected within each of these finite elements. The set of ordinary differential equations is approximated at each collocation point, leading to a set of

um1,1,k ) umin,k ∀m,k

(28)

uNmfe,Ncp,k ) u˜min,k ∀m,k

(29)

xn0,1,k ) xnin,k ∀n,k

(30)

Finally, eqs 31 and 32 constrain the values of both the system states and manipulated variables, respectively, to lie within acceptable lower and upper bounds.


∀n,f,c,k

(31)

ummin e umfck e ummax ∀m,f,c,k

(32)

xnmin

e

xnfck

e

xnmax

The following notations used in the above set of equations are n m , uf,c,k explained here for further clarity. xf,c,k , jxnk , ujmk , tf,c,k, tek, tsk represent the nth state in finite element f, collocation point c, slot k, mth manipulated variable in finite element f, collocation point c, slot k, the state desired value at the end of slot k, the desired value of the manipulated variable at the end of slot k, the time value inside finite element f, collocation point c, slot k, the final time at slot k, and the starting time at a slot k, respectively. Here, yi,k represents the binary variable to denote ′ is an auxiliary binary if a product i is assigned to slot k, yi,k variable, and zi,p,k denotes if product i is followed by product p in slot k. A few other important decision variables are processing time at slot k (pk), production rate (Gi), total production wheel time (Tc), amount produced for each product (Wi), processing time of product i at slot k (θik), transition time at slot k (θtk), total processing time of product i (Θi), and conversion (Xi). Some other parameters are number of products (Np), number of slots (Ns), number of finite elements (Nfe), number of collocation points (Ncp), number of system states (Nx), number of manipulated variables (Nu), demand rate (Di), length of finite element (hf,k), matrix of Radau quadrature weights (ΩNcp,Ncp), upper bound on processing time (θ max), nth steady state value n ), mth manipulated variable value of product i of product i (xss,i m (uss,i), feed stream volumetric flow rate (F0), and roots of Lagrange orthogonal polynomial (γNcp), whereas various unit cost parameters are inventory holding cost (Csi ), revenue (Cpi ), and raw material price (Cr). Additionally, the indices i, k, fe, co (or l), n, and m are used for denoting products, slots, finite elements, collocation points, system states, and manipulated variables, respectively. From the formulation, it is clear that the above ISC formulation is mixed integer nonlinear programming (MINLP) in nature. Interested readers are referred to the original work15 for details on the explanations of model equations. Assuming the uncertainty in demand, the way of adopting the aforementioned DMIDO model, using various preventive uncertainty handling frameworks, e.g. CCP, FMP, and RO, is shown in the following few sections. Integrated Scheduling and Control under Uncertainty CCP Approach. In CCP,37,38 constraints that are associated with random parameters, are ensured of getting satisfied with some user defined probability. In this framework, a standard optimization formulation with uncertainty parameter ξ min{f(x)|h(x,ξ) g 0}

(33)

can be expressed as min{f(x)|P(hj(x,ξ) g 0) g p} (j ) 1,...,u))

(34)

where f(x) is the objective function, x is the set of decision variables, ξ is the set of random parameters, P is the probability measure of the given probability space of the uncertain parameter, and p ∈ (0,1] is the probability level with which each of the u constraints hj(x, ξ) g 0 of the entire constraint set h(x, ξ) g 0 needs to be satisfied. The higher the value of p, the more reliable is the system with respect to the uncertain parameter. On the other hand, the set of feasible x is more and more shrunk with value of p close to unity. On the basis of the requirements of several constraints getting satisfied either

225

individually or together, the methodology is called individual or joint chance constrained programming, respectively. These two different concepts can be represented as eqs 34 and 35, respectively. min{f(x)|P[(hj(x,ξ) g 0) (j ) 1,...,u)] g p}

(35)

Assuming (i) a normal distribution for the random parameter, ξ, and (ii) that random parameters are separable from the decision variables, the constraints in eq 34 can be transformed into an equivalent deterministic form by the following coordinate transformation, where we correct the nominal requirement, ξjj, by quantile value multiplied by standard deviation of the uncertain variable and provide robustness of the generated optimal operating conditions under uncertain situations (see eq 36): P(hj(x,ξ) g 0) g p S P(h˜j(x) g ξj) g p Sh˜j(x) g ξˆ j: ) ξ¯ j + qpσξj

(36)

where ξj is the random parameter associated with the jth constraint, ξjj and σξj are the mean and standard deviation values of the corresponding random parameters, respectively, and qp is the p-quantile of the standard normal distribution with zero mean and unit standard deviation (e.g., 97% probability corresponds to qp ) q0.97 ) 2.0). A value of quantile closer to 2 represents model reliability of 97%; this would imply that an uncertain parameter would lie within the second quantile 97% of the time and, therefore, inclusion of such a description in the constraint satisfaction would represent the uncertainty to the above extent. In this case, the uncertainty terms are not separable from the decision variables, based on whether the uncertain terms are independent or dependent on each other, the corresponding mean and variance terms are incorporated into the equation in the same fashion and the deterministic equivalent form is derived. This deterministic equivalent generally becomes nonlinear as the computation of mean and variance of coefficients of decision variable leads to nonlinearity in terms of decision variables. However, the results for derivation of deterministic equivalent in CCP formulations are only known for limited cases, e.g. (i) uncertain variables are separable from the decision variables and (ii) constraints having uncertain variables are linear. The methods of solving chance constrained optimization problems for nonlinear processes under uncertainty have not been well developed. This limitation restricts us not to consider any uncertain parameter that has a nonlinear relationship with constraints in this study, e.g. uncertainty in a kinetic parameter. As uncertainty is considered in the demand, the demand equation (eq 6) of the DMIDO model is modified assuming a normal distribution for the uncertain parameter. The deterministic equivalent of the demand equation using the CCP approach (considering the right-hand side of the normal curve) is added to it, replacing the original demand equation. This enables us to simulate the uncertain situation under growing demand. The other uncertain situation of diminishing demands can also be simulated considering the left side of the normal curve simply by changing the sign of coordinate transformation. But we restrict ourselves to the right side of the normal curve. Upon modifying the equations, the following multiobjective uncertain model using CCP (MIDOCCP) formulation is derived.

226


max profit max qD j i + qDσD)Tc g 0 ∀i Wi - (D eqs 1-5, 7-32

(37)

As the reliability (qD) of the model has an inherent Pareto tradeoff with the profit of the model, we finally obtain the biobjective optimization formulation for the ISC model. This multiobjective optimization formulation is solved by the ε-constraint approach.59 FMP Approach A brief overview of FMP20,42,43 is presented in the following section. Assuming uncertain parameters as fuzzy numbers and the constraints associated with the uncertain parameters as fuzzy sets, a membership function is defined for each of the uncertain constraints in FMP, where the value of this membership function signifies the extent of constraint violation. For example, while solving the following problem, minimize f(x) such that gj(x) e b,

j ) 1, 2, ..., m

(38)

If we assume a jth linear constraint gj(x) e b in terms of decisions x and an uncertain parameter b that can vary in the range b to (b + ∆b), the membership function, µgj(x)(b), of this constraint, can be defined as follows:

{

if gj(x)bj+∆bj 1

(39)

This linearity assumption is taken to facilitate the linear propagation of the effect of the uncertain parameter with respect to the constraints. In the above constraint, it may often be desirable to achieve constraint satisfaction for higher values of b; a practical analogy to this requirement can be taken in terms of meeting more demand under uncertain conditions, where b is the crisply known value of the demand under certainly known conditions and the same demand b is allowed to vary over (b + ∆b) under uncertain conditions. FMP does not need the uncertain parameter to follow some distribution with known central tendency properties; it only needs the extent of margin that can be allowed on the uncertain parameters to be specified. The nature of the membership function can be of varied nature, however, has been assumed as linear here. In addition to defining membership function for constraints, similar membership functions are also formulated for objective functions with the lower and upper bounds on the objective function, specifying the expectation of the decision maker. For the above example, the objective function bounds can be calculated by solving two instances of the same optimization problem, the first with bj on the right-hand side and the other with bj + ∆bj. We denote these lower and upper bound values of the objective function as f max and f min, respectively. The membership function for the objective function in this case can be defined as

{

1

µf(x) ) 1 0

if f(x) > f min f(x) - f min if f max e f(x) e f min f max - f min if f(x) < f max

(40)

Similarly, the uncertainty in decision variable coefficients can also be captured by allowing the constraint to be flexible in a (non)linear fashion (as shown in eqs 39 and 40) between the two bounds, e.g. one bound when no uncertainty is considered and the other bound when uncertainty is considered in the decision variable coefficients. The final step in FMP approach is to construct an optimization model that attempts to balance the conflicting requirements of minimizing the objective function, while maximizing the satisfactions of the constraints. Bellman and Zadeh18 suggested that this can be done by maximizing the smallest membership function of all constraints, including the same for objective function. Using this max-min criterion with the above membership functions, the optimization problem under uncertainty can, therefore, be transformed in the following manner: max λ

0

Resiliency Issues in Integration of Scheduling and Control

Recommend Documents