Formal synthesis of closed-form sampled-data controllers for nonlinear continuous-time systems under STL specifications

We propose a counterexample-guided inductive synthesis framework for the formal synthesis of closed-form sampled-data controllers for nonlinear systems to meet general STL specifications. Rather than stating the STL specification for a single initial condition, we consider an (infinite) set of initial conditions. Candidate solutions are proposed using genetic programming, which evolves controllers based on a finite number of simulations. Subsequently, the best candidate is verified using reachability analysis; if the candidate solution does not satisfy the specification, an initial condition violating the specification is extracted as a counterexample. Based on this counterexample, candidate solutions are refined until eventually a solution is found. The resulting sampled-data controller is expressed as a closed-form expression, enabling the implementation in embedded hardware with limited memory and computation power. The effectiveness of our approach is demonstrated for multiple systems.


Introduction
Recent years have seen a surge in interest in controller synthesis for temporal logic specifications, realizing complex behavior beyond traditional stability requirements, see, e.g., the recent literature survey in [1]. Originally stemming from the field of computer science, temporal logic has been used to describe the correctness of complex behaviors of computer systems [2]. As it originally dealt with finite systems, (bi-)simulation approaches have been proposed to abstract infinite systems to finite systems [3,4]. However, as a downside, these approaches (e.g., [5][6][7][8]) typically suffer from the curse of dimensionality and return controllers in the form of enormous lookup tables [9].
Where certain temporal logics such as linear temporal logic reason over traces of finite systems, signal temporal logic (STL) reasons over continuous signals [10]. Besides a Boolean answer to whether the formula is satisfied, quantitative semantics of STL has been introduced [11,12], providing a quantititive measure on how robustly a formula is satisfied. These robustness measures enable optimization-based methods for temporal logic, such as model predictive control (MPC) [13][14][15][16][17], optimal trajectory planning [18], reinforcement learning [19], and neural networks [20]. Apart from optimization-based methods, other proposed approaches for STL specifications rely on control barrier functions (CBF) [21,22]. While the work in [21] does not optimize a robustness measure of the STL specification, the computation of the control input for every time step relies on online quadratic optimization. Alternatively, in [23,24] the synthesis for a fragment of STL is reformulated to a prescribed performance control problem, resulting in a continuous state feedback control law.
While (bi-)simulation approaches provide feedback strategies for all (admissible) initial conditions, only a limited number of optimization-based approaches consider a set of initial conditions [1], including [13,15,25]. In [16], tube MPC is used, in which a tube around a nominal initial condition is found for which the robustness measure is guaranteed. Similarly, the control barrier functions in [21] provide a forward invariant set around the initial condition.
In this work, we utilize genetic programming (GP) [26] and reachability analysis [27] to synthesize controllers. The benefit of genetic programming is that it is able to automatically find a structure for the controller, as the right structure is typically unknown beforehand [1]. Genetic programming has been used for formal synthesis for reach-avoid problems in [28,29], in which controllers and Lyapunov-like functions are automatically synthesized for nonlinear and hybrid systems. Also, reachability analysis has been used in formal controller synthesis for reach-avoid problems, e.g., in [25], MPC is combined with reachability analysis, whereas in [30][31][32] synthesizes a sequence of optimal control inputs [32] or linear controllers [30,31] for a sequence of time intervals.
Regardless, to the best of our knowledge, there are no closed-form controller synthesis methods which guarantee general STL specifications for a set of initial conditions. The goal of this work is to synthesize correct-by-construction closed-form controllers for nonlinear continuous-time systems subject to bounded disturbances for STL specifications. Moreover, we consider a sampled-data implementation of the controller, i.e., the controller output is only updated periodically and is held constant between sampling times. To this end, we propose a framework based on counterexampleguided inductive synthesis (CEGIS) (see e.g. [13,[33][34][35]), combining model checking for STL [36], the recent development of counterexample generation using reachability analysis [37], and genetic programming (GP) [26]. This CEGIS approach combines a learning step with a formal verification step, in this case GP and reachability analysis, respectively. Within this framework, violations obtained during verification are used to improve the learning process, until a controller which formally satisfies the desired specification is found. The resulting closed form and the sampled-data nature of the controller both enable its digital implementation.
The main contributions of this work are twofold: first of all, we propose a CEGIS framework combining genetic programming with reachability analysis for the synthesis of closed-form sampleddata controllers for STL specifications. To enable reasoning over reachable sets as opposed to singular trajectories, [36] introduced reachset temporal logic (RTL) and proposed a sound transformation from STL to RTL. Our second contribution is the definition of quantitative semantics for RTL, and proving that the quantitative semantics is sound and complete. Similar to the quantitative semantics of STL, these quantitative semantics provide a measure of how robustly a formula is satisfied.

Preliminaries
The set of real positive numbers is denoted by R ≥0 . The power set of a set S is denoted by 2 S . Finally, an n-dimensional zero vector is denoted by 0 n .

Signal temporal logic
We consider specifications expressed in signal temporal logic (STL) [10], using the following grammar: where ϕ, ϕ 1 , ϕ 2 are STL formula, and h(s) ≥ 0 is a predicate over a signal s : R ≥0 → R n and a function h : R n → R. The Boolean operators ¬ and ∧ denote negation and conjunction, respectively, and U [a,b] denotes the bounded until operator, i.e. until between a and b. We can also define other standard (temporal) operators from (1), such as disjunction ϕ 1 ∨ ϕ 2 := ¬(¬ϕ 1 ∧ ¬ϕ 2 ), next Given a set Y ⊂ R n which can be expressed as we denote the logic function indicating set membership by ϕ Y = i j h ij (x) ∼ 0. The satisfaction relation (s, t) |= ϕ indicates that the signal s starting at t satisfies ϕ. We consider the same definition of the semantics as in [36], which slightly deviates from e.g. [10] w.r.t. the until operator 1 . Since we build upon the results in [36], we have adopted the corresponding definition. STL is equipped with quantitative semantics ρ(s, ϕ, t) that provides a robustness measure of how well a signal s starting at time t satisfies or violates the STL specification [11,12]. If ρ(s, ϕ, t) is negative, lower values imply that ϕ is more strongly violated. Conversely, if ρ(s, ϕ, t) is positive, higher values imply that ϕ is satisfied more robustly.

Reachset temporal logic
Consider a closed-loop system described by the following differential inclusion: where ξ(t) ∈ R n denotes the state, F : R × R n ⇒ R n is a set-valued function, and I ⊂ R n is the set of initial conditions. In this work, we are not only interested in the STL performance of a singular trajectory, but rather of the set of all trajectories satisfying system Σ, defined by Now let us define the reachable set: Definition 2.1 (Reachable set). Given a system Σ, a mapping R e : R ≥0 → 2 R n is an exact reachable set if and only if: That is, a reachable set satisfies that ∀t ∈ R ≥0 , ∀ξ ∈ S(Σ) : ξ(t) ∈ R(t). We use the reachability analysis tool CORA [38], which returns a sequence of sets forming a reachable set given by The STL semantics over singular trajectories does not directly translate to the evaluation over reachable sets. To be able to reason directly over a reachable set, [36] introduced reachset temporal logic (RTL). The RTL fragment relevant for this work is given by: Here, ψ are propositional formulae over states x and φ formulae over a reachable set R : The semantics is defined as follows: Consider the following assumption on an STL formula ϕ, used throughout this work:  (2), let ϕ be an STL formula satisfying Assumption 2.1 for some value c, and R(t) be the reachable set of Σ in the form of (5) with t i+1 − t i = c. The transformation Υ from [36], bringing the STL formula ϕ into an RTL formula φ = Υ(ϕ), is sound, i.e.: This transformation Υ yields RTL formulae of the form where I, J i , K ij are finite index sets and ψ ijk are non-temporal subformulae. As can be seen, j closely relates to a time step c/2, whereas i and k relate to the number of conjunctions and disjunctions. The reachable set in (5) is formed by the reachable sequence R, which partitions time into an alternating sequence of points and open intervals. Similarly, the transformation from STL to RTL transforms the reasoning over an infinite set of time instances to reasoning over an alternating sequence of points and intervals. In this partition, the value c/2 can be seen as the time step between the points and a time interval. Due to this time partition, the transformation Υ is a sound transformation, but in general not complete, i.e., the converse of (6) does generally not hold. Therefore, the transformation Υ is subjected to some conservatism 2 , which can be reduced by taking smaller values of c.

Genetic programming
The controllers in this work are synthesized using genetic programming (GP) [26], a variant of genetic algorithms (GA) [39], which evolves entire programs rather than optimizing parameters. In our case, the evolved program is a controller based on elementary building blocks consisting of state variables and basic functions such as addition and multiplication. Within genetic programming, a candidate solution, called an individual, is represented by a data structure enabling easy manipulation, such as an expression tree. This data structure is called the genotype, whereas the individual itself, e.g., an analytic function, is referred to as the phenotype. A pool of individuals, called the population, is evolved based on a cost function, called the fitness function, which assigns a fitness score to all individuals. Depending on the fitness score, individuals can be selected to be recombined or modified using genetic operators, such as crossover and mutation. In the former, two subtrees of individuals are interchanged, whereas in the latter, a random subtree is replaced by a new random subtree. Each genetic operator has a user-defined rate, which determines the probability of the operator being applied to the selected individuals. A number of individuals are selected based on tournament selection: a fixed number of individuals are randomly selected from the population, and the individual with the highest fitness is returned. The process of selection and modification through genetic operators is repeated until a new population is created. The underlying hypothesis is that the average fitness of the population increases over many of these cycles, which are referred to as generations. The algorithm is terminated after a satisfying solution is found or a maximum number of generations is met. We use the variant grammar-guided genetic programming (GGGP) [28,40], which utilizes a grammar to which all individuals adhere: the population is initialized by creating random individuals adhering to the grammar and the used genetic operators are defined such that the resulting individuals also adhere to the grammar. The grammar is defined by the tuple (N , S, P), where N is a set of nonterminals, S a starting tree, and P a set of production rules, which relate nonterminals to possible expressions. An example of a grammar is shown in Figure 1a. In this grammar, the nonterminals correspond to polynomials pol , monomials mon over time t, and constants const . The starting tree S restricts the class of controllers to time-varying state feedback laws, linear in the state x ∈ R 2 . The genotype of an individual is constructed using the grammar as follows: starting with the starting tree, for all leaf nodes containing a nonterminal, a subtree is randomly selected from the corresponding production rules and put under the leaf node. This procedure is repeated until all leaf nodes are free of nonterminals. To prevent infinite depth trees due to recursive production rules, all recursive rules are omitted from the production rules after a fixed number of expansions of the nonterminals. To obtain the phenotype, all nonterminal nodes are replaced with their underlying subtrees, and the resulting expression tree is expressed as an analytic expression. Given the grammar in Figure 1a, an example of a fully expanded genotype is shown in Figure 1b, which has the corresponding phenotype of 9.5x 2 + 4.2tx 2 . The grammar-aware genetic operators are defined as follows: In the crossover operator, given two individuals, for each individual a random subtree with the same nonterminal root is selected and these are interchanged. In the mutation operator, a random subtree is selected and replaced with a newly grown subtree with the same nonterminal root.

Problem definition and solution approach
We consider nonlinear systems subject to disturbances of the form: Note that under this model, the disturbance ω ∈ Ω can change at every time instant. In this work we consider disturbance realizations w : R ≥0 → Ω, which are timedependent realizations of this uncertainty parameter ω. We consider sampled-data time-varying state-feedback controllers κ : where t k denotes the k-th sampling instant, t 0 = 0, and η is the sampling time. This results in a closed-loop system of the form (2) with, ∀t ∈ [t k , t k + η), The goal of this paper is formalized in the following: In Theorem 1, we have seen that (10) can be proven by translating the STL formula ϕ to the RTL formula φ using the transformation Υ, and subsequently proving (R, 0) |= φ. In this work, we propose a counterexample-guided inductive synthesis (CEGIS) framework to synthesize a controller such that (R, 0) |= φ, thereby solving Problem 3.1. The proposed framework consists of iteratively proposing a controller obtained through GGGP 3 and then formally verifying the RTL formula using reachability analysis. The proposed controller is designed based on a set of simulated trajectories, which correspond to pairs of initial conditions and disturbance realizations. The underlying idea is that these simulations are relatively fast to compute and provide a sensible search direction for the synthesis, whereas the reachability analysis verifies the resulting controller. For a given open-loop system Σ ol , STL formula ϕ, and grammar (N , S, P), the algorithm is initialized as follows: I2) The set of pairs of initial conditions and disturbance realizations I is initialized by randomly choosing n s initial conditions {x 1 , . . . , x ns } ⊂ I, with random disturbance realizations w i : Given the initialized data, the algorithm goes through the following cycle, illustrated in Figure 2, where each cycle is referred to as a refinement: A1) A candidate solution is proposed using GGGP, based on simulation trajectories corresponding to the set I.
A2) For the given candidate controller, the reachable set is computed.
(b) (R, t) |= φ, and a counterexample is extracted in the form of an initial condition for which there exists a disturbance realization s.t. the RTL specification is violated. For this initial condition, a disturbance realization is optimized. This pair of initial condition and disturbance realization is added to I and the algorithm returns to step A1).
(c) (R, t) |= φ and a maximum of refinements is reached, therefore the algorithm is terminated.
To quantify the violation or satisfaction of an RTL formula, we introduce quantitative semantics for RTL in the next section. The proposal of a candidate controller in step A1) is discussed in Section 5. The verification and counterexample generation in step A3) is discussed in Section 6.

Quantitative semantics
Inspired by the quantitative semantics of STL [11,12], we define quantitative semantics for RTL in this section. These quantitative semantics provide a robustness measure on how well the formula is satisfied. For an RTL formula φ with propositional subformulae ψ, the quantitative semantics is given by functions P (R, φ, t) and (x, ψ), respectively, recursively defined as: The quantitative semantics of STL are sound and complete [12,41]. The quantitative semantics of RTL also have these properties: Theorem 2 (Soundness and completeness). Let φ be an RTL formula, R a reachable set, and t a time instance, then: This is because on the boundary of an inequality, the distinction between inclusion or exclusion is lost within the quantitative semantics. That is, if (x, ψ) = 0, we also have (x, ¬ψ) = 0, hence the quantitative semantics of two mutually exclusive logic formulae evaluate to the same value.
The proof of Theorem 2 can be found in Appendix A. Consider an STL formula ϕ satisfying Assumption 2.1 and the corresponding RTL formula φ = Υ(ϕ) in the form of (7). Using the equivalences a (φ 1 ∧ φ 2 ) = a φ 1 ∧ a φ 2 and rewriting ψ ijk in disjunctive normal form, we can express the RTL formula as: where A ijk and B ijk a denote finite index sets, ∼∈ {≥, >}, and h ijk ab (x) ∼ 0 is a predicate over x. Using the quantitative semantics defined in Section 4, the robustness measure of this RTL formula is given by

Candidate controller synthesis
In this section, we detail step A1) of the proposed algorithm in Section 3, i.e., the proposal of a candidate controller. The candidate controller is synthesized using GGGP, by maximizing an approximation of the robustness measure, which is based on a finite number of simulated trajectories. The sampling time is equal to c/2 to coincide with the time instances at which the robustness measure P (R, φ, 0) is evaluated. For an RTL formula of the form (7), the first and the final time instances of relevance τ 0 and τ f , are given by τ 0 = 0 and τ f = c 2 max i∈I |J i |, respectively. Given a candidate controller κ : R ≥0 × R n → R m , a set of pairs of initial conditions and disturbance realizations I, and a time instance τ q , we consider an approximated reachable setR κ I (τ q ) formed by all corresponding simulated trajectories x : R ≥0 → R n : Provided this setR κ I (τ q ), we approximate the robustness measure by P (R κ I , φ, 0).

Outline of the candidate controller synthesis
The proposal of a candidate controller in step A1) undergoes the following steps, which are also illustrated in Figure 3: A1.a) We synthesize an analytic expression κ : R × R n → R m by using GGGP to solve: If for the resulting controller κ the robustness measure approximation P (R κ I , φ, 0) is negative, this optimization step in (13) is repeated. Otherwise, the algorithm continues to the next step.
A1.b) For each initial condition x i in I, an analytic expression for a disturbance realization w i : R → Ω is synthesized using GGGP, in which the robustness measure approximation is minimized, i.e.: arg max If the corresponding robustness degree approximation P (R κ I , φ, 0) is negative, the algorithm returns to step A1.a). Otherwise, if for all updated disturbance realizations the robustness measure approximation is positive, i.e., ∀i P (R κ {(x i ,ω i )} , φ, 0) > 0, the algorithm returns a candidate controller.

Reference-tracking controllers
To speed up the synthesis, we impose a structure to the solution, based on a nominal reference trajectory x ref (t) and a corresponding feedforward input u ff (t). That is, we consider a time-varying reference-tracking controller of the form: where κ fb : R ≥0 × R n → R m is a time-varying feedback controller. The feedforward input and reference trajectory can be computed beforehand as follows: i.e., maximizing: arg max Using the synthesized pair (u ff (t), x ref (t)), the user-defined grammar used within GGGP can be used to enforce the structure of a time-varying reference controller in (15)

Reachability analysis and verification
In this section we detail step A3) of the algorithm. In this work we consider polynomial zonotopes [42] as the set representation of the reachable set, motivated by its useful properties for counterexample generation, discussed in [37]. Intuitively, polynomial zonotopes maintain the dependencies between points in subsequent reachable sets under the reachability analysis operations. This enables the extraction of an initial condition corresponding to a point for which the specification is violated.
Using this method, we construct a counterexample in the form of a pair of initial condition and disturbance realization (x, w), such that the corresponding trajectory results in a violation of the RTL formula. After the reachability analysis, the algorithm undergoes the following steps: B1) For all subformulae φ ijk in (11b), the corresponding robustness sub-score (12b) is computed by solving the following nonlinear optimization problem 4 over the corresponding set R(jc/2): p * ijk = min B2) Given the robustness sub-scores p * ijk , compute the full robustness measure (12a): B3) As we rely on nonlinear optimization, we cannot guarantee to find the global optimum p * , but rather an upperboundp, such that P (R, φ, 0) = p * ≤p. Givenp, either: (a)p < 0, hence the RTL specification is violated. In this case, given the argument (ijk) * solving (17), we extract an initial condition x 0 corresponding to x (ijk) * , as described in [37]. For this initial condition x 0 , a disturbance realization w is synthesized similarly to step A1.b), i.e., GGGP is used to solve: The pair (x 0 , w) is subsequently added to I. This new set I is then used to improve upon the synthesized controller in step A1).
(b)p ≥ 0, hence the RTL specification is potentially satisfied. However, to guarantee this, we will perform an additional verification step, based on Satisfiability Modulo Theories (SMT) solvers [43], which are capable of verifying first-order logic formulae. The subformula (11b) holds if the following first-order logic formula holds: ∀x ∈ R j c 2 : where again ∼∈ {≥, >}. Suitable SMT solvers to verify (18) include Z3 [44] when R(jc/2) and h ijk ab are expressed as polynomials, and dReal [45] when these are expressed as general nonlinear expressions 5 . Given the Boolean answers to the subformulae in (11b) for all ijk, it is trivial to compute the Boolean answer to (11a).

Dealing with conservatism
Due to both the conservatism in the reachability analysis and the transformation from STL to RTL, it is possible that (R, 0) |= φ, whereas ∀ξ(0) ∈ I, (ξ, 0) |= ϕ, i.e., the desired STL specification holds for all initial conditions, whereas based on the reachability set, the RTL specification is not met. To counter this, the reachability analysis can be made less conservative by refining settings such as the time steps or Taylor order (see [47]). Secondly, the transformation Υ could be performed for a smaller time-discretization parameter c to obtain a less conservative RTL formula φ.
Issues due to conservatism can also be dealt with within the synthesis of a candidate controller in step A1). For example, the population of controllers within GGGP could be further optimized w.r.t. the robustness measure approximation, such that the added robustness could potentially compensate for the conservatism within the reachability analysis. Additionally, controllers can be optimized with respect to both robustness measure and complexity, as less complex controllers  Car 7  14  14 14 14 30  10  3  3  20  10 10  3  Path-planning 10 28  28 14 14 30  50  3  3  40  40 10  3  Aircraft 5  28  42 14 14 50  50  5  5  40  60 10  3   Table 2: Production rules P. might result in less conservatism within the reachability analysis. This results in a multi-objective optimization problem. In this work we consider the fitness criteria for complexity to be defined as the the number of nonterminals of an individual. We use the non-dominated sorting algorithm NSGA-II [48], a Pareto-optimal-aware sorting algorithm, which ranks candidate controllers based on the Pareto-optimality of both fitness criteria. This rank is then used as fitness value within the selection. To make sure the controller with the best robustness measure is always maintained within the population, this controller is always directly copied into the new generation. Finally, there is a gap between the approximated reachable setR κ I and the reachability analysis. There are two sources that can cause significant mismatches between this approximation and actual robustness measure. The first source is truncation errors of the integration scheme. Secondly, due to the added conservatism within reachability analysis, the reachable set can contain additional trajectories which are not admitted by the original system. To bridge this mismatch, we can consider an optional error signal ε added to the simulated trajectory x(τ q ), which is co-synthesized with the disturbance realizations, as will be shown in the case studies in the next section.

Case studies
In this section we demonstrate the effectiveness of the proposed framework on a car benchmark, path-planning problem and aircraft landing manoeuvre. The case studies are performed using an Intel Xeon CPU E5-1660 v3 3.00GHz using 14 parallel CPU cores. The GGGP algorithm is implemented in Mathematica 12 and the reachability is performed using CORA in MATLAB. For the nonlinear optimization and verification in Section 6, we use particle swarm optimization of the global optimization toolbox in MATLAB, and the SMT solver dReal with δ = 0.001, respectively. Across all benchmarks, the probability rate of the crossover and mutation operators being applied on a selected individual are 0.2 and 0.8, respectively. Each generation, parameters within an individual are optimized using Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [49], based on the same fitness function as used for GGGP. More specifically, we use the variant sep-CMA-ES [50], due to its linear space and time complexity. Benchmark-specific settings are shown in Table  1, which include the number of simulations n s , number of individuals, and the number of GGGP and CMA-ES generations. Note that the number of GGGP generations for κ and w i is the number of generations per step A1.a) and A1.b), and not the total of GGGP generations per proposal of a controller in step A1), which depends on the number of times step A1.a) and A1.b) are repeated. For each case-study, we use a grammar with nonterminals and production rules as shown in Table 2. These nonterminals correspond to general expressions expr , trigonometric functions trig , polynomial expression pol , monomials t and constants const . The expressions are formed by polynomials, a product of a polynomial and trigonometric functions, and a sum of two expressions. The trigonometric functions are restricted to hyperbolic tangents, sines and cosines with polynomial arguments. The polynomials are restricted to polynomials over time t. Note that per case study, different starting trees are used, such that potentially only a subset of the grammar is available. E.g., if the starting tree is pol , candidate solutions are restricted to polynomial solutions.
We use Runge-Kutta as numerical integration scheme. To keep a constant number of initial conditions in I, counterexamples are added using a first-in, first-out principle. To compensate for the gap between the simulation and the reachability analysis (as discussed in Section 7), we consider an added error signal bounded by the scaled vector field of the dynamics f , parameterized by where δ is a constant and σ : R ≥0 → [−1, 1] n×n a time-varying diagonal matrix which determines the sign and magnitude of the error signal. The constant δ is optimized after each reachability analysis such that the mismatch between the robustness measure and the approximated robustness measure is minimized, i.e.: arg min where {(x, w)} is the counterexample pair computed in Section 6. Finally, in reporting the synthesized controllers, its parameters are rounded from six to three significant numbers for space considerations.

Car benchmark
Let us consider a kinematic model of a car from [30]: where the states x 1 , x 2 , x 3 , x 4 denote the velocity, orientation, and x and y position of the car, respectively. Furthermore, u 1 and u 2 denote the inputs and w 1 and w 2 disturbances. The sampling time of the sampled-data controller is set to be 0.025 seconds. Similarly to [30], we consider a "turn left" maneuver over a time interval T = [0, 1], where within T , the trajectories stay within the safe set S and at the final time instant, the system is in the goal set, captured by the STL specification: We consider the following safe set S and goal set O: To guide the synthesis, we impose the reference-tracking controller structure from Section 5.2 and therefore we first design a feedforward signal and reference trajectory using GGGP. For u ff , x ref , we use polynomial expressions as a function of time t, for the feedback law κ we restrict the search space to reference-tracking controllers which are linear in the tracking error and polynomial in time: and for w i we consider saturated polynomials in time. This is done using the grammar with starting trees: Here, sat (ω i ,ωi) denotes a saturation function such that w i (t) ∈ Ω, where sat (ω i ,ωi) (x) = max(ω i , min(x, ω i )). Finally, for each disturbance realization, we co-evolve the error signal ε i in (19), which is dependent on the candidate controller κ and disturbance realization w i : where diag denotes a diagonal matrix. For the simulations and reachability analysis, we use a sampling time of 0.025 seconds and 0.0125 seconds, respectively. First, a feedforward control input and reference trajectory for a nominal initial condition are synthesized as described in Section 5.2. An example of a found feedforward controller and corresponding reference trajectory are shown in Table 5. For 10 independent runs, the average synthesis time of u ff and the reference trajectory per dimension x ref,i is shown in Table 4. Using these u ff and x ref as building blocks for the controller, κ is synthesized as described in step A1). An example of a synthesized K(t) in (22)  The corresponding reachable set is shown in Figure 4. We observe that the final reachable set is not within the goal set. The red dots represent the violation and the corresponding initial condition. After refining the controller iteratively, an example of a controller satisfying ϕ 1 after 3 refinements is shown in Table 5. The corresponding reachability analysis is shown in Figure 5 and it shows that for this controller the controller specification is formally met. For 10 independent synthesis runs of κ, statistics on the number of generations, number of refinements, complexity in terms of number of non-terminals, and computation time is shown in Tables 3 and 4 and Figure 9. In most cases, a solution was obtained around 3 refinements. However, due to the stochastic nature of the approach, in one case it took 20 refinements before a solution was found.

Input saturation
In our general framework, we do not canonically consider input saturation. Input saturation can be considered in multiple ways, such as restricting the grammar of the controller to include a saturation    function, or even a continuous approximation using e.g. a sigmoid function. However, the downside of such an approach is that the reachability analysis under these functions is typically challenging for state-of-the-art reachability tools, due to the strong nonlinearity or hybrid nature. Instead, for illustrative purposes, we incorporate the constraint within the STL specification, such that for all states in the reachable set the saturation bounds are not exceeded. Let us revisit the car benchmark, where we consider the same input constraints as in [30], namely u ∈ U = [−9.81, 9.81 The STL specification is extended to: with The synthesis statistics are shown in Tables 3 and 4 and Figure 9. An example of a synthesized  Figure 4. (22) is given by The corresponding reachability set is shown in Figure 6. In most cases, a solution was found in around 4 to 5 refinements, with the exceptions of two runs with 20 and 40 refinements, respectively.

Path planning for simple robot
Let us consider the path planning problem for a simple robot adopted from [17]. We deviate from [17] in considering the system in continuous time and consider bounded disturbances. The system is described by: where the state vector represents the x-velocity, y-velocity, x-position and y-position, respectively. The sampling time of the sampled-data controller is set to be 0.5 seconds. Similar to [17], we consider the specification in which the system needs to remain in a safe set S and eventually visit regions P 1 , P 2 and P 3 : [5,25] φ P1 ∧ ♦ [5,25] φ P2 ∧ ♦ [5,25] In [17], the input is constrained s.t. u ∈ U = [−1, 1] 2 . Similar to Section 8.2, we impose this constraint through the STL specification, yielding the following STL specification: where U is given by (24). We consider the same controller structure and grammar as the previous benchmark, with the exception of the grammar of the feedfoward input and reference trajectory. For Reachable set of the input over time. Yellow: initial set, gray: reachable set, red: safe set S and input constraints, green: target sets P 1 , P 2 , and P 3 , black: selection of simulated trajectories, blue: reachable sets at certain time instances within one of the target sets.
these elements, we extend the grammar to expressions which can include trigonometric functions, by using the grammar in Table 2 and the starting trees S u ff = ( expr , expr ) and S x ref,i = expr . For the simulations and reachability analysis, we use a sampling time of 0.5 seconds. The statistics on the synthesis is again shown in Tables 3 and 4 and Figure 9. An example of the controller elements u ff , x ref and K(t) of a synthesized controller are shown in Table 5. The corresponding reachable set of the state and input is shown in Figure 7. Across 10 independent runs, commonly in 1 to 2 refinements a solution was found, with one run requiring 8 refinement.

Landing maneuver
Let us consider the landing aircraft maneuver, adopted from [5]. The system model is given by D(ν 2 , x 1 ) = (2.7 + 3.08(1.15 + 4.2ν 2 ) 2 )x 2 1 , L(ν 2 , x 1 ) = (68.6(1.25 + 4.2ν 2 ))x 2 1 , where the states x 1 , x 2 , x 3 denote the velocity, flight path angle and the altitude of the aircraft, ν i denotes a disturbed input, where u 1 denotes the thrust of the engines and u 2 the angle of attack. Finally, D(ν, x 1 ) and L(ν, x 1 ) denote the lift and drag, respectively, and m = 60 · 10 3 kg, g = 9.81m/s 2 . The sampling time of the sampled-data controller is set at η = 0.25 seconds. Compared to [5], we do not consider measurement errors, but the proposed framework can be adapted arbitrarily to accommodate this type of disturbance. We define the following safe set, goal set and input bounds: and consider the following specification: [18,20] where the set U is given by (24). That is, trajectories are always within the safe set and satisfy the input constraints, until between 18 and 20 seconds the goal set is reached. We use the same controller structure and grammar as the path-planning problem. For the simulations and reachability analysis, we use a sampling time of 0.25 seconds. The algorithm settings are shown in Table 1. The statistics of 10 independent synthesis runs are again shown in Tables 3 and 4 and Figure 9. An example of the controller elements u ff , x ref and K(t) of a synthesized controller are shown in Table 5. The corresponding reachable set of the altitude over time, as well as the reachable sets of the pitch angles at multiple time instances are shown in Figure 8.

Discussion
In this section we discuss the main results from Section 8 and compare them to the results in the literature. Recall that a GGGP generation is the cycle of creating a new population through fitness evaluation, selection and applying genetic operators. A refinement is defined as the cycle of proposing a candidate solution based on GGGP, validation using reachability analysis, and extracting counterexamples. Therefore, in each refinement, there are one or multiple GGGP generations. First of all, Figure 9a shows a polynomial relation between the number of refinements and the total number  of GGGP generations. Secondly, Figure 9b shows a polynomial relation between the number of refinements versus the total computation time. Finally, Figure 9c illustrates that more refinements does not imply that complexity of the controller increases. However, the complexity of the found controller does seem to be dependent on the system and STL specification. While the computation time is related to the number of refinements, this relationship depends on the STL specification and the dynamics. For the car benchmark without and with input constraints, we observe that the added constraints within the STL specification increased the required number of generations, and typically required more time per refinement. Hence, the total computation time heavily depends on the STL specification, as expected. Additionally, we observe an increase in the median of the complexity of the resulting controllers. The input-constraint car and path-planning benchmarks are both four-dimensional systems, where the STL specification of the latter is more involved. Regardless, the path-planning problem has a lower computation time and requires less generations and number of refinements, indicating a dependency between the computation time and the dynamics of the system, which is also as expected.
In [30], the synthesis time for the car benchmark is around 10 seconds, which is significantly shorter than the synthesis time of the proposed framework. The resulting controller consists of a linear controller for each sampling time, resulting in 10 controllers in total. For longer time horizons and/or finer time discretization, this number of controllers increases. On the other hand, our method is able to find a single controller, independent of the sampling time.
Comparing the proposed framework to MPC approaches, such as the approach used for the path-planning problem in [17], we obtain a closed-form controller for which the STL specification is guaranteed, whereas MPC-based approaches require online optimization to compute the controller input.
For the aircraft benchmark, the abstraction-based method in [5] yields a controller which can be seen as a look-up table, in which the state-space is partitioned into over 2.26 million states. For each region within this partition, a finite set of admissible inputs are stored, resulting in a nondeterministic controller. This controller is stored as a binary decision diagram (BDD) with the size of 2.87 MB and can be further reduced to 0,15 MB by removing nondeterminism (see [9]). Additionally, to parse the BDD format, special libraries are required. On the other hand, our controller can be stored in a text file within 916 bytes and is simply expressed as single analytic function. The synthesis time in [5] is 674 seconds for the abstraction and 26 seconds for the controller synthesis, which is again significantly shorter than the presented framework. However, one could leverage that storage space in embedded hardware is finite, whereas bounds on offline computation time are typically less restrictive. Moreover, the abstraction yields a finite transition system with 9.38 · 10 9 transitions. For higher dimensional systems, due to the curse of dimensionality, the platform used to synthesize the controller can run into memory constraints. While the proposed framework in this work is computationally intensive, it does not suffer from the same curse of dimensionality w.r.t. memory constraints.
Relying on GP, the proposed framework is not a complete method. That is, the method is not guaranteed to find a solution in a finite number of iterations, regardless of its existence. Nevertheless, for the presented case-studies, in 10 independent runs a solution was always found. Since the search space is navigated nondeterministically, we observed that the number of GGGP generations, number of refinements and computation time can vary significantly for each run.
Across all benchmarks, the offline computation time for the proposed method is significantly larger than corresponding references. However, there are several elements in which the computation time can be improved. First of all, GGGP and the reachability tool are implemented in Mathematica and Matlab. By implementing these elements in lower-level programming languages, a speed-up is expected. Analyzing the contribution to the total computation time in Table 4, we observe that, in general, GGGP takes up the majority of the computation time. GGGP is highly parallelizable and is in this work not fully exploited, as we only consider 14 individuals, matching the number of used processor cores. By further exploiting the parallelazible nature of GGGP and therefore exploring a larger part of the search space each generation, a significant speed-up is expected. Thirdly, by limiting to a fragment of STL, e.g., by restricting h(s) to be linear, computing the robustness degree can be simplified and therefore improve the computation time of counterexamples. If additionally the robustness measure is upper bounded in a non-conservative manner, the usage of SMT solvers becomes redundant. This would significantly reduce the computation time for benchmarks such as the path-planning problem. Finally, we imposed input constraints through the STL specification. By using saturation functions in our grammar, the input constraints are satisfied by definition, simplifying the synthesis. However, as caveat, discontinuous functions such as saturation functions significantly complicate the reachability analysis.

Conclusion
We have proposed a framework for CEGIS-based correct-by-design controller synthesis for STL specifications based on reachability analysis and GGGP. The effectiveness has been demonstrated based on a selection of case-studies. While the synthesis time is outmatched by methods solving similar problems, the proposed method results in a compact closed-form analytic controller which is provably correct when implemented in a sampled-data fashion. This enables the implementation in embedded hardware with limited memory and computation resources.