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

Authors: Qianqian Duan; Genke Yang; Guanglin Xu; Changchun Pan
Journal: The Scientific World Journal
 (2014)

---

## Abstract

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.


---

## Body

## 1. Introduction

In recent years refineries have to explore all potential cost-saving 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–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 mixed-integer linear programming (MILP) model to solve a short-term 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) 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–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 Section2. 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.


## 2. 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). Figure1 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.Figure 1
Typical flow process for the refinery crude oil operations.A basic deterministic model assumption was presented in [5]. However, in this paper we will extend it to a 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.


### 2.2. 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.


## 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). Figure1 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.Figure 1
Typical flow process for the refinery crude oil operations.A basic deterministic model assumption was presented in [5]. However, in this paper we will extend it to a 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.


## 2.2. 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.


## 3. Mathematical Model

Indices and Sets
i
∈
{
1
,
…
,
N
S
}= crude oil storage tank and the crude oil in it.
j
∈
{
1
,
…
,
N
B
}= crude oil charging tank and the crude oil mix in the charging tank.
k
∈
{
1
,
…
,
N
C
}= key component of crude oil.
l
∈
{
1
,
…
,
N
CDU
}= crude distillation unit.
m
∈
{
1
,
…
,
N
M
} = macroperiod.
t
∈
{
1
,
…
,
N
SCH
}= time interval, N
SCH, denotes scheduling horizon.
T
m= the starting time of the macroperiodm.
v
∈
{
1
,
…
,
N
V
}= crude vessels.
Variables
D
j
,
l
,
t= 0-1 variable to denote if the crude oil mix in charging tank j charges CDU l at timet.
X
v
,
t
F= 0-1 variable to denote if vessel v starts unloading at timet.
X
v
,
t
L= 0-1 variable to denote if vessel v just completes unloading at timet.
X
v
,
t
W = 0-1 continuous variable to denote if vesselv is unloading its crude oil at timet.
Z
j
,
j
′
,
l
,
t = 0-1 continuous variable to denote if transition from crude mix (or charging tank)j toj
′at timet in CDUl.
f
i
,
j
,
k
,
t
SB= volumetric flow rate of component k from storage tank i to charging tank j at timet.
f
j
,
l
,
k
,
t
BC= volumetric flow rate of component k from charging tank j to CDU l at timet.
F
v
,
i
,
t
VS= volumetric flow rate of crude oil from vessel v to storage tank i at timet.
F
i
,
j
,
t
SB= volumetric flow rate of crude oil from storage tank i to charging tank j at timet.
F
j
,
l
,
t
BC= volumetric flow rate of crude oil mix from charging tank j to CDU l at timet.
T
v
F = vesselv unloading initiation time.
T
v
L = vesselv unloading completion and departure time.
v
j
,
k
,
t
B = volume of componentk in charging tankj at timet.
V
v
,
t
V = volume of crude oil in crude vesselv at timet.
V
i
,
t
S = volume of crude oil in storage tanki at timet.
V
j
,
t
B= volume of mixed oil in charging tankj at timet.
Parameters
C
V
UNLOAD= unloading cost of vessel v per time interval.
C
V
SEA= sea waiting cost of vessel v per time interval.
C
i
INVST= inventory cost of storage tank i per time per volume.
C
j
INVBL
= inventory cost of charging tank j per time per volume.
C
j
,
j
′
SETUP = changeover cost for transition from crude mixj toj
′ in CDU.
C
j
PEN = penalty cost for crude mix j shortfalls.
D
j
,
m
M= stochastic demand of crude mix j by CDUs during the macroperiodm.
F
v
,
i
VS
,
min
⁡= minimum crude oil transfer rate from vessel v to storage tanki.
F
v
,
i
VS
,
max
⁡= maximum crude oil transfer rate from vessel v to storage tanki.
F
i
,
j
SB
,
min
⁡= minimum crude oil transfer rate from storage tank i to charging tankj.
F
i
,
j
SB
,
max
⁡= maximum crude oil transfer rate from storage tank i to charging tankj.
F
j
,
l
BC
,
min
⁡= minimum crude oil charging rate from charging tank j to CDUl.
F
j
,
l
BC
,
max
⁡= maximum crude oil charging rate from charging tank j to CDUl.
T
v
ARR= crude vessel v arrival time around the docking station.
V
v
,
0
V= initial volume of crude oil in crude vesselv.
V
i
S
,
min
⁡ = minimum crude oil volume of storage tanki.
V
i
S
,
max
⁡= maximum crude oil volume of storage tanki.
V
i
,
0
S= initial crude oil volume of storage tanki.
V
j
B
,
min
⁡ = minimum mixed crude oil volume of charging tankj.
V
j
B
,
max
⁡ = maximum mixed crude oil volume of charging tankj.
ε
i
,
k
S = concentration of component k in the crude oil of storage tanki.
ε
j
,
k
,
0
B= initial concentration of component k in the crude mix of charging tankj.
ε
j
B
,
min
⁡= minimum concentration of component k in the crude mix of charging tankj.
ε
j
B
,
max
⁡= maximum concentration of component k in the crude mix of charging tankj.In this paper, a two-level time structure is used to model a general dynamic production system. We consider the crude mixj to be scheduled over a finite planning horizon consisting of macroperiods m
=
1
⋯
N
M. Each macroperiod m is divided into time intervals with fixed length, whereT
mrepresents the starting time of macroperiod m. LetD
j
,
m
M, a random variable, denote the demand of crude mix j in macroperiod m. We approximate the crude mix demand as a multivariate normal distribution with a specific correlation structureD
M
~
N
(
μ
,
Σ
)with covariance matrixΣ, whereD
M
≝
[
D
j
,
m
M
]
,
μ
≝
[
μ
j
,
m
]. 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:
Subject to:(2) Vessel Arrival and Departure Operation Rules.Each vessel arrives at the docking station for unloading only once throughout the scheduling horizon:(1a)
∑
t
=
1
N
SCH
X
v
,
t
F
=
1
∀
v
.
Each vessel leaves the docking station only once throughout the scheduling horizon. Consider
(1b)
∑
t
=
1
N
SCH
X
v
,
t
L
=
1
∀
v
.
Equation for unloading initiation time:
(1c)
T
v
F
=
∑
t
=
1
N
SCH
t
X
v
,
t
F
∀
v
.
Equation for unloading completion time:
(1d)
T
v
L
=
∑
t
=
1
N
SCH
t
X
v
,
t
L
∀
v
.
Each crude vessel should start unloading after arrival time set in the planning level:
(1e)
T
v
F
≥
T
v
ARR
∀
v
.
Duration of the vessel unloading is bounded by the initial volume of oil in the vessel divided by the maximum unloading rate:
(1f)
T
v
L
-
T
v
F
≥
⌈
V
v
,
0
V
max
⁡
⌈
F
v
,
i
V
S
,
max
⁡
⌉
⌉
∀
v
.
⌈
⌉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:
(1g)
T
v
+
1
F
≥
T
v
L
∀
v
.
Unloading is possible between timeT
v
F andT
v
L:
(1h)
X
v
,
t
W
≤
∑
τ
=
1
t
X
v
,
τ
F
,
X
v
,
t
W
≤
∑
τ
=
t
SCH
X
v
,
τ
L
∀
v
,
t
.(3) Material Balance Equations for the Vessel.For vesselv, the oil volume at time t equals the initial volume minus transferred volume from vessel v to storage tanks up to timet:(2a)
V
v
,
t
V
=
V
v
,
0
V
-
∑
i
=
1
N
ST
∑
τ
=
1
t
F
v
,
i
,
τ
V
S
∀
v
,
t
.
Operating constraints on crude oil transfer rate from vessel v to storage tank i at timet
(2b)
F
v
,
i
VS
,
min
⁡
X
v
,
t
W
≤
F
v
,
i
,
t
VS
≤
F
v
,
i
VS
,
max
⁡
X
v
,
t
W
∀
v
,
i
,
t
.
The volume of crude oil transferred from vessel v to storage tanks during the scheduling horizon equals the initial crude oil volume of vesselv:
(2c)
∑
i
=
1
N
ST
∑
t
=
1
N
SCH
F
v
,
i
,
t
VS
=
V
v
,
0
V
∀
v
,
i
,
t
.(4) Material Balance Equations for Storage Tanks.The oil volume in storage tank i at time t equals the initial oil volume in tank i plus the oil volume transferred from vessels to tank i up to the time t and minus oil volume removed from tank i to charging tanks up to the timet:(3a)
V
i
,
t
S
=
V
i
,
0
S
+
∑
v
=
1
N
V
∑
τ
=
1
t
F
v
,
i
,
τ
VS
-
∑
j
=
1
N
BT
∑
τ
=
1
t
F
i
,
j
,
τ
SB
∀
i
,
t
Operating constraints on crude oil transfer rate from storage tank i to charging tank j at timet, the term 1
-
∑
l
=
1
NCDU
D
j
,
l
,
t denotes that if charging tank j is charging any CDU, there is no oil transfer from storage tank i to charging tankj. Consider
(3b)
F
i
,
j
SB
,
min
⁡
(
1
-
∑
l
=
1
N
CDU
D
j
,
l
,
t
)
≤
F
i
,
j
,
t
SB
≤
F
i
,
j
SB
,
max
⁡
(
1
-
∑
l
=
1
N
CDU
D
j
,
l
,
t
)
∀
i
,
j
,
t
.
Volume capacity limitations for storage tank i at timet
(3c)
V
i
S
,
min
⁡
≤
V
i
,
t
S
≤
V
i
S
,
max
⁡
∀
i
,
t
.(5) Material Balance Equations for Charging Tanks.The crude oil mix in charging tank j at time t equals the initial oil volume in charging tank j plus the crude oil transferred from storage tanks to charging tankj up to timet and minus the crude oil mix j charged into CDUs up to timet:(4a)
V
j
,
t
B
=
V
j
,
0
B
+
∑
i
=
1
N
ST
∑
τ
=
1
t
F
i
,
j
,
τ
SB
-
∑
l
=
1
N
CDU
∑
τ
=
1
t
F
j
,
l
,
τ
BC
∀
j
,
t
.
Operating constraints on mixed oil transfer rate from charging tank j to CDU l at timet
(4b)
F
j
,
l
BC
,
min
⁡
D
j
,
l
,
t
≤
F
j
,
l
,
t
BC
≤
F
j
,
l
BC
,
max
⁡
D
j
,
l
,
t
∀
j
,
l
,
t
.
Volume capacity limitations for charging tank j at timet
(4c)
V
j
B
,
min
⁡
≤
V
j
,
t
B
≤
V
j
B
,
max
⁡
∀
j
,
t
.(6) Material Balance Equations for Component
k
in Charging Tanks.The volume of component k in charging tank j at time t equals the initial component k in charging tank j plus component k in crude oil transferred from storage tanks to charging tank j up to the time t and minus component k in crude oil mix j transferred to CDUs up to the timet:(5a)
v
j
,
t
B
=
v
j
,
0
B
+
∑
τ
=
1
t
∑
i
=
1
N
ST
f
i
,
j
,
τ
SB
-
∑
τ
=
1
t
∑
l
=
1
N
CDU
f
j
,
l
,
τ
BC
∀
j
,
k
,
t
.
Operating constraints on volumetric flow rate of component k from storage tank i to charging tankj
(5b)
f
i
,
j
,
k
,
t
SB
=
F
i
,
j
,
t
SB
ε
i
,
k
S
∀
i
,
j
,
k
,
t
.
Operating constraints on volumetric flow rate of component k from charging tank j to CDUl
(5c)
F
j
,
l
,
t
BC
ε
j
,
k
B
,
min
⁡
≤
f
j
,
l
,
k
,
t
BC
≤
F
j
,
l
,
t
BC
ε
j
,
k
B
,
max
⁡
∀
j
,
l
,
k
,
t
.
Volume capacity limitations for component k in charging tank j at timet
(5d)
V
j
,
t
B
ε
j
,
k
B
,
min
⁡
≤
v
j
,
k
,
t
B
≤
V
j
,
t
B
ε
j
,
k
B
,
max
⁡
∀
j
,
l
,
k
,
t
.(7) Operating Rules for Crude Oil Charging.Charging tank j can charge at most one CDU at any timet:
(6a)
∑
l
=
1
N
CDU
D
j
,
l
,
t
≤
1
∀
j
,
t
.
CDU l can be charged only by one charging tank at any timet:
(6b)
∑
j
=
1
N
BT
D
j
,
l
,
t
=
1
∀
l
,
t
.
If CDU l is charged by crude oil mix j at timet− 1 and charged byj
′at timet, changeover cost is involved:
(6c)
Z
j
,
j
′
,
l
,
t
≥
D
j
′
,
l
,
t
+
D
j
,
l
,
t
-
1
-
1
j
,
j
′
(
j
≠
j
′
)
=
1
,
…
,
N
B
∀
l
,
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.


## 4. Chance Constraint Based Deterministic Transformation

### 4.1. Expectation of Penalty Cost

The expectation of the underproduction penalty cost term for crude mixj in macroperiod m is equal to
(7)
E
[
∑
j
=
1
N
BT
C
j
PEN
∑
m
=
1
N
M
max
⁡
(
0
,
D
j
,
m
M
-
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
)
]
,
where
(8)
C
j
,
m
≜
max
⁡
(
0
,
D
j
,
m
M
-
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
)
=
-
min
⁡
⁡
(
0
,
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
-
D
j
,
m
M
)
=
-
[
min
⁡
⁡
(
D
j
,
m
M
,
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
)
-
D
j
,
m
M
]
.
To facilitate the calculation of the expectation, the standardization of the normally distributed variables D
j
,
m
M and variables∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BCis 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
(9)
x
j
,
m
=
D
j
,
m
M
-
D
j
,
m
M
^
σ
j
,
m
,
whereD
j
,
m
M
^denotes the mean of D
j
,
m
M and σ
j
,
m the square root of its variance. In the same way, the “standardization” of the variables∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BCdefines
(10)
K
j
,
m
=
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
-
D
j
,
m
M
^
σ
j
,
m
.Using the method in paper [17], the underproduction penalty termC
j
,
mis charging into the following form:
(11)
C
j
,
m
=
-
[
-
f
(
K
j
,
m
)
+
(
1
-
Φ
(
K
j
,
m
)
)
K
j
,
m
]
,
wheref is the standardized normal distribution function and Φ denotes the cumulative probability function of a standard normal random variable.


### 4.2. 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 parameterC
j
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 j in macroperiod m assumes the following form:
(12)
Pr
[
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
≥
D
j
,
m
M
]
≥
β
j
,
m
.
This constraint, known as achance constraint, imposes a lower boundβ
j
,
mon the probability that the crude mix demand realization will be greater than the planned salesD
j
,
m
Mfor crude mix j in macroperiod m. 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 ofD
j
,
m
M, the chance constraint can be written equivalently as
(13)
Pr
[
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
-
D
j
,
m
M
^
σ
j
,
m
≥
D
j
,
m
M
-
D
j
,
m
M
^
σ
j
,
m
]
≥
β
j
,
m
,
Pr
⌈
K
j
,
m
≥
x
j
,
m
⌉
≥
β
j
,
m
.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 constraint can be replaced with the following deterministic equivalent expression:(14)
ϕ
(
K
j
,
m
)
≥
β
j
,
m
.The application of the inverse of the normal cumulative distribution functionϕ
-
1, which is a monotonically increasing function, yields
(15)
K
j
,
m
≥
ϕ
-
1
(
β
j
,
m
)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
(
β
j
,
m
).


### 4.3. 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 single-product 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:(16)
Pr
[
⋂
(
j
,
m
)
∈
I
m
(
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
≥
D
j
,
m
M
)
]
≥
β
m
I
m
=
{
(
j
,
m
′
)
∣
1
≤
j
≤
N
BT
,
m
≤
m
′
≤
M
}
,
whereI
mis the set of product-period (
j
,
m
) combinations whose simultaneous demand satisfaction with probability of at leastβ
mis sought andI
m
=
{
(
j
,
m
′
)
∣
1
≤
j
≤
N
BT
,
m
≤
m
′
≤
M
} is the set of all joint chance constraint. The joint probability,β
m, 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 (16) can be written as separated equivalent deterministic equations. Initially we choose the T constraints in a manner such that together they are more stringent than (16). The initial set of individual linear constraints equation (17) replacing (16) was obtained using the following argument. First we denote the event ∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
B
C
≥
D
j
,
m
MbyA
j
,
m and its complement event∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
<
D
j
,
m
M by A
j
,
m
C. From Boole’s inequality of probability theory [22] it is well known that
(17)
P
[
⋃
(
j
,
m
)
∈
I
m
A
j
,
m
]
≤
∑
j
=
1
J
∑
m
=
1
M
P
[
A
j
,
m
]
.
Thus it follows that ifP
[
A
t
C
]
≤
(
β
m
/
T
), m
=
1,2
,
…
,
M, then
(18)
P
[
⋂
(
j
,
m
)
∈
I
m
A
j
,
m
]
=
1
-
P
[
⋃
(
j
,
m
)
∈
I
m
A
j
,
m
]
≥
1
-
∑
j
=
1
J
∑
m
=
1
M
P
[
A
j
,
m
C
]
≥
1
-
(
1
-
β
m
)
.
Now, becauseD
j
,
m
Mis normally distributed with mean μ
j
,
m and standard deviation σ
j
,
m
,
P
[
A
j
,
m
c
]
≤
(
1
-
β
m
)
/
T is equivalent to P
[
D
j
,
m
M
>
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
]
≤
(
1
-
β
m
)
/
T, which in its turn is equivalent to
(19)
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
B
C
≥
μ
j
,
m
+
z
α
/
T
σ
j
,
m
,
wherez
α
/
Tis the 100
(
1
-
(
(
1
-
β
m
)
/
T
)
th percentile of the standard normal distribution. Subsequently we will suppress the subscript forz and replace (19) by the following equation:
(20)
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
B
C
≥
μ
j
,
m
+
z
σ
j
,
m
.
The initial value ofzwill be chosen to be z
=
z
α
/
T.Using the individual constraints equation (20) results in a more conservative solution than the joint chance constraint equation (16). Therefore, it is required to update the value of𝓏in order to find a solution which satisfies the required confidence level. The method in updating𝓏proposed in [23] is used, which is based on the interpolation of𝓏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𝓏value. The algorithm is examined in the next section.


## 4.1. Expectation of Penalty Cost

The expectation of the underproduction penalty cost term for crude mixj in macroperiod m is equal to
(7)
E
[
∑
j
=
1
N
BT
C
j
PEN
∑
m
=
1
N
M
max
⁡
(
0
,
D
j
,
m
M
-
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
)
]
,
where
(8)
C
j
,
m
≜
max
⁡
(
0
,
D
j
,
m
M
-
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
)
=
-
min
⁡
⁡
(
0
,
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
-
D
j
,
m
M
)
=
-
[
min
⁡
⁡
(
D
j
,
m
M
,
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
)
-
D
j
,
m
M
]
.
To facilitate the calculation of the expectation, the standardization of the normally distributed variables D
j
,
m
M and variables∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BCis 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
(9)
x
j
,
m
=
D
j
,
m
M
-
D
j
,
m
M
^
σ
j
,
m
,
whereD
j
,
m
M
^denotes the mean of D
j
,
m
M and σ
j
,
m the square root of its variance. In the same way, the “standardization” of the variables∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BCdefines
(10)
K
j
,
m
=
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
-
D
j
,
m
M
^
σ
j
,
m
.Using the method in paper [17], the underproduction penalty termC
j
,
mis charging into the following form:
(11)
C
j
,
m
=
-
[
-
f
(
K
j
,
m
)
+
(
1
-
Φ
(
K
j
,
m
)
)
K
j
,
m
]
,
wheref is the standardized normal distribution function and Φ denotes the cumulative probability function of a standard normal random variable.


## 4.2. 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 parameterC
j
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 j in macroperiod m assumes the following form:
(12)
Pr
[
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
≥
D
j
,
m
M
]
≥
β
j
,
m
.
This constraint, known as achance constraint, imposes a lower boundβ
j
,
mon the probability that the crude mix demand realization will be greater than the planned salesD
j
,
m
Mfor crude mix j in macroperiod m. 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 ofD
j
,
m
M, the chance constraint can be written equivalently as
(13)
Pr
[
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
-
D
j
,
m
M
^
σ
j
,
m
≥
D
j
,
m
M
-
D
j
,
m
M
^
σ
j
,
m
]
≥
β
j
,
m
,
Pr
⌈
K
j
,
m
≥
x
j
,
m
⌉
≥
β
j
,
m
.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 constraint can be replaced with the following deterministic equivalent expression:(14)
ϕ
(
K
j
,
m
)
≥
β
j
,
m
.The application of the inverse of the normal cumulative distribution functionϕ
-
1, which is a monotonically increasing function, yields
(15)
K
j
,
m
≥
ϕ
-
1
(
β
j
,
m
)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
(
β
j
,
m
).


## 4.3. 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 single-product 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:(16)
Pr
[
⋂
(
j
,
m
)
∈
I
m
(
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
≥
D
j
,
m
M
)
]
≥
β
m
I
m
=
{
(
j
,
m
′
)
∣
1
≤
j
≤
N
BT
,
m
≤
m
′
≤
M
}
,
whereI
mis the set of product-period (
j
,
m
) combinations whose simultaneous demand satisfaction with probability of at leastβ
mis sought andI
m
=
{
(
j
,
m
′
)
∣
1
≤
j
≤
N
BT
,
m
≤
m
′
≤
M
} is the set of all joint chance constraint. The joint probability,β
m, 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 (16) can be written as separated equivalent deterministic equations. Initially we choose the T constraints in a manner such that together they are more stringent than (16). The initial set of individual linear constraints equation (17) replacing (16) was obtained using the following argument. First we denote the event ∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
B
C
≥
D
j
,
m
MbyA
j
,
m and its complement event∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
<
D
j
,
m
M by A
j
,
m
C. From Boole’s inequality of probability theory [22] it is well known that
(17)
P
[
⋃
(
j
,
m
)
∈
I
m
A
j
,
m
]
≤
∑
j
=
1
J
∑
m
=
1
M
P
[
A
j
,
m
]
.
Thus it follows that ifP
[
A
t
C
]
≤
(
β
m
/
T
), m
=
1,2
,
…
,
M, then
(18)
P
[
⋂
(
j
,
m
)
∈
I
m
A
j
,
m
]
=
1
-
P
[
⋃
(
j
,
m
)
∈
I
m
A
j
,
m
]
≥
1
-
∑
j
=
1
J
∑
m
=
1
M
P
[
A
j
,
m
C
]
≥
1
-
(
1
-
β
m
)
.
Now, becauseD
j
,
m
Mis normally distributed with mean μ
j
,
m and standard deviation σ
j
,
m
,
P
[
A
j
,
m
c
]
≤
(
1
-
β
m
)
/
T is equivalent to P
[
D
j
,
m
M
>
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
]
≤
(
1
-
β
m
)
/
T, which in its turn is equivalent to
(19)
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
B
C
≥
μ
j
,
m
+
z
α
/
T
σ
j
,
m
,
wherez
α
/
Tis the 100
(
1
-
(
(
1
-
β
m
)
/
T
)
th percentile of the standard normal distribution. Subsequently we will suppress the subscript forz and replace (19) by the following equation:
(20)
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
B
C
≥
μ
j
,
m
+
z
σ
j
,
m
.
The initial value ofzwill be chosen to be z
=
z
α
/
T.Using the individual constraints equation (20) results in a more conservative solution than the joint chance constraint equation (16). Therefore, it is required to update the value of𝓏in order to find a solution which satisfies the required confidence level. The method in updating𝓏proposed in [23] is used, which is based on the interpolation of𝓏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𝓏value. The algorithm is examined in the next section.


## 5. 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.


### 5.1. 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 (20) is replaced with (16) for a fixed scalar 𝓏 to determine the set of new controlled variables. The algorithm starts by choosing a high value for the initial z-value as in (20), which makes the corresponding solution satisfy the demand with a probability level higher than𝒫
target
=
β
m.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 (16). A subregion adaptive algorithm proposed by Genz and Bretz [24] is employed to carry out multivariate integration which makes this evaluation feasible.Step 4 (update𝓏-value).
If the evaluated probability level is in theγneighborhood of𝒫
target, the algorithm terminates since the goal of finding a schedule that satisfies the demand with a probability of 𝒫
target is accomplished; otherwise the z value is updated and the previous steps are repeated to obtain another schedule.To update the𝓏-value the following algorithm is used. The goal is to find a 𝓏-value in (20) that provides a schedule such that the demand can be satisfied with a probability of 𝒫
target
=
β
m over the entire time horizon. The 𝓏-value needs to be obtained iteratively.Step 1 (obtain the upper and lower bounds of𝓏).
First we choose𝓏
=
𝓏
α in (20). Obviously it yields a lower bound for the correct 𝓏-value. We call it 𝓏
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 𝒫
lower. Next we choose an arbitrary large value for 𝓏. We denote it by 𝓏
upper and obtain a similar estimate of the probability of demand satisfaction. We denote this probability by 𝒫
upper.Step 2 (obtain the upper and lower bounds of𝒫
target).
In this step we obtain the upper percentiles of the univariate standard normal distribution for these probabilities and denote𝒫
upperand𝒫
lower by 𝓏
2and 𝓏
1, respectively. We also denote the corresponding percentile for the𝒫
target value by𝓏
target.Step 3 (update𝓏-value).
Based on these values the updated𝓏-value is obtained using the following linear interpolation formula:
(21)
𝓏
new
=
𝓏
lower
+
(
𝓏
target
-
𝓏
1
𝓏
2
-
𝓏
1
)
(
𝒵
upper
-
𝒵
lower
)
.
If the 𝒵
new value is lower than 𝒵
2 and higher than 𝒵
target, then we replace𝒵
2by𝒵
new. If it is lower than𝒵
target and higher than 𝒵
1, we replace 𝒵
1 by 𝒵
new. We repeat this process using (16) until𝒫
targetis reached.


### 5.2. Update of the Planning and Scheduling Policy

As discussed in Section3, the demandD
M
~
N
(
μ
,
Σ
)follows a multivariate normal distribution with a specific correlation structure. In order to decrease the uncertainty associated with the random variables D
j
,
m
M 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 macroperiodmis described in the following steps.Step 1 (update statistics).
Update future demand statistics based on present and previous product demand realizations.Step 2 (resolve the problem).
Resolve the scheduling problem form
′
=
m
+
1
,
…
,
M.Step 3 (updatem).
Setm
←
m
+
1. If m
≤
M, return to Step 1.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. Step2 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 variableD
j
,
m
Mis partitioned into two sets. The first set
(22)
P
=
{
D
j
,
m
′
M
:
1
≤
m
′
≤
m
,
m
=
1
,
…
,
M
}
contains the crude mix demands realized in the past periods m
′
=
1
,
…
,
m. The second set
(23)
F
=
{
D
j
,
m
′
M
:
m
+
1
≤
m
′
≤
M
,
i
=
1
,
…
,
M
}
contains the uncertain crude mix demands for the future periods. Based on this partitioning, the variance-covariance matrix can be expressed as follows:
(24)
Σ
=
(
Σ
P
P
Σ
P
F
Σ
F
P
Σ
F
F
)
,
whereΣ
PPandΣ
FFare the variance-covariance submatrices of the crude mix demands belonging to the sets P andF, respectively. The elements of submatrices Σ
PF
=
(
Σ
FP
)
T are the covariances between the elements of P and F
. Let
(25)
Ξ
=
(
ξ
i
1
,
ξ
i
2
,
…
,
ξ
i
t
)
denote the realizations (outcomes) of the crude mix demands of set P associated with past periods. The conditional means of the uncertain future crude mix demands include
(26)
μ
F
/
P
=
μ
F
+
Σ
F
P
Σ
P
P
-
1
(
Ξ
-
μ
P
)
,
where μ
F
/
P denotes the conditional demand expectations andμ
Pandμ
F are the past and future mean values. The new (conditional) variance-covariance matrix is equal to
(27)
Σ
F
/
A
=
Σ
F
F
-
Σ
F
P
Σ
P
P
-
1
Σ
F
P
.
A detailed treatise of conditional multivariate normal probability functions can be found in Tong, 1990 [18]. Note that while the conditional mean valueμ
F
/
A depends on the demand realizations in past periods, the conditional covariance matrix Σ
F
/
A 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.


## 5.1. 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 (20) is replaced with (16) for a fixed scalar 𝓏 to determine the set of new controlled variables. The algorithm starts by choosing a high value for the initial z-value as in (20), which makes the corresponding solution satisfy the demand with a probability level higher than𝒫
target
=
β
m.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 (16). A subregion adaptive algorithm proposed by Genz and Bretz [24] is employed to carry out multivariate integration which makes this evaluation feasible.Step 4 (update𝓏-value).
If the evaluated probability level is in theγneighborhood of𝒫
target, the algorithm terminates since the goal of finding a schedule that satisfies the demand with a probability of 𝒫
target is accomplished; otherwise the z value is updated and the previous steps are repeated to obtain another schedule.To update the𝓏-value the following algorithm is used. The goal is to find a 𝓏-value in (20) that provides a schedule such that the demand can be satisfied with a probability of 𝒫
target
=
β
m over the entire time horizon. The 𝓏-value needs to be obtained iteratively.Step 1 (obtain the upper and lower bounds of𝓏).
First we choose𝓏
=
𝓏
α in (20). Obviously it yields a lower bound for the correct 𝓏-value. We call it 𝓏
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 𝒫
lower. Next we choose an arbitrary large value for 𝓏. We denote it by 𝓏
upper and obtain a similar estimate of the probability of demand satisfaction. We denote this probability by 𝒫
upper.Step 2 (obtain the upper and lower bounds of𝒫
target).
In this step we obtain the upper percentiles of the univariate standard normal distribution for these probabilities and denote𝒫
upperand𝒫
lower by 𝓏
2and 𝓏
1, respectively. We also denote the corresponding percentile for the𝒫
target value by𝓏
target.Step 3 (update𝓏-value).
Based on these values the updated𝓏-value is obtained using the following linear interpolation formula:
(21)
𝓏
new
=
𝓏
lower
+
(
𝓏
target
-
𝓏
1
𝓏
2
-
𝓏
1
)
(
𝒵
upper
-
𝒵
lower
)
.
If the 𝒵
new value is lower than 𝒵
2 and higher than 𝒵
target, then we replace𝒵
2by𝒵
new. If it is lower than𝒵
target and higher than 𝒵
1, we replace 𝒵
1 by 𝒵
new. We repeat this process using (16) until𝒫
targetis reached.


## 5.2. Update of the Planning and Scheduling Policy

As discussed in Section3, the demandD
M
~
N
(
μ
,
Σ
)follows a multivariate normal distribution with a specific correlation structure. In order to decrease the uncertainty associated with the random variables D
j
,
m
M 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 macroperiodmis described in the following steps.Step 1 (update statistics).
Update future demand statistics based on present and previous product demand realizations.Step 2 (resolve the problem).
Resolve the scheduling problem form
′
=
m
+
1
,
…
,
M.Step 3 (updatem).
Setm
←
m
+
1. If m
≤
M, return to Step 1.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. Step2 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 variableD
j
,
m
Mis partitioned into two sets. The first set
(22)
P
=
{
D
j
,
m
′
M
:
1
≤
m
′
≤
m
,
m
=
1
,
…
,
M
}
contains the crude mix demands realized in the past periods m
′
=
1
,
…
,
m. The second set
(23)
F
=
{
D
j
,
m
′
M
:
m
+
1
≤
m
′
≤
M
,
i
=
1
,
…
,
M
}
contains the uncertain crude mix demands for the future periods. Based on this partitioning, the variance-covariance matrix can be expressed as follows:
(24)
Σ
=
(
Σ
P
P
Σ
P
F
Σ
F
P
Σ
F
F
)
,
whereΣ
PPandΣ
FFare the variance-covariance submatrices of the crude mix demands belonging to the sets P andF, respectively. The elements of submatrices Σ
PF
=
(
Σ
FP
)
T are the covariances between the elements of P and F
. Let
(25)
Ξ
=
(
ξ
i
1
,
ξ
i
2
,
…
,
ξ
i
t
)
denote the realizations (outcomes) of the crude mix demands of set P associated with past periods. The conditional means of the uncertain future crude mix demands include
(26)
μ
F
/
P
=
μ
F
+
Σ
F
P
Σ
P
P
-
1
(
Ξ
-
μ
P
)
,
where μ
F
/
P denotes the conditional demand expectations andμ
Pandμ
F are the past and future mean values. The new (conditional) variance-covariance matrix is equal to
(27)
Σ
F
/
A
=
Σ
F
F
-
Σ
F
P
Σ
P
P
-
1
Σ
F
P
.
A detailed treatise of conditional multivariate normal probability functions can be found in Tong, 1990 [18]. Note that while the conditional mean valueμ
F
/
A depends on the demand realizations in past periods, the conditional covariance matrix Σ
F
/
A 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.


## 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 Table1. 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.Table 1
System information.
Scheduling horizon
30
Number of vessel arrivals
9
Vessel
Arrival time
Amount of crude
Concentration of key component 1
Concentration of key component 2
Vessel 1
1
100
0.01
0.04
Vessel 2
4
100
0.03
0.02
Vessel 3
7
100
0.05
0.01
Vessel 4
11
100
0.01
0.04
Vessel 5
14
100
0.03
0.02
Vessel 6
17
100
0.05
0.01
Vessel 7
21
100
0.01
0.04
Vessel 8
24
100
0.03
0.02
Vessel 9
27
100
0.05
0.01
Number of storage tanks
3
Storage tank
Capacity
Initial oil amount
Initial concentration of key component 1
Initial concentration of key component 2
Storage tank 1
100
20
0.0167
0.0333
Storage tank 2
100
50
0.03
0.0225
Storage tank 3
100
70
0.0433
0.0133
Number of charging tanks
3
Charging tank
Capacity
Initial oil amount
Initial concentration of key component 1
Initial concentration of key component 2
Charging tank 1
100
30
0.0167
0.0333
Charging tank 2
100
50
0.03
0.0225
Charging tank 3
100
30
0.0433
0.0133
Number of CDUs
2
Costs involved in vessel operation
Unloading cost: 8, sea waiting cost: 5
Tank inventory costs
Storage tank: 0.05; charging tank: 0.08
Changeover cost for charged oil switch
50
The penalty cost for crude mix shortfalls
0.5Table 2
Mean values of the product demand.
Macroperiod
Product 1
Product 2
Product 3
1
90
98
95
2
100
80
98
3
100
90
98The alternative model formulation is defined and solved using the proposed solution procedure.Mode:(28)
C
COST
=
∑
V
=
1
N
V
C
V
UNLOAD
(
T
v
L
-
T
v
F
)
+
∑
v
=
1
N
V
C
V
SEA
(
T
v
F
-
T
v
ARR
)
+
∑
i
=
1
N
ST
C
i
INVST
∑
t
=
1
N
SCH
(
V
i
,
t
S
+
V
i
,
t
-
1
S
2
)
+
∑
j
=
1
N
B
T
C
j
INVBL
∑
t
=
1
N
SCH
(
V
j
,
t
B
+
V
j
,
t
-
1
B
2
)
+
∑
t
=
1
N
SCH
∑
j
=
1
N
BT
∑
j
′
=
1
N
BT
C
j
,
j
′
SETUP
∑
l
=
1
N
CDU
Z
j
,
j
′
,
l
,
t
+
∑
j
=
1
N
BT
C
j
PEN
×
∑
m
=
1
N
M
max
⁡
(
0
,
D
j
,
m
M
-
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
B
C
)
Subject to:  Pr
[
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
≥
D
M
j
,
m
]
≥
β
j
,
m
h
i
n
d
a
w
w
i
i
i
i
i
[
D
M
j
,
m
Pr
⋂
(
j
,
m
)
∈
I
m
(
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
hindawiiiiiiiiiiii
hindawiiiiiiiiii
≥
D
M
j
,
m
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
B
C
)
]
≥
β
m
h
i
n
d
a
w
w
i
i
i
i
i
i
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 (29). 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 variance-covariance matrices. Consider
(29)
Σ
=
(
1.0
0.0
-
0.2
0.2
0.0
0.3
0.0
0.0
0.0
0.0
1.0
0.0
0.0
0.3
0.0
0.0
0.6
0.2
-
0.2
0.0
1.0
0.0
0.1
0.0
0.0
0.2
0.0
0.2
0.0
0.0
1.0
0.0
0.0
0.0
0.2
0.0
0.0
0.3
0.1
0.0
1.0
0.2
0.2
0.3
0.0
0.3
0.0
0.0
0.0
0.2
1.0
0.0
0.2
0.1
0.0
0.0
0.0
0.0
0.2
0.0
1.0
0.2
0.2
0.0
0.6
0.2
0.2
0.3
0.2
0.2
1.0
0.0
0.0
0.2
0.0
0.0
0.0
0.1
0.2
0.0
1.0
)
.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β
m, for the set of all 9 crude oil demands, is set to 50%. For the value ofβ
m 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.Figure 2
Comparison for the total operating cost with and without consideration of demand uncertainty.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.


### 6.2. Revision of Planning Policy

Using the correlation matrix in (29), 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:
(30)
Pr
[
⋂
(
j
,
m
)
∈
I
m
(
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
≥
D
j
,
m
M
)
]
≥
0.5
.
The resulting optimization problem is solved in the entire time horizon. Next, based on (i) the demand realizations in 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 macroperiodm
∈
(
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.Table 3
Employed realizations of the random product demands.
Macroperiod
Product 1
Product 2
Product 3
1
92.5504
107.0773
83.6901
2
104.2967
82.7651
92.1586
3
98.5225
93.5261
114.0730Table 4
Updates of the product demand expectations during the process of policy revision.
Policy revision
MacroperiodM
Product 1
Product 2
Product 3
1
1
90
98
95
2
100
80
98
3
100
90
98
2
2
100.0438
81.0850
98.0650
3
100
92.3
99.6155
3
3
100.7265
92.1073
98.6984Table 5
Updates of the product demand standard deviations during the process of policy revision.
Policy revision
MacroperiodM
Product 1
Product 2
Product 3
1
1
4.7434
4.9497
4.8734
2
5
4.4721
4.9497
3
5
4.7434
4.9497
2
2
4.8947
4.2417
4.7120
3
5.0
3.6691
4.8497
3
3
4.8821
3.3937
4.8020Table 6
Updates of the product demand correlation matrix.
Policy revision
1
1.0
0.0
−0.2
0.2
0.0
0.3
0.0
0.0
0.0
0.0
1.0
0.0
0.0
0.3
0.0
0.0
0.6
0.2
−0.2
0.0
1.0
0.0
0.1
0.0
0.0
0.2
0.0
0.2
0.0
0.0
1.0
0.0
0.0
0.0
0.2
0.0
0.0
0.3
0.1
0.0
1.0
0.2
0.2
0.3
0.0
0.3
0.0
0.0
0.0
0.2
1.0
0.0
0.2
0.1
0.0
0.0
0.0
0.0
0.2
0.0
1.0
0.2
0.2
0.0
0.6
0.2
0.2
0.3
0.2
0.2
1.0
0.0
0.0
0.2
0.0
0.0
0.0
0.1
0.2
0.0
1.0
2
1
−0.0045
−0.0671
0
0.2531
0
−0.0045
1
0.2146
0.2109
0.1352
−0.0646
−0.0671
0.2146
1
0
0.2546
0.1072
0
0.2109
0
1
0.2586
0.2041
0.2531
0.1352
0.2546
0.2586
1
−0.1583
0
−0.0646
0.1072
0.2041
−0.1583
1
3
1.0000
0.2672
0.2312
0.2672
1.0000
−0.1970
0.2312
−0.1970
1.0000Two 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 shown in Figure3. Examination of the data reveals that the actual production amount using the updating policy is close to the actual realized demand.Figure 3
Comparison for the realized demand and the actual production for crude mix under three circumstances.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.


## 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 (29). 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 variance-covariance matrices. Consider
(29)
Σ
=
(
1.0
0.0
-
0.2
0.2
0.0
0.3
0.0
0.0
0.0
0.0
1.0
0.0
0.0
0.3
0.0
0.0
0.6
0.2
-
0.2
0.0
1.0
0.0
0.1
0.0
0.0
0.2
0.0
0.2
0.0
0.0
1.0
0.0
0.0
0.0
0.2
0.0
0.0
0.3
0.1
0.0
1.0
0.2
0.2
0.3
0.0
0.3
0.0
0.0
0.0
0.2
1.0
0.0
0.2
0.1
0.0
0.0
0.0
0.0
0.2
0.0
1.0
0.2
0.2
0.0
0.6
0.2
0.2
0.3
0.2
0.2
1.0
0.0
0.0
0.2
0.0
0.0
0.0
0.1
0.2
0.0
1.0
)
.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β
m, for the set of all 9 crude oil demands, is set to 50%. For the value ofβ
m 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.Figure 2
Comparison for the total operating cost with and without consideration of demand uncertainty.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.


## 6.2. Revision of Planning Policy

Using the correlation matrix in (29), 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:
(30)
Pr
[
⋂
(
j
,
m
)
∈
I
m
(
∑
l
=
1
N
CDU
∑
T
m
T
m
+
1
-
1
F
j
,
l
,
t
BC
≥
D
j
,
m
M
)
]
≥
0.5
.
The resulting optimization problem is solved in the entire time horizon. Next, based on (i) the demand realizations in 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 macroperiodm
∈
(
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.Table 3
Employed realizations of the random product demands.
Macroperiod
Product 1
Product 2
Product 3
1
92.5504
107.0773
83.6901
2
104.2967
82.7651
92.1586
3
98.5225
93.5261
114.0730Table 4
Updates of the product demand expectations during the process of policy revision.
Policy revision
MacroperiodM
Product 1
Product 2
Product 3
1
1
90
98
95
2
100
80
98
3
100
90
98
2
2
100.0438
81.0850
98.0650
3
100
92.3
99.6155
3
3
100.7265
92.1073
98.6984Table 5
Updates of the product demand standard deviations during the process of policy revision.
Policy revision
MacroperiodM
Product 1
Product 2
Product 3
1
1
4.7434
4.9497
4.8734
2
5
4.4721
4.9497
3
5
4.7434
4.9497
2
2
4.8947
4.2417
4.7120
3
5.0
3.6691
4.8497
3
3
4.8821
3.3937
4.8020Table 6
Updates of the product demand correlation matrix.
Policy revision
1
1.0
0.0
−0.2
0.2
0.0
0.3
0.0
0.0
0.0
0.0
1.0
0.0
0.0
0.3
0.0
0.0
0.6
0.2
−0.2
0.0
1.0
0.0
0.1
0.0
0.0
0.2
0.0
0.2
0.0
0.0
1.0
0.0
0.0
0.0
0.2
0.0
0.0
0.3
0.1
0.0
1.0
0.2
0.2
0.3
0.0
0.3
0.0
0.0
0.0
0.2
1.0
0.0
0.2
0.1
0.0
0.0
0.0
0.0
0.2
0.0
1.0
0.2
0.2
0.0
0.6
0.2
0.2
0.3
0.2
0.2
1.0
0.0
0.0
0.2
0.0
0.0
0.0
0.1
0.2
0.0
1.0
2
1
−0.0045
−0.0671
0
0.2531
0
−0.0045
1
0.2146
0.2109
0.1352
−0.0646
−0.0671
0.2146
1
0
0.2546
0.1072
0
0.2109
0
1
0.2586
0.2041
0.2531
0.1352
0.2546
0.2586
1
−0.1583
0
−0.0646
0.1072
0.2041
−0.1583
1
3
1.0000
0.2672
0.2312
0.2672
1.0000
−0.1970
0.2312
−0.1970
1.0000Two 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 shown in Figure3. Examination of the data reveals that the actual production amount using the updating policy is close to the actual realized demand.Figure 3
Comparison for the realized demand and the actual production for crude mix under three circumstances.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.


## 7. 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 multiple-crude 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.



---

*Source: 730314-2014-03-16.xml*