• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Channel Capacity Optimization Based on Riemannian Trust Region Algorithm in IRS-Aided Communication System

    2022-10-27 04:42:02JinzhiLiuDanWangJiaminLiangZhiqiangMei
    China Communications 2022年10期

    Jinzhi Liu,Dan Wang,Jiamin Liang,Zhiqiang Mei

    College of Communication and Information Engineering,Chongqing University of Posts and Telecommunications,Chongqing 400065,China

    *The corresponding author,email:81559234@qq.com

    Abstract:Intelligent Reflecting Surface(IRS)can offer unprecedented channel capacity gains since it can reconfigure the signal propagation environment.We decide to maximize the channel capacity by jointly optimizing the transmit-power-constrained precoding matrix at the base station and the unit-modulusconstrained phase shift vector at the IRS in IRSassisted multi-user downlink communication.We first convert the resulting non-convex problem into an equivalent problem,then use the alternate optimization algorithm.While fixing the phase shift vector,we can obtain the optimal precoding matrix directly by adopting standard optimization packages.While fixing the precoding matrix,we propose the Riemannian Trust-Region(RTR)algorithm to solve this optimization problem.And the key of the RTR algorithm is the solution of the trust-region sub-problem.We first adopt the accurate solution based on Newton’s(ASNT)method to solve this sub-problem,which can obtain the global solution but cannot guarantee that the solution is optimal since the initial iteration point is difficult to choose.Then,we propose the Improved-Polyline(IPL)method,which can avoid the difficulty of the ASNT method and improve convergence speed and calculation efficiency.The numerical results show that the RTR algorithm has more significant performance gains and faster convergence speed compared with the existing approaches.

    Keywords:intelligent reflecting surface(IRS);multiuser;MIMO;channel capacity;unit modulus constraint;Riemannian manifold optimization

    I.INTRODUCTION

    As the 5-th Generation(5G)wireless communication networks putting into commercialization,individuals are eager for faster data transmission speed.The key technologies of 5G,such as massive multiple-input multiple-output(MIMO)and millimeter wave(mm-Wave)communications,have made practical contributions to realize this function.However,additional high hardware costs,massive power/energy consumption,and the site selection of 5G base stations are the main hindrances to their implementation in practice[1–7].Therefore,to realize green and sustainable development in 5G and beyond wireless networks,research on finding both spectrum-efficient and energy-efficient techniques is vital for sustainable capacity growth[8–12].

    To solve the challenges mentioned above,intelligent reflecting surface(IRS)is considered a promising green,high cost-effectiveness,and energy-efficient technique in the 5G and next-generation mobile communication system.IRS is an artificial surface composed of electromagnetic materials[3].Specifically,IRS can smartly tune the wireless propagation environment by integrating numerous low-cost passive reflecting elements on the surface[4].Furthermore,IRS possesses appealing advantages such as low profile and light-weight,which provide high flexibility for actual implementation.For instance,IRS can be easily installed/removed on/from the wall,ceiling,advertising panel,and even clothing[13–17][18].Moreover,IRS can provide additional control signal paths by optimizing its reflection coefficient,which can effectively eliminate unexpected channel interference between users[19].Motivated by the above merits,IRS has a bright prospect in future wireless networks.

    Figure 1.IRS-aided multi-user MIMO communication system model.

    Recently,there are several works on IRS reflection optimization in IRS-assisted wireless communication systems[7–12].Reference[7]utilized Semi-Definite Relaxation(SDR)method and alternate optimization technology to maximize the received power in IRSassisted multi-user multi-input single-output(MISO)communication system.However,the SDR method only can obtain approximate solutions with high complexity.In[10],the author adopted the gradient descent method and Sequential Fraction Programming(SFP)while designing resource allocation to maximize energy or spectrum efficiency at the same time.However,the unit modulus constraints caused by IRS were solved by a series of approximate sub-problems in the SFP method,which may lead to some errors and performance loss.The author in[11]studied the resource allocation design for secure communication in IRS-assisted communication,and adopted continuous convex approximation to enhance the security of the physical layer.In[12],the author proposed a complex circle manifold(CCM)method based on the gradient descent algorithm to solve the unit modulus constraints.However,the CCM method is a Riemannian first-order algorithm.When the solution is closer to the target value,the step size will be smaller,and the convergence speed will be slower,which causes multiple iterations.In summary,no method can maintain excellent performance while keeping low complexity and fast convergence speed in the IRS-assisted communication system.

    Motivated by the above,we consider an IRS-aided multi-user MIMO communication system in a single cell,as shown in Figure 1 in this work.In this system,a multi-antenna base station(BS)serves multiantenna users with the help of IRS with discrete phase shifts[20,21].Therefore,it allows us to formulate the channel capacity maximization problem subject to the maximum transmit power constraint of BS and the unit modulus constraints of IRS.Specifically,we jointly optimize the transmit beamformer at the BS and the phase shifter at the IRS to maximize the channel capacity of all users.However,the formulated optimization problem is hard to address since optimization variables are highly coupled,and the unit modulus constraints are highly non-convex.

    The contributions of this work are summarized as follows:

    1)For the non-convex optimization problem,we convert the original problem into a more tractable problem by using the equivalent relationship between the channel capacity and the weighted minimum mean square error(WMMSE)[12].Then we use the alternate optimization algorithm to optimize the precoding matrix at the BS and the phase shifts vector of IRS one by one.

    2)The resultant optimization problem becomes a convex optimization problem when fixing the phase shift vector.We can regard it as a second-order cone programming(SOCP)problem due to the power constraints.Then the standard optimization package is used to obtain the optimal precoding matrix[22].Given the precoding matrix,the unit modulus constraints of IRS can form a product manifold and embed it into the search space[23],and the constraint problem is reformulated into an unconstrained problem.Finally,we propose the Riemannian trust-region(RTR)algorithm to obtain the local optimal solution.And the solution of the trust-region sub-problem is the key to the RTR algorithm.We first adopt the accurate solution based on Newton(ASNT)method to solve the trust-region sub-problem,which can obtain the global solution.However,this method is difficult to select the initial iteration point and this point must be near the optimal point to ensure convergence.Thus,the ASNT method cannot guarantee that the solution is the global optimal solution.Then,we propose the Improved-Polyline(IPL)method,which can avoid the difficulty of initial value selection in the ASNTmethod and improve the convergence speed and calculation efficiency.

    Table 1.Symbol description.

    3)Numerical results show that the RTR with the IPL method converges faster than the RTR with the ASNT method.Moreover,it is shown that,compared with the existing benchmark schemes,the RTR algorithm has excellent performance gains while maintaining low complexity and fast convergence speed.Furthermore,numerical results also show that the IRSassisted communication system should make the best tradeoff between the number of reflecting elements of the IRS and the transmit power of the BS.

    The symbol description is shown in Table 1 below.

    II.SYSTEM MODEL AND PROBLEM FORMULATION

    2.1 System Model

    As shown in Figure 1,we investigate an IRS-aided multi-user MIMO communication system model in this paper.The IRS,composed of N reflecting elements,assiststhe downlink communication between the BS with M1multi-antennas K and M2users with multi-antennas.We only consider the path reflected by the IRS for the first time,and other paths are ignored owing to the impact of path loss.In addition,we assume that all channels involved are quasi-static flatfading channels and the BS can perfectly know the channel state information(CSI)of them1.Specifically,the IRS configuration in Figure 1 is the semi-passive structure2[24],where additional sensing devices are interlaced with IRS reflecting elements.Thus,the CSI can be obtained by utilizing the deep learning and compressive sensing tools in[25].

    Then,the received baseband signalykat the userkcan be written as:

    where HAI,k∈CM1×Ndenotes the baseband channel matrix between the BS and IRS;HIU,k∈CM2×Ndenotes the baseband channel matrix between the IRS and user;HAU,k∈CM1×M2denotes the baseband channel matrix between the BS and userk;wk∈CN(0,σ2I)is the additional white Gaussian noise at userk;φ=diag(αφ1,···,αφN)∈CN×Nis a diagonal matrix at the IRS,whereφn=ejθn,n=1,2,···,N,α∈(0,1]is the amplitude coefficient andθnrepresents the phase shift induced by then-th reflecting element.xk=Gksk∈CM1×1denotes the transmitted signal for userkat the BS where Gk∈CM1×qand sk∈Cq×1represents the linear precoding matrix and the information-bearing symbol data at userkrespectively,and q is the data stream of each user.Besides,skis independent,which should satisfyandwhenk/=j.

    For ease of implementation,we consider a practical IRS model where we assume that the amplitude coefficientα=1,and the phase shiftθninduced by then-th reflecting element is discrete as in[3,13,26].Specifically,we letddenote the number of bits,then the number of phase shift levelsLps=2d,and assume thatθnis uniformly quantized in the interval[0,2π).Thus,the set of discrete phase-shift values at each reflecting element is given byθ∈F={0,Δθ,...,(Lps-1)Δθ},whereΔθ=2π/Lps[13].

    2.2 Problem Formulation

    Then,the channel capacity for userkis expressed as:

    where Hk=HIU,kφHHAI,k+HHAU,kis the combined channel matrix.We have known that the transmit signal of the userkisxk=Gksk.Then,the sum of the transmitted signal power of all users should be less than or equal to the transmit power threshold:

    where Pmax>0 is the maximum transmit power.

    Then,we have known thatφn=ejθnin the IRS diagonal matrixφand discrete phase shiftθnis hard to solve directly.Thus,we denote the collection of diagonal elements ofφbyζ=[ζ1,...,ζN]H,whereζn=ejθn,?n∈N.To make the targeted problem more tractable,the discrete phase shiftθnmay be relaxed to their continuous counterparts.Specifically,we assume that the phase-shift resolution is infinite,i.e.,2d?1,then directly quantize each of the obtained continuous phase shifts to its nearest discrete value inF[1,3].Thus,the discrete phase shift constraintθn∈Fcan be equivalent to the unit modulus constraint[13]as follows:

    In summary,the channel capacity maximum problem can be written as:

    whereωkdenotes weighted factor at the userk.Constraint(6)ensures that the BS transmit power of all users is kept below the maximum transmit power Pmax,while constraint(7)is the unit modulus constraint at the IRS.

    III.CHANNEL CAPACITY MAXIMIZATION

    3.1 Problem Reformulation

    From the problem(P0),we can easily observe that the precoding matrix and the phase shift vector are coupled,and the constraint(7)is highly non-convex,so the problem is NP-hard.

    To solve this problem(P0),we transform it into a more tractable way by exploiting the relationship between the mean square error(MSE)of the optimal decoding matrix and the channel capacity expression[14].For ease of implementation,we use the linear decoding matrix,and then the estimated signal vector?skis expressed as:

    where DHk∈Cq×M2denotes the linear decoding matrix for the user.We have known that the transmitted signal and the noise are independent.Thus thek-th user’s MSE matrix can be given by

    Then,by introducing a set of auxiliary matrices W=and setting a set of decoding matrices asD={Dk,?k∈K},we define the similar functionrk(W,D,G,ζ)as[15]:

    Thus,we can build the relationship between the problem(P0)and this function in(10)by exploiting the relationship between the channel capacity and the meansquare error as[14,15].And the problem(P0)can be reformulated as:

    We can observe that the problem(P1)introduces two optimization variables compared with the problem(P0).However,the problem(P1)is more tractable to solve.Since,for any given G orζ,when the auxiliary matrices W and the decoding matrices D are also fixed,rk(W,D,G,ζ)is a concave function for each set of optimization matrices.

    In the following,we decide to solve the problem(P1)by adopting the alternate optimization(AO)algorithm to optimize each set of optimization matrices alternately.

    We first consider the decoding matrix Dand the auxiliary matrixW by exploiting the AO algorithm since these two optimization variables are only related to the functionrk(W,D,G,ζ).Thus,we can derive the optimal solution for the decoding matrix Dkby fixing W,G,ζ.Specifically,we find the first-order derivative of the functionrk(W,D,G,ζ)with respect to Dkand make it equal to zero.The optimal decoding matrix can be given by:

    Similarly,by fixing W,G,D,the optimal auxiliary matrix Wkcan be expressed as:

    In the following,we concentrate on optimizing the precoding matrix and the phase shift vector.

    3.2 Optimization with Respect to the Precoding MatricesG

    In this subsection,we optimize the transmit precoding matrix Gwhile fixingW,D,ζ.Specifically,we substitute Ekinto the problem(P1)and ignore the constant terms.Thus,the precoding matrix optimization problem can be expressed as:

    whereAk=ωkHHkDkWkDHkHk.

    We can easily observe that the precoding matrix optimization in(14)is a convex optimization problem,and the constraint(15)is a second-order cone constraint.Therefore,this problem can be regarded as a second-order cone programming(SOCP)problem,and the existing standard optimization tool such as CVX[22]can be used to obtain the optimal solution of the transmit precoding matrix.

    3.3 Optimization with Respect to the Phase Shift Vectorζ

    Similarly,we substitute Ekinto the problem(P1)while fixingW,D,G.Then,by ignoring the constant terms,the phase shift vector optimization problem can be expressed as:

    Then,we substitute the combined channel matrix Hk=HIU,kφHHAI,k+HHAU,kinto(18),the optimization problem in(18)can be reformulated as follows:

    Note thatζis a collection vector of diagonal elements ofφ,then according to the matrix equation in[16],we have:

    Similarly,we denote diagonal elements of the matrix V as,and arrive that Tr.Then,the optimization problem in(20)can be equivalently rewritten as:

    where U=B⊙C.

    Obviously,the unit modulus constraint in(19)is still the main obstacle to solve this optimization problem.Aiming at this non-convex constraint,we decide to use the Riemannian manifold optimization tool to embed the constraint into the search space[17].Then,this problem can be an unconstrained optimization problem based on the search space.We have known that the CCM method proposed in[12]is a first-order Riemannian manifold optimization method,which mainly uses the Riemannian gradient descent method.This method can obtain the global solution,but it can’t guarantee global optimal solution and convergences slowly.Therefore,we propose a second-order Riemannian manifold optimization method—Riemannian trust-region(RTR)algorithm.Specifically,the RTR algorithm uses second-order geometric shapes in the selected trust-region to avoid saddle points and obtains more accurate result while having low complexity and fast convergence speed.What’s more,the solution of the trust-region sub-problem to obtain search direction is the key of the RTR algorithm.ASNT method is first used to solve this sub-problem[27].However,this method is hard to choose the initial iterative point,and the point must be near the optimal point to guarantee convergence.Therefore,we propose the IPL method based on the accurate solution method in[27]and the piecewise linear interpolation method in[28].

    Figure 2.Riemannian manifold geometric interpretation.

    3.3.1 RTR Algorithm

    First,the unit modulus constraint of the problem(22)can be regarded as the product of N complex circles.Then,we define a complex circle product manifoldMas shown in Figure 2:

    Before algorithms are proposed,we first provide background knowledge about manifold optimization[23]:

    1)Tangent space

    According to Figure 2,we assume thatζ(i)is the current iteration point.For the product manifoldM,the tangent spaceTζ(i)Mis of great significance in determining the search direction of the algorithm.For any pointζ(i)∈M,the tangent space is composed of all tangent vectors passing through the point[23].Then,the tangent spaceTζ(i)Mat the pointζ(i)∈Mis as follows:

    where z is the tangent vector at the pointζ(i);0Nis the N-dimensional zero vector.

    2)Euclidean gradient and Riemannian gradient

    In Euclidean space,the Euclidean gradient represents the direction in which the objective function changes fastest at that point[23].Then,at the pointζ(i)∈M,the Euclidean gradient Gradζ(i)fof the objective function in(22)is as follows:

    Similar to the Euclidean gradient,the Riemannian gradient is the tangent direction where the objective function grows fastest.And the Riemannian gradient gradf(ζ(i))of the objective function is defined as a unique element in the tangent space,which satisfies the following conditions[23]:

    3)Retraction

    As shown in Figure 2,moving along the tangent vector direction of the pointζ(i)∈Mto find the target point on the manifoldMis defined as retraction,which is the mapping from the tangent space to the manifold[23].Riemannian exponential retraction is the most commonly used contraction in Riemannian manifold optimization,but its computational cost is high.Therefore,we decide to use second-order retraction replacing exponential retraction,because it not only has lower computational cost,but also is the approximation of exponential mapping[27].Therefore,the retractionRζ(i)(ηi)at the pointζ(i)∈Mis defined as:

    4)Riemannian-Hessian

    The Riemannian-Hessian is the covariant derivative of the gradient vector field with respect to the Levi-Civita connectionζ(i).Specifically,The Riemannian-Hessian Hessf(ζ(i))at the pointζ(i)is a linear mapping from the tangent spaceTζ(i)Mto itself,which is expressed as:

    Equation(29)can be then calculated by the directional derivative of the Riemannian gradient and Euclidean gradient,as shown below:

    In the following,we proposed the RTR algorithm to solve the unconstraint optimization problem in(22).The RTR algorithm is similar to the trust-region algorithm in Euclidean space[28].The basic idea of the trust-region algorithm is as follows:First,the given trust region radius is as the upper bound of the displacement length.Then,a closed sphere called trust region is determined with the current iteration point as the center and the upper bound as the radius.Next,the key step of the RTR algorithm is to approximate the objective function with a simple function and construct it as a trust-region sub-problem.We decide to address this sub-problem by adopting the ASNT method or the IPL method to obtain the trial step(search direction).Finally,we use an evaluation function to decide whether to accept the trial step and determine the trust-region radius for the next iteration.

    Specifically,we first use retraction to promote the objective functionon the manifoldMto the cost functionon the tangent spa ceTζ(i)M.Since the tangent space is Euclidean spa ce,we can define a qua dratic model of(ζi),and use the ASNT method and IPL method on the Euclidean space to calculate the extreme value of the trust-region around the point 0ζi∈TζiM.The choice around the point 0ζiis based on the retraction,we have(0ζi)=f(ζi)[27].Then,the extreme value is retracted from the tangent space to the point on the manifold,and this point is used as a candidate for the next iteration.

    Currently,we first need to construct a trust-region sub-problem.

    In thei-th iteration,it is sufficient to use the quadratic approximation to define a second-order Taylor expansion of the objective functionf(ζ(i))on the manifold to define a quadratic model[23].Then,the trust-region sub-problem can be given by the quadratic model at the pointζ(i)∈M:

    whereηi∈Tζ(i)Mis the search direction of the current iteration;Δiis the trust-region radius of the current iteration.Then,whether to adjust the trust-region radiusΔiand update the next iteration pointζ(i+1)is determined by the following evaluation function[23]:

    We can see that the definition ofρiis the ratio of th(e actual de)cline of the objective functionfRζ(i)(ηi)to the predicted decline of the quadratic model functionmζ(i)(0)-mζ(i)(ηi)in thei-th iteration.Ifρi<1/4,this indicates that the quadratic model is bad and we need to reduce the trust-region radius and recalculate the search directionηi.Ifρi→1,this indicates that the quadratic model is good and the objective function has a good approximation,thus the trust-region radius should be expanded.

    The critical steps of the Riemann trust region(RTR)algorithm are summarized in Algorithm 1.

    Lemm(a 1.)[29]ζ(*i)is the m(inim)um of f(ζ(i)),(assu)methat fζ(i)∈C2,gradfζ(i)=0,Hessfζ(i)ispositive defnite,then the sequence generated by theRTR algorithm will converge to ζ*(i),with a secondorder convergence rate.

    Algorithm 1.RTR algorithm.Require:maximum trust-re trust-region radiusΔ0 threshold,initial iteratio gorithm iteration stop th gion radiusˉΔ>0,initial∈(0,ˉΔ),step acceptance n point ζ(0)∈M and alreshold ε,i=0.1:while■■gradf(ζ(i))■■2>ε do 2:Solve the trust-region sub-problem in(32)to obtain the search direction ηi 3:Calculate ρi in(34)4:Adjust the trust-region radius:min(2Δi,ˉΔ),ρζi>3 4and‖ηi‖=1 Δi+1=■■■■■■■■ ■■■■■■■4Δi,ρζi<1 4 Δi,others 5:Calculate the next iteration point ζ(i+1)6:if ρi>ρ′then 7: ζ(i+1)=Rζ(i)(ηi)8:else 9: ζ(i+1)=ζ(i)10:end if 11:i←i+1 12:end while 13:return current iteration point ζ(i)Δi

    Thus,we can observe that the RTR algorithm has second-order convergence locally.

    Recall that the solution of the trust-region subproblem in step 2 is the crucial difficulty of the RTR algorithm.In the following,we propose two efficient methods to solve it.

    3.3.2 ASNT Method

    ASNT is the accurate solution method based on Newton’s method.First,we introduce the following lemma proposed in[28]

    Lemma 2.η*i is the solution of the trust-region subproblem in(32),if and only if the inequality μ*i≥0exists,which makes:

    where Hessf(ζi)+μ*i I is a positive semi-defnite matrix.

    According to Lemma 1,we can observe that the main idea of the accurate solution method is to solve the following equations:

    Whenμi=0,the solution of trust-region sub-problem isη*i=-Hessf(ζ(i))-1gradf(ζ(i));whenμi>0,the solution of trust-region sub-problem is obtained by solving the equalityto getμ*iand substituting it into the first equation of the equations in(36).

    Therefore,whenμi>0,we define the function as follows:

    We can easily see that the above functions(μi)is a monotonically increasing concave function.Then we decide to use Newton’s method in[28]to find the approximate rootμ*iofs(μi)=0,and then we can get the solution of the trust-region sub-problem as follows:

    The key steps of ASNT method are summarized in Algorithm 2.When the algorithm converges,we can obtain the optimal solutionμ*iand substitute it into(38)to obtain the optimal solution of the trust-region subproblemη*i.

    Algorithm 2.ASNT method.Require:Hessf(ζ(i)),initial iteration point μmi=μ0,first-order derivative ψ(μi)=s′(μi),and iteration stop threshold ε.1:calculate η*i 2:for m=0,1,2,···do 3:μm+1),gradf(ζ(i)i =μmi-s(μmi)ψ(μmi)4:if μm+1i -μmi≤ε then 5: μ*i=μm+1i 6: calculate η*i in(38)7: break;8:else 9: m←m+1 10:end if 11:end for

    3.3.3 IPL Method

    As the number of IRS reflecting elements N increases,the problem scale becomes bigger.We have known that the ASNT method is not suitable for large-scale problems because the initial iteration point must be controlled near the best point to ensure convergence[28].Therefore,we propose the IPL method to handle the trust-region sub-problem,which can avoid the difficulty of initial iteration point selection in the ASNT method and improve the calculation efficiency and accuracy.

    We define the function by the optimal solutionη*iof the ASNT method in(38):

    whereμi≥0.

    Recall the constraint of the trust-region sub-problem in(33),we decide to solve the equationg(μi)-Δi=0 to get the optimal solution whenμi>0 and then obtain the optimal solution of this sub-problemη*i.However,g(μi)-Δi=0 is a complex non-linear equation,which is difficult to handle in general.Thus,we adopted the following steps to simplify it.

    First,we readily observe thatg(μi)is a continuous monotonic decreasing function.Therefore,for any given trust-region radiusΔi,ifΔi≥,theng(μi)<0 for anyμi≥0,thus we letμi=0 to obtain the optimal solutionη*i;ifΔi≤,we letμi=0 and haveg(0)-Δi>0,and there must be a real numberμMletg(μM)-Δi<0.According to the zero-point theorem in[30],the rooted interval ofg(μi)-Δi=0 is[0,μM].Then we divide the rooted interval[0,μM]into M equal parts,and adopt the linear interpolation method in[30],which can convert the original complex non-linear equationg(μi)-Δi=0 into a simple linear equation to obtain the root in each cell[μm-1,μm],m=1,...,M,and then getη*i.

    The key steps of the IPL method are summarized in Algorithm 3.

    3.4 Alternate Algorithm to Solve Problem(P1)

    Based on the above optimization analysis of the precoding matrix and the phase shift vector,the key steps to solve the problem(P1)are summarized in Algorithm 4.

    3.5 Complexity Analysis

    We now analyze the complexity of the AO algorithm.The complexity of calculating the decoding matrix D(i)in step 2 isO(KM32).The complexity of calculating the auxiliary matrix W(i)in step 3 isO(Kq3).And in step 4,the complexity of calculating the precoding matrix G(i+1)by using the second-order cone programming in CVX is given byO(KM31).When utilizing the RTR algorithm to obtain the phase shift vector in step 5,the complexity of calculating projection and retraction are both linear complexity[12],and the complexity mainly depends on the calculation of the Riemannian Hessian while solving the trust-region sub-problem by adopting the ASNT or IPL method,which is given by O(M1N2+N3).We use Titto represent the total number of iterations of the RTR algorithm,then the complexity of the RTR algorithm isO(Tit(N3+M1N2+N)).

    Therefore,the total complexity of the alternate optimization algorithm in Algorithm 4 is:

    Algorithm 3.IPL method.Require:Hessf(ζ(i)),Δi,gradf(ζ(i)).1:Select appropriate μM randomly,and satisfy g(μM)<Δi 2:Divide the rooted interval[0,μM]into M equal parts,and denote the step size as τ=μm M,m=0,1,...,M.Calculate the node μm=mτ and the function ym at this node,and mark the point(μm,ym)as the interpolation node.3:Utilize the linear interpolation method to construct M straight lines.Then the straight-line equation between two adjacent points μm-1 and μm is as follows:lm(μ)=ym-ym-1 μm-μm-1·μ,m=1,...,M 4:calculate η*i 5:for m=0,1,...,M do 6:if y0≤Δi then 7: η*i=-(Hessf(ζ(i)))-1gradf(ζ(i))8:else if ym≤Δi≤ym+1 then 9: use the straight line lm+1(μ)to replace g(μi)and solve the equation lm+1(μ)-Δi=0 to get μ*i 10: η*i=-(Hessf(ζ(i))+μ*i I)-1gradf(ζ(i))11:else 12: m←m+1 13:end if 14:end for 15:return the solution of the trust-region sub-Problem η*i

    In addition,the CCM algorithm proposed in[12],which is compared with the RTR algorithm proposed in this paper,has a computational complexity ofO(TCCMN2+N3),where TCCMrepresents the total number of iterations of the CCM algorithm.The numerical results in the next section will compare their convergence speed.

    IV.NUMERICAL RESULTS

    Algorithm 4.AO algorithm.Require:the current number of iterations i=0,the maximum number of iterations imax,the initial precoding matrix G(0)and the initial phase shift vector ζ(0),the error tolerance value ε,the value of the objective function in the problem(P0)1:Given G(i)and ζ(i),calculate the optimal decoding matrix D(i)in(14)2:Given G(i),ζ(i),D(i),calculate the optimal decoding matrix W(i)in(15)3:Given W(i),ζ(i)and D(i),use CVX tool to solve the problem(16)to obtain the optimal precoding matrix G(i+1)4:Given G(i),W(i)and D(i),calculate the optimal phase shift vector ζ(i+1)by solving problem(22)with RTR algorithm of ASNT method or IPL method 5:if|R(G(i+1),ζ(i+1))-R(G(i),ζ(i))|R(G(i),ζ(i)) <ε or i>imax then 6:break;7:else 8:i←i+1 9:return step1;10:end if

    In this section,we provide numerical results to prove that employing IRS can dramatically improve the channel capacity and the proposed algorithm has excellent performance.Refer to[7],we consider an IRSassisted downlink communication system where BS consists of M1=4 antennas,K=4 users consists of M2=2 antennas,the number of data streams is q=2,and the number of IRS reflecting elements is N=20.And we assume that all the involved channels obey the independent Rayleigh distribution and exhibit flat fading.The path loss model in dB is expressed as[12]:

    where PL0is the path loss at the reference distanced0=1mand set to-30dB;αis the path loss exponent and set to 3 in our simulation results.Except for some specific parameters specified later,other parameters are summarized in Table 2 below.

    According to Table 2,we have known the fixed positions of the BS and IRS,then assume that users are randomly distributed in a circle centered at(xu,0)with radius 10m,as shown in Figure 3.For ease ofimplementation,it is assumed that the user is located near the IRS.In step 2 of the RTR algorithm,if the ASNT method is used,the RTR algorithm is denoted as RTR-ASNT,and if the IPL method is used,the RTR algorithm is denoted as RTR-IPL.

    Table 2.Simulation parameters.

    Figure 3.IRS-aided MIMO multi-user scenario(top view).

    4.1 Convergence of the AO Algorithm

    We first investigate the convergence behavior of the AO algorithm.Figure 4 shows the channel capacity versus the number of iterations for the various number of reflecting elements,i.e.,for N=10,30,and 50.We can observe that as N increases,the convergence speed becomes slightly slower since more optimization variables are involved,and then more iterations are required for convergence.However,for different values of N,the proposed algorithms still converge within 150 iterations,which confirms the effectiveness of the proposed algorithms.Furthermore,both the RTR-ASNT and RTR-IPL algorithms in the AO algorithm are tested.We can also observe that the RTRASNT and RTR-IPL algorithm converge to almost the same value.Still,the RTR-IPL algorithm converges faster than the RTR-ASNT algorithm as the number of reflecting elements increases.This is expected because as N increases,the problem-scale becomes bigger,the ASNT method is not suitable for the largescale problem since the initial iteration point must be selected near the best point to ensure convergence,and the IPL method can avoid selecting the initial iteration point in the ASNT method.We can conclude that the optimization performance of the RTR-IPL algorithm is better than that of the RTR-ASNT algorithm.

    Figure 4.Convergence behavior of the AO algorithm.

    Therefore,in the following numerical results,we always use the IPL method to solve the trust-region subproblem when adopting the RTR algorithm.And we compare the following benchmark schemes with the RTR algorithm:

    Figure 5.Comparison of convergence behavior between the RTR and CCM algorithm.

    1)CCM algorithm:Use the CCM algorithm in[12]to solve the problem(22)at step 4 of the AO algorithm.

    2)Random phase shift(RandPhase):We assumed that the phase shift of each reflecting element is independently and randomly generated from the setθn∈F={0,Δθ,...,(Lps-1)Δθ}.Then in the AO algorithm,we can skip step 5 directly.

    3)Without-IRS:Set a system that does not employ IRS.

    4.2 Comparison of Convergence Behavior between the RTR and CCM Algorithm

    Figure 5 shows the convergence performance of the RTR algorithm at the first iteration of the AO algorithm when N=10 and the comparison of convergence behavior between the RTR and CCM algorithm.From Figure 5,we can see that under the previous settings,both the RTR algorithm and CCM algorithm can achieve convergence within 100 iterations.However,we can obviously observe that the convergence speed of the RTR algorithm enjoys a quadratic local convergence rate and converges significantly faster than that of the CCM algorithm,and the convergence channel capacity value of the RTR algorithms is larger than that of the CCM algorithm.We have known that the complexity of both is lower based on the complexity analysis in the previous section,thus we can conclude that the RTR algorithm possesses better optimization performance compared to the CCM algorithm from Figure 5.

    Figure 6.The relationship between the channel capacity and IRS reflecting elements.

    Figure 7.The relationship between the channel capacity and IRS’s location.

    4.3 Impact of IRS Reflecting Elements

    In Figure 6,we can see the relationship between the channel capacity values of all the schemes and the number of reflecting elements at the IRS.We can observe that,as the number of reflecting elements increases,the performance gains of the proposed algorithm become significant compared with the other three benchmark schemes.For example,compared with the Without-IRS case,the performance gain of the RTR algorithm is only 14.2%when N=20,but the performance gain is increased to 43.8% when N=70.We can see that as N increases,all With-IRS schemes have a higher channel capacity value than the Without-IRS scheme since IRS is passive and has no additional radio frequency(RF)chain and power consumption.

    Therefore,we can conclude that the deployment of IRS in a wireless communication system is practical and feasible.

    Figure 8.The relationship between the channel capacity and user’s location.

    4.4 Impact of IRS’s Location

    In Figure 7,we compare the influence of the location of the IRS fromxIRS=20m toxIRS=100m on the channel capacity.We can observe that the channel capacity values of all schemes decrease asxIRSincreases when 20

    4.5 Impact of User’s Location

    In Figure 8,we compare the relationship between the channel capacity and the user’s horizontal distancexu.It can be observed that,compared with the other three benchmark schemes,asxuincreases,when the user is far away from the BS and closer to the IRS,the performance gain of the RTR algorithm is significant.This is because users receive strong reflected signals from the IRS,which can reduce interference among users in the cell.

    Figure 9.The relationship between the channel capacity and maximum transmit power.

    4.6 Impact of Maximum Transmit Power

    The relationship between channel capacity performance and maximum transmit power is plotted in Figure 9.We can see that when the maximum transmit power is too small to be practically realized,i.e.,if the transmit power from the BS to the user is insufficient to meet the user’s rate requirement,then the performance of all schemes corresponds to very low values from Figure 9.Then,as the maximum transmit power value increases,the slope of the curve will increase obviously,and the performance of the RTR algorithm is still the best.Therefore,it can be concluded by combining the analysis of Figure 6 and Figure 8 that the IRS-assisted communication system should make the best tradeoff between the number of reflecting elements of the IRS and the transmit power of the BS.

    4.7 Impact of the Number of BS Antennas

    Recall to Figure 8,we have already known that the RTR algorithm provides a higher channel capacity value than the benchmark scheme Without IRS by varying the user’s horizontal distancexu.However,we have never considered the scenario where the number of BS antennas exceeds(or equals)that of total data streams.In Figure 10,we compare the relationship between the channel capacity and the number of BS antennas when using an IRS with N=20.We can observe that,as compared to the Without-IRS scheme,the gap in performance gains of the proposed algorithm is getting smaller when the number of BS antennas M1becomes larger.This is because the direct channel between the BS and the user becomes stronger as the antennas become larger.However,additional RF chains and power amplifiers need to be deployed while increasing the number of antenna elements,which may lead to huge energy consumption.Therefore,we suggest that the large-scale passive IRS can be introduced into the future wireless system.

    Figure 10.The impact of the number of BS antennas.

    Figure 11.The impact of weights.

    4.8 Impact of Weights

    As mentioned in problem(P0),weights are used to control fairness among users.We assume that the locations of the four users are in the coordinates(30,0),(50,0),(70,0),and(90,0)respectively in Figure 3.Among them,the first user is closer to the BS than the second user,and the fourth user is closer to the IRS than the third user.Two sets of weight values were tested:1)ωk=0.25,?k;2)ω1=0.15,ω2=0.35,ω3=0.3,ω4=0.2.It can be found from Figure 11 that,when the weights are equal,the data rate of user 1 is the highest since user 1 is closest to the BS;the data rate of user 4 is slightly lower than that of user 1and higher than the other two,since user 4 is far from the BS but closest to the IRS,the channel gain of user 4 increases.In addition,the data rate difference between the four users is not significant when the weights are not equal.Therefore,in order to ensure fairness between users,the user with low channel gain should be assigned a higher weight to make the user’s data rate distribution more balanced.

    V.CONCLUSION

    We aim to maximize the channel capacity in the IRSassisted multi-user communication system in this paper.First,the resulting optimization problem is a nonconvex NP-hard problem since the optimization variables are highly coupled.Therefore,we convert the original problem into an easy-to-solve form.Then we use the alternate optimization algorithm to optimize the precoding matrix and the phase shift vector alternately.Under the unit modulus constraint,we proposed the RTR-ASNT algorithm and RTR-ILP algorithm,where the constrained problem becomes an unconstrained problem.Then the local optimal solution is obtained.Numerical results show that the proposed algorithms achieve significant performance gains compared to the benchmark algorithm while maintaining low complexity and fast convergence speed.What’s more,we can conclude that the deployment of IRS in a wireless communication system is effective and feasible.

    ACKNOWLEDGEMENT

    This work was supported by the General Program of Natural Science Foudation of Chongqing Province of China(cstc2021jcyj-msxmX0454).

    NOTES

    1The proposed algorithms in this paper are also applicable in the imperfect CSI case while utilizing the stochastic successive convex approximation technique in[26]

    2In general,the acquirement of the CSI is a challenging problem.there are two main approaches for the IRS-involved channel acquisition,depending on whether the IRS elements possess any active RF chains or not,termed as semi-passive IRS and(fully)passive IRS,respectively[24].The first approach with active RF chains at the IRS,conventional channel estimation methods can be adopted for the IRS to estimate the channels of the BS-IRS and IRS-user links,respectively.In contrast,for the second approach without receive RF chains at the IRS,the IRS reflection patterns can be designed together with the uplink pilots to estimate the concatenated BS-IRS-user channels[31][32]

    丁香六月欧美| а√天堂www在线а√下载| 一个人观看的视频www高清免费观看| 可以在线观看的亚洲视频| 亚洲欧美精品综合久久99| 色av中文字幕| 国产精品乱码一区二三区的特点| 国产精品久久久久久久电影 | 欧美xxxx黑人xx丫x性爽| 国产69精品久久久久777片| 丁香六月欧美| 婷婷精品国产亚洲av在线| 少妇的逼好多水| 18禁裸乳无遮挡免费网站照片| 成人亚洲精品av一区二区| 亚洲无线观看免费| 无遮挡黄片免费观看| 亚洲av电影在线进入| 麻豆国产av国片精品| 亚洲不卡免费看| 成人高潮视频无遮挡免费网站| 亚洲乱码一区二区免费版| 久久香蕉国产精品| 老司机午夜十八禁免费视频| 久久精品国产亚洲av涩爱 | 听说在线观看完整版免费高清| 久9热在线精品视频| 亚洲人成网站高清观看| 一区二区三区高清视频在线| 欧美一区二区国产精品久久精品| 51国产日韩欧美| 成人一区二区视频在线观看| 日韩精品青青久久久久久| 免费高清视频大片| 国产美女午夜福利| 啦啦啦观看免费观看视频高清| 亚洲成人免费电影在线观看| 国产国拍精品亚洲av在线观看 | 国产精品1区2区在线观看.| 啦啦啦免费观看视频1| 日韩欧美免费精品| 熟女少妇亚洲综合色aaa.| 亚洲欧美日韩东京热| 日韩欧美精品v在线| 欧美性感艳星| 欧美日韩亚洲国产一区二区在线观看| 国产老妇女一区| 99久久精品热视频| 成年人黄色毛片网站| 男女那种视频在线观看| 91字幕亚洲| 欧美一区二区国产精品久久精品| 少妇的丰满在线观看| 国产成+人综合+亚洲专区| 亚洲美女黄片视频| 女人十人毛片免费观看3o分钟| 亚洲欧美精品综合久久99| 亚洲最大成人中文| 香蕉av资源在线| www.www免费av| 级片在线观看| 亚洲中文字幕日韩| 国产亚洲欧美在线一区二区| 夜夜夜夜夜久久久久| 国产亚洲精品一区二区www| 免费看a级黄色片| 欧美性感艳星| 精品日产1卡2卡| 少妇熟女aⅴ在线视频| 亚洲成人中文字幕在线播放| a在线观看视频网站| 中文在线观看免费www的网站| 少妇高潮的动态图| 97超视频在线观看视频| 精品一区二区三区视频在线 | 最近最新中文字幕大全免费视频| 日韩欧美在线二视频| 国产精品99久久99久久久不卡| 一区福利在线观看| 国产精品久久视频播放| 99久久成人亚洲精品观看| 国产精品亚洲一级av第二区| 国产三级黄色录像| 国产精品亚洲av一区麻豆| 熟妇人妻久久中文字幕3abv| 久久久久久久久大av| 高清毛片免费观看视频网站| 日韩大尺度精品在线看网址| 国产av在哪里看| 日韩欧美免费精品| 尤物成人国产欧美一区二区三区| 12—13女人毛片做爰片一| 国产三级黄色录像| 女人被狂操c到高潮| 伊人久久精品亚洲午夜| 国产精品久久久人人做人人爽| 高清毛片免费观看视频网站| 人妻丰满熟妇av一区二区三区| 国产免费一级a男人的天堂| 国产精品三级大全| 色综合欧美亚洲国产小说| 欧美激情久久久久久爽电影| 12—13女人毛片做爰片一| 亚洲精品美女久久久久99蜜臀| 高清毛片免费观看视频网站| 色综合婷婷激情| 村上凉子中文字幕在线| 狠狠狠狠99中文字幕| 亚洲久久久久久中文字幕| 99视频精品全部免费 在线| 五月玫瑰六月丁香| 性欧美人与动物交配| 757午夜福利合集在线观看| avwww免费| 99riav亚洲国产免费| 久久久久精品国产欧美久久久| 久久精品91蜜桃| 久久精品影院6| 精品无人区乱码1区二区| 久久精品亚洲精品国产色婷小说| 午夜福利欧美成人| 久久性视频一级片| 成人永久免费在线观看视频| 国产91精品成人一区二区三区| 青草久久国产| 麻豆一二三区av精品| 午夜影院日韩av| 老鸭窝网址在线观看| 久久人妻av系列| 一进一出抽搐gif免费好疼| 久久久国产精品麻豆| 美女高潮的动态| 日韩精品中文字幕看吧| 长腿黑丝高跟| 色综合欧美亚洲国产小说| 欧美乱色亚洲激情| 每晚都被弄得嗷嗷叫到高潮| 亚洲美女视频黄频| 亚洲国产欧美人成| 99riav亚洲国产免费| 亚洲人成伊人成综合网2020| 级片在线观看| 欧美+日韩+精品| 1000部很黄的大片| 国产精品av视频在线免费观看| 国产成人福利小说| 欧美成人a在线观看| aaaaa片日本免费| 亚洲人成网站高清观看| 在线十欧美十亚洲十日本专区| 国产野战对白在线观看| 欧美成人一区二区免费高清观看| 真人做人爱边吃奶动态| 日本 av在线| 1024手机看黄色片| 亚洲无线观看免费| 中文字幕人妻丝袜一区二区| 国产精品久久久久久久久免 | 欧美又色又爽又黄视频| 最近最新中文字幕大全免费视频| 欧美最黄视频在线播放免费| 国产主播在线观看一区二区| 免费在线观看亚洲国产| 日韩中文字幕欧美一区二区| 国产伦人伦偷精品视频| 成人一区二区视频在线观看| 草草在线视频免费看| 亚洲成av人片免费观看| 人人妻人人看人人澡| 国产一区在线观看成人免费| 亚洲人成电影免费在线| 色老头精品视频在线观看| 国产高清视频在线观看网站| 亚洲av不卡在线观看| 亚洲人成网站在线播放欧美日韩| 欧美日韩乱码在线| 免费在线观看日本一区| avwww免费| 麻豆国产av国片精品| 内射极品少妇av片p| 国产91精品成人一区二区三区| 少妇人妻精品综合一区二区 | 国产乱人视频| 久久久久国内视频| 国产亚洲精品av在线| 一a级毛片在线观看| 亚洲最大成人手机在线| 少妇熟女aⅴ在线视频| 中文字幕人妻丝袜一区二区| 国产精品久久久久久久久免 | 亚洲自拍偷在线| 欧美成人一区二区免费高清观看| 好看av亚洲va欧美ⅴa在| 久久久久久久久大av| 99热6这里只有精品| 国产精品一区二区三区四区免费观看 | 亚洲国产精品成人综合色| 亚洲av成人精品一区久久| 亚洲乱码一区二区免费版| 日韩高清综合在线| 在线国产一区二区在线| 欧美黑人欧美精品刺激| 99久久综合精品五月天人人| 国内少妇人妻偷人精品xxx网站| 狠狠狠狠99中文字幕| 亚洲欧美日韩高清在线视频| 黄色女人牲交| 欧美中文日本在线观看视频| 老熟妇仑乱视频hdxx| 成人午夜高清在线视频| 国产亚洲精品av在线| 国产私拍福利视频在线观看| 中出人妻视频一区二区| 波多野结衣高清无吗| 嫩草影院入口| 在线播放无遮挡| 色综合亚洲欧美另类图片| 很黄的视频免费| 小说图片视频综合网站| 欧美成人性av电影在线观看| 亚洲国产欧美网| 日韩欧美一区二区三区在线观看| 国产单亲对白刺激| 亚洲成人久久性| 久久6这里有精品| 国产欧美日韩精品一区二区| 久久中文看片网| 桃色一区二区三区在线观看| 日本黄色视频三级网站网址| 免费观看人在逋| 尤物成人国产欧美一区二区三区| 久久精品国产自在天天线| 久久久国产成人免费| 色老头精品视频在线观看| 国产一区二区亚洲精品在线观看| 精品久久久久久成人av| 国产午夜精品久久久久久一区二区三区 | 国内久久婷婷六月综合欲色啪| 大型黄色视频在线免费观看| 午夜免费激情av| 国产精品99久久久久久久久| 亚洲无线在线观看| 国产精品一区二区免费欧美| 美女免费视频网站| 99久久精品国产亚洲精品| 精品人妻一区二区三区麻豆 | 欧美xxxx黑人xx丫x性爽| 在线观看舔阴道视频| 国内揄拍国产精品人妻在线| 久久草成人影院| 久久久久久九九精品二区国产| 亚洲精华国产精华精| 香蕉久久夜色| 久久香蕉国产精品| 亚洲在线观看片| 日本 av在线| 最新在线观看一区二区三区| 91字幕亚洲| 日韩欧美精品免费久久 | 亚洲欧美日韩无卡精品| 精品无人区乱码1区二区| 欧美三级亚洲精品| 欧美激情在线99| 99久久99久久久精品蜜桃| 99热只有精品国产| 精品人妻一区二区三区麻豆 | 国产视频内射| 身体一侧抽搐| 午夜免费激情av| 欧美激情久久久久久爽电影| 亚洲性夜色夜夜综合| 两个人看的免费小视频| 国产精品一及| 亚洲第一电影网av| 亚洲美女视频黄频| 精品久久久久久成人av| eeuss影院久久| 久9热在线精品视频| 三级男女做爰猛烈吃奶摸视频| 我要搜黄色片| 亚洲一区二区三区不卡视频| 亚洲狠狠婷婷综合久久图片| 国产激情偷乱视频一区二区| 19禁男女啪啪无遮挡网站| 尤物成人国产欧美一区二区三区| 婷婷六月久久综合丁香| 特大巨黑吊av在线直播| 婷婷精品国产亚洲av在线| 久9热在线精品视频| 麻豆成人av在线观看| 少妇裸体淫交视频免费看高清| 国产av麻豆久久久久久久| 国产91精品成人一区二区三区| 少妇的丰满在线观看| 波多野结衣高清无吗| 99在线视频只有这里精品首页| 欧美区成人在线视频| 久久久久久久精品吃奶| 国产精品国产高清国产av| av在线天堂中文字幕| 成人欧美大片| 成年女人永久免费观看视频| 国产探花在线观看一区二区| 亚洲18禁久久av| 在线播放国产精品三级| 丰满的人妻完整版| 丰满人妻熟妇乱又伦精品不卡| 成人三级黄色视频| 国产单亲对白刺激| 97人妻精品一区二区三区麻豆| x7x7x7水蜜桃| 久久久久亚洲av毛片大全| 成人性生交大片免费视频hd| 悠悠久久av| 男插女下体视频免费在线播放| 亚洲av二区三区四区| 性色av乱码一区二区三区2| 亚洲中文日韩欧美视频| 午夜福利免费观看在线| 亚洲国产欧美网| 国产一区二区亚洲精品在线观看| 麻豆一二三区av精品| 免费高清视频大片| 欧美最黄视频在线播放免费| 国产一区二区在线av高清观看| 男女那种视频在线观看| 欧美黑人欧美精品刺激| 国内精品久久久久精免费| 高清日韩中文字幕在线| 精品久久久久久久久久久久久| 一a级毛片在线观看| 亚洲第一欧美日韩一区二区三区| www日本黄色视频网| 美女高潮的动态| 精品人妻1区二区| 国产淫片久久久久久久久 | 制服人妻中文乱码| 人妻夜夜爽99麻豆av| 制服丝袜大香蕉在线| 欧美日本视频| 欧美性猛交╳xxx乱大交人| 久久久久久久久中文| 有码 亚洲区| 亚洲欧美日韩高清在线视频| 九色成人免费人妻av| 久久国产精品人妻蜜桃| 久久久久久久久久黄片| 日韩大尺度精品在线看网址| 国产一区二区在线av高清观看| 看免费av毛片| 床上黄色一级片| 免费人成在线观看视频色| 波多野结衣高清作品| 男女那种视频在线观看| 长腿黑丝高跟| 一级毛片高清免费大全| 久久久久久大精品| 在线观看舔阴道视频| 中文字幕熟女人妻在线| 欧美日韩国产亚洲二区| 最近视频中文字幕2019在线8| 国产激情偷乱视频一区二区| 噜噜噜噜噜久久久久久91| 久久久久久久久中文| 成人特级av手机在线观看| 成人精品一区二区免费| 最近视频中文字幕2019在线8| 9191精品国产免费久久| 我要搜黄色片| 内射极品少妇av片p| 一个人免费在线观看电影| 国产午夜精品论理片| 国内毛片毛片毛片毛片毛片| 日韩国内少妇激情av| а√天堂www在线а√下载| 亚洲欧美日韩东京热| 两个人视频免费观看高清| 亚洲国产精品合色在线| 午夜福利在线观看吧| 在线a可以看的网站| www.www免费av| 在线天堂最新版资源| 午夜免费激情av| 亚洲 欧美 日韩 在线 免费| 99在线视频只有这里精品首页| 国产高潮美女av| 操出白浆在线播放| 免费看a级黄色片| 成人鲁丝片一二三区免费| 久久人人精品亚洲av| 国产成人欧美在线观看| 精品国产三级普通话版| 真实男女啪啪啪动态图| 国产国拍精品亚洲av在线观看 | 黄片小视频在线播放| 九色国产91popny在线| 亚洲精品色激情综合| 久久久精品大字幕| 九九热线精品视视频播放| 亚洲av中文字字幕乱码综合| 久久久精品欧美日韩精品| 香蕉久久夜色| 免费搜索国产男女视频| 熟女电影av网| av女优亚洲男人天堂| 成人国产一区最新在线观看| 成人鲁丝片一二三区免费| 18禁黄网站禁片免费观看直播| 亚洲avbb在线观看| 色吧在线观看| 免费一级毛片在线播放高清视频| 国产精品精品国产色婷婷| 久久精品亚洲精品国产色婷小说| 啪啪无遮挡十八禁网站| 夜夜看夜夜爽夜夜摸| 久久久色成人| av天堂在线播放| 久久久国产成人精品二区| 亚洲成人中文字幕在线播放| 精品久久久久久久久久久久久| 日韩欧美三级三区| 麻豆一二三区av精品| 乱人视频在线观看| 老熟妇乱子伦视频在线观看| a在线观看视频网站| 国产免费一级a男人的天堂| 亚洲人与动物交配视频| 国产真实伦视频高清在线观看 | 国产精品亚洲av一区麻豆| 少妇丰满av| 好看av亚洲va欧美ⅴa在| 一个人免费在线观看电影| 日本熟妇午夜| 久久久久国内视频| 国产精品综合久久久久久久免费| 桃色一区二区三区在线观看| 搡女人真爽免费视频火全软件 | 久久婷婷人人爽人人干人人爱| 九色国产91popny在线| 国产黄a三级三级三级人| 色综合婷婷激情| 性色av乱码一区二区三区2| 老司机在亚洲福利影院| 一本久久中文字幕| 国产免费一级a男人的天堂| 99热6这里只有精品| 美女免费视频网站| 观看免费一级毛片| a级毛片a级免费在线| 老司机深夜福利视频在线观看| 欧美黑人巨大hd| 动漫黄色视频在线观看| 三级毛片av免费| 成人三级黄色视频| 国产又黄又爽又无遮挡在线| 亚洲av日韩精品久久久久久密| 99热6这里只有精品| 午夜福利18| 午夜福利成人在线免费观看| 99热只有精品国产| 色综合婷婷激情| 亚洲七黄色美女视频| 欧美最新免费一区二区三区 | 女人被狂操c到高潮| 午夜a级毛片| 真实男女啪啪啪动态图| or卡值多少钱| 国产黄片美女视频| 久久久国产成人免费| 亚洲乱码一区二区免费版| 亚洲精品乱码久久久v下载方式 | 伊人久久精品亚洲午夜| 十八禁人妻一区二区| ponron亚洲| 免费看日本二区| 亚洲激情在线av| av专区在线播放| 男女之事视频高清在线观看| 成人特级黄色片久久久久久久| 女警被强在线播放| 啦啦啦免费观看视频1| 国产 一区 欧美 日韩| 国产黄a三级三级三级人| 又黄又粗又硬又大视频| 十八禁人妻一区二区| 久久精品人妻少妇| 老司机福利观看| 美女高潮喷水抽搐中文字幕| 亚洲成av人片在线播放无| 亚洲黑人精品在线| 一进一出抽搐gif免费好疼| 久久精品国产99精品国产亚洲性色| 黄色女人牲交| 一级毛片高清免费大全| 99久久无色码亚洲精品果冻| 91久久精品电影网| 成熟少妇高潮喷水视频| 久久精品影院6| 欧美绝顶高潮抽搐喷水| 毛片女人毛片| 国产真实伦视频高清在线观看 | 精品99又大又爽又粗少妇毛片 | 色吧在线观看| 中文资源天堂在线| 精品久久久久久久末码| 亚洲精品亚洲一区二区| 久久伊人香网站| 哪里可以看免费的av片| 欧美激情久久久久久爽电影| 国产淫片久久久久久久久 | 日本精品一区二区三区蜜桃| 国产伦一二天堂av在线观看| 亚洲五月婷婷丁香| 亚洲美女黄片视频| 国产精品一区二区免费欧美| 制服丝袜大香蕉在线| 夜夜看夜夜爽夜夜摸| 亚洲成a人片在线一区二区| 精品国产美女av久久久久小说| 一区二区三区国产精品乱码| 老鸭窝网址在线观看| 日本免费一区二区三区高清不卡| 99热这里只有精品一区| 国产一区二区在线av高清观看| 久久久久国产精品人妻aⅴ院| 国产精品嫩草影院av在线观看 | 色老头精品视频在线观看| 欧美色视频一区免费| 五月伊人婷婷丁香| 99精品在免费线老司机午夜| 欧美成人性av电影在线观看| 成人av一区二区三区在线看| 欧美极品一区二区三区四区| 国产成人福利小说| 精品福利观看| 一进一出抽搐动态| 国产伦精品一区二区三区视频9 | 亚洲18禁久久av| 内射极品少妇av片p| 成年人黄色毛片网站| 亚洲精品日韩av片在线观看 | 亚洲激情在线av| 亚洲国产日韩欧美精品在线观看 | 中文字幕av成人在线电影| 国产精品久久电影中文字幕| 成人特级av手机在线观看| 欧美最新免费一区二区三区 | www.熟女人妻精品国产| 精品久久久久久久毛片微露脸| 欧美日韩综合久久久久久 | 免费看美女性在线毛片视频| 亚洲精品在线美女| 久久久国产成人精品二区| 天堂动漫精品| 成人国产综合亚洲| 热99re8久久精品国产| 国产精品久久久人人做人人爽| netflix在线观看网站| 久久久久国内视频| 亚洲自拍偷在线| 亚洲欧美日韩高清在线视频| 久久久久久久精品吃奶| 国产三级中文精品| 国产aⅴ精品一区二区三区波| 听说在线观看完整版免费高清| 成人特级黄色片久久久久久久| 又爽又黄无遮挡网站| 长腿黑丝高跟| 亚洲av二区三区四区| 免费在线观看成人毛片| 亚洲av美国av| 国产精品国产高清国产av| 黄色片一级片一级黄色片| 日韩欧美一区二区三区在线观看| 国产熟女xx| 久久性视频一级片| eeuss影院久久| 成人精品一区二区免费| a在线观看视频网站| 亚洲国产精品999在线| 99热精品在线国产| 久久天躁狠狠躁夜夜2o2o| 国产男靠女视频免费网站| 亚洲在线观看片| 色播亚洲综合网| av天堂在线播放| 欧美一级a爱片免费观看看| 岛国在线免费视频观看| 大型黄色视频在线免费观看| 欧美大码av| а√天堂www在线а√下载| 中文字幕av成人在线电影| 国产aⅴ精品一区二区三区波| 精品免费久久久久久久清纯| 精品人妻偷拍中文字幕| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 99久久无色码亚洲精品果冻| 精品人妻1区二区| 高潮久久久久久久久久久不卡| 夜夜爽天天搞| 在线观看av片永久免费下载| 亚洲一区高清亚洲精品| 国产伦精品一区二区三区视频9 | av天堂在线播放| 色综合欧美亚洲国产小说| 色吧在线观看| 日日摸夜夜添夜夜添小说| 舔av片在线| 亚洲av二区三区四区| 亚洲国产色片| 国产午夜精品久久久久久一区二区三区 | 午夜福利在线观看吧| 国产高清视频在线播放一区| 色视频www国产| 丰满的人妻完整版|