Program synthesis is the problem of computing correct-by-design programs from high-level specifications. Algorithms for this task have made substantial progress over recent years, in particular the architecture known as CEGIS, which is a recent approach to inductive synthesis.

Program synthesizers are an ideal fit for the synthesis of digital controllers, since the semantics of programs precisely capture the effects of finite-precision arithmetic. Surprisingly, the control literature has been oblivious to this relevant connection. A relevant exception is [40], which employs CEGIS on the synthesis of switching controllers for stabilizing continuous-time plants. This work hinges on and is limited by the capacity of the state-of-the-art SMT solvers to reason over linear arithmetic problems. Since this contribution employs finite switching actions for the digital controller, it avoids problems related to finite-precision arithmetic, but potentially suffers from the state-space explosion. Moreover, in [41] 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 both a set of goal states and safe states (this specification is also known as 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 controllers that we seek to synthesize, which are not switching in nature but rather continuous (as further detailed later).

The work in [2] synthesizes stabilizing controllers for continuous plants given as transfer functions, by exploiting bit-accurate verification of software-implemented digital controllers [9]. While this work also uses CEGIS, the approach is restricted to digital controllers for stable closed-loop systems expressed 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 the errors related to finite-precision arithmetic.

Furthermore, unlike the approach in [2], this work reasons over state-space models. A state-space model has well known advantages over the transfer function representation [20], as it allows synthesis of controllers with guarantees on the internal dynamics, e.g., safety. Our work indeed focuses on the safety of internal states (which we assume to fully observe), which is by and large overlooked in the standard (digital) control literature, by default focussed on stability/regulation/tracking properties. Moreover, our work integrates an abstraction-refinement (CEGAR) step inside the main CEGIS loop.

Beyond CEGIS-based architectures, there is an important line of research on provably correct control synthesis for dynamical models, which leverages formal abstractions. The approach underpinning tools such as Pessoa [32], related extensions [4, 30] and applications [51] synthesizes correct-by-design embedded control software by formally abstracting the model to a finite-state machine, and on the formal control synthesis over safety- and reachability-based temporal properties thereon. Whilst in principle this formal abstraction can account for the errors that we can deal with, it is expected to provide a controller of different nature than the one we generate. The obtained finite controller software can then be implemented (refined) over the concrete LTI model. However, relying on state-space discretization, this class of approaches is likely to incur scalability limitations.

Recent results in digital control have focused on separate aspects of discretization, e.g. delayed response [17] and finite-precision arithmetic, with the goal either to verify [16] the correctness of the implementation or to optimize [36] its design.

There are two different problems that arise from finite-precision arithmetic in digital controllers. The first is the error caused by the inability to represent the exact state of the physical system, while the second relates to rounding and saturation errors during mathematical operations. In [19], a stability measure based on the error of the digital dynamics ensures that the deviation introduced by finite-precision arithmetic does not lead to instability. Wu et al. [50] uses μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}-calculus to synthesise directly a digital controller, so that selected parameters result in stable model dynamics. The analyses in [42, 49] rely on an invariant computation on the discrete dynamics using semi-definite programming (SDP): while the former contribution uses a bounded-input and bounded-output (BIBO) notion of stability, the latter employs Lyapunov-based quadratic invariants. In both cases, the SDP solver uses floating-point arithmetic and soundness is checked by bounding the obtained error. An alternative is [37], where the verification of given control code is performed against a known model by extracting an LTI model of the code via symbolic execution: to account for rounding errors, upper bounds of their values are introduced in the verification phase. The work in [38] introduces invariant sets as a mechanism to bound the effect of quantization error on stabilization. Similarly, [29] evaluates the quantization error dynamics and calculates upper and lower bounds for the possible trajectory of the system, up to a finite time. Considering the problem of multi-modal dynamics, [25] uses numerical optimization techniques to learn optimal switching logic for hybrid systems. The last three approaches can be placed within the research area known as “hybrid systems theory.” The present contribution distances itself from all these cognate results.

A large body of work exists on evaluating fixed-point errors and bridging the gap between real value and fixed-point values in synthesis in other application areas and using other algorithmic techniques: Genetic programming has been used to minimise the error in synthesised fixed-point programs [15]; Smoothed proof search reduces the problem of parameter synthesis under boolean and quantitative objectives to a sequence of optimisation problems [12]; Synthesis of optimal fixed-point implementations of floating-point numerical software can be done using testing and induction [24]. All of these works present potential future applications for CEGIS in synthesising implementations in fixed-point arithmetic.

We consider models of physical plants expressed as ordinary differential equations, which we assume to be controllable [5]:1x˙(t)=Ax(t)+Bu(t),\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\dot{x}}(t) = {A}{x}(t)+ {B} {u}(t), \end{aligned}$$\end{document}where x∈Rn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x} \in {\mathbb {R}}^n$$\end{document}, u∈Rp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${u} \in {\mathbb {R}}^p$$\end{document}, A∈Rn×n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A} \in {\mathbb {R}}^{n \times n}$$\end{document}, B∈Rn×p\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${B} \in {\mathbb {R}}^{n \times p}$$\end{document}, and t∈R0+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t \in {\mathbb {R}}_0^+$$\end{document} denotes continuous time. We denote with x(0) the model initial condition, which can be non-deterministic. We assume full observability of the states of the model, namely the output states correspond with the model variables, which we can thus fully access.

Equation (1) is discretized in time [33, 48] with constant sampling intervals, each of duration Ts\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_s$$\end{document} (the sample time), into the difference equation2xk+1=Adxk+Bduk,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {x}_{k+1} =&{A}_d {x}_k+ {B}_d {u}_k, \end{aligned}$$\end{document}where Ad=eATs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A}_d=e^{{A}T_s}$$\end{document} and Bd=∫t=0TseAtdtB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${B}_d = \int _{t = 0}^{T_s} e^{{A} t} dt\ {B}$$\end{document}, and where k∈N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k \in {\mathbb {N}}$$\end{document} is a discrete counter and x0=x(0)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x}_{0}={x}(0)$$\end{document} denotes the initial state. We assume that specifications concern the model in (2), and we plan to devise controllers uk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${u}_k$$\end{document} to meet them. (The more general problem of synthesising controllers u(t) for (1) fall outside the scope of the present work.)

Models (1) and (2) depend on external non-determinism in the form of input signals u(t) and uk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_k$$\end{document}, respectively. Feedback architectures can be employed to manipulate properties and behaviors of the plant: we are interested in the synthesis of digital feedback controllers uk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_k$$\end{document}, as in Fig. 3, as practically implemented on field-programmable gate arrays or on digital signal processors, and as classically studied in [7].

Recall that the states of the model are fully accessible (namely, observable). We consider state feedback control architectures, where uk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_k$$\end{document} (notice we work with the discretized signal) is uk=rk-Kxk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_k = r_{k} - K x_k$$\end{document}. Here K∈Rp×n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K \in {\mathbb {R}}^{p \times n}$$\end{document} is a state-feedback gain matrix, and rk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{k}$$\end{document} is a reference signal (again digital). We will assume rk=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_k=0$$\end{document} (meaning that the reference signal is just a zero signal), thus obtaining the closed-loop model xk+1=(Ad-BdK)xk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x}_{k+1}=({A}_d-{B}_d{K}){x}_k$$\end{document} with the origin as its equilibrium point.

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 dynamical behavior of interest [7]. As argued later in this work, we will target a less standard objective, namely a quantitative safety requirement defined over a convex set around the origin, which opens up to more complex specifications [8, 45]. This is not typical in the digital control literature. We will further precisely account for the digital nature of the controller, which manipulates quantised signals as discrete quantities represented with finite precision. The new approach is fully automated and leverages an approach based on CEGIS, to be introduced shortly.Fig. 3

Closed-loop digital control setup, comprising an analogue model of the underlying real system, alongside a digital controller

In this work we employ the notion of asymptotic stability within the CEGIS loop, as a means for reducing the search space of possible safe controllers, where the notion of a safe controller is defined in the following section. As discussed later, for linear models a safe controller is necessarily asymptotically stable, although the reverse is not true. Qualitatively, (local) asymptotic stability is a property denoting the convergence of the model executions towards an equilibrium point, starting from any states in a neighborhood of the point. In the case of linear systems considered with a zero reference signal (as assumed above), the equilibrium point of interest is the origin (see Fig. 2 for the portrait of an asymptotically stable execution, converging to the origin).

Of interest for this work, it can be shown that a discrete-time LTI model is asymptotically stable if all the roots of its characteristic polynomial (i.e., the eigenvalues of the closed-loop matrix Ad-BdK\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_d - B_d K$$\end{document}) are inside the unit circle on the complex plane, i.e., if their absolute values are strictly less than one [7]. Whilst this simple sufficient condition can be either generalized or strengthen to be necessary, this is not necessary in the context of this work. What is technically key is that in this paper we shall express this asymptotic stability as a specification ϕstability\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _ stability $$\end{document}, and encode it in terms of a check known as Jury’s criterion [18]: this is an easy algebraic formula to check the entries of matrix K, so that the closed-loop dynamics are shaped as desired. We refer the interested reader to [18] for further details about Jury’s criterion, which we omit for the sake of space and since it is a standard result in the control literature that we just leverage for our overall approach.

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 model meets a given safety specification that is characterised by a given set around the origin. A safety specification gives rise to a requirement on the states of the model, namely that they remain within the safe set at all times (that is, over an infinite number of time steps). So the feedback controller (namely the choice of the gain matrix K) must ensure that the state never violates the requirement. Note that an asymptotically 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:3ϕsafety=∀k≥0,⋀i=1nxi̲≤xi,k≤xi¯,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \phi _ safety = \left\{ \forall k\ge 0,\, \bigwedge _{i=1}^{n}{\underline{x_{i}} \le x_{i,k} \le \overline{x_{i}}}\right\} , \end{aligned}$$\end{document}where xi̲\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\underline{x_{i}}$$\end{document} and xi¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{x_{i}}$$\end{document} are lower and upper bounds for the ith coordinate xi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{i}$$\end{document} of state x∈Rn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x\in {\mathbb {R}}^n$$\end{document} at the kth instant, respectively. This requires that the states will always be within an n-dimensional hyper-box.

Beyond the main requirement on safety, it is practically relevant to consider the constraints ϕinput\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _ input $$\end{document} on the input signal uk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_{k}$$\end{document} and ϕinit\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _ init $$\end{document} on the initial states x0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_0$$\end{document}, which we assume have given bounds: ϕinput={∀k≥0,⋀i=1pu̲≤ui≤u¯}\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _ input = \{\forall k \ge 0,\bigwedge _{i=1}^{p} {\underline{u}} \le u_{i} \le {\overline{u}}\} $$\end{document}, and ϕinit={⋀i=1nxi,0̲≤xi,0≤xi,0¯}\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _ init = \{ \bigwedge _{i=1}^{n} \underline{x_{i,0}} \le x_{i,0} \le \overline{x_{i,0}}\}$$\end{document}. The former constraint expresses that the control input, possibly shaped via state-feedback, might saturate in view of physical constraints. Notice that we will assume that the set of initial states lies within the safe set, since the contrary leads to dynamics that are trivially unsafe. Furthermore, whilst the problem can handle general hyper-boxes, we will in practice work with sets that contain the origin, particularly in view of the choice of the reference signal rk=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_k = 0$$\end{document}.

We employ the bounded model checking software tool CBMC [28] to model the controller behavior on digital signals with finite-bit precision. We express the controller semantics as C programs, which is CBMC’s preferred input language. CBMC symbolically executes every instruction in an input program by translating it to a Boolean satisfiability problem, which is satisfiable iff a certain property about the program holds. CBMC’s API allows to specify these properties in the form of assertions in the input program: we set these properties to express that our controllers generate (asymptotically stable and) safe dynamics.

A key contribution of this work is that it precisely replicates the finite-precision arithmetic within the digital controller it synthesizes, thus guaranteeing that controllers implemented with finite-precision arithmetic are safe. The specific components of the model we are concerned with, shown in Fig. 3, are the ADC, the digital controller and the DAC. Details of how we model the behaviour of the finite-precision arithmetic in these components are in Sect. 5. We must model the semantics precisely, in order to capture the full behaviour of the model. More specifically, we encompass the following features:

The ADC converts the analog signal x(t) into a digital signal xk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_k$$\end{document}, which is then fed into the controller, converting a real value to a finite-precision value.

In addition to precisely replicating the finite-precision arithmetic of the digital controller, we must consider that our model itself in (2) employs finite-precision arithmetic to represent the behaviour of the real system. In order to guarantee soundness, we therefore encompass the error that is due to modelling (as opposed to the nature of the digital controller): the representations used in the plant model and its arithmetic operations are carried out at finite precision. More specifically:

We account for the error introduced by finite-precision arithmetic applied over the model variables xk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_k$$\end{document} and uk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_k$$\end{document}, which are actually continuous quantities. We guarantee that the precision we use to represent the model variables is at least as precise as the precision used in the digital controller, and we use interval arithmetic to bound the incurred errors, as further detailed in Sect. 5.

In this paper we will use F⟨I,F⟩(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle I,F \rangle }(x)$$\end{document} to denote a real number x expressed with a fixed-point precision, using I bits to represent the integer part of the number and F bits to represent its decimal part. In particular, F⟨Ic,Fc⟩(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle I_c,F_c \rangle }(x)$$\end{document} denotes a real number x represented at the fixed-point precision of the controller, and F⟨Ip,Fp⟩(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle I_p,F_p \rangle }(x)$$\end{document} denotes a real number x represented at the fixed-point precision of the plant model. Ic\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_c$$\end{document} and Fc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_c$$\end{document} are determined by the controller. We pick Ip\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_p$$\end{document} and Fp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_p$$\end{document} for our synthesis such that Ip≥Ic\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_p \ge I_c$$\end{document} and Fp≥Fc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_p \ge F_c$$\end{document}, so that our controller can be represented accurately in the model domain. Thus any mathematical operations in our modelled digital controller will be in the range of F⟨Ic,Fc⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle I_c,F_c \rangle }$$\end{document}, and all other calculations in our model will be carried out in the range of F⟨Ip,Fp⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle I_p,F_p \rangle }$$\end{document}.

We further employ F⟨E,M⟩(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle E,M \rangle }(x)$$\end{document} to denote a real number x represented in a floating-point domain, with E bits representing the exponent part and M bits representing the mantissa part. In particular, we use F⟨Ec,Mc⟩(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle E_c,M_c \rangle }(x)$$\end{document} to denote a real number represented at the floating-point precision of the controller, whereas F⟨Ep,Mp⟩(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle E_p,M_p \rangle }(x)$$\end{document} denotes a real number represented at the floating-point precision of the plant model.Fig. 4

The general CEGIS framework

In this section, we give a brief description of the CEGIS framework [23, 43], which is illustrated in Fig. 4. CEGIS has been recently developed for the automated synthesis of software programs, and its setup is naturally that of mathematical logic. We consider an input specification of the form∃P.∀a.ϕ(a,P),\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \exists P .\, \forall a.\, \phi (a, P), \end{aligned}$$\end{document}where P ranges over functions (where a function is represented by the program computing it), a ranges over ground terms, and ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} is a quantifier-free logical formula. We interpret the ground terms over some domain D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {D}}$$\end{document}. This is a design problem, where the objective is to synthesise a valid P that satisfies ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} over all the a terms.

CEGIS breaks down this generally hard synthesis problem into two easier parts: an inductive synthesis phase (denoted by Synthesize in Fig. 4) and a verification phase (denoted by Verify in Fig. 4), which interacts via a finite set of tests inputs that is updated incrementally. Given the specification ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document}, the inductive synthesis procedure tries to find an existential witness P satisfying the specification ϕ(a,P)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi (a, P)$$\end{document} for all a in inputs (as opposed to all a∈D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a \in {\mathcal {D}}$$\end{document}). 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 verification phase, which checks whether it is a full solution (i.e., P satisfies the specification ϕ(a,P)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi (a, P)$$\end{document} for all a∈D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a\in {\mathcal {D}}$$\end{document}). 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. If the solution space is finite then the CEGIS loop is guaranteed to terminate by either finding a solution or showing that no solution exists.

In the context of the formal synthesis of safe controllers, of interest for this work, the set of possible inputs corresponds to the set of possible initial states and the candidate program P is a candidate controller K. The synthesis block generates a candidate controller that works for a subset of the possible initial states, and the verifier checks whether the controller works for all possible initial states.

There are a number of well known procedures to perform stability analysis of dynamical models [5]. Here we select the classical Jury’s stability criterion [7], in view of its efficiency and ease of integration within our implementation. This method checks the stability of a model working over the complex domain of its characteristic polynomial S(z), considered in its general form asS(z)=a0zN+a1zN-1+⋯+aN-1z+aN,a0≠0.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} S(z) = a_0z^N+a_1z^{N-1}+\cdots +a_{N-1}z+a_N,{\quad } a_0\ne 0. \end{aligned}$$\end{document}A standard result in Control theory [5] states that this polynomial can be obtained as a function of the state-space matrices Ad,Bd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_d,B_d$$\end{document}, and in particular its order N corresponds to the dimensions of the state variables (above, n). A sufficient condition for asymptotic stability of the closed-loop LTI model [7] is when all the roots of its characteristic polynomial S(z) (which correspond to the eigenvalues of the matrix Ad-BdK\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_d-B_dK$$\end{document}) are inside the unit circle in the complex plane, i.e., when the absolute values of the roots are less than one.

Skipping the full algebraic derivation for the sake of space (this can be found in [7]), the following matrix M with dimension (2N-2)×N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(2N-2)\times N$$\end{document} and elements m(·),(·)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{(\cdot ),(\cdot )}$$\end{document} is built from the coefficients of S(z) as:M=V(0)V(1)⋮V(N-2),\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} M=\left( \begin{array}{c} V^{(0)}\\ V^{(1)}\\ \vdots \\ V^{(N-2)} \end{array} \right) , \end{aligned}$$\end{document}where V(k)=[vij(k)]2×N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V^{(k)} = [v^{(k)}_{ij} ]_{2\times N}$$\end{document} is such that:vij(0)=aj-1,ifi=1v(1)(N-j+1)0,ifi=2vij(k)=0,ifj>n-kv1j(k-1)-v2j(k-1).v11(k-1)v21(k-1),ifj≤n-kandi=1v(1)(N-j+1)k,ifj≤n-kandi=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} v_{ij}^{(0)}= & {} \left\{ \begin{array}{l@{\quad }l} a_{j-1}, &{} \text{ if }~i=1\\ v_{(1)(N-j+1)}^{0},&{}\text{ if }~i=2 \end{array} \right. \\ v_{ij}^{(k)}= & {} \left\{ \begin{array}{l@{\quad }l} 0,&{}\text{ if }~j>n-k\\ v_{1j}^{(k-1)}-v_{2j}^{(k-1)} . \frac{v_{11}^{(k-1)}}{v_{21}^{(k-1)}}, &{} \text{ if }~j\le n-k ~\text{ and }~i=1\\ v_{(1)(N-j+1)}^{k},&{} \text{ if }~j\le n-k ~\text{ and }~i=2\\ \end{array} \right. \end{aligned}$$\end{document}and where k∈Z\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k \in {\mathbb {Z}}$$\end{document} is such that 0<k<N-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0< k < N - 2$$\end{document}.

We have that S(z) is the characteristic polynomial of an asymptotically stable system if and only if the following four conditions Ri,i=1,2,3,4,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_i, i = 1,2,3,4,$$\end{document} hold [7]: R1:S(1)>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_1: S(1) > 0$$\end{document}; R2:(-1)NS(-1)>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_2: (-1)^N S(-1) > 0$$\end{document}; R3:|a0|<aN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_3: |a_0| < a_N$$\end{document}; R4:m11>0∧m31>0∧m51>0∧…∧m(2N-3)(1)>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_4: m_{11}> 0 \wedge m_{31}>0 \wedge m_{51}>0 \wedge \ldots \wedge m_{(2N{-}3)(1)}>0$$\end{document}, where mij\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{ij}$$\end{document} denotes the element in position (i, j) of the matrix M, as defined previously.

Finally, the asymptotic stability property is finally encoded by a constraint expressed as the following formula: ϕstability={R1∧R2∧R3∧R4}.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _ stability = \{R_1 \wedge R_2 \wedge R_3 \wedge R_4\}. $$\end{document}

As described in Sect. 3.5, we use CBMC, a bit-precise bounded model checker, to synthesise and verify candidate controllers. CBMC manipulates precisely the fixed- or floating-point arithmetic used in the controller, as well as the ADC/DAC conversions, according to the IEEE standards.

We use finite-precision arithmetic to model the plant. This is an approximation that speeds up each CEGIS iteration, however it necessitates a further stage where we verify that the errors introduced by the approximation have not resulted in a controller that is unsafe when executed on a model expressed over real numbers. In this stage, we represent the plant model using double-precision floating-point numbers and we use the Boost interval arithmetic library [11] to bound the error in this representation. We employ a compositional numerical library to model the fixed-point arithmetic for the controller1 within double-precision floating-point numbers. We check that the controller is safe starting from each vertex of the set of initial states, and show that this is sufficient to prove safety from any state in this set (see Theorem 1).

We describe here the mathematics behind bounding the errors on the double-precision floating-point numbers. Recall we use F⟨E,M⟩(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle E,M \rangle }(x)$$\end{document} denote a real number x represented in a floating-point domain, with E bits representing the exponent part, and M bits representing the mantissa. In general the representation of a real number using the floating-point domain introduces an error, for which an upper bound can be given [10]. For each number x represented in the floating-point domain as F⟨E,M⟩(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle E,M \rangle }(x)$$\end{document}, we store an interval that encompasses this error. Further mathematical operations performed at the precision F⟨E,M⟩(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle E,M \rangle }(x)$$\end{document} will propagate this error, leading to further errors for which bounds can be derived [10].

The fixed-point arithmetic of the digital controller is performed on the upper and lower bound of the intervals from above independently, and the upper and lower bound of the result is taken as the interval result. For example, consider the conversion from the model precision to controller precision performed by the ADC on a single state value. The state value is represented as an interval {x.high,x.low}\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{x.high, x.low\}$$\end{document}, and the result of the conversion is an interval where the upper bound is the conversion of x.high and the lower bound is the conversion of x.low. Since the precision of the floating-point domain is greater than the precision of the controller, this is guaranteed to bound the real behaviour of the controller.

In this section we will quantify how the finite-precision arithmetic in a digital controller affects the safety and stability properties of an LTI model.

The inductive synthesis phase (synthesize) uses BMC to compute a candidate solution K that satisfies both the stability requirement and the safety specification, within a finite-precision model. In order to synthesize a controller that satisfies the stability requirement, we need the characteristic polynomial of the closed-loop model to satisfy Jury’s criterion [18] (see Sect. 4.1).

We fix an index k, and we synthesize a safe controller by unfolding the transition system (i.e., the closed-loop model) k steps and by selecting a controller K and a single initial state, such that the states at each step do not violate the safety requirement (see Sect. 3.4). That is, we ask the bounded model checker [14] if there exists a K that is safe for at least one x0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_0$$\end{document} in our set of all possible initial states, and a given fixed-point precision for the controller, and where the input signal remains within the specified bounds. The bounded model checker selects the controller K values using a satisfiability solver, i.e., by constraint solving. As such, we do not need to consider traditional control approaches such as pole placement. This approach is sound, i.e., the controller produced is guaranteed to be safe, if the current k is greater than the completeness threshold (see later step). We also assume some finite precision for the model and a given time discretisation, as described in Sect. 3.1. We use a fixed-point precision in this description, given by ⟨Ip,Fp⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle I_p,F_p\rangle $$\end{document}, but if we are considering a floating-point controller we will instead model the plant with a floating-point precision, ⟨Ep,Mp⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle E_p,M_p\rangle $$\end{document}. The checks that these assumptions hold are performed by the subsequent verify stages.

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

The bioreactor benchmark is a linear model of the cell mass concentration controlled through the dilution rate of a bioreactor [39]. The Chen benchmark correspond to a higher-order control system model employed as a case study for model-order reduction techniques [13]. The benchmarks Cruise 1 [20] and Cruise 2 [5] deal with automotive cruise control models, where the control is used to maintain the speed of an automobile at a constant value (after tracking a desired speed reference), and compensating disturbances and uncertainties. The models cstr and cstrtmp describe the pH dynamics of a reaction of an aqueous solution of sodium acetate with hydrochloric acid [46] and the temperature of a reaction [6] in a tank reactor. The DC motor plant describes the velocity dynamics of a direct-current electrical machine. The helicopter benchmark plant describes the transitional and rotational dynamics of a coaxial helicopter. The inverted pendulum benchmark describes a model for the cart position and for the angle of an inverted pendulum placed on the cart, which moves over a track by means of a DC motor. The magnetic pointer benchmark describes a magnetic pointer, whose angular position dynamics is controlled by a magnetic field. The magsuspension describes the dynamics of a magnetic car suspension system. The pendulum plant consists of a swinging point mass suspended from a frictionless pivot by means of a rod with negligible mass. The Regulator consists of a linear model of a synchronous electrical machine [26]. The satellite attitude control system plant describes the dynamics of a satellite attitude, i.e., the orientation angles. An attitude control system must maintain the desired orientation of a satellite with respect to an inertial frame. The springer-mass damper system plant is a standard model for the dynamics of several mechanical systems. The steam drum benchmark describes a linear model for the level dynamics of a boiler steam drum [31]. The suspension models the single-wheel suspension system of a car, namely the relative motion dynamics of a mass-spring-damper model, which connects the car to one of its wheels. The tapedriver benchmark models a computer tape driver, that is a computer storage device used to read and write data on magnetic tapes. The USCG cutter tampa heading angle plant describes the heading angle dynamics of a Coast Guard cutter.

Using the state-space models described in Sect. 8.1, the evaluation study has the following overall experimental goal:

Show that the multi-staged CEGIS approach is able to generate finite-precision digital controllers using fixed-point and floating-point arithmetic in a reasonable amount of time.

We provide the results in Table 1, where: Benchmark is the name of the corresponding benchmark; Order is the number of continuous variables of the model; F⟨Ip,Fp⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle I_p,F_p \rangle }$$\end{document}, where Ip\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_p$$\end{document} and Fp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_p$$\end{document} indicate the integer and fractional parts, and F⟨Ep,Mp⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle E_p,M_p \rangle }$$\end{document}, where Ep\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_p$$\end{document} and Mp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_p$$\end{document} indicate the exponent and mantissa part, are the fixed- and floating-point precisions used to model the plant, respectively; and Time is the total time required to synthesize a controller for the given model. Time-outs are indicated by ✗, where the time-out used is 5 h. The precision for the controller, F⟨Ic,Fc⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle I_c,F_c \rangle }$$\end{document}, is chosen to be Ic=8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_c = 8$$\end{document}, Fc=8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_c = 8$$\end{document} for fixed-point, whereas F⟨Ec,Mc⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}_{\langle E_c, M_c \rangle }$$\end{document} is chosen to be Wc=6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W_c = 6$$\end{document} and Fc=10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_c = 10$$\end{document} for half-precision floating-point format.

We separate our evaluation into two sets of results: fixed- and floating-point. We present a selection of case-studies that we are able to solve either for floating-point or fixed-point. Given a time discretisation, the floating-point controllers in general take longer to synthesise than the fixed-point ones. However, the fixed-point algorithm is often forced to try more time discretisations, because the fixed-point controller lacks the precision to control the system at the coarser time discretisation.Table 1

Synthesis times for fixed- and floating-point controllers

The mean run-time for the successful benchmarks is 225 s for fixed-point and 125 s for floating-point implementations, respectively. The median run-time for the successful benchmarks is 15.7 s for fixed-point and 17.8 s for floating-point implementations, respectively. We consider these times to be short enough to be of practical use to control engineers, and thus assert the success of the overall objective of the study that we have raised above. The completeness threshold depends upon the controller selected, and the SAT/SMT solver picks a controller with a low completeness threshold for all the time discretisations that we solve. The completeness threshold is zero where the controller synthesised results in a system with real valued eigenvalues.

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).2 The provided experimental environment runs multiple discretisations for each benchmark, and lists the fastest as the resulting synthesis time.

In this article we have presented a selection of case studies for which we have been able to automatically synthesise safe controllers. In the full experimental environment we have included further case studies that we have been unable to synthesise controllers for. This is in many cases due to timeouts, especially when the set completeness threshold is too large, or when a controller simply may not exist for a benchmark at a given time discretisation and required controller precision. Yet another source of incompleteness is the inability of the synthesize phase to employ a large-enough precision for the plant model.

We identify the following factors as potential limits to the generality of the obtained results.

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.