An Approximation Solution to Refinery Crude Oil Scheduling Problem with Demand Uncertainty Using Joint Constrained Programming

This paper is devoted to develop an approximation method for scheduling refinery crude oil operations by taking into consideration the demand uncertainty. In the stochastic model the demand uncertainty is modeled as random variables which follow a joint multivariate distribution with a specific correlation structure. Compared to deterministic models in existing works, the stochastic model can be more practical for optimizing crude oil operations. Using joint chance constraints, the demand uncertainty is treated by specifying proximity level on the satisfaction of product demands. However, the joint chance constraints usually hold strong nonlinearity and consequently, it is still hard to handle it directly. In this paper, an approximation method combines a relax-and-tight technique to approximately transform the joint chance constraints to a serial of parameterized linear constraints so that the complicated problem can be attacked iteratively. The basic idea behind this approach is to approximate, as much as possible, nonlinear constraints by a lot of easily handled linear constraints which will lead to a well balance between the problem complexity and tractability. Case studies are conducted to demonstrate the proposed methods. Results show that the operation cost can be reduced effectively compared with the case without considering the demand correlation.


Introduction
In recent years refineries have to explore all potential costsaving strategies due to intense competition arising from fluctuating product demands and ever-changing crude prices. Scheduling of crude oil operations is a critical task in the overall refinery operations [1][2][3].
Basically, the optimization of crude oil scheduling operations consists of three parts [4]. The first part involves the crude oil unloading, mixing, transferring, and multilevel crude oil inventory control process. The second part deals with fractionation, reaction scheduling, and a variety of intermediate product tanks control. The third part involves the finished product blending and distributing process. In this paper, we focus on the first part, as it is a critical component for refinery scheduling operations.
The crude oil scheduling problem has received considerable attention from researchers and different models have been developed on the basis of deterministic mathematical programming techniques. Lee et al. [5] developed a mixedinteger linear programming (MILP) model to solve a shortterm crude oil scheduling problem, in which the linearity of the bilinear constraints is maintained by replacing bilinear terms with individual component flows, but it can lead to composition discrepancy. To overcome the composition discrepancy problem, Wenkai et al. [6] and Reddy et al. [7] proposed an iterative MIP-NLP model and an iterative discrete-time MIP model, respectively. However the models rely on time discretization representations. Recently, most mathematical models of crude oil scheduling operations put emphasis on continuous-time formulation so as to shorten the gap between theoretical research and real-world operation. Chryssolouris et al. [8] studied the similar problem as Lee et al. [5] and took the temperature cut-points into consideration for each distillation. Jia and Kelly [4] formulated the same problem by a state-task-network (STN) 2 The Scientific World Journal continuous-time representation. Hu and Zhu [9] extended the event-based model of Jia et al., [10,11] to the slot model, by eliminating the redundant event points on others, reducing the size of the model, and hence the solution time of the problem.
Most of the current plant planning and scheduling models are based on deterministic programming. However, due to the volatile raw material prices, fluctuating products demands, and other changing market conditions, many parameters in a planning and scheduling model are usually uncertain. Neiro and Pinto [12] constructed a corporate planning model for multiple refineries using scenario based approach. Neiro and Pinto [12] developed a multiperiod MINLP model to deal with uncertainty in product price and crude price. However the scenario based approach provides no obvious information on the relation between reliability and profitability, which is crucial for decision makers. Several recent papers applied chance constrained programming models to the refinery short-term crude oil scheduling problem [13][14][15][16].
All of the aforementioned models in the domestic refinery optimization are either deterministic [5] or stochastic with independent demand distribution [14,15]. Cao et al. [15] considered the demand correlation for different crude mix in the same time period. However, they did not consider the correlation for the same crude mix in different time periods. In this paper we will consider the crude mix demand correlation in demands not only for different crude mix in the same time period but also for the same crude mix indifferent time periods. These considerations have practical significance. For example, if two crude mixes are predominantly used as raw materials in another process, their demands will be positively correlated in each time period. Alternatively, unusually high demand for a crude mix in one time period more often than not is followed by lower than normal demand in the next period, implying negative correlation. Taking into account such information, whenever it is available, enables a more efficient allocation of crude mix capacity to minimize cost and meet certain marketing objectives.
In this paper, we will propose a stochastic multiperiod model with considering the uncertain crude mix demand correlations. The model employs two-level time structure formulation in which the entire scheduling horizon is divided into several interrelated macroperiods. Each macroperiod with fluctuation demand consists of time intervals with fixed length. A chance constrained programming formulation is developed for solving the problem. The deterministic form of the stochastic constraint is used in solving the problem iteratively. In real-world situations, the future demand always changes as time rolls forward. To deal with the uncertainty of the model, it is important to adjust the planning policy and update the corresponding schedule and the correlation structure of the demand at the end of each time period based on the real sales [17].
The rest of the paper is organized as follows. The scheduling problem is specified in Section 2. In Section 3, a stochastic multiperiod model with two-level time structure formulation is formulated. In Section 4, the deterministic representation of the model is presented. In Section 5, an approximation method combining relax-and-tight technique is developed to solve the joint chance constrained problem above. By relax, we mean that the ΣA approach [18] is used to approximate the original-covariance matrix with a new one of simplified covariance structure. By tight, we mean that the joint chance constrains are transformed into several linear constraints with parameterized dependent which are more stringent than the original constraints. Moreover, an update policy upon the realization of the random demands of crude mix is described. A test problem involving correlated random crude mix demands is solved in Section 6, highlighting various modeling and algorithmic issues. Section 7 summarizes the work and provides some concluding remarks.

Problem Statements and Operation Rules
2.1. Problem Statements. The problem studied involves crude oil unloading process from vessels to storage tanks, transferring process from storage tanks to charging tanks (where several crude oils are mixed) and charging process from charging tanks to crude oil distillations (CDUs). Figure 1 shows the typical processes. During a given scheduling horizon, crude vessels arrive in the vicinity of the refinery docking station and, according to FCFS, wait for unloading of the preceding vessel in the docking station. At the docking station, crude oil is unloaded into storage tanks. Crude oil is then transferred from storage tanks to charging tanks which are buffers to produce a crude mix, of which component compositions were determined at the planning level. The crude oil mix in each charging tank is then charged into a CDU. Given the configuration of the multistage system as well as the uncertain arrival times of vessels, equipment capacity limitations, and key component concentration ranges, the problem is then to determine the following operating variables to minimize operating costs: (a) waiting time of each vessel at sea, (b) unloading time of each vessel, (c) crude unloading rate from vessels to storage tanks, (d) crude oil transfer and mixing rate from storage tank to charging tanks, (e) inventory levels of storage and charging tanks, (f) CDU charging rates, and (g) sequence for charging mixed crude into each CDU.
A basic deterministic model assumption was presented in [5]. However, in this paper we will extend it to a The Scientific World Journal 3 new mixed-integer nonlinear stochastic programming model with chance constraints. The demands for different crude mix are uncertain and possibly correlated, reflecting changing market conditions and periodic variation in customer orders. An optimal planning policy for minimizing expected cost is developed. This is achieved by letting crude mix demands be satisfied with at least a prespecified probability level.
The proposed stochastic model involves the following features and assumptions.
(1) A two-level time structure introduced by Fleischmann and Meyr [19] is adopted for modeling a general stochastic system. The entire planning interval is divided into macroperiods each with fluctuation demand and each macroperiod consists of time intervals with fixed length.
(2) The demands for different crude mixes are modeled as multivariate normally distributed random variables.
The normality assumption has been widely invoked in the literature because it captures the essential features of demand uncertainty and it is convenient to use. The use of more "complex" probability distributions is hindered by the fact that statistical information apart from mean and covariance estimates of product demands is rarely available.
(3) The demand correlation, in the problem formulation, stands for correlation of demands for different crude mixes in the same time period and demands for the same crude mix in different time periods.
(4) The penalty for some crude mix shortfalls is proportional to the amount of underproduction.

Operation Rules.
The operation rules are shown as follows: (1) the refinery uses only one docking station, and a new arriving vessel has to wait at sea until the anterior vessel leaves the docking station; (2) while a charging tank is charging CDU, crude oil from the storage tanks cannot be fed into the charging tank and vice versa; (3) each charging tank can only charge one CDU at each time interval; (4) each CDU can only be charged by one charging tank at each time interval; (5) CDUs must be operated continuously throughout the scheduling time horizon.

Mathematical Model
Indices and Sets  In this paper, a two-level time structure is used to model a general dynamic production system. We consider the crude mix to be scheduled over a finite planning horizon consisting of macroperiods = 1 ⋅ ⋅ ⋅ . Each macroperiod is divided into time intervals with fixed length, where represents the starting time of macroperiod . Let , , a random variable, denote the demand of crude mix in macroperiod . We approximate the crude mix demand as a multivariate normal distribution with a specific correlation structure . It means that the demands for different mix crudes not only are correlated but also are the demands for the same mixed crude in different macroperiods.
The mathematical model is formulated as follows.

(1) Operating Cost
Minimize: Operating cost = unloading cost for the crude vessel + cost for vessel waiting in the sea + inventory cost for storage and charging tanks + changeover cost + underproduction penalty cost Subject to: (2) Vessel Arrival and Departure Operation Rules. Each vessel arrives at the docking station for unloading only once throughout the scheduling horizon: Each vessel leaves the docking station only once throughout the scheduling horizon. Consider Equation for unloading initiation time: Equation for unloading completion time: Each crude vessel should start unloading after arrival time set in the planning level: The Scientific World Journal 5 Duration of the vessel unloading is bounded by the initial volume of oil in the vessel divided by the maximum unloading rate: ⌈⌉ corresponds to round-up of the next highest integer value. Vessel in the sea cannot arrive at the docking station for unloading unless the preceding vessel leaves: Unloading is possible between time V and V : (3) Material Balance Equations for the Vessel. For vessel v, the oil volume at time equals the initial volume minus transferred volume from vessel V to storage tanks up to time t: Operating constraints on crude oil transfer rate from vessel V to storage tank at time t The volume of crude oil transferred from vessel V to storage tanks during the scheduling horizon equals the initial crude oil volume of vessel v: (4) Material Balance Equations for Storage Tanks. The oil volume in storage tank at time equals the initial oil volume in tank plus the oil volume transferred from vessels to tank up to the time and minus oil volume removed from tank to charging tanks up to the time t: Operating constraints on mixed oil transfer rate from charging tank to CDU at time t Volume capacity limitations for charging tank at time t The proposed formulation involves stochastic expressions, which requires a different course of action for transforming it into an equivalent deterministic form. The deterministic equivalent representation of the expectation of the objective function is examined in the next section.

Expectation of Penalty
Cost. The expectation of the underproduction penalty cost term for crude mix in macroperiod is equal to where To facilitate the calculation of the expectation, the standardization of the normally distributed variables , and variables ∑ CDU =1 ∑ +1 −1 BC , , is performed first. Normal random variables can be recast into the standardized normal form, with a mean of zero and a variance of 1, by subtracting their mean and dividing by their standard deviation (square root of variance). This defines the standardized normal variables where, denotes the mean of , and , the square root of its variance. In the same way, the "standardization" of the variables ∑ Using the method in paper [17], the underproduction penalty term , is charging into the following form: where is the standardized normal distribution function and Φ denotes the cumulative probability function of a standard normal random variable.

Satisfaction Level for Single Crude Mix
Demand. The minimization of the objective function (1), as defined above, establishes the production and planned sales policy which most appropriately balances profits with inventory costs and underproduction shortfalls. A crude mix demand satisfaction level is not explicitly specified, but rather it is the outcome of the minimization of the profit function. While higher values of the parameter PEN conceptually increase the probability of demand satisfaction, this strategy may still lead to unacceptably low probabilities of satisfying certain crude mix demands (see examples). Therefore, the setting of explicit probability targets on crude mix demand satisfaction is much more desirable. A systematic way to accomplish this is to impose explicit lower bounds on the probabilities of satisfying a single crude mix demand. This requirement for crude mix in macroperiod assumes the following form: This constraint, known as a chance constraint, imposes a lower bound , on the probability that the crude mix demand realization will be greater than the planned sales , for crude mix in macroperiod . The deterministic equivalent representation can be obtained based on the concepts introduced by Charnes and Cooper [20]. Specifically, by subtracting the mean and dividing by the standard deviation of , , the chance constraint can be written equivalently as The right-hand side of the inequality within the probability sign is a normally distributed random variable with a mean of zero and a variance of 1. This implies that the chance The Scientific World Journal 7 constraint can be replaced with the following deterministic equivalent expression: The application of the inverse of the normal cumulative distribution function −1 , which is a monotonically increasing function, yields Inspection of the deterministic equivalent constraint reveals that it is linear in the deterministic variables and is composed of the mean of the original constraint augmented by the squared root of its variance times −1 ( , ).

Satisfaction Level for All Crude Mix Demand.
In some cases a probability target is desired for the demand satisfaction of a group of crude mix in different macroperiod. For example, when only up to 10%, unsatisfied crude mix demand can be tolerated throughout the entire horizon without distinguishing between different crude mixes. While singleproduct chance constraints are unaffected by correlations between crude mix demands, this is not always the case for joint chance constraints.
This gives rise to joint chance constraints. Joint chance constraints impose a probability target of simultaneously satisfying the demands for a group of crude mixes in different macroperiods as follows: where is the set of product-period ( , ) combinations whose simultaneous demand satisfaction with probability of at least is sought and = {( , ) | 1 ≤ ≤ BT , ≤ ≤ } is the set of all joint chance constraint. The joint probability, , is the probability of intersection of individual constraints to be satisfied.
We use the approach in [21] to solve the joint chance constrained problem. The joint chance constraint equation (17) can be written as separated equivalent deterministic equations. Initially we choose the constraints in a manner such that together they are more stringent than (17). The initial set of individual linear constraints equation (18) replacing (17) was obtained using the following argument. First we denote the event ∑ CDU =1 ∑ +1 −1 , , ≥ , by , and its complement event ∑ CDU =1 ∑ +1 −1 BC , , < , by , . From Boole's inequality of probability theory [22] it is well known that Now, because , is normally distributed with mean , and standard deviation , , where / is the 100(1 − ((1 − )/ ))th percentile of the standard normal distribution. Subsequently we will suppress the subscript for z and replace (20) by the following equation: The initial value of will be chosen to be = / . Using the individual constraints equation (21) results in a more conservative solution than the joint chance constraint equation (17). Therefore, it is required to update the value of z in order to find a solution which satisfies the required confidence level. The method in updating z proposed in [23] is used, which is based on the interpolation of z values. The deterministic forms of the objective function and stochastic constraint are used in solving the joint chance constrained problem iteratively by using a different z value. The algorithm is examined in the next section.

The Algorithm and Update Policy for the Joint Chance Constrained Problem
In this section, a method is developed to solve the joint chance constrained problem above, and then an update policy upon the realization of the random demands of crude mix is described.

Algorithm for the Joint Chance Constrained Problem.
The following three steps are performed iteratively.
Step 1 (initialization). Problem is solved with the equivalent deterministic objective function and the constraint equation (21) is replaced with (17) for a fixed scalar z to determine the set of new controlled variables. The algorithm starts by choosing a high value for the initial -value as in (21), which makes the corresponding solution satisfy the demand with a probability level higher than P target = . 8 The Scientific World Journal Step 2 (MINLP solution). The problem is implemented with GAMS (General Algebraic Modeling System). The CPLEX is used to solve this deterministic problem.
Step 3 (evaluate probability). Evaluate the multivariate normal probability (17). A subregion adaptive algorithm proposed by Genz and Bretz [24] is employed to carry out multivariate integration which makes this evaluation feasible.
Step 4 (updatez-value). If the evaluated probability level is in the neighborhood of P target , the algorithm terminates since the goal of finding a schedule that satisfies the demand with a probability of P target is accomplished; otherwise the value is updated and the previous steps are repeated to obtain another schedule.
To update the z-value the following algorithm is used. The goal is to find a z-value in (21) that provides a schedule such that the demand can be satisfied with a probability of P target = over the entire time horizon. The z-value needs to be obtained iteratively.
Step 1 (obtain the upper and lower bounds of z). First we choose z = z in (21). Obviously it yields a lower bound for the correct z-value. We call it z lower . We now run Step 2 of the algorithm for this lower bound and obtain an estimate of the probability with which the demand is being met. We call this probability P lower . Next we choose an arbitrary large value for z. We denote it by z upper and obtain a similar estimate of the probability of demand satisfaction. We denote this probability by P upper .
Step 2 (obtain the upper and lower bounds of P target ). In this step we obtain the upper percentiles of the univariate standard normal distribution for these probabilities and denote P upper and P lower by z 2 and z 1 , respectively. We also denote the corresponding percentile for the P target value by z target .
Step 3 (update z-value). Based on these values the updated z-value is obtained using the following linear interpolation formula: If the Z new value is lower than Z 2 and higher than Z target , then we replace Z 2 byZ new . If it is lower thanZ target and higher than Z 1 , we replace Z 1 by Z new . We repeat this process using (17) until P target is reached.

Update of the Planning and Scheduling Policy.
As discussed in Section 3, the demand ∼ ( , Σ) follows a multivariate normal distribution with a specific correlation structure. In order to decrease the uncertainty associated with the random variables , and thus increase the predictive power of the stochastic model, it is proposed to update the corresponding schedule at the end of each macroperiod based on the realized production [17]. The update of the schedule is accomplished by recalculating the conditional multivariate probability function given the demand realizations in previous periods.
After solving the scheduling problem, only the decision variables associated with the immediately following period may be taken as final. The decision variables referring to successive periods can be used for planning and activities related to the operation.
The update of the planning policy and schedule at the end of macroperiod is described in the following steps.
Step 1 (update statistics). Update future demand statistics based on present and previous product demand realizations.
The first step allows the use of updated demand forecasting information. If there are no new demand forecasts available, the conditional multivariate probability density function based on the realizations of the demands in previous periods can be utilized. The procedure for finding the conditional distribution function is outlined below.
Step 2 involves the solution of the problem, which, in turn, determines the new planning and scheduling policies.
The derivation of the conditional probability distribution based on the realization of previous random variables can be accomplished as follows. First, the random variable , is partitioned into two sets. The first set contains the uncertain crude mix demands for the future periods. Based on this partitioning, the variance-covariance matrix can be expressed as follows: where Σ PP and Σ FF are the variance-covariance submatrices of the crude mix demands belonging to the sets and F, respectively. The elements of submatrices Σ PF = (Σ FP ) are the covariances between the elements of and . Let denote the realizations (outcomes) of the crude mix demands of set associated with past periods. The conditional means of the uncertain future crude mix demands include where / denotes the conditional demand expectations and and are the past and future mean values. The new (conditional) variance-covariance matrix is equal to The Scientific World Journal 9 A detailed treatise of conditional multivariate normal probability functions can be found in Tong, 1990 [18]. Note that while the conditional mean value / depends on the demand realizations in past periods, the conditional covariance matrix Σ / is independent of the demand realizations Ξ. Updating the probability distribution at the end of each period decreases the uncertainty (variances) associated with the remaining random variables, thus increasing the predictive power of the stochastic model. This is observed in the case study in Section 6.

Case Study
In the case study, three different crude mixes are to be produced. The detailed data of the case study is given in Table 1. The time horizon of 30 days is divided into 30 equal times. The schedule involves three macroperiods with 10 days. The product demands in each macroperiod are described by normal multivariate probability distributions. The expected (mean) values of the demands are given in Table 2, and their standard deviation is assumed to be 25% of their mean values.  The alternative model formulation is defined and solved using the proposed solution procedure. Mode: ≥ others (1-7) .
Model incorporates not only single but also joint chance constraints for setting probability targets for the satisfaction of a group of products. First the effect of correlation on the economic parameters and the optimal product mix is evaluated and discussed, and next the proposed update policy is applied on the example problem.
The problem is implemented with GAMS. The CONOPT and CPLEX 4.0 solvers are used to solve the formulated models above.
6.1. Use of Joint Probability Constraints. Two separate cases are considered for model involving (i) uncorrelated and (ii) arbitrarily correlated crude mix demands. For the correlated crude mix demands, a variance-covariance matrix is constructed and shown in (30). The matrix has sparsity of 60%, and its off-diagonal elements vary between −0.3 and 0.6. The bias toward positive covariance is introduced to maintain semipositive definiteness which is a property of all variancecovariance matrices. Consider ) .

(30)
The resulting MILP problem involves the objective function and a set of linear constraints representing the deterministic equivalent representation of the individual and joint probability constraints. The joint probability index , for the set of all 9 crude oil demands, is set to 50%. For the value of is identified, the expectation of the cost is compared with those without considering demand correlation. Figure 2 shows the minimization expected cost for the two cases. These results indicate that the expectation of the total cost in the model considering the correlation decreases.
The lower expectation of the total cost results from the fact that the joint probability target can be met with smaller production levels than the ones ignoring correlations.

Revision of Planning Policy.
Using the correlation matrix in (30), a demand realization is randomly generated using the method presented by Tong (1980) [18]. The crude mix demands realizations are given in Table 3; a joint probability target of 50% is set to be the demand satisfaction of each crude mix in all time periods. This introduces joint probability constraints in the following form: The resulting optimization problem is solved in the entire time horizon. Next, based on (i) the demand realizations in Table 6: Updates of the product demand correlation matrix.  the first period, (ii) the correlation matrix, (iii) the mean values of the remaining random demands, and (iv) the remaining inventory at the beginning of the second period, a new planning revision problem is solved for macroperiod ∈ (2-3). This is repeated for all remaining macroperiods. The updated values of the mean vectors, standard deviations, and covariance matrix after each time period are given in Tables  4, 5, and 6.
Two separate cases are considered for the models involving (i) correlated crude mix demands using the proposed policy and (ii) correlated crude mix demands without the policy. The actual production amount for crude mix 1 is The realized demand Policy with updating Policy without updating shown in Figure 3. Examination of the data reveals that the actual production amount using the updating policy is close to the actual realized demand.
Clearly, when the solution is updated, the planning policy is much more effective because overproduction is decreased. The proposed scheduling updating policy by updating the probability distribution at the end of each period decreased uncertainty associated with the remaining random variables and thus increased the predictive power of the stochastic model for the subsequent periods.

Conclusion
In this paper, the multiperiod crude oil scheduling problem under demand uncertainty was addressed. A stochastic two-level time structure model considering the uncertain crude mix demand correlations is formulated. The model involves the minimization of the expected cost subject to constraints for the satisfaction of single-and/or multiplecrude mix demands with a prespecified level of probability. A chance constrained programming formulation is developed for solving the problem, and the deterministic form of the stochastic constraint is used to solve the problem in an iterative way. Moreover, a revision strategy for the update of the corresponding schedule at the end of each time period based on the realized production was also presented. This was accomplished by recalculating the conditional multivariate probability function given the demand realizations in previous periods.
Results on the correlated case study demonstrate that the expected profit and in particular the corresponding planning policy and schedule are strongly affected by the presence of correlations. The proposed planning revision and scheduling updating policy by updating the probability distribution at the end of each period decreased uncertainty associated with the remaining random variables and thus increased the predictive power of the stochastic model for the subsequent periods.