2021年五一數學建模比賽C題

2021年五?數學建模?賽C題

?章?錄

本?專挑數據挖掘、機器學習和NLP類型的題?做，有興趣也可以逛逛我的數據挖掘競賽專欄。

如果本篇博?對您有所幫助，請不要吝嗇您的點贊??思路

這道題的突破?就在于如何解決第?問，若能夠在第?問便提出?種量化評價風險的?法，下?的題就好解決了。

如何量化呢？這??統(tǒng)計的?法，對每?個感知器的所有時序數據（共5519個），從沖擊性風險、獨?性和偶發(fā)性風險、以及?相關風險三個??，來評價感知器的風險值。

?于第?問，我們可以沿?第?問的?法，以30分鐘為間隔，算出每個感應器的異常得分（得分?第?問?法求出），求這段時間內感應器異常得分的總和，找出前5個異常得分總和最?的時刻即可。并把異常總得分，作為該時刻的異常得分。之后，在每個最異常的時刻上，找出5個異常得分最?的感應器。

?于第三問，我們根據移動平均算法，求出未來15分鐘的風險值即可。

?于第四問，都算出風險了，安全評分還會遠嗎？第?問

從題?可知，數據存在波動，波動分成兩種：

1.正常波動：外界溫度或者產量變化的波動，傳感器誤報，規(guī)律性、獨?性、偶發(fā)性

2.?正常波動：?產過程中的不穩(wěn)定因素，持續(xù)性、聯(lián)動性

注意，正常波動中，有?個規(guī)律性，所以，正常波動是?噪聲嗎？

我們把沒有持續(xù)性和聯(lián)動性的參數，就叫?正常波動。

風險分析

我們來看看題?對風險的定義：

這些異常性波動的出現是?產過程中的不穩(wěn)定因素造成的，預?著可能存在安全隱患，我們視為風險性異常

所以，感應器誤識別，和因溫度等因素引起的獨?性和偶發(fā)性波動，都不是風險。獨?性、偶發(fā)性誤差（?風險也）

這種獨?性、偶發(fā)性的誤差，將原本相對平穩(wěn)的變量，添加上了?噪聲，造成數據帶有這種?噪聲誤差的原因可能是感應器誤識別，溫度、振動等外部環(huán)境：

風險

根據題?，所謂風險應該是有持續(xù)性的、聯(lián)動性的，如下所?：

灰?地帶

規(guī)律性波動

有?些傳感器變量，它的變量波動很?，如下所?：

從題?可知，有些誤差是具有規(guī)律性的，但規(guī)律性，?意味著聯(lián)動性和持續(xù)性。所以，對傳感器2來說，我們可以說是因為溫度等因素，也可以歸咎于?產過程。

但對本題??，我們將類似傳感器2的規(guī)律性波動，視為風險。

沖擊性誤差

如下所?，個?理解認為，感應器10和感應器5中的沖擊性波動，不應該理解為感應器的誤差，?應該視為?產過程中的沖擊性?擾。

1.483×med∣xi?x?∣xi，有：??xi,δ=1.5s?，再次計算：∑(x?i?x?)2/(p1.483×med∣xi?x?∣xi，有：??xi,δ=1.5s?，再次計算：∑(x?i?x?)2/(p?1)??xi,s?x?x?+δx??δx?im，于是沖擊性風險為m/n若xi<x??δx?+δ,其他若xi<x??δx?+δ,其他若xi>x?+δ若xi>x?+δ

算法A來?GB/T6379.5，應?此算法計算得到數據平均值和標準差的穩(wěn)健值，但也可以據此統(tǒng)計沖擊性風險的數據含量。

算法A的計算過程如下所?：?先是計算初始值：

x?=medxis?=

對每個?x?δ,?

x?i=

其中：

x?=∑x?i/ps?=1.134

重復：?x?δ,?

x?i=

直到的第三位有效數字和的對應數字在連續(xù)兩次迭代中不變。

因此，對感應器的數據，我們通過算法A后，得到最終的x?i，記取值為或的的數量為。

+δx??δ的數連續(xù)出現，我們可以給其添加+2，連續(xù)出現3次，記+δx??δ的數連續(xù)出現，我們可以給其添加+2，連續(xù)出現3次，記為1+2+31+1，以此類m′。權重。

如出現以及，記為1分，連續(xù)出現2次，記為1…，若間隔出現2次，則記為推，最終求和得到

于是，沖擊性風險即為：

m′risk1=n

獨?性和偶發(fā)性誤差

獨?性和偶發(fā)性，讓?個原本平滑的時間序列，成為了?組?噪聲。

理想狀態(tài)下，傳感器的數值應為：

加??噪聲后，應為：

μ,σ

或：μ,σ

因此，若?個時間序列，在打亂其時序后，滿?正態(tài)分布，則證明時間序列有?噪聲。?或者，時間序列滿?均勻分布，如上圖傳感器8，則亦可以證明時間序列的風險低。正態(tài)分布檢驗——Anderson–Darling檢驗（原理部分別看）

實踐表明，AD檢驗檢驗正態(tài)分布，似乎?KS檢驗還要差！，所以本??KS檢驗了

Anderson-Darling檢驗是Kolmogorov-Smirnov檢驗的改良，他能夠利?正態(tài)分布的分布函數的良好的數學特性，故?起Kolmogorov檢驗?對所有分布的檢驗，Anderson-Darling檢驗能針對正態(tài)分布，故其檢驗性質更加準確。

AD檢驗的原假設為：樣本服從正態(tài)分布N(μ,σ2)，若未知，則可以?樣本的估計量代替。

∑

i=1(2iYi?σ^=

n?1XiA2A2∞A2∑∞znμ^{?Xˉσ,2=,?

∑?X,∑

i=1(2iYi?σ^=

n?1XiA2A2∞A2∑∞znμ^{?Xˉσ,2=,?

∑?X,otherwise.(i=0μ,ifthemeanisknown.

ni?∑ni=1(Xi?μ)ifthevarianceisnotknown,butthemeanis.?12j1

1n∞)(4i1)e?(4i+1)∑i=1Xotherwise.2n+n22ifthevarianceisknown.π/(8z)∫0e81+w2)?w2(4i+1)2π2/(8z)dw(z

1A2=?n?n?1)(lnΦ(Yi)+ln(1?Φ(Yn?1?i)))

其中為：

XiYi=

均值和?差：

μ^

σ^2=

?為次序統(tǒng)計量，即將原始數據根據升序的?式進?排序。檢驗統(tǒng)計量的臨界值如表所?：

若?于臨界值，則拒絕原假設，即認為樣本不服從正態(tài)分布。

當然，A2在樣本接近的時候，如本例的5000個樣本，則服從：

2πPr(A2<z)=

?exp(?13.436+101.14AD?223.73AD2)A2檢驗supxFn(x)F(x)∑

i=1為指?函數，若有樣本落在范圍內，則為1。經驗函數根據樣本分布，模擬了總體的分布函數（CDF）。是上確界之xiF(x)。F(x)，?exp(?13.436+101.14AD?223.73AD2)A2檢驗supxFn(x)F(x)∑

i=1為指?函數，若有樣本落在范圍內，則為1。經驗函數根據樣本分布，模擬了總體的分布函數（CDF）。是上確界之xiF(x)。F(x)，則Dn

πn2πk=1(?1)k?1e??otherwisen?→∞∣B(F(t))∣∑?2k2x2exp(0.9177?4.279AD?1.38AD2)exp(?8.318+42.796AD?59.938AD2)supt∞=AD≥0.34AD≥0.2xe?(2kk=1?1)2π2/(8x2)?

p?value=

其中AD為的調整：

AD=A2×(1+(.75/n)+2.25/(n2))

因此我們可以計算出相應的p-value。若p-value越接近于1，意味著樣本越有可能是正態(tài)分布，也即風險越?。

參考?獻：

[A^2的分布函數]：EvaluatingtheAnderson-DarlingDistribution作者：GeorgeMarsaglia。

（上述?獻我放到代碼?件?了）均勻分布檢驗——Kolmogorov–Smirnov

Kolmogorov檢驗可以檢驗某個樣本的總體，是否服從給定的分布（任意分布）。當然，?于為什么不?Kolmogorov檢驗來判定正態(tài)性，是因為Kolmogorov檢驗不夠specification，不能完全運?正態(tài)分布的良好數學特性，所以較之Shapiro-Wilk檢驗來說，有效性較低。

但我們可以?Kolmogorov檢驗，來檢驗樣本的總體是否滿?均勻分布。定義檢驗統(tǒng)計量如下：

Dn=∣Fn(x)?F(x)∣

其中為經驗分布，為實際分布，如下：

1Fn(x)=nI[?∞,x](Xi)

I[?∞,x](Xi)supx意。

Kolmogorov檢驗的原假設為：樣本來源于總體分布

根據Glivenko–Cantelli原理，若樣本xi來源于總體分布將收斂于0?；颍?/p>

nD

其中，B(t)服從柯爾莫格羅夫分布，其分布函數如下：

∑∞Pr(K≤x)=1?2

于是，我們也可以根據上述分布，計算出p-value，p-value越?，原假設越不容易被拒絕，也即數據服從均勻分布，風險越?。

?r～tn?2∈[1,8]，也就是0.15s到2min之內的?相關系數，記為[r1,r2,?,r20]，并進?相關i

i=1?r～tn?2∈[1,8]，也就是0.15s到2min之內的?相關系數，記為[r1,r2,?,r20]，并進?相關i

i=1

通過Kolmogorov-Smirnov檢驗，計算出兩個p-value，記為p1,p2。當然，為了避免離群值對檢驗的影響，在均勻分布檢驗時，我們只取四分位數之間的值進?檢驗。

于是可以評價時間序列的??噪聲風險系數為：

risk′2=1?max(p′sw,p′ks)

?相關性系數和?相關檢驗

由于風險信號具有聯(lián)動性和持續(xù)性，所以，我們可以??相關系數來評判，我們的傳感器數據，是否也具備聯(lián)動性和持續(xù)性。?相關系數

?相關系數旨在計算當前時刻，和前k個時刻的Pearson相關系數其計算公式如下：

∑ni=k+1(xt?xˉt)(xt?k?xˉt?k)∑ni=k+1(xt?xˉt)2(xt?k?xˉt?k)2r=

若r越接近于1或-1，意味著數據之間的線性相關性越強。相關性檢驗

同理，我們也可以對數據的?相關性進?統(tǒng)計檢驗，原假設為H0：兩組數據的相關系數為0，即不存在相關性。

在原假設成?的情況下，可以證明的是，檢驗統(tǒng)計量：

n?21?r2

換句話說，只要：

tn?2(α/2)n?2+tn2?2(α/2)∣r∣≤

便可不拒絕原假設，反之拒絕原假設。

對于上述兩個?法，我們考慮計算出數據對k性檢驗檢驗，若原假設不被拒絕，則將對應的r置零，否則保留。

然后根據以下公式計算數據總體的聯(lián)動性、持續(xù)性風險：

∑8risk3=ri×1+2+?+8

上式的系數部分，表??相關系數的持續(xù)時間越長，則權重越?。

參考：《概率論與數理統(tǒng)計》陳希孺，第六章6.4節(jié)注意點

=0

求?相關系數時，有可能會遇到，即使是?條平直的線，也會具有很強的相關性，?如：=0

所以，我們要結合算法A，若計算出來的4分位數相同，則需要將r置0。另外，還要結合ADF檢驗，來判斷時序數據是否穩(wěn)定：

穩(wěn)定性檢驗——AugmentedDickey–FullerTest

也可以同時檢驗上述兩個特征：也即時序數據是否是?噪聲，或者是均勻分布的穩(wěn)定數據。換句話說，時序數據的均值和?差，不會隨著時間的變化?改變。

為此，本?采?：AugmentedDickey–Fuller檢驗

具體細節(jié)由于過于復雜，筆者??也有很多沒搞清楚，所以只能簡要介紹：ADF檢驗的原假設是：時序數據不是穩(wěn)定的，也即數據的均值和?差會隨著時間?變化。

因此，若拒絕了ADF降壓，則讓r風險計算總結

為了量化評判風險，我們?到了：

1.算法A：沖擊性風險

2.正態(tài)、均勻分布檢驗：量化誤差是否屬于風險

3.?相關系數：持續(xù)性和聯(lián)動性風險

結果展?

第?問

第?問是對0:00到23:00進?分析的。為了分析出某個時刻感應器的風險，我們需要以30分鐘為單位，進?分析，從?得出總風險。

如上，以30分鐘為單位，算出每個感應器的異常得分（得分?第?問?法求出），求這段時間內感應器異常得分的總和，找出前5個異常得分總和最?的時刻即可。并把異?？偟梅?，作為該時刻的異常得分。

在每個時刻上，找出5個異常得分最?的感應器：

當前時刻為：0days03:00:00異常的感應器有['感應器2','感應器5','感應器6','感應器7','感應器11']異常得分為156.45667125866223

當前時刻為：0days16:30:00異常的感應器有['感應器6','感應器7','感應器11','感應器12','感應器13']異常得分為156.29759895209494

當前時刻為：0days03:30:00異常的感應器有['感應器2','感應器5','感應器6','感應器7','感應器11']異常得分為156.1732779754605

當前時刻為：0days02:30:00異常的感應器有['感應器2','感應器5','感應器6','感應器7','感應器11']異常得分為155.8447347080813

當前時刻為：0days01:00:00異常的感應器有['感應器57','感應器5','感應器6','感應器7','感應器11']異常得分為155.66803974065186第三問

我們可以考慮兩種?法：

1.?移動平均的?法，以前5個歷史數據的均值，作為當前時刻的異常得分

2.?機器學習的?法，先擬合?個機器學習模型，再?來預測

這?就考慮第?種?法吧（好累…）

這?還要將采樣頻率轉為15分鐘，不過最后博主還是克服難關啦：

當前時刻為：23:15異常的感應器有['感應器6','感應器12','感應器14','感應器20','感應器46']異常得分為140.87487118928732

當前時刻為：23:30異常的感應器有['感應器6','感應器12','感應器14','感應器20','感應器46']異常得分為140.51674777337848

當前時刻為：23:45異常的感應器有['感應器6','感應器12','感應器14','感應器20','感應器46']異常得分為140.48243399368258

當前時刻為：24:00異常的感應器有['感應器6','感應器12','感應器14','感應器20','感應器46']異常得分為140.52817104974736第四問

max(xi)maxmax(xi)max(si)公式如下：xix′i=100×

最后安全性得分為：

si=100?x′i

并再次標準化：sis′i=100×

最后得結果：

時間0days00:30:00安全得分:27.577947203229396時間0days01:00:00安全得分:25.836655421584666時間0days01:30:00安全得分:4.943973968727239時間0days02:00:00安全得分:14.130602320269498時間0days02:30:00安全得分:6.268057686488378時間0days03:00:00安全得分:3.8362635876100994時間0days03:30:00安全得分:0.0時間0days04:00:00安全得分:1.7766079380091397時間0days04:30:00安全得分:21.946804063556208時間0days05:00:00安全得分:13.284571837058333時間0days05:30:00安全得分:16.635549520871987時間0days06:00:00安全得分:35.44008528065004時間0days06:30:00安全得分:18.169844626792514時間0days07:00:00安全得分:7.220945101360403時間0days07:30:00安全得分:6.284431776175831時間0days08:00:00安全得分:23.198847074161215時間0days08:30:00安全得分:18.68607186443374時間0days09:00:00安全得分:12.531824761538987時間0days09:30:00安全得分:10.149786238806326時間0days10:00:00安全得分:10.202274482709573時間0days10:30:00安全得分:34.877504464067854時間0days11:00:00

時間0days11:00:00安全得分:27.78782832066013時間0days11:30:00安全得分:34.685365726453085時間0days12:00:00安全得分:33.379030884846614時間0days12:30:00安全得分:37.14737579778459時間0days13:00:00安全得分:29.62661767882518時間0days13:30:00安全得分:27.52891244607456時間0days14:00:00安全得分:23.757506965479664時間0days14:30:00安全得分:24.369572633810925時間0days15:00:00安全得分:15.526304239330543時間0days15:30:00安全得分:22.087026222737826時間0days16:00:00安全得分:16.80616936