工程最佳化方法中的“最速下降法”和“DFP擬牛頓法”的 C 語言實現

標籤：des   style   blog   http   使用   strong   

這個小程式是研一上學期的“工程最佳化”課程的大作業。其實這題本可以用 MATLAB 實現，但是我為了鍛煉自己薄弱的編碼能力，改為用 C 語言實現。這樣，就得自己實現矩陣的運算（加減乘除、求逆、拷貝）；痛點是求偏導，通過查資料，發現可以通過導數定義，即取極限的方法，來逐步逼近求得梯度；另外，沒法做到輸入任意公式，只能將公式寫入程式碼為函數，而求導函數需要傳入公式，就直接傳入函數指標了。思考、編碼、調試、測試共耗費兩周左右時間，完成於 2013/01/10。雖然為了認真做這個大作業而耽誤了期末考試的複習，但我不後悔做出的選擇，因為我學到了我覺得真正有用的東西。

源碼託管在 Github 上：點此開啟連結

一、題目

用最速下降法和DFP擬牛頓法求解以下函數的最小值點以及最小值：

1.1 ，其中，，，

1.2 ，其中，，，

二、演算法

2.1最速下降法（steepest descent method）

演算法步驟：

（1）取初始點，精度,令；

（2）計算,若,則停，；否則轉（3）；

（3）一維搜尋：，

令,轉（2）。

2.2擬牛頓法（DFP）

演算法步驟：

（1）取初始點，允許誤差；

（2）求，若，令，演算法停止；否則轉（3）；

（3）令；

（4）令；

（5）求：，令；

（6）求，若，令，演算法停止；否則轉（6）；

（7）若，則令，，轉（3）；

否則令，，

計算；

令，轉（4）。

2.3成功—失敗法（用於一維搜尋）

演算法步驟：

（1）取初始點，初始步長和精度,計算；

（2）計算；

（3）若（搜尋成功），令；

若（搜尋失敗），若，令，停止迭代；

否者，令，轉（2）；

三、語言及演算法實現說明

3.1演算法實現語言及平台：

C語言+VC6.0（Debug模式）。

3.2幾個部分的思考：

（1）由於實現即時輸入函數多項式比較困難，本程式將函數多項式寫成模組，存入程式檔案中，由於程式使用函數指標，故可以陸續添加函數多項式而不必修改核心演算法的代碼；

（2）由於函數不同，取值範圍不同，則演算法需要不同的精度和步長，才能求得精確的結果，故本程式提供介面讓使用者指定；

（3）為了實現即時輸入變數維度，本程式使用動態記憶體分配，建立多維陣列，類比矩陣，用於儲存多維變數；

3.3演算法實現的重痛點分析：

（1）偏導數的求解：本程式使用偏導數的定義，即極限方法，求解指定點的函數值；

（2）DFP演算法中的計算：本程式用多維陣列來類比矩陣進行運算。

四、程式中的主要模組說明（完整程式及注釋見附錄）

4.1待求解的兩個函數：

其中vars為多維變數，n代表維度，這兩個模組返回函數在指定點的值。

/* 求函數1在指定點的值 */

double fun1(double **vars, int n)；

/* 求函數2在指定點的值 */

double fun2(double **vars, int n)；

4.2利用偏導的定義求某個點的偏導數：

其中f為指定函數，vars為多維變數，grads為梯度，n為維度，prec為使用者指定的精度；該模組求出函數的偏導存入矩陣grads中。

/* 用極限方法求指定點的偏導/梯度 */

void differ(double (*f)(double **vars, int n), double **vars, double **grads, int n, double prec)；

4.3成功—失敗法，用於一維搜尋：

其中f代表指定函數，vars為多維變數，d為二維搜尋的方向，n為維度，prec為使用者指定的進度，h為使用者指定的步長；

該模組將搜尋到的所對應的多維變數存入矩陣vars。

/* 成功失敗法，用於一維搜尋 */

void suc_fail(double (*f)(double **vars, int n), double **vars, double **d, int n, double prec, double h)；

4.4兩個核心演算法：

其中fun為待解函數的標號，n為維度，prec為使用者指定的精度，h為使用者指定的用於一維搜尋的步長；

這裡這兩個模組求出指定函數的最小值點和最小值並輸出。

/* 最速下降法（Speedest Descent Method）*/

void SD(int fun, int n, double prec, double h)；

/* DFP擬牛頓法 */

void DFP(int fun, int n, double prec, double h)；

五、程式使用說明

本程式將最速下降法和DFP法整合在一起，精度、步長、維度可由使用者指定：

（1）選擇方法（只輸入序號，‘0’退出）；

（2）選擇函數（只輸入序號）；

（3）輸入精度值（）；

（4）輸入一維搜尋的步長；

（5）輸入變數維度；

（6）輸入變數的每個分量；

斷行符號後程式開始使用指定方法對指定函數進行計算，計算過程中輸出迭代次數；

最後輸出結果：最小值點和最小值。

如所示（下一頁）：

六、運行結果及分析

6.1精度選擇：

（1）如下用最速下降法求函數1，精度取，步長取1，初值取（5,5,5），求解時陷入了無限迭代：

……

（2）對於（1）的輸入，僅修改精度為，僅迭代3次就求出了結果，且達到很高的精度，變數的三個分量和最優值都約等於0：

6.1.1小結

當精度值選擇太小，雖然可能得到更精確的結果，但會陷入死迴圈。當精度要求放鬆了一點，反而快速求出了精確結果，可見精度要選著適當，不可太大，也不可太小。以下實驗就選擇為精度值。

6.2一維搜尋的步長選擇：

（1）如下用最速下降法求函數1，精度取，步長取0.1，初值取（3,3,3），迭代3次求出結果，但是誤差很大：

（2）針對（1），僅將步長改為0.5，迭代4次求出結果，精度很高：

（4）如下用最速下降法求函數2，精度取，步長取0.5，

初值取（300,300,300），迭代9次求出結果，但是誤差很大：

（5）針對（4），僅將步長改為30，迭代17次求出結果，雖然結果與理想值0還是有一些誤差，但比（4）的結果精確了很多：

6.2.1小結

一維搜尋的步長也要選擇適當，否者求出的結果誤差很大。從以上對比可以看出，步長的選取要根據自變數的取值進行相應的調整：函數F1的，變數取3，步長h取0.5時誤差較小；函數F2的，變數取300，步長h取30時誤差較小，步長h取值為變數x取值的10%左右時誤差較小。

6.3比較最速下降法和DFP法：

6.3.1求解函數F1：精度取，步長取0.5，變數分別取（-5,-5,-5）、（5,5,5）

（1） 最速下降法

（2） DFP

6.3.2求解函數F2：精度取，步長取50，變數取（500,500,500）

（1）最速下降法

（2）DFP

6.3.3小結

由以上兩組對比可看出：

（1） 對於函數F1和F2，DFP演算法都比最速下降法迭代次數多；

（2） 對於函數F1和F2，DFP演算法都比最速下降法結果精確；