数学建模机器人避障

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

承诺书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则。

我们完全明白，在竞赛开始后参赛队员不能以任何形式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导老师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的，如果引用别人的成果获其他公开的资料（包括网上查询到的自资料），必须俺不找规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严重遵守竞赛规则，以保证竞赛的公正性、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。

我们授权全国大学生数学建模竞赛组委会，可将我们的论文一任何形式公开（包括进行网上公示，在书籍，期刊和其他媒体进行正式或非正式发表等）。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）：D题

我们参赛报名号为（如果赛区设置报名号的话）：13277001

所属学校（请填写完整的全名）：长江职业学院

参赛队员（打印并签名）：1. 何宝宝

2. 舒鑫鑫

3. 黎成硕

指导老师或知道教师组：陶燕芳胡芬侯谦民

日期：2012年9 月7 日赛区评阅编号（有赛区组委会评阅前进行编号）：

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

编号专用页

赛区评阅编号（由赛区组委会评阅前进行编号）：

赛区评阅记录（可供赛区评阅时使用）：

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）：

机器人避障问题的研究

摘要

本文主要研究了机器人在一定的区域中，且存在多个障碍物的情况下，由出发点到达目标点，以及由出发点经过途中的多个目标点最后回到出发点的两种情况。由于最短路径一定是由圆弧和与之相切的直线组成，因此我们将路径划分为若干个这种线圆结构来求解。

问题一：

对存在障碍物的两点之间，利用平面几何原理建立最短路程的模型，并用Matlab计算出两点之间可能存在的所有最短距离，然后整理这些距离，最后利用图论中的最短路的算法求得最优解。其次通过对几何原理的分析，同时结合CAD软件计算出机器人从起点出发经过多个目标点最后回到起点的最短路径。结果如下：O→A 最短路径为：471.0372

O→B 最短路径为：853.1121

O→C 最短路径为：1130.281

O→A→B→C→O 最短路径为：1884.149

问题二：

对于问题二，我们讨论了较简单的一种情形，即当圆心固定，半径变化时。在这种情形下，我们利用解析几何法建立了时间与半径之间的函数关系，并使用两种方法对函数进行求解。（1）直接利用Matlab软件绘出了函数图形，并求出了最优解，即当ρ=11.5038时，时间取到最小值T=94.5649。（2）尝试用鱼群算法求出了函数的最优解，最优解为ρ=11.5042T=94.9106，由比较可知两种算法的结果非常接近。

关键词：最短路径、穷举法、解析几何、最优化模型、鱼群算法、超越函数

1

一、问题重述

表1是一个800×800的平面场景图，在原点O(0, 0)点处有一个机器人，它只能在该平面场景范围内活动。图中有12个不同形状的区域是机器人不能与之发生碰撞的障碍物，障碍物的数学描述如下表。

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

机器人直线行走的最大速度为=5个单位/秒。机器人转弯时，最大转弯速度为=(ρ)=

，其中ρ是转弯半径。如果超过该速度，机器人将发生侧翻， .

无法完成行走。

请建立机器人从区域中一点到达另一点的避障最短路径和最短时间路径的数学模型。对场景图中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

二、问题的分析

问题一：要求求定点O（0,0）按照一定的行走规则绕过障碍物到达目标点的最短路径，我们先可以画出机器人行走的危险区域，这样的话，拐角处就是一个半径为10的圆弧。通过采用模拟作图的方法找出可能的最短路径（比如求O和A之间的最短路径，我们可以连接O和A之间的一段绳子，以拐角处的圆弧为支撑点拉紧，那么这段绳子的长度便是O到A的一条可能的最短路径），然后采用穷举法列出O到每个目标点的可能的最短路径。可以将原图中的几何图形单独拿出来转化为线圆模型，将机器人转弯圆弧以最小半径行走，并将直线距离转化为切点与已知点之间的距离。要求出最短路径需要求出起始点，拐点圆心，拐点圆心之间的中点和终点。根据CAD画图确定各个切点的坐标，再根据Matlab编程算出距离总长。最后比较其大小便可得出O到目标点的最短路径。

问题二：此问要求在时间最短的情况下求从出发点到达A点的路径。此问我们讨论了较简单的情形，即当圆心固定在障碍物5的顶点上，而拐弯半径发生变动的情形。通过问题一的分析，我们将圆心选定在正方形的左上顶点(80,210)处。通过平面几何原理，我们建立了时间与半径的函数关系。随后，利用Matlab直接对函数进行作图分析并求出了最优解；同时，我们又选用鱼群算法对函数进行迭代，验证了最优解。

三、模型的假设

（1）假设机器人可以抽象成一个点。

（2）假设机器人内部的控制系统没有出什么问题，能正常工作。

（3）假设假设机器人能迅速反应加速到最大速度和减速到最小速度，蓄势时间可以忽略不计。

（4）假设行走路线不受外界因素干扰，均为理想状态。

四、符号的说明

L 路径的总长度

a i圆心的横坐标

b j圆心的纵坐标

ρ转弯半径

T 最小总时间

五、理论准备

5.1理论与证明

猜想：具有圆形限定区域的最短路径是由两部分组成的：一部分是平面上的自然最短路径（即直线段），另一部分是限定区域的部分边界，这两部分是相切的，互相连接的。

证明：假设在平面中有A（a，0）和B（?a，0）两点，中间有一个半圆形的障碍物，证明从A到B的最路径为A EF?B。

3

4

图1

平面上连接两点最短的路径是通过这两点的直线段，但是连接两点的线段于障碍物相交，所以设法尝试折线路径。在y 轴上取一点C （0，y ），若y 适当大，则折线ACB 与障碍物不相交，折线ACB 的长度为：

|ACB |=2√a 2+y 2

显然|ACB |随着y 的减小而减小，减小y →y 得，即C →C ,使得AC 与C B 与障物相

切，切点分别为E 和F,显然AC B 是这种折线路径中最短的。由于满足0<α<π

2角

满足α

?

下面在考察一条不穿过障碍物的任何一条路径，设其分别于OE 和OF 的延长线交

与P 、Q 两点，记A 和P 之间的路径长度为AP ?，显然AP

?>|AP|,又由AE EO,所以|AP|>AE,从而AP

?>AE,同理可得BQ ?>BF 。 再来比较PQ 之间路径长度PQ

?和圆弧EF 的长度的大小。若PQ 之间的路径可有极坐标方程ρ=ρ(θ),则有ρ>10,可得：

PQ

?=∫√ρ2+ρ2dθ≥∫dθ?π3

?EF ? 结论：路径APQB 的长度超过路径AEFB 的长度，所以满足从A 到B 的最短路径为AEFB 的路径。 5.2模型的准备

在上面的证明中，我们就可以这样认为，起点到目标点无论中间障碍物有多少，最短路径都应该是若干个线圆组成的结构。在本题中存在障碍物的状况，且障碍物在拐点处的危险区域是一个半径为10的圆弧。所以结合上述证明，我们得知，

求两

5

点之间的最短路径中的转弯半径我们应该按照最小的转弯半径来算才能达到最优。

图2

如图2是将题目中原图的一个模型转化成上图所示，设A (x ,y )为原图中的原点O 也就是出发点,B(x 2,y 2) 是图中的终点，C(x 3,y 3)和D(x 4,y 4)是机器人在转弯的时候在限制的危险区域内（半径为10的圆弧）的切点，圆心是O(x 5,y 5) ，圆的半径为r ，AB 的长度是m ，AO 的长度为n ，BO 的长度为c ，由于要求出角DOC 所以，∠DOC=θ ，∠AOB=α ， ∠AOC=β ， ∠BOD=γ 于是有：

{

a =√(x 2?x )2+(y 2?y )2

b =√(x ?x 5)2+(y ?y 5)2

c =√

(x 2?x ５)2+(y ２?y 5)

2

在三角形AOB 中

α=arccos

ｎ2+ｃ2?ｍ

2

2ｎｃ

在直角三角形AOC 中:

β=arccos r

n ;

在直角三角形BOD 中：

γ=arccos r

c

;

6

所以根据以上的公式关系可以求出：θ=2π?α?β?γ

以及算出L =2?r 2+√c 2?r 2θr （ｍ，ｎ，ｃ可以根据障碍物的数据描述）求得。

六、模型的建立和求解

6.1问题一建模与求解 6.1.1建模过程

第一步：起点到目标点无论中间障碍物有多少，最短路径都应该是若干个线圆组成的结构。在本题中存在障碍物的状况，且障碍物在拐点处的危险区域是一个半径为10的圆弧，所以结合上述证明，我们得知，求两点之间的最短路径中的转弯半径应该按照最小的转弯半径来计算才能达到最优。

第二步：画出可能存在的路线，解决的就是O 到目标点A 的最短路径问题。很显然的一个问题是机器人从m 走的最短路径肯定是优于机器人从n 走的最短路径，给出了三条可能的最短路径，我们可以分别计算出三条可能路径的最短路径的长度，然后进行比较，取最小者就o 是到目标点B 的最优路径。

第三步：解决的是O 到目标点C 的最短路径问题。图中给出了两条可能路径的最短路径，我们同样可以分别计算出两条可能的最短路径，取最小者就是O 到目标点C 的最优路径。

图3

6.1.2模型的求解

（1）原点O 到目标点A 的最短路径为L =471.0372求解方法见(附录1) 过程见(附录

2)

（2）原点O到目标点B的可能路径有三条，即就有三条可能的最短路径。而机器人经过中间一条路径到达B,这条路径由三条线段和两段圆弧组成，直接用三中的解法是结不出来的。于是我们做了如下变换，先求出中间一条直线的中点设为F(180，370)，那么可以采用三中的解法，分别求O到F和F到B的最短路径，然后两段相加，便可求出O到B的最短路径。求解结果为L=853.1121

(3)机器人经过最下边一条路径，同理这条路径由四条直线和三个圆弧组成，同样可以采取②中的变换，分三部分求解。

综合①②③所述，O到目标点B的最短路径为L=853.1121。原点O到目标点C的可能路径由两条，和2中的方法一样，最终求解结果O到目标点C的最短路径为L=1130.281

(4)从原点O→A→B→C→O就是O→A ,A→B,B→C,C→O各段距离的最小路径值的和，解得L=1884.149

6.2模型建立与求解

6.2.1问题二模型

将题目中的原图的一个模型的转化成上图2所示，设A(x,y)为原图中的原点O也就是出发点,B(x2,y2)是图中的终点，C(x3,y3)和D(x4,y4)是机器人在转弯的时候在限制的危险区域内（半径为10的圆弧）的切点，圆心是O(x5,y5)，圆的半径为r，AB的长度是m，AO的长度为n，BO的长度为c，由于要求出∠DOC，所以，∠DOC=θ，∠AOB=α，∠AOC=β，∠BOD=γ

于是有：

{a=√(2)(2) b=√(5)(5) c=√(x2?x5)2+(y2?y5)2

在三角形AOB中

α=arccos n2 c2?m2

2nc

;

8

在直角三角形AOC中:

β=arccos r

n

;

在直角三角形BOD中：

γ=arccos r

c

;

所以根据以上的公式关系可以求出：θ=2π?α?β?γ

以及算出L=√n2?r2+√c2?r2+θr（m,n,c可以根据障碍物的数据描述）将O A的坐标代入上述式子解得θ的值

T=

√m2?r2+√n2?r2(1+e ? . ρ)

到圆心在m点时方程：

T=√802+2102?ρ2+√2202+902?m2+θ?ρ?arccos(√505

5050ρ)(1+e ? . ρ

2)

（1）解法一用Matlab画图，求极限（附录3）解得ρ=11.5得到T=94.5649 到圆心在n点时的方程：

T=

√2302+602?ρ2+√702+2402?ρ2+mρ(1+e ? . ρ2)

ρ=11.5 T=104.2

圆心在mn直线的方程：

T=

√x2+(?x?150)2?ρ2+(300?x)2+[300?(?x?150)2?y2]+θρ(1+10?0.1ρ2)

超越方程图像见附 5当x不变时，y增大时，f也在增大，当y不变，x增大时，f也增大，f与x，y存在递增的关系，图得当ρ=11.5 T=94.5649，

（2）解法二

鱼群模式不同于传统的问题解决方法，它提出了一种新的优化模型，对于状态多或无限状态的环境也不必全部遍历，允许一定的不确定性对于摆脱局部最优，从而寻找全局最优时有帮助的人工鱼群算法，这一模式具备分布处理参数和处置的鲁棒性强等能力。问题的解决时通过自治体在自主活动过程中以某种形式表现出来的。在寻优过程中通常会有两种方式表现出来：一种形式是通过人工鱼最终的分布情况来确定最优解的分布，通常随着寻优过程的进展，人工鱼往往会聚集在极值点的周围，而且全局最优解的即致电周围通常会聚集较多的人工鱼，另一种是在人工鱼的个体状态中表现出来的，即在寻优的过程中，跟踪记录最优个体的状态，就类似于遗传算法采用鱼群算法原理。

鱼群算法步骤如下：

9

第一步：确定种群规模N，在变量可行域内随机生成N个体，设定人工鱼可视为域Visual，步长step,拥挤度因子δ，尝试次数trynumber

第二步：计算初始鱼群各个体适应性值，去最优人工鱼状态机器赋值给公告板。

第三步：个体通过觅食、聚群、追尾行为更新自己，生成鱼群。

第四步：评价所有个体。某个个体优于公告板，则将公告板钢芯为该个体。

第五步：当公告板上最优解达到满意误差界内，算法结束，否则转第三步。于是给出函数式利用鱼群算法如下：

T=√802+2102?ρ2+√2202+902?m2+θ?ρ?arccos(√505

5050ρ)(1+e ? . ρ

2)

s.t x∈[0.300]

y∈[10,300]

fishnum=50; %生成50只人工鱼

MAXGEN=30; %最多迭代次数

trynumber=100;%最多试探次数

visual=1; %感知距离

delta=0.618; %是拥挤度因子

step=0.1; %步长

根据鱼群算法（见附录六）解得ρ=11.5 T=94.9106

在比较以上两种算法，时间T比较接近。算法二更精确，但是程序复杂，所以综合判断可以得知算法一更容易接受。

七、模型的优缺点

模型优点:

1、利用了线圆结构作图，运用CAD软件进行穷举法例举，作图不会出现人为的误差，会更加客观、清晰

2、以最简单的例子圆心在顶点上，最小半径为10个单位的方程，逐次深入建立了机器人行走的时间与路径的长度的关系式。

3、利用了鱼群模型计算机器人运行轨迹。

模型缺点：

1、以O点直接连线A点算出最短时间T=84.85281，得出必有最短时间很接近T=84.85281，未能建立完整的最短时间路径模型，而是取巧用了AUTOCAD软件画出了两弧连接的最短距离，并直接用AUTOCAD测量计算了最短时间T=87.91892。

2、因为实际情况限制速度不可能突变，所以以上模型是在较理想条件下建立的。

八、模型的推广与改进

10

本题模型建立了机器人在避障问题的方法，当障碍物足够多时，利用穷举法则不太实用了，AUTOCAD软件也不现实。因此我们想到另一个模型—最优化模型（迪杰斯特拉算法）。当障碍物为n时，迪杰斯特拉算法是典型最短路径算法，用于计算图或网中某个特定顶点到其他所有顶点的最短路径。主要特点是以起始点为中心向外，层层扩展直到扩展覆盖所有顶点。此算法需算出权值，数据一旦复杂计算麻烦，所以matlab编程好的人便能尝试一下。当半径足够大时(ρ>11.5)，半径越大弧就越接近直线，考虑有两条弧三条等就可以避过障碍物，多条弧时就不能有上述解得，

转弯圆心不知，半径越大，两点之间的弧的长度就近似直线，从而建立弧的模型

九、参考文献

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

[2] 梁国业，廖健平，数学建模，北京：冶金工业出版社，2004

[3] 李晓磊，邵之江，钱积新，鱼群算法，https://www.360docs.net/doc/fd4726286.html,/view/6053f68884868762caaed52e.html ，2012年09月 08日

[4] 谭永基，数学模型，上海：复旦大学出版社，2011

[5] 周培德，计算几何-算法与设计，北京：清华大学出版社，2005

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

[7] 袁玉军，陈婷，韩仁江，超越函数积分的五种解法，https://www.360docs.net/doc/fd4726286.html,/view/eca4e11a10a6f524ccbf8554.htm，2012年09月08日

11

附录：

程序一（问题一）

一function result=zongchang(T,W,V,r)

TV=sqrt((T(1)-V(1))^2+(T(2)-V(2))^2);

TW=sqrt((T(1)-W(1))^2+(T(2)-W(2))^2);

VW=sqrt((V(1)-W(1))^2+(V(2)-W(2))^2);

alpha1=acos((TV^2+VW^2-TW^2)/(2*TV*VW));

alpha2=acos(r/TV);

alpha3=acos(r/VW);

alpha4=2*pi-alpha1-alpha2-alpha3;%alpha4?a×aí??2D???TS1=sqrt(TV^2-r^2);%TS1,TS2?ù?a?2???D??%

S2W=sqrt(VW^2-r^2);

S1S2hu=r*alpha4;

result=TS1+S1S2hu+S2W

程序二（问题一）

二zongchang([0 0],[300 300],[230 60],10)

zongchang([0 0],[300 300],[80 210],10)

zongchang([0 0],[232.5 355],[230 60],10)

zongchang([232.5 355],[227.5 415],[235 300],10) zongchang([227.5 415 ],[185 565],[220 530],10) zongchang([185 565],[100 700],[150 600],10)

zongchang([0 0],[250 490],[230 60],10)

zongchang([250 490],[100 700],[270 680],10) zongchang([0 0],[157.5 355],[80 210],10)

Zongchang([175.5 535.5 ],[227.5 415 ],[235 300] ,10) zongchang([227.5 415],[185 565],[220 530],10) zongchang([185 565],[100 700],[150 600],10)

zongchang([0 0],[157.5 355],[80 210],10)

Zongchang([175.5 535.5 ],[227.5 415 ],[235 300] ,10) Zongchang([227.5 415 ],[100 700 ],[270 680] ,10)

zongchang([300 300],[285 540],[300 400],10) zongchang([285 540 ],[100 700] ,[270 680],10)

12

zongchang([300 300], [260 465] ,[220 530],10) zongchang([260 465],[100 700] ,[220 530],10)

zongchang([0 0],[320 80], [230 60] ,10) zongchang ([320 80],[455 150], [410,100],10) zongchang([445,150],[610 360], [500 200],10) zongchang([610 360],[720 560],[720 520] ,10) zongchang([720 560 ],[700 640] ,[720 640],10)

zongchang([0 0],[320 80], [230 60] ,10) zongchang([320 80],[455 170] ,[410 100],10) zongchang([455 170],[560 380], [500 240],10) zongchang([560 380],[720 560],[720 520 ],10) zongchang([720 560],[700,640], [720 600],10)

zongchang([0 0],[240 220],[80 210],10) zongchang([240 220 ],[470 330] ,[400 330],10) zongchang([470 330 ],[630 435],[540 300],10) zongchang([630 425], [720 560,], [720 520],10) zongchang([720 560],[700,640], [720 600],10)

zongchang([80 680],[385,640],[270,680],10) zongchang([570 660],[690,520],[640,520],10) zongchang([690,520],[720,560],[720,520],10) zongchang([720,560],[700,640],[720,640],10)

zongchang([385 640],[570 ,660],[ 500 600],10) zongchang([80,680],[385,640],[270,680],10) zongchang([385,640],[520,675],[500,600],10) zongchang([520,675],[605,730],[540,730],10) zongchang([605,730],[700,640],[670,730],10)

zongchang([100,700],[320,680],[270,680],10) zongchang([320,680],[495,705],[430,680],10) zongchang([495,705],[605,730],[540,730],10) zongchang([605,730],[700,640],[670,730],10)

zongchang([0 0],[320 80], [230 60] ,10) zongchang ([320 80],[455 150], [410,100],10) zongchang([445,150],[610 360], [500 200],10) zongchang([610 360],[720 560],[720 520] ,10) zongchang([720 560 ],[700 640] ,[720 640],10)

13

zongchang([0 0],[320 80], [230 60] ,10)

zongchang([320 80],[455 170] ,[410 100],10)

zongchang([455 170],[560 380], [500 240],10)

zongchang([560 380],[720 560],[720 520 ],10)

zongchang([720 560], [700 640],[720 600],10)

zongchang([0 0],[290 205],[80 210],10)

zongchang([290 205],[560 380],[500 240],10)

zongchang([560 380],[720 560],[720 520],10)

zongchang([720 560], [700 640],[720 600],10)

zongchang([0 0],[240 220],[80 210],10)

zongchang([240 220 ],[470 330] ,[400 330],10)

zongchang([470 330 ],[630 435],[540 300],10)

zongchang([630 425], [720 560,], [720 520],10)

zongchang([720 560],[700,640], [720 600],10)

zongchang([80 680],[385,640],[270,680],10)

zongchang([385 640],[570 ,660],[ 500 600],10)

zongchang([570 660],[690,520],[640,520],10)

zongchang([690,520],[720,560],[720,520],10)

zongchang([720,560],[700,640],[720,640],10)

zongchang([80,680],[385,640],[270,680],10)

zongchang([385,640],[520,675],[500,600],10)

zongchang([520,675],[605,730],[540,730],10)

zongchang([605,730],[700,640],[670,730],10)

zongchang([100,700],[320,680],[270,680],10)

zongchang([320,680],[495,705],[430,680],10)

zongchang([495,705],[605,730],[540,730],10)

zongchang([605,730],[700,640],[670,730],10)

程序三问题（二）

三clear

syms x y z m n

x=acos((220^2+90^2+80^2+210^2-300^2*2)/(2*sqrt((300-80)^2+90^2)*sqrt(80^ 2+210^2)));

y=acos(m/(sqrt((300-80)^2+90^2)));

z=acos(m/(sqrt(80^2+210^2)));

n=2*pi-x-y-z;

fplot('(sqrt(80^2+210^2-m^2)+sqrt(220^2+90^2-m^2)+m*(8917332864357639/22

14

51799813685248 - acos(1/5650*565^(1/2)*m) -

acos(1/5050*505^(1/2)*m))*(1+exp(10-0.1*m^2)))/5',[10,50])

syms y m

y=(sqrt(80^2+210^2-m^2)+sqrt(220^2+90^2-m^2)+m*(2*pi -

acos(1/5050*505^(1/2)*m) - acos(1/5650*565^(1/2)*m) - (5502580100633845*i)/4503599627370496)*(1+exp(10-0.1*m^2)))/5;

diff(y)

solve(ans)

y=(sqrt(80^2+210^2-11.5^2)+sqrt(220^2+90^2-11.5^2))/5+11.5*(891733286435 7639/2251799813685248 - acos(1/5650*565^(1/2)*11.5) -

acos(1/5050*505^(1/2)*11.5))*(1+exp(10-0.1*11.5^2))/5

程序四问题（二）

四

syms x y z m n

x=acos((230^2+60^2+70^2+240^2-2*300^2)/(2*sqrt(270^2+60^2)*sqrt(70^2+240 ^2)));

y=acos(m/(sqrt((230)^2+60^2)));

z=acos(m/(sqrt(270^2+240^2)));

n=2*pi-x-y-z;

fplot('((sqrt(230^2+60^2-m^2)+sqrt(70^2+240^2-m^2)+m*(598918793813111/14 0737488355328 - acos(1/5650*565^(1/2)*m) -

acos(17592186044416/6355139163046355*m)))*(1+exp(10-0.1*m^2)))/5',[10,50 ])

syms y m

y=((sqrt(230^2+60^2-m^2)+sqrt(70^2+240^2-m^2))+m*(598918793813111/140737 488355328 - acos(1/5650*565^(1/2)*m) -

acos(17592186044416/6355139163046355*m))*(1+exp(10-0.1*m^2)))/5;

diff(y)

solve(ans)

y=((sqrt(230^2+60^2-11.5^2)+sqrt(70^2+240^2-11.5^2)+11.5*(59891879381311 1/140737488355328 - acos(1/5650*565^(1/2)*11.5) -

acos(17592186044416/6355139163046355*11.5)))*(1+exp(10-0.1*11.5^2)))/5;程序六问题（二）

附录六

15

clc

clear all

close all

tic

figure(1);hold on

for i=1:1000

X=8+(16-8)*i/1000;

f=-(sqrt(50500-X^2)+sqrt(56500-X^2)+X*(2*pi-acos(X/224.7220505)-acos(X/2 37.6972865)-acos((250^2+237.6972865^2-424.2640687^2)/(2*224.7220505*237. 6972865)))*(1+exp(10-0.1*X^2)))/5;

plot(X,f);

hold on

end

%% 2?êyéè??

fishnum=50; %éú3é50??è?1¤ó?

MAXGEN=30; %×??àμü′ú′?êy

try_number=100;%×??àê?ì?′?êy

visual=1; %?D?a?àà?

delta=0.618; %óμ?·?èòò×ó

step=0.1; %2?3¤

%% 3?ê??ˉó?èo

lb_ub=[8,16,1];

X=AF_init(fishnum,lb_ub);

LBUB=[];

for i=1:size(lb_ub,1)

LBUB=[LBUB;repmat(lb_ub(i,1:2),lb_ub(i,3),1)];

end

gen=1;

BestY=-1*ones(1,MAXGEN); %??2??D×?ó?μ?oˉêy?μ

BestX=-1*ones(1,MAXGEN); %??2??D×?ó?μ?×?±?á?

besty=-100; %×?ó?oˉêy?μ

Y=AF_foodconsistence(X);

while gen<=MAXGEN

fprintf(1,'%d\n',gen)

for i=1:fishnum

%% ??èoDD?a

[Xi1,Yi1]=AF_swarm(X,i,visual,step,delta,try_number,LBUB,Y);

%% ×·?2DD?a

[Xi2,Yi2]=AF_follow(X,i,visual,step,delta,try_number,LBUB,Y);

if Yi1>Yi2

16

X(:,i)=Xi1;

Y(1,i)=Yi1;

else

X(:,i)=Xi2;

Y(1,i)=Yi2;

end

end

[Ymax,index]=max(Y);

figure(1);

plot(X(1,index),Ymax,'.','color',[gen/MAXGEN,0,0]) if Ymax>besty

besty=Ymax;

bestx=X(:,index);

BestY(gen)=Ymax;

[BestX(:,gen)]=X(:,index);

else

BestY(gen)=BestY(gen-1);

[BestX(:,gen)]=BestX(:,gen-1);

end

gen=gen+1;

end

plot(bestx(1),besty,'ro','MarkerSize',100)

xlabel('x')

ylabel('y')

title('ó?èo??·¨μü′ú1y3ì?D×?ó?×?±êò??ˉ')

%% ó??ˉ1y3ìí?

figure

plot(1:MAXGEN,BestY)

xlabel('μü′ú′?êy')

ylabel('ó??ˉ?μ')

title('ó?èo??·¨μü′ú1y3ì')

disp(['×?ó??aX￡o',num2str(bestx,'%1.4f')])

disp(['×?ó??aY￡o',num2str(besty,'%1.4f')])

toc

17