999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于改進原子軌道搜索算法的多工藝路線柔性作業車間問題研究

2024-01-31 07:04:10李佳蓉晁永生李純艷袁逸萍
機床與液壓 2024年1期
關鍵詞:模型

李佳蓉,晁永生,李純艷,袁逸萍

(新疆大學智能制造現代產業學院,新疆烏魯木齊 830017)

0 前言

在“中國制造2025” 的背景下,我國制造產業的生產模式正逐步由大批量單一產品向多品種、小批量定制化產品方向轉變。定制化產品在生產過程中,由于其結構復雜、多變等原因,通常有多條不同的生產加工路線。同時,日益多樣化的訂單需求,迫使制造單元更具柔性。

在傳統作業車間問題中,各工件的工序、加工機器以及對應的加工時間是固定的,而柔性作業車間問題(Flexible Job Shop Problem,FJSP)是在作業車間問題的基礎上兼顧考慮工件各工序機器選擇的不確定以及機器加工不同工序的時間不確定,這更符合實際生產的情況。

近年來,FJSP 已被國內外諸多學者廣泛研究,如,在模型搭建和求解方面,JALILVAND-NEJAD、FATTAHI[1]基于加工循環下的FJSP 問題提出一種混合整數線性規劃模型,并采用遺傳算法進行求解驗證。WU 等[2]研究了考慮夾具裝卸時間的雙資源約束FJSP 并建立了數學模型。VITAL-SOTO 等[3]建立了具有排序靈活性的FJSP 的數學模型,并提出一種混合細菌覓食優化算法求解該問題。PENG 等[4]構建了以加工時間、能耗和噪聲為優化目標的雙約束多目標的FJSP 模型,并提出一種混合離散多目標帝國競爭算法來進行求解。NOURI 等[5]提出一種基于混合元啟發式算法的聚類整體多智能體方法來求解FJSP。張宇嘉、宋威[6]以完工時間最優為目標,采用構建精英檔案和進步檔案的方法,提出一種雙檔案粒子群算法對FJSP 進行求解。雖然他們從模型和求解方法上對FJSP 進行了深入研究,但目前絕大多數的研究只考慮了工件的單個工藝方案,即工件的每道工序只考慮機器的可替代性。而在實際生產過程中工件可以有多條工藝路線,前人對考慮多工藝路線FJSP 的研究相當有限。

在柔性作業車間的綜合調度過程中,往往忽略多工藝路線這一關鍵因素,造成設備利用率較低、零件在某一機器堆積、單件產品加工時間過長等問題。因此,將多工藝路線下的柔性作業車間調度應用于產品制造對提高企業綜合競爭能力尤為重要。同時,由于加入多約束,模型在求解時的搜索空間范圍增大,一般算法容易在求解過程中陷入局部最優,導致搜索魯棒性差、效率低下。因此,本文作者建立一種多工藝路線柔性作業車間問題(Multiprocess Route Flexible Job Shop Problem,MRFJSP)模型,并提出一種改進原子軌道搜索(Improved Atomic Orbital Search,IAOS)算法,用于避免模型求解過程中效率低下的問題。

1 多工藝路線柔性作業車間調度問題

1.1 問題描述

MRFJSP 模型作為企業在生產過程中為各工件選擇合理加工路線、工序排序和機器分配,以實現調度目標的有效手段,可描述如下:有多個工件需要加工,每個工件的加工特征不同,不同特征可由多種加工方式實現,不同的加工方式可組合出多條加工路線,各路線包含的工序數不同,各工序可選的機器柔性,各機器加工同道工序的時間不同。由于在生產調度過程中會出現一種產品可按照不同工藝路線進行加工、一臺機器同一時刻只能加工一件產品等情況,為明確MRFJSP 模型中的約束關系,做出如下假設:

(1)零件可選擇不同工藝路線且不同工藝路線下的工序之間不存在優先級關系;

(2)禁止關閉閑置機器且不考慮機器故障;

(3)工序一旦開始就不允許中斷;

(4)零件的一道工序完成后即可到達加工的下一個機器。

文中模型相關參數及其含義見表1。

表1 相關參數及其含義Tab.1 Relevant parameters and their meanings

1.2 問題模型搭建

為提高多工藝路線柔性作業車間效率,以完工時間最小為調度優化目標函數,表示為式(1):

根據基準先行、先主后次、先面后孔的工藝規則,在保證有效和完整的前提下滿足以下約束:

(1)工藝路線約束。每個工件每次只能選擇一條工藝路線,表示為式(2):

(2)機器選擇約束。各工序只允許被一臺機器進行加工,表示為式(3):

(3)機器約束。在任意時刻每臺機器只能加工一道工序,表示為式(4):

(4)先序關系約束。各工件的不同工序不允許同時加工,表示為式(5):

(5)工件加工時間約束。加工工件的完成時間非負,表示為式(6):

(6)工序時間約束。工序的開始時間均不大于工序的完工時間,表示為式(7):

(7)加工路線決策變量約束。工件在第l條加工路線被選中,表示為式(8):

(8)機器決策變量約束。工件的第l條路線的第j道工序是否在機床k上加工,表示為式(9):

2 改進原子軌道搜索算法

原子軌道搜索(Atomic Orbital Search,AOS)算法是由AZIZI[7]在2021 年提出的一種高效的元啟發式算法,其設計源于量子力學原理和量子原子模型。AZIZI 等[8]結合AOS 研究了多個領域著名工程問題,驗證了AOS 適用于求解不同規模組合優化問題,具有較優的魯棒性、群體性。但考慮加入加工路線柔性和機器柔性等約束,問題的復雜性指數式增長,搜索空間急速擴大,此時若使用AOS 進行求解,反而會降低算法的搜索效率,影響全局最優值的獲得,因而本文作者提出一種改進原子軌道搜索(Improved Atomic Orbital Search,IAOS)算法。

2.1 原子軌道搜索算法基本原理

AOS 算法通過效仿原子核外部電子密度變化、電子能量狀態以及受各類因素干擾電子發生躍遷,從電子的選擇、搜索和更新這3 個方面建立數學模型,模擬算法的尋優過程[7]。

AOS 將原子核外的電子云模擬成薄的、同心的層作為算法的搜索空間,如圖1 所示。電子層表示為Xk,k表示電子層數。搜索空間中的電子用候選解Xi表示,候選解為電子在原子軌道模型中的映射,Xi的位置用一組決策變量[,,…,]表示,電子層Xk中的一組候選解如式(10)所示。候選解在搜索空間中的層數是由電子的概率密度確定的,在數學模型中使用概率密度函數進行確定。

圖1 電子躍遷示意Fig.1 Electron transition schematics

候選解在搜索空間中的初始位置由式(11)隨機確定,候選解的能量值表示為,即函數值。在求解最小化優化問題時,將候選解的能量值升序排序,能量值越小,候選解的概率密度值越高,以代表具有較低能量值的電子。

其中:i表示電子層中的第i個候選解,i∈{1,2,…,m};d表示維度,j∈{1,2,…,d}。

候選解進行位置更新時需考慮光或其他影響因素,如粒子、磁場等對電子的作用。電子受這些因素影響時會從基態躍遷至激發態,離開所在電子層,電子位置隨之變化,即候選解位置發生更新。候選解位置更新的方式如下:

(1)當光對電子產生作用(?≥PR)時:

①當電子吸收光子時(≥BE),電子能量升高,電子從低能級向高能級躍遷,電子位置更新,即候選解位置更新,可用式(12)表示:

②當電子發射光子時(<BE),電子能量降低,電子從高能級向低能級躍遷,電子位置更新,即候選解位置更新,可用式(13)表示:

(2)除了受光子的影響外,電子的能量還受其他因素作用(?<PR),此時電子位置也會更新,可用式(14)表示:

其中:?表示電子受影響的概率;PR表示光對電子作用的概率;分別表示更新前、后的候選解位置;LE表示所在層最低能量的電子位置;BS表示候選解所對應的結合態,用平均候選解位置向量表示;BE表示候選解所在層所對應的結合能,用層最低能量表示;ri表示一個隨機生成的電子位置;αi、βi和γi均為隨機數。

2.2 改進算法的編碼與解碼

基于建立的MRFJSP 模型,將候選解位置向量設為三段式向量,三段向量分別對應路線選擇、工序排序與機器分配。候選解位置向量維度d=N+2G,第一段向量長度為N,第二段、第三段向量長度均為G,N為工件數,G為加工所有工件的總工序數。候選解位置向量中各元素值均使用隨機函數均勻生成實現全局初始化,候選解位置向量中各元素取值范圍均為[-ε,ε],ε=N。

假設有一個MRFJSP,問題中有2 個工件,每個工件分別有2 條加工路線,每個工件的1、2 條路線分別有2、3 道工序。該問題的候選解由3 段向量構成,第一段向量長為2,后兩段的長度需根據第一段向量中的信息進行計算。若候選解第一段隨機位置向量為[1.8,-1.4],則通過編碼操作得到第一段向量為[2,1],即工件1 選擇路線2,工件2 選擇路線1。此時可知該加工任務的總工序數為5,則該候選解位置向量總長為12,圖2 為滿足第一段位置向量[1.8,-1.4]的一個候選解。

圖2 候選解位置向量Fig.2 Candidate solution position vector

2.2.1 編碼

候選解位置向量的編碼包含路線選擇、工序排序和機器分配3 個部分。對第一段路線選擇位置向量進行編碼時,編碼的長度對應總工件數,向量中的每個位置依次與工件號對應。根據各工件的可選路線數將候選解位置元素的取值范圍[-ε,ε]均勻分區,判斷候選解位置元素值所屬的分區,得到候選解位置元素對應的路線。以圖2 為例,候選解位置元素的取值范圍是[-2,2],工件1 的可選路線兩條,將范圍均勻分區為[-2,0)和[0,2],判斷元素所屬區間,則工件1 選擇路線2,工件2 選擇路線1。完成路線選擇位置向量的編碼、獲得總工序數后,采用最大位置規則對工序排序位置向量進行編碼[9],首先對工序排序位置向量進行降序排序,再依照所得的排序位置對工序順序進行重排列,得到編碼結果。編碼轉換過程如圖3所示。

圖3 工序排序位置向量編碼過程Fig.3 Process sorting position vector encoding process

完成前兩段編碼后才可開始第三段機器分配位置向量編碼。在該編碼過程中需考慮候選解位置元素與對應工件的可選機器之間的關系。文中采用文獻[10]的編碼方式,如公式(15),它將候選解位置向量映射到對應可選機器的機器序號,實現機器分配編碼,即公式對應計算出的值就是選擇的機器號,當計算的機器號不存在時,設置該機器號為1。

其中:為機器分配編碼;為候選解位置元素;s(d)為可選擇的機器集合;ε為總工件數。

由式(15)可對機器分配位置向量進行編碼,編碼的轉換過程如圖4 所示。

圖4 機器分配位置向量編碼過程Fig.4 Machine allocation position vector encoding process

2.2.2 解碼

解碼也包含路線選擇、工序排序和機器分配3 個部分。第一段路線選擇解碼需根據各工件的可選路線數將候選解位置元素的取值范圍[-ε,ε]均勻分區,判斷候選解編碼元素值所屬的分區范圍,在該范圍內隨機取值,從而實現路線選擇解碼。

完成第一段解碼后采用最大位置規則中的轉換公式(16),獲得第二段工序排序位置向量對應的解碼向量[9],然后將原工序排序位置向量進行升序排序,根據升序排序后的位置值對解碼向量重新排序,最終完成工序排序位置向量解碼。

完成第二段解碼后對第三段機器分配向量解碼,該解碼過程相當于機器分配編碼過程的逆轉換,轉換公式如式(17)所示,通過該公式可直接計算出對應的機器分配位置向量。若計算時出現可選機器數為1,則隨機生成該位置的解碼元素。

2.2.3 候選解位置向量的增碼與減碼

求解時候選解位置向量元素值會隨位置更新變化,當路線選擇位置向量中的某個元素值改變時,可能會導致工件加工路線變化,任務的總工序數隨之變化,候選解維度改變。當總工序數變化時,原始編碼長度將不符合編碼規則,無法求解,所以在解碼時要考慮候選解位置向量的增碼與減碼。

增碼與減碼的操作方法為:每次編碼前測量候選解原始維度對應的工序總數以及候選解路線選擇段對應的實際工序總數。若原始工序總數小于實際工序總數,則進行補碼,補碼數為實際工序總數與原始工序總數之差,使用隨機函數生成相應數量的位置元素補碼。反之,采用位置向量末端減碼,減碼數為實際工序總數與原始工序總數之差。

2.3 自體交叉

在作業車間調度方案中,由于工序在不同機器上的排布順序不同,完成加工任務的時間也不同,所以工件在機器上的排布往往有較多空閑時間,造成時間上的浪費。因此可以調整工序在機器上的排布,實現完工時間的優化。

自體交叉的操作過程如圖5 所示。自體交叉操作通過選取候選解單個個體工序排序段上的兩個隨機位置進行元素值互換,搜索排布中的空閑時間進行插補,縮短完工時間。

圖5 候選解自體交叉過程Fig.5 Candidate solution self-crossover process

2.4 變鄰域變異

作業車間調度問題中的鄰域結構直接影響局部最優解、全局最優解的獲得。變鄰域搜索不但能避免大量無效搜索過程,還能擴展當前的搜索空間,避免陷入局部最優,所以對MRFJSP 采用變鄰域搜索方法對候選解進行變異搜索。根據變鄰域搜索方法中鄰域范圍特性(范圍介于[2,4]時搜索效率和穩定性較好,范圍上限超過5 時搜索效率逐漸降低)[11],結合問題需要,合理設置鄰域范圍為[2,4]。

變鄰域變異的操作過程。先判斷候選解是否發生變鄰域變異;若發生變異,在鄰域范圍內隨機搜索得到一組工件位置;對工件位置數據全排列獲得多個可行候選解;計算出這些候選解的能量值再與原候選解進行比較;將最優候選解作為變鄰域變異后的新候選解。變鄰域變異的操作過程如圖6所示。

圖6 候選解變鄰域變異過程Fig.6 Candidate solution variation neighborhood mutation process

2.5 變工序數候選解精英保留策略

由于MRFJSP 中各工件有多條加工路線且生成的候選解維度不同,需考慮維度變化對搜索空間的影響,簡單地保留能量值較優的個體是不合理的。為保證子代搜索空間中候選解的多樣性以及避免算法陷入局部搜索空間,提出一種變工序數候選解精英保留策略。在算法生成初始候選解時對候選解標記維度,在候選解進行位置更新、自體交叉和變鄰域變異時更新維度標簽,完成每次迭代時保留能量值較優的候選解和各維度標簽中最優候選解,作為下一次迭代的初始搜索空間。

2.6 算法流程

具體算法流程如圖7 所示。

圖7 改進原子軌道搜索算法流程Fig.7 Improved atomic orbital search algorithm process

改進原子軌道搜索算法首先對候選解位置進行初始化并計算能量值,然后確定當前原子的結合態、結合能及最低能級,并使用概率密度函數分配候選解到各虛擬層,計算虛擬層的結合態、結合能及層能級,更新對應層候選解位置,使用自體交叉、變鄰域變異策略進行增強搜索,通過變工序數精英保留策略保留最優個體,重復上述過程直至滿足結束條件。

3 案例驗證與分析

本文作者在Windows 10、CPU 2.40 GHz、內存8 GB 的系統上采用MATLAB 2020b 編程實現IAOS 算法運行。IAOS 算法參數設置為:迭代次數為300,候選解個數為200,原子軌道數為5;光子對電子的作用率0.1;交叉概率0.8;變異概率0.2。

為能夠合理驗證前期建立的MRFJSP 模型,選取3 種規模的實例數據驗證,以文獻[12]中某內燃機車柔性生產作業車間的部分數據作為基礎數據[13]。這3 組 數據 分別 對應3、6、10 工 件在小、中、大規模問題上的加工任務且工藝路線可選。其中六工件六機器中規模問題的工藝路線信息如表2所示。

表2 中規模問題的工藝路線信息Tab.2 Process route information for medium scale problems

表中Ril和Oij分別表示工件的可選工藝路線和工序,i(i=1,2,…,6)表示工件號,l(l=1,2,…,6)表示路線號,j(j=1,2,…,6)表示工序號。表2 中工序Oij對應的可選機器和加工時間如表3 所示。

表3 表2 中工序Oij對應的加工信息Tab.3 Processing information corresponding to process Oij in Tab.2

3.1 模型求解

在FJSP 中各工件的加工路線固定、機器選擇柔性,無需考慮多個加工路線。在MRFJSP 中需要考慮加工路線柔性和機器選擇柔性,兩個模型的求解目標均為完工時間最小化。為避免隨機性,采用3 種不同規模問題對FJSP 模型和MRFJSP 模型進行測試,各運行10 次。中規模問題中FJSP 模型選擇各工件的原始車間加工路線進行求解,即參照表2 中的Ri1所對應的6 條加工路線,其中i=1,2,…,6。

經典FJSP 模型和MRFJSP 模型對3 種不同規模問題數據對進行求解,解得的最優調度結果如表4 所示,tbest、tavg和t分別表示最優完工時間、平均完工時間和算法求解的平均運行時間。兩模型中規模問題調度甘特圖如圖8、9 所示。

圖8 FJSP 中規模問題的最優調度甘特圖Fig.8 The optimal scheduling Gantt chart for scale problems in FJSP

表4 模型對比實驗結果Tab.4 Model comparison test results

3.2 求解結果對比分析

由表4 中兩模型實驗結果可看出,對于不同規模問題MRFJSP 模型求解得到的完工時間總體小于FJSP 模型。從小、中、大3 個規模問題上進行分析,使用MRFJSP 模型相較于FJSP 模型在完工時間上分別縮短了25.95%、25.6%和29.53%。對比圖8 和圖9 所示的最優調度方案甘特圖可看出MRFJSP 模型的調度方案相較于FJSP 模型的排程時間更緊湊,排布更合理。結合實驗結果和甘特圖綜合分析驗證了MRFJSP 模型求解方法能有效求解多工藝路線的FJSP,相對于經典的FJSP,文中提出的模型能有效縮短完工時間,提高生產效率。

圖9 MRFJSP 中規模問題的最優調度甘特圖Fig.9 The optimal scheduling Gantt chart for scale problems in MRFJSP

3.3 算法評價

將IAOS 與其他的元啟發式算法進行實驗對比,對比算法包括原子軌道搜索算法(AOS)[7]、遺傳算法(GA)[14]、變鄰域遺傳算法(GA-VNS)[12]。為避免隨機性,算法各運行10 次,參考文獻[15]的對比方法加入提升率R,用來表示IAOS 相較于其他元啟發式算法的提升水平,其計算公式為(tIAOS,Best-tAOS,Best)/tIAOS,Best。各算法在不同規模下的求解結果如表5 所示,各算法求解中規模問題的迭代曲線如圖10 所示。

圖10 中規模問題算法迭代曲線Fig.10 Algorithm iteration curves for medium scale problems

表5 算法實驗結果Tab.5 Algorithm experiment results

由表5 可知:從求解效果來看,對于小規模問題,除GA 算法的求解結果不佳外,其他3 種算法都取得了較好的解,IAOS 算法相較于AOS、GA 和GA-VNS算法的平均提升率分別提升了1.23%、1.64%和0.81%;對于中規模問題的求解,IAOS 搜索到的結果最好且與其他搜索算法的最優解拉開了一定的差距,相較于AOS、GA 和GA-VNS 算法的平均提升率分別提升了16.2%、27.12%和13.79%,由此體現出IAOS 改進的有效性以及在求解多工藝路線FJSP 上的優越性;對于大規模問題,AOS 和GA 算法的求解結果不理想,求解能力明顯下降,GA-VNS算法雖然在大規模問題上仍有一定的求解能力,但求解結果不如IAOS 優秀,證明了IAOS 算法在求解復雜度較高問題時的高效性。

由圖10 可知:對于中規模問題各算法在運算初期就能快速尋優,但在迭代50 次后明顯能看出AOS、GA 和GA-VNS 算法尋優能力均有所下降,IAOS 算法在迭代50 次后雖尋優速度減慢,但仍然具有一定的搜索能力。從收斂性上分析,IAOS、AOS、GA 和GA-VNS 算法分別在迭代175、38、57 和206 次時收斂。IAOS 雖然收斂速度較慢,但求得了優于其他算法的解;AOS 收斂最快,獲得了比較好的求解結果;GA 收斂也較快,但求得的結果最差;GA-VNS 收斂速度最慢,但其對GA 的改進有效,也取得了比較好的解。綜上可看出:雖然IAOS 收斂速度慢于AOS 和GA,時間復雜度略高于其他算法,但IAOS 的改進措施有效,在算法搜索的開始階段就能快速找到優于其他算法的解,能有效跳出局部最優,找到更好的解。

綜合來看,IAOS 算法在整體求解過程中尋優能力較好,應用于不同規模問題時的適用性較強,尤其是在求解復雜度較高的大規模問題時IAOS 算法比其他算法表現出更好的尋優能力。因此采用IAOS 算法進行多工藝路線柔性作業車間調度優化是可行的。

4 結論

(1)本文作者分析了工件的多工藝路線與柔性作業車間調度之間的聯系以及特點,建立了多工藝路線柔性作業車間調度模型,通過實驗驗證了多工藝路線柔性作業車間調度模型的有效性。

(2)設計了一種改進原子軌道搜索算法,引入一種三層編碼方式使其適用于求解MRFJSP;引入自體交叉和變鄰域變異增強局部搜索并避免陷入局部最優;迭代過程中設計了變工序數精英保留策略,擴大了搜索空間,使算法的求解質量和效率得到提高。

(3)采用3 組不同規模MRFJSP 對改進原子軌道搜索算法和其他3 種元啟發式算法的求解結果進行比較分析,驗證了改進原子軌道搜索算法的可行性和適用性。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 日韩A∨精品日韩精品无码| 91无码视频在线观看| 国内99精品激情视频精品| 久久亚洲黄色视频| 秋霞午夜国产精品成人片| 欧美一区二区三区不卡免费| 久久无码av三级| 欧美一区二区福利视频| 午夜一区二区三区| 国产91特黄特色A级毛片| 午夜视频在线观看区二区| 2020最新国产精品视频| 日本不卡在线| 乱人伦中文视频在线观看免费| 日韩欧美色综合| 特级aaaaaaaaa毛片免费视频 | 在线视频亚洲欧美| 超薄丝袜足j国产在线视频| 女同久久精品国产99国| 欧美一区福利| 国产激爽大片在线播放| 综合久久五月天| 亚洲AV无码久久精品色欲| 国产精品无码一二三视频| 亚洲国语自产一区第二页| 欧美19综合中文字幕| 亚洲国产精品不卡在线| 亚洲欧美一区在线| 看国产一级毛片| 欧美成人精品一级在线观看| 一本大道在线一本久道| 国产精品55夜色66夜色| 国产jizz| 亚洲一区无码在线| 国产视频一二三区| 国模沟沟一区二区三区| 国产另类乱子伦精品免费女| 无码又爽又刺激的高潮视频| 亚洲综合二区| 亚洲成A人V欧美综合| 丁香亚洲综合五月天婷婷| 无码丝袜人妻| 欧美激情视频二区三区| av一区二区三区高清久久| 久久人搡人人玩人妻精品| 国产一在线观看| а∨天堂一区中文字幕| 啪啪啪亚洲无码| 日韩亚洲综合在线| 久久婷婷综合色一区二区| 国产午夜福利在线小视频| 乱人伦99久久| 欧美视频二区| 国产精品免费露脸视频| 午夜精品国产自在| 暴力调教一区二区三区| 五月婷婷导航| 国产成人精品视频一区二区电影| 国内精品自在欧美一区| 国产污视频在线观看| 亚洲黄色激情网站| 日韩成人免费网站| 亚洲第一福利视频导航| 国产99精品视频| 激情无码字幕综合| 99久久人妻精品免费二区| 日韩二区三区| 亚洲欧美日韩天堂| 国产乱人伦偷精品视频AAA| 国产麻豆永久视频| 亚洲看片网| 精品久久久久久久久久久| 国产在线拍偷自揄观看视频网站| 亚洲精品成人7777在线观看| 亚洲性一区| 国产精品漂亮美女在线观看| 国产高清自拍视频| 国产小视频免费观看| 伊大人香蕉久久网欧美| 久久精品国产91久久综合麻豆自制| 国产高清在线观看91精品| 91无码网站|