Vathana Ly Vath
email: ly-vath@math.jussieu.fr

Mohamed Mnif
email: mohamed.mnif@enit.rnu.tn

Numerical approximation for an impulse control problem arising in portfolio selection under liquidity risk *

Keywords: Impulse control problem, Optimal transaction strategy, Monte Carlo method, Malliavin calculus. JEL Classification : G11 MSC Classification (2000) : 93E20, 65C05, 91B28, 60H30

, in which we suggest a model of liquidity risk and price impact and formulate the problem as an impulse control problem under state constraint. We show that our impulse control problem could be reduced to an iterative sequence of optimal stopping problems. Given the dimension of our problem and the complexity of its solvency region, we use Monte Carlo methods instead of finite difference methods to calculate the value function, the transaction and no-transaction regions. We provide a numerical approximation algorithm as well as numerical results for the optimal transaction strategy.

Introduction

In this article, we investigate numerical aspects of a portfolio selection problem studied in [START_REF] Vath | A model of optimal portfolio selection under liquidity risk and price impact[END_REF], in which we suggest a model of liquidity risk and price impact. Transactions are allowed only in discrete times and incur some fixed costs. Under the impact of liquidity risk, prices are pushed up when buying stock shares and moved down when selling shares. The investor maximizes his expected utility of terminal liquidation wealth, under a solvency constraint. This problem is formulated as an impulse control problem under state constraint. In [START_REF] Vath | A model of optimal portfolio selection under liquidity risk and price impact[END_REF], we characterize the value function as the unique constrained viscosity solution to the associated Hamilton-Jacobi-Bellman Quasi-Variational Inequality (HJBQVI). Our associated HJBQVI has, in addition to time variable, three variables : x, y, and p, respectively the cash holding, the stock holding, and the stock share price.

Hamilton-Jacobi-Bellman equations are usually solved by using numerical methods based on finite difference methods. The Howard algorithm, which consists in computing two sequences: the optimal strategy and the value function, is known to be efficient for the resolution of these types of equation. From Barles and Souganidis [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF], we know that a monotone, stable and consistant scheme insures the convergence of the algorithm to the unique viscosity solution of the HJBQVI. This algorithm has a complexity in O(N n ) where N is the number of points of the grid in one axis and n the dimension of the equation. Chancelier, Øksendal, and Sulem [START_REF] Chancelier | Combined stochastic control and optimal stopping, and application to numerical approximation of combined stochastic and impulse control[END_REF] used the Howard algorithm to solve numerically a bidimensional HJBQVI related to a problem of optimal consumption and portfolio with both fixed and proportional transaction costs. They solved the problem in a bounded domain and they assumed zero Neumann boundary conditions on the localized boundary. The disadvantage of the finite difference method is its suitability to only solve HJB equations when the solvency region has a simple shape such IR n + or when its boundaries are straight. In the latter paper, the solvency region presents some corners. The authors omitted the points of the domain where either the number of shares or the amount of money in the portfolio is non-positive.

Korn [START_REF] Korn | Portfolio optimization with strictly positive transaction costs and impulse control[END_REF] studied the problem of portfolio optimization with strictly positive transaction costs and impulse control. He presented a sequence of optimal stopping problems where the reward function is expressed in terms of the impulse operator. He proved the convergence of the sequence of optimal stopping problems towards the value function of the initial problem. Chancelier, Øksendal and Sulem [START_REF] Chancelier | Combined stochastic control and optimal stopping, and application to numerical approximation of combined stochastic and impulse control[END_REF] suggested an iterative method to solve the impulse control problem. They considered an auxiliary value function where the transactions number is bounded by a positive number.

In this article, we prove that both iterative methods coincide. We study numerically our problem by reducing the impulse control problem to an iterative sequence of optimal stopping problems. Then, we introduce a numerical approximation algorithm for every optimal stopping problem based on ideas of Monte Carlo numerical procedure which requires the computation of many conditional expectations. Several methods can be used for the valuation of these regression functions. We choose the Malliavin Calculus based Method suggested in [START_REF] Fournié | Applications of Malliavin calculus to Monte Carlo methods[END_REF] and then developed in [START_REF] Bouchard | On the Malliavin approach to Monte Carlo approximation of conditional expectations[END_REF]. Our numerical approach named value-iteration algorithm could be adapted to every shape of the solvency region and we don't need to assume some artificial boundary conditions.

The paper is organized as follows. We first show that the value function could be obtained as the limit of an iterative procedure when each step is an optimal stopping problem and the reward function is related to the impulse operator. We then provide a numerical method based on Malliavin calculus and give numerical results for the optimal transaction strategy.

Problem formulation

This section presents the details of the model. Let (Ω, F, IP ) be a probability space equipped with a filtration (F t ) 0≤t≤T supporting an one-dimensional Brownian motion W on a finite horizon [0, T ], T < ∞. We consider a continuous time financial market model consisting of a money market account yielding a constant interest rate r ≥ 0 and a risky asset (or stock) of price process P = (P t ). We denote by X t the amount of money (or cash holdings) and by Y t the number of shares in the stock held by the investor at time t.

Liquidity constraints. We assume that the investor can only trade discretely on [0, T ). This is modelled through an impulse control strategy α = (τ n , ζ n ) n≥1 : τ 1 ≤ . . . τ n ≤ . . . < T are stopping times representing the intervention times of the investor and ζ n , n ≥ 1, are F τn -measurable random variables valued in IR and giving the number of stock purchased if ζ n ≥ 0 or sold if ζ n < 0 at these times. The sequence (τ n , ζ n ) may be a priori finite or infinite. The dynamics of Y is then given by :

Y s = Y τn , τ n ≤ s < τ n+1 (2.1) Y τ n+1 = Y τn + ζ n+1 . (2.2)
Notice that we do not allow trade at the terminal date T , which is the liquidation date.

Price impact. The large investor affects the price of the risky stock P by his purchases and sales : the stock price goes up when the trader buys and goes down when he sells and the impact is increasing with the size of the order. We then introduce a price impact positive function Q(ζ, p) which indicates the post-trade price when the large investor trades a position of ζ shares of stock at a pre-trade price p. In absence of price impact, we have Q(ζ, p) = p. Here, we have Q(0, p) = p meaning that no trading incurs no impact and Q is nondecreasing in ζ with Q(ζ, p) ≥ (resp. ≤) p for ζ ≥ (resp. ≤) 0. Actually, in the rest of the paper, we consider a price impact function in the form

Q(ζ, p) = pe λζ ,
where λ > 0.

(2.

3)

The proportionality factor e λζ represents the price increase (resp. discount) due to the ζ shares bought (resp. sold). The positive constant λ measures the fact that larger trades generate larger quantity impact, everything else constant. This form of price impact function is consistent with both the asymmetric information and inventory motives in the market microstructure literature (see [START_REF] Kyle | Continuous auctions and insider trading[END_REF]).

We then model the dynamics of the price impact as follows. In the absence of trading, the price process is governed by

dP s = P s (bds + σdW s ), τ n ≤ s < τ n+1 , (2.4)
where b, σ are constants with σ > 0. When a discrete trading ∆Y s := Y s -Y s -= ζ n+1 occurs at time s = τ n+1 , the price jumps to P s = Q(∆Y s , P s -), i.e.

P τ n+1 = Q(ζ n+1 , P τ - n+1
).

(2.5)

Notice that with this modelling of price impact, the price process P is always strictly positive, i.e. valued in IR * + = (0, ∞). Cash holdings. We denote by θ(ζ, p) the cost function, which indicates the amount for a (large) investor to buy or sell ζ shares of stock when the pre-trade price is p :

θ(ζ, p) = ζQ(ζ, p).
In absence of transaction, the process X grows deterministically at exponential rate r : ), there is a fixed cost k > 0 to be paid. This results in a variation of cash holdings by ∆X s

dX s = rX s ds, τ n ≤ s < τ n+1 . ( 2
:= X s -X s -= -θ(∆Y s , P s -) -k, i.e. X τ n+1 = X τ - n+1 -θ(ζ n+1 , P τ - n+1 ) -k.
(

The assumption that any trading incurs a fixed cost of money to be paid will rule out continuous trading, i.e. optimally, the sequence (τ n , ζ n ) is not degenerate in the sense that for all n, τ n < τ n+1 and ζ n = 0 a.s. A similar modelling of fixed transaction costs is considered in [START_REF] Pliska | Optimal portfolio management with fixed transaction costs[END_REF] and [START_REF] Korn | Portfolio optimization with strictly positive transaction costs and impulse control[END_REF].

Liquidation value and solvency constraint. The solvency constraint is a key issue in portfolio/consumption choice problem. The point is to define in an economically meaningful way what is the portfolio value of a position in cash and stocks. In our context, we introduce the liquidation function (y, p) representing the value that an investor would obtained by liquidating immediately his stock position y by a single block trade, when the pre-trade price is p. It is given by : (y, p) = -θ(-y, p).

If the agent has the amount x in the bank account, the number of shares y of stocks at the pre-trade price p, i.e. a state value z = (x, y, p), his net wealth or liquidation value is given by :

L(z) = max[L 0 (z), L 1 (z)]1 y≥0 + L 0 (z)1 y<0 , (2.8)
where

L 0 (z) = x + (y, p) -k, L 1 (z) = x.
The interpretation is the following. L 0 (z) corresponds to the net wealth of the agent when he liquidates his position in stock. Moreover, if he has a long position in stock, i.e. y ≥ 0, he can also choose to bin his stock shares, by keeping only his cash amount, which leads to a net wealth L 1 (z). This last possibility may be advantageous, i.e. L 1 (z) ≥ L 0 (z), due to the fixed cost k. Hence, globally, his net wealth is given by (2.8). In the absence of liquidity risk, i.e λ = 0, and fixed transaction cost, i.e. k = 0, we recover the usual definition of wealth L(z) = x + py. Our definition (2.8) of liquidation value is also consistent with the one in transaction cost models where portfolio value is measured after stock position is liquidated and rebalanced in cash, see e.g. [START_REF] Cvitanic | Hedging and portfolio optimization under transaction costs : a martingale approach[END_REF] and [START_REF] Oksendal | Optimal consumption and portfolio with both fixed and proportional transaction costs[END_REF]. Another alternative would be to measure the portfolio value separately in cash and stock as in [START_REF] Deelstra | Dual formulation of the utility maximization problem under transaction costs[END_REF] for transaction cost models. This study would lead to multidimensional utility functions and is left for future research.

We then naturally introduce the liquidation solvency region (see Figure 1) :

S = z = (x, y, p) ∈ IR × IR × IR * + : L(z)
> 0 , and we denote its boundary and its closure by

∂S = z = (x, y, p) ∈ IR × IR × IR * + : L(z) = 0 and S = S ∪ ∂S.
The boundary of the solvency region may then be explicited as follows (see Figures 2 and3) :

∂S = ∂ -S ∪ ∂ y S ∪ ∂ x 0 S ∪ ∂ x 1 S ∪ ∂ x 2 S ∪ ∂ + S, where
∂ -S = z = (x, y, p) ∈ IR × IR × IR * + : x + (y, p) = k, y ≤ 0 ∂ y S = z = (x, y, p) ∈ IR × IR × IR * + : 0 ≤ x < k, y = 0 ∂ x 0 S = z = (x, y, p) ∈ IR × IR × IR * + : x = 0, y > 0, p < kλe ∂ x 1 S = z = (x, y, p) ∈ IR × IR × IR * + : x = 0, 0 < y < y 1 (p)), p ≥ kλe ∂ x 2 S = z = (x, y, p) ∈ IR × IR × IR * + : x = 0, y > y 2 (p), p ≥ kλe ∂ + S = z = (x, y, p) ∈ IR × IR × IR * + : x + (y, p) = k, y 1 (p) ≤ y ≤ y 2 (p), p ≥ kλe .
In the sequel, we also introduce the corner lines in ∂S :

D 0 = {(0, 0)} × IR * + ⊂ ∂ y S, D k = {(k, 0)} × IR * + ⊂ ∂ -S C 1 = {(0, y 1 (p), p) : p ∈ IR * + } ⊂ ∂ + S, C 2 = {(0, y 2 (p), p) : p ∈ IR * + } ⊂ ∂ + S.
Admissible controls. Given t ∈ [0, T ], z = (x, y, p) ∈ S and an initial state Z t -= z, we say that the impulse control strategy α

= (τ n , ζ n ) n≥1 is admissible if the process Z s = (X s , Y s , P s ) given by (2.1)-(2.2)-(2.4)-(2.5)-(2.6)-(2.7
) (with the convention τ 0 = t) lies in S for all s ∈ [t, T ]. We denote by A(t, z) the set of all such policies. We shall see later that this set of admissible controls is nonempty for all (t, z) ∈ [0, T ] × S. In the sequel, for t ∈ [0, T ], z = (x, y, p) ∈ S, we also denote Z 0,t,z s = (X 0,t,x s , y, P 0,t,p s ), t ≤ s ≤ T , the state process when no transaction (i.e. no impulse control) is applied between t and T , i.e. the solution to : starting from z at time t.

dZ 0 s =    rX 0 s 0 bP 0 s    ds +    0 0 σP 0 s    dW s , (2.9)
Investment problem. We consider an utility function U from IR + into IR, strictly increasing, concave and w.l.o.g. U (0) = 0, and s.t. there exist K ≥ 0, γ ∈ [0, 1) :

U (w) ≤ Kw γ , ∀w ≥ 0, (2.10)
We denote U L the function defined on S by :

U L (z) = U (L(z)).
We study the problem of maximizing the expected utility from terminal liquidation wealth and we then consider the value function :

v(t, z) = sup α∈A(t,z) IE e -r(T -t) U L (Z T ) , (t, z) ∈ [0, T ] × S. (2.11)
3 Viscosity solution of the associated Quasi-variational Hamilton-Jacobi-Bellman inequality

In Ly Vath, Mnif and Pham [START_REF] Vath | A model of optimal portfolio selection under liquidity risk and price impact[END_REF], we derived the HJB quasi-variational inequality satisfied by the value function (2.11)

min - ∂v ∂t -Lv , v -Hv = 0, on [0, T ) × S, (3.1)
where L as the infinitesimal generator associated to the system (2.9) corresponding to a no-trading period :

Lϕ = rx ∂ϕ ∂x + bp ∂ϕ ∂p + 1 2 σ 2 p 2 ∂ 2 ϕ ∂p 2 -rϕ,
H is the impulse operator defined by

Hϕ(t, z) = sup ζ∈C(z) ϕ(t, Γ(z, ζ)), (t, z) ∈ [0, T ] × S, Γ is the impulse transaction function defined from S × IR into IR × IR × IR * + : Γ(z, ζ) = (x -θ(ζ, p) -k, y + ζ, Q(ζ, p)), z = (x, y, p) ∈ S, ζ ∈ IR,
and C(z) the set of admissible transactions :

C(z) = ζ ∈ IR : Γ(z, ζ) ∈ S = {ζ ∈ IR : L(Γ(z, ζ)) ≥ 0} .
We related the value function (2.11) and the associated HJB quasi-variational inequality (3.1) by means of constrained viscosity solutions. The definition of viscosity solutions is given as follows:

Definition 3.1 (i) Let O ⊂ S. A locally bounded function u on [0, T ) × S is a viscosity subsolution (resp. supersolution) of (3.1) in [0, T ) × O if for all ( t, z) ∈ [0, T ) × O and ϕ ∈ C 1,2 ([0, T ) × S) s.t. (u * -ϕ)( t, z) = 0 (resp. (u * -ϕ)( t, z) = 0) and ( t, z) is a maximum of u * -ϕ (resp. minimum of u * -ϕ) on [0, T ) × O, we have min - ∂ϕ ∂t ( t, z) -Lϕ( t, z), u * ( t, z) -Hu * ( t, z) ≤ 0 (3.2) ( resp. ≥ 0). ( 3
lim (t , z ) → (t, z) z ∈ S v(t , z ) = 0, ∀(t, z) ∈ [0, T ) × D 0 (3.4) lim (t, z ) → (T, z) t < T, z ∈ S v(t, z ) = max[U L (z), HU L (z)], ∀z ∈ S, (3.5)
and the growth condition :

|v(t, z)| ≤ K 1 + x + p λ γ , ∀(t, z) ∈ [0, T ) × S (3.6)
for some positive constant K < ∞.

Convergence of the iterative scheme

We first introduce the following subsets of A(t, z), the set of the admissible impulse control strategies :

A n (t, z) := {α = (τ k , ξ k ) k=0,...,n ∈ A(t, z)}
and the corresponding value function v n , which describes the value function when the investor is required to trade at most n times:

v n (t, z) := sup α∈An(t,z) IE[e -r(T -t) U L (Z T )] (t, z) ∈ [0, T ] × S. (4.1)
For t ∈ [0, T ] and z = (x, y, p) ∈ S, if x, y are both nonnegative, we clearly have L(Z 0,t,z s ) ≥ 0, and so A 0 (t, z) is nonempty. Otherwise, if x < 0, y ≥ 0 or x ≥ 0, y < 0, due to the diffusion term P 0,t,z , it is clear that the probability for L(Z 0,t,z s ) to be negative before time T , is strictly positive, so that A 0 (t, z) is empty. Hence, the value function for n = 0 is initialized to:

v 0 (t, z) = IE e -r(T -t) U L (Z 0,t,z T ) if x ≥ 0, y ≥ 0 -∞ otherwise
We now show the convergence of the sequence of the value functions v n towards our initial value function v.

Lemma 4.1 For all (t, z) ∈ S lim n→∞ v n (t, z) = v(t, z).
Proof. From the definition of A n (t, z), we have:

A n (t, z) ⊂ A n+1 (t, z) ⊂ A(t, z). As such, v n (t, z) ≤ v n+1 (t, z) ≤ v(t, z),
which gives the existence of the limit and the first inequality:

lim n→∞ v n (t, z) ≤ v(t, z). (4.2)
Given ε > 0, from the definition of v, there exists an impulse control α = (τ 1 , τ 2 , ...; ξ 1 , ξ 2 , ...) ∈ A(t, z) such that

IE[e -r(T -t) U L (Z α T )] ≥ v(t, z) -ε, (4.3)
with Z α diffusing under the impulse control α.

We now set the control

α n := (τ 1 , τ 2 , ..., τ n-1 , τ ; ξ 1 , ξ 2 , ..., ξ n-1 , y τ n-1 ),
where τ n-1 < τ < min{τ n , T }. We see that α n ∈ A n (t, z) and consider the corresponding process Z (αn) . Using Fatou lemma, we obtain:

lim inf n→∞ IE[e -r(T -t) U L (Z (αn) T )] ≥ IE[lim inf n→∞ e -r(T -t) U L (Z (αn) T )] = IE[e -r(T -t) U L (Z α T )](4.4)
Using (4.3) and (4.4), we obtain

lim inf n→∞ v n (t, z) ≥ lim inf n→∞ IE[e -r(T -t) U L (Z (αn) T )] ≥ v(t, z) -ε.
As we obtain the latter inequality with an arbitrary ε > 0, and combining with the relation (4.2), we obtain the desired result:

lim n→∞ v n (t, z) = v(t, z). 2
Theorem 4.1 We define ϕ n (t, z) iteratively as a sequence of optimal stopping problems:

ϕ n+1 (t, z) = sup τ ∈S t,T IE e -r(τ -t) Hϕ n (τ, Z 0,t,z τ ) , ϕ 0 (t, z) = v 0 (t,
ϕ n (t, z) = v(t, z), (t, z) ∈ [0, T ] × S.
so that the iteration scheme for ϕ n provides an approximation for v.

Remark 4.2

The value function ϕ n satisfies the system of variational inequalities, which can be solved by induction starting from ϕ 0 :

min - ∂ϕ n+1 ∂t -Lϕ n+1 , ϕ n+1 -Hϕ n = 0, (t, z) ∈ [0, T ) × S,
together with the terminal condition:

ϕ n+1 (T, z) = Hϕ n (T, z).
Proof of Theorem 4.1. We show by induction that v n (t, z) = ϕ n (t, z), for all n. First, we have v 0 = ϕ 0 . Considering an impulse control strategy α 1 = (τ, ξ) ∈ A 1 (t, z), we clearly have

ϕ 1 (t, z) ≥ IE[e -r(τ -t) Hϕ 0 (τ, Z 0,t,z τ )],
≥ IE[e -r(τ -t) Hv 0 (τ, Z 0,t,z τ )].

From the definition of the operator H, we obtain

ϕ 1 (t, z) ≥ IE[e -r(τ -t) v 0 (τ, Γ(Z 0,t,z τ , ξ))], ∀α 1 = (τ, ξ) ∈ A 1 (t, z). (4.5)
Let Z (α 1 ) be the diffusion of Z, starting at time t, with Z (α 1 ) t = z, and evolving under the impulse control α 1 . Relation (4.5) becomes:

ϕ 1 (t, z) ≥ E[e -r(τ -t) v 0 (τ, Z (α 1 ) τ )], ∀ α 1 = (τ, ξ) ∈ A 1 (t, z). (4.6)
Given the arbitrariness of α 1 and by using the dynamic programming principle applied to v 1 (t, z), we obtain

ϕ 1 (t, z) ≥ v 1 (t, z).
From the definition of ϕ 1 , for a given ε > 0, there exists τ * such that

ϕ 1 (t, z) -ε ≤ IE[e -r(τ * -t) Hϕ 0 (τ * , Z 0,t,z τ * )]. (4.7)
From the compactness of the set of admissible transactions, there exists ξ * such that

ϕ 1 (t, z) -ε ≤ IE[e -r(τ * -t) v 0 (τ * , Γ(Z 0,t,z τ * , ξ * ))], ≤ E[e -r(τ * -t) v 0 (τ * , Z ( * ) τ * )],
where Z ( * ) is the processus starting at time t, with Z ( * ) t = z, and evolving under the impulse control α * := (τ * , ξ * ). Using the dynamic programming principle applied on v 1 (t, z), we obtain

ϕ 1 (t, z) -ε ≤ v 1 (t, z).
The latter inequality is satisfied for any value of ε > 0, as such, we have

ϕ 1 (t, z) ≤ v 1 (t, z), which leads to ϕ 1 (t, z) = v 1 (t, z), for all (t, z) ∈ [0, T ) × S.
By induction, assuming that for a given n, we have ϕ n (t, z) = v n (t, z), we will prove that ϕ n+1 (t, z) = v n+1 (t, z). By definition, we have for any

α n+1 = (τ 1 , ..., τ n+1 , ξ 1 , ..., ξ n+1 ) ∈ A n+1 (t, z), ϕ n+1 (t, z) ≥ IE[e -r(τ 1 -t) Hϕ n (τ 1 , Z 0,t,z τ 1 )], ≥ IE[e -r(τ 1 -t) v n (τ 1 , Γ(Z 0,t,z τ 1 , ξ 1 ))], ≥ IE[e -r(τ 1 -t) v n (τ 1 , Z (n+1) τ 1 )], (4.8)
where Z (n+1) is the diffusion starting at time t, with Z (n+1) t = z and evolves under the control α n+1 . Given the arbitrariness of the control α n+1 and by using the dynamic programming principle applied to v n+1 , relation (4.8) becomes:

ϕ n+1 (t, z) ≥ v n+1 (t, z).
To prove the opposite inequality, we use the definition of ϕ n+1 . For any ε > 0, there exists τ * such that

ϕ n+1 (t, z) -ε ≤ IE[e -r(τ * -t) Hϕ n (τ * , Z 0,t,z τ * )], (4.9)
≤ IE[e -r(τ * -t) Hv n (τ * , Z 0,t,z τ * )]. (4.10)

From the compactness of the set of admissible transactions, there also exists ξ * such that

Hv n (τ * , Z 0,t,z τ * ) = v n (τ * , Z (α * ) τ * ),
where Z (α * ) , the processus starting at time t, with Z t = z, evolves under the impulse control α * := (τ * , ξ * ). Using the dynamic programming principle applied on v n+1 , the relation (4.10) becomes

ϕ n+1 (t, z) -ε ≤ IE[e -r(τ * -t) v n (τ * , Z (α * ) τ * )], ≤ v n+1 (t, z).
The inequality is obtained for any given ε, this leads to the required inequality

ϕ n+1 (t, z) = v n+1 (t, z).

5 Numerical study

The objective of this section is the computation of a sequence of optimal stopping problem:

v n+1 (t, z) = sup τ ∈S t,T
E e -r(τ -t) Hv n (τ, X 0,t,x τ , y, P 0,t,p τ ) , z ∈ S and the associated trade region and the no-trade region. We choose the Monte Carlo numerical procedure for the implementation.

The Monte Carlo method

Let T m = {t l = lT /m} 0≤l≤m be the partition of the time interval T = [0, T ]. We denote by h t the time step T m , and by S m,t,T the subset of S t,T defined by S m,t,T = {τ ∈ S t,T ; τ ∈ T m }.

Let h z := (h x , h y , h p ) = (1/M 1 , 1/M 2 , 1/M 3 ), where (M 1 , M 2 , M 3 ) ∈ IN * 3
denotes the finite difference step in the state coordinate z = (x, y, p). Since the liquidation solvency region is unbounded, we localize S to

D = {z ∈ S s.t. -L 1 ≤ x ≤ L 1 , -L 2 ≤ y ≤ L 2 , 0 ≤ z ≤ L 3 },
where L 1 , L 2 and L 3 are positive constants. We define the grid

Ω hz = {z = (ih x , jh y , kh p ) ∈ D, -M 1 L 1 ≤ i ≤ M 1 L 1 , -M 2 L 2 ≤ j ≤ M 2 L 2 , 0 ≤ k ≤ M 3 L 3 } .
For the implementation, we simulate N independent Brownian motions as follows :

W t l+1 -W t l ∼ N (0, h t ).
Then, the price path is given by

P 0 t l+1 = P 0 t l e (b-σ 2 2 )ht+σ(Wt l+1 -Wt l )
For the approximation of the value function v n at the point (t, Z 0 t ) where Z 0 t is the random vector (X 0 t , y, Z 0 t ) (the randomness is only in the third component of this vector), t ∈ T m , two cases are possible :

Case 1: If Z t ∈ [-M 1 , M 1 ] × [-M 2 , M 2 ] × [-M 3 , M 3 ], then Ẑ0 t = N (Ω hz ) i=1 z i 1 A i (Z 0 t ),
where

N (Ω hz ) := Card{z s.t. z ∈ Ω hz } and (A i ) 1≤i≤N (Ω hz ) is

a Borel partition of S defined by

A i = z ∈ S s.t. |z i -z| = min 1≤j≤N (Ω hz ) |z j -z| .
|.| denotes the canonical Euclidean norm, and we take v n (t,

Z 0 t ) ≈ v n (t, Ẑ0 t ). Case 2: If Z t / ∈ [-M 1 , M 1 ] × [-M 2 , M 2 ] × [-M 3 , M 3 ], |Z 0 t -Ẑ0 t | could be large. To approximate v n (t, Z 0
t ), we use the growth condition of the value function

v(t, z) ≤ e ρ(T -t) γ (x + p λ (1 -e -λy )) γ (5.11)
where ρ is a positive constant s.

t. ρ > γ 1-γ b 2 +r 2 +σ 2 r(1-γ) σ 2
(See Propopsition 4.1 in Ly Vath, Mnif and Pham [START_REF] Vath | A model of optimal portfolio selection under liquidity risk and price impact[END_REF]). The approximation of v n (t, Z 0 t ) is given by

v n (t, Z 0 t ) ≈ v n (t, Ẑ0 t ) X 0 t + P 0 t (1-e -λy ) λ γ X0 t + P 0 t (1-e -λy ) λ γ (5.12)
The discrete time approximation for the value function v n is given by :

v n+1 (t, z) = sup τ ∈S m,t,T E e -r(T -τ ) Hv n (τ, X0,t,x τ , y, P 0,t,p τ ) , (t, z) ∈ T m × Ω hz .
The Snell envelop is computed by backward induction :

v n+1 (t m , z) = Hv n (t m , z)
and

v n+1 (t l-1 , z) = max Hv n (t l-1 , z); e -rht E[v n+1 (t l , Z 0 l )|F t l-1 ] , 1 ≤ l ≤ m,
where F t l-1 = σ(P t j , j ≤ l -1) is the discrete-time filtration. Hence :

E[v n+1 (t l , Z 0 l )|F t l-1 ] = E[v n+1 (t l , Z 0 l )|P t l-1 ] =: ρ(t l-1 , P 0 t l-1 ), 1 ≤ l ≤ m.

Estimation of the conditional expectation using Malliavin Calculus

Here, we are interested in computing the conditional expectation E[v n (t + h, Z 0 t+h )|P t ]. From the definition of v n , we have v n ≤ v. The main idea of the Malliavin method consists in using the Malliavin integration by part formula in order to get rid of the Dirac point masses in the following expression :

E[v n (t + h, Z 0 t+h )|P t = p] = E[v n (t + h, Z 0 t+h )δ p (P 0 t )] E[δ p (P 0 t )] . (5.13)
We focus on the calculation of E[v n (t+h, Z 0 t+h )δ p (P 0 t )]. We recall that P 0 t = p 0 e (b-σ 2 2 )t+σWt . We now define vn,t+h,x,y (Br) := v n (t + h, x, y, e (b-σ 2

2 )(t+h)+σBr ), and

pt := 1 σ (ln p p 0 -(b - σ 2 2 )t).
We obtain :

E[v n (t + h, Z 0 t+h )δ p (P 0 t )] = E[v n,t+h,X 0 t+h ,y (W t+h )δ pt (W t )]
. By the independence of Brownian motion's increments, we have :

E[v n,t+h,X t+h ,y (W t+h )δ pt (W t )] = vn,t+h,X t+h ,y (w 1 + w 2 )δ pt (w 1 )ϕ( w 1 √ t )ϕ( w 2 √ h ) dw 1 √ t dw 2 √ h ,
where ϕ is the density of standard one dimensional normal distribution. Using the growth condition of the value function v (5.11) we obtain

E[|v n (t + h, Z 0 t+h )| 2 ] < ∞.
(5.14)

Recalling that δ x (w 1 ) is a derivative of 1 w 1 ≥x , using (5.14) and by integration by parts formula with respect to w 1 variable and then with respect to variable w 2 , we get :

E[v n,t+h,X t+h ,y (W t+h )δ pt (W t )] = E vn,t+h,X t+h ,y (W t+h )1 [pt,∞) (W t ) - W t t + W t+h -W t h .
(5.15)

By denoting A h := Wt t - W t+h -Wt h
, it follows that :

E[v n (t + h, Z t+h )δ p (P t )] = E v n (t + h, Z t+h )1 [p,∞) (P t )A h .

Variance reduction by localization

By using Monte Carlo Method, we recover a convergence rate of the order √ N for the conditional expectation estimator where N is the simulation number. However, the variance of the estimator explodes as h tends to zero since lim sup

h-→0 A h = ∞ and lim inf h-→0 A h = -∞.
To find a remedy to this problem, we introduce localizing functions. Such functions catch the idea that the relevant information for the computation of E[g(S t+h )|S t = x] is located in the neighborhood of x. Let ϕ be an arbitrary localizing function. By definition, ϕ is smooth, bounded and it satisfies ϕ(0) = 1. Recalling the same arguments as in (5.15) and using (5.14), we obtain a family of alternative representations of the conditional expectation given by (5.13) :

E[v n (t + h, Z 0 t+h )δ p (P 0 t )] = E [v n (t + h, Z t+h )δ pt (W t )ϕ(W t -pt )] = E 1 Wt>pt v n (t + h, Z 0 t+h )(ϕ(W t -pt )A h -ϕ (W t -pt )) .
Moreover, it is possible to reduce the Monte Carlo estimator variance by a convenient choice of the localizing function. We consider the integrated mean square error :

J(ϕ) := IR E 1 Wt>pt v 2 n (t + h, Z 0 t+h )A 2 h,ϕ dx, (5.16)
where we adopted the following notation : A h,ϕ := ϕ(W t -pt )A h -ϕ (W t -pt ) and we are interested in minimizing J respect to the subset { ϕ smooth, bounded and ϕ(0) = 1}. Following [START_REF] Bouchard | On the Malliavin approach to Monte Carlo approximation of conditional expectations[END_REF], we prove that the optimal localizing function is given by :

ϕ(x) = e ν h x where ν h := E[v 2 n (t + h, Z 0 t+h )A 2 h ] E[v 2 n (t + h, Z 0 t+h )] 1 2
.

In conclusion, we obtain

E[v n (t + h, Z t+h )δ p (P t )] = E v n (t + h, Z t+h )1 [p,∞) e ν h (Wt-pt) (A h -ν h ) .

Algorithm and discrete value function formula

The algorithm computes two sequences {v n , ζ n } n≥1 by performing the following steps. Parameters: , λ, k, L 1 , L 2 , L 3 , N the number simulation, T , M 1 , M 2 , M 3 and m. Initialisation: v 0 = (v 0 (t, z)) (t,z)∈Tm×Ω hz , n = 0.

Step 1: Compute Hv n and ζ n on T m × Ω hz defined by

Hv n (t, z) = sup ζ∈ Ĉ(z) v n (t, Γ(z, ζ)), (t, z) ∈ T m × Ω hz , where Γ(z, ζ) = (x, ŷ, p) = N (Ω hz ) i=1 z i 1 A i (x -ζpe λy -k, y + ζ, pe λζ ), z i ∈ Ω hz and Ĉ(z) = {ζ ∈ IR s.t. L( Γ(z, ζ)) := max L0 (z), L 1 (z) 1 y≥0 + L0 (z)1 y<0 ≥ 0}, L0 (z) is the closest point of the grid (ih x ) -M 1 L 1 ≤i≤M 1 L 1 to the point x -l(y, p) -k.
Step 2: According to the previous section, we are able to calculate the value function :

v n+1 (t l , z) = max Hv n (t l , z); e -rht ρn (t l , z) , 0 ≤ l ≤ m -1, z ∈ Ω hz ,
Let us denote P (i) the i-th price simulation such that 1 ≤ i ≤ N , where N is the simulation number. Then, we define the estimators of ρ n by :

ρn (t l , z) = 1 N N i=1 v n (t l+1 , Z 0(i) t l+1 )1 [p,∞) e ν h (W (i) t l+1 -p (i) t l+1 ) (A (i) h -ν h ) 1 N N i=1 1 [p,∞) e ν h (W (i) t l+1 -p (i) t l+1 ) (A (i) h -ν h )
, where z = (x, y, p), A

h := W (i) t l t l - W (i) t l+1 -W (i) t l+1 h , pi t l = 1 σ (ln p p 0 -(b -σ (i)
2 )t l ) and W (i) i-th simulation of W . Taking into account the growth condition of the value function, we truncate the estimator ρn :

ρn (t l , z) := ρn (t l , z) ∧ 1 γ e ρ(T -t l ) (x + p λ ) γ ,
which improves the algorithm.

Step 3: Stopping test: If ||v n+1 -v n || ∞ ≤ , stop, otherwise go to step 1.

Numerical results

The computation is achieved with a cluster of 13 Intel Xeon Processors running at 2.8 Ghz with 2 Giga Bytes of RAM. Numerical tests are performed with the following numerical constants γ = 0.5, r = 0.1, α = 0.12, σ = 0.3. According to Bouchard and Touzi [START_REF] Bouchard | Discrete time approximation and Monte Carlo simulation of Backward Stochastic Differential Equations[END_REF], in order to achieve an error estimate of the order of n -1/2 , we have to choose a number N of simulated trajectories such that N = O(n 7/2 ). We choose N = 10000. Contrary to the policy-iteration algorithm (named Howard algorithm), the value-iteration algorithm needs more iterations to converge. At each iteration and each step of time we must estimate 1267 conditional expectations (it is the number of points in our domain). This explains that we need several days to achieve the whole computation. We equally mention that by using the probabilistic approach, we do not need to assume any boundary condition as in Chancelier, Øksendal, and Sulem [START_REF] Chancelier | Combined stochastic control and optimal stopping, and application to numerical approximation of combined stochastic and impulse control[END_REF]. However when the trajectories are outside our bounded domain, we approximate the value function by taking into account the growth condition (see (5.12)).

A partition of the solvency region S is displayed in figures ( 1)-( 2)-( 3) for different values of P and λ. It consists of three regions: Buy (B), Sell (S), and No-Trade (NT) regions. The domain between R1 and R2 corresponds to the region reached by the state variable after a purchase or a sale of risky asset, dictated by the optimal strategy. Due to the presence of fixed costs, the lines R1 and R2 do not coincide with D1 and D2 boundaries of the no-transaction region.

We equally try to see the sensitivity of different parameters and variables. First, there is a reduction in the No-Trade region when the price of the risky asset P increases, i.e. the line D1 moves downwards while the line D2 marginally moves upwards (see Figures (1)-( 2)). The interpretation of this observation is the following :

• in the case where the investor has a significant long position in the risky asset, he is required to reduce his risky asset position when the share price goes up. This phenomenon has also been observed in the Merton model [START_REF] Merton | Optimum consumption and portfolio rules in a continuous-time model[END_REF].

• in the case where the investor has a significant short position in the risky asset, he is required to buy back shares in order to reduce the risk when the share price goes up.

Second, we look at the impact of the coefficient of the impact price λ. We notice that when λ increases, the NT region widens (see Figures (2)-( 3)). In particular, the line D1 significantly moves upwards. Economically, it means that when the liquidity impact increases, the investor should trade less frequently.

Figure 2 : 1 CFigure 1 :

 211 Figure 2: The solvency region when p < k e Figure 3: The solvency region when p > k e

. 3 )

 3 (ii) A locally bounded function u on [0, T ) × S is a constrained viscosity solution of (3.1) in [0, T ) × S if u is a viscosity subsolution of (3.1) in [0, T ) × S and a viscosity supersolution of (3.1) in [0, T ) × S. In Ly Vath, Mnif and Pham [10], we obtained the following characterization Theorem 3.1 The value function v is continuous on [0, T )×S and is the unique (in [0, T )× S) constrained viscosity solution to (3.1) satisfying the boundary and terminal condition :

Figure 1 :Figure 2 :

 12 Figure 1: The optimal transaction policy for p=2, λ=0.5 and k=1

Figure 3 :

 3 Figure 3: The optimal transaction policy for p=3, λ=0.6 and k=1

  .6) When a discrete trading ∆Y s = ζ n+1 occurs at time s = τ n+1 with pretrade price P s -= P τ - n+1, we assume that in addition to the amount of stocks θ(∆Y s , P s -) = θ(ζ n+1 , P τ -

	n+1


Keywords: Glideslope approach, impulsive control, circular rendezvous, Hill-Clohessy-Wiltshire equations, linear programming

This paper presents minimum-fuel glideslope autonomous guidance algorithms for approaching a target evolving on a circular orbit. In the context of a chemical propulsion, the classical multipulse glideslope algorithm of Hablani is revisited and it is shown that when considering specific common directions for the glideslope such as V-bar and R-bar directions, a linear formulation of the circular minimum-fuel linearized rendezvous problem may be deduced. Unlike the classical glideslope algorithm for which there is no direct control on the fuel consumption, additional degrees of freedom and relevant decision variables may be identified by combining an analytical expression for the maximal guidance error and Hill-Clohessy-Wiltshire relative equations of motion. For a fixed-time rendezvous with a pre-assigned number of maneuvers, a fuel-optimal solution with a bounded guidance error is obtained by solving a linear programming problem. Numerical examples demonstrate the usefulness of the approach with respect to the classical ones when the approach corridor has to fulfill stringent geometrical restrictions such as line-of-sight constraints.

V-bar and R-bar Glideslope Guidance Algorithms for Fixed-Time Rendezvous:

A Linear Programming Approach Yassine Ariba * , * * , Denis Arzelier * * , Laura Sofia Urbina * * , Christophe Louembet * * * ,

INTRODUCTION

During the last fifteen years, there has been an increasing demand to perform autonomous rendezvous, proximity or On-Orbit Servicing operations between an active chaser spacecraft and a passive target spacecraft. Autonomy means that new challenges have to be met when designing appropriate guidance schemes. In particular, simplicity of onboard implementation while preserving optimality in terms of fuel consumption, is fundamental. For the proximal rendezvous phase, the assumption that the distance between the chaser and the target is small enough compared to the distance between the target and the center of attraction a linearized model for relative dynamics may be used to simplify the design of guidance schemes. In addition, if the maneuvers carried out on-board provide short high thrust pulses approximated as impulsive maneuvers, the design of efficient guidance scheme for the general rendezvous problem may be greatly simplified, Fehse (2003). The impulsive approximation for the thrust means that instantaneous velocity increments are applied to the chaser whereas its position remains unchanged.

This paper focuses on the fixed-time linearized minimumfuel impulsive proximal rendezvous problem. Depending on various operational and safety constraints, different approaching trajectories may be envisioned and proposed to realize the proximal rendezvous. In the literature, common approach strategies and directions of approaching a target in the close range phase of the rendezvous mission are known as V-bar (curvilinear orbit direction as a straight line) and R-bar(direction of the center of at-traction) approaches. This is mainly due to observability (line-of-sight constraints) and safety reasons imposing the requirement of a trajectory belonging to a cone-shaped approach corridor. One simple reference approaching trajectory complying to safety restrictions is known as the glideslope approach. A glideslope trajectory is a straight path connecting the current location of the chaser to its final destination which is a position of interest in space near the target. The glideslope approach has been first defined in the past for rendezvous and proximity operations involving the space shuttle, [START_REF] Pearson | The glideslope approach[END_REF]. This preliminary study has been extended and generalized later for any direction in space and circular reference orbit in [START_REF] Hablani | Guidance and relative navigation for autonomous rendezvous in a circular orbit[END_REF], [START_REF] Wang | Guidance algorithms for the near-distance rendezvous of on-orbitservicing spacecraft[END_REF] and for elliptic reference orbit in [START_REF] Okasha | Guidance, navigation and control for satellite proximity operations using tschauner-hempel equations[END_REF]. Indeed, the results presented in [START_REF] Hablani | Guidance and relative navigation for autonomous rendezvous in a circular orbit[END_REF] are well-known and define the so-called classical glideslope algorithm.

Our goal is to revisit the classical glideslope algorithm of Hablani in specific cases (V-bar and R-bar approaches) for which additional degrees of freedom helping to derive an optimization formulation may be defined. The inherent guidance error produced by the glideslope approach at each hump is bounded allowing to control it in a linear programming formulation of the guidance problem. Thanks to simple analytical derivations, two algorithms based on a linear programming formulation are designed and permit to obtain a minimum-fuel solution to the glideslope guidance problem while controlling the guidance error. The two proposed algorithms are compared to the classical glideslope algorithm and to analytical transfers from Fehse (2003) on two numerical examples (one for the V-bar and one for the R-bar) illustrating the interest of this new approach. Note that the trajectories are considered in an open-loop setup in this paper whereas real trajectories undergo orbital disturbances and therefore, have to be closed-loop controlled.

Notations: A ⊗ B stands for the usual Kronecker product between matrices A and B. O p×m , I m and 1 m denote respectively the null matrix of dimensions p×m, the identity matrix of dimension m and the m-vector composed of ones. M m and L m are respectively defined as a (m+1)×m matrix and a m square matrix defined by:

M m =         1 0 • • • 0 -1 1 . . . . . . 0 -1 . . . 0 . . . . . . . . . 1 0 • • • 0 -1         , L m =         0 0 • • • • • • 0 1 0 . . . . . . 0 . . . . . . . . . 0 . . . . . . . . . . . . 0 0 • • • 0 1 0         (1) ω = µ a 3
is the mean motion of the circular reference orbit where a is its radius and µ is the standard gravitational parameter.

CLASSICAL GLIDESLOPE APPROACH FOR RENDEZVOUS

Relative motion dynamics

This work focuses on the close range phase of the spacecraft rendezvous mission, Fehse (2003). When relative navigation information is available to the chaser, the relative motion expressed in the Local-Vertical-Local-Horizontal (LVLH) frame centered at the target spacecraft position and depicted on Figure 1 is used to characterize the dynamics of the chaser. The vector r defines the relative position of the chaser with respect to the target. Without any further assumptions, the relative dynamics are described by a set of nonlinear differential equations, Alfriend (2010), which are hard to handle for control purposes. Considering that the two spacecraft are sufficiently close to each other (proximity assumption) and under Keplerian hypothesis, linear relative dynamics may be derived as shown in the seminal work of [START_REF] Lawden | Optimal trajectories for space navigation[END_REF]. This set of equations are also known as the Tschauner-Hempel equations, [START_REF] Tschauner | Optimale beschleunigungs-programme fur des rendezvous manover[END_REF], that boil down to the well known Hill-Clohessy-Wiltshire equations, [START_REF] Clohessy | Terminal guidance system for satellite rendezvous[END_REF], when the target spacecraft orbit is supposed to be circular. Defining the state vector of the relative dynamics by X(ν) = [r(ν), v(ν)] T = [x, y, z, dx dν , dy dν , dz dν ] T where time is replaced by the true anomaly ν as the independent variable, the autonomous relative motion of the chaser has the following Linear Time-Invariant state space representation, Alfriend (2010):

dX dν (ν) =      0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 2 0 -1 0 0 0 0 0 0 3 -2 0 0      X(ν). (2)
The transition matrix Φ(ν, ν 0 ) of ( 2) is readily available such that X(ν) = Φ(ν, ν 0 )X(ν 0 ) for ν ≥ ν 0 where:

Φ(ν, ν 0 ) =     1 0 6(δν -sν ) (4sν -3δν)/ω 0 2(1 -cν )/ω 0 cν 0 0 sν /ω 0 0 0 4 -3cν 2(cν -1)/ω 0 sν /ω 0 0 6ω(1 -cν ) 4cν -3 0 2sν 0 -ωsν 0 0 cν 0 0 0 3ωsν -2sν 0 cν     (3)
where δν = ν -ν 0 , c ν = cos(ν -ν 0 ) and s ν = sin(ν -ν 0 ). It is assumed that only the chaser is cooperative using 6 ungimbaled identical chemical thrusters. The use of chemical propulsion allows to idealize possible thrusts as impulsive maneuvers providing instantaneous velocity changes while the relative position remains unchanged during firing. Thus, a controlled trajectory composed of N impulses is described by the following equation:

X(ν) = Φ(ν, ν 0 )X(ν 0 ) + N i=1 Φ(ν, ν i )B∆v i , (4)
where

ν 1 < ν 2 < • • • < ν N ≤ ν and ∆v i denotes the impulsive control applied at ν i . B = [O 3 I 3 ] T is the input matrix.
It is well-known and it may be observed from equations ( 2) and ( 3) that the in-plane motion (x -z) and the out-of-track motion are decoupled. This work will be focused on the in-plane motion for which fueloptimal V-bar and R-bar glideslope guidance schemes will be derived. Hereafter, the following notation describing the in-plane free motion is adopted:

r(ν) v(ν) = Φ rr (ν, ν 0 ) Φ rv (ν, ν 0 ) Φ vr (ν, ν 0 ) Φ vv (ν, ν 0 ) r 0 v 0 =    1 6(δν -s ν ) (4s ν -3δν)/ω 2(1 -c ν )/ω 0 4 -3c ν 2(c ν -1)/ω s ν /ω 0 6ω(1 -c ν ) 4c ν -3 2s ν 0 3ωs ν -2s ν c ν    (5) where r = [x z] T and v = [v x v z ] T .

Hablani's classical glideslope approach for rendezvous

As mentioned in the introduction, the most cited reference on this topic is the paper [START_REF] Hablani | Guidance and relative navigation for autonomous rendezvous in a circular orbit[END_REF] in which the so called classical inbound and outbound glideslope approaches for circular reference are presented in a general setup. This paper is only concerned by inbound decelerating glideslope trajectories and the main features of the classical glideslope algorithm are now recalled in this particular case. In the general case, this guidance trajectory is characterized by a straight line and its associated vector ρ(ν) = r c (ν) -r T , defining the commanded path as illustrated in Figure 2. Defining ρ 0 = r 0 -r T , the unit vector u defines the direction of the straight path where:

u = x T -x 0 ρ 0 z T -z 0 ρ 0 T .
The chaser is commanded to reach r T from r 0 following a specific commanded profile for the distance to go ρ 0 = ρ 0 :

x z T arget Chaser r T r0 u A r r c ρ O Fig. 2. Glideslope approach. ρ(ν) = ρ 0 e ( ρ0 -ρT )ν ωρ 0 + ρT ρ 0 ρ0 -ρT   e ( ρ0 -ρT )ν ωρ 0 -1   
(6) where ρ0 < 0, ρT < 0 are respectively the initial and final commanded velocities and ρ 0 is the initial distance to go. These quantities are specified by the designer and inputs for the classical algorithm. Note that ρ0 > ρT and that ρ T = 0 by definition of the trajectory. For a given set of these parameters, the basic principle of the classical algorithm is then to analytically compute a fixed number of impulses equally spaced in time over the transfer duration T . Each computed incremental velocity at r k is obtained as ∆v

k = v + k -v - k where v + k
is the departure velocity needed to go from r k to r k+1 and v - k is the arrival velocity at r k . Both quantities are simply obtained by solving the autonomous Hill-Clohessy-Wiltshire (5) at each step k.

The classical glideslope algorithm is straightforward and rapid to implement but suffers from key shortcomings. In particular, it is important to notice that the actual trajectory of the chaser will not be strictly along the commanded straight line path but will exhibit humps between the N points where an impulsive maneuver is performed and located on the commanded path (cf. Figure 2). These humps coming from the natural relative motion of the chaser driven by the Hill-Clohessy-Wiltshire equations are nothing but lateral guidance position errors possibly occurring in and out-of-plane that cannot be directly controlled in the classical glideslope algorithm. In addition, if the initial and final commanded velocities of the glideslope profile are a priori given, there is no degree of freedom left to control the transfer time and the consumption when X(0) and X(ν T ) are fixed. Indeed, the transfer time T is not fixed a priori but deduced from the initial and final commanded velocities ( ρ0 , ρT ) and from the initial distance to go ρ 0 :

T = ρ 0 ρ0 -ρT ln ρT ρ0 . ( 7
)
The consumption itself is computed a posteriori from the velocity increments without any possibility to optimize it for given side conditions of the rendezvous. The objective of the paper is then to propose a new optimization algorithm for particular and common glideslope approaching directions as V-bar and R-bar allowing to take these important features into account.

-Minimize the fuel-consumption for a given set of initial and final rendezvous conditions and an a priori fixed time of transfer; -Control the maximum guidance error by defining constraints on the humps profile.

Derivation of guidance error in orbital plane

Keeping in mind the requirements of the previous section, finding an analytical expression for the guidance error is a key technical step for the development of an efficient algorithm. Due to the decoupling property between inplane and out-of-plane keplerian relative motions, this study will be limited to glideslope approaches defined in the orbital plane.

Let us define the point A in the x -z plane (see Figure 2) as the orthogonal projection of the chaser position r on the glideslope straight line. We can thus write OA = (r -r 0 ) T u u + r 0 , and hence, guidance error is given by

T = r -OA.
The guidance error norm is then deduced as:

T 2 = δr -[δr T u] u T δr -[δr T u] u = δr T δr -[δr T u] 2 = δr 2 x (1 -u 2 x ) + δr 2 z (1 -u 2 z ) -2u x u z δr x δr z (8)
with δr = rr 0 , δr x = x -x 0 and δr z = z -z 0 .

MINIMUM-FUEL V-BAR GLIDESLOPE APPROACH VIA LINEAR PROGRAMMING

Constraint on guidance error

Considering a glideslope along a line parallel to the xaxis, the unitary direction vector u is [1 0] T . Hence, the guidance error norm ( 8) is reduced to

T 2 = δr 2 z , = (4 -3 cos δν)z 0 + 2 ω (cos δν -1)x 0 + 1 ω z 0 sin δν -z 0 2 , = 3(1 -cos δν)z 0 + 2 ω (cos δν -1)x 0 + 1 ω z 0 sin δν 2 . ( 9
)
When the guidance error distance is maximal, the velocity vector v is parallel to the glideslope line. A maximum condition is then given by z = 0 ⇔ 3ωz 0 sin δν -2x 0 sin δν + z 0 cos δν = 0. (10) Lemma 1. Let us define δν m the value of the true anomaly when guidance error is maximal and δν 0 the value of the true anomaly when guidance error is null, that is, when the satellite is back on the glideslope commanded path. Thus, δν 0 is the true anomaly interval for a hump. Considering a single hump, we have the relationship: δν 0 = 2 δν m . That is, the maximum occurs in the middle (in terms of true anomaly) of the hump.

From equation ( 9), guidance error distance can be readily upper-bounded with

| T | = (3 -3 cos δν)z 0 + 2 ω (cos δν -1)x 0 + 1 ω z 0 sin δν , ≤ (3 -3 cos δνm)z 0 + 2 ω (cos δνm -1)x 0 + 1 ω z 0 sin δνm .

Problem formulation

We aim now at combining motion equations, guidance error equations and identify relevant decision variables in order to formulate an optimization problem. The analysis of the previous paragraph provides an inequality that will be used to specify maximum conditions on guidance error for each hump, k = 0, . . . , N -1:

(3 -3 cos δνm)z 0 + 2 ω (cos δνm -1)x k + 1 ω z k sin δνm ≤ m k are equivalently    (3 -3 cos δνm)z 0 + 2 ω (cos δνm -1)x k + 1 ω z k sin δνm ≤ m k -(3 -3 cos δνm)z 0 - 2 ω (cos δνm -1)x k - 1 ω z k sin δνm ≤ m k (11)
where m k is given and specifies the maximal allowable guidance error during the (k + 1) th maneuver. We recall that assuming the interval δν 0 is constant, the parameter δν m is constant too.

After each maneuver, the chaser must be back on the glideslope. Such a requirement can be enforced by N equations of the form r k+1 = Φ rr r k + Φ rv v k , k = 0, . . . , N -1 ( 12) where

r k = x k z 0 , v k = x k z k , r 0 = x 0 z 0 , r N = x T z 0 .
r 0 and r T are given and impose the initial and final positions. The z-component of r k remains constant and equals z 0 since the glideslope is assumed to be parallel to the V-bar direction. The interval δν 0 between two impulses being constant, the transition matrix {Φ rr , Φ rv , Φ vr , Φ vv } is constant. We consider v k as the velocity vector resulting from the (k + 1) th control impulse. Its components together with x k form our decision variables. Hence, all positions on the glideslope, except the initial and final ones, are set free. Moreover, an additional equation may be specified so as to impose a final velocity condition v T :

v T = v N = Φ vr r N -1 + Φ vv v N -1 + ∆v N .
(13) v N is thus the desired final velocity vector (and not a variable) while ∆v N is an extra variable representing the last impulse.

The optimization problem consists in finding the values of

• velocity vectors v k for k = 0, . . . , N -1, • positions x k between two manoeuvres for k = 1, . . . , N -1, • the final impulse ∆v N , minimizing the overall consumption. This latter is based on impulses:

min N k=0 ∆v k 1 , with ∆v k = v k -v k-.
This criterion can be expressed with respect to the decision variables: (14) where v 0-is the initial velocity vector. This cost function involving absolute values can be transformed into a linear function with the introduction of new variables and inequality constraints, [START_REF] Bertsimas | Introduction to Linear Optimization[END_REF].

min v 0 -v 0-1 + N -1 k=1 v k -Φ vr r k-1 -Φ vv v k-1 1 + ∆v N 1 ,

Minimum-fuel LP algorithm

In this paragraph, the optimization to be addressed is stated and an algorithm is proposed. This algorithm computes appropriate impulses to transfer the satellite from an initial position r 0 to a final position r T while minimizing the consumption. Considering the decision variables previously defined, the cost criterion ( 14) and gathering constraints on guidance errors (11), on positions (12) and on the final velocity (13), a linear programming problem can be built. min c T X s.t. AX ≤ b, A eq X = b eq (15)

The vector X puts together all the decision variables,

X = x 1 . . . x N -1 v T 0 . . . v T N -1 ∆v T N d T 0 . . . d T N T , (16
) where d i ∈ R 2 are extra variables used to deal with the absolute values in the cost function. This cost criterion is then reduced to the sum of d i , i = 0, . . . , N , and the vector c is defined as

c T = [ O 1×3N +1 1 . . . 1 ] .
The equality constraint is built on equations ( 12) and ( 13), while the inequality constraint is built on equation ( 11) and additional inequalities introduced to cope with the absolute values in the cost function, [START_REF] Bertsimas | Introduction to Linear Optimization[END_REF]. The different matrices involved in the linear programming problem are defined as:

Aeq = M N -1 ⊗ C1 -I N ⊗ Φrv O 2N ×2 O 2N ×(2N +2) O 2×(N -1) O 2×(2N -2) Φvv I2×2 O 2×(2N +2) (17) and beq =   Φrrr 0 -C 2 z 0 1 N -2 ⊗ (Φrr -I 2 )C 2 z 0 ΦrrC 2 z 0 -r T v T -ΦvrC 2 z 0   , (18)
where 19) where

C 1 = 1 0 and C 2 = 0 1 . A =    O 4(N +1)×N -1 A O 4N ×2 I 2 -I 2 I N +1 ⊗ -I 2 -I 2 O 2N ×N -1 I N ⊗ M O 2N ×2 O 2N ×2(N +1)    (
A =   L N ⊗ -Φvv Φvv + I N ⊗ I 2 -I 2 O 2×2N O 2×2N   and b =              v 0 - -v 0 - Φvrr 0 -Φvrr 0 1 N -2 ⊗ ΦvrC 2 z 0 -ΦvrC 2 z 0 O 2×1 O 2×1 M 0 . . . M N -1              (20)
with

M =   2 ω cos δνm -1 1 ω sin δνm - 2 ω cos δνm -1 - 1 ω sin δνm   , M k = m k -3 -3 cos δνm z 0 m k + 3 -3 cos δνm z 0 . (21)
Remark 1. Note that expressions of above matrices have been slightly simplified because we have Φ rr C 1 = C 1 and Φ vr C 1 = 0 0 .

The procedure to run the glideslope method proposed in this paper is summarized in Algorithm 1. Given an initial {r 0 , v 0-} and a final {r T , v T } configuration, a transfer time T and a number N of maneuvers, the algorithm computes the optimal impulse sequence ∆v k , k = 0, . . . , N , minimizing the consumption while ensuring a maximal guidance error m k for each manoeuvre.

Algorithm 1: Minimum-fuel LP glideslope

Data: r0, r T , v - 0 , v T , N , T , ω, m = [m0, . . . , m N -1 ] Initialization: δt = T /N ; δν0 = ω δt ; δνmax = δν0/2;
Φrr, Φrv, Φvr, Φvv ← compute transition matrix for δt and ω;

for k = 0 to N -1 do M k ← compute matrix with m k , z0 and δνmax;
A ← compute matrix with Φvv, δνmax, ω and N ; b ← compute vector with Φvr, v0 -, M k , r0 and N ; Aeq ← compute matrix with Φrv, Φvv and N ; beq ← compute vector with Φrr, Φvr, r0, r T , v T and N ; c ← compute vector with N ; Solve: the linear programming problem, and store the optimal solution in X if a solution is found then

v k ← extracted from X, k = 0, . . . , N -1; ∆v N ← extracted from X; for k = 0 to N -1 do r k+1 ← Φrrr k + Φrvv k ; v k+1 -← Φvrr k + Φvvv k ; ∆v k ← v k -v k -;
else increase some m k , update M k and b, and re-run the solver;

EXTENSION TO R-BAR APPROACH

Constraint on excursion

In that case, the unitary direction vector u is [0 1] T . Hence, the guidance error norm (8) reduces to:

T 2 = δr 2 x , = 6z 0 (δν -sin δν) + 1 ω (4 sin δν -3δν)x 0 + 2 ω (1 -cos δν)z 0 2 . ( 22
)
The following lemma will be used in the next paragraph to provide an inequality constraint on the excursion for each hump. Unlike the V-bar case, a simple relationship between the true anomaly δν m when the guidance error is maximal and the maneuver interval δν 0 cannot be established. So, a conservative bound is introduced. Lemma 2. Let us define δν 0 the value of the true anomaly when guidance error is 0, that is, when the satellite is back on the glideslope commanded path. Considering a single hump such that δν 0 ∈ [0 arccos 3 4 ], a conservative upperbound of the guidance error is given by:

| T | ≤ α 1 |x 0 | + α 2 |z 0 | + α 3 |z 0 |, with α 1 = | 1 ω (4 sin δν 0 -3δν 0 )|, α 2 = | 2 ω (1 -cos δν 0 )| and α 3 = |6(δν 0 -sin δν 0 )|.
Although in the above lemma the range of δν 0 is restricted to [0 arccos 3 4 ], this result can be easily extended to a full revolution.

Formulation and algorithm

The problem formulation for the R-bar approach follows the same process as in Subsection 3.3. The set of N motion equations ( 12) with

r k = x 0 z k , v k = x k z k , r 0 = x 0 z 0 , r N = x 0 z T .
allows the satellite to go along the glideslope line. For the R-bar approach, the z-component is changing while the xcomponent remains constant. Thus, equations ( 12), ( 13) and ( 14) are considered again and combined with the new appropriate vector of decision variables:

X = z 1 . . . z N -1 v T 0 . . . v T N -1 ∆v T N d T 0 . . . d T N T .
Again, d i ∈ R 2 are extra variables used to deal with the absolute values in the cost function. Eventually, a minimum-fuel R-bar glideslope approach is formulated as a linear programming problem (15) with matrices defined as follows.

Aeq = Ã -I N ⊗ Φrv O 2N ×2 O 2N ×2N +2 O 2×N -2 ΦvrC 2 O 2×2N -2 Φvv I 2×2 O 2×2N +2 , (23)
with

à = I N -1 ⊗ C 2 O 2×N -1 + O 2×N -1 I N -1 ⊗ -ΦrrC 2 and beq =   Φrrr 0 -C 1 x 0 O 2(N -2)×1 C 1 x 0 -r T v T   , (24)
A =     Â A O 4N ×2 I 2 -I 2 I N +1 ⊗ -I 2 -I 2 O 8×N -1 I N -1 ⊗ M 1 I N ⊗ M 2 O 8N ×2 O 8N ×2(N +1)     (25)
where

A =   L N ⊗ -Φvv Φvv + I N ⊗ I 2 -I 2 O 2×2N O 2×2N   , Â =   O 8×N -2 O 8×1 I N -2 ⊗ -ΦvrC 2 ΦvrC 2 O 4(N -2)×1 O 4×N -2 O 4×1   , M 1 = 1 4 -1 4 α 3 , M 2 = 1 2 ⊗   α 1 α 2 α 1 -α 2 -α 1 α 2 -α 1 -α 2   , and b =           v 0 - -v 0 - Φvrr 0 -Φvrr 0 O 4(N -1)×1 1 8 m 0 -M 1 z 0 1 8 m 1 . . . 1 8 m N -1           . (26)

NUMERICAL EXAMPLES

This scenario, extracted from [START_REF] Hablani | Guidance and relative navigation for autonomous rendezvous in a circular orbit[END_REF], is a tangential impulse transfer along V-bar direction and considers a target flying an orbit at 400km altitude, which corresponds to an orbital rate of ω ≈ 0.001 rad/s. The standard glideslope algorithm is used to transfer the chaser from x = -500 m to x = -100 m in 10 maneuvers (N ) along 9 min (T = 540 s). This example's parameters are set as follows:

r 0 = -500 0 -20 , r T = -100 0 -20 , v 0 -= 0 0 0 .
These parameters have been injected in Algorithm 1 to compute the 10 pulses, and a final velocity has been specified as well (v T = 0). Figure 3 shows the chaser trajectories for three different methods in the x -z plane: the standard glideslope defined in [START_REF] Hablani | Guidance and relative navigation for autonomous rendezvous in a circular orbit[END_REF], the two impulses approach from Fehse (2003) and our Algorithm 1. In this simulation, the m k has been chosen identically equal to 1m. It is interesting to note in that configuration that the algorithm adjusts pulses to be regularly spaced. We remind that x k , 0 < k < N , are not constrained (and not ordered) and could be anywhere on the glideslope line. Regarding the consumption, Algorithm 1 outperforms the standard glideslope, reducing the cost from 8.09m/s to 2.31m/s. Note that we have included in the consumption of the standard approach an extra impulse to end up the manoeuvre with the desired final velocity vector v T . The consumption of the two tangential impulses approach is very low, 0.12m/s, but it results in a very large excursion (≈ 125m). Table 1 compares numerical simulations of the standard glideslope and Algorithm 1 for various number of impulses. For every test in Table 1, we set m k identically equal to 20m.

A second scenario, inspired from Di Cairano et al.

(2012), is considered. In this case, from an initial position (x 0 , z 0 ) = (250, 0), the chaser must be transferred to a target platform while remaining within a Line-of-Sight (LOS) cone. The docking port is located at the position coordinates (2.5, 0). The half angle of the LOS cone is 2deg and its vertex is 0.5m behind the docking port. The simulations of the standard glideslope and Algorithm 1, for N = 5 and T = 480s, are shown in Figure 4. Specifying the following upper-bounds on guidance error allows to maintain the spacecraft within the LOS cone. m = [ 5 1.5 0.4 0.1 0.03 ] Moreover, the resulting consumption is 3.87m/s, lower than the standard glideslope 4.90m/s.

CONCLUSION

In this paper, a revisited solution to the classical Hablani glideslope algorithm, [START_REF] Hablani | Guidance and relative navigation for autonomous rendezvous in a circular orbit[END_REF], for the impulsive close range rendezvous in a circular orbit framework is proposed for the the V-bar and R-bar approaches. It is aimed at tackling the two main drawbacks of this classical algorithm which are the lack of control on the bounds over the inherent guidance errors and the impossibility of minimizing the fuel consumption.

In response to these issues, a new problem formulation is developed, where the relative dynamics, the guidance error expressions and the identification of relevant decision variables are combined in order to derive a linear programming formulation of the guidance problem. This new proposed glideslope algorithm includes important design features like the specification of maximum conditions on the guidance error for each hump.

Two simulations are performed and compared to the classical glideslope algorithm given by [START_REF] Hablani | Guidance and relative navigation for autonomous rendezvous in a circular orbit[END_REF] and to the two-impulses approach by (Fehse, 2003). The key features of the exposed results are a significant improved fuel consumption with respect to Hablani's algorithm and a user-defined bound profile on the maximum guidance error, which turns out to be very useful when dealing when visibility constraints.

Fig. 1 .

 1 Fig. 1. LVLH frame for spacecraft rendezvous.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Three glideslope transfers along V-bar.

Table 1 .

 1

		N	2	3	4	10	20
	cons. m/s	Std Alg 1	3.61 2.26	4.76 2.29	5.63 2.30	8.09 2.31	9.32 2.31
	exc. m	Std Alg 1	26.2 13.8	16.2 6.2	11.1 3.5	2.67 0.56	0.78 0.15

ACKNOWLEDGEMENTS

The Authors would like to thank Jean-Claude Berges from CNES and Damiana Losa from Thales Alenia Space for the grants that partly supports this activity.


Jinlong Yang

Jean-Sébastien Franco

Franck Hétroy-Wheeler

Stefanie Wuhrer

Estimation of Human Body Shape in Motion with Wide Clothing

Keywords: human body modeling, shape and motion estimation, statistical shape space

Estimating 3D human body shape in motion from a sequence of unstructured oriented 3D point clouds is important for many applications. We propose the first automatic method to solve this problem that works in the presence of loose clothing. The problem is formulated as an optimization problem that solves for identity and posture parameters in a shape space capturing likely body shape variations. The automation is achieved by leveraging a recent robust pose detection method [1]. To account for clothing, we take advantage of motion cues by encouraging the estimated body shape to be inside the observations. The method is evaluated on a new benchmark containing different subjects, motions, and clothing styles that allows to quantitatively measure the accuracy of body shape estimates. Furthermore, we compare our results to existing methods that require manual input and demonstrate that results of similar visual quality can be obtained.

Introduction

Estimating 3D human body shape in motion is important for applications ranging from virtual change rooms to security. While it is currently possible to effectively track the surface of the clothing of dressed humans in motion [START_REF] Newcombe | Dynamicfusion: Reconstruction and tracking of non-rigid scenes in real-time[END_REF] or to accurately track body shape and posture of humans dressed in tight clothing [START_REF] Bogo | Detailed full-body reconstructions of moving people from monocular rgb-d sequences[END_REF], it remains impossible to automatically estimate the 3D body shape in motion for humans captured in loose clothing.

Given an input motion sequence of raw 3D meshes or oriented point clouds (with unknown correspondence information) showing a dressed person, the goal of this work is to estimate the body shape and motion of this person. Existing techniques to solve this problem are either not designed to work in the presence of loose clothing [START_REF] Weiss | Home 3D body scans from noisy image and range data[END_REF][START_REF] Helten | Personalization and evaluation of a real-time depth-based full body scanner[END_REF] or require manual initialization for the pose [START_REF] Wuhrer | Estimation of human body shape and posture under clothing[END_REF][START_REF] Neophytou | A layered model of human body and garment deformation[END_REF], which limits their use in general scenarios. The reason is that wide clothing leads to strong variations of the acquired surface that is challenging to handle automatically. We propose an automatic framework that allows to estimate the human body shape and motion that is robust to the presence of loose clothing.

Existing methods that estimate human body shape based on an input motion sequence of 3D meshes or oriented point clouds use a shape space that models human body shape variations caused by different identities and postures as prior. Such a prior allows to reduce the search space to likely body shapes and postures. Prior works fall into two lines of work. On the one hand, there are human body shape estimation methods specifically designed to work in the presence of loose clothing [START_REF] Wuhrer | Estimation of human body shape and posture under clothing[END_REF][START_REF] Neophytou | A layered model of human body and garment deformation[END_REF]. These techniques take advantage of the fact that observations of a dressed human in motion provides important cues about the underlying body shape as different parts of the clothing are close to the body shape in different frames. However, these methods require manually placed markers to initialize the posture. On the other hand, there are human body shape estimation methods designed to robustly and automatically compute the shape and posture estimate over time [START_REF] Weiss | Home 3D body scans from noisy image and range data[END_REF][START_REF] Helten | Personalization and evaluation of a real-time depth-based full body scanner[END_REF]. However, these methods use strong priors of the true human body shape to track the posture over time and to fit the shape to the input point cloud, and may therefore fail in the presence of loose clothing.

In this work, we combine the advantages of these two lines of work by proposing an automatic framework that is designed for body shape estimation under loose clothing. Like previous works, our method restricts the shape estimate to likely body shapes and postures, as defined by a shape space. We use a shape space that models variations caused by different identities and variations caused by different postures as linear factors [START_REF] Pishchulin | Building statistical shape spaces for 3d human modeling[END_REF]. This simple model allows for the development of an efficient fitting approach. To develop an automatic method, we employ a robust pose detection method that accounts for different identities [START_REF] Zuffi | The Stitched Puppet: A graphical model of 3d human shape and pose[END_REF] and use the detected pose to guide our model fitting. To account for clothing, we take advantage of motion cues by encouraging the estimated body shape to be located inside the acquired observation at each frame. This constraint, which is expressed as a simple energy that is optimized over all input frames jointly, allows to account for clothing without the need to explicitly detect skin regions on all frames as is the case for previous methods [START_REF] Neophytou | A layered model of human body and garment deformation[END_REF][START_REF] Balan | The naked truth: Estimating body shape under clothing[END_REF].

To the best of our knowledge, existing datasets in this research area do not provide 3D sequences of both body shape as ground truth and dressed scans for estimation. Therefore, visual quality is the only evaluation choice. To quantitatively evaluate our framework and allow for future comparisons, we propose the first dataset consisting of synchronized acquisitions of dense unstructured geometric motion data and sparse motion capture data of 6 subjects with 3 clothing styles (tight, layered, wide) under 3 representative motions, where the capture in tight clothing serves as ground truth body shape.

The main contributions of this work are the following.

-An automatic approach to estimate 3D human body shape in motion in the presence of loose clothing.

-A new benchmark consisting of 6 subjects captured in 3 motions and 3 clothing styles each that allows to quantitatively compare human body shape estimates.

Related work

Many works estimate human posture without aiming to estimate body shape, or track a known body shape over time. As our goal is to simultaneously estimate body shape and motion automatically and in the presence of loose clothing, we will focus our discussion on this scenario.

Statistical shape spaces. To model human body shape variations caused by different identities, postures, and motions, statistical shape spaces are commonly used. These shape spaces represent a single frame of a motion sequence using a low-dimensional parameter space that typically models shape variations caused by different identities and caused by different postures using separate sets of parameters. Such shape spaces can be used as prior when the goal is to predict a likely body shape under loose clothing. Anguelov et al. [START_REF] Anguelov | SCAPE: shape completion and animation of people[END_REF] proposed a statistical shape space called SCAPE that combines an identity model computed by performing principal component analysis (PCA) on a population of 3D models in standard posture with a posture model computed by analyzing near-rigid body parts corresponding to bones. This model performs statistics on triangle transformations, which allows to model non-rigid deformations caused by posture changes. Achieving this accuracy requires solving an optimization problem to reconstruct a 3D mesh from its representation in shape space. To improve the accuracy of the SCAPE space, Chen et al. [START_REF] Chen | Tensor-based human body modeling[END_REF] propose to combine the SCAPE model with localized multilinear models for each body part. To model the correlation of the shape changes caused by identity and posture changes, Hasler et al. [START_REF] Hasler | A statistical model of human pose and body shape[END_REF] perform PCA on a rotation-invariant encoding of the model's triangles. These models may be used as priors when estimating human body shape in motion, but none of them allow to efficiently reconstruct a 3D human model from the shape space.

To speed up the reconstruction time from the SCAPE representation, Jain et al. [START_REF] Jain | MovieReshape: tracking and reshaping of humans in videos[END_REF] propose a simplified SCAPE model, denoted by S-SCAPE in the following, that computes the body shape by performing PCA on the vertex coordinates of a training set in standard posture and combines this with a linear blend skinning (LBS) to model posture changes. Any posture variations present in the training data cause posture variation to be modeled in identity space, which is known to cause counter-intuitive deformations [START_REF] Pishchulin | Building statistical shape spaces for 3d human modeling[END_REF]. To remedy this, recently proposed shape spaces start by normalizing the posture of the training data before performing statistics and model shape changes caused by different factors such as identity and posture as multilinear factors [START_REF] Neophytou | Shape and pose space deformation for subject specific animation[END_REF][START_REF] Pishchulin | Building statistical shape spaces for 3d human modeling[END_REF][START_REF] Loper | SMPL: A skinned multi-person linear model[END_REF]. We use the normalized S-SCAPE model [START_REF] Pishchulin | Building statistical shape spaces for 3d human modeling[END_REF] in this work; however, any of these shape spaces could be used within our framework.

Recently, Pons-Moll et al. [START_REF] Pons-Moll | DYNA: a model of dynamic human shape in motion[END_REF] proposed a statistical model that captures fine-scale dynamic shape variation of the naked body shape. We do not model dynamic geometry in this work, as detailed shape changes are typically not observable under loose clothing.

Estimation of static body shape under clothing. To estimate human body shape based on a static acquisition in loose clothing and in arbitrary posture, the following two approaches have been proposed. Balan et al. [START_REF] Balan | The naked truth: Estimating body shape under clothing[END_REF] use a SCAPE model to estimate the body shape under clothing based on a set of calibrated multi-view images. This work is evaluated on a static dataset of different subjects captured in different postures and clothing styles. Our evaluation on 3D motion sequences of different subjects captured in different motions and clothing styles is inspired by this work. Hasler et al. [START_REF] Hasler | Estimating body shape of dressed humans[END_REF] use a rotation-invariant encoding to estimate the body shape under clothing based on a 3D input scan. While this method leads to accurate results, it cannot easily be extended to motion sequences, as identity and posture parameters are not separated in this encoding.

Both of these methods require manual input for posture initialization. In this work, we propose an automatic method to estimate body shape in motion.

Estimation of body shape in motion. The static techniques have been extended to motion sequences with the help of shape spaces that separate shape changes caused by identity and posture. Several methods have been proposed to fit a SCAPE or S-SCAPE model to Kinect data by fixing the parameters controlling identity over the sequence [START_REF] Weiss | Home 3D body scans from noisy image and range data[END_REF][START_REF] Helten | Personalization and evaluation of a real-time depth-based full body scanner[END_REF]. These methods are not designed to work with clothing, and it is assumed that only tight clothing is present.

Two more recent methods are designed to account for the presence of clothing. The key idea of these methods is to take advantage of temporal motion cues to obtain a better identity estimate than would be possible based on a single frame. Our method also takes advantage of motion cues.

Wuhrer et al. [START_REF] Wuhrer | Estimation of human body shape and posture under clothing[END_REF] use a shape space that learns local information around each vertex to estimate human body shape for a 3D motion sequence. The final identity estimate is obtained by averaging the identity estimates over all frames. While this shape space leads to results of high quality, the fitting is computationally expensive, as the reconstruction of a 3D model from shape space requires solving an optimization problem. Our method uses a simpler shape space while preserving a similar level of accuracy by using an S-SCAPE model that prenormalizes the training shapes with the help of localized information.

Neophytou and Hilton [START_REF] Neophytou | A layered model of human body and garment deformation[END_REF] propose a faster method based on a shape space that models identity and posture as linear factors and learns shape variations on a posture-normalized training database. To constrain the estimate to reliable regions, the method detects areas that are close to the body surface. In contrast, our method constrains the estimate to be located inside the observed clothing at every input frame, which results in an optimization problem that does not require a detection.

Both of these methods require manual input for posture initialization on the first frame. Additionally, a temporal alignment is required by Neophytou and Hilton. Computing temporal alignments is a difficult problem, and manual annotation is tedious when considering larger sets of motion sequences. In contrast, our method is fully automatic and addresses both aspects.

S-SCAPE model

In this work, we use the S-SCAPE model as prior for human body shape changes caused by different identities and postures. While we choose this shape space, any shape space that models identity and posture as multilinear factors could be used [START_REF] Neophytou | Shape and pose space deformation for subject specific animation[END_REF][START_REF] Loper | SMPL: A skinned multi-person linear model[END_REF]. Although such a simple shape space does not accurately model correlated shape changes, such as muscle bulging, it allows to effectively separate the different variations and can be fitted efficiently to input scans. This section briefly reviews the S-SCAPE model introduced by Jain et al. [START_REF] Jain | MovieReshape: tracking and reshaping of humans in videos[END_REF] that allows to separate the influence of parameters controlling identity and parameters controlling posture of a human body shape. In the following, we denote by β and Θ the parameter vectors that influence shape changes caused by identity and posture changes, respectively. In this work, we use the publicly available posture-normalized S-SCAPE model [START_REF] Pishchulin | Building statistical shape spaces for 3d human modeling[END_REF], where each training shape was normalized with the help of localized coordinates [START_REF] Wuhrer | Posture-invariant statistical shape analysis using Laplace operator[END_REF].

In the following, let N v denote the number of vertices on the S-SCAPE model, let s (β, Θ) ∈ R 3N v denote the vector containing the vertex coordinates of identity β in posture Θ, and let s (β, Θ) ∈ R 4N v denote the vector containing the corresponding homogeneous vertex coordinates. For the fixed posture Θ 0 that was used to train the identity space, S-SCAPE models the shape change caused by identity using a PCA model as

s (β, Θ 0 ) = Aβ + µ, (1)
where µ ∈ R 4N v contains the homogeneous coordinates of the mean body shape, A ∈ R 4N v ×d id is the matrix found by PCA, and d id is the dimensionality of the identity shape space. For a fixed identity β 0 , S-SCAPE models the shape change caused by posture using LBS as

s i (β 0 , Θ) = N b j=1 ω ij T j (Θ) s i (β 0 , Θ 0 ) , (2)
where s i and s i denote the standard and homogenous coordinate vector of the i-th vertex of s, N b denotes the number of bones used for LBS, T j (Θ) ∈ R 3×4 denotes the transformation matrix applied to the j-th bone, and ω ij denotes the rigging weight binding the i-th vertex to the j-th bone. Combining Eq. 1 and 2 in matrix notation leads to

s (β, Θ) = T (Θ) Aβ + T (Θ) µ, (3)
where T (Θ) ∈ R 3N v ×4N v is a sparse matrix containing the per-vertex transformations. Using this notation, it is easy to see that S-SCAPE is linear in both β and T (Θ), which allows for a simple optimization w.r.t. β and Θ.

Estimating model parameters for a motion sequence

We start by providing an overview of the proposed method. Fig. 1 shows the different parts of the algorithm visually. Given as input a trained S-SCAPE model and a motion sequence consisting of N f frames F i represented by triangle meshes with unknown correspondence, we aim to compute a single parameter vector β controlling the shape of the identity (as the identity of the person is fixed during motion) along with N f parameter vectors Θ i controlling the postures in each frame, such that

s i (β, Θ i ) is close to F i .
To fit the S-SCAPE model to a single frame F , we aim to minimize

E (F ,β ,Θ )=ω lnd E lnd (F ,β ,Θ )+ω data E data (F ,β ,Θ )+ω cloth E cloth (F ,β ,Θ ) (4)
w.r.t. β and Θ subject to constraints that keep β in the learned probability distribution of parameter values. Here, ω lnd , ω data , and ω cloth are weights that trade off the influence of the different energy terms. The energy E lnd measures the distance between a sparse set of provided landmarks, which correspond to distinctive positions on the human body, to their corresponding locations on s (β, Θ). The provided landmarks are computed automatically in the following.

The energy E data measures the distance between s (β, Θ) and F using a nearest neighbor cost. The energy E cloth is designed to account for loose clothing by encouraging s (β, Θ) to be located inside the observation F . For a motion sequence of N f frames, our goal is then to minimize

E F 1:N f ,β ,Θ 1:N f = N f i=1 E (F i ,β ,Θ i ) (5)
w.r.t. β and Θ 1:N f subject to constraints that keep β in the learned probability distribution of parameter values. Here,

F 1:N f = {F 1 , . . . , F N f } is the set of frames and Θ 1:N f = {Θ 1 , . . . , Θ N f } is the set of posture parameters.
The energy E cloth allows to take advantage of motion cues in this formulation as it encourages the body shape to lie inside all observed frames.

In the following sections, we detail the prior that is used to constrain β as well as the different energy terms. Optimizing Eq. 5 w.r.t. all parameters jointly results in a high-dimensional optimization problem that is inefficient to solve and prone to get stuck in undesirable local minima. After introducing all energy terms, we discuss how this problem can be divided into smaller problems that can be solved in order, thereby allowing to find a good minimum in practice.

Prior model for β

A prior model is used to ensure that the body shape stays within the learned shape space that represents plausible human shapes. The identity shape space is learned using PCA, and has zero mean and standard deviation σ i along the i-th principal component. Similarly to previous work [START_REF] Balan | The naked truth: Estimating body shape under clothing[END_REF], we do not penalize values of β that stay within 3σ i of the mean to avoid introducing a bias towards the mean shape. However, rather than penalizing a larger distance from the mean, we constrain the solution to lie inside the hyperbox ±3σ i using a constrained optimization framework. This constraint can be handled by standard constrained optimizers since the hyperbox is axis-aligned, and using this hard constraint removes the need to appropriately weigh a prior energy w.r.t. other energy terms.

Landmark energy

The landmark energy helps to guide the solution towards the desired local minimum with the help of distinctive anatomical landmarks. This energy is especially important during the early stages of the optimization as it allows to find a good initialization for the identity and posture parameters. In the following, we consider the use of N lnd landmarks and assume without loss of generality that the vertices corresponding to landmarks are the first N lnd vertices of s. The landmark term is defined as

E lnd (F , β, Θ) = N lnd i=1 s i (β, Θ) -l i (F ) 2 , (6)
where l i (F ) denotes the i-th landmark of frame F , s i (β, Θ) denotes the vertex corresponding to the i-th landmark of s (β, Θ), and • denotes the 2 norm. The landmarks l i (F ) are computed automatically with the help of the state of the art Stitched Puppet [START_REF] Zuffi | The Stitched Puppet: A graphical model of 3d human shape and pose[END_REF], which allows to robustly fit a human body model to a single scan using a particle-based optimization. Specifically, we once manually select a set of vertex indices to be used as landmarks on the Stitched Puppet model, which is then fixed for all experiments. To fit the Stitched Puppet to a single frame, randomly distributed particles are used to avoid getting stuck in undesirable local minima. We fit the Stitched Puppet model to frame F , and report the 3D positions of the pre-selected indices after fitting as landmarks l i (F ). While the Stitched Puppet aims to fit the body shape and posture of F , only the coordinates l i (F ) are used by our framework. Note that our method does not require accurate l i (F ), since l i (F ) are only used to initialize the optimization.

Using many particles on each frame of a motion sequence is inefficient. Furthermore, since the Stitched Puppet is trained on a database of minimally dressed subjects, using many particles to fit to a frame in wide clothing may lead to overfitting problems. This is illustrated in Fig. 2. To remedy this, we choose to use a relatively small number of particles which is set to 30. Starting at the second frame, we initialize the particle optimization to the result of the previous frame to guide the optimization towards the desired optimum.

Data energy

The data energy pulls the S-SCAPE model towards the observation F using a nearest neighbor term. This energy, which unlike the landmark energy considers all vertices of s, is crucial to fit the identity and posture of s to the input F as

E data (F , β, Θ) = N v i=1 δ N N s i (β, Θ) -N N (s i (β, Θ) , F ) 2 , (7)
where N v denotes the number of vertices of s and N N (s i (β, Θ) , F ) denotes the nearest neighbour of vertex s i (β, Θ) on F . To remove the influence of outliers and reduce the possibility of nearest neighbour mismatching, we use a binary weight δ N N that is set to one if the distance between s i and its nearest neighbor on F is below 200mm and the angle between their outer normal vectors is below 60 • , and to zero otherwise.

Clothing energy

The clothing energy is designed to encourage the predicted body shape s to be located entirely inside the observation F . This energy is particularly important when considering motion sequences acquired with loose clothing. In such cases, merely using E lnd and E data leads to results that overestimate the circumferences of the body shape because β is estimated to fit to F rather than to fit inside of F , see Fig. 3. To remedy this, we define the clothing energy as

E cloth (F ,β ,Θ )= N v i=1 δ out δ N N s i (β ,Θ )-N N (s i (β ,Θ ),F ) 2 +ω r β -β 0 2 , ( 8
)
where δ out is used to identify vertices of s located outside of F . This is achieved by setting δ out to one if the angle between the outer normal of N N (s i (β, Θ) , F ) and the vector s i (β, Θ) -N N (s i (β, Θ) , F ) is below 90 • , and to zero otherwise. Furthermore, ω r is a weight used for the regularization term, and β 0 is an initialization of the identity parameters used to constrain β.

When observing a human body dressed in loose clothing in motion, different frames can provide valuable cues about the true body shape. The energy E cloth is designed to exploit motion cues when optimizing E cloth w.r.t. all available observations F i . This allows to account for clothing using a simple optimization without the need to find skin and non-skin regions as in previous work [START_REF] Balan | The naked truth: Estimating body shape under clothing[END_REF][START_REF] Stoll | Video-based reconstruction of animatable human characters[END_REF][START_REF] Neophytou | A layered model of human body and garment deformation[END_REF]. The regularization ββ 0 2 used in Eq. 8 is required to avoid excessive thinning of limbs due to small misalignments in posture.

Fig. 3 shows the influence of E cloth on the result of a walking sequence in layered clothing. The left side shows overlays of the input and the result for ω cloth = 0 and ω cloth = 1. Note that while circumferences are overestimated when ω cloth = 0, a body shape located inside the input frame is found for ω cloth = 1. The comparison to the ground truth body shape computed as discussed in Sec. 6 is visualized in the middle and the right of Fig. 3, and shows that E cloth leads to a significant improvement of the accuracy of β.

Optimization schedule

Minimizing E F 1:N f , β, Θ 1:N f defined in Eq. 5 over all N f frames w.r.t. β and Θ i jointly is not feasible when considering motion sequences containing hundreds of frames as this is a high-dimensional optimization problem. To solve this problem without getting stuck in undesirable local minima, we optimize three smaller problems in order. Initial identity estimation. We start by computing an initial estimate β 0 based on the first N k frames of the sequence by optimizing E F 1:N k , β, Θ 1:N k w.r.t. β and Θ i . For increased efficiency, we start by computing optimal β i and Θ i for each frame using Eq. 4 by alternating the optimization of Θ i for fixed β i with the optimization of β i for fixed Θ i . This is repeated for N it iterations. Temporal consistency is achieved by initializing Θ i+1 as Θ i and β i+1 as β i starting at the second frame. As it suffices for the identity parameters to roughly estimate the true body shape at this stage, we set ω cloth = 0. In the first iterations, E lnd is essential to guide the fitting towards the correct local optimum, while in later iterations E data gains in importance. We therefore set ω data = 1 -ω lnd and initialize ω lnd to one. We linearly reduce ω lnd to zero in the last two iterations. We then initialize the posture parameters to the computed Θ i , and the identity parameters to the mean of the computed β i and iteratively

minimize E F 1:N k , β, Θ 1:N k w.r.t. Θ 1:
N k and β. This leads to stable estimates for Θ 1:N k and an initial estimate of the identity parameter, which we denote by β 0 in the following.

Posture estimation. During the next stage of our framework, we compute the posture parameters Θ N k +1:N f for all remaining frames by sequentially minimizing Eq. 4 w.r.t. Θ i . As before, Θ i+1 is initialized to the result of Θ i . As the identity parameters are not accurate at this stage, we set ω cloth = 0. For each frame, the energy is optimized N it times while reducing the influence of ω lnd in each iteration, using the same weight schedule as before. This results in posture parameters Θ i for each frame.

Identity refinement. In a final step, we refine the identity parameters to be located inside all observed frames F 1:N f . To this end, we initialize the identity parameters to β 0 , fix all posture parameters to the computed Θ i , and minimize

E F 1:N f , β, Θ 1:N f w.r.t. β.
As the landmarks and observations are already fitted adequately, we set ω lnd = ω data = 0 at this stage of the optimization.

Implementation details

The S-SCAPE model used in this work consists of N v = 6449 vertices, and uses d id = 100 parameters to control identity and d pose = 30 parameters to control posture by rotating the N b = 15 bones. The bones, posture parameters, and rigging weights are set as in the published model [START_REF] Pishchulin | Building statistical shape spaces for 3d human modeling[END_REF].

For the Stitched Puppet, we use 60 particles for the first frame, and 30 particles for subsequent frames. We use a total of N lnd = 14 landmarks that have been shown sufficient for the initialization of posture fitting [START_REF] Wuhrer | Estimation of human body shape and posture under clothing[END_REF], and are located at forehead, shoulders, elbows, wrists, knees, toes, heels, and abdomen. Fig. 1 shows the chosen landmarks on the Stitched Puppet model. During the optimization, we set N it = 6 and N k = 25 . The optimization w.r.t. β uses analytic gradients, and we use Matlab L-BFGS-B to optimize the energy. The setting of the regularization weight ω r depends on the clothing style. The looser the clothing, the smaller ω r , as this allows for more corrections of the parameters. In our experiments, we use ω r = 1 for all the sequences with layered and wide clothing in our dataset.

Evaluation

Dataset

This section introduces the new dataset we acquired to allow quantitative evaluation of human body shape estimation from dynamic data. The dataset consists of synchronized acquisitions of dense unstructured geometric motion data and sparse motion capture (MoCap) data of 6 subjects (3 female and 3 male) captured in 3 different motions and 3 clothing styles each. The geometric motion data are sequences of meshes obtained by applying a visual hull reconstruction to a 68-color-camera (4M pixels) system at 30FPS. The basic motions that were captured are walk, rotating the body, and pulling the knees up. The captured clothing styles are very tight, layered (long-sleeved layered clothing on upper body), and wide (wide pants for men and dress for women). The body shapes of 6 subjects vary significantly. Fig. 4 shows some frames of the database.

To evaluate algorithms using this dataset, we can compare the body shapes estimated under loose clothing with the tight clothing baseline. The comparison is done per vertex on the two body shapes under the same normalized posture. Cumulative plots are used to show the results.

Evaluation of posture and shape fitting

We applied our method to all sequences in the database. For one sequence of a female subject captured while rotating the body in wide clothing, Stitched Puppet fails to find the correct posture, which leads to a failure case of our method (see Fig. 2). We exclude this sequence from the evaluation. To evaluate the accuracy of the posture parameters Θ, we compare the 3D locations of a sparse set of landmarks captured using a MoCap system with the corresponding model vertices of our estimate. This evaluation is performed in very tight clothing, as no accurate MoCap markers are available for the remaining clothing styles. Fig. 5 summarizes the per-marker errors over the walking sequences of all subjects. The results show that most of the estimated landmarks are within 35mm of the ground truth and that our method does not suffer from drift for long sequences. As the markers on the Stitched Puppet and the MoCap markers were placed by non-experts, the landmark placement is not fully repeatable, and errors of up to 35mm are considered fairly accurate. To evaluate the accuracy of the identity parameters β, we use for each subject the walking sequence captured in very tight clothing to establish a ground truth identity β 0 by applying our shape estimation method. Applying our method to sequences in looser clothing styles of the same subject leads identity parameters β, whose accuracy can be evaluated by comparing the 3D geometry of s (β 0 , Θ 0 ) and s (β, Θ 0 ) for a standard posture Θ 0 . Fig. 6 summarizes the per-vertex errors over all motion sequences captured in layered and wide clothing, respectively. The left side shows the cumulative errors, and the right side shows the color-coded mean per-vertex error. The color coding is visualized on the mean identity of the training data. The result shows that our method is robust to loose clothing with more than 50% of all the vertices having less than 10mm error for both layered and wide clothing. The right side shows that as expected, larger errors occur in areas where the shape variability across different identities is high. Fig. 7 shows some qualitative results for all three types of motions and two clothing styles. Note that accurate body shape estimates are obtained for all frames. Consider the frame that shows a female subject performing a rotating motion in layered clothing. Computing a posture or shape estimate based on this frame is extremely challenging as the geometry of the layered cloth locally resembles the geometry of an arm, and as large portions of the body shape are occluded. Our method successfully leverages temporal consistency and motion cues to find reliable posture and body shape estimates.

Comparative evaluation

As we do not have results on motion sequences with ground truth for existing methods, this section presents visual comparisons, shown in Fig. 8. We compare to Wuhrer et al. [START_REF] Wuhrer | Estimation of human body shape and posture under clothing[END_REF] on the dancer sequence [START_REF] De Aguiar | Performance capture from sparse multi-view video[END_REF] presented in their work. Note that unlike the results of Wuhrer et al., our shape estimate does not suffer from unrealistic bending at the legs even in the presence of wide clothing. Furthermore, we compare to Neophytou and Hilton [START_REF] Neophytou | A layered model of human body and garment deformation[END_REF] on the swing sequence [START_REF] Vlasic | Articulated mesh animation from multi-view silhouettes[END_REF] presented in their work. Note that we obtain results of similar visual quality without the need for manual initializations and pre-aligned motion sequences. In summary, we present the first fully automatic method for body shape and motion estimation, and show that this method achieves state of the art results.

Comparison to

Comparison to Wuhrer et al. [START_REF] Wuhrer | Estimation of human body shape and posture under clothing[END_REF] Neophytou and Hilton [START_REF] Neophytou | A layered model of human body and garment deformation[END_REF] Fig. 8. Per comparison from left to right: input, result of prior works, our result.

Conclusion

We presented an approach to automatically estimate the human body shape under motion based on a 3D input sequence showing a dressed person in possibly loose clothing. The accuracy of our method was evaluated on a newly developed benchmark1 containing 6 different subjects performing 3 motions in 3 different styles each. We have shown that, although being fully automatic, our posture and shape estimation achieves state of the art performance. In the future, the body shape and motion estimated by our algorithm have the potential to aid in a variety of tasks including virtual change rooms and security applications.

Fig. 1 .

 1 Fig. 1. Overview of the proposed pipeline. From left to right: input frame, result of Stitched Puppet [1] with annotated landmarks, result after estimation of initial identity and posture, final result, and overlay of input and final result.

Fig. 2 .

 2 Fig. 2. Left: overfitting problem of Stitched Puppet in the presence of clothing. Input frame, Stitched Puppet result with 160 particles, and Stitched Puppet result with 30 particles are shown in order. Right: the failure case from our database caused by mismatching of Stitched Puppet.

Fig. 3 .

 3 Fig. 3. Influence of E cloth on walking sequence. Left: input data overlayed with result with ω cloth = 0 (left) and ω cloth = 1 (right). Middle: cumulative per-vertex error of estimated body shape with ω cloth = 0 and ω cloth = 1. Right: color-coded per-vertex error with ω cloth = 0 (left) and ω cloth = 1 (right).

Fig. 4 .

 4 Fig. 4. Six representative examples of frames of our motion database. From left to right, a female and male subject is shown for tight, layered, and wide clothing each.

Fig. 5 .

 5 Fig. 5. Accuracy of posture estimation over the walking sequences in tight clothing. Left: cumulative landmark errors. Right: average landmark error for each sequence.

Fig. 6 .

 6 Fig. 6. Summary of shape accuracy computed over the frames of all motion sequences of all subjects captured in layered and wide clothing. Left: cumulative plots showing the per-vertex error. Right: mean per-vertex error color-coded from blue to red.

Fig. 7 .

 7 Fig. 7. Overlay of input data and our result.

The benchmark can be downloaded at http://dressedhuman.gforge.inria.fr/.

Acknowledgments

Funded by France National Research grant ANR-14-CE24-0030 ACHMOV. We thank Yannick Marion for help with code to efficiently fit an S-SCAPE model to a single frame, Leonid Pishchulin for helpful discussions, Alexandros Neophytou and Adrian Hilton for providing comparison data, and Mickaël Heudre, Julien Pansiot and volunteer subjects for help acquiring the database. The database was acquired using the Kinovis platform: http://kinovis.inrialpes.fr/.