西电数模大作业

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

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

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

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

5、间、 非南沙和南沙之间其距离认定 为其平面欧几里得距离，南沙群岛各岛礁距离利用 ESPO 算法生成，由于南沙距离非南沙相对于其 尺度很大，故近似比较准确。方案形成过程南海岛礁经纬度采集和部分岛礁之间的测距借助百度地图及其提供的坐标拾取系统和测距工具，获取了我方23个岛礁、他方18个岛礁的经纬度信息和用于连通南沙我方部分岛礁（障碍外部）与中南暗沙、大环礁、中建岛、盘石屿等岛礁的 距离。并利用MATLAB计算了我方非南沙各岛礁之间的欧氏距离，从而生成了一部分最短路径完全图。由于表格较大，不在此呈现，详见表格。z二少T盘石帖浪花礁中沙环礪北1寻中廳細f |SsS苗左卫1昆丄比：f * #图1南海部分

6、岛礁华阳礁南沙群岛有障碍区域我方岛礁最短路径的形成由假设知，各敏感区域是由以外国侵占岛礁几何中心为圆心半径为12海里的区域，为凸集。若可行区域的边界是光滑曲面。则最短路径必由下列弧组成，它们或者是空间中的自然最短曲线，或者是 可行区域的边界弧。而且，组成最短路径的各段弧在连接点处必定相切。上述定理由J.W.Craggs在1973年证明。由于光滑曲线在计算机中不便于存储和处理，本文将敏感区域近似为其内接正多边形（本文中最终使用以 12边形近似的数据），并基于上述定理近似推断在障碍区两点之间的最短路径 必有以下线段组成，他们或是障碍多边形的边，或是障碍多边形各顶点及起始点之间的合法线段组 成，其中

7、合法线段为不与任何障碍多边形相交或在其外部的线段。因此，在上述分析的支持下，试 图用计算机解决这个平面简单凸多边形障碍区域最短路径问题首先要解决如何判断平面上点、线段 与多边形的位置关系。这是一个基本的计算几何问题，对于线段和多边形位置关系的判定，首先要 解决线段与线段之间、点与多边形之间的位置关系判定方法。对于点与多边形，可以利用射线法， 即由该点为端点任意引一条射线，若射线与多边形各边相交奇数次，则在内部，其他在外部（在边 上也认为在外部）；对于线段与线段之间，可利用正交变换（便于处理）将其中一条变换到X轴且一段与原点重合，进而判断另一条线段与X轴的交点与第一条线段的关系，若在其上则相交，

8、否则不想交（在本问题中如此判断合适）。所以对于线段与多边形，若线段与多边形任一条边相交，则 线段与多边形相交，否则进而判断线段中点与多边形的位置关系（由于多边形凸，在这种情况下， 中点与多边形的位置关系和线段与多边形关系等价），若在外，则在外，反之易得。上述算法可以 判定某条边的合法性，若不合法，则距离为无穷大，否则为欧氏距离。在上述算法的支撑下，可以 求得所有障碍多边形顶点及起始点互连而成的所有距离。再借助Dijkstra算法生成我方7个岛礁的最短路径。Dijkstra算法的思想是若从某点A到连通图上其他点最短路径生成是递增顺序的的，则到下一点X待生成的最短路径或者是A-X，或者经过已生成最

9、短路径的点。下面是结果呈现，数据在表格中给出。图26边形近似图312边形近似图4南沙所有最短路径目賂3图5永暑礁到东门礁三、测量中沙、西沙部分岛礁到南沙7个岛礁的距离由于中沙、西沙基本被我方控制，且距离南沙岛礁较远，因此，中沙、西沙到南沙部分岛礁的距离 可以用欧氏距离代替，下面分别测量了从中建岛、盘石屿、中南暗沙、中沙大环礁、黄岩岛到美济 礁、渚碧礁、永暑礁的距离，用于连通整个赋值图。然后利用Floyd算法生成最短路径完全图Dist2323和路径中继前驱 Path2323,其中，Distij为岛礁i与岛礁j的最短距离，它或是欧氏 距离，或是ESPO算法生成的距离，或是经过中继获得的最短距离。P

10、athij表是从岛礁i到岛礁j最短路径中的中级前驱，如果没有中继则为i。借助Path矩阵可以在生成 TSP路径后提供详细路径信息。II占起点总长：匸公里I 总W； 565,2里Ls*：屈舄长；&忆吃里总长已能亂吆里起点图5用于连通南沙与西沙、中沙的路径一、利用蚁群算法求解dist矩阵的最优TSP解蚁群算法是求解 TSP问题的经典算法，是一种元启发式搜索算法，借助群体智能来实现组合优化。 但是参数的选取并不容易，算法陷入局部最优的可能性很大，影响全局搜索能力。几个重要的参数 包括信息素重要程度指标alpha、可见度(距离倒数)重要程度指标beta、信息素挥发因子dispersefactor、蚂蚁

11、数量antnum和蚂蚁圈最大信息素Q。其中alpha=O时，算法退化为贪心算法，beta=O时，算法退化为随机搜索。在实验过程中试验了大量参数，发现蚂蚁数量过多越容易陷入局 部最优。最终选取 alpha=0.8, beta=2， antnum=34， dispersefactor=0.5，并引入 dispersefactor 的陷入 局部最优的自适应调节机制即dispersefactor: o oho” 自耳芦出於二tjp Q w -结论本文完满的解决了既定问题，具有良好的适用性。参考文献王万良 人工智能及其应用（第三版）群智能算法及其应用 Blog ：线与多边形位置关系附录 （源程序）、 F

12、loyd 算法path=zeros(23,23);dist=zeros(23,23);for i=1:23for j=1:23dist(i,j)=original(i,j); 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);endendendEnd二 :ESPO 算法含 Dijkstra 算法、几何判定算法、可视化(运行平台需要支持opencv)#include#include#define pi 3.1415926#define inf 9999999999usi

13、ng 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 = name;class Linepublic:point start, end;double length;Line() Line(point start, point end)this-start = start;this-end = end;lengt

14、h = 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) *(start.y - end.y); ;template class Polygon/whose methods appropriate for sample convex polygonspublic:int edgen;point vertaxedgenu

15、m;Line edgeedgenum;Polygon()void setPolygon(point vertax)edgen = edgenum;for (int i = 0; i vertaxi = vertaxi;for (int i = 0; i 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.e

16、nd, line1.start); subtract(line2.start, line1.start);subtract(line1.start, line1.start);double theta = -arctan(line1.end.x, line1.end.y);rotate(line1.end, theta);rotate(line2.start, theta);rotate(line2.end, theta);if (line2.start.y*line2.end.y 0)/ (line2.end.y -double x0 = line2.start.x - line2.star

17、t.y*(line2.end.x - line2.start.x) line2.start.y);if (x0 - line1.end.x)*x0 0)return 1;if (fabs(x0) 0.000001 | 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;te

18、mp.y = (p1.y + p2.y) / 2;return temp;bool interconnect(Line line)for (int i = 0; i edgen; i+)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

19、;point auxi; auxi.x = -500;auxi.y = -300;Line ray(p, auxi);for (int i = 0; i edgen; 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)

20、/on th edgereturn false;if (value = 1|value=3) cnt+;else if (value = 2) if (lock = 0) cnt+; lock = 1; else lock = 0;if (cnt % 2 = 1) return true;else return false;double arctan(double x, double y) if (x = 0 & y = 0) cout error 0) return atan(y / x);else if (x = 0)if (y 0) return pi / 2;else if (y =

21、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 class enemypublic:point

22、 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 radius = radius;for (int i = 0; i edgen; i+) vertaxi.x = center.x + radius*cos(2*pi /

23、edgenum * i); vertaxi.y = center.y + radius*sin(2*pi / edgenum * i); void drawarea() 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(p

24、1.x - p2.x)*(p1.x - p2.x) +(p1.y - p2.y)*(p1.y - p2.y); void espo(enemy*ene, int enenum, point china,int chinanum)Polygon *poly = new Polygonenenum;for (int i = 0; i enenum; i+) polyi.setPolygon(enei.vertax); enei.drawarea();int vertaxnum = enenum * ene0.edgen + chinanum; point*vertax = new pointver

25、taxnum;for (int i = 0; i chinanum; i+)vertaxi = 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), cou

26、t 已完成 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 = euclidean_metric(vertaxi, vertaxj);for

27、 (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 inf)prei = cnc;elseprei = -1;mindisti = distcnci;mindistcnc

28、 = 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 mindisttemp + disttempj)mindistj = mindisttemp + disttempj;prej = temp;for (int i = cnc+1; i chinanum; i

29、+)/cv:Mat mappp = mapaa.clone();int pathnum = 0;int *path = new intvertaxnum;cout 从 到 距离是 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

30、)vertaxpathij.x, (int)vertaxpathij.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

31、 str= os.str();/cv:imshow(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 = 中业岛,大现礁,小现礁,福禄礁,鸿庥岛 ,九章群礁 ,安达礁,舶兰礁 ,三角礁 ,毕生礁 ,双黄沙洲 ,库归礁,仁爱礁 ,仙娥礁,马欢岛,五方礁 ,火艾

32、礁,蒙自礁 ;double enelong = 114.29,113.85,113.99,113.62,114.36,114.43,114.71,114.58,115.34,113.67,114.42,114.61,115.86,115.45,11 5.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 Jstring china = 赤瓜礁 ,美济礁,永暑礁,

33、华阳礁 ,南薰礁,渚碧礁,东门礁 ; 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 enearea18;mapaa = cv:imread( 南沙诸岛 .png);for (int i = 0; i 18; i+)eneareai.setenemy(point(enelongi - baselong) * 111 / length*mapl,

34、 (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);三、蚁群算法含可视化#inclu

35、de#include#include#includeusing 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 inf 10.#define original 100cv:Mat background;double distsizesize;double pheromonesizesize;int presizesize

36、;string island = 永兴岛,西沙洲,银砾滩, 东岛, 甘泉岛 , 滨湄滩, 湛涵滩 , 玉琢礁 , 华光 礁, 盘石屿, 中建岛, 浪花礁, 黄岩岛, 中沙环礁北 , 中沙环礁南 ,中南暗沙 , 永暑礁, 华 阳礁, 赤瓜礁 , 美济礁, 东门礁 , 南薰礁 ,渚碧礁 ;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,115.44,112.96,112.83,114.28,115.5

37、4,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 y)this-name = name; this-x = x;

38、 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;bool haspassed(int nextpos)for (int i = 0; i ti

39、me; i+)if (nextpos = pathi)return true;return false;int nextpos()double pm = 0;for (int i = 0; i size; i+)if (haspassed(i) = false)alpha)*pow(distcurrentposi,pm += possibilityi=pow(pheromonecurrentposi, -beta);else possibilityi = 0;int posssize;int sum = 0;for (int i = 0; i size; i+) possibilityi /=

40、 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();pathtime = next;pathlength += distcurrentposnext; currentpos = ne

41、xt;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&num, int path) num = 0; int index = y; path0 = y;donum+;index = prexindex; pathnum =i

42、ndex; 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()void init()flag = 0; cnt = 0; optimal = inf;n0 = 5;void findoptimal()double min =

43、 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;void 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 alp