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

多變數函數邊界最小值搜尋

這一章中主要在討論解有限制條件最佳化模型的數值演算法,也就是一般所謂「非線性規劃(non-linear programming)」的範疇。先前提到,求多變數函數內部最小值大部分演算法都是屬於所謂直線搜尋法,不同演算法的迭代定義,只是在定義不同的搜尋方向。然而解有限制條件最佳化模型的數值演算法的發展則有相當大的分歧,主要可以分成兩類,一是序列近似法,二是直接搜尋法。

本章第一節起介紹的序列近似法的基本精神是,“用一系列簡單易解的子問題(sub-problem),去近似原先複雜難解的問題,而希望這一系列子問題的近似解,能夠逐漸逼近原複雜問題的真解”。什麼是“簡單易解的子問題”呢?先前介紹過解無限制條件的最佳化問題的演算法,了解如何求解這類問題後,相對於有限制條件的非線性最佳化問題,無限制條件的最佳化問題自然是較為“簡單易解”的問題。因此本章第一節中所介紹的第一類序列近似法,就是把原先有限制條件的非線性最佳化問題,轉化成一系列的無限制條件求最小值的問題。

在介紹序列近似法的第二大類之前,本章第二節先介紹一類在最佳化領域有非常重要地位的問題-「線性規劃(linear programming)」問題。線性規劃問題顧名思義,其目標函數及限制條件均為線性的函數,這一類問題有非常強健而有效率的演算法,可以迅速解出其最佳解,相對於有限制條件的非線性最佳化問題來說,線性規劃問題又成了“簡單易解”的問題。因此本章第四節中所介紹的第二類序列近似法,便是把原先有限制條件的非線性最佳化問題,轉化成一系列的線性規劃問題。序列近似法的概念,也經常是我們設計演算法求解工程最佳化問題時的基本想法。

本章第三節中所介紹的直接搜尋法,基本精神則與求多變數函數內部最小值的直線搜尋演算法相似,也是定義一個可行且下降的搜尋方向,並在所定義的搜尋方向求最小值,作為下一個迭代點。本章中將介紹可行方向法與梯度投影法兩種不同的演算法。

非線性規劃事實上是一個相當廣泛,且非常理論性的領域,各種非線性規劃的演算法也都已經有商業軟體可供設計者使用。本章對各個演算法的介紹,目的不在對各種演算法背後艱澀的數學理論作深入探討,主要重點是在對各數值演算法發展的系統脈絡、基本精神作一整理,使讀者在使用商業軟體時能對其演算法的特性有概括了解。

事實上絕大部分數值演算法假設所有函數是連續的、精確的、有解析形式而易於計算的,在處理問題的精神上,是把最佳化問題視為純粹的“數學問題”,然而這與實際的工程最佳化問題有相當差距,工程最佳化問題中的函數通常是不連續、不精確、沒有解析形式的。因此設計者在處理一般工程最佳化問題時,往往很難找到現成的演算法或商用電腦程式,可以立即直接套用,通常還是要根據問題的特性修改演算法,甚至重新設計演算法。如果能對各個數值演算法發展的系統脈絡、基本精神有充分了解,對設計者日後自行修改、設計特殊演算法,以處理其特定工程最佳化問題極有幫助。

1.     序列近似法(I)-懲罰函數法

這一節裡將要討論「懲罰函數法(penalty function method)」,也是前面提到第一類序列近似法,把原先有限制條件的非線性最佳化問題,轉化成一系列的無限制條件求最小值的問題。

首先考慮一個簡單的二變數求最小值問題,作為本章討論各演算法的“原型問題(prototype problem)”:

        min. 

        s.t.   

                                                                                        (1)

這個簡單的二變數求最小值問題如圖1所示,斜線部分為線性目標函數之等高線,左下角為兩個限制條件所形成的可行區域,最小點為兩個限制條件的交點(4, 3)

1. 二變數求最小值的原型問題

懲罰函數法的基本概念,就是組合目標函數與限制條件成為一個單一的「虛擬目標函數(pseudo-objective function),如式(2)所示。虛擬目標函數的特性,就是在設計點違反限制條件時給予“懲罰(penalty)”,增大虛擬目標函數的值,因此在求此虛擬目標函數的最小值的過程中,除了讓目標函數值變小之外,懲罰函數的存在使所得設計點不至違反限制條件。如式(2)中的函數P(x)便稱作「懲罰函數(penalty function)」。

                                                                      (2)

一個典型的懲罰函數定義如下:

                                                          (3)

(3)懲罰函數第一項中為不等式限制條件,其大於零的部分要給予“懲罰”,第二項中則為等式限制條件,若其不等於零便會得到“懲罰”。而(3)對限制條件違反越大,懲罰函數的數值也越大,與違反限制條件值的平方成正比。

因此懲罰函數法是以一系列求虛擬目標函數最小值的問題,逐步逼近最佳點。式(2)中係數的設定十分重要,如果設得太小,懲罰函數對不可行解“懲罰”的效果不夠,仍然可能產生不可行解。但如果設得太大,在虛擬目標函數中f(x)的值相對而言影響力太小,求的最小值時僅在求不違反限制條件,反而沒有在降低f(x)的值。

在迭代過程中,每次迭代時的設定也有所不同,而設定的策略則和懲罰函數的形式是所謂「外懲罰函數法(exterior penalty function method)」或者是「內懲罰函數法(interior penalty function method)」有關。外懲罰函數法每次迭代時的設定由小而大,也就是先求降低目標函數值再使其解不違反限制條件,從可行區間外逐步向可行區間逼近內懲罰函數法每次迭代時的設定則由大而小,也就是先求出可行解再漸漸降低目標函數值,從可行區間內逐步向邊界逼近。每次迭代得到的設計點則用作下一次迭代時的起始點。

2是式(1)中的原型問題取時的虛擬目標函數等高線圖。從圖2中可以看出虛擬目標函數的最小點與原來原型問題圖1中的最小點中十分相近。然而虛擬目標函數在原模型中可行區域邊界的附近不平滑,例如圖2中,在原問題的可行區域內,懲罰函數沒有作用,虛擬目標函數就是原先的線性函數,而越過可行區域邊界後虛擬目標函數又因為加了懲罰函數而突然變成非線性函數。事實上可以證明出在這個可行區域的邊界,虛擬目標函數的二次微分不連續,求其最小值時容易發生數值問題。此外太大時虛擬目標函數在原先不可行區域函數值變化的梯度極大,求其最小值時也可能會發生數值問題,因此運算上不夠強健是懲罰函數法主要的缺點。

◇作業1

Matlab將此函數繪出一系列類似圖2之等高線圖,將分別設定為0.001, 0.01, 0.1, 1, 10, 100, 1000,並標注每張圖中最小值之大略位置。以此討論內懲罰法和外懲罰法之精神,及設定之方式,及懲罰函數法的缺點。◇

2. 懲罰函數法之虛擬函數等高線圖

另外兩種常見的懲罰函數定義如下:

                                                                                             (4)

                                                                                  (5)

這兩種懲罰函數法都屬於內懲罰函數法,每次迭代時的設定由大而小,接近限制條件邊界時,虛擬目標,兩個懲罰函數值皆趨近無限大,像是在可行區域邊界築了一道高牆障蔽,使迭代過程中設計點不會進入不可行區域,因此這樣的懲罰函數也叫作「障蔽函數(barrier function)」。式(4)(5)兩個懲罰函數的不同,則在於式(5)的函數在可行區域邊界函數值變化的梯度較小,較不容易發生數值問題。

最後懲罰函數法程序可以用以下虛擬程式表示,注意在這個演算法中,如果是內懲罰法由大到小,因此;如果是外懲罰法由小到大,因此

  algorithm 1  Penalty Function Method;

      begin

           input: , , g;

           ;

           while TERMINATION = false do

                begin

                        Form pseudo-objective function ;

                        Solve min.  for ;

                        ;

                        ;

                end;

           output , ;

      end.

2.     線性規劃

許多學者把「線性規劃(Linear Programming, LP)」的發展列為二十世紀中期最重要的科學進步之一。線性規劃應用的領域非常廣泛,從經濟、財務金融規劃、到工廠排程、作業研究、以及其他工程或科學應用問題上,至今已經成為一個標準工具。

解線性規劃的演算法-由美國史丹福大學教授Dantzig1947年發明的「簡算法(simplex method)」,是一個非常成熟而強健的演算法,可以有效率地解出龐大的線性規劃問題,簡算法的發明也是線性規劃得到廣泛重視與應用的重大原因。這個演算法的概念,與本章前面討論解有限制條件最佳化模型演算法的兩個基本原則都不相同,簡算法既不是直接搜尋法也不是序列近似法,其搜尋方式是沿著可行區域邊界搜尋。本節的目的在概略介紹線性規劃與簡算法,而線性規劃問題將成為下面幾節中介紹的許多非線性規劃法的根本。

2.1 線性規劃模型

線性規劃在工程上的應用,最常見的應該是作業研究、工廠排程等工業工程相關的問題。以下便舉一個一般線性規劃教科書常見,簡單的工廠排程的例子[Hillier and Lieberman, 1990]

◇例題1

某公司生產玻璃窗門,該公司有三個工廠,工廠一生產鋁門窗架,工廠二生產木門窗架,工廠三生產玻璃並組合門窗。今三個工廠均有剩餘產能,而該公司計劃利用此剩餘產能生產兩個新產品,產品一為新型鋁門窗,產品二為新型木門窗,三個工廠可利用產能及兩個產品所需各工廠產能以及兩個產品的單位利潤如下表。

 

產品一

產品二

可利用產能

工廠一

1

0

4

工廠二

0

2

12

工廠三

3

2

18

單位利潤

3

5

 

產品一、產品二應採用何種生產組合,公司才能得到最大利潤?

首先令產品一、產品二的產量分別為,這個生產規劃問題可以很輕易寫成如下的最佳化模型:

        max.        

        s.t.           

                       

                       

                                                                                                        (6)

線性規劃顧名思義,其目標函數與限制條件均為線性函數,式(6)便是一個線性規劃模型。包含有n個設計變數及m個限制條件的線性規劃問題,其數學標準型式可以表示如下:

        max.        

        s.t.           

                                                                                                     (7)

注意此處標準形式中為“最大化”而非“最小化”,這主要是為了配合下一小節將要介紹的簡算法的形式,實質意義上並無不同。許多書中也會把式(7)寫成向量、矩陣相乘的形式,形式上更為簡潔。

2.2 簡算法

(6)的二變數線性規劃模型,可以繪成如圖3。從這個圖中你也許會回憶到,在中學時代就學過類似的問題,將代表目標函數的等高線直線描在圖上,可以從圖上很輕易地看出,這個問題的最佳解發生在可行區域的邊界的角點=(2, 6)

我們為什麼會有這樣的直覺呢?線性規劃問題中所有變數在目標函數、限制條件中都具有單調性,對於每一個變數一定有一個相對應的有效限制條件。因此線性規劃問題一定是「以限制條件界定」的問題,也就是說在最佳解時有效限制條件的個數等於設計變數的個數,對於一個n變數的線性規劃問題來說,最佳解一定發生在n個有效限制條件的交點上。簡算法正是根據線性規劃問題的這個特性而設計的,有系統地在可行區間邊界,逐個搜尋限制條件的交點,判斷哪些是有效的限制條件。

3. 二變數的線性規劃模型

簡算法的幾個重要步驟可以歸納如下:

(1)       起始步驟:自一可行區域的角點上開始搜尋。

(2)       迭代步驟:移動至相鄰更佳的角點上。

(3)       最佳測試:相鄰角點已經沒有比現在更佳時,即中止搜尋。

簡算法的這三個步驟觀念上十分簡單,以下便以式(6)的線性規劃模型為例,說明其數學演算過程。

(1)   起始步驟

這裡首先要介紹線性規劃中三個專有名詞:「基本變數(basic variable)」、「非基本變數(non-basic variable)」、及「差額變數(slack variable)。首先在式(6)中,我們可以把所有不等式限制條件左側與右側之間的“差額”補滿,改寫成如下等式限制條件的形式

        max.        

        s.t.           

                       

                       

                                                                                         (8)

其中新增加的變數便稱作「差額變數」,分別代表不等式限制條件的不等號左側與右側常數項41218的差額。差額變數為零時,表示其所對應之限制條件為有效限制條件,差額變數為正值時,表示其所對應之限制條件為無效限制條件。

簡算法中並沒有把差額變數和一般變數作特別區分,觀察式(8)與圖3,事實上每一個變數等於零時都代表了可行區域的一個邊界,差額變數為零時,分別代表三個限制條件的邊界,而一般變數為零時,則分別代表橫軸及縱軸的邊界。

在簡算法起始步驟中要選擇一個可行區域的角點作為搜尋出發點,最自然的選擇便是座標原點(0, 0),此時變數為零,而變數為正值。這裡設為零的變數便稱作「非基本變數」,設為正值的變數則為「基本變數」

在進行簡算法時,首先將式(8)改寫為如下的「表列形式(tabular form)

基本變數

f

 

 

1

-3

-5

0

0

0

0

0

1

0

1

0

0

4

0

0

2

0

1

0

12

0

3

2

0

0

1

18

注意基本變數在目標函數與三個限制條件中各只出現一次,且其對應的係數為1,同一式中僅有一個基本變數,而其他非基本變數又均為零。這樣形式的安排,使得目標函數值與基本變數值可直接在此表右端讀出,即

(2a) 迭代步驟1

迭代步驟中要由現在的角點移動至相鄰更佳的角點上,這其中便牽涉到兩個問題,一是在簡算法程序中如何判斷兩個角點為相鄰?其次是如何知道此相鄰角點比現在的角點更佳

移動至相鄰的角點,從幾何圖形上來說便是離開了某一條可行區域邊界直線(相對應原先為零的非基本變數成為基本變數),而接觸到另一條可行區域邊界直線(相對應的基本變數改設為零成為非基本變數)。如圖3中,起始點(0, 0)有兩個相鄰的角點,分別是(4, 0)(0, 6),由(0, 0)移動至(4, 0),碰到了限制條件,代表其對應的基本變數變成零而成為非基本變數,而離開了的邊界,代表原先的非基本變數此時大於零,成為基本變數。

至於起始點(0, 0)這兩個相鄰角點(4, 0)(0, 6)何者較佳呢?我們要檢查目標函數中變數的係數,這裡目標函數中變數的係數是+3,而變數的係數是+5,因此將增加顯然能更快速的增加目標函數值,因此我們的移動方向便選擇增加變數的方向。

決定往增加變數的方向移動後,下一個問題是,變數出現在兩個限制條件中,變數增加,應該停在哪一個角點?變數持續增加,最後將停在首先碰到的限制條件上,從表列形式中可輕易看出,在限制條件,而在限制條件,因此在這個迭代步驟中,變數增加將先碰到限制條件,設計點由(0, 0)移動至(0, 6),此時相對應的基本變數變成零而成為非基本變數,而原先的非基本變數成為基本變數。

(2-b) 迭代步驟2

移動至點(0, 6)後,(限制條件),代入其他各式消去變數後,這個線性規劃問題的表列形式可以更新如下:

基本變數

f

 

 

1

-3

0

0

2.5

0

30

0

1

0

1

0

0

4

0

0

1

0

0.5

0

6

0

3

0

0

-1

1

6

注意在這個消去過程中要達成的形式,是使基本變數僅在其相對應的限制條件中出現,且係數為1,如此才可直接從表右端目標函數值與基本變數值,在這裡

這個解是不是最佳解呢?此時的目標函數仍然有非基本變數的係數為正值,如果把這個非基本變數轉換成基本變數(也就是從零變成正值),目標函數有繼續增加的可能。接下來的處理方式和前面相同,變數出現在兩個限制條件中,變數持續增加,在碰到限制條件,而在碰到限制條件,因此在這個迭代步驟中,變數將停在限制條件上,設計點由(0, 6)移動至(2, 6),此時限制條件相對應的基本變數變成零而成為非基本變數,而離開的邊界後,原先的非基本變數便成為基本變數。

(3)   最佳測試

移動至點(2, 6)後,(由限制條件),代入其他各式消去基本變數,使基本變數僅在其相對應的限制條件中出現,且係數化簡為1之後,這個線性規劃問題的表列形式又可以更新如下:

基本變數

f

 

 

1

0

0

0

3/2

1

36

0

0

0

1

1/3

-1/3

2

0

0

1

0

1/2

0

6

0

1

0

0

-1/3

1/3

2

表右端可以讀出目標函數值與基本變數值分別為。檢查此時目標函數為目標函數中所有非基本變數係數均為負,可知由此點繼續移動無法進一步改進目標函數值,此點必為最佳點。此時差額變數為零,表示其所對應限制條件為有效限制條件,而差額變數大於零,表示其所對應限制條件不是有效限制條件。

◇作業2

某工廠有車床、銑床、磨床三種機具,每種機具的可使用時間車床為每週350小時,銑床每週500小時,磨床每週150小時。該工廠所生產的三個產品中每單位產量需此三種機具加工時間如下:

 

產品一

產品二

產品三

銑床

9

3

5

車床

5

4

0

磨床

3

0

2

銷售部門表示,產品三每週最大銷售量是20單位,產品一、產品二則無限制。產品一、產品二、產品三的單位利潤分別為50元、20元、25元。此工廠要決定如各產品各需生產多少單位,才能得到最大利潤。

將此問題建立一個線性規劃模型,並仿照前例,用簡算法解出,請清楚標明每一步驟。發揮耐心,一生之中只有這一次用手算解簡算法。◇

◇作業3

閱讀“最佳化電腦數值分析軟體GAMS使用簡介”,將前述線性規劃問題以GAMSWebOpt、或其他解最佳化問題之數值軟體解出。試解讀GAMS之輸出檔,並將結果與手算結果比較。◇

3.     直接搜尋法

除了序列近似法外,解有限制條件最佳化模型的數值演算法的第二類方法,是延續求多變數函數內部最小值的直線搜尋演算法精神,定義一個搜尋方向,並在所定義的搜尋方向求最小值,作為下一個迭代點。所不同的,是在無限制條件時,這個搜尋方向僅必須是下降方向,而在有限制條件情況的最小值搜尋過程中如果碰到可行區域邊界,這個搜尋方向必須同時是滿足限制條件的可行方向

3.1 可行方向法

先前曾經討論過下降方向與可行方向的定義,滿足式(9)的搜尋方向即為下降方向,而設計點落在限制條件的邊界上時,滿足式(10)的搜尋方向為可行方向。

                                                                                                    (9)

                                                                                                  (10)

但如圖4所示,在限制條件為非線性時,滿足的搜尋方向s很容易違反限制條件走入不可行區域,因此「可行方向法(feasible direction method)中可行方向的定義式(10)中,通常會加入一個與相關的正值係數,使搜尋方向稍微遠離可行邊界,如式(11)所示。

4. 加入係數q使搜尋方向稍微遠離可行邊界

                                                                              (11)

而我們要選擇的搜尋方向s,必須要使目標函數下降越多越好,即要使的值為最小,換一個表達方式,也就是使下式中的b值最大:

                                                                                              (12)

整理式(9)至式(12),求這個可行方向s可以轉化成以下線性規劃問題:

        max.         b

        s.t.           

                        ,

                        , i = 1, ..., n                                                            (13)

其中J為在設計點上有效限制條件的集合,對非線性限制條件建議值為1,對線性限制條件則為0,即搜尋方向保持在線性限制條件邊界上。在這個線性規劃問題中的變數數值並沒有大小限制,式(13)中最後的範圍限制主要在使此限制規劃問題不至於限制條件不足。

可行方向法在求出這個搜尋方向s之後,便如同搜尋內部最小值之直線搜尋演算法作直線搜尋,解出下一次迭代的設計點,但此時步長的決定除了要使目標函數在搜尋方向s上為最小之外,尚須注意不能使設計點違反限制條件。

◇作業4

考慮圖5中這個簡單的靜不定三桿件衍架結構最佳化問題,桿件123的截面積分別為,在三桿件的接點處受一大小為P的橫向作用力,及受一大小為8P方向向下的縱向作用力;Er分別為此靜不定衍架材料的彈性模數及比重。這是一個對稱結構,橫向作用力P的方向可以向左或向右,因此=。這個問題的目標是求取最佳的三桿件截面積值,使得此結構的重量為最輕,限制條件則為三桿件的最大應力值不可超過容許應力

5. 三桿件衍架結構

由基本的結構力學,可以推導出這個最佳化問題的數學模型:

        min. 

        s.t.   

               

               

其中限制條件為桿件13的最大應力值之限制條件,為桿件2的最大應力值之限制條件。為了討論上的簡便、明確,這裡先假設下列各個參數的數值:,同時,代入前式化簡,可以得到下式:

        min. 

        s.t.   

               

Matlab畫出這個最佳化設計模型的目標函數和限制條件。由點(100, 175)出發,求出在此點上可行方向法的搜尋方向。注意你可能需要如GAMS之類的軟體,來幫助你解式(13)的線性規劃問題。◇

3.2 梯度投影法

「梯度投影法(gradient projection method)也是一種直接搜尋法,其定義搜尋方向的基本精神與線性規劃法類似,也是希望沿者可行區域的邊界作搜尋。回顧搜尋內部最小值的演算法,在目前設計點上沒有有效限制條件時,梯度法是以為搜尋方向,而在目前設計點上有有效的限制條件時,梯度投影法的搜尋方向則定義為目標函數的梯度在該有效限制條件上的投影,也就是說希望定義一個沿著有效限制條件邊界移動的下降搜尋方向

假設在某迭代點有k個有效限制條件,則在此點上梯度投影法中的搜尋方向定義如下:

                                                                                        (14)

上式中第二項在修正目標函數的梯度,使得此搜尋方向能夠滿足

                                                                                                      (15)

將式(14)代入式(15)k個拉氏乘數可由下式中求得:

                                                      (16)

與前一節的可行方向法相同,對於非線性的限制條件,在有效限制條件上的直接投影,也就是滿足的搜尋方向很容易會違反限制條件進入不可行區域,因此對於每一個有效限制條件,這裡也要引入一個相對應的正值參數,使得搜尋方向s稍微遠離不可行區域邊界:

                                                                                                (17)

因此式(16)便成為

                                           (18)

其中參數的選定則和現在設計點對該有效限制條件是否為可行解有關,基本上如果越不可行,應該取得越大。

梯度投影法中從某一設計點之搜尋方向的定義,便是計算該點上的,設定後先解出滿足線性聯立方程式(18)中的,再代入式(14)求得搜尋方向s

◇作業5

在三桿件問題中,由點(100, 175)出發,求出在此點上梯度投影法的搜尋方向,並用Matlab畫出來,與可行方向法之搜尋方向比較之。

4.     序列近似法(II)

前一節中討論的線性規劃法,有非常強健而有效率的演算法及商業電腦程式,可以輕易地解出大規模的線性規劃問題,然而大多數科學及工程的最佳化問題中,不論是限制條件或是目標函數對於設計變數通常都是非線性的。從序列近似法的思考模式,解非線性問題演算法設計上很自然的發展方向,便是以一系列簡單易解的線性規劃子問題,來近似原來複雜難解的非線性最佳化問題,這種演算法稱作「序列線性規劃法(sequential linear programming, SLP)。本節中將繼續介紹的其他兩種演算法:「序列非線性規劃法(sequential non-linear programming)」、以及「中心法(method of centers)」,在發展的概念上都是屬於相同系統的序列近似法。

4.1 序列線性規劃法

考慮一非線性最佳化問題:

        min. 

        s.t.    , j = 1, ..., m                                                                     (19)

線性序列規劃法首先將目標函數及限制條件在起始設計點以一階泰勒級數線性化,展開成如下形式:

        min. 

        s.t.    , j=1, ...,m

                                                                                           (20)

(20)是一個線性規劃子問題,可以用線性規劃中的簡算法求解,得到下一次迭代的設計點。一般來說這個新設計點應比起始設計點更為接近原始非線性規劃問題的最小點,而在下一次迭代中,便將式(20)在設計點以一階泰勒級數線性化展開。如此反覆以線性規劃子問題去近似原先的非線性規劃問題,希望每一次迭代得到的新設計點,都比前一個設計點更接近真正的最小點,而在新設計點上的線性近似子問題,也越來越近似原非線性問題最小點的附近區域。最後線性規劃子問題的解,終於可以趨近原先非線性問題的最小點。

(20)中的最後一項限制條件稱為「移動限制(move limit),移動限制的給定,主要是因為利用線性序列規劃法解非線性最佳化問題時,函數是以一階泰勒級數方式近似展開,只有在設計點的附近才會較近似於原來的非線性曲線,超過一定範圍線性近似便相當不準確。因此在每次迭代過程中設計變數的改變量必須以移動限制加以侷限,否則演算過程中很可能會發生振盪及發散的情形。

此外在線性化之後的子問題中,所有變數都是單調性變數,目標函數中很可能有單調變數在線性化之後的子問題中沒有適當的限制條件可以界限住,造成子問題限制條件不足,因此對設計變數的移動限制,也有提供單調變數上下界的功能。

線性序列規劃法的進行流程,可以表示成如下虛擬程式,注意在迭代過程中每次迭代的移動限制大小,都以一個縮減比例係數b來修正。

  algorithm 2 Sequential Linear Programming;

      begin

           input: NLP problem Eq.(19), feasible initial design point ;

                        move limit , reduction factor b;

           ;

           while TERMINATION=false do

                begin

                        Form LP subproblem Eq. (20);

                        Solve LP subproblem for ;

                        while HIGHLY_INFEASIBLE()=true do

                        begin

                                ;

                                Solve LP subproblem for ;

                        end;

                        ;

                end;

           output , ;

      end.

序列線性規劃法每次迭代過程中,設計變數移動限制的定義是非常重要的。序列線性規劃法有一個惱人的特性,是其迭代過程中產生的設計點如果位於不可行區間時,則並不保證此次迭代線性化後的線性規劃模型有一可行解,即使有可行解,也極有可能會造成發散的情況出現。因此在序列線性規劃法中某次迭代出現過度不可行的解時,通常便放棄此解而將變數的移動限制縮小(如在虛擬程式中,移動限制乘上一個縮減比率),重新計算一次。然而序列線性規劃法的中止要件,通常是移動限制小於某一數值,導致每次迭代設計變數改變極小而中止,移動限制值縮小太快,也會使得收斂速率減慢亦或是停留在一非最佳點上。

因此移動限制值大小的改變,實在是序列線性規劃法中關係著控制收斂速率最重要的因素,然而如何適當地決定移動限制值大小而達成較佳的收斂效果,在實際問題的應用上卻是極為困難的一件事。相關文獻中有提到各種決定移動限制值的方法[Arora, 1989; Haftka and Gurdal, 1990; Chen, 1993],但並未有一確定規則來決定移動限制。移動限制決定的困難和其所造成的不穩定性,也是研究最佳化理論的學者向來不推薦使用序列線性規劃法的原因。

以式(1)中的原型問題為例,線性化後的線性規劃子問題如下式:

        min. 

        s.t.  

                                     (21)

設定各變數起始移動限制均為1,縮減比率0.9,從三個不同的起始點(1, 1)(5, 5)(10, 10)出發,可以得到如圖6的迭代過程圖。從圖中可以看出,序列線性規劃法從可行點(1, 1)或接近可行邊界的不可行點(5, 5)出發,都可以成功地收斂至最小點(4, 3),而從非常不可行點(10, 10)出發,序列線性規劃法則很快就發散了。

6. 序列線性規劃法迭代過程圖

◇作業6

嘗試利用序列線性規劃法解三桿件問題,從幾個不同的起始點出發做35次迭代,並繪出類似圖6之迭代歷程。討論你的觀察。注意你可能需要如GAMS之類的軟體,來幫助你解迭代過程中的線性規劃子問題。◇

4.2 序列非線性近似法

序列線性規劃法的原理簡單,在工程最佳化問題實際應用上也相當容易,因此一直是一個很受歡迎的演算法,然而卻有移動限制決定困難,迭代過程收斂性不穩定等缺點。究其失敗的根本原因,與內部最小值搜尋中的梯度法十分類似,是一次線性近似無法充分模擬、近似目標函數與限制條件中的非線性特徵

許多改進的演算法也發展出來試圖對函數非線性部份做比較好的近似,「倒數近似法(reciprocal approximation)便是其中之一。倒數近似法是對原先線性近似中的變數作一變數變換,令,式(20)中的線性化限制條件以一次泰勒式展開式對變數展開,便成為

                                                                (22)

代入上式可得

                                                           (23)

另外一種非線性近似法則是混合了線性近似與倒數近似,稱作「保守近似法(conservative approximation)。保守近似法中所謂“保守”,是在同一迭代點對限制條件的線性近似與倒數近似中,取其數值較大者,因為限制條件的形式是,所以用近似之限制條件數值較小者評估原始限制條件是否不可行,是較為“樂觀”的評估,而用近似之限制條件數值較大者評估原始限制條件是否不可行,是較為“保守”的評估。採用較保守的近似方式,可以減少近似後子問題之最佳解對原始限制條件為不可行的情況。

判斷線性近似與倒數近似何者較“保守”,可將兩者相減判斷何者較大,如-<0則採用,反之則採用

                                                             (24)

因此由的符號可以判知何者較大,我們可以得到保守近似法的近似通式如下:

                                                             (25)

其中如果<0表示倒數近似數值較大(較“保守”),故(即採用保守近似),否則採用線性近似,Gi = 1。如果將目標函數也以保守近似法作近似,近似子問題的可行區間將是一個「凸區間(convex domain)」,因此這種方法被稱作「凸線性化(convex linearization)」,是結構最佳化設計問題中常用的演算法[Braibant and Fleury, 1984]

倒數近似法和保守近似法有倒數項可以較成功地模擬原函數中的非線性,因此其收斂特性也較序列線性規劃法為佳。不過式(23)(24)均非線性,所以倒數近似法和保守近似法的近似子問題都已經不是可用簡算法輕易解出的線性規劃問題了。

◇作業7

選擇一個起始點,利用倒數近似法及保守近似法對三桿件問題做35次迭代,試與作業6的結果作一比較。注意你可能需要如GAMS之類的軟體,來幫助你解迭代過程中的倒數近似或保守近似子問題。◇

這些近似方法在工程最佳化問題上,目標函數或限制條件為內隱式函數時非常有用。從工程最佳化問題的觀點來看,有內隱式限制條件的工程最佳化問題是一個複雜、難解的問題,而有外顯形式的近似問題則有許多商業最佳化電腦軟體可以輕易解出。從這個觀點,在工程最佳化問題上一個很普遍的做法,是僅將問題中的內隱式限制條件作線性、倒數、保守近似,使成為一可用商業最佳化電腦軟體輕易解出的外顯形式問題。如此用一系列“簡單、易解”的外顯形式(但不一定是線性)問題,去模擬原先複雜、難解的內隱形式問題,便稱作「序列非線性規劃法」

4.3 序列二次規劃法

序列線性規劃法中移動限制的形式,是給每一次迭代的設計點一個“方形”的移動區域,而另一種可能的移動限制形式則是給每次迭代的設計點一個“圓形”的移動區域如下式

              (26)

其中附加係數0.5是便於後來微分計算中消去常數之用。

改用圓形的移動區域後,序列線性規劃法中的線性規劃子問題式(20)成為

        min. 

        s.t.    , j=1, ..., m

                                                                                  (27)

(27)也已經不再是線性規劃問題了。如何定義每次迭代時移動限制的大小,是序列線性規劃法中的大問題,改成圓形區域定義移動限制,還是並沒有解決序列線性規劃法的困難,移動步長x訂得太小的話,式(27)仍然可能沒有可行解。然而這個數學模型還可以進一步變形,把移動步長的限制條件併入目標函數中如式(28)所示,便無需實際去訂定移動步長x了:

        min. 

        s.t.            , j = 1, ..., m                     (28)

(28)是一個「二次規劃(quadratic programming)」問題,也就是其目標函數為二次函數,而所有限制條件均為線性。列出式(27)KKT最佳化條件且設其移動步長限制條件的拉氏乘數為1,你會發現與式(28)KKT最佳化條件完全相同(而且你會了解為何要乘上係數0.5來消去微分產生的常數),也就是說在此狀況下兩者的解完全相同,式(27)和式(28)是相等的最佳化模型,但在式(28)的形式中等於自動選擇了移動步長x的大小。

以式(28)的二次規劃問題為每次迭代的子問題之序列近似法,便稱作「序列二次規劃法(sequential quadratic programming, SQP)。這個演算法不用自行定義移動限制,而每次迭代的子問題式(28)也可用所謂“修正簡算法(modified simplex method)”解出,如果其線性化後的限制條件有可行區域存在,便不用擔心因移動限制訂得太小而造成迭代過程的子問題不可行。序列二次規劃法是一個頗為強健的演算法,也是目前在解大型非線性最佳化問題相當流行且重要的演算法。

◇作業8

選擇一個起始點,利用序列二次規劃法對三桿件問題做35次迭代,試與先前作業的結果作一比較。注意你可能需要如GAMS之類的軟體,來幫助你解迭代過程中的二次規劃子問題。◇

4.4 中心法

序列線性規劃法中每次迭代的線性規劃子問題之最佳解,對原非線性的限制條件如果為不可行時,常易造成整個迭代程序震盪、發散,而這個不穩定性也是序列線性規劃法嚴重的缺點。本節所要介紹的「中心法」,便是針對改進線性序列規劃法此項缺點所設計的演算法。

線性規劃問題的最佳解一定發生在限制條件的邊界的交點上,而在線性序列規劃法中,原限制條件有非線性的情況下,限制條件線性化後限制條件的交點通常會落在原限制條件的不可行區域內。因此中心法的重要特性,便是每次迭代所採用的設計點,不用前一次迭代線性規劃子問題的最佳點,而採用前一次迭代線性規劃子問題中,線性化後的目標函數與限制條件所建構區域中,所能包含最大“圓”的圓心,如此,希望每一次迭代所得到的設計點,仍然落在限制條件定義的可行區域內。

這個所謂的“圓”,在二變數平面問題中是一個圓,在三變數問題中便是一個“球”,而在變數數目超過三個時,變成了所謂“超球(hyper sphere)”。而求線性化後的目標函數與限制條件所建構區域中,所能包含最大“圓”的圓心,基本上也是解一個線性規劃問題。以本章的原型問題為例,以點(2, 2)為起始點,在點(2, 2)的目標函數與限制條件的線性近似如圖7所示。注意在中心法中起始點的選擇,必須是一個可行解。

7. 中心法

圖中線性化後的目標函數與限制條件所建構區域中,所能包含最大圓的圓心為A點,而點(2, 2)A點間的向量s即為本次迭代設計點移動的向量,求出s後即可求出下次的迭代點A

首先我們要求點A到直線f的距離。由圖6中的幾何關係可知,中心點A至目標函數的線性近似直線之距離,為向量s在目標函數在點(2, 2)梯度方向的相反方向之投影,即

                                                                                          (29)

而中心點A至限制條件的線性近似直線之距離則為向量在限制條件在點(2, 2)梯度方向的相反方向之投影,到任一限制條件的距離的通式可以寫成

                                                                              (30)

但限制條件的近似直線是原限制條件的一次泰勒氏展開式,因此在

                                     (31)

代入式(30)中可以得到

        , j = 1, ..., m                                                  (32)

得到圓心A與目標函數的近似直線之距離,以及圓心A與各限制條件的近似直線之距離後,我們可以推論這些近似直線所包含最大圓的半徑r應該小於,即

                                                                                     (33)

        , j = 1, ..., m                                                        (34)

最後整理式(33)(34),可得求此最大圓圓心的線性規劃子問題如下:

        max.         r

        s.t.   

                , j = 1, ..., m

                , i = 1, ..., n                                                                  (35)

注意式(35)中對於向量s的每一個元素,仍然訂有移動限制,但此處移動限制的主要功能是防止此線性規劃子問題限制條件不足。

中心法每次迭代求出的這個圓心,離線性化後的限制條件有一段距離,因此中心法每次迭代所選取的設計點,都應該是可行的設計點,這也是中心法最吸引人的地方。如此反覆以線性規劃子問題去求此圓心位置,希望每一次迭代得到的新設計點,都比前一個設計點更接近真正的最小點,而在新設計點上的線性化後的目標函數與限制條件所建構區域也越來越小。最後半徑r將趨近於零,而線性規劃子問題的解,終於可以趨近原先非線性問題的最小值。

中心法解決了線性序列規劃法迭代過程中因不可行解產生造成的不穩定性,然而中心法在接近最佳點時收斂的速度非常慢,因為中心法每次迭代不直接走到限制條件邊界,理論上來說要經過無限多次迭代才能收斂到理論最佳點。但是中心法每次迭代產生之設計點均為可行設計點,在幾次迭代中便能不斷改進設計點,求得一個接近理論最佳點的設計,這個特性在工程最佳化問題上,倒是十分實用的。

◇作業9

選擇一個起始點,利用中心法對三桿件問題做35次迭代,試與先前作業的結果作一比較。注意你可能需要如GAMS之類的軟體,來幫助你解如式(35)的線性規劃子問題。◇

參考資料

Arora, J. S. 1989. Introduction to Optimum Design, McGraw-Hill, New York.

Braibant, V. and Fleury, C., 1984. “Shape Optimal Design Using B-Spline.” Computer Methods in Applied Mechanics and Engineering, Vol. 44, pp. 247-267.

Chen, T.-Y., 1993. “Calculation of Move Limits for the Sequential Linear Programming Method.” International Journal for Numerical Methods in Engineering, Vol. 36, pp. 2661-2679.

Haftka, R. T., Gurdal, Z., Kamat, M. P., 1990. Element of Structural Optimization, 2nd Ed., Kluwer Academic Publishers, pp. 223-225, 251-281.

Hillier, F. S., and Lieberman, G. J., 1990. Introduction to Operation Research, McGraw-Hill, New York.