p_center解决方案.doc

问题：P-Center问题P-Center问题摘要问题：从一个点集合中选择p个点作为中心点，并把其他分配到某个选择的点，使得所有点到其对应的中心点的距离加起来最小。针对问题，分析得出p-center问题实质为选址问题。其研究背景为工厂的选址，属于运筹学中经典问题之一。运用智能算法中模拟退火随机局域搜索算法进行求解最小值。具体步骤为：由于题目中未提供点集合，所以首先通过文献查阅1和生活实际得到可靠的二维数据点分布，即表示一个点集合，存储方式为文件存储（data.txt）；其次加载点集合数据，采用模拟退火算法随机局域搜索算法2进行处理：1） 初始化：中心点个数，温度，温度缩减系数，最大迭代次数。解释：由于P-Center问题是以工厂的选址问题，加上编写的二维数据的分布情况，所以建立工厂的数量为事先已知条件，即；初试温度的设置是影响搜索性能重要因素之一。初始化温度高，则搜索到全局最优值的可能性大，但计算时间大幅度增加；反之，缩短了计算时间，但性能并不优越。所以采用多次实验的方式确定温度；为了保证较大的搜索空间，所以让系数接近于1；经过多次实验，确定迭代次数，此时结果较为理想。2） 迭代次，循环第三步到第四步。3）从临域中选择最佳的解，即确定最优值：产生新值；增量；若，则接受 作为新解，否者以概率接受作为新解。4）逐渐减少，且，最后跳至第二步。5）得到距离最小值然后通过模拟退火局部搜索算法，得迭代情况为：最后通过模拟退火局部搜索算法，得出分配图为：得出四个粗五角星为各自的中心点，其中颜色相同的属于各自颜色的中心点，即相同颜色距离各自中心点最短。通过Python得出最近距离为：102.401974373问题扩充：针对P-Center问题，还可以通过k-means聚类算法3进行解决，得到与最近搜索算法同样的结果。关键词： P-Center选址问题 模拟退火随机局域搜索算法 K-Means聚类算法目录P-Center问题2摘要21问题重述42数据预处理42.1数据来源42.2数据预处理方法42.3数据选取参考原则43问题分析43.1问题44 问题假设55 符号说明56 模型的建立与求解56.1解法一56.1.1模拟退火随机局部搜索算法56.1.2算法求解76.2解法二86.2.1K-Means聚类算法86.2.2算法求解87模型的评价98 参考文献9附录10附录一 模拟退火随机局域搜索算法Python代码10附录二 聚类算法Python代码12附录三 迭代次数与目标函数关系15附录四 结论图201问题重述选址问题是运筹学中经典的问题之一。经典的选址问题包括连续选址问题、离散选址问题、P-Median问题以及P-Center问题等。该问题属于P-Center问题，从一个点集合中选择p个点作为中心点，并把其他点分配到某个选择的点，使得所有点到其对应中心点的距离加起来最小。2数据预处理注：数据存储形式为：data.txt，可在附件一中查看2.1数据来源（1） 文献查阅（2） 生活实际2.2数据预处理方法（1） 数据选取：去除无效数据以及噪声数据。（2） 数据整合：将若干个数据库整合成一个数据存储结构。（3） 数据替代：将杂乱数据替代成易处理数据。（4） 数据变换：将原始数据转换为适合任务定价的形式。2.3数据选取参考原则（1） 统一数据源编码。（2） 清除唯一属性。（3） 清除重复属性。（4） 清除可忽略字段。（5） 进一步处理：清除异常数据，去除附带噪声数据，填充空缺值、丢失值。3问题分析3.1问题首先，通过文献查阅和生活实际得到可靠的二维数据点分布，将此二维分布作为数据点集合，然后通过模拟退火最近搜索算法4进行迭代优化，最终得到所有点到其中心点的距离。4 问题假设1. 假设根据数据特征或者政策（例如工厂政策）确定，即已经确定P-Center中的p中心个数。2.假设点集合为二维集合，不包括任何三维或者多维信息。5 符号说明符号说明可行解迭代更新后的解概率条件下更新优化目标函数模拟退火温度值中心点个数温度缩减系数最大迭代次数6 模型的建立与求解选址问题目前有多种求解方法，大致分为定性和定量两类：定性的方法主要是结合层次分析法和模糊综合法对各方案进行指标评价，最终得出最优选址；定量的方法主要是松弛算法和启发式算法以及两者的结合。而本题解决的是P-Center问题，即使用启发式算法-智能算法模拟退火随机局部搜索算法5进行解决。6.1解法一6.1.1模拟退火随机局部搜索算法模拟退火算法是一种贪心算法，但是其搜索过程引入了随机因素。在迭代更新可行解时，以一定概率来接受一个比当前解要差的解，因此有可能会跳出局部最优解，达到全局最优解，其搜索原理如下图：图1 模拟退火随机局部搜索算法示意图假定当前可行解为，迭代更新后的解为，那么对应的“能量差”定义为：以相应概率接受新值，此概率为：模拟退火随机局域搜索算法思路为：图2 模拟退火随机局域搜索算法思路图6.1.2算法求解注：详细Python程序见附录一与附件一根据算法思想，通过Python得到迭代与目标函数关系为：图3 模拟退火训练曲线图可以看出经过200次迭代，优化目标函数处于相对稳定状态。具体的目标函数与迭代次数对应关系见附录，以下为部分对应关系：图4 部分迭代次数与目标函数关系最终得出点集合分配图为：图5 点集合分配图得出四个粗五角星为各自的中心点，其中颜色相同的属于各自颜色的中心点，即相同颜色距离各自中心点最短。通过Python得出最近距离为：102.4019743736.2解法二通过智能算法中的模拟退火随机局域搜索算法可以得出相应的结论，为了检验以及去探索更多的解决方式，使用了聚类算法，以下为模型以及过程。具体的Python代码见附录二以及附件一。6.2.1K-Means聚类算法算法过程如下：1） 从N个点中随机选取p个点作为中心点2） 对剩余的点测量到其每个中心点的距离，并把它归到中心点的类3） 重新计算已经得到的各个类的中心点4） 迭代第二步、第三步直至新的中心点与原中心点相等或者小于指定阈值，算法结束6.2.2算法求解通过Python程序将通过文献查找的点集合数据data.txt进行聚类分析，得到与局域搜索算法类似的分配图，如下图所示：图6 聚类分析分配图此图的解释与模拟退火随机局域搜索算法相同，不再重复解释。同样可以得出所有点到对应中心点的距离最小为：102.401974373。7模型的评价7.1模型的优缺点模型的优点：（1）通过模拟退火随机局域搜索算法应用十分广泛，可以高效的解决NP完全问题（2）模拟退火随机搜索算法与K-Means聚类算法相互结合，相互印证，使得数据准确性得以保证。（3）通过模拟退火算法与K-Means算法得出的点集合分配图可以形象生动的得出二维的数据关系。模型的缺点：（1）模拟退火随机局域搜索算法初始化设置必须准确，要经过多次实验才能得到合适的初始化值。8 参考文献1CCEHC: An efficient local search algorithm for weighted partial maximum satisfiability,Artificial Intelligence,2017,2Local Search for Minimum Weight Dominating Set with Two-Level Configuration Checking and Frequency Based Scoring Function,Journal of Artificial Intelligence Research (JAIR),2017,3王守强. 多中心点聚类问题的随机算法D.山东大学,2010.4朱泓丞. 设施选址问题的研究与应用D.中国科学技术大学,2009.5蒋建林,徐进澎,文杰.基于单亲遗传模拟退火算法的顶点p-中心问题J.系统工程学报,2011,26(03):414-420.附录注：具体代码见附件附录一 模拟退火随机局域搜索算法Python代码import numpy as npimport matplotlib.pyplot as pltimport matplotlibimport scipy.spatial.distance as DSimport randomzandaoguangfont = matplotlib.font_manager.FontProperties(fname=a.TTF)def f(x): 目标函数，所有点到中心点的距离和 :param x: list 表示中心点的序号 :param dataSet: 数据集 :return: distances=DS.cdist(dataSet,dataSetx,:)#用于计算两个输入集合的距离 M=np.min(distances,axis=1)#得出附近点到中心点的最小值 return sum(M)#将最小值累加得到最小距离def find_next(x): 从邻域中选择最佳的解 :param x: :param dataSet: :return:最优值 ind=random.sample(range(p),1) x_new=list.copy(x) tmp=random.sample(range(dataSet.shape0),1)#dataSet.shape0返回行数 while tmp0 in x: tmp = random.sample(range(dataSet.shape0), 1) x_newint(ind0)=tmp0 if f(x_new)f(x): return x_new else: if random.random()np.exp(-(f(x_new)-f(x)/T):#模拟退火算法核心 return x_new else: return x#得到数据print (第一步：加载数据 ) dataSet=fileIn = open(data.txt)for line in fileIn.readlines(): lineArr = line.strip().split( ) dataSet.append(float(lineArr0), float(lineArr-1)dataSet=np.array(dataSet)p=4#设定问题的p中心参数的值T=1000#模拟退火的初始温度a=0.999#温度减少的系数，为了保证较大的搜索空间，所以非常接近于1Iter=200#最大迭代次数#模拟退火#param:history记录距离和最优值，方便比较最小值#param:best_x记录当前距离和最小值x=random.sample(range(dataSet.shape0),p)#随机产生初始解best_x=xhistory=for i in range(Iter): print(第%d次迭代：目标函数=%s%(i,f(x) x=find_next(x) T*=a if f(x)f(best_x): best_x=x history.append(f(best_x)print(距离和最小为%s%f(best_x)#模拟退火局部搜索训练曲线图#resul:t得出模拟退火的迭代次数与目标函数值的关系plt.figure(0)plt.plot(range(1,Iter+1),history,r-)plt.xlabel(迭代次数,fontproperties=zandaoguangfont)plt.ylabel(目标函数值,fontproperties=zandaoguangfont)plt.title(模拟退火训练曲线图,fontproperties=zandaoguangfont)#数据分配图#用红蓝黄绿等颜色代表p中心中的p类#其中经过搜索后的p个中心点用markersize=10的五边形表示#相同颜色的点相比较与相同颜色的加粗markersize=10的五边形，即相同颜色以相同颜色的加粗五边形为中心，此时距离最近。plt.figure(1)distances=DS.cdist(dataSet,dataSetbest_x,:)sorted_dist=np.argsort(distances,axis=1)color=r.,b.,y.,g.color0=rp,bp,yp,gpfor i in range(p): ind=sorted_dist:,0=i plt.plot(dataSetind,0,dataSetind,1,colori) plt.plot(dataSetbest_xi,0,dataSetbest_xi,1,color0i,markersize=10)plt.title(分配图,fontproperties=zandaoguangfont)plt.show()附录二 聚类算法Python代码from numpy import * import matplotlib.pyplot as plt import KMeans # 得到数据.print (第一步：加载数据. ) dataSet = fileIn = open(data.txt) for line in fileIn.readlines(): temp=lineArr = line.strip().split(t) temp.append(float(lineArr0)temp.append(float(lineArr1)dataSet.append(temp)fileIn.close() # 进行聚类.print (第二步：进行聚类. )dataSet = mat(dataSet) #设置为4个中心点k = 4 #调用KMeans文件中定义的kmeans方法。centroids, clusterAssment = KMeans.kmeans(dataSet, k) # 第三步：显示结果.#数据分配图#用不同种颜色代表p中心中的p类#其中经过搜索后的p个中心点用markersize=10的五边形表示#相同颜色的点相比较与相同颜色的加粗markersize=10的五边形，即相同颜色以相同颜色的加粗五边形为中心，此时距离最近。print (第三步：显示结果. )KMeans.showCluster(dataSet, k, centroids, clusterAssment)Kmeans.py# # kmeans: k-means cluster # Author :ZanDaoguang# Date :2018/03/09# HomePage : # Email : # from numpy import * import time import matplotlib.pyplot as plt # calculate Euclidean distance def euclDistance(vector1, vector2): return sqrt(sum(power(vector2 - vector1, 2) #求这两个矩阵的距离，vector1、2均为矩阵 # init centroids with random samples #在样本集中随机选取k个样本点作为初始质心def initCentroids(dataSet, k): numSamples, dim = dataSet.shape #矩阵的行数、列数 centroids = zeros(k, dim) #感觉要不要你都可以for i in range(k): index = int(random.uniform(0, numSamples) #随机产生一个浮点数，然后将其转化为int型centroidsi, : = dataSetindex, : return centroids # k-means cluster #dataSet为一个矩阵#k为将dataSet矩阵中的样本分成k个类 def kmeans(dataSet, k): numSamples = dataSet.shape0 #读取矩阵dataSet的第一维度的长度,即获得有多少个样本数据 # first column stores which cluster this sample belongs to, # second column stores the error between this sample and its centroid clusterAssment = mat(zeros(numSamples, 2) #得到一个N*2的零矩阵clusterChanged = True # step 1: init centroids centroids = initCentroids(dataSet, k) #在样本集中随机选取k个样本点作为初始质心 while clusterChanged: clusterChanged = False # for each sample for i in range(numSamples): #rangeminDist = 100000.0 minIndex = 0 # for each centroid # step 2: find the centroid who is closest #计算每个样本点与质点之间的距离，将其归内到距离最小的那一簇for j in range(k): distance = euclDistance(centroidsj, :, dataSeti, :) if distance minDist: minDist = distance minIndex = j # step 3: update its cluster #k个簇里面与第i个样本距离最小的的标号和距离保存在clusterAssment中#若所有的样本不在变化，则退出while循环if clusterAssmenti, 0 != minIndex: clusterChanged = True clusterAssmenti, : = minIndex, minDist*2 #两个*表示的是minDist的平方 # step 4: update centroids for j in range(k): #clusterAssment:,0.A=j是找出矩阵clusterAssment中第一列元素中等于j的行的下标，返回的是一个以array的列表，第一个array为等于j的下标pointsInCluster = dataSetnonzero(clusterAssment:, 0.A = j)0 #将dataSet矩阵中相对应的样本提取出来 centroidsj, : = mean(pointsInCluster, axis = 0) #计算标注为j的所有样本的平均值 print ( 聚类完成！) return centroids, clusterAssment # show your cluster only available with 2-D data #centroids为k个类别，其中保存着每个类别的质心#clusterAssment为样本的标记，第一列为此样本的类别号，第二列为到此类别质心的距离 def showCluster(dataSet, k, centroids, clusterAssment): numSamples, dim = dataSet.shape if dim != 2: print (Sorry! I can not draw because the dimension of your data is not 2!) return 1 mark = or, ob, og, ok, r, +r, sr, dr, len(mark): print (对不起，您输入的k太大了，请联系管理员山东科技大学昝道广) return 1 # draw all samples for i in range(numSamples): markIndex = int(clusterAssmenti, 0) #为样本指定颜色plt.plot(dataSeti, 0, dataSeti, 1, markmarkIndex) mark = Dr, Db, Dg, Dk, b, +b, sb, db, b, pb # draw the centroids for i in range(k): plt.plot(centroidsi, 0, centroidsi, 1, marki, markersize = 10)plt.title(Distribution diagram) plt.show()附录三 迭代次数与目标函数关系第0次迭代：目标函数=179.115046893第1次迭代：目标函数=205.447603458第2次迭代：目标函数=179.714990497第3次迭代：目标函数=190.540746471第4次迭代：目标函数=186.205430169第5次迭代：目标函数=255.113966307第6次迭代：目标函数=237.264836994第7次迭代：目标函数=360.246463689第8次迭代：目标函数=235.13584001第9次迭代：目标函数=245.721007671第10次迭代：目标函数=247.472895118第11次迭代：目标函数=273.485021319第12次迭代：目标函数=280.851747779第13次迭代：目标函数=280.308274938第14次迭代：目标函数=198.743816875第15次迭代：目标函数=195.036904207第16次迭代：目标函数=277.824032017第17次迭代：目标函数=292.981152147第18次迭代：目标函数=190.984905364第19次迭代：目标函数=191.252223971第20次迭代：目标函数=123.059185807第21次迭代：目标函数=190.533091518第22次迭代：目标函数=189.302660621第23次迭代：目标函数=187.360273683第24次迭代：目标函数=187.360273683第25次迭代：目标函数=216.753632115第26次迭代：目标函数=264.196159321第27次迭代：目标函数=264.296840503第28次迭代：目标函数=181.742051739第29次迭代：目标函数=122.674228117第30次迭代：目标函数=186.925795284第31次迭代：目标函数=256.518821535第32次迭代：目标函数=237.353645794第33次迭代：目标函数=191.156290007第34次迭代：目标函数=216.105125507第35次迭代：目标函数=147.578774857第36次迭代：目标函数=222.828821333第37次迭代：目标函数=215.525456567第38次迭代：目标函数=204.837254851第39次迭代：目标函数=197.185551202第40次迭代：目标函数=188.580114581第41次迭代：目标函数=265.711905553第42次迭代：目标函数=163.812778837第43次迭代：目标函数=162.068177064第44次迭代：目标函数=164.368700208第45次迭代：目标函数=161.816009847第46次迭代：目标函数=204.73869345第47次迭代：目标函数=162.769100846第48次迭代：目标函数=178.51942545第49次迭代：目标函数=173.561807727第50次迭代：目标函数=169.799603435第51次迭代：目标函数=168.651134516第52次迭代：目标函数=245.443837937第53次迭代：目标函数=169.449644817第54次迭代：目标函数=116.628068242第55次迭代：目标函数=171.256461075第56次迭代：目标函数=167.716292084第57次迭代：目标函数=163.959218108第58次迭代：目标函数=167.828227585第59次迭代：目标函数=186.81606983第60次迭代：目标函数=173.885593243第61次迭代：目标函数=171.185562964第62次迭代：目标函数=165.620351515第63次迭代：目标函数=171.009592894第64次迭代：目标函数=258.202349371第65次迭代：目标函数=300.310149921第66次迭代：目标函数=245.514351498第67次迭代：目标函数=158.696269997第68次迭代：目标函数=181.751581539第69次迭代：目标函数=190.838797732第70次迭代：目标函数=245.514351498第71次迭代：目标函数=187.41644075第72次迭代：目标函数=185.529583413第73次迭代：目标函数=192.972377185第74次迭代：目标函数=203.792334694第75次迭代：目标函数=179.422114776第76次迭代：目标函数=176.565275563第77次迭代：目标函数=177.060338213第78次迭代：目标函数=261.856371227第79次迭代：目标函数=186.265589017第80次迭代：目标函数=111.762469877第81次迭代：目标函数=164.693111801第82次迭代：目标函数=199.237534931第83次迭代：目标函数=201.605152809第84次迭代：目标函数=200.117346554第85次迭代：目标函数=226.507916216第86次迭代：目标函数=162.047130232第87次迭代：目标函数=180.016930471第88次迭代：目标函数=186.524458902第89次迭代：目标函数=184.228403531第90次迭代：目标函数=175.548638845第91次迭代：目标函数=134.209706505第92次迭代：目标函数=188.793411141第93次迭代：目标函数=178.706438929第94次迭代：目标函数=179.685911322第95次迭代：目标函数=244.590852758第96次迭代：目标函数=170.495284118第97次迭代：目标函数=251.458815477第98次迭代：目标函数=170.702193198第99次迭代：目标函数=263.478662358第100次迭代：目标函数=266.230685027第101次迭代：目标函数=233.806726474第102次迭代：目标函数=236.472566702第103次迭代：目标函数=225.158342676第104次迭代：目标函数=256.998374944第105次迭代：目标函数=253.877494761第106次迭代：目标函数=192.000533947第107次迭代：目标函数=178.526999848第108次迭代：目标函数=181.81273764第109次迭代：目标函数=209.9163263第110次迭代：目标函数=213.281909621第111次迭代：目标函数=188.821519052第112次迭代：目标函数=123.867798998第113次迭代：目标函数=105.140919713第114次迭代：目标函数=170.167350969第115次迭代：目标函数=162.045708686第116次迭代：目标函数=202.011426338第117次迭代：目标函数=208.780665419第118次迭代：目标函数=199.756181363第119次迭代：目标函数=147.2914882第120次迭代：目标函数=209.490628934第121次迭代：目标函数=188.906272556第122次迭代：目标函数=254.106389312第123次迭代：目标函数=182.700159056第124次迭代：目标函数=232.691441966第125次迭代：目标函数=169.767173494第126次迭代：目标函数=170.767353772第127次迭代：目标函数=188.919334801第128次迭代：目标函数=102.401974373第129次迭代：目标函数=186.11050464第130次迭代：目标函数=188.60993634第131次迭代：目标函数=275.973278089第132次迭代：目标函数=362.569042856第133次迭代：目标函数=270.404395791第134次迭代：目标函数=250.264903731第135次迭代：目标函数=240.202115223第136次迭代：目标函数=194.996586644第137次迭代：目标函数=159.9543548第138次迭代：目标函数=113.189269431第139次迭代：目标函数=197.900262404第140次迭代：目标函数=188.860310258第141次迭代：目标函数=209.126474986第142次迭代：目标函数=177.066658993第143次迭代：目标函数=177.066658993第144次迭代：目标函数=177.760102523第145次迭代：目标函数=177.549106799第146次迭代：目标函数=123.595979237第147次迭代：目标函数=178.527727413第148次迭代：目标函数=181.015296862第149次迭代：目标函数=180.151607634第150次迭代：目标函数=179.029042954第151次迭代：目标函数=179.118700201第152次迭代：目标函数=178.492873802第153次迭代：目标函数=267.259531299第154次迭代：目标函数=177.968871038第155次迭代：目标函数=184.081387786第156次迭代：目标函数=255.078061285第157次迭代：目标函数=249.847542458第158次迭代：目标函数=307.417266779第159次迭代：目标函数=250.027959865第160次迭代：目标函数=196.970093641第161次迭代：目标函数=196.970093641第162次迭代：目标函数=209.33031066第163次迭代：目标函数=204.843807568第164次迭代：目标函数=211.044681303第165次迭代：目标函数=149.340167212第166次迭代：目标函数=192.132343418第167次迭代：目标函数=140.343956644第168次迭代：目标函数=181.94000428第169次迭代：目标函数=160.148458356第170次迭代：目标函数=160.242174974第171次迭代：目标函数=110.389926868第172次迭代：目标函数=115.708666994第173次迭代：目标函数=114.032119783第174次迭代：目标函数=114.032119783第175次迭代：目标函数=196.53638158第176次迭代：目标函数=212.763689632第177次迭代：目标函数=191.777288702第178次迭代：目标函数=202.166598824第179次迭代：目标函数=180.713134792第180次迭代：目标函数=280.696308732第181次迭代：目标函数=177.432402563第182次迭代：目标函数=277.773816476第183次迭代：目标函数=284.328371408第184次迭代：目标函数=178.922381458第185次迭代：目标函数=186.85188453第186次迭代：目标函数=186.193903261第187次迭代：目标函数=215.784785022第188次迭代：目标函数=186.998732666第189次迭代：目标函数=166.417149078第190次迭代：目标函数=224.515369868第191次迭代：目标函数=240.503650114第192次迭代：目标函数=224.599534187第193次迭代：目标函数=245.57430458第194次迭代：目标函数=331.210864758第195次迭代：目标函数=236.657697786第196次迭代：目标函数=244.844816837第197次迭代：目标函数=238.634981095第198次迭代：目标函数=267.60921475第199次迭代：目标函数=231.441025132距离和最小为102.4019743附录四 结论图图7 模拟退火迭代训练曲线图图8 模拟退火随机局域搜索算法分配图图9 聚类算法分配图C+代码#include #include #include #include #include #include #include #include #include #define MAX 1e9using namespace std;struct pointdouble x;double y;int p = 4, T = 1000, Iter = 200;double a = 0.999;vector dataset;double getdis(point p1, point p2)return sqrt(pow(p1.x - p2.x), 2) + pow(p1.y - p2.y), 2);double f(vector x)double sum = 0;vectorvector dis;dis.resize(4);for (int i