查看: 2514|回复: 0

数学建模--图与网络（2）

字体大小: 正常放大

杨利霞

5273 主题	82 听众	17万积分

TA的每日心情

	开心 2021-8-11 17:59

签到天数: 17 天

[LV.4]偶尔看看III

网络挑战赛参赛者

自我介绍: 本人女，毕业于内蒙古科技大学，担任文职专业，毕业专业英语。

群组: 2018美赛大象算法课程

群组: 2018美赛护航培训课程

群组: 2019年数学中国站长建

群组: 2019年数据分析师课程

群组: 2018年大象老师国赛优

电梯直达

1^#

发表于 2018-10-30 11:25 |只看该作者 |正序浏览

|招呼Ta 关注Ta

最大流：

注意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代码：

最小费用求解) E8 _; N* G( [" n/ M8 V: c
# w2 f* h' Y  s$ T0 i
Lingo:3 s# Y8 p$ O, \

& R5 y% \2 J! p9 v6 Cmodel:
- G$ E! q# w0 ?9 _5 jsets:" W9 b  y) c! }7 t1 N7 e+ ~
nodes/s,1,2,3,t/:d;
0 ~( J* f. i; _8 @% larcs(nodes,nodes)/s 1,s 2,1 3,1 t,2 1,2 3,3 t/:b,c,f;
5 Q5 e3 }' R9 i$ d7 w# G( `endsets
; C9 B; W; g1 z2 `- @data:% d6 E: l! w9 u2 p8 Z
b = 4 1 6 1 2 3 2;: ^. g9 ^; I' h0 I# h5 H+ l
c = 10 8 2 7 5 10 4;3 |' N  {" F0 y( h7 ]0 W
d = 11 0 0 0 -11;
" P8 X8 Y8 W5 u  I$ z+ fenddata  E0 a% z# `$ f7 r: I* F8 h7 k, n$ |3 ^5 `
n = @size(nodes);! x- l' ?/ v! o; j5 j$ ~; J2 a
min = @sum(arcs:b*f);# P' E  z6 y( e6 f  [
@for(nodes(i)sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i)) = d(i));
7 E6 d' n% h1 j6 w( N7 l- W. _2 }@for(arcsbnd(0,f,c));
$ ~% _* h9 z' g; S, |end
6 p* w( H7 d0 T4 S( u& O) M1# ^6 p  a& l( p! F. x
2
) A* w& s1 \7 q/ [- W3
% @# I3 o& @1 G% ^4. D9 w% y  O, N9 W1 ~
5
9 q* G  J3 G7 o/ r6
4 ~4 y1 K" Z$ Z7
4 y6 r+ r2 r3 P) _- j9 W8: p6 r6 \, ?  u1 J9 n
9
1 g. n2 N# r, l* J% i10% m8 R% X' X& |
11
! W9 Y3 d; ^. W2 h12/ Y: L$ h; |, {1 w6 N
13+ m2 L: `  C9 s" E
14
) a$ i1 B% O, m- U' q% Y, \15" E9 P+ f+ L; E/ Z+ \8 w+ a
Matlab实现：: c: ~) }+ T, \' v# W

& l. B+ P- T7 n- c  Z% ~
8 V) `/ O4 z+ c* h( k( tn = 5;4 [  I  G) w1 k6 C0 o
%弧容量
# X0 x; J7 I5 h) Ua = zeros(5);
$ g" `  E! L& ^4 C0 ~0 F& D* na(1,[2 3]) = [10 8];
; W5 N) o4 X" u7 sa(2,[4,5]) = [2 7];% r) {" n$ q8 t" S: Y
a(3,[2 4]) = [5 10];
8 m3 J; O% \4 e  Ma(4,5) = 4;
/ P; F5 R( Y% ]" s  a4 X) _C = a;/ z+ }! \- B# ^8 v. E/ _
%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];
, @+ z* L7 O" A7 ]%弧上单元的费用5 B& A* F- |; R# V& k' l
a(1,[2 3]) = [4 1];8 _  n  l2 F' q- v" O, R
a(2,[4,5]) = [6 1];
* ]  E+ P; n! J7 [; Y4 ra(3,[2 4]) = [2 3];
6 v5 I5 |$ u  Q5 T0 [a(4,5) = 2;
8 t: F1 t2 ]" B0 F# ?% V# x7 p1 cb = a;5 F$ {0 n$ \2 {) Z8 S. c
%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];  j* f9 ~. ^  K% `4 E/ p4 ~
%wf表示最大流量。wf0表示预定的流量值
5 K4 P7 ~" W5 d$ q4 ?" |4 [wf = 0;5 F$ I$ ^7 {5 I: h/ w; K9 U
wf0 = Inf;
4 w# {( |. i; r3 j%取初始可行流f为零流. f% S) v$ K; ~
for i = 1:n  }- _- x% E* d( R, o( R
for j = 1:n
: |4 R1 _/ h9 E6 [1 s9 R       f(i,j) = 0;
1 J8 q& \, R! k9 {+ A! l% `# ? end! _+ g' E+ R4 l9 C8 {' v" e
end6 S' e2 H# ^6 |$ @  v/ M0 o
while (1)
8 P  ^7 B- h3 k/ h4 s( [! ]# ? %构造有向赋权图5 {9 `# t3 O( A+ T/ j9 |
for i = 1:n
+ a  I  f: k& j! D5 L) Z       for j = 1:n) ?" C4 e/ E% f# f, ^( R
         if j~=i
$ U0 H& W( j0 k6 t             a(i,j) = Inf;
3 _5 B- `! O$ h/ z% \          end
7 n' W0 \+ _; Y6 L. s, E4 a- K+ O       end
3 h, C6 N/ S( w0 M end9 N7 L! D  C/ D  }0 r
for i = 1:n; _) ~3 x" K' Q- R' Z
      for j = 1:n
8 ?4 p7 x, x- J: X          if C(i,j) > 0 && f(i,j) == 0
0 ~+ a) h7 y5 m; i  d" k             a(i,j) = b(i,j);
: }- d' P( }4 {# Y          elseif C(i,j) > 0 && f(i,j) == C(i,j)
) i0 S9 M, G" h. W. @7 A             a(j,i) = -b(i,j);
/ H; r/ U, T, C9 G" [          elseif C(i,j) > 0- D9 z! W4 B' [6 \) n1 z
            a(i,j) = b(i,j);+ v  k& G5 `( }$ J/ m4 d/ G
            a(j,i) = -b(i,j);
3 S9 k1 E0 p4 u7 g          end
1 ^, K  K0 R5 f* C, u$ V       end* c- H4 A, w& r* M" O
end3 v! v4 t' y* A5 Y9 G1 u4 `
%使用ford算法求最短路，赋初值
: M2 p! a1 R/ i( c7 E; _ for i = 2:n
7 ~. G: m" C6 G+ ~       p(i) = Inf;, ~9 Z1 E0 E9 Y8 T7 s: o0 I/ v
      s(i) = i;* i. t5 W- c7 Z7 U" Y+ V$ `
end1 S! u* |- }+ F) E5 H+ G
%求有向赋权图vs到vt的最短路,赋初值
0 ?8 Y+ [7 _5 x# k4 Q& I1 L for k = 1:n
* w( l( x; t1 n& X' L       pd = 1;
% m1 U6 ^0 H, h       for i = 2:n; ]# e; h! P9 y4 x$ T8 P
         for j = 1:n
% E7 G* m& R: K. H: q7 N, X# c* n             if p(i) > p(j) + a(j,i)( b$ T2 O4 p" [% h! m  X
                  p(i) = p(j) + a(j,i);3 u5 {8 n6 ~% h) m3 R
                  s(i) = j;
3 H  i) n2 \5 ~# O                   pd = 0;9 F/ Y3 c/ l: p6 a
            end
. n5 ?& M, L' }! x! q/ B8 }          end
% \' z0 z% T, ]: P5 D! `; @5 g/ ?- Q6 E       end
8 V8 \. x2 S+ [% D0 m- k) R& m       %求最短路的Ford算法结束
7 p% A4 h+ V+ z7 W( m8 g9 Z: B       if pd5 X! B; s! p, v9 m/ M# f# F
         break;. d5 c2 ^! D3 g, ^+ K) U
      end% S0 d# R0 {" b: _; c
end
& I3 i+ t) s) ` %不存在vs到vt的最短路，算法终止，注意在求最小费用最大流时构造有向赋权图中不会含有负权回路，不会出现k = n
- M: d% x* B" s$ K# A) ] if p(n) == Inf
, Y1 ^' q- H; n0 g/ w       break;
- Z) f) X: ~( z end
8 w/ @! z% z2 M2 {, r# l %进入调整过程,dvt表示调整量' Z3 B& b) [+ \2 [  l' H
dvt = Inf;
( x  e( @0 W4 F dvtt = Inf;
8 S2 v/ v. }, S5 L! g9 Z t = n;
: h- e1 C% l+ s( A while(1)7 Q: q/ Y/ W5 a3 p0 n
      if a(s(t),t) > 0
3 [0 l( l5 M" p2 f          %前向弧调整量
8 _9 Z: h3 ]5 c$ N4 h. x0 u: f          dvtt = C(s(t),t)-f(s(t),t);
8 ~) ^% G0 I3 C5 J# a. I0 V          %后向弧调整量, _  ^# f' B4 q  P) ~1 S
      elseif a(s(t),t) < 03 A* P' Q% N5 d" i
         dvtt = f(t,s(t));6 u2 a9 W, n, [0 w  G
      end
$ x: z9 E- `2 ?6 E# v+ p* ^) E       if dvt > dvtt+ P- y; \$ M3 N0 K8 ~9 d! l7 v
         dvt = dvtt;
& m. V3 u* @3 ]' M& t; o       end
+ j4 u- Q4 l) p* X( q4 _       %当t的标号为vs时，终止计算调整量
5 M. f0 w/ X& S/ h. `+ ?: f: }       if s(t) == 1
6 f4 n( Y- s9 \          break;
% V8 b- n, y: R       end
* A; ?% m+ E4 I8 p; U7 l       %继续调整前一段弧上的流f
; A/ Q+ s5 n" B9 N       t = s(t);% F, b/ F9 o. E8 m% \. T. Q9 Z
end
3 k- g6 l8 K/ a7 d+ ^8 ?) H; j pd = 0;! C) X% {& ?3 y
%如果最大流量大于或等于预定的流量值& ?( e. [0 x8 K( T! h- S& F5 {
if wf + dvt >= wf0
! l5 n& L7 u4 U$ M* P; o       dvt = wf0 - wf;* ^' I. ?6 B  O7 j, H
      pd = 1;. _6 h8 f1 n3 V+ d
end6 U7 Q& j2 ~! x2 z  q0 D
t = n;' R! s$ g, w; i; A5 F# b1 ?- A( K
%调整过程
3 \6 D* O- Y1 J% ? while(1)& l3 f0 X( I! |" c( K
if a(s(t),t) > 0/ d% a/ C- T7 ?# k! y# k9 i
      %前向弧调整
% X0 Y' K9 D3 x5 A       f(s(t),t) = f(s(t),t) + dvt;
3 s( F1 j: }( f! y" W elseif a(s(t),t) < 0
  A5 P: o) G* L) z/ v; L1 Q       %后向弧调整7 c; F5 z7 A* y2 I
      f(t,s(t)) = f(t,s(t)) - dvt;3 o- v: v4 l1 h
end 5 T, ]+ t* c3 V& Q
%当t的标号为vs时，终止调整过程, Q) X% I3 A4 T1 g  V4 F0 D
if s(t) == 1
7 U$ i7 [& P* `6 ~       break;2 t1 i( q* ~9 ]5 W( C) `: L$ q  n8 P+ N
end
. {5 w: W9 M$ }& L5 [ t = s(t);
, a, u7 F4 o) M! q- k1 e end
0 ?7 T+ X8 m: W9 I %如果最大流量达到预定的流量值. ^0 w6 w5 T+ {$ x# u) \% e# v% L
if pd8 d$ G9 Q3 F5 X1 n
      break;
4 n( m+ g& r3 c" @1 ]. f9 @ end, Q9 E3 f; E4 p! _
%计算最大流量
4 D, M% ^( |7 d2 F( Q% i9 E! t/ M wf = 0;1 p1 G" G# n& e4 h! \7 b
for j = 1:n! b6 X0 C) I; M$ P3 N& g( i6 B. s3 j
      wf = wf + f(1,j);
1 w8 `+ _! @5 ^4 K) `" ]9 o- K' C end
1 q9 x% Z' A* ]0 Send" e9 C% v" m; k' Z9 P
%计算最小费用
7 B6 o3 S. z' @+ v& gzwf = 0;7 M5 S5 i! `$ m! s( S* i+ q8 X
for i = 1:n/ j3 Q- o' z9 ^' @8 F
for j = 1:n6 ?( G' y# J& x: x6 j+ U- l5 X
      zwf = zwf + b(i,j)*f(i,j);
9 [) }- e7 E" u3 ~* i end
3 \$ p: a6 I+ f7 D" i' ]; v. hend
' \5 E' P# U8 e7 x: ?0 {6 i4 ?%最小费用最大流
: q) m4 J' w2 k& |+ k: Qf! w/ J1 p) x/ ^! Q. _+ }
%最小费用最大流量
5 y$ K" e  A) I2 d. P! rwf
6 D' c3 l3 J. S%显示最小费用+ v# l+ y" B- H4 f& J
zwf
' e& q5 @$ F& U8 [1
, u& u2 t# b2 o0 L+ j2
& s0 d" m: b) n6 Y3
, O+ q# C; k9 A3 R' U; d  Z) x4
2 v; {" E! K  z- G# I6 l56 k4 @2 ^- t1 g- y
69 \) G& t  ?) x9 S% ]7 Q/ x) x  s
7
- {0 ]5 Y3 S( b# w/ f- O! J8
6 \" P6 O: o' F* z5 Q9' H) L- e+ J8 d- b- ^4 t( M# ^) {
10
. B( L7 v* d" K7 O; [+ X11" V, i$ c3 x4 a+ y" @; S; ?5 K- V
12" ^( c0 M4 r4 g7 m' L. O
13
# N4 I- n. O2 i9 |' ^5 l. j14
  _( o0 {8 ]9 |1 M15& s/ i7 a4 m; p  B9 ^3 i& |* Q
16' K! l# ?* R' _* s4 a$ ]& c1 K. {
17/ u- y; V" B6 _! Y! K: g
18
, `7 i  Y3 r! D& t5 D' p1 p5 A' F19
: [% C( |/ d+ L/ B$ s' Y7 q20' d6 a: r, z& G3 a" w8 E5 H" C6 [$ N
21  q# ^7 K) |1 {5 `- g/ L
22
1 ^. Q+ l9 v3 J23
5 u5 f& |5 ?$ C# u1 ?; I24
( S' d& {" O( `+ g& A25
! ]  c. e& T3 g  h9 M8 w# }26
4 q/ v8 {" v1 Y5 I# Y$ C9 g0 M27, T4 b+ B& z$ g7 E
280 I4 \! `% {, C- X0 }" T% ]0 c$ V: ?
297 n4 l; V! N$ I
30
9 Y  {9 x& |, l% {% `6 f31' _  ?6 v  n% X" I
32
. o) M" e3 }# J6 R! P" K33
1 A( w1 v% [' p7 I& d$ Z( T2 @34% `+ j+ T+ j3 _: ]
35
& K% k/ e4 y, S1 f3 x+ I; d4 I$ s36( t  g+ E9 `) l+ K/ s( N( T
37: b7 z$ a- R% z& P' b( X
38
3 J# J& d+ w" X1 y1 e39* V1 ^, O+ u8 V' d: Y
40
41
& K" G! Q0 E' i! l, q5 b42$ N! u! \4 z" p" D4 P# h8 M8 y
43
: o. I# u5 {" i! w& _! U3 h! f3 P2 X44
0 i4 i/ t; s: B& O$ s5 T45
5 ~% g$ w( O4 _  {8 M2 p46& R7 ^0 a( Q0 g
47
" Q  r1 y) m: L1 s. J48
1 K# ^0 U. P3 @8 @* V49
: }; e4 U. X4 T# _" ^- v. P50
' L+ i2 |2 Q2 j% T& C# q* t% M' m5 n, z515 v! q0 U3 t$ r( b& }
52
& X1 N% K2 @' Z53
( m) v& n/ d, y: F1 v54' d* V7 _+ n  h  S3 |
55( N0 w1 g; n, u# p
56
. B( m. Z- U2 ?6 s5 n$ ^57
& @, C9 [: ~. h. I" m58
5 S' h. }7 J7 b; [- L$ M59( V, ~' W. ^! G$ a
60, v6 E6 _4 g9 e9 _# s0 R
61: }0 S2 p4 l: N5 V/ p$ y3 C
621 A5 f# L. J! b. x! M# p. d
63
- R; D( M2 B" U! _% c64( Z( @) V% e/ U2 {% [
65
665 I" T+ m) n8 V1 q- n* j! W
67
, `! [6 S4 p! a! c68, y+ D: V# q+ M, c
69& L% S2 c1 l6 ~! F
705 x9 I8 l- y. V4 U4 \: L
71. n. L1 ?4 q# @. |& i
72
9 T9 b& {) g4 c9 \73
& B6 z- p" s2 q; |! ^744 j& b5 W1 I  j
75
+ P, d8 b1 Z( D4 P( n  W/ T762 R0 p, _) e% H+ j* f" A2 @9 \
77; M  C" Y) {5 M* U. D
789 P$ A# B* O( X9 E- z
79
* R; M1 ]: y  b9 }2 `% o7 ?80
9 v4 ]0 R' E3 _+ s- R81
) j# m9 H3 F* N2 R, p9 f5 ~82
. ?, f! O" A8 F9 g6 V3 B& s) |! m3 @831 p! n+ m, a+ H3 Z; U7 J, o
84. Y) W2 {: w( [% `7 j# s5 Q0 p; h
85' p8 e  {* S; v, c$ ^* q: `& ~
86
0 N$ R+ F6 D1 ]6 T2 L8 f/ U87
5 I1 f2 A. {  D" U! Q3 K' T8 e88
1 ^# E/ P" P! J8 A0 h3 x. b+ [5 S- h89
# M& R# D5 {) B8 w" E, c! R90
, h: q$ L* d: K6 I3 G/ ~91% Y( S) C, |) ^# u* t! i
92
+ g, a6 R4 L2 i# G8 F! l93- h) z5 h) v9 ?
94
2 }% \* j2 l+ o6 |. i95
5 P8 n4 w" ~' k6 w5 D& {& h964 ]% _) m; q* b# ], J4 A5 b
97, H/ ~$ [% E% x3 T$ K, d
984 b) w$ l, ]  B3 O5 C9 F7 _% c
99
: |. ?' E2 ]+ I6 T' `100
! O, h, L+ B6 Y8 d5 T101% v. Z% a: G. c8 F$ h
102
& ], W. `7 w! D. H, `" H103
; H! ~5 l0 i- u: q& j. o( R; Y( S104
" n  e1 n& v9 J/ B105
  U& b/ k: b  W% R7 Z106
0 n6 s5 Q5 j) l2 m) q1 `* w( O107. u) u. u1 d; t: K$ O
108
/ Z; R; L: }- M$ l1091 z- V: b( h; O0 S6 l
110% ~" W  E3 @. E' }: j
111
( L" v3 l4 l3 U# s5 N" X+ G1 S112; T9 t' n8 i! ~* x7 ?
113
7 ], {- [) a# D4 d2 I114! C. U% @, F  j& Z" v% Q6 b4 ^3 ]: i% C
115
! h2 i9 N5 Q! ]2 V116: w0 E& i+ G. H% X
117
3 `1 |, X# B) N) E8 `118
) q$ a/ C: s* G/ n) H119: R% M8 j8 [3 ?, |2 Y
120
  G- A1 f! V0 }) ^* A121
7 q4 ]- ^) t/ p/ o122
1 \% x+ R9 ]0 g8 D( J6 f123* X5 G4 C: B' h  c7 n" G! }' h
1247 h9 v' x5 Q* M+ ]2 K
125* u7 @/ J8 A2 j3 J% q
126" o" s) X2 j$ Z* n! N' z# F
127# T! G  G% N/ }1 c' E# Z9 d
128
+ k7 x3 J8 x( G7 R3 \0 t. W129
) V/ {0 U* \- o6 f130
0 G6 a/ P. ?  j9 V" M* h131
* ?7 A1 e* H) ~; o9 ?. w1325 {7 v. \6 r4 |5 g- R$ k
133
; ?  _) ~2 Z2 f, u: D134# M) w4 p( V) Q& h
1351 R( G2 _2 B$ s' e, T) P0 N- @
136
0 N7 T$ T' ]- B8 Q8 o# E9 G1371 A+ E% s! Q* E2 U/ F
138
# a/ T) q1 \) v139: o& }, m6 j0 p6 `0 @( M
140
- c4 w* m, N9 B( K5 W* A匹配问题：
% u* p4 ?0 w, }0 N5 Y, f, J5 |7 m7 z2 {# N# n* h2 Z& D
最大匹配：
9 m) x! C" Y" Z0 m6 q/ K3 ]5 g# N6 G: q5 B

代码实现：

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

zan