Italo Atzeni
email: italo.atzeni@huawei.com

Marco Maso
email: marco.maso@huawei.com

Imène Ghamnia
email: imene.ghamnia@centralesupelec.fr

Ejder Baştug

Mérouane Debbah
email: merouane.debbah@huawei.com

Flexible Cache-Aided Networks with Backhauling

Caching at the edge is a promising technique to cope with the increasing data demand in wireless networks. This paper analyzes the performance of cellular networks consisting of a tier macro-cell wireless backhaul nodes overlaid with a tier of cache-aided small cells. We consider both static and dynamic association policies for content delivery to the user terminals and analyze their performance. In particular, we derive closedform expressions for the area spectral efficiency and the energy efficiency, which are used to optimize relevant design parameters such as the density of cache-aided small cells and the storage size. By means of this approach, we are able to draw useful design insights for the deployment of highly performing cache-aided tiered networks.

I. INTRODUCTION

A tremendous increase in smartphone usage during the last decade resulted in a 4000-fold mobile data traffic explosion in cellular networks [START_REF]Cisco Visual Networking Index: Global Mobile Data Traffic Forecast[END_REF]. In this setting, mobile operators constantly look for innovative solutions to satisfy performance indicators such as higher spectral and energy efficiency, improved coverage, and lower end-to-end delays: this calls for the development of innovative techniques at both device and network level in the coming years [START_REF] Andrews | What will 5G be?[END_REF]. Caching at the edge is regarded as one of the most promising approaches to cope with the increasing data demand [START_REF] Baştug | Cache-enabled small cell networks: modeling and tradeoffs[END_REF]. This is particularly true for heterogeneous and tiered network layouts, for which the extent of the performance gains brought by caching has been extensively assessed in terms of several performance metrics of practical relevance. For instance, [START_REF] Li | Optimization and analysis of probabilistic caching in n-tier heterogeneous networks[END_REF] and [START_REF] Liu | Caching policy toward maximal success probability and area spectral efficiency of cache-enabled HetNets[END_REF] focus on optimal probabilistic caching policies for enhanced content delivery, [START_REF] Abd-Elmagid | Cache-aided heterogeneous networks: Coverage and delay analysis[END_REF] studies the impact of caching on the coverage and on the delay experienced by the user terminals (UTs), and [START_REF] Khreishah | Joint caching, routing, and channel assignment for collaborative small-cell cellular networks[END_REF] jointly considers caching, routing, and channel assignment strategies for collaborative small cells (SCs).

Along similar lines, this work analyzes the system performance of cache-aided tiered networks comprised of a tier of macro-cell backhaul (BH) nodes overlaid with a tier of non-cooperative SCs [START_REF] Maso | Cacheaided full-duplex small cells[END_REF]. Each SC in the considered network has (limited) storage capabilities and pre-fetches popular files before the downlink transmission to its associated UTs. In this context, and differently from what has been previously proposed in the literature, our goal is to optimize relevant performance metrics for network planners and operators, and draw design guidelines for the deployment of highly performing cache-aided tiered networks. Using consolidated random This research has been supported in part by the ERC Starting Grant 305123 MORE (Advanced Mathematical Tools for Complex Network Engineering), the U.S. National Science Foundation under Grant CCF-1409228. spatial models to characterize the location of SCs, BHs and UTs [START_REF] Atzeni | Full-duplex MIMO small-cell networks with interference cancellation[END_REF], we investigate the area spectral efficiency (ASE) gains brought by caching popular files in the considered network using tools from stochastic geometry. To increase the generality and relevance of our study, we examine both static and dynamic association policies for content delivery to the UTs. In particular, in the static case, UTs can only be associated with SCs, regardless of the occurrence of cache hit/miss events at the latter; conversely, in the dynamic case, UTs are associated with SCs in case of cache hit and with wireless BHs otherwise.

The first contribution of this paper is the derivation of closed-form expressions and bounds for the ASE and the energy efficiency (EE) with the aforementioned UT association policies. Such expressions are subsequently used to optimize crucial design parameters such as SC density and storage size. Our findings highlight that the impact of caching on the ASE is minor as compared to that of the SC density, i.e., densifying the network is always more convenient than expanding the storage size at the SCs. Conversely, we show that increasing the storage size rather than the SC density is more effective in terms of EE of the network. In practice, the choice of investing on additional SCs or on more storage strongly depends on the goal of the network operator. Finally, it is worth observing that the performance of the considered network heavily hinges on the adopted UT association policy: a dynamic policy is always preferable whenever the network can support the larger complexity necessary for its effective implementation.

II. SYSTEM MODEL A. Network Model

Let us consider a set of mobile UTs in a two-tier ultradense network, where the lower tier comprises a set of noncooperative SCs ensuring network coverage and the upper tier consists in a set of BHs connected to the internet. In our model, each SC is associated with only one UT and one BH. We assume that the spatial distribution of the network nodes follows the stationary, independently marked Poisson point processes

(PPP) Φ x, u(x), b(x) ⊂ R 2 × R 2 × R 2 .
We use Φ SC {x} to denote the PPP of the SCs with spatial density λ (measured in SCs/m 2 ), with isotropic marks Φ UT u(Φ SC ) = {u(x)} x∈ΦSC and Φ BH b(Φ SC ) = {b(x)} x∈ΦSC representing the UTs and the BHs, respectively; evidently, Φ UT and Φ BH are dependent on Φ SC and have the same spatial density λ. In this setting, the employment of PPPs allows to capture the randomness of practical ultra-dense network deployments and, at the same time, obtain precise and tractable expressions for the system-level performance metrics [START_REF] Atzeni | Full-duplex MIMO small-cell networks with interference cancellation[END_REF]- [START_REF] Atzeni | Downlink cellular network analysis with LOS/NLOS propagation and elevated base stations[END_REF]. Let r yz y-z denote the distance between nodes y, z ∈ Φ: the distances of the UTs and the BHs from their associated SCs are assumed fixed and are denoted by R UT r xu(x) and R BH r b(x)x , ∀x ∈ Φ SC , respectively, with R BH > R UT .

B. Caching Model and UT Association

Let F {f i } F i=1 be a F -sized subset of all the files available in the internet, entirely accessible by each BH. Without loss of generality, and for simplicity in the notation, we assume that all F files have identical lengths. In the following, we refer to F as file catalog. Each UT requests files from F with probability denoted by P {p 1 , p 2 , . . . , p F }, with F i=1 p i = 1, which hinges on the files popularity over the whole network.

To offload the overlaying BH infrastructure, each noncooperative SC x is equipped with a storage unit ∆ x of size S ≤ F (measured in files/SC). Suppose that, at a given time instant, UT u(x) is interested in downloading file f i ∈ F: if f i is cached at its associated SC x, i.e., f i ∈ ∆ x , we have a cache hit event; conversely, if f i / ∈ ∆ x , we have a cache miss event. In this regard, we use / ∆ x to denote a cache miss event at x. Accordingly, we introduce the indicator function

1 / ∆ x 1, if / ∆ x 0, otherwise (1)
where its logical complement is denoted by 1 ∆x .

The effectiveness of the adopted cache-aided approach depends on the probability that any file requested by a given UT is cached at its serving SC, which is referred to as cache hit probability and is denoted by P hit . This quantity typically depends on the adopted caching policy, content placement, and distribution of the UT requests. Further details on these aspects are unnecessary for the scope of our work and we refer to [START_REF] Baştug | Cache-enabled small cell networks: modeling and tradeoffs[END_REF] and references therein for further details. Now, let R = T B ∈ R be the time-frequency resource used by the network for delivering one file from the BH tier to the UTs, with T ∈ R (resp. B ∈ R) defined as the amount of necessary time (resp. frequency) resources to perform this operation. Two UT association policies are considered in this work, i.e., static and dynamic. When the static UT association is adopted, operations in the two tiers occur on orthogonal resources. In this case, the UTs can only be served by their associated SCs, and the BH-to-SC/SC-to-UT links occupy two R 2 -sized resources: for instance, each link may use the entirety of B for half of the time (i.e., T 2 ) or the entirety of T for half of the bandwidth (i.e., B 2 ). Conversely, under dynamic UT association, operations in the two tiers occur on the same resource. In this case, the UTs are served by either the SCs or the BHs depending on the availability of the requested file at the SCs, and the SC-to-UT/BH-to-UT links occupy a R-sized resource. In practice, the second approach leads to a more efficient use of the time-frequency resource at the cost of more sophisticated inter-tier coordination and possibly

Static Dynamic

Cache hit Cache miss Cache hit Cache miss signaling. For notational simplicity, in the rest of the paper we differentiate between static and dynamic UT association by means of the super-indices S and D, respectively. The adopted UT association policy has a strong impact on the network response to a cache hit/miss event when UT u(x) requests file f i . In particular:

R 2 R 2 R 2 R R } BH } SC } UT
1) Static: If a cache hit occurs, then SC x transmits f i to u(x). Conversely, if a cache miss occurs, first x retrieves f i from its associated BH b(x), then transmits it to UT u(x). In both cases, an R 2 -sized resource is used by the SC-to-UT link due to the static UT association; 2) Dynamic: If a cache hit occurs, then x transmits f i to u(x). Conversely, if a cache miss occurs, u(x) can bypass x and connect directly to b(x) to download f i . In both cases, the whole resource R is used by either the SC-to-UT or the BH-to-UT link. A graphical representation of these operations, and their impact on the network resource usage, is provided in Figure 1.

C. Channel Model

In this work, all nodes have a single transmit/receive antenna. In addition, we assume that the SCs and the BHs transmit with powers ρ SC and ρ BH (measured in W), respectively. The propagation through the wireless channel combines pathloss attenuation and small-scale fading. For the former, we adopt the standard power-law pathloss model and define the pathloss function (y, z) r -α yz between nodes y, z ∈ Φ, with pathloss exponent α > 2. For the latter, we use h yz to denote the channel power fading gain between nodes y, z ∈ Φ: all channels are assumed to be subject to Rayleigh fading and, hence, h yz ∼ exp(1). 1 The signal-to-interference-plus-noise ratio (SINR) at SC x reads as

SINR x ρ BH R -α BH h b(x)x I x + σ 2 (2)
where I x is the overall interference power at SC x and σ 2 is the additive noise power. Likewise, the SINR at UT u(x) reads as

SINR u(x) ρ SC R -α UT h xu(x) I u(x) + σ 2 (3)
where I u(x) is the overall interference power at u(x). For static UT association, we have

I x = I (S) x and I u(x) = I (S)
u(x) , with

I (S)
x y∈ΦSC\{x}

ρ BH r -α b(y)x h b(y)x 1 / ∆ y (4)
I (S) u(x) y∈ΦSC\{x} ρ SC r -α yu(x) h yu(x) (5)
respectively. On the other hand, in case of dynamic UT association, we are only interested in the interference at the UT and we thus have

I u(x) = I (D)
u(x) , with

I (D) u(x) y∈ΦSC\{x} ρ SC r -α yu(x) h yu(x) 1 ∆y + ρ BH r -α b(y)u(x) h b(y)u(x) 1 / ∆ y . (6)

III. PERFORMANCE ANALYSIS

The system performance is closely related to the probability of the UTs successfully receiving the requested content, which is termed as success probability. This depends on two factors, i.e., the cache hit probability (cf. Section II-B) and the SINRs of the file transmissions (cf. Section II-C). Thus, in the following, we derive the success probability using both the considered UT association policies.

Our analysis focuses on a randomly chosen SC x, and its marks u(x) and b(x), referred to as typical SC, typical UT and typical BH, respectively: due to Slivnyak's theorem and to the stationarity of Φ, these nodes are representative of the whole network [START_REF] Atzeni | Full-duplex MIMO small-cell networks with interference cancellation[END_REF]. For a given SINR threshold θ, we assume that a file is successfully received by the typical UT if SINR x > θ∧SINR u(x) > θ and SINR u(x) > θ for static and dynamic UT association, respectively. The success probability is formalized in the following theorem, where we use the notation

Υ(z) πz 2 α csc 2π α α . (7)
Theorem 1. For static UT association, the success probability is given by

P (S) suc (θ) exp -θ σ 2 ρSC R α UT L (S) I u(x) θρ -1 SC R α UT × P hit + (1 -P hit ) exp -θ σ 2 ρBH R α BH L (S) Ix θρ -1 BH R α BH ( 8
)
where

L (S) Ix (s) exp -2πλ(1 -P hit )Υ(ρ BH s) (9)
L (S)

I u(x) (s) exp -2πλΥ(ρ SC s) (10)
are the Laplace transforms of the interference terms I (S) x and I (S) u(x) in ( 4)-( 5), respectively. For dynamic UT association, the success probability is given by

P (D) suc (θ) P hit exp -θ σ 2 ρSC R α UT L (D) I u(x) θρ -1 SC R α UT + (1 -P hit ) × exp -θ σ 2 ρBH R α BH 2π 0 L (D) I u(x) θρ -1 BH Ω(R UT , R BH , φ) dφ 2π (11)
with

Ω(R UT , R BH , φ) R 2 UT +R 2 BH +2R UT R BH cos φ α 2
and where

L (D) I u(x) (s) exp -2πλP hit Υ(ρ SC s) × exp -2πλ(1 -P hit )Υ(ρ BH s) (12)
is the Laplace transform of the interference term I (D) u(x) in ( 6). Proof: See Appendix I.

Since a closed-form expression for the success probability P (D) suc (θ) in ( 11) is not available, we provide a useful lower bound in the following corollary.

Corollary 1. For dynamic UT association, the success probability can be lower-bounded as

P (D) suc (θ) ≥ P hit exp -θ σ 2 ρSC R α UT L (D) I u(x) θρ -1 SC R α UT + (1 -P hit ) exp -θ σ 2 ρBH R α BH L (D) Ix θρ -1 BH R BH + R UT α . ( 13
)
Proof: The lower bound is obtained by considering the worst-case distance between BH and UT, i.e., R BH + R UT (cf. [START_REF] Atzeni | Full-duplex MIMO small-cell networks with interference cancellation[END_REF]).

The expressions above can be used to obtain other useful performance metrics. In this paper, we focus on the achievable ASE (measured in bps/Hz/m 2 ), which is readily obtained as

ASE (S) (θ, λ) 1 2 λP (S) suc (θ) log 2 (1 + θ) (14)
ASE (D) (θ, λ) λP (D) suc (θ) log 2 (1 + θ) (15)
for static and dynamic UT association, respectively. At this stage, it is worth observing that, in practice, the distances R UT and R BH depend on the SC density λ. Thus, in the following we set

R UT = βUT 2 √ λ and R BH = βBH 2 √
λ , with β BH > β UT ; observe that 1 2 √ λ represents the average distance between nodes in a PPP with spatial density λ. Figure 2 illustrates the achievable ASEs as defined in ( 14)-(15) against the cache hit probability P hit , with λ = 10 -2 SCs/m 2 , α = 4, θ = 1, σ 2 = 0 (i.e., interference limited case), ρ SC = 0.5 W, ρ BH = 1 W, β SC = 0.5, and β BH = 1 (the same parameters will be used in Section IV-A). The beneficial impact of caching is evident: as compared with the cache-free setup (which corresponds to P hit = 0), the achievable ASE with P hit = 0.25 (resp. P hit = 0.5) is increased by 188% (resp. 264%) for dynamic UT association and by 158% (resp. 218%) for static UT association. Furthermore, we observe that the dynamic UT association always outperforms the static UT association. Lastly, the lower bound on (15), derived using Corollary 1, is increasingly tight as P hit grows large.

IV. CACHE-AIDED SC NETWORK DEPLOYMENT

In this section, we look at possible design problems of interest to network operators and we aim at providing simple design guidelines for the deployment of cache-aided SC networks. As design parameters, we focus on the SC density λ ∈ [λ (min) , λ (max) ] and the storage size S ∈ [0, S (max) ], with S (max) ≤ F : therefore, we make the dependence on these parameters explicit in our performance metrics. For instance, the minimum SC density λ (min) can be chosen to ensure an interference-limited system performance (i.e., the background noise becomes negligible), whereas the maximum SC density λ (max) can be due to spatial limitations; similarly, the maximum storage size S (max) may arise from practical constraints regarding space or power consumption at the SCs.

In our analysis, we consider the simplest case of uniform files popularity, i.e., p i = 1 F , i = 1, . . . , F , with P hit (S) = S F , where we have made the dependence of the cache hit probability on S explicit. Note that, for the more realistic case of decreasing files popularity, i.e., p 1 > p 2 > . . . > p F , expanding the storage will have an increasingly less significant impact on the system performance.

A. Deployment Cost

Suppose that a network operator has a fixed monetary budget per m 2 c (measured in $/m 2 ) to deploy a cache-aided SC network. Let p λ denote the SC price (measured in $/SC) and let p S denote the storage price (measured in $/file). The following simple question naturally arises in this context: is it better to invest in deploying SCs or storage? The answer hinges on the choice of the performance metric and can be found by means of an appropriate optimization. For instance, to maximize the achievable ASE in ( 14)-( 15), the answer is given by the solution of the following optimization problem: (max) ].

max λ,S ASE(λ, S) s.t. p λ λ + p S λS ≤ c λ ∈ [λ (min) , λ (max) ] S ∈ [0, S
(P1) Alternatively, the system performance can be measured in terms of its EE. Before formalizing the corresponding problem, some definitions are in order. Let E tot (S) P hit (S)E hit + 1 -P hit (S) E miss (16) be the total energy consumed to provide a content to a UT (measured in J), where E hit is the energy needed by the SC to retrieve a cached file from its storage unit and transmit it to the UT in case of cache hit, and E miss is the energy needed by the BH to retrieve a non-cached file from the internet and transmit it to the UT in case of cache miss (via the SC for static UT association and directly for dynamic UT association). Moreover, we define the area energy consumption (measured in J/m 2 ) as AEC(λ, S) λE tot (S) and the EE (measured in bit/J) as

EE(λ, S) ASE(λ, S) AEC(λ, S) . (17)
Hence, the corresponding optimization problem is given by: max λ,S EE(λ, S) (max) ].

s.t. p λ λ + p S λS ≤ c λ ∈ [λ (min) , λ (max) ] S ∈ [0, S

(P2)

Using the same parameters as in Section III, we additionally set F = 10 7 files, λ (min) = 10 -4 SCs/m 2 , λ (max) = 10 -2 SCs/m 2 , S (max) = 5 × 10 6 files/SC, p λ = 250 $/SC, p λ = 0.005 $/file, E hit = 1 J, and E miss = 10 J. Assuming p λ λ + p S λ S = c, which excludes the trivial solution {λ , S } = {λ (max) , S (max) }, Figure 3 illustrates the feasible solution set of the considered problems with c = {1, 2.5, 5} $/m 2 .

Figures 4 and5 plot the achievable ASE and the EE, respectively, against the SC density λ, where each value of λ yields a value of S (and thus of P hit ) according to the constraint on the monetary budget (cf. Figure 3). The achievable ASE in Figure 4 is monotonically increasing in λ: this means that, to maximize the network ASE, one should invest in deploying more SCs (even if that means S = 0), whereas S should be increased only if there is remaining budget. Therefore, the solution of problem (P1) is given by

λ = min c p λ , λ (max) , S = 1 p S c λ -p λ . ( 18
)
On the contrary, the EE in Figure 5 is monotonically decreasing in λ: this means that, to maximize the network EE, one should invest in expanding the storage as much as possible (even if that means λ = λ (min) ), whereas λ should be increased only if there is remaining budget. Hence, the solution of problem (P2) reads as

S = min 1 p S c λ (min) -p λ , S (max) , λ = c p λ + p S S .
(19) In practice, increasing the SCs caching capabilities seems to be very useful for increasing the EE. Conversely, the impact of additional storage on the ASE is minor if the operator disposes of additional resources to increase the SC density. These results suggest that the answer to the question initially posed in this section, i.e., the choice of investing on additional SCs or on more storage, is not unique but strongly hinges on the ultimate goal of the network operator. Finally, we observe that a dynamic UT association policy is always preferable in terms of network performance whenever the larger complexity necessary to implement it can be afforded.

V. CONCLUSIONS

In this work, we study the performance of a tiered network where a tier of macro-cell wireless backhaul nodes is overlaid with a tier of cache-aided SCs. To frame a more realistic scenario, we consider both static and dynamic UT association policies. Building on random spatial models and using tools from stochastic geometry, we derive closed-form expressions and bounds for the achievable ASE for both UT association policies. Then, we study suitable optimization problems to identify deployment strategies for maximizing the ASE and EE of the network. Our findings show that the ASE is maximized by first increasing the SC density and then the storage size at the SCs, whereas an opposite approach should be adopted for maximizing the EE.

APPENDIX I PROOF OF THEOREM 1

Due to space limitations, we provide only a sketch of the proof. Considering static UT association, we have P (S) suc (θ) = P hit P[SINR u(x) > θ] + (1 -P hit )P[SINR u(x) > θ]P[SINR x > θ] and, building on [START_REF] Maso | Cacheaided full-duplex small cells[END_REF], [START_REF] Andrews | A tractable approach to coverage and rate in cellular networks[END_REF], we obtain Ix (s) defined in ( 9)- [START_REF] Andrews | A tractable approach to coverage and rate in cellular networks[END_REF]. Likewise, for dynamic UT association, we have P where (cf. [START_REF] Maso | Cacheaided full-duplex small cells[END_REF], [START_REF] Andrews | A tractable approach to coverage and rate in cellular networks[END_REF])

P[SINR u(x) > θ] = exp -θ σ 2 ρSC R α UT L (S) I u(x) θρ -1 SC R α UT 1 ∆x + exp -θ σ 2 ρBH R α BH L (S) I u(x) θρ -1 BH R α BH 1 / ∆ x (22) with L (D)
I u(x) (s) defined in (12). REFERENCES

Figure 1 .

 1 Figure 1. System model for static and dynamic UT association.

Figure 2 .

 2 Figure 2. Achievable ASE as a function of the cache hit probability.

2 Figure 3 .

 23 Figure 3. Feasible solution set of problems (P1)-(P2) for different values of the monetary budget c.

1 Figure 4 .

 14 Figure 4. Achievable ASE as a function of the SC density.

2 Figure 5 .

 25 Figure 5. EE as a function of the SC density.P[SINR u(x) > θ] = exp -θ σ 2 ρSC R α UT L (S) Ix θρ -1 SC R α UT

  (θ) = P[SINR u(x) > θ]

For more involved types of small-scale fading, e.g., Nakagami-m fading, we refer to[START_REF] Atzeni | Downlink cellular network analysis with LOS/NLOS propagation and elevated base stations[END_REF].


Sophie Arborio

Emmanuelle Simon

Stéphane Auvin

Anne De Saint

FAM-WEST. West Syndrom: building knowledge and singularity of the families experiences

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

FAM-WEST

West Syndrom: building knowledge and singularity of the families experiences S. Arborio (MCU HDR, Université de Lorraine, Crem) -sophie.arborio@univ-lorraine.fr E. Simon (MCU, Université de Lorraine, Crem) -emmanuelle.simon@univ-lorraine.fr Pr. S. Auvin (PU-PH, Hôpital universitaire Robert-Debré, Inserm) -stephane.auvin@rdb.

Preliminary Results

Knowledge as a necessary exchange A process as the interface of knowing and "knowing how" Knowing how to do with uncertainties ? Experiencing the knowledge

Emphasize the importance of experience in building mothers' knowledge of care.

Conclusion

How to build patient educational programs?

Emphazing through a diachronic process the families experiences ("knowing how") instead of only the medical "knowing that" help to adapt educational programs to families daily care.

Enfants de West, ASTB, Injeno Purpose Methods A multidisciplinary research

	How knowledge and "know-how"
	are constructed by families in the case of West's syndrome (WS)?	• 30 families with children ages from 11 months to 15 years
		Analysis		• 60 semi-structured interviews with the mothers
	Path / disability of illness	Complex routes therapeutic and diagnostic	Life story of families	• Participating observations at the hospitals and at the families (Facebook Enfants de West) through the discussion forum • Exchanges analysis homes (18 months)

aphp.fr Pr. A. De Saint Martin (PH, Hôpital de Hautepierre) -anne.desaintmartin@chru-strasbourg.fr Associations Les


Kevin Molloy

Laurent Denarie

Marc Vaisset

Thierry Simeon

Juan Cortés
email: juan.cortes@laas.fr

Thierry Sim

Juan Cort És

Juan Cort

Simultaneous System Design and Path Planning: A Sampling-based Algorithm

Keywords: robot motion planning, sampling-based algorithms, computational biology, protein design

HAL is

Introduction

System design and path-planning problems are usually treated independently. In robotics, criteria such as workspace volume, workload, accuracy, robustness, stiffness, and other performance indexes are treated as part of the system design [START_REF] Vijaykumar | Geometric optimization of serial chain manipulator structures for working volume and dexterity[END_REF][START_REF] Gosselin | A global performance index for the kinematic optimization of robotic manipulators[END_REF][START_REF] Merlet | Optimal design of robots[END_REF]. Path planning algorithms are typically applied to systems with completely fixed geometric and kinematic features. In this work, we propose an extension of the path planning problem, in which some features of the mobile system are not fixed a priori, and can be selected among a finite number of combinations. The goal is to find the best design (i.e. values for the variable features) to optimize the motion between given configurations.

A brute-force approach to solve this problem would involve individually solving motion planning problems for all possible designs, and then selecting the design providing the best result for the (path-dependent) objective function. However, because of the combinatorial explosion, only simple problems involving a small number of variable design features can be treated using this naïve approach. We propose a more sophisticated approach that simultaneously considers system design and path planning. A related problem is the optimization of geometric and kinematic parameters of a robot to achieve a given end-effector trajectory, usually referred to as kinematic synthesis [START_REF] Mccarthy | Kinematic synthesis[END_REF], or the reachability of desirable goal regions in cluttered workspaces [START_REF] Baykal | Optimizing design parameters for sets of concentric tube robots Prepared using sagej.cls using sampling-based motion planning[END_REF][START_REF] Baykal | Asymptotically optimal design of piecewise cylindrical robots using motion planning[END_REF]. Nonetheless, the problem we address in this work (see Section 2.1 for details) is significantly different, since we assume that all kinematic parameters and part of the geometry of the mobile system are provided as input. The design concerns a discrete set of features that can be associated to the bodies of the mobile system, such as shape or electrostatic charge, aiming to find the best possible path between two given configurations provided a path cost function. Very few works have considered such a hybrid design and path planning problem. One of the rare examples is a recently proposed method for aerial vehicle path planning [START_REF] Rudnick-Cohen | Multiobjective design and path planning optimization of unmanned aerial vehicles[END_REF] where the optimal path planning algorithm considers several possible flying speeds and wing reference areas to minimize path risk and time. Since the considered configuration space is twodimensional, the proposed solution is based on an extension of Dijkstra's algorithm working on a discrete representation of the search-space. This type of approach cannot be applied in practice to higher-dimensional problems, as the ones we address.

Sampling-based algorithms have been developed since the late 90s for path planning in high-dimensional spaces [START_REF] Kavraki | Probabilistic roadmaps for path planning in high dimensional configuration spaces[END_REF][START_REF] Lavalle | Planning Algorithms[END_REF], which are out of reach for deterministic, complete algorithms. Our work builds on this family of algorithms, which we extend to treat a combinatorial component in the search-space, associated to the systems design, while searching for the solution path.

Our approach presents some similarities with methods that extend sampling-based path planning algorithms to solve more complex problems such as manipulation planning [START_REF] Siméon | Manipulation planning with probabilistic roadmaps[END_REF]), multi-modal motion planning [START_REF] Hauser | Multi-modal motion planning in non-expansive spaces[END_REF] or minimum constraint removal planning [START_REF] Hauser | The minimum constraint removal problem with three robotics applications[END_REF], which also involve search spaces with hybrid structures. As in these other works, the proposed algorithm simultaneously explores multiple sub-spaces, aiming to find solutions more efficiently. Nevertheless, the hybrid design problem addressed here is different.

This paper presents the Simultaneous Design and Path-planning algorithm (SDAP), which is based on the Transition-based Rapidly-exploring Random Tree (T-RRT) algorithm [START_REF] Jaillet | Samplingbased path planning on configuration-space costmaps[END_REF]. As explained in Section 2.3, the choice of T-RRT as a baseline is guided by the type of cost function we apply for the evaluation of path quality. Nevertheless, other sampling-based algorithms can be extended following a similar approach. In particular, algorithms with optimality properties such as other variants of Rapidly-exploring Random Tree (RRT*) [START_REF] Karaman | Sampling-based Algorithms for Optimal Motion Planning[END_REF], Stable Sparse-RRT* (SST*) [START_REF] Li | Asymptotically optimal sampling-based kinodynamic planning[END_REF] or Fast Marching Tree (FMT*) [START_REF] Janson | Fast marching tree[END_REF]) could be possible alternatives for some classes of problem.

We demonstrate the good performance of the method on relatively simple, academic examples (Section 4). These simple examples allow us to apply the naïve exhaustive method, whose results can be used as a reference to evaluate the performance and the quality of the solutions produced by the SDAP algorithm. Results show that SDAP is able to find the best path-design pairs, requiring much less computing time than the naïve method. This advantage increases with the complexity of the problem.

Although the application of the proposed approach to problems of practical interest is out of the scope of this paper, we note that our motivation comes from problems in computational biology. For more than a decade, robotics-inspired algorithms have been applied to this area [START_REF] Moll | Roadmap Methods for Protein Folding[END_REF][START_REF] Al-Bluwi | Motion planning algorithms for molecular simulations: A survey[END_REF][START_REF] Gipson | Computational models of protein kinematics and dynamics: Beyond simulation[END_REF][START_REF] Shehu | Probabilistic search and optimization for protein energy landscapes[END_REF]. The method presented in this paper aspires to solve problems related to computational protein design [START_REF] Keating | Methods in protein design[END_REF][START_REF] Donald | Algorithms in Structural Molecular Biology[END_REF], which aims to create or modify proteins to exhibit some desired properties. Progress in this field promises great advances in pharmacology, biotechnology, and nanotechnology. Protein design is extremely challenging, and although there have been some considerable strides in the last years [START_REF] Privett | Iterative approach to computational enzyme design[END_REF][START_REF] Tinberg | Computational design of ligand-binding proteins with high affinity and selectivity[END_REF], the problem remains largely open. Current approaches focus on a static picture (i.e. search the amino-acid sequence that stabilizes a given structure), whereas dynamic aspects related to protein function are rarely considered. Our goal behind this work is to develop new methods to optimize functional protein motions. Flexible regions of a protein known as loops play instrumental roles in protein motion, and therefore, protein function. Section 5 presents a preliminary investigation on the application of the proposed method to designing a loop region in order to facilitate the transition of the protein between two stable states. In addition to computational protein design, applications of the proposed approach in robotics can be envisioned and are briefly mentioned in the conclusion (Section 6).

Problem Formulation and Approach

This section defines the problem addressed in this work, along with some notation, and presents an overview of the proposed approach.

Problem Definition

Let us consider an articulated linkage A consisting of n rigid bodies, A 1 ..A n . The kinematic parameters of A are static and are supplied as input. The geometry of the rigid bodies A i can admit some variability, as well as other physical properties (mass, electrostatic charge, ...). More precisely, a discrete set of m design features, f 1 ..f m , is defined and each body A i ∈ A is assigned a design feature f j ∈ F. We denote d as a vector of length n that represents the design features assigned to all the rigid bodies in A, i.e. d defines a particular design. D denotes the set of possible combinations of assignments of features for A, i.e. D defines all possible designs. D is referred to as the design space, which is a discrete space containing m n elements. A given configuration of A is denoted by q. Let C denote the configuration space. Note that for each q ∈ C, only a subset of the possible designs D can be assigned, since some designs are not compatible with some configurations due to self-constraints or environment constraints.

The workspace of A is constrained by a set of obstacles O i ∈ O. Note that the problem formulation could be extended to the design of the workspace, by associating a discrete set of design features to each obstacle. Nevertheless, for the sake of clarity, we consider in this paper that the obstacles have fixed features.

C d

f ree denotes all valid collision-free configurations of A for a given vector d of design features. A path P connecting two configurations q init and q goal of A with design d is defined as a continuous function P : [0, 1] → C, such that P (0) = q init and P (1) = q goal . The path is said to be collision-free if ∀t ∈ [0, 1], P (t) ∈ C d f ree . C f ree is the union of all individual C d f ree :

C f ree = d∈D C d f ree .
P f ree denotes the set of all feasible, collision-free paths connecting q init to q goal , considering all possible designs (∀d ∈ D).

A cost function c : C f ree × D → R associates to each pair {q, d} a cost value, ∀q ∈ C f ree and ∀d ∈ D. We define a costmap as the representation of the cost function on the configuration space of A for each design. Another cost function c P : P f ree × D → R is defined to evaluate the quality of paths. In this work, the path cost function c P is itself a function of the configuration cost function c, i.e. c P is a functional. More precisely, we consider the mechanical work criterion as defined in [START_REF] Jaillet | Samplingbased path planning on configuration-space costmaps[END_REF][START_REF] Devaurs | Optimal path planning in complex cost spaces with sampling-based algorithms[END_REF] to evaluate paths, which aims to minimize the variation of the configuration cost c along the path. This criterion is a suitable choice to evaluate path quality in many situations [START_REF] Jaillet | Samplingbased path planning on configuration-space costmaps[END_REF], and is particularly relevant in the context of molecular modeling. Nevertheless, other cost functions can be considered, such as the integral of c along the path. A discrete approximation, with constant step size δ = 1/l, of the mechanical work (MW) cost of a path P for a system design d can be defined as:

c P (P, d) = l k=1 max 0 , c P k l , d -c P k -1 l , d .
(1) The goal of our method is to find the best pair {P * , d * } such that: (2)

Illustrative examples

Figure 1 presents a simple example to help understand the problem formulation. In this case, the articulated mechanism A involves three bodies, in addition to the base and the end effector. The are three possible choices (design features f 1 , f 2 , f 3 ) for the shape of each body A i . Thus, the number of possible designs is 3 3 = 27. A cost function c is defined as the inverse of the minimum distance between the robot bodies and the obstacles. With the aim of favoring motions of the robot far from obstacles, c P can be defined as the integral of c along the path between q init and q goal . In this simple example, it is easy to figure out that the optimal design is d

* = [f 1 , f 2 , f 1 ]. f 3 q q goal 1 f init f 2 Figure 1.
Simple illustrative example of a simultaneous system design and path planning problem. There are three possible choices (design features) for each body of the robot. The best design, which optimizes a clearance-based cost function, is represented in blue.

Another simple example, presented in Section 4.1, will be used for a first proof of concept. In this other example, the design features do not involve the shape of the mobile system, but the electrostatic charge of each body. The objective is to minimize the cumulative variation of a potential energy function along the solution path.

Approach

A naïve approach to solve the problem would be to compute the optimal cost path for each design d ∈ D, and then choose the optimal design d * that minimizes c P . Such a brute-force approach can be applied in practice to simple problems involving a small number n of variable bodies and/or a few m design features (recall that the design space is size m n ). The method proposed below aims to solve the problem much more efficiently by combining both the discrete (design) and continuous (path) optimization in a single stage.

We assume that, for most problems of interest, the configuration space C is high-dimensional, so that exact/complete algorithms cannot be applied in practice to solve the path-planning part of the problem. For this, we build on sampling-based algorithms (LaValle and Kuffner 2001; [START_REF] Kavraki | Probabilistic roadmaps for path planning in high dimensional configuration spaces[END_REF], which have been very successful in the robotics community since the late 90s, and which have also been applied in other areas such as computational biology [START_REF] Moll | Roadmap Methods for Protein Folding[END_REF][START_REF] Al-Bluwi | Motion planning algorithms for molecular simulations: A survey[END_REF][START_REF] Gipson | Computational models of protein kinematics and dynamics: Beyond simulation[END_REF][START_REF] Shehu | Probabilistic search and optimization for protein energy landscapes[END_REF]. We also assume that the cardinality of the design space D is moderately high, such that a relatively simple combinatorial approach can be applied to treat the design part of the problem.

The idea is to explore C f ree to find paths between q init and q goal simultaneously considering all possible deigns d ∈ D. To reduce the number of configuration-design pairs {q, d} to be evaluated during the exploration, it is important to apply an effective filtering strategy. The choice of the particular sampling-based path planning algorithm and filtering strategy mainly depend on the type of objective function c P being considered. The approach described below has been developed to find good-quality solutions with respect to the MW path evaluation criterion (1). In this work, we extend the T-RRT algorithm [START_REF] Jaillet | Samplingbased path planning on configuration-space costmaps[END_REF], which finds paths that tend to minimize cost variation by filtering during the exploration tree nodes that would produce a steep cost increase. Following a similar approach, alternative algorithms and the associated filtering strategies could be developed to optimize other path cost functions. For instance, variants of RRT* [START_REF] Karaman | Sampling-based Algorithms for Optimal Motion Planning[END_REF] or FMT* [START_REF] Janson | Fast marching tree[END_REF]) could be considered for optimizing other types of monotonically increasing cost functions.

Algorithm

This section presents the SDAP algorithm, building on a single-tree version of T-RRT. However, the approach is directly applicable to multi-tree variants [START_REF] Devaurs | A multitree extension of the transition-based RRT: Application to ordering-and-pathfinding problems in continuous cost spaces[END_REF]. First, the basic algorithm is introduced, followed by additional explanation of the tree extension strategy and a brief theoretical analysis.

SDAP Algorithm

The SDAP pseudo-code is shown in Algorithm 1. A search tree, T , is created with q init as the root node. The tree is grown in configuration space through a series of expansion operations. Each node s in T encodes a configuration q and a set of designs D ⊆ D for which the configuration is valid, and which have not been filtered out in previous iterations for the construction of corresponding branches of the tree. In other words, each node's set of designs D is a subset of its parent's designs, i.e. Designs(s) ⊆ Designs(P arent(s)). This is essential to reduce the number of designs evaluated during the exploration. The algorithm keeps iterating until the goal configuration q goal is reached for at least one of the designs, or a maximum number of iterations is performed.

At each iteration, a random configuration q rand is generated (line 3). In T-RRT, a new node q new is created by expanding the nearest node in T q near in the direction of q rand for a distance δ. q new is then conditionally added to T based on a transition test. A common heuristic for this transition test is the Metropolis criterion, traditionally Algorithm 1: SDAP Algorithm input : the configuration space C; the design space D; the cost function c; the start state qinit; the goal state q goal ; maximum number of iterations MaxIter output: the tree T 1 T ← InitTree(qinit, D) 2 while not StoppingCriterion (T , q goal ,MaxIter) do applied in Monte Carlo methods [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF]. Transitions to lower cost nodes are always accepted and moves to higher costs nodes are probabilistically accepted. The probability to transition to a higher cost node is controlled by a temperature variable T . T-RRT dynamically controls T , as explained in the next section.

SDAP modifies the expansion and transition test functions of the standard T-RRT algorithm in order to address the design and path planning problems simultaneously. At each iteration, SDAP attempts to expand at least one node per design in D. The process is shown in Figure 2, where each design d ∈ D is encoded as a color on each node of the tree. During each iteration, SDAP expands a set of nodes that covers all designs D. In other words, the NearestNeighbors function (line 4) returns a set, Neighbors, containing the closest node to q rand for each design. In Figure 2, q rand is shown in black and the 3 nodes in the set Neighbors are circled in red. Each node in Neighbors is extended towards q rand (lines 6 -10) creating new candidate nodes which are labeled s 1 , s 2 , and s 3 in Figure 2. All 3 designs in s 1 fail the transition test, so the new node is not added to T . For s 2 the blue design passes the transition test and the node is added to T . Finally for s 3 , 1 of the 3 designs (yellow) fails the transition test, resulting in a node with 2 designs being added to T . Such a filtering of designs during the construction of the exploration tree is the key to SDAP's good performance.

Controlling Tree Expansion

The transition test is governed by the temperature parameter T . T-RRT automatically adjusts T during the exploration and has been shown highly effective in balancing tree exploration and tree refinement [START_REF] Jaillet | Samplingbased path planning on configuration-space costmaps[END_REF]. At each iteration, T-RRT adjusts T by monitoring the acceptance rate of new nodes. SDAP extends this idea by maintaining a separate temperature variable T (d) for each design d ∈ D. A given design d can appear in multiple nodes in Neighbors.

For each design d, the node in Neighbors closest to q rand is identified. The temperature T (d) is adjusted based on the success or failure of the extension operation from this node for design d.

The pseudocode for the transition test function is shown in Algorithm 2. The Neighbors set is processed in ascending order of the distance of each node from q rand (line 6 of Algorithm 1). For the node being expanded (s near ), each design d has its cost evaluated (line 4). Transitions to lower cost nodes are always accepted (line 8). Transitions to higher cost nodes are subjected to probabilistic acceptance (line 10). The set V (lines 12 and 18) tracks designs which have had their temperature adjusted during this iteration. The function returns the set D of designs that pass the transition test.

Theoretical Analysis

In this section we provide some theoretical analysis of the SDAP algorithm's completeness and path optimality. A theoretical analysis of the complexity of SDAP with respect to the brute-force approach is difficult, since both are stochastic processes. In this work, we instead provide empirical results in Section 4 that clearly show SDAP's efficiency versus an exhaustive search of paths for all possible designs.

Probabilistic

Completeness: SDAP's probabilistic completeness directly derives from that of RRT (LaValle and Kuffner 2001), which is inherited by T-RRT under the condition to guarantee a strictly positive probability of passing the transition test as explained in [START_REF] Jaillet | Samplingbased path planning on configuration-space costmaps[END_REF]. Since SDAP maintains this property by incorporating temperatures in the transition test for each given design d ∈ D, it also ensures the positive transition probability and that each C d f ree will be completely sampled,

T (d) ← T (d) / 2 (∆c / energyRange(T ,d))
15 else 16

T (d) ← T (d) • 2 T rate 17 if success then D ← D ∪ d 18 V ← V ∪ d 19 return (D)
thus maintaining the probabilistic completeness of the algorithm.

Path Optimality:

The current SDAP implementation is based on T-RRT, which has been empirically shown to compute paths that tend to minimize cost with respect to the MW criterion [START_REF] Jaillet | Samplingbased path planning on configuration-space costmaps[END_REF]), but without theoretical guarantee of optimality. Using anytime variants of T-RRT (AT-RRT or T-RRT * ) [START_REF] Devaurs | Optimal path planning in complex cost spaces with sampling-based algorithms[END_REF] would provide an asymptotic convergence guarantee. Implementing these within SDAP remains as future work.

Empirical Analysis and Results

As a proof of concept, we apply SDAP to a set of academic problems. SDAP is implemented as an adaptation of the Multi-T-RRT algorithm [START_REF] Devaurs | A multitree extension of the transition-based RRT: Application to ordering-and-pathfinding problems in continuous cost spaces[END_REF], with two trees growing from the initial and goal configurations. The search stops when the algorithm is able to join the two trees. For each problem, SDAP is compared against a naïve approach consisting of multiple independent runs of Multi-T-RRT on each design d ∈ D.

Test System Description

The test system is a 2D articulated mechanism with a fixed geometry surrounded with fixed obstacles. The bodies

c(q, d) = LJ(q, d) + ES(q, d) (3)
with:

LJ(q, d) = |A|-2 i=1   |A| j=i+2 2 • R A i A j 12 - 2 • R A i A j 6   + |A| i=1   |O| j=1 2 • R A i O j 12 - 2 • R A i O j 6   (4) ES(q, d) = |A|-2 i=1   |A| j=i+2 f i • f j A i A j   + |A| i=1   |O| j=1 f i • g j A i O j   (5)
where X i X j represents the Euclidean distance between the centers of the bodies/obstacles X i and X j . SDAP is empirically tested using several scenarios with the objective to find the optimal path-design pair (P * , d * ). First, a simple 4-body system is used to illustrate the effect of the design choice and to discuss the optimality of the solution. Then, two more complex systems, one involving a variable number of bodies to be designed and the other involving a variable number of degrees of freedom (DOF), are used to analyze the effect of the complexity of the problem on the performance of the algorithm.

4-Body System:

The first system consists of four bodies and five obstacles as shown in Figure 3 (top). q init and q goal correspond to fully stretched configurations, to the left and to the right, represented with solid and dashed outlines respectively. The design space consists of 3 4 = 81 possible combinations and the configuration space is 3 dimensional. Owing to the position of the positively and negatively charged obstacles, this scenario favors designs involving negative charges. However, the uncharged obstacles at the top and bottom of the workspace create a narrow passage that enforces the articulated chain to fold up, which penalizes designs with negative charges everywhere, owing to their repulsion. Therefore, this is a non-trivial design problem.

4.1.2 10-Body System: A larger system with 10 bodies and six obstacles is shown in Figure 3 (middle). The design space D contains 3 10 = 59049 possibilities, which cannot be exhaustively explored within a reasonable computing time, and is also challenging for SDAP because of memory issues (further discussed below). For that reason, three simplified versions of this scenario are constructed. The first one fixes the design for the first seven bodies A 1 ..A 7 as UUUUUUU. The remaining bodies (A 8 , A 9 , and A 10 ) can be designed, resulting in a design space of 3 3 = 27 designs. The second version expands the design space to the last 4 bodies (A 7 ..A 10 ) , resulting in a design space of 3 4 = 81 designs. The last case expands the design space to the last 5 bodies (A 6 ..A 10 ), resulting in a design space of 3 5 = 243. In all cases, the configuration space is 9 dimensional. All versions of the 10-body system are constrained with the same obstacles. They were chosen so that designs with strongly positively or negatively charged end-effectors will be trapped at local minima resulting from attractive or repulsive forces generated by the bottom obstacles.

30-Body System:

A third scenario, presented in Figure 3 (bottom), involves a mobile system with up to 30 bodies. This example is used to analyze the influence of the dimension of the configuration space on the performance of SDAP. For this, a variable number of articulated joints is considered. More precisely, rotational joints are considered for all the bodies, or only for one every 2, 3 and 6 bodies. Thus, the dimension of C is 30, 15, 10 and 5, respectively. As for the 10-body system, the design of the first bodies in the chain is fixed to U. Only the last 4 bodies (A 27 ..A 30 ) are designed.

Benchmark Results

We compare SDAP to a naïve approach (solving individual problems for each design d ∈ D) using the same Multi-T-RRT implementation. In other words, we compare one run of the SDAP algorithm against |D| runs of a singledesign path search. Multiple runs are performed (100 for the 4-body scenario, 50 for the 10-body scenario with 3 designed bodies, 20 for the 10-body scenario with 4 and 5 designed bodies, and 10 for the variants of the 30-body scenario) to reduce statistical variance inherent with stochastic methods. The single-design explorations can spend time trying to escape local minima associated with the costmaps of unfavorable designs, causing very long execution times and high-cost paths. As we are not interested in finding a solution path for every possible design but only for the designs with low-cost paths, a timeout is enforced for the single-design explorations of 5 minutes for the 4-body scenario, and 15 minutes for the 10body scenarios. For the 30-body problem, the differences in the dimensionality of the configuration space mandated different timeouts of 15, 30, 30, and 35 minutes for the 5, 10, 15, and 30 DOF scenarios respectively. No timeout was considered for SDAP (i.e. the algorithm always provided a solution). In all the cases, the T-RRT parameter T rate was set to 0.1, being the initial temperature T = 1 and the Boltzmann constant K = 0.00198721 (note that T and K do not have much physical meaning in the present context). During the exploration, the maximal values reached by T were around 100, being slightly lower for SDAP compared to the naïve approach. All the runs were performed in a single threaded process on a Intel(R) Xeon(R) CPU E5-2650 @ 2.00GHz processor with 32GB of memory.

Results for the 4-Body Scenario:

We did a simple experiment aiming to show how strongly the choice of the design may affect configuration-space costmap. Figure 4 for two designs along with a solution path. One design has a negatively charged end-effector (UUUN) and one has a positively charged end-effector (UUUP). Angles 1 and 2 are projected onto the x and y axis respectively, with angle 3 being set to minimize the cost function. For the UUUN design, the costmap is highly favorable to the desired motion, starting at a high-cost area and proceeding downhill to a low cost area. The costmap associated with the UUUP design shows a non-favorable motion between the two states.

Figure 5 presents results concerning the optimal design obtained by the naïve approach and by SDAP. Both algorithms provide very similar results (see the enlarged region in the figure). The best designs mainly involve negative and uncharged bodies, the best one being UNNU. As expected, owing to repulsive forces between the bodies of the mobile system, designs involving more than two nonconsecutive negative bodies are worst than those containing only one or two negative charges. In terms of computational performance, the naïve approach required 2 ± 0.1 minutes to compute the paths for the 81 designs individually, whereas SDAP required 1 ± 0.6 minutes. This represents a speed-up of 2 in favor of SDAP. As shown for the more complex examples below, the performance gain significantly increases when the dimensions of the design space and the configuration space increase. Regarding the number of nodes in the T-RRT trees, the 81 runs for the single-design searches generated approximately 77,000 nodes in total (i.e. around 1,000 nodes per design on average). SDAP generated approximately 39,000 nodes for the simultaneous exploration of all the designs, which is directly reflected with the increase in computational efficiency. Note that the high density of nodes created by the SDAP algorithm could be well exploited to improve the path cost by incrementally adding cycles in an anytime fashion, as proposed by [START_REF] Devaurs | Optimal path planning in complex cost spaces with sampling-based algorithms[END_REF]. Such an extension of the algorithm, which remains for future work (see Section 6), would guarantee asymptotic convergence to the global optimum.

Influence of the cardinality of the design space:

We used the 3 variants of the 10-body scenario to analyze the influence of the cardinality of the design space on the performance of the SDAP algorithm. Results are summarized in Table 1 and Figure 6. The speedup reported in the right column of the table clearly show that the difference in computing time between SDAP and the naïve approach increases significantly compared to the 4-body scenario. In the 3-designed-body problem, SDAP is 26 times faster than the single-design search on average. In the 4 and 5-designed-body cases, SDAP is 47 and 61 times faster respectively. While the cardinality of design space D is multiplied by 3 between each of the versions, the average execution time of SDAP only grows by a factor of 2.5 and 2.2 respectively. SDAP achieves this result by limiting the growth of the T-RRT search tree with respect to the size of D. The numbers in Table 1 show that the naïve approach yields a search tree scaling closely with the size of D, which directly correlates with higher computational requirements. D), and the size of the design space (D) for experiments allowing 3, 4, and 5 design positions. The average size of the search tree (expressed in thousands) and the wall clock execution time (minutes) are shown along with their standard deviations. The design space growth factor (GF) is shown along with the GF for each average tree sizes and run times (growth from the 3 design position problem). The tree size for the naïve approach follows the growth in the design space, where the SDAP tree grows at a much slower rate. The final column shows the speed up factor between the naïve approach and SDAP. Importantly, as in the case of the 4-body problem, SDAP successfully identifies the designs yielding the best path costs. Figure 7 compares the solutions found by SDAP to those obtained by the naïve approach.

|D| |D|

Note that our initial implementation of SDAP employs a brute force nearest neighbor search with O(n) complexity. The number of average nodes shown in Table 1 for the naïve approach represents the total of a single set of runs to cover each design in D. However, each individual execution produces a search tree containing 16,000 nodes on average for the 5-designed-body problem. The SDAP search tree for the same problem contains 432,000 nodes on average. Thus, a more efficient nearest neighbor search method would certainly further improve SDAP's relative performance.

The memory requirement for SDAP is directly correlated to the size of the search tree. Each node in the search tree takes on average 37 kB of memory for the 10body scenarios. For the 5-designed-body problem, SDAP required approximately 16 GB of memory (this also involves base data structures used for fast collision detection and distance computation, in addition to the system representation). Extrapolating from the data in Table 1, each additional body in the design space doubles the size of the search tree, and thus, doubles the memory requirements. We further discuss this in Section 6, and suggest some ideas to circumvent this issue.

Influence of the dimension of the configuration space:

The 4 variants of the 30-body scenario, involving a variable number of DOF, were used to analyze the influence of the dimensionality of the configuration space on the performance of the SDAP algorithm. Results are summarized in Table 2. The numbers in this table (see also Figure 8) show that performance gain of SDAP with respect to the naïve approach increases with the complexity of the problem. The improvement is very significant at the beginning, when changing from 5 DOF to Cost of the paths for the best designs found with the naïve approach (dotted blue and dashed green lines) and SDAP (red line) for the three variants of the 10-body problem: 3-designed-body problem (top), 4-designed-body problem (middle), 5-designed-body problem (bottom). SDAP solutions shown only for discovered paths (does not exhaustively search). SDAP discovers the same low-cost designs as the exhaustive single-design searches.

10 DOF, and tends to decrease for larger dimensions. Note however that the speedups reported in higher dimensions are underestimated due to use of a timeout for the singledesign runs of the naïve approach. As with the previous experiments, the best designs obtained by SDAP remain very similar to those obtained with the naïve approach.

As illustrated in Figure 9, which summarizes results about the influence of the problem complexity on the SDAP performance versus the naïve approach, the speedup tends to increase with both the cardinality of the design space and the dimension of the configuration space. It is important to Table 2. For the 30-Body scenario, the table shows results for different dimensions of the configuration-space, involving systems with 5, 10, 15 and 30 DOF. The average size of the search tree (expressed in thousands) and the wall clock execution time (minutes) are shown along with their standard deviations. The configuration space growth factor (GF) is shown along with the GF for each average tree sizes and run times (growth from the 5 DOF problem). The run times for SDAP grow at a slower rate than for the naïve approach, as shown in the speed up factors of the last column. The green line shows this with respect to the size of design space using the 10-body problem as benchmark. The blue line showcases the performance varying the number of DOF using the 30-body scenario.

mention that, for the 10-body problem, on average the naïve approach reached the timeout 15%, 23%, 30% of the time for the 3, 4, and 5-designed body problems respectively. For the 30-body problem, the timeout was reached 1%, 31%, 36% and 53% of the time for the 5, 10, 15, and 30 DOF experiments. This large number of failures when computing paths with the naïve approach means that many runs would have needed a longer time to eventually find a solution. Thus, larger timeout values would have certainly led to an even more impressive speedup in favor of SDAP.

A First Application to Protein Motion Design

Results presented in the previous section show the good performance of SDAP to solve relatively simple problems. This section deals with a more difficult problem related to protein motion design. First, the problem is formulated as a simultaneous design and path planning problem. Some simplifications are required in order to enable the application of the current version of the SDAP algorithm. A first application in this context is attempted for the design of protein region (a loop) aiming to facilitate its transition between two stable states.

Problem Definition

The protein motion design problem can be formulated as follows: given two protein scaffold conformations * q 1 and q 2 , find a sequence of n amino acids with (meta-)stable states for both structures, and favoring the conformational change (i.e. the motion) between them. Although existing protein design methods may identify sequences for which the two given conformations (q 1 and q 2 ) correspond to (meta-)stable states, these methods do not consider the feasibility of the transition path between these two states. In order to design a candidate protein for this motion, it is necessary to find a tuple (path between q 1 and q 2 , sequence of amino acids) that minimizes energy variation (i.e. the mechanical work) along the transition path.

For this preliminary application of the SDAP algorithm to protein motion design, we consider the coarse-grained model proposed by [START_REF] Brown | Coarsegrained sequences for protein folding and design[END_REF] (usually called BLN model) to limit the combinatorial complexity of the design space. More precisely, the 20 natural amino acids are grouped into 3 classes according to their chemical properties: B for the hydrophobic, L for the hydrophilic, and N for the neutral (see Table 3). This simplification reduces the cardinality of the design space to 3 n possible designs (instead of 20 n ). Furthermore, each amino acid (i.e. each body A 1 , . . . , A n ) is modeled as a single bead centered on its C α atom. Consecutive beads are linked by virtual bonds. This simplification allows us to neglect additional degrees of freedom associated with the protein side-chains. The configuration cost function corresponds to the potential energy, defined as:

c(q, d) = θ 1 2 k θ (θ -θ 0 ) 2 + φ A(1 + cos φ) + B(1 -cosφ) + C (1 + cos 3φ) + D 1 + cos φ + π 4 + i,j≥i+3 4S 1 σ r ij 12 -S 2 σ r ij 6 (6)
In this unitless equation, the first term is a bond angle energy term, where θ is the angle formed by two consecutive virtual bonds, k θ is a force constant of 20 rad -2 , and θ 0 is a reference angle of 105 • . The second term is a dihedral angle energy term, where φ is the dihedral angle formed by three consecutive virtual bonds, σ is the distance unit (σ = 1 Å), and A, B, C, D are parameters varying as a function of the local secondary structure (computed using DSSP [START_REF] Kabsch | Dictionary of protein secondary structure: pattern recognition of hydrogen-bonded and geometrical features[END_REF][START_REF] Touw | A series of PDBrelated databanks for everyday needs[END_REF])):

• A = 0, B = C = D = 1.
2 for an helical structure,

• A = 0.9, B = D = 0, C = 1.2 for a strand,

• and A = B = D = 0, C = 0.2 for other structures.

The last term considers pair-wise non bonded interactions, where r ij is the distance between the i th and the j th C α atoms in the sequence, and S 1 and S 2 are parameters that depend on the designed features associated to the pair of amino acid residues: S 1 = S 2 = 1 for B-B interactions, S 1 = 1 3 and S 2 = -1 for L-L and L-B interactions, and S 1 = 1 and S 2 = 0 for N-L, N-B, and N-N interactions.

Experiments

5.2.1 Test System: For this preliminary work, we have investigated Escherichia coli dihydrofolate reductase (ecDHFR). Prior work showed that the loop formed by residues 9-24 (called the Met20 loop) fluctuates between two conformations (called closed and occluded in related work) at a rate comparable with the speed of the reaction catalyzed by this enzyme [START_REF] Falzone | Dynamics of a flexible loop in dihydrofolate reductase from Escherichia coli and its implication for catalysis[END_REF], indicating that this flexible loop is implicated in this reaction. Other studies also demonstrated that the replacement of central residues of the Met20 loop (Met16, Glu17, Asn18, and Ala19) by glycine residues results in a striking decrease in the enzyme activity [START_REF] Li | Functional role of a mobile loop of Escherichia coli dihydrofolate reductase in transition-state stabilization[END_REF][START_REF] Osborne | Backbone Dynamics in Dihydrofolate Reductase Complexes: Role of Loop Flexibility in the Catalytic Mechanism †[END_REF]. Our aim is to study the effect of these mutations on the protein loop motion using SDAP.

Experiment Setup:

Structures of ecDHFR corresponding to the closed and occluded conformational states of the Met20 loop have been solved experimentally and are available in the Protein Data Bank with IDs 1RX2 and 1RX4, respectively. These two structures, which are the starting point for our studies, are illustrated in Figure 10. A comparison of these two structures shows that only the Met20 loop moves significantly, while the rest of the protein backbone is very rigid. Therefore, to reduce the dimension of the configuration space to be explored, all the protein amino acid residues excepting those of the Met20 loop (residues 9-24) are considered as rigid bodies and are fixed in the space. This simplified model involves 32 DOFs. To keep the two ends of the loop connected to the rest of the protein, kinematic loop closure constraints need to be satisfied during the exploration performed by SDAP. For this, we have applied the approach proposed by [START_REF] Cortés | Sampling-based motion planning under kinematic loop-closure constraints[END_REF] to extend sampling-based planners to closed-chain mechanisms. This approach, based on a decomposition of the degrees of freedom together with the combination of geometrically-constrained random sampling and inverse kinematics, was successfully applied to protein loops in previous work [START_REF] Cortés | Geometric algorithms for the conformational analysis of long protein loops[END_REF].

Concerning the design space, only the four central residues (16)(17)(18)(19) are associated to variable design features. The other residues are each constrained to keep the design feature B, L, or N, corresponding to its original amino acid type. This reduces the design space to a size of 3 4 = 81 possible designs. SDAP was applied to find paths from the occluded conformation to the closed one, considering an initial temperature of 300 K, and a maximum variation of bond rotations of 0.015 radians at each step.

Results:

For this difficult problem, SDAP was able to find paths connecting the closed and occluded conformations of the Met20 loop. The algorithm required about 45 minutes of CPU time, and generated around 50,000 nodes. However, it failed to identify the best designs for this motion. Indeed, all the nodes in the solution path were accepted by SDAP for all the designs. Moreover, 93% of the nodes generated by the algorithm were labeled as valid for the 81 possible designs. A more accurate analysis of the energy along the solution path revealed that the design-dependent component of the energy is negligible (about 5 orders of magnitude smaller) compared with the total energy variation induced by the loop motion. These results clearly indicate that the chosen coarse-grained energy function is not well adapted to protein motion design. A more detailed energy function, able to differentiate the quality of the sequences for a conformational transition, is required. However, using more accurate (all-atom) energy function is not possible in practice with the current implementation of SDAP, which does not scale to a sequence space of size 20 n , mainly due to memory issues. As mentioned below, the implementation of an enhanced version of SDAP able to deal with such a complex problem remains a challenging objective for future work.

Conclusion

In this paper, we have presented an original formulation of a challenging problem combining system design and path planning, and have proposed a new approach to solve it, building on sampling-based algorithms. Our first implementation of SDAP already shows significant gains in efficiency and accuracy compared with a brute-force approach.

While our work was primarily motivated by protein motion design problems, several applications in robotics can also be envisioned. In addition to the design of some robot's features to optimize its motion in a given workspace, it would also be possible to apply the proposed method to optimize the workspace layout for a given robot. One can also imagine applications for helping to the design of modular self-reconfigurable robots. We aim to implement SDAP within robot motion planning software to investigate these potential applications further.

For future work, in addition to investigating a variant of SDAP with asymptotic optimality guarantees based on the any-time T-RRT algorithm (AT-RRT) [START_REF] Devaurs | Optimal path planning in complex cost spaces with sampling-based algorithms[END_REF], we also aim to introduce further improvements. First, the exploration of very-high-dimensional configuration and design spaces implies computer memory issues (the resulting trees or graphs are very large). A solution to this problem would be to introduce pruning stages during the exploration, as is done in the SST* algorithm [START_REF] Li | Asymptotically optimal sampling-based kinodynamic planning[END_REF]. Larger design spaces also require SDAP to employ more sophisticated filters and heuristics. A promising direction is to perform statistical learning of the structure of the space, aiming to explore it more efficiently and to control the size of the search tree. To tackle the extreme combinatorial complexity of real-size problems in protein design, the algorithmic improvements described here clearly need to be combined with efficient parallelization strategies of the SDAP algorithm, similar to those proposed for T-RRT (Estaña et al. n.d.) for running on (possibly large) computer clusters.

c

  P (P * , d * ) = min{c P (P, d) | P ∈ P f ree , d ∈ D}.

Figure 2 .

 2 Figure 2. An expansion operation for SDAP. Designs are encoded as colors within each node. The nodes being expanded are circled in red. The expansion towards node s1 fails for all 3 designs, The expansion to s2 succeeds, and the expansion to s3 succeeds for 2 of the 3 designs.

Algorithm 2 :

 2 TransitionTest(T , s near , q new , c) input : the input tree T ; vector of temperatures T; parent node snear; new node qnew; the cost function c; temperature adjustment rate Trate; Boltzmann constant K internal: set of designs V with adjusted temperatures output : vector of designs D that pass the transition test 1 S ← φ 2 for d ∈ Designs(snear) do 3 if CollisionTest(qnew, d) == False then 4 cnear = c(Config(snear), d); cnew = c(qnew, d) exp(-∆c / (K • T (d))) > UniformRand()) then 11 success ← true 12 if d / ∈ V then 13 if success then 14

A

  1 ..A n are circles with radius R. The first body A 1 is a fixed base. The other bodies A 2 ..A n are articulated by a rotational joint centered on the previous rigid body that can move in the interval [0, 2π). A configuration q is described by a vector of n -1 angles corresponding to the value of each rotational joint. The features f 1 ..f n assigned to each body are electrostatic charges in F = {-1, 0, 1} (i.e. m = 3). The design vector d contains n charges f 1 ..f n associated to each rigid body A 1 ..A n of the mechanism. In the following, d will be written as a string, with each charge {-1, 0, 1} corresponding to N, U, and P respectively. For example, the design NPUN corresponds to the vector d = [-1, 1, 0, -1]. Obstacles O 1 ..O k are circles of radius R and have electrostatic charges with predefined values g i ∈ F.The cost function is inspired by a simple expression of the potential energy of a molecular system. It contains two terms, one corresponding to the Lennard-Jones potential and the other to the electrostatic potential. It is defined as:

Figure 3 .

 3 Figure 3. A 4-body (top), 10-body (middle) and 30-body (bottom) scenarios. Obstacles with positive charges (labeled P) are shown in solid red, negative (labeled N) in solid blue, neutral in gray (unlabeled). The initial state is shown in green with a solid line, a transition state shown in blue, and the goal state is shown in red with a dashed line.

Figure 4 .

 4 Figure 4. Configuration-space costmap of the 4-body system projected onto the first two DOF of the system expressed in radians. The initial configuration is indicated by the red dot on the left, and the goal by the green dot on the right. The plot on the left is for the UUUP design, the one to the right for UUUN. Each cell's cost is computed by finding the value of the third angle that minimizes the cost.

Figure 5 .

 5 Figure5. Cost of the paths for the best designs found with the naïve approach (dotted blue and dashed green lines) and SDAP (red line) for the 4-body scenario. SDAP solutions are shown only for discovered paths (does not exhaustively search). SDAP discovers the same low-cost designs as the exhaustive single-design searches.

Figure 6 .

 6 Figure 6. Evolution of the computing time (in logarithmic scale) with respect of the cardinality of the design space for the naïve approach and for SDAP.

  Figure7. Cost of the paths for the best designs found with the naïve approach (dotted blue and dashed green lines) and SDAP (red line) for the three variants of the 10-body problem: 3-designed-body problem (top), 4-designed-body problem (middle), 5-designed-body problem (bottom). SDAP solutions shown only for discovered paths (does not exhaustively search). SDAP discovers the same low-cost designs as the exhaustive single-design searches.

Figure 8 .

 8 Figure 8. Evolution of the computing time (in logarithmic scale) with respect of the dimension of the configuration space for the naïve approach and for SDAP.

Figure 9 .

 9 Figure9. The relative performance of SDAP over the naïve approach. The green line shows this with respect to the size of design space using the 10-body problem as benchmark. The blue line showcases the performance varying the number of DOF using the 30-body scenario.

  B Gly N Asn L Cys B Val B Ser N His L Leu B Trp B Thr L Gln L Ile B Tyr B Glu L Lys L Phe B Pro N Asp L Arg L Table 3. Correspondence of the 20 natural amino acids (columns 20) with the 3 coarse-grained design features (columns 3).

Figure 10 .

 10 Figure 10. Representation of the closed (left) and occluded (left) conformations of the Met20 loop in ecDHFR protein.

Table 1 .

 1 For each variant of the 10-body scenario, the table shows number of design bodies (

			Design GF	Nodes (x 10 3 )	Naïve Node GF	Time (min)	Time GF	Nodes (x 10 3 )	SDAP Node GF	Time (min)	Time GF	SDAP Speedup
	3	27	-	359 ± 34	-	147 ± 24	-	79 ± 30	-	6 ± 4	-	24.5
	4	81	3	1,204 ± 63	3.4	658 ± 62	4.5	205 ± 52	2.6	14 ± 8	2.3	47.0
	5 243	9 3,798 ± 119	10.5 2,157 ± 112	14.7 432 ± 103	5.5 35 ± 12	5.8	61.6

Prepared using sagej.cls

* In this paper, the term conformation is used to define a spatial arrangement of the atoms of the molecule. It is the equivalent to the term configuration for a robot.Prepared using sagej.cls

Funding

This work was partially supported by the French National Research Agency (ANR) under project ProtiCAD (grant number ANR-12-MONU-0015), and by the European Union's Horizon 2020 research and innovation program under (grant number 644271 AEROARMS).