Luca Greco

Hugues Mounier

Maria Bekcheva

An approximate characterisation of the set of feasible trajectories for constrained flat systems 1

Keywords: Constraints, Differentially flat systems, Bézier curves, Trajectory generation, Semialgebraic sets

We study differentially flat systems, subject to constraints on the input, the state and their derivatives. We are interested in providing an approximate characterisation of the set of trajectories satisfying both the equations of the system and the constraints. The flatness property and the choice of a polynomial parametrisation in terms of Bézier functions enable a finite dimension formulation of the problem in the parameter space of the Bézier control points. Such formulation provides a unified method to express a large variety of relevant constraint types such as, for instance, physical obstacles, actuator limitations, tubes around nominal trajectories, performance requirements and energy restrictions. We present a simple and effective inner approximation of the aforementioned set of trajectories obtained by simple algebraic manipulations of the original constraints formulation. The symbolic nature of the proposed characterisation allows a certain degree of adaptation to changes in the constraints without the need for any re-computation. The entire method (formulation an characterisation) is illustrated via a classical benchmark mechanical system.

Introduction

Steering the state of a system to a desired position or through a defined path is a fundamental control problem. It becomes notoriously difficult if the system is nonlinear and the tracking has to be performed while fulfilling physical constraints on the actuators and/or the state. To alleviate such complexity, a 2-dof scheme, separating the feedforward and feedback parts (see, e.g. [START_REF] Hagenmeyer | Robustness analysis with respect to exogenous perturbations for flatnessbased exact feedforward linearization[END_REF]), is often adopted. The feedback stage has to take place entirely online, while the feedforward part can be pre-computed offline, thus enabling the generation of trajectories satisfying the constraints.

A particularly skilful means to simplify the trajectory generation problem is to exploit some kind of system inversion or parametrisation property such as the differential flatness [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and examples[END_REF]), the orbital flatness Email addresses: luca.greco@centralesupelec.fr (Luca Greco), hugues.mounier@universite-paris-saclay.fr (Hugues Mounier), maria@inagral.com (Maria Bekcheva). 1 This paper was not presented at any IFAC meeting. 2 Corresponding author, Tel. +33 1 69 85 17 33. [START_REF] Fliess | A Lie-Backlund approach to equivalence and flatness of nonlinear systems[END_REF]) or the Liouvillian character [START_REF] Chelouah | Extensions of differential flat fields and Liouvillian systems[END_REF]). Under these structural properties, all system variables come as a direct (i.e. functional) expression of a particular (possibly fictitious) output, named a flat (or Liouvillian) output, and of a finite number of its derivatives. This parametrisation has a number of advantages, among which we mention: 1) An automatic satisfaction of the system's equations by the state and input induced by the parametrisation. Their trajectories will be named feasible in the sequel. 2) A direct mapping of the constraints in the flat output space. 3) A natural fulfilment of such constraints that is built in the flat output trajectory design. 4) An intrinsically unified formulation of any kind of constraints (input, output and state).

Most of the literature on steering for both differentially flat and general systems with constraints makes use of optimisation procedures, sometimes with collocation and/or discretisation. Moreover, the envisioned constraints are quite often in the state and input, but not in their successive derivatives. With the exception of [START_REF] Suryawan | Splines and polynomial tools for flatness-based constrained motion planning[END_REF] in the case of polynomial systems, when such collocation or discretisation procedures are involved, usually no guarantee of constraint fulfilment between the collocation or discretisation points can be given (see, e.g. [START_REF] Faiz | Trajectory planning of differentially flat systems with dynamics and inequalities[END_REF], [START_REF] Petit | Inversion based constrained trajectory optimization[END_REF], [START_REF] Oldenburg | Flatness and higher order differential model representations in dynamic optimization[END_REF], [START_REF] Joubert | An interactive tool for designing quadrotor camera shots[END_REF], [START_REF] Roberts | Generating dynamically feasible trajectories for quadrotor cameras[END_REF]).

Several authors sought to ensure the satisfaction of constraints over the entire time horizon without any post elaboration. A simple way to avoid any recourse to collocation or discretisation is to specialise the flat outputs to a linear combination of known functions. While the initial problem is infinite dimensional by nature (search in a general functional space), the one after specialisation is a finite dimensional one (e.g. in the weights of the linear combination).

For linear systems, [START_REF] Henrion | LMIs for constrained polynomial interpolation with application in trajectory planning[END_REF] parametrise the flat output with polynomial functions, which are naturally closed under differentiation; the corresponding constraint problem then amounts to a polynomial nonnegativity one, which is coped with a sum-of-squares decomposition and solved with semidefinite programming. [START_REF] Louembet | Motion planning for flat systems using positive B-splines: An LMI approach[END_REF] also make use of piecewise polynomial representations (B-splines) and a sum-of-squares formulation, while [START_REF] Suryawan | Splines and polynomial tools for flatness-based constrained motion planning[END_REF] resort to the convex hull property of B-splines, which introduces some conservatism. This convex hull property is also used by [START_REF] Flores | Trajectory generation for differentially flat systems via NURBS basis functions with obstacle avoidance[END_REF] together with NURBS, the optimal problem being solved through sequential quadratic programming. For nonlinear systems, some restrictions or simplifications are usually adopted. [START_REF] Louembet | Motion planning for flat systems using positive B-splines: An LMI approach[END_REF], for instance, consider an inner polytopic approximation of the constraints (actually, only linear constraints can be solved); [START_REF] Suryawan | Splines and polynomial tools for flatness-based constrained motion planning[END_REF] focus on a special class of polynomial systems (or polynomial approximations) where only linear and cubic monomial terms are present; [START_REF] Van Loock | Bspline parameterized optimal motion trajectories for robotic systems with guaranteed constraint satisfaction[END_REF] use a similar approximation to [START_REF] Louembet | Motion planning for flat systems using positive B-splines: An LMI approach[END_REF], but extend to polynomial upper and lower bounds for the constraints.

Another noticeable contribution [START_REF] Cichella | Optimal motion planning for differentially flat systems using Bernstein approximation[END_REF]) uses a Bernstein polynomial approximation (Bézier curves) of an optimal solution for a constraint minimisation problem of a differentially flat system. This approximation also minimises the original cost function in the Bernstein space of order N ; the approximant sequence then converges to the original optimal solution. Other interesting optimisation approaches, although more specifically focused on particular steering problems, are represented by [START_REF] Mellinger | Minimum snap trajectory generation and control for quadrotors[END_REF], [START_REF] Roberts | Generating dynamically feasible trajectories for quadrotor cameras[END_REF], [START_REF] Stoican | Obstacle avoidance via B-spline parametrizations of flat trajectories[END_REF].

Particularly relevant for the present paper is the work of [START_REF] Faulwasser | Constrained reachability and trajectory generation for flat systems[END_REF]. Instead of looking for a specific trajectory, solution of an optimisation problem, the authors formulate a constrained reachability problem. More in depth, they provide some formal conditions for the existence of a suitable control function capable of steering the state from an initial set-point to a final setpoint, while satisfying constraints on both the state and the input.

In this paper, we study differentially flat systems, subject to constraints on the input, the state and a finite (arbitrary large) number of their derivatives, without the need to explicitly specify starting and ending point conditions. We do not focus on the generation of one particular trajectory resulting from an optimisation procedure. Our aim is to provide an approximate characterisation of the whole set of trajectories satisfying both the equations of the system and the constraints. The states and the inputs are written in terms of the flat outputs (and their derivatives) and the constraints are, therefore, directly remapped on the desired outputs. We assume those mappings to be expressed (natively or through approximation) in polynomial form, but some special non polynomial forms are also discussed. In order to reduce the problem to a finite dimensional one, we choose a polynomial parametrisation of the desired outputs as Bézier functions. The convex hull property of Bézier functions allows the satisfaction of the constraints on the entire time horizon without resorting to discretisation methods and post-elaborations. The characterisation of the set of feasible and constraint-compliant trajectories, is then performed in the parameter space of the control points defining the Bézier curves.

The main contribution of the paper is twofold: 1) a unified formulation of the feasible constrained trajectories encompassing a large variety of constraint types; 2) a simple and effective inner approximation of the aforementioned set of trajectories.

The present formulation allows to specify, in a unified fashion, many relevant constraints in the states (e.g. initial and final points or regions, banned regions for obstacle avoidance, corridor constraints, tubes around nominal trajectories), in the inputs (e.g. saturations, upper and lower bounds due to actuator limitations, timevarying bounds) and in arbitrary polynomial functions of states, inputs and their derivatives (e.g. limitations on the snap, jerk, curvature, energy or power consumption).

The characterisation of the set of trajectories can be performed through standard formal (like the Cylindrical Algebraic Decomposition) or numerical (Sum of Squares formulations) methods. We discuss merits and limitations of both approaches and we present a simple characterisation which, despite some restrictions, offers many interesting advantages. More in depth, our characterisation:

• provides a family of approximating sets of increasing precision converging to the whole set of feasible constrained trajectories; • can be produced by simple algebraic manipulation of the general formulation; • retains the constraint structure. If some of the bounds in the constraints change, the approximation keeps its form and just adapts some coefficients without the need for any re-computation; • allows for a simple robustness analysis w.r.t. some parameters of the constraints; • allows the extraction of a specific trajectory from the set and a fast on-line re-computation of a new solution in case of changes in the constraints.

The entire method (formulation an characterisation) is thoroughly illustrated via two interesting test cases based on a classical benchmark mechanical system. Some preliminary results have been elaborated in the PhD thesis of the third author, whose relevant chapter is available at the following arXiv preprint [START_REF] Bekcheva | Constraints on nonlinear finite dimensional flat systems[END_REF].

Notation

The set of non-negative integers is denoted by N.

Given a scalar function z ∈ C κ (R, R) and a constant α ∈ N, we denote by z α the derivative tuple of z up to the order α κ: z α = z, ż, z, . . . , z (α) . Given the vector function

v = (v 1 , . . . , v q ), v i ∈ C κ (R, R)
and the tuple α = (α 1 , . . . , α q ), α i ∈ N, we denote by v α the tuple of derivatives of each component v i of v up to its respective order α i κ:

v α = v 1 , . . . , v (α1) 1 , v 2 , . . . , v (α2) 2 , . . . , v q , . . . , v (αq) q
. We denote by R[a 1 , . . . , a n ] (or by R[a] with a = (a 1 , . . . , a n ) ∈ R n ) the ring of polynomials in the variables a 1 , . . . , a n and with coefficients in R.

Problem statement

General problem formulation

Consider the nonlinear system

ẋ(τ ) = f (x(τ ), u(τ )) (1)
with state vector x = (x 1 , . . . , x n ) and control input u = (u 1 , . . . , u m ), x i , u j ∈ C κ ([0, +∞), R) for a suitable κ ∈ N. Systems with discontinuous or non differentiable functions are beyond the scope of this paper. One can apply techniques involving Filippov solutions or the functions algebra of [START_REF] Kaldmäe | Feedback linearization of possibly non-smooth systems[END_REF] to extend our results to these cases. We assume the state, the input and their derivatives to be subject to both inequality and equality scalar constraints of the form

G i (x α x i (τ ), u α u i (τ )) 0 ∀τ ∈ I i , ∀i ∈ {1, . . . , ν in } (2a) D j (x β x j (τ ), u β u j (τ )) = 0 ∀τ ∈ I j , ∀j ∈ {1, . . . , ν eq } (2b)
with each I j being either [0, T ] (continuous constraint) or a discrete set {t 1 , . . . , t γ }, 0 t 1 • • • t γ T < +∞ (discrete constraint), and α x i , β x j ∈ N n , α u i , β u j ∈ N m . We stress that the relations (2) specify objectives (setpoints, constraints, etc.) on the finite interval [0, T ]. Objectives can be also formulated on larger intervals as a concatenation of sub-objectives on a union of subintervals, provided that some continuity and/or regularity constraints are imposed on the boundaries of each sub-interval. We here focus on just one of such intervals.

Our aim is to provide a characterisation of the set of input and state trajectories (x, u) satisfying the system dynamics (1) (usually referred to as feasible trajectories) and the constraints (2). More formally we state the following problem.

Problem 1 (Constrained trajectory set) Let C be a subalgebra of C κ ([0, T ], R). Constructively characterise the set C cons ⊆ C n+m of all extended trajectories (x, u) satisfying the system equations (1) and the constraints (2).

Problem 1 can be considered as a generalisation of a constrained reachability problem (see for instance [START_REF] Faulwasser | Constrained reachability and trajectory generation for flat systems[END_REF]). In a reachability problem the stress is usually on initial and final set-points and the goal is to find a suitable input to steer the state from the initial to the final point while, possibly, respecting the constraints. Here we wish to give a functional characterisation of the overall set of extended trajectories (x, u) satisfying some given differential constraints. A classical constrained reachability problem can be cast in the present formalism by limiting the constraints G i and D j to x and u (not their derivatives) and by forcing two of the equality constraints to coincide with the initial and final set-points.

Problem 1 is difficult to be addressed in its general setting, thus, in the following, we make some restrictions to the class of systems and to the functional space C . As a first assumption we limit the analysis to differentially flat systems [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and examples[END_REF]).

Constraints in the flat output space

Let us assume that system (1) is differentially flat with flat output

3 y = (y 1 , . . . , y m ) = h(x, u ρ u ) , (3)
with ρ u ∈ N m . The parameterisation associated to y is

x = ψ(y η x ) (4a) u = ζ(y η u ) , (4b)
with η x ∈ N n and η u ∈ N m .

Through the dynamical extension algorithm [START_REF] Li | Decoupling theory of general multivariable analytic non-linear systems[END_REF]), we get the flat output dynamics

       y (θ1) 1 = φ 1 (y µ y 1 , u µ u 1 ) . . . y (θm) m = φ m (y µ y m , u µ u m ) , (5)
with µ y i = (µ y i1 , . . . , µ y im ) ∈ N m , µ u i = (µ u i1 , . . . , µ u im ) ∈ N m and θ i > max j µ y ji . The original n-dimensional dynamics (1) and the Θ-dimensional flat output dynamics (5) (Θ = i θ i ) are in one-to-one correspondence through (3) and (4). Therefore, the constraints (2) can be re-written as

Γ i (y ω in i (τ )) 0 ∀τ ∈ I i , ∀i ∈ {1, . . . , ν in } (6a) ∆ j (y ω eq j (τ )) = 0 ∀τ ∈ I j , ∀j ∈ {1, . . . , ν eq } (6b) with Γ i (y ω in i ) = G i ((ψ(y ηx )) α x i , ζ(y η u ) α u i ), ∆ j (y ω eq j ) = D j ((ψ(y η x ) β x j , ζ(y η u ) β u j ) and ω in i , ω eq j ∈ N m .
Thus, the Problem 1 can be transformed in terms of the flat output dynamics (5) and the constraints (6) as follows.

Problem 2 (Constrained flat output set)4 Let C y be a subspace of C σ ([0, T ], R) with σ = max((θ 1 , . . . , θ m ), ω in 1 , . . . , ω in ν in , ω eq 1 , . . . , ω eq ν eq ). Constructively characterise the set C cons y ⊆ C m y of all flat outputs satisfying the dynamics (5) and the constraints (6).

Working with differentially flat systems allows us to translate, in a unified fashion, all the state and input constraints as constraints in the flat outputs and their derivatives (see ( 6)). We remark that ψ and ζ in (4) are such that ψ(y ηx ) and ζ(y ηu ) satisfy the dynamics of system (1) by construction and lead to intrinsically feasible trajectories. In other words, the extended trajectories (x, u) of ( 1) are in one-to-one correspondence with y ∈ C m y given by (3). Hence, choosing y as a solution of Problem 2 ensures that x and u given by (4) are solutions of Problem 1.

Problem specialisation -Bézier curves

For any practical purpose, one has to choose the functional space C y to which all components of the flat output belong. In order to reduce the infinite dimensional complexity, one is led to choose an approximation space C app that is dense in C σ ([0, T ], R). A possible choice is to work with parametric functions expressed in terms of basis functions like, for instance, Bernstein-Bézier, Chebychev or Spline polynomials [START_REF] Prautzsch | Bézier and B-Spline Techniques[END_REF]).

Definition 1 A Bézier curve of degree N ∈ N in the Euclidean space R r is defined as

P (τ ) = N j=0 c j B jN (τ ), τ ∈ [0, 1]
where c j ∈ R r are control points and B jN (τ ) = N j (1τ ) N -j τ j are the Bernstein polynomials5 .

We stress that a generic Bézier curve P (t) with t ∈ [a, b] and a = b can be easily written in terms of the previously defined P (τ ), τ ∈ [0, 1] by means of the affine transformation τ = (t -a)/(b -a). For the sake of simplicity, we here set T = 1 and we chose as functional space

C app N =        g ∈ C σ ([0, 1], R) | g(τ ) = N j=0 c j B jN (τ ), N ∈ N, c j ∈ R r        . (7)
The set of Bézier functions of generic degree has the very useful property of being closed with respect to addition, multiplication, degree elevation, derivation and integration operations [START_REF] Farouki | Algorithms for polynomials in Bernstein form[END_REF]). As a consequence, any integro-differential polynomial operator applied to a Bézier curve, still produces a Bézier curve (in general of different degree). Therefore, if the flat outputs y are chosen in C app N and the operators Γ i (•) and ∆ j (•) in ( 6) are integro-differential polynomials, then such constraints can be expressed, after a reordering of the coefficients and without any added conservatism, in terms of Bézier curves in

C app N for some N ∈ N which is, in general, different from N .
The main difficulty in ensuring a polynomial form for the constraints ( 6) is represented by the mappings ψ(•) and ζ(•) in the flat parametrisation (4). Indeed, while the functions G i (•) and D j (•) in (2) are defined by the designer, the mappings ψ(•) and ζ(•) are imposed by the system. However, not the entire system needs to be expressed in polynomial form w.r.t. the flat outputs, but only the states and inputs that actually enter the constraints. Moreover, in some cases this is even not necessary: the trigonometric function in the example in Section 4 is managed without resorting to a polynomial representation. Even functions expressed as ratio of polynomials can be directly tackled without any change, provided that the polynomial at the denominator is always positive (see again the example in Section 4 for the variable ω).

In the other cases, polynomial constraints can be enforced by adding some conservativeness. A common technique consists in providing a polynomial upper and lower bound for the non polynomial function (see for instance [START_REF] Van Loock | Bspline parameterized optimal motion trajectories for robotic systems with guaranteed constraint satisfaction[END_REF]). If an estimate of the flat outputs is available (or can be imposed) in the interval [0, 1] in the form of a tube around a nominal trajectory or of maximum and minimum values, then the mappings ψ(•) and ζ(•) can be approximated up to a desired precision. There exists a wide literature on approximation techniques [START_REF] Powell | Approximation Theory and Methods[END_REF]) and many efficient methods based on Chebyshev polynomials [START_REF] Rivlin | The Chebyshev Polynomials: From Approximation Theory to Algebra and Number Theory[END_REF], [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]).

In the light of the previous analysis, we assume the following.

Assumption 1 Considering y ∈ C app N y(τ ) = N j=0 a j B jN (τ ) (8)
with a j := (a j1 , . . . , a jm ) ∈ R m , the constraint in ( 6) can be written as

Γ i (y ω in i (τ )) = N in i k=0 r in ik (a)B kN in i (τ ) 0 ∀τ ∈ I i , ∀i ∈ {1, . . . , ν in } (9a) ∆ j (y ω eq j (τ )) = N eq j k=0 r eq jk (a)B kN eq j (τ ) = 0 ∀τ ∈ I j , ∀j ∈ {1, . . . , ν eq } (9b)
where

N in i , N eq j ∈ N and r in ik , r eq jk ∈ R[a] polynomials in a := (a 1 , . . . , a N ) ∈ R m(N +1) .
Remark 1 (Time-varying and tube constraints) The use of Bézier functions, or more in general of basis functions, allows an interesting extension to the constraints formulation. Instead of the constant constraint Γ i (y ω in i (τ )) 0 one can impose the more general time-

varying constraint Γ i (y ω in i (τ )) N in i k=0 c ik B kN (τ )
for some c ik ∈ R. Indeed, by defining Γi (y

ω in i (τ )) = N in i k=0 r in ik (a) -c ik B kN (τ )
, the new constraint can be cast in the form (9a) with respect to the polynomials rin ik (a) = r in ik (a) -c ik .

Tube constraints can be imposed on states and inputs in a similar fashion. Such constraints are particularly relevant in obstacle avoidance tasks, where a non feasible nominal trajectory avoiding the obstacles can be easily available via geometric planning. The goal, in this case, is to compute a feasible trajectory, still avoiding obstacles, which is not too far from the nominal one. Let us consider, for instance, a vector v = (v 1 , . . . , v p ) of state variables (or inputs) expressed as polynomials (see previous discussion) of the flat outputs and their derivatives according to (4). If a nominal trajectory v(τ ) = (v 1 (τ ), . . . , vp (τ )) is given in terms of Bézier functions as vi (τ ) =

N in i k=0 vik B kN (τ )
, then a tube constraint can be imposed as follows

p i=1 (v i (τ ) -vi (τ )) 2 c 2
, where c is the radius of the tube. It is easy to verify that, after rearranging all the terms w.r.t. the Bernstein polynomials, the previous expression is in the form (9a).

Remark 2 (Polynomial ratio constraints) The adoption of Bézier functions enables another useful extension, namely the definition of constraints described in terms of polynomial ratios. Suppose to have the following constraints: For the sake of simplicity we consider here N n = N d = N , the general case can be put in this form by adding identically zero polynomials. We can split the previous constraints as follows:

b N k=0 d k (a)B kN (τ ) - N k=0 n k (a)B kN (τ ) 0 N k=0 n k (a)B kN (τ ) -b N k=0 d k (a)B kN (τ ) 0 .
These two constraints can then be cast in the form (9a) with respect to the polynomials rin ik (a) = bd k (a) -n k (a) and rin jk (a) = n k (a) -bd k (a) for some index i, j.

In order to reduce the infinite dimensional Problem 2 to a finite dimensional one, we exploit the well known fact (see for instance [START_REF] Prautzsch | Bézier and B-Spline Techniques[END_REF]) that any Bézier curve lies within the convex hull of the points (t * j , c j ), where t * j = j N are the Greville abscissae and c j the control points. By adding some conservatism, we can replace the constraints (9), which have to hold for any τ ∈ [0, 1], with a set of constraints expressed in terms of the (finitely many) control points of (9).

Definition 2 Let

r in = (r in 10 , . . . , r in 1N in 1 , . . . , r in ν in 0 , . . . , r in ν in N in ν in )
r eq = (r eq 10 , . . . , r eq 1N eq 1 , . . . , r eq ν eq 0 , . . . , r eq ν eq N eq ν eq

)

and r = (r in , r eq ), where r in ij and r eq kh are the polynomials occurring in (9). We define the semi-algebraic set associated to the constraints (9) as

S (r, A) = a ∈ A | r in (a) 0, r eq (a) = 0 (10)
for any parallelotope

A = [a 0 , a 0 ] × • • • × [a N , a N ] (11) with a i , a i ∈ (R ∪ {-∞, +∞}) m . The sets [a i , a i ] = [a i1 , a i1 ] × • • • × [a im , a im ],
with a ij < a ij , are themselves parallelotopes bounding the elements of the vectors

a i ∈ R m .
The parallelotope A, endowed with the usual R N +1 topology, is the topological space associated to the sheaf of available trajectories among which the user is allowed to choose a reference, through the representation (8).

A is defined by the designer in terms of the upper and lower bounds of the control points of the Bézier functions defining the flat outputs. The flat outputs have almost always a physical meaning and it is not difficult to find bounds for them. The semi-algebraic set S (r, A) represents the subset of A of those trajectories fulfilling the constraints (6). Picking a in S (r, A) automatically ensures that the output y in (8) satisfies the constraints (9). Unfortunately, however, there is no way to ensure a priori that such a solution exists even in the whole parallelotope A. The infinite dimensional Problem 2 has become a finite dimensional one in the vector a, and can be reformulated as follows.

Problem 3 (Semi-algebraic set) For any fixed parallelotope A as in (11), constructively characterise the semi-algebraic set S (r, A) in (10).

Remark 3 (Reducing the conservatism) The conservatism introduced by the latter formulation of the problem is due to the, potentially large, distance of a Bézier curve from its control points. A sharp estimate of the maximum distance of a Bézier curve from its defining control polygon, i.e. the piecewise linear function connecting the points (t * j , c j ), is provided in [START_REF] Nairn | Sharp, quantitative bounds on the distance between a Bezier curve and its control polygon[END_REF]. A way to reduce such a conservatism, without increasing the number of the parameters, is to apply a suitable degree elevation of the Bézier curve. It has been proved, see e.g. [START_REF] Prautzsch | Bézier and B-Spline Techniques[END_REF], that repeated degree elevations produce a sequence of control polygons converging to the underlying Bézier curve. An estimate of the reduction of the maximum distance is again given in [START_REF] Nairn | Sharp, quantitative bounds on the distance between a Bezier curve and its control polygon[END_REF].

3 Characterisation of the semi-algebraic set

Existing approaches

A full characterisation of the semi-algebraic set S (r, A) mentioned in Problem 3 would consist of a symbolic representation of each connected component defining the set and of, at least, one internal, sample point for each component. This goal can be achieved with exact, symbolic methods like the Cylindrical Algebraic Decomposition (CAD) [START_REF] Arnon | Cylindrical algebraic decomposition I: The basic algorithm[END_REF]). The CAD is a powerful method, originally developed for the existential quantifiers elimination problem, that provides a partition of the parameter space in cells described as simpler semialgebraic sets. On each cell the polynomials defining equalities and inequalities are sign-invariant (negative, positive or zero). Among all the cells, those where all the conditions are verified are returned and for each cell a sample point is also provided.

The CAD algorithm has two main drawbacks. The first is related to its computational complexity: it is doubly exponential in the number of variables. More precisely, for s polynomials of maximum degree d in n variables, the complexity is bounded by (sd) 2 O(n) [START_REF] Davenport | Real quantifier elimination is doubly exponential[END_REF], [START_REF] England | The complexity of cylindrical algebraic decomposition with respect to polynomial degree[END_REF]). There exist many improvements and optimisations of the basic algorithm (see for instance [START_REF] Wilson | Cylindrical algebraic subdecompositions[END_REF] for sub-CAD and [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF] for a general overview of the field), which increase the practical usability of the CAD, but none of them changes its theoretical complexity. There also exists a family of algorithms based on the critical points method (see [START_REF] Grigoriev | Solving systems of polynomial inequalities in subexponential time[END_REF], [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]) whose complexity is singly exponential.

The second drawback is the huge number of cells often produced by the partition. Even if their expression is potentially much simpler than that of the original semialgebraic set, their number negatively affect the CAD usability.

To overcome these limitations, many numerical, not exact, approaches have been proposed in the literature. They aim to solve a large variety of problems including ascertaining the non-negativity of a polynomial on a basic semi-algebraic set 6 , testing emptiness of semialgebraic sets and computing inner or outer approxima-tions of basic semi-algebraic sets. They all consist in converting the exact problem in an approximated one, which can be efficiently solved by exploiting the semidefinite programming theory and Sum of Squares (SoS) polynomials [START_REF] Laurent | Sums of Squares, Moment Matrices and Optimization Over Polynomials[END_REF], [START_REF] Lasserre | A sum of squares approximation of nonnegative polynomials[END_REF], [START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF]).

In this vein, an algorithm that can be exploited to solve, at least partially, the characterisation of Problem 3 has been proposed in [START_REF] Dabbene | Simple approximations of semialgebraic sets and their applications to control[END_REF]. It is able to provide a numerical approximation (inner or outer) of a basic semi-algebraic set. Unfortunately, despite its theoretically soundness and simplicity of implementation, it often suffers from numerical problems. A way to limit such shortcomings is to add further ad-hoc constraints forcing the correct solution, but there is no systematic procedure to guide in their choice.

Moreover, all the aforementioned numerical methods have the tendency to loose the link with the symbolic constraints. That is, if a bound on a variable or a parameter on a constraint is changed, a new numerical problem has to be solved and the previous solution is of scarce (or null) use.

A simple partial characterisation

We propose here a simple numerical method to provide an approximate characterisation of the semi-algebraic set in Problem 3, which has the advantage of retaining its link with the symbolic formulation and to require few computations.

The p-norm ellipsoids

As in [START_REF] Dabbene | Simple approximations of semialgebraic sets and their applications to control[END_REF], we look for an inner approximation of a basic semi-algebraic set instead of the general set S (r, A) in Problem 3. More precisely, we assume that the equality constraints can be solved via symbolic variable substitutions or can be converted, up to a desired precision, to a couple of inequalities. In many practical problems the use of equality constraints is limited to the specification of some waypoints (such as the departure and arrival points). As a consequence, we do not deem their potential replacement by couples of inequalities a significant limitation. We also assume that the inequalities come naturally, or can be put, in the form r i r i (a) r i for some r i , r i ∈ R and r i ∈ R[a]. We stress that, in practical cases, the inequalities often come in that form and, in general, the polynomials r i always admit a finite upper and lower bound when the parallelotope A is bounded, namely when a i , a i ∈ R m . Thus we make the following assumption.

Assumption 2 There exist M polynomials r i ∈ R[a] such that the semi-algebraic set S (r, A) in (10) can be transformed in the following basic semi-algebraic set defined in terms of 2M inequalities

S (A) = {a ∈ A | r i r i (a) r i , r i , r i ∈ R, i = 1, . . . , M } . ( 12
)
We want to compute an inner approximation A ⊂ S (A). To this aim, let us define the parallelotope

P = [r 1 , r 1 ] × • • • × [r M , r M ] (13)
with r i , r i ∈ R as in (12), and the vector polynomial function R : a → (r 1 (a), . . . , r M (a)). The inverse of R is a multivalued function defined as

R -1 : ρ → {a ∈ R m(N +1) | R(a) = ρ}.
We can easily verify that S (A) = R -1 (P) ∩ A, so, for any M ⊂ P we have

S (A) ⊃ R -1 (M ) ∩ A.
The idea is to build a suitable set M ⊂ P (or a family of them) and then to define the sought inner approximation as A = R -1 (M ) ∩ A, which automatically satisfies A ⊂ S (A). In particular, if M is a semi-algebraic set defined by only one inequality, then also R -1 (M ) is a semi-algebraic set defined by one inequality.

Definition 3 Let us define the unit

M -dimensional el- lipsoid in p-norm centred in c = (c 1 , . . . , c M ) ∈ R m as E M,p = q ∈ R M M i=1 |q i -c i | p b p i -1 0 (14)
for some constants b i > 0 representing the semi-axes. If p = 2d, d ∈ N we can drop the absolute value and the ellipsoid E M,p becomes a semi-algebraic set.

In the following we impose p = 2d, d ∈ N. Moreover, if we set

2b i = r i -r i 2c i = r i + r i (15)
we have E M,p ⊂ P. In the light of the previous discussion, we can state the following proposition.

Proposition 1 Given the semi-algebraic set

R -1 (E M,p ) = a ∈ R m(N +1) M i=1 r i (a) -c i b i p -1 0 (16
) where r i are the polynomials defining S (A) in ( 12), the set A M,p = R -1 (E M,p ) ∩ A satisfies A M,p ⊂ S (A). It is thus an approximation of S (A) parametrised by M and p.

The present approximation has the advantage of being directly generated with few symbolic manipulations, once the order p has been fixed. Moreover, the positions (15) show clearly the link with the symbolic constraints.

If the bounds r i , r i in (15) change, the approximation keeps its form and just adapts some coefficients without the need for any re-computation.

Remark 4 The use of the ellipsoidal approximation provides an alternative formulation to the inequality constraints expressed as ratio of polynomials. In the Remark 2 such ratios where converted in two sets of inequalities in terms of the polynomials rin ik (a) = bd k (a) -n k (a) and rin jk (a) = n k (a)-bd k (a). Working with the control points, and accepting the ensuing conservatism, leads to the inequalities bd k (a) -n k (a) 0 and n k (a) -bd k (a) 0. These inequalities can be put back in rational form as follows:

b n k (a) d k (a) b , (17)
with the further constraints 0 d k (a) to ensure the positivity of d k . We notice that such positivity constraints can be dropped whenever the polynomials d k are squares or sum of squares. This way we have two-boundary inequalities that can be directly included in the semi-algebraic set ( 16) by setting (with a slight abuse of notation)

r k (a) = n k (a) d k (a) 2b i = b -b 2c i = b + b . ( 18
)
If a purely polynomial representation of the set R -1 (E M,p ) is desired, the inequality in ( 16) can be put in rational form through the computation of the LCD and the denominator discarded due to its positivity.

Volume coverage

Proposition 2 Given M ∈ N, the sets A M,p = R -1 (E M,p ) ∩ A constitute a family of approximating sets for S (A) parametrised in p. In other words, A M,p1 ⊆ A M,p2 for any 2 p 1 p 2 ∈ N and

lim p→∞ A M,p = S (A) .
PROOF. By definition, the ellipsoids E M,p for fixed M are such that E M,p1 ⊆ E M,p2 for any 2 p 1 p 2 ∈ N.

Moreover lim p→∞ E M,p = P. By continuity of R -1 (•) we have A M,p1 ⊆ A M,p2 for any 2 p 1 p 2 ∈ N and

lim p→∞ A M,p = lim p→∞ R -1 (E M,p )∩A = R -1 (P)∩A = S (A) .
The previous proposition ensures that, if S (A) is not empty, then there exists p 2 such that for any p p, the set A M,p is not empty too, hence the approximation A M,p provides a set of feasible solutions. The previous assertions can be made more precise by considering firstly the volume of the ellipsoid E M,p (see for instance [START_REF] Wang | Volumes of generalized unit balls[END_REF])

V (E M,p ) = 2 M Γ M 1 p + 1 Γ M p + 1 M i=1 c i
where Γ(•) is the gamma function and secondly the one of the parallelotope P

V (P) = 2 M M i=1 c i .
The ratio between these two volumes is given by

(M, p) = Γ M 1 p + 1 Γ M p + 1
and, for fixed M , clearly lim p→∞ (M, p) = 1. Unfortunately, if we fix p and we compute lim M →∞ (M, p) = 0 we verify that, for increasing dimension, the volume of a parallelotope is concentrated at the corners. Hence, the higher the dimension M , the higher p must be to capture a relevant part of the volume of P. To give an idea of how much p must increase when M increases, let us as consider a linear relation between the two: p = kM for some k ∈ N and compute the limit

(k) = lim M →∞ (M, kM ) = e -γ k Γ 1 k + 1
where γ = -Γ (1) is the Euler-Mascheroni constant.

Being (M, kM ) (k) for any M , (k) provides a lower bound on the percentage of coverage of the volume of P. For instance, (1) ≈ 0.56, (2) ≈ 0.84 and (3) ≈ 0.92.

Internal sample points

The problem with this approximation is that, even if there exists a p 2 such that for any p p, the set A M,p is not empty, there is no general rule to compute such p. Generally speaking, A M,p is an approximation of S (A), but there is no assurance of A M,p being not empty.

To complete the characterisation of the semi-algebraic set, hence to solve the Problem 3, we need to establish whether the set is not empty and, in this case, to provide at least an internal sample point. The CAD method is of no use now as the number of variables is not changed, but the fact of having the set R -1 (E M,p ) defined by only one polynomial allows us to use simpler numerical methods. More precisely, we can pose the following minimisation problem min a∈A M i=1 ri(a)-ci bi p -1 and look for a local minimiser a * . If we can find a negative minimum value, then the set is not empty and the corresponding a * is the required sample point. The aforementioned minimisation problem has constraints of a very simple form, namely only upper and lower bounds on the variables. Therefore, a simple algorithm can be employed, whose core is a standard line search method based, for instance, on the gradient descent (Nocedal and Wright [Chapt. 3, 2006]). The constraints can be accounted for by simply modifying the search direction according to the gradient projection method (Nocedal and Wright [Chapt. 16, 2006]). It is important to remark that such a minimum point can be computed quite efficiently since the gradient (and Hessian if needed) can be symbolically provided to the minimisation algorithm. Unfortunately, if the minimum value is positive, there is no assurance of the emptiness of the set as the defining polynomial has in general many minima. Multi-start approaches are beneficial in this case to improve the chance of finding the global minimum. If a negative minimum cannot be found, the order p of the polynomial can be increased and the search repeated.

Remark 5 (On-line trajectory re-computation)

The advantage of retaining the link with the symbolic constraints becomes particularly relevant when an online re-computation of a trajectory is necessary cope with a change in the constraints. If some of the bounds r i , r i in ( 15) change, the previously computed trajectory could be no longer valid; but the formal expression of the set A M,p is not changed and it can still provide feasible trajectories. A quick on-line execution of the gradient descent algorithm can provide a new internal sample point, thus a valid trajectory.

Local approximations

As mentioned before, with the generic ellipsoidal approximation there is no assurance that for a fixed p the set A M,p is not empty. This lack is due to the fact that the centre c of the ellipsoid defined in ( 15) or ( 18) is in the middle of the parallelotope P, but there is no reason for this point to be the image of a valid solution. In fact, quite often the opposite is true: R -1 (c) / ∈ S (A). On the contrary, if an internal sample point s ∈ S (A) is known, it is possible to build a non-empty set A M,p for any order p (even just 2). The trick is to centre the ellipsoidal approximation in the image of the point s along the map R. More precisely, we define the new centre and the semi-axes as 17)). The semi-axes are defined so as to ensure that E M,p ⊂ P. Actually E M,p is now a local ellipsoidal approximation of a portion of P. Notice that, the point s belongs to the ensuing set R -1 (E M,p ) by construction, hence A M,p is not empty whatever p is.

c = (c 1 , . . . , c M ) = R(s) b i = min(c i -r i , r i -c i ) (or b i = min(c i -b, b -c i ) if we consider inequalities of the type (

Robustness of solutions

An interesting aspect of the present approximation is related to the possibility of performing a robustness analysis of the solutions with respect to individual constraints. Assume that a * ∈ A M,p is a solution to Problem 3 and r is one of the constraint bounds r k , r k for any fixed k. The constants b k and c k in (15) are the only ones to be functions of r. The defining polynomial in ( 16) computed in a * is now a mono-variate function (actually a ratio of polynomials or a polynomial) of the variable r:

l a * (r) = M i=1 r i (a * ) -c i (r) b i (r) p -1 (19)
The function l a * (•) can be easily exploited to evaluate the validity of the solution a * with respect to variations of r.

It is important to remark, however, the conservativeness of such a robustness analysis. The solution a * can still be valid for a certain r even if it fails the robustness check for that r. The robustness analysis is indeed related to the set A M,p , namely the set produced by the chosen ellipsoidal approximation, not to the set S (A). If the test fails for some fixed M and p, then a * / ∈ A M,p , but it is still possible that a * ∈ S (A).

Numerical test cases

In the present section we will illustrate the steps involved in the solution of Problem 3 via two different mechanical systems.

The unicycle

A kinematic model of the unicycle mobile robot is given by

ẋ = v cos θ ẏ = v sin θ (20) θ = ω
where (x, y) represents the position of the cart with respect to an inertial frame, θ is the orientation, that is the angle between the x-axis and the heading direction measured counterclockwise. The control inputs are the linear velocity v > 0 and the angular velocity ω. The model is flat with respect to the outputs x and y. The flat parameterisation (4) is given by

θ = ψ(x, ẋ, y, ẏ) = arctan ẏ ẋ (21a) v = ζ 1 (x, ẋ, y, ẏ) = ẋ2 + ẏ2 (21b) ω = ζ 2 (x, ẋ, ẍ, y, ẏ, ÿ) = ÿ ẋ -ẍ ẏ ẋ2 + ẏ2 . (21c)
We remark that, being v > 0, both θ and ω are well defined.

Test case 1

In the first case we consider a unicycle moving from a given initial position and orientation at time t = 0 to a given final position and orientation at time t = 1. We assume limited control authorities both in v and ω. The equality and inequality constraints (2) thus read

x(0) = x 0 x(1) = x f y(0) = y 0 y(1) = y f θ(0) = θ 0 θ(1) = θ f v v v ω ω ω .
For the flat outputs we choose Bézier curves

x(t) = 3 j=0 X j B j3 (t), y(t) = 3 j=0 Y j B j3 (t) .
Bézier curves are often represented in terms of the vectors of their control points:

X = [X 0 , X 1 , X 2 , X 3 ] and Y = [Y 0 , Y 1 , Y 2 , Y 3 ].
Using this notation, we can give an idea of the transformation induced by the derivative operations on the control points:

Ẋ = 3[X 1 -X 0 , X 2 -X 1 , X 3 -X 2 ] Ẍ = 6[X 0 -2X 1 + X 2 , X 1 -2X 2 + X 3 ]
where with Ẋ and Ẍ we denoted the vectors of control points of ẋ and ẍ respectively.

The equality constraints on the initial and final positions directly fix some control points as:

X 0 = x 0 , X 3 = x f , Y 0 = y 0 and Y 3 = y f .
Imposing the constraints on the orientation is slightly more involved since the parameterisation (21a) is not polynomial. Let us consider the case θ 0 (θ f is similar). We translate the original equality constraint in an equality and an inequality constraints in terms of the control points of the flat outputs. If θ 0 = π 2 we pose ẋ(0) = 0 and ẏ(0) > 0. If θ 0 = 3π 2 we pose ẋ(0) = 0 and ẏ(0

) < 0. If θ 0 ∈ π 2 , 3π 2
we pose ẏ(0) = ẋ(0) tan θ 0 and ẋ(0) < 0. Finally, if

θ 0 ∈ 0, π 2 ∪ 3π
2 , 2π we pose ẏ(0) = ẋ(0) tan θ 0 and ẋ(0) > 0. As an example let us consider θ 0 = 0. In this case we have ẏ(0) = 0 and ẋ(0) > 0, which directly translates in terms of control points to: 3(Y 1 -Y 0 ) = 0 and 3(X 1 -X 0 ) > 0. If a two-boundary inequality of the kind considered in ( 12) is required, an upper bound on 3(X 1 -X 0 ) can be easily added either assuming the parallelotope A in ( 11) is bounded or by recalling that | ẋ(0)| < |v|.

The inequality constraints for v can be easily written as (recall v has to be positive by hypothesis) v 2 ẋ2 + ẏ2 v 2 , which has to be satisfied for each control points of v.

As a result we have a set of five quadratic two-boundary inequalities. For example, the one for the first control point of v 2 is given by:

v 2 9(X 0 -X 1 ) 2 +9(Y 0 -Y 1 ) 2 v 2 .
The inequality constraints for ω are expressed as a ratio of polynomials, so we make reference to Remark 4 to produce a set of five two-boundary inequalities of the form ω Ni Di ω, where N i and D i are the control points associated to the Bézier functions at the numerator and denominator of the ratio (21c). As an example we report here the first control points

N 0 = 18(Y 0 -Y 1 )(X 0 -2X 1 +X 2 )-18(X 0 -X 1 )(Y 0 -2Y 1 +Y 2 ) and D 0 = 9(X 0 -X 1 ) 2 + 9(Y 0 -Y 1 ) 2 .
In order to make more concrete the present test case, let us fix the following values for the constraints:

x 0 = 0 x f = 1 y 0 = 0 y f = 1/2 θ 0 = 0 θ f = 0 v = 0 v = √ 10 ω = -10 ω = 10
and the following parallelotope for the control points: A = [-1, 1] 8 . Once the equality constraints have been removed by symbolic variable substitutions, the (basic) semi-algebraic set S (r, A) in ( 10) is made up of 22 polynomial inequalities of maximal degree 2 in the two variables X 1 and X 2 .

In this simple case the CAD is able to provide an elegant solution (see the blue triangle in figure 1) to Problem 3 given by:

1 10 < X 1 < 17 20 , 1 20 + X 1 < X 2 < 9 10 .
Concerning the ellipsoidal approximation described in Section 3.2, the set S (A) in ( 12) is made up of M = 12 two-boundary inequalities (i.e. 24 single inequalities), 5 of which expressed as ratios of polynomials and 7 expressed as polynomials, all in the two variables X 1 and In order to improve the approximation of S (A) without increasing the order of the ellipsoid, one can look for a local approximation as described in Section 3.2.4. Defining the sample point s as the new centre and the semiaxes accordingly (see again Section 3.2.4), a suitable approximation of order p = 4 is found (see the green set in figure 2). It is worth noting that, despite its low order, such local approximation is even better than the general one of order p = 8.

We can perform a robustness analysis of the solution s ≈ [0.34, 0.65] with respect to variations of the bounds v, ω and ω according to the method described in section 3.2.5. We use the ellipsoidal approximation of order p = 8 and we consider two independent variations in v and in ω > 0 such that ω = -ω and ω = ω. The functions l s (v) and l s (ω) in ( 19) turn out to be a ratio of polynomials with an always positive denominator. We focus now only on the numerators n ls (v) and n ls (ω) of the two functions. n ls (ω) takes on the very simple form n ls (ω) ≈ -0.2987ω can be used, without any re-computation, even if one of the bounds v, ω and ω varies according to the computed validity margins.

Test case 2

For the second test case we show how the present approach allows to impose constraints to the state variables and control inputs in a unified fashion. This is particularly relevant when performing, for instance, obstacle avoidance.

We choose quintic Bézier curves for the x and y variables and we impose them to stay in the interior of the region bounded by two Bézier curves as depicted in figure 3. Such constraints directly translate in terms of control points for the x and y variables as follows:

X = 0, 1 10 , 1 5 , 4 5 , 9 10 , 1 Y Y Y Y = 0, - 1 10 3 , 9 10 , 1 2 , 3 8 , 1 2 Y = 0, 1 10 , 6 5 , 11 10 , 3 4 , 1 2 .
Notice that, with the present constraints, the initial and final positions of the vehicle are the same of test case 1.

Concerning the initial and final orientation we strictly fix the former and we provide a bound for the latter:

θ 0 = 0 - π 3 θ f π 3 .
The constraints on the control inputs are:

v = 0 v = 10 ω = -50 ω = 50.
Since the control points for x are completely fixed, the parallelotope A has the form:

A = 5 i=0 [Y i , Y i ].
Once the equality constraints have been removed by symbolic variable substitutions, the (basic) semialgebraic set S (r, A) in ( 10) is made up of 47 polynomial inequalities (7 of degree 0, 12 of degree 1 and 28 of degree 2) in the three variables Y 2 , Y 3 and Y 4 . In this case we did not look for a solution with the CAD algorithm, but we directly computed an ellipsoidal approximation. The set S (A) in ( 12) is made up of M = 20 twoboundary inequalities (once those of degree zero have been removed being automatically satisfied), 7 of which expressed as ratios of polynomials and 13 expressed as polynomials. All the polynomials are of maximal degree 2 in the three variables Y 2 , Y 3 and Y 4 . A numerical representation of S (A), obtained by fine gridding the parameter space, is shown in blue in figure 4.

In order to compute the approximation A M,p = R -1 (E Despite the high order (p = 80) and the high coverage index of the parallelotope P, the approximation A M,p is appreciably smaller than S (A). The reason is that the solution set is manly located in a corner of P. In order to significantly lower the degree of the approximation, one

The planar manipulator

The dynamics of a planar manipulator can be written as (see [START_REF] Murray | A mathematical introduction to robotic manipulation[END_REF]) D(q)q + C(q, q) q = τ where τ = [τ 1 , τ 2 ] are the couples at the joints and q = [θ 1 , θ 2 ] with θ 1 the angle of the first link w.r.t. the horizontal axis measured CCW and θ 2 is the angle of the second link w.r.t. the axis of the first link measured CCW. The matrices are

D(q) = α + 2βc 2 δ + βc 2 δ + βc 2 δ C(q, q) = -βs 2 θ2 -βs 2 ( θ1 + θ2 )
βs 2 θ1 0 where c 2 = cos θ 2 and s 2 = sin θ 2 . The position of the end-effector is given by

x 2 = 1 cos θ 1 + 2 cos(θ 1 + θ 2 ) y 2 = 1 sin θ 1 + 2 sin(θ 1 + θ 2 ) .
For the simulations we choose: 1 = 0.5 m, 2 = 0.3 m, α ≈ 0.18428 kg m 2 , β = 0.027 kg m 2 , δ ≈ 0.01082 kg m 2 . The system is flat w.r.t. the outputs θ 1 , θ 2 .

This second example show how the present approach allows to impose constraints to the state variables and control inputs in a unified fashion. More specifically, we want the manipulator to move from the initial position θ 1 (0) = θ1 (0) = θ1 (0) = θ2 (0) = θ2 (0) = 0, θ 2 (0) = π 6 to the final position θ 1 (1) = π 2 , θ 2 (1) = 11π 6 , θ1 (1) = θ1 (1) = θ2 (1) = θ2 (1) = 0, with the endeffector avoiding a circular obstacle of radius R = 0.15 m centred in [x o , y o ] = ( 1 + 2 )[cos π 4 , sin π 4 ]. We also impose upper and lower bounds on θ i and τ i , i = 1, 2. The inequality constraints are hence:

0 ≤ θ 1 ≤ π 2 0 ≤ θ 2 ≤ 2π |τ 1 | ≤ 30 |τ 2 | ≤ 10 and R2 ≤ (x 2 -x o ) 2 + (y 2 -y o ) 2 ≤ ( 1 + 2 ) 2 .
The expressions for x 2 , y 2 , τ 1 and τ 2 are clearly not polynomials, so we approximate the trigonometric functions with polynomials whose coefficients are found with a least square method. More precisely, we found for sin θ 1 , cos θ 1 , sin θ 2 , cos θ 2 an approximation in terms of polynomials of order 3, 4, 5, 6 respectively, with maximum absolute error 0.00236, 0.00019, 0.01436, 0.00302 respectively. We remark that the smaller the interval of the angle, the smaller the order of the approximating polynomials.

We choose a Bézier of order 5 for θ 1 and of order 8 for θ 2 . The equality constraints associated to the initial and final points produce the following vectors of control points:

Θ 1 = 0, 0, 0, π 2 , π 2 , π 2
Θ 2 = π 6 , π 6 , π 6 , X 3 , X 4 , X 5 , 11π 6 , 11π 6 , 11π 6
where X 3 , X 4 and X 5 are free parameters. As suggested in Remark 3, we reduce the conservatism of the Bézier for θ 2 by a degree elevation from 8 to 10, which produces the control vector By the constraints on θ 2 and the analysis of the vector Θ2 , the parallelotope for the free control points is chosen as A = -4π 21 , 86π 21 × -688π 105 , 131π 21 × -44π 21 , 46π 21 .

The (basic) semi-algebraic set S (r, A) in ( 10) is made up of 500 polynomial inequalities (22 of degree 0, the other 478 ranging from degree 1 to 12) in the three variables X 3 , X 4 and X 5 . In this case we did not look for a solution with the CAD algorithm, but we directly computed an ellipsoidal approximation. The set S (A) in ( 12) is made up of M = 239 two-boundary inequalities, once those of degree zero have been removed being automatically satisfied. A numerical representation of S (A), obtained by fine gridding the parameter space, is shown in blue in figure 6. In order to compute the approximation A M,p = R -1 (E M,p ) ∩ A we chose here p = 100 ( (M, p) = 0.09) obtaining the set represented in red in figure 6. An approximation of order p = 300 ( (M, p) = 0.68, not shown here) produces an almost perfect coverage of S (A). The gradient descent algorithm initialised with starting point s 0 = [1, 1, 1] produces a valid internal sample point: s ≈ [7.16, 3.05, -0.90]. It is worth noting that the initial point s 0 ∈ A but is not a solution: s 0 / ∈ S (A).

In order to significantly lower the degree of the approximation, one can look for a local approximation as discussed in Section 3.2.4. Defining s as the new centre and the semi-axes accordingly (see again Section 3.2.4), a suitable approximation of order p = 10 is found (see the green surface in figure 6).

Conclusions

n this paper we investigate the problem of characterising the set of feasible trajectories satisfying constraints on the input, the state and their derivatives. We exploit the flatness property and a parametrisation based on Bézier functions to reduce the complexity of the problem to a finite dimension one. Other than providing a unified formulation encompassing a wide variety of different types of constraints, we present a simple symbolic approximation of the set of feasible constrained trajectories.

Future works will investigate how the present characterisation can be exploited to build other relevant sets, such as for instance the actuated energy set, and to analyse some system properties (e.g. sensitivity to parameters and constraint bounds, cross-influence of physical quantities like curvature or snap on energy consumption).

b

  Nn k=0 n k (a)B kN (τ ) N d k=0 d k (a)B kN (τ ) b where b < b ∈ R and n k , d k ∈ R[a 01 , . . . , a 0m , . . . , a N 1 , . . . , a N m ]. Let us also assume N d k=0 d k (a)B kN (τ ) > 0 for every τ ∈ [0, 1].

X 2 .

 2 Two inner approximations A M,p = R -1 (E M,p ) ∩ A of S (A) are shown in figure1, one for p = 8 ( (M, p) ≈ 0.37) in orange, the other for p = 20 ( (M, p) ≈ 0.81) in red. The gradient descent algorithm initialised with starting point s 0 = [1, 10 -5 ] produces a valid internal sample point: s ≈ [0.34, 0.65] (see the black dot in figure1). It is worth noting that the initial point s 0 ∈ A but is not a solution: s 0 / ∈ S (A).

Figure 3 .

 3 Figure 3. Bounds for the x and y states in test case 2.

Figure 4 .

 4 Figure 4. Solution domains for the test case 2. The basic semi-algebraic set S (A) is in blue, the generic ellipsoidal approximation of order p = 80 is in red, the local approximation of order p = 4 computed with respect to the sample point s = [0.92, 0.93, 0.66] is in green.

Figure 5 .

 5 Figure 5. Solution domains for the test case 2. The basic semi-algebraic set S (A) is in blue, the local approximation of order p = 8 computed with respect to the sample point s = [0.95, 0.95, 0.6] is in green.

Figure 6 .

 6 Figure 6. Solution domains for the planar manipulator. The basic semi-algebraic set S (A) is in blue, the generic ellipsoidal approximation of order p = 100 is in red, the local approximation of order p = 10 computed with respect to the sample point s ≈ [7.16, 3.05, -0.90] is in green.

  Figure1. Solution domains for the test case 1. The basic semi-algebraic set S (A) is in blue, the generic ellipsoidal approximation of order p = 20 is in red and that of order p = 8 is in orange. The black dot represents the sample point s ≈ [0.34, 0.65] computed by the gradient descent algorithm.

		1
		0.9
		0.8
		0.7
		0.6
	X 2	0.5
		0.4
		0.3
		0.2
		0.1
		0
		0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
							X 1

  Figure 2. Solution domains for the test case 1. The basic semi-algebraic set S (A) is in blue, the local approximation of order p = 4 computed with respect to the sample point s = [0.34, 0.65] is in green.

		1
		0.9
		0.8
		0.7
		0.6
	X 2	0.5
		0.4
		0.3
		0.2
		0.1
		0
		0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
							X 1
	8 + 6255.66, which leads to the valid-
	ity interval ω > 3.4684. n ls (v) takes on a more com-
	plicated form of a polynomial of order 8 in the variable
	v 2 . By numerically computing the roots, one finds the
	following validity interval 1.5269 < v < 3.5055. The ro-
	bustness test ensures that the solution s ≈ [0.34, 0.65]

  M,p ) ∩ A we chose here p = 80 ( (M, p) ≈ 0.96) obtaining the set represented in red in figure4. The gradient descent algorithm initialised with starting points 0 = 1/2[Y 2 + Y 2 , Y 3 + Y 3 , Y 4 + Y 4 ] produces a valid internal sample point: s ≈ [0.92, 0.93, 0.66].It is worth noting that the initial point s 0 ∈ A but is not a solution: s 0 / ∈ S (A).

We recall that the flat output y has the same dimension m as the input vector u.

Here the max operator is applied elementwise on each vector.

For all j = 0 and j = N , and by definition of the Bernstein polynomials, BjN (0) = BjN (1) = 0. For j = 0, by convention τ j = τ 0 = 1, whatever the value of τ is; likewise, for j = N , τ N -j = 1. Thus B0N (0) = N 0 , BNN (0) = N N .

A basic semi-algebraic set is built in terms of only (weak) polynomial inequalities.