当前位置:首页 >> 理学 >>

2012高教社杯全国大学生数学建模竞赛D题全国一等奖论文


2012 高教社杯全国大学生数学建模竞赛







我们仔细阅读了中国大学生数学建模竞赛的竞赛规则. 我们完全明白,在竞赛开始后参赛队员不能以任何方式(包括电话、电子邮件、网 上咨询等)与队外的任何人(包括指导教师)研究、讨论与赛题有关的问题。 我们知道,抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的 资料(包括网上查到的资料) ,必须按照规定的参考文献的表述方式在正文引用处和参 考文献中明确列出。 我们郑重承诺,严格遵守竞赛规则,以保证竞赛的公正、公平性。如有违反竞赛规 则的行为,我们将受到严肃处理。

我们参赛选择的题号是(从 A/B/C/D 中选择一项填写) : 我们的参赛报名号为(如果赛区设置报名号的话) : 所属学校(请填写完整的全名) : 参赛队员 (打印并签名) :1. 2. 3. 指导教师或指导教师组负责人 (打印并签名):

D

(此部分内容不便公开,见谅)

日期: 2012 年 9 月 10 日

赛区评阅编号(由赛区组委会评阅前进行编号) :

2012 高教社杯全国大学生数学建模竞赛

编 号 专 用 页

赛区评阅编号(由赛区组委会评阅前进行编号) :

赛区评阅记录(可供赛区评阅时使用) : 评 阅 人 评 分 备 注

全国统一编号(由赛区组委会送交全国前编号) :

全国评阅编号(由全国组委会评阅前进行编号) :

机器人避障问题
摘要
针对机器人避障问题, 本文分别建立了机器人从区域中一点到达另一点的避障的最 短路径、最短时间路径的非线性 0-1 整数规划模型。同时,本文为求带有 NP 属性的非 线性 0-1 整数规划模型, 构建了有效启发式算法, 利用 MATLAB 软件编程, 求得了 O→A、 O→B、O→C、O→A→B→A→C 的最短路径,同时得到了 O→A 的最短时间路径,求得的各 类最短路径均是全局最优。 针对区域中一点到达另一点的避障的最短路径问题,首先,本文证明了圆弧位置设 定在需要绕过障碍物的顶角上,且圆弧半径为 10 个单位时,能够使得机器人从区域中 一点到达另一点的行进路径最短;其次,本文将最短路径选择问题转化成了最短路径的 优选问题,根据避障条件,建立了具有较高普适性的避障最短路径的优化模型。为便于 求解,本文巧妙地将此优化模型转化成了以可行路径不与障碍物边界相交、不与圆弧相 交为约束条件, 以机器人从区域中一点达到另一点避障路径最短为目标的 0-1 规划模型; 再次,本文构建了两种有效的启发式算法,利用 MATLAB 软件编程求得了 O→A、O→B、O →C、 O→A→B→A→C 的最短路径, 最短路径长分别为 471.0372、 853.7001、 1088.1952、 2725.1596, 其中 O-->A 的最短路径为(0,0)→(70.5063,213.1405) →(75.975,219.1542) →(300,300),对应圆弧的圆心坐标为(80,210),O→B 的最短路径,对应圆弧的圆心坐 标:(60,300)、(150,435)、 (220、470) 、(220,530)、(150,600), O→C 经过的圆心: (410,100)、(230,60)、(720,520), (720,600) ,(500,200), O→A→B→C→O 经过的圆 心: (410,100), (230,60), (80,210), (220,530), (150,600), (270,680), (370,680), (430,680),(670,730),(540,730),(720,520),(720,600),(500,200)。 针对最短时间路径问题, 我们建立了从 o 点出发到任意目标点的 0-1 非线性整数规 划模型,同时针对题意要求,具体构建了从 o 点出发到 A 的最短时间路径的 0-1 非线性 整数规划模型,利用 LINGO 软件求解,获得了机器人从 o 点出发,到达 A 的最短时间路 径,求得最短时间路径下转弯半径为 12.9885 ,同时最短时间路径时间长为 94.2283 个单位。 相应圆弧的圆心坐标为 (82.1414, 207.9153),两切点坐标分别为(69.8045, 211.9779)、(77.7492,220.1387)。 本文确定路线思路循序渐进,先建立了有计算避障约束公式的普适性模型,再建 立了以不取相交点来简化 0-1 变量取值关系的简化模型;给出了二种启发式算法,最短 路径即最短时间路径具有一定可信度。同时第一个启发算法可以求得全局最优解,第二 个启发算法是针对问题的 NP 属性减少求解时间而构建的,两个算法都具有较重要的意 义。

【关键词】 机器人避障

最短路径

启发算法

0-1 规划模型

1

一、问题重述
在一个800× 800的平面场景图,在原点O(0, 0)点处有一个机器人,它只能在该平面 场景范围内活动。图中有12个不同形状的区域是机器人不能与之发生碰撞的障碍物,障 碍物的数学描述如下表:
编号 1 2 3 4 5 6 7 8 9 10 11 12 障碍物名称 正方形 圆形 平行四边形 三角形 正方形 三角形 长方形 平行四边形 长方形 正方形 正方形 长方形 (360, 240) (280, 100) (80, 60) (60, 300) (0, 470) (150, 600) (370, 680) (540, 600) (640, 520) (500, 140) 左下顶点坐标 (300, 400) 边长200 圆心坐标(550, 450),半径70 底边长140,左上顶点坐标(400, 330) 上顶点坐标(345, 210),右下顶点坐标(410, 100) 边长150 上顶点坐标(150, 435),右下顶点坐标(235, 300) 长220,宽60 底边长90,左上顶点坐标(180, 680) 长60,宽120 边长130 边长80 长300,宽60 其它特性描述

在图1的平面场景中,障碍物外指定一点为机器人要到达的目标点(要求目标点与 障碍物的距离至少超过10个单位)。规定机器人的行走路径由直线段和圆弧组成,其中 圆弧是机器人转弯路径。机器人不能折线转弯,转弯路径由与直线路径相切的一段圆弧 组成,也可以由两个或多个相切的圆弧路径组成,但每个圆弧的半径最小为10个单位。 为了不与障碍物发生碰撞, 同时要求机器人行走线路与障碍物间的最近距离为10个单位, 否则将发生碰撞,若碰撞发生,则机器人无法完成行走。

v ? v( ? ) ?

机器人直线行走的最大速度为 v0 ? 5 个单位/秒。机器人转弯时,最大转弯速度为

v0

法完成行走。

1 ? e10?0.1?

2

,其中 ? 是转弯半径。如果超过该速度,机器人将发生侧翻,无

请建立机器人从区域中一点到达另一点的避障最短路径和最短时间路径的数学模 型。对场景图中4个点O(0, 0),A(300, 300),B(100, 700),C(700, 640),具体计算: (1) 机器人从O(0, 0)出发,O→A、O→B、O→C和O→A→B→C→O的最短路径。 (2) 机器人从O (0, 0)出发,到达A的最短时间路径。 注:要给出路径中每段直线段或圆弧的起点和终点坐标、圆弧的圆心坐标以及机器 人行走的总距离和总时间。

二、问题分析
2.1 求取最短路径的分析 本问题要求机器人从区域中一点到达另一点的避障最短路径。 机器人只要做到转弯 时的圆弧半径最小为 10 个单位、与障碍物最近距离单时刻保持大于 10 个单位,那么可 行走的路径就有无数条,若想求得机器人从区域中一点到达另一点的避障最短路径,则 需要建立避障的最短路径模型,而建立避障的最短路径模型有一定难度。根据对问题的 分析,我们认为可以从简单做起,先确定小范围内最短路径条件,如圆弧位置的影响, 圆弧半径的大小,避免与障碍物碰撞条件等,通过确定最短路径条件来建立避障的最短 路径模型。对于最短路径的求取,我们可以通过确定穷举原则,利用穷举法来求解,当
2

然也可以通过构建启发式算法的进行求解。 2.2 最短时间路径的分析 对于要建立最短时间路径模型来说,我们容易知道影响的因素有直线行走速度、转 弯速度,同时还需要考虑使得最短时间路径条件,如圆弧位置(坐标)的影响,圆弧半 径的大小,避免与障碍物碰撞条件等。对于直线行进,我们希望行进速度越大越好,对 于机器人转弯时,转弯速度要有约束,要保证机器人不能发生侧翻。我们发现圆弧半径 的大小与转弯速度紧密相连,从转弯速度公式来分析,当转弯半径增大时,最大转弯速 度也增大,为在更短时间内行进到目标点,我们希望转弯速度为机器人的最大转弯速度 较好, 但有很大的可能是行进的路径不是最短的, 即行进路径有很大可能在增加。 于是, 我们需要做的工作是,在满足最短时间路径条件时,找到一个圆弧的坐标位置,同时确 定半径的大小,以求得最短时间路径。

三、模型假设
1.假设将机器人看成一个质点; 2.假设半径不变时,机器人在行进、转弯过程中能一直保持最大的速度; 3.假设启发算法是针对问题的 NP 属性减少求解时间而构建的。

四、符号说明
符号
Z
aij

含 避障最短路径



单位

备注

圆弧 i 切点到圆弧 j 切点的直线距离,即机器 人从圆弧 i 切点到圆弧 j 切点的直线路径

i, j ? 1...n, n ? 36, i ? j

bij

机器人从圆弧 i 行进至圆弧 j 切点时的转弯 路径

i, j ? 1...n, n ? 36, i ? j

xi yi

圆弧 i 的横坐标 圆弧 i 的纵坐标 圆弧对应圆心角 圆弧的半径 度 度

i ? 1..36
j ? 1..36

?
?

3

五、最短路径模型建立与求解
5.1 模型准备 5.1.1 确定圆弧位置与转弯半径 在建立机器人从区域中一点到达另一点的避障最短路径数学模型之前, 我们需要考 虑两个问题: 问题一:机器人从区域中一点到达另一点过程中,若中间有障碍物,则需要通过转 弯来绕过障碍物,那么,在转弯半径一定的情况下,怎样设定最佳圆弧位置,使得绕行 路径最短? 问题二:绕行路径是最短时,转弯半径的大小为多少? 针对考虑的问题一,我们取机器人从 O 到 A 点的行走过程来说明问题。在行走过程 中要求机器人行走线路与障碍物的最短距离为 10 个单位,圆弧(转弯)半径最小为 10 个单位。我们先令机器人转弯半径为 10 个单位,根据机器人行走过程中的要求,我们 易得两条极端的行走路径, 如图 1。 将路线 II 中圆弧 3 两切点线延长, 两延长交路线 I, 两交点处分别作半径为 10 个单位的圆弧,由此我们可得机器人从 O 到 A 点的行走时转 弯中心坐标的范围,如图 2 中四边形 abcd 。

图 1 两条极端路径 图 2 转弯中心坐标的范围 图 1 中路线 I 是理想化路线,机器人不能沿 800 ? 800 平面区域边界行走,800 ? 800 平面区域边界也可以看成是一个障碍物, 且有要求机器人行走线路与障碍物的最短距离 为 10 个单位,实际上作这样的处理并不会影响我们说明问题。 我们假设在平面中有 A(0, a) 和 O(? A,0) 两点,中间有一正方形的障碍物,将图 2 进 行转化,如图 3.

4

图 3 最短路径证明图 图 3 中, B, C... I 为切点, a, b, c, d 为圆弧圆心,四边形 abcd 为圆弧中心点的范围。 对于最佳圆弧位置确定,我们采用“覆盖法” 。我们容易知道,若路线 II 与 OA 构 成的区域 II 能够完整覆盖线 I 与 OA 构成的区域 I,即区域 I 属于区域 II,那么区域 II 的周长一定大于区域 I,否则。 图 3 中路线 I 与 OA 构成的区域 I 周长为直线段 OB ? CA 长度、圆弧 BC 长、 OA 长 之和,区域 I 周长 l1 为

l1 ? OB ? CA ? BC ? OA
机器人沿路线 I 的路径长 c1 可表示为

c1 ? OB ? CA ? BC
路线 II 与 OA 构成的区域 II 周长为直线段 OH ? FA 长度、 圆弧 FG 长、 长之和, OA 区域 II 周长 l 2 为

l2 ? OF ? GA ? FG ? OA
机器人沿路线 II 的路径长 c2 可表示为

c2 ? OF ? GA ? FG
显然我们知道区域 II 能够覆盖区域 I,即可得 l2 ? l1 ,进而可得到

c1 ? c2
同理,在圆弧中心点的范围任意取一点作为机器人转弯圆弧中点,并作路线 i ,再 将路线 i 与路线 I 做比较,可得到
c1 ? ci
5

由此,我们可得出结论:机器人从区域中一点到达另一点过程中,当圆弧位置设定 障碍物顶角上时,绕行路径最短,此时圆弧中心点坐标为障碍物顶角坐标。 针对考虑的问题二,为了更清晰说明绕行路径是最短时,转弯半径的大小为多少, 我们基于最小圆弧半径条件下使圆弧半径增大。为了保证机器人与障碍物不发生碰撞, 所以,需要保证大圆弧能够覆盖小圆弧对应圆的 1/4 圆弧。在设定好最佳圆弧位置情况 下,增加圆弧半径,比较最短路径的变化。假设圆弧半径为 R ( R ? 10 ) ,对应最短路 线如图 4。

(1)

(2)

图 4 圆弧半径为 R 最短路径 图 4(1)中 B 为两圆弧公共切点, C 为小圆弧切点( r ? 10 ) A 为大圆弧切点。 , 图 4(2)中 A 、 D 为大圆弧切点, B 、 C 为小圆弧切点。 针对图 4(1) ,根据“覆盖”思想与得出的结论,我们容易发现,当圆弧半径为 (2) R ( R ? 10 )时,行进路线与 OA 构成的区域显然是能够完全覆盖圆弧半径为 r 时构成的 区域,由此,可说明在圆弧位置设定为最佳的条件下,圆弧半径越小行进路径越短,而 圆弧半径最小为 10 个单位,进而说明圆弧的半径为 10 个单位时,绕行路径最短。 根据考虑的两个问题与证明结果,我们可得出结论:要使得机器人从区域中一点到 达另一点的行进路径最短,应使圆弧位置设定在需要绕过障碍物的顶角上最佳,此时圆 弧中心点坐标为障碍物顶角坐标,并且此时圆弧的半径为 10 个单位。 5.1.2 问题的转化 在模型准备中我们已得出要使得机器人从区域中一点到达另一点的行进路径最短, 应使圆弧位置设定在需要绕过障碍物的顶角上最佳, 此时圆弧中心点坐标为障碍物顶角 坐标,并且此时圆弧的半径为 10 个单位。因此,我们将 800 ? 800 的平面场景图进行处 理,处理原则有: 1.每个障碍物的顶角都设定一个圆弧; 2.圆弧坐标为障碍物顶点坐标; 3.圆弧的半径设定为 10 个单位; 4.给每一段圆弧从 2... n 标号,O 点标记为 1、36。 根据处理原则,得图 5。

6

图 5 处理后 800 ? 800 的平面场景图 在原问题中,若没进行确定圆弧位置与转弯半径以及平面场景的处理,原问题求解 将会很难,并且所有的转弯点均是未知,经处理后,我们将问题转化为在已知转弯点, 寻找合适的转弯点,使得路径最短,即我们将问题转化为了最短路径的优化问题。 5.2 避障最短路径模型的建立 问题转化为最短路径的优化问题后, 我们易知优化的目标函数为机器人在行进过程 中最短路径,通过 0-1 变量来选取经过的转弯点,因此可建立 0-1 规划模型。 5.2.1 目标函数 目标函数为机器人出区域中一点到达另一点的避障最短路径, 避障最短路径 Z 为行 进过程中直线路径与转弯路径之和,于是有
min Z ? ?? X ij (aij ? bij )
i ?1 j ?1 n n

(1)

其中, aij 为圆弧 i 切点到圆弧 j 切点的直线距离,即机器人从圆弧 i 切点到圆弧 j 切 点的直线路径; bij 为机器人从圆弧 i 行进至圆弧 j 切点时的转弯路径; X ij 为 0-1 变量, 表示若选择行走圆弧 i 后行走圆弧 j , X ij 为 1,否则为 0, i, j ? 1...n, n ? 36, i ? j 。

7

我们设 ( xik , yik ) 为圆弧 i 的切点坐标, 为各圆弧的编码, ? 1...36 , 为切点的次序, i k i
k ? 1,2 ,如 ( xi1 , yi1 ) 表示为圆弧 i 的第一个切点坐标。

机器人从圆弧 i 切点到圆弧 j 切点的直线路径 aij 为
aij ? ( xi 2 ? x j 2 ) 2 ? ( y i 2 ? y j1 ) 2 i, j ? 1..36, i ? j

(2)

机器人从圆弧 i 行进至圆弧 j 切点时的 bij 为
bij ? ?? i, j ? 1..36, i ? j

其中
( xi1 ? xi 2 ) 2 ? ( yi1 ? yi 2 ) 2 ? ? 2 ? ? 2 ? 2 ? 2 cos?

(3)

? 为转弯半径, ? ? 10 。
5.2.2 约束条件 1.避障条件的约束 在本问题中,障碍物边界可分为直线段与圆弧两种情况,针对障碍物边界不同的两 种情况,我们列出避障条件的约束: (1)避障约束条件 1 任意一条可行路径与所有障碍物的边界线段间的最短距离大于 10 个单位. 设平面内有两条线段 AB 和 CD , 点的坐标为 ( xi 2 , yi 2 ) , B 点的坐标为 ( x j1 , y j1 ) , A

C 点的坐标为 ( x p , y p ) , D 点的坐标为 ( xq , yq ) ,其中 A 、 B 可视为圆弧的切点坐标, C 、
D 为圆弧的切点坐标 AB 线段两圆弧的切点的连线,CD 视为障碍物的某一条边界线段。
m m 设 Ppq 是直线 AB 上的一点, Ppq 点的坐标 ( x m , y m ) 可以表示为 pq pq

? X ? xi 2 ? s ( x j1 ? xi 2 ) ? ? ?Y ? yi 2 ? s ( y j1 ? yi 2 ) ?
m m 当参数 0 ? s ? 1 时, Ppq 是线段 AB 上的点;当参数 s ? 0 时, Ppq 是 BA 延长线上的 m 点;当参数 s ? 1时, Ppq 是 AB 延长线上的点。

设 Qij 是直线 CD 上的一点, Qij 点的坐标 (U ,V ) 可以表示为
?U ? x P ? t ( xq ? x p ) ? ? ?Y ? y p ? s ( y q ? y p ) ?
Q Q 当参数 0 ? t ? 1时, ij 是线段 CD 上的点; 当参数 t ? 0 时, ij 是 DC 延长线上的点;

8

当参数 t ? 1 时, Qij 是 CD 延长线上的点。
m Ppq , Qij 两点之间的距离为
m Ppq Qij ? ( X ? U ) 2 ? (Y ? V ) 2

距离的平方为
m f ( s, t ) ? Ppq Qij ? ( X ? U ) 2 ? (Y ? V ) 2 2

? [( xi 2 ? x p ) ? s( x j1 ? xi 2 ) ? t ( xq ? x p )]2 ? [( yi 2 ? y p ) ? s( y j1 ? yi 2 ) ? t ( yq ? y p )]2
要求直线 AB , CD 之间的最短距离,也就是要求函数 f (s, t ) 的最小值。对 f (s, t ) 分 别求关于 s , t 的偏导数,并令偏导数为零: ? ?f ( s, t ) ? ?s ? 0 ? ?f ( s, t ) ? ?0 ? ?t 展开并整理后,得到下列方程组:9 ? [( x j1 ? xi 2 ) 2 ? ( y j1 ? yi 2 ) 2 ]s ? 2 ??[( x j1 ? xi 2 )( xq ? x p ) ? ( y j1 ? yi 2 )( yq ? y p )]t ? ?? ( xi 2 ? x j1 )( xi 2 ? x p ) ? ( yi 2 ? y j1 )( yi 2 ? y p ) ? ??[( x j1 ? xi 2 )( xq ? x p ) ? ( y j1 ? yi 2 )( yq ? y p )]t ? 2 2 ??[( xq ? x p ) ? ( yq ? y p ) ]t ?? ( x ? x )( x ? x ) ? ( y ? y )( y ? y ) i2 p q p i2 p q p ?
m 如果从这个方程组求出的参数 s , t 的值满足 0 ? s ? 1 , 0 ? t ? 1,说明 Ppq 点落在线 m 段 AB 上, Qij 点落在线段 CD 上,这时 Ppq Qij 的长度为
m Ppq Qij ? ( X ? U ) 2 ? (Y ? V ) 2

m 此时 Ppq Qij 就是线段 AB 与 CD 的最短距离。

如果从方程组求出的参数 s , t 的值不满足 0 ? s ? 1 , 0 ? t ? 1 ,说明不可能在线段
m m AB 内部找到一点 Ppq ,在线段 CD 内部找到一点 Qij ,使得 PpqQij 的长度就是线段 AB 与

CD 的最短距离。这时,还需要进行考虑的是:平面中一个点到一条线段的最短距离。
m 设编号为 m 的障碍物的圆心坐标 Ppq 和一条线段 AB , Ppq 点的坐标为 ( x m , y m ) , 设 m A pq pq

点的坐标为 ( xi 2 , yi 2 ) , B 点的坐标为 ( x j1 , y j1 ) 。 此时,直线 AB 的参数形式的方程为
9

? x ? xi 2 ? t ( x j1 ? xi 2 ) ? ? ? y ? yi 2 ? t ( y j1 ? yi 2 ) ?

直线上的点 ( x, y ) ,当参数 0 ? t ? 1时,是线段 AB 上的点;当参数 t ? 0 时,是 BA 延 长线上的点;当参数 t ? 1 时,是 AB 延长线上的点。
m 通过 Ppq 点,与直线 AB 垂直的平面方程为

( x j1 ? xi 2 )( x ? x m ) ? ( y j1 ? yi 2 )( y ? y m ) ? 0 pq pq
m m 下面求这个平面与直线 AB 的交点 Qij 。显然, Ppq Qij ? AB ,所以 Qij 点也是从 Ppq 点

向直线 AB 作垂线的垂足点。 将
? x ? xi 2 ? t ( x j1 ? xi 2 ) ? ? ? y ? yi 2 ? t ( y j1 ? yi 2 ) ?

代入平面方程,化简后解得
t? ( xi 2 ? x m )( xi 2 ? x j1 ) ? ( yi 2 ? y m )( yi 2 ? y j1 ) pq pq ( xi 2 ? x j1 ) 2 ? ( yi 2 ? y j1 ) 2
? X ? xi 2 ? t ( x j1 ? xi 2 ) ? ? ?Y ? yi 2 ? t ( y j1 ? yi 2 ) ?

然后,将上面得到的 t 的值代入直线方程,得到

其中, ( X , Y ) 为垂足点 Qij 的坐标。
m m 此时线段 Ppq Qij 的长度,也就是 Ppq 点到直线 AB 的垂直距离为
m Ppq Qij ? ( X ? x m ) 2 ? (Y ? y m ) 2 pq pq

m 如果前面求出的参数 t 的值满足 0 ? t ? 1, 说明垂足 Qij 点落在线段 AB 上, 这时 Ppq Qij m 的长度就是 Ppq 点到线段 AB 的最短距离。

如果前面求出的参数 t 的值满足 t ? 0 ,说明垂足 Qij 点不落在线段 AB 上,而是落在
m m m BA 的延长线上,这时 Ppq 点到线段 AB 的最短距离,就是 Ppq 点到 A 点的距离,即 Ppq A

? ( xi 2 ? x m ) 2 ? ( yi 2 ? y m ) 2 。 pq pq

如果前面求出的参数 t 的值满足 t ? 1 ,说明垂足 Qij 点不落在线段 AB 上,而是落在
m m 这时,Ppq 点到线段 AB 的最短距离, 就是 Ppq 点到 B 点的距离, Ppq B 即 m AB 的延长线上,

10

? ( x j1 ? x m ) 2 ? ( y j1 ? y m ) 2 。 pq pq

综上所述,即有平面内有两条线段最短距离为
d
ij qpm m ? qqpQij 0 ? s ? 1, 0 ? t ? 1 ? ?? m m ?min{ qqp A , qqp B } 否则 ?

因此, 我们可得到任意一条可行路径与所有障碍物的边界线段间的最短距离大于 10 个单位的约束条件

min ?d ij ?? 10 pqm
其中, i, j ? 1,2...n, n ? 36 , p, q ? 1,2...k , k ? 4 , m ? 1,2...m, m ? 12

(4)

(2)避障约束条件 2 对于标号为 2 的圆形障碍物,它与机器人行走路线的最近距离为 10 个单位,为使 满足与机器人行走路线的最近距离为 10 个单位要求,我们可以转化为圆形的圆心与障 碍物间的最近为 80 个单位的问题,即把圆形与直线线段的最短路径问题转化为点到直 线线段的的问题。 已知编号为 2 半径为 R ( R ? 70 )的圆形障碍物,边界并不是直线段,且该圆形障 碍物的圆形为 (550,450) 。 根据平面中一定点到一条线段的最短距离计算办法,我们可得

? O2Qij 0 ? s ? 1, 0 ? t ? 1 ? d ij ? ? ?min{ O2 A , O2 B } 否则 ? 因此,我们可得到圆形障碍物与平面内最近约束条件为
min ?dij ? ? R ? 10 i, j ? 1,2...n, n ? 36

(5)

2.各坐标点的约束 所有假设的坐标点都应在 800 ? 800 平面场景内,于是有
xi1, xi 2 , x j1 , yi1 , xi1 , yi 2 , y j1 , x p , x, y p , y p ? {0,800}

3.各圆弧点行进先后次序的约束 机器人从某个圆弧出发行至下一圆弧,先后次序的约束有
xi1 ? xi 2

yi1 ? yi 2
x j1 ? x j 2
y j2 ? y j2

(6)

xi 2 ? x j1

11

y i 2 ? x j1

其中, i, j ? 1..36, i ? j 。 5.2.3 避障最短路径模型 综合上述,建立避障最短路径模型:
min Z ? ?? X ij (aij ? bij )
i ?1 j ?1 n n

? ?Yij ? X ij ? ? xi1 ? xi 2 ? yi1 ? yi 2 ? ? y j2 ? y j2 ?x ? x j1 ? i2 ? y i 2 ? x j1 ? ? s.t ?bij ? ?? i, j ? 1..36, i ? j ? 2 2 2 2 2 ?( xi1 ? xi 2 ) ? ( yi1 ? yi 2 ) ? ? ? ? ? 2 ? cos? ? x , x , x , y , x , y , y , x , x, y , y ? (0,800) p p ? i1 i 2 j1 i1 i1 i 2 j1 p ?a ? ( x ? x ) 2 ? ( y ? y ) 2 i, j ? 1..36, i ? j i2 j2 i2 j1 ? ij ? X ij ? 0或1(i, j ? 1,2,...n, n ? 36, i ? j) ? ?min d ij ? 10 pqm ? ?d ij ? R ? 10 i, j ? 1,2...n, n ? 36 ?

? ?

5.3 基于避障条件转化下的最短路径简化模型 5.3.1 避障条件的转化 我们将任意一条可行路径与所有障碍物的边界线段间的最短距离大于 10 个单位转 化为两直线不相交、与 2 号障碍物圆周均不相交、切点范围的的控制。 出于简化避障条件的 0-1 变量取值关系,我们进行符号补充:
?1 选择走圆弧i后行进圆弧j i, j ? 1,2,..., n, n ? 36 X ij 为 0-1 变量,表示为 X ij ? ? ?0 否则
' ' y ij 表示圆弧 i 与圆弧 j 的圆心( x j , y j )所构成的线段的纵坐标, yij ?

xi ? x j yi ? y j

x,

x ? ( xi , x j ) 。

12

y m 表示编号为 m 的区域内一顶点( x p , y p )与相邻顶点( xq , y q ) ,所构成的线段的 pq
m 纵坐标, y pq ?

x p ? xq y p ? yq

x , x ? ( x p , xq ) 。当 m ? 2 时,顶点表示以(550,450)为圆点,80

个单位为半径的圆上的点,相邻顶点坐标则表示为以圆点对称的点。
? ?0 Yij 为 0-1 变量,表示为 Yij ? ? ?1 ?
' yij ? y m 存在唯一解 x ? ( x p , xq ) pq

否则



我们将避障约束条件替换为
Yij ? X ij

(7)

5.3.2 基于避障条件转化下的最短路径简化模型 于是,我们可得避障最短路径简化模型
min Z ? ?? X ij (aij ? bij )
i ?1 j ?1 n n

? ?Yij ? X ij ?x ? x i2 ? i1 ? y i1 ? y i 2 ? ? y j2 ? y j2 ?x ? x j1 ? i2 ? s.t ? y i 2 ? x j1 ? i, j ? 1..36, i ? j ?bij ? ?? ?( x ? x ) 2 ? ( y ? y ) 2 ? ? 2 ? ? 2 ? 2 ? 2 cos? i2 i1 i2 ? i1 ? x i1 , x i 2 , x j1 , y i1 , x i1 , y i 2 , y j1 , x p , x, y p , y p ? (0,800) ? ?a ij ? ( x i 2 ? x j 2 ) 2 ? ( y i 2 ? y j1 ) 2 i, j ? 1..36, i ? j ? ? X ij , Yij ? 0或1(i, j ? 1,2,...n, n ? 36, i ? j) ?

5.4 避障最短路径简化模型的求解 5.4.1 基于避障条件转化之下的最短路径的启发算法 对于避障最短路径的数学模型,是非线性 0-1 整数规划,具有 NP 属性。对于具体 的区域与障碍物,通过巧妙地将避障条件转化为相交约束,同时结合任意两点间最短距 离的 floyd 算法,我们构建出了能够对此问题规模求得全局最优解的较好算法。具体算 法思想是:针对此模型的 NP 属性特征,即此问题搜索可行域的规模越发,运行时间倍 增的不足,我们采用逐步缩减可行域的办法,将达到避障条件要求的所有可行路径搜索 出。 同时, 通过用两切点标识圆弧路径的办法, 将一个圆弧路径扩充为两个切点的表示, 建立出具有扩充点的邻接矩阵,于是将问题转化为赋权图上两点的最小距离问题,采用 Floyd 算法,即可得出避障最短路径。具体算法步骤为:
13

Setp1:得出所有满足避障条件的可行路径,分三步逐步进行。 (1)对于 n 个圆弧路径,任意两个圆弧间构造切线段及出发点到各圆弧的切线段,由 此形成初始的路径。判断每条路径中的直线段与各障碍物的边界线段是否相交,判断方 法为将两线段进行矢量跨立。若相交,去除该条路径中的此直线段。 (2)在(1)的基础上去除各条路径中的直线段与各圆弧路径对应圆周及 2 号圆形障 碍物边界相交的部分。判断各条路径中的直线段与圆周是否相关的方法为:判断圆心到 直线段的距离是否小于半径,若距离小于半径,则看垂足是否在线段上,垂足在此线段 上则相交,否则不相交。若圆心到直线段的距离大于半径,则不相交。 (3)依据避障条件,选取各切点在圆弧上的可取点范围,如图所示。

切点只能在 AB 上选取 判断条件为: 切点与圆心连线的向量与过圆心的两边界向量的夹角是否大于等于 90 度,若是,则在允许的取点范围。 基于以上(1) (2) (3)的剔除,剩余的路径即为满足避障条件的可行路径,而避 障最短路径必是这些可行路径中的一条。 Step2:构造扩充的赋权邻接矩阵。将每一个圆弧路径,转化为由两个切点标识,给 此两个切点间赋予边权为对应圆弧路径的弧长。对两切点间的是有直线段相连的,直接 将直线段长度赋予给此两切点间的边权。由此对所有可行路径,构造有扩充点的赋权图 的邻接矩阵。 Setp3:利用 Floyd 算法求扩充赋权图上任意两点间的最短路径,即可转化为任意两 点间的避障最短路径。 5.4.2 避障最短路径简化模型的求解过程与结果 根据构建的启发式算法步骤,利用 MATLAB 软件,我们可求解如下结果: (1)寻找可行路径 根据避障约束条件 1,任意一条可行路径与所有障碍物的边界线段间均不相交,求 得只满足路径线段与边界线段不相交时的可行路径如图 6。

14

图 6 可行路径图 只满足路径中的各直线段与各障碍物边界均不相交的可行路径既满足各直线段与 各障碍物边界均不相交的要求, 又同时满足与各圆弧路径对应的圆周及第 2 号障碍物圆 周均不相交的条件的可行路径如图 7。

图 7 满足三条件的可行路径图 在满足算法 step1 中(1) (2)的基础上,剔除掉切点不在允许范围之内的路径线 段,得出所有可行路径示意图,如图 8。

15

图 8 所有可行路径示意图 (2)对于每条可行路径的直线段和圆弧段,构造扩充点,对任意两点间按要求赋 权,则问题转化为在赋权图上求两点之间的最短距离。 (3)调用 Floyd 算法,先后求出 O-->A 的最短路径及相应最短路长,如图 9.

图 9 O-->A 的最短路径 0-->A
路线序号 1 圆弧起点坐标 70.5063,213.1405 圆弧终点坐标 75.975,219.1542 圆弧圆心坐标 80,210 最短路长 471.0372

同理,我们可得到 O-->B 的最短路径及相应最短路长、O-->C 的最短路径及相应最 短路长、O-->A-->B-->C-->O 的最短路径及相应最短路长.

16

图 10 O-->B 的最短路径 O-->B 的最短路径为: (0,0)-->(50.1354,301.6396)-->(51.6795,305.547)-->(141.6795,440.547)-->(147. 9621,444.79)-->(222.038,460.2096)-->(230,470)-->(230,530)-->(229.9563,530.9 242)-->(229.5746,532.8954)-->(229.2564,533.7791)-->(229.1263,534.0782)-->(2 25.4971,538.3544)-->(144.5036,591.6465)-->(140.6922,596.346)-->(100,700) 对应圆弧的圆心坐标:(60,300),(150,435), (220,470),(220,530),(150,600) 最短路长为:853.7001

图 11 O-->C 的最短路径图 O-->C 的最短路径为: (0,0)-->(232.1147,50.2273)-->(232.1694,50.2377)-->(412.1695,90.2373)-->(418 .3435,94.4905)-->(491.6536,205.5113)-->(492.0569,206.0862)-->(727.9298,513. 924)-->(730,520)-->(730,600)-->(727.7179,606.359)-->(700,640) 经过的圆心:(410,100), (230,60), (720,520), 圆心:(720,600),半径:10 圆
17

心:(500,200),半径:10 最短路长为:1088.1952

图 12 O-->A-->B-->C-->O 的最短路径图 O-->A-->B-->C-->O 的最短路径为: (0,0)-->(70.5063,213.1405)-->(75.975,219.1542)-->(76.2776,219.2811)-->(76.6 064,219.4067)-->(300,300)-->(229.5746,532.8954)-->(229.2564,533.7791)-->(22 9.1263,534.0782)-->(225.4971,538.3544)-->(144.5036,591.6465)-->(140.6922,59 6.346)-->(100,700)-->(270.5862,689.9828)-->(272.0002,689.799)-->(368.0003,6 70.2035)-->(370,670)-->(430,670)-->(434.0793,670.8716)-->(435.5886,671.7056 )-->(534.4132,738.2917)-->(540,740)-->(670,740)-->(679.7741,732.1462)-->(70 0,640)-->(727.7179,606.359)-->(730,600)-->(730,520)-->(727.9298,513.924)--> (492.0569,206.0862)-->(491.6536,205.5113)-->(418.3435,94.4905)-->(412.1695, 90.2373)-->(232.1694,50.2377)-->(232.1147,50.2273)-->(0,0) 经过的圆心: (410,100), (230,60), (80,210), (220,530), (150,600), (270,680), (370,680), (430,680), (670,730), (540,730), (720,520), (720,600), (500,200) 最短路长为:2725.1596。 本次求解,在搜索可行路径上程序运行时间为 10 分钟左右。 5.5 基于数值思想的改进算法 对每条路径,直接根据避障条件,判断它是否可行路径的方法——基于数值近似的 线段间最短距离判别法。 线段间最短距离判别法:给定路径端点坐标,所有障碍物边界线段的坐标,分别在 路径和某一边界线段上取一系列的点, 求得这系列点之间的最短距离即为线段间最短距 离的近似值,如果该近似值大于 10,则认为路径与该边界线段的距离满足要求,否则不 满足。直到该路径与所有的边界线段都满足距离要求,则该路径为可行路径并记录。易 知,若线段间的一系列点取值越细密,结果越精确。 按照该算法,我们编写 MATLAB 程序得到的所有可行路径如下图:

18

对此可行路径集合, 仍然使用基于避障条件转化之下的最短路径的启发算法的 Step2、 Step3 步寻找最短路径,最短路径结果完全等同基于避障条件转化之下的最短路径的启 发算法的结果,但是求解时间却大大缩减,说明该算法具有一定的有效性。

六、最短时间路径模型建立与求解
6.1 到任意目标点的避障最短时间路径模型的建立 基于问题一转化为最短路径的优化问题后, 我们易知优化的目标函数为机器人在行 进过程中最短时间路径,通过 0-1 变量来选取经过的转弯点的圆心坐标,两切点坐标及 转弯半径,此时将转弯点的圆心坐标,两切点坐标及转弯半径均作为决策变量。因此, 可建立 0-1 规划模型如下。 6.1.1 目标函数 目标函数为机器人出区域中一点到达另一点的避障最短时间路径, 避障最短路径 Z 为行进过程中直线路径与转弯路径的时间取和,故有
min Z ? ?? X ij (
i ?1 j ?1 n n

aij v0

?

bij v?

)

(8)

其中, aij 为切点 i 到切点 j 的直线距离;bij 为机器人从切点 i 行进至 j 切点时的转弯 路径; X ij 为 0-1 变量,表示若选择行走切点 ( xi , yi ) 后行走切点 ( x j , y j ) , X ij 为 1,否则 为 0, i, j ? 1...n, n ? 36, i ? j 。机器人从切点 i 到切点 j 的直线路径 aij 为

19

aij ? xi ? yi
2

2

? ( x j ? 300)2 ? ( y j ? 300)2

i, j ? 1..36, i ? j

(9)

机器人从切点 i 行进至 j 切点时的 bij 为
bij ? ??k i, j ? 1..36, i ? j

其中
( xi ? x j ) 2 ? ( yi ? y j ) 2 ? ?i ? ?i ? 2 ?i cos?
2

(10)

?k 为转弯半径。
6.1.2 约束条件 1.避障条件的约束 原理与最短路径的避障条件的原理保持一致。 2.各坐标点的约束 所有假设的坐标点都应在 800 ? 800 平面场景内,于是有
xi1, xi 2 , x j1 , yi1 , xi1 , yi 2 , y j1 , x p , x, y p , y p ? {0,800}

3.各圆弧点行进先后次序的约束 机器人从某个圆弧出发行至下一圆弧,先后次序的约束有
xi1 ? xi 2

yi1 ? yi 2
x j1 ? x j 2
y j2 ? y j2

(11)

xi 2 ? x j1 y i 2 ? x j1

其中, i, j ? 1..36, i ? j 。 4.转弯半径、圆心、与切点的坐标关系。 转弯半径等于圆心 ( x0 , y0 ) 与两个切点 ( x i , yi ) , ( x j , y j ) 坐标的直线距离。

?k ?
?k ?
5. 转弯半径 ? i 取值为:使
min(

( x0 ? xi ) 2 ? ( y0 ? yi ) 2
( x0 ? x j ) 2 ? ( y0 ? y j ) 2

aij v0

?

) v? 满足的 ?

bij

20

6.1.3 避障最短时间路径模型 综合上述,建立避障最短时间路径模型是在:
min Z ? ?? X ij (
i ?1 j ?1 n n

aij v0

?

bij v?

)

? ? ?Yij ? X ij ? ? xi1 ? xi 2 ? yi1 ? yi 2 ? ? y j2 ? y j2 ?x ? x j1 ? i2 ? yi 2 ? x j1 ? ? s.t ?bij ? ??k i, j ? 1..36, i ? j ? 2 2 2 ?( xi ? x j ) ? ( yi ? y j ) ? ? k ? ? k ? 2 ? k cos? ? x , x , x , y , x , y , y , x , x, y , y ? (0,800) p p ? i1 i 2 j1 i1 i1 i 2 j1 p ? 2 2 2 2 ?aij ? xi ? yi ? ( x j ? 300 ) ? ( y j ? 300 ) ? X , Y ? 0或1(i, j ? 1,2,...n, n ? 36, i ? j) ? ij ij ?? ? ( x ? x )2 ? ( y ? y )2 0 i 0 i ? k ? 2 2 ? ? k ? ( x0 ? x j ) ? ( y0 ? y j ) ?

6.2 求从 0-->A 的最短时间路径的简化模型 针对只求 0-->A 的最短时间路径, 我们具体问题具体分析, 建立了相应的简化模型。 障碍物 5 左上角的顶点与圆弧的距离大于等于 10 的避障约束可以表述如下:
d ? ? i ? (80 ? x0 ) 2 ? ( 210 ? y0 ) 2 ? 10

于是,我们可得一条路径时间最短下的最优 ? 的简化模型
min Z ? ?? X ij (
i ?1 j ?1 n n

aij v0

?

bij v?

)

?d ? ? ? (80 ? x ) 2 ? (210 ? y ) 2 ? 10 i 0 0 ? i, j ? 1..36, i ? j ?bij ? ??i ? 2 2 2 ?( xi ? x j ) ? ( yi ? y j ) ? ? i ? ? i ? 2 ? i cos? s.t ? ? xi1, xi 2 , x j1 , yi1 , xi1 , yi 2 , y j1 , x p , x, y p , y p ? (0,800) ? 2 2 2 2 i, j ? 1..36, i ? j ?aij ? xi ? x j ? ( x j ? 300 ) ? ( y j ? 300 ) ? ? X ij ? 0或1(i, j ? 1,2,...n, n ? 36, i ? j)

21

针对每一条可行路径,求出各路径对应的最优转弯半径 ?i (i ? 1,2) 则,最短时间路径的优化模型为:

?n n aij bij n n aij bij ? ? ? Z ? min ??? X ij ( ? ), ?? X ij ( ? )? v0 v?1 i ?1 j ?1 v0 v?2 ? ? i ?1 j ?1 ? ?
6.3 求从 0-->A 的最短时间路径的模型求解 根据 0-->A 的最短时间路径的模型,利用 LINGO 软件编程,比较两条可行路径对应 的最短时间, 我们求得最短时间路径下转弯半径为 12.9885 , 同时最短时间路径时间长 为 94.2283 个单位。相应圆弧的圆心坐标为(82.1414,207.9153),两切点坐标分别 为(69.8045,211.9779)、(77.7492,220.1387)。

七、模型评价
7.1 模型优点 (1)确定路线思路循序渐进,本文先建立了有计算避障约束公式的普适性模型, 再建立了以不取相交点来简化 0-1 变量取值关系的简化模型; (2)给出了二种启发式算法,最短路径即最短时间路径具有一定可信度。同时第 一个启发算法可以求得全局最优解,第二个启发算法是针对问题的 NP 属性减少求解时 间而构建的,两个算法都具有较重要的意义。 7.2 模型缺点 (1)本文将机器人看成一个质点,这将使机器人出现走区域边界的可能,可能会 出现与实际不符合的情况; (2)模型将假设机器人在行进、转弯过程中能一直保持最大的速度,诚然,现实 并非如此,所以我们得到的问题二的结果与实际最短时间路径会存在一定的误差。

参考文献
[1] 吴建国.数学建模案例精编[M], 北京:中国水电出版社,2006,210 [2] 杨秀月等.系统建模[M],北京:国防工业出版社,2006.5 [3] 张志涌等.精通 MATLAB6.5 版[M].北京: 北京航空航天大学出版社,2003,3

22

附录

题目平面图 MATLAB 程序清单:clear distancebetweenlines.m
Ax=A(1);Ay=A(2); Bx=B(1);By=B(2); Cx=C(1);Cy=C(2); Dx=D(1);Dy=D(2); if (A(1)-B(1))~=0 k=(By-Ay)/(Bx-Ax); b=By-k*Bx; ABXXL=linspace(A(1),B(1),M); ABYXL=k.*ABXXL+b; else ABXXL=linspace(A(1),B(1),M);
23

…………………………计算两线段间近似最短路径

function mind=distancebetweenlines(A,B,C,D,M)

ABYXL=linspace(A(2),B(2),M); end if (C(1)-D(1))~=0 k=(Dy-Cy)/(Dx-Cx); b=Dy-k*Dx; CDXXL=linspace(C(1),D(1),M); CDYXL=k.*CDXXL+b; else CDXXL=linspace(C(1),D(1),M); CDYXL=linspace(C(2),D(2),M); end mind=100000; for i=1:M for j=1:M if sqrt((ABXXL(i)-CDXXL(j))^2+(ABYXL(i)-CDYXL(j))^2)<=mind mind=sqrt((ABXXL(i)-CDXXL(j))^2+(ABYXL(i)-CDYXL(j))^2); end end end

FLOYD1.m
n=size(w,1); D=w; path=zeros(n,n);

?????????????Folyd 算法求最短路径和最短路程

function [L,R]=FLOYD1(w,s,t)

%??????±ê×?floyd??·¨ for i=1:n for j=1:n if D(i,j)~=inf path(i,j)=j; end end end for k=1:n for i=1:n for j=1:n if D(i,k)+D(k,j)<D(i,j) D(i,j)=D(i,k)+D(k,j); path(i,j)=path(i,k); end end end end L=zeros(0,0); R=s;
24

while 1 if s==t L=fliplr(L); L=[0,L]; return end L=[L,D(s,t)]; R=[R,path(s,t)]; s=path(s,t); end

huchang.m
%C??????×?±ê

?????????????已知圆心和圆上 2 点,计算弧长

function z=huchang(A,B,C,r) %A,B??????????×?±ê,r??°??? D=dot(A-C,B-C)/(norm(A-C)*norm(B-C));%???????ò?????? theta=acos(D);%????±í?? z=theta*r;

iscycleIntersect.m
if (line(3)-line(1))~=0

?????????????判断线段是否和圆相交

function z=iscycleIntersect(line,A,r) k=(line(4)-line(2))/(line(3)-line(1)); b=line(2)-k*line(1); d=abs(k*A(1)-A(2)+b)/sqrt(k^2+1); [x y]=solve(['x*(' num2str(k) ')-y+(' num2str(b) ')=0'],... ['(y-(' num2str(A(2)) '))*(' num2str(k) ')=(' num2str(A(1)) ')-x']); else x=line(3);y=A(2); d=abs(line(3)-A(1)); end x=double(x);y=double(y); if d>=r-0.5 z=0; elseif (y-line(4))*(y-line(2))>0 z=0; else z=1; end

islineIntersect.m
Ax=A(1);Ay=A(2); Bx=B(1);By=B(2); Cx=C(1);Cy=C(2); Dx=D(1);Dy=D(2);

?????????????判断两线段是否相交

function z=islineIntersect(A,B,C,D)

25

if ((Bx-Ax)*(Dy-Cy)-(By-Ay)*(Dx-Cx))*((Bx-Ax)*(Dy-Cy)-(By-Ay)*(Dx-Cx))~=0 r=((Ay-Cy)*(Dx-Cx)-(Ax-Cx)*(Dy-Cy))/((Bx-Ax)*(Dy-Cy)-(By-Ay)*(Dx-Cx)); s=((Ay-Cy)*(Bx-Ax)-(Ax-Cx)*(By-Ay))/((Bx-Ax)*(Dy-Cy)-(By-Ay)*(Dx-Cx)); if r>0&&r<=1&&s>0&&s<=1 z=1; else z=0; end else z=0; end

linedistance.m
a=B(1:2)-B(3:4); b=B(5:6)-B(3:4); c=A-B(3:4);

?????????????判断切点是否位于可行圆弧上

function z=linedistance(A,B)

theta1=acos(dot(a,c)/(norm(a)*norm(c)))*180/pi; theta2=acos(dot(b,c)/(norm(b)*norm(c)))*180/pi; if theta1>=90-0.5&&theta2>=90-0.5 z=0; else z=1; end

method1.m
clear clc close all theta=0:pi/100:2*pi;

?????????????启发式算法计算所有可行路径

zb{1}=[300 400;500 400;500 600;300 600]; zb{2}=[550 450;80 80]; zb{3}=[360 240;500 240;540 330;400 330;]; zb{4}=[280 100;410 100;345 210;]; zb{5}=[80 60;230 60;230 210;80 210]; zb{6}=[60 300;235 300;150 435;]; zb{7}=[0 470;220 470;220 530;0 530]; zb{8}=[150 600;240 600;270 680;180 680]; zb{9}=[370 680;430 680;430 800;370 800]; zb{10}=[540 600;670 600;670 730;540 730]; zb{11}=[640 520;720 520;720 600;640 600]; zb{12}=[500 140;800 140;800 200;500 200]; zb{13}=[0 0;800 0;800 800;0 800;]; zb{14}=[0 0 79]; zb{15}=[300 300 65]; zb{16}=[100 700 66];
26

zb{17}=[700 640 67]; %***********************???????ò???ù?????±??????·???*********************** ** lines=[];kxyh=[]; for i=1:length(zb) temp=zb{i}; if size(zb{i},1)==2 x=temp(2,1)*cos(theta)+temp(1,1); y=temp(2,2)*sin(theta)+temp(1,2); plot(x,y,'b-');hold on elseif size(temp,1)==3 temp0=temp+[-10 -10;10 -10;0 10]; for j=1:size(temp,1) if j==size(temp,1) end plot([temp(:,1);temp(1,1)],[temp(:,2);temp(1,2)],'b-');hold on elseif size(temp,1)==4 if i==13 temp0=temp+[10 10;-10 10;-10 -10;10 -10]; else temp0=temp+[-10 -10;10 -10;10 10;-10 10]; end for j=1:size(temp,1) if j==size(temp,1) tt=1; else tt=j+1;end lines=[lines [temp(j,1);temp(j,2);temp(tt,1);temp(tt,2)]]; end plot([temp(:,1);temp(1,1)],[temp(:,2);temp(1,2)],'b-');hold on elseif size(temp,1)==1 plot(temp(1),temp(2),'r.','MarkerSize',12);text(temp(1)+10,temp(2)+20,char( temp(3)));hold on end if i<=12&i~=2 for j=1:size(temp,1) temp1=temp(j,:); x=10*cos(theta)+temp1(1); y=10*sin(theta)+temp1(2); plot(x,y,'b-');hold on end end end axis([0 800 0 800]);grid on %***********************???????ò???ù????????????·???*********************** **
27

tt=1;

else

tt=j+1; end

lines=[lines [temp(j,1);temp(j,2);temp(tt,1);temp(tt,2)]];

kkk=1; for i=1:length(zb)-4 temp5=zb{i};if i==2 pp=1;else pp=size(temp5);end for j=1:pp if ~ismember(temp5(j,1),[0,800]) & ~ismember(temp5(j,2),[0,800]) if j==1 ppp=size(temp5,1); else ppp=j-1; end if j==size(temp5,1) tt=1; else tt=j+1;end kxyh(:,kkk)=[temp5(ppp,1);temp5(ppp,2);temp5(j,1);temp5(j,2);temp5(tt,1);te mp5(tt,2)]; yxzb(1:2,kkk)=temp5(j,:); if i==2 yxzb(3,kkk)=80;else yxzb(3,kkk)=10;end kkk=kkk+1; end end end % ******************??????????????O,A,B,C???÷?????????????·??**************** yxbh=0; for i=1:4 %**********************?????????????÷?°??????×?±ê*********************** * x=[];y=[];z=[];kk=1; temp=zb{13+i};A=temp(1:2); for j=1:12 temp2=zb{j}; if j==2 r=70;temp2=temp2(1,:); else r=10; end for k=1:size(temp2,1) B=temp2(k,:); if ~ismember(B(1),[0,800]) & ~ismember(B(2),[0,800]) [x(:,kk) y(:,kk)]=qiedian(A,B,r); z=[z kk]; kk=kk+1; end end end qdzb=[reshape(x,1,2*length(x));reshape(y,1,2*length(y));reshape([z;z],1,2*l ength(z))]; %*****%********************??????·????????·??*************************** for p=1:length(qdzb)
28

%*********************???????±???????à??**************************** for q=1:length(lines) temp3(p,q)=islineIntersect(temp,qdzb(1:2,p),lines(1:2,q),lines(3:4,q)); if temp3(p,q)>0 break;end end %*********************???????????????à??**************************** for q=(length(lines)+1):(length(lines)+length(yxzb)) temp3(p,q)=iscycleIntersect([A';qdzb(1:2,p)],yxzb(1:2,q-length(lines)),yxzb (3,q-length(lines))); if temp3(p,q)>0 break;end end %****************??????????·???????????????************************* if linedistance(qdzb(1:2,p),kxyh(:,qdzb(3,p)))==0 temp3(p,(length(lines)+length(yxzb)+1):(length(lines)+length(yxzb)+length(k xyh)))=0; else temp3(p,(length(lines)+length(yxzb)+1):(length(lines)+length(yxzb)+length(k xyh)))=1; end end ind{i}=find(sum(temp3,2)==0);temp4=ind{i}; for p=1:length(temp4) plot([A(1) qdzb(1,temp4(p))],[A(2) qdzb(2,temp4(p))],'r-');hold on end dykxlj{i,1}=[repmat(A',1,length(ind{i}));[qdzb(1,ind{i});qdzb(2,ind{i})];re pmat([length(yxzb)+i],1,length(ind{i}));qdzb(3,ind{i})]; end %%%******************?????????÷???????????????????·??********************** * yykxlj=[]; for i=1:length(yxzb) for j=i+1:length(yxzb) temp6=yuanyuanqiexian(yxzb(1:2,i),yxzb(3,i),yxzb(1:2,j),yxzb(3,j)); temp7=[]; for k=1:size(temp6,2) %********************???????±???????à??************************** * for L=1:length(lines) temp7(k,L)=islineIntersect(temp6(1:2,k),temp6(3:4,k),lines(1:2,L),lines(3:4
29

,L)); end %********************???????????????à??************************** ** for L=(length(lines)+1):(length(lines)+length(yxzb)) temp7(k,L)=iscycleIntersect(temp6(:,k),yxzb(1:2,L-length(lines)),yxzb(3,L-l ength(lines))); end %****************??????????·???????????????********************** *** if linedistance(temp6(1:2,k),kxyh(:,i))==0&&linedistance(temp6(3:4,k),kxyh(:,j ))==0 temp7(k,(length(lines)+length(yxzb)+1):(length(lines)+length(yxzb)+length(k xyh)))=0; else temp7(k,(length(lines)+length(yxzb)+1):(length(lines)+length(yxzb)+length(k xyh)))=1; end end ind2=find(sum(temp7,2)==0); if ~isempty(ind2) yykxlj=[yykxlj [temp6(:,ind2);repmat([i;j],1,length(ind2))]]; end end end for i=1:length(yykxlj) plot(yykxlj([1 3],i),yykxlj([2 4],i),'k-');hold on end save D11

method2.m
clear clc close all

?????????????线段距离判别法计算所有可行路径

theta=0:pi/100:2*pi; zb{1}=[300 400;500 400;500 600;300 600]; zb{2}=[550 450;80 80]; zb{3}=[360 240;500 240;540 330;400 330;]; zb{4}=[280 100;410 100;345 210;]; zb{5}=[80 60;230 60;230 210;80 210]; zb{6}=[60 300;235 300;150 435;];
30

zb{7}=[0 470;220 470;220 530;0 530]; zb{8}=[150 600;240 600;270 680;180 680]; zb{9}=[370 680;430 680;430 800;370 800]; zb{10}=[540 600;670 600;670 730;540 730]; zb{11}=[640 520;720 520;720 600;640 600]; zb{12}=[500 140;800 140;800 200;500 200]; zb{13}=[0 0 79]; zb{14}=[0 0 79]; zb{15}=[300 300 65]; zb{16}=[100 700 66]; zb{17}=[700 640 67]; %***********************???????ò???ù?????±??????·???*********************** ** lines=[];kxyh=[]; for i=1:length(zb) temp=zb{i}; if size(zb{i},1)==2 x=temp(2,1)*cos(theta)+temp(1,1); y=temp(2,2)*sin(theta)+temp(1,2); plot(x,y,'b-');hold on elseif size(temp,1)==3 temp0=temp+[-10 -10;10 -10;0 10]; for j=1:size(temp,1) if j==size(temp,1) end plot([temp(:,1);temp(1,1)],[temp(:,2);temp(1,2)],'b-');hold on elseif size(temp,1)==4 if i==13 temp0=temp+[10 10;-10 10;-10 -10;10 -10]; else temp0=temp+[-10 -10;10 -10;10 10;-10 10]; end for j=1:size(temp,1) if j==size(temp,1) tt=1; else tt=j+1;end lines=[lines [temp(j,1);temp(j,2);temp(tt,1);temp(tt,2)]]; end plot([temp(:,1);temp(1,1)],[temp(:,2);temp(1,2)],'b-');hold on elseif size(temp,1)==1 plot(temp(1),temp(2),'r.','MarkerSize',12);text(temp(1)+10,temp(2)+20,char( temp(3)));hold on end if i<=12&i~=2 for j=1:size(temp,1)
31

tt=1;

else

tt=j+1; end

lines=[lines [temp(j,1);temp(j,2);temp(tt,1);temp(tt,2)]];

temp1=temp(j,:); x=10*cos(theta)+temp1(1); y=10*sin(theta)+temp1(2); plot(x,y,'b-');hold on end end end axis([0 800 0 800]);grid on %***********************???????ò???ù????????????·???*********************** ** kkk=1; for i=1:length(zb)-4 temp5=zb{i};if i==2 pp=1;else pp=size(temp5);end for j=1:pp if ~ismember(temp5(j,1),[0,800]) & ~ismember(temp5(j,2),[0,800]) if j==1 ppp=size(temp5,1); else ppp=j-1; end if j==size(temp5,1) tt=1; else tt=j+1;end kxyh(:,kkk)=[temp5(ppp,1);temp5(ppp,2);temp5(j,1);temp5(j,2);temp5(tt,1);te mp5(tt,2)]; yxzb(1:2,kkk)=temp5(j,:); if i==2 yxzb(3,kkk)=80;else yxzb(3,kkk)=10;end kkk=kkk+1; end end end % ******************??????????????O,A,B,C???÷?????????????·??**************** yxbh=0; for i=1:4 %**********************?????????????÷?°??????×?±ê*********************** * x=[];y=[];z=[];kk=1; temp=zb{13+i};A=temp(1:2); for j=1:12 temp2=zb{j}; if j==2 r=70;temp2=temp2(1,:); else r=10; end for k=1:size(temp2,1) B=temp2(k,:); if ~ismember(B(1),[0,800]) & ~ismember(B(2),[0,800]) [x(:,kk) y(:,kk)]=qiedian(A,B,r);
32

z=[z kk]; kk=kk+1; end end end qdzb=[reshape(x,1,2*length(x));reshape(y,1,2*length(y));reshape([z;z],1,2*l ength(z))]; %*****%********************??????·????????·??*************************** temp3=ones(length(qdzb),length(lines)); for p=1:length(qdzb) %******************???????±???????à??????10????********************* ***** for q=1:length(lines) if distancebetweenlines(A,qdzb(1:2,p),lines(1:2,q),lines(3:4,q),100)>=10-0.01 temp3(p,q)=0; else temp3(p,q)=1; end end end ind{i}=find(sum(temp3,2)==0);temp4=ind{i}; for p=1:length(temp4) plot([A(1) qdzb(1,temp4(p))],[A(2) qdzb(2,temp4(p))],'r-');hold on end dykxlj{i,1}=[repmat(A',1,length(ind{i}));[qdzb(1,ind{i});qdzb(2,ind{i})];re pmat([length(yxzb)+i],1,length(ind{i}));qdzb(3,ind{i})]; end %%%******************?????????÷???????????????????·??********************** * yykxlj=[]; for i=1:length(yxzb) for j=i+1:length(yxzb) temp6=yuanyuanqiexian(yxzb(1:2,i),yxzb(3,i),yxzb(1:2,j),yxzb(3,j)); temp7=[]; for k=1:size(temp6,2) %******************???????±???????à??????10????******************* ***** for L=1:length(lines) temp7(k,L)=distancebetweenlines(temp6(1:2,k),temp6(3:4,k),lines(1:2,L),line s(3:4,L),100); if
33

distancebetweenlines(temp6(1:2,k),temp6(3:4,k),lines(1:2,L),lines(3:4,L),10 0)>=10-.01 temp7(k,L)=0; else temp7(k,L)=1; end end end ind2=find(sum(temp7,2)==0); if ~isempty(ind2) yykxlj=[yykxlj [temp6(:,ind2);repmat([i;j],1,length(ind2))]]; end end end for i=1:length(yykxlj) plot(yykxlj([1 3],i),yykxlj([2 4],i),'k-');hold on end

qiedian.m

?????????????计算圆外一点与圆的切点坐标

function [x,y]=qiedian(A,B,r) [b k]=solve(['(k*' num2str(B(1)) '-' num2str(B(2)) '+b)^2=' num2str(r) '^2*(1+k^2)'],... ['(k*' num2str(A(1)) '-' num2str(A(2)) '+b)=0']); ['(k*' num2str(B(1)) '-' num2str(B(2)) '+b)^2=' num2str(r) '^2*(1+k^2)']; ['(k*' num2str(A(1)) '-' num2str(A(2)) '+b)=0']; k=double(k); b=double(b); for i=1:length(k) ['x*' num2str(k(i)) '-y+' num2str(b(i)) '=0']; ['(y-' num2str(B(2)) ')*' num2str(k(i)) '=' num2str(B(1)) '-x']; [x(i) y(i)]=solve(['x*(' num2str(k(i)) ')-y+(' num2str(b(i)) ')=0'],... ['(y-(' num2str(B(2)) '))*(' num2str(k(i)) ')=(' num2str(B(1)) ')-x']); end x=double(x);y=double(y);

xslj.m

?????????????显示最短路径

function xslj(yxzb,sykxd,ROA) OA=['(' num2str(sykxd(1,ROA(1))) ',' num2str(sykxd(1,ROA(1))) ')']; for i=2:length(ROA) OA=[OA ['-->(' num2str(sykxd(1,ROA(i))) ',' num2str(sykxd(2,ROA(i))) ')']]; B(i)=sykxd(3,ROA(i)); end BB=unique(B); s=hist(B,BB); t=yxzb(:,BB(find(s>=2))); YX=[];
34

for i=1:size(t,2) YX=[YX '??????(' num2str(t(1,i)) ',' num2str(t(2,i)) ')??°?????' num2str(t(3,i)) ' end disp([OA]); disp(['?-??????????°?????' YX]); '];

yuanyuanqiexian.m

?????????????计算圆和圆的四条切线和切点

function [x]=yuanyuanqiexian(A,r1,B,r2) [b k]=solve(['(k*' num2str(A(1)) '-' num2str(A(2)) '+b)^2=' num2str(r1) '^2*(1+k^2)'],... ['(k*' num2str(B(1)) '-' num2str(B(2)) '+b)^2=' num2str(r2) '^2*(1+k^2)']); k=double(real(k)); b=double(real(b)); for i=1:length(k) temp5=solve(['x*(' num2str(k(i)) ')-y+(' num2str(b(i)) ')=0'],... ['(y-(' num2str(A(2)) '))*(' num2str(k(i)) ')=(' num2str(A(1)) ')-x']); x(1,i)=double(temp5.x);x(2,i)=double(temp5.y); temp5=solve(['x*(' num2str(k(i)) ')-y+(' num2str(b(i)) ')=0'],... ['(y-(' num2str(B(2)) '))*(' num2str(k(i)) ')=(' num2str(B(1)) ')-x']); x(3,i)=double(temp5.x);x(4,i)=double(temp5.y); end x=double(x); if length(k)==2 x(1,3)=A(1)-r1;x(2,3)=A(2);x(3,3)=B(1)-r1;x(4,3)=B(2); x(1,4)=A(1)+r1;x(2,4)=A(2);x(3,4)=B(1)+r1;x(4,4)=B(2); end

zdlj.m
clear clc load D11

?????????????构造距离矩阵, 并调用 Folyd 算法求最短路径

%*****************************?????ù?????????·???°??????????*************** *** sykxlj=[yykxlj]; for i=1:size(dykxlj,1) sykxlj=[sykxlj dykxlj{i}]; end ddzb=[yxzb [0;0;10] [300;300;10] [100;700;10] [700;640;10]]; %????×?±ê??O??A??B??C sykxd0=[sykxlj([1:2 5],:) sykxlj([3:4 6],:)]; sykxd=unique(sykxd0','rows');sykxd=sykxd'; %*********************Folyd??·¨????OA??OB??OC??×????·?????·??************** **** dis=ones(length(sykxd),length(sykxd))*inf; for i=1:length(sykxlj) C=find((ismember(sykxd',(sykxlj([1:2 5],i))','rows'))');
35

D=find((ismember(sykxd',(sykxlj([3:4 6],i))','rows'))'); dis(C,D)=sqrt((sykxlj(1,i)-sykxlj(3,i))^2+(sykxlj(2,i)-sykxlj(4,i))^2); dis(D,C)=dis(C,D); end for i=1:length(sykxd) for j=i+1:length(sykxd) if sykxd(3,i)==sykxd(3,j) dis(i,j)=huchang(sykxd(1:2,i),sykxd(1:2,j),ddzb(1:2,sykxd(3,i)),ddzb(3,sykx d(3,i))); end dis(j,i)=dis(i,j); end end dis=real(dis); [LOA,ROA]=FLOYD1(dis,1,107); [LOB,ROB]=FLOYD1(dis,1,20); [LOC,ROC]=FLOYD1(dis,1,198); [LAB,RAB]=FLOYD1(dis,107,20); [LBC,RBC]=FLOYD1(dis,20,198); disp('O-->A??×????·??????') xslj(yxzb,sykxd,ROA) disp(['×??¤??????' num2str(LOA(end))]) disp('O-->B??×????·??????') xslj(yxzb,sykxd,ROB) disp(['×??¤??????' num2str(LOB(end))]) disp('O-->C??×????·??????') xslj(yxzb,sykxd,ROC) disp(['×??¤??????' num2str(LOC(end))]) disp('O-->A-->B-->C-->O??×????·??????') xslj(yxzb,sykxd,[ROA(1:end-1) RAB(1:end-1) RBC(1:end-1) rot90(ROC,2)]) disp(['×??¤??????' num2str(LOA(end)+LAB(end)+LBC(end)+LOC(end))]) hold on plot(sykxd(1,ROA),sykxd(2,ROA),'b-') ; plot(sykxd(1,ROB),sykxd(2,ROB),'r-') ; plot(sykxd(1,ROC),sykxd(2,ROC),'g-') ; plot(sykxd(1,[ROA RAB RBC]),sykxd(2,[ROA RAB RBC]),'k-') ;hold off

LINGO 程序(请用 lingo9.0 以后版本) LINGO1.lg4 ?????????????问题 2 优化模型 1 障碍 5 左上角转弯) ( data: x1=0;y1=0; x2=300;y2=300;
36

dx=80;dy=210; v0=5; enddata min=OD+EA+DE; !满足(x3,y3)和(x4,y4)为切点; (x3-x1)*(y3-y0)-(y3-y1)*(x3-x0)=0; (x4-x2)*(y4-y0)-(y4-y2)*(x4-x0)=0; !半径减去圆心到顶点的距离必须大于等于 10 个单位, 即保证顶点与圆弧的最短距离大 于等于 10; d=((dx-x0)^2+(dy-y0)^2)^(1/2); r=((x3-x0)^2+(y3-y0)^2)^(1/2); r-d>=10; !切点到边界线的距离必须大于 10; x3<=dx; y3>=dy; x4<=dx; y4>=dy; !直线段所需时间; OD=((x3-x1)^2+(y3-y1)^2)^(1/2)/v0; EA=((x4-x2)^2+(y4-y2)^2)^(1/2)/v0; !圆弧段所需时间; v=v0/(1+@exp(10-0.1*r^2)); theta=@acos(((x4-x0)*(x3-x0)+(y4-y0)*(y3-y0))/r); DE=theta*r/v; LINGO2.lg4 ?????????????问题 2 优化模型 1 障碍 5 右下角转弯) ( model: data: x1=0;y1=0; x2=300;y2=300; dx=230;dy=60; v0=5; enddata min=OD+EA+DE; !满足(x3,y3)和(x4,y4)为切点; (x3-x1)*(x3-x0)+(y3-y0)*(y3-y1)=0; (x4-x2)*(x4-x0)+(y4-y2)*(y4-y0)=0;

37

!半径减去圆心到顶点的距离必须大于等于 10 个单位, 即保证顶点与圆弧的最短距离大 于等于 10; d=((dx-x0)^2+(dy-y0)^2)^(1/2); r=((x3-x0)^2+(y3-y0)^2)^(1/2); r=((x4-x0)^2+(y4-y0)^2)^(1/2); r-d>=10; !切点到边界线的距离必须大于 10; x3>=dx; y3<=dy; x4>=dx; y4<=dy; !直线段所需时间; OD=((x3-x1)^2+(y3-y1)^2)^(1/2)/v0; EA=((x4-x2)^2+(y4-y2)^2)^(1/2)/v0; !圆弧段所需时间; v=v0/(1+@exp(10-0.1*r^2)); theta=@acos(((x4-x0)*(x3-x0)+(y4-y0)*(y3-y0))/r^2); DE=theta*r/v;

38


相关文章:
...葡萄酒的评价2012高教社杯全国大学生数学建模竞赛
2012全国大学生数学建模全国一等奖 优秀论文 葡萄酒的评价2012高教社杯全国大学生数学建模竞赛_数学_自然科学_专业资料。采用方差检验、灰色关联、数据样本统计分析、...
2014高教社杯全国大学生数学建模竞赛D题获奖论文_图文
2014高教社杯全国大学生数学建模竞赛D题获奖论文_数学_自然科学_专业资料。2014高教社杯全国大学生数学建模竞赛获奖论文,为广大学子和数学建模人士提供依据和参考。...
2012高教社杯全国大学生数学建模竞赛B题获奖论文
2012高教社杯全国大学生数学建模竞赛B题获奖论文_其它_高等教育_教育专区。写的很好,注意格式今日推荐 50份文档 2014年注册会计师考试 ...
2012高教社杯全国大学生数学建模竞赛一等奖论文
2012 高教社杯全国大学生数学建模竞赛 承 诺 书 我们仔细阅读了中国大学生数学...我们参赛选择的题号是(从 A/B/C/D 中选择一项填写) : 我们的参赛报名号为...
2012高教社杯全国大学生数学建模竞赛A题论文
2012高教社杯全国大学生数学建模竞赛A题论文_理学_高等教育_教育专区。2012 高教社杯全国大学生数学建模竞赛 承 诺 书 我们仔细阅读了中国大学生数学建模竞赛的...
2013高教社杯全国大学生数学建模竞赛-B题论文_图文
2013高教社杯全国大学生数学建模竞赛-B题论文_教育学_高等教育_教育专区。碎...张艳. 图像拼接技术在文档图像扭曲识别中的应用与研究[D]. 北京:北方工业大学...
2012高教社杯全国大学生数学建模竞赛A题论文最新的
2012高教社杯全国大学生数学建模竞赛A题论文最新的_韩语学习_外语学习_教育专区...) 我们参赛选择的题号是(从 A/B/C/D 中选择一项填写) : A 我们的参赛报名...
2014高教社杯全国大学生数学建模竞赛题目 D
2014 高教社杯全国大学生数学建模竞赛题目(请先阅读“全国大学生数学建模竞赛论文格式规范”) D题 储药柜的设计 储药柜的结构类似于书橱,通常由若干个横向隔板和...
2012高教社杯全国大学生数学建模竞赛_图文
2012高教社杯全国大学生数学建模竞赛 - 2012 高教社杯全国大学生数学建模竞赛 联合赛区评奖结果(初稿) 全国大学生数学建模竞赛组委会,2012 年 10 月 27 日注: ...
2012年“高教社杯”全国大学生数学建模竞赛(CUMCM)国家...
2012年“高教社杯全国大学生数学建模竞赛(CUMCM)国家一等奖优秀论文C题目_数学...我们参赛选择的题号是(从 A/B/C/D 中选择一项填写) : C 我们的参赛报名...
更多相关标签: