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

    Online complex nonlinear industrial process operating optimality assessment using modi fied robust total kernel partial M-regression☆

    2018-05-25 11:26:24FeiChuWeiDaiJianShenXiaopingMaFuliWang

    Fei Chu *,Wei DaiJian Shen Xiaoping Ma Fuli Wang

    1 School of Information and Control Engineering,China University of Mining and Technology,Xuzhou 221116,China

    2 State Key Laboratory of Integrated Automation for Process Industries,Northeastern University,Shenyang 110819,China

    3 College of Information Science and Engineering,Northeastern University,Shenyang 110819,China

    1.Introduction

    In modern industrial processes,optimization of process operating performance has become a crucial issue and has received tremendous attention from both academia[1–4]and practice[5]recently because,in the face of growing competition,it represents the natural advantage for maximizing the economic profits.However,it is well known that the industrial process operating performance may deviate from the preliminary designed optimal operating point due to the disturbances and uncertainties.Therefore,process operating optimality assessment becomes important for industrial processes.

    Operating optimality assessment is an issue to answer the question“How optimal the current operating performance is?”[6],in other words it is to decide whether the operating performance is optimal or non-optimal and to get a measure on how far the current operating performance is from the optimum under the normal operating conditions[7].Namely,based on the assessment results,managers and operators can further adjust the industrial production and improve process operating performance.

    In recent years,some works have been published on the operating optimality assessment gradually[5–8].Ye et al.[5]proposed a probabilistic framework of operating safety and optimality assessment for multi-mode industrial processes,where an optimality assessment index is de fined and used to evaluate the process operating optimality.However,the assessment index is constructed based on the value of optimized objective function,such as plantcost,profits,productquality,the difficulty of calculating plant cost/profits online,and delay of production quality test,which makes the approach not suitable for the practical industrial application online,especially for the timeconsuming and large-scale complex industrial processes.Recognizing that the output-related process variation information is different between different performance grades but are maintained consistent in one performance grade,Liu et al.[7,8]proposed a performancesimilarity-based online assessment approach.They divided the training data into several steady performance grade based on the levels of an economic index,and then the process variation information of each steady performance grade was extracted by the featured extraction methods,such as principal component analysis(PCA)and total projection to latentstructures(T-PLS)[10,11],and used for operating optimality assessment.Although the PCA-based and T-PLS-based assessment methods are similar to some extent,the latter can overcome the low anti-interference ability and weak sensitivity drawbacks presented by the former,by removing the unrelated process variation with respect to the operating performance.Further in,Liu et al.proposed several improved versions of the performance-similarity and T-PLS-based online assessment approach to handle the nonlinearity issue[12],non-Gaussian multimode processes[13],and the misalignment between the actual process measurements and corresponding economic index[5].

    However,for real industrial process,the obtained process data may be subjected to outliers[14],which originates from sensor faults,missing data,process disturbances,heavy tails of noise distribution functions,malfunction of instruments,and so on;since PLS is known to be very sensible to outliers[15],the existence of a very small amount of outliers may lead to deceptive results of T-PLS-based operating optimality assessment methods.The approaches for operating optimality assessment listed above did not consider this issue,despite Liu et al.mentioning the outliers in reference[7].Based on the assumption that the data corresponding to a performance grade mainly distribute in an ellipse center area of the dataset,while the edge area of the ellipse is occupied by the outliers,the outliers were merely detected and removed from the modeling datasets by calculating the sample-similarity between each sample and the ellipse center of its data set before the assessment model for each performance grade were established.Nevertheless,for process with strong nonlinearities,the above assumption is usually not valid,and outliers are hard to be detected absolutely especially for processes with heavy-tail noise.Further,if tight sample-similarity threshold is used,too much “outliers”will be detected and removed,leading to insuf ficient samples especially for some performance grades which naturally lack samples.Meanwhile,outliers may contain some important process variations related to the operating performance,which should be extracted to establish exhaustive assessment models for each performance grade.

    To minimize the adverse effects of outliers,several robust PLS methods have been proposed[14].Among different robust PLS methods,a version of iteratively reweighted least squares algorithm for PLS has been obtained by Cumms and Andrew,called iteratively reweighted PLS.This method uses regression residuals to examine whether an object is an outlier or not,and calculates weight values to suppress the outliers in the model building step.Compared with other robust PLS methods,reweighted least squares algorithms possess low computational effort[15].Yet,this method is not resistant to bad leverage points(outliers in the space of input variables)[16].To solve this problem,partial robust M-regression(PRM)is proposed.With an appropriate choosing–weighting scheme,bad leverage points are also down-weighted at the same time.Further,a nonlinear version PRM[15],kernel Partial Robust M-regression,KPRM,is proposed which can cope with nonlinear robust problem.To the best of the author's knowledge,however,there have been few reports concerning robust nonlinear T-PLS method despite that the nonlinear robust problem is much more important in practical industrials.

    In this paper,we explore the issue of operating optimality assessment for complex industrial process based on performance-similarity considering nonlinearities and outliers simultaneously,and a general enforced online performance assessment framework is proposed,which consists of of fline training,online assessment and non-optimal cause identification.In practical,the process operating optimality performance is closely associated with optimal index(such as plant cost,profits or product quality),and in the of fline training,the training samples can be divided into different data sets according to the level of optimal index,and each data set roughly represents one performance grade;then,total projection to latent structures algorithm will be used to extract the process variations with respect to the optimal index to establish the performance assessment models for each performance grade.Considering the nonlinearities and robustness simultaneously,the primary version of T-PLS is not a good choice.A more powerful T-PLS algorithm is needed,which can cope with nonlinearities and minimize the adverse effects of outliers concurrently,establishing exhaustive assessment models for each performance grade.In the present study,we propose a novel modi fied robust nonlinear T-PLS method,called Total Kernel Partial M-regression(T-KPRM),by extending the choosing–weighting scheme to the kernelbased nonlinear T-PLS algorithm.T-KPRMalgorithm can also be treated as a post-processing method to decompose the KPRM kernel principal space and residual space,in which,with an appropriate choosing–weighting scheme,both bad leverage points and high residual points have been down-weighted.

    In the online assessment,based on the understanding that the optimalindex related process variation information differbetween different performance grades,but maintain consistent among the same performance grades,which is also known as performance-similarity[8],the current process operating performance can be evaluated by calculating the performance-similarities between online data and each performance grade modeling data.For the conversion processes between performance grades,the optimal index related process variation would change gradually from one performance grade to another,thus performance-similarities with an ascending rule can also be used to distinct the performance grade conversion.However,the existence of process disturbances and outliers online often make the assessment results not credible.Thus,an enforced online assessment procedure is proposed,and for the non-optimal operating performance,the responsible variables can be identified by the variable contribution-based nonoptimal cause identification method.

    The contribution of the present paper are summarized as follows:(1)a robust and nonlinear version of total projection to latent structure is proposed,and based on it,a novel and robust operating optimal assessment approach is developed which considers nonlinearities and outliers simultaneously;(2)an enforced framework for online performance assessment considering the existence of process disturbance and outliers is proposed,and the online assessment results contain both steady performance grades and the conversations between performance grades.

    The remainder ofthis paper is organized as follows:Section 2 gives a brief introduction of the algorithms of T-PLS and KPRM,and then the details of the proposed method based on T-KPRM for robust nonlinear process operating performance assessment that are presented in Section 3.The effect of the proposed method is demonstrated by a real case dense-medium coal preparation process,compared with the method proposed by Liu et al.[8].Finally,the paper ends with concluding remarks in Section 5 and some acknowledgements.

    2.Review of T-PLS and KPRM

    2.1.T-PLS algorithm

    T-PLS proposed by Zhou[9]is considered as an improved version of standard PLS.Through further decomposition,the output orthogonal and correlated part can be separated in the principal space,and large process variation can be separated from noise in residual space.Generally,compared with other improved PLS algorithms,e.g.,orthogonal signalcorrection based PLS and orthogonalPLS,T-PLS have been proven to be more effective in providing more accurate process information for those who are more concerned with certain aspects of the whole process information[17].

    Let X∈Rn×mbe the input matrix consisting of n samples with m process variables per sample,and Y∈Rn×vbe the output matrix with v single quality variables per sample.It should be noted that in this paper only the case with single output quality variable is considered,i.e.,the ash of the product coal,and T-PLS algorithm for a single output y∈Rn×lbe the outputmatrix with v single quality variables per sample.It should be noted that in this paper only the case with single output quality variable is considered,i.e.,the ash of the product coal,and T-PLS algorithm for a single output y∈Rn×1is used;for the multiple output case,one just need to treat all of the output variables and use T-PLS algorithm for multiple outputs of Y simultaneously.It is assumed that X and y are centered to zero and scaled to unitvariance.PLSdecomposes X and y as follows[9],

    where T∈Rn×Ais the score matrix,and P ∈ Rm×Aand q ∈ Rp×Aare the loading matrices for X and y,respectively.A is the number ofPLS components,which is usually determined by cross validation.E and f are the residual matrices of X and y,respectively.In the PLS procedure,to obtain the matrix T,the weight matrix W is used.Let R=W(PT W)-1,then T=XR.

    The PLS algorithm extracts the scores T to maximize the covariance between X and y.However,the PLS scores contain components orthogonal to y,which are output-irrelevant parts,and the X-residuals in PLS still contain a large variability,as PLS does not extract the X-variance in descending order.In addition,significant variations in the X-residuals may be affecting the product quality as well.This subspace is not modeled in the PLS model simple because it is not excited in the data that are used to build the PLS model.In summary,variations in both PLS scores and residuals can be relevantto the output y.The T-PLS algorithm is proposed to further decomposes the PLS scores and residuals to sort out components that are relevant to the output y.

    Then,based on the T-PLS algorithm,the X-scores and X-residuals can be further decomposed as follows,

    where ty∈ Rn×1,T0∈ Rn×(A-1)and Tr∈ Rn×Arare three score matrices and Py∈ Rm×1,Po∈ Rm×(A-1)and Pr∈ Rm×Arare the corresponding loading matrices.Eris the new residual matrix after performing PCA on E and can be expressed as Er=E(I-).Aris the number of output unrelated components.For a single output y,T-PLS needs only one score vector tyto predict the quality variable y.

    In Eq.(2),tyrepresents the variations only related to y in original T of PLS model,i.e.,the output-related variations in T,to represent the variations orthogonal to y in T,i.e.,the output-unrelated variations in T,Tris the major part of original X-residual E,and Eris the residual part of E after Tris removed.The T-PLS model decomposes input space X into four subspaces.They are the range spaces of Py,Po,Pr,and()(I-PRT)with dimensionsof1,A-1,Ar,m-A-Ar.A new sample vector xnewis partitioned into four portions,

    and the scores of xnewcan be calculated as follows,

    From the perspective of process operating performance assessment,only the score tyobtained from the of fline modeling and ty,newof online new sample are needed for they represent the output related process variations which can fully re flect the in fluence of process variation on the output,and it is appropriate to evaluate the process operating performance based on them[18].For the detailed procedures for single output T-PLS algorithm please refer to[19].

    2.2.KPRM algorithm

    The partial robust M-regression(PRM),developed from Iterative weighted PLS algorithm,is a robust version of the PLS1 approach[15].Through iterative computation,the algorithm can adaptively assign different weights to different sample data(the outliers are assigned close to 0 weight,the normal data are assigned close to 1 weight).The weighted method is used to overcome the in fluence of outlier on model and retain the part of useful information from outliers.Then,the outliers will be divided into two types by PRM algorithm.One is a high leverage point,which is far from the input data center of the sample points,and the other is the high residual point,which is the sample point of the large difference between the output value and the actual value.The algorithm uses different weighting methods to calculate the two types ofoutliers.The leverage weightofthe i th sample data wixlatent is computed as

    where c=4 represent Euclidean distance,medL1represent the median estimate of L1.tiis the score vector of the i th sample data,c is a tuning constant,taken as f(?).the weightfunction f(?)is called the “Fair”function,and is one of several possible weight functions cited in the original IRPLS.

    The residual weight of the i th sample datais computed as

    where riis the residual part between the actual value and predictive value of the i th sample data.~r is a robust estimator of the sample data.

    Using both types of outliers,we can obtain the comprehensive weights of the sample data,

    the weights of the sample data are updated continuously until the algorithm satis fies the convergence conditions[20].

    KPRM algorithm is developed on the basis of the PRM algorithm.The output data matrix y discussed in this paper is a single output variable.In orderto solve the nonlinearproblemofprocess data,a kernelfunction?is introduced to map the process data from original low dimensional space to high dimensionalfeature space F.In the high dimensionalfeature space F,the relationship of the process variables is linear.After the kernel function mapping,the sample data matrix will be changed to the characteristic data matrixΦ(X)=[?(x1),?(x2),…,?(xN)]T∈RN×M.The dimension M of the feature space can be any value,and even can be set to in finity.De fine kernel matrix K=Φ(X)ΦT(X).When the kernel function is introduced to calculate,we do notneed to know the concrete formofthe nonlinear mapping function Φ(?),and also there is no need for inner product computation in high dimensional space[12].The kernel matrix of high dimensional space can be obtained by calculating the kernel function Ki,jof the original low dimensional space.

    The steps of the KPRM algorithm are as follows:

    Step 1:The input data matrix X is projected to the high dimensional feature space F by nonlinear mapping function Φ(?).

    Step 2:Standard kernel matrix K,then the input kernel matrix K,and output matrix y are processed by PRM algorithm.

    Step 3:Initialize the weight Wiby using the formulas(11),(13),and(15).

    Step 4:Compute weighted input kernel matrix KWand output data matrix yWwhich obtained by weighting the original input kernel matrix X and output data yW.Then the PLS1 regression model is established after the weighted Input and output data,and the score vector is modi fied.

    Step 5:Calculate the residual riof each sample data,and update the weight Wiby using the formulas(11),(13),and(15).

    Step 6:Through consecutive calculation twice,the relative difference between q is less than a certain threshold(10-2),then the convergence condition is satis fied,and go to step 5,Otherwise go back to step 4.

    3.Proposed Approach for Robust Nonlinear Online Operating Optimal Assessment

    3.1.T-KPRM algorithm

    In this paper,the T-KPRM algorithm is proposed by combining the advantages of T-PLS and KPRM.T-KPPM algorithm can further decompose the high dimensionalprincipalcomponentsubspace and the residual subspace of the KPLS model.In the principal component subspace,the output related part and the output independent part will be separated.Meanwhile,in the residual subspace,the larger residual and the final noise will be separated.Thus,T-KPRM can accurately extract the variable information related to the output from process data and overcome the in fluence of outlier by adaptively assigning different weights to each sample data,obtaining exhaustive assessment models for each performance grade.

    Assuming that the process data matrix X∈RN×Jconsists of N samples with J process variables,and the output data matrix y∈RNonly contains one output process variable.The radial basis function is selected as the kernel functionwhere σ is the standard deviation parameter.The steps of the T-KPRM algorithm are summarized as follows,

    Step 1 The columns of the input data matrix X are processed by zero mean and unit variance.Similarly,the output data matrix y is standardized.

    Step 2 The inputdata matrix X is projected to the high dimensionalfeature space F by nonlinear mapping function Ф:xi∈RN→Φ(xi)∈F,and obtain the kernel matrix K,which is standard.

    Step 3 Input kernel matrix K and output matrixyare processed by PRM algorithm.Initialize the weight Wiby using the formulas(11),(13),and(15)and compute weighted input kernel matrix KWand output data matrix yW.Then the PLS1 regression model is established after the weighted KWand yW,and the score vector is modi fied.Calculate the residual riof each sample data,and update the weight Wiby using the formulas(11),(13)and(15)until the convergence condition is satis fied.

    Step 4 The inputdata matrix is KW,the outputdata matrix is yw.Extract convergence of uifromoutputmatrix yW,where i=1,KWi=KW,yWi=yW.(i)Make any column in the yWiequal to ui.(ii)Calculate the score vector of KW,ti=KWiui,ti←ti/‖ti‖.(iii)qi=(iv)Calculate the score vector of yWi,ui=yWiqi,ui←ui/‖ui‖.(v)Determine whether uimeet the convergence condition,if the condition is satis fied,then turn to step 5,otherwise go back to(i);

    Step 5 Calculate the load matrix of KWi,pi=Extract all the main elements,compute T,P,U,Q.(i)Make KWi+1=KWi-yWi+1=yWi-.(ii)Command i=i+1,repeat steps 5 and 6 until the main element A is extracted[4],the number of principal components can be determined by cross validation method.(iii)T=[t1,…,tA],P=[p1,…,pA],U=[u1,…,uA],Q=[q1,…,qA].Then,KW=TPT+E,yW=UQT+f.

    Step 6 The main element TPTis processed by PRM algorithm:TPT=+and the residual E is processed by PRM algorithm:

    If the output data matrix y is a single output variable,the T-KPRM model is as follows,

    where the score tyrepresents the direct parts related to yWin T.tyis also the most important variable in the evaluation model.The score tyrepresents the orthogonal parts related to yW,while Tris part ofa larger variance in E,and Eris the final residual part.

    When a new sample is introduced,the score matrix will be calculated as

    where knewis the kernel function of a new sample data xnewand calculated as

    and the mean centering of knewis

    with 1t=1/n?[1,1,…,1]T∈Rn.

    The score vector tyobtained by the T-KPRMalgorithmcan be used to establish an of fline assessment model of industrial process[7],which contains the information related to the comprehensive economic bene fit.

    3.2.Off-line modeling of performance grade

    Before establishing the evaluation model,the historical data can be divided into severaldata sets according to the comprehensive economic bene fit,and each data setwillrepresentthe corresponding performance grade.For complex industrialprocesses,it is appropriate to set up 3 to 4 performance grades[5].Setting too many performance grades,will make the performance grades very complicated,while setting too few performance grades,will reduce the accuracy of the assessment results.In this paper,four kinds of performance grades,i.e.,poor,general,suboptimal,and optimal,are used and numbered as 1,2,3,and 4 for simplicity,and they correspond to the levelofcomprehensive economic bene fit from low to high respectively.The modeling data sets for each performance grade are denoted as(Xc,yc),c=1,2,3,4,where Xc=[xc(1),xc(2),…,xc(Nc)]T∈RNc×J,yc∈RNc,and Ncis the number of sample data.The of fline evaluation model,and c=1,2,3,4.The score vector form of fline evaluation model calculated by the T-KPRM algorithm is expressed as follows,,i=1,2,…,Nc,where c is the number of performance grades.

    3.3.Enforced online performance assessment strategy

    In the online evaluation part,taking into account that the single sample data is not suf ficient to characterize the production process performance and is vulnerable to interference of disturbances and noises,thus an online sliding data window with width H is introduced as the basic analysis unit[8].The performance-similaritiesbetween the sliding data window and the corresponding performance grade are calculated at the current moment,and the similarity will be used to evaluate the process operating performance.In order to distinguish the deterministic performance grade and the performance grade conversion,the performance-similarity threshold ε(0.5<ε<1)will be used.When the maximum value of the performance-similarityis greater than the thresholdε,the process performance grade can be determined as p.The performance grade conversion is the performancesimilarity change gradually from the previous performance grade to the next one,and the value of the similarity will show a gradual increasing trend.Then,if the similarity is less than the threshold ε and its value meets the trend of increasing gradually,we can determine the process operates on a performance grade conversion process.

    However,despite the sliding window is used,the online assessmentresults maybe inexplicable or even non-actual due to process disturbances and uncertainties.Therefore,the online assessment strategy should be exhaustive and robust to the process disturbances and uncertainties.In this paper,we proposed a modi fied and enforced online operating performance assessment strategy based on three predetermined rules.The main idea is based on the continuity of for the process industries,e.g.chemical process,the industries process operation performance may change fromone performance grade to another performance grade,e.g.change from the optimal performance grade to sub-optimal performance grade.However,the change cannot be finished in an instant,generally,the industrial process performance changes gradually from one grade to another one,the industrial process will go through a conversion performance grade,thus the adjacent online performance assessment results cannot be two different certain performance grades.Ofcourse,we can'tcompletely rule outindividual extreme cases,which rarely appear in the practical industrial process and not be considered in this paper.Generally speaking,we can take all of the assessment results of moment k-1,k and k+1 to determine the ultimate assessment result of the current moment k.

    The enforced online operating performance assessment strategy is summered as follows,

    Step 1 Construct the online sliding data window Xon,k=[xon,k-H+1,…,xon,k]T,where H is the width of data window,for the initial running moment the data window can be constructed by copying the latest sample data.

    Step 2 Normalize xon,kwith the mean and standard deviation of each performance grade.Calculate the standard online window data matrix,c=1,2,3,4 by using of fline modeling.

    Step 3 Extract the scoresfrom,whereis the mean centering of

    Step 4 Calculate the Euclidean distance between the online data window and the corresponding performance grade c bywhereis the mean of scoreaccording to T-KPRM algorithm,and useto compute the similarity between the sliding data window and the corresponding performance grade c,suppose

    where

    Step 5 Then the online assessment result is determined by the following three predetermined rules.

    Rule 1 Firstly,determine whetherthe discriminate formulais satis fied,ifsatis fied,the currentprocess at moment k can be temporarily determined to belong to a certain process performance grade p.then we need to determine whetherthe assessmentresultofmoment k consist with that of moment k-1,if the assessment result of moment k is consist with that moment k-1,it can be if nally determined that the current process is belong to a certain process performance grade p.Otherwise,we need further determine whether the assessment results of moment k+1 and k-1 are the same,ifthe assessmentresult of moment k+1 is the same with moment k-1,it can be determined thatthe assessmentresultofmoment k is unreliable,and the assessment result of moment k should be consistwith those ofmoment k-1and k+1.Ifthe assessment result of moment k+1 is different with that of moment k-1,but consist with that of moment k,it can be determined that the assessment result of moment k and k+1 is reliable,and the performance grade of the process is changing,from certain grade to grade conversion or from grade conversion to certain grade.In the rare instances,the process may jump from one certain grade to another for the existence of strong disturbance,which is not considered in this study.

    Rule 2 If rule 1 is unsatis fied and the condition p=arg maxis satisif ed,it can be temporarily determined that the current process belongs to the performance grade conversion which is the performance-similarity change gradually from the previous performance grade to the next one.l is a positive integer and can be de fined according to production experience.Then we need to determine whether the assessment result of moment k consist with that of moment k-1,ifthe assessmentresult ofmoment k consists with thatmoment k-1,itcan be finally determined that the above assessment result is reliable.Otherwise,we need further determine whether the assessment results of moment k+1 and k-1 are same,if the assessmentresultofmoment k+1 is same with moment k-1,it can be determined that the assessment result of moment k is unreliable,and the assessmentresultof moment k should consist with that of moment k-1.

    Rule 3 If both rule 1 and 2 are unsatis fied due to the in fluence of process disturbances and uncertainties,e.g.the assessment results of moment k-1,k and k+1 are all different,the if nal assessment result of moment k and k+1 is considered unchanged and consist with the previous assessment results of moment k-1.

    Through the above three predetermined rules,we can effectively determine the certain performance grade or the performance grade conversion of the current moment.The enforced online assessment procedure is shown in Fig.1.

    3.4.Non-optimal cause identification

    When the process operating performance is non-optimal,itis meaningfulto find the corresponding causes variables.In the data-based fault diagnosis,the application of the contribution plots[21]method is common.Especially for the linear model,the method can be used to find out the relevant variables which are responsible for fault.However,the traditional contribution plots method is not suitable for the model established by kernel function.Recently,Peng et al.[22]proposed a contribution plots method based on T-KPLS to identify the faulty variables.In this section,the method of non-optimal cause identification based on T-KPRM is used.By calculating the partial derivative of the Euclidean distance which is between the online sliding data window and optimal performance grade,we can obtain the variable contribution,and the variable with the largest contribution is chosen as the cause variable for non-optimal operating performance.

    Fig.1.Procedure of the online assessment and non-optimal cause identification.

    The distancebetween the online sliding data window and optimal performance grade can be re-expressed as follow,

    whereandare the of fline evaluation model matrix and mean centered kernel vector of online sliding data window with respect to optimaloperating performance grade,respectively.Then the contribution of variable j tois calculated as

    whereis the j th variable of the sample data windowrepresents the trace of a matrix.

    4.Simulation Study:Dense Medium Coal Preparation Process

    Dense medium coal preparation process,also known as heavy medium coal preparation process has been widely used to upgrade run-of-mine coal in the modern coal industry by separating gangue from product coal to produces metallurgical coal or power station coal[23,24].In order to ensure the stability of product quality and maximum the comprehensive economic bene fit,the dense medium coal preparation process is expected to under long-term stable optimal operating performance.However,dense medium coal preparation process always suffers from some process disturbances and uncertainties,such as fluctuation of the raw coal property and deviate from the preliminary designed optimal operating point.Thus,developed suitable online operating performance assessment and nonoptimal cause identification methods are crucial for the optimal operation of dense medium coal preparation process.In this study,a dense medium coal preparation process is introduced as the research background,and the actual process data collected from this process is used to illustrate the feasibility and effectiveness of the proposed approach.

    4.1.Plant description

    The coal preparation plant utilizes the heavy medium cyclone to carry out coal sorting,raw coal,and circulating medium are fed into the mixing barrel at the same time,and the mixture is sorted by two products dense medium cyclone,then the over flow and under flow respectively are removed medium,dehydrated and related treatment.Finally,clean coal,middings and gangue products are obtained.While the medium is recovered,concentrated,and ultimately be recycled by medium system.The whole sorting process is divided into the coal flow process and the media flow process.The sensors used in the process include the density meter,the magnetic material,the liquid level meter,and so on.The density meter is used to monitor dielectric barrel suspension density,the liquid level gauge is used to monitor the liquid level of dielectric medium barrel,the magnetic content meter is used to monitor the coal slime content of sorting process.In addition,the whole system with all kinds of solenoid valve controls medium,water,and other flow.The density of heavy medium suspension is the main control object in the process of heavy medium separation.The control of the medium density contains density feedback and ash feedback.The density feedback is used to closed-loop control the density of heavy medium suspension,and the stability of the suspension density is ensured by the control method of the water supply,the diversion,and the adding media.Ash feedback is used to monitor the quality of the product,adjustthe given value ofheavy mediumsuspension density based on ash results and ensure the stability of the product quality[25,26].Fast ash test is the original coal ash indexing which is the basic index of coal quality in coal separating plant.The smaller the ash index is,the less the coal dust is,and the better the quality of coal is.Fig.2 shows the schematic diagram of dense medium coal preparation process.

    4.2.Of fline modeling for performance assessment

    We select five process variables and one output variable for online operating performance assessment of the dense medium coal preparation process.These five process variables are main density/g·L-1,main liquid level/m,main pressure/MPa,magnetic concentration/g·L-1,coal slime content/%,respectively.Then,the one output variable is:ash index of product coal/%.A total of 9900 samples are collected for establishing operating performance model.According to the practical production process,the dense medium coal preparation process is divided into four kinds ofperformance grade and corresponding performance grade conversion.Table 1 lists the division criterion for assessment in dense medium coal preparation process.

    Table 1 The division criterion for assessment

    Fig.2.Process flow diagram of heavy medium coal preparation.

    Itisworth to note thatthe ‘Optimal’in Table 1 isde fined according to the ash index 6.0–6.5 for the realprocess under study,which is different with the optimal de fined by an optimization problem.The ‘Optimal’in the table contains the optimal point and near optimal points,thus the‘Optimal’in Table 1 can also be considered as near optimal.In this section,the historical data of performance grade poor,general,suboptimal and optimal are collected from the actual process and constitute the modeling datasets(X1,y1),(X2,y2),(X3,y3)and(X4,y4)respectively,and the evaluation modelof each stable performance grade is established by the proposed T-KPRM algorithm.

    4.3.Online assessment without outliers

    In online assessment part,another 640 samples including four performance grades and the performance grade conversions are collected to verify the effectiveness of the proposed method.Then,the four performance grades:poor,general,suboptimal,and optimal respectively correspond to the number of test samples 103,230,130,and 112.The three performance grade conversions from poor to general,from general to suboptimal,and from suboptimal to optimal respectively correspond to the number of test samples 23,23,and 19.Set the parameters as H=4,l=1,and ε=0.8.In the simulation,the process operating performances in the test samples are poor→conversion→general→conversion→suboptimal→conversion→optimal.The online assessment results are shown in Fig.3.

    Fig.3(a)–(d)shows the performance-similarities between the online sliding data window and each of the performance grades.For the sample moment 1–103,the performance-similarities between the online data window and performance grade poor are greater than 0.8,which indicating that the current process operating performance is poor based on the assessment rule 1.Similarly,for the sampling moment 127–355,380–509,and 529–640,the performancesimilarities between the online data window and performance grade general,suboptimal,and optimal respectively are greater than 0.8,thus these illustrate that the corresponding process operating performance is general,suboptimal,and optimal separately.

    In addition to the four stable performance grades,at the left sampling moment the process runs into performance grade conversions.For instance,for the sampling moment 104–126,the performancesimilarities between the online data window and performance grade poor are no greater than 0.8,and the successive ascending con-are satis fied.Therefore,it can be identified that the process is operating on the performance grade conversion from poor to general based on the assessment rule 2.Similarly,atthe sampling moment356–378,510–528,itcan be determined that the process operating performance is conversing from general to suboptimal and from suboptimal to optimal respectively.Fig.3(e)shows the assessment results based on the similarities,the vertical axis 1,2,3,4 indicate the performance grade poor,general,suboptimal,and optimal in turn 1.5,2.5,3.5 represent three performance grade conversions respectively.

    4.4.Online assessment with outliers

    Fig.3.Similarities and online assessment results of the proposed method.

    Process data is often subjected to a variety of disturbances(such as sensor failures,process disturbances,etc.),which leads to the existence of outliers.The existence of outliers may make the modeling analysis method lose its proper generalization ability,leading to the decrease ofthe accuracy ofthe assessmentmodels.Thus the robustonline assessmentmethod is needed.In this section,in order to simulate the outliers,we choose to introduce the perturbation to the modeling data sets(X1,y1),(X2,y2),(X3,y3),and(X4,y4).Taking into account that,in the actual production,the number of modeling data sets of optimal and suboptimal performance grades are less than that of poor and general performance grades,5%data of the former modeling data sets are randomly selected to be introduced with disturbances.The way ofintroducing disturbance is that the 5%samples(‘outliers’)selected are divided into two equal parts,and the input of the first part are introduced with 30%disturbance,and the output of the other part are introduced with added 30%disturbance.Meanwhile,10%data of the modeling data of poor and general performance grades are randomly selected to be introduced with disturbance in the similar way as former respectively.The evaluation model( c=1,2,3,4)of each stable performance grade is established by the proposed method(noted as T-KPRM_based).Meanwhile,for comparison,the evaluation model( c=1,2,3,4)of each stable performance grade is also established by the existing method(noted as T-KPLS_based)mentioned in[8].

    Fig.4.Similarities and online assessment results with outliers of the existing method(T-KPLS_based).

    Fig.5.Similarities and online assessment results with outliers of the proposed method(T-KPRM_based).

    Another 603 samples with outliers including four performance grades and the performance grade conversions are collected to verify the robustness of the proposed method.The four performance grades are poor,general,suboptimal,and optimal respectively with the number of test samples 103,211,130,and 96.And the three performance grade conversions are frompoor to general,from generalto suboptimal,from suboptimal to optimal respectively with the number of test samples 23,21,19.Figs.4 and 5 show the evaluation results of the proposed method and the existing method.

    It can be seen that the assessment result of the existing method is not available for actual application due to the effect of outliers.On the other hand,the method proposed in this paper can overcome the effect of outliers and shown better robustness.The assessment result of the proposed method is similar with the assessment results without outliers mentioned in Section 4.3,even though there is time delay owingto the small value of the similarities threshold and use of data window,which is acceptable for actual application.

    Table 2 The assessment identification accuracy rate of the existing method(T-KPLS_based)

    Table 3 The assessment identification accuracy rate of the proposed method(T-KPRM_based)

    For further comparison,Tables 2 and 3 presentthe assessment identification accuracy rate ofthe existing method and the proposed method for certain performance grades respectively.The value of the performance-similarities threshold ε is chosen as 0.6,0.7,and 0.8.It can be seen that the assessment identification accuracy rate of the proposed method is much better than the existing method,which indicates that the proposed method is robust to outliers.It can be also seen from Table 3 that bigger similarities threshold may lead to lower identification accuracy of the proposed performance assessment method for the similarities between the online sliding data window and performance grades emerge obvious fluctuations due to the effect of outliers and noises.In the contrary,smaller similarities threshold may lead to time delay as shown in Fig.5.Therefore,the value of the similarities threshold should be reasonably determined according to actual production process,based on an overall analysis of identification accuracy and time delay.For the case used in this study,we suggest the suitable value of performance-similarities threshold is between 0.6 and 0.8.

    4.5.Non-optimal cause identification

    The non-optimal performance grade means that the industrial process is operating on the other stable performance grade or the performance grade conversion except optimal.The non-optimal cause is usually unchanging when the process operates on a non-optimal performance grade ora performance grade conversion.Thus,the identification result will not be updated until the assessment result changed.

    The corresponding non-optimal cause identification results for all the non-optimal operating performances are shown in Fig.6.The Z axis represents the variable contribution rate,the X axis indicate five process variables main density,main liquid level,main pressure,magnetic concentration and coal slime content and the six sampling moments on the Y axis represent respectively correspond to the nonoptimal operating performances.For example,at 515 sample moment,the variable contribution rate of main pressure and main liquid level is higher than the other process variables,thus the major non-optimal cause is main pressure and main liquid level at this moment.For 472 sample moment,the major non-optimal cause is main pressure and main density.In the same way we can find major non-optimal cause at sample moment 50,115,310,and 365 respectively.In summary,the results of online assessment and non-optimal cause identification can guide the operator timely adjust the control strategy of industrial production,improved the performance of the dense medium coal preparation process,and ensure the product coal quality.

    Fig.6.Non-optimal cause identification results.

    5.Conclusions

    In this paper,we proposed a novel approach for operating performance assessment of complex nonlinear industrial process based on modi fied robust total kernel partial M-regression,T-KPRM.The proposed approach T-KPRM can cope with nonlinearities and outliers simultaneously by combing the advantage of T-PLS and KPRM algorithm,and effective extract the optimal-index-related process variation information from process data to establish assessment models for each performance grade.In the online performance assessment part,we propose an enforced online identification framework taking account of the effects of process disturbances and uncertainties,and the process operational performance can be evaluated accurately based on the predetermined assessment rules.Further,the possible responsible cause for non-optimal operating performance can be found by variable contribution-based identification strategy.Finally,the proposed approach is illustrated with a real industrial case of dense medium coal preparation process,and the results show the efficiency ofthe proposed method comparing to the existing method,which have been confirmed by the actual process data from the dense medium coal preparation process without outliers.

    Acknowledgements

    The authors are grateful to Qi Wu from China University of Mining and Technology for his kind work on the simulation and associate professor Ruhai Lei from China University of Mining and Technology for sharing the real process data used in the case study.

    References

    [1]S.Pashah,A.Moinuddin,S.M.Zubair,Thermal performance and optimization of hyperbolic annular fins under dehumidifying operating conditions–analytical and numerical solutions,Int.J.Refrig.65(2016)42–54.

    [2]P.Wang,C.Yang,X.Tian,et al.,Adaptive nonlinear model predictive control using an on-line support vector regression updating strategy,Chin.J.Chem.Eng.22(7)(2014)774–781.

    [3]Z.Fei,K.Liu,B.Hu,et al.,An efficient latent variable optimization approach with stochastic constraints for complex industrial process,Chin.J.Chem.Eng.23(10)(2015)1670–1678.

    [4]F.Xu,H.Jiang,R.Wang,et al.,In fluence of design margin on operation optimization and control performance of chemical processes,Chin.J.Chem.Eng.22(1)(2014)51–58.

    [5]Y.Liu,F.Wang,Y.Chang,Operating optimality assessment based on optimality related variations and nonoptimal cause identification for industrial processes,J.Process Control 39(2016)11–20.

    [6]L.Ye,Y.Liu,Z.Fei,et al.,Online probabilistic assessment of operating performance based on safety and optimality indices for multimode industrial processes,Ind.Eng.Chem.Res.48(24)(2009)10912–10923.

    [7]Y.Liu,Y.Chang,F.Wang,Online process operating performance assessment and nonoptimal cause identification for industrial processes,J.Process Control 24(10)(2014)1548–1555.

    [8]Y.Liu,Y.Chang,F.Wang,et al.,Complex process operating optimality assessment and nonoptimal cause identification using modi fied total kernel PLS,The 26th Chinese Control and Decision Conference.IEEE 2014,pp.1221–1227.

    [9]D.Zhou,G.Li,S.J.Qin,Total projection to latent structures for process monitoring,AIChE J.56(1)(2010)168–178.

    [10]H.Li,D.Y.Xiao,Survey on data driven fault diagnosis methods,Control Decis.26(1)(2011)1–9.

    [11]Y.Wang,H.Cao,Y.Zhou,et al.,Nonlinear partial least squares regressions for spectral quantitative analysis,Chemom.Intell.Lab.Syst.148(2015)32–50.

    [12]Y.Liu,F.Wang,Y.Chang,et al.,Operating optimality assessment and nonoptimal cause identification for non-Gaussian multimode processes with transitions,Chem.Eng.Sci.137(2015)106–118.

    [13]Y.Xu,X.Deng,Fault detection of multimode non-Gaussian dynamic process using dynamic Bayesian independent component analysis,Neurocomputing 200(2016)70–75.

    [14]Z.Mao,Y.Chang,L.Zhao,Soft-sensor for copper extraction process in cobalt hydrometallurgy based on adaptive hybrid model,Chem.Eng.Res.Des.89(6)(2011)722–728.

    [15]Runda Jia,Zhizhong Mao,Yuqing Chang,Shuning Zhang,Kernel partial robust M-regression as a flexible robust nonlinear modeling technique,Chemom.Intell.Lab.Syst.100(2)(2010)91–98.

    [16]I.Hoffmann,S.Serneels,P.Filzmoser,et al.,Sparse partial robust M regression,Chemom.Intell.Lab.Syst.149(2015)50–59.

    [17]L.I.Gang,Q.I.N.Si-Zhao,J.I.Yin-Dong,et al.,Total PLS based contribution plots for fault diagnosis,Acta Automat.Sin.35(6)(2009)759–765.

    [18]J.J.Lee,C.H.Lim,H.S.Son,et al.,In vitro evaluation of the performance of Korean pulsatile ECLS(T-PLS)using precise quantification of pressure- flow waveforms,ASAIO J.51(5)(2005)604–608.

    [19]S.Yin,X.Zhu,O.Kaynak,Improved PLS focused on key-performance-indicatorrelated fault diagnosis,IEEE Trans.Ind.Electron.62(3)(2015)1651–1658.

    [20]J.González,D.Pena,R.Romera,A robust partial least squares regression method with applications,J.Chemom.23(2)(2009)78–90.

    [21]J.Liu,D.S.Chen,Faultisolation using modi fied contribution plots,Comput.Chem.Eng.61(2014)9–19.

    [22]K.Peng,K.Zhang,G.Li,et al.,Contribution rate plot for nonlinear quality-related fault diagnosis with application to the hot strip mill process,Control.Eng.Pract.21(4)(2013)360–369.

    [23]J.Bai,Y.Xing,Y.Chen,Application of three products heavy medium coal preparation process in Tangshan Chun'ao coal preparation plant,Clean Coal Technol.19(3)(2013)26–29.

    [24]M.Wang,L.Tang,The introduction and selection of the wear-resistant material used for the heavy medium pipeline in coal preparation plant,Coal Qual.Technol.1(2011)69–72.

    [25]D.Dou,J.Yang,J.Liu,et al.,A novel distribution rate predicting method of dense medium cyclone in theTaixi coal preparation plant,Int.J.Miner.Process.142(2015)51–55.

    [26]J.Chen,K.W.Chu,R.P.Zou,et al.,Prediction of the performance of dense medium cyclones in coal preparation,Miner.Eng.31(5)(2012)59–70.

    亚洲精品456在线播放app| 纯流量卡能插随身wifi吗| 久久综合国产亚洲精品| 精品午夜福利在线看| 中文字幕另类日韩欧美亚洲嫩草| 久久精品aⅴ一区二区三区四区 | 99国产综合亚洲精品| 久久精品aⅴ一区二区三区四区 | 啦啦啦在线观看免费高清www| 久久国产精品大桥未久av| 国产一区有黄有色的免费视频| 老司机亚洲免费影院| 激情五月婷婷亚洲| 亚洲天堂av无毛| 精品国产一区二区三区久久久樱花| 成人综合一区亚洲| 在线免费观看不下载黄p国产| 天堂8中文在线网| 国产成人91sexporn| 大话2 男鬼变身卡| 久久鲁丝午夜福利片| 亚洲性久久影院| 国产精品.久久久| 亚洲第一区二区三区不卡| av.在线天堂| 国产极品天堂在线| 哪个播放器可以免费观看大片| 亚洲国产精品专区欧美| 99久久中文字幕三级久久日本| 你懂的网址亚洲精品在线观看| 少妇的逼水好多| 欧美日韩亚洲高清精品| 国产精品蜜桃在线观看| 熟妇人妻不卡中文字幕| 午夜福利影视在线免费观看| 日本猛色少妇xxxxx猛交久久| 男人操女人黄网站| 日本黄大片高清| 午夜福利,免费看| 97在线人人人人妻| 欧美丝袜亚洲另类| 我要看黄色一级片免费的| 男女高潮啪啪啪动态图| 免费在线观看黄色视频的| 制服丝袜香蕉在线| 国产成人一区二区在线| 免费黄频网站在线观看国产| 国产日韩欧美视频二区| 9191精品国产免费久久| 午夜久久久在线观看| 午夜av观看不卡| 成人毛片60女人毛片免费| 午夜免费男女啪啪视频观看| 国产乱人偷精品视频| 一边摸一边做爽爽视频免费| 国语对白做爰xxxⅹ性视频网站| 亚洲国产精品一区三区| 国产男人的电影天堂91| 最近最新中文字幕免费大全7| 国产精品偷伦视频观看了| 久久久久国产网址| 国产毛片在线视频| 97超碰精品成人国产| 精品一区二区三区视频在线| 久久 成人 亚洲| 中文字幕最新亚洲高清| 国产亚洲最大av| 我的女老师完整版在线观看| 街头女战士在线观看网站| 肉色欧美久久久久久久蜜桃| 亚洲国产日韩一区二区| 久久毛片免费看一区二区三区| 精品亚洲乱码少妇综合久久| 三级国产精品片| 亚洲精品日本国产第一区| 日日摸夜夜添夜夜爱| 亚洲精品乱久久久久久| 伊人亚洲综合成人网| 日本-黄色视频高清免费观看| 九草在线视频观看| 亚洲性久久影院| kizo精华| 国产成人精品一,二区| 免费黄网站久久成人精品| 国产精品熟女久久久久浪| 亚洲成人手机| 久久午夜综合久久蜜桃| 国产成人免费观看mmmm| 一区在线观看完整版| 美女大奶头黄色视频| 国产av精品麻豆| 赤兔流量卡办理| h视频一区二区三区| videos熟女内射| 美女视频免费永久观看网站| 亚洲五月色婷婷综合| 亚洲精品成人av观看孕妇| 国产精品久久久久成人av| 国产av国产精品国产| 国产视频首页在线观看| 黑人猛操日本美女一级片| 男女边吃奶边做爰视频| 久久人人爽av亚洲精品天堂| 免费av中文字幕在线| www.熟女人妻精品国产 | 一二三四在线观看免费中文在 | 精品一区二区三卡| 老女人水多毛片| 丁香六月天网| 自线自在国产av| 免费看av在线观看网站| 国产黄色视频一区二区在线观看| 国产精品熟女久久久久浪| 好男人视频免费观看在线| 啦啦啦视频在线资源免费观看| 久久影院123| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲av成人精品一二三区| 欧美精品一区二区免费开放| 欧美精品一区二区大全| 亚洲欧美一区二区三区黑人 | 国产精品人妻久久久久久| 日韩视频在线欧美| 国产日韩欧美在线精品| 精品亚洲成a人片在线观看| 国产爽快片一区二区三区| 亚洲精品中文字幕在线视频| 一边摸一边做爽爽视频免费| 女性被躁到高潮视频| 久久国内精品自在自线图片| 嫩草影院入口| 美女福利国产在线| 91精品国产国语对白视频| 视频中文字幕在线观看| 青春草亚洲视频在线观看| 亚洲精品国产色婷婷电影| 美女xxoo啪啪120秒动态图| 丰满乱子伦码专区| 精品第一国产精品| 亚洲欧洲精品一区二区精品久久久 | 免费人妻精品一区二区三区视频| 成人毛片60女人毛片免费| 制服丝袜香蕉在线| 天堂中文最新版在线下载| 成人亚洲欧美一区二区av| 成人亚洲欧美一区二区av| 五月开心婷婷网| 久久鲁丝午夜福利片| 欧美最新免费一区二区三区| 日韩av在线免费看完整版不卡| 视频在线观看一区二区三区| 亚洲精品久久久久久婷婷小说| 亚洲精品视频女| 最新的欧美精品一区二区| 国产精品秋霞免费鲁丝片| 多毛熟女@视频| 国国产精品蜜臀av免费| 青春草亚洲视频在线观看| 一区二区三区四区激情视频| 美国免费a级毛片| 侵犯人妻中文字幕一二三四区| 欧美3d第一页| 亚洲国产欧美日韩在线播放| 久久久久久人妻| 免费黄频网站在线观看国产| 日韩电影二区| a级片在线免费高清观看视频| 岛国毛片在线播放| 国产av精品麻豆| 1024视频免费在线观看| 亚洲精品自拍成人| 99热国产这里只有精品6| 少妇的丰满在线观看| 国产高清国产精品国产三级| 中国国产av一级| 美女内射精品一级片tv| 国语对白做爰xxxⅹ性视频网站| 日韩 亚洲 欧美在线| 国产无遮挡羞羞视频在线观看| 久久久久久久大尺度免费视频| 久久国产精品男人的天堂亚洲 | 欧美亚洲 丝袜 人妻 在线| 街头女战士在线观看网站| 桃花免费在线播放| 日韩av免费高清视频| 中国国产av一级| 国产成人免费无遮挡视频| 又黄又粗又硬又大视频| a 毛片基地| 国产欧美日韩一区二区三区在线| 九九爱精品视频在线观看| 国产亚洲av片在线观看秒播厂| 亚洲精华国产精华液的使用体验| 欧美人与善性xxx| 亚洲成色77777| 亚洲第一区二区三区不卡| 国产精品 国内视频| 只有这里有精品99| 男人操女人黄网站| 一本色道久久久久久精品综合| 91久久精品国产一区二区三区| 精品少妇久久久久久888优播| 91精品国产国语对白视频| 亚洲国产最新在线播放| 97人妻天天添夜夜摸| a级毛色黄片| 亚洲国产看品久久| 国产免费现黄频在线看| 亚洲精品久久久久久婷婷小说| 91午夜精品亚洲一区二区三区| av免费观看日本| 久久久久视频综合| 青青草视频在线视频观看| av女优亚洲男人天堂| 国产精品99久久99久久久不卡 | 亚洲欧美一区二区三区黑人 | 香蕉丝袜av| 成年美女黄网站色视频大全免费| 永久网站在线| 免费人妻精品一区二区三区视频| av福利片在线| 一级毛片 在线播放| 一级a做视频免费观看| 日本av手机在线免费观看| 日韩免费高清中文字幕av| 最近的中文字幕免费完整| 国产成人欧美| 午夜影院在线不卡| 亚洲久久久国产精品| 国产高清不卡午夜福利| 精品福利永久在线观看| 免费人妻精品一区二区三区视频| 精品国产一区二区久久| 啦啦啦在线观看免费高清www| 成人漫画全彩无遮挡| 亚洲成人手机| 国产精品一区www在线观看| 丝袜在线中文字幕| 久久久国产欧美日韩av| 欧美变态另类bdsm刘玥| 亚洲中文av在线| 边亲边吃奶的免费视频| 香蕉国产在线看| 一级毛片 在线播放| 亚洲欧洲精品一区二区精品久久久 | 成人黄色视频免费在线看| 桃花免费在线播放| 狂野欧美激情性xxxx在线观看| av卡一久久| 最近的中文字幕免费完整| 免费黄网站久久成人精品| 少妇被粗大的猛进出69影院 | 哪个播放器可以免费观看大片| 亚洲av.av天堂| 91久久精品国产一区二区三区| 亚洲婷婷狠狠爱综合网| 韩国av在线不卡| 久久久精品94久久精品| 最近最新中文字幕免费大全7| 中国国产av一级| 日韩欧美精品免费久久| 日韩一区二区三区影片| 在线观看人妻少妇| 午夜福利,免费看| 九九在线视频观看精品| 亚洲精品视频女| 国产精品秋霞免费鲁丝片| 国产激情久久老熟女| 亚洲色图综合在线观看| 男人爽女人下面视频在线观看| 欧美3d第一页| 欧美97在线视频| 亚洲国产日韩一区二区| 午夜福利在线观看免费完整高清在| 日本av免费视频播放| 国产一区亚洲一区在线观看| 视频中文字幕在线观看| 老司机亚洲免费影院| 三上悠亚av全集在线观看| 亚洲成色77777| 日本vs欧美在线观看视频| 大香蕉97超碰在线| 亚洲av综合色区一区| 青春草视频在线免费观看| 国产一区二区三区综合在线观看 | 国产成人av激情在线播放| 免费日韩欧美在线观看| 蜜桃国产av成人99| 久久久精品免费免费高清| 99精国产麻豆久久婷婷| 日本av免费视频播放| 欧美国产精品一级二级三级| videossex国产| 18禁动态无遮挡网站| 欧美老熟妇乱子伦牲交| 天天操日日干夜夜撸| 18禁裸乳无遮挡动漫免费视频| 黄色毛片三级朝国网站| 一边摸一边做爽爽视频免费| 免费人成在线观看视频色| 99国产综合亚洲精品| 人人妻人人澡人人看| 精品少妇久久久久久888优播| 国产精品一区www在线观看| 97在线人人人人妻| 亚洲精品美女久久久久99蜜臀 | 久久久久久久久久人人人人人人| 免费高清在线观看视频在线观看| 在线观看三级黄色| 一区二区三区乱码不卡18| 男人舔女人的私密视频| 十八禁网站网址无遮挡| 亚洲激情五月婷婷啪啪| 精品午夜福利在线看| 夜夜骑夜夜射夜夜干| 亚洲图色成人| 性高湖久久久久久久久免费观看| 亚洲天堂av无毛| 最近中文字幕高清免费大全6| 国产无遮挡羞羞视频在线观看| 一区二区日韩欧美中文字幕 | 婷婷色综合大香蕉| 黄片播放在线免费| 妹子高潮喷水视频| 亚洲精华国产精华液的使用体验| 1024视频免费在线观看| 国产国拍精品亚洲av在线观看| 9191精品国产免费久久| 下体分泌物呈黄色| 蜜桃国产av成人99| 久久亚洲国产成人精品v| 亚洲激情五月婷婷啪啪| 久久人人爽人人爽人人片va| 少妇猛男粗大的猛烈进出视频| 99热全是精品| 久久婷婷青草| 在线看a的网站| 亚洲精品乱久久久久久| 日韩欧美精品免费久久| 精品亚洲乱码少妇综合久久| 久久久久国产网址| 69精品国产乱码久久久| 国精品久久久久久国模美| 久久精品夜色国产| 精品一品国产午夜福利视频| 人体艺术视频欧美日本| 搡女人真爽免费视频火全软件| 免费观看性生交大片5| 香蕉精品网在线| 亚洲性久久影院| 免费高清在线观看视频在线观看| 搡老乐熟女国产| 久久人人爽人人片av| 考比视频在线观看| 91午夜精品亚洲一区二区三区| 有码 亚洲区| 久久久久精品性色| 午夜福利乱码中文字幕| 纵有疾风起免费观看全集完整版| 九色成人免费人妻av| 久久久久网色| 成年动漫av网址| 超色免费av| 美女xxoo啪啪120秒动态图| 韩国av在线不卡| 欧美激情 高清一区二区三区| 亚洲精品日韩在线中文字幕| 激情五月婷婷亚洲| 亚洲精品成人av观看孕妇| 成人午夜精彩视频在线观看| 免费高清在线观看视频在线观看| 中文字幕av电影在线播放| 亚洲精品456在线播放app| 久久亚洲国产成人精品v| 成人手机av| 在线观看美女被高潮喷水网站| 999精品在线视频| 91精品三级在线观看| 黄色毛片三级朝国网站| 国产精品蜜桃在线观看| 国产一区二区激情短视频 | 曰老女人黄片| 日韩电影二区| 午夜日本视频在线| 国产精品女同一区二区软件| 妹子高潮喷水视频| 精品少妇黑人巨大在线播放| 免费观看a级毛片全部| 男男h啪啪无遮挡| 男女午夜视频在线观看 | 精品少妇久久久久久888优播| 久久精品久久久久久久性| 熟女av电影| 男人操女人黄网站| 99久久综合免费| 日日啪夜夜爽| 国产熟女欧美一区二区| 大码成人一级视频| 人妻少妇偷人精品九色| 国产又色又爽无遮挡免| 成人手机av| 国产成人免费无遮挡视频| 亚洲少妇的诱惑av| 成年美女黄网站色视频大全免费| 日韩大片免费观看网站| 国产精品久久久久成人av| 国产免费又黄又爽又色| 精品人妻偷拍中文字幕| 美女xxoo啪啪120秒动态图| 黄色 视频免费看| 久久精品人人爽人人爽视色| 人人妻人人澡人人看| 国语对白做爰xxxⅹ性视频网站| tube8黄色片| 高清黄色对白视频在线免费看| 黑人猛操日本美女一级片| 欧美成人午夜精品| 亚洲在久久综合| 日本wwww免费看| 女性被躁到高潮视频| 久久毛片免费看一区二区三区| 久久99蜜桃精品久久| 国产成人免费无遮挡视频| 极品少妇高潮喷水抽搐| 男女下面插进去视频免费观看 | 自线自在国产av| 插逼视频在线观看| 欧美另类一区| 欧美日韩亚洲高清精品| 老司机影院毛片| 午夜激情久久久久久久| 99热这里只有是精品在线观看| 五月开心婷婷网| 在线观看三级黄色| 少妇 在线观看| 免费大片黄手机在线观看| 国产成人91sexporn| 久久久久网色| 伦精品一区二区三区| 极品人妻少妇av视频| 下体分泌物呈黄色| 国产探花极品一区二区| 亚洲av.av天堂| 最近的中文字幕免费完整| 欧美人与善性xxx| 一级毛片 在线播放| 国产精品嫩草影院av在线观看| 午夜免费男女啪啪视频观看| 欧美日韩成人在线一区二区| 一二三四中文在线观看免费高清| 久久久国产欧美日韩av| 精品人妻偷拍中文字幕| 人人澡人人妻人| 黑人巨大精品欧美一区二区蜜桃 | 天美传媒精品一区二区| 乱人伦中国视频| 久久人妻熟女aⅴ| 成人影院久久| 如何舔出高潮| 最近最新中文字幕大全免费视频 | 亚洲av男天堂| 午夜激情av网站| 日本黄色日本黄色录像| 日韩大片免费观看网站| 精品国产一区二区三区久久久樱花| 一区二区av电影网| 欧美精品亚洲一区二区| 黄色配什么色好看| 国产亚洲欧美精品永久| 国产精品久久久久久精品古装| 99视频精品全部免费 在线| 亚洲欧洲精品一区二区精品久久久 | 我要看黄色一级片免费的| 国产伦理片在线播放av一区| av免费在线看不卡| 建设人人有责人人尽责人人享有的| 亚洲中文av在线| 欧美日韩精品成人综合77777| 久久精品国产综合久久久 | 边亲边吃奶的免费视频| 久久久久久久精品精品| 国产免费又黄又爽又色| 五月伊人婷婷丁香| 九色亚洲精品在线播放| 国产精品一区二区在线观看99| 国产片内射在线| 丁香六月天网| 午夜av观看不卡| av女优亚洲男人天堂| 成人免费观看视频高清| 9色porny在线观看| 久久久欧美国产精品| 精品久久国产蜜桃| 亚洲国产看品久久| 26uuu在线亚洲综合色| 免费黄色在线免费观看| 99国产综合亚洲精品| 久久人人爽av亚洲精品天堂| 少妇精品久久久久久久| 欧美人与性动交α欧美精品济南到 | 另类亚洲欧美激情| 亚洲欧洲国产日韩| 99热6这里只有精品| 波野结衣二区三区在线| av天堂久久9| 国产色婷婷99| 搡女人真爽免费视频火全软件| 欧美xxxx性猛交bbbb| 色网站视频免费| 在线观看美女被高潮喷水网站| 亚洲丝袜综合中文字幕| 色婷婷av一区二区三区视频| 久久精品久久久久久久性| 午夜福利在线观看免费完整高清在| 免费观看性生交大片5| 秋霞伦理黄片| 七月丁香在线播放| 丝袜在线中文字幕| 一级片免费观看大全| 午夜福利影视在线免费观看| 久久久久精品人妻al黑| 一级毛片电影观看| 丝袜喷水一区| 亚洲欧洲国产日韩| 欧美成人午夜免费资源| av免费观看日本| 亚洲成人一二三区av| 国产精品国产三级专区第一集| 曰老女人黄片| 一级毛片黄色毛片免费观看视频| 在线观看免费高清a一片| 国产探花极品一区二区| 国产视频首页在线观看| 久久久久久久亚洲中文字幕| 亚洲国产精品专区欧美| 晚上一个人看的免费电影| 一级,二级,三级黄色视频| 国产亚洲精品久久久com| 精品酒店卫生间| 亚洲 欧美一区二区三区| 久久久久久久久久成人| www.熟女人妻精品国产 | 欧美日韩视频精品一区| 人妻系列 视频| 精品福利永久在线观看| 春色校园在线视频观看| 最近中文字幕2019免费版| 中文字幕av电影在线播放| a级片在线免费高清观看视频| 大码成人一级视频| 97在线人人人人妻| 久久精品久久久久久噜噜老黄| 欧美亚洲日本最大视频资源| 汤姆久久久久久久影院中文字幕| a级毛片黄视频| 久久人妻熟女aⅴ| 人妻人人澡人人爽人人| 亚洲经典国产精华液单| 国产成人免费无遮挡视频| 99re6热这里在线精品视频| 久久久久人妻精品一区果冻| 男女无遮挡免费网站观看| 99久久中文字幕三级久久日本| 亚洲精品国产av蜜桃| av在线观看视频网站免费| 侵犯人妻中文字幕一二三四区| 国产福利在线免费观看视频| 亚洲综合精品二区| 香蕉丝袜av| 美女内射精品一级片tv| 国产男女超爽视频在线观看| 国产亚洲精品久久久com| 中文字幕另类日韩欧美亚洲嫩草| 精品国产国语对白av| 亚洲成人av在线免费| 亚洲欧美清纯卡通| 一级毛片黄色毛片免费观看视频| 一本久久精品| 亚洲av日韩在线播放| 久久久久精品性色| a级毛片在线看网站| 97人妻天天添夜夜摸| 男女边摸边吃奶| 少妇熟女欧美另类| 亚洲精品日韩在线中文字幕| 日韩 亚洲 欧美在线| 国产免费一级a男人的天堂| av不卡在线播放| 在线看a的网站| 曰老女人黄片| 交换朋友夫妻互换小说| 久久国内精品自在自线图片| 人体艺术视频欧美日本| 亚洲 欧美一区二区三区| 五月开心婷婷网| 王馨瑶露胸无遮挡在线观看| 国产精品久久久久久久电影| 亚洲国产最新在线播放| 日本-黄色视频高清免费观看| 国产精品熟女久久久久浪| 最近中文字幕2019免费版| 嫩草影院入口| 最近中文字幕2019免费版| 99久久精品国产国产毛片| 国产免费一级a男人的天堂| 大香蕉97超碰在线| 边亲边吃奶的免费视频| 人妻少妇偷人精品九色| 国产精品 国内视频| 少妇 在线观看| 免费高清在线观看视频在线观看| 性色av一级| 波多野结衣一区麻豆| 久热久热在线精品观看| 免费播放大片免费观看视频在线观看|