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

面向多最優解組合優化問題的決策求解算法

2022-06-08 09:19:56胡振震袁唯淋羅俊仁鄒明我
國防科技大學學報 2022年3期
關鍵詞:規劃

胡振震,袁唯淋,羅俊仁,鄒明我,陳 璟

(國防科技大學 智能科學學院, 湖南 長沙 410073)

最優化是應用很廣的一個數學分支,很多問題都可以看成是某種形式的最優化問題。最優化問題根據變量是否連續分為兩類,其中離散變量的優化問題稱為組合優化問題。組合優化通常是在一定的約束條件下,根據目標函數要求,從一個有限的集合里尋找一個最優的對象,如一個數、一個集合,或者一個排列[1-3]。在一定約束條件下求目標函數的極值問題也稱為數學規劃,因此組合優化問題也是規劃問題。組合優化的很多理論也基于規劃,如線性規劃、動態規劃、整數規劃等[4-6]。

組合優化的概念源于計算科學,與運籌、決策、管理、經濟等學科緊密相關,科學研究和工程應用十分廣泛。當前組合優化已經發展成一個涵蓋許多規劃問題的學科,比如:旅行商問題、背包問題、劃分問題、指派問題、最短路徑問題、網絡最大流問題、最小費用流問題、最小頂點覆蓋問題、最小支配集問題、最小生成樹問題、斯坦納最小樹問題等。盡管組合優化研究已經取得了豐碩的成果,但因為現實問題的復雜性,仍然還有很多問題需要研究:一是針對新應用問題的模型和方法研究,比如線性約束下的平行機排序問題[7]、“新高改”高中“走班”排課問題[8]、連續背包問題[9]、矩形背包問題[10]等;二是針對已有問題的新方法研究,比如尋路問題的規范排序 A 星算法[11]、圖匹配問題的深度強化學習方法[12]、旅行商問題的圖卷積網絡方法[13]、背包問題的強化學習方法[14]、二次背包問題的剪枝算法[15]等。

本文是針對一類新問題開展研究。現實中很多組合優化問題只有一個最優解,只要找到它就可完成求解,但也存在一些特殊的問題具有多個最優解,而且決策者必須獲取全部最優解才能據此適時做出決策,比如財務預算要在預算總值固定情況下得到所有的可能分配方案。這類問題具有兩個顯著特征:一是最優解是多個且需全部求出;二是目標物品或數值的總和是固定的。兩個特征使該類問題有別于單最優解問題,可稱之為多最優解組合優化問題。目前專門針對多最優解組合優化問題的研究少有報道,本文將其作為研究對象并以固定總和實數子集問題和費城中國餐館購買雞翅問題為例開展建模和分析。

1 相關工作

組合優化是離散變量的最優化問題,往往可以寫成線性規劃的形式。一個求最小目標的問題可用如下數學模型描述:

(1)

其中,x為決策變量,D表示問題或者決策變量的定義域(有限個值組成的集合),F為目標函數,G(x)≤0為約束條件。當求解的變量是整數時,轉變為整數規劃問題。

線性規劃問題常用的求解算法主要是單純形法及其相關擴展算法。整數規劃中的兩種主要算法——分枝限界法和割平面法與線性規劃的算法具有密切聯系[5]。對于變量取值是0或1的整數規劃問題(稱為0-1整數規劃問題),有兩種更為直觀的求解方法:一是枚舉法,即枚舉所有變量等于0或1的所有組合,然后判斷所有組合是否滿足約束條件并使目標函數最優; 第二種是隱枚舉法,同樣是枚舉所有的組合,但引入一些過濾條件,過濾掉局部不滿足約束或目標函數值非優的組合,從而減少計算量[16]。0-1整數規劃問題,也可以看作是多階段決策問題,從這一角度建模可以運用動態規劃的方法求解。動態規劃不是一種特殊的算法,而是一種考察問題的途徑和思路。針對不同的問題,往往需要建立針對性的模型進行求解[5,17]。

除了確定性算法外,組合優化問題還可采用啟發式算法。無論是人為的構造啟發式,還是基于一些最優化方法(如禁忌搜索、模擬退火、遺傳算法等)設計啟發式,主要目的是解決大規模問題中確定性算法無法在有限時間內求得精確解的問題,即利用啟發式算法在有限時間內獲得可接受的解(如近似解)。啟發式往往依賴于問題的性質,需要針對不同類型的問題修改,所以通常情況下啟發式設計并不容易。近年來深度強化學習被發現具有一定的潛力用來求解組合優化問題,可以通過學習來獲得好的啟發式[18-22]。但這些方法主要是針對單最優問題。從需要獲得多最優解的角度看,上述方法中枚舉法、隱枚舉法、基于動態規劃原理的算法最具有獲得問題的全部解的潛力,所以本文的研究主要集中于此。

2 問題定義與求解

2.1 問題形式化描述

多最優解組合優化問題可以定義為:從一個有限集合中尋找多個對象、數、集合或排列,使目標達到最優。它要求將多個最優解全部求出,且約束為等式。在數學形式上可以表示為:

(2)

固定總和實數子集問題和費城中國餐館購買雞翅問題是兩個典型示例。固定總和實數子集問題源于實際預算問題:數據集X包含一系列數據{8.05, 6.98, 6.19, 5, 22.96, 4.71, 4.74, 4.25, 6.34, 2.77, 7.31, 3.59, 18.28, 19.55},希望從中找到多組數(多個子集)使其累加總數為84.03。如果把問題的求解變量看作是取值為{0,1}的變量,即把問題看作是在數據集X中對每個數據做挑選的決策,選中則變量取為1,不選中則取為0,如式(3)所示:

(3)

購買雞翅問題源于一家位于美國費城的中餐館,其雞翅菜單與眾不同,不標注雞翅單價而是標注不同數量雞翅組合的價格,如表1所示。問題是求解購買指定數量雞翅時的全部最優惠購買方案。類似地,可以把問題建模為0-1整數規劃問題(0-1決策問題),如式(4)所示:

Find allx=arg(min(vx))

(4)

其中,b為要購買的雞翅總數。

表1 雞翅菜單

嚴格地說,固定總和實數子集問題是一個多最優解的約束滿足問題,但也可以轉換為最優化問題;購買雞翅問題則是一個標準的優化問題,要求購買固定數量雞翅的同時使得花費的代價最小。

2.2 枚舉和隱枚舉法

枚舉和隱枚舉法是多最優解組合優化問題的直觀解法。以固定總和實數子集問題為例,引入一對相反的不等式,式(3)的約束問題就能轉換為一個優化問題,如式(5)所示。這是典型的整數規劃形式,可以理解為n=14個變量(階段)的0-1整數規劃(決策)問題。

(5)

首先考慮枚舉法,由于數據集中總共包含14個數據,對應于14個變量,因此任意數量(1~14)的數據都可以進行組合,枚舉組合總數為214-1。枚舉過程可以利用字典序循環實現。由于目標函數是所選數據的和,所以每個組合都需要累加計算,運算量接近14乘以枚舉組合數。因此,對于數據集規模為n的問題,枚舉方法的時間復雜度約為O(n2n);枚舉過程不額外增加存儲空間,所以空間復雜度為O(n)。顯然時間復雜度的指數級增長源自枚舉組合數隨數據集規模增大產生的指數級增長。

隱枚舉法思想有些類似于分支定界法,基本思路是:當檢測到某些分支組合不符合約束要求,則避免進一步往下分支,從而減少總的計算量。隱枚舉過程可以理解為搜索樹的擴展,數據集中的每一個數作為一層的決策,以選擇和不選擇兩種情況擴展分支,因此搜索樹是一個不斷擴展的二叉樹,如圖1所示。圖中“×”表示當前節點已經不滿足約束條件,因此可以剪掉。如果不考慮剪枝,那么隱枚舉法的決策樹從根節點到第14層的葉子節點是完整的,葉子節點數量為214,因此時間復雜度為O(n2n),與枚舉法相同。若考慮剪枝則可以使其計算復雜度與枚舉法相比有明顯的下降,但下降的程度與問題的求解特征和方式相關,也與剪枝約束條件相關,不是一個固定值。

2.3 基于動態規劃的方法

若將問題看作多階段決策問題,可以自然地運用動態規劃原理來建模求解。對于一個n階段的決策問題sk+1=f(sk,xk),k=1,…,n,其中sk是第k階段的系統狀態,xk是第k階段的決策變量,優化問題可以看作求決策向量x1,x2,…,xn使得目標函數(指標)達到極值:

(6)

式中,Vk為各階段的目標函數。 根據動態規劃的最優性定理,可以將求指標極值問題轉換為分階段的遞推過程:

(7)

或者:

(8)

式(7)和式(8)分別稱為逆序和順序的動態規劃方程[4]。由此多階段的決策問題可以分解為在各個階段的狀態和目標函數遞推過程中做出最優的決策。

為解決固定總和實數子集問題,首先用一個簡單的示例來考察動態規劃的求解結構。問題是從數據集[1,2,3,4]中選擇數據子集使其和為7。這是一個完全類似的問題,可表示為整數規劃的形式,并轉換為多階段決策問題:

(9)

先考慮順序遞推的方法,定義狀態為決策(選擇該階段決策變量數據)前的數據總和。實際狀態的轉移和指標的遞推為:

(10)

根據各階段的狀態轉移式(10)可知,求決策變量需要從最后一個階段開始,若f4的所有可能值已知,則使f5最優可以將x4求出,進而知道f4的最優值,進一步將x3求出來,不斷遞推直到求出所有階段的決策變量。對于這樣一個求解過程,可以考慮兩種結構:一種是循環,一種是遞歸。循環的結構是先將k-1階段的所有狀態s和指標f求出,再求k階段的狀態和指標,直到最后一個階段,然后再求解決策變量。通俗講就是向前建表(遞推狀態和指標),向后查表(求決策變量)。而遞歸的結構是以遞歸的方式從最后一個階段的狀態開始遞歸,當存在前一階段的未知狀態和指標時,則遞歸地進入前一個階段的計算中,直至第一個階段的狀態和指標計算完成,這樣等價于完成了建表的過程,接著通過查表可以求解決策變量。

圖2給出了順序遞推形式下利用循環方法構造的求解結構。圖中的曲線表示遞推的過程。從狀態s1和指標f1開始遞推到s5和f5,根據式(10)可知f5=7為最優的目標。根據遞推的過程可以計算得到x4=1,s4=3,根據s4=3狀態可以知道有兩條遞推路徑到它。如果是只有一個解,那么查詢過程只要找到一條解路徑即可,但對于多最優解問題,必須獲取全部路徑,可以發現x4=1,x3=1,x2=0,x1=0和x4=1,x3=0,x2=1,x1=1是兩個最優解,圖中兩條藍色的路徑就是最優解的路徑。

圖3給出了順序遞推形式下利用遞歸方法構造的求解結構。圖中的曲線表示遞歸調用時的過程。根據約束條件狀態s5=7已知,但f5是未知的,要求f5則需遞歸的求f4,接著f3,直到f1。因為狀態s定義為決策前的和,因此f1只有一個值為0,其他狀態下的f值都定義為負無窮。那么遞歸完畢同樣得到所需的f值后,也可以根據查表得到解。

圖2 順序遞推循環求解結構Fig.2 Loop solving structure for sequential recursion

圖3 順序遞推遞歸求解結構Fig.3 Recursive solving structure for sequential recursion

從上述分析可知,順序遞推形式下,通過循環或者遞歸都能得到解。類似地,逆序遞推形式同樣也可以使用循環和遞歸兩種求解結構,因此一個動態規劃問題求解在實現上可以有四種結構,而且可以用圖證明逆序循環的結構和順序遞歸的結構是一致的,而逆序遞歸和順序循環的結構是一致的。總的來說,順序和逆序形式與指標推導式相關,順序遞推是指標從前往后推,決策變量從后往前計算;逆序遞推是指標從后往前推,決策變量從前往后計算。循環和遞歸則與指標遞推式的具體實現相關,循環是先計算指標推導式的等號右側項,而遞歸是先調用指標推導式的等號左側項。另外狀態的定義也不是固定的,也可以定義為各階段決策完成后的狀態,對應的則需修改遞推式的下標形式。

2.4 基于0-1決策動態規劃多最優解算法

根據前述分析,在動態規劃遞推建表過程中,多個最優解能夠用遞推路徑表示,與單最優解問題的差異主要體現在查表求決策變量過程中。一般的簡單遞推式查表,只能得到一條最優解路徑,但是多最優解問題必須得到多條最優解路徑。

以順序遞推循環求解結構為例,單個最優解問題可以用一個很簡單的邏輯進行求解:從最后一個階段開始,最優的fn已知,如果fn!=fn-1,那么必有xn=1,如果fn==fn-1,那么必有xn=0,然后根據最優的fn-1在前一個階段繼續求解。但若要得到多條解路徑則不同。

一種直觀解決思路是考慮建表過程中,當遞推某階段出現多種不同決策可以得到相同指標時,記錄多個決策,并在查表時根據這一信息向后搜索。然而這種方式會帶來內存空間的增大,因此若不希望增大內存空間,那么只能在建表的信息里實現對路徑搜索。觀察圖2中第3和第4階段,目標狀態f4=3可以從f3=0和f3=3得到,這是兩條路徑。不同于單最優解的邏輯,在查表過程中已知最優fn情況下,應判斷fn==fn-1+cn-1xn和fn==fn-1兩種情況,當都滿足條件時則采用分支的形式繼續下一階段。由于已經在建表過程中確認解的存在,因此完全可以采用深度優先搜索的方式得到不同路徑的解,而且搜索分支不用深入第1階段,只要達到fn==0(狀態也等于0)就可以終止分支,如算法1所示。

算法1 整數狀態的多最優解動態規劃算法

分析算法的結構可知,在狀態是整數的情況下,利用循環結構來實現狀態遞推,可以方便地利用數組(或列表)數據結構來記錄指標(目標函數)信息,并可同時利用數組的索引信息來表示階段和狀態。但是對于如式(3)這樣的實數問題,階段是整數可方便表示,但狀態卻是實數的,不便于表示。因此基于動態規劃的基本原理,考慮順序遞推的循環結構,針對性地提出一種狀態實數表示的算法,即采用列表加字典的數據結構算法,核心是利用實數表示的狀態作為字典的key值,如算法2所示。

算法2 實數狀態的多最優解動態規劃算法

實數狀態的多最優解算法與整數狀態的算法基本結構一致,差異主要由使用實數作為字典的key而產生。要注意由于實數在計算機中只是高精度的近似表示,所以在作為key以及運算時可能產生的問題可以使用固定模式的實數表達來解決。當然如果從另外一個角度看,實數狀態的問題也可以直接轉換為整數狀態的問題來解決,比如:式(5)可以將c和b同乘以100,那么整個問題變為整數問題,而且解是相同的。分析算法1可以知道,基于0-1決策的動態規劃法的空間復雜度主要由指標f的矩陣決定,為O(nb),b為整數表示的狀態總數(對于上述實數問題,b=84.03×100=8 403)。對于時間復雜度,計算主要體現在建表和查表過程中,建表過程中計算的區域約為nb/2,加上一些比較運算,可以認為在O(nb)量級。而查表過程由解的數量m和表的階段數n決定,最小的情況是nm,而最極端的情況是mn,因此時間復雜度在O(nb+nm)和O(nb+mn)之間的范圍內。

比較算法2與隱枚舉算法,建表過程中的s≤b判斷可以看作隱枚舉法的剪枝約束,而動態規劃相比隱枚舉法的主要優勢在于建表過程中對狀態的規約,相同的狀態合并到一起考慮,而隱枚舉法一直都是以叉樹的形式在擴展。基于動態規劃的算法正是在對狀態的不斷規約中使得其分支數不斷地縮減從而減小計算量。

3 兩種改進的多最優解算法

前面提出的基于0-1決策的算法有效解決了本文多最優解問題的兩個特征帶來的問題。固定物品(數值)總和的問題通過固定目標終狀態解決,多最優解的問題通過考慮多條解路徑查找來解決。同時由于狀態的壓縮使其相比枚舉和隱枚舉算法在計算復雜度上具有明顯的優勢。

購買雞翅問題,要求購買目標數量的物品,求代價最小的多個最優購買方案,仍采用前述方法開展分析。最直觀的求解思路仍然是采用枚舉,枚舉所有的購買方案組合,從中得到全部的最優惠方案。枚舉時需要考慮不同數量組合的雞翅可以購買多次的問題,比如要購買12只雞翅,可以購買3次4只,或2次6只,或1次4只1次8只等。因此購買b(>200)只雞翅,需要在b/4+b/5+…+b/200=n個數據的組合中尋找,比如購買256只雞翅,需要在587個數據中尋找最優組合,那么枚舉總數為2587-1,顯然這種規模已無法有效求解。可考慮用隱枚舉法和動態規劃法,即將問題轉換為n個階段的0-1決策問題。然而,隱枚舉法在購買數量上升到一定的程度后計算時間也變得不可接受,圖4表明采用隱枚舉法在購買數量為28時計算時間已經呈現指數級上升。

圖4 最優解數量較少情況下的計算時間Fig.4 Computation time of cases with small number of optimal solutions

基于0-1決策的多最優解動態規劃算法結果如圖4中紅色虛線所示,圖4中藍色實線是根據量級O(nb)估計的計算時間。結果表明在多數情況下計算時間會接近正比于O(nb),但在某些特殊情況下,比如圖中購買120和216只雞翅的情況,計算時間有跳躍式的上升,而這正是最優解數量較多的情況。由于圖4只是為說明問題,所以只選取了有限情況(曲線上星號標記)進行計算。

3.1 基于相同決策路徑合并的改進算法

顯然前述算法在最優解數量較多的特殊情況下,時間復雜度可能會往差的mn方向發展。這是因為當解的數量很多且階段數量很大時,使用遞歸的查表可能會近似于搜索一個n層m叉樹問題,這顯然是非常耗時的,例如:購買188只雞翅時,最優解的數量為69,利用前述算法計算時間為67.1 s;購買192只雞翅時,最優解數量為300,計算時間已達2 874.7 s,將變得不可接受。

顯然對于要求出全部最優解的給定問題,解的數量m是固定的,那么只能從問題的決策階段上來考慮降低極端條件下的計算復雜性。在前面的問題建模中,采用了0-1決策思路,其中一些階段的決策是購買相同數量的雞翅組合,因此這些階段的作用是相同的。如果能夠對其進行處理和簡化,即將相同作用的決策路徑合并,那么搜索空間將會大幅度減小。基于這種相同決策路徑合并的改進思路,可采用寬度優先搜索實現,如算法3所示。

算法3核心改進在于寬度優先搜索節點擴展過程中對相同決策效果的節點進行了刪減。購買相同雞翅數量情況下的改進算法與原算法計算時間比較如圖5所示,可以看到改進算法相比原來算法具有明顯的改進,在最優解數量比較多的情況下也沒有出現與原算法類似的跳躍性增長。

圖6給出了購買4~300只雞翅問題的改進算法計算時間結果。很明顯看到,計算時間與最優解數量的關系,計算時間總體符合O(nb+nm)的正比估計(圖中藍色線為根據該量級估計的計算時間),但隨著解數量的增加有一定的偏離,且解的數量越多則偏離越大,但總體性能仍可接受。

算法3 相同決策路徑合并的改進算法

圖5 改進算法與原始算法的計算時間比較Fig.5 Computation time comparison between improved algorithm and original algorithm

購買188只雞翅時,改進算法的時間為2.3 s,是原算法的1/29;購買192只雞翅時改進算法計算時間為5.3 s,是原算法的1/542。

圖6 改進算法計算時間與解數量的關系Fig.6 Relation of computation time to solution number for improved algorithm

3.2 基于0-x決策的改進算法

對原算法的改進除了進行相同決策路徑合并的思路外,其實也可以從建模的角度進行改進。如果能夠將0-1決策問題轉變為0-x決策問題,那么決策階段數也可大幅下降,從而降低極端情況下的時間復雜度。因此考慮將問題建模為[0,b/cj]決策問題,其中每個階段求解變量的最大值為b/cj,數學模型如式(11)所示。

(11)

這種改進思路是利用0-x決策建模減少總的決策階段數,使得表的整體搜索深度明顯下降,可采用深度優先類搜索算法實現,計算時間也將有明顯的改善。基于0-x決策的改進方法實現如算法4所示。

算法4 基于0-x決策的改進算法

與前述改進算法最大差別是,前述算法采用0-1決策,而算法4采用[0,b/cj]決策,在建表過程中同一階段的遞推需要比較b/cj種選擇,在查表過程中的遞歸則由原來的2個分支變為b/cj個分支,原理是通過階段內更多的比較和計算來換取決策階段數的下降。

算法4的復雜度仍然由建表和查表過程決定,表的結構為n行b列,只是n由原來0-1決策的b/4+b/5+…+b/200變為count(4,…,cj,…,200)。建表過程的時間復雜度仍為O(nb),查表過程沒有重復分支所以時間復雜度為O(nm)。圖7給出了0-x決策查表和0-1決策查表過程的局部差異,其中j階段的0-x決策對應了i到i-2階段的0-1決策。對應于xj=1,0-1決策有3種情況xi=1或者xi-1=1或者xi-2=1,而這三種情況其實是一個相同的決策,只是由于0-1決策結構導致分成3個階段來決策。所以,0-x決策的建模避免了0-1決策中的很多重復搜索,前面基于相同決策路徑合并的改進算法實質也是解決這一重復問題。

(a) 0-1決策(a) 0-1 Decision

(b) 0-x決策(b) 0-x Decision圖7 0-1決策和0-x決策的查表路徑差異Fig.7 Searching path difference between 0-1 and 0-x decision

圖8給出了購買4~300只雞翅利用基于0-x決策改進算法的運行時間情況,顯然計算時間基本符合時間復雜度O(nb+nm)的正比估計。與圖6相比可知,基于0-x決策的改進算法的性能相比基于0-1決策相同決策路徑合并的改進算法也有明顯的提升,例如最優解數量最多的購買298只雞翅的情況,計算時間也僅為0.437 5 s。由此得出結論,對于存在多次相同決策的問題,更適合建模為0-x決策問題,因為即便采用相同決策路徑合并,基于0-1決策的算法性能也達不到基于0-x決策算法的水平,而每次決策都不相同的問題,則更適合建模成0-1決策問題。(本文所有實驗代碼見:https://github.com/hushidong/multi-optimal-solution-combinatorial-optimization-problem)

圖8 基于0-x決策改進算法計算時間與解數量的關系Fig.8 Relation of computation time to solution number for improved algorithm based on 0-x decision

4 結論

本文提出了一類具有固定物品(數值)總和及多最優解特征的組合優化問題。該問題不同于一般的單最優解組合優化問題,必須求解所有最優解。以固定總和實數子集問題和費城中國餐館購買雞翅問題為例,比較分析了枚舉、隱枚舉和動態規劃三類不同方法。通過動態規劃求解結構分析,明確了狀態定義、指標遞推、決策變量求解對于算法實現的影響。提出基于0-1決策的多最優解的動態規劃算法,并針對實數狀態問題提出了基于字典數據結構的實數狀態表示求解算法。該算法在最優解數量較少時能夠獲得較好的性能,但最優解數量較多的情況下,計算時間呈現跳躍式上升。基于降低時間復雜度考慮,提出了壓縮決策階段的改進思路,并實現了兩種改進算法:基于相同決策路徑合并的0-1決策求解算法和基于0-x決策的求解算法。結果表明通過減小決策的階段數可使算法性能得到明顯提升,兩種改進算法的對比表明具有多次相同決策的問題更適合于建模成0-x決策問題。本文將一類多最優解的特殊的組合優化問題作為一類獨立問題專門進行研究,能夠幫助解決現實中一些多最優解的決策問題,且對開發新的建模和求解方法也有促進作用,下一步將持續推進算法優化和實際應用。

猜你喜歡
規劃
我們的規劃與設計,正從新出發!
房地產導刊(2021年6期)2021-07-22 09:12:46
“十四五”規劃開門紅
“十四五”規劃建議解讀
發揮人大在五年規劃編制中的積極作用
規劃計劃
規劃引領把握未來
快遞業十三五規劃發布
商周刊(2017年5期)2017-08-22 03:35:26
基于蟻群算法的3D打印批次規劃
多管齊下落實規劃
中國衛生(2016年2期)2016-11-12 13:22:16
十三五規劃
華東科技(2016年10期)2016-11-11 06:17:41
主站蜘蛛池模板: 国产农村妇女精品一二区| 中文字幕人成乱码熟女免费| 97综合久久| 国产真实自在自线免费精品| 波多野结衣无码中文字幕在线观看一区二区| 九色国产在线| 干中文字幕| 国产欧美日韩va另类在线播放| 亚洲成人网在线观看| 2019国产在线| 国产伦精品一区二区三区视频优播| 日韩最新中文字幕| 日本国产精品一区久久久| 91香蕉视频下载网站| 国产精品视频白浆免费视频| 日韩成人午夜| 美女被操91视频| 97精品久久久大香线焦| 国产成人久久777777| 国产女人喷水视频| 久久一本精品久久久ー99| 911亚洲精品| 91蝌蚪视频在线观看| 国产第二十一页| 国产香蕉在线视频| 99热这里只有免费国产精品| 国产日韩精品欧美一区灰| 亚洲色图综合在线| 久久中文字幕不卡一二区| 久久精品无码一区二区国产区| 欧美另类图片视频无弹跳第一页| 天天综合网色| 国产精品无码制服丝袜| 成人国产小视频| 国产福利影院在线观看| 欧美激情首页| 亚洲清纯自偷自拍另类专区| 国产69精品久久久久孕妇大杂乱| 丰满人妻久久中文字幕| 亚洲中字无码AV电影在线观看| 在线无码av一区二区三区| 国产无码精品在线| 亚洲欧美日韩另类| 午夜国产理论| 色综合久久久久8天国| 久久久久人妻一区精品| 日韩无码黄色网站| 国产欧美日韩综合一区在线播放| 欧美日韩精品综合在线一区| 在线视频亚洲色图| 国产人前露出系列视频| 四虎亚洲国产成人久久精品| 欧美亚洲欧美区| 热99精品视频| 色老头综合网| 国产高清无码第一十页在线观看| 91免费片| 亚洲国产亚洲综合在线尤物| 女人一级毛片| 尤物成AV人片在线观看| 国产精品美女自慰喷水| 成人永久免费A∨一级在线播放| 精品无码国产一区二区三区AV| 国产视频资源在线观看| 久久无码av三级| 亚洲愉拍一区二区精品| 99热国产这里只有精品9九| 国产精品人成在线播放| 一区二区三区四区日韩| 国产黄色爱视频| 五月六月伊人狠狠丁香网| 欧美激情成人网| 真人免费一级毛片一区二区| 色综合成人| 亚洲成a人片| 日韩天堂网| 全午夜免费一级毛片| www成人国产在线观看网站| 欧美一区二区啪啪| 五月婷婷综合网| 日韩欧美国产另类| 欧美第二区|