# A framework to measure the robustness of programs in the unpredictable
environment

Valentina Castiglioni  Reykjavik University, Iceland<EMAIL_ADDRESS><EMAIL_ADDRESS>,  Michele Loreti  University of Camerino, Italy
<EMAIL_ADDRESS>and  Simone Tini  University of Insubria, Italy
<EMAIL_ADDRESS>
###### Abstract.

Due to the diffusion of IoT, modern software systems are often thought to
control and coordinate _smart devices_ in order to manage assets and
resources, and to guarantee efficient behaviours. For this class of systems,
which interact extensively with humans and with their environment, it is thus
crucial to guarantee their “correct” behaviour in order to avoid unexpected
and possibly dangerous situations. In this paper we will present a framework
that allows us to _measure the robustness_ of systems. This is the ability of
a program to _tolerate_ changes in the environmental conditions and preserving
the original behaviour. In the proposed framework, the interaction of a
program with its environment is represented as a sequence of random variables
describing how both evolve in time. For this reason, the considered measures
will be defined among probability distributions of observed data. The proposed
framework will be then used to define the notions of _adaptability_ and
_reliability_. The former indicates the ability of a program to absorb
perturbation on environmental conditions after a given amount of time. The
latter expresses the ability of a program to maintain its intended behaviour
(up-to some reasonable tolerance) despite the presence of perturbations in the
environment. Moreover, an _algorithm_ , based on statistical inference, it
proposed to evaluate the proposed metric and the aforementioned properties.
Throughout the paper, two case studies are used to the describe and evaluate
the proposed approach.

###### Key words and phrases:

Adaptability, Reliability, Robustness, Hemimetrics, Unknown environment,
Probabilistic behaviour

##  1\. Introduction

With the advent of _IoT era_ [Kop11] we witnessed to the diffusion of a
plethora of _smart interconnected devices_ that permit _sensing data_ from the
_environment_ and to _operate_ on it. This revolution stimulated the design of
increasingly pervasive systems, often named _Cyber-Physical Systems_ (CPSs)
[RLSS10], that are thought to control and coordinate these devices in order to
manage assets and resources, and to guarantee efficient behaviours.

For this class of systems, which interact extensively with humans and with
their environment, it is thus crucial to guarantee their “correct” behaviour
in order to avoid unexpected and possibly dangerous situations. We can observe
that to evaluate this behaviour of a system, specific assumptions on the
environmental conditions are needed since the beginning of the design phase.
These assumptions can be expressed via formal tools that are able to describe
how the environment evolves and, if needed, how it reacts to the actions
performed by the program on it. For instance, while programming a drone that
is autonomously setting its trajectory to avoid obstacles, we have to consider
the possible wind scenarios describing the conditions under which our device
will operate. However, these scenarios cannot be fully deterministic and
should measure the likelihood of events that could constitute an hazard for
our device, like for instance a gust of wind. For this reason, _the evolution
of the environment_ can be described in terms of a _stochastic process_ that
gives a probabilistic measure of the _observed dynamics_ of the environment.
Similarly, the program controlling our system should incorporate _random
behaviours_ to tame and handle the uncertainty in the environmental
conditions. As a consequence, the interaction between our program and the
environment will be represented in terms of a _stochastic process_ whose
states describe how the environmental conditions and the state of our program
evolve along a computation.

In the literature, we can find a wealth of proposals of stochastic and
probabilistic models, as, e.g., Probabilistic Automata [Seg95], Stochastic
Process Algebras [HHH+94, CH08, BNL13], Labelled Markov Chains and Stochastic
Hybrid Systems [CLe07, HLS00], _Markov Decision Processes_ [Put05], and _ad
hoc_ solutions for specific application contexts, as, e.g., establishing
safety guarantees for drones flying under particular circumstances [VNGV16,
HJS+14]. Yet, in these studies, either the environment is not explicitly taken
into account or it is modelled only deterministically. In addition to that,
due to the variety of applications and heterogeneity of systems, no general
formal framework to deal with these challenges has been proposed so far. The
lack of concise abstractions and of an automatic support makes the analysis
and verification of the considered systems difficult, laborious, and error
prone.

We can observe that one of the main challenges in the analysis of CPSs is to
quantify how _variations in the environmental conditions_ can affect the
_expected behaviour_ of our system.

In this paper we will present a framework that allows us to _measure the
robustness_ of systems. This is the ability of a program to _tolerate_ changes
in the environmental conditions and preserving the original behaviour. Since
the interaction of the program with the enclosing environment is represented
as a sequence of random variables, the considered measures will be defined
among probability distributions of observed data. The proposed framework will
be then used to define the notions of _adaptability_ and _reliability_. The
former indicates the ability of a program to absorb perturbation on
environmental conditions after a given amount of time. The latter expresses
the ability of a program to maintain its intended behaviour (up-to some
reasonable tolerance) despite the presence of perturbations in the
environment. As an example, if we consider the above mentioned drone scenario,
a perturbation can be caused by a gust of wind that moves the drone out of its
trajectory. The controlling program is _adaptable_ if it can return on the
initial trajectory within a suitable amount of time while it is _reliable_ if
the drone is able to detect the presence of a gust of wind and can oppose to
it, being only slightly moved from its initial trajectory.

### Detailed outline of paper contributions

In this paper we first propose a semantic framework that permits modelling the
behaviour of programs operating in an evolving environment. In the proposed
approach, the interaction among these two elements is represented in a purely
_data-driven_ fashion. In fact, while the environmental conditions are
(partially) available to the program as a set of data, allowing it to adjust
its behaviour to the current situation, the program is also able to use data
to (partially) control the environment and fulfil its tasks. The program-
environment interplay is modelled in terms of the changes they induce on a set
of application-relevant data, henceforth called _data space_. This approach
will allow for a significant simplification in modelling the behaviour of the
program, which can be _isolated_ from that of the environment. Moreover, as
common to favour _computational tractability_ [ADB11, AKLP10], a _discrete
time_ approach is adopted.

In our model, a _system_ consists in _three distinct components_ :

  1. (1)

a _process_ $P$ describing the behaviour of the program,

  2. (2)

a _data state_ $\mathbf{d}$ describing the current state of the data space,
and

  3. (3)

an _environment evolution_ $\mathcal{E}$ describing the effect of the
environment on $\mathbf{d}$.

As we focus on the interaction with the environment, we abstract from the
internal computation of the program and model only its activity on
$\mathbf{d}$. At each step, a process can _read/update_ values in $\mathbf{d}$
and $\mathcal{E}$ applies on the resulting data state, providing a new data
state at the next step. To deal with the uncertainties, we introduce
_probability_ at two levels:

  1. (1)

we use the _discrete_ _generative probabilistic model_ [vGSS95] to define
processes, and

  2. (2)

$\mathcal{E}$ induces a _continuous_ _distribution_ over data states.

The behaviour of the system is then entirely expressed by its _evolution
sequence_ , i.e., the sequence of distributions over data states obtained at
each step. Given the novelties of our model, as a side contribution we show
that this behaviour defines a _Markov process_.

The second contribution is the definition of a _metric semantics_ in terms of
a (time-dependent) distance on the evolution sequences of systems, which we
call the _evolution metric_. The _evolution metric_ will allow us to:

  1. (1)

verify how well a program is fulfilling its tasks by comparing it with its
specification,

  2. (2)

compare the activity of different programs in the same environment,

  3. (3)

compare the behaviour of one program with respect to different environments
and changes in the initial conditions.

The evolution metric will consist of two components: a _metric on data states_
and the _Wasserstein metric_ [Vas69]. The former is defined in terms of a
(time-dependent) _penalty function_ allowing us to compare two data states
only on the base of the objectives of the program. The latter lifts the metric
on data states to a metric on distributions on data states. We then obtain a
metric on evolution sequences by considering the maximal of the Wasserstein
distances over time. The proposed metric is then used to formalise the notions
of _robustness_ , _adaptability_ , and _reliability_.

A _randomised algorithm_ that permits estimating the evolution sequences of
systems and thus for the evaluation of the evolution metric is then
introduced. Following [TK10], the Wasserstein metric is evaluated in time
$O(N\log N)$, where $N$ is the (maximum) number of samples. We already adopted
this approach in [CLT20a] in the context of finite-states self-organising
collective systems, without any notion of environment or data space. Moreover,
we will also show how the proposed algorithm can be used to estimate the
_robustness_ , _adaptability_ , and _reliability_ , of a system.

### Organisation of contents

After reviewing some background mathematical notions in Section 2, we devote
Section 3 to the presentation of our _model_. In Section 4 we discuss the
_evolution metric_ and all the ingredients necessary to define it. Then, in
Section 5 we provide the _algorithm_ , based on statistical inference, for the
evaluation of the metric. The notions of _adaptability_ and _reliability_ of
programs are formally introduced in Section 6 together with an example of
application of our algorithm to their evaluation. In these sections a running
scenario is used to clarify the proposed approach. This consists of a
stochastic variant of the three-tanks laboratory experiment described in
[RKO+97]. In Section 7 our methodology is used to analyse the _engine system_
proposed in [LMMT21]. This consists of two supervised self-coordinating
refrigerated engine systems that are subject to cyber-physical attacks aiming
to inflict _overstress of equipment_ [GGI+15]. Finally, Section 8 concludes
the paper with a discussion on related work and future research directions.

##  2\. Background

In this section we present the mathematical background on which we build our
contribution.

### Measurable spaces and measurable functions

A _$\sigma$ -algebra_ over a set $\Omega$ is a family $\Sigma$ of subsets of
$\Omega$ s.t.:

$\Omega\in\Sigma$,

$\Sigma$ is closed under complementation; and

$\Sigma$ is closed under countable union. The pair $(\Omega,\Sigma)$ is called
a _measurable space_ and the sets in $\Sigma$ are called _measurable sets_ ,
ranged over by $\mathbb{A},\mathbb{B},\dots$.

For an arbitrary family $\Phi$ of subsets of $\Omega$, the $\sigma$-algebra
_generated_ by $\Phi$ is the smallest $\sigma$-algebra over $\Omega$
containing $\Phi$. In particular, we recall that given a topology $T$ over
$\Omega$, the _Borel $\sigma$-algebra_ over $\Omega$, denoted
${\mathcal{B}}(\Omega)$, is the $\sigma$-algebra generated by the open sets in
$T$. For instance, given $n\in\mathbb{N}^{+}$, ${\mathcal{B}}(\mathbb{R}^{n})$
is the $\sigma$-algebra generated by the open intervals in $\mathbb{R}^{n}$.

Given two measurable spaces $(\Omega_{i},\Sigma_{i})$, $i=1,2$, the _product
$\sigma$-algebra_ $\Sigma_{1}\otimes\Sigma_{2}$ is the $\sigma$-algebra on
$\Omega_{1}\times\Omega_{2}$ generated by the sets
$\\{\mathbb{A}_{1}\times\mathbb{A}_{2}\mid\mathbb{A}_{i}\in\Sigma_{i}\\}$.

Given measurable spaces $(\Omega_{1},\Sigma_{1}),(\Omega_{2},\Sigma_{2})$, a
function $f\colon\Omega_{1}\to\Omega_{2}$ is said to be
$\Sigma_{1}$-_measurable_ if $f^{-1}\\!(\mathbb{A}_{2})\\!\in\\!\Sigma_{1}$
for all $\mathbb{A}_{2}\\!\in\\!\Sigma_{2}$, with
$f^{-1}(\mathbb{A}_{2})\\!=\\!\\{\omega\\!\in\\!\Omega_{1}\mid
f(\omega)\in\mathbb{A}_{2}\\}$.

### Probability spaces and random variables

A _probability measure_ on a measurable space $(\Omega,\Sigma)$ is a function
$\mu\colon\Sigma\to[0,1]$ s.t. $\mu(\Omega)=1$, $\mu(\mathbb{A})\geq 0$ for
all $\mathbb{A}\in\Sigma$ and $\mu(\bigcup_{i\in I}\mathbb{A}_{i})=\sum_{i\in
I}\mu(\mathbb{A}_{i})$ for every countable family of pairwise disjoint
measurable sets $\\{\mathbb{A}_{i}\\}_{i\in I}\subseteq\Sigma$. Then
$(\Omega,\Sigma,\mu)$ is called a _probability space_. We let
$\Pi{(\Omega,\Sigma)}$ denote the set of all probability measures over
$(\Omega,\Sigma)$.

For $\omega\in\Omega$, the _Dirac measure_ $\mathbf{\delta}_{\omega}$ is
defined by $\mathbf{\delta}_{\omega}(\mathbb{A})=1$, if $\omega\in\mathbb{A}$,
and $\mathbf{\delta}_{\omega}(\mathbb{A})=0$, otherwise, for all
$\mathbb{A}\in\Sigma$. Given a countable set of reals $(p_{i})_{i\in I}$ with
$p_{i}\geq 0$ and $\sum_{i\in I}p_{i}=1$, the _convex combination_ of the
probability measures $\\{\mu_{i}\\}_{i\in I}\subseteq\Pi{(\Omega,\Sigma)}$ is
the probability measure $\sum_{i\in I}p_{i}\cdot\mu_{i}$ in
$\Pi{(\Omega,\Sigma)}$ defined by $(\sum_{i\in
I}p_{i}\cdot\mu_{i})(\mathbb{A})=\sum_{i\in I}p_{i}\mu_{i}(\mathbb{A})$, for
all $\mathbb{A}\in\Sigma$. A probability measure $\mu\in\Pi{(\Omega,\Sigma)}$
is called _discrete_ if $\mu=\sum_{i\in
I}p_{i}\cdot\mathbf{\delta}_{\omega_{i}}$, with $\omega_{i}\in\Omega$, for
some countable set of indexes $I$. The _support_ of a discrete probability
measure $\mu$ is defined as
$\mathsf{supp}(\mu)=\\{\omega\in\Omega\mid\mu(\omega)>0\\}$.

Assume a probability space $(\Omega,\Sigma,\mu)$ and a measurable space
$(\Omega^{\prime},\Sigma^{\prime})$. A function
$X\colon\Omega\to\Omega^{\prime}$ is called a _random variable_ if it is
$\Sigma$-measurable. The _distribution measure_ , or _cumulative distribution
function_ (cdf), of $X$ is the probability measure $\mu_{X}$ on
$(\Omega^{\prime},\Sigma^{\prime})$ defined by
$\mu_{X}(\mathbb{A})=\mu(X^{-1}(\mathbb{A}))$ for all
$\mathbb{A}\in\Sigma^{\prime}$. Given random variables $X_{i}$ from
$(\Omega_{i},\Sigma_{i},\mu_{i})$ to
$(\Omega^{\prime}_{i},\Sigma^{\prime}_{i})$, $i=1,\ldots,n$, the collection
$\textbf{X}=[X_{1},\ldots,X_{n}]$ is called a _random vector_. The cdf of a
random vector X is given by the _joint_ distribution of the random variables
in it.

With a slight abuse of notation, in the examples and explanations throughout
the paper, we will sometimes use directly the cdf of a random variable rather
than formally introducing the probability measure defined on the domain space.
In fact, since we will consider Borel sets over n, for some $n\geq 1$, it is
usually easier to define the cdfs than the probability measures on them.
Similarly, when the cdf is absolutely continuous with respect to the Lebesgue
measure, then we shall reason directly on the _probability density function_
(pdf) of the random variable, namely the Radon-Nikodym derivative of the cdf
with respect to the Lebesgue measure.

Consequently, we shall also henceforth use the more suggestive term
_distribution_ in place of the terms probability measure, cdf and pdf.

### The Wasserstein hemimetric

A _metric_ on a set $\Omega$ is a function
$m\colon\Omega\times\Omega\to{}^{\geq 0}$ s.t. $m(\omega_{1},\omega_{2})=0$
iff $\omega_{1}=\omega_{2}$,
$m(\omega_{1},\omega_{2})=m(\omega_{2},\omega_{1})$, and
$m(\omega_{1},\omega_{2})\leq
m(\omega_{1},\omega_{3})+m(\omega_{3},\omega_{2})$, for all
$\omega_{1},\omega_{2},\omega_{3}\in\Omega$. We obtain a _pseudometric_ by
relaxing the first property to $m(\omega_{1},\omega_{2})=0$ if
$\omega_{1}=\omega_{2}$. Then, an _hemimetric_ is a pseudometric that is not
necessarily symmetric. A (pseudo-, hemi-)metric $m$ is $l$-_bounded_ if
$m(\omega_{1},\omega_{2})\leq l$ for all $\omega_{1},\omega_{2}\in\Omega$. For
a (pseudo-, hemi-)metric on $\Omega$, the pair $(\Omega,m)$ is a (_pseudo_ -,
_hemi-_)_metric space_.

Given a (pseudo-, hemi-)metric space $(\Omega,m)$, the (pseudo-, hemi-)metric
$m$ induces a natural topology over $\Omega$, namely the topology generated by
the open $\varepsilon$-balls, for $\varepsilon>0$,
$B_{m}(\omega,\varepsilon)=\\{\omega^{\prime}\in\Omega\mid
m(\omega,\omega^{\prime})<\varepsilon\\}$. We can then naturally obtain the
Borel $\sigma$-algebra over $\Omega$ from this topology.

In this paper we are interested in defining a _hemimetric on distributions_.
To this end we will make use of the Wasserstein lifting [Vas69] whose
definition is based on the following notions and results. Given a set $\Omega$
and a topology $T$ on $\Omega$, the topological space $(\Omega,T)$ is said to
be _completely metrisable_ if there exists at least one metric $m$ on $\Omega$
such that $(\Omega,m)$ is a complete metric space and $m$ induces the topology
$T$. A _Polish space_ is a separable completely metrisable topological space.
In particular, we recall that: [(i)]

is a Polish space; and

every closed subset of a Polish space is in turn a Polish space. Moreover, for
any $n\in N$, if $\Omega_{1},\dots,\Omega_{n}$ are Polish spaces, then the
Borel $\sigma$-algebra on their product coincides with the product
$\sigma$-algebra generated by their Borel $\sigma$-algebras, namely

${\mathcal{B}}(\bigtimes_{i=1}^{n}\Omega_{i})=\bigotimes_{i=1}^{n}{\mathcal{B}}(\Omega_{i}).$

(This is proved, e.g., in [Bog07] as Lemma 6.4.2 whose hypothesis are
satisfied by Polish spaces since they are second countable.) These properties
of Polish spaces are interesting for us since they guarantee that all the
distributions we consider in this paper are Radon measures and, thus, the
Wasserstein lifting is well-defined on them. For this reason, we also directly
present the Wasserstein hemimetric by considering only distributions on Borel
sets.

[Wasserstein hemimetric] Consider a Polish space $\Omega$ and let $m$ be a
hemimetric on $\Omega$. For any two distributions $\mu$ and $\nu$ on
$(\Omega,{\mathcal{B}}(\Omega))$, the _Wasserstein lifting_ of $m$ to a
distance between $\mu$ and $\nu$ is defined by

$\operatorname{\mathbf{W}}(m)(\mu,\nu)=\inf_{\mathfrak{w}\in\mathfrak{W}(\mu,\nu)}\int_{\Omega\times\Omega}m(\omega,\omega^{\prime})\mathrm{d}\mathfrak{w}(\omega,\omega^{\prime})$

where $\mathfrak{W}(\mu,\nu)$ is the set of the _couplings of $\mu$ and
$\nu$_, namely the set of joint distributions $\mathfrak{w}$ over the product
space $(\Omega\times\Omega,{\mathcal{B}}(\Omega\times\Omega))$ having $\mu$
and $\nu$ as left and right marginal, respectively, namely
$\mathfrak{w}(\mathbb{A}\times\Omega)=\mu(\mathbb{A})$ and
$\mathfrak{w}(\Omega\times\mathbb{A})=\nu(\mathbb{A})$, for all
$\mathbb{A}\in{\mathcal{B}}(\Omega)$.

Despite the original version of the Wasserstein distance being defined on a
metric on $\Omega$, the Wasserstein hemimetric given above is well-defined. We
refer the interested reader to [FR18] and the references therein for a formal
proof of this fact. In particular, the Wasserstein hemimetric is given in
[FR18] as Definition 7 (considering the compound risk excess metric defined in
Equation (31) of that paper), and Proposition 4 in [FR18] guarantees that it
is indeed a well-defined hemimetric on $\Pi(\Omega,{\mathcal{B}}(\Omega))$.
Moreover, Proposition 6 in [FR18] guarantees that the same result holds for
the hemimetric $m(x,y)=\max\\{y-x,0\\}$ which, as we will see, plays an
important role in our work (cf. Definition 4 below).

######  Remark 1.

As elsewhere in the literature, for simplicity and brevity, we shall
henceforth use the term _metric_ in place of the term hemimetric.

##  3\. The model

We devote this section to introduce the three components of our systems,
namely the _data space_ $\mathcal{D}$, the _process_ $P$ describing the
behaviour of the program, and the _environment evolution_ $\mathcal{E}$
describing the effects of the environment. As already discussed in the
Introduction, our choice of modelling program and environment separately will
allow us to isolate programs from the environment and favour their analysis.

In order to help the reader grasp the details of our approach, we introduce a
stochastic variant of the three-tanks laboratory experiment, that embodies the
kind of program-environment interactions we are interested in. Several
variants (with different number of tanks) of the $n$-tanks experiment have
been widely used in control program education (see, among the others, [AL92,
RKO+97, Joh00, ALGG+06]). Moreover, some recently proposed cyber-physical
security testbeds, like _SWaT_ [MT16, AGA+17], can be considered as an
evolution of the tanks experiment. Here, we consider a variant of the three-
tanks laboratory experiment described in [RKO+97].

Please notice that this example is meant to serve as a toy example allowing us
to showcase the various features of our tool, without carrying out a heavy
mathematical formulation. We delay to Section 7 the presentation of a case-
study showing how our framework can be applied to the analysis of (more
complex) real-world systems.

In the upcoming examples we will slightly abuse of notation and use a variable
name $x$ to denote all: the variable $x$, the (possible) function describing
the evolution in time of the values assumed by $x$, and the (possible) random
variable describing the distribution of the values that can be assumed by $x$
at a given time. The role of the variable name $x$ will always be clear from
the context.

Tank 1Tank 2Tank
3$q_{12}$$q_{23}$$q_{3}$$l_{max}$$l_{min}$$l_{goal}$$l_{1}$$l_{2}$$l_{3}$pump$q_{1}$

Flow rate under

the control of

the software

$q_{2}$

Flow rate under

the control of

the environment

pump

Flow rate under

the control of

the software

Figure 1. Schema of the three-tanks scenario.

As schematised in Figure 1, there are three identical tanks connected by two
pipes. Water enters in the first and in the last tank by means of a pump and
an incoming pipe, respectively. The last tank is equipped with an outlet pump.
We assume that water flows through the incoming pipe with a rate $q_{2}$ that
is determined by the environment, whereas the flow rates through the two pumps
(respectively, $q_{1}$ and $q_{3}$) are under the control of the program
controlling the experiment. The rates $q_{1},q_{2},q_{3}$ can assume values in
the range $[0,q_{max}]$, for a given a maximal flow rate $q_{max}$. The task
of the program consists in guaranteeing that the levels of water in the three
tanks fulfil some given requirements.

The level of water in tank $i$ at time $\tau$ is denoted by $l_{i}(\tau)$, for
$i=1,2,3$, and is always in the range $[l_{min},l_{max}]$, for suitable
$l_{min}$ and $l_{max}$ giving, respectively, the minimum and maximum level of
water in the tanks.

The principal difference between our version of the three-tanks and its
classical formulation is the flow rate $q_{2}$. In fact, we assume that
$q_{2}$ cannot be controlled by the program, that can only indirectly observe
the effects it has on $l_{3}$ and react to those. In our setting, this is
equivalent to say that $q_{2}$ is under the control of the environment.
Moreover, to obtain a more complex scenario, we assume that the value of
$q_{2}$ (at a given time) can be described only probabilistically. This can be
interpreted either as the effect of uncertainties (like, e.g., measurement
errors of the sensors having to determine $q_{2}$), or as the consequence of
an unpredictable behaviour of the environment. Consequently, the exact value
of $q_{2}(\tau)$ is unknown and will be evident only at execution time. Still,
we can consider different scenarios that render the assumptions we might have
on the environment. For instance, we can assume that the flow rate of the
incoming pipe is normally distributed with mean $q_{med}\in[0,q_{max}]$ and
variance $\delta_{q}>0$:

$\begin{array}[]{rcl}q_{2}(\tau+\Delta\tau)&\sim&N(q_{med},\delta_{q})\enspace.\end{array}$
(1)

In a more elaborated scenario, we could assume that $q_{2}$ varies at each
step by a value $v$ that is normally distributed with mean $0$ and variance
$1$. In this case, we have:

$\begin{array}[]{rcl}v(\tau)&\sim&N(0,1)\\\
q_{2}(\tau+\Delta\tau)&=&\min\left\\{\max\left\\{0,q_{2}(\tau)+v(\tau)\right\\},q_{max}\right\\}\enspace.\end{array}$
(2)

Conversely, the two pumps are controlled by the program that, by reading the
values of $l_{1}(\tau)$ and $l_{3}(\tau)$, can select the value of
$q_{1}(\tau+\Delta\tau)$ and $q_{3}(\tau+\Delta\tau)$. In detail, the
equations describing the behaviour of the program are the following:

$q_{1}(\tau+\Delta\tau)\\!=\\!\left\\{\\!\\!\begin{array}[]{ll}\max\\{0,q_{1}(\tau)-q_{step}\\}&\\!\text{if
}l_{1}(\tau)>l_{goal}+\delta_{l}\\\
\min\\{q_{max},q_{1}(\tau)+q_{step}\\}&\\!\text{if
}l_{1}(\tau)<l_{goal}-\delta_{l}\\\
q_{1}(\tau)&\\!\text{otherwise;}\end{array}\right.$ (3)
$q_{3}(\tau+\Delta\tau)\\!=\\!\left\\{\\!\\!\begin{array}[]{ll}\min\\{q_{max},q_{3}(\tau)+q_{step}\\}&\\!\text{if
}l_{3}(\tau)>l_{goal}+\delta_{l}\\\
\max\\{0,q_{3}(\tau)-q_{step}\\}&\\!\text{if
}l_{3}(\tau)<l_{goal}-\delta_{l}\\\
q_{3}(\tau)&\\!\text{otherwise.}\end{array}\right.$ (4)

Above, $l_{goal}$ is the desired level of water in the tanks, while
$q_{step}\in[0,q_{max})$ is the variation of the flow rate that is
controllable by the pump in one time step $\Delta\tau$. The idea behind
Equation (3) is that when $l_{1}$ is greater than $l_{goal}+\delta_{l}$, the
flow rate of the pump is decreased by $q_{step}$. Similarly, when $l_{1}$ is
less than $l_{goal}-\delta_{l}$, that rate is increased by $q_{step}$. The
reasoning for Equation (4) is symmetrical. In both cases, we use the threshold
$\delta_{l}>0$ to avoid continuous contrasting updates.

From now on, we will say that we are in _scenario 1_ if $q_{2}$ is determined
by Equation (1), and that we are in _scenario 2_ in case $q_{2}$ is determined
by Equation (2). In both scenarios, $q_{1}$ is determined by Equation (3) and
$q_{3}$ by Equation (4).

The dynamics of $l_{i}(\tau)$ can be modelled via the following set of
stochastic difference equations, with sampling time interval $\Delta\tau$:

$\begin{array}[]{rcl}l_{1}(\tau+\Delta\tau)&=&l_{1}(\tau)+\Delta\tau\cdot(q_{1}(\tau)-q_{12}(\tau))\\\
l_{2}(\tau+\Delta\tau)&=&l_{2}(\tau)+\Delta\tau\cdot(q_{12}(\tau)-q_{23}(\tau))\\\
l_{3}(\tau+\Delta\tau)&=&l_{3}(\tau)+\Delta\tau\cdot(q_{2}(\tau)+q_{23}(\tau)-q_{3}(\tau))\end{array}$
(5)

where $q_{1},q_{2},q_{3}$ are as above (i.e., the flow rates through,
respectively, the pump connected to the first tank, the incoming pipe, and the
outlet pump), and $q_{ij}$ denotes the flow rate from tank $i$ to tank $j$.
The evaluation of $q_{ij}(\tau)$ depends on $l_{i}(\tau),l_{j}(\tau)$, and on
the physical dimensions of the tanks, and it is obtained as follows. Let $A$
be the cross sectional area of each tank, and $a$ be the cross sectional area
of the connecting and outlet pipes. The volume balance difference equations of
each tank are then the following:

$\begin{array}[]{rcl}\displaystyle\frac{A}{\Delta\tau}\cdot(l_{1}(\tau+\Delta\tau)-l_{1}(\tau))&=&q_{1}(\tau)-q_{12}(\tau)\\\\[8.5359pt]
\displaystyle\frac{A}{\Delta\tau}\cdot(l_{2}(\tau+\Delta\tau)-l_{2}(\tau))&=&q_{12}(\tau)-q_{23}(\tau)\\\\[8.5359pt]
\displaystyle\frac{A}{\Delta\tau}\cdot(l_{3}(\tau+\Delta\tau)-l_{3}(\tau))&=&q_{2}(\tau)+q_{23}(\tau)-q_{3}(\tau).\\\
\end{array}$ (6)

We then apply Torricelli’s law to the equations in (6) to obtain the flow
rates $q_{12}$ and $q_{23}$:

$\begin{split}&q_{12}(\tau)=\begin{cases}a_{12}(\tau)\cdot
a\cdot\sqrt{2g\big{(}l_{1}(\tau)-l_{2}(\tau)\big{)}}&\text{ if
}l_{1}(\tau)\geq l_{2}(\tau)\\\\[2.84544pt] -a_{12}(\tau)\cdot
a\cdot\sqrt{2g\big{(}l_{2}(\tau)-l_{1}(\tau)\big{)}}&\text{
otherwise;}\end{cases}\\\ &q_{23}(\tau)=\begin{cases}a_{23}(\tau)\cdot
a\cdot\sqrt{2g\big{(}l_{2}(\tau)-l_{3}(\tau)\big{)}}&\text{ if
}l_{2}(\tau)\geq l_{3}(\tau)\\\\[2.84544pt] -a_{23}(\tau)\cdot
a\cdot\sqrt{2g\big{(}l_{3}(\tau)-l_{2}(\tau)\big{)}}&\text{
otherwise;}\end{cases}\end{split}$ (7)

where $g$ is the gravitational constant, and $a_{12},a_{23}$ are the loss
coefficients of the respective pipes. These coefficients are represented as
time dependent functions since they depend on the geometry of the pipes and on
the water level in the tanks. In our experiments, we use the approximation
$g=9.81$, and the following values for the aforementioned coefficients: [(i)]

$a=0.5$,

$a_{12}=0.75$,

$a_{23}=0.75$.  $\LHD$

We now proceed to formally introduce the three components of our systems and
their interaction.

### Modelling the data space.

We define the data space by means of a _finite_ set of _variables_
$\mathrm{Var}$ representing: [i)]

_environmental conditions_ , such as pressure, temperature, humidity, etc.,

_values perceived by sensors_ , which depend on the value of environmental
conditions and are unavoidably affected by imprecision and approximations
introduced by sensors, and

_state of actuators_ , which are usually elements in a discrete domain, like
$\\{\mathit{on},\mathit{off}\\}$. For each $x\in\mathrm{Var}$ we assume a
measurable space $(\mathcal{D}_{x},{\mathcal{B}}_{x})$, with
$\mathcal{D}_{x}\subseteq\real$ the domain of $x$ and ${\mathcal{B}}_{x}$ the
Borel $\sigma$-algebra on $\mathcal{D}_{x}$. Without loosing generality, we
can assume that each domain $\mathcal{D}_{x}$ is either a _finite set_ or a
_compact_ subset of . In particular, this means that each $\mathcal{D}_{x}$ is
a Polish space.

As $\mathrm{Var}$ is a finite set, we can always assume it to be ordered,
namely $\mathrm{Var}=\\{x_{1},\dots,x_{n}\\}$ for a suitable $n\in\mathbb{N}$.

[Data space] We define the _data space_ over $\mathrm{Var}$, notation
$\mathcal{D}_{\mathrm{Var}}$, as the Cartesian product of the variables
domains, namely
$\mathcal{D}_{\mathrm{Var}}=\bigtimes_{i=1}^{n}\mathcal{D}_{x_{i}}$. Then, as
a $\sigma$-algebra on $\mathcal{D}_{\mathrm{Var}}$ we consider the the product
$\sigma$-algebra
${\mathcal{B}}_{\mathcal{D}_{\mathrm{Var}}}=\bigotimes_{i=1}^{n}{\mathcal{B}}_{x_{i}}$.

When no confusion arises, we will use $\mathcal{D}$ and
${\mathcal{B}}_{\mathcal{D}}$ in place of $\mathcal{D}_{\mathrm{Var}}$ and
${\mathcal{B}}_{\mathcal{D}_{\mathrm{Var}}}$, respectively.

The data space $\mathcal{D}_{3T}$ for the considered three-tanks system is
defined on the set of variables
$\mathrm{Var}=\\{q_{1},q_{2},q_{3},l_{1},l_{2},l_{3}\\}$, with domains
$q_{i}\in\mathcal{D}_{q}=[0,q_{max}]$ and
$l_{i}\in\mathcal{D}_{l}=[l_{min},l_{max}]$, for $i=1,2,3$. We remark that
$\delta_{l}$ is not included in the data space. This is due to the fact that,
in our setting, $\delta_{l}$ is never modified by the environment (it is also
constant in time), and it is an inherent property of the program: different
programs might have different values of $\delta_{l}$, according to how capable
they are of controlling and modifying the flow rates, or of reading the levels
of water in the tanks.  $\LHD$

Elements in $\mathcal{D}$ are $n$-ples of the form $(v_{1},\dots,v_{n})$, with
$v_{i}\in\mathcal{D}_{x_{i}}$, which can be also identified by means of
functions $\mathbf{d}\colon\mathrm{Var}\to\real$ from variables to values,
with $\mathbf{d}(x)\in\mathcal{D}_{x}$ for all $x\in\mathrm{Var}$. Each
function $\mathbf{d}$ identifies a particular configuration of the data in the
data space, and it is thus called a _data state_.

[Data state] A _data state_ is a mapping
$\mathbf{d}\colon\mathrm{Var}\to\real$ from state variables to values, with
$\mathbf{d}(x)\in\mathcal{D}_{x}$ for all $x\in\mathrm{Var}$.

For simplicity, we shall write $\mathbf{d}\in\mathcal{D}$ in place of
$(\mathbf{d}(x_{1}),\dots,\mathbf{d}(x_{n}))\in\mathcal{D}$. Since program and
environment interact on the basis of the _current_ values of data, we have
that at each step there is a data state $\mathbf{d}$ that identifies the
_current state of the data space_ on which the next computation step is built.
Given a data state $\mathbf{d}$, we let $\mathbf{d}[x=v]$ denote the data
state $\mathbf{d}^{\prime}$ associating $v$ with $x$, and $\mathbf{d}(y)$ with
$y$ for any $y\neq x$. For $\overline{x}=(x_{i_{1}},\ldots,x_{i_{k}})$ and
$\overline{v}=(v_{i_{1}},\ldots,v_{i_{k}})$, we let
$\mathbf{d}[\overline{x}=\overline{v}]$ denote
$\mathbf{d}[x_{i_{1}}=v_{i_{1}}]\cdots[x_{i_{k}}=v_{i_{k}}]$.

### Modelling processes.

We introduce a simple process calculus allowing us to specify programs that
interact with a data state $\mathbf{d}$ in a given environment. We assume that
the action performed by a process at a given computation step is determined
probabilistically, according to the _generative_ probabilistic model [vGSS95].

[Syntax of processes] We let $\mathcal{P}$ be the set of _processes_ $P$
defined by:

$\begin{array}[]{rcl}P&::=&(\overline{e}\rightarrow\overline{x}).P^{\prime}\
\;\;|\;\;\mathrm{if}\;[e]\;P_{1}\;\mathrm{else}\;P_{2}\ \;\;|\;\;\sum_{i\in
I}p_{i}\cdot P_{i}\
\;\;|\;\;P_{1}\mathbin{\|}_{p}P_{2}\;\;|\;\;P_{1}~{}|~{}P_{2}\;\;|\;\;A\\\
e&::=&v\in\mathbb{R}\;\;|\;\;x\in\mathrm{Var}\;\;|\;\;op_{k}(e_{1},\ldots,e_{k})\end{array}$

where $p,p_{1},\ldots$ range over _probability weights_ in $[0,1]$, $I$ is
finite, $A$ ranges over _process variables_ , $op_{k}$ indicates a _measurable
operator_ $\mathbb{R}^{k}\rightarrow\mathbb{R}$, and $\overline{\cdot}$
denotes a finite sequence of elements. We assume that for each process
variable $A$ we have a single definition
$A\stackrel{{\scriptstyle\mathit{def}}}{{=}}P$. Moreover, we require that
$\sum_{i\in I}p_{i}=1$ for any process $\sum_{i\in I}p_{i}\cdot P_{i}$, and
that $\mathrm{var}(P_{1})\cap\mathrm{var}(P_{2})=\emptyset$ for any process
$P_{1}~{}|~{}P_{2}$, where $\mathrm{var}(P)=\\{x\in\mathrm{Var}\mid
x\in\overline{x}\text{ for any }(\overline{e}\rightarrow\overline{x})\text{
occurring in }P\\}$.

As usual, the expression $(\overline{e}\rightarrow\overline{x})$ in
$(\overline{e}\rightarrow\overline{x}).P^{\prime}$ is called a _prefix_. We
may use $\surd$ to denote the prefix $(\emptyset\rightarrow\emptyset)$. In a
single action a process can read and update a set of state variables. This is
done by process $(\overline{e}\rightarrow\overline{x}).P$ that first evaluates
the sequence of expressions $\overline{e}$ in the current data state and then
assigns the results to sequence of variables $\overline{x}$, which are stored
in the data state and evolve according to $\mathcal{E}$ before being available
at next step, when the process will behave as $P$. Process
$\mathrm{if}\;[e]\;P_{1}\;\mathrm{else}\;P_{2}$ behaves either as $P_{1}$, if
$e$ evaluates to $\top$, or as $P_{2}$, if $e$ evaluates to $\bot$. Then,
$\sum_{i\in I}p_{i}\cdot P_{i}$ is the _generative probabilistic choice_ : it
has probability $p_{i}$ to behave as $P_{i}$. We may write
$\sum_{i=1}^{n}p_{i}\cdot P_{i}$ as $p_{1}\cdot P_{1}+\ldots+p_{n}\cdot
P_{n}$. The _generative probabilistic interleaving_ construct
$P_{1}\mathbin{\|}_{p}P_{2}$ lets the two argument processes to interleave
their actions, where at each step the first process moves with probability $p$
and the second with probability $(1-p)$. We also provide a _synchronous
parallel composition_ construct $P_{1}~{}|~{}P_{2}$ that lets the two
arguments perform their actions in a synchronous fashion, provided $P_{1}$ and
$P_{2}$ do not modify the same variables. Finally, process variable $A$ allows
us to specify recursive behaviours by equations of the form
$A\stackrel{{\scriptstyle\mathit{def}}}{{=}}P$. To avoid Zeno behaviours we
assume that all occurrences of process variables appear _guarded_ by prefixing
constructs in $P$. We assume the standard notion of _free_ and _bound_ process
variable and we only consider _closed process terms_ , that is terms without
_free_ variables.

Actions performed by a process can be abstracted in terms of the _effects_
they have on the data state, i.e., via _substitutions_ of the form
$\theta=[x_{i_{1}}\leftarrow v_{i_{1}},\ldots,x_{i_{k}}\leftarrow v_{i_{k}}]$,
also denoted by $\overline{x}\leftarrow\overline{v}$ for
$\overline{x}=x_{i_{1}},\ldots,x_{i_{k}}$ and
$\overline{v}=v_{i_{1}},\ldots,v_{i_{k}}$. Since in Definition 3 process
operations on data and expressions are assumed to be measurable, we can model
the effects as ${\mathcal{B}}_{\mathcal{D}}$-measurable functions
$\theta\colon\mathcal{D}\to\mathcal{D}$ s.t.
$\theta(\mathbf{d}):=\mathbf{d}[\overline{x}=\overline{v}]$ whenever
$\theta=\overline{x}\leftarrow\overline{v}$. We denote by $\Theta$ the set of
effects. The behaviour of a process can then be defined by means of a function
$\mathsf{pstep}\colon\mathcal{P}\times\mathcal{D}\to\Pi(\Theta\times\mathcal{P})$
that given a process $P$ and a data state $\mathbf{d}$ yields a _discrete_
distribution over $\Theta\times\mathcal{P}$. The distributions in
$\Pi(\Theta\times\mathcal{P})$ are ranged over by $\pi,\pi^{\prime},\dots$.
Function $\mathsf{pstep}$ is formally defined in Table 1.

$\begin{array}[]{lcl}(\scalebox{0.9}{PR1})&&\mathsf{pstep}((\overline{e}\rightarrow\overline{x}).P^{\prime},\mathbf{d})=\delta_{(\overline{x}\leftarrow\llbracket\overline{e}\rrbracket_{\mathbf{d}},P^{\prime})}\\\\[2.84544pt]
(\scalebox{0.9}{PR2})&&\mathsf{pstep}(\mathrm{if}\;[e]\;P_{1}\;\mathrm{else}\;P_{2},\mathbf{d})=\begin{cases}\mathsf{pstep}(P_{1},\mathbf{d})&\text{
if }\llbracket e\rrbracket_{\mathbf{d}}=1\\\
\mathsf{pstep}(P_{2},\mathbf{d})&\text{ if }\llbracket
e\rrbracket_{\mathbf{d}}=0\end{cases}\\\\[2.84544pt]
(\scalebox{0.9}{PR3})&&\mathsf{pstep}(\sum_{i}p_{i}\cdot
P_{i},\mathbf{d})=\sum_{i}p_{i}\cdot\mathsf{pstep}(P_{i},\mathbf{d})\\\\[2.84544pt]
(\scalebox{0.9}{PR4})&&\mathsf{pstep}(P_{1}\mathbin{\|}_{p}P_{2},\mathbf{d})=p\cdot(\mathsf{pstep}(P_{1},\mathbf{d})\mathbin{\|}_{p}P_{2})+(1-p)\cdot(P_{1}\mathbin{\|}_{p}\mathsf{pstep}(P_{2},\mathbf{d}))\\\\[2.84544pt]
(\scalebox{0.9}{PR5})&&\mathsf{pstep}(P_{1}~{}|~{}P_{2},\mathbf{d})=\mathsf{pstep}(P_{1},\mathbf{d})~{}|~{}\mathsf{pstep}(P_{2},\mathbf{d})\\\\[2.84544pt]
(\scalebox{0.9}{PR6})&&\mathsf{pstep}(A,\mathbf{d})=\mathsf{pstep}(P,\mathbf{d})\qquad(\text{if
}A\stackrel{{\scriptstyle\mathit{def}}}{{=}}P)\\\ \end{array}$ Table 1.
Process Semantics

In detail, rule (PR1) expresses that in a single action a process can read and
update a set of variables: process $(\overline{e}\rightarrow\overline{x}).P$
first evaluates the sequence of expressions $\overline{e}$ in the data state
$\mathbf{d}$ and then stores the results to $\overline{x}$. These _new_ data
will evolve according to $\mathcal{E}$ before being available at next step,
when the process will behave as $P$. The evaluation of an expression $e$ in
$\mathbf{d}$, denoted by $\llbracket e\rrbracket_{\mathbf{d}}$, is obtained by
replacing each variable $x$ in $e$ with $\mathbf{d}(x)$ and by evaluating the
resulting ground expression. Then, process
$(\overline{e}\rightarrow\overline{x}).P$ evolves to the Dirac’s distribution
on the pair
$(\overline{x}\leftarrow\llbracket\overline{e}\rrbracket_{\mathbf{d}},P)$.
Rule (PR2) models the behaviour of
$\mathrm{if}\;[e]\;P_{1}\;\mathrm{else}\;P_{2}$. We assume that boolean
operations, like those in the $\mathrm{if}$-guard, evaluate to $1$ or $0$,
standing respectively for $\top$ and $\bot$. Hence,
$\mathrm{if}\;[e]\;P_{1}\;\mathrm{else}\;P_{2}$ behaves as $P_{1}$ if
$\llbracket e\rrbracket_{\mathbf{d}}=1$, and it behaves as $P_{2}$ otherwise.
Rule (PR3) follows the _generative_ approach to probabilistic choice: the
behaviour of process $P_{i}$ is selected with probability $p_{i}$. Rule (PR4)
states that process $P_{1}\mathbin{\|}_{p}P_{2}$ interleaves the moves by
$P_{1}$ and $P_{2}$, where at each step $P_{1}$ moves with probability $p$ and
$P_{2}$ with probability $(1-p)$. Given $\pi\in\Pi(\Theta\times\mathcal{P})$,
we let $\pi\mathbin{\|}_{p}P$ (resp. $P\mathbin{\|}_{p}\pi$) denote the
probability distribution $\pi^{\prime}$ over $(\Theta\times\mathcal{P})$ such
that: $\pi^{\prime}(\theta,P^{\prime})=\pi(\theta,P^{\prime\prime})$, whenever
$P^{\prime}=P^{\prime\prime}\mathbin{\|}_{p}P$ (resp.
$P^{\prime}=P\mathbin{\|}_{p}P^{\prime\prime}$), and $0$, otherwise. Rule
(PR5) gives us the semantics of the _synchronous_ parallel composition: for
$\pi_{1},\pi_{2}\in\Pi(\Theta\times\mathcal{P})$, we let
$\pi_{1}~{}|~{}\pi_{2}$ denote the probability distribution defined, for all
$\theta_{i}\in\Theta$ and $P_{i}\in\mathcal{P}$, $i=1,2$, by:

$(\pi_{1}~{}|~{}\pi_{2})(\theta_{1}\theta_{2},P_{1}~{}|~{}P_{2})=\pi_{1}(\theta_{1},P_{1})\cdot\pi_{2}(\theta_{2},P_{2})$

where the effect $\theta_{1}\theta_{2}$ is given by the concatenation of the
effects $\theta_{1}$ and $\theta_{2}$. Notice that the concatenation is well-
defined and it does not introduce contrasting effects on variables since we
are assuming that $\mathrm{var}(P_{1})\cap\mathrm{var}(P_{2})=\emptyset$.
Finally, rule (PR6) states that $A$ behaves like $P$ whenever
$A\stackrel{{\scriptstyle\mathit{def}}}{{=}}P$.

We can observe that $\mathsf{pstep}$ returns a discrete distribution over
$\Theta\times\mathcal{P}$.

######  Proposition 2 (Properties of process semantics).

Let $P\in\mathcal{P}$ and $\mathbf{d}\in\mathcal{D}$. Then
$\mathsf{pstep}(P,\mathbf{d})$ is a distribution with finite support, namely
the set
$\mathsf{supp}(\mathsf{pstep}(P,\mathbf{d}))=\\{(\theta,P^{\prime})\mid(\mathsf{pstep}(P,\mathbf{d}))(\theta,P^{\prime})>0\\}$
is finite and
$\sum_{(\theta,P^{\prime})\in\mathsf{supp}(\mathsf{pstep}(P,\mathbf{d}))}(\mathsf{pstep}(P,\mathbf{d}))(\theta,P^{\prime})=1$.

######  Proof 3.1.

The proof follows by a simple induction on the construction of
$\mathsf{pstep}(P,\mathbf{d})$.

We proceed to define the program $\mathsf{P_{Tanks}}$ responsible for
controlling the two pumps in order to achieve (and maintain) the desired level
of water in the three tanks. Intuitively, $\mathsf{P_{Tanks}}$ will consist of
a process $\mathsf{P_{in}}$ controlling the flow rate through the pump
attached to the first tank, and of a process $\mathsf{P_{out}}$ controlling
the flow rate through the outlet pump. The behaviour of $\mathsf{P_{in}}$ and
$\mathsf{P_{out}}$ can be obtained as a straightforward translation of,
respectively, Equation (3) and Equation (4) into our simple process calculus.
Formally:

$\begin{array}[]{lll}\mathsf{P_{in}}&\stackrel{{\scriptstyle\mathit{def}}}{{=}}&\mathrm{if}\;[l_{1}>l_{goal}+\delta_{l}]\;(\max\\{0,q_{1}-q_{step}\\}\rightarrow
q_{1}).\mathsf{P_{in}}\\\\[2.84544pt]
&&\mathrm{else}\quad\mathrm{if}\;[l_{1}<l_{goal}-\delta_{l}]\;(\min\\{q_{max},q_{1}+q_{step}\\}\rightarrow
q_{1}).\mathsf{P_{in}}\\\\[2.84544pt]
&&\phantom{\mathrm{else}}\quad\mathrm{else}\;\surd.\mathsf{P_{in}}\\\\[8.5359pt]
\mathsf{P_{out}}&\stackrel{{\scriptstyle\mathit{def}}}{{=}}&\mathrm{if}\;[l_{3}>l_{goal}+\delta_{l}]\;(\min\\{q_{max},q_{3}+q_{step}\\}\rightarrow
q_{3}).\mathsf{P_{out}}\\\\[2.84544pt]
&&\mathrm{else}\quad\mathrm{if}\;[l_{3}<l_{goal}-\delta_{l}]\;(\max\\{0,q_{3}-q_{step}\\}\rightarrow
q_{3}).\mathsf{P_{out}}\\\\[2.84544pt]
&&\phantom{\mathrm{else}}\quad\mathrm{else}\;\surd.\mathsf{P_{out}}\end{array}$

However, for $\mathsf{P_{Tanks}}$ to work properly, we need to allow
$\mathsf{P_{in}}$ and $\mathsf{P_{out}}$ to execute simultaneously. Since the
effects of read/update operators in $\mathsf{P_{in}}$ and $\mathsf{P_{out}}$
are applied to distinct variables (respectively, $q_{1}$ and $q_{3}$), we can
obtain the desired execution by defining $\mathsf{P_{Tanks}}$ as the
synchronous parallel composition of $\mathsf{P_{in}}$ and $\mathsf{P_{out}}$:

$\mathsf{P_{Tanks}}\stackrel{{\scriptstyle\mathit{def}}}{{=}}\mathsf{P_{in}}~{}|~{}\mathsf{P_{out}}$

$\LHD$

### Modelling the environment.

To model the action of the environment on data we use a mapping $\mathcal{E}$,
called _environment evolution_ , taking a data state to a distribution over
data states.

[Environment evolution] An _environment evolution_ is a function
$\mathcal{E}\colon\mathcal{D}\to\Pi(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$
s.t. for each $\mathbb{D}\in{\mathcal{B}}_{\mathcal{D}}$ the mapping
$\mathbf{d}\mapsto\mathcal{E}(\mathbf{d})(\mathbb{D})$ is
${\mathcal{B}}_{\mathcal{D}}$-measurable.

We notice that, due to the interaction with the program, in our model the
probability induced by $\mathcal{E}$ at the next time step _depends only_ on
the current state of the data space. It is then natural to assume that the
behaviour of the environment is modelled as a discrete time Markov process,
and to define function $\mathcal{E}$ we only need the initial distribution on
$\mathcal{D}$ (which in our setting will be the Dirac’s distribution on a data
state in $\mathcal{D}$) and the random vector $\mathbf{X}=[X_{1},\dots,X_{n}]$
in which each random variable $X_{i}$ models the changes induced by the
environment on the value of variable $x_{i}$.

To conclude the encoding of the three-tanks experiment from Example 3 into our
model, we need to define a suitable environment evolution. This can be derived
from Equations (5), (1) (or (2)), and (7) in the obvious way. For sake of
completeness, below we present functions $\mathcal{E}_{1}$ and
$\mathcal{E}_{2}$ modelling the evolution of the environment, respectively, in
scenario 1 (with Equation (1)) and scenario 2 (with Equation (2)). We remark
that:

  * •

The only stochastic variable in our setting is $q_{2}$, giving thus that the
distribution induced by $\mathcal{E}_{i}$ is determined by the distribution of
$q_{2}$, for $i=1,2$.

  * •

The values of $q_{1}$ and $q_{3}$ are modified by $\mathsf{P_{in}}$ and
$\mathsf{P_{out}}$, respectively. In particular, according to Definition 3,
this means that the environment does not affect the values of $q_{1}$ and
$q_{3}$ and, moreover, when we evaluate $\mathcal{E}_{i}(\mathbf{d})$, the
effects of the two process over $q_{1}$ and $q_{3}$ have already been taken
into account in $\mathbf{d}$, for $i=1,2$.

Define $f_{3T}\colon[0,q_{max}]\times\mathcal{D}_{3T}\to\mathcal{D}_{3T}$ as
the function such that, for any $x\in[0,q_{max}]$ and data state $\mathbf{d}$,
gives us the data state $f_{3T}(x,\mathbf{d})=\mathbf{d}^{\prime}$ obtained
from the solution of the following system of equations:

$\begin{cases}q_{12}=\begin{cases}a_{12}\cdot
a\cdot\sqrt{2g\big{(}\mathbf{d}(l_{1})-\mathbf{d}(l_{2})\big{)}}&\text{ if
}\mathbf{d}(l_{1})\geq\mathbf{d}(l_{2})\\\\[2.84544pt] -a_{12}\cdot
a\cdot\sqrt{2g\big{(}\mathbf{d}(l_{2})-\mathbf{d}(l_{1})\big{)}}&\text{
otherwise;}\end{cases}\\\ q_{23}=\begin{cases}a_{23}\cdot
a\cdot\sqrt{2g\big{(}\mathbf{d}(l_{2})-\mathbf{d}(l_{3})\big{)}}&\text{ if
}\mathbf{d}(l_{2})\geq\mathbf{d}(l_{3})\\\\[2.84544pt] -a_{23}\cdot
a\cdot\sqrt{2g\big{(}\mathbf{d}(l_{3})-\mathbf{d}(l_{2})\big{)}}&\text{
otherwise;}\end{cases}\\\
\mathbf{d}^{\prime}(l_{1})=\mathbf{d}(l_{1})+\Delta\tau\cdot(\mathbf{d}(q_{1})-q_{12})\\\
\mathbf{d}^{\prime}(l_{2})=\mathbf{d}(l_{2})+\Delta\tau\cdot(q_{12}-q_{23})\\\
\mathbf{d}^{\prime}(l_{3})=\mathbf{d}(l_{3})+\Delta\tau\cdot(x+q_{23}-\mathbf{d}(q_{3}))\\\
\mathbf{d}^{\prime}(q_{1})=\mathbf{d}(q_{1})\\\
\mathbf{d}^{\prime}(q_{2})=x\\\
\mathbf{d}^{\prime}(q_{3})=\mathbf{d}(q_{3})\end{cases}$

where $a,a_{12},a_{23},g,\Delta\tau$ are constants. Then, consider the random
variables $X\sim N(q_{med},\delta_{q})$ and $Y\sim N(0,1)$. From them, we
define the random variables $Q_{2}^{1}(\mathbf{d})=f_{3T}(X,\mathbf{d})$ and
$Q_{2}^{2}(\mathbf{d})=f_{3T}(\max\\{q_{max},\mathbf{d}(q_{2})+Y\\},\mathbf{d})$.
Then, for $i=1,2$ and for each data state $\mathbf{d}$, the distribution
$\mathcal{E}_{i}(\mathbf{d})$ is obtained as the distribution of $Q_{2}^{i}$.
$\LHD$

### Modelling system’s behaviour.

We use the notion of _configuration_ to model the state of the system at each
time step.

[Configuration] A _configuration_ is a triple $c=\langle
P,\mathbf{d}\rangle_{\mathcal{E}}$, where $P$ is a process, $\mathbf{d}$ is a
data state and $\mathcal{E}$ is an environment evolution. We denote by
$\mathcal{C}_{\mathcal{P},\mathcal{D},\mathcal{E}}$ the set of configurations
defined over $\mathcal{P},\mathcal{D}$ and $\mathcal{E}$. When no confusion
shall arise, we shall write $\mathcal{C}$ in place of
$\mathcal{C}_{\mathcal{P},\mathcal{D},\mathcal{E}}$.

Let $(\mathcal{P},\Sigma_{\mathcal{P}})$ be the measurable space of processes,
where $\Sigma_{\mathcal{P}}$ is the power set of $\mathcal{P}$, and
$(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$ be the measurable space of data
states. As $\mathcal{E}$ is fixed, we can identify $\mathcal{C}$ with
$\mathcal{P}\times\mathcal{D}$ and equip it with the product $\sigma$-algebra
$\Sigma_{\mathcal{C}}=\Sigma_{\mathcal{P}}\otimes{\mathcal{B}}_{\mathcal{D}}$:
$\Sigma_{\mathcal{C}}$ is generated by the sets
$\\{\langle\mathbb{P},\mathbb{D}\rangle_{\mathcal{E}}\mid\mathbb{P}\in\Sigma_{\mathcal{P}},\mathbb{D}\in{\mathcal{B}}_{\mathcal{D}}\\}$,
where $\langle\mathbb{P},\mathbb{D}\rangle_{\mathcal{E}}=\\{\langle
P,\mathbf{d}\rangle_{\mathcal{E}}\mid
P\in\mathbb{P},\mathbf{d}\in\mathbb{D}\\}$.

For $\mu_{\mathcal{P}}\in\Pi(\mathcal{P},\Sigma_{\mathcal{P}})$ and
$\mu_{\mathcal{D}}\in\Pi(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$ we let
$\mu=\langle\mu_{\mathcal{P}},\mu_{\mathcal{D}}\rangle_{\mathcal{E}}$ denote
the product distribution on $(\mathcal{C},\Sigma_{\mathcal{C}})$, i.e.,
$\mu(\langle\mathbb{P},\mathbb{D}\rangle_{\mathcal{E}})=\mu_{\mathcal{P}}(\mathbb{P})\cdot\mu_{\mathcal{D}}(\mathbb{D})$
for all $\mathbb{P}\in\Sigma_{\mathcal{P}}$ and
$\mathbb{D}\in{\mathcal{B}}_{\mathcal{D}}$. If $\mu_{\mathcal{P}}=\delta_{P}$
for some $P\in\mathcal{P}$, we shall denote
$\langle\delta_{P},\mu_{\mathcal{D}}\rangle_{\mathcal{E}}$ simply by $\langle
P,\mu_{\mathcal{D}}\rangle_{\mathcal{E}}$.

Our aim is to express the behaviour of a system in terms of the changes on
data. We start with the _one-step_ behaviour of a configuration, in which we
combine the effects on the data state induced by the activity of the process
(given by $\mathsf{pstep}$) and the subsequent action by the environment.
Formally, we define a function $\mathsf{cstep}$ that, given a configuration
$c$, yields a distribution on $(\mathcal{C},\Sigma_{\mathcal{C}})$ (Definition
3 below). Then, we use $\mathsf{cstep}$ to define the _multi-step_ behaviour
of configuration $c$ as a sequence
$\mathcal{S}_{c,0}^{\mathcal{C}},\mathcal{S}_{c,1}^{\mathcal{C}},\ldots$ of
distributions on $(\mathcal{C},\Sigma_{\mathcal{C}})$. To this end, we show
that $\mathsf{cstep}$ is a Markov kernel (Proposition 4 below). Finally, to
abstract from processes and focus only on data, from the sequence
$\mathcal{S}_{c,0}^{\mathcal{C}},\mathcal{S}_{c,1}^{\mathcal{C}},\ldots$, we
obtain a sequence of distributions
$\mathcal{S}_{c,0}^{\mathcal{D}},\mathcal{S}_{c,1}^{\mathcal{D}},\ldots$ on
$(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$ called the _evolution sequence_ of
the system (Definition 3.4 below).

[One-step semantics of configurations] Function
$\mathsf{cstep}\colon\mathcal{C}\to\Pi(\mathcal{C},\Sigma_{\mathcal{C}})$ is
defined for all configurations $\langle
P,\mathbf{d}\rangle_{\mathcal{E}}\in\mathcal{C}$ by

$\mathsf{cstep}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}})=\sum_{(\theta,P^{\prime})\in\mathsf{supp}(\mathsf{pstep}(P,\mathbf{d}))}\mathsf{pstep}(P,\mathbf{d})(\theta,P^{\prime})\cdot\langle
P^{\prime},\mathcal{E}(\theta(\mathbf{d}))\rangle_{\mathcal{E}}.$ (8)

The next result follows by
$\mathcal{E}(\theta(\mathbf{d}))\in\Pi(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$
and Proposition 2.

######  Proposition 3 (Properties of one-step configuration semantics).

Assume any configuration $\langle
P,\mathbf{d}\rangle_{\mathcal{E}}\in\mathcal{C}$. Then $\mathsf{cstep}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}})$ is a distribution on
$(\mathcal{C},\Sigma_{\mathcal{C}})$.

######  Proof 3.2.

Being $\mathbf{d}$ a data state in $\mathcal{D}$ and $\theta$ an effect,
$\theta(\mathbf{d})$ is a data state. Then, $\mathcal{E}(\theta(\mathbf{d}))$
is a distribution on $(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$, thus
implying that $\langle
P^{\prime},\mathcal{E}(\theta(\mathbf{d}))\rangle_{\mathcal{E}}$ is a
distribution on $(\mathcal{C},\Sigma_{\mathcal{C}})$ (see Notation 3).
Finally, by Proposition 2 we get that $\mathsf{cstep}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}})$ is a convex combination of distributions
on $(\mathcal{C},\Sigma_{\mathcal{C}})$ whose weights sum up to 1, which gives
the thesis.

Since $\mathsf{cstep}(c)\in\Pi(\mathcal{C},\Sigma_{\mathcal{C}})$ for each
$c\in\mathcal{C}$, we can rewrite
$\mathsf{cstep}\colon\mathcal{C}\times\Sigma_{\mathcal{C}}\to[0,1]$, so that
for each configuration $c\in\mathcal{C}$ and measurable set
$\mathbb{C}\in\Sigma_{\mathcal{C}}$, $\mathsf{cstep}(c)(\mathbb{C})$ denotes
the probability of reaching in one step a configuration in $\mathbb{C}$
starting from $c$. We can prove that $\mathsf{cstep}$ is the Markov kernel of
the Markov process modelling our system. This follows by Proposition 3 and by
proving that for each $\mathbb{C}\in\Sigma_{\mathcal{C}}$, the mapping
$c\mapsto\mathsf{cstep}(c)(\mathbb{C})$ is $\Sigma_{\mathcal{C}}$-measurable.

######  Proposition 4.

The function $\mathsf{cstep}$ is a Markov kernel.

######  Proof 3.3.

We need to show that $\mathsf{cstep}$ satisfies the two properties of Markov
kernels, namely

  1. (1)

For each configuration $c\in\mathcal{C}$, the mapping
$\mathbb{C}\mapsto\mathsf{cstep}(c)(\mathbb{C})$ is a distribution on
$(\mathcal{C},\Sigma_{\mathcal{C}})$.

  2. (2)

For each measurable set $\mathbb{C}\in\Sigma_{\mathcal{C}}$, the mapping
$c\mapsto\mathsf{cstep}(c)(\mathbb{C})$ is $\Sigma_{\mathcal{C}}$-measurable.

Item 1 follows directly by Proposition 3.

Let us focus on item 2. By Definition 3, for each $\langle
P,\mathbf{d}\rangle_{\mathcal{E}}\in\mathcal{C}$ we have that

$\mathsf{cstep}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}})(\mathbb{C})=\sum_{(\theta,P^{\prime})\in\mathsf{supp}(\mathsf{pstep}(P,\mathbf{d}))}\mathsf{pstep}(P,\mathbf{d})(\theta,P^{\prime})\cdot\langle
P^{\prime},\mathcal{E}(\theta(\mathbf{d}))\rangle_{\mathcal{E}}(\mathbb{C})\enspace.$

As, by definition, each $\theta\in\Theta$ and $\mathcal{E}$ are
${\mathcal{B}}_{\mathcal{D}}$-measurable functions, we can infer that also
their composition $\mathcal{E}(\theta(\cdot))$ is a
${\mathcal{B}}_{\mathcal{D}}$-measurable function. Since, moreover,
$\Sigma_{\mathcal{C}}$ is the (smallest) $\sigma$-algebra generated by
$\Sigma_{\mathcal{P}}\times{\mathcal{B}}_{\mathcal{D}}$ and every subset of
$\mathcal{P}$ is a measurable set in $\Sigma_{\mathcal{P}}$, we can also infer
that $\langle(\cdot),\mathcal{E}(\theta(\cdot))\rangle_{\mathcal{E}}$ is a
$\Sigma_{\mathcal{C}}$-measurable function. Finally, we recall that by
Proposition 2 we have that $\mathsf{supp}(\mathsf{pstep}(P,\mathbf{d}))$ is
finite. Therefore, we can conclude that $\mathsf{cstep}$ is a
$\Sigma_{\mathcal{C}}$-measurable function as linear combination of a finite
collection of $\Sigma_{\mathcal{C}}$-measurable functions (see, e.g., [Roy88,
Chapter 3.5, Proposition 19]).

Hence, the multi-step behaviour of configuration $c$ can be defined as a time
homogeneous Markov process having $\mathsf{cstep}$ as Markov kernel and
$\delta_{c}$ as initial distribution.

[Multi-step semantics of configurations] Let $c\in\mathcal{C}$ be a
configuration. The _multi-step behaviour_ of $c$ is the sequence of
distributions
$\mathcal{S}_{c,0}^{\mathcal{C}},\mathcal{S}_{c,1}^{\mathcal{C}},\ldots$ on
$(\mathcal{C},\Sigma_{\mathcal{C}})$ defined inductively as follows:

$\displaystyle\mathcal{S}_{c,0}^{\mathcal{C}}(\mathbb{C})=$
$\displaystyle\;\delta_{c}(\mathbb{C})\text{, for all
}\mathbb{C}\in\Sigma_{\mathcal{C}}$
$\displaystyle\mathcal{S}_{c,i+1}^{\mathcal{C}}(\mathbb{C})=$
$\displaystyle\int_{\mathcal{C}}\mathsf{cstep}(b)(\mathbb{C})\;\mathrm{d}(\mathcal{S}_{c,i}^{\mathcal{C}}(b))\text{,
for all }\mathbb{C}\in\Sigma_{\mathcal{C}}.$

We can prove that
$\mathcal{S}_{c,0}^{\mathcal{C}},\mathcal{S}_{c,1}^{\mathcal{C}},\ldots$ are
well defined, namely they are distributions on
$(\mathcal{C},\Sigma_{\mathcal{C}})$.

######  Proposition 5 (Properties of configuration multi-step semantics).

For any $c\in\mathcal{C}$, all
$\mathcal{S}_{c,0}^{\mathcal{C}},\mathcal{S}_{c,1}^{\mathcal{C}},\ldots$ are
distributions on $(\mathcal{C},\Sigma_{\mathcal{C}})$.

######  Proof 3.4.

The proof follows by induction. The base case for
$\mathcal{S}_{c,0}^{\mathcal{C}}$ is immediate. The inductive step follows by
Proposition 4.

As the program-environment interplay can be observed only in the changes they
induce on the data states, we define the _evolution sequence_ of a
configuration as the sequence of distributions over data states that are
reached by it, step-by-step.

[Evolution sequence] The _evolution sequence_ of a configuration $c=\langle
P,\mathbf{d}\rangle_{\mathcal{E}}$ is a sequence
$\mathcal{S}^{\mathcal{D}}_{c}\in\Pi(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})^{\omega}$
of distributions on $(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$ s.t.
$\mathcal{S}^{\mathcal{D}}_{c}=\mathcal{S}^{\mathcal{D}}_{c,0}\dots\mathcal{S}^{\mathcal{D}}_{c,n}\dots$
if and only if for all $i\geq 0$ and for all
$\mathbb{D}\in{\mathcal{B}}_{\mathcal{D}}$,

$\mathcal{S}^{\mathcal{D}}_{c,i}(\mathbb{D})=\mathcal{S}_{c,i}^{\mathcal{C}}(\langle\mathcal{P},\mathbb{D}\rangle_{\mathcal{E}}).$

##  4\. Towards a metric for systems

As outlined in the Introduction, we aim at defining a _distance_ over the
systems described in the previous section, called the _evolution metric_ ,
allowing us to do the following:

  1. (1)

Measure how well a program is fulfilling its tasks.

  2. (2)

Establish whether one program behaves better than another one in an
environment.

  3. (3)

Compare the interactions of a program with different environments.

As we will see, these three objectives can be naturally obtained thanks to the
possibility of modelling the program in isolation from the environment typical
of our model, and to our purely data-driven system semantics. Intuitively,
since the behaviour of a system is entirely described by its evolution
sequence, the evolution metric $\mathfrak{m}$ will indeed be defined as a
distance on the evolution sequences of systems. However, in order to obtain
the proper technical definition of $\mathfrak{m}$, some considerations are
due.

Firstly, we notice that in most applications, not only the program-environment
interplay, but also the tasks of the program can be expressed in a purely
data-driven fashion. For instance, if we are monitoring a particular feature
of the system then we can identify a (time-dependent) _specific value of
interest_ of a state variable, like an upper bound on the energy consumption
at a given time step. Similarly, if we aim at controlling several interacting
features or there is a range of parameters that can be considered acceptable,
like in the case of our running example, then we can specify a (time-
dependent) _set of values of interest_. At any time step, any difference
between them and the data actually obtained can be interpreted as a flaw in
systems behaviour. We use a _penalty function_ $\rho$ to quantify these
differences. From the penalty function we can obtain a _distance on data
states_ , namely a $1$-bounded _hemimetric_ $m^{\mathcal{D}}$ expressing how
much a data state $\mathbf{d}_{2}$ is worse than a data state $\mathbf{d}_{1}$
according to parameters of interests. The idea is that if $\mathbf{d}_{1}$ is
obtained by the interaction of $P_{1}$ with $\mathcal{E}$ and $\mathbf{d}_{2}$
is obtained from the one of $P_{2}$ with $\mathcal{E}$, $\mathbf{d}_{2}$ being
worse than $\mathbf{d}_{1}$ means that $P_{2}$ is performing worse than
$P_{1}$ in this case. Similarly, if $\mathbf{d}_{1}$ and $\mathbf{d}_{2}$ are
obtained from the interaction of $P$ with $\mathcal{E}_{1}$ and
$\mathcal{E}_{2}$, respectively, we can express whether in this particular
configuration and time step $P$ is behaving worse with one environment with
respect to the other.

Secondly, we recall that the evolution sequence of a system consists in a
sequence of _distributions_ over data states. Hence, we use the _Wasserstein
metric_ to lift $m^{\mathcal{D}}$ to a distance
$\operatorname{\mathbf{W}}(m^{\mathcal{D}})$ over distributions over data
states. Informally, the Wasserstein metric gives the expected value of the
ground distance between the elements in the support of the distributions.
Thus, in our case, the lifting expresses how much worse a configuration is
expected to behave with respect to another one at a given time.

Finally, we need to lift $\operatorname{\mathbf{W}}(m^{\mathcal{D}})$ to a
distance on the entire evolution sequences of systems. For our purposes, a
reasonable choice is to take the maximum over time of the pointwise (with
respect to time) Wasserstein distances (see Remark 9 below for further details
on this choice). Actually, to favour computational tractability, our metric
will _not_ be evaluated on _all_ the distributions in the evolution sequences,
but _only_ on those that are reached at certain time steps, called the
_observation times_ ($\mathrm{OT}$).

We devote the remaining of this section to a formalisation of the intuitions
above.

### A metric on data states.

We start by proposing a metric on data states, seen as _static components_ in
isolation from processes and environment. To this end, we introduce a _penalty
function_ $\rho\colon\mathcal{D}\to[0,1]$, a continuous function that assigns
to each data state $\mathbf{d}$ a penalty in $[0,1]$ expressing how far the
values of the parameters of interest in $\mathbf{d}$ are from their desired
ones (hence $\rho(\mathbf{d})=0$ if $\mathbf{d}$ respects all the parameters).
Since some parameters can be time-dependent, so is $\rho$: at any time step
$\tau$, the $\tau$-penalty function $\rho_{\tau}$ compares the data states
with respect to the values of the parameters expected at time $\tau$.

A requirement on the behaviour of the three-tanks system from Example 3 can be
that the level of water $l_{i}$ should be at the level $l_{goal}$, for
$i=1,2,3$. This can be easily rendered in terms of the penalty functions
$\rho^{l_{i}}$, for $i=1,2,3$, defined as the normalised distance between the
current level of water $\mathbf{d}(l_{i})$ and $l_{goal}$, namely:

$\rho^{l_{i}}(\mathbf{d})=\frac{|\mathbf{d}(l_{i})-l_{goal}|}{\max\\{l_{max}-l_{goal},l_{goal}-l_{min}\\}}\enspace.$
(9)

Clearly, we can also consider as a requirement that _all_ the $l_{i}$ are at
level $l_{goal}$. The penalty function thus becomes

$\rho^{\max}(\mathbf{d})=\max_{i=1,2,3}\rho^{l_{i}}(\mathbf{d})\enspace.$

The ones proposed above are just a few simple examples of requirements, and
related penalty functions, for the three-tanks experiment. Yet, in their
simplicity, they will allow us to showcase a number of different applications
of our framework while remaining in a setting where the reader can easily
foresee the expected results.  $\LHD$

A formal definition of the penalty function is beyond the purpose of this
paper, also due to its context-dependent nature. Besides, notice that we can
assume that $\rho$ already includes some tolerances with respect to the exact
values of the parameters in its evaluation, and thus we do not consider them.
The (_timed_) _metric on data states_ is then defined as the asymmetric
difference between the penalties assigned to them by the penalty function.

[Metric on data states] For any time step $\tau$, let
$\rho_{\tau}\colon\mathcal{D}\rightarrow[0,1]$ be the $\tau$-penalty function
on $\mathcal{D}$. The $\tau$-_metric on data states_ in $\mathcal{D}$,
$m_{\rho,\tau}^{\mathcal{D}}\colon\mathcal{D}\times\mathcal{D}\to[0,1]$, is
defined, for all $\mathbf{d}_{1},\mathbf{d}_{2}\in\mathcal{D}$, by

$m_{\rho,\tau}^{\mathcal{D}}(\mathbf{d}_{1},\mathbf{d}_{2})=\max\\{\rho_{\tau}(\mathbf{d}_{2})-\rho_{\tau}(\mathbf{d}_{1}),0\\}.$

Notice that $m_{\rho,\tau}^{\mathcal{D}}(\mathbf{d}_{1},\mathbf{d}_{2})>0$ iff
$\rho_{\tau}(\mathbf{d}_{2})>\rho_{\tau}(\mathbf{d}_{1})$, i.e., the penalty
assigned to $\mathbf{d}_{2}$ is higher than that assigned to $\mathbf{d}_{1}$.
For this reason, we say that
$m_{\rho,\tau}^{\mathcal{D}}(\mathbf{d}_{1},\mathbf{d}_{2})$ expresses _how
worse_ $\mathbf{d}_{2}$ is than $\mathbf{d}_{1}$ with respect to the
objectives of the system. It is not hard to see that for all
$\mathbf{d}_{1},\mathbf{d}_{2},\mathbf{d}_{3}\in\mathcal{D}$ we have
$m_{\rho,\tau}^{\mathcal{D}}(\mathbf{d}_{1},\mathbf{d}_{2})\leq 1$,
$m_{\rho,\tau}^{\mathcal{D}}(\mathbf{d}_{1},\mathbf{d}_{1})=0$, and
$m_{\rho,\tau}^{\mathcal{D}}(\mathbf{d}_{1},\mathbf{d}_{2})\leq
m_{\rho,\tau}^{\mathcal{D}}(\mathbf{d}_{1},\mathbf{d}_{3})+m_{\rho,\tau}^{\mathcal{D}}(\mathbf{d}_{3},\mathbf{d}_{2})$,
thus ensuring that $m_{\rho,\tau}^{\mathcal{D}}$ is a 1-bounded hemimetric.

######  Proposition 6.

Function $m_{\rho,\tau}^{\mathcal{D}}$ is a 1-bounded hemimetric on
$\mathcal{D}$.

When no confusion shall arise, we shall drop the $\rho,\tau$ subscripts. We
remark that penalty functions allow us to define the distance between two data
states, which are elements of n, in terms of a distance on . As we discuss in
Section 5, this feature significantly lowers the complexity of the evaluation
the evolution metric.

### Lifting $m^{\mathcal{D}}$ to distributions.

The second step to obtain the evolution metric consists in lifting
$m^{\mathcal{D}}$ to a metric on distributions on data states. In the
literature, we can find a wealth of notions of distances on distributions (see
[RKSF13] for a survey). For our purposes, the most suitable one is the
_Wasserstein metric_ [Vas69].

According to Definition 2 (given in Section 2), for any two distributions
$\mu,\nu$ on $(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$, the Wasserstein
lifting of $m^{\mathcal{D}}$ to a distance between $\mu$ and $\nu$ is defined
by

$\operatorname{\mathbf{W}}(m^{\mathcal{D}})(\mu,\nu)=\inf_{\mathfrak{w}\in\mathfrak{W}(\mu,\nu)}\int_{\mathcal{D}\times\mathcal{D}}m^{\mathcal{D}}(\mathbf{d},\mathbf{d}^{\prime})\;\mathrm{d}\mathfrak{w}(\mathbf{d},\mathbf{d}^{\prime})$

where $\mathfrak{W}(\mu,\nu)$ is the set of the _couplings of $\mu$ and
$\nu$_, namely the set of joint distributions $\mathfrak{w}$ over the product
space
$(\mathcal{D}\times\mathcal{D},{\mathcal{B}}_{\mathcal{D}}\otimes{\mathcal{B}}_{\mathcal{D}})$
having $\mu$ and $\nu$ as left and right marginal.

As outlined in Section 2, [FR18, Proposition 4, Proposition 6] guarantee that
$\operatorname{\mathbf{W}}(m^{\mathcal{D}})$ is a well-defined hemimetric on
$\Pi(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$.

######  Proposition 7.

Function $\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,\tau})$ is a
$1$-bounded hemimetric on $\Pi(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$.

### The evolution metric.

We now need to lift $\operatorname{\mathbf{W}}(m^{\mathcal{D}})$ to a distance
on evolution sequences. To this end, we observe that the evolution sequence of
a configuration includes the distributions over data states induced after
_each_ computation step. However, it could be the case that the changes on
data induced by the environment can be appreciated only along wider time
intervals. To deal with this kind of situations, we introduce the notion of
_observation times_ , namely a _discrete_ set $\mathrm{OT}$ of time steps at
which the modifications induced by the program-environment interplay give us
useful information on the evolution of the system (with respect to the
considered objectives). Hence, a comparison of the evolution sequences based
on the differences in the distributions reached at the times in $\mathrm{OT}$
can be considered meaningful. Moreover, considering only the differences at
the observation times will favour the computational tractability of the
evolution metric.

To compare evolution sequences, we propose a sort of _weighted infinity norm_
of the tuple of the Wasserstein distances between the distributions in them.
As weight we consider a non-increasing function
$\lambda\colon\mathrm{OT}\to(0,1]$ allowing us to express how much the
distance at time $\tau$, namely
$\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,\tau})(\mathcal{S}^{\mathcal{D}}_{c_{1},\tau},\mathcal{S}^{\mathcal{D}}_{c_{2},\tau})$,
affects the overall distance between configurations $c_{1}$ and $c_{2}$. The
idea is that, in certain application contexts, the differences in the
behaviour of systems that are detected after wide time intervals, can be less
relevant than the differences in the initial behaviour. Hence, we can use the
weight $\lambda(\tau)$ to _mitigate_ the distance of events at time $\tau$.
Following the terminology introduced in the context of behavioural metrics
[dAHM03, DGJP04], we refer to $\lambda$ as to the _discount function_ , and to
$\lambda(\tau)$ as to the _discount factor at time_ $\tau$. Clearly, the
constant function $\lambda(\tau)=1$, for all $\tau\in\mathrm{OT}$, means that
no discount is applied.

[Evolution metric] Assume a set $\mathrm{OT}$ of observation times and a
discount function $\lambda$. For each $\tau\in\mathrm{OT}$, let $\rho_{\tau}$
be a $\tau$-penalty function and let $m^{\mathcal{D}}_{\rho,\tau}$ be the
$\tau$-metric on data states defined on it. Then, the $\lambda$-_evolution
metric_ over $\rho$ and $\mathrm{OT}$ is the mapping
$\mathfrak{m}^{\lambda}_{\rho,\mathrm{OT}}\colon\mathcal{C}\times\mathcal{C}\to[0,1]$
defined, for all configurations $c_{1},c_{2}\in\mathcal{C}$, by

$\mathfrak{m}^{\lambda}_{\rho,\mathrm{OT}}(c_{1},c_{2})=\sup_{\tau\in\mathrm{OT}}\,\lambda(\tau)\cdot\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,\tau})\Big{(}\mathcal{S}^{\mathcal{D}}_{c_{1},\tau},\mathcal{S}^{\mathcal{D}}_{c_{2},\tau}\Big{)}.$

Since $\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,\tau})$ is a 1-bounded
hemimetric (Proposition 7) we can easily derive the same property for
$\mathfrak{m}^{\lambda}_{\rho,\mathrm{OT}}$.

######  Proposition 8.

Function $\mathfrak{m}^{\lambda}_{\rho,\mathrm{OT}}$ is a 1-bounded hemimetric
on $\mathcal{C}$.

Notice that if $\lambda$ is a _strictly_ non-increasing function, then to
obtain upper bounds on the evolution metric only a _finite_ number of
observations is needed.

######  Remark 9.

Usually, due to the presence of uncertainties, the behaviour of a system can
be considered acceptable even if it differs from its intended one _up-to a
certain tolerance_. Similarly, the properties of adaptability, reliability,
and robustness that we aim to study will check whether a program is able to
perform well in a perturbed environment _up-to a given tolerance_. In this
setting, the choice of defining the evolution metric as the pointwise maximal
distance in time between the evolution sequences of systems is natural and
reasonable: if in the worst case (the maximal distance) the program keeps the
parameters of interest within the given tolerance, then its entire behaviour
can be considered acceptable. However, with this approach we have that a
program is only as good as its worst performance, and one could argue that
there are application contexts in which our evolution metric would be less
meaningful. For these reasons, we remark that we could have given a
_parametric_ version of Definition 4 and defining the evolution metric in
terms of a generic _aggregation function_ $f$ over the tuple of Wasserstein
distances. Then, one could choose the best instantiation for $f$ according to
the chosen application context. The use of a parametric definition would have
not affected the technical development of our paper. However, in order to keep
the notation and presentation as simple as possible, we opted to define
$\mathfrak{m}^{\lambda}_{\rho,\mathrm{OT}}$ directly in the weighted infinity
norm form. A similar reasoning applies to the definition of the penalty
function that we gave in Example 4, and to the definition of the metric over
data states (Definition 4).

##  5\. Estimating the evolution metric

In this section we show how the evolution metric can be estimated via
statistical techniques. Firstly, in Section 5.1 we show how we can estimate
the evolution sequence of a given configuration $c$. Then, in Section 5.2 we
evaluate the distance between two configurations $c_{1}$ and $c_{2}$ on their
estimated evolution sequences.

###  5.1. Computing empirical evolution sequences

Given a configuration $\langle P,\mathbf{d}\rangle_{\mathcal{E}}$ and an
integer $k$ we can use function Simulate, defined in Figure 2, to sample a
sequence of configurations of the form

$\langle P_{0},\mathbf{d}_{0}\rangle_{\mathcal{E}},\langle
P_{1},\mathbf{d}_{1}\rangle_{\mathcal{E}},\ldots,\langle
P_{k},\mathbf{d}_{k}\rangle_{\mathcal{E}}.$

This sequence represents $k$-steps of a computation starting from $\langle
P,\mathbf{d}\rangle_{\mathcal{E}}=\langle
P_{0},\mathbf{d}_{0}\rangle_{\mathcal{E}}$. Each step of the sequence is
computed by using function SimStep, also defined in Figure 2. There we let
rand be a function that allows us to get a uniform random number in $(0,1]$
while $\textsc{sample}(\mathcal{E},\mathbf{d})$ is used to sample a data state
in the distribution $\mathcal{E}(\mathbf{d})$. We assume that for any
measurable set of data states $\mathbb{D}\in{\mathcal{B}}_{\mathcal{D}}$ the
probability of $\textsc{sample}(\mathcal{E},\mathbf{d})$ to be in $\mathbb{D}$
is equal to the probability of $\mathbb{D}$ induced by
$\mathcal{E}(\mathbf{d})$, namely:

$\mathit{Pr}(\textsc{sample}(\mathcal{E},\mathbf{d})\in\mathbb{D})=\mathcal{E}(\mathbf{d})(\mathbb{D})\enspace.$
(10)

If we consider the environment evolution from Examples 3 and 3, function
$\textsc{sample}(\mathcal{E},\mathbf{d})$ would sample a random value $x$ for
$q_{2}$ from the distribution chosen for it (according to the considered
scenario), and then solve the system of equations identified by
$f_{3T}(x,\mathbf{d})$.  $\LHD$

1:function Simulate($\langle P,\mathbf{d}\rangle_{\mathcal{E}},k$)

2: $i\leftarrow 1$

3: $c\leftarrow\langle P,\mathbf{d}\rangle_{\mathcal{E}}$

4: $a\leftarrow c$

5: while $i\leq k$ do

6: $c\leftarrow\textsc{SimStep}(c)$

7: $a\leftarrow a,c$

8: $i\leftarrow i+1$

9: end while

10: return $a$

11:end function

1:function SimStep($\langle P,\mathbf{d}\rangle_{\mathcal{E}}$)

2:
$\sum_{i=1}^{n}p_{i}(\theta_{i},P_{i})\leftarrow\mathsf{pstep}(P,\mathbf{d})$

3: $u\leftarrow\textsc{rand}()$

4: let $i$ s.t. $\sum_{j=1}^{i-1}p_{j}<u\leq\sum_{j=1}^{i}p_{j}$

5:
$\mathbf{d}_{i}^{\prime}\leftarrow\textsc{sample}(\mathcal{E},\theta_{i}(\mathbf{d}_{i}))$

6: return $\langle P_{i},\mathbf{d}_{i}^{\prime}\rangle_{\mathcal{E}}$

7:end function

Figure 2. Functions used to simulate behaviour of a configuration.

We can then extend this property to configurations and function SimStep:

######  Lemma 10.

For any configuration $\langle
P,\mathbf{d}\rangle_{\mathcal{E}}\in\mathcal{C}$, and for any measurable set
$\mathbb{C}\in\Sigma_{\mathcal{C}}$:

$\mathit{Pr}(\textsc{SimStep}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}})\in\mathbb{C})=\mathsf{cstep}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}})(\mathbb{C}).$

######  Proof 5.1.

To prove the thesis it is enough to show that the property

$\mathit{Pr}(\textsc{SimStep}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}})\in\mathbb{C})=\mathsf{cstep}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}})(\mathbb{C})$

holds on the generators of the $\sigma$-algebra $\Sigma_{\mathcal{C}}$. Hence,
assume that $\mathbb{C}=\langle\mathbb{P},\mathbb{D}\rangle_{\mathcal{E}}$ for
some $\mathbb{P}\in\Sigma_{\mathcal{P}}$ and
$\mathbb{D}\in{\mathcal{B}}_{\mathcal{D}}$. Then we have

$\displaystyle\mathsf{cstep}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}})(\mathbb{C})=\quad$
$\displaystyle\sum_{(\theta,P^{\prime})\in\mathsf{supp}(\mathsf{pstep}(P,\mathbf{d}))}\left(\mathsf{pstep}(P,\mathbf{d})(\theta,P^{\prime})\cdot\langle
P^{\prime},\mathcal{E}(\theta(\mathbf{d}))\rangle_{\mathcal{E}}\right)(\mathbb{C})$
$\displaystyle=\quad$ $\displaystyle\sum_{i=1}^{n}\left(p_{i}\cdot\langle
P_{i},\mathcal{E}(\theta_{i}(\mathbf{d}))\rangle_{\mathcal{E}}\right)(\mathbb{C})$
$\displaystyle=\quad$
$\displaystyle\sum_{i=1}^{n}p_{i}\cdot\mathds{1}_{\mathbb{P}}\\{P_{i}\\}\cdot\mathcal{E}(\theta_{i}(\mathbf{d}))(\mathbb{D})$
$\displaystyle=\quad$
$\displaystyle\sum_{\\{i|P_{i}\in\mathbb{P}\\}}p_{i}\cdot\mathcal{E}(\theta_{i}(\mathbf{d}))(\mathbb{D})$
$\displaystyle=\quad$
$\displaystyle\sum_{\\{i|P_{i}\in\mathbb{P}\\}}p_{i}\cdot\mathit{Pr}(\textsc{Sample}(\mathcal{E},\theta_{i}(\mathbf{d}))\in\mathbb{D})$
$\displaystyle=\quad$ $\displaystyle\mathit{Pr}(\textsc{SimStep}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}})\in\mathbb{C})$

where:

  * •

The first step follows by the definition of $\mathsf{cstep}$ (Definition 3).

  * •

The second step follows by
$\mathsf{pstep}(P,\mathbf{d})=p_{1}\cdot(\theta_{1},P_{1})+\ldots+p_{n}\cdot(\theta_{n},P_{n})$,
for some $p_{i}$, $\theta_{i}$ and $P_{i}$ (see Proposition 2).

  * •

The third and fourth steps follow by the definition of the distribution
$\langle\mu_{\mathcal{P}},\mu_{\mathcal{D}}\rangle_{\mathcal{E}}$ (Notation
3), in which $\mathds{1}_{\mathbb{P}}$ denotes the characteristic function of
the set $\mathbb{P}$, i.e., $\mathds{1}_{\mathbb{P}}\\{P^{\prime}\\}=1$ if
$P^{\prime}\in\mathbb{P}$ and $\mathds{1}_{\mathbb{P}}\\{P^{\prime}\\}=0$
otherwise.

  * •

The fifth step follows by the assumption on the function
$\textsc{Sample}(\\_)$ in Equation (10).

  * •

The last step follows by the definition of function $\textsc{SimpStep}(\\_)$.

1:function Estimate($\langle P,\mathbf{d}\rangle_{\mathcal{E}},k,N$)

2: $\forall i:(0\leq i\leq k):E_{i}\leftarrow\emptyset$

3: $counter\leftarrow 0$

4: while $counter<N$ do

5: $(c_{0},\ldots,c_{k})\leftarrow\textsc{Simulate}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}},k)$

6: $\forall i:E_{i}\leftarrow E_{i},c_{i}$

7: $counter\leftarrow counter+1$

8: end while

9: return $E_{0},\ldots,E_{k}$

10:end function

Figure 3.  Function used to obtain $N$ samples of the evolution sequence of a
configuration.

To compute the empirical evolution sequence of a configuration $c$ the
function Estimate in Figure 3 can be used.

Function $\textsc{Estimate}(c,k,N)$ invokes $N$ times function Simulate in
Figure 2 to sample a sequence of configurations $c_{0},\dots,c_{k}$, modelling
$k$ steps of a computation from $c=c_{0}$. Then, a sequence of observations
$E_{0},\ldots,E_{k}$ is computed, where each $E_{i}$ is the tuple
$c_{i}^{1},\ldots,c_{i}^{N}$ of configurations observed at time $i$ in each of
the $N$ sampled computations.

Notice that, for each $i\in\\{0,\dots,k\\}$, the samples
$c_{i}^{1},\dots,c_{i}^{N}$ are independent and moreover, by Lemma 10, they
are identically distributed. Each $E_{i}$ can be used to estimate the
distribution $\mathcal{S}^{\mathcal{C}}_{c,i}$. For any $i$, with $0\leq i\leq
k$, we let $\hat{\mathcal{S}}_{c,i}^{\mathcal{C},N}$ be the distribution such
that for any measurable set of configurations
$\mathbb{C}\in\Sigma_{\mathcal{C}}$ we have

$\hat{\mathcal{S}}_{c,i}^{\mathcal{C},N}(\mathbb{C})=\frac{|E_{i}\cap\mathbb{C}|}{N}.$

Finally, we let
$\hat{\mathcal{S}}_{c}^{\mathcal{D},N}=\hat{\mathcal{S}}_{c,0}^{\mathcal{D},N}\ldots\hat{\mathcal{S}}_{c,k}^{\mathcal{D},N}$
be the _empirical evolution sequence_ such that for any measurable set of data
states $\mathbb{D}\in{\mathcal{B}}_{\mathcal{D}}$ we have

$\hat{\mathcal{S}}_{c,i}^{\mathcal{D},N}(\mathbb{D})=\hat{\mathcal{S}}_{c,i}^{\mathcal{C},N}(\langle\mathcal{P},\mathbb{D}\rangle_{\mathcal{E}}).$

Then, by applying the weak law of large numbers to the i.i.d samples, we get
that when $N$ goes to infinite both $\hat{\mathcal{S}}_{c,i}^{\mathcal{C},N}$
and $\hat{\mathcal{S}}_{c,i}^{\mathcal{D},N}$ converge weakly to
$\mathcal{S}^{\mathcal{C}}_{c,i}$ and $\mathcal{S}_{c,i}^{\mathcal{D}}$
respectively:

$\lim_{N\rightarrow\infty}\hat{\mathcal{S}}_{c,i}^{\mathcal{C},N}=\mathcal{S}_{c,i}^{\mathcal{C}}\qquad\qquad\lim_{N\rightarrow\infty}\hat{\mathcal{S}}_{c,i}^{\mathcal{D},N}=\mathcal{S}_{c,i}^{\mathcal{D}}.$
(11)

The algorithms in Figures 2 and 3 have been implemented in Python to support
analysis of systems. The tool is available at
https://github.com/quasylab/spear.

Figure 4.  Simulation results. The dashed red line corresponds to
$l_{goal}=10$.

Figure 5.  Simulation of the variation of the level of water in each tank in
the two scenarios. The dashed green line corresponds to $l_{goal}=10$.

Figure 6.  Estimated distributions of the water levels in the tanks in the two
scenarios.

A simulation of the three-tanks laboratory experiment is given in Figure 4,
with the following setting: [(i)]

$l_{min}=0$,

$l_{max}=20$,

$l_{goal}=10$,

$\delta_{l}=0.5$,

$q_{max}=6$,

$q_{step}=q_{max}/5$,

$\delta_{q}=0.5$,

$\Delta\tau=0.1$. Moreover, we assume an initial data state $\mathbf{d}_{0}$
with $l_{i}=l_{min}$ and $q_{i}=0$, for $i=1,2,3$. The plot on the left hand
side of Figure 4 depicts a simulation run of the system $\mathbf{s}_{1}$,
namely the system in which the environment evolution follows scenario 1 and
having
$c_{1}=\langle\mathsf{P_{Tanks}},\mathbf{d}_{0}\rangle_{\mathcal{E}_{1}}$ as
initial configuration. On the right hand side we have the corresponding plot
for system $\mathbf{s}_{2}$, whose environment evolution follows scenario 2
and the initial configuration is thus
$c_{2}=\langle\mathsf{P_{Tanks}},\mathbf{d}_{0}\rangle_{\mathcal{E}_{2}}$.

In Figure 5 we consider the same simulations, but we compare the variations of
the level of water, in time, in the same tank in the two scenarios. This can
help us to highlight the potential differences in the evolution of the system
with respect to the two environments.

In Figure 6 we give an estimation of the distributions of $l_{1},l_{2},l_{3}$
obtained from our simulations of $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$. In
detail, the three plots in the top row report the comparison of the
distributions over $l_{1}$ obtained in the two scenarios at time step,
respectively, $50,100$ and $150$ when using $N=1000$ samples. The three plots
in the central row are the corresponding ones for $l_{2}$, and in the bottom
row we have the plots related to $l_{3}$. In both scenarios, we obtain
Gaussian-like distributions, as expected given the probabilistic behaviour of
the environment defined in Examples 3 and 3. However, we observe that while in
both cases the mean of the distributions is close to $l_{goal}=10$, there is a
significant difference in their variance: the variance of the distributions
obtained in scenario 1 is generally smaller than that of those in scenario 2
(the curves for scenario 1 are “slimmer” than those for scenario 2). We will
see how this disparity is captured by the evolution metric.  $\LHD$

As usual in the related literature, we can apply the classic _standard error
approach_ to analyse the approximation error of our statistical estimation of
the evolution sequences. Briefly, we let $x\in\\{l_{1},l_{2},l_{3}\\}$ and we
focus on the distribution of the means of our samples (in particular we
consider the cases $N=1000,5000,10000$) for each variable $x$. In each case,
for each $i\in\\{0,\dots,k\\}$ we compute the mean
$\overline{E_{i}}(x)=\frac{1}{N}\sum_{j=1}^{N}\mathbf{d}_{i}^{j}(x)$ of the
sampled data (recall that $E_{i}$ is the tuple of the sampled configurations
$c_{i}^{1},\dots,c_{i}^{N}$, with $c_{i}^{j}=\langle
P_{i}^{j},\mathbf{d}_{i}^{j}\rangle_{\mathcal{E}}$), and we evaluate their
_standard deviation_
$\tilde{\sigma}_{i,N}(x)=\sqrt{\frac{\sum_{j=1}^{N}(\mathbf{d}_{i}^{j}(x)-\overline{E_{i}}(x))^{2}}{N-1}}$
(see Figure 7(a) for the variation in time of the standard deviation of the
distribution of $x=l_{3}$). From $\tilde{\sigma}_{i,N}(x)$ we obtain the
_standard error of the mean_
$\overline{\sigma}_{i,N}(x)=\frac{\tilde{\sigma}_{i,N}(x)}{\sqrt{N}}$ (see
Figure 7(b) for the variation in time of the standard error for the
distribution of $x=l_{3}$). Finally, we proceed to compute the $z$-_score_ of
our sampled distribution as follows:
$z_{i,N}(x)=\frac{\overline{E_{i}}(x)-\mathbb{E}(x)}{\overline{\sigma}_{i,N}(x)}$,
where $\mathbb{E}(x)$ is the mean (or expected value) of the real distribution
over $x$. In Figure 7(c) we report the variation in time of the $z$-score of
the distribution over $x=l_{3}$: the dashed red lines correspond to $z=\pm
1.96$, namely the value of the $z$-score corresponding to a confidence
interval of the $95\%$. We can see that our results can already be given with
a $95\%$ confidence in the case of $N=1000$. Please notice that the
oscillation in time of the values of the $z$-scores is due to the
perturbations introduced by the environment in the simulations and by the
natural oscillation in the interval
$[l_{goal}-\delta_{l},l_{goal}+\delta_{l}]$ of the water levels in the
considered experiment (see Figure 4). A similar analysis, with analogous
results, can be carried out for the distributions of $l_{1}$ and $l_{2}$. In
Figure 7(d) we report the variation in time of the $z$-scores of the
distributions of the three variables, in the case $N=1000$.

(a) Standard deviation for $l_{3}$.

(b) Standard error for $l_{3}$.

(c) $z$-score for $l_{3}$.

(d) $z$-scores for $l_{1},l_{2},l_{3}$, $N=1000$.

Figure 7. Analysis of the approximation error, over time, of the distributions
of $l_{1},l_{2},l_{3}$, in scenario 2, for $N=1000,5000,10000$.

###  5.2. Computing distance between two configurations

Function Estimate allows us to collect independent samples at each time step
$i$ from $0$ to a deadline $k$. These samples can be used to estimate the
distance between two configurations $c_{1}$ and $c_{2}$. Following an approach
similar to the one presented in [TK10], to estimate the Wasserstein distance
$\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,i})$ between two (unknown)
distributions $\mathcal{S}^{\mathcal{D}}_{c_{1},i}$ and
$\mathcal{S}^{\mathcal{D}}_{c_{2},i}$ we can use $N$ independent samples
$\\{c^{1}_{1},\ldots,c^{N}_{1}\\}$ taken from
$\mathcal{S}^{\mathcal{C}}_{c_{1},i}$ and $\ell\cdot N$ independent samples
$\\{c^{1}_{2},\ldots,c_{2}^{\ell\cdot N}\\}$ taken from
$\mathcal{S}^{\mathcal{C}}_{c_{2},i}$. After that, we exploit the $i$-penalty
function $\rho_{i}$ and we consider the two sequences of values
$\omega_{1},\dots,\omega_{N}$ and $\nu_{1},\dots,\nu_{\ell\cdot N}$ defined as
follows:

$\\{\omega_{j}=\rho_{i}(\mathbf{d}_{1}^{j})\mid\langle
P_{1}^{j},\mathbf{d}_{1}^{j}\rangle_{\mathcal{E}_{1}}=c_{1}^{j}\\}\qquad\\{\nu_{h}=\rho_{i}(\mathbf{d}_{2}^{h})\mid\langle
P^{h}_{2},\mathbf{d}^{h}_{2}\rangle_{\mathcal{E}_{2}}=c_{2}^{h}\\}.$

We can assume, without loss of generality that these sequences are ordered,
i.e., $\omega_{j}\leq\omega_{j+1}$ and $\nu_{h}\leq\nu_{h+1}$. The value
$\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,i})(\mathcal{S}^{\mathcal{D}}_{c_{1},i},\mathcal{S}^{\mathcal{D}}_{c_{2},i})$
can be approximated as $\frac{1}{\ell N}\sum_{h=1}^{\ell
N}\max\\{\nu_{h}-\omega_{\lceil\frac{h}{\ell}\rceil},0\\}$. The next theorem,
based on results in [TK10, Vil08, Val74], ensures that the larger the number
of samplings the closer the gap between the estimated value and the exact one.

######  Theorem 11.

Let
$\mathcal{S}^{\mathcal{C}}_{c_{1},i},\mathcal{S}^{\mathcal{C}}_{c_{2},i}\in\Pi(\mathcal{C},\Sigma_{\mathcal{C}})$
be unknown. Let $\\{c^{1}_{1},\ldots,c^{N}_{1}\\}$ be independent samples
taken from $\mathcal{S}^{\mathcal{C}}_{c_{1},i}$, and
$\\{c^{1}_{2},\ldots,c_{2}^{\ell\cdot N}\\}$ independent samples taken from
$\mathcal{S}^{\mathcal{C}}_{c_{2},i}$. Let
$\\{\omega_{j}=\rho_{i}(\mathbf{d}^{j}_{1})\\}$ and
$\\{\nu_{h}=\rho_{i}(\mathbf{d}^{h}_{2})\\}$ be the ordered sequences obtained
from the samples and the $i$-penalty function. Then, the Wasserstein distance
between $\mathcal{S}^{\mathcal{D}}_{c_{1},i}$ and
$\mathcal{S}^{\mathcal{D}}_{c_{2},i}$ is equal, a.s., to

$\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,i})(\mathcal{S}^{\mathcal{D}}_{c_{1},i},\mathcal{S}^{\mathcal{D}}_{c_{2},i})=\lim_{N\rightarrow\infty}\frac{1}{\ell
N}\sum_{h=1}^{\ell
N}\max\left\\{\nu_{h}-\omega_{\lceil\frac{h}{\ell}\rceil},0\right\\}.$

######  Proof 5.2.

Let $c_{1}^{j}=\langle P^{j}_{1},\mathbf{d}^{j}_{1}\rangle_{\mathcal{E}}$, for
all $j=1,\dots,N$, and $c_{2}^{j}=\langle
P^{j}_{2},\mathbf{d}^{j}_{2}\rangle_{\mathcal{E}}$, for all
$j=1,\dots,\ell\cdot N$. We split the proof into two parts showing
respectively:

$\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,i})(\mathcal{S}^{\mathcal{D}}_{c_{1},i},\mathcal{S}^{\mathcal{D}}_{c_{2},i})=\lim_{N\to\infty}\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,i})(\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i},\hat{\mathcal{S}}^{\mathcal{D},\ell
N}_{c_{2},i})\enspace.$ (12)
$\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,i})(\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i},\hat{\mathcal{S}}^{\mathcal{D},\ell
N}_{c_{2},i})=\frac{1}{\ell N}\sum_{h=1}^{\ell
N}\max\left\\{\nu_{h}-\omega_{\lceil\frac{h}{\ell}\rceil},0\right\\}\enspace.$
(13)

  * •

Proof of Equation (12).

We recall that the sequence
$\\{\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{l},i}\\}$ converges weakly to
$\mathcal{S}^{\mathcal{D}}_{c_{l},i}$ for $l\in\\{1,2\\}$ (see Equation (11)).
Moreover, we can prove that these sequences converge weakly in
$\Pi(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$ in the sense of [Vil08,
Definition 6.8]. In fact, given the $i$-ranking function $\rho_{i}$, the
existence of a data state $\tilde{\mathbf{d}}$ such that
$\rho_{i}(\tilde{\mathbf{d}})=0$ is guaranteed (remember that the constraints
used to define $\rho_{i}$ are on the possible values of state variables and a
data state fulfilling all the requirements is assigned value $0$). Thus, for
any $\mathbf{d}\in\mathcal{D}$ we have that

$m^{\mathcal{D}}_{\rho,i}(\tilde{\mathbf{d}},\mathbf{d})=\max\\{\rho_{i}(\mathbf{d})-\rho_{i}(\tilde{\mathbf{d}}),0\\}=\rho_{i}(\mathbf{d})\enspace.$

Since, moreover, by definition $\rho_{i}$ is continuous and bounded, the weak
convergence of the distributions gives

$\displaystyle\int_{\mathcal{D}}\rho_{i}(\mathbf{d})\mathrm{d}(\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i}(\mathbf{d}))\to\int_{\mathcal{D}}\rho_{i}(\mathbf{d})\mathrm{d}(\mathcal{S}^{\mathcal{D}}_{c_{1},i}(\mathbf{d}))$
$\displaystyle\int_{\mathcal{D}}\rho_{i}(\mathbf{d})\mathrm{d}(\hat{\mathcal{S}}^{\mathcal{D},\ell
N}_{c_{2},i}(\mathbf{d}))\to\int_{\mathcal{D}}\rho_{i}(\mathbf{d})\mathrm{d}(\mathcal{S}^{\mathcal{D}}_{c_{2},i}(\mathbf{d}))$

and thus Definition 6.8.(i) of [Vil08] is satisfied. As $\mathcal{D}$ is a
Polish space, by [Vil08, Theorem 6.9] we obtain that

$\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{l},i}\to\mathcal{S}^{\mathcal{D}}_{c_{l},i}\text{
implies
}\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,i})(\mathcal{S}^{\mathcal{D}}_{c_{1},i},\mathcal{S}^{\mathcal{D}}_{c_{2},i})=\lim_{N\to\infty}\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,i})(\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i},\hat{\mathcal{S}}^{\mathcal{D},\ell
N}_{c_{2},i})\enspace.$

  * •

Proof of Equation (13).

For this part of the proof we follow [TK10]. Since the ranking function is
continuous, it is in particular ${\mathcal{B}}_{\mathcal{D}}$ measurable and
therefore for any distribution $\mu$ on
$(\mathcal{D},{\mathcal{B}}_{\mathcal{D}})$ we obtain that

$F_{\mu,\rho_{i}}(r):=\mu(\\{\rho_{i}(\mathbf{d})<r\\})$

is a well defined cumulative distribution function. In particular, for
$\mu=\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i}$ we have that

$F_{\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i},\rho_{i}}(r)=\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i}(\\{\rho_{i}(\mathbf{d})<r\\})=\frac{\left|\\{\langle
P_{1}^{j},\mathbf{d}_{1}^{j}\rangle_{\mathcal{E}_{1}}\in
E_{1,i}\mid\rho_{i}(\mathbf{d}_{1}^{j})<r\\}\right|}{N}\enspace.$

Since, moreover, we can always assume that the values
$\rho_{i}(\mathbf{d}_{1}^{j})$ are sorted, so that
$\rho_{i}(\mathbf{d}_{i}^{j})\leq\rho_{i}(\mathbf{d}_{1}^{j+1})$ for each
$j=1,\dots,N-1$, we can express the counter image of the cumulative
distribution function as

$F^{-1}_{\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i},\rho_{i}}(r)=\rho_{i}(\mathbf{d}_{1}^{j})\text{
whenever }\frac{j-1}{N}<r\leq\frac{j}{N}\enspace.$ (14)

A similar reasoning holds for $F^{-1}_{\hat{\mathcal{S}}^{\mathcal{D},\ell
N}_{c_{2},i},\rho_{i}}(r)$.

Then, by [FR18, Proposition 6.2], for each $N$ we have that

$\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,i})(\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i},\hat{\mathcal{S}}^{\mathcal{D},\ell
N}_{c_{2},i})=\int_{0}^{1}\max\left\\{F^{-1}_{\hat{\mathcal{S}}^{\mathcal{D},\ell
N}_{c_{2},i},\rho_{i}}(r)-F^{-1}_{\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i},\rho_{i}}(r),0\right\\}\mathrm{d}r\enspace.$

Let us now partition the interval $[0,1]$ into $\ell N$ intervals of size
$\frac{1}{\ell N}$, thus obtaining

$\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,i})(\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i},\hat{\mathcal{S}}^{\mathcal{D},\ell
N}_{c_{2},i})=\sum_{h=1}^{\ell N}\left(\int_{\frac{h-1}{\ell
N}}^{\frac{h}{\ell N}}\max\left\\{F^{-1}_{\hat{\mathcal{S}}^{\mathcal{D},\ell
N}_{c_{2},i},\rho_{i}}(r)-F^{-1}_{\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i},\rho_{i}}(r),0\right\\}\mathrm{d}r\right)\enspace.$

From Equation (14), on each interval $(\frac{h-1}{\ell N},\frac{h}{\ell N}]$
it holds that
$F^{-1}_{\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i},\rho_{i}}(r)=\rho_{i}(\mathbf{d}_{1}^{\lceil\frac{h}{\ell}\rceil})$
and $F^{-1}_{\hat{\mathcal{S}}^{\mathcal{D},\ell
N}_{c_{2},i},\rho_{i}}(r)=\rho_{i}(\mathbf{d}_{2}^{h})$. Moreover, both
functions are constant on each the interval so that the value of the integral
is given by the difference multiplied by the length of the interval:

$\displaystyle\operatorname{\mathbf{W}}(m^{\mathcal{D}}_{\rho,i})(\hat{\mathcal{S}}^{\mathcal{D},N}_{c_{1},i},\hat{\mathcal{S}}^{\mathcal{D},\ell
N}_{c_{2},i})={}$ $\displaystyle\sum_{h=1}^{\ell N}\frac{1}{\ell
N}\max\left\\{\rho_{i}(\mathbf{d}_{2}^{h})-\rho_{i}(\mathbf{d}_{1}^{\lceil\frac{h}{\ell}\rceil}),0\right\\}$
$\displaystyle=$ $\displaystyle\sum_{h=1}^{\ell N}\frac{1}{\ell
N}\max\left\\{\nu_{h}-\omega_{\lceil\frac{h}{\ell}\rceil},0\right\\}\enspace.$

By substituting the last equality into Equation (12) we obtain the thesis.

1:function Distance($c_{1},c_{2},\rho,\lambda,\mathrm{OT},N,\ell$)

2: $k\leftarrow\max\\{i\in\mathrm{OT}\\}$

3: $E_{1,1},\ldots,E_{1,k}\leftarrow\textsc{Estimate}(c_{1},k,N)$

4: $E_{2,1},\ldots,E_{2,k}\leftarrow\textsc{Estimate}(c_{2},k,\ell N)$

5: $m\leftarrow 0$

6: for all $i\in\mathrm{OT}$ do

7: $m_{i}\leftarrow\textsc{ComputeW}(E_{1,i},E_{2,i},\rho_{i})$

8: $m\leftarrow\max\\{m,\lambda(i)\cdot m_{i}\\}$

9: end for

10: return $m$

11:end function

1:function ComputeW($E_{1},E_{2},\rho$)

2: $(\langle
P^{1}_{1},\mathbf{d}^{1}_{1}\rangle_{\mathcal{E}_{1}},\ldots,\langle
P^{N}_{1},\mathbf{d}^{N}_{1}\rangle_{\mathcal{E}_{1}})\leftarrow E_{1}$

3: $(\langle
P^{1}_{2},\mathbf{d}^{1}_{2}\rangle_{\mathcal{E}_{2}},\ldots,\langle
P_{2}^{\ell N},\mathbf{d}_{2}^{\ell N}\rangle_{\mathcal{E}_{2}})\leftarrow
E_{2}$

4: $\forall j:(1\leq j\leq N):\omega_{j}\leftarrow\rho(\mathbf{d}^{j}_{1})$

5: $\forall h:(1\leq h\leq\ell N):\nu_{h}\leftarrow\rho(\mathbf{d}^{h}_{2})$

6: re index $\\{\omega_{j}\\}$ s.t. $\omega_{j}\leq\omega_{j+1}$

7: re index $\\{\nu_{h}\\}$ s.t. $\nu_{h}\leq\nu_{h+1}$

8: return $\frac{1}{\ell N}\sum_{h=1}^{\ell
N}\max\\{\nu_{h}-\omega_{\lceil\frac{h}{\ell}\rceil},0\\}$

9:end function

Figure 8. Functions used to estimate the evolution metric on systems.

The procedure outlined above is realised by functions Distance and ComputeW in
Figure 8. The former takes as input the two configurations to compare, the
penalty function (seen as the sequence of the $i$-penalty functions), the
discount function $\lambda$, the set $\mathrm{OT}$ of observation times which
we assume to be bounded, and the parameters $N$ and $\ell$ used to obtain the
samplings of computation.

Function Distance collects the samples $E_{i}$ of possible computations during
the observation period $[0,\max_{\mathrm{OT}}]$, where $\max_{\mathrm{OT}}$
denotes the last observation time.

Then, for each observation time $i\in\mathrm{OT}$, the distance at time $i$ is
computed via the function $\textsc{ComputeW}(E_{1,i},E_{2,i},\rho_{i})$.

Since the penalty function allows us to reduce the evaluation of the
Wasserstein distance in n to its evaluation on , due to the sorting of
$\\{\nu_{h}\mid h\in[1,\dots,\ell N]\\}$ the complexity of function ComputeW
is $O(\ell N\log(\ell N))$ (cf. [TK10]). We refer the interested reader to
[SFG+21, Corollary 3.5, Equation (3.10)] for an estimation of the
approximation error given by the evaluation of the Wasserstein distance over
$N$ samples.

###  5.3. Analysis of the three-tanks experiment

Our aim is now to use the evolution metric, and the algorithms introduced
above, to study various properties of the behaviour of the three-tanks system.
In particular, we consider the following two objectives:

  1. (1)

Comparison of the behaviour of two systems generated from the interaction of
the same program with two different environments, under the same initial
conditions.

  2. (2)

Comparison of the behaviour of two systems generated from the interaction of
two different programs with the same environment, under the same initial
conditions.

####  5.3.1. Same program, different environments

Consider the systems $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ introduced in
Example 5.1. We recall that the initial configuration of $\mathbf{s}_{i}$ is
$c_{i}=\langle\mathsf{P_{Tanks}},\mathbf{d}_{0}\rangle_{\mathcal{E}_{i}}$, for
$i=1,2$. Notice that since $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ are
distinguished only by the environment function, a comparison of their
behaviour by means of the evolution metric allows us to establish in which
scenario the program $\mathsf{P_{Tanks}}$ performs better with respect to the
aforementioned tasks.

(a) Evaluation of
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c_{1},c_{2})(\sim
0.15)$.

(b) Evaluation of
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c_{2},c_{1})(\sim
0.02)$.

Figure 9. Evaluation of the distance between $\mathbf{s}_{1}$ and
$\mathbf{s}_{2}$ (and $\mathbf{s}_{2},\mathbf{s}_{1}$) with respect to the
target $l_{3}=l_{goal}$. A red dot at coordinates $(x,y)$ means that
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}_{x}}(c_{1},c_{2})=y$,
where $\mathrm{OT}_{x}=\\{\tau\in\mathrm{OT}\mid\tau\geq x\\}$. The blue line
is the pointwise Wasserstein distance between the distributions in the
evolution sequences of $c_{i}$ and $c_{3-i}$.

Firstly, we focus on a single tank by considering the target $l_{3}=l_{goal}$,
and thus the penalty function $\rho^{l_{3}}$. We do not consider any discount
(so $\lambda$ is the constant function $\boldsymbol{1}$, i.e.,
$\lambda(\tau)=1$ for all $\tau$) and we let $\mathrm{OT}={\rm
Nature}\cap[0,150]$. In Figure 9 we show the evaluation of
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}$ between
$\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ (Figure 9(a)) and between
$\mathbf{s}_{2}$ and $\mathbf{s}_{1}$ (Figure 9(b)) over a simulation with
$5000$ samples ($N=1000$, $l=5$). Notice that the first distance is
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c_{1},c_{2})\sim
0.15$, whereas the second one is
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c_{2},c_{1})\sim
0.02$. The different scale of the two distances (the latter is less than 1/7
of the former) already gives us an intuition of the advantages the use of a
hemimetric gives us over a (pseudo)metric in terms of expressiveness.

Consider Figure 9(a). Recall that, given the definition of our hemimetric on
data states, the value of
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c_{1},c_{2})$
expresses how much the probabilistic evolution of the data in the evolution
sequence of $c_{2}$ is worse than that of $c_{1}$, with respect to the task
identified by $\rho^{3}$. In other words, the higher the value of
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c_{1},c_{2})$, the
worse the performance of $c_{2}$ with respect to $c_{1}$. Yet, for those who
are not familiar with the Wasserstein distance, the plot alone can be a little
foggy, so let us shed some light. Notice that, although after the first $40$
steps the distance decreases considerably, it never reaches $0$. This means
that at each time step, the distribution in the evolution sequence of $c_{2}$
assigns positive probability to (at least) one measurable set of data states
in which $l_{3}$ is farther from $l_{goal}$ than in all the measurable sets of
data states in the support of the distribution reached at the same time by
$c_{1}$. This is perfectly in line with our observations on the differences
between the variances of the distributions in Figure 6 (cf. the bottom row
pertaining $l_{3}$).

A natural question is then why the distance
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c_{2},c_{1})$, given
in Figure 9(b), is not equal to the constant $0$. This is due to the
combination of the Wasserstein distance with the hemimetric
$m^{\mathcal{D}_{3T}}_{\rho^{l_{3}}}$. Since
$m^{\mathcal{D}_{3T}}_{\rho^{l_{3}}}(\mathbf{d}_{1},\mathbf{d}_{2})=0$
whenever $\rho^{l_{3}}(\mathbf{d}_{2})<\rho^{l_{3}}(\mathbf{d}_{1})$ and the
Wasserstein lifting is defined as the infimum distance over the couplings, the
(measurable sets of) data states reached by $c_{2}$ that are better, with
respect to $\rho^{l_{3}}$, than those reached by $c_{1}$, are “hidden” in the
evaluation of
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c_{1},c_{2})$.
However, they do contribute to the evaluation of
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c_{2},c_{1})$.
Clearly, once we compute both
$\mathfrak{m}^{1}_{\rho^{l_{3}},\mathrm{OT}}(c_{1},c_{2})$ and
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c_{2},c_{1})$, we
can observe that the latter distance is less than $1/7$ of the former and we
can conclude thus that $\mathbf{s}_{1}$ shows a better behaviour than
$\mathbf{s}_{2}$.

As a final observation on Figure 9(b), we notice that, for instance, in the
time interval $I=[50,70]$ the pointwise distance between the evolution
sequences is $0$. This means that for each (measurable set) data state
$\mathbf{d}_{2}$ in the support of
$\mathcal{S}^{\mathcal{D}_{3T}}_{c_{2},\tau}$, $\tau\in I$, there is one
(measurable set) data state $\mathbf{d}_{1}$ in the support of
$\mathcal{S}^{\mathcal{D}_{3T}}_{c_{1},\tau}$, such that
$m^{\mathcal{D}_{3T}}_{\rho^{l_{3}}}(\mathbf{d}_{2},\mathbf{d}_{1})=0$.

(a) Evaluation of
$\operatorname{\mathbf{W}}(m^{\mathcal{D}_{3T}}_{\rho})(\mathcal{S}^{\mathcal{D}_{3T}}_{c_{1},\tau},\mathcal{S}^{\mathcal{D}_{3T}}_{c_{2},\tau})$
for $\rho\in\\{\rho^{l_{1}},\rho^{l_{2}},\rho^{l_{3}},\rho^{\max}\\}$.

(b) Evaluation of
$\operatorname{\mathbf{W}}(m^{\mathcal{D}_{3T}}_{\rho})(\mathcal{S}^{\mathcal{D}_{3T}}_{c_{2},\tau},\mathcal{S}^{\mathcal{D}_{3T}}_{c_{1},\tau})$
for $\rho\in\\{\rho^{l_{1}},\rho^{l_{2}},\rho^{l_{3}},\rho^{\max}\\}$.

Figure 10. Evaluation of the evolution metric with respect to different
targets.

Clearly, one can repeat a similar analysis for the other tasks of the system.
To give a general idea of the difference in the behaviour of the system in the
two scenarios, in Figure 10 we report the evaluation of the pointwise distance
between the evolution sequences of $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$
(Figure 10(a)), and viceversa (Figure 10(b)).

####  5.3.2. Different programs, same environment

We introduce two new programs: $\mathsf{P_{Tanks}^{+}}$ and
$\mathsf{P_{Tanks}^{-}}$. They are defined exactly as $\mathsf{P_{Tanks}}$,
the only difference being the value of $\delta_{l}$ used in the if/else guards
in $\mathsf{P_{in}}$ and $\mathsf{P_{out}}$. For $\mathsf{P_{Tanks}}$ we used
$\delta_{l}=0.5$; we use $\delta_{l}^{+}=0.7$ for $\mathsf{P_{Tanks}^{+}}$,
and $\delta_{l}^{-}=0.3$ for $\mathsf{P_{Tanks}^{-}}$. Our aim is to compare
the behaviour of $\mathbf{s}_{1}$ with that of systems $\mathbf{s}^{+}$ and
$\mathbf{s}^{-}$ having as initial configurations, respectively,
$c^{+}=\langle\mathsf{P_{Tanks}^{+}},\mathbf{d}_{0}\rangle_{\mathcal{E}_{1}}$
and
$c^{-}=\langle\mathsf{P_{Tanks}^{-}},\mathbf{d}_{0}\rangle_{\mathcal{E}_{1}}$.
Notice that the three systems are characterised by the same environment
evolution $\mathcal{E}_{1}$ and the same initial data state $\mathbf{d}_{0}$.
Moreover, as outlined above, we see the value of $\delta_{l}$ as an inherent
property of the program.

So, while the outcome of the experiment will be exactly the one the reader
expects, this example shows that with our framework we can compare the ability
of different programs to fulfil their objectives when operating on equal
footing, thus offering us the possibility of choosing the most suitable one
according to (possibly) external factors. For instance, we can imagine that
while $\delta_{l}^{-}$ allows for a finer control, it also causes more
frequent modifications to the flow rate $q_{1}$ than $\delta_{l}$ does. These
may in turn entail a higher risk of a breakdown of the pump attached to the
first tank. In this case, we can use the evolution metric to decide whether
the difference in the performance of the programs justifies the higher risk,
or if $\mathsf{P_{Tanks}}$ is still “good enough” for our purposes.

(a) Evolution metrics.

(b) Pointwise distances.

Figure 11. Evaluation of the distances between
$\mathsf{P_{Tanks}^{-}},\mathsf{P_{Tanks}},\mathsf{P_{Tanks}^{+}}$ with
respect to $\rho^{l_{3}}$.

In Figure 11 we present the results of the evaluations of
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c^{-},c_{1})$,
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c_{1},c^{+})$, and
$\mathfrak{m}^{\boldsymbol{1}}_{\rho^{l_{3}},\mathrm{OT}}(c^{-},c^{+})$, where
the target of the programs is $l_{3}=l_{goal}$ (i.e., the penalty function
$\rho^{l_{3}}$ is considered). In detail, in Figure 11(a) we report the values
of the three evolution metrics (see Figure 9 for an explanation of how to read
the graph), while in Figure 11(b) we report the pointwise (with respect to
time) Wasserstein distances between the distributions in the evolution
sequences of the systems.

##  6\. Robustness of programs

We devote this section to the analysis of the _robustness_ of programs with
respect to a data state and an environment. We also introduce two sensitivity
properties of programs, which we call _adaptability_ and _reliability_ ,
entailing the ability of the program to induce a similar behaviour in systems
that start operating from similar initial conditions.

###  6.1. Robustness

A system is said to be _robust_ if it is able to function correctly even in
the presence of uncertainties. In our setting, we formalise the _robustness_
of programs by _measuring_ their capability to _tolerate_ perturbations in the
environmental conditions.

First of all, we characterise the perturbations with respect to which we want
the program to be robust. Given a program $P$, a data state $\mathbf{d}$, and
an environment $\mathcal{E}$, we consider all data states
$\mathbf{d}^{\prime}$ such that the behaviour of the original system $\langle
P,\mathbf{d}\rangle_{\mathcal{E}}$ is at distance at most $\eta_{1}$, for a
chosen $\eta_{1}\in[0,1)$, from the behaviour of $\langle
P,\mathbf{d}^{\prime}\rangle_{\mathcal{E}}$. The distance is evaluated as the
evolution metric with respect to a chosen penalty function $\rho$ and a time
interval $I$. The idea is that $\rho$ and $I$ can be used to characterise any
particular situation we might want to study. For instance, if we think of a
drone autonomously setting its trajectory, then $I$ can be the time window
within which a gust of wind occurs, and $\rho$ can capture the distance from
the intended trajectory.

We can then proceed to measure the robustness. Let $\rho^{\prime}$ be the
penalty function related to the (principal) target of the system, and let
$\tilde{\tau}$ be an observable time instant. We evaluate the evolution
metric, with respect to $\rho^{\prime}$, between the behaviour of $\langle
P,\mathbf{d}\rangle_{\mathcal{E}}$ shown after time $\tilde{\tau}$, and that
of $\langle P,\mathbf{d}^{\prime}\rangle_{\mathcal{E}}$ for all the
perturbations $\mathbf{d}^{\prime}$ obtained above. We say that the robustness
of $P$ with respect to $\mathbf{d}$ and $\mathcal{E}$ is $\eta_{2}\in[0,1)$,
if $\eta_{2}$ is an upper bound for all those distances.

Notice that we use two penalty functions, one to identify the perturbations
($\rho$), and one to measure the robustness ($\rho^{\prime}$), which can be
potentially related to different targets of the system. For instance, in the
aforementioned case of the drone, $\rho^{\prime}$ can be related to the
distance from the final location, or to battery consumption. This allows for
an analysis of programs robustness in a general setting, where the
perturbations and the (long term) system behaviour take into account different
sets of variables. Moreover, the time instant $\tilde{\tau}$ plays the role of
a _time bonus_ given to the program to counter the effects of perturbations:
differences detectable within time $\tilde{\tau}$ are not considered in the
evaluation of the robustness threshold.

Definition 6.1 formalises the intuitions given above.

[Robustness] Let $\rho,\rho^{\prime}$ be two penalty functions, $I$ a time
interval, $\tilde{\tau}\in\mathrm{OT}$, and $\eta_{1},\eta_{2}\in[0,1)$. We
say that $P$ is
$(\rho,\rho^{\prime},I,\tilde{\tau},\eta_{1},\eta_{2})$-_robust_ with respect
to the data state $\mathbf{d}$ and the environment evolution $\mathcal{E}$ if

$\forall\,\mathbf{d}^{\prime}\in\mathcal{D}\text{ with
}\mathfrak{m}^{\lambda}_{\rho,\mathrm{OT}\cap I}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}},\langle
P,\mathbf{d}^{\prime}\rangle_{\mathcal{E}})\leq\eta_{1}$

it holds that

$\mathfrak{m}^{\lambda}_{\rho^{\prime},\\{\tau\in\mathrm{OT}\mid\tau\geq\tilde{\tau}\\}}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}},\langle
P,\mathbf{d}^{\prime}\rangle_{\mathcal{E}})\leq\eta_{2}\enspace.$

We can use our algorithm to measure the robustness of a given program. Given a
configuration $\langle P,\mathbf{d}\rangle_{\mathcal{E}}$, a set $\mathrm{OT}$
of observation times, a given threshold $\eta_{1}\geq 0$, a penalty function
$\rho$, and a time interval $I$, we can sample $M$ variations
$\\{\mathbf{d}_{1},\ldots,\mathbf{d}_{M}\\}$ of $\mathbf{d}$, such that for
any $i\in\\{1,\dots,M\\}$, $\mathfrak{m}^{\lambda}_{\rho,\mathrm{OT}\cap
I}(\mathbf{d},\mathbf{d}_{i})\leq\eta_{1}$. Then, for each sampled data state
$\mathbf{d}_{i}$ we can estimate the distance between $c=\langle
P,\mathbf{d}\rangle_{\mathcal{E}}$ and $c_{i}=\langle
P,\mathbf{d}_{i}\rangle_{\mathcal{E}}$, with respect to $\rho^{\prime}$, after
time $\tilde{\tau}$, namely
$\mathfrak{m}^{\lambda}_{\rho^{\prime},\\{\tau\in\mathrm{OT}\mid\tau\geq\tilde{\tau}\\}}(c,c_{i})$
for any $\tilde{\tau}\in\mathrm{OT}$. Finally, for each
$\tilde{\tau}\in\mathrm{OT}$, we let

$\xi_{\tilde{\tau}}=\max_{i\in\\{1,\dots,M\\}}\\{\mathfrak{m}^{\lambda}_{\rho^{\prime},\\{\tau\in\mathrm{OT}\mid\tau\geq\tilde{\tau}\\}}(c,c_{i})\\}\enspace.$

For the chosen $\eta_{1},\rho,I,\rho^{\prime}$ and $\tilde{\tau}$, each
$\xi_{\tilde{\tau}}$ gives us a lower bound to $\eta_{2}$, and thus the
robustness of the program.

Figure 12. Robustness of $\mathsf{P_{Tanks}}$ with respect to $\mathbf{d}_{s}$
and $\mathcal{E}_{1}$, for $M=100$, $\eta_{1}=0.3$, and $I=[0,50]$.

We can use our definition of robustness to show that whenever the program is
able to (initially) keep the difference with respect to the target
$l_{3}=l_{goal}$ bounded, then also the distance with respect to the target
$l_{1}=l_{goal}=l_{2}$ will be bounded. Formally, we instantiate the
parameters of Definition 6.1 as follows: $\rho=\rho^{l_{3}}$,
$\rho^{\prime}=\frac{1}{2}\rho^{l_{1}}+\frac{1}{2}\rho^{l_{2}}$, $I=[0,50]$,
and $\eta_{1}=0.3$. Figure 12 reports the evaluation of $\xi_{50}$ with
respect to $M=100$ variations $\mathbf{d}^{\prime}$ satisfying
$\mathfrak{m}^{\boldsymbol{1}}_{\rho,\mathrm{OT}\cap
I}(\langle\mathsf{P_{Tanks}},\mathbf{d}_{s}\rangle_{\mathcal{E}_{1}},\langle\mathsf{P_{Tanks}},\mathbf{d}^{\prime}\rangle_{\mathcal{E}_{1}})\leq
0.3$.  $\LHD$

###  6.2. Adaptability and reliability

We define the sensitivity properties of _adaptability_ and _reliability_ as
particular instances of the notion of robustness. Briefly, adaptability
considers only one penalty function, and reduces the time interval $I$ to the
sole initial time step $0$. Then, reliability also fixes $\tilde{\tau}=0$. Let
us present them in detail.

The notion of adaptability imposes some constraints on the _long term_
behaviour of systems, disregarding their possible initial dissimilarities.
Given the thresholds $\eta_{1},\eta_{2}\in[0,1)$ and an observable time
$\tilde{\tau}$, we say that a program $P$ is adaptable with respect to a data
state $\mathbf{d}$ and an environment evolution $\mathcal{E}$ if whenever $P$
starts its computation from a data state $\mathbf{d}^{\prime}$ that differs
from $\mathbf{d}$ for at most $\eta_{1}$, then we are guaranteed that the
distance between the evolution sequences of the two systems after time
$\tilde{\tau}$ is bounded by $\eta_{2}$.

[Adaptability] Let $\rho$ be a penalty function, $\tilde{\tau}\in\mathrm{OT}$
and $\eta_{1},\eta_{2}\in[0,1)$. We say that $P$ is
$(\tilde{\tau},\eta_{1},\eta_{2})$-_adaptable_ with respect to the data state
$\mathbf{d}$ and the environment evolution $\mathcal{E}$ if
$\forall\,\mathbf{d}^{\prime}\in\mathcal{D}$ with
$m^{\mathcal{D}}_{\rho,0}(\mathbf{d},\mathbf{d}^{\prime})\leq\eta_{1}$ it
holds
$\mathfrak{m}^{\lambda}_{\rho,\\{\tau\in\mathrm{OT}\mid\tau\geq\tilde{\tau}\\}}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}},\langle
P,\mathbf{d}^{\prime}\rangle_{\mathcal{E}})\leq\eta_{2}.$

We remark that one can always consider the data state $\mathbf{d}$ as the
ideal model of the world used for the specification of $P$, and the data state
$\mathbf{d}^{\prime}$ as the real world in which $P$ has to execute. Hence,
the idea behind adaptability is that even if the initial behaviour of the two
systems is quite different, $P$ is able to reduce the gap between the real
evolution and the desired one within the time threshold $\tilde{\tau}$. Notice
that being $(\tilde{\tau},\eta_{1},\eta_{2})$-adaptable for
$\tilde{\tau}=\min\\{\tau\mid\tau\in\mathrm{OT}\\}$ is equivalent to being
$(\tau,\eta_{1},\eta_{2})$-adaptable for all $\tau\in\mathrm{OT}$.

The notion of reliability strengthens that of adaptability by bounding the
distance on the evolution sequences from the beginning. A program is reliable
if it guarantees that small variations in the initial conditions cause only
bounded variations in its evolution.

[Reliability] Let $\rho$ be a penalty function and
$\eta_{1},\eta_{2}\in[0,1)$. We say that $P$ is
$(\eta_{1},\eta_{2})$-_reliable_ with respect to the data state $\mathbf{d}$
and the environment evolution $\mathcal{E}$ if
$\forall\,\mathbf{d}^{\prime}\in\mathcal{D}$ with
$m^{\mathcal{D}}_{\rho,0}(\mathbf{d},\mathbf{d}^{\prime})\leq\eta_{1}$ it
holds $\mathfrak{m}^{\lambda}_{\rho,\mathrm{OT}}(\langle
P,\mathbf{d}\rangle_{\mathcal{E}},\langle
P,\mathbf{d}^{\prime}\rangle_{\mathcal{E}})\leq\eta_{2}$.

The same strategy depicted above for the evaluation of the robustness can be
applied also in the case of adaptability and reliability. Briefly, we simulate
the behaviour of $M$ variations satisfying the requirement on the initial
distance between data states
$m^{\mathcal{D}}_{\rho,0}(\mathbf{d},\mathbf{d}_{i})\leq\eta_{1}$, and use the
bounds $\xi_{\tilde{\tau}}$ to estimate the adaptability bound $\eta_{2}$.
Similarly, for $\tau_{\min}=\min_{\mathrm{OT}}\tau$, $\xi_{\tau_{\min}}$ gives
a lower bound for its reliability.

(a) Adaptability modulo $\rho^{l_{3}}$.

(b) Comparison of adaptability modulo different penalty functions.

Figure 13. Adaptability of $\mathsf{P_{Tanks}}$ with respect to
$\mathbf{d}_{s}$ and $\mathcal{E}_{1}$, for $M=100$ and $\eta_{1}=0.3$.

We provide an example of the evaluation of the adaptability of
$\mathsf{P_{Tanks}}$. We consider the environment evolution $\mathcal{E}_{1}$,
i.e. the one corresponding to scenario 1, and, as initial data state, the data
state $\mathbf{d}_{s}$ such that $\mathbf{d}_{s}(q_{i})=0$ and
$\mathbf{d}_{s}(l_{i})=5$, for $i=1,2,3$. The rationale to consider
$\mathbf{d}_{s}$ in place of $\mathbf{d}_{0}$ is that $\rho(\mathbf{d}_{0})=1$
and thus $m^{\mathcal{D}_{3T}}_{\rho}(\mathbf{d}_{0},\mathbf{d}^{\prime})=0$
for all data states $\mathbf{d}^{\prime}$, for any penalty function
$\rho\in\\{\rho^{l_{1}},\rho^{l_{2}},\rho^{l_{3}},\rho^{\max}\\}$. Conversely,
there are several data states $\mathbf{d}^{\prime}$ such that
$m^{\mathcal{D}_{3T}}_{\rho}(\mathbf{d}_{s},\mathbf{d}^{\prime})>0$, so that
we can take into account the behaviour of the system when starting from
“worse” initial conditions. Let us stress that since the penalty functions
$\rho$ are evaluated only on (some of) the values $l_{i}$, for $i=1,2,3$, the
data states $\mathbf{d}^{\prime}$ that are within distance $\eta_{1}$ from
$\mathbf{d}_{s}$ can have any value in $\mathcal{D}_{q}$ initially assigned to
$q_{j}$, $j=1,2,3$ (and any value in $\mathcal{D}_{l}$ initially assigned to
the $l_{j}$ that are possibly not taken into account in the evaluation of
$\rho$).

To obtain the results reported in Figure 13 we considered, for each
$\rho\in\\{\rho^{l_{1}},\rho^{l_{2}},\rho^{l_{3}},\rho^{\max}\\}$, $M=100$
variations of the data state $\mathbf{d}_{s}$ within the (closed) ball of
radius $\eta_{1}=0.3$ with respect to $m^{\mathcal{D}_{3T}}_{\rho}$, and
performed each simulation over $5000$ samples ($N=1000$, $\ell=5$), keeping
$k=150$ as time horizon. Figure 13(a) represents the adaptability of
$\mathsf{P_{Tanks}}$ with respect to $\mathbf{d}$ and $\mathcal{E}_{1}$ for
$\rho=\rho^{l_{3}}$. We can then infer111Naturally, one should always take
into consideration the approximation error related to the computation of the
Wasserstein distance [SFG+21, Corollary 3.5, Equation (3.10)] for an exact
evaluation. Yet, in our presentation we are more interested in showing the
method rather than making exact calculations. that $\mathsf{P_{Tanks}}$ is
$(30,0.3,0.2)$-adaptable and $(50,0.3,0.05)$-adaptable with respect to
$\mathbf{d}_{s}$ and $\mathcal{E}_{1}$, when considering the target
$l_{3}=l_{goal}$. Figure 13(b) offers a comparison of the adaptability modulo
the various penalty functions defined for the three-tanks system.  $\LHD$

##  7\. Case-study: engine system

In this section we analyse a case study proposed in [LMMT21] and consisting in
two supervised self-coordinating refrigerated engine systems that are subject
to cyber-physical attacks aiming to inflict _overstress of equipment_
[GGI+15].

Since the impact of the attacks can be viewed as a flaw on system behaviour,
in our context it will be quantified by employing suitable penalty functions.
In particular, the impact of attacks aiming to overstress the equipment can be
quantified by adopting two different approaches: by simply measuring the level
of equipment’s stress or by comparing such a value with the alarm signals
generated by the system when attempting to perform attack detection. In the
former case one focuses only on damages inflicted by attacks, in the latter
case one is mainly interested in false positives and false negatives raised by
the detection system. In the former case we will use penalty functions
quantifying the level of stress of system’s equipment, in the latter case
penalty functions will quantify false negatives and false positives that are
produced when attempting to detect attacks and to raise alarm events. In both
cases, our simulations will allow us to _estimate the distance between the
genuine system and the system under attack_. Essentially, these simulations
will allow us to estimate the impact of attacks, thus quantifying how worse
the system under attack behaves with respect to the genuine one.

###  7.1. Modeling of the engine system

We start with describing a single engine. Then, we will describe the system
obtained by combining two identical engines with a supervisor. Finally, we
will introduce cyber-physical attacks tampering with a single engine or with
the whole system. Interestingly, our choice of modelling the program and the
environment separately well supports the specifications of these attacks:
since the attacks consist in pieces of malicious code that are somehow
injected in the logic of the system, an attack specified by a process $A$ can
be easily injected in a genuine system $\langle
P,\mathbf{d}\rangle_{\mathcal{E}}$ by composing $P$ and $A$ in parallel, thus
giving $\langle P~{}|~{}A,\mathbf{d}\rangle_{\mathcal{E}}$.

The engine we consider is a CPS whose logic aims to perform three tasks: [(i)]

regulate the speed,

maintain the temperature within a specific range by means of a cooling system,
and

detect anomalies. The first two tasks are on charge of a controller, the last
one is took over by an intrusion detection system, henceforth IDS.

The data space consists of the following variables:

  1. (1)

$\mathit{temp}$, representing a sensor detecting the temperature of the
engine;

  2. (2)

$\mathit{cool}$, representing an actuator to turn $\mathsf{on}$/$\mathsf{off}$
the cooling system;

  3. (3)

$\mathit{speed}$, representing an actuator to regulate the engine speed, with
values in the set $\\{\mathsf{slow},\mathsf{half},\mathsf{full}\\}$, where
$\mathsf{half}$ is the normal setting;

  4. (4)

$\mathit{ch\\_speed}$, representing a channel used by the IDS to command the
controller to set the actuator $\mathit{speed}$ at $\mathsf{slow}$, when an
anomaly is detected, or at $\mathsf{half}$, otherwise;

  5. (5)

$\mathit{ch\\_in}$ and $\mathit{ch\\_out}$, representing two channels used for
communication between the two engines: when an IDS commands to the controller
of its own engine to operate at $\mathsf{slow}$ speed, in order to compensate
the consequent lack of performance it asks to the controller of the other
engine to operate at $\mathsf{full}$ speed;

  6. (6)

$\mathit{stress}$, denoting the level of stress of the engine, due to its
operating conditions;

  7. (7)

six variables $\mathit{p_{k}}$, for $1\leq k\leq 6$, recording the
temperatures in the last six time instants;

  8. (8)

$\mathit{temp\\_fake}$, representing a noise introduced by malicious activity
of attackers tempering with sensor $\mathit{temp}$: we assume that the IDS can
access the genuine value $\mathit{temp}$, whereas the controller receives the
value of $\mathit{temp}$ through an insecure channel $\mathit{ch\\_temp}$ that
may be compromised and whose value is obtained by summing up the noise
$\mathit{temp\\_fake}$ and $\mathit{temp}$;

  9. (9)

$\mathit{ch\\_warning}$, representing a channel used by the IDS to raise
anomalies.

The environment evolution $\mathcal{E}$ affects these variables as follows:

  1. (1)

$\mathit{temp}$ changes by a value determined probabilistically according to a
uniform distribution ${\mathcal{U}}(I)$, where $I$ is the interval
$[-1.2,-0.8]$ if the cooling system is active ($\mathit{cool}=\mathsf{on}$),
whereas if the cooling system is inactive ($\mathit{cool}=\mathsf{off}$) then
$I$ depends on the speed: if $\mathit{speed}$ is $\mathsf{slow}$ then
$I=[0.1,0.3]$, if $\mathit{speed}$ is $\mathsf{half}$ then $I=[0.3,0.7]$, if
$\mathit{speed}$ is $\mathsf{full}$ then $I=[0.7,1.2]$;

  2. (2)

the variables $\mathit{p_{k}}$, for $k\in 1\dots 6$, are updated as expected
to record the last six temperatures detected by sensor $\mathit{temp}$;

  3. (3)

$\mathit{stress}$ remains unchanged if the temperature was below a constant
$\mathsf{max}$, set to 100 in our experiments, for at least $3$ of the last
$6$ time instants; otherwise, the variable is increased (reaching at most the
value $1$) by a constant $\mathsf{stressincr}$ that we set at $0.02$;

  4. (4)

all variables representing channels or actuators are not modified by
$\mathcal{E}$, since they are under the control of the program.

Summarising, the dynamics of $\mathit{p_{k}}$ and $\mathit{stress}$ can be
modelled via the following set of stochastic difference equations, with
sampling time interval $\Delta\tau=1$:

$\begin{array}[]{rcl}\mathit{p_{k}}(\tau+1)&=&\begin{cases}\begin{array}[]{cl}\mathit{temp}(\tau)&\text{
if }k=1\\\ \mathit{p_{k-1}}(\tau)&\text{ if
}k=2,\dots,6\end{array}\end{cases}\\\
\mathit{stress}(\tau+1)&=&\begin{cases}\begin{array}[]{cl}\mathit{stress}(\tau)+\mathsf{stressincr}&\text{
if }|\\{k\mid\mathit{p_{k}}(\tau)\geq\mathsf{max}\\}|>3\\\
\mathit{stress}(\tau)&\text{ otherwise}\end{array}\end{cases}\end{array}$ (15)

and $\mathit{temp}$ varies by a value that is uniformly distributed in an
interval depending on the state of actuators $\mathit{cool}$ and
$\mathit{speed}$:

$\begin{array}[]{rcl}\mathit{temp}(\tau+1)&=&\mathit{temp}(\tau)+v\\\
v&\sim&\begin{cases}\begin{array}[]{cl}{\mathcal{U}}[-1.2,-0.8]&\text{ if
}\mathit{cool}(\tau)=\mathsf{on}\\\ {\mathcal{U}}[0.1,0.3]&\text{ if
}\mathit{cool}(\tau)=\mathsf{off}\text{ and
}\mathit{speed}(\tau)=\mathsf{slow}\\\ {\mathcal{U}}[0.3,0.7]&\text{ if
}\mathit{cool}(\tau)=\mathsf{off}\text{ and
}\mathit{speed}(\tau)=\mathsf{half}\\\ {\mathcal{U}}[0.7,1.2]&\text{ if
}\mathit{cool}(\tau)=\mathsf{off}\text{ and
}\mathit{speed}(\tau)=\mathsf{full}.\end{array}\end{cases}\end{array}$ (16)

The whole engine, the controller and the IDS are modelled by the following
processes $\mathsf{Eng}$, $\mathsf{Ctrl}$ and $\mathsf{IDS}$, respectively:

$\begin{array}[]{lll}\mathsf{Eng}&\stackrel{{\scriptstyle\mathit{def}}}{{=}}&\mathsf{Ctrl}~{}|~{}\mathsf{IDS}\\\
\mathsf{Ctrl}&\stackrel{{\scriptstyle\mathit{def}}}{{=}}&\mathrm{if}\;[\mathit{ch\\_temp}<\mathsf{threshold}]\;\surd.\mathsf{Check}\;\mathrm{else}\;(\mathsf{on}\rightarrow\mathit{cool}).\mathsf{Cooling}\\\
\mathsf{Cooling}&\stackrel{{\scriptstyle\mathit{def}}}{{=}}&\surd.\surd.\surd.\surd.\mathsf{Check}\\\
\mathsf{Check}&\stackrel{{\scriptstyle\mathit{def}}}{{=}}&\mathrm{if}\;[\mathit{ch\\_speed}=\mathsf{slow}]\;((\mathsf{slow}\rightarrow\mathit{speed}),(\mathsf{off}\rightarrow\mathit{cool})).\mathsf{Ctrl}\\\
&&\mathrm{else}\;((\mathit{ch\\_in}\rightarrow\mathit{speed}),(\mathsf{off}\rightarrow\mathit{cool})).\mathsf{Ctrl}\\\
\mathsf{IDS}&\stackrel{{\scriptstyle\mathit{def}}}{{=}}&\mathrm{if}\;[\mathit{temp}>\mathsf{max}\wedge\mathsf{cool}=\mathsf{off}]\\\
&&((\mathsf{hot}\rightarrow\mathit{ch\\_warning}),(\mathsf{low}\rightarrow\mathit{ch\\_speed}),(\mathsf{full}\rightarrow\mathit{ch\\_out})).\mathsf{IDS}\\\
&&\mathrm{else}\;((\mathsf{ok}\rightarrow\mathit{ch\\_warning}),(\mathsf{half}\rightarrow\mathit{ch\\_speed}),(\mathsf{half}\rightarrow\mathit{ch\\_out})).\mathsf{IDS}\end{array}$

where we write $((e_{1}\to x_{1}),\dots,(e_{n}\to x_{n}))$ to denote the
prefix $(e_{1}\dots e_{n}\to x_{1}\dots x_{n})$.