# Tailoring Generative Adversarial Networks for Smooth Airfoil Design

Joyjit Chattoraj  Institute of High Performance Computing
Agency for Science, Technology and Research
Singapore 138632, Republic of Singapore
<EMAIL_ADDRESS>Jian Cheng Wong  Institute of High Performance
Computing
Agency for Science, Technology and Research
Singapore 138632, Republic of Singapore
<EMAIL_ADDRESS>Zhang Zexuan  Manna Dai  School of Continuing and
Lifelong Education
National University of Singapore
Singapore 119260, Republic of Singapore
<EMAIL_ADDRESS>Institute of High Performance Computing
Agency for Science, Technology and Research
Singapore 138632, Republic of Singapore
<EMAIL_ADDRESS>Xia Yingzhi  Institute of High Performance
Computing
Agency for Science, Technology and Research
Singapore 138632, Republic of Singapore
<EMAIL_ADDRESS>Li Jichao  Xu Xinxing  Institute of High
Performance Computing
Agency for Science, Technology and Research
Singapore 138632, Republic of Singapore
<EMAIL_ADDRESS>Institute of High Performance Computing
Agency for Science, Technology and Research
Singapore 138632, Republic of Singapore
<EMAIL_ADDRESS>Ooi Chin Chun  Institute of High Performance
Computing
Agency for Science, Technology and Research
Singapore 138632, Republic of Singapore
<EMAIL_ADDRESS>Yang Feng  Dao My Ha  Institute of High Performance
Computing
Agency for Science, Technology and Research
Singapore 138632, Republic of Singapore
<EMAIL_ADDRESS>Institute of High Performance Computing
Agency for Science, Technology and Research
Singapore 138632, Republic of Singapore
<EMAIL_ADDRESS>Liu Yong  Institute of High Performance Computing
Agency for Science, Technology and Research
Singapore 138632, Republic of Singapore
<EMAIL_ADDRESS>
###### Abstract

In the realm of aerospace design, achieving smooth curves is paramount,
particularly when crafting objects such as airfoils. Generative Adversarial
Network (GAN), a widely employed generative AI technique, has proven
instrumental in synthesizing airfoil designs. However, a common limitation of
GAN is the inherent lack of smoothness in the generated airfoil surfaces. To
address this issue, we present a GAN model featuring a customized loss
function built to produce seamlessly contoured airfoil designs. Additionally,
our model demonstrates a substantial increase in design diversity compared to
a conventional GAN augmented with a post-processing smoothing filter.

###### Index Terms:

GAN, generative AI, custom loss, airfoil design

##  I Introduction

Generative AI models, such as Generative Adversarial Network (GAN) and
Variational Autoencoder (VAE), have demonstrated promising outcomes in
engineering design. Applications span a range of fields, encompassing
structural optimization, materials design, and shape synthesis [1]. In recent
years, particularly in the domain of airfoil design [2], numerous variations
of GAN models have been investigated to discover novel shapes capable of
achieving optimal aerodynamic performance across a diverse set of operating
conditions [3, 4, 5, 6, 7].

However, it is frequently reported in the literature that GAN has limitations
generating smooth airfoil curves [8, 9, 10, 11]. A commonly employed strategy
to address this issue involves integrating a smoothing filter as a post-
processing technique [8].

Chen et al [12, 13, 14] proposed a novel model, BézierGAN, to avoid post-
processing of airfoil curves, instead, they employed a smoothing method, in
this case, Bézier curve, as an integral part of the GAN generator. During GAN
training, the model learns the optimized values for Bézier parameters. In
addition, a set of regularization terms has been imposed that enhances the
complexity of the model, which could potentially limit the model’s
generalizability across various applications.

Yonekura et al [15] used Conditional Wasserstein GAN (CWGAN-GP) and found that
the model can generate smooth airfoil curves. Notably, WGAN was originally
proposed to tackle the vanishing gradient problem in GAN models [16]. Hence,
whether the primary reason for non-smoothness lies in a vanishing gradient
remains an open question. Tan et al [9] with a similar CWGAN-GP model could
not generate smooth airfoil curves consistently. One noticeable difference
between these two models is that Tan et al [9] considered more than one
parameter for airfoil design including a shape-determined parameter (area).

Wang et al [17] achieved smooth airfoil shapes by integrating a VAE and a GAN,
referred to as VAEGAN, originally proposed in [18]. Here, VAE acts as a GAN
generator, and it can generate smooth curves without imposing any explicit
methodology for curve smoothness.

Figure 1: (a) The generator (G) of a trained conditional-GAN takes a noise
vector and a label vector as inputs and the generator outputs a 2D airfoil
curve. The generator is trained with a loss function $L^{(\text{G})}$ which
includes a binary cross entropy loss ($L_{\text{CE}}$) and a mean square error
loss $L_{\text{S}}$ to impose curve smoothing. (b) The generator outputs non-
smooth curves if $L^{(\text{G})}$ only includes $L_{\text{CE}}$, i.e., without
$L_{\text{S}}$. (c) The generator outputs desired smooth curves when
$L_{\text{S}}$ is included in the generator loss in addition to
$L_{\text{CE}}$.

In our proposed methodology, we attempt to tackle the non-smoothness problem
of GAN explicitly. We propose a custom loss function for the generator model
of conditional GAN (Fig. 1a). This loss function calculates a mean square
error if a generated curve deviates from its moving average position. The
custom loss function improves the ability of the GAN model to transition from
generating non-smooth curves to creating smooth airfoil curves, as illustrated
in Fig. 1b,c. Moreover, we find that the custom loss facilitates the GAN model
to generate more diverse airfoil shapes than the original GAN without the
custom loss function and incorporated with a smoothing filter as a post-
processing technique. As the custom loss requires computing only the moving
average, this approach is generalizable and can be easily adopted in other
engineering designs to impose smoothness.

The remaining sections of this article are organized as follows: (i) Dataset:
Here, we provide an overview of the airfoil data utilized for training the GAN
models. (ii) Methodology: This section delves into the specifics of the GAN
architecture, custom loss, and the training procedure. (iii) Results: We
analyze the effectiveness of the custom loss function in this segment. (iv)
Conclusions: Finally, we summarize the main findings and discuss potential
avenues for future research.

##  II Dataset

Our dataset consists of 1399 airfoils which are obtained from the UIUC Airfoil
Data Site [19] using an open-source pre-processing code [20, 21]. These
airfoils are essentially 2D curves with a set of X and Y coordinates, where X
coordinates strictly range from 0 to 1. The Y coordinates are flexible and
typically vary between $-$0.3 to 0.4 (Fig. 2). The maximum thickness ($\tau$),
i.e., measured perpendicular to the camber line of an airfoil, varies between
0.04 and 0.38 with a median value of 0.12.

Subsequently, these airfoils underwent assessment under three distinct
operating conditions: Reynolds number ($Re$), Mach number ($M$), and angle of
attack ($\alpha$). The airfoil’s aerodynamic performance was evaluated by
calculating the lift coefficient ($c_{l}$) and drag coefficient ($c_{d}$)
utilizing XFOIL [22]. The numerical simulation generated a comprehensive
dataset comprising 654,963 samples, averaging 468 samples per airfoil,
encompassing different operating conditions, i.e., $Re\ (2000000-12000000)$,
$M\ (0-0.3)$, and $\alpha\ (0-33)$, with their calculated $c_{l}$ and $c_{d}$.

Figure 2: An original airfoil (solid line) after reducing the number of points to 38 (symbols). TABLE I: 8 Classes and respective number of airfoils for training  Class | $c_{l}/c_{d}$ | $\alpha$ | $\tau$ | Airfoils
---|---|---|---|---
000 | low | low | low | 740
001 | low | low | high | 627
010 | low | high | low | 742
011 | low | high | high | 649
100 | high | low | low | 739
101 | high | low | high | 633
110 | high | high | low | 717
111 | high | high | high | 566

In essence, our GAN design challenge is framed as a classification problem,
where the objective is to create airfoil shapes that align with a particular
class. We define 8 classes to create airfoil shapes from 3 parameters: $\tau$
representative of airfoil geometry, $\alpha$ representative of operating
conditions, and the ratio $c_{l}/c_{d}$ representative of aerodynamic
performance (Table I). A class is characterized by combinations of low and
high values of $\tau$, $\alpha$, and $c_{l}/c_{d}$. The thresholds for
distinguishing between low and high values are established based on the median
values of these parameters, which are 0.12, 10, and 100, respectively. We
proceed to create a training dataset by determining the class to which each
airfoil belongs. Note that while an airfoil may be associated with multiple
classes, its presence in a class is considered only once. We find that the
number of airfoils varies between 566 and 742 in 8 classes, a total of 5413
samples, implying a reasonably balanced dataset for GAN training.

##  III Methodology

Figure 3: The architecture of conditional-GAN consists of a generator (G) and
a discriminator (D). Both G and D are fully connected deep neural networks.

###  III-A Pre-processing

We employ a pre-processing methodology to standardize airfoil coordinates for
GAN training. This process ensures uniform X coordinates across all airfoils,
maintaining consistency in their sequence, while strategically reducing the
number of points. Emphasis is placed on critical airfoil segments, such as the
leading and trailing edges. Employing a cosine distribution, we set 20 X
coordinates; excluding the two endpoints (X=0 and X=1), the remaining 18 X
coordinates each have 2 corresponding Y coordinates for the upper and lower
surfaces (Fig. 2). Consequently, a total of 38 Y coordinates distinguish one
airfoil from another, and these 38 Y coordinates serve as inputs for GAN
training.

###  III-B GAN architecture

We adopt the conditional GAN technique [23] which comprises a generator (G)
and a discriminator (D), depicted in Fig. 3. G is a fully connected deep
neural network with the rectified linear unit (ReLU) as an activated function.
G receives an input vector formed by concatenating a random noise vector (128
dimensions, Gaussian noise with mean 0 and standard deviation 1) and a label
vector (6 dimensions). Label vector represents the 8 classes (Table I), with
an example being 100101 label vector for class-011. G creates 38 Y coordinates
belonging to the 2D curve of an airfoil.

D receives an input vector formed by concatenating the real coordinates (38
dimensions) and generated coordinates (38 dimensions) from G associated with
the label vector (6 dimensions). D determines the probability of the generated
coordinates to be real. D is also a fully connected deep neural network; the
last layer of D uses the sigmoid activation function, and all the other
neurons use ReLU.

###  III-C Custom loss function

To enforce smoothness on the 2D curve of an airfoil, we tailor the generator
training by customizing the loss function $L^{\text{(G)}}$ as

$L^{\text{(G)}}=L_{\text{CE}}+L_{\text{S}},$ (1)

where $L_{\text{CE}}$ is the standard binary cross-entropy loss, whereas
$L_{\text{S}}$ introduces a penalty for any deviation of the generated Y
coordinates (with dimensions $n_{c}=38$) from their respective moving average
points. We compute the moving average $\bar{Y}_{i}$ for all Y coordinates with
a window of 3 points as

$\bar{Y}_{i}=\frac{1}{3}\sum_{j=-1}^{1}Y_{i+j},$ (2)

where $i=1,2,\ldots,n_{c}$, and the index $i+j$ iterates through 1 and $n_{c}$
in a cyclic manner. The coordinates are arranged in a loop where $i=1$
corresponds to the endpoint at $X_{i}=1$, $i=2$ corresponds to the upper
surface point at $X_{i}=0.996$, and $i=n_{c}$ corresponds to the lower surface
point at $X_{i}=0.996$. The deviation $\Delta_{i}$ between a Y coordinate
$Y_{i}$ and its moving average point is then

$\Delta_{i}=Y_{i}-\bar{Y}_{i}.$ (3)

We define the smoothing loss $L_{\text{S}}$ as the mean square deviation of
$\Delta_{i}$ as shown below

$L_{\text{S}}=\frac{\omega}{n_{c}}\sum_{i=1}^{n_{c}}\Delta_{i}^{2},$ (4)

where the coefficient $\omega$ determines the strength of $L_{\text{S}}$ with
respect to $L_{\text{CE}}$. For the discriminator, we use the binary cross-
entropy loss function only.

###  III-D GAN training

The generator (G) and discriminator (D) are trained simultaneously in a
competitive manner incorporating the custom loss function (1). In a training
loop, G creates 2D airfoil curves, and D evaluates both real and generated
airfoils without knowing the source. For both G and D, we use ADAM optimizer
and the learning rate $10^{-4}$. We set epochs at $30,000$. During training we
observed that beyond the first several thousand epochs, both the G-loss and
D-loss converged to approximately the same value, maintaining a steady-state
until the end of the training process. We experimented with three distinct
values for $\omega$ in (4): 1, 10, and 100. Our findings indicated that
setting $\omega$ to 10 resulted in optimal smoothing for the airfoils.

##  IV Results

Figure 4: Generated airfoil samples from the GAN trained without the smoothing
loss and augmented with a post-processing smoothing filter. Each panel
represents a class defined by low/high (0/1) $c_{l}/c_{d}$, $\alpha$, and
$\tau$. It displays 600 airfoil samples and their mean shape (dashed curve).

To evaluate the efficacy of our proposed custom loss function (1), we
conducted training experiments with two GAN models: one without a smoothing
loss and another incorporating the smoothing loss function, to be referred to
as smoothGAN, henceforth. We then compared their performances in generating
airfoils. A total of 600 airfoil samples were generated for each of the 8
classes, and we assessed these samples using three statistical metrics:

  1. 1.

$\text{ACC}^{(\tau)}$: Accuracy that measures the percentage of the total
airfoils possessing the correct (low/high) thickness at a given class.

  2. 2.

$\sigma^{(\tau)}$: Diversity in thickness, calculated as the standard
deviation of thickness after being rescaled by their maximum observed value.

  3. 3.

$S$: Diversity in shape, following [24], it is calculated as the average
distance of the mean shape $Y^{\text{(m)}}$ from all the generated samples
$\left\langle\sqrt{\sum_{i=1}^{n_{c}}(Y_{i}-Y_{i}^{\text{(m)}})^{2}}\right\rangle_{\text{samples}}$.

Figure 5: Generated airfoil samples from the GAN model trained with the
smoothing loss. Each panel represents a class defined by low/high (0/1)
$c_{l}/c_{d}$, $\alpha$, and $\tau$. It displays 600 airfoil samples and their
mean shape (dashed curve).

###  IV-A GAN without smoothing loss function

We found that GAN which is trained without the smoothing loss function always
generated non-smooth 2D curves, an example can be seen in Fig. 1b. A
Savitzky‐Golay filter [25] is further augmented with the GAN model as a post-
processing technique that smooths out the GAN-generated samples, which are
displayed in Fig. 4. We illustrate the mean shape of airfoils for each class
in Fig. 6a to showcase the shape diversity among classes. Additionally, the
distribution of $\tau$ is depicted in Fig. 6c,e as an alternative measure of
geometric diversity.

###  IV-B smoothGAN: GAN with smoothing loss function

The smoothGAN model trained with the custom loss function generates smooth 2D
curves, which do not require any post-processing. 600 samples for each of 8
classes, a total of 4800, plotted in Fig. 5, all exhibits smooth shapes. It is
evident in Fig. 5 and also in Fig. 6b,d,f that in comparison to GAN, smoothGAN
achieves more diversity in generated shapes. We further quantitatively
demonstrate the diversity through two metrics: $\sigma^{(\tau)}$ and $S$.

Figure 6: Comparison of airfoils which are generated from GAN trained without smoothing loss and augmented with a post-processing smoothing filter (3 panels in left column), and GAN trained with smoothing loss (3 panels in right column). (a) and (b) show the mean shape of 600 airfoils, a total of 8 shapes corresponding to 8 classes. Distribution of $\tau$ of 600 airfoils for a fixed class: (c) and (d) for class-000, and (e) and (f) for class-111. TABLE II: Performance of GAN Class | w/o Smoothing Loss | with Smoothing Loss
---|---|---
| $\text{ACC}^{(\tau)}(\%)$ | $\sigma^{(\tau)}$ | $S$ | $\text{ACC}^{(\tau)}(\%)$ | $\sigma^{(\tau)}$ | $S$
000 | 100 | 0.007 | 0.003 | 100 | 0.057 | 0.017
001 | 100 | 0.024 | 0.009 | 100 | 0.044 | 0.017
010 | 100 | 0.007 | 0.003 | 100 | 0.061 | 0.018
011 | 100 | 0.028 | 0.01 | 100 | 0.036 | 0.02
100 | 100 | 0.005 | 0.003 | 100 | 0.05 | 0.017
101 | 100 | 0.027 | 0.009 | 100 | 0.044 | 0.02
110 | 100 | 0.008 | 0.003 | 100 | 0.057 | 0.018
111 | 100 | 0.03 | 0.009 | 100 | 0.043 | 0.02

The comparative performance of these two GAN models is summarized in Table II.
Both models achieve a 100$\%$ accuracy ($\text{ACC}^{(\tau)}$), indicating
that the generated airfoils consistently exhibit the correct (high/low)
thickness corresponding to their respective classes. However, smoothGAN
exhibits notably elevated values for both $\sigma^{(\tau)}$ and $S$. As seen
in Fig. 6, distribution of $\tau$ is broader for smoothGAN than the former GAN
leading to approximately 1.5 to 10 times higher $\sigma^{(\tau)}$ values.
Similarly, smoothGAN exhibits approximately 2 to 6 times higher $S$ values
than the former GAN. It is noteworthy that the minimum value of S is zero
implying that the model generates only one airfoil shape. This situation may
occur in the presence of mode collapse, a common drawback associated with GANs
[26].

##  V Conclusions

In this study, we developed a methodology to generate smooth airfoils using
GAN. We showed that conditional GAN has limitations in generating smooth
curves, supporting the previous findings reported in the literature.
Typically, this shortcoming has been addressed by augmenting a smoothing
filter with GAN as a post-processing technique. We demonstrated explicitly
that employing such a post-processing technique can attain high accuracy.
Nevertheless, this approach may generate a restricted range of airfoil shapes,
and the absence of diversity in shape can constrain the GAN’s ability to
explore novel designs.

Our proposed methodology integrates a custom loss function with GAN,
comprising a smoothing loss that penalizes deviations of the generated shape
from the reference shape assessed through a moving average. This approach
enhances the GAN’s capability to consistently generate smooth airfoils, with 2
to 10 times greater diversity compared to the original GAN.

In future applications, we plan to implement our proposed GAN methodology in
various scenarios to generate smooth curves and surfaces. One potential
application could be in the realm of 3D airfoil design.

## Acknowledgment

This research is supported by Agency for Science Technology and Research
(A*STAR) under the AME Programmatic project: Explainable Physics-based AI for
Engineering Modelling & Design (ePAI) [Award No. A20H5b0142].

## References

  * [1] L. Regenwetter, A. H. Nobari, and F. Ahmed, “Deep generative models in engineering design: A review,” _Journal of Mechanical Design_ , vol. 144, no. 7, p. 071704, 2022.
  * [2] J. Li, X. Du, and J. R. Martins, “Machine learning in aerodynamic shape optimization,” _Progress in Aerospace Sciences_ , vol. 134, p. 100849, Oct. 2022.
  * [3] J. Li, M. Zhang, J. R. Martins, and C. Shu, “Efficient aerodynamic shape optimization with deep-learning-based geometric filtering,” _AIAA journal_ , vol. 58, no. 10, pp. 4243–4259, 2020.
  * [4] E. Yilmaz and B. German, “Conditional generative adversarial network framework for airfoil inverse design,” in _AIAA aviation 2020 forum_ , 2020, p. 3185.
  * [5] A. Heyrani Nobari, W. Chen, and F. Ahmed, “Pcdgan: A continuous conditional diverse generative adversarial network for inverse design,” in _Proceedings of the 27th ACM SIGKDD conference on knowledge discovery & data mining_, 2021, pp. 606–616.
  * [6] J. Li and M. Zhang, “On deep-learning-based geometric filtering in aerodynamic shape optimization,” _Aerospace Science and Technology_ , vol. 112, p. 106603, May 2021.
  * [7] Q. Du, T. Liu, L. Yang, L. Li, D. Zhang, and Y. Xie, “Airfoil design and surrogate modeling for performance prediction based on deep learning method,” _Physics of Fluids_ , vol. 34, no. 1, 2022.
  * [8] G. Achour, W. J. Sung, O. J. Pinon-Fischer, and D. N. Mavris, “Development of a conditional generative adversarial network for airfoil shape optimization,” in _AIAA Scitech 2020 Forum_ , 2020, p. 2261.
  * [9] X. Tan, D. Manna, J. Chattoraj, L. Yong, X. Xinxing, D. M. Ha, and Y. Feng, “Airfoil inverse design using conditional generative adversarial networks,” in _2022 17th International Conference on Control, Automation, Robotics and Vision (ICARCV)_. IEEE, 2022, pp. 143–148.
  * [10] G. B. Santos, A. V. Pantaleão, and L. O. Salviano, “Using deep generative adversarial network to explore novel airfoil designs for vertical-axis wind turbines,” _Energy Conversion and Management_ , vol. 282, p. 116849, 2023.
  * [11] K. Wada, K. Suzuki, and K. Yonekura, “Physics-guided training of gan to improve accuracy in airfoil design synthesis,” _arXiv preprint arXiv:2308.10038_ , 2023.
  * [12] W. Chen and M. Fuge, “B$\backslash$’eziergan: Automatic generation of smooth curves from interpretable low-dimensional parameters,” _arXiv preprint arXiv:1808.08871_ , 2018.
  * [13] W. Chen, K. Chiu, and M. Fuge, “Aerodynamic design optimization and shape exploration using generative adversarial networks,” in _AIAA Scitech 2019 Forum_ , 2019, p. 2351.
  * [14] W. Chen, K. Chiu, and M. D. Fuge, “Airfoil design parameterization and optimization using bézier generative adversarial networks,” _AIAA journal_ , vol. 58, no. 11, pp. 4723–4735, 2020.
  * [15] K. Yonekura, N. Miyamoto, and K. Suzuki, “Inverse airfoil design method for generating varieties of smooth airfoils using conditional wgan-gp,” _Structural and Multidisciplinary Optimization_ , vol. 65, no. 6, p. 173, 2022.
  * [16] M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein generative adversarial networks,” in _International conference on machine learning_. PMLR, 2017, pp. 214–223.
  * [17] Y. Wang, K. Shimada, and A. B. Farimani, “Airfoil gan: Encoding and synthesizing airfoils for aerodynamic shape optimization,” _Journal of Computational Design and Engineering_ , p. qwad046, 2023.
  * [18] A. B. L. Larsen, S. K. Sønderby, H. Larochelle, and O. Winther, “Autoencoding beyond pixels using a learned similarity metric,” in _International conference on machine learning_. PMLR, 2016, pp. 1558–1566.
  * [19] M. S. Selig. (1996) Uiuc airfoil data site. [Online]. Available: https://m-selig.ae.illinois.edu/ads/coord_database.html
  * [20] J. Li, S. He, and J. R. Martins, “Data-driven constraint approach to ensure low-speed performance in transonic aerodynamic shape optimization,” _Aerospace Science and Technology_ , vol. 92, p. 536–550, Sep. 2019.
  * [21] J. Li, M. A. Bouhlel, and J. R. Martins, “Data-based approach for fast airfoil analysis and optimization,” _AIAA Journal_ , vol. 57, no. 2, pp. 581–596, 2019.
  * [22] M. Drela, “Xfoil: An analysis and design system for low reynolds number airfoils,” in _Low Reynolds Number Aerodynamics: Proceedings of the Conference Notre Dame, Indiana, USA, 5–7 June 1989_. Springer, 1989, pp. 1–12.
  * [23] M. Mirza and S. Osindero, “Conditional generative adversarial nets,” _arXiv preprint arXiv:1411.1784_ , 2014.
  * [24] N. C. Brown and C. T. Mueller, “Quantifying diversity in parametric design: a comparison of possible metrics,” _AI EDAM_ , vol. 33, no. 1, pp. 40–53, 2019.
  * [25] W. H. Press and S. A. Teukolsky, “Savitzky-golay smoothing filters,” _Computers in Physics_ , vol. 4, no. 6, pp. 669–672, 1990.
  * [26] D. Yang, S. Hong, Y. Jang, T. Zhao, and H. Lee, “Diversity-sensitive conditional generative adversarial networks,” _arXiv preprint arXiv:1901.09024_ , 2019.


# On the continuity and smoothness of the value function in reinforcement
learning and optimal control

Hans Harder1 and Sebastian Peitz1 1Department of Computer Science, Paderborn
University, Paderborn, Germany<EMAIL_ADDRESS>
###### Abstract

The _value function_ plays a crucial role as a measure for the cumulative
future reward an agent receives in both reinforcement learning and optimal
control. It is therefore of interest to study how similar the values of
neighboring states are, i.e., to investigate the continuity of the value
function. We do so by providing and verifying upper bounds on the value
function’s modulus of continuity. Additionally, we show that the value
function is always Hölder continuous under relatively weak assumptions on the
underlying system and that non-differentiable value functions can be made
differentiable by slightly “disturbing” the system.

##  I Introduction

In reinforcement learning, an agent is put into some environment, acts in it,
accumulates traces of rewards, and tries to find actions that will give it the
best possible outcome. This is achieved by learning a policy (or feedback law)
that maximizes the cumulative future reward it receives, measured by the value
function.

A well-known two-phase approach to this issue is _policy iteration_ (see
[14]): In the _policy evaluation_ phase, the agent approximates the value
(i.e., cumulative reward) for each state by collecting trajectories following
the current policy. In the _policy improvement_ phase, the agent adapts its
policy by acting greedily with respect to the value function.

We study properties of the value function for a fixed policy, as in the policy
evaluation phase. (In optimal control, it is convention to only call the cost
functional evaluated for the optimal policy the “value function”. In our case,
every policy corresponds to a value function.) Essentially, this allows us to
discard the policy in our analysis by viewing it as part of the environment.
This setup is common in the analysis of temporal-difference learning, see,
e.g., [16, 15, 5].

For illustration purposes, let us introduce the formulation for the case where
both environment and policy are deterministic. We are given a set of states
$S$ and actions $A$, a model of the environment $\Phi:S\times A\rightarrow S$,
a reward function $r:S\times A\rightarrow\mathbb{R}$ and a discount factor
$\gamma\in(0,1)$. When a policy $\pi:S\rightarrow A$ is fixed, the model
reduces to a dynamical system $\Phi_{\pi}:S\rightarrow S$ defined by
$\Phi_{\pi}(x)=\Phi(x,\pi(x))$, and the reward function is simplified by
introducing $r_{\pi}(x)=r(x,\pi(x))$. The _value function_ is

$\displaystyle
v_{\pi}(x)=\sum_{n=0}^{\infty}\gamma^{n}r_{\pi}(\Phi_{\pi}^{n}(x)),$ (1)

where $(\cdot)^{n}$ denotes the $n$-fold application of $\Phi_{\pi}$.
Accordingly, one can define the _state-action value_ as

$\displaystyle q_{\pi}(x,a)=r(x,a)+\gamma v_{\pi}(\Phi(x,a)).$ (2)

Equation 1 is more or less the form that we are interested in, and we will
discard the subscript $(\cdot)_{\pi}$ in the forthcoming investigations by
assuming that $\pi$ is fixed. In other words, we will think of $r$ and $\Phi$
as functions from $S$ to $\mathbb{R}$ and $S$ respectively, independent of any
action.

Works such as [12] or [10] approximate (2) using neural networks and rely on
taking derivatives with respect to $a$. In practice, this is not a problem
since these models are differentiable by definition. However, when it comes to
the _true_ state-action value function, differentiability of $v_{\pi}$ is
required. But this is an assumption that needs good justification; there are
value functions that are _nowhere_ differentiable:

######  Proposition 1 (see [17]).

Let $\Phi(x)=4x(1-x)$ be the logistic map on $S=[0,1]$, put $r(x)=x$ and let
$\gamma\in[\frac{1}{2},1)$. Then
$v(x)=\sum_{n=0}^{\infty}\gamma^{n}r(\Phi^{n}(x))$ is nowhere differentiable.

######  Remark 1.

One can obtain the system of Proposition 1 by defining $r(x,a)=x$ and
$\Phi(x,a)=a\cdot 4x(1-x)$ for $A=S=[0,1]$. The policy $\pi(x)=1$ induces the
logistic map.

On the other hand, we _can_ get differentiability by “disturbing” the system.
Informally, if we define a new random dynamical system via “original system +
noise”, then the resulting value function is differentiable. We have
illustrated this in Figure 1 using the logistic map from Proposition 1.

Figure 1: The value function $v$ from Proposition 1 for the discount factor
$\gamma=0.8$. The “smoothed” version $w$ is the value function that one
obtains when disturbing the same system using Gaussian noise with standard
deviation $\sigma=0.01$, cf. Proposition 4.

In this paper, we will study the general case, i.e., we also analyze nowhere
differentiable value functions. It turns out that these functions are Hölder
continuous under weak assumptions, which is useful since this can be used to
bound the variance of $v$ if evaluated at random states:

######  Proposition 2.

Suppose $v:S\rightarrow\mathbb{R}$ for $S\subseteq\mathbb{R}^{N}$ is Hölder
continuous, i.e., there are $A\geq 0$ and $0<\beta\leq 1$ such that
$|v(x)-v(x^{\prime})|\leq A||x-x^{\prime}||_{2}^{\beta}$ for all
$x,x^{\prime}\in S$. Let $X$ be an $S$-valued random variable. Then

$\mathrm{Var}[v(X)]\leq\frac{A^{2}}{2^{1-\beta}}(\sum_{i=1}^{N}\mathrm{Var}[X_{i}])^{\beta}.$
(3)

###### Proof.

Any real-valued random variable $Y$ satisfies
$\mathrm{Var}[Y]=\frac{1}{2}\mathbb{E}[(Y-Y^{\prime})^{2}]$, where
$Y^{\prime}$ is an independent copy of $Y$. Hence, by Jensen’s inequality:

$\displaystyle\mathrm{Var}[v(X)]$
$\displaystyle=\textstyle\frac{1}{2}\mathbb{E}[(v(X)-v(X^{\prime}))^{2}]$
$\displaystyle\leq\textstyle\frac{1}{2}A^{2}\mathbb{E}[||X-X^{\prime}||^{2\beta}_{2}]$
$\displaystyle\leq\textstyle\frac{1}{2}A^{2}\mathbb{E}[||X-X^{\prime}||^{2}_{2}]^{\beta}$
$\displaystyle=\textstyle\frac{1}{2^{1-\beta}}A^{2}(\sum_{i=1}^{N}\frac{1}{2}\mathbb{E}[(X_{i}-X_{i}^{\prime})^{2}])^{\beta},$

which is just (3). ∎

We study the continuity of $v$ by categorizing it as “Lipschitz continuous”,
“Hölder continuous” or as a mix thereof, depending on the relationship between
the Lipschitz constant of the system and the discount factor. Hölder
continuity of $v$ is known under weak assumptions for the deterministic and
continuous-time setting, see [2], and it is known for the discrete-time
setting under similar assumptions as we do, see [4]. We improve these results
by giving a closed-form solution to an upper bound of $v$’s modulus of
continuity and show that it cannot be improved for linear reward functions.
Additionally, we consider both random and deterministic systems, as well as
continuous-time and discrete-time systems. Finally, we show how “disturbed”
systems can yield differentiable value functions. Besides [2] and [4], there
have been investigations under varying assumptions by [13, 11, 9, 3] and [8].
Proposition 1 is a result of [17], who studied the equation in a seemingly
unrelated setting.

#### Organization

In Section III, we derive continuity bounds for the value function in both
continuous- and discrete-time settings. We show how these bounds improve on
and entail known results about the Lipschitz resp. Hölder continuity of the
value function. The derived upper bound from Section III is verified
experimentally and shown to be sharp in Section IV. Finally, we show how to
obtain differentiability as soon as we add a little bit of Gaussian noise in
Section V.

##  II Definitions

We fix a probability triple $(\Omega,\mathcal{F},P)$ and a metric space
$(S,d)$. We call $S$ the _state space_ and suppose that $d$ is bounded;
$\sup_{x,x^{\prime}\in S}d(x,x^{\prime})\leq D$. For technical reasons, we
also have to assume that $S$ is part of a measurable space $(S,\Sigma)$ so
that measurable functions from $S$ to $\mathbb{R}$ can be composed or
integrated. Forthcoming, we will frequently use the _modulus (of continuity)_
of a uniformly continuous function $f:S\rightarrow\mathbb{R}$. This is the
function defined on $\mathbb{R}_{\geq 0}$ by

$\displaystyle W_{f}(d_{0})=\sup_{x,x^{\prime}\in
S}\\{|f(x)-f(x^{\prime})|:d(x,x^{\prime})\leq d_{0}\\}.$

Intuitively, $W_{f}$ is the best function that one can use when bounding
$|f(x)-f(x^{\prime})|$ in terms of $d(x,x^{\prime})$.

Next, we define “random systems”. Such systems model random transitions
between states over given periods. We also introduce the notion of _Lipschitz
continuity in Expectation_ , or LE-continuity, which is used to bound the
divergence between nearby states:

######  Definition 1.

Let $T\in\\{\mathbb{R}_{\geq 0},\mathbb{N}_{0}\\}$ be the time domain. A
measurable function $\Phi:T\times S\times\Omega\rightarrow S$, written
$\Phi_{t}(x)=\Phi_{t}(x,\cdot)=\Phi(t,x,\cdot)$, is called a random system (so
that $\Phi_{t}(x)$ is a random variable). If $T=\mathbb{R}_{\geq 0}$
($T=\mathbb{N}_{0}$), we call $\Phi$ a continuous-time (discrete-time) system.
We say that $\Phi$ is LE-continuous with constant $L\geq 0$ when

$\displaystyle\mathbb{E}[d(\Phi_{t}(x),\Phi_{t}(x^{\prime}))]\leq
L^{t}d(x,x^{\prime})\quad\forall t\in T,\,x,x^{\prime}\in X.$

This definition is very general and not strong enough to obtain the well-known
recursive representation of the value function. For this one needs stochastic
flows [6, 7], random dynamical systems [1] or Markov chains, but the latter
does not allow for a simple characterization of LE-continuity (see [7]).
However, the required properties from Definition 1 already suffice to get all
desired results since we study value functions in their infinite-sum (or
infinite-integral) form, which does not require the Markov property.

######  Example 1.

Let $\xi_{0},\xi_{1},\dots:\Omega\rightarrow\mathbb{R}^{n}$ be random vectors
and suppose $L\geq 0$ and $\varphi:S\times\mathbb{R}^{n}\rightarrow S$ satisfy

$\displaystyle d(\varphi(x,y),\varphi(x^{\prime},y))\leq
Ld(x,x^{\prime}),\quad\forall x,x^{\prime}\in S,y\in\mathbb{R}^{n}.$

Then the system defined by $\Phi_{n+1}(x)=\varphi(\Phi_{n}(x),\xi_{n})$ and
$\Phi_{0}(x)=x$ is LE-continuous with constant $L$.

######  Remark 2.

If the discrete metric $d(x,x^{\prime})=1_{x\neq x^{\prime}}$ is used, e.g.,
if $S$ is countable, then LE-continuity is trivially satisfied for $L=1$. LE-
continuity is therefore only necessary when $S$ is uncountable, e.g. because
it is an open or closed subset of $\mathbb{R}^{N}$. Intuitively, it requires
that trajectories cannot diverge too fast if the distance of the initial
states is small.

######  Remark 3.

The at worst exponential divergence of trajectories is a common assumption,
see, e.g., [2]. We alter it slightly by introducing the expectation. This does
not make a difference for deterministic systems.

If a random system is equipped with a reward function and a discount factor,
one obtains a corresponding value function. We discard the policy and the set
of actions:

######  Setting 1.

Let $r:S\rightarrow\mathbb{R}$ be a reward function for which there is a
concave and monotonically increasing $R:\mathbb{R}\rightarrow\mathbb{R}$ with
$R(0)=0$ such that $|r(x)-r(x^{\prime})|\leq R(d(x,x^{\prime}))$ for all
$x,x^{\prime}\in S$. Moreover, let $\Phi$ be a random system that is LE-
continuous with constant $L>0$. Let $\gamma\in(0,1)$ be a discount factor.
Define the value function as

$\displaystyle
v(x)=\begin{cases}\int_{0}^{\infty}\gamma^{t}\mathbb{E}[r(\Phi_{t}(x))]dt&T=\mathbb{R}_{\geq
0},\\\
\sum_{n=0}^{\infty}\gamma^{n}\mathbb{E}[r(\Phi_{n}(x))]&T=\mathbb{N}_{0}.\end{cases}$

##  III Continuity of the value function

Assuming that 1 holds, the overall goal of this section is to bound the
difference $|v(x)-v(x^{\prime})|$ by a quantity $H(d(x,x^{\prime}))$, where
$H$ is a function that is to be derived. The derived $H$ depends on multiple
factors: the Lipschitz constant $L$ and the bound $R$, the discount factor
$\gamma$, the maximal distance between two states $D$, and finally the chosen
time domain. In Theorem 1, we derive $H$ for the continuous case since the
discrete case can be reduced to it. We assume $L>1$ since $v$ is trivially
Lipschitz continuous otherwise (see Section III-A for the case $L\leq
1<1/\gamma$).

Figure 2: Visual depiction of the idea in Theorem 1.

######  Theorem 1.

Suppose 1 holds with $T=\mathbb{R}_{\geq 0}$ and $L>1$. Define for $d_{0}>0$:

$\displaystyle K(d_{0})=\begin{dcases}\frac{\ln L}{\ln(\gamma
L)}D^{1-\beta}d_{0}^{\beta}+\frac{\ln\gamma}{\ln(\gamma L)}d_{0}&L\neq
1/\gamma,\\\ (\ln(D/d_{0})+1)d_{0}&L=1/\gamma,\end{dcases}$

with $K(0):=0$ and $\beta=\ln(1/\gamma)/\ln L$. Setting

$\displaystyle H(d_{0})=\frac{R(K(d_{0}))}{\ln(1/\gamma)},$ (4)

we have

$\displaystyle|v(x)-v(x^{\prime})|\leq H(d(x,x^{\prime})),\quad\forall
x,x^{\prime}\in S.$ (5)

###### Proof.

Take $x,x^{\prime}\in S$ and define $d_{0}:=d(x,x^{\prime})$. For $d_{0}=0$,
use that $d$ is a metric to obtain $x=x^{\prime}$ and thus
$|v(x)-v(x^{\prime})|=0$. We have $K(0)=0$ and $R(0)=0$ by definition, hence
(5) applies. Forthcoming, suppose $d_{0}>0$. Using the concavity of $R$ with
Jensen’s inequality, we get

$\displaystyle|v(x)-v(x^{\prime})|\leq\int_{0}^{\infty}\gamma^{t}R(\mathbb{E}[d(\Phi_{t}(x),\Phi_{t}(x^{\prime}))])dt.$
(6)

We introduce the factor $\ln(1/\gamma)/\ln(1/\gamma)$ and notice that
$\int_{0}^{\infty}\ln(1/\gamma)\gamma^{t}dt=1$. Thus, we may apply Jensen’s
inequality a second time to bound $|v(x)-v(x^{\prime})|$ by:

$\displaystyle\frac{R\big{(}\ln(1/\gamma)\int_{0}^{\infty}\gamma^{t}\mathbb{E}[d(\Phi_{t}(x),\Phi_{t}(x^{\prime}))]dt\big{)}}{\ln(1/\gamma)}.$
(7)

Introducing

$\displaystyle K_{\delta}(d_{0})=\ln(1/\gamma)\big{(}\int_{0}^{\delta}(\gamma
L)^{t}d_{0}dt+\int_{\delta}^{\infty}\gamma^{t}Ddt\big{)},$ (8)

we find that for every $\delta\geq 0$:

$\displaystyle\ln(1/\gamma)\int_{0}^{\infty}\gamma^{t}\mathbb{E}[d(\Phi_{t}(x),\Phi_{t}(x^{\prime}))]dt\leq
K_{\delta}(d_{0})$ (9)

by bounding $\mathbb{E}[d(\Phi_{t}(x),\Phi_{t}(x^{\prime}))]$ by $L^{t}d_{0}$
in the interval from $0$ to $\delta$ and by $D$ in the interval from $\delta$
to $\infty$.

Since (8) is a bound for arbitrary $\delta$, we can find a $\delta$ that is
optimal, i.e., minimizes $K_{\delta}(d_{0})$: The term $L^{t}d_{0}$, as a
function of $t$, is strongly monotonically increasing if $L>1$, and thus it
must eventually exceed $D$. It is thus advantageous to bound
$\mathbb{E}[d(\Phi_{t}(x),\Phi_{t}(x^{\prime}))]$ by $D$ after this “crossing”
occurs and better to bound it by $L^{t}d_{0}$ before, and hence the optimal
choice for $\delta$ is

$L^{\delta}d_{0}=D\implies\delta=\ln(D/d_{0})/\ln L.$

Figure 2 shows the intuitive idea of this choice. The quantity $\delta$ must
be finite and non-negative since $D>d_{0}>0$ and $L>1$. Indeed, for this
particular choice,

$\displaystyle
L^{\delta}=e^{\ln(L)\delta}=D/d_{0}\quad\text{and}\quad\gamma^{\delta}=e^{\ln(\gamma)\delta}=D^{-\beta}d_{0}^{\beta}.$

Solving the integral for the case $L\neq 1/\gamma$:

$\displaystyle K_{\delta}(d_{0})$
$\displaystyle=\ln(1/\gamma)\frac{d_{0}((\gamma L)^{\delta}-1)}{\ln(\gamma
L)}+D\gamma^{\delta}$
$\displaystyle=\ln(1/\gamma)\frac{d_{0}(D^{1-\beta}d_{0}^{\beta-1}-1)}{\ln(\gamma
L)}+D^{1-\beta}d_{0}^{\beta}$
$\displaystyle=D^{1-\beta}d_{0}^{\beta}\big{(}1-\frac{\ln\gamma}{\ln(\gamma
L)}\big{)}+\frac{\ln\gamma}{\ln(\gamma L)}d_{0}$ $\displaystyle=\frac{\ln
L}{\ln(\gamma L)}D^{1-\beta}d_{0}^{\beta}+\frac{\ln\gamma}{\ln(\gamma
L)}d_{0}=K(d_{0}),$

where we have used $1-\ln\gamma/\ln(\gamma L)=\ln L/\ln(\gamma L)$. In the
case $L=1/\gamma$, the integral becomes:

$\displaystyle K_{\delta}(d_{0})$
$\displaystyle=\ln(1/\gamma)d_{0}\delta+D\gamma^{\delta}$
$\displaystyle=(\ln(D/d_{0})+1)d_{0}=K(d_{0}),$

using $\beta=1$. Substitute $K_{\delta}(d_{0})=K(d_{0})$ in (9) and apply
monotonicity of $R$ in (7) to obtain (5). ∎

######  Theorem 2.

Suppose 1 holds with $T=\mathbb{N}_{0}$ and $L>1$, and let $H$ be as in
Theorem 1. Then

$\displaystyle|v(x)-v(x^{\prime})|\leq\frac{\ln(1/\gamma)}{1-\gamma}H(d(x,x^{\prime})),\quad\forall
x,x^{\prime}\in X.$

###### Proof.

We reduce the problem to a continuous setting. Extend $\Phi$ from
$\mathbb{N}_{0}$ to $\mathbb{R}_{\geq 0}$ by defining
$\Phi_{t}^{c}(x):=\Phi_{\lfloor t\rfloor}(x)$ for $t\in\mathbb{R}_{\geq 0}$,
and suppose that the value function for $\Phi^{c}$ is given by $v^{c}$. Now,
for any $n\in\mathbb{N}_{0}$,

$\displaystyle\int_{n}^{n+1}\gamma^{t}\mathbb{E}[r(\Phi_{t}^{c}(x))]dt$
$\displaystyle=\mathbb{E}[r(\Phi_{n}(x))]\int_{n}^{n+1}\gamma^{t}dt$
$\displaystyle=\gamma^{n}\mathbb{E}[r(\Phi_{n}(x))]\frac{1-\gamma}{\ln(1/\gamma)}.$

Therefore,

$\displaystyle
v^{c}(x)=\sum_{n=0}^{\infty}\int_{n}^{n+1}\gamma^{t}\mathbb{E}[r(\Phi_{t}^{c}(x))]dt=\frac{1-\gamma}{\ln(1/\gamma)}v(x).$
(10)

Moreover, $\Phi^{c}$ has the same LE-constant as $\Phi$: Since $L>1$,

$\displaystyle\mathbb{E}[d(\Phi_{t}^{c}(x),\Phi_{t}^{c}(x^{\prime}))]$
$\displaystyle=\mathbb{E}[d(\Phi_{\lfloor t\rfloor}(x),\Phi_{\lfloor
t\rfloor}(x^{\prime}))]$ $\displaystyle\leq L^{\lfloor
t\rfloor}d(x,x^{\prime})\leq L^{t}d(x,x^{\prime}).$

Then use (10) and Theorem 1 to bound $|v(x)-v(x^{\prime})|$ by:

$\displaystyle\frac{\ln(1/\gamma)}{1-\gamma}|v^{c}(x)-v^{c}(x^{\prime})|\leq\frac{\ln(1/\gamma)}{1-\gamma}H(d(x,x^{\prime})).$

∎

###  III-A Hölder and Lipschitz continuity

The results from [2] and [4] are entailed by Theorems 1 and 2, but we can also
extend them into the stochastic setting: Under the assumptions of 1 and if the
reward function is Lipschitz continuous, i.e., $R(d_{0})=Cd_{0}$ for a $C\geq
0$, then there exist $\beta\in(0,1]$ and $A\geq 0$ such that:

$\displaystyle|v(x)-v(x^{\prime})|\leq Ad(x,x^{\prime})^{\beta},\quad\forall
x,x^{\prime}\in S.$ (11)

In fact, if

  * $L<1/\gamma$,

then $\beta=1$ (Corollary 1),

  * $L=1/\gamma$,

then any $\beta\in(0,1)$ is allowed (Corollary 2),

  * $L>1/\gamma$,

then $\beta=\ln(1/\gamma)/\ln L$ (Corollary 3).

The first case shows Lipschitz continuity, the second case Hölder continuity
for any exponent in $(0,1)$, and the last case Hölder continuity for the
exponent $\ln(1/\gamma)/\ln L$. The results follow from Theorems 1 and 2,
except for the straight-forward case $L<1/\gamma$. For $L=1/\gamma$, one can
take $\beta$ as close to $1$ as one desires, but the corresponding $A$ becomes
arbitrarily large. For the proofs, we refer to the appendix.

##  IV Sharpness and experiments

While Theorems 1 and 2 provide upper bounds on the modulus of $v$, it is still
unclear how good the bounds are. Since they only depend on $D$, $L$, $\gamma$
and $R$, our objective is now to construct a random system that is LE-
continuous with constant $L$ and whose corresponding value function has the
largest possible modulus. The difference between the resulting modulus and the
proposed bounds then indicates how much the latter could be decreased. We
start by looking at the deterministic map $x\mapsto x^{L}$ defined on $[0,1]$.
Afterward, we show that one can construct a random system whose value function
has a modulus that exactly matches the one proposed in Theorem 1.

###  IV-A An example for a steep value function

Take an arbitrary $L>1$ and define a dynamical system
$\Phi:[0,1]\rightarrow[0,1]$ via $\Phi(x)=x^{L}.$ For $n\in\mathbb{N}_{0}$,
set

$\Phi(n,x):=\Phi_{n}(x):=\Phi^{n}(x)=x^{(L^{n})}.$

The system $\Phi$ has a stable fixed point at $x=0$ and an unstable one at
$x=1$. We also have

$\partial_{x}\Phi_{n}(x)=L^{n}x^{L^{n}-1}\leq L^{n},$

which shows that $\Phi_{n}$ is Lipschitz continuous with constant $L^{n}$ for
the Euclidean metric $d(x,x^{\prime})=|x-x^{\prime}|$. If interpreted as a
random system, i.e., when $\Phi_{n}(x,\omega):=\Phi_{n}(x)$ for
$\omega\in\Omega$, it must be LE-continuous with constant $L$. For $r(x)=x$
and any discount factor $\gamma\in(0,1)$, the value function is

$\displaystyle v(x)=\sum_{n=0}^{\infty}\gamma^{n}x^{(L^{n})}.$

It is maximal for $x=1$, where it takes the value $1/(1-\gamma)$, and minimal
for $x=0$, where it is $0$. All other $x\in(0,1)$ quickly approach $0$ as $n$
increases. One can see the value functions for varying discount factors on the
left in Figure 3, depicting a steep increase as $x$ approaches $1$ for the
case $L\geq 1/\gamma$. If that happens, $v$’s derivative is unbounded at
$x=1$. The modulus takes the form $W_{v}(d_{0})=v(1)-v(1-d_{0}),$ see Lemma 1
in the appendix for details. As observable on the right of Figure 3, the
bounds provided by Theorem 2 are close to $W_{v}$ for varying values of
$\gamma$. Bounds that only assume knowledge about $L$ and $\gamma$ can hence
not be much better than Theorem 2 – at least for linear reward functions.

Figure 3: Left: The value functions corresponding to the example in Section IV
for $L=1.5$ and discount factors
$(\gamma_{0},\gamma_{1},\gamma_{2})=(0.5,0.9,0.99)$ when normalized to a
maximal value of $1$. Right: The moduli of continuity for the value functions
in comparison to the bounds given by Theorem 2, visualized by dashed lines in
the same color. The bounds from [4] are visualized using dotted lines (also
the same color).

###  IV-B Sharpness

In the time-continuous setting and for linear rewards, we cannot improve the
bound:

######  Proposition 3.

Let $D>0$, $L>1$, $\gamma\in(0,1)$ and suppose $r(x)=Ax$. Then there exists a
continuous-time random system with LE-constant $L$ satisfying 1 such that
$v$’s modulus is exactly $H$ from Theorem 1.

###### Proof.

Define a deterministic system on $S=[0,D]$ by
$\Phi_{t}(x,\omega)=\min\\{D,L^{t}x\\}$ and let
$d(x,x^{\prime})=|x-x^{\prime}|$. One may verify that $\Phi$ is indeed LE-
continuous with constant $L$. The value function takes the form

$\displaystyle v(x)=A\big{(}\int^{\delta}_{0}(\gamma
L)^{t}xdt+\int_{\delta}^{\infty}\gamma^{t}Ddt\big{)},$

where $\delta$ solves $L^{\delta}x=D$, i.e. $\delta=\ln(D/x)/\ln L$ for $x>0$
and $\delta=\infty$ otherwise. This is exactly

$v(x)=AK(x)/\ln(1/\gamma),$

where $K$ is from Theorem 1, see (8) and (9) with $K_{\delta}(d_{0})=K(d_{0})$
for $\delta=\ln(D/d_{0})/\ln L$. Let now $x,x^{\prime}\in S$. Then

$\displaystyle|v(x)-v(x^{\prime})|=\frac{|A|}{\ln(1/\gamma)}|K(x)-K(x^{\prime})|.$

Moreover, notice for $L\neq 1/\gamma$ and $D\geq x>0$:

$\displaystyle K^{\prime}(x)$ $\displaystyle=\frac{\ln\gamma}{\ln(\gamma
L)}(1-D^{1-\beta}x^{\beta-1})\geq 0,$ $\displaystyle K^{\prime\prime}(x)$
$\displaystyle=-\frac{\ln\gamma}{\ln(\gamma
L)}D^{1-\beta}(\beta-1)x^{\beta-2}\leq 0,$

which is due to

$\displaystyle\mathrm{sign}(\frac{\ln\gamma}{\ln(\gamma
L)})=\mathrm{sign}(\beta-1)=\mathrm{sign}(1-(\frac{x}{D})^{\beta-1}),$

for both $L<1/\gamma$ and $L>1/\gamma$. Thus, $K$ is monotonically increasing,
while $K^{\prime}$ is monotonically decreasing. Without loss of generality for
$x\geq x^{\prime}$,

$\displaystyle|K(x)-K(x^{\prime})|$ $\displaystyle=K(x)-K(x^{\prime})$
$\displaystyle=\textstyle\int_{x^{\prime}}^{x}K^{\prime}(y)dy$
$\displaystyle\leq\textstyle\int_{x^{\prime}}^{x}K^{\prime}(y-x^{\prime})dy$
$\displaystyle=\textstyle\int_{0}^{x-x^{\prime}}K^{\prime}(y)dy$
$\displaystyle=K(x-x^{\prime})-K(0)=K(x-x^{\prime}).$

Since the modulus is defined as the supremum over all $x,x^{\prime}$ and since
$K$ is monotonically increasing, the best way to maximize $K(x-x^{\prime})$
subject to $0\leq x-x^{\prime}\leq d_{0}$ is to put $x=d_{0}$ and
$x^{\prime}=0$:

$\displaystyle W_{v}(d_{0})=\sup_{0\leq x-x^{\prime}\leq
d_{0}}\frac{|A|K(x-x^{\prime})}{\ln(1/\gamma)}=H(d_{0}),$

which completes the proof for $L\neq 1/\gamma$. For the case $L=1/\gamma$, we
have for $x>0$:

$\displaystyle K^{\prime}(x)=\ln(D/x)\geq 0,\quad K^{\prime\prime}(x)=-1/x\leq
0.$

The remainder of the argument is analogous. ∎

##  V Disturbance implies Differentiability

If one stays in the deterministic realm, then one has to live with systems
that give rise to non-differentiable value functions. However: One can
_disturb_ such systems by adding a little bit of Gaussian noise in every step.
This per-step smoothing procedure suffices to obtain a differentiable value
function; we illustrate this for the case $S=\mathbb{R}^{n}$:

######  Proposition 4.

Let $\Phi:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$,
$r:\mathbb{R}^{n}\rightarrow\mathbb{R}$ be a bounded and differentiable reward
function, $\gamma\in(0,1)$ be a discount factor and $\sigma^{2}>0$ be a
variance. Define a random system

$\displaystyle\psi_{n+1}(x)=\psi_{n}(\Phi(x+\xi_{n})),\quad\psi_{0}(x)=x$

with $\xi_{0},\xi_{1},\dots\sim\mathcal{N}(0,\sigma^{2}I)$ iid and $I$ being
the $n\times n$ identity matrix. Then the value function $w$ of $\psi$ is
differentiable, and its gradient is given for
$\xi\sim\mathcal{N}(0,\sigma^{2}I)$ by

$\displaystyle\nabla w(x)=\nabla
r(x)+\frac{\gamma}{\sigma^{2}}\mathbb{E}[w(\Phi(x+\xi))\cdot\xi].$ (12)

###### Proof.

Because of the independence of $\xi_{0},\xi_{1},\dots$, the new value function
satisfies

$w(x)=r(x)+\gamma\mathbb{E}[w(\Phi(x+\xi))].$

Writing the expectation as a convolution,

$\mathbb{E}[w(\Phi(x+\xi))]=\mathbb{E}[w(\Phi(x-\xi))]=((w\circ\Phi)*p)(x),$

where $p:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is the density of
$\mathcal{N}(0,\sigma^{2}I)$. We can take the gradient with respect to $p$
_instead_ of $w\circ\Phi$,

$\displaystyle\nabla((w\circ\Phi)*p)(x)=((w\circ\Phi)*\nabla p)(x).$

Using $\nabla p(\xi)=-\frac{1}{\sigma^{2}}p(\xi)\cdot\xi$,

$\displaystyle\nabla\mathbb{E}[w(\Phi(x+\xi))]=\frac{1}{\sigma^{2}}\mathbb{E}[w(\Phi(x-\xi))\cdot(-\xi)],$

which gives (12) after changing variables from $-\xi$ to $\xi$. The gradient
$\nabla w(x)$ must also be finite, which follows from the fact that $\xi$ has
finite variance and from the boundedness of $w$ (by
$\sup_{x}|r(x)|/(1-\gamma)$, to be precise). ∎

In Figure 1, we have extended the logistic map’s domain from $[0,1]$ to
$\mathbb{R}$ by projecting onto $[0,1]$ in every step. We then added Gaussian
noise and obtained the value function. Proposition 4 shows that it must be
differentiable.

##  VI Conclusion

We investigated the continuity of the value function in reinforcement learning
and optimal control. In particular, we have shown that the value function can
be nowhere differentiable, as suggested by Proposition 1, even if the reward
function and the underlying system are both relatively well-behaved. Moreover,
we have derived an upper bound to the value function’s modulus of continuity
on bounded state spaces in both discrete and continuous time and have shown
sharpness of this bound for linear reward functions. As a consequence, we
argued that the value function is Hölder continuous, in line with [2, 4],
which is useful for bounding variances. Finally, we have shown how to obtain
differentiable value functions by introducing “disturbances”.

## Acknowledgements

H.H. and S.P. acknowledge financial support from the project “SAIL:
SustAInable Life-cycle of Intelligent Socio-Technical Systems” (Grant ID
NW21-059D), which is funded by the program “Netzwerke 2021” of the Ministry of
Culture and Science of the State of Northrhine Westphalia, Germany.

######  Lemma 1.

Consider the example from Section IV. Then $W_{v}(d_{0})=v(1)-v(1-d_{0}).$

###### Proof.

Note that $v$ is monotonically increasing, which follows from the non-
negativity of $v^{\prime}(x)$. Thus,

$\displaystyle W_{v}(d_{0})=\sup_{d_{0}\geq x-x^{\prime}\geq
0}\\{v(x)-v(x^{\prime})\\}.$

Moreover, $v^{\prime}(x)$ is monotonically increasing, since

$\displaystyle v^{\prime\prime}(x)=\sum_{n=0}^{\infty}(\gamma
L)^{n}(L^{n}-1)x^{L^{n}-2}\geq 0.$

Hence, for $1\geq x\geq x^{\prime}\geq 0$,

$\displaystyle v(x)-v(x^{\prime})$
$\displaystyle=\textstyle\int_{x^{\prime}}^{x}v^{\prime}(y)dy$
$\displaystyle\leq\textstyle\int_{x^{\prime}}^{x}v^{\prime}(y+(1-x))dy$
$\displaystyle=\textstyle\int_{x^{\prime}+(1-x)}^{1}v^{\prime}(y)dy$
$\displaystyle=v(1)-v(x^{\prime}+(1-x))$
$\displaystyle=v(1)-v(1-(x-x^{\prime})).$

Introducing $\delta=x-x^{\prime}$, the modulus becomes

$W_{v}(d_{0})=v(1)-\inf_{d_{0}\geq\delta\geq 0}v(1-\delta).$

Since $v$ is monotonically increasing, the $\delta$ that minimizes
$v(1-\delta)$ should be as large as possible. Since it is bounded by $d_{0}$,
the optimal choice becomes $\delta=d_{0}$. ∎

######  Corollary 1 (Case $<$).

Suppose 1 holds with $L<1/\gamma$ and $R(d_{0})=Cd_{0}$. Then (11) is
satisfied for $\beta=1$ and

$\displaystyle A=C\cdot\begin{cases}\frac{1}{\ln(1/(\gamma
L))}&T=\mathbb{R}_{\geq 0},\\\ \frac{1}{1-\gamma
L}&T=\mathbb{N}_{0}.\end{cases}$

###### Proof.

Set $d_{0}=d(x,x^{\prime})$. For $T=\mathbb{N}_{0}$,

$\displaystyle|v(x)-v(x^{\prime})|$
$\displaystyle\leq\sum_{n=0}^{\infty}\gamma^{n}CL^{n}d_{0}=\frac{C}{1-\gamma
L}d_{0}.$

Similarly, for $T=\mathbb{R}_{\geq 0}$,

$\displaystyle|v(x)-v(x^{\prime})|$
$\displaystyle\leq\int_{0}^{\infty}\gamma^{t}CL^{t}d_{0}dt=\frac{C}{\ln(1/(\gamma
L))}d_{0}.$

∎

######  Corollary 2 (Case $=$).

Suppose 1 holds with $L=1/\gamma$ and $R(d_{0})=Cd_{0}$. For $\beta\in(0,1)$,
(11) is satisfied by

$\displaystyle A=C\big{(}\textstyle\frac{1}{e(1-\beta)}+D^{1-\beta}(\ln
D+1)\big{)}\cdot\begin{cases}\frac{1}{\ln(1/\gamma)}&T=\mathbb{R}_{\geq 0},\\\
\frac{1}{1-\gamma}&T=\mathbb{N}_{0}.\end{cases}$

###### Proof.

For $T=\mathbb{R}_{\geq 0}$, rewrite $H$ from Theorem 1 as

$H(d_{0})=\frac{C}{\ln(1/\gamma)}d_{0}\ln(1/d_{0})+\frac{C(\ln
D+1)}{\ln(1/\gamma)}d_{0}.$ (13)

In the last term and for any $\beta\in(0,1)$, we can bound $d_{0}$ by
$d_{0}=d_{0}^{1-\beta}d_{0}^{\beta}\leq D^{1-\beta}d_{0}^{\beta}$. The
challenging term is $d_{0}\ln(1/d_{0})$, for which we would like to find a
$B\geq 0$ such that $Bd_{0}^{\beta}\geq d_{0}\ln(1/d_{0})$ for every
$d_{0}\in(0,\infty)$. Divide by $d_{0}^{\beta}$ to bound $B$ from below by

$h(d_{0}):=d_{0}^{1-\beta}\ln(1/d_{0}).$

Notice that

$\displaystyle h^{\prime}(d_{0})=((1-\beta)\ln(1/d_{0})-1)/d_{0}^{\beta},$

with

$\displaystyle h^{\prime}(d_{0})\geq 0$ $\displaystyle\iff-\ln d_{0}\geq
1/(1-\beta)$ $\displaystyle\iff d_{0}\leq e^{-1/(1-\beta)}$

and $h^{\prime}(d_{0}^{*})=0$ for $d_{0}^{*}=e^{-1/(1-\beta)}$. Indeed,
$h^{\prime}(d_{0})\geq 0$ for $d_{0}\in(0,d_{0}^{*}]$ and
$h^{\prime}(d_{0})\leq 0$ for $d_{0}\in[d_{0}^{*},\infty)$. Hence, $d_{0}^{*}$
is a global maximum of $h$, and thus

$\displaystyle B:=\frac{1}{e(1-\beta)}=h(d_{0}^{*})\geq h(d_{0})\quad\forall
d_{0}\in(0,\infty).$

In (13), we can thus bound $d_{0}\ln(1/d_{0})$ by $Bd_{0}^{\beta}$ and $d_{0}$
by $D^{1-\beta}d_{0}^{\beta}$. After factoring out $d_{0}^{\beta}$, we obtain

$H(d_{0})\leq\frac{C}{\ln(1/\gamma)}\big{(}\frac{1}{e(1-\beta)}+D^{1-\beta}(\ln
D+1)\big{)}d_{0}^{\beta},$

which shows the case for $T=\mathbb{R}_{\geq 0}$. For $T=\mathbb{N}_{0}$, use
Theorem 2 and then bound $H$ in the same way. ∎

######  Corollary 3 (Case $>$).

Suppose 1 holds with $L>1/\gamma$ and $R(d_{0})=Cd_{0}$. Then (11) is
satisfied when setting $\beta=\ln(1/\gamma)/\ln L$ and

$\displaystyle A=\frac{CD^{1-\beta}\ln L}{\ln(\gamma
L)}\cdot\begin{cases}\frac{1}{\ln(1/\gamma)}&T=\mathbb{R}_{\geq 0},\\\
\frac{1}{1-\gamma}&T=\mathbb{N}_{0}.\end{cases}$

###### Proof.

For $T=\mathbb{R}_{\geq 0}$, apply Theorem 1 and notice that the term
$(\ln\gamma/\ln(\gamma L))d_{0}$ in $K(d_{0})$ is non-positive. Relaxing the
bound by removing this term gives $A$. Additionally, note that $\beta$ is
strictly between 0 and 1 since $L>1/\gamma>1$ implies $\ln L>\ln(1/\gamma)>0$.
For $T=\mathbb{N}_{0}$, we use the same argument but factor in
$\ln(1/\gamma)/(1-\gamma)$ using Theorem 2. ∎

## References

  * [1] Ludwig Arnold.  Random Dynamical Systems, volume 1609, pages 1–43.  Springer Berlin Heidelberg, Berlin, Heidelberg, 1995.
  * [2] Martino Bardi and Italo Capuzzo-Dolcetta.  Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations.  Birkhäuser Boston, Boston, MA, 1997.
  * [3] Vincenzo Basco and Hélène Frankowska.  Lipschitz Continuity of the Value Function for the Infinite Horizon Optimal Control Problem Under State Constraints.  In Fatiha Alabau-Boussouira, Fabio Ancona, Alessio Porretta, and Carlo Sinestrari, editors, Trends in Control Theory and Partial Differential Equations, volume 32, pages 17–38. Springer International Publishing, Cham, 2019.
  * [4] Andrey Bernstein and Nahum Shimkin.  Adaptive-resolution reinforcement learning with polynomial exploration in deterministic domains.  Machine Learning, 81(3):359–397, December 2010.
  * [5] Justin A. Boyan.  Technical Update: Least-Squares Temporal Difference Learning.  Machine Learning, 49(2/3):233–246, 2002.
  * [6] A A Dorogovtsev and I. I. Nishchenko.  An analysis of stochastic flows.  Communications on Stochastic Analysis, 8(3), September 2014.
  * [7] Andrey A. Dorogovtsev.  Measure-Valued Processes and Stochastic Flows.  De Gruyter, October 2023.
  * [8] Minyi Huang.  Uniqueness of Constrained Viscosity Solutions in Hybrid Control Systems.  SIAM Journal on Control and Optimization, 46(1):332–355, January 2007.
  * [9] Kaito Ito, Takuya Ikeda, and Kenji Kashima.  Continuity of the Value Function for Stochastic Sparse Optimal Control.  IFAC-PapersOnLine, 53(2):7179–7184, 2020.
  * [10] Timothy P. Lillicrap, Jonathan J. Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa, David Silver, and Daan Wierstra.  Continuous control with deep reinforcement learning.  In Yoshua Bengio and Yann LeCun, editors, 4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings, 2016.
  * [11] Juan Pablo Rincón-Zapatero and Manuel S. Santos.  Differentiability of the value function without interiority assumptions.  Journal of Economic Theory, 144(5):1948–1964, September 2009.
  * [12] David Silver, Guy Lever, Nicolas Heess, Thomas Degris, Daan Wierstra, and Martin Riedmiller.  Deterministic Policy Gradient Algorithms.  In Proceedings of the 31st International Conference on International Conference on Machine Learning, volume 32, pages I–387–I–395, Beijing, China, 2014. JMLR.org.
  * [13] Bruno Strulovici and Martin Szydlowski.  On the smoothness of value functions and the existence of optimal strategies in diffusion models.  Journal of Economic Theory, 159:1016–1055, September 2015.
  * [14] Richard S. Sutton and Andrew G. Barto.  Reinforcement Learning: An Introduction.  Adaptive Computation and Machine Learning. MIT Press, Cambridge, Mass, 1998.
  * [15] Richard S. Sutton, Hamid Reza Maei, Doina Precup, Shalabh Bhatnagar, David Silver, Csaba Szepesvári, and Eric Wiewiora.  Fast gradient-descent methods for temporal-difference learning with linear function approximation.  In Proceedings of the 26th Annual International Conference on Machine Learning, pages 993–1000, Montreal Quebec Canada, June 2009. ACM.
  * [16] J.N. Tsitsiklis and B. Van Roy.  An analysis of temporal-difference learning with function approximation.  IEEE Transactions on Automatic Control, 42(5):674–690, May 1997\.
  * [17] Masaya Yamaguti and Masayoshi Hata.  Weierstrass’s function and chaos.  Hokkaido Mathematical Journal, 12(3), October 1983.


# Mapping illegal waste dumping sites with neural-network classification of
satellite imagery

Maria Roberta Devesa  Dymaxion LabsBuenos AiresArgentina
<EMAIL_ADDRESS>and  Antonio Vazquez Brust  Fundación Bunge y
BornBuenos AiresArgentina<EMAIL_ADDRESS>
(2021; 2021)

###### Abstract.

Public health and habitat quality are crucial goals of urban planning. In
recent years, the severe social and environmental impact of illegal waste
dumping sites has made them one of the most serious problems faced by cities
in the Global South, in a context of scarce information available for decision
making. To help identify the location of dumping sites and track their
evolution over time we adopt a data-driven model from the machine learning
domain, analyzing satellite images. This allows us to take advantage of the
increasing availability of geo-spatial open-data, high-resolution satellite
imagery, and open source tools to train machine learning algorithms with a
small set of known waste dumping sites in Buenos Aires, and then predict the
location of other sites over vast areas at high speed and low cost. This case
study shows the results of a collaboration between Dymaxion Labs and Fundación
Bunge y Born to harness this technique in order to create a comprehensive map
of potential locations of illegal waste dumping sites in the region.

Neural networks, dumping sites, urban waste management, transfer learning, GIS

††copyright: acmcopyright††copyright: rightsretained††journalyear: 2021††isbn:
††journalyear: 2021††conference: ACM SIGKDD Conference on Knowledge Discovery
and Data Mining ; August 14-18, 2021; Virtual Conference††booktitle: KDD
Humanitarian Mapping Wokshop ’21: ACM SIGKDD Conference on Knowledge Discovery
and Data Mining August 14-18, 2021; Virtual Conference

##  1\. Introduction

Waste management is a a critical utility, essential for urban life. It is also
a significant challenge for municipal governments all over the world. Often
rated among the top priorities faced by cities in developing countries, it can
also be taken for granted, missing from national or international political
agendas (Wilson and Velis, 2015). The problem is aggravated in developing
regions such as Latin America where the rapidly increasing production of
waste, coupled with burdened municipal budgets and scarce national attention,
result in poor waste management (Hoornweg and Giannelli, 2007). Dumping sites,
areas where intentionally and illegally abandoned waste accumulates, are a
physical manifestation of this inadequacy. Dumping sites generate severe
environmental degradation that disproportionately affects the most vulnerable
population (Pellow, 2004), specially those who -without access to better
options- occupy land near or even inside dumping sites. This close proximity
of human settlement to informal waste disposal areas has been associated to
severe health risks such as respiratory disease (Ibrahim et al., 2021) and
heavy metal poisoning (Cittadino et al., 2020). Due to their unplanned and
illicit nature, both the appearance and the growth of dumping sites are hard
to monitor and can evolve unnoticed by governments. This lack of visibility is
part of a larger issue, as precise data on the many aspects of waste
management is severely lacking all over the world. In 2015 an international
team tasked by the UN to measure the societal and environmental impact of poor
waste management found that ”availability and reliability of waste and
resource data is dire, and urgently needs attention” (Wilson and Velis, 2015).
As with other pressing developing world issues where a chronic absence of
timely information hinders decision-making, the exploration of remote sensing
approaches has the potential to yield useful and frequently updated results at
a fraction of the cost associated with on-the-ground data collection.

##  2\. Related Work

As far back as 1998 Algarni and Elsadiq (Ali Algarni and Elsadiq Ali, 1998)
reported their experience identifying waste disposal sites in the City of
Riyadh (Saudi Arabia) and its surroundings, applying digital processing
techniques to SPOT imagery to facilitate visual detection by experts using GIS
software. GIS-aided classification methods, taking advantage of spatial
modeling features from widely available GUI based software, continue to be
developed. Gill et al. (Gill et al., 2019) used Landsat thermal data to detect
unauthorised landfills in Kuwait by their heat generation, measuring ground
surface temperature and isolating suspiciously hot areas; Vambol et al.
(Vambol et al., 2019) experimented with a GIS-based workflow to detect illegal
waste dumps in the Kharkiv region (Ukraine) by analyzing the pixel brightness
values of high resolution satellite imagery obtained via Google Earth. On the
other hand, as expected given the fast development of general-purpose
classification algorithms during the last decades, recent experiences on
dumping detection aim to scale beyond traditional spatial modeling by relying
on machine learning techniques. Akinina et al. (Akinina et al., 2017) compared
the performance of three different algorithms applied to illegal dump
detection using Landsat data for the Ryazan Region in central Russia: Parzen
window, AdaBoost and neural network classification, with the latter yielding
the best results both in accuracy and speed. Skogsmo (Skogsmo, 2020) proposed
a support vector machine (SVM) approach for large-scale analysis of Sentinel
data, along with a case study of automated detection of waste dumps in
Kampala, Uganda. Due to our use of machine learning rather than a GIS
approach, our work most closely resembles the last two examples. Compared to
those, a key difference is our implementation of a deep learning model instead
of alternative classification algorithms like SVM or decision trees.

##  3\. Materials

###  3.1. Study area

This study focuses on the Buenos Aires Province in Argentina. The province has
an area of 307,571 km², although most of its population is concentrated on its
largest conurbation, the Buenos Aires Metropolitan Area (occupying 3,830 km²).
During the last decades urban solid waste production grew rapidly in the
region, overflowing the capacity of publicly run waste disposal centers and
propitiating the multiplication of illegal open air waste dumping sites. The
environmental and societal impact associated with the proliferation of dumping
sites made them one of the most conflictive patterns of land use change in the
region, and a significant risk for public health (Cittadino et al., 2012).
This situation makes Buenos Aires a relevant area for experimentation with
large-scale, high-frequency detection of waste dumping sites.

###  3.2. Data Sources

Satellite images used in this research are from the Sentinel-2 Earth
observation mission, retrieved via Google Earth Engine. The image resolution
used is 10 meters, and they were captured in the 2021-02-01 - 2021-02-28
period. After bench-marking results with different combinations, the best
performance was obtained using the three color bands R-G-B, the near-infrared
(IR), and the 2 near-infrared short wave bands (SWIR-1 and SWIR-2).

Two data sets with georeferenced locations of known illegal waste dumps were
available: A survey from the Autoridad de Cuenca Matanza Riachuelo (a public
agency tasked with improving the environmental quality of the region’s main
river basin), and a report from the Buenos Aires Province ombudsman. Together
they provide 81 annotations for open air waste dumps, which were used to train
the supervised machine learning algorithm presented in this study.

##  4\. Method

This study implements a Convolutional Neural Network (CNN) model by a U-Net
architecture developed for a segmentation process that identifies the
characteristic patterns of dumping sites. The output of the model will be an
image of the same size as the input, in which each pixel is the probability of
belonging to a waste dump. This means that the model will find not only their
position but also their boundaries. This feature is of great importance for
planning applications, as it allows measuring the growth of sites over time,
as well as their potential intersection with other land uses (watercourses,
inhabited areas, etc). The model needs an image set and annotations for the
training process. Each image has an associated binary mask that delimits the
object of interest. The image can have many bands (besides just RGB) so
additional information can be included as a band. Figure 1 shows a diagram of
the training process. This study uses 6-band images: three color bands (RGB),
a near-infrared (NIR) band, a short wave infrared (SWIR-1) band, and the
normalized difference of the two SWIR. Previous studies of methane emission
detection inspired the use of combined SWIR-1 and SWIR-2 bands, on the
assumption that open-air waste dumps may emit a significant amount of methane,
particularly from organic waste. NIR and SWIR bands are used because they are
more sensitive to green areas. As large dumping sites are usually located in
rural areas, these bands may help to classify the background better. Table 1
shows the improvement when these bands are included.

The following subsection details each step of the process.

Figure 1. Diagram of the training process

###  4.1. Pre-processing

The preprocessing stage aims to create the dataset used to train the model and
predict with it. A combination of tools from GDAL and pysatproc were used to
create image sets with the chosen bands and masks for each one. In addition,
the pysatproc tool sets the image size for the model and allows to perform a
sliding window over an annotation to increase the number of images available
for training. As an example, the Figure 2 shows an image with its mask.

Figure 2. A example of an image (left) and its respective mask (right) used to
train the model

###  4.2. The model

As it was briefly mentioned before, the U-Net is a type of CNN model. It was
initially created for biomedical image segmentation (Olaf Ronneberger, [n.d.])
although it also works for the segmentation of natural images. The U-Net is
constructed by multiple down-sampling and up-sampling layers, which are
connected by concatenating layers. This way it combines the high precision
activation nodes with nodes of the previous layer with better resolution.

###  4.3. The weighted binary cross entropy method

Under imbalanced data, the model tends to recognize background more than
objects of interest, thus learning to predict background terrain more than
objects. To deal with this condition the model uses weighted binary cross-
entropy (et al., [n.d.]). It allows us to consider different weight values for
each class (or one class, in this case), to force the loss function to focus
more on the cases that are less frequent on the dataset.

###  4.4. Post-processing

The prediction process produces a set of images with multiple bands as
classes, where each one has the probability of finding the object of interest
in the image, pixel by pixel. Our case only considers one class, so the output
consists of one band images.

To remove cases with low probability a filter process is implemented over all
predicted results, applying per-pixel probability threshold. Additionally, the
resulting images were polygonized to remove predicted objects based on surface
area, as the smaller ones are more likely to be false positives. The threshold
was set at 100 m², as waste dumping sites are usually larger than that.

##  5\. Results

A U-net model was performed to segment pixels showing dumping sites from
satellite images. The U-net architecture was built with four down-sampling
layers and four up-sampling layers, with six channels in the first
convolutional layer, and a 3×3 kernel size in all of the model’s convolutional
layers.

The dataset contains 1917 images with a size of 100x100 pixels per 6-band of
depth images and their respective masks. It was then divided into training,
validation, and test. The first two were used to perform the model and the
last one to obtain the metrics over untouched data. A split ratio of 0.1 for
testing, and a split ratio of 0.7 for training and validation. The metric
analyzed in this study is the Intersection over Union (IoU), which is mainly
used in applications related to object detection and object segmentation. It
is essentially a method to quantify the percent overlap between the target
mask and our prediction output, measuring the number of pixels common between
the target and prediction masks divided by the total number of pixels present
across both masks. The model was set with batch 16 and epoch 30, but after 18
epochs the model reaches the plateau; with a 0.0458 loss and a mean IoU of
0.675.

Table-1 shows the comparison of the model with different approaches. The loss
and the IoU values reached over the test sample are compared for a different
combination of bands on the input images. It shows that the optimal
performance is obtained using RGB-NIR-SWIR-NDSW

Image information | Loss | IoU
---|---|---
RGB | 0.0613 | 0.593
RGB-NIR | 0.0712 | 0.660
RGB-NIR-SWIR | 0.0652 | 0.650
RGB-NIR-SWIR-NDSW | 0.0458 | 0.675
Table 1. Table of the metrics for each approach.

As an example Figure 3 shows some model result examples, obtained after
processing images from 594 urban centers across the Buenos Aires province. The
images were originally captured between 2021-02-01 and 2021-02-28.

Figure 3. Results for different locations around Buenos Aires city. The red
line shows the area detected by the model as a potential illegal waste dumping
site.

Some of the main limitations found after running the prediction process over
the different zones of Buenos Aires were that despite that the model detected
several waste dumps, also other places like home constructions on the initial
state were incorrectly tagged.

In addition, small waste dumps (size less than 1 hectare) are hard to detect
due to the image resolution. And in some areas may be hard to get images
without clouds.

##  6\. Conclusions and future work

This research proposes a supervised machine learning method to detect illegal
waste dumping sites using satellite imagery. It consists of a U-Net
architecture based on CNN for classification and segmentation and weighted
cross-entropy. The model’s performance in predicting IoU is 0.6304 and it
reaches a loss of 0.0413.

To get further insights into model performance, future work can include images
from earlier years and additional annotations to train the model and predict
with it.

From a replication perspective, working with Sentinel- 2 satellite images is
critical not only for the current study but also for future iterations because
of their public availability and ease of acquisition. Free access to 10 meter
resolution satellite imagery paired with state-of-the-art open source machine
learning algorithms enable the fast analysis of extensive areas at a minimal
cost.

To facilitate monitoring growth of illegal dumping sites, we are currently
packaging our methodology and algorithms as cloud-based programming notebook,
and hosting hands-on training sessions aimed at municipal personnel.

## References

  * (1)
  * Akinina et al. (2017) Natalia V. Akinina, Maxim V. Akinin, Alexandr I. Taganov, and Michael B. Nikiforov. 2017\.  Methods of detection in satellite images of illegal dumps by using a method based on tree classifier. In _2017 6th Mediterranean Conference on Embedded Computing (MECO)_. IEEE, Bar, 1–3.  https://doi.org/10.1109/MECO.2017.7977179
  * Ali Algarni and Elsadiq Ali (1998) Dafer Ali Algarni and Abdalla Elsadiq Ali. 1998.  Mapping Waste-disposal Sites Using SPOT Remote Sensor Data: Riyadh Case.  _Journal of King Saud University - Engineering Sciences_ 10, 1 (1998), 15–29.  https://doi.org/10.1016/S1018-3639(18)30685-8
  * Cittadino et al. (2012) Alejandro Cittadino, V. D’hers, M.A. Igarzábal de Nistal, M. V. Majul, N. Ocello, and J. Zamorano. 2012.  _Atlas de la basura del Área Metropolitana de Buenos Aires_.  Wolkowicz Editores.
  * Cittadino et al. (2020) Alejandro Cittadino, Natalia Ocello, Maria Victoria Majul, Raquel Ajhuacho, Patricia Dietrich, and Maria Adela Igarzabal. 2020.  Heavy metal pollution and health risk assessment of soils from open dumps in the Metropolitan Area of Buenos Aires, Argentina.  _Environmental Monitoring and Assessment_ 192, 5 (May 2020), 291.  https://doi.org/10.1007/s10661-020-8246-x
  * et al. ([n.d.]) Fidel A. Guerrero-Pena et al. [n.d.].  Multiclass Weighted Loss for Instance Segmentation of Cluttered Cells.  _arxiv_ ([n. d.]).  https://arxiv.org/abs/1802.07465
  * Gill et al. (2019) Jasravia Gill, Kamil Faisal, Ahmed Shaker, and Wai Yeung Yan. 2019.  Detection of waste dumping locations in landfill using multi-temporal Landsat thermal images.  _Waste Management & Research: The Journal for a Sustainable Circular Economy_ 37, 4 (April 2019), 386–393.  https://doi.org/10.1177/0734242X18821808
  * Hoornweg and Giannelli (2007) Daniel Hoornweg and Natalie Giannelli. 2007.  _Managing Municipal Solid Waste in Latin America and the Caribbean : Integrating the Private Sector, Harnessing Incentives_.  Technical Report 28. World Bank.  https://openknowledge.worldbank.org/handle/10986/10639
  * Ibrahim et al. (2021) Mohd Faiz Ibrahim, Rozita Hod, Haidar Rizal Toha, Azmawati Mohammed Nawi, Idayu Badilla Idris, Hanizah Mohd Yusoff, and Mazrura Sahani. 2021.  The Impacts of Illegal Toxic Waste Dumping on Children’s Health: A Review and Case Study from Pasir Gudang, Malaysia.  _International Journal of Environmental Research and Public Health_ 18, 5 (Feb. 2021), 2221.  https://doi.org/10.3390/ijerph18052221
  * Olaf Ronneberger ([n.d.]) Thomas Brox Olaf Ronneberger, Philipp Fischer. [n.d.].  U-Net: Convolutional Networks for Biomedical Image Segmentation.  ([n. d.]).  https://arxiv.org/abs/1505.04597
  * Pellow (2004) David N. Pellow. 2004\.  The Politics of Illegal Dumping: An Environmental Justice Framework.  _Qualitative Sociology_ 27, 4 (2004), 511–525.  https://doi.org/10.1023/B:QUAS.0000049245.55208.4b
  * Skogsmo (2020) Markus Skogsmo. 2020\.  _A Scalable Approach for Detecting Dumpsites using Automatic Target Recognition with Feature Selection and SVM through Satellite Imagery_.  Ph.D. Dissertation. Uppsala University, Disciplinary Domain of Science and Technology, Mathematics and Computer Science, Department of Information Technology, Division of Visual Information and Interaction.  http://urn.kb.se/resolve?urn=urn:nbn:se:uu:diva-418792
  * Vambol et al. (2019) Sergij Vambol, Viola Vambol, Muniyan Sundararajan, and Iqbal Ansari. 2019.  The nature and detection of unauthorized waste dump sites using remote sensing.  _Ecological Questions_ 30, 3 (May 2019), 1\.  https://doi.org/10.12775/EQ.2019.018
  * Wilson and Velis (2015) David C Wilson and Costas A Velis. 2015.  Waste management – still a global challenge in the 21st century: An evidence-based call for action.  _Waste Management & Research: The Journal for a Sustainable Circular Economy_ 33, 12 (Dec. 2015), 1049–1051.  https://doi.org/10.1177/0734242X15616055

##  7\. About the authors

Maria Roberta Devesa is PhD. in high energy physics. Currently, works in
Computer Vision and Machine Learning since more than 4 years and has more than
7 years of experience in data analysis in different fields.

Antonio Vazquez Brust holds a MS degree in Urban Informatics. His area of
interest is the application of computational analysis to study urban
phenomena. He collaborates with NGOs and international agencies to better
understand and govern territorial processes.