Over-Approximation of Fluid Models

In the last decades, fluid (or mean-field) models underlying biochemical and computer systems have gained a lot of momentum. Possible examples are chemical reaction networks [1], optical switches [2] and layered queueing networks [3]. The main idea is to approximate the original stochastic model which is usually given in terms of a large-scale continuous time Markov chain (CTMC) by means of a compact system of ordinary differential equations (ODEs). When the number of agents (molecules, jobs, nodes etc.) present in the system tends to infinity, the simulation runs of a suitably scaled version of the CTMC can be shown to converge in probability to the deterministic solution of the underlying fluid model [4, 5]. The law of mass action from chemistry [6], for instance, has been shown to be the fluid model of a CTMC semantics stated on the molecule level [4].

Unfortunately, a precise parameterization of a fluid model is often not possible due to finite-precision measurements or stochastic noise [7, 8]. Hence, in order to verify that a nonlinear fluid model satisfies some given property in the presence of parameter functions that are subject to uncertainty, it becomes necessary to estimate the reachable set. This is because a finite set of possible ODE solutions (i.e., a proper subset of the reachable set) can only be used establish the presence of property violations but does not suffice to exclude their existence in general [9, 10]. Another reason is that closed-form expressions for reachable sets of nonlinear ODE systems are not known in general [11].

The estimation of reachable sets of continuous ODE systems is of crucial importance in the field of control engineering and has received a lot of attention over the decades. Linear ODE systems with uncertainties (alternatively, disturbances) are well-understood because in this case the reachable set can be shown to be convex. The situation where also matrix coefficients are uncertain [12], however, is more challenging than the standard control theoretical setting of additive uncertainties [13, 14, 15, 16]. Bounding the reachable set of nonlinear ODE systems is more difficult and there is a number of different techniques which complement each other. The abstraction approach approximates a nonlinear ODE system locally by an affine mapping [17, 18, 19] or a multivariate polynomial [20, 21]. The error can then be estimated using Taylor approximation and interval arithmetic [19, 22]. While abstraction techniques can cover many practical models, in general it is computationally prohibitive to obtain tight over-approximations for larger nonlinear systems [19, 23]. Lyapunov-like functions [9, 24, 25, 26, 27] known from the stability theory of ODE systems provide an alternative to abstraction techniques. Despite the fact that they often lead to tight bounds, their automatic computation is only possible in special cases [26]. In [28, 29, 30] it has been observed that over-approximation can be encoded as an optimal control problem. While theoretically appealing, the approach relies on the Hamilton-Jacobi equation, a partial differential equation which can only be solved for dynamical systems with few variables [31]. In a similar vein of research, [7, 32] used the necessary optimality conditions of Pontryagin’s principle [31] to derive heuristic estimations on reachable sets of nonlinear ODE systems.

Contributions. In the present paper, we introduce an over-approximation technique for the reachable sets of fluid models. The main idea is to exploit the fact that nonlinear fluid models can be related to the linear Kolmogorov equations of CTMCs. More specifically, the technical novelty of the present work is to prove that $i)$ a nonlinear fluid model can be over-approximated by solving a family inhomogeneous continuous time Markov decision processes (ICTMDPs) [33] with continuous action spaces and; $ii)$ to show that the family of ICTMDPs can be solved efficiently by modifying the strict version of Pontryagin’s principle [34] which is sufficient for optimality. This allows one to estimate the reachable set of an, in general, nonlinear ODE system by studying the reachable sets of a family of linear ODE systems.

For nonlinear fluid models, the proposed approach a) is efficient; b) induces bounds that can be expected to be tight and; c) allows for an algorithmic treatment in the case where the ODE system is given by multivariate polynomials, thus covering in particular biochemical models. A comparison with the state-of-the-art tool for reachability analysis CORA [35] in the context of the well-known SIRS model from epidemiology [7] confirms the potential of the proposed technique.

Related work on CTMDPs. With efficient solution techniques dating back to the sixties, CTMDPs [33] belong to one of the best studied classes of optimization problems. While there exists a large body of literature on homogeneous CTMDPs, however, much less is known about the inhomogeneous case. Moreover, most works on CTMDPs interpret controls as policies, meaning that only a subclass of uncertainties is admissible [36, Section 8.3]. Additionally, the cost function of interest is often either the discounted or the average cost [33] which cannot be used for the over-approximation of Kolmogorov equations. To the best of our knowledge, the only work which has studied the case of inhomogeneous CTMDPs featuring continuous action spaces and time dependent policies with respect to a cost which can be used for over-approximation is [37]. In this work, three concrete queueing systems were analyzed using Pontryagin’s principle [31]. For each of the three models, the underlying necessary conditions were shown to be already sufficient for optimality.

Paper outline. Section II provides a high-level discussion of our approach using a concrete example. Section III continues by introducing agent networks, a rich class of ODE systems that can be covered by our technique. In Section IV we first relate the reachable set of the original nonlinear ODE system underlying an agent network to the solution of a family of ICTMDPs. Afterwards, we present in Sections IV-A–IV-D an efficient solution approach to ICTMDPs. Section V compares a prototype implementation of the approach with CORA, while Section VI discusses how the approach complements existing approximation techniques. Section VII concludes the paper.

Notation. For nonempty sets $A$ and $\mathcal{I}$, let $A^{\mathcal{I}}$ denote the set of all functions from $\mathcal{I}$ to $A$. Note that elements of $A^{\mathcal{I}}$ can be interpreted as vectors with values in $A$ and coordinates in $\mathcal{I}$. We write $x\leq x^{\prime}$ for $x,x^{\prime}\in\mathbb{R}^{\mathcal{I}}$ whenever $x_{i}\leq x^{\prime}_{i}$ for all $i\in\mathcal{I}$. The equality of two functions $f$ and $g$, instead, is denoted by $f\equiv g$. By $\mathcal{S}$ we refer to the finite set of (agent) states; elements $V\in\mathbb{R}^{\mathcal{S}}$ of the reachable set of an ODE system are called concentrations instead. The derivative with respect to time of a function $V\in[0;T]\to\mathbb{R}^{\mathcal{S}}$ is denoted by $\dot{V}$. Instead, $\mathds{1}$ denotes the characteristic function.

We first discuss the problem and the proposed solution on the example of the SIRS model from epidemiology [38] that is given by the nonlinear ODE system

where $V^{\kappa_{\beta}}_{S},V^{\kappa_{\beta}}_{I}$ and $V^{\kappa_{\beta}}_{R}$ refers to the concentration of susceptible, infected and recovered agents, respectively, and $\kappa_{\beta}$ denotes the positive time-varying recovery rate parameter. We are interested in the case where the parameter function $\kappa_{\beta}$ is uncertain. More specifically, we assume that $\kappa_{\beta}\equiv\hat{\kappa}_{\beta}+u_{\beta}$, where $\hat{\kappa}_{\beta}$ is a known function resembling the nominal (or average) recovery parameter function, while $u_{\beta}$ is an unknown uncertainty which satisfies $|u_{\beta}(\cdot)|\leq\delta_{\beta}(\cdot)$ for some known function $\delta_{\beta}$. With this, the above ODE system rewrites to

The nominal solution $V^{0}$ corresponds to the case where $\kappa_{\beta}\equiv\hat{\kappa}_{\beta}$, i.e., when $u_{\beta}\equiv 0$. In practice, nominal parameter functions arise from finite-precision measurements, average behavior etc., while uncertainties account for the precision of measurements, conservative parameter estimations and stochastic noise.

Problem to solve. For a given time horizon $T>0$, we seek to provide, for each $0\leq t\leq T$, a superset which contains the reachable set $\mathcal{R}(t)=\{V^{u_{\beta}}(t)\mid|u_{\beta}(\cdot)|\leq\delta_{\beta}(\cdot)\}$. To this end, we bound the maximal deviation of $V^{u_{\beta}}$ from the nominal trajectory $V^{0}$, i.e., for each $B\in\{S,I,R\}$, we formally estimate the function

Since $V_{B}^{u_{\beta}}(t)=V_{B}^{0}(t)+V_{B}^{u_{\beta}}(t)-V_{B}^{0}(t)$, we infer

for all $0\leq t\leq T$ and $B\in\{S,I,R\}$. With this, it holds that

i.e., $\mathcal{R}(t)$ can be estimated by bounding the positive function $\mathcal{E}=(\mathcal{E}_{S}(\cdot),\mathcal{E}_{I}(\cdot),\mathcal{E}_{R}(\cdot))$. In what follows, we present a technique addressing this task.

First Step: Decoupling. Since the formal estimation of nonlinear dynamical systems is difficult, we relate the solution of (1) to that of a special linear ODE system. More specifically, we relate (1) to the linear Kolmogorov equations of a suitable CTMC. To this end, we first note that (1) is induced by the law of mass action [39] and the chemical reactions

The first reaction of (3) states that an infected agent can infect a susceptible one, while the second reaction implies that an infected agent eventually recovers. Instead, the third reaction expresses the fact that a recovered agent eventually loses its immunity and becomes susceptible again.

Apart from inducing the ODE system (1), the chemical reactions (3) induce also a probabilistic model. Intuitively, given a large group of agents interacting according to (3), the stochastic behavior of a single agent in the group is given in terms of a CTMC with the states $S,I,R$ such that at time $t$ the transition rate

The transition rate from state $S$ into state $I$ accounts for the fact that the probability of being infected is directly proportional to the concentration of infected agents.

The transition rates provided above imply that the transient probabilities of the CTMC satisfy the Kolmogorov equations

where $\pi^{u_{\beta}}_{S}(t)$, $\pi^{u_{\beta}}_{I}(t)$ and $\pi^{u_{\beta}}_{R}(t)$ denotes the probability that the fixed agent is susceptible, infected and recovered at time $t$, respectively.

We now make the pivotal observation that the solution $\pi^{u_{\beta}}$ of (5) with the initial condition given by $\pi^{u_{\beta}}(0)=V^{u_{\beta}}(0)$ is also a solution of (1), i.e., $\pi^{u_{\beta}}\equiv V^{u_{\beta}}$. (To see this, replace each $\pi^{u_{\beta}}_{A}$ with $V_{A}^{u_{\beta}}$ in (5).) Hence, if we are given $V^{u_{\beta}}_{I}$, the nonlinear ODE system (1) can be expressed in terms of the linear Kolmogorov equations (5).

Unfortunately, we cannot use (5) directly to estimate $\mathcal{E}$ from (2) because of the term $V^{u_{\beta}}_{I}$. We tackle this problem by replacing $V^{u_{\beta}}_{I}$ by $V^{0}_{I}+u_{I}$, where $V^{0}$ is the nominal trajectory of (1) in the case of $u_{\beta}\equiv 0$ and $u_{I}$ is a new uncertainty function that satisfies $|u_{I}(\cdot)|\leq\varepsilon_{I}(\cdot)$ for some positive function $\varepsilon_{I}$. This yields the linear ODE system

The key observation is that for any $u_{\beta}$, the uncertainty $u_{I}:=V^{u_{\beta}}_{I}-V^{0}_{I}$ induces a solution of (6) which coincides with the solution of (5), meaning that $\pi^{u_{\beta},u_{I}}\equiv V^{u_{\beta}}$ whenever $\pi^{u_{\beta},u_{I}}(0)=V^{u_{\beta}}(0)$.

Moreover, (6) is a linear ODE system that is decoupled from (1). Thus, instead of considering the maximal deviation of $V^{u_{\beta}}$ from $V^{0}$, $\mathcal{E}$, the above discussion motivates to focus on the maximal deviation of $\pi^{u_{\beta},u_{I}}$ from $\pi^{0,0}$, i.e.,

where $\varepsilon=(\varepsilon_{S}(\cdot),\varepsilon_{I}(\cdot),\varepsilon_{R}(\cdot))$ is a positive function, $B\in\{S,I,R\}$ and $\pi^{0,0}$ denotes the nominal solution of (6) when $u_{\beta}\equiv 0$ and $u_{I}\equiv 0$.

Intuitively, $\Phi$ takes a guess $\varepsilon=(\varepsilon_{S},\varepsilon_{I},\varepsilon_{R})$ for $\mathcal{E}$ as input and provides the new guess $\Phi(\varepsilon)=(\Phi_{S}(\varepsilon),\Phi_{I}(\varepsilon),\Phi_{R}(\varepsilon))$ for $\mathcal{E}$. Our goal is to find an $\varepsilon$ that satisfies $(\Phi_{B}(\varepsilon))(t)\leq\varepsilon_{B}(t)$ for all $B\in\{S,I,R\}$ and $0\leq t\leq T$ (or $\Phi(\varepsilon)\leq\varepsilon$ for short). This is because $\varepsilon_{I}$ is a guess for a bound on $|V_{I}^{u_{\beta}}-V_{I}^{0}|$, while $\Phi(\varepsilon)\leq\varepsilon$ implies that $|\pi_{I}^{u_{\beta},u_{I}}-\pi_{I}^{0,0}|$ is bounded by $\varepsilon_{I}$ whenever $|u_{I}|$ itself is bounded by $\varepsilon_{I}$. Building on this intuition, we will prove that $\Phi(\varepsilon)\leq\varepsilon$ implies $\mathcal{E}\leq\varepsilon$ and provide an algorithm that computes, whenever possible, the smallest such positive function $\varepsilon$.

Second Step: Approximation of Kolmogorov equations. The above discussion shows that an estimation of $\mathcal{E}$ requires one to evaluate the function $\Phi$. Thanks to the fact that (6) arises from (5) by decoupling, it can be seen that (6) describes the Kolmogorov equations of a CTMC with time-varying uncertain transition rates that are not coupled to the ODE system (1). This, in turn, allows one to compute any value of $\Phi$ by solving a family of tractable optimization problems. More specifically, the value $(\Phi_{B}(\varepsilon))(\hat{t})$ can be computed by determining two uncertainty functions $u^{\ast}_{\beta}$ and $u^{\ast}_{I}$ such that

with $\pi^{u_{\beta},u_{I}}(0)=V^{u_{\beta}}(0)$ and $\texttt{opt}\in\{\inf,\sup\}$. By interpreting the uncertainty functions $u^{\ast}_{\beta}$ and $u^{\ast}_{I}$ as optimal controls, (7) defines an optimal control problem with cost $\pi_{B}^{u_{\beta},u_{I}}(\hat{t})$.

A major result of the paper shows that uncertainty functions $u^{\ast}_{\beta}$ and $u^{\ast}_{I}$ can be efficiently computed by relying on a strict version of Pontryagin’s principle which is sufficient for optimality. Apart from solving (7) exactly, this allows one to devise an efficient procedure for the formal estimation of $\mathcal{E}$ whose bounds can be expected to be tight.

In this section, we introduce agent networks (ANs), a class of nonlinear ODE systems to which our over-approximation technique can be applied. ANs are, essentially, chemical reactions networks whose reaction rate functions are not restricted to the law of mass action. The distinctive feature of ANs is that their dynamics can be related to the linear Kolmogorov equations of CTMCs.

From a multiset $R_{j}$, we can extract two integer valued $|\mathcal{S}|$-vectors $d_{j}$ and $c_{j}$, counting how many agents in each state are transformed during a reaction (respectively produced and consumed). Specifically, for each $1\leq j\leq|\mathcal{F}|$, let $c_{jl},d_{jl}\in\mathbb{N}_{0}$ be such that

With these vectors, we can express the $j$-th reaction in the chemical reaction style [4] as follows:

We next introduce the ODE semantics of an AN.

The following example demonstrates Definition 1 and 2 in the context of the SIRS model from Section II. In particular, we remark the following.

In the following, we assume that an AN $(\mathcal{S},\mathcal{K},\mathcal{F})$ is accompanied by a finite time horizon $T>0$, a positive initial condition $V(0)\in\mathbb{R}^{\mathcal{S}}_{>0}$ and a Lipschitz continuous parameter function $\hat{\kappa}\in[0;T]\to\mathbb{R}_{>0}^{\mathcal{K}}$. Moreover, we require that each function $\Theta_{j}$

1. $i)$

   is analytic in $(V,\kappa)$ and linear in $\kappa$, i.e., it holds that $\Theta_{j}(V,c\kappa+c^{\prime}\kappa^{\prime})=c\Theta_{j}(V,\kappa)+c^{\prime}\Theta_{j}(V,\kappa^{\prime})$;

2. $ii)$

   satisfies $\Theta_{j}(V)=0$ whenever $V_{A_{l}}=0$ and $c_{j,l}>0$, where $c_{j,l}$ is as in (8).

Condition $i)$ enforces the existence of a unique solution (9) and allows us to apply Pontryagin’s principle in Section IV-A, while condition $ii)$ says essentially that the $j$-th reaction (8) can only take place when all its reactants have a positive concentration.

We wish to point out that $i)$ and $ii)$ can be easily checked because analytic functions are closed under summation, multiplication and composition. Additionally, functions $\Theta_{j}$ often enjoy a simple form in practical models (in the case of biochemistry, for instance, they are given in terms of monomials).

With $i)$ and $ii)$ in place, the following can be proven.

It can be seen that atomic transitions enforce conservation of mass, i.e., the creation and destruction of agents is ruled out at the first sight. This problem, however, can be alleviated by the introduction of artificial agent states, see [41].

Kolmogorov Equations of Agent Networks. Thanks to the fact that the dynamics of an AN arise from atomic transitions, it is possible to define a CTMC underlying a given AN which Kolmogorov equations are closely connected to the ODE system (9).

In the context of the SIRS example, Definition 3 gives rise to the transition rates (4), the uncertainty function $u=(u_{\beta})$ and the Kolmogorov equations (5). This is because the atomic transitions $S\to I$, $I\to R$ and $R\to S$ appear only in $R_{1}$, $R_{2}$ and $R_{3}$ of Example 1, respectively, thus yielding

where $\Theta_{1},\Theta_{2}$ and $\Theta_{3}$ are as in Example 1.

The next pivotal observation establishes a relation between the ODE system (9) and the Kolmogorov equations (11).

In the context of the SIRS example, Proposition 2 states that the solutions of (1) and (5) coincide whenever $\pi^{u}(0)=V(0)$.

It is possible to prove that (9) and (11) are the fluid limits of certain CTMC sequences in the case where $\pi^{u}(0)=V(0)/\lVert V(0)\rVert_{1}$ and the number of agents in the system tends to infinity, see [5, 41, 42] for details. We will not elaborate on this relation further because it is not required for the understanding of our over-approximation technique.

As anticipated in Section II and III, we estimate the reachable set of an AN with respect to an uncertainty set $\mathcal{U}_{\mathcal{K}}^{\delta}$, i.e., we bound $\mathcal{R}(t)=\{V^{u}(t)\mid u\in\mathcal{U}_{\mathcal{K}}^{\delta}\}$ for each $0\leq t\leq T$. To this end, we study the maximal deviation from the nominal trajectory $V^{0}$ attainable across $\mathcal{U}_{\mathcal{K}}^{\delta}$.

By Proposition 2, any trajectory $V^{u}$ of (9) coincides with the trajectory $\pi^{u}$ of (11) if $\pi^{u}(0)=V(0)$. Even though this allows one to relate the reachable set of a nonlinear system to that of a linear one, the transition rates of the coupled CTMC $(X^{u}(t))_{t\geq 0}$ depend on $V^{u}$. We address this by decoupling the transition rates of the coupled CTMC from $V^{u}$.

In the context of the AN from Example 1, the decoupled Kolmogorov equations (14) are given by (6) with $\mathfrak{u}\equiv(u_{\mathcal{K}},u_{\mathcal{S}})\equiv((u_{\beta}),(u_{I}))\in\mathcal{U}_{\mathcal{K}}^{\delta}\times\mathcal{U}_{\mathcal{S}}^{\varepsilon}$ $=\mathcal{U}_{\{\beta\}}^{\delta}\times\mathcal{U}_{\{S,I,R\}}^{\varepsilon}$. This is because the transition rates of the decoupled CTMC are

A direct comparison with the transition rates of the coupled CTMC given in (4) reveals that the original transition rate from $S$ into $I$, $V_{I}^{u_{\beta}}(t)$, is replaced with $V_{I}^{0}(t)+u_{I}(t)$.

The next result relates the original ODE system (9) to the decoupled Kolmogorov equations (14).

To provide an estimation of $\mathcal{E}$ using the Kolmogorov equations (14), we next define $\Phi(\varepsilon)$ as the maximal deviation from the nominal trajectory $\pi^{0}$ that can be attained across the uncertainties $u_{\mathcal{K}}\in\mathcal{U}_{\mathcal{K}}^{\delta}$ and $u_{\mathcal{S}}\in\mathcal{U}_{\mathcal{S}}^{\varepsilon}$.

As discussed in Section II, the goal is to find a positive function $\varepsilon$ such that $\Phi(\varepsilon)\leq\varepsilon$. This ensures that $|\pi^{\mathfrak{u}}-\pi^{0}|\leq\varepsilon$ for any $\mathfrak{u}=(u_{\mathcal{K}},u_{\mathcal{S}})\in\mathcal{U}_{\mathcal{K}}^{\delta}\times\mathcal{U}_{\mathcal{S}}^{\varepsilon}$ and implies, as stated in the next important result, that $\mathcal{E}\leq\varepsilon$.

A direct consequence of Theorem 1 is that a fixed point $\varepsilon^{\ast}$ of $\varepsilon\mapsto\Phi(\varepsilon)$ estimates $\mathcal{E}$ from above whenever $\varepsilon^{\ast}<V^{0}$.

The next result describes an algorithm for the computation of the least fixed point $\varepsilon^{\ast}$.