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

最小费用求解
( G4 n1 M$ Y5 A0 V9 i- ?: X- L/ q4 m% \# p
Lingo:
+ w  i; W7 e. `7 D7 [: `2 {% g
: z, [, U9 z  Z# Rmodel:! e3 l1 s6 M4 C5 h4 v' W
sets:
3 n& F6 `: B" V$ h- V; o( X- unodes/s,1,2,3,t/:d;4 N: |% @( e1 ?  ]3 H3 P
arcs(nodes,nodes)/s 1,s 2,1 3,1 t,2 1,2 3,3 t/:b,c,f;; X; b! {: Q5 R9 \8 K& S$ t
endsets" ?# K; r( ?" B: u
data:
  B5 E; w& ]& @: Ob = 4 1 6 1 2 3 2;
; ~, r- Y& t, l2 }- X* uc = 10 8 2 7 5 10 4;% ^! r$ b4 w2 W; G* g
d = 11 0 0 0 -11;
9 V# e' p8 _9 u* {8 venddata* \+ L" h1 o/ s
n = @size(nodes);* C9 J& v! I, o" h
min = @sum(arcs:b*f);0 B$ ~- x) B, |3 }, R( J- g' h
@for(nodes(i)sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i)) = d(i));9 D+ {4 o, h" W
@for(arcsbnd(0,f,c));
) w+ V1 M+ u% u0 U. x9 Z  h/ g# qend
0 c% o9 x1 |- B19 \7 @4 M7 i& s/ p% Z
20 c6 C6 _) A5 L8 P4 l0 u0 P
3
/ {/ }" N( n3 I8 X8 f- X/ f4" ?7 g) _6 J0 I1 c
5+ z! w1 _' l* ~4 J) u0 \4 f) C
6" g% m6 Q. J  J& [6 y1 L) `& J+ Q3 B
7) c3 c) \, e0 G9 M2 q
8
2 a' J; i. q5 T  |9$ l8 e. R1 ^, S3 l
10
( n3 F+ n5 b* S: |& d' D2 J11
- Q1 Z: m6 u6 p. i8 E124 j, C" Q8 h+ k7 v0 U! b
136 \8 y  F" ]" |* _4 v; J4 H
14& m5 X7 j% y" @4 w
15
6 p, E, g9 k4 ~% @Matlab实现：
+ h9 ~$ f( R9 b9 D, N) P: z6 r' L; m
3 h9 r7 `. G6 T6 E
( r5 Y. a+ G% q- E% l7 Vn = 5;; N. e: R. W& k/ X1 T7 C+ h# M! t
%弧容量
, K$ u0 E7 d7 w: p2 I( Ka = zeros(5);' l9 }7 l, o, h5 A" l
a(1,[2 3]) = [10 8];
9 G) N) k# s. s' r$ [a(2,[4,5]) = [2 7];3 [+ ^% |* c5 I
a(3,[2 4]) = [5 10];& k* K& u0 W# w* [: b& h
a(4,5) = 4;
& ~. t1 ?' M$ X' ~5 rC = a;( g& S  k* y5 _' [( A
%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];
4 X: ~' j/ }5 {/ i; _- k$ [% L) A1 E%弧上单元的费用, C5 \. j; [' D7 C! a
a(1,[2 3]) = [4 1];. _8 O8 F( [- S1 E7 G# ]2 D
a(2,[4,5]) = [6 1];/ V/ `9 j: ]: r3 ]  O4 V$ [
a(3,[2 4]) = [2 3];* Y0 r  N. |* ^0 ?! d- V& S& t
a(4,5) = 2;5 F3 K* e2 d' ?: B2 [. R3 b
b = a;8 N8 V- v( _4 h' w2 a, C% A. D4 y
%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];
+ `* q/ l: B) d7 E# M; r+ r%wf表示最大流量。wf0表示预定的流量值. G/ J7 X2 G( g
wf = 0;' O7 k+ V% F- l" }- s8 r
wf0 = Inf;
- U' i9 A+ \. B6 O%取初始可行流f为零流
" ?7 A; N7 k8 ?& U* h+ e, afor i = 1:n0 R7 y# ~8 r- g, j; _
for j = 1:n
: |& e7 @! }, r! Y       f(i,j) = 0;
. A. i' f, Q. j- [4 s8 l end
  H7 x( }. \/ f% k; X3 Pend
( A3 z7 @3 b) u9 q5 k  Mwhile (1)
7 \- `8 c( l4 X7 ~9 ^ %构造有向赋权图$ `# u, I! ^* q- J5 \% }
for i = 1:n1 [( F! h* F6 F
      for j = 1:n8 O2 z+ [8 D- N* _9 M" }% `, `
         if j~=i9 X3 B0 Z' e) G/ t& n
            a(i,j) = Inf;
9 P9 B( g9 s1 y1 P  Y) d" _/ m          end1 c3 M. s3 s; \3 l  E, f" I
      end
+ `2 d+ y& P  E( Q6 v, q- T' M% ]5 W end. r3 w: {( b3 b- _* }, `
for i = 1:n3 Y- H6 T7 X7 k' E/ }
      for j = 1:n7 p6 k2 v) x+ N" ?  N
         if C(i,j) > 0 && f(i,j) == 0
* n' d% }2 K/ a, o* i             a(i,j) = b(i,j);. d; u+ x! N7 E5 J& U% k% _
         elseif C(i,j) > 0 && f(i,j) == C(i,j): ~! Z2 k$ I5 A* P( B9 f
            a(j,i) = -b(i,j);
" A4 F# O! x9 ]- `0 N' i, U          elseif C(i,j) > 0
! h. V( j0 Y9 T+ e8 u& J             a(i,j) = b(i,j);6 ]9 y2 |! V% p# O) T6 Q( _- @
            a(j,i) = -b(i,j);
7 J( u2 L0 @0 ]+ U; v3 m( Y          end+ E; _4 f/ ^$ d% V2 G0 A
      end
: b$ Z$ {3 X: r% {) g end$ ~" `9 D: @+ A
%使用ford算法求最短路，赋初值
& l7 m9 |1 ~# H# o+ h& I! s0 T, k for i = 2:n1 }! z7 p7 H7 t4 S; x
      p(i) = Inf;
+ O6 e2 |: Z3 A3 }' Q       s(i) = i;( h' A, b$ i9 l( W7 n8 l
end& n: H( ?+ U& p6 ~& B- X
%求有向赋权图vs到vt的最短路,赋初值
. q0 N' s; [* ]& X for k = 1:n
# m) Q2 Q+ U& }6 {& V1 L       pd = 1;
% W& z3 P' Q0 H! f2 Z$ q% h       for i = 2:n$ ?0 I8 B6 k4 h" }( M
         for j = 1:n0 `. G9 q6 O* y1 l1 y
            if p(i) > p(j) + a(j,i)
4 W* Z7 }) W/ d; J                   p(i) = p(j) + a(j,i);& A. d# I; W. ]5 a
                  s(i) = j;  Q5 ^3 u4 E$ N  B
                  pd = 0;
; Q( ?9 b, A  C" @; s             end2 g. d" M2 g5 [" J; k& h, t: W$ Q
         end
6 b0 j$ ~# Z" l0 D9 Q       end' I' q# r( o+ Y' V" X2 a
      %求最短路的Ford算法结束
) w3 N& B% F3 k# Z* i) w: O  {       if pd
1 x( [. q; S4 j8 b' c4 I# c          break;
+ L- p% K) h0 x       end
( T6 c& b- f9 K end
' g3 v" \8 S6 X; c5 ?: s- H %不存在vs到vt的最短路，算法终止，注意在求最小费用最大流时构造有向赋权图中不会含有负权回路，不会出现k = n
/ q. k9 o" S! V- z2 D if p(n) == Inf
7 @" ^! I) j8 D# @. w       break;
$ x" Y# U' `+ U# |. |, Y0 Y end
3 p& z* A9 J) M4 ` %进入调整过程,dvt表示调整量
' _& K9 |8 ^% d4 o2 \ dvt = Inf;
" }+ c' a4 J4 U; O7 ^& R dvtt = Inf;8 K3 f& N& I2 g% Q6 R4 ?& U
t = n;* C% }9 f. c8 R& e! `- q5 E
while(1)$ l" }( j& |9 K" f9 m' {3 d
      if a(s(t),t) > 0
0 z% H3 J! B# x- q5 \  r  b          %前向弧调整量
& E# D  e" |+ O- ~          dvtt = C(s(t),t)-f(s(t),t);' k6 Q& M1 H  f' v
         %后向弧调整量5 U4 B$ q0 }/ [$ a# y$ r
      elseif a(s(t),t) < 0
) M, A( W" \! Q% Z# E          dvtt = f(t,s(t));! n& e% L8 a7 A/ {, H
      end
: |2 B* C, t/ Y       if dvt > dvtt1 ?2 ^- M) H( F* e9 }6 S/ Z
         dvt = dvtt;
) c8 h1 V4 N4 {# \8 }' P' g  h       end- z8 Q& K) r. u  \
      %当t的标号为vs时，终止计算调整量7 i7 S# s0 W  p; q9 _" N
      if s(t) == 1
7 S* x7 X5 `) |) p7 [+ O          break;
" t$ m9 r) a! @$ P- r/ E       end
: \! ^2 F5 r2 Q/ y  {% T$ B# b       %继续调整前一段弧上的流f
! q- _7 i7 \( U& A. D. u9 |' `7 W       t = s(t);# `, e& ~0 b, N  E- Q" A3 h! N% P
end
- R9 Y8 p1 n$ B# n pd = 0;$ _8 X( B5 D3 N# K, O0 C
%如果最大流量大于或等于预定的流量值; i9 h$ Y& b+ x. ^- E# U+ `* k- w
if wf + dvt >= wf0
8 |  o9 B2 ?/ y8 R! l" N) u& U       dvt = wf0 - wf;9 S  J  H; N( p( g: s3 f4 R  Y
      pd = 1;0 U+ c" R8 D6 y" L
end( n1 a& [9 }. Z* U, l0 Q
t = n;% ?# e* X( O) K
%调整过程
6 s3 T+ v: R' ~7 Z/ \( P1 K9 | while(1)  Q( I/ }' N% D! |
if a(s(t),t) > 0! H- r' Y; \5 s  g% @5 s
      %前向弧调整3 j7 o. |. K6 n+ k' E( _" `8 k
      f(s(t),t) = f(s(t),t) + dvt; - T2 y8 J3 }  B2 }  l* y) T: n" E
elseif a(s(t),t) < 0
0 J5 Z5 [: [# Q2 J) b       %后向弧调整, p" V" [6 @$ m& y% v$ P
      f(t,s(t)) = f(t,s(t)) - dvt;
6 c: L3 X' z0 h! W end
9 |: L; r2 b5 ^5 W %当t的标号为vs时，终止调整过程$ a- j: i+ X4 t; `( d: _$ v
if s(t) == 1$ Q+ c; G* A: Z% g. S
      break;9 ?/ n( R. D2 [* L
end: ^& W/ f8 C: q
t = s(t);, f6 t  e- [: @- i, ~9 d
end
7 @7 M. N) \/ f- D7 l8 J+ _ %如果最大流量达到预定的流量值
0 [3 ^$ G! Q! @ if pd2 m2 l1 `. M4 o2 b2 Q; {) T
      break;
/ G* f" r' K% { end
8 p$ q) d5 R/ `. j %计算最大流量5 B1 \- Q( J3 I/ c+ _# ]( E
wf = 0;* @  G* ?% T; w! C, o1 l9 N2 i3 j
for j = 1:n
6 P+ R" i/ v2 b2 b       wf = wf + f(1,j);! T( t  V1 X4 f& ^5 ~- |$ M
end' n! e* A7 S: q' X# H! Z0 g) `+ n
end8 `/ t/ L, V9 k/ ^
%计算最小费用6 F' t6 Q/ [; S
zwf = 0;
' s7 \; ?/ O9 L4 s7 Yfor i = 1:n
8 A4 G  k, w3 o( u" `4 s for j = 1:n
6 A- e* G' u5 T9 H6 ^       zwf = zwf + b(i,j)*f(i,j);2 k+ [/ f% k5 ?9 b# h' `- V
end+ z7 t7 |8 _9 t/ n' P  d
end
* B# \$ h* O( h! X- v7 t( B4 n/ e%最小费用最大流0 y) V- z( W: m# M- ~
f
5 X3 f$ J! b% }+ y/ p0 c$ ]%最小费用最大流量
, a7 H6 o  ^) u( twf& A* C8 a7 v8 l' O. {
%显示最小费用2 T! X  T4 ]& t, K
zwf
4 v0 S! ], q( t" f; P2 O1( p. T- h" f, w$ ]% F
2
) K5 s; {4 ?. e2 [8 a* E3
" V( R+ Z% O3 Z" z  W49 ?" v$ \0 w+ z! Z$ P9 ^
5
* t& \" r' G. Z7 U6
7 s2 b8 ]5 r5 U/ t; Y0 q6 S7- Z# f( o3 F2 X" w+ I5 F
8
$ g! B* g( e. G  z5 b/ V5 H2 M9; p( s" r" E# W; }& ^
101 b! I2 t; }$ V; [
11
9 ~" j) _/ O4 ^. N12
, ]- e7 S! S: `. F' e( `" O13# _( T5 \" }' I+ P7 E
146 l3 q7 S. b4 }
15
/ _& {6 Q8 r$ \. T" [# j16: @' O) ]5 ?7 C  l& v  Y. N2 }" Y
17) ^; F& U# G$ X( n1 d
18
( @: }9 R8 J! N0 S( C" Q$ B19
1 o4 B( U6 H/ R20
& o4 x( F1 }4 B5 [3 E& o21" {; T0 E: d- o4 v' O# z# r
22& V! r: q+ S2 b; d# d
23, ^7 P& ?" s0 [, ?. `# u, t9 Q
24
# N; B1 Q0 E$ o5 b# D$ S4 u25; ]1 V( O  l5 O4 r% h7 I/ _! t) }
26) w4 b# v) O, m) b# r6 E& m
27
8 r% h! ^* Y/ ]' i& k( P6 s28! c( V, J8 a+ [. U
29) R1 y: o. M0 f
30
6 V1 e2 [- }/ S8 L$ W& Z1 ~9 F+ a31
( X' Z$ {9 G8 E  o) c& F6 F& M. h32! A# i1 L0 n8 Q! K3 S
33
. k- b) w1 |4 b, u8 a34  m* b$ i* g0 O/ {
35" K" s) W; R/ f- h# T& n
36
1 K- C; r! c5 g" ^7 W$ x& J, N37
" {. Y- L) A3 h  B- }* S. \38
+ U9 e9 N- {+ o! k9 P' i, A39  I' S/ T+ G2 D% ]
40
# t& _  l" x, \: @41# X% `' }6 H, ]
428 w: l, C1 Z/ f4 |. N5 S
43
# {: @1 |) ?2 V3 P44) ~1 F; K# Z1 p2 b3 s! j
45$ S9 q; v( ^$ K! U% a
46' F- _# p( X; \
47& w8 I4 O% [% t
481 q* p; ], C3 \2 i+ Q
49
& T- F9 n  F  W7 X1 W50
* n% |2 T. [8 g5 s51; Q2 M3 i- Y8 P
52
9 X" w# k! P' Z& T# X0 B53, {& G% y% u- {. I3 o3 r
54
6 t6 U0 j5 A  ^' o8 N+ O2 ?0 r  ~553 M" {: Y% I( e+ N
56
+ W4 b$ ~3 S$ z  }4 J574 x! B# w7 f; y
58
! k8 J0 p9 M! e" _  o0 i59+ s" L7 e9 ?; U
60
/ Y8 A1 r9 H2 `1 c61( b8 h! T6 O5 Q( `( x( F
62" F; T3 e) d& {$ I
63) }9 P3 @1 y' C& h, B
64
6 s0 P) u1 I" S7 ]* \5 G65
# i( K. e  G( B; e/ e66
+ ^: p& C, V6 ?" C671 L" n5 n& v9 q2 m2 d
68, Q8 f( i8 l/ |& P& K
69# b- J" S" [( Q9 {. N% r
70% x6 P+ ~0 o; H- v
71
5 B/ ~  o% {* y& L7 z72
, |: h, b2 _6 \  w, k9 Q$ |73+ P3 C& |$ L. X/ x( |
74
750 b2 w/ z9 b2 R: W/ B& |
76$ R8 B& I( Q7 V" Z& ^+ y/ w9 B
77
. ^" Q  R4 \/ C- P78
1 L+ Z. {4 j5 {% `' A/ Q& |79
. X/ t" h( J; K, r; \3 B80. X. m- S* @* u0 e- Y+ {+ g1 \9 E
81/ f, v6 M# A( n' Q7 F( X
822 H* X1 Z. E( [+ ?& a
83
% A/ a' n: W' M& {7 t84- y7 H2 }- K: X2 j5 q
85, }, W1 L# \; g2 J, E4 V6 m
862 N0 @6 k9 y0 P( S1 J1 v
87
  c  b7 c# u- q/ D& \" i88  h1 \$ [9 C. U( J
89+ T# X# ?& g  z8 ?5 u
90
4 E: S. R+ X, W) u. a% z  q91! y- ~3 V& G6 |
92
6 {9 j" ?" @0 K93
& n" P  k8 k! y! \; i0 G1 A9 G# a94+ C' y" m- V" y4 ^) j- u
95
. A* Q# \( ?$ y. r1 i" b7 j8 c96
! P. z9 f( m: {1 q" h9 E97
- C' c6 ]# i: e98
5 j, D  ~. V* I# H0 R+ W; v99- J; M4 @# t; q
100
! [: V! S$ Y1 C) j1 ?3 P101
+ Z/ T6 N; b, M& {% {" i102& F! \  L3 L3 N0 C( ~, R. s, I
103) C0 I! _! V2 p
104
  Q9 ^( s. K7 n( w+ \: d5 f7 ~105
& l. ~' h3 n! h8 e0 j$ x  d) {106( J) K- e  |8 T) Q7 j
107
1 O1 X/ i$ Q9 H+ x108. ^. P+ @# |1 ]8 f, g: @8 A/ \
1096 T5 K; {& o! R; n
110
6 V- _) Z. m( E5 P) A  s111' |( x- i9 f, G' X% i! G
112
# B, M2 ~* d; C9 ]" ~( I) F3 V113
7 S* f- A* X* R" Q114: e* F4 m) Y: C) A3 |, q
115* l) y+ Z, q+ O1 B
116
  V1 \" y0 [* J" {6 `117
* d+ ?  n( B8 n; q% E$ w. H118
$ c- z% J3 v+ U* d: Z, S) F119& U! H& H3 ?5 b5 }
1201 @  m6 ~) H0 r* z/ Y+ e
121
9 O5 O; a9 P: Q" W( K) u' @9 m7 N, t122
1 K3 J. d; K4 p123
1 n( T/ j/ ^2 T6 a. [3 K, V124
; W; h% d/ m1 s0 Q1 i# n125
8 V! b) w# ?) A- R1 L+ m126
6 D" r0 O% d! H5 F4 Y127" f$ o% G1 T( I4 @6 z
128% B+ n$ Q/ k0 P. P3 p! J  ^
129+ c5 ]' L) D; t( g
130( f/ p: |  a& X3 f  y) p; x  O
131( N/ Z  K9 F  P2 l. ^
132) @; p+ ?- T3 _
1336 K) U5 n  ^0 L7 A1 ?
134
( R; v9 Q. [! N- \4 [135
9 z( t4 C" C2 x  [- D9 F3 ?1365 q' n7 |: e2 r, s3 P7 }6 I
137
) D6 _' T3 r- X& x3 R138
- T5 I# x; U2 L& H% @' z139+ }. J1 H2 d5 [- b5 f/ H
140
3 s4 c6 n* b& t7 y4 P匹配问题：
, V. D0 I" j8 s2 G  C, O7 E
* K, }3 R  u% O4 j2 l最大匹配：
' [' H3 C! A0 Q4 m# S3 g* [( g& G/ t
& o; E6 K7 e, L$ Q# M
代码实现：

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