劉福國 國欽光,殷炳毅,王守恩
(1.國網山東省電力公司電力科學研究院,山東 濟南 250003;2.山東電力研究院,山東 濟南 250003;3.華電濰坊發電有限公司,山東 濰坊 261000)
檢測發電機組效率時,鍋爐燃煤發熱量還無法連續測量,只能通過采集入爐煤樣品,在實驗室化驗分析得到[1]。用少量樣品煤的化驗數據表示檢測周期內入爐煤發熱量,其不確定性來自樣品采集和實驗室分析等2個環節[2],采樣環節引起的不確定度占95%以上[3]。文獻[4]研究了鍋爐采樣樣品煤發熱量的概率分布特性,對單個樣品發熱量的不確定度進行了評定。在實際生產中,機組效率常用供電煤耗表示,電廠日常供電煤耗檢測周期有時長達1個月、1個季度甚至1年。檢測周期內的燃煤發熱量一般采用多個樣品煤發熱量的平均值[5]。因此,多次采樣樣品發熱量均值的不確定度評定對機組效率檢測具有重要意義。
μx=μ
(1)
(2)
當總體N為正態分布或樣本容量n>30時,樣本均值的抽樣分布可作為正態分布,正態分布的參數μx、σx根據式(1)、式(2)計算,故可對n個樣品發熱量均值的不確定度進行評定。
但當總體N為非正態分布且樣本容量n<30時,樣本均值的抽樣分布為非正態分布,此時,如何對n個樣品發熱量均值的不確定度進行評定,是本文要解決的問題。
非正態總體N通常來自實際采樣實驗,為已知數據。要對n個樣品發熱量均值的不確定度進行評定,首先獲取樣本均值的抽樣樣本,傳統方法[9]是從總體的N個發熱量中隨機抽取一個容量為n的樣本,計算該樣本均值,采用同樣方式進行多次抽樣,得到樣本均值的抽樣樣本,利用該樣本對均值的不確定度進行評定。從總體的N個發熱量中抽取容量為n的樣本,重復抽樣總計可得到Nn個樣本,因此樣本均值抽樣樣本的最大容量為Nn。
獲得樣本均值的抽樣樣本的另一種方法是基于蒙特卡洛法[10,11]的隨機采樣實驗,與傳統抽樣方法相比,隨機采樣試驗獲得的虛擬樣本容量是任意的,采用隨機采樣試驗對n個樣品發熱量均值的不確定度進行評定的步驟為:
(1) 計算總體N的發熱量概率密度。
(2) 將總體N發熱量概率密度表示成多項式函數。
(3) 建立虛擬采樣系統。按照總體N的概率密度產生隨機數,作為虛擬采樣樣品的發熱量。這種利用給定概率密度函數進行隨機樣抽的系統,看作是樣品發熱量的虛擬采樣系統。
(4) 利用虛擬采樣系統進行n次采樣,得到n個虛擬樣品的發熱量數據,計算發熱量的平均值,得到第1個發熱量均值;采用同樣的方法得到m個發熱量均值,作為發熱量均值的抽樣樣本。
(5) 計算發熱量均值樣本的平均值、標準差及概率密度分布,在給定置信概率下求取包含系數,得到n個樣品發熱量均值擴展不確定度。
以某電廠鍋爐為例,對入爐煤發熱量均值不確定度進行評定。
表1為某電廠2017年5~7月入爐煤實際采樣試驗樣本N,該樣本包含276個樣品煤的發熱量[12]。這些發熱量數據既包含實驗室化驗分析引起的不確定性,又包含采制樣本環節引起的不確定性。研究表明,實驗室分析環節引起的不確定度約為0.05 MJ/kg[13~15],而采制樣本環節引起的不確定度為0.5~2.5 MJ/kg[2,4,12]。因此,忽略實驗室分析環節的不確定度,可認為表1中發熱量數據的不確定性是由采制樣本引起的。利用表1的數據對n次采樣發熱量均值的不確定度進行評定。
表1給出的樣本平均值μN=19.686 MJ/kg,標準差σN=0.933 MJ/kg。
表1 總體樣本N的發熱量Tab.1 Calorific value of sample N MJ/kg
總樣本N中,發熱量最小值為16.894 MJ/kg,最大值為21.217 MJ/kg。將此區間等分成若干個子區間,統計每個子區間內的樣品個數,計算落入區間的概率,概率密度是概率的變化率。利用Matlab統計工具箱中的hist函數可以統計落入每個子區間內的樣品數,利用ksdensity函數可以對樣本N的概率密度進行核心平滑估計[16]。由圖1給出的通過區間統計得到的概率密度的直方圖以及概率密度平滑估計可以看出,樣本發熱量不服從正態分布。概率密度平滑估計結果可用8階多項式函數f表示為:
圖1 總體樣本N的概率密度Fig.1 Probability density of population sample N
(3)
對于總體樣本N,式(3)所示的多項式函數由分段函數f1和f2組成,在不同區間上,函數f1和f2的系數見表2。
表2 多項式函數f的系數Tab.2 Coefficient of polynomial function f
根據概率密度函數式(3),采用挑選抽樣方法產生發熱量樣本,這種方法由馮·諾依曼在1951年提出[17]。如圖2所示,當已知概率密度函數f(Q)時,利用挑選抽樣法產生發熱量隨機數的主要步驟為:
(1) 在發熱量取值區間[a,b]內,產生均勻隨機數Q。
Q=(b-a)r1+a
(4)
式中:r1為區間[0,1]均勻隨機數,r1∈[0,1]。
(2) 在[0,fmax]區間內,產生均勻隨機數f。
f=r2fmax
(5)
式中:r2為區間[0,1]均勻隨機數,r2∈[0,1];fmax為函數f(Q)在區間[a,b]上的最大值。
(3) 當f≤f(Q)時,接受Q為所選的隨機數,否則,返回第(1)步,重新抽取1對(Q,f)。
根據圖2,挑選樣抽法的幾何解釋為:隨機選取位于矩形abcd內的點(Q,f),選擇位于曲線f(Q)下面的點,它們對應的發熱量Q服從概率密度為f(Q)的分布。
圖2 根據給定概率密度函數f(Q)產生隨機數QFig.2 Generating random numbers Q based on a given probability density function f(Q)
要對n次采樣發熱量均值的不確定度進行評定,首先利用第2.3節的虛擬采樣技術獲得均值樣本。按式(3)給定的概率密度函數,產生n個發熱量隨機數,計算它們的平均值,作為發熱量均值的第1個采樣數據。用相同方式得到m個發熱量均值的采樣數據,作為發熱量均值樣本。
取n=4和n=50,采用虛擬采樣方法分別得到發熱量均值樣本,樣本容量m取為m=10 000。對樣本進行統計計算,n=4時,樣本平均值μ4=19.689 MJ/kg,標準差σ4=0.480 MJ/kg;n=50時,樣本平均值μ50=19.689 MJ/kg,標準差σ50=0.135 MJ/kg。利用ksdensity函數對2個樣本的概率密度進行核心平滑估計,n=4時,結果見圖3,其概率密度可擬合成式(3)所示的多項式函數f(Q),限于篇幅,這里不再給出多項式系數;n=50時,區間統計得到的概率密度直方圖以及函數ksdensity的平滑估計結果見圖4,發熱量均值樣本非常接近正態分布,圖4中還給出了正態分布N(19.686,0.1352)的概率密度,可以看出,它和平滑估計結果吻合良好。
圖3 n=4時發熱量平均值的概率密度分布Fig.3 Probabilistic density distribution of the mean calorific value at n=4
圖4 n=50時發熱量平均值的概率密度分布與正態分布的對比Fig.4 Comparison of probability density of average calorific value with normal distribution at n=50
如圖3,對n=4時發熱量均值的不確定度進行評定,就是求給定的置信概率ε下的包含因子k,根據置信概率的定義[18]得到關于包含因子k的方程如式(6)。
求解式(6),可得到包含因子k。求解時取置信概率ε=0.95,式(6)中μ4=19.689,σ4=0.480。采用迭代法求解式(6),得到k=1.934。因此,n=4時發熱量均值的擴展不確定度為kσ4=0.928 MJ/kg,此時發熱量均值結果可表示為19.689±0.928 MJ/kg。
(6)
根據虛擬采樣方法,同樣可對n=50時發熱量均值的不確定度進行評定。利用Matlab的jbtest函數對n=50時發熱量均值樣本進行正態性檢驗[19]。結果表明,均值樣本服從正態分布,置信概率ε=0.95對應的包含因子k=1.96[20,21],發熱量均值的擴展不確定度為kσ50=0.265 MJ/kg。因此,n=50時發熱量均值的采樣化驗結果可表示為19.689±0.265 MJ/kg。
(1) 在發電機組效率檢測時,鍋爐燃煤多次采樣樣品發熱量均值的不確定度評定對于機組效率檢測具有重要意義。
(2) 已知燃煤發熱量樣本,要對n個樣品發熱量均值的不確定度進行評定,當發熱量樣本正態分布或n>30時,則n個樣品發熱量均值服從正態分布,該正態分布的參數、根據式(1)和式(2)計算,據此可對n個樣品發熱量均值的不確定度進行評定。
(3) 當發熱量樣本不服從正態分布且n<30時,要對n個樣品發熱量均值的不確定度進行評定,可采用挑選樣抽法進行虛擬采樣,得到n個樣品發熱量的均值樣本,然后根據該均值樣本的概率密度分布對發熱量均值的不確定度評定。
(4) 采用虛擬采樣均值樣本對n個樣品發熱量均值的不確定度評定結果與根據中心極限定理得到的結論一致。