當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

2021 年五一数学建模比赛 C 题

發布時間：2023/12/9 编程问答 29 豆豆

生活随笔收集整理的這篇文章主要介紹了 2021 年五一数学建模比赛 C 题小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

文章目錄

思路
第一問
- 風險分析
- - 獨立性、偶發性誤差（非風險也）
  - 風險
  - 灰色地帶
  - - 規律性波動
    - 沖擊性誤差
- 沖擊性風險判定與算法 A
- 獨立性和偶發性誤差
- - 正態分布檢驗—— Anderson–Darling 檢驗（原理部分別看）
  - 均勻分布檢驗——Kolmogorov–Smirnov 檢驗
  - 非正常波動系數
- 自相關性系數和自相關檢驗
- - 自相關系數
  - 相關性檢驗
  - 注意點
  - - 穩定性檢驗——Augmented Dickey–Fuller Test
- 風險計算總結
- 結果展示
第二問
第三問
第四問
代碼與提問
- 編程收錄

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

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

賽題官網：http://51mcm.cumt.edu.cn/

思路

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

如何量化呢？這里用統計的方法，對每一個感知器的所有時序數據（共 5519 個），從沖擊性風險、獨立性和偶發性風險、以及自相關風險三個方面，來評價感知器的風險值。

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

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

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

第一問

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

正常波動：外界溫度或者產量變化的波動，傳感器誤報，規律性、獨立性、偶發性

非正常波動：生產過程中的不穩定因素，持續性、聯動性

注意，正常波動中，有一個規律性，所以，正常波動是白噪聲嗎？

我們把沒有持續性和聯動性的參數，就叫非正常波動。

風險分析

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

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

所以，感應器誤識別，和因溫度等因素引起的獨立性和偶發性波動，都不是風險。

獨立性、偶發性誤差（非風險也）

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

風險

根據題目，所謂風險應該是有持續性的、聯動性的，如下所示：

灰色地帶

規律性波動

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

從題目可知，有些誤差是具有規律性的，但規律性，又意味著聯動性和持續性。所以，對傳感器 2 來說，我們可以說是因為溫度等因素，也可以歸咎于生產過程。

但對本題而言，我們將類似傳感器 2 的規律性波動，視為風險。

沖擊性誤差

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

沖擊性風險判定與算法 A

算法 A 來自 GB/T 6379.5，應用此算法計算得到數據平均值和標準差的穩健值，但也可以據此統計沖擊性風險的數據含量。

算法 A 的計算過程如下所示：
首先是計算初始值：
$x?=med?xis?=1.483×med?∣xi?x?∣\begin{array}{c} x^{*}=\operatorname{med} x_{i} \\ s^{*}=1.483 \times \operatorname{med}\left|x_{i}-x^{*}\right| \end{array}$
對每個 $x_i$ ，有：
$xi?={x??δ,若?xi<x??δx?+δ,若?xi>x?+δxi,其他?x_{i}^{*}=\left\{\begin{array}{cc} x^{*}-\delta, & \text { 若 } x_{i}<x^{*}-\delta \\ x^{*}+\delta, & \text { 若 } x_{i}>x^{*}+\delta \\ x_{i}, & \text { 其他 } \end{array}\right.$
其中： $δ=1.5s?\delta=1.5 s^{*}$ ，再次計算：
$x?=∑xi?/ps?=1.134∑(xi??x?)2/(p?1)\begin{array}{c} x^{*}=\sum x_{\mathrm{i}}^{*} / p \\ s^{*}=1.134 \sqrt{\sum\left(x_{\mathrm{i}}^{*}-x^{*}\right)^{2} /(p-1)} \end{array}$

重復：
$xi?={x??δ,若?xi<x??δx?+δ,若?xi>x?+δxi,其他?x_{i}^{*}=\left\{\begin{array}{cc} x^{*}-\delta, & \text { 若 } x_{i}<x^{*}-\delta \\ x^{*}+\delta, & \text { 若 } x_{i}>x^{*}+\delta \\ x_{i}, & \text { 其他 } \end{array}\right.$

直到 $s^*$ 的第三位有效數字和 $x^*$ 的對應數字在連續兩次迭代中不變。

因此，對感應器的數據，我們通過算法 A 后，得到最終的 $x_i^*$ ，記取值為 $x?+δx^*+\delta$ 或 $x??δx^*-\delta$ 的 $x_i^*$ 的數量為 $m$ ，于是沖擊性風險為 $m / n$ 。

我們也可以搞復雜一點，因為沖擊性持續時間越長，意味著風險越大。因此，若取值為 $x?+δx^*+\delta$ 或 $x??δx^*-\delta$ 的數連續出現，我們可以給其添加權重。

如出現以及，記為 1 分，連續出現 2 次，記為 $1 + 2$ ，連續出現 3 次，記為 $1 + 2 + 3$ …，若間隔出現 2 次，則記為 $1 + 1$ ，以此類推，最終求和得到 $m′m^\prime$ 。

于是，沖擊性風險即為：
$risk1=m′nrisk_1 =\frac {m^\prime}{n}$

獨立性和偶發性誤差

獨立性和偶發性，讓一個原本平滑的時間序列，成為了一組白噪聲。

理想狀態下，傳感器的數值應為：

加入白噪聲后，應為：

或：

因此，若一個時間序列，在打亂其時序后，滿足正態分布，則證明時間序列有白噪聲。又或者，時間序列滿足均勻分布，如上圖傳感器 8 ，則亦可以證明時間序列的風險低。

正態分布檢驗—— Anderson–Darling 檢驗（原理部分別看）

實踐表明，AD 檢驗檢驗正態分布，似乎比 KS 檢驗還要差！，所以本文用 KS 檢驗了

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

另外，比起 Shapiro-Wilk 檢驗，也是針對正態分布的特定的檢驗，但在樣本容量大于 5000 時，則準確度不高。具體見：Shapiro Wilk 檢驗維基百科：https://en.wikipedia.org/wiki/Shapiro%E2%80%93Wilk_test

AD 檢驗的原假設為：樣本服從正態分布 $N(μ,σ2)N(\mu, \sigma^2)$ ，若 $μ,σ\mu,\sigma$ 未知，則可以用樣本的估計量代替。

AD 檢驗的檢驗統計量為：
$A2=?n?1n∑i=1n(2i?1)(ln?Φ(Yi)+ln?(1?Φ(Yn?1?i)))A^{2}=-n-\frac{1}{n} \sum_{i=1}^{n}(2 i-1)\left(\ln \Phi\left(Y_{i}\right)+\ln \left(1-\Phi\left(Y_{n-1-i}\right)\right)\right)$
其中 $Y_i$ 為：
$Yi=Xi?μ^σ^Y_{i}=\frac{X_{i}-\hat{\mu}}{\hat{\sigma}}$
均值和方差：
$μ^={μ,if?the?mean?is?known.?Xˉ,=1n∑i=1nXiotherwise.?σ^2={σ2,if?the?variance?is?known.?1n∑i=1n(Xi?μ)2,if?the?variance?is?not?known,?but?the?mean?is.?1n?1∑i=1n(Xi?Xˉ)2,otherwise.?\begin{array}{l} \hat{\mu}=\left\{\begin{array}{ll} \mu, & \text { if the mean is known. } \\ \bar{X},=\frac{1}{n} \sum_{i=1}^{n} X_{i} & \text { otherwise. } \end{array}\right. \\ \hat{\sigma}^{2}=\left\{\begin{array}{ll} \sigma^{2}, & \text { if the variance is known. } \\ \frac{1}{n} \sum_{i=1}^{n}\left(X_{i}-\mu\right)^{2}, & \text { if the variance is not known, but the mean is. } \\ \frac{1}{n-1} \sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)^{2}, & \text { otherwise. } \end{array}\right. \end{array}$
而 $X_i$ 為次序統計量，即將原始數據根據升序的方式進行排序。檢驗統計量 $A^2$ 的臨界值如表所示：

若 $A^2$ 大于臨界值，則拒絕原假設，即認為樣本不服從正態分布。

當然， $A^2$ 在樣本接近 $∞\infin$ 的時候，如本例的 5000 個樣本，則 $A^2$ 服從：
$Pr?(A2<z)=2πz∑i=0∞(?12j)(4i+1)e?(4i+1)2π2/(8z)∫0∞ez8(1+w2)?w2(4i+1)2π2/(8z)dw\operatorname{Pr}\left(A^2<z\right)=\frac{\sqrt{2 \pi}}{z} \sum_{i=0}^{\infty}\left(\begin{array}{c} -\frac{1}{2} \\ j \end{array}\right)(4 i+1) e^{-(4i+1)^{2} \pi^{2} /(8 z)} \int_{0}^{\infty} e^{\frac{z}{8\left(1+w^{2}\right)}-w^{2}(4 i+1)^{2} \pi^{2} /(8 z)} d w$

當然，也可以根據近似公式計算 p-value，如下：
$p?value={exp?(1.2937?5.709AD?.0186AD2),AD≥0.6exp?(0.9177?4.279AD?1.38AD2)AD≥0.34exp?(?8.318+42.796AD?59.938AD2)AD≥0.2exp?(?13.436+101.14AD?223.73AD2)otherwise\begin{array}{l} p-value =\left\{\begin{array}{ll} \exp(1.2937 - 5.709AD - .0186AD^2), & AD \geq 0.6 \\ \exp(0.9177 - 4.279AD - 1.38AD^2) & AD \geq 0.34 \\ \exp(-8.318 + 42.796AD - 59.938AD^2) & AD \geq 0.2\\ \exp(-13.436 + 101.14AD - 223.73AD^2) & \text{otherwise} \\ \end{array} \right. \end{array}$
其中 AD 為 $A^2$ 的調整：
$A^2 \times (1 + (.75/n) + 2.25/(n^2))$
因此我們可以計算出相應的 p-value。若 p-value 越接近于 1，意味著樣本越有可能是正態分布，也即風險越小。

參考文獻：
A^2 臨界值表
[A^2 的分布函數]：Evaluating the Anderson-Darling Distribution 作者：George Marsaglia。
AD 檢驗維基百科
（上述文獻我放到代碼文件里了）
AD p-value 近似計算

均勻分布檢驗——Kolmogorov–Smirnov 檢驗

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

但我們可以用 Kolmogorov 檢驗，來檢驗樣本的總體是否滿足均勻分布。定義檢驗統計量如下：
$Dn=sup?x∣Fn(x)?F(x)∣D_{n}=\sup _{x}\left|F_{n}(x)-F(x)\right|$
其中 $F_n(x)$ 為經驗分布， $F (x)$ 為實際分布，如下：
$Fn(x)=1n∑i=1nI[?∞,x](Xi)F_{n}(x)=\frac{1}{n} \sum_{i=1}^{n} I_{[-\infty, x]}\left(X_{i}\right)$
$I[?∞,x](Xi)I_{[-\infty, x]}\left(X_{i}\right)$ 為指示函數，若有樣本落在范圍內，則為 1。經驗函數根據樣本分布，模擬了總體的分布函數（CDF）。 $sup_{x}$ 是上確界之意。

Kolmogorov 檢驗的原假設為：樣本 $x_i$ 來源于總體分布 $F (x)$ 。

根據 Glivenko–Cantelli 原理，若樣本 $x_i$ 來源于總體分布 $F (x)$ ，則 $D_n$ 將收斂于 0。或：
$nDn?π→∞sup?t∣B(F(t))∣\sqrt{n} D_{n} \stackrel{\pi \rightarrow \infty}{\longrightarrow} \sup _{t}|B(F(t))|$
其中， $B (t)$ 服從柯爾莫格羅夫分布，其分布函數如下：
$Pr?(K≤x)=1?2∑k=1∞(?1)k?1e?2k2x2=2πx∑k=1∞e?(2k?1)2π2/(8x2)\operatorname{Pr}(K \leq x)=1-2 \sum_{k=1}^{\infty}(-1)^{k-1} e^{-2 k^{2} x^{2}}=\frac{\sqrt{2 \pi}}{x} \sum_{k=1}^{\infty} e^{-(2 k-1)^{2} \pi^{2} /\left(8 x^{2}\right)}$

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

Kolmogorov 檢驗：https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test

非正常波動系數

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

于是可以評價時間序列的非白噪聲風險系數為：
$risk2′=1?max?(psw′,pks′)risk_2^\prime = 1 - \max(p_{sw}^\prime, p_{ks}^\prime)$

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

由于風險信號具有聯動性和持續性，所以，我們可以用自相關系數來評判，我們的傳感器數據，是否也具備聯動性和持續性。

自相關系數

自相關系數旨在計算當前時刻，和前 k 個時刻的 Pearson 相關系數其計算公式如下：
$\frac{\sum_{i=k+1}^n (x_t - \bar{x}_t) (x_{t-k} - \bar{x}_{t-k})} {\sqrt{\sum_{i=k+1}^n (x_t - \bar{x}_t)^2 (x_{t-k}-\bar{x}_{t-k})^2}}$

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

參考：https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.autocorr.html 和 Pearson 相關系數的定義寫出。

注意點

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

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

穩定性檢驗——Augmented Dickey–Fuller Test

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

為此，本文采用：Augmented Dickey–Fuller 檢驗

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

因此，若拒絕了 ADF 降壓，則讓 $r = 0$

風險計算總結

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

算法 A ：沖擊性風險

正態、均勻分布檢驗：量化誤差是否屬于風險

自相關系數：持續性和聯動性風險

結果展示

測得幾個風險最大的感應器，以及其圖片如下所示：

第二問

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

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

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

當前時刻為：0 days 03:00:00 異常的感應器有： ['感應器2', '感應器5', '感應器6', '感應器7', '感應器11'] 異常得分為： 156.45667125866223 當前時刻為：0 days 16:30:00 異常的感應器有： ['感應器6', '感應器7', '感應器11', '感應器12', '感應器13'] 異常得分為： 156.29759895209494 當前時刻為：0 days 03:30:00 異常的感應器有： ['感應器2', '感應器5', '感應器6', '感應器7', '感應器11'] 異常得分為： 156.1732779754605 當前時刻為：0 days 02:30:00 異常的感應器有： ['感應器2', '感應器5', '感應器6', '感應器7', '感應器11'] 異常得分為： 155.8447347080813 當前時刻為：0 days 01:00:00 異常的感應器有： ['感應器57', '感應器5', '感應器6', '感應器7', '感應器11'] 異常得分為： 155.66803974065186

第三問

我們可以考慮兩種方法：

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

用機器學習的方法，先擬合一個機器學習模型，再用來預測

這里就考慮第一種方法吧（好累…）

這里還要將采樣頻率轉為 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

第四問

至于第四問，其實就比較容易了。根據第二問求出的異常總得分，將其標準化為 100 即可。設 30 分鐘為間隔的異常總得分為 $x_i$ ，標準化公式如下：
$xi′=100×ximax?(xi)x_i^\prime = 100\times \frac{x_i}{\max{(x_i)}}$
最后安全性得分為：
$si=100?xi′s_i = 100-x_i^\prime$
并再次標準化：
$si′=100×simax?(si)s_i^\prime = 100\times \frac{s_i}{\max{(s_i)}}$
最后得結果：

時間： 0 days 00:30:00 安全得分: 27.577947203229396 時間： 0 days 01:00:00 安全得分: 25.836655421584666 時間： 0 days 01:30:00 安全得分: 4.943973968727239 時間： 0 days 02:00:00 安全得分: 14.130602320269498 時間： 0 days 02:30:00 安全得分: 6.268057686488378 時間： 0 days 03:00:00 安全得分: 3.8362635876100994 時間： 0 days 03:30:00 安全得分: 0.0 時間： 0 days 04:00:00 安全得分: 1.7766079380091397 時間： 0 days 04:30:00 安全得分: 21.946804063556208 時間： 0 days 05:00:00 安全得分: 13.284571837058333 時間： 0 days 05:30:00 安全得分: 16.635549520871987 時間： 0 days 06:00:00 安全得分: 35.44008528065004 時間： 0 days 06:30:00 安全得分: 18.169844626792514 時間： 0 days 07:00:00 安全得分: 7.220945101360403 時間： 0 days 07:30:00 安全得分: 6.284431776175831 時間： 0 days 08:00:00 安全得分: 23.198847074161215 時間： 0 days 08:30:00 安全得分: 18.68607186443374 時間： 0 days 09:00:00 安全得分: 12.531824761538987 時間： 0 days 09:30:00 安全得分: 10.149786238806326 時間： 0 days 10:00:00 安全得分: 10.202274482709573 時間： 0 days 10:30:00 安全得分: 34.877504464067854 時間： 0 days 11:00:00 安全得分: 27.78782832066013 時間： 0 days 11:30:00 安全得分: 34.685365726453085 時間： 0 days 12:00:00 安全得分: 33.379030884846614 時間： 0 days 12:30:00 安全得分: 37.14737579778459 時間： 0 days 13:00:00 安全得分: 29.62661767882518 時間： 0 days 13:30:00 安全得分: 27.52891244607456 時間： 0 days 14:00:00 安全得分: 23.757506965479664 時間： 0 days 14:30:00 安全得分: 24.369572633810925 時間： 0 days 15:00:00 安全得分: 15.526304239330543 時間： 0 days 15:30:00 安全得分: 22.087026222737826 時間： 0 days 16:00:00 安全得分: 16.806169363179187 時間： 0 days 16:30:00 安全得分: 27.19324868501274 時間： 0 days 17:00:00 安全得分: 0.9972329596947986 時間： 0 days 17:30:00 安全得分: 9.87443946696844 時間： 0 days 18:00:00 安全得分: 7.979879130013369 時間： 0 days 18:30:00 安全得分: 23.5299343143491 時間： 0 days 19:00:00 安全得分: 14.035391935319838 時間： 0 days 19:30:00 安全得分: 31.66537146146771 時間： 0 days 20:00:00 安全得分: 34.64135831503899 時間： 0 days 20:30:00 安全得分: 7.402710838145929 時間： 0 days 21:00:00 安全得分: 15.454310630953227 時間： 0 days 21:30:00 安全得分: 18.77425423612869 時間： 0 days 22:00:00 安全得分: 24.84119865998032 時間： 0 days 22:30:00 安全得分: 20.35897364747012 時間： 0 days 23:00:00 安全得分: 6.220522144385195 時間： 0 days 23:30:00 安全得分: 98.80570148706612 時間： 1 days 00:00:00 安全得分: 100.0

代碼與提問

若需要代碼，請點贊，關注、私信、說明題目和年份

如果有其他問題，請到評論區留言，私信提問，概不回答。也在此鼓勵大家獨立思考。

本人不會回訪，不互關，不互吹，以及謝絕諸如此類事

編程收錄

編程之精髓，固在復用。每做題之際，偶有靈光一現，或得謄錄網絡妙筆，觀之不忍釋手，是故獨辟一卷，畢錄于茲，以供覽閱。時在拋磚引玉，虛左待教。或偶逢旮旯罅隙，探索漫漫，復而柳暗花明，何其喜洋洋矣，然忖后來者亦常陷類似之困囿者何其蕓蕓，榫卯之機，舉手之勞，此間不可不錄也。

算法 A 代碼

總結

以上是生活随笔為你收集整理的2021 年五一数学建模比赛 C 题的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。

上一篇：利用 VMware vRealize -
下一篇： 3732: Network