Siamak Niroomandi

David Gonzalez

Icíar Alfaro

Felipe Bordeu Weldt

Adrien Leygue

Elías Cueto

Francisco Chinesta

D González

F Bordeu

Real Time Simulation of Biological Soft Tissues: A PGD Approach

Keywords: real time, surgery, model order reduction, proper generalized decomposition

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

INTRODUCTION

Real-time simulation is one of the most challenging scenarios for simulation-based engineering sciences. The term real time strongly depends on the particular pursued application, but surgery simulation is among the most restrictive ones. Haptic surgery simulators compute the response of biological soft tissues and give it back to the peripherals at, at least, 25 Hz of feedback rate if visual realism is needed and, notably, 500 Hz-1 kHz if haptic (force) response is desired [START_REF] Delingette | Hepatic surgery simulation[END_REF][START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF]. But biological soft tissues are known to be highly nonlinear [START_REF] Fung | Biomechanics: Mechanical Propeties of Living Tissues[END_REF][START_REF] Alastrue | Biomechanical modeling of refractive corneal surgery[END_REF][START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF], very frequently modeled in a hyperelastic framework. It is well-known, in addition, that at least nonlinear strain measures should be incorporated into the simulation, otherwise performing arbitrarily bad in terms of visual perception, spurious gain in volume, and others [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF].

This feedback rate for nonlinear problems constitutes indeed a great challenge for nowaday simulation techniques based upon FEMs. This is perhaps the ultimate reason for the lack, up to our knowledge, of surgery simulators of the second generation [START_REF] Delingette | Hepatic surgery simulation[END_REF], that is, those that incorporate stateof-the-art constitutive modeling of soft tissues to the simulation. A few references can be cited that incorporate nonlinear tissue behavior (mainly Kirchhoff-Saint Venant hyperelasticity) (see [START_REF] Taylor | On modelling of anisotropic viscoelasticity for soft tissue simulation: numerical solution and GPU execution[END_REF][START_REF] Taylor | High-speed nonlinear finite element analysis for surgical simulation using graphics processing units[END_REF][START_REF] Lim | Real time simulation of nonlinear tissue response in virtual surgery using the point collocation-based method of finite spheres[END_REF]). Many of them are indeed based upon explicit finite elements, possibly implemented in graphics processing units. This is due to the very astringent conditions dictated by the haptic feedback rates that prevent the use, in nowaday computers, of standard Newton-like methods for the solution of nonlinear systems of equations.

Recently, model order reduction has been seen as a powerful means to achieve real-time performance in the simulation of nonlinear solids. For instance, Barbič and James [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF] proposed recently a method for model order reduction of Kirchhoff-Saint Venant solids at haptic feedback rates. For more sophisticated constitutive laws, recent work on POD methods showed that visual feedback performance can be easily obtained, but this is not the case very often if we need haptic feedback [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF][START_REF] Niroomandi | Real-time simulation of surgery by reduced-order modeling and X-FEM techniques[END_REF][START_REF] Taylor | A reduced order explicit dynamic finite element algorithm for surgical simulation[END_REF]. In this case, linearization schemes for the resulting nonlinear systems of equations is mandatory, although not straightforward to accomplish. Very often, Newton strategies are out of reach, due to the inherent lack of time employed during stiffness matrix updates. In [START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF], a different strategy was developed that employed a combination of POD techniques with asymptotic numerical methods [START_REF] Cochelin | Asymptotic-numerical methods and Padé approximants for non-linear elastic structures[END_REF][START_REF] Abichou | Asymptotic numerical method for problems coupling several nonlinearities[END_REF][START_REF] Cao | An improved iterative method for large strain viscoplastic problems[END_REF][START_REF] Cochelin | The asymptotic numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF]. In that case, no tangent stiffness matrix updates are necessary, and only a set of linear problems, all with the same stiffness matrix, must be solved.

Proper generalized decomposition (PGD) methods, on the contrary, arose recently as a means to overcome the so-called curse of dimensionality associated to problems defined in high-dimensional spaces [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: transient simulation using space-time separated representations[END_REF][START_REF] Chinesta | Recent advances in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Ladeveze | The LATIN multiscale computational method and the proper generalized decomposition[END_REF]. Although the origins of the technique can be traced back to the large time increment method and its associated radial loading [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF], PGD is nowadays seen as a powerful method of model order reduction that generalizes POD (and hence its name). Its field of application has gone far beyond the initial objectives, and it is now being applied in a variety of fields: simulation-based control of processes [START_REF] Ghnatios | Methodological approach to efficient modeling and optimization of thermal processes taking place in a die: application to pultrusion[END_REF][START_REF] Ghnatios | Proper generalized decomposition based dynamic data-driven control of thermal processes[END_REF], efficient simulation of plates and shells [START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF], and simulation of gene regulatory networks [START_REF] Ammar | Reduction of the chemical master equation for gene regulatory networks using proper generalized decompositions[END_REF], to cite a few.

The main idea for the extension of a technique initially developed to deal with high-dimensional problems to a more general setting lies precisely in its ability to treat standard models as if they were multidimensional. The key idea was initially put forth in [START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the proper generalized decomposition[END_REF]. In it, parametric equations were cast into a multidimensional setting, thus taking advantage of the PGD solution structure. Going one step forward, PGD can advantageously be employed into real-time simulation frameworks by simply considering not only all parameters but also all possible boundary conditions (including initial boundary conditions, see [START_REF] Gonzalez | Proper generalized decomposition based dynamic data driven inverse identification[END_REF]) as new dimensions of the problem.

The resulting PGD solution to the problem is expressed as a finite sum of separable functions that provides actually a meta-model for the problem, for which no prior computer experiment (also known as snapshots in the model order reduction community) is necessary. This meta-model can then be successfully applied in real time to obtain the response of the system at kilohertz rates, as will be demonstrated in subsequent sections of this paper. This approach allows even to solve models on handheld devices such as tablets and smartphones, thus opening the range of possible applications of the technique.

In this paper, we present a novel technique for the simulation of biological soft tissues under hyperelasticity assumptions at haptic feedback rates. It is based on the use of the aforementioned PGD approach and an explicit linearization of the weak form generated by nonlinear strain measures. The use of PGD therefore allows not only for a complete generalization of previous works in the field (see [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF][START_REF] Niroomandi | Real-time simulation of surgery by reduced-order modeling and X-FEM techniques[END_REF]) but also to a completely new formulation of the problem. For the problem at hand, POD techniques compute the solution (the so-called snapshots) of different complete models for different contact positions between a surgical tool-scalpel-and an organ. By performing a statistical analysis on top of these results, POD techniques extract the so-called modes, that is, those displacement fields that best represent the solution of the complete problems. These modes, taken as Ritz-like, global basis, are then used to approximate problems different to the original ones, that is, for contact positions not initially considered. Therefore, the design of an appropriate simulation, or computer experiment, campaign is the key aspect of the method. Also, how to efficiently and accurately interpolate among reduced basis is another crucial aspect not fully resolved.

On the contrary, PGD methods consider a parametric problem. In this case, the parameter is the contact position. By formulating the parametric problem as a high-dimensional one, PGD methods allow for the efficient solution and ulterior storage of the model in the form of a sum of separable functions. Therefore, rather than creating reduced models for particular positions of the tool, PGD computes a general solution for any position of the tools, so that no subsequent interpolation of reduced models is necessary, nor the computation of snapshots. Thus, the proposed method is based upon an offline phase in which this general solution is computed and an online one in which the solution is only evaluated at impressive feedback rates, here on the order of kilohertz.

The paper is organized as follows. In Section 2, we introduce the basics of PGD applied to the problem of a hyperelastic solid under moving punctual loads, which is the most frequent case in surgery simulation. In Section 3, a very simple linearization of the nonlinear problem is introduced that allows for a simple yet effective computation of the PGD approach to the problem. Although this simple linearization is by no means the only possible one, its performance is analyzed in Section 4, through a series of benchmark problems. It is shown how the PGD approach to the problem of real-time simulation of soft tissue deformation opens new insights on how the problem can be attacked.

A PROPER GENERALIZED DECOMPOSITION APPROACH TO VIRTUAL SURGERY

As already mentioned in Section 1, the key issue in the use of PGD approaches for real-time simulation, and the one that makes it completely different in spirit from POD, lies in the formulation of the original problem as a parametric one. This parametric problem is then reformulated as a highdimensional problem by considering each parameter as a new coordinate in the state space. The PGD method then looks for an effective solution in the form of a finite sum of separable functions, so as to be able to avoid the curse of dimensionality associated to high-dimensional problems and mesh-based discretization techniques.

In this framework, the problem of determining the response of an organ to the load transmitted by the contact with a surgical tool could be formulated to determine the displacement at any point of the model, u.x, y, ´/, for any load position s and for any force vector orientation and module, t, thus rendering a problem defined in the physical space .R 3 /, plus a six-dimensional state space .R 6 /.

For the sake of simplicity in the following development, and without loss of generality, we assume a load vector t with unit module and oriented in the vertical direction. This renders a problem defined in R 6 .u D u.x, s//, with all the characteristics of the aforementioned one.

Let us consider the weak form of the equilibrium equations (balance of linear momentum). Again, for the sake of simplicity, we omit inertia terms. The interested reader could consult PGD-ODE [START_REF] Gonzalez | Proper generalized decomposition based dynamic data driven inverse identification[END_REF] for the treatment of ODEs in the framework of PGD. Under these assumptions, the weak form of the problem, extended to the whole geometry of the organ, and the portion of its boundary that is accessible to the surgeon, N t , ‡ consists in finding the displacement u 2 H 1 such that for all u 2 H 1 0 :

Z N Z r s u W d d N D Z N Z t2 u tdd N (1)
where D u [ t represents the boundary of the organ, divided into essential and natural regions, and where t D t1 [ t2 , that is, regions of homogeneous and nonhomogeneous, respectively, natural boundary conditions. Here, t D e k ı.x s/, where ı represents the Dirac delta function and e k the unit vector along the ´-coordinate axis (we consider here, for the ease of exposition, a unit load directed toward the negative ´axis of reference).

Once regularized, the Dirac delta term is approximated by a truncated series of separable functions in the spirit of the PGD method, that is,

t j m X iD1 f i j .x/g i j .s/ (2)
where m represents the order of truncation and f i j , g i j represent the j th component of vectorial functions in space and boundary positions, respectively. ‡ Typically, of all of the natural regions of the boundary, t -which in this case coincides with virtually all the boundary of organs, D @ , because they are not fixed or clamped but in contact to other organs-only a portion is accessible to the surgeon. In minimally invasive surgery, surgeons operate through a small incision on the skin of the patient, having access only to a limited portion of the boundary, here termed N .

The PGD approach to the problem is characterized by the construction, in an iterative way, of an approximation to the solution in the form of a finite sum of separable functions. Assume that we have converged to a solution, at iteration n of this procedure,

u n j .x, s/ D n X kD1 X k j .x/ Y k j .s/, ( 3 )
where the term u j refers to the j th component of the displacement vector, j D 1, 2, 3, and functions X k and Y k represent the separated functions used to approximate the unknown field, obtained in previous iterations of the PGD algorithm.

If we look for an improvement of this approximation, the (n C 1)th term will look like

u nC1 j .x, s/ D u n j .x, s/ C R j .x/ S j .s/, ( 4 )
where R.x/ and S .s/ are the sought functions that improve the approximation. In this framework, the admissible variation of the displacement will be given by

u j .x, s/ D R j .x/ S j .s/ C R j .x/ S j .s/. ( 5 )
At this point, several options are at hand so as to determine the new pair of functions R and S . The most frequently used, due to both its ease of implementation and good convergence properties, in general, is a fixed-point algorithm in which functions R and S are sought iteratively. We describe briefly the implementation of this algorithm.

Computation of S.s/ assuming R.x/ is known

In this case, following standard assumptions of variational calculus, we have

u j .x, s/ D R j .x/ S j .s/, ( 6 )
or, equivalently, u .x, s/ D R ı S , where the symbol 'ı' denotes the so-called entry-wise Hadamard or Schur multiplication for vectors. Once substituted into Equation (1), it gives

Z N Z r s .R ı S / W C W r s n X kD1 X k ı Y k C R ı S ! d d N D Z N Z t2 .R ı S / m X kD1 f k ı g k ! dd N , ( 7 )
or equivalently (we omit obvious functional dependencies)

Z N Z r s .R ı S / W C W r s .R ı S /d d N D Z N Z t2 .R ı S / m X kD1 f k ı g k ! dd N Z N Z r s .R ı S / R n d d N , (8)
where

R n represents R n D C W r s u n . ( 9
)
Because the symmetric gradient operates on spatial variables only, we have

Z N Z .r s R ı S / W C W .r s R ı S /d d N D Z N Z t2 .R ı S / m X kD1 f k ı g k ! dd N Z N Z .r s R ı S / R n d d N (10)
where all the terms depending on x are known, and hence, we can compute all integrals over and t2 (support of the regularization of the initially punctual load) to derive an equation to compute S .s/.

Computation of R.x/ assuming S.s/ is known

Equivalently, in this case, we have

u j .x, s/ D R j .x/ S j .s/, (11)
which, once substituted into Equation ( 1), gives

Z N Z r s .R ı S / W C W r s n X kD1 X k ı Y k C R ı S ! d d N D Z N Z t2 .R ı S / m X kD1 f k ı g k ! dd N . ( 12
)
In this case, all the terms depending on s (load position) can be integrated over N , leading to a generalized elastic problem to compute function R.x/.

This simple algorithm renders, in general, excellent convergence properties (see [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] and references therein).

ONE POSSIBLE EXPLICIT LINEARIZATION OF THE FORMULATION

The formulation introduced in Section 2 assumes implicitly small strains. But this assumption has been found to be very insufficient for virtual surgery. Strains are large very often, and when solved in a small strain setting, organs appear to suffer an unphysical gain in volume that renders the simulations clearly nonrealistic [START_REF] Delingette | Hepatic surgery simulation[END_REF]. Although real-time simulation of surgery does not look nowadays for accuracy levels similar to those common in usual engineering practice (quoting Cotin and Bro-Nielsen, [START_REF] Bro-Nielsen | Real-time volumetric deformable models for surgery simulation using finite elements and condensation[END_REF], '... the model may be physically correct if it looks right'), the inclusion of nonlinear strain measures seems to be crucial.

Soft tissues are frequently formulated under hyperelasticity assumptions [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF]. Again, for the sake of simplicity, we refer ourselves to a Kirchhoff-Saint Venant constitutive framework. Despite being very limited (and even unstable under compression due to the lack of polyconvexity of the strain energy functional), Kirchhoff-Saint Venant constitutive equations are widely used at this moment for real-time simulation of soft tissues (see [START_REF] Taylor | High-speed nonlinear finite element analysis for surgical simulation using graphics processing units[END_REF][START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Taylor | A reduced order finite element algorithm for surgical simulation[END_REF] among others).

The Kirchhoff-Saint Venant model is characterized by the energy density functional given by

‰ D 2 .tr.E // 2 C E W E (13)
where and are Lame's constants. The Green-Lagrange strain tensor, E , has the form

E D 1 2 F T F 1 Á D r s u C 1 2 r u r u T (14)
where F D r u C I is the gradient of deformation tensor. The second Piola-Kirchhoff stress tensor can be obtained by

S D @‰.E / @E D C W E ( 15
)
in which C is the fourth-order constitutive (elastic) tensor. Very few works have been written about PGD approximations for nonlinear solid mechanics problems, other than the Ladeveze's works in complex thermomechanical nonlinear models [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF]. Here, we focus on the nonlinear and parametric case within a fully separated representation. In this case, we restrict ourselves to quasi-static problems, for the sake of simplicity, and therefore introduce a pseudo-time t 2 OE0, 1 to perform the linearization.

Consistent linearizations of the resulting set of equations in the framework of PGD approximations are far from being trivial, so here we keep the formulation as simple as possible by performing a simple explicit linearization of the weak form of the problem.

Thus, load is applied along a series of time increments t, provoking increments in the displacement u.x, s/. At each time increment, a PGD fixed-point alternating directions algorithm similar to those introduced in Section 2 is employed. So, if we introduce the nonlinear strain measure given by Equation [START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF], into this incremental framework, we have (we omit obvious functional dependencies for clarity)

E t Ct D r s .u t C u/ C 1 2 r .u t C u/ r T .u t C u/ . (16)
Similarly, admissible variation of strain reads

E D r s .u / C 1 2 .r .u // r T .u t C u/ C 1 2 r .u t C u/ r T .u / D r s .u / C r .u / r T .u t C u/. ( 17
)
Once substituted into the weak form of the equilibrium equation, Equations ( 16) and ( 17) provide, for the left hand side term of Equation (1)-strain energy term-,

Z N Z .t / E W C W E d d N D Z N Z .t / r s .u / C r .u / r T .u t C u/ W C W Â r s u t C u C 1 2 r .u t C u/ r T .u t C u/ Ã d d N . ( 18
)
The simplest linearization of Equation ( 18) consists of keeping in the formulation only constant terms and those linear in u. We thus arrive at a weak form composed of 10 terms:

Z N Z .t / E W C W E d d N D Z N Z .t / r s .u / W C W r s u t d d N " ƒ' … T 1 C Z N Z .t / r s .u / W C W r s .u/d d N " ƒ' … T 2 C Z N Z .t / r s .u / W C W 1 2 r u t r T u t d d N " ƒ' … T 3 C Z N Z .t / r s .u / W C W r u t r T .u/d d N " ƒ' … T 4 C Z N Z .t / r .u / r T u t W C W r s u t d d N " ƒ' … T 5 C Z N Z .t / r .u / r T u t W C W r s .u/d d N " ƒ' … T 6 C Z N Z .t / r .u / r T u t W C W 1 2 r u t r T u t d d N " ƒ' … T 7 C Z N Z .t / r .u / r T u t W C W r u t r T .u/d d N " ƒ' … T 8 C Z N Z .t / r .u / r T .u/ W C W r s u t d d N " ƒ' … T 9 C Z N Z .t / r .u / r T .u/ W C W 1 2 r u t r T u t d d N " ƒ' … T 10 . ( 19
)
This renders a very simple scheme that has revealed, however, for judicious choice of the time step t, reasonable convergence properties, as will be demonstrated in Section 4.

Remark 1

The original work of P. Ladeveze on the large time increment method [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF] combined a space-time separated representation and thus produces a nonincremental solution of the problem. Generalized to this case, the displacement would be sought in the form u D u.x, s, t/. We have preferred, for simplicity of exposition, to keep the formulation as simple as possible, but the explicit linearization proposed in Equation ( 19) is by no means the only possible one.

Remark 2

Another possible choice for the aforementioned linearization is the standard forward Euler scheme. It has been noted, however, that instabilities in the results appear at zones subjected to compression, a typical characteristic of Kirchhoff-Saint Venant models [START_REF] Bonet | Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations[END_REF]. These are analyzed in Section 4. However, no spurious deformation modes have been observed by employing the simple explicit algorithm stated in Equation ( 19).

NUMERICAL RESULTS

Validation: nonlinear rod

To validate the explicit PGD approach introduced in Section 2 earlier, we considered the simple case of a Kirchhoff-Saint Venant beam subjected to a pure traction force F . In this simple case, the model can be simplified to a one-dimensional one, and an analytical solution for the axial displacement u at the bar tip can be obtained as

u D  ˇ 1 C 1 3ˇÃ , ( 20
)
where

ˇD K C r K 2 1 27 ! 1 3 , (21)
and, in turn, K D F=EA, with E as the Young's modulus of the material and A the area of the cross section at the undeformed configuration.

A model was thus constructed by considering a bar of length L D 400 mm, A D 40 40 mm 2 , E D 1.0 MPa, and F D 320 N. Load F , however, is considered to be applied at any point along the bar axis. We therefore compute a two-dimensional solution u D u.x, s/, where x represents the position along the bar axis and s the point of application of the load.

This simple example served to know the importance of the chosen time step in the overall convergence properties of the proposed method. As can be noticed from Figure 1, the error for a time step of 10 3 is O.10 4 /. Note that one single PGD approach is used, which is enriched at each time step, not a different PGD approximation within each time step.

It is important to remark that the real-time strategy here introduced is based upon the computation (only once and off-line) of a general, high-dimensional solution of the problem once for life. This general solution is then evaluated at real-time feedback rates but not recomputed. That is why the time taken in the offline computation of this general solution is not so important, because it will be performed only once.

Kirchhoff-Saint Venant beam bending

Further validation of the proposed strategy is obtained if we consider the problem of a beam bending under a transverse load (assumed vertical, for simplicity) applied at any point of its boundary. The problem has no analytical solution, up to our knowledge, and therefore, the general, multidimensional solution has been compared with a reference one obtained by standard finite element models (one for each considered load position) with consistent linearization and a Newton-Raphson iterative scheme.

The mesh is composed of tetrahedral elements, with only 9 9 nodes in the 40 40 mm 2 cross section and 21 nodes in the longitudinal direction, 400 mm long. Material parameters were Young's module E D 1.0 MPa and Poisson's coefficient D 0.25. With such a poor discretization, and employing tetrahedral elements, it is expected that a high error with respect to the exact solution will be obtained. However, it is not the purpose of this paper to obtain an accurate enough solution of the problem (which of course could be obtained by employing a more refined mesh and a finer time stepping).

The finite element model is shown in Figure 2, and the load is assumed to be applied at any of the points of the upper surface of the beam.

The obtained displacement for a particular location of the load (beam tip in this case, for comparison purposes) is depicted in Figure 3. Noteworthy, the simple explicit linearization algorithm here proposed does not imply a nonphysical gain of volume in the deformed model, which is the case for purely linear elastic models, nevertheless frequently employed in real-time simulation of surgery [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF].

Standard forward Euler algorithms for the linearization of the weak form of the problem rendered, in our experiments, spurious deformation modes. Although they are intrinsic of the Kirchhoff-Saint Venant constitutive model, it continues to be popular among the virtual surgery community because it constitutes the fastest way to avoid spurious gain in the volume typical of linear elastic models, when large deformations are being considered [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF]. An example of the result given by a forward Euler scheme for this same problem is shown in Figure 4. Although no unphysical gain in volume-typical of linear elastic approaches to the problem-is seen, the obtained displacement at beam tip is much higher than the reference one, obtained by FEMs. This is due to some well-known instabilities of the Kirchhoff-Saint Venant model under compression.

In general, the number n of functions employed in the approximation depends on the desired level of accuracy. Increasing the number of separated functions in the approximation leads, of course, to higher computational costs in the offline part of the method. But we highlight that this computation is carried out only once for life and stored in memory. The online part of the simulation is virtually not affected, because only some vector multiplications should be performed in addition, which do not alter the overall efficiency of the proposed method.

Palpation of the liver

The liver is the biggest gland in the human body, after the skin. Liver geometry has been obtained from the SOFA project [START_REF] Allard | SOFA an open source framework for medical simulation[END_REF] and post-processed to obtain a mesh composed of 8559 nodes and 10,519 tetrahedra (Figure 5). The liver is connected to the diaphragm by the coronary ligament so it seems reasonable to assume it to be constrained at the posterior face by the rest of the organs, whereas the anterior face is accessible to the surgeon. The inferior vena cava travels along the posterior surface, and the liver is frequently assumed clamped at that location. Although the assumed boundary conditions are not strictly correct from a physiological point of view, our main interest is to show that the model can be solved under real-time constraints with reasonable accuracy. Although the literature on the mechanical properties of the liver is not very detailed, we have assumed a Young's modulus of 160 kPa, and a Poisson coefficient of 0.48, thus nearly incompressible [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF].

The N surface, where the load can be located, has been defined as the whole boundary of the domain, even if in this case, only the frontal part of the organ is usually accessible to the surgeon. This region includes 2009 of the 8559 nodes of the model.

Model's solution was composed by a total of n D 167 functional pairs X k j .x/ Y k j .s/ (see Equation ( 3)). The third component (thus j D 3) of the first six modes X k 3 .x/ is depicted in Figure 6. The same is carried out in Figure 7 for functions Y , although in this case, they are defined only on the boundary of the domain, that is, N D @ . It is noteworthy that both X and Y sets of functions present a structure similar to that generated by POD methods, despite the fact that they are not, in general, optimal. Note how the frequency content of each pair of functions increases as we increase the number of the function, k.

X Y Z F 6 5 4 3 2 1 0 -1 -2 -3 -4 -5 (a)
X Y Z F 7 6 5 4 3 2 1 0 -1 -2 -3 (b)
The solution provided by the method agrees well with reference finite element solutions obtained by employing full-Newton-Raphson iterative schemes (following the same tendency than that shown for the beam bending problem). But, notably, the computed solution can be stored in a so compact form that an implementation of the method is possible on handheld devices such as smartphones and tablets. For instance, for Android-operated devices, an application has been developed (we call it iPGD and is freely downloadable from [START_REF] Bordeu | iPGD, an interactive PGD application for Android[END_REF]) that runs the model on a Motorola Xoom tablet running Android 3.0 without problems (only the surface of the model is represented for simplicity, given the limitations of the Android OS) (Figure 8). The 25-Hz feedback rate necessary for continuous visual perception is achieved without problems.

For more sophisticated requirements, such as those dictated by haptic peripherals, a simple laptop (in our case a MacBook pro running MAC OSX 10.7.4, equipped with 4-Gb RAM and an Intel core i7 processor at 2.66 GHz) is enough to achieve this performance. Feedback rates in the order of kilohertz are obtained without problems.

CONCLUSIONS

Model order reduction seems to play an important role in real-time simulation of soft biological tissues. In the last times, there have been a number of publications on the use of POD techniques to this class of problems. However, POD-based approaches seem to over-simplify models, and approaches other than those used in [START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF] do not reproduce properly the nonlinearity of soft tissues. In this paper, a new approach to the problem has been introduced. It is based on the use of PGD methods. This implies a complete change of paradigm, because PGD (in contrast to POD) do not need prior computer experiments to generate the snapshots needed to construct the optimal basis functions. The reduced approximation bases, in a separated form, are constructed on the fly.

PGD, on the contrary, operate in a two-stage approach. Firstly, a general meta-model is computed a priori once for life. During this intensive computation phase, the solution to the high-dimensional model is computed as a finite sum of separable functions. This compact solution, although not optimal, in general, provides a very light format to store the solution in the form of a metamodel that provides the solution to the problem for any parameter value. In this particular problem, parameters are chosen as the position of the contact force between organ and surgical tool (scalpel) and orientation of the load (thus rendering a problem defined in R 9 ). This meta-model is then evaluated under real-time restrictions very efficiently (e.g., reaching more than 1 kHZ in a MacBook pro laptop). In this paper, some benchmark examples have been given to justify the accuracy of the proposed approach. In particular, two different explicit algorithms for the time integration of the resulting equations have been proposed. These algorithms have shown to work well, although the development of more robust strategies of linearization of the weak form of the problem, based on the use of asymptotic numerical methods, is now being sought. This is part of our current effort in this research.

Figure 1 .

 1 Figure 1. Convergence of the tip displacement toward the reference solution as a function of the chosen time step.

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 2. Model for the beam bending problem.

Figure 5 .

 5 Figure 5. Finite element model for the human liver.

Figure 6 .

 6 Figure 6. Six first functions X k 3 .x/, k D 1, : : : 6, for the simulation of the liver.

Figure 7 .

 7 Figure 7. Six first functions Y k 3 .s/, k D 1, : : : 6, for the simulation of the liver. Note that, in this case, functions Y k .s/ are defined on the boundary of the liver only.

Figure 8 .

 8 Figure 8. An example of the implementation of the iPGD application for the liver problem.

ACKNOWLEDGEMENTS

This study was supported by the Spanish Ministry of Economy and Innovation (contract/grant number: CICYT DPI2011-27778-C02-01) and the European Union, Seventh Framework Program (Vision Advanced Infrastructure for Research (VISIONAIR) project, contract/grant number: EU-FP7-STREP-262044).


Achal Agrawal

Samson Lasaulce

Olivier Beaude

Raphael Visoz

A framework for optimal decentralized power control with partial CSI

A Framework for Optimal Decentralized Power Control with Partial CSI Achal Agrawal 1 , Samson Lasaulce 1 , Olivier Beaude 1 and Raphael Visoz 2 1 L2S (CNRS -CentraleSupélec-Univ. Paris Sud 11), France 2 Orange Labs, Issy-les-Moulineux, France 1 {agrawal, lasaulce, beaude}@lss.supelec.fr, 2 raphael.visoz@orange-ftgroup.com Abstract-In this paper, we use a recent information theoretical result to develop a general framework for finding optimal power control policies in the case of interference channels. The aforementioned result characterizes the achievable payoffs for an N -Agent (transmitters in our application) coordination problem with a certain information structure. We then provide an algorithm which exploits the characterization of achievable payoffs by conditional probability distributions to find optimal decision functions for the transmitters. Due to its general nature, the developed framework is conducive for applications to diverse scenarios in wireless communications. In this article, we restrict our attention to the case of decentralized power control in interference channels for different utility functions namely sum-rate, sum-energy and sum-goodput. Our approach has the following salient features: 1) The method proposes optimal power control functions for any given utility function as opposed to ad-hoc solutions for different utilities proposed in the literature, and 2) Noise in the channel estimation is taken into account, thus providing robust optimal solutions.

I. INTRODUCTION

One of the key contributions of this paper is to exploit the recent theorem derived in [START_REF] Larrousse | Coordinating partiallyinformed agents over state-dependent networks[END_REF] to find power control functions which may exploit the available knowledge optimally (in the long run). As an example, we treat the well studied problem of optimal power control in interference channels ( [START_REF] Zander | Performance of optimum transmitter power control in cellular radio systems[END_REF], [START_REF] Yates | A framework for uplink power control in cellular radio systems[END_REF], [START_REF] Stanczak | Fundamentals of resource allocation in wireless networks: theory and algorithms[END_REF]) and show that not only does our framework attain various payoffs identical to the state of the art, it also provides the optimal power control functions in diverse scenarios. The scenarios considered resemble those of [START_REF] Gohary | Robust IWFA for open-spectrum communications[END_REF], [START_REF] Anandkumar | Robust rate maximization game under bounded channel uncertainty[END_REF]. These works propose ad hoc solutions for specific utilities but not a generic framework which directly provides a power control function that aims at exploiting the available and arbitrary information as well as possible. Indeed our approach is very general and can be used to treat more complex and interesting scenarios by considering different kinds of channel state informations available at transmitters, different common objectives, robust power control taking the noisy communications into account, vectorial optimization etc.

The generality of the approach is by virtue of exploitation of some recent information theoretical results [START_REF] Larrousse | Coordinating partiallyinformed agents over state-dependent networks[END_REF]. A similar approach, albeit with a different information structure and for only 2 transmitter-reciever pairs was proposed in [START_REF] Larrousse | Implicit coordination in 2-agent team problems. application to distributed power allocation[END_REF]. The information structure considered in this paper is more adapted for the application proposed under realistic assumptions, the important assumptions being 1) No knowledge of future realizations of the nature state w.r.t. which the agents are coordinating, 2) Any number of agents with only a partial and noisy information about the current nature state, and 3) Decision functions must be found in a decentralized manner for having more robustness leading to simpler design architectures.

Indeed, for the case of power control in wireless interference channels, these assumptions have been paid some attention. The necessity of minimizing the information required at the transmitters for coordination led to approaches proposed in [START_REF] Paul De Kerret | Best-response team power control for the interference channel with local CSI[END_REF]), [START_REF] Gjendemsjø | Binary power control for multi-cell capacity maximization[END_REF], [START_REF] Gjendemsj | Binary power control for sum rate maximization over multiple interfering links[END_REF]. They consider only local channel state information at the transmitters (CSIT). This is also important for reducing the complexity of the optimizations to be performed. For sum-rate maximization, the scenario considered by us is comparable to a single-carrier version of the iterative water filling algorithm (IWFA) [START_REF] Scutari | Asynchronous iterative water-filling for gaussian frequency-selective interference channels[END_REF]. In the case of sum-energy utility, channel inversion strategy was proposed by [START_REF] Goodman | Power control for wireless data[END_REF], whereas the multi-carrier version of the problem was considered by [START_REF] Meshkati | A game-theoretic approach to energy-efficient power control in multicarrier CDMA systems[END_REF]. Decentralization of the power control, or for that matter for any wireless network design, is the reason why game theory has been extensively applied to such problems [START_REF] Lasaulce | Game theory and learning for wireless networks: fundamentals and applications[END_REF]. We thus compare our results to game theoretic equilibria analysis as most of the decentralized optimization literature only guarantees convergence to Nash equilibria.

The article follows the structure from more general results to specific application of the framework proposed. To this effect, Sec. II describes the performance limits of a general coordination problem characterizing the optimal achievable payoffs. The characterization lays the foundation for an algorithm proposed in Sec. III for finding the optimal decision policies of the coordinating transmitters. While we continuously illustrate what the developed framework corresponds to for wireless applications, in Sec. IV-A we describe the specific problem of power control in interference channels. A numerical analysis is conducted for diverse scenarios and the results are compared to existing literature in Sec. IV-B. Concluding remarks are provided in V II. PERFORMANCE LIMITS OF AN N -TRANSMITTERS INTERFERENCE COORDINATION PROBLEM Consider a coordination problem where N transmitterreceiver pairs are trying to coordinate the power emitted by them X i ∈ X i , i ∈ {1, ..., N } with respect to a nature state represented by X 0 ∈ X 0 to optimize a common payoff function w(x 0 , x 1 , x 2 ) over a long time period T . The nature state X 0 for applications to wireless communications is the channel gain matrix g ij = |h ij | 2 , with g ij being the channel gain of the link between transmitter i and reciever j. X 0 = (g 11 , ..., g N N ) is a random state which affects the common payoff function for the system and is not controlled by the coordinating transmitters. The realizations of the nature state X 0,t at each time instant t are i.i.d. and follow a distribution ρ. In wireless communications ρ is typically an exponential distribution for each channel gain g ij . Our aim is to characterize all the achievable expected payoffs under a certain information structure over the large time period T , T → ∞.

Firstly, we need to define the information structure under consideration. At every instant t, all the transmitters are assumed to have an image or a partial observation S i,t of the nature state X 0,t with respect to which they are coordinating. In the case of our application, this could be knowledge of only local CSIT, i.e. transmitter i observes a noisy version (in general) of only the direct link channel gain g ii . One could imagine other kinds of information available at the transmitters: for example, transmitter i observes all the links g ji ∀j. The observations S i,t are assumed to be generated by a memoryless channel whose transition probability is denoted by i (S i,t |X 0,t ). All transmitters have to choose the power emitted by them X i,t based on this information received. Formally, the sequences of decision functions for each agent i, σ i,t is defined as:

σ i,t : S t i × U → X i (1) (s i (1), s i (2)...s i (t), u(t)) → x i (t) (2)
where S t i = S i (1) × S i (2)... × S i (t) is the discrete observation alphabet till the instant t and X i is the power emitted for transmitter i with |S i |, |X i | < ∞. U t is the alphabet of the auxiliary variable U which is discussed in more detail later. s i (t), u(t) and x i (t) are the realizations of the corresponding variables at instant t.

The problem is said to be decentralized as each transmitter chooses its power independently based on the information received by it. Since the channel gains are generated by a random process, and thus in general also the corresponding chosen power levels, the quantity to be optimized is the expected objective function:

E Q [w(X)] = x∈X w(x 0 , x 1 , ..., x N )Q(x 0 , x 1 , ..., x N ) (3)
where X = (X 0 , X 1 , ..., X N ), x = (x 0 , x 1 , ..., x N ), and X = X 0 × N i=1 X i . An important point to note is that since expectation is a linear operator, optimizing the expected payoff is equivalent to finding the optimal distribution Q(x 0 , x 1 , ..., x N ). However, the optimization problem is not so trivial as indeed there are certain restrictions on the distributions Q that are implementable given the imposed information structure. For doing that, one first needs to define the notion of an implementable distribution.

Definition 1 (Implementability). Let the information structure be as defined in [START_REF] Larrousse | Coordinating partiallyinformed agents over state-dependent networks[END_REF]. The probability distribution Q(x 0 , x 1 , ..., x N ) is implementable if there exist decision functions σ i,t such that as T → +∞, we have for all x ∈ X ,

1 T T t=1 Q X0...X N ,t (x 0 , ..., x N , ) → Q(x 0 , ..., x N ) (4) where Q X0X1...X N ,t = Q X1,...X N |X0,i × ρ is the joint distribution induced by σ i,t at stage t.
As seen before, the expected payoff is characterized by the probability distribution Q over all the variables that intervene in the payoff function. Thus, the time averaged expected payoff w is said to be achievable, if and only if the corresponding distribution Q is implementable. The following theorem characterizes precisely the distributions that are implementable under the information structure (1).

Theorem 1.

[1] Assume the random process X 0,t to be i.i.d. following a probability distribution ρ and the available information to the transmitters S i,t to be the output of a discrete memoryless channel obtained by marginalizing the conditional probability (s 1 , ..., s N |x 0 ) . An expected payoff w is achievable in the limit T → ∞ if and only if it can be written as:

w = x0,x1,...x N ,u,s1,...s N ρ(x 0 )P u (u)× (s 1 , ...s N |x 0 )× N i=1 P Xi|Si,U (x i |s i , u) w(x 0 , x 1 , ...x N ) (5)
where U is an auxiliary variable which can be optimized.

The auxiliary variable U is an external lottery known to the transmitters beforehand which can be used to achieve better coordination. An example for how this variable U can help coordinate better in the case of power control is shown in Sec. IV-A.

Theorem 1 helps us find the best achievable payoffs on an average with arbitrary information to the transmitters as long as it follows the structure (in terms of when the information is available to whom) defined in equation [START_REF] Larrousse | Coordinating partiallyinformed agents over state-dependent networks[END_REF]. However, this theorem does not provide optimal sequences of decision functions σ i,t for the time averaged expected payoff. In the following section, we shall give an algorithm which helps find a suboptimal solution for the decision functions. The suboptimality is due to the high complexity of a multilinear optimisation problem. Nonetheless, In Sec.IV-B we see that the solutions obtained using this algorithm match the state of the art solutions proposed for different utilities discussed.

III. ALGORITHM FOR FINDING SUB-OPTIMAL LOW-COMPLEXITY DECISION FUNCTIONS

As seen from Theorem 1, optimal performance for the coordination problem depends only on the vector of conditional probabilities (P X1|S1,U , ..., P X N |S N ,U ). Thus, it suffices to find an optimum vector of lotteries for every action X i possible to obtain the best achievable payoff. However, one can see that this task is computationally very demanding. Instead, we apply a suboptimal approach by searching for the conditional probabilities in a distributed manner, thus reducing the complexity of the search. To this effect, we apply sequential best response dynamics [START_REF] Lasaulce | Game theory and learning for wireless networks: fundamentals and applications[END_REF]. This consists of each transmitter choosing the best corresponding power given that others keep their powers constant, and each transmitter doing so sequentially within an iteration. Note that by doing so, the optimization problem becomes linear, and hence its solution lies at a vertex, i.e. either (P Xi|S1,U ) = 0 or (P Xi|S1,U ) = 1. Thus the search for conditional probability distributions simplifies to search for decision functions

f i : S i × U → X i (s i (t), u(t)) → x i (t)
An important point to note is that we search for stationary strategy, i.e. f i does not depend on the time instant t. This simplification is practical as it leads to simpler design structures and easier to find in terms of complexity. Also, the auxiliary variable U , which serves a coordination key, is exchanged offline beforehand and thus does not entail any signalling cost.

The procedure outlined below is performed many times till convergence is achieved, typically in a few iterations. The number of iterations required for convergence of course depends on the number of agents coordinating, but it scales up very slowly as best response dynamics normally converge very fast. Furthermore, since we are considering a common 'team' payoff w(x 0 , x 1 , ..., x N ), the convergence of best response dynamics is guaranteed [START_REF] Monderer | Potential Games[END_REF]. Proposition 1. Algorithm 1 converges always for a common performance criteria w(x 0 , x 1 , ..., x N ).

Proof: The result can be proved by induction or more generally by calling for an exact potential game property (the latter argument may hold in the more general case in which the transmitter have different performance criteria).

For further justification of the simplification made by finding functions f i instead of conditional probability distributions, note that the optimization problem is linear w.r.t. each component of the vector (P Xi|S1,U ) that we are searching. The problem is therefore multilinear and Algo 1 will converge to the vertex of the polytop defined by the constraints. Note however that we might lose out on finding the global optimum by doing so. Nonetheless, the power control functions obtained by the simplification are sufficient to compare with the existing literature which typically provide optimal control functions and not mixed randomized strategies.

To apply best response dynamics, we rewrite the expected utility given by Theorem 1 in the following way:

w = x0,x1,...x N ,u,s1,...s N ρ(x 0 )P U (u)× (6)
(s 1 , ..., s N |x 0 )× (7)

N i=1 P Xi|Si,U (x i |s i , u) w(x 0 , x 1 , ...x N ) = ai,bi,u θ ai,bi,u P Xi|Si,U (x ai |s bi , u) (8)
where a i , b i , u are the respective indices of x i , s i , u and

θ ai,bi,u   a0 ρ(x a0 ) i (s i |x a0 ) a-i w(x a0 , x a1 , ...x a N )× b-i i =i i (s b i |x a0 ) i =i P X i |S i ,U (x i |s b i , u)   P U (u) (9)
where a -i , b -i are the indices which represent a i , b i being constant, while all the other indices are summed over. Written in this form, the linearity of the problem w.r.t. any one of the variables P Xi|Si,U (x i |s i , u) is even more explicit. Note that we have also assumed the independence of the observation channels of all the agents, i.e.

(s|x 0 ) = 1 (s 1 |x 0 ) × . . . × N (s N |x 0 ).
For interference channels, this assumption is reasonable as all transmitters receive independent feedback of the channel gains.

An important point to note is that the transmitters can run this algorithm offline as all they need to know are the channel statistics ρ. This helps in using the optimal functions for the entire timeslot where the channel statistics remain constant right from the start as opposed to taking some time initially to find the optimal functions online. The latter is the case for some decentralized solutions proposed such as the algorithms based on water-filling techniques proposed in [START_REF] Gohary | Robust IWFA for open-spectrum communications[END_REF], [START_REF] Anandkumar | Robust rate maximization game under bounded channel uncertainty[END_REF].

IV. APPLICATION TO POWER CONTROL IN INTERFERENCE CHANNELS

In the previous sections, we developed the framework, giving intuition as to what the defined quantities might represent in diverse optimization problems of wireless communications. In this section, we concentrate on the problem of power control in interference channels [START_REF] Lasaulce | Game theory and learning for wireless networks: fundamentals and applications[END_REF] and provide numerical analysis of some special cases for the above problem.

A. System Model

Consider N single-antenna Transmitter-Receiver pairs Transmitter i , Reciever i , i ∈ {1, 2..., N } with a single band interference channel with the channel gains being g ij , i, j ∈ {1, 2..., N } and g ij ∈ G. G is considered to be a discrete set and represents the alphabet of the possible channel gains. A well accepted model of statistics for the channel gains g ij = |h ij | 2 is Rayleigh fading. In this model, due to central limit theorem, the real and the imaginary components of h ij Algorithm 1: Proposed decentralized Algorithm for finding a control function inputs : X i ∀i ∈ {0, ..., N }, w(x 0 , x 1 , ...

x N ) ∀x, ρ(x 0 ), S|X0 (s|x 0 ) ∀x 0 , f init i ∀i ∈ {1...N }, P U (u) output: f * i (s i , u i ) ∀i ∈ {1...N } Initialization -f 0 i = f init i , iter = 0, iter max = 100 while ∀i f (iter-1) i -f iter i 2 ≥ AND iter ≤ iter max OR iter=0do iter = iter+1;
forall the i ∈ {1, . . . , N } do forall the s i ∈ S i do forall the u ∈ U do forall the x i ∈ X i do

Find the minimum (or maximum) coefficient θ ai,bi,u for minimization (or maximization) of payoff using (9); Update the function

f iter i (s i , u i ) ∈ arg min ai θ ai,bi,u ;
end end end end end follow a normal distribution, and thus the channel gain g ij follows exponential distribution (|h ij | follows Rayleigh distribution). We also assume that all channel gain distributions are independent of each other and that their realizations for each timeslot t is i.i.d. The transmitters transmit at discrete power levels P i ∈ P i quantized uniformly (in dB) between [0, P max ]. The utility for the pair Transmitter i , Reciever i depends on their signal to interference plus noise ratio (SINR) γ i (P , g) where P is the vector with each component P i being the power emitted by the i th transmitter. More precisely,

γ i (P , g) = P i g ii σ 2 + i =j P j g ji (10)
where σ 2 is the noise variance. Without loss of generality, we'll consider σ 2 = 1, and change P max to regulate the signal to noise ratio (SNR). SNR is defined to be the ratio

P max /σ 2 .
We shall consider only local CSIT, i.e. the direct channel gain g ii available at each transmitter i, based on which they have to choose the power level to emit at. For evaluation of the solutions found using our method, we shall consider 3 different payoff functions, which have vastly different optimal solutions and show that our method gives the same results as the state of the art which employ different methods for each payoff. The following common payoffs considered are:

• Sum-rate: w R (P , g)

= i log(1 + γ i ) • Sum-goodput: w G (P , g) = i Ω(γ i ) • Sum-energy: w E (P , g) = i Ω(γ i ) P i
Typical functions for Ω are Ω(x) = e -c x [START_REF] Belmega | Energy-efficient precoding for multiple-antenna terminals[END_REF], where c > 0 or Ω(x) = (1 -e -x ) M where M ≥ 1 is the packet length [START_REF] Goodman | Power control for wireless data[END_REF]. In simulations, we supposed c = 1 as it only changes the optimal solutions by a multiplicative constant. Note that proposition 1 holds for the above payoffs as we are considering sum of individual payoffs which trivially satisfies the required condition for convergence.

The following table explicits the relation between the input variables for the algorithm 1 and the corresponding quantities for the application to power control defined above.

General model Power Control Application

Nature State -X 0 Channel state -(g 11 , ..., g N N ) Decision of Transmitter i -X i Power emitted P i ∈ P i w(x 0 , x 1 ...x N ) w R (P , g), w G (P , g), w E (P , g) ρ(x 0 )

N i=1 N j=1 e -g ij , g ij = E(g ij ) (s i |x 0 ) s i = g ii or s i = g ii f init i f i = Pmax∀s ∈ S

TABLE I CORRESPONDENCE BETWEEN THE GENERAL FRAMEWORK AND ITS APPLICATION TO POWER CONTROL

In Theorem 1, we introduced an auxiliary variable U which could help in achieving better coordination. For better intuition, we show how a simply constructed external lottery helps transmitters coordinate in the case of power control.

Consider only N = 2 transmitter-receiver pairs communicating over a single band interference channel while maximizing the sum-rate w R (P , g). From [START_REF] Gjendemsjø | Binary power control for multi-cell capacity maximization[END_REF], we know that the optimal solutions (P * 1 , P * 2 ) ∈ {(P max , 0), (0, P max ), (P max , P max )}. However, if transmitters try to maximise their individual utility log(1 + γ i ), the Nash equilibrium is (P max , P max ) [START_REF] Lasaulce | Game theory and learning for wireless networks: fundamentals and applications[END_REF]. Now consider the following probability distribution for a random variable U with the realizations u ∈ U, |U| = 3.

P U (u) =    P r(u = u 1 ) = α 1 P r(u = u 2 ) = α 2 P r(u = u 3 ) = α 3 = 1 -(α 1 + α 2 ) (11)
A large sequence of realizations for U are drawn with the probability distribution specified in [START_REF] Scutari | Asynchronous iterative water-filling for gaussian frequency-selective interference channels[END_REF] and all transmitters share this drawn sequence. Now consider a trivial function common to all transmitters ∆ : U → P which maps every realization of U , u i to one of three possible optimal power control vectors {(P max , 0), (0, P max )(P max , P max )}. Under this setup, one can do monte-carlo simulations for a large number of channel gain realizations, and find the optimal probabilities (α * 1 , α * 2 ) by doing an exhaustive search. This is done by searching over a discrete space α i ∈ [0, 1], i ∈ {1, 2. For the discretization, the step size is chosen to be 0.01 for obtaining reasonably precise 'optimal' probability values. In Fig. 1, we plot the average payoff for the coordinated policy using the pre-drawn lottery U in terms of SNR (which is equal to P max due to normalization of channel noise variance σ 2 ). We compare the average payoffs for the coordinated policy as well as for nash equilibrium and see that at higher SNR, it is indeed preferable to coordinate using U . This is not surprising at all, as at higher SNRs, the interference generated is too high, whereas time-sharing helps avoid this.

In the following section, we provide the simulation setup and discuss the power control functions given by our approach. We do not consider the pre-drawn lottery U in our analysis, although in general U will only help improve the coordination performance. This is because the solution where transmitters ignore the coordination key is always possible.

B. Numerical Analysis

1) Simulation Setup: We implement Algo. 1 for the case of power control in a single band interference channel [START_REF] Lasaulce | Game theory and learning for wireless networks: fundamentals and applications[END_REF]. To this effect, we use the correspondence established in Sec. IV-A. We consider the case when N = 2. Note that exponential distributions for channel gain statistics g ij are characterized only by their mean g ij . We shall consider standard scenarios where (g ii , g ij ) = (1, 0.3), i = j.

As mentioned before, the noise variance of the channel has been normalized to one. Thus, the SNR is regulated via P max = 20db. The possible power levels emitted by In particular, it is seen that for maximizing sum-energy, the obtained power control function exhibits a threshold under which transmission should not occur for a given transmitter.

transmitters (P i ) are uniformly quantized between P i ∈ [-20, 20] db , with 50 quantization points.

Since the setup considered is symmetric for all Transmitter -Reciever pairs, we plot the power control functions f i (s i ), where S i = G ii , for the transmitter i. Note that since we do not consider the auxiliary variable u from here on, the function f i is represented with only one argument s i .

2) Results: In Fig. 2, we plot the power control functions f i obtained using the Algo 1 for the case N = 2 against g ii the direct channel gain between the i th Transmitter-Receiver pair. We do not consider noise in the channel estimation for the moment, and thus in this case s i = g ii . Three different utilities are considered, sum-rate, sum-energy, and sum-goodput.

For the case of sum-rate, it is known that binary power control P i ∈ {P min , P max } is optimal for 2 Transmitter-Reciever pairs [START_REF] Gjendemsjø | Binary power control for multi-cell capacity maximization[END_REF]. Moreover, as shown in [START_REF] Paul De Kerret | Best-response team power control for the interference channel with local CSI[END_REF], optimal power control functions with only local CSIT (g ii ) amounts to P i (g ii ) = P min if g ii ≤ g * and P i (g ii ) = P max otherwise.It is reassuring to find that our results verify this. Thus, in the case of sum-rate with local CSIT, we find exactly the same results as state of the art.

In the case of sum-energy, we see that all the available power P max is not used. Indeed, even in the case of only 1 Transmitter-Reciever pair, the optimal power control function is c g11 . We see that there is also a threshold value of g ii below which, P i = P min . Above the threshold value, the function is similar to the optimal solution obtained in the case of only 1 Transmitter-Receiver pair. The threshold function is also seen as solution in the case of Sum packet rate. This is because the utility function is not monotonous w.r.t. SIN R.

To the best of our knowledge, the power control functions found for sum-energy are new. Solutions in literature do not propose thresholding of power control functions to help reduce interference in case of bad direct channel gains. Admittedly, this requires greater analysis and comparisons with the state of the art. N=2 Fig. 3. Our algorithm reveals the shape of good power control functions in the presence of interference. In contrast with related works on energyefficiency, our work shows that thresholding is required to manage interference efficiently.

To investigate the function for sum-energy further, in 3 we compare the function obtained using our algorithm with that proposed by [START_REF] Goodman | Power control for wireless data[END_REF]. We see that unlike [START_REF] Goodman | Power control for wireless data[END_REF], we find a threshold below which emitting no power is more optimal. We also see that the presence of noise is makes the power control function more uniform. Our approach is shown to be much more robust to noise than one proposed by [START_REF] Goodman | Power control for wireless data[END_REF] as the function 1/g ii is very sensitive to noise at low values of g ii . Relative performance gain (w f -w N ash /w N ash in % ) of our algorithm when compared to average payoffs obtained by Nash equilibrium for sum-rate and sum-energy. The curve for sum-energy saturates as at high SNR, Pmax is not utilized as power emitted is much lower.

In Fig. 4, we compare the performance of the power control function with that of Nash equilibrium for different SNRs in the case of sum-energy and sum-rate. Thus, on the y -axis we see the relative performance gain (w fw N ash /w N ash in % ) of our algorithm when compared to average payoffs obtained by Nash equilibrium. Current state of the art propose solutions which converge to the Nash equilibrium. Thus, the solutions are generally not optimal while maximizing common utilities that are sum of the individual utilities. The evaluations for obtaining average expected payoffs was done for 10 6 channel realizations. They were drawn according to the distribution N i=1 N j=1 e -g ij . The channel gain statistics g ij were the same as the ones used for running the algorithm.

We see that for sum-energy, there is no performance gain at higher SNRs. This is because as seen from Fig. 2, the optimal power used is very low for all g ii . However, we see that when compared to nash equilibrium, we obtain a significant gain of around 55 %. In the case of sumrate, performance obtained matches with [START_REF] Paul De Kerret | Best-response team power control for the interference channel with local CSI[END_REF] since the power control functions obtained are identical. Not suprisingly, at higher SNRs, Algo 1 does much better as Nash equilibrium, which is all transmitters transmitting at P max generates too much interference.

One of the advantages of our framework is that it provides optimal power control functions even for noisy channel estimates. We illustrate this in Fig. 5 where power control functions for different levels of noise in channel estimation are plotted. This noise simulates the error in estimation due to noise in feedback transmission, or just simply statistical estimation error. The noise for the simulations is modeled to be gaussian, i.e. s i = g ii + z where Z ∼ N(0, σ 2 z ). As expected, the power control functions become more uniform at higher noise levels, as the information received is less reliable, and thus transmitters emit at a power level which maximizes the utility for an average case. Influence of channel estimation noise on the power control functions. We see that as noise increases, the optimal power emitted becomes more uniform. This is intuitive as higher the noise, higher the uncertainty in the observation leading to same power emitted for all channel gains observed in the asymptotic case.

V. CONCLUSIONS

The proposed framework was shown to be relevant in diverse scenarios of single band interference channel for finding optimal power control functions. Moreover, the power control functions depended only on local CSIT, thus having the merit of being implementable in a completely decentralized manner. Also, the solutions obtained take noise in the estimation of the channel gain into account. All the above features illustrate the generality of our approach in tackling problems of power control for maximizing sumutility functions.

However, the framework can be exploited further for tackling other problems in wireless communications as well. For example, one could consider the problem of power allocation in a multi-band interference channel. Also, the auxiliary variable U was not exploited, which in general will only make the solution more optimal. One could also consider different information available at transmitter i: g ij ∀j. In future, we hope to analyze more complex problems and show the merit of the proposed framework in finding decentralized strategies for transmitters to cooperate for optimizing common utilities.

Fig. 1 .

 1 Fig.1. The figure shows the performance in terms of sum-rate with : the proposed coordination mechanism using a public key shared beforehand; the single-carrier counterpart of the iterative water-filling algorithm. Without additional signalling, a significant gain is obtained.

Fig. 2 .

 2 Fig.2. The figure represents the power control functions provided by Algorithm 1 for the three performance metrics under consideration in this paper. In particular, it is seen that for maximizing sum-energy, the obtained power control function exhibits a threshold under which transmission should not occur for a given transmitter.

  Fig.[START_REF] Stanczak | Fundamentals of resource allocation in wireless networks: theory and algorithms[END_REF]. Relative performance gain (w f -w N ash /w N ash in % ) of our algorithm when compared to average payoffs obtained by Nash equilibrium for sum-rate and sum-energy. The curve for sum-energy saturates as at high SNR, Pmax is not utilized as power emitted is much lower.

  Fig. 5.Influence of channel estimation noise on the power control functions. We see that as noise increases, the optimal power emitted becomes more uniform. This is intuitive as higher the noise, higher the uncertainty in the observation leading to same power emitted for all channel gains observed in the asymptotic case.

ACKNOWLEDGMENT

Authors of this paper and particularly A. Agrawal would like to thank project LIMICOS -ANR-12-BS03-0005 for financing the project.