配煤動力學特性回歸模型的R語言實現

之前在研究生畢業論文當中用SPSS採用多元線性及非線性模型用於預測配煤動力學特性，本節主要結合R語言提供的函數和工具重現一下。

1 幾個概念：

單煤：就是從一個礦區開採出來的煤炭，一般性質相近；

配煤：為了取得特定的特性，由不同種單煤按照不同的比例混合而成，由於煤炭性質的複雜性，通常配煤的動力學特性無法由各種單煤加權平均求得；

煤的工業分析（煤的工業分析_百度百科）：煤的工業分析是指包括煤的水分（M ）、灰分（A ）、揮發分（V ）和固定碳（Fc )四個分析項目指標的測定的總稱。煤的工業分析是了解煤質特性的主要指標，也是評價煤質的基本依據。通常煤的水分、灰分、揮發分是直接測出的，而固定碳是用差減法計算出來的。廣義上講，煤的工業分析還包括煤的全硫分和發熱量的測定， 又叫煤的全工業分析。本文的實驗數據介於狹義與廣義之間，測量了發熱量，但是全硫在元素分析中進行。

煤的元素分析（煤的元素分析_百度百科）：煤的元素分析是對煤中的元素含量進行檢測和分析（一般用質量百分數表示），包括常規的C、H、O、N、S、Al、Si、Fe、Ca等元素含量，還可檢測煤中的痕量元素包括Ti、Na、K等。本文的實驗數據主要測量了最主要的C、H、O、N、S五種。

煤的著火溫度：煤的著火點是隨著熱力條件變化而變化的，並不是一個物理常數，只是在一定條件下得出的相對特徵值。通常把由緩慢氧化狀態轉變為高速燃燒狀態的瞬間過程稱為著火，此時的溫度即稱為著火溫度（Ignition Temperature）。

煤的活化能（Activation Energy）：即讓某一種煤開始發生化學反應（通常指燃燒）所需要克服的能量障礙。

2 數據描述：

有16種單煤及48種這些單煤按不同比例混合而成配煤的工業分析以及元素分析數據，並通過熱重分析、熱重微分分析及差熱分析得出的曲線結合Arrhenius方程及Doyle單曲線法求得各個不同煤樣的著火溫度（）和活化能。

這樣一共64組數據。

考慮到從實驗的角度求取著火溫度及活化能需要特殊的儀器並耗費較長的時間，希望能基於上述64組數據找到一個用於預測工業分析、元素分析已知煤樣的動力學特性，包括著火溫度、活化能等，特別是在混配煤動力學特性預測方面進行應用。

3 數據導入：

> library(xlsx)n載入需要的程輯包：rJavan載入需要的程輯包：xlsxjarsnn載入程輯包：『xlsx』nnThe following objects are masked from 『package:openxlsx』:nn createWorkbook, loadWorkbook, read.xlsx, saveWorkbook,n write.xlsxnn> filepath <- "G:/DataCruiser/workspace/Coal Analysis/Coal Analysis data.xlsx"n> CoalData <- read.xlsx(filepath, 1, encoding = "UTF-8")n> class(CoalData)n[1] "data.frame"n> head(CoalData)n 煤樣編號 煤樣名稱 Ignition_Temp Activation_Ene Mad. Aad. Vad. FCad.n1 1 黃陵 480 216035 3.85 23.68 26.98 45.49n2 2 劉橋 551 301993 1.00 14.62 10.44 73.94n3 3 大同 480 223548 2.50 14.42 26.06 57.02n4 4 義馬 403 160840 9.47 14.19 30.02 46.32n5 5 兗州新 466 226362 2.21 18.26 29.39 50.14n6 6 淮南 468 231235 1.62 25.47 23.32 49.59n Qnet.ad.J.g. Cad. Had. Nad. St.add. Oad.n1 24493.9 58.93 3.42 1.12 0.85 8.15n2 30646.5 75.90 3.14 1.15 0.42 3.81n3 28751.1 68.55 3.92 0.81 0.69 9.11n4 23802.3 58.44 2.92 0.83 0.44 13.71n5 27101.0 63.13 3.74 1.18 0.71 10.77n6 25702.6 61.26 3.32 1.09 0.52 6.72n

4 數據重命名：

著火溫度：ignitionTemperature

活化能：activationEnergy

水分：moisture

灰分：ash

揮發分：volatileMatter

固定碳：fixedCarbon

發熱量：LHV

碳：carbon

氫：hydrogen

氮：nitrogen

硫：sulfur

> names(CoalData) <- c("sampleNO","sampleName","ignitionTemperature","activationEnergy","moistrue",n+ "ash","volatileMatter","fixedCarbon","LHV","carbon","hydrogen","nitrogen","sulfur","oxygen")n> head(CoalData)n sampleNO sampleName ignitionTemperature activationEnergy moistrue ashn1 1 黃陵 480 216035 3.85 23.68n2 2 劉橋 551 301993 1.00 14.62n3 3 大同 480 223548 2.50 14.42n4 4 義馬 403 160840 9.47 14.19n5 5 兗州新 466 226362 2.21 18.26n6 6 淮南 468 231235 1.62 25.47n volatileMatter fixedCarbon LHV carbon hydrogen nitrogen sulfurn1 26.98 45.49 24493.9 58.93 3.42 1.12 0.85n2 10.44 73.94 30646.5 75.90 3.14 1.15 0.42n3 26.06 57.02 28751.1 68.55 3.92 0.81 0.69n4 30.02 46.32 23802.3 58.44 2.92 0.83 0.44n5 29.39 50.14 27101.0 63.13 3.74 1.18 0.71n6 23.32 49.59 25702.6 61.26 3.32 1.09 0.52n oxygenn1 8.15n2 3.81n3 9.11n4 13.71n5 10.77n6 6.72n

4 數據分組：

按不同因變數（ignitionTemperature和activationEnergy）進行數據分組。

> IT_Estimate <- as.data.frame(CoalData[,c("ignitionTemperature","moistrue",n+ "ash","volatileMatter","fixedCarbon","LHV","carbon","hydrogen","nitrogen","sulfur","oxygen")])n> head(IT_Estimate)n ignitionTemperature moistrue ash volatileMatter fixedCarbon LHVn1 480 3.85 23.68 26.98 45.49 24493.9n2 551 1.00 14.62 10.44 73.94 30646.5n3 480 2.50 14.42 26.06 57.02 28751.1n4 403 9.47 14.19 30.02 46.32 23802.3n5 466 2.21 18.26 29.39 50.14 27101.0n6 468 1.62 25.47 23.32 49.59 25702.6n carbon hydrogen nitrogen sulfur oxygenn1 58.93 3.42 1.12 0.85 8.15n2 75.90 3.14 1.15 0.42 3.81n3 68.55 3.92 0.81 0.69 9.11n4 58.44 2.92 0.83 0.44 13.71n5 63.13 3.74 1.18 0.71 10.77n6 61.26 3.32 1.09 0.52 6.72n> AE_Estimate <- as.data.frame(CoalData[,c("activationEnergy","moistrue",n+ "ash","volatileMatter","fixedCarbon","LHV","carbon","hydrogen","nitrogen","sulfur","oxygen")])n> head(AE_Estimate)n activationEnergy moistrue ash volatileMatter fixedCarbon LHVn1 216035 3.85 23.68 26.98 45.49 24493.9n2 301993 1.00 14.62 10.44 73.94 30646.5n3 223548 2.50 14.42 26.06 57.02 28751.1n4 160840 9.47 14.19 30.02 46.32 23802.3n5 226362 2.21 18.26 29.39 50.14 27101.0n6 231235 1.62 25.47 23.32 49.59 25702.6n carbon hydrogen nitrogen sulfur oxygenn1 58.93 3.42 1.12 0.85 8.15n2 75.90 3.14 1.15 0.42 3.81n3 68.55 3.92 0.81 0.69 9.11n4 58.44 2.92 0.83 0.44 13.71n5 63.13 3.74 1.18 0.71 10.77n6 61.26 3.32 1.09 0.52 6.72n

5 數據相關性檢驗：

這裡我們採用Pearson相關性分析，需要調用psych包中的corr.test()函數，結果如下：

> corr.test(IT_Estimate, use = "complete")nCall:corr.test(x = IT_Estimate, use = "complete")nCorrelation matrix n ignitionTemperature moistrue ash volatileMatternignitionTemperature 1.00 -0.58 0.09 -0.65nmoistrue -0.58 1.00 -0.42 0.31nash 0.09 -0.42 1.00 -0.20nvolatileMatter -0.65 0.31 -0.20 1.00nfixedCarbon 0.55 -0.12 -0.59 -0.65nLHV 0.39 -0.14 -0.79 -0.14ncarbon 0.33 -0.06 -0.82 -0.22nhydrogen -0.19 -0.33 -0.25 0.60nnitrogen 0.22 -0.49 -0.23 0.00nsulfur -0.07 -0.11 0.29 0.30noxygen -0.57 0.61 -0.29 0.67n fixedCarbon LHV carbon hydrogen nitrogen sulfurnignitionTemperature 0.55 0.39 0.33 -0.19 0.22 -0.07nmoistrue -0.12 -0.14 -0.06 -0.33 -0.49 -0.11nash -0.59 -0.79 -0.82 -0.25 -0.23 0.29nvolatileMatter -0.65 -0.14 -0.22 0.60 0.00 0.30nfixedCarbon 1.00 0.80 0.86 -0.16 0.32 -0.44nLHV 0.80 1.00 0.94 0.26 0.46 -0.31ncarbon 0.86 0.94 1.00 0.24 0.51 -0.39nhydrogen -0.16 0.26 0.24 1.00 0.58 0.26nnitrogen 0.32 0.46 0.51 0.58 1.00 -0.06nsulfur -0.44 -0.31 -0.39 0.26 -0.06 1.00noxygen -0.40 -0.16 -0.27 0.04 -0.40 -0.04n oxygennignitionTemperature -0.57nmoistrue 0.61nash -0.29nvolatileMatter 0.67nfixedCarbon -0.40nLHV -0.16ncarbon -0.27nhydrogen 0.04nnitrogen -0.40nsulfur -0.04noxygen 1.00nSample Size n[1] 64nProbability values (Entries above the diagonal are adjusted for multiple tests.) n ignitionTemperature moistrue ash volatileMatternignitionTemperature 0.00 0.00 1.00 0.00nmoistrue 0.00 0.00 0.02 0.38nash 0.47 0.00 0.00 1.00nvolatileMatter 0.00 0.01 0.11 0.00nfixedCarbon 0.00 0.36 0.00 0.00nLHV 0.00 0.28 0.00 0.27ncarbon 0.01 0.62 0.00 0.08nhydrogen 0.13 0.01 0.05 0.00nnitrogen 0.08 0.00 0.06 0.99nsulfur 0.56 0.39 0.02 0.02noxygen 0.00 0.00 0.02 0.00n fixedCarbon LHV carbon hydrogen nitrogen sulfurnignitionTemperature 0.00 0.04 0.25 1.00 1.00 1.00nmoistrue 1.00 1.00 1.00 0.25 0.00 1.00nash 0.00 0.00 0.00 0.94 1.00 0.46nvolatileMatter 0.00 1.00 1.00 0.00 1.00 0.40nfixedCarbon 0.00 0.00 0.00 1.00 0.31 0.01nLHV 0.00 0.00 0.00 0.79 0.01 0.33ncarbon 0.00 0.00 0.00 1.00 0.00 0.05nhydrogen 0.22 0.04 0.05 0.00 0.00 0.80nnitrogen 0.01 0.00 0.00 0.00 0.00 1.00nsulfur 0.00 0.01 0.00 0.04 0.62 0.00noxygen 0.00 0.22 0.03 0.77 0.00 0.74n oxygennignitionTemperature 0.00nmoistrue 0.00nash 0.51nvolatileMatter 0.00nfixedCarbon 0.03nLHV 1.00ncarbon 0.68nhydrogen 1.00nnitrogen 0.03nsulfur 1.00noxygen 0.00nn To see confidence intervals of the correlations, print with the short=FALSE optionn> corr.test(AE_Estimate, use = "complete")nCall:corr.test(x = AE_Estimate, use = "complete")nCorrelation matrix n activationEnergy moistrue ash volatileMatternactivationEnergy 1.00 -0.40 -0.03 -0.56nmoistrue -0.40 1.00 -0.42 0.31nash -0.03 -0.42 1.00 -0.20nvolatileMatter -0.56 0.31 -0.20 1.00nfixedCarbon 0.54 -0.12 -0.59 -0.65nLHV 0.39 -0.14 -0.79 -0.14ncarbon 0.37 -0.06 -0.82 -0.22nhydrogen -0.22 -0.33 -0.25 0.60nnitrogen 0.21 -0.49 -0.23 0.00nsulfur -0.17 -0.11 0.29 0.30noxygen -0.46 0.61 -0.29 0.67n fixedCarbon LHV carbon hydrogen nitrogen sulfur oxygennactivationEnergy 0.54 0.39 0.37 -0.22 0.21 -0.17 -0.46nmoistrue -0.12 -0.14 -0.06 -0.33 -0.49 -0.11 0.61nash -0.59 -0.79 -0.82 -0.25 -0.23 0.29 -0.29nvolatileMatter -0.65 -0.14 -0.22 0.60 0.00 0.30 0.67nfixedCarbon 1.00 0.80 0.86 -0.16 0.32 -0.44 -0.40nLHV 0.80 1.00 0.94 0.26 0.46 -0.31 -0.16ncarbon 0.86 0.94 1.00 0.24 0.51 -0.39 -0.27nhydrogen -0.16 0.26 0.24 1.00 0.58 0.26 0.04nnitrogen 0.32 0.46 0.51 0.58 1.00 -0.06 -0.40nsulfur -0.44 -0.31 -0.39 0.26 -0.06 1.00 -0.04noxygen -0.40 -0.16 -0.27 0.04 -0.40 -0.04 1.00nSample Size n[1] 64nProbability values (Entries above the diagonal are adjusted for multiple tests.) n activationEnergy moistrue ash volatileMatter fixedCarbonnactivationEnergy 0.00 0.04 1.00 0.00 0.00nmoistrue 0.00 0.00 0.02 0.38 1.00nash 0.80 0.00 0.00 1.00 0.00nvolatileMatter 0.00 0.01 0.11 0.00 0.00nfixedCarbon 0.00 0.36 0.00 0.00 0.00nLHV 0.00 0.28 0.00 0.27 0.00ncarbon 0.00 0.62 0.00 0.08 0.00nhydrogen 0.08 0.01 0.05 0.00 0.22nnitrogen 0.10 0.00 0.06 0.99 0.01nsulfur 0.19 0.39 0.02 0.02 0.00noxygen 0.00 0.00 0.02 0.00 0.00n LHV carbon hydrogen nitrogen sulfur oxygennactivationEnergy 0.04 0.09 1.00 1.00 1.00 0.01nmoistrue 1.00 1.00 0.25 0.00 1.00 0.00nash 0.00 0.00 0.94 1.00 0.46 0.51nvolatileMatter 1.00 1.00 0.00 1.00 0.40 0.00nfixedCarbon 0.00 0.00 1.00 0.31 0.01 0.03nLHV 0.00 0.00 0.79 0.01 0.33 1.00ncarbon 0.00 0.00 1.00 0.00 0.05 0.68nhydrogen 0.04 0.05 0.00 0.00 0.80 1.00nnitrogen 0.00 0.00 0.00 0.00 1.00 0.03nsulfur 0.01 0.00 0.04 0.62 0.00 1.00noxygen 0.22 0.03 0.77 0.00 0.74 0.00nn To see confidence intervals of the correlations, print with the short=FALSE optionn

從相關性檢驗中可以得出無論是著火溫度還是活化能，與水分、揮發分、固定碳、發熱量以及氧這五個變數相關較大，因此，以下考慮以此五個自變數採用多元線性回歸模型分別對著火溫度和活化能進行模擬，首先剔除不相關數據，僅包含我們感興趣的變數。

> AE_Estimate <- as.data.frame(AE_Estimate[,c("activationEnergy","moistrue", "volatileMatter","fixedCarbon","LHV","oxygen")])n> IT_Estimate <- as.data.frame(IT_Estimate[,c("ignitionTemperature","moistrue", "volatileMatter","fixedCarbon","LHV","oxygen")])n> n

用car包中的scatterplotMatrix()函數生成散點圖矩陣。

> library(car)nn載入程輯包：『car』nnThe following object is masked from 『package:psych』:nn logitnn> scatterplotMatrix(IT_Estimate, spread = FALSE, smoother.args = list(lty = 2), main = "Scatter Plot Matrix ignitionTemperature")> scatterplotMatrix(AE_Estimate, spread = FALSE, smoother.args = list(lty = 2), main = "Scatter Plot Matrix of activationEnergy")n

從上圖我們可以看出，著火溫度是一個單峰曲線，隨著水分、揮發分、氧含量的升高而降低，隨著固定碳以及發熱量升高而升高。活化能也有類似的規律。

5 回歸分析

下面使用lm()函數分別擬合各自的多元回歸模型，代碼如下：

> fit_AE <- lm(activationEnergy ~ moistrue + volatileMatter + fixedCarbon + LHV + oxygen, data = AE_Estimate)n> summary(fit_AE)nnCall:nlm(formula = activationEnergy ~ moistrue + volatileMatter + fixedCarbon + n LHV + oxygen, data = AE_Estimate)nnResiduals:n Min 1Q Median 3Q Max n-61685 -16171 613 10830 141303 nnCoefficients:n Estimate Std. Error t value Pr(>|t|) n(Intercept) 223703.974 63279.918 3.535 0.000808 ***nmoistrue -3938.243 3613.322 -1.090 0.280255 nvolatileMatter -6687.855 2682.473 -2.493 0.015536 * nfixedCarbon -3304.308 3005.542 -1.099 0.276134 nLHV 13.360 7.124 1.875 0.065797 . noxygen 1267.353 3105.625 0.408 0.684715 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nnResidual standard error: 33440 on 58 degrees of freedomnMultiple R-squared: 0.465,tAdjusted R-squared: 0.4188 nF-statistic: 10.08 on 5 and 58 DF, p-value: 5.577e-07nn> fit_IT <- lm(ignitionTemperature ~ moistrue + volatileMatter + fixedCarbon + LHV + oxygen, data = IT_Estimate)n> summary(fit_IT)nnCall:nlm(formula = ignitionTemperature ~ moistrue + volatileMatter + n fixedCarbon + LHV + oxygen, data = IT_Estimate)nnResiduals:n Min 1Q Median 3Q Max n-57.360 -8.439 -0.636 7.929 57.860 nnCoefficients:n Estimate Std. Error t value Pr(>|t|) n(Intercept) 499.767993 35.206960 14.195 < 2e-16 ***nmoistrue -5.759033 2.010339 -2.865 0.00580 ** nvolatileMatter -6.467892 1.492444 -4.334 5.9e-05 ***nfixedCarbon -3.753964 1.672189 -2.245 0.02860 * nLHV 0.012422 0.003964 3.134 0.00271 ** noxygen 1.795935 1.727872 1.039 0.30294 n---nSignif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1nnResidual standard error: 18.61 on 58 degrees of freedomnMultiple R-squared: 0.6926,tAdjusted R-squared: 0.6661 nF-statistic: 26.13 on 5 and 58 DF, p-value: 1.03e-13n

採用標準方法對模型進行診斷：

> confint(fit_AE)n 2.5 % 97.5 %n(Intercept) 9.703546e+04 350372.49266nmoistrue -1.117109e+04 3294.60612nvolatileMatter -1.205741e+04 -1318.30177nfixedCarbon -9.320555e+03 2711.93838nLHV -9.011095e-01 27.62129noxygen -4.949232e+03 7483.93717n> confint(fit_IT)n 2.5 % 97.5 %n(Intercept) 429.293606517 570.24237950nmoistrue -9.783163513 -1.73490270nvolatileMatter -9.455342988 -3.48044097nfixedCarbon -7.101214856 -0.40671337nLHV 0.004487498 0.02035647noxygen -1.662777356 5.25464759n> par(mfrow=c(2,2))n> plot(fit_AE)n> plot(fit_AE,main = "活化能診斷圖")n> plot(fit_IT,main = "著火溫度診斷圖")n

我們可以看出總體上數據符合正態性、線性的假設，在活化能的模型當中可能需要適當的加入一些二次項，另外，兩個模型都有部分離群點，可能是實驗的誤差導致。