Shenhong Li

Mahsa Derakhshani
email: m.derakhshani@lboro.ac.uk

Chung Shue Chen

Chung Shue

Chen

Sangarapillai Lambotharan
email: s.lambotharan@lboro.ac.uk

Outage Probability Analysis for Two-antennas MISO-NOMA Downlink with Statistical CSI

Keywords: NOMA, MISO, Outage probability

In this paper, we analyze the outage probability of the multi-user multiple-input single-output (MISO) downlink system by combining the non-orthogonal multiple access (NOMA) scheme. We derive tractable closed-form outage expressions given a minimum target rate for the individual users for the case of two antennas, by modeling cumulative distribution function (CDF) of received signal-to interference plus noise ratio (SINR). Simulation results illustrate the outage performance for different power allocation scenarios and verify the accuracy of our outage probability analysis.

I. INTRODUCTION

Currently, non-orthogonal multiple access (NOMA) has received significant attention for future generation networks that can achieve low latency, massive connectivity, high reliability, and high throughput [START_REF] Dai | Non-orthogonal multiple access for 5G: Solutions, challenges, opportunities, and future research trends[END_REF]. Multiple-input and singleoutput (MISO) [START_REF] Liu | Multiple-antenna-assisted non-orthogonal multiple access[END_REF] is one of the smart antenna technologies for wireless communications in which antennas are used at the transmitter side, and the received side only has one antenna. From the application viewpoint, the combination of NOMA and multiple antenna techniques is effective to enhance system performance. The essential intention of NOMA scheme in this work is to simultaneously serve multiple users over the same radio resources by applying successive interference cancellation (SIC) at the receiver side to cancel the inter-user interference.

Zero-forcing beamforming (ZFBF) is a practical multi-user transmission strategy. By designing one user's beamforming weight vector to be orthogonal to other selected user's channel vectors, ZFBF can completely eliminate multiuser interference [START_REF] Kimy | Non-orthogonal multiple access in a downlink multiuser beamforming system[END_REF]. Unlike conventional multi-user beamforming, NOMAbased beamforming involves designing a single beamforming vector to support multiple users. Hence, analyzing the MISO-NOMA systems becomes a prominent research trend. In terms of system performance, the authors in [START_REF] Liu | On the capacity comparison between MIMO-NOMA and MIMO-OMA[END_REF] demonstrated that the sum rate of MIMO-NOMA system is greater than that of MIMO-OMA system. Moreover, in order to further improve the system performance, the authors in [START_REF] Hanif | A minorization-maximization method for optimizing sum rate in the downlink of non-orthogonal multiple access systems[END_REF] studied the sum rate maximization problem for MISO-NOMA downlink system and proposed an iterative algorithm to solve a secondorder cone program (SOCP) with polynomial computational complexity. [START_REF] Liu | Joint beamforming and power optimization with iterative user clustering for MISO-NOMA systems[END_REF] examined iterative user clustering, and joint optimization with beamforming design and power allocation for MISO-NOMA downlink systems. The authors in [START_REF] Sun | Precoder design in statistical CSI aided non-orthogonal multiple access[END_REF] investigated the ergodic rate maximization problem for MISO-NOMA system considering statistical CSI under a constraint on the total transmission power. On the other hand, the authors in [START_REF] Mao | Power minimization resource allocation for underlay MISO-NOMA SWIPT systems[END_REF] proposed a minimum total power problem in MISO-NOMA system. A dynamic clustering optimization problem from a fairness perspective was investigated in [START_REF] Liu | Fairness of user clustering in MIMO non-orthogonal multiple access systems[END_REF] for MIMO-NOMA systems. Since outage probability is another evaluation criterion for system performance, it attracts many researchers to study it.

In conventional NOMA systems, the outage probability is widely used in the existing works [START_REF] Liu | Outage analysis and power allocation in uplink non-orthogonal multiple access systems[END_REF], [START_REF] Li | Outage-constrained robust power allocation for downlink MC-NOMA with imperfect SIC[END_REF] and [START_REF] Cui | Rethinking outage constraints for resource management in NOMA networks[END_REF] for performance analyses when the channel is characterized by random parameters. Furthermore, for MIMO-NOMA systems, the rate outage probability has been analyzed in [START_REF] Cui | Outage probability constrained MIMO-NOMA designs under imperfect CSI[END_REF] considering imperfect scenario, where channel gains are modelled as the estimated values plus an estimation error. Moreover, the outage probability expression in the case of channel estimation error is approximated by applying central limit theorem (CLT) in the case of a large number of antennas.

In this paper, the downlink of single cell MISO-NOMA system with one multi-antenna base station serving singleantenna users is considered where users are grouped to multiple clusters. NOMA is implemented for the users within each cluster, while ZFBF is used to cancel inter-cluster interference. We also assume that maximum two users can be in one cluster, to keep the SIC complexity low, which is a common assumption in the NOMA literature. We concentrate on the analytical study of the behavior of outage probability for individual users. Specifically, we derive a tractable closedform outage expression, which is defined as the probability that the instantaneous data rate is below a target value, for the case of two antennas. The outage probability of the users with more favorable channel conditions (called near users) can be obtained based on the cumulative distribution function (CDF) of exponential distribution for Rayleigh channels. On the other hand, the analysis of outage probability of the users with less favorable channel conditions (far users) is more challenging.

We derive approximated closed form expressions for far users based on the CDF of the ratio between two exponentially distributed random variables. This work only focuses on the case of two antennas at the BS since it is the first step to provide insights on rate outage analysis for MISO-NOMA systems under statistical CSI and as far as we know it has not been investigated in the MISO-NOMA literature previously. This analysis will be extended for a general case with N antennas in the future works.

The rest of this paper is organized as follows. Section II outlines the system model. Section III introduces the techniques of beamforming design and NOMA scheme. In Section IV, we formulate the received signal to interference plus noise ratio (SINR) of each user. Then, we analyse the outage probability of near and far users in Section V. Simulation results are shown in Section VI, followed by conclusions in Section VII.

We use boldface capital and lower case letters to denote matrices and vectors, respectively. (.) T to denote the transpose, E(.) to denotes the expectation.

II. SYSTEM MODEL

Consider a multiuser downlink MISO system consisting of one base station (BS) with two antennas, serving two clusters of users, each equipped with single receive antennas. The belonging users to one cluster share a beamforming, and thus, NOMA is implemented as the multiple access technique in each cluster. As shown in Fig. 1, we further assume that each cluster has only two users to keep the SIC complexity low. Based on the distance from BS, the two users in each cluster are called either near and far users. In practice, it is difficult to acquire the accurate channel state information due to channel estimation and quantization. In this work, we assume that only the statistical knowledge of CSI is available, hence the channel (row) vector for each user at cluster n is modeled as a random vector

h i,n = [h m i,n ] ∈ C 1×2
, where i represents the index for users, n denotes the index for clusters, and m represents the index for antennas. All of them can be equal to either 1 or 2. Assuming all wireless channels in the network exhibit Rayleigh block flat fading characteristics, the entries in h i,n can be modeled as independent and identically distributed (i.i.d.) complex Gaussian random variables with zero-mean and unit variance.

We adopt linear transmit precoding at the BS, where each cluster is assigned with one dedicated information beam. The weights for all antennas for the near user at cluster n can be collected into a 2×1 beamforming vector w n = [w 1,n , w 2,n ] T , where (.) T denotes the transpose of vector. Accordingly, the received signals in the cluster n for two users are given by

y i,n = h i,n √ d -λ i,n 2 ∑ j=1 √ P j w j x j + e i,n , (1)
(i = 1, 2; n = 1, 2)
where the transmitted signal x j is expressed as the scalar value

x j = √ α 1,j s 1,j + √ α 2,j s 2,j
, in which α 1,j and α 2,j are the transmit power allocation coefficients for two users in cluster j such that α 1,j + α 2,j = 1; s 1,j and s 2,j are the message signals for both users in cluster j; e i,n ∼ CN (0, σ 2 ) denotes the i.i.d. complex Gaussian noise value. Moreover P j is the transmit power for cluster j, d i,n is the distance of user i in cluster n from the BS, and λ is the path loss exponent.

III. PROPOSED TECHNIQUES A. Beamforming Design

Transmit beamforming uses multiple antennas and channel state information at the transmitter (CSIT) to steer radiated power towards directions of interest, limiting leakage in other directions. The traditional approach for multi user MIMO processing is to derive the sum capacity by dirty paper coding (DPC) [START_REF] Jindal | Dirty-paper coding versus TDMA for MIMO broadcast channels[END_REF]. However, such coding scheme is difficult to implement in practical system because of brute-forcing method. The simple zero-forcing method is utilized in this work to cancel the inter-cluster interference. Unlike conventional multi-user beamforming, NOMA-based multiuser beamforming involves designing a single beamforming vector to support multiple users in a cluster. In this system, the channel gain vectors of users are proposed as

h 1,1 = [h 1 1,1 , h 2 1,1 ], h 2,1 = [h 1 2,1 , h 2 2,1 ], (2)
h 1,2 = [h 1 1,2 , h 2 1,2 ], h 2,2 = [h 1 2,2 , h 2 2,2 ]
, where the superscript of h 1 1,1 is the index of antenna, the subscript {1, 1} represents the near user located at cluster 1, and the subscript {2, 1} represents the far user at cluster 1. Thus, the beamforming vector w n for cluster n is chosen to be orthogonal to all other clusters' channel vector h 1,n ′ , satisfying the following condition

h 1,n ′ |h 1,n ′ | w n = { 0, for n ′ ̸ = n 1, for n ′ = n. ( 3
)
The beamforming matrix in such MISO-NOMA system is designed to cancel out the inter-cluster interference for near users of all clusters, so that they can effectively perform SIC using NOMA.

When the normalized channel gain matrix for the near user is represented as

H 1 =   h 1 1,1 |h1,1| h 2 1,1 |h1,1| h 1 1,2 |h1,2| h 2 1,2 |h1,2|   , ( 4
)
where each antenna of channel gain follows the Rayleigh fading, the corresponding zero-beamforming matrix of the square channel gain matrix becomes

W = H 1 -1 = |h 1,1 ||h 1,2 | h 1 1,1 h 2 1,2 -h 1 1,2 h 2 1,1 [ h 2 1,2 |h 1,1 | -h 2 1,1 |h 1,2 | -h 1 1,2 |h 1,1 | h 1 1,1 |h 1,2 |
] .

B. NOMA scheme

In this work, NOMA is proposed to cancel the inter-user interference in each cluster. The SIC is utilized at the received side to separate signals in each cluster. However, it is hard to determine the decoding order since we only know the statistical knowledge of distribution of CSI at the transmitter. The SIC receiver first decodes the near user's signal by treating the far user's signal as the interference. After this detection stage, the near user's signal is removed from the composition signal. The far user treats the signal of the near user as interference without cancellation.

IV. RECEIVED SINR

In this section, we derive the signal to interference plus noise ratio (SINR) of the near and far users for the proposed MISO-NOMA system. The received signal of the near user in cluster 1 can be formulated as

y 1,1 = h 1,1 √ d -λ 1,1 P 1 w 1 x 1 + h 1,1 √ d -λ 1,1 P 2 w 2 x 2 + e 1,1 (5)
= h 1,1 √ d -λ 1,1 P 1 w 1 ( √ α 1,1 s 1,1 + √ α 2,1 s 2,1 )
+ h 1,1 √ d -λ 1,1 P 2 w 2 x 2 + e 1,1 ,
where h 1,1 √ d -λ 1,1 P 2 w 2 x 2 is the inter-cluster interference from the cluster 2, which can be removed by using ZF beamforming method and h 1,1

√ d -λ 1,1 P 1 w 1 √ α 2,1 s 2,1
is the inter-user interference from the far user in the same cluster, which will be canceled out with perfect SIC. Thus, y 1,1 will be reduced to

y 1,1 = h 1,1 √ d -λ 1,1 P 1 w 1 √ α 1,1 s 1,1 + e 1,1 . (6)
The received SINR of the near user in cluster 1 is then

SINR 1,1 = h 1,1 √ d -λ 1,1 P 1 w 1 √ α 1,1 2 σ 2 = |h 1,1 | 2 d -λ 1,1 α 1,1 P 1 σ 2 = χ 1,1 α 1,1 d -λ 1,1 P 1 σ 2 , ( 7
)
where χ 1,1 follows a gamma distribution with Gamma(2, 1).

Similarly, the received SINR of the near user in the cluster 2

will be SINR 1,2 = χ1,2α1,2d -λ 1,2 P2 σ 2 .
On the other hand, the received signal of the far user in cluster 1 is

y 2,1 = h 2,1 √ d -λ 2,1 P 1 w 1 x 1 + h 2,1 √ d -λ 2,1 P 2 w 2 x 2 + e 2,1 (8)
= h 2,1 √ d -λ 2,1 P 1 w 1 ( √ α 1,1 s 1,1 + √ α 2,1 s 2,1 )
+ h 2,1 √ d -λ 2,1 P 2 w 2 x 2 + e 2,1 .
Unlike the near user, the ZFBF vector does not satisfy (3) for the far user in cluster n since we only have one beamforming vector, and the far user does not perform SIC. Therefore, the received SINR of far user in cluster 1 is

SINR 2,1 = |h 2,1 w 1 | 2 d -λ 2,1 α 2,1 P 1 |h 2,1 w 1 | 2 d -λ 2,1 α 1,1 P 1 + |h 2,1 w 2 | 2 d -λ 2,1 P 2 + σ 2 . ( 9
)
The expression of received SINR of the far user in cluster 2 also can be written as

SINR 2,2 = |h 2,2 w 2 | 2 d -λ 2,2 α 2,2 P 2 |h 2,2 w 2 | 2 d -λ 2,2 α 1,2 P 2 + |h 2,2 w 1 | 2 d -λ 2,2 P 1 + σ 2 , ( 10
)
Since the inter-user interference and inter-cluster interference cannot be removed for the far user, it is required to calculate the distribution of the squared-norm of the product between channel gain and the beamforming vector for further analysis.

V. OUTAGE PROBABILITY ANALYSIS

We define the outage probability as the probability that the achieved data rate of an individual user is below a predefined value. According to the Shannon capacity and the SINR expressions in ( 7) and ( 9), the data rate of the near user and the far user in cluster 1 can be expressed as

R 1,1 = log(1 + SINR 1,1 ) = log ( 1 + χ 1,1 α 1,1 d -λ 1,1 P 1 σ 2 ) , (11) R 2,1 = log(1 + SINR 2,1 ) = log ( 1 + |h 2,1 w 1 | 2 d -λ 2,1 α 2,1 P 1 |h 2,1 w 1 | 2 d -λ 2,1 α 1,1 P 1 + |h 2,1 w 2 | 2 d -λ 2,1 P 2 + σ 2 )
.

A. Far user

The outage probability of the far user can be directly defined as the probability that the data rate falls below the target rate.

Pr(R 2,1 ≤ γ 2,1 ) = (12)
Pr

( |h 2,1 w 1 | 2 d -λ 2,1 α 2,1 P 1 |h 2,1 w 1 | 2 d -λ 2,1 α 1,1 P 1 + |h 2,1 w 2 | 2 d -λ 2,1 P 2 + σ 2 ≤ e γ2,1 -1 )
Considering high SNR regimes, the noise power is negligible, and hence, the outage probability of the far user in cluster 1 can be approximated by a lower bound ignoring the noise power. Consequently, after some simple mathematical manipulations, the outage probability can be summarized as

Pr(R 2,1 ≤ γ 2,1 ) ≈ Pr ( α 2,1 P 1 e γ2,1 -1 -α 1,1 P 1 ≤ |h 2,1 w 2 | 2 |h 2,1 w 1 | 2 P 2 ) . ( 13
)
To derive the closed-form outage expression of the far user, it is necessary to know the ratio distribution between |h

(b) h1,1| 2 |h 1 2,1 h 2 1,2 -h 2 2,1 h 1 1,2 2
Fig. 2: Curve fitting method and |h 2,1 w 2 | 2 , where h 2,1 is a row vector with complex normal random variables for Rayleigh channels, but the distribution of the column vector w n is unknown. Hence, we will analyze the expressions, both h 2,1 w 1 and h 2,1 w 2 in the following. First, according to the beamforming matrix, the beamforming vectors w 1 and w 2 can be obtained as

w 1 =   h 2 1,2 |h1,1| h 1 1,1 h 2 1,2 -h 1 1,2 h 2 1,1 -h 1 1,2 |h1,1| h 1 1,1 h 2 1,2 -h 1 1,2 h 2 1,1   , w 2 =   -h 2 1,1 |h1,2| h 1 1,1 h 2 1,2 -h 1 1,2 h 2 1,1 h 1 1,1 |h1,2| h 1 1,1 h 2 1,2 -h 1 1,2 h 2 1,1   . ( 14
)
Then, the ratio between |h 2,1 w 1 | 2 and |h 2,1 w 1 | 2 can be described as

|h 2,1 w 2 | 2 |h 2,1 w 1 | 2 = |h 1,2 | 2 |h 2 2,1 h 1 1,1 -h 1 2,1 h 2 1,1 | 2 |h 1,1 | 2 |h 1 2,1 h 2 1,2 -h 2 2,1 h 1 1,2 | 2 , ( 15
)
where

h 1 1,1 , h 2 1,1 , h 1 1,2 , h 2 1,2 , h 1 2,1
and h 2 2,1 are i.i.d complex, zero mean random normal distributed variables. The distribution of product of two complex normal distributed random variables follows the product normal distribution, however this distribution is not helpful to find tractable expressions for outage probabilities. To calculate [START_REF] Cui | Outage probability constrained MIMO-NOMA designs under imperfect CSI[END_REF], we need to find the ratio distribution between

h 1,2 | 2 |h 2 2,1 h 1 1,1 -h 1 2,1 h 2 1,1 2 and h 1,1 | 2 |h 1 2,1 h 2 1,2 -h 2 2,1 h 1 1,2
2 . Based on the study in [START_REF] Romero-Jerez | On the distribution of the squared norm of non-circular complex gaussian random variables with applications[END_REF],

the squared norm of complex normal random variables can be approximated as exponential distributed random variables when real and imaginary parts are independent and they have the same variance values. We also verify this by using a curve fitting method as shown in Fig. 2. It can be clearly observed that the PDFs of numerator and denominator of |h2,1w2| 2 |h2,1w1| 2 can be accurately approximated by an exponential distribution. Therefore, the ratio of two squared norm expressions can be approximated as the ratio between two independent exponential distributions as X Y , where X ∼ Exp(λ 1 ) and Y ∼ Exp(λ 2 ). The corresponding probability density functions (PDF) of X and Y are

f X (x) = λ 1 e -λ1x , ( 16
)
f Y (y) = λ 2 e -λ2y .
The main challenge of this step is to derive the PDF and CDF of this ratio distribution. First, assuming Z = X Y as another random variable representing the ratio of the random variable X and Y . Using the method of transformation [START_REF] Mathai | Jacobians of matrix transformations and functions of matrix arguments[END_REF], the PDF of Z is derived as

Z = X Y , X = ZY, ( 17
, f X,Y (x, y) is derived as f X (x)f Y (y) = λ 1 λ 2 e -λ1x+λ2y , hence for X = Y Z, f X,Y (yz, y) = λ 1 λ 2 e -(λ1z+λ2)y .
Thus, the marginal PDF of random variable Z is given by

f Z (z) = ∫ ∞ -∞ f (yz, y)|J|dy = ∫ ∞ 0 f (yz, y)|J|dy (18) = λ 1 λ 2 ∫ ∞ 0 e -(λ1z+λ2)y ydy. Let θ = λ 1 z + λ 2 , f Z (z) = λ 1 λ 2 ∫ ∞ 0 ye -θy dy. ( 19
)
By definition, E(z) = ∫ ∞ 0 θze -θz dz = 1 θ which is the mean of an exponential distribution with parameter θ, f Z (z) = λ1λ2 θ 2 . Therefore, the PDF of the ratio distribution is

f Z (z) = λ 1 λ 2 (λ 1 z + λ 2 ) 2 . ( 20
)
The CDF of the ratio distribution will be

F Z (z) = ∫ z 0 f (t)dt = ∫ z 0 λ 1 λ 2 (λ 1 t + λ 2 ) 2 dt = λ 1 z λ 1 z + λ 2 . ( 21
)
According to the PDF and CDF of the ratio distribution, the closed-form outage expression is derived as

Pr(R 2,1 ≤ γ 2,1 ) = Pr(SINR 2,1 ≤ e γ2,1 -1) (22) = 1 -Pr ( |h 2,1 w 2 | 2 |h 2,1 w 1 | 2 ≤ α 2,1 P 1 (e γ2,1 -1)P 2 - α 1,1 P 1 P 2 ) = 1 - λ 1 λ 1 + λ2 ( α 2,1 P 1 (e γ 2,1 -1)P 2 - α 1,1 P 1 P 2 ) .
Subsequently, based on the condition of transmit power allocation coefficient that α 2,1 = 1 -α 1,1 and introducing a new variable, power ratio coefficient ρ = P1 P2 , we have

Pr(R 2,1 ≤ γ 2,1 ) = 1 - λ 1 λ 1 + λ2(e γ 2,1 -1) ρ(1-α1,1e γ 2,1 ) , ( 23
)
where

λ 1 = λ 2 , since |h 1,2 | 2 |h 2 2,1 h 1 1,1 -h 1 2,1 h 2 1,1 | 2 and |h 1,1 | 2 |h 1 2,1 h 2 1,2 -h 2 2,1 h 1 1,2 | 2
have the same distribution. Therefore, the final outage expression of the far user in cluster 1 will be

Pr(R 2,1 ≤ γ 2,1 ) = e γ2,1 -1 e γ2,1 (1 -α 1,1 ρ) -(1 -ρ) . ( 24
)
The outage probability of the far user in cluster 2 has the similar closed-form expression as

e γ 2,2 -1 e γ 2,2 (1- α 2,1 ρ )-(1-1 ρ ) .

B. Near user

As SIC is performed at the near user, the outage probability of the near user also depends on the successful SIC of the far user's signal. Therefore, the outage probability of the near user in cluster n can be defined as Pr out 1,n = 1 -Pr suc 1,n × Pr suc 2,n .

(25)

(n = 1, 2)
First, the success probability of the near user can be obtained based on the cumulative distribution function (CDF) of an exponential random variable, which is represented as

Pr suc 1,n (R 1,n ≥ γ 1,n ) = Pr ( χ 1,n ≥ (e γ1,n -1)σ 2 α 1,n d -λ 1,n P n ) (26) = 1 -γ ( 2, (e γ1,n -1)σ 2 α 1,n d -λ 1,n P n )
where γ 1,n is the minimum reserved rate, and P n is the transmit power in cluster n.

Then, the probability of successful detection of the far user's signal by the near user can be formulated as

Pr suc 2,n (R ′ 2,n ≥ γ ′ 2,n ) (27) =      1 -γ ( 2, (e γ ′ 2,n -1)σ 2 d -λ 1,n Pn ( 1-α1,ne γ ′ 2,n ) ) , f ≥ 0 0, f < 0 where R ′ 2,n = log ( 1 + |h1,n| 2 d -λ 1,n α2,nPn |h1,n| 2 d -λ 1,n α1,nPn+σ 2
) and f = 1e γ ′ 2,n α 1,n . Finally, the outage probability of the near user in cluster n can be obtained.

VI. SIMULATION RESULTS

To validate the accuracy of the analytical outage expressions, Monte Carlo simulations are provided to assess the precision of the closed-form outage expressions for both users. The transmit power allocation coefficient α 1,1 is smaller than α 1,2 due to the properties of power control for downlink NOMA system. The value of the transmit power for cluster 1 is fixed to P 1 = 10.

Fig. 3 shows the relationship between the outage probability of the far user in cluster 1 and the target rate. The power ratio coefficient ρ = P1 P2 acts a key factor in the closed-form outage expression for far user. Thus, we consider two cases, ρ = 1 and ρ = 10. For ρ = 1 the power levels of two clusters are the same which also implies that both clusters are at the same distance from BS, while ρ = 10 represents a case where the power level of cluster 1 is larger than the one for cluster 2. We also assume that d 2,1 = 10 and d 1,1 = 1 unless otherwise stated. It is observed that the outage probabilities of far user increase with the target rate. For comparison, the outage probability with same power level ρ = 1 is greater than in the case of ρ = 10. It is also shown that the analytical results match well with the numerical results. To further study the effects of power ratio ρ on outage probabilities, in Fig. 4, we compare the outage probabilities of far users for both clusters. We observe that the outage probability decreases with ρ for cluster 1, but it increases with ρ for cluster 2. In particular, the outage probabilities cross together when power level is the same (i.e., ρ = 1). It also confirms that the analytical results for outage probabilities have as high accuracy compared with the numerical results.

Fig. 5 illustrates the outage probability of the near user versus target rate with various values of α 1,n . The SNR of cluster n is fixed to 20dB. As can be observed, the outage probability of the near user improves as the target rate increases. Moreover, the outage probability decreases with the growth of α 1,n and the outage performance becomes worse as the distance of the near user from BS increases. Fig. 6 demonstrates the effects of SNR (which is defined as P n /σ 2 ) on the outage probability of the near user in cluster n. For γ 1,n = 1, we can see that the outage probability decreases with the SNR. Moreover, the results indicate that the outage probability of the near user decreases with the growth of α 1,n since the allocated power to the far user becomes smaller. In addition, the outage probability increases with the distance from BS to near user. The outage performance of the near user has a fine agreement between the numerical and analytical results, and verifies the validity of our outage probability analysis.

VII. CONCLUSIONS

In this paper, the closed-form outage expressions of individual users have been derived for downlink MISO-NOMA system in case of two antennas. The outage probability of far users can be approximated as the CDF of the ratio of two exponentially distributed random variables, which can be solved by applying Jacobian transformation method. On the other hand, the outage probability of the near users for both clusters is proven to be in the form of the CDF of an exponential distribution. This work is the first step to derive tractable expressions for the outage probability of MISO NOMA systems under assumption of statistical CSI and hence it only investigates the case of two antennas. In the future, we will extend this work for a general case of N antennas.

Fig. 1 :

 1 Fig. 1: Downlink MISO-NOMA system

  According to the properties of joint PDF of x and y

Fig. 3 :Fig. 4 :

 34 Fig. 3: The outage probability of the far user vs target rate with different α 1,1

Fig. 5 :Fig. 6 :

 56 Fig. 5: The outage probability of the near user vs target rate with different α 1,n (SNR = 20dB)

  2,1 w 1 | 2

	4500												4500
	4000					Histogram PDF of numerator Exponential approximation		4000					Histogram PDF of denominator Exponential approximation
	3500												3500
	3000												3000
	2500												2500
	2000												2000
	1500												1500
	1000												1000
	500												500
	0	0	2	4	6	8	10	12	14	16	18	20	0	0	2	4	6	8	10	12	14	16	18	20
	(a) h1,2| 2 |h 2 2,1 h 1 1,1 -h 1 2,1 h 2 1,1	2

VIII. ACKNOWLEDGEMENT

The authors would like to thank Jakob Hoydis (Bell Lab) and Prof Francois Baccelli (INRIA) for their helpful discussion and comments.

This work has been supported by the Engineering and Physical Sciences Research Council through grant EP/R006385/1. A part of the work was carried out at LINCS (www.lincs.fr)


Nicolas Bilger
email: bilger@lms.polytechnique.fr

François Auslender
email: auslende@lms.polytechnique.fr

Michel Bornert
email: bornert@lms.polytechnique.fr

Renaud Masson
email: renaud.masson@edf.fr

New bounds and estimates for porous media with rigid perfectly plastic matrix

We derive new rigorous bounds and self-consistent estimates for the effective yield surface of porous media with a rigid perfectly plastic matrix and a microstructure similar to Hashin's composite spheres assemblage. These results arise from a homogenisation technique that combines a pattern-based modelling for linear composite materials and a variational formulation for nonlinear media.

solids and structures / porous media / bounds / perfect plasticity

Nouvelles bornes et estimations pour les milieux poreux à matrice rigide parfaitement plastique

Résumé

On présente de nouvelles bornes rigoureuses et estimations autocohérentes pour la surface de plasticité macroscopique de milieux poreux à matrice rigide parfaitement plastique dont la microstructure s'apparente à l'assemblage des sphères composites de Hashin. Elles sont déduites d'une technique d'homogénéisation qui combine une approche par motifs pour les composites linéaires et une formulation variationnelle pour les milieux non linéaires.

solides et structures / milieux poreux / bornes / plasticité parfaite

Version française abrégée

, qui permet une description améliorée des interactions dans l'environnement immédiat d'un pore.

La démarche est la suivante. Les équations du problème local non linéaire sont traitées selon les principes de la méthode sécante modifiée [START_REF] Suquet | Overall properties of nonlinear composites: a modified secant moduli theory and its link with Ponte Castañeda's nonlinear variational procedure[END_REF], qui rejoint l'approche variationnelle de Ponte Castañeda [START_REF] Castañeda | New variational principles in plasticity and their application to composite materials[END_REF]. Contrairement aux mises en oeuvre usuelles de ces principes, on fait ici appel, comme dans [START_REF] Bornert | Self-consistent approaches and strain heterogeneities in two-phase elastoplastic materials[END_REF], à un composite linéaire de comparaison multiphasé, avec plusieurs modules sécants dans la matrice, et non pas un seul module constant. Plus précisément, la microstructure est classiquement [START_REF] Perrin | Contribution à l'étude théorique et numérique de la rupture ductile des métaux[END_REF][START_REF] Gȃrȃjeu | Contribution à l'étude du comportement non linéaire de milieux poreux avec ou sans renforts[END_REF] supposée s'apparenter à l'assemblage des sphères composites de Hashin [START_REF] Hashin | The elastic moduli of heterogeneous materials[END_REF], formées d'un pore sphérique entouré d'une gaine de matrice, de taille variable mais de fraction volumique égale à la porosité globale et réparties aléatoirement dans un volume élémentaire représentatif de façon à couvrir tout l'espace (Fig. 1a). Dans le but de mieux représenter l'hétérogénéité radiale des champs mécaniques induite par la plasticité, les couronnes de matrice sont discrétisées en N couches concentriques (Fig. 1b) et l'on affecte à chaque couche r un module sécant différent µ r . Le milieu linéaire de comparaison ainsi obtenu admet une description par motifs morphologiques représentatifs [START_REF] Stolz | Analyse morphologique et approches variationnelles du comportement d'un milieu élastique hétérogène[END_REF], les motifs étant les sphères composites multicouches ainsi définies. L'homogénéisation de son comportement peut être abordée analytiquement au travers de l'approche variationnelle présentée dans [START_REF] Stolz | Analyse morphologique et approches variationnelles du comportement d'un milieu élastique hétérogène[END_REF][START_REF] Bornert | Morphologically representative pattern-based bounding in elasticity[END_REF], qui conduit, sous l'hypothèse réaliste que les motifs sont répartis de manière isotrope, au problème d'une d'inclusion composite multicouche noyée dans un milieu de référence homogène sous chargement homogène à l'infini (Fig. 1c), dont la solution analytique est connue [START_REF] Hervé | N -layered inclusion-based micromechanical modelling[END_REF]. Suivant le choix du milieux de référence, une borne supérieure [START_REF] Stolz | Analyse morphologique et approches variationnelles du comportement d'un milieu élastique hétérogène[END_REF][START_REF] Bornert | Morphologically representative pattern-based bounding in elasticity[END_REF] ou une estimation autocohérente [START_REF] Bornert | Self-consistent approaches and strain heterogeneities in two-phase elastoplastic materials[END_REF][START_REF] Bornert | Morphologically representative pattern-based modelling in elasticity[END_REF] des propriétés du composite linéaire, et par suite de la surface de plasticité du milieu poreux, sont obtenues. Les équations à traiter sont résumées par le système (3), où σ 0 est la contrainte d'écoulement de von Mises de la matrice, T le taux de triaxialité macroscopique (rapport entre contrainte moyenne m et équivalente eq ) et la fonction F fait référence à l'homogénéisation du milieu linéaire de comparaison. Les dérivées partielles de cette dernière donnent accès au moment d'ordre deux ε 2 eq r de la déformation dans les couches selon (2) et sont calculées numériquement.

La résolution de (3) par une méthode de point fixe pour plusieurs valeurs de T permet de construire la fonction seuil point par point. Les résultats pour une porosité f de 1 % et 5 % sont donnés en Fig. 2. Avec un nombre de couches suffisamment grand (N 10) et pour un chargement purement hydrostatique, on converge vers la solution exacte (Fig. 2b), également donnée par le critère de Gurson. Avec une seule couche, on retrouve pour tout T la borne de Hashin et Shtrikman non linéaire (NLHS+). Sous chargement purement déviatorique, la discrétisation n'a qu'un faible effet (Fig. 2a), mais elle permet d'améliorer les bornes et estimations existantes (Fig. 2c) dont en particulier celles à caractère analytique dues à Gȃrȃjeu et Suquet [START_REF] Gȃrȃjeu | Contribution à l'étude du comportement non linéaire de milieux poreux avec ou sans renforts[END_REF][START_REF] Gȃrȃjeu | Effective properties of porous ideally plastic or viscoplastic materials containing rigid particles[END_REF]. Le gain est plus net pour des triaxialités intermédiaires. L'estimation autocohérente est très proche de la borne supérieure, mais s'en écarterait pour une porosité plus élevée.

Introduction

During the last 25 years, noticeable progress has been made in the modelling of nonlinear porous media. Gurson, in his pioneering work on a hollow rigid plastic sphere [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth. Part 1 -Yeld Criteria and flow rules for porous ductile media[END_REF], exhibited a macroscopic yield function which involves, besides classical parameters such as the equivalent von Mises stress eq and the yield stress σ 0 , the porosity f and the hydrostatic stress m :

eq σ 0 2 + 2f cosh 3 m 2σ 0 -1 -f 2 = 0 ( 1 )
This criterion has the noticeable property to be an upper bound [START_REF] Perrin | Contribution à l'étude théorique et numérique de la rupture ductile des métaux[END_REF][START_REF] Gȃrȃjeu | Contribution à l'étude du comportement non linéaire de milieux poreux avec ou sans renforts[END_REF] for the particular class of microstructures described by Hashin's [START_REF] Hashin | The elastic moduli of heterogeneous materials[END_REF] composite spheres assemblage (CSA), made up of spherical pores surrounded by concentric shells of matrix, in different sizes but with the same volume fraction, and distributed in such a way that they fill the whole space (Fig. 1a). Moreover, for purely hydrostatic loads, it gives the exact solution. Another rigorous upper bound for the macroscopic plastic yield surface has been proposed later on by Ponte Castañeda [START_REF] Castañeda | New variational principles in plasticity and their application to composite materials[END_REF] and Michel and Suquet [START_REF] Michel | The constitutive law of nonlinear viscous and porous materials[END_REF]; it is based on an extension to nonlinear behaviour of the linear Hashin and Shtrikman upper bound and applies to any microstructure in which the spatial distribution of the pores is isotropic. It will be referred to by NLHS+. When the triaxiality factor T = m / eq goes to zero, this new bound improves on Gurson's criterion. However, for infinite triaxiality, Gurson's criterion keeps predominance. More recently, Gȃrȃjeu [START_REF] Gȃrȃjeu | Contribution à l'étude du comportement non linéaire de milieux poreux avec ou sans renforts[END_REF] has derived another analytical upper bound that applies to Hashin's assemblage, using more refined admissible strain fields in the composite spheres than according to Gurson's analysis. His criterion is very similar to (1), except that the first term has to be multiplied by some complex coefficient τ (f ). Since τ (f ) < 1, his result, referred to as GSH+ in the following, improves on Gurson's criterion for non purely hydrostatic loads and gives the exact result for T = ∞. However, it lies slightly above NLHS+ for T = 0. Several authors [START_REF] Leblond | Exact results and approximate models for porous viscoplastic solids[END_REF][START_REF] Gȃrȃjeu | Effective properties of porous ideally plastic or viscoplastic materials containing rigid particles[END_REF] have suggested to replace τ (f ) by 1 + 2 3 f ; the obtained estimate (GS in the following) of the yield function coincides then with NLHS+ for T = 0 and with (1) for T = ∞, but there is no proof that such a relation yields an upper bound. In the present paper, new bounds and estimates for the yield function of porous media with a microstructure similar to Hashin's CSA are derived; they are consistent with the above mentioned best bounds for T = 0 and T = +∞ and improve on them for any other triaxiality.

We consider a representative volume element of a such a porous medium. The homogeneous and isotropic matrix is rigid and perfectly plastic, follows the von Mises criterion, with yield stress σ 0 , and obeys the normality rule. Under these assumptions and when is submitted to macrohomogeneous boundary conditions, an effective behaviour of the porous material can be defined; it is rigid perfectly plastic and obeys the normality rule [START_REF] Castañeda | Nonlinear composites[END_REF] as well. When, in addition, the spatial distribution of the pores is isotropic, the overall yield criterion depends only on the three invariants of the overall stress . The proposed models, alike those mentioned above, are neither sensitive to the third invariant nor to the sign of the pressure, so that the overall yield function depends only on eq and | m |. To construct it, we determine, for any imposed value of the triaxiality T , the overall stress that leads to some non vanishing overall strain E, with arbitrary magnitude. The full yield function in the (| m |, eq ) plane is then obtained when T varies from 0 to +∞.

We make use of the "modified secant procedure" proposed by Suquet [START_REF] Suquet | Overall properties of nonlinear composites: a modified secant moduli theory and its link with Ponte Castañeda's nonlinear variational procedure[END_REF] which is equivalent to Ponte Castañeda's variational procedure [START_REF] Castañeda | New variational principles in plasticity and their application to composite materials[END_REF]. But unlike these authors, we are dealing with a "linear comparison composite" (LCC), in the sense of Ponte Castañeda, with a piecewise constant but non uniform distribution of shear moduli in the matrix. This enriched formulation of the comparison medium leads to a better description of the plastic flow in the nonlinear medium and is the key of the obtained improvements. The multiphase LCC is handled by an extension of the linear Hashin and Shtrikman procedure to microstructures described by means of morphologically representative patterns [START_REF] Stolz | Analyse morphologique et approches variationnelles du comportement d'un milieu élastique hétérogène[END_REF][START_REF] Bornert | Morphologically representative pattern-based bounding in elasticity[END_REF], which leads to an (almost) analytical solution of the linear problem. The detail of the procedure is given in the next section and the results are commented and compared to above mentioned earlier models in the third one.

Nonlinear variational procedure based on multilayered composite spheres

2.1. Modified secant approach. -If we assume for simplicity that the local deviatoric stress evolves in a proportional and monotonic way, the deformation theory of plasticity applies. The matrix being isotropic and incompressible, its local behaviour may be described by a secant constitutive law involving a local secant shear modulus µ(ε eq ) which depends at each point x of on the local equivalent strain ε eq (x). The nonlinear problem can be simplified by replacing this distribution of local moduli by a piecewise uniform one, the secant moduli being equal to µ r in the sub-domain r of the matrix. This defines the so-called LCC, which is an N + 1 phase composite made up of the pores and the N sub-phases r , r ∈ [1, N], in the matrix, with volume fractions c r . According to the modified secant approach, the moduli µ r should be related, through the local secant constitutive law of the matrix, to some reference strain equal to the secondorder moment of the local strain field in the LCC over the domain r : ε r eq = ε 2 eq r where a r is the volume average of a over phase r. With such a choice, the effective properties of the LCC lead to an upper bound for the nonlinear composite [START_REF] Castañeda | New variational principles in plasticity and their application to composite materials[END_REF], in terms of overall potentials. In the present case, one gets an outer bound for the yield surface. When the overall properties of the LCC are isotropic, with µ hom and k hom the effective shear und bulk moduli, when an overall triaxiality T is prescribed and when an overall equivalent strain E eq is imposed, the general relation for the second order moments given in [START_REF] Suquet | Overall properties of nonlinear composites: a modified secant moduli theory and its link with Ponte Castañeda's nonlinear variational procedure[END_REF] specialises into:

ε 2 eq r = E 2 eq c r ∂µ hom ∂µ r + 3 T µ hom k hom 2 ∂k hom ∂µ r (2)
The nonlinear set of equations to determine the overall plastic limit is then given by ( 3), where F provides the linear isotropic effective properties of the multiphase LCC. This system has to be solved in the µ r by, for instance, an iterative procedure such as the classical fixed point method,

         ∀r ∈ [1, N], µ r = σ 0 3 ε 2 eq r
, ε 2 eq r given by ( 2)

k hom , µ hom = F c s , µ s , s ∈ [1, N] eq = 3µ hom E eq and m = T eq (3)
2.2. Multilayered composite spheres assemblage. -The full implementation of this algorithm requires the computation of the function F , which depends on the phase distribution. Alike other approaches mentioned above, we assume that the microstructure of the porous medium belongs to the class of composite spheres assemblages (Fig. 1a), made up of hollow spheres with a pore volume fraction equal to the macroscopic porosity, but different radii so that the whole domain can be filled with such spheres. When this microstructure is submitted to plastic deformation, strong deformation gradients may appear near the pores. Under hydrostatic pressure, these gradients are mostly along the radii of the composite spheres. This observation suggests that, in order to describe more accurately the local plastic flow around the pores, it makes sense to discretise the matrix into concentric shells, as suggested in (Fig. 1b). As in [START_REF] Bornert | Self-consistent approaches and strain heterogeneities in two-phase elastoplastic materials[END_REF], the multiphase LCC the present approach relies on is made of multilayered composite spheres, consisting of a pore surrounded by N concentric layers. For simplicity, we assume that the geometry of the composite spheres of different sizes are homothetic. Thus the rth phase, r ∈ [1, N], of our multiphase LCC, with uniform shear modulus µ r , will be the union of all rth layers of the composite spheres. Such a phase distribution can be described by means of "morphologically representative patterns" [START_REF] Stolz | Analyse morphologique et approches variationnelles du comportement d'un milieu élastique hétérogène[END_REF], a pattern being, in the present case, a multilayered sphere of a given size. The derivations presented in [START_REF] Stolz | Analyse morphologique et approches variationnelles du comportement d'un milieu élastique hétérogène[END_REF][START_REF] Bornert | Morphologically representative pattern-based bounding in elasticity[END_REF], based on the Hashin and Shtrikman variational formulation, allow us to construct bounds and generalised self-consistent estimates [START_REF] Bornert | Self-consistent approaches and strain heterogeneities in two-phase elastoplastic materials[END_REF][START_REF] Bornert | Morphologically representative pattern-based modelling in elasticity[END_REF] for the effective properties of such composites. It is shown that, under the assumption of an isotropic distribution of the centres of the patterns, the effective properties can be deduced from the local fields that arise in such a multilayered inclusion embedded in an infinite linear isotropic reference medium (Fig. 1c), when it is submitted to homogeneous strain or stress conditions at infinity. Estimates of the effective shear (resp. bulk) modulus are obtained from the ratio between the average deviatoric (resp. hydrostatic) stress and the average deviatoric strain (resp. dilatation) in the composite inclusion. When the infinite medium has the properties of the stiffest layer, these estimates are indeed upper bounds for µ hom (resp. k hom ), whereas when it has the estimated effective properties, so-called generalised self-consistent estimates are generated.

The multilayered composite inclusion problem can be solved analytically. The relations given in [START_REF] Hervé | N -layered inclusion-based micromechanical modelling[END_REF] provide exact and easy to implement expressions for the local and average stresses and strains in the layers, which can be extended to incompressible phases or to voids, as required for the present application. The combination of this solution and the system (3) finally leads to an outer bound and a self-consistent estimate for the effective yield surface of the nonlinear porous material. Note first that the derivatives in (2) are computed numerically because the analytical relations available from [START_REF] Hervé | N -layered inclusion-based micromechanical modelling[END_REF] are too complex. Second, the initial assumption of proportional loading does not restrict the generality of the results, since the obtained local and overall stresses do not depend on the magnitude of the arbitrarily small overall strain. This property could also be justified more rigorously within a formalism using dissipation potentials [START_REF] Castañeda | Nonlinear composites[END_REF], not to be addressed here.

Results

The new bound (Nφ+) obtained for 1% porosity with N = 1, 2, 5 and 20 layers of same thickness as well as the new self-consistent estimate obtained with 20 layers (20φsc) are compared on Fig. 2 This result can easily be predicted from theoretical arguments since the local fields in the CSA are known exactly for a purely hydrostatic load: within any composite sphere, the displacement fluctuation is radial so that the stress and strain components in a local spherical coordinate system depend only on the distance r to the centre of the corresponding pore. For an incompressible matrix, the displacement fluctuation in the matrix varies like 1/r 2 , the local equivalent strain like 1/r 3 and the local secant modulus like σ 0 r 3 . Such a distribution can be approached exactly by the LCC when the discretisation is refined, so that the fields in the nonlinear composite can be described exactly within the secant approach. One can check that the numerical solution of (3) does indeed converge to this distribution of secant moduli. For intermediate triaxialities, the new bound derived with 20 layers is always somewhat lower than Gurson's criterion. For a purely deviatoric load, the latter is definitively too stiff; the new bound is very close to NLHS+ and seems to be independent of the discretisation. In fact, it is not: one can plot in this case the evolution of µ s with respect to r and observe a non constant evolution which becomes more and more significant when the porosity is increased. Thus the new bound is always softer than NLHS+ even if the difference can be neglected in practice for small porosity. As expected, the self-consistent estimate is always within the bound for a similar discretisation, but very close to it for small porosity. The plot on Fig. 2c provides, for a porosity of 5%, a comparison between the new bound (20φ+) and estimate (20φsc) calculated for 20 layers, and the best up to now available rigorous bound (GSH+) and estimate (GS) mentioned in the introduction. In order to emphasise the improvements, it gives, as a function of the polar angle θ = arctan(T ) in the (| m |, eq ) plane, the difference between the 'radii'

ρ = 2 m + 2
eq predicted by these models and that given by Gurson's criterion. Except in a small zone at a very large triaxiality where the number of layers is too small, we observe that 20φ+ improves on GSH+ and 20φsc improves on GS. Moreover, it is worth noting that the 20φ+ bound is tighter than the GS estimate. When the triaxiality tends to infinity, we observe that the new bound and estimate tend to the Gurson bound but are less accurate than the other models. This discrepancy can be reduced with an increased discretisation, but at a significant numerical cost. The GSH+ remains the best fully analytical bound valid for any CSA. When the composite spheres are isotropically distributed, Nφ+ is more restrictive. However, in this case, the less accurate but fully analytical and very simple GS criterion seems to be a reasonable model for more practical applications.

Figure 1 .

 1 Figure 1. Composite spheres assemblage (a); discretisation into layers (b); multilayered inclusion (c).

Figure 1 .

 1 Figure 1. Assemblage de sphères composites (a) ; discrétisation en couches (b) ; inclusion multicouche (c).

  (a and b) to the Gurson model and the Hashin and Shtrikman upper bound NLHS+. The fact that 1φ+ and NLHS+ coincide clearly results from the fact that the linear upper bound involved by both approaches are identical. With more layers, the new bound always improves on NLHS+. For a purely hydrostatic load, (Nφ+) converges to the Gurson model when N is increased. With 10 layers the discrepancy is for instance lower than 1%.

Figure 2 .

 2 Figure 2. (b) Yield surfaces for f = 1%. (a) zoom for small T . (c) Deviation of yield surfaces with respect to Gurson's criterion for f = 5% in (ρ, θ) polar coordinates of the ( m , eq ) plane.

Figure 2 .

 2 Figure 2. (b) Surfaces seuils pour f = 1 %. (a) zoom pour T petit. (c) Écarts entre surfaces seuils et critère de Gurson pour f = 5 % en coordonnées polaires (ρ, θ) du plan ( m , eq ).

Acknowledgements. The authors are grateful to P. Suquet and A. Zaoui for their comments and advices.


T Karimzada

H Maigre

Modelisation of dynamic crack propagation criteria

The aim of this study is the analysis of the dynamic crack propagation criteria in materials under impact loading. We want to define which mechanical quantities monitor the propagation, then use these parameters in numerical simulations and finally compare the results with dynamic propagation experiments. First, we propose a ID semi-analytical model. In this model several mechanical quantities can be used to define the propagation criterion, like the energy release rate ... The effect of the applied loading and the reflections of waves is also very important, leading sometimes to a stick-slip propagation. Next, we present a 2D modelisation. In this case, it is not possible to study analytically the effect each criterion and we need to use a numerical simulation. We propose a numerical method based on finite element. The propagation is simulated by releasing nodes with a the release criterion defined from the propagation criterion. The numerical scheme is chosen to respect the equivalence between the total dissipation and the fracture dissipation. In parallel, we have performed dynamic propagation experiments on PMMA. The samples are DCB specimens loading with a wedge in a Kolsky set-up. During each test, the location of the crack has been recorded 4 times with high speed cameras (I microsecond shutter).

INTRODUCTION

Cracks in structures are always sources of failure particularly when the loading is transient. It is know that dynamic loading affects the toughness of materials but experiments giving precise measurement of the dynamic stress intensity factor at the onset of the propagation are difficult to achieve [1]. Moreover it is not clear that the dynamic toughness is directly related to the dynamic stress intensity factor and one can imagine a relation depending on the whole evolution of the stress at the crack tip till the initiation. After initiation, the study of the dynamic propagation is much more complicated for both aspects theoretical and experimental. For example, the mechanical state at the crack-tip is controlled by two intensity factors: a stress one K1" and a kinematical one K1" [2]. The parameters governing the propagation must also depend on the thermal heating appearing at the crack-tip because of the large amount of dissipated energy released to fracture the material in a very short time [3,[START_REF] Rittel | /nt[END_REF].

The aim of the present study is to propose different tools to analyse the dynamic propagation of cracks especially the interaction between the loading and the crack-tip, the dependence of the criteria on the crack speed and the thermo-mechanical coupling. The first one consists in a !-Dimensional model [5,6]. With this model the relation between the applied lm i ding and the mechanical state at the crack-tip is obtained analytically and the propagation is a solution of a scalar differential equation depending only on the criteria. The second one is 2-Dimensional approach with dynamic propagation experiments coupled with numerical simulations.

I-DIMENSIONAL MODEL

Modelisation

We consider an elastic DCB structure subjected only to transversal displacements. We assume also that the stress state is defined by a shear force only. For example this modelisation can be applied to fiber composite where axial deformation is negligible compared to transverse distortion. Due to symmetry we just consider the upper half of the DBC. We note u(x, t) the transverse displacement, T(x, t) the internal shear, a(t) the length of the crack and F(t) the applied loading at the extremity x = 0 of the DCB (figure 1).

F(t) t X a(t) L Figure 1. Upper half DCB model.
The equations governing this problem are the elastic behaviour and the equation of motion :

T(x, t) = fl S u .• (x, t) (1)
T(x, t) .• = pS u.n(x, t)

flS and pS are respectively the elastic shearing modulus and the density of the material considered for this structure.

Initial and boundary conditions

Before solving the problem we need to fixe initial and boundary condition. Initially we suppose the DCB at rest so we have :

T(x, t) = 0 ; u(x, t) = 0 Vx E [0; L], 'lit< 0 (3)
In x = 0 we apply a given force F(t) which is supposed to be inactive in the past and it is related to the shear force by the following condition:

T(O, t) = -F(t) 'lit
Note that by convention when the crack is opening the shear force T is negative.

(4)

Normally we should prescribe a boundary condition in x = L . In our case the DCB is free along the crack and totally fixed after the crack-tip:

u(x, t) = 0 'l:tx E [a(t) ; L], Vt
The second boundary condition is then written like this: u(a(t), t) = 0 'lit (5) A special attention have to be paid on this condition because this a moving condition with the crack propagation depending on the crack length which is also an unknown variable of the problem.

Criteria of propagation

The evolution of the crack length depends on the criterion of propagation chosen at the crack-tip. With this 1-D model there is no stress singularity and therefore no stress intensity factors to quantify the loading of the crack. To define this loading we calculate the energy release rate G(t) which is defined whatever the mechanical model of fracture is. In our case we finally obtained the following expression:

G(t) = JJ� r(t)2 (1-�) (6)
Where T-(t) = T(a-(t), t)) is the shear force just before the crack-tip and c = ..J!JSipS the celerity of the elastic shear waves in the DCB.

The shear force is not continuous at the.moving crack-tip and the jump is obtained using shock wave relation which leads to the definition of the shear force in front of the crack-tip:

r(t) = T(a+(t), t) = T(a-(t), t) -� T(a-(t), t) = r(t) (1 -�) (7)
Putting this result in the definition of G we obtain :

(8)
This expression is very similar to the expression given in the 2-dimensional model with two stress intensity factors, the stress intensity factor and the kinematic intensity factor.

So in our model a criterion of propagation will appear as a function of these four quantities (G, r, T+, a) in which only two of them are really independent.

Equation of propagation

The general solution of equation ( 1) and ( 2) is a set of two unknown scalar functions f and g:

u(x, t) = f(t-x/c) + g(t + x/c) (9)
Putting this expression in the initial and the boundary conditions we find a explicit relation between the applied loading and the criterion of propagation. This relation gives a non linear differential equation of the first order on a(t). In the first stage (t < t1) when we consider the arrival of the first incident wave this equation becomes:

_ F(t-a(t)/c) T (t) = T(a(t), t) = -2 l +ale Vt < t1 ( 10
)
After a certain time (t1 < t) reflected waves at the crack-tip and at the extremity x = 0 will arrive again at the crack-tip and perturb the propagation. The new differential equation is:

T( () ) __ 2 F(t -a(t)/c)-F1(t) -... .. a t , t - 1 + ale v t > t1 (11)
Where F1(t) is a known function depending on the previously calculated evolution of the crack propagation.

Numerical illustrations

We present two examples of propagation obtained with this model. In these examples the load applied is linear increasing force and we compare two criteria of propagation. The first one is defined by a critical release rate value Gc as:

a = 0 if G(t) < Gc and G(t) = Gc if a > 0
The second one is defined by a critical shear force value T c as:

(12) a=O if -T\t) <Tc and -T'"(t)=Tc if a>O (13)

We consider this second criterion because in the 2-dimensional model it will correspond to a critical stress intensity factor.

All dimensions are taken to I, so the initial crack length is 1, the elastic celerity is 1, ... The critical values are also taken to 1. The applied force is F(t) = 0.5 + 0.4 t so that the propagation just initiates at the arrival of the first wave at the crack-tip. The main result is plotted on the graph figure 2 . We can see in both case that the celerity is always lower than 1 which the maximum theoretical celerity. Even if the applied load is increasing continuously reflections at the crack-tip produce a sudden decrease of the crack velocity. The difference between Gc criterion and Tc criterion is not very important, the first one just leads to faster crack propagation.

Dynamic propagation experiments

Experimental setup

Our experimental setup consists of DCB specimens placed between split Hopkinson pressure bars (SHPB). The opening of the crack is obtained with a 90° wedge (figure 3). Bars are made of maraging and specimen of PMMA. The specimens are transparent and we used 4 electronic high speed camera to take 4 micro-second pictures during test at 4 different moments. The initial crack length is 28 mm so the crack can dynamically propagate over 70 mm. We did 4 experiments with roughly the same impact loading. The velocity of the striker was between 9rnls and 9.7rnls and the force applied to the wedge was about 8kN at the maximum for a duration of 200 11s. We put grease between the wedge and the specimen and we assume a non frictional contact.

Experimental results

Theoretically the cracks should propagate along symmetric plane of the specimen. In fact we observe the symmetry in the beginning of the propagation (-25 mm) and then a slight and smooth deviation occurs. This is due mainly to the non perfectly symmetric contact between the specimen and the bars (figure 4). As we did 4 tests under the same loading we obtained 16 pictures taken at 16 different moments. If we assume the reproductibility of the experiments we have 16 different positions of the crack (figure 5). We can observe a regular propagation of the crack with a slight acceleration at 220 11s. The average celerity is 370 rnls. The loading entering the specimen at 120 J..l S, propagation starts 30 J..I. S after and lasts 200 J..I. S.

These durations are to be compared to 40 flS which is the time for the elastic waves to go from the incident bar to the transmitter bar. So it seems that the crack propagation is not affected by the reflections on the boundaries of the specimen. The examination of the crack surfaces shows the commonly observed feature of conical markings. More precisely we can see a rather mirror aspect at the beginning and a rather rough feature after which is in agreement with the observed acceleration. As we did not have photoelasticity or caustic facilities to measure directly the dynamic stress intensity factors we need to simulate numerically the previous experiments to investigate the a representative criterion of propagation. A first numerical simulation has been done using an elasto-dynamic finite element method. In these simulations the input forces are the experimental data got from the SHPB records and the propagation was simulated by releasing the nodal force at the crack-tip according to the experimental observation. The numerical result is the energy release rate G as function of time computed using path invariant integral. We have obtained a lot of oscillations around a roughly constant average value during the propagation. These oscillations come essentially from the numerical technique used to perform the propagation. In the lD model we obtained regular propagation with some jumps of celerity corresponding to particular reflected waves. In the experiments there should be some sudden variation but not at each time step. The next step will consist in simulating the propagation with a given criterion and find the one giving a numerical crack propagation close to the experimental one.

Conclusions

We have presented some preliminary results on dynamic propagation of cracks.

MFigure 2

 2 Figure 2 Evolution of the crack velocity for two different criteria of propagation

Figure 3 .

 3 Figure 3. SHPB setup with DCB specimen loaded with a wedge

Figure 4 .

 4 Figure 4. Four J.IS pictures showing the position of the crack during a dynamic propagation.

Figure 5 .

 5 Figure 5. Crack position as function of time in 4 dynamic propagation tests 3.4 Numerical simulation

  Theses results are based on a semi-analytical !-Dimensional model and on 2-Dimensional experiments. The ID model gives exact results on the condition of propagation as a function of the criterion. Theses ID results are very useful to find the right way to perform 2D simulations which are the only way to analyze carefully propagation experiments.


Mathilde Bigot

Fabien Chauveau

Ralph Sinkus

Olivier Beuf

Simon A Lambert

Fabien Chauveau²

Characterization of fibril aggregates in ex vivo rat brain with multi-frequency MR Elastography -Preliminary results

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Accurate and early diagnosis of dementia represents a major challenge in an ageing society. Pathological hallmarks of these diseases are fibrils, formed by the aggregation of endogenous proteins misshaped to β-sheet conformation. Accurate localization of these fibrils could improve diagnosis given that diverse dementias protein accumulation starts in different locations.

We suggest using Multi-frequency Magnetic Resonance Elastography (MRE) to characterize those aggregates. MRE has already been used to characterize the spatial distribution of sub-voxel obstacles [START_REF] Jugé | Microvasculature alters the dispersion properties of shear waves--a multi-frequency MR elastography study[END_REF][START_REF] Lambert | Bridging Three Orders of Magnitude: Multiple Scattered Waves Sense Fractal Microscopic Structures via Dispersion[END_REF] . Other groups have found stiffness changes in transgenic mice [START_REF] Majumdar | Early-stage analysis of murine models of Familial Alzheimer's disease: Preliminary results. 1 st MRE Workshop[END_REF][START_REF] Murphy | Magnetic resonance elastography of the brain in a mouse model of Alzheimer's disease: initial results[END_REF][START_REF] Munder | MR elastography detection of early viscoelastic response of the murine hippocampus to amyloid β accumulation and neuronal cell loss due to Alzheimer's disease[END_REF] ; however the substrate of these changes is not well characterized: effects of fibrils, neuronal loss, vascularization or inflammation are not discriminated. The purpose of this work is to characterize the local changes induced by injected fibrils on the brain mechanical properties.

Methods

Four Sprague Dawley rats were injected with 10µL of α-synuclein fibrils at 200µM in the striatum area. The opposite hemisphere received a control injection of Phosphate Buffered Saline (PBS). Animals were sacrificed 2 to 4 weeks after injections. The brains were immediately embedded in 1% agarose gel (SIGMAF8630) with artificial cerebrospinal fluid [START_REF]Preparation of Artificial CSF[END_REF] . The brain was then imaged on a Bruker 4.7T scanner with a home-designed volume coil dedicated to small sample MRE. RARE sequence was used as anatomical image. MRE measurements were performed at 800, 900, 1000 and 1200Hz. 1000Hz acquisition was done first and repeated before the last frequency was acquired, to confirm that brain mechanical properties were not degraded during acquisition time (90 minutes). MRE sequence was derived from RARE sequence [START_REF] Lefebvre | Comparison of viscoelastic property characterization of plastisol phantoms with magnetic resonance elastography and high-frequency rheometry[END_REF] . Motion was encoded in the three directions with the following parameters: 8 snapshot of the waves, 7 1-mm-thick slices, FOV of 2.5x2.5cm², in plane matrix of 64×64, leading to voxel size of 0.391x0.391x1mm [START_REF] Majumdar | Early-stage analysis of murine models of Familial Alzheimer's disease: Preliminary results. 1 st MRE Workshop[END_REF] .

Processing involves a 3D direct curl inversion [START_REF] Sinkus | Elasticity reconstruction: Beyond the assumption of local homogeneity[END_REF] . Reconstruction was performed with the application of a Gaussian filter (σ = 0.39x0.39x1mm 3 , kernel = 3x3x3 voxels) on the unwrapped phase images. Real part k r (inversely proportional to wave speed) of the complex wave number was estimated in manually defined ROIs corresponding to the whole brain, the α-synuclein injection localization and the contralateral injection. Those ROIs are reported on anatomical and k r reconstructed images (Fig. 1). The exponent y of the frequency power law demonstrated by Lambert et al. [START_REF] Lambert | Bridging Three Orders of Magnitude: Multiple Scattered Waves Sense Fractal Microscopic Structures via Dispersion[END_REF] was derived from k r maps reconstructed at different frequencies for each ROI. Allometric fitting (a*frequency y ) was realized using Origin software (OriginPro 2016 64-bit, http://www.OriginLab.com). Paired sample Wilcoxon tests were performed to compare the dispersion in the ROIs.

Results

Mechanical parameters of the brain remained stable during whole imaging duration: for the two 1000Hz acquisitions, reconstructed k r varied in average by 9 ± 6% for the whole brain, 13 ± 9% for the fibril injection and 4 ± 3% for the PBS injection. For the whole brain ROI, this variation is inferior to the standard deviations of each k r measured at 1000 Hz.

No significant difference was found between k r estimated from the three ROIs whatever the frequency investigated (Fig. 2). Interestingly, y was superior in average in the α-synuclein ROI (0.80 ± 0.09) compared to both the whole brain (0.65 ± 0.07) and the contralateral PBS injection (0.58 ± 0.04) (Fig. 3). In two cases, y in the injection ROI was close to y in the whole brain. For two rats, y in PBS seemed close to y in the whole brain, but less for two others. Wilcoxon tests did not qualify this augmentation as significant.

Discussion and conclusions

Single-frequency MRE for fibril aggregates detection does not appear relevant as no change in k r was detected whatever the frequency investigated. These results support recent findings in Alzheimer's transgenic animals suggesting that fibrils alone are not responsible for viscoelasticity changes at 1000Hz 3 .

This work shows for the first time that α-synuclein fibrils, injected in an important quantity without over toxicity and inflammation at the sub-acute stage, can hardly be detected using multi-frequency MRE. An increase of the parameter y is observed at the fibril location. This parameter could be sensitive to microarchitecture alteration, as indicated by Schregel et al. [START_REF] Schregel | Demyelination reduces brain parenchymal stiffness quantified in vivo by magnetic resonance elastography[END_REF] . A decrease of this same parameter is observed in a contralateral PBS injection. These changes are not significant, due to the small number of imaged animals, but also to some values, very close to whole brain dispersion. These values could be explained by a wrong estimation of the injection area when placing the ROIs, which will need validation from ongoing histological analysis.

Additional acquisitions are needed to confirm these results. The sensitivity of the method to detect the presence of fibrils must be studied by decreasing the concentration of injected fibrils.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: a) Central slice of anatomical RARE image of an ex vivo brain embedded in agarose and b) reconstructed k r map of the same brain, derived from MRE acquisitions. On both images are represented the ROIs used to calculate k r in the whole brain (green), an α-synuclein injection (blue) and a control PBS injection (yellow). Needle traces from the injection are visible in both blue and yellow ROIs on the anatomical image.

Acknowledgement

This work was supported by the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR) and by PEPS CNRS "Balanced".

We want to acknowledge the PILoT imaging platform, member of France Life Imaging network (grant ANR-11-INBS-0006).for the support provided on image acquisition.


Minh Quan Dao

Davide Lanza

Vincent Frémont

End-to-End Deep Neural Network Design for Short-term Path Planning

Early attempts on imitating human driving behavior with deep learning have been implemented in an reactive navigation scheme which is to directly map the sensory measurement to control signal. Although this approach has successfully delivered the first half of the driving task -predicting steering angles, learning vehicle speed in an end-to-end setting requires significantly large and complex networks as well as the accompanying dataset. Motivated by the rich literature in trajectory planning which timestamps a geometrical path under some dynamic constraints to provide the corresponding velocity profile, we propose an end-to-end architecture for generating a non-parametric path given an image of the environment in front of a vehicle. The level of accuracy of the resulting path is 70%. The first and foremost benefit of our approach is the ability of incorporating deep learning into the navigation pipeline. This is desirable because the neural network can ease the hardness of developing the see-think-act scheme, while the trajectory planning at the end adds a level of safety to the final output by ensuring it obeys static and dynamic constraint.

I. INTRODUCTION

Motion planning methods for autonomous vehicles are classically developed to sequentially perform path planning, obstacles avoidance, and trajectory optimization [START_REF] Katrakazas | Real-time motion planning methods for autonomous on-road driving: State-ofthe-art and future research directions[END_REF], [START_REF] Gu | Focused trajectory planning for autonomous on-road driving[END_REF], [START_REF] Xu | A real-time motion planner with trajectory optimization for autonomous vehicles[END_REF]. Path planning module, which can be implemented by RRT [START_REF] Kuwata | Motion planning for urban driving using rrt[END_REF], A* [START_REF] Montemerlo | Junior: The stanford entry in the urban challenge[END_REF], Lattices and motion primitives [START_REF] Ziegler | Spatiotemporal state lattices for fast trajectory planning in dynamic on-road driving scenarios[END_REF], or function optimization [START_REF] Ziegler | Trajectory planning for bertha a local, continuous method[END_REF], takes into account the geometry characteristic of the environment to produce a collisionfree (with regard to static obstacle) non-parametric path. On the other hand, the trajectory optimization module, which subsumes the obstacles avoidance, aims to optimize both way points and the corresponding velocity profile so that they respect vehicle's kinematics and dynamic obstacles [START_REF] Liu | Speed profile planning in dynamic environments via temporal optimization[END_REF], [START_REF] Ziegler | Trajectory planning for bertha a local, continuous method[END_REF].

Deep learning can help simplify the path planning and trajectory optimization pipeline above with just a deep neural network which learns human driving behavior in a supervised manner. The first successful example traces back to ALVINN [START_REF] Pomerleau | Alvinn: An autonomous land vehicle in a neural network[END_REF], where a shallow network was used to calculate the steering angle directly from images. This work is revised in the deep learning era in [START_REF] Bojarski | End to end learning for selfdriving cars[END_REF], the authors went beyond a mere pattern recognition, learning the entire steering angle prediction pipeline for autonomous cars by building a mapping from images obtained by a forward camera to steering angle with a deep CNN. Taking the similar approach but with different input, [START_REF] Maqueda | Event-based vision meets deep learning on steering prediction for self-driving cars[END_REF] used visual information obtained by an event-based camera to train their steering model. A more global approach to motion planning using deep learning is to 1 Authors are with Centrale Nantes, LS2N -UMR 6004 Nantes, France email: vincent.fremont@ec-nantes.fr learn a spatial traversability maps [START_REF] Wulfmeier | Watch this: Scalable costfunction learning for path planning in urban environments[END_REF], [START_REF] Wulfmeier | Large-scale cost function learning for path planning using deep inverse reinforcement learning[END_REF]. In these works, rather than just predict a single steering angle, the model is designed to learn a cost function with can be later used for path planning. The common shortage of these works is their inability of addressing vehicle speed.

To solve this problem by the end-to-end learning, [START_REF] Xu | End-to-end learning of driving models from large-scale video datasets[END_REF] uses recurrent layers together with CNN network to learn a complete driving policy (steering angle and vehicle speed). Though this work has showed its capability and robustness, it comes with the cost of an extremely large dataset and fairly complex network as well as the training process. Taking a different approach, [START_REF] Hubschneider | Integrating end-to-end learned steering into probabilistic autonomous driving[END_REF] proposed an integrated solution where a CNN is trained on monocular image data (as in [START_REF] Bojarski | End to end learning for selfdriving cars[END_REF]) to output a path, this path is then used as the initial guess for a Particle Swarm Optimization algorithm [START_REF] Kennedy | A discrete binary version of the particle swarm algorithm[END_REF] which in turn transforms the path to a complete trajectory.

Another method to infer vehicle speed is to calculate it proportional to the collision probability [START_REF] Loquercio | Dronet: Learning to fly by driving[END_REF]. The network designed in this work is made of three consecutive ResNet blocks following by two parallel fully-connected layers, which respectively output the steering angle and a collision probability. The steering prediction is learned through a regression problem, while the collision prediction is addressed as a binary classification problem.

In this paper, we propose and evaluate a deep neural network architecture inspired by DroNet for short-term path planning which is to predict a sequence of steering angles directly from an image obtained by forward camera, hence an end-to-end model. The main difference between [START_REF] Loquercio | Dronet: Learning to fly by driving[END_REF] and our work is the statement of steering angle prediction problem. In fact, as mentioned in [START_REF] Penkov | Applying techniques in supervised deep learning to steering angle prediction in autonomous vehicles[END_REF], it can be transformed from a regression problem of continuous values to a classification problem where the steering angle range is tessellated into discrete spans with width of 0.01 radians. Such choice of span width is justified by the jitter of steering angle applied by a human driver in straight road. Moreover, the calculated steering sequence is mapped into a non-parameterized path and, once the path is output, any motion planning algorithm can be implemented to timestamp the path.

The rest of this paper is organized as follows. Sec.II adapts the DroNet architecture for learning a single steering angle through a classification problem. In Sec.III, we further modify the resulted architecture and the used dataset to enable the network to learn a geometrical path in an endto-end setting. The performance of the path planning model derived in the Sec.III is evaluated in Sec.IV. Then Sec.V describes the short-term path generation using the steering angles sequence and car-like vehicles motion constraints.

Concluding remarks are made in Sec. in VI.

II. LEARNING STEERING ANGLES

Taking a different approach compared to the majority of researches in end-to-end learning for autonomous driving which formulates the prediction of steering angles as a regression problem [START_REF] Pomerleau | Alvinn: An autonomous land vehicle in a neural network[END_REF], [START_REF] Bojarski | End to end learning for selfdriving cars[END_REF], [START_REF] Loquercio | Dronet: Learning to fly by driving[END_REF], our model is designed to be a classifier. The reason of our choice stems from the fact that for a regression-based network, there is a continuous range of angles to infer, but only a finite number of samples with which to train against. By tessellating the range of steering angles into discrete bins, the requirement of infinite training samples to fully cover such continuous range is no longer effective. Moreover, Sec.IV shows that with sufficiently small bins, our classifier can outperform the DroNet -a regression-based model.

A. Network Architecture [START_REF] Loquercio | Dronet: Learning to fly by driving[END_REF] has showed that their the architecture responds strongly to "line-like" features in the forward images which has a strong relation with the resulted steering angles. Motivate by this work, we design our model out of their ResNetmade body and put a classifier on top of it. This classifier is made of 2 dense layers. The first has 800 neurons activated by ReLU, while the second has 227 neurons (equal to the number of classes) and is activated by Softmax function. The conceptual architecture is shown in Fig. 1. Each ResNet block in this figure is comprised of 3 convolutional layers: 2 on the main path and 1 on the shortcut (see Fig. 2).

The hyperparameters of each ResNet block are shown in Tab. I.

TABLE I MODIFIED RESNET CNN BODY PARAMETERS

The intuition behind our model is that the ResNet-made body should learn better how to output useful feature maps which probably contains roads shape and drivable area, while the classifier on its top should learn how to output the steering angle given the provided feature map.

B. Data Preparation

The dataset used to train our model is Udacity dataset challenge 2 1 . This dataset contains several hours of driving on suburban road in good weather and lightning condition.

After processing, the dataset is organized as a time-order list. An element of this list contains an image captured by frontfacing camera, the associated steering angle, GPS coordinate of the vehicle at this time step and other information.

To decide the class of any steering angle, a histogram of all steering angles in the dataset ranging from -2.051 to 1.903 radian is built (see Fig. 4). The width of every bin of this histogram is 1 degree. The class of an angle is the index of the bin it belongs to. As can be seen in Fig. 4, the distribution of collected steering angles among these classes is imbalance. To ensure the network does not overlook the less frequent classes during the training process, class i is assigned a weight w(i) based on the median frequency balancing method in [START_REF] Eigen | Predicting depth, surface normals and semantic labels with a common multi-scale convolutional architecture[END_REF].

w(i) =
median frequency of the dataset frequency of this class (1)

Here, the numerator is the ratio between the number of samples in class i and the total number of samples in the dataset. The denominator is the median of the all frequencies. The resulted weight w(i) is later used to scale the contribution of every sample in class i to the total loss function.

C. Model Training

The model's weights are initialized randomly and the network is trained by minimizing the cross entropy loss function using Adam optimizer with default parameters. After training for 100 epochs with 1200 batch size, our model's accuracy on validation set peaks at 66% before dropping due to over fitting. With this level of accuracy, its qualitative performance measured by Root Mean Square Error (RMSE) and Explained Variance score (EVA) are respectively 0.1083 and 0.8338. These values are competitive, compared to recent development in end-to-end learning model: DroNet [START_REF] Loquercio | Dronet: Learning to fly by driving[END_REF] and [START_REF] Maqueda | Event-based vision meets deep learning on steering prediction for self-driving cars[END_REF]. The details comparison is carried out in Sec.IV.

III. LEARNING A GEOMETRICAL PATH

As shown in the previous section, a single steering angle can be learned through a classification problem. Nevertheless, knowing the steering angle is just half of the autonomous driving task. The other control signal to be provided, is the vehicle's speed. There are no means of inferring a vehicle speed given a single steering angle at the same time instance. However, the rich trajectory planning literature suggests that a geometrical path can be timestamped to generate a velocity profile [START_REF] Katrakazas | Real-time motion planning methods for autonomous on-road driving: State-of-the-art and future research directions[END_REF] such that it can satisfy some dynamic constraint. Therefore, in this section, we modify the resulted architecture from the previous section and the used dataset to enable the network to learn a geometrical path.

A. Network Architecture

Inspired by [START_REF] Hubschneider | Integrating end-to-end learned steering into probabilistic autonomous driving[END_REF], a path can be encoded as a sequence of steering angles, each of which is applied to a predefined traveling distance. Based on this insight, a network can learn a path by learning a sequence of steering angles. This leads to the replacement of a single classifier on top of the

B. Dataset Preparation

The new model's output is interpreted as a sequence of 5 steering angles. Each angle is applied to 2 meters of traveling distance. As a result, a training sample is prepared as following:

• X: an image of the environment in front of the vehicle • y: a list of one-hot vectors. This first vector denotes the class of the steering angle associates with the frame represented by X. The i -th vector represents the steering angle of the frame i × 2 meters away from X.

Such definition of y suggests that the data set used to train path planning model needs to explicitly contain the distance information between two adjacent labels. This distance can be retrieved from the GPS coordinate of each frame provided in the original Udacity dataset. Upon completely being generated, the whole dataset is divided into training set and validation set with respectively 19000 and 2350 samples.

C. Model Training

Since the ResNet-made body is the same the model which predicts a single steering angle in Sec.II, the weight of the ResNet-made body of the path planning model in this section is trained from the best weight obtained in Sec.II. Conversely, all classifiers' weight are initialized randomly. The loss function is chosen to be the cross-entropy and is minimized by the Adam optimizer. The choice of hyperparameters is the same as in Sec. II.

After training for 50 epochs, each of 5 five classifiers in our model achieves at least 70% of accuracy.

IV. MODEL PERFORMANCE

A. Quantitative Performance

Using the same approach of [START_REF] Loquercio | Dronet: Learning to fly by driving[END_REF], [START_REF] Maqueda | Event-based vision meets deep learning on steering prediction for self-driving cars[END_REF], the quantitative performance is measured by Root Mean Square Error (RMSE) and the Explained Variance score (EVA). The performance over these metrics of 5 classifiers in our path planning model compared to a constant estimator, which always predicts 0 as steering angle, a random one, the DroNet [START_REF] Loquercio | Dronet: Learning to fly by driving[END_REF], and the model taking input from an event-based camera [START_REF] Maqueda | Event-based vision meets deep learning on steering prediction for self-driving cars[END_REF] is shown in Tab. II. All classifiers in our model outperform the DroNet in both RMSE and EVA, while fall behind eventbased model with a small margin in RMSE. This shows by tessellating the range of steering angle into sufficiently small intervals, a classification-based model can deliver a better result, compared to a regression-based one.

B. Qualitative Performance

The comparison between the histograms of predicted angle classes and their ground truth on validation set are shown in Fig. 4. This figure indicates a relative match between the predicted distribution and the true distribution.

In addition, the normalized confusion matrix of the first classifier is shown in Fig. 5. This matrix features a clear, large magnitude main diagonal. This means the majority of predicted angle classes is actually the true class. Nevertheless, there are a few strong cells in the bottom of Fig.

C. Layer Activation Visualization

In an attempt to understand how our model produces its prediction, the outputs of each ResNet block in the path planning model are displayed in Fig. 6. This figure shows that the first block recognizes lane mark and vehicles, while the second block segments the drivable area. The last block learns a down-sampled mapping. Together, these three ResNet blocks learn useful feature maps which contains roads shape and drivable area, while the classifiers on the top learn to calculate the steering angle given those feature maps.

V. INTERPRETING A STEERING SEQUENCE AS A PATH

Since the motion of car-like vehicles is constrained to be circular around its Instantaneous Center of Rotation (ICR) (see Fig. 7), a sequence steering angles can be interpreted as a geometrical path (i.e. sequence of way points) by applying each angle in the sequence to a predefined traveling distance of s meters.

In Fig. 7, L is the distance between the front and rear axle. δ i (i = 0 . . . N ) is a steering angle (N + 1 is the length of steering angles sequence). At each time instance, the pose of the vehicle is represented by the pose of the local frame attached to the center of its rear axle -O i x i y i z i (i = 0 . . . N ).

x i goes from the rear axle to the front axle, and is perpendicular to these axles. z i is orthogonal to the plane of motion, and pointing outward. y i is defined such that O i x i y i z i is right-handed. The target is to calculate the position of the center of the front axle relative to the local body frame at the presence -O 0 x 0 y 0 z 0 . Assuming that the vehicle's motion is planar, the transformation from frame O i x i y i z i to frame O i+1 x i+1 y i+1 z i+1 is described by:

i T i+1 = Rot z,φi i t i+1 0 1×2 1 (2)
Here, i t i+1 is the coordinate of O i+1 in frame i, and Rot z,φi is the matrix represents the rotation around the zaxis by an angle φ i . The radius of the circular motion around ICR i is calculated as

R i = L tan δ i (3)
Eq.3 implies that when the steering angle is close to zero (i.e. the steering wheel is kept at the neutral position), the radius of motion approach infinity, hence a straight motion. R i , the coordinate of O i+1 in frame i is

i t i+1 = 2R i sin φ i 2   cos φi 2 sin φi 2   ( 4
)
φ i is the angle between x i and x i+1 . As shown in Fig. 7, this angle can be calculated by:

φ i = s R i = s tan δ i L (5)
With Eq. 5 and Eq. 4. The transformation from frame i to frame i + 1 in Eq. 2 is now fully defined. The local position of the center of the front axle in homogeneous form is:

i L i = L, 0, 1 T (6)
This position is transformed into the local frame at the present time by the following equation

0 L i = 0 T i i L i = i j=1 j-1 T j i L i (7)
To test the quality of the predicted path and its ground truth, the trained path planning model is used to infer path from forward images taken from Udacity dataset. The inference process is implemented on a laptop equipped with an NVIDIA GeForce MX130, an Intel Core i7-8650U (1.90GHz), and 16GB of RAM. The inference time for 1 sample (i.e. 1 image) is 1 millisecond. Examples of path generated by both true and predicted sequence of steering are shown in Fig. 8, while the extended video of path planning model's output compared to ground truth can be found in this link: https://youtu.be/X2fi2xVr2jE. Fig. 8 as well as the video shows a good match between the predicted path and its ground truth, which in turn proves the quality of the prediction of our model.

VI. CONCLUSIONS

In this paper, we explored an end-to-end learning approach to path planning for autonomous vehicles. In details, a neural network made of three ResNet blocks and an array of The advantage of this approach is that it pays the way for the integration of deep neural network into the motion planning framework. Specifically, the geometrical path learned by the network is then timestamped using trajectory optimization technique to finally produce a parameterized path (a path with a velocity profile). Beside providing two crucial control signals of the driving task (steering angle and velocity), this integrated approach enhances the reliability of the predicted trajectory while ensures it is natural and consistent with vehicles' kinematic.

For the future work, the accuracy of the array of classifiers needs to be improved. Furthermore, the network generalization should be evaluated on other autonomous vehicles datasets with different camera parameters. Since a sequence of steering angles implies a time order, it might be helpful if the network can learn a temporal relation among the steering angles. This can be done by exploring the application of recurrent layers such as LSTM cells to enable the model learning such time relation.

The code used in this paper is hosted on GitHub in Minh-Quan Dao's ECN-E2E repository: https://github. com/quan-dao/ECN-E2E.

Fig. 1 .Fig. 2 .

 12 Fig.1. Modified DroNet architecture[START_REF] Loquercio | Dronet: Learning to fly by driving[END_REF]

5

 5

Fig. 3 .Fig. 4 .

 34 Fig. 3. Path planning architecture (CNN part)

Fig. 5 .

 5 Fig. 5. Normalized confusion matrix of the first classifier.

Fig. 6 .

 6 Fig. 6. Output of each ResNet block with respect to different input images. From left to right, the images are respectively the input image, output of the first, second, and third block.

Fig. 7 .

 7 Fig. 7. Car-like vehicle's motion diagram with two different value of steering angles. Each of these angles is applied to an arc distance of s meters.

Fig. 8 .

 8 Fig. 8. Display of the predicted path. Left: Images captured by front-facing camera. Right: Path in cartesian coordinates

https://github.com/udacity/self-driving-car/ tree/master/datasets/CH2

ACKNOWLEDGMENT

This research has been conducted as part of the HI-ANIC (Human Inspired Autonomous Navigation In Crowds) project, funded by the French Ministry of Education and Research and the French National Research Agency (ANR-17-CE22-0010).