1、利用 Matlab 解決數學問題一、線性規劃求解線性規劃的 Matlab 解法單純形法是求解線性規劃問題的最常用、最有效的算法之一。單純形法是首先由George Dantzig 于 1947 年提出的，近 60 年來，雖有許多變形體已被開發，但卻保持著同樣 的基本觀念。 由于有如下結論： 若線性規劃問題有有限最優解， 則一定有某個最優解是可行 區域的一個極點。基于此，單純形法的基本思路是：先找出可行域的一個極點，據一定規則 判斷其是否最優；若否，則轉換到與之相鄰的另一極點，并使目標函數值更優；如此下去， 直到找到某一最優解為止。 這里我們不再詳細介紹單純形法， 有興趣的讀者可以參看其它線 性規

2、劃書籍。下面我們介紹線性規劃的 Matlab 解法。中線性規劃的標準型為min cT x such that Ax bx基本函數形式為 linprog(c,A,b) ，它的返回值是向量 x 的值。還有其它的一些函數調用形式 (在Matlab 指令窗運行 help linprog 可以看到所有的函數調用形式) ，如：x,fval=linprog(c,A,b,Aeq,beq,LB,UB,X 0,OPTIONS)這里 fval 返回目標函數的值， Aeq 和 beq 對應等式約束 Aeq* x beq ， LB和 UB 分別是變 量 x 的下界和上界， x0 是 x 的初始值， OPTIONS是控制

3、參數。例 2 求解下列線性規劃問題max z 2x1 3x2 5x3x1 x2 x3 72x1 5x2 x3 10 x1,x2,x3 0解 （ i）編寫 M 文件c=2;3;-5;a=-2,5,-1; b=-10;aeq=1,1,1;beq=7;x=linprog（-c,a,b,aeq,beq,zeros（3,1）value=c'*x（ii ）將 M 文件存盤，并命名為。（iii ）在 Matlab 指令窗運行 example1 即可得所求結果。例 3 求解線性規劃問題min z 2x1 3x2 x3x1 4x2 2x3 83x1 2x2 6 x1,x2,x3 0解 編寫 Matlab

4、 程序如下：c=2;3;1;b=8;6;a=1,4,2;3,2,0;x,y=linprog(c,-a,-b,zeros(3,1)二、整數規劃整數規劃問題的求解可以使用 Lingo 等專用軟件。對于一般的整數規劃規劃問題，無法 直接利用 Matlab 的函數， 必須利用 Matlab 編程實現分枝定界解法和割平面解法。 但對于指 派問題等特殊的 0 1 整數規劃問題或約束矩陣 A是幺模矩陣時，有時可以直接利用 Matlab 的函數 linprog 。例 8 求解下列指派問題，已知指派矩陣為3 8 2 10 38729764275842359 10 6 9 10解：編寫 Matlab 程序如下：c

5、=3 8 2 10 3;8 7 2 9 7;6 4 2 7 5 8 4 2 3 5;9 10 6 9 10;c=c(:);a=zeros(10,25);for i=1:5a(5+i,i:5:25)=1;a(i,(i-1)*5+1:5*i)=1;end b=ones(10,1);x,y=linprog(c,a,b,zeros(25,1),ones(25,1)求得最優指派方案為 x15 x 23 x32 x44 x 51 1 ，最優值為 21。三、非線性規劃Matlab 中非線性規劃的數學模型寫成以下形式min f ( x)Ax BAeq x BeqC( x) 0 ,Ceq (x) 0其中 f (

6、x) 是標量函數， A, B, Aeq, Beq是相應維數的矩陣和向量， C(x),Ceq(x) 是非線性 向量函數。Matlab 中的命令是X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)它的返回值是向量 x，其中 FUN 是用 M 文件定義的函數 f (x)；X0 是 x的初始值； A,B,Aeq,Beq 定義了線性約束 A* X B,Aeq* X Beq，如果沒有等式約束，則 A=,B=,Aeq=,Beq= ；LB和 UB是變量 x的下界和上界，如果上界和下界沒有約束， 則 LB=，UB=，如果 x無下界，則 LB=-inf，如果 x

7、 無上界，則 UB=inf；NONLCON是用 M 文件定義的 非線性向量函數 C(x),Ceq(x) ；OPTIONS定義了優化參數，可以使用 Matlab 缺省的參數設置。例 2 求下列非線性規劃問題22min f (x) x1x2 8x12 x2 02x1 x2 2 0x1,x2 0.( i)編寫 M 文件function f=fun1(x);f=x(1)2+x(2)2+8;和M文件function g,h=fun2(x);g=-x(1)2+x(2);h=-x(1)-x(2)2+2; %等式約束( ii)在 Matlab 的命令窗口依次輸入options=optimset;x,y=fmi

8、ncon( 'fun1' ,rand(2,1),zeros(2,1),'fun2' , options)就可以求得當 x1 1,x2 1 時，最小值 y四、圖論兩個指定頂點之間的最短路徑10。問題如下： 給出了一個連接若干個城鎮的鐵路網絡，在這個網絡的兩個指定城鎮間， 找一條最短鐵路線。以各城鎮為圖 G 的頂點，兩城鎮間的直通鐵路為圖G 相應兩頂點間的邊，得圖 G 。對稱為 e的權， 得到賦權圖 G 。G 的G 的每一邊 e，賦以一個實數 w(e) 直通鐵路的長度， 子圖的權是指子圖的各邊的權和。 問題就是求賦權圖 G 中指定的兩個頂點 u0,v0 間的具最小

9、權的軌。這條軌叫做 u 0, v0間的最短路，它的權叫做 u0,v0間的距離，亦記作 d(u0,v0 )。求最短路已有成熟的算法：迪克斯特拉( Dijkstra )算法，其基本思想是按距 u0 從近到 遠為順序，依次求得 u0到G 的各頂點的最短路和距離，直至 v0 (或直至 G 的所有頂點) ， 算法結束。為避免重復并保留每一步的計算信息，采用了標號算法。下面是該算法。(i) 令l(u0) 0，對 v u0，令 l(v)，S0 u0，i 0。(ii) 對每個 v Si ( Si V Si )，用minl(v),l(u) w(uv)u Si代替l(v) 。計算minl(v) ，把達到這個最小值

10、的一個頂點記為 ui 1，令Si 1 Si ui 1。v Si(iii). 若i |V| 1，停止；若 i |V | 1，用 i 1代替i ，轉(ii)。算法結束時，從 u0到各頂點 v的距離由 v的最后一次的標號 l(v)給出。在 v進入 Si 之 前的標號 l(v)叫 T標號， v進入Si時的標號 l(v)叫P標號。算法就是不斷修改各項點的T標號，直至獲得 P 標號。若在算法運行過程中，將每一頂點獲得 P 標號所由來的邊在圖上 標明，則算法結束時， u0 至各項點的最短路也在圖上標示出來了。例 9 某公司在六個城市 c1,c2, , c6中有分公司，從 ci到 cj的直接航程票價記在下述矩

11、陣的 (i , j )位置上。( 表示無直接航路) ，請幫助該公司設計一張城市 c1到其它城市間的票價最便宜的路線圖。050402510500152025150102040201001025252010055102525550用矩陣 an n( n為頂點個數)存放各邊權的鄰接矩陣，行向量pb、 index1 、 index2 、d 分別用來存放 P 標號信息、標號頂點順序、標號頂點索引、最短通路的值。其中分量1 當第 i 頂點已標號pb(i )；0 當第 i 頂點未標號index2( i ) 存放始點到第 i 點最短通路中第 i 頂點前一頂點的序號；d(i) 存放由始點到第 i 點最短通路的值

12、。 求第一個城市到其它城市的最短路徑的 Matlab 程序如下： clear;clc;M=10000; a(1,:)=0,50,M,40,25,10; a(2,:)=zeros(1,2),15,20,M,25; a(3,:)=zeros(1,3),10,20,M; a(4,:)=zeros(1,4),10,25;a(5,:)=zeros(1,5),55;a(6,:)=zeros(1,6);a=a+a'pb(1:length(a)=0;pb(1)=1;index1=1;index2=ones(1,length(a);d(1:length(a)=M;d(1)=0;temp=1;while

13、sum(pb)<length(a)tb=find(pb=0);d(tb)=min(d(tb),d(temp)+a(temp,tb);tmpb=find(d(tb)=min(d(tb);temp=tb(tmpb(1);pb(temp)=1;index1=index1,temp;index=index1(find(d(index1)=d(temp)-a(temp,index1);if length(index)>=2index=index(1);endindex2(temp)=index;endd, index1, index2每對頂點之間的最短路徑計算賦權圖中各對頂點之間最短路徑，顯

14、然可以調用Dijkstra 算法。具體方法是：每次以不同的頂點作為起點， 用 Dijkstra 算法求出從該起點到其余頂點的最短路徑， 反復執行 n 次 這樣的操作，就可得到從每一個頂點到其它頂點的最短路徑。這種算法的時間復雜度為3O(n3) 。第二種解決這一問題的方法是由Floyd R W提出的算法，稱之為 Floyd 算法。假設圖 G 權的鄰接矩陣為A0 ，a11a12a1na21a22a2nA0an1an2ann來存放各邊長度，其中：aii 0 i 1,2, ,n；aiji, j 之間沒有邊，在程序中以各邊都不可能達到的充分大的數代替；aij wij wij 是 i, j 之間邊的長度，

15、 i, j 1,2, ,n 。對于無向圖， A0 是對稱矩陣， aij aji 。Floyd算法的基本思想是： 遞推產生一個矩陣序列 A0,A1, ,Ak, , An ，其中 Ak(i, j)表 示從頂點 vi 到頂點 vj 的路徑上所經過的頂點序號不大于k的最短路徑長度。計算時用迭代公式：Ak(i,j) min( Ak 1(i,j),Ak 1(i,k) Ak 1(k,j)k 是迭代次數， i, j,k 1,2, ,n 。最后，當 k n 時， An 即是各頂點之間的最短通路值。例10 用 Floyd算法求解例 1。矩陣 path用來存放每對頂點之間最短路徑上所經過的頂點的序號。Floyd算法

16、的 Matlab程序如下：clear;clc;M=10000;a(1,:)=0,50,M,40,25,10;a(2,:)=zeros(1,2),15,20,M,25;a(3,:)=zeros(1,3),10,20,M;a(4,:)=zeros(1,4),10,25;a(5,:)=zeros(1,5),55;a(6,:)=zeros(1,6);b=a+a'path=zeros(length(b);for k=1:6for i=1:6for j=1:6if b(i,j)>b(i,k)+b(k,j)b(i,j)=b(i,k)+b(k,j);path(i,j)=k;endendend e

17、nd b, path4.2.1 prim 算法構造最小生成樹設置兩個集合 P和Q,其中 P用于存放 G 的最小生成樹中的頂點，集合 Q存放 G 的最 小生成樹中的邊。令集合 P的初值為 P v1 （假設構造最小生成樹時，從頂點 v1出發）， 集合 Q的初值為 Q。prim 算法的思想是，從所有 p P，v V P的邊中， 選取具有最小權值的邊 pv ，將頂點 v 加入集合 P中，將邊 pv加入集合 Q中，如此不斷重復，直到 P V 時，最小生成樹構造完畢，這時集合 Q 中包含了最小生成樹的所有邊。prim 算法如下：（i） P v1 ，Q（ ii） while P Vpv min（ wpv,

18、p P,v V PP P vQ Q pvend例 11 用 prim 算法求右圖的最小生成樹。Matlab我們用 result 3 n 的第一、二、三行分別表示生成樹邊的起點、終點、權集合。 程序如下：clc;clear;M=1000;a(1,2)=50; a(1,3)=60;a(2,4)=65; a(2,5)=40;a(3,4)=52;a(3,7)=45;a(4,5)=50; a(4,6)=30;a(4,7)=42;a(5,6)=70;a=a;zeros(2,7);a=a+a'a(find(a=0)=M; result=;p=1;tb=2:length(a);while length

19、(result)=length(a)-1temp=a(p,tb);temp=temp(:);d=min(temp);jb,kb=find(a(p,tb)=d);result=result,j;k;d;p=p,k;tb(find(tb=k)=;j=p(jb(1);k=tb(kb(1);end result4.2.1 Kruskal 算法構造最小生成樹科茹斯克爾( Kruskal)算法是一個好算法。 Kruskal 算法如下 :(i)選 e1 E(G) ，使得 w(e1) min 。(ii)若e1,e2, ,ei 已選好，則從 E(G) e1,e2, ,ei 中選取 ei 1，使得 G e1,e2

20、, ,ei ,ei 1 中無圈，且 w(ei 1) min 。(iii)直到選得 e 1 為止。例 12 用 Kruskal 算法構造例 3 的最小生成樹。我們用 index2 n 存放各邊端點的信息， 當選中某一邊之后， 就將此邊對應的頂點序號中 較大序號 u改記為此邊的另一序號 v，同時把后面邊中所有序號為 u 的改記為 v 。此方法的 幾何意義是：將序號 u 的這個頂點收縮到 v 頂點， u 頂點不復存在。后面繼續尋查時，發現 某邊的兩個頂點序號相同時，認為已被收縮掉，失去了被選取的資格。Matlab 程序如下：clc;clear;M=1000;a(1,2)=50; a(1,3)=60;

21、a(2,4)=65; a(2,5)=40;a(3,4)=52;a(3,7)=45;a(4,5)=50; a(4,6)=30;a(4,7)=42;a(5,6)=70;i,j=find(a=0)&(a=M);b=a(find(a=0)&(a=M);data=i'j'b'index=data(1:2,:);loop=max(size(a)-1;result=;while length(result)<looptemp=min(data(3,:);flag=find(data(3,:)=temp);flag=flag(1);v1=data(1,flag);

22、v2=data(2,flag);if index(1,flag)=index(2,flag)result=result,data(:,flag);endif v1>v2index(find(index=v1)=v2;elseindex(find(index=v2)=v1;enddata(:,flag)=;index(:,flag)=;endresult旅行商( TSP)問題一名推銷員準備前往若干城市推銷產品， 然后回到他的出發地。 如何為他設計一條最短 的旅行路線 (從駐地出發，經過每個城市恰好一次，最后返回駐地)這個問題稱為旅行商問 題。用圖論的術語說，就是在一個賦權完全圖中，找出一個

23、有最小權的 Hamilton 圈。稱這 種圈為最優圈。 與最短路問題及連線問題相反， 目前還沒有求解旅行商問題的有效算法。 所 以希望有一個方法以獲得相當好(但不一定最優)的解。一個可行的辦法是首先求一個 Hamilton 圈 C ，然后適當修改 C 以得到具有較小權的另 一個 Hamilton 圈。修改的方法叫做改良圈算法。設初始圈 C v1v2 vnv1 。( i)對于 1 i 1 j n ，構造新的 Hamilton 圈 :Cij v1v2 vivjvj 1vj 2 vi 1vj 1vj 2 vnv1 ,它 是 由 C 中 刪 去 邊 vivi 1 和 vjvj 1 ， 添 加邊 viv

24、j 和 vi 1vj 1 而得 到 的 。 若 w(vivj) w(vi 1vj 1) w(vivi 1) w(vjvj 1)，則以 Cij代替C，Cij叫做 C的改良圈。ii)轉( i)，直至無法改進，停止。用改良圈算法得到的結果幾乎可以肯定不是最優的。 為了得到更高的精確度， 可以選擇 不同的初始圈，重復進行幾次算法，以求得較精確的結果。這個算法的優劣程度有時能用 Kruskal算法加以說明。 假設 C 是 G 中的最優圈。 則對于 任何頂點 v ，C v是在 G v中的 Hamilton 軌，因而也是 G v的生成樹。 由此推知： 若 T 是 G v中的最優樹，同時 e和 f 是和 v關聯的兩條邊，并使得 w(e) w(f) 盡可能小，則 w(T) w(e) w( f) 將是 w( C)的一個下界。這里介紹的方法已被進一步發展。 圈的修改過程一次替換三條邊比一次僅替換兩條邊更 為有效；然而，有點奇怪的是，進一步推廣這一想法，就不利了。例 13 從北京( Pe)乘飛機到東京 (T)、紐約 (N)、墨西哥城 (M) 、倫敦 (L)、巴黎 (Pa)五城 市做旅游， 每城市恰去一次再回北京， 應如何安排旅游線， 使旅程