From Data to Action: Accelerating Refinery Optimization with AI

Highlights From Data to Action: Accelerating Refinery Optimization with AI

arXiv:2605.15085v1 [stat.ML] 14 May 2026

Dániel Pfeifer, Ábrahám Papp, Tibor Bernáth, Tamás Zoltán Varga, Márk Czifra, Botond Szilágyi, Edith Alice Kovács • Accelerating Big Data Refinery Optimization and decision making • Understanding business opportunities by analyzing types of anomalies • New multivariate Anomaly Detection methodologies • Enhancing LP-based decision making with Anomaly Detection tools

From Data to Action: Accelerating Refinery Optimization with AI Dániel Pfeifera , Ábrahám Pappa , Tibor Bernáthb , Tamás Zoltán Vargab , Márk Czifrab , Botond Szilágyia , Edith Alice Kovácsa a

Budapest University of Technology and Economics, Műegyetem rkp. 3., Budapest, 1111, Hungary b MOL Group, Dombóvári út 28., Budapest, 1117, Hungary

Abstract Nowadays refinery optimization utilizes sheer amounts of data, which can be handled with modern Linear Programming (LP) software, but the interpreting and applying the results remains challenging. Large petrochemical companies use massive models, with hundreds of thousands of input matrix elements. The LP solution is mathematically correct, but simplifications are made in the model, and data supply errors may occur. Therefore, further insight is needed to trust the results. The LP solver does not have a memory, so additional understanding could be gained by analyzing historical data and comparing it to the current plan. As such, machine learning approaches were suggested to support decision making based on the LP solution. Among these, Anomaly Detection tools are proposed to be used in tandem with the LP output. A transformed version of the popular ECOD methodology is applied. New methods are proposed to handle high-dimensional data: choosing the most informative pairs. Then, this is used alongside two 2D Anomaly Detection algorithms, revealing several business opportunities and data supply errors in the MOL refinery scheduling and planning architecture. Keywords: Anomaly Detection, Big Data, Linear Programming, Marginal Values, Business Opportunities 1. Introduction Modern technology applications generate vast amounts of data. Discovering hidden patterns within existing data and leveraging them for business insights is becoming increasingly important.

In general, Big Data refers to large and complex datasets that are often used for predictive analytics. It is a term that relates to massive, heterogeneous, and often unstructured digital data that are difficult to handle with traditional data management tools and techniques [1]. The data analysis-based technologies that various organizations and companies have recently implemented are called Business intelligence. In M. Hamzehi et al. [2], the data set used is a data set related to the sales system of a pharmaceutical company, which enters the data into a data warehouse online. The paper aims to provide an efficient model for optimizing the products sales system in a pharmaceutical company using clustering methods based on machine learning indicators and algorithms. Rath et al. [3] introduced a structure to build a Business Intelligence framework. They show that machine learning approaches play an important role in the tasks of Business Intelligence in large commercial organizations. Another type of investigation is presented in a recent paper by F. Ridzuan et al. [4] related to the problem of outliers. It includes recommendations to improve the quality of data and data collection systems are provided by finding several factors that greatly influence the presence of outliers. Machine learning approaches are intensively applied to address real-world problems in the fourth industrial revolution. These kinds of approaches are now an integral part of the operations of most oil and gas companies, which they use to translate datasets, including large volumes of real-time information into actionable insights. A good overview of these methods is given in [5] and [6]. Aspects related to big data in this field are discussed in [7]. A clear literature overview of the methodologies for addressing scheduling, planning, and supply chain management of oil refinery operations can be found in [8]. A paper by Grossmann [9] establishes enterprise-wide optimization as the framework for integrating supply, manufacturing, and distribution decisions across multiple levels and sites. Another perspective on surveying the history and current state of AI and machine learning applications in chemical engineering can be found in [10]. The research presented in this paper is motivated by the problems which occur during the application of an efficient optimization solver in a Big Data environment. Linear Programming (LP) models have been a staple in Refinery Planning and Optimization since the late 20th century. These large mathematical models take hundreds of thousands of inputs and produce hundreds of thousands of outputs, to find the most economical way to operate one or more 2

refineries, petrochemical plants, their feedstock supply, market demands all while considering every relevant logistic or operational constraint. In the past decades, these models have become necessary, as finding the optimum and maximizing profits have become a must among shrinking margins, rising energy prices, volatile external environment and ever stricter regulatory pressure on the industry. An even more developed idea in a recent study to apply an AI safety system that operates in parallel with operational AI, continuously supervising system safety behavior can be found in [11]. MOL Group is a Central European integrated oil company, that has also been using an LP optimizer, AspenTech PIMS (Process Industry Modelling System) [12] for more than 30 years to optimize their monthly plans and evaluate the refinery economics. During this time a wide range of supporting tools were developed internally with the sole purpose of helping to understand the result of the LP optimization. Due to the sheer amount of data and the complex technological and optimization knowledge needed to interpret the results it remains a daunting task. While LP gives a mathematically optimal solution it is the optimization experts’ responsibility together with refinery scheduling, logistics and other disciplines to make sure that these plans make sense not just mathematically but economically and operationally as well. Analyzing the entirety of the input & output data is not possible as of now, since the size and complexity of the problem is in the world of big data. So far key process and economic indicators were defined and monitored, however everchanging external factors such as Covid, NatGas price volatility and other side effects of the RUS-UA war can fundamentally change main drivers of optimization. To tackle these obstacles, AI solutions may be used to help our understanding, find business opportunities and maximize margins. Anomaly Detection tools are being implemented that compare the inputs & outputs of new plans to old ones which have already been checked by all stakeholders and accepted by the management. This study will focus on three new Anomaly Detection methodologies that have been fine-tuned by applying them to real data from monthly plans and illustrate how LP and Anomaly Detection can go hand in hand to enhance each other’s, and the company’s performance. The remainder of the paper is the following. In Section 2, we introduce the Linear Programming methodology used to solve Refinery Optimization problems. We detail how Marginal and Incremental Values help experts find business opportunities, and introduce the LP software used. In Section 3, 3

we introduce Anomaly Detection as a tool for finding such business opportunities. Through several real life examples, we show kinds of anomalies that may occur. Finally, we introduce the input data that the Anomaly Detection will use. In Section 4, we detail the exact mathematical methodology of Anomaly Detection, and introduce two 2D algorithms and one 1D modified algorithm for detecting values or value pairs that do not conform the historical distribution. Additionally, we introduce a pair selection methodology that can vastly reduce the number of possible pairs to look through. In Section 5, we illustrate some real life 1D and 2D anomalies we have found with this methodology. Finally, Section 6 summarizes the paper. 2. Refinery Planning & Optimization 2.1. Linear Programming methodology A central methodology of oil refinery planning is Linear Programming. An early formulation of an LP/NLP refinery planning model, covering processing units, intermediate streams, and blending relations with real industrial case studies can be found in [13]. Here, we will introduce the primal and dual Linear Programming problems and show how they are related to each other. Let A ∈ Rm×n be the technological matrix obtained from data supply containing known parameters (e.g. prices and costs, capacities, blending proportions etc.), x ∈ Rn the decision variables of the primal problem (e.g. amounts bought and sold, amounts produced, amounts transferred, etc.), b ∈ Rm the constraint limits of the primal problem, (e.g. upper bounds of the available resources), and c ∈ Rn the weights of the objective function. The primal linear programming problem has the following form: (P)

maxn cT x x∈R

subject to Ax ≤ b, x≥0

(1)

In the primal problem we maximize a linear objective function, the solutions (decision variables) are the elements of x vector. Based on the primal problem (Formula 1), we can formulate the dual problem in the following way: (D)

bT y

min

y∈Rm

subject to AT y ≥ c, y≥0 4

(2)

The decision variables of the dual problem are contained in y ∈ Rm . Its components yi are the dual decision variables, one for each primal constraint. By the strong duality theorem [14, p. 146-155] if both problems have solutions, then the maximum of the primal objective and the minimum of dual objective functions are equal. In the refinery case, there are more unknowns than equations, therefore n > m. Marginal Values are an important concept in LP. Their values are generated in each LP run. They can be extracted for all constraints. Non-zero Marginal values appear in the cases where variables are non basis. Marginal Values indicate what would be the increment of the objective (gross profit) function in case of an infinitesimal change on the given bound. The more the given bound constrains the profit generation the higher the Marginal Value will be. If the given bound is not constraining, the profit generation there is no (or 0) MV on that constraint. The Marginal Value of a row of constraints aTi x ≤ bi is ∂ max{cT x|Ax ≤ b} ∂bi From the dual problem, the Marginal Value of a column can also be defined. Since the dual problem has a transposed A matrix, each constraint in the dual problem can be interpreted as a column of the original problem. So the Marginal Value of a column j is Pi =

Dj =

∂ min{bT y|AT y ≥ c} ∂cj

A constraint aTi x ≤ bi is called basis, if ai x̂ = bi holds where x̂ is the solution of the LP. 2.2. Refinery Economics LP optimization is particularly valuable in the complex environment of refineries and petrochemical plants, where numerous variables and constraints must be managed simultaneously. Marginal Values (MV) Marginal Values are rates of change of the Objective Function with a change in the activity of a variable. As partial derivatives of the Objective Function – they are only valid at the increments around the optimal activity. 5

Dj : Applies to column activities (variables), and represents rate of change in objective function as the bound (min or max) of a variable is increased. Pi : Applies to rows (equations), and represents the rate of change in the objective function as the constraint (min or max) of a row is decreased.

Figure 1: Example: How much we

would earn (USD/t) on selling the Incremental Values (IV) Incremental Values are internal values of 40001st ton? Marginal Values are a material within the planning premises. For only valid for the increments around the original 40 000 tons, while BEV any active material and utility in the LP an considers average for the next 40 kt. incremental value is defined as well. These refer to the substitution value of the infinitesimal amount of the given product. Dj of a Basis Column: All basis columns will have a Dj equal to zero because the cost of producing the next increment and its incremental product value balance each other and the slope is zero. If a column has an entry (price or cost) in the objective function row, the Dj will be the difference between this entry and the incremental value of the column in a solution. Break-even Value: A break-even activity is a value which does not alter a refiner’s economics. In other words, they do not lose (nor make) money in the process. The break-even value of an activity is the unit price (or cost) of the activity that preserves the original objective function of the refinery LP model before the specified changes have taken place. The break-even position provides a benchmark for negotiating a discount or premium of a transaction.

2.3. Detecting Business Opportunities with Incremental and Marginal Values Marginal values (MV) are a key concept in LP, and support understanding the optimization drivers and result in a complex system. These values are generated in each LP run. Incremental Values are obtained from Marginal Values. They can only be calculated for variables xi for which ci ̸= 0 (or in other words, variables that appear in the objective function). In this case, these variables represent some physically ready material that can be sold, and, as such, we are interested in how its sell price is assembled. Incremental Value is then equal to the buy 6

costs of this product, plus its production costs, transfer costs, storage costs, etc., minus its Marginal Value. At every stage of production, this represents how valuable the given product is inside the refinery.

Figure 2: Simple representation of Marginal and Incremental Value dynamics. Any transformation in the LP (transfer between sites, inventorying between periods or a processing unit) can be described with the same methodology. Incremental and Marginal values are connected and can be calculated from each other.

Marginal and Incremental Values are real treasures for an optimizer, and most of the time they are the key for understanding complex optimization dynamics. Therefore the MOL Group Optimization team developed multiple reports and concepts to understand the LP results through MVs: • Reports for ranked MVs to capacities, sales/purchases, other constraining variables are especially useful to indicate the most important bottlenecks and optimization drivers, often revealing model stability and improper settings. • Graphs showing values of same materials across multiple periods and sites – explaining market sharing, transfer and inventory dynamics. • A simple diagrams explaining MV dynamics for transfers, inventory and unit utilization. • Graphical reports of unit feed and product MVs revealing bottleneck products and asymmetries in operation. • Multi-case (direction) runs specifically focusing on MV dynamics of a given material and downstream inflection points. Although these reports and developments already give strong support for result analysis, the dashboards are only created from a fraction of the LP 7

Figure 3: A product value diagram for one of MOL Group’s steam crackers, during the times of the historically high natural gas prices. In these months the historically lowest value Fuel gas had been the most valuable product of the Steam crackers. Cracker operation differences among the periods can be understood easier by observing Ethylene value differences among the periods.

results; and all of them focus on a single LP run, with no memory and not using relationships of variables from historical optimization results. This necessitates the introduction of a more sophisticated method that does have a memory and can use and connect historical information. 2.4. Linear Programming software - AspenTech PIMS It is a flexible and easy to use productivity tool for economic planning in the process industries. PIMS (Process Industry Modelling System) [12] employs Linear Programming technique to optimize the operation and economics of refineries, petrochemical plants, or other industry facilities. The A matrix representing a PIMS model is a collection of linear equations. Each of these equations represents one aspect of the refinery or petrochemical plant operation, logistics or commercial constraint. It may be used for a wide variety of short term and strategic planning purposes, such as the evaluation of alternative feedstocks, optimization of product slates, evaluation of grass root opportunities or expansions, and many others. Although PIMS is called a linear programming system, it has many nonlinear elements. 8

Figure 4: A product value diagram from multiple LP runs – in this case a crude evaluation. These graphs can reveal how different straight run products’ (direct products of crude distillation units) values change in case of different crudes processed, showing the product shortage (when high value is observed) or product surplus (low value observed – in extreme cases even negative ones)

Nowadays, the use of PIMS is common in the oil industry, as about 75% of the world’s refineries and about 60% of petrochemical companies use it for production planning and scheduling [12]. MOL Group has been using the PIMS software, by AspenTech since 1993. 3. Goal of Anomaly Detection 3.1. Real life challenges of refinery optimization MOL Group is an integrated oil and gas company, operating three refineries and two petrochemical sites in Central Europe in relative proximity (∼ 500km), with the majority of the sites being landlocked. This makes logistics and inventory management key optimization drivers, which require multi-site, multi-period (3-4 month) models. For this, MOL Group is using one of the largest and most complex XPIMS models in the world. The current model matrix consists of ∼ 45 000 rows, ∼ 60 000 columns and over 600 000 non-zero elements. The size of the mathematical model in itself represents a significant challenge regarding feasibility, stability and solving time. 9

Figure 5: The Initial Plan (IP) process. It is the planning and optimization process in MOL Group, done on a monthly basis, to provide optimized plans for the whole Downstream business for the upcoming 3-4 months. The first step in the planning is the data supply, based on which operational and external (economic) data is updated in the model. Based on this “raw” data the first mathematically feasible case (C00) is created. This is followed by an iterative process, together with scheduling, logistic and other stakeholders the plan is finetuned to reflect the operational reality. The Anomaly Detector aims to support this iterative process.

On top of this the data management and business understanding needed for the optimization and economic calculations require years of experience in the industry and a deep understanding in the field of optimization with linear programming. The key challenges of data management is the data supply process, which is the first step of every planning cycle and includes dozens of stakeholders providing their corresponding data from external environment; as well as operational parameters and planned inventory levels to be set as drivers or constraints in the model. Due to the sheer amount of data in this process, mistakes (like missing or wrong values) and inconsistencies (illogical values) or discrepancies (e.g. MIN>MAX)) are unavoidable, and we need to handle them. This is a classic case of “garbage in garbage out”. 10

Another key challenge of optimization stems from the complexity of interpreting the outputs (results) of the model. We need to make sure the result is not just mathematically optimal, but it is operationally viable as well, and makes sense economically. Interpreting results, understanding relations and the economic drivers require deep domain expertise, analytical thinking, and years of experience. In refinery optimization, monthly plans are being created, which means one month is one period. What this means is when the monthly balance is being calculated PIMS only sees the beginning and end point of the month, which in this case is the opening and closing inventory of a given material. It does not check the rundown of the intra-month inventories. Just because PIMS found a plan in the model feasible, it does not mean it will be feasible in real life (scheduling) as well. A detailed explanation covering the gap between planning and scheduling models in process industries can be found in [15]. This gap between planning and scheduling can be illustrated by the following example: We have two units, Unit A and Unit B. Unit A is producing the feedstock of Unit B. We have an inventory tank for this material, with minimum of 5 kt and maximum of 60 kt storage capacity. We are going to examine several scenarios in a 30 day long period.

Base scenario: Inventory change is +15 kt

Scenario A: The shutdown occurs on the 20th of the month.

Base scenario: Unit A and B are operating on maximum capacity: Unit A (feed producer) = 150 kt/month (5 kt/day), Unit B (consumer) = 135 kt/month (4.5 kt/day), Inventory change = 15 kt/month (0.5 kt/day). 11

If we would like to analyze a 10 day long shut down of Unit A, the monthly balance would look like this from optimization point of view: Unit A = 100 kt/month, Unit B = 115 kt/month, Inventory change = −35 kt/month. In the actual scheduling, we can have several different scenarios effecting the feasibility of the plan: Scenario A: Unit A has a planned shut down for 10 days, from the 20th of the month: Inventory change = −35 kt/month: 0.5 kt/day until the 20th, −4.5 kt/day from the 20th. This scenario is still feasible, but reaches the inventory minimum by the end of the period.

Scenario B: The shutdown occurs on the 10th of the month (infeasible, because we go below the inventory minimum).

Scenario C: The shutdown occurs on the start of the month (infeasible, because we go below the inventory minimum).

Scenario B: Unit A has a planned shut down for 10 days, from the 10th of the month: Inventory change = −35 kt/month 0.5 kt/day until the 10th. −4.5 kt/day between the 10th and the 20th, 0.5 kt/day from the 20th. The scenario is infeasible, because we go below the inventory minimum. Scenario C: Unit A has a planned shut down for 10 days, from the 1st Figure 6: Optimal Scenario: The shutdown of the month: Inventory change = does not cause an inventory deficit or over−35 kt/month, −4.5 kt/day until flow. 12

the 10th, 0.5 kt/day from the 10th. The scenario is infeasible, because we go below the inventory minimum. To avoid this, we will need to set a 10 kt lower throughput on Unit B, and a 10 kt higher closing inventory compared to the Base scenario. These, and similar settings that are overarching the discrepancy between planning and scheduling have to manually be input into the LP optimization in an iterative manner and require experience and a close cooperation between the corresponding teams. Interpreting results, understanding connections and the economic drivers require deep domain expertise, analytical thinking, and years of experience. The aim of the Anomaly Detection tool is to aid us in this process, by evaluating past cases and recognizing similar patterns, just like a seasoned expert. 3.2. Detecting unusual values Any value or pattern its own could never be considered an "anomaly" without knowing how this value usually behaves. This is why considering historical data is necessary for Anomaly Detection. For example, if a price of a product in the current period is 500$, but it used to move in the range 100$ − 200$, then we can be almost sure that it is an anomaly. Though this also depends on how many historical cases we collected. If the 100$ − 200$ range originates from only a couple datapoints, then the 500$ value may not be considered as anomalous as if we had many historical examples. This leads to further considerations in the actual methodology (See Section 4). 3.3. About the input data In AspenTech PIMS we are differentiating between two types of data in the A matrix and the b vector: structural and transactional. Structural data describes the complete asset setup from the refinery configuration, technology, to logistics and market points. It is also referred to as “the model” itself. These structures are relatively fix and not changing as often, only when new units are added, existing ones revamped or corrected due to differences compared to reality are found. Since the plans’ structures are constantly changing as units are revamped, this may lead to new variables introduced or removed. Therefore the size of

13

A and b, as well as x are subject to change, which leads to further difficulties when applying the anomaly detector on historical data (See the last 3 paragraphs of Section 4.2.1). Transactional data on the other hand describes the actual operational parameters of the refinery, from available unit capacities, inventory levels, logistical constraints and technological parameters to crude supply, contracted sales and purchases and external environment. These are called “input data” as they are changing from month-to-month and we require data supply to keep them updated. The anomaly detection tool is aimed at these. In MOL Group’s PIMS model we are investigating the following types of input data with the anomaly detection tool: • Sales – materials sold (from the from refinery) • Purchase – materials purchased (to the refinery) • Capacity – refinery unit capacities • Proclim – refinery process limitations • Bounds – bounded variables to represent various constraints • Transfer – material transfer between plants and depos • Blending – fuel blending numbers • Inventory – refinery and depo inventories • Material Balance – material balance of all streams Some of the constraints required for these categories are non-linear by nature, and additional sequential linearizations with piecewise linear functions are required. 3.4. Business opportunities Although a comprehensive methodology and several reports were developed in MOL Group’s optimization team to support the understanding of the drivers and opportunities of this extremely complex optimization dilemma, these face two challenges:

14

1. Time-series: the previous methods were focusing on single case data – there was no connection provided to historical results, therefore any pattern or anomaly recognition could only be performed by the optimizer expert, and as such was limited by the expert’s time, experience and analysis skills. 2. Limited number of KPIs monitored: although the most important results were checked and monitored in simple dashboards, many relations remained hidden from the tens of thousands of data points. 3.5. Technical challenges Any unsupervised machine learning algorithm, just like Anomaly Detection, expects clean data, with all input values clearly available for all historical cases. However, the data listed in Section 3.3 was available though multiple different SQL data tables, which not only needed to be collected, but also joined through multiple types of key columns. This procedure was encoded, with an additional interface for modelers to change which SQL tables data is extracted from, and on which keys they are joined. 4. Anomaly Detection methodology Multiple Machine learning algorithms have been developed for anomaly detection, grouped into supervised and unsupervised learning. See a detailed overview of process fault detection and diagnosis in chemical engineering in [16]. We are now interested in the unsupervised anomaly detection case, since we do not have labeled data. DBSCAN and the Clustering-based Local Outlier Factor are two clustering-based techniques for outlier detection. Neighbor-based techniques consider outliers as lonely observations "far" from others observations. K-nearest neighbors (KNN) is a simple technique in which the average distance to k nearest neighbors is computed based on distance metrics such as Euclidean, Manhattan or other distances. The points with the biggest average distances are marked as outliers [17]. Additionally, dimensionality reduction methods, such as PCA; as well as statistical monitoring methods are suggested in [18]. Density-based algorithms cluster the points based on the distance between a data point and a cluster. Singular points, or points in the lower density regions, are marked as outliers. Local Outlier Factor and Connectivity-Based Outlier Factor are two widely used density-based algorithms that are very 15

efficient at finding high-density regions (normal instances) and low-density regions (outliers). Finally, angle-based methods consider not only the distance between points but also the directions of the distance vectors. An angle-based outlier factor is calculated, and data with small values are labeled as anomalies. [19]. Proximity measures are easier to calculate than to fit statistical distributions, but the high dimension makes the distance-based methods unusable and also the distribution fitting hard. This is why in the next subsection we tackle the high-dimensional problem by using the most informative pairs. 4.1. Univariate Anomaly Detector - ECOD The ECOD Anomaly Detector was originally developed by Z. Li et al. in 2022 [20]. The idea was to calculate the (one-dimensional) cumulative distribution functions (c.d.f.’s) of each input variable in the historical (train) data Fi and compare each input variable of the new monthly plan’s data to these Fi ’s. If the value of the c.d.f. Fi (x) at the new datapoint x is extraordinarily high or extraordinarily low, then it must be an anomaly.

Figure 7: The location of anomalies expected, depending on the skewness of the historical data’s CDF. Skewness captures which side the long tail is located on, which is where anomalistic values could show up (dark red).

If the difference between 0 and Fi (x) (on the left side of the c.d.f. function), or the difference between 1 and Fi (x) (on the right side of the c.d.f. function) is extremely small, then the given value is an anomaly. In fact, this difference could be used as an “Anomaly Score”. This is useful because sorting by this Anomaly Score allows us to find the largest anomalies. The smaller this value, the farther away the new input value is from its old distribution, so the more likely that it is an anomaly. Moreover, in the multivariate case, 16

these scores are comparable to each other, which makes it possible to select the most anomalous input indices. More precisely, [20] suggests that we find the skewness of Fi . In the case of right skewness, the anomalies accumulate on the left of Fi while in the case of left skewness, the anomalies will be on the right. Therefore, depending on the skewness of Fi the following value could be used as an anomaly score: ( Fi (x) in the case of right skewness 1 − Fi (x) in the case of left skewness However, in case of anomalies, the listed Fi (x) and 1 − Fi (x) values get exponentially small in the tails of the c.d.f., so the use of − log2 was suggested to yield the final ECOD anomaly score of ( − log2 (Fi (x)) in the case of right skewness ECOD(x) = − log2 (1 − Fi (x)) in the case of left skewness This way, large ECOD Anomaly Scores mean large anomalies, small scores mean regular values (non-anomalies). This is as far as the original ECOD methodology goes, however, when actually applying this method on real data, we have found the following shortcomings: 1. Very often, the value of log(0) = −∞ would be encountered. That means that (in the case of right skewness), Fi (x) = 0, meaning that in the historical data, every value of the given input was greater than x or in other words, we have never seen such a small x value in the past. Or in the case of left skewness, we have never seen such a large x value in the past. So for example, when the entire system is flawed, e.g. it encounters a loophole of being able to transfer/buy/sell infinitely, such large values can occur very frequently. All of these cases resulted in many ECOD Anomaly Scores of −∞, which could no longer be compared to each other. 2. Quite often, the Fi functions were calculated from only a few data points, which meant that the ECOD Anomaly Scores weren’t telling the whole truth. This has also yielded cases of log(0) = −∞ with the only difference that sometimes entirely correct values would also gain an infinitely large anomaly score. 17

To solve these issues, we have applied the following modifications: 1. We have redefined the ECOD Anomaly Score to be ( c · miny∈Train Data (|x − y|) if ECOD(x) is undefined ECOD′ (x) := ECOD(x) otherwise where c is a scaling factor that indicates how much more important it is for the anomaly score to indicate that the given value does not exist in the train data. By default, it is set to c = 10, but we have also made it possible for the user to change it in the future. Moreover, if the first case is chosen, the output also indicates that this x value was never seen in the input dataset, which we have called “AA type anomalies”. A cutoff was chosen for the rest of the cases. The first 5000 anomalies found this way (sorted by the ECOD Anomaly Score) were called “A type anomalies”. 2. We have introduced three more cutoffs to the system. Firstly, constant c.d.f.s were discarded. Secondly, the train data would only contain c.d.f.s that are calculated from at least k historical samples (where we have chosen k = 5 for this parameter). Thirdly, the proportion between the datapoints that the given c.d.f. is calculated from and the number of monthly plans the historical data contains must be at least p (where we have chosen p = 0.05 for this parameter). These two empirical cutoffs avoid calculating the c.d.f. from too few data points, and we will not be checking for anomalies in these cases. 4.2. Multivariate Anomaly Detection Anomaly Detection in high dimensions is a multifaceted problem, and the choice of detection techniques is strongly influenced by how anomalies are defined, as well as by the nature of the input data and the expected outputs. These factors lead to a wide range of problem formulations, each requiring different analytical and computational approaches. The presence of anomalies, which may be data flaws or measurement errors, can lead to model misspecification and misleading results in the case of classical models such as Linear Programming. In contrast, there are cases when anomalies are the main carriers of important information, which may lead to significant economic opportunities. 18

In the following subsections, our anomaly detection methods are introduced that help the experts improve the quality of data and enhance the work of the decision-makers. 4.2.1. A pair-selection method for bivariate Anomaly Detection In unsupervised Anomaly Detection in high dimensions, a main problem is caused by the sparsity of the data. To overcome this, we apply the idea introduced in [21] to use a subset of the bivariate marginal distributions. At the core of this method is the goal of finding an approximating multivariate probability distribution, defined by only bivariate probability distributions, chosen from the marginals of the approximated one, with the property that the Kullback-Leibler divergence is minimized. In [22] it is shown that finding the maximum weight spanning tree in a complete graph, in which the vertices are associated with the random variables and the weights are defined by the mutual information of connected variables, determines the bivariate marginals included in the approximation. The Anomaly Detection method presented in [21] employs a parametric approach for modeling the bivariate marginals by separately fitting the bivariate copulas and the marginals. Since in our case we deal with an extreme large amount of data and variables, we decide to use Gaussian bivariate marginals, which are faster to fit and are more robust. In this paper, we introduce two techniques to identify anomalies in the bivariate marginals, as described in the following subsections. In the Section 4.1, we have seen how one-dimensional c.d.f.s Fi can be used to find rare, unusual values, which we have called anomalies. In the following two subsections, we will be using certain two-dimensional marginal random vectors (Xi , Xj ), where i and j are two different indices of the column vectors of the LP model. Now due to the high-dimensionality of the data we will explain how we select the important pairs (Xi , Xj ); then explain how to apply two different types of anomaly detectors. The used methodology goes the following way: 1. Pre-define groups of input pairs we would like to calculate anomalies for. In the case of large oil refinery models, unfortunately, calculating an Anomaly Score for all possible pairs seems to not be feasible, not to mention the interpretation of a high Anomaly Score between certain variables becomes fuzzy. For example, there is no point in comparing the Marginal Value of a purchase to the capacity of a completely different plant – even if such a pair has unusual values, the planning team 19

wouldn’t be able to use this information. Therefore, a manually defined pre-selection of variable groups is necessary, which in our case, reduced the number of pairs to look through from ∼ 7 billion to ∼ 4 million. This was done by the expert modelers, then afterwards the Anomaly Detector automatically does this filtering every time it is run. 2. Further filter the potential pairs using Kendall’s Tau. Calculate a Kendall correlation matrix between all pairs defined in the first step. To do this, we create a complete graph with nodes for each variable, and edge weights equal to the absolute value of Kendall’s Tau between each connected pair. Then we calculate a maximum weight spanning tree, and only keep edges (pairs) that are inside this spanning tree. This avoids listing the same anomaly source twice, and only lists the most significant anomalies. For example, if the values of all of i, j and k are extreme, then all pairs (i, j), (j, k), and (j, i) would be listed as anomalies. Using a spanning tree however, one of these pairs – the one with the lowest two-dimensional anomaly score – will be discarded. Finally, we obtained pairs with mostly high Kendall correlation. After this procedure, some low-Kendall correlation pairs have still remained, so we decided to implement a cutoff, K, and remove pairs that have an absolute Kendall correlation below K (in our case, K = 0.4). This makes sure we are really only dealing with existing relationships. In our case, this process reduced the number of potential pairs to look through from ∼ 4 million to ∼ 20 thousand. 3. Then, we are preparing for the two multidimensional Anomaly Detection methods described in Sections 4.2.2 and 4.2.3. For all of the (x, y) pairs obtained, calculate the best fitting line y = ax + b, an R2 linear fit score of this line, an average e and a variances s2e over a fit line (see Formula 3). Then, for all the pairs obtained, also calculate µ, the 2dimensional historical mean of the pair’s values, and a 2 × 2 covariance matrix V from the same historical data. In our experience, calculating such a large correlation matrix in Step 2 takes several hours. To avoid having to wait this much every time the Anomaly Detector runs, we have applied the train – retrain methodology. This means that to “fit” this Anomaly Detector, we first need to train (or retrain) it, where Steps 1-3 will be calculated. When applied to new data (tested), only the steps described in Sections 4.2.2 and 4.2.3 will be applied. There is however one caveat in Step 2 of this methodology. In reality, there 20

are a huge number of missing values. For one, not all input values are present in all monthly plans, in fact, every plan is usually missing a fair number of input values. However, a correlation matrix calculation expects input vectors of the same length, and we are dealing with non-equal length input vectors (the amount of missing values is different between each input). Matching up the correct values of each input variable’s historical range for every element of the correlation matrix turned out to be an incredibly demanding task from a computational point of view. To resolve this, we are using a penalty function method. This method calculates the correlation matrix normally: It fills in actual missing values with the average of the given input’s values. Then it applies a multiplier to each row and each column of the correlation matrix. This multiplier ranges between 0 and 1, and is closer to 1 the more missing values we have. It makes sense to linearly penalize the correlation matrix, namely if for example 30% of the values were missing, then we multiply the given correlation by 1 − 0.3 = 0.7. However, each row and column may need to be multiplied, so we instead multiply by the root of this, so that two roots will cancel out to obtain a linear penalization. q , where More specifically, we multiply each row and each column by n−s n s is the number of values considered to be “missing”, and n is the number of historical cases. This is the actual correlation matrix we use in Step 2. This reduced the calculation time of the correlation matrix from more than 10 hours to around 2 − 3 hours. 4.2.2. Linear Regression Anomaly Detector Anomalies do not just come alone. Sometimes two different LP model values can individually be in the expected range, but together still be far from the usual. See Figure 8. The two coordinates of the red dot represent the current values of the x and the y variables, respectively. When projected to the x axis, the red dot will be lost in the range of all the other black dots, and same for the y direction. Yet the rest of the black Figure 8: Example of an 2dots all lie near the blue line, while the red dimensional anomaly that would dot lies far away. In this case, the y value not be detected in one dimension 21

was historically approximately some linear function of x, let’s say y ≈ ax + b, except in the red case.

Figure 9: The Anomaly Score must depend on the variance of the train data. While the red dot is the same distance away from the blue line in both cases, it is a relatively realistic value on the left, and an anomaly on the right.

The goal of the Linear Regression Anomaly Detector is first to find such variable pairs for which historical data would suggest a strong linear correlation (a.k.a., the y variable is usually very close to ax + b for some values of a and b), and for these pairs, give the best estimation for a and b, and a linear fit score. The linear fit score defines how strong this correlation is – see Figure 9: In the left case, x and y have a weak correlation, so low linear fit score, and in the right case, this correlation is high. Even though the red dot is the same distance away from the blue line in both cases, since the linear fit score on the left is smaller, the red dot is less of an anomaly than on the right. The linear fit score mentioned is actually the Pearson-correlation (R2 ) score. When applying the Linear Regression Anomaly Detector, we can read the new input data, run through all 20 thousand pairs filtered in Section 4.2.1, and see if any of them are far outside of the expected y = ax + b line. To quantify this, we calculate a pairwise Anomaly Score, in the following way: First we draw a normal distribution over the y = ax + b line with variance equal to the train data’s variance around this line. Then, the Anomaly Score will be the value taken on by this normal distribution at the new pair’s coordinates. This means that a lower Anomaly Score actually describes a greater anomaly. In terms of formulas, the Anomaly Score of a new pair of 22

values (xn+1 , yn+1 ) will be following:   1 (yn+1 − (axn+1 + b) − e)2 AnomalyScore(xn+1 , yn+1 ) = p exp − 2s2e 2πs2e where n

1X e= (yi − (axi + b)) n i=1 s2e =

n 1 X (yi − (axi + b) − e)2 n − 1 i=1

(3)

y = ax + b is the best fitting line on the (x, y) data The parameters a, b, e and s2e were calculated and saved in Step 3 of Section 4.2.1. 4.2.3. Multivariate Sampling Anomaly Detector In the Multivariate Sampling methodology, we defined scores based on ellipses, obtained from the level sets of a 2-dimensional multivariate distribution, calculated from historical data. On Figure 10, the red ellipse is the best fitting ellipse to the data, such that it has the lowest area that still surrounds 99% of the historical data, and only exactly 1% of the historical data is outside. A 2-dimensional normal distribution is defined by a mean vector µ (2 values, the in- Figure 10: The 99% level set of a dividual averages of the x and y directions), 2D multivariate normal distribution fitted to a variable pair of historical and a 2×2 covariance matrix V . These were values. already calculated in Step 3 of Section 4.2.1. Afterwards, the Multivariate Sampling algorithm will be the following: We calculate a pairwise Anomaly Score. A fitted normal distribution already exists, so the Anomaly Score will simply be the value this function takes on at the new data’s coordinates. A p-cutoff-value can also be calculated. This value expresses how rare such a pair is. If the p-cutoff-value is 0.01 or 1%, then the pair of values is outside the above shown ellipse. If the 23

p-cutoff-value is 0.05 or 5%, then the pair is outside the ellipse defined by the 5% cutoff. These two cutoffs are calculated by Multivariate Sampling (the namesake of the method). Since the data is sparse, we sample 1500 normally distributed two-dimensional values using a Monte-Carlo simulation, with the same mean µ and covariance V as the historical data. (This 1500 can be changed according to the size of the original problem.) In terms of formulas, the Anomaly Score of a new pair of values (xn+1 , yn+1 ) will be following: AnomalyScore(xn+1 , yn+1 ) =   1 1 T −1 exp − ((xn+1 , yn+1 ) − µ) V ((xn+1 , yn+1 ) − µ) =p 2 (2π)2 det V where µ and V are calculated and saved in Step 3 of Section 4.2.1. 4.3. Assembling the full Anomaly Detector Figure 11 contains an example historical range of a pair of inputs considered. We then apply two Anomaly Detection methods introduced in Section 4.2.2 and Section 4.2.3. Both of these methods have a 1% cutoff based on historical data. (Outside of this cutoff, historically, only 1% of the data was found.) The intersections of these regions are interpreted as different types of anomalies. • Non-anomalous region (green): The region that neither method considers an anomaly. • Significant anomalies (orange): The region that both methods consider an anomaly. Values found here are very likely to be anomalous, therefore the "significant" marker. • Disproportionate anomalies (purple): The region that Linear Regression considers an anomaly, but Multivariate sampling does not. Since value pairs "close" to the Linear Regression line have approximately the same proportion, values "far away" from this line can be considered disproportionate. • Suez-type anomalies (purple): The region that Linear Regression does not consider an anomaly, but Multivariate Sampling does. In these regions, value pairs are both either unexpectedly large or unexpectedly small (e.g. two prices skyrocketed), but compared to each 24

other, they are still about the same. The name comes from the fact that in 2021, famously, a ship got stuck in the Suez-canal, which lead to such anomalies - prices went up, but all by approximately the same amount. [23]

Figure 11: The regions determined by the Anomaly Detector for each pair of historical values considered. • The red ellipse is determined by Multivariate Sampling, and it is the best fitting ellipse to the datapoints such that 1% of the historical data is outside. • The green line is determined by Linear Regression and it is the best fitting line to the data. The dashed green line is chosen such that 1% of the historical data is outside of it. • The central green region is the "non-anomalous" region. • The two purple regions are where Disproportionate Anomalies are found. • The two blue regions are where Suez-type Anomalies are found. • The two large orange regions are where Significant Anomalies are found.

As mentioned in Section 4.2.1, the program needs to be periodically retrained. We chose a period of 3 months for this. Figure 12 shows what the (re)training algorithm does, and how it affects the main Anomaly Detection algorithm. 25

26 Figure 12: The outline of the Anomaly Detection process. First the Anomaly Detector is fit using historical PIMS Solutions. This calculates the 1D CDFs, important & highly-correlating pairs to look through in the multidimensional case, and the parameters a, b, e, s2e , µ, V for each considered pair (x, y). Then the Anomaly Detector uses these pre-calculated values to find the one- and two dimensional Anomaly Scores for each considered value or pair. The resulting 2D anomalies are transformed via the Venn diagram shown for an easier interpretation of the Anomaly types.

5. Illustrative Results In this section, we will showcase some of the Anomaly Detection findings. Due to the data being sensitive, we will multiply all numerical values by random number between [0.9, 1.1]. We will present three examples of anomalies we have found. In general, the LP solver outputs a result, which, among other data, contains the current plan’s Marginal Values. Sorting them in a decreasing order of their absolute value, allows us to filter out the particularly outstanding values of the current plan. However, this does not consider the historical values of them. The Anomaly Detection methodology does use this historical information, and sorts Marginal Values by an Anomaly Score, which indicates how different they are from the historical distribution. If the current plan contains an order of magnitude jump Marginal Value, that will be noticed. Marginal Values are usually of varying order of magnitude. If such an order of magnitude jump occurs for typically small Marginal Values, those go unnoticed. Large Marginal Values of a constraint may indicate "mathematical strain". This means that if we just had 1 more unit of material to satisfy the constraint, then our profit would increase in accordance. This is also called a "tight" system. Our goal, in all cases, is first of all to correct all unrealistic values and typos; and secondly, to gain a better business understanding and a trust in the LP solution. The planners need to be confident that their plan is without issues when handing it over to the scheduler colleagues. Example 1 This is an example of such a large Marginal Value. A Blend constraint in one of the plans had a Marginal Value of 80143.31. This was significant, because of the 5004 cases that we trained on, 3501 contained an entry for this Marginal Value. However, historically, this Marginal Value was always in the range [−1.02, 39.41], and most of the time it was below 2. This was an internal material flow, whose Marginal Value is expected to be 0, or at least, very close to 0, which historically held. A price of a barrel of oil is in the range of $50 − $80. The barrel factor is around 7.5 − 8. This ratio compares the price of products sold that were produced from the given barrel to the purchase cost of the barrel of material. This means that the materials bought for the above $50 − $80, after chemical 27

transformations, are sold for around $350 − $700. This does not mean that the profit is $300 − $600, because the cost of refining takes up the majority of this price. In fact, a healthy margin on a given barrel is around $3 − $4. Marginal Values represent a similar idea: they indicate how much profit we could generate by adding infinitesimally more input material. The value of this infinitesimally small material should realistically not be greater than the margin on a barrel of material. However, in this case, the obtained Marginal Value was much larger than expected, 80143.31. The potential causes of this may be one or more of the following: 1. The output of the LP was almost infeasible. 2. The constraint corresponding to this Marginal Value is "tight". 3. Tight specifications on the recipe in the current situation. 4. The system cannot use up the materials included in this blend. 5. etc. The potential solutions for 1 − 4 are the following: 1-2. Along with other large Marginal Values, find related capacities / maximum values that are possible to increase to "loosen" the system. 3. If the specification is tight, release it. 4. Find material balance limits that can be released. Or provide the input materials from other inventories as well. An example for altering the amount of material used for this blend from a given inventory is to increase it from by about 1 t, as seen on Figure 13. This requires modifying other constraints as well, since materials from the same inventory may be used elsewhere. Compared to the total material balance, this change is minor, but mathematically it releases the constraint, obtaining a realistic Marginal Value. Through this process, we can iterate the model to find realistic values, in order to hand over a trustworthy LP solution. Example 2 This is an example of a pair-anomaly we found. We are doing a weightproportional blending of several materials. Two components of this blend, component proportions A and B used to historically always have a relation of B = 1.5133 · A, however, in the currently tested case, B = 0.8787 · A. 28

Figure 13: If the given inventory maximum capacity is increased by 1 ("loosening" the model), the corresponding Marginal Value may decrease to realistic values. The magnitude of the jump depends on the magnitude of the material balance value. Of course, we do not know the shape of the above curve, we can only find its values by reevaluating the model.

This is why the Linear Regression Amomaly Detector output this pair as an anomaly. Moreover, the Multivariate Sampling Anomaly Detector also found this to be an anomaly, because historically both A and B were small, below 0.08, but in the currently tested case, both, A and B were above 0.3. Therefore, this pair is in the orange zone of Figure 11, so it is a Significant Anomaly. There are two types of blending: formula blend, and specification blend. Formula blends always follow the same proportions in every blend. In the case of specification blends, the system can set these proportions (between given limits), with the goal of outputting a blend with a given quality. Whenever experts find such an anomaly, they need to check which type of blend this was. This was a formula blend, so planners need to be aware of the change. In this case, if this showed up as an anomaly, we need to disregard it, since this originated from a model change. In case of a specification blend, it is possible that the external environment changes, e.g. the price structure is different; or that the internal environment changes, e.g. which plants are under maintenance, what kinds of materials 29

are available. In this case, achieving the same quality may require unusual blend proportions. If we believe that these were not causing the issue, then data supply for the limits or capacities can be incorrect. If the model ran without infeasibility, and with correct data supply, then it was capable of fulfilling the quality constraints. So in this case, even if this blend is unusual compared to historical blends, then it should not cause any issues. Example 3 This is an example of two inventory capacity maximums, A and B. Historically we have found that B = 0.5916 · A − 0.6159 holds with an R2 = 0.8225. However in this case, A and B were extremely far away from this line. The potential causes for this anomaly may be data supply issues, or an actual inventory capacity increase or decrease. Sometimes, due to operational requirements, schedulers can add or remove extra tanks, which could be the cause of this anomaly. If this was not an operative change, then experts will need to discuss with the colleagues responsible for this inventory whether or not this was an actual inventory capacity increase (physically more barrels are available or additional warehouses are rented), or if it was an error. 6. Summary This paper emphasizes the need to reconsider the concept of anomalies. Besides their well-known applications, they can serve as a tool for expert reasoning, enhancing more effective business strategies. Data gathered from various sources, such as external datasets, internal records, or software may overlook key insights relevant to experts. We have shown how the traditional LP methodology can be enhanced with AI methods such as Anomaly Detection. This new approach is not only capable of finding unusual, unrealistic and unattainable values that hinder planning, but also uncover hidden potentials and bottlenecks in each plan. By using not only the input and output values of the LP problem, but also each constraint’s Marginal Value, the found anomalies may indicate new business opportunities that could not have been considered without Anomaly Detection due to the sheer size of the data. We have split this issue into two sub-problems: one- and multidimensional Anomaly Detection. For the one-dimensional case, we have applied 30

and modified the well-known ECOD methodology to fit our problem. For the multidimensional case, we chose the most informative pairs, first using grouping based on expert knowledge, and secondly, automatic Kendall-correlation and spanning tree-based filtering. We have introduced two new methodologies for finding 2D anomalies in this filtered dataset. Using historical data, we have trained a Linear Regression-based and a Multivariate Sampling-based algorithm, which can then be used to detect 3 types of 2D anomalies. By choosing fitting names for these groups, they can then be easily interpreted by experts. Using these introduced methods, we have found several 1D and 2D anomalies. Due to the sheer amount of data available, it would be impossible to manually find meaningful anomalies. Once the Anomaly Detector found the largest anomalies, and sorted by Anomaly Score, then experts can take a look at them, and using their experience, they can interpret these results with improved business understanding. We believe that this work opens a new avenue for gaining insight into business opportunities hidden within Big Data, by applying the Anomaly Detection methodologies introduced. In the future, one may extend this methodology into finding 3D anomalies, and making them interpretable. Alternatively, one may consider using distributions other than the Gaussian distribution.

Acknowledgements The research reported in this paper is part of project no. BME-NVA02, implemented with the support provided by the Ministry of Innovation and Technology of Hungary from the National Research, Development and Innovation Fund, financed under the TKP2021 funding scheme. Additional funding was provided by NKFIH for the Cooperative Doctoral Program.

31

References [1] L. Rodríguez-Mazahua, C.-A. Rodríguez-Enríquez, J. L. SánchezCervantes, J. Cervantes, J. L. García-Alcaraz, G. Alor-Hernández, A general perspective of big data: applications, tools, challenges and trends, The Journal of Supercomputing 72 (8) (2016) 3073–3113. [2] M. Hamzehi, S. Hosseini, Business intelligence using machine learning algorithms, Multimedia tools and applications 81 (23) (2022) 33233– 33251. [3] M. Rath, Realization of business intelligence using machine learning, Internet of Things in Business Transformation: Developing an Engineering and Business Strategy for Industry 5.0 (2021) 169–184. [4] F. Ridzuan, W. M. N. W. Zainon, Diagnostic analysis for outlier detection in big data analytics, Procedia Computer Science 197 (2022) 685–692. [5] P. Larrañaga, D. Atienza, J. Diaz-Rozo, A. Ogbechie, C. E. PuertoSantana, C. Bielza, Industrial applications of machine learning, CRC press, 2018. [6] D. Carou, A. Sartal, J. P. Davim, Machine learning and artificial intelligence with industrial applications, Springer. doi 10 (2022) 978–3. [7] D. Klipa, I. Ristić, A. Radonjić, I. Scepanović, et al., Big data and artificial intelligence, International Journal of Management Trends: Key Concepts and Research 1 (1) (2022) 3–14. [8] N. K. Shah, Z. Li, M. G. Ierapetritou, Petroleum refining operations: key issues, advances, and opportunities, Industrial & Engineering Chemistry Research 50 (3) (2011) 1161–1170. [9] I. Grossmann, Enterprise-wide optimization: A new frontier in process systems engineering, AIChE Journal 51 (7) (2005) 1846–1857. [10] V. Venkatasubramanian, The promise of artificial intelligence in chemical engineering: Is it here, finally?, AIChE Journal 65 (1) (2019).

32

[11] R. Thakur, Z. Sajid, F. Khan, Artificial intelligence (ai) safety system for safe & trustworthy autonomy, Digital Chemical Engineering (2026) 100308. [12] Aspen unified pims, https://www.aspentech.com/en/products/msc/ aspen-unified-pims. [13] J. M. Pinto, L. F. L. Moro, A planning model for petroleum refineries, Brazilian Journal of Chemical Engineering 17 (2000) 575–586. [14] D. Bertsimas, J. N. Tsitsiklis, Introduction to linear optimization, Vol. 6, Athena scientific Belmont, MA, 1997. [15] I. Harjunkoski, C. T. Maravelias, P. Bongers, P. M. Castro, S. Engell, I. E. Grossmann, J. Hooker, C. Méndez, G. Sand, J. Wassick, Scope for industrial applications of production scheduling models and solution methods, Computers & Chemical Engineering 62 (2014) 161–193. [16] V. Venkatasubramanian, R. Rengaswamy, K. Yin, S. N. Kavuri, A review of process fault detection and diagnosis: Part i: Quantitative model-based methods, Computers & chemical engineering 27 (3) (2003) 293–311. [17] T. T. Dang, H. Y. Ngan, W. Liu, Distance-based k-nearest neighbors outlier detection method in large-scale traffic data, in: 2015 IEEE International Conference on Digital Signal Processing (DSP), IEEE, 2015, pp. 507–510. [18] S. J. Qin, Survey on data-driven industrial process monitoring and diagnosis, Annual reviews in control 36 (2) (2012) 220–234. [19] H.-P. Kriegel, M. Schubert, A. Zimek, Angle-based outlier detection in high-dimensional data, in: Proceedings of the 14th ACM SIGKDD international conference on Knowledge discovery and data mining, 2008, pp. 444–452. [20] Z. Li, Y. Zhao, X. Hu, N. Botta, C. Ionescu, G. H. Chen, Ecod: Unsupervised outlier detection using empirical cumulative distribution functions, IEEE Transactions on Knowledge and Data Engineering 35 (12) (2022) 12181–12193. 33

[21] G. Horváth, E. Kovács, R. Molontay, S. Nováczki, Copula-based anomaly scoring and localization for large-scale, high-dimensional continuous data, ACM Transactions on Intelligent Systems and Technology (TIST) 11 (3) (2020) 1–26. [22] C. Chow, C. Liu, Approximating discrete probability distributions with dependence trees, IEEE transactions on Information Theory 14 (3) (1968) 462–467. [23] 2021 suez canal obstruction, world-middle-east-56505413.

34

https://www.bbc.com/news/