實(shí)驗(yàn)一工業(yè)生產(chǎn)的產(chǎn)量問題

1、實(shí)驗(yàn)七 生產(chǎn)計(jì)劃中的產(chǎn)量問題實(shí)驗(yàn)?zāi)康摹?介紹無約束最優(yōu)化方法的一些基本概念。2了解幾種常見的無約束優(yōu)化問題的求解方法，如迭代算法、最速下降法(梯度法)、牛頓法(Newt on)、擬牛頓法。3學(xué)習(xí)掌握用MATLAB優(yōu)化工具箱中的命令來求解無約束優(yōu)化問題。實(shí)驗(yàn)內(nèi)容】某公司生產(chǎn)一種產(chǎn)品有甲、 乙兩個(gè)品牌， 試討論產(chǎn)銷平衡下的最大利潤。 所謂產(chǎn)銷平衡 指公司的產(chǎn)量等于市場上的銷量。 利潤既取決于銷量和單件價(jià)格， 也依賴于產(chǎn)量和單件成本。按照市場規(guī)律，甲種品牌的價(jià)格 pi固然會(huì)隨其銷量Xi的增長而降低；同時(shí)乙品牌銷量X2的 增長也會(huì)使甲的價(jià)格有稍微下降，根據(jù)實(shí)際情況，可以確定價(jià)格與銷量成線性關(guān)系，即p1

2、 = 300 2.35 x1 0.09 x2乙的價(jià)格卩2遵循同樣的規(guī)律，有p2 = 480 0.14 x1 2.98 x2甲品牌的成本會(huì)隨著其產(chǎn)量的增長而降低，按實(shí)際情況可假設(shè)為負(fù)指數(shù)關(guān)系，即有q1 = 38 e 0.023X1 + 1 16乙品牌的成本q2遵循同樣的規(guī)律，有0.018 x2q2 = 94 e2 + 145試確定甲、乙兩種品牌的產(chǎn)量，使公司獲得的總利潤最大。【實(shí)驗(yàn)準(zhǔn)備】 無約束最優(yōu)化方法是指在沒有約束條件限制下， 求多變量實(shí)值函數(shù)極值的方法。 無約束 最優(yōu)化問題的數(shù)學(xué)表達(dá)式為min f (x) , x =( x1, x2，Xn) Rn(1)一般f為非線性函數(shù)，x是n維實(shí)變量，實(shí)

3、際上這是一個(gè)多元函數(shù)無條件極值 問題。由于一個(gè)求極大值問題， 可以添加負(fù)號(hào)的方式轉(zhuǎn)化為求極小值問題， 因此通常只討論求極小 值問題。應(yīng)該注意的是，極值問題的解，即極值點(diǎn)，都是 局部最優(yōu)解 ，全局最優(yōu)解 只能從局 部最優(yōu)解的比較中得到。如何求解無約束最優(yōu)化問題( 1)的最優(yōu)解呢？一般是采用迭代方法，即先選擇一個(gè)初 始點(diǎn)， 再尋找該點(diǎn)處的下降方向，我們稱為搜索方向。在該方向上求極小點(diǎn)，得到一個(gè)新的 點(diǎn)。這個(gè)新的點(diǎn)要優(yōu)于原來的點(diǎn)， 即新點(diǎn)處的目標(biāo)函數(shù)值小于原來點(diǎn)處的目標(biāo)函數(shù)值。然后在新點(diǎn)處再尋找下降方向和該方向上求極小點(diǎn)，如此下去，最終得到最優(yōu)點(diǎn)。因此， 求解無約束最優(yōu)化問題需解決兩個(gè)問題： 一是在

4、某些方向上的一維極小點(diǎn)，我們也稱為一維搜索； 另一個(gè)是尋找某些點(diǎn)處的下降方向， 這是無約束最優(yōu)化方法中最重要的一 個(gè)問題。我們先了解第一個(gè)問題最常用的搜索方法。1. 求一維極小的二點(diǎn)二次插值方法設(shè)d(k)是點(diǎn)x(k)處的一個(gè)搜索方向，要在該方向上尋優(yōu)問題，轉(zhuǎn)化為求一維函數(shù)()=f ( x(k) + d(k)(2)的求極值問題。最常用的一維搜索方法是插值法，即用某些點(diǎn)的函數(shù)值(或?qū)?shù)值)構(gòu)造插值函數(shù)，用插值函數(shù)的極小點(diǎn)來近似函數(shù)()的極小點(diǎn)。這里介紹一種有效的插值方法，稱為二點(diǎn)二次插值方法，即用二點(diǎn)處的函數(shù)值和一個(gè)點(diǎn)處的導(dǎo)數(shù)值構(gòu)造二次函數(shù)，反復(fù)用二次函數(shù)的極小點(diǎn)來逼近函數(shù)()的極小點(diǎn)。算法1二點(diǎn)

5、二次插值法(1) 取初始點(diǎn) 1 (如(如 =0.1 )，置精度要求(2)則置(3)(4)計(jì)算如果(5)計(jì)算(6)，則停止計(jì)算(1 = 0),初始步長 (如，并計(jì)算1 =( 1),;否則置 =2 = ( 2 )2 1 ),則置=21( 2 J'1 ( 2 1)作為極小點(diǎn))；和步長縮減因子(1)，轉(zhuǎn)(3)否則置(0, 1)= ，轉(zhuǎn)2.最速下降法前面介紹了一維搜索的二點(diǎn)二次插值方法, 來看看兩個(gè)概念。2)F面討論如何選擇搜索方向的問題,我們先定義1稱n維向量(X1X2)T為函數(shù)f(x)在x處的梯度，記為 f(x)Xn f (X)=(X1定義2設(shè)d是任意的單位向量，若極限 limn 0Xnf

6、(X ")一空存在，則稱該極限為(3)函數(shù)f (X)在X處沿方向d的一階方向?qū)?shù)，簡稱為方向?qū)?shù)，記為f(X)，即d_f(x) = lim f(X ad) f(X)(4dn 0a最速下降法的基本思想：選取一點(diǎn)x1作為初始點(diǎn)，計(jì)算該點(diǎn)的梯度f(x)，求該點(diǎn)處的最速下降方向，即令 d1 = f (X1)，再沿d1方向前進(jìn)，尋找該方向上的極小點(diǎn)，得 到點(diǎn)x2，再計(jì)算 f(x2)，令d2 = f(x2)，沿d2方向前進(jìn)，得到點(diǎn)x3，如此下去 具體算法如下：算法2最速下降法(1) 取初始點(diǎn)x1，置精度要求為，置k = 1(2) 若| f (xk) | v ，則停止計(jì)算(xk為無約束問題的最優(yōu)解

7、)；否則置dk = f (xk)(3) 一維搜索，求一維問題min (a) = f ( xk ad k)ak 1kk得 ak，置 x = x ad 。(4) 置k = k +1，轉(zhuǎn)步驟(2)在算法中，為精度要求，即當(dāng)梯度接近于0時(shí)，我們就認(rèn)為達(dá)到極小點(diǎn)，終止計(jì)算。這樣做的目的是避免算法產(chǎn)生死循環(huán)，算法中的一維搜索可用算法1來求解。最速下降法是一種最基本的算法，它在最優(yōu)化方法中占有重要地位。其優(yōu)點(diǎn)是工作量小， 存儲(chǔ)變量較少，初始點(diǎn)要求不高；缺點(diǎn)是收斂慢，效率不高，有時(shí)達(dá)不到最優(yōu)解。下面介紹一種簡單而直觀的方法 Newton法。3. 牛頓法22 f將f (x)的二階導(dǎo)數(shù)記作f(X)=()是一個(gè)n

8、n矩陣，稱黑塞(Hessian )XiXj矩陣(簡記H(X)。當(dāng)f (x)二階偏導(dǎo)數(shù)連續(xù)時(shí)，H(X)是對(duì)稱矩陣。Newton法的基本思想是用一個(gè)二次函數(shù)去近似目標(biāo)函數(shù)f(x)，然后精確地求出這個(gè)二次函數(shù)的極小點(diǎn)，并把這個(gè)極小點(diǎn)作為所示函數(shù)的極小點(diǎn)x的近似值。k 1*k設(shè)x 是f(x)極小點(diǎn)x的第k 1次近似，將f (X)在x點(diǎn)作二階Taylor展開，略去 高于二階的項(xiàng)，用dk = xk 1 xk得到f (xk 1) = f (xk dk)f(xk) + f (xk) dk + 1 d 2 f (xk) d k(5)為了求dk使f(xk dk)最小，只需將(5)式右端對(duì)dk求導(dǎo)，并令其為零，可得

9、 f (xk) + 2 f(xk) d k = 0(6)稱為牛頓方程。當(dāng)黑塞矩陣2 f (xk)正定時(shí)，可解得dk 一 ( 2f(xk)-1 f(xk)( 7)稱為牛頓方向。其迭代公式為k 1 k2 k -1kx = x - ( f(x ) f(x )(8)4. 擬牛頓法牛頓法的優(yōu)點(diǎn)是在接近最優(yōu)解時(shí)具有2階收斂性，缺點(diǎn)是計(jì)算復(fù)雜，而且會(huì)出現(xiàn)病態(tài)、不正定等情況。為克服這些問題，同時(shí)又保持較快的收斂，禾U用第k和第k 1步得到的xk、k 1kk 1k 12kx 、 f(x ) . f (x )根據(jù)擬牛頓條件構(gòu)造一正定矩陣G近似代替f(x )或H k 1近似代替(2 f (xk 1)-1，得到的迭代

10、算法稱為擬牛頓法。其中有兩個(gè)常用的著名的迭代公式為BFGS公式和DFP公式，可參見其他最優(yōu)化方法書籍或MATLAB幫助。5. MATLAB求解無約束優(yōu)化問題的主要命令(1) MATLAB5.2及以下版本使用的命令x = fmin( ' fun ' , x1 , x2 )求一元函數(shù) y = f( x )在x1 , x2 內(nèi)的極小值x = fmin( ' fun ' , x1 , x2 , optio ns )同上，參數(shù) optio ns 的定義由表 1給出x , opti ons = fmin( .)同上，同時(shí)返回參數(shù)optio ns 的值fmin函數(shù)采用黃金分割

11、法和拋物線插值法，fun可直接用函數(shù)表達(dá)式表示，也可以是用M文件定義的函數(shù)名。建立函數(shù)形式的M-文件格式如下，文件名fun.m (般M函數(shù)文件名取函數(shù)名)function f = fun( x )f = f( x )x = fmins( ' fun ' , x0 )求多元函數(shù)(1)以x0為迭代初值的局部極小值x = fmins( ' fun ' , x0 , options )同上，參數(shù) options 的定義由表 1 給出x , optio ns = fmin s(.)同上，同時(shí)返回參數(shù) optio ns 的值fmins函數(shù)采用Nelder - Meade單純

12、形搜索法，fun可直接用函數(shù)表達(dá)式表示，也 可以是用M文件定義的函數(shù)名。x = fminu( ' fun ' , xO )x = fminu( ' fun ' , x0 , opti ons )x , opti ons = fminu( ' fun ' , x0 ,.)的值fminu函數(shù)為無約束優(yōu)化提供了三種算法，由 提供了二種算法，由optio ns (7)控制。這里,求多元函數(shù)(1)以x0為迭代初值的局部極小值同上，參數(shù)options 的定義由表1給出同時(shí)返回參數(shù)optio ns同上,opti ons (6) fminu必須先用控制，為步長一

13、維搜索 M-文件定義函數(shù) fun。(2) MATLAB5.3以上版本使用的命令fminbnd替代fmin求解一元函數(shù)極值，使用格式、搜索算法與之相同，x , Fval =fminbnd( ' fun ' , x0 ,.)同時(shí)返回解x處的函數(shù)值，而不是參數(shù)options的返回值fminsearch替代fmins求解多元函數(shù)(1)極值，使用格式、搜索算法與之相同fminunc替代fminu求解多元函數(shù)(1)極值，使用格式、搜索算法與之相同有關(guān)上述命令的詳細(xì)信息和使用方式可在幫助文件中了解，或在命令框里輸入help“命令名”查閱。在MATLAB5.3以上的各種最新版本中fmin、fm

14、ins和fminu命令仍然有用，但在MATLAB勺未來版本將可能刪除這些命令。(3) 命令中參數(shù)optio ns的有關(guān)定義在大多數(shù)MATLAB化命令函數(shù)中有一個(gè)控制參數(shù)options，它是一個(gè)有18個(gè)分量的向量，包含了在優(yōu)化程序中要用到的參數(shù)，以便在計(jì)算最優(yōu)值時(shí)控制精度要求、輸出形式、搜索算法、迭代次數(shù)、步長等等是。在對(duì)優(yōu)化命令函數(shù)的第一次調(diào)用時(shí)，向量options會(huì)自動(dòng)使用缺省值；當(dāng)然在調(diào)用前，可以對(duì)optio ns的某些分量進(jìn)行賦值，以達(dá)到控制要求。表1 參數(shù)options各分量說明序號(hào)功能定義分量說明缺省值含義1輸出形式options(1)=1 ，有中間結(jié)果輸出 options(1)=

15、1,給出警告信息0，無中間結(jié)果輸出2解x(k)的精度(k)options(2)設(shè)置x 的精度|x(k1)(k) xx(k)|< 1e 43函數(shù)f (k)的精度options(3)設(shè)置f (k)的精度1)f (k)f (k)-< 1e 44函數(shù)g(k)的精度options(4)設(shè)置g (k)的精度，由 命令 constr、minmax使用終止判 據(jù)g(k1) g(k)-W 彳a 7(k) g5主要算法options(5)=1 ，高斯-牛頓法0，LM方法6搜索方向算法options(6)=1 ,擬牛頓法DFP公式 options(6)=2，最速下降法0，擬牛頓法的BFGS公式7步長一維

16、搜索算法options=1，三次多項(xiàng)式插值0,混合的二次和三次多項(xiàng)式8函數(shù)值輸岀options©)輸出極值點(diǎn)處的函數(shù) 值9梯度檢查options(9)=1，最初幾步，梯度將與有限微分計(jì)算的結(jié)果比較0,不作檢查10函數(shù)計(jì)算次數(shù)option(10)輸出函數(shù)計(jì)算次數(shù)11梯度計(jì)算次數(shù)option(11)輸出梯度計(jì)算次數(shù)12約束梯度計(jì)算次數(shù)option(12)輸出約束梯度計(jì)算次 數(shù)13等式約束的個(gè)數(shù)option(13)確定等式約束數(shù)目0，等式約束個(gè)數(shù)為014最大迭代次數(shù)option(14)輸入迭代的最大次數(shù)100 n，n為自變量個(gè)數(shù)200 n (fmins )500 n (fmin )15目標(biāo)

17、優(yōu)化由 attgoal 命令使用，option(15) 輸入須達(dá)目的的目標(biāo)個(gè)數(shù)016差分步長最小值用差分計(jì)算梯度時(shí)步長的下限1e 817差分步長最大值用差分計(jì)算梯度時(shí)步長的上限0.118步長參數(shù)，一般取1或更小【實(shí)驗(yàn)方法與步驟】1. MATLAB化工具箱中解無約束問題的命令和參數(shù)options的基本用法下面用例題來予以說明例 1 求解 min f(x) = ex 5 x , 1< x< 2輸入命令：>> x=fmin( 'exp(x) - 5*x' , 1 , 2 ) , f=exp( x ) - 5*xx =1.6094f =-3.0472同時(shí)，該問題

18、可以用 fminbnd求解，得到的解答一樣輸入命令：>> x , fval=fmi nbn d( 'exp( x ) - 5*x' , 1 , 2 )x =1.6094fval =-3.0472例2用擬牛頓法的DFP公式求解min f (x) = 4x" + x； x3x2，初值為(1, 2)，要 求精度為1e 7,給出函數(shù)計(jì)算次數(shù)以及函數(shù)值首先建立M-文件，文件名取函數(shù)名fun.mfunction f = fun( x )f = 4 * x( 1 )A2 + x( 2 )A2 - x( 1 )A3 * x( 2 )輸入命令：>> x0=1,2

19、;>> opt( 2 )=1e-6;%設(shè)置自變量要求的精度>> opt( 3 )=1e-6;%設(shè)置函數(shù)所要求的精度>> opt( 6 )=1;%搜索算法用擬牛頓法的DFP公式求解>> x, opt=fminu( ' fun ', x0, opt);%返回參數(shù) options 的值給向量 opt可得到結(jié)果x =1.0e-008 *-0.34930.7566參數(shù)options的值返回給向量opt，全部元素略。>> n=opt( 10 );%函數(shù)計(jì)算次數(shù)賦值給nn =35>> y=opt( 8 );%在極值點(diǎn)處函

20、數(shù)值賦值給yy =1.0606e-0162. 引例問題的分析與求解由題設(shè)可知，甲品牌產(chǎn)品單件獲利為p1 q1，乙品牌產(chǎn)品單件獲利為P2 q2，由產(chǎn) 銷平衡原理，所有產(chǎn)品的銷量即為產(chǎn)量，則甲、乙兩種產(chǎn)品總獲利為Z(X1 ,X2)=( Pi qi) Xi +( P2 q2) x?容易看出，原問題實(shí)際上轉(zhuǎn)化為求二元函數(shù)z(x1 ,X2)的極大值，為用MATLAB優(yōu)化工具箱中的fminunc求解，需將其轉(zhuǎn)化為求函數(shù)y(x1, x2) = z(x1, x2)的最小值。為確定初始值，先忽略成本，并令p1價(jià)格中x2項(xiàng)的較小系數(shù)0.09和p2中x1項(xiàng)較小的 系數(shù)0.14等于零(因?yàn)樗鼈儗?duì)價(jià)格的作用比較微弱，暫時(shí)可忽略不計(jì))，則確定初值問題轉(zhuǎn)化為求Z(Xi,X2)=( 300 2.35 Xi) Xi +( 480 2.98 X?) X：的極值，很容易可以求得xi =300 = 63.83 , x2 = 480= 80.54，我們用它作為原2 2.3522.98問