Assalé Adjé
email: assale.adje@univ-perp.fr

Pierre-Loïc Garoche
email: pierre-loic.garoche@onera.fr

Victor Magron
email: victor.magron@imag.fr

A Sums-of-Squares Extension of Policy Iterations

In order to address the imprecision often introduced by widening operators in static analysis, policy iteration based on min-computations amounts to considering the characterization of reachable value set of a program as an iterative computation of policies, starting from a postfixpoint. Computing each policy and the associated invariant relies on a sequence of numerical optimizations. While the early research efforts relied on linear programming (LP) to address linear properties of linear programs, the current state of the art is still limited to the analysis of linear programs with at most quadratic invariants, relying on semidefinite programming (SDP) solvers to compute policies, and LP solvers to refine invariants.

We propose here to extend the class of programs considered through the use of Sums-of-Squares (SOS) based optimization. Our approach enables the precise analysis of switched systems with polynomial updates and guards. The analysis presented has been implemented in Matlab and applied on existing programs coming from the system control literature, improving both the range of analyzable systems and the precision of previously handled ones.

Introduction

A wide set of critical systems software including controller systems, rely on numerical computations. Those systems range from aircraft controllers, car engine control, anti-collision systems for aircrafts or UAVs, to nuclear powerplant monitors and medical devices such a pacemakers or insulin pumps.

In all cases, the software part implements the execution of an endless loop that reads the sensor inputs, updates its internal states and controls actuators. However the analysis of such software is hardly feasible with classical static analysis tools based on linear abstractions. In fact, according to early results in control theory from Lyapunov in the 19th century, a linear system is defined as asymptotically stable iff it satisfies the Lyapunov criterion, i.e. the existence of a quadratic invariant. In this view, it is important to develop new analysis tools able to support quadratic or richer polynomial invariants.

x ∈ X in ; w h i l e ( r 0 1 ( x ) 10 and . . . and r 0 n 0 ( x ) 1 0 ) { c a s e ( r1 1 ( x ) 10 and . . . and r 1 n 1 ( x ) 1 0 ) : x = T 1 ( x ) ; c a s e . . . c a s e ( r i 1 ( x ) 10 and . . . and r i n i ( x ) 1 0 ) : x = T i ( x ) ; } (a) A one-loop program with a switch-case loop body.

Let X 0 = {y | r 0 1 (y) 1 0 and . . . and r 0 n 0 (y) 1 0} and for all cases i, X i = {y | r i 1 (y) 1 0 and . . . and r i n i (y) 1 0}. We can define the discrete-time switched system:

x 0 ∈ X in , ∀ k ∈ N, x k+1 = T (x k ), where T (x) = T i (x) if x ∈ X i ∩ X 0 (b)
The discrete-time switched system correspondence of the program.

Figure 1: One-loop programs with switch-case loop body and its representation as a switched system.

While most controllers are linear, their actual implementation is always more complex. e.g. in order to cope with safety, additional constructs such as antiwindups or saturations are introduced. Another classical approach is the use of linear parameter varying controllers (LPV) were the gains of the linear controller vary depending on conditions: this becomes piecewise polynomial controllers at the implementation level.

We are interested here in bounding the set of reachable values of such controllers using sound analyses, that is computing a sound over approximation of reachable states. We specifically focus on a class of programs larger than linear systems: constrained piecewise polynomial systems.

A loop is performed while a conjunction of polynomial inequalities 1 is satisfied. Within this loop, different polynomial updates are performed depending on conjunctions of polynomial constraints. This class of programs is represented in Fig. 1a. The program can be analyzed through its switched system representation (see Fig. 1b). In the obtained system, the conditions to switch are only governed by the state variable: at each time k, we consider the dynamics T i such that x k ∈ X i . So, the switch conditions do not depend on the time (we do not consider the time when we reach X i ) and are not defined from random variables.

From the point-of-view of the dynamical systems or control theory, to find an over-approximation of reachable values set of a program of the class described in Fig. 1a can be reduced to find a positive invariant of its switched system representation (see Fig. 1b).

Moreover, the class of switched systems where the switch conditions only depend on the state variable is classical in control theory and includes the nonlinear systems with dead-zone, saturation, resets or hysteresis.

Related Works

Reachability analysis is a long-standing problem in dynamical systems theory, especially when the system dynamics is nonlinear. In the particular case of polynomial systems, this problem has recently attracted several research efforts.

In [START_REF] Wang | Polynomial Level-Set Method for Polynomial System Reachable Set Estimation[END_REF], the authors use a method based on sublevel sets of polynomials to analyze the reachable set of a continuous time polynomial system with initial/general constraints being encoded by semialgebraic sets. Their method relies on a so-called iterative advection algorithm based on Sums-of-Squares (SOS) and semidefinite programming (SDP) to compute either inner of outer approximation of the backward reachable set, also named domain of attraction (DoA). The stability analysis of continuous-time hybrid systems with SOS certificates was investigated in [START_REF] Papachristodoulou | Robust stability analysis of nonlinear hybrid systems[END_REF]. [START_REF] Posa | Stability analysis and control of rigid-body systems with impacts and friction[END_REF] have recently applied analogous techniques to perform stability analysis and controller synthesis in the context of robotics. Other studies rely on SOS reinforcement and moment relaxations to obtain hierarchies of approximations converging to sets of interest such as the DoA in continuous-time, either from outside in [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF] or from inside in [START_REF] Korda | Inner approximations of the region of attraction for polynomial dynamical systems[END_REF]. This approach relies on a linearization of the ordinary differential equation involved in the polynomial system, based on Liouville's Equation satisfied by adequate occupation measures. This framework has been extended to hybrid systems in [START_REF] Shia | Convex computation of the reachable set for controlled polynomial hybrid systems[END_REF], as well as to synthesis of feedback controllers in [START_REF] Majumdar | Convex optimization of nonlinear feedback controllers via occupation measures[END_REF]. Liouville's Equation can also be used to approximate other sets of interest, such as the maximal controlled invariant for either discrete-or continuous-time systems (see [START_REF] Korda | Convex computation of the maximum controlled invariant set for polynomial control systems[END_REF]).

By contrast with these prior works, our approach focuses on computing over approximation of the forward reachable set of a discrete-time polynomial system with finitely many guards, thanks to an algorithm relying on SOS and template policy iterations.

Template abstractions were introduced in [SSM04] as a way to define an abstraction based on an a-priori known vector of templates, i.e. numerical expressions over the program variables. An abstract element is then defined as a vector of reals defining bounds b i over the templates p i :

p i (x 1 , . . . , x d ) ≤ b i .
Initially templates were used in the classical abstract interpretation setting to compute Kleene fixpoints with linear functions p i . Typically, the values of the bound b i increases during the fixpoint computation until convergence to a post-fixpoint.

Later in [START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF] the authors proposed to consider generalized templates including quadratic forms and compute directly the fixpoint of these template-based abstractions using numerical optimization. When considering the sub-class of linear programs, fixpoints can be computed by using a finite sequence of linear (LP) and semidefinite (SDP) optimization problems.

Two dual approaches, respectively called Max-policies and Min-policies, can be applied. Maxpolicies [START_REF] Gawlitza | Precise fixpoint computation through strategy iteration[END_REF] iterate from initial states and compute policies as relaxations through rewriting of an optimization problem (forgetting about rank conditions). Min-policies [GGTZ07, AGG10] rely on duality principle and determine a policy through the computation of a Lagrange multiplier.

precision of the analysis of boundedness for such systems. Here, we develop a policy iteration on template domains based on polynomials. The current approach improves the previous frameworks developed in [START_REF] Gaubert | Static analysis by policy iteration on relational domains[END_REF][START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF] to handle affine (or very specific) piecewise affine systems with quadratic templates and semidefinite programming. This improvement comes from the use of SOS programming to develop an SOS extension of a relaxed functional, sharing the same properties as the one defined in [START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF]. In particular, we show, as in [START_REF] Adjé | Policy iteration in finite templates domain[END_REF], that this policy iteration has the desired convergence guarantees.

Organization of the paper

The paper is structured as follows: we characterize the class of programs considered -constrained piecewise discrete-time polynomial systems -together with their collecting semantics as a least fixpoint. Then, in Section 3, we recall definitions of generalized templates, their expression as an abstract domain and the definition of the abstract transfer function. Section 4 proposes an abstraction of the transfer function using an SOS reinforcement, while Section 5 relies on this abstraction to perform policy iteration. Finally Section 6 presents experiments.

Constrained Piecewise Discrete-Time Polynomial Systems: Definition and Collecting Semantics

In this paper, we are interested in proving automatically that the set of values taken by the variables of the analyzed program is bounded. In the rest of the paper, we analyze the program by decomposing it using its dynamic system representation. The boundedness problem is thus reduced to prove that the set of all the possible trajectories of a dynamical system is bounded. Since the analyzed program has the form depicted by Figure 1a, we consider the special class of discrete-time dynamical systems introduced in Figure 1b that is:

(i) their dynamic law T is a piecewise polynomial function, and

(ii) the state variable x is constrained to live in some given basic semi-algebraic set2 .

We recall that a set is a basic semi-algebraic set if and only if it can be represented as a conjunction of strict or weak polynomial inequalities ("basic" means that no disjunction occurs). Note that T is a piecewise polynomial function with respect to a given partition, meaning that if we restrict T to be an element of the partition then T is a polynomial function.

We now give a formal definition of constrained piecewise discrete-time polynomial system (PPS for short).

First to define a PPS, we need a partition. Let I be a finite set of partition labels and X =

{X i ⊆ R d | i ∈ I} be
a partitioning, that is a given family of basic semi-algebraic sets satisfying the following:

i∈I X i = R d , ∀ i, j ∈ I, i = j ⇒ X i ∩ X j = ∅ . (1)
Each set X i of the partition corresponds to a case in the loop body of the program given in Figure 1a as the cases are assumed to be disjoint. By definition of basic semi-algebraic sets, it follows that for all i ∈ I, there exists a family of n i polynomials {r i j , j ∈ [n i ]} such that:

X i = x ∈ R d | r i j (x) 1 0 ∀ j ∈ [n i ] . ( 2
)
where 1 is either < or ≤ and [n i ] denotes the set {1, . . . , n i }. Eq. (2) matches with the program given in Figure 1a. Guards are defined as conjunctions of polynomial inequalities and thus are basic semi-algebraic sets.

The second tool needed to define a PPS is the piecewise polynomial dynamic relative to the partition. Let T : R d → R d be a piecewise polynomial function w.r.t. to the partitioning X . By definition, there exists a family of polynomials {T i : R d → R d , i ∈ I} such that for all i ∈ I:

T (x) = T i (x), ∀x ∈ X i .
(3)

Eq. (3) matches with the polynomial updates in Figure 1a. Finally, it remains to define the initialization and the set where the state variable lies. Let X in and X 0 be two basic semi-algebraic sets of R d , X in supposed to be compact, i.e. closed and bounded. The two sets can be represented as in Eq. (2) using their respective family of n in and n 0 polynomials:

X in = x ∈ R d | r in j (x) 1 0 ∀ j ∈ [n in ] and X 0 = x ∈ R d | r 0 j (x) 1 0 ∀ j ∈ [n 0 ]
, where for all j ∈ [n in ], r in j : R d → R and for all k ∈ [n 0 ], r 0 k : R d → R is a polynomial. The set X in and X 0 respectively denote the set of initial states of the program and the set which defines the loop condition in Figure 1a.

Let X be the family of sets {X i , i ∈ I} satisfying Eq. (1) and T be the family of functions {T i , i ∈ I} satisfying Eq. (3). We define the PPS associated to the quadruple (X in , X 0 , X , T ) as the system satisfying the following dynamic:

x 0 ∈ X in , and ∀ k ∈ N, if x k ∈ X 0 , x k+1 = T (x k ) . ( 4
)
In the rest of the paper, a PPS dynamical system is identified by a quadruple (X in , X 0 , X , T ).

Example 2.1 (Running example). We consider the following running example corresponding to a slightly modified version of [START_REF] Amir | Switched stability of nonlinear systems via sos-convex lyapunov functions and semidefinite programming[END_REF]Example 3]. By comparison with [AJ13, Example 3], the semialgebraic set X 1 (resp. X 2 ) is introduced to represent conditions under which we use the polynomial update T 1 (resp. T 2 ). The PPS is the quadruple (X in , X 0 , {X 1 , X 2 }, {T 1 , T 2 }), where:

X in = [-1, 1] × [-1, 1] X 0 = R 2 and X 1 = {x ∈ R 2 | -x 2 1 + 1 ≤ 0} X 2 = {x ∈ R 2 | x 2 1 -1 < 0}
and the family of functions {T 1 , T 2 }, defined as follows:

T 1 (x 1 , x 2 ) = 0.687x 1 + 0.558x 2 -0.0001x 1 x 2 -0.292x 1 + 0.773x 2 and T 2 (x 1 , x 2 ) = 0.369x 1 + 0.532x 2 -0.0001x 2 1 -1.27x 1 + 0.12x 2 -0.0001x 1 x 2
Its (discretized) reachable value set is simulated and depicted at Figure 2.

-1

-0.5 0 .5 1 -1 -0.5 0 .5 1 X 2 X 1 X 1 Figure 2: Running example simulation
We recall that our objective is to prove automatically that the set of the possible trajectories of the system is bounded. This set, also called the reachable values set or the collecting semantics of the system, is defined as follows:

R = k∈N T (k) | X 0 (X in ) (5)
where T | X 0 is the restriction of T over X 0 .The notation T (k)

| X 0 stands for composing k-times the map T | X 0 .

To prove this boundedness property, we can compute R and do some analysis to prove that R is bounded. Nevertheless, the computation of R cannot be done in general and instead, we have to compute an over-approximation of R and show that this approximation is bounded.

The usual abstract interpretation methodology to characterize and to construct this overapproximation relies on the representation of R as the smallest fixed point of a monotone map over a complete lattice. In other words, R satisfies:

R = T (R ∩ X 0 ) ∪ X in = i∈I T i (R ∩ X 0 ∩ X i ) ∪ X in .
Let us define ℘(R d ) the set of subsets of R d and introduce the function

F : ℘(R d ) → ℘(R d ) defined for all C ∈ ℘(R d ) by: F (C) = T (C ∩ X 0 ) ∪ X in = i∈I T i (C ∩ X 0 ∩ X i ) ∪ X in . (6)
Thus, R is the smallest fixed point of F and from Tarski's Theorem, since F is monotone and

℘(R d ) is a complete lattice: R = min{C ∈ ℘(R d ) | F (C) ⊆ C} . ( 7
)
Finally to compute an over-approximation of R it suffices to compute a set C such that

F (C) ⊆ C. A set C which satisfies F (C) ⊆ C is called an inductive invariant 3 .
The rest of the paper addresses the computation of a sound over-approximation of R using its definition as the smallest fixpoint of F (Eq. ( 7)).

Basic Semi-algebraic Inductive Invariants Set

An easy way to over-approximate the set of reachable values is to restrict the set of inductive invariants that we consider. We propose to restrict the class of such invariants to basic semialgebraic sets using template abstractions. A template abstraction consists of representing a given set as the intersection of sublevel sets of a-priori fixed functions depending on the state variables. Such functions are called templates. Then computing an inductive invariant in the templates domain boils down to providing, for each template p, a bound w(p) such that the intersection over the templates p of sublevel sets {x ∈ R d | p(x) ≤ w(p)} is an inductive invariant. In our context, a template is simply an a-priori fixed multivariate polynomial.

The overall method is not new and corresponds to a specialization of the templates abstraction (see [START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF][START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF]) to polynomial templates. However, in practice, the method developed in [AGG10, AGG12] is restricted to template polynomials of degree 2 (quadratic forms) and affine systems or a very restricted class of piecewise affine systems.

Next we give formal details about the polynomial template abstraction and the equations that must satisfy the template bounds vector w to generate an inductive invariant. From now on, we denote by P the set of templates and by F P, R the set of functions from P to R = R∪{-∞, +∞}.

We equip F P, R with the functional partial order

≤ F i.e. v ≤ F w iff v(p) ≤ w(p) for all p ∈ P.
Let w ∈ F P, R . The sets that we consider are of the form:

w = {x ∈ R d | p(x) ≤ w(p), ∀ p ∈ P} . ( 8
) Example 3.1. Let us define q 1 (x) = q 1 (x 1 , x 2 ) = x 2 1 , q 2 (x) = q 2 (x 1 , x 2 ) = x 2
2 and let us consider a well-chosen polynomial p of degree 6. We will explain in Subsection 5.3 how to generate automatically this template p. Let us define P run := {q 1 , q 2 , p}.

Consider the function w 0 over P run , w 0 (q 1 ) = w 0 (q 2 ) = 2.1391 and w 0 (p) = 0, the set 3. Now let us take the function w 1 over P run defined by w 1 (q 1 ) = 1.5503, w 1 (q 2 ) = 1.9501 and w 1 (p) = 0, the set

w 0 = {(x 1 , x 2 ) ∈ R 2 | x 2 1 ≤ w 0 (q 1 ), x 2 2 ≤ w 0 (q 2 ), p(x) ≤ w 0 (p)} is presented in Figure
w 1 = {(x 1 , x 2 ) ∈ R 2 | x 2 1 ≤ w 1 (q 1 ), x 2 2 ≤ w 1 (q 2 ), p(x) ≤ w 1 (p)} is also presented in Figure 3.
Since the set of black dots in Figure 3 belongs to w 0 and w 1 , we guess that w 0 and w 1 contain both the reachability set of the system from Example 2.1. To formally prove it, one way is to show that they are also inductive invariants for this system.

We restrict the class of inductive invariants to those of the form (8) and characterize the inductiveness for such sets. Since each polynomial template p is fixed, the considered variables that we handle are the template bounds w ∈ F P, R . Therefore, we need to translate the

-1.5 -1 -0.5 0 .5 1 1.5 1.46 1.24 -1.5 -1 -0.5 0 .5 1 1.5 1.46 1.39
The semialgebraic set of gray denotes the set w 0 of Example 3.1 since p(x) ≤ 0 is included in

x 2 i ≤ 2.1391, for i ∈ [1, 2].
The dark gray region denotes the semialgebraic set w 0 . Black dots are actual reachable states of R obtained by simulation.

i∈I T i (w ∩ X 0 ∩ X i ) ∪ X in ⊆ w .
By definition, w is an inductive invariant iff:

∀ p ∈ P, ∀x ∈ i∈I T i (w ∩ X 0 ∩ X i ) ∪ X in , p(x) ≤ w(p) .
Using the definition of the supremum, w is an inductive invariant iff:

∀ p ∈ P, sup x ∈ i∈I T i (w ∩ X 0 ∩ X i ) ∪ X in p(x) ≤ w(p) .
Now, let us consider p ∈ P. Using the fact that for all A, B ⊆ R d and for all functions f , sup

A∪B f = sup{sup A f, sup B f }: sup x∈ i∈I T i (w ∩X 0 ∩X i )∪X in p(x) = sup sup i∈I sup x∈T i (w ∩X 0 ∩X i ) p(x), sup x∈X in p(x) .
By definition of the image:

sup

x∈ i∈I T i (w ∩X 0 ∩X i )∪X in p(x) = sup sup i∈I sup y∈w ∩X 0 ∩X i p(T i (y)), sup x∈X in p(x) .
Next, we introduce the following notation, for all p ∈ P:

F i (w)(p) := sup x∈w ∩X i ∩X 0 p(T i (x)) and X in † (p) := sup x∈X in p(x) .
Finally, we define the function from F P, R to itself, for all w ∈ F P, R :

F (w) := sup sup i∈I F i (w), X in † .
Note that , † correspond exactly to the notations used in [START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF]. By construction we obtain the following proposition:

Proposition 3.1. Let w ∈ F P, R . Then w is an inductive invariant (i.e. F (w ) ⊆ w ) iff F (w) ≤ F w. From Prop. 3.1, inf{w ∈ F P, R n | F (w) ≤ F w} identifies the smallest inductive invariant w of the form (8).
Example 3.2. Let us consider the system defined at Example 2.1. Let us consider the same templates basis P run from Example 3.1 i.e. P run = {q 1 , q 2 , p} where q 1 (x) = x 2 1 , q 2 (x) = x 2 2 and p is a well-chosen polynomial of degree 6. Let w ∈ F P run , R . For i = 1 and the templates q 1 , we have:

F 1 (w)(q 1 ) = sup -x 2 1 +1≤0 x 2 1 ≤w(q 1 ), x 2 2 ≤w(q 2 ), p(x)≤w(p) (0.687x 1 + 0.558x 2 -0.0001x 1 x 2 ) 2 . Indeed, X 1 = {x ∈ R 2 | -x 2
1 + 1 ≤ 0} and X 0 = R 2 and the dynamics associated with X 1 is the polynomial function T 1 defined for all x ∈ R 2 by: T 1 (x) = 0.687x 1 +0.558x 2 -0.0001x 1 x 2 -0.292x 1 +0.773x 2 . Since q 1 computes the square of the first coordinates, this yields q 1 (T 1 (x)) = (0.687x

1 + 0.558x 2 - 0.0001x 1 x 2 ) 2 .
With w ∈ F P, R , computing F (w) boils down to solving a finite number of nonconvex polynomial optimization problems. General methods do not exist to solve such problems. In Section 4, we propose a method based on Sums-of-Squares (SOS) to over-approximate F (w).

SOS-based Relaxed Semantics

In this section, we introduce the relaxed functional on which we will compute a fixpoint, yielding a further over-approximation of the set R of reachable values. This relaxed functional is constructed from a Lagrange relaxation of maximization problems involved in the evaluation of F and Sums-of-Squares strengthening of polynomial nonnegativity constraints. First, we recall mandatory background related to Sums-of-Squares and their application in polynomial optimization. The interested reader is referred to [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] for more details.

Sums-of-Squares Programming

Let R[x] 2m stands for the set of polynomials of degree at most 2m and Σ[x] ⊂ R[x] be the cone of Sums-of-Squares (SOS) polynomials, that is Σ[x] := { i q 2 i , with q i ∈ R[x] }. Our work will use the simple fact that for all q ∈ Σ[x], then q(x) ≥ 0 for all x ∈ R d as the set Σ[x] contains only nonnegative polynomials. In other words, for any given polynomial q, we can strengthen the constraint of q being nonnegative into the existence of an SOS decomposition of q.

For q ∈ R[x] 2m , finding an SOS decomposition i q 2 i = q valid over R d is equivalent to solve the following matrix linear feasibility problem:

q(x) = b m (x) T Q b m (x) , ∀x ∈ R d , ( 9
)
where b m (x) := (1, x 1 , . . . , x d , x 2 1 , x 1 x 2 , . . . , x m d ) (the vector of all monomials in x up to degree m) and the Gram matrix Q, being a semidefinite positive matrix (i.e. all the eigenvalues of Q are nonnegative). The size of Q (as well as the length of b m ) is d+m d . Example 4.1. Consider the bi-variate polynomial q(x)

:= 1 + x 2 1 -2x 1 x 2 + x 2 2 . With b 1 (x) = (1, x 1 , x 2 ), one looks for a semidefinite positive matrix Q such that the polynomial equality q(x) = b 1 (x) T Q b 1 (x) holds for all x ∈ R 2 . The matrix Q =    1 0 0 0 1 -1 0 -1 1   
satisfies this equality and has three nonnegative eigenvalues, which are 0, 1, and 2, respectively associated to the three eigenvectors e 0 := (0, 1/ √ 2, 1/ √ 2) , e 1 := (1, 0, 0) and e 2 := (0, 1/ √ 2, -1/ √ 2) .

Defining the matrices L := (e 1 e 2 e 0 ) =

1 0 0 0 1 √ 2 1 √ 2 0 -1 √ 2 1 √ 2
and D = 1 0 0 0 2 0 0 0 0 , one obtains the decompo-

sition Q = L D L and the equality q(x) = (L b 1 (x)) T D (L b 1 (x)) = σ(x) = 1 + (x 1 -x 2 ) 2 , for all x ∈ R 2 .
The polynomial σ is called an SOS certificate and guarantees that q is nonnegative.

In practice, one can solve the general problem (9) by using semidefinite programming (SDP) solvers (e.g. Mosek [AA00], SDPA [YFN + 10], CSDP [START_REF] Borchers | Csdp, a c library for semidefinite programming[END_REF]). For more details about SDP, we refer the interested reader to [START_REF] Vandenberghe | Semidefinite programming[END_REF].

The SOS reinforcement of polynomial optimization problems consists of restricting polynomial nonnegativity to being an element of Σ[x]. In case of polynomial maximization problems, the SOS reinforcement boils down to computing an upper bound of the real optimal value. For example let p ∈ R[x] and consider the unconstrained polynomial maximization problem sup{p(x), x ∈ R d }.

Applying SOS reinforcement, we obtain:

sup{p(x), x ∈ R d } = inf{η | η -p(x) ≥ 0} ≤ inf{η | η -p(x) ∈ Σ[x]} . ( 10
)
Now, let p, q ∈ R[x] and consider the constrained polynomial maximization problem: sup{p

(x)|q(x) ≤ 0, x ∈ R d }. Let λ ∈ Σ[x], then: sup q(x)≤0, x∈R d p(x) ≤ sup x∈R d p(x) -λ(x) • q(x) .
Indeed, suppose q(x) ≤ 0, then -λ(x)q(x) ≥ 0 and p(x) ≤ p(x) -λ(x)q(x). Finally taking the supremum over {x ∈ R d | q(x) ≤ 0} provides the above inequality. Since sup{p(x) -λ(x) • q(x), x ∈ R d } is an unconstrained polynomial maximization problem then we apply an SOS reinforcement (as in Eq. ( 10)) and we obtain:

sup q(x)≤0, x∈R d p(x) ≤ sup x∈R d p(x) -λ(x) • q(x) ≤ inf{η | η -p -λq ∈ Σ[x]} .
Finally, note that this latter inequality is valid whatever λ ∈ Σ[x] and so we can take the infimum over λ ∈ Σ[x] which leads to:

sup q(x)≤0, x∈R d p(x) ≤ inf λ∈Σ[x] sup x∈R d p(x) -λ(x) • q(x) ≤ inf η-p-λq∈Σ[x] λ∈Σ[x] η . ( 11
)
In Eq. ( 11), λ is an SOS polynomial but to exploit linear programming solvers in policy iterations (see the fourth assertion of Prop. 5.2) we restrict λ to be a nonnegative scalar and in this case, since positive scalars are sum-of-squares polynomials of degree 0, we obtain a safe over-approximation of the right-hand-side of Eq. ( 11).

In presence of several constraints, we assign to each constraint an element σ ∈ Σ[x], and we consider the product of σ with its associated constraint and then the sum of all such products. This sum is finally added to the objective function.

The use of such SOS polynomials for constrained polynomial optimization problem can be seen as a generalization of the S-procedure from [START_REF] Yakubovich | S-procedure in nonlinear control theory[END_REF]. We refer to [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF] or [START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF] for applications in control. Note that the existence of SOS decompositions of positive polynomials over compact sets is ensured by the Putinar Positivstellensatz from [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF].

Relaxed semantics

The computation of F as a polynomial maximization problem cannot be directly performed using numerical solvers. We use the SOS reinforcement mechanisms described above to relax the computation and characterize an abstraction of F .

We still assume the knowledge of the template basis P, involving polynomials of degree at most 2m. Let us define F (P, R + ) the set of nonnegative functions over P i.e. g ∈ F (P, R + ) iff for all p ∈ P, g(p) ∈ R + . Let p ∈ P and w ∈ F P, R . Starting from the definition of F i , one obtains the following:

F i (w) (p) = sup q(x)≤w(q), ∀q∈P r i j (x)≤0, ∀ j∈[n i ] r 0 k (x)≤0, ∀ k∈[n 0 ] p(T i (x)) ≤ inf λ∈F (P,R + ) σ∈Σ[x],µ l ∈Σ[x],γ l ∈Σ[x] deg(σ)≤2m deg T i deg(µ l r i l )≤2m deg T i deg(γ l r 0 l )≤2m deg T i sup x∈R d p(T i (x)) + q∈P λ(q)(w(q) -q(x))
-

n i l=1 µ l (x)r i l (x) - n 0 l=1 γ l (x)r 0 l (x) ≤ inf λ,σ,µ l ,γ l ,η η s. t.                η -p • T i - q∈P λ(q)(w(q) -q) + n i l=1 µ l r i l + n 0 l=1 γ l r 0 l = σ , λ ∈ F (P, R + ) , σ ∈ Σ[x], µ l ∈ Σ[x], γ l ∈ Σ[x], η ∈ R , deg(σ) ≤ 2m deg T i , deg(µ l r i l ) ≤ 2m deg T i , deg(γ l r 0 l ) ≤ 2m deg T i .
(using an SOS reinforcement to remove the sup)

We denote by Σ[x] n the set of n-tuples of SOS polynomials. For clarity purpose, the dependency on i is omitted within the notations of the multipliers µ l and γ l . Moreover, let us write n i l=1 µ l r i l (resp. n 0 l=1 γ l r 0 l ) as µ, r i (resp. γ, r 0 ). Finally, we write F R i (w) (p) the over-approximation of F i (w) (p), defined as follows:

F R i (w) (p) = inf λ,σ,µ,γ,η η s. t.        η -p • T i - q∈P λ(q)(w(q) -q) + µ, r i + γ, r 0 = σ λ ∈ F (P, R + ) , σ ∈ Σ[x], µ ∈ Σ[x] n i , γ ∈ Σ[x] n 0 , η ∈ R , deg(σ) ≤ 2m deg T i , deg( µ, r i + γ, r 0 ) ≤ 2m deg T i . ( 12
)
In Equation ( 12), the notation λ is a vector of Lagrange multipliers. Each multiplier is associated with a constraint constructed from a template i.e. a constraint q(x) -w(q) ≤ 0. We also introduce the vector of SOS polynomials µ and γ. Their role is to take into account the presence of the constraints x ∈ X i and x ∈ X 0 in the computation of F R i (w) (p). Recall that X i and X 0 are basic semi-algebraic sets, then the size of the vectors µ and γ are equal to the number of polynomials defining X i and X 0 .

We conclude that, for all i ∈ I, the evaluation of F R i can be done using SOS programming, since it is reduced to solve a minimization problem with a linear objective function and linear combination of polynomials constrained to be sum-of-squares.

Note that F R i defined at Eq. ( 12) is the SOS extension of the relaxed function defined in [START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF]. Indeed, considering the special case where T i is affine, the templates p, q and the test functions r i , r 0 are quadratic, the vectors µ i and γ i are restricted to be nonnegative scalars, then F R i corresponds to the relaxed function defined in [START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF] at Eq. (3.12).

Example 4.2. We still consider the running example defined at Example 2.1 and take again the same templates basis P run of Example 3.1 composed of q 1 : x → x 2 1 and q 2 : x → x 2 2 . and a well-chosen polynomial p of degree 6. For the index of the partition i = 1. Recall that T 1 (x) =

0.687x 1 +0.558x 2 -0.0001x 1 x 2 -0.292x 1 +0.773x 2 and X 1 = {x ∈ R 2 | -x 2 1 + 1 ≤ 0} and thus r 1 1 (x) = -x 2 1 + 1. Let w ∈ F P run , R , then: F R 1 (w) (q 1 ) = inf λ,σ,µ,η η s. t.          η -(0.687x 1 + 0.558x 2 -0.0001x 1 x 2 ) 2 -λ(q 1 )(w(q 1 ) -x 2 1 ) -λ(q 2 )(w(q 2 ) -x 2 2 ) -λ(p)(w(p) -p(x)) + µ(x)(1 -x 2 1 ) = σ(x) λ ∈ F (P, R + ) , σ ∈ Σ[x], µ ∈ Σ[x], η ∈ R , deg(σ) ≤ 6, deg(µ) ≤ 6 .
In practice, one cannot find any feasible solution of degree less than 6, thus we replace the degree constraint by the more restrictive one: deg(σ) ≤ 6, deg(µ) ≤ 6.

The computation of F requires the approximation of X in † := sup{p(x), x ∈ X in }. Since X in is a basic semi-algebraic set and each template p is a polynomial, then the evaluation of X in † boils down to solving a polynomial maximization problem. Next, we use SOS reinforcement described above to over-approximate X in † with the set X in R , defined as follows:

X in R (p) := inf      η η -p + ν n in , r n in = σ 0 , η ∈ R, σ 0 ∈ Σ[x], ν in ∈ Σ[x] n in , deg(σ 0 ) ≤ 2m, deg( ν n in , r n in ) ≤ 2m      .
Thus, the value of X in R (p) is obtained by solving an SOS optimization problem. Since X in is a nonempty compact basic semi-algebraic set, this problem has a feasible solution (see the proof of [Las01, Th. 4.2]), ensuring that X in R (p) is finite valued.

Example 4.3. The initialization set

X in of Example 2.1 is [-1, 1] × [-1, 1].
It can be written as:

{(x 1 , x 2 ) ∈ R 2 | x 2 1 -1 ≤ 0, x 2 2 -1 ≤ 0}.
Then, considering the same template basis of Example 4.2 and the template q 1 :

X in R (q 1 ) := inf      η η -x 2 1 + ν n in 1 (x)(x 2 1 -1) + ν n in 2 (x)(x 2 2 -1) = σ 0 (x), η ∈ R, σ 0 ∈ Σ[x], ν in 1 , ν in 2 ∈ Σ[x], deg(σ 0 ) ≤ 6, deg( ν n in 1 ) ≤ 6, deg( ν n in 2 ) ≤ 6     
.

It is easy to see that taking for all

x ∈ R 2 , ν n in 1 (x) = 1 and for all x ∈ R 2 , ν n in 2 (x) = 0 leads to η -x 2 1 + ν n in 1 (x)(x 2 1 -1) + ν n in 2 (x)(x 2 2 -1) = η -1 = σ 0 (x)
. Thus for η = 1 and for all x ∈ R 2 , σ 0 (x) = 0, we obtain X in R (q 1 ) ≤ 1. We will see at

Prop. 4.1, that X in † ≤ F X in R . Thus, since X in † (q 1 ) = sup{x 2 1 | (x 1 , x 2 ) ∈ [-1, 1] × [-1, 1]} = 1, we conclude that 1 ≤ X in R (q 1 ) and
X in R (q 1 ) = 1.
Finally, we define the relaxed functional F R for all w ∈ F P, R for all p ∈ P as follows:

F R (w) (p) = sup sup i∈I F R i (w) (p), X in R (p) . ( 13
)
As we followed the construction proposed in Section 4.1, the relaxed functional F R provides a safe over-approximation of the abstract semantics F .

Proposition 4.1 (Safety).

The following statements hold:

1. X in † ≤ F X in R ;
2. For all i ∈ I, for all w ∈ F P, R , F i (w) ≤ F F R i (w);

3. For all w ∈ F P, R , F (w) ≤ F F R (w).

An important property that we will use to prove some results on policy iteration algorithm is the monotonicity of the relaxed functional.

Proposition 4.2 (Monotonicity).

1. For all i ∈ I,

w → F R i (w) is monotone on F P, R ; 2. The function w → F R (w) is monotone on F P, R . Proof. Let p ∈ P. For v ∈ F P, R , λ ∈ F (P, R + ), µ ∈ Σ[x] n i , γ ∈ Σ[x] n 0 and η ∈ R such that deg( µ, r i + γ, r 0 ) ≤ 2m deg T i , we define the polynomial in x, ψ λ,µ,γ,η (v) := η -p • T i - q∈P λ(q)(v(q) -q) + µ, r i + γ, r 0 . We define for v ∈ F P, R the set R(v) = {(λ, µ, γ, η) | ψ λ,µ,γ,η (v) ∈ Σ[x], λ ∈ F (P, R + ) , µ ∈ Σ[x] n i , γ ∈ Σ[x] n 0 , η ∈ R}. Now,
let us take w, w ∈ F P, R such that w ≤ F w . We have:

ψ λ,µ,γ,η (w) = η -p • T i -q∈P λ(q)(w(q) -q) + µ, r i + γ, r 0 = η -p • T i -q∈P λ(q)(w(q) -w (q) + w (q) -q) + µ, r i + γ, r 0 = η -p • T i -q∈P λ(q)(w (q) -q) + µ, r i + γ, r 0 -q∈P λ(q)(w(q) -w (q)) = ψ λ,µ,γ,η (w ) + q∈P λ(q)(w (q) -w(q)) .
Then, from w ≤ F w and the fact that λ(q) are nonnegative scalars, if ψ λ,µ,γ,η (w ) is an SOS polynomial, so is ψ λ,µ,γ,η (w) as a sum of a SOS polynomial and a nonnegative scalar. Hence, we have R(w ) ⊆ R(w). Finally, we recall that if A ⊆ B, then inf B ≤ inf A . We conclude that

F R i (w) (p) ≤ F R i (w ) (p). 2.
The mapping F R is monotone as supremum of monotone maps.

From the third assertion of Prop. 4.1, if w satisfies F R (w) ≤ F w then F (w) ≤ F w and from Prop. 3.1, w is an inductive invariant and thus R ⊆ w . This result is formulated as the following corollary.

Corollary 4.1 (Over-approximation). For all w ∈ F P, R such that F R (w) ≤ F w then R ⊆ w .

Policy Iteration in Polynomial Templates Abstract Domains

We are interested in computing the least fixpoint R R of F R , R R being an over-approximation of R (least fixpoint of F ). As for the definition of R, it can be reformulated using Tarski's theorem as the minimal post-fixpoint:

R R = min{w ∈ F P, R |F R (w) ≤ F w} .
The idea behind policy iteration is to over-approximate R R using successive iterations which are composed of

• the computation of polynomial template bounds using linear programming,

• the determination of new policies using SOS programming, until a fixpoint is reached. Policy iteration navigates in the set of post-fixpoints of F R and needs to start from a post-fixpoint w 0 known a-priori. It acts like a narrowing operator and can be interrupted at any time. For further information on policy iteration, the interested reader can consult [CGG + 05, GGTZ07].

Policies

Policy iteration can be used to compute a fixpoint of a monotone self-map defined as an infimum of a family of affine monotone self-maps. In this paper, we propose to design a policy iteration algorithm to compute a fixpoint of F R . In this subsection, we give the formal definition of policies in the context of polynomial templates and define the family of affine monotone self-maps. We do not apply the concept of policies on F R but on the functions F R i exploiting the fact that for all i ∈ I, F R i is the optimal value of a minimization problem. Policy iteration needs a selection property, that is, when an element w ∈ F P, R is given, there exists a policy which achieves the infimum. In our context, since we apply the concept of policies to F R i , it means that the minimization problem involved in the computation of F R i has an optimal solution. In our case, for w ∈ F P, R and p ∈ P, an optimal solution is a vector

(λ, σ, µ, γ) ∈ F (P, R + ) × Σ[x] × Σ[x] n i × Σ[x] n 0
such that, using (12), we obtain:

F R i (w) (p) = p • T i + q∈P λ(q)(w(q) -q) -µ, r i -γ, r 0 + σ and deg(σ) ≤ 2m deg T i , deg( µ, r i + γ, r 0 ) ≤ 2m deg T i . ( 14
)
Observe that in Eq. ( 14), F R i (w) (p) is a scalar whereas the right-hand-side is a polynomial. The equality in this equation means that this polynomial is a constant polynomial. Then we introduce the set of feasible solutions for the SOS problem

F R i (w) (p): Sol(w, i, p) = {(λ, σ, µ, γ) ∈ F (P, R + ) × Σ[x] × Σ[x] n i × Σ[x] n 0 | Eq. (14) holds} . ( 15
)
Since policy iteration algorithm can be stopped at any step and still provides a sound overapproximation, we stop the iteration when Sol(w, i, p) = ∅. Now, we are interested in the elements w ∈ F (P, R) such that Sol(w, i, p) is non-empty:

FS P, R = {w ∈ F P, R | ∀ i ∈ I, ∀ p ∈ P, Sol(w, i, p) = ∅} . ( 16
)
The notation FS P, R was introduced in [AGG12] to define the elements w ∈ F P, R satisfying Sol(w, i, p) = ∅. In [AGG12, Section 4.3], we could ensure that Sol(w, i, p) = ∅ using Slater's constraint qualification condition. In the current nonlinear setting, we cannot use the same condition, which yields a more complicated definition for FS P, R .

Finally, we can define a policy as a map which selects, for all w ∈ F S P, R , for all i ∈ I and for all p ∈ P a vector of Sol(w, i, p). More formally, we have the following definition: Definition 5.1 (Policies in the policy iteration SOS based setting). A policy is a map π : FS P, R → ((

I × P) → F (P, R + ) × Σ[x] × Σ[x] n i × Σ[x] n 0 ) such that: ∀ w ∈ FS P, R , ∀ i ∈ I, ∀ p ∈ P, π(w)(i, p) ∈ Sol(w, i, p).
We denote by Π the set of policies. For π ∈ Π, let us define π λ as the map from FS P, R to (I × P) → F (P, R + ) which associates with w ∈ F S P, R and (i, p) ∈ I × P the first element of π(w)(i, p) i.e. if π(w)(i, p) = (λ, σ, µ, γ) then π λ (w)(i, p) = λ. The equality π λ (w)(i, p) = λ means that when we perform the policy iterations algorithm, we select the vector of Lagrange multipliers λ associated with the constraints of the form q(x) ≤ w(q). The purpose of this selection is to update the value of w using the direction λ. The other coordinates composing π(w)(i, p) that is σ, µ, γ do not serve the policy iterations algorithm but are only used to take in consideration the sets X i and X 0 in the computation of F R i (w)(p). As said before, policy iteration exploits the linearity of maps when a policy is fixed. We have to define the affine maps we will use in a policy iteration step. With π ∈ Π, w ∈ FS P, R , i ∈ I and p ∈ P and λ = π λ (w)(i, p), let us define the map ϕ λ w,i,p : F P, R → R as follows:

v → ϕ λ w,i,p (v) = q∈P λ(q)v(q) + F R i (w) (p) - q∈P λ(q)w(q) . ( 17
)
Then, for π ∈ Π, we define for all w ∈ F S P, R , the map Φ π(w) w from F P, R → F P, R . Let v ∈ F P, R and p ∈ P:

Φ π(w) w (v)(p) = sup sup i∈I ϕ λ w,i,p (v), X in R (p) . ( 18
)
Example 5.1. Let us consider Example 4.2 and the function w 0 (q 1 ) = w 0 (q 2 ) = 2.1391 and w 0 (p) = 0. Then there exists two SOS polynomials µ and σ such that, for all x ∈ R d :

F R 1 (w) (q 1 ) = (0.687x 1 + 0.558x 2 -0.0001x 1 x 2 ) 2 + λ(q 1 )(2.1391 -x 2 1 ) +λ(q 2 )(2.1391 -x 2 2 ) -λ(p)p(x) -µ(x)(1 -x 2 1 ) + σ(x) = 1.5503 ,
with λ(q 1 ) = λ(q 2 ) = 0 and λ(p) = 2.0331. It means that λ, µ and σ are computed such that

(0.687x 1 + 0.558x 2 -0.0001x 1 x 2 ) 2 + λ(q 1 )(2.1391 -x 2 1 ) + λ(q 2 )(2.1391 -x 2 2 ) -λ(p)p(x) -µ(x)(1 - x 2 1 ) + σ(x) is actually a constant polynomial.
Then (λ, µ, σ) ∈ Sol(w 0 , 1, q 1 ) and we can define a policy π(w 0 ) such that π(w 0 )(1, q 1 ) = (λ, µ, σ) and thus π λ (w 0 )(1, q 1 ) = (0, 0, 2.0331). We can thus define for v ∈ F (P run , R), the affine mapping:

ϕ λ w 0 ,1,q 1 (v) = λ(q 1 )v(q 1 )+λ(q 2 )v(q 2 )+λ(p)v(p)+ F R 1 (w) (q 1 )-λ(q 1 )w(q 1 )-λ(q 2 )w(q 2 )-λ(p)w(p) = 2.1391v(p) + 1.5503.
Let us denote by F (P, R) the set of finite valued function on P i.e g ∈ F (P, R) iff g(p) ∈ R for all p ∈ P.

Proposition 5.1 (Properties of ϕ λ i,w,p ). Let π ∈ Π, w ∈ FS P, R and (i, p) ∈ I × P. Let us write λ = π λ (w)(i, p). The following properties are true:

1. ϕ λ w,i,p is affine on F (P, R) ; 2. ϕ λ w,i,p is monotone on F P, R ; 3. ∀ v ∈ F P, R , F R i (v)(p) ≤ ϕ λ w,i,p (v) ; 4. ϕ λ w,i,p (w) = F R i (w)(p) .
Proof. Let w ∈ FS P, R , i ∈ I, p ∈ P and π ∈ Π. 1. The fact that ϕ

π λ (w)(i,p) w,i,p
is affine follows readily from the definition (Eq. ( 17)).

The monotonicity of ϕ

π λ (w)(i,p) w,i,p
follows from the nonnegativity of π λ (w)(i, p).

3. Let v ∈ F P, R . Since w ∈ FS P, R , there exists (λ, σ, µ, γ) ∈ F (P, R + ) × Σ[x] × Σ[x] n i × Σ[x] n 0 such that deg(σ) ≤ 2m deg T i , deg( µ, r i + γ, r 0 ) ≤ 2m deg T i : F R i (w) (p) = p • T i + q∈P λ(q)(w(q) -q) -µ, r i -γ, r 0 + σ .
Writing λ = π λ (w)(i, p), we get:

ϕ λ w,i,p (v) = q∈P λ(q)v(q) - q∈P λ(q)w(q) + p • T i + q∈P λ(q)(w(q) -q) -µ, r i -γ, r 0 + σ = p • T i + q∈P λ(q)(v(q) -q) -µ, r i -γ, r 0 + σ .
Finally,

ϕ λ w,i,p (v) -p • T i - q∈P λ(q)(v(q) -q) + µ, r i γ, r 0 = σ , (19)
and recall that (Eq. ( 14))

F R i (w) (p) = inf λ,σ,µ,γ,η η s. t.        η -p • T i - q∈P λ(q)(w(q) -q) + µ, r i + γ, r 0 = σ λ ∈ F (P, R + ) , σ ∈ Σ[x], µ ∈ Σ[x] n i , γ ∈ Σ[x] n 0 , η ∈ R , deg(σ) ≤ 2m deg T i , deg( µ, r i + γ, r 0 ) ≤ 2m deg T i .
From Eq. ( 19), (λ, σ, µ, γ, ϕ λ w,i,p (v)) is a feasible solution of the latter minimization problem and we conclude that

F R i (v) (p) ≤ ϕ π λ (w)(i,p) w,i,p (v). 4. ϕ π λ (w)(i,p) w,i,p (w) = q∈P λ(q)w(q) + F R i (w) (p) - q∈P λ(q)w(q) = F R i (w) (p) .
The properties presented in Prop. 5.1 imply some useful properties for the maps Φ π(w) w .

Proposition 5.2 (Properties of Φ π(w) w

). Let π ∈ Π and w ∈ F S P, R . The following properties are true:

1. Φ π(w) w is monotone on F P, R ; 2. F R ≤ F Φ π(w) w ; 3. Φ π(w) w (w) = F R (w) ; 4. Suppose that the least fixpoint of Φ π(w) w is L ∈ F (P, R).
Then L can be computed as the unique optimal solution of the linear program:

inf    p ∈P v(p ) | ∀ (i, p) ∈ I × P, ϕ π λ (w)(i,p) i,w,p (v) ≤ v(p), ∀q ∈ P, X in R (q) ≤ v(q)    . ( 20
)
LP problem (20) corresponds exactly to the linear program presented in the case of quadratic templates [AGG12, Eq. 4.4].

Proof. Let π ∈ Π and w ∈ FS P, R .

1. The map Φ π(w) w is monotone as the map ϕ

π λ (w)(i,p) w,i,p
is monotone for all i ∈ I and for all p ∈ P, and the the fact that the point-wise supremum of monotone maps is also monotone.

2. Let v ∈ F P, R and let p ∈ P. Recall that:

F R (v) (p) = sup sup i∈I F R i (v) (p), X in R (p) ,
and from the third assertion of Prop. 5.1, we have for all i ∈ I,

F R i (v)(p) ≤ ϕ π λ (w)(i,p) w,i,p
, by taking the supremum over I and then the supremum with X in R (p), we obtain that F R (v)(p) ≤ Φ π(w) w (v)(p), yielding the desired result.

3. This result follows readily from the fourth assertion of Prop. 5.1 and the definition of Φ π(w) w (Eq. ( 18)). 4. By Tarski's theorem and as Φ π(w) w is monotone, Φ π(w) w has a least fixpoint in F P, R . Let L be this least fixpoint supposed to be finite valued. Now, from Tarski's theorem and the definition of Φ π(w) w , we have:

L = inf{v | Φ π(w) w (v) ≤ F v} = inf v | ∀ (i, p) ∈ I × P, ϕ π λ (w)(i,p) i,w,p (v) ≤ v(p), ∀q ∈ P, X in R (q) ≤ v(q) .
Let us suppose that there exists a feasible solution v such that q∈P v(q) < q∈P L(q). Note that since (inf{v, L}) ≤ inf{v, L}. This contradicts the minimality of L. We conclude that L is the optimal solution of Linear Program (20).

X in R ≤ F v,
Remark 1. We recall that the linear constraints in Problem (20) come from the use of the function defined at Equation (17) which is affine on the variable v. The linear forms are defined from the vector of Lagrange multipliers λ found when we solve the minimization problem involved in Equation (12). If we had allowed a vector of SOS polynomials λ as vector of Lagrange multipliers, we would obtain a set of polynomial inequalities that we would solve using SOS programming. The resulted problem would not have a feasible solution.

For example, let us consider an SOS polynomial template p, an SOS (non scalar) polynomial λ and a scalar c. Then, in this case, an analog of Problem (12) would be:

min{v(p) ∈ R | λ(x)v(p) + c ≤ v(p), ∀ x ∈ R, v(p) ≥ X in † (p)}
We assumed that p is a SOS polynomial template, implying that X in † (p) is strictly positive. Since λ(x) is a non scalar SOS polynomial and v(p) > 0, then v(p)(1 -λ(x)) -c is negative for some sufficiently large x. This proves the infeasibility of the problem.

Recall that a function g : R d → R is upper-semicontinuous at x iff for all (x n ) n∈N converging to

x, then lim sup n→+∞ g(x n ) ≤ g(x).

Proposition 5.3. Let p ∈ P. Then w → F R (w)(p) is upper-semicontinuous on FS P, R ∩ F (P, R).

Proof. Let π ∈ Π, w ∈ FS P, R ∩ F (P, R) and p ∈ P. Let i ∈ I. Let (w n ) n∈N be a sequence of elements of F (P, R) converging to w. Let λ = π λ (w)(i, p). Since ϕ λ i,w,p is affine on F (P, R), then ϕ λ i,w,p is continuous on F (P, R) and finally v → Φ π(w) w (v)(p) is continuous on F (P, R) as a finite supremum of continuous functions on F (P, R). Then from the second point of Prop. 5.2, for all

n ∈ N, F R (w n )(p) ≤ Φ π(w) w (w n )(p)
. By taking the lim sup, we obtain:

lim sup n→+∞ F R (w n )(p) ≤ lim sup n→+∞ Φ π(w) w (w n )(p) = Φ π(w) w (w)(p) = F R (v)(p).

Policy Iteration

Next, we describe the policy iteration algorithm. We suppose that we have a post-fixpoint w 0 of F R in F (P, R).

We detail step by step the algorithm presented in Figure 4. At Line 1, the algorithm is initialized and thus k = 0. At Line 4, we compute F R (w k ) using Eq. ( 13) and solve the SOS problem involved in Eq. ( 12). At Line 6, if for all i ∈ I and for p ∈ P, the SOS problem involved in Eq. ( 12) has an optimal solution, then a policy π is available and we can choose any optimal solution of SOS problem involved in Eq. (12) as policy. If an optimal solution does not exist then the algorithm stops and return w k . Now, if a policy π has been defined, the algorithm goes to Line 12 and we can define Φ π(w k ) w k following Eq. ( 18). Then, we solve LP problem (20) and define the new bound on templates w k+1 as the smallest fixpoint of Φ

π(w k ) w k
. Finally, at Line 13, k is incremented.

If for some k ∈ N, w k / ∈ FS P, R and w k-1 ∈ FS P, R then the algorithm stops and returns w k . Hence, we set for all l ≥ k, w l = w k .

Theorem 5.1 (Convergence result of the algorithm presented in Figure 4). The following statements hold:

1. For all k ∈ N, w k ∈ F (P, R) and F R (w k ) ≤ w k ;

2. The sequence (w k ) k≥0 generated by Algorithm 4 is decreasing and converges;

input : w 0 ∈ F (P, R), a post-fixpoint of F R output: a fixpoint w = F R (w) if ∀ k ∈ N, w k ∈ FS P, R

Let w

∞ = lim k→+∞ w k , then F R (w ∞ ) ≤ w ∞ . Furthermore, if for all k ∈ N, w k ∈ FS P, R and if w ∞ ∈ FS P, R then F R (w ∞ ) = w ∞ .
Proof. 1. We reason by induction. We have F R (w 0 ) ≤ w 0 and w 0 ∈ F (P, R) by assumption. Now suppose that for some k ∈ N, F R (w k ) ≤ w k and w k ∈ F (P, R). If w k / ∈ F S P, R then w l = w k for all l ≥ k and then we have proved the result. Now suppose that w k ∈ FS P, R and let us take π ∈ Π such that Φ

π(w k ) w k (w k ) = F R (w k ). From induction property Φ π(w k ) w k (w k ) ≤ w k and thus w k is a post-fixpoint of Φ π(w k ) w k belonging to F (P, R). Since every post-fixpoint of Φ π(w k ) w k is greater than X in R then least fixpoint of Φ π(w k ) w k
is finite valued and thus it is the optimal solution w k+1 of Problem (20). Moreover from the second point of Prop. 5.2, F R (w k+1 ) ≤ F Φ π(w k ) w k (w k+1 ) and since

w k+1 is the least fixpoint of Φ π(w k ) w k (w k+1 ) then F R (w k+1 ) ≤ w k+1
. This completes the proof and for all k ∈ N, w k ∈ F (P, R) and

F R (w k ) ≤ w k . 2. Let k ∈ N. If w k / ∈ FS P, R then w k+1 = w k ≤ w k . Now suppose that w k ∈ FS P, R and let π ∈ Π such that Φ π(w k ) w k (w k ) = F R (w k ), then from the third point of Prop. 5.2, Φ π(w k ) w k (w k ) = F R (w k ) ≤ w k ; the inequality results from the first assertion. Then w k is a post-fixpoint of Φ π(w k ) w k . Since w k+1 is the least fixpoint of Φ π(w k ) w k and Φ π(w k ) w k
is monotone then from Tarski's theorem w k+1 ≤ w k . From the first point, for all k ∈ N, w k ∈ F (P, R). Moreover by definition of F R , X in R ≤ F R (w k ) for all k ∈ N, then from the first point (w k ) k≥0 is lower bounded then it converges to some w ∞ .

3. If for some k, w k / ∈ FS P, R and w k-1 ∈ FS P, R , then w ∞ = w k and we have F R (w ∞ ) ≤ w ∞ from the first point. Now suppose that for all k ∈ N, w k ∈ FS P, R . Since F R is monotone then for all k ∈ N, F R (w ∞ ) ≤ F R (w k ) ≤ w k from the first point. Now taking the limit of the right-hand side, we get

F (w ∞ ) ≤ w ∞ . Now, let k ∈ N and let π ∈ Π such that Φ π(w k ) w k (w k ) = F R (w k ).
From the second point, w k+1 ≤ w k and from the monotonicity of

Φ π(w k ) w k , we have w k+1 = Φ π(w k ) w k (w k+1 ) ≤ Φ π(w k ) w k (w k ) = F R (w k ).
By taking the lim sup on k, we get w ∞ ≤ lim sup k→+∞ F R (w k ). As F R is upper-semicontinuous on FS P, R ∩ F (P, R) then, if

w ∞ ∈ FS P, R , w ∞ ≤ lim sup k→+∞ F R (w k ) ≤ F R (w ∞ ) and so w ∞ = F R (w ∞ ).

Initialization and templates choice

In Section 3, we have made the assumption that the template basis was given by an oracle. Moreover, in Algorithm 4, we suppose that we have a post-fixpoint w 0 ∈ F (P, R) of F R . Now, we give details about the templates basis choice and the computation of a post-fixpoint w 0 ≥ w ∞ . The templates basis choice relies on the computation of a template basis composed of one element. This single template is constructed by the method developed in [START_REF] Adjé | Property-based polynomial invariant generation using sums-of-squares optimization[END_REF] and is then completed using the strategy proposed in [AGM15, Ex. 9]. The single template computation also permits us to compute w 0 . Actually, the method developed in [START_REF] Adjé | Property-based polynomial invariant generation using sums-of-squares optimization[END_REF] is constructed by using the definition of being a post-fixpoint of F R . Indeed, suppose that the templates basis is constituted of one template p then w 0 is a post-fixpoint F R if and only if F R (w 0 )(p) ≤ w 0 (p). This is equivalent to:

X in R = inf{η | η -p + n in j=1 ν in j r in j ∈ Σ[x], ν in ∈ Σ[x] n in } ≤ w 0 ,
and for all i ∈ I:

F R i (w 0 ) (p) = inf λ,µ,γ,η η ≤ w 0 . s. t. η -p • T i -λ i (w 0 -p) + µ, r i + γ, r 0 ∈ Σ[x] λ i ≥ 0, µ ∈ Σ[x] n i , γ ∈ Σ[x] n 0 , η ∈ R
By definition of the infimum, it is equivalent to the existence of ν in ∈ Σ[x] and for all i ∈ I of

λ i ≥ 0, µ i ∈ Σ[x] n i , γ i ∈ Σ[x] n 0 such that: w 0 -p + n in j=1 ν in j r in j ∈ Σ[x] w 0 -p • T i -λ i (w 0 -p) + µ, r i + γ, r 0 ∈ Σ[x] . ( 21
)
Now to find a template, it suffices to find p such that Eq. ( 21) holds. However, the following two issues remain. First, without an objective function, p = 0 is a solution of Eq. (21). A workaround to avoid this trivial solution consists of optimizing a certain objective function under the constraints given in Eq. (21). In [START_REF] Adjé | Property-based polynomial invariant generation using sums-of-squares optimization[END_REF], a similar optimization procedure (Problem (13) of [START_REF] Adjé | Property-based polynomial invariant generation using sums-of-squares optimization[END_REF]) is used to prove a property of the form R ⊆ {x ∈ R d | κ(x) ≤ α}, for a given real-valued function κ. Here, we are interested in proving the boundedness of the reachable value set, which corresponds to minimize α with κ = • 2 2 . Second, finding λ i and p satisfying Eq. (21) boils down to solving a bilinear SOS problem, which is not easy to handle in practice. Thus, we fix λ i = 1 as in Lyapunov equations. We also take The semialgebraic sets denotes templates of degree 7 (blue), 8 (red), 9 (green) and 10 (black). The first box denotes the initial bounds as obtained in figure 3. The second one is the one obtained after 1 iterations. Except degree 7 that converged in 4 iterations, all others converged in 1. Degree 10 faced numerical issues and did not allow to refine the bounds without errors. The presented analysis has been applied to available examples of the control community literature: piecewise linear systems, polynomial systems, etc. We gathered the examples matching our criteria: discrete systems, possibly piecewise, at most polynomial. In all the considered cases, no common quadratic Lyapunov existed. In other words, not only the existing linear abstractions such as intervals or polyhedra would fail in computing a non trivial post-fixpoint, but also the existing analyses dedicated to digital filters such as [Fer04, GS07, AGG12, RJGF12]. The analysis has been implemented in Matlab and relies on the Mosek SDP solver [AA00], through the Yalmip [L 04] SOS front-end. Without outstanding performances, all experiments are performed within a few seconds per iteration, which makes us believe that a more serious implementation would perform better. We recall that the analysis could be interrupted at any point, still providing a safe upper bound.

We next present the examples handled by our SOS policy iteration algorithm:

Example 6.1. The following example corresponds to [START_REF] Feng | Stability analysis of piecewise discrete-time linear systems[END_REF]Ex. 2.1] and represents a piecewise linear system with 2 cases handling 3 variables. The initial set is:

X in = [-1, 1] 3 .
The set where the state-variable lies is:

X 0 = R 3 .
The sets defining the partition of the state-space are:

X 1 = {(x, y, z) ∈ R 3 |x ≤ 0}, X 2 = {(x, y, z) ∈ R 3 |x > 0} .
Finally the dynamics associated to the partition are:

T 1 (x, y, z) =   
x + 0.5y -0.3x + 0.8y 0.4z Example 6.2. We consider the example [START_REF] Feng | Stability analysis of piecewise discrete-time linear systems[END_REF]Ex. 3.3] which describes a piecewise linear system with 4 cases handling 2 variables. The initial set is:

X in = [-1, 1] 2 .
The set where the state-variable lies is:

X 0 = R 2 .
The sets defining the partition of the state-space are:

X 1 = {(x, y) ∈ R 2 |x ≤ -1}, X 2 = {(x, y) ∈ R 2 |x ∈] -1, 1] ∧ y > 0}, X 3 = {(x, y) ∈ R 2 |x ∈] -1, 1] ∧ y ≤ 0}, X 4 = {(x, y) ∈ R 2 |x > 1} .
Finally the dynamics associated to the partition are:

T 1 (x, y) = 0.9x -0.01y 0.1x + y -0.02 , T 4 (x, y) = 0.9x -0.01y 0.1x + y + 0.02 , T 2 (x, y) = T 3 (x, y) = x -0.02y 0.02x + 0.9y .

Example 6.3. The following example is the piecewise quadratic system with 2 cases handling 2 variables [START_REF] Amir | Switched stability of nonlinear systems via sos-convex lyapunov functions and semidefinite programming[END_REF]Ex. 3]. The initial set is:

X in = [-1, 1] 2 .
The set where the state-variable lies is:

X 0 = R 2 .
The sets defining the partition of the state-space are:

X 1 = {(x, y) ∈ R 2 | -x 4 + x 2 -1 ≤ 0}, X 2 = {(x, y) ∈ R 2 |x 4 -x 2 + 1 < 0} .
Finally the dynamics associated to the partition are:

T 1 (x, y) = 0.687x + 0.558y -0.0001xy -0.292x + 0.773y , T 2 (x, y) = 0.369x + 0.532y -0.0001x 2 -1.27x + 0.12y -0.0001xy .

Conclusion

We proposed an extension of policy iteration algorithms, using Sum-of-Squares programming. This extension allows to consider the wider class of disjunctive polynomial programs. In this new setting, we showed that we keep the advantage of policy iteration algorithms, while producing a sequence of increasingly safe over-approximations of the reachability set.

As future work, we plan to generalize this algorithm to programs involving non-polynomial updates, including square roots, divisions as well as transcendental functions. The computational method developed in the paper could be also generalized to other classes of nonlinear switched systems involving either random or temporal switching.

Figure 3 :

 3 Figure 3: Semialgebraic sets w for Example 3.1

  q∈P v is finite. Then we have inf{v, L} ≤ L and inf{v, L} = L. As Φ π(w) w is monotone and as v and L are feasible, we have Φ π(w) w

Figure 4 :

 4 Figure 4: SOS-based policy iteration algorithm for PPS programs.

Figure 5 :

 5 Figure 5: Different templates and associated bounds computed with Policy Iterations

1 is either the strict < or the weak (≤) comparison operator over reals.

For instance membership of the sub-level set {x ∈ R d | 1 -x 2 2 ≤ 0}, thus this does not entail boundedness of variable values.

In the dynamical systems theory, the inductive invariant sets are called positive invariant.

Acknowledgments

The research leading to these results has partly received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement nr. 306595 "STATOR", the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program "Investissement d'avenir", the RTRA/STAE BRIEFCASE project grant, the ANR projects INS-2012-007 CAFEIN, and ASTRID VORACE.

* The work was done when the author was at IRIT, Université Paul Sabatier at Toulouse and was supported by the CIMI (Centre International de Mathématiques et d'Informatique) Excellence program ANR-11-LABX-0040-CIMI within the program ANR-11-IDEX-0002-02 during a postdoctoral fellowship.

Contributions

The present paper is a followup of [START_REF] Adjé | Property-based polynomial invariant generation using sums-of-squares optimization[END_REF], in which Sums-of-Squares (SOS) programming is used to analyze properties (such as boundedness or safety) of piecewise discrete polynomial systems. The main contribution is the extension of the Min-policy iteration algorithm to improve the w 0 = 0 since p has a constant part. Finally, to obtain a template p, we solve the following SOS problem:

Let (p, w) be a solution of Problem (22). In [AGM15, Prop. 1], we proved that the set {x ∈ R d | p(x) ≤ 0} defines an inductive invariant. To complete the template basis, we use the strategy proposed in [AGM15, Ex. 9], that is, we work with the templates basis {x → x 2 i , i ∈ [d]} ∪ {p}. We thus use the inductive invariant set {p(x) ≤ 0, x 2 i ≤ w} as initialization i.e. the initial bound is w 0 (q) = w if q = p and w (0) (q) = 0 if q = p. As opposed to the approach of [AGM15], we avoid increasing the degree of polynomial p to obtain better bounds on the reachable values set.

Computational considerations

The number of (a-priori unknown) coefficients of the polynomial p (of degree 2m and d variables) appearing in Problem 22 is 2m+d d . Similarly, the number of coefficients of each σ i (resp. µ i j and γ i j ) is

Therefore, our framework is expected to be tractable when either d or m is small. As mentioned in [AGM15, Section 4], one could address bigger instances while exploiting sparsity properties of the initial system, as in [START_REF] Waki | Sums of squares and semidefinite programming relaxations for polynomial optimization problems with structured sparsity[END_REF].

Experiments

Details of the running Example.

Recall that our running example is given by the following PPS: (X in , X 0 , {X 1 , X 2 }, {T 1 , T 2 }), where:

and the functions relative to the partition {X 1 , X 2 } are:

The first step consists in constructing the template basis and compute the template p and bound w on the reachable values as a solution of Problem (22). We fix the degree of p to 6. The template p generated from Matlab is of degree 6 and is equal to

The upper bound w is equal to 2.1343. As suggested in Section 5.3, we can take the template basis P run = {p, x → x 2 1 , x → x 2 2 }. We write q 1 for x → x 2 1 and q 2 for x → x 2 2 . The basic semialgebraic {x ∈ R 2 | p(x) ≤ 0, q 1 (x) ≤ 2.1343, q 2 (x) ≤ 2.1343} is an inductive invariant and the corresponding bounds function is w 0 = (w 0 (q 1 ), w 0 (q 2 ), w 0 (p)) = (2.1343, 2.1343, 0).

As in Line 4 of Algorithm 4, we compute the image of w 0 by F R using SOS (Eq. ( 12)). We found that F R (w 0 )(q 1 ) = 1.5503, F R (w 0 )(q 2 ) = 1.9501 and F R (w 0 )(p) = 0 . Since w 0 ∈ F S P run , R , Algorithm 4 goes to Line 6 and the computation of F R (w 0 ) permits to determine a new policy π(w 0 ). The important data is the vector λ. For example, for i = 1 and the template q 1 , the vector λ is (0, 0, 2.0331). It means that we associate for each template q a weight λ(q). In the case of λ = (0, 0, 2.0331), λ(q 1 ) = 0, λ(q 2 ) = 0 and λ(p) = 2.0332. For i = 1, the template q 1 and the bound vector w 0 , the function ϕ λ w 0 ,1,q 1 (v) = 2.0331v(p) + 1.5503. To get the new invariant, Algorithm 4 goes to Line 12 and we compute a bound vector w 1 solution of Linear Program (20). In this case, it corresponds to the following LP problem: min v(q 1 ) + v(q 2 ) + v(p) 2.0331v(p)+1.5503≤v(q 1 ), 1.0429v(p)+1.2235≤v(q 2 ), 0.9535v(p)-0.0248≤v(p), (i=2) 0.4578v(p)+0.8843≤v(q 1 ), 0.2048v(p)+1.9501≤v(q 2 ), 0.9985v(p)-3.4691e-7≤v(p), (i=1) 1≤v(q 1 ), 1≤v(q 2 ), 0≤v(p),

We obtain:

w 1 (q 1 ) = 1.5503, w 1 (q 2 ) = 1.9501 and w 1 (p) = 0 . Then, we come back to Line 4 of Algorithm 4 and we compute F R (w 1 ) using the SOS program Eq. ( 12). The implemented stopping rule is F R (w k ) -w k ∞ ≤ 1e-6 and since F R (w 1 ) -w 1 ∞ ≤ 1e-6, Algorithm 4 terminates. The computed sets are presented in Figure 3, page 8. Figure 5 presents the semialgebraic sets obtained with higher dimensional templates, up to degree 10. Results are similar but could lead to different numbers of iterations depending on the degree. In case of multiple iterations, the final value is also reached at iterate 1 and is slightly modified by following iterations.

Example 6.4. The following example is the hand-crafted piecewise polynomial of degree 3 with 2 cases developed in [START_REF] Adjé | Property-based polynomial invariant generation using sums-of-squares optimization[END_REF]. The initial set is:

The set where the state-variable lies is:

The sets defining the partition of the state-space are:

Finally the dynamics associated to the partition are:

The table 1 summarizes the examples considered, the bounds obtained, the degree of the polynomial templates and the number of iterations performed before reaching the fixpoint. 12 max (10) "No good invariant" occurs when the template synthesis fails, i.e. does not provide a sound postfixpoint or some numerical issues occurs during the policy iterations phase. It seems to be due to the large size of the SOS problems together with numerical issues related to the interior point methods implemented in the relying solvers.


Damien Chablat
email: damien.chablat@irccyn.ec-nantes.fr

Stéphane Caro

Emilie Bouyer

É Bouyer

The Optimization of a Novel Prismatic Drive

Keywords: Cam design, Transmission, Hertz pressure, Slide-o-Cam Conjugate cams

Introduction

In robotic and mechatronic applications, whereby motion is controlled using a piece of software, the conversion from rotational to translational motions is usually realized by means of ball-screws or linear actuators. The both are gaining popularity. However they present some drawbacks. On the one hand, ball-screws comprise a high number of moving parts, their performance depending on the number of balls rolling in the shaft groove. Moreover, they have a low load-carrying capacity, due to the punctual contact between the balls and the groove. On the other hand, linear bearings are composed of roller-bearings to figure out the previous issue, but these devices rely on a form of direct-drive motor, which makes them expensive to produce and maintain.

A novel transmission, called Slide-o-Cam is depicted in Fig. [START_REF] Gonz A′ Lez-Palacios | The design of a novel pure-rolling transmission to convert rotational into translational motion[END_REF] and was introduced in [START_REF] Gonz A′ Lez-Palacios | The design of a novel pure-rolling transmission to convert rotational into translational motion[END_REF] to transform a rotational motion to a translational one. Slide-o-Cam is composed of four main elements: ( i ) the frame; ( ii ) the cam; ( iii ) the follower; and ( iv ) the rollers. The input axis on which the cams are mounted, named camshaft, is driven at a constant angular velocity by means of an actuator under computer-control. Power is transmitted to the output, the translating follower, which is the roller-carrying slider, by means of pure-rolling contact between the cams and the rollers. The roller comprises two components, the pin and the bearing. The bearing is mounted to one end of the pin, while the other end is press-fit into the roller-carrying slider. Consequently, the contact between the cams and rollers occurs at the outer surface of the bearing. The mechanism uses two conjugate cam-follower pairs, which alternately take over the motion transmission to ensure a positive action; the rollers are thus driven by the cams throughout a complete cycle. Therefore, the main advantages of cam-follower mechanisms with respect to the other transmissions, which transform rotation into translation are: ( i ) the lower friction; ( ii ) the higher stiffness; and ( iii ) the reduction of wear. Many researchers have devoted their efforts to study contact stress distribution and predict surface fatigue life in machine parts under different types of loading. Indeed, when two bodies with curved surfaces, for example, a cam and a roller, are pressed together, the contact is not linear but a surface. The stress occurred may generate failures such as cracks, pits, or flaking in the material. Heinrich Rudolf Hertz (1857-1894) came up with a formula to evaluate the amount of surface deformation when two surfaces (spherical, cylindrical, or planar) are pressed each other under a certain force and within their limit of elasticity.

Synthesis of Planar Cam Mechanisms

Let the x -y frame be fixed to the machine frame and the uv frame be attached to the cam, as depicted in Fig. 2. [START_REF] Gonz A′ Lez-Palacios | The design of a novel pure-rolling transmission to convert rotational into translational motion[END_REF] O is the origin of both frames, 2 O is the center of the roller, and C is the contact point between the cam and the roller. The geometric parameters are illustrated in the same figure. The notation used in this figure is based on the general notation introduced in [START_REF] González-Palacios | Cam Synthesis[END_REF], namely, ( i ) p is the pitch, i.e., the distance between the center of two rollers on the same side of the follower; ( ii ) e is the distance between the axis of the cam and the line of centers of the rollers; ( iii ) 4 a is the radius of the roller-bearing, i.e., the radius of the roller; ( iv ) ψ is the angle of rotation of the cam, the input of the mechanism;

( v ) s is the position of the center of the roller, i.e, the displacement of the follower, which is the output of the mechanism; ( vi ) μ is the pressure angle; and ( vii ) f is the force transmitted from the cam to the roller. The above parameters as well as the surface of contact on the cam, are determined by the geometric relations derived from the Aronhold-Kennedy Theorem [START_REF] Waldron | Kinematics, Dynamics, and Design of Machinery[END_REF]. As a matter of fact, when the cam makes a complete turn ( = 2 ψ π Δ ), the displacement of the roller is equal to p , the distance between two rollers on the same side of the roller-carrying slider ( = s p Δ ). Furthermore, if we consider that Fig. 3 illustrates the home configuration of the roller, the latter is below the x -axis when = 0 ψ . Therefore, (0) = /2 s p and the input-output function s is defined as follows:

( ) = 2 2 p p s ψ ψ π - (1)
The expressions of the first and second derivatives of ( ) s ψ are given by:

c c u b b a nd v b b a ψ ψ δ ψ ψ ψ δ ψ + - - - + - - (3
ψ α ψ α ψ α δ ψ α ′ - ′ + + ⎛ ⎞ - ⎜ ⎟ ′ + ⎝ ⎠ (4)
where 1

α is the directed angle between the axis of the cam and the translating direction of the follower. 1

α is positive in the counterclockwise (ccw) direction. Considering the orientation adopted for the input angle ψ and the direction defined for the output s , as depicted in Fig. 3,

1 = /2 α π - (5)
The nondimensional design parameter η is defined below and will be used extensively in what remains. = / e p η (6) The expressions of 2 b , 3 b and δ can be simplified using Eqs. ( 1), ( 2), (4a-c), ( 5) and ( 6):

2 2 2 3 = 2 = ( 2 1) ( ) 2
= arctan 2 1 p b p b π πη ψ π π ψ π δ πη -+ - ⎛ ⎞ - ⎜ ⎟ - ⎝ ⎠ (7)
From Eq. ( 7), η cannot be equal to 1/(2 ) π . That is the first constraint on η . An extended angle Δ was introduced in [START_REF] Lee | Design for Manufacturability of Speed-Reduction Cam Mechanisms[END_REF] to know whether the cam profile can be closed or not. Angle Δ is defined as a root of Eq. ( 3). In the case of Slide-o-Cam, Δ is negative, as shown in Fig. 4. Consequently, the cam profile is closed if and only if 2

ψ π Δ ≤ ≤ -Δ . O 1 x y u v C ψ (a) O 1 x y C ψ -u -v (b) O 1 x y C ψ u v (c)

Curvature of the Cam Profile

The curvature of any planar parametric curve, in terms of the Cartesian coordinates u and v , and parameterized with parameter ψ , is given by [START_REF] López-Cajún | Optimization of Cam Mechanisms[END_REF]:

2 23 / 2 ( ) ( ) ( ) ( ) = [ ( ) ( ) ] v u u v u v ψ ψ ψ ψ κ ψ ψ ′ ′′ ′ ′′ - ′ ′ + (9)
The curvature p κ of the pitch curve is given in [START_REF] Renotte | The Design of a Novel Prismatic Drive for a Three-DOF Parallel-Kinematics Machine[END_REF] as

2 2 2 3 / 2 2 [( ) 2(2 1)( 1)] = [( ) (2 1) ] p p π ψ π πη πη κ ψ π πη - + - - - + - (10)
provided that the denominator never vanishes for any value of ψ , i.e., provided that 1/(2 )

η π ≠ (
In [START_REF] Chablat | Design Strategies of Slide-o-Cam Transmission[END_REF], the authors assume that the cam profile must be fully convex. The consequence on the design parameters was > 1/ η π .

To increase the design parameter space, we accept now to have a non convex cam. But, when the cam push the roller, the sign of the local radius c , the cam profile is not valid. In the former condition is not satisfied, we have a lap-back movement in the curve as depicted in Fig. 5.

Pressure Angle

The pressure angle of cam-roller-follower mechanisms is defined as the angle between the common normal at the cam-roller contact point C and the velocity of C as a point of the follower Chablat D., Caro S., et Bouyer E., "The Optimization of a Novel Prismatic Drive", Problems of Mechanics, No 1(26), pp. 32-42, 2007. [START_REF] López-Cajún | Optimization of Cam Mechanisms[END_REF], as depicted in Fig. 2, where the pressure angle is denoted by named μ . This angle plays an important role in cam design. The smaller | | μ , the better the force transmission. In the case of high-speed operations, i.e., angular velocities of cams exceeding 50 rpm, the pressure-angle is recommended to lie advised to be smaller than 30 o .

For the case at hand, the expression for the pressure-angle of μ is given in [START_REF] López-Cajún | Optimization of Cam Mechanisms[END_REF] as:

( ) tan = ( ) s e s ψ μ ψ ′ -
Considering the expressions for s and s′ , and using the parameter η given in Eqs.( 1), (2)a and ( 6), respectively, the expression for the pressure angle becomes the pressure angle turns to be:

2 tan = n n n πη μ ψ π - -

Conjugate Cams

To reduce the pressure angle, several cams can be assembled in the same cam-shaft. We note m the number of cams and β the angle of rotation between two adjacent cams, i.e.,

2 = nm π β (14)
On the Slide-o-Cam mechanism designed in [START_REF] Gonz A′ Lez-Palacios | The design of a novel pure-rolling transmission to convert rotational into translational motion[END_REF], two conjugate cams with one lobe each and, of which = β π , were used. Figure 6 shows two cam profiles with one and two lobes. In [START_REF] Chablat | Design Strategies of Slide-o-Cam Transmission[END_REF], the authors has studies the influence of design parameters η , 4 a and n on the values of the pressure angle, the cam driving the roller and the influence of the number of conjugate cams on the maximum value of pressure angle. Here is a sum up of the corresponding results::

• Influence of parameter η : The lower η , the lower the absolute value of the pressure angle, with 1/ η π ≥ . • Influence of the radius of the roller 4 a : The lower 4 a , the lower the absolute value of the pressure angle.

• Influence of the number of lobes n : The lower n , the lower the absolute value of the pressure angle.

• Influence of the Number of Conjugate Cams: The higher the number of conjugate cams, the lower the absolute value of the pressure angle.

These results are the same if we consider a single cam with several lobes or a two-or threeconjugate-cam mechanism and stay valid for single or conjugate-cam mechanisms and is independent of the number of lobes. However, the contact stress issues and the application were not taken into account in that research work.

Physical constraints

Let us assume that the surfaces of contact are ideal, smooth and dry, with negligible friction. Two relations follow from the strength of materials. Besides, the bearing shafts are subject to shearing, whereas the camshafts are subject to shearing and bending. Consequently, we come up with the following relations: is the torque applied to the camshaft; is the maximum allowable stress inside the cam axis which cannot be exceeded in the camshaft;

is the maximum stress inside the bearing's axis which cannot be exceeded in the bearing shaft;

Contact Stresses

When two bodies with curved surfaces, for example, a cam and a roller, are pressed together, the contact is not linear but along a surface, due to the inherent material compliance, the stresses developed in the two bodies being three-dimensional. Those contact stresses generate typical failures as cracks, pits, or flaking in the surface material. Heinrich Rudolf Hertz (1857-1894) proposed some formulas to evaluate the width of the band of contact and the maximum pressure, for the case of loaded contact between two cylinders, as depicted in Fig. 7. On the one hand, the width B of the band of contact is defined below:

1 2 1 2 1 2 16 ( ) = ( ) F K K R R B L R R + + (17)
where, F is the axial load, 1 R and 2 R are the radius of the two cylinders in contact, L is the width of the cylinders, 1 K and 2 K are coefficients which characterize the materials of the two cylinders.

2 2 1 2 1 2 1 2 1 1 = , = K K E E ν ν π π - - (18)
where, In our case, the two bodies in contact are the roller and the cam. The roller is a cylinder, the cam is not. However, we can approximate the cam locally by a cylinder. Consequently, it is possible to compute the Hertz pressure with respect to ψ , which is the angle of rotation of the cam, i.e, the input of the mechanism.

For a given cam profile, the maximum value of the Hertz pressure is obtained for the minimum radius of curvature of the cam. Obviously, the Hertz pressure is a maximum when B is a minimum and F , the magnitude of the force f r transmitted by the cam, is a maximum. Consequently, the higher B and the smaller F , the lower the Hertz pressure.

B depends on several parameters, amongst them, the equivalent radius of the contact, eq R , which is R is a minimum. Therefore, to compute the maximum value of the Hertz pressure, we have to consider the lowest value of 2 R with respect to ψ . Here, we consider the active part of the cam only.

1 2 1 2 = eq R R R R R + ( 20
) 1 R is constant, since it
The minimum value of 2 R , for a two-conjugate cam mechanism, is obtained when = /n ψ π -Δ , with n the number of lobes of the cam, and Δ the extended angle. The load F is computed by

2 2 = x y F f f + (21)
where,

0 0 2 = = , = = tan tan y m y x f C F f F f p π θ θ (22)
with m C the torque of the motor; p , the pitch of the follower; and θ , the angle depicted in Fig. 2. The value of F depends on the mechanism input, ψ whereas x f , attains a minimum when = /n ψ π -Δ . Consequently, the Hertz pressure, which is proportional to / F B reaches a maximum when = /n ψ π -Δ . Table 1 presents the maximum Hertz pressure allowed for some common materials. The values are given in MPa, with stat P , the allowable pressure for a static load. It is not advised to apply more than 40% of stat P to reach an infinite fatigue life. Obviously, the maximum allowable pressure depends also on the shape of the different parts in contact. A thin part is less stiff than a thick one. For example, in our case, we can assume that a Chablat D., Caro S., et Bouyer E., "The Optimization of a Novel Prismatic Drive", Problems of Mechanics, No 1(26), pp. 32-42, 2007. multilobe cam is less stiff than a single-lobe cam. However, we will only consider here the material for the determination of the allowable pressures values.

3.

Optimization of a Slide-o-Cam

Influence of the Different Parameters on the Hertz Pressure

The maximum value of the pressure depends on several parameters, namely, the number of conjugate cams, the material of the parts in contact, the geometry of the cam, and the load applied. Therefore, we have different ways to minimize the Hertz pressure.

• Increase the number of conjugate cams, raise the number of conjugate cams, reduce the length of the active part of the cam, and decrease the minimum value of 2 R . Nevertheless, when we increase the number of conjugate cams, we only consider one contact point. As a matter of fact, we assume that when two cams can drive the rollers, the cam with the smaller absolute value of pressure angle effectively drives the follower;

• Decrease the axial load, which is possible by minimizing the motor torque, or by increasing the pitch. Moreover, since the component of F on the X -axis is low when compared with the one of the Y -axis, it is more convenient to minimize y f ;

• Choose a material with a lower Young modulus, i.e., a more compliant material, thus increasing the surface of contact, hence, decreasing the pressure. However, when the material is more compliant, its plastic domain occurs for smaller stresses.

• Decrease the radius of the cam (defined by p and e ) and the radius of the roller 4 a .

Case study

A motivation of this research work is to design a Slide-o-Cam transmission for high-speed machines. As mentioned in the introduction, this mechanism should be suitable for the Orthoglide, which is a low power machine tool, as shown in Fig. 8 [START_REF] Chablat | The Design of a Novel Prismatic Drive for a Three-DOF Parallel-Kinematics Machine[END_REF]. Here is a list of its features:

• ball screw engine torque = 1.2 N.m;

• ball screw engine velocity = 0 to 3000 rpm;

• ball screw pitch = 20 mm/turn; • axial static load = 376 N;

• stiffness = 130 N/ μ m.

We choose here to design Slide-o-Cam with only one lob on each cam and only two cams to reduce the siez of the transmission. Let us assume that the maximum stress that the shafts can support is equal to 150 MPa. The minimum diameter of the bearing shaft bear φ to transmit the load is equal to 1.8 mm and the cam shaft cam φ is equal to 3.75 mm when the pitch is equal to 20 mm/turn.

Optimization of Slide-o-Cam

When the pitch and the torque of the transmission are fixed, four parameters can be changed, cam φ , bear φ , L and the material. Here, we use either steel (respectively aluminum), of which Young modulus is equal to 2 210 000 M / Pa mm (respectively 2 69 000 M / Pa mm ). To reduce the size of the transmission, we use only two cams even if the pressure angle can be smaller by using three cams [START_REF] Chablat | Design Strategies of Slide-o-Cam Transmission[END_REF].

The pressure angle does not depend on L . So, we can represent the isovalues of μ as a function of bear φ and cam φ as depicted in Fig. 9 taking into account the constraint on the maximal value of μ , i.e. < 30 μ o and the curvature of the cam profile. The optimal design with respect to the Hertz pressure is obtained with the maximal value of L . However, if we consider the size of the transmission such a parameter has to be bounded. Figure 12 illustrates a possible shape of Slide-o-Cam to replace the ball-screws of the Orthoglide.

Summary of results and conclusions

This paper dealt with the optimization of Slide-o-Cam mechanism based on the Pressure angle and Hertz pressure. The kinematic constraints of the Orhtoglide was used to defined to pitch and the input torque but the result was general. For a given input torque and pitch, the optimal design parameters of the cam, i.e. cam φ and bear φ are computed when we consider the maximal admissible pressure angle max μ as a constraint and the strength of the material.

Conversely to previous works, the optimization of the cam parameters based on the Hertz pressure does not converge toward the results obtained by an optimization based on the pressure angle. The optimal values are not the ones obtained for a minimum radius of the roller 4 a and a minimum diameter of the cam cam φ . Further research on the fatigue of the cam and the bearing are currently carried on.

5.

Figure 1 :

 1 Figure 1: Layout of Slide-o-Cam

Figure 2 :

 2 Figure 2: Parameterization of Slide-o-Cam Figure 3: Home configuration of the mechanism

Figure 4 : 2 . 1 . 2 O

 4212 Figure 4: Orientations of the cam found when = 0 c v : (a) = ψ Δ ; (b) = ψ π ; and (c) = 2 ψ π -Δ

  11) Let c ρ and p ρ be the radii of curvature of both the cam profile and the pitch curve, respectively, and c κ the curvature of the cam profile. Since the curvature is the reciprocal of the radius of curvature, we have = , due to the definition of the pitch curve, it is apparent that Eq. (12) in terms of c κ and p κ , we obtain the curvature of the cam profile as

ρFigure 5 :

 5 Figure 5: Cam profile and local curvature of the cam In [10], the minimal value of c ρ is obtain for 2 2 2 2 2 4 4 = min

Figure 6 :

 6 Figure 6: Assembly of two cams with = 50 p , 4 = 10 a and = 9 e : (a) one lobe; (b) two lobes

  Chablat D., Caro S., et Bouyer E., "The Optimization of a Novel Prismatic Drive", Problems of Mechanics, No 1(26), pp.32-42, 2007.

Figure 7 :

 7 Figure 7: Parallel cylinders in contact and heavily loaded

Figure 8 :

 8 Figure 8: The Orthoglide (© CNRS Photothèque / CARLSON Leif)

φ

  bear[mm] φ cam[mm]

Figure 9 :τ

 9 Figure 9: Isovalues of max μ for the Orthoglide constraints

Figure 11 :

 11 Figure 11: Variation of h max P as a function of L for = 6.7 bear

  The cam profile is determined by the displacement of the contact point C around the cam. The Cartesian coordinates of this point in the uv frame take the form[START_REF] González-Palacios | Cam Synthesis[END_REF]

	Chablat D., Caro S., et Bouyer E., "The Optimization of a Novel Prismatic Drive", Problems of Mechanics, No 1(26),
	pp. 32-42, 2007.				s	( ) = /(2 ) a p nd s ψ π ′	( ) = 0 ψ ′ ′				(2)
	2 ( ) = cos	(	3	4	)cos(	) a	( ) =	2	sin	(	3	4	)sin(	)

  Chablat D., Caro S., et Bouyer E., "The Optimization of a Novel Prismatic Drive", Problems of Mechanics, No 1(26), pp.32-42, 2007.

	1 E and 2 E are the corresponding Young moduli.
	On the other hand, the Hertz pressure is defined as follows:
	h P	=	4 L B F π	(19)

1 ν , 2 ν are the Poisson ratios of the materials of the cylinders 1 and 2,

  is the radius of the roller, i.e. 1

	4 R a and 2 = R is the local radius of the R . Finally, for a given cam profile, B is a R ρ . Therefore, B depends only on 2 cam, i.e. 2 = c
	minimum when 2

Table 1 : Allowable pressures [MPa]

 1

	Material	stat P maximum	Recommended value of max P
	Stainless steel	650	260
	Improved steel	1600 to 2000	640 to 800
	Grey cast iron	400 to 700	60 to 280
	Aluminum	62.5	25 to 150
	Polyamide	25	10

Table 2 ,

 2 the maximal Hertz pressures are smaller for aluminum. However, we have to choose steel parts for the Slide-o-Cam Chablat D., Caro S., et Bouyer E., "The Optimization of a Novel Prismatic Drive", Problems of Mechanics, No 1(26), pp. 32-42, 2007.because the advised value for h max P is 150 MPa for aluminum and 800 MPa for steel.

Table 2 : Maximal Hertz pressure for steel and aluminum parts [Mpa]

 2

	L	Steel	Aluminum
	10	974	558
	20	689	394
	30	562	322
	40	487	279
	50	435	249
	60	397	228

IRCCyN: UMR n o 6597 CNRS, École Centrale de Nantes, Université de Nantes, École des Mines de Nantes

Acknowledgments This research work was mainly made by Émilie Bouyer during a internship at McGill University from April, the 15th to August, the 31st 2006. The authors would like to acknowledge the Prof. Angeles from McGill University for making the collaborated study possible. Chablat D., Caro S., et Bouyer E., "The Optimization of a Novel Prismatic Drive", Problems of Mechanics, No 1(26), pp. 32-42, 2007.

6.