西电数模大作业

1、 关于中国南海岛礁巡航路线方案的探讨 姓名：刘发强 2016年11月8日 摘要南海岛礁巡航路线方案的确定实质上是一个求取平面有障碍区域两点之间最短路径问题和TSP问题的混合。针对前者，本文借助ESPO算法的思想，将外国非法侵占岛礁所确定的12海里敏感区域抽象为正多边形，在对所有顶点及起始点之间任一条路径的合法性判定（借助计算几何中点、线段、多边形的位置关系）之后，利用Dijkstra算法生成障碍密集区（南沙群岛）我方任意两岛礁之间的最短路径，后对我方全体（西沙、南沙、中沙）岛礁利用Floyd算法生成最短路径完全图的邻接矩阵，为后续处理奠定基础。针对后者，利用自适应的蚁群算法求取最优解，经检验效

2、果较好。同时借助openCV提供的图像处理和画图函数对算法的运行和结果进行了可视化呈现，易于理解调试。理论上，本文的解决方案对于这类问题有很广泛的适用性。关键词：TSP问题；ESPO算法；Dijkstra算法；Floyd算法；openCV；蚁群算法；算法可视化 绪论 随着中国海权意识的逐步觉醒和军事力量的日益强大以及美国重返亚太战略的实施，南海态势在多方的博弈下波谲云诡，不安定不稳定性极大。为了维护我国的主权和领土完整，对我国实际控制的岛礁进行周期性的补给和巡航势在必行。因此，在周围情况极其复杂的南海确定一个最优的巡航路径非常关键。本文主要借助ESPO算法和蚁群算法对我国南海实际控制的包括西沙

3、、中沙、部分南沙在内的23个岛礁进行了路径规划，比较圆满的解决了上述问题，且具有良好的可扩展性和适用性。本文的基本假设：1、 地球表面近似为球面，在经线上，每纬度约111千米，在纬线上，每经度约为111*cos千米，其中为纬度。由于所考查区域（纬度从东经111°到东经117°，北纬8°到北纬18°）面积相对较小且在低纬度，故近似为平面即每经度111千米。利用百度地图的测距工具和坐标拾取系统得出的结果与对平面上两点欧几里得距离相比最大相对误差低于1%，这对于本问题的解决足够准确。2、 岛礁抽象为其几何中心处一点。3、 本文考察西沙群岛、中沙群岛的所有岛礁和

4、南沙群岛的7个岛礁。由于中沙群岛大部分在中部呈一紧凑环装分布，为了简化且不失准确性，将大环礁抽象为环北一点和环南一点。4、 本文所规划的路径除不可避免之外不进入外国侵占岛礁所形成的12海里敏感区域，由于东门礁、南薰礁距离鸿庥岛和九章群礁过近，在此两处将12海里放宽为6海里。5、 由于部分岛礁之间平面欧几里得距离远大于其面积尺度，故在无障碍区即西沙、中沙我方岛礁的最短距离可以由平面欧几里得距离代替。并且最优解中，可以预见路径穿越岛礁的可能性极小，再次印证上述假定的合理性。6、 在构造最短路径完全图的过程中，对于非南沙与非南沙之间、非南沙和南沙之间其距离认定为其平面欧几里得距离，南沙群岛各岛礁距离

5、利用ESPO算法生成，由于南沙距离非南沙相对于其尺度很大，故近似比较准确。方案形成过程南海岛礁经纬度采集和部分岛礁之间的测距借助百度地图及其提供的坐标拾取系统和测距工具，获取了我方23个岛礁、他方18个岛礁的经纬度信息和用于连通南沙我方部分岛礁（障碍外部）与中南暗沙、大环礁、中建岛、盘石屿等岛礁的距离。并利用MATLAB计算了我方非南沙各岛礁之间的欧氏距离，从而生成了一部分最短路径完全图。由于表格较大，不在此呈现，详见表格。图 1南海部分岛礁南沙群岛有障碍区域我方岛礁最短路径的形成由假设知，各敏感区域是由以外国侵占岛礁几何中心为圆心半径为12海里的区域，为凸集。若可行区域的边界是光滑曲面。则最

6、短路径必由下列弧组成，它们或者是空间中的自然最短曲线，或者是可行区域的边界弧。而且，组成最短路径的各段弧在连接点处必定相切。上述定理由J.W.Craggs在1973年证明。由于光滑曲线在计算机中不便于存储和处理，本文将敏感区域近似为其内接正多边形（本文中最终使用以12边形近似的数据），并基于上述定理近似推断在障碍区两点之间的最短路径必有以下线段组成，他们或是障碍多边形的边，或是障碍多边形各顶点及起始点之间的合法线段组成，其中合法线段为不与任何障碍多边形相交或在其外部的线段。因此，在上述分析的支持下，试图用计算机解决这个平面简单凸多边形障碍区域最短路径问题首先要解决如何判断平面上点、线段与多边形

7、的位置关系。这是一个基本的计算几何问题，对于线段和多边形位置关系的判定，首先要解决线段与线段之间、点与多边形之间的位置关系判定方法。对于点与多边形，可以利用射线法，即由该点为端点任意引一条射线，若射线与多边形各边相交奇数次，则在内部，其他在外部（在边上也认为在外部）；对于线段与线段之间，可利用正交变换（便于处理）将其中一条变换到X轴且一段与原点重合，进而判断另一条线段与X轴的交点与第一条线段的关系，若在其上则相交，否则不想交（在本问题中如此判断合适）。所以对于线段与多边形，若线段与多边形任一条边相交，则线段与多边形相交，否则进而判断线段中点与多边形的位置关系（由于多边形凸，在这种情况下，中点与

8、多边形的位置关系和线段与多边形关系等价），若在外，则在外，反之易得。上述算法可以判定某条边的合法性，若不合法，则距离为无穷大，否则为欧氏距离。在上述算法的支撑下，可以求得所有障碍多边形顶点及起始点互连而成的所有距离。再借助Dijkstra算法生成我方7个岛礁的最短路径。Dijkstra算法的思想是若从某点A到连通图上其他点最短路径生成是递增顺序的的，则到下一点X待生成的最短路径或者是A-X，或者经过已生成最短路径的点。下面是结果呈现，数据在表格中给出。图 2 6边形近似图 3 12边形近似图 4南沙所有最短路径 图 5永暑礁到东门礁三、测量中沙、西沙部分岛礁到南沙7个岛礁的距离由于中沙、西沙基

9、本被我方控制，且距离南沙岛礁较远，因此，中沙、西沙到南沙部分岛礁的距离可以用欧氏距离代替，下面分别测量了从中建岛、盘石屿、中南暗沙、中沙大环礁、黄岩岛到美济礁、渚碧礁、永暑礁的距离，用于连通整个赋值图。然后利用Floyd算法生成最短路径完全图Dist2323和路径中继前驱Path2323,其中，Distij为岛礁i与岛礁j的最短距离，它或是欧氏距离，或是ESPO算法生成的距离，或是经过中继获得的最短距离。Pathij表是从岛礁i到岛礁j最短路径中的中级前驱，如果没有中继则为i。借助Path矩阵可以在生成TSP路径后提供详细路径信息。图 5用于连通南沙与西沙、中沙的路径1、 利用蚁群算法求解di

10、st矩阵的最优TSP解蚁群算法是求解TSP问题的经典算法，是一种元启发式搜索算法，借助群体智能来实现组合优化。但是参数的选取并不容易，算法陷入局部最优的可能性很大，影响全局搜索能力。几个重要的参数包括信息素重要程度指标alpha、可见度（距离倒数）重要程度指标beta、信息素挥发因子dispersefactor、蚂蚁数量antnum和蚂蚁圈最大信息素Q。其中alpha=0时，算法退化为贪心算法，beta=0时，算法退化为随机搜索。在实验过程中试验了大量参数，发现蚂蚁数量过多越容易陷入局部最优。最终选取alpha=0.8，beta=2，antnum=34，dispersefactor=0.5，并

11、引入dispersefactor的陷入局部最优的自适应调节机制即dispersefactor<=(dispersefactor*limit)0.5，该迭代修正将最终使dispersefactor趋于limit不至于修正过度，取limit为0.4。下面是最优解，为3235.49千米。所有图片都在附件中给出。 结论 本文完满的解决了既定问题，具有良好的适用性。 参考文献 王万良 人工智能及其应用（第三版）群智能算法及其应用 Blog：线与多边形位置关系 附录(源程序）1、 Floyd算法 path=zeros(23,23); dist=zeros(23,23); for i=1:23 for

12、 j=1:23 dist(i,j)=original(i,j); if dist(i,j)<inf path(i,j)=i; else path(i,j)=-1; end if i=j dist(i,j)=0; end end end for k=1:23 for i=1:23 for j=1:23 if dist(i,j)>dist(i,k) + dist(k,j) dist(i,j) = dist(i,k) + dist(k,j); path(i,j) = path(k,j); end end end End二:ESPO算法含Dijkstra算法、几何判定算法、可视化（运行平台

13、需要支持opencv）#include<iostream>#include<opencv2opencv.hpp>#define pi 3.1415926#define inf 9999999999using namespace std;cv:Mat mapaa;class pointpublic:double x, y;string name;point()x = y = 0;point(double x, double y)this->x = x;this->y = y;void setname(string name)this->name = na

14、me;class Linepublic:point start, end;double length;Line()Line(point start, point end)this->start = start;this->end = end;length = sqrt(start.x - end.x)*(start.x - end.x) + (start.y - end.y) *(start.y - end.y);void setlength()length = sqrt(start.x - end.x)*(start.x - end.x) + (start.y - end.y)

15、*(start.y - end.y);template <int edgenum>class Polygon/whose methods appropriate for sample convex polygonspublic:int edgen;point vertaxedgenum;Line edgeedgenum;Polygon()void setPolygon(point vertax)edgen = edgenum;for (int i = 0; i < edgenum; i+)this->vertaxi = vertaxi;for (int i = 0; i

16、 < edgenum; i+)edgei.start = vertaxi;edgei.end = vertax(i + 1) % edgen;edgei.setlength();int intersect(Line line1, Line line2)subtract(line1.end, line1.start);subtract(line2.end, line1.start);subtract(line2.start, line1.start);subtract(line1.start, line1.start);double theta = -arctan(line1.end.x,

17、 line1.end.y);rotate(line1.end, theta);rotate(line2.start, theta);rotate(line2.end, theta);if (line2.start.y*line2.end.y < 0)double x0 = line2.start.x - line2.start.y*(line2.end.x - line2.start.x) / (line2.end.y - line2.start.y);if (x0 - line1.end.x)*x0 < 0)return 1;if (fabs(x0) < 0.000001

18、| fabs(x0 - line1.end.x) < 0.000001)return 2;else if (fabs(line2.start.y*line2.end.y) < 0.000001)return 3;elsereturn 0;point intermediate(point p1, point p2)point temp;temp.x = (p1.x + p2.x) / 2;temp.y = (p1.y + p2.y) / 2;return temp;bool interconnect(Line line)for (int i = 0; i < edgen; i+

19、)if (line.length < 0.00001)return false;if (intersect(line, edgei) = 1)return true;if (inside(intermediate(line.start, line.end) = false)return false;elsereturn true;bool inside(point p)int cnt = 0;int lock = 0;point auxi;auxi.x = -500;auxi.y = -300;Line ray(p, auxi);for (int i = 0; i < edgen;

20、 i+)Line temp = edgei;int value = intersect(ray, temp);if (fabs(temp.length) < 0.00001)return false;if (fabs(p.y - temp.start.y)*(temp.end.x - temp.start.x) - (temp.end.y - temp.start.y)*(p.x - temp.start.x) <= 0.00001)/on th edgereturn false;if (value = 1|value=3)cnt+;else if (value = 2)if (l

21、ock = 0)cnt+;lock = 1;elselock = 0;if (cnt % 2 = 1)return true;elsereturn false;double arctan(double x, double y)if (x = 0 && y = 0)cout << "error" << endl;elseif (x > 0)return atan(y / x);else if (x = 0)if (y > 0)return pi / 2;else if (y < 0)return -pi / 2;els

22、e if (y >= 0)return atan(y / x) + pi;else if (y < 0)return atan(y / x) - pi;void rotate(point &p1, double theta)point temp = p1;p1.x = temp.x*cos(theta) - temp.y*sin(theta);p1.y = temp.y*cos(theta) + temp.x*sin(theta);void subtract(point& p1, point p2)p1.x -= p2.x;p1.y -= p2.y;template

23、<int edgenum>class enemypublic:point center;double radius;int edgen;string name;point vertaxedgenum;void setenemy(point enemy,double radius,string name)this->center = enemy;this->radius = radius;this->name = name;edgen = edgenum;for (int i = 0; i < edgen; i+)vertaxi.x = center.x +

24、radius*cos(2*pi / edgenum* i);vertaxi.y = center.y + radius*sin(2*pi / edgenum * i);enemy()void setradius(double radius)this->radius = radius;for (int i = 0; i < edgen; i+)vertaxi.x = center.x + radius*cos(2*pi / edgenum * i);vertaxi.y = center.y + radius*sin(2*pi / edgenum * i);void drawarea(

25、)for (int i = 0; i < edgen; i+)line(mapaa, cv:Point(vertaxi.x, vertaxi.y), cv:Point(vertax(i + 1) % edgen.x, vertax(i + 1) % edgen.y), cv:Scalar(052, 25, 65);double euclidean_metric(point p1, point p2)return sqrt(p1.x - p2.x)*(p1.x - p2.x) +(p1.y - p2.y)*(p1.y - p2.y);voidespo(enemy<12>*ene

26、, int enenum, point china,int chinanum)Polygon<12> *poly = new Polygon<12>enenum;for (int i = 0; i < enenum; i+)polyi.setPolygon(enei.vertax);enei.drawarea();int vertaxnum = enenum * ene0.edgen + chinanum;point*vertax = new pointvertaxnum;for (int i = 0; i < chinanum; i+)vertaxi =

27、chinai;for (int i = 0; i < enenum * ene0.edgen; i+)vertaxi + chinanum = enei / ene0.edgen.vertaxi % ene0.edgen;double *dist = new double*vertaxnum;for (int i = 0; i < vertaxnum; i+)disti = new doublevertaxnum;for (int i = 0; i < vertaxnum; system("cls"), cout<<"已完成"

28、;<<i<<"/"<<vertaxnum,i+)for (int j = 0; j < vertaxnum;j+)if (j = i)distij = inf;continue;Line temp(vertaxi, vertaxj);int lock = 0;for (int t = 0; t < enenum; t+)if (erconnect(temp)=true)lock = 1;break;if (lock = 1)distij = distji = inf;else distij = distji

29、= euclidean_metric(vertaxi, vertaxj);for (int cnc = 0; cnc < chinanum; cnc+)int *generated = new intvertaxnum;for (int i = 0; i < vertaxnum; i+)generatedi = 0;int *pre = new intvertaxnum;double *mindist = new doublevertaxnum;for (int i = 0; i < vertaxnum; i+)generatedi = 0;if (distcnci <

30、 inf)prei = cnc;elseprei = -1;mindisti = distcnci;mindistcnc = 0;for (int i = 0; i < vertaxnum; i+)double min = inf;int temp = -1;for (int j = 0; j < vertaxnum; j+)if (generatedj = 0 && mindistj <= min)min = mindistj;temp = j;generatedtemp = 1;for (int j = 0; j <vertaxnum; j+)if

31、(!generatedj && mindistj>mindisttemp + disttempj)mindistj = mindisttemp + disttempj;prej = temp;for (int i = cnc+1; i < chinanum; i+)/cv:Mat mappp = mapaa.clone();int pathnum = 0;int *path = new intvertaxnum;cout << "从 " << << " 到 " &

32、lt;< << "距离是" << mindisti/800*397.7 <<"千米"<< endl;int index = prei;pathpathnum = i;while (index != -1)pathnum+;pathpathnum = index;index = preindex;for (int ij = 0; ij < pathnum; ij+)circle(mapaa, cv:Point(int)vertaxpathij.x, (int)vertaxp

33、athij.y), 5, cv:Scalar(i*10, i, i*12);line(mapaa, cv:Point(int)vertaxpathij.x, (int)vertaxpathij.y), cv:Point(int)vertaxpath(ij + 1) % (pathnum + 1).x, (int)vertaxpath(ij + 1) % (pathnum + 1).y), cv:Scalar(i*10, 45, 255);ostringstream os;os << mindisti / 800 * 397.7;string str= os.str();/cv:im

34、show("wall1", mapaa);cv:imwrite("完整.bmp", mapaa);/getchar();int main()double length = 394.8, width = 267.7, baselong = 112.28, baselati = 11.20, radius = 44.11, mapl = 800, maph = 536;string obstacle = "中业岛","大现礁","小现礁","福禄礁","鸿庥岛"

35、;,"九章群礁","安达礁","舶兰礁","三角礁","毕生礁","双黄沙洲","库归礁","仁爱礁","仙娥礁","马欢岛","五方礁","火艾礁","蒙自礁" ;double enelong = 114.29,113.85,113.99,113.62,114.36,114.43,114.71,114.58,115.34,113.67,

36、114.42,114.61,115.86,115.45,115.86,115.72,114.90,114.78 ;double enelati = 11.03,10.05,9.94,10.18,10.15,9.82,10.30,10.37,10.18,8.94,10.64,10.81,9.7,9.35,10.66,10.47,10.83,11.09 ;string china = "赤瓜礁","美济礁","永暑礁","华阳礁","南薰礁","渚碧礁","东门礁&qu

37、ot; ;double cnlong = 114.28,115.54,112.96,112.83,114.22,114.08,114.57 ;double cnlati = 9.68,9.86,9.59,8.85,10.17,10.88,9.92 ;point*chin = new point7;enemy<12> enearea18;mapaa = cv:imread("南沙诸岛.png");for (int i = 0; i < 18; i+)eneareai.setenemy(point(enelongi - baselong) * 111 / le

38、ngth*mapl, (baselati - enelatii) * 111 / width*maph),radius,obstaclei);enearea4.setradius(20);enearea5.setradius(20);for (int i = 0; i < 7; i+)chini.x = (cnlongi - baselong) * 111 / length*mapl;chini.y = (baselati - cnlatii) * 111 / width*maph;chini.setname(chinai);espo(enearea, 18, chin, 7);三、蚁群

39、算法含可视化#include<iostream>#include<stdlib.h>#include<time.h>#include<opencv2opencv.hpp>using namespace std;#define size 23int Q = 1000;double Q0 = 100;#define antnum 34double alpha = 0.8;double beta = 2;#define cir 10000double dispersefactor = 0.5;#define original 100cv:Mat bac

40、kground;double distsizesize;double pheromonesizesize;int presizesize;string island = "永兴岛","西沙洲","银砾滩", "东岛", "甘泉岛", "滨湄滩", "湛涵滩", "玉琢礁", "华光礁", "盘石屿", "中建岛", "浪花礁", "黄岩岛"

41、;, "中沙环礁北", "中沙环礁南","中南暗沙", "永暑礁", "华阳礁", "赤瓜礁", "美济礁", "东门礁", "南薰礁","渚碧礁" ;double longitude = 112.35,112.25,112.24,112.73,111.68,112.47,112.56,112.03,111.71,111.79,111.2,112.52,117.76,114.71,114.06,1

42、15.44,112.96,112.83,114.28,115.54,114.57,114.22,114.08 ;double latitude = 16.84,16.98,16.76,16.63,16.51,16.33,16.4,16.33,16.22,16.03,15.77,16.03,15.15,16.07,15.51,13.94,9.59,8.85,9.68,9.86,9.92,10.17,10.88 ;class pointpublic:double x, y;string name;point()void setinfor(string name, double x, double

43、y)this->name = name;this->x = x;this->y = y;point islansize;class Antpublic:double possibilitysize;int currentpos;int time;double pathlength;int pathsize;Ant()void init()currentpos = rand() % size;time = 0;pathlength = 0;path0 = currentpos;for (int i = 0; i < size; i+)possibilityi = 0;bo

44、ol haspassed(int nextpos)for (int i = 0; i < time; i+)if (nextpos = pathi)return true;return false;int nextpos()double pm = 0;for (int i = 0; i < size; i+)if (haspassed(i) = false)pm += possibilityi=pow(pheromonecurrentposi, alpha)*pow(distcurrentposi, -beta);else possibilityi = 0;int posssize

45、;int sum = 0;for (int i = 0; i < size; i+)possibilityi /= pm;sum+=possi = (int)(possibilityi*cir);int x;int randnum = rand() % sum;sum = 0;for (int i = 0; i < size; i+)sum += possi;if (randnum < sum)return i;int i = 0;while (haspassed(i) = true)i+;return i;void move()time+;int next= nextpos

46、();pathtime = next;pathlength += distcurrentposnext;currentpos = next;if (time = size - 1)pathlength += distnextpath0;time = 0;bool hasgonethrough(int x, int y)for (int i = 0; i < time; i+)if (pathi = x&&path(i + 1) % size = y)return true;return false;void getpath(int x, int y, int&nu

47、m, int path)num = 0;int index = y;path0 = y;donum+;index = prexindex;pathnum =index; while (index!= x);for (int i = 0; i <=num/ 2; i+)int temp = pathi;pathi = pathnum - i;pathnum - i = temp;class antColonypublic:Ant antantnum;double optimal;int x;int n0;int opathsize;int cnt;int flag;antColony()v

48、oid init()flag = 0;cnt = 0;optimal = inf;n0 = 5;void findoptimal()double min = optimal;for (int i = 0; i < antnum; i+)if (anti.pathlength < min)min = anti.pathlength;x = i;for (int j = 0; j < size; j+)opathj = antx.pathj;flag = 0;if (fabs(min - optimal) < 0.000001)flag+;optimal = min;voi

49、d run(int n)for (int i = 0; i < n; i+)round(1000);show();draw();void show()cout << optimal << "千米" << " "for (int i = 0; i < size; i+)cout << << " "cout << endl;cout << alpha << " " << beta << " " << dispersefactor << endl;void update()for (int i = 0; i < size; i+)for (int j = 0; j < size; j+)double delta = 0;for (int k = 0; k < antnum; k+)if (antk.hasgone