Formal Synthesis of Analytic Controllers for Sampled-Data Systems via Genetic Programming

TL;DR

This paper introduces an automated method using genetic programming to synthesize formal controllers for nonlinear sampled-data systems, ensuring safety and reachability, verified through SMT solvers, applicable beyond polynomial systems.

Contribution

It develops a novel, generalizable genetic programming approach for synthesizing safety-critical controllers with formal guarantees for nonlinear sampled-data systems.

Findings

Successfully synthesized controllers for multiple systems

Verified correctness using SMT solvers

Applicable to non-polynomial systems

Abstract

This paper presents an automatic formal controller synthesis method for nonlinear sampled-data systems with safety and reachability specifications. Fundamentally, the presented method is not restricted to polynomial systems and controllers. We consider periodically switched controllers based on a Control Lyapunov Barrier-like functions. The proposed method utilizes genetic programming to synthesize these functions as well as the controller modes. Correctness of the controller are subsequently verified by means of a Satisfiability Modulo Theories solver. Effectiveness of the proposed methodology is demonstrated on multiple systems.

Tables

Table 1. Table 2: Polynomial systems and results for 10 runs. t 𝑡 t : total time, μ 𝜇 \mu : mean, σ 𝜎 \sigma : standard deviation.

System	Linear	2nd-order	3rd-order
$f(x,u)$	$\begin{array}{c}{\dot{x}}_1=x_2\\ {\dot{x}}_2=-x_1+u\end{array}$		$\begin{array}{c}{\dot{x}}_1=-10x_1+10x_2+u\\ {\dot{x}}_2=28x_1-x_2-x_1{x}_3\\ {\dot{x}}_3=x_1{x}_2-2.6667x_3\end{array}$
$S$	${[-1,1]}^2$	${[-1,1]}^2$	${[-5,5]}^2$
$I$	${[-0.5,0.5]}^2$	${[-0.5,0.5]}^2$	${[-1.2,1.2]}^3$
$G$	${[-0.1,0.1]}^2$	${[-0.05,0.05]}^2$	${[-0.3,0.3]}^3$
$𝒢$	$\{-1,0,1\}$	$\{-1,0,1\}$	$\{-100,-50,-5,0,5,50,100\}$
$\epsilon ,h$	$(2,1),0.01s$	$(7,0),0.01s$	$(3800,6800,1900),0.001s$
$\begin{array}{c}\\ \text{\# gen.}\\ t\left[s\right]\end{array}$

Table 2. Table 3: Nonpolynomial systems and results for 10 runs. t 𝑡 t : total time, μ 𝜇 \mu : mean, σ 𝜎 \sigma : standard deviation

System	Pendulum	Pendulum on a cart
$f(x,u)$	$\begin{array}{c}{\dot{x}}_1=x_2\\ {\dot{x}}_2=-\left(\frac{b}{J}+\frac{K^2}{J{R}_a}\right)x_2-\frac{mlg}{J}\sin \left(x_1\right)+\frac{K}{J{R}_a}u\\ \\ m=5.50\cdot {10}^{-2}\text{ kg},l=4.20\cdot {10}^{-2}\text{ m},\\ J=1.91\cdot {10}^{-4}\text{ kg}{\text{m}}^2,g=9.81\text{ m}/{\text{s}}^2,\\ K=5.36\cdot {10}^{-2}\text{ Nm/A},R_a=9.50\Omega ,\\ b=3.0\cdot {10}^{-6}\text{Nms}\end{array}$	$\begin{array}{c}{\dot{x}}_1=x_2\\ {\dot{x}}_2=\frac{g}{l}\sin \left(x_1\right)-\frac{b}{m{l}^2}{x}_2+\frac{1}{ml}\cos \left(x_1\right)u\\ \\ g=9.8\text{ m}/{\text{s}}^2,b=2\text{Nms}\\ l=0.5\text{ m},m=0.5\text{ kg}.\end{array}$
$S$	$[-2\pi ,2\pi ]\times [-100,100]$	$[-2\pi ,2\pi ]\times [-10,10]$
$I$	$[-\pi ,\pi ]\times [-10,10]$	${[-0.5,0.5]}^2$
$G$	$[-1.0,-0.5]\times [-1.0,1.0]$	$[-0.25,0.25]\times [-1,1]$
$𝒢$	$\{-10,-5,0,5,10\}$	$\{-6,-2,0,2,6\}$
$\epsilon ,h$	$(600,12700),0.001s$	$(200,3200),0.001s$
$\begin{array}{c}\\ \text{\# gen.}\\ t\left[s\right]\end{array}$

Table 3. Table 4: Production rules 𝒫 𝒫 \mathcal{P}

$𝒩$	Rules
	$::=x_1-x_{c,1}\left|\dots \right|x_n-x_{c,n}$
	$::=$ Random Real $\in [-10,10]$
	,

Table 4. Table 5: Results for 10 runs for the Pendulum on a cart system without pre-defining 𝒢 𝒢 \mathcal{G} .

	Min	Max	$\mu$	$\sigma$
# gen.	4	14	7.6	2.84
$t\left[s\right]$	135.59	536.12	267.82	124.48

Equations

$$\dot{\xi }(t)=f(\xi (t),u(t)),$$

$$\begin{array}[]{r}\forall\xi(t_{0})\in I,\exists T,\forall t\in[t_{0},T]:\xi(t)\in S\wedge\xi(T)\in G.\end{array}$$

$$\begin{array}[]{r}\!\!\!\forall\xi(t_{0})\in I,\exists T,\forall t\geq t_{0},\forall\tau\geq T\!:\!\xi(t)\in S\wedge\xi(\tau)\in\!G.\!\!\!\end{array}$$

$$\mathcal{G}=\{g_q:\mathcal{X}\to \mathcal{U}\mid q\in Q\}.$$

$$\forall x\in I:$$

$$\forall x\in \partial S:$$

$$\forall x\in A\backslash G,\exists q\in Q,\forall z\in R_q(x):$$

$$\left\{\begin{array}[]{rl}u(t)&=g_{q_{k}}(\xi(t_{k})).\\ q_{k}(t_{k})&=\operatorname*{arg\,min}\limits_{q\in Q}\max\limits_{z\in R_{q}(\xi(t_{k}))}\dot{V}_{q}(\xi(t_{k}),z)\end{array}\right.$$

$$\forall x\in \partial G:V(x)>\beta$$

$$\forall x\in G\backslash int(B),\exists q\in Q,\forall z\in R_q(x):{\dot{V}}_q(x,z)\le -\gamma$$

$$r_q(x,\tau ,e)=x+\tau f(x,g_q(x))+\frac{1}{2}{\tau }^{2}e,$$

$$R_q(s)=\bigcup _{(\tau ,e)\in E}{r}_q(s,\tau ,e)$$

$${\epsilon }_i=\underset{(x,u)\in \mathcal{X}\times \mathcal{U}}{\max }\left|\frac{\partial f_i(x,u)}{\partial x}{f}_i(x,u)\right|.$$

$$(\forall s\in C_i){\varphi }_i(s)\ge 0,i=1,\dots ,5.$$

$$e_{{\varphi }_i}$$

$$f_{\mathrm{samp},{\varphi }_i}:=(1+e_{{\varphi }_i}{)}^{-1},i=1,...,5.$$

$$f:=\sum _{i=1}^j{w}_i{f}_{\mathrm{samp},{\varphi }_i}+\sum _{i=1}^j{f}_{\mathrm{SMT},{\varphi }_i},j=3\vee j=5.$$

$$V^{\prime}(x)=V(x)-\max\Big{(}\max_{x\in I_{\mathrm{samp}}}(V(x)),0\Big{)},$$

$$\forall x\in S\backslash G:V(x)\ge V(x_c)$$

$$\begin{array}[]{r}\forall x\in D,\forall q\in Q:\max\limits_{z\in R_{q}(x)}\dot{V}_{q}(x,z)>-\gamma\implies\\ \min\limits_{p\in Q}(\dot{V}_{p}(x,x)+\alpha_{p}(x))<\dot{V}_{q}(x,x)+\alpha_{q}(x),\end{array}$$

$$q_k(t_k)=\underset{q\in Q}{\mathrm{argmin}}({\dot{V}}_q(\xi (t_k),\xi (t_k))+{\alpha }_q(\xi (t_k))).$$

$$V(x)=-4015.83+10.8526x_1^{\prime }+199.048x_1^{\prime 2}+0.311673x_2+18.8116x_1^{\prime }{x}_2+2.23916x_2^2,$$

$$u(t_k)=\max (-6,\min (6,g_{q_k}(\xi (t_k)))).$$

$$V$$

$$\mathcal{G}$$

Full text

Supported by NWO Domain TTW under the CADUSY project #13852.

The original version of this article has been accepted to CDC 2018. This version contains minor corrections.

1 Introduction

Modern controller design for nonlinear continuous systems often involves both reachability and safety specifications. Furthermore, digital controller implementations typically impose that states are measured periodically and that control signals are held constant in between sampling. This paper proposes an approach to automatically synthesize periodically switched state feedback controllers for a special subclass of safety and reachability specifications for nonlinear sampled-data systems.

Two popular paradigms for automatic controller synthesis for reachability and safety specifications are: 1) abstraction and simulation, and 2) Control Lyapunov functions (CLF) and Control Barrier Functions (CBF).

The first approach abstracts the infinite system to a finite one, which simplifies the formal controller synthesis for temporal logic specifications [1]. For nonlinear systems, tools implementing this approach include PESSOA [2], SCOTS [3] and CoSyMa [4]. The second approach deals with the system as an infinite system. Control Lyapunov functions [5] and Control Barrier Functions [6] are design tools for stabilization and safety specifications respectively. In [7] and [8] attempts are made to combine both CLFs and CBFs. Automatic synthesis of these functions is often done by posing the problem as a sum of squares (SOS) problem, which can be solved through convex optimization, see e.g. [9] and [10]. Drawbacks of the abstraction and (bi-)simulation approach are that it requires discretization of the state space and that the resulting controller is often an enormous look-up table in the form of a sparse matrix or a binary decision diagram (BDD). On the other hand, the SOS programming paradigm is limited to polynomial systems. Although reformulation of some nonpolynomial systems to an SOS formulation exists, e.g. [9, 11] and references therein, polynomial Lyapunov functions can be too restrictive, as global asymptotical stability of a polynomial system does not imply the existence of a polynomial Lyapunov function [12].

To overcome these limitations, we propose a framework which uses genetic programming (GP) in combination with a Satisfiability Modulo Theories (SMT) solver. GP is an evolutionary algorithm which evolves encoded representations of symbolic functions [13], rather than just fitting optimal parameters given a predefined structure. An SMT solver is a tool which uses a combination of background theories to determine whether a first-order logic formula can be satisfied [14]. Our approach uses a CLF-like function and a predefined switching law that infers a reachability and safety specification. The proposed framework uses GP to automatically generate both candidate CLFs and optionally the controller modes of a periodically switched state feedback controller. The SMT solver is subsequently used to formally verify the candidate solutions. By using GP, we allow ourselves to search for solutions that include nonpolynomial functions. Furthermore, the synthesized controllers are expressed as analytic expressions that are significantly more compact than BDDs returned by abstraction-based methods.

This work is a follow-up to [15], in which also a combination of GP and SMT solvers is used. The main contributions of this paper are: 1) synthesis w.r.t. a predefined periodic sampling time, rather than arbitrary switching with a (more conservative) minimum dwell-time and 2) the use of a different and less conservative CLBF. Additionally, more benchmark examples are provided. Other related work is found in [16], in which robust CLFs for switched systems with reach-while-stay (RWS) specifications are synthesized using a counterexample-guided synthesis. However, in [16] the controller modes are pre-specified, while in our approach these modes can also be discovered automatically, eliminating the need for prior input space discretization. Furthermore, this paper extends the set of specifications to include invariance of the goal set. Finally, similar to [15], the theoretical lower bounds on the minimum dwell-times reported in [16] are often very conservative.

Notation

Let $\mathbb{Z}_{\geq 0}=\left\{0,1,2\dots\right\}$. Let us denote the boundary and the interior of a set $D$ with $\partial D$ and $int(D)$ respectively. The image and inverse image of set $A$ under $f$ are denoted by $f[A]$ and $f^{-1}[A]$. Finally, the Euclidean norm is denoted by $\|\cdot\|$.

2 Problem definition

In this paper we design sampled-data state feedback controllers for nonlinear continuous-time systems described by

[TABLE]

where the variables $\xi(t)\in\mathcal{X}\subseteq\mathbb{R}^{n}$ and $u(t)\in\mathcal{U}\subseteq\mathbb{R}^{m}$ denote the state and input respectively. Due to the sampled-data nature of the controller, $u(t)=g(x(t_{k}))$, $\forall t\in[t_{k},t_{k}+h)$, where $h>0$ denotes a constant sampling time.

2.1 Control specification

Given a compact safe set $S\subseteq\mathcal{X}$, compact initial set $I\subset S$ and compact goal set $G\subset S$, we consider the following specifications:

1. $\mathtt{CS}_{1}$

Reach while stay (RWS): all trajectories starting in $I$ eventually reach $G$, while staying within $S$:

[TABLE]

1. $\mathtt{CS}_{2}$

Reach and stay while stay (RSWS): all trajectories starting from $I$ eventually reach and stay in $G$, while staying within $S$:

[TABLE]

This paper addresses the following problem:

Problem 2.1**.**

Given the compact sets $(S,I,G)$ and system (1), synthesize a sampled-data state feedback controller $u(t)=g(x(t_{k}))$ such that the closed-loop system satisfies specification $\mathtt{CS}_{1}$ or $\mathtt{CS}_{2}$.

We propose to solve Problem 2.1 by using a periodically switched controller based on a CLBF, as will be established in the next section. The CLBFs and controller modes are synthesized using grammar-guided genetic programming (introduced in Section 4) and verified by means of an SMT solver. The overall algorithm is described in Section 5.

3 Control strategy

In this section we discuss the used control strategy and establish how it solves problem 2.1 by means of Theorem 1 and Corollary 1.

3.1 Control Lyapunov Barrier Function

Consider a set of controller modes with index set $Q\!\subset\!\mathbb{Z}_{\geq 0}$:

[TABLE]

Given the system (1), an initial state $x=\xi(t_{k})$, let us denote the (over-approximated) reachable set for $t\in[t_{k},t_{k}+h]$ under a controller mode $q$ as $R_{q}(x)$ s.t. given a $q$, $\forall t\in[t_{k},t_{k}+h]:\xi(t)\in R_{q}(\xi(t_{k}))$. The construction of $R_{q}$ is discussed in Section 5.1. We consider a switching controller based on a CLBF defined as follows.

Definition 3.1** (Control Lyapunov Barrier Function).**

A function $V\in\mathcal{C}^{1}(S,\mathbb{R})$ is a Control Lyapunov Barrier Function (CLBF) w.r.t. the compact sets $(S,I,G)$, $S\subseteq\mathcal{X}$, $I,G\subseteq int(S)$, system (1), and controller modes (4) if there exists a scalar $\gamma>0$ such that

[TABLE]

where $A:=\{x\in S\mid V(x)\leq 0\}$ and $\dot{V}_{q}(x,z)=\langle\nabla V(z),f(z,g_{q}(x))\rangle$.

Remark 1**.**

The choice of $\gamma$ is arbitrary, because if a solution $V^{*}$ exists for $\gamma^{*}$, there always exists a linear transformation of $V^{*}$ such that the inequalities in (5) are satisfied for any $\gamma$.

Proposition 1**.**

Given a CLBF $V$, $\exists e\in\mathbb{R}$ s.t. $e=\inf_{x\in S\backslash G}V(x)$. Furthermore, the sublevel set $L_{c}:=\{x\in S\mid V(x)\leq c\}$ is compact.

Proof.

Since $V(x)$ is continuous and $S$ is compact, $V[S]\subset\mathbb{R}$ is compact and hence $V[S\backslash G]\subseteq V[S]$ is bounded, i.e. $\exists e\in\mathbb{R}$ s.t. $e=\inf_{x\in S\backslash G}V(x)$. Moreover, $Y:=\{y\in V[S]~{}|~{}y\leq c\}$ and its inverse image $V^{-1}[Y]=L_{c}$ are compact. ∎

3.2 Control policy

Given a CLBF $V$, we consider periodically switching controllers of the form

[TABLE]

where $t_{k+1}=t_{k}+h$, $t\in[t_{k},t_{k}+h)$.

3.3 Reach while stay

The presented controller strategy based on the CLBF enforces specification $\mathtt{CS}_{1}$, as shown in the following theorem.

Theorem 1**.**

Given a system (1), CLBF $V$ w.r.t. compact sets $(S,I,G)$ and controller (6), then (2) holds.

Proof.

For $\xi(t_{0})\in I$ it follows from (5a) and the definition of $A$ that $V(\xi(t_{0}))\in A$. From (5c) it follows that for all $\xi(t_{k})\in A\backslash G$ there exists a $q\in Q$ such that $\forall t\in[t_{k},t_{k}+h]:\dot{V}_{q}(\xi(t_{k}),\xi(t))\leq-\gamma$. Selecting such a mode using controller (6), applying the comparison theorem (see e.g. [17]), and using $\forall x\in A$, $V(x)\leq 0$, it follows that $\forall k\in\mathbb{Z}_{\geq 0}$, $\forall t\in[t_{k},t_{k}+h]$, $\forall\xi(t_{k})\in A\backslash G$: $V(\xi(t))\leq V(\xi(t_{k}))-\gamma h\leq-\gamma h$ Therefore, $\xi(t_{k})\in A\backslash G$ implies $\forall t\in[t_{k},t_{k}+h]$, $V(\xi(t))$ will decrease and thus cannot reach $\partial S$, as from (5b) we have $\forall x\in\partial S:V(x)>0$. Since from proposition 1 it follows that $V(x)$ on $A\backslash G\subseteq S\backslash G$ is lower bounded, $V(\xi(t))$ will decrease until in finite time $\xi(t)$ leaves $A\backslash G$ and can only enter $G$, therefore (2) holds. ∎

3.4 Reach and stay while stay

The conditions in (5) are not sufficient for forward invariance of (a subset of) the goal set, as they do not impose that $\forall x\in\partial G,V(x)<\inf_{y\in S\backslash G}V(y)$. Therefore some trajectories starting in $\partial G$ might enter $S\backslash G$ before entering $G$ again. The following corollary establishes sufficient conditions for specification $\mathtt{CS}_{2}$.

Corollary 1**.**

Given a system (1), CLBF $V$ w.r.t. compact sets $(S,I,G)$, and a controller (6), if $\exists\beta\in\mathbb{R}$ such that

[TABLE]

where $B:=\{x\in S\mid V(x)\leq\beta\}$, then (3) holds.

Proof.

From Theorem 1 we have that there exists a time $t_{K}\geq t_{0}$ such that $\xi(t_{K})\in G$. Analogous to the proof of Theorem 1, from (7b) it follows that $\forall\xi(t_{K})\in G$, $\xi(t)$ with $t\geq t_{K}$ enters in finite time $G\cap B$. From the definition of $B$ and Proposition 1 it follows that $B$ is compact and thus $G\cap B$ is compact. From (7b) and controller (6) we have that $\forall x\in\partial(G\cap B),z\in R_{q}(x):\dot{V}_{q}(x,z)\leq-\gamma$. Combining this with (7a), we have that all states $\xi(t)\in\partial(G\cap B)$ cannot reach $\partial G$ and $V(\xi(t))$ decreases, thus these trajectories will remain within $G\cap B$. Therefore it follows that $G\cap B\subseteq G$ is forward invariant. As $G\subseteq int(S)$, we have that (3) holds. ∎

Remark 2**.**

Comparing the CLBF for RSWS to the CLBF in [15], in this work the condition on the derivative of $V$ is only imposed for the sublevel set $A\subset S$, rather than the entire safe set $S$. Secondly, the CLBF in this work involves only 2 parameters $y$ and $\beta$, as opposed to 5 in [15].

4 Grammar-guided genetic programming

Genetic programming (GP) is an evolutionary algorithm capable of synthesizing entire functions, in our case a CLBF and controller modes, that minimize a cost function, without pre-specifying a fixed structure [13]. The algorithm is initialized with a random population of candidate solutions (individuals). Each individual is assigned a fitness using the fitness function, which reflects how well the design goal is satisfied. Individuals are then selected based on this fitness to undergo genetic operations. The resulting individuals form a new generation. This cycle is repeated for many generations under the expectation that the average fitness increases, until a solution is found or a maximum number of generations is met.

In GP, solutions (the phenotypes) are encoded in a certain representation (the genotypes) that allows for easy modification. In this work we use a grammar-guided genetic programming (GGGP) algorithm, similar to the work of [18], which uses a tree representation that is constructed based on a Backus-Naur form (BNF) grammar [19]. The BNF grammar consists of the tuple $\{\mathcal{N},\mathcal{S},\mathcal{P},\mathcal{P}^{*}\}$, where $\mathcal{N}$ denotes the set of nonterminals, $\mathcal{S}\in\mathcal{N}$ the start symbol, $\mathcal{P}$ the set of production rules, and $\mathcal{P}^{*}$ the set of terminal production rules, which contains no recursive rules. An example of a simple grammar to construct monomials is given by $\mathcal{N}=\{\mathrm{mon,var}\},~{}\mathcal{S}=\left<\mathrm{mon}\right>$, $\mathcal{P}$ in Table 1, and $\mathcal{P}^{*}$ obtained by omitting the recursive rules from $\mathcal{P}$. Here $\left<\mathrm{mon}\right>$ denotes monomials and $\left<\mathrm{var}\right>$ scalar variables. Using $\mathcal{P}$, $\left<\mathrm{mon}\right>$ can be mapped to either $\left<\mathrm{var}\right>$ or $\left<\mathrm{var}\right>\times\left<\mathrm{mon}\right>$.

A parse tree is constructed using the BNF grammar as follows. Starting with the start symbol nonterminal $\mathcal{S}$, a random corresponding rule is chosen from the production rules $\mathcal{P}$. This rule forms a subtree that is put under the nonterminal. Subsequently, all nonterminals in the leaves of the resulting tree are similarly expanded, until all leave nodes contain no nonterminals anymore. To limit the tree depth, $\mathcal{P}^{*}$ is used if a predefined depth is reached, such that the number of recursive rules is limited. The final parse tree is transformed into the phenotype by replacing all nonterminals with their underlying subtrees, yielding a new parse tree corresponding directly to a function. Figure 1 shows a fully grown genotype synthesized using the example grammar, as well as the transformation to its phenotype.

We use the genetic operators crossover and mutation, which take the role of exploitation and exploration of genotypes respectively. The crossover operator takes two individuals and switches two random subtrees with the same nonterminal root. The mutation operator takes a single individual and replaces a subtree corresponding to a random nonterminal with a new subtree grown from that nonterminal.

As stated before, we aim to synthesize the pair $(V,\mathcal{G})$. For both the CLBF and controller modes we use a separate parse tree, which we refer to as a gene. In this case, the genotype is formed by two genes.

5 Automatic CLBF and controller synthesis

In this section the overall algorithm is described.

5.1 One-step ahead reachable set

In this work the reachable set is constructed by using Euler’s forward method and bounding the local truncation error (LTE). This yields the following analytic expression

[TABLE]

such that the over-approximated reachable set is given by

[TABLE]

with $E:=[0,h]\times\Pi_{i=1}^{n}[-\varepsilon_{i},\varepsilon_{i}]$ and

[TABLE]

While this construction can be quite conservative, it allows for relatively simple analytic expressions.

5.2 Fitness

For each inequality (5a)-(5c) (and optionally (7a)-(7b)), an independent fitness value is constructed, consisting of a sample-based and an SMT solver-based fitness value. The former gives a measure of how much the inequalities in (10) are violated, whereas the SMT solver is used to provide a formal guarantee on whether these inequalities are satisfied. In this work we use the SMT solver dReal [20], which is able to verify nonlinear inequalities over the reals. Furthermore, in case a formula is not satisfied, the SMT solver can be used to provide a counterexample, which can again be used for the sample-based verification.

Inequalities (5a)-(5c) and (7a)-(7b) can be rewritten as 111By taking $C_{1}=I$, $C_{2}=\partial S$, $C_{3}=S\backslash G\times\Pi_{q\in Q}E$, $C_{4}=\partial G$, $C_{5}=G\times\Pi_{q\in Q}E$, $\phi_{1}(s)=-V(s)$, $\phi_{2}(s)=V(s)-c$, $\phi_{3}(s)=\chi_{A}(s_{x})(-\gamma-\min_{q\in Q}(\dot{V}(s_{x},r_{q}(s_{x},s_{\tau,q},s_{e,q}))$, $\phi_{4}(s)=V(s)-\beta-c$, $\phi_{5}(s)=\chi_{B}(s_{x})(-\gamma-\min_{q\in Q}(\dot{V}(s_{x},r_{q}(s_{x},s_{\tau,q},s_{e,q})))$, where for inequalities 3 and 5 the state $s$ is partitioned as $s=(s_{x},s_{\tau,1},s_{e,1},\dots,s_{\tau,\bar{q}},s_{e,\bar{q}})$ with $s_{x}\in\mathbb{R}^{n}$, $s_{\tau,q}\in\mathbb{R}$, $s_{e,q}\in\mathbb{R}^{n}$, $\bar{q}=|Q|$, $c$ is an arbitrary small real number to make the strict inequality non-strict, and $\chi_{D}(s)$ denotes a membership function of set $D$, i.e. $\chi_{D}(s)=1$ if $s\in D$ and zero otherwise.

[TABLE]

Given a finite set $C_{i,\mathrm{samp}}=\{x_{1},\dots,x_{n}\}$, with $C_{i,\mathrm{samp}}\subset C_{i}$, the sample-based fitness is based on an error measure w.r.t. $\phi_{i}$ defined as

[TABLE]

Using this measure, the sampled-based fitness is given by

[TABLE]

The SMT-based fitness $f_{\mathrm{SMT},\phi_{i}}$ is $1$ if it follows from dReal that the inequality is satisfied and [math] otherwise. This fitness value is only computed if an individual satisfies $f_{\mathrm{samp},\phi_{i}}=1$ for all $i\in\{1,2,3\}$ (or $i\in\{1,...,5\}$). Otherwise, $f_{\mathrm{SMT},\phi_{i}}$ are set to 0 for all conditions.

To prioritize finding a $V$ that satisfies (5a) and (5b) before checking the condition on its derivative $\dot{V}$ in (5c), (and similarly for the additional conditions (7a)-(7b)), we use the weights $w_{i}=\lfloor w_{i-1}f_{\mathrm{samp},\phi_{i-1}}\rfloor$ for $i=2,\dots 5$ and $w_{1}=1$. The overall fitness is given by

[TABLE]

Finally, to promote the selection of equivalent, but less complex individuals, candidates with the same fitness (12) are ranked according to the number of their parameters. If this is still not decisive, they are subsequently ranked based on their lowest maximum parameter.

5.3 Numerical optimization

In order to speed up the convergence of the fitness, each generation the parameters of the individuals are optimized using Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [21], which is an evolutionary optimization algorithm that is regarded to be robust with regard to discontinuous fitness functions. We use the variant sep-CMA-ES [22], because of its linear time and space complexity. In our grammar, we have the rule ‘const’ which creates a random constant. In every generation, these constants are then optimized using sep-CMA-ES, where their initial values are the current parameter values.

5.4 Algorithm outline

Provided a system, specification sets $(S,I,G)$, a grammar, and sample sets $C_{i,\mathrm{samp}}$, the proposed approach consists of the following steps:

A random population of individuals is generated as described in Section 4. 2. 2.

The parameters of the individuals are optimized using sep-CMA-ES based on the sample fitness. 3. 3.

The full fitness (12) is computed. 4. 4.

Counterexamples returned by the SMT solver are added to $C_{i,\mathrm{samp}}$. 5. 5.

The best individuals are copied to the next generation. 6. 6.

A new population is created by selecting individuals using tournament selection [13] and modifying them using genetic operators. 7. 7.

Steps 3 to 5 are repeated until all conditions are satisfied (i.e. $f_{\mathrm{SMT},\phi_{i}}=1$ for all $i$) or a maximum number of generations is reached.

Remark 3**.**

If the maximum number of generations is reached before an individual achieves $f_{\mathrm{SMT},\phi_{i}}=1$ for all $i$, no guarantees on the specification are provided.

Remark 4**.**

It is possible to pre-define the controller modes in $\mathcal{G}$, such that only the CLBF is synthesized.

5.5 Additional operations

To aid finding the correct bias $s$ of $V(x)$ such that (5a) is satisfied, the following biasing is performed before each fitness evaluation within CMA-ES:

[TABLE]

where $I_{\mathrm{samp}}$ denotes a subsampled set of $I$. To guide the search further, we impose the additional condition

[TABLE]

where $x_{c}$ denotes the center of the goal set.

6 Implementation

The switching law in (6) is computationally intensive to check online. By offline designing $\alpha_{q}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ for all $q\in Q$ such that

[TABLE]

allows us to replace the switching law with:

[TABLE]

Intuitively, when at a point $x$ a mode $q^{\prime}$ is not viable under the reachable set $R_{q}(x)$, the nominal system $\dot{V}_{q^{\prime}}(x,x)$ plus buffer $\alpha_{q^{\prime}}(x)$ should not minimize the set $\bigcup_{q\in Q}\dot{V}_{q}(x,x)+\alpha_{q}(x)$, such that it is not selected by the switching law.

Theorem 2**.**

Given a CLBF, if $\forall q\in Q$, $\alpha_{q}(x)$ satisfies (15) for $D=A\backslash G$, switching law (16) yields that (2) holds.

Proof.

This proof is by contradiction. By definition of the CLBF, for all $x\in A\backslash G$, there always exists a $q$ such that $\max_{z\in R_{q}(x)}\dot{V}_{q}(x,z)\leq-\gamma$. Assume that when using switching law (16), we have $\max_{z\in R_{q_{k}}(\xi(t))}\dot{V}_{q_{k}}(\xi(t),z)>-\gamma$. It then follows from (15) that $\min_{q\in Q}(\dot{V}_{q}(x,x)+\alpha_{q}(x))<\dot{V}_{q_{k}}(x,x)+\alpha_{q_{k}}(x)$. This directly contradicts the switching law $q_{k}=\min_{q\in Q}(\dot{V}_{q}(x,x)+\alpha_{q}(x))$. Hence switching law (16) can only select a $q_{k}$ such that $\max_{z\in R_{q_{k}}(\xi(t_{k}))}\dot{V}_{q_{k}}(\xi(t_{k}),z)\leq-\gamma$, which is guaranteed to exist by the design of the CLBF. The remainder of the proof is analogous to the proof of Theorem (1). ∎

Corollary 2**.**

Given a CLBF satisfying (7), if $\forall q$, $\alpha_{q}(x)$ satisfies (15) for $D=A\backslash int(B)$, using switching law (16) yields that (3) holds.

Proof.

This proof is analogous to the proof of Theorem 2 and Corollary 1. ∎

The functions $\alpha_{q}(x)$ can be designed and verified offline using again an SMT solver.

7 Case studies

In this section the effectiveness of the approach is demonstrated for a simple linear system, polynomial systems of second and third order (see Table 2), and two nonpolynomial systems (see Table 3). The systems and specifications are adopted from [23] and references therein, with the exception of the Pendulum system, adopted from [15].

For these case studies, we fixed the control mode vector field $\mathcal{G}$ and synthesized controllers for the reach-while-stay specification $\mathtt{CS}_{1}$. Across all these case studies, we used a population of 16 individuals, a maximum of 50 generations, and a maximum of 30 generations within CMA-ES. The mutation and crossover rates were both chosen to be 0.5. The number of test samples and maximum number of additional counterexamples were set to 100 and 300 respectively. For the counterexamples, a first-in-first-out principle was used. The (arbitrary) $\gamma$ of the CLBF was set to $\gamma=0.1$ and the precision parameter of dReal set to $\delta=0.001$. The values of $\varepsilon_{i}$ are obtained using bisection and the SMT solver. The GGGP algorithm and CMA-ES are implemented in Mathematica, running on an Intel Xeon CPU E5-1660 v3 3.00GHz using 8 CPU cores.

The used grammar is defined by $\mathcal{S}_{V}=\left<\mathrm{const}\right>+\left<\mathrm{expr}\right>$, $\mathcal{N}$ and $\mathcal{P}$ as shown in Table 4, and $\mathcal{P}^{*}$ is obtained by removing all recursive rules from $\mathcal{P}$. While this grammar restricts to polynomial CLBFs, the proposed approach can also be used for nonpolynomial CLBFs. Finally, the maximum recursive rule depth was set to be 7.

To show repeatability, the synthesis was repeated 10 times for each benchmark. Statistics on the number of generations and the total synthesis time are shown in Table 2 and 3. With the exception of the third-order polynomial system, in all 10 runs a solution was found for each benchmark. For the third-order system only a single run did not find a solution within 50 generations. Given the found solutions, we used again bisection and the SMT solver to find a $\beta$ such that the conditions in Corollary 1 hold. We found a $\beta$ such that the conditions in Corollary 1 hold for: 4 solutions of the linear system, 1 of the 2nd-order system, 8 of the 3rd-order system, 0 of the pendulum system and 6 of the pendulum on cart system. Hence for these solutions the stronger specification $\mathtt{CS}_{2}$ is guaranteed.

One of the found solutions for the pendulum system is

[TABLE]

where $x_{1}^{\prime}=(0.75+x_{1})$. We manually designed $\alpha(x)=[\alpha_{1}(x),\dots\alpha_{5}(x)]^{T}$ to be $[100,0,0,0,500]$ such that (15) holds. Figure 2 shows the phase plot of the closed-loop system for $\xi(t_{0})\in\{(-\pi,10),(-2,-5),(1.5,0),(\pi,10)$}. It can be seen that indeed all trajectories satisfy $\mathtt{CS}_{1}$.

For this solution, we could not find a $\beta$ such that Corollary 1 holds. Nevertheless, by increasing the goal set to $G$ to $[-1,-0.5]\times[-1.5,1.5]$, it can be verified that for $\beta=-4012.3$ the conditions in Corollary 1 hold.

The used sampling times $h$ are significant larger than the minimal dwell-times reported in [15] and [16]. For example, for the pendulum on cart system benchmark we used $h=0.001$ seconds, whereas [16] reports a theoretical minimum dwell-time of $2\cdot 10^{-6}$ seconds.

7.1 Evolving $\mathcal{G}$

Let us reconsider the pendulum on a cart from Table 3 and specification $\mathtt{CS}_{1}$, but without pre-specifying $\mathcal{G}$. We saturate the input with

[TABLE]

A separate gene for the controller modes $\mathcal{G}$ is used with start symbol $\mathcal{S}_{G}=\left<\mathrm{G}\right>$ and the product rules in Table 4. The results for 10 runs are shown in Table 5. Comparing Table 3 with 5 we observe a comparable number of generations required to find a solution, although a longer computation time per generation is observed. However, the benefit is that no discretization of the input space is required. One of the found solutions is given by

[TABLE]

Note that $\mathcal{G}$ consists of only a single mode, hence no switching law is required when implementing this controller. Finally, for $\beta=-19.5313$, $V$ satisfies (7), hence using this controller also guarantees $\mathtt{CS}_{2}$.

8 Discussion

This paper presented a method for automatic synthesis of a periodically switched state feedback controller for nonlinear sampled-data systems with reachability and safety specifications. Preliminary results have been shown for several nonlinear systems up to the third-order. It was shown that the framework was able to synthesize CLBFs given pre-defined controller modes, but was also capable of synthesizing controller modes automatically, eliminating the need to discretize the input space beforehand. Moreover, it is possible to find a single controller automatically, removing the need for the switching law entirely.

For almost all benchmarks runs, solutions were found within 50 generations. Nevertheless, a drawback of the proposed methodology is that there is no guarantee that a solution will be found within a number of generations.

A straightforward improvement to obtain less conservative sampling times is by using a less conservative over-approximation of the reachable set, e.g. by using higher order Taylor series approximations or using local bounds rather than for the entire domain.

Finally, the functions $\alpha_{q}(x)$ that simplify the switching condition are currently synthesized by hand. In future work, the aim is to automate this synthesis as well, for example by again using the combination of GP with SMT solvers.