💺


敲到码穷处,坐看云起时。☁️ ☁️ ☁️

数学建模系列文章——总结篇《数模美一国一退役选手的经验分享[2021纪念版]》.


🍬

一、前驱知识准备

1、Lingo简介

  LINGO是 Linear Interactive and General Optimizer 的缩写,即“交互式的线性和通用优化求解器”。

  Lingo超强的优化计算能力在很多方面(线性规划、非线性规划、线性整数规划、非线性整数规划、非线性混合规划、二次规划等)比matlab、maple等强得多。

  Lingo编程简洁明了,数学模型不用做太大的改动,便可以直接用Lingo语言编程呈现,十分直观


2、01规划、整数规划、(非)线性规划等规划的区别

☁️ ☁️ ☁️

  ①0-1规划 是决策变量仅取值0或1的一类 特殊的整数规划

   补充:0-1 变量可以数量化地描述诸如开与关、取与弃、有与无等现象所反映的逻辑关系。因此 0-1 规划非常适合解决如工厂选址 、生产计划安排、背包问题、人员安排等多种问题。

  ②整数规划 是指规划中的变量(全部或部分)限制为整数。

  补充:若在线性模型中,变量限制为整数,则称为整数线性规划。所流行的求解整数规划的方法往往只适用于整数线性规划。从约束条件的构成又可细分为线性,二次和非线性的整数规划。

  ③线性规划(Linear programming,简称 LP ),是运筹学中的一个重要分支,它是一种数学方法。研究 线性约束条件线性目标函数极值问题 的数学理论和方法。

  ④非线性规划 是一种求解目标函数或约束条件中有一个或几个 非线性函数 的最优化问题的方法。


3、优化(/规划)模型的组成

(1)目标函数

  (1)目标函数:一般表示成求 某个数学表达式 的 最小值最大值

(2)决策变量

  (2)决策变量:目标函数值 取决于 哪些变量

(3)约束条件

  (3)约束条件:对变量附加一些条件限制(通常用等式或不等式来表示)。


4、Lingo模型的基本组成

(1)集合定义部分

(2)数据段

(3)初始化段

(4)计算段

(5)目标函数和约束条件

  :详细说明可参考本文最后的 参考附录【2】 📚。


二、样例:国赛数模11年B题《交巡警服务平台的设置与调度》的第一问(前2小问)

1、问题描述:

  该市中心城区A的交通网络和现有的20个交巡警服务平台的设置情况如下图所示。(相关数据可在本文最后 参考附录【1】下载)

  ①现在请为各交巡警服务平台分配管辖区域,使其在所管辖的范围出现突发事件时,尽量能在3分钟内有交巡警(警车时速为60km/h)到达事发地。

  ②对于重大突发事件,需要调度全区20个交巡警服务平台的警力资源,对进出该去的13条交通要道实现快速全封锁。实际中一个平台的警力最多封锁一个结点,请给出该区交巡警服务平台合理的调度方案。【注:13条交通要道即图中13个红点。】

在这里插入图片描述




2、解题思路


  ①根据实际情况,我们首先要对数据进行处理,利用Folyd算法,借助 matlab 软件对其进行求解,求出20个交巡警服务平台到各个结点之间的实际距离
  然后建立优化模型,找出目标函数,并找出相应的约束条件,借助 lingo 对其进行求解,得出结论。


  ②由上面那张图,我们明显地可以看出,有些结点到所有的交巡警平台的距离都不会小于3km,所以显然不能满足 “3分钟” 那个条件(因为车速60km/h→1km/min)。故我们考虑到交巡警平台的重新铺设,将20个交巡警平台重新在交通网络中铺设,得到最优结果。


  ③我们设定一个正向型指标,使增设平台的目标就是:让该 指标 分数最大。



3、数据预处理


Folyd算法如下:

  Step1:构造 0-1 矩阵 A0
在这里插入图片描述
  来存放各节点之间的关系,其中,aij = 0 表示 i,j 两节点之间 没有 直接连接。 a ij = 1 表示 i,j 两节点之间直接连接。

  Step2:递推产生一个矩阵序列 A0、A1、A2、… 、Ak、Ak+1、… 、An-1、An。其中 Ak ( i , j ) 表示从节点 i 到结点 j 的路径所经过的节点序号不大于 k 的最短路径长度。
  计算时用的迭代公式为:

A k ( i , j ) = m i n ( A k − 1 ( i , j ) , A k − 1 ( i , k ) + A k − 1 ( k , j ) ) A_k(i,j)=min(A_{k-1}(i,j),A_{k-1}(i,k)+A_{k-1}(k,j)) Ak(i,j)=min(Ak1(i,j),Ak1(i,k)+Ak1(k,j))

  其中,k 为迭代次数, i, j, k = 1, 2, 3, … , n
  最后,当 k=n 时,An 即是各节点之间的最短距离。

  注: Table是个 n*m 的结构体表 ,封装了题目里面的所有数据。
  其构造的内容详见 这一篇文章《数学建模——matlab绘制 地图 散点图连线图 (运用plot、scatter、struct、xlsread等函数)》
  链接: https://blog.csdn.net/Wang_Dou_Dou_/article/details/119007126?spm=1001.2014.3001.5501.
  Folyd 的原理可以参考这篇《算法篇——详解 弗洛伊德(Folyd)算法【用 C/C++ 和 matlab 实现】》
  链接:https://blog.csdn.net/Wang_Dou_Dou_/article/details/119485893?spm=1001.2014.3001.5501.

  通过 matlab 编程求解如下:


A_0 = zeros(92,92);   % 初始化 注:A区域有92个节点
for i=1:92
    next_point = Table(i).NextJieDian;   % 注:Table(i).NextJieDian 是一个 1*n 的矩阵
    for j=1:length(next_point)
        if 1 <= next_point(j) && next_point(j) <= 92
            A_0(i,next_point(j)) = 1;
        end
    end
end

zuo_x_A = zeros(92,1);
zuo_y_A = zeros(92,1);
for i=1:92    % 筛选出 A 区域里面的数据点
    zuo_x_A(i) = Table(i).X;
    zuo_y_A(i) = Table(i).Y;
end

X_0 = zeros(92,92);
for i=1:92    % 计算A区域中 两两结点之间的距离
    for j=1:92
        X_0(i,j) = sqrt( ( zuo_x_A(i) - zuo_x_A(j) ).^2 + ( zuo_y_A(i) - zuo_y_A(j) ).^2 );   
    end
end
Fal_X = X_0 .* A_0;   % 只保留 两个结点有连线的 "距离值"

for i=1:92
    for j=1:92
        if A_0(i,j) == 0    % 距离等于0,表示网络中两节点之间不相连,则表示距离无穷大
            Fal_X(i,j) = 99999;   
        end
    end
end

for i=1:92
    Fal_X(i,i) = 0;   % 自身到自身的节点距离为0
end

% Folyd算法:
path = zeros(92);
for k=1:92
    for i=1:92
        for j=1:92
            if Fal_X(i,j) > Fal_X(i,k) + Fal_X(k,j)
                Fal_X(i,j) = Fal_X(i,k) + Fal_X(k,j);   % 修改最短距离
                path(i,j) = path(k,j);                          % 记录最短途径
            end
        end
    end
end

  运行结果如下

在这里插入图片描述
  Fal_X 表示 节点 i 到节点 j 之间的最短距离。path 记录着 节点 i 到节点 j 的最短途径上, 结点 j 的上一结点 k。




4、第①小问 模型的建立和求解

  根据上述分析和数据处理,明显有些节点到所有交巡警服务平台的距离都不会小于 3km 。所以不能满足那个 “3分钟”(因为车速60km/h→1km/min) 的要求。

  所以,我们建立以下模型。并且在满足 “3分钟之内尽量到达” 的条件上,再使得对于少数不能感赶到的,我们尽量使其趋于最短时间。
在这里插入图片描述

说明:

  <1> x i j x_{ij} xij 是 0/1 变量,当 x i j = 1 x_{ij}=1 xij=1 表示第 j 个节点分配在第 i 个交巡警服务平台的管辖内;当 x i j = 0 x_{ij}=0 xij=0 表示第 j 个节点 没有 分配在第 i 个交巡警服务平台的管辖内。

  <2> L i j L_{ij} Lij 是 结点 i 到 结点 j 的最短距离。

  <3>目标函数为 ∑ i = 1 20 x i j L i j , ( j = 1 , 2 , . . . , 92 ) \sum \limits_{i=1}^{20} x_{ij}L_{ij} ,(j=1,2,...,92) i=120xijLij,(j=1,2,...,92) 表示第 j 个节点距离第 i 个交巡警服务平台距离最近时,就将第 j 个节点分配在第 i 个交巡警服务平台的管辖区域内。并 尽量 使交巡警平台 能满足 “3分钟” 的那个条件。

  <4> ∑ i = 1 20 x i j = 1 , ( j = 1 , 2 , . . . , 92 ) \sum \limits_{i=1}^{20} x_{ij}=1 ,(j=1,2,...,92) i=120xij=1,(j=1,2,...,92) 表示每个节点(总共92个) 只能归属于一个交巡警服务平台的管辖。

  <5> ∑ j = 1 92 x i j ≥ 1 , ( i = 1 , 2 , . . . , 20 ) \sum \limits_{j=1}^{92} x_{ij} ≥1,(i=1,2,...,20) j=192xij1,(i=1,2,...,20) 表示每个交巡警服务平台至少管辖一个节点。【也就是大家也都得干干活】

  用 lingo 软件对其进行求解:


model:
sets:  ! 集合段;
R/1..92/;
S/1..20/;  
T/1..92/;
JuZhen(S,T): x,L; 
BigJuZhen(R,T):Fal_X;
endsets

data:  ! 数据段;
Fal_X = @ole('C:\Users\ASUS\Desktop\CSDN\6交巡警\cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls','data_1');    ! 向Excel读数据
@ole('C:\Users\ASUS\Desktop\CSDN\6交巡警\cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls','结果') = X;         ! 向Excel写数据
enddata

calc:  ! 计算段;
@for( R(i)| (i #LE# 20): @for( T(j): L(i,j) = Fal_X(i,j) ) );
endcalc

min = @sum( JuZhen: x*L );      ! 目标函数;

! 以下3个函数是约束条件;
@for( T(j): @sum( S(i):x(i,j) ) = 1 );    

@for( S(i): @sum( T(j):x(i,j) ) >= 1 );  

@for( JuZhen: @bin(x)  );      ! 01变量
end

  运行结果如下:(已转到 Excel )

在这里插入图片描述


  后续的可视化处理,即用 matlab 软件对其进行画图:


x = xlsread('cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls',7,'A1:CN20')
table_2 = struct([]);   % 结构体初始化
for i=1:20
    table_2(i).PingTai = [];
end

for i=1:20
    for j=1:92
        if x(i,j) == 1
            table_2(i).PingTai = [ table_2(i).PingTai,j ];  % 将 x 的数据读到结构体
        end
    end  
end
Save = zeros(20,2);  % 初始化
for i=1:20
    Z = cell2mat( struct2cell( table_2(i) ) );  % 转换操作:结构体 → 元胞数组 → 矩阵
    biggest_dis = 0 ;
    x0 = Table(i).X;
    y0 = Table(i).Y;
    for j=1:length(Z)
        x1 = Table(Z(j)).X;
        y1 = Table(Z(j)).Y;
        dis = sqrt( ( x0 - x1 ).^2 + ( y0 - y1 ).^2 );
        if dis > biggest_dis
            Save(i,1) = dis;   % 记录最长时间(也就是最短距离,因为车速是1km/h)
            Save(i,2) = j;     % 记录交通路口编号
        end
    end
end
xlswrite('cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls',Save,9,'M2:N22');
hold on;
new_line = cell2mat( struct2cell( table_2(1) ) );  % 结构体 → 元胞数组 → 矩阵
for j=1:length(new_line)   % 给第一个交巡警服务平台画管辖区域
    xx = [ Table(1).X,Table(new_line(j)).X ] ;
    yy = [ Table(1).Y,Table(new_line(j)).Y ] ;
    scatter(xx(:,2),yy(:,2),10,'filled','m');
    plot(xx,yy,'m-.','LineWidth',1)
end

new_line = cell2mat( struct2cell( table_2(2) ) );
for j=1:length(new_line)   % 给第二个交巡警服务平台画管辖区域
    xx = [ Table(2).X,Table(new_line(j)).X ] ;
    yy = [ Table(2).Y,Table(new_line(j)).Y ] ;
    scatter(xx(:,2),yy(:,2),10,'filled','r');
    plot(xx,yy,'r-.','LineWidth',1)
end

  运行结果如下:(注:只画出两个交巡警服务平台的管辖区域,即 粉色 和 红色)

在这里插入图片描述




5、第②小问 模型的建立和求解

  ☁️ ☁️ ☁️

  首先,我们用 matlab ,统计出13个节点(A区的13个出入口)和20个交巡警服务平台的距离表。


Model_2 = zeros(20,13);
t = 1;
ans = [];
for i=1:92
    if Table(i).A_Qu_ChuRuKou == 1
        for j = 1:20
            Model_2(j,t) = Fal_X(j,i);  % 获取对应的最短距离
        end
        t = t + 1;
        ans = [ans , i];    % 获取编号
    end
end
Model_2
ans
xlswrite('cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls',Model_2,9);

  运行结果如下:(已用 Excel 作处理)

在这里插入图片描述
  :①表中的单位是mm。(比例尺1:100000)

    ②最上面一排项目中,括号外是打表的序号,括号内是A区13个出入口的编号。

  接着,考虑到题目中所说的 “快速全面封锁” 这点,而时间的长短是由 所有距离中的最长距离 决定的,故不能建立 所有路程之和最短的模型,而要建立一个 时间最短模型。因为制约时间的并非是所有路程之和而是其中的最长路径。

  所以建立的模型如下:

在这里插入图片描述

说明:
  <1> L0 为 13个A区出入口分别到20个服务平台的 13 个最短距离中的 最大值。

  <2> ∑ j = 1 13 x i j ≤ 1 , ( i = 1 , 2 , . . . , 20 ) \sum \limits_{j=1}^{13} x_{ij}≤1 ,(i=1,2,...,20) j=113xij1,(i=1,2,...,20) 表每个交巡警平台只能封锁一个节点。

  <3> ∑ i = 1 20 x i j = 1 , ( j = 1 , 2 , . . . , 13 ) \sum \limits_{i=1}^{20} x_{ij}=1 ,(j=1,2,...,13) i=120xij=1,(j=1,2,...,13) 表示每个A区出入口(总共13个) 都必须有交巡警进行封锁。

  然后,我们用 lingo 编程求解。


model:

sets:  
S/1..20/;
T/1..13/;
link(S,T): x,L;
endsets

data:
L = @ole('C:\Users\ASUS\Desktop\CSDN\6交巡警\cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls','dis');  !读数据;
@ole('C:\Users\ASUS\Desktop\CSDN\6交巡警\cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls','TRY') = x;  !写数据;
L0 = 57.01;    !从Excel直接获取的数据;
enddata 

min = L1;

L1 >= L0;

@for( S(i): @SUM( T(j): x(i,j) ) <= 1 );

@for( T(j): @SUM( S(i): x(i,j) ) = 1 );

@for( link: x*L <= L1 );

@for( link: @bin(x) );

end

  运行结果如下:(已用 Excel 作处理)

在这里插入图片描述
  :lingo解得的 L1 = 80.15,故最长的路程为第29个节点,但只有 8015m ,并且满足将各个节点全部封锁的条件。所以从 时间的角度 来看,该模型是值得肯定的。【其实这里有一个 对比方案 ,但我没写,那个对比方案是 以13个A区出入口20个服务平台的距离和最小 为规划目标。但这个对比方案的模型计算结果没 80.15 好。】




三、总结:

  在做 优化(/规划)模型 的题目时,一般用 lingo 要方便一点。但在此之前,我们一般要进行一下 数据预处理

  之后,我们要 找到并建立 以下三个东西:

    (1)目标函数
    (2)决策变量
    (3)约束条件

  最后,我没有详细地阐述其原理,只阐述了有什么用、怎么用。详细原理可以参考本文最后的 参考附录 📚 📚 📚。


四、参考附录:

[1] 《2011高教社杯全国大学生数学建模竞赛赛题 B题 数据》
链接: http://www.mcm.edu.cn/html_cn/node/a1ffc4c5587c8a6f96eacefb8dbcc34e.html.

[2] 《数学建模之Lingo基础知识与应用》
链接: https://blog.csdn.net/sunyueqinghit/article/details/81708836.

[3] 《LINGO学习笔记01》
链接: https://blog.csdn.net/Temmie1024/article/details/108863568.

[4] 《算法篇——详解 弗洛伊德(Folyd)算法【用 C/C++ 和 matlab 实现】》
链接: https://blog.csdn.net/Wang_Dou_Dou_/article/details/119485893?spm=1001.2014.3001.5501.

数学建模系列文章——总结篇《数模美一国一退役选手的经验分享[2021纪念版]》.


纯手码字 码图不易,多多支持 🙈🙈

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐