Nageswaran Praveenganesh

Emilien Demory

Christine Gamon

Véronique Blandin

Pierre Y Chavant
email: yves.chavant@ujf-grenoble.fr

-1,3,2-dioxaborinane)

Keywords: boron, palladium, catalysis, arylboronic, Suzuki Efficient borylation of reactive aryl halides with MPBH (4, 4, 6-trimethyl

The combination of 4,4,6-trimethyl-1,3,2dioxaborinane, a particularly stable and inexpensive borylation reagent, and Buchwald's palladium catalyst provides a simple, fast, cost-effective borylation of electron-rich, reactive iodides, bromides, and triflates, to produce stable, easily purified boronic esters.

The increasing number of applications of the Suzuki-Miyaura coupling prompts the design of efficient and cost-effective methods for the preparation of the arylboronic precursors. Among these, the borylation of aryl halides or triflates with pinacolborane 1 (PinBH) is an excellent alternate to the more classical addition of organometallics to trialkoxyborane. Murata's team 2 and ours 3 proposed the replacement of PinBH in this reaction by the easily accessible, much less expensive and more conveniently prepared MPBH (4,4,6-trimethyl-1,3,2dioxaborinane 4 1). In addition, Billingsley and Buchwald 1i completed in 2008 the original study of Baudoin et al. 1c on the very positive role of bulky phosphine ligands in the borylation with PinBH. We wondered if the advantages of both implementations could be combined to provide a convenient and cost-effective access to aryl boronic esters from aryl halides; we discuss herein the advantages and limitations of this system. MPBH is a particularly stable dioxaborinane, 4 readily accessible from hexyleneglycol and borane. We described earlier 2,3 convenient preparations using BH 3 -DMS or B 2 H 6 generated from NaBH 4 and MeSO 3 H. We considered the use of N,N-diethylaniline-borane, a safe 5 source of borane for large-scale applications, but the reaction of hexylene glycol with this reagent was sluggish. Pd/C is known 6 to accelerate the decomplexation and solvolysis of amine-boranes in methanol. In an analogous manner (Scheme 1), we found that reacting Then we studied the palladium-catalyzed borylation of various aryl halides and triflates with MPBH and a catalyst prepared from Pd 2 (dba) 3 and CyJohnPhos 9 according to Baudoin's procedure 1c (Scheme 2, Table 1). 10 In preliminary experiments, we found that toluene and dioxane were equally suitable solvents. With PinBH, a dramatic feature of this catalyst system compared to others 1-3 is the extreme rapidity of the process, when applied to electron-rich aryl iodides 11 or bromides. 12 The same is true with MPBH. This high reactivity allows complete conversion of reactive aryl iodides and bromides (entries 1-4, 7, 8) between 20 and 40°C within a few hours. To our knowledge, operating at 20°C has not been described for this reaction yet. 13 The process remained efficient at low levels of catalyst: the boronate ester was produced from p-bromo-N,Ndimethylaniline in 90% yield in 6 h at 20°C when 0.3 mol% Pd was used (10 mmol run, entry 7). When less reactive iodides were used, some hydrodehalogenation was observed. This side-reaction is always met in analogous borylation with PinBH, and is favored by electro-withdrawing substituents on the aryl ring. In the present work as well as in the literature, 1 it is the only side-reaction: the homocoupling product was never observed. We indicate whenever available the Ar-B(OR) 2 /Ar-H ratio, but low-boiling arenes were seldom accurately measurable. Note that all reactions of Table 1 were run to complete conversion of the starting material. Thus, lower isolated yields reflect a higher level of hydrodehalogenation. With the reactive substrates of entries 1-4, 7 and 8, the reduced product was below detec-tion. In all cases, it was found to be less present than with the PdCl 2 (TPP) 2 catalyst in our former work. 3 It can be concluded that in the case of reactive, electron-rich aryl iodides and bromides, MPBH is a cost-effective substitute to pinacolborane in a fast and efficient palladium-catalyzed borylation. Disappointingly, with less reactive aryl halides, a clearcut difference in reactivity between PinBH and MPBH appeared, with MPBH favouring the unwanted hydrodehalogenation. The difference remained limited with electron-depleted iodides (entries 5 and 6). Entry 5 (69% yield) should be compared with Masuda's 1b result: borylation of ethyl 4-iodobenzoate with PinBH in the presence of PdCl 2 dppf yielded 79% Ar-Pin (hydrodehalogenated product 17%). More dramatic differences appeared with electron-depleted bromides, for which hydrodehalogenation became deleterious (entries 9-12). In the same conditions, PinBH was much more efficient: compare entries 9 and 10, 11 and 12. The borylation of chlorides was practically limited to the most reactive, electron-rich ones (entries 15-17). [START_REF] Very Recently | borylated aryl chlorides with both PinBH and MPBH in the presence of excess Bu 4 NI[END_REF] Remarkably, aryl triflates 1b,15 were borylated efficiently (entries 13-14). The extent of the hydrodetriflation seemed less dependent of the electro-donating (entry 13) or electro-withdrawing nature of the substituents (entry 14, compare also with entry 11). So far, all reported optimizations of this reaction were done with electron-rich substrates. [1][2][3] We decided to check the influence of reaction parameters on the borylation-hydrodehalogenation ratio when the substrate is electron-depleted; 2-bromonaphtalene was our model of "unreactive" substrate. 16 Reactions were monitored by GC with internal standard. Although we could not increase the yield in boronic ester 2g above 50%, some interesting trends appeared. Higher temperatures favor borylation over hydrodehalogenation (37% yield in 2g at 80°C, 50% at 150°C), making this reaction a good candidate for closed-vessel, microwaves-heated reactions. Nevertheless, at identical temperatures, microwaves heating did not lead to significantly faster reaction than conventional heating. With 2-bromonaphtalene and MPBH, other ligand systems (dppf, 1a,b dpePhos, 1f excess TPP, 3 N-(dicyclohexylphosphino)-2-(2'-tolyl)indole 17 or di-tert-butylbipyridyl) led to slower reactions and increased hydrodehalogenation. Monitoring the course of the reaction showed that the borylation/hydrodehalogenation ratio is lower at early stages of the reaction. Examples have been reported with PinBH 18 where prolonged heating of the reaction mixtures caused an increase of the hydrodehalogenation. In the present conditions, both reactions are concurrent. Prolonged heating after the end of the conversion of the bromide never caused any significant loss of borylated product. A larger excess of MPBH (3 eq.) slightly favored hydrodehalogenation: best results were obtained with only 1.5 eq. dialkoxyborane.

Our hypothesis is that MPBH is more prone than PinBH to some disproportionation [START_REF] Rose | The disproportionation of dialkoxyboranes yields trialkoxyboranes and BH 3 . See[END_REF] in the reaction conditions, releasing small amounts of BH 3 that would favor Pd-catalyzed reduction of the substrate. Indeed, replacing MPBH by BH 3 :Et 3 N caused extensive hydrodehalogenation. As a conclusion, the combination of MPBH, a particularly stable and inexpensive reagent, and Buchwald's palladium catalyst provides a simple, very fast, costeffective borylation of reactive, electron-rich aryl iodides, bromides, and triflates, to produce stable, easily purified boronic esters, that can be readily used as such in Suzuki couplings. 1,2 Scheme: Borylation of 2-bromonaphtalene

MPBH 2a-j Scheme 2

  Borylation of aryl halides and triflates.

	presence of Pd/C (5 % by weight of 10 % Pd/C reagent)
	resulted in the complete consumption of the diol in 30
	min at room temperature. Distillation of the reaction mixture in vacuo provided neat MPBH in 75% yield. 7,8
		Cy 2 P
	X X = I, Br, Cl, OTf R +	20-100 o C Et 3 N Pd 2 (dba) 3	R	B	O O

N, N-diethylaniline-borane and hexyleneglycol in the

Table 1

 1 Borylation of various aryl halides and triflates with MPBH catalyzed by Pd/CyJohnPhos. a,10

	Entry substrate		conditions	ArB(O	Product
					R)2/Ar	yield
					H b	(%) c
	1	Me 2 N		3% Pd, 20°C,		2a
				<0.5 h	> 99:1	90
			I
	2	MeO		1% Pd, 20°C,		2b
				1 h	> 95:5	91
			I
	3	H 2 N		0.5% Pd,		2c
				20°C,	> 97:3	79
			I	0.5 h
	4	NH 2		1% Pd, 20°C,		2d
				3 h	> 97:3	88
		I
	5	MeOOC		1% Pd, 80°C,		2e
				3 h	70:30	69
			I
	6			1% Pd, 40°C,		2f
				16 h	n.d. d	60
		O 2 N	I
	7	Me 2 N		0.3% Pd, 20°C,	> 99:1	2a 90 e
		Br	6 h
	8	MeO		0.5% Pd,		2b
				40°C,	>95:5	81
		Br	6 h
	9			MPBH (3 eq.)		2g
				1% Pd, 100°C	30:70	_
		Br	Tol. f , 30 h
	10			PinBH (3 eq.)	92:8	3g
				1% Pd, 100°C		_
				Tol. f , 15 h
		O		MPBH		2h
	11	Ph		Tol. f , 15min g 1% Pd, 150°C	29:71	_
			Br	PinBH
	12			1% Pd, 150°C	70:30	3h
				Tol. f , 15min g		litt. 1i 70
	13	Me		1% Pd, 80°C,		2i
				6 h	n.d.	82
		OTf
	14	O		1% Pd, 80°C,		2i
		Me		16h	n.d.	73
		OTf
	15	MeO		1% Pd,		2b
				100°C, 16h	>95:5	90
		Cl
	16	OMe		2% Pd, 80°C,		2j
				16 h	51:49	45
		Cl
	17	Me		1 % Pd, 80°C,		2i
				24 h	n.d.	62
		Cl
	a Unless otherwise stated: MPBH 1.5 eq., Et3N 3 eq., Pd (as
	Pd2(dba)3), CyJohnPhos 2 eq. per Pd, in dioxane.

b

Measured by GC, except entries 3 and 4: 1 H-NMR. c Isolated yields after chromatography. d n.d.: not determined. f run in toluene. g sealed vessel, microwaves heating.

Acknowledgment

We gratefully thank the Université J. Fourier and the CNRS for financial support.

Int. Ed. 1985, 24, 878; (f) PraveenGanesh, N.; d'Hondt, S.; Chavant, P. Y. J. Org. Chem. 2007, 72, 4510. (5) Wilkinson, H. S.; Tanoury, G. J.; Wald, S. A.; Senanayake, C. H. Org. Process Res. Dev. 2002, 6, 146. (6) (a) Couturier, M.; Tucker, J. L.; Andresen, B. M.; Dube, P.; Negri, J. T. Org. Lett. 2001, 3, 465; (b) Couturier, M.; Tucker, J. L.; Andresen, B. M.; Dube, P.; Brenek, S. J.; Negri, J. T. Tetrahedron Lett. 2001, 42, 2285. (7) Preparation of MPBH: under a well-ventilated fume hood (H 2 evolution), a 100-mL flask equipped with a Claisen distillation head was charged with 320 mg of 10 % Pd/C and anhydrous dioxane (10 mL). The flask was flushed with nitrogen and freshly distilled 2-methyl-2,4pentanediol (50 mmol, 6.4 g) in dioxane (10 mL) was added. A solution of N, N-diethylaniline-borane (60 mmol, 9.78 g) in dioxane (10 mL) was added at 20°C over 30 min. The reaction mixture was stirred for an additional 30 min at 20°C. Low pressure distillation gave neat MPBH (4.8 g, 75%, b.p. 45 °C / 35 Torr). (8) Our attempts to use N, N-diethylaniline as the base in subsequent borylations failed. (9) [1,1'-biphenyl]-2-yl-dicyclohexyl-phosphine; Wolfe, J. P.; Singer, R. A.; Yang, B. H.; Buchwald, S. L. J. Am. Chem. Soc. 1999, 121, 9550. (10) Typical procedure for the borylation: an oven-dried Schlenk vessel (or a 10-mL microwaves vial) was charged with Pd 2 (dba) 3 (2.3 mg, 2.5 µmol, described as 0.5 % in Table 1) and CyJohnPhos (3.5 mg, 10 µmol, always 2 eq./Pd) and placed under an atmosphere of Argon. Anhydrous dioxane (0.6 mL), the aryl halide (0.5 mmol), triethylamine (152 mg, 1.5 mmol) and MPBH (96 mg, 0.75 mmol) were introduced (solid aryl halides were added along with the other solid reagents). The reaction mixture was then heated at the indicated temperature until the aryl halide has been completely consumed as determined by gas chromatography. The reaction was allowed to cool to room temperature, and filtered through a thin pad of Celite (eluent ether). The eluate was concentrated and the crude material purified via flash chromatography on silica gel. Products 2a-j 2,3 , 3g 1b and 3h 1i have been previously described. 4,4,6-trimethyl-2-naphthalen-2-yl-[1, 3,2]dioxaborinane 2g:

1 H NMR (400 MHz, CDCl 3 ) δ ppm 8.34 (s, 1H), 7. 89-7.85 (m, 2H), 7.82-7.76 (m, 2H), 7.48-7.40 (m, 2H), 4.38 (dqd, J = 11.8, 6.2, 2.9 Hz, 1H), 1.87 (dd, J = 13.7, 2.9 Hz, 1H), 1.62 (dd, J = 13.7, 11.8 Hz, 1H), 1.41 (s, 3H), 1.39 (s, 3H), 1.38 (d, J = 6.2 Hz, 3H); 13 C NMR (75.5 MHz, CDCl 3 ) δ ppm 134.8 (x2), 132.99, 130.14, 128.68, 127.67, 129.5 (very broad, C-B), 126.62, 126.39, 125.44, 71.16, 65.14, 46.15, 31.38, 28.27, 23.30; 11 B NMR (128 MHz, CDCl 3 ) δ ppm 27.10; IR (neat, cm -1 ): 3053, 2972, 2932, 2911, 1630, 1598, 1503, 861, 823, 765, 746, 685; MS (EI, 70eV) 255 (26), 254 (100), 253 (28). (11) (a) Baudoin, O.; Decor, A.; Cesario, M.; Gueritte, F. Synlett 2003, 2009; (b) Joncour, A.; Decor, A.; Thoret, S.; Chiaroni, A.; Baudoin, O. Angew. Chem. Int. Ed. 2006, 45, 4149; (c) Poriel, C.; Lachia, M.; Wilson, C.; Davies, J. R.; Moody, C. J. J. Org. Chem. 2007, 72, 2978. (12) (a) Combs, A.; Zhu, W.; Crawley, M.; Glass, B.; Polam, P.; Sparks, R. B.; Modi, D.; Takvorian, A.; McLaughlin, E.; Yue,E.; Wasserman, Z.; Bower, M.; Wei, M.; Rupar, M.; Ala, P. J.; Reid, B. M.; Ellis, D.; Gonneville, L.; Emm, T.; Taylor, N.; Yeleswaram, S.; Li, Y.; Wynn, R.; Burn, T. C.; Hollis, G.; Liu, P. C. C.; Metcalf, B. J. Med. Chem. 2006, 49, 3774; (b) Mentzel, U. V.; Tanner, D.; Tonder, J. E. J. Org. Chem. 2006, 71, 5807. (13) For related room-temperature Suzuki couplings, see: (a) Wolfe, J. P.; Tomori, H Sadighi, J. P.; Yin, J.; Buchwald, S. L. J. Org. Chem. 2000, 65, 1158; (b) Billingsley, K. L.;

Supporting Information for this article is available online at http://www.thieme-connect.com/ejournals/toc/synlett.

Supporting Information to

Efficient borylation of reactive aryl halides with MPBH (4,4,6-trimethyl-1,3,2-dioxaborinane).

DOI:

Nageswaran PraveenGanesh, Emilien Demory, Christine Gamon, Véronique Blandin, Pierre Y. This article describes our study on borylation of aryl halides and triflates with 4,4,6-trimethyl-1,3,2-dioxaborinane (MPBH), a stable and easily available reagent. The combination of MPBH and Buchwald's palladium catalyst provides a simple, fast, cost-effective borylation of electron-rich, reactive aryl iodides, bromides, and triflates, to produce stable, easily purified boronic esters.

In the course of this work we found that the desired borylation can be hampered by a concurrent hydrodehalogenation of the substrate. We studied the influence of various parameters on the extent of this side-reaction in the case of 2bromonaphtalene as the model. The results of this study are summarized in the table below.

Standard procedure for the borylation of 2-bromonaphtalene

An oven-dried Schlenk vessel (or a 10-mL microwaves vial) was charged with 2-bromonaphtalene (0.5 mmol), Pd 2 (dba) 3 (2.3 mg, 0.5 %) and CyJohnPhos (3.5 mg, 2 %) and placed under an atmosphere of Argon. Anhydrous toluene (1 mL), triethylamine (152 mg, 1.5 mmol), MPBH (96 mg, 0.75 mmol) and tetradecane (as an internal standard) were introduced with a syringe. The reaction mixture was then heated at the indicated temperature for the indicated time. When microwaves were used the reaction was hydrolysed after this time and analyzed by GC (relative response factors of the components were determined from mixtures of pure samples). Reactions performed in an oil bath at atmospheric pressure were monitored by repeated analysis of small aliquots.

Typical microwaves experiment: heating profiles

The profiles below (entry 31) are typical. This reaction was performed on a 0.5 mmol scale in a 10 mL vessel; solvent was toluene.


Konstantin Avrachenkov
email: k.avrachenkov@sophia.inria.fr

Bruno Ribeiro
email: ribeiro@cs.umass.edu

Don Towsley
email: towsley@cs.umass.edu

Improving Random Walk Estimation Accuracy with Uniform Restarts

Keywords: Sampling, Random Walk, Spectral Gap, PageRank, Online Social Network Echantillonnage, Marche Alátoire, Ecart Spectral, PageRank, Réseaux Sociaux

This work proposes and studies the properties of a hybrid sampling scheme that mixes independent uniform node sampling and random walk (RW)-based crawling. We show that our sampling method combines the strengths of both uniform and RW sampling while minimizing their drawbacks. In particular, our method increases the spectral gap of the random walk, and hence, accelerates convergence to the stationary distribution. The proposed method resembles PageRank but unlike PageRank preserves time-reversibility. Applying our hybrid RW to the problem of estimating degree distributions of graphs shows promising results.

Introduction

Many networks, including on-line social networks (OSNs) and peer-to-peer (P2P) networks, exist for which it is impossible to obtain a complete picture of the network. This leaves researchers with the need to develop sampling techniques for characterizing and searching large networks. Sampling methods can be classified as based on independent uniform sampling or crawling. These two classes of sampling methods have their advantages and drawbacks. Our work proposes and studies the properties of a hybrid sampling scheme that mixes independent uniform node sampling and random walk (RW)-based crawling. We show that our sampling method combines the strengths of both uniform and RW sampling while minimizing their drawbacks.

Within the class of uniform sampling methods, uniform node sampling is widely popular and has the advantage of sampling disconnected graphs. In an online social network (OSN) where users are associated with unique numeric IDs, uniform node sampling is performed by querying randomly generated IDs. In a P2P network like Bittorrent, uniform node sampling is performed by querying a tracker server [START_REF] Konrath | Attacking a swarm with a band of liars: evaluating the impact of attacks on bittorrent[END_REF]. In practice, however, these samples are expensive (resource-wise) operations (the ID space in an OSN, such as Facebook and MySpace, is large and sparse and tracker queries can be rate-limited [START_REF] Konrath | Attacking a swarm with a band of liars: evaluating the impact of attacks on bittorrent[END_REF]). For instance, in MySpace we expect only 10% the IDs to belong to valid users [START_REF] Gauvin | Measurement and genderspecific analysis of user publishing characteristics on myspace[END_REF], i.e., only one in every ten queries successfully finds a valid MySpace account.

Within crawl-based sampling methods, random walk (RW) sampling is among the most popular methods [START_REF] Bisnik | Optimizing random walk search algorithms in p2p networks[END_REF][START_REF] Gkantsidis | Hybrid search schemes for unstructured peer-to-peer networks[END_REF][START_REF] Gkantsidis | Random walks in peer-to-peer networks: algorithms and evaluation[END_REF][START_REF] Lv | Search and replication in unstructured peerto-peer networks[END_REF][START_REF] Ribeiro | Estimating and sampling graphs with multidimensional random walks[END_REF][START_REF] Volz | Probability based estimation theory for Respondent-Driven Sampling[END_REF]. Let G = (V, E) be an undirected, non-bipartite graph with n nodes. RW sampling is preferred because it requires few resources and, when G is connected, can be shown to produce asymptotically unbiased estimates of

f (G) = ∀v∈V h(v).
(

Moreover, when G is connected, a RW visits all nodes in G in O(n 3 ) [START_REF] Lovász | Random walks on graphs: a survey[END_REF] steps w.h.p., a useful property when searching unstructured networks (such as P2P networks). Note that the above formal RW guarantees require G to be connected. In the real-world, however, networks may consist of several disconnected components, e.g. Twitter [22] and Livejournal [START_REF] Ribeiro | Estimating and sampling graphs with multidimensional random walks[END_REF], to cite two known examples. Moreover, the performance of such methods are closely tied to the difference between the largest and the second largest eigenvalues of the associated RW transition probability matrix. This difference is also known as the spectral gap which we denote as δ. More precisely, let (X t : t = 0, 1, 2, . . . ) be a discrete-time Markov chain associated with a random walk over G with transition probability matrix P = [p ij ], ∀i, j ∈ V , where p ij = 1/d i and d i is the degree of node i ∈ V . The eigenvalues of P are

1 = λ 1 ≥ λ 2 ≥ λ 3 ≥ • • • ≥ λ n ≥ -1. The spectral gap is defined as δ = (1 -λ 2 ). When G is disconnected δ = 0.
Typically, even a connected real complex network has small δ, explained by the clustered structure of these networks [START_REF] Leskovec | Statistical properties of community structure in large social and information networks[END_REF], and as a consequence the RW tends to get "trapped" inside subgraphs.

Unfortunately the accuracy of estimates of f (G) (eq. ( 1)) obtained with a RW is inversely proportional to the spectral gap δ (revisited later in this section). The main contribution of this work is to increase the spectral gap, δ, of a RW by combining it with independent uniform node sampling. Although node sampling can be very expensive, when combined with RW sampling, the resulting algorithm can significantly reduce the estimation error and search time at a negligible increase in overhead. The idea is as follows, add auxiliary edges among all the nodes with weight α/n > 0. The hybrid sampling scheme, which we call RWuR, corresponds to a random walk on this modified graph. Our results show a significant increase in the spectral gap δ when our hybrid scheme is used (even when α is small). Note that when α = 0 we RWuR is a regular RW and in the limit when α → ∞ RWuR performs independent uniform sampling.

In what follows we revisit the connection between the spectral gap, δ, and estimation errors; connecting δ with the Mean Squared Error (MSE) of the RW estimates.

Spectral Gap and Estimation Error

For now we assume G is connected. Suppose we wish to estimate f (G) of eq. ( 1) from B sampled nodes obtained by a stationary RW, (X 0 , X 1 , . . . , X B-1 ). A widely used estimator of f (G) is [START_REF] Ribeiro | Estimating and sampling graphs with multidimensional random walks[END_REF][START_REF] Volz | Probability based estimation theory for Respondent-Driven Sampling[END_REF]

f = 1 B B-1 t=0 h ′ (X t ) , X t ∈ V (2)
where

h ′ (v) = h(v)/π v , v ∈ V ,
and π v is the stationary distribution of the random walk. We note that if the graph is undirected, the RW is time reversible and its stationary distribution is given by

π i = d i 2|E| , ∀i ∈ V. (3)
To simplify the notation we drop the dependence of f on G. The MSE of f is given by

E[( f -f ) 2
]. Now we explore the connection between the spectral gap, δ, and the MSE of f . In a stationary RW, f is an unbiased estimator of f [START_REF] Volz | Probability based estimation theory for Respondent-Driven Sampling[END_REF]. Thus, the Mean Squared Error (MSE) is also the variance of f . Let var π ( f ) denote the variance of f in a stationary RW and

E[( f -f ) 2 ] = var( f ) + bias( f ), where bias( f ) = (E[ f ] -f ) 2 .
The assymptotic ratio between (var( f ) + bias( f )) and var π ( f ) is a function of the spectral gap δ [1, Chapter 4.1]

sup f lim B→∞ var( f ) + bias( f ) var π ( f ) = 1 + λ 2 1 -λ 2 = 2 -δ δ .
Note that δ also determines the mixing time of the RW [START_REF] Sinclair | Improved bounds for mixing rates of Markov chains and multicommodity flow[END_REF], which means that when δ ≪ 1 (i.e., λ 2 ≈ 1) it takes many steps for the random walk to converge to the stationary distribution (a potential source of bias in f when the RW does not start in steady state).

INRIA

We now turn our attention to finding the relationship between var π ( f ) and δ. Consider h from eq. ( 1). Assume that ∀v∈V h(v) = 0 (this assumption is used to simplify our calculations) and h

′ 2 2 ∆ = ∀v∈V (h ′ (v)) 2 π v > 0. Let σ 2 ∆ = lim B→∞ B (var( f ) + bias( f )).
From the inequality [1, Chapter 4, Proposition 29]

σ 2 B 1 - 2 δB ≤ var π ( f ) ≤ 2 h ′ 2 2 δB 1 + δ 2B (4)
we have a relationship between the MSE of ĝ and the spectral gap δ. Eq. ( 4) shows that the MSE upper bound decreases as δ increases. One important class of f (G) functions is the class that measures θ k , the fraction of nodes with degree k,

f k (G) ∆ = ∀v∈V h k (v) ∆ = ∀v∈V 1(d v = k)/n, ,
where

1(x = y) = 1 if x = y and 1(x = y) = 0, otherwise. Note that h k 2 2 = θ 2 k /Π k , where Π k = ∀v∈V π v 1(d v = k), which yields var π ( fk ) ≤ 2θ 2 k Π k δB 1 + δ 2B . (5)
Thus, eq. ( 5) shows that the error in estimating the fraction of nodes with degree k is uper bounded by the inverse of: (1) the spectral gap δ and (2) the probability that the RW finds a node with degree k, Π k . In Section 2 we see that increasing parameter α of our hybrid RW increases δ but also decreases Π k when k is larger than the average degree. This tradeoff can be seen in the experiments of Section 3, where the MSE of the fraction of high degree nodes, at first, decreases as we increase α until a certain (unknown optimal) point where a further increase in α increases the MSE. As future work we will investigate this optimal value of α.

Reducing mixing time by restart

As we have observed in the previous section, many complex networks have small spectral gaps and hence random walks on such networks can have a negative impact on the accuracy (variance) and bias of the estimates.

In this work we are interested in methods for accelerating the rate of convergence to the stationary distribution based on the addition of auxiliary transitions. Thus, the natural first method to investigate is PageRank [START_REF] Brin | The anatomy of a large-scale hypertextual Web search engine[END_REF]. Namely, a random walk follows some outgoing link with probability c and with probability 1 -c it jumps to an arbitrary node of the network chosen according to the uniform distribution. Then, the modified random walk is described by the following transition matrix

P = cP + (1 -c) 1 n 1 1 T , ( 6
)
where 1 is a vector of ones with an appropriate dimension. Then, the stationary distribution of the modified random walk π(c) is a unique solution of the following equations

π(c) P = π(c), π(c)1 = 1.
It is known [START_REF] Haveliwala | The second eigenvalue of the Google matrix[END_REF] that the second largest eigenvalue of matrix P is equal to c. Thus, by choosing c not close to one, we can significantly increase the spectral grap δ. However, as was observed for example in [START_REF] Litvak | Characterization of tail dependence for in-degree and pagerank[END_REF], PageRank's steady state distribution can be just weakly correlated with node degree. Furthermore, there are cases when ranking of nodes according to PageRank is sensitive to the change of parameter c [START_REF] Bressan | Choose the damping, choose the ranking?[END_REF]. Thus, the stationary distribution of the modified random walk can be significantly distorted. To mitigate the latter problem, we suggest a variation where we connect all the nodes in the graph with a weight α/n. The difference with PageRank is that in our variation the uniform restart occurs not with a fixed probability but with a probability depending on the node degree. Specifically, the transition probability in our modification is given by pij

= α/n+1 di+α , if i has a link to j, α/n di+α , if i does not have a link to j. ( 7
)
The advantage of such a modification is that the new random walk is also reversible with the following stationary distribution πi (α) =

d i + α 2|E| + nα ∀i ∈ V, ( 8
)
from which the original stationary distribution (3) can easily be retrieved.

Let us now show that this second method also improves algebraic connectivity. First, we consider a regular graph where all nodes have degree d. Theorem 2.1 Let G = (V, E) be an undirected regular graph with degree d. Let λk (α) and λ k be the eigenvalues of the Markov chains associated with the modified and original random walks on G, respectively. Then, all the eigenvalues corresponding to the modified random walk except the unit eigenvalue are scaled as follows:

λk (α) = d d + α λ k , k = 2, ..., n. ( 9
)
Proof. The transition probabilities of the modified random walk ( 7) can be written in the following matrix form

P = d d + α P + α dn 11 T ,

INRIA

where P is the transition matrix corresponding to the original random walk. Since 1 is an eigenvector corresponding to the unit eigenvalue, we can apply Brauer's Theorem [START_REF] Brauer | Limits for the characteristic roots of a matrix, iv: Applications to stochastic matrices[END_REF]. Brauer's Theorem says that if λ is an eigenvalue and x the corresponding eigenvector of matrix A then λ + v T x is an eigenvalue of matrix A + xv T for any vector v and the other eigenvalues of A + xv T coincide with the eigenvalues of A. Thus, matrix P + α dn 11 T has the following eigenvalues: 1+α/d, λ 2 , ..., λ n . Since the eigenvalues of a matrix multiplied by a scalar are the eigenvalues of that matrix multiplied by the same scalar, the eigenvalues of matrix P are: 1, d/(d + α)λ 2 , ..., d/(d + α)λ n .

If we expand d/(d + α) in ( 9) into a power series with respect to α we can rewrite (9) as follows:

λ2 (α) = 1 - α d λ 2 + o(α).
Thus, for small values of α the spectral gap can be approximated as follows:

δ ≈ α d λ 2 . ( 10
)
Let us now consider the case of a general undirected graph.

Theorem 2.2 Let G = (V, E) be a general undirected graph. Then, the second largest eigenvalue λ2 (α) of the modified random walk has the following connection with the second largest eigenvalue λ 2 of the original random walk

λ2 (α) = 1 - n k=1 1 d k u 2k v 2k n k=1 u 2k v 2k α λ 2 + n k=1 1 d k u 2k n j=1 v 2j n n k=1 u 2k v 2k α + o(α), ( 11
)
where u 2 and v 2 are respectively left and right Fiedler eigenvectors of the original graph.

Proof. Let us analyze the equation

P (α)v 2 (α) = λ2 (α)v 2 (α) (12)
with the help of perturbation theory techniques [START_REF] Avrachenkov | Analytic perturbation theory and its applications[END_REF][START_REF] Baumg Artel | Analytic perturbation theory for matrices and operators[END_REF]. We expand P (α) as a power series with respect to α. Namely, we write P (α) = P + αΓ (1) + α 2 Γ (2) + ... , [START_REF] Haveliwala | The second eigenvalue of the Google matrix[END_REF] where in particular the coefficient of the first order term is given by

Γ (1) = diag 1 d k 1 n 11 T -P , ( 14
)
where diag 1 d k is a diagonal matrix with the elements 1 d k on the diagonal. We also expand λ2 (α) and v2 (α) as power series λ2 (α) = λ (0) + αλ (1) + α 2 λ (2) + ... [START_REF] Leskovec | Statistical properties of community structure in large social and information networks[END_REF] and v2 (α) = v (0) + αv (1) + α 2 v (2) + ... . [START_REF] Litvak | Characterization of tail dependence for in-degree and pagerank[END_REF] Next we substitute P (α), λ2 (α) and v2 (α) in the form of power series ( 13), ( 15) and ( 16) into equation [START_REF] Gkantsidis | Random walks in peer-to-peer networks: algorithms and evaluation[END_REF]. Thus, we have (P + αΓ (1) + α 2 Γ (2) + ...)(v (0) + αv (1) + α 2 v (2) + ...) = (λ (0) + αλ (1) + α 2 λ (2) + ...)(v (0) + αv (1) + α 2 v (2) + ...).

Equating terms with the same powers of α yields

P v (0) = λ (0) v (0) . ( 17
)
Since we are interested in the second largest eigenvalue of P (α), we conclude that v (0) = v 2 and λ (0) = λ 2 , where v 2 is the eigenvector of P corresponding to the second largest eigenvalue λ 2 of P . Then, collecting terms with α, we get

Γ (1) v 2 + P v (1) = λ (1) v 2 + λ 2 v (1) . ( 18
)
Premultiplication of equation ( 18) by the left eigenvector u T 2 corresponding to the second largest eigenvalue λ 2 (u

T 2 P = λ 2 u T 2 ) leads to u T 2 Γ (1) v 2 = λ (1) 2 u T 2 v 2 .
or

λ (1) = u T 2 Γ (1) v 2 u T 2 v 2 λ 2
Let us consider in more detail the expression u T 2 Γ (1) v 2 .

u T 2 Γ (1) v 2 = 1 n u T 2 diag 1 d k 1 1 T v 2 -u T 2 diag 1 d k P v 2 = 1 n n k=1 1 d k u 2k n j=1 v 2j - n k=1 1 d k u 2k v 2k λ 2
In the latter equality we use the fact that P v 2 = λ 2 v 2 . Thus, we obtain formula [START_REF] Gkantsidis | Hybrid search schemes for unstructured peer-to-peer networks[END_REF].

Similar to the above analysis we can derive any number of terms in the power series expressions for λ(α) and v(α). For example, equating the terms with α 2 in the perturbed equation yields (1) , where v (1) can be found from equation [START_REF] Lv | Search and replication in unstructured peerto-peer networks[END_REF]. Namely, v (1) is a solution of the equation

λ (2) = 1 u T 2 v 2 u T 2 Γ 2 v 2 λ 2 -u T 2 (Γ 1 P + λ (1) I)v
(P -λ 2 I)v (1) = (Γ 1 P + λ (1) I)v 2 ,
INRIA in the orthogonal complement to the subspaces corresponding to the eigenvalues 1 and λ 2 .

We can obtain a solution with the help of the group reduced resolvent

H {1,λ2} = 1 2πi Γ 1 ζ (P -ζI) -1 dζ,
where Γ is a contour in the complex plane enclosing all the eigenvalues of P except 1 and λ 2 [START_REF] Avrachenkov | Analytic perturbation theory and its applications[END_REF]. Then, we have v (1) = H {1,λ2} (Γ

1 P + λ (1) I)v 2 .
Even though Theorem 2.2 provides a connection between the second largest eigenvalues of the original and modified graphs, we cannot readily deduce from expression [START_REF] Gkantsidis | Hybrid search schemes for unstructured peer-to-peer networks[END_REF] if the spectral gap actually decreases. Therefore, next we analyse a typical case where we can obtain more insight from formula [START_REF] Gkantsidis | Hybrid search schemes for unstructured peer-to-peer networks[END_REF].

We note that v 2 and u T 2 are Fiedler vectors. The Fiedler vectors indicate principal clusters of the original graph. Let us represent the transition matrix for the original graph in the following form

P = P (0) + εC = P 1 0 0 P 2 + εC, ( 19
)
where P 1 , P 2 represent the transitions inside the principal clusters and εC represents transitions between the principal clusters. We choose the blocks P 1 and P 2 to be transition matrices, which means that some elements of εC corresponding to the blocks P 1 and P 2 are negative. Of course, all elements of the sum are non-negative. Now we are ready to state the next result.

Theorem 2.3

Given that the original graph has two principal Fiedler components with the same average node degree d, the following connection between the eigenvalues of the modified and original graphs take place

λ2 (α) = 1 - α d λ 2 + o(α) + O(ε, |E 1 | -1 , |E 2 | -1 ). ( 20
)
Proof. The proof is postponed to the Appendix.

In particular, we conclude from [START_REF] Ribeiro | Estimating and sampling graphs with multidimensional random walks[END_REF] that for small values of α the value of the spectral gap can be approximated as follows:

δ ≈ α d λ 2 . ( 21
)
Comparing ( 10) and ( 21), it is curious to observe that in the case of the general undirected graph the parameter d in ( 10) is replaced by the average node degree d.

The expression [START_REF] Sinclair | Improved bounds for mixing rates of Markov chains and multicommodity flow[END_REF] provides simple guidelines for the choice of α. For instance, if the original graph has spectral gap very close to zero and we would like the spectral gap of the modified graph to be approximately equal to 0.1 we choose α = 0.1 d.

Numerical Results

In the following preliminary experiments we use one real-world graph and one random graph. The real-world graph has 5, 204, 176 nodes and 77, 402, 652 edges and was collected in a nearly complete crawl of the Livejournal social blog network [START_REF] Mislove | Measurement and Analysis of Online Social Networks[END_REF]. The Livejournal graph has average degree 14.6 and a giant strongly connected component with 5, 189, 809 nodes and the remaining nodes form a number of small connected components. The random graph is created by connecting, with one edge, two Barabási-Albert graphs [START_REF] Barabási | Emergence of scaling in random networks[END_REF], G 1 and G 2 , with average degrees 2 and 10, respectively. We call the former the BA 2 graph. A RW over the BA 2 graph resembles the RW transition probability matrix described in eq. ( 19), where P i , i = 1, 2, are the transition probability matrices of the RW on each Barabási-Albert graph (G 1 and G 2 , respectively) and ǫ is small. Different from the example shown in Section 2, the average degrees of G 1 and G 2 are different.

Our goal is to compare Random Walks (RWs) against Random Walks with uniform Restarts (RWuRs) in estimating the degree distribution of the graph, i.e., we seek to estimate Θ k , the fraction of nodes with degree greater than k. Let Θk be the estimated value of Θ k .We use

NMSE k = E ( Θk -Θ k ) 2 /Θ k (22)
to measure the estimation accuracy. Parameters: Our experiments have the following parameters. The sampling budget B, which is used in both RW and RWuR. The sampling budget of a RW determines the number of steps. The sampling budget of RWuR does not directly determine the number of steps. This is because there is a sampling budget penalty, c, associated with each restart. For instance, a RWuR that performs m uniform restarts walks ⌊B -mc⌋ steps and gathers ⌊B -mc⌋ + m observations. In all our experiments we use B = n/100 (i.e., the budget is 1% of the total number of nodes in the graph).

Initial RW states: All experiments initialize RW in steady state while RWuR is initialized from uniformly sampled nodes (i.e., RWuR does not start in steady state). This initialization favors RW over our RWuR algorithm. Still, as seen next, our RWuR algorithm outperforms RW in all scenarios (for a given choice of α).

Our first experiment is based on an undirected version of the Livejournal graph. Figure 1(a) shows the empirical NMSE, eq. ( 22), on the Livejournal graph obtained from 20, 000 runs, RWuR restart weights α = 0.01, 10, and RWuR restart penalty c = 10. We choose c = 10 to match the 1/10 hit ratio of MySpace's uniform node sampling [START_REF] Gauvin | Measurement and genderspecific analysis of user publishing characteristics on myspace[END_REF]. We observe that estimates obtained with RWuR α = 0.01 are more precise than the estimates obtained with RW (particularly for small degree nodes and with almost no difference for high degree nodes). Thus, RWuR is able to reduce the NMSE even when restarts are rare (i.e., α is small). We perform the same experiment with increased restart weight α = 10 (Figure 1(a)) and observe that increasing α also increases the accuracy of RWuR for estimating the head of the distribution but decreases the accuracy at estimating its tail than both RW and RWuR with α = 0.01. Note that as we increase α RWuR gets closer to perfoming independent uniform node sampling. We also perform the same experiment over the BA 2 graph (with the only difference that the smallest restart weight is now α = 0.1). Note that the BA 2 graph has a clear RW bottleneck (the edge that connects the two otherwise disconnected components). Figure 1(b) shows the empirical NMSE. Unsurprisingly, the improvement in NMSE obtained by RWuR against RWs is even more pronounced than the improvement observed in our Livejournal experiments. Note that here, similar to the Livejournal experiment, increasing the restart weight from α = 0.1 to α = 10 also increases the accuracy of RWuR for estimating the head of the distribution but decreases its accuracy at estimating the tail.

The above empirical observations, on the relationship between the NMSE of the degree distribution tail and α, prompt us to revisit the analysis performed at the end of Section 1. Putting together eqs. ( 5) and (22) yield

NMSE k ∝ 1 Π k δ , ( 23
)
where Π k is as defined at the end of Section 1. From Section 2 we know that

Π k = ∀v∈V 1(d v = k) k + α 2|E| + nα .
Thus, when k > d ( d is the average degree), Π k decreases with α which implies (eq. ( 23)) that the NMSE increases with α. Similarly, when k < d the NMSE decreases with α. Now let's look at the spectral gap δ. Section 2 shows that δ = 1 -d/(d + α)λ 2 for a d-regular graph and δ ≈ αλ 2 / d for graphs with two quasi-disconnected components that have same average degree, assuming α to be small. These results indicate that δ increases with α.

As one increases α, when k > d the tradeoff between Π k increasing and δ decreasing the NMSE can explain the behavior observed in our numerical results. Similarly when k < d, the NMSE monotonically decreases with α, which can also be observed in our numerical results.

Effect of the restart cost c:

The uniform restarts required in our RWuR algorithm can be expensive. In MySpace [START_REF] Gauvin | Measurement and genderspecific analysis of user publishing characteristics on myspace[END_REF] one needs to perform, in average, 10 queries in order to obtain an user that can be used to restart the RWuR. In the following experiments we compare RW and RWuR with different restart costs c. Figure 2(a) shows the empirical NMSE of RW and RWuR (c = 1, 100; α = 0.01) on the Livejournal graph. The curves in Figure 2(a) are all on top of each other for degrees higher than 120. We notice no loss in NMSE when increasing the restart cost. This is expected as restarts are rare when α = 0.01.

Conclusions and Future Work

Our work proposed and studied the properties of a hybrid sampling scheme (RWuR) that mixes independent uniform node sampling and random walk (RW)-based crawling. Our sampling method combines the strengths of both uniform and RW sampling while minimizing their drawbacks. RWuR can be used in any OSN and P2P networks that allow uniform node sampling (usually at a premium cost), such as MySpace, Facebook, and Bittorrent. We have formally shown under two scenarios that, when compared to a regular RW, RWuR INRIA has larger spectral gap, consequently reducing the mixing time. We also observe that RWuR has a positive impact on reducing the estimation error when compared to regular RWs. As part of our future work we plan to investigate the use of RWuRs to reduce the time to search for a file in a P2P network or a user in an OSN.

We note that d i,k = d k except for the nodes which connect the principal Fiedler clusters. Therefore, we can write

a = 1 - α µ 1 + µ 2 µ 1 |C 1 | 2|E 1 | + µ 2 |C 2 | 2|E 2 | λ 2 + O(ε, |E 1 | -1 , |E 2 | -1 ) = 1 - α d λ 2 + O(ε, |E 1 | -1 , |E 2 | -1 ),
where the last equality follows from the assumption that the principal Fiedler clusters have the same average node degree d = 2|E

1 d k u 2k = k∈C1 1 d k π 1,k + k∈C2 1 d k π 2,k + O(ε) = k∈C1 1 d k d 1,k 2|E 1 | + k∈C2 1 d k d 2,k 2|E 2 | + O(ε) = |C 1 | 2|E 1 | + |C 2 | 2|E 2 | + O(ε, |E 1 | -1 , |E 2 | -1 ) = O(ε, |E 1 | -1 , |E 2 | -1 )
In the latter equality we again use the fact that the principal Fiedler clusters have the same average node degree

d = 2|E 1 |/|C 1 | = 2|E 2 |/|C 2 |.
Thus, only term a provides significant contribution to the final formula [START_REF] Ribeiro | Estimating and sampling graphs with multidimensional random walks[END_REF].

Figure 1 :

 1 Figure 1: Estimation error of RW and RWuR with varying α (larger is worst).

Figure 2

 2 (b) shows the empirical NMSE of RW and RWuR (c = 1, 10, 100, 1000; α = 0.1) on the BA 2 graph. Observe that RWuR outperforms RW when c = 1, 10, 100; the exception happens when c = 1000. Thus, in a graph with a strong RW bottleneck, RWuR is the estimation method of choice even if the cost of restart is high (e.g., c = 100).

Figure 2 :

 2 Figure 2: Estimation error of RWuR with varying restart cost c (larger is worst).

  1 |/|C 1 | = 2|E 2 |/|C 2 |. Next we consider term b. Let us analyze the sum

	n k=1	1 d k u 2k .

n k=1

RR n°7394

Acknowledgements

This research was sponsored in part by the U.S. Army Research Laboratory under Cooperative Agreement Number W911NF-09-2-0053. The views and conclusions contained in this document are those of the authors and do not represent the official policies, either expressed or implied, of the U.S. Army Research Laboratory or the U.S. Government. The U.S Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation hereon. This research was also sponsored in part by the European Commission through the ECODE project (INFSO-ICT-223936) of the European Seventh Framework Programme (FP7).

Appendix: Proof of Theorem 2.3

By construction, the matrix P (0) has two left eigenvectors corresponding to the eigenvalue one u 1,1 = π 1 0 , u 1,2 = 0 π 2 , and two right eigenvectors corresponding to the eigenvalue one

Let us introduce the following matrices

Since the left and right eigenvectors of matrix P (0) form a bi-orthogonal system, we have

Let us also introduce the aggregated transition generator

Next we use the perturbation technique to determine a form of the Fiedler vectors. The perturbed eigenvalue problem can be written as follows:

(P (0) + εC)(v (0) + εv (1) + ...) = (1 -ελ (1) + ...)(v (0) + εv (1)

We set λ (0) = 1, since the eigenvalue one of P (0) splits into two eigenvalues [START_REF] Avrachenkov | Analytic perturbation theory and its applications[END_REF][START_REF] Delebecque | A reduction process for perturbed markov chains[END_REF]. The first perturbed eigenvalue remains equal to one and corresponds to the stationary distribution, the second perturbed eigenvalue corresponds to the Fiedler vectors. We are interested in this second eigenvalue. By equating the terms with the zero power of ε in (26) we obtain

We conclude from the above equation that v (0) is given by some linear combination of the right eigenvectors of P (0) corresponding to the eigenvalue one. Thus, we have

with d ∈ R 2×1 . Equating the terms with the power one of ε we obtain (1) .

Premultiplying the above equation by U 1 and noticing that U 1 (P (0) -I) = 0, we get

RR n°7394

Substituting (27) into the above equation, we obtain the eigenvalue problem for the aggregated generator M d = -λ (1) d.

Denote the entries of the aggregated generator as follows:

The above eigenvalue problem for the aggregated generator has two eigenvalues 0 and µ 1 +µ 2 .

The zero eigenvalue corresponds to the perturbed stationary distribution. Thus we take the second eigenvalue λ (1) = µ 1 + µ 2 corresponding to the Fiedler vectors. As vector d we can take d = [µ 1 -µ 2 ] T which leads to

Similar considerations can be performed for the perturbed eigenvalue problem corresponding to the left Fiedler vector (u (0) + εu (1) + ...)(P (0) + εC) = (1 -ελ (1) + ...)(u (0) + εu (1)

which results in the left eigenvector problem for the aggregated generator eM = -λ (1) e.

and, consequently, with e = [1 -1] we have

Let us now substitute the expressions for the Fiedler vectors (28) and ( 29) into [START_REF] Gkantsidis | Hybrid search schemes for unstructured peer-to-peer networks[END_REF]. Towards this end, we consider separately the two terms in ( 11)

Let us first consider term a.

where the summations are taken over the two principal Fiedler clusters C 1 and C 2 . Next we use the formula for the stationary distribution of the random walk.

INRIA


Carlos Rengifo

Yannick Aoustin

Christine Chevallereau

Franck Plestan

A Penalty-based Approach for Contact Forces Computation in Bipedal Robots

Computation of contact forces is essential for the simulation of mechanical systems with unilateral constraints, like bipedal robots. Most methods are based on the rigid body assumption. They can be categorized into constraint-based and penalty-based approaches. In the former, contact forces are computed by solving an optimization problem based on linear or nonlinear complementarity conditions. Unfortunately, these methods cannot be directly applied to articulated systems described in generalized coordinates. In the second approach, spring-damper models are used to minimize interpenetration between the surfaces in contact. The main criticism to penalty approaches are parameter tunning, static friction handling, and the difficulties to treat multiple simultaneous unilateral contacts. In this work we present a new compliant approach based on input-output feedback linearization. The main advantages of the proposed approach are, the spring-damper parameters are independent of the parameters of the system (i.e masses, inertias), no a priori-defined velocity thresholds are required to distinguish between dynamic and static friction, multiple simultaneous unilateral contacts are naturally handled. The proposition has been succesfully applied to the simulation of a 3D bipedal walking robot.

I. INTRODUCTION

Bipedal robots are mechanical systems with contacts intermittence, unilateral constraints and impacts. These systems are traditionally described as hybrid dynamical systems and by consequence biped walking simulators commonly uses an event-driven approach: impacts are regarded as the events producing transitions between biped phases. In the case of bipedal robots with point-feet there are only four possible phases: right single support, left single support, double support and fly. In robots with no trivial feet the number of states and events dramatically increases. For example, in a rectangular footed biped, each foot may be contacting in one of 10 ways (at one of four vertices, one of four edges, with its face, or not in contact at all) and the foot can impact the ground in 9 different ways. Thus, the resulting system is a finite state machine with 100 states and 81 possibles events. In [START_REF] Hale | Robot simulation, collisions and contacts[END_REF], as contacts with static friction are handled separately from those with dynamic friction, their software platform for humanoid robot simulation considers 361 possible contact states. The contact enumeration problem is discussed in [START_REF] Yamane | Dynamic filter -concept and implementation of online motion generator for human figures[END_REF]. Given the exploding complexity of the event-driven approach as the number of contact points increases, techniques stemmed from computer graphics animations are currently applied to humanoid robot simulators. These techniques can be categorized into constraint-based and penalty-based approaches. In the former, contact forces are computed by solving an optimization problem based on linear or nonlinear complementarity conditions. The most natural way to formulate these conditions is in terms of accelerations and forces. However it leads to problems without solution in the classical sense. Even with the use of impulsive forces, the existence of solutions can not be guaranteed [START_REF]Rigid-body dynamics with friction and impact[END_REF]. Instead of acceleration and forces, complementarity conditions can be formulated in terms of velocities and time integrals of the forces, provided the transformation of the differential equations of the system into finite differences. Many variants of this approach, called time-stepping, have been proposed. Anitescu [START_REF] Anitescu | A time-stepping method for stiff multibody dynamics with contact and friction[END_REF] and Stewart [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with coulomb friction[END_REF] use a linear complementarity formulation requiring the approximation of the friction cone by a polyhedron. The advantage of Anitescu and Stewart methods is the existence and uniqueness of solutions can be guaranteed [START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of painlevé's problem[END_REF]. Unfortunately, these methods are formulated using maximal coordinates and by consequence cannot be directly applied to articulated systems described in generalized coordinates. When maximal coordinates are used, joints must be treated as bilateral constraints. The major disadvantage of this is interpenetration due to numerical drift. Interpenetration must be corrected using post stabilization methods like [START_REF] Cline | Post-stabilization for rigid body simulation with contact and constraints[END_REF]. The proposition of Liu and Wang [START_REF] Liu | Computation of three-dimensional rigid-body dynamics with multiple unilateral contacts using time-stepping and Gauss-Seidel methods[END_REF] uses maximal coordinates also, but it does not require some friction cone approximation. However, bilateral constraints are not considered, so its applicability to robotic systems is very limited. All timestepping approaches have the advantage that collisions and contact forces are computed in a unified way. On the other hand, they have many limitations from the computational point of view as stated in [START_REF] Drumwright | A fast and stable penalty method for rigid body simulation[END_REF], and they do not have been extended to articulated bodies described with generalized coordinates. In [START_REF] Chardonnet | Dynamic simulator for humanoids using constraint-based method with static friction[END_REF] a simulator for the robot HRP2 using a combination of time-stepping and the acceleration-force formulation proposed in [START_REF] Ruspini | Collision/contact models for the dynamic simulation of complex environments[END_REF] is presented. In the penalty based approach spring-damper models are used to minimize interpenetration between the surfaces in contact [START_REF]Stable penalty-based model of frictional contacts[END_REF]. Penalty methods are used in most biped simulators. Some recent examples are Saika-3 [START_REF] Hwang | An order n dynamic simulator for a humanoid robot with a virtual spring-damper contact model[END_REF], JOHNNIE [START_REF] Buschmann | Dynamics simulation for a biped robot: modeling and experimental verification[END_REF] and OpenHRP [START_REF] Hirukawa | Experimental evaluation of the dynamic simulation of biped walking of humanoid robots[END_REF]. In these biped simulators spring-damper models are used to compute normal and tangential contact forces. In [START_REF] Hwang | An order n dynamic simulator for a humanoid robot with a virtual spring-damper contact model[END_REF], a classical spring-damper model is used, the problem is that sticking normal forces appear just before contact separation. In [START_REF] Buschmann | Dynamics simulation for a biped robot: modeling and experimental verification[END_REF], a Heaviside function is used to prevent this problem. In [START_REF] Hirukawa | Experimental evaluation of the dynamic simulation of biped walking of humanoid robots[END_REF] the computation of normal contact forces is based on the nonlinear spring-damper model proposed in [START_REF] Marhefka | Simulation of contact using a nonlinear damping model[END_REF]. This model assures continuous normal contact forces. The problem is that sticking forces arise when an external force separates the objects at high velocity. That is the case in biped walking simulations when the stance foot takes off after the impact. The three above mentioned methods do not require velocity thresholds to distinguish between static and dynamic friction. The main criticisms to these penalty-based approaches are parameter tunning, and the difficulties to treat multiple simultaneous unilateral contacts. In this paper a new penalty-based approach based on inputoutput feedback linearization is proposed. The main advantages of the proposed approach are: the parameters of the contact model are independent of those of the system (i.e masses, inertias), no a priori-defined velocity threshold are required to distinguish between dynamic and static friction and multiple simultaneous unilateral contacts are naturally handled. This paper is organized as follows. Section II describes the proposed methodology for the contact force computation using two simple examples. In section III the methodology is applied to a bipedal robot with rectangular feet. In Section V, the numerical problem of contact stability is introduced. In Section IV simulation results are presented. The last section of this paper is devoted to conclusions and perspectives.

II. THE PROPOSED APPROACH

In order to illustrate the general idea of the proposed methodology, two simple examples are considered. The first one is devoted to the normal force computation, and the second one to friction forces. The incidence of the contact model parameters on the resulting reaction forces will be also illustrated.

A. Bouncing ball

Consider a bouncing ball of mass m as illustrated in Figure 1. The differential equation for the position is:

mp z = -mg + f n (1)
The normal contact force f n is zero when the position p z is positive, and f n is greater than zero when p z = 0. If we want to simulate the system (1) then raises the question how to compute the reaction force between the ball and the ground?.

If a rigid impact with restitution coefficient zero is supposed, then the resulting f n has two components as illustrated in Figure 2. A Dirac's delta function of weight -m ṗz t - i ; and constant value mg for t > t i . The impulse, which occurs at the impact instant t = t i , is required to discontinuously change the velocity from a negative value to zero. The constant component of f n guarantees pz (t > t i ) = 0. The function f n of the Figure 2 assures zero velocity and acceleration after the impact. This example clearly shows the In the penalty methods used in the simulators for the bipedal

fn t t = t i mg -m • ṗz t - i • δ (t -t i )
Fig. 2. Theoretical reaction force for the bouncing ball example when a rigid impact with zero restitution coefficient is supposed.

robots Saika-3 [START_REF] Hwang | An order n dynamic simulator for a humanoid robot with a virtual spring-damper contact model[END_REF], Johnnie [START_REF] Buschmann | Dynamics simulation for a biped robot: modeling and experimental verification[END_REF] and OpenHRP [START_REF] Hirukawa | Experimental evaluation of the dynamic simulation of biped walking of humanoid robots[END_REF], impact and resting contact forces are computed in a unified way. When a typical penalty-based method like [START_REF] Hwang | An order n dynamic simulator for a humanoid robot with a virtual spring-damper contact model[END_REF], is applied to the bouncing ball system, f n is computed using a springdamper model:

f n =    0 p z > 0 -k p p z -k v ṗz p z ≤ 0 (2)
As noted in [START_REF] Marhefka | Simulation of contact using a nonlinear damping model[END_REF], the main inconvenient of expressions like [START_REF] Buschmann | Dynamics simulation for a biped robot: modeling and experimental verification[END_REF], is that before contact separation forces tending to stick the ball to the ground always appear . The expressions used in the biped simulators described in [START_REF] Hirukawa | Experimental evaluation of the dynamic simulation of biped walking of humanoid robots[END_REF] and [START_REF] Buschmann | Dynamics simulation for a biped robot: modeling and experimental verification[END_REF], avoid this problem using nonlinear functions.

In the approach proposed in this paper, when p z ≤ 0, we consider that the acceleration of the ball is imposed by the contact model. This acceleration is opposed to interpenetration

pd z = -k p p z -k v min (0, ṗz ) (3)
Hereafter pd z will be called the desired acceleration. If k p > 0 and k v > 0, the term min (0, ṗz ) assures a desired positive acceleration for all p z < 0. Once we have pd z , the following step consists to compute the required f n to obtain pz = pd z . For the model (1) such force is

f n = mg + mp d z .
f n in the general case is given by

f n =    0 p z > 0 mg + mp d z p z ≤ 0 (4)
The procedure for the computation of f n when p z ≤ 0 is summarized in Figure 3: first, the desired acceleration is computed using the equation ( 3) and then the corresponding contact force is computed using (4). In the general case a saturation function is applied to assure f n ≥ 0. In this simple case pd z ≥ 0 implies f n ≥ 0. As the equation ( 3) can be rewritten as

pd z = k p p d z -p z + k v max 0, ṗd z -ṗz f n
is the required force to obtain a desired position p d z = 0 and a desired velocity ṗd z = 0. A fundamental difference between (2) and ( 4) is steady state interpenetration depth. In the former case, this value is equal to -g.m/k p and for the second expression is zero. However, in the Section V the case of a nonzero desired interpenetration is discussed. Given ṗz (t i ) < 0 and p z (t i ) = 0, the velocity and the position of the ball at the impact instant t = t i , the dynamical behaviour in contact is determined by the following differential equation:

pz = -k p p z -k v min (0, ṗz ) , t ≥ t i
Parameter selection for the contact model concerns only the equation above and does not depend on the mass of the ball.

In more complex systems the proposed procedure allows a parameters selection independently of the mathematical model of the systems and the number of points in contact with the ground. In

B. Brick on a frictional ramp

Consider a brick sliding down a ramp as illustrated in Figure 6. The differential equation for the position is:

m px = mg sin α -f r (5)
Where f r is the frictional force. If the initial velocity is strictly positive ( ṗx (0) > 0), and the frictional force is greater than the gravitational one (µ cos α > sin α), then the velocity of the brick goes to zero and the brick will be stopped. In such case, if Coulomb's friction model is applied, f r is a discontinuous function of the velocity: Now we will apply a procedure similar to that of the last example. The desired tangential acceleration is chosen as:

f r =    µmg cos α sign ( ṗx ) | ṗx | > 0 mg sin α ṗx = 0 (6
pd x = -k v ṗx (7)
The required friction force f r to obtain px = pd x from ( 5) is given by f r = mg sin α -mp d

x . However, the magnitude of the frictional force is limited by the normal force and the friction coefficient. Thus, a saturation function must be applied in order to assure an f r such that |f r | ≤ µf n . The normal force exerted by the ramp on the brick is

f n = mg cos α f r = sign (f ′ r ) • min (|f ′ r |, µmg cos α) (8)
with f ′ r = mg sin α + mk v ṗx The equation ( 8) can be rewritten in a more classical form:

f r =    µmg cos α • sign ( ṗx ) | ṗx | > v s mg sin α + mk v ṗx | ṗx | ≤ v s (9)
With 5) becomes px = -k v ṗx . In such case, the velocity and the force respectively converge to zero and to mg sin α.

v s = (g/k v )(µ cos α -sin α). When | ṗx | ≤ v s the model (
The procedure for the computation of f r is summarized in Figure 7: first, the desired acceleration is computed using the equation ( 7) and then the corresponding contact force is computed using [START_REF] Liu | Computation of three-dimensional rigid-body dynamics with multiple unilateral contacts using time-stepping and Gauss-Seidel methods[END_REF]; second a saturation function is applied to this force in order to assure |f r | ≤ µf n . As it can be seen from the Figure 7, f r is the required force to obtain a desired velocity ṗd x = 0, given the constraint |f r | ≤ µf n .

force acceleration

Desired a ramp

Brick on Friction

III. COMPUTATION OF CONTACT FORCES IN A BIPEDAL ROBOT USING INPUT-OUTPUT FEEDBACK LINEARIZATION

The bipedal robot used in this paper is depicted in Figure 8. It consists of five rigid links and two rectangular feet connected by fourteen motorized joints to form a serial structure. Eight possible contact points are considered, one for each vertice of the foot. Impacts on edges are treated as a simultaneous collision with the ground of two vertices. Flat foot impacts are treated as a simultaneous collision with the ground of four vertices. The mathematical model of the biped is given by: A (q) q a + H (q p , q v ) = Γ e + D (q p ) F,

With

q a =   0 V0 0 ω0 q   , q v =   0 V 0 0 ω 0 q   Inertial frame 3 DoF 3 DoF 3 DoF 1 DoF 1 DoF Zero frame 3 DoF v 2 v 4 v 3 v 1 v 5 v 6 v 7 v 8
p =   X 0 α q   , Γ e =   0 0 Γ  
where X 0 and α are respectively the position and orientation of the right foot in the inertial frame. q ∈ IR 14 are joint positions. 0 V 0 ∈ IR 3 and 0 ω 0 ∈ IR 3 are respectively the linear and angular velocity of the zero-frame. 0 V0 ∈ IR 3 and 0 ω0 ∈ IR 3 are respectively the linear and angular acceleration of the zero-frame. These velocities and accelerations are defined as:

0 V 0 ≡ 0 R I • I V 0 0 ω 0 ≡ 0 R I • I V 0 (11) 0 V0 ≡ 0 R I • d I V 0 dt 0 ω0 ≡ 0 R I • d I ω 0 dt
where I V 0 ∈ IR 3 and I ω 0 ∈ IR 3 are respectively the linear and angular velocity of the zero-frame with respect to the inertial frame. The rotation matrix 0 R I ∈ IR 3×3 expresses the orientation of the inertial frame with respect to the zero-frame. As the right foot moves, the matrix 0 R I is not constant. It can be shown that [START_REF] Ruspini | Collision/contact models for the dynamic simulation of complex environments[END_REF] implies:

d 0 V0 dt = 0 V0 -0 ω 0 × 0 V 0 d 0 ω0 dt = 0 ω0 ,
as consequence the time derivative of q v is not q a . Concerning the equation ( 10), A ∈ IR 20×20 is the inertia matrix, H ∈ IR 20 is the vector including Coriolis, centrifugal and gravitational terms, Γ ∈ IR 14 is the control input, D ∈ IR 20×24 , is the transposed of the Jacobian matrix relating the Cartesian velocities of the vertices of the foot with the joint velocities. The vector F = [f T 1 . . . f T 8 ] T ∈ IR 24 contains all contact force vectors. Each contact force f i is a three dimensional vector containing the tangential forces f it ∈ IR 2 and the normal force

f n ∈ IR f i = f it f in , f it = f ix f iy .
Contact forces are subject to the following constraints:

f in ≥ 0, |f it | ≤ µf in i = 1 . . . 8 (12)
Every contact point has an associated scalar constraint φ i (q, α) ≥ 0. Being φ i (q, α) the distance between the vertice i of the foot and the ground. When φ i (q, α) > 0, the vertice is not in contact with the ground and the corresponding f i is zero. Conversely, when φ i (q, α) = 0, the magnitude of the corresponding f i is not zero. Given n c (t k ), the set of contact points satisfying the condition φ i (q, α) = 0 at the instant time t = t k , the term D (q, α) F of the equation ( 10) at such instant can be written in terms of the nonzero contact forces

D (q, α) F = i∈nc D i (q, α) f i ( 13
)
D i ∈ IR 20×3 is composed by the three columns of D corresponding to the forces f in , f ix and f iy . For example, if at a given instant the vertices v 1 and v 4 (see Figure 8) are contact with the ground, then n c = {1, 4}. In such case D 1 is composed of the columns of the matrix D corresponding to the forces f 1n , f 1x , f 1y , and D 4 of those corresponding to f 4n , f 4x , f 4y . The equation ( 13) can be rewritten in a compact form as:

D (q, α) F = D c (q, α) F c ( 14
)
D c is composed by the matrix D i , i ∈ n c , and F c by the vectors f i , i ∈ n c . Now consider P = [p T 1 . . . p T 8 ] T ∈ IR 24 , the vector containing the Cartesian position of the vertices of the foot in the inertial frame. With p i (q) = [p ix , p iy , p iz ] T ∈ IR 3 . The relationship between Ṗ and q v is given by the Jacobian matrix J = D T :

Ṗ = D T q v
Taking the derivative of the kinematic equation above, the Cartesian accelerations P can be exprimed as a linear function of F ,

P = D T q a + B = D T A -1 D • F + D T A -1 (Γ e -H) + B (15)
with

B ≡ -D T   0 ω 0 × 0 V 0 0 0   + ḊT q v
If instead of the eight Cartesian positions, we consider only the Cartesian positions of the vertices in contact with the ground (P c (q) = [ p T i ] T , i ∈ n c ), the equation ( 15) becomes:

Pc = D T c A -1 D c • F c + D T c A -1 (Γ e -H) + B c (16) with B c ≡ -D T c   0 ω 0 × 0 V 0 0 0   + ḊT c q v
The main point of our approach is to consider that Pc is imposed by the contact model and then solve for F c from [START_REF]Stable penalty-based model of frictional contacts[END_REF]. The accelerations given by the contact model will be denoted by P d . The component pd i of the vector P d is the required acceleration to asymptotically stops the motion of the contact point i ( ṗi → 0) and removes its interpenetration (p iz → 0)

pd i = -k p   0 0 p iz   -k v   ṗix ṗiy min (0, ṗiz )   , i ∈ n c (17)
Once the desired accelerations are defined, the following step is to compute the required F ′ c to obtain

P = P d D T c A -1 D c • F ′ c = P d -D T c A -1 (Γ e -H) + B c (18)
The system of linear equations (18) may have an infinity of solutions. For example, in a bipedal robot, when the four vertices of the foot are in contact with the ground during a single support phase, the contact forces cannot be uniquely determined. In such cases, the solution with minimum Euclidean norm is computed using the Moore-Penrose pseudo-inverse. Once F ′ c is computed, the following step consists to build the vector F ′ composed of F ′ c and zeros. After that, saturations are applied to F ′ in order to obtain contact force vector F satisfying [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with coulomb friction[END_REF]. The following saturation function assures a nonnegative normal contact force:

f in = max 0, f ′ in (19)
As the maximal value of |f it | is µf in , the maximal values for |f ix | and |f iy | are respectively µf in cos θ and µf in sin θ.

With θ = atan f ′ iy , f ′ ix . Thus, the components of the tangential force vector are computed using the procedure depicted in Figure [START_REF] Liu | Computation of three-dimensional rigid-body dynamics with multiple unilateral contacts using time-stepping and Gauss-Seidel methods[END_REF]. This procedure can be exprimed mathematically as:

f ix = sign f ′ ix • min |f ′ ix |, µf in cos θ f iy = sign f ′ iy • min |f ′ iy |, µf in sin θ (20)
As it can be seen from this equation, the transition between dynamic and static friction does not produce discontinuous tangential forces and it does not require a priori defined velocity thresholds. If F ′ c satisfies the physical constraints [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with coulomb friction[END_REF], then Pc = P d c In such case, contact forces and accelerations are decoupled. In this approach, inspired from input-output feedback linearization, parameter selection for the contact P.D Direct Geometric Model Direct Kinematic Model

P d = 0 Ṗ d = 0 P d qv, qp F ′ F Pc, Ṗc
Robot Dynamics Linearisation and Decoupling Saturation Fig. 10. Contact force computation for the bipedal robot described by [START_REF] Marhefka | Simulation of contact using a nonlinear damping model[END_REF].

model is independent on the number of points in contact with the ground and the parameters of the robot (i.e mass, inertias). The contact force computation for the system ( 10) is described in Figure 10.

IV. NUMERICAL TESTS

The procedure described in this paper for contact force computation has been applied to the simulation of a walking gait consisting of five steps. The selected parameters for the Equation (17) were k p = 500 and k v = 250. In Figure 11 the distances between the vertices of the foot and the ground are traced as a function of the time. In Figure 12, these distances are traced when the right foot is in contact with the ground after its first impact. Figure 13 shows the corresponding normal contact forces. V. CONTACT STABILITY In this section the numerical stability of the contact between the foot and the ground will be introduced. It has important consequences on normal forces computation. In the absence of external forces, the contact between an object and the ground is said in equilibrium if after a given instant T c , the normal force exerted by the ground is equal to the weight of the object. In such case, the application of a penaltybased approach leads to p z = ṗz = pz = 0, ∀t > T c . If an external force is added to the normal force exerted by the body on the ground, the equilibrium is broken and interpenetration immediately appears. However, after a transient time, position, velocity and acceleration return to zero and the normal force converges toward a new steady state value. If the external force mentioned above disappears, the contact is broken and a positive distance appears. As consequence, the normal contact force becomes zero. As no opposition to the gravity force exists, the object will be again in contact with the ground after a short instant generating a impact. These impacts are not physical ones, they occur only in simulation. For this reason, thereafter they will be called false impacts. This phenomena occurs in bipedal robot simulation. During the single support phases, the foot in contact with the ground is subject to a time varying forces generated by the motion of the other bodies. As we can see from Figure 13, false impacts produce discontinuities (vertical lines) in the computed normal forces. Contact stability can be improved if steady state interpenetration is allowed. Instead of converge towards p iz = 0, the Equation (17) can be modified for obtaining convergence towards p iz = ǫ.

VI. CONCLUSION

In this paper we have presented a new penalty approach based on input-output feedback linearization. This idea coming from control theory, allows to treat multiple simultaneous contacts independently. Another is parameters of the contact model are independent of the parameters of the system. The computation of contact forces is composed of two steps. First, the required forces to avoid interpenetration and stop the motion of the contact points are computed. Second, nonlinearities are applied to these forces in order to satisfy physical constraints. The proposition was successfully applied to the simulation of a 3D bipedal walking robot. A simple solution was proposed to improve the numerical stability of the contacts between the foot and the ground.

Fig. 1 .

 1 Fig. 1. Bouncing ball

Fig. 3 .

 3 Fig. 3. Bouncing ball. Computation of the normal force when pz ≤ 0.

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. Bouncing ball vertical position. Height (mm) versus time (sec): kp = 0, kv = 5000 (dashed line), kp = 1, kv = 5000 (solid line), kp = 5, kv = 5000 (dotted line).

Fig. 7 .

 7 Fig. 7. Brick on frictional ramp. Computation of the friction force.

Fig. 8 .

 8 Fig. 8. General diagram of the bipedal robot and notation for the vertices.

  q

Fig. 9 .

 9 Fig. 9. Nonlinearities used to get tangential contact forces inside of the friction cone.

Vertices v 3 , v 7 (Vertices v 4 , v 8 (Fig. 11 .

 374811 Fig. 11. Vertical position of the feet. Height (m) versus time (sec) for a walking gait consisting of four steps. Right foot (solid line), Left foot (dotted line)

Vertice v 4 (Fig. 12 .

 412 Fig. 12. Vertical position of the feet around zero height. Height (mm) versus time (sec) for the vertices of the right foot after its first impact.

Vertice v 4 (Fig. 13 .

 413 Fig.[START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of painlevé's problem[END_REF]. Normal contact forces. Normal contact forces (N ) versus time (sec). These forces correspond to the simulations results presented in the Figure12.

  The normal contact forces corresponding to ǫ = 1 × 10 -6 m and the same values of k p and k v (k p = 500, k v = 250) are showed in Figure14.

Vertice v 4 (Fig. 14 .

 414 Fig.[START_REF]Rigid-body dynamics with friction and impact[END_REF]. Normal contact forces. Normal contact forces (N ) versus time (sec) when steady-state interpenetration is allowed. These results must be compared of those presented in Figure13.

  )

		2
		1.5
		1
	Height [millimeters]	0.5
		0
		-0.5
		0.4 -1	0.6	0.8	1	1.2	1.4	1.6	1.8	2	2.2	2.4
							Time [sec]