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

分批補料反應過程的非固定終端經濟優化控制

2021-08-31 06:59:02湯舒淇薄翠梅俞輝李俊張登峰張泉靈金曉明
化工學報 2021年8期
關鍵詞:優化生產

湯舒淇,薄翠梅,俞輝,李俊,張登峰,張泉靈,金曉明

(1南京工業大學電氣工程與控制科學學院,江蘇南京 211816;2南京工業大學智能制造研究院,江蘇南京 210009;3浙江大學控制科學與工程學院,浙江杭州 310027)

引 言

隨著現代社會對多品種、高質量、精細化大宗化工產品的迫切需求,當前我國石化行業的生產正向著原料多元化、產品高值化發展。由于間歇、半間歇的生產模式更適合生產小批量、多品種、高附加值的產品,因此更受青睞[1-2]。間歇、半間歇生產與連續生產的區別主要在于它的生產操作范圍不局限于傳統意義上的穩態點,同時因為間歇生產的時間有限性,必須對整個生產過程進行有效的排產,這給生產帶來了靈活性,但也給過程優化控制帶來了困難[3]。

目前大部分優化與控制理論方法都是以連續過程為研究對象,通過將批次反應的終端狀態固定,終端狀態作為優化的終端約束,可將用于連續過程的優化控制方法應用到批次生產中,從而解決批次過程經濟優化問題[4]。批次過程優化問題的求解通常采取直接法或間接法。間接法指利用變分法等最優化原理得到問題的數學解析解。例如采用基于龐特里亞金最大原理的優化方法,以求取在批次生產周期終點固定時間下產物濃度最大化的最佳操作曲線軌跡[5]。直接法是將狀態變量或輸入變量離散化,把原動態優化問題轉化為非線性規劃問題(nonlinear programming,NLP)。例如因非線性約束導致的解析解計算困難,可對控制輸入變量進行參數化以找到最佳的批次過程操作條件[6],或用有限元正交配置法徹底離散化優化問題[7]。而離散化后的NLP則可以用各類基于梯度或隨機搜索算法計算求解,如文獻中將改進的蟻群算法用于發酵過程中基質流加率控制軌跡的動態優化[8]。

由于不確定的干擾會導致反應持續時間不確定,因此固定的批次終端優化方法在很多場合不再適用。在不確定干擾下固定終端約束對于批次補料過程是不經濟的[9]。這是因為預先設定生產目標的條件并不適用于外部條件改變后的生產。例如,在不確定的系統參數下,將分批發酵過程的自由終端時間和反應物初始濃度作為優化調節變量,根據發酵的實時狀態變化隨時調控操作參數,獲取最佳的發酵品質[10]。對于多批次迭代學習進行計算的情況,不固定終端將帶來批次不等長的情況,通過搜索相似的模式分別構造它們對應的置信區域[11-12]。抑或通過混合時間標度策略可將自由終端時間的優化問題轉化為固定終端時間的優化問題[13],為解決不等長批次的批間學習提供了較好的參考。

在沒有終端約束的情況下,對于間歇生產的操縱優化而言,操縱變量可行域通常會增加。對于計算量過大且無法通過相關約束縮小可行域簡化計算的大型化工生產過程來說,通常采用貫序優化或雙層優化的結構[14]。雙層優化是指上層通過經濟優化計算操縱變量的最佳操作軌跡[15];底層以上層的優化結果作為固定參數,通過基于優化的最佳操作區域內尋底層局部最優控制,例如用模型預測控制(MPC)實現快速跟蹤控制[16]。該結構成功地將經濟優化問題和控制問題分離開,以降低解決問題的復雜性,已廣泛應用于復雜化工過程優化系統。然而雙層優化結構的潛在缺點是上下層所采用的不同模型會導致優化效果偏離預期[17]。因此,近年提出將操作優化和最優控制集成到同一層的單層優化控制方法。為了提高單層優化控制的魯棒性,設計了魯棒的狀態估計和退避約束的方法[18];另外一些研究者采用了一種復合優化目標函數,同時考慮了過程產量和目標函數敏感性[19]。文獻[20]還提出了一種基于NMPC和DRTO的批次協同優化控制方法,使NMPC更加適用于批次生產過程,可以同時考慮批次生產長度對經濟的影響。在單層的優化控制中,應用了滾動時域控制(receding horizon control,RHC)的單層優化控制又可以視為一種經濟模型預測控制(economic model predictive control,EMPC),應用于反應過程中的研究[21-24]。例如文獻[25]通過為同一過程制定不同的EMPC經濟目標函數,以實現不同的經濟優化效果。也有文獻研究了非固定終端的收斂性問題[26-27],但利用EMPC目標函數的多樣性,將非固定終端約束應用于批次經濟優化控制研究成果并不夠充分。

因此在上述研究的基礎上,提出了一種應用于分批補料的非固定終端經濟優化控制方法,用于考慮批次生產長度對經濟的影響,解決批次周期不確定約束優化控制問題。本文基于經濟模型預測控制方法,提出了非固定終端的經濟優化問題;提出控制變量差異參數化的方法提高優化計算效率;采用內點罰函數法求解帶非線性等式、不等式約束的優化問題;給出了基于滾動時域的優化控制方法及其流程圖;在苯胺加氫分批補料工藝中進行了實驗測試,并與多回路復雜PI控制進行了對比研究,以驗證該方法的有效性。

1 非固定終端優化控制

1.1 優化問題的提出

任意被控對象動態過程,反應過程可以寫成如下常微分方程組(ordinary differential equations,ODEs)形式:

令優化求解的初值x(ts)=xs,u(ts)=us,ts為優化時的當前時刻,定義一個包含(xs,us)的不變集Ω,且Ω?X×U,則x(t),u(t)∈Ω代表在約束內u(t)的作用下x(t)也處于約束內。

在連續過程中,EMPC通過在滾動時域算法中求解含經濟目標函數的最優化問題,實時更新系統的運行條件。該優化問題形式通常如下:

但將EMPC應用到間歇生產過程時,則需要考慮如下問題:(1)間歇生產無穩態點,也因此不存在預測范圍內需要達到的設定點;(2)間歇生產的生產時間有限,如果預測范圍PEMPC保持不變,隨著時間的滾動,預測范圍將超過生產時間范圍。

因此考慮將針對連續過程提出的EMPC應用于間歇生產過程,并且對間歇生產模式提出有針對性的改進。EMPC策略在連續過程中的收斂性通常通過增加預測時域范圍來保證。根據文獻[26],隨著優化范圍的增長,系統性能會收斂到最佳狀態,相較于無限時域優化有指數級誤差衰減。因此優化時域應當在可計算范圍內足夠長。對于批次反應過程來說,最長的有效預測時域即批次的終端,則優化預測時域P=Tf-ts,以獲得足夠的收斂性能。假定間歇生產的一個批次周期為Tf,用批次結束時間的產量約束x(Tf)=xf代替預測范圍內的設定點,以達到預先設定的生產目標。因此在每個采樣點要解決一個固定終端的最優控制問題,如式(5)所示:

終端時間放寬后,對于經濟目標函數的定義也能不局限于最小化控制輸入。原本的生產計劃指標是通過一個終端等式約束達到的,給求解計算帶來了很大的復雜度。由于實際生產中的經濟目標通常是利益最大化。

通過將目標函數以最大化利益的方式設置,即可將原本的生產計劃指標的等式約束轉化為產品收入項;將原本的最小化控制輸入目標函數,轉化為生產成本項。最終達到最大化利益生產的目標。而根據間歇過程的生產特點,產品收益通常在一個批次生產完成之后才會產生,即在目標函數中只有終端項。而生產成本則是在生產前有一次性成本,在整個生產過程中有連續累積的成本,因此既有終端項又有積分項。完整的經濟目標函數形式應當如下:

其中,Veconomic、V1、V2、V3分別為過程經濟利益函數、累積成本函數、產品收益函數、生產前一次性成本去量綱后的函數,V3通常是一個常數。經過改進后非固定終端的優化控制問題如式(7)所示。

在t時刻的最優控制輸入序列集為Φ,它是一系列控制變量的集合。這是因為當終端不固定的情況下,不同的操縱變量有不同的離散化規則,具體在1.2節中介紹。

可以觀察到,上述表達式計算的是一個動態經濟優化問題,且通過計算動態經濟優化和在每個采樣時刻更新系統狀態,將開環優化轉化為反饋控制策略。當不考慮過程干擾以及模型失配等不確定性因素,該優化問題在最初一次求解得到的最優終端時間與后續過程中計算得到的最優終端時間應當是一致的。而發生擾動后,經濟最優的情況會使Tf趨向更有利的區間。而反應最終的產物x(Tf)將由優化問題根據過程情況與經濟條件調整。

1.2 控制軌跡與狀態變量平穩性

式(7)的動態優化問題通過優化計算的目標函數來優化控制輸入變量,可以最大程度地降低過程成本。但是,由于該優化中存在較大的可行域區間,所得的控制輸入的最佳軌跡并不一定適合化工生產的實際過程。操縱變量的反復大跳動與引起的狀態變量急劇增減可能會在運行過程中損壞設備,或使具有高敏感度的系統進入失控狀態。為限制狀態變量或操縱變量的跳變,可在約束中加入導數或累積上限。但是如果遇到確實需要較大跳變才能抵消的大擾動,該約束會限制控制動作的行為,以至于控制可行性不高。因此,在原目標函數加入限制重要系統狀態變量的導數和操縱變量的累積差項,這樣當經濟項大到可以忽略控制平穩項時,意味著需更新控制動作。更新的目標函數如下:

其中,α和β分別是維度為nx與nu的加權因子向量,α的第k個和β的第i個分量分別對應系統狀態變量平緩變化的程度和控制輸入變量的穩定變化強度。

2 控制變量的參數化

在解決式(7)的優化問題中,使用控制變量參數化法(control variables parameterization,CVP)描述控制變量的軌跡。在控制變量參數化方法中,只有控制變量被離散化,而其他變量則仍然按照連續的考慮。控制變量的參數化給無限維的優化問題轉化為有限維的NLP提供可能,而保留其他向量的連續特性又大幅度降低了求解規模。

面對批次反應過程的反應周期不定的情況,使用傳統的固定網格化離散方法固然是一種可行的方法,然而該方法下,并不能根據過程動態特性自主選擇不同控制變量的微調精度。另外,若要更精確地逼近最優控制軌跡,就要將網格寬度設置得短一點,也增加了網格數量,且加大了NLP的維數和計算量。因此對于不同控制變量的參數化定義方法非常重要,這里選擇了根據控制變量動態特性合理確定網格,并且隨著迭代實時更新的控制序列分段方法,如圖1所示。

圖1 控制變量的獨立參數化示意圖Fig.1 Schematic diagram of control vectors respectively parameterized

由于批次終端時間并不事先可知,因此當優化所得的批次終端時間變化時,需要在每個優化問題求解之前(即當前時刻ts)時更新分割參數,如式(9)和式(10)所示。

這兩個更新公式分別對應于控制范圍超過預期的終止時間和控制范圍與預期終止時間的比例小于ωi。其中ωi是一個無量綱參數,表示第i個控制變量的控制時域在一個批次周期中合適的占比。

3 優化問題的求解

為解決優化問題,本文將采用內點最優化法(interior point optimizer,IPOPT)進行求解,它通過在目標函數中添加障礙項消去不等式約束,并通過迭代減小障礙權重,不斷逼近最優解。通過控制變量參數化方法將式(8)中的優化問題轉化為NLP后即可以直接求解,但所需求解的問題是一個含有非線性等式約束、不等式約束的較為復雜的約束優化問題。為了降低求解難度,采用了ODE求解集成和障礙函數法分別解決等式、不等式約束,將有約束的NLP轉化為無約束的NLP來求解。

式(7)問題中的等式約束即過程模型的微分方程,將微分方程的求解記作算子D(·),則t>ts階段的狀態變量的預測值,通過ODE求解器所求的解可以表達成如下形式:

其中,邊界條件x(ts)=xs,且ts<σts)時刻狀態變量的數值解x?(t)。在目標函數式(8)中將所有的狀態變量用式(11)替代,即將微分方程的數值解集成到目標函數的計算中,則自動滿足等式約束。則新的優化問題轉化為如式(12)所示的NLP。

考慮到一些優化問題的極小值點可能位于可行域的邊界上,為得到這樣的極小值點,內點罰函數需要在迭代過程中不斷調整罰因子來逐步削弱障礙函數對最優值點的影響,從而使算法產生的點列在可行域內部逐步逼近原規劃問題的極小值點。本方法適合處理求解最優解在可行域邊界上但在靠近可行域邊界時有次最優解的問題。考慮如式(12)所示的優化問題,通過控制變量參數化后優化操縱變量為m=[u'1,…,u'nu,Tf]',不等式約束條件用fi表示,ncon是約束條件的個數,即可行域內點集為:

建立懲罰項

其中,π是懲罰障礙權重,fbarrier為障礙函數,為了獲得連續梯度,選取對數形式:

給定初始π0>0,m0∈X,權重衰減系數γ>1,允許誤差δ>0。求解P(mk,πk)的近似最優解,即|?m P(mk,πk)|≤δ,如果fbarrier(mk)/πk>δ則令πk+1=γπk,否則結束迭代計算,而近似最優解即為mk。

假設存在最優解m*,則同時存在對偶最優變量λ*與原始操縱變量一同滿足Karush-Kuhn-Tucker(KKT)條件:

在這個問題中,λi可相當于

則當π→∞時,互補松弛條件成立,即內點法產生的最優解相當于滿足KKT條件。

4 非固定終端經濟優化控制流程

非固定終端經濟優化控制遵循滾動時域優化控制,直到滿足終止條件為止。如圖2所示,在初始化系統后,根據式(8)的目標函數和式(7)中的等式、不等式約束,形成一個動態優化問題。根據當前的測量信息[x(ts),u(ts),d(ts)],用式(11)將ODE求解器集成到目標函數的計算中,用以根據當前時刻狀態估計未來時刻狀態變量的值。通過CVP離散化操縱變量,將動態優化問題轉化為NLP。通過IPOPT方法求解帶有邊界約束的NLP,從而獲得最優批次生產時間Tf*與最優操縱變量序列uij*,以ΔtOC執行操縱變量,依據式(9)、式(10)更新Nu,并且向前推進ts。通過判斷Nu的個數判斷所剩時間是否不足以進行一次優化控制,如果小于1則不再進行優化求解,再根據時間ts是否達到Tf結束一批次的生產,否則,重新測量當前狀態信息,在下一時刻重新求解優化問題。

圖2 非固定終端經濟優化控制流程Fig.2 Flow chart of economic optimization control with unfixed terminal

根據流程圖,非固定終端經濟優化控制的具體步驟如下:

(1)初始化參數設置與初始操作條件;

(2)從時間ts開始,根據優化問題式(12)求解出最優控制參數集(Φ*)和最佳批次終端時間(Tf*),并計算終端狀態與該時刻下的經濟目標函數;

5 實驗測試

5.1 分批補料反應過程機理模型

本算法使用的測試案例來自文獻[28]的苯胺(aniline,A)加氫制環己胺(cyclohexylamine,CHA)模型,其反應方程式如下:

苯胺加氫通常采取批次生產的模式,在連續攪拌釜式反應器內一次性加入苯胺,然后以批次或半批次的方式加入氫氣,反應是放熱的,需要對反應器降溫以帶走多余的熱量。反應所需的夾套換熱反應器如圖3所示。

圖3 分批進料反應器示意圖Fig.3 Schematic diagramof fed-batch reactor

描述該系統的動態數學模型由總質量平衡、3個成分的平衡、反應器液體上的能量平衡和夾套能量平衡組成。

假設密度恒定,質量平衡如下:

其中,Fin是氫氣進料的體積流量,VR是容器內反應物的體積。

成分平衡如下:

反應器的能量平衡如下:

夾套的能量平衡如下:

反應器壓力由溫度和液體組成計算得出,苯胺和環己胺的蒸氣壓常數以及氫的亨利定律常數是使用Chao-Seader物理性能數據包從Aspen Plus獲得的數據計算得出。

其中,Xa、Xb、Xc分別為苯胺、氫氣、環己胺在反應器中的摩爾分數。

5.2 目標函數的設置

苯胺加氫過程的經濟目標函數設置如下:

其中,上角標“f”代表t=Tf時各種變量的值,從其物理意義中可以看出,經濟目標函數的值即過程加料、冷卻水消耗和一次性投料的成本減去產品收益。然而由于實際上的目標函數還有關于穩定運行的項,因此會將經濟函數的數量級縮放到0~1之間,即在Veconomic外乘系數ε=1/Vmax,Vmax是理論上能達到的最大收益,即不需要任何成本,將所有一次性投料都轉化為產品所能獲得的收益。然而實際上這是不存在的,所以經濟項永遠小于1。

在約束條件的問題上,各個操縱變量、狀態變量的邊界約束不僅包括生產時的物理約束、安全約束,還包括根據實際化工過程經驗所得的可行約束。即根據所需產品的質量,估計推導出操縱變量的操縱范圍,并作為過程約束加在優化問題中。由于無終端約束,收緊了過程約束,因此求解的速度大幅度提升。在無擾動發生時,生產目標的可達性有了新的保障。而大的擾動發生時,生產指標達到的概率急劇降低,此時基于最大化收益的優化方法也能極力止損。

進料流量的變化會引起苯胺加氫反應過程系統狀態發生更大的變化,而系統狀態變量中需要進行安全約束的主要就是氫氣濃度和反應器溫度,反應器壓力和冷卻劑出口溫度也幾乎完全分別取決于這兩個變量,因此為了限制溫度上升的速度需要對溫度經行斜率限制。除此以外,開環下整個批次的時間大約在300 min以內,取控制范圍在預測范圍的1/10左右,為了對過程速率有更加細致的控制,將氫氣進料的控制變量離散化參數設置得比冷卻劑進料參數更細膩。因此針對α、β、Nu、ΔTu四個參數的整定結果見表2。

表2 優化過程的參數整定Table 2 Tuning parameter settings of the optimized control

5.3 仿真測試

基于上述苯胺加氫反應模型,采用C++語言進行批次反應過程控制與優化仿真,實現環己胺制備過程的單批次周期的最優化。為了測試本方法在變批次長度上的優勢,將優化控制與選擇與分程結構的多回路PI控制進行了對比。PI控制的條件設定參照文獻[29]用離線優化得到的最佳參數運行,分別用冷卻液進料與氫氣進料進行獨立回路控制,再為氫氣進料增加低選控制結構,避免進料導致壓力過大,而溫度控制采用兩種輸入的分程控制結構。PI控制的參數調節依據文獻[30-31],為KcTC=0.05,TiTC=10 min,KcPC=0.3,TiPC=30 min。在50 min時對冷卻劑入口溫度加入了擾動,從350 K升高到400 K,分別測試了有無擾動狀況下的兩種方法控制輸出情況對比,見圖4~圖9。

圖4 有擾動下冷卻水流量曲線對比Fig.4 Comparison of coolant flow with disturbance

控制穩定性如圖4和圖5所示。在PI控制下所需的冷卻劑用量明顯多于優化控制下,擾動發生時控制變量也發生了振蕩。相比本文提出的優化控制方法的控制曲線一直處在可接受的變化速率中。從圖6和圖7中在擾動發生時溫度、壓力這兩個重要被控變量曲線來看,PI控制下的壓力曲線在擾動發生時總是會出現振蕩。實際上,在反復調試PID參數后發現擾動帶來的振蕩是與反應機理模型相關的,是單回路閉環控制不可避免的。而優化控制采用一種滾動優化控制的策略,通過求解一連串控制序列和預測狀態變量可以很好地避免振蕩的發生。

圖5 有擾動下進料流量曲線對比Fig.5 Comparison of feed flow with disturbance

圖6 溫度控制抗擾動Fig.6 Anti-disturbance comparison of temperature control

圖7 壓力控制抗擾動Fig.7 Anti-disturbance comparison of pressure control

反應終止時間如圖8和圖9所示。在擾動的影響下批次的反應明顯變緩慢了,這是因為溫度升高后必須降低送料速度,以避免反應器內溫度超標。因此兩種控制方法下批次終止時間都變長了。除此以外,可以觀察到PI控制下的反應在接近反應結束時,反應速率急劇下降。根據式(22)可以推斷出,此刻反應器中的苯胺濃度較低,則反應速度會下降,為了促使反應的繼續,則需要提高另一種反應物的濃度,即需要向反應器中加入大量的氫,促使反應速率保持在較高的水平,這也是反應器中壓力在這段時間驟然升高的原因。在PI選擇控制的策略下,壓力過高時選擇犧牲溫度控制。對于優化控制而言所有變量都在同一個優化中求解,因此壓力過高的情況在預測、優化中被避免了。

圖8 反應器內苯胺濃度變化Fig.8 Changes of aniline concentration

圖9 反應器內環己胺濃度變化Fig.9 Changesof cyclohexylamine concentration

兩種控制下的反應經濟情況如表3所示。其中的凈收入計算方法是:凈收入=反應產物獲益-反應過程消耗(包括一次性投入物料、加料、冷卻等消耗)。可以看到無擾動情況下,PI控制的凈收入甚至比優化控制要高,因為在無擾動的情況下PI控制的反饋控制足夠將反應穩定在較優的操作區間。但是算上時間成本后,PI控制的平均收益則略低于優化控制。在有擾動的情況下,兩種控制下的經濟收益均低于無擾動,且PI控制的凈收入下降了45%,相比之下優化控制凈收入降低了25%的情況較好。

表3 兩種控制下的經濟情況對比Table 3 Economic profitability of two kinds of control

除此以外,對于分批補料過程的特性,當容器內的一次性投料的反應物濃度低時,需要反應速率大幅下降,大量的投入消耗和時間成本對反應物產量的提升不大。因此在終端不固定的情況下可以獲得更多方面的收益。

6 結 論

本文提出了一種分批補料過程的非固定終端的經濟優化控制方法,該方法將非固定終端經濟優化與多變量協同控制集成在滾動時域控制結構中,可以更好地實現單批次內生產利益最大化。由于將批次生產過程的關鍵操縱變量和批生產周期的持續時間都視為優化和控制問題的自由度,因此可以較為準確地管理批次生產的終止時間,在反應條件允許范圍內最大程度壓縮低效率的生產階段。基于目標函數去預測評估優化最佳操縱變量和操作時間,能夠穩定地將過程中發生的不確定擾動引起的偏移盡快反饋調控現有條件下經濟更優的操作區間,不斷更新關鍵操縱變量的控制分段函數的分割數及其寬度,從而可靈活優化操縱變量和操作時間的軌跡。單層的經濟優化控制算法避免了雙層優化控制的滯后性,并且基于懲罰函數內點法的實時優化在有效范圍能最快地求解最優軌跡線,解決單批次周期的優化和控制問題,與復雜多回路PI控制對比分析發現,本方法明顯更靈活有效,測試結果明顯優于基于選擇-分程的復雜控制方法,通過優化計算每個批次生產的最佳操作條件及其生產周期,實現批次反應過程生產時間與經濟效益的最優化管理。

符號說明

——分別為反應器內苯胺、氫氣、環己胺的初始濃度,kmol/m3

——冷卻液的比熱容,kJ/(kmol?K)

DR——反應器底部直徑,m

Fin,0——加料流量初值,m3/s

Fj0——冷卻液流量初值,m3/s

Fjmax,Fjmin——分別為冷卻液流量的最大、最小值,m3/s

Fin,max,Fin,min——分別為加料流量的最大、最小值,m3/s

ΔHR——反應放熱,kJ/kmol

k0——反應速率常數,m3/(s?kmol)

Tin——反應物進料溫度,K

Tjin——冷卻液入口溫度,K

Tjout,0——冷卻液出口初始溫度,K

TRmax,TRmin——分別為反應器溫度的最大、最小值,K

TR0——反應器溫度的初值,K

Vj——冷卻液體積,m3

VRmax,VRmin——分別為反應器內反應物體積的最大、最小值,m3

——反應器內反應物初始體積,m3

θa,θb,θc,θj——分別為苯胺、氫氣、環己胺、冷卻液的經濟參數,CNY/mol

ρj——冷卻液的密度,kg/m3

猜你喜歡
優化生產
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
用舊的生產新的!
“三夏”生產 如火如荼
S-76D在華首架機實現生產交付
中國軍轉民(2017年6期)2018-01-31 02:22:28
基于低碳物流的公路運輸優化
現代企業(2015年2期)2015-02-28 18:45:09
安全生產重于泰山
主站蜘蛛池模板: 91无码人妻精品一区| 天堂岛国av无码免费无禁网站| 国产精彩视频在线观看| 色悠久久综合| 91 九色视频丝袜| 亚洲va在线∨a天堂va欧美va| 99re经典视频在线| 久久精品无码一区二区国产区| 免费日韩在线视频| 日本欧美在线观看| 国产精品久久久久久久久久98| 久草国产在线观看| 精品国产福利在线| 国内精品视频在线| 午夜a级毛片| 欧美五月婷婷| 亚洲成a人在线播放www| 伊人久久大香线蕉aⅴ色| 免费在线看黄网址| 一本大道无码日韩精品影视| 福利一区三区| 国产成人精品免费视频大全五级| 99在线视频网站| julia中文字幕久久亚洲| 欧美色图第一页| 成人福利在线看| 99re在线免费视频| 亚洲男人天堂2018| 91高清在线视频| 欧美福利在线观看| 免费国产无遮挡又黄又爽| 国产一区二区色淫影院| 丁香婷婷激情网| 国产区福利小视频在线观看尤物| 色综合天天综合中文网| 亚洲一区二区三区中文字幕5566| 91精品国产自产91精品资源| 亚洲男人的天堂在线观看| 午夜色综合| 亚洲中文无码av永久伊人| 2024av在线无码中文最新| 人妻丝袜无码视频| 日本影院一区| 亚洲日韩AV无码一区二区三区人| 国产极品粉嫩小泬免费看| 夜夜操天天摸| 久久77777| 日韩亚洲综合在线| 国产拍在线| 久久美女精品| 女人18毛片水真多国产| 国产欧美视频综合二区| 97国产成人无码精品久久久| 久久久久久久久亚洲精品| 国产一区二区网站| 国产视频入口| www.日韩三级| 玖玖精品视频在线观看| 亚洲精品色AV无码看| 国模沟沟一区二区三区 | 91外围女在线观看| 欧美成人午夜视频免看| 亚洲国产精品VA在线看黑人| swag国产精品| 美女扒开下面流白浆在线试听| 亚洲高清在线天堂精品| 国产日韩欧美在线视频免费观看| 色综合久久综合网| 波多野吉衣一区二区三区av| 国产一区二区三区视频| 国产一级毛片高清完整视频版| 日本一区中文字幕最新在线| 精品久久久无码专区中文字幕| 亚洲手机在线| 99re在线观看视频| 国产va欧美va在线观看| 亚洲手机在线| 国产精品亚洲专区一区| 国产亚洲精品无码专| 久久男人资源站| 日韩高清欧美| 国产精品欧美激情|