用鲍威尔法求解最优化问题

chlh_jd · 发表于 2013-10-6 16:12

本帖最后由 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
(setq mid (nth i lst)
mid (ch-lst new j mid))
(ch-lst mid i lst))
(ch-lst new i lst)))))

;;(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 
    (setq  mid (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_limit nil
eps 1e-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_limit nil
eps 1e-6)
(Powell val con arg arg_limit eps);;->((8.0 6.0) 8.0)

自贡黄明儒 · 发表于 2013-10-8 08:34

这么深的东西，实际工作中从来没有用到过，今后可能也用不到。帮你顶一下

crtrccrt · 发表于 2013-10-8 12:56

一个字
伪原码

chlh_jd · 发表于 2013-10-8 17:36

哦，编辑的时候，前面的被我删除了

mahuan1279 · 发表于 2017-11-11 10:41

(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_limit  nil
   eps  1e-6)
(Powell val con arg arg_limit eps)

为什么以上运行不出来呢？

landsat99 · 发表于 2022-5-27 14:42

Powell用lisp实现很特别。赞一个

landsat99 · 发表于 2022-5-27 15:09

楼主函数可以运行。改变函数后迭代有时导致acad死机

主程序迭代是否增加约束条件，，否则收敛慢的函数计算时间过长，容易死机

账号		自动登录	找回密码
密码			注册

[自我挑战] 用鲍威尔法求解最优化问题

评分

点评