Jean-Bastien Grill
email: jean-bastien.grill@inria.frmichal.valko@inria.fr

Michal Valko

Rémi Munos
email: munos@google.com

Blazing the trails before beating the path: Sample-efficient Monte-Carlo planning

You are a robot and you live in a Markov decision process (MDP) with a finite or an infinite number of transitions from state-action to next states. You got brains and so you plan before you act. Luckily, your roboparents equipped you with a generative model to do some Monte-Carlo planning. The world is waiting for you and you have no time to waste. You want your planning to be efficient. Sample-efficient. Indeed, you want to exploit the possible structure of the MDP by exploring only a subset of states reachable by following near-optimal policies. You want guarantees on sample complexity that depend on a measure of the quantity of near-optimal states. You want something, that is an extension of Monte-Carlo sampling (for estimating an expectation) to problems that alternate maximization (over actions) and expectation (over next states). But you do not want to StOP with exponential running time, you want something simple to implement and computationally efficient. You want it all and you want it now. You want TrailBlazer.

Introduction

We consider the problem of sampling-based planning in a Markov decision process (MDP) when a generative model (oracle) is available. This approach, also called Monte-Carlo planning or Monte-Carlo tree search (see e.g., [START_REF] Kocsis | Bandit-based Monte-Carlo planning[END_REF]), has been popularized in the game of computer Go [START_REF] Coulom | Efficient selectivity and backup operators in Monte-Carlo tree search[END_REF][START_REF] Gelly | Modification of UCT with patterns in Monte-Carlo Go[END_REF][START_REF] Silver | Mastering the game of Go with deep neural networks and tree search[END_REF] and shown impressive performance in many other high dimensional control and game problems [START_REF] Browne | A survey of Monte Carlo tree search methods[END_REF]. In the present paper, we provide a sample complexity analysis of a new algorithm called TrailBlazer.

Our assumption about the MDP is that we possess a generative model which can be called from any state-action pair to generate rewards and transition samples. Since making a call to this generative model has a cost, be it a numerical cost expressed in CPU time (in simulated environments) or a financial cost (in real domains), our goal is to use this model as parsimoniously as possible.

Following dynamic programming [START_REF] Bellman | Dynamic Programming[END_REF], planning can be reduced to an approximation of the (optimal) value function, defined as the maximum of the expected sum of discounted rewards: E t≥0 γ t r t , where γ ∈ [0, 1) is a known discount factor. Indeed, if an ε-optimal approximation of the value function at any state-action pair is available, then the policy corresponding to selecting in each state the action with the highest approximated value will be O (ε/ (1 -γ))-optimal [START_REF] Bertsekas | Neuro-Dynamic Programming[END_REF].

Consequently, in this paper, we focus on a near-optimal approximation of the value function for a single given state (or state-action pair). In order to assess the performance of our algorithm we measure its sample complexity defined as the number of oracle calls, given that we guarantee its consistency, i.e., that with probability at least 1 -δ, TrailBlazer returns an ε-approximation of the value function as required by the probably approximately correct (PAC) framework.

We use a tree representation to represent the set of states that are reachable from any initial state. This tree alternates maximum (MAX) nodes (corresponding to actions) and average (AVG) nodes (corresponding to the random transition to next states). We assume the number K of actions is finite. However, the number N of possible next states is either finite or infinite (which may be the case when the state space is infinite), and we will report results in both the finite N and the infinite case. The root node of this planning tree represents the current state (or a state-action) of the MDP and its value is the maximum (over all policies defined at MAX nodes) of the corresponding expected sum of discounted rewards. Notice that by using a tree representation, we do not use the property that some state of the MDP can be reached by different paths (sequences of states-actions). Therefore, this state will be represented by different nodes in the tree. We could potentially merge such duplicates to form a graph instead. However, for simplicity, we choose not to merge these duplicates and keep a tree, which could make the planning problem harder. To sum up, our goal is to return, with probability 1 -δ, an ε-accurate value of the root node of this planning tree while using as low number of calls to the oracle as possible. Our contribution is an algorithm called TrailBlazer whose sampling strategy depends on the specific structure of the MDP and for which we provide sample complexity bounds in terms of a new problem-dependent measure of the quantity of near-optimal nodes. Before describing our contribution in more detail we first relate our setting to what has been around.

Related work

In this section we focus on the dependency between ε and the sample complexity and all bound of the style 1/ε c are up to a poly-logarithmic multiplicative factor not indicated for clarity. Kocsis and Szepesvári [START_REF] Kocsis | Bandit-based Monte-Carlo planning[END_REF] introduced the UCT algorithm (upper-confidence bounds for trees). UCT is efficient in computer Go [START_REF] Coulom | Efficient selectivity and backup operators in Monte-Carlo tree search[END_REF][START_REF] Gelly | Modification of UCT with patterns in Monte-Carlo Go[END_REF][START_REF] Silver | Mastering the game of Go with deep neural networks and tree search[END_REF] and a number of other control and game problems [START_REF] Browne | A survey of Monte Carlo tree search methods[END_REF]. UCT is based on generating trajectories by selecting in each MAX node the action that has the highest upper-confidence bound (computed according to the UCB algorithm of Auer et al. [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF]). UCT converges asymptotically to the optimal solution, but its sample complexity can be worst than doubly-exponential in (1/ε) for some MDPs [START_REF] Munos | From bandits to Monte-Carlo tree search: The optimistic principle applied to optimization and planning[END_REF]. One reason for this is that the algorithm can expand very deeply the apparently best branches but may lack sufficient exploration, especially when a narrow optimal path is hidden in a suboptimal branch. As a result, this approach works well in some problems with a specific structure but may be much worse than a uniform sampling in other problems.

On the other hand, a uniform planning approach is safe for all problems. Kearns et al. [START_REF] Kearns | A sparse sampling algorithm for nearoptimal planning in large Markov decision processes[END_REF] generate a sparse look-ahead tree based on expanding all MAX nodes and sampling a finite number of children from AVG nodes up to a fixed depth that depends on the desired accuracy ε. Their sample complexity is2 of the order of (1/ε) log(1/ε) , which is non-polynomial in 1/ε. This bound is better than that for UCT in a worst-case sense. However, as their look-ahead tree is built in a uniform and non-adaptive way, this algorithm fails to benefit from a potentially favorable structure of the MDP.

An improved version of this sparse-sampling algorithm by Walsh et al. [START_REF] Thomas | Integrating sample-based planning and model-based reinforcement learning[END_REF] cuts suboptimal branches in an adaptive way but unfortunately does not come with an improved bound and stays non-polynomial even in the simple Monte Carlo setting for which K = 1.

Although the sample complexity is certainly non-polynomial in the worst case, it can be polynomial in some specific problems. First, for the case of finite N , the sample complexity is polynomial and Szörényi et al. [START_REF] Szörényi | Optimistic planning in Markov decision processes using a generative model[END_REF] show that a uniform sampling algorithm has complexity at most (1/ε) 2+log(KN )/(log(1/γ)) . Notice that the product KN represents the branching factor of the lookahead planning tree. This bound could be improved for problems with specific reward structure or transition smoothness. In order to do this, we need to design non-uniform, adaptive algorithm that captures the possible structure of the MDP when available, while making sure that in the worst case, we do not perform worse than a uniform sampling algorithm.

The case of deterministic dynamics (N = 1) and rewards considered by Hren and Munos [START_REF] Hren | Optimistic Planning of Deterministic Systems[END_REF] has a complexity of order (1/ε) (log κ)/(log(1/γ)) , where κ ∈ [1, K] is the branching factor of the subset of near-optimal nodes. 3 The case of stochastic rewards has been considered by Bubeck and Munos [START_REF] Bubeck | Open-loop optimistic planning[END_REF] but with the difference that the goal was not to approximate the optimal value function but the value of the best open-loop policy which consists in a sequence of actions independent of states. Their sample complexity is (1/ε) max(2,(log κ)/(log 1/γ)) .

In the case of general MDPs, Buşoniu and Munos [START_REF] Buşoniu | Optimistic planning for Markov decision processes[END_REF] consider the case of a fully known model of the MDP. For any state-action, the model returns the expected reward and the set of all next states (assuming N is finite) with their corresponding transition probabilities. In that case, the complexity is (1/ε) log κ/(log(1/γ)) , where κ ∈ [0, KN ] can again be interpreted as a branching factor of the subset of near-optimal nodes. These approaches use the optimism in the face of uncertainty principle whose applications to planning have been have been studied by Munos [START_REF] Munos | From bandits to Monte-Carlo tree search: The optimistic principle applied to optimization and planning[END_REF]. TrailBlazer is different. It is not optimistic by design: To avoid voracious demand for samples it does not balance the upperconfidence bounds of all possible actions. This is crucial for polynomial sample complexity in the infinite case. The whole Section 3 shines many rays of intuitive light on this single and powerful idea.

The work that is most related to ours is StOP by Szörényi et al. [START_REF] Szörényi | Optimistic planning in Markov decision processes using a generative model[END_REF] which considers the planning problem in MDPs with a generative model. Their complexity bound is of the order of (1/ε) 2+log κ/(log(1/γ))+o (1) , where κ ∈ [0, KN ] is a problem-dependent quantity. However, their κ defined as lim ε→0 max(κ 1 , κ 2 ) (in their Theorem 2) is somehow difficult to interpret as a measure of the quantity of near-optimal nodes. Moreover, StOP is not computationally efficient as it requires to identify the optimistic policy which requires computing an upper bound on the value of any possible policy, whose number is exponential in the number of MAX nodes, which itself is exponential in the planning horizon. Although they suggest (in their Appendix F) a computational improvement, this version is not analyzed. Finally, unlike in the present paper, StOP does not consider the case N = ∞ of an unbounded number of states.

Our contributions

Our main result is TrailBlazer, an algorithm with a bound on the number of samples required to return a high-probability ε-approximation of the root node whether the number of next states N is finite or infinite. The bounds use a problem-dependent quantity (κ or d) that measures the quantity of near-optimal nodes. We now summarize the results.

Finite number of next states (N < ∞):

The sample complexity of TrailBlazer is of the order of 4 (1/ε) max(2,log(N κ)/ log(1/γ)+o( 1)) , where κ ∈ [1, K] is related to the branching factor of the set of near-optimal nodes (precisely defined later).

Infinite number of next states (N = ∞):

The complexity of TrailBlazer is (1/ε) 2+d , where d is a measure of the difficulty to identify the near-optimal nodes. Notice that d can be finite even if the planning problem is very challenging. 5 We also state our contributions in specific settings in comparison to previous work.

• For the case N < ∞, we improve over the best-known previous worst-case bound with an exponent (to 1/ε) of max(2, log(N K)/ log(1/γ)) instead of 2 + log(N K)/ log(1/γ) reported by Szörényi et al. [START_REF] Szörényi | Optimistic planning in Markov decision processes using a generative model[END_REF]. • For the case N = ∞, we identify properties of the MDP (when d = 0) under which the sample complexity is of order (in 1/ε 2 ). This is the case when there are non-vanishing actiongaps 6 from any state along near-optimal policies or when the probability of transitionning to nodes with gap ∆ is upper bounded by ∆ 2 . This complexity bound is as good as Monte-Carlo sampling and for this reason TrailBlazer is a natural extension of Monte-Carlo sampling (where all nodes are AVG) to stochastic control problems (where MAX and AVG nodes alternate). Also, no previous algorithm reported a polynomial bound when N = ∞. • In MDPs with deterministic transitions (N = 1) but stochastic rewards our bound is

(1/ε) max(2,log κ/(log 1/γ)) which is similar to the bound achieved by Bubeck and Munos [START_REF] Bubeck | Open-loop optimistic planning[END_REF] in a similar setting (open-loop policies). • In the evaluation case without control (K = 1) TrailBlazer behaves exactly as Monte-Carlo sampling (thus achieves a complexity of 1/ε 2 ), even in the case N = ∞. • Finally TrailBlazer is easy to implement and is numerically efficient. 4 neglecting logarithmic terms in ε and δ 5 since when N = ∞ the actual branching factor of the set of reachable nodes is infinite 6 defined as the difference in values of best and second-best actions 3 2 Monte-Carlo planning with a generative model Objective For any node s, we define the value function V [s] as the optimum over policies π (giving a successor to all MAX nodes) of the sum of discounted expected rewards playing policy π,

1: Input: δ, ε 2: Set: η ← γ 1/ max(2,log(1/ε)) 3: Set: λ ← 2 log(ε(1 -γ)) 2 log log(K) (1-η) log(γ/η) 4: Set: m ← (log(1/δ) + λ)/((1 -γ)
V [s] = sup π E t≥0 γ t r st s 0 = s, π ,
where γ ∈ (0, 1) is the discount factor. If s is an AVG node, V satisfies the following Bellman equation,

V [s] = E [r s ] + γ s ∈C[s] p(s |s)V [s ] . If s is a MAX node, then V [s] = max s ∈C[s] V [s ] .
The planner has access to the oracle which can be called for any AVG node s to either get a reward r or a transition τ which are two independent random variables identically distributed as r s and τ s respectively. With the notation above, our goal is to estimate the value V [s 0 ] of the root node s 0 using the smallest possible number of oracle calls. More precisely, given any δ and ε, we want to output a value µ ε,δ such that P [|µ ε,δ -V [s 0 ]| > ε] ≤ δ using the smallest possible number of oracle calls n ε,δ . The number of calls is the sample complexity of the algorithm.

Blazing the trails with TrailBlazer

To fulfill the above objective, our TrailBlazer constructs a planning tree T which is, at any time, a finite subset of the potentially infinite tree. Only the already visited nodes are in T and explicitly represented in memory. Taking the object-oriented paradigm, each node of T is a persistent object with its own memory which can receive and perform calls respectively from and to other nodes. A node can potentially be called several times (with different parameters) during the run of TrailBlazer and may reuse (some of) its stored (transition and reward) samples. In particular, after node s receives a call from its parent, node s may perform internal computation by calling its own children in order to return a real value to its parent.

Pseudocode of TrailBlazer is in Figure 1 along with the subroutines for MAX nodes in Figure 3 and AVG nodes in Figure 2. A node (MAX or AVG) is called with two parameters m and ε, which represent some requested properties of the returned value: m controls the desired variance and ε the desired maximum bias. We now describe the MAX and AVG node subroutines.

1: Input: m, ε 2: L ← all children of the node 3: ← 1 4: while |L| > 1 and U ≥ (1 -η)ε do 5: U ← 2 1-γ log(K /(δε))+γ/(η-γ)+λ+1 6: for b ∈ L do 7: µ b ← call b with ( , Uη/(1 -η)) 8: end for 9: L ← b : µ b + 2U 1-η ≥ sup j µj -2U 1-η 10: ← + 1 11: end while 12: if |L| > 1 then 13: Output: µ ← max b∈L µ b 14: else { L = {b } } 15: b ← arg max b∈L µ b 16:
µ ← call b with (m, ηε) 17:

Output: µ 18: end if Figure 3: MAX node MAX nodes A MAX node s keeps a lower and an upper bound of its children values which with high probability simultaneously hold at all times. It sequentially calls its children with different parameters in order to get more and more precise estimates of their values. Whenever the upper bound of one child becomes lower than the maximum lower bound, this child is discarded. This process can stop in two ways: 1) The set L of the remaining children shrunk enough such that there is a single child b left. In this case, s calls b with the same parameters that s received and uses the output of b as its own output.

2) The precision we have on the value of the remaining children is high enough. In this case, s returns the highest estimate of the children in L. Note that the MAX node is eliminating actions to identify the best. Any other best-arm identification algorithm for bandits can be adapted instead.

AVG nodes Every AVG node s keeps a list of all the children that it already sampled and a reward estimate r ∈ R. Note that the list may contain the same child multiple times (this is particularly true for N < ∞). After receiving a call with parameters (m, ε), s checks if ε ≥ 1/(1 -γ). If this condition is verified, then it returns zero. If not, s considers the first m sampled children and potentially samples more children from the generative model if needed. For every child s in this list, s calls it with parameters (k, ε/γ), where k is the number of times a transition toward this child was sampled. It returns r + γµ, where µ is the average of all the children estimates.

Anytime algorithm

TrailBlazer is naturally anytime. It can be called with slowly decreasing ε, such that m is always increased only by 1, without having to throw away any previously collected samples. Executing TrailBlazer with ε and then with ε < ε leads to the same amount of computation as immediately running TrailBlazer with ε.

Practical considerations

The parameter λ exists so the behavior depends only on the randomness of oracle calls and the parameters (m, ε) that the node has been called with. This is a desirable property because it opens the possibility to extend the algorithm to more general settings, for instance if we have also MIN nodes. However, for practical purposes, we may set λ = 0 and modify the definition of U in Figure 3 by replacing K with the number of oracle calls made so far globally.

Cogs whirring behind

Before diving into the analysis we explain the ideas behind TrailBlazer and the choices made.

Tree-based algorithm

The number of policies the planner can consider is exponential in the number of states. This leads to two major challenges. First, reducing the problem to multi-arm bandits on the set of the policies would hurt. When a reward is collected from a state, all the policies which could reach that state are affected. Therefore, it is useful to share the information between the policies. The second challenge is computational as it is infeasible to keep all policies in memory.

These two problems immediately vanish with just how TrailBlazer is formulated. Contrary to Szörényi et al. [START_REF] Szörényi | Optimistic planning in Markov decision processes using a generative model[END_REF], we do not represent the policies explicitly or update them simultaneously to share the information, but we store all the information directly in the planning tree we construct. Indeed, by having all the nodes being separate entities that store their own information, we can share information between policies without explicitly having to enforce it.

We steel ourselves for the detailed understanding with the following two arguments. They shed light from two different angles on the very same key point: Do not refine more paths than you need to! Delicate treatment of uncertainty First, we give intuition about the two parameters which measure the requested precision of a call. The output estimate µ of any call with parameters (m, ε) verifies the following property (conditioned on a high-probability event),

∀λ E e λ(µ-V[s]) ≤ exp α + ε|λ| + σ 2 λ 2 2 , with σ 2 = O (1/m) and constant α. (1)
This awfully looks like the definition of µ being uncentered sub-Gaussian, except that instead of λ in the exponential function, there is |λ| and there is a λ-independent constant α. Inequality 1 implies that the absolute value of the bias of the output estimate µ is bounded by ε,

E [µ] -V [s] ≤ ε.
As in the sub-Gaussian case, the second term 1 2 σ 2 λ 2 is a variance term. Therefore, ε controls the maximum bias of µ and 1/m control its sub-variance. In some cases, getting high-variance or low-variance estimate matters less as it is going to be averaged later with other independent estimates by an ancestor AVG node. In this case we prefer to query for high variance rather than a low one, in order to decrease sample complexity.

From σ and ε it is possible to deduce a confidence bounds on |µ -V [s]| by typically summing the bias ε and a term proportional to the standard deviation σ = O (1/ √ m). Previous approaches [START_REF] Szörényi | Optimistic planning in Markov decision processes using a generative model[END_REF][START_REF] Bubeck | Open-loop optimistic planning[END_REF] consider a single parameter, representing the width of this high-probability confidence interval. TrailBlazer is different. In TrailBlazer, the nodes can perform high-variance and low-bias queries but can also query for both low-variance and low-bias. TrailBlazer treats these two types of queries differently. This is the whetstone of TrailBlazer and the reason why it is not optimistic.

Refining few paths

In this part we explain the condition |SampledNodes| > m in Figure 2, which is crucial for our approach and results. First notice, that as long as TrailBlazer encounters only AVG nodes, it behaves just like Monte-Carlo sampling -without the MAX nodes we would be just doing a simple averaging of trajectories. However, when TrailBlazer encounters a MAX node it locally uses more samples around this MAX node, temporally moving away from a Monte-Carlo behavior. This enables TrailBlazer to compute the best action at this MAX node. Nevertheless, once this best action is identified with high probability, the algorithm should behave again like Monte-Carlo sampling. Therefore, TrailBlazer forgets the additional nodes, sampled just because of the MAX node, and only keeps in memory the first m ones. This is done with the following line in Figure 2,

ActiveNodes ← SampledNodes(1 : m).
Again, while additional transitions were useful for some MAX node parents to decide which action to pick, they are discarded once this choice is made. Note that they can become useful again if an ancestor becomes unsure about which action to pick and needs more precision to make a choice. This is an important difference between TrailBlazer and some previous approaches like UCT where all the already sampled transitions are equally refined. This treatment enables us to provide polynomial bounds on the sample complexity for some special cases even in the infinite case (N = ∞).

TrailBlazer is good and cheap -consistency and sample complexity

In this section, we start by our consistency result, stating that TrailBlazer outputs a correct value in a PAC (probably approximately correct) sense. Later, we define a measure of the problem difficulty which we use to state our sample-complexity results. We remark that the following consistency result holds whether the state space is finite or infinite. Theorem 1. For all ε and δ, the output µ ε,δ of TrailBlazer called on the root s 0 with (ε, δ) verifies

P [|µ ε,δ -V [s 0 ]| > ε] < δ.

Definition of the problem difficulty

We now define a measure of problem difficulty that we use to provide our sample complexity guarantees. We define a set of near-optimal nodes such that exploring only this set is enough to compute an optimal policy. Let s be a MAX node of tree T . For any of its descendants s, let c →s (s ) ∈ C [s ] be the child of s in the path between s and s. For any MAX node s, we define

∆ →s (s ) = max x∈C[s ] V [x] -V [c →s (s )] .
∆ →s (s ) is the difference of the sum of discounted rewards stating from s between an agent playing optimally and one playing first the action toward s and then optimally.

Definition 1 (near-optimality). We say that a node s of depth h is near-optimal, if for any even depth h ,

∆ →s (s h ) ≤ 16 γ (h-h )/2 γ(1 -γ)
with s h the ancestor of s of even depth h . Let N h be the set of all near-optimal nodes of depth h.

Remark 1. Notice that the subset of near-optimal nodes contains all required information to get the value of the root. In the case N = ∞, when p(s|s ) = 0 for all s and s , then our definition of near-optimality nodes leads to the smallest subset in a sense we precise in Appendix C. We prove that with probability 1 -δ, TrailBlazer only explores near-optimal nodes. Therefore, the size of the subset of near-optimal nodes directly reflects the sample complexity of TrailBlazer.

In Appendix C, we discuss the negatives of other potential definitions of near-optimality.

Sample complexity in the finite case

We first state our result where the set of the AVG children nodes is finite and bounded by N .

Definition 2. We define κ ∈ [1, K] as the smallest number such that

∃C ∀h, |N 2h | ≤ CN h κ h .
Notice that since the total number of nodes of depth 2h is bounded by (KN ) h , κ is upper-bounded by K, the maximum number of MAX's children. However κ can be as low as 1 in cases when the set of near-optimal nodes is small.

Theorem 2. There exists C > 0 and K such that for all ε > 0 and δ > 0, with probability 1 -δ, the sample-complexity of TrailBlazer (the number of calls to the generative model before the algorithm terminates) is

n(ε, δ) ≤ C(1/ε) max(2, log(N κ) log(1/γ) +o(1)) (log(1/δ) + log(1/ε)) α ,
where α = 5 when log(N κ)/ log(1/γ) ≥ 2 and α = 3 otherwise.

This provides a problem-dependent sample-complexity bound, which already in the worst case (κ = K) improves over the best-known worst-case bound O (1/ε) 2+log(KN )/ log(1/γ) [START_REF] Szörényi | Optimistic planning in Markov decision processes using a generative model[END_REF]. This bound gets better as κ gets smaller and is minimal when κ = 1. This is, for example, the case when the gap (see definition given in Equation 2) at MAX nodes is uniformly lower-bounded by some ∆ > 0.

In this case, this theorem provides a bound of order (1/ε) max(2,log(N )/ log(1/γ)) . However, we will show in Remark 2 that we can further improve this bound to (1/ε) 2 .

Sample complexity in the infinite case

Since the previous bound depends on N , it does not apply to the infinite case with N = ∞. We now provide a sample complexity result in the case N = ∞. However, notice that when N is bounded, then both results apply.

We first define gap ∆(s) for any MAX node s as the difference between the best and second best arm,

∆(s) = V [i ] -max i∈C[s],i =i V [i] with i = arg max i∈C[s] V [i] . (2)
For any even integer h, we define a random variable S h taking values among MAX nodes of depth h, in the following way. First, from every AVG nodes from the root to nodes of depth h, we draw a single transition to one of its children according to the corresponding transition probabilities. This defines a subtree with K h/2 nodes of depth h and we choose S h to be one of them uniformly at random. Furthermore, for any even integer h < h we note S h h the MAX node ancestor of S h of depth h . Definition 3. We define d ≥ 0 as the smallest d such that for all ξ there exists a > 0 for which for all even h > 0,

E     K h/2 1 S h ∈ N h 0≤h <h h ≡0(mod 2) ξ γ h-h 1 ∆(S h h )<16 γ (h-h )/2 γ(1-γ)     ≤ aγ -dh If no such d exists, we set d = ∞.
This definition of d takes into account the size of the near-optimality set (just like κ) but unlike κ it also takes into account the difficulty to identify the near-optimal paths.

Intuitively, the expected number of oracle calls performed by a given AVG node s is proportional to: (1/ε 2 ) × (the product of the inverted squared gaps of the set of MAX nodes in the path from the root to s) × (the probability of reaching s by following a policy which always tries to reach s).

Therefore, a near-optimal path with a larger number of small MAX node gaps can be considered difficult. By assigning a larger weight to difficult nodes, we are able to give a better characterization of the actual complexity of the problem and provide polynomial guarantees on the sample complexity for N = ∞ when d is finite. Theorem 3. If d is finite then there exists C > 0 such that for all ε > 0 and δ > 0, the expected sample complexity of TrailBlazer satisfies

E [n(ε, δ)] ≤ C (log(1/δ) + log(1/ε)) 3 ε 2+d •
Note that this result holds in expectation only, contrary to Theorem 2 which holds in high probability.

We now give an example for which d = 0, followed by a special case of it. Lemma 1. If there exists c > 0 and b > 2 such that for any near-optimal AVG node s,

P [∆ (τ s ) ≤ x] ≤ cx b ,
where the random variable τ s is a successor state from s drawn from the MDP's transition probabilities, then d = 0 and consequently the sample complexity is of order 1/ε 2 . Remark 2. If there exists ∆ min such that for any near-optimal MAX node s, ∆(s) ≥ ∆ min then d = 0 and the sample complexity is of order 1/ε 2 . Indeed, in this case as

P [∆ s ≤ x] ≤ (x/∆ min ) b
for any b > 2 for which d = 0 by Lemma 1.

Conclusion

We provide a new Monte-Carlo planning algorithm TrailBlazer that works for MDPs where the number of next states N can be either finite or infinite. TrailBlazer is easy to implement and is numerically efficient. It comes packaged with a PAC consistency and two problem-dependent sample-complexity guarantees expressed in terms of a measure (defined by κ) of the quantity of near-optimal nodes or a measure (defined by d) of the difficulty to identify the near-optimal paths. The sample complexity of TrailBlazer improves over previous worst-case guarantees. What's more, TrailBlazer exploits MDPs with specific structure by exploring only a fraction of the whole search space when either κ or d is small. In particular, we showed that if the set of near-optimal nodes have non-vanishing action-gaps, then the sample complexity is O(1/ε 2 ), which is the same rate as Monte-Carlo sampling. This is a pretty decent evidence that TrailBlazer is a natural extension of Monte-Carlo sampling to stochastic control problems.

A Consistency

First, we comment on the behavior of the nodes when called multiple times. Remark 3. Notice that the output of the call to a node s with parameters (m, ε) and global δ does not depend on when this call was performed. That means if this call is performed again later in the execution with the same parameters (m, ε) and δ, the output would be the same. This is because, first, the MAX nodes does not store anything so that their behavior only depend on s, m, ε, and δ. Second, the AVG nodes store the already sampled rewards, transitions and the order they were sampled but when called with parameters (m, ε) and δ, it returns the average of the m first sampled children only.

From now we consider any δ, η, and γ, such that 0 < γ < η < 1 and 0 < δ < 1/2.

We define a random variable µ s (m, ε) to be the output of the node s called with parameters (m, ε) and global parameters δ and η. Notice that following Remark 3, µ s (m, ε) does depend on s, m, ε, δ, η and all the transitions and rewards of itself and its descendants but does not depend on the current state of the algorithm.

We define the function U(m, ε) such that the value U(m, ε) + ε is the width of a high-probability confidence interval for µ s (m, ε) which is used in MAX node,

U(m, ε) def = 2 1 -γ log(K /(δε)) + γ/(η -γ) + 1 m • We define max (ε) def = max{ : U( , ε) ≥ (1 -η)ε}.
By definition of TrailBlazer, this is the maximum value that the children of a MAX node would be ever called with. We also define (m, ε)

def = max(m, max (ε))
We now construct event B s for which Lemma 2 proves that it holds with high probability and that on B s , TrailBlazer works well. Definition 4. For any node s, ε > 0, and m ∈ N + ,

• If ε ≥ 1/(1 -γ), B s (m, ε) def = Ω. • Else, if s is a MAX node, B s (m, ε) def = ∩ i∈C[s] ∩ ≤ max(ε) C i, (ε) ∩ B i, max (ε) with    C i, (ε) def = |µ i , η 1-η U ( , ε) -V [i]| ≤ 1 + η 1-η U( , ε) B i, (ε) def = B i , η 1-η U ( , ε) . • Else, if s is an AVG node, B s (m, ε) def = ∩ k∈{1,..,m} B τ s,k (ν s (k, m), ε/γ), with ν s (k, m) def = |{j : τ s,k = τ s,j and j ≤ m}|.
Intuitively, C is the event where the confidence bounds are satisfied. Using events C, Definition 4 recursively constructs the events B needed for Lemma 2.

We remark that there is no loop in this recursive definition as the definition of B s only uses {B i } i where i is a child of s. Moreover, this way we define B s in a finite number of steps as going down through MAX and AVG nodes increases ε to εη/γ and reaches 1/(1 -γ) in a finite number of steps.

Moreover, note that B s is not an intersection of a finite number of C i, because τ s,k is a random variable and therefore in the AVG case, B s is not just the intersection of m events. Remark 4. Notice that, by induction, for any δ > 0, η > 0, ε 1 , ε 2 , m 1 , m 2 , and any node s,

ε 1 ≤ ε 2 and m 1 ≥ m 2 =⇒ B s (m 1 , ε 1 ) ⊆ B s (m 2 , ε 2 ).
Lemma 2. For any k ∈ N + and any node s, integers m ≥ m, and ε > 0, such that δ < (m, ηε) -hs(ε) and

• either k is odd, s is an AVG node, and ε > (γ/η) (k-1)/2

1-γ ,

• or k is even, s is a MAX node, and ε > (γ/η) (k/2)-1

η(1-γ)

, then the two following properties are verified

Property A: P [B s (m, ε)] ≥ 1 -δ (m, ε) K/η 2 hs(ε) -1 K -1 Property B: ∀λ, E e λ(µs(m,ε)-V[s]) B s (m , ε) ≤ exp 2δ (m, ε) η 2 hs(ε) × 1 -εθ(s) η/γ -1 + ε|λ| + λ 2 2(1 -γ) 2 m
,

where θ(s) = 1 -γ if s is a MAX node (1 -γ)/γ if s is an AVG node and h s (ε) = max 0, log(θ(s)γε) log(γ/η)
Lemma 2 claims that on a high-probability event B, the sub-variance of µ s is of order at most 1/m and its bias of order at most ε up to the term independent of λ. We prove it by induction on k.

Proof.

1 • Let k = 1.
Then, k is odd and s is an AVG node. Since ε > 1/(1 -γ), s outputs zero by the definition. For any m ≥ m we have

B s (m, ε) = B s (m , ε) = Ω, since ε > 1/(1 -γ).
Property A: By definition h s (ε) ≥ 0, therefore δ

( (m,ε)K/η 2 ) hs(ε) -1 K-1 ≥ 0. Then, P [Ω] = 1 ≥ 1 -δ (m, ε) K/η 2 hs(ε) -1 K -1 •
Property B: First, as the rewards are bounded from above by one, we have

V [s] ≤ h γ h ≤ 1/(1-γ).
Second, V [s] ≥ 0 as the rewards are bounded from below by zero. Finally, εθ(s) ≤ 1, hence

2δ (m, ε) η 2 hs(ε) × 1 -εθ(s) η/γ -1 ≥ 0.
We combine these bounds to get

∀λ, E e λ(µs(m,ε)-V[s]) B s (m, ε) = E e λ(µs(m,ε)-V[s]) = E e λ(0-V[s]) ≤ exp |λ| 1 -γ ≤ exp (ε|λ|) ≤ exp 2δ (m, ε) η 2 hs(ε) × 1 -εθ(s) η/γ -1 + ε|λ| + λ 2 2(1 -γ) 2 m • 2 • Let k ∈ N + .
Assume that Lemma 2 is true for any k ≤ k -1 and let us prove it for k. We prove it separately for k odd and even.

k is even: Let s be a MAX node, m any integer and ε such that

ε > (γ/η) (k/2)-1 η(1 -γ) and 1 (m, ε) > δ.
We deduce that

ηε > η (γ/η) (k/2)-1 η(1 -γ) = (γ/η) (k/2)-1 1 -γ •
Also, any child of s is an AVG node as nodes alternate between AVG and MAX. We can thus apply the induction assumption to any child of s with parameters (m , ε ) with ε < ηε and m ≤ (m, ε)).

For any child i ∈ C [s] of s and any integer > 0, we remind the reader the definition of B i, ,

B i, = B i , η 1-η U( , ε) . If U( , ε) ≥ (1 -η)ε, then η 1-η U( , ε
) ≥ ηε and ≤ max (ε). We can thus apply our induction assumption to get that for any ≤ max ,

A P [B i, ] ≥ 1 -δ , η 1-η U( , ε) K/η 2 hi η 1-η U( ,ε) -1 K -1 (3) B ∀λ, E e λ(µi( , η 1-η U( ,ε))-V[i]) B i, (4)
≤ exp     2δ   , η 1-η U( , ε) η 2   hi η 1-η U( ,ε) × 1 -η 1-η U( , ε)θ(i) η/γ -1 + η 1-η U( , ε)|λ| + λ 2 2(1 -γ) 2    
We first prove Property A.

Combining Inequality 4 with Lemma 6, we get for all i ∈ C [s] and ≤ max ,

P [C i, c |B i, ] ≤ 2 exp 2δ max (ηε) hi(ηε) × 1 -εθ(i) η/γ -1 exp   - U( , ε) 1 + η 1-η -η 1-η U( , ε) 2 2 (1-γ) 2    ≤ 2 exp 2 1 -εθ(i) η/γ -1 exp - (1 -γ) 2 2 U( , ε) 2 as (δ < 1/ ) ≤ 2 exp 2γ η -γ exp - (1 -γ) 2 2 1 (1 -γ) 2 4 log (K /(δε(1 -γ))) + 4γ/(η -γ) + 4 ≤ 1 2 δε(1 -γ) K 2 ≤ 1 2 δε(1 -γ) K 2 •
Using a union bound over all ≤ max we get

P [∪ ≤ max C i, c |B i, ] ≤ l 1 2 δε(1 -γ) K 2 ≤ π 2 12 δε(1 -γ) K ≤ δε(1 -γ) K •
Now, we can finally bound the event for node s as

P [B s (m, ε)] = P ∩ i∈C[s] (∩ ≤ max C i, ∩ B i, max ) = 1 -P ∪ i∈C[s] (∪ ≤ max C i, c ∪ B i, max c ) ≥ 1 - i∈C[s] P [∪ ≤ max C i, c ∪ B i, max c ] = 1 - i∈C[s] (1 -P [∩ ≤ max C i, ∩ B i, max ]) = 1 - i∈C[s] (1 -P [∩ ≤ max C i, |B i, max ] P [B i, max ]) .
We use Inequality 3 to get

P [B i, max ] ≥ 1 -δ ( max , ηε) K/η 2 hi(ηε) -1 K -1 •
To relate h s (ε) to its child, we realize that i is an AVG node and hence,

h i (ηε) = max 0, log (θ(s)γηε(1 -γ)) log (γ/η) = log (θ(s)γηε(1 -γ)) log (γ/η) = log (ηε(1 -γ)) log (γ/η) = log(η/γ) + log ((1/γ)ε(1 -γ)) log (γ/η) = h s (ε) -1.
We go back to bounding P [B s (m, ε)],

P [B s (m, ε)] ≥ 1 - i∈C[s] 1 - δε(1 -γ) K 1 -δ ( max , ηε) K/η 2 hi(ηε) -1 K -1 ≥ 1 - i∈C[s] δ K + δ ( max , ε) K hs(ε)-1 -1 K -1 (as ε ≤ 1/(1 -γ)) = 1 -δ 1 + K ( max , ε) K hs(ε)-1 -1 K -1 ≥ 1 -δ (m, ε) K/η 2 hs(ε) -1 K -1 •
To prove Property B, we introduce additional notation. We define

µ i, def = µ i , η 1-η U( , ε)
• We also define the empirically best and the best child as

ı def = arg max i∈C[s] µ i, and i def = arg max i∈C[s] V [i] .
We define the random variable L as the set of all children that would be kept in Line 9 in Figure 3 if they are called with parameter as

L def = i ∈ C [s] : µ i, + 2U( , ε) 1 -η ≥ sup j µ j, - 2U( , ε) 1 -η .
Notice that the MAX node is not constructing sets {L } iteratively: The reason is that some of the children might have been eliminated in the earlier rounds. However, at the end of the while loop in Figure 3, we get set L that is the intersection of all {L } ,

L = ∩ ≤ max L .
We define G (which is not a random variable) as

G def = i ∈ C [s] : V [i] + U( , ε) 1 -η ≥ V [i ] - U( , ε) 1 -η • (5)
For simplicity, we also define

ı def = ı max , G def = G max , and for all i ∈ C [s], µ i def = µ i, max .
We can now prove Property B. By definition of C i, and induction assumption, we have

∩ i∈C[s], ≤ max C i, ⊆ {∀i, : |µ i, -V [i]| ≤ U( , ε)/(1 -η)}
and we set

C def = ∩ i∈C[s], ≤ max C i, .
We distinguish two cases depending on the cardinality of the deterministic set G and assume C.

First case: |G| = 1. If |L| > 1 then for all i ∈ L such that i ∈ G = {i } µ i ≤ V [i] + ε (because C i, max holds) < V (i ) -ε (because i / ∈ G s = {i }) ≤ µ s,t i (because C i , max holds)
Therefore arg max i∈L µ i = i which implies that i is called with parameters (m , ε ) with ε ≤ ηε ≤ ε and m ≥ m.

If |L| = 1, then L = {i } and consequently i is called with parameters (m, ε).

In both cases, on event C, µ s (m, ε) is the output of a call to i with parameters (m , ε ) with m ≥ m and ε ≤ ε and m ≤ max . We use Property B of induction assumption as max ≥ m and

U( max , ε)η/(1 -η) ≤ ε. ∀λ, E e λ(µ i (m ,ε )-V[i ]) B i , max = E e B i , max ≤ exp   2δ (m , ε ) η 2 h i (ε ) × 1 -ε θ(i ) η/γ -1 + ε |λ| + λ 2 2(1 -γ) 2 m   ≤ exp 2δ (m, ηε) η 2 hs(ηε) × 1 -ηεθ(s) η/γ -1 + ε|λ| + λ 2 2(1 -γ) 2 m
Also, because for all i, j ∈ C [s] with i = j, B i and µ j are independent, we can claim that

∀λ, E e λ(µs(m,ε)-V[s]) ∩ i∈C[s] B i, max ≤ exp 2δ (m, ηε) η 2 hs(ηε) × 1 -ηεθ(s) η/γ -1 + ε|λ| + λ 2 2(1 -γ) 2 m
• Using Lemma 7 we finaly get for all λ

E e λ(µs(m,ε)-V[s]) B s (m, ε) = E e λ(µs(m,ε)-V[s]) ∩ i∈C[s] (∩ ≤ max C i, ∩ B i, max ) ≤ E e λ(µs(m,ε)-V[s]) ∩ i∈C[s] B i, max P ∩ i∈C[s], ≤ max C i, ∩ i∈C[s] B i, max ≤ exp 2δ (m,ηε) η 2 hs(ηε) × 1-ηεθ(s) η/γ-1 + ε|λ| + λ 2 2(1 -γ) 2 m 1 -δε(1 -γ)
To simplify the expression, we note that for 0 < x ≤ 1/2, we have that 1/(1-x) ≤ 1+2x ≤ exp(2x). We also use that δ < 1/2. exp 2δ (m,ηε)

η 2 hs(ηε) × 1-ηεθ(s) η/γ-1 1 -δε(1 -γ) ≤ exp 2δ (m, ηε) η 2 hs(ηε) × 1 -ηεθ(s) η/γ -1 + 2δε(1 -γ) ≤ exp 2δ (m, ε) η 2 hs(ε) 1 -ηε(1 -γ)/γ η/γ -1 + ε(1 -γ) = exp 2δ (m, ε) η 2 hs(ε) 1 -ε(1 -γ) η/γ -1 = exp 2δ (m, ε) η 2 hs(ε) 1 -εθ(s) η/γ -1 Second case: |G| > 1. For all i ∈ G l , µ i, ≥ V [i] -U( , ε)/(1 -η) (because C holds) ≥ V (i ) -3U( , ε)/(1 -η) (because i ∈ G l ) ≥ µ s,t i -4U( , ε)/(1 -η) (because C holds)
Therefore for all ≤ max , G ⊂ L .

As a result |L| > |G| > 1. The output is then the maximum of the estimates in L. The best estimate µ ı in L verifies

µ ı ≥ µ i (by definition of ı) ≥ V (i ) -U( max , ε)/(1 -η) (because C holds).
It also verifies that

µ ı ≤ V ( ı) + U( max , ε)/(1 -η) (because C ı, max holds) ≤ V (i ) + U( max , ε)/(1 -η)
(by definition of i ).

Since µ s (m, ε) = µ ı , we have

µ s (m, ε) -V [s] ≤ U( max , ε) 1 -η ≤ ε. It follows that E e λ(µs(m,ε)-V[s]) |B s (m, ε) ≤ exp (|λ|ε) ≤ exp 2δ (m, ε, ) η 2 hs(ε,) × 1 -ε, θ(s) η/γ -1 + |λ|ε + λ 2 2(1 -γ) 2 m
• k is odd: Let s be a AVG node, m ∈ N + 0 , and ε such that

ε > (γ/η) (k-1)/2 (1 -γ) and 1 max(m, max (ε)) > δ.
Hence, we can deduce

ε γ > (γ/η) (k-1)/2 γ(1 -γ) = (γ/η) k-1 2 -1 η(1 -γ) •
Also, any child of s is a MAX node as nodes alternate between AVG and MAX. We can thus apply induction assumption to any child of s with parameters (m , ε ) with ε ≥ ε/γ and m ≤ max(m, max (ε)).

For any child i ∈ C [s] of s and any integer ≤ m we can thus apply our induction assumption to get

P [B i (l, ε/γ)] ≥ 1 -δ (l, ε/γ) K/η 2 hi(ε/γ) -1 K -1 ∀λ, E e λ(µi( ,ε/γ))-V[i]) B i ( , ε/γ) ≤ exp 2δ ( , ε/γ) η 2 hi(ε/γ) × 1 -ε/γθ(i) η/γ -1 + |λ|ε/γ + λ 2 2(1 -γ) 2
We state again the definition of B s

B s (m, ε) def = ∩ i∈C[s] B i (ν i (m), ε/γ) with ν i (m) def = |{j : τ s,j = i and j ≤ m}|.
We first prove property A.

P ∪ i∈C[s] B i (ν i (m), ε/γ) c ≤ i∈C[s] P [B i (ν i (m), ε/γ) c ] = k i∈C[s] P [B i (ν i (m), ε/γ) c |ν i (m) = k] dπ(k) = k i∈C[s] P [B i (k i , ε/γ) c |ν i (m) = k i ] dπ(k) = k i∈C[s] P [B i (k i , ε/γ) c ] dπ(k) ≤ k i∈C[s];ki>0 δ (k i , ε/γ) K/η 2 hi(ε/γ) -1 K -1 π(k) ≤ δ (k i , ε) K/η 2 h(ε) -1 K -1 π(k) = δ (k i , ε) K/η 2 h(ε) -1 K -1
We now prove property B.

E   exp   (λ/m) i∈C[s] ν i (m)µ i (ν i (m), ε/γ)   B s (m, ε)   = k E   exp   (λ/m) i∈C[s] ν i (m)µ i (ν i (m), ε/γ)   B s (m, ε), (ν i (m)) = (k i )   dπ(k) = k E   exp   (λ/m) i∈C[s] k i µ i (k i , ε/γ)   B s (m, ε)   dπ(k) = k E   i∈C[s] exp ((λ/m)k i µ i (k i , ε/γ)) B s (m, ε)   dπ(k) ≤ k exp   i∈C[s] k i 2δ (k i , ε/γ) η 2 hi(ε/γ) × 1 -ε/γθ(i) η/γ -1 + |(λ/m)|ε/γ + (λ/m) 2 2(1 -γ) 2 k i   dπ(k) = k exp m2δ (m, ε/γ) η 2 hi(ε/γ) × 1 -ε/γθ(i) η/γ -1 + |λ|ε/γ + λ 2 2(1 -γ) 2 m dπ(k) ≤ exp 2δ (m, ε) η 2 hs(ε) × 1 -εθ(s) η/γ -1 + |λ|ε/γ + λ 2 2(1 -γ) 2 m
We can finally bound

E   exp   λr + γ(λ/m) i∈C[s] ν i (m)µ i (ν i (m), ε/γ)   B s (m, ε)   ≤ exp λ 2 2m + 2δ (m, ε) η 2 hs(ε) × 1 -εθ(s) η/γ -1 + |λ|ε + λ 2 γ 2 2(1 -γ) 2 m ≤ exp 2δ (m, ε) η 2 hs(ε) × 1 -εθ(s) η/γ -1 + |λ|ε + λγ √ 2m(1 -γ) + λ √ 2m 2 = exp 2δ (m, ε) η 2 hs(ε) × 1 -εθ(s) η/γ -1 + |λ|ε + λ 2 2(1 -γ) 2 m •
Now we stat and prove Theorem 4 that is equivalent to Theorem 1. Theorem 4. Let µ be the output of TrailBlazer.For all δ > 0, ε > 0, there exists an event B such that

P [B] ≥ 1 -δ and on B, |µ -V [s 0 ]| ≤ ε.
Proof. We apply Lemma 2 to the root node s 0 . Let us note the following:

m = (log(1/δ) + λ)/((1 -γ) 2 ε 2 ) λ = 2 log(ε(1 -γ)) 2 log log(K) (1-η) log(γ/η)
We set δ = exp(λ)δ . We apply Lemma 2 on s 0 with (m, ε/2). We take

B = B s0 (log(1/δ) + λ)/((1 -γ) 2 ε 2 ).
We check that

δ (m, ε) K/η 2 hs(ε) -1 K -1 ≤ δ
and apply Lemma 6 with property B of Lemma 2 to get

|µ s0 (m, ε/2) -V [s 0 ]| ≤ ε.

B Sample complexity

Definition 5 (node classification). A node s of depth h is said to be

• suboptimal if there exists h < h, such that ∆ →s (s h ) > 16 γ (h-h )/2
γ(1-γ) ,

• optimal if for all h < h, ∆ →s (s h ) = 0,

• near-optimal in all other cases.

Lemma 3. For any node s of depth h, integer m, and ε > 0, a call to s with (m, ε) does not generate call to nodes of depth h > h + h max (ε) with If B s0 (m, ε) holds then for all nodes s of depth h and h < h, the following holds,

h max (ε) ≤ 2 log (ε(1 -γ)) log (γ/
s is suboptimal =⇒ Γ s = 0 ∆(s h ) > 16 γ (h-h )/2 γ(1 -γ)
=⇒ any call to s was generated by a call from s h with m ≤ Γ s h .

Proof. For any node s of depth h, we assume that Γ s > 0 and we want to prove that s is not suboptimal. Let h be some integer with h < h. We note (m s h , ε s h ) the parameters s h has been called with. If Γ s > 0, then, by Lemma 3, that means that s h performs a call to its child i s in direction to s with parameters (m, ε) such that h max (ε) ≥ h -h .

For any m , ε , if s h performs a call to i s with parameters (m , ε ) with ε ≤ ε then 1. either i s ∈ L max(ms h ,εs h ) and ηε s h ≤ ε, 2. or i s ∈ L for some such that η 1-η U ≤ ε.

We use the definition of G (5) and the result of Equation 6.

Case 1. For all we have L ⊂ G therefore

V [i ] -V [i s ] ≤ 4U 1-η ≤ 4 ε η • Case 2. By definition of max , U max+1 ≤ (1 -η)ε s h . We also have that ∀ ∃a ≥ 0 U +1 U = + 1 log( + 1) + a log( ) + a ≥ + 1 ≥ 1/ √ 2.
We put it together using i s ∈ L max ⊂ H max ,

V [i ] -V [i s ] ≤ 4U max 1 -η ≤ U max ≤ 4 √ 2 1 -η U max+1 ≤ 4 √ 2ε s ≤ 4 √ 2ε/η. Either V [i ] -V [i s ] ≤ 4ε η or V [i ] -V [i s ] ≤ 4ε √ 2 η
from which we deduce that

V [i ] -V [i s ] ≤ 4ε √ 2 η • We use h -h ≤ h max (ε) ≤ 2 log(ε(1-γ)) log(γ/η) + 1 to get (γ/η) (h-h )/2 ≥ ε(1 -γ)γ/η and using η ≤ 1 ε ≤ 1 γ(1 -γ) (1/η) (h-h )/2 γ (h-h )/2 . As h -h ≤ h max (ε 0 ), log (1/η) (h-h )/2 ≤ log (1/η) hmax(ε0)/2 ≤ log 1/γ log(1/ε) log(1/ε) log (η/γ) ≤ 1.
We combine

ε ≤ e γ(1 -γ) γ (h-h )/2 and V [i ] -V [i s ] ≤ 4ε √ 2 η , to finally get ∆ →s (s h ) = V [i ] -V [i s ] ≤ 4e √ 2 γ(1 -γ) γ (h-h )/2 ≤ 16 γ (h-h )/2 γ(1 -γ) .
This is true for any h < h which proves that s is not suboptimal.

We now assume that ∆(s h ) > 16 γ (h-h )/2 γ(1-γ) and prove that if a call to i later generated calls to s, then this call to i was done with (m, ε) with m ≤ Γ s h . By Lemma 3, h max (ε) ≥ h -h . We distinguish three cases, 1. either s h calls i within the while loop, 2. or s h calls i outside the while loop with |L| > 1, 3. or s h calls i outside the while loop with |L| = 1.

Case 1. If this call is made to i = i with ε = η 1-η U , with the same reasoning as above we conclude that V

[i ] -V [i s ] ≤ 4e √ 2
γ(1-γ) γ h-h which contradicts the assumption. Case 2. We use again the same reasoning as above with U max to conclude that if ε s h ≥ ηε then

V [i ] -V [i s ] ≤ 4e √ 2 γ(1-γ) γ h-h . Case 3.
In this last case, the call is done with (m s h , ε s h ) and m s h ≤ Γ s h .

B.1 Finite case: Proof of Theorem 2

Proof. Any near-optimal node of depth h called with parameters (m , ε ) generates at most

log(1/(εδ)) (1 -η) 2 (1 -γ) 2 ε 2 samples.
This follows from the same reasoning as in previous theorem. Then we can use that ε ≤ ε(η/γ) h/2 . Therefore, the maximum number of samples T is

T ≤ h≤hmax N h κ h log(1/(εδ)) (1 -η) 2 (1 -γ) 2 ε 2 η h/2 /γ h/2 • If N κ/γ 2 > 1, then, T ≤ O h max (N κη/γ) hmax log(1/(εδ)) (1 -η) 2 (1 -γ) 2 ε 2 = O h max (1/ε) (log(N κ)+2 log(η/γ))/ log(η/γ) log(1/(εδ)) (1 -η) 2 (1 -γ) 2 ε 2 = O log(1/ε) log(η/γ) (1/ε) log(N κ)/ log(η/γ) log(1/(εδ)) (1 -η) 2 (1 -γ) 2 = O (1/ε) log(N κ)/ log(1/γ)+o(1) log(1/(εδ)) log(1/ε) 4 (1 -γ) 5 Setting η = γ 1/ max(2,log(1/ε)) = O (1/ε) log(N κ)/ log(1/γ)+o(1) log(1/δ) + log(1/ε) 1 -γ 5 If N κ/γ 2 < 1, then, T ≤ O (1/ε) 2 1 1 -N κ/γ 2 log(1/(εδ)) (1 -η) 2 (1 -γ) 2 = O (1/ε) 2 (log(1/δ) + log(1/ε)) 3 (1 -N κ/γ 2 )(1 -γ) 4 • Setting η = γ 1/ max(2,log(1/ε)) If N κ/γ 2 = 1, then, T ≤ O h max (1/η) hmax log(1/(εδ)) (1 -η) 2 (1 -γ) 2 ε 2 T ≤ O log(1/(εδ)) log(1/ε) 4 (1 -γ) 5 ε 2 (1/η) hmax = O(1) T ≤ O (1/ε) 2 log(1/δ) + log(1/ε) 1 -γ 5 • B.

Infinite case

For any MAX node s of depth h and any h > h , we define Γ h s as the maximum m the node s has been called with; such that this call generated calls up to depth h. Lemma 5. There exists ξ such that for any node s of depth h and any ancestor s h of s of depth h ≤ h,

E [Γ s ] ≤ 1 (s is not suboptimal) p(s|s 0 ) 0≤h <h h ≡0(mod 2) 1 ξγ h-h 1 ∆(s h )<16 γ (h-h )/2 γ(1-γ) Γ s0 • Proof. We prove that E Γ h s h +1 ≤ 1 (s h +1 is not suboptimal) E Γ s h h p(s|s ) 1 ξγ h-h 1 ∆(s h )<16 γ (h-h )/2 γ(1-γ) • If ∆ s h ≥ 16 γ (h-h )/2
γ(1-γ) then by Lemma 4, either s is suboptimal and therefore Γ s = 0 or the calls of s h towards s which generate calls up to depth h are done with parameters m ≤ E Γ h s .

In the other case, a node called with (m, ε) generates at most O 1/ε 2 samples. Moreover, since s h is of depth h , we know that it is called with ε ≥ O (γ/η) (h-h )/2 ≥ ξγ (h-h )/2 /K for some ξ as

(1/η) hmax = O(1).
Also, if no sample from s h reaches s h +1 then Γ s = 0. Together, we get that

Γ s h +1 ≤ 1   One of the 1 ξγ h-h 1 ∆(s h )<16 γ (h-h )/2 γ(1-γ) samples from s h goes to s h +1   Γ s h E Γ s h +1 ≤ P   One of the 1 ξγ h-h 1 ∆(s h )<16 γ (h-h )/2 γ(1-γ) samples from s h goes to s h +1   E Γ s h ≤ E Γ s h h p(s|s ) 1 ξγ h-h 1 ∆(s h )<16 γ (h-h )/2 γ(1-γ) • Now E Γ s h +1 = 0 and therefore s h +1 is h-suboptimal. Moreover, if one of the s h is h-sub-
optimal then s is suboptimal. To finish the proof, we multiply all these inequalities over h .

B.2.1 Proof of Theorem 3

Proof. We bound the expected number of oracle calls as

E [n(ε, δ)] ≤ E   hmax(ε) h=1 s of depth h Γ s   = hmax(ε) h=1 s of depth h E [Γ s ] ≤ hmax(ε) h=1 K h/2 E     1 S h ∈ N h 0≤h <h h ≡0(mod 2) 1 ξγ 2(h-h ) 1 ∆(S h h )<16 γ (h-h )/2 γ(1-γ)     Γ s0 ≤ hmax(ε) h=1 aγ -dh Γ s0 ≤ ah max (ε) O (1/ε) d Γ s0 ≤ O log(1/ε)(1/ε) d (log(1/δ) + log(1/ε)) 2 ε 2 ≤ O (log(1/δ) + log(1/ε)) 3 ε 2+d •
This is true on B s (m, ε) which holds with probability 1 -δ but as δ O 1/ε 3 , the result also holds on the whole probability space.

B.2.2 Proof of Lemma 1

Proof. We note

x h def = 1 ξγ h-h and y h def = 16 γ (h-h )/2 γ(1-γ) so that E     1 S h ∈ N h 0≤h <h h ≡0(mod 2) 1 ξγ h-h 1 ∆(S h h )<16 γ (h-h )/2 γ(1-γ)     = E     1 S h ∈ N h 0≤h <h h ≡0(mod 2) x 1(∆(S h h )<y h ) h     .
Since at all AVG nodes transitions are taken independently:

E     1 S h ∈ N h 0≤h <h h ≡0(mod 2)
ξ 2 γ (b-2)h   ≤ P S h ∈ N h exp ∞ h =0 ξ 2 γ (b-2)h ≤ P S h ∈ N h ξ 3
We then conclude with

P S h ∈ N h ≤

C On the choice of the near-optimality definition

In this part, we discuss the possible alternative choices to the near-optimality set N h and explain the choice we made in Definition 1. From this definition, the knowledge of the rewards of the near-optimal nodes is enough to compute an optimal policy. Therefore, whatever the rewards associated with the non-near-optimal nodes are, they do not change the optimal policy. A good adaptive algorithm, would not exploring the whole tree, but would with high probability only explore a set not significantly larger than a set N h of near-optimal nodes. There are other possible definitions of the near-optimality set satisfying this desired property that we could possibly consider for the definition of N h . Ideally, this set would be as small as possible.

Alternative definition 1 An idea could be to consider for any node s and an ancestor s , the value of the policy starting from s choosing the action leading to s every time it can and playing optimally otherwise. Let V s (s ) be this value. A natural approach would be to define N h only based on V s (s 0 ) with s 0 being the root. To ensure that exploring the near-optimal nodes is enough to compute the optimal policy we need any near-optimal nodes s of depth h to verify

(V [s 0 ] -V s (s 0 )) p(s|s 0 ) ≤ γ h 1 -γ •
When there is no AVG node, p(s) = 1 and this definition coincides with the near-optimal set defined in OLOP [START_REF] Bubeck | Open-loop optimistic planning[END_REF]. Nevertheless, in our case, p(s) can be low and even 0 in the infinite case. When p(s|s 0 ) = 0 for all s, any node would be near-optimal which is bad because then the near-optimality set is large and the sample complexity of the resulting algorithm suffers.

Alternative definition 2 There is a smarter way to define a near-optimality set. We could define it as all the nodes s of depth h such that ∀h ≤ h,

h i=h p(s i |s h )γ i (V [s i ] -V [c(s i , s)]) ≤ γ h 1 -γ •
Taking h = 0, we recover the alternative definition 1. Alternative definition 2 defines a smaller nearoptimality set, as the condition needs to be verified for all h and not just for h = 0. This definition would also lead to a smaller near-optimality set than our Definition 1. Indeed, in Definition 1 we consider only the first term of the sum above. This enables the algorithm to compute its strategy from a node s only based on the value of its descendants. This is an ingredient leading to better computational efficiency of TrailBlazer. Indeed, the computational complexity of our algorithm is linear in the number of calls to the generative model. In addition, when the transition probabilities are low, both definitions are close and coincide in the infinite case, when the probability transitions are continuous and therefore p(s|s ) = 0. In this case (p(s|s ) = 0), for any set of nodes of not null probability mass, it is possible to change the value of the rewards on this subset such that the value of the root change. In this sense, this definition of near-optimal nodes is minimal when p(s|s ) = 0.

D Auxiliary material

Figure 2 :

 2 Figure 2: AVG node

  h ) h = P S h ∈ N h 0≤h <h h ≡0(mod 2) 1 + x h P ∆(S h h ) < y h ≤ P S h ∈ N h 0≤h <h h ≡0(mod 2) 1 + ax h y b h ≤ P S h ∈ N h exp   0≤h <h,h ≡0(mod 2) ax h y b h   ≤ P S h ∈ N h exp   0≤h <h,h ≡0(mod 2)

  Ka 16 γ (h-h )/2 γ(1 -γ) b ≤ ξ 4 /K h/2 .

Lemma 6 .β 2 λ 2 2 , 2 •β 2 λ 2 2 ≤ 2 • 2 • 2 • 7 .E

 62222227 If for some α > 0, β, M , a random variable X verifies that E e λX ≤ M exp |λ|α + then for all u > 0P [|X| ≥ u] ≤ 2M exp -(u -α) 22βProof. Using Markov inequality for u > α and setting λ = u-α β 2 , we getP [X ≥ u] ≤ E e λX e λu ≤ M exp |λ|α -λu + M exp -(u -α) 22βSimilarly, for u < -α,P [X ≤ u] ≤ M exp -(u + α) 2 2βTaking a union bound, we have that for all u > α,P [|X| ≥ u] ≤ 2M exp -(u -α) 22βLemma For any positive random variable X and event A,[X] = E [X|A] P [A] + E [X|A c ] P [A c ] ≥ E [X|A] P [A]

  η) •MAX and AVG nodes alternate. First, Any AVG node called with ε can only perform calls with ε ≥ ε/γ. Second, MAX nodes called with ε can only perform calls with ε ≥ ηε. Finally, an AVG node called with ε = 1/(1 -γ) cannot perform a call to deeper node. For any node s called with different (m t , ε t ) during the algorithm execution, we let Γ s = max t m t be the maximum m t that s has been called with. Lemma 4. We denote the parameters the root has been called with by (m 0 , ε 0 ).

neglecting exponential dependence in γ

nodes that need to be considered in order to return a near-optimal approximation of the value at the root

Acknowledgements The research presented in this paper was supported by French Ministry of Higher Education and Research, Nord-Pas-de-Calais Regional Council, a doctoral grant of École Normale Supérieure in Paris, Inria and Carnegie Mellon University associated-team project EduBand, and French National Research Agency projects ExTra-Learn (n.ANR-14-CE24-0010-01) and BoB (n.ANR-16-CE23-0003)


Keywords: Antennas, antenna surrogate model, statistical antenna modeling, spherical harmonics, polynomial chaos

A new methodology of statistical modeling of the far field (FF) radiated by antennas undergoing random disturbances is presented. Firstly, the radiated FF is transformed into a parsimonious form using the Spherical Modes Expansion Method (SMEM); then a surrogate model relating the parsimonious field with the input random parameters is constructed using the Polynomial Chaos Expansion Method (PCEM). The combination of the SMEM and the PCEM allows to develop a compact and precise model with a minimized experimental design cost. The obtained model is computationally costless for generating statistical samples of disturbed antennas easily usable as surrogate models in various types of analyses. In order to demonstrate its performance, the proposed methodology is validated with a deformable canonical antenna -a dipole undergoing three independent random deformations (stretching, bending and torsion), deriving a compact and precise surrogate model.

I. INTRODUCTION

OWADAYS antennas are increasingly used in complex and fluctuating environment where random disturbances from the nearby surroundings can have non-negligible impact on antennas performance. Besides, antennas may be also subject to other strong randomness due to their intrinsic variability. This is the case for wearable, deformable or reconfigurable antennas [START_REF] Rais | A review of wearable antenna[END_REF], which are rapidly emerging in the communication and medical fields. Therefore, quantitatively characterize the performance of antennas undergoing random disturbances is of significant importance.

In the past few years, efforts have been made to interpret random effects on antenna global performance indicators (GPIs) such as the total efficiency, bandwidth, resonance frequency, radiation patterns [START_REF] Mellah | Statistical Modeling of the "Antenna-Head" Interaction[END_REF] [START_REF] Boeykens | Cylindrical bending of Deformable textile rectangular patch antennas[END_REF], etc., or on diffuse channels [START_REF] Kersaudy | Stochastic analysis of scattered field by building facades using polynomial chaos[END_REF].

To our knowledge, it is the first time a methodology is proposed for modeling statistically the FF radiated by antennas in uncertain configurations or environments. Such overall statistical modeling provides a complete and efficient substitute for simulated (with EM solvers) or real antennas: any FF antenna response in the considered fluctuating conditions can be computed with a negligible computational cost. This kind of surrogate model can be useful for various types of analyses beyond This paper is submitted on April 22, 2016 for review.

J. Du and C. Roblin are with the Department of Communication and Electronics, Télécom ParisTech (Institut Mines-Télécom and LTCI), Paris 75634, France (e-mail: jinxin.du@telecom-paristech.fr). GPI-related ones, such as combined modeling approaches, e.g. hybrid statistical antenna and directional channel models analyses [START_REF] Sibille | A Channel-Based Statistical Approach to Antenna Performance in UWB Communications[END_REF] (joint antenna/channel analyses can be driven with a low computational cost) or complete antenna models plugged in asymptotic simulators (e.g. ray tracing), etc., notably when the knowledge of both polarization and phase is mandatory (e.g. for MIMO systems). The only real drawback of the proposed method compared to the statistical modeling of specific GPIs is an increased complexity of the modeling.

In this letter, we present a methodology for the overall modeling of "variable antennas", in which the procedure is split into two steps. The first one is to achieve a parsimonious but complete representation of the radiated FF, and the second consist in developing a statistical model which can quantify the effects of antennas random variability on their radiation properties. Several techniques such as the Spherical Mode Expansion Method (SMEM) [START_REF] Hansen | Spherical near-field antenna measurements[END_REF], the Singularity expansion method [START_REF] Baum | On the singularity expansion method for the solution of Electromagnetic interaction problems[END_REF], [START_REF] Ch | Ultra compressed parametric modelling of UWB antenna measurements[END_REF], and the Fourier series expansion [START_REF] Mhanna | Parametric statistical modeling of power gain patterns for RFID backscattering channels[END_REF] can be used to efficiently represent the field. The underlying common principle is to represent the EMF in a more suitable space in order to achieve a considerable compression of information, prior to consider randomness effects with an achievable statistical approach. The SMEM was chosen for its efficiency as regards parsimony, in particular for small to moderate size antennas (for terminals, sensors, RFID tags, etc.). The Polynomial Chaos Expansion Method (PCEM) was chosen for its efficiency to model stochastic multivariate input/output problems, and its ability to deduce compact and accurate surrogate models with considerably lower Experimental Design (ED) cost than classical methods such as e.g. Monte Carlo. Once the model is available, one can easily generate antennas corresponding to any values of the random input parameters.

II. PROPOSED METHODOLOGY

A. Parsimony of the EMF with the SMEM

According to [START_REF] Ch | Characterization tools of antennasin the Time Domain[END_REF], the electromagnetic far field of any passive antenna can be fully represented by its Antenna Transfer Function (ATF) H (f, r ). Based on the Spherical Modes Expansion theory, the ATF can be expanded as:

1 1 0 4 ˆˆ( , ) ( , ) ( ) ( ) ( ) jkr p p p re f f H f a f        r E r r  H ∞   (1)
where p = 2{n(n + 1) + m 1} + s is the condensed global index related to the modality s (= 1 or 2), degree n, and order m [START_REF] Hansen | Spherical near-field antenna measurements[END_REF]. In the spherical coordinate system, the Vector Spherical Harmonics (VSHs)  p (r ) are orthonormal bases in the space of vector spherical functions and H p (f ) are complex coefficients which are obtained by projection H p (f ) H (f, r ),  p * (r ) .

Statistical Modeling of Disturbed Antennas Based on the Polynomial Chaos Expansion

Jinxin Du, Christophe Roblin, Member, IEEE

N

In practice the SME is truncated to a maximal degree Nmax (i.e. to Pmax = 2Nmax (Nmax + 2)) corresponding to the most energetically dominant VSH modes (N max ≈ kr 0 +  ( kr 0 3), where r 0 is the radius of the antenna minimal sphere [START_REF] Jensen | On the number of modes in spherical wave expansion[END_REF]), which is rather low for small/moderate size antennas. Moreover, as the real and imaginary parts of VSH modes have independent contributions to the total energy, we can keep only most dominant real and/or imaginary parts of H p (f ), as in (2):

mod( , 2) /2 ˆˆ( , ) ( ) ( ) ( ) d q q q q I f j j f         r r  H   (2)
where

q = 2(p  1) + u with u  {1, 2},  q (f ) =Re H p (f ) or  q (f ) =Im H p (f ) if u = 1 or 2; "mod(a, n)" means "a modulo n" and "   
  " is the ceiling round function; Id is the index set of the selected most dominant  q (f ) such that the collected energy reaches a predefined relative threshold Ec_th (e.g. 99 %) compared to the total energy corresponding to degree Nmax. Note that the cardinal of Id is expected to be much smaller than 2Pmax, achieving a sparser field representation which is beneficial for improving the model's compactness. By reversing the SME transform we are able to reconstruct the ATF from the set of selected dominant  q (f ). To evaluate the accuracy of the truncated expansion, an error indicator named Relative Root Mean Square Error of the Realized Gain (RRMSE of Gr) is defined in (

ˆ[ ( , )

/

m m r init rcs r r m init r m RRMSE G f G f G f d G f d           r r r (3)
where G r init and G r rcs denote respectively the realized gain of the initial and reconstructed fields, d is the infinitesimal solid angle,

 m = {(θ, φ) | G r (θ, φ) ≥ Max θ, φ
G r θ, φ 10 dB} is the main lobe(s) solid angle and reminding that G r (f, r )=ǁH (f, r )ǁ 2 [START_REF] Ch | Characterization tools of antennasin the Time Domain[END_REF]. This indicator is also used to assess the field variations due to disturbances, G r rcs and G r init being replaced by the Gr of the fields with and without disturbances.

B. Polynomial Chaos Expansions (PCE) of VSH coefficients

The PCE is an efficient method for describing the uncertainty propagation characteristics of an uncertain system. The input/output relation representing such system is denoted Y = (X) , where X = (X 1 , X 2 … X M ) is the random vector (spanning the input stochastic space) with known Probability Density Function (PDF) f X , and Y is a scalar output which can be expanded as the weighted sum of polynomials (4) [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF]:

= ( )= ( ) M Y y      X X    (4)
where α is a multi-index, Φ α (X) form a multivariate polynomial basis in the Hilbert space

L P X 2 (R M , R) with the inner product <u, v> = u(x)v(x)f X (x)dx R M
, and y α are the PCE coefficients. The PCE may also have sparse expression where only desired polynomials are conserved, as in ( 5):

( ) Y y      X    (5)
where  ⊂ N M is a subset of index of "dominant polynomials" of total degree less than a prescribed value.

The prediction accuracy of the PCEM will be evaluated and controlled by the cross-validation indicator known as the Leave-One-Out error (e LOO ) [START_REF] Molinaro | Prediction error estimation: A comparison of resampling methods[END_REF] defined in (6):

( ) \ ( ) 2 ( ) 2 1 1 [ ( ) ( )] / ( ) N N n P C n n n LOO Y n n e Y           x x ( 6
)
where  = {x (1) ,…, x (N) } T is the Experimental Design (ED) of N realizations of X throughout which the system is pre-evaluated, with responses noted as  = {Y (1) ,…, Y (N) } T ;  PC\n is the PCE model derived from the whole set ,  except {x (n) , Y (n) }, and μ Y is the mean of the responses.

Taking the random disturbances as inputs, and the dominant  q (f ) as outputs, we can expand each  q (f ) into the form of ( 5) to obtain a compact surrogate model [START_REF] Hansen | Spherical near-field antenna measurements[END_REF] for the ATFs:

mod( , 2) / 2 ˆˆ( , ) ( ) ( ) ( ) d q q q q q I j j y            r r H X X      (7)
Note that the parameter "f " is omitted in [START_REF] Hansen | Spherical near-field antenna measurements[END_REF] because it can be either set at a fixed value or considered as a component of X.

There are a variety of methods for estimating the PCE coefficients. Without trying to be exhaustive, we can mention the Monte-Carlo (MC), the Quadrature Projection (QP), the Ordinary Least Square Error (OLSE), and the Hybrid Least Angle Regression (HLARS) methods [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF], etc. Among all these methods, the HLARS has the advantage of being more compact as it conserves only polynomials of the most "correlation-dependent" importance. Moreover, the HLARS is also computationally efficient as thanks to its compact expression, the ED can be designed much tighter, reminding that the cost for the antennas pre-evaluation is proportional to the cardinal of the ED (NED = card(ED)). Additionally, the HLARS can be used in an adaptive way in which the PCE is applied iteratively for increasing polynomial degrees until the expansion accuracy cannot be further improved. In this letter, the adaptive HLARS method is adopted to estimate the PCE coefficients [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF].

C. Statistical Modeling procedure

The overall modeling of a given antenna undergoing random disturbances X comprises the following steps:

1. Initialize the ED  and the starting value of Ec_th; impose the RRMSE threshold (RRMSEth) for the SME truncation and the Leave-One-Out error threshold (e LOO_th ) for the PCE-HLARS algorithm; 2. Pre-evaluate antenna responses throughout the ED set.

Precisely, for each value x (n) , n = 1, …, NED, design the corresponding antenna, carry out the EM simulations or measurements and compute the ATF H (r , x (n) ); 3. Apply the SMEM to each ATF under the energy threshold Ec_th criterion and obtain the dominant index set I d (n) ; 4. Define the aggregated modes index set

I D = ⋃ I d (n) N ED n=1 ; 5.
Re-expand all ATFs on modes { q/2 (r ), q ∈ I D } and obtain corresponding VSH coefficients  q , q ∈ I D . Note that the steps 4 -5 are to guarantee that all expansions are based on the same set of modes, whereas the energy criterion Ec_th is still satisfied; 6. Reconstruct the ATFs from the truncated modes { q ,  q/2 (r ), q ∈ I D } and assess the error of reconstruction. If RRMSE ≤ RRMSEth is satisfied then pass to step 7, otherwise increase Ec_th and go back to the step 3;

7. Apply the PCE-HLARS method to each  q and finally get PCE models { q = ∑ y α q α∈ q Φ α X , q ∈ I D }; 8. Assess e LOO . If e LOO e LOO_th is satisfied for all PCE models, or if NED exceeds a predefined acceptable value NED,max (i.e. too costly to continue) then terminate the modeling procedure; otherwise enrich the ED (i.e.  new =[;  add ]) and go back to the step 2. Eventually, a surrogate model relating quantitatively and explicitly the random disturbances X and the ATFs via the whole set of polynomial coefficients {y α q , α ∈  q , q ∈ I D }is established, allowing the ATF generation for any value of X.

III. APPLICATION TO A DEFORMABLE DIPOLE

The proposed methodology is, among others, applied to a lossless half-wave wire dipole which is a simple and very well-known canonical antenna (so an ideal "test antenna") to verify its relevance, and assess its accuracy and robustness. Besides, a few results are briefly given for a deformable patch.

A. Design of the deformable dipole

The original antenna (Fig. 1 (a)) is designed around 2.3 GHz and well matched (|S11| < 10 dB) over a frequency band of [2.1, 2.4] GHz. The dipole is assumed to be deformable using three independent parameters: stretching (with stretch ratio ), bending (with curvature  or bending radius Rb = 1/), and torsion (with torsion rate α), see Fig. 1 (b)-(d). We suppose that the input parameters (, , α) are uniformly distributed over respectively [1, 1.05], [0, 0.1] mm -1 , and [0, 0.1] rad•mm -1 . Fig. 2 shows the combined effects of disturbances on the reflection coefficient |S11| and on the radiation patterns (at f = 2.3 GHz). Compared to the original dipole, variations of more than 15 dB in |S11| and 13.5 dB in Gr, i.e. 92.6 % in terms of RRMSE are observed; these large variations justify to resort to an accurate quantitative approach. In the following, the considered frequency is set to f = 2.3 GHz.

B. Modeling description and results

The modeling approach described in section II.C is applied to the dipole. Firstly, the ED is initialized with 100 samples using the Latin-Hypercube-Sampling (LHS) algorithm [START_REF] Mckay | A comparison of Three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]; the desired thresholds RRMSEth and eLOO_th are fixed at 2 % and 4 %; the starting value for Ec_th is set at 80 %. Secondly, the pre-evaluation of the dipole is carried out in CST ® MWS (with the FEM solver) to obtain radiation fields which are sampled every 5º in both elevation and azimuth. The full truncation degree of the SME is set at Nmax = 3 (i.e. Pmax = 30), for which the collected energy is proved to be very close to the total energy of the EMF in reality. The "Maximin" technique is used for enriching the initial ED, the principle of which is to add new points sequentially to the ED while maximizing the minimum Euclidean distance between the new points and the existing ones. At last, we use the UQLab, a toolbox for the uncertainty quantification developed by the team of B. Sudret [START_REF] Marelli | UQLab user manual -Polynomial Chaos Expansions[END_REF], to implement the PCE-HLARS algorithm.

As shown in Fig. 3, while increasing the collected energy threshold Ec_th, the RRMSE of the SME decreases gradually but meanwhile the number of selected dominant  q increases rapidly, which degrades the compactness of the model. The RRMSEth is ultimately satisfied with Ec_th = 99.9 %, with 35 (<< 2Pmax) selected dominant  q . Fig. 4 shows the e LOO of the PCE models for the selected dominant  q . As the mode indexing is chosen once, there is no way to guarantee a monotonic decrease of the modal energy with the global index q for all samples in the ED (as energy dominance can vary from one sample to another). It is why e LOO does not increase monotonically with q in Fig. 4 (although of course a less energetic mode is more sensitive to noise). Finally it can be seen that when NED increases the error decreases for all PCE models and the e LOO_th is satisfied for all models for NED = 260, for which we get a surrogate model with about 2000 polynomial coefficients.

In order to verify its accuracy and robustness, the constructed model has been used to generate the most deformed dipole (i.e. (, , α) = (1.05, 0.1, 0.1)) and a comparison of radiation patterns between the initial field, the reconstructed field (with dominant VSH modes), and the generated field (with the proposed surrogate model) is presented in Fig. 5(a). It is noted that the maximal difference in Gr is only 0.7 dB and the RRMSE is less than 8 % in this extreme deformation case; furthermore, the model has been used to generate dipoles for 500 random values of (, , α), which are compared to their simulated (with CST ® MWS) counterparts, as shown in Fig. 5(b). We see that the RRMSE of the generated antennas remains below 2 % for 88 % of the cases and the extreme error is still less than 10 %, which confirms the accuracy and the robustness of the model. It should be noted that there is always a trade-off between the model accuracy and the cost to achieve that precision. One can choose the most appropriate trade-off according to each practical needs. Another conclusion we can draw is that the model also verifies high parsimony. Compared to the initial antennas representation by a huge amount of sampled EMF data (~10 6 ), the surrogate model requires only a limited number of parameters (i.e. ~10 3 ), reaching a data compression rate of about 10 3 . It is actually much higher as input data (parameters or angle variables) can be interpolated as desired. Besides, the modeling was applied to other antennas such as e.g. a microstrip-fed patch with variable length and substrate curvature, thickness and r to show its feasibility for any type of antennas or deformations. For NED = 280, a model of ~ 2600 PCE coefficients and RRMSE < 15 % was achieved.

IV. CONCLUSION

The proposed methodology allows to develop compact and precise surrogate models of the FF radiated by "random antennas" with a minimized cost. This methodology has been applied to the modeling of a deformable dipole to prove its feasibility, efficiency and high accuracy. The derived model is computationally costless for generating deformable dipoles usable as surrogate models in various types of analyses requiring the complete field knowledge (including polarization and phase). The proposed modeling methodology is applicable to any other types of random disturbances or antennas. Future and already on-going work will address more realistic applications considering e.g. flexible printed antennas for terminals or sensors and/or perturbations due to the nearby surrounding environment, including the modeling of the frequency behavior (S11( f ), etc.).  q e LOO (%)

Fig. 1 .

 1 Fig. 1. Dipole undergoing different deformations: (a) original dipole; (b) stretched dipole; (c) curved dipole (d) curved and twisted dipole.

Fig. 2 .

 2 Fig. 2. Combined effect of stretching, bending, and torsion (, , α) on the dipole: (a) reflection coefficient |S11|; (b) radiation patterns (in E/H plane) at f = 2.3 GHz.

Fig. 3 .

 3 Fig. 3. (a) RRMSE of reconstructed fields; (b) number of selected dominant q: for different collected energy threshold Ec_th for NED = 100.

Fig. 4 .

 4 Fig. 4. e LOO of the PCE models for selected dominant q for different ED.

Fig. 5 .

 5 Fig. 5. (a) Radiation patterns (in E/H plane) of the initial, reconstructed and generated fields for the most deformed case at f = 2.3 GHz; (b) RRMSE of 500 generated dipoles for random values of (, , α).

ACKNOWLEDGMENT

The authors would like to thank Pr. B. Sudret and the UQLab team for providing PCE-HLARS freewares for academics, and Dr. J. Wiart for fruitful discussions and useful advices.


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

Philippe Wenger

Félix Majou

The Optimal Design of Three Degree-of-Freedom Parallel Mechanisms for Machining Applications

Keywords: Parallel Machine Tool, Isotropic Design, and Singularity

The subject of this paper is the optimal design of a parallel mechanism intended for three-axis machining applications. Parallel mechanisms are interesting alternative designs in this context but most of them are designed for three-or six-axis machining applications. In the last case, the position and the orientation of the tool are coupled and the shape of the workspace is complex. The aim of this paper is to use a simple parallel mechanism with two-degree-of-freedom (dof) for translational motions and to add one leg to have one-dof rotational motion. The kinematics and singular configurations are studied as well as an optimization method. The three-degree-of-freedom mechanisms analyzed in this paper can be extended to four-axis machines by adding a fourth axis in series with the first two.

Introduction

Parallel kinematic machines (PKM) are commonly claimed to offer several advantages over their serial counterparts, like high structural rigidity, high dynamic capacities and high accuracy [START_REF] Treib | Similarity laws of serial and parallel manipulators for machine tools[END_REF]. Thus, PKM are interesting alternative designs for high-speed machining applications. The first industrial application of PKMs was the Gough platform, designed in 1957 to test tyres [START_REF] Gough | Contribution to discussion of papers on research in automobile stability, control an tyre performance[END_REF]. PKMs have then been used for many years in flight simulators and robotic applications [START_REF] Merlet | The parallel robot[END_REF] because of their low moving mass and high dynamic performances [START_REF] Treib | Similarity laws of serial and parallel manipulators for machine tools[END_REF]. This is why parallel kinematic machine tools attract the interest of most researchers and companies. Since the first prototype presented in 1994 during the IMTS in Chicago by Gidding&Lewis (the VARIAX), many other prototypes have appeared.

To design a parallel mechanism, two important problems must be solved. The first one is the identification of singular configurations, which can be located inside the workspace. For a six-dof parallel mechanism, like the Gough-Stewart platform, the location of the singular configurations is very difficult to characterize and can change under small variations in the design parameters [START_REF] Merlet | The parallel robot[END_REF]. The second problem is the non-homogeneity of the performance indices (condition number, stiffness...) throughout the workspace [START_REF] Treib | Similarity laws of serial and parallel manipulators for machine tools[END_REF]. To the authors' knowledge, only one parallel mechanism is isotropic throughout the workspace [START_REF] Kong | Generation of parallel manipulators with three translational degrees of freedom based on screw theory[END_REF] but the legs are subject to bending. Moreover, this concept is limited to three-dof mechanisms and cannot be extended to four or five-dof parallel mechanisms. Numerous papers deal with the design of parallel mechanisms [START_REF] Kong | Generation of parallel manipulators with three translational degrees of freedom based on screw theory[END_REF][START_REF] Hervé | Structural Synthesis of Parallel Robots Generating Spatial Translation[END_REF]. However, there is a lack of four-or five-dof parallel mechanisms, which are especially required for machining applications [START_REF] Gao | New kinematic structures for 2-, 3-, 4-, and 5-DOF parallel manipulator designs[END_REF]. To decrease the cost of industrialization of new PKM and to reduce the problems of design, a modular strategy can be applied. The translational and rotational motions can be divided into two separated parts to produce a mechanism where the direct kinematic problem is decoupled. This simplification yields also some simplifications in the definition of the singular configurations. The organization of this paper is as follows. Next section presents design problems of parallel mechanisms. The kinematic description and singularity analysis of the parallel mechanism used, are reported in sections 3.1 and 3.2. Sections 3.3 and 3.4 are devoted to design and the optimization.

About parallel kinematic machines

General remarks

In a PKM, the tool is connected to the base through several kinematic chains or legs that are mounted in parallel. The legs are generally made of telescopic struts with fixed foot points (Figure 1a), or fixed length struts with moveable foot points (Figure 1b).

P P

Figure 1a: A bipod PKM Figure 1b: A biglide PKM For machining applications, the second architecture is more appropriate because the masses in motion are lower.

The linear joints can be actuated by means of linear motors or by conventional rotary motors with ball screws. A classification of the legs suitable to produce motions for parallel kinematic machines is provided by [START_REF] Gao | New kinematic structures for 2-, 3-, 4-, and 5-DOF parallel manipulator designs[END_REF] with their degrees of freedom and constraints. The connection of identical or different kinematic legs permits the authors to define two-, three-, four-and five-dof parallel mechanisms. However, it is not possible to remove one leg from a four-dof to produce a three-dof mechanism because no modular approach is used.

Singularities

The singular configurations (also called singularities) of a PKM may appear inside the workspace or at its boundaries. There are two types of singularities [START_REF] Gosselin | Singularity analysis of closedloop kinematic chains[END_REF]. A configuration where a finite tool velocity requires infinite joint rates is called a serial singularity. These configurations are located at the boundary of the workspace. A configuration where the tool cannot resist any effort and in turn, becomes uncontrollable is called a parallel singularity. Parallel singularities are particularly undesirable because they induce the following problems (i) a high increase in forces in joints and links, that may damage the structure, and (ii) a decrease of the mechanism stiffness that can lead to uncontrolled motions of the tool though actuated joints are locked. Figures 2a and2b show the singularities for the biglide mechanism of Fig. 1b. In Fig. 2a, we have a serial singularity. The velocity amplification factor along the vertical direction is null and the force amplification factor is infinite.

Figure 2b shows a parallel singularity. The velocity amplification factor is infinite along the vertical direction and the force amplification factor is close to zero. Note that a high velocity amplification factor is not necessarily desirable because the actuator encoder resolution is amplified and thus the accuracy is lower. The determination of the singular configurations for twodof mechanisms is very simple; conversely, for a six-dof mechanism like the Gough-Stewart platform, a mechanism with six-dof, the problem is very difficult [START_REF] Merlet | The parallel robot[END_REF]. With a modular architecture, when the position and the orientation of the mobile platform are decoupled, the determination of the singularities is easier.

P P

Kinetostatic performance of parallel mechanism

Various performance indices have been devised to assess the kinetostatic performances of serial and parallel mechanisms. The literature on performance indices is extremely rich to fit in the limits of this paper (service angle, dexterous workspace and manipulability…) [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems[END_REF].

The main problem of these performance indices is that they do not take into account the location of the tool frame. However, the Jacobian determinant depends on this location [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems[END_REF] and this location depends on the tool used.

To the authors' knowledge there is no parallel mechanism, suitable for machining, for which the kinetostatic performance indices are constant throughout the workspace (like the condition number or the stiffness…).

For a serial three-axis machine tool, a motion of an actuated joint yields the same motion of the tool (the transmission factors are equal to one). For a parallel machine, these motions are generally not equivalent.

When the mechanism is close to a parallel singularity, a small joint rate can generate a large velocity of the tool. This means that the positioning accuracy of the tool is lower in some directions for some configurations close to parallel singularities because the encoder resolution is amplified. In addition, a high velocity amplification factor in one direction is equivalent to a loss of stiffness in this direction. The manipulability ellipsoid of the Jacobian matrix of robotic manipulators was defined two decades ago [START_REF] Salisbury | Articulated Hands: Force Control and Kinematic Issues[END_REF]. The JJ -1 eigenvalues square roots, γ 1 and γ 2 , are the lengths of the semi-axes of the ellipse that define the two velocity amplification factors between the actuated joints velocities and the velocity vector t & (λ 1 = 1/γ 1 and

λ 2 = 1/γ 2 )
. For parallel mechanisms with only pure translation or pure rotation motions, the variations of these factors inside the Cartesian workspace can be limited by the following constraints

min max i i i λ λ λ ≤ ≤
Unfortunately, this concept is quite difficult to apply when the tool frame can produce both rotational and translational motions. In this case, indeed the Jacobian matrix is not homogeneous [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems[END_REF]. A first way to solve this problem is its normalization by computing its characteristic length [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems[END_REF][START_REF] Chablat | On the Kinetostatic Optimization of Revolute-Coupled Planar Manipulators[END_REF]. The second approach is to limit the values of all terms of the Jacobian matrix to avoid singular configuration and to associate these values to a physical measurement (See section 3.4).

3 Kinematics of mechanisms studied

Kinematics of a parallel mechanism for translational motions

The aim of this section is to define the kinematics and the singular configuration of a two-dof translational mechanism (Figure 3), which can be extended to threeaxis machines by adding a third axis in series with the first two. The output body is connected to the linear joints through a set of two parallelograms of equal lengths

i i L A B =
, so that it can move only in translation. The two legs are PPa identical chains, where P and Pa stand for Prismatic and Parallelogram joints, respectively. This mechanism can be optimized to have a workspace whose shape is close to a square workspace and the velocity amplification factors are bounded [START_REF] Chablat | Conception Isotropique d'une morphologie parallèle: Application à l'usinage[END_REF]. The joint variables 1 ρ and 2 ρ are associated with the two prismatic joints. The output variables are the Cartesian coordinates of the tool center point

T y x P ] [ = .
To control the orientation of the reference frame attached to P, two parallelograms can be used, which also increase the rigidity of the structure, Figure 3.

P B A P B A - in two possible directions, we obtain ) ( 1 1 1 1 a b i a p - × θ + = & & & (1a) ) ( 2 2 2 2 a b i a p - × θ + = & & & (
ρ = & & e a
, respectively. For an isotropic configuration to exist where the velocity amplification factors are equal to one, we must have 0 . 2 1 = e e [START_REF] Chablat | Conception Isotropique d'une morphologie parallèle: Application à l'usinage[END_REF] (Figure 5). We would like to eliminate the two passive joint rates ⊥e b a or ⊥e b a . To avoid these two singularities, the range limits are defined in using suitable bounds on the velocity factor amplification (See section 3.3).

Kinematics of a spatial parallel mechanism with one-dof of rotation

The aim of this section is to define the kinematics of a simple mechanism with two-dof of translation and onedof of rotation. To be modular, the direct kinematic problem must be decoupled between position and orientation equations. A decoupled version of the Gough-Stewart Platform exists but it is very difficult to build because three spherical joints must coincide [START_REF] Khalil | Kinematic Analysis and Singular configurations of a class of parallel robots[END_REF]. Thus, it cannot be used to perform milling applications. The main idea of the proposed architecture is to attach a new body with the tool frame to the mobile platform of the two-dof mechanism defined in the previous section. The new joint admits one or two-dofs according to the prescribed tasks.

To add one-dof on the mechanism defined in section 3, we introduce one revolute joint between the previous mobile platform and the tool frame. Only one leg is necessary to hold the tool frame in position. Figure 4 shows the mechanism obtained with two translational dofs and one rotational dof. The architecture of the leg added is PUU where P and U stand for Prismatic and Universal joints, respectively [START_REF] Gao | New kinematic structures for 2-, 3-, 4-, and 5-DOF parallel manipulator designs[END_REF].

The new prismatic joint is located orthogonaly to the first two prismatic joints. This location can be easily justified because on this configuration, i.e. when 1 1 -⊥b a b p , the third leg is far away from serial and parallel singularities. Let ρ & be referred to as the vector of actuated joint rates and p & as the velocity vector of point P,

1 2 2 [ ] T ρ ρ ρ = & & & & ρ and [ ] T x y = p & & &
Due to the architecture of the two-dof mechanism and the location of P, its velocity on the z-axis is equal to zero. p & can be written in three different ways by traversing the three chains A are given by We want to eliminate the passive joint rates i θ & and α & from Eqs. (3a-c), which we do upon dot-multiplying Eqs. (3a-c) by

P B A i i , 1 1 1 1 ( ) θ = + × - p a i b a & & & (3a) 2 2 2 2 ( ) θ = + × - p a i b a & & & (3b) 3 3 3 3 3 ( ) ( ) ( ) θ α β = + + × - + × - p a
1 1 1 ρ = & & e a ,
i i a b -, 1 1 1 1 1 1 ( ) ( ) T T ρ - = - b a p b a e & & (4a) 2 2 2 2 2 2 ( ) ( )
T T ρ - = - b a p b a e & & (4b) 3 3 3 3 3 3 3 3 3 ( ) ( ) ( ) ( )
T T T ρ β - = - + - × - b a p b a e b a j p b & & & (4c)
Equations (4a-c) can be cast in vector form, namely,

= t J& ρ with B A J 1 - = and T y x ] [ β = & & & t
where A and B are the parallel and serial Jacobian matrices, respectively, , i.e.,

T T T T ⎡ ⎤ - ⎢ ⎥ ≡ - ⎢ ⎥ ⎢ ⎥ - -- × - ⎣ ⎦ b a A b a b a b a j p b ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - - - ≡ 3 3 3 2 2 2 1 1 2 2 3 3 3 3 3 ( ) 0 ( ) 0 ( ) ( ) ( )
-⊥ a b e . However, these singular configurations are simple and can be avoided by proper limits on the actuated joints.

Optimization of the useful workspace for translational motions

Two types of workspaces can be defined, (i) the Cartesian workspace is the manipulator's workspace defined in the Cartesian space, and (ii) the useful workspace is defined as a subset of the Cartesian workspace. Workspace and size are prescribed where some performance indices are prescribed.

For parallel mechanism, the useful workspace shape should be similar to the one of classical serial machine tools, which is parallelepipedic if the machine has three translational degrees of freedom for instance. So, a square useful workspace is prescribed here where the velocity amplification factors remain under the prescribed values. Two square useful workspaces can be used, (i) he first one has horizontal and vertical sides (Figure 5a) and (ii) the second one has oblique sides but its size is higher (Figure 5b). (Figure 6) and the velocity amplification factors are computed for each configuration. This method was developed in [START_REF] Majou | A Novel method for the design of 2-DOF Parallel mechanisms for machining applications[END_REF]. In each case, velocity amplification factor extrema are located along the sides i j P P : they start from 1 at point S, then they vary until they reach prescribed boundaries (1/ 3 3 i λ ≤ ≤ ). Then computation (which analytical expressions λ i ) has been obtained with Maple along the four sides of the square and for the same length of the leg equal to one. For the first mechanism, solution (a), the size of the optimal surface is equal to 0,89 m 2 and for the second mechanism, solution (b), the size is equal to 0,62 m 2 . The result obtained for the solution (b) is smaller than for the solution (a) but is more appropriate for the extension three-axis mechanism of Figure 4. In effect, we want the axis of rotation to be parallel to one of the side of the useful workspace. In the next section, the lengths of the third leg will optimize to achieve this square useful workspace without singularity.

Optimization of the third-axis for rotational motions

The aim of this section is to define the two lengths of the leg, 1

L B P = and 2 3 3 L A B = 3
such that it is possible to achieve the maximum range variation of third axis β without meeting a singular configuration throughout the square workspace with the size defined in the previous section. The fist step of this optimization is to find the location of the prismatic joint. When P is on the center of the square workspace, we chose to place the third leg furthest away from singular configuration, i.e. when As it is defined in Section 3.2, the third leg is in a singular configuration whenever (Figure 7d). In the optimization function, we set: With a suitable chose of lengths, it is easy to obtain a range variation of β higher than 180°. So, this value can be reduced if we increase the constraint defined in Eqs. (4). However, when 2 1 1,8 L L = , the range of variation is optimal.

Conclusions

In this paper, a parallel mechanism with two degrees of position and one degree of rotation is studied. All the actuated joints are fixed prismatic joints, which can be actuated by means of linear motors or by conventional rotary motors with ball screws. Only three types of joints are used, i.e., prismatic, revolute and universal joints. All the singularities are characterized easily because position and orientation are decoupled for the direct kinematic problem and can be avoided by proper design. The lengths of the legs as well as their positions is optimized, to take into account the velocity amplification factors for the translational motions and to avoid the singular configuration for the rotational motions.

Figure

  Figure 2a: A serial singularity

Figure 3 :

 3 Figure 3: Parallel mechanism with two-dof To produce the third translational motion, it is possible to place orthogonally a third prismatic joint. The velocity p & of point P can be expressed in two different ways. By traversing the closed loop ) ( 2 2 1 1

1 θ & and 2 θ

 2 & from Eqs. (1a-b), which we do upon dot--b) can be cast in vector form, namelyρ & & B p A =, with A and B denoted, respectively, as the parallel and serial Jacobian matrices, is defined as the vector of actuated joint rates and p & is the velocity of point P, i.e., When A and B are not singular, we obtain the relations,

2 2 A

 22 B are colinear, i.e.

Figure 4 :

 4 Figure 4: Parallel mechanism with two-dof of translation and one-dof of rotation

T

  There are two new singularities when one leg is added. The first one is a parallel singularity when

Figure 5 :

 5 Figure 5: Cartesian workspace and isotropic configuration To find the best useful workspace (center locus and size), we can shift the useful workspace along x-axis ( ) x Δ and

Figure 6 :

 6 Figure 6: Looking for the best useful workspace (center locus and size)

Figure 7 :

 7 Figure 7: Optimal, serial and parallel configuration of the third leg

1 L 2 L 2 =1,8L 1 Figure 8 :

 12218 Figure 8: Range of variation in degrees of β as a function of 1 L and 2 L

& & & & & (3c)

Acknowledgments

This research was partially supported by the CNRS (Project ROBEA "Machine à Architecture compleXe").