# The RITAS algorithm: a constructive yield monitor data processing algorithm

Luis Damiano  Department of Statistics, Iowa State University, Ames, IA, US
Jarad Niemi  Department of Statistics, Iowa State University, Ames, IA, US

(December 11,2020)

###### Abstract

Yield monitor datasets are known to contain a high percentage of unreliable
records. The current tool set is mostly limited to observation cleaning
procedures based on heuristic or empirically-motivated statistical rules for
extreme value identification and removal. We propose a constructive algorithm
for handling well-documented yield monitor data artifacts without resorting to
data deletion. The four-step Rectangle creation, Intersection assignment and
Tessellation, Apportioning, and Smoothing (RITAS) algorithm models sample
observations as overlapping, unequally-shaped, irregularly-sized, time-
ordered, areal spatial units to better replicate the nature of the destructive
sampling process. Positional data is used to create rectangular areal spatial
units. Time-ordered intersecting area tessellation and harvested mass
apportioning generate regularly -shaped and -sized polygons partitioning the
entire harvested area. Finally, smoothing via a Gaussian process is used to
provide map users with spatial-trend visualization. The intermediate steps as
well as the algorithm output are illustrated in maize and soybean grain yield
maps for five years of yield monitor data collected at a research agricultural
site located in the US Fish and Wildlife Service Neal Smith National Wildlife
Refuge.

††footnotetext: Creative component for the Master of Science in Statistics
https://dr.lib.iastate.edu/handle/20.500.12876/93760

##  1 Introduction

As defined by the National Research Council (1997), precision agriculture is a
data-centered discipline comprising data acquisition at an appropriate scale,
data interpretation and analysis, and management response at an appropriate
scale and time. Miller et al. (1988) were one of the firsts to employ
geostatistics explicitly in precision agriculture, a discipline that started
in the late 80’s when agriculture shifted from farm-level to site-specific
crop management (Oliver, 2010). Their study, among other things, introduced
geostatistics for mapping patterns in soil phosphorus or potassium via
interpolation. Computer mapping of yield and soil is one of the main uses of
this technology, typically to help customize crop management across and within
fields by identifying less productive areas at a sub-field scale (Lowenberg-
DeBoer and Erickson, 2019). Yield data are now recorded automatically for a
wide variety of crops including cereal grains, oilseeds, fiber, forage,
biomass, fruits and vegetables. These data are known to be accurate at a
global scale, yet there exist nuances at a local scale that affect
visualization and downstream analysis, and jeopardizes the credibility and
validity of management decisions downstream. The advent of big data has been
impacting this field largely, especially in terms of data acquisition,
interpretation, and analysis.

Numerous works have established that yield datasets contain a high percentage
of unreliable data and have proposed cleaning procedures (Blackmore and
Marshall, 1996; Moore, 1998; Blackmore, 1999; Thylén et al., 2000; P et al.,
2003; Simbahan et al., 2004; Ping and Dobermann, 2005; Sudduth and Drummond,
2007; Sudduth et al., 2012; Mark et al., 2013; Leroux et al., 2018, 2019; Vega
et al., 2019). Mechanical measurement errors have been studied in great depth
(Arslan and Colvin, 1999, 2002; Grisso et al., 2002; Burks et al., 2004;
Hemming and Chaplin, 2005; Fulton et al., 2009; Schuster et al., 2017), yet we
could not find data adjustment protocols based on mechanistic models.
Typically, systematic errors are identified using either heuristics or
statistical rules with acceptable empirical results but limited probabilistic
motivation. With almost unanimous preference for systematic error detection
and removal, we have found little development of alternative strategies for
error correction or integration such as those proposed by Bachmaier (2007,
2010); in fact, some methodologies discard as much as one third of the
original dataset as summarized by Lyle et al. (2013).

In many cases, procedures such as (Vega et al., 2019) require the user to set
values for tuning parameters such as thresholds. These constants are not
learned from the data but set arbitrarily by the user, potentially hindering
comparisons across users, locations, and years. Moreover, manually-set
thresholds become a hurdle in the context of big data: both the amount and the
heterogeneity of yield monitor datasets have been increasing as raw data
become more abundant and diverse. Additionally, the methods reviewed in (Lyle
et al., 2013) have a general preference for working with crop yield as the
main input of the cleaning procedures. Being a quantity fusing data collected
from more than one sensor, each with its own propagating measurement error, it
typically has a lower signal-to-noise ratio than the mass random variable.
Also, scaling mass to homogenize unequally-sized observational units may
potentiate extreme values resulting in a new variable with higher variability.

Some processing rules, for example the positional error removal in (Blackmore,
1999), are based on the marginal distribution of yield and fail to account for
the spatial nature of the data. Others such as (Leroux et al., 2018; Vega et
al., 2019) collapse all the information in one point on a 2-dimensional plane,
thus failing to recognize that the recorded data are in fact associated with
overlapping, unequally-shaped, and irregularly-sized areal units. Each
observation is in reality a realization from a continuous spatial stochastic
process: time order and spatial superposition are intrinsic characteristics of
the destructive sampling scheme worth modeling.

The many variables logged by the yield monitor equipment include harvested
mass, geocordinates, distance traveled, and swath width, which can be used to
create an areal representation of the data. Depending on the hardware, the
reported travel distance and speed may be measured by a speed sensor, a GPS
receiver, a radar, or a ultranosic sensor (Mulla, 2013). When distance is not
available, it can be linearly estimated using speed and cycle length or less
preferably approximated using the euclidean distance between two subsequent
coordinates. The swath width, also known as gathering or cutting width, can be
time dependent and the combine’s path that may not be correctly reflected in
the data logged by the monitor (Ross et al., 2008).

Blackmore (1999) observed that there is a discrepancy between the theoretical
and the effective harvested area. One of the reasons is that the recorded
header width differs from the width of the header that is truly full of crop
at any given time step. While in practice it can be set equal to a fixed
proportion of the cutter bar width, e.g. 95%, land finishes and areas close to
voids require special treatment as the header could even be empty. In order to
avoid any use of the recorded width, and circumvent the uncertainty associated
with it, the authors introduced potential mapping. In this technique, the
recorded mass of data points belonging to circular neighbourhoods is
aggregated and assigned to a new spatial point whose location corresponds to
the circle geocenter. The aggregated mass is then spatially re-normalized by
the area of the polygons in the Voronoi diagram formed by these new points.
This approach has two limitations. First, it transforms data into spatial
points and thus incurs in information loss about the shape and size of the
areal units (e.g. different shapes could have the same centroid). Second, an
additional gridding technique needs to be applied to analyze temporal trends
in yield maps (Blackmore et al., 2003), effectively displacing the yield
measurements twice: all the polygon information is first collapsed into a
single point which is displaced afterwards.

Previously, Han et al. (1997) introduced a bitmap approach to determine the
actual combine cut width, compute the effective area size, and approximate the
new centroid of each spatial unit after removing the area covered by previous
observations. To track the covered area at each time step, the author
superimposes a regular grid where each cell functions as a bitmask indicating
whether the pixel has been harvested before or not, hence the name of the
methodology. The article does not show how this methodology performs in the
presence of sharp turns, where overlaps tend to be larger during the pivoting
motion, as these seem to have been removed from the analysis. It is worth
noting that, even after processing, the intermediate spatial representations
have some overlap and skips left unexplained by the author. Again, as these
intermediate spatial representations are collapsed into spatial points, some
information is lost. Furthermore, these new data points may not be spatially
aligned if temporal analysis was pursued next. Finally, the authors warn that
the bitmap initialization process is more complicated when there exists non-
crop features such as rivers, canals, or roadways.

Expanding on the bitmap concept, Drummond et al. (1999) developed a method
where observational units are represented by polygons constructed from
position and trajectory information. These are processed in chronological
harvest order by Boolean subtraction to compute the actual harvested area
during each time step, effectively retaining full information about the areal
units. As drawbacks, the authors mention the potential overestimation at the
boundaries as well as the computational complexity of the algorithm. In the
general case, its complexity is order $O(N^{2})$, yet several strategies for
specific cases are discussed there.

The principal contribution of this work is to propose a new constructive
algorithm for yield monitor data processing that accounts for overlapping,
unequally-shaped, irregularly-sized, time-ordered, areal spatial units without
resorting to data deletion. Data collection specifics are firstly described.
Then, the RITAS algorithm is motivated and detailed on a step-by-step basis.
Finally, our methodology is illustrated with a real dataset application.

##  2 Data collection

For our purposes, we define crop yield as crop grain mass per unit area. Our
quantity of interest is a ratio of two quantities that involve a set of
typically more than six sensors, antennas, and devices attached to combine
harvesters that altogether calculate and record grain yield in real time as
the machine operator harvests the productive field. At the end of each cycle,
which lasts a pre-specified number of seconds, the monitor logs measurements
from several on-board sensors for more than 15 variables related to
geolocation, crop characteristics, and machine diagnostics. Dataset size could
range from one to more than a hundred thousand or more observations depending
on factors such as field size, cycle length, and land features.

Grain yield mapping using geo-referenced measurements from a combine-mounted
grain yield monitor are prevalent in agriculture. Reported accuracy of
continuous yield monitoring ranges 93-99.5% depending on equipment type and
brand, calibration regime, flow rate, and environmental conditions at harvest
(Birrell et al., 1996; Fulton et al., 2009; Lyle et al., 2013). Arslan and
Colvin (2002) investigates the specific impact of yield monitor calibration in
yield estimation accuracy.

Previous agronomical studies offer insight on the physicalities of the natural
system whose sampling process we replicate. Ross et al. (2008) is a synoptic
work focused on yield estimation from yield monitor data covering conceptual,
modeling, and mechanical aspects. Grain mass measurement involves a series of
operations to separate clean grain from other materials: gathering, cutting,
pickup, and feeding, threshing, separation, and cleaning. As a result, there
is a time delay in flow rate measurement as well as several mass losses
associated with plant gathering, processing, and leaking until the grain is
measured on or at the exit of the clean grain elevator. At the end of each
cycle, a receiver estimates the GPS antenna position, whose placement in the
combine varies across models. Their capabilities vary largely, typically being
more accurate at measuring velocity rather than position and having a root
mean square horizontal uncertainty of 0.5-3.

##  3 The RITAS algorithm

Our algorithm constructs yield maps through the following steps: Rectangle
creation, Intersection assignment, Tessellation, Apportioning, and Smoothing
(RITAS). The overarching goal of this process is to mimic the real world
harvesting processes. Figure 1 provides an illustration of a subset of these
steps.

Figure 1: Close-up illustration of selected processing steps. Left:
Construction of the vehicle polygons from the GPS location data. Black dots
mark the location at the end of each logging cycle, green dots correspond to
the vertices computed according to the displacement vector implied by two
consecutive spatial points. The distance traveled could be reported by the
yield monitor or estimated from the vector length. Center: Spatial polygons
reveal overlap (darker areas) due to driving maneuvers. Right: Tessellation
eliminates the overlap by assigning intersecting areas to the first polygon in
time.

Each precision yield data set is assumed to have time-ordered rows containing
the following information: mass harvested $m_{t}$, 2-dimensional spatial
coordinate $(x_{t},y_{t})$, and swath half-width $w_{t}$ for $t=0,\ldots,T,\
T\in\mathbb{N}$. We assume any lag time has been effectively pre-processed by
appropriately matching the mass harvested to its 2-dimensional location. See,
for example, Burks et al. (2004); Hemming and Chaplin (2005); Sudduth and
Drummond (2007); Sudduth et al. (2012).

### Step 1: Rectangle creation

Figure 1 (left) illustrates the construction of a rectangular polygon
representing the harvested area between each sequential pair of spatial
coordinates. The position vector $\mathbf{s_{t}}=(x_{t},y_{t})$ represents the
location of the tracked device in a 2-dimensional plane at the time step $t$,
which we assume to be the midpoint of the combine harvester head. The
rectangle is then uniquely identified by the position of its four vertices,
two representing the beginning of the harvested area at time step $t-1$ and
two representing the end of the harvested area at time step $t$. The linear
displacement vector is equal to the vector difference between the position
vectors at two subsequent time steps $\mathbf{s_{t}}-\mathbf{s_{t-1}}$. The
first two vertices are computed as the endpoints of a line segment
perpendicular to the displacement vector with midpoint at $\mathbf{s_{t-1}}$
and length equal to the swath width $2w_{t}$. The remaining two vertices are
found using the same procedure but pivoting on the midpoint $\mathbf{s_{t}}$
instead. Given a pair of position vectors $\mathbf{s_{t}},\mathbf{s_{t-1}}$,
denote $b_{t}=(y_{t}-y_{t-1})/(x_{t}-x_{t-1})$ the slope of the connecting
line, and define $dx_{t}=w_{t}(1+b_{t}^{-2})^{-\frac{1}{2}}$ and
$dy_{t}=-dx_{t}b_{t}^{-1}$. Then, the rectangle associated with mass $m_{t}$
has the four vertices coordinates given by $\\{(x_{i}\pm dx_{t},y_{i}\pm
dy_{t}):i=t,t-1\\}$.

As the first spatial point has no displacement, this processing yields $T$
rectangles with vertices collected in the set $\mathcal{P}=\\{P_{\tau}$:
$\tau\in\\{1,\dots,T\\}\\}$.

### Step 2: Intersection assignment and Tessellation

Figure 1 (middle) shows end result of this construction which produces
rectangles that overlap, an area called the _intersection_ , due to adjacent
harvester paths. Geometries in $\mathcal{P}$ represent the area over which the
combine harvester passed, which may differ from the effectively harvested area
at each time step since the header may not be full of crop when harvesting the
field boundaries (e.g. on outer edges, around voids), or within them. As a
byproduct of the destructive sampling scheme, the discrepancy is linked to the
spatial superposition of the observational units that arise due to local
harvesting dynamics such as turns, wedges, parallel lines, and traveling over
harvested areas. Globally within a field, this phenomenon can be exacerbated
by factors such as field characteristics, or narrow row crops Ross et al.
(2008). Generally, overlapping cannot be assumed symmetrical nor time-
invariant. For example, pivoting motions overlap on the inner side and the
proportion of duplicated area varies at each time step depending on factors
such as the sharpness of the turns, the angle of the wedges, or the obstacles
faced by the operator.

As a general framework to model the time-varying effectively harvested area,
we run a time-ordered apportioning procedure over the rectangles in
$\mathcal{P}$. Let
$\tilde{P}_{\tau}=P_{t}\setminus\left(\bigcup_{i=1}^{t-1}P_{i}\right)$ be the
time-ordered relative complement of the previously harvested area in the
rectangle corresponding to the time step $t$. By doing this, we effectively
map the set of overlapping rectangles $\mathcal{P}$ into a set of non-
overlapping polygons
$\tilde{\mathcal{P}}=\\{\tilde{P}_{\tau}:\tau\in\\{1,\dots,T\\}\\}$. When the
original dataset has no voids (e.g. unplanted, flooded, or more generally
inaccessible sub areas), $\tilde{\mathcal{P}}$ form a flat plane covered by
tiles with no overlaps and no gaps (non-periodic tessellation). Some of the
computational considerations discussed in Drummond et al. (1999), such as the
time scaling and appropriate techniques to reduce algorithm complexity, are
relevant for implementing this step.

This step creates $T$ polygons partitioning the harvested area, each
associated with the same mass it had in Step 1 but now the area may be smaller
(due to the intersection removal) and thus yield is more accurately captured.

### Step 3: Apportioning

The elements in $\tilde{\mathcal{P}}$ are unsuitable for spatial techniques
that do not accommodate areas with irregular shapes and heterogeneous area
sizes. Two equally-sized long rectangles, one in a vertical and the other in a
horizontal position, would have the same centroid yet they convey different
information about yield to the north or west of the centroids. Alternatively,
two centered equally-shaped polygons with different sizes carry different
information about the spatial coverage of the collected data.

To normalize the areal representation in terms of both shape and size, we
superimpose a regular grid of square pixels, assign portions of the polygons
into grid pixels, and apportion the harvested mass associated with the
polygons to the pixels. The first two steps involve topology operations on
geometries whereas the last one involves manipulation of the numerical data.
Accounting for these two spatial features is by itself a methodological
improvement over the surveyed algorithms, which simply reduce data to spatial
points such as the displacement vector endpoint or, less commonly, the polygon
centroid.

Let $\mathcal{P}^{*}$ be a set of $N\in\mathbb{N}$ equally-sized, non-
overlapping, and contiguous squared pixels forming a partition covering all
the elements in $\tilde{\mathcal{P}}$. The constant $N$ reflects the user
preference in terms of resolution, selected either by the total number of
pixels, more intuitively by the pixel length in meters, or the target
computational time investment. Although pixels size can be set arbitrarily,
its choice should consider the accuracy of the position system. We take the
pairwise intersection among the elements in $\tilde{\mathcal{P}}$ and
$\mathcal{P}^{*}$ and compute $\pi_{\tau,n}\in[0,1]$ the proportion of the
area in the $\tau$-th non-overlapping polygon $\tilde{P}_{\tau}$ that
intersects with the $n$-th pixel in the grid for $\tau\in\\{1,\dots,T\\}$ and
$n\in\\{1,\dots,N\\}$. The corresponding proportion of the harvested mass
$m_{\tau}$ associated with $\tilde{P}_{\tau}$ is assigned to $P^{*}_{n}$. The
total harvested mass associated with the $n$-th pixel is given by
$m^{*}_{n}=\sum_{\tau=1}^{T}\pi_{\tau,n}\ m_{\tau}$. The resulting polygons in
$\mathcal{P}^{*}$ resemble much the Basic Areal Units (BAUs) as defined by
Nguyen et al. (2012) in the context of massive data fusion: fine-scale,
nonoverlapping, areal regions representing the smallest resolution at which
data is aggregated.

Two key aspects behind the gridding strategy are worth mentioning. First, for
the purpose of apportioning, we assume that the harvested mass associated with
the spatial polygons is distributed uniformly within each unit. This is a
sensible assumption in this context as combine harvester log data in short
cycles and the spatial polygons represent small areas for the scale of the
underlying crop growth process.

Second, if one imposes regularity and equal-shape conditions on the grid
elements and also forces it to cover all the elements in
$\tilde{\mathcal{P}}$, the sum of the pixels area may exceed the sum of the
tessellated polygons area. The excess, found at both the harvested region
outer borders and the inner voids boundaries, can be diagnosed visually with
ease. It decreases as the grid resolution, or equivalently as the total number
of pixels $N$, increases and so can be controlled at the cost of additional
computational time for topology operations and smoothing. Concretely, when a
Gaussian Process is applied for smoothing as described below, exact inference
has time complexity in the order of $O(N^{3})$ and storage demands of of
$O(N^{2})$. In other words, as we double the map resolution, we octuple the
number of operations and quadruple the need for storage. Alternatively, one
could exclude the pixels with less than an arbitrary proportion of area
effectively covered by the tessellated observations. A sensible choice, such
as a minimum coverage of 50%, will tend to balance under and over-covered
polygons. Due care should be taken so that apportioning is still applied
validly: mass should be allocated to pixels only in proportion of the actual
overlapping area, discarding any part of the tessellated polygons that is not
covered by any pixel. Note that the grid resolution serves only for the
purpose of spatial aggregation, and need not be the same resolution used for
the visualization that will be introduced in the following step.

This step creates $N$ polygons, determined by the user based on tesselation
resolution, partitioning the harvested area each with an associated mass of
harvested crop. Figure 2 (bottom left) shows the regular grid with apportioned
mass. The regular tesselation provides constant areas and meaningful
centroids, and therefore we can use standard spatial smoothing techniques
directly on mass, as opposed to yield.

### Step 4: Smoothing

McCullagh and Clifford (2006) provided empirical evidence suggesting that the
non-anthropogenic spatial variation in yield, defined as patterns that cannot
be explained by topography or human intervention, matches the characteristics
of the de Wijs process plus white noise. We smooth using a Gaussian Process
(GP) with a Matérn covariance on the logarithm of mass (Handcock and Stein,
1993; Guttorp and Gneiting, 2006), which becomes similar to a de Wijs process
as the length scale tends to zero. Compared with the more common powered
exponential covariance functions, e.g. Gaussian or exponential, the Matérn
adds an additional parameter that controls local smoothness, i.e.
differentiability, and therefore is often more accurate for real world
processes.

Covariance parameters are estimated via Maximum Likelihood, and smoothed
values are found for each tile following Cressie (1993). Specifically, for
each pixel we have a predicted mean $\hat{\mu}_{\ell}\in\mathbb{R}$ and
variance $\hat{\sigma}^{2}_{\ell}>0$ for the logarithm of mass. We using the
following formulas to convert back to the mean and variance of mass

$\hat{\mu}_{m}=\exp\left(\hat{\mu}_{\ell}+\hat{\sigma}^{2}_{\ell}/2\right),\quad\mbox{and}\quad\hat{\sigma}^{2}_{m}=\exp\left(2\hat{\mu}_{\ell}+\hat{\sigma}^{2}_{\ell}\right)\left[\exp\left(\hat{\sigma}^{2}_{\ell}\right)-1\right].$

Finally, yield is calculated by dividing the tile mass by the tile area.

We have implemented our algorithm in the R programming language, using the R
packages doParallel, foreach, gstat, rgeos, and sp (Pebesma, 2004; Pebesma and
Bivand, 2005; Bivand et al., 2013; Gräler et al., 2016; Microsoft and Weston,
2017; Corporation and Weston, 2018; R Core Team, 2019; Bivand and Rundel,
2019).

##  4 STRIPS yield mapping

In this section, we illustrate the functioning and the results of our
methodology when applied to yield monitor data collected from the same
agricultural site over the years. We start with a detailed discussion of the
production of the yield map for one specific year, and we finally show some
resulting visualizations for the same fields but different crops and years.

The data arises from a study conducted at the Neal Smith National Wildlife
Refuge in central Iowa to quantify the impact of grassland-to-cropland
conversion on nitrate-nitrogen (NO3–N) concentrations in soil and shallow
groundwater and to assess the potential for perennial filter strips to
mitigate increases in NO3–N levels, run-off reduction, and soil nutrient loss
(Zhou et al., 2010). The experiment was run in different study sites within
the 3000-ha area managed by the U.S. National Fish and Wildlife Service,
located in the Walnut Creek watershed in Jasper County, Iowa. In this case
study, we focus on one specific site named Basswood that is situated west to
the Basswood Trailhead.

Basswood, located at WGS84 15 N 0477097E 4598644N, has a total area of
approximately 13 Ha. Nearly 81% of the surface is cropland; most of the
remaining proportion is reconstructed prairie vegetation planted as part of
the experimental design. For the purpose of our yield maps, these areas are
treated as voids because no data were collected from that surface. Cropland in
the experiment is in a maize–soybean rotation using standard no-till soil and
weed-management techniques. Geographic coordinates, sample time, moisture
content, and maize (2008, 2010, 2012, 2014) and soybean (2009, 2011, 2013,
2015) flow rate were reported by a Case IH AFS Pro-600 combine-mounted yield
monitor every 1-3 s during crop harvest, resulting in a fine-scale spatially
referenced dataset of crop yields across the study area. Schulte et al. (2017)
provides a detailed account of the experiment protocol and the resulting
improvements in the biodiversity and the delivery of multiple ecosystem
services.

The yield monitor dataset for the year 2012 has 4,239 observations logged
every three seconds starting from the southwest corner of the site. The swath
width was reported to be constant at 6.10m. The distance traveled during each
cycle, with an overall mean of 3.7 m, has three modes with centers at 1.9,
4.0, and 6.0 m. The distribution of the yield reported by the monitor, with
its median located at 6.3 mg/ha and the 90% of observations being within 1.4
and 10.8 mg/ha, is symmetric and platykurtic. Visual inspection suggests that
there are approximately 20 extreme values on the right tail. Small areas
within the field boundaries without data, visualized as voids in the maps,
correspond to small portions of soil allocated for nonproductive purposes
(e.g. perennial crops, research equipment).

Figure 2: Visualization of selected intermediate steps involved in the
algorithm for the harvested grain yield of Basswood on year 2012 (maize). Top
left: Point map where each observation is marked with equally-sized points
placed at the each logging location. Information relevant to spatial trends
are hard to observe, such as area coverage and overlaps. Top right: Map of the
clipped polygons. Highly intertiment coloring in noisy areas hinder the
visualization of the spatial trends, hence the need for smoothing. Bottom
left: Map of the aggregated grid at a 5 m resolution. This step produces
equally-sized, regularly-shaped polygons suitable for spatial interpolation.
Bottom right: Map to be shown to the user. By increasing color homogeneity at
a local level, low- and high-yield areas have a larger contrast and become
easier to identify.

Figure 2 (top left) displays the simplest form of a yield map. Data points are
visualized as equally-sized symbols with colors indicating the yield value at
each recorded location. A single-hue shade of green is chosen to ease the
interpretation with dark shades being intuitively associated with higher crop
yield. The characteristics of the yield distribution vary largely across
sites, crop types, and years. Since the main objective of yield map analysis
is to identify sub-field areas with different performance levels, the points
are colored according to a measure of relative yield within the dataset as
opposed to an absolute scale. As a general rule, we define the color gradient
as a function of the the empirical quartiles in order to guarantee that low
and high yield measurements are uniformly represented at the same time that
the interpretation guidelines remain consistent across maps.

The proposed coloring scheme helps capturing the large-scale spatial trends:
the east area, both above and below the perennial crop strip, and the
northwest area show the best relative performance whereas the southwest area
is underperforming. Another evident pattern is outer borders and borders
around inner voids displaying lower yield, which could be explained by soil
fertility properties or could simply be a byproduct of how data is collected.
Data artifacts due to narrow finishes are well documented: if the header cut
was not full when harvesting the boundaries and the operator failed to
manually flag it, yield would be underestimated. Careful consideration should
be given to neighboring points that are not on the borders. Light green points
on the northern borderline are adjacent to dark green points, suggesting that
this might very well be an artifact due to deficiencies in the the data
collection procedure. On the other hand, data points on the southwest zones
are consistently underperforming, thus suggesting the existence of an actual
spatial trend.

Using points to visualize the data provides no information about the shape of
the areal observational unit and hides the overlaps in the harvested areas,
which should be considered appropriately when computing the estimated yield at
a given spatial location. In fact, from the figure it is not evident that
there is 9.0% of areal overlap, defined as the percent excess of the sum of
the individual rectangles area over the polygons union area.

We apply the RITAS algorithm. Figure 4 (top right), displaying the constructed
spatial polygons, makes overlapping patterns more evident: (i) subsequent
samples overlap during turns, especially on the inner side; (ii) near voids,
where the landscape requires more maneuvering; (iii) wedges formed by
perpendicular passing, for example on the western part of the field; (iv)
driving from one part to another; (v) between parallel passings and narrow
segments.

Because overlapping produces a systematic overestimation of the effective area
size, thus biasing down the estimated yield, its correct treatment can uncover
underperforming areas. In all these cases, the coloring suggests that highly
overlapping polygons are associated with lower yields. We note, however, that
yield was computed using the theoretical polygon area which overestimates the
effective harvest area. In Figure 4 (middle left), which shows the reshaped
polygons and the yield computed with the new effective area, we notice that
some of the sub field areas with low-yield polygons now display a better
performance suggesting that this visual artifact was indeed caused by the
overlapping. As a computational note, when producing the reshaped polygons, 10
spatial polygons whose area had been fully harvested in previous time steps
were dropped from the dataset causing an minor leakage of 0.1% of the total
harvest mass; in case of major situations, aggregating the mass of fully
nested geometries would be appropriate.

As the ultimate goal of the visualization work is to support the user’s
decision making process, the clipped map is not adequate. Crop management
decisions at a sub-field scale are based on spatial trends whereas the
measurements are highly noisy due to a combination of at least 10 possible
types of data collection error as discussed in Lyle et al. (2013). For
example, in Figure 4 there are zones of predominating high and low yield
contaminated with scattered observations with the opposite color, and there
are also zones with a combination where it is difficult to identify local
trends.

Smoothing is thus typically applied. As discussed in the previous section,
off-the-shelf smoothing techniques are not suitable for unequally-sized
polygons. We superimpose a grid with 4,194 squares at a 5-meter resolution and
apportion the harvested mass of the reshaped polygons into the corresponding
pixels. As we retain those pixels with at least 50% of their area covered by
observations only, we find local zones with under or over coverage. Overall,
the whole grid covers 2.6% less of the sampled surface, a discrepancy that can
be diminished by increasing the grid resolution.

The results are seen in Figure 2 (bottom right). The effect of smoothing on
the signal to noise ratio is evident from the figure. The richness of
smoothing is not only on the visualization, but also on the interpretation of
the estimated parameters. We then propose that yield maps should not only
include the spatial polygons, but also statistical information useful for the
map user. The range, for instance, is informative for map users to better
understand the scale of the spatial effect and adjust the scale of their
decisions accordingly.

The smooth map is consistent with the main spatial trends observed so far.
Clear patterns become more evident now; in the mixed areas, smoothing helps
not only to identify the overall local trend but also to make internal
breaks/borders more visible. Additional features can help with interpretation
include contour lines, e.g. each 1 mg/ha similar to Blackmore (1999), or color
schemes based on spatial clusters.

Figure 3: Visualization of the harvested grain yield for Basswood. Rows: 2010
(maize), 2014 (maize), 2009 (soybean), 2011 (soybean). Columns: apportioned
(left), smoothed (right). Despite soybean larger variability, the smooth map
offer a clear visualization of the spatial trends.

##  5 Discussion

We have proposed an algorithm for data processing and visualization that is
designed for situations where the destructive sampling scheme incorporates
other non-trivial features, namely time order and unequally-sized,
irregularly-shaped, and possibly overlapping areal observational units.
Reproducing the sampling process, the rules unpack the georeferenced
coordinates, recorded as points in a 2-dimensional space, into polygons that
are next tiled. The resulting geometries are aggregated using basic areal
units and apportioned assuming that the variable of interest follows a uniform
distribution at the local scale. The final output of the processing rules is
in the form of a regular tessellation of the harvested area, equipped with
both their associated observed and smoothed yield; the output is meant for
both data analysis and visualization. As a byproduct, because the location of
the pixels in both the aggregation and smoothing grids can be chosen freely
without affecting the validity of the algorithm, spatial registration of the
output is readily available.

To illustrate the main characteristics of our procedure, we have presented a
case study of grain yield data processing and visualization. Whereas the
previous work in precision agriculture presents an almost unanimous preference
for extreme values detection and removal, with up to 32% of the data excluded
Lyle et al. (2013), our physically principled approach generates a smooth
visualization while retaining as much information as collected. This
improvement is partly explained by our reproducing of the sampling scheme, but
also by a twist on how the problem has been traditionally framed: recognizing
that local yield is a highly volatile measurement unit due to scaling and
error propagation, and thus effectively shifting the baseline variable from
yield to mass. Also in contrast with previously proposed methodologies, most
of which require the map producer to set arbitrary thresholds, our algorithm
is autonomous. There are two main advantages associated with this: (i) large
amount of grain yield datasets can be processed automatically, a relevant
feature in the current times where data collection has became accessible and
widespread even at a small farming scale; and (ii) the final outputs become
more consistent across different users, fields, and years.

When doing exact inference for a Gaussian Process, smoothing requires the
inversion of a square matrix with size equal to the number of observations in
the aggregation step, which in this case study is 3,738 pixels. Unless
resorting to approximate interpolation methods, significant time is expected
to be consumed for this purpose. As a palliative, our implementation
optionally divides the prediction space into smaller subsets to be computed in
parallel. Although there is no direct gain in the matrix inversion time, it
provides with some marginal improvements as out prediction space is large with
more than 90,000 pixels.

Future improvements of this algorithm could come from three perspectives:
assumptions, visualization, and computations. One of the main assumptions is
that, at a local scale, the random variable of interest follows a uniform
distribution. Although this is a reasonable idealization for the case study
given the small area of the observational units relative to the rate of change
of the underlying process, applications for highly-distanced observations or
rapidly-changing underlying processes may require the extension to a different
distribution (e.g. exponential decay from the centroids of the polygons to its
boundaries). Yield monitors may record additional variables that are currently
not exploited by our algorithm, which could be enhanced for example by
smoothing via universal kriging with automatic relevance determination for
feature selection.

To improve visualizations, the smooth map could profit from better display
techniques to signal the contrast between low and high valued areas (e.g.
contour lines, spatial clustering techniques). In the specific case of grain
yield maps, we still observe that some of the well-known sources of yield data
error transpire through the algorithm. Some of these can be treated previously
to running the algorithm, for example the time lag effect, the harvester fill
mode error, and position errors as described in Blackmore (1999).

On the computational side, the case study suggests that the most time
consuming steps are smoothing and tiling, respectively. The former could be
reduced by resourcing to approximation methods for Gaussian Process spatial
interpolation (Shi and Cressie, 2007; Cressie and Johannesson, 2008; Katzfuss
and Cressie, 2011; Nguyen et al., 2012, 2014), approximate linear algebra
routines, or a more efficient interpolation technique. Although the intrisic
sequential nature of the sampling scheme limits the potential of
parallelization in the tiling step, some performance improvements could be
potentially achieved by subviding the spatial domain into disjoint blocks that
should after be accordingly recoupled (divide and conquer). Alternatively, the
bitmap matrix of Han et al. (1997) could be revisited as an approximation to
our polygon approach. These gains could turn the current work into a near
real-time algorithm.

## References

  * Arslan and Colvin (1999) Arslan, S. and Colvin, T. S. (1999).  Laboratory performance of yield monitor.  Applied Engineering in Agriculture, 15(3):189–195.  Publisher: American Society of Agricultural and Biological Engineers (ASABE).
  * Arslan and Colvin (2002) Arslan, S. and Colvin, T. S. (2002).  Grain yield mapping: yield sensing, yield reconstruction, and errors.  Precision Agriculture, 3(2):135–154.  Publisher: Springer Nature.
  * Bachmaier (2007) Bachmaier, M. (2007).  Using a robust variogram to find an adequate butterfly neighborhood size for one-step yield mapping using robust fitting paraboloid cones.  Precision Agriculture, 8(1):75–93.  Publisher: Springer Science and Business Media LLC.
  * Bachmaier (2010) Bachmaier, M. (2010).  A yield mapping procedure based on robust fitting paraboloid cones on moving elliptical neighborhoods and the determination of their size using a robust variogram.  Positioning, 01(1):27–41.  Publisher: Scientific Research Publishing, Inc,.
  * Birrell et al. (1996) Birrell, S. J., Sudduth, K. A., and Borgelt, S. C. (1996).  Comparison of sensors and techniques for crop yield mapping.  Computers and Electronics in Agriculture, 14(2):215–233.
  * Bivand and Rundel (2019) Bivand, R. and Rundel, C. (2019).  rgeos: interface to geometry engine - open source (’GEOS’).
  * Bivand et al. (2013) Bivand, R. S., Pebesma, E., and Gomez-Rubio, V. (2013).  Applied spatial data analysis with R, second edition.  Springer, NY.
  * Blackmore and Marshall (1996) Blackmore, B. S. and Marshall, C. J. (1996).  Yield mapping: errors and algorithms.  American Society of Agronomy, Crop Science Society of America, Soil Science Society of America.
  * Blackmore (1999) Blackmore, S. (1999).  Remedial correction of yield map data.  Precision Agriculture, 1(1):53–66.  Publisher: Springer Nature.
  * Blackmore et al. (2003) Blackmore, S., Godwin, R. J., and Fountas, S. (2003).  The analysis of spatial and temporal trends in yield map data over six years.  Biosystems Engineering, 84(4):455–466.  Publisher: Elsevier BV.
  * Burks et al. (2004) Burks, T. F., Shearer, S. A., Fulton, J. P., and Sobolik, C. J. (2004).  Effects of time-varying inflow rates on combine yield monitor accuracy.  Applied Engineering in Agriculture, 20(3):269–275.  Publisher: American Society of Agricultural and Biological Engineers (ASABE).
  * Corporation and Weston (2018) Corporation, M. and Weston, S. (2018).  doParallel: foreach parallel adaptor for the ’parallel’ package.
  * Cressie and Johannesson (2008) Cressie, N. and Johannesson, G. (2008).  Fixed rank kriging for very large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(1):209–226.  Publisher: Wiley.
  * Cressie (1993) Cressie, N. A. C. (1993).  Statistics for spatial data.  John Wiley & Sons, Inc.
  * Drummond et al. (1999) Drummond, S. T., Fraisse, C. W., and Sudduth, K. A. (1999).  Combine harvest area determination by vector processing of GPS position data.  Transactions of the ASAE, 42(5):1221–1228.  Publisher: American Society of Agricultural and Biological Engineers (ASABE).
  * Fulton et al. (2009) Fulton, J. P., Sobolik, C. J., Shearer, S. A., Higgins, S. F., and Burks, T. F. (2009).  Grain yield monitor flow sensor accuracy for simulated varying field slopes.  Applied Engineering in Agriculture, 25(1):15–21.  Publisher: American Society of Agricultural and Biological Engineers (ASABE).
  * Grisso et al. (2002) Grisso, R. D., Jasa, P. J., Schroeder, M. A., and Wilcox, J. C. (2002).  Yield monitor accuracy: successful farming magazine case study.  Applied Engineering in Agriculture, 18(2).  Publisher: American Society of Agricultural and Biological Engineers (ASABE).
  * Gräler et al. (2016) Gräler, B., Pebesma, E., and Heuvelink, G. (2016).  Spatio-temporal interpolation using gstat.  The R Journal, 8(1):204–218.
  * Guttorp and Gneiting (2006) Guttorp, P. and Gneiting, T. (2006).  Studies in the history of probability and statistics XLIX on the matérn correlation family.  Biometrika, 93(4):989–995.
  * Han et al. (1997) Han, S., Schneider, S. M., Rawlins, S. L., and Evans, R. G. (1997).  A bitmap method for determining effective combine cut width in yield mapping.  Transactions of the ASAE, 40(2):485–490.  Publisher: American Society of Agricultural and Biological Engineers (ASABE).
  * Handcock and Stein (1993) Handcock, M. S. and Stein, M. L. (1993).  A bayesian analysis of kriging.  Technometrics, 35(4):403–410.  Publisher: Taylor & Francis Group.
  * Hemming and Chaplin (2005) Hemming, N. and Chaplin, J. (2005).  Determining lag time for mass flow in a combine harvester.  Transactions of the ASAE, 48(2):823–829.  Publisher: American Society of Agricultural and Biological Engineers (ASABE).
  * Katzfuss and Cressie (2011) Katzfuss, M. and Cressie, N. (2011).  Spatio-temporal smoothing and EM estimation for massive remote-sensing data sets.  Journal of Time Series Analysis, 32(4):430–446.  Publisher: Wiley.
  * Leroux et al. (2018) Leroux, C., Jones, H., Clenet, A., Dreux, B., Becu, M., and Tisseyre, B. (2018).  A general method to filter out defective spatial observations from yield mapping datasets.  Precision Agriculture, 19(5):789–808.  Publisher: Springer Nature.
  * Leroux et al. (2019) Leroux, C., Jones, H., Clenet, A., and Tisseyre, B. (2019).  Knowledge discovery and unsupervised detection of within-field yield defective observations.  Computers and Electronics in Agriculture, 156:645–659.  Publisher: Elsevier BV.
  * Lowenberg-DeBoer and Erickson (2019) Lowenberg-DeBoer, J. and Erickson, B. (2019).  Setting the record straight on precision agriculture adoption.  Agronomy Journal, 111(4):1552.  Publisher: American Society of Agronomy.
  * Lyle et al. (2013) Lyle, G., Bryan, B. A., and Ostendorf, B. (2013).  Post-processing methods to eliminate erroneous grain yield measurements: review and directions for future development.  Precision Agriculture, 15(4):377–402.  Publisher: Springer Science and Business Media LLC.
  * Mark et al. (2013) Mark, S., Adriano, A., and Jose, M. (2013).  A simple method for filtering spatial data.  In Precision agriculture13 (pp. 259-266), Precision Agriculture 2013 - Papers Presented at the 9th European Conference on Precision Agriculture, ECPA 2013.
  * McCullagh and Clifford (2006) McCullagh, P. and Clifford, D. (2006).  Evidence for conformal invariance of crop yields.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 462(2071):2119–2143.  Publisher: The Royal Society.
  * Microsoft and Weston (2017) Microsoft and Weston, S. (2017).  foreach: provides foreach looping construct for r.
  * Miller et al. (1988) Miller, M. P., Singer, M. J., and Nielsen, D. R. (1988).  Spatial variability of wheat yield and soil properties on complex hills.  Soil Science Society of America Journal, 52(4):1133.  Publisher: Soil Science Society of America.
  * Moore (1998) Moore, M. (1998).  An investigation into the accuracy of yield maps and their subsequent use in crop management.
  * Mulla (2013) Mulla, D. J. (2013).  Twenty five years of remote sensing in precision agriculture: key advances and remaining knowledge gaps.  Biosystems Engineering, 114(4):358–371.  Publisher: Elsevier BV.
  * National Research Council (1997) National Research Council (1997).  Precision agriculture in the 21st century.  National Academies Press.
  * Nguyen et al. (2012) Nguyen, H., Cressie, N., and Braverman, A. (2012).  Spatial statistical data fusion for remote sensing applications.  Journal of the American Statistical Association, 107(499):1004–1018.  Publisher: Informa UK Limited.
  * Nguyen et al. (2014) Nguyen, H., Katzfuss, M., Cressie, N., and Braverman, A. (2014).  Spatio-temporal data fusion for very large remote sensing datasets.  Technometrics, 56(2):174–185.  Publisher: Informa UK Limited.
  * Oliver (2010) Oliver, M., editor (2010).  Geostatistical applications for precision agriculture.  Springer Netherlands.
  * P et al. (2003) P, N., T, M., and M, D. (2003).  An algorithm for automatic detection and elimination of defective yield data. Wageningen Academic Publishers.
  * Pebesma (2004) Pebesma, E. J. (2004).  Multivariable geostatistics in s: the gstat package.  Computers & Geosciences, 30:683–691.
  * Pebesma and Bivand (2005) Pebesma, E. J. and Bivand, R. S. (2005).  Classes and methods for spatial data in r.  R News, 5(2):9–13.
  * Ping and Dobermann (2005) Ping, J. L. and Dobermann, A. (2005).  Processing of yield map data.  Precision Agriculture, 6(2):193–212.  Publisher: Springer Science and Business Media LLC.
  * R Core Team (2019) R Core Team (2019).  R: a language and environment for statistical computing.
  * Ross et al. (2008) Ross, K. W., Morris, D. K., and Johannsen, C. J. (2008).  A review of intra-field yield estimation from yield monitor data.  Applied Engineering in Agriculture, 24(3):309–317.  Publisher: American Society of Agricultural and Biological Engineers (ASABE).
  * Schulte et al. (2017) Schulte, L. A., Niemi, J., Helmers, M. J., Liebman, M., Arbuckle, J. G., James, D. E., Kolka, R. K., O’Neal, M. E., Tomer, M. D., Tyndall, J. C., Asbjornsen, H., Drobney, P., Neal, J., Ryswyk, G. V., and Witte, C. (2017).  Prairie strips improve biodiversity and the delivery of multiple ecosystem services from corn–soybean croplands.  Proceedings of the National Academy of Sciences, 114(42):11247–11252.  Publisher: Proceedings of the National Academy of Sciences.
  * Schuster et al. (2017) Schuster, J. N., Darr, M. J., and McNaull, R. P. (2017).  Performance benchmark of yield monitors for mechanical and environmental influences.  In 2017 ASABE annual international meeting, page 1.  American Society of Agricultural and Biological Engineers.
  * Shi and Cressie (2007) Shi, T. and Cressie, N. (2007).  Global statistical analysis of MISR aerosol data: a massive data product from NASAs terra satellite.  Environmetrics, 18(7):665–680.  Publisher: Wiley.
  * Simbahan et al. (2004) Simbahan, G. C., Dobermann, A., and Ping, J. L. (2004).  Screening yield monitor data improves grain yield maps.  Agronomy Journal, 96(4):1091.  Publisher: American Society of Agronomy.
  * Sudduth and Drummond (2007) Sudduth, K. A. and Drummond, S. T. (2007).  Yield editor.  Agronomy Journal, 99(6):1471.  Publisher: American Society of Agronomy.
  * Sudduth et al. (2012) Sudduth, K. A., Drummond, S. T., and Myers, D. B. (2012).  Yield editor 2.0: software for automated removal of yield map errors.  In 2012 dallas, texas, july 29 - august 1, 2012. American Society of Agricultural and Biological Engineers.
  * Thylén et al. (2000) Thylén, L., Algerbo, P. A., and Giebel, A. (2000).  An expert filter removing erroneous yield data.  In Robert, P. C., Rust, R. H., and Larson, W. E., editors, Proceedings of the 5<span class=”nocase”>${}^{\textrm{th}}$</span> international conference on precision agriculture, bloomington, minnesota, USA, 16-19 july, 2000., pages 1–9. American Society of Agronomy.
  * Vega et al. (2019) Vega, A., Córdoba, M., Castro-Franco, M., and Balzarini, M. (2019).  Protocol for automating error removal from yield maps.  Precision Agriculture.  Publisher: Springer Nature.
  * Zhou et al. (2010) Zhou, X., Helmers, M. J., Asbjornsen, H., Kolka, R., and Tomer, M. D. (2010).  Perennial filter strips reduce nitrate levels in soil and shallow groundwater after grassland-to-cropland conversion.  Journal of Environment Quality, 39(6):2006.  Publisher: American Society of Agronomy.

##  6 Acknowledgements

Funding was provided by Iowa State University through the Presidential
Interdisciplinary Research Initiative on C-CHANGE: Science for a Changing
Agriculture and the Foundation for Food and Agriculture Research (FFAR). The
dataset used to illustrate our work was collected as part of the Science-Based
Trials of Rowcrops Integrated with Prairie Strips (STRIPS) project, which
involved the work of the US Fish and Wildlife Service Neal Smith National
Wildlife Refuge, numerous technicians and student researchers, and Heartland
Cooperative.

(a) Observations as point on a 2D-space

(b) Rectangle creation step output

(c) Tessellation step output

(d) Tessellated plane partition

(e) Apportioning step output

(f) Smoothing step output

Figure 4: Step-by-step progression of Basswood (2012).


# LQ Control of Traffic Flow Models via Variable Speed Limits

Brian Block and Stephanie Stockar This material is based upon work supported
by the National Science Foundation Graduate Research Fellowship Program under
Grant No. DGE-1343012 and NSF CAREER Award 2042354.The authors are with the
Department of Mechanical and Aerospace Engineering and the Center for
Automotive Research, The Ohio State University, 930 Kinnear Road, Columbus, OH
43212 USA<EMAIL_ADDRESS>
###### Abstract

In this paper, an extension of a linear control design for hyperbolic linear
partial differential equations is presented for a first-order traffic flow
model. Starting from the Lighthill-Whitham-Richards (LWR) model, variable
speed limit control (VSL) is applied through a modification of Greenshield’s
equilibrium flow model. Then, an optimal linear quadratic (LQ) controller is
designed on the linear LWR model. The LQ state feedback function is found via
the solution of a Riccati differential equation. Unlike previous studies, the
control input is the rate of change of the input, not the input itself. The
proposed controller is then verified on both the linear and nonlinear models.
In both cases, the controller is able to drive the system to a desired density
profile. In the nonlinear application, a higher control gain is needed to
achieve similar results as in the linear case.

##  I INTRODUCTION

Congestion is a growing concern in major cities as it leads to increased
energy usage and pollution [1]. In 2022, the average commuter in the Unites
States spent around 73 extra hours in traffic. This is up over 40% from 2016
[2, 3]. The modeling of traffic flow and its control offer opportunities to
mitigate congestion, improve travel time, and reduce overall energy usage
associated with the transportation sector.

One approach to modeling traffic flow in cities and other large scale systems
like highways, is through a macroscopic model that looks at the density and
flow of traffic [4]. A common macroscopic traffic flow model is the Lighthill-
Whitham-Richards (LWR) model, which consists of a hyperbolic conservation
partial differential equation (PDE) in density [5, 6]. As a first-order model,
the LWR model relies on an equilibrium flow equation to determine the velocity
of vehicles [4]. The equilibrium flow-density relationship can be modeled
various ways including using a triangular relationship, Greenshield’s model,
or an exponential model [7, 8, 9]. The LWR model is commonly used for large
scale traffic congestion control as it is both simple and accurate [10, 11,
12, 7].

Macroscopic traffic congestion mitigation strategies rely on either boundary
control, such as ramp metering, or in-domain control, such as using variable
speed limits (VSL). For example, ramp metering is used in [13, 14] where
backstepping control is used to drive a continuous traffic flow model to a
density setpoint. For the same purpose, a model free reinforcement learning
approach is also used in [15]. Similarly, proportional-integral control is
used for ramp metering [16] for a discrete traffic flow model. In addition to
tracking a density setpoint, time and space dependent density trajectories are
tracked via boundary control in [7, 17] with the goal of removing excess cars
from the road and matching desired inflows and outflows.

Alternatively to ramp metering, VSL control relies on changing the speed
limits within a stretch of highway. The speed limit can be changed such that
the flow entering a traffic jam can be lessened or the flow leaving a traffic
jam can be increased so that the jam is mitigated faster. Both ramp metering
and VSL control are used in [18] to reduce on-ramp queues and total travel
time of vehicles within highway links. Model predictive control is also used
to reduce total travel time through both ramp metering and VSL control [19].
Both of these approaches use discrete traffic models. VSL control on a
continuous PDE is used in [20] by developing feedback control laws that
stabilize a desired density profile.

A possible way to realize VSL control that has not been explored is through a
linear quadratic regulator (LQR). In previous traffic applications, LQR has
been used in conjunction with non-cooperative Nash games to balance the flow
of vehicles [21]. This application, though, uses boundary feedback and is
developed on a discrete model, not a continuous traffic flow model. In other
hyperbolic PDE applications, LQR has been used as well. In [22], an infinite-
dimensional LQ controller is designed for boundary control of a system
containing a continuous stirred tank and a plug flow reactor. In [23], a state
feedback operator that controls the distributed jacket temperature of a
reactor is computed by solving a matrix Riccati differential equation. Both
applications have the same purpose of achieving a desired chemical
concentration.

Most in-domain VSL control of traffic flow has been studied on discretized
models [24, 9, 18, 19]. The approach in [20] was the first to apply VSL
continuously in space and time. While LQR for infinite dimensional systems has
been proven effective in other areas, its application to in-domain traffic
control poses a challenge because the rate of change of the control input is
paramount. The speed limit cannot change too quickly over space as that
creates unrealistic driving behavior such as large changes in the acceleration
of vehicles. This paper presents a novel LQR formulation where instead of
controlling the input, the rate of change is controlled for a continuous
hyperbolic PDE model. This approach also differs from previous continuous VSL
approaches [20] as, instead of controlling the VSL rate, the rate of change of
the speed limit over the length of highway is controlled.

##  II MODEL DESCRIPTION

###  II-A Lighthill-Whitham-Richards Model

The LWR model [5, 6] is a conservation of mass equation given by

$\frac{\partial\rho}{\partial t}+\frac{\partial q}{\partial z}=0$ (1)

where $\rho(z,t)$ is the traffic density and $q(z,t)$ is a fundamental diagram
which gives the equilibrium relationship between traffic density and flow. The
fundamental diagram is defined as

$q=\rho u$ (2)

where $u=U(\rho)$ is the equilibrium relationship between density and traffic
velocity. A common equilibrium relationship is Greenshield’s model which is a
decreasing function of density and is given as

$U(\rho)=U_{max}\bigg{(}1-\frac{\rho}{\rho_{max}}\bigg{)}$ (3)

where $U_{max}$ is the maximum velocity and $\rho_{max}$ is the maximum
density. The relationship between traffic flow and density using Greenshield’s
model is shown in Fig. 1. In Greenshield’s model, the critical density,
$\rho_{crit}$, which is where traffic goes from free flow to congested, is
half of $\rho_{max}$.

Figure 1: Fundamental diagram showing free flow region (green) and congested
region (red).

###  II-B Linearized LWR Model

The linearized LWR model is obtained by defining the perturbation variables as

$\rho=\rho_{0}+\Delta\rho$ (4) $u=u_{0}+\Delta
u=U_{max}\bigg{(}1-\frac{\rho_{0}+\Delta\rho}{\rho_{max}}\bigg{)}$ (5)

Then, by substituting (4),(5) into (1) the linear LWR model can be written as

$\frac{\partial\Delta\rho}{\partial t}+u_{0}\frac{\partial\Delta\rho}{\partial
z}+\rho_{0}\frac{\partial\Delta u}{\partial z}=0$ (6)

Expanding all the terms gives

$\frac{\partial\Delta\rho}{\partial
t}+U_{max}\bigg{(}1-\frac{\rho_{0}}{\rho_{max}}\bigg{)}\frac{\partial\Delta\rho}{\partial
z}+\\\ \rho_{0}\frac{\partial}{\partial
z}U_{max}\bigg{(}1-\frac{\Delta\rho}{\rho_{max}}\bigg{)}=0$ (7)

which simplifies to

$\frac{\partial\Delta\rho}{\partial
t}+U_{max}\bigg{(}1-2\frac{\rho_{0}}{\rho_{max}}\bigg{)}\frac{\partial\Delta\rho}{\partial
z}=0$ (8)

##  III CONTROL FORMULATION

###  III-A Variable Speed Limit Control

Variable speed limit control involves changing the speed limit in order to
control the flow of traffic [18, 9, 24, 20]. In this work, VSL is added to
Greenshield’s model to affect $U_{max}$ and thus $q$. The equilibrium velocity
(3) then becomes

$U_{VSL}(\rho)=bU_{max}\bigg{(}1-\frac{\rho}{\rho_{max}}\bigg{)}$ (9)

where $b$ is the VSL rate. The maximum speed $U_{max}$ acts as the speed limit
and the VSL rate $b$ works as the control input to either decrease or increase
the maximum speed. The definition of the control input $b$ introduces a new
nonlinearity into the model. Following the same procedure as before, perturbed
quantities for the rate of change of the speed limit and velocity are
introduced

$b=b_{0}+\Delta b$ (10) $u=(b_{0}+\Delta
b)U_{max}\bigg{(}1-\frac{\rho_{0}+\Delta\rho}{\rho_{max}}\bigg{)}$ (11)

The perturbation of density stays unchanged from (4). So, (7) now becomes

$\frac{\partial\Delta\rho}{\partial
t}+b_{0}U_{max}\bigg{(}1-\frac{\rho_{0}}{\rho_{max}}\bigg{)}\frac{\partial\Delta\rho}{\partial
z}+\rho_{0}\frac{\partial\Delta u}{\partial z}=0$ (12)

Inserting (11) into (12) gives

$\frac{\partial\Delta\rho}{\partial
t}+b_{0}U_{max}\bigg{(}1-\frac{\rho_{0}}{\rho_{max}}\bigg{)}\frac{\partial\Delta\rho}{\partial
z}+\\\
\rho_{0}U_{max}\bigg{(}1-\frac{\rho_{0}}{\rho_{max}}\bigg{)}\frac{\partial\Delta
b}{\partial
z}-b_{0}U_{max}\frac{\rho_{0}}{\rho_{max}}\frac{\partial\Delta\rho}{\partial
z}=0$ (13)

which simplifies to

$\frac{\partial\Delta\rho}{\partial
t}+b_{0}U_{max}\bigg{(}1-2\frac{\rho_{0}}{\rho_{max}}\bigg{)}\frac{\partial\Delta\rho}{\partial
z}+\\\
\rho_{0}U_{max}\bigg{(}1-\frac{\rho_{0}}{\rho_{max}}\bigg{)}\frac{\partial\Delta
b}{\partial z}=0$ (14)

The variable $\frac{\partial\Delta b}{\partial z}$ is the change in the
perturbation of the VSL rate over the length of road. When there is no VSL,
$\frac{\partial\Delta b}{\partial z}=0$ and $b_{0}=1$, and (14) returns to
(8).

###  III-B Optimal Control Design for Hyperbolic PDE model

The first step to derive the state feedback control policy is to write the PDE
in its equivalent state space formulation on a Hilbert space $\mathcal{H}$
[23]. For a linear hyperbolic PDE of the form

$\displaystyle\frac{\partial x}{\partial t}$ $\displaystyle=V\frac{\partial
x}{\partial z}+Mx+B_{0}u$ (15) $\displaystyle y$ $\displaystyle=C_{0}x$ (16)

this results in

$\displaystyle\dot{x}(t)$ $\displaystyle=Ax(t)+Bu(t)$ (17) $\displaystyle
y(t)$ $\displaystyle=Cx(t)$ (18)

where $A$ is the linear operator defined as

$A=V\cdot\frac{d.}{dz}+M\cdot I$ (19)

on the domain

$D(A)=\\{x\in\mathcal{H}:\frac{dx}{dz}\in\mathcal{H}\textrm{ and }x(0)=0\\}$
(20)

with $V\in\mathbb{R}^{n\times n}$ being a symmetric matrix and $M$, $B_{0}$,
and $C_{0}$ being real continuous space-varying matrix functions. In addition,
$x$ is assumed to be absolutely continuous. It is proven in [23] that if $V<0$
then $A$ generates an exponentially stable $C_{0}$-semigroup.

The linear-quadratic optimal control problem on an infinite time interval
finds a control input $u_{opt}$ such that the functional $J(x_{0},u)$ is
minimized. The cost functional is given as

$J(x_{0},u)=\int_{0}^{\infty}\langle Cx(t),QCx(t)\rangle+\langle
u(t),Ru(t)\rangle\,dt$ (21)

The solution of (21) can be found by finding the non-negative self-adjoint
operator $P$ which solves

$[A^{*}P+PA+C^{*}QC-PBR^{-1}B^{*}P]x=0$ (22)

When $(A,B)$ is exponentially stabilizable and $(Q^{\frac{1}{2}}C,A)$ is
exponentially detectable (22) has a unique, non-negative solution $P$ and (21)
is minimized by the unique control input $u_{opt}$ given by

$u_{opt}(t)=K_{0}x(t)$ (23)

where the optimal feedback is

$K_{0}=-R^{-1}B^{*}P$ (24)

The proof that $(A,B)$ is exponentially stabilizable and
$(P^{\frac{1}{2}}C,A)$ is exponentially detectable given $V<0$ is shown in
[23]. Furthermore,

$P:=\Phi(z)I$ (25)

where $\Phi(z)$ is the solution of the matrix Riccati differential equation
given by

$\displaystyle\begin{split}V\frac{d\Phi}{dz}&=M^{*}\Phi+\Phi
M+C_{0}^{*}Q_{0}C_{0}-\Phi B_{0}R_{0}^{-1}B_{0}^{*}\Phi\\\
\Phi(L)&=0\end{split}$ (26)

where $Q_{0}$ is a positive matrix and $R_{0}$ is a self-adjoint positive
matrix. Thus, (24) becomes

$K_{0}=-R_{0}^{-1}(z)B^{*}_{0}(z)\Phi(z)I$ (27)

###  III-C LQ Control Deisgn for LWR Model with VSL

In order to find $\Phi(z)$, we take (14) and transform it into

$\frac{\partial\Delta\rho}{\partial
t}=-b_{0}U_{max}\bigg{(}1-2\frac{\rho_{0}}{\rho_{max}}\bigg{)}\frac{\partial\Delta\rho}{\partial
z}-\\\
\rho_{0}U_{max}\bigg{(}1-\frac{\rho_{0}}{\rho_{max}}\bigg{)}\frac{\partial\Delta
b}{\partial z}$ (28)

so that it is now in the form of (15). It follows that

$\displaystyle\begin{split}V&=-b_{0}U_{max}\bigg{(}1-2\frac{\rho_{0}}{\rho_{max}}\bigg{)}\\\
M&=0\\\
B_{0}&=-\rho_{0}U_{max}\bigg{(}1-\frac{\rho_{0}}{\rho_{max}}\bigg{)}\\\
C_{0}&=I\end{split}$ (29)

The condition of $V<0$ is met if $\rho_{0}<\frac{1}{2}\rho_{max}$ in (29).
This means that $\rho_{0}$ must always fall in the free flow regime, which is
a valid assumption since the equilibrium profile should never fall in the
congested regime to avoid traffic jams. The control input for the LWR model
with VSL is the change of the speed limit over the stretch of highway. Using
(26) and (29) the matrix Riccati equation becomes

$\displaystyle\begin{split}V\frac{d\Phi}{dz}&=Q_{0}-\Phi
B_{0}R_{0}^{-1}B_{0}^{*}\Phi\\\ \Phi(L)&=0\end{split}$ (30)

As in [23], $R_{0}=1$. The values for $V$ and $B_{0}$ will be inserted into
the equation at the end of the derivation for $\Phi$ for simplicity. Using
separation of variables, (30) becomes

$\int_{0}^{\Phi}\frac{V}{Q_{0}-B_{0}^{2}\Phi^{2}}\,d\Phi=\int_{0}^{z}\,dz=z+C$
(31)

The left hand side of (31) can be solved by a combination of integration by
partial fractions and substitution. After applying both methods, (31) becomes

$\frac{V}{2B_{0}\sqrt{Q_{0}}}(\ln|B_{0}\Phi+\sqrt{Q_{0}}|-\ln|B_{0}\Phi-\sqrt{Q_{0}}|)=z+C$
(32)

Then,

$\displaystyle\begin{split}\ln\bigg{|}\frac{B_{0}\Phi+\sqrt{Q_{0}}}{B_{0}\Phi-\sqrt{Q_{0}}}\bigg{|}&=\frac{2B_{0}\sqrt{Q_{0}}}{V}(z+C)\\\
\bigg{|}\frac{B_{0}\Phi+\sqrt{Q_{0}}}{B_{0}\Phi-\sqrt{Q_{0}}}\bigg{|}&=\exp\bigg{(}\frac{2B_{0}\sqrt{Q_{0}}}{V}(z+C)\bigg{)}\\\
\big{|}B_{0}\Phi\\!+\\!\sqrt{Q_{0}}\big{|}\\!&=\\!\big{|}B_{0}\Phi\\!-\\!\sqrt{Q_{0}}\big{|}\exp\bigg{(}\frac{2B_{0}\sqrt{Q_{0}}}{V}(z\\!+\\!C)\bigg{)}\end{split}$
(33)

The solution to (33) can be either

$B_{0}\Phi+\sqrt{Q_{0}}=(B_{0}\Phi-\sqrt{Q_{0}})\exp\bigg{(}\frac{2B_{0}\sqrt{Q_{0}}}{V}(z+C)\bigg{)}$
(34)

or

$B_{0}\Phi+\sqrt{Q_{0}}=-(B_{0}\Phi-\sqrt{Q_{0}})\exp\bigg{(}\frac{2B_{0}\sqrt{Q_{0}}}{V}(z+C)\bigg{)}$
(35)

Because of the condition that $\Phi(L)=0$, only (35) can be the solution for
$\Phi$. The final solution for the state feedback function $\Phi$ using the
condition $\Phi(L)=0$ is

$\Phi=\frac{\sqrt{Q_{0}}\bigg{(}\exp(\frac{2B_{0}\sqrt{Q_{0}}}{V}(z-L))-1\bigg{)}}{B_{0}\bigg{(}\exp(\frac{2B_{0}\sqrt{Q_{0}}}{V}(z-L))+1\bigg{)}}$
(36)

Combining (23), (27), and (36) with $x(t)=\Delta\rho$ for the linear model
leaves us with the final solution for the optimal control input

$u_{opt}(z,t)\\!=\\!\frac{\partial b}{\partial
z}\\!=\\!-\sqrt{Q_{0}}\frac{\bigg{(}\exp(\frac{2B_{0}\sqrt{Q_{0}}}{V}(z-L))-1\bigg{)}}{\bigg{(}\exp(\frac{2B_{0}\sqrt{Q_{0}}}{V}(z-L))+1\bigg{)}}\Delta\rho$
(37)

The parameter $Q_{0}$ is left as a tuning parameter and its influence on the
solution of $\Phi$ is shown in Fig. 2.

Figure 2: Influence of $Q_{0}$ on the state feedback function, $\Phi$.

##  IV SIMULATION RESULTS

###  IV-A Uncontrolled Scenario

The case study investigated in this paper simulates a stretch of highway with
an oscillating and increasing flow entering it. This situation could occur
because of both on-ramp flow and traffic light operation. Because the base
model used for the LQ control is linear, the simulation is kept within the
free flow region so that no traffic jams, or shocks in the system, occur. The
initial density on the stretch of highway is

$\rho(0,z)=10\sin\bigg{(}\frac{\pi z}{L}\bigg{)}+\rho_{0}$ (38)

and the upstream boundary condition is

$\rho(t,0)=5\exp(-L\cdot 10^{-6}t)\sin\bigg{(}\frac{\pi
t}{20}\bigg{)}+\rho_{0}+\frac{t}{4L}$ (39)

The model parameters for the simulation are given in Table I. Based off of the
fundamental diagram shown in Fig. 1 and parameters in Table I, this means that
the the density across the highway must be kept below $\rho_{crit}=80$
cars/km. The resulting density evolution and velocity of the highway section
is shown in Fig. 3.

TABLE I: Model Parameters. Parameter | Value | Unit
---|---|---
Maximum density $\rho_{max}$ | 160 | [cars/km]
Maximum speed $U_{max}$ | 115 | [kph]
Average density $\rho_{0}$ | 50 | [cars/km]
Road length $L$ | 2000 | [m]
Simulation time $T$ | 120 | [s]

The results of the baseline simulation using the linear model are shown in
Fig. 3(a),(b), and the results of the baseline nonlinear simulation are shown
in Fig. 3(c),(d). For both models, the baseline scenario results stay within
the free flow regime as the maximum density reached is less than
$\rho_{crit}$. There is good agreement between the results of the linear and
nonlinear model. As the density at the beginning of the length of road
approaches $\rho_{crit}$ at around 100 s, the solutions between the two models
start to differ.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 3: Baseline results for linear model (a), (b) and nonlinear model (c),
(d) given the initial and boundary conditions in (38), (39).

###  IV-B Linear LQR Results

The performance of the LQ controller is first investigated using the linear
LWR model. In order to evaluate the performance, the total number of cars in
the highway section is calculated as

$\\#\textrm{ of cars}=\int_{0}^{L}\rho\,dz$ (40)

As the average desired density is 50 cars/km, the desired amount of cars left
in the highway section is 100 cars. The results of applying the developed LQ
controller to the linear model are shown in Fig. 4. Four different values of
the tuning parameter $Q_{0}$ are investigated. For $Q_{0}=1e-6$, the control
is not enough to lower the density by a reasonable amount. Towards the end of
the simulation, at around 100 s, the density even starts to rise. For both
$Q_{0}=1e-5$ and $Q_{0}=5e-5$ the control is enough to lower the density, but
at a certain point it no longer drives the density to the desired amount of
cars left on the road. For $Q_{0}=1e-5$ this stagnation happens at 80 s and
for $Q_{0}=5e-5$ this happens at around 40 s. One thing to note, though, is
that when $Q_{0}=5e-5$ the amount of cars does approach the desired amount,
but never reaches it. The controller is able drive the linear system to the
desired density when $Q_{0}=5e-4$, and it reaches it at around 20 s.

Figure 4: Total cars on highway section for varying values of $Q_{0}$ when
control is applied to linear model.

###  IV-C Nonlinear LQR Results

For the nonlinear model, the control input is no longer $\frac{\partial
b}{\partial z}$, instead it is the actual VSL rate, $b$. So, in order to apply
the control to the nonlinear model, first the linear model is used to
determine the optimal $\frac{\partial b}{\partial z}$, then the nonlinear
control is found by

$b_{opt}=\int_{0}^{L}\frac{\partial b}{\partial z}\,dz$ (41)

The nonlinear model (1) then uses (9) for the equilibrium velocity equation
where $b$ is defined as in (41). The results for the same values of $Q_{0}$ as
used before, but this time on the nonlinear model, are shown in Fig. 5. For
$Q_{0}=1e-6$, the system is mostly unaffected. When $Q_{0}=1e-5$ and
$Q_{0}=5e-5$ the traffic is reduced, but it does not reach the desired
density. With $Q_{0}=5e-4$, the controller is able to get close to the desired
density by 40 s, but does not achieve it.

Figure 5: Total cars on highway section for varying values of $Q_{0}$ when control is applied to nonlinear model. |  |
---|---|---
(a) | (b) | (c)
Figure 6: Density (a), velocity (b), and variable speed limit (c) of linear model with $Q_{0}=5e-5$. |  |
---|---|---
(a) | (b) | (c)
Figure 7: Density (a), velocity (b), and variable speed limit (c) of nonlinear
model with $Q_{0}=5e-5$.

###  IV-D Comparison

As shown in both Fig. 4 and Fig. 5, the performance of the controller differs
between the linear and nonlinear implementation. In each case, for the same
value of $Q_{0}$, the linear model is affected more than the nonlinear model.
For example, when $Q_{0}=5e-4$ the controller is able to drive the system to
the desired traffic density in the linear model, but cannot do it for the
nonlinear model. Similarly, the linear model reaches the desired density at 20
s, whereas in the nonlinear model it does not get close until 40 s.

To further investigate the performance of the designed controller on the two
models Fig. 6 and Fig. 7 show the resulting density, velocity, and speed limit
for the linear and nonlinear model, respectively. The results shown are for
$Q_{0}=5e-5$. Looking at Fig. 6(a) and Fig. 7(a), it can be seen that the
nonlinear controller is not as effective. The high density region resulting
from the initial condition at $x=500$ m to $x=1500$ m is dissipated in the
linear model, and it never reaches the downstream boundary. On the other hand,
in the nonlinear model, the high density region continues to move until it
reaches the downstream boundary at around 20 s. As well, the upstream boundary
conditions propagate further into the road in the nonlinear simulation, but in
the linear simulation the controller is able to mitigate them quickly.

Another difference can be seen in the speed limits for the linear and
nonlinear model, Fig. 6(c) and Fig. 7(c), respectively. In the nonlinear
simulation, the controller needs to keep a higher speed limit for longer at
the end of the length of road. The high speed limit lasts until just before 20
s in the nonlinear model whereas in the linear model the high speed limit
region ends after around 5 s.

As explained earlier, the control input for the linear model is the rate of
change of the VSL rate over space $\frac{\partial b}{\partial z}$, but in the
nonlinear model it is the actual VSL rate, $b$. This could explain the lag in
response of the nonlinear model when compared to the linear model.

##  V CONCLUSIONS

In this paper, a LQ controller was designed for the LWR model using variable
speed limit control. First, the linearized LWR model was defined and a
variable speed limit was applied to the model as a modification of
Greenshield’s model. Then, a state feedback function was computed via the
solution of a Riccati differential equation. In this application, instead of
the control input being just a distributed parameter, the control input was
the rate of change of the distributed parameter. The designed controller was
then verified on both the linear and nonlinear LWR model, and the performance
of the two controllers was compared. Future work will investigate how to tune
$Q_{0}$ as well as how to develop a controller that can be used in either a
mixed traffic or congested traffic scenario to mitigate traffic jams in the
presence of shocks.

## References

  * [1] D. Schrank, T. Lomax, and B. Eisele, “2021 urban mobility report,” Texas Transportation Institute, 2021. [Online]. Available: http://mobility.tamu.edu/ums/report
  * [2] B. Pishue, “2022 INRIX Global Traffic Scorecard,” in _INRIX (January 2023)_ , 2023.
  * [3] G. Cookson and B. Pishue, “2016 INRIX Global Traffic Scorecard,” in _INRIX (January 2017)_ , 2017.
  * [4] M. Treiber and A. Kesting, _Traffic Flow Dynamics: Data, Models and Simulation_. Springer Science & Business Media, 2012.
  * [5] M. J. Lighthill and G. B. Whitham, “On kinematic waves ii. a theory of traffic flow on long crowded roads,” _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ , vol. 229, no. 1178, pp. 317–345, 1955.
  * [6] P. I. Richards, “Shock waves on the highway,” _Operations research_ , vol. 4, no. 1, pp. 42–51, 1956.
  * [7] L. Tumash, C. Canudas-de Wit, and M. L. Delle Monache, “Robust tracking control design for fluid traffic dynamics,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_. IEEE, 2019, pp. 4085–4090.
  * [8] V. Knoop, _Introduction to Traffic Flow Theory: An introduction with exercises_. TU Delft Open, 2017.
  * [9] R. C. Carlson, I. Papamichail, and M. Papageorgiou, “Local feedback-based mainstream traffic flow control on motorways using variable speed limits,” _IEEE Transactions on intelligent transportation systems_ , vol. 12, no. 4, pp. 1261–1276, 2011.
  * [10] A. Spiliopoulou, M. Kontorinaki, M. Papageorgiou, and P. Kopelias, “Macroscopic traffic flow model validation at congested freeway off-ramp areas,” _Transportation Research Part C: Emerging Technologies_ , vol. 41, pp. 18–29, 2014.
  * [11] H. Yu, M. Diagne, L. Zhang, and M. Krstic, “Bilateral boundary control of moving shockwave in LWR model of congested traffic,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 3, pp. 1429–1436, 2020.
  * [12] H. Liu, C. Claudel, and R. Machemehl, “Robust traffic control using a first order macroscopic traffic flow model,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 23, no. 7, pp. 8048–8062, 2021.
  * [13] H. Yu and M. Krstic, “Adaptive output feedback for aw-rascle-zhang traffic model in congested regime,” in _2018 Annual American Control Conference (ACC)_. IEEE, 2018, pp. 3281–3286.
  * [14] ——, “Traffic congestion control for Aw-Rascle-Zhang model,” _Automatica_ , vol. 100, pp. 38–51, 2019.
  * [15] H. Yu, S. Park, A. Bayen, S. Moura, and M. Krstic, “Reinforcement learning versus PDE backstepping and PI control for congested freeway traffic,” _arXiv preprint arXiv:1904.12957_ , 2019.
  * [16] C. Pasquale, S. Sacone, and S. Siri, “Closed-loop stability of freeway traffic systems with ramp metering control,” in _2018 IEEE Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 223–228.
  * [17] L. Tumash, C. C. de Wit, and M. L. Delle Monache, “Boundary control design for traffic with nonlinear dynamics,” _IEEE Transactions on Automatic Control_ , 2021.
  * [18] I. Papamichail, K. Kampitaki, M. Papageorgiou, and A. Messmer, “Integrated ramp metering and variable speed limit control of motorway traffic flow,” _IFAC Proceedings Volumes_ , vol. 41, no. 2, pp. 14 084–14 089, 2008.
  * [19] S. Liu, H. Hellendoorn, and B. De Schutter, “Model predictive control for freeway networks based on multi-class traffic flow and emission models,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 18, no. 2, pp. 306–320, 2016.
  * [20] I. Karafyllis and M. Papageorgiou, “Feedback control of scalar conservation laws with application to density control in freeways by means of variable speed limits,” _Automatica_ , vol. 105, pp. 228–236, 2019.
  * [21] D. Pisarski and C. Canudas-de Wit, “Nash game-based distributed control design for balancing traffic density over freeway networks,” _IEEE Transactions on Control of Network Systems_ , vol. 3, no. 2, pp. 149–161, 2015\.
  * [22] A. A. Moghadam, I. Aksikas, S. Dubljevic, and J. F. Forbes, “Lq control of coupled hyperbolic pdes and odes: Application to a cstr-pfr system,” _IFAC Proceedings Volumes_ , vol. 43, no. 5, pp. 721–726, 2010.
  * [23] I. Aksikas, A. Fuxman, J. F. Forbes, and J. J. Winkin, “Lq control design of a class of hyperbolic pde systems: Application to fixed-bed reactor,” _Automatica_ , vol. 45, no. 6, pp. 1542–1548, 2009.
  * [24] R. Carlson, “Mainstream traffic flow control on motorways,” Ph.D. dissertation, Ph. D. dissertation, Tech. Univ. Crete, Chania, Greece, 2011.