用鲍威尔法求解最优化问题
本帖最后由 chlh_jd 于 2013-10-8 17:44 编辑大家好！
很久以前就在想：用LISP写一个求解最优化的通用函数，因为LISP语言本身就比较适合；基于对数学的弱弱了解，写了个粗糙的函数和大家讨论，希望您能多参与这个话题，感谢在先！
;;(Powell val con arg arg_limit eps)
(defun Powell(val con arg
arg_limit eps / ss_powell_arg0
ss_powell_1dmin ss_powell_main_function
ss_powell_a0 ss_powell_dfm_lst ss_powell_dfm
ss_powell_x0 ss_powell_f0 ss_powell_x1
ss_powell_f1 ss_powell_x2 ss_powell_f2)
;; --------------------------------------------------------;;
;; Use Powell method to solve solve optimization problems
;;
;; Args:
;; val -- Expressions calculation function
;; con -- Restriction function
;; arg -- A string list of the Parameters's name , such as '( "x1" "x2" ...)
;; arg_limit -- parameters's limits , like '( (-1e14 1e14) (0 1e14) ...)
;; eps -- accuracy
;;
;; by GSLS(SS) 2013.10.6
;;---------------------------------------------------------;;
(defun ss_powell_1dmin(x -xm x0 +xm / f0 fa fb -fm +fm)
;; One-dimensional search
;; This function has not the Powell method .
;; How to change it connect with Gradient direction ?
(set (read x) x0)
(setq f0 (if (con)
(val)
1e99))
(set (read x) -xm)
(setq fa (if (con)
(val)
1e99))
(set (read x) (* 0.5 (+ x0 -xm)))
(setq -fm (if (con)
(val)
1e99))
(set (read x) +xm)
(setq fb (if (con)
(val)
1e99))
(set (read x) (* 0.5 (+ x0 +xm)))
(setq +fm (if (con)
(val)
1e99))
(cond
((or (equal x0 -xm eps) (equal x0 +xm eps))
x0)
((< fa -fm f0 +fm fb)
-xm)
((> fa -fm f0 +fm fb)
+xm)
((and (> fa -fm f0) (< f0 +fm fb))
(ss_powell_1dmin x (*0.5 (+ x0 -xm)) x0 (* 0.5 (+ x0 +xm))))
((and (> fa -fm) (< -fm f0 +fm fb))
(ss_powell_1dmin x -xm (* 0.5 (+ x0 -xm)) x0))
((and (> fb +fm) (> fa -fm f0 +fm))
(ss_powell_1dmin x x0 (* 0.5 (+ x0 +xm)) +xm))
((and (equal x0 -xm (sqrt eps)) (equal x0 +xm (sqrt eps)))
(/ (+ x0 -xm +xm) 3.))
))
(setq ss_powell_arg0 arg)
(if arg_limit
(mapcar
(function
(lambda (x y)
(cond
((or (not y) (apply (function or) y))
(set (read (strcat x "_limit")) (list -1e14 1e14)))
((not (car y))
(set (read (strcat x "_limit")) (list -1e14 (cadr y))))
((not (cadr y))
(set (read (strcat x "_limit")) (list (car y) 1e14)))
(t (set (read (strcat x "_limit")) y))
)))
arg
arg_limit)
(foreach x arg
(set (read (strcat x "_limit")) (list -1e14 1e14)))
) ;_ x1_limit
(setq ss_powell_x0
(mapcar
(function
(lambda (x)
(*
0.5
(+
(car (eval (read (strcat x "_limit"))))
(cadr (eval (read (strcat x "_limit"))))))))
arg))
(defun ss_powell_main_function(/ ss_powell_i)
(ss-set arg ss_powell_x0)
(setq ss_powell_f0
(if (con)
(val)
1e99))
(setq ss_powell_a0
(mapcar (function (lambda (x)
0.0))
arg)
ss_powell_dfm_lst
nil
ss_powell_i 0)
(mapcar
(function
(lambda (x/
ss_powell_a1_f0 ss_powell_a1_x0
ss_powell_a1_l ss_powell_a1_-xm
ss_powell_a1_+xm ss_powell_a1_ai
ss_powell_a1_f1 )
(setq ss_powell_a1_f0
(if(con)(val)1e99))
(setq ss_powell_a1_x0 (eval (read x)))
(setq ss_powell_a1_l (eval (read (strcat x "_limit"))))
(setq ss_powell_a1_-xm (car ss_powell_a1_l)
ss_powell_a1_+xm (cadr ss_powell_a1_l))
(setq ss_powell_a1_ai
(- (ss_powell_1dmin
x
ss_powell_a1_-xm
ss_powell_a1_x0
ss_powell_a1_+xm)
ss_powell_a1_x0))
(setq ss_powell_a0 (ch-lst ss_powell_a1_ai ss_powell_i ss_powell_a0))
(ss-set arg
(mapcar (function +) ss_powell_x0 ss_powell_a0))
(setq ss_powell_a1_f1
(if (con)
(val)
1e99))
(setq ss_powell_dfm_lst
(cons (- ss_powell_a1_f0 ss_powell_a1_f1) ss_powell_dfm_lst))
(setq ss_powell_i (1+ ss_powell_i))
))
arg)
(setq ss_powell_dfm_lst (reverse ss_powell_dfm_lst))
(setq ss_powell_x1
(mapcar (function +) ss_powell_x0 ss_powell_a0))
(setq ss_powell_x2
(mapcar (function (lambda (x y)
(- (* 2. x) y)))
ss_powell_x1
ss_powell_x0))
(ss-set arg ss_powell_x1)
(setq ss_powell_f1
(if (con)
(val)
1e99))
(ss-set arg ss_powell_x2)
(setq ss_powell_f2
(if (con)
(val)
1e99)))
(ss_powell_main_function)
(while (not (equal ss_powell_x1 ss_powell_x0 eps))
(setq ss_powell_dfm (apply (function max) ss_powell_dfm_lst))
(if
(or (>= ss_powell_f2 ss_powell_f0)
(>= (* (+ ss_powell_f0 ss_powell_f2 (* -2. ss_powell_f1))
(- ss_powell_f0 ss_powell_f1 ss_powell_dfm))
(* 0.5
ss_powell_dfm
(expt (- ss_powell_f0 ss_powell_f2) 2))))
(progn
(if (< ss_powell_f1 ss_powell_f2)
(setq ss_powell_x0 ss_powell_x1)
(setq ss_powell_x0 ss_powell_x2))
;_ (setq ss_powell_dfm_i (vl-sort-i ss_powell_dfm_lst (function >)))
)
(progn;_this is not
(setq ss_powell_i
(car
(vl-sort-i ss_powell_dfm_lst (function >))))
(setq ss_powell_a0
((lambda(/ i)
(setq i -1)
(mapcar
(function (lambda (y)
(setq i (1+ i))
(if (= i ss_powell_i)
y
(* -2.y))))
ss_powell_a0))))
(setq ss_powell_x0 (mapcar (function +) ss_powell_x0 ss_powell_a0))
;_(setq ss_powell_dfm_i (shift(vl-sort-i ss_powell_dfm_lst (function >))))
)
;_(setq ss_powell_x0 ss_powell_x1)
)
;|
(setq arg (mapcar (function (lambda (i)
(nth i arg)))
ss_powell_dfm_i))
(setq ss_powell_x0 (mapcar (function (lambda (i)
(nth i ss_powell_x0)))
ss_powell_dfm_i))|;
(ss_powell_main_function))
(ss-set arg ss_powell_x1)
(setq ss_powell_f1
(if (con)
(val)
1e99))
(list (mapcar (function (lambda (x) (eval(read x)))) ss_powell_arg0) ss_powell_f1))
;;; Used function
(defun ss-set(ss_set_str_lst ss_set_lst / ss_set_i)
(setq ss_set_i 0)
(repeat (length ss_set_lst)
(set (read (nth ss_set_i ss_set_str_lst))
(nth ss_set_i ss_set_lst))
(setq ss_set_i (1+ ss_set_i))))
;;-----------------------------
(defun ch-lst(new i lst / j len fst mid)
(if (/= (type i) (quote list))
(cond
((not (listp lst))
lst)
((minusp i) lst)
((> i (setq len (length lst))) lst)
((> i (/ len 2))
(reverse (ch-lst new (1- (- 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 (vl-list* (cadddr lst)
(caddr lst)
(cadr lst)
(car lst)
fst)
lst (cddddr lst)))
(reverse fst))
(list new)
(cdr lst))))
(progn
(setq j (cadr i)
i (car i))
(if j
(progn
(setqmid (nth i lst)
mid (ch-lst new j mid))
(ch-lst mid i lst))
(ch-lst new i lst)))))
应用例子1：求10(x1+x2-5)^2+(x1-x2)^2的极小点和极小值，理论解： ((2.5 2.5) 0)
(defun val()
;; 10(x1+x2-5)^2+(x1-x2)^2
(+ (* 10. (expt (+ x1 x2 -5.) 2)) (expt (- x1 x2) 2)))
(defun con () T)
(setq arg(list "x1" "x2")
arg_limitnil
eps1e-6)
(Powell val con arg arg_limit eps);;->((2.5 2.5) 1.32989e-011)
应用例子2：求60-10x1-4x2+x1^2+x2^2-x1x2的极小点和极小值，理论解：((8 6) 8)
(defun val()
;;60-10x1-4x2+x1^2+x2^2-x1x2
(+ 60. (* -10. x1) (* -4. x2) (* x1 x1) (* x2 x2) (* -1. x1 x2)))
(defun con () T)
(setq arg(list "x1" "x2")
arg_limitnil
eps1e-6)
(Powell val con arg arg_limit eps);;->((8.0 6.0) 8.0) 这么深的东西，实际工作中从来没有用到过，今后可能也用不到。帮你顶一下 一个字
伪原码
哦，编辑的时候，前面的被我删除了 (defun val()
;;x1^2+2x2^2-2x1x2-4x1
(+ (* x1 x1) (* 2. x2 x2) (* -4. x1)(* -2. x1 x2)))
(defun con () T)
(setq arg(list "x1" "x2")
arg_limitnil
eps1e-6)
(Powell val con arg arg_limit eps)
为什么以上运行不出来呢？ Powell用lisp实现很特别。赞一个
楼主函数可以运行。改变函数后迭代 有时导致acad死机
主程序迭代是否增加约束条件，，否则收敛慢的函数计算时间过长，容易死机
页:
[1]