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

最小费用求解
$ p) h3 W! j+ @4 v
' o/ C/ R! s2 R7 F2 h- u, ZLingo:. a% ^% }! R& q; m4 h
1 h2 o" \# q+ k1 ~. ~) i
model:+ k: a0 Z' p2 g* ~
sets:
8 W/ W7 t4 D8 ^3 z8 y5 [nodes/s,1,2,3,t/:d;
( p( D" d8 U# T9 s* ^+ L' Darcs(nodes,nodes)/s 1,s 2,1 3,1 t,2 1,2 3,3 t/:b,c,f;
& a# F" M, ]& C: m  x7 u: dendsets9 y0 i/ v+ o# Z9 k# ?) s
data:
! ], Q% a/ E- p( i8 J3 eb = 4 1 6 1 2 3 2;5 H2 K6 ]- q) i
c = 10 8 2 7 5 10 4;
4 }6 q) L* J( p1 I+ Fd = 11 0 0 0 -11;
8 {9 ^) T- [4 d: o# n# ^enddata1 p5 k9 o' Z" V1 l* _* J0 g
n = @size(nodes);7 e3 d( ^! O- l+ E% w& F
min = @sum(arcs:b*f);
. y; [8 [' R: i! K' D: G  w! B@for(nodes(i)sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i)) = d(i));5 u  `9 M9 i8 x# r
@for(arcsbnd(0,f,c));) X% R* `2 B, t
end; o$ q8 n# p0 H
1# _, s) T1 K$ `/ y% F
2
- {2 d) a5 Q9 H3 z3
% K1 B" j! v8 N- o8 G$ g; w4
# q; l2 b* v  L! \; g5# m5 s5 S$ A; U3 w8 s& Z
6. ]3 i6 r3 h7 C: X: g
7
" ^( ~% f6 v4 d# Y8$ Y- H  n4 }$ S3 D3 S
9
' e9 H7 d; }6 X: b/ W. ~, v10) s- n6 v4 f7 L# `5 h" G
11& X/ }9 e: B* N- E( f$ h8 N& ]
12
; i3 U1 V; R3 s5 R13
- y6 @# F- H6 T14' x  ~8 _( y; k0 H6 |( {- ?- r
154 Q  G$ K% E4 f2 P% z/ Y
Matlab实现：; |( D' x; [* h7 ^+ V' ]! V
  J) w& z* k: F9 W, V

1 ^7 W& M* O% R* q( R' un = 5;
! o4 ?7 M1 m, r. q2 W* G%弧容量
5 ^: H; O- u9 Y5 d6 Xa = zeros(5);
( l- o% T6 }+ O6 [! |5 ya(1,[2 3]) = [10 8];
8 n2 L; w9 s1 ^. b: Ma(2,[4,5]) = [2 7];" s/ C! w- A- }& ^+ m
a(3,[2 4]) = [5 10];. t1 y% S  h; j7 R: p
a(4,5) = 4;
: z* n/ ~+ B& p& DC = a;
' {: B( a& M( Q1 ^" P%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];+ @( O  b4 ]$ ~( Y! ^; F
%弧上单元的费用
3 c$ s, ?$ I7 A" L( ea(1,[2 3]) = [4 1];0 k: S4 J# v4 N. a/ }$ e% b
a(2,[4,5]) = [6 1];6 u& D) F' n4 ^3 M3 a9 Y
a(3,[2 4]) = [2 3];# A2 m; }/ K( i* _$ _% ^
a(4,5) = 2;
% T/ _1 d  B" _9 }# a  Wb = a;6 h8 A6 d6 x- ]+ j2 d: O
%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' C9 o3 v$ C3 k- \
%wf表示最大流量。wf0表示预定的流量值
7 p( g0 U/ m0 m! |: z( }7 k7 p6 ], hwf = 0;
7 A, V) t1 |0 o! i+ owf0 = Inf;
5 \" l' @& h4 k! z% L%取初始可行流f为零流
/ W3 N7 e: R1 I! t  |2 U& qfor i = 1:n
5 W$ r7 W5 U; B6 ^, _- P for j = 1:n0 l( I; k6 y9 y' z- O
      f(i,j) = 0;" O8 s( [. h$ B! j+ @
end6 x( w8 M" I' t' x
end
& d) ^1 b" z, x! }6 h( Y2 \while (1)
" \) L# h0 H/ D6 A! V9 f %构造有向赋权图0 F4 o# W+ ]. x0 z/ b" J+ Z* D# \
for i = 1:n
# E3 c$ U/ m  L6 M5 _       for j = 1:n
5 ?; M9 D3 i( k7 S3 [          if j~=i1 k$ g5 a( a1 _3 g  \
            a(i,j) = Inf;- I4 m+ q3 J" ]- j% |9 G1 p
         end
! f+ ]; {4 r5 F/ K       end$ |1 E$ [+ W' K$ C1 ^
end
; q6 h3 w$ r! r' `) d for i = 1:n! ?) H$ W8 k8 y, Z: `
      for j = 1:n
! [) E: Z6 \# E0 L          if C(i,j) > 0 && f(i,j) == 0$ _, F" C! w; D
            a(i,j) = b(i,j);5 A; L8 t( ~7 B( ~$ r' T. {1 w
         elseif C(i,j) > 0 && f(i,j) == C(i,j)
; p1 d! E1 [) {" G# t             a(j,i) = -b(i,j);" ?3 t2 e' {# E# R
         elseif C(i,j) > 0. b0 [$ z# Q: e( V4 ]: B
            a(i,j) = b(i,j);
/ ?( Y% K: q% a: [; g             a(j,i) = -b(i,j);
3 p) \( J9 l2 ?          end5 S# T9 e' P& N. ^
      end
6 p6 m- w  q( i4 F+ ]2 ] end
* v6 I- P: Z* d4 i& V: Y %使用ford算法求最短路，赋初值7 t5 v( {* Y% z& v
for i = 2:n
; }4 \% O2 q6 Q& w8 B% U3 ^8 T; m- k, W       p(i) = Inf;
/ l5 q$ o. m$ g7 Z! E+ M       s(i) = i;$ O5 P: \  t, `9 U
end
$ |. K. Y9 Y5 O %求有向赋权图vs到vt的最短路,赋初值
, D$ C# s! a& S6 Y2 w7 a  u for k = 1:n
8 @/ x( {2 c  l/ d2 n       pd = 1;; X+ b+ K, v. n8 A1 i( L8 {6 d9 s* p
      for i = 2:n
! p2 C% a/ \+ n8 J. N( L+ u          for j = 1:n
  U4 i/ A' i$ o+ t             if p(i) > p(j) + a(j,i)( Y" m3 [$ u7 ?) e" q
                  p(i) = p(j) + a(j,i);
4 u$ F9 F" t& V0 A" n3 {                   s(i) = j;
( f; M6 G" q: J' b3 x6 q4 ]                   pd = 0;/ f* o$ R* H* n+ O  d- K" r, O
            end0 m( a2 R8 d: ?
         end
" v+ D0 [1 H# A' t       end
& {& {$ \/ h" h" J/ x" e% b& T3 H       %求最短路的Ford算法结束
4 k( Z1 X( N; u* u       if pd
' {3 f) _1 v, `* ]          break;
3 @2 J. w+ E4 E% l       end
$ Q6 @9 P# J- |/ F3 g4 o end
0 P9 l) H' q9 A4 D+ O- B0 ? %不存在vs到vt的最短路，算法终止，注意在求最小费用最大流时构造有向赋权图中不会含有负权回路，不会出现k = n
. C4 f/ p9 U5 v* x6 J& M  f if p(n) == Inf) F- t" `/ H, @4 @! B' X, M
      break;: @: x1 K- X4 R3 i6 M: H3 y9 o
end
7 O7 a2 H1 R% W; d %进入调整过程,dvt表示调整量: H+ i) u& A) U  [3 i) a( b
dvt = Inf;! o- {/ A8 c+ e% [( W, k- ~
dvtt = Inf;
+ Z2 i" f( w& N4 I/ @2 X4 R t = n;2 {" f/ Z$ K- M3 t# X
while(1)
7 B0 k( Y! q' [7 Z' P  h2 Q       if a(s(t),t) > 0
! ]" \2 s, C4 @9 M7 H: `- g          %前向弧调整量6 W2 _6 N4 u' f: L8 T
         dvtt = C(s(t),t)-f(s(t),t);9 d4 B0 Y: y$ c: |  O8 C2 p9 n
         %后向弧调整量# B$ h% g  M' V: O2 ]
      elseif a(s(t),t) < 0
. Y' B9 s+ o* i          dvtt = f(t,s(t));
% v1 {' n" P$ p- A) P0 [2 B& Y       end
" A* t' k  O# `$ b* Q1 X: j       if dvt > dvtt3 O3 j% k# c" `( J' i
         dvt = dvtt;
' \- N" w. l/ o0 @: F: S! k& G       end
, P4 ~: U8 z# z, V       %当t的标号为vs时，终止计算调整量8 V$ B% Q! x5 Y+ L  A1 a0 R
      if s(t) == 1
) S" r' ~& E7 x4 }' ~/ X! h          break;: g& s, r' W' |& Y
      end
9 F  a+ T6 F8 @$ y1 F$ y0 N; M/ V       %继续调整前一段弧上的流f" b! ?8 }, F$ k6 T% i: Y& S" S
      t = s(t);$ K( d" M9 s0 L1 b
end
6 K9 W$ ]6 r/ s9 ?. n0 S5 T pd = 0;
5 {- C9 o1 d) |$ F! H" k) Y %如果最大流量大于或等于预定的流量值
5 h5 K% W$ I1 F& J if wf + dvt >= wf0
1 E" `% N3 S) H8 N/ x* X* \       dvt = wf0 - wf;0 Y4 l2 _* |! l
      pd = 1;
* J! b2 I% V: O/ x, f- S. S end* Z& X$ t4 G7 ^% R
t = n;4 U9 o0 z+ [% u' e6 C, c/ E" Y
%调整过程
9 j3 X: `+ C7 v while(1): B- ^- c& W; x8 j' C/ A* M
if a(s(t),t) > 0
. @7 [6 N) F* R5 Q7 V3 H       %前向弧调整- K/ O' u; W8 N1 x& I- ]
      f(s(t),t) = f(s(t),t) + dvt; & K9 W3 I0 K1 M
elseif a(s(t),t) < 0: t$ g- V$ H  [: p* q
      %后向弧调整# P9 G3 G+ p; C5 X  H; q! S  {$ v
      f(t,s(t)) = f(t,s(t)) - dvt;
! p) r  ]2 E' ^ end % T' w9 y, u- B. u5 K( S8 W
%当t的标号为vs时，终止调整过程
* F$ t  y% S+ K6 i# H if s(t) == 1
8 |, p* Y0 O3 M. m       break;- `5 m' c# G+ I6 |5 m5 H
end
6 d0 c: ]3 k$ D; ]8 ?# c t = s(t);# _8 Z2 `6 h* h
end
! E6 A" L; B) Z% {- ^8 K, ?& k %如果最大流量达到预定的流量值
+ ?" r: X. o1 [) `( w" B( z if pd% Q2 G) s; P) B) G1 Z. D
      break;
  j. r* ^7 z' ]* u end8 V& N4 V; m$ U! u8 e3 e6 o) p
%计算最大流量' W1 R7 f$ R: B- R0 H, U) e; P
wf = 0;
# L: V* }0 `1 X) |: _3 q5 [4 m for j = 1:n
* H; i; e( g. J8 S% \. y       wf = wf + f(1,j);# o/ ?$ e# O" t1 c" E7 X
end
# t, Z! D0 _2 Iend8 M+ t7 I9 \; u8 X" |
%计算最小费用# I1 \" U* x8 K" ^8 j0 P# I
zwf = 0;
/ J# X5 F9 _0 T' z  j: O% z# rfor i = 1:n+ M$ u' B7 ]4 w9 g1 m3 [% p
for j = 1:n- E+ a8 ?: i" o( I/ z. N+ r
      zwf = zwf + b(i,j)*f(i,j);- K' d8 r+ i( o4 R4 }5 o; t
end
1 D& |9 J( w0 F# R4 {* e& z- Fend. ~" n/ T. Q1 S& W/ y# O* X# R, f
%最小费用最大流- e: S. U5 d- V7 T
f# s3 y3 R+ r+ L+ z( e) \. r5 |
%最小费用最大流量
; J5 R% n3 p9 Owf; P9 l- I8 H8 y  y
%显示最小费用
3 t2 E; P  |; L: ?& t' o- s8 c7 Wzwf
7 C8 D- u& ]8 l% _3 [7 k: A1; d$ j9 t9 @% _: j
2
- g! P! C+ j$ v5 Z* i$ v& I34 D  R. G! `) Q5 D. o
4
7 m1 ?" d7 |& _/ W1 f, D5
& i, N5 @% P! Z$ g0 X) p5 o6# P- R. w* H4 D" ~
7+ d0 i+ C: k  f/ X5 A0 t/ ~1 M! [
8
$ U* e# X* P3 F* R! `8 ^- p$ H! c9
' @/ w# ]* |% k0 e6 {' Z- I10
/ H; @, M/ B# a- z6 n9 S116 G- c$ w$ G( a- [0 j9 @* U. P; O
12
: g. d7 Q% p0 }3 m8 A13
; @9 B+ l) Y( ~% Z14
/ t2 [% [0 x- k% M15/ q2 l: S# W% \
16, `1 O- W& q5 T( O# C+ y
17: g# y5 ?& N7 S3 T8 l, F' [
182 |' e$ x  C' D! ~6 l& k; H
19
) q4 n: g7 P3 T' {) V3 S9 k. i208 l+ K4 f! w$ C/ Y% ]/ P, A* K! o
21
4 l7 L7 c. L6 j2 f0 ~" D# u3 b& m, i22
" g, S, P' X( O0 }  @; f237 n0 G. `- w% C. W7 Z7 R  f1 h- h
243 a7 C9 M7 v) C% m2 X
252 u5 u' v0 W. Q9 r! R& G- h# z
26
, u$ b1 v/ r( x3 G# O, h7 W278 g* A$ X0 k% }0 h# j
28
* l; _1 }+ v. U1 l294 h% Z( c0 {" h' N" W: @9 e1 D: P+ N
306 B$ V7 o! I3 w0 t7 a7 ]) Q
314 n& l2 S+ c: [) f/ `9 \7 S
32
1 |/ u+ u& G, X2 @337 P1 f/ q: y- }. J) Y2 \
34
  k6 `% e4 p- g9 z8 x35
$ {; v0 t: b! D' q36
1 E' k$ d5 m( w1 d37, ~( l* g5 K2 u3 N
38
% p* C6 G3 J5 Y4 T* Z- T9 ~: x391 B* j' R" H  M/ E+ K  e& n( h- j8 ^
40( i3 k7 u8 P1 E! C+ R/ ]1 _& t4 U6 z
41
5 u+ t2 _7 A; X( S42) ?0 q3 Y& W5 J* w6 L# L2 y/ Y( |
43& X) T. J  Y0 k/ V) j
44
8 Z1 [) N# R  E* s5 ^! Y# t45
8 N0 @5 U( M* j9 ]  }46, a, l3 M0 p/ _) ?' K8 o( T
47
/ n  {: j, U: M# g% g484 @. m, ]. D: l" n! w4 O  B
49' _( W& \% c- v( j
50# U$ k  T: F9 O# U
51
! f* {3 h) |! w52
8 L' x2 L  s, `53) _8 x& V  L2 z: X) p0 w  @
549 F7 p: B% Q+ Z" O
55
( ~) V0 M, J0 V  e" U560 X  C4 _/ G9 U3 r- ^2 B
57
8 l8 F- m% m  X5 x0 w: z# _58" r5 J7 ^7 k! Y3 x/ U
59
, s/ u3 \3 p" m6 x0 e60
) h! w7 L$ J1 s0 o61
+ r& G5 |8 W) x. c  L62
' r4 Q; ^; L; P0 N* \6 f' N63
1 n/ d( X) n+ z8 i64
6 E/ T; r1 }5 j* O0 x. g1 r1 j65
* C* \3 G) L9 j3 b6 B  S66! l7 g8 \$ u0 U) d, M/ K; G+ g( |
67
* W/ P" c1 P  h  a# `8 B68
" r+ S  X- g! ~$ V9 Q3 f1 i4 M694 G6 H% I+ }& [" g
70( i& A% }$ e  E8 P
71* k. }' O; Z' R  _& C5 F9 [
72
( I0 J# F  [" B9 G+ B73, ~: j5 Y, x+ T: I, D0 c, `
74
; O! M1 B( b1 @- i75
. J% N4 \  C( j/ @* m76
) L# }! w9 U/ i: F1 y- F6 a; t77
2 d6 Z( ^# @0 _! i) V* R/ e78
, j+ @+ ?, ~8 o79$ Q) }; J! ^4 z3 {  [! A
80
4 G2 K! Z% j0 ~5 N) z0 _81" Z1 U5 Y, b! Y/ @6 n
821 D* f7 I/ \3 |$ w6 F
83! ]( B: y6 a/ u2 H0 r
84
* ?- x& F& z! \* h85
8 m6 T5 g+ |% D! }4 d7 v; H+ W( Z864 D) ^0 B: O, [! Z( G& X+ y5 p3 H
87% r0 T3 q( w$ D' T7 z
885 V* F3 d9 F: O; W( E. G7 a
89/ C: ~4 y1 \  I$ }5 `
90
5 P1 J$ N4 Y# H% [91
/ `9 g/ P$ }& O2 Z# d- {$ S92) ]* e# [7 W3 q9 v6 t7 T$ |8 W
93
# R/ v" C% Q1 r  U$ `4 ^94
& b9 S' D( f+ C6 o1 S95* [/ o6 e/ G  a; n9 D! o% G1 N
968 L- @( o; z8 V! R% d5 ]8 z6 }
97( j& D9 X& c7 S+ m% I. o5 _
988 I) m; J/ r- x$ z' U6 ~
99
& J/ I) }; q) t6 v100! W6 t. T/ I3 b8 _
101
  V% _) v: a' ~6 V5 U; b% e102
. _0 G/ ^' u; U7 l  h: `1035 D4 x+ T: b/ N9 e
1045 [- e& N3 Y4 q* @+ V7 r" d
105( p  m' X9 |( `4 S, v
106; v3 S7 C' v$ F7 H2 `! k' k
107
3 N1 H* Y3 ?2 B& M/ `/ N: c) D108+ N, r3 u+ n& t
1093 ^$ i0 ^' k! c  t- {1 B7 K9 W  @8 A
1104 p+ l) n) m+ F: ~& e
111
/ q& K+ R8 R8 ?. P9 C* i112
8 [* x: r; {3 P9 v, ~4 G9 {113; ^6 r/ h3 J- g* q0 U
114# \! W) q" a  v+ B6 ]
1159 @+ s/ u( x* P& T; m
1167 T4 Y/ O0 f2 V+ q
117  x1 ?( I* c& X  Q  I& f
118" V6 y3 G- J8 u% Q! Y2 t: V* ^- c" W
119
  _$ k, s# E2 I1204 ~( I' Z. Y% h. H3 X
121
0 L; v# Z& A+ `7 ~122
4 T5 t6 F8 H" Z+ q. j123+ d0 ^. @# o1 K+ ?3 V" m' z; B- s3 ?
124* c/ k* E% M5 G4 x% T$ E% }
125
4 ^' {7 w3 k& v6 \$ v: h126
' ^- X' C  T1 h, z. F7 ~127
6 \7 H# J, S+ r+ C0 F128* k7 p+ R) h# E5 j
129
4 G9 [$ J4 {7 p/ [( s, s130
: D' A1 V0 D$ V1 g( S$ \* a131
. T3 F2 x) M5 y( L132* \% l8 L6 C3 W& ^* p
133
6 n  d3 O$ V/ I1 m1 Y134
4 }/ C0 A) i3 v135
7 ]5 s* m( m' b6 m136
- b, l$ ?1 r2 a( r: y) [- H: h1 p3 [137
6 @9 i5 f. V+ V, L+ K( |2 F138
, K+ |/ H3 A! @: m" k9 D139
6 C) {! {2 \7 k0 T( A, O" V140
, ]9 O' D3 P$ M: {7 N5 {. }# i匹配问题：
: g/ E1 z2 j8 L) A9 l9 s0 x; b. c/ S0 H5 n5 k
最大匹配：
0 L9 e7 F' z1 G, k5 d3 C8 x+ U% y! e' |

代码实现：

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