//Logo Image
作者:徐業良 (1995-09-15);推薦:徐業良 (2001-06-19);最近更新:徐業良(2004-10-29)
附註:本文為元智大學機械工程研究所最佳化設計課程教材,僅限於教學上學生個人使用,原書初版由宏明書局印行。

多變數函數內部最小值搜尋

在討論完單變數函數最小值搜尋之後,這一章要擴展到多變數函數內部最小值的搜尋,也就是在最佳化設計模型中沒有限制條件,或者所有限制條件均非有效的狀況。

先前提到,一般迭代形式求最小值的數值演算法,都包括了初始值給定、迭代定義、以及中止要件三個部分,這一章中我們可以看到,求多變數函數內部最小值的各種演算法,也只是迭代定義部分有不同的變化。而大部分求多變數函數內部最小值的演算法都是屬於所謂「直線搜尋(line search)」法,不同演算法的迭代定義,只是在定義不同的搜尋方向,並且在所定義的搜尋方向求得該直線上目標函數的最小值,作為下一個迭代點。因此在每一次迭代中的直線搜尋,事實上都是單變數函數求最小值的問題,先前所討論單變數函數求最小值的方法都可以應用上。

這一章對多變數函數內部最小值搜尋演算法的討論,將先由工程上最直覺的「函數值比較法(function comparison method)」出發,接下來便將繼續討論直線搜尋法,和其所涵蓋的各種二次、一次、零次的演算法。

1.     函數值比較法

考慮一個簡單的雙彈簧系統平衡問題,如圖1所示[Vanderplaats, 1984]

1. 簡單的雙彈簧系統平衡問題

彈簧受到P1的水平力及P2垂直力後,系統彈力位能PE可以表示成如下式:

        (1)

其中為點A受力後的水平與垂直位移。這個簡單的雙彈簧系統受力之平衡狀態,即是系統位能為最小值時的位移狀態,因此求此系統平衡時點A的水平與垂直位移,即是在求PE之最小值。

要求PE最小值可能有許多種方法,當然我們可以用解析形式的解法,令,兩式聯立求解。但是在目標函數形式較複雜,或者設計變數的數目較多時,解析形式的解法相當困難。例如把式(1)作一次微分,所得式子將相當複雜,聯立求解析形式解幾乎不可能,因此這裡的重點還是放在迭代形式的數值方法

要求式(1)的最小值,工程上最直覺、最直接的方法可能是不斷調整設計變數的值,比較在不同變數值組合之下,何者目標函數值最小,這也就是所謂的目標函數值比較法。至於如何調整設計變數的值,則有三種可能的作法,這裡首先討論所謂「隨機搜尋法(random search)

如以下虛擬程式所示,隨機搜尋法中設計者首先對各個設計變數給予一個合理的上界與下界,在每一次迭代中,便對每一變數產生一個介於零與一之間的亂數,而此次迭代中的變數值即設為

                                                                                        (2)

這裡k是迭代次數,則必定在此變數的上界與下界之間。得到此次迭代的設計點後,接下來檢查在此設計點上目標函數值是否較原先儲存的最小值為小,如果是,便必須更新此函數最小值,否則繼續迭代下去,直到超過最大容許迭代次數為止。

   algorithm 1 Random Search

        begin

             input: , , i = 1, ..., n; ;

             for i = 1 to n do

             begin

                 generate random number ;

                 ;

             end;

             ; ;

             ;

             while  do

             begin

                 for i = 1 to n do

                 begin

                      generate random number ;

                      ;

                 end;

                 if  then

                 begin

                      ;

                      ;

                 end;

                 ;

             end;

             output , ;

        end.

◇作業1

Matlab寫一隨機搜尋法的電腦程式,以第一節中的雙彈簧系統位能問題作測試,為設計變數設定合理的上下界,列印出其經過第100200300、……、1000次迭代之後的最佳設計變數值與最佳目標函數值。重新執行一次,再列印出其經過第100200300、……、1000次迭代之後的最佳設計變數值與最佳目標函數值。本章作業中寫程式的目的,僅在實際了解各個演算法的特性,以及可能發生的數值問題,因此程式只要能處理兩變數問題即可。◇

隨機搜尋法顯然相當沒有效率,它需要非常多次目標函數值的運算,而大部分的函數值運算卻都沒有用而被浪費掉了,且不管迭代多少次,隨機搜尋法都完全沒有任何保證可以收斂至理論最小值。但是這個方法在許多應用場合還是有用的,它是一個「零次方法」,當目標函數是不平滑函數(non-smooth function),有許多不連續點,或者一次微分不存在時,一般直線搜尋演算法幾乎都無法應用,隨機搜尋法卻仍然有效。此外原理簡單、容易寫成電腦程式、且運算上也比較沒有數值問題,也是隨機搜尋法的優點。

對於一般平滑函數,自然絕對沒有必要使用隨機搜尋法,第二種調整設計變數值的方式,則是應用現有設計點上的數值資訊,來決定如何修改現有設計。平滑函數都有幾何上的「連貫性(coherence)」,從現有設計點上的數值資料,如一次、二次微分,往往便可判定目標函數在此設計點附近的整體趨勢,進而決定應該如何改變設計點以改進目標函數值。本章下一節將要介紹的直線搜尋數值方法,都是屬於這種運用設計點上的數值資訊,來決定下一個迭代點的方法。

但就工程最佳化問題來說,運用設計點上的數值資訊,來決定下一個迭代點的這類演算法最大的盲點,卻也正是其僅使用每一迭代設計點上的數值資訊。這類演算法假設我們所處理的問題都只是單純的數學函數,而忽略了設計者對此問題所具有工程上的經驗或知識,可能對最佳化的過程很有幫助,這是十分可惜的。

因此第三種調整設計變數值的作法,對於一般設計者、工程師可能是最直覺、最常用的方法,不是靠隨機搜尋猜測,也不僅只靠每一次迭代設計點上的數值資訊來決定搜尋方向,而是以設計者自己的經驗與對此系統的知識,配合前幾次改變設計點所得結果所顯現的趨勢,作一個綜合判斷來決定如何調整設計變數組合,改變現有的設計。當然對一般設計者來說,這種改變設計的方法通常近乎是「嘗試錯誤」,最大問題是很難如最佳化數值演算法一般,有系統地決定下一個迭代點,逐步找到最佳解,這也是在設計程序之中,我們需要使用最佳化數值演算法(或數值最佳化演算的概念)的重要原因。

事實上如何將設計者對某些特定工程問題的知識、經驗模擬下來,與最佳化的數值方法加以整合,使其成為有系統的最佳化演算法的一部份,是工程最佳化設計上值得研究的一個領域。

◇作業2

觀察作業1隨機搜尋的迭代歷程,是否可以得到一些「經驗法則」如何可以讓此演算法更“聰明”一些?敘述你的觀察,並嘗試將這些經驗法則加入隨機搜尋演算法中,與作業1比較看看,新的演算法搜尋效率是否提高了?◇

2.     直線搜尋法的觀念

搜尋函數最小值常常被比喻作在濃霧中下山,你的目的是要走到最低的山谷中,但是卻看不到整個山的形狀,只能看到你所站立位置附近的地形,這時候你應該採取怎麼樣的策略,才能走到山谷?

一個直覺而簡單的策略是,你環顧四周,選定一個下坡的方向,沿著這個方向一直走,直到又要開始上坡了,你便停下來,再次環顧四周,選定另一個下坡方向,朝著這個方向前進。如果你停下來,環顧四周,發現周圍全是上坡,沒有下坡方向,那麼你很有理由判定,你應該在某一個山谷了。這個簡單的濃霧下山策略,便是多變數函數內部最小值直線搜尋演算法的基本原則。

2.1 下降的搜尋方向

前一節中提到,一般平滑函數都有幾何上的連貫性,從現有設計點上的數值資料往往便可判定目標函數在此設計點上的整體趨勢,進而決定如何改變設計點以改進目標函數值。在直線搜尋演算法中,這些現有設計點上的數值資料,便是在用來決定搜尋方向。

以圖1中的雙彈簧系統為例,式(1)中的目標函數PE可以繪成如圖2的等高線圖。從圖中A點出發,首先我們必須要決定一個搜尋方向s。而這個搜尋方向s必須是“下坡”方向,也就是在這個方向上的微小擾動,必須能夠降低目標函數的值。很顯然的,在第k次迭代中的下降方向必須能夠滿足下式:

                                                            (3)

設計點上的下降方向可能有無限多個,但一個很明顯的下降方向是梯度向量的相反方向

       

            (4)

採用負的梯度向量(gradient)作為搜尋方向的演算法,就稱作「梯度法(gradient method)。由式(4)可知,在設計點,負梯度方向也是「最陡下降方向(steepest descent direction)」,所以梯度法也稱作「最陡下降法(steepest decent method)。至於其他演算法採用不同的搜尋方向,在下面兩節會進一步說明之。

2. 雙彈簧系統位能的等高線圖

2.2 搜尋方向上步長的決定

又如圖2所示,從點A出發,選擇負梯度向量為搜尋方向,接下來便要計算在這個搜尋方向上,什麼時候由下坡轉成上坡而應該要停下來,也就是要決定在這個搜尋方向上的「步長(step length)」大小。在最佳化演算法中,步長通常都是以符號a來表示,決定這個a值,事實上是一個單變數求最小值的問題。如圖2中,如果將此等高線圖由A-B截面垂直切下,可以得到如圖3所示曲線,在第k次迭代時,我們的目的便是在求出,步長是多少,f(a)有最小值,而,其中都是已知,所以在這個平面上的曲線是以a為變數的單變數函數。例如在點Aa=0,函數值為而決定最適當的步長,基本上就是在解圖3中的單變數函數最小值問題min. f(a),求出在最小點B的步長,此時函數值為

3. 決定步長是一個解單變數函數最小值的問題

因此求多變數函數內部最小值的直線搜尋演算法可以由以下虛擬程式表示,不同的直線搜尋演算法僅是虛擬程式中的副程式SEARCH_DIRECTION(x)有所不同。各種不同演算法對搜尋方向的定義,將在第三節中繼續討論之。

   algorithm 2 Line Search

        begin

             input: ;

             ;

             ;

             while TERMINATION(x, k) = false do

             begin

                 SEARCH_DIRECTION(x);

                 STEP_LENGTH(x, s);

                 ;

                 ;

             end;

             output , x;

        end.

◇作業3

將上述Line Search演算法繪成流程圖,逐步解釋其意義,特別注意解釋三個副程式的意義及其輸入、輸出。以Matlab建構這個主程式,並選擇一個中止要件,完成TERMINNATION副程式。◇

副程式STEP_LENGTH(x, s)基本上是在解圖3中的單變數函數最小值問題 min. f(a)。解這個問題可能有許多種方法,例如可以用類似先前提到的求函數零點的方法,如二分法或割線法,求f'(a) = 0之解。f'(a)可以下式表示:

                                                                                      (5)

其中是設計點對步長的變化率,因為搜尋方向為直線,所以(5)中求,從圖2中來看,也就是求搜尋方向和目標函數梯度向量互相垂直的點。實際計算上,則是求時的值。

其次可以用多項式近似法來求函數f(a)的最小值,首先計算f(a)上幾個點的函數值,用一條多項式曲線逼近,再計算此多項式曲線的最小值。當然這條近似多項式函數曲線的最小值,並不保證是原函數f(a)的最小值,且如不確定最小值的位置,求得的最小值是由近似曲線外差得到而非內差得到,這個近似最小值的準確性更是無法保證。

另外也可以用黃金切割法來解f(a)的最小值,然而如先前的討論,以黃金切割法解min. f(a)時必須輸入搜尋區間的兩個邊界值,滿足,也就是輸入滿足所謂“高-低-高”模式搜尋的不確定區間。這個搜尋區間起始點可設定在,另一邊界值則必須作進一步搜尋,以下便是可以找出滿足此“高-低-高”模式邊界值的虛擬程式。

   algorithm 3 Find Interval

        begin

             input: d;                                                       : d is a small step length;

             ;

             ;

             ;

             while  do

             begin

                 ;

                 ;

                 ;

                 ;

             end;

             output , ;

        end.

在下一節會討論到的牛頓法搜尋中,理想的步長a是等於1,因此前面虛擬程式中的起始輸入值d應該小於等於0.5。求出後,便可利用黃金切割法,求出f(a)的最小值,決定此次迭代移動步長。

4是雙彈簧系統位能問題中由(-5, -4)出發,梯度法迭代的過程圖。從圖中可以看到,越接近最小值時,梯度法每次迭代所移動的距離越短,收斂速度很慢,這也是梯度法實用上並不受歡迎的原因。

4. 雙彈簧系統位能問題梯度法的迭代過程

◇作業4

選擇任何一個方法,完成上題中STEP_LENGTH的副程式,特別注意其輸入、輸出。以負梯度方向為搜尋方向,完成SEARCH_DIRECTION副程式。以雙彈簧平衡問題為例,由不同起始點出發,可否產生如圖4之搜尋結果?◇

除了迭代次數之外,最佳化演算法效率也需要考慮每次迭代所需函數計算次數。由以上討論可知,在每一次迭代中要作精準的直線搜尋,需要計算目標函數值很多次,例如在黃金切割法中,首先為黃金切割法找尋滿足“高-低-高”模式的邊界值須作多次函數值計算,接下來用黃金切割法求f(a)的最小值又需要多次函數值運算。如何減少每次迭代中這些函數值計算的次數,實為求多變數函數內部最小值演算法效率高低的關鍵。下一節裡即在介紹所謂「不精準的直線搜尋(inexact line search)」,也就是在每次迭代中不須找到函數精準的最小值,以降低每次迭代中目標函數值計算的次數

2.3 不精準的直線搜尋

前一節中提到,要精準地搜尋函數f(a)的最小值需要多次函數值運算,然而沿著下降的搜尋方向行走,是否一定要在最低點才停下來呢?並不一定,事實上只要能夠找到步長a滿足「足夠下降(sufficient decrease)」的要求即可[Gill, Murray, and Wright, 1989]。下式中的Goldstein-Armijo原則,即是一個對第k次迭代的步長足夠下降要求的例子:

                               (6)

其中訂定後,在搜尋方向作直線搜尋時,每一次作新的函數值計算,都應代入式(7)作檢查,如果滿足式(6),這次迭代便可以停止而尋找新的搜尋方向亦即只要找到一個a值能夠滿足式(6),該次迭代便可停在此點。

Goldstein-Armijo原則其實有相當明確的幾何意義,首先的線性近似式可以表示如下,

                                                                   (7)

如圖5所示,點A的座標,直線為通過點A的切線(注意此處變數為,如式(5)),是的狀況,通過A點的水平線,則是的狀況,選定,式(6)即是限定步長的選擇必須使此次迭代目標函數值必須落在圖5兩條直線之間,直線限制步長不會太大,而直線則限制步長不會太小。

5. 不精準的直線搜尋

在不精準的直線搜尋中,的訂定顯然是十分重要的,如果範圍定得太小,每一次迭代中要找到滿足此要求的步長,所需函數值計算的次數勢必增多,但如果範圍定得太大,每一次迭代所得到的新設計點目標函數值的下降不夠,反而會造成總迭代次數增多。

1是雙彈簧系統位能問題,採用梯度法精準的直線搜尋與不精準的直線搜尋兩種處理方式,每次迭代的迭代點,迭代後的目標函數值,與所需目標函數值計算的次數。由此表中可以看出,在不精準的直線搜尋()中每次迭代所需函數值計算次數減少很多。一般來說不精準的直線搜尋每次迭代函數值的下降的速度可能較慢,但兩相比較起來,使用不精準的直線搜尋,計算效率上還是較精準的直線搜尋為佳。

1. 精準的直線搜尋與不精準的直線搜尋迭代過程比較

 

迭代次數

0

1

2

3

4

5

6

精準直線搜尋

迭代點

2.00, 2.00

5.56, 1.06

5.79, 2.36

6.74, 2.10

6.86, 3.04

7.45, 2.87

7.52, 3.55

函數值

-5.35

-29.17

-33.42

-318

-38.07

-39.34

-40.21

計算次數

22

22

25

22

29

22

25

不精準直線搜尋

迭代點

2.00, 2.00

6.83, 0.73

6.90, 3.63

7.04, 2.79

7.29, 3.75

7.43, 3.01

7.54, 3.78

函數值

-5.35

-24.93

-37.22

-38.56

-39.07

-39.62

-40.11

計算次數

1

3

3

3

1

1

2

◇作業5

修改上題中的程式:(1)每次迭代中,計數並輸出函數值計算次數,(2)STEP_LENGTH的副程式中增加不精準直線搜尋之選項。再以雙彈簧平衡問題為例,由不同起始點出發,可否產生如表1之比較表?◇

3.     二次與一次直線搜尋演算法

這一節裡將介紹其他幾個重要的求多變數函數最小值的直線搜尋演算法。如前所述,大部分直線搜尋演算法基本觀念都是相同的,只是在定義下降的搜尋方向時有所差異,因此這一節裡將只討論各個演算法對搜尋方向的定義。與先前討論單變數函數最小值搜尋演算法時類似,在定義搜尋方向時只需要目標函數函數值的演算法,可歸類為零次方法,在定義搜尋方向時需要目標函數梯度向量的演算法,可歸類為一次方法,需要目標函數的Hessian矩陣的演算法,則被歸類為二次方法

3.1 二次方法-牛頓法

先前討論單變數函數最小值搜尋時,曾經提到牛頓法,而此處求多變數函數最小值的二次方法,牛頓法也是最基本而重要的演算法。一個多變數函數的泰勒氏展開式可表示如下:

                                                   (8)

與單變數函數的牛頓法相同,這裡也是以二次函數式(8)去近似原函數f(x),而以此二次函數的最小點為下一次迭代點。解式(8)最小點的一次必要條件可得:

       

                                                                           (9)

因此牛頓法在第k次迭代的搜尋方向為

                                                                                   (10)

這個搜尋方向稱作「牛頓方向(Newton’s direction)。在這個推導中牛頓方向理想的步長a應該為1,但實際演算法進行時步長固定為1容易造成發散,每次迭代時通常還是依照前一節中討論的方法去計算步長,這種作法稱作「修正牛頓法(Modified Newton’s Method)。圖6是由(-5, -4)出發,以修正牛頓法解雙彈簧系統位能問題的迭代過程。牛頓法的搜尋方向雖然並非迭代點上下降最快的方向,但牛頓法整體的收斂速度遠較梯度法為快,如果原來的目標函數本身就是二次式,那麼式(8)便不只是近似式,而就是目標函數完整的泰勒展開式,牛頓法可在一次迭代便找到函數的最小點。

6. 以修正牛頓法解雙彈簧系統位能問題的迭代過程

牛頓法在實際應用上的一個重要的問題,是如果式(8)中的Hessian矩陣不是正定矩陣,解式(8)之一次必要條件所求出的搜尋方向可能不是一個下降的搜尋方向。

一個一般的下降的搜尋方向可用“Metric”表示成如下形式:

                                                                                                 (11)

其中矩陣便稱作Metric,例如,式(11)定義的方向為最陡下降方向,當,式(11)定義的方向為牛頓方向

如果是一個下降方向,則

                                                                             (12)

矩陣必須是一個正定矩陣。但是牛頓法中每次迭代的Hessian矩陣的反矩陣不一定是正定矩陣,此時牛頓方向不是一個下降的搜尋方向,碰到這種情況可以有以下幾種修正方法:

(1)   令此次迭代,也就是在此次迭代採用最陡下降方向而不用牛頓方向,這種方式稱作「交替梯度與牛頓法(Alternating gradient and Newton's method)

(2)   令此次迭代,嘗試幾個正數,直到為正定矩陣為止。

(3)   修正Cholesky分解,將矩陣一般Cholesky分解中的各分解矩陣修正,使得到,其中是一個元素均為正數或零的對角矩陣(diagonal matrix),而修正後的矩陣則為一個正定矩陣。

(4)   Hessian矩陣作「spectral 分解」,如的形式,其中的特徵向量組合成的矩陣,而則為對角線是的特徵值的對角矩陣。作此分解後,令中負的特徵值一律為正而成為,也就是將可能的上升方向顛倒符號,再計算

這幾種處理方法在此不作數學上的細節討論,但了解其中(2)(3)(4)項的基本想法,是在滿足下降方向的前提下,去建構一個儘量近似於Hessian矩陣的反矩陣的“Metric”,也就是盡量產生一個接近牛頓方向的搜尋方向,這個觀念與下一節將要討論的「準牛頓法(Quasi Newton’s method)」的觀念十分類似。

另外式(10)中,定義牛頓方向需要計算Hessian矩陣的反矩陣,如果Hessian矩陣是奇異矩陣(singular matrix)或狀況不佳的矩陣(ill-conditioned matrix),求反矩陣時一定會產生數值問題,應用上應該加以注意。

◇作業6

在上題程式SEARCH_DIRECTION副程式中,增加牛頓方向的選項,用修正牛頓法的方式解雙彈簧問題。注意此函數Hessian矩陣的解析解非常難求,試用有限差分法得到近似的矩陣,且每次迭代時輸出Hessian矩陣的特徵值。迭代過程中Hessian矩陣是否一直是正定矩陣?◇

3.2 一次方法-準牛頓法

本章第二節中提到梯度法的收斂非常慢,甚至一般研究最佳化數值方法理論的學者都建議不要用梯度法。而相對的牛頓法是二次近似,能對目標函數中非線性的特性作適當模擬,收斂速度就快得多。不過收斂速度快,相對付出的代價便是牛頓法中每次迭代都需要目標函數的Hessian矩陣,組成這個矩陣要計算目標函數對設計變數的二次微分,運算上是相當昂貴的。因此這一節中所要介紹的一次演算法「準牛頓法」,便是著眼於如何只用一次微分資料,來產生近似牛頓方向的搜尋方向

準牛頓法也是屬於所謂「變動metric(variable metric method)的一種,其每次迭代所用的metric ,都是由前一次迭代的metric加上一個「更新矩陣(updating matrix)」而得的:

        ,                                                                      (13)

其中是一個對稱的更新矩陣,定義如下

                     (14)

pyst分別定義為

                                                                                                   (15)

                                                                                     (16)

                                                                                                        (17)

                                                                                                     (18)

(13)(18)的數學細節也不在此討論,但從對這幾個式子的觀察可以了解準牛頓法的基本概念,也是嘗試用前幾次迭代的一次微分資訊去建構一個近似於原目標函數的Hessian矩陣之反矩陣的metric,來定義一個接近牛頓方向的搜尋方向

在式(14)q是一個自訂的參數,如果q= 0,此式便成為DFP更新(Davidon-Fleecher-Powell update)[Davidon, 1959; Fletcher and Powell, 1963];如果q= 1,此式則稱作BFGS更新(Broydon-Fletcher-Goldfarb-Shanno update)[Broydon, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970]。圖7是由(-5, -4)出發,以BFGS更新解雙彈簧系統位能問題的迭代過程,比較圖467,可以看出準牛頓法的搜尋方向較近似於牛頓法。

7. 雙彈簧系統位能問題BFGS更新的迭代過程

◇作業7

在上題程式SEARCH_DIRECTION副程式中,增加BFGS更新的選項,解雙彈簧問題。◇

3.3 共軛方向法

如果兩個搜尋方向滿足下式,則稱作共軛方向。

                                                                                                   (19)

其中H是目標函數的Hessian矩陣。

前面提到以牛頓方向為搜尋方向,對於一個二次目標函數,可以在一次搜尋找到其最小值。而數學上可以證明出,以共軛方向為搜尋方向,對於一個n變數的二次目標函數,可以用n或更少的共軛搜尋方向,找到其最小值(Zangwill, 1969)。共軛方向的數學意義在此不作深入探討,而共軛方向法的基本概念,就是每次迭代所定義的搜尋方向,要和前n-1次迭代的搜尋方向為共軛方向,如此以n個共軛方向搜尋,可以模擬一次牛頓方向搜尋,卻不須計算其Hessian矩陣

數學上可以證明,對於一個二次目標函數來說,如果第k次迭代的搜尋方向,要和前n-1次迭代的搜尋方向為共軛方向,應定義如下:[Gill, Murray, Wright, 1989]

        ,

        ,                                                                     (20)

(20)的搜尋方向僅是梯度法的搜尋方向加以修正,電腦程式應用上十分簡單。圖8是由(-5, -4)出發,以共軛方向法解雙彈簧系統位能問題的迭代過程,注意其搜尋方向與梯度法、牛頓法、準牛頓法的BFGS更新等方法的差異。

共軛方向法在實際應用上可能發生的問題,是以式(20)定義之共軛方向,有可能不是目標函數的下降方向,因此共軛方向求出之後,應先計算的值,判斷其是否函數的下降方向,如不是則應令,再讓整個演算法重新開始。

數學家Powell曾經發展出一種零次的共軛方向法[Powell, 1964],其原理是在一個n變數目標函數中,可以用n+1個零次搜尋方向來模擬一個式(20)中的共軛方向而不須目標函數的一次微分值。照此推算,需要n(n+1)個這樣的零次搜尋方向才能模擬一個牛頓方向。Powell的方法雖然不需要計算函數的一次微分,但當變數的數目增加時,其所需要的搜尋次數成平方增加,是其最主要的缺點。

8. 以共軛方向法解雙彈簧系統位能問題的迭代過程

◇作業8

在上題程式SEARCH_DIRECTION副程式中,增加共軛方向的選項,解雙彈簧問題。注意在共軛方向求出後,須先檢查其是否為下降方向,如否,當次迭代下降方向須先回復為負梯度方向。◇

◇作業9

進一步增加你程式中的搜尋方向,使其可以彈性選擇最陡下降方向、牛頓方向、BFGS修正、共軛方向等四種方向。

(a)    以函數測試你的程式,這個函數稱作Rosenbrock的香蕉函數,常用來測試各種最佳化數值演算法,如圖9所示,這個函數有一個香蕉形的山谷,最小點為(1, 1)T。用四種不同的搜尋方向作迭代,試由不同起始點出發,將迭代過程中的設計點繪於圖9等高線圖中。討論你的觀察。

(b)   以函數測試你的程式,這個函數稱如圖10所示,有4個區域最小值。用4種不同的搜尋方向作迭代,試由不同起始點出發,將迭代過程中的設計點繪於下頁等高線圖中。討論你的觀察。◇

9. 作業9(a)函數等高線圖

10. 作業9(b)函數等高線圖

參考資料

Broydon, C. G., 1970. “The Convergence of a Class of Double Rank Minimization Algorithms, part I and II.” J. Inst. Math. Appl., Vol. 6, pp. 76-90, 222-231.

Davidon, W. C., 1959. “Variable Metric Method for Minimization.” Argone National Laboratory, ANL-5990 Rev., University of Chicago.

Fletcher, R., 1970. “A New Approach to Variable Metric Algorithms.” Computer Journal, Vol. 13, pp. 317-322.

Fletcher, R. and Powell, M. J. D., 1963. “A Rapidly Convergent Method for Minimization.” Computer Journal, Vol. 6, No. 2, pp. 163-168.

Gill, P. E., Murray, W., and Wright, M. H., 1989. Practical Optimization, Eighth Printing, Academic Press, San Diego, California.

Goldfarb, D., 1970. “A Family of Variable Metric Methods Derived by Variational Means.” Math. Comput., Vol. 24, pp. 23-36.

Powell, M. J. D., 1964. “An Efficient Method for Finding the Minimum of a Function of Several Variables without Calculating Derivatives.” Computer Journal, Vol. 7, No. 4, pp. 155-162.

Shanno, D. F., 1970. “Conditioning of Quasi-Newton Methods for Function Minimization.” Math. Comput., Vol.24, pp. 647-656.

Vanderplaats, G. N., 1984. Numerical Optimization Techniques for Engineering Design With Applications, McGraw-Hill, New York.

Zangwill, W. I., 1969. Nonlinear Programming: A Unified Approach, Prentice-Hall, Englewood Cliffs, N.J.