多目標規(guī)劃MATLAB2012年wgx

1、MATLAB求解多目標規(guī)劃一、0-1規(guī)劃的MATLAB求解數(shù)學模型：MIN fx S.T. Ax=b Aeqx=beq x=0,1命令格式：x = bintprog(f)x = bintprog(f, A, b)x = bintprog(f, A, b, Aeq, beq)x = bintprog(f, A, b, Aeq, beq, x0)x = bintprog(f, A, b, Aeq, beq, x0, options)x, fval = bintprog(.)x,fval, exitflag = bintprog(.)x, fval, exitflag, output = bintp

2、rog(.)數(shù)學模型：MIN lambda S.T. F(x)-weight* lambda =goal(達到目標） Ax=b（線性不等式約束） Aeqx=beq（線性等式約束） C(x)=0（非線性不等式約束） Ceq(x)=0（非線性等式約束） lb=x 1 % Two output arguments G = . % Gradients evaluated at xEndThe gradient consists of the partial derivative dF/dx of each F at the point x.二、多目標規(guī)劃的MATLAB求解二、多目標規(guī)劃的MATLAB求

3、解有關(guān)優(yōu)化參數(shù)設置：options = optimset(GradConstr,on）約束條件的梯度方向參數(shù)設置為on時，用下列函數(shù)定義：function c,ceq,GC,GCeq = mycon(x)c = . % Nonlinear inequalities at xceq = . % Nonlinear equalities at xif nargout 2 % Nonlcon called with 4 outputs GC = . % Gradients of the inequalities GCeq = . % Gradients of the equalitiesEnd注意：

4、一般 weight=abs(goal)模型：x=(A+BKC)x+Bu，設計K滿足目標： Y=Cx1）循環(huán)系統(tǒng)的特征值（由命令eig(A+B*K*C)確定）的目標為goal = -5,-3,-12）K中元素均在-4,4中; 設特征值的weight= abs(goal),定義目標函數(shù)F如下： function F = eigfun(K,A,B,C)F = sort(eig(A+B*K*C); % Evaluate objectives，由小到大排列優(yōu)化程序為：A = -0.5 0 0; 0 -2 10; 0 1 -2;B = 1 0; -2 2; 0 1;C = 1 0 0; 0 0 1; K0

5、 = -1 -1; -1 -1; % Initialize controller matrixgoal = -5 -3 -1; % Set goal values for the eigenvaluesweight = abs(goal) % Set weight for same percentagelb = -4*ones(size(K0); % Set lower bounds on the controllerub = 4*ones(size(K0); % Set upper bounds on the controlleroptions = optimset(Display,iter

6、); % Set display parameterK,fval,attainfactor = fgoalattain(K)eigfun(K,A,B,C). goal,weight,lb,ub,options)二、舉例-有關(guān)循環(huán)控制系統(tǒng)優(yōu)化問題運行結(jié)果如下Active constraints: 1 2 4 9 10K = -4.0000 -0.2564 -4.0000 -4.0000fval = -6.9313 -4.1588 -1.4099attainfactor = -0.3863二、舉例-有關(guān)循環(huán)控制系統(tǒng)優(yōu)化問題如果至少保證38.63%的目標精確匹配，設置GoalsExactAchiev

7、e參數(shù)值為3options = optimset(GoalsExactAchieve,3);K,fval,attainfactor = fgoalattain(. (K)eigfun(K,A,B,C),K0,goal,weight,lb,ub,.options)After about seven iterations, a solution is K = -1.5954 1.2040 -0.4201 -2.9046fval = -5.0000 -3.0000 -1.0000attainfactor = 1.0859e-20表明目標已完全匹配二、舉例-有關(guān)循環(huán)控制系統(tǒng)優(yōu)化問題謝謝！ 初等模型舉例

8、 常見類型定性模型經(jīng)驗公式（擬合、插值）量綱分析比例模型 2.1 崖高的估算假如你站在崖頂且身上帶著一只具有跑表功 能的計算器，你也許會出于好奇心想用扔下 一塊石頭聽回聲的方法來估計山崖的高度， 假定你能準確地測定時間，你又怎樣來推算 山崖的高度呢，請你分析一下這一問題。我有一只具有跑 表功能的計算器。方法一假定空氣阻力不計，可以直接利用自由落體運動的公式來計算。例如， 設t=4秒，g=9.81米/秒2，則可求得h78.5米。 我學過微積分，我可以做 得更好，呵呵。 除去地球吸引力外，對石塊下落影響最大的當 屬空氣阻力。根據(jù)流體力學知識，此時可設空氣阻力正比于石塊下落的速度，阻力系 數(shù)K為常數(shù)

9、，因而，由牛頓第二定律可得： 令k=K/m，解得 代入初始條件 v(0)=0，得c=g/k，故有 再積分一次，得： 若設k=0.05并仍設 t=4秒，則可求 得h73.6米。 聽到回聲再按跑表，計算得到的時間中包含了 反應時間 進一步深入考慮不妨設平均反應時間 為0.1秒 ，假如仍 設t=4秒，扣除反應時間后應 為3.9秒，代入 式，求得h69.9米。 多測幾次，取平均值再一步深入考慮代入初始條 件h(0)=0，得到計算山崖高度的公式： 將e-kt用泰勒公式展開并 令k 0+ ，即可得出前面不考慮空氣阻力時的結(jié)果。還應考慮回聲傳回來所需要的時間。為此，令石塊下落 的真正時間 為t1，聲音傳回來

10、的時間記 為t2，還得解一個方程組： 這一方程組是非線性的，求解不太容易，為了估算崖高竟要去解一個非線性主程組似乎不合情理 相對于石塊速度，聲音速度要快得多，我們可 用方法二先求一次 h，令t2=h/340，校正t，求石塊下落時間 t1t-t2將t1代入式再算一次，得出崖高的近似值。例如， 若h=69.9米，則 t20.21秒，故 t13.69秒，求得 h62.3米。 2.2 錄像帶還能錄多長時間錄像機上有一個四位計數(shù)器，一盤 180分鐘的錄像帶在開始計數(shù)時為 0000，到結(jié)束時計數(shù)為1849，實際走時為185分20秒。我們從0084觀察到0147共用時間3分21秒。若錄像機目前的計數(shù)為142

11、8，問是否還能錄下一個 60分鐘的節(jié)目？rRl由得到又 因和 得 積分得到即從而有我們希望建立一個錄像帶已錄像時 間t與計數(shù)器計 數(shù)n之間的函數(shù)關(guān)系。為建立一個正確的模型，首 先必須搞清哪些量是常量，哪些量是變量。首先，錄像 帶的厚 度W是常量，它被繞在一個半徑 為r的園盤上，見圖。磁帶轉(zhuǎn)動中線速 度v顯然也是常數(shù)，否則圖象聲音必然會失真。此外，計數(shù)器的讀 數(shù)n與轉(zhuǎn)過的圈數(shù)有關(guān)，從而與轉(zhuǎn)過的角 度成正比。 rRl 此式中的三個參數(shù)W、v和r均不易精確測得，雖然我們可以從上式解出t與n的函數(shù)關(guān)系，但效果不佳，故令 則可將上式簡化為： 故令上式又可化簡記成 t= an2+bn t= an2+bn

12、rRl上式以a、b為參數(shù)顯然是一個十分明智的做法，它為公式的最終確立即參數(shù)求解提供了方便。將已知條件代入，得方程組： 從后兩式中消 去t1，解得a=0.0000291, b=0.04646,故t=0.0000291 n2+0.04646n，令n=1428，得到t=125.69（分）由于一盒錄像帶實際可錄像時間為185.33分，故尚可錄像時間 為59.64分，已不能再錄下一個60分鐘的節(jié)目了。 在解決實際問題時，注意觀察和善于想象是十分重要的，觀察與想象不僅能發(fā)現(xiàn)問題隱含的某些屬性，有時還能順理成章地找到解決實際問題的鑰匙。本節(jié)的幾個例子說明，猜測也是一種想象力。沒有合理而又大膽的猜測，很難做出

13、具有創(chuàng)新性的結(jié)果。開普勒的三大定律（尤其是后兩條）并非一眼就能看出的，它們隱含在行星運動的軌跡之中，隱含在第谷記錄下來的一大堆數(shù)據(jù)之中。歷史上這樣的例子實在太多了。在獲得了一定數(shù)量的資料數(shù)據(jù)后，人們常常會先去猜測某些結(jié)果，然后試圖去證明它。猜測一經(jīng)證明就成了定理，而定理一旦插上想象的翅膀，又常常會被推廣出許多更為廣泛的結(jié)果。即使猜測被證明是錯誤的，結(jié)果也決不是一無所獲的失敗而常常是對問題的更為深入的了解。 2.3最短路徑與最速方案問題 例5（最短路徑問題） 設有一個半徑為 r 的圓形湖，圓心為 O。A、B 位于湖的兩側(cè)，AB連線過O，見圖?，F(xiàn)擬從A點步行到B點，在不得進入湖中的限 制下，問怎樣

14、的路徑最近。 ABOr將湖想象成凸出地面的木樁， 在AB間拉一根軟線，當線被拉緊時將得到最短路徑。根據(jù)這樣的想象，猜測 可以如下得到最短路徑： 過A作圓的切線切圓于E，過B作圓的切線切圓 于F。最短路徑為由線 段AE、弧EF和線段FB連接而成的連續(xù)曲線（根據(jù)對稱性，AE，弧EF，F(xiàn)B連接而成的連續(xù)曲線也是）。EFEF以上只是一種猜測，現(xiàn)在來證明這一猜測是正確的。為此，先介紹一下凸集與凸集的性質(zhì)。定義2.1（凸集）稱集合 R為凸集，若x1、x2R及0，1，總有x1+（1+）x2R。即若x1、x2R，則x1、x2的連線必整個地落 在R中。定理2.2（分離定理）對平面中的凸 集R與R外的一點K，存在

15、直線 l , l 分離R與K，即R與K分別位于 l 的兩側(cè)（注：對一般的凸 集R與R外的一點K，則存在超平面分 離R與K），見圖。klR下面證明猜想猜測證明如下：（方法一）顯然， 由AE、EF、FB及AE，EF，F(xiàn)B圍成的區(qū)域 R是一凸集。利用分離定理易證最短徑不可能經(jīng)過R外的點，若不然，設 為最短路徑，過R外的一點M，則必存在直 線l分離M與R，由于路徑是連續(xù)曲線，由A沿到M，必交l于M1，由M沿到B又必交l于M2。這樣，直線 段M1M2的長度必小于路 徑M1MM2的長度，與是A到B的最短路徑矛盾，至此，我們已證明最短路徑必在凸集R內(nèi)。不妨設路徑經(jīng)湖的上方到達B點，則弧EF必在路徑F上，又直

16、線段AE是由A至E的最短路徑，直線FB是由F到B的最短路徑，猜測得證。ABOrEFEFM1M2Ml還可用微積分方法求弧長，根據(jù)計算證明滿足限止條件的其他連續(xù)曲線必具有更大的長度；此外，本猜測也可用平面幾何知識加以證明等。 根據(jù)猜測不難看出， 例5中的條件可以大大放松，可以不必 設AB過圓心，甚至可不必設湖是圓形的。例如對 下圖，我們可斷定由A至B的最短路徑必 為l1與l2之一，其證明也不難類似給出。 ABl1l2D到此為止，我們的研討還只局限于平面之中，其實上述猜測可十分自然地推廣到一般空間中去。1973年，J.W.Craggs證明了以上結(jié)果：若可行區(qū)域的邊界是光滑曲面。則最短路徑必由下列弧組

17、成，它們或者是空間中的自然最短曲線，或者是可行區(qū)域的邊界弧。而且，組成最短路徑的各段弧在連接點處必定相切。例6 一輛汽車停于 A處并垂直于AB方向，此汽車可轉(zhuǎn)的最小圓半徑為 R，求不倒車而由 A到B的最短路徑。解（情況1）若|AB|2R，最短路徑由 弧AC與切線BC組成（見圖 ）。（情況2）若|AB|0為推力，fS，故由連續(xù)函數(shù)的性質(zhì)存在 某TT，S（T）=S但這一結(jié)果與=(t)是最優(yōu)方案下的車速的假設矛盾，因為用我們猜測的推車方法推車，只 需T時間即可將車推到修車處， 而T0可推出 0的置換矩陣P步2 確定 步3 取 ，用 代替步4 若 =0，停；否則，返回步1。例2. 為方便起見，我們來分

18、解一個元素均為非負整數(shù)的3階雙隨機矩陣，（由Birkhoff定理，r5）解：取 ，=min 1, 3, 3 =1分解成，再取 因min 5, 5, 3 = 3，又有，取 于是又有 易得分解結(jié)果為：尚需解決的問題是如何求P，使得Pij0必有 。讀者不難發(fā)現(xiàn)，此問題可以通過求解一個兩分圖上的最大流（或最大匹配）來實現(xiàn)，計算量為O(n4)，是多項式時間可解的。具體方法為：作一兩分圖，若 ，則作邊（i, j），令邊容量為1，這樣，可作出P的充要條件是該最大流問題的最大流量為n。對例9.33，n=3。由于所有 ，先取 ， P1為 相應的兩分圖為：于是又可求得，相應的兩分圖為： 又可得 ，如此下去，直到作

19、不出P為至，由于 的特殊性質(zhì)及Birkhoff定理，上述分解必能在不超過r= (n1)2 + 1步內(nèi)終止。上述開關(guān)設計方法要求在通訊衛(wèi)星上設置(n1)2 + 1種不同的開關(guān)模式（即Pk），當n稍大時，(n1)2 + 1仍顯得太大而使得使用時不便。例如，當n=41時，(n1)2 + 1=| 60 |。為實用方便，人們研究了限止開關(guān)模式個數(shù)的相應問題。若要求rn，即要求通訊衛(wèi)星上至多設置n種開關(guān)模式，則問題化為令rn，求不超過n個置換矩陣Pk及k，使之滿足： min S.t 為了使任意一對發(fā)射法與接收站之間的傳送均為可能實現(xiàn)的，自然應要求Pk滿足（5.1） （5.2） （右面的矩陣有n2個值為1的

20、分量，每一Pk 恰有n個1分量）故r=n。容易看出，（5.1）隱含著T的每一元素只能被唯一的P復蓋，即T的元素在分解中是不可分割的，這當然是一個好性質(zhì)，使實際操作時較為方便，但可惜的是對一般的雙隨機矩陣，分解很可能無解。例3 若取， ， （注意：T已是雙隨機矩陣，行和列和均為10） 則min S.t 的解為1=3，2=4，3=5。（大于10）而 但等號經(jīng)常并不成立。1985年，F(xiàn)Rendel證明，在給定滿足（5.2）的置換矩陣P1,Pn后，求解問題（5.1）是NP難的，從而不可能存在多項式時間算法，除非P=NP。現(xiàn)要求r2n一種自然而方便的開關(guān)設置為引入兩組各有n個開關(guān)模式的置換矩陣P1,Pn

21、，Q1,Qn，滿足下面的（5.3)式：例如，當n=3時，可令：（注：這種設置方法保持了其內(nèi)在的對稱性，不失為一種明智的做法。）現(xiàn)在，我們來分解例9.33中的雙隨機矩陣 ，令 =，得方程組求出各對角線與反對角線上的三個元素之和，并作一些簡單的消去運算；將矩陣的所有元素相加，可得下面的方程組：注意到（5.3），易證空間 的維數(shù)為 5，故 之一可任取，（稍加注意即可保持非負性），例如，令3=0，求得 ，故有讀者不難驗證，上述方法可推廣到n是奇數(shù)的一般情況。事實，由各對角線元素之和可導出n1個方程，由各反對角線元素之和又可導出n1個方程，加上矩陣所有元素之和導出的等式，共計可導出2 n1個方程，并易知它們是獨立的。另一方面空間 的維數(shù)恰為2 n1，故 之一可任取，而通過方程組解得所有的 ，（只須注意保持其非負性即可）但當n為偶數(shù)時，情況就不大相同了。讓我們先來觀察一下n=4的情況。當n=4時，易見， 具有非常特殊的結(jié)構(gòu)，一般的偶數(shù)階雙隨機矩陣，即使其元素是非負整數(shù)，也無法用Pk、Qk來分解。當 具有上述結(jié)構(gòu)時，能否用Pk和Qk來分解呢？易見，由各對角線元素之和可導出：另外，由