2011年数学建模国赛B题的交巡警平台问题是一个经典的整数规划问题，常规的方法是用整数规划求解。但为了练习一下智能算法，我还是决定用模拟退火针对该问题的第二问写一些求解心得。一方面作为学习的记录，另一方面也给有同样问题的朋友们一些启发。

欧式距离计算

首先，初始化数据。各个节点之间有些有路线有些没有路线，根据给出的data计算出有路线的节点之间的距离。该距离将在下文迪杰斯特拉最短路径法的有向图作为权重使用。

%% 数据初始化
pos = readmatrix('data.xls','Sheet','全市交通路口节点数据','Range','B2:C93'); % A区各点坐标
point = readmatrix('data.xls','Sheet','全市交通路口的路线','Range','A2:B144'); % 路径节点
p_start = []; % 开始点
p_end = []; % 结束点
for i = 1:143if point(i,2) <= 92p_start = [p_start, point(i,1)];p_end = [p_end, point(i,2)]; % 筛选A市区范围内的路径end
end
d = sqrt((pos(p_start,1) - pos(p_end,1)).^2 + (pos(p_start,2) - pos(p_end,2)).^2); % 路径距离

迪杰斯特拉最短路径法

根据迪杰斯特拉最短路径法，画出各节点之间路径的有向图，并根据有向图得到最短路径矩阵D。

%% 迪杰斯特拉最短路径法
s = p_start; % 开始点
t = p_end; % 结束点
w = d; % 权重，本题代表路径距离
G = graph(s,t,w); % 生成有向图
plot(G, 'EdgeLabel', G.Edges.Weight, 'linewidth', 1); % 做出有向图
D = distances(G);   % 最短路径矩阵

图其实不重要，对我们真正有用的，是那个最短路径矩阵D。根据D，我们才有方法建立优化模型。

该问题简化后可以如下表述：节点1-20（代表交巡警平台）中派13个到13个指定节点（全市区出入口位置）的最优方案。

主要有两种主要思路，第一种是不同时出发，让总路程最短，该思路省油；第二种思路是同时出发，让总时间最短，该思路省时间。第一种思路建模比较容易，优化函数怎么写比较直观。第二种思路更加贴合实际，但是模型的优化函数建立比较麻烦，且假设警车速度都相同本质上也不是那么合理。因此我采用第一种思路求解。

优化目标设置

%% 围堵方案目标：总路径最短
vital_point = readmatrix('data.xls','Sheet','全市区出入口的位置','Range','C2:C14'); % 交通要道位置
police_point = 1:1:20; police_point = police_point.'; % 交巡警平台位置
% distance_per = D(police_point, vital_point); % 每个交警平台围堵特点节点的距离
% distance_total = sum(distance_per)  % 构造总路径函数obj_fun

其中vital_point就是全市出入口位置。police_point代表交巡警平台。

优化函数

优化函数中，将police_point设置为决策变量x，因为police_point是在20个中选13个，即对象不固定，而vital_point虽然在顺序上不固定，但是对象是固定的。

此时的x是一个序列，每一个值对应的交巡警平台即去向vital_point中相同索引的节点。该序列在1-20中挑选13个数字随机排序。

function f = obj_fun(x,D,vital_point) % 目标函数f = 0;for i = 1:13f = f + D(vital_point(i),x(i));end
end

罚函数

罚函数是处理约束条件的函数，本题中没有约束函数，可以让CV=0。遇到有约束条件的问题，罚函数可以按照如下形式书写。

function CV = obj_fun_CV(x)  % 约束条件函数G1 = (g1(x) > 0).*(g1(x)); % 约束条件1的罚函数G2 = (g2(x) ~= 0).*(g2(x)); % 约束条件2的罚函数CV = G1+G2; % 罚函数累加
end

模拟退火算法过程

第一步，设置参数

设置句柄函数和退火过程的物理量参数。

%% 设置求解问题的参数
problem.numVar = 13;       %变量个数
problem.fun = @(x,D,vital_point)obj_fun(x,D,vital_point); %优化目标函数名称
problem.fun_CV = @(x)obj_fun_CV(x);  %约束条件%% 模拟退火的参数
SAParameters.temperature = 1000; % 初始温度 设置的足够大的话，可以在初始拥有更好的性能
SAParameters.kb = 0.3; % 温度系数
SAParameters.alpha = 0.99; % 降温系数
SAParameters.penalty = 1.5; % 惩罚系数
SAParameters.num = 1000; % 马尔可夫链长度
SAParameters.Tmin = 1; % 结束温度

第二步：初始解的生成

在写优化函数时候已经说过，就是将1-20中挑选13个数随机排列。具体实现时可以先将1-20随机排列，然后挑选前13个，代入优化函数中计算出适应度。

%% 解的初始化
while(1)variables = randperm(20);variables = variables(1:13);variables = variables';  % 初始化一条可行的路径amount = 13;CV = problem.fun_CV(variables);if(CV==0)break;end
end
var_final = variables; % 初始化最终最优解
T = SAParameters.temperature; % 初始化温度
E0_OBJ = problem.fun(variables,D,vital_point); % 初始化目标函数值
E0_CV = problem.fun_CV(variables); % 初始化CV值
E0 = E0_OBJ+SAParameters.penalty*E0_CV; % 最终目标值
E_OBJ_f = E0; % 初始化最佳温度

第三步：退火过程

交换法产生新解

每个温度下，都要进行一个马尔科夫链循环的新解适应度的计算。新解产生的方法采用二交换和三交换法，即随机交换两个或三个序列中的值。

%% 通过二交换或三交换产生一个新解if (rand < 0.5) % 随机决定是进行两交换还是三交换% 两交换法产生新的解ind1 = 0; ind2 = 0;while (ind1 == ind2)ind1 = ceil(rand.*amount);ind2 = ceil(rand.*amount);endtmp1 = variables(ind1);variables(ind1) = variables(ind2);variables(ind2) = tmp1; % 随机交换2个变量的位置else% 三交换法产生新的解ind1 = 0; ind2 = 0; ind3 = 0;while (ind1 == ind2) || (ind1 == ind3) ...|| (ind2 == ind3) || (abs(ind1-ind2) == 1)ind1 = ceil(rand.*amount);ind2 = ceil(rand.*amount);ind3 = ceil(rand.*amount);endtmp1 = ind1;tmp2 = ind2;tmp3 = ind3;% 确保ind1 < ind2 < ind3if (ind1 < ind2) && (ind2 < ind3)elseif (ind1 < ind3) && (ind3 < ind2)ind2 = tmp3;ind3 = tmp2;elseif (ind2 < ind1) && (ind1 < ind3)ind1 = tmp2;ind2 = tmp1;elseif (ind2 < ind3) && (ind3 < ind1) ind1 = tmp2;ind2 = tmp3; ind3 = tmp1;elseif (ind3 < ind1) && (ind1 < ind2)ind1 = tmp3;ind2 = tmp1; ind3 = tmp2;elseif (ind3 < ind2) && (ind2 < ind1)ind1 = tmp3;ind2 = tmp2; ind3 = tmp1;endtmplist1 = variables((ind1+1):(ind2-1));variables((ind1+1):(ind1+ind3-ind2+1)) = variables((ind2):(ind3));variables((ind1+ind3-ind2+2):ind3) = tmplist1;end

Metropolis准则淘汰劣解

计算出新解的适应度后，与原来的解进行比较，倘若比原来小，直接替换；倘若比原来大，基于概率P接受新解。该解的迭代法称为Metropolis准则。

%% 移动后的目标值计算
E_OBJ = problem.fun(variables,D,vital_point); % 移动后的目标函数值
E_CV = problem.fun_CV(variables); % 移动后的CV值
E = E_OBJ+SAParameters.penalty*E_CV;
dE = E-E0;
if (dE<0 && E_CV==0)var_final = variables; % 优解迭代E_OBJ_f=E_OBJ;
end
prob=exp(-dE/SAParameters.kb/T); % P值计算
if(dE>0 && rand()>prob)var_final = variables_temp; % 按P概率接受劣解,否则还原
end

退火过程总代码如下：

while (T>=SAParameters.Tmin) % 开始降温for i = 1:SAParameters.num % 马尔科夫链temp_variables = variables;%% 通过二交换或三交换产生一个新解if (rand < 0.5)   % 随机决定是进行两交换还是三交换% 两交换法产生新的解ind1 = 0; ind2 = 0;while (ind1 == ind2)ind1 = ceil(rand.*amount);ind2 = ceil(rand.*amount);endtmp1 = variables(ind1);variables(ind1) = variables(ind2);variables(ind2) = tmp1; % 随机交换2个变量的位置else% 三交换法产生新的解ind1 = 0; ind2 = 0; ind3 = 0;while (ind1 == ind2) || (ind1 == ind3) ...|| (ind2 == ind3) || (abs(ind1-ind2) == 1)ind1 = ceil(rand.*amount);ind2 = ceil(rand.*amount);ind3 = ceil(rand.*amount);endtmp1 = ind1;tmp2 = ind2;tmp3 = ind3;% 确保ind1 < ind2 < ind3if (ind1 < ind2) && (ind2 < ind3)elseif (ind1 < ind3) && (ind3 < ind2)ind2 = tmp3;ind3 = tmp2;elseif (ind2 < ind1) && (ind1 < ind3)ind1 = tmp2;ind2 = tmp1;elseif (ind2 < ind3) && (ind3 < ind1) ind1 = tmp2;ind2 = tmp3; ind3 = tmp1;elseif (ind3 < ind1) && (ind1 < ind2)ind1 = tmp3;ind2 = tmp1; ind3 = tmp2;elseif (ind3 < ind2) && (ind2 < ind1)ind1 = tmp3;ind2 = tmp2; ind3 = tmp1;endtmplist1 = variables((ind1+1):(ind2-1));variables((ind1+1):(ind1+ind3-ind2+1)) = variables((ind2):(ind3));variables((ind1+ind3-ind2+2):ind3) = tmplist1;end%% 移动后的目标值计算E_OBJ = problem.fun(variables,D,vital_point); % 移动后的目标函数值E_CV = problem.fun_CV(variables); % 移动后的CV值E = E_OBJ+SAParameters.penalty*E_CV;dE = E-E0;if (E_OBJ<=E_OBJ_f && E_CV==0)var_final = variables; % 保留新解E_OBJ_f=E_OBJ;endprob=exp(-dE/SAParameters.kb/T);if(dE>0 && rand()>prob)var_final = variables_temp; % 还原旧解endE0_OBJ=problem.fun(variables,D,vital_point); %初始目标函数值E0_CV=problem.fun_CV(variables); %初始CV值E0=E0_OBJ+SAParameters.penalty*E0_CV;end
T = T*SAParameters.alpha; % 降温
end

计算结果

我最后计算出来的最短总距离为1038.75，vital_point对于police_point方案如下：

第一列为交巡警平台节点标号，第二列为对应前往的全市出入口节点标号。

总结

用智能算法处理整数规划类问题还是比较有意义的。很多连续性问题其实有许多解法，用智能算法往往是大材小用，但是像这种离散类问题，倘若可以用智能算法求解，可以更加灵活，当问题呈几何级数复杂后，也会有时间复杂度和空间复杂度上的优势。

基于模拟退火的交巡警平台调度问题相关推荐

matlab交巡警平台,数学建模2011年交巡警平台题目matlab进行数据练习
题目: 题目给出的数据: 还有城市中平台和道路设置的图示. 本次练习的目的: 主要熟悉matlab导入excel数据,编写matlab程序对数据进行处理练习需要完成的目标: 通过编写matlab程序 ...
交巡警求全市距离MATLAB程序,数学建模国赛优秀论文-交巡警服务平台的设置与调度.doc...
PAGE 第 PAGE 2页,共 NUMPAGES 1页 2011高教社杯全国大学生数学建模竞赛承诺书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则. 我们完全明白,在竞赛开始后参赛队员不能以 ...
交巡警求全市距离MATLAB程序,交巡警服务平台的设置与调度论文
交巡警服务平台的设置与调度摘要合理安排交巡警管辖范围具有现实的意义,既提高巡警的办事效率又可以更方便的为百姓解决问题. 针对问题一中的第一小问,建立路口搜索模型.编写MATLAB程序,对A区每个平 ...
[数学建模] 2. 交巡警服务平台的设置与调度
必看: 数据与司守奎老师的lingo书籍---<LINGO软件及应用>PDF版.数据及代码,百度云链接,自取,仅供共享学习.评论区的链接没有书籍的PDF,论文找不到了...,不然也给大家直 ...
全国大学生数学建模竞赛2011B题交巡警服务平台的设置与调度MATLAB程序与LINGO程序
目录一.第1问二.第2问三.第3问四.第4问五.第5问一.第1问 %第1问MATLAB程序 %% 对应计算两点间距离(A区) a=xlsread('2011B.xls','全市交通路口节点 ...
交巡警求全市距离MATLAB程序,高教社杯全国大学生数学建模竞赛B题参考答案
交巡警服务平台的设置与调度优化分析摘要本文以实现警察的刑事执法.治安管理.交通管理.服务群众四大职能为宗旨,利用有限的警务资源,根据城市的实际情况与需求合理地设置了交巡警服务平台.分配各平台的管 ...
2011年数学建模国赛B题(交巡警服务平台的设置与调度模型)论文.doc
本博客下载链接有可以修改的word版本, 可免费下载阅览学习, 也可作为数学建模相关课程作业修改上交. 链接:https://pan.baidu.com/s/1xJKV-MeA1qJAevXpfLpk ...
交巡警求全市距离MATLAB程序,数学建模论文交巡警服务平台的设置与调度
内容简介: 数学建模论文交巡警服务平台的设置与调度本文针对如何设置交巡警服务平台.各平台的管辖范围.警务资源调度等相关问题建立模型并进行研究.利用MATLAB软件对数据进行处理,计算出任意两节点的 ...
交巡警求全市距离MATLAB程序,交巡警服务平台的设置与调度数学建模论文_
第二类:? ?????>∈=∈30min 2ij p j c V i A ,即节点i 到任何交巡警服务平台节点的最短距离大于mm 30,这样的节点有92,61,39,38,29,28,将其转化为 ...

基于模拟退火的交巡警平台调度问题