Didier Henrion

Christophe Louembet

Convex inner approximations of nonconvex semialgebraic sets applied to fixed-order controller design *

Keywords: polynomials, nonconvex optimization, LMI, fixed-order controller design

We describe an elementary algorithm to build convex inner approximations of nonconvex sets. Both input and output sets are basic semialgebraic sets given as lists of defining multivariate polynomials. Even though no optimality guarantees can be given (e.g. in terms of volume maximization for bounded sets), the algorithm is designed to preserve convex boundaries as much as possible, while removing regions with concave boundaries. In particular, the algorithm leaves invariant a given convex set. The algorithm is based on Gloptipoly 3, a public-domain Matlab package solving nonconvex polynomial optimization problems with the help of convex semidefinite programming (optimization over linear matrix inequalities, or LMIs). We illustrate how the algorithm can be used to design fixed-order controllers for linear systems, following a polynomial approach.

Introduction

The set of controllers stabilizing a linear system is generally nonconvex in the parameter space, and this is an essential difficulty faced by numerical algorithms of computer-aided control system design, see e.g. [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF] and references therein. It follows from the derivation of the Routh-Hurwitz stability criterion (or its discrete-time counterpart) that the set of stabilizing controllers is real basic semialgebraic, i.e. it is the intersection of sublevel sets of given multivariate polynomials. A convex inner approximation of this nonconvex semialgebraic stability region was obtained in [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF] in the form of linear matrix inequalities (LMI) obtained from univariate polynomial positivity conditions, see also [START_REF] Karimi | Robust control of polytopic systems by convex optimization[END_REF]. Convex polytopic inner approximations were also obtained in [START_REF] Nurges | Robust pole assignment via reflection coefficients of polynomials[END_REF], for discrete-time stability, using reflection coefficients. Convex inner approximations make it possible to design stabilizing controllers with the help of convex optimization techniques, at the price of loosing optimality w.r.t. closed-loop performance criteria (H 2 norm, H ∞ norm or alike).

Generally speaking, the technical literature abounds of convex outer approximations of nonconvex semialgebraic sets. In particular, such approximations form the basis of many branch-and-bound global optimization algorithms [START_REF] Neumaier | Complete search in continuous global optimization and constraint satisfaction[END_REF]. By construction, Lasserre's hierarchy of LMI relaxations for polynomial programming is a sequence of embedded convex outer approximations which are semidefinite representable, i.e. which are obtained by projecting affine sections of the convex cone of positive semidefinite matrices, at the price of introducing lifting variables [START_REF] Henrion | Solving nonconvex optimization problems -How Glop-tiPoly is applied to problems in robust and nonlinear control[END_REF].

After some literature search, we could not locate any systematic constructive procedure to generate convex inner approximations of nonconvex semialgebraic sets, contrasting sharply with the many convex outer approximations mentioned above. In the context of fixed-order controller design, inner approximations correspond to a guarantee of stability, at the price of loosing optimality. No such stability guarantee can be ensured with outer approximations.

The main contribution of this paper is therefore an elementary algorithm, readily implementable in Matlab, that generates convex inner approximations of nonconvex sets. Both input and output sets are basic semialgebraic sets given as lists of defining multivariate polynomials. Even though no optimality guarantees can be given in terms of volume maximization for bounded sets, the algorithm is designed to preserve convex boundaries as much as possible, while removing regions with concave boundaries. In particular, the algorithm leaves invariant a given convex set. The algorithm is based on Gloptipoly 3, a public-domain Matlab package solving nonconvex polynomial optimization problems with the help of convex LMIs [START_REF] Henrion | GloptiPoly 3: moments, optimization and semidefinite programming[END_REF]. Even though the algorithm can be useful on its own, e.g. for testing convexity of semialgebraic sets, we illustrate how it can be used to design fixed-order controllers for linear systems, following a polynomial approach.

Convex inner approximation

Given a basic closed semialgebraic set

S = {x ∈ R n : p 1 (x) ≤ 0 . . . p m (x) ≤ 0} (1)
where p i are multivariate polynomials, we are interested in computing another basic closed semialgebraic set

S = {x ∈ R n : p1 (x) ≤ 0 . . . p m(x) ≤ 0} (2)
which is a valid convex inner approximation of S, in the sense that S ⊂ S.

Ideally, we would like to find the tightest possible approximation, in the sense that the complement set S\ S = {x ∈ S : x / ∈ S} is as small as possible. Mathematically we may formulate the problem as the volume minimization problem inf

S S\ S dx
but since set S is not necessarily bounded we should make sure that this integral makes sense. Moreover, computing the volume of a given semialgebraic set is a difficult task in general [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF], so we expect that optimizing such a quantity is as much as difficult. In practice, in this paper, we will content ourselves of an inner approximation that removes the nonconvex parts of the boundary and keeps the convex parts as much as possible.

Detecting nonconvexity

Before describing the method, let us recall some basics definitions on polynomials and differential geometry.

Let x ∈ R n → p i (x) ∈ R[x]
be a multivariate polynomial of total degree d. Let

g i (x) = ∂p i (x) ∂x j j=1...n ∈ R n [x]
be its gradient vector and

H i (x) = ∂ 2 p i (x) ∂x j ∂x k j,k=1...n ∈ R n×n [x]
its (symmetric) Hessian polynomial matrix. Define the optimization problem

q i = min x,y y T H i (x, y)y s.t. p i (x) = 0 p j (x) ≤ 0, j = 1 . . . m, j = i y T g i (x) = 0 y T y = 1 (3)
with global minimizers {x 1 . . . x k i } and {y 1 . . . y k i }.

Let us make the following nondegeneracy assumption on defining polynomials p i (x):

Assumption 1 There is no point x such that p i (x) and g i (x) vanish simultaneously while satisfying p j (x) ≤ 0 for j = 1, . . . , m, j = i. Since the polynomial system p i (x) = 0, g i (x) = 0, involves n + 1 equations for n unknowns, Assumption 1 is satisfied generically. In other words, in the Euclidean space of coefficients of polynomials p i (x), instances violating Assumption 1 belong to a variety of Lebesgue measure zero, and an arbitrarily small perturbation on the coefficients generates a perturbed set S ǫ satisfying Assumption 1.

Theorem 1 Under Assumption 1, polynomial level set (1) is convex if and only if q i ≥ 0 for all i = 1, . . . , m.

Proof: The boundary of set S consists of points x such that p i (x) = 0 for some i, and p j (x) ≤ 0 for j = i. In the neighborhood of such a point, consider the Taylor series

p i (x + y) = p i (x) + y T g i (x) + y T H i (x)y + O(y 3 ) (4)
where O(y 3 ) denotes terms of degree 3 or higher in entries of vector y, the local coordinates. By Assumption 1, the gradient g i (x) does not vanish along the boundary, and hence convexity of the boundary is inferred from the quadratic term in expression (4). More specifically, when y T g i (x) = 0, vector y belongs to the hyperplane tangent to S at point x. Let V be a matrix spanning this linear subspace of dimension n -1 so that y = V ŷ for some ŷ. The quadratic form y T H i (x)y = ŷT V T H i V ŷ can be diagonalised with the congruence transformation ŷ = U ȳ (Schur decomposition), and hence y

T H i (x)y = ȳT U T V T H i V U ȳT = n-1 i=1 h i (x)ȳ 2 i .
The eigenvalues h i (x), i = 1, . . . , n -1 are reciprocals of the principal curvatures of the surface. Problem (3) then amounts to finding the minimum curvature, which is non-negative when the surface is locally convex around x.

In the case of three-dimensional surfaces (n = 3), the ideas of tangent plane, local coordinates and principal curvatures used in the proof of Theorem 1 are standard notions of differential geometry, see e.g. Section 3.3. in [START_REF] Carmo | Differential geometry of curves and surfaces[END_REF] for connections between principal curvatures and eigenvalues of the local Hessian form (called the second fundamental form, once suitably normalized).

Figure 1: Hyperboloid of one sheet (white), with tangent plane (gray) at the origin, a saddle point with a tangent convex parabola (thick black) and a tangent concave hyperbola (thick black).

As an example illustrating the proof of Theorem 1, consider the hyperboloid of one sheet

S = {x ∈ R 3 : p 1 (x) = x 2
1 -x 2 2 -x 3 ≤ 0} with gradient and Hessian

g 1 (x) =   2x 1 -2x 2 -1   , H 1 (x) =   2 0 0 0 -2 0 0 0 0   .
At the origin x = 0, the tangent plane is T = {y ∈ R 3 : y 3 = 0} and p 1 (y) = 2y 2 1 -2y 2 2 is a bivariate quadratic form with eigenvalues 2 and -2, corresponding respectively to the convex parabola {x : x 2 2 + x 3 = 0} (positive curvature) and concave hyperbola {x : x 2 1 -x 3 = 0} (negative curvature), see Figure 1. Theorem 1 can be exploited in an algorithmic way to generate a convex inner approximation of a semialgebraic set.

Algorithm 1 (Convex inner approximation)

Input: Polynomials p i , i = 1 . . . m defining set S as in [START_REF] Ackermann | Robust control: the parameter space approach[END_REF]. Small nonnegative scalar ǫ.

Output: Polynomials pi , i = 1 . . . m defining set S as in (2).

Step 1: Let i = 1.

Step 2: If deg p i ≤ 1 then go to Step 5.

Step 3: If p i (x) ∈ S, solve optimization problem (3) for optimum q i and minimizers

{x 1 . . . x k }. If p i (x) / ∈ S, go to Step 5.
Step 4: If q i < 0, then select one of the minimizers x j , j = 1 . . . k i , let p m+1 = g i (x j )(xx j ) + ǫ. Then let m = m + 1, and go to step 3.

Step 5:

Let i = i + 1. If i ≤ m then go to Step 2.
Step 6:

Return pi = p i , i = 1, . . . m.
The idea behind the algorithm is as follows. At Step 3, by solving the polynomial optimization problem of Theorem 1 we identify a point of minimal curvature along algebraic varieties defining the boundary of S. If the minimal curvature is negative, then we separate the point from the set with a gradient hyperplane, and we iterate on the resulting semialgebraic set. At the end, we obtain a valid inner approximation.

Note that

Step 2 checks if the boundary is affine, in which case the minimum curvature is zero and there is no optimization problem to be solved.

The key parameter of the algorithm is the small positive scalar ǫ used at Step 4 for separating strictly a point of minimal curvature, so that the algorithm does not identify it again at the next iteration. Moreover, in Step 4, one must elect arbitrarily a minimizer. We will discuss this issue later in this paper.

Finally, as pointed out to us by a referee, the ordering of the sequence of input polynomials p i has an impact on the sequence of output polynomials pi , and especially on the size of the convex inner approximation S. However, it seems very difficult to design a priori an optimal ordering policy.

Matlab code and geometric examples

At each step of Algorithm 1 we have to solve a potentially nonconvex polynomial optimization problem. For that purpose, we use Gloptipoly 3, a public-domain Matlab package [START_REF] Henrion | GloptiPoly 3: moments, optimization and semidefinite programming[END_REF].

The methodology consists in building and solving a hierarchy of embedded linear matrix inequality (LMI) relaxations of the polynomial optimization problem, see the survey [START_REF] Laurent | Sums of squares, moment matrices and optimization over polynomials[END_REF].

The LMI problems are solved numerically with the help of any semidefinite programming solver (by default Gloptipoly 3 uses SeDuMi). Under the assumption that our original semi-algebraic set is compact, the sequence of minimizers obtained by solving the LMI relaxations is ensured to converge mononotically to the global minimum. Under the additional assumption that the global optima live on a zero-dimensional variety (i.e. there is a finite number of them), Gloptipoly 3 eventually extracts some of them (not necessarily all, but at least one) using numerical linear algebra. The LMI problems in the hierarchy have a growing number of variables and constraints, and the main issue is that we cannot predict in advance how large has to be the LMI problem to guarantee global optimality.

In practice however we observe that it is not necessary to go very deep in the hierarchy to have a numerical certificate of global optimality.

Hyperbola

Let us first with the elementary example of an unbounded nonconvex hyperbolic region

S = {x ∈ R 2 : p 1 (x) ≤ 0} with p 1 (x) = -1 + x 1 x 2 , for which optimization problem (3) reads min 2y 1 y 2 s.t. x 2 y 1 + x 1 y 2 = 0 -1 + x 1 x 2 = 0 y 2 1 + y 2 2 = 1. Necessary optimality conditions yield immediately k 1 = 2 global minimizers x 1 = √ 2 2 (1, 1), y 1 = √ 2 2 (1, -1) and x 2 = √ 2 2 (-1, -1), y 2 = √ 2
2 (-1, 1), and hence two additional (normalized) affine constraints

p 2 (x) = -2 + x 1 + x 2 and p 3 (x) = -2 -x 1 -x 2 defining the slab S = {x : p i (x) ≤ 0, i = 1, 2, 3} = {x : -2 ≤ x 1 + x 2 ≤ 2} which is indeed a valid inner approximation of S.

Egg quartic

Now we show that Algorithm 1 can be used to detect convexity of a semialgebraic set. Consider the smooth quartic sublevel set

S = {x ∈ R 2 : p 1 (x) = x 4 1 + x 4 2 + x 2 1 + x 2 ≤ 0} represented on Figure 2. Assumption 1 is ensured since the gradient g 1 (x) = [2x 1 (x 2
1 + 2) 4x 3 2 + 1] cannot vanish for real x. A Matlab implementation of the first steps of the algorithm can be easily written using Gloptipoly 3: Notice that we solve the LMI relaxation of order 3 (e.g. moments of degree 6) of problem [START_REF] Helton | Sufficient and necessary conditions for semidefinite representability of convex hulls and sets[END_REF]. In GloptiPoly, an LMI relaxation is solved with the command msol which returns two output arguments: status and obj. Argument status can take the following values:

% problem data x 1 x 2 -1 -0.
• -1 if the LMI relaxation is infeasible or could not be solved for numerical reasons;

• 0 if the LMI relaxation could be solved but it is impossible to detect global optimality and to extract global optimizers, in which case obj is a lower (resp. upper) bound on the global minimum (resp. maximum) of the original optimization problem;

• +1 if the LMI relaxation could be solved, global optimality is certified and global minimizers are extracted, in which case obj is the global optimum of the original optimization problem.

Running the above script, Gloptipoly returns obj = 2.0000 and status = 1, certifying that the minimal curvature is strictly positive, and hence that the polynomial sublevel set is convex.

Note that in this simple case, convexity of set S follows directly from positive semidefiniteness of the Hessian

H 1 (x) = diag (12x 2 1 + 2, 12x 2 2
), yet Algorithm 1 can systematically detect convexity in more complicated cases.

Waterdrop quartic

Consider the quartic S = {x ∈ R 2 : p 1 (x) = x 4 1 + x 4 2 + x 2 1 + x 3 2 ≤
0} which has a singular point at the origin, hence violating Assumption 1.

Applying Algorithm 1, the LMI relaxation of order 4 (moments of degree 8) yields a globally minimal curvature of -0.094159 achieved at the 2 points x 1 = (-0.048892, -0.14076) and x 2 = (0.048896, -0.14076). With the two additional affine constraints p k (x) = g 1 (x k )(x -x k ) ≤ 0, k = 2, 3, the resulting set S has a globally minimal curvature of 1 certified at the LMI relaxation of order 4, and therefore it is a valid convex inner approximation of S, see Figure 3. This example illustrates that Algorithm 1 can work even when Assumption 1 is violated.

Here the singularity is removed by the additional affine constraints. This example also shows that symmetry of the problem can be exploited, since two global minimizers are found (distinct points with the same minimal curvature) to remove two nonconvex parts of the boundary simultaneously.

Singular quartic

Consider the quartic S = {x ∈ R 2 : p 1 (x) = x 4 1 + x 4 2 + x 3 2 ≤ 0} which has a singular point at the origin, hence violating Assumption 1.

Running Algorithm 1, we obtain the following sequence of bounds on the minimum curvature, for increasing LMI relaxation orders: GloptiPoly is not able to certify global optimality, so we can only speculate that the global minimum is zero and hence that set S is convex, see Figure 4. We may say that set S is numerically convex. Indeed if we strenghten the constraint p 1 (x) ≤ 0 into p 1 (x)+ǫ ≤ 0 for a small positive ǫ, say 10 -3 , then GloptiPoly 3 certifies global optimality and convexity with obj = -4.0627e-7 at the 4th LMI relaxation. On the other hand, if we relax the constraint into p 1 (x) + ǫ ≤ 0 with a negative ǫ = -10 -3 , then GloptiPoly 3 certifies global optimality and nonconvexity with obj = -0.22313 at the 4th LMI relaxation. We can conclude that the optimum of problem 3 is sensitive, or ill-conditioned, with respect to the problem data, the coefficients of p 1 (x). The reason behind this ill-conditioning is the singularity of S at the origin, see Figure 5 which represents the effect of perturbing the constraint p 1 (x) ≤ 0 around the singularity.

Control applications

In this section we focus on control applications of Algorithm 1, which is used to generate convex inner approximation of stability regions in the parameter space.

Third-order discrete-time stability region

Algorithm 1 can lend insight into the (nonconvex) geometry of the stability region. Consider the simplest non-trivial case of a third-order discrete-time polynomial x 1 + x 2 z + x 3 z 2 + z 3 which is stable (roots within the open unit disk) if and only if parameter x = (x 1 , x 2 , x 3 ) lies within the interior of compact region S = {x ∈ R 3 :

p 1 (x) = -x 1 -x 2 -x 3 -1 ≤ 0, p 2 (x) = x 1 -x 2 + x 3 -1 ≤ 0, p 3 (x) = x 2 1 -x 1 x 3 + x 2 -1 ≤ 0}.
Stability region S is nonconvex, delimited by two planes p 1 (x) = 0, p 2 (x) = 0 and a hyperbolic paraboloid p 3 (x) = 0 see e.g. [START_REF] Ackermann | Robust control: the parameter space approach[END_REF]Example 11.4].

Optimization problem (3) corresponding to convexity check of the hyperbolic paraboloid reads as follows:

min -2y 2 1 + 2y 1 y 3 s.t. x 2 1 -x 1 x 3 + x 2 -1 = 0 -x 1 -x 2 -x 3 -1 ≤ 0 x 1 -x 2 + x 3 -1 ≤ 0 (2x 1 -x 3 )y 1 + y 2 + x 3 y 3 = 0 y 2 1 + y 2 2 + y 2 3 = 1.
(

) 5
The objective function and the last constraint depend only on y, and necessary optimality conditions obtained by differentiating the Lagrangian -2y 2 1 + 2y 1 y 3 + t(y 2 1 + y 2 2 + y 2 3 -1) with respect to y yield the symmetric pencil equation

  -4 + 2t 0 2 0 2t 0 2 0 2t     y 1 y 2 y 3   = 0
From the determinant of the above 3-by-3 matrix, equal to t(t 2 -2t -1), we conclude that multiplier t can be equal to 1 -√ 2,0 or 1 + √ 2. The choice t = 0 implies y 1 = 0, y 2 = 1, y 3 = 0 which is inconsistent with the last but one constraint in [START_REF] Henrion | Ellipsoidal approximation of the stability domain of a polynomial[END_REF]. The choice

t = 1 - √ 2 yields y 1 = ±(1 + √ 2)α, y 2 = 0, y 3 = ±α with α = 1/ 4 -2 √
2 and the objective function -2y 2 1 + 2y

1 y 3 = -1 + √ 2. The choice t = 1 + √ 2 yields y 1 = ±α, y 2 = 0, y 3 = ±(-1 - √ 2)
α and the objective function -1 -√ 2, a negative minimum curvature. Therefore region S is indeed nonconvex.

From the remaining constraints in [START_REF] Henrion | Ellipsoidal approximation of the stability domain of a polynomial[END_REF], we conclude that the minimal curvature points x can be found along the portion of parabola √ 2x 2 1 -x 2 + 1 = 0 included in the half-planes (2+ √ 2)x 1 +x 2 +1 ≥ 0 and -(2+ √ 2)x 1 +x 2 +1 ≥ 0. Any plane tangent to the hyperbolic paraboloid p 3 (x) = 0 at a point along the parabola √ 2x 2 1 -x 2 + 1 = 0 can be used to generate a valid inner approximation of the stability region. For example, with the choice x 1 = (0, 1, 0), we generate the gradient half-plane p 4 (x) = g 3 (x 1 )(x -x 1 ) = x 2 -1 ≤ 0.

More generally, for discrete-time polynomials of degree n ≥ 3, stability region S is the image of the box B = [-1, 1] n (of so-called reflection coefficients) though a multiaffine mapping, see e.g. [START_REF] Nurges | Robust pole assignment via reflection coefficients of polynomials[END_REF] and references therein. The boundary of S consists of ruled surfaces, and the convex hull of S is generated by the images of the vertices of B through the multiaffine mapping. It would be interesting to investigate whether this particular geometry can be exploited to generate systematically a convex inner approximation of maximum volume of the stability region S.

Fixed-order controller design

Consider the open-loop discrete-time system (-2z 2 + 1)/(z 3 + z 2 + a), parametrized by a ∈ R, in negative closed-loop configuration with the controller (x 1 z + x 2 )/(z + 1). The characteristic polynomial is equal to q(z)

= 4 k=0 q k z k = z 4 + 2(1 -x 1 )z 3 + (1 -2x 2 )z 2 + (a + x 1 )z + a + x 2 ,
     p 1 p 2 p 3 p 4 p 5       =       -1 1 -1 1 -1 4 -2 0 2 -4 -6 0 2 0 -6 4 2 0 -2 -4 -1 -1 -1 -1 -1             q 0 q 1 q 2 q 3 q 4       .
The affine inequalities p k (x 1 , x 2 ) < 0, k = 1, 2, . . . , 5 define a polytope in the controller parameter plane (x 1 , x 2 ) ∈ R 2 , and the inequality p 6 (x 1 , x 2 ) < 0 defines a cubic region.

In the case a = 0, with the following Gloptipoly 3 implementation of Steps 1-3 of Algorithm 1: 2); g6 = diff(p6,x); % gradient H6 = diff(g6,x); % Hessian % LMI relaxation order order = input('LMI relaxation order = '); % build LMI relaxation P = msdp(min(y'*H6*y), p6==0, p1<=0, p2<=0, p3<=0, p4<=0, p5<=0, ... g6*y==0, y'*y==1, order); % solve LMI relaxation [status,obj] = msol(P) we obtain a negative lower bound obj = -3.5583 at the 2nd LMI relaxation, which is inconclusive. At the 3rd LMI relaxation, we obtain a positive lower bound obj = 0.8973 which certifies convexity of the stability region, see Figure 6. Since the stability region S = {x ∈ R 2 : p k (x) ≤ 0, k = 1, 2, . . . , 6} is convex, we can optimize over it with standard techniques of convex optimization. More specifically, a recent result in [START_REF] Lasserre | On convex optimization without convex representation[END_REF] indicates that any limit point of any sequence of admissible stationary points of the logarithmic barrier function f (x) = -6 k=1 log p k (x) is a Karush-Kuhn-Tucker point satisfying first order optimality condition. In particular, the gradient of f (x) vanishes at the analytic center of the set. Using Maple (or a numerical local optimization method) we can readily obtain the analytic center x * 1 ≈ 0.57975, x * 2 ≈ 0.13657 (five-digit approximations of algebraic coefficients of degree 17) corresponding to a controller well inside the stability region. Such a controller can be considered as non-fragile, in the sense that some uncertainty on its coefficients will not threaten closed-loop stability. Now for the choice a = -3/4 we carry on again our study of convexity of the stability region with the help of a similar GloptiPoly script. At the 2nd LMI relaxation we obtain a negative lower bound obj = -385.14 which is inconclusive. At the 3rd LMI relaxation, we obtain a negative lower bound obj = -380.88 which is also inconclusive. Eventually, at the 4th LMI relaxation, we obtain a negative lower bound obj = -380.87 which is certified to be the global minimum with status = 1. The point x 1 at which the minimum curvature is achieved is a vertex of the stability region, and the tangent at this point of the nonconvex part of the boundary is used to generate a valid inner approximation S, see Figure 7. Any point chosen in this triangular region corresponds to a stabilizing controller. We see that here the choice of the point of minimum curvature is not optimal in terms of maximizing the surface of S. A point chosen elsewhere along the negatively curved part of the boundary would be likely to generate a larger convex inner approximation.

mpol x y 2 p1 = -x(1)+x(2); p2 = -6*x(1)+4*x(2); p3 = -10*x(2)-4; p4 = -8+4*x(2)+6*x(1); p5 = -4+x(1)+x(2); p6 = 6*x(1)^2*x(2)+3*x(1)^2-10*x(1)*x(2)-2*x(1)-3*x(2)^3+6*x(2)^2+x(

Optimal control with semialgebraic constraints

In Model Predictive Control (MPC), an optimal control problem is solved recursively. This resolution is usually based on direct methods that consist of deriving a nonlinear program from the optimal control problem by discretization of the dynamics and the path constraints. Since the embedded software has strict specification on algorithm complexity and realtime computation, convexity of the program is a key feature [START_REF] Ross | Issues in the real-time computation of optimal control[END_REF]. Indeed, in this context, our convex inner approximation of the admissible space become valuable to speed up the computation even at the price of some conservatism.

In open-loop control design, convexity of the problem is a matter of concern especially when the optimal control problem is part of an MPC procedure. In this case, the optimal control problem is solved mostly using direct methods that transfom it into a parametric optimization problem. Convexity permits to limits the complexity of the resolution and so reduces the computation time of an optimal solution. Unfortunately, in dynamic inversion techniques based on differential flatness, the generally convex constraints on the states and inputs are replaced by nonconvex admissible sets in the flat output space, see [START_REF] Ross | Issues in the real-time computation of optimal control[END_REF] and reference therein for details. Thus, in such a method, it is necessary to design inner convex approximation of the admissible subset to develop a tractable algorithm [START_REF] Louembet | Path planning for satellite slew manoeuvers: a combined flatness and collocation based approach[END_REF].

Consider the following optimal control problem min x,u

t f t 0 u 2 (t)dt s.t. ẋ = 0 1 0 0 x + 0 1 u x(t 0 ) = x 0 , x(t f ) = x f p 1 (x) ≤ 0.
The objective of this problem is to steer the linear system from an initial state to a final state in a fixed time inside the admissible state subset S defined e.g. by the waterdrop quartic defined in section 4.3:

S = {x ∈ R 2 : p 1 (x) = x 4 1 + x 4 2 + x 2 1 + x 3 2 ≤ 0}. (6)
We describe thereafter a classical methodology for solving the previous optimal control problem using flatness-based dynamic inversion, see [START_REF] Milam | A new computational approach to realtime trajectory generation for constrained mechanical systems[END_REF][START_REF] Petit | Inversion based constrained trajectory optimization[END_REF][START_REF] Lévine | Flat output characterisation for linear systems using polynomial matrices[END_REF] for other examples. As the dynamics are linear and fully actuated, dynamic inversion can be used to develop an efficient algorithm for the considered problem [START_REF] Lévine | Flat output characterisation for linear systems using polynomial matrices[END_REF]. Thus, the system trajectory can be parametrized by a user-specified sufficiently smooth function x 1 (t) = f (t) so that x 2 (t) = ḟ (t) and u(t) = f (t). The function f (t) is classically described by a chosen basis b(t) and the associated vector of weighting coefficient α such that

f (t) = k α k b k (t).
In order to derive a finite dimensional program, the admissible set constraint is discretized and enforced at a finite number of time instants

{t i } i=1,...,N such that t 0 ≤ t 1 < t 2 < • • • < t N ≤ t f . Since S is nonconvex, we obtain a finite-dimensional nonlinear nonconvex program: min α i u 2 (α, t i ) s.t. x(α, t 0 ) = x 0 , x(α, t f ) = x f p 1 (x(α, t i )) ≤ 0, i = 1, . . . , N.
The inner approximation S calculated previously in section 4.3 is given by S = {x ∈ R 2 :

p 1 (x) ≤ 0, p 2 (x) = g 1 (x 1 )(x -x 1 ) ≤ 0, p 3 (x) = g 1 (x 2 )(x -x 2
) ≤ 0} where g 1 (x 1 ) and g 1 (x 2 ) is the gradient of p 1 (x) evaluated at x = x 1 and x = x 2 , respectively. The use of the inner approximation S as admissible subset leads to the following convex program:

min α i u 2 (α, t i ) s.t. x(α, t 0 ) = x 0 , x(α, t f ) = x f p 1 (x(α, t i )) ≤ 0, p 2 (x(α, t i )) ≤ 0, p 3 (x(α, t i )) ≤ 0, i = 1, . . . , N.
In the following, we set t 0 = 0, x 0 = [0.3000, -0.8000] and t f = 2.5, x f = [-0.3000, -0.8000].

The time function f (t) is a 5-segment-piecewise polynomial of the 4th order (degree 3) defined on a B-spline basis. We run both programs for different values of N. In table 1 we compare the computation times and optimal costs. See Figure 8 for the state trajectories.

For this example, we observe the positive effect that convexity has on the reduction of 1: Computation times and optimal costs of the nonconvex and convexified optimal control problems, as functions of the number N of discretization points. the computational burden, balanced by the relatively small loss of performance.

Conclusion

We have presented a general-purpose computational algorithm to generate a convex inner approximation of a given basic semialgebraic set. The inner approximation is not guaranteed to be of maximum volume, but the algorithm has the favorable features of leaving invariant a convex set, and preserving convex boundaries while removing nonconvex regions by enforcing linear constraints at points of minimum curvature. Even though our initial motivation was to construct convex inner approximations of stability regions for fixed-order controller design, our algorithm can be used on its own for checking convexity of semialgebraic sets.

Each step of the algorithm consists in solving a potentially nonconvex polynomial optimization problem with the help of a hierarchy of convex LMI relaxations. For this we use Gloptipoly 3, unfortunately with no guarantee of a priori computational burden, even though in practice it is observed that global optimality is ensured at a moderate cost, as soon as the dimension of the ambient space is small. Numerical experiments indicate that the approach may be practical for ambient dimensions up to 4 or 5. For larger problems, we can rely on more sophisticated nonlinear or global optimization codes [START_REF] Neumaier | Complete search in continuous global optimization and constraint satisfaction[END_REF], even though this possibility has not been investigated in this paper. Indeed, our main driving force is to contribute with a readily available Matlab implementation.

Our algorithm returns a sequence of polynomials such that the intersection of their sublevel sets is geometrically convex. However, the individual polynomials (of degree two or more) are not necessarily convex functions. One may therefore question the relevance of applying a relatively complex algorithm to obtain a convex inner approximation in the form of a list of defining polynomials which are not necessary individually convex.

A recent result of [START_REF] Lasserre | On convex optimization without convex representation[END_REF] indicates however that any local optimization method based on standard first-order optimality conditions for logarithmic barrier functions will generate a sequence of iterates converging to the global minimum of a convex function over convex sets. In other words, geometric convexity seems to be more important that convexity of the individual defining polynomials.

Indeed, if convexity of the inner approximation is guaranteed in the presented work, convexity of the defining polynomials would allow the use of constant multipliers to certificate optimality in a nonlinear optimization framework. Instead, with no guarantee of convexity of the defining polynomials, the geometric proprety of convexity of the sets is more delicate to exploit efficiently by optimization algorithms.

Finally, let us emphasized that it is conjectured that all convex semialgebraic sets are semidefinite representable in [START_REF] Helton | Sufficient and necessary conditions for semidefinite representability of convex hulls and sets[END_REF], see also [START_REF] Lasserre | Convex sets with semidefinite representation[END_REF]. It may then become possible to fully exploit the geometric convexity of our inner convex through an explicit representation as a projection of an affine section of the semidefinite cone. For example, in our target application domain, this would allow to use semidefinite programming to find a suboptimal stabilizing fixed-order controller.

Figure 3 :

 3 Figure 3: Nonconvex waterdrop quartic (light gray) and its convex inner approximation (dark gray) obtained by adding affine constraints at two points x 1 and x 2 of minimal curvature.

  5000 • 10 -1 -7.7502 • 10 -2 -8.5855 • 10 -3 -4.9525 • 10 -3

Figure 4 :

 4 Figure 4: Numerically convex singular quartic.

Figure 5 :

 5 Figure 5: Perturbed quartic p 1 (x) + ǫ ≤ 0 (bold line) can be convex (ǫ = 10 -3 ) or nonconvex (ǫ = -10 -3 ) near singularity of original quartic level set p 1 (x) = 0 (light line).

Figure 6 :

 6 Figure 6: Convex stability region (dark gray), with analytic center (cross) corresponding to a fixed-order controller.

Figure 7 :

 7 Figure 7: Convex inner approximation (dark gray) of nonconvex fourth-order discretetime stability region (light gray).

Figure 8 :

 8 Figure 8: Optimal trajectories (bold) in nonconvex admissible set (left) and in convex inner approximation (right).

  and it is Schur stable (all roots in the open unit disk) if and only if p k < 0, k = 1, 2, . . . , 4 and p 6 = -p 2 p 3 p 4 + p 2 2 p 5 + p 1 p 2 4 < 0 where 

Acknowledgements

The first author is grateful to J. W. Helton for many discussions and ideas leading to Theorem 1.


Jean-Baptiste Fournier

Comment on "Are stress-free membranes really 'tensionless' ?" by

Keywords: 87, 16, dj -dynamics and fluctuations 68, 03, Kn -Dynamics (capillary waves) 68, 35, Md -Surface thermodynamics, surface energies

In recent years, it has been argued that the tension parameter driving the fluctuations of fluid membranes should differ from the imposed lateral stress, the "frame tension" [1,2]. In particular, stress-free membranes were predicted to have a residual fluctuation tension. Some experimental evidence of this property has been reported in [3]. In a recent paper [4], Schmid argued that the reasoning published in [2] (where the frame tension is calculated by averaging the stress tensor and by differentiating the free energy), is inherently inconsistent-in the sense that a linearized theory, the Monge model, was used to predict a non-linear effect. In the present comment, we show that the criticism argued in [4] does not hold.

The correct way to proceed in the calculations is to be consistent at first-loop order, i.e., at first order in the small parameter k B T /κ, where k B T is the temperature in energy units and κ the membrane bending rigidity. In [4], however, all energies were given in units of k B T , and thus all tracks of the small parameter in which the expansion should be made were lost.

Restoring the k B T factors, eq. ( 3) of [4] reads

G(N, σ int , A p ) = σ int A p + k B T N -1 2 ln 2σ int λ 2 k B T -2 + 1 + σ int A p 4πκ(N -1) ln 1 + 4πκ(N -1) σ int A p . (1)
Here, G is the total free-energy of the membrane described at the Gaussian level in the Monge gauge, N the total number of fluctuating degrees of freedom (proportional to the total number of lipids), σ int the "internal" tension appearing in the Gaussian Hamiltonian

H = 1 2 d 2 r ⊥ [σ int (∇h) 2 + κ(∇ 2 h) 2
] functional of the membrane height function h(r ⊥ ), A p the fixed projected area above which the membrane stands, and λ a small distance cutoff in the direction orthogonal to the membrane.

Differentiating the free-energy with respect to the projected area, one obtains the frame tension σ f = ∂G/∂A p , yielding

σ f σ int = 1 + k B T 8πκ ln 1 + 4πκ(N -1) σ int A p - k B T 8πκ 4πκ(N -1) σ int A p + O(ǫ 2 ), (2)
Comment which corresponds to eq. ( 4) of [4]. Note that this result is only valid at first-order in the small parameter ǫ = k B T /κ because H is only the Gaussian approximation of the complete Helfrich Hamiltonian. The excess area (A -A p )/A p can be obtained from the relation A = ∂G/∂σ int , yielding

A -A p A p = k B T 8πκ ln 1 + 4πκ(N -1) σ int A p + O(ǫ 2 ), (3)
which corresponds to eq. ( 5) of [4]. Thus (A -A p )/A p is also a first-order quantity. The argumentation in [4] begins by defining

y = ln 1 + 4πκ(N -1) σ int A p ≃ 8πκ k B T A -A p A p = O(1). (4)
Schmid, however, fails to notice that this is a zeroth-order quantity (in ǫ): omitting k B T in the denominator, Schmid considers y as a first-order quantity, like (A -A p )/A p . However,

(A -A p )/A p = O(ǫ) is small, but 8πκ/(k B T ) = O(ǫ -1
) is large, and so the product is O(ǫ 0 = 1). Then, from eqs. ( 2) and ( 4), Schmid writes

σ f σ int = 1 - k B T 8πκ (e y -y -1) ≈ 1 - 4πκ k B T A -A p A p 2 + . . . . (5)
and says: "Hence the frame tension is found to differ from the internal tension to second order in (A -A p )/A p . On the other hand, the Monge model approximates planar interfaces only up to the first order in (∇h) 2 or (A -A p )/A p . Thus the results (6) and (4) 1 are not rigorous". However, since y = O(1), using the power series expansion e yy -1 = 1 2 y 2 + O(y 3 ) in eq. ( 5) is not justified. Even the result is in contradiction with the conclusion drawn: since 4πκ/(k B T ) = O(ǫ -1 ) and [(A -A p )/A p ] 2 = O(ǫ 2 ), the correction (σ f -σ int )/σ int appearing in eq. ( 5) would be a first-order correction.

The correct answer, at first-loop order (at first-order in k B T /κ), is thus the one given in eq. ( 2). In conclusion, the reasoning in [4] arguing that the results of Fournier and Barbeta are erroneous does not hold.


Anne-Florence Bitbol

Kévin Sin

Jean-Baptiste Fournier

Universal amplitudes of the Casimir-like interactions between four types of rods in fluid membranes

Keywords: 05.40.-a -Fluctuation phenomena, random processes, noise, and Brownian motion 87.16.dj -Membranes, bilayers, and vesicles: dynamics and fluctuations 87.16.Ka -Filaments, microtubules, their networks, and supramolecular assemblies

The fluctuation-induced, Casimir-like interaction between two parallel rods of length L adsorbed on a fluid membrane is calculated analytically at short separations d L. The rods are modeled as constraints imposed on the membrane curvature along a straight line. This allows to define four types of rods, according to whether the membrane can twist along the rod and/or curve across it. For stiff constraints, all the interaction potentials between the different types of rods are attractive and proportional to L/d. Two of the four types of rods are then equivalent, which yields six universal Casimir amplitudes. Repulsion can occur between different rods for soft constraints. Numerical results obtained for all ranges of d/L show that the attraction potential reaches kBT for d/L 0.2. At separations smaller than dc ≈ L(L/ p) 1/3 , where p is the rod persistence length, two rods with fixed ends will bend toward each other and finally come into contact because of the Casimir interaction.

When a field with long-range correlations is confined between objects that constrain its fluctuations, generic long-range forces arise. The most famous example is the Casimir interaction between two uncharged metal plates, which originates from the constraints they impose on the quantum electrodynamical vacuum fluctuations [1]. Similar interactions arise between objects immersed in a classical fluid, due to reduced thermal fluctuations. Fluids with long-range correlations produce long-range Casimirlike forces with universal characteristics. These interactions have been predicted to occur in binary mixtures close to a critical point, superfluids, liquid crystals and fluid membranes or interfaces [2,3]. Direct measurements of thermal Casimir-like forces have been achieved recently between a colloid and a surface in a critical binary mixture [4].

In biomimetic membranes, which are fluid interfaces characterized by their bending rigidity [5], Casimirlike forces have been predicted to arise between phaseseparated domains [6], integral proteins [6], adsorbed stiff rods [7] and semi-flexible polymers [8]. All these inclusions impose constraints on the thermal fluctuations of the membrane shape. The interactions between inclusions within biomembranes, and their organization, plays a central role in many cellular regulatory processes and is of fundamental importance [START_REF] Alberts | Molecular Biology of the Cell[END_REF]. Biologically relevant stiff polymers adsorbed on membranes include actin filaments [START_REF] Limozin | [END_REF]11], FtsZ rings [12] and the β-Amyloid fibrils involved in Alzheimer's disease [13,14]. Adsorption of goethite (α-FeOOH) nanorods on membranes has also been reported [15].

In this paper, we study the fluctuation-induced Casimirlike interaction free energy F C between two parallel rods of length L adsorbed on a membrane1 , at short and intermediate separations d < ∼ L. It can be calculated by studying rods whose only effect is to restrict the membrane fluctuations: then, F C corresponds to the d-dependent part of the free energy of the system. The Casimir-like force can be obtained by direct differentiation of F C . We focus on membranes dominated by their bending rigidity κ, i.e., we assume that the length scale (κ/σ) 1/2 , which compares its effect to that of the membrane tension σ, is much larger than all other relevant lengths. Other cases are briefly discussed at the end of this paper. In the asymptotic limit d L, previous analytical calculations showed that F C is attractive and decays as k B T (L/d) 4 [7], where k B is Botzmann's constant, and T is the temperature. Thus, even for d 3L, its magnitude ≈ 10 -4 k B T is far too small to effectively bring rods together. At short separations d L, one expects a constant scale-free Casimir interaction per unit length, i.e., F C ∝ -k B T L/d. Effective attraction at d < ∼ L will thus result if the numerical prefactor is of order unity. Golestanian [8] showed that the finite bending rigidity of semiflexible polymers, or polymer rods, embedded in membranes, has dramatic effects on their Casimir interaction at short separations: shape fluctuations, which locally reduce the separation of the polymer rods, become unstable. Precisely, the free energy eigenvalues of these modes were shown to be λ(q) (κ r /k B T )q 4 -(2A/π)/d 3 , where κ r is the bending rigidity of the polymer rods and A 2.91. Instabilities occur if λ(q min ) < 0 with q min ≈ π/L, i.e., for d < d c , where

d c L ≈ 2A π 5 1/3 L p 1/3 . ( 1
)
Here, p = κ r /(k B T ) is the persistence length of the rod.

For actin filaments [START_REF] Limozin | [END_REF]11] of length L 100 nm, with p 17 µm, this yields d c /L ≈ 5×10 -2 . For L/ p 10 -4 , which may be realistic for actin bundles, the undulation instability threshold goes down to d c /L ≈ 10 -2 . It is thus realistic to consider stiff rods at separations in the range 0.01 < ∼ d/L < ∼ 1. In this paper we shall give analytical results for d L, and numerical results for all ranges of d/L.

Let us describe the shape of the membrane by its height z = h(r) above a reference plane parametrized by r = (x, y). We shall refer to the latter coordinates as "inplane" and to z as "out-of-plane". In the Gaussian limit of small deformations around the flat shape, the effective Hamiltonian of the membrane reads [5]:

H = dr κ 2 ∇ 2 h 2 . ( 2
)
Note that the Gaussian curvature does not contribute, due to the Gauss-Bonnet theorem [START_REF] Chaikin | Principles of condensed matter physics[END_REF], and that the membrane tension σ has been neglected as previously discussed. The presence of a rigid rod lying within the membrane (or absorbed on it) is traditionally described by imposing the constraint [6,7]:

h(r) = a + b • r (3) for r in a narrow rectangle [x 0 -, x 0 + ] × [0, L],
where a and b describe translational and tilt degrees of freedom, respectively. This is fully equivalent to imposing in the same domain the constraints:

h yy = 0 , h xx = 0 , h xy = 0 , ( 4
)
where h ij ≡ ∂ i ∂ j h. An infinitely narrow rod can be described by imposing conditions (4) for x = x 0 and y ∈ [0, L]. The single condition h yy = 0, however, is sufficient to specify a rod. In the absence of the other conditions, the membrane can curve across the rod and twist along it. It follows that we may consider four types of rods (fig. 1):

h yy = h xx = h xy = 0.
The actual constraints imposed on the coarse-grained field h by a given rod depend on its microscopic properties, such as its width and its interactions with the lipids that constitute the membrane. Let us introduce the natural short length cutoff a of the theory, which is of order of the membrane thickness (it cannot be smaller since the membrane is considered as a thin surface): the field h is constructed through a coarse-graining, for instance a simple moving average, of the height of each lipid molecule over the lengthscale a. First, let us consider a rod with width smaller than a that interacts with the lipids through short-range (e.g. hydrophobic) interactions, which can be neglected at a distance a. Such a rod cannot affect h xx or h xy , and it will thus be of the "ct" type. On the contrary, a rod with a width equal to, or a little bit larger than a will be of the "c t" type. Indeed, this width is sufficient for the thin flat region arising from the presence of the rod to remain after the coarse-graining. A rod which has strong anisotropic attractive interactions with the lipids, with a range at least of order a, would behave as a "c t" rod too. In the case of intermediate ranges of interactions, the rod could constrain the slope but not the curvature of the membrane, in which case one would obtain the "c t" type. Finally, the "ct" type could correspond to a rod with width of order a and with short-range interactions, but composed of units free to twist relative to one another. In our coarse-grained description, all such rods can be considered as one-dimensional since their width is of order a or smaller. Hence, the microscopic properties of a rod (width, interactions) will appear only indirectly in our work, through the rod type. Let us first calculate the Casimir interaction between two parallel rods of the simplest "ct" type, in the d L regime. Let us assume periodic boundary conditions with period L along the y direction. Generalizing the constraint h yy = 0 to rods with a finite rigidity, we account for the presence of the rods by adding to H the following interaction energies between each rod and the membrane:

H i = L 0 dy κ ⊥ 2 h 2 yy (x i , y), (5)
for i ∈ {1, 2}, with x 1 = 0 and x 2 = d > 0. Thus, we allow out-of-plane fluctuations of the rods with a bending rigidity κ ⊥ = k B T ⊥ . The hard constraint h yy = 0 is recovered in the limit κ ⊥ → ∞. Note that we implicitly assume in any case an infinite in-plane rigidity κ → ∞, which is justified if d/L is not too small, as discussed above. Adapting the field-theoretical method of Li and Kardar [START_REF] Li | [END_REF] to these general constraints, we introduce two auxiliary fields ψ i (y), and we make the substitution

e -βH i = Dψ i e -1 2 dy -1 ⊥ ψ 2 i +i dr h(r)ψ i (y)δ(x-xi) (6)
in the partition function

Z = Dh exp[-β(H + H 1 + H 2 )]
of the membrane containing the two rods, where β = (k B T ) -1 . Performing the Gaussian integral over h(r) yields

Z = Dψ 1 Dψ 2 e -1 2L q Ψ t (-q)M (d,q)Ψ(q) , (7)
where

ψ i (q) = L 0 ψ i (y) exp(-iqy), Ψ(q) = (ψ 1 (q) ψ 2 (q)) t , M (d, q) = -1 ⊥ + Γ(0, q) Γ(-d, q) Γ(d, q) -1 ⊥ + Γ(0, q) (8) and Γ(x i , q) = k B T κ q 4 ∞ -∞ dp 2π e ipxi (p 2 + q 2 ) 2 = k B T 4κ |q| (1 + |q|x i ) e -|q|xi ( 9
)
is the Fourier transform along y of Γ(x i , y) = ∂ 4 G(x i , y)/∂y 4 , where G(r -r ) = r|(βκ∇ 4 ) -1 |r is the correlation function of the free membrane. Performing the Gaussian integrals over the fields ψ i (q) yields the free energy -k B T ln Z for two parallel "ct" rods. Discarding d-independent terms, we obtain the Casimir interaction

F ct|ct C (d) = 1 2 k B T q ln det [M (d, q)] = k B T L 2 dq 2π ln 1 - Γ(d, q) -1 ⊥ + Γ(0, q) 2 = k B T L d φ ct|ct (d) , ( 10
)
where

φ ct|ct (d) = ∞ 0 dx 2π ln 1 - 1 + x 1 + d/x 2 e -2x , (11) d = 4κ k B T d ⊥ . ( 12
)
It can be shown that this integral expression of the Casimir-like interaction is consistent with eqs. ( 4) and ( 8) of ref. [8]. Our method for calculating the Casimir-like interaction, however, is more direct, and we shall extract explicit forms for the asymptotic behaviors and extend it to other types of rods. Two asymptotic Casimir interaction regimes occur (fig. 2). First, in the "soft" limit where d tends to infinity (κ ⊥ → 0), the integral in eq. ( 11) can be performed analytically:

F ct|ct C,soft (d) - 7 8π d2 k B T L d = - 7k B T 128π (κ ⊥ /κ) 2 L d 3 . ( 13
)
This interaction, which scales as d -3 , is quite small, however, since it is proportional to d-2 .

Conversely, in the "hard" limit where d tends to zero (κ ⊥ → ∞), φ ct|ct reaches a constant limit, denoted by A ct|ct :

F ct|ct C (d) = A ct|ct k B T L d , ( 14
)
A ct|ct = ∞ 0 dx 2π ln 1 -(1 + x) 2 e -2x
≈ -0.46237. (15) This "hard" limit corresponds to rods imposing h yy = 0. Note that for ⊥ = βκ ⊥ → ∞, the q 4 factors in Γ(d, q) and Γ(0, q) cancel each other out in eq. [START_REF] Limozin | [END_REF]. Hence, imposing h yy = 0 on the rod is mathematically equivalent to setting h = 0 for x = x 0 and y ∈ [0, L] (by means of external forces and torques). This indicates that the relative tilt of the rods is effectively frozen in this regime. Physically, this comes from the fact that we are in the d L regime: a relative tilt of these long rods would cost a lot of energy. More formally, as we use periodic boundary conditions with period L in the y direction in our calculation, which amounts to neglecting the effect of the edges of the rods, h yy (x i , y) = 0 for y ∈ [0, L] entails that h(x i , y) is constant. Setting these constants to zero is legitimate since their presence simply corresponds to a possible change of reference plane.

For ordinary rods with κ r = κ ≈ κ ⊥ , which corresponds to p = ≈ ⊥ , it is the "hard" limit that is the physical one. Indeed,

d L gives d 4βκL/ ⊥ < 1, since the condition L < ⊥ k B T /(4κ) ≈ 10 -2
⊥ is required to avoid the in-plane undulation instability for small d/L, as previously discussed. Hence, the "hard" interaction, eq. ( 14), is the physical one for hard rods at separations d L. This asymptotic regime is attained slowly, however, so that one might rather expect to see an intermediate effective power law (see fig. 2). In principle, the "soft" 1 is violated, and in the left shaded area the in-plane undulation instability is expected to occur. The dashed lines show the asymptotic regimes ( 13) and ( 14). In the white region, one would measure βF ct|ct C ≈ -0.2 (L/d) 1.2 . As discussed in the text, the regime d 1 could be compatible with d/L 1 for anisotropic rods with a very soft out-of-plane bending rigidity.

limit can be reached by strongly anisotropic rods with a large κ and a small κ ⊥ , such that ⊥ k B T /(4κ) d L. Such rods should rather be called ribbons.

Let us now calculate the Casimir interactions between the other pairs of rods in the d L regime. Since the "hard" limit is more realistic, we shall consider the situation where the conditions h yy = 0, h xy = 0 and h xx = 0 are strictly imposed. Defining a vectorial differential operator D α for each rod type α ∈ {ct, c t, ct, ct }:

D ct = (∂ 2 y ) , D c t = (∂ x ∂ y , ∂ 2 y ) t , D ct = (∂ 2 x , ∂ 2 y ) t , D ct = (∂ 2 x , ∂ x ∂ y , ∂ 2 y ) t , (16)
the constraints imposed by a rod of type α can be written as D α h(x, y) = 0, for x ∈ {0, d} and y ∈ [0, L]. Following the steps detailed above in the case of two "ct" rods, we obtain:

F α|β C (d) = 1 2 k B T q ln det M α|β (d, q) , ( 17
)
where M α|β is the Hermitian block matrix:

M α|β (d, q) = G α|α (0, q) G α|β (-d, q) G β|α (d, q) G β|β (0, q) , ( 18
) in which (G α|β ) ij = D α i D β j G. Note that when D α i D β j con- tains ∂ 4
x (e.g., for two "ct" rods), the corresponding element of G α|α (x, q) includes a term proportional to δ(x). This term is discarded in x = 0, as in Ref. [START_REF] Li | [END_REF]. It is shown further on, on an example, that this regularization procedure is in agreement with zeta-function regularization. We thus obtain the Casimir interaction between a rod of type α and a rod of type β:

F α|β C (d) = A α|β k B T L d , ( 19
)
with A α|β (= A β|α ) given by

A α|β = ∞ 0 dx 2π ln 1 -f α|β (x) . ( 20
)
The function f α|β depends on the types of rods at stake:

f ct|ct (x) = (1 + x) 2 e -2x , ( 21
)
f ct|c t(x) = 1 + 2 x + 2 x 2 e -2x , ( 22
)
f ct|ct (x) = (1 + 2 x) e -2x , ( 23
) f c t|c t(x) = 2 1 + 2 x 2 e -2x -e -4x , ( 24
)
f c t|ct (x) = 4 x e -2x + e -4x , ( 25
)
f ct|ct (x) = 2 e -2x -e -4x . ( 26
)
We find that the results involving "c t" rods are strictly identical to those involving "ct" rods:

f ct |β ≡ f ct|β for all β.
The numerical values of the A α|β coefficients are given in table 1. They are all negative, which means that all the rods attract each other. Casimir interactions are ordinarily attractive for like boundary conditions and repulsive for unlike boundary conditions [3]. Here, however, the situation is not so simple, since all the rods share the h yy = 0 condition, while they differ on the other ones. We obtain increasing interactions upon going from "ct" to "ct/c t" then to "c t". Besides, as mentioned previously, "ct" and "c t" rods behave equivalently. However, this equivalence holds only in the "hard" limit. Indeed, if we put a finite rigidity κ ⊥ on h yy , as in eq. ( 5), the "ct" rod vanishes altogether in the limit κ ⊥ → 0 while the "c t" rod becomes a different object: a flexible "ribbon" setting only h xx = 0. We have calculated the Casimir interaction between this type of ribbon and a "ct" rod imposing h yy = 0. We find A rib|ct 0.012 > 0, i.e., these antagonistic boundary conditions (h yy = 0 vs. h xx = 0) indeed yield a repulsion.

Exact calculations of the fluctuation-induced interactions between rods at short separations d L were made possible by neglecting boundary effects, which results into translational invariance. Formally, this was implemented through periodic boundary conditions. At larger separations, the effect of the edges of the rods cannot be neglected anymore. However, the Casimir interaction energy can be computed numerically for all d/L. Numerical schemes based on the Li and Kardar method have been developed in Refs. [18,19]. The idea is to discretize the functional integral over the auxiliary fields ψ i that enforce the constraints. Using a discretization scheme that will be discussed elsewhere [20], we have calculated the Casmir interaction between two rigid "ct" rods at any separation (fig. 3). Although the microscopic cutoff wavelength is explicitly involved in our calculations, we obtain universal results that only depend on d/L. At large separations, we recover the (L/d) 4 power-law obtained by Golestanian et al. [7]. At small separations, we recover our result (14).

As can be seen on fig. 3, the Casimir interaction potential reaches -k B T for d/L 0.2 (while using the asymptotic form would yield the overestimated value d/L 2). Thus, two rods approaching each other at separations less than half their length will effectively attract each other due to the Casimir effect.

In order to check the regularization procedure mentioned after Eq. ( 18), we have recalculated the interaction between two rigid rods of the "ct" type by using zeta-function regularization [START_REF] Elizalde | Ten physical applications of spectral zeta functions[END_REF], which is rather straightforward for such rods.

To obtain the Casimir energy F ct|ct C (d), we calculate the partition function Z = Dh exp[-βH] of the membrane patch P = [0, d] × [0, L] situated between the rods (the free energy of the rest of the membrane is independent of d). The functional integral in Z runs over the functions h(r) that comply with the boundary conditions imposed by the rods, i.e., D ct h(x, y) = 0 for x ∈ {0, d}, and y ∈ [0, L]. Let us as-sume, in addition, that h(x, y) = 0 on the rods (it will not change the result since we have seen previously that the relative tilt of the rods is effectively frozen). Finally, let us take periodic boundary conditions with period L in the y direction. Then, H = 1 2 κ P dr hOh, where O = ∇ 4 is the biharmonic operator on P with the above-mentioned boundary conditions. Since O is positive definite Hermitian, we can write Z ∝ 1/ det(O). Thus, discarding dindependent terms, we obtain

F ct|ct C (d) = 1 2 k B T Tr ln(O). The eigenvalues of O are λ nm = π 4 n 2 d 2 + 4 m 2 L 2 2 , n ∈ N * , m ∈ Z, (27)
and the associated eigenfunctions are sin(nπx/d) × exp(im2πy/L). Hence, we obtain

F ct|ct C (d) = 1 2 k B T ∞ n=1 ∞ -∞ dm ln(λ nm ) . ( 28
)
While Eq. ( 28) is divergent, its finite Casimir part can be extracted by calculating

ζ O (s) = ∞ n=1
∞ -∞ dm λ -s nm , which is finite for s > 1 4 , and yields by analytical continuation:

F ct|ct C (d) = - 1 2 k B T lim s→0 ζ O (s) = -k B T π 12 L d . ( 29
)
This result is in agreement with the value A ct|ct = -π/12 obtained above with the regularization of Li and Kardar [START_REF] Li | [END_REF]. This agreement also implies that a half-plane imposing the same boundary conditions as a rod will yield the same Casimir interaction (in the regime L d). Physically, this equivalence is due to the fact that the Casimir interaction arises from the confinement of the fluctuation modes between the rods (or the half-planes), i.e., in the patch P.

Finally, we have checked that the undulation instability discussed in Ref. [8] can be recovered from the Casimir interaction calculated in this paper. Let us consider the fundamental deformation mode:

δx = u(y) = U sin(πy/L ) , y ∈ [0, L ] , ( 30
)
where δx is the in-plane displacement of each rod in the direction of the other rod. Inextensibility implies L = L -π 2 U 2 /(4L) + O(U 3 ). In the limit d L, we may use the proximity-force (or Derjaguin) approximation [START_REF] Derjaguin | [END_REF] to calculate the total free energy of the two rods:

F def k B T L 0 dy A α|β d -2u + p u 2 - A α|β L d + O(U 3 ) = A α|β L d 4U πd 1 + πU 2d + π 4 p 2L U 2 L 2 + O(U 3 ). ( 31
)
Minimizing F def with respect to U , we find that the equilibrium deformation amplitude U eq reaches d/2, i.e., that the two rods touch each other when d < d c , where

d c L = -4(2 + π)A α|β π 5 1/3 L p 1/3 . ( 32
)
We thus obtain an instability threshold very similar to that of eq. ( 1). For instance, for two "ct" or "c t" rods, with A α|β -0.46, the prefactor of (L/ p ) 1/3 is 0.31, while in eq. ( 1) it is 0.27. In this paper, we have assumed that the effect of the membrane tension σ was negligible with respect to that of its bending rigidity κ. This hypothesis is valid in usual flaccid membranes (σ ≈ 10 -7 N/m) for rods shorter than (κ/σ) 1/2 ≈ 1 µm. However, σ can be increased up to about 10 -3 N/m by applying an external tension to the membrane: in this case, the effect of tension will become important. The methods presented here in the bending rigidity-dominated regime can also be applied to the tension-dominated limit, where L and d are much larger than (κ/σ) 1/2 . For instance, the interaction between two "ct" rods imposing h yy = 0 is then still given by Eq. 19, but with A ct|ct = -π/24. This result can be obtained both through the method based on Ref. [START_REF] Li | [END_REF] and through zeta-function regularization. This tension-dominated interaction is relevant in liquid-vapor interfaces [23]. A general calculation of the Casimir-like interaction between two rods, including both κ and σ, is beyond the scope of this paper. This problem has been addressed very recently in the simpler case of two small circular inclusions [24].

In summary, we have studied the Casimir interaction between two parallel rods adsorbed on a fluid membrane characterized by its bending energy and subject to thermal fluctuations. We have considered four types of boundary conditions, depending on whether the membrane can twist along the rod and/or curve across it. Whatever their type, rigid rods of length L at short separations d L, experience an attractive Casimir interaction scaling as k B T L/d. In this regime, two of the four types of rods are equivalent, which yields six universal Casimir amplitudes. When taking into account the finite persistence length p of the rod, we find, in agreement with Ref. [8], that at separations less than d c ∼ L(L/ p ) 1/3 two rods with fixed ends will bend toward each other and finally come into contact because of the Casimir interaction. Besides, we have shown that rods with a very soft out-of-plane bending rigidity, but a large in-plane bending rigidity, still experience a powerlaw Casimir interaction (albeit small), which scales as d -3 . For a large, but not infinite out-of-plane bending rigidity, an effective power law in d -α with α close to 1 would be observed. This nuances the conclusion of Ref. [8] according to which a finite out-of-plane rigidity, no matter how large, should destroy the long-range Casimir interaction. We have also calculated numerically the Casimir interaction between two rods at any separations, finding at d/L ≈ 0.5 a crossover between the short separation d -1 behavior and the large separation d -4 behavior. We have shown, in the case of the simplest boundary condi-tions, that at separations d/L < ∼ 0.2 the Casimir attraction reaches k B T . In conclusion, due to the very long range of the Casimir interaction, long rods adsorbed on an effectively tensionless membrane will actually attract each other-by the Casimir effect-when their separation becomes significantly smaller than their length.

Fig. 1 :

 1 Fig. 1: Four types of rods for which the membrane may or may not curve (c) across the rod or twist (t) along it.

Fig. 2 :

 2 Fig. 2: Logarithmic plot of the Casimir interaction (10) between two "ct" rods of length L and separation d, calculated analytically under the assumption d/L 1, for L = 10 -3 p ≡ 10 -3 ⊥ and κ = 25 kBT (yielding d = 0.1d/L). In the right shaded area the condition d/L1 is violated, and in the left shaded area the in-plane undulation instability is expected to occur. The dashed lines show the asymptotic regimes (13) and (14). In the white region, one would measure βF

Fig. 3 :

 3 Fig. 3: Logarithmic plot of the universal Casimir interaction between two rigid "ct" rods of length L and separation d, calculated numerically in all regimes of d/L. Circles (resp. squares) corresponds to a length L of the rods 75 (resp. 50) times the microscopic cutoff. The straight line shows the analytical asymptotic behavior F ct|ct C /(kBT ) = A ct|ct d/L, valid for d/L 1.

Table 1 :

 1 Universal amplitudes A α|β of the Casimir interactions between the various types of rods, rounded off to the third decimal place. The exact value of A ct|ct is -π/12.

	H α H H H H β ct	ct -0.262 -0.357 -0.549 ct / ct c t
	ct / ct	-0.462 -0.672
	c t	-0.924

Our model is valid both for embedded and adsorbed rods.