【建模算法】matlab基于求解器intlinprog求解TSP问题

最新推荐文章于 2024-04-25 16:09:41 发布
最新推荐文章于 2024-04-25 16:09:41 发布
阅读量2.4k 收藏 26
点赞数 4
分类专栏： 数学建模 文章标签： matlab 算法 图论
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/baidu/article/details/124990349
版权

【建模算法】matlab基于求解器intlinprog求解TSP问题

本案例说明如何使用二元整数规划来求解经典的TSP问题。此问题涉及找到一条历经一系列停留点（城市）的最短回路（路径）。在本例中有 52 个停留点，但你可以很轻松地更改 nStops 变量以得到不同规模的问题。对最初的问题进行求解后得到的解会包含子回路。这意味着找到的最优解并没有给出一条穿过所有点的连续路径，而是有几个独立的环路。然后，你可以使用迭代过程来确定子回路，添加约束，并重新运行优化，直到消除子回路。

一、问题描述

​ 本案例以52个城市为例，假定52个城市的位置坐标如表1所列。寻找出一条最短的遍历52个城市的路径。

城市编号X坐标Y坐标
1565575
225185
3345750
4945685
5845655
6880660
725230
85251000
95801175
106501130
111605620
121220580
131465200
1415305
15845680
16725370
17145665
18415635
19510875
20560365
21300465
22520585
23480415
24835625
25975580
261215245
271320315
281250400
29660180
30410250
31420555
32575665
3311501160
34700580
35685595
36685610
37770610
38795645
39720635
40760650
41475960
4295260
43875920
44700500
45555815
46830485
47117065
48830610
49605625
50595360
511340725
521740245

二、求解步骤

2.1 问题表示

将要进行整数线性规划的TSP问题表示如下：

  • 生成所有可能的行程，意味着经过所有不同停留点对组。
  • 计算每次行程的距离。
  • 要最小化的代价函数是行程中每次旅行的旅行距离之和。
  • 决策变量是与每个行程相关联的二元变量，其中每个 1 表示回路中存在的一次行程，每个 0 表示不在回路中的一次行程。
  • 为确保回路包括每个经停留点，应加入这样一个线性约束：每个停留点都正好涉及两段行程。这意味着一段行程到达该停留点，一段行程离开该停留点。

2.2 生成停留点

读取数据，生成停留点

data=importdata('city052.xlsx');
X=data.data.Sheet1(:,2:3);   %各城市坐标

nStops = 52; %52城市，可修改, 问题规模为N^2
stopsLon = X(:, 1);   %城市x坐标
stopsLat = X(:, 2);   %城市y坐标

2.3 计算停留点之间的距离

由于存在52个停留点，因此有 1326 次行程，这意味着有 1326 个二元变量（# 变量 = 52 选择 2）。

生成所有行程，这意味着生成所有停留点对组。

idxs = nchoosek(1:nStops,2);

计算所有行程距离，假设地球表面是平的以便使用毕达哥拉斯法则。

dist = hypot(stopsLat(idxs(:,1)) - stopsLat(idxs(:,2)), ...
             stopsLon(idxs(:,1)) - stopsLon(idxs(:,2)));
lendist = length(dist);

根据 dist 向量的定义，一次回路的长度为

dist'*x_tsp

其中 x_tsp 是二元解向量。这是您尝试最小化的旅行距离。

2.4 创建图和绘制地图

将问题表示为图。创建一个图，其中停留点是节点，行程是边。

G = graph(idxs(:,1),idxs(:,2));

使用图论图显示停留点。绘制没有图边的节点。

%绘制停留点
figure;
hGraph = plot(G,'XData',stopsLon,'YData',stopsLat,'LineStyle','none','NodeLabel',{});
for i = 1: size(X, 1)
    text(X(i, 1), X(i, 2), [' ',num2str(i)], 'color', [1, 0, 0]);
end
hold on

在这里插入图片描述

2.5 约束

创建线性约束，使每个停留点都有两个关联行程，因为每个停留点上必须有一次到达该停留点的行程和一次离开该停留点的行程。

%创建线性约束，使每个停留点都有两个关联行程，因为每个停留点上必须有一次到达该停留点的行程和一次离开该停留点的行程
Aeq = spalloc(nStops,length(idxs),nStops*(nStops-1)); %稀疏矩阵
for ii = 1:nStops
    whichIdxs = (idxs == ii); %查找包含站点ii的行程
    whichIdxs = sparse(sum(whichIdxs,2)); %包括ii位于两端的行程
    Aeq(ii,:) = whichIdxs'; %包含在约束矩阵中
end
beq = 2*ones(nStops,1);

2.6 二元边界

所有决策变量均为二元变量。现在，将 intcon 参数设置为决策变量的数目，设置每个变量的下界为 0，上界为 1。

intcon = 1:lendist;
lb = zeros(lendist,1);
ub = ones(lendist,1);

2.7 使用 intlinprog 进行优化

此问题已可以求解。要隐藏迭代输出，请关闭默认的迭代输出。

opts = optimoptions('intlinprog','Display','off');
[x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts);

创建一个新图，将解行程作为边。为此，在某些值不是精确整数的情况下对解进行舍入，并将生成的值转换为 logical。

x_tsp = logical(round(x_tsp));
Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2),[],numnodes(G));

2.8 可视化解

%可视化解
highlight(hGraph,Gsol,'LineStyle','-')
title('具有子回路的解')

从图上可以看出，解有几个子回路。到当前为止指定的约束没有阻止这些子回路的发生。要防止出现任何可能的子回路，您将需要极大数量的不等式约束。

2.9 子回路约束

由于无法添加所有子回路约束，需要采用迭代方法。检测当前解中的子回路，然后添加不等式约束以防止发生这些特定子回路。通过执行此操作，您只需几次迭代即可找到合适的回路。

使用不等式约束消除子回路。具体的工作原理如下：如果您在一个子回路中有五个点，则您有五条连接这些点的线来创建该子回路。可实现一个不等式约束来消除此子回路，即，这五个停留点之间存在的线必须小于或等于四条。

还要找到这五个停留点之间的所有线，并限制解中这些线不超过四条。这是正确的约束，因为如果解中存在五条或更多条线，则该解将有子回路（具有 n 个节点和 n 条边的图始终包含一个环路）。

通过识别 Gsol 中的连通分量来检测子回路，该图是用当前解中的边构建的。conncomp 返回一个向量，该向量包含每条边所属的子回路的编号。

tourIdxs = conncomp(Gsol);
numtours = max(tourIdxs); %子回路数目
fprintf('当前子回路数: %d\n',numtours);

加入线性不等式约束以消除子回路，并反复调用求解器求解，直到只剩下一个子回路。

A = spalloc(0,lendist,0); %分配稀疏线性不等式约束矩阵
b = [];

pic_num = 1;
while numtours > 1 %当子回路大于1时进入循环 
    %增加子回路约束
    b = [b;zeros(numtours,1)]; %分配b
    A = [A;spalloc(numtours,lendist,nStops)]; %猜测要分配多少个非零
    for ii = 1:numtours
        rowIdx = size(A,1) + 1; %索引计数器 
        subTourIdx = find(tourIdxs == ii); %提取当前子回路
        variations = nchoosek(1:length(subTourIdx),2);
        for jj = 1:length(variations)
            whichVar = (sum(idxs==subTourIdx(variations(jj,1)),2)) & ...
                       (sum(idxs==subTourIdx(variations(jj,2)),2));
            A(rowIdx,whichVar) = 1;
        end
        b(rowIdx) = length(subTourIdx) - 1;
    end

    %再次尝试优化
    [x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,A,b,Aeq,beq,lb,ub,opts);
    x_tsp = logical(round(x_tsp));
    Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2),[],numnodes(G));
    
    %可视化结果
    hGraph.LineStyle = 'none'; %删除之前高亮显示的路径
    highlight(hGraph,Gsol,'LineStyle','-')
    drawnow;
    pause(1);
    F=getframe(gcf);
    I=frame2im(F);
    [I,map]=rgb2ind(I,256);
    if pic_num == 1
        imwrite(I,map,'TSP052.gif','gif', 'Loopcount',inf,'DelayTime',0.5);
    else
        imwrite(I,map,'TSP052.gif','gif','WriteMode','append','DelayTime',0.5);
    end
    
    %本次多少子回路?
    tourIdxs = conncomp(Gsol);
    numtours = max(tourIdxs); %子回路数目
    fprintf('当前子回路数: %d\n',numtours)
    pic_num = pic_num + 1;
end

title('最优解');
hold off

2.10 解质量

此解表示一个可行的回路，因为它是单一环路。但它是否为一次代价最低的回路？找出答案的一种方法是检查输出结构体。

disp(output.absolutegap)

绝对间隙小意味着该解或者是最优的，或者总长度接近最优。

三、求解结果

示例动画（此处为抓取的32城市求解动画）。
在这里插入图片描述
当前子回路数: 7
当前子回路数: 1

总距离：7544.3659

最优解轨迹图：
在这里插入图片描述

matlab完整源代码

%52城市TSP问题求解
clear,clc;
data=importdata('city052.xlsx');
X=data.data.Sheet1(:,2:3);   %各城市坐标

nStops = 52; %52城市，可修改, 问题规模为N^2
stopsLon = X(:, 1);   %城市x坐标
stopsLat = X(:, 2);   %城市y坐标


%生成所有停留点对组
idxs = nchoosek(1:nStops,2);

%计算所有行程距离，假设地球表面是平的以便使用毕达哥拉斯法则
dist = hypot(stopsLat(idxs(:,1)) - stopsLat(idxs(:,2)), ...
             stopsLon(idxs(:,1)) - stopsLon(idxs(:,2)));
lendist = length(dist);

%创建一个图，其中停留点是节点，行程是边。
G = graph(idxs(:,1),idxs(:,2));

%绘制停留点
figure;
hGraph = plot(G,'XData',stopsLon,'YData',stopsLat,'LineStyle','none','NodeLabel',{});
for i = 1: size(X, 1)
    text(X(i, 1), X(i, 2), [' ',num2str(i)], 'color', [1, 0, 0]);
end
hold on

%创建线性约束，使每个停留点都有两个关联行程，因为每个停留点上必须有一次到达该停留点的行程和一次离开该停留点的行程
Aeq = spalloc(nStops,length(idxs),nStops*(nStops-1)); %稀疏矩阵
for ii = 1:nStops
    whichIdxs = (idxs == ii); %查找包含站点ii的行程
    whichIdxs = sparse(sum(whichIdxs,2)); %包括ii位于两端的行程
    Aeq(ii,:) = whichIdxs'; %包含在约束矩阵中
end
beq = 2*ones(nStops,1);

%所有决策变量均为二元变量，将intcon参数设置为决策变量的数目，设置每个变量的下界为 0，上界为 1。
intcon = 1:lendist;
lb = zeros(lendist,1);
ub = ones(lendist,1);

%使用 intlinprog 进行优化求解
opts = optimoptions('intlinprog','Display','off');
[x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts);

%创建一个新图，将解行程作为边，在某些值不是精确整数的情况下对解进行舍入，并将生成的值转换为 `logical`。
x_tsp = logical(round(x_tsp));
Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2),[],numnodes(G));

%可视化解
highlight(hGraph,Gsol,'LineStyle','-')
title('具有子回路的解')

tourIdxs = conncomp(Gsol);
numtours = max(tourIdxs); %子回路数目
fprintf('当前子回路数: %d\n',numtours);

A = spalloc(0,lendist,0); %分配稀疏线性不等式约束矩阵
b = [];

pic_num = 1;
while numtours > 1 %当子回路大于1时进入循环 
    %增加子回路约束
    b = [b;zeros(numtours,1)]; %分配b
    A = [A;spalloc(numtours,lendist,nStops)]; %猜测要分配多少个非零
    for ii = 1:numtours
        rowIdx = size(A,1) + 1; %索引计数器 
        subTourIdx = find(tourIdxs == ii); %提取当前子回路
        variations = nchoosek(1:length(subTourIdx),2);
        for jj = 1:length(variations)
            whichVar = (sum(idxs==subTourIdx(variations(jj,1)),2)) & ...
                       (sum(idxs==subTourIdx(variations(jj,2)),2));
            A(rowIdx,whichVar) = 1;
        end
        b(rowIdx) = length(subTourIdx) - 1;
    end

    %再次尝试优化
    [x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,A,b,Aeq,beq,lb,ub,opts);
    x_tsp = logical(round(x_tsp));
    Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2),[],numnodes(G));
    
    %可视化结果
    hGraph.LineStyle = 'none'; %删除之前高亮显示的路径
    highlight(hGraph,Gsol,'LineStyle','-')
    drawnow;
    pause(1);
    F=getframe(gcf);
    I=frame2im(F);
    [I,map]=rgb2ind(I,256);
    if pic_num == 1
        imwrite(I,map,'TSP052.gif','gif', 'Loopcount',inf,'DelayTime',0.5);
    else
        imwrite(I,map,'TSP052.gif','gif','WriteMode','append','DelayTime',0.5);
    end
    
    %本次多少子回路?
    tourIdxs = conncomp(Gsol);
    numtours = max(tourIdxs); %子回路数目
    fprintf('当前子回路数: %d\n',numtours)
    pic_num = pic_num + 1;
end

title('最优解');
hold off

%disp(output.absolutegap);  %输出解的质量
disp(['总距离：', num2str(costopt)]);
disp('-------------------------');

本案例数据与源代码下载地址：https://download.csdn.net/download/baidu/85467610

果州做题家
关注
  • 4
    点赞
  • 踩
  • 26
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
【MATLAB模型】线性规划linprog 0-1规划intlinprog，非线性规划fmincon 罚函数法fminunc
什么是快乐星球的博客
07-30 2781
线性、非线性的区别： matlab中线性规划的标准形式： 线性规划的例子： 求最大值，则 -c matlab中线性规划的代码： %max 线性规划 clc;clear; c=[2;3;-5]; a=[-2,5,-1;1,3,1]; b=[-10;12]; aeq=[1,1,1]; beq=7; [x,y]=linprog(-c,a,b,aeq,beq,zeros(3,1)) %min 线性规划 clc;clear all; c=[2;3;1]; a=[-...
基于遗传算法的TSP问题求解（matlab源码）
04-26
TSP (traveling salesman problem，旅行商问题)是典型的NP完全问题，即其最坏情况下的时间复杂度随着问题规模的增大按指数方式增长，到目前为止还未找到一个多项式时间的有效算法。本资源给出了“基于遗传算法的一个TSP问题求解”的matlab完整源代码。
matlab基于求解器intlinprog求解TSP问题
05-26
matlab基于求解器intlinprog求解52城市TSP问题完整数据与代码。本案例说明如何使用二元整数规划来求解经典的TSP问题。此问题涉及找到一条历经一系列停留点（城市）的最短回路（路径）。在本例中有 52 个停留点，但你可以很轻松地更改 `nStops` 变量以得到不同规模的问题。对最初的问题进行求解后得到的解会包含子回路。这意味着找到的最优解并没有给出一条穿过所有点的连续路径，而是有几个独立的环路。然后，你可以使用迭代过程来确定子回路，添加约束，并重新运行优化，直到消除子回路。
基于MATLAB和YALMIP的CPLEX求解器的配电网故障重构模型研究——以33节点为例
最新发布
lxQGfFrIn的博客
04-25 399
配电网故障重构是一项关键的电力系统运维任务，旨在快速准确地定位和修复配电网中的故障节点，保障电网的连通性和供电稳定性。基本内容：以33节点为研究对象，编制配电网故障重构模型，采用图论知识保证配电网的连通性和辐射性，以网损和负荷损失作为目标函数，包括潮流约束、电压电流约束、sop约束、辐射性约束等，程序运行稳定。首先，我们建立了一套完整的重构模型，具体步骤如下。在故障重构过程中，我们需要考虑潮流约束、电压电流约束、SOP约束以及辐射性约束等多个约束条件，以确保重构后的电网能够满足安全可靠的运行要求。
MATLAB线形规划函数linprog、intlinprog与二次规划函数quadprog
STM89C56的博客
04-26 6336
一、目标函数简介 线形规划问题的数学模型: s.t 线形规划，顾名思义就是目标函数与约束条件均为线形函数（一次函数）。 二、线形规划函数linprog linprog函数的用法大致分为以下几种用法： 1、不等式约束 x = linprog(f,A,b) ，f为目标函数，A为不等式约束的系数矩阵，b为不等式右侧的值，返回值x为...
Matlab绘制双纽线、莫比乌斯环，双纽线公式、重积分
Wayne的CSDN博客
08-13 8572
Matlab绘制双纽线、莫比乌斯环；双纽线的二重积分
27 用linprog、fmincon求 解线性规划问题（matlab程序）
m0_57943157的博客
07-30 1554
与linprog相比，多了参数intcon,代表了整数决策变量所在的位置。即使赛题中有非线性目标/约束，第一想法也应是将其。A=[6,-1,1;在比赛时，遇到非线性形式是家常便饭。不同资产配置下收益率/风险不同，总资金有限；求解此类问题的一般步骤，后几期会逐步增加用。求解线性规划问题，求目标函数的最小值，f=[-5,4,2];，更多细节参见帮助文档。l=[-1,0,0];求最大值时，c加上负号：-c。%%%%用linprog求解。%%%%用fmincon求解。
Matlab中intlinprog函数的用法总结
weixin_44092113的博客
12-01 1437
intlinprog是matlab中用于求解混合整数线性规划(Mixed-integer linear programming)的一个函数，用法基本和linprog差不多Matlab中，该模型的标注写法如下\[f,x,b,beq,lb,ub,intcon是向量；A和Aeq是矩阵 \]
遗传算法_matlab_TSP问题求解
04-10
资源名：遗传算法_matlab_TSP问题求解 资源类型：matlab项目全套源码 源码说明： 全部项目源码都是经过测试校正后百分百成功运行的，如果您下载后不能运行可联系我进行指导或者更换。 适合人群：新手及有一定经验的...
数学建模混合粒子群算法求解TSP问题代码
07-26
数学建模：混合粒子群算法求解TSP问题的matlab代码，可学习参考。
matlab遗传算法应用TSP问题求解含源代码.zip
11-16
matlab遗传算法应用TSP问题求解含源代码.zip
数学建模-模拟退火算法求解TSP问题Matlab源码.zip
01-20
数学建模-模拟退火算法求解TSP问题Matlab源码.zip
matlab loopcount,求助一个数值积分问题，用matlab的quadgk函数来计算，谢谢！
weixin_36266554的博客
03-22 615
引用回帖:wurongjun at 2017-07-31 11:49:49你好!我的Matlab是旧版!麻烦你把你的quadgk函数代码贴一个上来,谢谢啦!用命令 type quadgk就可以看到源代码,在复制粘贴即可!function [q,errbnd] = quadgk(FUN,a,b,varargin)%QUADGKNumerically evaluate integral, adap...
模拟退火算法解决TSP（python实现 110+行代码）【gif生成】
肥宅Sean
12-19 3951
简述 代码我是基于我之前写的两篇，一篇是遗传算法TSP的Python实现，一篇是模拟退火算法的解决TSP的C++实现。 模拟退火算法理论+Python解决函数极值+C++实现解决TSP问题 遗传算法解决TSP问题 Python实现【160行以内代码】 效果演示 对比 相比于遗传算法来说没有保持历史中的较优数据，但是通过退火的算法思维，很邻近点搜索的想法，任然能保持较为正确的收敛结果 代码 i...
Matlab中的优化问题求解方法与示例分析
vipfanxu的博客
12-12 572
除了线性规划问题外，Matlab还提供了非线性规划问题和整数规划问题的求解方法，读者可以根据具体问题的需求选择适合的函数进行求解。lb和ub为变量的下界和上界。其中，f为目标函数的系数向量，A和b为不等式约束的系数矩阵和常数向量，Aeq和beq为等式约束的系数矩阵和常数向量，lb和ub为变量的下界和上界。其中，f为目标函数的系数向量，intcon为整数变量的索引向量，其他参数的含义与之前介绍的函数相对应。根据上述代码，可以得到最大化利润为22500元，此时产品A的生产数量为4个，产品B的生产数量为2个。
遗传算法解决TSP问题MATLAB实现（详细）
热门推荐
徐奕的专栏
02-01 8万+
问题定义：巡回旅行商问题 给定一组n个城市和俩俩之间的直达距离，寻找一条闭合的旅程，使得每个城市刚好经过一次且总的旅行距离最短。 TSP问题也称为货郎担问题，是一个古老的问题。最早可以追溯到1759年Euler提出的骑士旅行的问题。1948年，由美国兰德公司推动，TSP成为近代组合优化领域的典型难题。 TSP是一个具有广泛的应用背景和重要理论价值的组合优化问题。 近年来，有很多解决该问题的较为有效...
Matlab中 intlinprog函数用法简介
LiDong的博客
11-07 6802
Matlab中 intlinprog函数用法简介 本来想要自己亲手写一遍的，发现了一优质博文基本上跟我做过的例题大差不差，所以就直接放上链接。 参考链接 https://www.cnblogs.com/goodtwo/p/11145699.html ...
运筹学1——整数规划之分支定界法与MATLAB的intlinprog函数
奇迹的博客
06-25 3894
1 背景 直接用fmincon求解最优化问题时，得到的结果可能是小数，但针对某一些问题，要求结果必须是整数，称为整数规划问题。 分支定界法是一种求解整数规划的算法，具体介绍可看 分支定界法MATLAB实现 具体问题： 2 用MATLAB直接求解整数规划 intlinprog是一个求解整数规划的函数，具体使用方法和fmincon类似，不清楚fmincon的可看我之前的文章MATLAB优化函数fmincon。简单来说，就是把优化问题化为标准形式。 令x=[x11,x12,x13,…,x54,x55,x56]
TSP算法（0-1规划） ---用Matlab和Lingo联合实现
Hello_leiyuanyi的博客
08-20 3014
  MATLAB处理数据的程序：   处理之前需要进行数据的处理，把数据转换为0和1，转化为相关的矩阵；用i表示行，用j表示列，增加一行用0表示。   文件名根据需要修改即可。 load 'C:\Users\Administrator\Desktop\table2.txt'; a = table2; [m,n]=size(a); d=zeros(n+1,n+1); %定义距离矩阵;...
基于matlab灰狼算法求解多旅行商问题
02-08
综上所述，基于matlab的灰狼算法可以用来求解多旅行商问题。通过灰狼算法的迭代优化策略，可以得到近似最优解。这种方法具有较高的搜索能力和全局优化能力，在实际应用中具有一定的实用性和效果。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值