研究生数学建模比赛感触
建模 数学 比赛 研究生
2023-06-13 09:13:16 时间
研究生数学建模比赛感触
一.题目
本次研究生建模比赛，我深有感触，我们队友一起选的B题，题目是关于空气质量的二次建模问题，此题明显是考察数据分析的，我做两问，以及第三问借鉴一下别人的想法。总体感觉难度有点大，因为是第一次参加，对模型不够熟练，只做一部分，不过也收获了不少，比如数据预处理，表的合并，写入文件，以及缺失值处理，还有统计方法，除此外还用到模型。主要是python的pandas,numpy以及matplotlib库的使用，还有聚类（KMeans），逻辑回归，支持向量机的应用.
题目
二.建模过程以及代码、图片、数据等
问题一 这个应该就是一个公式计算问题
代码如下（使用 jupyter notebook）:
#!/usr/bin/env python
# coding: utf-8
# In[1]:
import pandas as pd
import numpy as np
import math
# In[2]:
data = pd.read_excel("附件1 监测点A空气质量预报基础数据.xlsx",sheet_name=2)
data1 =data[497:501]
data1
# In[3]:
#得到各个时间段的浓度值
so2 = data1.iloc[:,2].tolist()
no2 = data1.iloc[:,3].tolist()
pm10 = data1.iloc[:,4].tolist()
pm25 = data1.iloc[:,5].tolist()
o3 = data1.iloc[:,6].tolist()
co = data1.iloc[:,7].tolist()
# In[4]:
#co的AOI计算函数
def co_js(co):
if co >= 0 and co <2:
BPhi=2
BPlo=0
IAQIhi=50
IAQIlo=0
Cp=co
elif co >= 2 and co <4:
BPhi=4
BPlo=2
IAQIhi=100
IAQIlo=50
Cp=co
elif co >= 4 and co < 14:
BPhi=14
BPlo=4
IAQIhi=100
IAQIlo=150
Cp=co
elif co >= 14 and co < 24:
BPhi=24
BPlo=14
IAQIhi=200
IAQIlo=150
Cp=co
elif co >= 24 and co < 36:
BPhi=36
BPlo=24
IAQIhi=300
IAQIlo=200
Cp=co
elif co >=36 and co < 48:
BPhi=48
BPlo=36
IAQIhi=400
IAQIlo=300
Cp=co
elif co >= 48 and co < 60:
BPhi=60
BPlo=48
IAQIhi=500
IAQIlo=400
Cp=co
else:
IAQIp=0
IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo)+IAQIlo))
return IAQIp
# In[5]:
#co的AOI计算结果
co_IAQI=[]
for i in co:
co_IAQI.append(co_js(i))
co_IAQI
# In[6]:
def so2_js(so2):
if so2 >= 0 and so2 < 50:
BPhi=50
BPlo=0
IAQIhi=50
IAQIlo=0
Cp=so2
elif so2 >= 50 and so2 < 150:
BPhi=150
BPlo=50
IAQIhi=100
IAQIlo=50
Cp=so2
elif so2 >= 150 and so2 < 475:
BPhi=475
BPlo=150
IAQIhi=150
IAQIlo=100
Cp=so2
elif so2 >= 475 and so2 < 800:
BPhi=800
BPlo=475
IAQIhi=200
IAQIlo=150
Cp=so2
elif so2 >= 800 and so2 < 1600:
BPhi=1600
BPlo=800
IAQIhi=300
IAQIlo=200
Cp=so2
elif so2 >= 1600 and so2 < 2100:
BPhi=2100
BPlo=1600
IAQIhi=400
IAQIlo=300
Cp=so2
elif co >= 2100 and co < 2620:
BPhi=2620
BPlo=2100
IAQIhi=500
IAQIlo=400
Cp=so2
else:
IAQIp=0
IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo)+IAQIlo))
return IAQIp
# In[7]:
so2_IAQI=[]
for i in so2:
so2_IAQI.append(so2_js(i))
so2_IAQI
# In[8]:
def no2_js(no2):
if no2 >= 0 and no2 < 40:
BPhi=40
BPlo=0
IAQIhi=50
IAQIlo=0
Cp=no2
elif no2 >= 40 and no2 < 80:
BPhi=80
BPlo=40
IAQIhi=100
IAQIlo=50
Cp=no2
elif no2 >= 80 and no2 < 180:
BPhi=180
BPlo=80
IAQIhi=150
IAQIlo=100
Cp=no2
elif no2 >= 180 and so2 < 280:
BPhi=280
BPlo=180
IAQIhi=200
IAQIlo=150
Cp=no2
elif no2 >= 280 and no2 < 565:
BPhi=565
BPlo=280
IAQIhi=300
IAQIlo=200
Cp=no2
elif no2 >= 565 and no2 < 750:
BPhi=750
BPlo=565
IAQIhi=400
IAQIlo=300
Cp=no2
elif co >= 750 and co < 940:
BPhi=940
BPlo=750
IAQIhi=500
IAQIlo=400
Cp=no2
else:
IAQIp=0
IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo)+IAQIlo))
return IAQIp
# In[9]:
no2_IAQI=[]
for i in so2:
no2_IAQI.append(no2_js(i))
no2_IAQI
# In[10]:
def pm10_js(pm10):
if pm10>= 0 and pm10 < 50:
BPhi=50
BPlo=0
IAQIhi=50
IAQIlo=0
Cp=pm10
elif pm10 >= 50 and pm10 < 150:
BPhi=150
BPlo=50
IAQIhi=100
IAQIlo=50
Cp=pm10
elif pm10 >= 150 and pm10 < 250:
BPhi=250
BPlo=150
IAQIhi=150
IAQIlo=100
Cp=pm10
elif pm10 >= 250 and pm10 < 350:
BPhi=350
BPlo=250
IAQIhi=200
IAQIlo=150
Cp=pm10
elif pm10 >= 350 and pm10 < 420:
BPhi=420
BPlo=360
IAQIhi=300
IAQIlo=200
Cp=pm10
elif pm10 >= 420 and pm10< 500:
BPhi=500
BPlo=420
IAQIhi=400
IAQIlo=300
Cp=pm25
elif pm10 >= 500 and pm10 < 600:
BPhi=600
BPlo=500
IAQIhi=500
IAQIlo=400
Cp=pm10
else:
IAQIp=0
IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo)+IAQIlo))
return IAQIp
# In[11]:
pm10_IAQI=[]
for i in pm10:
pm10_IAQI.append(pm10_js(i))
pm10_IAQI
# In[12]:
def pm25_js(pm25):
if pm25 >= 0 and pm25 < 35:
BPhi=35
BPlo=0
IAQIhi=50
IAQIlo=0
Cp=pm25
elif pm25 >= 35 and pm25 < 75:
BPhi=75
BPlo=35
IAQIhi=100
IAQIlo=50
Cp=pm25
elif pm25 >= 75 and pm25 < 115:
BPhi=115
BPlo=75
IAQIhi=150
IAQIlo=100
Cp=pm25
elif pm25 >= 115 and pm25 < 150:
BPhi=350
BPlo=250
IAQIhi=200
IAQIlo=150
Cp=pm25
elif no2 >= 150 and no2 < 250:
BPhi=250
BPlo=150
IAQIhi=300
IAQIlo=200
Cp=pm25
elif no2 >= 250 and no2 < 350:
BPhi=350
BPlo=250
IAQIhi=400
IAQIlo=300
Cp=pm25
elif co >= 350 and co < 500:
BPhi=500
BPlo=350
IAQIhi=500
IAQIlo=400
Cp=pm25
else:
IAQIp=0
IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo)+IAQIlo))
return IAQIp
# In[13]:
pm25_IAQI=[]
for i in pm25:
pm25_IAQI.append(pm25_js(i))
pm25_IAQI
# In[14]:
def o3_js(o3):
if o3 >= 0 and o3 < 100:
BPhi=100
BPlo=0
IAQIhi=50
IAQIlo=0
Cp=o3
elif o3 >= 100 and o3 < 160:
BPhi=160
BPlo=100
IAQIhi=100
IAQIlo=50
Cp=o3
elif o3 >= 160 and o3 < 215:
BPhi=115
BPlo=75
IAQIhi=150
IAQIlo=100
Cp=o3
elif o3 >= 215 and o3 < 265:
BPhi=350
BPlo=250
IAQIhi=200
IAQIlo=150
Cp=o3
elif no2 >= 265 and no2 < 800:
BPhi=800
BPlo=265
IAQIhi=300
IAQIlo=200
Cp=o3
else:
IAQIp=0
IAQIp=int(math.ceil((IAQIhi-IAQIlo)*(Cp-BPlo)/(BPhi-BPlo)+IAQIlo))
return IAQIp
# In[15]:
o3_IAQI=[]
for i in o3:
o3_IAQI.append(o3_js(i))
o3_IAQI
# In[16]:
#整合数据
AQIA = np.transpose(np.array([co_IAQI,so2_IAQI,no2_IAQI,o3_IAQI,pm10_IAQI,pm25_IAQI]))
AQIA
# In[17]:
#统计浓度最大值
AQIA_end = AQIA.max(axis=1)
AQIA_end
# In[18]:
AQIA_index = np.argmax(AQIA,axis=1)
AQIA_index
# In[19]:
#首要污染物名称
AQIA_columns = []
columns = ["co","so2","no2","o3","pm10","pm25"]
for i in AQIA_index:
AQIA_columns .append(columns[i])
AQIA_columns
# In[20]:
#首要污染物
wrw = []
for i in range(4):
if AQIA_end[i] <= 50:
wrw.append("无")
else:
wrw.append(AQIA_columns[i])
wrw
# In[21]:
sj = data1.iloc[:,0]
sj = np.array(sj,dtype="datetime64[D]")
dd = data1.iloc[:,1]
# In[22]:
sj
# In[23]:
#写入文件
columns =["监测日期","地点","AQI","首要污染物"]
DATA = pd.DataFrame(data=np.array([sj,dd,AQIA_end,wrw]).T,columns=columns)
DATA.to_excel("AQI.xlsx")
运行结果如图：
监测日期 | 地点 | AQI | 首要污染物 | |
|---|---|---|---|---|
0 | 2020-08-25 | 监测点A | 60 | |
1 | 2020-08-26 | 监测点A | 46 | 无 |
2 | 2020-08-27 | 监测点A | 218 | |
3 | 2020-08-28 | 监测点A | 258 |
问题二：做的分类，不知道是否正确，分成六类
#!/usr/bin/env python
# coding: utf-8
# In[1]:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
# In[2]:
data = pd.read_excel("附件1 监测点A空气质量预报基础数据.xlsx",sheet_name=1)
data = data.replace("—",np.nan)
data1 = data.dropna(axis=0,how="any")
# In[3]:
data1
# In[4]:
so2 = np.array(data1.iloc[:,2])
no2 = np.array(data1.iloc[:,3])
pm10 = np.array(data1.iloc[:,4])
pm25 = np.array(data1.iloc[:,5])
o3 = np.array(data1.iloc[:,6])
co = np.array(data1.iloc[:,7])
wd = np.array(data1.iloc[:,8])
sd = np.array(data1.iloc[:,9])
qy = np.array(data1.iloc[:,10])
fs = np.array(data1.iloc[:,11])
fx = np.array(data1.iloc[:,12])
# In[5]:
y = [so2,no2,pm10,pm25,o3,co]
colors =["red","black","orange","blue","green","purple"]
list = ["so2","no2","pm10","pm25","o3","co"]
# In[6]:
for i in range(6):
plt.figure(figsize=(20,8))
plt.subplot(3,2,i+1)
plt.scatter(wd,y[i],c=colors[i])
plt.xlabel("wd")
plt.ylabel(list[i])
plt.savefig("./wd/"+list[i]+"_wd"+".jpg")
plt.show()
# In[7]:
for i in range(6):
plt.figure(figsize=(20,8))
plt.subplot(3,2,i+1)
plt.scatter(sd,y[i],c=colors[i])
plt.xlabel("sd")
plt.ylabel(list[i])
plt.savefig("./sd/"+list[i]+"_sd"+".jpg")
plt.show()
# In[8]:
for i in range(6):
plt.figure(figsize=(20,8))
plt.subplot(3,2,i+1)
plt.scatter(qy,y[i],c=colors[i])
plt.xlabel("qy")
plt.ylabel(list[i])
plt.savefig("./qy/"+list[i]+"_qy"+".jpg")
plt.show()
# In[9]:
for i in range(6):
plt.figure(figsize=(20,8))
plt.subplot(3,2,i+1)
plt.scatter(fx,y[i],c=colors[i])
plt.xlabel("fx")
plt.ylabel(list[i])
plt.savefig("./fx/"+list[i]+"_fx"+".jpg")
plt.show()
# In[10]:
for i in range(6):
plt.figure(figsize=(20,8))
plt.subplot(3,2,i+1)
plt.scatter(fs,y[i],c=colors[i])
plt.xlabel("fs")
plt.ylabel(list[i])
plt.savefig("./fs/"+list[i]+"_fs"+".jpg")
plt.show()
# In[11]:
from sklearn.decomposition import PCA
# In[12]:
columns = ["wd","sd","qy","fs","fx"]
tzdata = np.array([wd,sd,qy,fs,fx]).T
type(tzdata)
# In[13]:
pac_sum = []
for i in range(1,6):
pca = PCA(i).fit(tzdata)
pac_sum.append(pca.explained_variance_ratio_.sum())
pac_sum
# In[14]:
wd1 = pac_sum[0]
sd1 = pac_sum[1]-pac_sum[0]
qy1 = pac_sum[2]-pac_sum[1]
fs1 = pac_sum[3]-pac_sum[2]
fx1 = pac_sum[4]-pac_sum[3]
list = [wd1,sd1,qy1,fs1,fx1]
list
# In[15]:
import numpy as np
plt.rcParams['font.family']=['sans-serif']
plt.rcParams['font.sans-serif']=['SimHei']
plt.figure(figsize=(20,8))
plt.plot(range(1,6),pac_sum)
plt.xticks([1,2,3,4,5])
plt.ylim(0.9,1.01)
plt.xlabel("选择特征个数")
plt.ylabel("方差贡献率误差总和")
plt.savefig("tz_"+"var_sum"+".jpg")
plt.show()
# In[16]:
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score,silhouette_samples
# In[17]:
#查看样本标签
x = KMeans(5).fit(co.reshape(-1,1))
x.labels_
# In[18]:
#查看质心
x.inertia_
# In[19]:
x.cluster_centers_
# In[20]:
#测试打分情况
silhouette_score(co.reshape(-1,1),x.labels_)
# In[21]:
co.reshape(-1,1)
# In[22]:
silhouette_samples(co.reshape(-1,1),x.labels_)
# In[23]:
#确定分类的个数，为5-6类最好
list_inertia = []
for i in range(2,10):
x = KMeans(n_clusters=i).fit(co.reshape(-1,1))
list_inertia.append(x.inertia_)
plt.figure(figsize=(25,8))
plt.plot(range(2,10),list_inertia)
plt.xlabel("n")
plt.ylabel("inertia")
plt.savefig("km_"+"co"+".jpg")
plt.show()
# In[24]:
def cs_ct(sp):
x = KMeans(6).fit(sp.reshape(-1,1)).labels_
return x
# In[25]:
cs_ct(co)
# In[26]:
cl = []
for i in cs_ct(co):
cl.append(str(i))
cl
# In[27]:
m = cs_ct(co)
m.shape
# In[28]:
co
# In[29]:
#测试区间最大值以及最小值
data = pd.Series(data=co,index=cs_ct(co))
data[1].max()
data[1].min()
# In[30]:
np.array(data[0])
# In[31]:
type(data)
# In[32]:
#获取分类区间
def qj_lb(sp):
x = KMeans(6).fit(sp.reshape(-1,1)).labels_
data = pd.Series(data=sp,index=x)
list1 = []
list2 = []
for i in range(6):
list1.append(data[i].min())
list2.append(data[i].max())
list3 = []
for i in range(6):
print("第{}个类别的范围为：{}到{}".format(i+1,list1[i],list2[i]))
#list3.append("第{}个类别的范围为：{}到{}".format(i+1,list1[i],list2[i]))
# for i in range(6):
# plt.figure(figsize=(20,8))
# plt.subplot(3,2,i+1)
# plt.plot(np_sj,sp)
# plt.xlabel("lb")
# plt.ylabel("fw")
# plt.show()
print("-"*50)
# In[33]:
#遍历循环，打印分类数据
sj =[so2,no2,pm10,pm25,o3,co,wd,sd,qy,fs,fx]
bq =["so2","no2","pm10","pm25","o3","co","wd","sd","qy","fs","fx"]
for i in range(len(bq)):
print(bq[i]+"浓度影响")
qj_lb(sj[i])
#写入文件格式不对，未写入，故注释
#f = open(r"./nd.txt","w")
#f.write(str(qj_lb(sj[i])))
# In[112]:
def zx_wz(sp):
x = KMeans(6).fit(sp.reshape(-1,1)).labels_
data = pd.Series(data=sp,index=x)
list1 = []
list2 = []
list3=[]
for i in range(6):
list1.append(data[i].min())
list2.append(data[i].max())
return (list1,list2)
# In[ ]:
# In[84]:
zx_wz(co)
# In[79]:
def zx_d(sj):
list=[]
for i in range(6):
list.append(zx_wz(sj)[1][i]-zx_lb(co)[0][i])
return list
# In[87]:
zx_d(co)
# In[111]:
colors = ["red","blue","green","yellow","black","purple","orange","pink","coral","gold"]
for i in range(10):
plt.figure(figsize=(60,25))
plt.subplot(2,5,i+1)
plt.scatter(range(1,7),zx_d(sj[i]),c=colors[i])
plt.xlabel(bq[i])
plt.ylabel("分类质心位置")
plt.title(bq[i]+"分类图")
plt.savefig("./分类/"+bq[i]+"分类.jpg")
plt.grid()
plt.show()
图片展示：
方差
特征
so2浓度
第1个类别的范围为：7.0到9.0
第2个类别的范围为：-2.0到4.0
第3个类别的范围为：18.0到47.0
第4个类别的范围为：13.0到17.0
第5个类别的范围为：5.0到6.0
第6个类别的范围为：10.0到12.0
--------------------------
no2浓度
第1个类别的范围为：33.0到51.0
第2个类别的范围为：77.0到112.0
第3个类别的范围为：2.0到18.0
第4个类别的范围为：52.0到76.0
第5个类别的范围为：113.0到211.0
第6个类别的范围为：19.0到32.0
--------------------------
pm10浓度
第1个类别的范围为：82.0到116.0
第2个类别的范围为：25.0到40.0
第3个类别的范围为：-4.0到24.0
第4个类别的范围为：41.0到58.0
第5个类别的范围为：59.0到81.0
第6个类别的范围为：117.0到217.0
--------------------------
pm25浓度
第1个类别的范围为：25.0到37.0
第2个类别的范围为：-5.0到12.0
第3个类别的范围为：54.0到79.0
第4个类别的范围为：80.0到163.0
第5个类别的范围为：13.0到24.0
第6个类别的范围为：38.0到53.0
--------------------------
o3浓度
第1个类别的范围为：188.0到405.0
第2个类别的范围为：-2.0到22.0
第3个类别的范围为：84.0到128.0
第4个类别的范围为：51.0到83.0
第5个类别的范围为：129.0到187.0
第6个类别的范围为：23.0到50.0
--------------------------
co浓度
第1个类别的范围为：0.7到0.8
第2个类别的范围为：1.1到1.4
第3个类别的范围为：0.1到0.5
第4个类别的范围为：0.9到1.0
第5个类别的范围为：1.5到2.5
第6个类别的范围为：0.6到0.6
--------------------------
wd浓度
第1个类别的范围为：23.9到27.6
第2个类别的范围为：14.8到19.6
第3个类别的范围为：31.2到38.2
第4个类别的范围为：19.7到23.8
第5个类别的范围为：5.8到14.7
第6个类别的范围为：27.7到31.1
--------------------------
sd浓度
第1个类别的范围为：75.0到84.0
第2个类别的范围为：42.0到56.0
第3个类别的范围为：85.0到99.0
第4个类别的范围为：14.0到41.0
第5个类别的范围为：57.0到66.0
第6个类别的范围为：67.0到74.0
--------------------------
qy浓度
第1个类别的范围为：1003.0到1006.9
第2个类别的范围为：1011.3到1015.5
第3个类别的范围为：1007.0到1011.2
第4个类别的范围为：1020.1到1029.2
第5个类别的范围为：1015.6到1020.0
第6个类别的范围为：993.5到1002.9
--------------------------
fs浓度
第1个类别的范围为：0.7到1.0
第2个类别的范围为：2.1到2.6
第3个类别的范围为：1.1到1.5
第4个类别的范围为：2.7到5.8
第5个类别的范围为：0.1到0.6
第6个类别的范围为：1.6到2.0
--------------------------
fx浓度
第1个类别的范围为：40.4到84.0
第2个类别的范围为：261.7到311.7
第3个类别的范围为：174.1到261.6
第4个类别的范围为：84.1到172.8
第5个类别的范围为：0.0到40.3
第6个类别的范围为：311.8到360.0
--------------------------
后面的没做，就不展示。
不管咋样，这次竞赛收获还是很多的，有付出才有收获，加油。
相关文章
- 3D游戏建模真的很累吗？前景怎么样？又需要哪样的基础
- 数学建模：绘图经验总结
- 【Manning新书】Python时间序列预测，手把手教你实战时序建模
- 数据建模方法及步骤图_comsol建模步骤教程
- 数学建模中的选址问题_数学建模停车场规划问题
- BricsCAD 23 for Mac(CAD建模软件)v23.1.05永久激活版
- 数学建模的一些方法_对数学建模的认识
- 数学建模——熵权法
- 数据建模的精华：很少有人真正理解数据模型的形态
- [Medical Image Analysis | 论文简读] Factorizer: 医学图像分割的可扩展可解释上下文建模方法
- 深度网络CTR建模
- 用户多兴趣建模MIND
- 数学建模常用模型02：插值与拟合
- 数学建模常用模型03：层次分析法
- 2023年（美赛）美国大学生数学建模竞赛ABCDEF题思路资料汇总【全网最全】
- DAZ Studio Pro Edition 4.15数字人物三维建模与动画制作佼佼者
- ReCap Pro 2019三维建模和扫描数据转换的专业工具+全版本安装包
- 3DSMax如何下载安装？三维动画建模软件 3DSMax2024中文版下载
- SketchUp建筑建模软件，Sketch Up软件下载安装
- R语言使用马尔可夫链对营销中的渠道归因建模|附代码数据
- AI浪潮下，语音识别建模技术的演进 | 雷锋网公开课
- Art of Illusion是一个用Java语言写成的基于GPL v2的开源三维建模和渲染软件，可以创建高质量的具有真实感的三维模型。它使用方便，可以对纹理和材料等效果进行编辑，可以产生具有真实感的图像。