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, outp

2、ut = bintprog(.),數(shù)學模型：MIN lambda S.T. F(x)-weight* lambda =goal(達到目標） Ax=b（線性不等式約束） Aeqx=beq（線性等式約束） C(x)=0（非線性不等式約束） Ceq(x)=0（非線性等式約束） lb=x=ub F=f1(x),f2(x),為多目標的目標函數(shù)； F與C(x),Ceq(x)都是通過function來定義； 命令格式： x = fgoalattain(fun,x0,goal,weight) x = fgoalattain(fun,x0,goal,weight,A,b) x = fgoalattain(fun

3、,x0,goal,weight,A,b,Aeq,beq) x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub),二、多目標規(guī)劃的MATLAB求解,命令格式： x = fgoalattain(fun,x0,goal,weight) x = fgoalattain(fun,x0,goal,weight,A,b) x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq) x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub) x = fgoalattain(

4、fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon) x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,. lb,ub,nonlcon,options) x,fval = fgoalattain(.) x,fval,attainfactor = fgoalattain(.) x,fval,attainfactor,exitflag = fgoalattain(.) x,fval,attainfactor,exitflag,output = fgoalattain(.) x,fval,attainfactor,ex

5、itflag,output,lambda = fgoalattain(.),二、多目標規(guī)劃的MATLAB求解,x = fgoalattain(myfun,x0,goal,weight,A,b,Aeq,beq,. lb,ub,mycon) where mycon is a MATLAB function such as function c,ceq = mycon(x) c = . % compute nonlinear inequalities at x. ceq = . % compute nonlinear equalities at,二、多目標規(guī)劃的MATLAB求解,x = fgoala

6、ttain(myfun,x0,goal,weight) where myfun is a MATLAB function such as function F = myfun(x) F = . % Compute function values at x.,有關優(yōu)化參數(shù)設置： options = optimset(GradObj,on)目標函數(shù)的梯度方向參數(shù)設置為on時，用下列函數(shù)定義： function F,G = myfun(x) F = . % Compute the function values at x if nargout 1 % Two output arguments G =

7、 . % Gradients evaluated at x End The gradient consists of the partial derivative dF/dx of each F at the point x.,二、多目標規(guī)劃的MATLAB求解,二、多目標規(guī)劃的MATLAB求解,有關優(yōu)化參數(shù)設置： options = optimset(GradConstr,on）約束條件的梯度方向參數(shù)設置為on時，用下列函數(shù)定義： function c,ceq,GC,GCeq = mycon(x) c = . % Nonlinear inequalities at x ceq = . % No

8、nlinear equalities at x if nargout 2 % Nonlcon called with 4 outputs GC = . % Gradients of the inequalities GCeq = . % Gradients of the equalities End 注意：一般 weight=abs(goal),模型：x=(A+BKC)x+Bu，設計K滿足目標： Y=Cx 1）循環(huán)系統(tǒng)的特征值（由命令eig(A+B*K*C)確定）的目標為goal = -5,-3,-1 2）K中元素均在-4,4中; 設特征值的weight= abs(goal),定義目標函數(shù)F如

9、下： 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 = -1 -1; -1 -1; % Initialize controller matrix goal = -5 -3 -1; % Set goal values for the eigenvalues weight = abs(goal) % Set weight for

10、same percentage lb = -4*ones(size(K0); % Set lower bounds on the controller ub = 4*ones(size(K0); % Set upper bounds on the controller options = optimset(Display,iter); % Set display parameter K,fval,attainfactor = fgoalattain(K)eigfun(K,A,B,C). goal,weight,lb,ub,options),二、舉例-有關循環(huán)控制系統(tǒng)優(yōu)化問題,運行結(jié)果如下 Ac

11、tive constraints: 1 2 4 9 10 K = -4.0000 -0.2564 -4.0000 -4.0000 fval = -6.9313 -4.1588 -1.4099 attainfactor = -0.3863,二、舉例-有關循環(huán)控制系統(tǒng)優(yōu)化問題,如果至少保證38.63%的目標精確匹配，設置GoalsExactAchieve參數(shù)值為3 options = optimset(GoalsExactAchieve,3); K,fval,attainfactor = fgoalattain(. (K)eigfun(K,A,B,C),K0,goal,weight,lb,ub,.

12、options) After about seven iterations, a solution is K = -1.5954 1.2040 -0.4201 -2.9046 fval = -5.0000 -3.0000 -1.0000 attainfactor = 1.0859e-20表明目標已完全匹配,二、舉例-有關循環(huán)控制系統(tǒng)優(yōu)化問題,謝謝！,初等模型舉例,常見類型,定性模型 經(jīng)驗公式（擬合、插值） 量綱分析 比例模型,2.1 崖高的估算,假如你站在崖頂且身上帶著一只具有跑表功 能的計算器，你也許會出于好奇心想用扔下 一塊石頭聽回聲的方法來估計山崖的高度， 假定你能準確地測定時間，你又怎

13、樣來推算 山崖的高度呢，請你分析一下這一問題。,方法一,我學過微積分，我可以做 得更好，呵呵。,令k=K/m，解得,代入初始條件 v(0)=0，得c=g/k，故有,再積分一次，得：,若設k=0.05并仍設 t=4秒，則可求 得h73.6米。,聽到回聲再按跑表，計算得到的時間中包含了 反應時間,進一步深入考慮,不妨設平均反應時間 為0.1秒 ，假如仍 設t=4秒，扣除反應時間后應 為3.9秒，代入 式，求得h69.9米。,多測幾次，取平均值,再一步深入考慮,2.2 錄像帶還能錄多長時間,錄像機上有一個四位計數(shù)器，一盤 180分鐘 的錄像帶在開始計數(shù)時為 0000，到結(jié)束時計 數(shù)為1849，實際走

14、時為185分20秒。我們從 0084觀察到0147共用時間3分21秒。若錄像 機目前的計數(shù)為1428，問是否還能錄下一個 60分鐘的節(jié)目？,又 因和 得,積分得到,即,從而有,此式中的三個參數(shù)W、v和r均不易精確測得，雖然我們可以從上式解出t與n的函數(shù)關系，但效果不佳，故令 則可將上式簡化為：,故,t= an2+bn,上式以a、b為參數(shù)顯然是一個十分明智的做法，它為公式的最終確立即參數(shù)求解提供了方便。將已知條件代入，得方程組：,從后兩式中消 去t1，解得a=0.0000291, b=0.04646,故t=0.0000291 n2+0.04646n，令n=1428，得到t=125.69（分）由于

15、一盒錄像帶實際可錄像時間為185.33分，故尚可錄像時間 為59.64分，已不能再錄下一個60分鐘的節(jié)目了。,2.3最短路徑與最速方案問題,例5（最短路徑問題）,設有一個半徑為 r 的圓形湖，圓心為 O。A、 B 位于湖的兩側(cè)，AB連線過O，見圖。 現(xiàn)擬從A點步行到B點，在不得進入湖中的限 制下，問怎樣的路徑最近。,以上只是一種猜測，現(xiàn)在來證明這一猜測是正確的。為此，先介紹一下凸集與凸集的性質(zhì)。,下面證明猜想,猜測證明如下：,還可用微積分方法求弧長，根據(jù)計算證明滿足限止條件的其他連續(xù)曲線必具有更大的長度；此外，本猜測也可用平面幾何知識加以證明等。,到此為止，我們的研討還只局限于平面之中，其實上

16、述猜測可十分自然地推廣到一般空間中去。1973年，J.W.Craggs證明了以上結(jié)果：,例6 一輛汽車停于 A處并垂直于AB方向，此 汽車可轉(zhuǎn)的最小圓半徑為 R，求不倒車而由 A到B的最短路徑。,例7 駕駛一輛停于A處與AB成1角度的汽 車到B處去，已知B處要求的停車方向必須 與 AB成2角，試找出最短路徑（除可轉(zhuǎn) 的最小圓半徑為R外，不受其他限止）。,最速方案問題,例8 將一輛急待修理的汽車由靜止開始沿一 直線方向推至相隔 S米的修車處，設阻力不 計 ，推車人能使車得到的推力 f 滿足: -BfA , f0為推力，f0為拉力。問怎樣 推車可使車最快停于修車處。,此問題為一泛函極值問題，求解十

17、分困難，為得出一個最速方案。我們作如下猜測：,邏輯模型,例1 擬將一批尺寸為124的的商品裝入尺寸為666的正方體包裝箱中，問是否存在一種裝法，使裝入的該商品正好充滿包裝箱。,解 將正方體剖分成27個222的小正方體，并 按下圖所示黑白相間地染色。,再將每一222的小正方體剖分成111的小正方體。 易見，27個222的正方體中，有14個是黑的，13個是白的（或13黑14白），故經(jīng)兩次剖分，共計有112個111的黑色小正方體和104個111的白色小正方體。 雖然包裝箱的體積恰好是商品體積的27倍，但容易看到，不論將商品放置在何處，它都將占據(jù)4個黑色和4個白色的111小正方體的位置，故商品不可能充

18、滿包裝箱。,德國著名的藝術家Albrecht Drer(1471-1521)于1514年曾鑄造了一枚名為“Melencotia I”的銅幣。令人奇怪的是在這枚銅幣的畫面上充滿了數(shù)學符號、數(shù)字及幾何圖形。這里，我們僅研究銅幣右上角的數(shù)字問題,所謂的魔方是指由1n2這n2個正整數(shù)按一定規(guī)則排列成的一個n行n列的正方形 。n稱為此魔方的階 。,Drer魔方：4階，每一行之和為34，每一列之和為34，對角線（或反對角線）之和是34，每個小方塊中的數(shù)字之和是34，四個角上的數(shù)字加起來也是34,什么是Drer魔方,多么奇妙的魔方！,銅幣鑄造時間：1514年,構造魔方是一個古老的數(shù)學游戲，起初它還和神靈聯(lián)系

19、在一起，帶有深厚的迷信色彩。傳說三千二百多年前（公元前2200年），因治水出名皇帝大禹就構造了三階魔方（被人們稱“洛書”），至今還有人把它當作符咒用于某些迷信活動，大約在十五世紀時，魔方傳到了西方，著名的科尼利厄斯阿格里帕（1486-1535）先后構造出了39階的魔方 。,如何構造魔方,奇數(shù)（不妨n=5）階的情況,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,魔方數(shù)量隨階數(shù)n增長的速度實在是太驚人了！,同階魔方的個數(shù),允許構成魔方的數(shù)取任意實數(shù),允許取實數(shù),n階魔方A、B，任意實數(shù)、,A+B是n階魔方,具有指

20、定性質(zhì)的魔方全體構成一個線性空間,問題已發(fā)生了實質(zhì)性變化,注：刻畫一個線性空間只需指出它的維數(shù)并求出此線性空間的一組基底,松馳問題的討論,1在第一行中共有4種取法，為保持上述性質(zhì)的成立，第二行中的1還有兩種取法。當?shù)诙械?也取定后，第三行與第四行的1就完全定位了，故一共可作出8個不同的最簡方陣，稱之為基本魔方并記之為Q1， ，Q8,仍以4階方陣為例。 令R為行和，C為列和，D為對角線和，S為小方塊和 定義0-方：R=C=D=S=0 定義1-方：R=C=D=S=4,R=C=D=S=1的方陣構成的線性空間具有什么樣的性質(zhì)？,類似于構造n維歐氏空間的標準基，利用0和1我們來構造一些R=C=D=S=

21、1的最簡單的方陣。,顯然， Drer空間（簡稱D空間）中任何一個元素都可以用Q1，Q2，Q8來線性表示，但它們能否構成D空間的一組基呢？,等號兩邊對應元素相比較，得r1=r2=r7=0， 所以 是線性無關,是D空間的最小生成集。,令D 即解方程組：,=,解得 D=,研究Albrecht Drer鑄造的銅幣,2004年浙江大學數(shù)學建模競賽 （B題）通訊衛(wèi)星上的開關設置 地面上存在著n個接收站與n個發(fā)送站，而在通訊衛(wèi)星上則設置了若干種開關模式。開關模式可用矩陣P=(pij)來表示，若衛(wèi)星可接收發(fā)送站i發(fā)射的信息并將信息傳送回地面的接收站j時，矩陣中的元素pij =1，否則pij =0。通訊衛(wèi)星上的

22、接收發(fā)送任務也可以用一個矩陣T=（tij）來表示，其元素tij為需經(jīng)通訊衛(wèi)星傳遞的由i發(fā)點發(fā)送到j接收點的信息量的傳送時間長度。由于技術上的原因，當發(fā)送站i在發(fā)送給接收站j信息時，它不能同時發(fā)送給別的接收站信息；同樣，當接收站j在接收發(fā)送站i的信息時，也不能同時接收其他發(fā)送站發(fā)送的信息。你的任務是：,設計一組開關模式，k=1, ,r（注：r應當盡可能?。?，使得對任意給定的任務矩陣T，衛(wèi)星開關設置均能完成要求的發(fā)接收任務。 設計一個算法，在發(fā)接收任務T給出后，可根據(jù)你設計的開關模式（k=1,r）求出各開關的使用時間k，使得在完成預定傳送任務的前提下使用各開關模式的總時間最短。 同樣由于技術上的原

23、因，開關模式的總數(shù)r有一個上限。當需要傳送的任務數(shù)數(shù)量較大時，仍無法分派任務。請你想一些辦法來解決這一困難，（當然，這時你可能要作出一些犧牲，即傳送時間可能會增加一些）。,問題及模型,問題的標準形式為：在地面上存在著n個收站與n個發(fā)戰(zhàn)，而在通訊衛(wèi)星上則設置了若干種開關模式。開關模式可用矩陣P=(pij)來表示，若衛(wèi)星可接收發(fā)射站i發(fā)射的信息并將信息傳送回地面的接收站j時，矩陣元素pij =1，否則pij =0。通訊衛(wèi)星的接發(fā)任務也可用一矩陣T=（tij）來表示，其元素tij為需經(jīng)通訊衛(wèi)星傳遞的由i發(fā)點發(fā)送到j接受點的信息量的傳送時間長度。問題要求求r并設計一組開關模式Pk，k=1, ,r及模式

24、Pk的使用時間k，使得在完成預定傳送任務的前提下各開關模式使用的總時間最短，即要求求解下面的問題：,例1 設,這是一個有3個發(fā)送站與3個接收站的實例，tij在矩陣中已給出，例如由發(fā)站1傳送到收站1的通訊量為3單位時間等。,分析 容易看出，三個發(fā)站需傳送的時間分別為6、5、5；而三個收站需接收的時間分別為6、3、7。為完成全部傳送任務，通訊衛(wèi)星總傳送時間至少應為7單位時間，即的下界為7。,由于技術上的原因，當發(fā)站i在發(fā)送給收站j信息時，它不能同時發(fā)送給別的收站信息；同樣，當收站j在接收發(fā)站i的信息時，也不能同時接收其他發(fā)站發(fā)送的信息。這一要求說明，任一開關模式Pk應具有以下性質(zhì)：（1）Pk的每一

25、行中有且只有一個1，每一列中也有且只有一個1；（2）所有的1均位于不同的行列中。,滿足（1）、（2）的矩陣 被稱為置換矩陣，n階置換矩陣Pk共有n!個，當n較大時，我們不可能在通訊衛(wèi)星上設置這么多種不同的開關模式。因而，為了設計出切實可行的開關模式，我們還得另想辦法。,（設計方法1）,注意到Pk每行（或列）元素之和均為1，故不管如何指派開關的使用時間（即不論如何取k），矩陣,均具有某些特殊的性質(zhì)，例如其行和（及列和）均為同一常數(shù)。這樣的矩陣構成一個線性空間（參見邏輯模型第一節(jié) Drer魔方），為減少開關模式的種類，可取此空間的一組基底作為開關模式。在使用這種開關模式時，無論T的元素tij怎么取

26、，通訊衛(wèi)星對每一發(fā)（收）點的開通時間總和是恒定的。在這種開關模式下，可按如下方式指派各開關模式的使用時間：,步1 先將T改變?yōu)?， 滿足：,將T化為 的方法一般有無窮多種，如可如下化法：,令 事實上， ，（即通訊衛(wèi)星傳送總時間的下界）。,令,其中,用這種方法化例中的T，得到,的任一行（或列）中元素之和均為7。,定義1 稱行和、列%和均相等的矩陣為雙隨機矩陣（Doubly stochastic matrix）,定理1 （Birkhoff定理，1944）任一n階雙隨機矩陣均可寫成至多 （n1）2+1個置換矩陣的非線性組合。,的分解可如下進行：,步1 選取由Pij0可推出 0的置換矩陣P,步2 確定

27、,步3 取 ，用 代替,步4 若 =0，停；否則，返回步1。,例2. 為方便起見，我們來分解一個元素均為非負整數(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。由于所有 ，先取,

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

29、之間的傳送均為可能實現(xiàn)的，自然應要求 Pk滿足,（5.1）,（5.2）,（右面的矩陣有n2個值為1的分量，每一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,一種自然而方便的開關設置為引入兩組各有n個開關模式的置換矩陣P1,Pn，Q1,Qn，滿足下面的（5.3)式：,例如，當n=3時，可令：,（注：這種設置方法保持了其內(nèi)在的對稱性，不失為一種明智的做法。）,現(xiàn)在，我們來分解例9.33中的雙隨機矩陣 ，令 = ，得方程組,求出各對角線與反對角線上的三個元素之和，并作一些簡單的消去運算； 將矩陣的所有元素相加，可得下面的方程組：,注意到（5.3），易證空間 的維數(shù)為 5， 故 之一可任取，（稍加注意即可保持非負性）， 例如，令3=0，求得 ，故有,讀者不難驗證，上述方法可推廣到n是奇數(shù)的一般情況。事實，由