# Federated Learning over Next-Generation Ethernet Passive Optical Networks

Oscar J. Ciceri, Carlos A. Astudillo, Zuqing Zhu, and Nelson L. S. da Fonseca
Institute of Computing - University of Campinas
Campinas 13089-971, SP, Brazil
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Oscar J. Ciceri, Carlos A. Astudillo, Zuqing Zhu, and
Nelson L. S. da Fonseca O. J. Ciceri, C. A. Astudillo and N. L. S da Fonseca
are with the Institute of Computing, University of Campinas 13083-852, Brazil
(emails: oscar<EMAIL_ADDRESS>nfonseca@ic.unicamp.br).Z. Zhu is with the School of Information Science and
Technology, University of Science and Technology of China, Hefei, Anhui
230027, P. R. China (email: zqzhu@ieee.org).

###### Abstract

Federated Learning (FL) is a distributed machine learning (ML) type of
processing that preserves the privacy of user data, sharing only the
parameters of ML models with a common server. The processing of FL requires
specific latency and bandwidth demands that need to be fulfilled by the
operation of the communication network. This paper introduces a Dynamic
Wavelength and Bandwidth Allocation algorithm for Quality of Service (QoS)
provisioning for FL traffic over 50 Gb/s Ethernet Passive Optical Networks.
The proposed algorithm prioritizes FL traffic and reduces the delay of FL and
delay-critical applications supported on the same infrastructure.

###### Index Terms:

Dynamic Wavelength and Bandwidth Allocation, Federated Learning, Ethernet
Passive Optical Networks, Quality of Service, Multi-services, Future Access
Networks.

##  I Introduction

Artificial Intelligence (AI) applications are widely employed in various
different sectors of our daily life, such as in smart cities, smart homes,
smart transportation, smart health, and finance. Moreover, the upcoming 6th
Generation mobile networks (6G) will support an even wider range of
applications based on AI, such as holographic communications, super-smart
homes, and cooperative autonomous robots. The 6G system will change the
network vision from ”connected things” to ”connected intelligence” [1], in
which advanced and ubiquitous AI will empower numerous applications.

Traditional use of Machine Learning (ML) models relies on batch processing in
a central server and the employment of datasets containing user data. With the
worldwide adoption of data protection and privacy legislation, the creation of
datasets and applications based on ML has been considerably limited. One way
of coping with such restrictions is the adoption of Federated Learning (FL),
which is a distributed way of processing machine learning algorithms that does
not disclose private data. In FL, clients train a local ML model using a
private dataset, and the parameters of these local models are then sent to a
central server. The server produces a global model on the basis of the
numerous parameter values received and distributes this global model to the
clients for further training. This round of processing is repeated until the
global model produces results with an acceptable level of accuracy. In this
way, user privacy is preserved. The most common approaches for the
consolidation of the parameters sent by the clients to produce the global
model (FedAvg) rely on the assumption that clients are synchronized and that
local datasets are independent and identically distributed [2]. Such FL can be
used for the training of large ML models involving thousands or even millions
of parameters.

The processing of Federated Learning models has brought challenges to
communication networks. FL clients may produce highly bursty traffic when
uploading their model parameters to the server. For example, clients training
a Convolutional Neural Networks model with few convolution layers and
thousands of parameters may need to send hundreds of megabytes. When millions
of parameters are involved, the amount of bytes sent can be of the order of
gigabytes [3]. Moreover, FL may impose stringent communication delays for the
uploading of client parameters to enhance fast convergence to the global
model, especially when the federation involves numerous clients. To cope with
diverse communication delays, the server may either wait the arrival of the
local parameters from all clients, increasing convergence time, or exclude the
late arriving data from the parameter consolidation step, which reduces the
accuracy of the model [4]. In addition, FL may also require a very large
number of training rounds to produce accurate global models [5]. These
challenges calls for efficient resource allocation mechanisms to meet the FL
requirements.

Passive Optical Network (PON) is a cost-efficient access network technology
for delivering broadband services [6]. Operators have already deployed 10 Gb/s
Time Division Multiplexing (TDM) PONs during the past two decades. In recent
years, the ITU and IEEE standardization groups have proposed next-generation
PONs based on Time and Wavelength Division Multiplexing (TWDM) to increase the
network capacity for supporting demanding applications and services. TWDM
allows allocation in various wavelength channels of 25 Gb/s (50G-EPON) and 10
Gb/s (40G-XPON) [7].

Infrastructure service Providers owning a PON can increase network
utilization, as well as their profits, by offering a diverse variety of
services to a variety of different customers, such as residential users,
enterprises, and mobile network operators (MNO). The capacity to support the
Quality of Service (QoS) requirements of emerging 5G/6G applications, such as
distributed Machine Learning, Tactile Internet, or Mobile Fronthauling (MFH)
calls for efficient use of network resources. Nonetheless, the aforementioned
challenges of FL applications make QoS provisioning challenging.

A few approaches have been proposed to deal with FL processing over PONs ([4,
8]). One of these is the Bandwidth Slicing (BS) approach for TDM-PONs, that
reserves portions of the PON capacity for FL clients [4]. The bandwidth for
each of the slices granted per cycle is orders of magnitude less than that
required for transmitting a model update, which implies that several
scheduling cycles will be required to fully upload the parameters of the
client models. An architecture for scalable federated learning involving two-
step of aggregation over PONs was introduced in [8]. In this proposal, the
parameters of local models are first aggregated at clients connected to
Optical Network Unit (ONU) and then aggregated on a server connected to the
Optical Line Terminator (OLT). As a consequence, the amount of upstream
traffic remains relatively constant regardless of the number of clients in the
federation. However, the high volume of traffic generated by the FL clients
can create network bottlenecks, which impacts on the time to upload of the
local parameters.

This paper discusses the main issues for supporting FL applications over PONs
and introduces a Dynamic Wavelength and Bandwidth Allocation (DWBA) algorithm
for 50G-EPONs that dynamically prioritizes FL traffic while maintaining the
traditional guaranteed bandwidth scheme for PON customers. Two prioritization
policies have been proposed to reduce the delay of FL traffic and delay-
critical applications. In the first, the intra-ONU scheduler strictly
prioritizes the FL traffic over that from other types of applications (FL-
first policy). In the second, the delay-critical traffic is prioritized over
the FL traffic (DC-first policy). The BS algorithm introduced by Li et. al
(2020) [4] was adapted to employ multiple wavelengths, as well as an adaptive
polling cycle as required in 50G-EPON networks with dynamic resource
allocation for comparison purpose.

Results show that the DC-first policy increases the FL model accuracy and
reduces the delay of federated learning and delay-critical applications when
compared to the BS approach and the FL-first policy.

##  II Resource Allocation in Passive Optical Networks

PON is a network access technology that offers larger capacity, greater cost-
efficiency, and more energy savings than do other network access technologies.
There are two main PON standards: Ethernet PON (EPON) and Gigabit Capable
(GPON), with EPON being less expensive. GPON transmission system ofemploys
synchronous frames issued at every 125 $\mu$s, while those of EPON use
Ethernet frames asynchronously for transmissions based on granted cycles of
variable duration. While traditional PON standards allow bit rates of 1
$\nicefrac{{Gb}}{{s}}$ and 10 $\nicefrac{{Gb}}{{s}}$, the next-generation PON
standards allow those of 40 $\nicefrac{{Gb}}{{s}}$ to 100
$\nicefrac{{Gb}}{{s}}$.

The 50 Gb/s optical access network standardized in IEEE 50G-EPON 802.3ca-2020
[9] is a promising technology for adoption by InP to support emerging services
with strigent latency and bandwidth requirements. This 50G-EPON technology
employs the Time and Wavelength Division Multiple Access (TWDMA) technique for
controlling uplink transmissions between the ONUs and the OLT. There are three
main TWDM-PON-based access architectures for the connectivity between the OLT
and ONUs: Multiple-Scheduling Domain (MSD), Single-Scheduling Domain (SSD),
and Wavelength Agile (WA). In the first, ONUs transmit on a single wavelength
at a time. In the second, ONUs can transmit simultaneously on all wavelengths,
and in the third, more than one wavelength can be granted to a single ONU.

In this technology, the signaling protocol Multipoint Control Protocol (MPCP)
is employed for resource allocation. This protocol uses Report and Gate
messages for this allocation. Report messages are sent on the upstream to the
OLT by the ONUs to request bandwidth, while Gate messages are sent on the
downstream by the OLT to the ONUs to inform the grated wavelength(s) and
transmission windows, as well as the starting time of the next transmission
window. Resource allocation is carried out in two steps, one for wavelength
allocation and the other for bandwidth allocation. The use of different
schemes for transmission on multiple wavelengths can be defined on the basis
of conventional Dynamic Bandwidth Allocation (DBA) algorithms for TDM-PONs.

For dynamic bandwidth allocation over EPONs, the Interleaved Polling with
Adaptive Cycle Time (IPACT) algorithm has been adopted to complement the MPCP
protocol. The IPACT algorithm employs an interleaved polling and statistical
multiplexing technique that leads to efficient upstream channel usage. The
limited policy has been used to assure bandwidth to ONUs according to pre-
defined Service Level Agreements. Moreover, the original IPACT algorithm
employs a single wavelength channel for scheduling. It has been modified to
operate with multiple wavelengths in [10], [11] and [12]. The modified IPACT
algorithm was proposed for the SSD and MSD architectures [10]. Additional
algorithms have been proposed: the Water-Fill (WF) [11] to promote fairness in
the wavelength utilization and First-Fit (FF) [12] to provide less delay.
Moreover, when there is no scheduler for Quality of Service (QoS) provisioning
in the PON, the First-Come-First-Served (FCFS) queuing policy is employed.
However, this strategy does not consider the priority or required
bandwidth/delay of the applications.

Figure 1: An overview of federated learning over passive optical networks.

The performance of diverse applications in a PON is ensured by the adopted QoS
mechanism, it controls the way frames are queued, prioritized, and scheduled.
Such assurance of QoS can be provided by either the ONU or OLT. In the single-
level architecture, the ONUs reports individual queue sizes, while the OLT
distributes the bandwidth for each type of traffic. In the hierarchical
architecture, the OLT allocates bandwidth for each ONU, and the ONUs manage
the amount of bandwidth to be allocated for each traffic.

The most straightforward method of facilitating QoS provisioning is the
Differentiated Service approach. It classifies network traffic and delivers
different services to different applications. The simplest way to implement
Differentiated Services is to employ strict priority scheduling. The ONUs
categorize the incoming traffic and put it in the buffer, imposing the
prioritization of different traffic classes. With the employment of
Differentiated Services, however, the PON can support packetized voice and
video with strict bandwidth and latency constraints, as well as best-effort
traffic [13]. However, none of the exiting mechanisms have been specially
designed to support the QoS requirements of FL applications.

##  III Resource Allocation for Federated Learning

Fig. 1 illustrates the scenario of FL processing over a PON; where clients are
connected to the ONUs and the server is attached to the OLT.

The FL training process can be either asynchronous or synchronous. In the
former, the global model parameters are computed as soon as the server
receives updates of the parameters of the local models from a certain number
of clients. In the latter, the server aggregates the local parameters that
arrive in a period of constant duration excluding the parameters from the late
arriving stragglers. This exclusion, however, reduces the accuracy of the
model as well as increases the time required for convergence to a final global
model.

The synchronization time per round includes the downstream, computing,
network, propagation, and aggregation delays, as shown in Fig. 2a. The
downstream delay includes the propagation and transmission delays of the
parameters of the global model from the server to the clients. The computing
time is the time taken to train the local model at the client in each round.
The network delay is the time spent in communication the local model
parameters from the clients to the server. The propagation delay is the time
taken by the bits transmitted to travel from the client to the server on the
network medium. The aggregation delay is the processing time of the
aggregation algorithm.

The network delay depends on the network load and the resource allocation
mechanism employed for the PON for allocating bandwidth and wavelength(s) to
the ONUs. The computing time depends on the capacity of the client and the
size of the training dataset, while the propagation delay depends on the
distances between the server and the clients. Moreover, long network delays
may increase the number of straggler clients, decreasing the model accuracy
and increasing the convergence time.

The time taken to transmit the local model parameters to the server depends on
the bandwidth allocated to the FL traffic. In general, PON customers receive a
portion of the total available bandwidth in the PON due to the shared nature
of the upstream channel. Residential and business customers usually have
guaranteed bandwidth from tens to hundreds of Mbps, and other PON customers
can require on demand up to tens of Gbps. However, the large size of the local
model parameters, which can be in the order of gigabytes, may demand several
seconds to be fully transmitted, even with guaranteed bandwidths on the order
of Gbps, as shown in Figure 2b.

The unique characteristics of FL applications introduce challenges for the
management of the available bandwidth in scenarios with limited bandwidth and
diversity of customers, services and applications. The recently proposed BS
algorithm to the support of QoS for FL applications assures a slice of
bandwidth for the FL traffic. This is then allocated according to the
ascending order of downstream client delay and computing time. However, the
approach is not efficient for PON scenarios when loads are high and traffic is
bursty. The large bandwidth slice required to serve FL traffic properly
reduces the available bandwidth for other PON customers. It also reduces the
statistical multiplexing gain due to the static allocation of the network
resources. On the other hand, if the slice is much smaller than the total
available bandwidth, the granted bandwidth is likely to be insufficient to
serve timely the FL processing due to the large size of the FL packets. As a
consequence, clients send a small portion of the FL frames per cycle,
requiring numerous scheduling cycles to be fully served. If two or more
clients have FL frames in the queue at a given time, then only one client can
use the slice, while the others will have to wait for that slice to become
available, as shown in Figure 2c.

Furthermore, the BS approach is not compatible with traditional PON business
models, in which customers rent portions of the PON capacity from the InP to
support their applications and services. Moreover, it is unclear who pays for
a shared bandwidth slice.

In summary, even though the BS approach reduces the latency for FL
applications in relation to the traditional First Come First Served approach,
the static bandwidth allocation and the incompatibility of the approach with
traditional business models may lead to issue of QoS support and deployment.

(a) Delays of Federated Learning

(b) Large FL Packet Size

(c) Bandwidth Slice for FL

Figure 2: Resource allocation issues for Federated Learning over PONs

##  IV DWBA Schemes for Supporting FL traffic over 50G-EPON networks

This section introduces the proposed scheme for providing QoS for FL
applications over 50G-EPON networks.

###  IV-A Adaptation of the Bandwidth Slicing Approach to TWDM-EPON

The BS algorithm [4] serves FL traffic by dynamically allocating bandwidth on
a single wavelength for FL clients. It calculates the number of cycles an FL
client requires to be completely served based on the required bandwidth and
the fixed polling cycle length of 125 $\mu s$ employed in the GPON technology.
However, next-generation PONs employ multiple wavelengths and the cycle
duration is unknown a priori when an adaptive polling cycle mechanism is
employed, such as in the EPON technology. Therefore, the proposed BS algorithm
cannot be directly employed in TWDM-EPON networks.

We proposed an adaptation of the BS approach for TWDM-EPONs called multi-
wavelength BS algorithm (MW-BS). It deals with multiple wavelengths and
employs an adaptive polling cycle for dynamic resource allocation. A portion
of the PON capacity (bandwidth slice) is still reserved for the FL traffic in
each scheduling cycle, but instead of using the polling cycle information to
grant the transmission windows for FL traffic and then share the remaining
slice capacity with other traffic types, MW-BS reserves the total slice to the
FL traffic as long as a bandwidth request from any FL client exists. The use
of a dynamic polling cycle reduces the FL traffic delays and avoids the need
for information about the duration of the unknown upcoming cycles. Three
variations of the MW-BS algorithm were proposed for different TWDM wavelength
allocation policies, namely MW-BS-SSD, MW-BS-MSD, and MW-BS-FF.

Figure 3: Flow chart of the bandwidth slicing algorithm adapted for TWDM-
EPONs.

The flow diagram in Figure 3 summarizes the proposed DWBA scheme residing on
the OLT. The ONUs sends the Report message requesting bandwidth for Federated
Learning as well as conventional applications. When a Report message arrives
from a ONU containing a bandwidth request for FL traffic, the OLT first grants
the bandwidth from the reserved slice, if available (Block 1); otherwise, the
OLT allocates bandwidth for the conventional traffic (Block 2).

For allocation of bandwidth for the slice, The OLT also reserves bandwidth for
the ONU for upcoming cycles. It selects the wavelength(s) as a function of the
TWDM wavelength allocation policy, and calculates the next starting time for
the FL transmission. For the SSD policy, the OLT grants all wavelengths. For
the MSD policy, the OLT grants a predetermined-fixed wavelength. For the FF
policy, the OLT grants the first available wavelength, and then calculates the
granted transmission window for each allocated wavelength, depending on the
number of wavelengths and the portion of the PON capacity designed for FL use.
If the granted window is equal to the requested windows, the FL traffic will
be fully served, and the OLT will make the bandwidth slice available for the
next cycle.

The OLT also calculates the granted bandwidth for the convectional
applications. If the OLT has previously allocated the bandwidth for the slice,
it selects these wavelength(s) for the FL traffic. Otherwise, the OLT selects
the wavelength(s) for the FL traffic and calculates the next starting time for
that FL transmission as a function of the TWDM wavelength allocation policies
involved. The transmission window for the conventional applications is
calculated according to the limited policy. Finally, the OLT issues and send a
Gate message with the granted bandwidths for both FL and conventional
applications.

###  IV-B DWBA for Federated Learning

To address the issues of resource allocation raised here, we introduce a DWBA
algorithm that supports QoS provisioning for Federated Learning applications
while meeting the requirements of delay-critical applications in TWDM-EPON
networks. The algorithm is called DWBA for Federated Learning (DWBA-FL).

The idea behind our proposal is to allow PON customers to employ their
guaranteed bandwidth for the scheduling of the FL application, without
jeopardizing the QoS provisioning of other delay-critical applications. To
achieve this, the proposed mechanism adopts the widely-used Differentiated
Service approach to tackle the QoS provisioning problem of Federated Learning
applications over Ethernet PON. The prioritization of FL traffic is used to
comply with the traditional business model, as well as to improve statistical
multiplexing gain.

The proposal employs an intra-ONU scheduler with a strict priority queuing
policy to serve the ONU queues. The ONUs arbitrates the transmission demands
of the different applications. Upon the arrival of a Report message, the OLT
calculates the transmission window according to the conventional limited
policy and selects the wavelength(s) on the basis of the TWDM wavelength
allocation policies. The OLT then sends a Gate message containing the resource
allocation decision. Upon the receipt of that Gate message, the intra-ONU
scheduler distributes the received bandwidth among the queues in the ONU. In
our model, traffic is classified as Federated Learning, delay-critical, delay-
sensitive, or best-effort. The ONUs mantains four different queues for
buffering frames for these types of traffic.

We propose two prioritization policies. The FL-first policy defines the FL
traffic as that of the $highest$ priority and the delay-critical, delay-
sensitive, and best-effort being of $high$, $medium$ and $low$ priority,
respectively. On the one hand, this strict prioritization of FL frames can
reduce synchronization time for FL processing. It can also increase the delay
of delay-critical application because the FL traffic requires a large
bandwidth per cycle. To help alleviate this problem, we propose the DC-first
policy, in which delay-critical traffic also has the highest priority and
Federated Learning traffic the $high$ priority. Moreover, the proposal was
defined for all TWDM architectures.

##  V Performance Evaluation

The performance of the proposed DWBA scheme was evaluated using an EPON
simulator (EPON-Sim), previously validated in [14]. This extension was
extended to support the three architectures, SSD, MSD and WA, proposed for
50G-EPON networks. Moreover, our proposal and the BS approach were introduced
in the simulator.

###  V-A Simulation Model and Setup

The simulation scenarios include a 50G-EPON network with $1$ OLT serving 32
ONUs on an optical distribution network with a tree topology. Two wavelength
channels of 25 Gbit/s were employed for upstream transmission, giving a total
capacity of 50 Gbit/s. The total available bandwidth in the PON was equally
distributed among the ONUs, so that each ONU has the same guaranteed bandwidth
$b_{i}$, while the aggregated offered load per ONU $l_{i}$ varied from
$0.6\cdot b_{i}$ to $1.0\cdot b_{i}$ (for the sake of clearness and brevity,
herein after, $b_{i}$ is omitted from the offered load values of ONU $i$).

The aggregated load included the traffic of the four different types of
application: Federated Learning, delay-critical, delay-sensitive, and Best
Effort. The benchmarking framework for learning in federated settings LEAF [5]
was used to generate the FL traffic. The FEMNIST dataset and CNN with two 5×5
convolution layers were used for model training, while the FedAvg algorithm
was employed to aggregate the local parameters in the server. Other
configurations for the learning process, such as learning rate and batch size,
followed the settings defined in [15]. FL clients generated 26.4 MBytes of
data in each round of training. Moreover, the ONUs put the local parameters
into frames according to the Ethernet protocol, which has a Maximum
Transmission Unit of 1500 bytes and a header field for signaling (preamble) of
20 bytes.

The delay-critical applications were modeled employing a Constant Bit Rate
(CBR) flow. It was coded with a fixed-size packet of 70 bytes and an inter-
arrival time of 12.5 $\mu$s, which produces an offered load of 44.8 Mbps. The
rest of the offered load $l_{i}$ was evenly distributed between delay-
sensitive and Best Effort traffic. The traffic streams were generated
employing Pareto ON-OFF sources. The ON period time and packet-burst size
followed a Pareto and Bounded Pareto distributions, respectively. The
aggregated traffic at the ONU had a Hurst parameter of 0.8. Moreover, the
packet lengths are uniformly distributed between $64$ and $1518$ bytes.

A threshold value of $\theta=0.015$ was employed in the MW-BS algorithm, as in
[4]. This algorithm reduces the bandwidth for each ONU since it reserves
bandwidth for the slice. Moreover, we employ the same aggregated offered load
in the simulated algorithms to make a fair comparison. The duration of the
guard period was set to 0.624 $\mu s$ with a maximum cycle length of 1 ms.
Each simulation scenario lasted 100 s and was replicated $10$ times.

(a) Delay-critical application

(b) Federated Learning application

Figure 4: Average delay produced by the delay-critical and FL applications.
$P1$ and $P2$ mean, respectively, DWBA-FL algorithm with FL-first policy and
DC-first policy.

(a) Involved clients vs. synchronization time

(b) Learning accuracy (with 2000 rounds) vs. synchronization time

Figure 5: Percentage of the learning accuracy and involved clients depending
on the synchronization time. $P2$ means DWBA-FL algorithm with DC-first
policy. Figure 6:  Network delay performance for the Federated Learning
application. Each boxplot shows 10th, 30th, 50th, 80th, and 100th percentile
of network delay for FL clients; $P2$ means DWBA-FL algorithm with DC-first
policy..

###  V-B Simulation Results and Discussion

Mean delay values obtained by DWBA-FL were lower than 80 ms and 150 ms for the
FL traffic in both underloaded and overloaded conditions, respectively. The
delay values given by BS were at least twice as high as those given by our
proposal (Fig. 4b). This improvement is a consequence of the large windows
allocated for transmissions of FL traffic when our proposal is employed.

Moreover, the use of the DC-first policy produced lower delay values for the
delay-critical traffic lower than those given by both the FL-first policy and
the BS algorithm (Fig. 4a). This result occurs because the bandwidth slice is
statically allocated for the FL traffic. Furthermore, the strict
prioritization of FL traffic employing the FL-first policy and the huge amount
of traffic produced by the FL application leads to delay-critical application
to bandwidth starvation. The mean delay of the delay-critical traffic produced
by the FL-first policy was from 200 us to 1000 us and is more than that
produced by FL-first policy. Thus, the DWBA-FL with DC-first policy produce
mean delay values for the Federated Learning and delay-critical application
lower than the other algorithms.

Furthermore, the FF policy produces a slight decrease in delay values for both
traffics than the other wavelength allocation policies (i.e.,  SSD and MSD).
This results are a consequence of the wastage of bandwidth due to the
excessive uses of guard periods and poor multiplexing gain when employing the
SSD and MSD, respectively.

In Fig. 5a, the blue curve shows the proportion of clients involved as a
function of the computing time. It shows the minimal synchronization time per
round without any communication delay. The MW-BS algorithm requires a longer
synchronization time per round to produce the same percentage of the involved
clients than that required by the proposed scheme with DC-first policy. For
example, it is required a synchronization time of 1.9 s and 2.1 s to produce a
percentage of involved clients of 50 % with our proposal and the MW-BS
algorithm, respectively. To achieve a training accuracy of 76% (Fig. 5b), the
proposed scheme can reduce 9.5% of the training time compared to the BS
algorithm (i.e.,  0.2 s less for a synchronization time of 2.1 s), when the
total traffic load is 0.8.

Fig. 6 shows the network delay as a function of the ONU offered load. The MW-
BS produce delay values greater than 300 ms, whereas, with the DWBA-FA
algorithm, these values are reduced to less than 150 ms. Moreover, for 80 % of
the clients, which is the typical percentage of clients that produce accuracy
greater than 75 % (see Fig. 5), the MW-BS scheme imposes a network delay
greater than 200 ms, while the DWBA-FA imposes delay values lower than 100 ms,
under underloaded conditions (i.e.,  load ¡ 0.85). In summary, DWBA-FL reduces
the network delay when compared to the MW-BS scheme. This reduction in delay
may decrease the number of stragglers, which in the end leads to a faster
convergence and greater model accuracy.

##  VI Conclusion

This paper has introduced a resource allocation (RA) scheme for supporting
Federated Learning applications in TWDM-EPONs networks. Our proposal includes
a strict prioritization for federated learning and delay-critical traffic,
which maximizes the allocated bandwidth and reduce the delay for both types of
application. Our proposal reduces the synchronization time without
compromising the number of involved clients. It also reduces the delay of
time-critical and federated learning applications. Future research directions
are envisioned as follows. Various FL applications coexisting on a given
network infrastructure, which have different size of the local model
parameters, number of clients and synchronization time. Mechanisms are needed
to appropriately address the QoS provisioning for the diversity of FL
applications. These schemes may schedule the FL traffics based on required
bandwidth but also considering the number of straggler clients, the diverse FL
packet sizes, and synchronization time.

## Acknowledgments

This work was partially sponsored by grant #15/24494-8, São Paulo Research
Foundation (FAPESP), and CNPq.

## References

  * [1] K. B. Letaief, W. Chen, Y. Shi, J. Zhang, and Y.-J. A. Zhang, “The roadmap to 6g: Ai empowered wireless networks,” _IEEE Communications Magazine_ , vol. 57, no. 8, pp. 84–90, 2019.
  * [2] B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y Arcas, “Communication-efficient learning of deep networks from decentralized data,” in _Artificial Intelligence and Statistics_. PMLR, 2017, pp. 1273–1282.
  * [3] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 770–778.
  * [4] J. Li, X. Shen, L. Chen, and J. Chen, “Bandwidth slicing to boost federated learning over passive optical networks,” _IEEE Communications Letters_ , vol. 24, no. 7, pp. 1492–1495, 2020.
  * [5] S. Caldas, S. M. K. Duddu, P. Wu, T. Li, J. Konečnỳ, H. B. McMahan, V. Smith, and A. Talwalkar, “Leaf: A benchmark for federated settings,” _arXiv preprint arXiv:1812.01097_ , 2018.
  * [6] F. de Melo Pereira, N. L. S. da Fonseca, and D. S. Arantes, “A fair scheduling discipline for ethernet passive optical networks,” _Comput. Networks_ , vol. 53, no. 11, pp. 1859–1878, 2009.
  * [7] J. S. Wey, “The outlook for pon standardization: a tutorial,” _Journal of Lightwave Technology_ , vol. 38, no. 1, pp. 31–42, 2020.
  * [8] J. Li, L. Chen, and J. Chen, “Scalable federated learning over passive optical networks,” _arXiv preprint arXiv:2010.15454_ , 2020.
  * [9] “Ieee standard for ethernet amendment 9: Physical layer specifications and management parameters for 25 gb/s and 50 gb/s passive optical networks,” _IEEE Std 802.3ca-2020 (Amendment to IEEE Std 802.3-2018 as amended by IEEE 802.3cb-2018, IEEE 802.3bt-2018, IEEE 802.3cd-2018, IEEE 802.3cn-2019, IEEE 802.3cg-2019, IEEE 802.3cq-2020, IEEE 802.3cm-2020, and IEEE 802.3ch-2020)_ , pp. 1–267, 2020.
  * [10] M. P. McGarry, M. Reisslein, and M. Maier, “Wdm ethernet passive optical networks,” _IEEE Communications Magazine_ , vol. 44, no. 2, pp. 15–22, 2006\.
  * [11] L. Wang, X. Wang, M. Tornatore, H. S. Chung, H. H. Lee, S. Park, and B. Mukherjee, “Dynamic bandwidth and wavelength allocation scheme for next-generation wavelength-agile epon,” _Journal of Optical Communications and Networking_ , vol. 9, no. 3, pp. B33–B42, 2017.
  * [12] S. B. Hussain, W. Hu, H. Xin, and A. M. Mikaeil, “Low-latency dynamic wavelength and bandwidth allocation algorithm for ng-epon,” _Journal of Optical Communications and Networking_ , vol. 9, no. 12, pp. 1108–1115, 2017.
  * [13] M. P. McGarry, M. Reisslein, and M. Maier, “Ethernet passive optical network architectures and dynamic bandwidth allocation algorithms,” _IEEE Communications Surveys & Tutorials_, vol. 10, no. 3, pp. 46–60, 2008.
  * [14] O. J. Ciceri, C. A. Astudillo, G. B. Figueiredo, Z. Zhu, and N. L. S. da Fonseca, “Passive optical networking for 5g and beyond 5g low-latency mobile fronthauling services,” 2021.
  * [15] J. Li, L. Chen, and J. Chen, “Scalable federated learning over passive optical networks,” _arXiv preprint arXiv:2010.15454_ , 2020.

Oscar Ciceri<EMAIL_ADDRESS>completed his five-year degree in
Electronics and Telecommunications Engineering at Universidad del Cauca
(UNICAUCA), Colombia, in 2015, and his M.S. degree in Computer Science at the
University of Campinas (UNICAMP), Brazil, in 2019. He is a Ph.D. student and
researcher in the Network Computing Laboratory (LRC) at that university. His
research interests include Passive Optical Networks (PONs), 5G and 5G beyond
networks, Machine Learning, and virtualization. His current work is on the
low-latency PON technologies for 5G and 5G beyond mobile optical access
systems.
---
Carlos Astudillo<EMAIL_ADDRESS>received his bachelor’s degree
in Electronics and Telecommunications Engineering from the University of
Cauca, Colombia in 2009, and his MSc degree in Computer Science from the
University of Campinas (UNICAMP), Brazil, in 2015. Currently, he is a PhD
candidate in Computer Science at the Institute of Computing (IC), UNICAMP. His
master’s thesis received 1st place at the 29th Theses and Dissertations
Contest (CTD) of the Brazilian Computer Society (SBC), and the 2016 Best
Master’s Thesis Award from IC/UNICAMP. His research interests include radio
resource allocation and quality of service provisioning for machine to machine
communications and the Internet of Things, backhauling/fronthauling in 5G &
B5G cellular networks, and machine learning for networking.
---
Zuqing Zhu [M’07, SM’12]<EMAIL_ADDRESS>received his Ph.D. degree from the
Department of Electrical and Computer Engineering, University of California,
Davis, in 2007. He is currently a full professor at USTC. He has published
more than 200 papers in referreed journals and conferences. He was an IEEE
Communications Society Distinguished Lecturer (2018–2019) and a Senior Member
of OSA.
---
Nelson L. S. da Fonseca [M’88, SM’01]<EMAIL_ADDRESS>obtained his
Ph.D. degree from the University of Southern California in 1994. He is Full
Professor at the Institute of Computing, State University of Campinas
(UNICAMP). He published 400+ papers and supervised 70+ graduate students. He
is Senior Editor for the IEEE Communications Magazine and the IEEE Systems
Journal. He is also an Associate Editor for Peer-to-Peer Networking and
Applications and Computer Networks. He is past EiC of the IEEE Communications
Surveys and Tutorials. He served as the IEEE ComSoc VP Publications, VP
Technical and Educational Activities, and VP Member Relations.
---
  *[AI]: Artificial Intelligence
  *[6G]: 6th Generation mobile networks
  *[ML]: Machine Learning
  *[FL]: Federated Learning
  *[PON]: Passive Optical Network
  *[TDM]: Time Division Multiplexing
  *[PONs]: Passive Optical Network
  *[TWDM]: Time and Wavelength Division Multiplexing
  *[MFH]: Mobile Fronthauling
  *[BS]: Bandwidth Slicing
  *[ONU]: Optical Network Unit
  *[OLT]: Optical Line Terminator
  *[DWBA]: Dynamic Wavelength and Bandwidth Allocation
  *[EPON]: Ethernet PON
  *[]:
  *[GPON]: Gigabit Capable
  *[IEEE]: Institute of Electrical and Electronics Engineers
  *[InP]: Infrastructure service Provider
  *[TWDMA]: Time and Wavelength Division Multiple Access
  *[ONUs]: Optical Network Unit
  *[MSD]: Multiple-Scheduling Domain
  *[SSD]: Single-Scheduling Domain
  *[WA]: Wavelength Agile
  *[MPCP]: Multipoint Control Protocol
  *[DBA]: Dynamic Bandwidth Allocation
  *[IPACT]: Interleaved Polling with Adaptive Cycle Time
  *[WF]: Water-Fill
  *[FF]: First-Fit
  *[QoS]: Quality of Service
  *[CNN]: Convolutional Neural Network
  *[CBR]: Constant Bit Rate
  *[EPONs]: Ethernet PON


# Weakly coupled systems of eikonal equations for autonomous vehicles

Maria Teresa Chiri1, Kenneth D Czuprynski2 and Ludmil T Zikatanov3 *The
research of Ludmil Zikatanov is based upon work supported by and while serving
at the National Science Foundation. Any opinion, findings, and conclusions or
recommendations expressed in this material are those of the author and do not
necessarily reflect the views of the National Science Foundation.1 Department
of Mathematics and Statistics, Queen’s University, Kingston, ON, Canada
<EMAIL_ADDRESS>Applied Research Laboratory, The Pennsylvania State
University, University Park, PA, USA<EMAIL_ADDRESS>National Science
Foundation, Alexandria, VA, USA<EMAIL_ADDRESS>
###### Abstract

In this paper, we study solutions for a weakly coupled system of eikonal
equations arising in an optimal path-planning problem for autonomous vehicles.
The model considered takes into account the possibility of two types of
breakdown, partial and total, which happen at a known, spatially inhomogeneous
rate. In particular, we show how to bypass the delicate degenerate coupling
condition by using existing results on weakly coupled systems of Hamilton
Jacobi equations. Then we consider finite element method schemes built for
convection-diffusion problems to construct approximate solutions for this
system and produce some numerical simulations.

##  I Introduction

In the era of autonomous vehicles, there is an increasing need to develop
robust and widely applicable path-planning algorithms. The models introduced
so far are meant for a variety of autonomous vehicles which range from
planetary exploration rovers [1] to flying drones [2], underwater vehicles
[3], and of course cars [4]. Although the literature on path planning methods
used in robotics is already vast (see [5, 6, 7, 8, 9] for some examples) and
the number of models continues to increase over the years, there are still
important issues that need to be addressed. A new challenge, for instance, is
to derive models taking into account accidental events such as vehicle
breakdowns, even in the case of a completely known environment. In this
respect, in [10] the authors consider an optimal path model in which a vehicle
switches between two different modes, representing a partial and a total
breakdown condition. The switch happens stochastically at known rates and each
mode has its own deterministic dynamics and running cost. Their approach is
based on dynamic programming in continuous state and time and leads to
globally optimal trajectories by solving a system of two weakly coupled
eikonal equations. This represents a simple problem in the context of the
general family of weakly coupled systems of Hamilton-Jacobi equations which
appear frequently in the literature of optimal control problems with random
switching costs governed by Markov chains [11, 12, 13]. Recently weakly
coupled Hamilton-Jacobi systems were studied from a viewpoint of weak KAM
theory. For example, [14, 15] and [16] investigated the large-time behavior of
the solution for time-dependent problems. In [17], the authors studied
homogenization for weakly coupled systems and the rate of convergence to
matched solutions while [18] generalized the notion of Aubry sets for the case
of systems and proved comparison principle with respect to their boundary data
on Aubry sets.
We will study the system derived in [10] in the framework of viscosity
solutions, as previously done in [19, 20]. In particular, we will prove a
comparison principle in the classical sense [21] with the difference that the
boundary conditions will be prescribed not only along the perimeter of the
domain but also on the Aubry set of the system. On the other side, the
existence of solutions is a very delicate issue. We will discuss an example of
a one-dimensional weakly coupled system of eikonal equations showing how the
lack of boundary conditions on the Aubry set affects the uniqueness of
solutions and we will provide boundary conditions for which the system does
not admit solutions.
Numerically, a range of different methods have been applied to solve the
eikonal equations. Fast marching methods [22] and fast sweeping methods [23]
are efficient algorithms for the causal propagation of boundary information
into the domain. In [10] the coupled system of eikonal equations is solved
using FMMs within a value iteration scheme which uses policy evaluations to
speed up convergence. A number of finite element formulations have also been
considered for the eikonal equations (cf. [24, 25, 26]). In [25], an
artificial viscosity regularization with a homotopy scheme for the stable
reduction of the viscosity is introduced. High-order discontinuous Galerkin
formulations have also been developed for the eikonal equations [27] and are
amenable to parallelization. In this work, we consider an artificial viscosity
formulation similar to [25] but consider finite element schemes built for
convection-diffusion problems [28] as means of achieving stable numerical
solutions. We point out that numerical simulations will be performed in a case
in which the Aubry set is empty to avoid the theoretical complications
described above.

The paper is organized as follows: in Section II we will describe the model
introduced in [10], in Section III we will prove a comparison principle for
viscosity solutions of the weakly coupled system of eikonal equations and show
an example of non-uniqueness and non-existence; in Section IV we will derive a
numerical scheme based on a finite element method to approximate solutions for
this system, while simulations will be discussed in section V. Finally in
Section VI we will describe possible future research directions.

##  II The model

Inspired by [10], we consider the problem of optimal path planning for an
autonomous vehicle (AV) subject to two potential breakdown events:
1\. Partial breakdown or Mode 1, in which the vehicle is damaged but can keep
moving toward its target;
2\. Total breakdown or Mode 2, in which the vehicle is damaged and needs to
move toward a repair depot.
The classical formulation of the optimal path planning problem consists of an
agent that tries to optimally travel from one point to another while obeying
an ordinary differential equation describing its motion. In our analysis, we
will consider AVs moving on a stretch of road modeled by a rectangle
$\Omega=[0,L]\times[0,S]$ where $L$ represents the length and $S$ is the
width. The vehicles’ dynamic is described by the ODE

$\dot{y}(t)=d(t)f(y(t)),\quad y(0)=x_{0}$

where $d:[0,T]\to S^{1}$ is the direction of motion, $f$ is the speed and
$y(t)$ is the trajectory starting at $x_{0}\in\Omega$. Each breakdown type is
assumed to occur at some known rate and accrues some predetermined cost.
Denote with $T_{b}$ the time at which the first partial breakdown occurs,
$T_{G}$ the time at which the AV would reach the target if no breakdown
occurred, and with $T_{1}=\min\\{T_{b},T_{G}\\}$, then the expected cost in
Mode 1 is given by

$\displaystyle J_{1}(x,a(\cdot))=\mathbb{E}\Big{[}$
$\displaystyle\int_{0}^{T_{1}}K_{1}(y(s))+\lambda(y(s))R(y(s))ds$
$\displaystyle+\chi_{[T_{b}<T_{G}]}u_{2}(y(T_{b}))\Big{]}\,,$ (1)

where $u_{2}$ is the value function in Mode 2, $K_{1}$ is the running cost,
$\lambda>0$ is the rate at which the partial breakdown occurs (trajectory
dependent) and $R$ is the repair cost. If a partial breakdown occurs before
the robot reaches the goal, then the last term in (II) is nonzero and the AV
switches to Mode 2. Similarly, let $T_{B}$ be the time of the first total
breakdown, $T_{D}$ the time the AV would reach the depot if no breakdown
occurred, and denote with $T_{2}=\min\\{T_{B},T_{D}\\}$, then the expected
cost in Mode 2 is given by

$\displaystyle J_{2}(x,a(\cdot))=\mathbb{E}\Big{[}$
$\displaystyle\int_{0}^{T_{2}}K_{2}(y(s))ds$
$\displaystyle+\chi_{[T_{B}<T_{D}]}\Big{(}R(y(T_{B}))+u_{1}(y(T_{B}))\Big{)}$
$\displaystyle+\chi_{[T_{B}\geq T_{D}]}R(y(T_{D}))+u_{1}(y(T_{D}))\Big{]}\,,$
(2)

where $K_{2}$ is the running cost. The last two terms in (II) again encode the
optimal cost-to-go whenever a mode switch occurs. Let $G\subset\Omega$ denote
the goal locations for admissible AV trajectories and let $D\subset\Omega$
denote the location of the repair depots that the AV can navigate to after
undergoing a partial breakdown. According to the derivation described in [29],
the value functions $u_{1}$ and $u_{2}$ solve the following system of eikonal
equations:

$\begin{cases}|\nabla u_{1}|f_{1}+\phi_{1}(u_{1}-u_{2})=K_{1}+\lambda
R,\quad&x\in\Omega\\\ |\nabla
u_{2}|f_{2}+\phi_{2}(u_{2}-u_{1})=K_{2}+\phi_{2}R,\quad&x\in\Omega\\\
u_{1}=0,\quad&x\in G\\\ u_{2}=R+u_{1},\quad&x\in D\end{cases}$ (3)

Observe that $i$-th equation in the system depends only on $\nabla u_{i}$, but
not on $\nabla u_{j}$ for $i\neq j$, therefore we refer to $\eqref{eikonal-
system-original}$ as a weakly coupled system. The functions
$\lambda,\phi_{2}:\Omega\rightarrow(0,\infty)$ correspond to the rates of
exponential random variables modeling the time until another complete
breakdown and they are trajectory dependent; the function
$\phi_{1}(x):\Omega\rightarrow(0,\infty)$ denotes the rate at which partial
breakdowns occur. The asymmetry in the notation for the rates of the
exponential random variables and the rate for the partial breakdowns comes
from the necessity of giving a consistent notation to the coupling matrix. The
functions $K_{i}:\Omega\rightarrow[0,+\infty)$ denote the running cost
associated with the cost functional, while $f_{i}:\Omega\to[0,+\infty)$ are
the space dependent speeds.

##  III Viscosity solutions and comparison principle

In this section, we want to prove a comparison principle for viscosity
solutions of the weakly coupled system (3). We start by recalling the
definition of viscosity solution introduced in [30], in the context of the
system we are considering.

###### Definition 1

A continuous function $u:\Omega\to\mathbb{R}^{2}$, with
$u(x)=(u_{1}(x),u_{2}(x))$, is a viscosity subsolution (resp. supersolution)
of (3), if for each $i=1,2$ and test function $\varphi$ of class
$\mathcal{C}^{1}$, when $u_{i}-\varphi$ attains its local maximum (resp. local
minimum) at $x$, then

$\begin{aligned} |\nabla\varphi|f_{1}+\phi_{1}u_{1}&\leq K_{1}+\lambda
R+\phi_{1}u_{2},\\\ |\nabla\varphi|f_{2}+\phi_{2}u_{2}&\leq
K_{2}+\phi_{2}R+\phi_{2}u_{1},\end{aligned}\quad(\hbox{resp.}~{}\geq~{})$ (4)

A function is called a viscosity solution of (3) if it is both a viscosity
subsolution and supersolution.

A standard assumption in literature for weakly coupled systems of Hamilton-
Jacobi equations is the so-called monotonicity condition, i.e. the $i$-th
Hamiltonian is increasing in $u_{i}$ and decreasing in $u_{j}$ with $j\neq i$.
Such a condition is clearly satisfied by (3). Moreover, the Hamiltonians are
locally Lipschitz continuous, coercive, and strictly convex in the gradient of
$u$. On the other side, the coupling matrix for the system is given by

$D(x)=\begin{pmatrix}\phi_{1}&-\phi_{1}\\\ -\phi_{2}&\phi_{2}\end{pmatrix}$
(5)

which is degenerate (the sum of the element of each row is equal to zero) and
irreducible for every $x\in\Omega$, this simply means that coupling is not
trivial. Degenerate couplings have been avoided for a long time, since they
represented an obstruction to achieving uniqueness of the solution (cf. [19]).
This was until it was realized this type of coupling induces pathologies on
the Aubry set which behaves as a boundary in the inner part of the domain [31,
14] . This set, in the particular case of (3), is given by

$\mathcal{A}\doteq\\{x\in\Omega,K_{1}(x)+K_{2}(x)+(\lambda(x)+\phi_{2}(x))R(x)=0\\}$

(since the right-hand side of both the equations is non-negative) and we may
have $\mathcal{A}\setminus\partial\Omega\neq\emptyset$. Therefore to prove a
comparison principle and consequently to get the uniqueness of the solution,
it is essential to add boundary conditions on $\mathcal{A}\cup\partial\Omega$.
In the following, we will use the notation $USC(\Omega)$ and $LSC(\Omega)$ for
the set of upper semicontinuous and lower semicontinuous functions on
$\Omega$.

###### Theorem 1

[Comparison Principle] Let $\phi_{1},\phi_{2},\lambda,R:\Omega\to\mathbb{R}$
be continuous functions. Let $u\in USC(\Omega)$ and $v\in LSC(\Omega)$ be
respectively bounded subsolutions and supersolution of (3). Assume that
$u(x)\leq v(x)$ for every $x\in\mathcal{A}\cup\partial\Omega$, then $u(x)\leq
v(x)$ on $\bar{\Omega}$

###### Proof:

Consider $\delta\in(0,1)$ an set

$M_{\delta}\doteq\sup_{[0,L]\times[0,S]}\sup_{i=1,2}\\{\delta u_{i}-v_{i}\\}$

If $M_{\delta}\leq 0$ there is nothing to prove, therefore we assume it is
strictly positive. By compactness the $sup$ is a maximum achieved at a point
$x_{0}\in[0,L]\times[0,S]$ for some $i_{0}$. Denote with $\mathcal{I}$ the set
of the index for which the maximum is achieved. We can distinguish three
cases.
Case 1. $\mathcal{I}=\\{1,2\\}$ and $x_{0}\in\mathcal{A}\cup\partial\Omega$.
In this case we have

$\displaystyle M_{\delta}$ $\displaystyle=\frac{1}{2}\big{[}(\delta
u_{1}-v_{1})(x_{0})+(\delta u_{2}-v_{2})(x_{0})\big{]}$
$\displaystyle\leq\frac{1}{2}\leq(\delta-1)\big{(}v_{1}(x_{0})+v_{2}(x_{0})\big{)}$
$\displaystyle\leq(1-\delta)|v|_{\infty}.$

Case 2. $\mathcal{I}=\\{1,2\\}$ and
$x_{0}\notin\mathcal{A}\cup\partial\Omega$. For instance, we can assume
without loss of generality that $K_{1}(x_{0})+\lambda(x_{0})R(x_{0})\neq 0$.
In this case, we proceed classically by considering

$\sup_{\Omega\times\Omega}\Big{\\{}\delta
u_{1}(x)-v_{1}(y)-\frac{|x-y|^{2}}{2\varepsilon^{2}}-|x-x_{0}|^{2}\Big{\\}}.$

Let $(\bar{x},\bar{y})$ be the point at which this maximum is achieved. This
value is clearly greater than $M_{\delta}$, moreover, the following properties
are satisfied:

$\bar{x},\,\bar{y}\to
x_{0}\,\qquad\frac{|x-y|^{2}}{2\varepsilon^{2}},\,|\bar{x}-x_{0}|^{2}\to
0\quad\hbox{ as }\varepsilon\to 0.$

Let $p_{\varepsilon}=\frac{\bar{x}-\bar{y}}{\varepsilon^{2}}$. Since $u_{1}$
is viscosity subsolution of (3), we get

$\displaystyle\Big{|}{p_{\varepsilon}+2(\bar{x}-x_{0})}\Big{|}f_{1}(\bar{x})+\delta\phi_{1}(\bar{x})$
$\displaystyle(u_{1}(\bar{x})-u_{2}(\bar{x}))$
$\displaystyle\leq\delta\big{(}K_{1}(\bar{x})+\lambda(\bar{x})R(\bar{x})\big{)}.$

On the other side, since $v_{1}$ is viscosity supersolution

$|p_{\varepsilon}|f_{1}(\bar{y})+\phi_{1}(\bar{y})(v_{1}(\bar{y})-v_{2}(\bar{y}))\geq\big{(}K_{1}(\bar{y})+\lambda(\bar{y})R(\bar{y})\big{)}.$

By subtracting the two inequalities we first get

$\displaystyle\Big{|}{p_{\varepsilon}+2(\bar{x}-x_{0})}\Big{|}f_{1}(\bar{x})-|p_{\varepsilon}|f_{1}(\bar{y})$
$\displaystyle=\Big{(}\Big{|}{p_{\varepsilon}+2(\bar{x}-x_{0})}\Big{|}-|p_{\varepsilon}|\Big{)}f_{1}(\bar{x})+|p_{\varepsilon}|\big{(}f_{1}(\bar{x})-f_{1}(\bar{y})\big{)}$

which goes to $0$ as $\varepsilon\to 0$, and

$\displaystyle\delta\phi_{1}(\bar{x})$
$\displaystyle\big{(}u_{1}(\bar{x})-u_{2}(\bar{x})\big{)}-\phi_{1}(\bar{y})\big{(}v_{1}(\bar{y})-v_{2}(\bar{y})\big{)}$
$\displaystyle=\phi_{1}(\bar{x})\big{(}\delta u_{1}(\bar{x})-\delta
u_{2}(\bar{x})-v_{1}(\bar{y})+v_{2}(\bar{y})\big{)}$
$\displaystyle\qquad\qquad+(\phi_{1}(\bar{x})-\phi_{1}(\bar{y}))(v_{1}(\bar{y})-v_{2}(\bar{y}))$

where

$(\phi_{1}(\bar{x})-\phi_{1}(\bar{y}))(v_{1}(\bar{y})-v_{2}(\bar{y}))\to
0\qquad\hbox{as }\quad\varepsilon\to 0.$

Therefore we get

$\displaystyle\phi_{1}(\bar{x})$ $\displaystyle\big{(}\delta
u_{1}(\bar{x})-\delta u_{2}(\bar{x})-v_{1}(\bar{y})+v_{2}(\bar{y})\big{)}$
$\displaystyle\leq(\delta-1)\big{(}K_{1}(\bar{x})+\lambda(\bar{x})R(\bar{x})\big{)}+o_{\varepsilon}(1).$

Since $\delta u_{j}-v_{j}$ is upper-semicontinuous for $j=1,2$, it follows
that

$\limsup_{\varepsilon\to 0}(\delta u_{j}(\bar{x})-v_{j}(\bar{y}))\leq\delta
u_{j}(x_{0})-v_{j}(x_{0})\leq M_{\delta}.$ (6)

hence

$-\phi_{1}(\bar{x})(\delta
u_{2}(\bar{x})-v_{2}(\bar{y}))\geq-\phi_{1}(\bar{x})M_{\delta}+o_{\varepsilon}(1).$

Moreover we have $\delta u_{1}(\bar{x})-v_{1}(\bar{y})\geq M_{\delta}$,
therefore

$\phi_{1}(\bar{x})(\delta
u_{1}(\bar{x})-v_{1}(\bar{y}))\geq\phi_{1}(\bar{x})M_{\delta}.$

We can conclude that

$(\delta-1)\big{(}K_{1}(\bar{x})+\lambda(\bar{x})R(\bar{x})\big{)}\geq 0$

which leads to contraddiction when $\varepsilon\to 0$ since
$x_{0}\notin\mathcal{A}$.

Case 3. $\mathcal{I}\neq\\{1,2\\}$. We can argue almost like in Case 2, but
the estimate needs to be more rigorous on the index not in $\mathcal{I}$.
Indeed let’s assume that $\mathcal{I}=\\{{1\\}}$. The previous computations
still hold, but for (6), if $j=2$

$\limsup_{\varepsilon\to 0}(\delta u_{2}(\bar{x})-v_{2}(\bar{y}))\leq\delta
u_{2}(x_{0})-v_{2}(x_{0})\leq M_{\delta}-\eta$

for some $\eta>0$. It follows that

$\phi_{1}(x_{0})\eta\leq(\delta-1)\big{(}K_{1}(\bar{x})+\lambda(\bar{x})R(\bar{x})\big{)}+o_{\varepsilon}(1)$

which again leads to contradiction as $\varepsilon\to 0$.

We can conclude that the only possibility is
$M_{\delta}\leq(1-\delta)|v|_{\infty}$ which in the case of $\delta=1$ implies
$M_{1}\leq 0$. ∎

###### Lemma 1

Let $u\in USC(\Omega)$ be a bounded viscosity subsolution of (3), then $u$ is
Lipschitz continuous on $\Omega$.

###### Proof:

Since $\Omega$ is a compact set, we have that
$K_{1},K_{2},\phi_{1},\phi_{2},\lambda$ are bounded. Therefore

$|\nabla u_{i}|f_{i}(x)\leq C_{i}\qquad\hbox{i=1,2}$

which implies that $u$ is Lipschitz. ∎

The comparison principle proved in Theorem 1 implies that if a solution
exists, it must be unique. On the other hand, it is not clear which type of
boundary conditions on $\partial\Omega$ and $\mathcal{A}$ ensure the existence
of viscosity solutions. We borrow an example from [14] to illustrate how the
set $\mathcal{A}$ can affect the uniqueness of solutions for a weekly coupled
system. Moreover, we provide bundary conditions on $\mathcal{A}$ for which
such a system has no solution.

###### Example 1

Consider the one-dimensional problem

$\begin{cases}|u^{\prime}_{1}(x)|+u_{1}(x)-u_{2}(x)=F(x),\quad&x\in(-1,1)\\\
|u^{\prime}_{2}(x)|+u_{2}(x)-u_{1}(x)=F(x)\quad&x\in(-1,1)\\\ u_{i}(\pm
1)=0\quad&i=1,2\end{cases}$ (7)

with $F(x)=2|x|$. It is immediate to check that both
$u_{1}(x)=u_{2}(x)=1-x^{2}$ and $u_{1}(x)=u_{2}(x)=\min\\{1-x^{2},x^{2}-C\\}$,
$C\in(0,1)$, are viscosity solutions of (7). Here the set $\mathcal{A}$ is
given by $\\{0\\}$, hence without prescribing a boundary condition in $0$ we
do not have uniqueness. If we just add the condition $u_{i}(0)=\frac{1}{2}$,
this immediately selects the unique solution
$u_{1}(x)=u_{2}(x)=\min\\{1-x^{2},x^{2}-\frac{1}{2}\\}$. On the other side for
$u_{i}(0)=2$, no solution satisfies (7).

##  IV Numerical Solution

In this section the numerical solution of the weakly coupled system (3) using
the finite element method (FEM) is considered. The numerical scheme will solve
a sequence of systems augmented by artificial viscosity terms. Formally, we
have the problem: Find $u^{\varepsilon}_{1}$ and $u^{\varepsilon}_{2}$ such
that

$\displaystyle-\varepsilon\Delta u_{1}^{\varepsilon}+|\nabla
u^{\varepsilon}_{1}|f_{1}+\phi_{1}u^{\varepsilon}_{1}$
$\displaystyle=K_{1}+\lambda R+\phi_{1}u^{\varepsilon}_{2},$ (8)
$\displaystyle-\varepsilon\Delta u^{\varepsilon}_{2}+|\nabla
u^{\varepsilon}_{2}|f_{2}+\phi_{2}u^{\varepsilon}_{2}$
$\displaystyle=K_{2}+\phi_{2}R+\phi_{2}u^{\varepsilon}_{1},$

for $x\in\Omega$ with $u^{\varepsilon}_{1}=0$ for $x\in G$ and
$u^{\varepsilon}_{2}=R+u^{\varepsilon}_{1}$ for $x\in D$, as
$\varepsilon\rightarrow 0$. Although the viscosity term provides some
mechanism for control of the solution, we are still left to deal with the
original nonlinearity which is accomplished through linearization as in [25].
Let $u^{\varepsilon_{n}}_{k}$ denote the solution of (8) for the $n^{th}$
value of $\varepsilon$, where $k=1,2$. Linearization of the absolute value
term results in the approximation

$\displaystyle|\nabla u^{\varepsilon_{n}}_{k}|$ $\displaystyle\approx|\nabla
u^{\varepsilon_{n-1}}_{k}|+\frac{\nabla u^{\varepsilon_{n-1}}_{k}}{|\nabla
u^{\varepsilon_{n-1}}_{k}|}\left(\nabla u^{\varepsilon_{n}}_{k}-\nabla
u^{\varepsilon_{n-1}}_{k}\right)$ $\displaystyle\approx\frac{\nabla
u^{\varepsilon_{n-1}}_{k}}{|\nabla u^{\varepsilon_{n-1}}_{k}|}\nabla
u^{\varepsilon_{n}}_{k},$ (9)

for $\nabla u^{\varepsilon_{n-1}}_{k}$ not identically zero. Inserting the
approximation (IV) into (8) leads to the following iteration. Let
$u^{\varepsilon_{0}}_{1}$ and $u^{\varepsilon_{0}}_{2}$ denote initial guesses
and let $\beta^{\varepsilon_{n-1}}_{k}$ denote the normalization of $\nabla
u^{\varepsilon_{n-1}}_{k}$, then for $n=1,2,\dots$, we solve

$\displaystyle-\varepsilon_{n}\Delta
u^{\varepsilon_{n}}_{1}+\beta^{\varepsilon_{n-1}}_{1}\cdot\nabla
u^{\varepsilon_{n}}_{1}f_{1}+\phi_{1}u^{\varepsilon_{n}}_{1}$
$\displaystyle=K_{1}+\lambda R+\phi_{1}u^{\varepsilon_{n}}_{2}$ (10)
$\displaystyle-\varepsilon_{n}\Delta
u^{\varepsilon_{n}}_{2}+\beta^{\varepsilon_{n-1}}_{2}\cdot\nabla
u^{\varepsilon_{n}}_{2}f_{2}+\phi_{2}u^{\varepsilon_{n}}_{2}$
$\displaystyle=K_{2}+\phi_{2}R+\phi_{2}u^{\varepsilon_{n}}_{1}$

for $x\in\Omega$ with $u^{\varepsilon_{n}}_{1}=0$ for $x\in G$ and
$u^{\varepsilon_{n}}_{2}=R+u^{\varepsilon_{n}}_{1}$ for $x\in D$. Equation
(10) now resembles a convection-diffusion problem for which a number of
numerical schemes exist. In this work, we formulate the discrete system via a
streamline diffusion approximation.

###  IV-A Streamline Diffusion

For the discretization of equation (10) we use a streamline diffusion
approximation [32, 33] and in our derivation we follow the settings in [34].
We seek an approximate solution which is a piecewise linear, continuous
function on a given simplicial mesh $\mathcal{T}_{h}$ covering exactly the
computational domain $\Omega$. Namely, letting $T$ denote a generic element we
have: $\overline{\Omega}=\bigcup_{T\in\mathcal{T}_{h}}T$. Further, assume that
the sets $G$ and $D$ are partitioned by the mesh exactly. Then for $\Gamma$
denoting the set $G$ or $D$, the piecewise linear finite element space
containing our solutions is given by $V_{\Gamma}(h)$

$\displaystyle V_{\Gamma}(h):=\left\\{v\in C(\overline{\Omega})\;|\;v|_{T}\in
P_{1}(T),\;\forall T\in\mathcal{T}_{h},v|_{\Gamma}=g\right\\}$

for the appropriate function $g$, and where $P_{1}(T)$ denotes the space of
polynomials of degree less than or equal to one over $T$. Note that the
solution spaces for $u_{1}$ and $u_{2}$ differ slightly due to the sets $G$
and $D$, and are given by $V_{G}(h)$ and $V_{D}(h)$, respectively.

For ease of presentation in the derivation of the weak form of (10), we
suppress the superscript notation as well as the iteration variable $n$. The
test functions for the streamline diffusion formulation have the form
$v+\theta h\nu\beta\cdot\nabla v$ for $v\in V_{\Gamma}(h)$, where
$\nu=(|\beta|^{2}+1)^{-1/2}$. The variable $h$ denotes the mesh size and
$\theta$ is a positive constant. Consider the first equation in (10), then for
$v_{1}\in V_{D}(h)$ multiplying by a test function and integrating yields

$\displaystyle\left(-\varepsilon\Delta u_{1}+\beta_{1}\cdot\nabla
u_{1}f_{1}+\phi_{1}u_{1},v_{1}+\theta h\nu\beta_{1}\cdot\nabla v_{1}\right)$
(11) $\displaystyle=\left(\phi_{1}u_{2},v_{1}+\theta h\nu\beta_{1}\cdot\nabla
v_{1}\right)$ $\displaystyle+\left(K_{1}+\lambda R,v_{1}+\theta
h\nu\beta_{1}\cdot\nabla v_{1}\right)$

where $(u,v)=\int_{\Omega}uvd\omega$. We remark that because $v\in
V_{\Gamma}(h)$, the components of $\nabla v$ are constant on each element and
therefore integrations of the form $(w,\theta h\nu\beta\cdot\nabla v)$ are
considered elementwise. Working from right to left, define the linear
functional $l_{1}(\cdot)$ by

$\displaystyle l_{1}(v_{1}):=\left(K_{1}+\lambda R,v_{1}+\theta
h\nu\beta_{1}\cdot\nabla v_{1}\right)$ (12)

and the bilinear form $b_{1}(\cdot,\cdot)$ by

$\displaystyle b_{1}(u_{2},v_{1}):=\left(\phi_{1}u_{2},v_{1}+\theta
h\nu\beta_{1}\cdot\nabla v_{1}\right).$ (13)

The final term can be split into the sum of the viscosity term and its
remainder. Then for $\beta_{1}=([\beta_{1}]_{1},[\beta_{1}]_{2})^{T}$ the
introduction of the matrix

$\displaystyle{\bf
B}_{1}=\begin{pmatrix}[\beta_{1}]_{1}^{2}&[\beta_{1}]_{1}[\beta_{1}]_{2}\\\
[\beta_{1}]_{1}[\beta_{1}]_{2}&[\beta_{1}]_{2}^{2}\end{pmatrix}f_{1}$

allows us to represent the next term as

$\displaystyle\left(\beta_{1}\cdot\nabla u_{1}f_{1}+\phi_{1}u_{1},v_{1}+\theta
h\nu\beta_{1}\cdot\nabla v_{1}\right)$ (14) $\displaystyle=\theta
h\nu\left(\nabla u_{1},{\bf B}_{1}\nabla
v_{1}\right)+\left(\beta_{1}\cdot\nabla u_{1}f_{1},v_{1}\right)$
$\displaystyle+\left(\phi_{1}u_{1},v_{1}+\theta h\nu\beta_{1}\cdot\nabla
v_{1}\right).$

Performing integration by parts on the viscosity term yields

$\displaystyle-\varepsilon\left(\Delta u_{1},v_{1}+\theta
h\nu\beta_{1}\cdot\nabla v_{1}\right)$ (15) $\displaystyle=\varepsilon(\nabla
u_{1},\nabla v_{1})-\varepsilon\left(\nabla u_{1},v_{1}+\theta
h^{p}\nu\beta_{1}\cdot\nabla v_{1}\right)_{\partial\Omega}$

where $(\nabla u,v)_{\Gamma}=\int_{\Gamma}v\nabla u\cdot ndS$. We note that
due to the definition of $V_{\Gamma}(h)$, all higher derivatives of $\nabla
v_{1}$ and $\beta_{1}$ are zero for any given element. In this case,
integration by parts results in an additional term due to the set $D$ not
necessarily being the entire (or part) of the boundary of the domain. By
letting ${\bf K}_{1}={\bf I}\varepsilon+\theta h\nu{\bf B}_{1}$ we can define
the bilinear form $a_{1}(\cdot,\cdot)$ as

$\displaystyle a_{1}(u_{1},v_{1})$ $\displaystyle:=\left(\nabla u_{1},{\bf
K}_{1}\nabla v_{1}\right)+\left(\beta_{1}\cdot\nabla u_{1}f_{1},v_{1}\right)$
(16) $\displaystyle\hskip 18.06749pt+\left(\phi_{1}u_{1},v_{1}+\theta
h\nu\beta_{1}\cdot\nabla v_{1}\right)$ $\displaystyle\hskip
18.06749pt-\varepsilon\left(\nabla u_{1},v_{1}+\theta h\nu\beta_{1}\cdot\nabla
v_{1}\right)_{\partial\Omega}$

which represents the combination of (14) and (15). Using the linear and
bilinear forms defined above, equation (11) becomes

$\displaystyle a_{1}(u_{1},v_{1})-b_{1}(u_{2},v_{1})$
$\displaystyle=l_{1}(v_{1}),\quad\forall v_{1}\in V_{D}(h).$

This provides the weak form for the first equation in (10). By following the
same derivation and introducing analogous bilinear and linear forms we arrive
at the weak form of (10): find $(u_{1},u_{2})\in(V_{D},V_{G})$ such that

$\displaystyle a_{1}(u_{1},v_{1})-b_{1}(u_{2},v_{1})$
$\displaystyle=l_{1}(v_{1}),\quad\forall v_{1}\in V_{D}(h)$ (17)
$\displaystyle a_{2}(u_{2},v_{2})-b_{2}(u_{1},v_{2})$
$\displaystyle=l_{2}(v_{2}),\quad\forall v_{2}\in V_{G}(h).$

The system (17) then represents our finite element formulation for a given
$\varepsilon$. This results in a single linear system which solves for the
epsilon dependent $u_{1}$ and $u_{2}$ simultaneously using the algebraic
multilevel solver Multigraph 2.1 (cf. [35, 36]). We remark that this is in
contrast to the value iteration approach in [10] which iterates between each
equation in the coupled system. In our setting, the only iteration is due to
decreasing $\varepsilon$.

##  V Numerical Example

We consider the road scenario depicted in Figure 1. Here, a pair of disabled
vehicles are present on the shoulder of the road due to a collision.

Figure 1: Road scenario

The assumption is that a breakdown is more likely in regions near the pair of
disabled vehicles due to potential debris etc. which the AV must account for
in navigating to the goal location indicated by the green dot; in this case,
we assume the goal location and repair depot coincide. This increased
likelihood of a breakdown is represented in the model by the function
$\phi(x)$ given in Figure 2.

Figure 2: Counter plot of the spatially dependent likelihood of breakdown
function $\phi(x)$ associated with the road layout in Figure 1.

For the above scenario we define $\Omega=[0,2]\times[0,1]$, let
$G=D=(1.9,0.5)$, and uniformly refine a total of seven times. We set
$f_{1}=1$, $f_{2}=0.2$, $R=1$, $K_{1}=1$, $K_{2}=1$ and define
$\phi(x)=7e^{-5(x-1)^{2}-5*y^{2}}$ which is plotted above. Additionally, we
take $\lambda=1$ and $\phi_{2}(x)=3$. Figures 3 and 4 display the contour
plots associated with the value functions which correspond to vehicle
operation in modes 1 and 2, respectively. A trajectory from any starting point
$(x_{0},y_{0})$ is obtained by moving in the direction of the gradient. The
contours show that the optimal trajectories move away from the disabled
vehicles as expected.

Figure 3: Value function $u_{1}$ corresponding to normal operation. Figure 4:
Value function $u_{2}$ which corresponds to operation after a partial
breakdown.

##  VI Conclusions and future work

In this work, we considered the model of optimal path planning with random
breakdowns introduced in [10] and analyzed from a theoretical point of view.
We proved a comparison principle for viscosity solutions of the weakly coupled
system of eikonal equations of the value functions. In particular, we
described how to bypass the non-trivial degenerate coupling condition (5).
Then we showed an example of how the lack of boundary conditions on the Aubry
set can compromise the uniqueness of the solution. Finally, in the same
example, we have provided boundary conditions for which the system does not
have solutions. It still remains an open question how to choose boundary
conditions to ensure the existence of a solution. Indeed, the few existence
results proved for weakly coupled systems require that the coupling matrix is
non-degenerate [19], or that the Aubry set is empty [37], or that the domain
is a torus [14].

Future direction may also include the extension of the model in [10] to
domains topologically more challenging, like networks, assuming not only that
the vehicle is subject to breakdowns but also that the speed varies because of
heterogeneous road conditions or the presence of lights at the nodes.
Numerically, this work took advantage of finite element schemes used in
convection-diffusion problems by iteratively solving a linearized form of the
system augmented with artificial viscosity (10); specifically, the streamline
diffusion approximation was used as a means of obtaining stable numerical
solutions at each iteration. A number of schemes exist for such problems and
the exploration of their utility in the context of the coupled system
considered in this paper is one potential direction for future work.

## References

  * [1] P. Tompkins, _Mission-directed path planning for planetary rover exploration_. Carnegie Mellon University, 2005.
  * [2] D. Nieto-Hernández, C.-F. Méndez-Barrios, J. Escareno, V. Ramírez-Rivera, L. Torres, and H. Méndez-Azúa, “Non-holonomic flight modeling and control of a tilt-rotor mav,” in _2019 6th International Conference on Control, Decision and Information Technologies (CoDIT)_ , 2019, pp. 1947–1952.
  * [3] A. Aguiar and A. Pascoal, “Regulation of a nonholonomic autonomous underwater vehicle with parametric modeling uncertainty using lyapunov functions,” in _Proceedings of the 40th IEEE Conference on Decision and Control (Cat. No.01CH37228)_ , vol. 5, 2001, pp. 4178–4183 vol.5.
  * [4] J. Reeds and L. Shepp, “Optimal paths for a car that goes both forwards and backwards,” _Pacific journal of mathematics_ , vol. 145, no. 2, pp. 367–393, 1990.
  * [5] K. Alton and I. M. Mitchell, “Optimal path planning under different norms in continuous state spaces,” in _Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006. ICRA 2006._ IEEE, 2006, pp. 866–872.
  * [6] B. Chazelle, “Approximation and decomposition of shapes,” _Algorithmic and Geometric Aspects of Robotics_ , vol. 1, pp. 145–185, 1985.
  * [7] E. W. Dijkstra, “A note on two problems in connexion with graphs:(numerische mathematik, 1 (1959), p 269-271),” 1959.
  * [8] L. E. Kavraki, P. Svestka, J.-C. Latombe, and M. H. Overmars, “Probabilistic roadmaps for path planning in high-dimensional configuration spaces,” _IEEE Transactions on Robotics and Automation_ , vol. 12, no. 4, pp. 566–580, 1996.
  * [9] S. LaValle, “Rapidly-exploring random trees: A new tool for path planning,” _Research Report 9811_ , 1998.
  * [10] M. Gee and A. Vladimirsky, “Optimal path-planning with random breakdowns,” _IEEE Control Systems Letters_ , vol. 6, pp. 1658–1663, 2021.
  * [11] H. T. Davis, _Introduction to nonlinear differential and integral equations_. US Atomic Energy Commission, 1960.
  * [12] W. H. Fleming and H. M. Soner, _Controlled Markov processes and viscosity solutions_. Springer Science & Business Media, 2006, vol. 25.
  * [13] M. H. Protter and H. F. Weinberger, _Maximum principles in differential equations_. Prentice Hall, Englewood Cliffs, 1967.
  * [14] F. Camilli, O. Ley, P. Loreti, and V. D. Nguyen, “Large time behavior of weakly coupled systems of first-order hamilton–jacobi equations,” _Nonlinear Differential Equations and Applications NoDEA_ , vol. 19, pp. 719–749, 2012.
  * [15] F. Cagnetti, D. Gomes, and H. V. Tran, “Adjoint methods for obstacle problems and weakly coupled systems of pde,” _ESAIM: Control, Optimisation and Calculus of Variations_ , vol. 19, no. 3, pp. 754–779, 2013.
  * [16] H. Mitake and H. V. Tran, “Remarks on the large time behavior of viscosity solutions of quasi-monotone weakly coupled systems of hamilton–jacobi equations,” _Asymptotic Analysis_ , vol. 77, no. 1-2, pp. 43–70, 2012.
  * [17] H. Mitake, A. Siconolfi, H. V. Tran, and N. Yamada, “A lagrangian approach to weakly coupled hamilton–jacobi systems,” _SIAM Journal on Mathematical Analysis_ , vol. 48, no. 2, pp. 821–846, 2016.
  * [18] A. Davini and M. Zavidovique, “Aubry sets for weakly coupled systems of hamilton–jacobi equations,” _SIAM Journal on Mathematical Analysis_ , vol. 46, no. 5, pp. 3361–3389, 2014.
  * [19] H. Engler and S. M. Lenhart, “Viscosity solutions for weakly coupled systems of hamilton-jacobi equations,” _Proceedings of the London Mathematical Society_ , vol. 3, no. 1, pp. 212–240, 1991.
  * [20] H. Ishii and S. Koike, “Viscosity solutions for monotone systems of second–order elliptic pdes,” _Communications in partial differential equations_ , vol. 16, no. 6-7, pp. 1095–1128, 1991.
  * [21] M. Bardi and I. Capuzzo-Dolcetta, _Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations_ , ser. Systems & Control: Foundations & Applications. Birkhäuser Boston, Inc., Boston, MA, 1997, with appendices by Maurizio Falcone and Pierpaolo Soravia. [Online]. Available: https://doi.org/10.1007/978-0-8176-4755-1
  * [22] J. A. Sethian, “A fast marching level set method for monotonically advancing fronts.” _proceedings of the National Academy of Sciences_ , vol. 93, no. 4, pp. 1591–1595, 1996.
  * [23] C.-Y. Kao, S. Osher, and Y.-H. Tsai, “Fast sweeping methods for static hamilton–jacobi equations,” _SIAM journal on numerical analysis_ , vol. 42, no. 6, pp. 2612–2632, 2005.
  * [24] A. Caboussat, R. Glowinski, and T.-W. Pan, “On the numerical solution of some eikonal equations: An elliptic solver approach,” _Chinese Annals of Mathematics, Series B_ , vol. 36, no. 5, pp. 689–702, 2015.
  * [25] Y. Yang, W. Hao, and Y.-T. Zhang, “A continuous finite element method with homotopy vanishing viscosity for solving the static eikonal equation,” _Communications in Computational Physics_ , vol. 31, no. 5, pp. 1402–1433, 2022.
  * [26] J.-L. Guermond, F. Marpeau, and B. Popov, “A fast algorithm for solving first-order pdes by l1-minimization,” _Communications in Mathematical Sciences_ , vol. 6, no. 1, pp. 199–216, 2008.
  * [27] D. Flad, A. Pradhan, and S. Murman, “Arbitrary order solutions for the eikonal equation using a discontinuous galerkin method,” _arXiv preprint arXiv:2108.05950_ , 2021.
  * [28] K. W. Morton, _Revival: Numerical solution of convection-diffusion problems (1996)_. CRC Press, 2019.
  * [29] E. Cartee, A. Farah, A. Nellis, J. Van Hook, and A. Vladimirsky, “Quantifying and managing uncertainty in piecewise-deterministic markov processes,” _SIAM/ASA Journal on Uncertainty Quantification_ , vol. 11, no. 3, pp. 814–847, 2023.
  * [30] M. G. Crandall, L. C. Evans, and P.-L. Lions, “Some properties of viscosity solutions of hamilton-jacobi equations,” _Transactions of the American Mathematical Society_ , vol. 282, no. 2, pp. 487–502, 1984.
  * [31] G. Namah and J.-M. Roquejoffre, “Remarks on the long time behavior of the solutions of hamilton-jacobi equations,” _Communications in Partial Differential Equations_ , vol. 24, no. 5-6, pp. 883–893, 1999. [Online]. Available: https://doi.org/10.1080/03605309908821451
  * [32] A. N. Brooks and T. J. R. Hughes, “Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations,” _Comput. Methods Appl. Mech. Engrg._ , vol. 32, no. 1-3, pp. 199–259, 1982, fENOMECH ”81, Part I (Stuttgart, 1981). [Online]. Available: https://doi.org/10.1016/0045-7825(82)90071-8
  * [33] C. Johnson, _Numerical solution of partial differential equations by the finite element method_. Cambridge University Press, Cambridge, 1987.
  * [34] R. E. Bank, P. S. Vassilevski, and L. T. Zikatanov, “Arbitrary dimension convection–diffusion schemes for space-time discretizations,” _Journal of Computational and Applied Mathematics_ , vol. 310, pp. 19–31, 2017.
  * [35] R. E. Bank and R. K. Smith, “An algebraic multilevel multigraph algorithm,” _SIAM Journal on Scientific Computing_ , vol. 23, no. 5, pp. 1572–1592, 2002\.
  * [36] R. E. Bank. (2015) Multigraph 2.1. [Online]. Available: https://ccom.ucsd.edu/~reb/software.html
  * [37] F. Camilli and P. Loreti, “Comparison results for a class of weakly coupled systems of eikonal equations,” _Hokkaido Mathematical Journal_ , vol. 37, no. 2, pp. 349 – 362, 2008. [Online]. Available: https://doi.org/10.14492/hokmj/1253539559