查看: 2515|回复: 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代码：

最小费用求解
; I# q3 [7 y3 L3 S* ]  s& ?) m. x& e+ m5 G" f7 n
Lingo:1 {5 j9 H7 a6 S; f
: J7 v1 A" |; b
model:
3 I/ j* M* u+ ysets:
nodes/s,1,2,3,t/:d;
3 @2 x6 {, E1 B' ^7 Larcs(nodes,nodes)/s 1,s 2,1 3,1 t,2 1,2 3,3 t/:b,c,f;
! s6 ]4 D+ q$ s5 g6 S+ Yendsets4 f( j* ]$ d& G
data:
& }$ V2 v9 m' P$ `b = 4 1 6 1 2 3 2;2 ^9 T9 H. u, B' E# b
c = 10 8 2 7 5 10 4;
- j9 C$ f+ g  n. nd = 11 0 0 0 -11;" V+ X5 U+ Z3 w: g' l
enddata
! |/ t2 H4 I  Vn = @size(nodes);
4 T6 i$ F- {4 _0 P3 xmin = @sum(arcs:b*f);
" K5 h: T/ N. v6 l/ K6 T5 S9 l@for(nodes(i)sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i)) = d(i));
+ q& K' ~! \, R- p* g* N: Z@for(arcsbnd(0,f,c));
7 H* r6 s& k, `! x4 t, |' z2 ~end/ G9 |* w: q  s5 _* }) B7 D
1' F8 W7 Z# z. ]
2
3 _5 Y8 \- I; G( z( ]; n1 L3$ L: O& ~2 J1 y% {& f
4: b) Q6 j& V, O" t/ Y" j) k- h
5
2 z. p% ~' @7 v/ P' h6
6 U* x7 m& t( r7 h7
4 F: @8 c1 n; s: \8+ J, q% G& H' a. [8 q3 ~
9
  c# u: D& T# k5 Y, o$ |10' m7 X7 S( h5 [4 q" T; p5 n
11
! j8 y3 y1 w( Y' X12
8 [0 C2 T. |% }; n13: V' U. V9 C5 |" _" a
14
5 W: q# |$ ?8 S6 j# A3 `& v15
: S/ N/ R6 h* Z' c  l$ h  RMatlab实现：9 B3 h) A6 c% u7 G* c( w$ a

' a" q* s: {8 _/ X' ?) {+ h: A% ^' y& p, B3 f8 k3 R
n = 5;2 @3 f9 P" [+ f/ ?5 U9 Y
%弧容量
8 o8 p2 f- R- n6 z9 f; _a = zeros(5);9 }! {  y  ?# E0 L9 E
a(1,[2 3]) = [10 8];# w/ g& ^+ t8 B2 k- k
a(2,[4,5]) = [2 7];! L4 z' ]$ t" [6 f; P5 @0 `
a(3,[2 4]) = [5 10];: `! H: _6 ?* o
a(4,5) = 4;& D; t+ d1 O3 Q' q: x7 Y7 I0 a5 w3 t) C
C = a;2 `6 t0 B  o6 U! E. I1 Y7 N
%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];( j7 O) s2 I" h5 j
%弧上单元的费用' G3 h+ D. s6 q9 ^- u- V9 I7 B
a(1,[2 3]) = [4 1];
2 r! F1 v4 A! k7 A/ Z3 ia(2,[4,5]) = [6 1];! r. u  {9 J  P
a(3,[2 4]) = [2 3];
9 c8 z3 S9 K2 `* da(4,5) = 2;
3 A7 P3 h# k0 }/ b" A0 Sb = a;. ]0 M( N# K8 r' F: |
%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];" @8 ?4 n" T+ r& |( M
%wf表示最大流量。wf0表示预定的流量值, S& {0 `3 C! S5 k  T$ k9 k
wf = 0;4 s/ \3 T/ N) E8 n2 b/ B
wf0 = Inf;
* [: K0 `% K4 y+ L! K$ M7 ]7 |%取初始可行流f为零流& z: h( i: v- x1 }8 j$ x
for i = 1:n
: A6 ^  C9 ]# _6 F; l for j = 1:n; G$ d3 O* m5 o" r  R; i% z" z
      f(i,j) = 0;( W+ K% E6 t0 e3 Q$ l
end2 T" [+ e# e/ r6 @; ~- \. N/ V
end
3 ~8 M- `& o* s- N  _* K$ fwhile (1)
5 X2 g9 {) Z1 v; {: b9 P+ M %构造有向赋权图' z0 R# O. ^9 n- ^) p% v  {" ]4 O
for i = 1:n
. {& l( u* a" O1 x       for j = 1:n& i1 n* l# m* |* t) d
         if j~=i: H+ V- [% g6 g0 `& c$ ?
            a(i,j) = Inf;
) d8 H) k4 R4 ^4 [" j: K          end3 E$ d; _! U. ]: x
      end5 c6 ^& R6 p/ X) }- b. I8 f
end  k- z9 Y# T+ N- L2 x# k7 V
for i = 1:n% k3 W3 r* o# u' C+ @
      for j = 1:n! k1 s1 D- ]6 y5 z' U
         if C(i,j) > 0 && f(i,j) == 0
* a+ F/ `9 r9 {$ e3 m8 A7 Y             a(i,j) = b(i,j);" @+ e5 p. j3 s1 m
         elseif C(i,j) > 0 && f(i,j) == C(i,j)% P5 a9 k/ ?) P5 w" h% J
            a(j,i) = -b(i,j);
0 m/ O" W$ N  v, x' g          elseif C(i,j) > 0- d. J& h. b( X8 G, `/ O7 f& Z
            a(i,j) = b(i,j);
) B8 d4 J. L3 k" O. Q: U) {# N             a(j,i) = -b(i,j);5 T# a* F) m/ E7 W
         end  W/ |- i" L7 k+ ^8 e. ~) D
      end
# L' Q5 W! y+ Z. C: U end
1 E* Q# e' `! e$ s %使用ford算法求最短路，赋初值/ M  ]; d8 A" w- ~3 I/ Y. [; w
for i = 2:n* U+ _3 X: s+ G/ d3 o+ D% [+ e2 q$ N
      p(i) = Inf;
' m. v3 y! P% _: h. l       s(i) = i;1 V4 ?* B& O4 o' d3 R5 z
end  ~' U" n8 }8 s: H' i9 T
%求有向赋权图vs到vt的最短路,赋初值
3 o* z) I' |- @& Y for k = 1:n! T+ i0 @, S3 H0 t. G" p$ E
      pd = 1;/ D! q) p/ v+ y+ ^, k
      for i = 2:n4 x/ Y5 A% w( X2 A" ^
         for j = 1:n
& V" H( r) ]6 y0 m             if p(i) > p(j) + a(j,i)
% c$ A, p- u0 m6 D                   p(i) = p(j) + a(j,i);4 T$ X: F6 n; a" _* n4 \, _6 n
                  s(i) = j;6 ?- l4 a- Y8 m5 T3 i
                  pd = 0;
% ?/ M) z: \/ G             end) K& _( d- k  R* z( U+ u  Y  Y; O
         end+ R/ h+ W% @: n
      end- X9 B' j! X9 c  Q3 l
      %求最短路的Ford算法结束! p: ?2 ~* b+ k7 g1 t  h  k2 @' @( L
      if pd" ^6 \* L) u9 |& h9 x: Z. b* q
         break;# o3 o6 V3 G/ ~$ K" `$ v/ y
      end
" e1 b) e6 Y$ o4 n3 t" } end
$ W8 K5 d; c+ S! P %不存在vs到vt的最短路，算法终止，注意在求最小费用最大流时构造有向赋权图中不会含有负权回路，不会出现k = n
& _: K. S& k& ^7 R3 s- P0 { if p(n) == Inf
4 S1 M2 g; w3 T3 ^) V1 C       break;
8 w2 C1 Y; Z- p3 r- _/ D8 N( r& v end7 z' p( s8 l9 _: c" [( D+ {
%进入调整过程,dvt表示调整量5 ]' }) h. h& e' i- i4 f6 V
dvt = Inf;
# y8 o8 Z8 i0 c; [8 H dvtt = Inf;1 F4 I& {- h7 c0 T: T" q
t = n;& ]% T$ r: ?: x3 s
while(1)
: i5 L# ^2 X. g$ I       if a(s(t),t) > 07 N' D0 ?$ @5 X: T; Q
         %前向弧调整量  m+ o% t, a9 @
         dvtt = C(s(t),t)-f(s(t),t);
$ S1 n- Q' ?  d( L! p% W6 a0 L          %后向弧调整量
4 q) T7 e3 ~2 M/ i7 I7 t/ F       elseif a(s(t),t) < 00 {( A+ v2 S6 E$ h" U
         dvtt = f(t,s(t));- _; t- y! F- U/ z# \/ ~
      end/ s0 B- i8 X1 b$ L/ l
      if dvt > dvtt/ W# g% @4 E2 N) Y1 F
         dvt = dvtt;9 }  a( `; B9 C8 U) L# t
      end
; I: T* k5 R% H! p' G* @/ ^; L% d       %当t的标号为vs时，终止计算调整量
, D2 h% X& m5 T/ Y, V3 J7 z       if s(t) == 1
+ ]' x9 r$ R5 i% P          break;  ^( f# e* |% ^
      end5 E" P! J8 Z! D/ U  I
      %继续调整前一段弧上的流f3 T" `5 ^. w/ D: P* X8 }  ^" T+ m
      t = s(t);
: v. S, E. z; M- n end6 M: i( D7 }8 r1 ~# e" |* d
pd = 0;
: x' |0 |1 b* m$ k# a+ { %如果最大流量大于或等于预定的流量值7 ^- M  e( }4 c# F' p# T
if wf + dvt >= wf0/ g8 E- [5 E5 R* n
      dvt = wf0 - wf;, l$ `: \& t1 P7 `2 L. j
      pd = 1;& p: E9 w/ o; p0 B5 b6 ~
end
3 {- w# e( T( {7 E3 ~ t = n;
, t. r8 o  {* ~3 R5 D; q %调整过程
/ p) Z7 ^" a( F: B9 ` while(1)2 R' N! C# m9 _* `
if a(s(t),t) > 01 ]" ?' ?9 {+ o& D
      %前向弧调整% A+ ^! d$ T0 x: A! v
      f(s(t),t) = f(s(t),t) + dvt;
5 A2 c: ~2 H1 w5 P. O, R- U7 H  {. ^ elseif a(s(t),t) < 0; W( W7 P- |8 z0 J9 O
      %后向弧调整
) a3 I( ~- m! G4 w7 U7 C$ H$ E; y% G1 R       f(t,s(t)) = f(t,s(t)) - dvt;+ V/ M0 r) |9 A1 f9 ^
end
$ ]+ \2 c3 J# v) s9 v %当t的标号为vs时，终止调整过程
. ~7 U+ z7 g: |) d, v if s(t) == 1* Y/ c4 ^0 o# v6 x) e8 b
      break;% h6 L8 E' H) ?4 T; R
end
/ P, C2 c" z" N9 e- T t = s(t);
$ [, l  j3 ~8 X8 z end
3 S  [& }6 \4 f6 w$ q3 c9 d %如果最大流量达到预定的流量值
6 \2 U. {0 p; p$ ]% q* E if pd' H* [  V8 S- r6 d! M
      break;
. j/ P  g+ i  z. X end
/ h3 K+ J( i0 \% O %计算最大流量
6 T+ C2 S0 M. b: l wf = 0;2 B) t" B  K: D& w2 e- n# [6 |; \
for j = 1:n* c5 n3 i& g% B" K
      wf = wf + f(1,j);
  S8 T& N( s" E) m" r end* y7 X$ M* A- \: f
end
/ G5 h9 ?6 o3 U: u: {%计算最小费用' x: O3 ^. \9 _  T% A" [
zwf = 0;! b  O( R7 g" V/ g0 X3 t
for i = 1:n7 v9 v- A' C& g8 e; o: S0 {
for j = 1:n
0 d# \  s# c/ M7 H2 ~, w* V       zwf = zwf + b(i,j)*f(i,j);
- N7 @, q, t3 `8 _# O8 s" H end
0 K& ?% Q6 v( @6 G4 o3 jend0 ]4 A/ u5 z& X+ Z3 X$ Y% e9 X
%最小费用最大流
( @6 g) L/ k/ D" D: c$ g  D$ e8 tf
4 W' }% m( b+ T8 ^$ q%最小费用最大流量, C( p# n! X6 V8 {
wf4 G) H7 K1 ~9 F6 o9 }8 [% e
%显示最小费用! ?- ?2 G4 y+ [, y, b9 m8 Z" n
zwf
3 ]3 e% I6 w6 B& D1
4 n5 \. |3 s; D2 ]" o* N+ |2
+ T2 ~7 t$ i. e4 V3& \3 J: K( p( y, n- G& I/ O( M) e
4
, d6 `4 _% }- }% [5% @8 z6 J& u3 }% T3 o5 k
6
6 K# Q( L4 L$ X0 V6 l/ W3 a( W7/ h! k; X! }7 j  Q( {0 M6 x
8% E0 P  C- l6 h$ h* W. j
9
& I7 t" t/ E2 X) ~' \% ^105 `: s% C0 m% |( ]7 |
119 Q; J( e: H1 G4 _0 \0 |
12. R5 D' ^. Q# l: F& B& N! X
13# Z2 @( e5 Q4 ]: [( a+ E+ p1 y4 n
14
8 [! ~; S  B  ?* c15
! l# R$ ^. e. e+ a7 d3 Q3 L169 N6 d2 b" Z* j9 @( V" A  A  E$ w
17+ N# N1 }9 H" b7 B
18
. ?" H' X" R+ t* X+ J" q8 E3 c19
5 X4 Z* d. J4 h2 e; S3 u* n20
7 b7 h, N9 M6 f% s21
0 a7 a  p' J; B22$ o- \& z( k0 q6 l4 u% G" o6 S
23
/ \9 y) \$ W6 y' C; A- }24
* f+ k1 h# D0 L2 O8 o2 h& A25+ {' [8 Z: {6 p% x. h+ m
261 \+ I% ?  M0 P" P) i3 @
27
, N0 H8 N% O/ E# a4 [. E  z285 T5 j0 [# i" A
29
& X. n) `! d" H& O# Y% U30
% ^; j2 ]1 K* F% S" @* f31
/ ~. Z+ C" E# ?. l) n$ n325 [: g9 B6 m! r/ T7 x" s2 a: x
33
( k% |( O# x6 u. V9 C8 S0 ]34
; l# f4 k, x$ s, l& l; O  b3 _35+ g+ F! o" I  a1 j' _% q& m2 m
36( {' v( \4 c+ c0 g
37
: g3 w- e3 S- q; K# |38
. n+ _. y2 F+ w7 e7 b$ J6 a' N39
4 Z9 ~4 `6 b7 [' @2 F: `- e401 E2 [3 x" U! f; a4 }
411 d! b0 ?; O9 h$ s1 u  ]
42% [4 v8 B, }! x8 ?8 Z4 _
43
6 U4 q: n8 V  r1 p" I44
3 v! T! S( ^6 Y8 e5 ?1 i" D457 E! F1 X; R4 `% E8 t7 d. t
46% L0 Y' ]9 m( J- l" y
47  L7 `1 c$ k' c) C0 k& v- R
48
7 x# N; o5 y3 ]: s% r/ z4 c9 \5 a49
; u% |4 X) |) b" A: y' I% W50' i4 ^& e0 m, p1 z
51
4 F. ?7 x5 K! v( u% G52. T4 D- L5 i0 ]# X
53
5 n6 |/ }0 ^# f54
" @8 @: Y- l/ ~6 {% {55# x' ]/ n: \" ?  [$ ]
564 g& n/ P, j% _) W) {
57! v/ M2 |' D$ R& N; T: Z0 A9 a
58" n0 ~# ^" D6 E2 i7 ]+ T
59( ?! \5 N  N' v' o
60  l. b& y- ?; f4 x7 H, L4 N3 k; W
613 F. h* V7 H* u% m4 k- v' m5 O
625 _. `, v* _' E0 e) a+ `' ], o' l' H
63
2 B- ~# R, f" |2 h9 @+ E% u' n& ~64
. u3 O8 j" ~0 e- p/ A# b65
# I$ ?! x5 z, m0 E5 Y% {66
( H7 o# c, m8 J& t7 m% n6 _67
) I- C+ M/ ~( a8 O68
- ~( B8 X5 i3 p69" q6 i3 a# B3 i' S
70
/ y% z% F, G: u7 j1 u( S* D- a& d714 s& O7 d" v, Q& s
729 s8 g* y* q; G, E: G
730 T/ g  m- S* W2 p
74
* I' D4 d  t, [. a: {75/ S* U  z2 I) i; A# R: ?0 x! h3 s
76
. h2 Z2 p( E+ B: N- o77
$ @' q) ~* p. }* G' B78
7 T* {7 B% f6 T& `- H79( i: v8 c7 B" P
803 g, D- c* }3 U9 G( |$ ]& |
81
* }/ j# Z# b- M( f5 U6 I4 H- s2 ?82
- D* ?5 ~& ^' L% o* v. m836 P1 d. ?& b/ K, ?6 E5 K
84
+ `! T, L1 _# b. G85
: X' C: {/ L" d0 [86
( M% ]' t1 K( c87
6 ?' n) Q! u9 O# E. Q: N- S88
6 _( C1 |4 C; z, l0 _89
; e4 f4 Y4 W' W: ^& |90
  M! v1 X8 d" {! J912 O* \8 A2 A; n8 h2 x
92
8 z  W% h4 [7 y/ m3 u93
; t  k3 Z2 b+ ]4 l8 x) P) F" u1 x94
, N9 `: l& L! X; u9 T( d957 W# Q( @& s+ s0 S& V, B
96
2 G. `1 z0 H( c4 J3 u97% i. d4 G1 h6 F% e
98
/ ^  H! p4 F3 r- Y0 V99, b# O2 z2 m. N
1009 s  T- V( h) ]/ j# d  [7 }' N% g
101
4 U. }) M! n4 X$ ?: @; R: I" ?3 A. b1021 a" j- i- Q, U/ ?
103- m4 E$ a* t. R/ f2 i% s% y
104
  L7 n% f3 z% \: J! K- i105
8 f: @2 [- y8 B1066 L8 I3 V) g/ |% I- s9 c6 Q! j8 J; o$ T
107
( ]6 Z, h9 {! e2 ?6 v108
: v3 ?/ M+ U/ S. P109
/ }" }3 S# o) P9 P) c! L1102 a  m) t" p3 P! W! ^
111
; H* c! m" t" a" i8 u  w112, N: @! O) V% r# j6 L, |& m0 Q4 t
113( g% V" [2 c% n3 W8 `" @0 S
114
: N1 g2 u, }. P' s9 e  F8 V0 L115
+ o( O3 h" v) ^" S* z$ z116
+ z) z+ e. p  g6 P1 i# F117
: O: x; O0 u5 |1189 F' k& s2 Z" r: y. E
1193 |& ~% Q) D+ A& n9 b, B2 d+ c/ ]
1206 l( N2 C( O: o1 Z2 p
121
& x: |2 K  P# d, D7 [- C6 V122
% |4 R4 Y  R3 D2 T- \! n+ D- ^1234 I- ^+ ~: V& T
124
# Z- v! }- R' \* l& H: j3 v$ {1259 H9 K+ ], Z$ T1 ~
126" J3 O& X1 B- w0 f
1270 P- D7 i; K" Y7 Q; p  r! Y7 C0 _
1288 ^3 h6 x" j& d
129# I5 n  R' d; W8 z  f# Y9 ^# H
130
- n" \1 e# E, o7 N1 `' R8 A! A. e& z1316 L' L! m$ ^. ]8 V( H; }
132
2 K" q  a6 S1 ^( W, q& H1 {! C$ f0 q133& @0 [3 N& ]  `' n- [
134* ]7 }# m% K( ], Z
135
1 W$ n' k/ J5 W' n136
! X) g& m  \  a, B1372 m: E: `  ?6 h8 q2 Q+ D
1383 p4 o' \6 t) [, ]; _  o" t
139
& C$ ]+ N0 m7 o; m% w4 e+ W( h+ B* I; X140; p$ f" ?: Y  ]) t" W
匹配问题：
3 H" {* A% j% G7 H% d" g
3 G6 s& K" h* h" z: z最大匹配：2 q7 |- L! I  v8 E: L. o
1 ~. c4 I) e+ \6 T4 X/ H, q  N

代码实现：

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