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

最小费用求解$ r+ V  J4 u& d

, B7 M4 F! u# x' {$ ]3 P5 FLingo:
5 ], q6 H) ~4 L0 k
4 `% z  D$ T( D. W" j3 pmodel:
$ e3 t$ ^2 L+ Lsets:
9 i' b! D3 ^" I% R" rnodes/s,1,2,3,t/:d;
3 i  ?- L1 A3 `5 P% ]& tarcs(nodes,nodes)/s 1,s 2,1 3,1 t,2 1,2 3,3 t/:b,c,f;, U" D! n& q' j9 M5 p
endsets6 b# e( K- `: E9 W) N* f
data:+ I6 i6 @& E* L
b = 4 1 6 1 2 3 2;
9 ^, }. ~- T! V7 Z! t+ h( \c = 10 8 2 7 5 10 4;
2 V9 c5 f  U0 D) G$ U! u7 _d = 11 0 0 0 -11;0 D0 J4 G  X0 |2 ~
enddata
$ I4 L, d; u9 E8 o. x, gn = @size(nodes);( r: }! r- [4 v  u" {
min = @sum(arcs:b*f);
9 w" S0 @9 ?4 H9 U. }@for(nodes(i)sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i)) = d(i));
4 f: p( s0 X" ^5 s) k% q$ H; x" q@for(arcsbnd(0,f,c));
4 s2 t* P0 Y; h- ^9 K6 k" w9 P$ jend
; }- G4 k8 D/ n1( s9 n; K9 ]; I# w  |
28 Q0 T- `: P; ]( g
3
4. r# F6 \0 N( G3 Q# L, a
54 n- l/ T$ D/ X: O, D
6& E: J) E: A6 @/ K1 t9 r& R; P2 f
7
+ P+ }& t6 |( _8% z8 j) d4 [% B# F
9
- e8 N# h; H" b4 c, E$ M$ f( C6 {10
2 i% Y( u/ X6 ?. h: [& X11
: y/ F( P! O! s$ ?1 t12" w7 M* v" `+ e1 @* R
13
& W/ ^# {- S/ U, W2 f14
2 X, n/ g3 J2 P8 t* I" }15
+ a: s% \) V1 _; }) [" M8 RMatlab实现：( {  m& J+ F# M% x" j8 w; ]3 ^

, s; p$ w; y' n" ~; j0 R! L
+ S, W/ M# ~. q4 [n = 5;
; M3 K/ J# o1 _8 r%弧容量
2 S% f$ q6 \6 e) S6 C7 g+ V! \: u) Ja = zeros(5);; T6 \8 x2 J# j
a(1,[2 3]) = [10 8];
/ G) I' n- x: R- D5 w- n5 ka(2,[4,5]) = [2 7];
7 n& H& o6 N( b: u1 X: E# ^a(3,[2 4]) = [5 10];" V" r' x7 d; e
a(4,5) = 4;
& u- ~9 x9 c& _% X/ Z! ^C = a;8 ~0 g& y+ N0 x+ W' r' |$ Y
%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];3 _  g5 s* R# w$ |" i: ?- q
%弧上单元的费用
  }) S3 Q3 \7 W. v# F" d0 q! ~a(1,[2 3]) = [4 1];
- P0 y( A- [& k8 }; ka(2,[4,5]) = [6 1];+ V, f0 d; i% _/ g% D
a(3,[2 4]) = [2 3];
. t* Z1 m# V: y8 r* A2 ~. qa(4,5) = 2;% Y' D! `( g/ p5 q9 x
b = a;
& }0 `4 }5 M3 h' I+ o+ m0 L5 U0 ^3 z%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];
: n( s; L! f$ \: y  H( R1 A0 C%wf表示最大流量。wf0表示预定的流量值
9 q: v1 L, m! V- S2 S$ _) _wf = 0;
5 u  V  F( |8 U8 ~/ |wf0 = Inf;
4 C( `7 i( k% L. W%取初始可行流f为零流
& I% z* M7 F$ Ifor i = 1:n. t. f+ L' }( z  u( r' C
for j = 1:n
+ Z4 U  ~6 [7 _3 S       f(i,j) = 0;
( @" t) l" P8 O, a5 [1 r6 ^, W end
! [) a& G/ U8 }+ S0 Uend" Q, j1 t  b- B& P0 u. G
while (1)) v6 e9 D; [( h6 ~( \
%构造有向赋权图( ?, y( t7 F- {. N
for i = 1:n
: [2 M* r. k" R+ M; W- d. `. H! H       for j = 1:n( A7 L9 v/ Y. A  g
         if j~=i5 @, a3 ]( q+ C/ u) F1 h' M4 P
            a(i,j) = Inf;
5 U" _, \" e2 d) c9 a0 |          end
1 b  h; F6 r1 l" r. T* M" ~       end3 G" |" T7 m: I3 q. I
end- M0 ~# F2 R" n
for i = 1:n' h0 R* p! F" h/ G% g( g
      for j = 1:n
: [- A7 J) I# N* j$ A2 s          if C(i,j) > 0 && f(i,j) == 0: m; A+ [8 n6 S) B0 L$ C$ {
            a(i,j) = b(i,j);7 G+ w; X: j: E+ _
         elseif C(i,j) > 0 && f(i,j) == C(i,j)
$ S" n0 i- y* b             a(j,i) = -b(i,j);! b* o3 d9 G0 z; ?9 R+ ]% b/ a1 ^
         elseif C(i,j) > 0! k9 R4 }! w) ]4 u% u. s
            a(i,j) = b(i,j);7 o6 [6 G$ k' A: f8 X
            a(j,i) = -b(i,j);0 e7 @- m: U  G
         end
% w  X4 ], o6 \) h* r1 Q       end
# P& j  E- I4 d+ S' Y9 j# J end& L! Z4 \, p: E% Q5 h% y3 i
%使用ford算法求最短路，赋初值. K) J. W- s+ Z& ^
for i = 2:n
4 E/ E6 l4 G6 `5 u5 H3 j       p(i) = Inf;
! j- F8 A* r/ r9 L       s(i) = i;6 X: q0 e* n" \3 r  t3 @; z2 g
end
0 [( i' R* A8 Q) }) H" D %求有向赋权图vs到vt的最短路,赋初值
8 D. B! j6 P) `3 N for k = 1:n
; G! j' i. O+ R* A       pd = 1;: b& U' W9 X1 z' E" F
      for i = 2:n1 ]- w" h$ h8 p. [: b) e# k
         for j = 1:n# C- S8 l+ Q% `' o
            if p(i) > p(j) + a(j,i)5 Z, d$ |6 p: L5 C
                  p(i) = p(j) + a(j,i);
) l5 I1 X8 k2 }- B! w$ {) L                   s(i) = j;
* }* s3 A; N& w; |                   pd = 0;+ z. }. N2 J8 f, Y. j% o
            end
& P/ j; S% k' t8 e' w6 \$ n7 I& ?          end% n) z7 _- @: ?0 t
      end5 ?) i( m/ r$ v& I, j/ v
      %求最短路的Ford算法结束
- Y* I1 V  i7 c9 {' l( E) L3 m! z/ o       if pd/ a/ q. R2 b4 r! E2 k
         break;( ^: h! @+ u  j! X' X
      end" j. P7 ]+ e* ?8 R" @
end4 Q+ k4 j6 s: W+ g0 ^1 f
%不存在vs到vt的最短路，算法终止，注意在求最小费用最大流时构造有向赋权图中不会含有负权回路，不会出现k = n* h6 I' E7 C% d( g6 }# e
if p(n) == Inf
7 U2 _/ r& h2 b! ?8 R% s* P" E       break;& i: Y& F! C" s, R* P. R
end
/ `9 |' i* {6 K4 t2 s  Z %进入调整过程,dvt表示调整量
* P1 y5 Z* m$ e9 S( W. s dvt = Inf;6 {: j. h! d, Z2 Y3 ]1 U, E
dvtt = Inf;# X3 m$ d4 T* f- t5 U3 q% v
t = n;" `5 y) }1 m6 a0 T, Y' N2 _
while(1)
8 J5 E' [2 d# A7 n1 n1 h  W0 J7 ~       if a(s(t),t) > 07 f" e2 A3 f1 x. c  g
         %前向弧调整量
0 W. Y# ?/ f3 R4 I$ t          dvtt = C(s(t),t)-f(s(t),t);5 \; Z) |' ?) Q! I
         %后向弧调整量3 g7 d% X  Z  d# I
      elseif a(s(t),t) < 0& p, J7 v) ]. i0 X, a
         dvtt = f(t,s(t));
3 M( ?+ q: K; \$ ~9 n1 l       end' e! }+ k+ {! H; j& f6 ~  }
      if dvt > dvtt" D' D5 T3 o9 r( u4 a
         dvt = dvtt;5 ~1 N6 ~/ M& X+ B
      end" k$ U- @, `5 S# M: n
      %当t的标号为vs时，终止计算调整量
: z: J4 E+ t# `( G0 p- Y6 B       if s(t) == 1
8 p  X! s" u) x$ Z$ `          break;
/ ~* U  Q( e- z6 y' m       end
4 a' x+ ^8 u# D) V& `6 |- V' l       %继续调整前一段弧上的流f' V' Q8 _. u* H/ u3 a/ p# G3 ?( x" Z+ V
      t = s(t);
+ E/ B3 q+ x$ u9 O end' l& N# t* E4 H2 ?& M2 a, P
pd = 0;) v/ C" v6 c  w) n; E& _3 {
%如果最大流量大于或等于预定的流量值
+ K, W6 m6 W3 V: S if wf + dvt >= wf0
# D& \; [* `- r1 i: T6 {- P0 _       dvt = wf0 - wf;" b* t$ u9 L/ \- l1 M) e
      pd = 1;
: V0 F6 O5 _5 O# U! [# z0 Z end
8 I8 ^) p# N. C* l. n) {  v/ }( Y t = n;
% g3 A: \' p5 X8 S  n. p %调整过程, {! N. x0 \1 e  u' M
while(1)
0 Z8 B2 f  o4 j$ m if a(s(t),t) > 0  {' S( A& n8 I* a1 g
      %前向弧调整+ C& x6 w/ ]; ]  K5 q: b) \* j
      f(s(t),t) = f(s(t),t) + dvt;
5 I+ s! n0 j( t6 u elseif a(s(t),t) < 0  P+ M5 a1 Z/ V4 Y7 v
      %后向弧调整
8 i% @* A. ^$ C- k: l% K1 @       f(t,s(t)) = f(t,s(t)) - dvt;! s$ C0 N* {) ^' _4 W, V
end
# x) O$ \# }4 J1 L5 { %当t的标号为vs时，终止调整过程
6 C, O# l" w6 ] if s(t) == 1; y4 o$ [# \2 |3 [2 F
      break;
8 |) o9 ]4 O" b8 T. }6 A end
. c# }& J- _4 w, [" x t = s(t);
  g' l  g) N# R7 K& h+ z end8 W, C0 p6 K) G
%如果最大流量达到预定的流量值
1 R7 s$ C& O7 S7 q; R% m if pd
2 X6 c- p9 o* D: k; f" V       break;
# N2 e7 U) a& n  a6 Q" {) G end
: L/ M8 U0 O  c2 H %计算最大流量
5 M( x6 i. n! a' B7 o9 J wf = 0;' v+ V, g5 i& O- w! W/ M  G
for j = 1:n
9 t& J7 p. v% [; x* z       wf = wf + f(1,j);
( o" P. K0 C) ] end) H1 |( C3 M" [; ?! n
end0 B: k' e4 E* t( E3 Y1 k: ]' C
%计算最小费用
zwf = 0;
& `2 E7 N/ j  W' Dfor i = 1:n
' M2 E7 i$ K. }. v/ J for j = 1:n% U0 v& r1 h/ H/ }% |
      zwf = zwf + b(i,j)*f(i,j);
& y2 O# P* V: }* V" j1 s8 F end& T4 P6 p, U0 a5 B9 f
end
6 q& p/ d2 O' m9 c5 _. h%最小费用最大流
- ~  s- [( w6 Y- x# K1 ef; k* q  d& r" {! A- r% m
%最小费用最大流量: i. {1 V7 t0 x# ]6 \4 o  N, l
wf
1 ^( X6 K' N, _. g%显示最小费用, Y3 `& M9 ~* {8 ~' f: q
zwf3 G, v. C0 X' [: C7 M1 M
1
( Q5 J; {8 q9 S4 P1 L2
' ]! P/ G7 f) l/ ?6 m5 l# a7 l0 J35 J- b/ o$ S- A4 P! I8 }  ~- E3 X( l
4
- J! a# a8 I' b9 g1 [, i1 _52 R8 N, x  \+ x
6( }" i4 b- K1 j& B
7, V' J6 F- w0 {; }
8( U7 ^# |& f; c" e
9
2 y- l0 T' \, R. o( ~! e& o10+ u% T( p4 C$ {  M% L' n* G
11
/ N4 e; M* g# X. J# u: {12' P/ x, a# j  g( T
13
' s3 S/ h2 Q! r2 {* L3 q, Q141 O& o9 S1 q) I# p
15
/ M; W9 n/ _4 s  I6 e: E! ~5 g0 K' R16+ E0 S! K( h3 v  z# X2 L. g/ y' \
17+ x# D) e/ e; L, J
18
6 ]. ?% T& F$ d$ t19
- _; }, a4 t9 ~+ v$ O20+ s1 C' |' }4 L, D& q$ d
215 K! s6 n' ?/ b+ C% e; X
22
7 A* P' P9 v5 Y& R238 l6 H7 j9 I1 W- g# ?8 e4 `: k
24  z7 B$ \% l7 A5 J$ J( _9 f
25
+ p: D: [. x$ M* @- P260 n0 ^; O& P* a* a* O1 I
27/ A( u4 S! _8 z1 X8 s6 t3 t
28
' |  S; @4 U: o, d2 L3 w29# P9 z/ p* I& u& x( d, j' `6 e6 F) P
30( Y7 R& c* x  e% P
31
0 E# |  r/ }1 p/ S32
3 \6 s0 N7 u" [& X33
. t2 V! ]3 Y: ]% n" W34
8 f* {/ ?# q4 x" K. O0 B! J9 \2 X# O. M356 y5 P; i& C5 l, q% `8 u+ `: D
366 O5 @2 [7 J* ]% P
37
! m" b/ o5 ?; b! E38) {3 I+ Y% C5 Q+ a  W1 Q
39' L$ z% k" l0 U+ D5 T6 v! n# L
40+ d8 M5 V" S  \: U' U' B
41
' j/ l# x# G; s2 X42
, q; a' W6 N6 ^5 a, B43
, |  z; D' h# D' z9 G' \6 A44
8 ]% @) P7 {4 C: `458 Y1 O+ e1 ?; D
46
- E1 R5 c+ P7 |: N( h9 U47
0 i' M2 p1 f3 ?; c48" B- X. W" H1 I
49
) S9 Y7 s* R9 d50' _' o& o  F6 E# E% S' B8 r
51. E2 T+ U% y+ G4 q
52) S3 [( d6 A( `9 ?& T* T) r9 g6 m
53
, Q$ X$ f2 a- }" _$ C54+ U* W! T3 M% n  c8 A' l& P- z5 @
556 f, d! T& ^1 a2 C$ Q7 T
56
' N( {( I' e: b# r2 W! w57
& w& b" y9 \, _2 d; Q; g# Z58
5 ~3 [( L+ }1 U4 A" L593 i4 @, Y# ^1 _$ j' B* `: b
603 N! o9 S3 A" U( S) Z
61
) c/ U, K9 ?2 Q0 H, H62# N$ i% r% B9 i7 C
63
( b; [( E1 @8 f- \" v$ O! ?( Z649 W' b7 w& X. P# p
653 u+ X8 p9 K! J3 Y. l4 r8 _& {
66: D' F7 i; P$ {& f8 h/ I
67
5 m3 a& k! ?# o1 U683 C0 \! P! \3 v
69
0 V) |' e9 U( _4 F: ^$ r1 Q70% w% E8 m7 A8 q  {
71
% x: E$ b8 h/ b* o& ]722 G1 s6 ^+ ]3 ^1 k0 F/ [
73
7 @# E+ I! H8 N7 h5 B+ v74
  x- j- D" M8 C75
$ z( m  h$ B( u( D/ H' P/ k76) F. ]$ s! R7 H$ Q5 }8 |$ h
777 y6 i0 I/ x6 ^
789 D) p8 Z" e' |$ o) o' h; `
79
2 y  p5 x0 n+ |8 x& V2 _2 z80
5 h* J. z, {% I- p: K  R818 v: S7 _$ f% X1 S. y
82& H+ h9 K5 P. N+ x$ Q+ c% }
83
: a$ Y# G# B2 `: m4 w84" L+ }/ o/ A9 |  a5 S
85# c" F6 S  g* Q2 j9 U7 E
86
  q5 T3 g$ Q! h87
0 ]' Z7 b- @1 F0 I: `7 T/ x88' h1 v( I. b9 Q( @1 H- ?
89
+ }4 e; u- S) i. O" H90
  S, n! j3 ]' {8 I91
* c$ h' k8 X9 q92
/ j9 g7 l9 P0 y: j93
; U: \& K2 }+ N  q- S94
! @0 N# E: Y! a: _. n4 |95  H" L9 b" I/ m$ Z! a% r$ D$ `1 G9 ?2 D
96
- p4 Y$ r# t/ e5 n97
/ ]% b; H/ i2 A! f: ?981 d" P9 n4 M: g
996 U7 m; |+ n8 C% K) F0 _
100( m3 U! T- b* b- M6 R' g6 w$ K  y
101
0 O$ B! y( L% M. V" b1027 g4 F: `0 f7 Z' [5 R
103/ n5 @* v/ z# c
104
$ A1 R' o$ t) d( h8 n" {. \/ T) {1056 [- }& X& s8 ]* k( G, e
106
4 E" ?! n* {3 V% s8 F  S. T8 f107
) N. f# A, S9 R; T108
/ l0 Z/ @7 W% z8 F9 N109  _& {! z  Y3 |+ d. t
110
4 w* y# l( u4 n$ i+ o  ]1 e111
# ?2 x' d5 _- ^% u, t- x) p1123 d# V& E- d3 {/ G! c
113; A5 H+ L6 c% Q. X3 ]  e7 q
114
4 f5 T0 L4 U- }7 }115
/ L8 f* \5 f/ ~* \, {116$ m: E; E3 C. h3 s
1178 d" U, D; K, ^7 f3 L
118
" o7 W3 U1 X( e  O# n119
+ `5 x, \9 H1 G120
8 ]! `9 M! g/ T1217 J: C  f6 a: b, P/ J# t
122% P  y3 [3 u0 ~+ i  w" ]
123
" x( \, Q5 p+ J0 v, V1249 J. E, O$ L( b; s
125+ @6 j, p0 e# M, G! U3 X2 ?, O+ Q
126
; l# |, Z! K7 i$ o) b2 h4 F127
5 ]0 l" i+ E; b  Y! {- \! ?128" X6 N" e. w) j8 r0 D3 |
129# G. T" e9 X& A- V( Q
130
% E! o8 s# U# W/ [- p/ `& E131
2 J* J# B5 p+ U132
1 \$ ^0 Z! A. F  {5 \5 O/ s133
$ m' S3 E1 ^0 i1 ^6 r' d8 m134
9 r5 l* C8 C' P135
3 @3 c- i! c6 w0 O: `" O136' i0 C) N7 g+ |0 `- l* _& f
137" K' P/ a% O$ w+ M8 g
138  G8 |! ~) _" P* p( N# E% d3 s' p6 M
139
) R0 c4 X4 c  C% `) z: H140
5 ~' _) W0 {% @) `' b9 T匹配问题：
" i7 Q8 q; U$ Y6 N# [: A! Q/ I' L
最大匹配：  [% v' D$ c2 W/ v

  v7 e% l9 a: v6 O) F" d
代码实现：

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