Python
7.20
openpyxl
- 读取某一个指定行列的数值
from openpyxl import load_workbook
# 加载同一个文件夹下面的一个excel文件
wb = load_workbook("1.xlsx")
sheet = wb["AUTO3006-001"]
res = sheet.cell(row=1, column=8).value
for i in range(2, 20):
res = sheet.cell(row=i, column=8).value
print(res)
- 读取某一列的数，并存放到列表当中去
from openpyxl import load_workbook
# 加载同一个文件夹下面的一个excel文件
wb = load_workbook("1.xlsx")
# 读取表中的某一个sheet
ws = wb["AUTO3006-001"]
# 建立一个数组来放置数据
res = []
# 读取一列的值
for i in ws["H"]:
# 读取sheet返回的是一个列表，取值的时候用value去读取
res.append(i.value)
print(res)
number = []
for i in ws["A"]:
# 读取sheet返回的是一个列表，取值的时候用value去读取
number.append(i.value)
print(number)
seaborn
- 知乎可视化（一）：核密度估计图kdeplot
import matplotlib.pyplot as plt
import numpy as np
# 导入numpy包，用于生成数组
import seaborn as sns
# 习惯上简写成sns
x = np.random.randn(100)
# 随机生成100个符合正态分布的数
sns.kdeplot(x, cut=0.2)
# cut是什么意思啊，我还没看出来
plt.show()
- 知乎可视化（一）：distplot集合了matplotlib的hist()与核函数估计kdeplot的功能，增加了rugplot分布观测条显示与利用scipy库fit拟合参数分布的新颖用途。
- 知乎可视化（四）：stripplot(分布散点图)
7.22
Openpyxl
- 读取sheetname，修改sheetname
from openpyxl import load_workbook
# 读取sheetname
wb = load_workbook('1.xlsx')
# 或者不知道名字时
sheet_names = wb.sheetnames # 返回一个列表
ws1 = wb[sheet_names[0]] # 第一个sheet的名字
ws2 = wb[sheet_names[1]] # 第二个sheet的名字
ws3 = wb[sheet_names[2]] # 第三个sheet名字
print(ws1)
print(ws2)
print(ws3)
# 修改sheet的名字
ws1.title = 'AUTO'
wb.save("1.xlsx")
# 要在文件处于关闭状态才能修改之后覆盖原文件
print(ws1)
- 获取单元格信息
from openpyxl import load_workbook
# 打开excel文件
wb = load_workbook('1.xlsx')
# 获取活跃表对象
sheet = wb.active
# 获取单元格对应的 Cell 对象
a6 = sheet['A6'] # A1 表示A列中的第一行，这儿的列号采用的是从A开始的
print(a6)
# 获取单元格中的内容
content = a6.value
print(content) # 结果是: Rank
# 获取单元格的行和列信息
row = a6.row
print('行:', row) # 结果： 1
column = a6.column
print('列:', column) # 结果： 6
coordinate = a6.coordinate
print(coordinate) # 结果：A6
7 23
- 一次读取一个文件多个sheet，同一列的数据
from openpyxl import load_workbook
import time
# 记录开始时间
start = time.time()
# 加载文件
wb = load_workbook("1.xlsx")
# ===================================
def Mysleeptime(hour, min, sec):
return hour * 3600 + min * 60 + sec
second = Mysleeptime(0, 0, 1)
# 这个second用来指定延迟函数的时间长短，可延可不延的
# print(second)
# >>>>>>>>>>>>>>>>>>>>>>>>>>批量读取表格内容<<<<<<<<<<<<<<<<<<<<
# >>>>>>>>>>>>>>>>>>>>>>>>>>读取A1-A3的内容<<<<<<<<<<<<<<<<<<<<
for i in range(4):
if i != 0:
k = "AUTO%s" % i # 用k来表示当前的sheet
print(k) # 拼接字符
res = []
ws = wb[k]
flag = 1 # 标志当前是第几个sheet
# 读取当前sheet中指定列的值
for j in ws["A"]:
# 读取sheet返回的是一个列表，取值的时候用value去读取
res.append(j.value)
print("读取%s" % k, res, sep='=')
flag = flag + 1
time.sleep(second) # 延迟函数
else:
print("结束")
end = time.time() # 计算结束时间
print('Running time: %s Seconds' % (end - start))
# ===================================
7.24
字典 key
用于存放具有映射关系的数据
字典的键不能重复，且只能用不可变的数据作为字典的键
字典中的数据是没有先后顺序的，字典的存储是无序的
字典是一种可变的容器，可以存储任意类型的数据
-
使用字典的字面值创建字典，里面允许出现若干键值对
-
创建空字典 dic = {}
-
创建非空字典
dic = { 1: 'One', 2: 'Two', 'age1': 890, 'age9': 9890, }
-
-
使用字典的构造函数创建字典
- dic = dict() #创建一个空字典
- 使用可迭代对象初始化一个字典
dic = dict(
[
('name','tarena'),
('age',15),
['addr','BJ']
]
)
# 这里是几个意思？？？？
- 使用关键字传参形式创建一个字典
dic=(name='tarena',age=15)
- 修改健对应的值
a = {
'w': 1,
'h': 2,
3: 900990}
print(a[3])
a[3] = 'yyyyyy'
# 键存在的时候直接替换
a[9] = 'wwww'
print(a)
# 键不存在的时候直接添加
导入ID和箱数 代码
tip：19-11-25没有“签收汇总”就把这个表移出去了
from openpyxl import load_workbook
import os
import time
import pandas as pd
# 2020/7/24
# 记录开始时间
start = time.time()
# ---------------------------------------------------------------------------------
# 读取excel文件的相应列的值，存入列表
# 输入:文件名，表单名，列名称
# 输出：读出该列从第二行开始的值并存入列表，输出列表
# ---------------------------------------------------------------------------------
def ReadColumn(Document, Sheet, Column):
Wb = load_workbook(Document)
Ws = Wb[Sheet]
output = []
for i in Ws[Column]:
row = i.row
if row > 1:
output.append(i.value)
return output
# ----------------------------------------------------------------------------------
# ==============================19年10月数据导入 data19_10============================
data19_10 = {}
for i in range(1, 32):
x = "19年10月\广医市配2019-10-%s.xlsx" % i
y = "签收汇总"
o = "O"
g = "G"
# 有些天是没有相应文件的，判断是否存在该文件
if not os.path.exists(x):
continue
a = ReadColumn(x, y, o)
b = ReadColumn(x, y, g)
temp = {
i: {}
}
length = len(a)
for j in range(0, length):
q = a[j]
p = b[j]
# 判断原来是否有该ID，若有则累加，无则直接保存
if q not in temp[i]:
temp[i][q] = p
else:
temp[i][q] = p + temp[i][q]
# 判断是否q是否为空或0，是则结束遍历
if q is None:
break
# 将这一张表的key-value,以i为二级字典名，添入一级字典data
data19_10.update(temp)
print(temp)
end = time.time() # 计算结束时间
print('Running time: %s Seconds' % (end - start))
# ==============================19年11月数据导入 data19_11============================
data19_11 = {}
for i in range(1, 32):
x = "19年11月\广医市配2019-11-%s.xlsx" % i
y = "签收汇总"
o = "O"
g = "G"
# 有些天是没有相应文件的，判断是否存在该文件
if not os.path.exists(x):
continue
a = ReadColumn(x, y, o)
b = ReadColumn(x, y, g)
temp = {
i: {}
}
length = len(a)
for j in range(0, length):
q = a[j]
p = b[j]
# 判断原来是否有该ID，若有则累加，无则直接保存
if q not in temp[i]:
temp[i][q] = p
else:
temp[i][q] = p + temp[i][q]
# 判断是否q是否为空或0，是则结束遍历
if q is None:
break
# 将这一张表的key-value,以i为二级字典名，添入一级字典data
data19_11.update(temp)
print(temp)
end = time.time() # 计算结束时间
print('Running time: %s Seconds' % (end - start))
# ==============================19年12月数据导入data19_12============================
data19_12 = {}
for i in range(1, 32):
x = "19年12月\广医市配2019-12-%s.xlsx" % i
y = "签收汇总"
o = "O"
g = "G"
# 有些天是没有相应文件的，判断是否存在该文件
if not os.path.exists(x):
continue
a = ReadColumn(x, y, o)
b = ReadColumn(x, y, g)
temp = {
i: {}
}
length = len(a)
for j in range(0, length):
q = a[j]
p = b[j]
# 判断原来是否有该ID，若有则累加，无则直接保存
if q not in temp[i]:
temp[i][q] = p
else:
temp[i][q] = p + temp[i][q]
# 判断是否q是否为空或0，是则结束遍历
if q is None:
break
# 将这一张表的key-value,以i为二级字典名，添入一级字典data
data19_12.update(temp)
print(temp)
end = time.time() # 计算结束时间
print('Running time: %s Seconds' % (end - start))
差不多运行了1300s
天呐 21分钟
测试结果通过：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zi8Zb69X-1677154847151)(F:\2020Summer\暑期实习科研\实习\图片\导入数据成功.JPG)]
如上图，545433送货ID对应在表9中出现了3次，累加为165箱，实际字典中的value也为165箱
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-RiTs5yAr-1677154847152)(F:\2020Summer\暑期实习科研\实习\图片\导入数据成功2.JPG)]
如上图，确实是都存进去了哈，nb
7.27
经纬度转换距离
import math
# ===================================================================================
# 通过经纬度测算距离
# INPUTS: A点经度 A点纬度 B点经度 B点纬度
# OUTPUTS: 返回值为考虑球面的两点直线距离,单位：km
# ===================================================================================
def distance(lng_A, lat_A, lng_B, lat_B):
# 地球半径
R = 6371.004
pi = 3.141592654
Mlng_A = lng_A
Mlat_A = 90 - lat_A
Mlng_B = lng_B
Mlat_B = 90 - lat_B
C = math.sin(Mlat_A*pi/180) * math.sin(Mlat_B*pi/180) * math.cos((Mlng_A - Mlng_B)*pi/180) +math.cos(Mlat_A*pi/180) * math.cos(Mlat_B*pi/180)
if C >= 1:
C = 1
Distance = R * math.acos(C)
return Distance
d = distance(113.37935836, 31.71785761, 114.3933937, 30.63068048)
print(d)
向高德地图获取经纬度
# INPUTS：字符型的地址
# OUTPUTS：经度，纬度
def getcode(site):
parameters = {'address': site, 'key': '80115f36f89a4dd9fb3c4dd3911b149f'}
base = 'http://restapi.amap.com/v3/geocode/geo'
response = requests.get(base, parameters)
info_site = response.json()
lng = info_site['geocodes'][0]['location'].split(',')[0]
lat = info_site['geocodes'][0]['location'].split(',')[1]
lat = float(lat)
lng = float(lng)
return lng, lat
7.28
格式化处理文档一列内容
import requests
import math
from openpyxl import load_workbook
# 这个文档现在用来处理地址数据
import os
import time
# ---------------------------------------------------------------------------------
# 读取excel文件的相应列的值，存入列表
# 输入:文件名，表单名，列名称
# 输出：读出该列从第二行开始的值并存入列表，输出列表
# ---------------------------------------------------------------------------------
def ReadColumn(Document, Sheet, Column):
Wb = load_workbook(Document)
Ws = Wb[Sheet]
output = []
for i in Ws[Column]:
row = i.row
if row > 1:
output.append(i.value)
return output
# ----------------------------------------------------------------------------------
for i in range(1, 2):
wb = "19年11月\广医市配2019-11-%s.xlsx" % i
ws = "签收汇总"
o = "O"
g = "G"
k = "K"
# 有些天是没有相应文件的，判断是否存在该文件
if not os.path.exists(wb):
continue
a = ReadColumn(wb, ws, o)
b = ReadColumn(wb, ws, g)
place = ReadColumn(wb, ws, k)
book = load_workbook(wb)
sheet = book["签收汇总"]
for j in range(len(a)):
if place[j] is None:
break
if not ("佛山" in place[j]):
if not ("广州" in place[j]):
strlist = ['广州市', place[j]]
place[j] = ''.join(strlist)
print(place[j])
col = j+2
sheet.cell(row=col, column=11).value = place[j]
if not ("广东" in place[j]):
col = j+2
strlist = ['广东省', place[j]]
place[j] = ''.join(strlist)
print(place[j])
sheet.cell(row=col, column=11).value = place[j]
# 保存文件
book.save(wb)
json格式的读写
import json
# json.dump() 用于将dict类型的数据转成str，并写入到json文件中 , 相当于f.write(json.dumps())
dict = {'a': 'wo', 'b': 'zai', 'c': 'zhe', 'd': 'li'}
json.dump(dict,open(r'C:\Users\zy\Documents\GitHub\python3\searchTest\json.json','w'))
# json.load() 用于从json文件中读取数据，相当于f.read(json.loads())
filename = (r'C:\Users\zy\Documents\GitHub\python3\searchTest\json.json')
jsObj = json.load(open(filename))
print(jsObj)
# 字符类型
print(type(jsObj))
保存一个月的数据为三维字典
1级字典：19年x月
2级字典：1-31（代表某一天）
3级字典：{
‘amount’ :
‘lng’ :
‘lat’ :
}
import requests
import math
from openpyxl import load_workbook
# 这个文档现在用来处理地址数据
import os
import time
import json
# 记录开始时间
start = time.time()
# ---------------------------------------------------------------------------------
# 读取excel文件的相应列的值，存入列表
# 输入:文件名，表单名，列名称
# 输出：读出该列从第二行开始的值并存入列表，输出列表
# ---------------------------------------------------------------------------------
def ReadColumn(Document, Sheet, Column):
Wb = load_workbook(Document)
Ws = Wb[Sheet]
output = []
for i in Ws[Column]:
row = i.row
if row > 1:
output.append(i.value)
return output
# -------------------------------------------------------------------------------------
# INPUTS：字符型的地址
# OUTPUTS：经度，纬度
def getcode(site):
parameters = {'address': site, 'key': '80115f36f89a4dd9fb3c4dd3911b149f'}
base = 'http://restapi.amap.com/v3/geocode/geo'
response = requests.get(base, parameters)
info_site = response.json()
lng = info_site['geocodes'][0]['location'].split(',')[0]
lat = info_site['geocodes'][0]['location'].split(',')[1]
lat = float(lat)
lng = float(lng)
return lng, lat
# ===================================================================================
# 通过经纬度测算距离
# INPUTS: A点经度 A点纬度 B点经度 B点纬度
# OUTPUTS: 返回值为考虑球面的两点直线距离,单位：km
# ===================================================================================
def distance(lng_A, lat_A, lng_B, lat_B):
# 地球半径
R = 6371.004
pi = 3.141592654
Mlng_A = lng_A
Mlat_A = 90 - lat_A
Mlng_B = lng_B
Mlat_B = 90 - lat_B
C = math.sin(Mlat_A * pi / 180) * math.sin(Mlat_B * pi / 180) * math.cos((Mlng_A - Mlng_B) * pi / 180) + math.cos(
Mlat_A * pi / 180) * math.cos(Mlat_B * pi / 180)
Distance = R * math.acos(C)
return Distance
# ==================== =======19年10月数据导入 data19_10并保存为json==========================
data19_10 = {}
for i in range(1, 32):
wb = "19年10月\广医市配2019-10-%s.xlsx" % i
ws = "签收汇总"
o = "O"
g = "G"
k = "K"
# 有些天是没有相应文件的，判断是否存在该文件
if not os.path.exists(wb):
continue
a = ReadColumn(wb, ws, o)
b = ReadColumn(wb, ws, g)
place = ReadColumn(wb, ws, k)
temp = {
i: {}
}
length = len(a)
for j in range(0, length):
q = a[j]
# 判断是否q是否为空或0，是则结束遍历
if q is None:
break
p = b[j]
x, y = getcode(place[j])
# 判断原来是否有该ID，若有则累加，无则直接保存
if q not in temp[i]:
temp[i][q] = {
'amount': p,
'lng': x,
'lat': y
}
else:
temp[i][q]['amount'] = p + temp[i][q]['amount']
# 将这一张表的key-value,以i为二级字典名，添入一级字典data
data19_10.update(temp)
print(temp)
# json.dump() 用于将dict类型的数据转成str，并写入到json文件中 , 相当于f.write(json.dumps())
json.dump(data19_10, open(r'E:\PycharmProjects\Project1\19年10月\19_10.json', 'w'))
end = time.time() # 计算结束时间
print('Running time: %s Seconds' % (end - start))
# json.load() 用于从json文件中读取数据，相当于f.read(json.loads())
# filename = r'E:\PycharmProjects\Project1\19年10月\19_10.json'
# jsObj = json.load(open(filename))
# print(jsObj)
# print(type(jsObj))
Tips：索引的时候每一级都要加字符来索引，比如找第一个字典的‘1’日的送货ID317的总箱数
data19_10【‘1’】【‘317’】【‘amount’】才可以找到相应的数量
7.29
facility example
Model
$Minimize (Cy_j+\sum_i\sum_jD_{ij}x_{ij}) $
my_obj = C 价格系数向量
ub 上限
lb 下限
rhs 资源约束
senses 是什么
7.30
Cplex for Python – Day 1
（continuing to Friday）
CVRP
m i n i m i z e ∑ i , j ∈ A c i j x i j minimize\sum_{i,j\in A}c_{ij}x_{ij} minimize∑i,j∈Acijxij
s.t.
∑ j ∈ V , j ≠ i x i j = 1 , i ∈ N \sum_{j\in V,j \ne i}x_{ij}=1,i\in N ∑j∈V,j=ixij=1,i∈N
∑ i ∈ V , i ≠ j x i j = 1 , j ∈ N \sum_{i\in V,i \neq j}x_{ij}=1,j\in N ∑i∈V,i=jxij=1,j∈N
i f ( x i j = 1 ) , u i + q j = u j , i , j ∈ A : j ≠ 0 , i ≠ 0 if (x_{ij}=1),u_i+q_j=u_j,i,j\in A:j\neq 0,i\neq 0 if(xij=1),ui+qj=uj,i,j∈A:j=0,i=0
q i ≤ u i ≤ Q , i ∈ N q_i\le u_i\le Q,i\in N qi≤ui≤Q,i∈N
x i j ∈ 0 , 1 , i , j ∈ A x_{ij}\in{0,1},i,j\in A xij∈0,1,i,j∈A
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-2oZrfrJ7-1677154847152)(F:\2020Summer\暑期实习科研\实习\图片\121.JPG)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-6F14vxO6-1677154847153)(F:\2020Summer\暑期实习科研\实习\图片\212.JPG)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-b2IWX2FU-1677154847153)(F:\2020Summer\暑期实习科研\实习\图片\333.JPG)]
7.31
Cplex for Python – Day2
random
-
numpy.random.randint(lowb,upb,num) ==========随机生成指定范围内num个整数
-
random.randint(a,b) =========================随机生成指定范围内1个整数
-
================== 一次生成对应列表 N 里面的值的Key的键对 ================================
q = {i: np.random.randint(1, 10) for i in N}
- random.uniform()
用于生成一个指定范围内的随机浮点数，两个参数其中一个是上限，一个是下限，如果a>b则生成随机数
n,其中a<n<b。如果a<b则b<=n<=a - random.random()没有输入参数，生成0-1内的随机浮点数，输入参数 size 指定生成数组数目
matplotlib.pyplot
- 折线图
import matplotlib.pyplot as plt
#图形输入值
input_values = [1,2,3,4,5]
#图形输出值
squares = [1,4,9,16,25]
#plot根据列表绘制出有意义的图形，linewidth是图形线宽，可省略
plt.plot(input_values,squares,linewidth=5)
#设置图标标题
plt.title("Square Numbers",fontsize = 24)
#设置坐标轴标签
plt.xlabel("Value",fontsize = 14)
plt.ylabel("Square of Value",fontsize = 14)
#设置刻度标记的大小
plt.tick_params(axis='both',labelsize = 14)
#打开matplotlib查看器，并显示绘制图形
plt.show()
- 绘制一系列点
import matplotlib.pyplot as plt
x_values = [1,2,3,4,5]
y_values = [1,4,9,16,25]
plt.scatter(x_values[3:],y_values[3:],s=100)
# s指定绘制的散点的半径大小, [3:] x,y列表下标3及之后的数据
plt.show()
-
注释
plt.annotate('$q_{%d}=%d$' % (i, q[i]), (loc_x[i]+2, loc_y[i]))引号内是注释内容，第一个括号数据填充注释，第二个括号指定注释位置坐标
-
常见颜色参数
"colors"
============= ===============================
character color
============= ===============================
``'b'`` blue 蓝
``'g'`` green 绿
``'r'`` red 红
``'c'`` cyan 蓝绿
``'m'`` magenta 洋红
``'y'`` yellow 黄
``'k'`` black 黑
``'w'`` white 白
============= ===============================
'marker'
============= ===============================
character description
============= ===============================
``'.'`` point marker
``','`` pixel marker
``'o'`` circle marker
``'v'`` triangle_down marker
``'^'`` triangle_up marker
``'<'`` triangle_left marker
``'>'`` triangle_right marker
``'1'`` tri_down marker
``'2'`` tri_up marker
``'3'`` tri_left marker
``'4'`` tri_right marker
``'s'`` square marker
``'p'`` pentagon marker
``'*'`` star marker
``'h'`` hexagon1 marker
``'H'`` hexagon2 marker
``'+'`` plus marker
``'x'`` x marker
``'D'`` diamond marker
``'d'`` thin_diamond marker
``'|'`` vline marker
``'_'`` hline marker
============= ===============================
docplex
##########################################
# max 3 * x1 + 5 * x2 + 4 * x3
# s.t.
# 2 * x1 + 3 * x2 <= 1500
# 2 * x2 + 4 * x3 <= 800
# 3 * x1 + 2 * x2 + 5 * x3 <= 2000
# x1, x2, x3 >= 0
###########################################
from docplex.mp.model import Model #导出库，只用这一个就够了
model = Model() #创建模型
var_list = [i for i in range(0, 4)] #创建列表
X = model.continuous_var_list(var_list, lb=0, name='X') #创建变量列表
model.maximize(3 * X[0] + 5 * X[1] + 4 * X[2]) #设定目标函数
#添加约束条件
model.add_constraint(2 * X[0] + 3 * X[1] <= 1500)
model.add_constraint(2 * X[1] + 4 * X[2] <= 800)
model.add_constraint(3 * X[0] + 2 * X[1] + 5 * X[2] <= 2000)
sol = model.solve() #求解模型
print(sol) #打印结果
以下是用随机生成的数据作为输入来求解的代码
虽 老报错 如下：
但是最优解也解出来了！就很奇怪
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-9gnzAdIu-1677154847154)(F:\2020Summer\暑期实习科研\实习\图片\7311.JPG)]
from docplex.mp.model import Model
import numpy as np
import matplotlib.pyplot as plt
import json
n = 20
H = 1
V = [i for i in range(0, n)]
U = V
Q = np.random.randint(1, 20, n)
lng = np.random.random(n) * 100 # 0-100之间的随机数
lat = np.random.random(n) * 100 # 0-100之间的随机数
Z = [500 for i in V]
K = 0.001
C = [10 for j in V]
T = 100
plt.scatter(lng, lat, c='r')
for i in U:
plt.annotate('$q_{%d}=%d$' % (i, Q[i]), (lng[i]+2, lat[i]))
plt.show()
x = [(i, j) for i in U for j in V]
D = {(i, j): np.hypot(lng[i] - lng[j], lat[i] - lat[j]) for i, j in x} # 距离矩阵
# 创建模型
mdl = Model('model')
# 添加决策变量
x = mdl.binary_var_dict(x, name='x') # 创建binary变量列表
y = mdl.binary_var_dict(V, name='y') # 创建连续变量列表
# 约束函数
for i in U:
for j in V:
mdl.add_constraint(T >= D[i, j]*x[1, j])
mdl.add_constraints(mdl.sum(x[i, j] for j in V) == 1 for i in U)
mdl.add_constraints(Z[j] >= mdl.sum(Q[j]*x[i, j] for i in U) for j in V)
# mdl.add_constraints(K*Z[j] <= mdl.sum(x[i, j]*Q[j] for i in U) for j in V)
# 这个约束一直一直加不进去呢？？？？
mdl.add_indicator_constraints(mdl.indicator_constraint(x[i, j], y[j] == 1) for i in U for j in V)
# 目标函数
mdl.minimize(mdl.sum(H*D[i, j]*x[i, j] for j in V for i in U) + mdl.sum(C[j] * y[j] for j in V))
solution = mdl.solve(log_output=True)
print(solution)
print(solution.solve_details)
导入实际数据
数据容量有限，20个以上可能就会报错
import json
import math
import matplotlib.pyplot as plt
from docplex.mp.model import Model
# ===================================================================================
# 通过经纬度测算距离
# INPUTS: A点经度 A点纬度 B点经度 B点纬度
# OUTPUTS: 返回值为考虑球面的两点直线距离,单位：km
# ===================================================================================
def distance(lng_A, lat_A, lng_B, lat_B):
# 地球半径
R = 6371.004
pi = 3.141592654
Mlng_A = lng_A
Mlat_A = 90 - lat_A
Mlng_B = lng_B
Mlat_B = 90 - lat_B
C = math.sin(Mlat_A * pi / 180) * math.sin(Mlat_B * pi / 180) * math.cos((Mlng_A - Mlng_B) * pi / 180) + math.cos(
Mlat_A * pi / 180) * math.cos(Mlat_B * pi / 180)
if C >= 1:
C = 1
Distance = R * math.acos(C)
return Distance
# ============================================================================
filename = (r'19_10.json')
data = json.load(open(filename))
date = '8'
if date not in data:
print("%s 号当日无订货记录。" % date)
exit(0)
# 把data['x']的key值转化为列表输出
key = list(data[date])
data_date = data[date]
n = 8
key = key[50:n + 50]
# print(key)
lng = []
lat = []
amount = []
for i in range(0, n):
lng.append(data_date[key[i]]['lng'])
lat.append(data_date[key[i]]['lat'])
amount.append(data_date[key[i]]['amount'])
V = [i for i in range(0, n)]
U = V
Z = [300 for i in V]
K = 0.001
C = [10 for j in V]
T = 100
H = 1
plt.scatter(lng, lat, c='r')
for i in V:
plt.annotate('$q_{%d}=%d$' % (i, amount[i]), (lng[i]+2, lat[i]))
plt.show()
x = [(i, j) for i in U for j in V]
D = {(i, j): distance(lng[i], lat[i], lng[j], lat[j]) for i, j in x} # 距离矩阵
# 创建模型
mdl = Model('model')
# 添加决策变量
x = mdl.binary_var_dict(x, name='x') # 创建binary变量列表
y = mdl.binary_var_dict(V, name='y') # 创建binary变量列表
# 约束函数
for i in U:
for j in V:
mdl.add_constraint(T >= D[i, j] * x[1, j])
mdl.add_constraints(mdl.sum(x[i, j] for j in V) == 1 for i in U)
mdl.add_constraints(Z[j] >= mdl.sum(amount[i] * x[i, j] for i in U) for j in V)
# mdl.add_indicator_constraints(
# mdl.indicator_constraint(y[j] == 1, K*Z[j] <= mdl.sum(x[i, j]*amount[i] for i in U)) for j in V
# )
mdl.add_indicator_constraints(mdl.indicator_constraint(x[i, j], y[j] == 1) for i in U for j in V)
# 目标函数
mdl.minimize(mdl.sum(H * D[i, j] * x[i, j] for j in V for i in U) + mdl.sum(C[j] * y[j] for j in V))
solution = mdl.solve(log_output=True)
print("solution: ", solution)
print("solve_details: ", solution.solve_details)
8.6
8.10
将按照星期几分类存放为json文件
import json
import matplotlib.pyplot as plt
# =============================================================================================
def WeekData(month, day, data):
Tues = {}
Wed = {}
Thur = {}
Fri = {}
Sat = {}
Sun = {}
Mon = {}
if month == '10':
if day % 7 == 0:
Mon.update(data)
elif (day - 1) % 7 == 0:
Tues.update(data)
elif (day - 2) % 7 == 0:
Wed.update(data)
elif (day - 3) % 7 == 0:
Thur.update(data)
elif (day - 4) % 7 == 0:
Fri.update(data)
elif (day - 5) % 7 == 0:
Sat.update(data)
elif (day - 6) % 7 == 0:
Sun.update(data)
elif month == '11':
if day % 7 == 0:
Thur.update(data)
elif (day - 1) % 7 == 0:
Fri.update(data)
elif (day - 2) % 7 == 0:
Sat.update(data)
elif (day - 3) % 7 == 0:
Sun.update(data)
elif (day - 4) % 7 == 0:
Mon.update(data)
elif (day - 5) % 7 == 0:
Tues.update(data)
elif (day - 6) % 7 == 0:
Wed.update(data)
elif month == '12':
if day % 7 == 0:
Sat.update(data)
elif (day - 1) % 7 == 0:
Sun.update(data)
elif (day - 2) % 7 == 0:
Mon.update(data)
elif (day - 3) % 7 == 0:
Tues.update(data)
elif (day - 4) % 7 == 0:
Wed.update(data)
elif (day - 5) % 7 == 0:
Thur.update(data)
elif (day - 6) % 7 == 0:
Fri.update(data)
return Mon, Tues, Wed, Thur, Fri, Sat, Sun
# ===============================================================================================
filename = (r'19_10.json')
data10 = json.load(open(filename))
filename = (r'19_11.json')
data11 = json.load(open(filename))
filename = (r'19_12.json')
data12 = json.load(open(filename))
MonTotal = {}
TuesTotal = {}
WedTotal = {}
ThurTotal = {}
FriTotal = {}
SatTotal = {}
SunTotal = {}
Tues = {}
Wed = {}
Thur = {}
Fri = {}
Sat = {}
Sun = {}
Mon = {}
for i in range(1, 32):
day = str(i)
if day not in data10:
print("19年10月 %s 号当日无订货记录。" % day)
continue
Mon, Tues, Wed, Thur, Fri, Sat, Sun = WeekData('10', i, data10[day])
MonTotal.update(Mon)
TuesTotal.update(Tues)
WedTotal.update(Wed)
ThurTotal.update(Thur)
FriTotal.update(Fri)
SatTotal.update(Sat)
SunTotal.update(Sun)
for j in range(1, 32):
day = str(j)
if day not in data11:
print("19年11月 %s 号当日无订货记录。" % day)
continue
Mon, Tues, Wed, Thur, Fri, Sat, Sun = WeekData('11', j, data11[day])
MonTotal.update(Mon)
TuesTotal.update(Tues)
WedTotal.update(Wed)
ThurTotal.update(Thur)
FriTotal.update(Fri)
SatTotal.update(Sat)
SunTotal.update(Sun)
for k in range(1, 32):
day = str(k)
if day not in data12:
print("19年12月 %s 号当日无订货记录。" % day)
continue
Mon, Tues, Wed, Thur, Fri, Sat, Sun = WeekData('12', k, data12[day])
MonTotal.update(Mon)
TuesTotal.update(Tues)
WedTotal.update(Wed)
ThurTotal.update(Thur)
FriTotal.update(Fri)
SatTotal.update(Sat)
SunTotal.update(Sun)
# ==========================计算一周内每天的订货量总量并画个折线图=========================
amount = [0, 0, 0, 0, 0, 0, 0]
for num in range(1, 20000000):
ID = str(num)
if ID in MonTotal:
amount[0] += MonTotal[ID]['amount']
# print('1 : %d' % amount[0])
if ID in TuesTotal:
amount[1] += TuesTotal[ID]['amount']
if ID in WedTotal:
amount[2] += WedTotal[ID]['amount']
if ID in ThurTotal:
amount[3] += ThurTotal[ID]['amount']
if ID in FriTotal:
amount[4] += FriTotal[ID]['amount']
if ID in SatTotal:
amount[5] += SatTotal[ID]['amount']
if ID in SunTotal:
amount[6] += SunTotal[ID]['amount']
# print('7 : %d' % amount[6])
plt.plot(range(len(amount)), amount)
plt.show()
8.11
Gurobi建模基本概念：
-
Parameter （参数）控制优化器的行为， 需要在优化启动前设置。
-
Attributes （属性）控制模型（包括模型 、变量、约束、目标等对象）的特性。
-
TupleList和TupleDict，可以筛选对象来组成约束和变量
-
quicksum 相当于 ∑ \sum ∑
-
x.prod(cost) 变量和系数相乘后累加，x与cost有相同的键值
建立7个场景
它数据多了，MemoryError
import matplotlib.pyplot as plt
import numpy as np
import gurobipy as gp
from gurobipy import GRB
import json
import math
# ===================================================================================
# 通过经纬度测算距离
# INPUTS: A点经度 A点纬度 B点经度 B点纬度
# OUTPUTS: 返回值为考虑球面的两点直线距离,单位：km
# ===================================================================================
def distance(lng_A, lat_A, lng_B, lat_B):
# 地球半径
R = 6371.004
pi = 3.141592654
Mlng_A = lng_A
Mlat_A = 90 - lat_A
Mlng_B = lng_B
Mlat_B = 90 - lat_B
C = math.sin(Mlat_A * pi / 180) * math.sin(Mlat_B * pi / 180) * math.cos((Mlng_A - Mlng_B) * pi / 180) + math.cos(
Mlat_A * pi / 180) * math.cos(Mlat_B * pi / 180)
if C >= 1:
C = 1
Distance = R * math.acos(C)
return Distance
# =======================================================================================
filename1 = (r'MonList.json')
Mon = json.load(open(filename1))
filename2 = (r'TuesList.json')
Tues = json.load(open(filename2))
filename3 = (r'WedList.json')
Wed = json.load(open(filename3))
filename4 = (r'ThurList.json')
Thur = json.load(open(filename4))
filename5 = (r'FriList.json')
Fri = json.load(open(filename5))
filename6 = (r'SatList.json')
Sat = json.load(open(filename6))
filename7 = (r'SunList.json')
Sun = json.load(open(filename7))
U1, amount1, lng1, lat1 = gp.multidict(Mon)
U2, amount2, lng2, lat2 = gp.multidict(Tues)
U3, amount3, lng3, lat3 = gp.multidict(Wed)
U4, amount4, lng4, lat4 = gp.multidict(Thur)
U5, amount5, lng5, lat5 = gp.multidict(Fri)
U6, amount6, lng6, lat6 = gp.multidict(Sat)
U7, amount7, lng7, lat7 = gp.multidict(Sun)
D1 = {
(i, j): distance(lng1[i], lat1[i], lng1[j], lat1[j]) for i in U1 for j in U1 if i != j
}
D2 = {
(i, j): distance(lng2[i], lat2[i], lng2[j], lat2[j]) for i in U2 for j in U2 if i != j
}
D3 = {
(i, j): distance(lng3[i], lat3[i], lng3[j], lat3[j]) for i in U3 for j in U3 if i != j
}
D4 = {
(i, j): distance(lng4[i], lat4[i], lng4[j], lat4[j]) for i in U4 for j in U4 if i != j
}
D5 = {
(i, j): distance(lng5[i], lat5[i], lng5[j], lat5[j]) for i in U5 for j in U5 if i != j
}
D6 = {
(i, j): distance(lng6[i], lat6[i], lng6[j], lat6[j]) for i in U6 for j in U6 if i != j
}
D7 = {
(i, j): distance(lng7[i], lat7[i], lng7[j], lat7[j]) for i in U7 for j in U7 if i != j
}
for i in U1:
plt.scatter(lng1[i], lat1[i], c='r')
plt.annotate('$q_{%s}=%d$' % (i, amount1[i]), (lng1[i]+2, lat1[i]))
# plt.plot(loc_x[0], loc_y[0], c='b', marker='s')
plt.show()
gurobi里面facility.py的示例
# Facility location: a company currently ships its product from 5 plants
# to 4 warehouses. It is considering closing some plants to reduce
# costs. What plant(s) should the company close, in order to minimize
# transportation and fixed costs?
#
# Note that this example uses lists instead of dictionaries. Since
# it does not work with sparse data, lists are a reasonable option.
#
# Based on an example from Frontline Systems:
# http://www.solver.com/disfacility.htm
# Used with permission.
import gurobipy as gp
from gurobipy import GRB
# Warehouse demand in thousands of units
demand = [15, 18, 14, 20]
# Plant capacity in thousands of units
capacity = [20, 22, 17, 19, 18]
# Fixed costs for each plant
fixedCosts = [12000, 15000, 17000, 13000, 16000]
# Transportation costs per thousand units
transCosts = [[4000, 2000, 3000, 2500, 4500],
[2500, 2600, 3400, 3000, 4000],
[1200, 1800, 2600, 4100, 3000],
[2200, 2600, 3100, 3700, 3200]]
# Range of plants and warehouses
plants = range(len(capacity))
warehouses = range(len(demand))
# Model
m = gp.Model("facility")
# Plant open decision variables: open[p] == 1 if plant p is open.
open = m.addVars(plants,
vtype=GRB.BINARY,
obj=fixedCosts,
name="open")
# Transportation decision variables: transport[w,p] captures the
# optimal quantity to transport to warehouse w from plant p
transport = m.addVars(warehouses, plants, obj=transCosts, name="trans")
# You could use Python looping constructs and m.addVar() to create
# these decision variables instead. The following would be equivalent
# to the preceding two statements...
#
# open = []
# for p in plants:
# open.append(m.addVar(vtype=GRB.BINARY,
# obj=fixedCosts[p],
# name="open[%d]" % p))
#
# transport = []
# for w in warehouses:
# transport.append([])
# for p in plants:
# transport[w].append(m.addVar(obj=transCosts[w][p],
# name="trans[%d,%d]" % (w, p)))
# The objective is to minimize the total fixed and variable costs
m.modelSense = GRB.MINIMIZE
# Production constraints
# Note that the right-hand limit sets the production to zero if the plant
# is closed
m.addConstrs(
(transport.sum('*', p) <= capacity[p]*open[p] for p in plants), "Capacity")
# Using Python looping constructs, the preceding would be...
#
# for p in plants:
# m.addConstr(sum(transport[w][p] for w in warehouses)
# <= capacity[p] * open[p], "Capacity[%d]" % p)
# Demand constraints
m.addConstrs(
(transport.sum(w) == demand[w] for w in warehouses),
"Demand")
# ... and the preceding would be ...
# for w in warehouses:
# m.addConstr(sum(transport[w][p] for p in plants) == demand[w],
# "Demand[%d]" % w)
# Save model
m.write('facilityPY.lp')
# Guess at the starting point: close the plant with the highest fixed costs;
# open all others
# First open all plants
for p in plants:
open[p].start = 1.0
# Now close the plant with the highest fixed cost
print('Initial guess:')
maxFixed = max(fixedCosts)
for p in plants:
if fixedCosts[p] == maxFixed:
open[p].start = 0.0
print('Closing plant %s' % p)
break
print('')
# Use barrier to solve root relaxation
m.Params.method = 2
# Solve
m.optimize()
# Print solution
print('\nTOTAL COSTS: %g' % m.objVal)
print('SOLUTION:')
for p in plants:
if open[p].x > 0.99:
print('Plant %s open' % p)
for w in warehouses:
if transport[w, p].x > 0:
print(' Transport %g units to warehouse %s' %
(transport[w, p].x, w))
else:
print('Plant %s closed!' % p)
创建一个指示变量约束
MODEL.addGenConstrIndicator(binvar, binval, expression, name=“”)
指示变量 binvar 的值取 binval 时, 进行约束 expression
案例:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-gv2VebTW-1677154849773)(null)]
方法一: 构造指示变量 , 则上述约束转化为
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-lnTGuhgu-1677154849475)(null)]
方法二: 转为二次约束, 若矩阵非正定矩阵, 则无法求解
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-5X8qdcCl-1677154849188)(null)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-fFedjzID-1677154847155)(F:\2020Summer\暑期实习科研\实习\图片\4.JPG)]
8.13
基于python遗传算法的选址问题
更多详细内容~参见这篇博客——
https://blog.csdn.net/Jarrodche/article/details/93136237?ops_request_misc=%7B%22request%5Fid%22%3A%22159730579019724839261750%22%2C%22scm%22%3A%2220140713.130102334.pc%5Fall.%22%7D&request_id=159730579019724839261750&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2allfirst_rank_ecpm_v3~pc_rank_v2-1-93136237.first_rank_ecpm_v3_pc_rank_v2&utm_term=python+遗传算法+选址&spm=1018.2118.3001.4187)
模型如下：
给定某一地区所有需求点的地址集合，要求从备选地址中选出一定数目的地址建立工厂，从而建立一系列的配送区域，实现各个需求点的配送，使得在选出点建立的工厂与各需求点形成的配送系统总配送费用最少。为了便于建立模型，需作一定的假设，假设系统满足如下一些条件：
-
仅在一定的被选对象内考虑新的工厂建厂；
-
运输费用与运输量成正比；
-
需求点的需求按区域总计；
-
工厂生产能力可以满足需求；
-
各需求点的需求量一定且为己知；
-
总费用只考虑固定的工厂建址费用和运输费用。
该选址问题的混合整数规划模型为：
式中，m为工厂备选地点的个数；n为需求点的个数；p为可设置工厂的最大个数； x i j x_{ij} xij为从工厂 i 到需求点 j 的运输量； y i y_i yi为整数变量，当 y i = 1 y_i=1 yi=1 时表示 i 地被选中，当时表示地未被 y i = 0 y_i=0 yi=0 选中； b j b_j bj为第 j 地的需求量； a i a_i ai为被选工厂 i 的建设容量（供应能力）； c i j c_{ij} cij为从 i 地到 j 地的包括装卸运输费在内的单位运输费用； d i d_i di为被选工厂 i 的固定费用（包括基本投资费和固定经营费）。
目标函数(1)中，等式右边第一项是从工厂到需求点的运输成本；第二项是工厂设置的固定成本。式(2)表示从各工厂向某需求点供给的物资总和应满足该需求点的需求量；式(3)表示如果i地点被选中建厂，则从它发出的货物总量不超过它的建设容量（供应能力）；式(4)表示选上的工厂个数不得超过原定的最大限额。
import geatpy as ga # 解遗传算法的库
from pymprog import * # 解运输问题的库
import numpy as np
import xlrd
# -------------------解码函数---------------------------
def decode(Y, A, B, m):
# step1
List1 = np.multiply(np.array(Y), np.array(A)) # 计算 y*a
sum1 = 0 # 再把他们累加起来
for i in List1:
sum1 += i
sum2 = 0
for j in B:
sum2 += j
# step2
while sum1 < sum2:
k = np.random.randint(0, m)
if Y[k] == 0:
sum1 += A[k]
Y[k] = 1
else:
k = np.random.randint(1, m + 1)
continue
# else:
# print("算法终止")
return Y
# -----------------------适值函数-------------------------------
# 固定成本函数
def Fixed_cost(List2):
fixed_cost = 0 # 再进行累加计算 计算建厂的固定成本
for i in List2:
fixed_cost += i
return fixed_cost
# 运输成本函数
def report(Y_decode, A, B): # 返回(产地,销地):运价；产地:产量； 销地:销量
# 把目标函数简化成产大于销的运输问题
a = ("A1", "A2", "A3", "A4", "A5", "A6", "A7", "A8", "A9", "A10") # 工厂（产地）
b = ("B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8", "B9", "B10",
"B11", "B12", "B13", "B14", "B15", "B16", "B17", "B18", "B19", "B20", "B21") # 需求地（销地）
datas = dict() # 定义一个空字典 (产地,销地):运价
datac = dict() # 定义一个空字典 产地:产量
datax = dict() # 定义一个空字典 销地:销量
List3 = np.multiply(np.array(Y_decode), np.array(A)) # 计算 y*a 产量
# （产地，销地）：运价
Data = xlrd.open_workbook('单位运输费用.xlsx') # 读取Excel运价列表
table = Data.sheets()[0] # 读取Excel列表的sheet1
for i in range(0, 10): # 先从十行里面拿出一行来
x = np.array(table.row_values(i))
for j in range(0, 21): # 再从21列中取一列
data = x[j]
datas[a[i], b[j]] = data # (产地,销地):运价
# 产地：产量
y = List3
for i in range(0, 10):
datac[a[i]] = y[i] # 产地:产量
# 销地：销量
for j in range(0, 21):
datax[b[j]] = B[j] # 销地:销量
return (datas, datac, datax)
def Transport_cost(datas, datac, datax):
a = ("A1", "A2", "A3", "A4", "A5", "A6", "A7", "A8", "A9", "A10") # 工厂（产地）
b = ("B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8", "B9", "B10",
"B11", "B12", "B13", "B14", "B15", "B16", "B17", "B18", "B19", "B20", "B21") # 需求地（销地）
# 开始模型计算
X = np.zeros((10, 21)) # 生成10行21列的0矩阵
cost = 0.0
begin('Transport')
x = var('x', datas.keys()) # 调运方案
minimize(sum(datas[i, j] * x[i, j] for (i, j) in datas.keys()), 'Cost') # 总运费最少
for i in datac.keys(): # 产地产量约束
sum(x[i, j] for j in datax.keys()) == datac[i]
for j in datax.keys(): # 销地销量约束
sum(x[i, j] for i in datac.keys()) == datax[j]
solve()
for (i, j) in datas.keys():
# print('x[i, j].primal',x[i, j].primal)
if x[i, j].primal > 0 and datas[i, j] != 0:
# print("产地:%s -> 销地:%s 运输量:%-2d 运价:%2d" % (i, j, int(x[i, j].primal), int(datas[i, j])))
X[a.index(i)][b.index(j)] = int(x[i, j].primal)
cost = int(vobj())
# print("X", X)
# print("总费用:%d" % int(vobj()))
end()
return (X, cost)
# -----------------------主函数开始-------------------------------------
def main():
# 需求点的需求量
B = [6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 0]
# 工厂容量（生产能力）
A = [30, 40, 50, 55, 45, 48, 36, 50, 60, 35]
# 各待选工厂的固定费用
D = [65, 50, 40, 30, 20, 50, 45, 30, 35, 25]
# 从待选工厂到需求点的单位运输费用
Data = xlrd.open_workbook('单位运输费用.xlsx')
table = Data.sheets()[0]
Y = [] # 定义变量码
Y_decode = [] # 定义解码生成的变量码
m = 10 # 工厂个数
T = 100 # 迭代次数
decode_after_chrom2 = np.zeros((60, 10)) # 存储60条Y_decode
datas = dict() # 定义一个空字典 (产地,销地):运价
datac = dict() # 定义一个空字典 产地:产量
datax = dict() # 定义一个空字典 销地:销量
fixed_cost = 0.0 # 初始固定成本
X = np.zeros((10, 21)) # 运量矩阵
transport_cost = 0.0 # 初始运输费用
cost = 0.0 # 初始总费用
best_y = np.zeros((1, 10)) # 所有代中的最优选址方案
min_t = 0 # 最优的方案出现的代数
all_min_cost = 1000 # 所有代中最小总费用
all_min_X = np.zeros((10, 21)) # 所有代中最小总费用的运量矩阵
FitnV = np.zeros((1, 60)) # 适应度列表
NewChrlx = [] # 轮盘赌存储位置
New_chrom1 = np.zeros((60, 10))
New_chrom2 = np.zeros((60, 10))
# ---------------------遗传编码---------------------------
# 生成附加码
chrom1 = ga.crtpp(60, 10, 10) # 种群大小，染色体长度，种群最大元素
print(chrom1)
# 生成变量码
FieldDR = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])
chrom2 = ga.crtip(60, FieldDR) # 种群大小，上界下界区域描述器
print(chrom2)
add_code = np.zeros((1, 10)) # 定义父代单行附加码
new_add_code = np.zeros((1, 10)) # 定义子代单行附加码
variable_code = np.zeros((1, 10)) # 定义父代单行变量码
new_variable_code = np.zeros((1, 10)) # 定义子代单行变量码
index_v = 0 # 定义初始索引值
begin_New_chrom2 = np.zeros((60, 10)) # 定义经过交叉后的新的变量码种群
again_New_chrom1 = np.zeros((60, 10)) # 经过变异后的附加码种群
# tem_list = []
a = 0
for t in np.arange(0, T):
min_cost = 1000 # 初始最小总费用
min_X = np.zeros((10, 21)) # 初始最小总费用的运量矩阵
min_i = 400 # 每一代中最小的运价在种群对应的位置
better_y = np.zeros((1, 10)) # 每一代中的最有选址方案
print("==========================第%d代种群======================" % t)
for i in np.arange(0, 60):
# print("第%d个染色体" % i)
Y = chrom2[i]
# print("-------------%d-----------------" % i)
# print("Y",Y)
Y_decode = np.array(decode(Y, A, B, m))
decode_after_chrom2[i] = Y_decode
List2 = np.multiply(np.array(D), np.array(Y_decode)) # 计算 d*y
fixed_cost = Fixed_cost(List2)
print("fixed_cost", fixed_cost)
(datas, datac, datax) = report(Y_decode, A, B)
(X, transport_cost) = Transport_cost(datas, datac, datax)
cost = fixed_cost + transport_cost
print('第%d个染色体的cost' % i, cost)
FitnV[0][i] = 1 / cost
if cost < min_cost:
min_cost = cost
min_X = X
min_i = i
better_y = Y_decode
# print("第%d代种群中最小的cost"%t, min_cost)
# print("最优运输方案", min_X)
# print('对应的位置', min_i)
NewChrlx = ga.rws(np.array(FitnV).T, 59) # 轮盘赌
# print(NewChrlx)
for i in NewChrlx:
New_chrom1[i] = chrom1[i]
New_chrom2[i] = decode_after_chrom2[i]
for j in range(0, 60): # 把最大的适值得到染色体拿到新的种群中
if j not in NewChrlx:
New_chrom1[j] = chrom1[j]
New_chrom2[j] = decode_after_chrom2[j]
# 遗传操作
Pc = 0.6 # 交叉概率
Pm = 0.1 # 变异概率
New_chrom1 = ga.xovpm(New_chrom1, 0.6)
# print("交叉后的附加码New_chrom1",New_chrom1)
# print("chrom1",chrom1)
for change_1 in np.arange(0, 60):
# print('------------修改chrom2第%d次-----------'%change_1)
new_add_code = New_chrom1[change_1]
add_code = chrom1[change_1]
# print(type(add_code))
variable_code = New_chrom2[change_1]
for change_2 in np.arange(0, 10):
index_v = add_code.tolist().index(new_add_code[change_2])
# print("index_v",index_v)
new_variable_code[0, change_2] = variable_code[index_v]
begin_New_chrom2[change_1] = new_variable_code
# print("begin_New_chrom2[change_1]",begin_New_chrom2[change_1])
new_variable_code = np.zeros((1, 10))
# print("交叉后的变量码begin_New_chrom2",begin_New_chrom2)
# 变异操作
ret = np.random.random()
if ret > Pm:
for mutation_v in np.arange(0, 60):
new_add_code = New_chrom1[mutation_v]
# print('未翻转的附加码new_add_code',new_add_code)
random_1 = np.random.randint(0, 10)
random_2 = np.random.randint(0, 10)
a = np.abs(random_1 - random_2)
tem_list = np.zeros(a + 1)
if a != 0:
if random_2 > random_1:
for i, k in zip(np.arange(0, np.abs(random_2 - random_1) + 1), np.arange(random_1, random_2 + 1)):
tem_list[i] = new_add_code[k]
# print("未翻转的单行tem_list",tem_list)
# print("未翻转的tem_list类型",type(tem_list))
tem_list = tem_list.tolist()
# print("转换为列表的tem_list", tem_list)
tem_list.reverse()
# print("翻转后的tem_list",tem_list)
for i, k in zip(np.arange(0, np.abs(random_2 - random_1) + 1).tolist(),
np.arange(random_1, random_2 + 1).tolist()):
new_add_code[k] = tem_list[i]
# print("翻转后的附加码",new_add_code)
again_New_chrom1[mutation_v] = new_add_code
else:
for i, k in zip(np.arange(0, np.abs(random_2 - random_1) + 1), np.arange(random_2, random_1 + 1)):
tem_list[i] = new_add_code[k]
# print("未翻转的单行tem_list",tem_list)
# print("未翻转的tem_list类型",type(tem_list))
tem_list = tem_list.tolist()
# print("转换为列表的tem_list",tem_list)
tem_list.reverse()
# print("tem_list",tem_list)
for i, k in zip(np.arange(0, np.abs(random_2 - random_1) + 1).tolist(),
np.arange(random_2, random_1 + 1).tolist()):
new_add_code[k] = tem_list[i]
again_New_chrom1[mutation_v] = new_add_code
# print("循环中的again_New_chrom1", again_New_chrom1)
else:
again_New_chrom1[mutation_v] = new_add_code
# print("else循环中的again_New_chrom1", again_New_chrom1)
# print("New_chrom1",New_chrom1)
# print('变异后的附加码again_New_chrom1',again_New_chrom1)
chrom1 = again_New_chrom1
chrom2 = begin_New_chrom2
print('遗传操作后的chrom2', chrom2)
print("第%d代种群中最小的cost" % t, min_cost)
print("最优运输方案", min_X)
print('对应的位置', min_i)
if all_min_cost > min_cost:
all_min_cost = min_cost
all_min_X = min_X
best_y = better_y
min_t = t
# New_chrom1[60]= chrom1[min_i]
# New_chrom2[60] = decode_after_chrom2[min_i]
# print("FitnV",np.array(FitnV).T.shape)
print("所有代中种群中最小的cost", all_min_cost)
print("所有代中的最优运输方案", all_min_X)
print("最优的选址方案", best_y)
print("最优的选址方案出现在第%d代" % min_t)
if __name__ == "__main__":
main()
场景建模 Gurobi
import matplotlib.pyplot as plt
import gurobipy as gp
from gurobipy import *
import json
import math
# ===================================================================================
# 通过经纬度测算距离
# INPUTS: A点经度 A点纬度 B点经度 B点纬度
# OUTPUTS: 返回值为考虑球面的两点直线距离,单位：km
# ===================================================================================
def distance(lng_A, lat_A, lng_B, lat_B):
# 地球半径
R = 6371.004
pi = 3.141592654
Mlng_A = lng_A
Mlat_A = 90 - lat_A
Mlng_B = lng_B
Mlat_B = 90 - lat_B
C = math.sin(Mlat_A * pi / 180) * math.sin(Mlat_B * pi / 180) * math.cos((Mlng_A - Mlng_B) * pi / 180) + math.cos(
Mlat_A * pi / 180) * math.cos(Mlat_B * pi / 180)
if C >= 1:
C = 1
Distance = R * math.acos(C)
return Distance
# =================================================================================
# 将以str作为key值的字典转化为以int为key值的字典,并只保留前n个key-value对
# 输入：dict_str, n
# 输出：dict_int
# ===================================================================================
def Key_StrToInt(dict_str, n):
dict_int = {}
Keys = dict_str.keys()
Keys = list(Keys)
Keys = Keys[0:n]
for i in Keys:
i = int(i)
temp = {
i: dict_str[str(i)]
}
dict_int.update(temp)
return dict_int
# ==================读入以str为key值的字典=====================
filename1 = (r'MonList.json')
Mon_str = json.load(open(filename1))
filename2 = (r'TuesList.json')
Tues_str = json.load(open(filename2))
filename3 = (r'WedList.json')
Wed_str = json.load(open(filename3))
filename4 = (r'ThurList.json')
Thur_str = json.load(open(filename4))
filename5 = (r'FriList.json')
Fri_str = json.load(open(filename5))
filename6 = (r'SatList.json')
Sat_str = json.load(open(filename6))
filename7 = (r'SunList.json')
Sun_str = json.load(open(filename7))
# ===================================================================================
n = 100 # 定义取每个场景的多少个数据
# 把字典转化成以int数字为key值的字典
Mon_int = Key_StrToInt(Mon_str, n)
Tues_int = Key_StrToInt(Tues_str, n)
Wed_int = Key_StrToInt(Wed_str, n)
Thur_int = Key_StrToInt(Thur_str, n)
Fri_int = Key_StrToInt(Fri_str, n)
Sat_int = Key_StrToInt(Sat_str, n)
Sun_int = Key_StrToInt(Sun_str, n)
U1, amount1, lng1, lat1 = gp.multidict(Mon_int)
U2, amount2, lng2, lat2 = gp.multidict(Tues_int)
U3, amount3, lng3, lat3 = gp.multidict(Wed_int)
U4, amount4, lng4, lat4 = gp.multidict(Thur_int)
U5, amount5, lng5, lat5 = gp.multidict(Fri_int)
U6, amount6, lng6, lat6 = gp.multidict(Sat_int)
U7, amount7, lng7, lat7 = gp.multidict(Sun_int)
capacity = 3000
FixedCost = 100
T = 300
# 求两个集合的并集
U1 = list(U1)
U2 = list(U2)
U3 = list(U3)
U4 = list(U4)
U5 = list(U5)
U6 = list(U6)
U7 = list(U7)
V = list(set(U1).union(U2, U3, U4, U5, U6, U7))
# T = [T for h in range(len(V))]
Z = [capacity for i in V]
C = [FixedCost for c in V]
# 把并集之后的数据写入lng[],lat[]，作为仓库潜在选址集合坐标
lng = {}
lat = {}
for k in V:
if k in U1:
a = lng1[k]
b = lat1[k]
elif k in U2:
a = lng2[k]
b = lat2[k]
elif k in U3:
a = lng3[k]
b = lat3[k]
elif k in U4:
a = lng4[k]
b = lat4[k]
elif k in U5:
a = lng5[k]
b = lat5[k]
elif k in U6:
a = lng6[k]
b = lat6[k]
elif k in U7:
a = lng7[k]
b = lat7[k]
x = {
k: a
}
y = {
k: b
}
lng.update(x)
lat.update(y)
# 计算各场景下距离矩阵D_jk
D1 = {
(i, j): distance(lng1[i], lat1[i], lng[j], lat[j]) for i in U1 for j in V
}
D2 = {
(i, j): distance(lng2[i], lat2[i], lng[j], lat[j]) for i in U2 for j in V
}
D3 = {
(i, j): distance(lng3[i], lat3[i], lng[j], lat[j]) for i in U3 for j in V
}
D4 = {
(i, j): distance(lng4[i], lat4[i], lng[j], lat[j]) for i in U4 for j in V
}
D5 = {
(i, j): distance(lng5[i], lat5[i], lng[j], lat[j]) for i in U5 for j in V
}
D6 = {
(i, j): distance(lng6[i], lat6[i], lng[j], lat[j]) for i in U6 for j in V
}
D7 = {
(i, j): distance(lng7[i], lat7[i], lng[j], lat[j]) for i in U7 for j in V
}
# 创建解集下标
x1 = [(i, j) for i in U1 for j in V]
x2 = [(i, j) for i in U2 for j in V]
x3 = [(i, j) for i in U3 for j in V]
x4 = [(i, j) for i in U4 for j in V]
x5 = [(i, j) for i in U5 for j in V]
x6 = [(i, j) for i in U6 for j in V]
x7 = [(i, j) for i in U7 for j in V]
x1 = tuplelist(x1)
x2 = tuplelist(x2)
x3 = tuplelist(x3)
x4 = tuplelist(x4)
x5 = tuplelist(x5)
x6 = tuplelist(x6)
x7 = tuplelist(x7)
V = tuplelist(V)
# Model
SM = gp.Model("Scenario")
# 添加变量
y = SM.addVars(V,
vtype=GRB.BINARY,
obj=C,
name="y")
x1 = SM.addVars(x1, vtype=GRB.BINARY, name="x1")
x2 = SM.addVars(x2, vtype=GRB.BINARY, name="x2")
x3 = SM.addVars(x3, vtype=GRB.BINARY, name="x3")
x4 = SM.addVars(x4, vtype=GRB.BINARY, name="x4")
x5 = SM.addVars(x5, vtype=GRB.BINARY, name="x5")
x6 = SM.addVars(x6, vtype=GRB.BINARY, name="x6")
x7 = SM.addVars(x7, vtype=GRB.BINARY, name="x7")
SM.update()
# 设置目标函数
SM.setObjective(
sum(y[j] * FixedCost for j in V)
+ sum(x1[i, j] * D1[i, j] for i in U1 for j in V)
+ sum(x2[i, j] * D2[i, j] for i in U2 for j in V)
+ sum(x3[i, j] * D3[i, j] for i in U3 for j in V)
+ sum(x4[i, j] * D4[i, j] for i in U4 for j in V)
+ sum(x5[i, j] * D5[i, j] for i in U5 for j in V)
+ sum(x6[i, j] * D6[i, j] for i in U6 for j in V)
+ sum(x7[i, j] * D7[i, j] for i in U7 for j in V),
GRB.MINIMIZE)
# 添加约束
SM.addConstrs(T >= D1[i, j] * x1[i, j] for i in U1 for j in V)
SM.addConstrs(T >= D2[i, j] * x2[i, j] for i in U2 for j in V)
SM.addConstrs(T >= D3[i, j] * x3[i, j] for i in U3 for j in V)
SM.addConstrs(T >= D4[i, j] * x4[i, j] for i in U4 for j in V)
SM.addConstrs(T >= D5[i, j] * x5[i, j] for i in U5 for j in V)
SM.addConstrs(T >= D6[i, j] * x6[i, j] for i in U6 for j in V)
SM.addConstrs(T >= D7[i, j] * x7[i, j] for i in U7 for j in V)
SM.addConstrs(sum(x1[i, j] for j in V) == 1 for i in U1)
SM.addConstrs(sum(x2[i, j] for j in V) == 1 for i in U2)
SM.addConstrs(sum(x3[i, j] for j in V) == 1 for i in U3)
SM.addConstrs(sum(x4[i, j] for j in V) == 1 for i in U4)
SM.addConstrs(sum(x5[i, j] for j in V) == 1 for i in U5)
SM.addConstrs(sum(x6[i, j] for j in V) == 1 for i in U6)
SM.addConstrs(sum(x7[i, j] for j in V) == 1 for i in U7)
SM.addConstrs(sum(x1[i, j] * amount1[i] for i in U1) <= capacity for j in V)
SM.addConstrs(sum(x2[i, j] * amount2[i] for i in U2) <= capacity for j in V)
SM.addConstrs(sum(x3[i, j] * amount3[i] for i in U3) <= capacity for j in V)
SM.addConstrs(sum(x4[i, j] * amount4[i] for i in U4) <= capacity for j in V)
SM.addConstrs(sum(x5[i, j] * amount5[i] for i in U5) <= capacity for j in V)
SM.addConstrs(sum(x6[i, j] * amount6[i] for i in U6) <= capacity for j in V)
SM.addConstrs(sum(x7[i, j] * amount7[i] for i in U7) <= capacity for j in V)
for i in U1:
for j in V:
SM.addGenConstrIndicator(x1[i, j], 1, y[j] == 1, name="ind1")
for i in U2:
for j in V:
SM.addGenConstrIndicator(x2[i, j], 1, y[j] == 1, name="ind2")
for i in U3:
for j in V:
SM.addGenConstrIndicator(x3[i, j], 1, y[j] == 1, name="ind3")
for i in U4:
for j in V:
SM.addGenConstrIndicator(x4[i, j], 1, y[j] == 1, name="ind4")
for i in U5:
for j in V:
SM.addGenConstrIndicator(x5[i, j], 1, y[j] == 1, name="ind5")
for i in U6:
for j in V:
SM.addGenConstrIndicator(x6[i, j], 1, y[j] == 1, name="ind6")
for i in U7:
for j in V:
SM.addGenConstrIndicator(x7[i, j], 1, y[j] == 1, name="ind7")
# Save model
SM.write('facilityPY.lp')
# Solve
SM.optimize()
print(SM.objVal)
x1 = SM.getAttr('x', x1)
x2 = SM.getAttr('x', x2)
x3 = SM.getAttr('x', x3)
x4 = SM.getAttr('x', x4)
x5 = SM.getAttr('x', x5)
x6 = SM.getAttr('x', x6)
x7 = SM.getAttr('x', x7)
y = SM.getAttr('x', y)
for j in V:
if y[j] == 1:
print('y', j, y[j])
for i in U1:
for j in V:
if 0.9 <= x1[i, j] <= 1.1:
plt.plot([lng1[i], lng[j]], [lat1[i], lat[j]], 'b')
plt.scatter(lng1[i], lat1[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U2:
for j in V:
if 0.9 <= x2[i, j] <= 1.1:
plt.plot([lng2[i], lng[j]], [lat2[i], lat[j]], 'b')
plt.scatter(lng2[i], lat2[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U3:
for j in V:
if 0.9 <= x3[i, j] <= 1.1:
plt.plot([lng3[i], lng[j]], [lat3[i], lat[j]], 'b')
plt.scatter(lng3[i], lat3[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U4:
for j in V:
if 0.9 <= x4[i, j] <= 1.1:
plt.plot([lng4[i], lng[j]], [lat4[i], lat[j]], 'b')
plt.scatter(lng4[i], lat4[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U5:
for j in V:
if 0.9 <= x5[i, j] <= 1.1:
plt.plot([lng5[i], lng[j]], [lat5[i], lat[j]], 'b')
plt.scatter(lng5[i], lat5[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U6:
for j in V:
if 0.9 <= x6[i, j] <= 1.1:
plt.plot([lng6[i], lng[j]], [lat6[i], lat[j]], 'b')
plt.scatter(lng6[i], lat6[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U7:
for j in V:
if 0.9 <= x7[i, j] <= 1.1:
plt.plot([lng7[i], lng[j]], [lat7[i], lat[j]], color='b')
plt.scatter(lng7[i], lat7[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.show()
# 画出订货商户地点
for i in U1:
plt.scatter(lng1[i], lat1[i], c='r')
plt.xlim(113, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U2:
plt.scatter(lng2[i], lat2[i], c='r')
plt.xlim(113, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U3:
plt.scatter(lng3[i], lat3[i], c='r')
plt.xlim(113, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U4:
plt.scatter(lng4[i], lat4[i], c='r')
plt.xlim(113, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U5:
plt.scatter(lng5[i], lat5[i], c='r')
plt.xlim(113, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U6:
plt.scatter(lng6[i], lat6[i], c='r')
plt.xlim(113, 114)
plt.ylim(22.4, 24)
plt.show()
for i in U7:
plt.scatter(lng7[i], lat7[i], c='r')
plt.xlim(113, 114)
plt.ylim(22.4, 24)
plt.show()
8.14
遗传算法建模开始
Numpy
numpy教程参看此博客：
https://www.qikegu.com/docs/3418#numpy%E6%95%B0%E7%BB%84%E8%BD%B4
- 获取数组元素的大小：
#获取数组元素的大小
import numpy as np
a = np.array([[1,2,3]])
print("数组元素大小：", a.itemsize, "字节")
# 输出
# 数组元素大小： 8 字节
- 获取数组的大小和行列数
import numpy as np
a = np.array([[1,2,3,4,5,6,7]])
print("大小:",a.size)
print("形状:",a.shape)
# 输出
# 大小: 7
# 形状: (1, 7)
- 重构数组对象
import numpy as np
a = np.array([[1,2],[3,4],[5,6]])
print("原数组：")
print(a)
a=a.reshape(2,3)
print("改变后：")
print(a)
# 输出如下：
原数组：
[[1 2]
[3 4]
[5 6]]
改变后：
[[1 2 3]
[4 5 6]]
- 数组切片
NumPy中，数组切片可以从数组中提取指定范围的数组元素。NumPy中的数组切片方法与python中的列表切片方法类似。
切片语法：
arr_name[start: end: step]
-
- [:]表示复制源列表
-
- 负的index表示，从后往前。-1表示最后一个元素。
>>> a = np.array([1,2,3,4,5,6,7,8,9])
>>> print(a[:5])
[1 2 3 4 5]
>>> print(a[1:5:2])
[2 4]
- linspace(start, end, number)
linspace()函数的作用是: 返回给定区间内均匀分布的值。下面的示例，在给定的区间5-15内返回10个均匀分布的值
示例
import numpy as np
a = np.linspace(0, 15, 5) # 打印5个值，这些值在给定的区间0-15上均匀地分布
print(a)
# 输出：
[ 0. 3.75 7.5 11.25 15. ]
-
获取数组元素中的最大值、最小值以及元素的和
ndarray.max()、ndarray.min()和ndarray.sum()函数，分别用于获取数组元素中的最大值、最小值以及元素和。
import numpy as np
a = np.array([1,2,3,10,15,4])
print("数组:",a)
print("最大值:",a.max())
print("最小值:",a.min())
print("元素总和:",a.sum())
# 输出：
数组: [ 1 2 3 10 15 4]
最大值: 15
最小值: 1
元素总和: 35
-
Numpy 数组轴
NumPy多维数组由轴表示，其中axis-0表示列，axis-1表示行。我们可以通过轴对列或行进行计算。例如，求某行元素的和。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-odyqtbSP-1677154847156)(F:\2020Summer\暑期实习科研\实习\图片\np数轴.JPG)]
计算每一列中的最大元素、每一行中的最小元素，以及每一行的和：
import numpy as np
a = np.array([[1,2,30],[10,15,4]])
print("数组:",a)
print("每列的最大元素:",a.max(axis = 0))
print("每行的最小元素:",a.min(axis = 1))
print("每行的和:",a.sum(axis = 1))
# 输出：
数组: [[ 1 2 30]
[10 15 4]]
每列的最大元素: [10 15 30]
每行的最小元素: [1 4]
每行的和: [33 29]
-
数组间的算术运算
import numpy as np a = np.array([[1,2,30],[10,15,4]]) b = np.array([[1,2,3],[12, 19, 29]]) print("a+b\n",a+b) print("axb\n",a*b) # 同样位置的数字相乘之后，在同样位置构成的数组 print("a/b\n",a/b) # 同样位置的数字相除之后，在同样位置构成的数组 # 输出 a+b [[ 2 4 33] [22 34 33]] axb [[ 1 4 90] [120 285 116]] a/b [[ 1. 1. 10. ] [ 0.83333333 0.78947368 0.13793103]] -
数组拼接
numpy中，可以垂直或水平拼接2个数组。
import numpy as np a = np.array([[1,2,30],[10,15,4]]) b = np.array([[1,2,3],[12, 19, 29]]) print("垂直拼接\n",np.vstack((a,b))); print("水平拼接\n",np.hstack((a,b))) # 输出 垂直拼接 [[ 1 2 30] [10 15 4] [ 1 2 3] [12 19 29]] 水平拼接 [[ 1 2 30 1 2 3] [10 15 4 12 19 29]]
11.7
过滤商业客户
import json
filename8 = (r'Bsn_id.json')
Bsn_ID = json.load(open(filename8))
k = len(Bsn_ID['RECORDS'])
filename9 = (r'ZengCheng_XT_ID.json')
ZCXT = json.load(open(filename9))
filename1 = (r'MonList.json')
Mon_str = json.load(open(filename1))
filename2 = (r'TuesList.json')
Tues_str = json.load(open(filename2))
filename3 = (r'WedList.json')
Wed_str = json.load(open(filename3))
filename4 = (r'ThurList.json')
Thur_str = json.load(open(filename4))
filename5 = (r'FriList.json')
Fri_str = json.load(open(filename5))
filename6 = (r'SatList.json')
Sat_str = json.load(open(filename6))
filename7 = (r'SunList.json')
Sun_str = json.load(open(filename7))
Mon_Phar = {}
Tues_Phar = {}
Wed_Phar = {}
Thur_Phar = {}
Fri_Phar = {}
Sat_Phar = {}
Sun_Phar = {}
# 过滤商业客户
for i in range(0, k):
j = Bsn_ID['RECORDS'][i]['delivery_addr_id']
if j in Mon_str.keys():
Mon_Phar[j] = Mon_str[j]
if j in Tues_str.keys():
Tues_Phar[j] = Tues_str[j]
if j in Wed_str.keys():
Wed_Phar[j] = Wed_str[j]
if j in Thur_str.keys():
Thur_Phar[j] = Thur_str[j]
if j in Fri_str.keys():
Fri_Phar[j] = Fri_str[j]
if j in Sat_str.keys():
Sat_Phar[j] = Sat_str[j]
if j in Sun_str.keys():
Sun_Phar[j] = Sun_str[j]
json.dump(Mon_Phar, open(r'E:\PycharmProjects\Project1\AnotherTest\Mon_Phar.json', 'w'))
json.dump(Tues_Phar, open(r'E:\PycharmProjects\Project1\AnotherTest\Tues_Phar.json', 'w'))
json.dump(Wed_Phar, open(r'E:\PycharmProjects\Project1\AnotherTest\Wed_Phar.json', 'w'))
json.dump(Thur_Phar, open(r'E:\PycharmProjects\Project1\AnotherTest\Thur_Phar.json', 'w'))
json.dump(Fri_Phar, open(r'E:\PycharmProjects\Project1\AnotherTest\Fri_Phar.json', 'w'))
json.dump(Sat_Phar, open(r'E:\PycharmProjects\Project1\AnotherTest\Sat_Phar.json', 'w'))
json.dump(Sun_Phar, open(r'E:\PycharmProjects\Project1\AnotherTest\Sun_Phar.json', 'w'))
11.9
带画图例的测试（Another Test）
import matplotlib.pyplot as plt
import gurobipy as gp
from gurobipy import *
import json
import math
import numpy as np
# ===================================================================================
# 通过经纬度测算距离
# INPUTS: A点经度 A点纬度 B点经度 B点纬度
# OUTPUTS: 返回值为考虑球面的两点直线距离,单位：km
# ===================================================================================
def distance(lng_A, lat_A, lng_B, lat_B):
# 地球半径
R = 6371.004
pi = 3.141592654
Mlng_A = lng_A
Mlat_A = 90 - lat_A
Mlng_B = lng_B
Mlat_B = 90 - lat_B
C = math.sin(Mlat_A * pi / 180) * math.sin(Mlat_B * pi / 180) * math.cos((Mlng_A - Mlng_B) * pi / 180) + math.cos(
Mlat_A * pi / 180) * math.cos(Mlat_B * pi / 180)
if C >= 1:
C = 1
Distance = R * math.acos(C)
return Distance
# =================================================================================
# 将以str作为key值的字典转化为以int为key值的字典,并只保留前n个key-value对
# 输入：dict_str, n
# 输出：dict_int
# ===================================================================================
def Key_StrToInt(dict_str, n):
dict_int = {}
Keys = dict_str.keys()
Keys = list(Keys)
Keys = Keys[0:n]
for i in Keys:
i = int(i)
temp = {
i: dict_str[str(i)]
}
dict_int.update(temp)
return dict_int
# ==================读入以str为key值的字典=====================
filename1 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Mon.json')
Mon_str = json.load(open(filename1))
filename2 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Tues.json')
Tues_str = json.load(open(filename2))
filename3 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Wed.json')
Wed_str = json.load(open(filename3))
filename4 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Thur.json')
Thur_str = json.load(open(filename4))
filename5 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Fri.json')
Fri_str = json.load(open(filename5))
filename6 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Sat.json')
Sat_str = json.load(open(filename6))
filename7 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Sun.json')
Sun_str = json.load(open(filename7))
# ===================================================================================
n = 25 # 定义取每个场景的多少个数据
# 把字典转化成以int数字为key值的字典
Mon_int = Key_StrToInt(Mon_str, n)
Tues_int = Key_StrToInt(Tues_str, n)
Wed_int = Key_StrToInt(Wed_str, n)
Thur_int = Key_StrToInt(Thur_str, n)
Fri_int = Key_StrToInt(Fri_str, n)
Sat_int = Key_StrToInt(Sat_str, n)
Sun_int = Key_StrToInt(Sun_str, n)
U1, amount1, lng1, lat1 = gp.multidict(Mon_int)
U2, amount2, lng2, lat2 = gp.multidict(Tues_int)
U3, amount3, lng3, lat3 = gp.multidict(Wed_int)
U4, amount4, lng4, lat4 = gp.multidict(Thur_int)
U5, amount5, lng5, lat5 = gp.multidict(Fri_int)
U6, amount6, lng6, lat6 = gp.multidict(Sat_int)
U7, amount7, lng7, lat7 = gp.multidict(Sun_int)
capacity = 300000
FixedCost = 5000
T = 200
# 求两个集合的并集
U1 = list(U1)
U2 = list(U2)
U3 = list(U3)
U4 = list(U4)
U5 = list(U5)
U6 = list(U6)
U7 = list(U7)
V = list(set(U1).union(U2, U3, U4, U5, U6, U7))
# T = [T for h in range(len(V))]
Z = [capacity for i in V]
C = [FixedCost for c in V]
# 把并集之后的数据写入lng[],lat[]，作为仓库潜在选址集合坐标
lng = {}
lat = {}
for k in V:
if k in U1:
a = lng1[k]
b = lat1[k]
elif k in U2:
a = lng2[k]
b = lat2[k]
elif k in U3:
a = lng3[k]
b = lat3[k]
elif k in U4:
a = lng4[k]
b = lat4[k]
elif k in U5:
a = lng5[k]
b = lat5[k]
elif k in U6:
a = lng6[k]
b = lat6[k]
elif k in U7:
a = lng7[k]
b = lat7[k]
x = {
k: a
}
y = {
k: b
}
lng.update(x)
lat.update(y)
# 计算各场景下距离矩阵D_jk
D1 = {
(i, j): distance(lng1[i], lat1[i], lng[j], lat[j]) for i in U1 for j in V
}
D2 = {
(i, j): distance(lng2[i], lat2[i], lng[j], lat[j]) for i in U2 for j in V
}
D3 = {
(i, j): distance(lng3[i], lat3[i], lng[j], lat[j]) for i in U3 for j in V
}
D4 = {
(i, j): distance(lng4[i], lat4[i], lng[j], lat[j]) for i in U4 for j in V
}
D5 = {
(i, j): distance(lng5[i], lat5[i], lng[j], lat[j]) for i in U5 for j in V
}
D6 = {
(i, j): distance(lng6[i], lat6[i], lng[j], lat[j]) for i in U6 for j in V
}
D7 = {
(i, j): distance(lng7[i], lat7[i], lng[j], lat[j]) for i in U7 for j in V
}
# 创建解集下标
x1 = [(i, j) for i in U1 for j in V]
x2 = [(i, j) for i in U2 for j in V]
x3 = [(i, j) for i in U3 for j in V]
x4 = [(i, j) for i in U4 for j in V]
x5 = [(i, j) for i in U5 for j in V]
x6 = [(i, j) for i in U6 for j in V]
x7 = [(i, j) for i in U7 for j in V]
x1 = tuplelist(x1)
x2 = tuplelist(x2)
x3 = tuplelist(x3)
x4 = tuplelist(x4)
x5 = tuplelist(x5)
x6 = tuplelist(x6)
x7 = tuplelist(x7)
V = tuplelist(V)
# Model
SM = gp.Model("Scenario")
# 添加变量
y = SM.addVars(V,
vtype=GRB.BINARY,
obj=C,
name="y")
x1 = SM.addVars(x1, vtype=GRB.BINARY, name="x1")
x2 = SM.addVars(x2, vtype=GRB.BINARY, name="x2")
x3 = SM.addVars(x3, vtype=GRB.BINARY, name="x3")
x4 = SM.addVars(x4, vtype=GRB.BINARY, name="x4")
x5 = SM.addVars(x5, vtype=GRB.BINARY, name="x5")
x6 = SM.addVars(x6, vtype=GRB.BINARY, name="x6")
x7 = SM.addVars(x7, vtype=GRB.BINARY, name="x7")
SM.update()
# 设置目标函数
SM.setObjective(
sum(y[j] * FixedCost for j in V)
+ sum(x1[i, j] * D1[i, j] for i in U1 for j in V)
+ sum(x2[i, j] * D2[i, j] for i in U2 for j in V)
+ sum(x3[i, j] * D3[i, j] for i in U3 for j in V)
+ sum(x4[i, j] * D4[i, j] for i in U4 for j in V)
+ sum(x5[i, j] * D5[i, j] for i in U5 for j in V)
+ sum(x6[i, j] * D6[i, j] for i in U6 for j in V)
+ sum(x7[i, j] * D7[i, j] for i in U7 for j in V),
GRB.MINIMIZE)
# 添加约束
SM.addConstrs(T >= D1[i, j] * x1[i, j] for i in U1 for j in V)
SM.addConstrs(T >= D2[i, j] * x2[i, j] for i in U2 for j in V)
SM.addConstrs(T >= D3[i, j] * x3[i, j] for i in U3 for j in V)
SM.addConstrs(T >= D4[i, j] * x4[i, j] for i in U4 for j in V)
SM.addConstrs(T >= D5[i, j] * x5[i, j] for i in U5 for j in V)
SM.addConstrs(T >= D6[i, j] * x6[i, j] for i in U6 for j in V)
SM.addConstrs(T >= D7[i, j] * x7[i, j] for i in U7 for j in V)
SM.addConstrs(sum(x1[i, j] for j in V) == 1 for i in U1)
SM.addConstrs(sum(x2[i, j] for j in V) == 1 for i in U2)
SM.addConstrs(sum(x3[i, j] for j in V) == 1 for i in U3)
SM.addConstrs(sum(x4[i, j] for j in V) == 1 for i in U4)
SM.addConstrs(sum(x5[i, j] for j in V) == 1 for i in U5)
SM.addConstrs(sum(x6[i, j] for j in V) == 1 for i in U6)
SM.addConstrs(sum(x7[i, j] for j in V) == 1 for i in U7)
SM.addConstrs(sum(x1[i, j] * amount1[i] for i in U1) <= capacity for j in V)
SM.addConstrs(sum(x2[i, j] * amount2[i] for i in U2) <= capacity for j in V)
SM.addConstrs(sum(x3[i, j] * amount3[i] for i in U3) <= capacity for j in V)
SM.addConstrs(sum(x4[i, j] * amount4[i] for i in U4) <= capacity for j in V)
SM.addConstrs(sum(x5[i, j] * amount5[i] for i in U5) <= capacity for j in V)
SM.addConstrs(sum(x6[i, j] * amount6[i] for i in U6) <= capacity for j in V)
SM.addConstrs(sum(x7[i, j] * amount7[i] for i in U7) <= capacity for j in V)
# SM.addConstr(sum(y[j] for j in V) == 1)
for i in U1:
for j in V:
SM.addGenConstrIndicator(x1[i, j], 1, y[j] == 1, name="ind1")
for i in U2:
for j in V:
SM.addGenConstrIndicator(x2[i, j], 1, y[j] == 1, name="ind2")
for i in U3:
for j in V:
SM.addGenConstrIndicator(x3[i, j], 1, y[j] == 1, name="ind3")
for i in U4:
for j in V:
SM.addGenConstrIndicator(x4[i, j], 1, y[j] == 1, name="ind4")
for i in U5:
for j in V:
SM.addGenConstrIndicator(x5[i, j], 1, y[j] == 1, name="ind5")
for i in U6:
for j in V:
SM.addGenConstrIndicator(x6[i, j], 1, y[j] == 1, name="ind6")
for i in U7:
for j in V:
SM.addGenConstrIndicator(x7[i, j], 1, y[j] == 1, name="ind7")
# Save model
SM.write('facilityPY.lp')
# Solve
SM.optimize()
print(SM.objVal)
x1 = SM.getAttr('x', x1)
x2 = SM.getAttr('x', x2)
x3 = SM.getAttr('x', x3)
x4 = SM.getAttr('x', x4)
x5 = SM.getAttr('x', x5)
x6 = SM.getAttr('x', x6)
x7 = SM.getAttr('x', x7)
y = SM.getAttr('x', y)
for j in V:
if y[j] == 1:
print('y', j, y[j])
flag = 1
for i in U1:
for j in V:
if 0.9 <= x1[i, j] <= 1.1:
if 0.9 < flag < 1.1:
plt.plot([lng1[i], lng[j]], [lat1[i], lat[j]], 'b')
s1 = plt.scatter(lng1[i], lat1[i], color='r')
s2 = plt.scatter(lng[j], lat[j], color='g')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
else:
plt.plot([lng1[i], lng[j]], [lat1[i], lat[j]], 'b')
plt.scatter(lng1[i], lat1[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.title('(a)')
plt.show()
for i in U2:
for j in V:
if 0.9 <= x2[i, j] <= 1.1:
if 1.9 < flag < 2.1:
plt.plot([lng2[i], lng[j]], [lat2[i], lat[j]])
plt.scatter(lng2[i], lat2[i], color='r', label='merchant')
plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend(loc="upper left")
flag = flag + 1
else:
plt.plot([lng2[i], lng[j]], [lat2[i], lat[j]], 'b')
plt.scatter(lng2[i], lat2[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.title('(b)')
plt.show()
for i in U3:
for j in V:
if 0.9 <= x3[i, j] <= 1.1:
if 2.9 < flag < 3.1:
plt.plot([lng3[i], lng[j]], [lat3[i], lat[j]], 'b')
plt.scatter(lng3[i], lat3[i], color='r', label='merchant')
plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend(loc="upper left")
flag = flag + 1
else:
plt.plot([lng3[i], lng[j]], [lat3[i], lat[j]], 'b')
plt.scatter(lng3[i], lat3[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.title('(c)')
plt.show()
for i in U4:
for j in V:
if 0.9 <= x4[i, j] <= 1.1:
if 3.9 < flag < 4.1:
s1 = plt.scatter(lng4[i], lat4[i], color='r', label='merchant')
plt.plot([lng4[i], lng[j]], [lat4[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
else:
plt.plot([lng4[i], lng[j]], [lat4[i], lat[j]], 'b')
plt.scatter(lng4[i], lat4[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('(d)')
plt.show()
for i in U5:
for j in V:
if 0.9 <= x5[i, j] <= 1.1:
if 4.9 < flag < 5.1:
s1 = plt.scatter(lng5[i], lat5[i], color='r', label='merchant')
plt.plot([lng5[i], lng[j]], [lat5[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
else:
plt.plot([lng5[i], lng[j]], [lat5[i], lat[j]], 'b')
plt.scatter(lng5[i], lat5[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('(e)')
plt.show()
for i in U6:
for j in V:
if 0.9 <= x6[i, j] <= 1.1:
if 5.9 < flag < 6.1:
s1 = plt.scatter(lng6[i], lat6[i], color='r', label='merchant')
plt.plot([lng6[i], lng[j]], [lat6[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
else:
plt.plot([lng6[i], lng[j]], [lat6[i], lat[j]], 'b')
plt.scatter(lng6[i], lat6[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('(f)')
plt.show()
for i in U7:
for j in V:
if 0.9 <= x7[i, j] <= 1.1:
if 6.9 < flag < 7.1:
s1 = plt.scatter(lng7[i], lat7[i], color='r', label='merchant')
plt.plot([lng7[i], lng[j]], [lat7[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
else:
plt.plot([lng7[i], lng[j]], [lat7[i], lat[j]], color='b')
plt.scatter(lng7[i], lat7[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('(g)')
plt.show()
画到一张图上的测试
folium示例
import folium
import webbrowser
m=folium.Map(location=[40.009867,116.485994],zoom_start=10) # 绘制地图，确定聚焦点
folium.Marker([40.2,116.7],popup='<b>浮标上面的那个文字</b>').add_to(m) # 定一个点，放到地图m上
folium.Marker([40.22,116.72],popup='<b>浮标上面的那个文字</b>',icon=folium.Icon(color='red')).add_to(m)
# 把浮标变成红色
folium.Marker([40.24,116.74],popup='<b>浮标上面的那个文字</b>',icon=folium.Icon(color='green',icon='info-sign')).add_to(m)
# 浮标改图样
#标记一个空心的圈
folium.Circle(
location=[40.2,117.7],
radius=10000,
color='crimson',
popup='popup',
fill=False
).add_to(m)
#标记一个实心圆
folium.CircleMarker(
location=[39.2,117.7],
radius=100,
popup='popup',
color='#DC143C',#圈的颜色
fill=True,
fill_color='#6495ED' #填充颜色
).add_to(m)
m.save('f1.html'
webbrowser.open('f1.html')
11.9
标注到html网页上
# coding=utf-8
import matplotlib.pyplot as plt
import gurobipy as gp
from gurobipy import *
import json
import math
import folium
import webbrowser
# ===================================================================================
# 通过经纬度测算距离
# INPUTS: A点经度 A点纬度 B点经度 B点纬度
# OUTPUTS: 返回值为考虑球面的两点直线距离,单位：km
# ===================================================================================
def distance(lng_A, lat_A, lng_B, lat_B):
# 地球半径
R = 6371.004
pi = 3.141592654
Mlng_A = lng_A
Mlat_A = 90 - lat_A
Mlng_B = lng_B
Mlat_B = 90 - lat_B
C = math.sin(Mlat_A * pi / 180) * math.sin(Mlat_B * pi / 180) * math.cos((Mlng_A - Mlng_B) * pi / 180) + math.cos(
Mlat_A * pi / 180) * math.cos(Mlat_B * pi / 180)
if C >= 1:
C = 1
Distance = R * math.acos(C)
return Distance
# =================================================================================
# 将以str作为key值的字典转化为以int为key值的字典,并只保留前n个key-value对
# 输入：dict_str, n
# 输出：dict_int
# ===================================================================================
def Key_StrToInt(dict_str, n):
dict_int = {}
Keys = dict_str.keys()
Keys = list(Keys)
Keys = Keys[0:n]
for i in Keys:
i = int(i)
temp = {
i: dict_str[str(i)]
}
dict_int.update(temp)
return dict_int
# ==================读入以str为key值的字典=====================
filename1 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Mon.json')
Mon_str = json.load(open(filename1))
filename2 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Tues.json')
Tues_str = json.load(open(filename2))
filename3 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Wed.json')
Wed_str = json.load(open(filename3))
filename4 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Thur.json')
Thur_str = json.load(open(filename4))
filename5 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Fri.json')
Fri_str = json.load(open(filename5))
filename6 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Sat.json')
Sat_str = json.load(open(filename6))
filename7 = (r'E:\PycharmProjects\Project1\AnotherTest\HP_Sun.json')
Sun_str = json.load(open(filename7))
# ===================================================================================
n = 120 # 定义取每个场景的多少个数据
# 把字典转化成以int数字为key值的字典
Mon_int = Key_StrToInt(Mon_str, n)
Tues_int = Key_StrToInt(Tues_str, n)
Wed_int = Key_StrToInt(Wed_str, n)
Thur_int = Key_StrToInt(Thur_str, n)
Fri_int = Key_StrToInt(Fri_str, n)
Sat_int = Key_StrToInt(Sat_str, n)
Sun_int = Key_StrToInt(Sun_str, n)
U1, amount1, lng1, lat1 = gp.multidict(Mon_int)
U2, amount2, lng2, lat2 = gp.multidict(Tues_int)
U3, amount3, lng3, lat3 = gp.multidict(Wed_int)
U4, amount4, lng4, lat4 = gp.multidict(Thur_int)
U5, amount5, lng5, lat5 = gp.multidict(Fri_int)
U6, amount6, lng6, lat6 = gp.multidict(Sat_int)
U7, amount7, lng7, lat7 = gp.multidict(Sun_int)
capacity = 300000
FixedCost = 5000
T = 200
# 求两个集合的并集
U1 = list(U1)
U2 = list(U2)
U3 = list(U3)
U4 = list(U4)
U5 = list(U5)
U6 = list(U6)
U7 = list(U7)
V = list(set(U1).union(U2, U3, U4, U5, U6, U7))
# T = [T for h in range(len(V))]
Z = [capacity for i in V]
C = [FixedCost for c in V]
# 把并集之后的数据写入lng[],lat[]，作为仓库潜在选址集合坐标
lng = {}
lat = {}
for k in V:
if k in U1:
a = lng1[k]
b = lat1[k]
elif k in U2:
a = lng2[k]
b = lat2[k]
elif k in U3:
a = lng3[k]
b = lat3[k]
elif k in U4:
a = lng4[k]
b = lat4[k]
elif k in U5:
a = lng5[k]
b = lat5[k]
elif k in U6:
a = lng6[k]
b = lat6[k]
elif k in U7:
a = lng7[k]
b = lat7[k]
x = {
k: a
}
y = {
k: b
}
lng.update(x)
lat.update(y)
# 计算各场景下距离矩阵D_jk
D1 = {
(i, j): distance(lng1[i], lat1[i], lng[j], lat[j]) for i in U1 for j in V
}
D2 = {
(i, j): distance(lng2[i], lat2[i], lng[j], lat[j]) for i in U2 for j in V
}
D3 = {
(i, j): distance(lng3[i], lat3[i], lng[j], lat[j]) for i in U3 for j in V
}
D4 = {
(i, j): distance(lng4[i], lat4[i], lng[j], lat[j]) for i in U4 for j in V
}
D5 = {
(i, j): distance(lng5[i], lat5[i], lng[j], lat[j]) for i in U5 for j in V
}
D6 = {
(i, j): distance(lng6[i], lat6[i], lng[j], lat[j]) for i in U6 for j in V
}
D7 = {
(i, j): distance(lng7[i], lat7[i], lng[j], lat[j]) for i in U7 for j in V
}
# 创建解集下标
x1 = [(i, j) for i in U1 for j in V]
x2 = [(i, j) for i in U2 for j in V]
x3 = [(i, j) for i in U3 for j in V]
x4 = [(i, j) for i in U4 for j in V]
x5 = [(i, j) for i in U5 for j in V]
x6 = [(i, j) for i in U6 for j in V]
x7 = [(i, j) for i in U7 for j in V]
x1 = tuplelist(x1)
x2 = tuplelist(x2)
x3 = tuplelist(x3)
x4 = tuplelist(x4)
x5 = tuplelist(x5)
x6 = tuplelist(x6)
x7 = tuplelist(x7)
V = tuplelist(V)
# Model
SM = gp.Model("Scenario")
# 添加变量
y = SM.addVars(V,
vtype=GRB.BINARY,
obj=C,
name="y")
x1 = SM.addVars(x1, vtype=GRB.BINARY, name="x1")
x2 = SM.addVars(x2, vtype=GRB.BINARY, name="x2")
x3 = SM.addVars(x3, vtype=GRB.BINARY, name="x3")
x4 = SM.addVars(x4, vtype=GRB.BINARY, name="x4")
x5 = SM.addVars(x5, vtype=GRB.BINARY, name="x5")
x6 = SM.addVars(x6, vtype=GRB.BINARY, name="x6")
x7 = SM.addVars(x7, vtype=GRB.BINARY, name="x7")
SM.update()
# 设置目标函数
SM.setObjective(
sum(y[j] * FixedCost for j in V)
+ sum(x1[i, j] * D1[i, j] for i in U1 for j in V)
+ sum(x2[i, j] * D2[i, j] for i in U2 for j in V)
+ sum(x3[i, j] * D3[i, j] for i in U3 for j in V)
+ sum(x4[i, j] * D4[i, j] for i in U4 for j in V)
+ sum(x5[i, j] * D5[i, j] for i in U5 for j in V)
+ sum(x6[i, j] * D6[i, j] for i in U6 for j in V)
+ sum(x7[i, j] * D7[i, j] for i in U7 for j in V),
GRB.MINIMIZE)
# 添加约束
SM.addConstrs(T >= D1[i, j] * x1[i, j] for i in U1 for j in V)
SM.addConstrs(T >= D2[i, j] * x2[i, j] for i in U2 for j in V)
SM.addConstrs(T >= D3[i, j] * x3[i, j] for i in U3 for j in V)
SM.addConstrs(T >= D4[i, j] * x4[i, j] for i in U4 for j in V)
SM.addConstrs(T >= D5[i, j] * x5[i, j] for i in U5 for j in V)
SM.addConstrs(T >= D6[i, j] * x6[i, j] for i in U6 for j in V)
SM.addConstrs(T >= D7[i, j] * x7[i, j] for i in U7 for j in V)
SM.addConstrs(sum(x1[i, j] for j in V) == 1 for i in U1)
SM.addConstrs(sum(x2[i, j] for j in V) == 1 for i in U2)
SM.addConstrs(sum(x3[i, j] for j in V) == 1 for i in U3)
SM.addConstrs(sum(x4[i, j] for j in V) == 1 for i in U4)
SM.addConstrs(sum(x5[i, j] for j in V) == 1 for i in U5)
SM.addConstrs(sum(x6[i, j] for j in V) == 1 for i in U6)
SM.addConstrs(sum(x7[i, j] for j in V) == 1 for i in U7)
SM.addConstrs(sum(x1[i, j] * amount1[i] for i in U1) <= capacity for j in V)
SM.addConstrs(sum(x2[i, j] * amount2[i] for i in U2) <= capacity for j in V)
SM.addConstrs(sum(x3[i, j] * amount3[i] for i in U3) <= capacity for j in V)
SM.addConstrs(sum(x4[i, j] * amount4[i] for i in U4) <= capacity for j in V)
SM.addConstrs(sum(x5[i, j] * amount5[i] for i in U5) <= capacity for j in V)
SM.addConstrs(sum(x6[i, j] * amount6[i] for i in U6) <= capacity for j in V)
SM.addConstrs(sum(x7[i, j] * amount7[i] for i in U7) <= capacity for j in V)
SM.addConstr(sum(y[j] for j in V) == 1)
for i in U1:
for j in V:
SM.addGenConstrIndicator(x1[i, j], 1, y[j] == 1, name="ind1")
for i in U2:
for j in V:
SM.addGenConstrIndicator(x2[i, j], 1, y[j] == 1, name="ind2")
for i in U3:
for j in V:
SM.addGenConstrIndicator(x3[i, j], 1, y[j] == 1, name="ind3")
for i in U4:
for j in V:
SM.addGenConstrIndicator(x4[i, j], 1, y[j] == 1, name="ind4")
for i in U5:
for j in V:
SM.addGenConstrIndicator(x5[i, j], 1, y[j] == 1, name="ind5")
for i in U6:
for j in V:
SM.addGenConstrIndicator(x6[i, j], 1, y[j] == 1, name="ind6")
for i in U7:
for j in V:
SM.addGenConstrIndicator(x7[i, j], 1, y[j] == 1, name="ind7")
# Save model
SM.write('facilityPY.lp')
# Solve
SM.optimize()
print(SM.objVal)
x1 = SM.getAttr('x', x1)
x2 = SM.getAttr('x', x2)
x3 = SM.getAttr('x', x3)
x4 = SM.getAttr('x', x4)
x5 = SM.getAttr('x', x5)
x6 = SM.getAttr('x', x6)
x7 = SM.getAttr('x', x7)
y = SM.getAttr('x', y)
Draw_map = folium.Map(location=[23.20, 113.45], zoom_start=12)
for j in V:
if y[j] == 1:
print('y', j, y[j])
folium.Marker([lat[j], lng[j]], popup='<b>仓库</b>').add_to(Draw_map) # 定一个点，放到地图m上
flag = 1
for i in U1:
for j in V:
if 0.9 <= x1[i, j] <= 1.1:
if 0.9 < flag < 1.1:
# plt.plot([lng1[i], lng[j]], [lat1[i], lat[j]], 'b')
s1 = plt.scatter(lng1[i], lat1[i], color='r')
s2 = plt.scatter(lng[j], lat[j], color='g')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
folium.Circle(
radius=200,
location=[lat1[i], lng1[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
else:
# plt.plot([lng1[i], lng[j]], [lat1[i], lat[j]], 'b')
plt.scatter(lng1[i], lat1[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
folium.Marker([lat[j], lng[j]], popup='<b>ddd</b>').add_to(Draw_map) # 定一个点，放到地图m上
folium.Circle(
radius=200,
location=[lat1[i], lng1[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
for i in U2:
for j in V:
if 0.9 <= x2[i, j] <= 1.1:
if 1.9 < flag < 2.1:
# plt.plot([lng2[i], lng[j]], [lat2[i], lat[j]])
plt.scatter(lng2[i], lat2[i], color='r', label='merchant')
plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend(loc="upper left")
flag = flag + 1
folium.Circle(
radius=200,
location=[lat2[i], lng2[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
else:
# plt.plot([lng2[i], lng[j]], [lat2[i], lat[j]], 'b')
plt.scatter(lng2[i], lat2[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
folium.Circle(
radius=200,
location=[lat2[i], lng2[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
for i in U3:
for j in V:
if 0.9 <= x3[i, j] <= 1.1:
if 2.9 < flag < 3.1:
# plt.plot([lng3[i], lng[j]], [lat3[i], lat[j]], 'b')
plt.scatter(lng3[i], lat3[i], color='r', label='merchant')
plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend(loc="upper left")
flag = flag + 1
folium.Circle(
radius=200,
location=[lat3[i], lng3[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
else:
# plt.plot([lng3[i], lng[j]], [lat3[i], lat[j]], 'b')
plt.scatter(lng3[i], lat3[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
folium.Circle(
radius=200,
location=[lat3[i], lng3[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
for i in U4:
for j in V:
if 0.9 <= x4[i, j] <= 1.1:
if 3.9 < flag < 4.1:
s1 = plt.scatter(lng4[i], lat4[i], color='r', label='merchant')
# plt.plot([lng4[i], lng[j]], [lat4[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
folium.Circle(
radius=200,
location=[lat4[i], lng4[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
else:
# plt.plot([lng4[i], lng[j]], [lat4[i], lat[j]], 'b')
plt.scatter(lng4[i], lat4[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
folium.Circle(
radius=200,
location=[lat4[i], lng4[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
for i in U5:
for j in V:
if 0.9 <= x5[i, j] <= 1.1:
if 4.9 < flag < 5.1:
s1 = plt.scatter(lng5[i], lat5[i], color='r', label='merchant')
# plt.plot([lng5[i], lng[j]], [lat5[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
folium.Circle(
radius=200,
location=[lat5[i], lng5[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
else:
# plt.plot([lng5[i], lng[j]], [lat5[i], lat[j]], 'b')
plt.scatter(lng5[i], lat5[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
folium.Circle(
radius=200,
location=[lat5[i], lng5[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
for i in U6:
for j in V:
if 0.9 <= x6[i, j] <= 1.1:
if 5.9 < flag < 6.1:
s1 = plt.scatter(lng6[i], lat6[i], color='r', label='merchant')
# plt.plot([lng6[i], lng[j]], [lat6[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
folium.Circle(
radius=200,
location=[lat6[i], lng6[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
else:
# plt.plot([lng6[i], lng[j]], [lat6[i], lat[j]], 'b')
plt.scatter(lng6[i], lat6[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
folium.Circle(
radius=200,
location=[lat6[i], lng6[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
for i in U7:
for j in V:
if 0.9 <= x7[i, j] <= 1.1:
if 6.9 < flag < 7.1:
s1 = plt.scatter(lng7[i], lat7[i], color='r', label='merchant')
# plt.plot([lng7[i], lng[j]], [lat7[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
folium.Circle(
radius=200,
location=[lat7[i], lng7[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
else:
# plt.plot([lng7[i], lng[j]], [lat7[i], lat[j]], color='b')
plt.scatter(lng7[i], lat7[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
folium.Circle(
radius=200,
location=[lat7[i], lng7[i]],
popup='popup',
color='crimson',
fill=False,
).add_to(Draw_map)
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('HuangPu')
plt.show()
Draw_map.save('HuangPu.html')
webbrowser.open('HuangPu.html')
钟落潭数据筛选成json
from openpyxl import load_workbook
import os
import time
import json
filename1 = (r'E:\PycharmProjects\Project1\Mon.json')
Mon = json.load(open(filename1))
filename3 = (r'E:\PycharmProjects\Project1\1_送货地址id.json')
ZLT = json.load(open(filename3))
filename2 = (r'E:\PycharmProjects\Project1\Tues.json')
Tues = json.load(open(filename2))
filename4 = (r'E:\PycharmProjects\Project1\Wed.json')
Wed = json.load(open(filename4))
filename5 = (r'E:\PycharmProjects\Project1\Thur.json')
Thur = json.load(open(filename5))
filename6 = (r'E:\PycharmProjects\Project1\Fri.json')
Fri = json.load(open(filename6))
filename7 = (r'E:\PycharmProjects\Project1\Sat.json')
Sat = json.load(open(filename7))
filename8 = (r'E:\PycharmProjects\Project1\Sun.json')
Sun = json.load(open(filename8))
filename9 = (r'E:\PycharmProjects\Project1\AnotherTest\ZhongLT.json')
NONO = json.load(open(filename9))
ZLT_CK = {}
for i in range(1, len(ZLT['RECORDS'])):
temp = {}
if ZLT['RECORDS'][i]['delivery_addr_id'] in NONO.keys():
temp[ZLT['RECORDS'][i]['delivery_addr_id']] = NONO[ZLT['RECORDS'][i]['delivery_addr_id']]
ZLT_CK.update(temp)
print(len(ZLT_CK))
json.dump(ZLT_CK,open(r'E:\PycharmProjects\Project1\AnotherTest\ZLT_CK.json','w'))
for i in range(1, len(ZLT['RECORDS'])):
temp = {}
if ZLT['RECORDS'][i]['delivery_addr_id'] in Mon.keys():
temp[ZLT['RECORDS'][i]['delivery_addr_id']] = Mon[ZLT['RECORDS'][i]['delivery_addr_id']]
ZLT_CK.update(temp)
for i in range(1, len(ZLT['RECORDS'])):
temp = {}
if ZLT['RECORDS'][i]['delivery_addr_id'] not in ZLT_CK.keys():
if ZLT['RECORDS'][i]['delivery_addr_id'] in Tues:
temp[ZLT['RECORDS'][i]['delivery_addr_id']] = Tues[ZLT['RECORDS'][i]['delivery_addr_id']]
ZLT_CK.update(temp)
print(len(ZLT_CK))
for i in range(1, len(ZLT['RECORDS'])):
temp = {}
if ZLT['RECORDS'][i]['delivery_addr_id'] not in ZLT_CK.keys():
if ZLT['RECORDS'][i]['delivery_addr_id'] in Wed:
temp[ZLT['RECORDS'][i]['delivery_addr_id']] = Wed[ZLT['RECORDS'][i]['delivery_addr_id']]
ZLT_CK.update(temp)
for i in range(1, len(ZLT['RECORDS'])):
temp = {}
if ZLT['RECORDS'][i]['delivery_addr_id'] not in ZLT_CK.keys():
if ZLT['RECORDS'][i]['delivery_addr_id'] in Thur:
temp[ZLT['RECORDS'][i]['delivery_addr_id']] = Thur[ZLT['RECORDS'][i]['delivery_addr_id']]
ZLT_CK.update(temp)
for i in range(1, len(ZLT['RECORDS'])):
temp = {}
if ZLT['RECORDS'][i]['delivery_addr_id'] not in ZLT_CK.keys():
if ZLT['RECORDS'][i]['delivery_addr_id'] in Fri:
temp[ZLT['RECORDS'][i]['delivery_addr_id']] = Fri[ZLT['RECORDS'][i]['delivery_addr_id']]
ZLT_CK.update(temp)
for i in range(1, len(ZLT['RECORDS'])):
temp = {}
if ZLT['RECORDS'][i]['delivery_addr_id'] not in ZLT_CK.keys():
if i in Sat:
temp[ZLT['RECORDS'][i]['delivery_addr_id']] = Sat[ZLT['RECORDS'][i]['delivery_addr_id']]
ZLT_CK.update(temp)
for i in range(1, len(ZLT['RECORDS'])):
temp = {}
if ZLT['RECORDS'][i]['delivery_addr_id'] not in ZLT_CK.keys():
if ZLT['RECORDS'][i]['delivery_addr_id'] in Sun:
temp[ZLT['RECORDS'][i]['delivery_addr_id']] = Sun[ZLT['RECORDS'][i]['delivery_addr_id']]
ZLT_CK.update(temp)
print(len(ZLT_CK))
11.10
论文里面用到的代码 画图
# coding=utf-8
import matplotlib.pyplot as plt
import gurobipy as gp
from gurobipy import *
import json
import math
# ===================================================================================
# 通过经纬度测算距离
# INPUTS: A点经度 A点纬度 B点经度 B点纬度
# OUTPUTS: 返回值为考虑球面的两点直线距离,单位：km
# ===================================================================================
def distance(lng_A, lat_A, lng_B, lat_B):
# 地球半径
R = 6371.004
pi = 3.141592654
Mlng_A = lng_A
Mlat_A = 90 - lat_A
Mlng_B = lng_B
Mlat_B = 90 - lat_B
C = math.sin(Mlat_A * pi / 180) * math.sin(Mlat_B * pi / 180) * math.cos((Mlng_A - Mlng_B) * pi / 180) + math.cos(
Mlat_A * pi / 180) * math.cos(Mlat_B * pi / 180)
if C >= 1:
C = 1
Distance = R * math.acos(C)
return Distance
# =================================================================================
# 将以str作为key值的字典转化为以int为key值的字典,并只保留前n个key-value对
# 输入：dict_str, n
# 输出：dict_int
# ===================================================================================
def Key_StrToInt(dict_str, n):
dict_int = {}
Keys = dict_str.keys()
Keys = list(Keys)
Keys = Keys[0:n]
for i in Keys:
i = int(i)
temp = {
i: dict_str[str(i)]
}
dict_int.update(temp)
return dict_int
# ==================读入以str为key值的字典=====================
filename1 = (r'E:\PycharmProjects\Project1\MonList.json')
Mon_str = json.load(open(filename1))
filename2 = (r'E:\PycharmProjects\Project1\TuesList.json')
Tues_str = json.load(open(filename2))
filename3 = (r'E:\PycharmProjects\Project1\WedList.json')
Wed_str = json.load(open(filename3))
filename4 = (r'E:\PycharmProjects\Project1\ThurList.json')
Thur_str = json.load(open(filename4))
filename5 = (r'E:\PycharmProjects\Project1\FriList.json')
Fri_str = json.load(open(filename5))
filename6 = (r'E:\PycharmProjects\Project1\SatList.json')
Sat_str = json.load(open(filename6))
filename7 = (r'E:\PycharmProjects\Project1\SunList.json')
Sun_str = json.load(open(filename7))
# ===================================================================================
n = 110 # 定义取每个场景的多少个数据
# 把字典转化成以int数字为key值的字典
Mon_int = Key_StrToInt(Mon_str, n)
Tues_int = Key_StrToInt(Tues_str, n)
Wed_int = Key_StrToInt(Wed_str, n)
Thur_int = Key_StrToInt(Thur_str, n)
Fri_int = Key_StrToInt(Fri_str, n)
Sat_int = Key_StrToInt(Sat_str, n)
Sun_int = Key_StrToInt(Sun_str, n)
U1, amount1, lng1, lat1 = gp.multidict(Mon_int)
U2, amount2, lng2, lat2 = gp.multidict(Tues_int)
U3, amount3, lng3, lat3 = gp.multidict(Wed_int)
U4, amount4, lng4, lat4 = gp.multidict(Thur_int)
U5, amount5, lng5, lat5 = gp.multidict(Fri_int)
U6, amount6, lng6, lat6 = gp.multidict(Sat_int)
U7, amount7, lng7, lat7 = gp.multidict(Sun_int)
capacity = 3000
FixedCost = 500000
T = 20
H = 365*5
# 求两个集合的并集
U1 = list(U1)
U2 = list(U2)
U3 = list(U3)
U4 = list(U4)
U5 = list(U5)
U6 = list(U6)
U7 = list(U7)
V = list(set(U1).union(U2, U3, U4, U5, U6, U7))
# T = [T for h in range(len(V))]
Z = [capacity for i in V]
C = [FixedCost for c in V]
# 把并集之后的数据写入lng[],lat[]，作为仓库潜在选址集合坐标
lng = {}
lat = {}
for k in V:
if k in U1:
a = lng1[k]
b = lat1[k]
elif k in U2:
a = lng2[k]
b = lat2[k]
elif k in U3:
a = lng3[k]
b = lat3[k]
elif k in U4:
a = lng4[k]
b = lat4[k]
elif k in U5:
a = lng5[k]
b = lat5[k]
elif k in U6:
a = lng6[k]
b = lat6[k]
elif k in U7:
a = lng7[k]
b = lat7[k]
x = {
k: a
}
y = {
k: b
}
lng.update(x)
lat.update(y)
# 计算各场景下距离矩阵D_jk
D1 = {
(i, j): distance(lng1[i], lat1[i], lng[j], lat[j]) for i in U1 for j in V
}
D2 = {
(i, j): distance(lng2[i], lat2[i], lng[j], lat[j]) for i in U2 for j in V
}
D3 = {
(i, j): distance(lng3[i], lat3[i], lng[j], lat[j]) for i in U3 for j in V
}
D4 = {
(i, j): distance(lng4[i], lat4[i], lng[j], lat[j]) for i in U4 for j in V
}
D5 = {
(i, j): distance(lng5[i], lat5[i], lng[j], lat[j]) for i in U5 for j in V
}
D6 = {
(i, j): distance(lng6[i], lat6[i], lng[j], lat[j]) for i in U6 for j in V
}
D7 = {
(i, j): distance(lng7[i], lat7[i], lng[j], lat[j]) for i in U7 for j in V
}
# 创建解集下标
x1 = [(i, j) for i in U1 for j in V]
x2 = [(i, j) for i in U2 for j in V]
x3 = [(i, j) for i in U3 for j in V]
x4 = [(i, j) for i in U4 for j in V]
x5 = [(i, j) for i in U5 for j in V]
x6 = [(i, j) for i in U6 for j in V]
x7 = [(i, j) for i in U7 for j in V]
x1 = tuplelist(x1)
x2 = tuplelist(x2)
x3 = tuplelist(x3)
x4 = tuplelist(x4)
x5 = tuplelist(x5)
x6 = tuplelist(x6)
x7 = tuplelist(x7)
V = tuplelist(V)
# Model
SM = gp.Model("Scenario")
# 添加变量
y = SM.addVars(V,
vtype=GRB.BINARY,
obj=C,
name="y")
x1 = SM.addVars(x1, vtype=GRB.BINARY, name="x1")
x2 = SM.addVars(x2, vtype=GRB.BINARY, name="x2")
x3 = SM.addVars(x3, vtype=GRB.BINARY, name="x3")
x4 = SM.addVars(x4, vtype=GRB.BINARY, name="x4")
x5 = SM.addVars(x5, vtype=GRB.BINARY, name="x5")
x6 = SM.addVars(x6, vtype=GRB.BINARY, name="x6")
x7 = SM.addVars(x7, vtype=GRB.BINARY, name="x7")
SM.update()
# 设置目标函数
SM.setObjective(
sum(y[j] * FixedCost for j in V)
+ sum(H*x1[i, j] * D1[i, j] for i in U1 for j in V)
+ sum(H*x2[i, j] * D2[i, j] for i in U2 for j in V)
+ sum(H*x3[i, j] * D3[i, j] for i in U3 for j in V)
+ sum(H*x4[i, j] * D4[i, j] for i in U4 for j in V)
+ sum(H*x5[i, j] * D5[i, j] for i in U5 for j in V)
+ sum(H*x6[i, j] * D6[i, j] for i in U6 for j in V)
+ sum(H*x7[i, j] * D7[i, j] for i in U7 for j in V),
GRB.MINIMIZE)
# 添加约束
SM.addConstrs(T >= D1[i, j] * x1[i, j] for i in U1 for j in V)
SM.addConstrs(T >= D2[i, j] * x2[i, j] for i in U2 for j in V)
SM.addConstrs(T >= D3[i, j] * x3[i, j] for i in U3 for j in V)
SM.addConstrs(T >= D4[i, j] * x4[i, j] for i in U4 for j in V)
SM.addConstrs(T >= D5[i, j] * x5[i, j] for i in U5 for j in V)
SM.addConstrs(T >= D6[i, j] * x6[i, j] for i in U6 for j in V)
SM.addConstrs(T >= D7[i, j] * x7[i, j] for i in U7 for j in V)
SM.addConstrs(sum(x1[i, j] for j in V) == 1 for i in U1)
SM.addConstrs(sum(x2[i, j] for j in V) == 1 for i in U2)
SM.addConstrs(sum(x3[i, j] for j in V) == 1 for i in U3)
SM.addConstrs(sum(x4[i, j] for j in V) == 1 for i in U4)
SM.addConstrs(sum(x5[i, j] for j in V) == 1 for i in U5)
SM.addConstrs(sum(x6[i, j] for j in V) == 1 for i in U6)
SM.addConstrs(sum(x7[i, j] for j in V) == 1 for i in U7)
SM.addConstrs(sum(x1[i, j] * amount1[i] for i in U1) <= capacity for j in V)
SM.addConstrs(sum(x2[i, j] * amount2[i] for i in U2) <= capacity for j in V)
SM.addConstrs(sum(x3[i, j] * amount3[i] for i in U3) <= capacity for j in V)
SM.addConstrs(sum(x4[i, j] * amount4[i] for i in U4) <= capacity for j in V)
SM.addConstrs(sum(x5[i, j] * amount5[i] for i in U5) <= capacity for j in V)
SM.addConstrs(sum(x6[i, j] * amount6[i] for i in U6) <= capacity for j in V)
SM.addConstrs(sum(x7[i, j] * amount7[i] for i in U7) <= capacity for j in V)
# SM.addConstr(sum(y[j] for j in V) == 1)
for i in U1:
for j in V:
SM.addGenConstrIndicator(x1[i, j], 1, y[j] == 1, name="ind1")
for i in U2:
for j in V:
SM.addGenConstrIndicator(x2[i, j], 1, y[j] == 1, name="ind2")
for i in U3:
for j in V:
SM.addGenConstrIndicator(x3[i, j], 1, y[j] == 1, name="ind3")
for i in U4:
for j in V:
SM.addGenConstrIndicator(x4[i, j], 1, y[j] == 1, name="ind4")
for i in U5:
for j in V:
SM.addGenConstrIndicator(x5[i, j], 1, y[j] == 1, name="ind5")
for i in U6:
for j in V:
SM.addGenConstrIndicator(x6[i, j], 1, y[j] == 1, name="ind6")
for i in U7:
for j in V:
SM.addGenConstrIndicator(x7[i, j], 1, y[j] == 1, name="ind7")
# Save model
SM.write('facilityPY.lp')
# Solve
SM.optimize()
print(SM.objVal)
x1 = SM.getAttr('x', x1)
x2 = SM.getAttr('x', x2)
x3 = SM.getAttr('x', x3)
x4 = SM.getAttr('x', x4)
x5 = SM.getAttr('x', x5)
x6 = SM.getAttr('x', x6)
x7 = SM.getAttr('x', x7)
y = SM.getAttr('x', y)
for j in V:
if y[j] == 1:
print('y', j, y[j])
flag = 1
for i in U1:
for j in V:
if 0.9 <= x1[i, j] <= 1.1:
if 0.9 < flag < 1.1:
# plt.plot([lng1[i], lng[j]], [lat1[i], lat[j]], 'b')
s1 = plt.scatter(lng1[i], lat1[i], color='r')
s2 = plt.scatter(lng[j], lat[j], color='g')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
else:
# plt.plot([lng1[i], lng[j]], [lat1[i], lat[j]], 'b')
plt.scatter(lng1[i], lat1[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('(a)')
plt.show()
for i in U2:
for j in V:
if 0.9 <= x2[i, j] <= 1.1:
if 1.9 < flag < 2.1:
# plt.plot([lng2[i], lng[j]], [lat2[i], lat[j]])
plt.scatter(lng2[i], lat2[i], color='r', label='merchant')
plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend(loc="upper left")
flag = flag + 1
else:
# plt.plot([lng2[i], lng[j]], [lat2[i], lat[j]], 'b')
plt.scatter(lng2[i], lat2[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('(b)')
plt.show()
for i in U3:
for j in V:
if 0.9 <= x3[i, j] <= 1.1:
if 2.9 < flag < 3.1:
# plt.plot([lng3[i], lng[j]], [lat3[i], lat[j]], 'b')
plt.scatter(lng3[i], lat3[i], color='r', label='merchant')
plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend(loc="upper left")
flag = flag + 1
else:
# plt.plot([lng3[i], lng[j]], [lat3[i], lat[j]], 'b')
plt.scatter(lng3[i], lat3[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('(c)')
plt.show()
for i in U4:
for j in V:
if 0.9 <= x4[i, j] <= 1.1:
if 3.9 < flag < 4.1:
s1 = plt.scatter(lng4[i], lat4[i], color='r', label='merchant')
# plt.plot([lng4[i], lng[j]], [lat4[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
else:
# plt.plot([lng4[i], lng[j]], [lat4[i], lat[j]], 'b')
plt.scatter(lng4[i], lat4[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('(d)')
plt.show()
for i in U5:
for j in V:
if 0.9 <= x5[i, j] <= 1.1:
if 4.9 < flag < 5.1:
s1 = plt.scatter(lng5[i], lat5[i], color='r', label='merchant')
# plt.plot([lng5[i], lng[j]], [lat5[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
else:
# plt.plot([lng5[i], lng[j]], [lat5[i], lat[j]], 'b')
plt.scatter(lng5[i], lat5[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('(e)')
plt.show()
for i in U6:
for j in V:
if 0.9 <= x6[i, j] <= 1.1:
if 5.9 < flag < 6.1:
s1 = plt.scatter(lng6[i], lat6[i], color='r', label='merchant')
# plt.plot([lng6[i], lng[j]], [lat6[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
else:
# plt.plot([lng6[i], lng[j]], [lat6[i], lat[j]], 'b')
plt.scatter(lng6[i], lat6[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('(f)')
plt.show()
for i in U7:
for j in V:
if 0.9 <= x7[i, j] <= 1.1:
if 6.9 < flag < 7.1:
s1 = plt.scatter(lng7[i], lat7[i], color='r', label='merchant')
# plt.plot([lng7[i], lng[j]], [lat7[i], lat[j]], 'b')
s2 = plt.scatter(lng[j], lat[j], color='g', label='FDC')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
else:
# plt.plot([lng7[i], lng[j]], [lat7[i], lat[j]], color='b')
plt.scatter(lng7[i], lat7[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.title('(g)')
plt.show()
论文内500个点的代码
# coding=utf-8
import matplotlib.pyplot as plt
import gurobipy as gp
from gurobipy import *
import json
import math
# ===================================================================================
# 通过经纬度测算距离
# INPUTS: A点经度 A点纬度 B点经度 B点纬度
# OUTPUTS: 返回值为考虑球面的两点直线距离,单位：km
# ===================================================================================
def distance(lng_A, lat_A, lng_B, lat_B):
# 地球半径
R = 6371.004
pi = 3.141592654
Mlng_A = lng_A
Mlat_A = 90 - lat_A
Mlng_B = lng_B
Mlat_B = 90 - lat_B
C = math.sin(Mlat_A * pi / 180) * math.sin(Mlat_B * pi / 180) * math.cos((Mlng_A - Mlng_B) * pi / 180) + math.cos(
Mlat_A * pi / 180) * math.cos(Mlat_B * pi / 180)
if C >= 1:
C = 1
Distance = R * math.acos(C)
return Distance
# =================================================================================
# 将以str作为key值的字典转化为以int为key值的字典,并只保留前n个key-value对
# 输入：dict_str, n
# 输出：dict_int
# ===================================================================================
def Key_StrToInt(dict_str, n):
dict_int = {}
Keys = dict_str.keys()
Keys = list(Keys)
Keys = Keys[0:n]
for i in Keys:
i = int(i)
temp = {
i: dict_str[str(i)]
}
dict_int.update(temp)
return dict_int
# ==================读入以str为key值的字典=====================
filename1 = (r'E:\PycharmProjects\Project1\MonList.json')
Mon_str = json.load(open(filename1))
filename2 = (r'E:\PycharmProjects\Project1\TuesList.json')
Tues_str = json.load(open(filename2))
filename3 = (r'E:\PycharmProjects\Project1\WedList.json')
Wed_str = json.load(open(filename3))
filename4 = (r'E:\PycharmProjects\Project1\ThurList.json')
Thur_str = json.load(open(filename4))
filename5 = (r'E:\PycharmProjects\Project1\FriList.json')
Fri_str = json.load(open(filename5))
filename6 = (r'E:\PycharmProjects\Project1\SatList.json')
Sat_str = json.load(open(filename6))
filename7 = (r'E:\PycharmProjects\Project1\SunList.json')
Sun_str = json.load(open(filename7))
# ===================================================================================
n = 500 # 定义取每个场景的多少个数据
# 把字典转化成以int数字为key值的字典
Mon_int = Key_StrToInt(Mon_str, n)
U1, amount1, lng1, lat1 = gp.multidict(Mon_int)
capacity = 3000
FixedCost = 500000
T = 20
H = 365*5
# 求两个集合的并集
U1 = list(U1)
V = list(U1)
# T = [T for h in range(len(V))]
Z = [capacity for i in V]
C = [FixedCost for c in V]
# 把并集之后的数据写入lng[],lat[]，作为仓库潜在选址集合坐标
lng = {}
lat = {}
for k in V:
if k in U1:
a = lng1[k]
b = lat1[k]
x = {
k: a
}
y = {
k: b
}
lng.update(x)
lat.update(y)
# 计算各场景下距离矩阵D_jk
D1 = {
(i, j): distance(lng1[i], lat1[i], lng[j], lat[j]) for i in U1 for j in V
}
# 创建解集下标
x1 = [(i, j) for i in U1 for j in V]
x1 = tuplelist(x1)
V = tuplelist(V)
# Model
SM = gp.Model("Scenario")
# 添加变量
y = SM.addVars(V,
vtype=GRB.BINARY,
obj=C,
name="y")
x1 = SM.addVars(x1, vtype=GRB.BINARY, name="x1")
SM.update()
# 设置目标函数
SM.setObjective(
sum(y[j] * FixedCost for j in V)
+ sum(H*x1[i, j] * D1[i, j] for i in U1 for j in V),
GRB.MINIMIZE)
# 添加约束
SM.addConstrs(T >= D1[i, j] * x1[i, j] for i in U1 for j in V)
SM.addConstrs(sum(x1[i, j] for j in V) == 1 for i in U1)
SM.addConstrs(sum(x1[i, j] * amount1[i] for i in U1) <= capacity for j in V)
# SM.addConstr(sum(y[j] for j in V) == 1)
for i in U1:
for j in V:
SM.addGenConstrIndicator(x1[i, j], 1, y[j] == 1, name="ind1")
# Save model
SM.write('facilityPY.lp')
# Solve
SM.optimize()
print(SM.objVal)
x1 = SM.getAttr('x', x1)
y = SM.getAttr('x', y)
for j in V:
if y[j] == 1:
print('y', j, y[j])
flag = 1
for i in U1:
for j in V:
if 0.9 <= x1[i, j] <= 1.1:
if 0.9 < flag < 1.1:
# plt.plot([lng1[i], lng[j]], [lat1[i], lat[j]], 'b')
s1 = plt.scatter(lng1[i], lat1[i], color='r')
s2 = plt.scatter(lng[j], lat[j], color='g')
plt.legend((s1, s2), ('merchant', 'FDC'), loc='upper left')
flag = flag + 1
else:
# plt.plot([lng1[i], lng[j]], [lat1[i], lat[j]], 'b')
plt.scatter(lng1[i], lat1[i], color='r')
plt.scatter(lng[j], lat[j], color='g')
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(112.9, 114)
plt.ylim(22.4, 24)
plt.show()
关于html只显示一个点的问题
def file_content_alter(file_name, replace_str_dict):
file_data = ""
with open(file_name, "r", encoding="utf-8") as f:
for line in f:
for k, v in replace_str_dict.items():
line = line.replace(k, v)
file_data += line
with open(file_name, "w", encoding="utf-8") as f:
f.write(file_data)
file_content_alter(html_path, {"https://code.jquery.com/jquery-1.12.4.min.js": "https://libs.baidu.com/jquery/2.0.0/jquery.min.js",
"https://rawcdn.githack.com/python-visualization/folium/master/folium/templates": "."})
479
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
