//Logo Image
最近更新:徐業良(2005-10-31)﹔推薦:徐業良(2005-09-11)
附註:本文為元智大學機械所最佳化設計課程教材。

第五章 單變數函數最小值搜尋

建立了最佳化設計模型之後,接下來的工作便是要求這個最佳化設計模型的最佳解。前面章節中提到,當最佳化設計模型沒有限制條件、或限制條件都無效(inactive)時,最佳化設計模型的最佳解是一個內部最佳解,也就是求目標函數的最小值。

單變數函數最小值搜尋的演算法,是許多求多變數目標函數最小值的演算法發展的基礎,因此這裡對解最佳化設計模型的數值演算法的討論,便由單變數最小值搜尋的演算法出發。你也許認為求單變數函數最小值是一個非常簡單的問題,前面章節中提到,令目標函數一次微分為零,其解即為此目標函數的靜止點,這樣的求解方法是求目標函數最小值的「解析方法(analytical method)。然而在求目標函數一次微分為零的解時,往往還是需要用到迭代形式的數值方法,且在目標函數的一次微分很難求得、不存在、或不連續時,並無法由令一次微分為零的解析方法求目標函數最小值,而必須以迭代形式的「數值方法(numerical method),直接求目標函數的最小值。

這一章裡,首先介紹求函數最小值演算法迭代的觀念,第二小節則討論求函數「零點(zero)」,也就是求函數值為零之解的數值演算法,第三小節介紹求單變數函數最小值的數值演算法,最後第四節則在介紹多項式近似法求函數零點或最小值。

1.         求最佳解演算法迭代的概念

在工程設計過程中,設計者常常根據設計需求先做出第一代的初始設計,然後測試、評估這個初始設計,看看這個設計是否已經滿足所有設計需求,是否還有進一步改進的空間,接著設計者根據評估的結果、自己專業領域的知識、或者過去的設計經驗,產生第二版的改良設計,然後再重新測試、評估這個改良設計。如此持續作了好幾版的改良,直到設計者認為這個設計滿足了所有設計需求、已經足夠好,或者已經沒有額外時間和資源來發展新版的設計,設計者便決定停止。

求最佳解數值演算法的概念其實非常類似,一般迭代形式求最佳解演算法,都包括了初始值給定、迭代定義(iteration definition)、以及中止要件(termination condition)三個部分,這三個部分的關係,可大致由圖1簡單表示。在研究不同演算法時,只是這三個部分(特別是迭代定義的部分)的內容各有不同,而各個演算法的整體結構,則幾乎都是一致的。

1. 解最佳化設計模型演算法的基本迭代觀念

使用各種求最佳點數值演算法時,初始值的給定是很重要的。直覺地說,設計者給定的初始值離最佳點越近,迭代收斂的速率越快。除了收斂的速率之外,初始值的給定也要考慮到幾乎所有數值演算法的收斂性質都是屬於「區域收斂(local convergence)」,而非「全域收斂(global convergence),也就是說,如果設計者給定的初始點不在某個特定區間之內,或者沒有滿足某項特殊要求,則這個演算法可能根本不會收斂。

因此在實用上,設計者不應任意給定初始值,而應該考慮到所選擇的演算法有無特殊要求,同時就實際工程問題中各設計變數的物理意義考量,儘量給定一個合理且滿足所有限制條件的可行設計點,作為演算法的初始點

中止要件的部分,則可能有許多種不同的定義,端視演算法的性質,或某一類特殊問題的要求而定。在這裡首先要強調的是,“中止”和“收斂”並不相同,事實上在許多狀況下,演算法中止並不代表演算法已經收斂至理論最佳值。數值演算法在滿足中止要件時,有時表示現在的設計點在預設的容差範圍下、非常接近最佳點,或者是每次迭代之間的變化已經非常小了、繼續迭代下去也不會有明顯進步,或者是在可使用的計算資源下,已經無法或不值得繼續改進現在的設計。

具體來說,中止要件的定義通常有如下幾種形式:

(1)       ,其中L為先前所定義的拉氏函數,為在第k次迭代時的設計變數向量,是設計者預設的一個微小值。滿足這個中止要件時,設計點幾乎滿足KKT條件,非常接近最佳點。然而這個中止要件在一般工程最佳化問題上,並不十分實用,主要是因為在一般工程問題中,一次微分的計算昂貴且不精確,且除非值設得很大,否則一般很難找到滿足這個中止要件的設計點。

(2)       ,或,其中為在第k次迭代時的目標函數值,是設計者預設的兩個微小值。滿足這個中止要件,代表演算法在每次迭代的結果差異已經非常小,繼續迭代下去也不會有明顯進步。這個中止要件定義十分適合工程最佳化問題,但是要特別了解的是,滿足這個中止要件有時是因為演算法中移動的步伐長度太短,兩次迭代之間的移動速度太慢,而並不保證此時設計點已經接近理論最佳點。兩個微小值的設定,則應視設計者對目標函數及設計變數要求的精確度而訂。

(3)       迭代次數N。實用上來說,即使在流程圖中沒有標明出來,每一個數值演算法都會包括這一項中止要件,使得在最佳化設計模型有錯誤(例如具有單調性的設計變數無適當限制條件),或者演算法無法收斂的情況下,演算法不會無窮盡地迭代下去。通常在測試階段,設計者對所解的最佳化設計模型正確性並無絕對信心,或對所使用的演算法特性還無法完全掌握時,可把最大迭代次數N訂得低一些,以節省計算資源。演算法如果因為迭代次數超過最大迭代次數而中止的話,通常也代表這次運算是失敗的,所得結果很可能是發散的,甚至是不可行的設計點。

最後圖1中迭代定義的部分,則是每一種演算法的精華所在,各個演算法的不同通常也就是在迭代定義部分設計的不同。一般來說,迭代定義基本的形式與功能,都在根據前一次(或前幾次)迭代的設計點上的數值資料,來決定下一次迭代的設計點。而從迭代定義所需要的數值資料不同,我們又可以把演算法歸類為「零次方法(zero order method)」、「一次方法(first order method)」、或「二次方法(second order method)。例如接下來第二、第三小節中介紹的單變數函數零點或最小值搜尋方法,都只需要單變數函數的函數值,而不需要一次微分,且只要求函數值為連續,而不要求一次微分為連續,因此可歸類為零次方法。本章第四節介紹的多項式近似法求函數零點或最小值,則需要函數的一次或二次微分資料,可歸類為一次方法或者二次方法。

比較起來,零次方法每次迭代所需要的數值資料最少,但是一般來說收斂也較慢(所需迭代次數較多);一次或二次方法收斂所需迭代次數也許較少,但是每次迭代所需的數值資料卻也昂貴得多。

2.         求單變數函數的零點

本節首先介紹求單變數函數的零點的數值演算法,包括二分法和割線法。

2.1 二分法

中學時代的基礎數學應該就介紹過以下定理1的這個連續函數的基本觀念:

◇定理1

f(x)為一連續函數,若,且f(a)f(b)<0,則在區間[a, b]中必定存在一解x*,使得f(x*) = 0

這一章裡第一個要討論的求單變數函數零點的數值演算法,是所謂「二分法(bisection)」,也正是反覆應用定理1的概念。一個最基本的想法是,既然在區間[a, b]中,如果我們不斷減小區間[a, b],使得在多次迭代之後,那麼也就解出來了。因此二分法這個演算法中迭代定義的部分,就在如何有系統、有效率地減小這個「不確定區間(interval of uncertainty)

如圖2所示,f(a)f(b) < 0時,在區間[a, b]中有一零點,要有系統地縮小這個區間,一個直覺的選擇是取區間[a, b]的中點為下一個迭代點,也就是。而如圖2中,f(x1)f(b)<0,因此便取代成為a點,進入下一次迭代,同理,如果f(x1)f(a)<0,因此便取代成為下一次迭代中的b點。再如圖2中,取第二次迭代點,第三次迭代點,如此反覆迭代直到滿足中止要件,也就是f(x*)非常接近零為止。注意在這個方法中,每次迭代不確定區間便減小一半。

2. 二分法求函數零點

二分法可以寫成如下的虛擬程式,了解這個虛擬程式,應該有助於更嚴謹地了解這個演算法的邏輯性。注意在這個虛擬程式中,哪些部分是屬於初始值給定,哪些部分是中止要件、迭代定義。

algorithm 1 Bisection;

    begin

        input: a, b;                                       :

        ;

        ;

        while  do

            begin

                if  then

                else ;

                ;

                ;

            end;

                output , ;

    end.

2.2 割線法

另外一個求單變數函數零點的零次方法叫做「割線法(secant method),或者有的書稱此方法為「線性內差法(linear interpolation method)。如圖3所示,割線法顧名思義,是在函數上前兩次迭代點之間畫一條割線,來近似這個單變數函數,並求出這條割線與橫軸的交點,作為下一次迭代點。割線法的迭代定義如下式:

                                                                               (1)

其中fkfk-1便是第kk-1次迭代的函數值,xk為第k次迭代的設計點,xk+1是割線與橫軸的交點,也就是下一次迭代的新設計點。

3. 割線法求函數零點

割線法亦可寫成如下虛擬程式:

algorithm 2 Secant Method;

    begin

        input:

        ;

        while  do

            begin

               

                ;

            end;

                output , ;

    end.

割線法的「收斂速率(convergence rate)」較二分法為快(此處將不對演算法的收斂速率作正式、理論性的討論),但其主要問題是,如果xkxk-1在函數零點的同側,也就是xkxk-1形成的割線所求出與橫軸的交點xk+1是線性外差而非內差,如圖3中的點x4x2x3的線性外差,那麼割線法很可能會發散

割線法概念上屬於「序列近似法(sequential approximation method),不斷用兩點間的割線來近似原先的曲線。序列近似法的概念和二分法中不斷有系統地減小搜尋的不確定區間的想法,也正是下面介紹單變數函數最小值搜尋演算法兩種不同的體系。

◇作業1

在二分法和割線法的虛擬程式中,請標明出哪些部分是屬於初始值給定,哪些部分是中止要件、迭代定義。◇

◇作業2

Matlab寫一二分法與割線法的電腦程式,並用以下三個例題來作測試。何者收斂較快?如果有不收斂的情況,畫出函數圖形解釋其理由,並且說明如何可以補救這種不收斂的狀況。

(a)    ,起始點為13

(b)    ,起始點為-12

(c)    ,起始點為-101。◇

3.         求單變數函數的最小點

這一節將介紹兩個求單變數函數最小值的零次方法:「費邦那西搜尋(Fibonacci search)「黃金切割搜尋(golden section search)。這兩種方法迭代定義的基本概念,都和前一節中介紹的二分法十分類似,也就是在每次迭代過程中,不斷有系統地減小搜尋的不確定區間。

3.1 單模態函數

討論二分法搜尋函數零點時有一個基本的假設,是函數在起始搜尋區間中至少有一個零點。這裡對函數最小值的搜尋,則是假設所搜尋的函數在起始搜尋區間中是「單模態(unimodal),簡單的說,就是函數在起始搜尋區間中只有一個最小值。函數是單模態的正式定義如下[Gill, Marray, and Wright, 1989]

◇定義1

如果存在唯一的使得任何x1 < x2,滿足

        x2 < x*,則f(x1) > f(x2)

        x1 > x*,則f(x1) < f(x2)

則函數f(x)在區間[a, b]中是單模態。◇

注意這個定義的形式是單模態的充分條件。這個定義的內涵其實非常簡單,直接一點地說,這個定義告訴我們一個單模態函數只有一個最小點x*,在最小點x*的左邊,函數是單調遞減的,而在最小點x*的右邊,函數則成了單調遞增。定義1為這一節中搜尋單變數函數最小值的演算法提供了一個類似定理1、決定搜尋的不確定區間的準則。在定理1中搜尋單變數函數零點,我們要尋找兩個點函數值相乘為負值,函數在這兩點之間才至少存在一個零點;以下定理2是定義1的變形,搜尋單變數函數最小點,我們要尋找的是三個搜尋點其函數值呈“高-低-高”模式,這三個點構成之區間內才有最小點存在。

◇定理2

設函數在區間中是單模態,。若,則最小點(最小點在x1右邊才可能滿足定義1),反之若,最小點(最小點在x2左邊才可能滿足定義1)。◇

如圖4所示,在區間之間加入兩個搜尋點,如果b三點的函數值便呈一“高-低-高”模式,此時最小點,因此搜尋的不確定區間便由原先的縮減為。相反的,如果a三點的函數值便呈一“高-低-高”模式,此時最小點,因此搜尋的不確定區間便由原先的縮減為。了解縮小不確定區間的基本要求後,接下來要介紹的兩個演算法費邦那西搜尋和黃金切割搜尋,重點便在於如何定義區間中的比較點,以得到較快速的收斂。

4. 單變數函數最小值搜尋,尋找“高-低-高”模式。

3.2 費邦那西搜尋

從一個起始不確定區間出發,要如何決定中間搜尋點的位置,才能有效率地縮減這個不確定區間呢?例如在圖4中這兩個搜尋點要放在什麼位置,這次迭代完成後才能得到最小的不確定區間?直覺上來說,你也許會回答放在三等分的點,但這時不論最小點或者,下一次迭代的不確定區間都可縮減為原來的三分之二。然而根據費邦那西搜尋,這兩點應該都放在區間中點附近,相隔一微小距離d,如圖5(a),如此下一次迭代的不確定區間將縮減為原來的二分之一,是更有效率的分隔方式。

費邦那西搜尋的搜尋點位置配置,是基於一個特別的數列,叫做費邦那西數列,這個數列的通式如下:

        ,                                                                   (2)

由式(2)可知費邦那西數列頭幾個元素是1, 1, 2, 3, 5, 8, 13 .....。費邦那西搜尋中,給定N個搜尋點,便可依據費邦那西數列來決定這些搜尋點的位置。例如圖5(a)中,當N=2,即有兩個搜尋點時,便決定了這兩個搜尋點的位置分別在,兩點間隔一微小距離d,而根據這些搜尋點的高-低-高模式,最後不確定區間將縮減為

又如圖5(b)N=3,有三個搜尋點的情況,決定了前兩個搜尋點的位置為,不確定區間先縮減為2/3,第三個搜尋點放入之後,最後不確定區間縮減為,或

再看圖5(c)N=4,有四個搜尋點的情況,應該就能完全了解費邦那西搜尋的程序。便決定了頭兩個搜尋點的位置為2/53/5,不確定區間先縮減為3/5;第三個搜尋點放入後,不確定區域再縮減為;第四個搜尋點放入之後,最後不確定區間縮減為,或者

(a) N=2                                                   (b) N=3

(c) N=4

5. 費邦那西搜尋

依據這個程序,在有N個搜尋點的情況,費邦那西搜尋最後能將不確定區域縮減至1/FN。因此在決定變數的不確定區間需要的精確度後,可從費邦那西數列中查出需要多少個搜尋點才能夠滿足此精確度的需求,再依據費邦那西數列安排這些搜尋點的位置。

然而這種方式在電腦程式中勢必要儲存或計算產生一大串費邦那西數列,以決定搜尋點的數目或位置。一個自然的問題是,如果對搜尋點個數沒有限制,也就是說如果有無限多個搜尋點時,應該要如何安排這些搜尋點的位置?

下一節中介紹的黃金切割搜尋,是由費邦那西搜尋演變而來,但是搜尋點個數N可以趨近於無線大,也就是說不用事先決定搜尋點的數目,或最終不確定區間的大小,而縮小不確定區間的方式和費邦那西搜尋相同,是一個相當普遍的單變數函數搜尋法。

◇作業3

試繪出類似圖5N=5, 6費邦那西搜尋,不確定區間的切割模式。用N=5搜尋下列三個函數的最小值所在之不確定區間。

(a)    ,起始點為-13

(b)   ,起始點為-33

(c)    ,起始點為010

(d)   第一章籃球架設計樣本問題,,起始點為01。◇

3.3 黃金切割搜尋

從費邦那西搜尋的程序中應該可以體會到,要得到最有效率、最快速的收斂,基本的原則就是要儘量重複利用每一個搜尋點,而黃金切割搜尋中搜尋點位置的決定,也就是希望能夠儘量運用到前一次迭代的搜尋點。如圖6所示,當頭兩個搜尋點對稱地放在不確定區間[0,1]中,決定了“高-低-高”模式、把不確定區間縮減為之後,下一次迭代放入第三個搜尋點時,希望仍然能夠用到這兩個搜尋點,因此這兩個搜尋點的位置在縮減後的不確定區間中,希望仍然能保持和原先不確定區間相同的比例,對照圖6,即,也就是希望能夠滿足下式:

                                                                       (3)

                                                                                           (4)

        t為正數)

前述推導中1-t = 0.6180也就叫作「黃金切割數(golden section number)或者「黃金切割比(golden section ratio)。黃金切割搜尋中,不論搜尋點個數N為何,只要反覆把新的一個搜尋點放在原先不確定區間的t1-t處,便可把不確定區間縮減為原先的1-t倍。

6. 黃金切割搜尋

那麼黃金切割搜尋和費邦那西搜尋又有什麼關係呢?如果檢查費邦那西數構成的數列,其中是第k個費邦那西數,可以得到{1.000, 0.500, 0.667, 0.600, 0.625, 0.615, 0.619, ....}費邦那西數列的極限值正是黃金切割比1-t因此如果k趨近於無限大,也就是說作費邦那西搜尋時可以有無限多個搜尋點,每次新加入搜尋點位置,必定在原區間的t1-t處,而每次迭代後不確定區間縮小為原先的1-t倍。

除了用在單變數函數最小值搜尋之外,黃金切割比也有一些視覺、美學上的意義。例如藝術家達文西研究人體造型,便曾結論到一個“完美”的人體,頭頂到肚臍長度與肚臍至腳底長度的比例,就是黃金切割比。另外戲劇劇場中,最合適觀眾視覺感受的舞台鏡框高度與寬度比,也等於黃金切割比。甚至汽車造型設計中前後輪輪距和車長的比例,也多設在近於黃金切割比。

最後黃金切割搜尋可以寫成如下的虛擬程式:

algorithm 3 Golden Section;

    begin

      input: ;               :, , .

      ;

     

     

      while  do

        begin

            if  then

            begin

                ;

                ;

                ;

            end;

            else

            begin

                ;

                ;

                ;

            end;

            ;

        end;

      output , ;

    end.

◇作業4

寫一黃金切割法的電腦程式,並用作業3中的四個函數來作測試,演算法中的中止要件。◇

4.         多項式近似法

第二、第三小節中介紹的單變數函數零點或最小值搜尋方法,都屬於零次方法,只需要單變數函數的函數值,而不需要一次微分,且只要求函數值為連續,而不要求一次微分為連續。這一節則要介紹多項式近似法--「牛頓法(Newton's method)」求函數零點或最小值,多項式近似法屬於序列近似法,需要函數的一次或二次微分資料,可歸類為一次方法或者二次方法。

4.1 牛頓法求單變數函數零點

牛頓法解單變數函數零點的原理,和2.2節中討論的割線法十分類似,割線法中是以函數上兩點形成的割線作線性近似,牛頓法則是以函數某一點上的切線作線性近似。給定函數上一點,以及這一點的函數值和一次微分後,在這一點上的線性近似式可以寫成

                                                                                        (5)

(5)X軸交點則為

                                                                                                (6)

(6)即為牛頓法求單變數函數零點演算法中的迭代定義。將割線法稍作修改,便可以得到如下牛頓法的虛擬程式:

algorithm 4 Newton’s Method for Finding a Zero;

    begin

        input:

        ;

        while  do

            begin

               

                ;

            end;

                output , ;

    end.

◇作業5

寫一牛頓法求單變數函數零點的電腦程式,以作業1中的例題作測試。其收斂速率與二分法或割線法比較如何?試由不同的起始點出發,看是否有不收斂的情況。◇

4.2 牛頓法求單變數函數最小值

回顧第一章工程最佳化設計的樣本問題中,目標函數為內隱式函數,為求得其最小值,第一章中先以有限元素分析計算出變數a5個不同值時,的函數值(變數b的值由等式限制條件決定),再以曲線適應(Curve Fitting)求出之二次近似式,,利用此二次近似式,求出在a=0.47時,之最小值為30.60mm(如圖7)。

7. 樣本問題中以二次多項式近似內隱式限制條件

接下來再次修改有限元素分析模型,使a=0.47mb=0.88m,實際進行有限元素分析,所得最大變形量31.03mm,與近似方程式預測值相差1.41%,這個誤差尚可接受,便不再繼續做迭代。

以牛頓法求單變數函數最小值,原理上是以一元二次多項式去近似此單變數函數,和前面敘述的程序非常相似,但是樣本問題的處理中是以變數a散佈在在整個定義域中數個資料點,嘗試做「全域近似(global approximation)」,牛頓法則是以某一個設計點上,函數的二次微分資料做「區域近似(local approximation)」。

牛頓法以初始設計點上函數的二次微分資料得到原單變數函數的一元二次多項式近似式後,該一元二次多項式的最小點很容易求出,而每次迭代中便以此近似最小點為下次迭代點,逐步逼近原函數的最小點。如此以一系列簡單、已知的問題,去近似原本複雜、未知的問題,這也正是後續將正式介紹序列近似法的基本概念。

接續前一節牛頓法解單變數函數零點的演算法,牛頓法求單變數函數最小值演算法的推導其實十分容易。在式(5)(6)中,假設另一函數F,且,那麼求函數f的零點便相等於求函數F的最小點。此時式(6)的迭代定義成為

                                                                                               (7)

這個演算法要特別注意的是,如果在式(7)為負值,也就是說在設計點函數的一元二次多項式近似式為一上凸函數而非下凹函數,這時候點為一元二次多項式近似式的最大值而非最小值,則牛頓法也會發散。

◇作業6

修改作業5中的程式,使成為一個牛頓法求單變數函數最小點的電腦程式,以作業3中的例題作測試。其收斂速率與黃金切割法比較如何?試由不同的起始點出發,看是否有不收斂的情況。◇

5.         演算法體系整理

最後圖8將本章介紹的演算法做一個體系上的整理。本章介紹了6個求單變數函數零點和單變數函數最小值的演算法,就演算法的結構來說,可以分為“縮小不確定區間”(包括二分法、費邦那西搜尋、黃金切割搜尋)和“序列近似”(包括割線法和兩種牛頓法)兩個體系。縮小不確定區間的方法,在基於初始不確定區間包含函數零點或最小點的前提下,每次迭代縮小此不確定區間,演算法設計的重點在於如何最有效率、最快速地縮小不確定區間。序列近似法的精神則在“以一系列簡單、已知的問題,去近似原本複雜、未知的問題”,期待簡單、已知問題的解會逐步近似原本複雜、未知問題的解。此外圖8中每個演算法也標明是屬於零次方法、一次方法、或二次方法。

8. 演算法體系整理