方法1:精确重心法
import numpy as np
import pandas as pd
import math as m
from math import radians, cos, sin, asin, sqrt
#球面距离函数
def geodistance(lng1,lat1,lng2,lat2):
#lng1,lat1,lng2,lat2 = (120.12802999999997,30.28708,115.86572000000001,28.7427)
lng1, lat1, lng2, lat2 = map(radians, [float(lng1), float(lat1), float(lng2), float(lat2)]) # 经纬度转换成弧度
dlon=lng2-lng1
dlat=lat2-lat1
a=sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
distance=2*asin(sqrt(a))*6371*1000 # 地球平均半径,6371km
distance=round(distance/1000,3)
return distance
data = pd.read_excel("data1212.xls")#W,C是费率和货量,对距离进行加权,转换成运费成本。
#精确重心法迭代
WC = np.array(data['W']) * np.array(data['C'])
WCX = (np.array(data['X_c']) * WC).sum()
WCY = (np.array(data['Y_c']) * WC).sum()
x0 = WCX / WC.sum()#初始等效重心
y0 = WCY / WC.sum()#初始等效重心
d_j=[]
for m in range(len(np.array(data['X_c']))):
d_j.append(geodistance(np.array(data['X_c'])[m],np.array(data['Y_c'])[m],x0,y0))#球面距离
#d_j.append(((np.array(data['X_c'])[m] - x0)**2 + (np.array(data['Y_c'])[m]-y0)**2)**0.5) #欧式距离
T= (WC * d_j).sum()#算总成本
print('重心法初始选点位置:({},{})'.format(x0,y0))
print('总费用T0:{}'.format(T))
#选址迭代
for i in range(10):
WC_j = WC/d_j
WCX_j = ((np.array(data['X_c']) * WC)/d_j).sum()
WCY_j = ((np.array(data['Y_c']) * WC)/d_j).sum()
x = WCX_j / WC_j.sum()
y = WCY_j / WC_j.sum()
d_j=[]
for j in range(len(np.array(data['X_c']))):
d_j.append(geodistance(np.array(data['X_c'])[j],np.array(data['Y_c'])[j],x,y))#球面距离
#d_j.append(((np.array(data['X_c'])[j] - x)**2 + (np.array(data['Y_c'])[j]-y)**2)**0.5) #欧式距离
T = (WC * d_j).sum()
print('经{}次迭代后选址点位置:({},{})'.format(i+1,x,y))
print('总费用T{}:{}'.format(i+1,T))
输出结果:
data数据样式:最终选址确定在经纬度:113.28,23.02 总运输成本:42557.4182
方法2:遗传算法选址
#搭建遗传算法问题框架,单个设施选址
import numpy as np
import geatpy as ea
#选址目标函数
def km1(x,y):
K=0
for i in range(len(data)):
K+=geodistance(x,y,data.X_c.values[i],data.Y_c.values[i])*data.W.values[i]*data.C.values[i]
return K
class MyProblem(ea.Problem): # 继承Problem父类
def __init__(self):
name = 'MyProblem' # 初始化name(函数名称,可以随意设置)
M = 1 # 初始化M(目标维数)
maxormins = [1] # 初始化maxormins(目标最小最大化标记列表,1:最小化该目标;-1:最大化该目标)
Dim = 2 # 初始化Dim(决策变量维数)
varTypes = [0] * Dim # 初始化varTypes(决策变量的类型,元素为0表示对应的变量是连续的;1表示是离散的)
lb = [100, 20] # 决策变量下界
ub = [120, 26] # 决策变量上界
lbin = [1,1] # 决策变量下边界(0表示不包含该变量的下边界,1表示包含)
ubin = [1,1] # 决策变量上边界(0表示不包含该变量的上边界,1表示包含)
# 调用父类构造方法完成实例化
ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)
def evalVars(self, Vars): # 目标函数
x1 = Vars[:, [0]]
x2 = Vars[:, [1]]
f = km1(x1[0],x2[0])
for mm in range(1,10):#跟后面种群数相同
f=np.vstack((f,km1(x1[mm],x2[mm])))
###设置约束范围
CV = np.hstack([x1- 130,
x2 - 25])
return f, CV
def calReferObjV(self): # 设定目标数参考值(本问题目标函数参考值设定为理论最优值)可不设置,对结果无影响
referenceObjV = np.array([[300000]])
return referenceObjV
#遗传算法求解
if __name__ == '__main__':
# 实例化问题对象
problem = MyProblem()
# 构建算法
algorithm = ea.soea_DE_rand_1_bin_templet(problem,
ea.Population(Encoding='RI', NIND=10),
MAXGEN=100, # 最大进化代数。
logTras=10) # 表示每隔多少代记录一次日志信息,0表示不记录。
algorithm.mutOper.F = 0.5 # 差分进化中的参数F
algorithm.recOper.XOVR = 0.7 # 重组概率
# 求解
res = ea.optimize(algorithm, verbose=True, drawing=1, outputMsg=True, drawLog=False, saveFlag=True)
print(res)
微信公众号后台回复
加群:加入全球华人OR|AI|DS社区硕博微信学术群
资料:免费获得大量运筹学相关学习资料
人才库:加入运筹精英人才库,获得独家职位推荐
电子书:免费获取平台小编独家创作的优化理论、运筹实践和数据科学电子书,持续更新中ing...
加入我们:加入「运筹OR帷幄」,参与内容创作平台运营
知识星球:加入「运筹OR帷幄」数据算法社区,免费参与每周「领读计划」、「行业inTalk」、「OR会客厅」等直播活动,与数百位签约大V进行在线交流
文章须知
文章作者:用户007
微信编辑:疑疑
文章转载自『Python学习杂记』公众号,原文链接:选址问题(一)-精确重心法和遗传算法
关注我们
FOLLOW US