数学建模 铺路问题的最优化模型

铺路问题的最优化模型

摘要

本文采用了两种方法，一种是非线性规划从而得出最优解，另一种是将连续问题离散化利用计算机穷举取最优的方法。

根据A地与B地之间的不同地质有不同造价的特点，建立了非线性规划模型和穷举取最优解的模型，解决了管线铺设路线花费最小的难题。

问题一：在本问题中，我们首先利用非线性规划模型求解，我们用迭代法求出极小值（用Matlab实现），计算结果为总费用最小为748.6244万元，管线在各土层中在东西方向上的投影长度分别为15.6786km，3.1827 km，2.1839 km，5.8887km，13.0661km。然后，我们又用穷举法另外建立了一个模型，采用C语言实现，所得最优解为最小花费为748.625602万元，管线在各土层中在东西方向上的投影长度分别为15.70km,3.20km,2.20km,5.90km,13.00km。

问题二：本问题加进了一个非线性的约束条件来使转弯处的角度至少为160度，模型二也是如此。非线性规划模型所得计算结果为最小花费为750.6084万元，管线在各土层中在东西方向上的投影长度分别为14.4566km,4.3591km,2.5984km,6.5387km,12.0472km。遍历模型所得最优解为最小花费为750.821154万元，管线在各土层中在东西方向上的投影长度分别为14.10km,4.30km, 2.70km,6.70km,12.20km。

问题三：因为管线一定要经过一确定点P，我们将整个区域依据P点位置分成两部分，即以A点正东30km处为界，将沙土层分成两部分。非线性规划模型最小花费为752.6432万元，管线在各土层中在东西方向上的投影长度分别为21.2613km,3.3459km,2.2639km,3.1288km,2.4102km,7.5898km。遍历模型最小花费为752.649007万元，管线在各土层中在东西方向上的投影长度分别为21.30km,3.30km,2.30km,3.10km,2.40km,7.60km。

关键词：非线性规划逐点遍历穷举法

一．问题重述

准备在A地与B地之间修建一条地下管线，B地位于A地正南面26km和正东40km 交汇处，它们之间有东西走向岩石带。地下管线的造价与地质特点有关，下图给出了整个地区的大致地质情况，显示可分为三条沿东西方向的地质带，其宽度分别为：沙土地质带宽C1，C5；沙石地质带宽C2；沙石土地质带宽：C4；岩石地质带宽C3。

在给定三种地质条件上每千米的修建费用的情况如下：

试解决以下几个问题：

(1) 图中直线AB显然是路径最短的，但不一定最便宜；而路径ARSB过岩石和沙石的路径最短，但是否是最好的路径呢？试建立一个数学模型，确定最便宜的管线铺设路线。（若C1=6，C2=4，C3=5，C4=6，C5=5，确定最便宜的管线铺设路线。）

160，确定最便宜的管线铺设路线。

(2) 铺设管线时，如果要求管线转弯时，角度至少为0

(3) 铺设管线时，如果要求管线必须通过位于沙石地质带或岩石地质带中的某一已知点P（位于A地正南面18km和正东30km交汇处）时，确定最便宜的铺设路线。

二.模型假设

1、修建费用仅与管线长度和不同地质的造价有关，不含其他费用；

2、在无特殊要求情况下，管线可以向任意方向延伸；

3、不考虑管线宽度；

4、所有管线都铺设在同一水平面上；

三.符号说明

()x f 为修建总费用

1x 为管线与沙土层1c 中东西方向上的投影长度

2x 为管线与沙石层2c 中东西方向上的投影长度

3x 为管线与岩石层3c 中东西方向上的投影长度

4x 为管线与沙石土层4c 中东西方向上的投影长度（在问题三中指在过P 点的东西方

向的直线上的P 点以西的投影长度）

5x 为管线与沙土层5c 中东西方向上的投影长度（在问题三中指在过P 点的东西方向

的直线上的P 点以东的投影长度）

6x 为管线与沙土层5c 中东西方向上的投影长度 1p 为沙土层1c 每千米的修建费用 2p 为沙石层2c 每千米的修建费用

3p 为岩石层3c 每千米的修建费用 4p 为沙石土层4c 每千米的修建费用

5p 为沙土层5c 每千米的修建费用(在问题三中指在沙石土层每千米的修建费用) 6p 为问题三中沙土层6c 每千米的修建费用

4c 在问题一、二中指沙石土层的宽度，在问题三中指沙石土层P 点以上的半层的宽

度

5c 在问题一、二中指沙石土层的宽度，在问题三中指沙石土层P 点以下的半层的宽

度

6c 问题三中最下面的沙土层的宽度

四.问题分析

4．1 问题一:

本问题主要围绕由A 点到B 点铺设管线展开，要求花费最少。根据不同地质条件的花费，确定在某一土层中铺设管线的长度。我们采用了两种方法求得最少的花费，分别为非线性规划模型和逐点遍历模型。

4．1．1 方案一

我们首先利用非线性规划求解，可以得出一个关于工程总造价的目标函数f(x)，而且可知f(x)在整个区域连续且可微，f(x)符合在某一点有局部极小点的条件。因此我们用迭代法求出极小值（用Matlab 实现），我们分别选用了几组不同的初始值来保证所得到的极小值也是整个区域上的最小值。

4．1．2 方案二

我们又用穷举法另外建立了一个模型，用来确保模型一的结果是最小值，采用C 语言实现，我们先在每两种不同地质间的交界线上每隔0.1km 确定一个点，然后每条交界线都任取一点，连线，得出一条路径。之后将每一条可能的路径都遍历一遍，将最小值和对应的点保存，得出结果。

4．2 问题二

本问题与问题一相比，增加了约束条件“要求管线转弯时，角度至少为0

160”，我们在问题一所建立的两种模型的基础上均增加相应约束条件，通过求出管线转弯处的管线角度的正切值，并利用反正切函数得出管线角度，从而对管线的铺设方向加以限制，得出最少花费的管线铺设线路。

4．3 问题三

本问题要求铺设管线一定要经过一确定点P ，因此可以将此问题分为两步，即从A 到P 的路径为第一步，从P 到B 的路径为第二步。因为从A 到P 的路径选择及其花费与从P 到B 的路径选择及其花费无关，所以求出第一步从A 到P 的最优解，以及第二步求从P 到B 的最优解，这两的最优解之和便为整个管线铺设的最优解。

五.模型建立与求解

5．1 问题一

5．1．1 方案一.

根据题意，在第i 个土层中的管线长度为

22i i i c x d +=

所以，在该层中的修建花费为

i i i d p S ?=

则总花费为

∑=+?=5

1

22)(i i i i c x p x f

因此得到目标函数

∑=+?=5

122)(:min i i i i c x p x f

然后所要修建的地区为A 地正南面26km 和正东40km 所表示的区域，在每个土层中管线在东西方向的投影长度应大于0km 小于40km ，且所有土层中管线在东西方向上的投影长度之和小于40km ，因此可确定约束条件：

400≤≤i x

405

1

=∑=i i

x

运用MATLAB 软件编程，得到计算结果为总费用最小为748.6244万元，管线在各土层中在东西方向上的投影长度分别为15.6786km ，3.1827 km ，2.1839 km ，5.8887km ，13.0661km 。

5．1．2方案二

先在每两种不同土层的交界线上每隔0.1km 确定一个点，然后在每条交界线上都任取一点，并连线，得出一条可能路径。再将每一条可能的路径按公式

∑=+?=5

122)(i i i i c x p x f

逐一计算花费，找到花费的最小值和其对应的点，确定最优路径。在此方案中，采用C 语言编程进行遍历，所得最优解为最小花费为748.625602万元，管线在各土层中在东西方向上的投影长度分别为15.70km,3.20km,2.20km,5.90km,13.00km 。

5．2 问题二

本题也为确定最便宜的管线铺设路线，所以与问题一有相同的目标函数

∑=+?=5

122)(:min i i i i c x p x f

及约束条件：

400≤≤i x ;

405

1

=∑=i i

x

;

根据本题中所要求的管线转弯角度大于0

160，利用管线在各土层中在东西方向上的

投影长度与相应土层宽度得出管线转弯所形成的角的正切值，即i

i

c x ，再利用反正切函数算出具体角度。由此得到新的约束条件：

0arctan arctan 2

18070

11≤???

?

?????? ??+???? ??--?++i i i i c x c x ππ;

4,3,2,1=i

0180110

arctan arctan 211≤?-???

? ??+???? ??-++ππ

i i i i

c x c x ; 4,3,2,1=i

在问题一建立的模型的基础上，依据本题中新增非线性约束条件，建立新的模型，

利用MATLAB 编程，所得计算结果为最小花费为750.6084万元，管线在各土层中在东西方向上的投影长度分别为14.4566km,4.3591km,2.5984km,6.5387km,12.0472km 。

利用相同的约束条件，利用C 语言编程遍历，所得最优解为最小花费为750.821154万元，管线在各土层中在东西方向上的投影长度分别为14.10km,4.30km, 2.70km,6.70km,12.20km 。

5．3 问题三

根据本题中管线必须通过已知点P （位于A 地正南面18km 和正东30km 交汇处）的约束条件，我们将整个区域依据P 点位置分成两部分，即以A 点正东30km 处为界，将沙土层分成两部分，使整个修建区域变成6个土层。在问题一所建立的模型上加以改进，使目标函数变为：

∑=+?=6

122)(:min i i i i c x p x f

并将约束条件改为：

300≤≤i x ;

4,3,2,1=i ;100≤≤i x

6,5=i

304

1=∑=i i

x

;

106

5

=∑=i i

x

;

利用非线性规划模型，MATLAB 编程，所得最优解为:最小花费为752.6432万元，

管线在各土层中在东西方向上的投影长度分别为21.2613km,3.3459km,2.2639km, 3.1288km,2.4102km,7.5898km

利用遍历模型，C 语言编程，所得最优解为：最小花费为752.649007万元，管线在各土层中在东西方向上的投影长度分别为21.30km,3.30km,2.30km,3.10km,2.40km, 7.60km 。

六.模型的评价与改进

对于模型一，存在的缺点是用Matlab中的fmincon函数所求最优解可能只是局部最优解，必须代入几组不同的初始迭代值，来确定所求解为全局最优解，但仍有可能遗漏。对于模型二，缺点是精度不够小，当精度取到0.1时，计算机要用几分钟才能得出结果，精度更小时所需时间更长，而且在问题二中，模型二因为精度太低而跳过了最优路径，所以与模型一所得结果有一定差距。但模型二的结果可用于检验模型一所得结果是否为最优解。

七.参考文献

[1] 曹卫华郭正，最优化设计方法及 MATLAB的实现，北京：化学工业出版社，2005.

[2] 孙辉吴润秀，C语言程序设计，北京：中国铁道出版社，2007.

[3] 李尚志，数学建模竞赛教程，南京：江苏教育出版社，1996.

[4] 何文章，大学数学实验，哈尔滨：哈尔滨工程大学出版社，2000.

[5] 边馥萍，数学模型方法与算法，北京：高等教育出版社，2005.

[6] 姜启源，数学模型，北京：高等教育出版社，1993.

附件

附件一：

问题一的Matlab语言源代码：

先建立目标函数的M文件：

function f=fun(x);

f=12*sqrt(x(1)^2+36)+18*sqrt(x(2)^2+16)+28*sqrt(x(3)^2+25)+16*sqrt(x(4)^2

+36)+12*sqrt(x(5)^2+25)

主程序：

x0=[0;0;0;0;0];

A=[1 0 0 0 0;0 1 0 0 0;0 0 1 0 0;0 0 0 1 0;0 0 0 0 1];

b=[40 40 40 40 40];

Aeq=[1 1 1 1 1];beq=[40];

[x,z]=fmincon('fun',x0,A,b,Aeq,beq)

运算结果：

x =

15.6786

3.1827

2.1839

5.8887

13.0661

z =

748.6244

附件二：

问题一的C语言源代码：

#include

#include

#define JD 0.1 //取精度为0.1

main()

{

double x[6]={0,0,0,0,0,40},s=0,min=100000,best[5]; //s用来临时保存每一条路径的价格，x用来临时保存每一条路径

int i,j,k,l,n,c[5]={6,4,5,6,5},p[5]={12,18,28,16,12},D; //min 用来记录最小价格，best用来记录最佳路径

D=40/JD; //D为每一条交界线上可分的段数

for (i=0;i

{

for(j=i;j

{

for(k=j;k

{

for(l=k;l

{

for(n=1;n<6;n++)

{

s+=sqrt((x[n]-x[n-1])*(x[n]-x[n-1])+c[n-1]*c[n-1])*p[n-1];

}

if(s

{ //则将其存在min中，同时路径保存

min=s; //在best中。

for(n=0;n<5;n++)

best[n]=x[n];

}

s=0;

x[4]+=JD;

}

x[3]+=JD;

x[4]=x[3];

}

x[2]+=JD;

x[3]=x[2];

}

x[1]+=JD;

x[2]=x[1];

}

printf("最少花费为%f\n",min);

for(i=0;i<5;i++)

printf("在第%d条交界线上的位置为%f\n",i,best[i]);

}

运算结果：

最少花费为748.625602

在第0条交界线上的位置为0.000000

在第1条交界线上的位置为15.700000

在第2条交界线上的位置为18.900000

在第3条交界线上的位置为21.100000

在第4条交界线上的位置为27.000000

附件三：

问题二的Matlab语言源代码：

先建立目标函数的M文件：

function f=fun(x);

f=12*sqrt(x(1)^2+36)+18*sqrt(x(2)^2+16)+28*sqrt(x(3)^2+25)+16*sqrt(x(4)^2 +36)+12*sqrt(x(5)^2+25)

然后建立非线性约束条件的M文件：

function[c,ceq]=mycon(x)

c(1)=pi*70/180-pi*0.5+atan(x(1)/6)-atan(x(2)/4);

c(2)=pi*70/180-pi*0.5+atan(x(2)/4)-atan(x(3)/5)

c(3)=pi*70/180-pi*0.5+atan(x(3)/5)-atan(x(4)/6);

c(4)=pi*70/180-pi*0.5+atan(x(4)/6)-atan(x(5)/5);

c(5)=pi*0.5-atan(x(1)/6)+atan(x(2)/4)-pi*110/180;

c(6)=pi*0.5-atan(x(2)/4)+atan(x(3)/5)-pi*110/180;

c(7)=pi*0.5-atan(x(3)/5)+atan(x(4)/6)-pi*110/180;

c(8)=pi*0.5-atan(x(4)/6)+atan(x(5)/5)-pi*110/180;

ceq=[];

主程序：

x0=[0.1;0.1;0.1;0.1;0.1];

A=[1 0 0 0 0;0 1 0 0 0;0 0 1 0 0;0 0 0 1 0;0 0 0 0 1];

b=[40,40,40,40,40];

Aeq=[1 1 1 1 1];beq=[40];

[x,z]=fmincon('fun',x0,A,b,Aeq,beq,[],[],'mycon')

运算结果：

x =

14.4566

4.3591

2.5984

6.5387

12.0472

z =

750.6084

附件四

问题二的C语言源代码：

#include

#include

#define JD 0.1 //取精度为0.1

int tj(double a,double b,double c1,double d,double e) //满足转弯处至少为160度的条件

{

int n,f,c[5]={6,4,5,6,5};

double x[6]={0,a,b,c1,d,e};

for(n=1;n<5;n++)

{

if(3.1416*0.5-atan((x[n]-x[n-1])/c[n-1])+atan((x[n+1]-x[n])/c[n])>=1.22 1731&&3.1416*0.5-atan((x[n]-x[n-1])/c[n-1])+atan((x[n+1]-x[n])/c[n])<=1.9 19863)

f=1;

else

return 0;

}

return f;

}

main()

{

double x[6]={0,0,0,0,0,40},s=0,min=100000,best[5]; //s用来临时保存每一条路径的价格，x用来临时保存每一条路径

int i,j,k,l,n,c[5]={6,4,5,6,5},p[5]={12,18,28,16,12},D; //min 用来记录最

小价格，best用来记录最佳路径

D=40/JD; //D为每一条交界线上可分的段数

for (i=0;i

{

for(j=i;j

{

for(k=j;k

{

for(l=k;l

{

for(n=1;n<6;n++)

{

s+=sqrt((x[n]-x[n-1])*(x[n]-x[n-1])+c[n-1]*c[n-1])*p[n-1];

}

if(s

{ //且满足转弯处不小于160度的条件，

min=s; //则将其存在min中，同时路径保存

for(n=0;n<5;n++) //在best中。

best[n]=x[n];

}

s=0;

x[4]+=JD;

}

x[3]+=JD;

x[4]=x[3];

}

x[2]+=JD;

x[3]=x[2];

}

x[1]+=JD;

x[2]=x[1];

}

printf("最少花费为%f\n",min);

for(i=0;i<5;i++)

printf("在第%d条交界线上的位置为%f\n",i,best[i]);

}

运算结果：

最少花费为750.821154

在第0条交界线上的位置为0.000000

在第1条交界线上的位置为14.100000

在第2条交界线上的位置为18.400000

在第3条交界线上的位置为21.100000

在第4条交界线上的位置为27.800000

附件五：

问题三的Matlab语言源代码：

先建立目标函数的M文件：

function f=fun3(x);

f=12*sqrt(x(1)^2+36)+18*sqrt(x(2)^2+16)+28*sqrt(x(3)^2+25)+16*sqrt(x(4)^2

+9)+16*sqrt(x(5)^2+9)+12*sqrt(x(6)^2+25)

主程序：

x0=[0;0;0;0;0;0];

A=[1 0 0 0 0 0;0 1 0 0 0 0;0 0 1 0 0 0;0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1];

b=[30 30 30 30 10 10];

Aeq=[1 1 1 1 0 0;0 0 0 0 1 1];beq=[30;10];

[x,z]=fmincon('fun3',x0,A,b,Aeq,beq)

运算结果：

x =

21.2613

3.3459

2.2639

3.1288

2.4102

7.5898

z =

752.6432

附件六：

问题三的C语言源代码：

#include

#include

#define JD 0.1 //取精度为0.1

main()

{

double x[7]={0,0,0,0,30,30,40},s=0,min,min1=100000,min2=100000,best[6]; //s用来临时保存每一条路径的价格，x用来临时保存每一条路径

int i,j,k,n,c[6]={6,4,5,3,3,5},p[6]={12,18,28,16,16,12},D; //min 用来记

录最小价格，best用来记录最佳路径

D=30/JD; //D为每一条交界线上可分的段数

for (i=0;i

for(j=i;j

{

for(k=j;k

{

for(n=1;n<5;n++)

{

s+=sqrt((x[n]-x[n-1])*(x[n]-x[n-1])+c[n-1]*c[n-1])*p[n-1];

}

if(s

{ //则将其存在min中，同时路径保存

min1=s; //在best中。

for(n=0;n<5;n++)

best[n]=x[n];

}

s=0;

x[3]+=JD;

}

x[2]+=JD;

x[3]=x[2];

}

x[1]+=JD;

x[2]=x[1];

}

printf("从A到P最少花费为%f\n",min1);

for(i=0;i<4;i++)

printf("在第%d条交界线上的位置为%f\n",i,best[i]);

//第二部分

best[4]=30;

s=0;

for(i=0;i

{

for(n=5;n<7;n++)

{

s+=sqrt((x[n]-x[n-1])*(x[n]-x[n-1])+c[n-1]*c[n-1])*p[n-1];

}

if(s

{ //则将其存在min中，同时路径保存

min2=s; //在best中。

best[5]=x[5];

}

s=0;

x[5]+=JD;

}

printf("从P到B最少花费为%f\n",min2);

for(i=4;i<6;i++)

printf("在第%d条交界线上的位置为%f\n",i,best[i]);

min=min1+min2;

printf("总最少花费为%f\n",min);

}

运算结果：

从A到P最少花费为582.011987

在第0条交界线上的位置为0.000000

在第1条交界线上的位置为21.300000

在第2条交界线上的位置为24.600000

在第3条交界线上的位置为26.900000

从P到B最少花费为170.637021

在第4条交界线上的位置为30.000000

在第5条交界线上的位置为32.400000

总最少花费为752.649007