建模仿真與優(yōu)化設計.ppt

,建模仿真與優(yōu)化設計 第一部分 中北大學 張保成,掌握優(yōu)化問題的數學描述方法； 2. 熟練掌握常用優(yōu)化算法。,一 優(yōu)化模型的一般意義 二 非線性規(guī)劃建模無約束最優(yōu)化 非線性規(guī)劃建模有約束非線性規(guī)劃 四 連續(xù)結構體建模與優(yōu)化設計,學 習 內 容,目 的,（一）優(yōu)化模型的數學描述,一 優(yōu)化模型的一般意義,“受約束于”之意,（二）優(yōu)化模型的分類,1.根據是否存在約束條件 有約束問題和無約束問題。 2.根據設計變量的性質 靜態(tài)問題和動態(tài)問題。,3.根據目標函數和約束條件表達式的性質 線性規(guī)劃，非線性規(guī)劃，二次規(guī)劃，多目標規(guī)劃等。,（1）非線性規(guī)劃 目標函數和約束條件中，至少有一個非線性函數。,（2）線性規(guī)劃（LP） 目標函數和所有的約束條件都是設計變量 的線性函數。,（3）二次規(guī)劃問題 目標函數為二次函數，約束條件為線性約束,5. 根據變量具有確定值還是隨機值 確定規(guī)劃和隨機規(guī)劃。,4. 根據設計變量的允許值,整數規(guī)劃（0-1規(guī)劃）和實數規(guī)劃。,（三）建立優(yōu)化模型的一般步驟,1.確定設計變量和目標變量； 2.確定目標函數的表達式； 3.尋找約束條件。,二、優(yōu)化求解的數學基礎,函數的梯度,泰勒展開,二階導數矩陣,矢量的概念、運算和點積,矩陣的運算和逆矩陣,（一）方向導數,和 分別是 函數 f( x1 , x2 )在x0點處沿坐標軸x1和x2方向的變化率,故函數 f( x1 , x2 )在x0(x10 , x20)點處沿某一方向S的變化率為： 稱為該函數沿此方向的方向導數 偏導數可以看作是函數沿坐標軸 方向的方向導數，并有 （二）梯度 二元函數在點x0的梯度是由函數在該點的各一階偏導數組成的向量。即 :,二、優(yōu)化求解的數學基礎,設S方向單位向量 則有 函數的梯度具有以下性質： （1）函數在一點的梯度是一個向量。梯度的方向是該點函數值上升最快的方向，與梯度相反的方向是該點函數值下降的最快的方向，梯度的大小就是它的模長。 （2）一點的梯度方向是與過該點的等值線或等值面的切線或切平面相垂直的方向，或者說是該點等值線或等值面的法線方向。 （3）梯度是函數在一點鄰域內局部形態(tài)的描述。在一點上升得快的方向，離開該領域后就不一定上升得快，甚至可能下降。,二、優(yōu)化求解的數學基礎,（三）泰勒展開 為了便于數學問題的分析和求解，往往需要將一個復雜的非線性函數簡化成線性函數或二次函數。簡化的方法可以采用泰勒展開式。 由高等數學可知，一元函數 f(x)若在點x0的鄰域內n階可導，則函數可在該點鄰域內作泰勒展開： 二元函數 f(x)在點x0(x10 , x20)也可以作泰勒展開，展開式一般取前三項，即：,二、優(yōu)化求解的數學基礎,將上式寫為矩陣形式 其中，G(x0)稱為函數 f( x1 , x2 )在點x0處的二階導數矩陣或海賽矩陣。,二、優(yōu)化求解的數學基礎,在優(yōu)化計算中，當某點附近的函數值采用泰勒展開式作近似表達時，研究該點鄰域的極值問題需要分析二次型函數是否正定。當對任何非零向量x使 則二次型函數正定，G為正定矩陣,二、優(yōu)化求解的數學基礎,（四）二次函數 當將函數的泰勒展開式取到二次項時得到二次函數形式。優(yōu)化計算經常把目標函數表示為二次函數以使問題分析得以簡化。在線性代數中將二次齊次函數稱作二次型，其矩陣形式 在優(yōu)化計算中，當某點附近的函數值采用泰勒展開式作近似表達時，研究該點鄰域的極值問題需要分析二次型函數是否正定。當對任何非零向量x使 則二次型函數正定，G為正定矩陣,二、優(yōu)化求解的數學基礎,對于一般二次函數 矩陣有正定和負定之分。對于所有非零向量： （1）若有 ，則稱矩陣是正定的； （2）若有 ，則稱矩陣是半正定的； （3）若有 ，則稱矩陣是負定的； （4）若有 ，則稱矩陣是半負定的； （5）若有 ，則稱矩陣是不定的,二、優(yōu)化求解的數學基礎,可以證明，正定二次函數具有以下性質： （1）正定二次函數的等值線或等值面是一簇同心圓或同心橢球。橢圓簇或橢球簇的中心就是該二次函數的極小點。 （2）非正定二次函數在極小點附近的等值線或等值面近似于橢圓或橢球。,二、優(yōu)化求解的數學基礎,（五）無約束優(yōu)化問題的極值條件 二次函數 f(x1 , x2)在x0取得極值的必要條件為 充分條件為：該點處海賽矩陣正定，即,二、優(yōu)化求解的數學基礎,正定：各階主子式均大于零。,現代設計方法概論課程教案,（六） 下降迭代算法及其收斂性,無約束最優(yōu)化問題求優(yōu)過程的求解方法大致分為兩類。 （1）解析法,二、優(yōu)化求解的數學基礎,現代設計方法概論課程教案,（2）數值迭代計算,（六） 下降迭代算法及其收斂性,二、優(yōu)化求解的數學基礎,三、優(yōu)化設計的迭代計算,（一）優(yōu)化問題的求解方法,1、優(yōu)化問題的本質,優(yōu)化問題的本質求極值的數學問題。,2、優(yōu)化問題的求解方法,理論上,解析法,數值計算法,能求解嗎？,由于實際優(yōu)化數學模型的目標函數及約束函數往往是非線性的,解析法求解非常困難，甚至無法實現,（二）數值計算法的迭代方法,1、數值計算法的數學基礎,計算方法,2、數值計算法的迭代方法,數值計算法可以較好地解決這類問題,三、優(yōu)化設計的迭代計算,從目標函數出發(fā),構造一種使目標函數值逐次下降的數值計算方法,利用計算機進行反復迭代運算,一步步搜索、調優(yōu)逐步逼近函數極值點或最優(yōu)點,直到滿足一定的精度時終止迭代計算,最后所逼近的設計點即最優(yōu)點，所得到的解即一定精度下的近似解,三、優(yōu)化設計的迭代計算,（三）迭代計算的迭代過程,由選定的初始點x(0)出發(fā),使?jié)M足,由于各設計點的函數值依次下降，可見迭代點不斷向理論最優(yōu)點逼近，最后可得到一定精度下的近似最優(yōu)點，記作x*。,三、優(yōu)化設計的迭代計算,迭代過程圖示：,X(0),X(1),X(k),s(0),a(0),X(2),x*,三、優(yōu)化設計的迭代計算,（四）迭代計算的終止準則,由于數值迭代是逐步逼近最優(yōu)點而獲得近似解的，所以要考慮優(yōu)化問題解的收斂性及迭代過程的終止條件。,1、迭代的收斂性,指某種迭代程序產生的序列xk (k=0,1,2,)收斂于,2、通常采用的終止準則,（1）點距準則,|xk+1-xk| 1,相鄰兩個設計點的距離已達到充分小,或,兩迭代點的坐標分量之差,三、優(yōu)化設計的迭代計算,（2）函數下降量準則,相鄰迭代點的函數值下降量達到充分小,或,相鄰迭代點的函數值的相對值達到充分小,（3）梯度準則,無約束最優(yōu)化,非線性規(guī)劃建模,教學目的,教學內容,2、掌握用數學軟件包求解無約束最優(yōu)化問題。,1、了解無約束最優(yōu)化基本算法。,1.無約束優(yōu)化基本思想及基本算法。,3. 用MATLAB求解無約束優(yōu)化問題。,2. MATLAB優(yōu)化工具箱簡介,無約束最優(yōu)化問題,求解無約束最優(yōu)化問題的的基本思想,*無約束最優(yōu)化問題的基本算法,標準形式：,求解無約束最優(yōu)化問題的基本思想,求解的基本思想 ( 以二元函數為例 ),5,3,1,連續(xù)可微,多局部極小,唯一極小 (全局極小),搜索過程,最優(yōu)點 (1 1) 初始點 (-1 1),-1,1,4.00,-0.79,0.58,3.39,-0.53,0.23,2.60,-0.18,0.00,1.50,0.09,-0.03,0.98,0.37,0.11,0.47,0.59,0.33,0.20,0.80,0.63,0.05,0.95,0.90,0.003,0.99,0.99,1E-4,0.999,0.998,1E-5,0.9997,0.9998,1E-8,坐標輪換法,一、坐標輪換法的迭代過程 如圖，以二次函數為例。,x2,x0 x01,x02 x21,x11,x12,x1,x21,坐標輪換法,任取一初始點x0作為第一輪的始點x01,先沿第一坐標軸的方向e1作一維搜索，用一維優(yōu)化方法確定最優(yōu)步長11，得第一輪的第一個迭代點x11=x01+ 11 e1，然后以x11為新起點，沿第二坐標軸的方向e2作一維搜索，確定步長21 ，得第一輪的第二個迭代點x21=x11+ 11 e2 第二輪迭代，需要 x11x01 x12 x02+ 12 e1 x22=x12+ 22 e2 依次類推，不斷迭代，目標函數值不斷下降，最后逼近該目標函數的最優(yōu)點。,坐標輪換法,二、終止準則 可以采用點距準則或者其它準則。 注意： 若采用點距準則或函數值準則，其中采用的點應該是一輪迭代的始點和終點，而不是某搜索方向的前后迭代點。,坐標輪換法,三、坐標輪換法的流程圖,入口,給定：x0，,K=1,i=1,Xik=x0,沿ei方向一維搜索求i xik=xi-1k+ ikei x=xk f=f(x),i=n?,|xnk-x0k|?,x*=x f*=f(x*),出口,i=i+1,x0=x0k,k=k+1,N,Y,N,Y,坐標輪換法,四、例題 五、小結 坐標輪換法程序簡單，易于掌握。但是計算效率比較低，尤其是當優(yōu)化問題的維數較高時更為嚴重。一般把此種方法應用于維數小于10的低維優(yōu)化問題。 對于目標函數存在 “脊線”的情況，在脊線 的尖點處沒有一個坐標方 向可以使函數值下降，只 有在銳角所包含的范圍搜 索才可以達到函數值下降 的目的，故坐標輪換法對 此類函數會失效。,x2,x1,脊線,鮑威爾方法方向加速法,鮑威爾方法是直接利用函數值來構造共軛方向的一種共軛方向法。這種方法是在研究具有正定矩陣G的二次函數 的極小化問題時形成的。其基本思想是在不用導數的前提下，在迭代中逐次構造G的共軛方向。 一、共軛方向的概念 設G為一正定對稱矩陣，若有一組非零向量S1，S2， ，Sn滿足 SiTGSj=0 (i j )，則稱這組向量關于矩陣G共軛。 共軛方向對于構造一種有效的算法是很重要的。以正定二元二次函數為例，我們進行探討。,鮑威爾方法,正定的二元二次函數的等值線為一組橢圓，任選初始點x0沿某個下降方向S0作一維搜索，得x1 x1=x0+ 0S0此時，點x1的梯度必然與方向S0垂直，即有 f(x1)TS0=0S0與某一等值線相切于x1點。下一次的迭代，若選擇負梯度方向為搜索方向，將產生鋸齒現象。為避免鋸齒的產生，我們取迭代方向S1直指極小點x*，如圖所示。,x0,x*,x1,1S1,- f(x1),S1,0S0,鮑威爾方法,若選定這樣的方向，則對于二元二次函數只需進行S0 和S1兩次搜索就可以求得極小點x*，即 x*=x1+ 1S1由于f(x1)=Gx1+b，當x1 x*時 1S1 0，由于x*是f （x）的極小點，應滿足極值必要條件，故f(x*)=Gx*+b=0，即 f(x*)=G（x1+ 1S1）+b= f(x1)+ 1GS1 =0 等式兩端同乘以(S0)T,得（S0）TGS1=0 ，故兩個向量S0、S1是G的共軛向量。 因此，若要使第二個迭代點成為正定二元二次函數的極小點，只要使兩次一維搜索的方向S0、S1關于函數的二階導數矩陣G共軛就可以了。,鮑威爾方法,二、共軛 方向的生成 設xk、xk+1為從不同點出發(fā)，沿同一方向Sj 進行一維搜索得到的兩個極小點，如圖所示。根據梯度和等值面的性質， Sj 和xk、xk+1兩點處的梯度gk、gk+1之間存在如下關系 (Sj)Tgk=0 (Sj)Tgk+1=0又因為xk、xk+1兩點處的梯度 可表示為 gk=Gxk+b gk+1=Gxk+1+b 兩式相減，得 gk+1-gk=G(xk+1-xK),鮑威爾方法,因此有 (Sj)T (gk+1-gk)= (Sj)T G(xk+1-xK)=0 若取方向Sk= xk+1-xK，則Sk和Sj對G共軛。這說明只要沿方向Sj 分別對函數作兩次一維搜索，得到兩個極小點，則這兩點的連線方向就是與Sj 共軛的方向。,鮑威爾方法,三、鮑威爾基本算法 如圖所示，以三維二次目標函數的無約束優(yōu)化問題為例。,x1,x3,x2,e1,e2,e3,e2,e3,e3,x01,x11,x21,x31,x1,x12,x22,x32,x13,x23,x33,x2,x3,S3,-S1,S2,S2,S1,鮑威爾方法,鮑威爾基本算法的步驟： 第一環(huán)基本方向組取單位坐標矢量系e1、 e2、 e3 、 en，沿這些方向依次作一維搜索，然后將始末兩點相連作為新生方向，再沿新生方向作一維搜索，完成第一環(huán)的迭代。以后每環(huán)的基本方向組是將上環(huán)的第一個方向淘汰，上環(huán)的新生方向補入本環(huán)的最后而構成。n維目標函數完成n環(huán)的迭代過程稱為一輪。從這一輪的終點出發(fā)沿新生方向搜索所得到的極小點，作為下一輪迭代的始點。這樣就形成了算法的循環(huán)。,鮑威爾方法,鮑威爾基本算法的缺陷： 可能在某一環(huán)迭代中出現基本方向組為線性相關的矢量系的情況。如第k環(huán)中，產生新的方向： Sk=xnk-x0k=1kS1k+ 2kS2k + + nkSnk 式中， S1k、S2k 、 、 Snk為第k環(huán)基本方向組矢量， 1k 、 2k 、 、 nk為個方向的最優(yōu)步長。 若在第k環(huán)的優(yōu)化搜索過程中出現1k =0，則方向Sk表示為S2k 、 S3k 、 、 Snk的線性組合，以后的各次搜索將在降維的空間進行，無法得到n維空間的函數極小值，計算將失敗。,鮑威爾方法,鮑威爾基本算法的退化,如圖所示為一個三維優(yōu)化問題的示例，設第一環(huán)中1 =0 ，則新生方向與e2 、e3共面，隨后的各環(huán)方向組中，各矢量必在該平面內，使搜索局限于二維空間，不能得到最優(yōu)解。,鮑威爾方法,四、鮑威爾修正算法 在某環(huán)已經取得的n+1各方向中，選取n個線性無關的并且共軛程度盡可能高的方向作為下一環(huán)的基本方向組。 鮑威爾修正算法的搜索方向的構造： 在第k環(huán)的搜索中， x0k 為初始點，搜索方向為S1k、S2k 、 、 Snk，產生的新方向為Sk ，此方向的極小點為xk。點 xn+1k=2xnk-x0k ， 為x0k對xnk的映射點。 計算x0k 、 x1k 、 、 xnk 、 xk、 x n+1k 各點的函數值，記作： F1=F(x0k) F2=F(xnk) F3=F(xn+1k) = F(xmk) -F(xm-1k) 是第k環(huán)方向組中，依次沿各方向搜索函數值下降最大值，即Smk方向函數下降最大。,鮑威爾方法,為了構造第k+1環(huán)基本方向組，采用如下判別式： 按照一下兩種情況處理： 1、上式中至少一個不等式成立，則第k+1環(huán)的基本方向仍用老方向組S1k、S2k 、 、 Snk。 k+1的環(huán)初始點取 x0k+1=xnk F2<F3 x0k+1=xn+1k F2F3 2、兩式均不成立，則淘汰函數值下降最大的方向，并用第k環(huán)的新生方向補入k+1環(huán)基本方向組的最后，即k+1環(huán)的方向組為S1k、S2k 、 、 Sm-1k、Sm+1k 、 Snk 、 Sn+1k 。k+1環(huán)的初始點取 x0k+1=xk xk是第k環(huán)沿Sn+1k方向搜索的極小點。,鮑威爾方法,x0k,x1k,x2k,x3k,xm-1k,xmk,Snk,xnk,Smk,S3k,S2k,Sk,xk,xn+1k,F1,F2,F3,鮑威爾方法,鮑威爾算法的終止條件： | xk-x0k | 五、鮑威爾算法的迭代步驟及流程圖 六、例題,梯度法,優(yōu)化設計是追求目標函數值最小，因此，自然可以設想從某點出發(fā)，其搜索方向取該點的負梯度方向，使函數值在該點附近下降最快。這種方法也稱為最速下降法。 一、基本原理梯度法的迭代公式為： x(k+1)=x(k)-(k)g(k)其中g(k)是函數F(x)在迭代點x(k)處的梯度f(xk) ， (k)一般采用一維搜索的最優(yōu)步長，即 f(x(k+1)=f(x(k)-(k)g(k) =min f(x(k)-(k)g(k)=min()據一元函數極值條件和多元復合函數求導公式，得 ()= -( f(x(k)-(k)g(k)T g(k) =0即 （ f(x(k+1)T g(k) =0或 (g(k+1)Tg(k)=0,梯度法,此式表明，相鄰的兩個迭代點的梯度是彼此正交的。也即在梯度的迭代過程中，相鄰的搜索方向相互垂直。梯度法向極小點的逼近路徑是鋸齒形路線，越接近極小點，鋸齒越細，前進速度越慢。 這是因為，梯度是函數的局部性質，從局部上看，在該點附近函數的下降最快，但從總體上看則走了許多彎路，因此函數值的下降并不快。,梯度法,二、迭代終止條件 采用梯度準則： | g(k) | 三、迭代步驟（1）任選初始迭代點x(0)，選收斂精度 。（2）確定x(k)點的梯度（開始k=0)（3）判斷是否滿足終止條件| g(k) | ？若滿足輸出最優(yōu)解，結束計算。否則轉下步。（4）從x(k)點出發(fā)，沿-g(k)方向作一維搜索求最優(yōu)步長(k)。得下一迭代點 x(k+1)=x(k)-(k)g(k) ，令k=k+1 返回步驟（2）。,梯度法,四、梯度法流程圖,入口,給定： x(0)，,k=0,|g(k)| ？,x*=x(k) f*=f(x(k),出口,x(k)= x(0,計算：g(k),k=k+1,沿g(k)方向一維搜索， 求最優(yōu)步長(k)。,N,Y,共軛梯度法,共軛梯度法是共軛方向法的一種，因為該方法中每一個共軛向量都是依賴于迭代點處的負梯度而構造出來的，所以稱作共軛梯度法。 一、共軛梯度法的搜索方向 共軛梯度法的搜索方向采用梯度法基礎上的共軛方向，如圖所示，目標函數F(x)在迭代點xk+1處的負梯度為- f(xk+1)，該方向與前一搜索方向Sk互為正交，在此基礎上構造一種具有較高收斂速度的算法，該算法的搜索方向要滿足以下兩個條件： （1）以xk+1點出發(fā)的搜索方向 Sk+1是- f(xk+1)與Sk的線性組合。即,xk,x*,xk+1,- f(xk+1),Sk+1,Sk,共軛梯度法,Sk+1 = - f(xk+1) + kSk （2）以Sk與Sk+1為基底的子空間中，矢量Sk與Sk+1相共軛，即滿足 Sk+1T G Sk = 0 二、 k的確定 確定方法不作要求。記住 三、 共軛梯度法的算法 （1）選初始點x0和收斂精度。 （2）令k=0，計算S0 = - f(x0) 。 （3）沿Sk方向進行一維搜索求(k)，得 x(k+1)=x(k)+(k)S(k) （4）計算 f(xk+1) ，若| f(xk+1)| ，則終止迭代，取x*=xk+1；否則進行下一步。,共軛梯度法,（5）檢查搜索次數，若k=n，則令x0=xk+1，轉（2）， 否則，進行下一步。 （6）構造新的共軛方向 Sk+1 = - f(xk+1) + kSk 令k=k+1，轉（3）,共軛梯度法,四、共軛梯度法流程圖,入口,k=0，計算：- f(x0),| f(xk+1) | ？,出口,求(k) ，x(k+1)= x(k) +(k)S(k),計算： f(xk+1),x*=xk+1 f(x*)=f(xk+1),Y,N,給定： x(0)，,k n ？,x0=xk+1,N,Y,Sk+1 = - f(xk+1) + kSk,K=K+1,共軛梯度法,五、共軛梯度法的特點 共軛梯度法屬于解析法，其算法需求一階導數，所用公式及算法簡單，所需存儲量少。該方法以正定二次函數的共軛方向理論為基礎，對二次型函數可以經過有限步達到極小點，所以具有二次收斂性。但是對于非二次型函數，以及在實際計算中由于計算機舍入誤差的影響，雖然經過n次迭代，仍不能達到極小點，則通常以重置負梯度方向開始，搜索直至達到預定精度，其收斂速度也是較快的。 六、例題,牛頓法,牛頓法是求無約束最優(yōu)解的一種古典解析算法。牛頓法可以分為原始牛頓法和阻尼牛頓法兩種。實際中應用較多的是阻尼牛頓法。,原始牛頓法,一、原始牛頓法的基本思想 在第k次迭代的迭代點x(k)鄰域內，用一個二次函數去近似代替原目標函數f(x)，然后求出該二次函數的極小點作為對原目標函數求優(yōu)的下一個迭代點，依次類推，通過多次重復迭代，是迭代點逐步逼近原目標函數的極小點。如圖所示。,F(x),1(x),0(x),x0,x1,x2,x*,原始牛頓法,二、原始牛頓法的迭代公式 設目標函數f(x)具有連續(xù)的一、二階導數，在x(k)點鄰域內取f(x)的二次泰勒多項式作近似式，即取逼近函數(x)為設xk+1為(x)極小點，根據極值的必要條件，應有(xk+1)=0，即 f(xk)+ G(xk)x=0 又 x= xk+1 - xk 得 xk+1 = xk- G(xk)-1 f(xk) 即牛頓法迭代公式，方向- G(xk)-1 f(xk)稱為牛頓方向,原始牛頓法,三、原始牛頓法的特點 若用原始牛頓法求某二次目標函數的最優(yōu)解，則構造的逼近函數與原目標函數是完全相同的二次式，其等值線完全重合，故從任一點出發(fā)，一定可以一次達到目標函數的極小點。 因此，牛頓法是具有二次收斂性的算法。其優(yōu)點是：對于二次正定函數，迭代一次即可以得到最優(yōu)解，對于非二次函數，若函數二次性較強或迭代點已經進入最優(yōu)點的較小鄰域，則收斂速度也很快。 原始牛頓法的缺點是：由于迭代點的位置是按照極值條件確定的，并未沿函數值下降方向搜索，因此，對于非二次函數，有時會使函數值上升，即 f(xk+1) f(xk)，而使計算失敗。,阻尼牛頓法,一、對原始牛頓法的改進 為解決原始牛頓法的不足，加入搜索步長(k)因此，迭代公式變?yōu)椋?xk+1 = xk - (k)G(xk)-1 f(xk) 這就是阻尼牛頓法的迭代公式，最優(yōu)步長(k)也稱為阻尼因子，是沿牛頓方向一維搜索得到的最優(yōu)步長。 二、阻尼牛頓法的迭代步驟 （1）給定初始點和收斂精度 （2）計算 f(xk) 、 G(xk)、 G(xk)-1 （3）求xk+1 = xk - (k)G(xk)-1 f(xk) （4）檢查收斂精度，若| xk+1- xk < |則x*=xk+1，停止，否則 k=k+1，返回（2）繼續(xù),阻尼牛頓法,三、阻尼牛頓法的特點 優(yōu)點： 由于阻尼牛頓法每次迭代都在牛頓方向進行一維搜索，避免了迭代后函數值上升的現象，從而保持了牛頓法的二次收斂性，而對初始點的選擇沒有苛刻的要求。 缺點： 1、對目標函數要求苛刻，要求函數具有連續(xù)的一、二階導數；為保證函數的穩(wěn)定下降，海賽矩陣必須正定；為求逆陣要求海賽矩陣非奇異。 2、計算復雜且計算量大，存儲量大,變尺度法,Matlab優(yōu)化工具箱簡介,1.MATLAB求解優(yōu)化問題的主要函數,2. 優(yōu)化函數的輸入變量,使用優(yōu)化函數或優(yōu)化工具箱中其它優(yōu)化函數時, 輸入變量見下表:,3. 優(yōu)化函數的輸出變量下表:,4控制參數options的設置,(3) MaxIter: 允許進行迭代的最大次數,取值為正整數.,Options中常用的幾個參數的名稱、含義、取值如下:,(1)Display: 顯示水平.取值為off時,不顯示輸出; 取值為iter時,顯示每次迭代的信息;取值為final時,顯示最終結果.默認值為final.,(2)MaxFunEvals: 允許進行函數評價的最大次數,取值為正整數.,例：opts=optimset(Display,iter,TolFun,1e-8) 該語句創(chuàng)建一個稱為opts的優(yōu)化選項結構,其中顯示參數設為iter, TolFun參數設為1e-8.,控制參數options可以通過函數optimset創(chuàng)建或修改。命令的格式如下：,(1) options=optimset(optimfun) 創(chuàng)建一個含有所有參數名,并與優(yōu)化函數optimfun相關的默認值的選項結構options.,（2）options=optimset(param1,value1,param2,value2,.) 創(chuàng)建一個名稱為options的優(yōu)化選項參數,其中指定的參數具有指定值,所有未指定的參數取默認值.,(3)options=optimset(oldops,param1,value1,param2, value2,.) 創(chuàng)建名稱為oldops的參數的拷貝,用指定的參數值修改oldops中相應的參數.,用Matlab解無約束優(yōu)化問題,其中（3）、（4）、（5）的等式右邊可選用（1）或（2）的等式右邊。 函數fminbnd的算法基于黃金分割法和二次插值法，它要求目標函數必須是連續(xù)函數，并可能只給出局部最優(yōu)解。,常用格式如下： （1）x= fminbnd (fun,x1,x2) （2）x= fminbnd (fun,x1,x2 ，options) （3）x，fval= fminbnd（.） （4）x，fval，exitflag= fminbnd（.） （5）x，fval，exitflag，output= fminbnd（.）,主程序為wliti1.m: f=2*exp(-x).*sin(x); fplot(f,0,8); %作圖語句 xmin,ymin=fminbnd (f, 0,8) f1=-2*exp(-x).*sin(x); xmax,ymax=fminbnd (f1, 0,8),例2 對邊長為3米的正方形鐵板，在四個角剪去相等的正方形以制成方形無蓋水槽，問如何剪法使水槽的容積最大？,解,先編寫M文件fun0.m如下: function f=fun0(x) f=-(3-2*x).2*x;,主程序為wliti2.m: x,fval=fminbnd(fun0,0,1.5); xmax=x fmax=-fval,運算結果為: xmax = 0.5000,fmax =2.0000.即剪掉的正方形的邊長為0.5米時水槽的容積最大,最大容積為2立方米.,命令格式為: （1）x= fminunc（fun,X0 ）；或x=fminsearch（fun,X0 ） （2）x= fminunc（fun,X0 ，options）； 或x=fminsearch（fun,X0 ，options） （3）x，fval= fminunc（.）； 或x，fval= fminsearch（.） （4）x，fval，exitflag= fminunc（.）； 或x，fval，exitflag= fminsearch （5）x，fval，exitflag，output= fminunc（.）； 或x，fval，exitflag，output= fminsearch（.）,2、多元函數無約束優(yōu)化問題,標準型為：min F(X),3 fminunc為中型優(yōu)化算法的步長一維搜索提供了兩種算法， 由options中參數LineSearchType控制： LineSearchType=quadcubic(缺省值)，混合的二次和三 次多項式插值； LineSearchType=cubicpoly，三次多項式插,使用fminunc和 fminsearch可能會得到局部最優(yōu)解.,說明:,fminsearch是用單純形法尋優(yōu). fminunc的算法見以下幾點說明：,1 fminunc為無約束優(yōu)化提供了大型優(yōu)化和中型優(yōu)化算法。由options中的參數LargeScale控制： LargeScale=on(默認值),使用大型算法 LargeScale=off(默認值),使用中型算法,2 fminunc為中型優(yōu)化算法的搜索方向提供了4種算法，由 options中的參數HessUpdate控制： HessUpdate=bfgs（默認值），擬牛頓法的BFGS公式； HessUpdate=dfp，擬牛頓法的DFP公式； HessUpdate=steepdesc，最速下降法,例3 min f(x)=(4x12+2x22+4x1x2+2x2+1)*exp(x1),1、編寫M-文件 fun1.m: function f = fun1 (x) f = exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1); 2、輸入M文件wliti3.m如下: x0 = -1, 1; x=fminunc(fun1,x0); y=fun1(x),3、運行結果: x= 0.5000 -1.0000 y = 1.3029e-10,3.用fminsearch函數求解,輸入命令: f=100*(x(2)-x(1)2)2+(1-x(1)2; x,fval,exitflag,output=fminsearch(f, -1.2 2),運行結果: x =1.0000 1.0000 fval =1.9151e-010 exitflag = 1 output = iterations: 108 funcCount: 202 algorithm: Nelder-Mead simplex direct search,4. 用fminunc 函數,(1)建立M-文件fun2.m function f=fun2(x) f=100*(x(2)-x(1)2)2+(1-x(1)2,(2)主程序wliti44.m,Rosenbrock函數不同算法的計算結果,可以看出，最速下降法的結果最差.因為最速下降法特別不適合于從一狹長通道到達最優(yōu)解的情況.,例5 產銷量的最佳安排 某廠生產一種產品有甲、乙兩個牌號，討論在產銷平衡的情況下如何確定各自的產量，使總利潤最大. 所謂產銷平衡指工廠的產量等于市場上的銷量.,基本假設,1價格與銷量成線性關系,2成本與產量成負指數關系,模型建立,若根據大量的統(tǒng)計數據,求出系數b1=100,a11=1,a12=0.1,b2=280, a21=0.2,a22=2,r1=30,1=0.015,c1=20, r2=100,2=0.02,c2=30,則 問題轉化為無約束優(yōu)化問題：求甲,乙兩個牌號的產量x1，x2，使 總利潤z最大.,為簡化模型,先忽略成本,并令a12=0,a21=0,問題轉化為求: z1 = ( b1 - a11x1 ) x1 + ( b2 - a22x2 ) x2 的極值. 顯然其解為x1 = b1/2a11 = 50, x2 = b2/2a22 = 70, 我們把它作為原問題的初始值.,總利潤為： z(x1,x2)=(p1-q1)x1+(p2-q2)x2,模型求解,1.建立M-文件fun.m: function f = fun(x) y1=(100-x(1)- 0.1*x(2)-(30*exp(-0.015*x(1)+20)*x(1); y2=(280-0.2*x(1)- 2*x(2)-(100*exp(-0.02*x(2)+30)*x(2); f=-y1-y2;,2.輸入命令: x0=50,70; x=fminunc(fun,x0), z=fun(x),3.計算結果: x=23.9025, 62.4977, z=6.4135e+003 即甲的產量為23.9025,乙的產量為62.4977,最大利潤為6413.5.,非線性規(guī)劃建模,有約束非線性規(guī)劃,教學目的,教學內容,2、掌握用數學軟件求解優(yōu)化問題。,1、直觀了解非線性規(guī)劃的基本內容。,1.非線性規(guī)劃的基本理論。,2. 用數學軟件求解非線性規(guī)劃。,3. 鋼管訂購及運輸優(yōu)化模型,*非線性規(guī)劃的基本解法,非線性規(guī)劃,非線性規(guī)劃的基本概念,定義 如果目標函數或約束條件中至少有一個是非線性函數時的最優(yōu)化問題就叫做非線性規(guī)劃問題,非線性規(guī)劃的基本概念,一般形式: （1） 其中 ， 是定義在 En 上的實值函數，簡記:,其它情況: 求目標函數的最大值或約束條件為小于等于零的情況，都可通過取其相反數化為上述一般形式,定義1 把滿足問題（1）中條件的解 稱為可行解（或可行點），所有可行點的集合稱為可行集（或可行域）記為D即 問題(1)可簡記為 ,定義2 對于問題(1)，設 ， 若存在 ,使得對一切 ，且 ，都有 ，則稱X*是f(X)在D上的局部極小值點（局部最優(yōu)解）特別地當 時,若 ，則稱X*是f(X)在D上的嚴格局部極小值點（嚴格局部最優(yōu)解）,定義3 對于問題(1),設 ,對任意的 ,都有 則稱X*是f(X)在D上的全局極小值點（全局最優(yōu)解）特別地當 時，若 ，則稱X*是f(X)在D上的嚴格全局極小值點（嚴格全局最優(yōu)解）,非線性規(guī)劃的基本解法,SUTM外點法,SUTM內點法（障礙罰函數法）,1. 罰函數法,2. 近似規(guī)劃法,罰函數法,罰函數法基本思想是通過構造罰函數把約束問題轉化為一系列無約束最優(yōu)化問題，進而用無約束最優(yōu)化方法去求解這類方法稱為序列無約束最小化方法.簡稱為SUMT法 其一為SUMT外點法,其二為SUMT內點法,其中T(X,M)稱為罰函數，M稱為罰因子，帶M的項稱為罰項，這里的罰函數只對不滿足約束條件的點實行懲罰:當 時，滿足各 ，故罰項=0，不受懲罰當 時，必有 的約束條件，故罰項0，要受懲罰,SUTM外點法,罰函數法的缺點是：每個近似最優(yōu)解Xk往往不是容許解，而只能近似滿足約束，在實際問題中這種結果可能不能使用；在解一系列無約束問題中，計算量太大，特別是隨著Mk的增大，可能導致錯誤,1、任意給定初始點X0，取M11，給定允許誤差 ，令k=1； 2、求無約束極值問題 的最優(yōu)解，設為Xk=X(Mk)，即 ； 3、若存在 ，使 ，則取MkM( )令k=k+1返回（2），否則，停止迭代得最優(yōu)解 . 計算時也可將收斂性判別準則 改為 .,SUTM外點法(罰函數法)的迭代步驟,SUTM內點法（障礙函數法）,內點法的迭代步驟,近似規(guī)劃法的基本思想：將問題(3)中的目標函數 和約束條件 近似為線性函數，并對變量的取值范圍加以限制，從而得到一個近似線性規(guī)劃問題,再用單純形法求解之，把其符合原始條件的最優(yōu)解作為(3)的解的近似,近似規(guī)劃法,每得到一個近似解后，都從這點出發(fā)，重復以上步驟,這樣，通過求解一系列線性規(guī)劃問題，產生一個由線性規(guī)劃最優(yōu)解組成的序列，經驗表明，這樣的序列往往收斂于非線性規(guī)劃問題的解。,近似規(guī)劃法的算法步驟如下,用MATLAB軟件求解,其輸入格式如下: 1.x=quadprog(H,C,A,b); 2.x=quadprog(H,C,A,b,Aeq,beq); 3.x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB); 4.x=quadprog(H,C,A,b, Aeq,beq ,VLB,VUB,X0); 5.x=quadprog(H,C,A,b, Aeq,beq ,VLB,VUB,X0,options); 6.x,fval=quaprog(.); 7.x,fval,exitflag=quaprog(.); 8.x,fval,exitflag,output=quaprog(.);,1、二次規(guī)劃,例1 min f(x1,x2)=-2x1-6x2+x12-2x1x2+2x22 s.t. x1+x22 -x1+2x22 x10, x20,1、寫成標準形式：,2、 輸入命令： H=1 -1; -1 2; c=-2 ;-6;A=1 1; -1 2;b=2;2; Aeq=;beq=; VLB=0;0;VUB=; x,z=quadprog(H,c,A,b,Aeq,beq,VLB,VUB),3、運算結果為： x =0.6667 1.3333 z = -8.2222,s.t.,1. 首先建立M文件fun.m,定義目標函數F（X）: function f=fun(X); f=F(X);,2、一般非線性規(guī)劃,其中X為n維變元向量，G(X)與Ceq(X)均為非線性函數組成的向量，其它變量的含義與線性規(guī)劃、二次規(guī)劃中相同.用Matlab求解上述問題，基本步驟分三步：,3. 建立主程序.非線性規(guī)劃求解的函數是fmincon,命令的基本格式如下： (1) x=fmincon(fun,X0,A,b) (2) x=fmincon(fun,X0,A,b,Aeq,beq) (3) x=fmincon(fun,X0,A,b, Aeq,beq,VLB,VUB) (4) x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon) (5)x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon,options) (6) x,fval= fmincon(.) (7) x,fval,exitflag= fmincon(.) (8)x,fval,exitflag,output= fmincon(.),輸出極值點,M文件,迭代的初值,參數說明,變量上下限,注意： 1 fmincon函數提供了大型優(yōu)化算法和中型優(yōu)化算法。默認時，若在fun函數中提供了梯度（options參數的GradObj設置為on），并且只有上下界存在或只有等式約束，fmincon函數將選擇大型算法。當既有等式約束又有梯度約束時，使用中型算法。 2 fmincon函數的中型算法使用的是序列二次規(guī)劃法。在每一步迭代中求解二次規(guī)劃子問題，并用BFGS法更新拉格朗日Hessian矩陣。 3 fmincon函數可能會給出局部最優(yōu)解，這與初值X0的選取有關。,1、寫成標準形式： s.t.,2x1+3x2 6 s.t x1+4x2 5 x1,x2 0,例2,2、先建立M-文件 fun3.m: function f=fun3(x); f=-x(1)-2*x(2)+(1/2)*x(1)2+(1/2)*x(2)2,3、再建立主程序youh2.m： x0=1;1; A=2 3 ;1 4; b=6;5; Aeq=;beq=; VLB=0;0; VUB=; x,fval=fmincon(fun3,x0,A,b,Aeq,beq,VLB,VUB),4、運算結果為： x = 0.7647 1.0588 fval = -2.0294,1先建立M文件 fun4.m,定義目標函數: function f=fun4(x); f=exp(x(1) *(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1);,x1+x2=0 s.t. 1.5+x1x2 - x1 - x2 0 -x1x2 10 0,例3,2再建立M文件mycon.m定義非線性約束： function g,ceq=mycon(x) g=x(1)+x(2);1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10;,3主程序youh3.m為: x0=-1;1; A=;b=; Aeq=1 1;beq=0; vlb=;vub=; x,fval=fmincon(fun4,x0,A,b,Aeq,beq,vlb,vub,mycon),3. 運算結果為： x = -1.2250 1.2250 fval = 1.8951,例4,1先建立M-文件fun.m定義目標函數: function f=fun(x); f=-2*x(1)-x(2);,2再建立M文件mycon2.m定義非線性約束： function g,ceq=mycon2(x) g=x(1)2+x(2)2-25;x(1)2-x(2)2-7;,3. 主程序fxx.m為: x0=3;2.5; VLB=0 0;VUB=5 10; x,fval,exitflag,output =fmincon(fun,x0,VLB,VUB,mycon2),4. 運算結果為: x = 4.0000 3.0000 fval =-11.0000 exitflag = 1 output = iterations: 4 funcCount: 17 stepsize: 1 algorithm: 1x44 char firstorderopt: cgiterations: ,應用實例： 供應與選址,某公司有6個建筑工地要開工，每個工地的位置（用平面坐標系a，b表示，距離單位：千米 ）及水泥日用量d(噸)由下表給出。目前有兩個臨時料場位于A(5,1)，B(2,7)，日儲量各有20噸。假設從料場到工地之間均有直線道路相連。 （1）試制定每天的供應計劃，即從A，B兩料場分別向各工地運送多少噸水泥，使總的噸千米數最小。 （2）為了進一步減少噸千米數，打算舍棄兩個臨時料場，改建兩個新的，日儲量各為20噸，問應建在何處，節(jié)省的噸千米數有多大？,（一）、建立模型,記工地的位置為(ai，bi)，水泥日用量為di，i=1,6;料場位置為(xj，yj)，日儲量為ej，j=1,2；從料場j向工地i的運送量為Xij。,當用臨時料場時決策變量為：Xij， 當不用臨時料場時決策變量為：Xij，xj，yj。,（二）使用臨時料場的情形,使用兩個臨時料場A(5,1)，B(2,7).求從料場j向工地i的運送量為Xij，在各工地用量必須滿足和各料場運送量不超過日儲量的條件下，使總的噸千米數最小，這是線性規(guī)劃問題. 線性規(guī)劃模型為：,設X11=X1, X21= X 2, X31= X 3, X41= X 4, X51= X 5, X61= X 6 X12= X 7, X22= X 8, X32= X 9, X42= X 10, X52= X 11, X62= X 12 編寫程序gying1.m,計算結果為：,x = 3.0000 5.0000 0.0000 7.0000 0.0000 1.0000 0.0000 0.0000 4.0000 0.0000 6.0000 10.0000 fval = 136.2275,（三）改建兩個新料場的情形,改建兩個新料場,要同時確定料場的位置(xj,yj)和運送量Xij，在同樣條件下使總噸千米數最小.這是非線性規(guī)劃問題。非線性規(guī)劃模型為：,設 X11=X1, X21= X 2, X31= X 3, X41= X 4, X51= X 5, X61= X 6 X12= X 7, X22= X 8, X32= X 9, X42= X 10, X52= X 11, X62= X 12 x1=X13, y1=X14, x2=X15, y2=X16,（1）先編寫M文件liaoch.m定義目標函數。,(2) 取初值為線性規(guī)劃的計算結果及臨時料場的坐標: x0=3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7; 編寫主程序gying2.m.,(3) 計算結果為：,x= 3.0000 5.0000 0.0707 7.0000 0 0.9293 0 0 3.9293 0 6.0000 10.0707 6.3875 4.3943 5.7511 7.1867 fval = 105.4626 exitflag = 1,(4) 若修改主程序gying2.m, 取初值為上面的計算結果: x0= 3.0000 5.0000 0.0707 7.0000 0 0.9293 0 0 3.9293 0 6.0000 10.0707 6.3875 4.3943 5.7511 7.1867,得結果為: x=3.0000 5.0000 0.3094 7.0000 0.0108 0.6798 0 0 3.6906 0 5.9892 10.3202 5.5369 4.9194 5.8291 7.2852 fval =103.4760 exitflag = 1,總的噸千米數比上面結果略優(yōu).,(5) 若再取剛得出的結果為初值, 卻計算不出最優(yōu)解.,(6) 若取初值為: x0=3 5 4 7 1 0 0 0 0 0 5 11 5.6348 4.8687 7.2479 7.7499, 則計算結果為: x=3.0000 5.0000 4.0000 7.0000 1.0000 0 0 0 0 0 5.0000 11.0000 5.6959 4.9285 7.2500 7.7500 fval =89.8835 exitflag = 1 總的噸千米數89.8835比上面結果更好.,通過此例可看出fmincon函數在選取初值上的重要性.,連續(xù)結構體建模與優(yōu)化設計,連續(xù)體形狀優(yōu)化例一,連續(xù)體形狀優(yōu)化例二,注：并非結 構變形動畫， 而是結構受 力狀況下自 動尋優(yōu)過程 注意：網格 變化過程 （網格拓撲 不改變，類 似蜘蛛網受 力變化）,離散模型敏度分析方法,對E求導,？,對E求導,離散法主要為有限 元程序開發(fā)者應用,代入已知得所求,連續(xù)模型敏度分析方法,敏度表達為位移、 應力等和區(qū)域形狀 變化量的邊界積分 或區(qū)域積分的形式,圖 連續(xù)體變形過程的映射示意圖,網格均勻化(Mesh Smoothing),啟發(fā)式算法,啟發(fā)式算法(heuristic algorithm)是相對于最優(yōu)化算法提出的,啟發(fā)式算法定義：一種技術。這種技術使得在可接受的計算 的計算費用內去尋找最好的解，但不一定能保證所得解的可 行性和最優(yōu)性，甚至在多數情況下，無法闡述所得解同最優(yōu) 解的近似程度,結構優(yōu)化之滿應力法,滿應力法理論上受啟發(fā)于生物進化方式,滿應力法是一種經過實踐檢驗的較好的優(yōu)化方法,滿應力法數學模型可以看成是最優(yōu)問題的近似模型,連續(xù)體拓撲優(yōu)化例一,注意：圖片中的孤立塊和鋸齒邊界形狀,連續(xù)體拓撲優(yōu)化例二,注意：圖片中的棋盤圖案,連續(xù)體拓撲優(yōu)化例三,注意：拓撲優(yōu)化是 一種概念設計,映射反演方法例一,微分方程,代數方程,求解代數方程,方程解,拉氏變換,反拉氏變換,RMI方法使用條件：,映射須是兩類數學對象之間 的一一對應關系； 映射須是可定映的； 逆映射（反演）具有能行性,授課結束 謝 謝 大 家！,