2021華為杯數(shù)學建模 記錄

其實很早就想著記錄一下建模比賽的過程，但是大概是因為過程中似乎沒有令人記憶深刻的點，加之研究生第一個學期比較忙（課多+寫論文），所以給一直拖下去了。這學期論文第二稿提交后，也閑下來了一些，就趁著還有一些印象記錄下來。但畢竟也過了半年左右，難免有一些記不清楚的地方。
比賽官網(wǎng)鏈接: https://cpipc.acge.org.cn/.

組隊

三個人一個隊，研一的新生可以和其他學校的同學（比如本科或者高中的同學）一起組隊，研一以上的只能在本校組隊。我們當時就是宿舍的三個人，看通知有這個比賽而且學校報銷報名費就報上去了，當時沒有什么雄心壯志的。目標就是成功完賽交上論文退報名費。同時也沒有數(shù)學建模的基礎（師兄說這個比賽很簡單給了我們自信），也沒有提前看相關的內容，沒有分工。。。。
所以如果想沖刺好一些的名次的話，還是建議好好準備，聽別人說是上海落戶加分最簡單的比賽了。

選題

題目的壓縮包是很早就可以在網(wǎng)站上下載的，但是解壓需要密碼，到時間公布解壓的密碼，才可以看到里面的題目，這樣減輕了網(wǎng)站下載的壓力。
上午8點公布題，但我們宿舍睡到了11點半。。。起來看題：

• A題，華為專項題——“相關矩陣組的低復雜度計算和存儲建模”，要求是對矩陣儲存算法進行優(yōu)化，看了一下就pass了，屬于是題目都看不懂。
• B題，“空氣質量預報二次建模”，要求是已知一個實測數(shù)據(jù)和一次預報數(shù)據(jù)（一次預報是通過模型算出的，沒有考慮到實際數(shù)據(jù) ）利用實測數(shù)據(jù)對模型的一次預報數(shù)據(jù)進行建模優(yōu)化，得到更為準確的二次預報數(shù)據(jù)。看完題我和隊友說能做，因為我的本科畢設是做電離層預報的，兩者的思路很類似，做起來比較得心應手。
• C題，“帕金森病的腦深部電刺激治療建模研究”，要求是對腦部刺激后，神經(jīng)元的一系列反應和傳遞進行建模，看完題也沒太考慮，因為題目理解就感覺比較累。
• D題，“抗乳腺癌候選藥物的優(yōu)化建模”，主要是對化合物合成藥物的成分進行選擇，有一個同學選了這道題。
• E題，“信號干擾下的超寬帶（UWB）精確定位問題”，要求是對收集到的定位數(shù)據(jù)進行處理建模和應用場景的測試以及分類。這道題也比較想選，因為專業(yè)是和定位有關，相關的背景知識也很了解，而且看了一下題目給的數(shù)據(jù)，以及要求的分類模型，標簽也很清楚，我覺的用SVM去做是一定可以的，甚至效果會很好。
• F題，“航空公司機組優(yōu)化排班問題”，要求是對飛行員的排班優(yōu)化問題進行建模，因為題目太長加上前面有兩個候選題，就pass了。

看完題討論了一會就準備在B和E題中選了，然后出去吃飯的時候商量了一下，最后因為我說B題后面兩問我一個人就可以搞完，于是就決定做B題了。最后分工就是一個室友畫圖，一個室友做前兩問，我做后兩問。論文的話就各自寫各自做的部分。

問題重述

根據(jù)提供的監(jiān)測點長期空氣質量預報基礎數(shù)據(jù)，包括污染物濃度一次預報數(shù)據(jù)、氣象一次預報數(shù)據(jù)、氣象實測數(shù)據(jù)和污染物濃度實測數(shù)據(jù)進行分析建模，完成以下問題。

• 問題1：AQI計算

  計算監(jiān)測點A從2020年8月25日到8月28日每天實測的AQI和首要污染物。

• 問題2：氣象條件和AQI的相關性分析

  在污染物排放情況不變的條件下，某一地區(qū)的氣象條件有利于污染物擴散或沉降時，該地區(qū)的AQI會下降，反之會上升。根據(jù)對污染物濃度的影響程度，對氣象條件進行合理分類，并闡述各類氣象條件的特征。

• 問題3：單監(jiān)測點的空氣污染物二次預報模型

  建立一個同時適用于A、B、C三個監(jiān)測點的二次預報數(shù)學模型，用來預測未來三天6種常規(guī)污染物單日濃度值，要求二次預報模型預測結果中AQI預報值的最大相對誤差應盡量小，且首要污染物預測準確度盡量高。并預測監(jiān)測點A、B、C在2021年7月13日至7月15日6種常規(guī)污染物的單日濃度值，計算相應的AQI和首要污染物。

• 問題4：多監(jiān)測點的區(qū)域協(xié)同二次預報模型

  相鄰區(qū)域的污染物濃度往往具有一定的相關性，區(qū)域協(xié)同預報可能會提升空氣質量預報的準確度。監(jiān)測點A的臨近區(qū)域內存在監(jiān)測點A1、A2、A3，建立包含A、A1、A2、A3四個監(jiān)測點的協(xié)同預報模型，要求二次模型預測結果中AQI預報值的最大相對誤差應盡量小，且首要污染物預測準確度盡量高。使用該模型預測監(jiān)測點A、A1、A2、A3在2021年7月13日至7月15日6種常規(guī)污染物的單日濃度值，計算相應的AQI和首要污染物。并研究與問題3的模型相比，協(xié)同預報模型能否提升針對監(jiān)測點A的污染物濃度預報準確度？說明原因。

第一問解題

第一問很簡單，只要按照公式寫一些代碼計算出來就好了，注意一下臭氧是滑動平均值。

第二問解題

該問主要是做相關性分析，難度門檻也不高，所以該算的東西一定要算，然后根據(jù)算出來的相關性給到對應的解釋與分析，基本這一問就完成了。我覺的最重要的是把圖做好看，因為如果只給一堆數(shù)字，看上去就會很亂，用圖的話在論文里展示出來也會很豐富。
同時考慮到有些時間點缺失數(shù)據(jù)，后面兩問做的時候數(shù)據(jù)也需要補全，所以在做相關性分析之前，先把缺失數(shù)據(jù)的時間點給內插了出來，方法就是拉格朗日插值法。


計算過程：
首先使用監(jiān)測點A的每日污染物濃度數(shù)據(jù)計算出A點每日的AQI。之后根據(jù)每小時的氣象條件數(shù)據(jù)計算出日均的氣象條件參數(shù)，以二十四個整點時刻的均值代表當天的氣象條件參數(shù)水平。然后按全年十二個月將數(shù)據(jù)進行分類，計算每個月內AQI與各項氣象條件參數(shù)的相關系數(shù)，由此評估不同月份中氣象條件對AQI的不同影響，得到氣象條件與AQI變化的季節(jié)性特征。另外，計算了每種污染物與四種氣象條件參數(shù)共二十四組的相關系數(shù)，由此體現(xiàn)不同氣象條件對不同污染物濃度的影響的差異。
相關系數(shù)計算公式：

每個月的相關系數(shù)結果：

分析：

• 濕度：在一年中的絕大部分時間，空氣質量指數(shù)AQI與濕度呈負相關，相關系數(shù)為負，呈現(xiàn)出一定的季節(jié)性特征，在八、九月份相關系數(shù)絕對值達到最大，為-0.70，表明在這兩個月份濕度的上升會顯著降低AQI，在五月份和十二月份，AQI與濕度呈正相關，相關系數(shù)分別為0.10和0.20，在十一月份，AQI與濕度的相關系數(shù)為0，表明在該月內，AQI與濕度不相關。
• 壓強：AQI與其相關性呈現(xiàn)出一定的季節(jié)性特征，在春季兩者呈正相關，三月份和五月份兩者的相關系數(shù)達到最大，為0.30，在其它時間大都呈負相關，其中負相關系數(shù)絕對值最小的月份為一月和六月，相關系數(shù)為-0.20，負相關系數(shù)絕對值最大的月份為十二月，相關系數(shù)為-0.70，值得一提的是，從十月份到十二月份，AQI與壓強呈負相關，且相關系數(shù)絕對值從0.30一直增長到0.70，表明進入秋冬季后，AQI與壓強的相關性一直在增強，另外，在九月份，AQI與壓強呈現(xiàn)了正相關，這可能與及其它因素如監(jiān)測點A的地理位置和其它氣象因素等的影響有關。
• 風速：AQI與其相關系數(shù)在全年皆為負，在這種情況下，兩者的相關性仍呈現(xiàn)出一定的季節(jié)性特征，在春夏季，兩者的相關系數(shù)絕對值較小，絕對值最小出現(xiàn)在四月份和六月份，對應的相關系數(shù)為-0.20，絕對值最大也僅為0.5，出現(xiàn)在五月份，進入秋冬季后，AQI與風速的相關系數(shù)絕對值逐漸增大，相關性增強，最大出現(xiàn)在十二月份，最大絕對值為0.8。值得注意的是，根據(jù)[1]對成都市的污染物濃度變化研究，污染物的濃度均值隨著風速的增加而減小，AQI應該上升，但在本例的分析中，隨著風速的增加，AQI卻呈下降趨勢，即空氣質量降低，這可能是因為監(jiān)測點A的地理位置等因素導致風速增加利于污染物的形成或聚集，另外，分析時未加入風向對污染物濃度的影響，這可能也是導致AQI與風速呈負相關的原因之一。
• 溫度：在全年的大部分時間中，AQI與溫度呈正相關，且呈現(xiàn)出一定的季節(jié)性特征。在春季，兩者的相關性較低，三月份兩者的相關性最低，相關系數(shù)為0.0，表明兩者在三月份不相關，在春季兩者的相關系數(shù)絕對值最大也僅為0.30，最大出現(xiàn)在五月份，相關系數(shù)為-0.30。在夏季，兩者的相關性逐漸正向攀升，其相關系數(shù)從六月份的-0.10，到七月份的0.40，到八月份的0.70達到全年最高，此后一直維持在較高的水平，至一月份減小至0.20后，二月份又升高至0.50。 綜合考慮四種氣象條件參數(shù)對AQI的影響隨月份的變化可以發(fā)現(xiàn)，每種氣象條件參數(shù)與AQI的相關性都呈現(xiàn)出一定的季節(jié)性特征。
• 在春季，AQI與壓強呈正相關，與其它三種氣象條件參數(shù)呈負相關，但相關性最強的參數(shù)為風速，表明在春季，風速的變化對AQI的影響更顯著；
• 在夏季，AQI與溫度呈正相關，與其它三種氣象條件參數(shù)呈負相關，相關性最強的參數(shù)為溫度和濕度，另外兩種參數(shù)相較之下對AQI的影響更小；
• 在秋季，AQI與溫度呈正相關，與其它三種氣象條件參數(shù)呈負相關，而相關性最強的兩種參數(shù)依次為溫度和風速，表明在秋季，溫度和風速的變化對AQI的影響更顯著；
• 在冬季，AQI與溫度呈正相關，與其它三種氣象條件參數(shù)呈負相關，相關性最強的兩種參數(shù)依次為風速和溫度，表明在冬季，風速和溫度的變化對AQI的影響更顯著。
  為了呈現(xiàn)出氣象條件對每種污染物濃度各自的影響，本文對2019年4月16日 ~2021年7月13日二氧化硫（SO2）、二氧化氮（NO2）、粒徑小于10μm的顆粒物（PM10）、粒徑小于2.5μm的顆粒物（PM2.5）、臭氧（O3）、一氧化碳（CO）六種污染物濃度與四種氣象條件參數(shù)溫度(℃)、濕度(%)、氣壓(MBar)、風速(m/s)進行了線性回歸分析，結果如圖所示。
  
  另外，對于監(jiān)測站A的2020年7月23日~2021年7月13日的氣象一次預報數(shù)據(jù)，本文對其中所有的預報元素兩兩計算了相關系數(shù)，得到的相關性熱圖，如圖所示。
  
  可以看出，大部分預報元素之間相關性很弱，就本文所關心的六種污染物濃度而言，除臭氧之外的五種污染物，它們的每小時平均濃度與近地2米溫度，地表溫度，比濕，大氣壓，長波輻射相較別的元素有較強的相關性，這幾種因素涉及到的氣象條件參數(shù)有溫度，濕度，氣壓。另外，它們兩兩之間也有較強的相關性。對臭氧而言，它僅與邊界層高度、感熱通量、潛熱通量、短波輻射和地面太陽能輻射這四種元素有較強的相關性，而與其它元素的相關性很弱甚至沒有，這可能與的形成機制有關。

第三問解題

這一問我的思路就是用LSTM（長短期記憶模型）去做，因為之前對這一塊有比較多的了解和準備。解題的想法就是訓練一個模型，輸入是歷史數(shù)據(jù)序列（一次預報和實測數(shù)據(jù)），輸出是二次預報數(shù)據(jù)。而題目提供了兩年的數(shù)據(jù)，從數(shù)據(jù)量的角度看，完成一個模型的訓練也是足夠的。具體的步驟如下：

LSTM結構：

網(wǎng)絡模型：

第一個輸入?yún)?shù)為需要預報的空氣污染物的實測數(shù)據(jù)，對于每種污染物，分別輸入網(wǎng)絡處理。第二個參數(shù)為對應空氣污染物的一次預報數(shù)據(jù)。第三個參數(shù)輸入的則是上一節(jié)分析得到的與污染物相關性較強的實測數(shù)據(jù)，根據(jù)第二問中的分析，這里選擇為風速和濕度兩個氣象條件。由于臭氧O3是六種污染物中唯一的二次污染物，其并非來自污染源的直接排放，而在大氣中經(jīng)過一系列化學及光化學反應生成的。因此在對臭氧O3進行預報時，同時輸入了對應時間的NO2實測數(shù)據(jù)值。同時引入了日內時參數(shù) (hour in day, HD)，以使得網(wǎng)絡可以學習到日周期變化的特性，同時考慮參數(shù)的連續(xù)性，在處理時將 HD 參數(shù)分為 sin，cos 兩個分量 HDs,HDc 表示。

對于輸入的時間長度，由于輸入矩陣需要為標準矩陣，考慮了兩種設計：

由于需要預報未來三天64組數(shù)據(jù)（第一天8：00至第三天23：00，一小時一組數(shù)據(jù)），輸入數(shù)據(jù)中一次預報的數(shù)據(jù)時間即為未來這三天64小時對應的數(shù)據(jù)，而與實際觀測相關的數(shù)據(jù)則為預報時刻前的64組數(shù)據(jù)，這樣兩類數(shù)據(jù)可以組成一個二維矩陣輸入網(wǎng)絡。

將需要預報的未來64組數(shù)據(jù)分為N組，進行層層遞進的預報方式。如N=4時，每次預報64/4=16小時的數(shù)據(jù)，對于第一次預報，和1中類似，為構成二維矩陣，輸入數(shù)據(jù)中一次預報的數(shù)據(jù)時間即為未來16小時對應的數(shù)據(jù)，而與實際觀測相關的數(shù)據(jù)則為預報時刻前的16組數(shù)據(jù)。對于第二次預報，實際觀測相關的數(shù)據(jù)則替換為第一次預報的結果，以此類推，最后得到64小時的預報值。

在測試驗證過程中，我們發(fā)現(xiàn)兩種預報策略最后的精度相差不大，同時第二種策略訓練時間是第一種的數(shù)倍之多，權衡之下考慮到實際應用場景，我們最后選擇了第一種預報策略進行后續(xù)的研究和結果評估。
(關于這個點我還思考了挺久的，因為要讓輸入的數(shù)據(jù)是整齊的矩陣形式，而實測數(shù)據(jù)和一次預報的數(shù)據(jù)時間又不是對齊的，所以需要人為的把數(shù)據(jù)對齊，才能輸入網(wǎng)絡進行訓練)
后面訓練過程也沒什么好說的，這里放一個結果圖：

第四問解題

其實和第三問一樣，只是多了個需要考慮其他站的影響，對于我寫的LSTM模型來說，就是多加一個輸入的維度而已。但是我取了個很吊的名字——區(qū)域智能協(xié)同的空氣污染物二次預報模型（哈哈哈哈哈哈哈）

后面過程和第三問也基本一樣，放一個對比結果圖：

分析結果：
通過問題三和問題四的監(jiān)測站A在不同條件下進行的預報結果可知，兩種條件下的預報結果有一定的差異性和相似性。在兩種條件下的預報結果中，7月13日均無首要污染物， 并且在13日當天，各污染物濃度在兩種條件下相對變化率不大，因此AQI同樣差異不大但在7月14日與7月15日AQI差異較大，引起AQI差異的主要原因來自臭氧的變化，SO2、 PM10 、PM2.5、CO預報濃度發(fā)生的改變不明顯，NO2預報濃度差異較大，因此預報結果 的變化主要由O3與NO2引起。通過監(jiān)測站A、A1、A2、A3的聯(lián)測預報，A站O3與NO2在單 站預報的基礎上大大減小，而通過分析A1、A2、A3的數(shù)據(jù)可以看到，該三個站的O3與NO2 均比A站的O3與NO2要小，這是區(qū)域協(xié)同模型預報結果比單監(jiān)測站預報結果的根本原因。 通過加入測站周圍的測站，同時提供模型以地理信息相關因素，會使得監(jiān)測站的數(shù)據(jù)受到 周圍測站數(shù)據(jù)的修正，說明區(qū)域協(xié)同模型比單監(jiān)測站模型的預報結果更能體現(xiàn)該區(qū)域的未來污染物變化特性。

最后寫了一些模型的評價與推廣，就是硬吹，現(xiàn)在看尬的不行。

賽后總結

其實也沒啥好總結的，前面看題的時候太自信了，導致一共四天的比賽前面兩天大家都在上課開組會吃飯擺爛，結果到了最后兩天著急了，最后一天通宵搞到7點把論文提交了后吃了個飯回宿舍睡到了下午。也是因為沒什么經(jīng)驗，后面非常的趕，而我又要把兩問的代碼實現(xiàn)，把結果畫成圖，很多評估模型的部分我也沒能弄的比較詳細，一些后面覺的很棒的想法也沒能實現(xiàn)。但是也不能說后悔吧，怎樣的經(jīng)歷都是自己的選擇，自己的選擇就是最好的。以后爭取做的更好。
記憶最深刻的就是通宵那個晚上，我寫論文寫的有點累了，于是拿了瓶rio，去樓頂吹風，當時看著萬籟寂靜的珞珈山，靜靜的站著，感覺很奇怪，想起了在林芝的晚上，看著前面的南迦巴瓦峰，仿佛手中抓住了時間的流逝，山就在那里，可以遠遠的看著他，但更想去觸摸他。
我們都在朝著一個具體的方向向前奔走
12月份出結果，運氣很好，這份我覺的半完成度的論文拿了二等獎。

一些代碼
代碼可能有些亂，當時是ipynb格式寫的。其實很簡單，用tensorflow的keras庫直接搞的，搭建網(wǎng)絡很方便，訓練的代碼可能還沒有處理數(shù)據(jù)和畫圖寫的代碼多。

#!/usr/bin/env python # coding: utf-8from tensorflow import keras from tensorflow.keras.models import Sequential from tensorflow.keras import layers from tensorflow.keras.optimizers import RMSpropimport matplotlib.pyplot as plt import matplotlibimport numpy as np import csvfolder = r'C:\Users\yes\A站結果\\' real_csv = 'A點內插觀測值.csv' pre_csv = 'A點一次預報值.csv'with open(folder+real_csv) as f:data = csv.reader(f, delimiter=",")real_data = []for line in data:line = list(map(float,line))real_data.append((line)) real_data = np.array(real_data) print(real_data)with open(folder+pre_csv) as f:data = csv.reader(f, delimiter=",")pre_data = []for line in data:line = list(map(float,line))pre_data.append((line)) pre_data = np.array(pre_data) print(pre_data)pre_all_one = np.vstack((real_data[:,4],pre_data[:,4])) pre_all = np.hstack((pre_all_one.transpose(),real_data[:,10:12])) #針對O3特殊處理 pre_all = np.hstack((pre_all,real_data[:,1].reshape(len(pre_all),1)))data_mean = pre_all.mean(axis=0) data_std = pre_all.std(axis=0)deal_data = (pre_all - data_mean) / data_stddef generator(data, lookback, delay, min_index, max_index, target_num=0,shuffle=False, batch_size=64):if max_index is None:max_index = len(data) - delay - 1i = min_index + lookbackwhile 1:if shuffle:rows = np.random.randint(min_index + lookback, max_index, size=batch_size)else:if i + batch_size >= max_index:i = min_index + lookbackrows = np.arange(i, min(i + batch_size, max_index))i += len(rows)samples = np.zeros((len(rows),lookback,data.shape[-1]))targets = np.zeros((len(rows),delay))for j, row in enumerate(rows):indices_sam = range(rows[j] - lookback, rows[j])indices_tar = range(rows[j] + 1, rows[j] + 1 + delay)samples[j] = data[indices_sam]targets[j] = data[indices_tar][:,0]yield samples, targetslookback = 64 delay = 64 batch_size = 16train_gen = generator(deal_data,lookback=lookback,delay=delay,min_index=0,max_index=None,shuffle=True,batch_size=batch_size) val_gen = generator(deal_data,lookback=lookback,delay=delay,min_index=8001,max_index=None,batch_size=batch_size) test_gen = generator(deal_data,lookback=lookback,delay=delay,min_index=8001,max_index=None,batch_size=batch_size)# This is how many steps to draw from `val_gen` # in order to see the whole validation set: val_steps = (len(deal_data) - 8001 - lookback) // batch_size# This is how many steps to draw from `test_gen` # in order to see the whole test set: test_steps = (len(deal_data) - 8001 - lookback) // batch_sizemodel = Sequential() model.add(layers.LSTM(64,#return_sequences = True,input_shape = (None,deal_data.shape[-1])))#input_shape = (None,1))) #model.add(layers.LSTM(32)) model.add(layers.Dense(64))model.compile(optimizer=RMSprop(), loss='mae') history = model.fit(train_gen,steps_per_epoch=200,epochs=30,validation_data=val_gen,validation_steps=val_steps,verbose = 1)import matplotlib.pyplot as pltloss = history.history['loss'] val_loss = history.history['val_loss']epochs = range(len(loss))plt.figure() plt.subplots(figsize=(20,10)) plt.plot(epochs, loss, 'ro', label='Training loss') plt.plot(epochs, val_loss, 'b', label='Validation loss') plt.title('Training and validation loss') plt.legend(fontsize=24) plt.show()test_data = deal_data[-64:,].reshape(1,64,5) results=model.predict(test_data, batch_size=1) ypx = results * data_std[0] + data_mean[0]#設置標簽格式 matplotlib.rcParams['xtick.labelsize'] = 18 matplotlib.rcParams['ytick.labelsize'] = 18 # X、Y軸刻度標簽字體大小 matplotlib.rcParams['axes.titlesize'] = 24 matplotlib.rcParams['axes.labelsize'] = 24 fig, ax = plt.subplots(figsize=(10,6),) plt.ylim(0,200) plt.plot(pre_all[-64:,1] , '^', color='r', label='First_predict', markersize=8) plt.plot(ypx.reshape(64), 'o', color='b', label='Twice_predict', markersize=8) plt.legend(fontsize=20) plt.xlabel('Epoch(1h)') plt.ylabel('O3(μg/m3)') plt.savefig('O3.jpg', dpi=300)#保存圖片

總結

以上是生活随笔為你收集整理的“华为杯”中国研究生数学建模竞赛 记录的全部內容，希望文章能夠幫你解決所遇到的問題。

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