# Eco-Routing Using Open Street Maps

R. K. Ghosh  Department of EECS, IIT Bhilai, Raipur 492015, INDIA
<EMAIL_ADDRESS>Vinay R  Arnab Bhattacharya

###### Abstract

A vehicle’s fuel consumption depends on its type, the speed, the condition,
and the gradients of the road on which it is moving. We developed a Routing
Engine for finding an eco-route (one with low fuel consumption) between a
source and a destination. Open Street Maps has data on road conditions. We
used CGIAR-CSI road elevation data [4] to integrate the road gradients into
the proposed route-finding algorithm that modifies Open Street Routing Machine
(OSRM). It allowed us to dynamically predict a vehicle’s velocity, considering
both the conditions and the road segment’s slope. Using Highway EneRgy
Assessment (HERA) methodology, we calculated the fuel consumed by a vehicle
given its type and velocity. We have created both web and mobile interfaces
through which users can specify Geo coordinates or human-readable addresses of
a source and a destination. The user interface graphically displays the route
obtained from the proposed Routing Engine with a detailed travel itinerary.

##  I Introduction

The stiffest challenge faced by humankind today is to maintain a low ambient
Air Quality Index (AQI) [26]. Besides having an adverse impact on climate [6,
10, 30], poor air quality affects flora and fauna. The presence of chemicals,
biological, and physical particles in the atmosphere is the major reason for
air pollution. It is not possible to prevent naturogenic pollution. However,
it remains under the tolerance limit over time. The pollution due to
anthropogenic sources increases in proportion due to accelerated human-
activities [14]. The significant causes of pollution are toxic emissions from
industries, agricultural activities, waste accumulation, and fossil fuel
combustion in the generation of electricity and transportation service.

The emission due to the burning of fuel in the transportation sector is a
dominant factor affecting air quality, especially in large cities. It releases
Green House Gases (GHG) and Non-Methane Volatile Organic Compounds (NMVOC)
into the atmosphere. NMVOC mainly contributes to particulate matter.
Controlling particulate matter is a significant challenge especially in India.
Transportation sources lead to almost one-third of PM pollution in India. In
fact, ten most polluted cities in the world are located in India. Figure 1
[27] gives an idea of the severity levels of PM pollution in different cities
across northern India.

Figure 1: Contribution of the transport sector to total emissions of the main
air pollutants [1].

The composition of direct GHG is water vapor (H2O), Nitrous Oxide (NO), Carbon
Dioxide (C2O), Methane (CH4), Ozone (O3). Indirect GHG has NMVOC, Sulpher
Dioxide (S2O), and Carbon Monoxide (CO) in its composition. Without GHG, the
average surface temperature of Earth would be
-18$\mathrm{\SIUnitSymbolDegree}$ instead of the present average of
15$\mathrm{\SIUnitSymbolDegree}$ [21]. However, excessive emission of indirect
GHG may cause a significant rise in global temperature and adversely affecting
the Earth’s climate [33]. As one can observe from the Fig 2, the most
significant factor contributing to emissions of all the major air pollutants
except SO2 happens to be road transportation system.

Figure 2: Contribution of the transport sector to total emissions of the main
air pollutants [1].

The severity of pollution in Indian cities has created a near-disastrous
situation for human survival. The economic growth coupled with population
growth in cities has also led to increased mobility. The number of registered
vehicles in India has multiplied in the last decade and is still growing.
According to a report [37], during the eight years 2011-2018, the increase in
sales of two-wheelers alone is from 117 million to 210 million as indicated by
Fig. 3.

Figure 3: Two wheeler sales in India during the period 2011-2018.

Not only does the fuel used in the residential and transportation sectors
contribute to the air pollutants, but the movements of vehicles also cause re-
suspension of dust, keeping the density of the ambient concentrations of
particulate matter more or less at the same level.

The minimum travel time is the most important optimization in measuring the
efficiency of transportation services. Geographic Information systems (GIS)
help in achieving it. All the web/mobile GIS systems provide the shortest
distance or the shortest time routes by collecting real-time or cumulative
traffic patterns. The fuel consumption of a vehicle is not a linear function
of velocity. So, optimizing the route for the shortest possible time does not
necessarily optimize a vehicle’s fuel consumption. It is essential to work on
the tools that reduce emissions from transportation services. The cumulative
effects of creating such a tool-chain can be summarized as follows:

  * •

It lowers the cost of transport and contributes to the overall fiscal
management of a nation’s economy.

  * •

It reduces the emission improves AQI in cities, making them more livable.
Overall impacts lead to a better environment.

  * •

It leads to the conservation of fossil fuel, which is estimated to last only
between 50-100 years at the present rate of consumption.

This paper provides a method for estimating fuel consumed by the vehicle for a
given path. By incorporating the estimation of fuel with the technique of
finding the shortest route between a pair of source and destination, we can
estimate the fuel paths for traveling on vehicles on the roads. We have
developed a web and a mobile interface to find a complete travel itinerary by
considering the fuel consumption factor.

##  II Background and Related Work

Eco-friendly routing of vehicular movements is one of the significant aspects
of the overall approach for control of emission and conservation of fossil
fuel. Before addressing the problem, we need to understand the background of
appropriate models used to estimate emission and fuel consumption for
different vehicle types. The concern here is to integrate both fuel
consumption assessments and the road elevation data with widely use routing
machines such as Google Maps or Open Street Maps to generate an eco-friendly
route between a source and a destination.

###  II-A Eco-friendly Routing

Petroleum-based fuel accounts for almost 35% of the input cost in the
transportation industry [12, 23]. Any fuel-saving measures will be a boon for
the industry and cut down the emissions drastically. Eco-routing solutions
[35] allow for a prediction of the most fuel-economic routes.

Many factors affect vehicle emissions and fuel consumption. However, the three
most important factors are vehicle characteristics, road characteristics, and
traffic conditions. Each type of vehicle has different aerodynamic drags,
different engine technology, which determines its fuel efficiency. Fuel
consumption and vehicle emissions are sensitive to changes in travel speed,
traffic patterns, etc. A small portion of the route may account for a
substantial part of total fuel consumption depending on traffic intensity
[38]. The gradient of a road plays another important role in understanding a
vehicle’s fuel efficiency as real-world fuel consumption and emission models
are sensitive to the road grade estimates. Some research [40] indicates that
the road grade accounts for approximately 1-3% of fuel consumption used in
light vehicles. The eco-friendly route will be the least fuel-consuming route
among all the possible routes given a source and a destination. Therefore, the
estimate of fuel consumption requires, so we need the entire road curvature
data for the complete path from a source to a destination. Finally, it is
extremely convenient and a practical aid for a vehicle driver if Petroleum-
based fuel accounts for almost 35% of the input cost in the transportation
industry [12, 23]. Any fuel-saving measures will be a boon for the industry
and cut down the emissions drastically. Eco-routing solutions [35] allow for a
prediction of the most fuel-economic routes.

Many factors affect vehicle emissions and fuel consumption. However, the three
most important factors are vehicle characteristics, road characteristics, and
traffic conditions. Each type of vehicle has different aerodynamic drags,
different engine technology, which determines its fuel efficiency. Fuel
consumption and vehicle emissions are sensitive to changes in travel speed,
traffic patterns, etc. A small portion of the route may account for a
substantial part of total fuel consumption depending on traffic intensity
[38]. The gradient of a road plays another important role in understanding a
vehicle’s fuel efficiency as real-world fuel consumption and emission models
are sensitive to the road grade estimates. Some research [40] indicates that
the road grade accounts for approximately 1-3% of fuel consumption used in
light vehicles. The eco-friendly route will be the least fuel-consuming route
among all the possible routes, given a source and destination. Therefore, the
estimate of fuel consumption requires, so we need the entire road curvature
data for the complete path from a source to a destination. Finally, it is
extremely convenient and a practical approach for a vehicle driver if the
existing navigation tools may provide eco-route for commutation and
transportation. Finding eco-friendly routing between a source and a
destination requires three basic components:

  * •

A model for estimating emissions from different vehicle types and a
compilation method for calculating aggregated emission data under different
road surfaces and traffic conditions.

  * •

A database of roads providing a digital road map with historical data and data
fusion strategies for integrating crowdsourced real-time traffic.

  * •

A user interface for taking input and display the route on a road map, along
with fuel consumed, and travel itinerary

Developing an eco-friendly navigation tool is far too complex to build from
scratch. It requires a sustained long term large team effort, which can only
be possible by Google, Microsoft, etc. Therefore, in this paper, we employed
the following short cuts:

  1. 1.

For road network databases, we leveraged two well-established navigation and
route machines, namely, Google Maps [18] and Open Street Maps [5].

  2. 2.

For estimating fuel consumption, we used HERA methodology [36].

  3. 3.

For display of maps and least fuel cost routes, we used a third-party front-
end builder.

  4. 4.

For compiling travel itinerary, we used Google multipoint representation of
turns, etc.

  5. 5.

For integrating a real-time mobile navigation interface, we employed another
third party application.

###  II-B HERA Methodology

Many models exist for the estimation of emissions due to road traffic and
energy consumption. We can analyze the traffic flows from two general
perspectives,

  1. 1.

Microscopic.

  2. 2.

Macroscopic.

A Microscopic view of traffic takes into account the vehicle-related variables
and traffic flow characteristics. The vehicle-related variables, among others,
include the length of the vehicle, the speed, and the acceleration. The
traffic flow characteristics refer to the variability in traffic streams,
which depends on physical constraints and the complex behavior of driving. The
distance between two vehicles and time variability between different traffic
streams constitutes variability in a traffic stream. For macroscopic analysis
of traffic, we try to extrapolate an aggregate picture of vehicular movements
[24] instead of looking into each vehicle in a traffic stream.

Two prominent simulators widely used to model emission from vehicles are:

  * •

The COPERT [29] – COmputer Programme to calculate Emissions from Road
Transport.

  * •

MOVES [3] – Motor Vehicles Emission Simulator.

Both models use their databases to obtain the emission factors as a function
of the average cycle of speeds. COPERT is mainly used in Europe, while MOVES
is used in the USA. Other types of simulators take microscopic details into
account to estimate emission maps with speed and acceleration.

The energy consumption and emission modeling tools provide quantitative
evaluation for a given route. For our purposes, we have used models created
based on HERA [36] methodology. It uses an average speed consumption model
based on COPERT emission factors. These factors are adjusted with a correction
factor to incorporate the vehicles’ gradient effect in each road section.

The average speed is a key factor in fuel consumption and emissions of a
vehicle. Several studies [36] have found a relationship between fuel
consumption and vehicle speed. Therefore, it is very appropriate to use models
which are based on speed. Since the HERA consumption model is based on
consumption and emission factors from COPERT, it considers a set of vehicle-
centric parameters to calculate fuel consumption. For example, it uses the
fuel type and engine capacity, and engine technology for light vehicles
(passenger cars). It uses the maximum total weight, load state, engine
technology, and road gradient for heavy vehicles.

HERA [36] uses a slightly altered set of parameters. It take the average speed
as the leading factor and the road gradient as a correction factor. It
considers six vehicle types, namely, passenger cars, light-duty vehicles,
motorcycles, rigid trucks, articulated trucks, and buses. HERA methodology
provides an analytical formula to determine an approximates relationship
between the fuel Consumption and the average speed for each category of
vehicle. Table I summarizes the relationships.

Vehicle Type | Fuel Consumption Function for 0% Road Gradient (gfuel/veh-km)
---|---
Motor Cycle | 25.722 + (276.13/V) + (-0.254)*V + 0.00311*V2
Passenger Car | 54.7 + (496/V) + (-0.542)*V + 0.0042*V2
Light-Duty Vehicle | 146.27 + ((-0.0000106)/V) + (-2.596)*V + 0.01984*V2
Rigid Truck | 152.96 + (604.156/V) + (-2.295)*V + 0.0238*V2
Articulated Truck | 332.603 + (1680.879/V) + (-4.676)*V + 0.0311*V2
Bus | 281.735 + (4186.178/V) + (-3.457)*V + 0.0216*V2
TABLE I: Approximate relationships between speed and fuel consumption for each
vehicle type.

Overall Fuel Consumption,

$FC_{i}=FCZ_{i}*RGF_{i},$

where $FCZ$ is the fuel consumption on the road with zero gradients, $RGF$ is
the Road Gradient Factor, and $i$ represents a different vehicle type such as
a car, bus, etc.

HERA uses each vehicle category’s consumption factors and consumption curves
data from COPERT database to construct the above formulae. By applying the
non-linear regression techniques and parametrization of the consumption
curves, HERA derives the consumption model function for each vehicle type,
excluding the Road Gradient Factor RGF.

###  II-C Road Gradient Factor

Apart from average vehicle speeds, many other parameters affect fuel
consumption. These include a wear-and-tear of vehicle, the technology used,
and the conditions of the road, the elevation of the road (or Road Gradient),
Road Surface.

The gradient of a road increases or decreases the resistance of a vehicle to
traction. The increase or decrease in the engine load has a corresponding
effect on emission and fuel consumption rates. The effect of road gradient on
the vehicle is not symmetrical, i.e., the extra fuel consumed by traveling
uphill won’t be equal to the reduced fuel consumption when going downhill. In
general, the road gradient effect is more on the lighter vehicle than on the
heavy-duty vehicles.

For a given speed $V$, the fuel consumed on the flat road is different
compared to either uphill or downhill. HERA uses a correction factor which is
based on the results of the German Emission Factor Program [8]. Some of the
factors included in the correction are:

  1. 1.

The road gradient

  2. 2.

The pollutant

  3. 3.

The average speed of the vehicle

  4. 4.

The vehicle category

The road gradient factor $RGF_{i}$ for each vehicle type $i$ is found
according to following formula:

$\begin{split}RGF_{i}&=A6_{i,k}*V^{6}+A5_{i,k}*V^{5}+A4_{i,k}*V^{4}+A3_{i,k}*V^{3}\\\
&+A2_{i,k}*V^{2}+A1_{i,k}*V+A0,\end{split}$

where RGFi is the road gradient factor for the vehicle category $i$, $V$ is
the vehicle’s average speed in km/hr. $A0$ to $A6$ are the coefficients
depending on the gradient and the type of vehicle provided by Hickman et al.
[17].

##  III Google Maps and Open Street Map

Developing an eco-friendly vehicle routing from scratch is not feasible
without the support of a large team. We aim to assess problems likely to
appear in trying to enhance an existing navigation tool for eco-friendliness.
Therefore, we started by exploring the different APIs provided by Google Maps
and Open Street Map.

###  III-A Geo-Coding

Let us first start with a precise definition of the problem.

######  Definition III.1 (Problem statement).

Given a source $S$ and a destination $D$, find the most fuel-efficient route
from $S$ to $D$.

For navigation, a user may specify a source or a destination. The source or
the destination would be in the form of a human-readable address such as a
street address. A human-readable street address is converted into geographic
coordinates by a process known as Geo-coding. More precisely, given a street
address as input, Geo-coding outputs a (latitude, longitude) pair. The
geographic coordinates are convenient as they can be used to place markers on
a map and find the shortest distance or elevation between a source and a
destination.

###  III-B Google Map APIs

Google Maps API tool provides an API for extracting geo-coordinates for an
address by using the following url:

https://maps.googleapis.com/maps/api/geocode/json?address=<ADDRESS>&key=<KEY>

Along with latitude and longitude coordinates, it also provides other address
components such as state, country, pin code, etc. For example, using the
ADDRESS = ”cse IIT Kanpur” returns Geo-coordinates (26.5143, 80.2348).

####  III-B1 Places Nearby

Google API can be used to find important utility points near a pre-specified
location. For example, utility points of interest could be the nearby
restaurants, car service stations, hospitals, fuel pumps, convenience stores,
govt offices, etc. The API call given below extracts the location and other
information about the nearby restaurants.

https://maps.googleapis.com/maps/api/place/nearbysearch/json?location=26.5143,80.2348&rankby=distance&type=restaurant&key=<KEY>
The above example returns the restaurant details near the given coordinates
(CSE IIT Kanpur), sorted by the distance.

####  III-B2 Elevation

The elevation values represent the elevations relative to the sea level. For
finding the elevation at some Geo coordinates, we use the following GoogleMap
API.

https://maps.googleapis.com/maps/api/elevation/json?locations=26.5143,80.2348&key=<KEY>
It returns the elevation at CSE IIT Kanpur.

####  III-B3 Direction

A series of Geo coordinates along the path represents a direction. To extract
the Geo coordinate series for a given source-destination pair, we use the
following GoogleMap API. We need to go via these Geo coordinates while
driving.

https://maps.googleapis.com/maps/api/directions/json?origin=<GeoCordinates>&destination=<GeoCordinates>&key=<KEY>
This API returns polylines (a sequence of straight lines) from a source to a
destination. This polyline contains encrypted geo coordinates. We need to
decrypt them and extract a series of geo-coordinates from all those polylines.

Figure 4 shows the graphical display of the direction from CSE IIT Kanpur to
Lucknow Airport.

Figure 4: Direction service from Google for route between CSE IIT Kanpur to
Lucknow Airport

####  III-B4 Elevation Along a Path

By using ‘getElevationAlongPath’ API, We can find the road elevations along
the path from a source to a destination. For example the top half of Fig. 5
shows the elevations along the path from ‘CSE IIT Kanpur’ to ‘Lucknow
Airport’.

Figure 5: Elevation and normalized elevation along the path

To understand how the road is, we normalized the road elevations by talking
about the minimum raise/drop in the elevation between two consecutive points
$\triangle$ e=5m. The lower half of Fig. 5 shows the normalized elevations.

Using the direction service, we can even find the distance and the time
required to travel along the path. Therefore, it is possible to calculate the
average velocity along the route. By using elevation service, we can calculate
the gradient along a road. Once the road gradient is available, the equation
II-B can be used to find the fuel consumption for a given route.

####  III-B5 Fuel Estimation

We built a simple application that uses the above APIs and calculates Fuel
consumption for the route provided by Google Maps. Fig. 6 and 7 provide
screenshots of two different source-destination pairs obtained from our
application.

Figure 6: CSE IIT Kanpur to Lucknow airport

We divided the road into intervals of 256km each to get a clear understanding
of the road’s elevation.

Figure 7: IITK to Delhi

As shown in the Fig. 7 the first 256km has is a 20m rise in elevation
incurring a fuel consumption of 14.5 liters. However, for the next part 247km,
there is a 65m rise in the elevation, and the fuel consumption increases to
18.2 liters. So, the gradient of a road segment affects fuel consumption. A
vehicle requires more fuel to travel on a steep ascent than a relatively lower
ascent.

###  III-C Limitations of Google Map API

The limitation of Google Maps is that it only estimates the fuel consumption
for the given routes. In other words, it is possible to calculate fuel
consumption if we already know the path. However, we aim to estimate the least
fuel consumption route. So the routing algorithm has to incorporate fuel
consumption as a basis for finding the path and not the time or distance.
Unfortunately, Google Maps don’t have an API which uses fuel consumption as a
basis for route computation. Furthermore, we cannot apply our routing
algorithms on Google Maps because it does not provide us the maps’ raw data.

##  IV Open Street Maps (OSM)

Limitations of Google Map APIs forced us to search for the alternative to find
raw map data. We searched the projects where the users collect the raw data.
We found the open-source collaborative mapping [13] called ”Open Street Map”
(OSM) which not only competes with commercial services like Google Maps but
also provides the raw map data. OSM data is entirely crowdsourced, and anyone
can use/edit the map data [5] at no cost at all. Structurally, OSM consists of
two different types of elements:

  * •

Node: It contains places defined with it’s Geo Coordinates

  * •

Way: It contains the roads, boundaries of the building, places. The roads have
a detailed description like the road’s name, whether the road is
motorway/trunk/residential, etc.

There is a project called Open Source Routing Machine (OSRM). It is a Routing
Engine relying on the implementation of multi-level Dijkstra’s shortest path
algorithm on a graph built from the OSM data. The edge weights of the graph
are distances or duration for covering those distances. Therefore, the routing
based on the least distance or the least travel time. However, the advantage
of the above implementation is that it allows one to fix a velocity, depending
on the road type. Table II gives a specification of those values.

htb

Road Type | velocity(km/hr)
---|---
motorway | 90
motorway_link | 45
trunk | 85
trunk_link | 40
primary | 65
primary_link | 30
secondary | 55
secondary_link | 25
tertiary | 40
tertiary_link | 20
unclassified | 25
residential | 25
living_street | 10
service | 15

TABLE II: Velocity on various roads

Because OSRM provides raw map data and the fact that it has a specification of
vehicle velocity based on road type, we decided to implement our project on
Open Street Maps and OSRM.

##  V Optimal Fuel Routing Machine (OFRM)

The road network graph built from OSM by OSRM has edge weights assigned either
as distance or as the time to those distances considering the road types. For
computing the least fuel cost path, we need to change the edge weights to fuel
consumption parameters using equation II-B, which is straightforward. However,
in the equation, we still have to incorporate Road Gradient Factor. But the
OSM data set like Google Maps doesn’t have those details. So, the elevation
needs to be computed only from available data. Let us find a method to
determine the road gradient by integrating the elevation data.

###  V-A Integrating Elevation

The gradient of a road affects the velocity of a vehicle moving on the way.
Since speed affects a vehicle’s fuel consumption, we need to know the road
gradient to estimate fuel consumption when a vehicle moves either up or down a
road gradient. We assign a penalty coefficient to velocity, which varies
according to the known gradient of the road. The penalty function is
determined based on the time required to cover a distance over flat road
versus the time needed to cover the same distance on the road with known
elevation either in ascent or in descent. A road elevation has the effect of
slowing down the vehicle while moving up or down the gradient. Hence, the
estimation of fuel consumption can be made based on the modified velocity of
the car.

Before deciding on a penalty factor, it is important to understand the road
gradient’s calculation as indicated above. Ascent or descent that a vechicle
should make is equal to the difference in the heights between the starting and
the endpoint of a road segment, as shown in Fig. 8 [39].

Figure 8: Gradient of a road segment.

Given the altitude and the length (hypotenuse) of the road segment, we can
find the elevation. The measurement of all road gradient is with reference to
the sea level. If the elevation data is available for every point on the
globe, we can easily calculate the raise or drop in elevation between two
points.

Digital Elevation Model (DEM) [41] is a digital representation of the Earth’s
topography. 3D Computer Graphics (3D-CG) represents the land surface generated
from the terrain elevation data. DEM is created by collecting elevation from
equally or unequally spaced points [28]. The collected elevation information
can be used for continuous spatial representation of any terrain. DEM just
gives elevation information of terrain, including trees, or any human-made
structures with reference to Mean Sea Level (MSL).

Depending upon how the elevation information is structured and stored, we can
classify DEM into two types:

  1. 1.

Regular square grids,

  2. 2.

Triangulated irregular networks (TIN).

Grid DEM is stored as a rastor map in which each pixel has the information
about the elevation of terrain. TIN data also includes features like peaks,
slopes, conic pits. So, TIN-DEM considered being more accurate for extracting
structural information. It is possible to obtain elevation information using
field surveys from topographic contours, aerial photographs, or satellite
imagery using photogrammetric techniques. Photogrammetry is an of three-
dimensional measurements from two-dimensional data. Recently radar
interferometric techniques and Laser altimetry have also been used to generate
a DEM. Today very fine resolution DEMs at the near-global scale are readily
available from various sources [28].

Some of the sources of global elevation data set are.

  1. 1.

Global 30 Arc-Second Elevation (GTOPO30): GTOPO [2] a global digital elevation
model (DEM) with a horizontal grid spacing of 30 arc seconds (approximately 1
kilometer). GTOPO30, completed in late 1996, was developed over three years
through a collaborative effort led by staff at the US Geological Survey’s EROS
Data Center (EDC). It requires topographic information from several rasters
and vector sources.

  2. 2.

NOAA Global Land One-km Base Elevation (GLOBE): GLOBE [15] is an
internationally designed, developed, and independently reviewed global DEM, at
a latitude-longitude grid spacing of 30 arc seconds (approx. 1km).

  3. 3.

Shuttle Radar Topography Mission (SRTM): SRTM [19] is an international
research effort that obtained DEM on a near-global scale to generate the most
complete high-resolution digital topographic database of Earth. Using two
radar antennas and a single pass, it collected sufficient data to create a
digital elevation model. This 1-arc second global DEM has a spatial resolution
of about 30 meters in most areas. Also, it covers most of the world with an
absolute vertical height accuracy of 10 meters approximately.

  4. 4.

ASTER Global Digital Elevation Model: NASA and Japan’s joint operation was the
birth of Advanced Spaceborne Thermal Emission and Reflection Radiometer
(ASTER) [31]. ASTER data are used to create detailed maps of the surface
temperature of land, emissivity, reflectance, and elevation.

  5. 5.

Lidar DEM: It is an airplane laser system that sends light pulses. It reflects
the terrain and the objects on it. LiDAR [9] transmits and receives
electromagnetic radiation of the near IR band. Unlike photogrammetry, LiDAR
data collection is not affected by sun angle. LIDAR precision depends on GPS,
IMU, and Laser measurements.

DEM | DEM Resolution | DEM Accuracy | Area of missing data
---|---|---|---
GTOPO30 [11] | 900 mts | 30m | None
NOAA Globe [16] | 900 mts | 30m | None
SRTM [19] | 90 mts and 30 mts | 10m | Topographically steep areas
ASTER GDEM [31] | 30 mts | 14 mts | Areas with constant cloud cover
Lidar [9] | 1 mts | 15 cm | available at certain areas
TABLE III: Various DEM services

From the table III, we can see that DEMs GTOPO30 and NOAA GLOBE has poor
resolution and accuracy. LIDAR gives high accuracy and resolution, but it’s
just for specialized enterprise coverage specific to certain areas. So, the
choice of SRTM data is best for achieving good accuracy, and an acceptable
resolution.

The data collected by SRTM has voids [34]. Two groups of scientists worked to
fill the voids for different purposes, the CGIAR-CSI versions and the USGS
HydroSHEDS dataset. The CGIAR-CSI version 4 provides the best global coverage
for the full resolution SRTM dataset. The HydroSHEDS dataset was generated for
hydrological applications and is suitable for consistent drainage and water
flow information. In our investigation, we have used elevation information
created by Shuttle Radar Topography Mission (SRTM), which is provided by
CGIAR-CSI.

Using CGIAR-CSI GeoPortal, we can get the elevation data of terrain in the
intervals of 5${}^{\circ}\times 5^{\circ}$ tiles. The elevation produced by
SRTM has a resolution of 90m at the equator. The elevation data are available
in both ArcInfo ASCII and GeoTiff format [4], which is easy to process for our
purpose.

Figure 9 gives an illustration of a sample SRTM data points.

Figure 9: Information about dataset

The same vertical line points have the same longitude, but different latitudes
will differ with 0.000833 degrees. Similarly, the same horizontal line points
have the same latitude, but different longitudes differ by 0.000833 degrees.

If we want the elevation at the Geo co-ordinate $P(x,y)$, we can find the
nearest data point $Q$ in the dataset and may assume $Q$’s elevation as $P$’s
elevation. It will possibly work well for city areas that are typically flat
and have little elevations at a few places. However, it won’t be accurate for
mountainous roads. The elevation values differ a lot between two points on the
stretch of road passing through mountains. So we need to generate a new data
point using the dataset.

####  V-A1 Bilinear interpolation

Generally, interpolation [7] is used to find a new data point using the known
data points. Fig 10 illustrates how interpolation is used.

Figure 10: The four corners are the known data points and the point P at which
we want to interpolate.

Let the function $E:P\rightarrow R$, where $R$ is the set of points on the
real line, and $P$ is at any point on the road. Suppose, we know the values of
set of four given points $A_{11}$, $A_{21}$, $A_{22}$, $A_{12}$. First, we
need to carry out linear interpolation along the $X$-direction that results:

$E(x,b_{1})=\frac{a_{2}-x}{a_{2}-a_{1}}E(A_{11})+\frac{x-a_{1}}{a_{2}-a_{1}}E(A_{21})$
$E(x,b_{2})=\frac{a_{2}-x}{a_{2}-a_{1}}E(A_{12})+\frac{x-a_{1}}{a_{2}-a_{1}}E(A_{22})$

Next, we need to perform linear interpolation along the $Y$-direction, that
results:

$E(x,y)=\frac{b_{2}-y}{b_{2}-b_{1}}E(x,b_{1})+\frac{y-b_{1}}{b_{2}-b_{1}}E(x,b_{2})$

We can use non-linear interpolations also.

####  V-A2 Example

Now let’s consider following the mountainous road from Tirumala to Tirupati
and compare the elevations along the path with naive, interpolation, and
Google Maps data.

Figure 11: A stretch of 260m on a mountainous road.

We first divide the path of 258m into seven segments. In the figure, the
marker $A$ is the starting position, and the marker $H$ is the final position.

$A\xrightarrow{59m}B\xrightarrow{48m}C\xrightarrow{37m}D\xrightarrow{13m}E\xrightarrow{19m}F\xrightarrow{44m}G\xrightarrow{38m}H$

Position Elevation | Rounding of to
nearest point(m) | Interpolation | Google Maps
---|---|---|---
A | 817 | 822.33 | 821.5
B | 839 | 815.27 | 812.22
C | 789 | 809.46 | 806.35
D | 800 | 804.802 | 804.167
E | 800 | 806.41 | 805.41
F | 800 | 811.72 | 810.55
G | 827 | 822.99 | 822.4
H | 827 | 827.79 | 827.1
TABLE IV: Comparing SRTM results with Google Maps

The table IV have the comparison of the SRTM elevation data without
interpolation, with interpolation and Google Map elevation. From this, we can
observe that SRTM data with interpolation gives us good accuracy. Now we can
calculate elevation at any given Geo co-ordinate. Therefore, we can now
calculate the gradient of a road segment easily. Depending on the gradient, we
calculate Road Gradient Factor. Once the gradient of the road is known, we can
use Equation II-B to change the edge weights of the graph into fuel
consumption. Therefore, we can find the least fuel path for a given source and
a destination. The next section deals with the results.

##  VI Encoding the Polyline

The route is fetched in the JSON format, where all the data are represented in
the form of strings. Since the Geo coordinates are in float values, a
precision of 5-6 decimal places with the integer part may require a string of
hundreds of characters. For example, a single value of a latitude ”13.351534”
in JSON format requires nine bytes. However, any float can fit into 4 bytes.
It implies that sending hundreds of points in JSON string format would lead to
exchanging a lot of messages. The question is, can we apply some compression
technique in encoding?

###  VI-A A Simple Idea

One can think of a simple-minded approach as follows. Suppose a float value is
always considered as accurate up to five decimal places. Then we first convert
a float value into an integer by multiplying it with $10^{5}$. Next, we
convert the resulting integer into binary and then represent each byte as the
corresponding character. However, some of these characters sometimes turn out
to be control characters such as tab, space, etc. However, if we encode such a
way that every character ASCII value between 63-126 are alphabets, brackets,
etc., we can compress a polyline.

The Polyline Encoding provided in Algorithm 1 help us to achieve it [32].

Procedure  _Polyline_Encoding( $G$)_

// G-array of Geo coordinates (Polyline) ;

String output ;

factor $\leftarrow 10^{5}$;

end $\leftarrow$ length($G$);

for  _each $i\in$ [0, end)_ do

// Encoding latitude of $i^{th}$ point;

output += ENCODE($G$[i][0]*factor);

// Encoding longitude of $i^{th}$ point;

output += ENCODE($G$[i][1]*factor);

end for

return _$output$_ ;

end

Algorithm 1 Encoding Polyline

Procedure  _Encode( $A$)_

// A is an integer String output;

$B$ = $A<<1$ ; // Left shift by 1 bit

if  _A $<$ 0_ then

$B=\sim B$; // Complement

end if

while  _B $\geq$ 0x20_ do

// Concatenation

strcat(output, ASCII(0x20 $|$ ($B$ & 0x1f)) + 63);

$B>>5$; // Right shift by 5 bits

end while

return _output;_

end

Algorithm 2 Encode a number

The above polyline compression algorithm compresses JSON data of a route
around 90kB to 15kB. When a user requests a route to our engine, it gives the
response in JSON format. Since the response data is in compressed form, it
helps to reduce network latency. To evaluate the effect of compression, we
sent a hundred routes through the network to a user continuously. It took five
seconds for an end-to-end transmission for non-compressed route data, and
three seconds compressed where the number of hops between server and the user
was one.

The above polyline compression algorithm compresses JSON data of a route
around 90kB to 15kB. When a user requests a route to our engine, it gives the
response in JSON format. Since the response data is in compressed form, it
helps to reduce network latency. To evaluate the effect of compression, we
sent a hundred routes through the network to a user continuously. It took five
seconds for an end-to-end transmission for non-compressed route data, and
three seconds compressed where the number of hops between server and the user
was one.

##  VII Results and Analysis

In this section, we discuss the results from the Optimal Fuel Routing Machine
(OFRM). At present, OFRM web applicattion accepts Geo-coordinates or the
addresses as source and destination. The user can place markers for a source
and another for a destination on the map. OFRM displays a complete description
of the route between a given source-destination pair in terms of each road
segment’s distance, the associated directions, and the turns required. It also
displays three values, namely, the distance, the time, and the fuel
consumption needed to travel in that entire path. A few screenshots from OFRM
are given below for the reader to understand the advantage of navigational
outputs from OFRM.

Figure 12: IISC Bengaluru to Bengaluru Bus Station, input supplied as address

From Fig 12, we can observe that there is a small box on the top left where
the user specifies inputs for the source and the destination. A graphical
display of the route between src to dst appears marked in thick blue. In the
right part of the display window, the travel itinerary details are provided
with exact navigational instructions for distance to be traveled for a segment
of road and the type of turns (left, right or u-turns) to be executed. To show
the map tiles, we used Mapbox APIs [25], for the detailed traveling itinerary,
we used leaflet library [22]. The user may select the source and the
destination on the map instead of providing the input address. The markers are
associated with (lat, long) pairs. So, the marker-based inputs will be
interpreted as Geo-coordinates for the routing engine. The user also has an
option to supply the input directly as Geo-coordinates.

Figure 13: Bengaluru Airport to IISC Bengaluru, selecting inputs on the map.

Fig 13 illustrates an example, where the user selected the source as
Bengaluru’s International Airport and the destination as IISC on the map. So,
the input is treated as Geo-coordinates. As a result, we can see the graphical
route and travel itinerary between the source and destination.

If a user wants to find a route that goes via some chosen places, then it is
possible to request the routing engine to display a least fuel cost path by
specifying any finite number of way points.

Figure 14: Bengaluru Airport to IISC Bengaluru via Whitefield.

For example, suppose a user wants to travel from the source (airport) to the
destination(IISC) via Whitefield. Then the user gives the three inputs, the
source, the destination, and the way point. The routing engine will give a
route from the airport to IISC via Whitefield, which we can see in the Fig 14.

#### Comparing mountainous Roads

To investigate how the road gradient affects fuel consumption, we selected
some of the mountainous roads, namely, Tirupati-Tirumala, Aluva-Munnar, Kullu-
Manali, and Kalka-Shimla.

For each road, we compared the duration of the travel from Google Maps, the
output from OFRM (our routing engine) without considering the Road Gradient,
then the output of OFRM taking the gradient of the road from SRTM dataFactor.

one-way Rd | Elevation | Google Maps | OFRM-Zero | OFRM-Gradient
---|---|---|---|---
Tirumala Up Rd - 17km | 710m gain | 32min | 16min, 1.03 Litres | 29min, 1.58 Litres
Tirumala Down Rd - 17km | 710m drop | 37min | 21min, 1.13 Litres | 38min, 1.07 Litres

two-way Rd | Elevation | Google Maps | OFRM-Zero | OFRM-Gradient
---|---|---|---|---
Aluva - Munnar Rd 18km | 570m gain | 40min | 29min, 1.31 Litres | 51min, 1.94 Litres
Aluva - Munnar Rd 18km | 570m drop | 40min | 29min, 1.31 Litres | 51min, 1.94 Litres

two-way Rd | Elevation | Google Maps | OFRM-Zero | OFRM-Gradient
---|---|---|---|---
Kullu to Manali Rd - 38km | 680m gain | 1hr7min | 29min, 2.33 Litres | 43min, 2.85 Litres
Manali to Kullu Rd - 38km | 680m drop | 1hr05min | 29min, 2.33 Litres | 43min, 2.85 Litres

two-way Rd | Elevation | Google Maps | OFRM-Zero | OFRM-Gradient
---|---|---|---|---
Kalka - Shimla Rd - 25km | 520m gain | 40min | 18min, 1.51 Litres | 36min, 1.84 Litres
Kalka - Shimla Rd - 25km | 520m drop | 40min | 18min, 1.51 Litres | 36min, 1.84 Litres

TABLE V: Comparison of several mountainous roads using Google Maps, our
routing engine with and without considering the gradient.

The above polyline compression algorithm compresses JSON data of a route
around 90kB to 15kB. When a user requests a route to our engine, it gives the
response in JSON format. Since the response data is in compressed form, it
helps to reduce network latency. To evaluate the effect of compression, we
sent a hundred routes through the network to a user continuously. It took five
seconds for an end-to-end transmission for non-compressed route data, and
three seconds compressed where the number of hops between server and the user
was one.The different columns of Table V represent the following values:

  1. 1.

Google Maps \- it gives the output from the routing engine of Google Maps. It
gives the travel time without considering the traffic.

  2. 2.

OFRM-Zero \- it gives two outputs time and fuel consumed as estimated by the
proposed routing engine without considering the gradient of the road.

  3. 3.

OFRM-Gradient \- it also gives two outputs time and the amount of fuel spent
as estimated by the proposed routing engine when the road gradient is also
considered.

From the Table V we observe that output produced by OFRM-Gradient appear more
realistic in terms of the duration of travel time as well as the estimation of
fuel consumption. Since the road gradient either ascent or descent slows down
the vehicle speed and increases the total duration of travel time compared to
the corresponding outputs obtained from OFRM-zero.

###  VII-A Mobile Application

We also designed a mobile application using open source project osmbonuspack
[20], where the user can supply a source and a destination. The application
sends a request to the OFRM routing engine and using the response from OFRM.
The application gives a graphical display of the least fuel route between the
source and the destination and the itinerary.

The application can take the current location of the user from device’s GPS
chip as input. In this case, the user don’t need to provide a source. The
application displays a marker at the user’s current location. As he changes
the location, the marker location in the map also changes correspondingly. But
we maintained the display of the route from his/her position to the
destination. If the user goes outside the route that the routing engine
provided, the mobile application sends a new request with the user’s current
position and gets corresponding route updates.

(a) Route at current location.

(b) Travel itinerary.

Figure 15: Route and travel itinerary from current location of the user to
destination.

Figure 15 illustrate an example with GPS input. On the left of the figure we
can see the source as ”My position”, which is CSE Department and Destination
as Lucknow Airport, we can see graphical route between source and destination,
and also overall distance, duration and fuel required. In the right we have
detailed travel itinerary.

(a) Route at current location.

(b) Update with location.

Figure 16: As the current location changes the route gets updated.

Fig 16 illustrates the route update when the user moved away from the current
route. We can see the route from the current location of user (CSE Department)
as the new route gets updated.

###  VII-B Experiment on traffic trends in Indian cities

Until now, we have chosen the vehicle velocity from Table II for the
calculation of fuel consumption. The table specifies that the velocity is only
dependent on the type of road. In other words, we used a set of time-invariant
values for a vehicle’s velocity. However, in a real-world scenario, human
activities vary throughout the day, and traffic patterns vary accordingly.
Choosing a constant velocity results in either underestimating or
overestimating the fuel consumption.

Choosing an appropriate value for the velocity is an important factor for
computing fuel consumption. So we need to understand the trend of traffic
flows throughout the day and obtain a dynamic set of velocity coefficients for
the velocities for the road types from Table II. In other words, by
determining velocity coefficients it is possible to capture the dynamicity of
the traffic conditions. We approached the problem as follows in order to
search for the appropriate velocity coefficients. We generated hundreds of
source-destination pairs for different locations in five different cities
across in India. Then using Google Map API, we determined the duration of
travel for each pair of source-destination on a weekday at different times of
the day.

We compared the travel times between source-destination pairs at 9 AM, 5 PM,
and 11 PM with the travel time at 3 AM using. At 9 AM and 5 PM, the vehicular
movements reach peaks in almost all cities. At 11 PM, traffic is related to
the nightlife in cities, but the traffic is expected to be relatively low
compared to 9 AM or 5 PM. We have taken 3 AM outputs as the reference because
minimal vehicular movements is expected around the time as it coincides with
the normal cycle of human non-activity. Therefore, 3 AM is treated as zero
traffic time. The values are tabulated in Table VI, where tX represents the
average time taken to travel 1km distance at ”$X$:00” hours for
$X\in\\{3,9,17,23\\}$.

city | $\frac{t_{9}-t_{3}}{t_{3}}$ | $\frac{t_{17}-t_{3}}{t_{3}}$ | $\frac{t_{23}-t_{3}}{t_{3}}$
---|---|---|---
Bangaluru | 0.846 $\pm$ 0.232 | 0.934 $\pm$ 0.213 | 0.240 $\pm$ 0.035
Delhi | 0.698 $\pm$ 0.169 | 0.670 $\pm$ 0.160 | 0.242 $\pm$ 0.035
Hyderabad | 0.566 $\pm$ 0.153 | 0.745 $\pm$ 0.148 | 0.272 $\pm$ 0.055
Kolkata | 0.588 $\pm$ 0.140 | 0.627 $\pm$ 0.193 | 0.252 $\pm$ 0.042
Mumbai | 0.736 $\pm$ 0.278 | 0.832 $\pm$ 0.170 | 0.258 $\pm$ 0.044
TABLE VI: Additional time taken to travel 1km at different times of the day
w.r.t 3AM.

From Table VI, we find that for travel in Bangaluru city:

  * •

At 9AM on an average it takes an extra time of (0.846 $\pm$ 0.232)t3 mins to
travel 1km distance,

  * •

At 5PM on an average it takes an extra duration of (0.934 $\pm$ 0.213)t3 mins,
and

  * •

At 11 PM on average, it takes an extra duration of (0.240 $\pm$ 0.035)t3 mins.

Similarly, we can infer the factors for additional time requirements to
account for traffic in the other cities, namely Mumbai, Kolkata, Delhi and
Hyderabad.

To test the accuracy of the proposed routing engine, we have compared its
outputs for a set of sources-destination pairs with that of the Google Maps
API when traveled around 9 AM, and we computed the ratio (average time
provided by OFRM/average time provided by Google API). For Bangluru, the value
of the ratio is 0.577, and it is 0.533 for Delhi. The reason for the huge
deviation with respect to Google Maps is attributed to the static velocity
values used by OFRM without reference to variation in time of the day. Since
Google Maps collects live data, its outputs are different from our routing
engine outputs, which uses velocity values from the static Table II. It
implies that one single static velocity table for all cities and all traffic
trends cannot provide meaningful output. So, we have to choose different set
of velocity values for 9AM. Using statistics of average values from Table VI,
we find that the average velocity of the vehicle at 3AM is almost twice
compared to the velocity of the vehicle at 9AM. Using this relation we can
generate velocity Table VII for different roads for 9 AM, by taking Table II
as reference.

Road Type | velocity(km/hr)
---|---
motorway | 45
motorway_link | 22.5
trunk | 42.5
trunk_link | 20
primary | 32.5
primary_link | 15
secondary | 27.5
secondary_link | 12.5
tertiary | 20
tertiary_link | 10
unclassified | 12.5
residential | 12.5
living_street | 5
service | 7.5
TABLE VII: Velocity on various roads.

When we reevaluated our routing engine which uses velocity Table VII with the
Google Map at 9AM, the ratio (average time provided by OFRM/average time
provided by Google API) in Bangluru increased to 0.886. In Delhi, it rose to
0.941. So it is essential to generate velocity values for more time intervals
for different cities. OFRM can use these velocity values to produce results
specific to the city’s traffic trends as in Fig. 17. It clearly illustrates
that the vehicular traffic during 9 PM and 5 PM have more or less similar
trends. However, 11 AM traffic is a bit lower but not as low as 3 AM traffic.

Figure 17: Analyzing traffic in several cities at various times.

##  VIII Conclusion

We extended the Open Source Routing Machine (OSRM) to build OFRM. OFRM
provides the optimal fuel path between a given source-destination pair. It is
possible to specify any number of waypoints between a source-destination pair.
Using the proposed routing engine as the backend, one can design a pretty good
user interface for road navigating. We have built a simple web application and
mobile applications that use OSRM API in section VII. The OFRM routing
engine’s integration into a full-fledged navigation app like Google Maps or
Open Street Maps is not difficult but will take time.

In particular, some of the future work are:

  1. 1.

We will focus on improving the methods of estimating fuel consumption to
consider specific internal parameters of the car, namely vehicle radiator,
exhaust heat, wind resistance, rolling resistance, and the car’s acceleration,
instead of using a statistical model like HERA.

  2. 2.

Offline mobile application - Instead of requesting a route to an online
routing server, the routing engine should compute locally in the user’s
mobile.

  3. 3.

We can include technology to convert itinerary travel data to speech, and
also, we can improve the aesthetics of the user interface.

  4. 4.

We can add a feature, such that when a user selects an option of low in fuel,
our routing engine should be able to provide a route via gas station.

## Acknowledgements

This research is a part of the collaborative project on the eco-driving system
with Hanyang University, Seoul, Republic of South Korea. It is sponsored by a
grant from the Korean Evaluation Institute of Technology (KEIT). We are
thankful to KEIT for the grant.

## References

  * [1] European Environment Agency.  Data and maps indicators.  https://www.eea.europa.eu/data-and-maps/indicators.  Accessed on 31st July, 2020.
  * [2] EROS Center.  Usgs eros archive - digital elevation - global 30 arc-second elevation (gtopo30).  https://www.usgs.gov/.
  * [3] FHWA Resource Center.  Introduction to moves 2010.  https://www.epa.gov/sites/production/files/2016-06/documents/introduction-to-moves2010webinar.pdf, 2010\.  Accessed on 31st July, 2020 from EPA Office of Transportation and Air Quality.
  * [4] CGIAR CSI.  Srtm data.  http://srtm.csi.cgiar.org/srtmdata/.  Accessed on 31st July, 2020.
  * [5] Steve Coast.  User page in OSM.  https://www.openstreetmap.org.  Accessed on 31st July, 2020.
  * [6] Aaron J. Cohen, H. Ross Anderson, Bart Ostro, Kiran Dev Pandey, Michal Krzyzanowski, Nino Künzli, Kersten Gutschmidt, Arden Pope, Isabelle Romieu, Jonathan M. Samet, and Kirk Smith.  The global burden of disease due to outdoor air pollution.  Journal of Toxicology and Environmental Health, Part A, 68(13-14):1301–1307, 2005.
  * [7] Philip J Davis.  Interpolation and approximation.  Courier Corporation, 1975.
  * [8] Peter De Haan and Mario Keller.  Modelling fuel consumption and pollutant emissions based on real-world driving patterns: the hbefa approach.  International journal of environment and pollution, 22(3):240–258, 2004.
  * [9] Dean B Gesch.  Analysis of lidar elevation data for improved identification and delineation of lands vulnerable to sea-level rise.  Journal of Coastal Research, pages 49–58, 2009.
  * [10] A. Ghorani-Azam, B. Riahi-Zanjani, and M. Balali-Mood.  Effects of air pollution on human health and practical measures for prevention in iran.  Journal of Research in Medical Sciences, 21, 2016.
  * [11] USGS Global.  Arc-second elevation (gtopo30), 30.
  * [12] Adel Gohari, Nasir Matori, Khamaruzaman Wan Yusof, Iraj Toloue, and Kin Cho Myint.  Effects of the fuel price increase on the operating cost of freight transport vehicles.  In E3S Web of Conferences, volume 34, page 01022. EDP Sciences, 2018\.
  * [13] Michael F Goodchild.  Citizens as sensors: the world of volunteered geography.  GeoJournal, 69(4):211–221, 2007.
  * [14] Nancy B Grimm, Stanley H Faeth, Nancy E Golubiewski, Charles L Redman, Jianguo Wu, Xuemei Bai, and John M Briggs.  Global change and the ecology of cities.  science, 319(5864):756–760, 2008.
  * [15] David A. Hastings and Paula K. Dunbar.  Global land one-kilometer base elevation (globe) digital elevation model.  Technical report, National Oceanic and Atmospheric Administration, National Geophysical Data Center, Broadway, Boulder, Colorado 80303, U.S.A, 1999\.  Key to Geophysical Records Documentation (KGRD) 34.
  * [16] David A Hastings and Paula K Dunbar.  Global land one-kilometer base elevation (globe).  In Geodesy Beyond 2000, pages 101–106. Springer, 2000.
  * [17] J Hickman, D Hassel, R Joumard, Z Samaras, and S Sorenson.  Methodology for calculating transport emissions and energy consumption.  The National Academies of Sciences, Engineering, and Medicine, pages 75 – 81, 1999.
  * [18] Google Inc.  Google maps platform.  https://maps.google.com.  Accessed on 31st July, 2020.
  * [19] Andy Jarvis, Hannes Isaak Reuter, Andrew Nelson, Edward Guevara, et al.  Hole-filled srtm for the globe version 4.  available from the CGIAR-CSI SRTM 90m Database (http://srtm. csi. cgiar. org), 15:25–54, 2008.
  * [20] M Kergell.  About osmbonuspack.  https://github.com/MKergall/osmbonuspack/.  Accessed on 31st July, 2020.
  * [21] Hervé Le Treut, Ulrich Cubasch, and Myles Allen.  Historical overview of climate change science.  Notes, 16, 2005.
  * [22] Leaflet.  Api: Leaflet routing machine.  http://www.liedman.net/leaflet-routing-machine/api/.  Accessed on 31st July, 2020.
  * [23] Todd Litman.  Transportation cost and benefit analysis.  Victoria Transport Policy Institute, 31, 2009.
  * [24] Sven Maerivoet and Bart De Moor.  Traffic flow theory.  arXiv preprint physics/0507126, 2005.
  * [25] Mapbox.  Maps service.  https://docs.mapbox.com/api/maps/.  Accessed on 31st July, 2020.
  * [26] Mauro Masiol, Claudio Agostinelli, Gianni Formenton, Enzo Tarabotti, and Bruno Pavoni.  Thirteen years of air pollution hourly monitoring in a large city: Potential sources, trends, cycles and effects of car-free days.  Science of The Total Environment, 494-495:84 – 96, 2014.
  * [27] Niall McCarthy.  India has the most polluted cities on earth.  https://www.statista.com/chart/17239/average-level-of-particulate-matter-pollution/, 2019\.
  * [28] D. Nagesh.  Basic concept of remote sensing.  https://nptel.ac.in/courses/105108077/module7/lecture25.pdf, 2014\.  Accessed on 31st July, 2020.
  * [29] Leonidas Ntziachristos, Dimitrios Gkatzoflias, Chariton Kouridis, and Zissis Samaras.  Copert: a european road transport emission inventory model.  In Information technologies in environmental engineering, pages 491–504. Springer, 2009.
  * [30] Christina Nunze.  Climate 101: Air pollution.  National Geographic, 2, 2019.
  * [31] Frank Paul, Christian Huggel, and Andreas Kääb.  Combining satellite multispectral image data and a digital elevation model for mapping debris-covered glaciers.  Remote sensing of Environment, 89(4):510–518, 2004.
  * [32] Google Maps Platform.  Encoded polyline algorithm.  https://developers.google.com/maps/documentation/utilities/polylinealgorithm.  Accessed on 31st July, 2020 from Google Maps Platform.
  * [33] V. Ramanathan and Y. Feng.  Air pollution, greenhouse gases and climate change: Global and regional perspectives.  Atmospheric Environment, 43(1):37 – 50, 2009.  Atmospheric Environment - Fifty Years of Endeavour.
  * [34] Natalie Robinson, James Regetz, and Robert P Guralnick.  Earthenv-dem90: A nearly-global, void-free, multi-scale smoothed, 90m digital elevation model from fused aster and srtm data.  ISPRS Journal of Photogrammetry and Remote Sensing, 87:57–67, 2014\.
  * [35] Marc Schröder and Pedro Cabral.  Eco-friendly 3d-routing: A gis based 3d-routing-model to estimate and reduce co2-emissions of distribution transports.  Computers, Environment and Urban Systems, 73:40–55, 2019.
  * [36] Natalia Sobrino, Andres Monzon, and Sara Hernandez.  Reduced carbon and energy footprint in highway operations: the highway energy assessment (hera) methodology.  Networks and Spatial Economics, 16(1):395–414, 2016.
  * [37] I. Wagner.  India’s two-wheeler sales 2011-2018.  https://www.statista.com/statistics/318023/two-wheeler-sales-in-india/, June 2019.
  * [38] Jinghui Wang, Ahmed Elbery, and Hesham A Rakha.  A real-time vehicle-specific eco-routing model for on-board navigation applications capturing transient vehicle behavior.  Transportation Research Part C: Emerging Technologies, 104:1–21, 2019.
  * [39] Wikipedia.  Grade(slope).  https://en.wikipedia.org/wiki/Grade_(slope).  Accessed on 31st July, 2020.
  * [40] Enjian Yao and Yuanyuan Song.  Study on eco-route planning algorithm and environmental impact assessment.  Journal of Intelligent Transportation Systems, 17(1):42–53, 2013\.
  * [41] Weihua Zhang and David R Montgomery.  Digitaair pollution top 10 worst polluted cities in the world, according to iqair. l elevation model grid size, landscape representation, and hydrologic simulations.  Water resources research, 30:1019–1028, 1994.


# Synthesizing Controller for Safe Navigation using Control Density Function

Joseph Moyalan, Sriram S.K.S Narayanan, Andrew Zheng, and Umesh Vaidya
Financial support from of NSF CPS award 1932458 and NSF 2031573 is greatly
acknowledged. Joseph Moyalan, Sriram S.K.S Narayanan, Andrew Zheng, and Umesh
Vaidya are with the Department of Mechanical Engineering, Clemson University,
Clemson, SC<EMAIL_ADDRESS>
###### Abstract

We consider the problem of navigating a nonlinear dynamical system from some
initial set to some target set while avoiding collision with an unsafe set. We
extend the concept of density function to control density function (CDF) for
solving navigation problems with safety constraints. The occupancy-based
interpretation of the measure associated with the density function is
instrumental in imposing the safety constraints. The navigation problem with
safety constraints is formulated as a quadratic program (QP) using CDF. The
existing approach using the control barrier function (CBF) also formulates the
navigation problem with safety constraints as QP. One of the main advantages
of the proposed QP using CDF compared to QP formulated using CBF is that both
the convergence/stability and safety can be combined and imposed using the
CDF. Simulation results involving the Duffing oscillator and safe navigation
of Dubin car models are provided to verify the main findings of the paper.

##  I Introduction

Most control system applications in robotics and automotive engineering
include driving the nonlinear system dynamics from an initial set to a target
set while avoiding certain unsafe sets. Safe navigation is a well-known
problem in the robotics community, with applications extending to aerospace,
unmanned ground vehicles, manufacturing, power systems, etc. Existing
literature on the safe navigation problem involves jointly solving the safety
and stability problem using artificial potential methods [1]. However, using
attractive and repulsive potentials often leads to local minima, which is a
well-known problem [2, 3]. Another approach includes [4], which makes use of
the control Lyapunov function (CLF) and control barrier function (CBF) to
design a feedback law to achieve simultaneously both convergence to a target
set and avoidance of an unsafe set. However, the control design is problem-
specific and less intuitive, and if the safety and convergence objectives
conflict, then such a feedback law can’t be designed. Another approach uses
the combination of the CLF and CBF to construct a quadratic program (QP) to
solve for the desired safe control [5, 6]. However, finding CLF and CBF to
formulate the QP is not trivial.
The solution to the navigation problem is also achieved in the dual-density
space. The density function was first introduced in [7] as a dual to Lyapunov
function for stability analysis. This density function was later used as a
safety certificate using the sum of squares optimization method for the
analysis and safety of the nonlinear system [8]. In [9], the authors utilized
navigation measures for designing safe controllers using a convex formulation.
In [10], the authors formulate the safety problem as a co-design problem of
finding the density function and optimal safe controller. However, the
iterative approach of estimating the density function from the optimal control
can be computationally expensive, depending upon the complexity of the
problem. In contrast, we design a controller for a given known density
function, which is constructed based on the occupancy-based interpretation of
the density. Similarly, the use of linear operators such as Koopman and
Perron-Frobenius (P-F) operators for convex data-driven approaches for optimal
control and control with safety constraints have been explored in [11, 12, 13,
14]. The convex data-driven approaches heavily rely on observable functions to
lift the dynamics to the space of functions. The number of observables
required to lift the dynamics to function space increases with an increase in
the dimension or complexity of the underlying nonlinear system, which makes it
computationally expensive. As such, these approaches suffer the curse of
dimensionality. In [15], to avoid the curse of dimensionality, the authors
provide an analytical expression for the density-based feedback controller for
the navigation of single integrator dynamics, which can be viewed as a dual
construction of the navigation functions from [3].
The control design in [15] is only applicable to single integrator dynamics or
systems, which can be reformulated as a single integrator using a change of
coordinates and inverse dynamics approach. This paper extends the density-
based controller introduced in [15] to nonlinear systems with drift.
Therefore, unlike [15], our approach can be applied to any general nonlinear
system that contains drift. In this paper, we introduce the notion of control
density function (CDF). The CDF is an extension of the density function, just
like CLF is the extension of the Lyapunov function for control systems. The
construction of CDF is provided based on [15], which is then utilized to
formulate the problem statement as a QP with CDF-based constraints. The CDF-
based constraint ensures a nonlinear control system’s simultaneous convergence
and safety. This contrasts with the approaches given in [5, 6] where one needs
to augment CLF constraints with CBF to ensure convergence. Another difference
includes almost everywhere convergence of the system dynamics, which is a
weaker notion of convergence than control Lyapunov functions. As a result, the
control law obtained from CDF constraints ensures convergence for almost all
initial conditions. The set of all initial conditions not converging to the
target set will have a zero Lebesgue measure [16]. We also provide an example
of an underactuated system in the form of the Dubin car model, where obstacles
are only present in the subspace of the system dynamics. Finally, unlike [15],
we also show that the density-based safe controller can be combined with a
nominal controller for optimal performance by modifying the cost function to
include the norm of the difference between the desired safe control and the
nominal control. The rest of the paper is organized as follows. Section II
contains the problem statement and some preliminaries, and Section III
consists of the paper’s main results. Section IV discusses the computational
framework for the prescribed QP. In Section V, we provide some simulation
results followed by some conclusions.

##  II Preliminaries and Problem Statement

Consider the dynamical system of the form

$\displaystyle\dot{{\mathbf{x}}}={\mathbf{f}}({\mathbf{x}})+{\mathbf{g}}({\mathbf{x}}){\mathbf{u}}$
(1)

where ${\mathbf{x}}\in D\subseteq\mathds{R}^{n}$ and ${\mathbf{u}}\in
U\subseteq\mathds{R}^{m}$ are the states and the control inputs respectively.
We also assume that
${\mathbf{f}},{\mathbf{g}}\in\mathcal{C}^{1}(D,\mathds{R}^{n})$ are
continuously differentiable functions on $D$. The unsafe sets are represented
by $\mathcal{U}\subset D$. Next, ${\mathbf{X}}_{0},{\mathbf{X}}_{T}\subset D$
represents the initial set and target set, respectively. In this paper, we
assume the target is located at the origin, i.e., ${\mathbf{X}}_{T}={0}$. We
represent $s_{t}({\mathbf{x}})$ to be the solution of (1) with respect to some
control ${\mathbf{u}}$ at time $t$ starting from the initial condition
${\mathbf{x}}$. We also denote $\mathcal{B}(D)$ to be the Borel
$\sigma$-algebra on $D$ and ${\cal M}(D)$ as the vector space of real-valued
measures on ${\cal B}(D)$ and $m(\cdot)$ denotes Lebesgue measure. Also, we
represent $\bar{D}:=D\setminus\mathcal{N}_{\eta}$ where $\mathcal{N}_{\eta}$
represents a small neighborhood of $\eta$ radius around ${\mathbf{X}}_{T}$.
The notation $\nabla_{{\mathbf{x}}}$ denotes $[\frac{\partial}{\partial
x_{1}},\dots,\frac{\partial}{\partial x_{n}}]^{\top}$ where
${\mathbf{x}}\in\mathds{R}^{n}$. Also, $tr(A)$ represents the trace of the
matrix $A$.

######  Problem 1.

(Almost everywhere (a.e.) safe navigation) The primary objective of this paper
is to design a control ${\mathbf{u}}$ to drive the trajectories of the system
given by (1) from almost all initial conditions (w.r.t. Lebesgue measure) from
the initial set ${\mathbf{X}}_{0}$ to a target set ${\mathbf{X}}_{T}$ while
avoiding the unsafe set $\mathcal{U}$.

###  II-A Density function for safe navigation

The construction of the density-based safe control for the dynamics given in
(1) is inspired by [15]. We define the unsafe set $\mathcal{U}$. Let there be
$K$ number of obstacles. We define a continuous scalar-valued function
$c_{k}({\mathbf{x}})$ for $k=1,\dots,K$. Now, each obstacle set can be defined
as follows:

$\displaystyle C_{k}:=\\{{\mathbf{x}}\in D:c_{k}({\mathbf{x}})\leq 0\\}$ (2)

Therefore, the set defining the total unsafe region is given by
$\mathcal{U}:=\bigcup_{k=1}^{K}\;C_{k}$. Similarly, we will define another
continuous scalar-valued function $b_{k}({\mathbf{x}})$ to define the sensing
region for each obstacle set as follows:

$\displaystyle B_{k}:=\\{{\mathbf{x}}\in D:b_{k}({\mathbf{x}})\leq
0\\}\setminus C_{k}$ (3)

Now, we use the functions $c_{k}({\mathbf{x}})$ and $b_{k}({\mathbf{x}})$ to
formulate a smooth inverse bump function. First, we start by constructing the
following functions,

$\displaystyle m_{k}$
$\displaystyle=\frac{c_{k}({\mathbf{x}})}{c_{k}({\mathbf{x}})-b_{k}({\mathbf{x}})},$
$\displaystyle\psi_{k}({\mathbf{x}})$
$\displaystyle=\frac{\exp(\frac{-1}{m_{k}})}{\exp(\frac{-1}{m_{k}})+\exp(\frac{-1}{1-m_{k}})}.$

Using the functions $m_{k}({\mathbf{x}})$ and $\psi_{k}({\mathbf{x}})$, we
define a piece-wise smooth inverse bump function $\Psi_{k}({\mathbf{x}})$ as
follows:

$\displaystyle\Psi_{k}({\mathbf{x}})=\begin{cases}0,&{\mathbf{x}}\in C_{k}\\\
\psi_{k}({\mathbf{x}}),&{\mathbf{x}}\in B_{k}\\\
1,&\text{otherwise}\end{cases}.$ (4)

Note that $\prod_{k=1}^{K}{\mathbf{\Psi}}_{k}({\mathbf{x}})$ encodes the
unsafe set $\mathcal{U}$. To encode information about the target set
${\mathbf{X}}_{T}$, we use $V({\mathbf{x}})$ which acts as a distance function
from the current state ${\mathbf{x}}$ to the target state ${\mathbf{X}}_{T}$.
The $V({\mathbf{x}})$ can be modified to adjust to the geometry of the
underlying configuration space. In the Euclidean space with
${\mathbf{x}}\in\mathds{R}^{n}$, we use
$V({\mathbf{x}})=({\mathbf{x}}-{\mathbf{X}}_{T})^{\top}P({\mathbf{x}}-{\mathbf{X}}_{T})$
for some $P>0$. Finally, the density function $\rho({\mathbf{x}})$ for the
safe control is given as follows:

$\displaystyle\rho({\mathbf{x}})=\frac{\prod_{k=1}^{K}{\mathbf{\Psi}}_{k}({\mathbf{x}})}{V({\mathbf{x}})^{\alpha}}$
(5)

for some $\alpha>0$.

It can be observed that the density function $\rho({\mathbf{x}})$ given in (5)
is a smooth function for all ${\mathbf{x}}\in D$. Fig. 1a shows an environment
with one obstacle set $\mathcal{U}$ and target ${\mathbf{X}}_{T}$. Fig. 1b
shows the corresponding density function representation. Note that the density
function $\rho({\mathbf{x}})$ takes minimum value for
${\mathbf{x}}\in\mathcal{U}$ and max value for
${\mathbf{x}}\in{\mathbf{X}}_{T}$.

Figure 1: (a) Environment setup with unsafe set $\mathcal{U}$ and target
${\mathbf{X}}_{T}$, (b) density function $\rho({\mathbf{x}})$ for this
environment.

######  Assumption 1.

We assume that system dynamics given by (1) is locally linearizable inside the
region $\mathcal{N}_{\eta}$ containing ${\mathbf{X}}_{T}=0$. We also assume
that the linearized dynamics at the origin are stabilizable. Therefore the
local stability of ${\mathbf{X}}_{T}=0$ can be achieved by designing a local
linear feedback controller.

The construction of density function in (5) was proposed in [15] and used for
solving the safe navigation problem with simple integrator dynamics of the
form $\dot{\mathbf{x}}={\mathbf{u}}$. The feedback controller ${\mathbf{u}}$
for the safe navigation was shown to be of the form

$\displaystyle{\mathbf{u}}({\mathbf{x}})=\frac{\partial\rho}{\partial{\mathbf{x}}}$
(6)

This paper’s main contribution is introducing the CDF for solving the
navigation problem as stated in Problem 1 for nonlinear system (1) with drift.

##  III Control Density Function (CDF)

The density functions can be utilized to check the stability of any nonlinear
system [17, 14]. However, such density functions cannot be used to design a
safe control. By drawing inspiration from how Lyapunov functions were extended
to CLF [18, 19], and how barrier functions were extended to CBF [20, 6], we
propose in this section an extension of the density function and call it as
the CDF. Let us again consider the nonlinear system as follows:

$\displaystyle\dot{{\mathbf{x}}}={\mathbf{f}}({\mathbf{x}})+{\mathbf{g}}({\mathbf{x}}){\mathbf{u}}$
(7)

where ${\mathbf{u}}=[u_{1},\dots,u_{m}]^{\top}$ and
${\mathbf{g}}({\mathbf{x}})=[g_{1}({\mathbf{x}}),\dots,g_{m}({\mathbf{x}})]$.
All the vector fields are assumed to belong to
$\mathcal{C}^{1}(\bar{D},\mathds{R}^{n})$.

The following theorem provides the main result for the safe navigation control
for a nonlinear system using CDF.

######  Theorem 1.

Under Assumption 1, given system dynamics in (1) and density function given by
(5), the system trajectories can be driven from almost all initial conditions
to a target set ${\mathbf{X}}_{T}$ while avoiding unsafe set $\mathcal{U}$ if
there exists a control ${\mathbf{u}}\in U$ and $\lambda>0$ such that

$\displaystyle\nabla\cdot({\mathbf{f}}({\mathbf{x}})\rho+{\mathbf{g}}({\mathbf{x}}){\mathbf{u}}\rho)\geq
0,\;\;a.e.\;{\mathbf{x}}\in\bar{D}$
$\displaystyle\nabla\cdot({\mathbf{f}}({\mathbf{x}})\rho+{\mathbf{g}}({\mathbf{x}}){\mathbf{u}}\rho)\geq\lambda>0,\;\;\forall\;{\mathbf{x}}\in{\mathbf{X}}_{0}$
(8)

The proof is provided in the Appendix.

Based on the result of Theorem 1, we see that solving for the desired safe
control ${\mathbf{u}}$ is an infinite-dimensional feasibility problem.
However, if we assume that ${\mathbf{f}}({\mathbf{x}})$,
${\mathbf{g}}({\mathbf{x}})$ and $\rho({\mathbf{x}})$ are known to us, then we
can solve for ${\mathbf{u}}$ point-wise along the system trajectory by
formulating it as QP utilizing the CDF constraints given in (8). The QP-CDF
formulation is given below:

$\displaystyle\;\;\;\;\;\;\;\;\min_{{\mathbf{u}}}\;\|{\mathbf{u}}\|^{2}$ s.t.
$\displaystyle\nabla\cdot({\mathbf{f}}({\mathbf{x}})\rho+{\mathbf{g}}({\mathbf{x}}){\mathbf{u}}\rho)\geq
0,\;\;a.e.\;{\mathbf{x}}\in\bar{D}$
$\displaystyle\nabla\cdot({\mathbf{f}}({\mathbf{x}})\rho+{\mathbf{g}}({\mathbf{x}}){\mathbf{u}}\rho)\geq\lambda>0,\;\;\forall\;{\mathbf{x}}\in{\mathbf{X}}_{0}$
(9)

######  Remark 1.

If there exists a nominal control, ${\mathbf{u}}_{0}$, for the system given by
(7), we can reformulate (9) to enforce ${\mathbf{u}}_{0}$ in the absence of
unsafe sets by introducing a cost function which minimizes
$\|{\mathbf{u}}-{\mathbf{u}}_{0}\|^{2}$.

##  IV Computational framework

The inequality in (8) can be expanded as follows:

$\displaystyle\nabla\cdot({\mathbf{f}}({\mathbf{x}})\rho)+\nabla_{{\mathbf{x}}}^{\top}\left[{\mathbf{g}}({\mathbf{x}})\rho\right]{\mathbf{u}}+tr\left((\nabla_{{\mathbf{x}}}{\mathbf{u}}^{\top})^{\top}{\mathbf{g}}\rho\right)>0$
(10)

We observe that while trying to solve for ${\mathbf{u}}$ point-wise along the
system trajectory, the term $\nabla u_{i}$ is a spatial operator and requires
information on points in the neighborhood of the system trajectory. Therefore,
to calculate $\nabla u_{i}$, we perturb the point on the trajectory by
$\epsilon$ along the directional basis to obtain points $[z_{1},\dots,z_{n}]$
where $z_{i}\in\mathds{R}^{n}$. Here,
$z_{i}={\mathbf{x}}+\epsilon\;\textbf{e}_{i}$ where $\textbf{e}_{i}$ is a
column vector consisting of all zeros except at $i^{th}$ position where the
value is $1$. These points around the trajectory will be used to calculate the
$\nabla u_{i}$.

To make things easier, we can split the inequality given in (10) as follows:

$\displaystyle\nabla\cdot({\mathbf{f}}\rho)+\nabla_{{\mathbf{x}}}^{\top}\left[{\mathbf{g}}({\mathbf{x}})\rho\right]{\mathbf{u}}\geq\beta\rho$
$\displaystyle\left|tr\left((\nabla_{{\mathbf{x}}}{\mathbf{u}}^{\top})^{\top}{\mathbf{g}}\right)\right|<\beta$
(11)

Now we can write $\nabla u_{i}$ as follows:

$\displaystyle\nabla
u_{i}=\left[\frac{u_{i}^{1}-u_{i}}{\epsilon},\dots,\frac{u_{i}^{n}-u_{i}}{\epsilon}\right]^{\top}$
(12)

where ${\mathbf{u}}^{j}=[u^{j}_{1},\dots,u^{j}_{m}]^{\top}$ represents the
control satisfying (11) at point $z_{j}$. This can be written in matrix form
for all the $m$ control values as follows:

$\displaystyle(\nabla_{{\mathbf{x}}}{\mathbf{u}}^{\top})^{\top}=\begin{bmatrix}\frac{u_{1}^{1}-u_{1}}{\epsilon}&\dots&\frac{u_{1}^{n}-u_{1}}{\epsilon}\\\
&\ddots&\\\
\frac{u_{m}^{1}-u_{m}}{\epsilon}&\dots&\frac{u_{m}^{n}-u_{m}}{\epsilon}\end{bmatrix}$
(13)

Therefore, we can rewrite (11) as follows:

$\displaystyle\nabla\cdot({\mathbf{f}}({\mathbf{x}})\rho({\mathbf{x}}))+\nabla_{{\mathbf{x}}}^{\top}\left[{\mathbf{g}}({\mathbf{x}})\rho({\mathbf{x}})\right]{\mathbf{u}}\geq\beta\rho({\mathbf{x}})$
$\displaystyle\nabla\cdot({\mathbf{f}}(z_{1})\rho(z_{1}))+\nabla_{{\mathbf{x}}}^{\top}\left[{\mathbf{g}}(z_{1})\rho(z_{1})\right]{\mathbf{u}}^{1}\geq\beta\rho(z_{1})$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\vdots$
$\displaystyle\nabla\cdot({\mathbf{f}}(z_{n})\rho(z_{n}))+\nabla_{{\mathbf{x}}}^{\top}\left[{\mathbf{g}}(z_{n})\rho(z_{n})\right]{\mathbf{u}}^{n}\geq\beta\rho(z_{n})$
$\displaystyle\left|tr((\nabla_{{\mathbf{x}}}{\mathbf{u}}^{\top})^{\top}\;{\mathbf{g}}({\mathbf{x}}))\right|<\beta$
(14)

Here, the number of decision variables to solve the above linear inequalities
will be $m(n+1)$. The algorithm 1 summarizes the steps to solve for safe
navigation control using QP-CDF as given in (9).

######  Remark 2.

The value of $\alpha$ determines the rate of convergence towards the target
set. This stems from the fact that $\alpha$ appears in the distance function
$V({\mathbf{x}})$ which contains information on the target set. It can be
inferred that the convergence rate can be increased by increasing the value of
$\alpha$ and vice-versa. Similarly, increasing the sensing region around the
obstacle will lead to smoother avoidance control, and decreasing the sensing
region will lead to more aggressive avoidance control near the obstacle.

Input: ${\mathbf{f}},{\mathbf{g}},\rho,{\mathbf{x}}_{0},\beta,N$

for  _$k=1:N$_ do

$z_{j}={\mathbf{x}}_{k-1}+\textbf{e}_{j}\;\;\forall\;j=1,\dots,n$.

Solve for
_${\mathbf{u}}_{k},{\mathbf{u}}^{1}_{k},\dots,{\mathbf{u}}^{n}_{k}$_ do

$\min\;\|{\mathbf{u}}_{k}\|^{2}+\|{\mathbf{u}}^{1}_{k}\|^{2}+\dots+\|{\mathbf{u}}^{n}_{k}\|^{2}$

s.t.

$\nabla\cdot({\mathbf{f}}({\mathbf{x}}_{k-1})\rho({\mathbf{x}}_{k-1}))+$

$\nabla_{{\mathbf{x}}}^{\top}\left[{\mathbf{g}}({\mathbf{x}}_{k-1})\rho({\mathbf{x}}_{k-1})\right]{\mathbf{u}}_{k}>\beta\rho({\mathbf{x}}_{k-1})$,

$\nabla\cdot({\mathbf{f}}(z_{1})\rho(z_{1})+\nabla_{{\mathbf{x}}}^{\top}\left[{\mathbf{g}}(z_{1})\rho(z_{1})\right]{\mathbf{u}}^{1}_{k}>\beta\rho(z_{1})$,

$\;\;\;\;\;\;\;\vdots$

$\nabla\cdot({\mathbf{f}}(z_{n})\rho(z_{n})+\nabla_{{\mathbf{x}}}^{\top}\left[{\mathbf{g}}(z_{n})\rho(z_{n})\right]{\mathbf{u}}^{n}_{k}>\beta\rho(z_{n})$,

$\left|tr((\nabla_{{\mathbf{x}}}{\mathbf{u}}_{k}^{\top})^{\top}\;{\mathbf{g}}({\mathbf{x}}_{k-1}))\right|<\beta$

${\mathbf{x}}_{k}={\mathbf{x}}_{k-1}+\Delta
t({\mathbf{f}}({\mathbf{x}}_{k-1})+{\mathbf{g}}({\mathbf{x}}_{k-1}){\mathbf{u}}_{k})$

Algorithm 1 QP-CDF

##  V Simulation Results

In this section, we provide some navigation results for different types of
system dynamics starting with the Duffing oscillator.

###  V-A Duffing oscillator

Let us consider the Duffing oscillator dynamics:

$\displaystyle\dot{x}_{1}=x_{2}$
$\displaystyle\dot{x}_{2}=x_{1}-x_{1}^{3}-0.1x_{2}+u$ (15)

Figure 2: (a) Trajectories converging to the target (green) while avoiding the
unsafe set (gray), (b) Control plots for the controlled Duffing oscillator
dynamics obtained by solving QP-CDF.

Fig. 2 provides a safe navigation control trajectory for the Duffing
oscillator by solving QP-CDF. The density function is constructed based on
equations (4) and (5). The functions used to describe the obstacle set are
given by $c({\mathbf{x}}):=\|{\mathbf{x}}-o\|^{2}-r_{1}^{2}$ where $o$ and
$r_{1}$ are the center and radius of the circular obstacle respectively.
Similarly, the function used to define sensing region is given by
$b({\mathbf{x}}):=\|{\mathbf{x}}-o\|^{2}-r_{2}^{2}$ where $r_{2}$ is the
radius of the sensing region of the circular obstacle. Here, the obstacle
radius is 0.5 units and the sensing radius is 0.7 units with its center
located at $[0,0]$. The matrix $P$ used to construct the function
$V({\mathbf{x}})$ is obtained by linearizing the dynamics around the target
point and solving the algebraic Riccati equation with the identity matrix as
the state and control gain matrix. Here, the control bounds are $\pm 2$. It
can be observed that the control action obtained from the QP-CDF is minimally
invasive. The system trajectories follow along the vector field of the Duffing
oscillator when it is away from the obstacle. The control from the QP-CDF is
mainly active near the obstacles by making the system trajectories drive away
from the obstacles while trying to converge to the target set.

###  V-B Dubin car model

Let us consider the Dubin car model as follows:

$\displaystyle\dot{x}_{1}$ $\displaystyle=v\cos(\theta)$ (16a)
$\displaystyle\dot{x}_{2}$ $\displaystyle=v\sin(\theta)$ (16b)
$\displaystyle\dot{\theta}$ $\displaystyle=\omega$ (16c)

where $u$ and $\omega$ are the control inputs. We assume that the obstacles to
be avoided are only present in the $x_{1}-x_{2}$ space. In this example, we
first reformulate (16a)-(16b) in the form of single integrator dynamics as
given in (17) and solve for $u_{1}$ and $u_{2}$ using the QP-CDF given in (9).
These control values will be then utilized in the design of $v$ and $\omega$.
The single integrator dynamics in the $x_{1}-x_{2}$ space is as follows:

$\displaystyle\dot{x}_{1}$ $\displaystyle=u_{1}$ (17a)
$\displaystyle\dot{x}_{2}$ $\displaystyle=u_{2}$ (17b)

The safe control for (17) is obtained by solving the QP-CDF given in (9).
Here, $c_{k}({\mathbf{x}}):=\|{\mathbf{x}}-o_{k}\|^{2}-r_{1k}^{2}$ and
$b_{k}({\mathbf{x}}):=\|{\mathbf{x}}-o_{k}\|^{2}-r_{2k}^{2}$ where $o_{k}$,
$r_{1k}$ and $r_{2k}$ are the center, radius and sensing radius of the
$k^{th}$ circular obstacle. The control obtained as the solution of the QP-CDF
can be used to calculate $u$ and $\tilde{\theta}$ as follows:

$\displaystyle
v=\sqrt{u_{1}^{2}+u_{2}^{2}},\;\;\;\tilde{\theta}=\tan^{-1}\left(\frac{u_{2}}{u_{1}}\right)$
(18)

Now, the control $\omega$ needs to be designed such that
$\theta-\tilde{\theta}$ tends to zero. Therefore, we consider the Lyapunov
function given by $0.5(\theta-\tilde{\theta})^{2}$ which gives us the
following control law for $\omega$:

$\displaystyle\omega=\dot{\tilde{\theta}}-k(\theta-\tilde{\theta})$ (19)

for some $k>0$. Fig. 3 provides a safe navigation trajectory for the unicycle
model in $x_{1}-x_{2}$ space. The obstacle radius is 2 units and the sensing
radius is 2.5 units with the center for the two obstacles located at $[3,1]$
and $[7.5,-1]$ respectively. The value of gain $k$ is chosen to be 10. The
control action mainly drives the system trajectories away from the obstacle
set while converging on the target set.

Figure 3: (a) Trajectories in the $x_{1}-x_{2}$ space converging to the target
(green) while avoiding unsafe sets (gray), (b) state trajectory plot of
$\theta$ for safe navigation.

##  VI Conclusion

The problem of safe navigation of nonlinear dynamical systems is considered.
We introduce the notion of CDF as an extension to the density-based
controllers. The CDF-based safety constraints are inspired by the occupancy-
based interpretation of the measure associated with the density function. We
then formulate the navigation problem with safety constraints as a QP using
CDF. The existing approach for navigation problems includes CBF-based QP. The
advantage of CDF-based QP over CBF-based QP is that both the convergence and
safety constraints can be imposed using CDF. Finally, we provide simulation
results for safe navigation of the Dubin car model and Duffing oscillator to
showcase the validity of CDF-based QP.

##  VII Appendix

The proof of Theorem 1 relies on the following Lemma.

######  Lemma 1.

If

$\displaystyle\int_{0}^{\infty}\int_{{\mathbf{X}}_{0}}\mathds{1}_{\mathcal{U}}(s_{t}({\mathbf{x}}))d{\mathbf{x}}\;dt=0$
(20)

then

$\displaystyle\int_{{\mathbf{X}}_{0}}\mathds{1}_{\mathcal{U}}(s_{t}({\mathbf{x}}))d{\mathbf{x}}=0\;\;\;\forall
t\geq 0$ (21)

i.e., the amount of time system trajectories spend in set $\mathcal{U}$
starting from the positive measure set ${\mathbf{X}}_{0}$ is equal to zero.

Proof: The proof is done by the method of contradiction. Let us assume (21) is
not true. Then there exists a time $\bar{t}$ such that

$\int_{{\mathbf{X}}_{0}}\mathds{1}_{\mathcal{U}}(s_{\bar{t}}({\mathbf{x}}))d{\mathbf{x}}>0.$

Now, from the continuity of the solution of the differential equation, there
exists a $\Delta$ such that

$\int_{\bar{t}}^{\bar{t}+\Delta}\int_{{\mathbf{X}}_{0}}\mathds{1}_{\mathcal{U}}({\mathbf{s}}_{t}({\mathbf{x}}))d{\mathbf{x}}\;dt>0$

Therefore,

$\displaystyle 0<$
$\displaystyle\int_{\bar{t}}^{\bar{t}+\Delta}\int_{{\mathbf{X}}_{0}}\mathds{1}_{\mathcal{U}}({\mathbf{s}}_{t}({\mathbf{x}}))d{\mathbf{x}}\;dt$
$\displaystyle\leq$
$\displaystyle\int_{0}^{\infty}\int_{{\mathbf{X}}_{0}}\mathds{1}_{\mathcal{U}}(s_{t}({\mathbf{x}}))d{\mathbf{x}}\;dt=0.$

Hence, we arrive at a contradiction.

Proof of Theorem 1: Let us consider the following:

$\displaystyle\nabla\cdot({\mathbf{f}}({\mathbf{x}})\rho+{\mathbf{g}}({\mathbf{x}}){\mathbf{u}}\rho)=h({\mathbf{x}})$
(22)

where $h({\mathbf{x}})\geq 0$ and $h({\mathbf{x}})\geq\lambda>0$ for
${\mathbf{x}}\in{\mathbf{X}}_{0}$. Now, through the method of characteristics,
the function $\rho({\mathbf{x}})$ can be written as follows:

$\displaystyle\rho({\mathbf{x}})=\int^{\infty}_{0}h(s_{-t}({\mathbf{x}}))\left|\frac{\partial
s_{-t}({\mathbf{x}})}{\partial{\mathbf{x}}}\right|dt,$ (23)

where $|\cdot|$ represents the determinant. This can be easily verified by
simple substitution of (23) in (22) and using the fact that

$\displaystyle\lim_{t\rightarrow\infty}h(s_{-t}({\mathbf{x}}))\left|\frac{\partial
s_{-t}({\mathbf{x}})}{\partial{\mathbf{x}}}\right|=0.$ (24)

The limit in (24) is the consequence of $\rho({\mathbf{x}})$ being bounded in
$\bar{D}$ and using Barbalat’s Lemma. The term inside the integral in (23) can
be written using the linear Perron-Frobenius (P-F) operator as follows:

$\displaystyle[\mathbb{P}_{t}h]({\mathbf{x}})=h(s_{-t}({\mathbf{x}}))\left|\frac{\partial
s_{-t}({\mathbf{x}})}{\partial{\mathbf{x}}}\right|.$ (25)

Therefore,

$\displaystyle\rho({\mathbf{x}})=\int^{\infty}_{0}[\mathbb{P}_{t}h]({\mathbf{x}})dt.$
(26)

Now, utilizing (24), we can write

$\displaystyle\lim_{t\rightarrow\infty}[\mathbb{P}_{t}h]({\mathbf{x}})=0\implies\lim_{t\rightarrow\infty}[\mathbb{P}_{t}\mathds{1}_{{\mathbf{X}}_{0}}]({\mathbf{x}})=0$
(27)

which follows because
$h({\mathbf{x}})\geq\lambda>0\;\;\forall{\mathbf{x}}\in{\mathbf{X}}_{0}$ and
using dominated convergence theorem. Here, $\mathds{1}_{{\mathbf{X}}_{0}}$
represents the indicator function for ${\mathbf{X}}_{0}$. Now, for any
$A\subseteq\bar{D}$, we have

$\displaystyle\int_{A}[\mathbb{P}_{t}\mathds{1}_{{\mathbf{X}}_{0}}]({\mathbf{x}})d{\mathbf{x}}=$
$\displaystyle\int_{\bar{D}}[\mathbb{P}_{t}\mathds{1}_{{\mathbf{X}}_{0}}]({\mathbf{x}})\mathds{1}_{A}({\mathbf{x}})d{\mathbf{x}}$
$\displaystyle=$
$\displaystyle\int_{\bar{D}}\mathds{1}_{{\mathbf{X}}_{0}}({\mathbf{x}})\mathds{1}_{A}(s_{t}({\mathbf{x}}))d{\mathbf{x}}$
(28)

This can be observed by using the definition of P-F operator in (25) and doing
the change of variables in integration such as
${\mathbf{y}}=s_{-t}({\mathbf{x}})$ and $d{\mathbf{y}}=\left|\frac{\partial
s_{-t}({\mathbf{x}})}{\partial{\mathbf{x}}}\right|d{\mathbf{x}}$ and
relabelling. The right hand side of (28) can be seen as follows:

$\displaystyle\int_{A}[\mathbb{P}_{t}\mathds{1}_{{\mathbf{X}}_{0}}]({\mathbf{x}})d{\mathbf{x}}=m\\{{\mathbf{x}}\in{\mathbf{X}}_{0}:s_{t}({\mathbf{x}})\in
A\\}.$

Therefore, using (27), we observe that

$\displaystyle
0=\int_{A}[\lim_{t\rightarrow\infty}\mathbb{P}_{t}\mathds{1}_{{\mathbf{X}}_{0}}]({\mathbf{x}})d{\mathbf{x}}=m\\{{\mathbf{x}}\in{\mathbf{X}}_{0}:\lim_{t\rightarrow\infty}s_{t}({\mathbf{x}})\in
A\\}.$

The above statement can be generalized for any measurable Lebesgue set
$A\subseteq\bar{D}$. Therefore,

$\displaystyle
m\\{{\mathbf{x}}\in{\mathbf{X}}_{0}:\lim_{t\rightarrow\infty}s_{t}({\mathbf{x}})\neq
0\\}=0.$

Now, from the construction of the density function, we know that
$\rho({\mathbf{x}})=0\;\forall\;{\mathbf{x}}\in\mathcal{U}$. Therefore,

$\displaystyle\int_{\mathcal{U}}\int^{\infty}_{0}[\mathbb{P}_{t}\mathds{1}_{{\mathbf{X}}_{0}}]({\mathbf{x}})dtd{\mathbf{x}}\leq\int_{\mathcal{U}}\rho({\mathbf{x}})d{\mathbf{x}}=0.$
(29)

Utilizing the Markov property of the P-F operator and the fact that indicator
functions are non-negative functions, we can rewrite (29) as follows:

$\displaystyle\int_{\mathcal{U}}\int^{\infty}_{0}[\mathbb{P}_{t}\mathds{1}_{{\mathbf{X}}_{0}}]({\mathbf{x}})dtd{\mathbf{x}}=0.$
(30)

Now, doing the change of variables in integration such as
${\mathbf{y}}=s_{-t}({\mathbf{x}})$ and $d{\mathbf{y}}=\left|\frac{\partial
s_{-t}({\mathbf{x}})}{\partial{\mathbf{x}}}\right|d{\mathbf{x}}$ and
relabelling, the left-hand side of the (30) can be written as follows:

$\displaystyle\int_{\mathcal{U}}\int^{\infty}_{0}[\mathbb{P}_{t}\mathds{1}_{{\mathbf{X}}_{0}}]$
$\displaystyle({\mathbf{x}})dtd{\mathbf{x}}=\int_{\bar{D}}\int^{\infty}_{0}[\mathbb{P}_{t}\mathds{1}_{{\mathbf{X}}_{0}}]({\mathbf{x}})\mathds{1}_{\mathcal{U}}({\mathbf{x}})dtd{\mathbf{x}}$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\int_{\bar{D}}\mathds{1}_{{\mathbf{X}}_{0}}({\mathbf{x}})\mathds{1}_{\mathcal{U}}(s_{t}({\mathbf{x}}))d{\mathbf{x}}\;dt$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\int_{{\mathbf{X}}_{0}}\mathds{1}_{\mathcal{U}}(s_{t}({\mathbf{x}}))d{\mathbf{x}}\;dt=0.$
(31)

Now, from (31) and using Lemma 1, we can conclude the following:

$\displaystyle\int_{{\mathbf{X}}_{0}}\mathds{1}_{\mathcal{U}}(s_{t}({\mathbf{x}}))d{\mathbf{x}}=0\;\;\;\forall
t\geq 0$

## References

  * [1] O. Khatib, “Real-time obstacle avoidance for manipulators and mobile robots,” _The international journal of robotics research_ , vol. 5, no. 1, pp. 90–98, 1986.
  * [2] B. Krogh, “A generalized potential field approach to obstacle avoidance control,” in _Proc. SME Conf. on Robotics Research: The Next Five Years and Beyond, Bethlehem, PA, 1984_ , 1984, pp. 11–22.
  * [3] E. Rimon, _Exact robot navigation using artificial potential functions_. Yale University, 1990.
  * [4] M. Z. Romdlony and B. Jayawardhana, “Uniting control lyapunov and control barrier functions,” in _53rd IEEE Conference on Decision and Control_. IEEE, 2014, pp. 2293–2298.
  * [5] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs for safety critical systems,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 8, pp. 3861–3876, 2016.
  * [6] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control barrier functions: Theory and applications,” in _2019 18th European control conference (ECC)_. IEEE, 2019, pp. 3420–3431.
  * [7] A. Rantzer, “A dual to lyapunov’s stability theorem,” _Systems & Control Letters_, vol. 42, no. 3, pp. 161–168, 2001.
  * [8] A. Rantzer and S. Prajna, “On analysis and synthesis of safe control laws,” in _42nd Allerton Conference on Communication, Control, and Computing_. University of Illinois, 2004, pp. 1468–1476.
  * [9] U. Vaidya, “Optimal motion planning using navigation measure,” _International Journal of Control_ , vol. 91, no. 5, pp. 989–998, 2018.
  * [10] Y. Chen, M. Ahmadi, and A. D. Ames, “Optimal safe controller synthesis: A density function approach,” in _2020 American Control Conference (ACC)_. IEEE, 2020, pp. 5407–5412.
  * [11] X. Ma, B. Huang, and U. Vaidya, “Optimal quadratic regulation of nonlinear system using koopman operator,” in _2019 American Control Conference (ACC)_. IEEE, 2019, pp. 4911–4916.
  * [12] B. Huang and U. Vaidya, “A convex approach to data-driven optimal control via perron-frobenius and koopman operators,” _IEEE Transactions on Automatic Control_ , 2022.
  * [13] H. Yu, J. Moyalan, U. Vaidya, and Y. Chen, “Data-driven optimal control of nonlinear dynamics under safety constraints,” _IEEE Control Systems Letters_ , vol. 6, pp. 2240–2245, 2022.
  * [14] J. Moyalan, H. Choi, Y. Chen, and U. Vaidya, “Data-driven optimal control via linear transfer operators: A convex approach,” _Automatica_ , vol. 150, p. 110841, 2023.
  * [15] A. Zheng, S. S. Narayanan, and U. Vaidya, “Safe navigation using density functions,” _IEEE Robotics and Automation Letters_ , 2023.
  * [16] J. Moyalan, H. Choi, Y. Chen, and U. Vaidya, “Sum of squares based convex approach for optimal control synthesis,” in _2021 29th Mediterranean Conference on Control and Automation (MED)_. IEEE, 2021, pp. 1270–1275.
  * [17] R. Rajaram, U. Vaidya, M. Fardad, and B. Ganapathysubramanian, “Stability in the almost everywhere sense: A linear transfer operator approach,” _Journal of Mathematical analysis and applications_ , vol. 368, no. 1, pp. 144–156, 2010.
  * [18] Z. Artstein, “Stabilization with relaxed controls,” _Nonlinear Analysis: Theory, Methods & Applications_, vol. 7, no. 11, pp. 1163–1173, 1983.
  * [19] E. D. Sontag, “A lyapunov-like characterization of asymptotic controllability,” _SIAM journal on control and optimization_ , vol. 21, no. 3, pp. 462–471, 1983.
  * [20] P. Wieland and F. Allgöwer, “Constructive safety using control barrier functions,” _IFAC Proceedings Volumes_ , vol. 40, no. 12, pp. 462–467, 2007.