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:143
if point(i,2) <= 92
p_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:13
f = 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);
end
tmp1 = 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);
end
tmp1 = ind1;tmp2 = ind2;tmp3 = ind3;
% 确保ind1 < ind2 < ind3
if (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;
end
tmplist1 = variables((ind1+1):(ind2-1));
variables((ind1+1):(ind1+ind3-ind2+1)) = variables((ind2):(ind3));
variables((ind1+ind3-ind2+2):ind3) = tmplist1;
endMetropolis准则淘汰劣解
计算出新解的适应度后，与原来的解进行比较，倘若比原来小，直接替换；倘若比原来大，基于概率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);
end
tmp1 = 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);
end
tmp1 = ind1;tmp2 = ind2;tmp3 = ind3;
% 确保ind1 < ind2 < ind3
if (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;
end
tmplist1 = 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;
end
prob=exp(-dE/SAParameters.kb/T);
if(dE>0 && rand()>prob)
var_final = variables_temp; % 还原旧解
end
E0_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方案如下：
第一列为交巡警平台节点标号，第二列为对应前往的全市出入口节点标号。
总结
用智能算法处理整数规划类问题还是比较有意义的。很多连续性问题其实有许多解法，用智能算法往往是大材小用，但是像这种离散类问题，倘若可以用智能算法求解，可以更加灵活，当问题呈几何级数复杂后，也会有时间复杂度和空间复杂度上的优势。
936
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
