机器人避障问题国家一等奖论文-毕设论文.doc

D题 机器人避障问题
摘 要
本文综合运用分析法、图论方法、非线性规划方法，讨论了机器人避障最短路径和最短时间路径求解问题。
针对问题一，首先，通过分析，建立了靠近障碍物顶点处转弯得到的路径最短、转弯时圆弧的半径最小时和转弯圆弧的圆心为障碍物的顶点时路径最短、转弯在中间目标点附近时，中间目标点位于弧段中点有最短路径的三个原理，基于三个原理，其次对模型进行变换，对障碍物进行加工，扩充为符合条件的新的区域并在转弯处圆角化构成障碍图，并通过扩充的跨立实验，得到切线和圆弧是否在可避障区的算法，第三，计算起点、中间目标点和最终目标点和各圆弧及圆弧之间的所有可避障切线和圆弧路径，最后给这些定点赋一个等于切线长度或弧度的权值构成一个网络图，然后利用Dijkstra算法求出了O-A、O-B，O-C的最短路径为O-A：471.0372个单位，O-B：853.7001个单位，O-C：1086.0677个单位；对于需要经中间目标点的路径，可运用启发规则分别以相邻的目标点作为起点和终点计算，确定路径的大致情况，在进一步调整可得到O-A-B-C-O的最短路径为2748.699个单位。
针对问题二，主要研究的是由出发点到达目标点A点的最短时间路径，我们在第一问的基础上考虑路径尽可能短且圆弧转弯时的圆弧尽量靠近障碍物的顶点，即确定了圆弧半径最小时的圆弧内切于要确定的圆弧时存在最小时间路径，建立以总时间最短为目标函数，采用非线性规划模型通过Matlab编程求解出最短时间路径为最短时间路程为472.4822个单位，其中圆弧的圆心坐标为，最短时间为94.3332秒。圆弧两切点的坐标分别为（70.88,212.92）、（77.66,219.87）。
关键字： Dijkstra算法 跨立实验 分析法 非线性规划模型

一.问题的重述
图是一个800×800的平面场景图，在原点O(0, 0)点处有一个机器人，它只能在该平面场景范围内活动。图中有12个不同形状的区域是机器人不能与之发生碰撞的障碍物，障碍物的数学描述如下表：
编号
障碍物名称
左下顶点坐标
其它特性描述
1
正方形
(300, 400)
边长200
2
圆形

圆心坐标(550, 450)，半径70
3
平行四边形
(360, 240)
底边长140，左上顶点坐标(400, 330)
4
三角形
(280, 100)
上顶点坐标(345, 210)，右下顶点坐标(410, 100)
5
正方形
(80, 60)
边长150
6
三角形
(60, 300)
上顶点坐标(150, 435)，右下顶点坐标(235, 300)
7
长方形
(0, 470)
长220，宽60
8
平行四边形
(150, 600)
底边长90，左上顶点坐标(180, 680)
9
长方形
(370, 680)
长60，宽120
10
正方形
(540, 600)
边长130
11
正方形
(640, 520)
边长80
12
长方形
(500, 140)
长300，宽60
在图1的平面场景中，障碍物外指定一点为机器人要到达的目标点(要求目标点与障碍物的距离至少超过10个单位)。规定机器人的行走路径由直线段和圆弧组成，其中圆弧是机器人转弯路径。机器人不能折线转弯，转弯路径由与直线路径相切的一段圆弧组成，也可以由两个或多个相切的圆弧路径组成，但每个圆弧的半径最小为10个单位。为了不与障碍物发生碰撞，同时要求机器人行走线路与障碍物间的最近距离为10个单位，否则将发生碰撞，若碰撞发生，则机器人无法完成行走。
机器人直线行走的最大速度为个单位/秒。机器人转弯时，最大转弯速度为，其中是转弯半径。如果超过该速度，机器人将发生
侧翻，无法完成行走。
请建立机器人从区域中一点到达另一点的避障最短路径和最短时间路径的数学模型。对场景图中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的最短时间路径。
注：要给出路径中每段直线段或圆弧的起点和终点坐标、圆弧的圆心坐标以及机器人行走的总距离和总时间。

图 800×800平面场景图
二.问题的分析
本问题的难点在于机器人要到达指定的目标点需要满足以下两个约束条件：
1. 机器人的行走路径由直线段和圆弧组成，不能折线转弯，转弯路径由与直线路径相切的一段圆弧组成，且每个圆弧的半径最小为10个单位；
2. 要求机器人行走线路与障碍物间的最近距离为10个单位。
因此，我们在建立模型求解机器人到达目标点的最短路径时需优先考虑这两个约束条件。
首先我们可以根据约束条件将机器人行走的危险区域进行扩张，即所有的障碍物都向外扩张了10个单位。机器人所走的路径都是直线与圆弧的构成，故存在线和弧、弧与弧之间的切线，可以考虑在所有代表出发点与其它圆弧之间切线的顶点与源点连成一条边，权值均为0，同理在所有代表目标点到其它圆弧切线的顶点与终点连成一条边，权值均为0，这样题目就转化成了求源点到达终点之间的最短路径问题了。
对于问题二，要求最短时间路径，则要考虑的是要以最大速度行走。直线行走时就是最大速度，但在转弯圆弧处因为转弯半径越小，行走的速度就越慢，则需在第一问的前提下增大圆弧的半径，则圆弧的转弯半径和圆心都不确定，通过建立模型确定O-A之间的转弯时圆弧的半径和圆心。

三、模型的基本假设
根据对该问题的分析，本文对所建立的模型提出以下基本假设：
1.机器人的性能足够好，能准确地沿圆弧转弯
2.机器人行走过程中不会意外停止
3.图中所给数据障碍物都是真实的
4.机器人能够抽象成点来处理
四、符号说明
指每段路径的长度
：机器人直线行走时的最大速度
：机器人的最大转弯速度
：每段路径所用的时间

五、模型的建立与求解
5.1模型的准备
模型的准备一
首先给出以下三个前提及其证明：
1.靠近障碍物顶点处转弯得到的路径最短
如图1所示：

图1

设A为起点，B为终点，矩形的阴影部分是障碍区，C为障碍物的顶点，D为障碍物外任意一点，连接AD，BD，AC，BC延长交AD于E点，因为在三角形中两边之和大于第三边，所以有：
BD+DE>BE；BE=BC+CE；
AE+CE>AC；AD=AE+DE；
两个不等式相加，得：AE+DE+AE+CE>BE+AC；
即：BD+AD>AC+BC.
于是，得证由A到B在顶点C处转弯时为最短路径。
2.转弯时圆弧的半径最小时路径最短和转弯圆弧的圆心为障碍物的顶点时路径最短
要证明机器人转弯时圆弧的半径最小时路径最短，则可以把P到A的路径比拟成一条弹性绳，由于路径要绕过障碍物，故要拉长弹性绳以绕过障碍物，由于机器人转弯时只能是走圆弧路径且有圆弧半径约束，即可以在障碍物的顶点处加入圆环，将其圆心固定在顶点处，由于弹性绳的弹性势能和伸长量存在如下能量关系：，故弹性绳弹性势能最小时伸长量最小，此时正是路径最短的情况：
如图2所示：

图2
根据最小势能原理可知，当弹性体平衡时，系统势能最小。即弹性体在自由条件下，有由高势能向低势能转化的趋势。现在将圆环看成也有弹性，在如图所示的条件下为初始状态。圆环受力为F，此时圆环有缩小的趋势，随着圆环的缩小系统趋于平衡，弹性绳有最小势能。由能量守恒也可以说明，弹性绳的弹性势能转化为弹性圆环的弹性势能，于是弹性绳的弹性势能减小。
因此，随着圆环的半径的减小，弹性绳的势能减小，即最短路径变短。所以转弯时圆弧的半径最小时路径最短。转弯圆弧的圆心为障碍物的顶点时路径最短。
3.转弯在中间目标点附近时，中间目标点位于弧段中点有最短路径
当机器人从O点出发不能直接到达另一目标点，而需经过中间目标点A时，则必须在A点转弯处形成一个圆弧，中间目标点位于弧段中点时且半径最小时有最短路径，下面给予证明：由上述2的证明可仍然采用此证法，设从两转向处引
出一根弹性绳，要使其经过A 点，如图所示，在其自然伸长的前提下，必定会在A 点出弯折。由于机器人有转弯时圆弧的半径最小时路径最短，从两转向处引出一根弹性绳，要使其经过A 点，如图所示，在其自然伸长的前提下，必定会在A 点出弯折。由于机器人有最小转弯半径约束，故在A 处加一最小转弯半经的圆，即r=10个单位，使其受力平衡时，即为路径的最优解。
根据以上证明作出如下图形(图3)：

图3

作法：过A点分别作圆上M,N的切线，作的角平分线，在上距离A点10个单位处选取一点O，以O为圆心，10个单位为半径，作圆O，再过M,N分别作圆O的切线，即为路径。
模型准备二：
切线的剔除：
基于上述准备一的三个前提条件和题中给出的两约束条件即（1）机器人的行走路径由直线段和圆弧组成，不能折线转弯，转弯路径由与直线路径相切的一段圆弧组成，且每个圆弧的半径最小为10个单位（2）要求机器人行走线路与障碍物间的最近距离为10个单位。故存在一些不符合要求的切线，要将这些切线剔除，可做示意图如下：

图4

如图所示，阴影部分为障碍物，用圆来代表弧，AC、BC、AD、BD、EH、FG分别如图示的切线，由于AD、BC、BD、FG、EH经过障碍物，所以要剔除，AC为符合要求切线。(图中应该剔除的切线用虚线表示)。
模型准备三
由上准备一和准备三可知，对于机器人转弯的路径按照一个圆弧、两个圆弧、三个圆弧等多个圆弧枚举出转弯的情况：
转过了一个圆弧,如图5:

图5
如图，设()，O()，()，、的坐标分别为和，，，, 则: ,，，
又=，可知：
，
，
，
故B到A的路径长度为： 。
其中，，，，
即。
并且利用相切的关系可以得出下面两个方程：
可以解出切点，同理可以解出切点。
转过了两个圆弧，如图6：
情况一：

图6

假设两圆心坐标分别为和,半径均为，O点坐标为，则容易求得:

这样我们就可以利用只有一个圆弧中的方法，先求A到O，再求O到B，于是分两段就可以求解。同理如果有更多的转弯，同样可按照此种方法分解。
情况二：

图7
设圆心坐标分别为和,半径均为，，的长为,长为,长为，长为，长为；，，，，，；由这些点的几何关系可以得到：

其中：，，，， ，.
则，，，，，
所以： .
情况三:
当障碍物本身是个圆弧时：

图8
此种情况计算路径的方法同情况二。
综上几种情况的讨论，可以知道：机器人的避障路径都是由以上几种情况组成，故我们在求由O 点到各目标点时路径长度都可以将路径分解成以上几种情况来求解。

5.2模型的建立与求解
5.2.1问题一的模型建立与求解
模型的建立：
题中给出了12个多边形障碍物的区域，故我们根据图中的约束条件利用Matlab(程序见附录一)将障碍物扩充为如图9：
图9
即形成一个新的障碍区域，先求出所有的切线，包括出发点和目标点到所有圆弧的切线以及所有圆弧与圆弧之间的切线，给这些定点赋一个等于切线长度的权值，如果某两条切线有一个公切圆弧，则代表这两条曲线的顶点是一条直线的两个端点，边上的权值等于这两条切线之间的劣弧长度。那么在所有代表出发点与其它圆弧之间切线的顶点与源点连成一条边，权值均为0，同理在所有代表目标点到其它圆弧切线的顶点与终点连成一条边，权值均为0，这样题目就转化成了求源点到达终点之间的最短路径问题了，这里最短路径就是指经过所有顶点与边的权值之和最小。
对于最短路径的求解，有以下步骤：
1）画出出发点和目标点和各个圆弧的切线，以及圆弧与圆弧之间的切线，
但是切线有可能经过障碍物的内部或危险区域，也可能出现切线重复的状况，既有很多不合理的切线，再根据约束条件我们对模型做了以下修正：
l 检验切线两个端点是否在障碍物内部。
l 检验切线两切线是否相交
l 检验圆弧所对应的圆心，即障碍物的顶点到切线的距离是否小于10。
除上述三个指标外还需考虑一种特殊情况，即公切线在一个圆弧上只有一个切点(如图10)
因此我们作如下规定:如果某条切线与某段圆弧相切，且切点不在切线的端点上，则该切线为不合理。

图10
2）然后将不合理的切线剔除，假设剔除过后有m条合理切线，那么就有个顶点，设这些点的权值（），即第条合理曲线的长度。为边的权值，即第条弧的长度。
3）然后把路径确定下来，按照求得权值矩阵给图中的顶点及边长赋值。利用Dijkstra算法对其求解出最短路径。
根据上述的模型通过编程求解和几何求解求出了最短路径，并且求出了最短路径中每段直线段或圆弧的起点和终点坐标、圆弧的圆心坐标。

模型的结果
通过编程求解得出最短路径的线路坐标和最短路径：

表一 各线路的最短路径

路径

最短路径的坐标线路

最短路径的长度(个单位)

O-A

471.0372

O-B

853.7001

O-C

1086.0677

O-A-B-C-O

2748.699

其中：
O-A:
圆弧圆心坐标为（80,210）。
O-B:
圆弧圆心坐标依次为：（60,300），（345,210），（220,470），（220,530），（150,600）。
O-C:
圆弧圆心坐标依次为（230,60），（410,100），（500,200），（720,520），（720,600）。
O-A-B-C-O:
圆弧圆心的坐标：A（299.327,309.9773），B（108.0849,694.115），C（708.990,644.3794）（其他圆心在障碍物的顶点时未列出，只列出了在A,B,C三个目标点的圆心坐标）。
各路径最短线路图如下：

图11 O-A的最短路径

图12 O-B的最短路径

图13 O-C的最短路径

图14 O-A-B-C-O最短路径图
5.2.2问题二模型的建立与求解
模型分析
由第一问的模型可知，当转弯圆弧的圆心在顶点时，行走路径越靠近两点连线，且圆弧半径越小路径越短，所花的时间也就越少；但根据最大转弯速度为，其中是转弯半径，圆弧半径越小转弯的速度也就越慢，所花的时间又更长，基于这两者结合，将转弯路径的圆弧与问题一中O点到A
点圆弧相内切，如图15所示

图15
设圆弧的圆心为，问题一中的顶点圆心为（80,210），
,,,,,,;则有，，，，，，，， ，故可列出时间的计算表达式为，，
要确定出时间的最短路径，我们采用非线性规划模型求解，建立时间最短的目标函数。
模型的建立与求解
机器人从O点到A点，由第一问得最短路径，设第一段线段所用时间为,第二段圆弧段所用时间为，第三段线段所用时间为,为两圆心的距离，那么
目标函数可以表示为：

利用Matlab编程求得总时间最小为秒，最小时间总路程为472.4822个单位，其中转弯处圆弧圆心坐标为，此时圆弧的半径为11.197个单位，圆弧两切点的坐标分别为（70.88,212.92）、（77.66,219.87）。

六、模型的评价
模型的优点：
1. 模型简单易懂，对一般过路障问题实用。
2. 问题一中利用权值网络对问题的整体把握，有全局包揽局部的效果，避免了只考虑小部分的情况。
3. 对问题建立模型前，通过一些定理的证明及模型准备为模型的建立提供了理论的支撑，避免了出现泛泛而谈。
模型的缺点：
1. 此模型利用整体把握的考虑，但由于它遍历计算的节点很多，计算具体数值时可能存在一定的误差，效率也会更低。
2. 当障碍物为其他极其不规则的形状时，模型需考虑多几方面的情形，要对模型进行适当的修改。
七、参考文献
[1]姜启源，数学模型，北京，高等教育出版社 ，2003
[2] 尤承业，解析几何，北京，北京大学出版社，2004
[3]周培德，计算几何—算法与设计，北京清华大学出版社，2005
[4]王文波，数学建模及其基础知识详解，武汉，武汉大学出版社，2005
[5]祝红芳 王从庆，机器人路径规划的元胞自动机算法，江西，江西科学，第27卷第1期，36-40,2009.2
[6]龚劬，图论与网络最优化算法，重庆，重庆大学出版社，2009
[7]李明哲,图论及其算法,北京，机械工业出版社，2010
[8]王桂平,图论算法理论、实现及运用，北京，北京大学出版社，2011

附 录
附录一

路径图的Matlab程序:

运行的约束条件：
open fig1.fig
text(0,0,'O','fontsize',12);
text(300,300,'A','fontsize',12);
R=10;
O=[80 210];
P1=[70.5060 213.1406];
P2=[76.6054 219.4066];
line([0 70.5060],[0 ,213.1406],'LineWidth',2);

an1=angle(complex(P1(1)-O(1),P1(2)-O(2)));
an2=angle(complex(P2(1)-O(1),P2(2)-O(2)));
if an2<0
an2=an2+2*pi;
end
an=linspace(an1,an2,50);
plot(O(1)+R*cos(an),O(2)+R*sin(an),'LineWidth',2);

dot1=[76.6054,219.4066];
dot2=[300,300];
line([dot1(1) dot2(1)],[dot1(2) dot2(2)],'LineWidth',2);

OA
open fig1.fig
text(0,0,'O','fontsize',12);
text(300,300,'A','fontsize',12);
R=10;
O=[80 210];
P1=[70.5060 213.1406];
P2=[76.6054 219.4066];
line([0 70.5060],[0 ,213.1406],'LineWidth',2);

an1=angle(complex(P1(1)-O(1),P1(2)-O(2)));
an2=angle(complex(P2(1)-O(1),P2(2)-O(2)));
if an2<0
an2=an2+2*pi;
end
an=linspace(an1,an2,50);
plot(O(1)+R*cos(an),O(2)+R*sin(an),'LineWidth',2);

dot1=[76.6054,219.4066];
dot2=[300,300];
line([dot1(1) dot2(1)],[dot1(2) dot2(2)],'LineWidth',2);

OB
open fig1.fig
text(0,0,'O','fontsize',12);
text(100,700,'B','fontsize',12);
R=10;
O=[0 0];

dot1=O;
dot2=[50.1353,301.6396];
linePic(dot1,dot2);

O1=[60,300];
P1=[50.1353,301.6396];
P2=[50.1353,301.6396];
arcPic(O1,P1,P2,R)

dot1=[51.6797,305.5746];dot2=[141.6979,439.6089];
linePic(dot1,dot2);

P1=[141.6979,439.6089];
P2=[150,444.0343];
O1=[150,435];
arcPic(O1,P1,P2,R);

dott=[150,444.0343, 222.1970,460.2429];
dot1=dott(1:2);dot2=dott(3:4);
linePic(dot1,dot2);

arcc=[222.1970,460.2429, 230,470, 220,470];
P1=arcc(1:2);
P2=arcc(3:4);
O1=arcc(5:6);
arcPic(O1,P1,P2,R);

dott=[230 470 230 530];
dot1=dott(1:2);dot2=dott(3:4);
linePic(dot1,dot2);

arcc=[ 230 530 225.4967 538.3538 220,530 ];
P1=arcc(1:2);
P2=arcc(3:4);
O1=arcc(5:6);
arcPic(O1,P1,P2,R);

dott=[225.4967 538.3538 144.5033 591.6462];
dot1=dott(1:2);dot2=dott(3:4);
linePic(dot1,dot2);

arcc=[ 144.5033 591.6462 140.6916 596.3458 150 600 ];
P1=arcc(1:2);
P2=arcc(3:4);
O1=arcc(5:6);
arcPic(O1,P1,P2,R);

dott=[140.6916 596.3458 100 700];
dot1=dott(1:2);dot2=dott(3:4);
linePic(dot1,dot2);

% £¨0,0£© £¨50.1353,301.6396£©
% £¨50.1353,301.6396£© £¨51.6797,305.5746£© £¨60,300£© 10
% £¨51.6797,305.5746£© £¨141.6979,439.6089£©
% £¨141.6979,439.6089£© £¨150,444.0343£© £¨150,435£© 10
% £¨150,444.0343£© £¨222.1970,460.2429£©
% £¨222.1970,460.2429£© £¨230,470£© £¨220,470£© 10
% £¨230,470£© £¨230,530£©
% £¨230,530£© £¨225.4967,538.3538£© £¨220,530£© 10
% £¨225.4967,538.3538£© £¨144.5033,591.6462£©
% £¨144.5033,591.6462£© £¨140.6916,596.3458£© £¨150,600£© 10
% £¨140.6916,596.3458£© £¨100,700£©

OABCO
open fig1.fig
text(0,0,'O','fontsize',12);
text(700,640,'C','fontsize',12);
R=10;
O=[0 0];

data={ [0 0 70.506 213.1406 ]
[70.506,213.1406 76.6215,219.412 80,210 ]
[76.6215 219.412 294.4385,298.5652 ]
[294.4385,298.5652 289.9197,306.5856 290.8727 304.0857 ]
[289.9197 306.5856 229.4072,533.3617 ]
[229.4072,533.3617 225.4967,538.3538 220,530 ]
[225.4967,538.3538 144.5033,591.6462 ]
[144.5033,591.6462 140.865,595.9316 150,600 ]
[140.865,595.9316 98.9499,690.0466 ]
[98.9499,690.0466 108.9534,704.0772 108.0849,694.115 ]
[ 108.9534,704.0772 270.8686,689.9622 ]
[ 270.8686,689.9622 272 689.798 270,680 ]
[ 272 689.798 368,670.798 ]
[ 368,670.798 370,670 370,680]
[ 370,670 430,670 ]
[430,670 435.5878 671.8432 430,680 ]
[435.5878 671.8432 534.4122 738.2932 ]
[534.4122 738.2932 540,740 540,730 ]
[ 540,740 670,740 ]
[670,740 679.7738,732.1149 670,730 ]
[679.7738,732.1149 699.2162,642.2645 ]
[ 699.2162,642.2645 701.3151,637.9688 708.99,644.3794 ]
[ 701.3151,637.9688 710.2187 602.0799 ]
[710.2187 602.0799 730,600 720,600 ]
[ 730,600 730,520 ]
[730,530 727.9377 513.9178 720,520 ]
[ 727.9377 513.9178 492.0623 206.0822 ]
[ 492.0623 206.0822 491.6552,205.5103 550,450 ]
[ 491.6552,205.5103 418.3448,94.4897 ]
[418.3448,94.4897 412.1693,90.2381 400,330 ]
[ 412.1693,90.2381 232.1693,50.2381 ]
[ 232.1693,50.2381 232.1149,50.2262 230,60 ]
[ 232.1149,50.2262 0,0 ]}
for i=1:length(data)
if length(data{i})==4
dot1=data{i}(1:2);dot2=data{i}(3:4);
linePic(dot1,dot2);
else
P1=data{i}(1:2);
P2=data{i}(3:4);
O1=data{i}(5:6);
arcPic(O1,P1,P2,R);
end
end

运行的约束条件：
open fig1.fig
text(0,0,'O','fontsize',12);
text(300,300,'A','fontsize',12);
R=10;
O=[80 210];
P1=[70.5060 213.1406];
P2=[76.6054 219.4066];
line([0 70.5060],[0 ,213.1406],'LineWidth',2);

an1=angle(complex(P1(1)-O(1),P1(2)-O(2)));
an2=angle(complex(P2(1)-O(1),P2(2)-O(2)));
if an2<0
an2=an2+2*pi;
end
an=linspace(an1,an2,50);
plot(O(1)+R*cos(an),O(2)+R*sin(an),'LineWidth',2);

dot1=[76.6054,219.4066];
dot2=[300,300];
line([dot1(1) dot2(1)],[dot1(2) dot2(2)],'LineWidth',2);

OA
open fig1.fig
text(0,0,'O','fontsize',12);
text(300,300,'A','fontsize',12);
R=10;
O=[80 210];
P1=[70.5060 213.1406];
P2=[76.6054 219.4066];
line([0 70.5060],[0 ,213.1406],'LineWidth',2);

an1=angle(complex(P1(1)-O(1),P1(2)-O(2)));
an2=angle(complex(P2(1)-O(1),P2(2)-O(2)));
if an2<0
an2=an2+2*pi;
end
an=linspace(an1,an2,50);
plot(O(1)+R*cos(an),O(2)+R*sin(an),'LineWidth',2);

dot1=[76.6054,219.4066];
dot2=[300,300];
line([dot1(1) dot2(1)],[dot1(2) dot2(2)],'LineWidth',2);

OB
open fig1.fig
text(0,0,'O','fontsize',12);
text(100,700,'B','fontsize',12);
R=10;
O=[0 0];

dot1=O;
dot2=[50.1353,301.6396];
linePic(dot1,dot2);

O1=[60,300];
P1=[50.1353,301.6396];
P2=[50.1353,301.6396];
arcPic(O1,P1,P2,R)

dot1=[51.6797,305.5746];dot2=[141.6979,439.6089];
linePic(dot1,dot2);

P1=[141.6979,439.6089];
P2=[150,444.0343];
O1=[150,435];
arcPic(O1,P1,P2,R);

dott=[150,444.0343, 222.1970,460.2429];
dot1=dott(1:2);dot2=dott(3:4);
linePic(dot1,dot2);

arcc=[222.1970,460.2429, 230,470, 220,470];
P1=arcc(1:2);
P2=arcc(3:4);
O1=arcc(5:6);
arcPic(O1,P1,P2,R);

dott=[230 470 230 530];
dot1=dott(1:2);dot2=dott(3:4);
linePic(dot1,dot2);

arcc=[ 230 530 225.4967 538.3538 220,530 ];
P1=arcc(1:2);
P2=arcc(3:4);
O1=arcc(5:6);
arcPic(O1,P1,P2,R);

dott=[225.4967 538.3538 144.5033 591.6462];
dot1=dott(1:2);dot2=dott(3:4);
linePic(dot1,dot2);

arcc=[ 144.5033 591.6462 140.6916 596.3458 150 600 ];
P1=arcc(1:2);
P2=arcc(3:4);
O1=arcc(5:6);
arcPic(O1,P1,P2,R);

dott=[140.6916 596.3458 100 700];
dot1=dott(1:2);dot2=dott(3:4);
linePic(dot1,dot2);

% £¨0,0£© £¨50.1353,301.6396£©
% £¨50.1353,301.6396£© £¨51.6797,305.5746£© £¨60,300£© 10
% £¨51.6797,305.5746£© £¨141.6979,439.6089£©
% £¨141.6979,439.6089£© £¨150,444.0343£© £¨150,435£© 10
% £¨150,444.0343£© £¨222.1970,460.2429£©
% £¨222.1970,460.2429£© £¨230,470£© £¨220,470£© 10
% £¨230,470£© £¨230,530£©
% £¨230,530£© £¨225.4967,538.3538£© £¨220,530£© 10
% £¨225.4967,538.3538£© £¨144.5033,591.6462£©
% £¨144.5033,591.6462£© £¨140.6916,596.3458£© £¨150,600£© 10
% £¨140.6916,596.3458£© £¨100,700£©

OABCO
open fig1.fig
text(0,0,'O','fontsize',12);
text(700,640,'C','fontsize',12);
R=10;
O=[0 0];

data={ [0 0 70.506 213.14