Automated Synthesis of Safe Digital Controllers for Sampled-Data Stochastic Nonlinear Systems

TL;DR

This paper introduces an automated method for synthesizing digital controllers that guarantee safety in stochastic nonlinear systems, combining rapid candidate generation with rigorous verification to ensure probabilistic safety thresholds.

Contribution

The paper presents a novel synthesis approach that integrates Monte Carlo-based candidate generation with SMT-based verification for safety guarantees in nonlinear stochastic systems.

Findings

Successfully applied to artificial pancreas, powertrain, and quadruple-tank systems.

Achieved high safety probability with rapid candidate synthesis.

Validated safety guarantees using SMT-based confidence intervals.

Abstract

We present a new method for the automated synthesis of digital controllers with formal safety guarantees for systems with nonlinear dynamics, noisy output measurements, and stochastic disturbances. Our method derives digital controllers such that the corresponding closed-loop system, modeled as a sampled-data stochastic control system, satisfies a safety specification with probability above a given threshold. The proposed synthesis method alternates between two steps: generation of a candidate controller pc, and verification of the candidate. pc is found by maximizing a Monte Carlo estimate of the safety probability, and by using a non-validated ODE solver for simulating the system. Such a candidate is therefore sub-optimal but can be generated very rapidly. To rule out unstable candidate controllers, we prove and utilize Lyapunov's indirect method for instability of sampled-data nonlinear systems. In the subsequent verification step, we use a validated solver based on SMT (Satisfiability Modulo Theories) to compute a numerically and statistically valid confidence interval for the safety probability of pc. If the probability so obtained is not above the threshold, we expand the search space for candidates by increasing the controller degree. We evaluate our technique on three case studies: an artificial pancreas model, a powertrain control model, and a quadruple-tank process.

Tables

Table 1. Table 1. Controller synthesis for the artificial pancreas system. l 𝑙 l – controller degree, K P subscript 𝐾 𝑃 K_{P} , K I subscript 𝐾 𝐼 K_{I} , K D subscript 𝐾 𝐷 K_{D} – controller gains, [ a ∗ , b ∗ ] superscript 𝑎 superscript 𝑏 [a^{*},b^{*}] – confidence interval for safety probability (with c = 0.99 𝑐 0.99 c=0.99 ), # o ​ ( # o 0 ) subscript # 𝑜 subscript superscript # 0 𝑜 \#_{o}(\#^{0}_{o}) – number of points used by the non-verified ODE solver at the end (beginning) of iteration, # c ​ ( # u ​ n ) subscript # 𝑐 subscript # 𝑢 𝑛 \#_{c}(\#_{un}) – number of candidates (unstable) sampled by the optimization algorithm, CPU(opt) – total (only optimize procedure) runtime in minutes.

$l$	$K_P$ $\times {10}^3$	$K_I$ $\times {10}^7$	$K_D$ $\times 10$	$[a^{\ast },b^{\ast }]$	${\#}_o\left({\#}_o^0\right)$	${\#}_c\left({\#}_{un}\right)$	CPU(opt)
0	-5.006	-	-	[0.938,0.988]	8(1)	165(0)	1130(586)
1	-5.4	-2.179	-	[0.939,0.989]	8(8)	217(0)	1534(873)
2	-5.716	-1.88	-2.002	[0.942,0.992]	8(8)	301(0)	2100(1312)

Table 2. Table 2. Controller synthesis for the fuel control system. See caption of Table 1 .

$l$	$K_P$	$K_I$	$K_D\times {10}^3$	$[a^{\ast },b^{\ast }]$	${\#}_o\left({\#}_o^0\right)$	${\#}_c\left({\#}_{un}\right)$	CPU(opt)
0	0.2713	-	-	[0.783,0.834]	128(64)	74(7)	1068(944)
1	0.2004	0.0537	-	[0.954,1]	128(128)	134(22)	1838(1690)
2	0.2082	0.0759	-4.9551	[0.963,1]	128(128)	214(72)	2337(2165)

Table 3. Table 3. Controller synthesis for the quadruple-tank system. K P 1 subscript 𝐾 subscript 𝑃 1 K_{P_{1}} , K I 1 subscript 𝐾 subscript 𝐼 1 K_{I_{1}} , K D 1 subscript 𝐾 subscript 𝐷 1 K_{D_{1}} , K P 2 subscript 𝐾 subscript 𝑃 2 K_{P_{2}} , K I 2 subscript 𝐾 subscript 𝐼 2 K_{I_{2}} , K D 2 subscript 𝐾 subscript 𝐷 2 K_{D_{2}} – controller gains, [ a ∗ , b ∗ ] superscript 𝑎 superscript 𝑏 [a^{*},b^{*}] – confidence interval for safety probability (with c = 0.99 𝑐 0.99 c=0.99 ), # c ​ ( # u ​ n ) subscript # 𝑐 subscript # 𝑢 𝑛 \#_{c}(\#_{un}) – number of candidates (unstable) sampled by the optimization algorithm, CPU(opt) – total (only optimize procedure) runtime in minutes, ∗ – the number of points in discretisation was increased from 1 to 2.

$K_{P_1}$	$K_{I_1}$	$K_{D_1}$	$K_{P_2}$	$K_{I_2}$	$K_{D_2}$	$[a^{\ast },b^{\ast }]$	${\#}_c\left({\#}_{un}\right)$	CPU(opt)
10.916	-	-	13.085	-	-	[0.85,0.90]	58(1)	164(30)
8.463	0.650	-	12.749	-	-	[0.89,0.94]	146(1)	252(73)
6.555	1.233	-	10.057	1.359	-	[0.94,0.99]	251(1)	370(113)
6.576	1.144	1.737	9.019	1.048	-	[0.93,0.98]	507(2)∗	717(246)
6.422	1.057	-0.075	6.760	1.724	3.973	[0.90,0.95]	654(2)∗	917(346)

Table 4. Table 4. Parameter values for the glucose-insulin regulatory model. w 𝑤 w (kg) is the body weight.

par	value	par	value	par	value
$w$	100	$k_e$	0.138	$k_{12}$	0.066
$k_{a_1}$	0.006	$k_{a_2}$	0.06	$k_{a_3}$	0.03
$k_{b_1}$	0.0034	$k_{b_2}$	0.056	$k_{b_3}$	0.024
$t_{maxI}$	55	$V_I$	$0.12\cdot w$	$V_G$	$0.16\cdot w$
$F_{01}$	$0.0097\cdot w$	$t_{maxG}$	40	$F_R$	0
$EG{P}_0$	$0.0161\cdot w$	$A_G$	0.8	$k_{\mathrm{int}}$	0.025

Equations

$$\frac{d}{dt}x(t)=$$

$$y(t_k)=$$

$$u(t)=$$

$$e(k)=r(k)-y(k),k\in {\mathbb{Z}}^{\ge 0}$$

$$u(k)=-\sum _{i=1}^L{a}_{i}u(k-i)+\sum _{i=0}^L{b}_{i}e(k-i),$$

$$x^c(k+1)$$

$$u(k)$$

$$\frac{d}{dt}{x}_1(t)$$

$$y_1(t_k)$$

$$x_1^c(k+1)$$

$$u_1(k\tau )$$

$$\frac{d}{dt}{x}_2(t)$$

$$x_2^c(k+1)$$

$$u_2(k\tau )$$

$$\underset{\Vert (x_1,u_1)\Vert \to 0}{\lim }\frac{\Vert g(x_1,u_1)\Vert }{\Vert (x_1,u_1)\Vert }=0\text{and}\underset{\Vert x\Vert \to 0}{\lim }\frac{\Vert l(x_1)\Vert }{\Vert x_1\Vert }=0.$$

$$\Vert g(x_1,u_1)\Vert \le \gamma \Vert x_1\Vert +\gamma \Vert u_1\Vert \text{and}\Vert l(x_1)\Vert \le \gamma \Vert x_1\Vert ,$$

$$\{\begin{array}{l}\frac{d}{dt}{x}_1(t)=A{x}_1(t)+B{C}_c{x}_1^c(k)-B{D}_{c}C{x}_1(k\tau )+g(x_1(t),u_1(k\tau ))\\ x_1^c(k+1)=G_c{x}_1^c(k)-H_{c}C{x}_1(k\tau )-H_{c}l(x_1(k\tau )),\end{array}$$

$$\Vert x_1(t)\Vert \le h_1(t-k\tau ,\gamma )\Vert x_1(k\tau )\Vert +h_2(t-k\tau ,\gamma )\Vert x_1^c(k)\Vert$$

$$x_1((k+1)\tau )$$

$$x_1^c(k+1)$$

$$\hat{g}(k\tau )={\int }_0^{\tau }{e}^{A(\tau -\lambda )}g(x_1(k\tau +\lambda ),u_1(k\tau ))d\lambda .$$

$$\Vert \hat{g}(k\tau )\Vert \le \gamma {\hat{h}}_1(\tau )\Vert x_1(k\tau )\Vert +\gamma {\hat{h}}_2(\tau )\Vert x_1^c(k)\Vert ,$$

$$x_2((k+1)\tau )$$

$$x_2^c(k+1)$$

$$\hat{G}:=\left[\begin{array}[]{cc}G-HD_{c}C&-HC_{c}\\ G_{c}&-H_{c}C\end{array}\right],$$

$$T\hat{G}T^{-1}=\left[\begin{array}[]{cc}G_{1}&0\\ 0&G_{2}\end{array}\right],$$

$$y(k):=\left[\begin{array}[]{c}x_{1}(k\tau)\\ x_{1}^{c}(k)\end{array}\right],\quad z=Ty=\left[\begin{array}[]{c}z_{1}\\ z_{2}\end{array}\right]=\left[\begin{array}[]{c}T_{1}y\\ T_{2}y\end{array}\right],$$

$$z(k+1)=T\hat{G}T^{-1}z(k)+T\left[\begin{array}[]{c}\hat{g}(k\tau)\\ \hat{l}(k\tau)\end{array}\right]=\left[\begin{array}[]{c}G_{1}z_{1}(k)+\hat{g}_{1}(k\tau)\\ G_{2}z_{2}(k)+\hat{g}_{2}(k\tau)\end{array}\right].$$

$$1+2\delta =\underset{i}{\min }|{\lambda }_i(G_1)|.$$

$$(1+\delta {)}^2{G_1^T}^{-1}{M}_1{G}_1^{-1}-M_1=-Q_1,G_2^T{M}_2{G}_2/(1+\delta {)}^2-M_2=-Q_2.$$

$$G_1^T{M}_1{G}_1-M_1={\bar{Q}}_1+({\delta }^2+2\delta )M_1,$$

$$G_2^T{M}_2{G}_2-M_2=-{\bar{Q}}_2+({\delta }^2+2\delta )M_2,$$

$$V(k):=z_1(k{)}^T{M}_1{z}_1(k)-z_2(k{)}^T{M}_2{z}_2(k),$$

$$V$$

$$-z_2(k+1{)}^T{M}_2{z}_2(k+1)+z_2(k{)}^T{M}_2{z}_2(k)$$

$$=z_1(k{)}^T(G_1^T{M}_1{G}_1-M_1)z_1(k)+2{\hat{g}}_1(k\tau {)}^T{M}_1(G_1{z}_1(k)+{\hat{g}}_1(k\tau ))$$

Full text

††copyright: none††conference: ; 2019††journalyear: ;

1. Introduction

Digital control (Ogata, 1995) is essential in many cyber-physical and embedded systems applications, ranging from aircraft autopilots to biomedical devices, due to its superior flexibility and scalability, and lower cost compared to its analog counterpart. The synthesis of analog controllers for linear systems is well-studied (Nise, 2016), but its extension to nonlinear and stochastic systems has proven much more challenging. Furthermore, digital control adds extra layers of complexity, e.g., time discretization and signal quantization. A common problem in both digital and analog control is the lack of automated synthesis techniques with provable guarantees, especially for properties beyond stability (e.g., safety) for nonlinear stochastic systems.

In this paper we address this problem by introducing a new method for the synthesis of probabilistically safe digital controllers for a large class of stochastic nonlinear systems, viz. sampled-data stochastic control systems. In such systems, the plant is a set of nonlinear differential equations subject to random disturbances, and the digital controller samples the noisy plant output, generating the control input with a fixed frequency.

Controllers are usually designed to achieve stability of the closed-loop system. Lyapunov’s indirect method provides conditions under which the stability of an equilibrium point of a nonlinear system follows from the stability of that point for the linearized version of the system (Khalil, 2002). Lyapunov’s method for sampled-data nonlinear systems is much more involved. Previous work (Nesic and Teel, 2004) provides sufficient conditions on the sampled-data linearized system that ensure stability of the sampled-data nonlinear system. Unfortunately, it is difficult to verify these conditions algorithmically. In this paper, we instead prove necessary conditions for stability that are easy to verify, and use them to restrict the controller synthesis domain. However, a stable system is not necessarily safe, as during the transient the system might reach an unsafe, catastrophic state. The synthesis approach that we propose overcomes this issue by deriving controllers that are safe.

Given an invariant $\phi$ (i.e., a correctness specification), and a nonlinear plant with stochastic disturbances and noisy outputs, our method synthesizes a digital controller such that the corresponding closed-loop system satisfies $\phi$ with probability above a given threshold $\vartheta$. The synthesis algorithm (Algorithm 1 in Section 5) is illustrated in Figure 1. It works by alternating between two steps: generation of a candidate controller $\mathbf{p}_{c}$, and verification of the candidate. $\mathbf{p}_{c}$ is generated via the optimize procedure (see Algorithm 2), which maximizes a Monte Carlo estimate of the satisfaction probability by simulating a discrete-time approximation of the system with a non-validated ODE solver. Such a candidate is, therefore, sub-optimal but very rapid to generate. To rule out unstable controller candidates, we prove and utilize Lyapunov’s indirect method for instability of sampled-data nonlinear systems. Along with $\mathbf{p}_{c}$, optimize returns an approximate confidence interval (CI) $[a,b]$ for the satisfaction probability.

Next, in the verification step (procedure verify), we use a validated solver based on SMT (Satisfiability Modulo Theories) to compute a numerically and statistically valid CI $[a^{\prime},b^{\prime}]$ for the satisfaction probability of $\mathbf{p}_{c}$. If the deviation between the approximate CI $[a,b]$ and the precise CI $[a^{\prime},b^{\prime}]$ is too large, indicating that the candidates generated by optimize are not sufficiently accurate, we increase the precision of the non-validated, fast solver (procedure update_discretization). If instead the precise probability is not above the threshold $\vartheta$, we expand the search space for candidates by increasing the controller degree.

Summarizing, the novel contributions of this paper are:

• •

we synthesize digital controllers for nonlinear systems subject to stochastic disturbances and measurement noise, while state-of-the-art approaches consider linear systems only;

• •

we prove Lyapunov’s indirect method for instability of nonlinear systems in closed-loop with digital controllers;

• •

we present a novel algorithm that synthesizes digital controllers with guaranteed probabilistic safety properties.

2. Sampled-data Stochastic Systems

We consider sampled-data stochastic control systems (SDSS), a rich class of control systems where the plant is specified as a nonlinear system subject to random disturbances. The controller periodically samples the plant output subject to random noise generating, using the plant output’s history, a control input that is kept constant during the sampling period with a zero-order hold; see Figure 2. The controller is characterized by a number of unknown parameters, which are the target of our synthesis algorithm.

Definition 2.1 (Sampled-data Stochastic Control System).

An SDSS can be described in the following state-space notation:

[TABLE]

where $x(\cdot)\in\mathbb{R}^{n}$ is the state of the plant; $x_{0}$ is the initial state at time $t=0$; $\mathbf{d}(\cdot)\in\mathbb{R}^{q}$ is the disturbance; $y(\cdot)$ is the plant output, which is a function of the state with additive i.i.d. noise $\eta\sim\mathcal{N}(0,W)$ with covariance matrix $W$; $u(\cdot)\in\mathbb{R}^{m}$ is the control input, updated at every sampling period $\tau>0$ by the digital controller $h$ (defined in Section 3); and $\mathbf{p}\in\mathbb{P}\subset\mathbb{R}^{2L+1}$ is the vector of unknown controller parameters, where $\mathbb{P}$ is a hyperbox (i.e., a product of closed intervals). The dynamics of the plant is governed by the vector field $f$, which is assumed to be in $C^{1}$, hence Lipschitz-continuous. We also assume that the output map $o(\cdot)$ is in $C^{1}$.

We assume that there are no time lags for transmitting the plant output to the controller and the control input to the plant. The disturbance $\mathbf{d}(\cdot)$ is a piecewise-continuous function having, for a time horizon $T$, a finite number of discontinuities. The discontinuity points and the value of $\mathbf{d}(\cdot)$ at each sub-domain can be defined in terms of a finite number of random parameters drawn from arbitrary distributions. Note that these assumptions on $\mathbf{d}(\cdot)$ are reasonably mild and allow us to define very general classes of systems, which subsume, for instance, numerical solutions of stochastic differential equations (Rüemelin, 1982).111Such numerical solutions rely on computing the value of the Wiener process at discrete time points, which makes it a special case of our disturbances.

3. Digital Controllers

The operation of a digital controller is succinctly indicated in Equation (1). These computations are generally performed using current and past output samples and past input samples.

Definition 3.1 (Digital Controller for SDSS).

Given an SDSS, we denote $y(k)=y(t_{k})$ and $u(k)=u(t_{k})$ and define the tracking error as

[TABLE]

where $r(\cdot)$ is the reference signal. The output of the controller is

[TABLE]

where $u(j)=e(j)=0$ for $j<0$222Note that if the controller has been previously deployed, i.e., it starts from a non-empty history, then $u(j),e(j)$ may be nonzero for $j<0$. and $L$ is the controller degree.

Controller design amounts to finding a degree $L$ and coefficients $\{a_{i}\}_{i=1}^{L},\{b_{i}\}_{i=0}^{L}$ that ensure the desired behavior of the closed-loop system. The vector of parameters $\mathbf{p}$ defined in (1) is recovered by setting $\mathbf{p}=[b_{0},a_{1},b_{1},\ldots,a_{L},b_{L}]^{T}$. An alternative description of the controller is via the state-space representation

[TABLE]

where $x^{c}(k)$ is the state of the controller and matrices $(G_{c},H_{c},C_{c},D_{c})$ need to be designed. The above two representations are equivalent. Given a controller in the form of (3), one can transform it to the representation (4), for instance by taking states as memories that store previous values of inputs/outputs. Given matrices $(G_{c},H_{c},C_{c},D_{c})$, one can easily compute coefficients $\{a_{i},b_{i}\}$ in (3) using matrix multiplications (Ogata, 1995).

3.1. Stability of the closed-loop system

A necessary requirement of any controller is stability of the closed-loop system. In our setting, we have the influence of both external inputs $(\mathbf{d},\eta,r)$ and initial states $(x_{0},x_{0}^{c})$. A suitable notion is input-to-state stability (ISS) (Sontag, 2008), which implies that bounded input signals must result in bounded outputs and, at the same time, that the effect of initial states must disappear as time goes to infinity. A necessary requirement of ISS is Lyapunov stability of the system ‘without’ external inputs, stated in the next definition, a requirement that can be applied to both continuous- and discrete-time systems.

Definition 3.2.

Consider a dynamical system with state space $\mathcal{D}\subset\mathbb{R}^{n}$ and without any external inputs, where $x_{e}\mathord{\in}\mathcal{D}$ is an equilibrium point and $\mathcal{D}$ is open. Then $x_{e}$ is called Lyapunov stable if for every $\epsilon>0$ there exists a $\delta>0$ such that for all $x(0)\in\mathcal{D}$ with $\|x(0)-x_{e}\|\leq\delta$, we have $\|x(t)-x_{0}\|\leq\epsilon$ for all $t\geq 0$.

It is very easy to verify stability for discrete-time linear systems.

Proposition 3.3 (Stability of Linear Systems (Khalil, 2002)).

A linear discrete-time system $x(k+1)=Gx(k)$ is Lyapunov stable at $x_{e}=0$ and $\lim_{k\rightarrow\infty}x(k)=0$ if and only if all eigenvalues of $G$ are inside unit circle. This condition is equivalent to the existence of positive definite matrices $M,Q$ such that $G^{T}MG-M=-Q$.

In the remainder of this section, we consider a version of SDSS in Definition 2.1 controlled by (4) without any external input, i.e., when $(\mathbf{d},\eta,r)$ are identically zero. We study Lyapunov stability of the closed-loop system without external inputs, which is necessary for having input-to-state stability. Let us put $(\mathbf{d},\eta,r)\equiv 0$ and define $x_{1}:=x-x_{e}$, $u_{1}:=u-u_{e}$, with $(x_{e},u_{e})$ being the equilibrium point of SDSS (1), i.e., $f(x_{e},u_{e},0)=0$. Similarly, define $y_{1}:=y-y_{e}$ with $y_{e}=o(x_{e})$ and $x_{1}^{c}:=x_{c}-x_{e}^{c}$ with $x_{e}^{c}$ being the equilibrium point for the controller. We then denote the plant dynamics after eliminating external inputs based on shifted version of variables by

[TABLE]

where $\bar{f}(x_{1},u_{1})=f(x_{1}+x_{e},u_{1}+u_{e},0)$ and $\bar{o}(x_{1})=o(x_{1}+x_{e})-y_{e}.$ Thus $x_{1}=0,u_{1}=0$ is an equilibrium point for (5). The digital controller dynamics is likewise given by

[TABLE]

where the minus sign is due to $r(\cdot)=0$ and negative feedback in (2). The nonlinear system (5) is controlled with the digital controller (6) by setting $u_{1}(t)=u_{1}(k\tau)$ for all $t\in[k\tau,(k+1)\tau)$, $k\in\mathbb{Z}^{\geq 0}$.

Ensuring stability of the sampled-data nonlinear control system (5)-(6) is difficult in general. Sufficient conditions for preserving stability under linearization are provided in (Nešić et al., 1999; Nesic and Teel, 2004), but they are hard to verify automatically. Rather, we provide an easy-to-check necessary condition for Lyapunov stability to reject unsuitable controllers. This necessary condition is based on Lyapunov’s indirect method developed here for sampled-data nonlinear systems. In the following we prove that if the linearized closed-loop system has an eigenvalue outside the unit circle, the nonlinear closed-loop system (5)-(6) is not Lyapunov stable, thus the system (1)-(4) is not input-to-state stable.

We first consider the linearized version of the closed-loop system (5)-(6), which is

[TABLE]

where $A:=\frac{\partial\bar{f}}{\partial x_{1}}(0,0)$, $B:=\frac{\partial\bar{f}}{\partial u_{1}}(0,0)$, and $C:=\frac{\partial\bar{o}}{\partial x_{1}}(0)$, Define $g(x_{1},u_{1}):=\bar{f}(x_{1},u_{1})-Ax_{1}-Bu_{1}$ and $l(x_{1}):=\bar{o}(x_{1})-Cx_{1}$ i.e., the non-linear terms describing the deviation between non-linear and linearized functions. Thus, we have

[TABLE]

Then, for any $\gamma>0$ there exists an $r(\gamma)>0$ such that

[TABLE]

for all $x_{1}\in\mathbb{R}^{n},u_{1}\in\mathbb{R}^{m}$ with $\|(x_{1},u_{1})\|\leq r(\gamma)$. We now simplify the dynamics of the closed-loop nonlinear system as

[TABLE]

The next lemma establishes a bound on $x_{1}(t)$ for any $t\in[k\tau,k\tau+\tau]$, as a function of $x_{1}(k\tau)$ and $x^{c}_{1}(k)$. Due to space constraints we present the proof of this lemma in the appendix.

Lemma 3.4.

Under dynamics (8), for a given $t\in[k\tau,(k+1)\tau]$ and any $\gamma>0$, we have

[TABLE]

if $\|(x_{1}(t_{1}),u_{1}(k\tau))\|\leq r(\gamma)$ for all $t_{1}\in[k\tau,t]$, with $r(\gamma)$ satisfying property (7). Functions $h_{1},h_{2}$ are continuous and nonnegative with $h_{1}(0,\gamma)=1$ and $h_{2}(0,\gamma)=0$.

The upper bound (9) enables us to study the effect of the nonlinear terms $g(\cdot)$ and $l(\cdot)$ in the sampled version of the dynamics, which can be written as

[TABLE]

with $G=e^{A\tau}$, $H=\int_{0}^{\tau}e^{A\lambda}Bd\lambda$, $\hat{l}(k\tau)=-H_{c}l(x_{1}(k\tau))$, and

[TABLE]

Next, we derive a bound for $\hat{g}(\cdot)$ in terms of $x_{1}$ and $x_{1}^{c}$.

Lemma 3.5.

For any $\gamma>0$, there exist continuous functions $\hat{h}_{1},\hat{h}_{2}$ such that the following inequality holds for $\hat{g}(\cdot)$ defined in (10),

[TABLE]

if $\|(x_{1}(t),u_{1}(k\tau))\|\leq r(\gamma)$ for all $t\in[k\tau,(k+1)\tau]$, with $r(\gamma)$ satisfying property (7). Functions $\hat{h}_{1},\hat{h}_{2}$ are nonnegative with $\hat{h}_{1}(0)=0$ and $\hat{h}_{2}(0)=0$.

The explicit form of $\hat{h}_{1},\hat{h}_{2}$ is provided, along with the proof, in the appendix. We are now ready to state our main result of this section.

Theorem 3.6.

The continuous-time nonlinear system (5) controlled with the digital controller (6) is unstable if the linearized continuous-time system controlled by the same digital controller has a pole outside the unit circle.

Sketch of the proof. We prove the theorem by contradiction. We show that there is an $r>0$ such that for any $\delta>0$ we can find an initial state for the system and the controller with $\|(x_{1}(0),x_{1}^{c}(0))\|\leq\delta$ and time $t>0$ with $\|(x_{1}(t),x_{1}^{c}(t))\|>r$. We construct a Lyapunov function for the linear system that is strictly increasing on a suitable set of initial states. By a proper selection of $r$, we show that this Lyapunov function is also strictly increasing on the nonlinear system if the trajectory remains inside the ball with radius $r$, which is a contradiction.

Proof.

The closed-loop linearized system will have the following dynamics in discrete time

[TABLE]

with $G=e^{A\tau}$ and $H=\int_{0}^{\tau}e^{A\lambda}Bd\lambda$. This gives the following state transition matrix

[TABLE]

which is assumed to have at least one eigenvalue outside the unit circle. We cluster the eigenvalues of $\hat{G}$ into a group of eigenvalues outside the unit circle and a group of eigenvalues on or inside the unit circle. Then there is a nonsingular matrix $T$ such that

[TABLE]

where $G_{1}^{-1}$ is stable. In other words, $G_{1}$ contains all of the eigenvalues of $\hat{G}$ from the first group. (The matrix $T$ can be found for instance by transforming $\hat{G}$ into its real Jordan form.) Let us define

[TABLE]

where the partitions of $z$ and $T$ are compatible with the dimensions of $G_{1}$ and $G_{2}$. The dynamics of $z$ becomes

[TABLE]

Now define $\delta>0$ by

[TABLE]

Then both $G_{2}/(1+\delta)$ and $(1+\delta)G_{1}^{-1}$ are stable matrices. According to Proposition 3.3, there are positive definite matrices $M_{1},M_{2}$ and $Q_{1},Q_{2}$ such that the following matrix equalities hold

[TABLE]

and we get

[TABLE]

with $\bar{Q}_{1}:=G_{1}^{T}Q_{1}G_{1}$ and $\bar{Q}_{2}:=(1+\delta)^{2}Q_{2}$. For the function

[TABLE]

we have

[TABLE]

where the positive constant

[TABLE]

is obtained using definition of $\hat{g}_{1},\hat{g}_{2}$ as a function of $\hat{g}$ and Lemma 3.5 as

[TABLE]

Similarly, we use the Cauchy-Schwarz inequality for

[TABLE]

as

[TABLE]

Then, function $V(\cdot)$ satisfies

[TABLE]

with $c_{0}=\min_{i,j}\{\lambda_{i}(\bar{Q}_{1}),\lambda_{j}(\bar{Q}_{2})\}$.

Take $0<\gamma_{0}\leq 1$ sufficiently small such that $c_{0}-c_{1}\gamma_{0}-c_{2}\gamma_{0}^{2}\geq 0$ with its associated radius $r_{0}=r(\gamma_{0})$. Note that this is always possible since $\hat{h}_{1},\hat{h}_{2}$, thus $c_{1},c_{2}$, are bounded on the interval $\gamma\in(0,1]$. Then we have $V(k+1)\geq(1+\delta)^{2}V(k)$ as long as $\|(x_{1}(t),u_{1}(k\tau))\|\leq r_{0}$. For the proof of instability, take radius

[TABLE]

and any initial condition $y(0)=(x_{1}(0),x_{1}^{c}(0))$ such that

[TABLE]

We claim that the trajectory starting from $y(0)$ will always leave the ball with radius $r_{1}$. Suppose this is not true, i.e., $\|y(t)\|\leq r_{1}$ for all $t\geq 0$. Then $\|(x_{1}(t),u_{1}(k\tau))\|\leq r_{0}$ for all $t\in[k\tau,k\tau+\tau]$ and $k\in\mathbb{Z}^{\geq 0}$ and $V(k+1)\geq(1+\delta)^{2}V(k)\geq(1+\delta)^{2k}V(0)$. Then $\lim_{k\rightarrow\infty}V(k)=\infty$, which contradicts the boundedness of $y(t)$. ∎

Stability test

The characteristic polynomial $P(s)$ of the linearized system (11) is

[TABLE]

which is a polynomial whose coefficients depend on the choice of parameters $\mathbf{p}$ for the digital controller (1) in either of the representations (3)-(4). As we have shown in Theorem 3.6, if this polynomial has a root $s$ outside unit circle, i.e., $\|s\|>1$, then the closed-loop sampled-data nonlinear system is unstable, so we can eliminate that controller from the synthesis domain.

4. Digital Controller Synthesis

In this section we define the digital controller synthesis problem for SDSSs, using the same notation of Definition 2.1. Recall that we consider controller parameters $\mathbf{p}\in\mathbb{P}$ (where $\mathbb{P}$ is a hyperbox in $\mathbb{R}^{2L+1}$), a finite time horizon $T$, and a safety invariant $\phi(x)$ defined as a predicate over the SDSS state vector $x$ (a quantifier-free FOL formula over the theory of nonlinear real arithmetic). The synthesized controller should ensure safety with respect to the probability measure induced by the stochastic disturbance $\mathbf{d}(\cdot)$ and the measurement noise $\eta$. In particular, the probability of satisfying $\phi(x(t))$ for all $t\in[0,T]$ should be above a given, user-defined threshold $\vartheta$. We further constrain the controller search by limiting the maximum controller degree to $L$ (see (3) in Definition 2.1).

Below, we denote by $\theta$ the stochastic uncertainty due to the SDSS disturbances $\mathbf{d}$ and the measurement noise $\eta$ up to time $T$.

Definition 4.1 (Digital Controller Synthesis).

Given an SDSS, a time bound $T<\infty$, a maximum controller degree $L$, and a probability threshold $\vartheta\in(0,1)$, the digital controller synthesis problem is finding the degree $l^{*}$

[TABLE]

and controller parameters $\mathbf{p}^{*}\in C(l^{*})$ where

[TABLE]

and $\mathbb{P}_{l}$ is the parameter space of controllers of degree $l$. For clarity we indicate that $x$ depends (indirectly) on the controller parameters $\mathbf{p}$ and on the stochastic uncertainty $\theta$. If $C(l)=\emptyset$ for all $l\leq L$, then we say the problem is infeasible.

In general the above synthesis problem is very hard to solve exactly due to the presence of nonlinearities, ordinary differential equations (ODEs) introduced by the SDSS dynamics, and multi-dimensional integration for computing probabilities. In fact, a decision version of the digital controller synthesis problem (i.e., given $l\leq L$ decide whether $C(l)$ is nonempty), is easily shown to be undecidable. While Satisfiability Modulo Theory (SMT) approaches, e.g., (Gao et al., 2012), can now in principle handle nonlinear arithmetics via a sound numerical relaxation, their computational complexity is exponential in the number of variables. In particular, our stochastic optimization problem is high-dimensional and currently infeasible for fully formal approaches, but it can be tackled using the mixed SMT-statistical approach of (Shmarov and Zuliani, 2016), which computes statistically and numerically sound confidence intervals. As such, we replace (exact) probability with empirical mean, thereby obtaining a Monte Carlo version of Definition 4.1.

Definition 4.2 (Digital Controller Synthesis - Monte Carlo).

Given an SDSS $H$, a time bound $T<\infty$, a maximum controller degree $L$, a probability threshold $\vartheta\in(0,1)$, and a confidence value $c\in(0,1)$. The Monte Carlo digital controller synthesis problem is finding

[TABLE]

and controller parameters $\mathbf{\hat{p}}^{*}\in\hat{C}(\hat{l}^{*})$ where

[TABLE]

and $\boldsymbol{\theta}_{K}=(\theta_{1},\ldots,\theta_{K})$ is a finite-dimensional random vector in which each $\theta_{i}$ is independent and identically distributed from $\theta$, $\text{Prob}_{\boldsymbol{\theta}_{K}}$ is the product measure of $K$ copies of the probability measure of $\theta$, and $\bar{X}(\mathbf{p},\boldsymbol{\theta}_{K})=\frac{1}{K}\sum_{i=1}^{K}I\left\{\forall t\in[0,T]\;\;\phi(x(\mathbf{p},\theta_{i},t))\right\}$, with $I\{\cdot\}$ being the indicator function.

Note that by the law of large numbers when $K\rightarrow\infty$ and the confidence value $c$ is sufficiently close to one, we have that $\hat{C}(l)$ approximates $C(l)$ arbitrarily well. Also, note that the constraints in Definition 4.2 are simpler to solve than those in 4.1, as they do not involve absolutely precise integration of the probability measure — we only require to decide the constraints with some statistical confidence $c$. However, even with this simplification a decision version of the Monte Carlo digital controller synthesis problem (i.e., deciding whether $\hat{C}(l)=\emptyset$) remains undecidable when plants with nonlinear ODEs are involved. Intuitively, that is because evaluating the elements of $\hat{C}(l)$ amounts to solving reachability, which is well known to be an undecidable problem for general nonlinear systems. Hence, one can only solve the Monte Carlo controller synthesis problem approximately, and that is what we aim to do in the next Section. Finally, note that we can generate sample realizations of $\theta$: recall from Section 2 that the disturbance $\mathbf{d}(\cdot)$ is defined from a finite number of random variables. This implies that $\theta$ is finite-dimensional as the number of random noise variables is also finite due to the sampling period.

5. Synthesis Algorithm

In this section we present an algorithm for approximately solving the Monte Carlo controller synthesis problem of Definition 4.2. Our synthesis algorithm starts from controllers with degree $l=0$ and iteratively increase $l$ until the constraint $\hat{C}(l)\neq\emptyset$ is satisfied or $l$ reaches a maximum value.

The synthesis algorithm, summarized in Algorithm 1, consists of two nested loops. The inner loop (lines 1-9) consists of two main stages: optimization and verification. Procedure optimize (line 1) aims at finding controller parameters $\mathbf{p}$ that (approximately) maximizes the empirical probability that the closed-loop system with a discrete-time version of the plant satisfies property $\phi$ over the finite time horizon $[0,T]$; optimize also returns an approximate confidence interval (CI) $[a,b]$ for such probability. Optimization is based on the cross-entropy algorithm, but our approach could work with different black-box optimization algorithms too, such as Gaussian process optimization (Rasmussen, 2004) and particle swarm optimization (Kennedy, 2011). Then, procedure verify (line 1) checks the candidate controller $\mathbf{p}$ in closed-loop with the original (continuous-time) plant model and computes a precise CI $[a^{\prime},b^{\prime}]$ for $\text{Prob}_{\boldsymbol{\theta}_{K}}\{\bar{X}(\mathbf{p},\boldsymbol{\theta}_{K})\geq\vartheta\}$. The reason for using the continuous-time plant only in verify is due to the high computational complexity of validated numerical ODE solving compared to solving its discrete-time approximation. The interval returned by verify is compared against the current best verified interval, which is then updated accordingly (line 1).

The procedures optimize and verify are iterated until the approximate CI (for the discrete-time plant) $[a,b]$ and the verified CI (for the continuous-time plant) $[a^{\prime},b^{\prime}]$ overlap to a certain length, or a maximum number of iterations is reached (line 9). In the outer loop, if $a^{\prime}\geq\vartheta$ (i.e., the LHS of the verified confidence interval is larger than $\vartheta$) then $\mathbf{p}$ is a witness for $\hat{C}(l)\neq\emptyset$, with probability at least $c$. Therefore, $\mathbf{p}$ approximately solves the digital controller synthesis problem of Definition 4.2, and the algorithm terminates. (The approximation lies in the fact that we cannot guarantee that the synthesized controller has minimum degree.) Otherwise ($a^{\prime}<\vartheta$), we increase the controller degree $l$ up to a maximum degree $L$.

In the inner loop, line 1 improves the approximation of the closed-loop system used in optimize. This can be any adjustment to the ODE solver complexity (e.g., increasing the Taylor series order). In our case, it corresponds to increasing the number of time points used for ODE integration. We next explain both optimize and verify in more detail.

Procedure optimize

In Algorithm 2 we give the pseudocode for the optimize function (line 1 of Algorithm 1). It implements a modified cross-entropy (CE) optimization algorithm (Rubinstein, 1999) that repeatedly samples (from the CE distribution of) controller parameters, evaluates their performance, and guides the search towards parameters that increase the safety probability (i.e., probability of satisfying $\phi$ over $[0,T]$). Sample performance is computed first by the stability check using Theorem 3.6 (line 2 of Algorithm 2). If a controller does not pass the test, i.e., it is necessarily unstable, it is rejected. Otherwise we compute a CI for the probability (over $\theta$) of the closed-loop system to satisfy the invariant $\phi$ (line 2). For this purpose, we consider a discrete-time version of the plant, simulated via an approximate ODE solver based on the first term of the Taylor series expansion. The time steps $0=t_{1}<\cdots<t_{G}=T$ for ODE integration are obtained by discretizing the time between the controller sampling points (defined by $\tau$ — see Definition 2.1) using the discretization parameter $\mathbf{m}$.

To compute the CI, our implementation uses sequential Bayesian estimation for efficiency reasons (Shmarov and Zuliani, 2015), but other standard statistical techniques may also be employed (e.g., the Chernoff-Hoeffding bound). After an adequate number of controller parameters are sampled and evaluated, the best performing sample is chosen (line 2), and the CE distribution is updated accordingly (line 2, see (Rubinstein, 1999) for more details). This is repeated until a maximum number of iterations is reached.

Procedure verify

We use the ProbReach tool (Shmarov and Zuliani, 2015) to compute a CI for the probability that a candidate digital controller in closed-loop with the plant satisfies the time-bounded invariant $\phi$ (line 1 of Algorithm 1). This step is necessary since the candidate controller has been obtained using an approximate, discrete-time solver for simulating the (continuous) plant dynamics, while ProbReach uses instead an SMT solver (Gao et al., 2013) to handle the plant dynamics in a sound manner. In particular, ProbReach allows to derive a guaranteed confidence interval for $\text{Prob}_{\boldsymbol{\theta}_{K}}\{\bar{X}(\mathbf{p},\boldsymbol{\theta}_{K})\geq\vartheta\}$ with confidence $c$, where $K$ is the number of samples of the Monte Carlo synthesis problem (see Definition 4.2). (We note that in general $K$ will depend on the size of the interval and the confidence $c$.)

Procedure verify consists of two steps. The first step builds a hybrid system (the model format accepted by ProbReach) representing the closed-loop system under the candidate controller. In the second step, verify invokes ProbReach with three parameters: the hybrid system, and the required minimum size $\xi$ and confidence $c$ of the confidence interval to compute. We remark that the size of the confidence interval cannot be guaranteed in general (Shmarov and Zuliani, 2016) because of the undecidability of reasoning about nonlinear arithmetics. As such, the confidence interval returned by ProbReach via verify can be fully trusted from both the statistical and numerical viewpoints: while the interval size might be larger than $\xi$, the confidence is guaranteed to be at least $c$, as the sampled controllers are evaluated by SMT and verified numerical techniques.

Theorem 5.1.

Let $S$ be an SDSS for which the synthesis problem of Definition 4.1 is feasible for a given $\vartheta\in(0,1)$ and controller degree $l\leq L$. Suppose that Algorithm 1, with parameters $S,L,\vartheta$ and $c\in(0,1)$, returns a controller $\mathbf{p}$ and an interval $[a,b]$ such that $a\geq\vartheta$. Then $\mathbf{p}$ is a solution of the Monte Carlo problem of Definition 4.2 for $S,\vartheta$, and $c$.

Proof (sketch).

It suffices to note that Algorithm 1 terminates either by finding a controller of minimal degree whose safety probability is at least $\vartheta$, with confidence $c$, or by finding a controller of degree $l$ with the best confidence interval produced by verify. By hypothesis a controller of degree $l$ that satisfies $\phi$ with probability larger than $\vartheta$ exists, and Algorithm 1 returns an interval whose LHS is larger than $\vartheta$. Therefore, by verify, we know that this interval has confidence $c$. ∎

We remark that Algorithm 1 can output sub-optimal controllers: this is unavoidable when using stochastic optimization methods such as cross-entropy. However if the controller design problem is feasible, multiple restarts of Algorithm 1 will eventually find the minimal controller with probability 1.

6. Case studies and evaluation

We evaluate our approach on three case studies: a model of insulin control for Type 1 diabetes (T1D) (Hovorka, 2011), also known as the artificial pancreas (AP), a model of a powertrain control system (Jin et al., 2014), and a quadruple-tank process. For all case studies we use the following input parameters for Algorithm 1: $\xi=0.05$, $c=0.99$ and $\alpha=0.5$. The experiments were performed on a 32-core Intel 2.90GHz system running Ubuntu.

Digital PID controllers.

While our algorithm can synthesize any digital controller as per Definition 3.1, we here exemplify its use via proportional-integral-derivative (PID) controllers, one of the most popular control techniques. A PID controller output is the weighted sum of three terms: the error itself weighted with $K_{P}$, its rate of change weighted with $K_{D}$, and accumulated error weighted with $K_{I}$. The input/output equation of a digital PID controller is

[TABLE]

Essentially, the controller needs to store the previous value of the input and the previous two values of the error. In the following case studies, we focus on the synthesis of controllers in the PID form, hence we consider a maximum degree $L=2$.

6.1. Artificial pancreas

The AP is a system for the automated delivery of insulin therapy that is required to keep blood glucose (BG) levels of diabetic patients within safe ranges, typically between 4-11 mmol/L. A so-called continuous glucose monitor (CGM) sends BG measurements to a control algorithm that computes the adequate insulin input. PID control is one of the main techniques (Steil et al., 2011), and is also found in commercial devices (Kanderian Jr and Steil, 2014).

Meals are the major disturbance in insulin control, which make full closed-loop control challenging. Our approach is therefore well suited to solve this problem because it can synthesize controllers attaining arbitrary safety probability by minimizing the impact of such disturbances. To model insulin and glucose dynamics, we employ the nonlinear model of Hovorka et al. (Hovorka et al., 2004), considered as one of the most faithful models. The plant has nine state variables describing insulin and glucose concentration in different physiological compartments. We evaluate the system for a time bound of 24 hours.

In our SDSS model, we consider three meals (respectively represent breakfast, lunch and dinner) with random timing and random amount, expressed by the following normally-distributed parameters: the amount of carbohydrates (CHO) of each meal in grams, $D_{G_{0}}\sim\mathcal{N}(50,100)$, $D_{G_{1}}\sim\mathcal{N}(70,100)$ and $D_{G_{2}}\sim\mathcal{N}(60,100)$, and the waiting times between meals, $T_{1}\sim\mathcal{N}(300,100)$ and $T_{2}\sim\mathcal{N}(300,100)$. The corresponding disturbance input is given by:

[TABLE]

The system output $y(t)$ is the CGM measurement (performed every 5 minutes), given by the equation $y(t)=C(t)+\eta(t)$, where $C$ is the state variable for interstitial glucose and $\eta(t)$ is white Gaussian sensor noise with standard deviation $0.25$.

The control input $u(t)$ is the insulin infusion rate computed by the PID controller. The tracking error is defined as $e(t)=r(t)-y(t)$ with the constant reference signal $r(t)=6.11$ mmol/L. The total infusion rate is given by $u(t)+u_{b}$ where $u_{b}$ ($\approx 0.05548$) is the basal insulin, i.e., a low and continuous dose to regulate glucose outside meals. The value of $u_{b}$ is chosen to guarantee a steady-state BG value equals to $r(t)$ in absence of meals. This steady state is used as the initial state of the system.

Safety property

Insulin control seeks to prevent hyperglycemia (BG above 11 mmol/L) and hypoglycemia (BG below 4 mmol/L). Hypoglycemia happens when the controller overshoots its dose and has more severe health effects than hyperglycemia, which is tolerated to a small extent after meals. For this reason we consider a safe BG range of $[4,16]$ mmol/L, which strictly avoids hypoglycemia and allows for some post-meal hyperglycemia tolerance. In addition, we want that the glucose level stays close to the reference signal towards the end of the 24 hours (1440 minutes). Our invariant is given by:

[TABLE]

where $G$ is the state variable for the BG concentration.

In the synthesis algorithm, we use a probability threshold of $\vartheta=0.95$ (we want to satisfy the above invariant with probability not below 95%), confidence $c=0.99$, confidence interval size $\xi=0.05$, and $\alpha=0.5$.

Synthesis results

Table 1 shows the PID controllers synthesized at each iteration of the algorithm. The domain of controller parameters was chosen as follows: $K_{p}\in[-10^{-2},10^{-3}]$, $K_{i}\in[-10^{-5},10^{-6}]$ and $K_{d}\in[-1,10^{-1}]$. Even though none of the synthesized controllers achieves the probability threshold $\vartheta=0.95$, the degree-2 controller (PID) is very close to satisfying the property, with a $99\%$-confidence interval of $[0.94242,0.99242]$.

To better understand the performance of the controllers, we analyze their behavior on 1,000 Monte Carlo executions of the system. Results, reported in Figure 3, evidence that hyper- and hypo-glycemia episodes are never sustained.

6.2. Powertrain system

We consider the automotive air-fuel control system adapted from the powertrain control benchmark in (Jin et al., 2014). The plant model consists of a system of three nonlinear ODEs describing the dynamics of the engine in relation to the throttle air dynamics, intake manifold and air-fuel path.

The system has two exogenous inputs, captured by the disturbance vector $\mathbf{d}(t)=[\omega(t)\ \theta_{in}(t)]^{T}$, where $\omega$ (rad/s) is the engine speed ($\omega\sim\mathcal{N}(105,4)$), and $\theta_{in}$ (degrees) is the throttle angle. $\theta_{in}(t)$ is defined as a pulse train wave with random amplitude $a\sim\mathcal{N}(30.6,25)$ and period $\zeta=4$:

[TABLE]

The noisy plant output is $y(t)=\lambda(t)+\eta(t)$, where $\lambda(t)$ is the air/fuel ratio, and $\eta(t)\sim\mathcal{N}(0,0.0625)$.

The engine is controlled by a PID controller that seeks to maintain a constant air/fuel ratio equals to the stoichiometric value $\bar{\lambda}=14.7$, that is when the engine performs optimally. The tracking error is thus given by $e(t)=y(t)-\bar{\lambda}$. The control signal $u(t)$ determines the amount of fuel entering the system.

Safety property

We consider the following invariant

[TABLE]

where $\mu(t)=(\lambda(t)-\bar{\lambda})/\bar{\lambda}$ is the relative error from the setpoint. The first conjunct states that the air/fuel ratio should constantly be within $\pm 100\%$ of the ideal ratio $\bar{\lambda}$. The second conjunct states that whenever the input throttle angle $\theta_{in}$ rises (at time $t=0$) or falls ($t=\zeta/2$), the plant should settle within time $\zeta/8$ and remain in the settling region ($\pm 5\%$ around $\bar{\lambda}$) until the next rise or fall (happening after time $\zeta/2$). We set the probability threshold to $\vartheta=0.96$.

Synthesis results

Table 2 shows the PID controllers synthesized at each iteration of the algorithm. The domain of controller parameters was chosen as follows: $K_{p}\in[-0.1,0.5]$, $K_{i}\in[-0.05,0.2]$ and $K_{d}\in[-0.05,0.05]$. With our algorithm, we could synthesize a degree-2 controller (PID) satisfying the threshold. The optimal degree-1 controller has similar performance (both yield confidence intervals with RHS equals to 1), albeit below the threshold.

Compared to the AP case study, we observe that the powertrain model requires generating (and verifying) fewer candidate parameters. At the same time the dynamics of the powertrain system appear more challenging to control as the model requires more ODE integration steps (see column $\#_{o}(\#^{0}_{o})$ of Tables 2 and 1).

6.3. Quadruple-tank process

We consider a quadruple-tank process adapted from (Johansson, 2000), which consists of four interconnected water tanks. The process is illustrated in Figure 4. This model is an example of a multiple-input and multiple-output (MIMO) system with multivariable right half-plane zeros (Glad and Ljung, 2000) (such zeros bring performance limitations in control problems). We extended the deterministic model of (Johansson, 2000) to include uncertainties in the valve settings and random disturbances in the process that remove water from the tanks.

The process is controlled in a decentralized fashion, by which two digital controllers are designed for the input-output pairs $(u_{1},y_{1})$ and $(u_{2},y_{2})$, where $u_{1}$ and $u_{2}$ are the input voltages for the pumps, and $y_{1}$ and $y_{2}$ are the water level measurements obtained as $y_{1}=0.5\cdot h_{1}$ and $y_{2}=0.5\cdot h_{2}$, where $h_{1}$ and $h_{2}$ are the water levels in tanks $1$ and $2$, respectively. In this case study we assume that the pumps can only add water to the tanks (and cannot pump it out).

We consider a scenario where at time 0 and then twice after every minute we remove a random amount of water (which is model through reducing the corresponding water levels by a random value $\sim U(0,3)$) from tanks 1 and 2. Every time such a disturbance happens, the valves parameters are randomly reset to $\gamma_{1}\sim N(0.7,0.223)$ and $\gamma_{2}\sim N(0.6,0.223)$. The system is subject to a measurement noise modeled as a white Gaussian noise with variance 0.33.

Safety property

After each disturbance, we require that the system reaches the desired water levels in tanks 1 and 2 (within 1 centimeter above or below the corresponding set points $r_{1}=12.4$ and $r_{2}=12.7$) within 5 seconds, and that the water levels stay close to the setpoints for the remaining 55 seconds, before the next disturbance occurs. Also, all four water levels $h_{1},h_{2},h_{3},h_{4}$ must always stay in the interval $[0,20]$ and the input voltages $u_{1},u_{2}$ for both pumps must be in the range $[0,24]$.

Synthesis results

The domain of controller parameters was chosen as $K_{P_{1}}\in[-1,20]$, $K_{I_{1}}\in[-1,10]$, $K_{D_{1}}\in[-1,10]$, $K_{P_{2}}\in[-1,20]$, $K_{I_{2}}\in[-1,10]$, $K_{D_{2}}\in[-1,10]$. The controller synthesis results are presented in Table 3, which shows that we can obtain a confidence interval of up to $[0.94,0.99]$ for the safety probability by using two PI controllers (see third row of Table 3). Note that the performance of the controller is not improved by including the derivative terms $K_{D_{1}},K_{D_{2}}$ (see last two rows). This is due to the optimization algorithm which works by sampling a finite number of controller parameters and thus, might fail to explore parameter regions with better safety probability.

7. Related Work

Although recent papers (Abate et al., 2017b, a; Duggirala and Viswanathan, 2015) have addressed the synthesis of safe digital controllers for linear and deterministic systems, synthesis for the class of stochastic nonlinear systems that we consider has not yet been tackled.

While the approach proposed in this paper can be applied in principle to any digital controller, our case studies focus on PID controllers. Several methods have been proposed for the synthesis of PID controller for nonlinear and stochastic plants (Su et al., 2005; Fliess and Join, 2013; Guo and Wang, 2005; He and Liu, 2011; Duong and Lee, 2012; Shmarov et al., 2017; Alimguzhin et al., 2017). However, none of these methods can provide safety guarantees beside the work of (Shmarov et al., 2017) (discussed at the end of the section) .

We pursue a different direction with respect to the classical solutions proposed in the literature, by providing a framework for the efficient synthesis of discrete-time digital controller for nonlinear and stochastic plants that are safe and robust (with respect to a probabilistic reachability property) by construction.

The control synthesis problem considered here requires to solve a parameter synthesis problem over a closed-loop system modeled as a stochastic nonlinear system. In contrast with existing parameter synthesis techniques for stochastic and continuous nonlinear systems (Bartocci et al., 2015; Bortolussi and Sanguinetti, 2015; Haghighi et al., 2015), our approach, which targets specifically digital controllers, takes advantage of a computationally efficient stability check that rules out unstable controller candidates, thereby reducing the computational effort.

The problem of controller synthesis under safety requirements has been investigated mostly for Model Predictive Control (MPC) (Camacho and Alba, 2013) whose goal is to find the control input that optimizes the predicted performance of the closed-loop system up to a finite horizon. The work of (Karaman et al., 2008; Raman et al., 2014; Kim et al., 2017; Wongpiromsarn et al., 2012; Pant et al., 2017; Li et al., 2017; Sadigh and Kapoor, 2016) consider safety requirements expressed as temporal logic formulas, and synthesize MPC controllers that optimize the robust satisfaction of the formula (Donzé and Maler, 2010) (i.e., a quantitative measure of satisfaction). MPC is an online method that requires solving at runtime often computationally expensive optimization problems. In contrast, our approach performs controller synthesis at design time.

The closest paper to our own is the work of (Shmarov et al., 2017), where the authors have recently proposed a method to synthesize continuous-time PID controllers for nonlinear stochastic plants such that the resulting closed-loop system satisfies safety and performance requirements specified as bounded reachability properties. That method works under the assumptions that the system can measure the output of the plant continuously and without sensing noise. However, this is not realistic in the majority of embedded systems whose operations are governed by a discrete-time clock and where sensor noise is unavoidable.

8. Conclusions

The synthesis of digital controllers for cyber-physical systems with nonlinear and stochastic dynamics is a challenging problem, and for such systems, no automated methods currently exist for deriving controllers with rigorous and quantitative safety guarantees. In this paper, we have presented a solution to this problem based on two key contributions: a method to check the candidate controllers for the sampled-data nonlinear system with respect to stability; and a two-stage synthesis algorithm that alternates between a fast candidate generation phase (based on Monte-Carlo sampling and non-validated ODE solving) and a verification phase where we derive numerically and statistically valid confidence intervals on the safety probability of the closed-loop system. With this method, we managed to synthesize controllers for three nonlinear systems (artificial pancreas, powertrain, and quadruple-tank process) characterized by large stochastic disturbances and sensing noise. As future work, we plan to extend our method to hybrid systems and controllers with fixed-point precision.

Appendix A Stability Requirement

The following lemma provides a bound on exponent of a matrix and is is used in proving our main theorem.

Lemma A.1.

For any matrix $A$ and any $\epsilon>0$, there is a constant $C(\epsilon)$ such that

[TABLE]

with $a=\max_{i}Re(\lambda_{i}(A))$.

Proof.

We use the Jordan normal form $J$ of $A$ that satisfies $A=P^{-1}JP$ for some invertible matrix $P$:

[TABLE]

where $J_{i}$ is the $i^{\text{th}}$ block of $J$ with size $k_{0}\in\mathbb{N}$ and can be written as $J_{i}=\lambda_{i}I+N$. Matrix $N$ is a matrix of all zeros except identities on the first superdiagonal. Since $N^{m}=0$ for all $m>k_{0}$, we have the following for $J_{i}$:

[TABLE]

The claim is true by taking $C(\epsilon):=\|P^{-1}\|_{2}\|P\|_{2}\max_{i}C_{i}(\epsilon)$. ∎

Proof of Lemma 3.4.

Define the function $W(t):=\|x_{1}(t)\|^{2}$,

[TABLE]

We write this inequality in terms of $\sigma(t):=\sqrt{W(t)}$ as

[TABLE]

[TABLE]

Rename $L:=\|A\|+\gamma$, $L_{1}:=\|BC_{c}\|+\gamma\|C_{c}\|$ and $L_{2}:=\|BD_{c}C\|+\gamma\|D_{c}C\|+\gamma^{2}\|D_{c}\|$, to get

[TABLE]

with functions

[TABLE]

where $L,L_{1},L_{2}$ depend on $\gamma$ as defined above. ∎

Proof of Lemma 3.5.

Using the assumption of $\|(x_{1}(t),u_{1}(k\tau))\|\leq r(\gamma)$ for all $t\in[k\tau,(k+1)\tau]$, we get

[TABLE]

Then we employ Lemma 3.4 to get

[TABLE]

The claim holds with any $\hat{h}_{1},\hat{h}_{2}$ with

[TABLE]

which can be selected as the following according to Lemma A.1, with $\Gamma:=C(1)$ and $\alpha:=a+1$,

[TABLE]

∎

Appendix B Gluco-regulatory ODE model

The model consists of three subsystems:

• •

Glucose Subsystem: it tracks the masses of glucose (in mmol) in the accessible ($Q_{1}$) and non-accessible ($Q_{2}$) compartments, $G$ (mmol/L) represents the glucose concentration in plasma, $EGP_{0}$ (mmol/min) is the endogenous glucose production rate and $U_{G}(t)$ (mmol/min) is the glucose absorption rate from the gut.

• •

Gut absorption: this subsystem uses a chain of two compartments, $G_{1}$ and $G_{2}$ (mmol), to model the absorption dynamics of ingested food, given by the disturbance $D_{G}(t)$. $A_{g}$ is the CHO bio-availability. $t_{maxG}$ (min) is the time of maximum appearance rate of glucose.

• •

Interstitial glucose: $C$ is the subcutaneous glucose concentration (mmol/L) detected by the CGM sensor and has a delayed response w.r.t. the blood concentration $G$.

• •

Insulin Subsystem: it represents absorption of subcutaneously administered insulin. It is defined by a two-compartment chain, $S_{1}$ and $S_{2}$ measured in U (units of insulin), where $u(t)$ (U/min) is the administration of insulin computed by the PID controller, $u_{b}$ (U/min) is the basal insulin infusion rate and $I$ (U/L) indicates the insulin concentration in plasma.

• •

Insulin Action Subsystem: it models the action of insulin on glucose distribution/transport, $x_{1}$, glucose disposal, $x_{2}$, and endogenous glucose production, $x_{3}$ (unitless).

The model parameters are given in Table 4.

[TABLE]

Appendix C Fuel Control System Model

The dynamics of the engine (plant) are given by the following set of ODEs:

[TABLE]

where $p$ (bar) is the intake manifold pressure; $\theta$ (degrees) is the throttle angle; $\lambda$ is the air/fuel ratio; $\theta_{in}$ (degrees) is the throttle angle input disturbance; $\hat{\theta}$ is the throttle plate angle and is defined by:

[TABLE]

$\dot{m}_{c}$ (g/s) is the air inflow rate to cylinder and is defined by:

[TABLE]

$\omega$ (rad/s) is the engine speed disturbance; $\dot{m}_{af}$ is the inlet air mass flow rate, defined by:

[TABLE]

$F_{c}$ is the commanded fuel input defined as $F_{c}=(1+u(t))\dot{m}_{c}/\bar{\lambda}$, where $u(t)$ is the control input and $\bar{\lambda}$ is the ideal air/fuel ratio.

Parameter values are: $c_{1}=0.41328$, $c_{2}=-0.366$, $c_{3}=0.08979$, $c_{4}=-0.0337$, $c_{5}=0.0001$, $c_{6}=2.821$, $c_{7}=-0.05231$, $c_{8}=0.10299$, $c_{9}=-0.00063$, $c_{10}=1$, $c_{12}=0.9$, $c_{25}=1$, $c_{26}=4$.

Appendix D Model of the Quadruple-Tank Process

The dynamics of the Quadruple-Tank Process are given by the following set of ODEs (Johansson, 2000):

[TABLE]

where $h_{i},a_{i},A_{i}$, $i\in\{1,2,3,4\}$ are the water level, cross-section of the outlet hole, and cross-section of tank $i$, respectively. Inputs $u_{1},u_{2}$ indicate the voltages applied to the pumps and the corresponding flows are $k_{1}u_{1},k_{2}u_{2}$. The parameters $\gamma_{1},\gamma_{2}\in(0,1)$ show the settings of the valves. The flow to tank $1$ is $\gamma_{1}k_{1}u_{1}$ and the flow to tank $4$ is $(1-\gamma_{1})k_{1}u_{1}$ (similarly for the other two tanks). The acceleration of gravity is denoted by $g$. The water levels of tanks $1,2$ are measured by sensors as $k_{c}h_{1},k_{c}h_{2}$. The parameter values are: $A_{1}=A_{3}=28\,\text{cm}^{2}$, $A_{2}=A_{4}=32\,\text{cm}^{2}$, $a_{1}=a_{3}=0.071\,\text{cm}^{2}$, $a_{2}=a_{4}=0.057\,\text{cm}^{2}$, $k_{c}=0.5\,\text{V}/\text{cm}$, $g=9.81\,\text{m}/\text{s}^{2}$. We have chosen the steady state values $h_{1}^{0}=12.4\,\text{cm},h_{2}^{0}=12.7\,\text{cm},h_{3}^{0}=1.8\,\text{cm},h_{4}^{0}=1.4\,\text{cm},u_{1}^{0}=3.00\,\text{V},u_{2}^{0}=3.00\,\text{V},k_{1}=3.33\,\text{cm}^{3}/\text{Vs}$, and $k_{2}=3.35\,\text{cm}^{3}/\text{Vs}$.