数学建模社区-数学中国

标题: 数学建模--图与网络（2） [打印本页]

作者: 杨利霞 时间: 2018-10-30 11:25
标题: 数学建模--图与网络（2）

最大流：

注意Matlab 中的最大流问题是必须是单源和单汇问题，因此这里需要构建虚拟的源点S和汇点G。
clc,clear
a = zeros(9,9);
a(1,[2:4]) = [20,20,100];
a(2,[5 6 8]) = [30,10,40];
a(3,[7 8]) = [10,50];
a(4,[5:8]) = [20,10,40,5];
a([5:8],9) = [20,20,60,20];
a = sparse(a);
[b,c] = graphmaxflow(a,1,9)

最大流拓展：最小费用最大流

仅仅是在求得最大流之后进行对最小费用的求解：

clc,cleara = zeros(5);a(1,[2 3]) = [10 8;a(2,[4,5]) = [2 7;a(3,[2 4]) = [5 10;a(4,5) = 4;a = sparse(a);[b c] = graphmaxflow(a,1,5)

最大流量为11

自定义Matlab代码：

最小费用求解& Y+ n- W& y' k& H2 C
- f- m, r/ [% j8 z
Lingo:" i% ?3 W8 F. H  m5 \: X  k" V; N

4 @  ?: [  G2 H; o! R% b1 tmodel:3 m# ^7 i* ~4 P& |' i4 X
sets:( L2 j3 ^& p+ u( h+ \
nodes/s,1,2,3,t/:d;, G  ]. C: n; B" [! U$ L
arcs(nodes,nodes)/s 1,s 2,1 3,1 t,2 1,2 3,3 t/:b,c,f;% Q+ b7 j  F: `" U/ z# e3 J$ d: T
endsets
9 j; x& G2 [; L+ Q9 s8 f& R8 ~data:6 K4 g3 K3 ?/ `, D/ t& l
b = 4 1 6 1 2 3 2;
7 L/ F, }1 L6 K& i$ ?6 X( Pc = 10 8 2 7 5 10 4;
* s0 F2 b- `8 C9 Q6 s. M+ Zd = 11 0 0 0 -11;6 p, ^0 ~1 S+ q7 h! {( F
enddata
( p: f2 u  {+ W" d8 Qn = @size(nodes);7 R+ H* D! e" a  `3 V# N
min = @sum(arcs:b*f);
: m3 n/ l& |2 [6 h- s  a@for(nodes(i)sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i)) = d(i));+ Y" J/ |( |( u$ F( \2 G3 F9 m
@for(arcsbnd(0,f,c));
3 d: N; o1 x6 [; n4 L& cend
) A& j* K! |# z" u- J, h! r1
- R  ?- F# {7 C/ X. p2" T. [# k' K0 `; b% h2 K
3
! S1 n4 i! Q: f% X3 c4
6 |2 \8 }$ l$ m0 A  s9 A$ ^$ W7 |5* H5 V! p7 M; W) g7 P( L. o( a
6; m+ s- y8 @; T. y
75 T/ k" d, a% j; |0 B- s4 _
84 o% h  F9 I% N  ^6 C* Z
9
5 Z, ?9 {* T6 F# W  [2 t3 h10( f! X! b! ?/ z3 `' ]/ i' D
11
1 C$ T" z4 X; }1 k4 k8 O12) U( |$ D% q, F' |, y" M
13
) `! M% h- g' H& S& A9 t, k14
2 }# Q# b6 X: S2 }15
+ R! v$ o1 C6 _9 @8 }+ eMatlab实现：
1 N, l$ s2 E7 c6 r5 x$ r) ?/ E2 `5 W" d

, V- J3 O" ?! v' R- Gn = 5;
% s- T3 f, t3 ~- ?( T* H%弧容量
) V+ C; U3 j0 [7 U; p$ \9 D+ e3 t0 ?a = zeros(5);
% o( k0 K, p8 Y0 D. n6 I; b- Na(1,[2 3]) = [10 8];
5 @# N% c4 `0 q$ ba(2,[4,5]) = [2 7];/ _/ W' b& [: c3 |4 W' m: F
a(3,[2 4]) = [5 10];
7 D0 m! `8 R- ta(4,5) = 4;
3 G# y9 _* f8 W. G% |C = a;( M6 }( ~8 r! B# u6 z/ ~) V
%C = [0 15 16 0 0;0 0 0 13 14;0 11 0 17 0;0 0 0 0 8;0 0 0 0 0];
; B+ g/ e% b8 v0 p%弧上单元的费用( p2 x* j" ~$ W  M9 l
a(1,[2 3]) = [4 1];5 o: \8 w% F9 K0 o+ m
a(2,[4,5]) = [6 1];  Q1 t1 {6 O/ T3 A: t) ?5 H
a(3,[2 4]) = [2 3];( U7 v' G+ b# ^
a(4,5) = 2;* n& j8 d" S$ [) e1 V9 l1 X
b = a;
. W& j3 o$ O& ]1 n1 A%b = [0 4 1 0 0;0 0 0 6 1;0 2 0 3 0;0 0 0 0 2;0 0 0 0 0];
* s$ n( Z0 R5 z: l%wf表示最大流量。wf0表示预定的流量值
3 \0 J4 p/ R1 _3 T( wwf = 0;& y- K: |/ A% d* m% D/ S3 J% \
wf0 = Inf;
5 A+ i& l4 @9 l$ p%取初始可行流f为零流
' g, ^' }% i$ Q4 m3 |" ffor i = 1:n- s: r4 V$ s: d& ]7 U; F" `
for j = 1:n
2 L. g: @5 x: F0 I( O! x       f(i,j) = 0;
$ x* X" @5 m! @: w$ o end
# H4 S5 g0 g4 T* O* q' Dend2 k: ]5 g) F' q+ i
while (1)0 R' M, u% q, S* b) ?
%构造有向赋权图4 R8 d- f& H+ ?: H' R
for i = 1:n
! J, |7 j- P0 B; s4 ]! y& S       for j = 1:n
$ S- y- {# _3 F0 _) c* ]          if j~=i) o" ^1 X$ ]! E1 B, _
            a(i,j) = Inf;
- T1 ~3 d8 P* w& U% a4 m          end0 c2 b  U* X& Z* s
      end
5 L& o$ H( g' y, x& c# C3 i end
/ h' Q- T# h! } for i = 1:n6 g- W' `7 E3 H
      for j = 1:n/ R/ ^7 \- g2 y3 j
         if C(i,j) > 0 && f(i,j) == 0+ v& o( ^; Q9 |- L( X1 Z* @
            a(i,j) = b(i,j);
" s7 q" o/ ~+ E( r- E9 z: o  }          elseif C(i,j) > 0 && f(i,j) == C(i,j)  A5 {" D/ j' u# c& }: \) e7 M
            a(j,i) = -b(i,j);; F9 c  y* e8 C! _. g. _/ e
         elseif C(i,j) > 0( B+ ~- k4 b( U$ w0 P# \6 r
            a(i,j) = b(i,j);
5 D- G1 L" y6 c' V# c             a(j,i) = -b(i,j);% }$ x6 b% C5 }
         end+ f/ N- K& ^8 D
      end
$ J: L% {% s7 F. q end
) Z9 g8 h+ U0 C1 v2 a" y %使用ford算法求最短路，赋初值% o* X8 |3 a. k  r4 L# D
for i = 2:n
2 H3 r4 |8 M2 ]( n, x' K' J       p(i) = Inf;1 X, b# l# E5 r7 P7 w
      s(i) = i;
! w+ y& g2 d+ x+ ] end
8 i& G% v) b! T: h1 L %求有向赋权图vs到vt的最短路,赋初值9 D: ]8 ~- N* Z3 _2 N
for k = 1:n
* u7 Q5 c3 }% W- D/ n       pd = 1;& Z6 l/ b& [( }
      for i = 2:n( F9 f' w2 s3 _6 [# h. M# C9 `
         for j = 1:n% Y7 c$ x! W4 I6 Y6 a
            if p(i) > p(j) + a(j,i)" n" g- @- |) ~% u& y8 u" D  W
                  p(i) = p(j) + a(j,i);6 j* |) l6 n& K( U- i3 |$ Q1 [: @
                  s(i) = j;
" N# Z* A  l4 q                   pd = 0;
6 a) h+ }; r3 R7 s7 u             end9 S9 v8 N+ |4 t- S' h' b; S) X  s
         end* F8 @! p" r5 v# R3 d+ \
      end1 E  l7 Q' l# t3 a4 F4 {
      %求最短路的Ford算法结束
) M9 D. l  T3 M. d$ z, y2 n' C4 c       if pd
+ Q; z6 i) s, z7 W+ E5 \          break;
- k, X+ y3 t/ T# _/ A7 {       end7 S+ x9 x" l! m4 B9 u, ]
end5 P- w/ q6 H8 _' u2 h" d* E
%不存在vs到vt的最短路，算法终止，注意在求最小费用最大流时构造有向赋权图中不会含有负权回路，不会出现k = n8 C* n, M: d2 z" g& k4 y. _3 _( ?
if p(n) == Inf
- ~; r! ~/ B5 K* }( u5 D, ?       break;
' e8 x) q- n1 Z3 ]! J( M end
. U$ S$ I  ]6 T$ o( v %进入调整过程,dvt表示调整量
$ X' }$ q, M( n, z  K: j dvt = Inf;
! p' i- O  J1 }4 Y dvtt = Inf;
8 Q( p( _% F! @5 e0 W+ l) k" d t = n;
# S9 }+ [2 U  {; {5 J while(1)
. n' ^6 I8 D1 {- D3 n% A       if a(s(t),t) > 0, W9 e  Y' S( d5 M; f* Y6 d
         %前向弧调整量
) M; q7 \% Z7 ]! ^- x          dvtt = C(s(t),t)-f(s(t),t);
6 Z( ?7 V$ K# E0 y* [- u5 c          %后向弧调整量
6 e! v, J/ G6 g4 L. k7 b( i9 d# w       elseif a(s(t),t) < 0+ V$ m  M8 i1 c7 R# @
         dvtt = f(t,s(t));  M5 g1 e! E/ K: f* |" o
      end5 x, ~% @; H' `) Z3 P) u' ?, a
      if dvt > dvtt9 F" q; e4 X3 {) D  X/ y/ O- |4 o' y+ k
         dvt = dvtt;7 n8 u0 [8 A  A+ i0 U
      end3 k- `7 g9 y$ y
      %当t的标号为vs时，终止计算调整量! t1 c' E8 w& J) ?. t) _
      if s(t) == 1& G# ?: z, ~( o0 H! D
         break;: o/ g6 j2 Y7 U5 u* N2 S& B: T
      end
) X. w4 c2 P8 {2 U4 r( g0 _' U9 l       %继续调整前一段弧上的流f
: Q/ Z, Z9 n3 v! [       t = s(t);
- k( g# t% U9 [( b* `5 U end4 q' \& E5 f: R& K
pd = 0;
: C3 l$ F% H0 U/ ?0 ?9 M %如果最大流量大于或等于预定的流量值1 p4 ^; v$ t5 f* v3 }& s* H; M
if wf + dvt >= wf0
/ g: F- K' W% Y3 o, b, P       dvt = wf0 - wf;# A/ c, F7 Q' E: K: b
      pd = 1;
; y, r9 F9 b6 j: O% {! P end
& ^, V/ s' y" H6 x/ y" L! B1 |8 y t = n;' L) L8 x& u1 L5 |% U3 L
%调整过程% o8 r- R& }0 I0 W
while(1)" |: n& H- U* p# U; {+ {8 [
if a(s(t),t) > 0
' ?4 ~0 j2 a3 c+ [       %前向弧调整  [: r. C0 z; d, c  t* ~) O
      f(s(t),t) = f(s(t),t) + dvt; % f& t8 B+ ]  I8 z. \/ x
elseif a(s(t),t) < 0" f  w& M- Q$ D- i. a/ }
      %后向弧调整
9 H5 s# M5 u3 Q' n( L       f(t,s(t)) = f(t,s(t)) - dvt;/ P% [3 u- q! ?& F
end ) v9 c! G* T( t$ u1 T3 B& y. a
%当t的标号为vs时，终止调整过程
- f) |- N$ u& G! c( f, j( B if s(t) == 1
. e; T' B! [# A! U2 R) ~4 R. ~4 q& A       break;" l0 I# h4 X! i- J4 g  N1 u4 v. R6 y
end7 ~( {' Y' l8 w2 D4 N( t$ m
t = s(t);
- h% W7 w8 s) `0 J% e end
- a& b% }" R5 @8 J3 ~6 s %如果最大流量达到预定的流量值
2 I# V% P% v. |! B: \, w if pd
+ n) X  B- X( p' |( h$ l9 ^       break;
. q# h" I! }, q" a' P* O  f end5 N. g4 _. i% M6 i/ p  S
%计算最大流量3 j6 ^2 g# N2 |' F+ T
wf = 0;0 x' e# c* K1 I5 @( y- E
for j = 1:n
' e6 k3 t1 B# M* f' v" s       wf = wf + f(1,j);- m$ K' o; f. i( |# L
end: S& i- A8 ~* z2 D: A
end8 d5 U1 C8 W. |/ t
%计算最小费用* y$ ^1 o- Z3 _- S6 K
zwf = 0;
0 G) L2 I  f" ^6 o+ Zfor i = 1:n
6 t% d7 h" y& O for j = 1:n
) s, j) G" W% d       zwf = zwf + b(i,j)*f(i,j);% |$ V* j3 S. q
end5 Y0 z5 z5 x9 U
end
9 t1 y  B: Y4 c! n  q) W%最小费用最大流
7 t/ x/ o/ T. `6 j) O: Uf$ f: ^8 C" V) v) O) Z" O/ A
%最小费用最大流量8 N* Z* b. P) h: U3 t
wf" N# e. m9 y# r7 s, Z
%显示最小费用0 _  |8 d* ^# Y& G1 t
zwf
. @; `* T$ B8 x: H  c11 p3 w  f# R. I) y  e" m2 j! T  n
2. y( T# P5 K4 w3 |6 N* f7 A/ U) m
3$ x- ~1 H6 S# R! K9 s! x3 d
45 d, h/ F% k- _' Q- |7 S
5
# V5 R. v8 g8 M0 M1 ]) x6
+ p0 f- K& e: |% D* f8 E7
5 s- B  Q4 x  L/ R' Y  J8
6 G( M- ]" \6 R9
7 I5 d4 s& l$ p% I* ]/ T: U108 H, |1 x* q( Y% G& [
11
( B/ u, t2 y+ A# S0 ]12
8 r4 s% N4 C/ u. ~! i: _/ f  m13/ r  g3 h% S, }7 f4 i
14
9 t( l6 J  P$ a# S15
' y7 F# O8 A3 F8 P2 ]16
4 V9 S4 j$ }2 J  b; e. W17
, }/ _3 f- n& Z+ p4 M! e18
* e! [9 k6 H; z! S19+ D3 P% D6 c9 o$ q( t- {0 D1 B
205 `+ E/ ?2 M2 r! L) P4 j
211 t/ y/ F$ i5 J. Q, ^
22
+ n+ y, }7 R/ c23
& g$ w1 j  q( g4 M& R' F2 C$ M24
1 e3 c- \" R/ N1 ^2 e7 M% x25
$ Q. f, `, I$ o26
3 t6 @. A! A- f27: Z7 O& n! ?- N& n. T
28
: Y& a9 B* S) g' G. ~6 D9 J% l29/ L" d: I5 M0 b: D6 x; j: w
30
8 Q$ v% h6 a- D) z3 y( R* |" Y318 M; T% p& s+ e3 E/ i5 E$ l
32' v' Z) a4 q3 A2 U' [9 b+ R) r
33! K9 G5 ]5 p- E2 P  Y" i
34
0 g- L# ?6 Z' H0 t35
+ _% h. d; g, m* T2 z" g% N36
, i; Q' F% t4 `$ {, @0 L37' U) u9 h$ H9 u
38
5 c+ m- b' R0 o! {- X' z39/ K; @- D1 }9 R5 w& W! B! q8 G
408 A* b# @' ]/ A: o" ~- E. u8 y
419 e' F7 g) X! \: c8 B% y
42/ r1 L" u6 u# O3 V' o
43
2 }+ g* F( ~0 g$ M449 q8 Y8 [8 S. p- x
45
1 m" o/ w1 I1 y' o# f46. y2 G. t/ {, W& o% E
47( x# }% T) Y# G
48
( v7 h* a! [/ ?# B! j, s49
1 X8 F( f9 m2 q. E8 M509 r  ^4 N- [) f. |% G9 `* Z+ r4 l
511 s! e$ T1 g0 x8 @
527 r0 E9 l2 r/ m2 O* ^( m
53
* [1 s- t# u9 H: c/ e% L545 M7 f; Y- r8 z8 N) A; R  H
558 v; n# T% \: l2 q2 A6 M
56) {8 o: m; _! D- o% c0 w+ n
57
; W$ M# ]& @0 B2 A5 l" U58
$ F' G8 j" E7 B1 v1 S+ ^; y591 X, i4 F  P7 H! |
60) o5 p/ e+ E% m1 D' o, h
61: s% N7 U! o! h: J$ R  B0 i. j% a
627 a; C  r6 F: v: k
63! K4 G! w7 M6 X5 A) [: I$ @' x2 W
64
; X1 v& @$ O) Z" V  c65
' C- M2 \2 L1 H" D66: T' m, g# i6 [- B! I+ r
67
( o. `* b% c: u- D2 |  v68
8 E* x3 T7 w" r5 U$ N8 |( o- I69% T9 W) _/ j! K: J
70
. @& y3 F6 j& N' L# O71# [6 h) y! Q. |, q
72+ H  M$ j# `/ d$ E
73% b# c' Y# t$ @6 @9 l: a6 W
74
# M& [' _7 W2 t. s  O. v8 b, U75% P0 ^$ i+ x2 G8 k
761 ^$ _. w* F7 d: [# h
77
0 v, Z( h: J- V- _5 D7 R7 {78
2 m% B# j2 b' B79* D8 k. Q) t3 u" T7 }
809 Y7 }; ?& O# N4 p+ T! v
81. {+ {# _  `/ ?  U# [4 @
82
6 [% D' p( B, C! ]0 C8 s5 \83
: d2 I# C  p- b0 `84
( D0 M' V6 K2 S; ]' a% c85
9 O1 N$ t5 J, E0 `/ |3 ^866 n) m. _( z  i0 t% l: I
87- C; n  C( m+ a- ]9 I4 b+ p1 M
88
& t9 m$ Y8 }  p6 i89
- s+ ]. m; p# w* U7 ^% `4 b" l90
  D4 J7 L7 X- n91( Y" G9 R2 ?9 y& z7 v
92
2 [% _; v. Y. R3 _& o( \93: C8 u; q: N0 K/ K" B
94
/ k2 K' C) e: g5 F( ?9 y3 w95
- m% t8 h4 v5 c; `! ]6 H96
; O1 j- l$ s2 b97% Q, S- n! P/ c* Y6 F0 }8 M& K: E5 U
98
9 E/ H. l- t7 B$ Q99! D+ P- g" e6 R6 t# V& Z% `% N
100
( |0 w- W. N' J0 A8 {3 w) b101# p" _' b% e- @- k" G& X
102
$ x3 I9 o) a* e3 t2 j' Q1 |103
0 w: I( R( P, O+ D& p104" J* H" d) G6 ^$ r6 u2 x
1050 ]- A1 Q! Z: @
106
8 R' f  D1 [- `8 a3 X6 K107- O* \4 w* Y/ f- _. D
108
6 `! S; e. P* y7 ^: m- P) q5 I7 B1094 X( l4 F# P+ @2 x2 ]
110% p* z. s8 O/ E8 I0 e; Z0 r
111
4 ?3 |8 k* |4 e8 ^112
1 k6 h8 T2 `! D  U) N, ~5 G/ R1136 h$ A- @" H8 Y8 n2 V/ _6 {) Z% O. g
114
' y9 Q; y8 p/ g. k4 J3 v115* W; V4 z9 @0 m1 N0 ~( K
116' f9 A& |0 j: [/ r% n
117
4 V8 ~0 i/ a7 b: ~) c! F118) ]7 F. {2 E# O, J+ E
119
7 X: e$ o& R. u# y  g3 f120
9 `' b3 W8 e- V# F% v9 ~7 M121& b3 F4 Y" \& Y$ U7 q( o  s. y
122
% H3 D( P4 J+ G" S$ T2 b4 ?1232 g7 m1 x7 B% R- b/ ]/ Y3 y1 H
124
) B  i; P3 a; v; M/ Y1259 O- T4 d8 f- y5 R5 L5 q
126% J3 R6 Q! j0 F! \; \( F; A& U) J2 O
127, ~  O! D% g2 n. R, m, w
128
, L0 }* g! t$ \" K  g# _; t129' G! \+ E# Y. g' ?7 s% m( E
130
) x- h; T( E) U0 w$ ~6 l6 m: @/ @131  w! f2 [' C6 W) o- k6 Z
1329 T9 _% _% R1 S  R  {- h$ \
133+ B$ Y3 V' h4 Q$ f/ F+ m# d8 L! d
134, d2 C8 x2 H2 Q3 j# l
1352 t% t, f! |0 z' t" G7 ?; R; b6 q
1364 e# ~& [  y) l+ @8 u2 s4 M7 u
1370 [# [/ U2 F7 g( b- s3 ^6 F: D4 O! e$ n8 E
1386 X. f  Q+ l0 D2 e' R5 ~7 S1 G0 x
139
3 R: @2 G# D5 y. @, i1 N140. E# ~6 X* J1 ~0 v
匹配问题：% C3 S6 K) z6 O% _
1 e$ e; I7 h0 {. \4 q, J
最大匹配：- I2 C9 p6 ?# i) d" @, Q8 w

9 _" [- s' J/ V$ k0 v
代码实现：

m = 5;
n = 5;
A = [0 1 1 0 0;1 1 0 1 1;0 1 1 0 0;0 1 1 0 0;0 0 0 1 1];
M(m,n) = 0;
for i = 1:m
for j = 1:n
%求初始匹配M
      if A(i,j)
         M(i,j) = 1;
         break;
      end
end
%仅含一条边的初始匹配M
if M(i,j)
      break;
end
end

while (1)
%记录X中点的标号和标记
for i = 1:m
      x(i) = 0;
end
%记录Y中点的标号和标记
for i = 1:n
      y(i) = 0;
end
%寻找X中M的所有非饱和点
for i = 1:m
      pd = 1;
      for j = 1:n
         if M(i,j)
            pd = 0;
         end
      end
      if pd
         x(i) = -n-1;
      end
end
pd = 0;
while(1)
      xi = 0;
      for i = 1:m
         if x(i) < 0
            xi = i;
            break;
            end
         end
         if(xi == 0)
            pd = 1;
            break;
         end
         x(xi) = x(xi)*(-1);
         k = 1;
         for j = 1:n
            if A(xi,j)&&y(j)==0
                  y(j) = xi;
                  yy(k) = j;
                  k = k + 1;
            end
         end
         if k > 1
            k = k - 1;
            for j = 1:k
                  pdd = 1;
                  for i = 1:m
                     if M(i,yy(j))
                        x(i) = -yy(j);
                        pdd = 0;
                        break;
                     end
                  end
                  if pdd
                     break;
                  end
            end
            if pdd
                  k = 1;
                  j = yy(j);
                  while(1)
                     P(k,2) = j;
                     P(k,1) = y(j);
                     j = abs(x(y(j)));
                     if j == n+1
                        break;
                     end
                     k = k+1;
                  end
                  for i = 1:k
                     if M(P(i,1),P(i,2))
                        M(P(i,1),P(i,2)) = 0;
                        else
                              M(P(i,1),P(i,2)) = 1;
                        end
                     end
                     break;
                  end
            end
         end
         if pd
            break;
         end
      end
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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
最佳匹配

代码实现：

n = 4;
A = [4 5 5 1;2 2 4 6;4 2 3 3;5 0 2 1];
M(n,n) = 0;
for i = 1:n
L(i,1) = 0;
L(i,2) = 0;
end
%初始化可行点标记L
for i = 1:n
for j = 1:n
      if L(i,1) < A(i,j)
         L(i,1) = A(i,j);
      end
end
end
%生成子图Gl
for i = 1:n
for j = 1:n
      if L(i,1) + L(j,2) == A(i,j)
         Gl(i,j) = 1;
         else
            Gl(i,j) = 0;
      end
end
end
%获得仅含Gl的一条边的初始匹配M
ii = 0;
jj = 0;
for i = 1:n
for j = 1:n
      if Gl(i,j)
         ii = i;
         jj = j;
         break;
      end
end
if(ii)
      break;
end
end
M(ii,jj) = 1;
for i = 1:n
S(i) = 0;
T(i) = 0;
NIS(i) = 0;
end

while (1)
for i = 1:n
      k = 1;
      for j = 1:n
         if M(i,j)
            k = 0;
            break;
         end
      end
      if k
         break;
      end
end
%获得最佳匹配M,算法终止
if k == 0
      break;
end

%S = {xi}
S(1) = i;
jss = 1;
jst = 0;
while(1)
jsn = 0;
%选择NL的值
for i = 1:jss
      for j = 1:n
         if Gl(S(i),j)
            jsn = jsn + 1;
            NIS(jsn) = j;
            for k = 1:jsn-1
                  if NIS(k) == j
                     jsn = jsn - 1;
                  end
            end
         end
      end
end
%判断NL(S) = T ?
if jsn == jst
      pd = 1;
      for j = 1:jsn
         if NIS(j) ~= T(j)
            pd = 0;
            break;
         end
      end
end
%如果NL(S) = T  计算al的值
if (jsn == jst) && pd
      al = Inf;
      for i = 1:jss
         for j = 1:n
            pd = 1;
            for k = 1:jst
                  if T(k) == j
                     pd = 0;
                     break;
                  end
            end
            if pd && (al > L(S(i),1) + L(j,2) - A(S(i),j))
                  al = L(S(i),1) + L(j,2) - A(S(i),j);
            end
         end
      end
      %调整可行点标记
      for i = 1:jss
         L(S(i),1) = L(S(i),1) - al;
      end
      %调整可行点标记
      for j = 1:jst
         L(T(j),2) = L(T(j),2) + al;
      end
      %生成子图Gl
      for i = 1:n
         for j = 1:n
            if L(i,1) + L(j,2) == A(i,j)
                  Gl(i,j) = 1;
                  else
                     Gl(i,j) = 0;
            end
                  M(i,j) = 0;
                  k = 0;
         end
      end
      %获得仅含Gl的一条边的初始匹配M
      ii = 0;
      jj = 0;
      for i = 1:n
         for j = 1:n
            if Gl(i,j)
                  ii = i;
                  jj = j;
                  break;
            end
         end
         if(ii)
            break;
         end
      end
      M(ii,jj) = 1;
      break;
      else
         for j = 1:jsn
            pd = 1;
            for k = 1:jst
                  if T(k) == NIS(j)
                     pd =0
                     break;
                  end
            end
            if pd
                  jj = j;
                  break;
            end
         end
         %判断y是否是M的饱和点
         pd = 0;
         for i = 1:n
            if M(i,NIS(jj))
                  pd = 1;
                  ii = i;
                  break;
            end
         end
         if pd
            jss = jss + 1;
            S(jss) = ii;
            jst = jst + 1;
            T(jst) = NIS(jj);
            else
                  for k = 1:jst
                     M(S(k),T(k)) = 1;
                     M(S(k+1),T(k)) = 0;
                  end
                  if jst == 0
                     k = 0;
                  end
                  M(S(k+1),NIS(jj)) = 1;
                  break;
            end
         end
      end
end
MaxZjpp = 0;
for i = 1:n
      for j = 1:n
         if M(i,j)
            MaxZjpp = MaxZjpp + A(i,j);
         end
      end
end
M
MaxZjpp

欢迎光临数学建模社区-数学中国 (http://www.madio.net/)