Olivier Pinaud
email: olivier.pinaud@g2elab.grenoble-inp.fr

Olivier Chadebec
email: olivier.chadebec@g2elab.grenoble-inp.fr

Laure-Line L Rouve
email: line.rouve@g2elab.grenoble-inp.fr

Jean-Louis L Coulomb

Jean-Michel M Guichon
email: michel.guichon@g2elab.grenoble-inp.fr

Andrea Vassilev
email: andrea.vassilev@cea.fr

A Bayesian Approach for Spherical Harmonic Expansion Identification: Application to Magnetostatic Field Created by a Power Circuitry

Keywords: Bayes theorem, inverse problem theory, magnetostatics, random variables propagation, spherical harmonics

This paper deals with the use of the Bayesian approach to inverse an underdetermined magnetostatic problem based on spherical harmonic expansion. Identification of the spherical harmonic coefficients is helped thanks to some a priori information. This information comes from a numerical model statistically studied to define an average-state vector and a covariance matrix. The whole approach is applied for the study of the magnetostatic field inside an electric vehicle, created by its power circuitry. It demonstrates the strength of merging a priori information and measured information in order to obtain an efficient identification of magnetic sources created by a complex set of conductors.

I. INTRODUCTION

T

In the case where propagation phenomena can be neglected, when associated frequencies are extremely low, from dc to a few hundreds of kilohertz, different software can numerically compute the magnetic field created by complex systems [1], [2], or [3] for example. Even if the modeling remains necessary to correctly handle the physical phenomena and characterize the magnetic field, it is obvious, and unfortunately very often verified by experiment, that measured fields differ from computed ones. The main reason is that some information remains always unknown for a structure: part of the geometry or the sources, unwanted unbalances, defaults, etc. As a consequence, direct modeling cannot give the whole representation of the magnetic fields. Moreover, this method may become totally unrealistic when several complex structures must be modeled to take into account all the sources, onboard a vehicle for example.

The second way to evaluate the electromagnetic fields is mainly based on experiment. Many sensors are spread either around the studied structure or the studied area to assess the radiated field levels and frequencies. In order to limit the number of sensors, these measurements can be combined with models that are mathematical expressions of the fields, like spherical harmonic expansions [START_REF] Hobson | The Theory of Spherical Harmonics[END_REF]. An inverse problem must be solved to identify the parameters of the harmonic model thanks to measurements [START_REF] Legris | Identification de l'état magnétique d'un système ferromagnétique à partir de mesures du champ proche[END_REF]- [START_REF] Vincent | Identification of equivalent multipolar electromagnetic sources by spatial filtering[END_REF]. Then, the field can be estimated at any point of the studied zone. Nevertheless this approach remains invasive, because many sensors must be located preferably on sphere surfaces [START_REF] Rouve | Optimal magnetic sensors location for spherical harmonics identification applied to radiated electrical devices[END_REF], either all around the sources or the studied area. For the specific application of electric vehicles for example, magnetic field is studied inside an area surrounded by sources. Obviously, because of the passengers, it is clearly impossible to set as many sensors anywhere in the compartment.

It is why this paper proposes an original method that enables magnetic field evaluation by merging together experimental information (magnetic field measurements) and a priori information (direct modeling of the main knowledge of the sources). The developed approach tries to combine these information's to obtain a better evaluation of the magnetic field than using each of them separately. To do so, a Bayesian approach is used to solve the inverse problem based on the field harmonic expansion. All this development is based on previous work [START_REF] Vassilev | Magnetic field exposure assessment in electric vehicles[END_REF], where magnetic field of several electric cars has been characterized from dc up to 10 MHz. Despite a quite complex spectrum, for the sake of science progress, we have chosen to study only the magnetic field at very low frequency (< 1 Hz). This choice has led to the numerical modeling of the main ferromagnetic parts HANKS to the latest developments of power electronics leading to higher and higher compactness, the improvements of electric energy storage and its availability, increased by onboard production, electricity has been spread as a major energy vector in the transport domain, as well as on sea, air, or land. Electric vehicles like buses, trams, or cars, are more and more common in our everyday life. In order to carry and transform the electric power, sizeable currents circulate through the vehicles, in the presence of magnetic materials. As a consequence, the question of the level of the electromagnetic fields created inside such vehicles can arise. This question may be motivated for many reasons, the main ones being electromagnetic compatibility or human exposure.

To address this general problem of evaluation of electromagnetic fields created by electric systems, two approaches can be found. The first one consists in proposing a model of the studied structure, and once geometry and physical sources are described, in computing the associated fields. The method seems simple but requires in fact a lot of knowledge, research, and development. and of the dc power circuitry of an electric vehicle [START_REF] Pinaud | Forward model computation of magnetostatic fields inside electric vehicles[END_REF]. This has required several month of work and demonstrates the high complexity of such an approach. In this paper, only a simplified shape of the power circuitry is used to develop the whole proposed approach combining a priori modeling of the sources and few measurements of the magnetic field.

In a first part, the approximation of the magnetic field by a truncated spherical harmonic expansion is presented [START_REF] Hobson | The Theory of Spherical Harmonics[END_REF], [START_REF] Durand | Electrostatique, Tome 2, Problèmes généraux conducteurs[END_REF]. Then, the identification of the expansion coefficients by an inverse problem resolution [START_REF] Tikhonov | Méthodes de résolution de problèmes mal poses[END_REF]- [START_REF] Tarantola | Inverse Problem Theory[END_REF] is discussed. A focus is especially made on the Bayesian theory [START_REF] Tarantola | Inverse Problem Theory[END_REF], [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory (Signal Processing series)[END_REF] ,which merges a priori knowledge of the sources with measurements of the field taking into account their uncertainties in order to solve the associated inverse problem. In a second step, the Bayesian implementation is described to estimate the a priori information; it means an a priori set of spherical harmonic coefficient and its associated covariance matrix, obtained from a numerical modeling of physical sources. Then, the method is numerically tested to show how a priori information enables two important aspects: the sensors number reduction and the identification solution stabilization with respect to measurement noise. The whole approach is finally experimentally validated on a laboratory mockup and some conclusions are drawn.

II. WRITING THE INVERSE PROBLEM

When studying magnetic field in air region, whatever the sources, general solution of Maxwell equations [START_REF] Durand | Electrostatique, Tome 2, Problèmes généraux conducteurs[END_REF], [START_REF] Durand | Magnétostatique[END_REF] leads to the spherical harmonic expansion of this magnetic field. Inside an internal zone (blue sphere of radius r 0 on Fig. 1 for example) surrounded by magnetic sources, the general expression of the magnetic field is given by

B (r, θ, ϕ) = -μ 0 4π ∞ n =1 +n m =-n a n,m • ∇ (r n • Y m n (θ, ϕ)) (1)
where Y m n are the spherical harmonic functions (derived from the Legendre polynomials), (r, θ, ϕ) are the spherical coordinates, "a n,m " the spherical harmonic coefficients and " ∇" the gradient mathematical operator.

When radius "r" gets smaller, high order "n" contribution becomes negligible because of the specific law "r n ". So, for a given radius r ≤ r 0 (centered on the blue sphere for example), the decomposition can be truncated at the order n = N max . Now, if the number of triaxial sensor is K, the forward linear mathematical problem that links harmonic coefficients to the magnetic field measured on sensors is

G • A = B mes (2)
where G is a (3K, N 2 max + 2N max ) matrix that depends on the harmonic decomposition order, the sensor position (r, θ, ϕ) and their number, A is a (N 2 max + 2N max ) vector of "a n,m " harmonic coefficients and B mes is a (3K) vector of the measured magnetic field components.

The inverse problem consists in finding a set of coefficients that fulfils the linear expression (2) but above all, that permits to determine the magnetic field in an area where sensors cannot be placed (here, it is inside the blue sphere of Fig. 1). The inverse problem can be solved when it is well posed; i.e., enough sensors (as many as unknowns) correctly located (well spread on a measuring sphere surrounding the studied area). Unfortunately, it is generally not the case [START_REF] Tarantola | Inverse Problem Theory[END_REF] and this is particularly true in our application. We can set only a few sensors, less than the number of unknowns (harmonic coefficients); so the problem is then underdetermined. Moreover, because the sensor's position is strongly constrained, it is impossible to have a uniform distribution on the sphere surface (i.e., all around the passenger for example); the problem is then also ill posed. Trying to solve (2) under these conditions leads to a pseudosolution according to the least-square criterion (3) but this result generally hardly fits outside the sensors positions.

A LSC = G t • G -1 • G t • B mes . ( 3
)
A probabilistic view of the inverse problem may help to find the most likely solution [START_REF] Tarantola | Inverse Problem Theory[END_REF]. Thus, the linear matrix system (2) is rewritten as a density probability of having B mes given a coefficients vector A, also called the likelihood probability: p(B mes |A). For linear discrete problems, the Gaussian distribution is commonly used [START_REF] Tarantola | Inverse Problem Theory[END_REF], [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory (Signal Processing series)[END_REF] and gives the following expression:

p (B mes |A) ∝ exp - 1 2 (GA -B mes ) t S -1 m (GA -B mes ) . ( 4
) With S m a (3K, 3K) covariance matrix linked to uncertainties due to measurements and/or inherent to the chosen mathematical model.

The expression within the brackets of ( 4) is a quadratic form and when the derivative is cancelled, the corresponding vector of coefficients (vector A containing all the "a n,m ") gives the maximum of the probability distribution. This solution is given by the maximum of likelihood (ML) estimator as follow [START_REF] Legris | Identification de l'état magnétique d'un système ferromagnétique à partir de mesures du champ proche[END_REF]. This solution "A ML " maximizes the probability (4) and is equivalent to a weighted least-square criterion.

A ML = G t • S -1 m • G -1 • G t • S -1 m • B mes . ( 5
)
Now, the Bayes theorem [START_REF] Tarantola | Inverse Problem Theory[END_REF] permits to bring some a priori information about the behavior of the solution by expressing the so called a posteriori probability density p(A|B mes ).

p (A|B mes ) ∝ p (B mes |A) • p (A) (6)
where p(A) is the probability function of the a priori information (general behavior of the solution) also classically described by a Normal law distribution

p (A) ∝ exp - 1 2 (A -A 0 ) t S -1 0 (A -A 0 ) . ( 7
)
With A 0 a (N 2 max + 2N max ) vector of a priori harmonic coefficients, obtained from an a priori modeling of the studied system (illustrated here by Fig. 1). S 0 is the A 0 associated covariance matrix and comes from a statistical study of the a priori numerical model. Finally, the maximum a posteriori (MAP) estimator leads to "A MAP " and maximizes [START_REF] Rouve | Application of the spherical harmonic model to identify the magnetic state of a system and to extrapolate its signature[END_REF].

A MAP = A 0 + G t S -1 m G + S -1 0 -1 G t S -1 m (B mes -GA 0 ) . ( 8
) This harmonic solution is clearly centered on the a priori information A 0 . This a priori solution is then modified proportionally to the difference between the actual field measurements and the a priori field levels (estimated from the a priori information). The modification factor naturally takes into account the measurement/model uncertainties trough the covariance matrix S m and the a priori information uncertainties through the covariance matrix S 0 . Finally, the a posteriori solution ( 8) is a natural balance between the a priori information and the measured information thanks to the defined uncertainties.

This allows us to merge together modeled information and measured information into an identification problem. In Section IV, ML solution and MAP solution will be compared to illustrate the benefits of bringing some a priori information into an underdetermined inverse problem. Before, the estimation of the a priori solution A 0 and covariance matrices S m and S 0 must be computed in order to determine the solution A ML (5) and A MAP [START_REF] Rouve | Optimal magnetic sensors location for spherical harmonics identification applied to radiated electrical devices[END_REF].

III. MODELING THE A Priori INFORMATION

Only for illustration purpose, we will work with a simplified model (see Fig. 4) of the power circuitry of an electric vehicle studied in previous work [START_REF] Pinaud | Forward model computation of magnetostatic fields inside electric vehicles[END_REF]. Let us assume some a priori information about this circuitry that can be modeled (described) like on the Fig. 1. Red lines are power cables and black lines represent electric current distribution inside batteries packs. The difference between the studied circuitry (see Fig. 4) and the a priori modeled one (see Fig. 1) is deliberate. This is to clearly demonstrate the benefits of bringing some knowledge into an identification problem. From this model, we will first compute the corresponding harmonic coefficients "a nm " using some properties of the spherical harmonic functions. This will give us the a priori solution defined by the vector A 0 . Then, the associated covariance matrix S 0 will be computed thanks to the propagation of random variable.

A. Computing the A Priori Harmonic Solution

The Biot-Savart formula allows the exact computation of the magnetic field created by the electric current flowing inside conductors. Nevertheless, the a priori solution to bring into our Bayesian approach is the set of harmonic coefficients "a nm " (contained in A 0 ) that describes the average harmonic state of the magnetic field anywhere inside the studied area (inside the blue sphere for example).

To obtain these coefficients, a useful property of spherical harmonic functions is used [START_REF] Hobson | The Theory of Spherical Harmonics[END_REF], [START_REF] Durand | Magnétostatique[END_REF]. Because these functions describe mathematical orthogonal basis, their scalar product (on a sphere, whose radius is r 0 for example) is always zero except when functions are the same

Y m n (θ, ϕ) • Y m n (θ, ϕ) • dS = 4π • r 2 0 2n + 1 . ( 9
)
As a consequence, the scalar product [START_REF] Pinaud | Forward model computation of magnetostatic fields inside electric vehicles[END_REF] between the radial magnetic field component B r [START_REF] Durand | Electrostatique, Tome 2, Problèmes généraux conducteurs[END_REF] and a given spherical harmonic function Y m 0 n 0 , gives a value proportional to the corresponding spherical harmonic coefficient [START_REF] Durand | Magnétostatique[END_REF].

Int = B r • Y m 0 n 0 (θ 0 , ϕ 0 ) • r 2 0 sin (θ 0 ) dθ 0 dϕ 0 (10) B r (r 0 , θ, ϕ) = -μ 0 4π ∞ n =1 +n m =-n a n,m • n • r n -1 0 •Y m n (θ, ϕ) (11)
Int = -μ 0 • a n 0 ,m 0 • n 0 • r n 0 +1 0 2n 0 + 1 • δ n 0 • δ m 0 . (12)
Finally, the radial magnetic field component is exactly computed on the blue sphere surface, over a regular grid for example, by the Biot-Savart law. Then, the numerical integration of (10) allows us to determine the harmonic coefficients values thanks to [START_REF] Durand | Magnétostatique[END_REF]. This gives a set of harmonic coefficients allowing describing the magnetic field everywhere inside the defined blue sphere. If interested in another area, simply define another sphere, which will lead to another set of coefficients.

B. Propagation of Random Variable

Once the a priori solution is computed (assumed harmonic average state definition contained in A 0 ), the uncertainties must be defined to estimate their influence on the solution.

In our problem, we assume that uncertainties concerning the source mainly come from the power circuitry geometry, as electric current may be directly measured. As a consequence, each coordinates of every coil point describing the inductor geometry is treated as a random variable with a mean value "μ" (defining the presented a priori model of Fig. 1) and a standard deviation "σ μ ."

We must determine how geometric uncertainties affect the variations of the spherical harmonic coefficients. The difficulty lies on their nonlinear link with the physical description. Nevertheless, it is not absolutely necessary to know the mathematical expression between inputs (physical parameters like electric currents and their position) and outputs (spherical harmonic coefficients). To do so, the Monte Carlo approach (MC) is well known and allows to estimate the standard deviation of the outputs, but its main drawback is its computational cost. Several thousands of draws must be made that include the computation of B r (by the Biot-Savart law in this example) and the numerical integration [START_REF] Pinaud | Forward model computation of magnetostatic fields inside electric vehicles[END_REF]. If the radial component B r is not simply computed by a Biot-Savart law and must be obtained from a more complex model (including, for example, some magnetic and conductive materials), the MC computation time becomes far too high. To reduce the number of computations, the unscented transformation (UT) [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF]- [START_REF] Ferber | Conducted EMI of DC-DC converters with parametric uncertainties[END_REF] is an interesting alternative for the propagation of random variables through models.

This approach, developed by Hulmann and Julier, is based on the intuition that: "from a limited number of parameters, it should be easier to approximate the statistical moments than the probabilistic distribution itself." It is explained that the probability distribution function of a random variable "U" of dimension "p" (13), defined with a mean value "ū" ( 14) and a covariance matrix "S uu " [START_REF] Heinz | Regularization of Inverse Problems[END_REF], can be approximated by 2p + 1 symmetrical weighted value [START_REF] Tarantola | Inverse Problem Theory[END_REF].

U = [U 1 U 2 . . . U p ] T (13) ū = [μ 1 μ 2 . . . μ p ] T ( 14
)
S uu = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ σ μ 1 σ μ 1 μ 2 . . . σ μ 1 μ p σ μ 2 μ 1 σ μ 2 . . . σ μ 2 μ p . . . σ μ p μ 1 . . . σ μ p μ 2 . . . . . . . . . σ μ p ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (15)
U 0 = ū, W 0 = k (n + k) U i = ū + (n + k) S uu i , W i = 1 2(n + k) ( 16
)
U i+n = ū - (n + k) S uu i , W i+n = 1 2(n + k) . With i ∈ [1, p], (p + k) S uu i
is the ith line or column of the square root of the covariance matrix (p + k) S uu and W i the associated probabilistic weighting. The parameterk ∈ and allows a better precision for the evaluation of the output statistical moments whose order is above three. In this paper, only the first two statistical moments are looked for (the mean value and its covariance), and thus, k = 1.

As an example, Fig. 2 illustrates the approximation of the probability distribution function for a 1-D random variable (p = 1), defined with a mean value "μ x " and a standard deviation "σ x ". Now applied to our application case, let's assume that this p-dimensional random variable called "U" contains all the geometrical coordinate parameters defining the electric conductor geometry. Then, each single spherical harmonic coefficient "a nm " is linked to "U" through a nonlinear function "f nm " [a mathematical picture of equations ( 10) to [START_REF] Durand | Magnétostatique[END_REF]].

a n,m = f n,m [U] . ( 17
)
In consequences, each realization of "U" [defined according to [START_REF] Tarantola | Inverse Problem Theory[END_REF]] leads to an associated value for each studied "a nm " coefficient, and thus, the vector "A", containing all the a priori spherical harmonic coefficients [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF], becomes a q-dimensional random variable (q = N 2 max + 2N max ) with a mean value "A 0 = ā" and a covariance matrix "S 0 ."

A = [a 1,0 a 1,1 a 1,-1 . . . a N m a x ,-N m a x ] T . ( 18
)
The mean value and the covariance matrix of this qdimensional random variable are then computed as follow:

A 0 = ā = 2p i=0 W i • A i S 0 = 2p i=0 W i • (A i -ā) (A i -ā) T . ( 19
)
Finally, compared to the MC stochastic method, this approach deterministically selects computational points, called sigma points [START_REF] Tarantola | Inverse Problem Theory[END_REF], depending on the statistical moments of the input random variables, in order to compute the outputs values, and thus, to evaluate their statistical moments [START_REF] Julier | The scaled unscented transformation[END_REF]. For the studied case, because of Biot-Savart law (short time computation), both approaches (MC or UT) can be applied and results compared. But UT should more clearly demonstrate its strength for much time consuming computations: if model integrates some ferromagnetic parts for example.

C. Evaluation of the A Priori Solution Covariance Matrix

The whole a priori inductor geometry (see Fig. 1) is defined by 49 independent coil points. Let us assume that each of them is defined by three independent parameters: position toward x-axis, y-axis, and z-axis. Now, each parameter is treated as random with a mean value (μ) that defines the a priori "average state model" (i.e., Fig. 1) and a standard deviation (σ) that defines the possible range of variation of the coordinate. For information, the a priori (see Fig. 1) and the studied geometry (see Fig. 4) are about 1.7-m long, 0.5-m wide, and 0.1-m height. This corresponds approximatively to a real vehicle power circuitry at scale 1:2. Then, the standard deviation is defined of 0.03 m toward x-axis and y-axis and of 0.01 m toward z-axis.

As a consequence, the associated random variable 'U" (containing all the coil points coordinates) is of dimension p = 147 and its covariance matrix "S uu " is only diagonal because every parameters are independent to each other. Then UT algorithm deterministically chooses 295 a priori electric conductors geometry using [START_REF] Tarantola | Inverse Problem Theory[END_REF]. These configurations are then used to compute the associated spherical harmonic coefficients [using (10) - [START_REF] Durand | Magnétostatique[END_REF]] and this finally leads to the covariance matrix evaluation using the formula [START_REF] Julier | The scaled unscented transformation[END_REF]. In comparison, for MC computation, each parameter follows a Gaussian probability distribution around its corresponding mean value and defined standard deviation. The electric conductor geometry then randomly varies for each spherical harmonic coefficients evaluation using [START_REF] Pinaud | Forward model computation of magnetostatic fields inside electric vehicles[END_REF] to [START_REF] Durand | Magnétostatique[END_REF]. This is done thousands of times to finally compute the associated covariance matrix.

Both obtained covariance matrices are quite similar and lead to the same identification results: see the normalized harmonic coefficients distribution (see Fig. 3) obtained with the MAP estimator using the second nonnoised sensors configuration (reference to Section IV-B). This confirms the strength of using a limited number of selected computational points for evaluating the covariance matrix of the a priori spherical harmonic coefficients from the assumed physical uncertainties. On a standard PC (Intel Core i5-2.5Ghz, 4Go RAM, 32 bits), the UT algorithm required a computation of about 20 s to evaluate the covariance matrix while MC algorithm required about 100 s.

Nevertheless, only the covariance matrix evaluated by the MC approach will be used in the following parts because the main objective of this paper is to demonstrate the benefits of bringing a priori information into an underdetermined inverse problem.

IV. BAYES THEOREM NUMERICAL APPLICATION

In this section, ML and MAP estimators will be compared for illustrating the benefits of bringing an a priori solution into an underdetermined inverse problem. This is done through three key steps; a well-posed inverse problem (enough sensors compared to the number of unknowns, with a good distribution to have a good spatial magnetic field sampling), an ill-posed un- derdetermined inverse problem (not enough sensors with a bad distribution) and a noise sensitivity study is finally carried out for each of them.

It must be pointed out that all the results presented in this section are all obtained from numerical computations; experimental results will be discussed in Section V.

The study consists in using the identified spherical harmonic coefficients, either obtained from ML or MAP estimator, to compute a posteriori the magnetic field approximation on a validation path and compare the result with a reference computation obtained from Biot-Savart law.

The path of comparison runs from the bottom to the top on the surface of the yellow sphere (see Fig. 4). This latter is centered on the center of the blue sphere, but with a smaller radius, representing the whole studied area border.

The basic parts of the problem are as follows: a priori information (A 0 and S 0 ) comes from the a priori numerical model illustrated by Fig. 1 while the measured information comes from virtual sensors placed around the blue sphere on the model Fig. 4. The difference between the a priori geometric modeling (see Fig. 1) and the numerically measured one (see Fig. 4) is deliberate and quite clear; a large electric current loop is nearly cancelled on the right of the a priori model (see Fig. 1) and the parallel conductor are farther from each other on the measured model (see Fig. 4). The Biot-Savart computation of the magnetic field modulus on the validation path, running on the yellow sphere, confirms these differences (see Fig. 5). These differences are meant for illustrating the benefits of bringing a priori information into an underdetermined identification problem. The knowledge we have about a system is always partly wrong because of uncertainties; here, the a priori information wrongness is exaggerated.

For both a priori model and studied model (the virtually measured one), we know that the magnetic field on the blue sphere surface is very well approximated by a spherical harmonic expansion truncated at maximum order N max = 6.

This defines the maximum number of required spherical harmonic coefficients (i.e., 48 values) for having a good description of the magnetic field anywhere inside the studied area (inside Fig. 6. Nonnoised magnetic field modulus on the validation path after coefficients identification using the first sensors configuration. the blue sphere). Then, two nonnoised measurement configurations are presented first in order to evaluate the theoretical identification quality and a measurement noise sensitivity study is carried out to evaluate the solution stability.

A. Well-Posed Inverse Problem

The first test case is composed of 16 triaxes sensors well spread all over the blue sphere surface. So, the problem is well posed (as many measured information as unknowns to identify) and because measurements are nonnoised (exact computation with Biot-Savart law), the inverse problem is correctly solved. Fig. 6 compares the reference induction modulus (Biot-Savart computation on the validation path, yellow sphere surface) with the induction modulus approximation, obtained from the identified harmonic coefficients either after using ML estimator or MAP estimator.

This result demonstrates that the inverse problem succeeds as expected when enough sensors can be well placed within a low noised environment. Also, it must be pointed out here, the very little difference between the ML estimator and the MAP estimator: a quite natural result with this measurement configuration. Nevertheless, the noise sensitivity study will demonstrate the benefits of bringing some a priori information for the stability of the solution.

B. Ill-Posed Undetermined Inverse Problem

This second test case brings more complexity by decreasing the sensor number and by modifying their distribution because of some geometrical constraints for example. The problem is no longer well posed, it is now ill posed because of a nonoptimal sensors distribution, and is here underdetermined, because of a lack of measured information. Only six triaxes sensors are spread on the bottom of the blue sphere surface. This location would correspond to the passenger seat position for example. An ML estimator would hardly succeed because we know it is necessary to find all the 48 spherical harmonic coefficients while only 18 measurements are available. MAP estimator may give better results thanks to the a priori information brought, describing the assumed statistical behavior of the solution. As before, this configuration is nonnoised and the reference induc-Fig. 7. Nonnoised magnetic field modulus on the validation path after coefficients identification using the second sensors configuration. tion modulus (black curve) computed on the validation path is compared with the ones obtained either with ML or MAP coefficients estimator (see Fig. 7).

There is a clear difference between the ML estimator solution (only using measured information) and MAP estimator solution (using both a priori and measured information). An MAP estimator leads to a very good result all over the validation path (blue curve) contrarily to ML estimator (red curve), which only partly succeeds near the measurement sensors locations.

C. Noise Sensitivity Study

This part of the study is dedicated to the evaluation of the solution sensitiveness to noise. We assume perfect mathematical modeling (i.e., proper harmonic expansion truncation at order N max = 6) for the virtual measurement sensors location. Thus, only the diagonal of the covariance matrix S m is classically defined with the measurement noise variance. This noise is defined to follow the same probability density distribution on each sensor axis: S m is then proportional to identity matrix. The virtual measurement still comes from a Biot-Savart computation but the level is then randomly modified with a uniform probability density distribution over the range ±100 nT. This defines an average signal to noise ratio of about 20% over the blue sphere surface. Our objective here is to have a look on the statistical behavior of the inverse problem solution with respect to measurement noise. Then, the identification of harmonic coefficient is done a hundred times from randomly noised virtual measurements and the result variations are studied.

Let us remind that the observed result is the modulus of magnetic field, whose approximation is computed on the validation path and compared to the reference (Biot-Savart law).

Table I resumes the error for the well-posed inverse problem (first test case with 16 sensors, well spread all over the blue sphere, Section IV-A) while Table II concerns the ill-posed underdetermined inverse problem (second test case with only six sensors, gathered on the bottom of the blue sphere, Section IV-B). The error study consists in the computation of the root sum square (RSS) of the difference between the reference signal and the a posteriori harmonic expansion of the magnetic field modulus over the yellow sphere. Results clearly show the ML estimator solution sensitiveness to noise. As the S m matrix is proportional to identity matrix, it has no positive effect on the quality of the identification [see expression ( 5)]. In comparison, an MAP estimator gives far better results (better stability) because of the form of the solution (8). This estimation is a kind of balance between the a priori information about the unknowns (spherical harmonic coefficients) and measurements (magnetic field components) thanks to the use of the covariance matrices S 0 and S m . This observation is even more particularly true when the problem is ill posed and underdetermined; an ML estimator is highly noise sensitive while MAP estimator still leads to a very good solution and with an impressive stability.

V. BAYES THEOREM EXPERIMENTAL VALIDATION

For this last section, an experimental mock-up has been built (see Fig. 8) in order to confirm the previous results. Electric conductors have been set according to the numerically studied model (see Fig. 4). It must be pointed out that the real geometry naturally slightly differs from the model but this is not a problem in our approach because the a priori information inherently takes into account this difference. Compared to the numerical results discussed before, where sensors were placed over the identification sphere, experimental tests introduce a significant modification: sensors are now placed over two perpendicular planes. This configuration would correspond to a case where sensors are inserted into the passenger seat; this permits to fulfill the sensor localization constraint inside the vehicle (see Fig. 9). Two sets of Bartington fluxgate sensors are used: six triaxes sensors are placed on the two planes for the identification purpose (Mag-3IE) and other six triaxes sensors are placed inside the studied area (Mag-3MSES) to check the identification results (see Fig. 9). Table III gives the sensors characteristics. Only the magnetostatic field created by the electric power circuitry is studied (as in Section IV): so the measurement is done in three steps (with electric current OFF, ON, and OFF) in order to subtract the uniform Earth magnetic field. On each step, the time-mean value is computed in order to reduce the measurement noise (due to environment and sensor). Each measurement is done several times to confirm that the recorded level really corresponds to the magnetic field created by the electric current distribution. Then, the covariance matrix S m is arbitrarily set proportional to the identity matrix with a standard deviation of 1.10 -10 T.

Table IV gives the RSS error between the magnetic field modulus computed after the inverse problem and the measured one on the checking positions (blue points on Fig. 9).

First of all, it must be pointed out that experimental results are less good than the ones obtained before with the numerical tests (see Section IV-B), essentially because of a practical simplification linked to the sensors distribution. As sensors are placed on planes, they are a bit further from the blue sphere surface and the corresponding spherical harmonic expansion truncation (i.e., at a higher measuring radius) should be made with a higher order (N max > 6). Despite this fact, in order to keep the exact same a priori information as for the numerical illustration (see Section IV), the sixth-order harmonic expansion has been kept for identification. This introduces a well-known harmonic modeling misfit affecting the identification quality but let us remind that the objective here is to demonstrate the benefit of bringing some a priori information into an identification problem. To do so, an ML estimator which only uses measured information, is compared with the MAP estimator, which uses both a priori and measured information.

Despite this practical simplification, the presented experimental configuration allows confirming the benefits of merging a priori information and measured information. The identification result is far better when both sources of information are used than only the measured ones. The a priori magnetic field modulus has also been computed on checking sensors; as we can see, the a priori information presents a misfit of about 70% over these locations. Thanks to the MAP estimator, the misfit is clearly reduced to about 10% and is even better than ML estimator result, which misfit is about 37%. This permits to quantitatively measure how the Bayes approach is able to make the solution converge toward its real.

VI. CONCLUSION AND DISCUSSION

In this paper, the magnetostatic field has been studied inside an area. Thanks to the use of the spherical harmonic expansion, the whole magnetic field has been described, inside a sphere, with a very few parameters whatever the number of sources or their complexity. When neither direct modeling nor the measurements can succeed alone in describing the whole magnetic field distribution because of some practical constraints, mixing these two aspects into a statistical approach of the inverse problem allows better results.

The direct modeling (i.e., the a priori information) has been statistically studied to assess the confidence level to put on. To do so, the MC algorithm has been used quite naturally, but also a less known one called UT. Thanks to the selection of a few computational points, this approach is far more lightweight than an MC algorithm for the evaluation of the outputs covariance matrix. The UT approach underlies toward some good perspectives for the propagation of uncertainties through complex numerical models.

The a priori information has been modeled and inserted into the Bayesian approach of the inverse problem. This allows to drastically reducing the number of required sensors but not only. Numerical illustrations have also shown promising results about the quality and the stability of the identification with respect to noise. Even if the first experimental results presented are not so excellent (compared to the numerical ones), because of some simplifications due to practical constraints, they are still very good, which demonstrates the strength of the developed approach.

The whole method is applied here for the study of a magnetostatic field distribution created by an electric power circuitry but thanks to the use of the spherical harmonic expansion of the magnetic field, it could easily take into account some ferromagnetic material as their contribution can be simply seen as another source, linearly or not dependent of the first ones. Finally, the principle is deeply general and could be widely applied to a lot of other engineering domains, promising a lot for future development.

Olivier Pinaud received the electrical engineering degree from ENSE3, Grenoble Institute of Technology, Grenoble, France, in 2010 and the Ph.D. degree in electrical engineering from University of Grenoble, Grenoble, in 2014.

Since 2011, he has been working with Grenoble Electrical Engineering Laboratory (G2Elab), Grenoble.

Olivier Chadebec received the Diploma degree in electrical engineering and the Ph.D. degrees from the Grenoble Institute of Technology, Grenoble, France, in 1997 and 2001, respectively.

He is currently a CNRS Senior Researcher (directeur de recherche CNRS) and leads the "Models, Methods and Methodologies Applied to Electrical Engineering" research group of G2ELab, University of Grenoble, Grenoble. He has been a Visiting Scientist at the Technology Centre, Federal University of Santa Catarina, Florianópolis, Brazil, in 2012-2013. His research interests include computational electromagnetism (finite element and integral methods), equivalent magnetic sources identification by solving inverse problems, and low magnetic fields measurements.

Laure-Line Rouve received the Ph.D. degree from Grenoble Institute of Technology, graduate degree from the Electrical Engineering School of Grenoble, Grenoble Institute of Technology, Grenoble, France, in 1993, and the Ph.D. degree in electrical engineering in 1996.

She has been working as a Research Engineer in the team of Low Magnetic Field, Grenoble Electrical Engineering Laboratory, Grenoble, since 1997. She is a Specialist of characterization and reduction of low magnetic fields, thanks to modeling and measurement. A first part of her research is dedicated to ship magnetic silencing, with a special focus on magnetic treatment and magnetoelasticity phenomena. Another part of her research deals with the stray field generated by electrical machines with applications to stealth and fault diagnosis. Since 2008, she has been the in charge of the technical team of Low Magnetic Field and its specific facility of Magnetic Metrology Laboratory in Low Field.

Jean-Louis Coulomb received the electrical engineering degree in 1972, the "Docteur de 3°cycle" degree in 1975 and the "Docteur d'Etat" degree in 1981 from the Institut National Polytechnique de Grenoble, Grenoble, France.

He is currently a Professor Emeritus at the "Ecole Nationale Supérieure d'Ingénieurs pour l'Energie, l'Eau et l'Environnement," Grenoble Institute of Technology, Grenoble. He works also at the Grenoble Electrical Engineering Laboratory, where he was the former Head of the "Low Magnetic Fields" team.

Jean-Michel Guichon received the Diploma degree in electrical engineering and the Ph.D. degree from Grenoble Institute of Technology, Grenoble, France, in 1998 and 2001, respectively.

He is an Associate Professor at the Université Joseph Fourier, Grenoble. He has been with the Grenoble Electrical Engineering Laboratory (G2ELab), St. Martin d'Hères, France, since 2003, working in the field of power electronics. He is involved in the development of design tools for power electronics, in particular, electromagnetics modeling. Andrea Vassilev received the mechanical engineering degree, in 1996, from Ecole Centrale de Lyon, Écully, France.

From 1996 to 2001, he was with PSA Peugeot Citroen, as a Mechanical Simulation Engineer. Since 2001, he has been with CEA-LETI, Grenoble, France, a research institute for electronics and information technologies.

Fig. 1 .

 1 Fig. 1. A priori electric current distribution.

Fig. 2 .

 2 Fig. 2. Example of probability density approximation of a 1-D random variable.

Fig. 3 .Fig. 4 .

 34 Fig.3. Identified spherical harmonic coefficients with MAP estimator using either MC or UT method for S 0 evaluation.

Fig. 5 .

 5 Fig. 5. Biot-Savart computation of the magnetic field modulus on the yellow sphere surface from either a priori model or measured model.

Fig. 8 .

 8 Fig. 8. Visualization of the experimental mock-up.

Fig. 9 .

 9 Fig. 9. Visualization of experimental sensors distribution.

TABLE I RSS

 I ERROR BETWEEN REFERENCE SIGNAL AND COMPARED SIGNAL FOR THE WELL-POSED INVERSE PROBLEM

				±100 nT Noise Applied
	Identification Method	No Noise (Fig. 6)	Average (100 values)	Maximum (100 values)	Minimum (100 values)
	Maximum of Likelihood	0.021%	17.858%	171.187%	0.240%
	Maximum A Posteriori	0.015%	0.177%	0.712%	0.031%
	Given in percentage of the reference signal.

TABLE II RSS

 II ERROR BETWEEN REFERENCE SIGNAL AND COMPARED SIGNAL FOR THE UNDERDETERMINED INVERSE PROBLEM

				±100 nT Noise Applied
	Identification Method	No Noise (Fig. 7)	Average (100 values)	Maximum (100 values)	Minimum (100 values)
	Maximum of Likelihood	50.033%	94.858%	421.851%	4.036%
	Maximum A Posteriori	0.021%	0.470%	1.615%	0.086%

Given in percentage of the reference signal.


Huimei Han

Jun Zhao
email: 2junzhao@ntu.edu.sg

Dusit Niyato
email: dniyato@ntu.edu.sg

Marco Di Renzo
email: marco.di.renzo@gmail.com

Quoc-Viet Pham

Intelligent Reflecting Surface Aided Network: Power Control for Physical-Layer Broadcasting

Keywords: Intelligent reflecting surface, wireless communication, power control, quality of service

As a recently proposed idea for future wireless systems, intelligent reflecting surface (IRS) can assist communications between entities which do not have high-quality direct channels in between. Specifically, an IRS comprises many low-cost passive elements, each of which reflects the incident signal by incurring a phase change so that the reflected signals add coherently at the receiver. In this paper, for an IRS-aided wireless network, we study the problem of power control at the base station (BS) for physical-layer broadcasting under quality of service (QoS) constraints at mobile users, by jointly designing the transmit beamforming at the BS and the phase shifts of the IRS units. Furthermore, we derive a lower bound of the minimum transmit power at the BS to present the performance bound for optimization methods. Simulation results show that, the transmit power at the BS approaches the lower bound with the increase of the number of IRS units, and is much lower than that of the communication system without IRS.

I. INTRODUCTION

Benefiting from various advanced technologies, the fifthgeneration (5G) communications achieves great improvements in spectral efficiency, such as massive antennas deployment at the base station (BS) (i.e., massive multiple-input multipleoutput), non-orthogonal multiple access, millimeter wave communications, and ultra-dense hetnets. However, these advanced technologies will introduce great mounts of energy consumption, resulting in the high complexity of hardware implementation, which brings great challenges for practical implementations [1]. For example, the ultra-dense HetNets mean that there are lots of BSs in the network, and the energy consumption scales up with respect to the number of BSs. The massive antenna arrays consist of active elements and transmit/recieve data, thus consuming energy expensively.

Intelligent reflecting surface. To improve the spectral efficiency and reduce the energy consumption, researchers are exploring new ideas for future wireless systems [2]- [START_REF] Basar | Index modulation techniques for next-generation wireless networks[END_REF]. Among theses ideas, intelligent reflecting surface (IRS) has been considered in several studies [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF]- [START_REF] Nadeem | Asymptotic analysis of large intelligent surface assisted MIMO communication[END_REF]. An IRS is a planar array consisting of many reflecting and nearly passive units. Each IRS unit is controlled by the BS remotely to change the phase of the incident signal, so that the signals at the receiver can add coherently. In other words, IRS intelligently adjusts the propagation conditions to improve communication quality between the BS and mobile users (MUs). Since each IRS unit only reflects signals in a passive way, instead of transmitting/receiving signals in an active way, the energy consumption is very low. In addition, due to the characteristics of lightweight and low profile, the IRS can be deployed on walls/building facades, and the channel model between the BS and the IRS is usually characterized as line of sight (LoS) [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF].

Our contributions. To address the problem of power control at the BS for physical layer broadcasting under quality of service (QoS) constraints at the MUs in an IRS-aided network, we propose to employ the alternating optimization algorithm to jointly design the transmit beamforming at the BS and the IRS units. Furthermore, we derive the lower bound of the minimum transmit power for the broadcast setting to present the performance bound for optimization methods. Simulation results show that, for the broadcasting transmit pattern, the transmit power at the BS approaches the lower bound with the increase of the number of IRS units, and is much lower than that of the communication system without IRS.

Comparing this paper with [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF], [START_REF]Beamforming optimization for intelligent reflecting surface with discrete phase shifts[END_REF]. Recently, Wu and Zhang [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF], [START_REF]Beamforming optimization for intelligent reflecting surface with discrete phase shifts[END_REF] also considered downlink power control under QoS, with phase shifts of IRS units having continuous domains in [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF] and discrete domains in [START_REF]Beamforming optimization for intelligent reflecting surface with discrete phase shifts[END_REF]. The differences between our paper and [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF], [START_REF]Beamforming optimization for intelligent reflecting surface with discrete phase shifts[END_REF] are twofold. First, our paper considers the broadcast setting, while [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF], [START_REF]Beamforming optimization for intelligent reflecting surface with discrete phase shifts[END_REF] are for the unicast setting. Second, under the line of sight (LoS) channel model between the BS and IRS, we analyze a lower bound of the minimum transmit power in the general setting (i.e., an arbitrary number of antennas at the BS, an arbitrary number of IRS units, and an arbitrary number of MUs), while [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF], [START_REF]Beamforming optimization for intelligent reflecting surface with discrete phase shifts[END_REF] only derived the relationship between the transmit power and the received power in the special setting of considering single user and single-antenna BS, ignoring the channel between BS and MU, and assuming that the channel model between the BS and IRS is Rayleigh fading.

Other related work. For IRS-aided wireless communications, in addition to [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF], [START_REF]Beamforming optimization for intelligent reflecting surface with discrete phase shifts[END_REF] above and [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF]'s conference version [START_REF] Wu | Intelligent reflecting surface enhanced wireless network: Joint active and passive beamforming design[END_REF] for downlink power control under QoS, optimizing transmit power at BS is addressed in [START_REF] Nadeem | Asymptotic analysis of large intelligent surface assisted MIMO communication[END_REF] to maximize the minimum SINR among users and in [START_REF] Guo | Weighted sumrate optimization for intelligent reflecting surface enhanced wireless networks[END_REF], [START_REF] Pan | Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer[END_REF] to maximize the weighted sum of downlink rates. An earlier draft [START_REF] Zhao | Optimizations with intelligent reflecting surfaces (IRSs) in 6G wireless networks: Power control, quality of service, max-min fair beamforming for unicast, broadcast, and multicast with multi-antenna mobile users and multiple IRSs[END_REF]

... , b r H ,1 r h ,2 r h , r K h ,1 b h ,2 b h , b K h
Fig. 1: An IRS-aided communication system consisting of a base station (BS), multiple mobile users (MUs), and an IRS comprising many IRS units, where phase shifts incurred by the IRS units are remotely controlled by the BS.

under QoS in the unicast, multicast, and broadcast settings. As an updated version, this current paper adds simulation results and also derives a lower bound for the minimum transmit power at the BS. In the absence of IRS, downlink power control under QoS for the broadcast setting is studied in the seminal work [START_REF] Sidiropoulos | Transmit beamforming for physical-layer multicasting[END_REF].

Organization. The remainder of this paper is organized as follows. Section II presents the system model and formulates the problem of power control under QoS. In Section III, we describe an algorithm to solve the problem. A lower bound for the minimum transmit power is elaborated in Section IV. Sections V and VI give numerical results and the conclusion respectively.

Notation. We utilize italic letters, boldface lower-case and upper-case letters to denote the scalars, vectors and matrices respectively. (•) T and (•) H stand for the transpose and conjugate transpose of a matrix, respectively. We utilize D i,j and x i to stand for the element in the i th row and j th column of D and the i th element of x respectively. C denotes the set of all complex numbers. I is the identity matrix. CN (µ, σ 2 ) denotes a circularly-symmetric complex Gaussian distribution with mean µ and variance σ 2 . Let • and | • | denote the Euclidean norm of a vector and cardinality of a set respectively. diag(x) means a diagonal matrix with the element in the i th row and i th column being the i th element in x. arg(x) stands for the phase vector. E(•) and Var(•) are the expectation and variance operations, respectively. For a square matrix M , we use M -1 and M 0 to denote its inverse and positive semidefiniteness. respectively.

II. SYSTEM MODEL AND PROBLEM DEFINITION

A. System model

We consider IRS-aided communications in the broadcast setting, where there are a BS with M antennas and an IRS with N IRS units, and K single-antenna MUs, as shown in Fig. 1. We consider that the BS utilizes linear transmit precoding as the beamforming vector, denoted by w ∈ C M×1 , and thus the transmitted signal at the BS is x = ws where s is the broadcasted data. When BS broadcasts the signal x, it will arrive at each MU via indirect and direct channels, and the received signal at each MU is the superposed signal from the two channels. More specifically, for the indirect channel, the transmitted signal x travels from the BS to the IRS, reflected by the IRS, and finally travels from the IRS to these K MUs. For the direct channel, the transmitted signal x travels from BS to these K MUs directly.

Let Φ = diag(β 1 e jθ1 , . . . , β N e jθN ) denote the reflection coefficient matrix at the IRS, where β n and θ n denote the amplitude factor and phase shift respectively. In this paper, we assume that the IRS only changes the phase of the reflected signal, i.e., θ n ∈ [0, 2π) and

β n = 1. Let H b,r ∈ C N ×M , h H r,i ∈ C 1×N
, and h H b,i ∈ C 1×M be the BS-IRS channel, IRSi th MU channel, and BS-i th MU channel respectively. Then, the received signal at MU i is given by

y i = (h H r,i ΦH b,r + h H b,i )ws + n i , i = 1, . . . , K, (1
) where n i ∼ CN (0, σ 2 i ) denotes the additive white Gaussian noise at MU i.

We assume that the broadcasted data s is normalized to unit power. Then, the signal-to-interference-plus-noise ratio (SINR) at MU i can be written as

SINR i = |h H i (Φ)w| 2 σ 2 i , (2)
where

h H i (Φ) = h H r,i ΦH b,r + h H b,i means the overall downlink channel from the BS to MU i.

B. Problem definition

The problem of power control under QoS for broadcasting, is to minimize the transmitted power at the BS under QoS. Note that the transmitted power at the BS is w 2 , and that the QoS of MU i is usually characterized by its SINR. Then, this problem can be formulated as

(P1): min w,Φ w 2 (3a) s.t. |h H i (Φ)w| 2 σ 2 i ≥ γ i , ∀i == 1, . . . , K, (3b)
0 ≤ θ n < 2π, n = 1, . . . , N, (3c)
where γ i is the SINR target. Furthermore, without loss of generality, we assume that all MUs have the same SINR target and the same noise variance, i.e.,

γ k = γ, σ 2 k = σ 2 .

III. ALTERNATING OPTIMIZATION ALGORITHM

In this section, we utilize alternating optimization, which is used for multivariate optimization in an alternating manner, to solve Problem (P1), as described in Algorithm 1. More specifically, we first optimize w given Φ, and then optimize Φ given w, which is performed iteratively to obtain the desired w and Φ. In the following, we describe the details of the j th iteration to illustrate the alternating optimization algorithm.

Optimizing w given Φ (j-1) . Given Φ (j-1) obtained during the (j -1) th iteration, Problem (P1) becomes the conventional

s.t. |h H i (Φ (j-1) )w| 2 σ 2 ≥ γ, i = 1, . . . , K, (4b)
Note that Problem (P2) is non-convex because of the non-convex constraint, which can be solved by a relaxation of Problem (P2) based on semi-definite programming (SDP) [START_REF] Cramér | Random Variables and Probability Distributions[END_REF], (P3): min X trace(X) (5a)

s.t. trace(XH i (Φ j-1 )) ≥ γσ 2 , i = 1, . . . , K, (5b) X 0, (5c)
where X and H i (Φ (j-1) ) are defined as X := ww H and H i (Φ (j-1) ) := h i (Φ (j-1) )h i (Φ (j-1) ) H respectively.

Apparently, Problem (P3) is an SDP, and we can utilize the convex optimization solvers (e.g., CVX [START_REF] Grant | CVX: MATLAB software for disciplined convex programming[END_REF]) to solve this problem. After X is available, the Gaussian randomization [START_REF] So | On approximating complex quadratic optimization problems via semidefinite programming relaxations[END_REF] is applied to obtain solution to Problem (P2). Note that, when utilizing the Gaussian randomization, we can obtain many candidate solutions to Problem (P2), and we select the one with the minimum power as the value of w during the j th iteration, denoted by w (j) .

Finding Φ given w (j) . Given w (j) , Problem (P1) becomes the following feasibility check problem of finding Φ:

Algorithm 1 Alternating optimization to find w and Φ for Problem (P1).

1: Initialize Φ as Φ (0) := diag(e jθ (0) 1 , . . . , e jθ (0) N ), where θ Given Φ as Φ (j-1) , solve Problem (P3) to obtain w (j) ; 5:

Compute the object function value P (j) t ← w (j) 2 ; 6:

if 1 -

P (j) t P (j-1) t ≤ ε then 7:
break; {Comment: ε controls the number of executed iterations before termination. The algorithm terminates if the relative difference between the transmit power obtained during the j th iteration and the (j -1) th iteration is no greater than ε.}

8:

end if {Comment: Finding Φ given w:}

9:

Given w as w (j) , solve Problem (P6) to obtain Φ (j) ; 10:

if Problem (P6) is then 11:

break;

12:

end if 13: end while

(P4) : Find Φ (6a) s.t. |h H i (Φ)w (j) | 2 σ 2 ≥ γ, (6b) 0 ≤ θ n < 2π, n = 1, ..., N. (6c)
Let φ = [e jθ1 , . . . , e jθN ] H , a i = diag(h H r,i )H b,r w (j) , and b i = h H b,i w (j) . Then, the (P4) can be rewritten as (P5) : Find φ (7a)

s.t. φ H , 1 A i φ 1 + b i b H i ≥ γσ 2 . ( 7b
)
|φ n | = 1, n = 1, . . . , N. (7c)
where

A i = a i a H i , a i b H i b i a H i , 0 .
Note that, since the constraints (7c) are non-convex, Problem (P5) is a non-convex optimization problem. However, by introducing an auxiliary variable t satisfying |t| = 1, Problem (P5) can be converted to a homogeneous quadratically constrained quadratic program (QCQP). Specifically, define

v := t φ 1 = φt t , and
V := vv H . (8)
Then, a relaxation of Problem (P5) based on SDP is

(P6) : max V ,α K i=1 α i (9a) s.t. trace(A i V ) + b i b H i ≥ α i + γσ 2 , (9b) V n,n = 1, n = 1, . . . , N + 1 (9c) V 0, α i ≥ 0, i = 1, . . . , K, (9d)
where variable α i can be described as MU i's "SINR residual" in the phase shift optimization [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF].

Similar to Problem (P2), we can utilize the convex optimization solvers to solve problem (P6). After V is available, the Gaussian randomization is applied to obtain many candidate solutions to (P4), denoted by [Φ (j) 1 , . . . , Φ

c ] where c is the number of candidate solutions. The rule of selecting one as the value of Φ during the j th iteration, denoted by Φ (j) , is described as follows.

First, we define

f := min(|h H 1 (Φ)w| 2 , . . . , |h H K (Φ)w| 2 ), (10)
where w = w w denotes the transmit beamforming direction at the BS. Replacing Φ and w in Eq. ( 10) with Φ (j-1) and w (j) = w (j) w (j) respectively, we can obtain the value of f after optimizing w given Φ (j-1) , denoted by f (j) ow (the subscript "o" represents optimization).

Next, after replacing Φ and w in Eq. ( 10) with Φ (j) k (k = 1, . . . , c) and w (j) respectively, we can obtain the value of f corresponding to Φ

(j) k , denoted by f (j) k . If f (j) k satisfies f (j) k ≥ f (j)
ow , we incorporate it into a set G, and select the Φ corresponding to the maximum element in G as Φ (j) . We denote the maximum element in G as f of f after optimizing Φ given w (j) .

Proposition 1. The rule of selecting one as the value of Φ during the j th iteration, ensures the objective value in Problem (P2) is non-increasing over the iterations.

Proof: Let P t = w 2 denote the transmit power. Given Φ, Problem (P2) can be rewritten as

min w P t s.t. P t |h H i (Φ)w| 2 σ 2 ≥ γ ∀i ∈ {1, 2, . . . , K}. (11)
Apparently, the minimum value of P t is

P t = γσ 2 min(|h H 1 (Φ)w| 2 , . . . , |h H K (Φ)w| 2 ) = γσ 2 f . ( 12
)
Note that, if the value of f after optimizing w given Φ is non-decreasing over the iterations, then P t is non-increasing over the iterations; i.e., if f

(j+1) ow ≥ f (j)
ow , then we have

P (j+1) t ≤ P (j)
t . Based on the rule of selecting one as the value of Φ during the j th iteration, it is easily to derive f

(j) oΦ ≥ f (j) ow .
Then, if the w (j+1) is the optimal solution to Problem (P2) during the (j + 1) th iteration, we derive f

(j+1) ow ≥ f (j) oΦ ≥ f (j) ow .
Hence, we have P (j+1) t ≤ P (j) t , which means that P t is non-increasing over the iterations.

IV. A LOWER BOUND FOR MINIMUM TRANSMIT POWER

In this section, for the IRS-aided broadcast pattern, we derive a lower bound of the minimum transmit power.

We assume that the BS-MUs and IRS-MUs channel are Rayleigh fading, and that BS-IRS channel is LoS. We consider the uncorrelated Rayleigh fading channel model for IRSi th MU and BS-i th MU; i.e. h r,i ∼ CN (0, β 2 r,i I), h b,i ∼ CN (0, β 2 b,i I), where β 2 r,i and β 2 b,i account for the path loss of IRS-MUs and BS-MUs respectively. Let (x BS , y BS , z BS ) and (x IRS , y IRS , z IRS ) be the coordinate of BS and IRS respectively. Then, the channel between the BS and IRS is given by [START_REF] Nadeem | Asymptotic analysis of large intelligent surface assisted MIMO communication[END_REF]

H b,r = β h 2 2 sg T , (13)
where

s = [s 1 , . . . , s m , . . . , s M ] T , b = [g 1 , . . . , g n , . . . , g N ] T s m = exp j 2π λ d BS (m -1)sinφ LoS1 sinθ LoS1 , m = 1,. . . ,M g n = exp j 2π λ d IRS (n -1)sinφ LoS2 sinθ LoS2 , n = 1,. . . ,N θ LoS1 = tan -1 d BS-IRS z IRS -z BS , θ LoS2 = π -θ LoS1 φ LoS1 = π -tan -1 y IRS -y BS x IRS -x BS , φ LoS2 = π + φ LoS1 , (14)
where λ is wavelength, d BS and d IRS are the inter-antenna separation at the BS and IRS respectively, φ LoS1 and φ LoS2 are the LoS azimuth at BS and IRS respectively, θ LoS1 and θ LoS2 denote the elevation angle of departure at BS and elevation angle of arrival at IRS respectively, β 2 h accounts for the path loss of IRS-BS, and d BS-IRS represent the distance between the BS and IRS.

Next, we present the details of deriving the lower bound of transmit power P t with respect to the number of IRS units N , the number of MUs K, and the number of antennas M , considering the following two cases of parameter settings: 1) K = 1 and M > 1; 2) K > 1 and M > 1. In addition, when discussing the case of K = 1, we omit the subscript i of β b,i and β r,i for presentation simplicity.

Case 1): K = 1 and M > 1. Based on Eq. ( 11), given Φ, the minimum transmit power P t is P t = σ 2 γ |h H 1 (Φ)w| 2 . Furthermore, because h H 1 (Φ) is random variance, transmit power P t should be considered to be the average transmit power, which is more accurately written as

P t = σ 2 γ E(|h H 1 (Φ)w| 2 ) . ( 15
)
This means that minimizing the transmit power is equivalent to maximizing the term

E(|h H 1 (Φ)w| 2 ). For |h H 1 (Φ)w|, we have |h H 1 (Φ)w| = |h H r,1 ΦH b,r w + h H b,1 w| (a) ≤ |h H r,1 ΦH b,r w| + |h H b,1 w|. ( 16
)
Based on the triangle inequality, Eq. 16(a) holds if and only

if arg(h H r,1 ΦH b,r w) = arg(h H b,1 w) = ϕ 0 . Let A = |h H r,1 ΦH b,r w|, B = |h H b,1 w|.
Then, the maximum value of E(|h H 1 (Φ)w| 2 ) with respect to Φ and w, denoted by Q 1 , is given by

Q 1 = max(E(|h 1 H (Φ)w| 2 )) = E((A + B) 2 ) = E(A 2 ) + 2E(AB) + E(B 2 ). (17)
Next, we discuss how to derive each term in Eq. [START_REF] Grant | CVX: MATLAB software for disciplined convex programming[END_REF]. For E(A 2 ), we have

                                       E(A) (a) = E( N n=1 |h H r,1,n || M m=1 H b,r,m,n w m |) = E( N n=1 |h H r,1,n ||C n |) (b) = E(|h H r,1,1 |)(|C 1 | + |C 2 | + . . . + |C N |), E 2 (A) = (|C 1 | + |C 2 | + . . . + |C N |) 2 E 2 (|h H r,1,1 |) (c) = (N |C 1 |) 2 β 2 r π 4 (d) ≤ πN 2 β 2 h β 2 r M 8 , Var(A) = Var( N n=1 |h H r,1,n ||C n |) (e) ≤ β 2 r 2 (2 - π 2 ) × β 2 h M 2 × N, E(A 2 ) = E 2 (A) + Var(A) ≤ πN 2 β 2 h β 2 r M 8 + N β 2 r β 2 h M 4 (2 -π 2 ), (18)
β 2 r 2 (2 -π 2 ).
This paper appears in the Proceedings of IEEE International Conference on Communications (ICC) 2020.

Please see https://arxiv.org/abs/1912.03468 for the full journal version. Feel free to contact us for questions/remarks.

P t ≥ P L K=1,M>1 = σ 2 γ max(E(|h 1 H (Φ)ω| 2 )) = σ 2 γ Q 1 = σ 2 γ πN 2 β 2 b β 2 r M 8 + N β 2 r β 2 h M 4 (2 - π 2 ) + N πβ r β h β b M 2 √ 2 + β 2 b 2 (2 - π 2 ) + πβ 2 b M 4 . (22)
For E(B 2 ), we have For E(AB), we have

                               E(B) (a) = E( M m=1 |h b,m H w m |) = E( M m=1 |h b,m H ||w m |) = (|w 1 | + |w 2 | + . . . + |w M |)E(|h H b,1 |), E 2 (B) = (|w 1 | + |w 2 | + . . . + |w M |) 2 E 2 (|h H b,1 |) (b) ≤ (M |w 1 |) 2 β 2 b π 4 = πβ 2 b M 4 , Var(B) = Var( M m=1 |h b,m H ||w m |) (c) ≤ β 2 b 2 (2 - π 2 ), E(B 2 ) = E 2 (B) + Var(B) ≤ πβ 2 b M 4 + β 2 b 2 (2 - π 2 ), (19)
E(AB) = E 2 (A)E 2 (B) = N πβ r β h β b M 4 √ 2 . ( 20
)
Substituting Eq. ( 18)-Eq. ( 20) into Eq. ( 17), we have

Q 1 = E(A 2 ) + 2E(AB) + E(B 2 ) = πN 2 β 2 b β 2 r M 8 + N β 2 r β 2 h M 4 (2 - π 2 ) + N πβ r β h β b M 2 √ 2 + β 2 b 2 (2 - π 2 ) + πβ 2 b M 4 .
(21) Then, substituting Eq. (21) into Eq. ( 15), the lower bound of the minimum transmit power at the BS in the case of K = 1, M > 1 is obtained in Eq. ( 22).

Case 2): K > 1 and M > 1. The minimum value of P t satisfying the constrains in Eq. ( 11), is

P t = γσ 2 min(E(|h H 1 (Φ)w| 2 ), . . . , E(|h H K (Φ)w| 2 )) . ( 23
)
Based on min(E(|h H 1 (Φ)w| 2 ), . . . , E(|h H K (Φ)w| 2 )) ≤ min(Q 1 , Q 2 , . . . , Q K ), the lower bound of the minimum transmit power at the BS is given by

P t ≥ P L K>1,M>1 = rσ 2 min(Q 1 , Q 2 , . . . , Q K ) , (24)
where Q i = max(E(|h H i (Φ)w| 2 )). Eq. ( 21) only presents how to get the value of Q 1 , and we can use the same way to compute the other values of Q i (i = 1, . . . , K).

In addition, since Problem P(3) has K linear constrains and M 2 variables, the complexity of solving Problem P [START_REF] Hu | Beyond massive MIMO: The potential of data transmission with large intelligent surfaces[END_REF] is O((K + M 2 ) 3.5 ) for one iteration [START_REF] Sidiropoulos | Transmit beamforming for physical-layer multicasting[END_REF]. Similarly, the complexity of solving Problem P( 6) is O(2K + (N + 1) 2 ) 3.5 ) for one iteration. Hence, the complexity of the proposed alternating optimization is O((K +M 2 ) 3.5 )+O(2K +(N +1) 2 ) 3.5 ) for one iteration. A future direction for us is to reduce the computational complexity and one potential idea is to use manifold optimization [START_REF] Pan | Intelligent reflecting surface for multicell mimo communications[END_REF].

V. SIMULATION RESULTS

In this section, we utilize numerical results to validate the derived lower bound of the transmit power and the alternating optimization algorithm. We assume that the BS with uniform linear array of antennas is located at (0,0,0), and that the IRS with uniform linear array of IRS units is located at (0,50,0). The inter-antenna and inter-unit separation at BS and IRS are half wavelength. The purpose of deploying IRS is to improve the signal strength. To illustrate this benefiting, we assume that the MUs are uniformly located at the half circle centered at the IRS with radius 2 m as shown in Fig. 2 Figs. 3 and4 show the variance of transmit power at BS with the number of IRS units for M = 20, 30, 40 and K = 2 respectively. We can see that, the transmit power decreases with the increase of the number of IRS units and the number of antennas at the BS, which significantly lower than the baselines. This indicates that, deploying IRS can actually improve the signal strength and thus decrease the transmit power at the BS. Furthermore, the results from Figs. 3 and 4 also show that the transmit power at the BS approaches the lower bound with the increase of the the number of IRS units, which coincide with the analysis results. Notice that Fig. 3 is more obvious than Fig. 4, and actually the speed of approaching the lower bound in Fig. 4 is extremely slow. 50. The results show that, with the increase of the number of MUs, the transmit power increases, dramatically lower than the baselines, and the gap between the transmit power and the lower bound widens up. To obtaining a better bound which grows with K is our future direction.

VI. CONCLUSION

In this paper, we have proposed a solution to the power control under QoS for an IRS-aided wireless network. Specifically, we utilize the alternating optimization algorithm to jointly optimize the transmit beamforming at the BS and the passive IRS units at the IRS. Furthermore, we derived a lower bound of the minimum transmit power for the IRS-enhanced physical layer broadcasting. Simulation results show that, the transmit power at the BS approaches the lower bound with the number of IRS units, and is significantly lower than that of the communication system without IRS. Transmit power at the BS (dB)

  1, 2, . . . , N ) is chosen uniformly at random from [0, 2π); 2: Initialize the iteration number j ← 1; 3: while 1 do {Comment: Optimizing w given Φ:} 4:

  (j) oΦ , which is the value This paper appears in the Proceedings of IEEE International Conference on Communications (ICC) 2020. Please see https://arxiv.org/abs/1912.03468 for the full journal version. Feel free to contact us for questions/remarks.

  where C n = M m=1 H b,r,m,n w m , step (a) follows from the fact that arg(h H r,1 ΦH b,r w) = ϕ 0 , step (b) follows from the fact that h r,1 ∼ CN (0, β 2 r,1 I). For step (c), Since the element in H b,r,m,n has the same amplitude and [w 1 , . . . , w M ] is the normalized vector, it is easy to derive that |C 1 | = . . . = |C N |. Step (c) also follows from the fact that |h r,1,1 H | has distribution of Rayleigh with mean βr √ π 2 , and steps (d)(e) follow from the fact that term |C 1 | 2 ≤ Mβ 2 h 2 and that |h r,1,1 H | has distribution of Rayleigh with variance

  where step (a) follows from the fact that arg(h H b,1 w) = ϕ 0 , steps (b) and (c) follow from the fact that |h b,1 H | has distribution of Rayleigh with mean β b π 2 ), step (b) also follows from the fact that term (|w 1 | + |w 2 | + . . . + |w M |) 2 takes the maximum value if |w 1 | = |w 2 | = . . . = |w M | because of M m=1 |w m | 2 = 1.

Fig. 2 :

 2 Fig. 2: The location of IRS, BS and MUs in the simulation.Figs.3 and 4show the variance of transmit power at BS with the number of IRS units for M = 20, 30, 40 and K = 2 respectively. We can see that, the transmit power decreases with the increase of the number of IRS units and the number of antennas at the BS, which significantly lower than the baselines. This indicates that, deploying IRS can actually improve the signal strength and thus decrease the transmit power at the BS. Furthermore, the results from Figs.3 and4 also show that the transmit power at the BS approaches the lower bound with the increase of the the number of IRS units, which coincide with the analysis results. Notice that Fig.3is more obvious than Fig.4, and actually the speed of approaching the lower bound in Fig.4is extremely slow. Fig.5show the variance of transmit power at BS with the number of MUs K ranging from 1 to 10 for M = 40, N =

Fig. 5 Fig. 3 :Fig. 4 :

 534 Fig. 2: The location of IRS, BS and MUs in the simulation.Figs.3 and 4show the variance of transmit power at BS with the number of IRS units for M = 20, 30, 40 and K = 2 respectively. We can see that, the transmit power decreases with the increase of the number of IRS units and the number of antennas at the BS, which significantly lower than the baselines. This indicates that, deploying IRS can actually improve the signal strength and thus decrease the transmit power at the BS. Furthermore, the results from Figs.3 and4 also show that the transmit power at the BS approaches the lower bound with the increase of the the number of IRS units, which coincide with the analysis results. Notice that Fig.3is more obvious than Fig.4, and actually the speed of approaching the lower bound in Fig.4is extremely slow. Fig.5show the variance of transmit power at BS with the number of MUs K ranging from 1 to 10 for M = 40, N =

Fig. 5 :

 5 Fig. 5: Transmit power vs. the number of MUs.

  of our current paper summarizes problems of downlink power control This paper appears in the Proceedings of IEEE International Conference on Communications (ICC) 2020. Please see https://arxiv.org/abs/1912.03468 for the full journal version. Feel free to contact us for questions/remarks.

		IRS
	BS
		MU K
	MU1	MU2

  This paper appears in the Proceedings of IEEE International Conference on Communications (ICC) 2020. Please see https://arxiv.org/abs/1912.03468 for the full journal version. Feel free to contact us for questions/remarks.

	power control problem under QoS in the downlink broadcast
	channel without an IRS.
	(P2): min w	w 2	(4a)

  [1] S. Zhang, Q.Wu, S. Xu, and G. Y. Li, "Fundamental green tradeoffs: Progresses, challenges, and impacts on 5G networks," IEEE Communications Surveys Tutorials, vol. 19, no. 1, pp. 33-56, Firstquarter 2017. [2] S. Buzzi, C. I, T. E. Klein, H. V. Poor, C. Yang, and A. Zappone, "A survey of energy-efficient techniques for 5G networks and challenges ahead," IEEE Journal on Selected Areas in Communications, vol. 34, no. 4, pp. 697-709, April 2016.

	1	2	3	4	5	6	7	8	9	10
					Number of MUs K