点集最短回路（TSP问题）讨论
本帖最后由 chlh_jd 于 2012-10-9 23:09 编辑点集最短回路实际上是TSP问题的简化版，属于N-P难问题，目前针对这个问题国内已有很多研究。
这里有篇TSP问题算法综述 ：http://www.joces.org.cn/CN/article/downloadArticleFile.do?ttachType=PDF&id=9165
维基百科主要列了：遗传算法、模拟退火、群集智能（词条上没有解释，个人猜想是一种结合人工智能把点簇归到一起最佳出路放置在点簇凸包上的一种组合算法、算法关键应该是在点集归簇，或者说其实质是模拟退火归簇与遗传算法的结合）、局部搜索和多主体优化系统。http://zh.wikipedia.org/wiki/TSP
美国 Richard Johnsonbaugh, Marcus Schaefer 著作 《大学算法教程》里面提及了一些NP完全的基本算法有：蛮干法、随机方法、近似方法、参数化方法、启发式搜索 等。
二叉树描述法：有见过采用最优二叉树算法求解的TSP的，一般结合遗传算法作前期最优二叉树数据构造。
http://baike.soso.com/v6426781.htm
启发式搜索法：基于直观或经验构造的算法，现代优化算法中常常仅是其中的一步，单独用来求解TSP不太合适；应用实例还是有的：《一种改进的TSP问题启发式算法》
TSP 问题的脂肪计算复杂性与启发式算法设计
最邻近法：一种偏向于解决特定几何问题的弱算法，一般用于遗传算法的基因变异，也有人用来求解简单TSP：
单元划分法： 其结果与凸包塌陷方式接近，但凸包塌陷方式得全局性要比它好得多。
禁忌搜索法：常常容易陷入局部最优，需要设置一定的跳出关卡来扩散，以得到最大化的全局最优，
常用于解决排序问题（如车间作业排序、图案着色等）；因为代码较容易实现，也有很多人用来求解TSP，
举例： http://file.lw23.com/b/b4/b4a/b4a11e29-de82-4ccd-9ca1-3efc9c15c4a8.pdf
http://file.lw23.com/a/aa/aa2/aa230e4c-0bd3-4e8d-8eb2-258b576faa34.pdf
http://file.lw23.com/1/12/121/12126da9-4356-4a13-a16a-0c37f760bcf7.pdf
http://166.111.121.20:9080/mathjournal/XTGC803/xtgc803004.caj.pdf
神经网络法：清华大学数学系列教材的《现代优化计算方法》中有介绍一种叫弹性变形算法，其实质就是神经网络法，过程与凸包塌陷的逆过程相似。
百度文库 《SOFM神经网络最近插入法混合算法在TSP问题中应用研究》
《Hopfield神经网络在TSP问题中的应用 》
模拟退火法：基于概率分布的局部搜索算法的推广，常用于解决下料问题；因为建模相对简单，也不乏用来解TSP的实例；
http://www.vckbase.com/index.php/wv/1196
旅行商问题(TS P)的改进模拟退火算法
蚁群优化算法：常用于医疗数据挖掘；由于需要大量的信息积累，应用于TSP显得效率较低，常常需要结合其他手段来较好地实现。
百度文库有篇讲稿 《蚁群算法》
遗传算法：目前公认解决TSP等NP-Hard问题的最佳算法，但人们常常不满足于他误差和效率，也因此衍生出非常多的其他算法应用到基因选取、种群优化、交叉变异等组合遗传算法。
编程过程较为复杂（尤其是采用VLISP来实现），目前的编码多数采用二进制，也有采用十进制、十六进制的。前些年要找些类似的源代码学习很困难，现在就很多了。
有关研究非常多《基于混合遗传算法的TSP问题优化研究 》
组合算法： 模拟退火+遗传算法
基于模拟退火的遗传优化算法在TSP问题中的应用
个人认为采用组合遗传算法解决这个问题比较合适，采用VLISP实现最好加个类库处理（如矩形分布、点簇分布、离散分布...）。
早年收集了些遗传算法的文章也上传在这里
俄罗斯的VLisp顶级高手ElpanovEvgeniy曾在Theswamp上提出来讨论，现在看来他已经采用遗传算法解决了这个问题。
http://www.theswamp.org/index.php?topic=30434.75
这里引用Evgeniy提供的2个表供大家讨论用
chlh_jd 发表于 2013-1-6 03:38
@zdqwy19
我提供的TSP函数带有一定的遗传特性，并不一定能得到最优解，尤其对于规则点阵；该算法可 ...
你的程序每次运行结果都是一样的，没有看出遗传算法的随机性啊 mahuan1279 发表于 2017-11-13 21:21
测试图为什么会出现交叉现象？
用遗传算法得到一个近优解。
本帖最后由 chlh_jd 于 2012-10-9 11:51 编辑
贴出用 ElpanovEvgeniy 的方法改进 纯VLISP 版本 ， 比较适用于不规则离散点集
;;;------------------------TSP------------------------------------------------------------;;;
;;;---------------------------------------------------------------------------------------;;;
(defun c:test (/ foo f2 ptl lst l n i d0 l0 l1 d1)
;;by GSLS(SS)
;;refer ElpanovEvgeniy's method fromhttp://www.theswamp.org/index.php?topic=30434.75
;;2012-8-10
(defun foo (l / D D0 D1)
(setq l0 (mapcar (function list) (cons (last l) l) l)) ;_setq
;_defun
(setq d0 (get-closedpolygon-length l))
(while
(> d0
(progn
(foreach a l0
(setq d (get-closedpolygon-length l))
(setq l1 (vl-remove (car a) (vl-remove (cadr a) l)))
(setq l1 (f1 (car a) l1))
(setq l1 (f1 (cadr a) l1))
(if (> d
(setq d1 (get-closedpolygon-length l1))
)
(setq d d1
l l1
) ;_setq
) ;_if
(setq l1 (vl-remove (car a) (vl-remove (cadr a) l)))
(setq l1 (f1 (cadr a) l1))
(setq l1 (f1 (car a) l1))
(if (> d
(setq d1 (get-closedpolygon-length l1))
)
(setq d d1
l l1
)
)
)
d
) ;_progn
) ;_<
(setq d0 d)
) ;_while
(setq d (get-closedpolygon-length l))
l
)
(defun f1 (a l)
(ins-lst a (get-closest-i l a) l)
)
(defun f2 (lst)
(mapcar (function (lambda (p0 p p1 / a)
(setq a (- (angle p p0) (angle p p1)))
(if (< a (- pi))
(abs (+ a pi pi))
(if (> a pi)
(abs (- a pi pi))
(abs a)
)
)
)
)
(cons (last lst) lst)
lst
(reverse (cons (car lst) (reverse (cdr lst))))
)
)
(setq ptl (my-getpt)
ptl (mapcar (function (lambda (p) (list (car p) (cadr p)))) ptl)
)
(setq t1 (getvar "MilliSecs"))
(setq lst (Graham-scan ptl))
(foreach a lst
(setq ptl (vl-remove a ptl))
)
(while (and (> (length ptl) 2) (setq l (Graham-scan ptl)))
(foreach p l
(setq ptl (vl-remove p ptl))
(setq n (get-minadddist-i lst p))
(setq lst (ins-lst p n lst))
)
)
(if ptl
(foreach p ptl
(setq n (get-minadddist-i lst p))
(setq lst (ins-lst p n lst))
)
)
(setq lst (foo lst))
(setq l (f2 lst))
(setq i0
l0 lst
n(length lst)
d0 (get-closedpolygon-length lst)
)
(foreach a l
(if (and (< a _pi3) (= (setq p (nth i lst)) (nth i l0)))
(progn
(if (= i 0)
(setq p0 (last lst))
(setq p0 (nth (1- i) lst))
)
(if (= i (1- n))
(setq p1 (car lst))
(setq p1 (nth (1+ i) lst))
)
(setq m (list (list p0 p1 p)
(list p1 p p0)
(list p1 p0 p)
(list p p0 p1)
(list p p1 p0)
)
)
(setq l1
(car (vl-sort (mapcar (function (lambda (x)
(ch-para-lst x i lst)
)
)
m
)
(function (lambda (e1 e2)
(< (get-closedpolygon-length e1)
(get-closedpolygon-length e2)
)
)
)
)
)
)
(setq d1 (get-closedpolygon-length l1))
(if (< d1 d0)
(setq d0d1
lst l1
)
)
)
)
(setq i (1+ i))
)
(setq l (f2 lst))
(setq i0
l0 lst
d0 (get-closedpolygon-length lst)
)
(foreach a l
(if (and (< a _pi2) (setq p (nth i l0)))
(progn
(setq l1 (f1 p (vl-remove p lst)))
(setq d1 (get-closedpolygon-length l1))
(if (< d1 d0)
(setq d0d1
lst l1
)
)
)
)
(setq i (1+ i))
)
(entmake
(append (list '(0 . "LWPOLYLINE")
'(100 . "AcDbEntity")
'(8 . "temp")
'(62 . 1)
'(100 . "AcDbPolyline")
(cons 90 (length lst))
'(70 . 1)
)
(mapcar (function (lambda (p) (cons 10 p))) lst)
)
)
(setq t2 (getvar "MilliSecs"))
(princ (strcat "\nTSP Length :" (rtos d0 2 0) "."))
(princ (strcat "\nUse Time :" (rtos (- t2 t1) 2 0) "ms."))
(princ)
)
;;;Use Funtions
;;;--------------------------------------------------------------
;; Convex hull of pts , Graham scan method
;; by Highflybird
(defun Graham-scan (ptl / hPs rPs PsY Pt0 sPs P Q)
(if (< (length ptl) 4) ;3点以下
ptl ;是本集合
(progn
(setq rPs (mapcar (function (lambda (x)
(if (= (length x) 3)
(cdr x) x)))
(mapcar 'reverse ptl));_点表的X和Y交换
PsY (mapcar 'cadr ptl) ;_点表的Y值的表
Pt0 (reverse (assoc (apply 'min PsY) rPs)) ;_最下面的点
sPs (sort-ad ptl Pt0) ;_按角度距离排序点集
hPs (list (caddr sPs) (cadr sPs) Pt0) ;_开始的三点
)
(foreach n (cdddr sPs) ;从第4点开始
(setq hPs (cons n hPs) ;把Pi加入到凸集
P (cadr hPs) ;Pi-1
Q (caddr hPs) ;Pi-2
)
(while (and q (> (det n P Q) -1e-6)) ;如果左转
(setq hPs (cons n (cddr hPs)) ;删除Pi-1点
P (cadr hPs) ;得到新的Pi-1点
Q (caddr hPs) ;得到新的Pi-2点
)))
hPs ;返回凸集
))
)
;;;以最下面的点为基点，按照角度和距离分类点集
(defun sort-ad (pl pt)
(vl-sort pl
(function (lambda (e1 e2 / an1 an2)
(setq an1 (angle pt e1)
an2 (angle pt e2))
(if (equal an1 an2 1e-6);_这里降低误差，以适应工程需求
(< (distance pt e1) (distance pt e2))
(< an1 an2)
))))
)
;;定义三点的行列式,即三点之倍面积
(defun det (p1 p2 p3)
(- (* (- (car p2) (car p1)) (- (cadr p3) (cadr p1)))
(* (- (car p3) (car p1)) (- (cadr p2) (cadr p1)))
))
;;;
;;;------------------------
(defun my-getpt (/ ss i en l)
(setq ss (ssget '((0 . "point"))))
(setq i -1)
(while (setq en (ssname ss (setq i (1+ i))))
(setq l (cons (cdr (assoc 10 (entget en))) l))
)
)
;;;------------------------
;;;
;;(ins-lst 10 5 '(1 2 3 4 5))
;; i 为新插入元素的位置
(defun ins-lst (new i lst / len fst)
(cond
((minusp i)
lst
)
((> i (setq len (length lst)))
lst
)
((> i (/ len 2))
(reverse (ins-lst new (- len i) (reverse lst)))
)
(t
(append
(progn
(setq fst nil)
(repeat (rem i 4)
(setq fst (cons (car lst) fst)
lst (cdr lst)
)
)
(repeat (/ i 4)
(setq fst (cons (cadddr lst)
(cons (caddr lst)
(cons
(cadr lst)
(cons
(car lst)
fst
)
)
)
)
lst (cddddr lst)
)
)
(reverse fst)
)
(list new)
lst
)
)
)
)
;;;------------------------
;;
;;(ch-para-lst '(7 8 9) 3 '(1 2 3 4 5))
(defun ch-para-lst (para i lst / len fst)
(setq len (length lst))
(cond
((minusp i)
lst
)
((> i (1- len))
lst
)
((= i 0)
(cons (cadr para)
(cons (caddr para)
(reverse (cons (car para) (cdr (reverse (cddr lst)))))
)
)
)
((= i (1- len))
(reverse
(append (cdr (reverse para))
(cddr (reverse (cons (last para) (cdr lst))))
)
)
)
((> i (/ len 2))
(reverse
(ch-para-lst (reverse para) (- len i 1) (reverse lst))
)
)
(t
(append
(progn
(setq fst nil)
(repeat (rem i 4)
(setq fst (cons (car lst) fst)
lst (cdr lst)
)
)
(repeat (/ i 4)
(setq fst (cons (cadddr lst)
(cons (caddr lst)
(cons
(cadr lst)
(cons
(car lst)
fst
)
)
)
)
lst (cddddr lst)
)
)
(reverse
(cons (caddr para)
(cons (cadr para) (cons (car para) (cdr fst)))
)
)
)
(cdr lst)
)
)
)
)
;;;------------------------
;;
(defun get-minadddist-i (lst p)
(car
(vl-sort-i
(mapcar (function (lambda (p1 p2)
(- (+ (distance p p1) (distance p p2))
(distance p1 p2)
)
)
)
(cons (last lst) lst)
lst
)
'<
)
)
)
;;;------------------------
(defun get-closest-i (lst p)
(car
(vl-sort-i
(mapcar
(function
(lambda (p1 p2 / pt d d1 d2)
(setq pt (inters p
(polar p (+ (/ pi 2.) (angle p1 p2)) 1.)
p1
p2
nil
)
d(distance p1 p2)
d1 (distance p p1)
d2 (distance p p2)
)
(if pt
(if (equal (+ (distance pt p1) (distance pt p2)) d 1e-8)
(distance p pt)
d2
)
1e99
)
)
)
(cons (last lst) lst)
lst
)
'<
)
)
)
;;;------------------------
;;
(defun get-closedpolygon-length (l)
(apply (function +)
(mapcar (function (lambda (p1 p2)
(distance p1 p2)
)
)
(cons (last l) l)
l
)
)
) 比较深奥，不过很有用 强大！学习！谢谢！ 收藏 深奥,太深奥 好！效果如何，前一个和opendcl有冲突，新的不知道，回去试一下。 回去试了一下1300个点约2-3分钟，感觉有点慢，不过图形蛮理想。 @zdqwy19
感谢您的测试，希望能有更好的建议。