1、大學數(shù)學實驗大學數(shù)學實驗Mathematical Experiments 實驗實驗7 7 無約束優(yōu)化無約束優(yōu)化 最優(yōu)化最優(yōu)化(Optimization)是工程技術、經(jīng)濟管理、是工程技術、經(jīng)濟管理、科學研究、社會生活中經(jīng)常遇到的問題科學研究、社會生活中經(jīng)常遇到的問題, 如如:優(yōu)化模型和算法的重要意義結構設計結構設計資源分配資源分配生產(chǎn)計劃生產(chǎn)計劃運輸方案運輸方案解決優(yōu)化問題的手段解決優(yōu)化問題的手段 經(jīng)驗積累，主觀判斷經(jīng)驗積累，主觀判斷 作試驗，比優(yōu)劣作試驗，比優(yōu)劣 建立數(shù)學模型，求解最優(yōu)策略建立數(shù)學模型，求解最優(yōu)策略最優(yōu)化最優(yōu)化: : 在一定條件下，尋求使目標最大在一定條件下，尋求使目標最大(

2、(小小) )的決策的決策 運籌學運籌學(OR: Operations/Operational Research)管理科學管理科學(MS: Management Science)決策科學決策科學 (DS: Decision Science)(最)優(yōu)化理論是運籌學的重要內(nèi)容無無約約束束優(yōu)優(yōu)化化OR/MS/DS優(yōu)化優(yōu)化(Optimization), 規(guī)劃規(guī)劃(Programming)線線性性規(guī)規(guī)劃劃非非線線性性規(guī)規(guī)劃劃網(wǎng)網(wǎng)絡絡優(yōu)優(yōu)化化組組合合優(yōu)優(yōu)化化整整數(shù)數(shù)規(guī)規(guī)劃劃不不確確定定規(guī)規(guī)劃劃多多目目標標規(guī)規(guī)劃劃目目標標規(guī)規(guī)劃劃動動態(tài)態(tài)規(guī)規(guī)劃劃 優(yōu)化問題的數(shù)學模型可行域目標函數(shù)，決策變量，fxnxRxxfz

3、),(min 可行解可行解（只滿足（只滿足(2)）與）與最優(yōu)解最優(yōu)解（滿足（滿足(1),(2)） 無約束優(yōu)化無約束優(yōu)化（只有（只有(1)）與）與約束優(yōu)化約束優(yōu)化（(1),(2)）)2(, 2 , 1, 0)(. .)1 ()(minmixgtsxfzix 實際問題一般總有約束，何時可用無約束優(yōu)化處理？實際問題一般總有約束，何時可用無約束優(yōu)化處理？4. 優(yōu)化工具箱的使用優(yōu)化工具箱的使用2. 無約束優(yōu)化的基本方法：無約束優(yōu)化的基本方法：梯度法，牛頓法，擬牛頓法梯度法，牛頓法，擬牛頓法1. 優(yōu)化問題的最優(yōu)解條件；算法模式優(yōu)化問題的最優(yōu)解條件；算法模式無約束優(yōu)化的主要內(nèi)容3. 非線性最小二乘法非線性最

4、小二乘法5. 實際問題中的無約束優(yōu)化模型實際問題中的無約束優(yōu)化模型實例實例1 1 產(chǎn)銷量安排產(chǎn)銷量安排假設假設A產(chǎn)銷平衡產(chǎn)銷平衡假設假設Bp隨隨x (兩種牌號兩種牌號)增加而減小，呈線性關系增加而減小，呈線性關系12111211121211111, 0,aaaabxaxabp某廠生產(chǎn)兩個牌號的同一種產(chǎn)品，如何確定產(chǎn)量使利潤最大某廠生產(chǎn)兩個牌號的同一種產(chǎn)品，如何確定產(chǎn)量使利潤最大21222221222212122, 0,aaaabxaxabp22211121,)()(),(max21xqpxqpxxzxx目目 標標利潤最大利潤最大0,0,2222221111112211crcerqcrcerqx

5、xq隨隨x (本牌號本牌號)增加而減小，呈負指數(shù)關系增加而減小，呈負指數(shù)關系假設假設C實例實例1 1 產(chǎn)銷量安排產(chǎn)銷量安排無約束無約束(非線性非線性)規(guī)劃規(guī)劃x1, x2 0 ?0yxVOR2x=629, y=375309.00 (1.30)864.3(2.0)飛機飛機x=?, y=?VOR1x=764, y=1393161.20 (0.80)VOR3x=1571, y=25945.10 (0.60)北DMEx=155, y=987飛機與監(jiān)控臺（圖中坐標和測量距離的單位是“公里”）實例實例2 飛機精確定位問題飛機精確定位問題 飛機精確定位模型飛機精確定位模型42424)()(),(atan2d

6、yyxxyyxxiii不考慮誤差因素不考慮誤差因素超定方程組，超定方程組，非線性最小二乘！非線性最小二乘！）飛機位置坐標（要求計算：，距離誤差為記測量距離為，角度誤差為記測量角度為標分別為已知數(shù)據(jù)：設備位置坐yxdiiyxiiii, . ; 3,.,1, 1,.4;),(44量綱不符！量綱不符！？ 242424312)()(),(atan2 dyyxxyyxxMiniiiix,y飛機精確定位模型飛機精確定位模型44242444)()(),(atan2dyyxxdyyxxiiiiii考慮誤差因素考慮誤差因素Min x; Min y; Max x; Max y. 以距離為約束，優(yōu)化角度誤差之和（或

7、平方和）；以距離為約束，優(yōu)化角度誤差之和（或平方和）；或以角度為約束，優(yōu)化距離誤差或以角度為約束，優(yōu)化距離誤差. 非線性規(guī)劃非線性規(guī)劃誤差非均勻分布！誤差非均勻分布！ ？ ？僅部分考慮誤差僅部分考慮誤差! 角度與距離的角度與距離的“地位地位”不應不不應不同！同！有人也可能會采用其他目標，如：有人也可能會采用其他目標，如：飛機精確定位模型飛機精確定位模型2424244231)()(),(atan2),( yyxxdyyxxyxEMiniiiii誤差一般服從什么分布？誤差一般服從什么分布？正態(tài)分布！正態(tài)分布！不同的量綱如何處理？不同的量綱如何處理？無約束非線性最小二乘模型無約束非線性最小二乘模型歸

8、一化處理！歸一化處理！角度需要進行預處理，角度需要進行預處理，如利用如利用atan2函數(shù)函數(shù), 值域值域(-pi, pi)無約束優(yōu)化:最優(yōu)解的分類和條件)(xfMinx給定一個函數(shù)給定一個函數(shù) f( (x),),尋找尋找 x* 使得使得 f( (x*)最小，即最小，即nTnxxxx),(21其中其中局部最優(yōu)解局部最優(yōu)解全局最優(yōu)解全局最優(yōu)解必要條件必要條件0),()(1*Txxnffxfx*f(x)xlxgo充分條件充分條件0)(, 0)(*2*xfxfHessian陣陣nnjixxff22求解無約束優(yōu)化的基本思路：下山法基本思想基本思想在在n中某一點，確定一個中某一點，確定一個搜索方向搜索方向

9、及沿該方向的及沿該方向的移動步長移動步長，得到使目標函數(shù)下降的新的點，得到使目標函數(shù)下降的新的點迭代步驟迭代步驟Step 1 初始化：初始點初始化：初始點x0，終止準則等，終止準則等Step 2 迭代改進：方向迭代改進：方向d k，步長，步長 k)()(1kkxfxfkkkkdxx1Step 3 終止檢驗：得到近優(yōu)解或終止檢驗：得到近優(yōu)解或k+1k轉轉2選擇選擇d k , k 使使 f 下降更快下降更快 不同算法不同算法搜索方向的選擇搜索方向的選擇1. 1. 最速下降法最速下降法（梯度法）（梯度法） 下下 降降 方方 向向 最速下降方向最速下降方向迭代改進格式迭代改進格式算法特點算法特點初始階

10、段改進較快，最優(yōu)解附近改進較慢初始階段改進較快，最優(yōu)解附近改進較慢kkTkkkkdxfxfdxfxf)()()()(10)(kkTdxf)(kkxfd)(1kkkxfxx將)(1kxf在kx點作泰勒展開，只保留一階項，有暫不考慮搜索步長，暫不考慮搜索步長，可設可設 k=1( (負梯度方向負梯度方向) )2. Newton方法方法特點特點 局部局部2階收斂階收斂;需計算需計算Hessian陣陣,它可能病態(tài)或不正定它可能病態(tài)或不正定將將f(xk+1)在在xk點作泰勒展開至二階項，用點作泰勒展開至二階項，用d替代替代dk 牛牛 頓頓 方方 程程 牛牛 頓頓 方方 向向迭迭 代代 格格 式式)()(2

11、kkxfdxfdxfddxfxfdxfxfkTkTkkk)(21)()()()(21)()(12kkkxfxfd)()(121kkkkxfxfxx)()(11kkkkxFxFxx比較比較的牛頓法解0)(xF)()(xfxF求d使f(xk+1)極小右端對d導數(shù)為00)()(2dxfxfkk3. 擬擬Newton方法方法目的目的不計算不計算Hessian陣，克服病態(tài)、不正定、計陣，克服病態(tài)、不正定、計算復雜等缺陷，同時保持收斂較快的優(yōu)點算復雜等缺陷，同時保持收斂較快的優(yōu)點回顧解方程組回顧解方程組 F(x)=0的擬牛頓法的擬牛頓法)()()(11kkkkkkxFxFxxAA 滿足使kkkkAAAA計

12、算用迭代方法11思路思路)()(11kkkkxFxFxx)()(11kkkkxFAxx)(xfMinx優(yōu)化問題優(yōu)化問題)()(xFxf相當fGfxFxf222),()(代替陣不一定正定，構造正定相當3. 擬擬Newton方法方法( (續(xù)續(xù)) )kkkkkkHHHGGG11或迭代公式構造按照按照 擬牛頓條件擬牛頓條件：kkkkkkfHxfxG11或)()(,11kkkkkkxfxffxxx記)(1kkkkxfHxx設在第k步, Gk已得到, Hk=(Gk)-1,可計算)(1112kkkkxfHxx于是有3.1 Davidon-Fletcher-Powell(DFP)公式公式kkTkkTkkkkT

13、kTkkkfHfHffHfxxxH)()()()(3. 擬擬Newton方法方法(續(xù)續(xù))kTkTkkkkTkkkTkTkkkTkkkTkkfxfxGGxffxfffxxGxG)()()()()()()(1 (3.2 Broyden-Fletcher-Goldfarb-Shanno(BFGS)公式公式kkTkkTkkkkTkTkkkxGxGxxGxfffG)()()()(kTkTkkkkTkkkTkTkkkTkkkTkkfxxfHHfxfxxxfxfHfH)()()()()()()(1 (搜索步長的確定搜索步長的確定線性搜索線性搜索(Line Search)確定步長方法確定步長方法問題問題給定x

14、k和方向dk, 確定步長 k, 使得)(minkkdxf優(yōu)化優(yōu)化算法算法黃金分割黃金分割(0.618)法、法、Fibonacci法、法、Newton切線法、割線法、切線法、割線法、2次或次或3次插值法次插值法等等一維優(yōu)化問題0kddf0)()(kkkTkdxfd非線性最小二乘(Least Square)擬合問題問題給定給定(t i, y i), i=1,n, 擬合一個函數(shù)擬合一個函數(shù)y=f(t, x)， 其中其中x為待定的參數(shù)向量為待定的參數(shù)向量, f 對對x非線性。非線性。優(yōu)化模型優(yōu)化模型)()()(minxrxrxRTx),()(xtfyxriiiTnxrxrxr)(),()(1記誤差記誤

15、差（有時為記號方便，目標乘以（有時為記號方便，目標乘以1/2）根據(jù)目標函數(shù)是根據(jù)目標函數(shù)是 r(x) 的二次函數(shù)的特點的二次函數(shù)的特點構造簡單算法構造簡單算法非線性最小二乘擬合非線性最小二乘擬合)()(xrxJRTSxJxJRT2)()(22mmlkiixxrxr)()(22討論討論 牛頓法要計算牛頓法要計算Hessian矩陣，其中矩陣，其中S計算量大計算量大 若若 f 對對x線性線性, 則化為線性最小二乘擬合則化為線性最小二乘擬合, 此時此時S=0特定算法考慮如何特定算法考慮如何忽略或近似矩陣忽略或近似矩陣S SmnjixrxJ)/()(記)(xr的雅各比陣為)()()(xrxrxRT)()

16、(21xrxrSiniiGauss-Newton算法算法(GN)：忽略矩陣忽略矩陣S)()()()(kkkkTkxrxJdxJxJf用用R代替，下降方向代替，下降方向dk滿足滿足GN算法收斂性特點算法收斂性特點收斂性依賴收斂性依賴 f 對對 x 的線性程度的線性程度, 及偏差及偏差r的大小的大小(“小殘量問小殘量問題題”)非線性最小二乘擬合非線性最小二乘擬合)()(2kkxfdxf牛頓方程牛頓方程 )()(xrxJRT)()(22xJxJRT)()(21xrxrSiniiG-N算法修正算法修正)()()()(kkkkkTkxrxJdIxJxJ其中 k0為修正參數(shù).dk位于位于GN方向方向( k

17、很小很小)和和負梯度方向負梯度方向( k很大很大)之間之間LM算法特點算法特點非線性最小二乘擬合非線性最小二乘擬合)()()()(kkkkTkxrxJdxJxJ防止 JTJ 出現(xiàn)病態(tài) LM算法算法: Levenbery (1944), Marquardt (1963)若若GN方向滿足約束方向滿足約束(*), k=0; 否則否則 k0先確定最大步長先確定最大步長hk，再確定方向，再確定方向dk信賴域算法特點信賴域算法特點LM算法的另一種觀點算法的另一種觀點(*) |. .|)()(|min22kkkkkhxxtsxxxJxr)()()()(kkkkkTkxrxJdIxJxJ其解為其解為Mor (

18、1978) 給出給出LM算法一種可靠和有效的實現(xiàn)算法一種可靠和有效的實現(xiàn)信賴域（置信域，信賴域（置信域，Trust Region）優(yōu)化工具箱基本用法：基本用法：x=fminbnd(f,lb,ub)x=fminunc(f,x0)x=fminunc(f,x0,options,P1,P2,.)x=fminsearch(f,x0,options,P1,P2,.)f.m f(x)的的m文件名文件名x0 初始點初始點; x 最優(yōu)解最優(yōu)解P1,P2, 傳給傳給fun的參數(shù)的參數(shù)中間輸入項缺省用中間輸入項缺省用 占位占位nxRxxfMin),(模型：2min 222babyax其中：例8 , 1 x)x3si

19、nx(min 1：例非線性最小二乘法非線性最小二乘法)()()(minxrxrxRTubxlb),()(xtfyxriiiTnxrxrxr)(),()(1x,resnorm,res,exitf,out,lambda,jacob=lsqnonlin(fun,x0,lb,ub,options,P1,P2,)輸入的用法與輸入的用法與fminunc類似，但注意：類似，但注意：f.mf(x)的的m文件名文件名: function y=f(x,t)fun.mr(x)的的m文件名文件名: function r=fun(x,t,y)輸出輸出 resnorm=r(x)T*r(x), res=r(x)(誤差向量誤

20、差向量)x,resnorm,res,exitf,out,lambda,jacob=lsqcurvefit(f,x0,t,y,lb,ub,options,P1,P2,)非線性最小二乘法非線性最小二乘法)()()(minxrxrxRTubxlbt0.250.511.523468c19.21 18.15 15.36 14.10 12.899.327.455.243.01ktretc)(例例3 用下面一組數(shù)據(jù)擬合系數(shù)r, k : 控制精度，觀察中間結果，控制迭代次數(shù)等optimset /顯示控制參數(shù)顯示控制參數(shù)optimset optfun /顯示顯示optfun的控制參數(shù)的控制參數(shù)opt=optim

21、set /控制參數(shù)設為控制參數(shù)設為(即缺省值即缺省值)opt=optimset(optfun)/optfun控制參數(shù)設缺省值控制參數(shù)設缺省值opt=optimset(par1,val1,par2,val2,.)opt=optimset(oldopts,par1,val1,.)opt=optimset(oldopts,newopts)val=optimget(opt,par1,par2,)val=optimget(opt,par1,par2, default)options 控制參數(shù)設定控制參數(shù)設定/獲取獲取: optimset; optimgetDiagnostics on | off /是否

22、顯示診斷信息是否顯示診斷信息Display off | iter | final | notify /顯示信息的級別顯示信息的級別GradObj on | off/是否采用分析梯度是否采用分析梯度Jacobian on | off /采用分析采用分析Jacob陣（用于約束優(yōu)化中）陣（用于約束優(yōu)化中）LargeScale on | off /是否采用是否采用大規(guī)模算法大規(guī)模算法MaxFunEvals 最大函數(shù)調用次數(shù)最大函數(shù)調用次數(shù)MaxIter 最大迭代次數(shù)最大迭代次數(shù)TolCon 約束的控制精度（用于約束優(yōu)化中）約束的控制精度（用于約束優(yōu)化中）TolFun 函數(shù)值的控制精度函數(shù)值的控制精度T

23、olX 解的控制精度解的控制精度主要控制參數(shù)主要控制參數(shù)( (對大規(guī)模對大規(guī)模/ /中等規(guī)模算法均有效中等規(guī)模算法均有效) )最一般的輸出形式最一般的輸出形式x,f,exitflag,out,grad,hess=fminunc(.)更多輸出：最優(yōu)值等更多輸出：最優(yōu)值等f 目標函數(shù)值目標函數(shù)值exitflag 0收斂收斂,0達到函數(shù)或迭代次數(shù)達到函數(shù)或迭代次數(shù), 0不收斂不收斂Output iterations 實際迭代次數(shù)實際迭代次數(shù) funcCount 實際函數(shù)調用次數(shù)實際函數(shù)調用次數(shù) algorithm 實際采用的算法實際采用的算法 cgiterations 實際實際PCG迭代次數(shù)（大規(guī)模

24、算法用）迭代次數(shù)（大規(guī)模算法用） stepsize 最后迭代步長（中等規(guī)模算法用）最后迭代步長（中等規(guī)模算法用） firstorderopt 一階最優(yōu)條件（梯度的范數(shù)）一階最優(yōu)條件（梯度的范數(shù)） message 收斂信息等收斂信息等grad 目標函數(shù)的梯度目標函數(shù)的梯度hess 目標函數(shù)的目標函數(shù)的Hessian矩陣矩陣用用optimset控制參數(shù)控制參數(shù)選擇選擇算法選擇：算法選擇：optimset中參數(shù)控制中參數(shù)控制LargeScale on | off （on為缺省）為缺省） 搜索步長的算法選擇搜索步長的算法選擇(缺省：缺省：quadcubic)系統(tǒng)缺省采用大規(guī)模算法系統(tǒng)缺省采用大規(guī)模算法

25、(如置信域法，如果可能如置信域法，如果可能)當設為當設為off時缺省：時缺省：BFGS、混合、混合2, 3次多項式插次多項式插值值LineSearchType=quadcubic 混合混合2, 3次多項式插值次多項式插值LineSearchType=cubicpoly 3次多項式插值次多項式插值HessUpdate = dfp （DFP算法）算法）HessUpdate = steepdesc（最速下降算法）最速下降算法）搜索方向的算法選擇搜索方向的算法選擇(缺省：缺省：BFGS)計算結果計算結果 方向 步長 最優(yōu)解x 最優(yōu)值f n BFGS 2,3次 (9.9997e-001 9.9994e-

26、001) 1.0944e-009 165 DFP 2,3次 (9.9997e-001 9.9994e-001) 2.1173e-009 165 steep 2,3次 (-8.0263e-001 6.5064e-001) 3.2536e+000 1003 BFGS 3次 (1.0001e+000 1.0003e+000) 3.1432e-008 162 DFP 3次 (9.0468e-001 8.1572e-001) 9.8270e-003 495steep 3次 (-1.1831e+000 1.4056e+000) 4.7695e+000 1002222)1 ()(100),(min. 5xx

27、yyxf例精確解：x=y=1, f(x,y)=0無約束優(yōu)化無約束優(yōu)化nxRxxfMin),(模型：采用分析梯度：采用分析梯度：opt 中中 GradObj=on 時：時：x=fminunc(fun,x0,opt,)fun.m中還要有中還要有一般形式一般形式 function f,g=fun(x)(xf222)1 ()(100),(min. 5xxyyxf續(xù)例)(200)1 (2)(40022xyxxyxf算出 方向 步長 最優(yōu)解x 最優(yōu)值f n BFGS 2,3次 (1.0001e+000 1.0002e+000) 8.9857-009 105 DFP 2,3次 (1.0001e+000 1.

28、0003e+000) 2.2499e-008 109 BFGS 3次 (1.0000e+000 1.0001e+000) 2.7164e-009 53 DFP 3次 (7.4711e-001 5.5212e-001) 6.7609e-002 144計算結果計算結果與不用分析梯度的結果比較與不用分析梯度的結果比較 方向 步長 最優(yōu)解x 最優(yōu)值f n BFGS 2,3次 (9.9997e-001 9.9994e-001) 1.0944e-009 165 DFP 2,3次 (9.9997e-001 9.9994e-001) 2.1173e-009 165 BFGS 3次 (1.0001e+000 1

29、.0003e+000) 3.1432e-008 162 DFP 3次 (9.0468e-001 8.1572e-001) 9.8270e-003 495x0*O的圖形222)1 ()(100),(xxyyxfRosenbrock 函數(shù)（香蕉函數(shù)）算法選擇算法選擇：BFGS公式，混合公式，混合2，3次插值，一般較好。次插值，一般較好。幾個值得注意的問題幾個值得注意的問題精度控制精度控制：對迭代次數(shù)有重大影響，應適當選擇。：對迭代次數(shù)有重大影響，應適當選擇。梯度函數(shù)梯度函數(shù)：利用分析梯度：利用分析梯度可能可能改進算法的性能改進算法的性能改變初始值改變初始值 由一個初值出發(fā)通常得到局部最優(yōu)解，由一個

30、初值出發(fā)通常得到局部最優(yōu)解，如果函數(shù)存在多個局部最優(yōu)，只有改變初值，對局如果函數(shù)存在多個局部最優(yōu)，只有改變初值，對局部最優(yōu)進行比較，才有部最優(yōu)進行比較，才有可能可能得到全局最優(yōu)解。得到全局最優(yōu)解。其他算法選擇其他算法選擇： （詳細用法請查閱（詳細用法請查閱help文檔）文檔）高度非線性或不連續(xù)時可用程序高度非線性或不連續(xù)時可用程序 fminsearch(fun,x0)單變量時可用程序單變量時可用程序 fminbnd(fun,v1,v2)非線性最小二乘法：使用分析導數(shù)非線性最小二乘法：使用分析導數(shù))()()(minxrxrxRTubxlb),()(xtfyxriiiTnxrxrxr)(),()(

31、1x,resnorm,res,exitf,out,lambda,jacob=lsqnonlin(fun,x0,lb,ub,options,P1,P2,)輸入的用法與輸入的用法與fminunc類似，但注意：類似，但注意：fun.m r(x) 的的m文件名，文件名，Jacobian=on時含有導數(shù)信息時含有導數(shù)信息一般形式：一般形式：function r,J = fun(x)(xr缺省：大規(guī)模算法（缺省：大規(guī)模算法（LargeScale = on ）當當LargeScale = off ：缺省：缺省： Levenberg-Marquardt算法；算法；LevenbergMarquardt=off:

32、Gauss-Newton法法一維搜索（線搜索）步長選擇與一維搜索（線搜索）步長選擇與fminunc中類似中類似非線性最小二乘法：算法選擇非線性最小二乘法：算法選擇)()()(minxrxrxRTubxlb實例實例1 產(chǎn)銷量安排產(chǎn)銷量安排28010022 .01 .0122211211baaaa已知已知數(shù)據(jù)數(shù)據(jù))3020(,02. 0015. 0,10030TTTcr22222212121112121111)()()(min2211xcerxaxabxcerxaxabxfxx, 2 , 1, 0,212211iaabxaxabpiiiiiii, 2 , 1, 0,icrcerqiiiixiiii22211121,)()(),(max21xqpxqpxxzxx原問題原問題x = 23.9025 62.4977 y = 6.4135e+003即甲產(chǎn)量為23.9025，乙產(chǎn)量為62.4977，最大利潤為6413.5命令和最優(yōu)解命令和最優(yōu)解x=fminunc(f, x0)shili0701.m702/502/22221111abxabx，其