Bibek Kc

Hassan Alkhwaildi
email: hassan.alkhwaildi@siu.edu

Multi-Objective TLBO and GWO-based Optimization for Placement of Renewable Energy Resources in Distribution System

Keywords: Distributed Generation, Fuzzy Decision, Hybrid Teaching Learning -Gray Wolf Algorithm, Multi-Objective Optimization, Renewable Energy Sources

The use of renewable solar and wind resources as distributed generation sources in distribution networks has been welcomed by network operators. In order to exploit the maximum benefits of using these distributed products, the location of installation and their capacity should be determined optimally in the distribution network. In this paper, in order to optimize the placement of solar panels and wind turbines in the distribution network with the aim of reducing losses and improving reliability based on Energy Not Supplied subscribers (ENS), a multi-objective evolutionary algorithm based on fuzzy decision method, called the Multi-Objective Hybrid Training Learning Based Optimization-Grey Wolf Optimizer (MOHTLBOGWO) proposed that has a High optimization speed and not trapped at all in the optimal local. At first, the candidate buses are set for the installation of renewable resources using the Loss Sensitivity Factor (LSF). Then the proposed method is used to determine the location and optimal capacity of renewable resources through the candidate bases. Proposed issues have been implemented in a single-objective and multiobjective manner on a 33 bus IEEE radial distribution network. Also, in this paper, the effect of distributing renewable resources on the characteristics of the distribution network is evaluated. The results obtained from the proposed algorithm are compared with the results of other algorithms to demonstrate the superiority of the proposed method in reducing losses, improving reliability, and increasing the financial profit of the network. Simulation results show the better performance of the proposed method in comparison with Teaching-Learning Based Optimization (TLBO) and Grey Wolf Optimiser (GWO) methods and past studies to achieve optimal results. Also, the results show that distributing of the capacity and location of distributed renewable generation leads to a further reduction in losses and a better improvement of the reliability criterion.

Introduction

The growing consumption of electric energy, mainly produced by burning fossil fuels, leads to various issues, such as environmental and financial issues. On the other hand, nearly 15% of active power is produced in large power plants [START_REF] Holttinen | Design and operation of power systems with large amounts of wind power[END_REF][START_REF] Coello-Coello | An Introduction to Evolutionary Multi-Objective Optimization with Some Applications in Pattern Recognition[END_REF]. This power is lost in transmission lines and distribution networks. Reducing power losses is one of the most interesting and most important issues in power system studies. There are several ways to minimize losses in power systems at the level of transmission and distribution networks [START_REF] Hayerikhiyavi | Comprehensive Analysis of Continuously Variable Series Reactor Using GC Framework[END_REF]. DG placement is one of these methods. These methods have advantages such as improving reliability, improving the voltage profile, and reducing the power loss of the active. However, inappropriate placement and inappropriate DG planning can weaken the distribution network characteristics [START_REF] Mehta | Optimal selection of distributed generating units and its placement for voltage stability enhancement and energy loss minimization[END_REF][START_REF] Daud | A review: Optimal distributed generation planning and power quality issues[END_REF]. Several studies have been conducted in the field of exploitation of distribution networks by optimal utilizing the dispersed generation sources in the distribution network. Classical, analytical, and meta-exploratory methods are the methods that have been proposed in this field. Classical methods are methods that can determine optimal response over a short period of time, but with increasing problem dimensions they cannot determine the optimal response. Analysis methods are methods that require additional computations and may not be able to achieve absolute optimal points [START_REF] Jenkinson | Improved optimization for the robust and accurate linear registration and motion correction of brain images[END_REF]. Recently; Metaheuristic methods have been used to solve the problems of distribution network utilization. In [START_REF] Zamzam | Data-driven learning-based optimization for distribution system state estimation[END_REF], the Teachinglearning algorithm has been used and optimization has been achieved as a single-objective or two-objective. Placement of wind turbines and solar cells in the distribution network with the goal of reducing losses and improving the voltage stability conducted as a two-objective issue using the Particle Swarm Optimization algorithm [START_REF] Golestani | Distributed generation dispatch optimization by artificial neural network trained by particle swarm optimization algorithm[END_REF][START_REF] Behbahani | Safety Assessment on Merge and Diverge Areas using Fuzzy Inference System, Artificial Neural Network, and Particle Swarm Optimization[END_REF]. In [START_REF] Prabha | Optimal placement and sizing of multiple distributed generating units in distribution networks by invasive weed optimization algorithm[END_REF][START_REF] Yazdani | Distribution Simultaneous Sizing, Sitting, and Service Area Determination of Sub-transmission Substations and Distributed Generations Considering Load Uncertainty by GSO Algorithm[END_REF], in order to reduce the true power losses in distribution networks, the optimal placement of dispersed generation resources have been used in a combination and single-objective approach. In this study, a sensitivity analysis method was used to determine the best mode for locating dispersed generation sources and capacitors in the network. In [START_REF] Khatod | Evolutionary programming based optimal placement of renewable distributed generators[END_REF], the optimal placement of renewable distributed generation sources is presented using evolutionary planning and its results are compared with the genetic method. The results show that, compared to the proposed method, the genetic algorithm finds the optimal problem at lower speeds. In [START_REF] Sattianadan | Optimal placement of DG in distribution system using genetic algorithm[END_REF], the locating of dispersed generation in distribution networks is provided by considering the voltage stability criteria. In this study, the genetic algorithm is used to find the best answer. In [START_REF] Kaur | Effective Loss Minimization and Allocation of Unbalanced Distribution Network[END_REF][START_REF] Rao | Teaching-learning-based optimization algorithm[END_REF], a new and high-efficiency algorithm is proposed using a firefly algorithm to determine the location and capacity of dispersed generation sources in an unbalanced distribution network in a single-objective optimization. The goal of the problem is to reduce network losses.

The use of Genetic Algorithm (GA) and Simulated Annealing algorithm for optimal allocation of DG in distribution network [START_REF] Gandomkar | A combination of genetic algorithm and simulated annealing for optimal DG allocation in distribution networks[END_REF][START_REF] Tavana | An integrated intuitionistic fuzzy AHP and SWOT method for outsourcing reverse logistics[END_REF] are presented in the form of single-objective optimization with the aim of reducing of losses. In [START_REF] Moradi | A combination of genetic algorithm and particle swarm optimization for optimal DG location and sizing in distribution systems[END_REF][START_REF] Shirai | Configuration of fuzzy membership function circuits using current mode[END_REF], the installation location of the DG is determined by the GA algorithm and its capacity is determined by the Particle Swarm Optimization (PSO). This work has been done on the basis of the weight coefficient method and is presented with the aim of reducing power losses, decreasing voltage oscillation of the bus and improving the voltage stability of the distribution network. In [START_REF] Falaghi | ACO based algorithm for distributed generation sources allocation and sizing in distribution systems[END_REF][START_REF] Shakarami | A new approach for network reconfiguration problem in order to deviation bus voltage minimization with regard to probabilistic load model and DGs[END_REF], the Ant Colony Optimization (ACO) algorithm is proposed to determine the best installation location of the DG in the distribution network. This algorithm is presented in the form of single-objective optimization for reducing losses. Harmony Search algorithm (HS) is presented in [START_REF] Rao | Optimal Placement of Multiple Distributed Generator by Hs Algorithm[END_REF] as an optimization method for solving this problem with the goal of reducing losses based on single-objective optimization. The Artificial Colony Bee (ABC) algorithm has been introduced in [START_REF] Padma Lalitha | Application of fuzzy and ABC algorithm for DG placement for minimum loss in radial distribution system[END_REF][START_REF] Vosoogha | An Intelligent Power Prediction Method for Wind Energy Generation Based on Optimized Fuzzy System[END_REF] to determine the capacity, optimal power coefficient and DG installation location. This algorithm is presented in the form of single-objective optimization for reducing active losses. In [START_REF] Hadidian-Moghaddam | A multi-objective optimal sizing and siting of distributed generation using ant lion optimization technique[END_REF], the Differential Evolution Algorithm (DEA) has been used to optimize single-objective DGs and achieve minimum losses. In [START_REF] Devi | Application of modified bacterial foraging optimization algorithm for optimal placement and sizing of distributed generation[END_REF], the Bacterial Foraging Optimization (BFO) Search Algorithm was used to optimize the planning of DG units in the distribution system with the goal of reducing losses, reducing operating costs, and improving voltage stability. This work is based on the weighting coefficient method and the multi-objective optimization form. The Imperialistic Competitive Algorithm (ICA) is used in [START_REF] Soroudi | Imperialist competition algorithm for distributed generation connections[END_REF] as a singleobjective optimization to minimize losses in the system. The Plant Growth Simulation Algorithm (PGSA) has been used in [START_REF] Prabha | Determining the optimal location and sizing of distributed generation unit using Plant Growth Simulation Algorithm in a radial distribution network[END_REF][START_REF] Rajaram | Power system reconfiguration in a radial distribution network for reducing losses and to improve voltage profile using modified plant growth simulation algorithm with Distributed Generation (DG)[END_REF][START_REF] Khatibi | Comparative Analysis of Power System Model Reduction[END_REF] to determine the unit capacity and the loss sensitivity factor in choosing optimal DG location. Firefly Algorithm (FA) is used in [START_REF] Shakarami | A new approach for network reconfiguration problem in order to deviation bus voltage minimization with regard to probabilistic load model and DGs[END_REF][START_REF] Nadhir | Firefly algorithm for optimal allocation and sizing of Distributed Generation in radial distribution system for loss minimization[END_REF] to optimize singleobjective DGs with the goal of reducing power losses. In [START_REF] Ali | Ant Lion Optimization Algorithm for optimal location and sizing of renewable distributed generations[END_REF][START_REF] Reddy | Ant lion optimization algorithm for optimal sizing of renewable[END_REF], Ant Lion Optimizer (ALO) algorithm is used to determine optimal location of installation and DG size based on renewable resources with the goal of reducing losses and improving the profiles and voltage stability in a multiobjective optimization based on weight coefficient method. In [START_REF] Hasanien | A cuckoo search algorithm optimizer for steady-state analysis of self-excited induction generator[END_REF], the optimal placement of wind turbines in the distribution network is proposed using Cuckoo Search Algorithm (CSA) in a single-objective optimization with the goal of reducing power losses. In [START_REF] Choi | Multi-objective optimization of voltage-stability based on congestion management for integrating wind power into the electricity market[END_REF][START_REF] Shokouhandeh | Optimal Estimation of Capacity and Location of Wind, Solar and Fuel Cell Sources in Distribution Systems Considering Load Changes by Lightning Search Algorithm[END_REF], the locating of wind turbines and solar panels is proposed to reduce power losses and improve voltage stability using PSO algorithm based on multi-objective optimization and weight coefficient method.

As mentioned in most of the previous studies, the problem of positioning is presented as single or twoobjective optimization with the goal of reducing losses and improving the voltage profile, based on the weight coefficients method. In order to achieve a more accurate and realistic exploitation of the distribution network, the problem of the location of distributed generation should be considered as a multi-objective. One of the best ways for solving multiobjective problems is the use of Pareto levels to determine the optimal answer, which in this paper focuses on it and presents a multi-objective hybrid algorithm. On the other hand, the placement of renewable resources from the perspective of reliability, and considering the network failure rate and energy not-supplied network, has been less investigated. In this paper, these concepts are considered. In this paper, the proposed method is based on the combination of Teaching-learning and gray wolfs algorithms The Teaching-learning algorithm was presented in 2012 [START_REF] Rao | Teaching-learningbased optimization: an optimization method for continuous nonlinear large scale problems[END_REF] and the high convergence rate and the lack of control parameters were the features of this method. The gray wolf method is presented in 2014 [START_REF] Mirjalili | Grey wolf optimizer[END_REF] and has been greatly welcomed to solve the power engineering optimization problems. This method has features such as high convergence speed, high computational power, fast and large searches and Avoid Local Optimal. Therefore, in this paper, in order to locating solar panels and wind turbines in a 33-buses radial standard distribution network and reducing losses and improving reliability, Utilizing the benefits of both methods, a multiobjective hybrid algorithm based on a fuzzy decision method called MOTLBOGWO is proposed. The main contribution of this article are following:

• Use of a new cumulative smart method called the HTLBOGWO method • Solve multi-objective problem of locating renewable resources based on fuzzy decision method in MOHTLBOGWO algorithm • Study of the location of renewable energy sources from the perspective of reliability • Assess the dispersal of location and capacity of renewable energy sources over network characteristics The rest of this paper is organized as follows: In Section 2, Loss Sensitivity Factor Formulation is described for determining the candidate buses for the installation of renewable energy sources. In Section 3, the problem formulation including optimization target function and problem constraints is presented. In Section 4, the proposed algorithm, the fuzzy decision method and implementation of the proposed method in problem solving described. In Section 5, the simulation results are presented and the conclusions are presented in Section 6.

Loss Sensitivity Factor

In this study, the Loss-Sensitivity Factor (LSF) is presented to show buses that are sensitive to losses [START_REF] Reddy | Application of loss sensitivity factor and genetic algorithm for capacitor placement for minimum loss in radial distribution system[END_REF][START_REF] Asadi | Optimal placement and sizing of capacitor banks in harmonic polluted distribution network[END_REF]. Also, using this factor can reduce the search area and the time of the optimization process. In other words, the buses that are in trouble for losses are identified by this factor, and thus instead of searching all the network buses for the installation of dispersed generation by the optimization method, the search area is limited to these buses, thereby the time of the calculations is reduced. For a transmission line 'L' connected between the buses "i" and "k", as shown in Figure 1, the active power losses on this line are defined by RI2, defined as Eq. ( 1) [START_REF] Ali | Ant Lion Optimization Algorithm for optimal location and sizing of renewable distributed generations[END_REF][START_REF] Reddy | Application of loss sensitivity factor and genetic algorithm for capacitor placement for minimum loss in radial distribution system[END_REF].

(P Q )R P (V ) - + = (1)
The LSF can be calculated as follows [START_REF] Sardi | A loss sensitivity factor method for locating ES in a distribution system with PV units[END_REF]:

ik loss k ik 2 k k P 2Q R Q (V ) -  =  (2)
The normalized voltage is obtained by dividing the voltage in the base conditions by 0.95 [START_REF] Ali | Ant Lion Optimization Algorithm for optimal location and sizing of renewable distributed generations[END_REF][START_REF] Reddy | Application of loss sensitivity factor and genetic algorithm for capacitor placement for minimum loss in radial distribution system[END_REF]. Buses whose their voltages are less than 1.01 can be considered as the candidate bass for DG installation. It's worth noting that LSFs introduce a sequence of buses for the installation of DG.

Problem Formulation

In this study, multi-objective optimization of location and capacity of solar panels and wind turbines in the distribution network using the MOHTLBOGWO hybrid algorithm is proposed to reduce losses and improve reliability. In order to achieve Pareto set levels, Problem optimization variables, the optimal location and capacity of dispersed wind and solar generation sources are determined using intelligent singleobjective optimization (using the HTLBOGWO method) as well as multi-objective based on decision-making fuzzybased approach.

Objective Function

In this study, the minimization of power losses and improvement of Reliability is considered as an objective function and is described below.

• Line loss

Total network losses are equal to the losses of all network lines. To calculate network losses, you must first get the current of lines. To calculate the current of lines, it is also necessary to solve the load distribution problem and obtain the voltage of the buses [START_REF] Haggi | The effect of rescheduling power plants and optimal allocation of STATCOM in order to Improve power system static security using TLBO algorithm[END_REF]. By calculating the buses voltages, the line currents and network losses are calculated as follows:

(3) [START_REF] Mehta | Optimal selection of distributed generating units and its placement for voltage stability enhancement and energy loss minimization[END_REF] where k is between the buses i and j. Rk and Xk are the resistance and reactance of the k-th line, and Nb is the number of grid lines.

• Reliability

Reliability indicators are distribution system reliability assessment factors that were first introduced by the IEEE in 1998 to evaluate the distribution system. Basic reliability indicators are three load point indicators, a) The average rate of breakdown, b) the average withdrawal, and (c) the annual withdrawal period. A radial distribution system includes a set of series equipment including lines, brackets, fuses, and disconnected switches, and so on. A subsystem connected to the load point of a system requires that all equipment between it and the supply point be in operation. Therefore, the principle of serial systems can be considered directly for these systems as follows.

1

n p i i   = =  (5) 1 n p i i i ur  = =  (6) 1 1 n ii p i p n p i i r u r    = = ==   (7)
In the above equation, n is represent the number of outage events affecting load point p, i  denotes Equipment failure rate I (failure/year) and ri is repair time equip I (hours). Based on the basic indices of the load point and energy consumption at load points, energy not-supplied (ENS) (kWh/year) is calculated as follows [START_REF] Sabpayakom | Power losses reduction and reliability improvement in distribution system with Very Small Power Producers[END_REF].

() 1 nl a j j j ENS L u = =  (8)
where, nl represents the total number of load points, the unavailability of the load point j (hours/year), and the average load connected at the load point j (kW).

Constraints

In solving this optimization problem, a series of equal and unequal constraints is considered, which is expressed in sequence.

• Equilibrium power equations

The sum of the algebraic input and output power in the distribution system should be equal to:

ij k kk VV I R jX - = + 2 1 b N loss k k k P R I = =  DG DG N L N Swing DG Lineloss i 1 i 1 q 1 N L N Swing DG Lineloss i 1 i 1 q 1 P P (i) P (i) Pd(q) Q Q (i) Q (i) Qd(q) = = = = = = + = + + = +       (9)
• Voltage constraint

The value of the voltage in each bus should be between a minimum and a maximum value as follows:

min i max V V V  (10)
In the above relation, Vmin and Vmax are considered to be 0.95 and 1.05 p.u respectively [START_REF] Ali | Ant Lion Optimization Algorithm for optimal location and sizing of renewable distributed generations[END_REF].

• Constraints on the production of DG

To prevent reverse power flow, the installed DG capacity in the network is limited as follows, which should not exceed the power supplied by the post [START_REF] Ali | Ant Lion Optimization Algorithm for optimal location and sizing of renewable distributed generations[END_REF].

DG N L N DG Lineloss i 1 i 1 q 1 3 P (i) P (i) Pd(q) 4 = = =     +        (11) DG N L N DG Lineloss i 1 i 1 q 1 3 Q (i) Q (i) Qd(q) 4 = = =     +        (12)
On the other hand, DGs themselves have the minimum and the maximum production capacity is as follows [START_REF] Sabpayakom | Power losses reduction and reliability improvement in distribution system with Very Small Power Producers[END_REF].

min max DG DG DG P P (i) P   ( 13
) min max DG DG DG Q Q (i) Q   (14)
• Line capacity constraints

The Complex power transmission of each line should be less than the nominal value given in the following equation.

Li

Li(rated) SS 

Overview of Proposed Algorithm

TLBO

The TLBO algorithm is a smart optimization method that was introduced by Rao [START_REF] Rao | Teaching-learning-based optimization algorithm[END_REF] based on the influence of teacher to students to increase scientific level of class. This method is based on this principle that the teacher tries to close class level to himself and students, in addition to exploit the teacher's knowledge with regard to other classmates, use their knowledge to increase level of them. Because of the teacher can't bring level of individual students to himself, so tries to increase the average level of whole class and evaluates the class level based on the exams and students' scores. The mathematical expression of this approach is that first the population of problem variables (teacher and students) are defined randomly. All of these populations are compared together by the objective function and set of variables with best solution are considered as the teacher.

This approach is divided into two phases: teacher phase and student phase.

• Teacher phase

In this step teacher tries to bring class average to himself. But since it is very difficult, teacher tries to increase class average from Mi to M_new. Each set of problem variables are updated based on the difference of these two values. Difference of these two values can be saved by the parameter Diff_Mean as follows:

_ ( _ ) i i f i Diff Mean r M new T M =- ( 16
)
Where Tf is the teacher parameter that is selected randomly between 1 and 2. The ri is a random number between 0 and 1. Using the follow equation each set of variables are updated.

, ,

_ new i old i i X X Diff Mean =+ ( 17
)
• Student phase Students in addition to teacher's knowledge, benefit from each other's knowledge. The mathematical expression of this approach is that in each step and in each repetition each set of variable (student) selects one of students randomly. For example student i selects student j and this i is opposite of j. If the student j has more knowledge respects to student i then the student i updates his status based on the following equation: ,, ()

new i old i i i j X X r X X = + - (18)
The student status is varied as follows:

,, ()

new i old i i j i X X r X X = + - (19)
After the all students changed their status, their level is evaluated by the objective function. Under these conditions the best student is compared with the teacher of previous step and if a better result has, is replaced with previous iteration teacher. This process is continued to obtain convergence conditions.

GWO

One of the population-based intelligent and evolutionary algorithms is the Gray Wolf algorithm (GWO) which was first announced by Mir-Jalili in 2014 [START_REF] Mirjalili | Grey wolf optimizer[END_REF]. In this algorithm, the performance and behaviour of gray wolfs for hunting is simulated. Figure 2 shows that parameters such as α, β, δ, and ω represent the leaders of the group so that α directs the group as leader of the group and has important decisions about hunting, resting place and so on. The second group of leadership belongs to β. In addition to being able to help α due to having a good decision, β members are also the best substitutes for the α wolfs when they are old or dead; ω but is at the bottom of this group. They are the last wolves that are allowed to eat. The other group members that are not α, β and ω are called δ [START_REF] Mirjalili | Grey wolf optimizer[END_REF]. • Investigating, pursuing and following the hunt; • Pursuing, sieging and harassing the hunt until it stops; and • Attacking to the hunt.

In the simulation of the GWO, α is considered as a top answer. After that, the next two responses after α are considered β and δ. Finally, the rest answers are regarded ω.

To model a suitable situation that gray wolves encircle their hunt during hunting are shown as follows [START_REF] Mirjalili | Grey wolf optimizer[END_REF].

(20) [START_REF] Shakarami | A new approach for network reconfiguration problem in order to deviation bus voltage minimization with regard to probabilistic load model and DGs[END_REF] where t represents the repetition, A , C , and P X indicate the coefficients vectors, the position of the bait and the gray wolf respectively. The coefficients vector is obtained from the following equation [START_REF] Mirjalili | Grey wolf optimizer[END_REF]. According to the wolf ability due to finding and surrounding the hunt and in order to having a mathematical simulation of this kind of behaviour, it is assumed that the three top wolves in the group have more knowledge about the hunting place. By storing the position of the three above mentioned wolves, ω wolves have to change their location based on the position of the three top members of the group which can be expressed as follows [START_REF] Mirjalili | Grey wolf optimizer[END_REF]:

1 2 3 , , D C X X D C X X D C X X       = - = - = - ( 24
) 1 2 3 1 2 3 , , X X a D X X a D X X a D       = - = - = - (25) 1 2 3
( 1) 3

X X X Xt + + + = (26)

HTLBOGWO (TLRBO)

The TLBO algorithm is a robust algorithm capable of solving engineering problems. However, a new phase has been added to this algorithm to increase the global and local search capability as well as track the optimal answer. This phase is due to the choice of first, second and third students as the best answers after the student phase. After selecting the top three students, the remaining students will learn from these three students. The process of selection students and their ranking is considered to be the ranking phase, which is called the TLRBO algorithm.

Since the ranking phase and follow up of the students and their learning from the top three students are very similar to the behaviour of the gray wolf algorithm, So in the formulation of TLRBO, method of selection the wolf-alpha, beta, and delta and the method of updating the location of the wolfs of Omega to follow up the top three wolves have been used. Therefore, the proposed algorithm will be a combination of the TLBO and GWO algorithms.

HTLBOGWO (TLRBO) Multi-Objective Algorithm

Multi-objective issues have multiple goals, which are mostly contradictory. The answer to these problems is a set of solutions called Pareto optimal solutions [START_REF] Haggi | Multi-Objective PMU Allocation for Resilient Power System Monitoring[END_REF]. This set includes Pareto optimal solutions that represent the best balance between the objectives. Multi-objective optimization is considered as a minimization problem and formulated as follows [START_REF] Sabpayakom | Power losses reduction and reliability improvement in distribution system with Very Small Power Producers[END_REF]:

  1 2 0 min : F(x)= ( x), ( x),..., ( x), imize f f f ( 27
)
to: g (x) 0, i 1,2,...,

i subject m  = (28)
h (x) 0, i 1,2,...,

i p = = ( 29
) i i i L x U  (30)
where 0 represents the number of targets, m denotes the number of inequality constraints and p is the number of equality constraints. [Li, Ui] are the boundaries variable i. Given the nature of the multi-objective problems, various solutions cannot be compared through mathematical relations operators. In this case, Pareto optimal parental concepts allow us to compare two solutions in a multiobjective search space [START_REF] Sabpayakom | Power losses reduction and reliability improvement in distribution system with Very Small Power Producers[END_REF].

Two key points in finding a suitable set of optimal Pareto solutions for the given problem are convergence and coverage. Convergence refers to the ability of a multiobjective algorithm to determine the exact approximation of Pareto's optimal solution. The coverage reflects the distribution of Pareto's optimal response along the objectives. Since many of the current multi-objective algorithms are deductive, coverage and the number of answers for decision-making after the optimization process are very important. The ultimate goal for a multi-objective optimizer is to find the exact approximation of the optimal Pareto correct solution (convergence) with uniform distribution (coverage) for all purposes.

To solve multi-objective problems using the MOHTLBOGWO algorithm, it first equipped to an archive to store and restores the best approximation of the optimal Pareto solution. The position updating of the search factors for MOHTLBOGWO is the same, but the student's position is selected from the archive. In order to find a suitable expansion for the Pareto optimal frontier, a place from the low-density Pareto optimal frontier, similarly for multiobjective particle swarm optimization algorithm (MOPSO) [START_REF] Robles-Rodriguez | Multi-objective particle swarm optimization (MOPSO) of lipid accumulation in Fed-batch cultures[END_REF] is selected. To find the low-density area of Pareto optimal frontier, the search space should be divided. This is accomplished by finding the best and worst objectives of Pareto's optimal solution, the definition of the hyper-sphere to cover all the solution, and the division of the hyper-sphere into sub-hyper-sphere in each repetition. After the creation of the section, the choice is made by a roulette wheel mechanism with the following probability for each section proposed by Coello Coello [START_REF] Coello-Coello | An Introduction to Evolutionary Multi-Objective Optimization with Some Applications in Pattern Recognition[END_REF].

i i c P N = (31)
where c is a constant number greater than one, and Ni is the number of obtained Pareto's optimal response in the i-th section. This allows MOHTLBOGWO algorithm equations with higher probability to choose positions from low population sections.

Fuzzy Decision Making

The ultimate goal of a multi-objective optimization algorithm is to identify solving in the Pareto optimal set. However, identifying the entire Pareto optimal set is not possible due to its wide dimensions to prove the optimality, and therefore not recommended. Therefore, to investigate the feasibility of having Pareto optimal sets in multiple optimization problems, there is a practical approach. In the present study, a fuzzy approach is applied to select the best solution of the Pareto set. The j-objective function of a solution in Pareto set fj is defined by a membership function j as follows [START_REF] Taghavirashidizadeh | A Genetic algorithm for multi-objective reconfiguration of balanced and unbalanced distribution systems in fuzzy framework[END_REF]:

f f f f f f f f         =           - - (32)
where min j f and max j f are the minimum and maximum values of the objective function j. For each i solution, the membership function is calculated as follows [START_REF] Shirai | Configuration of fuzzy membership function circuits using current mode[END_REF]:

1 1 1 n i j j j m n i j i j    = = = =   ( 33
)
where n is the number of objective functions and m is the number of solutions. The answers have the maximum value of μi for the best compromise answer.

Implementation of MOHTLBOGWO

The steps to implement the proposed method in problem solving are as follows:

• Step 1: Random generation of the initial population from the set of variables including the location of the installation and the capacity of solar and wind units as well as the power factor of wind turbines.

• Step 2: the value of the target function is calculated for each set of variables and the best set is selected in terms of the value of the objective function as the representative of the entire population. • Step 3: Each set of variables is updated by the proposed algorithm, and if the new variables have better results, they will be replaced with the previous set. • Step 4: if the condition for convergence does not exist, go to step 1. Otherwise, go to step 5. • Step 5: Stop the algorithm.

Simulation Results and Discussion

Testing System

As shown in Figure 3, in this paper, the proposed method is implemented on the IEEE standard 33-buses radial network. In the 33-buses network, the total consumption of this network is 3720 kW and 2300 kVAR. The 33-buses Network has 37 branches. System information is presented in [START_REF] Etemad | Exploiting endusers capability in voltage support by means of power electronic transformer[END_REF].

Simulation Strategies

In this paper, the optimal placement of solar panels and wind turbines for the purpose of reducing losses and improving reliability indicators are presented using the proposed method HTLBOGWO. First, simulations were performed in a single-objective optimization. Then, multiobjective optimization was proposed based on fuzzy decision method and the results were compared and analyzed. In single-objective optimization, simulation results obtained from TLBO, GWO, and HTLBOGWO methods are evaluated based on indicators of loss reduction and energy not-supplied subscriber reduction (reliability improvement). In multi-objective optimization, the results of multiobjective methods MOGWO, MOTLBO and MOTLBOGWO are presented. In Figure 4, the LSF curve of the 33-buses network is presented. In this figure, the buses 6, 28, 29, 30, 9, 10, 13, 8, 27, 31 and 26 are based on the LSF curve considered as the candidate buses for installing wind turbines and solar panels Which are similar to Candida buses in reference [START_REF] Ali | Ant Lion Optimization Algorithm for optimal location and sizing of renewable distributed generations[END_REF].

Single-Objective Optimization Results

In this section, the results of single-objective positioning of solar panels and wind turbines are presented for the purpose of reducing losses and reducing ENS for 1, 2 and 3 DGs 1 megawatts. In Table 1, the results are presented with the objective of reducing losses. In the positioning of 3 solar panels, the buses 13, 24 and 30 were selected for installation with a capacity of 996, 938 and 852 kW respectively. The network losses before positioning were 202.68 kilowatts. After locating one, two, and three panels, this value has dropped to 127.28, 86.55 and 72.11 kW, respectively. In locating of 3 wind turbines, the buses 30, 13 and 24 were selected for installation with a capacity of 997, 1000 and 789 kW and with power factors 0.8659, 0.8122 and 0.8726 respectively. With locating 1, 2 and 3 wind turbines, the network losses dropped from 202.68 to 81.43, 32.17 and 13.68 kilowatts, respectively. With 1, 2 and 3 wind turbines, the network losses dropped from 202.68 to 81.43, 32.17 and 13.68 kilowatts, respectively. Also, the minimum voltage is also improved by the use of panels and turbines. The results also show that the ENS value before the locating of solar panels and wind turbines was 6.695 MW, and after locating 3 panels it reached 0.628 MW, and after locating 3 turbines reached 0.671 MWh. The results show that the cost of network losses with the optimal locating of solar panels and wind turbines has decreased and financial benefits from loss reduction have increased. In addition, the results have shown that increasing the number of DGs has led to a further reduction in losses, reduction cost of losses, and ENS, as well as further increase in the minimum voltage and net profit of the network. The results show that the performance of wind turbines was better than solar panels in improving the distribution network characteristics. In other words, in terms of reduction of losses and ENS, as well as the improvement of minimum voltage, it achieves better results due to the injection of reactive power in addition to the active power to the distribution network.

Table 1. Optimal sizing DG's for 33 bus IEEE system (active power loss objective function)

The convergence curve of the combined HTLBOGWO combination method, along with TLBO and GWO methods for single-objective optimization of solar panels and wind turbines with the goal of reducing losses for 3 DGs, is shown in Figure 5 and Figure 6, respectively. As can be seen, in locating solar panels and wind turbines, the HTLBOGWO method has less convergence fluctuations, and achieved a lower loss rate than the TLBO and GWO with a higher convergence rate and less repetition.

Single-objective optimization results with the goal of reducing losses using TLBO, GWO and HTLBOGWO methods for 3DG are presented in Table 2. As can be seen, the proposed method is better than other methods in terms of reducing losses, reducing cost of losses, reducing ENS, improving the minimum voltage and increasing the financial profit of the network.

Table 3 presents the results of single-objective positioning of solar panels and wind turbines with the goal of reducing ENS for 1, 2 and 3 DG 1 megawatts. As can be seen, with the locating of 1 and 2 solar panels, the network losses are decreasing, but in the positioning of the three panels, the losses compared to the positioning of the two turbines increased from 87. [START_REF] Tavana | An integrated intuitionistic fuzzy AHP and SWOT method for outsourcing reverse logistics[END_REF] contrast to the positioning of the three panels, the ENS value was strongly Dropped to 0.0025 MWh. In the positioning of wind turbines, with increasing its number, the amount of losses decreased, the minimum network voltage increased and the amount of ENS decreased, so that it reached 0.002 MW in the positioning of 3 turbines. Also, with the optimal locating of turbines, the cost of losses is reduced and the amount of financial gain in the locating of the three turbines has a maximum value. This table also shows that wind turbines have a better effect on the characteristics of the distribution network than solar panels. Single-objective optimization results with the aim of reducing the ENS obtained from the HTLBOGWO method are presented in Table 4 and compared with the TLBO and GWO methods. The amount of losses by the proposed method was 97.04 kilowatts and the losses obtained from the TLBO and GWO methods were 101.25 and 103.17 kW respectively. Also, the amount of ENS by the proposed method, 0.025 MWh and by TLBO and GWO methods were 0.048 and 0.044 MW, respectively, which indicates the superiority of the proposed method in solving the problem of positioning solar panels. Also, according to Table 4, in the locating of solar panels and wind turbines, the HTLBOGWO method performed better than other methods in terms of reducing losses, reducing ENS, improving the minimum voltage and financial profit of the network.

Assessment of the Renewable Resource Dispersion Effect

In this section, the effect of dispersing the capacity of solar panels and wind turbines has been evaluated in problem solving. In other words, the effect of using a 3 MW DG with the use of three DGs with a capacity of 1 MW on the power loss, reliability index, minimum network voltage, cost of losses and financial profit of the network has been evaluated. The results obtained in Table 5 and Table 6 are presented for the purpose of decreasing losses and decreasing the ENS in a single-objective optimization, respectively. As can be seen, when using three DGs with a capacity of 1 MW instead of using a 3 MW DG The amount of network losses and ENS will be lower and the network obtaining the higher minimum voltage , more financial gain and reduce cost of losses. For example, according to Table 5, in dispersal conditions, the network losses decreased from 103.96 to 72.11 kV and ENS values decreased from 4.379 to 0.628 megawatts / h. According to Table 6, under dispersal conditions, the amount of network losses decreased from 161.26 to 97.04 kilowatts and ENS values from decreased 0.11 to 0.025 MWh. The voltage profile curve of 33-buses network before locating and under dispersion conditions for solar panels and wind turbines is plotted in Figure 7. Also Figure 7 shows the voltage profile curve of 33-buses network before locating and under conditions where there is no dispersion condition. Figure 7 and Figure 8 showed that the distribution and distribution of power in the network buses leads to improvement of the network voltage profile.

889 [START_REF] Vosoogha | An Intelligent Power Prediction Method for Wind Energy Generation Based on Optimized Fuzzy System[END_REF], 903 [START_REF] Khatibi | Comparative Analysis of Power System Model Reduction[END_REF], 994 [START_REF] Prabha | Optimal placement and sizing of multiple distributed generating units in distribution networks by invasive weed optimization algorithm[END_REF] 996/1@ 13, 938/1@ 24, 852/1@ 30

Multi-Objective Optimization Results

In this section, multi-objective optimization of the locating of solar panels and wind turbines is presented based on the fuzzy decision approach using the HTLBOGWO method and the results are compared with the TLBO and GWO multi-objective methods. The Pareto optimal set of solutions for placement of 1 panel and 1 wind turbine with a capacity of 1 MW is shown in Figure 9 and Figure 10, respectively. As shown in Figure 9, the MOHTLBOGWO method, in comparison with the MOTLBO and MOGWO response scattering, offers better results in terms of achieving an optimal response. Also, in the terms of achieving to lower losses, lower ENS, minimum voltage, lower cost of losses and higher financial profit, the proposed method performance was better than the other two methods. Also, according to Table 7 and Table 8, in the multi-objective optimization, items such as the losses, ENS, the minimum voltage, the cost of losses, and the amount of financial gain are located between the single-objective values that aimed at reducing losses and reducing ENS. In other words, in the multi-objective optimization, there is a Compromise between the losses and ENS. In Figure 11, the multi-objective positioning of both panels and turbines simultaneously is presented using various methods. The numerical results are also given in Table 9. The results show that the proposed MOTLBOGWO method offers better Pareto-optimal solution than other methods and has obtained the optimal solution with less dispersion. Also, according to Table 9, can be seen that the simultaneous use of 1 panel and 1 turbine instead of use one panel or one turbine, has led to a further reduction of losses and ENS, a further increase in the minimum voltage and net profit of the network.

Conclusions

In this paper, the new MOHTLBOGWO method was used to determine the location and capacity of solar panels and wind turbines for the purpose of reducing losses and improving reliability in the distribution network. The proposed problem was evaluated based on single-objective and multi-objective optimization based on fuzzy decision making in optimizing the use of renewable resources in the 33 Buses distribution network. The effectiveness of the proposed method in various optimizations was presented and compared with the MOTLBO and MOGWO results. The results showed that in the optimal placement of the solar panels and wind turbines, unlike their single-objective optimization, there is a compromised between the reduction of losses and improved reliability, and the optimal solution was determined from the set of answers of the Pareto levels. In comparison with the MOTLBO and MOGWO methods, the MOHTLBOGWO method has a better convergence in Pareto levels and achieves optimal solutions. Also, the ability of the proposed method in single-objective optimization In comparison with previous studies in reducing losses, improving the minimum voltage and increasing network profits was confirmed. The results also showed that the dispersal of renewable resources in the distribution network would result in more reduction of losses, greater improvement on reliability, and more financial gain on the network.

Figure 1 .

 1 Figure 1. Distribution Network Equivalent Circuit [32] 22 k k ik ik loss 2 k

Figure 2 .

 2 Figure 2. Leader hierarchy in gray wolves group The principles of GWO function are as follows:

  vectors of random numbers at a distance are shown by and that are selected from [0, 1] area and vector a decreases from 2 to 0 during repetitions.

Figure 3 .

 3 Figure 3. IEEE Standard Radial 33-buses Distribution Network

Figure 4 .

 4 Figure 4. LSF curve of the 33-buses network

Figure 5 .

 5 Figure 5. Convergence curve of different methods with the aim of reducing the losses in the 3 PV positioning

Figure 6 .

 6 Figure 6. Convergence curve of different methods with the aim of reducing losses in the 3 WT positioning

Figure 7 .

 7 Figure 7. Network voltage profile of 33 buses before placing solar panels and under dispersion conditions

Figure 8 .

 8 Figure 8. Network voltage profile of 33 buses before placing wind turbines and under dispersion conditions

Figure 9 .Figure 10 .

 910 Figure 9. Pareto optimal answers set curve in the placement of the solar panels by various methods

Figure 11 .

 11 Figure 11. Pareto optimal answers set curve (PV + WT only with maximum capacity of 3 megawatts)

  to 97.04 kilowatt, but in

			Maximum Size 3MW (kW/pf/@Bus)
		No PV	1PV	2PV	3PV
	Power Losses (kW)	210.98	127.28	86.55	72.11
	ENS (MWh/yr)	6.695	4.158	1.923	0.628
	Minimum Voltage (pu)	0.91308	0.9285	0.9629	0.9667
					996/1@ 13,
				893/1@ 10,	938/1@ 24,
	Size and Location	--	1000/1@30	1000/1@ 30	852/1@ 30
	Total Size (kW)	--	1000	1893
	Losses Cost ($/yr)	110891	66899.8	45491.3	37903.5
	Net Saving ($/yr)	--	43991.2	65399.7	72987.5
		No WT	1WT	2WT	3WT
	Power Losses (kW)	210.98	81.43	32.17	13.68
	ENS (MWh/yr)	6.695	4.158	2.233	0.640
	Minimum Voltage (pu)	0.91308	0.9360	0.9796	0.9892
					997/0.8659 @30,
				861/ 0.8742 @10,	1000/0.8122 @13,
	Size and Location	--	1000/0.8011@30	1000/0.8091@30	789/ 0.8726@24
	Total Size (kW)	--
	Losses Cost ($/yr)	110891	42801	16913.5	7191.6
	Net Saving ($/yr)	--	68090	93977.5	103699.4

Table 2 .

 2 Optimal sizing DG's for 33-bus IEEE system (active power loss objective function-algorithms comparison)

			Maximum Size 3MW (kW/pf/@Bus)
		No PV	TLBO	GWO	HTLBOGWO
	Power Losses (kW)	210.98	74.10	73.29	72.11
	ENS (MWh/yr)	6.695	0.734	0.659	0.628
	Minimum Voltage (pu)	0.9130	0.9621	0.9648	0.9667
			872 (29), 858 (13), 946
	Size and Location	--

Table 3 .

 3 Optimal sizing DG's for 33 bus IEEE system (ENS objective function)

	Total Size (kW)	--
	Losses Cost ($/yr)	110891	38949.0	38525.7	37903.5
	Net Saving ($/yr)	--	71942	72365.3	72987.5
		No WT	TLBO	GWO	HTLBOGWO

Table 4 .

 4 Optimal sizing DG's for 33 bus IEEE system (ENS objective function-Algorithms comparison)

			Maximum Size 3MW (kW/pf/@Bus)
		No PV	TLBO	GWO	HTLBOGWO
	Power Losses (kW)	210.98	101.25	103.17	97.04
	ENS (kWh/yr)	6.695	0.048	0.044	0.035
	Minimum Voltage (pu)	0.9130	0.9694	0.9677	0.9705
			986 (9),	807 (29),	709 (30),
	Size and Location	--	789 (29),	830 (13),	810 (10),
			1000 (13)	1000 (10)	957 (13)
	Total Size (kW)	--
	Losses Cost ($/yr)	110891	53217.5	54228.8	51008.3
	Net Saving ($/yr)	--	57673.5	56662.2	59882.7
		No WT	TLBO	GWO	HTLBOGWO
	Power Losses (kW)	210.98	58.17	59.61	55.43
	ENS (kWh/yr)	6.695	0.018	0.016	0.011
	Minimum Voltage (pu)	0.9130	0.9808	0.9800	0.98248
			893	652	993
			(9) 0.8760, 762	(29) 0.9660, 978	(13) 0.8579, 802
	Size and Location	--	(30) 0.982,	(13) 0.9465, 920	(29) 0.8427, 762
			1000	(30) 0.9824	(30) 0.8209
			(13) 0.9519
	Total Size (kW)	--
	Losses Cost ($/yr)	110891	30577.0	31334.6	29136.6
	Net Saving ($/yr)	--	80314	79556.4	81754.4

Table 5 .

 5 Results of DGs distribufication (SOO-Loss minimization)

	Maximum Size 3MW (kW/pf/@Bus)

Table 6 .

 6 Results of DGs distribufication (SOO-ENS minimization)

	Maximum Size 3MW (kW/pf/@Bus)
	1DG (3000)	3DG (1000)
	(kVA/P.F)	(kVA/P.F)

Table 7 .

 7 Fuzzy multi-objective results for 1 PV placement Maximum Size 3MW (kW/pf/@Bus)

		No DG	1PV (Loss objective)	1PV (ENS objective)	1PV (Loss+ENS bjective)
	Power Losses (kW)	210.98	101.96	161.26	124.62
	ENS (kWh/yr)	6.695	4.379	1.337	1.602
	Minimum Voltage (pu)	0.9130	0.9510	0.9466	0.9486
	Size and Location	--	2585/1@29	2097/1@10	2007/1@30
	Total Size (kW)	--	2585	2097	2007
	Losses Cost ($/yr)	110891	53590.1	84761.3	65505.4
	Net Saving ($/yr)	--	57300.9	26129.7	45385.6

Table 8 .

 8 Fuzzy multi-objective results for 1 WT placement Maximum Size 3MW (kW/pf/@Bus)

		No DG	1WT (Loss objective)	1WT (ENS objective)	1WT (Loss+ENS objective)
	Power Losses (kW)	210.98	61.36	100.09	73.3816
	ENS (kWh/yr)	6.695	4.406	0	1422.4586
	Minimum Voltage (pu)	0.9130	0.9667	0.9564	0.95677

Table 9 .

 9 Fuzzy multi-objective results for 1 PV and 1 WT placement simultaneously Maximum Size 3MW (kW/pf/@Bus)

		No DG	1 PV	1 WT	1PV+1WT
	Power Losses (kW)	210.98	124.6299	73.3816	38.3211
	ENS (kWh/yr)	6.695	1.602	1.422	0.342
	Minimum Voltage (pu)	0.9130	0.9486	0.9567	0.9805
	Size and Location	--	2007/1@30	2078/0.833@30	873/1@13 1405/0.8@30
	Total Size (kW)	--	2007	2078	2278
	Losses Cost ($/yr)	110891	65505.4	38569.3	20141.5
	Net Saving ($/yr)	--	45385.6	72321.7	90749.5


Mérouane Debbah
email: merouane.debbah@supelec.fr

ON THE LIMITING MOMENTS OF VANDERMONDE RANDOM MATRICES

In this paper, analytical methods for finding moments of random Vandermonde matrices are developed. Vandermonde Matrices play an important role in signal processing and communication applications such as direction of arrival estimation, sampling theory or precoding. Within this framework, we extend classical freeness results on random matrices with i.i.d entries and show that Vandermonde structured matrices can be treated in the same vein with different tools. We focus on various types of Vandermonde matrices, namely Vandermonde matrices with or without uniformly distributed phases. In each case, we provide explicit expressions of the moments of the associated Gram matrix, as well as more advanced models involving the Vandermonde matrix. Comparisons with classical i.i.d. random matrix theory are provided and free deconvolution results are also discussed.

INTRODUCTION

We will consider Vandermonde matrices V of dimension N × L of the form

V = 1 √ N      1 • • • 1 e -jω1
• • • e -jωL . . . . . . . . . e -j(N -1)ω1 • • • e -j(N -1)ωL

     (1)
where ω 1 ,...,ω L are independent and identically distributed (phases) taking values on [0, 2π). Such matrices occur frequently in many applications, such as finance [START_REF] Norberg | On the Vandermonde matrix and its application in mathematical finance[END_REF], signal array processing [START_REF] Schmidt | Multiple emitter localization and signal parameter estimation[END_REF][START_REF] Wax | Detection of signals by information theoretic criteria[END_REF][START_REF] Johnson | Array Signal processing: Concepts and Techniques[END_REF][START_REF] Roy | ESPRIT-estimation of signal parameters via rotational invariance techniques[END_REF][START_REF] Porat | Analysis of the asymptotic relative efficiency of the MUSIC algorithm[END_REF], ARMA processes [START_REF] Klein | On Stein's equation, Vandermonde matrices and Fisher's information matrix of time series processes. part I: The autoregressive moving average process[END_REF], cognitive radio [START_REF] Sampaio | Vandermonde matrices for security applications[END_REF], security [START_REF] Sampaio | Vandermonde frequency division multiplexing[END_REF], wireless communications [START_REF] Wang | Vandermonde-Lagrange mutually orthogonal flexible transceivers for blind CDMA in unknown multipath[END_REF] and biology [START_REF] Waterfall | Sloppymodel universality class and the Vandermonde matrix[END_REF] and have been much studied. The main results are related to the distribution of the determinant of (1) [START_REF] Girko | Theory of Random Determinants[END_REF]. The large majority of known results on the eigenvalues of the associated Gram matrix concern Gaussian matrices [START_REF] Mehta | Random Matrices[END_REF] or matrices with independent entries. None have dealt with the Vandermonde case. For the Vandermonde case, the results depend heavily on the distribution of the entries, and do not give any hint on the asymptotic behaviour as the matrices become large. In the realm of wireless channel modelling, [START_REF] Muller | A random matrix model of communication via antenna arrays[END_REF] has provided some insight on the behaviour of the eigenvalues of random Vandermonde matrices for a specific case, without any formal proof. We prove here that the case is in fact more involved than what was claimed.

In many applications, N and L are quite large, and we may be interested in studying the case where both go to ∞ at a given ratio, with L N → c. Results in the literature say very little on the asymptotic behaviour of (1) under this growth condition. The results, however, are well known for other models. The factor 1 √ N , as well as the assumption that the Vandermonde entries e -jωi lie on the unit circle, are included in [START_REF] Norberg | On the Vandermonde matrix and its application in mathematical finance[END_REF] to ensure that our analysis will give limiting asymptotic behaviour. Without this assumption, the problem at hand is more involved, since the rows of the Vandermonde matrix with the highest powers would dominate in the calculations of the moments when the matrices grow large, and also grow faster to infinity than the 1 √ N factor in [START_REF] Norberg | On the Vandermonde matrix and its application in mathematical finance[END_REF], making asymptotic analysis difficult. In general, often the moments, not the moments of the determinants, are the quantities we seek. Results in the literature also say very little on the moments of Vandermonde matrices. The literature says very little on the mixed moments of Vandermonde matrices and matrices independent from them. This is in contrast to Gaussian matrices, where exact expressions [START_REF] Thorbjørnsen | Mixed moments of Voiculescu's Gaussian random matrices[END_REF] and their asymptotic behaviour [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF] are known using the concept of freeness [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF] which is central for describing the mixed moments.

The derivation of the moments are a useful basis for performing deconvolution. For Gaussian matrices, deconvolution has been handled in the literature [START_REF] Anderson | Asymptotic theory for principal component analysis[END_REF][START_REF] Abed-Meraim | A subspace algorithm for certain blind identification problems[END_REF][START_REF] Thorbjørnsen | Mixed moments of Voiculescu's Gaussian random matrices[END_REF][START_REF] Ryan | Free deconvolution for signal processing applications[END_REF]. Similar flavored results will here be proved for Vandermonde matrices. Concerning the moments, it will be the asymptotic moments of random matrices of the form V H V which will be studied, where (.) H denotes hermitian transpose. We will also consider mixed moments of the form DV H V, where D is a square diagonal matrix independent from V. While we provide the full computation of lower order moments, we also describe how the higher order moments can be computed. Tedious evaluation of many integrals is needed for this, but numerical methods can also be applied. Surprisingly, it turns out that the first three limit moments can be expressed in terms of the Marchenko Pastur law [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF][START_REF] Tulino | Random Matrix Theory and Wireless Communications[END_REF]. For higher order moments this is not the case, although we state an interesting inequality involving the Vandermonde limit moments and the moments of the classical Poisson distribution and the Marchenko Pastur law, also known as the free Poisson distribution [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF].

This paper is organized as follows: Section 2 contains a general result for the mixed moments of Vandermonde matrices and matrices independent from them. We will differ between the case where the phase ω in (1) are uniformly distributed on [0.2π), and the more general cases. The case of uniformly distributed phases is handled in section 3. In this case it turns out that one can have very nice expressions, for both the asymptotic moments, as well as for the lower order moments. Section 4 considers the more general case when ω has a continous density, and shows how the asymptotics can be described in terms of the case when ω is uniformly distributed. Section 5 discusses our results and puts them in a general deconvolution perspective, comparing with other deconvolution results, such as those for Gaussian deconvolution.

In the following, upper (lower boldface) symbols will be used for matrices (column vectors) whereas lower symbols will represent scalar values, (.) T will denote transpose operator, (.) ⋆ conjugation and (.) H = (.) T ⋆ hermitian transpose. I n will represent the n × n identity matrix. We let tr n be the normalized trace for matrices of order n × n, and T r the non-normalized trace. V will be used only to denote Vandermonde matrices with a given phase distribution. The dimensions of the Vandermonde matrices will always be N × L unless otherwise stated, and the phase distribution of the Vandermonde matrices will always be denoted by ω.

A GENERAL RESULT FOR THE MIXED MOMENTS OF VANDERMONDE MATRICES

We first state a general theorem applicable to Vandermonde matrices with any phase distribution. The proof for this theorem, as well as for theorems succeeding it, are based on calculations where partitions are highly involved. We denote by P(n) the set of all partitions of {1, ..., n}, and we will use ρ as notation for a partition in P(n). The set of partitions will be equipped with the refinement order ≤, i.e. ρ 1 ≤ ρ 2 if and only if any block of ρ 1 is contained within a block of ρ 2 . Also, we will write ρ = {ρ 1 , ..., ρ k }, where ρ j are the blocks of ρ, and let |ρ| denote the number of blocks in ρ. We denote by 0 n the partition with n blocks, and by 1 n the partition with 1 block.

In the following D r (N ), 1 ≤ r ≤ n are diagonal L × L matrices, and V is of the form (1). We will attempt to find

M n = lim N →∞ E[tr L ( D 1 (N )V H VD 2 (N )V H V • • • × D n (N )V H V)]
(2) for many types of Vandermonde matrices, under the assumption that L N → c, and under the assumption that the D r (N ) have a joint limit distribution as N → ∞ in the following sense:

Definition 1 We will say that the {D r (N )} 1≤r≤n have a joint limit distribution as N → ∞ if the limit

D i1,...,is = lim N →∞ tr L (D i1 (N ) • • • D is (N )) (3)
exists for all choices of i 1 , ..., i s . For ρ = {ρ 1 , ..., ρ k }, with ρ i = {ρ i1 , ..., ρ i|ρi| }, we also define

D ρi = D iρ i1 ,...,iρ i|ρ i | ,
and

D ρ = k i=1 D ρi .
Had we replaced Vandermonde matrices with Gaussian matrices, free deconvolution results [START_REF] Ryan | Free deconvolution for signal processing applications[END_REF] could help us compute the quantities D i1,...,is from M n . For this, the cumulants of the Gaussian matrices are needed, which asymptotically have a very nice form. For Vandermonde matrices, the role of cumulants is taken by the following quantites

Definition 2 Define K ρ,ω,N = 1 N n+1-|ρ| × (0,2π) |ρ| n k=1 1-e jN (ω b(k-1) -ω b(k) ) 1-e j(ω b(k-1) -ω b(k) ) , dω 1 • • • dω |ρ| , ( 4
)
where ω ρ1 , ..., ω ρ |ρ| are i.i.d. (indexed by the blocks of ρ), all with the same distribution as ω, and where b(k) is the block of ρ which contains k (where notation is cyclic, i.e. b(-1) = b(n)). If the limit

K ρ,ω = lim N →∞ K ρ,ω,N
exists, then K ρ,ω is called a Vandermonde mixed moment expansion coefficient. These coefficients will for Vandermonde matrices play the same role as the cumulants do for large Gaussian matrices. We will not call them cumulants, however, since they don't share the same multiplicative properties (embodied in what is called the moment cumulant formula).

The following is the main result of the paper. Different versions of it adapted to different Vandermonde matrices will be stated in the succeeding sections.

Theorem 1 Assume that the {D r (N )} 1≤r≤n have a joint limit distribution as N → ∞. Assume also that all Vandermonde mixed moment expansion coefficients K ρ,ω exist.

Then the limit

M n = lim N →∞ E[tr L ( D 1 (N )V H VD 2 (N )V H V • • • × D n (N )V H V)]
(5) also exists when L N → c, and equals

ρ∈P(n) K ρ,ω c |ρ|-1 D ρ . (6)
The proof of theorem 1 can be found in [START_REF] Ryan | Random Vandermonde matrices-part I: Fundamental results[END_REF]. Although the limit of K ρ,ω,N as N → ∞ may not exist, it will be clear from section 4 that it exists when the density of ω is continous. Theorem 1 explains how convolution with Vandermonde matrices can be performed, and also provides us an extension of the concept of free convolution to Vandermonde matrices. Note that when

D 1 (N ) = • • • = D n (N ) = I L , we have that M n = lim N →∞ E tr L V H V n ,
so that our our results also include the limit moments of the Vandermonde matrices themselves. M n corresponds also to the limit moments of the empirical eigenvalue distribution

F N V H V defined by F N V H V (λ) = #{i|λ i ≤ λ} N ,
(where λ i are the (random) eigenvalues of V H V), i.e.

M n = lim N →∞ E λ n dF N (λ) .
(6) will also be useful on the scaled form

cM n = ρ∈P(n) K ρ,ω (cD) ρ . (7)
When

D 1 (N ) = D 2 (N ) = • • • = D n (N )
, we denote their common value D(N ), and define the sequence

D = (D 1 , D 2 , ...) with D n = lim N →∞ tr L ((D(N )) n ).
In this case D ρ does only depend on the block cardinalities |ρ j |, so that we can group together the K ρ,ω for ρ with equal block cardinalities. If we group the blocks of ρ so that their cardinalities are in descending order, and set

P(n) r1,r2,...,r k = {ρ = {ρ 1 , ..., ρ k } ∈ P(n)||ρ i | = r i ∀i},
where

r 1 ≥ r 2 ≥ • • • ≥ r k ,
and also write

K r1,r2,...,r k = ρ∈P(n)r 1 ,r 2 ,...,r k K ρ,ω , (8)
then, after performing the substitutions

m n = (cM ) n = c lim N →∞ E tr L D(N )V H V n , d n = (cD) n = c lim N →∞ tr L (D n (N )) , (9)
(7) can be written

m n = r 1 ,...,r k r1+•••+r k =n K r1,r2,...,r k k j=1 d rj . (10)
For the first 5 moments this becomes

m 1 = K 1 d 1 m 2 = K 2 d 2 + K 1,1 d 2 1 m 3 = K 3 d 3 + K 2,1 d 2 d 2 1 + K 1,1,1 d 3 1 m 4 = K 4 d 4 + K 3,1 d 3 d 1 + K 2,2 d 2 2 + K 2,1,1 d 2 d 2 1 + K 1,1,1,1 d 4 1 m 5 = K 5 d 5 + K 4,1 d 4 d 1 + +K 3,2 d 3 d 2 + K 3,1,1 d 3 d 2 1 + K 2,2,1 d 2 2 d 1 + K 2,1,1,1 d 2 d 3 1 + K 1,1,1,1,1 d 5 1 . . . . . . (11
) Thus, the algorithm for computing the asymptotic mixed moments of Vandermonde matrices with matrices independent from them can be split in two:

• (9), which scales with the matrix aspect ratio c, and [START_REF] Waterfall | Sloppymodel universality class and the Vandermonde matrix[END_REF], which performs computations independent of the matrix aspect ratio c.

• (
Similar splitting of the algorithm for computing the asymptotic mixed moments of Wishart matrices and matrices independent from them was derived in [START_REF] Ryan | Free deconvolution for signal processing applications[END_REF]. Although the matrices D i (N ) are assumed to be determinstic matrices throughout the paper, all formulas extend naturally to the case when D i (N ) are random matrices independent from V. The only difference when the D i (N ) are random is that certain quantities are replaced with fluctuations. D 1 D 2 should for instance be replaced with

lim N →∞ E tr L (D(N )) tr L (D(N )) 2 when D i (N ) is random.
In the next sections, we will derive and analyze the Vandermonde mixed moment expansion coefficients K ρ,ω for various cases, which is essential for the the algorithm (11).

UNIFORMLY DISTRIBUTED ω

We will let u denote the uniform distribution on [0, 2π). We can write

K ρ,u,N = 1 (2π) |ρ| N n+1-|ρ| × (0,2π) |ρ| n k=1 1-e jN (x b(k-1) -x b(k) ) 1-e j(x b(k-1) -x b(k) ) dx 1 • • • dx |ρ| , ( 12
)
where integration is w.r.t. Lebesgue measure. In this case one particular class of partitions will be useful to us, the noncrossing partitions: Definition 3 A partition is said to be noncrossing if, whenever i < j < k < l, i and k are in the same block, and also j and l are in the same block, then all i, j, k, l are in the same block. The set of noncrossing partitions is denoted by N C(n).

The noncrossing partitions have already shown their usefulness in expressing the freeness relation in a particularly nice way [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF]. Their appearance here is somewhat different than in the case for the relation to freeness: Theorem 2 Assume that the {D r (N )} 1≤r≤n have a joint limit distribution as N → ∞, Then the Vandermonde mixed moment expansion coefficient

K ρ,u = lim N →∞ K ρ,u,N
exists for all ρ. Moreover, 0 < K ρ,u ≤ 1, the K ρ,u are rational numbers for all ρ, and K ρ,u = 1 if and only if ρ is noncrossing.

The proof of theorem 2 can be found in [START_REF] Ryan | Random Vandermonde matrices-part I: Fundamental results[END_REF]. Due to theorem 1, theorem 2 guarantees that the asymptotic mixed moments (5) exist when L N → c for uniform phase distribution, and are given by [START_REF] Porat | Analysis of the asymptotic relative efficiency of the MUSIC algorithm[END_REF]. The values K ρ,u are in general hard to compute for higher order ρ with crossings. We have performed some of these computations. It turns out that the following computations suffice to obtain the 7 first moments.

Lemma 1

The following holds:

K {{1,3},{2,4}},u = 2 3 K {{1,4},{2,5},{3,6}},u = 1 2 K {{1,4},{2,6},{3,5}},u = 1 2 K {{1,3,5},{2,4,6}},u = 11 20 K {{1,5},{3,7},{2,4,6}},u = 9 20 K {{1,6},{2,4},{3,5,7}},u = 9 20 .
The proof of lemma 1 is given in [START_REF] Ryan | Random Vandermonde matrices-part I: Fundamental results[END_REF]. Combining theorem 2 and lemma 1 into this form, we will prove the following:

Theorem 3 Assume D 1 (N ) = D 2 (N ) = • • • = D n (N ).
When ω = u, [START_REF] Waterfall | Sloppymodel universality class and the Vandermonde matrix[END_REF] takes the form

m 1 = d 1 m 2 = d 2 + d 2 1 m 3 = d 3 + 3d 2 d 1 + d 3 1 m 4 = d 4 + 4d 3 d 1 + 8 3 d 2 2 + 6d 2 d 2 1 + d 4 1 m 5 = d 5 + 5d 4 d 1 + 25 3 d 3 d 2 + 10d 3 d 2 1 + 40 3 d 2 2 d 1 + 10d 2 d 3 1 + d 5 1 m 6 = d 6 + 6d 5 d 1 + 12d 4 d 2 + 15d 4 d 2 1 + 151 20 d 2 3 + 50d 3 d 2 d 1 + 20d 3 d 3 1 + 11d 3 2 + 40d 2 2 d 2 1 + 15d 2 d 4 1 + d 6 1 m 7 = d 7 + 7d 6 d 1 + 49 3 d 5 d 2 + 21d 5 d 2 1 + 497 20 d 4 d 3 + 84d 4 d 2 d 1 + 35d 4 d 3 1 + 1057 20 d 2 3 d 1 + 693 10 d 3 d 2 2 + 175d 3 d 2 d 2 1 + 35d 3 d 4 1 + 77d 3 2 d 1 + 280 3 d 2 2 d 3 1 + 21d 2 d 5 1 + d 7 1 .
Theorem 2 and lemma 1 reduces the proof of theorem 3 to a simple count of partitions. Theorem 3 is proved in [START_REF] Ryan | Random Vandermonde matrices-part I: Fundamental results[END_REF]. To compute higher moments m k , K ρ,u must be computed for partitions of higher order.

Following the proof of theorem 2, we can also obtain formulas for the fluctuations of mixed moments of Vandermonde matrices. We will not go into details on this, but only state the following equations without proof:

lim N →∞ E tr L D(N )V H V n tr L D(N )V H V m = E tr L D(N )V H V n D m 1 c lim N →∞ E T r D(N )V H V 2 tr L D(N )V H V 2 = 4 3 d 2 2 + 4d 2 d 2 1 + 4d 3 d 1 + d 4 . (13)
Following the proof of theorem 2 again, we can also obtain exact expressions for moments of lower order random Vandermonde matrices with uniformly distributed phases, not only the limit. We state these only for the first four moments.

Theorem 4 Assume D 1 (N ) = D 2 (N ) = • • • = D n (N ).
When ω = u, [START_REF] Waterfall | Sloppymodel universality class and the Vandermonde matrix[END_REF] takes the exact form

m 1 = d 1 m 2 = 1 -N -1 d 2 + d 2 1 m 3 = 1 -3N -1 + 2N -2 d 3 +3 1 -N -1 d 1 d 2 + d 3 1 m 4 = 1 - 20 3 N -1 + 11N -2 - 37 6 N -3 d 4 + 4 -12N -1 + 8N -2 d 3 d 1 + 8 3 -5N -1 + 19 6 N -2 d 2 2 +6 1 -N -1 d 2 d 2 1 + d 4 1 .
Theorem 4 is proved in [START_REF] Ryan | Random Vandermonde matrices-part I: Fundamental results[END_REF]. Exact formulas for the higher order moments also exist, but they become increasingly complex, as entries for higher order terms L -k also enter the picture. These formulas are also harder to prove for higher order moments. In many cases, exact expressions are not what we need: First order approximations (i.e. expressions where only the L -1 -terms are included) can suffice for many purposes. In [START_REF] Ryan | Random Vandermonde matrices-part I: Fundamental results[END_REF], we explain how the simpler case of these first order approximations can be computed. It seems much harder to prove a similar result when the phases are not uniformly distributed.

ω WITH CONTINOUS DENSITY

The following result tells us that the limit K ρ,ω exists for many ω, and also gives a useful expression for them in terms of the density of ω, and K ρ,u .

Theorem 5 The Vandermonde mixed moment expansion coefficients K ρ,ω = lim N →∞ K ρ,ω,N exist whenever the density p ω of ω is continous on [0, 2π). If this is fulfilled, then

K ρ,ω = K ρ,u (2π) |ρ|-1 2π 0 p ω (x) |ρ| dx . ( 14
)
The proof is given in [START_REF] Ryan | Random Vandermonde matrices-part I: Fundamental results[END_REF].

Besides providing us with a deconvolution method for finding the mixed moments of the {D r (N )} 1≤r≤n , theorem 5 also provides us with a way of inspecting the phase distribution ω, by first finding the moments of the density, i.e. 2π 0 p ω (x) k dx. However, note that we can not expect to find the density of ω itself, only the density of the density of ω. To see this, define

Q ω (x) = µ ({x|p ω ≤ x})
for 0 ≤ x ≤ ∞, where µ is uniform measure on the unit circle. Write also q ω (x) as the corresponding density, so that q ω (x) is the density of the density of ω. Then it is clear that

2π 0 p ω (x) |ρ| dx = ∞ 0 x n q ω (x)dx. (15)
These quantities correspond to the moments of the measure with density q ω , which can help us obtain the density q ω itself (i.e. the density of the density of ω). However, the density p ω can not be obtained, since we see that any reorganization of its values which do not change its density q ω will provide the same values in [START_REF] Thorbjørnsen | Mixed moments of Voiculescu's Gaussian random matrices[END_REF]. Note also that theorem 5 gives a very special role to the uniform phase distribution, in the sense that it minimizes the moments of the Vandermonde matrices V H V. This follows from ( 14), since

2π 0 p u (x) |ρ| dx ≤ 2π 0 p ω (x) |ρ| dx
for any density p ω . In [START_REF] Ryan | Random Vandermonde matrices-part II: Applications[END_REF], several examples are provided where the integrals ( 14) are computed.

DISCUSSION

We have already explained that one can perform deconvolution with Vandermonde matrices in a similar way to how one can perform deconvolution for Gaussian matrices. We have, however, also seen that there are many differences.

Convergence rates

In [START_REF] Thorbjørnsen | Mixed moments of Voiculescu's Gaussian random matrices[END_REF], almost sure convergence of Gaussian matrices was shown by proving exact formulas for the distribution of lower order Gaussian matrices. These deviated from their limits by terms of the form 1/L 2 . In theorem 4, we see that terms of the form 1/L are involved, which indicates that we can not hope for almost sure convergence of Vandermonde matrices. There is no reason why Vandermonde matrices should have the almost sure convergence property, due to their very different degree of randomness when compared to Gaussian matrices. Figures 1,2 show the speed of convergence of the moments of Vandermonde matrices (with uniformly distributed phases) towards the asymptotic moments as the matrix dimensions grow, and as the number of samples grow. The differences between the asymptotic moments and the exact moments are also shown. To be more precise, the MSE of figures 1 and 2 is computed as follows:

1. K samples V i are independently generated using (1).

The 4 first sample moments

mji = 1 L tr n V H i V i j
(1 ≤ j ≤ 4) are computed from the samples. 4. The 4 first exact moments E j are computed using theorem 4.

5. The 4 first asymptotic moments A j are computed using theorem 3.

6. The mean squared error (MSE) of the first 4 estimated moments from the exact moments is computed as

4 j=1
Mj -E j 2 .

7. The MSE of the first 4 exact moments from the asymptotic moments is computed as

4 j=1 (E j -A j ) 2 .
Figures 1 and2 are in sharp contrast with Gaussian matrices, as shown in figure 3. First of all, it is seen that the asymptotic moments can be used just as well instead of the exact moments (for which expressions can be found in [START_REF] Ryan | Channel capacity estimation using free probability theory[END_REF]), due to the O(1/N 2 ) convergence of the moments. Secondly, it is seen that only 5 samples were needed to get a reliable estimate for the moments.

Inequalities between moments of Vandermonde matrices and moments of known distributions

We will state an inequality involving the moments of Vandermonde matrices, and the moments of known distributions from probability theory. The classical Poisson distribution with rate λ and jump size α is defined as the limit Matrices are on the form 1 N XX H with X a complex standard Gaussian matrix. The MSE of the moments and the asymptotic moments is also shown. of

1 - λ n δ 0 + λ n δ α * N
as n → ∞ [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF]. For our analysis, we will only need the classical Poisson distribution with rate c and jump size 1.

We will denote this quantity by ν c . The free Poisson distribution with rate λ and jump size α is defined similarly as the limit of

1 - λ n δ 0 + λ n δ α ⊞N
as n → ∞, where ⊞ is the free probability counterpart of classical additive convolution [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF][START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF]. For our analysis, we will only need the free Poisson distribution with rate 1 c and jump size c. We will denote this quantity by µ c . µ c is the same as the better known Marchenko Pastur law, i.e. it has the density [START_REF] Hiai | The Semicircle Law, Free Random Variables and Entropy[END_REF]

f µc (x) = (1 - 1 c ) + δ 0 (x) + (x -a) + (b -x) + 2πcx , (16)
where

(z) + = max(0, z), a = (1 - √ c) 2 , b = (1 + √ c) 2 .
Since the classical (free) cumulants of the classical (free) Poisson distribution are λα n [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF], we see that the (classical) cumulants of ν c are c, c, c, c, ..., and that the (free) cumulants of µ c are 1, c, c 2 , c 3 , .... In other words, if a 1 has the distribution µ c , then

φ(a n 1 ) = ρ∈N C(n) c n-|ρ| = ρ∈N C(n) c |K(ρ)|-1 = ρ∈N C(n) c |ρ|-1 . (17)
Here we have used the Kreweras complementation map, which is an order-reversing isomorphism of N C(n) which satisfies |ρ| + |K(ρ)| = n + 1 (here φ is the expectation in a non-commutative probability space). Also, if a 2 has the distribution ν c , then

E(a n 2 ) = ρ∈P(n) c |ρ| . (18)
We immediately recognize the c |ρ|-1 -entry of theorem 1 in ( 17) and ( 18) (except for an additional power of c in ( 18)).

Combining theorem 2 with

D 1 (N ) = • • • = D n (N ) = I N , (17)
, and ( 18), we thus get the following corollary to theorem 2:

Corollary 1 Assume that V has uniformly distributed phases. Then the limit moment

M n = lim N →∞ E tr L V H V n satsifies the inequality φ(a n 1 ) ≤ M n ≤ 1 c E(a n 2 ),
where a 1 has the distribution µ c of the Marchenko Pastur law, and a 2 has the Poisson distribution ν c . In particular, equality occurs for m = 1, 2, 3 and c = 1 (since all partitions are noncrossing for m = 1, 2, 3).

Corollary 1 thus states that the moments of Vandermonde matrices with uniformly distributed phases are bounded above and below by the moments of the classical and free Poisson distributions, respectively. The different Poisson distributions enter here because their (free and classical) cumulants resemble the c |ρ|-1 -entry in theorem 1, where we also can use that K ρ,u = 1 if and only if ρ is noncrossing to get a connection with the Marchenko Pastur law. To see how close the asymptotic Vandermonde moments are to these upper and lower bounds, the following corollary to theorem 3 contains the first moments:

Corollary 2 When c = 1, the limit moments The first three moments coincide for the three distributions, and are 1, 2, and 5, respectively.

M n = lim N →∞ E tr L V H V n ,
The numbers f p n and p n are simply the number of partitions in N C(n) and P(n), respectively. The number of partitions in N C(n) equals the Catalan number C n = 1 n+1 2n n [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF], so they are easily computed. The number of partitions of P(n) are also known as the Bell numbers B n [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF]. They can easily be computed from the recurrence relation

B n+1 = n k=0 B k n k .
It is not known whether the limiting distribution of our Vandermonde matrices has compact support. Corollary 2 does not help us in this respect, since the Marchenko Pastur law has compact support, and the classical Poisson distribution has not. In figure 4, the mean eigenvalue distribution of 640 samples of a 1600 × 1600 Vandermonde matrix with uniformly distributed phases is shown. While the Poisson distribution ν 1 is purely atomic and has masses at 0, 1, 2, and 3 which are e -1 , e -1 , e -1 /2, and e -1 /6 (the atoms consist of all integer multiples), the Vandermonde histogram shows a more continous eigenvalue ditribution, with the peaks which the Poisson distribution has at integer multiples clearly visible here as well (the peaks are not as sharp though). We remark that the support of V H V goes all the way up to N , but lies within [0, N ]. It is also unknown whether the peaks at integer multiples in the Vandermonde histogram grow to infinity as we let N → ∞. From the histogram, only the peak at 0 seems to be of atomic nature. In figures 5 and 6, the same histogram is shown for 1600×1200 (i.e. c = 0.75) and 1600 × 800 (i.e. c = 0.5) Vandermonde matrices, respectively. It should come as no surprise that the effect of decreasing c is stretching the eigenvalue density vertically, and compressing it horizontally. just as the case for the different Marchenko Pastur laws. Eigenvalue histograms for Gaussian matrices which in the limit give the corresponding (in the sense of corollary 1) Marchenko Pastur laws for figures 5 (i.e. µ 0.75 ) and 6 (i.e. µ 0.5 ), are shown in figures 7 and 8.

Deconvolution

Deconvolution with Vandermonde matrices (as stated in (6) in theorem 1) differs from the Gaussian deconvolution counterpart [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF] in the sense that there is no multiplicative [START_REF] Nica | Lectures on the Combinatorics of Free Probability[END_REF] structure involved, since K ρ,ω is not multiplicative in ρ. The Gaussian equivalent of theorem 3 (i.e. V H V replaced with 1 N XX H , with X an L × N complex, standard, Gaussian N XX H , with X an L × N = 1200 × 1600 complex, standard, Gaussian matrix. N XX H , with X an L × N = 800 × 1600 complex, standard, Gaussian matrix. matrix) is

m 1 = d 1 m 2 = d 2 + d 2 1 m 3 = d 3 + 3d 2 d 1 + d 3 1 m 4 = d 4 + 4d 3 d 1 + 3d 2 2 + 6d 2 d 2 1 + d 4 1 m 5 = d 5 + 5d 4 d 1 + 5d 3 d 2 + 10d 3 d 2 1 + 10d 2 2 d 1 + 10d 2 d 3 1 + d 5 1 m 6 = d 6 + 6d 5 d 1 + 6d 4 d 2 + 15d 4 d 2 1 + 3d 2 3 + 30d 3 d 2 d 1 + 20d 3 d 3 1 + 5d 3 2 + 10d 2 2 d 2 1 + 15d 2 d 4 1 + d 6 1 m 7 = d 7 + 7d 6 d 1 + 7d 5 d 2 + 21d 5 d 2 1 + 7d 4 d 3 + 42d 4 d 2 d 1 + 35d 4 d 3 1 + 21d 2 3 d 1 + 21d 3 d 2 2 + 105d 3 d 2 d 2 1 + 35d 3 d 4 1 + 35d 3 2 d 1 + 70d 2 2 d 3 1 + 21d 2 d 5 1 + d 7 1 , (19)
(where the m i and the d i are computed as in ( 9) by scaling the respective moments by c). This follows immediately from asymptotic freeness, and from the fact that 1 N XX H converges to the Marchenko Pastur law µ c . In particular, when all D i (N ) = I L and c = 1, we obtain the limit moments: 1,2,5,14,42,132,429, which also were listed in corollary 2. One can also write down a Gaussian equivalent to the fluctuations of Vandermonde matrices (13) (fluctuations of Gaussian matrices are handled more thoroughly in [START_REF] Mingo | Second order freeness and fluctuations of random matrices: I. Gaussian and Wishart matrices and cyclic Fock spaces[END_REF]). These are (20) These equations can be proved using the same combinatorical methods as in [START_REF] Ryan | Channel capacity estimation using free probability theory[END_REF]. Only the first equation is here stated as an exact expression. The second and third equations also have exact counterparts, but their computations are more involved. Similarly, one can write down a Gaussian equivalent to theorem 4 for the exact moments. For the first three moments (the fourth moment is dropped, since this is more involved), these are

m 1 = d 1 m 2 = d 2 + d 2 1 m 3 = 1 + N -2 d 3 + 3d 1 d 2 + d 3 1 .
This follows from a careful count of all possibilities after the matrices have been multiplied together (for this, see also [START_REF] Ryan | Channel capacity estimation using free probability theory[END_REF], where one can see that the restriction that the matrices D i (N ) are diagonal can be dropped in the Gaussian case). It is seen, contrary to theorem 4 for Vandermonde matrices, that the second exact moment equals the second asymptotic moment from [START_REF] Ryan | Free deconvolution for signal processing applications[END_REF], and also that the convergence is faster (i.e. O(n -2 )) for the third moment (this will also be the case for higher moments).

CONCLUSION AND FURTHER DIRECTIONS

We have shown how asymptotic moments of random Vandermonde matrices can be computed analytically, and treated many different cases. Vandermonde matrices with uniformly distributed phases proved to be the easiest case and was given separate treatment, and it was shown how the case with more general phases could be expressed in terms of the case of uniformly distributed phases. In addition to the general asymptotic expressions stated, exact expressions for the first moments of Vandermonde matrices with uniformly distributed phases were also stated.

Throughout the paper, we assumed that only diagonal matrices were involved in mixed moments of Vandermonde matrices. The case of non-diagonal matrices is harder to address, and should be addressed in future research. The analysis of the support of the eigenvalues is also of importance, as well as the behavior of the maximum and minimum eigenvalue. The methods presented in this paper can not be used directly to obtain explicit expressions for the asymptotic mean eigenvalue distribution, so this is also a case for future research. A way of attacking this problem could be to develop for Vandermonde matrices analytic counterparts to what one has in free probability (such as the Rand S-transform and their connection with the Stieltjes transform).

3. The 4 Fig. 1 .

 41 Fig.1. MSE of the first 4 estimated moments from the exact moments for 80 samples for varying matrix sizes, with N = L. Matrices are on the form V H V with V a Vandermonde matrix with uniformly distributed phases. The MSE of the first 4 exact moments from the asymptotic moments is also shown.

Fig. 2 .Fig. 3 .

 23 Fig.2. MSE of the first 4 moments from the actual moments for 320 samples for varying matrix sizes, with N = L. Matrices are on the form V H V with V a Vandermonde matrix with uniformly distributed phases. The MSE of the moments and the asymptotic moments is also shown.

  the moments f p n of the Marchenko Pastur law µ 1 , and the moments p n of the Poisson distribution ν 1 satisfy f p 4 = 14 ≤ M 4 = 44 3 ≈ 14.67 ≤ p 4 = 15 f p 5 = 42 ≤ M 5 = 146 3 ≈ 48.67 ≤ p 5 = 52 f p 6 = 132 ≤ M 6 = 3571 20 ≈ 178.55 ≤ p 6 = 203 f p 7 = 429 ≤ M 7 = 2141 3 ≈ 713.67 ≤ p 7 = 877.

Fig. 4 .

 4 Fig. 4. Histogram of the mean eigenvalue distribution of 640 samples of V H V, with V a 1600 × 1600 Vandermonde matrix with uniformly distributed phases.

Fig. 5 .

 5 Fig. 5. Histogram of the mean eigenvalue distribution of 640 samples of V H V, with V a 1600 × 1200 Vandermonde matrix with uniformly distributed phases.

Fig. 6 .

 6 Fig. 6. Histogram of the mean eigenvalue distribution of 640 samples of V H V, with V a 1600 × 800 Vandermonde matrix with uniformly distributed phases.

Fig. 7 .

 7 Fig. 7. Histogram of the mean eigenvalue distribution of 20 samples of 1N XX H , with X an L × N = 1200 × 1600 complex, standard, Gaussian matrix.

Fig. 8 .

 8 Fig. 8. Histogram of the mean eigenvalue distribution of 20 samples of 1N XX H , with X an L × N = 800 × 1600 complex, standard, Gaussian matrix.

E

  tr n D(N ) 1 N XX H 2 = (tr n (D(N )) 2 + 1 nN tr n (D(N ) 2 ) E tr n D(N ) 1 N XX H n = (tr n (D(N )) n + O(N -2 ) E tr n D(N ) 1 N XX H tr n D(N ) 1 N XX H 2 = tr n (D(N ))tr n (D(N ) 2 ) + O(N -2 ).

This work was supported by Alcatel-Lucent within the Alcatel-Lucent chair on flexible radio.