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

    AWK-TIS:An Improved AK-IS Based on Whale Optimization Algorithm and Truncated Importance Sampling for Reliability Analysis

    2023-02-26 10:18:06QiangQinXiaoleiCaoandShengpengZhang

    Qiang Qin,Xiaolei Cao and Shengpeng Zhang

    1Aerospace Science and Industry Corporation Defense Technology Research and Test Center,Beijing,100854,China

    2Key Laboratory for Thermal Science and Power Engineering of Ministry of Education,Department of Engineering Mechanics,Tsinghua University,Beijing,100084,China

    ABSTRACT In this work,an improved active kriging method based on the AK-IS and truncated importance sampling(TIS)method is proposed to efficiently evaluate structural reliability.The novel method called AWK-TIS is inspired by AK-IS and RBF-GA previously published in the literature.The innovation of the AWK-TIS is that TIS is adopted to lessen the sample pool size significantly,and the whale optimization algorithm(WOA)is employed to acquire the optimal Kriging model and the most probable point(MPP).To verify the performance of the AWK-TIS method for structural reliability,four numerical cases which are utilized as benchmarks in literature and one real engineering problem about a jet van manipulate mechanism are tested.The results indicate the accuracy and efficiency of the proposed method.

    KEYWORDS Structural reliability;active kriging;whale optimization algorithm;AK-IS

    Nomenclature

    AK Active kriging

    TIS Truncated importance sampling

    WOA Whale optimization algorithm

    MPP Most probable point

    MCS Monte Carlo simulation

    FORM First order reliability method

    SORM Second order reliability method

    LSF Limit state function

    RBF Radial basis function

    DOE Design of experiment

    ξ Regression coefficient vector

    z(x) Stationary Gaussian process

    σz2Process variance

    R(xi,xj) Correlation function

    θ?Optimal correlation parameter

    PfFailure probability

    δPfCoefficient of variation of the failure probability

    λPositive penalty coefficient

    ζSmall constant

    χThreshold

    x(t)ii-thsearch agent in thet-thiteration

    The best agent in thet-thiteration

    AiRandom value

    BiRandom value

    CiRandom value

    NISNumber of importance sampling

    βReliability index

    1 Introduction

    Structural reliability analysis is often carried out to quantitatively assess the probability of structures’safe or failure state under the impact of uncertain factors in the environmental,structural and load parameters.Various structural reliability analysis methods,which can be broadly categorized into three main types:simulation-based methods,analytical methods or approximation methods,and meta-modeling methods,have been proposed over the past three decades[1,2].

    Concerning simulation-based methods,the crude Monte Carlo simulation(MCS)method is the basic form of this classification,which evaluates the failure probability by dividing the number of samples falling into the failure domain by all the samples generated from the limit state evaluations.Although MCS is the most convenient implementation approach,it is extremely time-consuming when dealing with small failure probability and complex numerical models(e.g.,fidelity finite element models(FEM)).Meanwhile,tremendous variation techniques of the MCS such as subset simulation[3,4],importance sampling [5,6] and directional sampling [7,8] have been developed to obtain more accurate results with fewer samplings.However,the main drawback of insufficient still exists in those methods.

    As analytical methods,the traditional first-order and the second-order reliability methods(FORM and SORM)[9–12]respectively approximate the limit state function(LSF)around the most probable failure point(MPP)with first-order and incomplete second-order functions.The derivatives of the LSF with respect to the random variables have thus to be evaluated.Hence,these methods are not suitable for implicit limit state function and high nonlinearity cases.

    Meta-modeling methods have been developed to balance the computational accuracy and efficiency during the structural reliability analysis,such as response surfaces[13–15],polynomial chaos expansion [16,17],artificial neural networks [18,19],support vector machine [20,21],radial basis function(RBF)[22,23]and Kriging metamodel[24,25].In recent years,Kriging metamodel methods have obtained remarkable attention in the state-of-the-art literature.From the design of experiment(DOE) aspect of view,the strategy for constructing a Kriging model can broadly be classified into two sorts,the one is “one-shot”and the other is adaptive sampling or sequential sampling methods[26].Literally,“one-shot”means producing enough design points before generating a Kriging model without supplementing new sample points in succeeding courses.Nevertheless,the adaptive sampling methods generate a relatively small number of design points at the initial step,and then add one or more samples in each course of iteration according to certain regulations [27].Furthermore,the exclusive feature of Kriging that it can predict the local uncertainty along with the prediction value promotes the development of the adaptive Kriging[24,28].

    Two typical algorithms,the efficient global reliability analysis(EGRA)[28]and adaptive Kriging(AK) with MCS (AK-MCS) [29],are generously adopted and depict outstanding characters in diverse structural reliability problems.The amount of research has then further ameliorated the performance of active Kriging-based methods in recent years.These ameliorations resort to different and more efficient learning functions,sampling schemes and stopping criteria.Learning functions include the expected feasible function(EFF)[28],the learning functionU[29],the learning functionH[30],expected risk function (ERF) [31],least improvement function (LIF) [32],reliability-based expected improvement function(REIF)[33],Folded Normal based Expected Improvement Function(FNEIF) [34],K-means and weighted K-medoids algorithms-assisted learning process [35],and PDEM-oriented expected improvement function (PEIF) [36].With regard to stopping criteria,it is common to see that when the values of learning functions reach certain thresholds the learning process stopped.For sampling tactics,importance sampling,subset simulation,line sampling,etc.have been adopted to evaluate failure probability in active Kriging(AK)reliability analysis,and some AK-based sampling methods were presented,such as AK-IS[37],AK-SS[38],AK-DS[39],AKOIS[40],AK-ARBIS [41],AK-SESC [42] and AKSE [43].By combining these sampling methods with learning functions,the number of calls of limit state functions has reduced dramatically with acceptable accuracy of variation coefficient of failure probability.Despite the aforementioned learning functions,the correlation parameterθalso plays a vital role in constructing an accurate Kriging model.Thus,the pattern search algorithm is applied in DACE to calculate the optimum ofθ[44].However,as a local optimization algorithm,pattern search is sensitive to the initial value and may get trapped in local optimum.Then,various global optimization algorithms,such as genetic algorithm(GA)[33,45],DIRECT algorithm[46],particle swarm optimization(PSO)[47],and artificial bee colony(ABC)[48],are adopted to find the optimal value ofθby solving the maximum likelihood equation.Compared with GA,PSO,and the pattern search,the whale optimization algorithm(WOA)[49]proposed in the year of 2016 has its own advantages,such as simple principle,easy implementation,high convergence speed and calculation accuracy.

    This paper aims to improve the efficiency and accuracy of the AK-IS for structural reliability analysis.In this paper,AWK-TIS is proposed,which is inspired by AK-IS [34] and RBF-GA [50],coupling with WOA and Truncated importance sampling (TIS) [51,52],and the WOA in AWK-TIS is adopted to seek the optimum of correlation parameterθand the MPP in different phases of each course of iterations.Moreover,inside of the optimalβ-sphere treated as the safe state,the samples are eliminated by TIS.Therefore,the prediction of Kriging model and the selection of the best next sample are only implemented through the candidate samples outside the optimalβ-sphere,which significantly reduces the computational time and saves computer memory in comparison with the original AK-IS.

    The remainder of the manuscript is organized as follows.Section 2 briefly introduces the Kriging method,AK-IS method,TIS and WOA.In Section 3,the procedures and the main steps of the AWKTIS method are proposed.In Section 4,the applicability of the proposed technique is validated by four numerical cases and one real engineering problem of jet van manipulation mechanism.Finally,conclusions are summarized in Section 5.

    2 Backgrounds

    Before introducing the proposed AWK-TIS algorithm,some basic theories are briefly recalled,e.g.,the Kriging method,the main steps of the AK-IS algorithm and truncated importance sampling(TIS)method,and the WOA.

    2.1 Basic Theory about Kriging Method

    Kriging metamodel,which consists of a parametric linear regression model and a nonparametric stochastic process,is an interpolation technique based on statistical theory.It needs a DOEs to determine its stochastic parameters and then predictions of the response can be inferred on any unknown point.Give a set of initial DOEsX=[x1,x2,...,xm],withxi?Rn(i=1,2,...,m) theith experiment,andG=[G(x1),G(x2),...,G(xm)]is the corresponding response toX.The approximate relationship between any experimentxand the responseG(x)can be denoted as[41,42]:

    whereξT=[ξ1,...,ξp]is a regression coefficient vector,fT(x)=is the basic functions vector which makes a global simulation in design space.In the ordinary Kriging,F(ξ,x)is a scalar and always taken asF(ξ,x)=ξ;z(x)is a stationary Gaussian process with zero mean and covariance between two data pointsxiandxjis denoted as:

    whereσz2is the process variance,andxi,xjare data points from the whole samplesX.Ris the correlation function aboutxiandxjwith a correlation parameter vectorθthat has to be estimated by pattern search method in DACE.A diversity of models can be adopted to define correlation functionsuch as Gaussian correlative model,experimental and linear model,etc.The widely employed Gaussian model is accepted in the paper and can be defined as:

    wherenis the dimension of the sample,xki,andθkare thekth components of data samplexi,xjandθ,respectively.

    Define correlation matrixR=[R(xi,xj)]m×m,Fis am×1 unit vector,thenξandσz2are given by:

    whereR,andare functions ofθ.

    Then at an unknown pointx,the Best Linear Unbiased Predictor(BLUP)of the response(x)is shown to be a Gaussian random variant(x)~N((x),(x)),where

    wherer(x)=[R(x,x1),R(x,x2),...,R(x,xm)],v(x)=FTR?1r(x)?1.(x)is usually taken as the estimated(x)at pointx.That means Kriging is an exact interpolation method.

    The correlation parameterθcan be obtained through the maximum likelihood estimation(MLE)[47,48]:

    where |R| is the determinant ofR.The optimal Kriging model can be obtained if the optimal correlation parameterθ?is guaranteed.The Eq.(8)can be solved by pattern search method,genetic algorithm,PSO,etc.It should be noted that since the introduction of WOA,this algorithm has been used in many fields and showed its superiority over many other metaheuristic algorithms such as PSO,ABC,GA,etc.As a result,WOA is adopted to searchθ?in this study.

    2.2 AK-IS Algorithm

    AK-IS mainly consists of two steps.The first step is the FORM approximation,which is utilized to find the MPP that is not accurate enough to evaluate the probability.The second step is to establish an AK model with a special learning function and stopping criteria to reduce the calls of original limit state function.In this step,IS method is adopted to calculate the failure probability and its coefficient of variation.

    IS requires the definition of a joint PDF?nthat is taken as a standard Gausssian one centered at the MPP found in the first step for the new sampling.The failure probability is estimated as[37]:

    whereIF(u)is indicator function,IF(u))={0 ifG(u)>0 and 1 ifG(u)≤0}.

    A number ofNISsamples denotedare generated according to importance sampling density function?n.The integral is then expressed by:

    The variance of the failure probability estimator is given as:

    The coefficient of variationδISof the failure probability estimator is then:

    Besides,the learning functionUproposed in AK-MCS is still used in AK-IS:

    whereμ(x)andσ(x)are Kriging mean and standard deviation,respectively.Moreover,the sample that has the minimum ofU(x)is called the best next sample and it is iteratively enriched the DOE to update the Kriging model.The learning process stops if the minimum of the learning functionUis larger than its threshold.Besides,the active procedure continues andNISnew samples are enlarged until the coefficient of variationδISis smaller than a limit value which is set to 0.05 in AK-IS.

    2.3 Truncated Importance Sampling

    Truncated importance sampling(TIS)method is developed by introducing a hypersphere with the radial,which is indicated byβshown in Fig.1,equals to the minimum distance from the origin to the limit state surface inu-space[51,52].

    Figure 1:Sketch of TIS

    The red point in Fig.1 indicates the MPP.

    The indicator function of the outer of theβ-sphere is defined as follows:

    where the inner of theβ-sphere contains part of the safe domain or even all of the district,while other domains including all the failure domain and part of the safe domain are out of theβ-sphere.It can be seen that the indicator functionIFof the samples located in theβ-sphere is=0,and there is unnecessary to call the true model to evaluate whether safe or not.The integral is then estimated by:

    The variance of the failure probability estimator is estimated as:

    The coefficient of variation(C.O.V)δTISof the failure probability estimator is then:

    2.4 Whale Optimization Algorithm

    The Whale Optimization Algorithm (WOA) is a novel nature-inspired population-based metaheuristic algorithm,presented by Seyedali Mirjalili and Andrew Lewis from Griffith University in 2016 [49].It was abstracted from the special hunting behavior of humpback whales.This foraging behavior is named bubble-net feeding method.

    Two approaches namely,shrinking encircling mechanism and spiral updating position,are designed to mathematically model the exploitation phase of the WOA.The first approach is formulated as follows:

    wheretindicates the current iteration,AiandCiare coefficients.x(t)iandrepresenti-thsearch agent and the best agent in thet-thiteration,respectively.

    whereqis linearly reduced from 2 to 0 over the course of iterations,ri1,ri2are random numbers in[0,1].Aiis a random value in the interval [?2q,2q].When the random value ofAiis in [?1,1],the next position of a search agent can be in any position of the agent and the position of the current best agent.The parameterCiis a random value in [0,2].This component provides random weights for prey in order to stochastically emphasize(Ci>1)or deemphasize(Ci<1)the effect of prey in defining the distance between the search agent and the current best agent.Besides,the second approach,the process of spiral updating position,is expressed as follows:

    wherebis a constant for defining the‘9’-shaped path of distinctive bubbles created by the humpback whales,lis a random number in[?1,1].

    The humpback whales hunt the prey by using the shrinking circle mechanism and spiral updating method simultaneously,and the two approaches have the same possibility to be chosen,so the model of the exploitation phase is as follows:

    whereri3are random numbers in[0,1]to switch the different mathematical model of exploitation.

    However,the process of search for prey is defined as the exploration phase,and the model of this phase is described as follows:

    wherex(t)jis thej-thrandom search agent which is different from thex(t)i,andAiis also defined by Eq.(19),but it is absolute value,|Ai|,is greater than 1 to emphasize exploration and let the WOA to search the whole search space.

    3 AWK-TIS Algorithm

    In this section,the main steps of the proposed method are introduced.Firstly,the WOA-Kriging model is constructed based on the basic Kriging model.Secondly,the AWK-TIS method is developed.

    3.1 WOA-Kriging Model

    In the proposed WOA-Kriging model,the optimal value of correlation parameterθ?is calculated by integrating WOA into the basic Kriging model.The pattern search method is adopted to seekθ?in DACE.However,the pattern search method is a local optimization algorithm that needs an initial solution to search the optimal and it is easy to get trapped in the local optimum.Alternatively,while seekingθ?,as an efficient global optimization algorithm not depend on the initial solution.WOA is embedded in the DACE toolbox to solve Eq.(8),instead of pattern search algorithm.Then,in our implementations the fitness function for the MLE according to the correlation parameterθis to minimize the function below[45]:

    Given the fitness function above,the finding process ofθ?becomes a typical minimization problem.Thus,the population(whale positions)is firstly initialized in an-dimensional search space,and the population size is set to 25 in this study.Each position corresponds to a correlation parameter vector.Next,a set of mathematical rules in WOA such as encircling prey,spiral bubble-net feeding and search for prey are repeatedly conducted at each iteration.Then,the fitness values of every population are obtained,the fittest solution is ascertained.Finally,the global optimal correlation parameter is acquired until the satisfaction of the stopping criteria that is the maximum number of iterations equals to 100.The main WOA parameters are unchanged during the optimization procedure.Consequently,the WOA-Kriging model is constructed by the optimal correlation parameterθ?.

    3.2 The AWK-TIS Algorithm

    In this subsection,the AWK-TIS algorithm is introduced in detail.In AWK-TIS,the real limit state function is substituted by a WOA-Kriging model; the failure probability and its coefficient of variation are computed by TIS method in an active learning way.AWK-TIS mainly consists of three components.Firstly,a WOA-Kriging model is constructed according to the optimal parameterθ?searched by WOA and the DOE generated through Latin Hypercube Sampling method in the entire design space.Secondly,MPP is evaluated by WOA and the evaluation process is treated as solving a constrained optimization.Thirdly,the active learning method and TIS are adopted together to conduct the selection of the best next sample and the evaluation of the probability.At this stage,the WOAKriging model is updated according to the enriched DOE.The general sketch of the AWK-TIS is shown in Fig.2,and the detailed procedure is proposed as follows:

    Step 1:Generation of the initial DOE.Latin Hypercube Sampling(LHS)method is employed to generate the initial samples in the whole design space.Then,the structural response or the real limit state function of every initial sample is evaluated.Thus,the initial DOE is composed of the initial samples and their corresponding responses.The number of the initial samples has an effect on the accuracy of the WOA-Kriging model;however,there are no certain rules to certificate this number.

    Step 2:Construction of the WOA-Kriging model.Construct the WOA-Kriging based on initial DOE by embedding WOA to DACE toolbox to find the optimal parameterθ?.Calculatingθ?will spend a little more computing time,but it is worth since the time consumed for finding accurate parameters could be ignored compared with that for time-consuming simulations.The initial value ofθis unnecessary in this metamodel construction,and the correlation function is chosen to be the Gaussian function.

    Step 3:Evaluation of the MPP.The WOA is utilized to calculate the MPP.The notion of reliability indexβ,is the minimum distance from the origin to the limit state surface in the standard normal space.The reliability index can be therefore calculated as the following constrained optimization problem[9,47,48]:

    wherexdenotes the vector of random variants,which isx=(x1,x2,...,xn);denotes the standard normal variants;μxiandσxirepresent the mean and standard deviation of the random variantxi,respectively.(x)is the limit state function constructed by WOA-Kriging model.

    Figure 2:General sketch of the AWK-TIS algorithm

    The problem in Eq.(25)is a constrained nonlinear optimization problem.The aim of constraint optimization is to search for feasible solutions with better objective values.However,constrained optimization problems are difficult to solve than unconstrained optimization problems due to the presence of constraints and their interrelationship between the objective functions.As a result,the main assignment while solving the constrained optimization problem is to deal with the constraints.

    Step 3.1:Construction of unconstrained function

    To avoid the violation of constraints,unfeasible solutions should be adapted to feasible solutions.In the meta-heuristic community,the most common technique is to take advantage of penalty functions to handle these constraints.Thus,a penalty function is employed to convert the constrained optimization problem in Eq.(25)to the unconstrained one in Eq.(26)[21]:

    where,λis the positive penalty coefficient,and its initial value is set to 103in this paper.

    Step 3.2:Calculation of reliability index by WOA

    The WOA is adopted to solve the unconstrained function.The inner loop is to find the optimal solution of the Eq.(26) by WOA,while the outer loop is to make sure that the reliability indexβis accurate enough.Thus,two stop criteria should be satisfied simultaneously as follows.First,the relative error between the reliability indices in two subsequent outer loop iterations is acceptable[47]:

    wherekandk+1 indicate thek-th andk+1-th iteration in the outer loop,respectively,and relative errorζis a small constant(e.g.,0.01).

    Second,it is required that the value of the WOA-Kriging model at the“potential”MPP should be relatively small which means that the“potential”MPP is rather close to the LSF[48]:

    whereχis the threshold which is also a small constant(e.g.,0.0001).

    If the two criteria above are met simultaneously,the MPP is acquired;otherwise,letλ=10λand go back to Step 3.1.It is noted that in this step the parameters in WOA are remain the same as its original version.

    Step 4:Generation of NIS samples.NISsamples are generated according to importance sampling density function centered on the MPP found by Step 3.The response of those samples will be evaluated by WOA-Kriging if the active learning procedure requires it.

    Step 5:Conduction of the active learning procedure.Add{xu,G(xu)}into the DOE set to rebuild the WOA-Kriging model,and evaluate the learning functionUamongNISsamples.The learning process continues until the stop criterion,i.e.,min(U(xi))>2 is satisfied,which means that the probability of making wrong sign estimation should not exceed 2.3%.

    Step 6:Evaluation of the failure probability and the coefficient of variation δTIS.After the WOAKriging metamodel is constructed by the active learning method,the next step is to adopt TIS to analyze the reliability.The failure probability and the coefficient of variation are estimatedδTISby Eqs.(15) and (17),respectively,and make sure that theδTISsatisfies the predefined threshold.IfδTISdoes not satisfy the threshold,a set of newNISsamples are generated to enrich the originalNISsamples.In this paper,the threshold of theδTISis taken as 0.05.

    Step 7:Output the unbiased estimation of the failure probability.IfδTISsatisfies the threshold,i.e.,δTIS≤0.05 then AWK-TIS terminates,and the failure probability is acquired.

    4 Validation Cases

    In this section,the efficiency of AWK-TIS is illustrated by four numerical cases and one engineering case.Each of the cases in this study has only a single failure region,which means only one MPP in each case.All the case results are averaged over 10 different runs.The accuracy and efficiency of different algorithms is compared according to the number of calls to the LSF (Ncall),δPf(the coefficient of variation of failure probability)andεPf(the relative percentage error of failure probability in comparison with the reference value).

    4.1 Case 1:2D Application

    A 2D numerical case,which has two independent random variables with standard normal distribution,is chosen to illustrate the process of the proposed method.The limit state function reads

    In the proposed method,an initial Kriging model is constructed based on 12 samples within the domain of[±6]inuspace.The results of different methods are summarized in Table 1.The proposed method AWK-TIS is compared with AK-IS,MetaAK-IS2,AK-SS.The reference value evaluated by MCS with 5×107samples is taken from[37],and the value is 2.85×10?5.In order to test the robustness of AWK-TIS,10 different runs are performed withNIS=1×105initial samples centered on the MPP confirmed by WOA.According toNcall,i.e.,the number of calls to the LSF,the proposed AWK-TIS outperforms the majority of active learning Kriging-based methods except for AK-SS;however,theδPf,i.e.,the coefficient of variation of failure probability of AK-SS is higher than that of AWK-TIS.

    Table 1: Results of the 2D application:Case 1

    Furtherly,the approximation process of the limit state function is schematically illustrated in Fig.3,including LSF(green line),Kriging model(red dotted line),initial samples(black fork dots),added samples(pink circle)and MPP(blue hexagon).Moreover,the MPP searched by WOA is located almost in the same position after 5 iterations while the enriched DOE are updated sequentially byUlearning function.The convergence process of failure probabilityPfand the C.O.V ofPffor AWK-TIS are presented in Fig.4.

    It can be seen from Figs.3 and 4 that these newly computed samples are located in the vicinity of the limit state function and are used to improve the accuracy of Kriging metamodel.Moreover,AWK-TIS requires much less performance function computations than AK-IS.Indeed,AWK-TIS needs 9 computations of real limit state function after the construction of Kriging model by initial DOE.

    As described in Section 2.4,the WOA is adopted to acquire the correlation parameterθ.The number of search agents is 25,the number of iterations is set to 100,and other parameters in WOA are remaining as the original form.Fig.5 depicts the convergence of the Kriging correlation parameterθshown in Eq.(8)when the update process of the WOA-Kriging model is stopped.Moreover,Fig.5 shows that the maximum likelihood estimation ofθconvergence to the minimum value of 6.49×10?6at the 16thiteration,and the optimumθis found at this point.Compared with PSO and GA,the WOA shows its high convergence speed and calculation accuracy.

    Figure 3:The approximation process of the LSF for case 1

    Figure 4:Convergence process of Pf and C.O.V of Pf for case 1

    Figure 5:The process to find θ?using WOA for case 1

    4.2 Case 2:A Cantilever Beam Example

    This case has been widely used in literature [38].The cantilever beam,which has a rectangular cross-section,is subjected to a distributed uniform load.The limit state function which is defined by the maximum vertical displacement of the beam is given as

    wherex1andx2are assumed to be independent random variables following normal distribution.The statistical properties of the above two variables are shown in Table 2.The results obtained by the proposed AWK-TIS,MSC,AK-MCS,and AK-SS are listed in Table 3.

    Table 2: Statistical properties of the variables for case 2

    Table 3: Results of case 2

    From Table 3,it can be seen that theNcallof the proposed AWK-TIS is the same as that of AK-SS,but it is less than the number of samples applied in AK-MCS.Besides,the one can see that all thePfcalculated by the different methods have high accuracy comparing with the MCS solution.However,the coefficient of variation of failure probability of AWK-TIS is the smallest in the four methods.

    Similarly to Figs.3 and 4,Fig.6 shows the final approximation of AWK-TIS and Fig.7 shows the convergence process of failure probability and the C.O.V of failure probability.As shown in Fig.6,the added samples are almost in the vicinity of the real LSF,and the LSF constructed by Kriging can be accurately approximated in the domain of the sampling region.On the contrary,the Kriging model shows inexact outside the sampling region.Nevertheless,it should be indicated that the region outside the sampling area has little influence on the failure probability calculation due to the MPP is far away from this area.From Fig.7,the proposed method has converged to the final results at about 12 added samples,which indicates the efficiency and accuracy of the AWK-TIS.

    Figure 6:The approximation of the LSF for case 2

    Figure 7:Convergence process of Pf and C.O.V of Pf for case 2

    4.3 Case 3:A Non-Linear Oscillator

    A non-linear undamped single degree of freedom system depicted in Fig.8 is analyzed in this case.The limit state function is given as[54]

    whereω0=Six random input variables includingC1,C2,M,R,T1,andF1are listed in Table 4.The initial number of DOE is set as 20.The AWK-TIS is compared with several other existing methods (MCS,AK-MCS,PAK-Bn[54]) and their results are listed in Table 5.Same as the two cases above,the reference values to compare the reliability analysis results are thePf,,C.O.V ofPf,δPf,εPf,estimated by MCS.

    Figure 8:A non-linear oscillator

    Table 4: Statistical properties of the variables for case 3

    Table 4 (continued)VariableMeanStandard deviationDistribution T110.2Normal F10.60.2Normal

    Table 5: Results of case 3

    It can be seen from Table 5 that although the relative error of failure probability obtained by the proposed method is slightly larger than that of AK-MCS and PAK-Bn,theNcallof AWK-TIS is significantly reduced contrasted with that of the two methods above.Additionally,the C.O.V of failure probability is rather small compared with the other three methods.Hence,it indicates that AWK-TIS could solve the problem of small failure probability with a moderate number of random variables.

    Equally,Fig.9 presents the convergence process of failure probability and C.O.V of failure probability of this case.From Fig.9,one can see that the failure probability has converged with almost 50 added samples,which means that the proposed AWK-TIS is suitable for the calculation of the small failure probability.

    Figure 9:Convergence process of Pf and C.O.V of Pf for case 3

    4.4 Case 4:A Roof Truss Structure

    A roof truss structure,which is selected to verify the proposed method,is shown in Fig.10.The roof truss undertaking the uniformly distributed loadqthat can be converted to the nodal loadP=ql/4.The response of the roof truss is the perpendicular deflection △Cof nodeCis denoted as[55]:

    whereAC,AS,EC,ESrespectively are sectional area,elastic modulus concrete and steel bars,lis the length of the structure.The mean and variation coefficient of the six independent normal distribution random variables are reported in Table 6.The threshold of △Cshould not be larger than 0.03 m,and the LSF of this case is thus defined as

    Figure 10:The schematic diagram of a roof truss

    Table 6: Statistical properties of the variables for case 3

    The AWK-TIS is compared with other existing methods (MCS,AK-MCS,KAIS,RBF) and the results are summarized in Table 7.The results by MCS calculations with 1.06×106samples are regarded as reference values in this case.

    Table 7: Results of case 4

    Table 7 (continued)MethodNcallPfδPfεPf(%)RBF[25]12+1439.736×10?3/3.629 AWK-TIS20+799.357×10?30.010.404

    It can be seen from Table 7 that AWK-TIS needs only 99Ncallto acquire a satisfactory result while the other three approaches demand more than 150Ncall.The smallNcallhas highlighted the numerical efficiency of the proposed approach in estimating the roof truss structure.Furthermore,the AWK-TIS has the second relative error of 0.404%comparing with the reference value ofPfevaluated by MCS,followed by AK-MCS and RBF.

    By investigating the convergence process ofPfand C.O.V ofPfpresented in Fig.11,we can see that both of the curves converged at about 25 added samples,which means that the proposed AWK-TIS is suitable for dealing with small failure probabilities and relatively large number of random variables.Additionally,the convergence threshold can be advisably loosed to further improve computational efficiency.

    Figure 11:Convergence process of Pf and C.O.V of Pf for case 4

    4.5 Case 5:A Jet Vane Manipulation Mechanism

    Vertical launch missile is widely used on land air defense and sea air defense.The performance of a jet vane system is extremely vital to the flying control of vertical launch missile.As a significant component of the jet van system,the manipulation mechanism drives the rudder rotating in the engine gas flow that generates pitch moment to control the orientation of the missile,which is shown in Fig.12a.The manipulate mechanism convers the linear motion of the servo motor to the rotation of the jet van,and the schematic of this mechanism is shown in Fig.12b.

    Figure 12:A jet van manipulation mechanism

    In Fig.12b,AandBrepresent the positions of the two ends of the linkage respectively,andCis the center of the shaft.L1andL2are the length of linkage and rocker respectively.L3is the distance betweenAandCalong theX-axial.ais the distance betweenAandC.Sis the maximum displacement ofAalong the negative direction ofX-axial,andAl,Blare the left limit position of the linkage.θ2is the angle between theBCand theX-axial,and it is the function ofL1,L2,L3,aandS,which is shown in Eq.(34)

    whereθ1is the angle between the linkage and the X-axial.

    In this case,the linkage lengthL1,the length of the rockerL2,the distance betweenAandCalong the X-axialL3,and the maximum displacement of linear servo systemSare taken as random variables.Table 8 lists the details of random variables.To assure the normal operation of the manipulation mechanism,the maximum angle of the rudder shaft should not exceed a given threshold.Therefore,the limit state function is defined as

    whereθmis the threshold of the maximum angle of the rudder shaft which is taken as 118.2°,θ2(X)is the maximum angle of the rudder shaft which is calculated by solving Eq.(34),theXrepresents all the random variables.

    Table 8: Statistical properties of the variables for case 5

    The results of the different methods are shown in Table 9,in which results estimated by the AKMCS and AK-IS methods are also provided for the sake of comparison.Besides,the calculation times or CPU-time of different AK-based algorithm are also exhibited in this table.The sample pool of the MCS,AK-MCS and AK-IS are 5×107,5×106and 4×106,respectively.

    Table 9: Results of the case 5

    Table 9 shows that the predictive accuracy of the AWK-TIS notably outperforms the other two AK methods considering the relative error of failure probability.Besides,theNcallis less than that of AK-MCS and AK-IS,which indicates that the proposed AWK-TIS can solve real engineering problem efficiently.From the CPU-time results,it is concluded that although the process of optimizing Kriging model by WOA algorithm additionally increases the computational time,the utilization of the TIS method reduces the calculation time.

    The convergence of failure probability and coefficient of variation of failure probability is presented in Fig.13.From Fig.13,thePfand C.O.V ofPfhave converged to the final results at about the 10 added samples.Thus,the computational performance of the AWK-TIS when dealing with this case can be improved with a loose stopping criterion.

    Figure 13:Convergence process of Pf and C.O.V of Pf for case 5

    5 Conclusions

    The AWK-TIS method is proposed as a novel reliability method to deal with structural reliability problems,in which the WOA is performed to seek the optimum correlation parameter of Kriging model and the MPP,and then the added samples generated byUlearning function are enriched to DOE to refine the Kriging model.The MPP is calculated by WOA through solving the constrained nonlinear optimization problem until certain two thresholds are satisfied simultaneously.Afterward,the TIS method is applied to calculate the failure probability and coefficient of variation of the failure probability.In AWK-TIS,the sample pool is much smaller than that of AK-IS,which significantly expedites the learning efficiency.Five test cases including four numerical cases and one engineering case are adopted to verify the performance of the AWK-TIS method.Results prove that the proposed method can achieve high computational accuracy and efficiency.

    AWK-TIS method has two main drawbacks:(1)AWK-TIS,like AK-IS,is dependent on the MPP to ensure the next best point; however,the solitary sampling center has confined its employment in dealing with multi-MPP problems,(2)the process of acquiring the MPP is transformed to solve the constrained nonlinear optimization problem,which may affect the efficiency of the proposed method,(3) due to the restriction of IS-based algorithm,the proposed algorithm does have some defects in dealing with high-dimensional problems.Future works will focus on three aspects: (1) adopt other efficient sampling methods in lieu of TIS to solve multi-MPP problems,(2)improve the efficiency of WOA by hybridising other optimizers or employing other more effective optimization algorithms to solve the constrained optimization problem to obtain the MPP,(3)alternative modeling and analysis methods need to be incorporated for effective validation of AWK-TIS to high-dimensional problems.

    Funding Statement: This work is supported by the Technical Basic Scientific Research Projects of State Administration of Science,Technology and Industry for National Defence,PRC (Grant No.JSZL2019204C001).

    Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.

    亚洲片人在线观看| а√天堂www在线а√下载| 激情在线观看视频在线高清| 又爽又黄无遮挡网站| 亚洲精品影视一区二区三区av| 亚洲美女视频黄频| 国产亚洲精品久久久久久毛片| 淫妇啪啪啪对白视频| 中亚洲国语对白在线视频| 久久精品国产清高在天天线| 精品一区二区三区人妻视频| 制服人妻中文乱码| 一边摸一边抽搐一进一小说| 欧美午夜高清在线| 日韩中文字幕欧美一区二区| 午夜日韩欧美国产| 国产成人a区在线观看| 婷婷丁香在线五月| 国产精品久久久久久人妻精品电影| 午夜福利在线在线| 在线国产一区二区在线| 成人18禁在线播放| 欧美激情在线99| 综合色av麻豆| 亚洲无线在线观看| 深爱激情五月婷婷| 日日干狠狠操夜夜爽| 18禁在线播放成人免费| 国产伦人伦偷精品视频| 性色av乱码一区二区三区2| 欧美黄色片欧美黄色片| 亚洲精品粉嫩美女一区| 亚洲av日韩精品久久久久久密| 精品久久久久久久毛片微露脸| 九九久久精品国产亚洲av麻豆| 国产免费av片在线观看野外av| 18禁在线播放成人免费| 国产精华一区二区三区| 狂野欧美激情性xxxx| 制服人妻中文乱码| 国产黄片美女视频| 香蕉久久夜色| 国产激情欧美一区二区| 欧美丝袜亚洲另类 | 欧美激情在线99| 男女视频在线观看网站免费| 成人无遮挡网站| 一个人观看的视频www高清免费观看| 免费看十八禁软件| 欧美色视频一区免费| 一级作爱视频免费观看| 一级黄片播放器| 蜜桃亚洲精品一区二区三区| 国产色爽女视频免费观看| 亚洲av中文字字幕乱码综合| 精品国产亚洲在线| 村上凉子中文字幕在线| 国产三级黄色录像| 久久久久九九精品影院| 欧美一级毛片孕妇| 免费一级毛片在线播放高清视频| 亚洲中文字幕一区二区三区有码在线看| 国产一区二区在线观看日韩 | 欧美极品一区二区三区四区| 国产熟女xx| 日韩中文字幕欧美一区二区| 欧美bdsm另类| 18禁在线播放成人免费| 日本撒尿小便嘘嘘汇集6| 真实男女啪啪啪动态图| 女同久久另类99精品国产91| 久久99热这里只有精品18| 精品久久久久久久毛片微露脸| 在线观看免费午夜福利视频| 夜夜躁狠狠躁天天躁| 少妇熟女aⅴ在线视频| 我的老师免费观看完整版| 男人舔女人下体高潮全视频| 亚洲av熟女| 亚洲人与动物交配视频| 欧洲精品卡2卡3卡4卡5卡区| 国产午夜精品论理片| 欧美日韩黄片免| 国产蜜桃级精品一区二区三区| 偷拍熟女少妇极品色| 免费大片18禁| 99在线视频只有这里精品首页| 午夜精品在线福利| 免费看美女性在线毛片视频| 久9热在线精品视频| 国内久久婷婷六月综合欲色啪| 伊人久久精品亚洲午夜| 在线视频色国产色| 一个人免费在线观看的高清视频| 美女 人体艺术 gogo| 最新美女视频免费是黄的| 欧美绝顶高潮抽搐喷水| 成人av一区二区三区在线看| 久久精品影院6| 久久6这里有精品| 少妇人妻精品综合一区二区 | 成人精品一区二区免费| 日韩免费av在线播放| 国产91精品成人一区二区三区| 亚洲av第一区精品v没综合| 在线观看免费午夜福利视频| av专区在线播放| 夜夜看夜夜爽夜夜摸| 丁香欧美五月| 亚洲欧美一区二区三区黑人| 欧美中文综合在线视频| 女同久久另类99精品国产91| 亚洲,欧美精品.| 一个人免费在线观看电影| 黄色片一级片一级黄色片| 夜夜躁狠狠躁天天躁| 国产老妇女一区| 久久亚洲精品不卡| av在线天堂中文字幕| 成人av在线播放网站| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 亚洲无线观看免费| 国产男靠女视频免费网站| 久久香蕉国产精品| 51国产日韩欧美| 国产精品一区二区免费欧美| 久久久久久人人人人人| 欧美黑人欧美精品刺激| 日本一二三区视频观看| 免费人成在线观看视频色| 久久精品人妻少妇| 亚洲av熟女| 国产v大片淫在线免费观看| 一本综合久久免费| 久久亚洲精品不卡| 操出白浆在线播放| 天天一区二区日本电影三级| 国产精品一区二区三区四区久久| 很黄的视频免费| 国产探花在线观看一区二区| 91字幕亚洲| 90打野战视频偷拍视频| 99久国产av精品| a级一级毛片免费在线观看| 精品国产亚洲在线| 岛国视频午夜一区免费看| 欧美激情在线99| 亚洲欧美精品综合久久99| 一本综合久久免费| 国产伦精品一区二区三区视频9 | 亚洲成av人片在线播放无| 久久久久亚洲av毛片大全| 国产亚洲欧美在线一区二区| 午夜福利成人在线免费观看| 一区二区三区激情视频| 国产精品一区二区三区四区久久| 免费在线观看亚洲国产| 欧美精品啪啪一区二区三区| 久久久成人免费电影| 夜夜爽天天搞| 国产激情偷乱视频一区二区| 亚洲精品在线观看二区| 午夜福利在线观看吧| 最新美女视频免费是黄的| 12—13女人毛片做爰片一| 国产一区二区亚洲精品在线观看| 亚洲第一欧美日韩一区二区三区| 国产午夜福利久久久久久| 五月玫瑰六月丁香| 久久久国产成人精品二区| 99久久九九国产精品国产免费| 日韩欧美一区二区三区在线观看| aaaaa片日本免费| 国产成+人综合+亚洲专区| 在线观看一区二区三区| 国产成人影院久久av| 高潮久久久久久久久久久不卡| 亚洲 国产 在线| 久久久久久大精品| 欧美极品一区二区三区四区| 亚洲 欧美 日韩 在线 免费| 亚洲无线在线观看| 亚洲avbb在线观看| 小说图片视频综合网站| 天堂√8在线中文| 日韩有码中文字幕| 好看av亚洲va欧美ⅴa在| 欧美bdsm另类| 国产美女午夜福利| 国产熟女xx| 国产不卡一卡二| 中国美女看黄片| 久久国产精品人妻蜜桃| 久久人人精品亚洲av| 性色av乱码一区二区三区2| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 99国产极品粉嫩在线观看| 18禁美女被吸乳视频| 精品久久久久久成人av| 日本成人三级电影网站| 琪琪午夜伦伦电影理论片6080| 国产v大片淫在线免费观看| 免费看光身美女| 国产高清激情床上av| a在线观看视频网站| 丁香六月欧美| 国产免费av片在线观看野外av| 黄片大片在线免费观看| 无限看片的www在线观看| 男女下面进入的视频免费午夜| 日本一二三区视频观看| 欧美精品啪啪一区二区三区| 国产精品一区二区三区四区久久| 首页视频小说图片口味搜索| 国产老妇女一区| 九色成人免费人妻av| 狂野欧美激情性xxxx| 国产亚洲精品久久久久久毛片| 欧美性猛交黑人性爽| 成人亚洲精品av一区二区| 精品一区二区三区视频在线观看免费| 国产精品爽爽va在线观看网站| 老司机福利观看| 九九在线视频观看精品| 90打野战视频偷拍视频| 九色国产91popny在线| 国产免费一级a男人的天堂| 深爱激情五月婷婷| 啦啦啦观看免费观看视频高清| 少妇人妻一区二区三区视频| 好男人在线观看高清免费视频| 国产午夜精品论理片| 搡女人真爽免费视频火全软件 | 亚洲va日本ⅴa欧美va伊人久久| 久久久久久久久大av| 操出白浆在线播放| 国产精品98久久久久久宅男小说| 亚洲无线在线观看| 级片在线观看| 亚洲精品亚洲一区二区| 中亚洲国语对白在线视频| 一级毛片女人18水好多| 久久这里只有精品中国| www国产在线视频色| 夜夜看夜夜爽夜夜摸| 午夜亚洲福利在线播放| 国产探花在线观看一区二区| 久久精品国产自在天天线| 国产精品嫩草影院av在线观看 | e午夜精品久久久久久久| 成人国产综合亚洲| 午夜福利免费观看在线| 日本一本二区三区精品| 国产一区二区在线观看日韩 | 又紧又爽又黄一区二区| 久久久精品欧美日韩精品| 亚洲欧美精品综合久久99| 999久久久精品免费观看国产| 偷拍熟女少妇极品色| 国产又黄又爽又无遮挡在线| 日韩精品青青久久久久久| 亚洲av一区综合| 午夜久久久久精精品| 亚洲真实伦在线观看| 免费搜索国产男女视频| 亚洲无线观看免费| 亚洲18禁久久av| 日本黄色片子视频| 12—13女人毛片做爰片一| 日日干狠狠操夜夜爽| 精品一区二区三区av网在线观看| 亚洲精品日韩av片在线观看 | 两个人看的免费小视频| 午夜免费成人在线视频| 日本一二三区视频观看| 国产精品野战在线观看| 久久99热这里只有精品18| 亚洲精品久久国产高清桃花| 午夜免费激情av| 午夜久久久久精精品| 午夜免费成人在线视频| 99精品久久久久人妻精品| 国产毛片a区久久久久| 五月玫瑰六月丁香| 久久久久精品国产欧美久久久| 亚洲精华国产精华精| 特大巨黑吊av在线直播| 久久久久久人人人人人| 亚洲av免费高清在线观看| 老熟妇乱子伦视频在线观看| 99国产精品一区二区蜜桃av| 免费看光身美女| 久久久精品欧美日韩精品| 国产成人av激情在线播放| 国产精品av视频在线免费观看| 免费在线观看成人毛片| 男人舔奶头视频| 高清毛片免费观看视频网站| 欧美日韩精品网址| 免费搜索国产男女视频| 内射极品少妇av片p| 19禁男女啪啪无遮挡网站| 黄色丝袜av网址大全| 97人妻精品一区二区三区麻豆| 精品熟女少妇八av免费久了| 身体一侧抽搐| 久久性视频一级片| 国产亚洲精品久久久久久毛片| av视频在线观看入口| 9191精品国产免费久久| 国产黄a三级三级三级人| 国产精品 欧美亚洲| 国产真实伦视频高清在线观看 | 亚洲精品影视一区二区三区av| 久久精品亚洲精品国产色婷小说| 亚洲狠狠婷婷综合久久图片| 成人欧美大片| 久久中文看片网| 一卡2卡三卡四卡精品乱码亚洲| 91字幕亚洲| 日本一二三区视频观看| 一级黄色大片毛片| 亚洲五月天丁香| 日本免费一区二区三区高清不卡| 青草久久国产| 黄色日韩在线| 久久亚洲真实| 亚洲最大成人中文| 亚洲七黄色美女视频| 午夜免费男女啪啪视频观看 | 国产探花极品一区二区| 午夜福利免费观看在线| 中文字幕av在线有码专区| 亚洲国产欧洲综合997久久,| 男女那种视频在线观看| 亚洲av免费在线观看| 成人亚洲精品av一区二区| 日本与韩国留学比较| 别揉我奶头~嗯~啊~动态视频| 亚洲人成网站在线播放欧美日韩| 午夜老司机福利剧场| 成人永久免费在线观看视频| 天天添夜夜摸| 91九色精品人成在线观看| 久久久久久久午夜电影| 麻豆一二三区av精品| 麻豆国产97在线/欧美| 国产久久久一区二区三区| 亚洲国产精品sss在线观看| 国产日本99.免费观看| 一区二区三区高清视频在线| 十八禁人妻一区二区| www日本在线高清视频| 欧美一区二区亚洲| 国产成+人综合+亚洲专区| 中文字幕熟女人妻在线| 在线国产一区二区在线| 18禁在线播放成人免费| 成人永久免费在线观看视频| 一个人免费在线观看电影| 亚洲国产日韩欧美精品在线观看 | 国产精品1区2区在线观看.| 一夜夜www| 国产伦在线观看视频一区| 女生性感内裤真人,穿戴方法视频| 国产激情欧美一区二区| 成年女人永久免费观看视频| netflix在线观看网站| 国产精品三级大全| 啦啦啦观看免费观看视频高清| 国产精品,欧美在线| 国产在视频线在精品| 草草在线视频免费看| 亚洲男人的天堂狠狠| 啦啦啦免费观看视频1| 久久草成人影院| 国产午夜精品久久久久久一区二区三区 | 亚洲最大成人手机在线| 国产一级毛片七仙女欲春2| 国产精品久久视频播放| 熟女电影av网| 成年人黄色毛片网站| 成人特级黄色片久久久久久久| 国产精品 国内视频| 又黄又爽又免费观看的视频| 老鸭窝网址在线观看| 免费无遮挡裸体视频| 久久久久久久精品吃奶| 白带黄色成豆腐渣| 嫁个100分男人电影在线观看| 99视频精品全部免费 在线| 99riav亚洲国产免费| 国产黄片美女视频| 99riav亚洲国产免费| 国产亚洲精品久久久久久毛片| 免费看十八禁软件| 日韩中文字幕欧美一区二区| 欧美日韩一级在线毛片| 亚洲精品亚洲一区二区| 日本黄色片子视频| 精品久久久久久,| 婷婷丁香在线五月| 久久精品国产亚洲av香蕉五月| 精品99又大又爽又粗少妇毛片 | 亚洲自拍偷在线| 欧美黄色淫秽网站| eeuss影院久久| 亚洲精品乱码久久久v下载方式 | 免费av观看视频| 综合色av麻豆| 国产精品99久久久久久久久| 丝袜美腿在线中文| 女同久久另类99精品国产91| 久久午夜亚洲精品久久| 日韩人妻高清精品专区| 99久久精品国产亚洲精品| 无遮挡黄片免费观看| 狂野欧美激情性xxxx| 亚洲精品成人久久久久久| 蜜桃亚洲精品一区二区三区| 免费av毛片视频| 欧美大码av| 久久中文看片网| 淫妇啪啪啪对白视频| 99在线人妻在线中文字幕| 99热精品在线国产| 偷拍熟女少妇极品色| 免费看光身美女| 欧美在线一区亚洲| 亚洲成人免费电影在线观看| 最近在线观看免费完整版| 精品电影一区二区在线| 久久久国产成人精品二区| 3wmmmm亚洲av在线观看| 高清在线国产一区| 亚洲熟妇中文字幕五十中出| 在线观看舔阴道视频| 日韩欧美免费精品| 激情在线观看视频在线高清| 久久九九热精品免费| 三级国产精品欧美在线观看| 午夜激情欧美在线| 丁香欧美五月| 国产 一区 欧美 日韩| xxxwww97欧美| 无遮挡黄片免费观看| 成熟少妇高潮喷水视频| 国产精品久久久久久精品电影| 男女那种视频在线观看| 色综合站精品国产| 国产爱豆传媒在线观看| 99热6这里只有精品| 精品久久久久久,| 欧美日韩亚洲国产一区二区在线观看| 搡女人真爽免费视频火全软件 | 精品久久久久久,| 成熟少妇高潮喷水视频| 五月伊人婷婷丁香| 国内久久婷婷六月综合欲色啪| 国产精品免费一区二区三区在线| 国产欧美日韩精品一区二区| 欧美性感艳星| 一个人看的www免费观看视频| 亚洲最大成人手机在线| 神马国产精品三级电影在线观看| 国产97色在线日韩免费| 男人舔女人下体高潮全视频| 真实男女啪啪啪动态图| 国产一区二区三区在线臀色熟女| 又黄又爽又免费观看的视频| 黄片大片在线免费观看| 久久久久亚洲av毛片大全| 久久草成人影院| 国产蜜桃级精品一区二区三区| 免费在线观看日本一区| 欧美性感艳星| 国产精品亚洲一级av第二区| 女警被强在线播放| 国产成人系列免费观看| av黄色大香蕉| 有码 亚洲区| 国产探花极品一区二区| 99riav亚洲国产免费| 国产精品综合久久久久久久免费| 美女高潮的动态| 搡女人真爽免费视频火全软件 | 国产精品av视频在线免费观看| 国产高清三级在线| 亚洲av电影在线进入| 亚洲欧美日韩东京热| 免费搜索国产男女视频| 我的老师免费观看完整版| xxxwww97欧美| 亚洲激情在线av| 久久久久久久亚洲中文字幕 | 美女高潮的动态| 亚洲美女视频黄频| 看片在线看免费视频| 日韩欧美免费精品| 国产精品久久久久久久电影 | 国产精品香港三级国产av潘金莲| 久久国产精品人妻蜜桃| 国产一区二区三区在线臀色熟女| 国产免费男女视频| 成年人黄色毛片网站| av天堂在线播放| 欧美一级毛片孕妇| 亚洲av成人av| 久久精品国产亚洲av涩爱 | 国产色爽女视频免费观看| 久久久国产精品麻豆| 天天一区二区日本电影三级| 亚洲五月婷婷丁香| 久久欧美精品欧美久久欧美| 法律面前人人平等表现在哪些方面| 十八禁网站免费在线| 欧美日韩国产亚洲二区| 免费看十八禁软件| 波多野结衣高清无吗| 欧美极品一区二区三区四区| 亚洲熟妇中文字幕五十中出| 天天一区二区日本电影三级| 久久九九热精品免费| 国产av在哪里看| 日韩欧美国产在线观看| 一个人免费在线观看的高清视频| h日本视频在线播放| 精品熟女少妇八av免费久了| eeuss影院久久| 国内精品久久久久精免费| 天天添夜夜摸| 小说图片视频综合网站| 国产中年淑女户外野战色| 亚洲精品在线观看二区| 丰满人妻熟妇乱又伦精品不卡| 九九热线精品视视频播放| 婷婷亚洲欧美| 狂野欧美白嫩少妇大欣赏| 欧美bdsm另类| 窝窝影院91人妻| 91在线精品国自产拍蜜月 | 亚洲精品456在线播放app | 亚洲美女视频黄频| a级一级毛片免费在线观看| 丝袜美腿在线中文| 欧美一区二区精品小视频在线| 国产单亲对白刺激| 亚洲熟妇中文字幕五十中出| 亚洲va日本ⅴa欧美va伊人久久| 真人一进一出gif抽搐免费| 久久国产乱子伦精品免费另类| 99热这里只有是精品50| 欧美一区二区亚洲| 亚洲 欧美 日韩 在线 免费| 中出人妻视频一区二区| 一二三四社区在线视频社区8| 一个人看视频在线观看www免费 | 国产av在哪里看| 国产在视频线在精品| 舔av片在线| h日本视频在线播放| 18禁裸乳无遮挡免费网站照片| 国产午夜福利久久久久久| 两人在一起打扑克的视频| 高清日韩中文字幕在线| 亚洲在线观看片| 亚洲内射少妇av| 久久精品国产综合久久久| 中亚洲国语对白在线视频| 日韩精品青青久久久久久| 午夜久久久久精精品| 18禁在线播放成人免费| 欧美三级亚洲精品| 国产97色在线日韩免费| 精品乱码久久久久久99久播| 国产成+人综合+亚洲专区| 欧美色视频一区免费| 日韩 欧美 亚洲 中文字幕| 美女被艹到高潮喷水动态| 久久天躁狠狠躁夜夜2o2o| 成人国产一区最新在线观看| 国产av一区在线观看免费| 乱人视频在线观看| 国内精品美女久久久久久| 成人亚洲精品av一区二区| 九九热线精品视视频播放| 国产老妇女一区| 欧美丝袜亚洲另类 | 精品久久久久久,| 中文字幕久久专区| 亚洲精华国产精华精| 99国产极品粉嫩在线观看| 亚洲在线自拍视频| av视频在线观看入口| 日韩精品青青久久久久久| 中文字幕人妻熟人妻熟丝袜美 | 熟女电影av网| 成年女人毛片免费观看观看9| 女生性感内裤真人,穿戴方法视频| 国产精品98久久久久久宅男小说| 国产97色在线日韩免费| 一二三四社区在线视频社区8| 婷婷精品国产亚洲av| 97超级碰碰碰精品色视频在线观看| 亚洲国产欧洲综合997久久,| 国产精品国产高清国产av| 操出白浆在线播放| 亚洲,欧美精品.| 欧美日韩瑟瑟在线播放| 51国产日韩欧美| 久久中文看片网| 99久久九九国产精品国产免费| 99热6这里只有精品| 亚洲熟妇中文字幕五十中出| 精品福利观看|