2020年数维杯数学建模C题 垃圾转运优化模型设计求解全过程文档及程序

2020年数维杯数学建模

C题 垃圾转运优化模型设计

原题再现：

  随着我国人口的不断增加及城镇化进程的快速推进，城市面临了众多公共管理方面的难题。如生活垃圾、废气废水及排泄物等等的处理问题。截止2019年底我国拥有十多个千万规模以上的大型城市，城镇人口数量达到了8.48亿人。

  数据统计结果表明我国的人均日垃圾产量为1.0至1.2公斤。每天大型城市都有上万吨规模的垃圾需要得到妥善的处理，将这些垃圾统一运输到垃圾处理中心至少需要几百辆车。在难以降低城市人均垃圾产量的情况下，如果不对城市垃圾进行有效处理，那么很快就会有大量的城市周围土地被占用，同时因垃圾而导致城市空气质量的迅速下降。对此国内外对垃圾处理问题展开了众多的研究，也出现了众多垃圾处理技术。

  通常生活垃圾分为可回收垃圾和不可回收垃圾，也有干垃圾和湿垃圾的分类方法。在我国开始施行了诸如可回收和不可回收及干垃圾和湿垃圾分类投放措施。然而，在居民对垃圾分类知识的匮乏及垃圾分类投放的不便捷性，垃圾分类投放的措施依然没有得到全面的推广应用。

  请您利用数学建模的方法去解决如下四个问题：

  问题1：搜集整理某个城市的交通网络数据及现有垃圾转运中心数据，并结合此城市城区地理位置、人口等信息预测出城市关键节点处的垃圾每日产量数据；

  问题2：提供垃圾转运站的合理选址结果，并根据指定的单个或多个垃圾处理中心位置与处理能力，建立在不同车容量下的垃圾转运优化模型与算法；

  问题3：通常每天能够正常工作的垃圾运输车辆或多或少，请您提供在指定的垃圾运输车辆总数下的最优垃圾运输任务指派方案，并尽量按照垃圾运输司机的运输意愿完成任务指派；

  问题4：如果同时考虑垃圾的类型，请您重新对问题2和问题3进行建模，并举例说明。我们特别期待深入讨论大规模情形下的优化模型与算法的论文。

整体求解过程概述(摘要)

  随着我国人口的增加和城市化进程的加快，城市生活垃圾的数量呈现逐年较快增长的趋势，因此处理处理问题变得尤为重要，而收集运输环节是城市生活垃圾处理的重要一环，严重影响到垃圾处理的经济效益、社会效益和环境效益。如今城市生活垃圾的研究和管理较多停留在末端治理和源头控制方面，而对于中间收集运输环节的关注相对较少，垃圾收运车辆的调度和路径优化问题的研究难以进行实际应用。现有研究目标较为单一且对多车型的研究较少，并且很少有进行分类运输的研究。因此本文以垃圾收运路径问题的研究缺口入手，提出了基于城市生活垃圾分类的清运车辆调度。
  基于岭回归模型和GM（1，N）模型，建立了城市生活垃圾产量预测的多元组合模型，并以青岛市为例进行了生活垃圾的产量预测。基于岭回归模型和GM（1，N）模型建立的多元组合模型具有很好的预测效果，其预测精度和变量解释程度高于单个模型，同时通过权重系数的分配，降低了单一模型发展趋势的单一性、定向性，综合了单个模型的自身趋势，使得预测曲线在拟合和预测的节点上更加平滑。通过对城币生沽垃圾收运系统进行定性和定量分析，结合青岛城市生活垃圾清运现状，描述了现有收运路径并分析存在的问题，对中转站设备配置以及清运车辆调度进行研究。主要分两部分，一部分是对现有混装清运模式的车辆调度进行优化，另一部分是研究分类清运的车辆调度。首先结合排队论优化中转站卸料口的配置并进行经济效益分析;然后依据垃圾收运模式的特点，建立模型，将实际问题中的限制转换成约束条件，提出垃圾分类清运的车辆调度模式，模型求解方面，使用了粒子群算法在MATLAB软件上进行编程求解，得到不同模式下清运车辆的调度方案。

问题分析：

  青岛，别称岛城、琴岛、胶澳，是山东省副省级市、计划单列市，国务院批复确定的中国沿海重要中心城市和滨海度假旅游城市、国际性港口城市。截至2019年，全市下辖7个区、代管3个县级市，总面积11293平方千米，2019年全市常住总人口949.98万人。其中，市区常住人口645.20万人。

  青岛市政府于2019年11月18日出台《青岛市城市生活垃圾管理办法》，该办法明确提出:生活垃圾治理规划应当结合生活垃圾的产量预测和成分特点，确定生活垃圾管理事业的发展方向，统筹生活垃圾处置流量、流向，明确生活垃圾处置结构和设施总体布局。在生活垃圾处理方式上，垃圾清洁填埋目前占据了绝大部分比重;在垃圾的运输方式上，青岛推行清洁直运的垃圾运输政策。其正在运行的二期填埋场设计处理能力2671吨/日，而目前实际处理数量超过6000吨/日，为设计处理能力的两倍以上，处于超负荷运作。因此做好生活垃圾产量的预测工作对垃圾的运输、处置和规划具有重要意义。以2011-2019年青岛市的生活垃圾实际清运量作为生活垃圾的产量数据，数据来源于青岛环境保护局历年发布的青岛市环境状况公报。相关影响因子的历年数据来源于青岛市统计局历年发布的青岛统计年鉴中的人口与就业人员、交通运输、国内贸易、旅游、城市建设、环境保护、教育、物价等章节。

模型假设：

  目前的混装清运模式是使用配备移动式容器的清运车辆，属于垃圾收运模型中的拖曳模型。清运车辆路径是从停车场出发，去往垃圾楼，置换箱体后原路返回中转站卸载垃圾箱的过程，在允许的时间范围内，箱式车可以再次从中转站出发，去往下一个垃圾楼进行清运，重复作业多次最终在规定时间窗内返回停车场。清运模式见下图 3-3，虚线代表装载空集装箱的路径，实线代表装载垃圾集装箱的路径，其中 3-3(a)为停车场与中转站在不同位置的情况，3-3(b)表示停车场与中转站在同一位置的情况。考虑到模型的适用性，建立多个停车场、多个中转站、停车场与中转站位置不同的车辆调度模型，而在计算模型时，考虑到北京市的实际情况，转换为停车场与中转站位置相同的车辆调度模型。
  (1)各垃圾楼垃圾量己知并且固定，不会随着时间的变化而发生增减，为更加贴近实际，早班和晚班的垃圾量分别计算;
  (2)垃圾楼的开放时间与清运工作时间一致;
  (3)清运车辆每次清运一箱垃圾，每个垃圾楼的垃圾可以由多辆车或一辆车多次
清运;
  (4)清运车辆数量充足;
  (5)除高峰期时间段外，车辆在每条线路上的行驶速度一致;
  (6)清运车辆收集垃圾后原路返回中转站，往返行驶距离相同;
  (7)车辆根据到达中转站的时间先到先服务;
  (8)运输成本只与运输距离呈正相关关系。

模型的建立与求解

  对于诸多的影响因子，需要用一定的方式进行甄别和筛选。本文采用邓氏关联度理论对预选的变量进行关联度分析，保留关联度大的影响因子参与建模。邓氏关联分析属于色关联分析，其基本思想是根据序列曲线几何形状的相似程度来判断其联系是否紧密。线越接近，相应序列之间的关联度就越大。

岭回归模型

  为了保证模型数据源与模型自身的统一性，本文采用一维模型GM（1，1）和一元线性函数生成了相关影响因子未来5年的预测值，分别作为GM（1，N）和岭回归模型预测数据来源，预测值及拟合情况见表3、表4。


  以大型设备至所辖转运站的加权运输距离S最短(式2-1)各大型生活垃圾处理设备的每日垃圾处理量E的标准差最小（式2-2）及各大型设备所在地的恶臭影响程度最低(式2-3)为目标函数，以式(2-4)为约束条件，建立大型生活垃圾处理设备选址模型如下

  大型生活垃圾处理设备选址的数学模型是一个多目标优化问题，直接求解较复杂，因此通过无量纲化将量纲统一，并将fi加权求和转换为一个目标函数F，以便于模型的求解。

小型处理设备选址

  在主要考虑经济效益的前提下，大型设备能将每日城市生活垃圾总量M全部处理完，但较长的运输距离将增加垃圾恶臭影响的环保成本。为了量化垃圾恶臭影响的环保成本，假设单位运输时间，单位垃圾量的环保成本为EC元。下面将对大型生活垃圾处理设备的总成本BC与小型生活垃圾处理设备的总成本SC进行比较，得出小型设备选址的一般性结论。

  2) 小型设备选址的一般性结论
  当第i个转运站至第j个大型生活垃圾处理设备的加权距离S与转运站的生活垃圾量RL满足小型生活垃圾处理设备的总成本SC小于大型设备的总成本BC时, 即可在该转运站设置小型生活垃圾处理设备。

  分析表5中的数据,发现大、小型设备的处理能力相差过大。故在满足大小设备处理能力上下限的约束下,比较大、小型设备在不同处理能力下总成本差异,重新计算最优的设备处理能力,使其满足不同城市生活垃圾的收运与处理。

生活垃圾处理设备处理能力的函数关系

大小型生活垃圾处理设备的数量与选址

模型求解

模型求解

  模型以原宣武区的垃圾清运作业为研究对象，在实际问题中车场和中转站距离很近忽略不计，而中转站的数量为 1 个，即马家楼中转站。根据垃圾清运模式的特点，使用粒子群算法进行求解。粒子群算法容易陷入局部最优解，因此引入爬山算法，通过追踪鸟群群首路径中的最优解，调整整体走向。
  1)求解步骤:
  设置每个粒子的初始位置限制和粒子速度限制。粒子位置限制设为[0,1]，速度限制设为[-0.1,0.1]。根据粒子当前位置适应度函数，将每个粒子的目标函数值作为其适应度值。为了防止出现局部最优的情况，引入小型爬山算法，对粒子位置进行变异计算。对于每个粒子，适应度值与其经过的个体历史最优值位置比较，如果该值优于历史最优值则保留，该值取代个体历史最优值，进入下一步迭代。对于每个粒子，适应度值与粒子种群经历过的最优位置进行比较，如果该值优于粒子群的历史最优值则保留，该值取代个体历史最优值，否则仍保留历史最优值，进入下一步迭按位置和速度的更新对粒子种群中的每个粒子进行迭代更新从而产生新的进化种群。

表 6 垃圾楼相关数据

  本小节的模型结果首先对混装清运调度目前的经验方案进行优化研究，优化后仍存在一定的等待时间，因此考虑第三章中转站卸料口配置优化的结果对模型进行求解，同时也提出了清运车辆的调度方案。车辆调度可以参考下表 7 执行。

  本章在分析垃圾混装清运的现状的基础上，提出清运车调度存在的问题，针对问题进行混装清运车辆调度模型。混装清运模型是垃圾收运过程中使用移动式容器收运的 Rollon-Rolloff 问题，可以抽象成带时间窗、带中间处理设施、多个车场的 VRP 问题。模型设立四个目标函数，通过归一和加权得到总体目标函数;根据混装清运特点对模型进行约束。以马家楼中转站及其对应的原宣武区的垃圾清运作业为例，采用粒子群算法对模型进行求解，并提出车辆调度方案，优化原有车辆调度。

问题四模型假设

  (1)各垃圾楼垃圾量已知并且固定，不会随着时间的变化而发生增减，为更加贴近实际，早班和晚班的垃圾量分别计算;
  (2)将每个社区作为一个投放点，各垃圾种类的比例相同，个别放点根据特点进行微调;
  (3)可回收物使用箱式车清运，厨余垃圾和其他垃圾使用分仓清运车清运，车辆数量充足;
  (4)分仓清运车的各仓载重量固定;
  (5)在载重和时间允许情况下，车辆可去往多个地点进行清运;
  (6)道路没有转向、单行等限制
  (7)垃圾楼的开放时间与清运工作时间一致;
  (8)除高峰期时间段外，车辆在每条线路上的行驶速度一致;
  (9)节点之间的行驶距离固定;
  (10)车辆根据到达中转站的时间先到先服务;
  (11)运输成本只与运输距离呈正相关关系。

模型 MT1 的建立

  在分类联合清运的模式中，可以将整体模型看作是由两个模型构成，可回收物的清运是单独的路径，作为其中一个模型;另外一个模型是厨余垃圾和其他垃圾的联合清运的路径模型，建模如下:
  （1）目标函数

  (2)约束条件

论文缩略图：

程序代码：

import random as ra
import numpy as np
import adj_matrix_china as adj
import copy
num=31
#x0:经济 y0:禀赋 z0:初始效益
x0=np.array([126.79 ,123.88 ,77.94 ,89.61 ,107.80 ,\
 93.44 ,51.66 ,79.30 ,60.81 ,67.49 ,60.97\
 ,45.34 ,52.71 ,48.78 ,36.09 ,39.45 ,\
 35.02 ,46.81 ,39.78 ,31.86 ,44.70 ,41.86 \
 ,49.26 ,44.20 ,\
38.12 ,26.79 ,42.00 ,36.61 ,27.14 ,34.22 ,32.75])
"""
x0=np.array([100.00 ,95.91 ,55.10 ,65.62 ,75.37 ,57.42 ,34.87 ,65.33 ,47.19 ,\
 43.08 ,46.24 ,34.46 ,41.37 ,35.09 ,26.51 ,27.88 ,22.06 ,\
 35.67 ,28.26 ,26.20 ,34.38 ,32.42 ,40.64 ,34.01 ,33.10 ,\
 29.82 ,33.06 ,29.87 ,20.12 ,29.23 ,28.31])
 """
y0=np.array([0.21900552,0.150779763,0.224193101,0.143499445,0.169196543,0.178010654
,0.272209386,\
0.194657663,0.262827531,0.23860188,0.220629418,0.270156393,0.320706923,0.276874761
,\
0.317457365,0.32817422,0.391237114,0.303330628,0.331573028,0.298208533,0.27237412,\
0.303493275,0.26215307,0.348405972,0.329620889,1,0.279305443,0.44507755,0.55094282
4,\
0.324156628,0.352843656])
y0=y0*110
z0=x0/y0
x1=np.copy(x0)
y1=np.copy(y0)
z1=np.copy(z0)
#中间变量
x2=np.zeros(num)
y2=np.zeros(num)
z2=np.zeros(num)
#单位步长
h=5
#效益权重
w1=1/3
w2=0.2 
#合作矩阵 aij:i 对 j 提供效益帮助 1 单位：1 i 对 j 提供资源帮助 1 单位：1j
A=np.zeros((num,num),dtype=complex)
#邻接矩阵
B=np.zeros((num,num))
B=copy.deepcopy(adj.A)
print('测试邻接矩阵')
print('江西','湖北',adj.Is_adj('江西','湖北'))
print('江西','北京',adj.Is_adj('江西','北京'))
print('江苏','河南',adj.Is_adj('江苏','河南'))
print('江西','贵州',adj.Is_adj('江西','贵州'))
"""
#test num=3
B[0,1]=B[1,0]=1
B[1,2]=B[2,1]=1
B[0,2]=B[2,0]=1
"""
#test num=8
"""
B[0,1]=B[1,2]=B[2,4]=B[4,7]=1
B[6,7]=B[5,6]=B[3,5]=B[0,3]=1
"""
for i in range(0,num):
 for j in range(i+1,num):
 B[j,i]=B[i,j]
def cal_benkm(zk1,yk1,zk0,yk0):#m 对 k 带来的好处
 return zk1*yk1-zk0*yk0
def cal_zm1(zm,zk):#k 对 m 的效益提升 返回 zm'
 zm1=max(w1*zk,zm+w2*(zk-zm))
 return zm1
 
 
def judge1(k,m):
 zk0=z1[k]
 yk0=y1[k]
 zm0=z1[m]
 ym0=y1[m]
 if(y0[k]<y0[m]-h):#m 可以给 k 提供资源
 ym1=ym0-h
 zm1=cal_zm1(zm0,zk0)
 ben_m=cal_benkm(zm1,ym1,zm0,ym0)
 if(ben_m>0):
 return 1#双方有益
 if(zm0<zk0*w1):#k 帮扶 m,无需代价
 return -1#单方帮扶
 return 0#拒绝合作
 
 
def add_collab(k,m,s):#增强合作
 #参数 s:s=2 k 经济帮 m，m 资源帮 k;s=-2 k 与 m 调换
 # s=1 k 经济帮扶 m s=-1 m 经济帮扶 k
zk0=z1[k]
 yk0=y1[k]
 zm0=z1[m]
 ym0=y1[m]
 if(s==2):
 A[k,m]+=1
 A[m,k]+=1j
 y1[m]=ym0-h
 y1[k]=yk0+h
 z1[m]=cal_zm1(zm0,zk0)#k 帮 m 对 m 的效益提升
 elif(s==-2):
 A[k,m]+=1j
 A[m,k]+=1
 z1[k]=cal_zm1(zk0,zm0)
 y1[k]=yk0-h
 y1[m]=ym0+h
 elif(s==1):
 A[k,m]+=1
 z1[m]=cal_zm1(zm0,zk0)
 elif(s==-1):
 A[m,k]+=1
 z1[k]=cal_zm1(zk0,zm0)
#返回 k 对 m 的合作状态
# s=2 k 经济帮 m，m 资源帮 k;s=-2 k 与 m 调换
# s=1 k 经济帮扶 m s=1 s 经济帮扶 k) 
def stats(k,m):
 a=A[k,m]
 b=A[m,k]##?会不会 a b 都有实虚部?
 if(a.real>0 and b.imag>0):
 return 2
 elif(b.real>0 and a.imag>0):
 return -2
 elif(a.real>0 and abs(b)<1e-5):
 return 1
 elif(b.real>0 and abs(a)<1e-5):
 return -1
 return 0
 
#main program:
def find():
 i=ra.randint(0,num-1)#任选一省份
 zi0=z1[i]
 yi0=y1[i]
 max_benf_m=0
 flag=0
s=0
 #参数 s:s=2 i 经济帮 k，k 资源帮 i;s=-2 k 与 i 调换
 # s=1 i 经济帮扶 k s=-1 k 经济帮扶 i
 for k in range(0,num):
 
 if(B[i,k]!=0):
 yk0=y1[k]
 zk0=z1[k]
 benf_ik=0
 benf_ki=0
 if(zi0>zk0):#i 可经济帮 k
 status=judge1(i,k)
 if(status==1):
 zk1=cal_zm1(zk0,zi0)#后者对前者的效益提升 返回前者修改值
 yk1=yk0-h
yi1=yi0+h
benf_ik=cal_benkm(zi0,yi1,zi0,yi0)
benf_ki=cal_benkm(zk1,yk1,zk0,yk0)
ss=2
 elif(status==-1):
 zk1=cal_zm1(zk0,zi0)
benf_ki=cal_benkm(zk1,yk0,zk0,yk0)
ss=1
 else:#k 可经济帮 i
 status=judge1(k,i)
 if(status==1):
 zi1=cal_zm1(zi0,zk0)
yi1=yi0-h
yk1=yk0+h
benf_ik=cal_benkm(zi1,yi1,zi0,yi0)
benf_ki=cal_benkm(zk0,yk1,zk0,yk0)
ss=-2
 elif(status==-1):
 zi1=cal_zm1(zi0,zk0)
 benf_ik=cal_benkm(zi1,yi0,zi0,yi0)
 ss=-1
 if((benf_ik+benf_ki)>max_benf_m):
 max_benf_m=benf_ik+benf_ki
 flag=k
 s=ss
 #end if
 #end if of B[i,k]!=0
 if(flag!=0):
 add_collab(i,flag,s)
def cal_ben_improv(z0,z1):
 s=0
 for i in range(0,num):
 s+=cal_benkm(z1[i],y1[i],z0[i],y0[i])
print(s)
 
for i in range(0,10000):
 find()
cal_ben_improv(z0,z1)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188