# Event detection from novel data sources: Leveraging satellite imagery
alongside GPS traces

Ekin Ugurel  University of Washington, Seattle, WA
<EMAIL_ADDRESS>,  Steffen Coenen  DKS Associates, Seattle, WA
<EMAIL_ADDRESS>,  Minda Zhou Chen  University of Washington, Seattle, WA
<EMAIL_ADDRESS>and  Cynthia Chen  University of Washington, Seattle, WA
<EMAIL_ADDRESS>
## Abstract

Rapid identification and response to breaking events, particularly those that
pose a threat to human life such as natural disasters or conflicts, is of
paramount importance. The prevalence of mobile devices and the ubiquity of
network connectivity has generated a massive amount of temporally- and
spatially-stamped data. Numerous studies have used mobile data to derive
individual human mobility patterns for various applications. Similarly, the
increasing number of orbital satellites has made it easier to gather high-
resolution images capturing a snapshot of a geographical area in sub-daily
temporal frequency. We propose a novel data fusion methodology integrating
satellite imagery with privacy-enhanced mobile data to augment the event
inference task, whether in real-time or historical. In the absence of boots on
the ground, mobile data is able to give an approximation of human mobility,
proximity to one another, and the built environment. On the other hand,
satellite imagery can provide visual information on physical changes to the
built and natural environment. The expected use cases for our methodology
include small-scale disaster detection (i.e., tornadoes, wildfires, and
floods) in rural regions, search and rescue operation augmentation for lost
hikers in remote wilderness areas, and identification of active conflict areas
and population displacement in war-torn states. Our implementation is open-
source on GitHub: https://github.com/ekinugurel/SatMobFusion.

## Introduction

In an increasingly interconnected world, the ability to rapidly identify and
respond to breaking events, particularly those that pose a threat to human
life such as natural disasters or conflicts, is of paramount importance. The
effectiveness of response efforts can be significantly impacted by the speed
and accuracy of information gathering, analysis, and dissemination, with the
potential to save lives and mitigate damage. However, obtaining such
information poses significant challenges due to various factors. For example,
there have been instances where news agencies faced difficulties in promptly
or completely delivering information about ongoing events due to a scarcity of
reporters present in the affected areas Stelter and Cohen (2008). The ubiquity
of mobile devices and the near-constant network connectivity that
characterizes our modern era have resulted in the generation of a massive
amount of temporally- and spatially-stamped data (hereafter called mobile
data). These datasets have been the subject of numerous studies aiming to
derive individual human mobility patterns for various applications, including
quantifying urban vitality Sulis et al. (2018), modeling epidemic spread
Alessandretti (2022), and inferring the activities of individuals Liao et al.
(2007). Simultaneously, the increasing number of orbital satellites has
facilitated the collection of high-resolution images, capturing snapshots of
geographical areas with sub-daily temporal frequency. In this paper, we
propose a novel data fusion methodology that integrates satellite imagery with
mobile data to identify breaking events in real time. This approach leverages
the strengths of both data types: mobile data provides an approximation of
human mobility, proximity to one another, and the built environment, while
satellite imagery offers visual information on physical changes to the built
and natural environment. This method holds significant potential for three
critical applications. First, it can enhance the detection of small-scale
disasters such as tornadoes, wildfires, and floods in rural regions, where
sparse populations and limited infrastructure can impede timely detection and
response. Second, it can augment search and rescue operations for lost hikers
in remote wilderness areas, providing critical geospatial data to guide rescue
teams and improve the chances of a successful rescue. Before and after images
of the Aru Glacier avalanche in Figure 1 suggest that large-scale changes in
the natural environment can be inferred rather straightforwardly using
satellite imagery Tian et al. (2017). Finally, it can aid in the
identification of active conflict areas and population displacement in war-
torn states, offering invaluable insights to humanitarian organizations and
policy-makers. By leveraging the strengths of both satellite imagery and GPS
traces, this approach overcomes the limitations of both data types and
augments the ability to monitor, understand, and respond to dynamic events in
diverse contexts.

Figure 1. Images before and after the Aru Glacier collapse on 17 July 2016.
Images are from Sentinel-2 on 18 June 2016 (a) and 21 July 2016 (b) with 10 m
resolution. From (Tian et al., 2017)

We demonstrate an application of the proposed method on a small-scale tornado
that ripped through Muskogee, Oklahoma on May 15th, 2020. By combining mobile
data from the area as well as satellite images before and after the event, we
infer the occurrence of a natural disaster, identify damaged areas, and
highlight hot spots for individuals’ travel through an area (both spatially
and temporally).

### Our Contributions

  * •

We develop a practice-oriented data fusion methodology to augment the
inference task for a variety of use cases;

  * •

We outline expected use cases and demonstrate the case study of a small-scale
natural disaster for which the proposed method is appropriate, and;

  * •

We provide an open-source implementation of the discussed methods.

The rest of the paper is organized as follows: Related Work discusses previous
work on leveraging satellite images, GPS traces, social media posts, and the
fusion of such datasets in a variety of contexts. Dataset and Software
describes the datasets and the open-source software we employ. The Methodology
section details our methods, which include data preprocessing steps as well as
analysis specific to each data type. The Experiments section outlines our case
study on a tornado in Muskogee, Oklahoma, while the Use Cases section
describes other potential use cases for the proposed method. We conclude with
a brief discussion of the implications of our findings.

## Related Work

Historically, capturing changes in the built and natural environment often
involved on-the-ground data collection and analysis Haddad (2011). Surveys and
field observations have been the most traditional method to capture such
changes, and continue to be so today in many industries. Prior to the advent
of modern satellites, aerial photography taken from aircrafts was another
common method to document large-scale changes to the built environment Bewley
(2003). While these two methods of data collection are sufficient for many
studies, they (a) lack the immediacy called for in urgent real-time scenarios
and (b) tend to miss developments in the absence of the proper equipment or
manpower. In the last three decades, new methods of data collection in the
built and natural environment have become prominent, including the use of
satellite images, GPS traces from mobile devices, and user-generated content
from social media platforms Jongman et al. (2015); Ji et al. (2018); Yin et
al. (2021).

This literature review aims to provide a comprehensive understanding of
current and recent research in the fields of mobile data analysis, satellite
imagery, and crowdsourcing, highlighting works that have fused any two of them
to assist in the inference task. Table 1 summarizes some of the papers
discussed in the following paragraphs.

Table 1. Studies that have used mobile data/satellite imagery for various purposes. Author | Date |  Datasets Used |  Applications
---|---|---|---
Jongman et al. | 2015 |  Satellite imagery,
Twitter posts |  Early flood detection in the Philippines and Pakistan
Molinier et al. | 2016 |  Satellite imagery,
mobile app data |  Deriving biomass maps of forests
Fayne et al. | 2017 |  Satellite imagery |  Flood detection in the Lower Mekong Basin
Ji et al. | 2018 |  Satellite imagery |  Collapsed building detection in 2010 Haiti Earthquake
Sulis et al. | 2018 |  Transit smart card data, Twitter posts |  Quantifying diversity and vitality in London
Said et al. | 2019 |  Satellite imagery,
social media posts |  Natural disaster detection
Pokhryiyal et al. | 2020 |  Satellite imagery,
GPS traces |  Estimation of poverty levels in Haiti
Yin et al. | 2021 |  Satellite imagery,
GPS traces |  Road attribute detection in Singapore and Jakarta
Belcastro et al. | 2021 |  Social media posts |  Small-scale emergency detection

### Satellite Imagery and Remote Sensing

The use of satellite imagery for event detection and inference has been well-
documented. Satellite images provide a bird’s-eye view of the Earth’s surface,
enabling the detection of physical changes over time. Though the granularity
of satellite images alone may not provide the nuanced information needed in
certain scenarios, they are useful in understanding big-picture changes on the
built and natural environment. Particularly in the last five years, the
proliferation of commercial satellite imagery providers has enabled access to
this data type at a short notice and in real-time. This has been used in
various applications, such as disaster management Voigt et al. (2007), urban
planning Albert et al. (2017); Karunanithi et al. (2016), and environmental
monitoring Xu (2010).

Recent advancements in machine learning have further enhanced the capabilities
of satellite imagery analysis. For instance, convolutional neural networks
have been used to automatically detect and classify objects in satellite
images Zhang et al. (2016). Meta’s novel object classification method, Segment
Anything Kirillov et al. (2023), has been proposed for the geospatial domain
with great promise Wu and Osco (2023). These developments have opened up new
possibilities for event detection and inference, as it allows for the
automated analysis of large volumes of satellite imagery and greatly reduces
image processing time, allowing for subsequent immediate response on the
ground.

### GPS Traces from Mobile Devices

Parallel to the advancements in satellite imagery analysis, the proliferation
of GPS-enabled mobile devices has led to an explosion of mobile data,
providing detailed information about the movements of individuals and groups
Zheng et al. (2008). Research has delved into the use of mobile data for
tracking population movements Bengtsson et al. (2011), understanding social
behavior Hong et al. (2008); Calabrese et al. (2010), predicting and managing
traffic patterns Herrera et al. (2010), and more.

The richness of this data source has also led to the development of new
methods for event detection and inference. For instance, GPS traces have been
used to detect social events Castro et al. (2014); Pan et al. (2013), infer
the activities of individuals Liao et al. (2007), and understand where
individuals conduct their daily activities Chen et al. (2014). However, the
fusion of mobile data with other data sources for real-time event
identification and inference in the built and natural environment,
particularly in the context of emergencies and disasters, is an area that
deserves further exploration.

### Crowdsourcing and Social Media Data

Crowdsourcing has emerged as a powerful tool for gathering and analyzing data,
particularly in the context of event detection and response. The term
”crowdsourcing” refers to the practice of obtaining information or input into
a task or project by enlisting the services of a large number of people,
either paid or unpaid, typically via the internet Estellés-Arolas and L.
Guevara (2012). In the context of event detection, crowdsourcing can be used
to gather real-time information from individuals who are directly experiencing
or observing the event Goodchild (2007). One of the most notable examples of
crowdsourcing in event detection is the Ushahidi platform, which was used to
gather and map reports of violence in Kenya after the post-election violence
in 2008 Okolloh (2009). Since then, Ushahidi has been used in various other
contexts, including disaster response and environmental monitoring.

Crowdsourcing is closely linked to social media’s rise to prominence as a
crucial news outlet during disasters. A mountain of literature has been
published on the promises of social media data in augmenting the inference
task in a range of scenarios. Social media websites provide a platform for
individuals to share real-time information on developing events, including
their impacts and the needs of affected communities Palen et al. (2007). The
information circulated on social media tends to dissipate quicker and be more
specific than official news outlets, as anyone and everyone has the ability to
share and post Palen (2008). This information has been used by disaster
response organizations to better understand the situation on the ground and
coordinate their response efforts.

Several studies have explored the use of social media data in disaster
response. For instance, (Sakaki et al., 2010) developed a system that uses
Twitter data to detect earthquakes in real-time. Similarly, (Ashktorab et al.,
2014) developed a system that uses Twitter data to extract actionable
information during disasters, such as requests for help or reports of
infrastructure damage. (Guan and Chen, 2014) developed a metric to quantify
the evolution of disasters based on Twitter activities, demonstrating its use
for damage detection during Hurricane Sandy, affecting much of the eastern
United States and various Caribbean nations. More recently, social media data
has been combined with other data sources for disaster response. For example,
(Jongman et al., 2015) proposed a method for fusing social media data and
satellite imagery to detect floods. This approach leverages the strengths of
both data sources, using social media data to provide real-time information
about the flood and satellite imagery to provide a macro view of the flood’s
impacts. We refer the reader to (Said et al., 2019) for a broader review of
these approaches.

Despite the significant potential of crowdsourcing and social media data for
event detection and response, there are notable limitations that must be
acknowledged. The quality and reliability of crowdsourced data can vary
significantly. Misinformation or inaccurate reports can be disseminated,
either unintentionally due to confusion or panic or intentionally as a form of
manipulation or cyberattack Allcott et al. (2019). This can lead to false
positives in event detection or misdirected response efforts. Additionally,
crowdsourcing is inherently dependent on the availability and willingness of
individuals to contribute information. In certain scenarios, such as in remote
wilderness areas or war-torn regions, the number of individuals able to
contribute data may be limited. Similarly, during small-scale or rapidly
unfolding events, there may not be sufficient time for crowdsourced
information to be gathered and analyzed. These limitations highlight the need
for complementary approaches to event detection and response. Our proposed
method addresses these limitations by providing a reliable, real-time source
of information that is not dependent on individual contributions or vulnerable
to the same misinformation risks.

### Data Fusion for Event Inference

While there is a substantial body of research on the use of mobile data and
satellite imagery independently, there is a noticeable gap in the literature
when it comes to their combined use for real-time event identification and
response. This gap represents a significant opportunity for new research that
can contribute to both academic knowledge and practical applications in
disaster management, conflict monitoring, as well as search and rescue
operations.

The fusion of satellite imagery and GPS traces for event inference is a
relatively new area of research and is uniquely poised to overcome the
limitations placed on only analyzing each data source individually. Though GPS
traces have emerged as a prominent data source, it should be noted that such
datasets tend to exclude certain populations, such as the extreme poor
Pokhriyal et al. (2020). In developing countries, this population may
represent a significant share of the population. Thus, in the event of an
emergency, mobile data cannot be the sole indicator of irregular movement
patterns. Satellite imagery itself may also be inadequate in event inference
due to the lack of specific mobility information. However, when combined,
these two data sources can supplement each other and leverage the additional
information to better inform the inference task.

Though previous studies have explored the intersection of these two data
sources, the direction of application we intend to use our tool differs from
previous works. For instance, (Yin et al., 2021) proposed a method for fusing
satellite imagery and GPS traces to infer and fill in road attributes such as
speed limits and lane widths, with the goal of developing a more comprehensive
view of road networks globally. (Pokhriyal et al., 2020) fused satellite
images with mobile data to estimate poverty levels in Haiti, where a
nationwide census had not been conducted since 2012. In developing countries
where there may not be the economic resources necessary to complete annual
censuses, using a combination of data sources to estimate various baselines
may emerge as a crucial tool for financial aid allocation and future urban
planning. A combination of mobile data and satellite imagery has also been
used to study forests and obtain biomass measurements Molinier et al. (2016).
This study combined citizen science-powered mobile data collection and
satellite imagery to garner forest information and produce biomass maps in a
cost effective manner.

Notably, both the work of Yin et al. and Pokhriyal et al. provide information
that can be used for future planning. However, we hope to use GPS traces and
satellite imagery in real time for instantaneous data analysis and event
identification. Additionally, unlike in Molinier et al., where mobile data is
actively collected through an app that must be user-operated, we obtain GPS
traces passively through users’ interactions with various apps (not a specific
one). Thus, our event detection and inference capabilities are completely
powered by remote data, and without the need for active, real-time human
participation.

## Dataset and Software

For mobile data, we employ privacy-protected, passively generated GPS data
from an American data solution provider specializing in geospatial analytics.
The dataset contains discrete GPS points of 2,000 anonymous, opted-in
individuals in the Tulsa metropolitan area between December 2019 and July
2020. The locational data recorded includes timestamps, unique device
identifiers, latitudes and longitudes, and a measure of data precision (i.e.,
a spatial radius for which the provider has 95% confidence in the reported
coordinates). In addition to anonymizing the data, the data provider
obfuscates home locations to the census-block group level, and removes
sensitive points of interest from the dataset, in order to preserve privacy.
It’s also crucial to note that the temporal granularity of this data is
determined by user activity, meaning data points are generated irregularly and
not at fixed intervals. This non-uniformity of temporal data represents the
actual, inconsistent frequency of device usage among the sampled population.

For satellite imagery, we pull data from Planet Labs, a private American
company that specializes in Earth imaging. Planet operates a large
constellation of miniature satellites, known as Doves, which continually
capture and transmit images of Earth. The 3-band images we use have undergone
orthorectification, a process that corrects the geometric distortions in an
image to represent a flat surface, thereby ensuring that the scale is uniform
throughout. This makes the images reliable for distance measurements and for
overlaying with other geospatial data. Additionally, the images have been
color-corrected to ensure that the colors represented are as close as possible
to what the human eye would perceive if viewing the scene directly. This
enhances the clarity and interpretability of the images, making them a
valuable resource for our research.

We implemented our methodology in the Python programming language van Rossum
(1995). We used pandas Wes McKinney (2010) to read and wrangle mobile data and
Rasterio Gillies et al. (2013–) to manipulate satellite images. We also
leveraged NumPy Harris et al. (2020) for a variety of array operations,
GeoPandas Jordahl et al. (2019) for operations involving coordinates, and
scikit-mobility Pappalardo et al. (2019) to conduct mobile data cleaning and
preprocessing. Last but not least, we used matplotlib Hunter (2007) to produce
any visualizations of our results. Any other chunk of code we developed can be
found on our GitHub repository: https://github.com/ekinugurel/SatMobFusion.

## Methodology

Our methodology has two modules—one that deals with the analysis of mobile
data, and another with that of satellite imagery. In both modules, we discuss
relevant preprocessing methods as well as analysis steps to aid in the
inference task. Figure 2 shows the high-level steps in preprocessing and
analysis of both datasets while also outlining key geospatial analysis
processes for the three use cases discussed in the Introduction.

Figure 2. Methodological flowchart; the top bins describe the order of
processes employed in this paper, while the bottom outlines key geospatial
analysis metrics to carry out the inference task

### Mobile Data

For mobile data, the general preprocessing workflow is as follows:

  1. (1)

Spatial Partitioning

  2. (2)

Coordinate Reference System (CRS) Projection

  3. (3)

Data Filtering

To avoid wastefulness in computation, we begin by partitioning the mobile data
spatially, which we do by drawing a bounding box using four sets of
coordinates. In practice, there are several ways to achieve this at scale:
Leveraging a geocoding application protocol interface (API) to obtain the
coordinates based on an existing point of interest, partitioning by federal
and/or state designations such as census tracts, ZIP codes, or county borders,
manually typing in a bounding box (i.e., by defining the geographic
coordinates of the corners of the box), or obtaining a KMZ file from the
National Oceanic and Atmospheric Administration’s Damage Assessment Toolkit.
For this effort, we used the last option to derive the centroid coordinates of
the affected area around which we drew a $1\;\textrm{km}^{2}$ bounding square.

Once the data is spatially partitioned, we project coordinates to the
appropriate CRS. Different CRSs are designed to accurately represent different
aspects of the Earth’s surface. Some are designed to preserve area, while
others preserve shape, distance, or direction. In the context of event
inference, the optimal choice of CRS should minimize bias in the derived
mobility metrics, which may exist due to distortions of distances and areas.
Therefore, we use an equal-area projection around UTM Zone 15, which covers
the area around our case study.

We then filter out any erroneous data points based on a threshold of segment
velocity. Since this is done ”as the crow flies” and not with respect to the
underlying transportation network, this threshold should be set
conservatively. Velocity filtering discards infeasible jumps and oscillations
in trajectories, which may occur as a result of tall buildings blocking a
satellite’s line of sight to the mobile device (termed ”urban canyon effect”,
Ben-Moshe et al. (2011)). These oscillations may also occur due to low data
precision as a result of the ”cold-start problem”—when a signal dropped by
WiFi is not immediately picked up by a satellite, requiring some time before
being able to acquire navigation data and calculate a position Lehtinen et al.
(2008).

Once preprocessing is finished, we conduct inference-related analyses using
the following steps:

  1. (1)

Experimental Design Partitioning

  2. (2)

Exploratory Data Analysis (EDA)

  3. (3)

Statistical Testing

Specifically, we first define pre-, during-, and post-event periods that
temporally segment the data. The purpose of this partitioning is to derive key
metrics for each period, which would assist in detecting anomalies and enable
more formal statistical tests. The duration of these periods depends on the
nature of the analysis—for natural disasters like tornadoes or hurricanes, the
”during” period may be a few hours, while for applications in search & rescue
or in conflict areas, the same period may only last a few minutes.

In conducting EDA on time-series mobile data, the objectives are to understand
the structure of the data (i.e., identifying immediate trends, patterns, or
outliers), detect seasonality (i.e., periodic fluctuations), and prepare the
data for more formal anomaly detection techniques and statistical tests. The
EDA process can also assist in formulating hypotheses about possible causes
for anomalies by exploring relationships between the target variable and other
variables. In the context of mobile data within constrained geography, the
metrics we analyze include the radius of gyration $r_{g}$ (the characteristic
distance traveled by an individual in a given period, Equation 1) Song et al.
(2010), the number of extended stays at a particular location $c$ (Equation
3), and the number of visits per time unit $v$.

(1) $r_{g}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(r_{i}-r_{cm})^{2}}$

where $r_{i}$ represents the $i=1,\ldots,n$ locations recorded for the
individual and $r_{cm}=\frac{1}{n}\sum_{i=1}^{n}r_{i}$ is the center of mass
of the period’s trajectory.

(2) $\displaystyle c$ $\displaystyle=\sum_{i=1}^{n}c_{i},$ (3) $\displaystyle
c_{i}$ $\displaystyle=\begin{cases}1&\text{if }t_{i}=1,l_{i}=1\\\
0&\text{otherwise}\end{cases}$

where $t_{i}$ is an indicator variable denoting whether a device is observed
more than once within a 15-minute period and $l_{i}$ is an indicator variable
denoting whether multiple observations of the device are within a 100-meter
radius (only considered if $t_{i}=1$). We note that these thresholds are
robust and suggest adjusting them to suit the underlying domain.

(4) $v=\sum_{j=1}^{m}v_{j}$

where $j=1,\ldots,m$ indexes the subsections of the inspected geographical
area and $v_{j}$ is an indicator variable denoting whether the individual was
located in that subsection in the given period.

Finally, we conduct statistical tests to confirm the presence of anomalies. In
this paper, we employ the $Z$-score method (Equation 5), in which the
principle is to assume the data follows a normal distribution and then to
identify data points that are too far from the mean.

(5) $Z=\frac{x-\mu}{\sigma}$

where $x$ is a data point, $\mu$ is the mean of the dataset, and $\sigma$ is
the standard deviation. We consider Z-scores above 3 (or below -3) an anomaly.

### Satellite Imagery

The general workflow for processing and analyzing satellite images is as
follows:

  1. (1)

Specify tradeoff between region-of-interest (ROI) coverage and time gap to the
event of interest

  2. (2)

Pull imagery and project to CRS

  3. (3)

Image Analysis: RGB, greyscale, and NDVI differences

Providers like Planet Labs have multiple orbital satellites whose coverage
area may intersect—the goal of each satellite is to pass over a particular
geography at least once per day; however, various factors can influence the
availability of imagery and cause irregular gaps in observation, including
cloud cover, technical issues with onboard cameras or the image transmission
process, and regulatory restrictions imposed by nation-states. For these
reasons, we derive a utility metric (Equation 6) to determine the best
available before and after image for an event of interest. This metric weighs
the image’s area of coverage against its time gap with the event using a
tradeoff coefficient—the idea is to maximize the spatial overlap between the
image and the ROI while minimizing the number of days between the images.

(6) $u=\frac{A_{img}}{A_{evt}}-\frac{|t_{img}-t_{evt}|}{\phi}$

where the first term denotes the percent of ROI covered by the image (i.e.,
$A_{img}$ is the area of the image and $A_{evt}$ is the area affected by the
event), the numerator of the second term denotes the number of days between
the images, and $\phi$ is the coverage tradeoff coefficient. In this work, we
use $\phi=0.25$, meaning that 25% (percentage points) more coverage of the AOI
is defined to be worth as much as the image being taken 1 day before or after
the event.

As an additional filter to the above utility metric, we require all images to
have less than 50% cloud coverage. Once we have the satellite images, we
project them to the same CRS we used for the mobile data to avoid
misalignment.

We analyze the images to identify any physical changes to the built and
natural environment. This analysis involves looking at the regular visible
spectrum (RGB), a grayscale version (i.e., luminosity, Equation 7), and the
Normalized Difference Vegetation Index (NDVI, Equation 8), which measures the
greenness and the density of the vegetation captured in a satellite image.

(7) $\displaystyle NDVI$ $\displaystyle=\frac{NIR-R}{NIR+R}$ (8)
$\displaystyle GS$ $\displaystyle=0.299R+0.587G+0.114B$

where $NIR$, $R$, $G$, and $B$ are near-infrared, red, green, and blue
reflectances, respectively. NDVI is always between -1 and 1, where negative
values signify the presence of clouds, water, or snow, values close to zero
signify earthy rocks and bare soil, and positive values signify the existence
of some vegetation (with more positive meaning more green).

Figure 3 shows a decision tree outlining when the proposed method should be
employed based on a series of conditions, organized by the cognitive load
required to succeed in the inference task. Intuitively, crowdsourcing may be
the preferred option for most organizations and individuals but can be
unreliable due to the threat of misinformation or cybersecurity breaches.
Similarly, while cameras and other in-situ sensors could provide longer visual
feeds that are also easier to interpret, their availability in certain
settings tends to be sparse. If the above two conditions are true, and cloud-
free satellite imagery is available, then employing our method may be of
interest to organizations and individuals alike.

Figure 3. Usage decision tree for the proposed method

Relating to the above figure, we emphasize three points about the proposed
methodology:

  1. (1)

The process of event detection starts with mobile data, not satellite imagery,
due to the feasibility of passively monitoring geography based on mobility
metrics. However, the fusion of the two sources is what enables us to make an
inference.

  2. (2)

The proposed methodology is designed to be flexible, allowing the user to
define an event of interest using a configuration file.

  3. (3)

This approach can be used with both real-time and historical data, depending
on the context of the application, making it a potentially powerful tool for
improving emergency response times and information percolation speeds.

## Case Study

To validate our proposed methodology, we conducted a case study on the EF-1
tornado that hit near Muskogee, Oklahoma on May 15, 2020. This event provided
an opportunity to test our approach in a real-world scenario involving a
small-scale disaster in a rural region.

After spatially partitioning the GPS traces to include the subset within the
bounding box of the tornado, we projected all data to UTM Zone 15 and filtered
the data by the logic described in Section Methodology. Then, we defined the
two weeks before the day of the event (May 2 - May 15) as the before, the day
of the event (May 15) as the during, and the week after the event (May 15 -
May 22) as the after periods and began exploring the data. Figure 4 visualizes
the time series of two metrics, highlighting the periods by color; there are
no noticeable changes in the radius of gyration, while the number of extended
stays at a given location within the ROI increases sharply the day after the
event. However, Z-score tests on both metrics do not detect any anomalies.

(a) Radius of Gyration

(b) Number of Stays per day

Figure 4. EDA of mobility metrics in the before, during, and after periods

However, when we look at the third metric (visits per hour to the ROI, Figure
5), we see an unmistakable outlier; on the day of the tornado, visits to the
area increased sharply at 2 pm and again at 6 pm, substantially beyond the
usual fluctuations captured by the 95% confidence interval (light blue shade)
in the before and after periods. Furthermore, Z-score tests on this metric
confirm the hypothesis: Two anomalies exist in the data and they are at 2 pm
($Z=3.47$) and at 6 pm ($Z=4.87$) on May 15th.

Figure 5. Visits per hour to the ROI in the before, during, and after periods
sorted by hour of day

Following the anomaly detection, we pull in satellite imagery as described in
the previous section. We use the utility metric in Equation 6 to land on two
images, one from May 9 and another from May 18. These two images are free of
clouds and strike a balance between the percent of ROI covered and the number
of days separating them from the event.

We analyze the images based on visual inspection, the greyscale difference,
and the NDVI difference (Figure 6). We highlight two spots on the regular
color band images in which changes to the soil and a corner building are
apparent to the naked eye. Furthermore, while noisy, the greyscale difference
image shows changes to luminosity and the NDVI difference image shows loss of
vegetation throughout the captured area (note that 0.0 is shown in light
green).

The images show clear evidence of physical changes to the ROI. When combined
with insights from GPS traces, and even without hindsight bias, we would infer
that a mildly-destructive natural disaster occurred, piquing passerby’s
interest as well as landowners’ concern on the day of the event (resulting in
greater visits to the ROI).

(a) Before the Muskogee tornado on May 9, 2020

(b) After the Muskogee tornado on May 18, 2020

(c) Greyscale difference between the two images

(d) NDVI difference between the two images

Figure 6. Satellite imagery with 10m resolution documenting changes to the
built and natural environment after the Muskogee tornado

## Use Cases

The proposed methodology is expected to outperform existing methods for data
collection and information percolation for a range of use cases. We especially
anticipate broad applicability in the following settings:

  * •

Intelligent transportation systems: The likelihood of finding parking
approximate to a desired location and in a timely manner has been shown to be
a determinant of travelers’ mode choice Assemi et al. (2020). Therefore,
providing accurate parking-related information can augment an individual’s
decision-making process. Public and private entities can leverage the proposed
method to provide real-time information on parking availability. The temporal
sparsity of satellite imagery can be augmented with GPS traces to provide an
online learning framework for the likelihood of finding parking.

  * •

Events in rural areas: In rural areas where traditional data collection
methods may be sparse or non-existent, our methodology could provide valuable
insights into events as they unfold. The combination of satellite imagery and
mobile data could help identify and respond to emergencies such as wildfires,
floods, or lost hikers.

  * •

Extreme-weather events: During extreme-weather events like hurricanes,
tornadoes, or blizzards, our methodology could provide real-time information
about the event’s impact and aid in coordinating response efforts. The ability
to quickly identify areas of high impact could help prioritize resource
allocation and rescue operations.

  * •

War-torn states: In regions affected by conflict, obtaining reliable on-the-
ground information can be challenging. Satellite imagery combined with mobile
data could help identify areas of active conflict, population displacement, or
infrastructure damage.

  * •

Areas without network connectivity: In areas where network connectivity is
limited or non-existent, such as remote wilderness areas or during large-scale
network outages, satellite imagery can still provide valuable data. When
combined with historical or intermittent mobile data, this could help identify
and respond to events such as lost hikers or the spread of wildfires.

  * •

National security matters: Our methodology could also be used in the context
of national security, helping to identify potential threats or unusual
activities. The fusion of satellite imagery with mobile data could provide a
more comprehensive picture of activities across a wide geographical area.

In each of these use cases, our methodology aims to provide a more nuanced and
timely understanding of events as they unfold, potentially improving response
times and outcomes.

## Conclusion

We have presented a novel data fusion methodology combining satellite imagery
with GPS traces generated from mobile devices. Our approach leverages the
strengths of both data types: mobile data provides an approximation of human
mobility, proximity to one another, and the built environment, while satellite
imagery offers visual insights into physical changes to the built and natural
environment. Furthermore, we have demonstrated the case study of a small-scale
tornado, showing the promise of our method in identifying anomalies in
movement patterns using mobile data, changes to the built environment using
satellite imagery, and inferring the cause of both using our fusion
methodology.

Our methodology is viable due to the increasing accessibility of commercial
data providers specializing in remote sensing and mobile devices. This trend
is only expected to continue—the number of active orbital satellites has grown
exponentially since 2017 orb (2022), while the percentage of Americans that
own a smartphone has risen to 97% mob (2021). The continued commercialization
of such products will enable researchers and practitioners to monitor,
categorize, and understand the natural and built environments in ways
unimaginable just 50 years ago. We hope our method is a step in this
direction, enabling others to pursue endeavors made possible by the fusion of
these data sources.

## Data Access Statement

The authors note that only Ekin Ugurel had access to the privacy-protected GPS
dataset. Data supporting this study cannot be made available due to privacy-
preserving research agreements with the data provider.

## Author Contributions

The authors confirm contribution to the paper as follows: study conception and
design: E. Ugurel, S. Coenen; data collection: E. Ugurel, S. Coenen; analysis
and interpretation of results: E. Ugurel, S. Coenen; draft manuscript
preparation: E. Ugurel, M. Chen, S. Coenen, C. Chen. All authors reviewed the
results and approved the final version of the manuscript.

## References

  * mob [2021] Mobile fact sheet.  Technical report, Pew Reseearch Center, 2021.  URL https://www.pewresearch.org/internet/fact-sheet/mobile/.
  * orb [2022] Online index of objects launched into outer space.  Technical report, United Nations - Office of Outer Space Affairs, 2022\.
  * Albert et al. [2017] A. Albert, J. Kaur, and M. Gonzalez.  Using convolutional networks and satellite imagery to identify patterns in urban environments at a large scale.  pages 1357–1366, 08 2017.  doi: 10.1145/3097983.3098070.
  * Alessandretti [2022] L. Alessandretti.  What human mobility data tell us about COVID-19 spread.  _Nature Reviews Physics_ , 4(1):12–13, Jan. 2022\.  ISSN 2522-5820.  doi: 10.1038/s42254-021-00407-1.  URL https://www.nature.com/articles/s42254-021-00407-1.  Number: 1 Publisher: Nature Publishing Group.
  * Allcott et al. [2019] H. Allcott, M. Gentzkow, and C. Yu.  Trends in the diffusion of misinformation on social media.  _Research & Politics_, 6(2):2053168019848554, 2019.
  * Ashktorab et al. [2014] Z. Ashktorab, C. Brown, M. Nandi, and A. Culotta.  Tweedr: Mining twitter to inform disaster response.  In _International Conference on Information Systems for Crisis Response and Management_ , 2014.
  * Assemi et al. [2020] B. Assemi, D. Baker, and A. Paz.  Searching for on-street parking: An empirical investigation of the factors influencing cruise time.  _Transport Policy_ , 97:186–196, 2020.
  * Ben-Moshe et al. [2011] B. Ben-Moshe, E. Elkin, H. Levi, and A. Weissman.  Improving Accuracy of GNSS Devices in Urban Canyons.  page 6, 2011.
  * Bengtsson et al. [2011] L. Bengtsson, X. Lu, A. Thorson, R. Garfieldarfier, and J. Schreeb.  Improved response to disasters and outbreaks by tracking population movements with mobile phone network data: A post-earthquake geospatial study in haiti.  _PLoS medicine_ , 8:e1001083, 08 2011.  doi: 10.1371/journal.pmed.1001083.
  * Bewley [2003] R. Bewley.  Aerial survey for archaeology.  _The Photogrammetric Record_ , 18:273 – 292, 12 2003.  doi: 10.1046/j.0031-868X.2003.00023.x.
  * Calabrese et al. [2010] F. Calabrese, M. Colonna, P. Lovisolo, D. Parata, and C. Ratti.  Real-time urban monitoring using cell phones: A case study in rome.  _IEEE transactions on intelligent transportation systems_ , 12(1):141–151, 2010.
  * Castro et al. [2014] P. Castro, D. Zhang, C. Chen, S. Li, and G. Pan.  From taxi gps traces to social and community dynamics: A survey.  _ACM Computing Surveys_ , 46, 01 2014.  doi: 10.1145/2543581.2543584.
  * Chen et al. [2014] C. Chen, L. Bian, and J. Ma.  From traces to trajectories: How well can we guess activity locations from mobile phone traces?  _Transportation Research Part C: Emerging Technologies_ , 46:326–337, 2014.
  * Estellés-Arolas and L. Guevara [2012] E. Estellés-Arolas and F. G. L. Guevara.  Towards an integrated crowdsourcing definition.  _Journal of Information Science_ , 38, 04 2012.  doi: 10.1177/0165551512437638.
  * Gillies et al. [2013–] S. Gillies et al.  Rasterio: geospatial raster i/o for Python programmers, 2013–.  URL https://github.com/rasterio/rasterio.
  * Goodchild [2007] M. F. Goodchild.  Citizens as sensors: the world of volunteered geography.  _GeoJournal_ , 69:211–221, 2007.
  * Guan and Chen [2014] X. Guan and C. Chen.  Using social media data to understand and assess disasters.  _Natural hazards_ , 74:837–850, 2014.
  * Haddad [2011] N. Haddad.  From ground surveying to 3d laser scanner: A review of techniques used for spatial documentation of historic sites.  _Journal of King Saud University - Engineering Sciences_ , 23:109–118, 06 2011.  doi: 10.1016/j.jksues.2011.03.001.
  * Harris et al. [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, Sept. 2020.  doi: 10.1038/s41586-020-2649-2.  URL https://doi.org/10.1038/s41586-020-2649-2.
  * Herrera et al. [2010] J. C. Herrera, D. B. Work, R. Herring, X. J. Ban, Q. Jacobson, and A. M. Bayen.  Evaluation of traffic data obtained via gps-enabled mobile phones: The mobile century field experiment.  _Transportation Research Part C: Emerging Technologies_ , 18(4):568–583, 2010.
  * Hong et al. [2008] S.-J. Hong, J. Thong, J. Y. Moon, and K. Tam.  Understanding the behavior of mobile data services consumers.  _Information Systems Frontiers_ , 10:431–445, 09 2008.  doi: 10.1007/s10796-008-9096-1.
  * Hunter [2007] J. D. Hunter.  Matplotlib: A 2d graphics environment.  _Computing in Science & Engineering_, 9(3):90–95, 2007.  doi: 10.1109/MCSE.2007.55.
  * Ji et al. [2018] M. Ji, L. Liu, and M. Buchroithner.  Identifying collapsed buildings using post-earthquake satellite imagery and convolutional neural networks: A case study of the 2010 haiti earthquake.  _Remote Sensing_ , 10(11):1689, 2018.
  * Jongman et al. [2015] B. Jongman, J. Wagemaker, B. Revilla Romero, and E. Coughlan de Perez.  Early flood detection for rapid humanitarian response: harnessing near real-time satellite and twitter signals.  _ISPRS International Journal of Geo-Information_ , 4(4):2246–2266, 2015.
  * Jordahl et al. [2019] K. Jordahl, J. V. den Bossche, J. Wasserman, J. McBride, J. Gerard, M. Fleischmann, J. Tratner, M. Perry, C. Farmer, G. A. Hjelle, S. Gillies, M. Cochran, M. Bartos, L. Culbertson, N. Eubank, maxalbert, S. Rey, A. Bilogur, D. Arribas-Bel, C. Ren, J. Wilson, M. Journois, L. J. Wolf, L. Wasser, Ömer Özak, YuichiNotoya, F. Leblanc, Filipe, C. Holdgraf, and A. Greenhall.  geopandas/geopandas: v0.6.1, Oct. 2019.  URL https://doi.org/10.5281/zenodo.3483425.
  * Karunanithi et al. [2016] M. Karunanithi, S. Kumar, and P. Ponnusamy.  Use of high resolution google earth satellite imagery in landuse map preparation for urban related applications.  _Procedia Technology_ , 24:1835–1842, 12 2016.  doi: 10.1016/j.protcy.2016.05.231.
  * Kirillov et al. [2023] A. Kirillov, E. Mintun, N. Ravi, H. Mao, C. Rolland, L. Gustafson, T. Xiao, S. Whitehead, A. C. Berg, W.-Y. Lo, P. Dollár, and R. Girshick.  Segment anything.  _arXiv:2304.02643_ , 2023.
  * Lehtinen et al. [2008] M. Lehtinen, A. Happonen, and J. Ikonen.  Accuracy and time to first fix using consumer-grade gps receivers.  In _2008 16th International Conference on Software, Telecommunications and Computer Networks_ , pages 334–340, 2008.  doi: 10.1109/SOFTCOM.2008.4669506.
  * Liao et al. [2007] L. Liao, D. Fox, and H. Kautz.  Extracting places and activities from gps traces using hierarchical conditional random fields.  _The International Journal of Robotics Research_ , 26(1):119–134, 2007.
  * Molinier et al. [2016] M. Molinier, C. López-Sánchez, T. Toivanen, I. Korpela, J. J. Corral-Rivas, R. Tergujeff, and T. Häme.  Relasphone—mobile and participative in situ forest biomass measurements supporting satellite image mapping.  _Remote Sensing_ , 8:869–891, 10 2016.  doi: 10.3390/rs8100869.
  * Okolloh [2009] O. Okolloh.  Ushahidi, or ‘testimony’: Web 2.0 tools for crowdsourcing crisis information.  _Participatory learning and action_ , 59(1):65–70, 2009.
  * Palen [2008] L. Palen.  Online social media in crisis events.  _Educause Quarterly_ , 31:76–78, 01 2008.
  * Palen et al. [2007] L. Palen, S. Vieweg, J. Sutton, S. B. Liu, and A. Hughes.  Crisis informatics: Studying crisis in a networked world.  In _Proceedings of the Third International Conference on E-Social Science_ , pages 7–9, 2007.
  * Pan et al. [2013] G. Pan, G. Qi, Z. Wu, D. Zhang, and S. Li.  Land-use classification using taxi gps traces.  _Intelligent Transportation Systems, IEEE Transactions on_ , 14:113–123, 03 2013.  doi: 10.1109/TITS.2012.2209201.
  * Pappalardo et al. [2019] L. Pappalardo, F. Simini, G. Barlacchi, and R. Pellungrini.  scikit-mobility: A python library for the analysis, generation and risk assessment of mobility data.  _arXiv preprint arXiv:1907.07062_ , 2019.
  * Pokhriyal et al. [2020] N. Pokhriyal, O. Zambrano, J. Linares, and H. Hernández.  _Estimating and Forecasting Income Poverty and Inequality in Haiti Using Satellite Imagery and Mobile Phone Data_.  06 2020.  doi: 10.18235/0002466.
  * Said et al. [2019] N. Said, K. Ahmad, M. Riegler, K. Pogorelov, L. Hassan, N. Ahmad, and N. Conci.  Natural disasters detection in social media and satellite imagery: a survey.  _Multimedia Tools and Applications_ , 78:31267–31302, 2019\.
  * Sakaki et al. [2010] T. Sakaki, M. Okazaki, and Y. Matsuo.  Earthquake shakes twitter users: real-time event detection by social sensors.  In _Proceedings of the 19th international conference on World wide web_ , pages 851–860, 2010.
  * Song et al. [2010] C. Song, Z. Qu, N. Blumm, and A.-L. Barabási.  Limits of predictability in human mobility.  _Science_ , 327(5968):1018–1021, 2010.
  * Stelter and Cohen [2008] B. Stelter and N. Cohen.  Citizen journalists provided glimpses of mumbai attacks.  _The New York Times_ , 2008.  URL https://www.nytimes.com/2008/11/30/world/asia/30twitter.html#:~:text=From%20his%20terrace%20on%20Colaba,every%20move%20on%20the%20Internet.
  * Sulis et al. [2018] P. Sulis, E. Manley, C. Zhong, and M. Batty.  Using mobility data as proxy for measuring urban vitality.  _Journal of Spatial Information Science_ , 2018(16):137–162, 2018.
  * Tian et al. [2017] L. Tian, T. Yao, Y. Gao, L. Thompson, E. Mosley-Thompson, S. Muhammad, J. Zong, C. Wang, S. Jin, Z. Li, and et al.  Two glaciers collapse in western tibet.  _Journal of Glaciology_ , 63(237):194–197, 2017\.  doi: 10.1017/jog.2016.122.
  * van Rossum [1995] G. van Rossum.  Python tutorial.  Technical Report CS-R9526, Centrum voor Wiskunde en Informatica (CWI), Amsterdam, May 1995.
  * Voigt et al. [2007] S. Voigt, T. Kemper, T. Riedlinger, R. Kiefl, K. Scholte, and H. Mehl.  Satellite image analysis for disaster and crisis-management support.  _IEEE transactions on geoscience and remote sensing_ , 45(6):1520–1528, 2007.
  * Wes McKinney [2010] Wes McKinney.  Data Structures for Statistical Computing in Python.  In Stéfan van der Walt and Jarrod Millman, editors, _Proceedings of the 9th Python in Science Conference_ , pages 56 – 61, 2010.  doi: 10.25080/Majora-92bf1922-00a.
  * Wu and Osco [2023] Q. Wu and L. P. Osco.  samgeo: A Python package for segmenting geospatial data with the Segment Anything Model (SAM), May 2023.  URL https://doi.org/10.5281/zenodo.7966658.
  * Xu [2010] H. Xu.  Analysis of impervious surface and its impact on urban heat environment using the normalized difference impervious surface index (ndisi).  _Photogrammetric Engineering & Remote Sensing_, 76(5):557–565, 2010.
  * Yin et al. [2021] Y. Yin, A. Tran, Y. Zhang, W. Hu, G. Wang, J. Varadarajan, R. Zimmermann, and S.-K. Ng.  Multimodal fusion of satellite images and crowdsourced gps traces for robust road attribute detection.  In _Proceedings of the 29th International Conference on Advances in Geographic Information Systems_ , pages 107–116, 2021.
  * Zhang et al. [2016] L. Zhang, L. Zhang, and B. Du.  Deep learning for remote sensing data: A technical tutorial on the state of the art.  _IEEE Geoscience and remote sensing magazine_ , 4(2):22–40, 2016.
  * Zheng et al. [2008] Y. Zheng, Q. Li, Y. Chen, X. Xie, and W.-Y. Ma.  Understanding mobility based on gps data.  In _Proceedings of the 10th international conference on Ubiquitous computing_ , pages 312–321, 2008.


# Deep-learning-based Early Fixing for Gas-lifted Oil Production Optimization:
Supervised and Weakly-supervised Approaches

Bruno M. Pacheco  Laio O. Seman  Eduardo Camponogara  Department of Automation
and Systems Engineering, Federal University of Santa Catarina (UFSC),
Florianópolis, Brazil (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>eduardo.camponogara@ufsc.br).

###### Abstract

Maximizing oil production from gas-lifted oil wells entails solving Mixed-
Integer Linear Programs (MILPs). As the parameters of the wells, such as the
basic-sediment-to-water ratio and the gas-oil ratio, are updated, the problems
must be repeatedly solved. Instead of relying on costly exact methods or the
accuracy of general approximate methods, in this paper, we propose a tailor-
made heuristic solution based on deep learning models trained to provide
values to all integer variables given varying well parameters, early-fixing
the integer variables and, thus, reducing the original problem to a linear
program (LP). We propose two approaches for developing the learning-based
heuristic: a supervised learning approach, which requires the optimal integer
values for several instances of the original problem in the training set, and
a weakly-supervised learning approach, which requires only solutions for the
early-fixed linear problems with random assignments for the integer variables.
Our results show a runtime reduction of 71.11% Furthermore, the weakly-
supervised learning model provided significant values for early fixing,
despite never seeing the optimal values during training.

###### keywords:

Mixed-integer optimization, Deep learning, Weakly-supervised learning, Early
fixing, Oil production systems

††thanks: This research was partly funded by Petróleo Brasileiro S.A. under
grant 2018/00217-3, FAPESC under grant 2021TR2265, CNPq and CAPES.

##  1 Introduction

The maximization of oil production in an offshore platform is a challenging
problem due to the physical models’ complexity and various operational
constraints (Luguesi et al., 2023; Müller et al., 2022). Furthermore, the
variations in the oil wells and the multitude of technologies that can be
employed account for the many configurations this problem presents. Therefore,
it becomes necessary to rely on optimization models to reach optimal
conditions for the operation.

Due to the nonlinearities of the oil output stream from gas-lifted wells, the
formulation of the optimization model is not straightforward. The approach of
Müller et al. (2022) uses a piecewise linear model of the wells and formulates
the problem as a Mixed-Integer Linear Program (MILP). More specifically, the
relationship between the liquid production of the wells and the lift-gas
injection, together with the wellhead pressure, is modeled as a piecewise
linear function. The piecewise-linear functions are defined by studying the
behavior of the relations through simulation runs or field tests and defining
breakpoints between which the relation is considered to be linear. Special
Ordered Set of type 2 (SOS2) constraints are then used to define the
operational region for each constraint, which defines the active boundaries
for a region in the state space in which the function approximation is linear.
The final model is linear with binary variables and SOS2 variable sets.

Solving MILPs exactly requires algorithms such as branch-and-bound and branch-
and-cut (Lee and Mitchell, 2009). These algorithms usually require long
runtimes due to the many iterations of solving linear relaxations of the
original problem. One can use approximate methods instead of relying solely on
exact algorithms to optimize oil production. Such methods can improve the
efficiency of the optimization process but may provide suboptimal solutions.
Derivative-free methods like genetic algorithms, simulated annealing, and
particle swarm optimization can also be used to find near-optimal solutions
efficiently (Seman et al., 2020). Another approximate method is to early fix
the variables, reducing the dimension of the problem. In the case of MILPs,
one can develop a heuristic to fix all integer variables, reducing the problem
to a linear program, which can be solved very efficiently with algorithms such
as the simplex. As pointed out by Bengio et al. (2021), such heuristics can be
very hard to handcraft, which makes machine learning (ML) models natural
candidates for the task.

In this paper, we propose two deep-learning approaches to reduce the runtime
of the gas-lifted oil production MILP through early fixing. One is a
supervised learning approach that requires that the optimal set of binary
variables is known for each problem instance; instances of the problem must be
solved exactly to build the training data. With this data, the model is then
trained to provide the optimal binary variables given the parameters of the
MILP. The other is a weakly-supervised approach that requires just the
solutions to the early-fixed problem, i.e., after fixing the integer variables
to any (binary) value. In other words, the training data is generated from
solving linear problems, resulting from randomly fixing the integer variables
in the original MILP. To the best of our knowledge, this is the first work to
successfully implement learning-based heuristics to speed up the solution of
oil production optimization.

###  1.1 Related work

In this section, we provide a brief overview of related work in approximate
methods for MILP, focusing on deep-learning-based methods.

The use of heuristics in MILP solvers is common. For instance, the SCIP solver
(Vigerske and Gleixner, 2018) uses primal heuristics to find feasible
solutions. However, recent studies have highlighted the development of
heuristics based on ML techniques. Particularly, Bengio et al. (2021) have
suggested that the potential benefits of using ML include reduced
computational time and improved solution quality.

Ding et al. (2019) presented a learning-based heuristic to accelerate the
solution-finding process. The authors propose a Graph Neural Network (GNN)
model that predicts solution values for branching. The developed heuristic is
used to guide a branch-and-bound tree search.

Li and Wu (2022) train a learning-based heuristic for early fixing MILPs
within an Alternating Direction Method of Multipliers (ADMM) algorithm. By
formulating the early fixing as a Markov decision process, the authors use
reinforcement learning to train the heuristic. The authors showed that the
proposed heuristic significantly improves the solving speed and can even
improve the solution quality.

Pacheco et al. (2023) explore using GNNs as early-fixing heuristics for the
Offline Nanosatellite Task Scheduling (ONTS) problem. In this direction, the
authors implement a GNN-based solution that uses bipartite graphs, feature
aggregation, and message-passing techniques. The results suggest that GNNs can
be an effective method for aiding in the decision-making process of MILPs.

Finally, Anderson et al. (2022) proposed a weakly-supervised approach for
warm-starting gas network problems. The authors present a model that generates
feasible initial solutions which are used to accelerate the MILP solver’s
convergence. The results show a 60.5% decrease in the runtime.

In summary, the literature offers promising approaches for accelerating MILP
solvers using ML techniques, particularly through early fixing. However, to
the best of our knowledge, no one has applied such techniques to oil
production optimization. In this context, this paper is the first to explore
the use of supervised and weakly-supervised learning approaches to the oil
production maximization problem, leveraging surrogate models for the liquid
production and offering insights into this growing area that hybridizes
optimization and deep learning.

##  2 Problem Statement

In this section, we present the problem formulation, which is based on Müller
et al. (2022), with only gas–lifted oil wells connected to manifolds.
Nevertheless, it is easy to see that our methodological approach can be
applied to variations of this problem, e.g., platforms with satellite wells,
wells with electrical submersible pump systems, and subsea manifolds.

###  2.1 Well model

A single production platform can extract from multiple oil wells. Each well
$n\in\mathcal{N}$ has its liquid production $q_{\rm l}^{n}$ induced by the
wellhead pressure ${\rm whp}^{n}$ and a lift-gas injection rate $q_{\rm
gl}^{n}$. The relationship between $q_{\rm l}^{n}$, ${\rm whp}^{n}$, and
$q_{\rm gl}^{n}$ is modeled based on the natural characteristics of the well
and the gas-oil ratio (GOR, referred to as ${\rm gor}^{n}$) and basic-
sediment-to-water ratio (BSW, referred to as ${\rm bsw}^{n}$) of its liquid
production. Both ${\rm bsw}^{n}$ and ${\rm gor}^{n}$ are measured through
separation tests and are considered static during the optimization. However,
as they change with time, their updates drive new executions of the
optimization process to keep the results reliable.

###  2.2 Piecewise linearization

We use the Marlim simulator (Seman et al., 2020), a proprietary software from
Petrobras, to model the liquid output of each well,

$q_{\rm l}^{n}=\textrm{Marlim}(q_{\rm gl}^{n},{\rm whp}^{n},{\rm bsw}^{n},{\rm
gor}^{n}).$ (1)

An example of the liquid flow function of a real well can be seen in Figure 1
with fixed values for ${\rm bsw}^{n}$ and ${\rm gor}^{n}$. As both $q_{\rm
gl}^{n}$ and ${\rm whp}^{n}$ can be controlled, but have a nonlinear
relationship with the outcome, we apply piecewise linearization to $q_{\rm
l}^{n}$ as a function of both. More precisely, let $\mathcal{K}_{\rm
gl}=\\{1,\ldots,k_{\rm gl}\\}$ and $\mathcal{J}_{\rm whp}=\\{1,\ldots,j_{\rm
whp}\\}$ be sets of indices for lift-gas injection and wellhead pressure. Let
also $\mathcal{Q}_{\rm gl}^{n}=\\{q_{\rm
gl}^{n,k}\in\mathbb{R}_{+}:k\in\mathcal{K}_{\rm gl}\\}$ and
$\mathcal{W}^{n}_{\rm whp}=\\{{\rm
whp}^{n,j}\in\mathbb{R}_{+}:j\in\mathcal{J}_{\rm whp}\\}$ be the respective
breakpoint values for well $n$, and

$\mathcal{Q}_{\rm liq}^{n}=\left\\{q_{\rm l}^{n,k,j}:q_{\rm
l}^{n,k,j}=\hat{q}_{\rm l}^{n}(q_{\rm gl}^{n,k},{\rm whp}^{n,j},{\rm
bsw}^{n},{\rm gor}^{n}),\right.\\\ \left.{k\in\mathcal{K}_{\rm
gl},j\in\mathcal{J}_{\rm whp}}\right\\}$ (2)

be the liquid flow rates for the well at the breakpoints. $\mathcal{Q}_{\rm
liq}$ is obtained by using Marlim as in Eq. (1) with the adjusted parameters
${\rm gor}^{n}$ and ${\rm bsw}^{n}$.

The piecewise approximation is then given by

$\displaystyle q_{\rm gl}^{n}~{}$ $\displaystyle=~{}\sum_{k\in\mathcal{K}_{\rm
gl}}\sum_{j\in\mathcal{J}_{\rm whp}}\theta^{n}_{k,j}q_{\rm gl}^{n,k}$ (3a)
$\displaystyle{\rm whp}^{n}~{}$ $\displaystyle=~{}\sum_{k\in\mathcal{K}_{\rm
gl}}\sum_{j\in\mathcal{J}_{\rm whp}}\theta^{n}_{k,j}{\rm whp}^{n,j}$ (3b)
$\displaystyle q_{\rm l}^{n}~{}$ $\displaystyle=~{}\sum_{k\in\mathcal{K}_{\rm
gl}}\sum_{j\in\mathcal{J}_{\rm whp}}\theta^{n}_{k,j}q_{\rm l}^{n,k,j}$ (3c)
$\displaystyle 1~{}$ $\displaystyle=~{}\sum_{k\in\mathcal{K}_{\rm
gl}}\sum_{j\in\mathcal{J}_{\rm whp}}\theta^{n}_{k,j}$ (3d)
$\displaystyle\theta^{n}_{k,j}~{}$ $\displaystyle\geq
0,~{}k\in\mathcal{K}_{\rm gl},\,j\in\mathcal{J}_{\rm whp}.$ (3e)

To ensure piecewise linearization, SOS2 constraints are necessary for the
values of $\theta$ that correspond to $q_{\rm gl}$ and ${\rm whp}$, as follows

$\displaystyle\eta_{k}^{n}~{}$ $\displaystyle=~{}\sum_{j\in\mathcal{J}_{\rm
whp}}\theta^{n}_{k,j},\,k\in\mathcal{K}_{\rm gl}$ (4a)
$\displaystyle\eta_{j}^{n}~{}$ $\displaystyle=~{}\sum_{k\in\mathcal{K}_{\rm
gl}}\theta^{n}_{k,j},\,j\in\mathcal{J}_{\rm whp}$ (4b)
$\displaystyle\left\\{\eta_{k}^{n}\right\\}_{k\in\mathcal{K}_{\rm
gl}}~{}~{}\text{is SOS2}$ (4c)
$\displaystyle\left\\{\eta_{j}^{n}\right\\}_{j\in\mathcal{J}_{\rm
whp}}\text{is SOS2}.$ (4d)

SOS2 constraints imply that no more than two consecutive elements of the
ordered set of variables are nonzero (Beale and Tomlin, 1969).

###  2.3 Problem formulation

All oil wells are connected to a hub that directs their liquid productions
$q_{\rm l}^{n}$ to separators. The separators then split the liquid flow into
the oil, gas, and water phase flows. This separation depends on the BSW and
GOR of each well $n\in\mathcal{N}$,

$\displaystyle q_{\rm oil}^{n}$ $\displaystyle=q_{\rm l}^{n}\cdot(1-{\rm
bsw}^{n})$ (5a) $\displaystyle q_{\rm water}^{n}$ $\displaystyle=q_{\rm
l}^{n}\cdot{\rm bsw}^{n}$ (5b) $\displaystyle q_{\rm gas}^{n}$
$\displaystyle=q_{\rm l}^{n}\cdot(1-{\rm bsw}^{n})\cdot{\rm gor}^{n}.$ (5c)

The total gas flow is limited by the maximum lift-gas flow available
$\bar{q}_{\rm gl}$ though

$\sum_{n\in\mathcal{N}}q_{\rm gl}^{n}\leq\bar{q}_{\rm gl}.$ (6)

Finally, the objective is to maximize the total volume of oil extracted

$\sum_{n\in\mathcal{N}}q_{\rm oil}^{n}.$ (7)

We can express the problem as

$\begin{split}\max_{\boldsymbol{q}_{\rm gl},\boldsymbol{{\rm
whp}}}~{}&\eqref{eq:objective}\\\ \textrm{s.t.}~{}~{}&\eqref{eq:pwl-
approximation}\textrm{--}\eqref{eq:q-gl-max},\end{split}$ (8)

where $\boldsymbol{q}_{\rm gl},\boldsymbol{{\rm
whp}}\in\mathbb{R}^{|\mathcal{N}|}$ are the flow of lift-gas injected and the
well-head pressure of each well.

Note that the problem, in this case using piecewise linearization of the
liquid flow, is an MILP parameterized by ${\rm bsw}^{n}$, ${\rm gor}^{n}$,
$\bar{q}_{\rm gl}$, and $Q_{\rm liq}^{n}$. Therefore, let $\pi\in\Pi$ be the
vector of problem parameters. We can write the problem as

$\begin{split}P(\pi)=\max_{x,z}~{}&c^{T}x\\\
\text{s.t.}~{}&A(\pi)x+C(\pi)z\leq b(\pi)\\\ &x\in
X\subset\mathbb{R}^{n_{x}},z\in Z\subset\\{0,1\\}^{n_{z}}\end{split}$ (9)

where $x$ is the vector of continuous variables (e.g., $q_{\rm l}^{n}$,
$q_{\rm gl}^{n}$, $\theta^{n}_{k,j}$) and $z$ is the vector of binary
variables necessary for the SOS2 constraints. More precisely,
$z=(\boldsymbol{z}_{\rm gl},\boldsymbol{z}_{{\rm whp}})$ is such that, for all
$n\in\mathcal{N}$,

$\displaystyle\eta^{n}_{k}$ $\displaystyle\leq z_{\rm gl}^{n,k},k=1$ (10a)
$\displaystyle\eta^{n}_{j}$ $\displaystyle\leq z_{{\rm whp}}^{n,j},j=1$ (10b)
$\displaystyle\eta^{n}_{k}$ $\displaystyle\leq\max(z_{\rm gl}^{n,k},z_{\rm
gl}^{n,k-1}),k\in\mathcal{K}_{\rm gl}\setminus\\{1,k_{\rm gl}\\}$ (10c)
$\displaystyle\eta^{n}_{j}$ $\displaystyle\leq\max(z_{{\rm whp}}^{n,j},z_{{\rm
whp}}^{n,j-1}),j\in\mathcal{J}_{{\rm whp}}\setminus\\{1,j_{{\rm whp}}\\}$
(10d) $\displaystyle\eta^{n}_{k}$ $\displaystyle\leq z_{\rm
gl}^{n,k-1},k=k_{\rm gl}$ (10e) $\displaystyle\eta^{n}_{j}$ $\displaystyle\leq
z_{{\rm whp}}^{n,j-1},j=j_{{\rm whp}}$ (10f)
$\displaystyle\sum_{k\in\mathcal{K}_{\rm gl}\setminus\\{k_{\rm gl}\\}}z_{\rm
gl}^{n,k}=1,\sum_{j\in\mathcal{J}_{{\rm whp}}\setminus\\{j_{{\rm
whp}}\\}}z_{{\rm whp}}^{n,j}=1$ (10g) $\displaystyle z_{\rm
gl}^{n}\in\\{0,1\\}^{k_{\rm gl}-1}$ (10h) $\displaystyle z_{{\rm
whp}}^{n}\in\\{0,1\\}^{j_{{\rm whp}}-1}$ (10i)

in which $z_{\rm gl}^{n,k}=1$ if the $k$-th interval is selected, meaning that
$q_{\rm gl}^{n}\in[q_{\rm gl}^{n,k},q_{\rm gl}^{n,k+1}]$ and only
$\eta^{k}_{n}$ and $\eta^{k+1}_{n}$ can be nonzero, otherwise $z_{\rm
gl}^{n,k}$ assumes value $0$. The semantics of the binary variables $z_{\rm
whp}^{n,j}$ are analogous. Notice that, because $z_{\rm gl}^{n,k}$ is binary,
the operator $\max(z_{\rm gl}^{n,k},z_{\rm gl}^{n,k-1})$ in (10d) can be
equivalently represented by the linear form $(z_{\rm gl}^{n,k}+z_{\rm
gl}^{n,k-1})$. The formulation (10) is equivalent to the SOS2 constraints in
(4c) and (4d).

##  3 Methodology

###  3.1 Early fixing

Suppose we can determine which linearization region of $q_{\rm l}^{n}$ is to
be selected (i.e., which variables $\eta_{k}^{n}$ and $\eta_{j}^{n}$ can have
nonzero values). In that case, the SOS2 constraints can be removed; thus, the
problem becomes completely linear. This is equivalent to fixing the $z$
variables in the standard formulation (9). Therefore, let us write $P(\pi,z)$
as the problem (9) but with fixed $z$ values, i.e., with integer variables $z$
treated as parameters.

An early fixing heuristic provides an assignment $\hat{z}$ to the integer
variables. Ideally, the assignment will be such that $P(\pi,\hat{z})=P(\pi)$.
Since the early-fixed problem is an LP, it can be solved much faster than the
original MILP problem. Therefore, the total cost of solving the early-fixed
problem is the cost of solving the LP and the cost of running the heuristic.
In practice, however, a trade-off between the cost of running the heuristic
and the gap between $P(\pi,\hat{z})$ and $P(\pi)$ is expected.

Our proposed approach is to develop an early fixing heuristic based on a deep-
learning model. We want a deep learning model $N:\Pi\to[0,1]^{n_{z}}$ with
which we can compute $\hat{z}=\lfloor N(\pi)\rceil$. Two distinct approaches
to training such a model are proposed.

###  3.2 Supervised learning approach

We can train a model for early fixing by feeding it with instances of the MILP
problem and optimizing the model’s parameters such that its output
approximates the optimal binary assignment. Let us define a dataset

$D_{\rm sup}=\left\\{(\pi,z^{\star})\in\Pi\times
Z:P(\pi,z^{\star})=P(\pi)\right\\},$

which associates instances’ parameter vectors $\pi$ with the optimal binary
assignment $z^{\star}$. Note that this dataset requires us to solve to
optimality all MILP instances available.

Let us define a deep learning model

$\begin{split}N_{\rm sup}:\Pi\times\Theta_{\rm sup}&\to[0,1]^{n_{z}}\\\
\pi;\theta_{\rm sup}&\mapsto N_{\rm sup}(\pi;\theta_{\rm sup})\end{split}$
(11)

for which $\theta_{\rm sup}$ is the vector of model parameters that can be
trained. Then, it is possible to optimize the model’s parameters such that for
each vector of parameters in $D_{\rm sup}$, the model’s output approximates
$z^{\star}$. Namely,

$\min_{\theta_{\rm sup}}\sum_{(\pi,z^{\star})\in D_{\rm sup}}\mathcal{L}_{\rm
sup}(N_{\rm sup}(\pi;\theta_{\rm sup}),z^{\star})$ (12)

where $\mathcal{L}_{\rm sup}$ can be, for example, the binary cross entropy
between the elements of both vectors.

###  3.3 Weakly-supervised learning approach

We propose an alternative learning approach that does not require solving MILP
problems. First, we recall that the target of the deep learning model is to
provide $\hat{z}$ such that $P(\pi,\hat{z})$ is maximized111Since (9) is a
maximization problem.. Our proposed approach is to train a _surrogate_ model
that approximates $P(\cdot,\cdot)$, and differentiate through this surrogate
to train the early fixing heuristic using gradient descent methods.

This approach requires only a dataset of assignments for the integer variables
paired with the objective of the respective early-fixed problem. Let us define

$D_{\rm weak}=\left\\{(\pi,\hat{z},p)\in\Pi\times
Z\times\mathbb{R}:p=P(\pi,\hat{z})\right\\},$

which is built with (random) samples of $\hat{z}\in Z$. Then, we train a model

$\begin{split}S:\Pi\times[0,1]^{n_{z}}\times\Theta_{S}&\to\mathbb{R}\\\
\pi,\hat{z};\theta_{S}&\mapsto S(\pi,\hat{z};\theta_{S})\end{split}$ (13)

in a supervised manner, such that its output approximates $P(\pi,\hat{z})$,
that is,

$\min_{\theta_{S}}\sum_{(\pi,\hat{z},p)\in D_{\rm
weak}}\mathcal{L}_{S}(S(\pi,\hat{z};\theta_{S}),p),$ (14)

where $\mathcal{L}_{S}$ can be, for example, the mean squared error.

Now, let us define

$\begin{split}N_{\rm weak}:\Pi\times\Theta_{\rm weak}&\to[0,1]^{n_{z}}\\\
\pi;\theta_{\rm weak}&\mapsto N_{\rm weak}(\pi;\theta_{\rm weak}),\end{split}$
(15)

which can be trained in an unsupervised manner to maximize the surrogate
model’s output given the candidate values for fixing

$\max_{\theta_{\rm weak}}\sum_{\pi\in_{R}\Pi}S(\pi,N_{\rm
weak}(\pi;\theta_{\rm weak});\theta_{S}),$ (16)

where $\pi\in_{R}\Pi$ denotes that $\pi$ is chosen randomly from $\Pi$, and
the summation is iterated over as many random samples as desired. Note that
only $\theta_{\rm weak}$ is optimized during the training of $N_{\rm weak}$,
i.e., $S$ is unchanged.

##  4 Experiments and Results222Code available in github/brunompacheco/early-
fixing-oil-production.

###  4.1 Data

For our experiments, we consider the problem of optimizing oil production from
a single well ($|\mathcal{N}|=1$). We use data from a real subsea oil well,
provided by Petrobras444The data is not made available as it is an
intellectual property of Petrobras. . We set a target for the early fixing
models to generalize to different values of ${\rm bsw}$, ${\rm gor}$, and
$\bar{q}_{\rm gl}$. Therefore, we describe the parameter space of our problem
as

$\pi=({\rm bsw},{\rm gor},\bar{q}_{\rm
gl})\in\Pi=[0.5,1.0]\times[0,300]\times[4000,12500]$

assuming that the BSW can be no lower than 0.5, the GOR is always smaller than
300, and the maximum lift-gas flow is always larger than $10^{5}$. More
precisely, we assume that, in practice, ${\rm bsw}\sim U(0.5,1.0)$, ${\rm
gor}\sim N(100,25)$, and $\bar{q}_{\rm gl}\sim U(4000,12500)$. Note that we
omit the liquid flow function from the parameter vector, once we deal with a
single well and, thus, $\mathcal{Q}_{\rm liq}$ can be uniquely determined by
the other parameters.

The liquid flow function is always linearized with the same breakpoints for
both $q_{\rm gl}$ and ${\rm whp}$. This makes the domain of the binary
variables consistent across instances. The selected breakpoints are

$\displaystyle{\rm whp}^{j}$ $\displaystyle=14k-4,\quad$ $\displaystyle
k=1,\dots,6$ (17a) $\displaystyle q_{\rm gl}^{k}$
$\displaystyle=2500j-2500,\quad$ $\displaystyle j=1,\dots,6.$ (17b)

An example of $\mathcal{Q}_{\rm liq}$ with these breakpoints can be seen in
Figure 1. With the fixed breakpoints, we can describe the domain of the binary
variables as

$\begin{split}Z=\left\\{(z_{\rm gl},z_{{\rm
whp}})\in\\{0,1\\}^{5}\times\\{0,1\\}^{5}:\sum_{i=1}^{5}z_{\rm
gl}=1,\right.\\\ \left.\sum_{i=1}^{5}z_{{\rm whp}}=1\right\\},\end{split}$
(18)

in which $z_{\rm gl}$ indicates the pair of $\eta_{k}$ variables that can take
nonzero values, while $z_{{\rm whp}}$ indicates the pair of $\eta_{j}$
variables that can take nonzero values.

Figure 1: Simulated liquid flow $\mathcal{Q}_{\rm liq}$ (see Eq. (2)) based on
data of a real well using Marlin. The color scale is directly proportional to
the z-axis. We adopted $q_{\rm l}^{n,k,j}=-1$ for combinations of $q_{\rm
gl}^{n,k}$ and ${\rm whp}^{n,j}$ that were considered infeasible by the
simulator.

For the supervised learning task, we build $D_{\rm sup}$ with 500 instances of
the MILP problem from different combinations of ${\rm bsw}$, ${\rm gor}$, and
$\bar{q}_{\rm gl}$. Gurobi was used to solve the MILPs, upon which
$z^{\star}=(z^{\star}_{\rm gl},z^{\star}_{{\rm whp}})$, the optimal value for
the binary variables, was extracted. We use the same $\pi$ vectors as in
$D_{\rm sup}$ for the weakly-supervised learning task. For each $\pi$ we draw
12 random candidates $\hat{z}\in Z$, that is, we ensure that each $\hat{z}$
respects the constraints in (10). Then, $D_{\rm weak}$ is built by solving the
LPs using Gurobi to compute $p=P(\pi,\hat{z})$. Whenever $\pi$ and $\hat{z}$
resulted in an infeasible problem, we added them to $D_{\rm weak}$ with
$p=-1$. In total, $D_{\rm sup}$ contains 500 data points of the form $(({\rm
bsw},{\rm gor},\bar{q}_{\rm gl}),z^{\star})$, while $D_{\rm weak}$ contains
6000 data points of the form $(({\rm bsw},{\rm gor},\bar{q}_{\rm
gl}),\hat{z})$.

###  4.2 Supervised learning experiments

For supervised learning, using $D_{\rm sup}$, we choose

$\begin{split}N_{\rm sup}:\Pi\times\Theta_{\rm
sup}&\to\Delta^{5}\times\Delta^{5}\\\ {\rm bsw},{\rm gor},\bar{q}_{\rm
gl};\theta_{\rm sup}&\mapsto(\hat{z}_{\rm gl},\hat{z}_{{\rm whp}})=N_{\rm
sup}({\rm bsw},{\rm gor},\bar{q}_{\rm gl};\theta_{\rm sup}),\end{split}$ (19)

where $\Delta^{5}$ is the $5$-dimensional simplex set, as a neural network
with 2 hidden layers of $25$ neurons each and a $10$-dimensional output. The
inputs are normalized before the first layer. We use ReLU activation for all
layers except the last one. The last layer’s output is divided into two
vectors of dimension $5$, passing through the softmax function, thus mapping
both into 5-dimensional simplexes. The softmax in the output of the network
ensures that, after rounding, it respects the binary constraints, that is,
$\lfloor\hat{z}_{\rm gl}\rceil$ and $\lfloor\hat{z}_{{\rm whp}}\rceil$ respect
(10).

We randomly split $D_{\rm sup}$ into training and test sets with an 80-20
ratio. The training data was used to select the optimal network architecture
and hyperparameters (learning rate, batch size, and number of epochs). After
this adjustment, the entirety of the training set was used to train $5$
$N_{\rm sup}$ models with random initialization. Adam is used to optimizing
the parameters of the models on the training set such that $\mathcal{L}_{\rm
sup}$ is minimized (see equation (12)). We use

$\mathcal{L}_{\rm sup}(\hat{z}_{\rm gl},\hat{z}_{{\rm whp}},z_{\rm
gl}^{\star},z_{{\rm whp}}^{\star})=\sum_{i=1}^{k_{\rm gl}-1}{\rm
BCE}(\hat{z}_{{\rm gl},i},z_{{\rm gl},i}^{\star})\\\ +\sum_{j=1}^{k_{\rm
whp}-1}{\rm BCE}(\hat{z}_{{\rm whp},j},z_{{\rm whp},j}^{\star}),$

where BCE is the binary cross-entropy function. We use batches of $64$
elements and an initial learning rate of 0.001. Each model is trained for
$100$ epochs. The performance on the test set is summarized in Table 1.

Table 1: Early fixing performance on the test set. We measure the models’ accuracy in predicting the binary variables’ optimal values. _Infeasible_ is the ratio of instances that became infeasible after early fixing. _Objective gap_ is the mean relative decrease of objective value by performing early fixing, for all instances that remained feasible after early fixing. The values reported are the average of 5 runs (models with random initialization), except for _Baseline_.  Model | Accuracy | Infeasible | Objective gap
---|---|---|---
Supervised | 99.78% | 0.11% | 0.01%
Weakly-supervised | 32.31% | 4.19% | 3.64%
Baseline | 0% | 0% | 21.42%

###  4.3 Weakly-supervised learning experiments

We build the surrogate model

$\begin{split}S:\Pi\times[0,1]^{10}\times\Theta_{S}&\to\mathbb{R}\\\ {\rm
bsw},{\rm gor},\bar{q}_{\rm gl},z_{\rm gl},z_{{\rm
whp}};\theta_{S}&\mapsto\hat{p}=S({\rm bsw},{\rm gor},\bar{q}_{\rm
gl},z;\theta_{S})\end{split}$ (20)

as a neural network with $3$ hidden layers of $10$ neurons each. ReLU is used
as an activation function at each layer except the last one, which has no
activation function. Inputs are normalized and a factor of $2000$ scales the
output. Dropout is applied during training at each hidden layer with a
probability of $20$% for each neuron.

The dataset $D_{\rm weak}$ is split randomly into a training and test set
following an 80-20 ratio. The optimal architecture and hyperparameters
(learning rate, batch size, number of epochs) were determined in the same way
as in Sec. 4.2. We use the squared error as the loss function

$\mathcal{L}_{S}(\hat{p},p)=|\hat{p}-p|^{2}.$

We use Adam to minimize the loss function on the training set, with an initial
learning rate of $0.001$ and mini-batches of $1024$ samples. $5$ models with
random initialization are trained for $20$ epochs each. On the test set, the
surrogate models correctly predicted the feasibility (its output was lower
than $0$) an average of $78.23$% of the time. The average MAE on the instances
that were not infeasible was $60.05$. Figure 2 shows an example of the
surrogate model’s performance compared to the real objective values.

Figure 2: Comparison of the surrogate model $S(\pi,z)$ (in blue) and the
target function $P(\pi,z)$ (in red) for all early fixing possibilities
($\forall z\in Z$) for a sample instance of the problem. Values of -1 in _P_
indicate that the problem is infeasible for that given combination of the
variables.

The early fixing model $N_{\rm weak}$ has the same architecture as $N_{\rm
sup}$. $5$ models with random initialization are trained as described in
Section 3.3, but using the parameters $\pi$ from the training set of $D_{\rm
sup}$ described above. Each $N_{\rm weak}$ model is trained with a different
surrogate model $S$. The models are trained using Adam with an initial
learning rate of $0.01$ and a batch of $64$ for $100$ epochs each. The
performance on the $D_{\rm sup}$ test set is reported in Table 1.

###  4.4 Baseline Model

As a reference for the deep learning models’ performance, we compute the
results of always fixing the same values for the binary variables. In our
problem, the safest option in this approach is always to pick the region with
the smallest values possible for $q_{\rm gl}$ and ${\rm whp}$. This always
results in a feasible problem. The performance of this baseline approach on
the test set of $D_{\rm sup}$ can be seen in Table 1.

###  4.5 Early Fixing Impacts

To evaluate the impacts of early fixing in the optimization, we measure the
runtime of solving the original MILP, the runtime of the early-fixed problem
(which is an LP), and the runtime of the early fixing models. We perform these
experiments on all instances of $D_{\rm sup}$, i.e., with the problems defined
by the parameters $\pi$ in the $D_{\rm sup}$ dataset used in the experiments
above.

We found that the original problem can be solved, on average, in 0.90 ms,
while the early-fixed problem is solved in an average of 0.18 ms. Considering
the 0.08 ms the early-fixing models took, on average, during the experiments,
the early-fixing approach takes, on average, 0.26 ms, representing a 71.11%
runtime reduction.

##  5 Conclusions

Our experiments show that deep-learning-based early fixing models successfully
speed up the optimization of the offshore gas-lifted oil production problem,
with a 71.11% runtime reduction. Training in a supervised learning setting,
although with a significantly higher cost for collecting the training data, is
undoubtedly a superior approach concerning the weakly-supervised setting.
Nevertheless, the experiments with the weakly-supervised approach indicate
that it is possible to develop an early fixing heuristic when optimal
solutions to the MILPs are unavailable or too hard to obtain. Still, the
weakly-supervised approach needs further refinement to achieve competitive
results.

Further research is still necessary on the suitability of the deep-learning-
based early fixing for harder problems, e.g., multiple wells connected to
manifolds, with more operational constraints and more integer variables.
Moreover, the early fixing approaches presented in this paper theoretically
apply to any MILP. However, the performance of the models and the practical
viability of them is still an open theme for research.

## References

  * Anderson et al. (2022) Anderson, L., Turner, M., and Koch, T. (2022).  Generative deep learning for decision making in gas networks.  _Mathematical Methods of Operations Research_ , 95(3), 503–532.  10.1007/s00186-022-00777-x.
  * Beale and Tomlin (1969) Beale, E. and Tomlin, J. (1969).  Special facilities in a general mathematical programming system for nonconvex problems using ordered sets of variables.  _Operational Research_ , 69, 447–454.
  * Bengio et al. (2021) Bengio, Y., Lodi, A., and Prouvost, A. (2021).  Machine learning for combinatorial optimization: A methodological tour d’horizon.  _European Journal of Operational Research_ , 290(2), 405–421.  10.1016/j.ejor.2020.07.063.
  * Ding et al. (2019) Ding, J.Y., Zhang, C., Shen, L., Li, S., Wang, B., Xu, Y., and Song, L. (2019).  Accelerating primal solution findings for mixed integer programs based on solution prediction.  10.48550/arXiv.1906.09575.
  * Lee and Mitchell (2009) Lee, E.K. and Mitchell, J.E. (2009).  Integer Programming: Branch and Bound Methods.  In C.A. Floudas and P.M. Pardalos (eds.), _Encyclopedia of Optimization_ , 1634–1643. Springer US, Boston, MA.  10.1007/978-0-387-74759-0_286.
  * Li and Wu (2022) Li, L. and Wu, B. (2022).  Learning to accelerate approximate methods for solving integer programming via early fixing.  10.48550/arXiv.2207.02087.
  * Luguesi et al. (2023) Luguesi, C., Camponogara, E., Seman, L.O., González, J.T., and Leithardt, V.R.Q. (2023).  Derivative-free optimization with proxy models for oil production platforms sharing a subsea gas network.  _IEEE Access_ , 11, 8950–8967.  10.1109/ACCESS.2023.3239421.
  * Müller et al. (2022) Müller, E.R., Camponogara, E., Seman, L.O., Hülse, E.O., Vieira, B.F., Miyatake, L.K., and Teixeira, A.F. (2022).  Short-term steady-state production optimization of offshore oil platforms: wells with dual completion (gas-lift and ESP) and flow assurance.  _TOP_ , 30(1), 152–180.  10.1007/s11750-021-00604-2.
  * Pacheco et al. (2023) Pacheco, B.M., Seman, L.O., Rigo, C.A., Camponogara, E., Bezerra, E.A., and dos Santos Coelho, L. (2023).  Graph neural networks for the offline nanosatellite task scheduling problem.  ArXiv:2303.13773 [cs.LG].
  * Seman et al. (2020) Seman, L.O., Miyatake, L.K., Camponogara, E., Giuliani, C.M., and Vieira, B.F. (2020).  Derivative-free parameter tuning for a well multiphase flow simulator.  _Journal of Petroleum Science and Engineering_ , 192, 107288.  10.1016/j.petrol.2020.107288.
  * Vigerske and Gleixner (2018) Vigerske, S. and Gleixner, A. (2018).  SCIP: global optimization of mixed-integer nonlinear programs in a branch-and-cut framework.  _Optimization Methods and Software_ , 33(3), 563–593.  10.1080/10556788.2017.1335312.