?

基于COOT 優化算法和改進型重復PI 控制的三相LCL 型光伏并網系統

2024-02-06 11:22楊春輝屈莉莉戴宏躍
關鍵詞:跟隨者領導者阻尼

楊春輝,屈莉莉*,廖 慧,戴宏躍

(1.佛山科學技術學院機電工程與自動化學院,廣東 佛山 528225;2.廣東創電科技有限公司,廣東 佛山 528000;3.廣州科技貿易職業學院 智能制造學院,廣東 廣州 511442)

2020 年9 月,中國在聯合國大會上提出了“雙碳”目標:在2030 年之前實現“碳達峰”以及在2060年之前實現“碳中和”的目標[1]。在“雙碳”目標的引導下,光伏產業并網發電在政策扶持和鼓勵下迅速發展。光伏陣列的最大功率跟蹤(MPPT)和光伏并網逆變器是光伏產業并網發電系統的重要環節,光伏并網系統由光伏陣列和并網逆變器兩個部分組合而成[2]。兩級式光伏并網系統通常利用前級DC/AC 電路完成MPPT,實現光伏陣列最大功率輸出,然后再經過后級DC/AC 電路并入到公共電網[3]。作為太陽能光伏發電系統與電網互連的關鍵接口裝置,MPPT 和光伏并網逆變器控制策略直接決定系統并網運行的可靠性和安全性。

傳統MPPT 方法(如電導增量法INC、擾動觀察法P&O)雖然有很好的局部尋優能力,但是需要按照一定步長進行追蹤,步長過大會出現功率振蕩,步長過小追蹤到最大功率點的時間過長,而且在局部遮擋情況下,傳統MPPT 方法還很容易陷入到局部最優[4]。近年來出現的人工智能算法(如蟻群算法、灰狼算法、白冠雞優化算法等)以其超強的全局搜索能力成為學者們關注的熱點。

并網逆變器的逆變器橋輸出電壓中含有豐富的開關諧波,需要在逆變器中引入濾波器,LCL 型濾波器與L 型濾波器相比具有更好的諧波抑制能力和衰減特性,因此被廣泛應用于并網逆變器[5]??紤]到LCL 濾波器是一個三階系統,存在諧振尖峰,需要合適的阻尼策略來穩定系統。常用阻尼策略可分為:無源阻尼和有源阻尼,無源阻尼采用串聯或并聯阻尼的方式增加系統阻尼,但該方法會帶來額外的功率損耗;有源阻尼從控制策略的角度解決諧振尖峰,通過虛擬阻尼進行狀態反饋,實現諧振尖峰抑制,比如電容電流反饋阻尼控制,但是需要額外增加電流傳感器[6]。對于并網電流的控制,常見的控制策略有PI 控制、PR 控制、重復控制等。文獻[7]指出PI 控制具有結構簡單、算法成熟等優點,但對正弦量給定無法做到無靜差跟蹤;文獻[8]采用電流比例諧振(PR)控制,省去了坐標變換,但由于使用比例諧振控制器的個數較多,不僅增加了計算量,而且系統變得更復雜;文獻[9]指出重復控制有良好的穩態特性,能夠抑制周期性干擾,但在動態響應過程中由于延遲了一個基波周期,當系統給定發生變化時,動態性能較差;文獻[10]提出基于重復控制+PI 控制復合控制,彌補了重復控制動態性能較差的缺點,但建模過程忽略采樣計算延時帶來的影響。

針對傳統MPPT 算法追蹤到最大功率點所需時間長以及容易陷入局部最優的問題,本文提出一種基于COOT 優化算法。通過對白冠雞種群進行層級劃分,適應度值高者為種群領導者,其余為跟隨者。跟隨者具有兩種位置更新方式,分別為主動更新和被動更新。主動更新過程中,跟隨者根據鏈式運動或隨機運動更新位置,并不依賴領導者;在被動更新過程中,跟隨者需要向領導者方向發生靠攏,種群需要朝著最佳區域前進,因此在跟隨者向領導者靠攏的同時,領導者也要同時不斷調整自身位置向最優區域靠近。通過這兩種位置更新方式,可以使得在追蹤功率的時候避免陷入到局部最優。

針對LCL 在高頻處會帶來一個較大的諧振以及考慮到在實際生產過程中,廠家一般在逆變器側加入電流傳感器來實現過流保護,為了節約成本,本文提出僅采用逆變器側電流反饋的有源阻尼控制策略。為了避免求解三階系統特征根帶來的數學限制,通過分析三相LCL 等效電路,不同于用傳統伯德圖進行分析,本文提出基于阻尼因子的設計方法。在數學建模過程中,把PI 控制和采樣計算延時加入到有源阻尼高頻諧振抑制中一起建模,然后在這個基礎上再加入重復控制。為了減少奇次諧波,內模只選取周期采樣次數的1/2,搭建出基于奇次基波級聯式結構的重復PI 控制策略的并網逆變器??紤]到在實際中電網頻率會出現一定波動,利用拉格朗日插值法逼近內模的小數部分,得出基于冪級數展開的分數延遲濾波器,通過該濾波器可以使重復控制的內膜更貼近實際電網基波頻率,從而提高控制器的諧波抑制性能。

1 最大功率點跟蹤

圖1 為COOT 優化算法的MPPT 系統框圖,其由5 塊光伏陣列組件、Boost 升壓電路、COOT 的MPPT 控制算法、PWM等構成。

圖1 COOT 優化算法的MPPT 系統框圖

1.1 COOT 優化算法原理

COOT 優化算法主要模擬白冠雞在自然界中獲取食物的行為,從而實現算法尋優的目的[11]。

1.1.1 初始化

假設白冠雞種群大小為N,從白冠雞種群中選取10%的個體為白冠雞的領導者,記為LeaderPos,剩余的白冠雞為跟隨者,記為FollowPos。對白冠雞種群位置(包括領導者以及跟隨者)進行初始化處理,即

其中,CootPos(i)表示第i 只白冠雞的位置,ub 和lb 分別表示搜索空間的上界和下界,rand 為0 到1 的隨機數。

1.1.2 跟隨者位置主動更新

當一個取值范圍為[0,1]的隨機數rand1>0.5,白冠雞跟隨者將選擇主動更新策略,主動更新策略包括兩種運動方式,通過一個取值范圍為[0,1]的隨機數rand2 來決定采用哪一種運動方式。

(1)隨機運動。當rand2<0.5,跟隨者進行隨機運動,為了體現運動的隨機性,首先在位置空間生成一個隨機位置Q,其數學公式為

隨機運動有利于算法跳出局部最優解,使用這種方法更新白冠雞跟隨者位置的計算表達式為

其中,FollowPos(i)表示第i 個跟隨者的位置,R1是區間[0,1]中的隨機數,參數A 為[0,1]中線性遞減的因子,t 為當前迭代次數,Iteration 為最大迭代次數。

(2)鏈式運動。當rand2>0.5,執行鏈式運動,算法將相鄰兩只白冠雞中間位置實現鏈式運動,其數學公式為

其中,FollowPos(i-1)是上一只白冠雞跟隨者的位置。

1.1.3 跟隨者位置被動更新

當一個取值范圍為[0,1]的隨機數rand1<0.5,白冠雞跟隨者將選擇被動更新策略。通常情況下,白冠雞跟隨者根據白冠雞領導者的位置調整自已的位置并且朝著它們的方向移動,利用k 來控制白冠雞領導者的位置引導,即

其中,i 是跟隨者的索引數,k 是領導者的索引數,mod 是取余函數。跟隨者隨領導者的位置更新公式為

其中,LeaderPos(k)為第k 個領導者的位置,R 和R2為區間[0,1]的隨機數。

1.1.4 領導者帶領跟隨者走向最佳區域

白冠雞種群的領導者不斷更新它們朝著最佳區域目標方向的位置,從而帶領白冠雞跟隨者種群走向最佳區域,白冠雞領導者的位置更新計算表達式為

其中,LeaderPos(i)為第i 個領導者的位置,gBest 為種群最優位置,R3和R4為區間[0,1]的隨機數,參數B 為在[0,2]中線性遞減因子。

1.2 COOT 優化算法的MPPT 流程

(1)實時采集光伏陣列光照強度和環境溫度,然后得到光伏陣列輸出特性曲線(P-U 曲線);

(2)利用COOT 算法進行搜索,追蹤到光伏陣列工作在全局最大功率點處的功率,最終得到最佳電壓;

(3)通過PWM控制Boost 電路,使得光伏陣列工作在全局最大功率點處(Boost 電路通過控制占空比來調節電壓,由于電壓跟功率是有對應關系的光伏陣列曲線,得到最佳占空比即得到最佳電壓即得到最大功率);

(4)當使用COOT 算法追蹤GMPP 時,考慮到外部光照或者溫度會發生變化,光伏陣列輸出特性曲線會發生偏移,GMPP 的位置也可能發生變化,為降低功率損失需要重新啟動COOT 算法進行追蹤。本文設置功率變化情況滿足下式就重啟算法

本文算法的具體步驟為:

(1)初始化COOT 算法的參數,包括:設置白冠雞數量N=50,其中,白冠雞領導者的數量LeaderPos=5,白冠雞跟隨者的數量FollowPos=45,最大迭代次數Iteration=15,以及初始化白冠雞種群的位置;

(2)將采集到電壓、電流值進行采樣;

(3)判斷迭代次數是否滿足迭代跳轉次數15 次,若滿足則轉至步驟(6);

(4)計算一開始設置白冠雞領導者和白冠雞跟隨者對應的適應度值(適應度值即功率),找出白冠雞種群個體的最優位置gBest(位置即占空比);

(5)把第一輪得到的最佳適應度值和最佳位置代入到位置更新中,產生隨機數rand1 和rand2;若rand1>0.5,跟隨者進行主動更新,并不依賴領導者,當rand2>0.5,跟隨者進行鏈式運動,當rand2<0.5,跟隨者進行隨機運動;若rand1<0.5,跟隨者進行被動更新,跟隨者需要向領導者方向發生靠攏。種群需要朝著最佳區域前進,因此在跟隨者向領導者靠攏的同時,領導者也要同時不斷調整自身位置向最優區域靠近;

(6)記錄本次迭代中的最佳適應度值;

(7)計算本次迭代更新的適應度值與目前找到的最佳適應度值是否滿足式(8),若不滿足則轉回步驟(1);

(8)測量到最佳適應度值即為全局最大功率,結束算法。

基于COOT 的MPPT 控制算法的流程圖如圖2 所示。

圖2 基于COOT 的MPPT 控制算法的流程圖

2 LCL 高頻諧振抑制

為了抑制LCL 高頻諧振,在過去10 年里,基于狀態變量的主動阻尼控制策略被廣泛提出[12]。在實際生活中,工廠生產的并網逆變器通常將逆變器側電流傳感器與逆變器產品集成,以實現過流保護為目的。僅采用逆變器側電流反饋不僅可以節約成本,還可以避免三相LCL 在DQ 變換解耦時需要進行復雜的數學運算,同時為了避免求解三階系統特征根時的數學限制,提出基于阻尼因子的設計方法。在數學建模過程中加入采樣計算延時以及PI 控制,通過后續推導驗證設計的主動阻尼控制策略與實際被動阻尼控制策略的效果相同。

2.1 LCL 并網逆變器無源阻尼控制策略

無源阻尼控制策略較容易實現,而且不需要增加額外的控制手段,但需要引入阻尼電阻R,導致有功功率的損耗,通常情況下在逆變器電感上串聯R;考慮到三相LCL 的狀態反饋跟單相LCL 狀態反饋大致類似,為了簡化本文就以單相LCL 進行數學建模,該控制方法拓撲圖及其結構框圖如圖3 所示。

圖3 無源阻尼拓撲和結構框圖

按照圖3 所示,以逆變器側電壓Vinv為輸入,網側電流ig為輸出,則LCL 傳遞函數為

由式(9),可以得到阻尼比為

根據式(9)可知,R 值越大,其阻尼控制的效果越好。當R 值過大時,LCL 對高頻諧波的抑制效果明顯,同時R 產生的有功損耗會進一步增大;當R 值較小時,其抑制固有諧振的能力變弱。因而確定R 值需綜合考慮對濾波器的固有諧振抑制以及高次諧波的抑制要求。

2.2 LCL 并網逆變器有源阻尼控制策略

圖4 為三相LCL 并網逆變器有源阻尼控制圖,逆變器側電感為L1、電網側電感為L2、濾波電容為C以及電網阻抗為Lg,忽略所有的等效串聯電阻確保系統無被動阻尼,其中LT=L2+Lg,測量電網電壓以進行鎖相環同步,測量逆變器側電流并向有源阻尼提供一個附加阻尼GS的控制器,Gc(s)為PI 控制器以跟蹤參考電流,而且還考慮到采樣計算調制會帶來Gd(s)=e-1.5sTs的延時。

圖4 三相LCL 并網逆變器有源阻尼控制圖

根據圖4 畫出其有源阻尼控制框圖,如圖5 所示,考慮到d 軸跟q 軸一致,本文只拿d 軸進行數學分析,并考慮到逆變器增益KPWM。

圖5 三相LCL 并網逆變器有源阻尼控制框圖

按照圖5 所示,以逆變器側電壓Vinv為輸入,網側電流iL2為輸出,則GLCL1(s)傳遞函數為

從式(11)可以明顯看出,逆變器側電流的有源阻尼控制策略存在兩個阻尼源,一個是系統本身固有阻尼kp,另一個是本文加入的附加阻尼GS。

主動阻尼控制策略若想與實際被動阻尼控制策略產生諧振抑制效果一致,式(9)應與式(11)等價。為簡化計算,本文先忽略時間延遲,阻尼增益K=GS+kp和阻尼電阻R 之間的關系可用下式表示,即

當考慮延遲的時候,式(11)分母二次項引入指數e-1.5sTs將影響到開環極點的分布,為了將有源阻尼控制策略與無源阻尼控制策略等效,需要等同于在實際被動阻尼電阻R 上再串聯一個虛擬阻抗Z(s),如圖6 所示。

圖6 串聯阻抗后無源阻尼等效控制框圖

Z(s)可表示為

把s=jw 代入,則Z(s)化簡為

根據上述公式可知,虛擬阻抗Z(s)由虛擬電阻器Req和虛擬電抗器jXeq表示。將電感L1與jXeq合并,形成等效電感器Leq,即

串聯虛擬阻抗后,實際無源阻尼的LCL 諧振頻率wr'為

根據式(10)得出串聯阻抗后無源阻尼的阻尼因子為

通過前面分析得出阻尼因子應該是越大越好。

比例增益kp取決于串聯電感、交叉頻率wc(一般取值為0.3 wr)和直流母線電壓,kp和積分時間常數計算公式分別為

計算得到kp=0.06,ki=0.004 4,根據前面分析,當K=0.092 時,存在最佳阻尼因子,附加阻尼Gs=K-kp=0.032;把計算后的參數代入進行仿真;聯立式(13)、(15)、(16)、(17)、(18)、(19),作出阻尼增益K 與阻尼因子ζ 的函數關系圖,如圖7 所示。

圖7 阻尼增益與阻尼因子函數關系圖

由圖7 可看出,當K=0.092 時,電路處于最佳阻尼因子狀態,可以提供最佳的阻尼能力和更大的穩定裕度。隨著Lg變化,電路由剛性網絡變換成弱性網絡,K 的穩定區域擴大,最佳阻尼因子從0.141 提高到0.4,但最佳阻尼因子K 沒有變化。由圖8 可以看出,隨著電網阻抗增加,LCL 諧振峰值變更加平滑,且增益裕度一直大于0 dB,系統保持穩定;證明該有源阻尼控制策略在剛性網絡和弱性網絡都具有良好的魯棒性。

圖8 不同Lg 下頻率響應曲線

3 級聯式重復PI 控制設計

在電網外部,由于諧波的存在,將會嚴重影響系統的穩定,相比于偶次諧波,奇次諧波會對系統造成更大的危害[13]。為了有效地降低奇次諧波,縮短重復控制器開始修正跟蹤誤差時的延遲時間,筆者將內膜進行改進,內模只選取周期采樣次數的1/2。奇次諧波重復控制器是重復控制將一個基本單元放入回路中,該基本單元基于采樣周期延遲的正反饋,由N/2 個采樣周期延遲的負反饋組成,其中N=Ts/T0(N 為一個周期的采樣次數,Ts為系統的采樣周期,T0為輸出正弦電壓基波周期),可以在某個頻率及其所有諧波下引入無限增益。改進后的重復控制結構如圖9b 所示。

圖9 傳統和奇次諧波重復單元

PI 控制器中,比例環節可以實現動態響應速度快,積分環節可以實現無靜差跟蹤,但PI 控制很容易出現周期性干擾,沒有很好的抑制作用;重復控制由于內部存在延時環節,動態性能較差,但在周期信號控制中具有良好的穩態性能[14]。本文將二者結合在一起,并考慮到降低奇次諧波,采用基于奇次內膜級聯式復合重復控制結構,如圖10 所示。重復控制器與反饋級聯閉環,第1 個周期內通過閉環進行調節,從第2 個周期開始,重復控制再進行無靜差跟蹤,以提高系統的動態響應。

圖10 基于奇次諧波級聯式復合重復控制框圖

圖10 中,z-N/2為超前環節;內膜包括低通濾波器和延遲環節Q(z)z-N/2;C(z)為補償器,提供相位補償和幅值補償;Gc(z)為PI 控制器;Gp(z)為被控對象;r(z)為參考電流;d(z)為擾動;y(z)為輸出電壓;E(z)為跟蹤誤差。

對重復控制器的設計即對三個環節參數低通濾波器Q(z)、采樣次數N 和補償環節C(z)=krzMS(z)的設計。為了簡化設計過程,通常設Q(z)=0.95、kr=0.9,則N=10 000/50,即只用對濾波器補償S(z)補償和zM超前環節進行設計即可。

三相LCL 并網逆變器的系統傳遞函數為式(11),在重復控制器加入之前已經插入PI 控制器,經過前面計算,比例函數kp=0.06,ki=0.004 4。

將三相LCL 并網逆變器通過雙線性變換法離散化后,得到

S(z)一般設計為低通濾波器,二階濾波器起的主要是對高頻信號的衰減作用,二階濾波器僅提供高頻衰減,因此可以根據二階濾波器的截止頻率來設計,補償器C(z)中的濾波器S(z)一般采用的是二階低通濾波器的形式,如式(17)所示,取低通濾波器的截止頻率為受控對象的諧振頻率ω=31 400 rad/s,阻尼系數0.707,通過雙線性比變換法離散后,最后設計的補償器C(z)中的濾波器S(z)為:

從圖11 中可看出高頻段增益得到了有效抑制。因為濾波器和逆變器本身都存在滯后,它則由超前環節zM來補償,首先畫出S(z)GLCL1(z)的Bode 圖,再依次畫z-M的Bode 圖;通過伯德圖比較分析模塊z-M 與補償對象之間的相位滯后特性來完成M 參數的選擇,以達到最佳的相位補償效果。圖12 為不同階數的超前補償效果圖。

圖11 S(z)GLCL1(z)的伯德圖

圖12 超前環節z-M 的相頻特性

z-2的相頻特性曲線在中低頻段與被控對象的相頻特性曲線基本一致??梢哉J為提前2 拍的超前補償環節z2可以有效補償和糾正系統的相位滯后問題,并保證系統中低頻零相移的控制特性。

4 頻率自適應級聯式重復PI 控制

傳統重復控制的離散域內膜G(z)=1/(1-Q*z-N/2),其中N=fs/f0(fs為采樣頻率,f0為電網頻率)。在理想情況下,采樣頻率和電網頻率的比值N 為整數,此時重復控制能夠有效抑制各次諧波,但由于弱電網存在,電網頻率會出現一定程度的波動,N 中會出現小數部分,即

其中,Ni=int(N)表示整數部分,F=N-Ni表示小數(0

由拉格朗日插值原理可得

由式(25)可以得出基于冪級數展開的分數延遲濾波器,如圖13 所示。

圖13 M 階分數延遲濾波器

當F=d,d 為整數時,式(25)依然成立,即

用矩陣形式表示為

其中

利用拉格朗日插值法計算得到

當M 分別取1、2、3 時,濾波器的系數如表1 所示。

表1 分數延遲濾波器的系數

本文選取M=3 作為實際的插值階數。

在實際生活中,電網頻率需要在49.6 Hz~50.4 Hz 之間波動,以7 次諧波為例,當電網頻率為49.6 Hz 和50.4 Hz 時,傳統內膜(N=200)和加入分數延遲環節內膜(N=10 000/49.6=201.6 和10 000/50.4=198.4)的幅值特性對比如圖14 所示。頻率波動時加入分數延遲環節的RC 在7 次諧波處的增益由9.71 dB 增加到26 dB,諧波抑制效果更好。

圖14 7 次諧波處的幅頻特性

對應的控制框圖如圖15 所示。

圖15 頻率自適應級聯式重復PI 控制框圖

5 仿真驗證

考慮到電流紋波限制、無功功率要求和電網諧波標準,表2 列出了主要的電路參數。

表2 電路主要參數

對變化的光照模式下COOT 優化算法的追蹤能力進行仿真,設定在時間t=0 s 時,光伏組件受光照強度的情況1 為:1 000、900、600、800、400 W/m2,最大功率為Pmax=4 250 W;在t=0.6 s 時所受光照情況2為:600、500、200、400、200 W/m2,最大功率為Pmax=3 430 W;變化遮擋下的輸出有功功率、無功功率曲線如圖16 所示。

圖16 變化光照下有功、無功功率曲線

光照情況1 時,基于COOT 算法的MPPT 系統在0.42 s 后追蹤到全局GMPP,為4 250 W;t=0.6 s時,光照情況2 變化,COOT 算法重啟,并在0.82 s 時準確追蹤到新的全局GMPP,為3 430 W??芍贑OOT 方法的MPPT 系統在變化的光照模式下可以快速準確跟蹤到全局最大功率點處。

將頻率由原來的工頻50 Hz 變為49.6 Hz,采用冪級數展開的分數延遲濾波器下基于奇次內膜級聯式復合重復控制并網策略(即頻率自適應級聯式重復PI 控制并網策略),系統輸出電壓-電流波形和電流諧波分析圖如圖17 所示。

圖17 49.6 Hz 下改進型重復PI 控制電壓-電流波形和并網電流諧波含量

將頻率由原來的工頻50 Hz 變為50.4 Hz,采用冪級數展開的分數延遲濾波器下基于奇次內膜級聯式復合重復控制并網策略(即頻率自適應級聯式重復PI 控制并網策略),系統輸出電壓-電流波形和電流諧波分析圖如圖18 所示。

圖18 50.4 Hz 下改進型重復PI 控制電壓-電流波形和并網電流諧波含量

由圖18 可知,無論是光照強度還是電網頻率發生變化時,并網電流與電網電壓頻率和相位依舊一致;并且工頻由50 Hz 變為49.6 Hz,輸出電流的THD 僅為1.56%;工頻由50 Hz 變為50.4 Hz,輸出電流的THD 僅為1.31%;說明無論是光照強度發生變化或者頻率發生漂移時,采用頻率自適應級聯式重復PI 控制可使系統的輸出電流對參考電流的跟蹤誤差小,而且可以明顯看出奇次諧波含量也很少。

6 結語

本文將白冠雞優化算法結合到MPPT 上面,可以快速準確追蹤到全局最優功率值沒有陷入到局部最優;本文僅采用逆變器側電流反饋進行有源阻尼控制,不同于利用伯德圖進行分析,提出基于阻尼因子的設計方法,在數學建模過程中把采樣計算延時和PI 控制加入,驗證得到的數學模型與無源阻尼的數學模型是等效的,最后加入重復控制;在這個基礎上,為了減少奇次諧波,提出基于奇次內膜的級聯式復合重復控制策略降低奇次諧波的危害;考慮到實際電網頻率會出現一定波動,提出冪級數展開的分數延遲濾波器,提高重復內膜對電網頻率擾動的適應性;最后把它們運用到三相LCL 型光伏并網系統Simulink 仿真中,通過將光照強度和工頻頻率進行變化,驗證本文提出的基于COOT 優化算法和改進型重復PI 控制策略下,并網電壓和電流具有良好的正弦形、電流可以很快跟蹤到給定值、THD 值較小且奇次諧波含量低。

猜你喜歡
跟隨者領導者阻尼
N維不可壓無阻尼Oldroyd-B模型的最優衰減
關于具有阻尼項的擴散方程
具有非線性阻尼的Navier-Stokes-Voigt方程的拉回吸引子
閉目塞聽,才是領導者的第一大忌
由城市臺的“跟隨者”到縣域“三農”媒體的 “領導者”
從“跟隨者”到“引領者”
—— 甕福集團PPA項目成為攪動市場的“鯰魚”
真誠是領導者的最高境界
跟隨者
具阻尼項的Boussinesq型方程的長時間行為
出口跟隨者會受益于開拓者嗎?——來自中國工業企業的證據
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合