Over-Approximation of Fluid Models

TL;DR

This paper introduces a novel, efficient method for computing formal bounds on the reachable states of uncertain nonlinear fluid models, aiding in verification of biochemical and performance systems.

Contribution

It relates nonlinear fluid model reachability to Markov decision processes and offers optimal control-based solutions for tight over-approximations.

Findings

The method provides tight bounds on the reachable set.

It outperforms existing over-approximation approaches.

The technique is efficient for large-scale systems.

Abstract

Fluid models are a popular formalism in the quantitative modeling of biochemical systems and analytical performance models. The main idea is to approximate a large-scale Markov chain by a compact set of ordinary differential equations (ODEs). Even though it is often crucial for a fluid model under study to satisfy some given properties, a formal verification is usually challenging. This is because parameters are often not known precisely due to finite-precision measurements and stochastic noise. In this paper, we present a novel technique that allows one to efficiently compute formal bounds on the reachable set of time-varying nonlinear ODE systems that are subject to uncertainty. To this end, we a) relate the reachable set of a nonlinear fluid model to a family of inhomogeneous continuous time Markov decision processes and b) provide optimal and suboptimal solutions for the family by relying on optimal control theory. The proposed technique is efficient and can be expected to provide tight bounds. We demonstrate its potential by comparing it with a state-of-the-art over-approximation approach.

Tables

Table 1. TABLE I: Results obtained by applying CORA and the heuristic implementation of Algorithm  1 to the SIRS model ( V ). While CORA terminated for D ≤ 6 𝐷 6 D\leq 6 , no estimations could be obtained in the case of D ≥ 7 𝐷 7 D\geq 7 due to out-of-memory errors. The heuristic implementation of Algorithm  1 , instead, is slower than CORA but scales to larger systems and provides tight bounds.

			CORA		Algorithm 1 for $𝒯\left(0.04\right)$		Algorithm 1 for $𝒯\left(0.03\right)$
$D$	$\left|𝒮\right|$	$\left|𝒦\right|$	Run time	$\stackrel{\text{Bound on}}{{sup}_t{\parallel ℰ\left(t\right)\parallel }_{\infty }}$	Run time	$\stackrel{\text{Bound on}}{{sup}_t{\parallel ℰ\left(t\right)\parallel }_{\infty }}$	Run time	$\stackrel{\text{Bound on}}{{sup}_t{\parallel ℰ\left(t\right)\parallel }_{\infty }}$
1	3	3	0m	0.187	2m	0.158	3m	0.163
2	6	8	1m	0.238	4m	0.124	6m	0.130
3	9	15	1m	0.296	8m	0.109	11m	0.116
4	12	24	2m	0.377	13m	0.100	17m	0.109
5	15	35	4m	0.494	18m	0.096	24m	0.103
6	18	48	5m	0.694	19m	0.091	29m	0.101
7	21	63	—	—	30m	0.091	43m	0.094
8	24	80	—	—	40m	0.084	53m	0.096
9	27	99	—	—	50m	0.087	65m	0.096
10	30	120	—	—	61m	0.091	78m	0.095

Equations

$${\dot{V}}_S^{{\kappa }_{\beta }}$$

$${\dot{V}}_I^{{\kappa }_{\beta }}$$

$${\dot{V}}_R^{{\kappa }_{\beta }}$$

$${\dot{V}}_S^{u_{\beta }}$$

$${\dot{V}}_I^{u_{\beta }}$$

$${\dot{V}}_R^{u_{\beta }}$$

$${\mathcal{E}}_B(t)=\sup \{|V_B^{u_{\beta }}(t)-V_B^0(t)|\mid |u_{\beta }(\cdot )|\le {\delta }_{\beta }(\cdot )\}$$

$$V_B^0(t)-{\mathcal{E}}_B(t)\le V_B^{u_{\beta }}(t)\le V_B^0(t)+{\mathcal{E}}_B(t)$$

$$\mathcal{R}(t)\subseteq\prod_{B\in\{S,I,R\}}\big{[}V_{B}^{0}(t)-\mathcal{E}_{B}(t);V_{B}^{0}(t)+\mathcal{E}_{B}(t)\big{]},$$

$$S+I$$

$$\bullet \text{ from state S into state I is VIuβ(t);}$$

$$\bullet \text{ from state I into state R is κ^β(t)+uβ(t);}$$

$$\bullet \text{ from state R into state S is 1.}$$

$${\dot{\pi }}_S^{u_{\beta }}$$

$${\dot{\pi }}_I^{u_{\beta }}$$

$${\dot{\pi }}_R^{u_{\beta }}$$

$${\dot{\pi }}_S^{u_{\beta },u_I}$$

$${\dot{\pi }}_I^{u_{\beta },u_I}$$

$${\dot{\pi }}_R^{u_{\beta },u_I}$$

$$\begin{gathered} ({\Phi }_B(\epsilon ))(t)=\sup \{|{\pi }_B^{u_{\beta },u_I}(t)-{\pi }_B^{0,0}(t)| \\ |u_{\beta }(\cdot )|\le {\delta }_{\beta }(\cdot )\text{ and }|u_I(\cdot )|\le {\epsilon }_I(\cdot )\}, \end{gathered}$$

$$\pi_{B}^{u^{\ast}_{\beta},u^{\ast}_{I}}(\hat{t})=\texttt{opt}\{\pi_{B}^{u_{\beta},u_{I}}(\hat{t})\mid\pi^{u_{\beta},u_{I}}\text{ solves~{}(\ref{ex_eq_sir_atomic_uu}) and}\\ |u_{\beta}(\cdot)|\leq\delta_{\beta}(\cdot),|u_{I}(\cdot)|\leq\varepsilon_{I}(\cdot)\},$$

$$c_{j,l}=\sum _{A_l\to A_{l^{\prime }}\in R_j}1\text{and}{d}_{j,l^{\prime }}=\sum _{A_l\to A_{l^{\prime }}\in R_j}1.$$

$$c_{j,1}{A}_1+\dots +c_{j,|\mathcal{S}|}{A}_{|\mathcal{S}|}\text{xlongrightarrow}{\Theta }_j{d}_{j,1}{A}_1+\dots +d_{j,|\mathcal{S}|}{A}_{|\mathcal{S}|}$$

$$\begin{gathered} {\mathcal{U}}_{\mathcal{K}}^{\delta }:=\{u:[0;T]\to {\mathbb{R}}_{>0}^{\mathcal{K}}\mid |u_{\alpha }(\cdot )|\le {\delta }_{\alpha }(\cdot ) \\ \text{and }u\text{ is measurable}\} \end{gathered}$$

$${\dot{V}}_B^u(t)$$

$$\displaystyle:=\sum_{1\leq j\leq|\mathcal{F}|}(d_{j,B}-c_{j,B})\Theta_{j}\big{(}V^{u}(t),\hat{\kappa}(t)+u(t)\big{)}$$

$$R_1=\{S\to I,I\to I\},$$

$${\Theta }_1(V,\kappa )=V_S{V}_I,$$

$$S+I$$

$${\dot{W}}_B(t)$$

$$\displaystyle\Theta_{j}\big{(}W,\hat{\kappa}+u\big{)}$$

$$r_{B,C}(V,\kappa )=\sum _{1\le j\le |\mathcal{F}|\mid B\to C\in R_j}{\Theta }_j(V,\kappa )/V_B$$

$${\dot{\pi }}_B^u(t)$$

$$:=-\sum _{C:C\ne B}{r}_{B,C}(V^u(t),\hat{\kappa }(t)+u(t)){\pi }_B^u(t)$$

$$+\sum _{C:C\ne B}{r}_{C,B}(V^u(t),\hat{\kappa }(t)+u(t)){\pi }_C^u(t)$$

$$r_{S,I}(V^u(t),\hat{\kappa }(t)+u(t))$$

$$r_{I,R}(V^u(t),\hat{\kappa }(t)+u(t))$$

$$r_{R,S}(V^u(t),\hat{\kappa }(t)+u(t))$$

Full text

I Introduction

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.

II The Main Idea in a Nutshell

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

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.

III Technical Preliminaries

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.

Definition 1**.**

An agent network (AN) is a triple $(\mathcal{S},\mathcal{K},\mathcal{F})$ of a finite set of states $\mathcal{S}=\{A_{1},\ldots,A_{|\mathcal{S}|}\}$, parameters $\mathcal{K}$ and reaction rate functions $\mathcal{F}$. Each reaction rate function $\Theta_{j}:\mathbb{R}_{\geq 0}^{\mathcal{S}\cup\mathcal{K}}\rightarrow\mathbb{R}_{\geq 0}$

• •

describes the rate at which reaction $j$ occurs;

• •

takes concentration and parameter vectors $V\in\mathbb{R}_{>0}^{\mathcal{S}}$ and $\kappa\in\mathbb{R}_{>0}^{\mathcal{K}}$, respectively;

• •

is accompanied by a multiset $R_{j}$ of atomic transitions of the form $A_{l}\to A_{l^{\prime}}$, where $A_{l}\to A_{l^{\prime}}$ describes an agent in state $A_{l}$ interacting and changing state to $A_{l^{\prime}}$.

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

[TABLE]

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

[TABLE]

We next introduce the ODE semantics of an AN.

Definition 2**.**

For a given AN $(\mathcal{S},\mathcal{K},\mathcal{F})$, a continuous parameter function $\hat{\kappa}:[0;T]\to\mathbb{R}_{>0}^{\mathcal{K}}$ and a piecewise continuous function $\delta:[0;T]\to\mathbb{R}^{\mathcal{K}}_{>0}$ with $\delta_{\alpha}(\cdot)<\hat{\kappa}_{\alpha}(\cdot)$ with $\alpha\in\mathcal{K}$, let

[TABLE]

denote the set of admissible uncertainties. Then, the reachable set of $(\mathcal{S},\mathcal{K},\mathcal{F})$ with respect to $\mathcal{U}_{\mathcal{K}}^{\delta}$ is given by the solution set $\{V^{u}\mid u\in\mathcal{U}_{\mathcal{K}}^{\delta}\}$, where $V^{u}$ solves

[TABLE]

for all $B\in\mathcal{S}$. The reachable set at time $0\leq t\leq T$ is given by $\mathcal{R}(t)=\{V^{u}(t)\mid u\in\mathcal{U}_{\mathcal{K}}^{\delta}\}$.

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

Remark 1**.**

Throughout the paper, the SIRS model from Section II is used to explain definitions and statements. All example environments refer to it.

Example 1**.**

Consider the agent network $(\{S,I,R\},$ $\{\beta\},$ $\{\Theta_{1},\Theta_{2},\Theta_{3}\})$ given by

[TABLE]

where $V=(V_{S},V_{I},V_{R})$ and $\kappa=(\kappa_{\beta})$. Let the time-varying uncertain recovery rate parameter be given by $\kappa_{\beta}\equiv\hat{\kappa}_{\beta}+u_{\beta}$, where $\hat{\kappa}_{\beta}$ denotes the nominal trajectory and $u=(u_{\beta})\in\mathcal{U}^{\delta}_{\{\beta\}}$ is the uncertainty function for some positive $\delta=(\delta_{\beta})$ such that $\delta_{\beta}<\hat{\kappa}_{\beta}$. The AN induces the reactions

[TABLE]

while the ODE system (9) is given by (1).

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.

Proposition 1**.**

In the case $i)$ and $ii)$ hold true, (9) admits a unique solution $V^{u}$ on $[0;T]$ for any uncertainty function $u\in\mathcal{U}^{\delta}_{\mathcal{K}}$. Moreover, there exists an $\eta>0$ such that $V^{u}(t)\in\mathbb{R}^{\mathcal{S}}_{\geq\eta}$ for all $u\in\mathcal{U}^{\delta}_{\mathcal{K}}$ and $0\leq t\leq T$.

Proof.

Local existence and uniqueness are ensured by [31, Section 3.3.1]. Let us define $W(0):=V(0)$,

[TABLE]

for all $B\in\mathcal{S}$ and let $e(\alpha)\in\mathbb{R}^{\mathcal{K}}$ denote the vector with $e(\alpha)_{\alpha^{\prime}}=1$ if $\alpha=\alpha^{\prime}$ and $e(\alpha)_{\alpha^{\prime}}=0$ when $\alpha\neq\alpha^{\prime}$. With this, $ii)$ implies for all $1\leq j\leq|\mathcal{F}|$ and $u\in\mathcal{U}^{\delta}_{\mathcal{K}}$ that

[TABLE]

because the function $\Theta_{j}$ is nonnegative. Hence, $\dot{V}^{u}(t)=F(V^{u}(t),\hat{\kappa}(t)+u(t))\geq G(t,V^{u}(t))$, thus implying that $V^{u}\geq W$ for all $u\in\mathcal{U}^{\delta}_{\mathcal{K}}$. We next show that $W$ is positive on $[0;T]$. To this end, let us assume towards a contradiction that there is $0<\tau\leq T$ such that $W_{A}(\tau)=0$ for some $A\in\mathcal{S}$. Thanks to the continuity of $W$, we may assume without loss of generality that $W$ is positive on $[0;\tau)$. With $\mathcal{W}(s):=W(\tau-s)$, it holds that $\dot{\mathcal{W}}(s)=-G(\tau-s,\mathcal{W}(s))$. There exists a sufficiently small interval $[0;\tau^{\prime})$ on which Euler’s sequence given by $(\mathcal{W}^{l})_{l\geq 0}$, where $\mathcal{W}^{0}:=\mathcal{W}(0)$ and $\mathcal{W}^{l+1}:=\mathcal{W}^{l}-\Delta t\cdot G(\tau-l\Delta t,\mathcal{W}^{l})$, converges to a local solution of $\mathcal{W}$ [40]. By construction, the sequence has to converge to a positive function on $(0;\tau^{\prime})$ as $\Delta t\to 0$. However, thanks to $ii)$, $\mathcal{W}_{A}^{0}=0$ implies $\mathcal{W}_{A}^{k}=0$ for all $k\geq 0$ regardless how small $\Delta t>0$ is, thus yielding a contradiction. Moreover, since $W>0$ and $\sum_{B\in\mathcal{S}}F_{B}(V^{u}(t),\hat{\kappa}(t)+u(t))=0$ for all $t\geq 0$, we also infer the existence of $V^{u}$ on the whole $[0;T]$. ∎

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).

Definition 3**.**

For a given AN $(\mathcal{S},\mathcal{K},\mathcal{F})$, define

[TABLE]

for all $B,C\in\mathcal{S}$ with $B\neq C$, $V\in\mathbb{R}^{\mathcal{S}}_{>0}$ and $\kappa\in\mathbb{R}_{>0}^{\mathcal{K}}$. Then, the coupled CTMC $(X^{u}(t))_{t\geq 0}$ underlying $(\mathcal{S},\mathcal{K},\mathcal{F})$ and $u\in\mathcal{U}_{\mathcal{K}}^{\delta}$ has state space $\mathcal{S}$ and its transition rate from state $B$ into state $C$ at time $t$ is $r_{B,C}(V^{u}(t),\hat{\kappa}(t)+u(t))$. The coupled Kolmogorov equations of $(X^{u}(t))_{t\geq 0}$ are

[TABLE]

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

[TABLE]

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).

Proposition 2**.**

For any uncertainty $u\in\mathcal{U}_{\mathcal{K}}^{\delta}$ and

[TABLE]

the solution of (11) exists on $[0;T]$ and satisfies $\pi^{u}(t)=V^{u}(t)$ for all $0\leq t\leq T$.

Proof.

Note that (11) rewrites into (9) if $\pi^{u}_{B}$ and $\dot{\pi}^{u}_{B}$ is replaced with $V^{u}_{B}$ and $\dot{V}^{u}_{B}$ for all $B\in\mathcal{S}$, respectively. With this, the claim follows from Proposition 1. ∎

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.

IV 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}$.

Definition 4**.**

For a given AN $(\mathcal{S},\mathcal{K},\mathcal{F})$ with uncertainty set $\mathcal{U}_{\mathcal{K}}^{\delta}$, the maximal deviation at time $t$ of (9) from $V^{0}$ is

[TABLE]

with $B\in\mathcal{S}$ and $\mathcal{E}=(\mathcal{E}_{B})_{B\in\mathcal{S}}$. With this, it holds that

[TABLE]

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}$.

Definition 5**.**

For $\varepsilon<V^{0}$ and $\mathfrak{u}=(u_{\mathcal{K}},u_{\mathcal{S}})\in\mathcal{U}_{\mathcal{K}}^{\delta}\times\mathcal{U}_{\mathcal{S}}^{\varepsilon}$, let $(\mathcal{D}^{\mathfrak{u}}(t))_{t\geq 0}$ be the decoupled CTMC with transition rates \big{(}r_{B,C}(V^{0}(t)+u_{\mathcal{S}}(t),\hat{\kappa}(t)+u_{\mathcal{K}}(t))\big{)}_{B,C} and the decoupled Kolmogorov equations

[TABLE]

where $f$ is as in Definition 3 and $\mathcal{U}_{\mathcal{S}}^{\varepsilon}$ is defined similarly to $\mathcal{U}_{\mathcal{K}}^{\delta}$ from Definition 2.

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

[TABLE]

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)$.

Remark 2**.**

Note that $V^{0}$ can be efficiently computed using a numerical ODE solver and by setting $u$ in (9) to zero.

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

Proposition 3**.**

Assume that $\mathcal{E}<V^{0}$. Then, for any $u_{\mathcal{K}}\in\mathcal{U}_{\mathcal{K}}^{\delta}$, there exists some $u_{\mathcal{S}}\in\mathcal{U}_{\mathcal{S}}^{\mathcal{E}}$ such that the solution of (14) subject to the initial condition $V(0)$ satisfies $\pi^{\mathfrak{u}}(t)=V^{u_{\mathcal{K}}}(t)$ for all $0\leq t\leq T$.

Proof.

For $\varepsilon$ with $\mathcal{E}\leq\varepsilon<V^{0}$, the definition of $\mathcal{E}$ implies that $u_{\mathcal{S}}:=V^{u_{\mathcal{K}}}-V^{0}\in\mathcal{U}_{\mathcal{S}}^{\varepsilon}$ for any $u_{\mathcal{K}}\in\mathcal{U}_{\mathcal{K}}^{\delta}$. Since $\pi^{u_{\mathcal{K}},u_{\mathcal{S}}}$ from (14) coincides with $\pi^{u_{\mathcal{K}}}$ from (11), Proposition 2 yields the claim. ∎

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}$.

Definition 6**.**

For a piecewise continuous function $\varepsilon<V^{0}$, let $\Phi(\varepsilon)=(\Phi_{B}(\varepsilon))_{B\in\mathcal{S}}$ be given by

[TABLE]

$(\Phi_{B}(\varepsilon))(t)$* denotes the maximal deviation of $\pi^{\mathfrak{u}}_{B}(t)$ from $\pi_{B}^{0}(t)$, where $\pi^{0}$ arises from $\pi^{\mathfrak{u}}$ in (14) if $\mathfrak{u}=0$.*

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$.

Theorem 1**.**

If $\Phi(\varepsilon)\leq\varepsilon$, then $\mathcal{E}\leq\varepsilon$.

Remark 3**.**

For the benefit of presentation, we prove Theorem 1 in Section IV-B by invoking the strict version of Pontryagin’s principle presented in Section IV-A.

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}$.

Theorem 2**.**

Fix some small $\varepsilon^{(0)}>0$ and set

[TABLE]

for all $k\geq 0$. If $\lim_{k\to\infty}\varepsilon^{(k)}=\varepsilon$ such that $\varepsilon\neq\infty$, then $\varepsilon$ is the smallest fixed point of $\Phi$ which satisfies $\varepsilon\geq\varepsilon^{(0)}$.

Proof.

Obviously, $\Phi$ is monotonic increasing, i.e., $\varepsilon\leq\varepsilon^{\prime}$ implies $\Phi(\varepsilon)\leq\Phi(\varepsilon^{\prime})$. With this, Kleene’s fixed point theorem yields the claim. ∎

Note that the computation of the sequence $(\varepsilon^{(k)})_{k}$ can be terminated if $\varepsilon^{(k+1)}<V^{0}$ is violated because in such case no bound can be obtained.

IV-A Optimal Solutions for inhomogeneous CTMDPs

In each step of the fixed point iteration from Theorem 2, a new value of $\Phi$ has to be computed. To this end, for any $0\leq\hat{t}\leq T$ and $A\in\mathcal{S}$, we have to

[TABLE]

While the solution of such optimization problems is particulary challenging in the case of nonlinear dynamics, time-varying systems such as (16) are easier to come by. This is because (16) is a linear system with additive and multiplicative uncertainties. More formally, (16) is linear in concentrations variables if the parameter variables are fixed and linear in parameter variables when the concentration variables are fixed.

Remark 4**.**

It is worth noting that (16) can be rewritten in the case of minimization (maximization is similar) to

[TABLE]

This defines a CTMDP with finite state space $\mathcal{S}$ and action space \big{(}\prod_{\alpha\in\mathcal{K}}[-\delta_{\alpha}(t);\delta_{\alpha}(t)]\big{)}\times\big{(}\prod_{A\in\mathcal{S}}[-\varepsilon_{A}(t);\varepsilon_{A}(t)]\big{)} at time $t$. The CTMDP is inhomogeneous due to the presence of the function $V^{0}$ in the transition rates from Definition 5.

For the benefit of presentation, we write in that what follows $\mathfrak{u}\in\mathcal{U}_{\mathcal{K}\cup\mathcal{S}}^{\mathfrak{b}}$ instead of $(u_{\mathcal{K}},u_{\mathcal{S}})\in\mathcal{U}_{\mathcal{K}}^{\delta}\times\mathcal{U}_{\mathcal{S}}^{\varepsilon}$, where $\mathfrak{b}_{\alpha}=\delta_{\alpha}$ and $\mathfrak{b}_{A}=\varepsilon_{A}$ for all $\alpha\in\mathcal{K}$ and $A\in\mathcal{S}$, respectively. Moreover, we recall that a solution of a differential inclusion $\dot{z}\in G(z)$ is any absolutely continuous function $z$ which satisfies $\dot{z}\in G(z)$ almost everywhere.

We solve (16) by modifying the strict version of Pontryagin’s principle [34] which is sufficient for optimality. Our modification of [34] is less general than the original because it is stated for CTMCs but it makes weaker assumptions (the concavity of $\hat{H}$ is required on positive values only).

Theorem 3**.**

For any $p\in\mathbb{R}^{\mathcal{S}}$, let $H(t,\pi,(u_{\mathcal{K}},u_{\mathcal{S}}),p)=\sum_{A\in\mathcal{S}}p_{A}h_{A}(t,\pi,(u_{\mathcal{K}},u_{\mathcal{S}}))$ and assume that, for any $0\leq t\leq\hat{t}$ and $p\in\mathbb{R}_{\geq 0}^{\mathcal{S}}$, the function

[TABLE]

is concave on $\mathbb{R}_{>0}^{\mathcal{S}}$. Then, any solution of the differential inclusion

[TABLE]

subject to (12) and $p(\hat{t})\equiv-\mathds{1}_{\{A=\cdot\}}(\cdot)$ ($p(\hat{t})\equiv\mathds{1}_{\{A=\cdot\}}(\cdot)$) minimizes (maximizes) the value of $\pi_{A}(\hat{t})$.

Proof.

The proof follows the argumentation of [34]. Fix some $u_{\mathcal{K}}\in\mathcal{U}_{\mathcal{K}}^{\delta}\cap\mathcal{C}([0;\hat{t}])$ and $u_{\mathcal{S}}\in\mathcal{U}_{\mathcal{S}}^{\varepsilon}\cap\mathcal{C}([0;\hat{t}])$ and let $\pi$ denote the solution underlying $\dot{\pi}(t)=h(t,\pi(t),(u_{\mathcal{K}}(t),u_{\mathcal{S}}(t)))$. Note that it suffices to consider continuous uncertainties because standard results from ODE theory and functional analysis ensure that the maximal value $\Phi(\varepsilon)$ can be attained by continuous uncertainties, that is

[TABLE]

where $\mathcal{C}^{\delta}_{\mathcal{K}}=\mathcal{U}^{\delta}_{\mathcal{K}}\cap\mathcal{C}([0;\hat{t}])$ and $\mathcal{C}^{\varepsilon}_{\mathcal{S}}=\mathcal{U}^{\varepsilon}_{\mathcal{S}}\cap\mathcal{C}([0;\hat{t}])$. For the ease of notation, let $\pi^{\ast}$, $p$ and $u^{\ast}$ denote a solution of the differential inclusion and set

[TABLE]

where $\cdot$ denotes the dot product. Thanks to the fact that $\dot{p}(t)=-p(t)\cdot(\partial_{\pi}h)(t,\pi^{\ast}(t),u^{\ast}(t))$, integration by parts yields

[TABLE]

With this, it holds that

[TABLE]

where the first inequality is implied by the concavity of $\pi\mapsto p\cdot h^{\ast}$, while the second inequality follows from the definition of $p\cdot h^{\prime}$ and the choice of $u^{\ast}$. In the case where we seek to maximize the value of $\pi_{A}(\hat{t})$, we note that $p_{\cdot}(\hat{t})=\mathds{1}_{\{A=\cdot\}}$ yields $0\geq\pi_{A}(\hat{t})-\pi^{\ast}_{A}(\hat{t})$. Since the case of minimization is similar, the proof is complete. ∎

We next identify structural conditions on $(\mathcal{D}^{\mathfrak{u}}(t))_{t\geq 0}$ which can be easily checked and that imply the technical requirement of concavity of Theorem 3.

• (A1)

For any $B,C\in\mathcal{S}$ and $0\leq t\leq T$, there exist Lipschitz continuous $k^{B\to C},k^{B\to C}_{i}\in[0;T]\to\mathbb{R}_{\geq 0}$ such that the transition rate function $r_{B,C}$ from Definition 3 satisfies

[TABLE]

for all $u_{\mathcal{K}}\in\mathbb{R}^{\mathcal{K}}$ and $u_{\mathcal{S}}\in\mathbb{R}^{\mathcal{S}}$.

• (A2)

For each $i\in\mathcal{K}\cup\mathcal{S}$, there exist unique $B_{i},C_{i}\in\mathcal{S}$ such that $k_{i}^{B\to C}\not\equiv 0$ implies $B=B_{i}$, $C=C_{i}$ and $k_{i}^{B\to C}>0$.

Assumption (A1) requires, essentially, the transition rate functions to be linear in the uncertainties, while (A2) forbids the same uncertainty to affect more than one transition of the decoupled CTMC $(\mathcal{D}^{\mathfrak{u}}(t))_{t\geq 0}$.

The next example demonstrates that our running example satisfies condition (A1) and (A2).

Example 2**.**

Recall that the transition rates of the decoupled CTMC of Example 1 are given by (15). Hence, $k^{S\to I}\equiv V^{0}_{I}$, $k^{I\to R}\equiv\hat{\kappa}_{\beta}$ and $k^{R\to S}\equiv 1$ and (A1) holds true. Condition (A2), instead, follows with $B_{I}=S$, $C_{I}=I$, $k_{I}^{S\to I}\equiv 1$ and $B_{\beta}=I$, $C_{\beta}=R$, $k_{\beta}^{I\to R}\equiv 1$.

The following crucial theorem can be shown in the presence of $\textbf{(A1)}-\textbf{(A2)}$. We wish to stress that the result can be also applied to an ICTMDP which is not induced by an AN.

Theorem 4**.**

Assume that $\textbf{(A1)}-\textbf{(A2)}$ hold true and fix some $A\in\mathcal{S}$. Then, the differential inclusion

[TABLE]

subject to $p_{B}(\hat{t})=-\mathds{1}_{\{A=B\}}$ ($p_{B}(\hat{t})=\mathds{1}_{\{A=B\}}$), with $B\in\mathcal{S}$, $0\leq t\leq\hat{t}$ and

[TABLE]

for $i\in\mathcal{K}\cup\mathcal{S}$, has a solution. Moreover, for any solution $p$ of (18), (4), the underlying solution $\pi$ that satisfies (12) and

[TABLE]

minimizes (maximizes) the value $\pi_{A}(\hat{t})$.

Proof.

In the following, we verify that a solution of (18)-(20) is a solution of the differential inclusion from Theorem 3. To this end, we first observe that

[TABLE]

Hence, we infer that

[TABLE]

This and Theorem 3 show that an optimal control $u^{\ast}$ must satisfy (4). Moreover, it implies that

[TABLE]

for all $0\leq\lambda\leq 1$ and $\pi,\pi^{\prime}\in\mathbb{R}_{>0}^{\mathcal{S}}$, thus yielding linearity (and thus also concavity) of $\hat{H}$ on $\mathbb{R}_{>0}^{\mathcal{S}}$. The last statement follows by noting that

[TABLE]

for all $E\in\mathcal{S}$. ∎

We wish to stress that Theorem 4 ensures that any solution $p$ of the differential inclusion (18), (4) induces an ODE solution $\pi$ of (20) such that $\pi_{A}(\hat{t})=\pi_{A}^{\ast}(\hat{t})$, where $\pi^{\ast}_{A}(\hat{t})$ denotes the solution of (16). This stands in stark contrast to the standard version of Pontryagin’s principle [31] which provides only necessary conditions for optimality, meaning that the value $\pi_{A}(\hat{t})$ arising from the standard version [31] may fail to satisfy $\pi_{A}(\hat{t})=\pi^{\ast}_{A}(\hat{t})$.

Solving a differential inclusion is a challenging task and requires one to assume in practice that it does not exhibit sliding or gazing modes [43, 44]. Fortunately, the next crucial results states that it is possible to obtain a specific solution of the differential inclusion (18), (4) by solving a Lipschitz continuous ODE system.

Theorem 5**.**

By replacing (4) with

[TABLE]

the differential inclusion (18) becomes an ODE system which is Lipschitz continuous in $t$ and $p$. With this change in place, the statement of Theorem 4 remains valid.

Proof.

Let $\mathcal{P}$ denote the drift of the ODE system (18) which underlies (21), that is

[TABLE]

with $u^{\ast}$ being as in (21). Fix some $(\tilde{t},\tilde{p})\in[0;\hat{t}]\times\mathbb{R}^{\mathcal{S}}$ and pick further two sequences $(t^{l},p^{l})_{l}$ and $(\tau^{l},\wp^{l})_{l}$ in $[0;\hat{t}]\times\mathbb{R}^{\mathcal{S}}$ which converge both to $(\tilde{t},\tilde{p})$ as $l\to\infty$. We first show that $\mathcal{P}_{B}(t^{l},p^{l})-\mathcal{P}_{B}(\tau^{l},\wp^{l})\to 0$ as $l\to\infty$. To this end, it suffices to observe that any $i\in\mathcal{K}\cup\mathcal{S}$ with $\psi_{i}(\tilde{t},\tilde{p})=0$ implies $\tilde{p}_{B_{i}}-\tilde{p}_{C_{i}}=0$ (recall that $k_{i}^{B\to C}\not\equiv 0$ yields $k_{i}^{B\to C}>0$). Hence, it holds that

[TABLE]

as $l\to\infty$. This shows the continuity of $\mathcal{P}$. To see also the Lipschitzianity, define

[TABLE]

for each $i\in\mathcal{K}\cup\mathcal{S}$ with $k_{i}^{B_{i}\to C_{i}}\not\equiv 0$. Note that $\psi_{i}(t,p)>0$ if and only if $p_{C_{i}}-p_{B_{i}}>0$ because $k_{i}^{B_{i}\to C_{i}}>0$ whenever $k_{i}^{B_{i}\to C_{i}}\not\equiv 0$. Moreover, for any $s\in\{-1,+1\}^{\mathcal{K}\cup\mathcal{S}}$, $\mathcal{P}$ is Lipschitz continuous on any bounded subset of $\bigcap_{i}G^{s_{i}}_{i}$ because $k_{i}^{B_{i}\to C_{i}}$ and $k^{B_{i}\to C_{i}}$ are Lipschitz continuous on $[0;\hat{t}]$. This shows that $\mathcal{P}$ is Lipschitz continuous on any bounded subset of $\bigcup_{s}\bigcap_{i}G^{s_{i}}_{i}$. With this, the continuity of $\mathcal{P}$ implies that $\mathcal{P}$ is Lipschitz continuous on any bounded subset of $[0;\hat{t}]\times\mathbb{R}^{\mathcal{S}}$. ∎

In the remainder of the paper, we replace (4) by (21). Theorem 5 ensures that (18) admits a unique solution $p$ and that the underlying optimal uncertainty $u^{\ast}(\cdot,p(\cdot))$ induces the minimal (maximal) value $\pi_{A}^{\ast}(\hat{t})$ via (21) and (20).

The next example demonstrates Theorem 4 and 5 in the context of the SIRS model from Example 1.

Example 3**.**

*We have seen in Example 2 that our running example satisfies the requirements of Theorem 4. In particular, if $\hat{\kappa}\equiv 1$, then (18) and (21) rewrite to *

[TABLE]

*and *

[TABLE]

respectively. The minimal value of, say, $\pi^{\mathfrak{u}}_{I}(\hat{t})$ can be obtained as follows. First, solve the ODE system (22) where the boundary condition is given by $p_{I}(\hat{t})=-1$ and $p_{S}(\hat{t})=p_{R}(\hat{t})=0$. Afterwards, using the obtained solution $p$, solve the ODE system (20) using the controls $u^{\ast}_{I}(\cdot,p(\cdot))$ and $u^{\ast}_{\beta}(\cdot,p(\cdot))$. A possible solution is visualized in Figure 1.

While Theorem 5 solves the problem from a theoretical point of view, it has to be noted that a numerical solution $\underline{p}$ of the Lipschitz continuous ODE system (4), (21) is an approximation of the true solution $p$. Hence, for any $\tilde{t}$ with $\psi_{i}(\tilde{t},p(\tilde{t}))\approx 0$, the computation of the optimal uncertainty $u_{i}(\tilde{t})$ may be hindered by the numerical errors underlying the ODE solver. The next crucial result addresses this issue by stating that, essentially, for each such $\tilde{t}$ the choice of $u_{i}(\tilde{t})$ is not important.

Theorem 6**.**

For any $\xi>0$, it is possible to efficiently compute some $\zeta>0$ such that the following holds. If $\mathfrak{u}\in\mathcal{U}_{\mathcal{K}\cup\mathcal{S}}^{\mathfrak{b}}$ is such that for all $i\in\mathcal{K}\cup\mathcal{S}$ it holds that $\mathfrak{u}_{i}(t)=u^{\ast}_{i}(t,p(t))$ whenever $|\psi_{i}(t,p(t))|\geq\zeta$, then $|\pi_{A}^{\mathfrak{u}}(\hat{t})-\pi_{A}^{\ast}(\hat{t})|\leq\xi$, where $\pi_{A}^{\mathfrak{u}}$ and $\pi_{A}^{\ast}$ is the solution of (14) and (16), respectively.

Proof.

With $c_{1}:=2\sup\{\mathfrak{b}_{i}(t)\mid i\in\mathcal{K}\cup\mathcal{S},0\leq t\leq T\}$, $c_{2}:=\lVert V(0)\rVert_{1}\cdot\max\{|k^{B_{i}\to C_{i}}_{i}(t)|\mid i\in\mathcal{K}\cup\mathcal{S},0\leq t\leq T\}$ and $\zeta=\xi/(T\cdot|\mathcal{K}\cup\mathcal{S}|c_{1}c_{2})$, let the function $u\in\mathcal{U}_{\mathcal{K}\cup\mathcal{S}}^{\mathfrak{b}}$ be such that $u_{i}(t)=u_{i}^{\ast}(t,p(t))$ whenever $|p_{C_{i}}(t)-p_{B_{i}}(t)|\geq\zeta$. Using the same notation as in the proof of Theorem 3, we infer in the case when $p(\hat{t})\equiv\mathds{1}_{\{A=\cdot\}}(\cdot)$ the following:

[TABLE]

where the first identity holds true because $\pi\mapsto p\cdot h^{\ast}$ is linear (see proof of Theorem 4), while the other identities follow as in the proof of Theorem 3. The above calculation yields

[TABLE]

where the second equality follows from the proof of Theorem 4. ∎

The above theorem states, essentially, that the choice of $u_{i}(t)$ is irrelevant at all time points $t$ with $|\psi_{i}(t,p(t))|<\zeta$. Hence, an uncertainty $\mathfrak{u}$ that is induced by a numerical solution of (18) can be used to solve (16).

We end the section by mentioning the following generalization of Theorem 4.

Remark 5**.**

The statement of Theorem 4 extends to the case where one seeks to minimize (maximize) the linear combination $\sum_{A\in\mathcal{S}}\sigma_{A}\pi_{A}(\hat{t})$, where $\sigma\in\mathbb{R}^{\mathcal{S}}$. The corresponding boundary condition of (18) is given by $p(\hat{t})=-\sigma$ ($p(\hat{t})=\sigma$).

IV-B Proof of Theorem 1

Armed with Theorem 4 and Theorem 5, we are now in a position to prove Theorem 1 under the assumption that the decoupled CTMC from Definition 5 satisfies (A1) and (A2). The assumption will be dropped in Section IV-C.

Proof of Theorem 1.

Let us assume towards a contradiction that $\delta$ and $\varepsilon$ are positive piecewise constant functions, that $\hat{\kappa}$ is analytic and that there exists an analytic uncertainty function $u_{\mathcal{K}}\in\mathcal{U}_{\mathcal{K}}^{\delta}$, a time $0<\hat{t}\leq T$ and some $A\in\mathcal{S}$ such that $|V^{u_{\mathcal{K}}}_{A}(\hat{t})-V^{0}_{A}(\hat{t})|=\varepsilon_{A}(\hat{t})$ and $|V^{u_{\mathcal{K}}}_{B}(t)-V^{0}_{B}(t)|<\varepsilon_{B}(t)$ for all $B\in\mathcal{S}$ and $0\leq t<\hat{t}$. Since $\varepsilon_{A}(\hat{t})>0$, we may assume without loss of generality that $V_{A}^{u_{\mathcal{K}}}(\hat{t})>V_{A}^{0}(\hat{t})$. With this, we consider the optimization problem

[TABLE]

where $h$ is as in (14). Let $\pi^{\ast}_{A}(\hat{t})$ denote the solution of (23) and set $\mathfrak{u}:=(u_{\mathcal{K}},u_{\mathcal{S}})$ with $u_{\mathcal{S}}:=V^{u_{\mathcal{K}}}-V^{0}$. Since $\varepsilon_{A}(\hat{t})=|\pi_{A}^{\mathfrak{u}}(\hat{t})-\pi_{A}^{0}(\hat{t})|\leq|\pi_{A}^{\ast}(\hat{t})-\pi_{A}^{0}(\hat{t})|=(\Phi_{A}(\varepsilon))(\hat{t})\leq\varepsilon_{A}(\hat{t})$, we infer $\pi^{\ast}_{A}(\hat{t})=\pi_{A}^{\mathfrak{u}}(\hat{t})$. Hence, $\mathfrak{u}$ is an optimal control and Theorem 4 and 5 imply that $u_{i}(t)=u_{i}^{\ast}(t)$ whenever $\psi_{i}(t,p(t))\neq 0$, where $u_{i}^{\ast}(t)$ is as in (21), $i\in\mathcal{K}\cup\mathcal{S}$ and $p$ solves (18) and (21). At the same time, the analyticity of $u_{\mathcal{K}}$, $\hat{\kappa}$ and $(\Theta_{j})_{j}$ implies that $u_{\mathcal{S}}=V^{u_{\mathcal{K}}}-V^{0}$ is analytic as well [41]. Since $\delta$ and $\varepsilon$ are piecewise constant and none of the $u_{i}$ can be locally constant (otherwise the $u_{i}$ in question would be constant on the whole $[0;T]$ by the identity theorem), we infer that $\psi_{i}(\cdot,p(\cdot))\equiv 0$ for all $i\in\mathcal{K}\cup\mathcal{S}$. Recall from the proof of Theorem 6 that

[TABLE]

and that the Hamiltonian $H(t^{\prime\prime},\pi^{\prime\prime},(u^{\prime\prime}_{\mathcal{K}},u^{\prime\prime}_{\mathcal{S}}),p^{\prime\prime})$ is invariant with respect to the value of $u^{\prime\prime}_{i}$ when $\psi_{i}(t^{\prime\prime},p^{\prime\prime})=0$. This, the above discussion and (24) imply that $\pi^{\ast}_{A}(\hat{t})=\pi_{A}^{(u^{\prime}_{\mathcal{K}},u^{\prime}_{\mathcal{S}})}(\hat{t})$ for any uncertainty $(u^{\prime}_{\mathcal{K}},u^{\prime}_{\mathcal{S}})$. As this contradicts $V_{A}^{u_{\mathcal{K}}}(\hat{t})>V_{A}^{0}(\hat{t})$, we infer the statement of the theorem in the case where $\delta$ and $\varepsilon$ are piecewise constant and $u_{\mathcal{K}}$ and $\hat{\kappa}$ are analytic. Thanks to the fact that analytic and piecewise constant functions are dense in set of bounded measurable functions on $[0;T]$, this suffices the claim. ∎

IV-C Sub-Optimal Solutions for inhomogeneous CTMDPs

It may happen that the decoupled CTMC from Definition 5 violates (A1) or (A2). We next discuss a procedure which allows one to transform a CTMC violating (A1)-(A2) into one which satisfies (A1)-(A2). We convey the main ideas using concrete examples.

The extension of Example 1 discussed next induces a decoupled CTMC which violates (A1).

Example 4**.**

Consider the agent network $(\{S,I,R\},$ $\{\alpha,\beta\},$ $\{\Theta_{1},\Theta_{2},\Theta_{3}\})$ given by

[TABLE]

where $V=(V_{S},V_{I},V_{R})$ and $\kappa=(\kappa_{\alpha},\kappa_{\beta})$. Let the time-varying uncertain infection and recovery parameter functions be given by $\kappa_{\alpha}\equiv\hat{\kappa}_{\alpha}+u_{\alpha}$ and $\kappa_{\beta}\equiv\hat{\kappa}_{\beta}+u_{\beta}$, respectively, where $u=(u_{\alpha},u_{\beta})\in\mathcal{U}^{\delta}_{\{\alpha,\beta\}}$ and $\delta=(\delta_{\alpha},\delta_{\beta})$. Then, the AN induces the reactions

[TABLE]

and the transition rates of the decoupled CTMC from Definition 5 are given by

[TABLE]

Since

[TABLE]

leads to the nonlinear term $u_{\alpha}u_{I}$, the decoupled CTMC does not satisfy (A1).

The idea is to substitute any nonlinear expression of uncertainties by a new uncertainty that bounds the original nonlinear expression. For instance, in the case of Example 4, we substitute $u_{\alpha}u_{I}$ with the new uncertainty $u_{\alpha|I}$ and set $\mathfrak{b}_{\alpha|I}:=\mathfrak{b}_{\alpha}\mathfrak{b}_{I}$ because $|u_{\alpha}(\cdot)u_{I}(\cdot)|\leq\mathfrak{b}_{\alpha}(\cdot)\mathfrak{b}_{I}(\cdot)$.

This motivates the following concept.

Definition 7**.**

For $\varepsilon<V^{0}$ a family of transition rates $(\hat{r}_{B,C})_{B,C}$ is an envelope of the transition rates $(r_{B,C})_{B,C}$ from Definition 5 if there exist Lipschitz continuous functions $k^{B\to C},k^{B\to C}_{i}\in[0;T]\to\mathbb{R}_{\geq 0}$, an index set $\mathcal{I}$ with $\mathcal{I}\cap(\mathcal{K}\cup\mathcal{S})=\emptyset$ and a piecewise continuous function $b:[0;T]\to\mathbb{R}^{\mathcal{I}}_{>0}$ such that for all $(u_{\mathcal{K}},u_{\mathcal{S}})\in\mathcal{U}_{\mathcal{K}}^{\delta}\times\mathcal{U}_{\mathcal{S}}^{\varepsilon}$ one can pick a $u_{\mathcal{I}}\in\mathcal{U}_{\mathcal{I}}^{b}$ so that

[TABLE]

for all $0\leq t\leq T$.

A possible envelope of the decoupled CTMC from Example 4 is given by $\hat{r}_{I,R}:=r_{I,R}$, $\hat{r}_{R,S}:=r_{R,S}$ and

[TABLE]

with $\mathcal{I}=\{\alpha|I\}$ and $\mathfrak{b}_{\alpha|I}:=\mathfrak{b}_{\alpha}\mathfrak{b}_{I}$.

By construction, any envelope satisfies (A1). It may however happen that an envelope does not satisfy (A2). To see this on an example, we extend Example 4 to the multi-class SIRS model [7] in which the overall population of agents is partitioned into classes, thus providing a better picture of the actual spread dynamics [38].

Example 5**.**

With $N\geq 2$ being the number of classes, the multi-class SIRS agent network is given by the atomic reactions

[TABLE]

where $1\leq\nu,\mu\leq N$. In the case where all rates are subject to uncertainty, the reactions are

[TABLE]

The first reaction expresses the fact that a susceptible agent of class $\nu$ may be infected by an infected agent from class $\mu$. The transition rates of the decoupled CTMC are

[TABLE]

The nonlinear terms $u_{\alpha_{\nu,\mu}}u_{I_{\mu}}$ prevent the decoupled CTMC to satisfy (A1). Similarly to Example 4, we thus consider the envelope

[TABLE]

with $\mathfrak{b}_{\alpha_{\nu,\mu}|I_{\mu}}:=\mathfrak{b}_{\alpha_{\nu,\mu}}\mathfrak{b}_{I_{\mu}}$. Unfortunately, envelope (5) violates (A2) because each $u_{I_{\mu}}$ is contained in $\hat{r}_{S_{1},I_{1}}$, …, $\hat{r}_{S_{N},I_{N}}$.

We continue by observing that envelope (5) can be transformed into a set of transition rates which satisfies (A1) and (A2). Indeed, if we substitute in each $\hat{r}_{S_{\nu},I_{\nu}}$ from (5) the uncertainty $u_{I_{\mu}}$ with $u_{I_{\nu,\mu}}$, the transition rates

[TABLE]

define a CTMC which satisfies (A1) and (A2). This is because every uncertainty function $u_{i}$, where $i\in\mathcal{K}\cup\mathcal{S}\cup\mathcal{I}$ and

[TABLE]

with $\mathfrak{b}_{I_{\nu,\mu}}:=\mathfrak{b}_{I_{\mu}}$ and $\mathfrak{b}_{\alpha_{\nu,\mu}|I_{\mu}}:=\mathfrak{b}_{\alpha_{\nu,\mu}}\mathfrak{b}_{I_{\mu}}$, appears in exactly one transition rate $\tilde{r}_{B,C}$.

This above discussion motivates the following.

Definition 8**.**

Assume that $(\hat{r}_{B,C})_{B,C}$ is an envelope of $(r_{B,C})_{B,C}$ given by

[TABLE]

for all $B,C\in\mathcal{S}$ and let $\mathcal{I}_{0}\dot{\cup}\mathcal{I}_{1}=\mathcal{K}\cup\mathcal{S}\cup\mathcal{I}$ be such that $(\hat{r}_{B,C})_{B,C}$ violates (A2) for each $i\in\mathcal{I}_{1}$. Then, the transition rate $\tilde{r}_{B,C}$ arises from $\hat{r}_{B,C}$ by substituting each occurrence of $u_{i}$ in $\hat{r}_{B,C}$ with $u_{i|B\to C}$, where $i\in\mathcal{I}_{1}$. By setting $\mathfrak{b}_{i|B\to C}:=\mathfrak{b}_{i}$, the coarsening of $(\hat{r}_{B,C})_{B,C}$ is given by $(\tilde{r}_{B,C})_{B,C}$.

Remark 6**.**

Note that (IV-C) is, up to a renaming of indices, a coarsening of the envelope (5). This can be seen by substituting each $u_{I_{\nu,\mu}}$ with $u_{I_{\mu}|S_{\nu}\to I_{\nu}}$.

The next result states that the coarsening of an envelope of the decoupled CTMC from Definition 5 allows one to estimate $\mathcal{E}$ from Definition 4.

Theorem 7**.**

Given the decoupled CTMC from Definition 5, let us assume that $(\hat{r}_{B,C})_{B,C}$ is an envelope for $(r_{B,C})_{B,C}$. Let further $(\tilde{r}_{B,C})_{B,C}$ denote the coarsening of $(\hat{r}_{B,C})_{B,C}$ as given in Definition 8 and let

[TABLE]

denote the Kolmogorov equations underlying the coarsening $(\tilde{r}_{B,C})_{B,C}$. Set further $\tilde{\pi}^{\mathfrak{u}}(0)=V(0)$ and

[TABLE]

Then, applying the fixed point iteration algorithm of Theorem 2 to $\Psi$ instead of $\Phi$ yields a bound on $\mathcal{E}$. Moreover, Theorem 4 and Theorem 5 carry over to the extended set of uncertainties $\mathcal{U}^{\delta}_{\mathcal{K}}\times\mathcal{U}^{\varepsilon}_{\mathcal{S}}\times\mathcal{U}^{b}_{\mathcal{I}}$ and can be used to compute $\Psi$.

Proof.

To see that $\Psi(\varepsilon)\leq\varepsilon$ implies $\mathcal{E}\leq\varepsilon$, let $\varepsilon$, $\delta$, $\hat{\kappa}$, $u_{\mathcal{K}}$, $u_{\mathcal{S}}$, $\mathfrak{u}$, $\hat{t}$ and $A$ be as in the proof from Section IV-B. Then, it holds that $\varepsilon_{A}(\hat{t})=|\pi_{A}^{\mathfrak{u}}(\hat{t})-\pi_{A}^{0}(\hat{t})|\leq|\tilde{\pi}_{A}^{\ast}(\hat{t})-\tilde{\pi}_{A}^{0}(\hat{t})|=(\Psi_{A}(\varepsilon))(\hat{t})\leq\varepsilon_{A}(\hat{t})$, where $\tilde{\pi}_{A}^{\ast}(\hat{t})$ solves

[TABLE]

and the inequality $|\pi_{A}^{\mathfrak{u}}(\hat{t})-\pi_{A}^{0}(\hat{t})|\leq|\tilde{\pi}_{A}^{\ast}(\hat{t})-\tilde{\pi}_{A}^{0}(\hat{t})|$ holds true because the definition of the envelope and the coarsening of an envelope ensure that for any function $(u_{\mathcal{K}},u_{\mathcal{S}})\in\mathcal{U}_{\mathcal{K}}^{\delta}\times\mathcal{U}_{\mathcal{S}}^{\varepsilon}$ there exists some function $u_{\mathcal{I}}\in\mathcal{U}_{\mathcal{I}}^{b}$ such that

[TABLE]

for all $0\leq t\leq T$. This implies $\tilde{\pi}_{A}^{\ast}(\hat{t})=\pi_{A}^{\mathfrak{u}}(\hat{t})$. Moreover, it can be observed that Theorem 4 and Theorem 5 hold for any CTMC which transition rates satisfy (A1) and (A2). Hence, they can be used to compute $\tilde{\pi}_{A}^{\ast}(\hat{t})$ (just replace the index set $\mathcal{K}\cup\mathcal{S}$ with $\mathcal{K}\cup\mathcal{S}\cup I$). This and the definition of the envelope and the coarsening of an envelope ensure the existence of some $u_{\mathcal{I}}\in\mathcal{U}_{\mathcal{I}}^{b}$ such that $\mathfrak{\tilde{u}}=(u_{\mathcal{K}},u_{\mathcal{S}},u_{\mathcal{I}})$ is optimal, i.e., $\tilde{\pi}_{A}^{\ast}(\hat{t})=\tilde{\pi}_{A}^{\mathfrak{\tilde{u}}}(\hat{t})$. In the case $b_{i}(\cdot)$ is piecewise constant for all $i\in\mathcal{I}$, the argumentation from Section IV-B leads to the desired contradiction. With this, the density argument from Section IV-B implies that $\mathcal{E}\leq\varepsilon$. Since $\Psi$ is monotonic increasing, the proof is complete. ∎

It is interesting to note that Theorem 1 allows one to derive bounds on $\mathcal{E}$ from Definition 4 using any estimation technique that applies to (16). Put different, Theorem 4 and 5 can be replaced, in principle, by any over-approximation technique applicable to time-varying linear systems with uncertain additive and multiplicative uncertainties. Note, however, that the bounds obtained by Theorem 4 and 5 cannot be improved if the decoupled CTMC satisfies (A1)-(A2) and can be expected to be tight even if (A1)-(A2) are violated because Theorem 7 relies on optimal control theory.

IV-D Algorithm

The previous sections gives rise to Algorithm 1 which summarizes all steps of our approximation technique. Apart from line 4 that has to compute an envelope of the transition rates of the decoupled CTMC from Definition 5 (see Definition 7), all steps of Algorithm 1 can be automatized. Indeed, the computation of the coarsening in line 8 is the variable substitution introduced in Definition 8, while $\Psi(\varepsilon^{\mathit{old}})$ in line 11 can be obtained by applying, for any $A\in\mathcal{S}$ and $0\leq\hat{t}\leq T$, Theorem 4 in order to compute the maximal and minimal value of $\tilde{\pi}^{\mathfrak{\tilde{u}}}_{A}(\hat{t})$ across all $\mathfrak{\tilde{u}}\in\mathcal{U}^{\delta}_{\mathcal{K}}\times\mathcal{U}^{\varepsilon}_{\mathcal{S}}\times\mathcal{U}^{b}_{\mathcal{I}}$.

Computation of an envelope. In the case the reaction rate functions of the agent network are multivariate polynomials as in the case of our running example discussed in Example 1-5, the envelope from line 4 can be efficiently computed by Algorithm 2. This makes our approach particularly suited to models from the field of biochemistry. Note that Algorithm 2 replaces any product of uncertainties by a new uncertainty and bounds it by the maximal value of the replaced product similarly to the discussion following Example 4.

Computation of $\Psi$. We conclude the section by discussing a rigorous and a heuristic implementation of line 11. In the case of the former, one has to combine Theorem 4, 5 and 6 with a verified numerical ODE solver as [45] that provides apart from a numerical solution also an estimation of the underlying numerical error [40]. Additionally to the numerical error, one has to account for the discretization error arising in the computation of $\Psi(\varepsilon)$, where the idea is to evaluate $(\Psi(\varepsilon))(\cdot)$ only at grid points $\hat{t}_{l}$ from $\mathcal{T}(\Delta t)=\{0,\Delta t,2\Delta t,\ldots,T\}$ by computing the maximal and minimal value of $\tilde{\pi}^{\mathfrak{\tilde{u}}}_{A}(\hat{t}_{l})$ from (7) for all $A\in\mathcal{S}$ and $\hat{t}_{l}\in\mathcal{T}(\Delta t)$.

The following result can be proven.

Theorem 8**.**

For any positive function $\varepsilon<V^{0}$ and $\xi\in(0;1)$, the function $t\mapsto(\Psi(\varepsilon))(t)$ can be approximated with precision $\xi$ by solving $4|\mathcal{S}|\Lambda T\xi^{-1}$ ODE systems of size $|\mathcal{S}|$, where

[TABLE]

and $\tilde{h}$ denotes the Kolmogorov equations underlying the transition rates $(\hat{r}_{B,C})_{B,C}$ from Algorithm 1.

Proof.

We prove that we need to solve $2|\mathcal{S}|\Lambda T\xi^{-1}$ instances of (18) and (20), respectively. To this end, we first note that $\tilde{\pi}^{\mathfrak{\tilde{u}}}$ from (7) is absolutely continuous and has derivative $\tilde{h}(\cdot,\tilde{\pi}^{\mathfrak{\tilde{u}}},\mathfrak{\tilde{u}})$ almost everywhere. This yields

[TABLE]

for any $0\leq t_{1}\leq t_{2}\leq T$. Since $\lVert\tilde{\pi}^{\mathfrak{\tilde{u}}}(t)\rVert_{\infty}\leq\lVert V(0)\rVert_{1}$ for all $0\leq t\leq T$, we infer that $\lVert\tilde{\pi}^{\mathfrak{\tilde{u}}}(t_{2})-\tilde{\pi}^{\mathfrak{\tilde{u}}}(t_{1})\rVert_{\infty}\leq\Lambda|t_{2}-t_{1}|$. This implies that we miss the actual value of $(\Psi_{\cdot}(\varepsilon))(\cdot)$ by at most $\Lambda\Delta t$ if we compute the maximal and minimal value of $\tilde{\pi}^{\mathfrak{\tilde{u}}}_{A}(\hat{t}_{l})$ for all $A\in\mathcal{S}$ and $\hat{t}_{l}\in\mathcal{T}(\Delta t)=\{0,\Delta t,2\Delta t,\ldots,T\}$. With this, we note that $\Lambda\Delta t\leq\xi$ implies $\Delta t\leq\xi/\Lambda$ which, in turn, induces $T/\Delta t=\Lambda T\xi^{-1}$ grid points. Since we need to compute the minimum and maximum value of $\tilde{\pi}^{\mathfrak{u}}_{A}(\hat{t}_{l})$ for all $A\in\mathcal{S}$ and $\hat{t}_{l}\in\mathcal{T}(\Delta t)$, this yields the claim. ∎

Apart from the rigorous implementation, our approach allows for a heuristic implementation. Here, instead of using a verified numerical ODE solver, one invokes a standard ODE solver in which the numerical error is minimized heuristically by varying the integration step size [40]. Similarly, one accounts heuristically for the discretization underlying $\mathcal{T}(\Delta t)$ by gradually refining an initially coarse discretization of $[0;T]$ until the approximations of $\Psi$ are reasonably close.

We wish to point out that both implementations naturally apply to parallelization because each single ODE system can be solved independently from the others.

V Numerical Evaluation

In this section we study the potential of our technique by applying it on the multi-class SIRS model from Section IV-C. To this end, we implemented an experimental prototype of the heuristic version from Section IV-D in Matlab by relying on the (non-verified) numerical ODE solver provided by the Matlab command ode45s. The heuristic implementation was compared with the state-of-the-art reachability analysis tool CORA [35] that covers nonlinear ODE systems with multiplicative uncertainty functions.

All experiments were performed on a 3.2 GHz Intel Core i5 machine with 8 GB of RAM. The Matlab solver ode45s was invoked with its default settings. For CORA, instead, the time step was set to $0.004$, while the expert settings were chosen as in the nonlinear tank example from the CORA manual [46]. The main findings are as follows.

• •

The bounds obtained by the heuristic implementation are tight;

• •

CORA is faster than the heuristic implementation in the case of smaller systems;

• •

The heuristic implementation scales to models that cannot be covered by CORA.

The above confirms the discussion from Section IV-D concerning the complexity of our approach. Indeed, for smaller ODE systems, our approach is inferior to CORA because of the discretization of the time interval $[0;T]$. However, for abstraction approaches such as CORA it is computationally prohibitive to obtain tight over-approximations for larger nonlinear systems in general. Instead, our approach requires to solve $4|\mathcal{S}|\Lambda T\xi^{-1}$ ODE systems of size $|\mathcal{S}|$, see Theorem 8. Moreover, at least as far as the heuristic implementation is considered, we were able to obtain tight bounds. In summary, we argue that our technique has the potential to complement state-of-the-art over-approximation approaches.

Multi-class SIRS Model. The global dynamics underlying (25) is given by

[TABLE]

where $u_{\alpha_{\nu,\mu}}$, $u_{\beta_{\nu}}$, $u_{\gamma_{\nu}}$ are time-varying uncertain functions and $1\leq\nu,\mu\leq D$ for some $D\geq 1$. The system has $3D$ ODE variables and $D^{2}+2D$ uncertainties. As discussed in Section IV-C, the transition rates (IV-C) define a function $\Psi$ as in Theorem 7 that can be used to estimate the function $\mathcal{E}$ underlying (V).

In our experiments, we randomly chose $\hat{\kappa}_{\alpha_{\nu,\mu}}\equiv 1.00$, $\hat{\kappa}_{\beta_{\nu}}\equiv 2.00$, $\hat{\kappa}_{\gamma_{\nu}}\equiv 3.00$ and $V^{u_{\mathcal{K}}}_{S_{\nu}}(0)=4.00+0.10(\nu-1)$, $V^{u_{\mathcal{K}}}_{I_{\nu}}(0)=V^{u_{\mathcal{K}}}_{R_{\nu}}(0)=1.00$ for all $1\leq\nu,\mu\leq D$. The time horizon was set to $T=3.00$, while all parameters were subject to uncertainties with modulus not higher than $\zeta=0.03$, i.e., $\mathfrak{b}_{\theta}(t)=\zeta$ for all $\theta\in\mathcal{K}$ and $0\leq t\leq T$.

Table I and Figure 2 summarize our findings. With increasing $D$, the tightness of the bounds provided by CORA decreases while the corresponding running times increase. In principle, the tightness can be improved by using stricter parameters (e.g., by decreasing the step size). This, however, increases the time and space requirements. Likewise, the over-approximation of larger models requires more resources in general. On our machine, for instance, $D\geq 7$ or time steps below $0.004$ led to out-of-memory errors. The heuristic implementation of Algorithm 1, instead, scales to larger instances of the running example and provides tight bounds. Indeed, since the $\Lambda$ from (29) can be chosen as $6D$ in the case of the ODE system (25), Theorem 8 implies that one has to solve $4|\mathcal{S}|\Lambda T\xi^{-1}=216D^{2}\xi^{-1}$ ODE systems of size $3D$ in order to guarantee that a numerical approximation of $\Psi(\varepsilon)$ from Theorem 7 misses the actual value of $\Psi(\varepsilon)$ by at most $\xi>0$. We approximated the values of $\Psi$ from Algorithm 1 using discretizations $\mathcal{T}(0.04)$ and $\mathcal{T}(0.03)$, where

[TABLE]

The run times account for the computation of the sequence $(\varepsilon^{(k)})_{k}$ from Algorithm 1 with $\eta=10^{-4}$. In agreement with Theorem 8, the running times exhibit a polynomial growth. Moreover, discretizations $\mathcal{T}(0.04)$ and $\mathcal{T}(0.03)$ induce bounds that are reasonably close.

VI Discussion

While Pontryagin’s principle and its extensions to systems with uncertain parameters have been used in the context of reachability analysis [7, 32], to the best of our knowledge, the principle has not been applied in the context of formal over-approximation of a general class of nonlinear ODE systems. This is because the principle is in general only a necessary condition for optimality, while its strict versions [34, 47] require concavity or convexity which is rarely satisfied by nonlinear ODE models. Additionally, Pontryagin’s principle induces in general a differential inclusion which can only be solved under additional assumptions [43, 44]. The present work addresses those problems by a) approximating the original nonlinear ODE system by a family of linear Kolmogorov equations (14) with multiplicative and additive uncertainties and; by b) showing that each family member can be over-approximated tightly and efficiently using a modified version of the strict version of Pontryagin’s principle [34].

The proposed approach is complementary to existing approximation techniques. Indeed, while it is less efficient than approaches that are based on monotonic systems and differential inequalities [48, 49, 22, 8], it may provide tighter bounds because it relies on optimal control theory. Instead, for approaches based on abstraction [21, 35] and the Hamilton-Jacobi equation [29, 30], in general it becomes computationally prohibitive to obtain tight over-approximations for larger nonlinear systems [19, 23]. Another point worth stressing is that many approaches applicable to nonlinear ODE models assume time-invariant uncertain parameters and uncertain initial conditions, while the present technique focusses on nonlinear ODE systems with time-varying uncertain parameters and fixed initial conditions. Since the proposed approximation technique relies on the availability of a concrete nominal solution, a direct extension to sets of initial conditions seems not to be possible. This notwithstanding we wish stress that it is particulary suited to systems biology where initial concentrations can be measured while reaction rates are often difficult to obtain and may vary with time.

VII Conclusion

In this work we presented an over-approximation technique for nonlinear ODE systems with time-varying uncertain parameters. Our approach provides verifiable bounds in terms of a family of linear Kolmogorov equations with uncertain additive and multiplicative time-varying parameters. To ensure efficient computation and tight estimations, we have established, to the best of our knowledge, a novel efficiently computable solution technique for a class of inhomogeneous continuous time Markov decision processes.

The presented over-approximation technique is efficient and can be expected to provide tight bounds because it relies on optimal control theory and allows for an algorithmic treatment in the case where the ODE system is given by multivariate polynomials. This makes it particularly suited to models from (bio)chemistry.

By comparing our approach with a state-of-the-art over-approximation technique in the context of the multi-class SIRS model from epidemiology [7], we have provided numerical evidence for the potential of our approach. The most pressing line of future work is the development of a tool which provides a rigorous implementation of the technique.

Acknowledgement

The author thanks Mirco Tribastone for helpful discussions. Parts of the work have been conducted when the author was with IMT Lucca. The author is supported by a Lise Meitner Fellowship that is funded by the Austrian Science Fund (FWF) under grant number M 2393-N32 (COCO).