Automated Formal Synthesis of Digital Controllers for State-Space Physical Plants

TL;DR

This paper introduces an automated method for synthesizing safe digital controllers for linear time-invariant physical plants, combining CEGIS and reachability analysis to ensure stability and safety.

Contribution

It presents a novel, sound, and automated approach that integrates counterexample guided synthesis with reachability analysis for controller synthesis.

Findings

Successfully synthesizes controllers for complex physical models

Ensures system safety and stability through formal verification

Demonstrates effectiveness on digital control benchmarks

Abstract

We present a sound and automated approach to synthesize safe digital feedback controllers for physical plants represented as linear, time invariant models. Models are given as dynamical equations with inputs, evolving over a continuous state space and accounting for errors due to the digitalization of signals by the controller. Our approach has two stages, leveraging counterexample guided inductive synthesis (CEGIS) and reachability analysis. CEGIS synthesizes a static feedback controller that stabilizes the system under restrictions given by the safety of the reach space. Safety is verified either via BMC or abstract acceleration; if the verification step fails, we refine the controller by generalizing the counterexample. We synthesize stable and safe controllers for intricate physical plant models from the digital control literature.

Figures

Tables

Table 1. Table 1: Experimental results

#	Benchmark	Order	Multi-staged		Abstraction
			${ℱ}_{⟨I_p,F_p⟩}$	Time	${ℱ}_{⟨I_p,F_p⟩}$	Time
1	Cruise Control	1	8,16	8.40 s	16,16	2.17 s
2	DC Motor	2	8,16	9.45 s	20,20	2.06 s
3	Helicopter	3	✗	✗	16,16	1.37 s
4	Inverted Pendulum	4	8,16	9.65 s	16,16	0.56 s
5	Magnetic Pointer	2	✗	✗	28,28	44.14 s
6	Magnetic Suspension	2	12,20	10.41 s	16,16	0.61 s
7	Pendulum	2	8,16	14.02 s	16,16	0.60 s
8	Suspension	2	12,20	73.66 s	✗	✗
9	Tape Driver	3	8,16	10.10 s	16,16	68.24 s
10	Satellite	2	8,16	9.43 s	16,16	0.67 s

Equations

$$\dot{x}(t)=Ax(t)+Bu(t),x\in {\mathbb{R}}^n,u\in {\mathbb{R}}^m,A\in {\mathbb{R}}^{n\times n},B\in {\mathbb{R}}^{n\times m},$$

$$x_{k+1}=A_d{x}_k+B_d{u}_k$$

$$u_k=r_k-K{x}_k.$$

$$x_{k+1}=(A_d-B_{d}K)x_k+B_d{r}_k.$$

$${\varphi }_{safety}⟺\forall k\ge 0.\underset{i=1}{\overset{n}{\bigwedge }}\underline{x_i}\le x_{i,k}\le \overline{x_i},$$

$$u_k$$

$$x_{k+1}$$

$$x_k$$

$$N=\left\{{\nu }_1+{\nu }_2+{\nu }_3:{\nu }_1\in \left[-\frac{q_1}{2}\frac{q_1}{2}\right]\wedge {\nu }_2\in \left[-\frac{q_2}{2}\frac{q_2}{2}\right]\wedge {\nu }_3\in \left[-q_3{q}_3\right]\right\}$$

$$({\varphi }_{init}\wedge {\varphi }_{input}\wedge u_k=-K{x}_k)\Rightarrow {\varphi }_{init}^K$$

$$x=$$

$${\hat{X}}^{\#}=\mathcal{A}{X}_0+{\mathcal{B}}_{n}N,X_0=\left\{x:x\models {\varphi }_{init}\right\},$$

$$A_{d}=\left[\begin{array}[]{ccc}2.6207&-1.1793&0.65705\\ 2&0&0\\ 0&0.5&0\end{array}\right],\quad B_{d}=\left[\begin{array}[]{c}8\\ 0\\ 0\end{array}\right]$$

$$x_{k+1}$$

$$x_1$$

$$x_2$$

$$x_n$$

$$=(A_d-B_{d}K{)}^n{x}_0+\sum _{i=0}^{i=n-1}(A_d-B_{d}K{)}^i{B}_{d}k\delta .$$

$$0+(I-A_d+B_{d}K{)}^{-1}{B}_{d}k\delta .$$

$$\mathcal{F}⟨I_c,F_c⟩(y(t))=y_k:y_k\in \left[y(t)-\frac{q_c}{2}y(t)+\frac{q_c}{2}\right].$$

$$x_{k+1}=A_d{x}_k+B_d(u_k+{{\nu }_2}_k),u_k=-K{x}_k+{{\nu }_1}_k,$$

$$x_{k+1}$$

Full text

1 Introduction

Linear Time Invariant (LTI) models represent a broad class of dynamical systems with significant impact in numerous application areas such as life sciences, robotics, and engineering [4, 13]. The synthesis of controllers for LTI models is well understood, however the use of digital control architectures adds new challenges due to the effects of finite-precision arithmetic, time discretization, and quantization noise, which is typically introduced by Analogue-to-Digital (ADC) and Digital-to-Analogue (DAC) conversion. While research on digital control is well developed [4], automated and sound control synthesis is challenging, particularly when the synthesis objective goes beyond classical stability. There are recent methods for verifying reachability properties of a given controller [14]. However, these methods have not been generalized to control synthesis. Note that a synthesis algorithm that guarantees stability does not ensure safety: the system might transitively visit an unsafe state resulting in unrecoverable failure.

We propose a novel algorithm for the synthesis of control algorithms for LTI models that are guaranteed to be safe, considering both the continuous dynamics of the plant and the finite-precision discrete dynamics of the controller, as well as the hybrid elements that connect them. We account for the presence of errors originating from a number of sources: quantisation errors in ADC and DAC, representation errors (from the discretization introduced by finite-precision arithmetic), and roundoff and saturation errors in the verification process (from finite-precision operations). Due to the complexity of such systems, we focus on linear models with known implementation features (e.g., number of bits, fixed-point arithmetic). We expect a safety requirement given as a reachability property. Safety requirements are frequently overlooked in conventional feedback control synthesis, but play an important role in systems engineering.

We give two alternative approaches for synthesizing digital controllers for state-space physical plants, both based on CounterExample Guided Inductive Synthesis (CEGIS) [30]. We prove their soundness by quantifying errors caused by digitization and quantization effects that arise when the digital controller interacts with the continuous plant.

The first approach uses a naïve technique that starts by devising a digital controller that stabilizes the system while remaining safe for a pre-selected time horizon and a single initial state; then, it verifies unbounded safety by unfolding the dynamics of the system, considering the full hyper-cube of initial states, and checking a completeness threshold [17], i.e., the number of iterations required to sufficiently unwind the closed-loop state-space system such that the boundaries are not violated for any larger number of iterations. As it requires unfolding up to the completeness threshold, this approach is computationally expensive.

Instead of unfolding the dynamics, the second approach employs abstract acceleration [7] to evaluate all possible progressions of the system simultaneously. Additionally, the second approach uses abstraction refinement, enabling us to always start with a very simple description regardless of the dynamics complexity, and only expand to more complex models when a solution cannot be found.

We provide experimental results showing that both approaches are able to efficiently synthesize safe controllers for a set of intricate physical plant models taken from the digital control literature: the median run-time for our benchmark set is $7.9$ seconds, and most controllers can be synthesized in less than $17.2$ seconds. We further show that, in a direct comparison, the abstraction-based approach (i.e., the second approach) lowers the median run-time of our benchmarks by a factor of seven over the first approach based on the unfolding of the dynamics.

Contributions

We compute state-feedback controllers that guarantee a given safety property. Existing methods for controller synthesis rely on transfer function representations, which are inadequate to prove safety requirements. 2. 2.

We provide two novel algorithms: the first, naïve one, relies on an unfolding of the dynamics up to a completeness threshold, while the second one is abstraction-based and leverages abstraction refinement and acceleration to improve scalability while retaining soundness. Both approaches provide sound synthesis of state-feedback systems and consider the various sources of imprecision in the implementation of the control algorithm and in the modeling of the plant. 3. 3.

We develop a model for different sources of quantization errors and their effect on reachability properties. We give bounds that ensure the safety of our controllers in a hybrid continuous-digital domain.

2 Related Work

CEGIS -

Program synthesis is the problem of computing correct-by-design programs from high-level specifications. Algorithms for this problem have made substantial progress in recent years, for instance [16] to inductively synthesize invariants for the generation of desired programs.

Program synthesizers are an ideal fit for the synthesis of digital controllers, since the semantics of programs capture the effects of finite-precision arithmetic precisely. In [27], the authors use CEGIS for the synthesis of switching controllers for stabilizing continuous-time plants with polynomial dynamics. The work extends to affine systems, but is limited by the capacity of the state-of-the-art SMT solvers for solving linear arithmetic. Since this approach uses switching models instead of linear dynamics for the digital controller, it avoids problems related to finite precision arithmetic, but potentially suffers from state-space explosion. Moreover, in [28] the same authors use a CEGIS-based approach for synthesizing continuous-time switching controllers that guarantee reach-while-stay properties of closed-loop systems, i.e., properties that specify a set of goal states and safe states (constrained reachability). This solution is based on synthesizing control Lyapunov functions for switched systems that yield switching controllers with a guaranteed minimum dwell time in each mode. However, both approaches are unsuitable for the kind of control we seek to synthesize.

The work in [2] synthesizes stabilizing controllers for continuous plants given as transfer functions by exploiting bit-accurate verification of software-implemented digital controllers [5]. While this work also uses CEGIS, the approach is restricted to digital controllers for stable closed-loop systems given as transfer function models: this results in a static check on their coefficients. By contrast, in the current paper we consider a state-space representation of the physical system, which requires ensuring the specification over actual traces of the model, alongside the numerical soundness required by the effects of discretisation and finite-precision errors. A state-space model has known advantages over the transfer function representation [13]: it naturally generalizes to multivariate systems (i.e., with multiple inputs and outputs); and it allows synthesis of control systems with guarantees on the internal dynamics, e.g., to synthesize controllers that make the closed-loop system safe. Our work focuses on the safety of internal states, which is usually overlooked in the literature. Moreover, our work integrates an abstraction/refinement (CEGAR) step inside the main CEGIS loop.

The tool Pessoa [21] synthesizes correct-by-design embedded control software in a Matlab toolbox. It is based on the abstraction of a physical system to an equivalent finite-state machine and on the computation of reachability properties thereon. Based on this safety specification, Pessoa can synthesize embedded controller software for a range of properties. The embedded controller software can be more complicated than the state-feedback control we synthesize, and the properties available cover more detail. However, relying on state-space discretization Pessoa is likely to incur in scalability limitations. Along this research line, [3, 20] studies the synthesis of digital controllers for continuous dynamics, and [34] extends the approach to the recent setup of Network Control Systems.

Discretization Effects -

The classical approach to control synthesis has often disregarded digitalization effects, whereas more recently modern techniques have focused on different aspects of discretization, including delayed response [10] and finite word length (FWL) semantics, with the goal either to verify (e.g., [9]) or to optimize (e.g., [24]) given implementations.

There are two different problems that arise from FWL semantics. The first is the error in the dynamics caused by the inability to represent the exact state of the physical system, while the second relates to rounding and saturation errors during computation. In [12], a stability measure based on the error of the digital dynamics ensures that the deviation introduced by FWL does not make the digital system unstable. A more recent approach [33] uses $\mu$-calculus to directly model the digital controller so that the selected parameters are stable by design. The analyses in [29, 32] rely on an invariant computation on the discrete system dynamics using Semi-Definite Programming (SDP). While the former uses bounded-input and bounded-output (BIBO) properties to determine stability, the latter uses Lyapunov-based quadratic invariants. In both cases, the SDP solver uses floating-point arithmetic and soundness is checked by bounding the error. An alternative is [25], where the verification of given control code is performed against a known model by extracting an LTI model of the code by symbolic execution: to account for rounding errors, an upper bound is introduced in the verification phase. The work in [26] introduces invariant sets as a mechanism to bound the quantization error effect on stabilization as an invariant set that always converges toward the controllable set. Similarly, [19] evaluates the quantization error dynamics and bounds its trajectory to a known region over a finite time period. This technique works for both linear and non-linear systems.

3 Preliminaries

3.1 State-space representation of physical systems

We consider models of physical plants expressed as ordinary differential equations (ODEs), which we assume to be controllable and under full state information (i.e., we have access to all the model variables):

[TABLE]

where $t\in\mathbb{R}_{0}^{+}$, where $A$ and $B$ are matrices that fully specify the continuous plant, and with initial states set as $x(0)$. While ideally we intend to work on the continuous-time plant, in this work Eq. (1) is soundly discretized in time [11] into

[TABLE]

where $k\in\mathbb{N}$ and $x_{0}=x(0)$ is the initial state. $A_{d}$ and $B_{d}$ denote the matrices that describe the discretized plant dynamics, whereas $A$ and $B$ denote the continuous plant dynamics. We synthesize for requirements over this discrete-time domain. Later, we will address the issue of variable quantization, as introduced by the ADC/DAC conversion blocks (Fig. 1).

3.2 Controller synthesis via state feedback

Models (1) and (2) depend on external non-determinism in the form of input signals $u(t)$ and $u_{k}$, respectively. Feedback architectures can be employed to manipulate the properties and behaviors of the continuous process (the plant). We are interested in the synthesis of digital feedback control algorithms, as implemented on Field-Programmable Gate Arrays or Digital Signal Processors. The most basic feedback architecture is the state feedback one, where the control action $u_{k}$ (notice we work with the discretized signal) is computed by:

[TABLE]

Here, $K\in\mathbb{R}^{m\times n}$ is a state-feedback gain matrix, and $r_{k}$ is a reference signal (again digital). The closed-loop model then takes the form

[TABLE]

The gain matrix $K$ can be set so that the closed-loop discrete dynamics are shaped as desired, for instance according to a specific stability goal or around a specific dynamical behavior [4]. As argued later in this work, we will target more complex objectives, such as quantitative safety requirements, which are not typical in the digital control literature. Further, we will embrace the digital nature of the controller, which manipulates quantized signals as discrete quantities represented with finite precision.

3.3 Stability of closed-loop systems

In this work we employ asymptotic stability in the CEGIS loop, as an objective for guessing controllers that are later proven sound over safety requirements. Asymptotic stability is a property that amounts to convergence of the model executions to an equilibrium point, starting from any states in a neighborhood of the point (see Figure 3 for the portrait of a stable execution, converging to the origin). In the case of linear systems as in (4), considered with a zero reference signal, the equilibrium point of interest is the origin.

A discrete-time LTI system as (4) is asymptotically stable if all the roots of its characteristic polynomial (i.e., the eigenvalues of the closed-loop matrix $A_{d}-B_{d}K$) are inside the unity circle of the complex plane, i.e., their absolute values are strictly less than one [4] (this simple sufficient condition can be generalised, however this is not necessary in our work). In this paper, we express this stability specification $\phi_{\mathit{stability}}$ in terms of a check known as Jury’s criterion [11]: this is an easy algebraic formula to select the entries of matrix $K$ so that the closed-loop dynamics are shaped as desired.

3.4 Safety specifications for dynamical systems

We are not limited to the synthesis of digital stabilizing controllers – a well known task in the literature on digital control systems – but target safety requirements with an overall approach that is sound and automated. More specifically, we require that the closed-loop system (4) meets given safety specifications. A safety specification gives raise to a requirement on the states of the model, such that the feedback controller (namely the choice of the gains matrix $K$) must ensure that the state never violates the requirement. Note that a stable, closed-loop system is not necessarily a safe system: indeed, the state values may leave the safe part of the state space while they converge to the equilibrium, which is typical in the case of oscillatory dynamics. In this work, the safety property is expressed as:

[TABLE]

where $\underline{x_{i}}$ and $\overline{x_{i}}$ are lower and upper bounds for the $i$-th coordinate $x_{i}$ of state $x\in\mathbb{R}^{n}$ at the $k$-th instant, respectively. This means that the states will always be within an $n$-dimensional hyper-box.

Furthermore, it is practically relevant to consider the constraints $\phi_{\mathit{input}}$ on the input signal $u_{k}$ and $\phi_{\mathit{init}}$ on the initial states $x_{0}$, which we assume have given bounds: $\phi_{\mathit{input}}={\forall k.\underline{u}\leq u_{k}\leq\overline{u}}$, $\phi_{\mathit{init}}=\bigwedge_{i=1}^{n}\underline{x_{i,0}}\leq x_{i,0}\leq\overline{x_{i,0}}.$ For the former, this means that the control input might saturates in view of physical constraints.

3.5 Numerical representation and soundness

The models we consider have two sources of error that are due to numerical representation. The first is the numerical error introduced by the fixed-point numbers employed to model the plant, i.e., to represent the plant dynamics $A_{d}$, $B_{d}$ and $x_{k}$. The second is the quantization error introduced by the digital controller, which performs operations on fixed-point numbers. In this section we outline the notation for the fixed-point representation of numbers, and briefly describe the errors introduced. A formal discussion is in Appendix 0.B.1.

Let $\mathcal{F}_{\langle I,F\rangle}(x)$ denote a real number $x$ represented in a fixed point domain, with $I$ bits representing the integer part and $F$ bits representing the decimal part. The smallest number that can be represented in this domain is $c_{m}=2^{-F}$. Any mathematical operations performed at the precision $\mathcal{F}_{\langle I,F\rangle}(x)$ will introduce errors, for which an upper bound can be given [6].

We will use $\mathcal{F}_{\langle I_{c},F_{c}\rangle}(x)$ to denote a real number $x$ represented at the fixed-point precision of the controller, and $\mathcal{F}_{\langle I_{p},F_{p}\rangle}(x)$ to denote a real number $x$ represented at the fixed-point precision of the plant model ($I_{c}$ and $F_{c}$ are determined by the controller. We pick $I_{p}$ and $F_{p}$ for our synthesis such that $I_{p}\geq I_{c}$ and $\allowbreak F_{p}\geq F_{c}$). Thus any mathematical operations in our modelled digital controller will be in the range of $\mathcal{F}_{\langle I_{c},F_{c}\rangle}$, and all other calculations in our model will be carried out in the range of $\mathcal{F}_{\langle I_{p},F_{p}\rangle}$. The physical plant operates in the reals, which means our verification phase must also account for the numerical error and quantization errors caused by representing the physical plant at the finite precision $\mathcal{F}_{\langle I_{p},F_{p}\rangle}$.

3.5.1 Effect on safety specification and stability

Let us first consider the effect of the quantization errors on safety. Within the controller, state values are manipulated at low precision, alongside the vector multiplication $Kx$. The inputs are computed using the following equation:

[TABLE]

This induces two types of the errors detailed above: first, the truncation error due to representing $x_{k}$ as $\mathcal{F}{\langle I_{c},F_{c}\rangle}(x_{k})$; and second, the rounding error introduced by the multiplication operation. We represent these errors as non-deterministic additive noise.

An additional error is due to the representation of the plant dynamics, namely

[TABLE]

We address this error by use of interval arithmetic [22] in the verification phase.

Previous studies [18] show that the FWL affects the poles and zeros positions, degrading the closed-loop dynamics, causing steady-state errors (see Appendix 0.B for details) and eventually de-stabilizing the system [5]. However, since in this paper we require stability only as a precursor to safety, it is sufficient to check that the (perturbed, noisy) model converges to a neighborhood of the equilibrium within the safe set (see Appendix 0.A.1).

In the following, we shall disregard these steady-state errors (caused by FWL effects) when stability is ensured by synthesis, and then verify its safety accounting for the finite-precision errors.

4 CEGIS of Safe Controllers for LTI Systems

In this section, we describe our technique for synthesizing safe digital feedback controllers using CEGIS. For this purpose, we first provide the synthesizer’s general architecture, followed by describing our two approaches to synthesizing safe controllers: the first one is a baseline approach that relies on a naïve unfolding of the transition relation, whereas the second uses abstraction to evaluate all possible executions of the system.

4.1 General architecture of the program synthesizer

The input specification provided to the program synthesizer is of the form $\exists P.\,\forall a.\,\sigma(a,P)$, where $P$ ranges over functions (where a function is represented by the program computing it), $a$ ranges over ground terms, and $\sigma$ is a quantifier-free formula. We interpret the ground terms over some finite domain $\mathcal{D}$. The design of our synthesizer consists of two phases, an inductive synthesis phase and a validation phase, which interact via a finite set of test vectors inputs that is updated incrementally. Given the aforementioned specification $\sigma$, the inductive synthesis procedure tries to find an existential witness $P$ satisfying the specification $\sigma(a,P)$ for all $a$ in inputs (as opposed to all $a\in\mathcal{D}$). If the synthesis phase succeeds in finding a witness $P$, this witness is a candidate solution to the full synthesis formula. We pass this candidate solution to the validation phase, which checks whether it is a full solution (i.e., $P$ satisfies the specification $\sigma(a,P)$ for all $a\in\mathcal{D}$). If this is the case, then the algorithm terminates. Otherwise, additional information is provided to the inductive synthesis phase in the form of a new counterexample that is added to the inputs set and the loop iterates again. More details about the general architecture of the synthesizer can be found in [8].

4.2 Synthesis problem: statement (recap) and connection to program synthesis

At this point, we recall the the synthesis problem that we solve in this work: we seek a digital feedback controller $K$ (see Eq. 3) that makes the closed-loop plant model safe for initial state $x_{0}$, reference signal $r_{k}$ and input $u_{k}$ as defined in Sec. 3.4. We consider non-deterministic initial states within a specified range, the reference signal to be set to zero, saturation on the inputs, and account for digitization and quantization errors introduced by the controller.

When mapping back to the notation used for describing the general architecture of the program synthesizer, the controller $K$ denotes $P$, $(x_{0},u_{k})$ represents $a$ and $\phi_{\mathit{stability}}\wedge\phi_{\mathit{input}}\wedge\phi_{\mathit{init}}\wedge\phi_{\mathit{safety}}$ denotes the specification $\sigma$.

4.3 Naïve Approach: CEGIS with multi-staged verification

An overview of the algorithm for controller synthesis is given in Fig. 2. One important observation is that we verify and synthesize a controller over $k$ time steps. We then compute a completeness threshold $\overline{k}$ [17] for this controller, and verify correctness for $\overline{k}$ time steps. Essentially, $\overline{k}$ is the number of iterations required to sufficiently unwind the closed-loop state-space system, which ensures that the boundaries are not violated for any other $k{>}\overline{k}$.

Theorem 4.1

There exists a finite $\overline{k}$ such that it is sufficient to unwind the closed-loop state-space system up to $\overline{k}$ in order to ensure that $\phi_{\mathit{safety}}$ holds.

Proof

A stable control system is known to have converging dynamics. Assume the closed-loop matrix eigenvalues are not repeated (which is sensible to do, since we select them). The distance of the trajectory from the reference point (origin) decreases over time within subspaces related to real-valued eigenvalues; however, this is not the case in general when dealing with complex eigenvalues. Consider the closed-loop matrix that updates the states in every discrete time step, and select the eigenvalue $\vartheta$ with the smallest (non-trivial) imaginary value. Between every pair of consecutive time steps $k\,T_{s}$ and $(k+1)\,T_{s}$, the dynamics projected on the corresponding eigenspace rotate $\vartheta T_{s}$ radians. Thus, taking $\overline{k}$ as the ceiling of $\frac{2\pi}{\vartheta T_{s}}$, after $k{\geq}\overline{k}$ steps we have completed a full rotation, which results in a point closer to the origin. The synthesized $\overline{k}$ is the completeness threshold. ∎

Next, we describe the different phases in Fig. 2 (blocks 1 to 4) in detail.

The inductive synthesis phase (synthesize) uses BMC to compute a candidate solution $K$ that satisfies both the stability criteria (Sec. 3.3) and the safety specification (Sec. 3.4). To synthesize a controller that satisfies the stability criteria, we require that a computed polynomial satisfies Jury’s criterion [11]. The details of this calculation can be found in the Appendix.

Regarding the second requirement, we synthesize a safe controller by unfolding the transition system $k$ steps and by picking a controller $K$ and a single initial state, such that the states at each step do not violate the safety criteria. That is, we ask the bounded model checker if there exists a $K$ that is safe for at least one $x_{0}$ in our set of all possible initial states. This is sound if the current $k$ is greater than the completeness threshold. We also assume some precision $\langle I_{p},F_{p}\rangle$ for the plant and a sampling rate. The checks that these assumptions hold are performed by subsequent verify stages.

The first verify stage, safety, checks that the candidate solution $K$, which we synthesized to be safe for at least one initial state, is safe for all possible initial states, i.e., does not reach an unsafe state within $k$ steps where we assume $k$ to be under the completeness threshold. After unfolding the transition system corresponding to the previously synthesized controller $k$ steps, we check that the safety specification holds for any initial state. This is shown in Alg.  1. 3. 3.

The second verify stage, precision, restores soundness with respect to the plant’s precision by using interval arithmetic [22] to validate the operations performed by the previous stage. 4. 4.

The third verify stage, complete, checks that the current $k$ is large enough to ensure safety for any $k^{\prime}{>}k$. Here, we compute the completeness threshold $\overline{k}$ for the current candidate controller $K$ and check that $k{\geq}\overline{k}$. This is done according to the argument given above and illustrated in Fig. 3.

Checking that the safety specification holds for any initial state can be computationally expensive if the bounds on the allowed initial states are large.

Theorem 4.2

If a controller is safe for each of the corner cases of our hypercube of allowed initial states, it is safe for any initial state in the hypercube.

Thus we only need to check $2^{n}$ initial states, where $n$ is the dimension of the state space (number of continuous variables).

Proof

Consider the set of initial states, $X_{0}$, which we assume to be convex since it is a hypercube. Name $v_{i}$ its vertexes, where $i=1,\ldots,2^{n}$. Thus any point $x\in X_{0}$ can be expressed by convexity as $x=\sum_{i=1}^{2^{n}}\alpha_{i}v_{i}$, where $\sum_{i=1}^{2^{n}}\alpha_{i}=1$. Then if $x_{0}=x$, we obtain

[TABLE]

where $x_{k}^{i}$ denotes the trajectories obtained from the single vertex $v_{i}$. We conclude that any $k$-step trajectory is encompassed, within a convex set, by those generated from the vertices. ∎

4.3.1 Illustrative Example

We illustrate our approach with an example, extracted from [13]. Since we have not learned any information about the system yet, we pick an arbitrary candidate solution (we always choose $K=[0\ 0\ 0]^{T}$ in our experiments to simplify reproduction), and a precision of $I_{p}=13$, $F_{p}=3$. In the first verify stage, the safety check finds the counterexample $x_{0}=[-0.5\ 0.5\ 0.5]$. After adding the new counterexample to its sets of inputs, synthesize finds the candidate solution $K=[0\,0\,0.00048828125]^{T}$, which prompts the safety verifier to return $x_{0}=[-0.5\,-0.5\,-0.5]$ as the new counterexample.

In the subsequent iteration, the synthesizer is unable to find further suitable candidates and it returns UNSAT, meaning that the current precision is insufficient. Consequently, we increase the precision the plant is modelled with to $I_{p}=17$, $F_{p}=7$. We increase the precision by 8 bits each step in order to be compliant with the CBMC type API. Since the previous counterexamples were obtained at lower precision, we remove them from the set of counterexamples. Back in the synthesize phase, we re-start the process with a candidate solution with all coefficients [math]. Next, the safety verification stage provides the first counterexample at higher precision, $x_{0}=[-0.5\ 0.5\ 0.5]$ and synthesize finds $K=[0\ 0.01171875\ 0.015625]^{T}$ as a candidate that eliminates this counterexample. However, this candidate triggers the counterexample $x_{0}=[0.5\ -0.5\ -0.5]$ found again by the safety verification stage. In the next iteration, we get the candidate $K=[0\ 0\ -0.015625]$, followed by the counterexample $x_{0}=[0.5\ 0.5\ 0.5]$. Finally, synthesize finds the candidate $K=[0.01171875\ -0.013671875\ -0.013671875]^{T}$, which is validated as a final solution by all verification stages.

4.4 Abstraction-based CEGIS

The naïve approach described in Sec. 4.3 synthesizes a controller for an individual initial state and input with a bounded time horizon and, subsequently, it generalizes it to all reachable states, inputs, and time horizons during the verification phase. Essentially, this approach relies on the symbolic simulation over a bounded time horizon of individual initial states and inputs that form part of an uncountable space and tries to generalize it for an infinite space over an infinite time horizon.

Conversely, in this section, we find a controller for a continuous initial set of states and set of inputs, over an abstraction of the continuous dynamics [7] that conforms to witness proofs at specific times. Moreover, this approach uses abstraction refinement enabling us to always start with a very simple description regardless of the complexity of the overall dynamics, and only expand to more complex models when a solution cannot be found.

The CEGIS loop for this approach is illustrated in Fig. 4.

We start by doing some preprocessing:

1. (a)

Compute the characteristic polynomial of the matrix $(A_{d}-B_{d}K)$ as $P_{a}(z)=z^{n}+\sum_{i=1}^{n}{(a_{i}-k_{i})z^{n-i}}$. 2. (b)

Calculate the noise set $N$ from the quantizer resolutions and estimated round-off errors:

[TABLE]

where $q_{1}$ is the error introduced by the truncation in the ADC, $q_{2}$ is the error introduced by the DAC and $q_{3}$ is the maximum truncation and rounding error in $u_{k}=-K\cdot\mathcal{F}_{\langle I_{c},F_{c}\rangle}(x_{k})$ as discussed in Section 3.5. More details on how to model quantization as noise are given in Appendix 0.B.2. 3. (c)

Calculate a set of initial bounds on $K$, $\phi_{\mathit{init}}^{K}$, based on the input constraints

[TABLE]

Note that these bounds will be used by the synthesize phase to reduce the size of the solution space. 2. 2.

In the synthesize phase, we synthesize a candidate controller $K\in\mathbb{R}\langle I_{c},F_{c}\rangle^{n}$ that satisfies $\phi_{\mathit{stability}}\wedge\phi_{\mathit{safety}}\wedge\phi_{\mathit{init}}^{K}$ by invoking a SAT solver. If there is no candidate solution we return UNSAT and exit the loop. 3. 3.

Once we have a candidate solution, we perform a safety verification of the progression of the system from $\phi_{\mathit{init}}$ over time, $x_{k}\models\phi_{\mathit{safety}}$. In order to compute the progression of point $x_{0}$ at iteration $k$, we accelerate the dynamics of the closed-loop system and obtain:

[TABLE]

As this still requires us to verify the system for every $k$ up to infinity, we use abstract acceleration again to obtain the reach-tube, i.e., the set of all reachable states at all times given an initial set $\phi_{\mathit{init}}$:

[TABLE]

where $\mathcal{A}=\bigcup_{k=1}^{\infty}(A_{d}-B_{d}K)^{k},\mathcal{B}_{n}=\bigcup_{k=1}^{\infty}\sum_{i=0}^{k}(A_{d}-B_{d}K)^{i}B_{n}$ are abstract matrices for the closed-loop system [7], whereas the set $N$ is non-deterministically chosen.

We next evaluate $\hat{X}^{\#}\models\phi_{\mathit{safety}}$. If the verification holds we have a solution, and exit the loop. Otherwise, we find a counterexample iteration $k$ and corresponding initial point $x_{0}$ for which the property does not hold, which we use to locally refine the abstraction. When the abstraction cannot be further refined, we provide them to the abstract phase. 4. 4.

If we reach the abstract phase, it means that the candidate solution is not valid, in which case we must refine the abstraction used by the synthesizer.

1. (a)

Find the constraints that invalidate the property as a set of counterexamples for the eigenvalues, which we define as $\phi_{\Lambda}$. This is a constraint in the spectrum i.e., transfer function) of the closed loop dynamics. 2. (b)

We use $\phi_{\Lambda}$ to further constrain the characteristic polynomial $z^{n}+\sum_{i=1}^{n}(a_{i}-k_{i})z^{n-i}=\prod_{i=1}^{n}(z-\lambda_{i}):|\lambda_{i}|<1\wedge\lambda_{i}\models\phi_{\Lambda}$. These constraints correspond to specific iterations for which the system may be unsafe. 3. (c)

Pass the refined abstraction $\phi(K)$ with the new constraints and the list of iterations $k$ to the synthesize phase.

4.4.1 Illustrative Example

Let us consider the following example with discretized dynamics

[TABLE]

Using the initial state bounds $\underline{x_{0}}=-0.9$ and $\overline{x_{0}}=0.9$, the input bounds $\underline{u}=-10$ and $\overline{u}=10$, and safety specifications $\underline{x_{i}}=-0.92$ and $\overline{x_{i}}=0.92$, the synthesize phase in Fig. 4 generates an initial candidate controller K=[\begin{array}[]{ccc}0.24609375&-0.125&0.1484375\end{array}]. This candidate is chosen for its closed-loop stable dynamics, but the verify phase finds it to be unsafe and returns a list of iterations with an initial state that fails the safety specification (k,x_{0}){\in}\{(2,[\begin{array}[]{ccc}0.9&-0.9&0.9\end{array}]),\allowbreak(3,[\begin{array}[]{ccc}0.9&-0.9&-0.9\end{array}])\}. This allows the abstract phase to create a new safety specification that considers these iterations for these initial states to constrain the solution space. This refinement allows synthesize to find a new controller K=[\begin{array}[]{ccc}0.23828125&-0.17578125&0.109375\end{array}], which this time passes the verification phase, resulting in a safe system.

5 Experimental Evaluation

5.1 Description of the benchmarks

A set of state-space models for different classes of systems has been taken from the literature [13, 23, 31, 1] and employed for validating our methodology.

DC Motor Rate plants describes the angular velocity of a DC Motor, respectively. The Automotive Cruise System plant represents the speed of a motor vehicle. The Helicopter Longitudinal Motion plant provides the longitudinal motion model of a helicopter. The Inverted Pendulum plant describes a pendulum model with its center of mass above its pivot point. The Magnetic Suspension plant provides a physical model for which a given object is suspended via a magnetic field. The Magnetized Pointer plant describes a physical model employed in analogue gauges and indicators that is rotated through interaction with magnetic fields. The 1/4 Car Suspension plant presents a physical model that connects a car to its wheels and allows relative motion between the two parts. The Computer Tape Driver plant describes a system to read and write data on a storage device.

Our benchmarks are SISO models (Section 3). The Inverted Pendulum appears to be a two-output system, but it is treated as two SISO models during the experiments. All the state measurements are assumed to be available (current work targets the extension of our framework to observer-based synthesis).

All benchmarks are discretized with different sample times [11]. All experiments are performed considering $\underline{x_{i}}=-1$ and $\overline{x_{i}}=1$ and the reference inputs $r_{k}=0,\forall k>0$. We conduct the experimental evaluation on a 12-core 2.40 GHz Intel Xeon E5-2440 with 96 GB of RAM and Linux OS. We use the Linux times command to measure CPU time used for each benchmark. The runtime is limited to one hour per benchmark.

5.2 Objectives

Using the state-space models given in Section 5.1, our evaluation has the following two experimental goals:

1. EG1

(CEGIS) Show that both the multi-staged and the abstraction-based CEGIS approaches are able to generate FWL digital controllers in a reasonable amount of time. 2. EG2

(sanity check) Confirm the stability and safety of the synthesized controllers outside of our model.

5.3 Results

We provide the results in Table 1. Here Benchmark is the name of the respective benchmark, Order is the number of continuous variables, $\mathcal{F}_{\langle I_{p},F_{p}\rangle}$ is the fixed-point precision used to model the plant, while Time is the total time required to synthesize a controller for the given plant with one of the two methods. Timeouts are indicated by ✗. The precision for the controller, $\mathcal{F}_{\langle I_{c},F_{c}\rangle}$, is chosen to be $I_{c}=8$, $F_{c}=8$.

For the majority of our benchmarks, we observe that the abstraction-based back-end is faster than the basic multi-staged verification approach, and finds one solution more (9) than the multi-staged back-end (8). In direct comparison, the abstraction-based approach is on average able to find a solution in approximately 70% of the time required using the multi-staged back-end, and has a median run-time 1.4 s, which is seven times smaller than the multi-staged approach. The two back-ends complement each other in benchmark coverage and together solve all benchmarks in the set. On average our engine spent 52% in the synthesis and 48% in the verification phase.

The median run-time for our benchmark set is 9.4 s. Overall, the average synthesis time amounts to approximately 15.6 s. We consider these times short enough to be of practical use to control engineers, and thus affirm EG1.

There are a few instances for which the system fails to find a controller. For the naïve approach, the completeness threshold may be too large, thus causing a timeout. On the other hand, the abstraction-based approach may require a very precise abstraction, resulting in too many refinements and, consequently, in a timeout. Yet another source of incompleteness is the inability of the synthesize phase to use a large enough precision for the plant model.

The synthesized controllers are confirmed to be safe outside of our model representation using MATLAB, achieving EG2. A link to the full experimental environment, including scripts to reproduce the results, all benchmarks and the tool, is provided in the footnote as an Open Virtual Appliance (OVA).111www.cprover.org/DSSynth/controller-synthesis-cav-2017.tar.gz The provided experimental environment runs multiple discretisations for each benchmark, and lists the fastest as the result synthesis time.

5.4 Threats to validity

Benchmark selection: We report an assessment of both our approaches over a diverse set of real-world benchmarks. Nevertheless, this set of benchmarks is limited within the scope of this paper and the performance may not generalize to other benchmarks.

Plant precision and discretization heuristics: Our algorithm to select suitable FWL word widths to model the plant behavior increases the precision by 8 bits at each step in order to be compliant with the CBMC type API. Similarly, for discretization, we run multiple discretizations for each benchmark and retain the fastest run. This works sufficiently well for our benchmarks, but performance may suffer in some cases, for example if the completeness threshold is high.

Abstraction on other properties: The performance gain from abstract acceleration may not hold for more complex properties than safety, for instance “eventually reach and always remain in a given safe set”.

6 Conclusion

We have presented two automated approaches to synthesize digital state-feedback controllers that ensure both stability and safety over the state-space representation. The first approach relies on unfolding of the closed-loop model dynamics up to a completeness threshold, while the second one applies abstraction refinement and acceleration to increase speed whilst retaining soundness. Both approaches are novel within the control literature: they give a fully automated synthesis method that is algorithmically and numerically sound, considering various error sources in the implementation of the digital control algorithm and in the computational modeling of plant dynamics. Our experimental results show that both approaches are able to synthesize safe controllers for most benchmarks within a reasonable amount of time fully automatically. In particular, both approaches complement each other and together solve all benchmarks, which have been derived from the control literature.

Future work will focus the extension of these approaches to the continuous-time case, to models with output-based control architectures (with the use of observers), and to the consideration of more complex specifications.

Appendix 0.A Stability of Closed-loop Models

0.A.1 Stability of closed-loop models with fixed-point controller error

The proof of Jury’s criterion [11] relies on the fact that the relationship between states and next states is defined by $x_{k+1}=(A_{d}-B_{d}K)x_{k}$, all computed at infinite precision. When we employ a FWL digital controller, the operation becomes:

[TABLE]

where $\delta$ is the maximum error that can be introduced by the FWL controller in one step, i.e., by reading the states values once and multiplying by $K$ once. We derive the closed form expression for $x_{n}$ as follows:

[TABLE]

The definition of asymptotic stability is that the system converges to a reference signal, in this case we use no reference signal so an asymptotically stable system will converge to the origin. We know that the original system with an infinite-precision controller is stable, because we have synthesized it to meet Jury’s criterion. Hence, $(A_{d}-B_{d}K)^{n}x_{0}$ must converge to zero.

The power series of matrices converges [15] iff the eigenvalues of the matrix are less than 1 as follows: $\sum_{i=0}^{\infty}T^{i}=(I-T)^{-1}$, where $I$ is the identity matrix and $T$ is a square matrix. Thus, our system will converge to the value

[TABLE]

As a result, if the value $(I-A_{d}+B_{d}K)^{-1}B_{d}k\delta$ is within the safe space, then the synthesized fixed-point controller results in a safe closed-loop model. The convergence to a finite value, however, will not make it asymptomatically stable.

Appendix 0.B Errors in LTI models

0.B.1 Errors due to numerical representation

We have used $\mathcal{F}_{\langle I,F\rangle}(x)$ denote a real number $x$ represented in a fixed point domain, with $I$ bits representing the integer part and $F$ bits representing the decimal part. The smallest number that can be represented in this domain is $c_{m}=2^{-F}$. The following approximation errors will arise in mathematical operations and representation:

Truncation: Let $x$ be a real number, and $\mathcal{F}_{\langle I,F\rangle}(x)$ be the same number represented in a fixed-point domain as above. Then $\mathcal{F}_{\langle I,F\rangle}(x)=x-\delta_{T}$ where the error $\delta_{T}=x\ \%_{c_{m}}\ \tilde{x}$, and $\%_{c_{m}}$ is the modulus operation performed on the last bit. Thus, $\delta_{T}$ is the truncation error and it will propagate across operations. 2. 2.

Rounding: The following errors appear in basic operations. Let $c_{1},c_{2}$ and $c_{3}$ be real numbers, and $\delta_{T1}$ and $\delta_{T2}$ be the truncation errors caused by representing $c_{1}$ and $c_{2}$ in the fixed-point domain as above.

1. (a)

Addition/Subtraction: these operations only propagate errors coming from truncation of the operands, namely $\mathcal{F}_{\langle I,F\rangle}(c_{1})\pm\mathcal{F}_{\langle I,F\rangle}(c_{2})=c_{3}+\delta_{3}$ with $|\delta_{3}|\leq|\delta_{T1}|+|\delta_{T2}|$. 2. (b)

Multiplication: $\mathcal{F}_{\langle I,F\rangle}(c_{1})\cdot\mathcal{F}_{\langle I,F\rangle}(c_{2})=c_{3}+\delta_{3}$ with $|\delta_{3}|\leq|\delta_{T1}\cdot\mathcal{F}_{\langle I,F\rangle}(c_{2})|\allowbreak+|\delta_{T2}\cdot\mathcal{F}_{\langle I,F\rangle}(c_{1})|+c_{m}$, where $c_{m}=2^{-F}$ as above. 3. (c)

Division: the operations performed by our controllers in the FWL domain do not include division. However, we do use division in computations at the precision of the plant. Here the error depends on whether the divisor is greater or smaller than the dividend: $\mathcal{F}_{\langle I,F\rangle}(c_{1})/\mathcal{F}_{\langle I,F\rangle}(c_{2})=c_{3}+\delta_{T3}$ where $\delta_{T3}$ is $(\delta_{T2}\cdot c_{1}-\delta_{T1}\cdot c_{2})/(\delta_{T2}^{2}-\delta_{T2}c_{2})$, 3. 3.

Overflow: The maximum size of a real number $x$ that can be represented in a fixed point domain as $\mathcal{F}_{\langle I,F\rangle}(x)$ is $\pm(2^{I-1}+1-2^{-F})$. Numbers outside this range cannot be represented by the domain. We check that overflow does not occur.

0.B.2 Modeling quantization as noise

During any given ADC conversion, the continuous signal will be sampled in the real domain and transformed by $\mathcal{F}\langle I_{c},F_{c}\rangle(x)$ (assuming the ADC discretization is the same as the digital implementation). This sampling uses a threshold which is defined by the less significant bit ($q_{c}=c_{m_{c}}=2^{-F_{c}}$) of the ADC and some non-linearities of the circuitry. The overall conversion is

[TABLE]

If we denote the error in the conversion by $\nu_{k}=y_{k}-y(t)$ where $t=nk$, and $n$ is the sampling time and $k$ the number of steps, then we may define some bounds for it $\nu_{k}\in[-\frac{q_{c}}{2}\ \ \frac{q_{c}}{2}]$.

We will assume, for the purposes of this analysis, that the domain of the ADC is that of the digital controller (i.e, the quantizer includes any digital gain added in the code). The process of quantization in the DAC is similar except that it is calculating $\mathcal{F}\langle I_{dac},F_{dac}\rangle(\mathcal{F}\langle I_{c},F_{c}\rangle(x))$. If these domains are the same ($I_{c}=I_{dac},\allowbreak F_{c}=F_{dac}$), or if the DAC resolution in higher than the ADCs, then the DAC quantization error is equal to zero. From the above equations we can now define the ADC and DAC quantization noises ${\nu_{1}}_{k}\in[-\frac{q_{1}}{2}\ \ \frac{q_{1}}{2}]$ and ${\nu_{2}}_{k}\in[-\frac{q_{2}}{2}\ \ \frac{q_{2}}{2}]$, where $q_{1}=q_{c}$ and $q_{2}=q_{\mathit{dac}}$. This is illustrated in Fig. 3.2 where $Q_{1}$ is the quantizer of the ADC and $Q_{2}$ the quantizer for the DAC. These bounds hold irrespective of whether the noise is correlated, hence we may use them to over-approximate the noise effect on the state space progression over time. The resulting dynamics are

[TABLE]

which result in the following closed-loop dynamics:

[TABLE]