# Knowledge-Assisted Dual-Stage Evolutionary Optimization of Large-Scale Crude
Oil Scheduling

Wanting Zhang, Wei Du, Guo Yu, , Renchu He,
Wenli Du, and Yaochu Jin Manuscript received 8 March 2023; revised 9 August
2023; accepted 18 September 2023. This work was supported by National Natural
Science Foundation of China (Key Program: 62136003), Major Program of Qingyuan
Innovation Laboratory (Grant No. 00122002), National Natural Science
Foundation of China (62173144, 62373154) and the Shanghai Committee of Science
and Technology, China (Grant No.22DZ1101500). (Corresponding authors: Wei Du;
Wenli Du.)Wanting Zhang, Wei Du, Renchu He, Wenli Du, and Yaochu Jin are with
the Key Laboratory of Smart Manufacturing in Energy Chemical Process, Ministry
of Education, East China University of Science and Technology, Shanghai
200237, China (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>yaochu.jin@uni-bielefeld.de).Wenli
Du is also with the Qingyuan Innovation Laboratory, Quanzhou 362801,
China.Yaochu Jin is also with the Chair of Nature Inspired Computing and
Engineering, Faculty of Technology, Bielefeld University, 33619 Bielefeld,
Germany.Guo Yu is with the Institute of Intelligent Manufacturing, Nanjing
Tech University, Nanjing 211816, China (e-mail: gysearch@163.com).

###### Abstract

With the scaling up of crude oil scheduling in modern refineries, large-scale
crude oil scheduling problems (LSCOSPs) emerge with thousands of binary
variables and non-linear constraints, which are challenging to be optimized by
traditional optimization methods. To solve LSCOSPs, we take the practical
crude oil scheduling from a marine-access refinery as an example and start
with modeling LSCOSPs from crude unloading, transportation, crude distillation
unit processing, and inventory management of intermediate products. On the
basis of the proposed model, a dual-stage evolutionary algorithm driven by
heuristic rules (denoted by DSEA/HR) is developed, where the dual-stage search
mechanism consists of global search and local refinement. In the global search
stage, we devise several heuristic rules based on the empirical operating
knowledge to generate a well-performing initial population and accelerate
convergence in the mixed variables space. In the local refinement stage, a
repair strategy is proposed to move the infeasible solutions towards feasible
regions by further optimizing the local continuous variables. During the whole
evolutionary process, the proposed dual-stage framework plays a crucial role
in balancing exploration and exploitation. Experimental results have shown
that DSEA/HR outperforms the state-of-the-art and widely-used mathematical
programming methods and metaheuristic algorithms on LSCOSP instances within a
reasonable time.

###### Index Terms:

Large-scale crude oil scheduling, evolutionary optimization, heuristic rules,
local refinement

##  I Introduction

Production scheduling plays a vital role in improving the economic performance
and competitiveness of refineries [1]. Short-term crude oil scheduling is a
leading segment for the overall refinery scheduling. It aims to make the
short-term decision for a week or ten days, depending on the operation status
and material transfer. The scheduling quality directly affects the stability
of subsequent production. However, scheduling optimization for crude oil is
challenging with the interaction of discrete events and continuous processes.
Currently, crude scheduling still relies on the experience of schedulers due
to the lack of effective techniques or software tools [2]. Therefore, short-
term scheduling optimization for crude oil has attracted extensive attention
from both academics and industries.

Typically, mathematical programming (MP) models represented with discrete- [3]
or continuous-time [4] are widely formulated to describe short-term scheduling
problems of crude oil. The discrete models divide the whole scheduling horizon
into several periods with uniform time slots. The operating activities that
begin or end at the slots are represented by binary variables. It is
convenient for complex scheduling problems to model as discrete-time
representations. However, the commonly-used way to improve the accuracy by
using the small slots may lead to a considerable increase of binary variables,
which brings difficulties in computing. The continuous-time models using fewer
binary variables in single [5] or multiple time-grids models [6, 4] have been
reported to overcome this deficiency. Nevertheless, the number of events for
operations in these models is required to be given in advance and only
determined by trial and error, which is difficult to implement in real-world
applications. As a result, the practicality and optimality of the modeling
method are suggested to be balanced in large-scale complex problems.

In addition, the blending constraints involving bilinear terms in crude oil
scheduling optimization lead to mixed-integer nonlinear programming (MINLP).
As discussed in [7], current commercial solvers may fail to converge to an
optimal solution when solving MINLP problems. In order to solve this problem
efficiently, several strategies have been reported. The common way is to solve
mixed-integer linear programming (MILP) and nonlinear programming (NLP)
iteratively [8, 9, 10]. However, this approach may not find a feasible
solution due to the composition concentration discrepancy caused by the
relaxation of nonlinear terms [11]. Similarly, the work in [12] presents a
rolling-horizon framework, which eliminates the composition discrepancy, but
it may fail to get a feasible schedule since myopic behavior [6]. Therefore,
despite the fact that MP methods have been used extensively, it compromises
the originality of complex problems and the quality of solutions.

As introduced above, most research based on MP methods pays attention to
attaining a better schedule with fewer binary variables, nonlinear terms, and
less CPU time. However, with the increase in scheduling scale, the number of
discrete variables and nonlinear terms expands significantly, which leads to
an exponential growth of computational time [13]. In recent years,
metaheuristic algorithms as a powerful tool have been introduced to solve
nonlinear [14] and large-scale [15] optimization problems. These advantages
underscore the superiority of metaheuristics over traditional methods in
various production scheduling problems. Numerous studies have examined the use
of metaheuristics in scheduling, as reviewed in references [16, 17, 18]. More
specifically, Pan et al. [19] explored flexible job-shop scheduling and
developed a bi-population evolutionary algorithm with a feedback mechanism for
energy-efficient optimization. To address the exponential complexity in
distributed flow-shop scheduling optimization, Wang et al. [20] proposed a
metaheuristic algorithm based on a knowledge-based cooperative strategy, which
effectively tackles this strongly NP-hard problem. Furthermore, Du et al. [21]
introduced a decision variable classification to support evolutionary
optimization in solving high-dimensional order scheduling problems. Evidently,
these studies demonstrate that metaheuristics have emerged as effective
approaches for handling complex scheduling problems. In addition, another
attractive feature of metaheuristics compared to MP is that they can be easily
incorporated with sophisticated simulation models that introduce many
realistic details. Specifically, these details contain constraints that cannot
clearly be modeled by equations or would give rise to models with a large
number of variables but relatively few free decision variables.

Motivated by the advantages of metaheuristics, a number of evolutionary
algorithms (EAs) have been investigated on crude oil scheduling problems. For
instance, a structural adaptive genetic algorithm (SAGA) for crude scheduling
problems is proposed in [22] and further applied to the optimization for the
robust scheduling of crude oil [13, 23]. Though this graph-based
representation used in this algorithm provides a sparse problem
representation, the convergence of SAGA remains slower in large-scale problems
with a huge number of decision variables. The work in [24, 25] investigates
the multi-objective optimization problem for crude oil scheduling in an inland
refinery and applies GA to solve the assignment of charging tanks and
distillers. However, their methods are only implemented on condition that
storage tanks have sufficient inventory. The realistic scale scheduling
problems are far from being resolved. While crude oil scheduling problems have
attracted extensive EA-based research, there is a notable lack of emphasis on
the utilization of knowledge. In general, specific knowledge can be derived
from the characteristics of problems, and its application can significantly
enhance search efficiency for real-world problems [26]. In the context of
crude oil scheduling problems, empirical knowledge embedded within operations,
such as the preference for selecting crude types during the blending process,
is easily understandable and extractable. Consequently, more efforts should be
made to tackle large-scale complex scheduling based on metaheuristics and
improve the applicability of approaches for practical operations.

Based on the above discussions, it can be found that there is still a big gap
between the existing theoretical research and actual processes in large-scale
crude scheduling optimization. To this end, in this paper, we propose a novel
dual-stage evolutionary algorithm driven by heuristic rules (DSEA/HR) to
efficiently address the large-scale crude oil scheduling problems (LSCOSPs) in
the real-world refinery. The contributions of this work are listed as follows:

  1. 1.

In contrast to the small or medium-scale problems resolved in existing
studies, large-scale crude oil scheduling problems are modeled from crude
unloading, transportation, crude distillation unit processing, inventory
management of intermediate products to the practical processing and operating
constraints.

  2. 2.

In the DSEA/HR, the dual-stage search consists of global search and local
refinement. The former is to speed up the optimization process, while the
latter is to reform the infeasible solutions. Particularly, we extract
problem-specific knowledge by analyzing the impact of crude blending
operations on search space. This knowledge is utilized to formulate two
heuristic rules for population initialization, thereby improving the global
search efficiency within the mixed decision space.

  3. 3.

The proposed DSEA/HR has been successfully implemented on practical LSCOSP
instances, and the experimental results have shown that our approach is
superior to the compared methods in performance and computational efficiency.

The remainder of the paper is organized as follows. Section II formulates the
mathematical model of the LSCOSP. The proposed DSEA/HR approach based on
heuristic rules and a repair strategy is presented in Section III. Section IV
reports the experimental results on real-world cases, together with
performance analysis. Finally, the paper is concluded in Section V.

##  II Large-Scale Crude Oil Scheduling Problem

###  II-A Problem Description

Figure 1: Crude oil scheduling structure of marine-access refineries.

The scope of the studied LSCOSP is illustrated in Fig. 1. Firstly, crude oil
is transferred to tanks from very large crude carriers (VLCCs) or small-size
vessels. Then the crude from different tanks is blended and charged to crude
distillation units (CDUs) where various intermediate products such as residue,
diesel, and naphtha are produced. Finally, these intermediates are fed to
specific downstream units. Notably, we focus on a marine-access refinery in
this paper. The range of component concentration in the tank fluctuates widely
due to marine-access refineries without dedicated terminal storage tanks.
Compared with the scheduling of inland refineries, severe composition
discrepancy will arise due to the relaxation of nonlinear terms in
deterministic methods. In addition, the residue may exceed loads of secondary
units by overusing inexpensive heavy crudes [27]. To address this issue, we
have also taken into account the inventory of four types of residue tanks. The
residue types are related to the operating modes in CDUs and crude types.
Apparently, the above operational features undoubtedly make resource
allocation more complex.

###  II-B Problem Formulation

Based on the above problem description, we present the mathematical model of
the LSCOSP in this section, utilizing the notations listed in Table I. It is
worth noting that we chose the discrete-time representation of the LSCOSP for
the following reasons: (1) this representation is extensive and effective for
complex industrial problems on a large scale [28]; (2) the discrete-time
representation conforms with the conditions of the marine-access refineries
where the practical operations only occur at regular time points. Next, we
formulate the constraints and objective function as follows.

TABLE I: Nomenclature.

Notation |  Meaning
---|---
Indices and sets: |
$v$ |  Index of vessels, $v\in\\{1,...,V\\}$
$t$ |  Index of tanks, $t\in\\{1,...,T\\}$
$c$ |  Index of crude types, $c\in\\{1,...,C\\}$
$k$ |  Index of crude properties, $k\in\\{1,...,K\\}$
$u$ |  Index of CDUs, $u\in\\{1,...,U\\}$
$r$ |  Index of residue types, $r\in\\{1,...,R\\}$
$s$ |  Index of product types, $s\in\\{1,...,S\\}$
$n$ |  Index of periods, $t\in\\{1,...,N\\}$
$RC_{r}$ |  Set of crude allowed to product residue $r$
Parameters: |
$AD_{v}$ |  Arrival day of vessel $v$
${MT}_{u}$ |  The maximum number of the charging tank connected by CDU $u$ in each period
${CR}_{r,n}$ |  Consumption rate of residue type $r$ in period $n$
${FU}_{u}^{L}$ |  Lower bound of crude received by CDU $u$ in one period
${FU}_{u}^{U}$ |  Upper bound of crude received by CDU $u$ in one period
$P_{c,k}$ |  Key property $k$ of crude $c$
$P_{u,k}^{L}$ |  Lower bound of feed property $k$ of CDU $u$
$P_{u,k}^{U}$ |  Upper bound of feed property $k$ of CDU $u$
${TL}_{t}^{L}$ |  Lower bound of liquid level of crude tank $t$
${TL}_{t}^{U}$ |  Upper bound of liquid level of crude tank $t$
${IR}_{r}^{L}$ |  Lower bound of inventory of residue $r$
${IR}_{r}^{U}$ |  Upper bound of inventory of residue $r$
${TS}_{t}$ |  Area of bottom of tank $t$
${yield}_{u,c,s}$ |  Product $s$ yield of crude $c$ in CDU $u$
${FO}_{s,u}^{L}$ |  Lower bound of the outflow of product $s$ from CDU $u$
${FO}_{s,u}^{U}$ |  Upper bound of the outflow of product $s$ from CDU $u$
Binary variables:
${XC}_{v,c,t,n}$ |  =1, if vessel $v$ transfers crude $c$ to tank $t$ in period $n$, 0 otherwise
${XT}_{t,u,n}$ |  =1, if tank $t$ is charging to CDU $u$ in period $n$, 0 otherwise
${XU}_{u,r,n}$ |  =1, if CDU $u$ produces residue type $r$ in period $n$, 0 otherwise
Continuous variables:
${CTU}_{t,u,n}$ |  0-1 continuous variables; =1, if connection status of tank $t$ to CDU $u$ changes between period $n-1$ and period $n$, 0 otherwise
${CUR}_{u,r,n}$ |  0-1 continuous variables; =1, if connection status of residue type $r$ in CDU $u$ changes between period $n-1$ and period $n$, 0 otherwise
$CO_{u,n}$ |  0-1 continuous variables; =1, if connection status changes between period $n\\!-\\!1$ and period $n$, 0 otherwise
${IV}_{v,c,n}$ |  Amount of crude $c$ in vessel $v$ at the end of period $n$
${ITC}_{t,c,n}$ |  Amount of crude $c$ in tank $t$ at the end of period $n$
${FV}_{v,c,t,n}$ |  Amount of crude $c$ unloaded by vessel $v$ to tank $t$ at the end of period $n$
${FR}_{u,r,n}$ |  Weight of reside type $\\!r\\!$ produced by CDU $\\!u\\!$ at the end of period $n$
${IT}_{t,n}$ |  Inventory of crude tank $t$ at the end of period $n$
${IR}_{r,n}$ |  Inventory of residue $r$ at the end of period $n$
${FTC}_{t,c,u,n}$ |  Amount of crude $c$ charged by tank $t$ to CDU $u$ in period $n$
${FUC}_{u,c,n}$ |  Total amount of crude $c$ received by CDU $u$ in period $n$
${FU}_{u,n}$ |  Total feed of CDU $u$ in period $n$
${FT}_{t,u,n}$ |  Total amount of crude charged by tank $t$ to CDU $u$ in period $n$
$FTC_{t,c,u,n}$ |  Amount of crude $c$ received by CDU $u$ from tank $t$ in period $n$
${FO}_{u,s,n}$ |  Weight of product $s$ output from CDU $u$ in period $n$

####  II-B1 Constraints

In the LSCOSP, two groups of constraints based on discrete events and
continuous processes are concluded:

  * •

Operating constraints, which are denoted by the binary variables representing
the operating requirements.

  * •

Transfer constraints, which include material balance, capacity, and processing
conditions.

Given that a complete mathematical model involves numerous formulations, in
order to improve the readability of this paper, we have provided the primary
constraints, while the remaining constraints can be found in Section S.I of
the supplementary file.

  * a)

Operating constraints:

(1) (2) (3) (4) (5)

where Eq. (1) specifies that each crude tank is capable of receiving only one
type of crude from a vessel during each period. This constraint is in place to
prevent the mixing of different quality crudes in the same tank during
unloading. Eq. (2) ensures that the inlet and outlet operations of each tank
cannot be implemented simultaneously. Eq. (3) means that no more than two CDUs
can charge a single tank simultaneously. Eq. (4) defines the maximum number of
the charging tank that can be used. Eq. (5) indicates that each CDU can only
produce one type of residue during each period. The remaining operating
constraints are provided in Eqs. (S.1)-(S.11) of the supplementary file.

  * b)

Transfer constraints:

(6) (7)

where Eq. (6) assures that the concentration of the crude mix composition in
the tanks is equal to that in the CDU feed streams. Noteworthy, this equation
contains two non-convex bilinear terms, leading to poor convergence in large-
scale problems. Eq. (7) requires that the key properties (i.e., density,
sulfur content, and total acid number) of the feed should be satisfied in the
allowable range after the crude blending process. The remaining transfer
constraints are listed in Eqs. (S.12)-(S.23) of the supplementary file.

####  II-B2 Objective function

Minimizing the total operating costs is commonly adopted as the optimization
objective, which includes the vessel waiting cost, the inventory cost, and the
changeover cost [11]. To ensure efficient production capacity, the unloading
and processing rates are both regarded as maximum limits in the LSCOSP.
Accordingly, the objective function can be simplified to minimize the
changeover cost and is presented as follows:

(8)

where $\omega$ is the cost coefficient for the single changeover. Eq. (8)
indicates that the changeover cost is incurred by switching feed among
different tanks and switching operating modes in CDUs.

In this way, the MINLP formulation for the LSCOSP can be described as follows:

$\displaystyle\left\\{\begin{array}[]{l}\min\,\,\,F\\\
s.t.\quad\rm{Eqs.}\eqref{a-1}-\eqref{b-2},\\\
\qquad\,\,\rm{Eqs.}(S.1)-(S.23).\end{array}\right.$ (12)

###  II-C Illustrative Instance

We provide a simple example to illustrate the crude oil scheduling process.
This example involves a single vessel ($V_{1}$), four crude tanks
($T_{1},T_{2},T_{3},T_{4}$), three types of crude oil ($C_{1},C_{2},C_{3}$),
and one CDU ($CDU_{1}$) capable of producing two types of residue
($R_{1},R_{2}$) in a refinery. The scheduling horizon spans three days. Table
II lists crucial parameters, including the initial quantities of crude oil
($IV$) of the vessel and its arrival time ($AD$). The table also presents the
initial inventory of tanks ($ITC$), the sulfur content ($P_{sul}$), and the
residue yield ($Yld$) of each crude type. Additionally, the last row in this
table outlines the permissible range of crude oil types for processing with
each residue type. To be specific, vessel $V_{1}$ carries 60 kt of crude
$C_{2}$ and is scheduled to arrive on the second day. The initial tank
inventory comprises 50 kt of crude $C_{3}$, 30 kt of crude $C_{1}$, 20 kt of
crude $C_{2}$, and 20 kt of crude $C_{1}$ for tanks $T_{1}$-$T_{4}$,
respectively. The sulfur contents of crude $C_{1}$-$C_{3}$ are 0.5%, 1.1%, and
2.7%, respectively, with corresponding residue yields of 33.64%, 13.65%, and
29.65%. Due to the unique design regulations of each CDU, they can produce
different types and quantities of residue by blending various crude oils. In
the given example, CDU1 can produce residue $R_{1}$ by utilizing crude
$C_{1}$-$C_{3}$, and it can also generate residue type $R_{2}$ by blending
crude $C_{1}$ and $C_{2}$.

Fig. 2 illustrates a feasible solution with two changeovers over three days
for the above scenario, which meets all constraints described in Section II-B.
More specifically, Fig. 2(a) clearly demonstrates the feasibility of the
operating constraints. Figs. 2(b) and 2(c) reflect the alteration of the feed
status (e.g., flow rate and sulfur content) supplied to the CDU and the
corresponding residue outputs during each changeover. It is evident that
ensuring the feed status complies with specific boundary limitations and
adjusting the residue processing mode for the CDU to meet the production
demand are both critical for maintaining scheduling feasibility in terms of
materials transfer.

However, it is important to note that the number of discrete variables and
nonlinear terms significantly increases with the expansion of the scheduling
horizon and the number of resources, such as tanks, crude types, CDUs, and
product types. This expansion results in exponential time complexity and can
even exceed the solving capacity of MINLP solvers. Given the advantages of EAs
in solving large-scale and nonconvex problems, we propose a novel EA-based
approach called DSEA/HR in the following section for solving the presented
LSCOSP.

TABLE II: Numerical Example of the Studied Problem.

Vessel ($c|IV|AD$) |  $V_{1}$
---|---
|  ($C_{2}|$ 60 kt $|$ 2nd day)

Tank ($c|ITC$) |  $T_{1}$ |  $T_{2}$ |  $T_{3}$ |  $T_{4}$
---|---|---|---|---
|  ($C_{3}|$50 kt) |  ($C_{1}|$30 kt) |  ($C_{2}|$20 kt) |  ($C_{1}|$20 kt)

Crude ( $P_{sul}|Yld$) |  $C_{1}$ |  $C_{2}$ |  $C_{3}$
---|---|---|---
|  (0.5%$|$33.64%) |  (1.1%$|$13.65%) |  (2.7%$|$29.65%)

Residue ${c}$ |  $RC_{R1}$ |  $RC_{R2}$
---|---|---
|  $\\{C_{1},C_{2},C_{3}\\}$ |  $\\{C_{1},C_{2}\\}$

Figure 2: An illustration of an feasible schedule. (a) Gantt chart of a
feasible solution. (b) Flow rate and sulfur content of feed to the CDU. (c)
Outputs of two types of residue in the CDU.

##  III Dual-Stage Evolutionary Optimization for LSCOSPs

For large-scale problems, metaheuristic algorithms are popularly used since
the above MINLP formulations are difficult to converge. This section
introduces a dual-stage evolutionary algorithm called DSEA/HR to solve the
LSCOSP rather than directly solving the large-scale MINLP model. The core of
the proposed DSEA/HR lies in its dual-stage search mechanism, which includes a
global search and a local refinement stage. According to the diverse search
characteristics of these two stages, we employ the competitive swarm optimizer
(CSO) [29] and the composite differential evolution (CoDE) [30] optimizer as
the fundamental optimizers in global and local optimization, respectively. It
is worth noting that the key performance improvement stems from the dual-stage
search framework itself, not limited to any specific optimizer. The
implementation details of the proposed DSEA/HR are elaborated below.

###  III-A Encoding and decoding schemes

In solving the LSCOSP, the encoding scheme is essential for EAs. A potential
solution should reflect the operating target tanks and their flow rate. To
this end, the structure of the encoding scheme for one period is comprised of
four components: (1) receiving tanks, (2) receiving flow rate, (3) charging
tanks, and (4) charging flow rate. They are illustrated in Fig. 3. In
Component I, every two bits ${RT}_{v,1}$ and ${RT}_{v,2}$ denote receiving
tanks. Similarly, in Component II, every two bits ${RF}_{v,1}$ and
${RF}_{v,2}$ indicate receiving flow rate. The size of them is $2|V|$ in one
period, where $|V|$ is the total number of vessels. Moreover, in Components
III and IV, every single bit represents a charging tank and its charging flow,
respectively. The size of them is the sum of ${MT}_{u}$ for all CDUs.

Based on the above encoding, a complete schedule result can be mapped clearly.
In every period, the unloading tanks and their amounts for each vessel are
respectively determined by Component I and II from left to right one by one.
Note that the effectiveness of bits depends on whether the vessel $v$ is
unloading. The charging tanks and their amounts for $u$-th CDU are
respectively indicated by {$CT_{u,1}$, … ,$CT_{u,FT_{u}}$} and {$CF_{u,1}$, …
,$CF_{u,{MT}_{u}}$} represented in Component III and IV. To guarantee the
sustainability of activities in the entire scheduling horizon, each vessel can
unload from two tanks in a specific order during each period. Additionally,
the duration of charging operations must be in integer periods, implying
switches are only permitted at the beginning of the next period.

###  III-B Properties of LSCOSPs

In the LSCOSPs, each type of crude oil is limited to producing a specific
range of residue types. For instance, let $(C_{\Delta},R_{\Delta})$ denote the
set of crude types in a tank ($\Delta\\!=\\!t$) or a vessel ($\Delta\\!=\\!v$)
and their corresponding residue types that can be produced, respectively. If
$(C_{1},R_{1})$ and $(C_{2},R_{2})$ satisfy the conditions
$R_{1}\\!\\!=\\!\\!\\{1,2\\}$ and $R_{2}\\!\\!=\\!\\!\\{2,3\\}$, then the
blended crude $C_{blend}\\!=\\!\\{C_{1}\cup C_{2}\\}$ can only produce the
residue type $R_{blend}\\!=\\!\\{1,2\\}\cap\\{2,3\\}=\\{2\\}$. According to
this principle, two properties of the LSCOSPs can be derived and used as
valuable knowledge to enhance search efficiency.

Property 1: When assigning discharging tanks, it is important to avoid
significant disparities in the producible residue types between the crude oil
in the discharging tank and that in the vessel.

Property 2: After selecting charging tanks for CDUs, it is crucial to ensure
the diversity of residue types produced by other storage tanks.

Notably, formulating the above properties in terms of equations is
complicated. Therefore, we propose several heuristic rules that represent the
knowledge as probabilities and incorporate them into the search process.

Figure 3: Representation of encoding for each period.

###  III-C Population initialization with heuristic rules

In order to improve convergence for solving LSCOSPs, two heuristic rules based
on empirical knowledge are utilized in the population initialization of
DSEA/HR. Accordingly, Rules I and II are detailed below, followed by the
procedure of population initialization shown in Algorithm 1.

Rule I is devoted to the assignment of discharging tanks. According to
Property 1, randomly blending different types of crude may easily lead to an
infeasible schedule since different types of residue are only produced by
specific crudes. In this work, tanks are given different priorities for
discharging. Firstly, the highest priority is assigned to the empty tank. If
there is no empty tank, we define the similarity between the unloaded crude
and the mixture in the tank, formulated in Eq. (13). Then the tanks with the
highest similarity are chosen because the tanks with higher similarity own a
wider processable range of residue types than others after mixing the unloaded
crude. Finally, the tank with the largest available capacity is selected to
receive among the tanks with the highest similarity.

${R_{v,c,t}}=\left|{\bigcap\limits_{c\in\left\\{{\mathop{UJ}\nolimits_{v},\mathop{BJ}\nolimits_{t}}\right\\}}{CP}_{c}}\right|,\forall
v,t,$

(13)

where ${R_{v,c,t}}$ represents similarity between crude types in the tank t
and the crude c to be unloaded from vessel $v$; ${CP}_{c}$ denotes the range
of types in which the residue produced by crude c; ${UJ}_{v}$ and ${BJ}_{t}$
are the types of crude unloaded from vessel v and stored in tank t,
respectively.

Algorithm 1 Population initialization.

0: $M$ (population size), problem parameters

0: $P(M)$ (initialized population)

1: for $m=1\to M$ do

2: Randomly initialize $\\{CF_{1,1},\cdots,CF_{u,{MT}_{u}}\\}$ as charging
amount for each CDU;

3: Generate $\\{CT_{1,1},\cdots,CT_{u,{MT}_{u}}\\}$ as the charging tanks for
each CDU $\leftarrow$ Rule II and $\\{CF_{1,1},\cdots,CF_{u,{MT}_{u}}\\}$;

4: Generate $\\{RT_{1,1},\cdots,RT_{v,2}\\}$ as receiving tanks for each
vessel $\leftarrow$ Rule I and initial capacity of tanks;

5: Calculate $\\{RF_{1,1},\cdots,RF_{v,2}\\}$ as unloading amount based on
maximum capacity of $\\{RT_{1,1},\cdots,RT_{v,2}\\}$;

6: Add initial solution $m$ into population $P(M)$;

7: end for

Rule II aims to generate a group of charging tanks for each CDU. Typically,
the crudes in a tank can produce multiple types of residue. According to
Property 2, it is promising to satisfy the inventory of residue if the tanks
that can produce a variety of residue types are reserved. In this way, the
probability of selection for the charging tank is presented in Eq. (14). The
tank producing fewer types of residue is preferred to be selected, ensuring a
wide feasible region in the later periods.

$Prob_{t}=\frac{{\frac{1}{{\left|{\bigcap\limits_{c\in{BJ}_{t}}{CP}_{c}}\right|}}}}{{\sum\limits_{t^{\prime}\in
AT}{\frac{1}{{\left|{\bigcap\limits_{c\in{BJ}_{t^{\prime}}}{CP}_{c}}\right|}}}}},t\in
AT,$

(14)

where $AT$ indicates a set of crude tanks whose level is not fall below the
lower bound after they charge CDU with amount specified by $CF$.

###  III-D Constraints handling and fitness evaluation

In this paper, the LSCOSP is formulated as a large-scale constrained
optimization problem:

$\begin{aligned} \mathop{\min}\limits_{x}\quad&f(x)\\\
s.t.\quad&x=(x_{1},\cdots,x_{D})\in\Omega\\\ &{g_{i}}(x)\leq 0,i=1,\cdots,p\\\
&{h_{i}}(x)=0,i=p+1,\cdots,q\end{aligned},$ (15)

where $x$ is a solution including $D$ decision variables; $\Omega$ is the
decision space; ${g_{i}}(x)$ and ${h_{i}}(x)$ are the $i$-th inequality and
$(i-p)$-th equality constraint, respectively; $f(x)$ is the objective value of
solution $x$. Eqs. (16) and (17) respectively define the degree of violation
and the binary flag for $i$-th constraint.

$\displaystyle\begin{split}CV_{i}(x)&=\left\\{\\!\\!\\!{\begin{array}[]{*{20}{l}}\max(0,{g_{i}}(x)),&\\!\\!\\!\\!if\,i\leq
p\\\
\max(0,\left|{{h_{i}}(x)}\right|),&\\!\\!\\!\\!otherwise\end{array}}\right.\\!\\!\\!,i=1,\ldots,q,\end{split}$
(16)
$\displaystyle\begin{split}CVN_{i}(x)&=\left\\{{\begin{array}[]{*{20}{l}}1,&CV_{i}(x)>0\\\
0,&otherwise\end{array}}\right.,i=1,\ldots,q,\end{split}$ (17)

Accordingly, the number and the degree of constraint violations are calculated
by Eqs. (18) and (19), respectively. Meanwhile, the objective value of
individual $x$ is obtained from $f(x)$. In this way, the fitness of individual
$x$ can be denoted as $(CVN(x),CV(x),f(x))$.

$CVN(x)=\sum\nolimits_{i=1}^{q}{C{VN_{i}}(x)},$

(18)

$CV(x)=\sum\nolimits_{i=1}^{q}{C{V_{i}}(x)}.$

(19)

Further, the feasibility-based criterion [31] is used in the solution
selection to accelerate the convergence of the population to feasible regions.
To be specific, the fitness of solution $x$ is better than another solution
$y$ if the following conditions hold:

$\displaystyle\begin{split}\left\\{\\!\\!{\begin{array}[]{*{20}{l}}\scalebox{0.9}{$\\!f(x)\\!<\\!f(y),$}&\scalebox{0.9}{$\\!\\!if\,CVN(x),CVN(y)\\!=0$}\\\
\scalebox{0.9}{$\\!CV(x)\\!<\\!CV(y),$}&\scalebox{0.9}{$\\!\\!if\,CVN(x)\\!=CVN(y)$}\\\
\scalebox{0.9}{$\\!f(x)\\!<\\!f(y),$}&\scalebox{0.9}{$\\!\\!if\,CV(x)\\!=CV(y)$}\\\
\scalebox{0.9}{$\\!CVN(x)\\!<\\!CVN(y),$}&\scalebox{0.9}{$\\!\\!otherwise$}\end{array}}\right..\end{split}$
(20)

The above selection operation is adopted into both the global and local search
stages.

###  III-E Global search stage for LSCOSPs

In the dual-stage framework of the proposed DSEA/HR, the global search serves
as the first stage. Its goal is to evolve populations within a mixed decision
space that involves large-scale discrete and continuous variables. Given that
the performance of conventional optimizers, such as the particle swarm
optimizer (PSO) [32], seriously degrades as the dimension of the search space
increases. It is more appropriate to embed optimizers specifically designed
for large-scale global optimization problems into the global search stage of
DSEA/HR. To this end, we have adopted CSO [29] as the basic optimizer for the
global search stage in this paper. CSO is well-suited for solving widely used
high-dimension problems due to its novel competition mechanism.

In CSO, neither the particle’s best position nor the global best position from
the basic PSO is used. Instead, the particles are randomly grouped in pairs.
Then the particles from the same pair compete with each other. The particle
with better fitness will be selected to proceed to the next generation. The
other one will be updated by the modified velocity and position learned from
the winner in this competition. The updated strategies are presented in Eqs.
(21) and (22).

$\displaystyle\begin{split}\scalebox{0.9}{${V_{l,n}}(g+1)=$}&\scalebox{0.9}{${R_{1}}(n,g){V_{l,n}}(g)$}\\\
&\scalebox{0.9}{$+{R_{2}}(n,g)({X_{w,n}}(g)-{X_{l,n}}(g))$}\\\
&\scalebox{0.9}{$+\varphi{R_{3}}(n,g)(\overline{{X_{n}}}(g)-{X_{l,n}}(g)),$}\end{split}$
(21)

${X_{l,n}}(g+1)={X_{l,n}}(g)+{V_{l,n}}(g+1),$

(22)

where $V_{l,n}(g+1)$ and $X_{l,n}(g+1)$ denote the velocity and position of
the loser in the $n$-th round of the competition in generation $t$+1,
respectively; ${R_{1}}(n,g)$, ${R_{2}}(n,g)$, ${R_{3}}(n,g)$ are vectors
randomly generated between 0 to 1 after the n-th competition in generation g;
$\overline{{X_{n}}}(g)$ is the mean value vector of all particles in
generation g for each dimension; $\varphi$ is a control parameter, which is
set to zero as recommended in [29].

The above strategies for updating velocity and position are used to generate
the swarm of the next generation in the global search stage. Algorithm 2
presents the procedure of the global search with CSO.

Algorithm 2 Procedure of the global search.

0: $P_{0}(M)$ (initialized population), maximum evaluation number

0: $S_{CSO}(Z)$ (best solutions by CSO)

1: $g=0$;

2: /*Population evolution with CSO */

3: while the maximum evaluation number is not achieved do

4: Calculate the $fitness$ of each particle $X_{i}(g)$ in $P_{g}(M)$ using
(15)-(19);

5: $A=P_{g}(M)$, $P_{g+1}=\varnothing$;

6: while $A\neq\varnothing$ do

7: Randomly select two particles $X_{1}(g)$, $X_{2}(g)$ from $A$;

8: if $fitness(X_{1}(g))>fitness(X_{2}(g))$ then

9: $X_{w}(g)=X_{1}(g)$, $X_{l}(g)=X_{2}(g)$;

10: else

11: $X_{w}(g)=X_{2}(g)$, $X_{w}(g)=X_{1}(g)$;

12: end if

13: $P_{g+1}(M)\leftarrow$ $X_{w}(g)$;

14: Update $X_{l}(g)$ using (21) and (22);

15: $P_{g+1}(M)\leftarrow$ Updated $X_{l}(g)$;

16: Remove $X_{1}(g)$, $X_{2}(g)$ from $A$;

17: end while

18: $g=g+1$;

19: end while

20: $S_{CSO}(Z)\leftarrow$ Select solutions with the best fitness from
$P_{g}(M)$;

###### Remark 1

It is worth mentioning that the proposed dual-stage evolution framework is
flexible and can accommodate any optimizer. Thus, the choice of the employed
CSO for the proposed DSEA/HR is not unique. While the optimizer for the global
stage does not play a decisive role in terms of performance improvement, using
an optimizer tailored for large-scale global optimization is helpful in
improving the algorithm’s efficacy. This advantage is also demonstrated by the
detailed experimental comparison provided in Section IV-C2.

###  III-F Repair strategy based on local refinement

After optimizing the global search stage for the LSCOSPs, the discrete
variables representing discharging and charging tanks in the elite solutions
are expected to converge, aided by the guidance of heuristic rules. However,
finding a feasible solution in the global search stage remains challenging due
to the complexity of the mixed decision space. To mitigate this, we have
developed a local refinement stage with a problem-specific repair strategy for
handling infeasible solutions.

To enhance the exploitation capability of DSEA/HR, we have focused on
optimizing continuous variables. We have introduced a repair strategy that
fixes the discrete variables, effectively transforming the mixed-variable
problem into a continuous-variable one. This strategy, focusing on local
variables (continuous variables), is termed local refinement. Specifically,
during the local refinement stage, the values of discrete variables inherit
the results obtained from the global search, while the charging flow
represented by the continuous variable will be optimized in the EA-based
framework. Note that the EA-based optimizer cannot be replaced by any
deterministic nonlinear solver. This is because the preservation of the
diversity of feasible solutions is essential in solving the scheduling problem
effectively.

Algorithm 3 Procedure of the local refinement.

0: $S_{CSO}(Z)$ (best solutions by CSO), $H$ (population size), maximum
evaluation number

0: $S_{best}$ (final solution)

1: if $S_{CSO}(Z)$ violate processing constraints then

2: /* Population initialization for local refinement */

3: $i=1$, $n=0$, $g=0$;

4: for $i=1\to H$ do

5: Randomly initialize $\\{CF_{1,1},\cdots,CF_{u,{MT}_{u}}\\}$ as charging
amount for each CDU for individual $i$ ;

6: Duplicate charging tanks for each CDU
$\\{CT_{1,1},\cdots,CT_{u,{MT}_{u}}\\}$ from the particle $n$ to the
individual $i$;

7: Add individual $i$ into $Q_{g}$;

8: $n=n+1$;

9: if $n>Z$ then

10: $n=1$;

11: end if

12: end for

13: Calculate the $fitness$ of each individual $X(g)$ in $Q_{g}(H)$ using
(15)-(19);

14: /*Population evolution with CoDE */

15: while the maximum evaluation number is not achieved do

16: $Q_{g+1}=\varnothing$;

17: for $i=1\to n$ do

18: Generate trail vectors of $X(g)$ by the generation strategies and random
control parameters in CoDE;

19: Evaluate trail vectors and select them with $X(g)$ using (15)-(19);

20: $Q_{g+1}$ $\leftarrow$ Update charging amount for each CDU
$\\{CF_{1,1},\cdots,CF_{u,{MT}_{u}}\\}$ of $X(g)$;

21: end for

22: $g=g+1$;

23: end while

24: $S_{best}\leftarrow$ Select solutions with the best fitness from
$Q_{g}(H)$;

25: else

26: $S_{best}\leftarrow S_{CSO}(Z)$;

27: end if

With respect to the choice of the optimizer, it is hopeful for local search
operators to enhance the exploitation ability over continuous spaces. With
this information, we employ CoDE [30] as an exemplary optimizer to enhance the
exploitation ability in the local refinement stage. In CoDE, the control
parameters and trial vector generation strategies are randomly selected from a
candidate pool. The candidate pool of strategies and parameters can be
adjusted to suit different problem-solving scenarios. Algorithm 3 lists the
procedure of the local refinement using CoDE for repairing infeasible
solutions. In Lines 3-12, the initial population is composed of individuals
represented by continuous variables, while the discrete variables remain the
same as elite solutions from the global stage. The population for the local
refinement is evolved using CoDE operators to find feasible solutions, as
shown in Lines 15-23 of Algorithm 3.

Remark 2: It should be noted that the optimization scopes of the global search
stage and the local refinement stage are distinct. The global search stage
emphasizes the optimization of discrete variables, assisted by heuristic
rules, while the local refinement stage focuses on optimizing the continuous
variables of potential solutions. The integration of these two stages is
expected to enhance the performance of the proposed DSEA/HR in mixed decision
spaces. Furthermore, the selection of the optimizer for the local refinement
stage is not confined to a single option. An optimizer with a robust local
exploitation ability, serving the purpose of effective local refinement, would
facilitate the improvement of solutions towards more promising regions.

Figure 4: The flowchart of DSEA/HR for solving LSCOSPs.

###  III-G Framework of DSEA/HR

The flowchart of DSEA/HR is illustrated in Fig. 4, which consists of two
fundamental building modules: global search and local refinement. The initial
population is generated based on the heuristic rules. In the process of
evolution, two effective operators (i.e., CSO and CoDE) are adopted in global
search and local refinement, respectively. Individual fitness is evaluated by
the designed criterion before the selection. If the global search stage fails
to obtain a feasible solution, the population can be repaired by the local
refinement stage.

##  IV Case Study

###  IV-A Experimental Settings

The effectiveness of the proposed DSEA/HR on LSCOSPs is demonstrated by the
large-scale cases from a real-world marine-access refinery. The refinery has
two berths, twenty-one tanks (T1-T21), three CDUs (CDU1-CDU3), and four types
of residue tanks (R1-R4). The scheduling horizon is ten days (20 periods)
since the short-term scheduling of the refinery is carried out three times a
month. For the same refinery, the complexity of the scenarios is influenced by
certain conditions, such as the crude oil parcels to receive, initial tank
inventory, and the properties of crude oil. Three cases are tested in the case
study. Case 1 considers several new types of parcels that are different from
the initial component in each tank. The challenge of this case is that
unloading operations may significantly affect the stability of crude
components in tanks. Case 2 and Case 3 have a large number of low-quality
crudes. In particular, more than half of the crudes are high sulfur content or
heavy density in Case 3. Proper blending is crucial to satisfy the constraints
of properties and yields. The details of all the cases are listed in Section
S.II of the supplementary file.

Three groups of comparison experiments were designed as follows.

_Group 1:_ Comparison between DSEA/HR and state-of-the-art commercial MINLP
optimization solvers, including ANTIGONE [33], SCIP [34], and SHOT [35].

_Group 2:_ Comparison between DSEA/HR and state-of-the-art metaheuristics.
This includes DSEA/HR-PSO, where CSO is replaced by PSO [32] in the global
search stage of DSEA/HR, as well as CSO [29], the efficient recursive
differential grouping (ERDG) [36], the self-adaptive differential evolution
with multi-trajectory search (SaDE-MMTS) [37], and two metaheuristics for
crude oil scheduling mentioned in Section I, namely SAGA [22] and COSO-GA for
single objective [24].

_Group 3:_ Comparison between DSEA/HR and its two variants, including DSEA/HR
without the incorporation of heuristic rules (denoted by DSEA/HR-V1) and
DSEA/HR without local refinement (denoted by DSEA/HR-V2).

All metaheuristics were implemented with MATLAB 2020b, while the MINLP model
for LSCOSP was programmed and solved in GAMS version 33.1. In comparison with
deterministic approaches, the performance of algorithms was evaluated using
objective values and CPU times as metrics. The mean, standard deviation (std),
and feasible rate (FR) [38] were used as indicators to evaluate the optimality
and stability of experimental results compared with metaheuristics.

To ensure fair comparisons, each case was optimized using DSEA/HR or other
metaheuristics over 20 independent runs. The Wilcoxon’s rank-sum test at a
0.05 significance level was employed to statistically compare the algorithm
performance. In addition, all experiments were executed on Windows 10 with
Intel(R) Xeon(R) Gold 6256 CPU @ 3.60 GHz and 16.0 GB RAM.

###  IV-B Parameter Settings

####  IV-B1 Operators

The proposed DSEA/HR integrates the CSO in the global search stage and the
CoDE in the local refinement stage. To demonstrate the simplicity and
adaptability of our approach, the control parameters in CSO (i.e., $\varphi$)
and CoDE (i.e., scaling factor and crossover rate) are set based on their
original papers’ recommendations [29, 30].

####  IV-B2 Swarm and population size

The swarm size of CSO is set to 100 according to the recommendation in [29],
while the population size of CoDE requires determination, as setting details
are not provided in [30]. Table III presents a sensitivity analysis on
different population sizes ($PS$). The performance of DSEA/HR using
$PS\\!\\!=\\!\\!60$ outperforms than those using $PS\\!\\!=\\!\\!30$ and
$PS\\!\\!=\\!\\!90$ in three cases. Consequently, we set $PS$ as $60$ in this
paper.

TABLE III: Indicator Results Obtained by DSEA/HR that Adopts Different Setting
Values for $PS$ after 20 Runs.

Case | PS=30 | PS=60 | PS=90
---|---|---|---
mean$\pm$std | FR |  mean$\pm$std | FR |  mean$\pm$std | FR
Case 1 | $13.7500\pm 2.2913$ | $70\%$ |  $12.2500\pm 1.6819$* | $95\%$ | $14.5000\pm 2.7434$ | $80\%$
Case 2 | $14.9000\pm 4.3758$ | $60\%$ |  $11.3000\pm 1.5927$* | $90\%$ | $13.2000\pm 3.3182$ | $50\%$
Case 3 | $25.3000\pm 2.2501$ | $75\%$ |  $23.4500\pm 2.4382$* | $80\%$ | $26.1000\pm 2.7701$ | $60\%$

The best results are shown in boldface. ‘*’ indicates that the result is
significant better than the peer algorithm at a 0.05 level by the Wilcoxon’s
rank-sum test.

####  IV-B3 Number of function evaluations

The total number of fitness evaluations of the populations is used as the
termination criterion and set to a sufficiently large value to ensure
convergence. Specifically, the number of fitness evaluations for the two
search stages of DSEA/HR is set as $1\times 10^{5}$ and $3\times 10^{4}$,
respectively.

###  IV-C Comparisons with other state-of-the-art algorithms

In this subsection, we first verify the effectiveness of DSEA/HR with
traditional deterministic methods. Then we compare the performance of DSEA/HR
with metaheuristic algorithms tailored for large-scale optimization problems.

####  IV-C1 Comparisons with deterministic algorithms

The scale of variables in the LSCOSP is related to the following product of
several factors, i.e., the number of unloading crude types $\times$ tanks
$\times$ periods, the number of tanks $\times$ CDUs $\times$ periods, and the
number of CDUs $\times$ residue types $\times$ periods. Apparently, the scale
of the problem is considerable with the increase in the scheduling horizon. In
order to investigate the effectiveness of the proposed approach for large-
scale optimization, we extend each of the above cases to two examples with
different lengths of the horizon, in which one implements three-day
scheduling, while the other carries out the scheduling with ten days. Table IV
lists the problem scale of each example.

TABLE IV: Problems scale for MINLP formulation

Case | Example |  Parcels $\times$ Tanks $\times$ Crude types $\times$ CDUs $\times$ Residue types | Horizon (periods) | Total var. | Binary var. | Total cons. | Nonlinear cons.
---|---|---|---|---|---|---|---
Case 1 | Ex.1A | $1\times 19\times 14\times 3\times 4$ | 3 days (6) | 9948 | 538 | 19687 | 11172
Ex.1B | $6\times 19\times 14\times 3\times 4$ | 10 days (20) | 36960 | 3660 | 79942 | 40964
Case 2 | Ex.2A |  $1\times 21\times 13\times 3\times$ 4 | 3 days (6) | 10200 | 534 | 20062 | 11466
Ex.2B |  $4\times 21\times 13\times 3\times$ 4 | 10 days (20) | 35920 | 2740 | 75843 | 42042
Case 3 | Ex.3A | $2\times 21\times 11\times 3\times 4$ | 3 days (6) | 8928 | 546 | 17836 | 9702
Ex.3B | $8\times 21\times 11\times 3\times 4$ | 10 days (20) | 33760 | 3820 | 74412 | 35574

TABLE V: Comparison of MINLP Algorithms and DSEA/HR Results

Example | DSEA/HR | ANTIGONE | SCIP | SHOT
---|---|---|---|---
Objective value | CPU time (s) | Objective value | CPU time (s) | Objective value | CPU time (s) | Objective value | CPU time (s)
Ex.1A | 0 | 15 | 0 | 110 | 0 | 234 | $-$ | $-$
Ex.1B | 12 | 70 | $-$ | $-$ | $-$ | $-$ | $-$ | $-$
Ex.2A | 1 | 14 | 0 | 540 | 0 | 474 | 10 | 43
Ex.2B | 11 | 98 | $-$ | $-$ | $-$ | $-$ | $-$ | $-$
Ex.3A | 2 | 19 | 0 | 2454 | $-$ | $-$ | $-$ | $-$
Ex.3B | 23 | 286 | $-$ | $-$ | $-$ | $-$ | $-$ | $-$

The best results are shown in boldface. ‘$-$’ indicates that no feasible
solution was obtained within the allowable CPU time (3600 s).

The proposed DSEA/HR is compared with three commercial MINLP solvers (i.e.,
ANTIGONE [33], SCIP [34], and SHOT [35]) on the above examples. These solvers
integrate various state-of-the-art optimization technologies, such as
nonlinear branch-and-bound [39], extended supporting hyperplane [40], and
decomposition-based outer approximation [41].

Table V presents the numerical results, in which we set the cost coefficient
$\omega=1$ in the objective function for the convenience of calculation. For
scheduling problems within three days (i.e., Examples 1A, 2A, and 3A), DSEA/HR
and MINLP-based algorithms are able to find feasible solutions. Example 1A
considers a network consisting of 19 tanks, three CDUs, and four types of
residue tanks. In this example, initial inventory levels in crude tanks are
high, facilitating the charging operation. It can be observed that DSEA/HR and
two MINLP-based algorithms find the schedule without any changeover
operations, while DSEA/HR has better performance on computing time. This is
because MP-based methods require much time to deal with multiple equations,
even in a case with a simple initial situation. In contrast, with the guidance
of heuristic rules, the proposed EA-based approach can effortlessly search the
feasible region on the condition of sufficient initial crude inventory. For
Examples 2A and 3A, the number of tanks increases to 21, and more low-quality
crude is contained in tanks. In Example 2A, the properties of crudes in eight
tanks are high-sulfur or high-density, which exceeds the limits of the feed
properties allowed by CDUs. Similarly, the properties of crude in five tanks
are both high-sulfur and high-density in Example 3A. Apparently, the strict
quality properties constraints can lead to a small feasible region, which
affects the computing time and the accuracy of solutions. However, our
approach obtains the approximate optimal results in a very short time, as
shown in Examples 2A and 3A in Table V. It is worth mentioning that the CPU
time for DSEA/HR does not remarkably increase with the lower crude oil
properties.

For Examples 1B-3B, the length of the scheduling horizon goes up to ten days.
Except for more crude oil parcels needed to receive, the other conditions are
the same as Examples 1A-3A. As observed in Table V, the MINLP methods cannot
find any feasible solution for Examples 1B-3B within the CPU time limit. In
contrast, the proposed DSEA/HR identifies the feasible solution within an
acceptable computing time. Overall, though our approach is inferior to the MP-
based method in the optimality for general-scale problems, it is prospective
to provide a feasible solution for larger problems via the mechanisms of
heuristic rules and the repair strategy. In order to confirm this assertion,
we will discuss the performance of the two mechanisms in Section IV-D.

####  IV-C2 Comparisons with metaheuristic algorithms

The proposed DSEA/HR is compared with six state-of-the-art metaheuristics
(i.e., DSEA/HR-PSO, CSO, ERDG, SaDE-MMTS, SAGA, and COSO-GA) on LSCOSPs. Among
these compared algorithms, DSEA/HR-PSO is a variant of DSEA/HR that utilizes
PSO for the global search of the original DSEA/HR. CSO is a simple and
effective large-scale optimization algorithm that has also been adopted as the
optimizer of the proposed DSEA/HR in the global search stage. ERDG, designed
in the cooperative co-evolution framework, has been proposed to solve large-
scale optimization problems with an efficient variable grouping strategy.
SaDE-MMTS is a modification of the original SaDE, satisfying both global and
local search requirements and demonstrating remarkable performance on the
high-dimensional test functions up to 1000 dimensions. SAGA and COSO-GA have
both been developed to address crude oil scheduling problems: the former
provides a graph-based structure for encoding, while the latter applies GA
only to optimize integer variables.

TABLE VI: Indicator Results Obtained by DSEA/HR and State-of-the-Art EAs after
20 Runs.

Algorithm | Case 1 | Case 2 | Case 3
---|---|---|---
| mean$\pm$std | FR | mean$\pm$std | FR | mean$\pm$std | FR
DSEA/HR |  $12.2500\pm 1.6819$* | $95\%$ |  $11.3000\pm 1.5927$* | $90\%$ |  $23.4500\pm 2.4382$* | $80\%$
DSEA/HR-PSO | $12.8000\pm 1.5079$ | $95\%$ | $13.7500\pm 1.6162$ | $75\%$ | $24.5500\pm 2.5021$ | $70\%$
CSO | $14.3500\pm 2.7582$ | $20\%$ | $15.5000\pm 2.6852$ | $15\%$ | $26.1000\pm 2.2919$ | $5\%$
ERDG | $13.7000\pm 2.5152$ | $10\%$ | $15.8000\pm 2.2384$ | $10\%$ | $28.1000\pm 2.1250$ | $0\%$
SaDE-MMTS | $15.7000\pm 1.6575$ | $25\%$ | $14.4500\pm 2.1879$ | $20\%$ | $26.6500\pm 1.4318$ | $10\%$
SAGA | $14.1500\pm 2.5397$ | $45\%$ | $14.1000\pm 2.1002$ | $15\%$ | $26.4500\pm 1.9595$ | $10\%$
COSO-GA | $15.2500\pm 2.5726$ | $50\%$ | $13.9000\pm 1.6827$ | $20\%$ | $25.5000\pm 1.9057$ | $20\%$

The best results are shown in boldface. ‘*’ indicates that the result is
significant better than the peer algorithm at a 0.05 level by the Wilcoxon’s
rank-sum test.

From the results summarized in Table VI, it can be observed that DSEA/HR
significantly outperforms the other comparison algorithms in term of both
optimality and stability on all cases. Notably, three competitors (i.e., CSO,
ERDG, and SaDE-MMTS) struggled to find feasible solutions, and ERDG even
failed to find any feasible solution within 20 runs for the cases with the
relatively small feasible regions. These outcomes emphasize that relying on a
single optimizer is insufficient for real-world scheduling problems. The
results obtained from SAGA and COSO-GA show that these two approaches,
specifically designed for crude oil scheduling problems, do not exhibit
significant performance when applied to solving large-scale instances.
Additionally, the variant DSEA/HR-PSO is slightly worse than DSEA/HR in terms
of optimality, while DSEA/HR-PSO is comparable to DSEA/HR in terms of
stability. This comparison indicates that while the superior performance of
DSEA/HR is mainly due to the dual-stage framework, the incorporation of
powerful optimizers still contributes to identifying superior solutions.

The above experimental results demonstrate the significance of the dual-stage
search framework in the proposed DSEA/HR, thereby verifying the main
motivation of this article. Additionally, there exists an opportunity to
further enhance the efficacy of the algorithm by incorporating several
strategies. In the subsequent section, we will discuss the roles of kernel
strategies designed for each stage of DSEA/HR, namely heuristic rules and the
repair strategy.

###  IV-D Effectiveness of the Components in DSEA/HR

To investigate the effectiveness of two key components of DSEA/HR, namely the
heuristic rules and the repair strategy, we adapt the original DSEA/HR into
two variants. These variants are named DSEA/HR-V1 and DSEA/HR-V2. In
DSEA/HR-V1, the heuristic rules for population initialization in the DSEA/HR
algorithm are removed. As for DSEA/HR-V2, the repair strategy in the DSEA/HR
algorithm is discarded.

TABLE VII: Indicator Results Obtained by DSEA/HR and its Two Variants
Algorithms after 20 Runs.

Algorithm | Case 1 | Case 2 | Case 3
---|---|---|---
| mean$\pm$std | FR | mean$\pm$std | FR | mean$\pm$std | FR
DSEA/HR |  $12.2500\pm 1.6819$* | $95\%$ |  $11.3000\pm 1.5927$* | $90\%$ |  $23.4500\pm 2.4382$* | $80\%$
DSEA/HR-V1 | $13.6500\pm 1.8994$ | $50\%$ | $14.9500\pm 2.0641$ | $40\%$ | $26.2000\pm 1.8525$ | $30\%$
DSEA/HR-V2 | $13.1000\pm 1.6512$ | $65\%$ | $12.5500\pm 1.5381$ | $20\%$ | $25.2500\pm 1.6504$ | $15\%$

The best results are shown in boldface. ‘*’ indicates that the result is
significant better than the peer algorithm at a 0.05 level by the Wilcoxon’s
rank-sum test.

Figure 5: Box plots for the objective values by DSEA/HR, DSEA/HR-V1, and
DSEA/HR-V2 in Cases 1-3. (a) Case 1. (b) Case 2. (c) Case 3.

Table VII shows the statistical results obtained by DSEA/HR and the two
variants. It can be observed that DSEA/HR significantly outperforms its two
variants in all cases. Moreover, in the comparison of the two variants, due to
the high inventory levels and the proper crude quality in Case 1, the guidance
of heuristic rules is dominant. Thus, the performance of DSEA/HR-V2 with
heuristic rules is better than DSEA/HR-V1. As discussed in Section IV-C, more
low-quality crude in Cases 2 and 3 are considered, which leads to a small
feasible region. Thus, repairing infeasible solutions becomes indispensable.
As expected, Table VII demonstrates that DSEA/HR-V1 slightly improves
performance compared to DSEA/HR-V2 through the assistance of the repair
strategy. Similar observations are noted in the box plots in Fig. 5. The
heuristic rules contribute to assisting DSEA/HR in the search for better
solutions in case adequate computational resources are available. On the other
hand, the repair strategy improves the computing performance of DSEA/HR in
cases with smaller feasible regions. To gain insight into how these two
components respectively work well in addressing the LSCOSPs, a detailed
analysis of scheduling results is presented in the following section.

####  IV-D1 Effectiveness of the proposed heuristic rules in global search

As mentioned earlier, the proposed DSEA/HR comprises two heuristic rules: one
for assigning discharging tanks (denoted as Rule I) and the other for
selecting charging tanks for each CDU (denoted as Rule II). To facilitate a
better understanding and explanation of their effectiveness, we present a
feasible schedule solved by the proposed DSEA/HR for Case 2 in Fig. 6. The
horizontal axis represents the scheduling time, while the vertical axis
indicates crude oil tanks.

Figure 6: Gantt chart of Case 2 obtained by DSEA/HR. The blue numbers indicate
instances of reasonable unloading operations. Figure 7: Gantt chart of Case 2
obtained by DSEA/HR-V1. The red numbers indicate instances of unreasonable
unloading operations. Figure 8: Gantt chart of Case 3 obtained by DSEA/HR and
DSEA/HR-V2. The red numbers indicate the results obtained by DSEA/HR-V2, while
the blue numbers indicate the results obtained by DSEA/HR.

Regarding the effectiveness of Rule I on unloading operations, it can be
observed from Fig. 6 that the tanks (i.e., T09, T12, and T16) continuously
charge to CDU1 between periods 10 and 20 after receiving crude C11 (indicated
by blue numbers in Fig. 6). In contrast, Fig. 7, obtained by DSEA/HR-V1,
depicts a schedule that violates the property constraint due to the absence of
heuristic rules for unloading operations. Specifically, the feed densities of
CDU1 during periods 16-20 (i.e., 0.889, 0.9, 0.907, 0.909, 0.909 $g/cm^{3}$)
exceed its upper limit (0.885 $g/cm^{3}$). This occurs because the low-density
crude (i.e., C11) that satisfies the processing requirement is unloaded into
tanks already containing other incompatible crudes (i.e., C2, C4, and C9) for
CDU1. Apparently, the random unloading operations without the heuristic
guidance contribute significantly to infeasibility and increased changeovers.

Furthermore, the performance of heuristic rules is also evident in charging
operations. Fig. 7 shows an infeasible schedule solved by DSEA/HR-V1. A
comparison between Fig. 6 and Fig. 7 reveals that DSEA/HR achieves fewer
changeovers than DSEA/HR-V1. This improved performance can be attributed to
the heuristic rule designed for charging operations, which enhances the
diversity of the tanks charging to a CDU. In contrast, due to the absence of
this rule, several tanks (i.e., T14, T18, and T21) often charge a large amount
of crude to two CDUs simultaneously. Although it is allowed in this problem,
the drawback is that rapid consumption of inventory results in frequent
changeovers between tanks and CDUs.

Based on the above comparisons, it can be concluded the proposed two heuristic
rules are conducive to the performance improvement of the proposed DSEA/HR.
The mechanism of heuristic rules, derived from empirical knowledge, can also
be extended to more real-world instances.

####  IV-D2 Effectiveness of the repair strategy in local refinement

We have implemented a repair strategy during the local refinement stage to
enhance the algorithm’s search for feasible solutions. To elaborate on this,
we have used Case 3 as an example, which poses a significant challenge due to
the involvement of a large amount of low-quality crude types.

Our repair strategy, focusing on modifying the flow rate charging to CDUs,
effectively solves this issue. Fig. 8 demonstrates this with a schedule
involving all decisions of charging tanks, which is obtained by DSEA/HR-V2 and
then modified by the repair strategy. From the original results without repair
in Fig. 8, we can see that the sulfur content of the feed charging to CDU1
exceeds the upper limit during periods 1-6. In addition, the yield of naphtha
produced by CDU2 falls below the lower bound during periods 12-20. These
violations, however, can be completely eliminated through the local adjustment
of the flow rate. Specifically, with the help of the repair strategy in the
local refinement stage, the flow rate tied to decisions causing violations is
further optimized in the continuous decision variable space. Following the
flow rate optimization, we can observe from Fig. 8 that this schedule can be
modified into a feasible solution that satisfies all constraints. As a result,
implementing the repair strategy in the local refinement stage enhances
DSEA/HR’s capability to search for feasible regions, particularly in instances
that involve multiple low-quality crudes.

Based on the above discussions, it can be concluded that both the heuristic
rules and the repair strategy improve the optimization ability of the proposed
DSEA/HR in discrete and continuous variable spaces, respectively. The
integration of these strategies further enhances the solving capability of the
EA-based algorithm, enabling it to find feasible solutions within an
acceptable time.

##  V Conclusions

It is challenging to find the feasible solution of the large-scale crude oil
scheduling problem, especially in a reasonable time. This paper investigated
the LSCOSPs from a large marine-access refinery in practice. Firstly, the
LSCOSPs accommodating the practical operating features were modeled. Following
the model, we developed a novel EA-based framework called DSEA/HR to address
the LSCOSPs. Based on the empirical knowledge, we designed two heuristic rules
for population initialization. One heuristic rule guides the decisions of
unloading operations, while the other drives the decisions of charging
operations. Moreover, to improve the efficacy of the algorithm in solving
LSCOSPs with a small feasible region, we proposed a repair strategy to
intensify the search for the continuous decision space.

In the experiments, we have compared the proposed algorithm with the state-of-
the-art MINLP methods and metaheuristic algorithms on a number of LSCOSP
instances. The comparative experimental results have demonstrated that DSEA/HR
is superior to the competitors and is able to converge to the feasible region
for all instances in an acceptable time. Meanwhile, the efficacy of the
special designs is also further verified, including heuristic rules and the
repair strategy. To the best of our knowledge, it is the first attempt to
solve such large-scale crude oil scheduling problems that are very close to
the operational features of an actual refinery. The applicability of the
proposed method has also been verified in three practical cases.

A few issues remain to be further investigated in future work. First, while
problem-specific knowledge is conducive to improving the efficacy of the
algorithm, it cannot be directly applied to other scheduling problems. Second,
the experimental results indicate that the algorithm does not guarantee the
discovery of a feasible solution in each run. The feasible rate of algorithms
tends to decrease as the complexity of the cases increases. To address these
concerns, future research will focus on automatically extracting problem-
specific knowledge by integrating emerging knowledge transfer techniques and
incorporating efficient search strategies to enhance the stability of
scheduling approaches.

## References

  * [1] N. K. Shah, Z. Li, and M. G. Ierapetritou, “Petroleum refining operations: Key issues, advances, and opportunities,” _Ind. Eng. Chem. Res._ , vol. 50, no. 3, pp. 1161–1170, 2011.
  * [2] F. Li, M. Yang, W. Du, and X. Dai, “Development and challenges of planning and scheduling for petroleum and petrochemical production,” _Front. Eng. Manag._ , vol. 7, no. 3, pp. 373–383, 2020.
  * [3] G. K. Saharidis, M. Minoux, and Y. Dallery, “Scheduling of loading and unloading of crude oil in a refinery using event-based discrete time formulation,” _Comput. Chem. Eng._ , vol. 33, no. 8, pp. 1413–1426, 2009\.
  * [4] Y. Zhao, N. Wu, Z. Li, and T. Qu, “A novel solution approach to a priority-slot-based continuous-time mixed integer nonlinear programming formulation for a crude-oil scheduling problem,” _Ind. Eng. Chem. Res._ , vol. 55, no. 41, pp. 10 955–10 967, 2016.
  * [5] P. C. P. Reddy, I. A. Karimi, and R. Srinivasan, “A new continuous-time formulation for scheduling crude oil operations,” _Chem. Eng. Sci._ , vol. 59, no. 6, pp. 1325–1341, 2004.
  * [6] J. Cerdá, P. C. Pautasso, and D. C. Cafaro, “Efficient approach for scheduling crude oil operations in marine-access refineries,” _Ind. Eng. Chem. Res._ , vol. 54, no. 33, pp. 8219–8238, 2015.
  * [7] J. Li, W. Li, I. A. Karimi, and R. Srinivasan, “Improving the robustness and efficiency of crude scheduling algorithms,” _AIChE J._ , vol. 53, no. 10, pp. 2659–2680, 2007.
  * [8] S. Mouret, I. E. Grossmann, and P. Pestiaux, “A new Lagrangian decomposition approach applied to the integration of refinery planning and crude-oil scheduling,” _Comput. Chem. Eng._ , vol. 35, no. 12, pp. 2750–2766, 2011\.
  * [9] L. S. Assis, E. Camponogara, B. C. Menezes, and I. E. Grossmann, “An MINLP formulation for integrating the operational management of crude oil supply,” _Comput. Chem. Eng._ , vol. 123, pp. 110–125, 2019.
  * [10] L. S. Assis, E. Camponogara, and I. E. Grossmann, “A MILP-based clustering strategy for integrating the operational management of crude oil supply,” _Comput. Chem. Eng._ , vol. 145, p. 107161, 2021.
  * [11] S. Yadav and M. A. Shaik, “Short-term scheduling of refinery crude oil operations,” _Ind. Eng. Chem. Res._ , vol. 51, no. 27, pp. 9287–9299, 2012\.
  * [12] P. C. P. Reddy, I. A. Karimi, and R. Srinivasan, “Novel solution approach for optimizing crude oil operations,” _AIChE J._ , vol. 50, no. 6, pp. 1177–1197, 2004.
  * [13] D. Panda and M. Ramteke, “Reactive scheduling of crude oil using structure adapted genetic algorithm under multiple uncertainties,” _Comput. Chem. Eng._ , vol. 116, pp. 333–351, 2018.
  * [14] S. Das and P. N. Suganthan, “Differential evolution: A survey of the state-of-the-art,” _IEEE Trans. Evol. Comput._ , vol. 15, no. 1, pp. 4–31, 2011.
  * [15] Y. Tian, L. Si, X. Zhang, R. Cheng, C. He, K. C. Tan, and Y. Jin, “Evolutionary large-scale multi-objective optimization: A survey,” _ACM Comput. Surv._ , vol. 54, no. 8, 2021.
  * [16] M. Li and G. Wang, “A review of green shop scheduling problem,” _Inf. Sci._ , vol. 589, pp. 478–496, 2022.
  * [17] K. Gao, Z. Cao, L. Zhang, Z. Chen, Y. Han, and Q. Pan, “A review on swarm intelligence and evolutionary algorithms for solving flexible job shop scheduling problems,” _IEEE/CAA J. Automatica Sinica_ , vol. 6, no. 4, pp. 904–916, 2019.
  * [18] J. Branke, S. Nguyen, C. W. Pickardt, and M. Zhang, “Automated design of production scheduling heuristics: A review,” _IEEE Trans. Evol. Comput._ , vol. 20, no. 1, pp. 110–124, 2016.
  * [19] Z. Pan, D. Lei, and L. Wang, “A bi-population evolutionary algorithm with feedback for energy-efficient fuzzy flexible job shop scheduling,” _IEEE Trans. Syst., Man, Cybern., Syst._ , vol. 52, no. 8, pp. 5295–5307, 2022.
  * [20] J. Wang and L. Wang, “A bi-population cooperative memetic algorithm for distributed hybrid flow-shop scheduling,” _IEEE Trans. Emerg. Topics Comput. Intell._ , vol. 5, no. 6, pp. 947–961, 2021.
  * [21] W. Du, W. Zhong, Y. Tang, W. Du, and Y. Jin, “High-dimensional robust multi-objective optimization for order scheduling: A decision variable classification approach,” _IEEE Trans. Ind. Informat._ , vol. 15, no. 1, pp. 293–304, 2019.
  * [22] M. Ramteke and R. Srinivasan, “Large-scale refinery crude oil scheduling by integrating graph representation and genetic algorithm,” _Ind. Eng. Chem. Res._ , vol. 51, no. 14, pp. 5256–5272, 2012.
  * [23] D. Panda and M. Ramteke, “Preventive crude oil scheduling under demand uncertainty using structure adapted genetic algorithm,” _Appl. Energy_ , vol. 235, pp. 68–82, 2019.
  * [24] Y. Hou, N. Wu, M. Zhou, and Z. Li, “Pareto-optimization for scheduling of crude oil operations in refinery via genetic algorithm,” _IEEE Trans. Syst., Man, Cybern., Syst._ , vol. 47, no. 3, pp. 517–530, 2017.
  * [25] Y. Hou, N. Wu, Z. Li, Y. Zhang, T. Qu, and Q. Zhu, “Many-objective optimization for scheduling of crude oil operations based on NSGA-III with consideration of energy efficiency,” _Swarm Evol. Comput._ , vol. 57, p. 100714, 2020.
  * [26] Y. Jin, Ed., _Knowledge Incorporation in Evolutionary Computation_. Heidelberg, Germany: Springer, 2005.
  * [27] X. Dai, L. Zhao, Z. Li, W. Du, W. Zhong, R. He, and F. Qian, “A data-driven approach for crude oil scheduling optimization under product yield uncertainty,” _Chem. Eng. Sci._ , vol. 246, p. 116971, 2021.
  * [28] J. M. Wassick and J. Ferrio, “Extending the resource task network for industrial applications,” _Comput. Chem. Eng._ , vol. 35, no. 10, pp. 2124–2140, 2011.
  * [29] R. Cheng and Y. Jin, “A competitive swarm optimizer for large scale optimization,” _IEEE Trans. Cybern._ , vol. 45, no. 2, pp. 191–204, 2015\.
  * [30] Y. Wang, Z. Cai, and Q. Zhang, “Differential evolution with composite trial vector generation strategies and control parameters,” _IEEE Trans. Evol. Comput._ , vol. 15, no. 1, pp. 55–66, 2011.
  * [31] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: NSGA-II,” _IEEE Trans. Evol. Comput._ , vol. 6, no. 2, pp. 182–197, 2002.
  * [32] J. Kennedy and R. Eberhart, “Particle swarm optimization,” in _Proc. ICNN’95–Int. Conf. Neural Networks_ , vol. 4, 1995, pp. 1942–1948.
  * [33] R. Misener and C. A. Floudas, “ANTIGONE: algorithms for continuous/integer global optimization of nonlinear equations,” _J. Glob. Optim._ , vol. 59, no. 2, pp. 503–526, 2014.
  * [34] T. Achterberg, “SCIP: Solving constraint integer programs,” _Math. Program. Comput._ , vol. 1, no. 1, pp. 1–41, 2009.
  * [35] A. Lundell, J. Kronqvist, and T. Westerlund, “The supporting hyperplane optimization toolkit for convex MINLP,” _J. Glob. Optim._ , vol. 84, no. 1, pp. 1–41, 2022.
  * [36] M. Yang, A. Zhou, C. Li, and X. Yao, “An efficient recursive differential grouping for large-scale continuous problems,” _IEEE Trans. Evol. Comput._ , vol. 25, no. 1, pp. 159–171, 2021.
  * [37] S. Zhao, P. N. Suganthan, and S. Das, “Self-adaptive differential evolution with multi-trajectory search for large-scale optimization,” _Soft Comput._ , vol. 15, no. 11, pp. 2175–2185, 2011.
  * [38] Z. Liu and Y. Wang, “Handling constrained multiobjective optimization problems with constraints in both the decision and objective spaces,” _IEEE Trans. Evol. Comput._ , vol. 23, no. 5, pp. 870–884, 2019.
  * [39] J. N. Hooker, “Planning and scheduling by logic-based benders decomposition,” _Oper. Res._ , vol. 55, no. 3, pp. 588–602, 2007.
  * [40] J. Kronqvist, A. Lundell, and T. Westerlund, “The extended supporting hyperplane algorithm for convex mixed-integer nonlinear programming,” _J. Glob. Optim._ , vol. 64, no. 2, pp. 249–272, 2016.
  * [41] P. Muts, I. Nowak, and E. M. T. Hendrix, “The decomposition-based outer approximation algorithm for convex mixed-integer nonlinear programming,” _J. Glob. Optim._ , vol. 77, no. 1, pp. 75–96, 2020.