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

最小费用求解
) h3 }! m- O2 p4 N0 `( f  F  m1 b9 u) z+ Q2 r; s4 r% n3 M
Lingo:
: I! {; f$ q( I7 I+ j2 y& K1 L- x4 I1 ]
model:6 @2 _0 ]! |+ x+ O3 ^, I
sets:
4 V2 |8 R" [' Y2 i3 fnodes/s,1,2,3,t/:d;; ~$ x- n; {: [1 n
arcs(nodes,nodes)/s 1,s 2,1 3,1 t,2 1,2 3,3 t/:b,c,f;
( [- n8 I. B. N- F0 qendsets# Z3 q& |! ~! C8 o8 R. `
data:
# B6 X8 S+ r% k* g! xb = 4 1 6 1 2 3 2;
, q0 f$ R( L( o% G! T$ o  A% {  oc = 10 8 2 7 5 10 4;1 v  I$ O2 G! v3 K  R
d = 11 0 0 0 -11;6 W1 i8 ]( P& F' ]! m# @5 O
enddata; x. i' Y9 F6 {. s% h- Z6 `
n = @size(nodes);
! U1 c! p/ f9 a/ \" F0 |+ z2 omin = @sum(arcs:b*f);5 J6 q6 i1 p: Z- e: R' y4 k) A9 j3 ^
@for(nodes(i)sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i)) = d(i));
, T. M" {' i7 O+ V5 R@for(arcsbnd(0,f,c));1 [8 ^. p9 ^, s$ n: x! B
end0 J; Q( u. u4 c; J, d1 U# J& E
1( o. Y+ T) m5 u- J# H; l
22 b- p0 Y4 P' K1 \
3
6 S- C. t+ g/ |0 g; F4  ~5 c3 V+ ^* Z: B2 s3 P
5; h; W) ]1 k  P
6
8 v% A  v8 d1 X) F7) O3 {! Y: U. e) ~& x
8
* n, a) i7 e5 A% X& d9
5 S: p* O' r- N4 i% s10
  G7 F* r/ N3 O: B* j! I0 V, P& A* X11
9 E/ M. E2 w% F122 L, M# b& H5 j  M
134 B; F: w; x# Q* |5 d1 C
14  S6 {2 n) G4 c1 Q6 U7 o
153 F2 }( L0 e6 ~% \% `8 x
Matlab实现：% X$ Q( a: Z9 i0 n) n; R$ y2 c

! _9 x9 b. ]; h0 z" _: v/ p' ?( }
6 a9 k% I8 x$ ?" En = 5;
  ~8 `$ [) [: _1 r7 O! _%弧容量1 z2 |" ^6 \" e# s1 u% \
a = zeros(5);
/ L4 s9 |0 ~2 U* f& Ma(1,[2 3]) = [10 8];
0 g+ n9 l5 j$ m# ]: p+ Z+ g; p$ Oa(2,[4,5]) = [2 7];7 m: |$ k  ]' x% K2 ^4 v
a(3,[2 4]) = [5 10];
# _* L7 \  ~# i* R6 f% t$ G$ Ma(4,5) = 4;  ~/ f  t( _2 f- h
C = a;
9 z- b3 p2 z" O%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];
1 k) b" \! m* T! I3 P, w0 n%弧上单元的费用
. F6 ]; R9 N/ x0 o: `( n+ ya(1,[2 3]) = [4 1];
$ r& O4 ]1 J7 N1 ?a(2,[4,5]) = [6 1];2 l2 b; `$ Q6 O# @7 u
a(3,[2 4]) = [2 3];" A) u" m* c2 _
a(4,5) = 2;' I- L7 j+ l4 N' B, L6 E
b = a;
( G9 y) V5 @2 x$ r8 E- n& m  }%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  o8 t1 I! E+ ^$ H6 C%wf表示最大流量。wf0表示预定的流量值
7 h2 y" ~/ o$ A+ Owf = 0;
# x7 n' B$ E4 n% \4 d2 w5 P2 lwf0 = Inf;. B# v$ X4 Q% T, {/ N+ \0 p9 F4 r
%取初始可行流f为零流
1 A6 Y' ?: {# D7 T; G  Efor i = 1:n* D, ~9 @* [) K/ h% c& W
for j = 1:n  d; P9 U3 u9 ~" \7 Y
      f(i,j) = 0;/ E' ?+ b0 F2 {$ O( n- S  V  v3 E
end
+ s$ w  l$ y2 a* l4 b" ]% ?end9 a+ r' X  f( w" f$ m
while (1)+ \, B/ T# A! g' F/ P
%构造有向赋权图( A+ b: {8 f1 m- ^" u) t& s
for i = 1:n' @1 ~# L- t8 f/ G) s
      for j = 1:n
( f1 d- V" }2 q1 R1 v2 h          if j~=i
0 ]3 O% N7 V5 a0 I4 z             a(i,j) = Inf;; B# ]; n' \6 `- y0 x8 g
         end' G3 U, N) R# i% ?& h5 Q+ H+ a- i
      end  _: X  I) F, d
end
8 D- @3 m' A4 e0 R for i = 1:n  g: w- V# m: R" d  w- f0 V% j4 j
      for j = 1:n" v7 l# \/ Y  S+ ]4 s. o% A7 a0 x
         if C(i,j) > 0 && f(i,j) == 0' |7 {3 x- q. R+ z+ B: [
            a(i,j) = b(i,j);9 \! R7 a$ N- t( h# \0 m  k( d
         elseif C(i,j) > 0 && f(i,j) == C(i,j)
: t9 Q6 t1 j. n% T: J7 Z             a(j,i) = -b(i,j);: F! L' a! k% Y  ?6 u4 L
         elseif C(i,j) > 02 U+ k2 e" F- _! m
            a(i,j) = b(i,j);
; A$ w5 w6 f& b& ]8 I1 g; C             a(j,i) = -b(i,j);/ @7 P4 y) `- l  ?2 p1 W5 ^
         end  Y+ W0 l! v% V0 [3 F' V3 N  I0 E
      end
6 @8 O8 ^# k; F1 q5 I end
8 p2 S) r' j/ V  t %使用ford算法求最短路，赋初值) t5 W" ?! X' H$ g- q
for i = 2:n
# S3 W( ?* P- r4 m2 U' S       p(i) = Inf;+ ~( H: I, y5 u6 V
      s(i) = i;
" {' w' G0 Q) A- {% I end7 i) O: l6 s# ]2 D) S3 {) ^
%求有向赋权图vs到vt的最短路,赋初值% `* M7 g, w* y, v1 a+ t) K( Q
for k = 1:n, m1 N/ R) X9 t1 D
      pd = 1;3 g# x. V( i# Q
      for i = 2:n- N4 n7 ~. ?& w9 G+ X, c3 r) X
         for j = 1:n
9 \3 r, q: w8 `/ }1 G9 _" q0 c             if p(i) > p(j) + a(j,i)
: S1 ~2 U5 x. j1 Z' E                   p(i) = p(j) + a(j,i);. o1 ?8 h( N- k: B: b* ^
                  s(i) = j;
% ]# ^! x/ k* c# C3 d( _' `                   pd = 0;
, C$ E" a/ k) E2 n( e5 e* d             end* z( g2 D: z, }/ I# b* G
         end
2 r9 g$ h. J( M- f2 q$ J& R1 n       end5 Z2 g& a3 C, Z# e3 i7 `9 @' h3 T
      %求最短路的Ford算法结束0 e/ O7 X; N$ x- m6 y5 g
      if pd9 [2 T6 o/ [1 l8 w+ w. t
         break;! B7 [3 c) z" m; n1 l5 U
      end9 C" B) H2 X% @' X; a
end* G: c) k$ Y& K) ^* O
%不存在vs到vt的最短路，算法终止，注意在求最小费用最大流时构造有向赋权图中不会含有负权回路，不会出现k = n
/ f% M- g9 o( S if p(n) == Inf& u3 k, D! `* Y4 e9 R& V
      break;
+ }- o7 V  ^: l0 O# G4 b5 X3 U6 g4 b end( K( {7 S7 s) b( A% }2 g# h
%进入调整过程,dvt表示调整量3 ]) @3 }4 Q5 g5 ]( i6 i2 o2 j2 T
dvt = Inf;* G& q2 U" O) L( ~" n
dvtt = Inf;
7 A9 f. s& [6 \$ S t = n;
  M$ S. p1 C! _1 l2 J3 H5 x while(1)
. z5 X5 N1 ~* v1 T+ \7 {, s       if a(s(t),t) > 0
- g; A5 l: X7 C5 y" L          %前向弧调整量1 c' m3 Q/ y4 D: Y! U
         dvtt = C(s(t),t)-f(s(t),t);  m# Y2 u  N  i# k; B
         %后向弧调整量
5 b' P/ m$ b' `$ |& g! {: x       elseif a(s(t),t) < 0
8 @( {9 ~/ {8 f, a          dvtt = f(t,s(t));
1 j2 {6 J7 U+ i& ?8 v% L: x& c2 t       end0 L/ U( W& z8 O9 [& a
      if dvt > dvtt
8 q) N* D3 r/ u, C+ Q          dvt = dvtt;2 ^. J. X: W9 b) [& H5 F# g4 A* @
      end9 |; ?$ {1 {) E  S0 d
      %当t的标号为vs时，终止计算调整量
6 e' u& j5 T" O  R8 L" u$ G       if s(t) == 1# l) q% I  y* l3 K  J% }
         break;9 B: z6 L1 L8 K& l. W; Z1 E  l
      end
( e& }$ g0 X+ d' t' T9 g2 h. a       %继续调整前一段弧上的流f- s: i6 m. O$ r
      t = s(t);6 e4 ]- c* P/ w, Y9 C
end- ~, S" F8 g! f( j8 n
pd = 0;
9 d$ c) W: y7 e" i  S, h% [; X %如果最大流量大于或等于预定的流量值
8 t+ g5 ~# w4 l5 a8 T3 F if wf + dvt >= wf0  L3 [7 Z6 v  d1 ^: C
      dvt = wf0 - wf;
/ B; a3 p7 H! i% K- _$ Q7 f       pd = 1;
* g7 D5 S6 O. |+ {+ [ end
. N7 P" h( B5 F t = n;
- ~# v/ x+ {. y* v$ p3 |/ u' Z( ~ %调整过程  |+ [& |  @7 c1 o4 ?
while(1)
0 M& F( w2 v: A* A: I if a(s(t),t) > 0  m% q$ p  q9 R: {! q
      %前向弧调整
% L* p: ^. W/ @" w/ ~5 ?$ Y) C       f(s(t),t) = f(s(t),t) + dvt;
$ ^: {1 J" q8 W! B. w/ s! o2 ?& ` elseif a(s(t),t) < 0
' U: p4 W* G8 O3 v- @; c  a' W9 U       %后向弧调整
' d" d$ ^2 ?7 ?4 Q+ }       f(t,s(t)) = f(t,s(t)) - dvt;
0 G6 d; p% Q3 P end
" w2 V- h+ e  {* a %当t的标号为vs时，终止调整过程) x( |! s( t+ E: Z$ z) L
if s(t) == 1
      break;" ]1 G! n/ ~: G9 B
end
9 c3 G$ \8 |( K- x- \ t = s(t);1 i! ?+ _) \$ C9 o& @! {, m
end1 D, q7 {, ~& l5 s
%如果最大流量达到预定的流量值
) h( ^5 B! i6 h, r if pd
0 r# j/ |; y* t' z$ H& g7 R$ M       break;% t# W$ t6 d( I/ b1 ^) r* M1 o1 V
end
% @3 k* L" j) e) Q %计算最大流量
6 f! H1 P- k" p9 I" ?- X wf = 0;: _* Z/ v+ y: W2 x7 F* ~5 o. ?  X
for j = 1:n
. J0 V/ H% ]4 g  j0 |1 g% J       wf = wf + f(1,j);
1 C% J0 _: L7 \# P6 O end" }6 g. `; x: `  U  X5 i5 N# k
end
  e9 e0 k! v6 J. `" }9 {%计算最小费用
& h( Y2 ?3 N3 M: n3 D2 ^zwf = 0;
% B4 J2 C0 {/ R2 Z+ Cfor i = 1:n& N& u6 G9 \' T; M0 G0 o
for j = 1:n
+ r2 o/ [! v; Z       zwf = zwf + b(i,j)*f(i,j);
5 V, z- B& D0 L5 \ end" \: }$ J0 y: X% H5 Q/ C
end
$ r  C$ A) e8 J* a6 d# f% \7 j%最小费用最大流" _& a/ i6 {, m: {8 x6 F
f; R4 _, _$ @# L5 s( E& C
%最小费用最大流量, c; j# C; j/ ?- ~8 a6 K5 f
wf
# [0 l4 F5 g4 i1 }+ i* z%显示最小费用% A3 a  k- L% }& J0 y; u# ?$ q6 C5 n
zwf
) y/ I' Z' g' Z5 M# ~0 m7 a. M1
5 J* h. _" X) h7 L& T2
( F7 u9 v0 R  b2 i3
  o& ^" y. O8 S! Z  W4
' M" M9 L# w' H) L$ j1 t" x9 t$ z6 _6 {5
. F" {9 d& h' B7 q. c6
% |7 u( @3 P/ O+ o7
+ ]% x3 m& V$ V8
: w4 U* E, Y. i, Y1 o9
0 D4 c8 @, [& \9 Y7 K8 o10+ W7 ?% Z$ T" I
11
/ ^4 X) q2 j; o2 ]12
: \& }# H; P$ A& _) n+ l. l13
/ m. D- e6 k' x14
- Z. A4 ~- ~* j15
7 y& e& a6 t: o: `6 h) o16
& I! X. `+ z% n; z  E- c) ^17: P. Z8 _/ S/ Q: g% v
18* c1 B3 G- x- ~9 r
19
  e) t3 e5 k9 B" e( w( |/ a4 b. g; `206 Z% ~0 `$ a: z. W8 w; J
21
4 Q7 o1 Z' O) S, q22
+ C0 e) l+ M7 U- W" |23* ?4 H; K. v! E& V' S6 T8 h
247 _: Z2 o9 }/ e& g
25" o+ a% F  Z) G
268 m& _5 Q3 C# x. S% Q: K
27
& b9 b) f8 b3 H! H, D3 G7 V" I28# w0 t  W& d# u& y
292 f; o7 ^5 V+ R6 m2 w: k
302 Y, i2 k- k$ N
317 F& N7 V1 K) r$ N' |
32
$ ~! B1 p4 A! O( H& Q# J33
2 Y5 [  i; f9 v) b34( r* _) S1 d# J" p& Q
35
  y( b' M1 }2 m% b5 ^366 x$ K& {' t, i
37
+ ^6 N+ k# M6 [! i( K9 l38
* `  X0 I. ^) A' J+ m" e% O4 D9 h39# B$ i8 T, X) M& @2 i- }
40  q- y5 j$ [' i5 N2 \
41
/ G' t/ K; X% g# z, f. m, P425 r) W8 u  d$ }; a5 ]# k; N
43
* R3 A! o  s- p$ F5 R44
$ u9 t: ?% E4 o9 m& ~$ `; S9 |. T458 ^; _' j4 c5 y8 f) `2 O
46
! v5 ~3 h8 C5 i, r47
0 w/ O0 Y- I8 M3 d' y+ i4 t, z2 w48" Q  M9 w: j4 X0 |
492 C( ^- d$ Z/ W7 i% L+ P
50
& v* R/ m) x9 D% s& x+ J% D9 m51' m$ {& m* T5 ]0 h: C" N8 o6 i
52( v$ n) p, v4 `: C9 K$ I
53
( L1 G' F1 t" a+ R/ r! m547 k% u! J  W- |7 b
55
( e) @4 E8 V' E# L56) H% v/ ?; e6 e
57
, ~( c1 n7 F  h( ]) C! T58
# }+ J& a; N' o" ?59
0 O4 A/ f( Q9 m+ D: t- r0 B605 |, _2 T$ Z& w. T. r
61
. U7 B* Z7 c& G624 x  A& _7 T% K- ^* a
63
6 @4 p7 n5 c, y4 s$ g64" ~' A0 S: E# b7 e  ?& G
65
' G3 B. @6 I  _9 J% M) v5 w, @66
. T! B4 l8 ?0 m  l4 N67" N# {3 H4 A3 w$ W! C
68( V  @. X) w# c. S0 b
69
0 G% o2 s8 h0 p  v6 ]70$ ]! m* t' v* o$ M0 V( u( Z' C7 ^
712 s1 x9 i9 K6 s& e" h$ Q
72
6 Z; v+ ^! }9 w; \( P  y9 ]73
1 H0 v$ z  H' J/ w. o746 o% E4 H3 ~6 t7 ^- [
75
$ P/ W$ R3 h6 D# f, [- O76# Z" i! ]/ j% `2 ~2 m+ E
77
  I. z% ~6 C! }5 A78
791 L0 k- N, @  U; S
802 X+ `4 P$ c; b# b
81
$ y5 t  M4 s! ~* [9 Q6 Q7 q. E82- l6 y, {" J) s$ A& V% E
83
; S3 I; Q+ P, e" C3 G' L) G0 ~: D84
- j; ^( Z5 ]. A5 _3 }85
$ k& D+ m) m+ w; j9 ~- R86
  N9 n  o5 o2 `$ U8 j) L1 M873 n( h$ h& f, l3 Q9 f
881 s8 @& n- W! n' G% M6 [& s  `3 ?
89
3 T  F" E+ w4 [/ P90; T; B0 _' S/ G# |% u/ x& u1 O
91
. c- @' ?3 @$ ~. e92
- P  J0 i1 f3 w0 |) f6 Y93) k( w* _* }1 @: D* j
94- l- y+ y2 G4 w
95
9 F7 ]& I% u" O$ U& S' N0 f# B96, B# m7 R( Y* q) u. @
971 E+ O5 o& b& a( }- }/ v' x3 z
98
  c: {) b! ]" }& T) ^99
* R/ J( F9 r1 F100
# A: A9 G' R: H! |  x  w1017 p5 z% t9 f' N
102
" h, b0 _( z( u% Y3 K; N103
( M8 q* l* [( K104% D9 a9 H! s' D3 R
105- N2 f% N% v# R3 Y9 {3 V/ z$ t( }1 g
106# n7 v& X$ Y0 g5 a$ m
107" z0 \+ q8 ]  X/ J
108
0 n; E- O& g! ]/ D/ |1091 n9 @  n. G; ]0 M& B6 j+ X6 F
1102 J3 a% h# h/ s
111( D" C$ C0 g8 ]: d1 D
112' w% j5 _6 [2 c# g1 x8 _
1134 r: Y# K$ }1 w0 ?& y: Q
114' o6 c; q) F2 M4 P+ M  w3 D' x
115
! N  Y9 v* q0 d4 e( Z5 z5 d116
# K9 I& {, F! P1 W: [3 t. y7 B1178 y: U4 p6 o% D# U
118
( C" q; b7 C* X$ b119: V  |( l% n3 ?& t3 K, q
120& ~! q8 d  J- p9 I- Y0 a
121
5 ?2 k0 }5 H0 C6 f: y# V122
: a7 m9 u9 T4 U( d123% M( d% Y8 t& _# E" G3 Q
124
7 z+ I  _2 g( O2 e- G" p# A125' O! A0 e! y0 b1 `" b; P
1269 s/ @' B# R+ I5 F
127  W7 I8 x* J, E5 l
128
5 V( t) L. M5 Q; Q+ ^5 @129+ H1 q" Y3 u# n8 f
130
131: z; r9 K* g% A2 p' o6 P
132
( J/ Q. N% n$ Y$ B133
3 N- g; y; r- J1347 j: u7 P' m2 ]! I
1353 X5 a+ [" a; p1 F" q
136
. P0 J8 @. l4 K3 z- B- L5 Q137
% B4 n$ [+ [4 [' m4 I7 M138
" v, o( B  a% a$ A* V  j3 T139; x% `# n7 M1 C" b6 j) A  z+ R
140
0 s1 ]0 w9 C$ |6 {: t2 j匹配问题：/ Y: G9 z- f2 A( B- Z7 h
- r9 J) A' c' c0 i" ~( @
最大匹配：3 T$ H8 T6 s! Z/ L

% l9 D4 i. z# s0 V6 g2 d: D2 q
代码实现：

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