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

最小费用求解
8 K/ Y4 Z" ]- l/ l$ t3 n, N) {0 Y+ M! p
Lingo:
+ d+ ?- t9 c  Y# j& l6 N) Q! k! X9 n4 C
model:
$ U' }: g) {9 l' k7 Bsets:
" v6 L- }9 m2 Mnodes/s,1,2,3,t/:d;$ H: w4 @+ y6 r; E
arcs(nodes,nodes)/s 1,s 2,1 3,1 t,2 1,2 3,3 t/:b,c,f;
$ B) \* |5 h  s" B& Tendsets$ k+ O( H1 I, M$ m4 W& E) _
data:7 @  I. f& ?3 S, b! ?5 o, A
b = 4 1 6 1 2 3 2;
" @* q' }# F5 B3 U9 }/ H+ kc = 10 8 2 7 5 10 4;
* S# J5 h' r5 G. S: T) @( Q# ~d = 11 0 0 0 -11;. b6 M+ d, \( h5 x  {* {
enddata$ q' R1 ?' C  s* \/ h7 b
n = @size(nodes);# K  u# [/ }3 h7 [" T5 ~5 F1 ?' t
min = @sum(arcs:b*f);: }7 d5 L( R8 h- b# T" F
@for(nodes(i)sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i)) = d(i));
( B. H3 Q7 ?2 O' n' e9 J@for(arcsbnd(0,f,c));  J/ c; _. B3 W% l& W, X7 w. S7 E
end( l3 l+ l8 I& F0 s0 N$ g5 p
1& m+ x6 k7 c# J( L( M5 t* L; P
2/ ~+ z& O/ K1 ]9 q: Q0 }
3
' R6 J% `" o9 C( S4
- V7 q4 l4 W1 u0 }4 l  u& k8 Q; [* c5/ F  I. P; \8 A4 F& K: M# e$ W* S4 g
6
' X+ x; ~2 T5 g7 ^! f+ I3 X7
/ ?% e( W9 j! w1 G5 n5 ~7 g# Q8
0 f9 _0 Y! L9 ?6 o/ q, O+ b0 f( `9
- K  E! b2 J: ]4 z, U10
3 I' H1 [: }' q/ d! W115 R) n9 F' R# @
12" ^5 X2 K# w2 N' a3 C/ a
13: H8 Z. p. t9 l) |) h8 j, p
14
' w8 H0 Q' C) f- j15
% N* p  t3 B- N5 ?6 q" U( Z1 VMatlab实现：
+ c8 P; i: ]$ S% q
  w" K/ ~  Z* G- W5 f5 {1 S( N. I
n = 5;
8 @5 v& n; G& B8 G%弧容量5 j7 w2 Z# J- Z) [. S2 K
a = zeros(5);
( Z& C  b# X8 f% \9 Ya(1,[2 3]) = [10 8];" @3 T1 C2 D" A' c- ^+ W
a(2,[4,5]) = [2 7];
/ q/ L, s; v* u. I4 w; y7 _a(3,[2 4]) = [5 10];
& c1 ~" E* L6 j: V* C" @- ~  xa(4,5) = 4;4 ~0 }0 P- }5 ^5 v, v; {
C = a;
, r( G& o; v( l+ a/ T9 B%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];- ?5 @' }8 v0 ?$ ?8 ?9 P: Q
%弧上单元的费用
( D) G/ I3 s$ f0 y$ ba(1,[2 3]) = [4 1];4 _, p5 ~  ]" |" j
a(2,[4,5]) = [6 1];
$ ?$ n: i& G9 g% J0 oa(3,[2 4]) = [2 3];, n" L2 N$ \$ q7 s! X! W
a(4,5) = 2;' E- C3 y1 ^4 X7 j; m7 \) b& t
b = a;
/ S- G3 R' P, D/ 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];
1 v9 w' ]% d- q%wf表示最大流量。wf0表示预定的流量值+ v) B/ w. @$ D6 r& M
wf = 0;
, q7 S7 }( q. Y/ D. q9 W+ @" t6 Owf0 = Inf;6 a* P7 ^  x8 u/ C7 A  M9 @7 L
%取初始可行流f为零流
# U' r: g% o( B& t! ofor i = 1:n1 {: \$ H' i/ r3 V3 l! Q
for j = 1:n- ~3 I6 S  w% ~( @
      f(i,j) = 0;
. _9 Z+ @& i) i$ [; V end
* @/ b0 `/ x' B/ Vend
" \/ X9 M% K9 e- ^& y( B4 n. ]1 qwhile (1)
/ t3 i: C7 y* P2 n# Z. i' \& O2 c %构造有向赋权图, f" X, p. e- q4 C/ S1 Q
for i = 1:n
3 Q: o2 O! H' A+ `8 _. x$ m       for j = 1:n
% R! _, I& }- j2 m8 A- Q! p          if j~=i+ U. A. u7 }& P3 t- ]1 ?
            a(i,j) = Inf;9 B$ b$ z, B  n, u. n- [
         end, G  g' G# u% S5 Y
      end
7 W' W8 X2 c+ T; ]1 i9 m9 v) m$ ~% z end  ^# p8 Z6 o/ b) E6 V
for i = 1:n
* d. C1 s# `! b. p3 [       for j = 1:n
) l7 A/ f3 f4 N5 R1 C) u          if C(i,j) > 0 && f(i,j) == 0
% V$ J1 D2 f. l& v             a(i,j) = b(i,j);- H; }( Q7 Y! p3 e+ }# U, S
         elseif C(i,j) > 0 && f(i,j) == C(i,j)
$ H7 a3 C* Y, A3 B             a(j,i) = -b(i,j);
1 S% J( y& Q. {/ z( H+ a2 s          elseif C(i,j) > 0! s' ~  r. A, t! k% a
            a(i,j) = b(i,j);8 q5 ^7 V2 {  H! G& _
            a(j,i) = -b(i,j);
8 I0 z! _, S* {8 G: `          end. q% ~  |3 ~' p. L' _& v) m# ^: |
      end3 a+ Z$ t4 T* K) l: v' e: V
end
7 G2 q: w$ Y+ f4 B %使用ford算法求最短路，赋初值
5 L6 i, g2 ~) y for i = 2:n
      p(i) = Inf;  [- i+ Z. w! q$ G8 m
      s(i) = i;" a. X% P3 A$ W- g; P2 M- R
end
9 \% k7 u1 |% y" { %求有向赋权图vs到vt的最短路,赋初值
2 {6 }7 u6 C2 |, @1 Y for k = 1:n  P. K4 P+ ^. M+ U$ u
      pd = 1;6 Z/ p2 t3 b9 [. `5 }
      for i = 2:n& e) Y3 L- h) r! `
         for j = 1:n7 ?0 P, s5 C3 A+ D6 q
            if p(i) > p(j) + a(j,i)3 x( Y& T9 r: s2 w0 j% o
                  p(i) = p(j) + a(j,i);
. y- H2 L3 v$ n                   s(i) = j;
/ h* {. ?* T- j% |; y- F, Y6 a                   pd = 0;: E, s# a! m* i9 {% B
            end; L9 S. e: j3 H* q" G$ i& Q  n
         end! f! q/ Y4 m% P# ~7 i. B5 w" _
      end3 C$ I3 Z9 P6 R9 u
      %求最短路的Ford算法结束
- ^- n5 S+ ]# v4 D9 L! i# l& S       if pd
. _0 w, \4 K9 x: w. D' u* @+ O* @          break;
. t) e: h* e5 z  y! H) z& R/ K! ?       end% M+ w5 [: \0 _' ?; B" [5 ~8 ?: _
end
' K8 l3 ~% i; n1 H7 k$ h  | %不存在vs到vt的最短路，算法终止，注意在求最小费用最大流时构造有向赋权图中不会含有负权回路，不会出现k = n2 Y3 G) v$ H0 o/ H  e4 A4 x
if p(n) == Inf6 x! `5 N+ ?! ^' U1 Y0 x' l  V
      break;- b1 J, a3 c7 r# F, b( l# \
end) q( a1 e/ w( X: |8 e
%进入调整过程,dvt表示调整量
, w  U7 B+ Q" q7 T; j" l+ W) } dvt = Inf;, f6 X0 i8 `7 t" Y2 b
dvtt = Inf;* h0 z0 y. _) m+ d, a. K/ L
t = n;! X  o- f7 j  f* j5 w2 S
while(1)
) Z# F. p# X0 S1 _4 h       if a(s(t),t) > 0
0 d0 V0 C) ^2 h5 }( @# v" o( Q          %前向弧调整量8 G& X1 \! U& ?- Y, m. _2 s
         dvtt = C(s(t),t)-f(s(t),t);( r* Y- F9 b+ ^+ {
         %后向弧调整量
" i9 R( r9 E3 Z: q" d6 M9 I( s- F! O0 T       elseif a(s(t),t) < 0
& \" |- \' w- b  W          dvtt = f(t,s(t));4 L- O3 `+ i: x
      end6 B% b5 v& v% y6 j  A3 S
      if dvt > dvtt
5 |) n# l" Y$ {& c          dvt = dvtt;
$ s  O* D- Z7 v5 U$ o7 ]7 g       end3 y2 L& f' s! N0 ]! v# V* O- m
      %当t的标号为vs时，终止计算调整量
! m' t, k. m* B       if s(t) == 1
1 N6 l# a2 D) i& R# Q+ p! V          break;
1 |6 b5 i& ~8 W" ?" O       end
& r" y1 X. A6 K' t       %继续调整前一段弧上的流f- G' T) R! ^( y$ C4 z  n
      t = s(t);
( u  Z1 `1 @+ j& K- n end) n2 p: v) Y- U! k
pd = 0;) j4 u- F. R0 ~5 h5 w9 }7 Q# T( q
%如果最大流量大于或等于预定的流量值5 [6 R& H5 \. }0 z9 K: R) z$ H
if wf + dvt >= wf0
) r4 ?7 F/ \% b       dvt = wf0 - wf;' c8 x% ?5 }$ f' W3 |
      pd = 1;
( a; A, i3 w1 ~ end1 j5 h6 K" I0 T, O
t = n;$ Y$ q3 ~+ q) d, ]) L0 j
%调整过程
( X  O6 D2 M2 I& w  i; J2 P while(1)2 n1 X$ L0 H7 Z' v
if a(s(t),t) > 0: w% a+ B! P8 x6 f* m6 @
      %前向弧调整8 C9 d! y) e/ F3 c3 Y7 }0 e+ a
      f(s(t),t) = f(s(t),t) + dvt;
9 ?3 m) M0 y/ U elseif a(s(t),t) < 0
: c1 l* g+ A( `# Y/ Q       %后向弧调整
, }0 u# R8 ]7 w8 p       f(t,s(t)) = f(t,s(t)) - dvt;
4 C+ d$ c1 ]' s' u1 _8 O! E end / Z3 d$ |' M' Y8 `$ D8 C
%当t的标号为vs时，终止调整过程
$ }5 l: L; b7 D8 J9 } if s(t) == 13 T7 W; A! h3 @5 O- V
      break;* b3 F. m  L% i. v
end2 Y% _# U( K+ F# t; O& J
t = s(t);6 K0 f6 s2 b# @5 X  t4 J
end
; c! O( V! {+ Y  { %如果最大流量达到预定的流量值) E5 ~% z# x5 N6 Q  o
if pd
, _6 S# P, ^* K3 R% U9 r7 E4 U       break;; [* h; B& }: ]$ H* v* @
end# e* \4 ~4 w0 m- ?; X+ x
%计算最大流量
/ }% A9 Q2 b( _ wf = 0;
/ Y) ]- k$ Q( }# H- p8 W for j = 1:n* H+ s( b9 ]- _4 V1 o
      wf = wf + f(1,j);7 v% p) w  U9 _
end
; p$ }0 ]3 j2 B  Eend9 D! z' j4 ^4 l  [7 c* M9 ]
%计算最小费用; F1 x. j8 }6 N% V2 f/ I& h
zwf = 0;
; M3 G; p/ K2 ^9 S5 H- U, bfor i = 1:n- |; X/ q) [' q+ ^/ X* ~
for j = 1:n
  Y" m$ \, \# Y, u$ D& r2 d       zwf = zwf + b(i,j)*f(i,j);
/ r( T* g' A8 h% b end* W: H3 g% y9 M( _
end
4 i! u9 Z* _* o6 U$ V8 \9 a%最小费用最大流
  \7 ^. ]& f0 H8 xf2 r' N7 i" ]) K- Y# p3 i
%最小费用最大流量# U) o4 B) \* ?' }
wf
5 ?5 V0 A  O$ r, y: F3 Z%显示最小费用
( g+ g, c7 v) B, S" D( y" u( a' e/ hzwf
4 D7 k  D$ r' w1+ I9 w5 Y( J! T: Q4 u/ l9 C
2+ Z# J5 f. y8 Q4 Q: c, w/ a
3
% \) m1 O" x. @4
8 ~6 T+ F, S% V7 V5
( ^. x4 E& F) a% G+ n' }$ a6
$ I& i' n2 n, u5 S7% y9 I, e3 N7 z. K1 I7 J% ~& j) F
85 N8 h! G6 |7 f! a8 L: [. K; a
9  W# ?1 [* J4 }: S
10- s* D8 ~) m& B0 n6 t
11
+ i! Q, k% q: s5 F2 y. ]/ U2 W12
) y! W* e2 `/ ^  N! s4 n3 o135 c4 l9 p0 m# R! L5 \3 t' U
14% r% R- j+ A! P8 F
15  S6 t% k1 }: \1 {8 h5 N3 }
161 w1 z) o9 E4 f2 y/ a
17
/ b0 u) S* Z  A9 S18; H' ?7 H' y3 v
19
& h* T- z6 O3 ]4 o9 J: d1 E20
; q+ o* v- [$ @9 b; M7 e) E21
. f* q. p- c) n- D22) s" `# d0 O) }) v- x/ z
23+ [9 a8 u1 Q$ \9 p3 ~
24. o' [9 M3 T& B& w- [
25
0 t# p" T* a6 }% [+ }26
) \* t; ]. |% c6 b27
- _! r4 Z" h3 g, t, F28
' L* J5 A" W9 o! q/ [3 q( r. Y29
! h# \7 B5 Z, c, _30" R3 n/ M" Z. K8 f1 t" j, e
31- q* `6 V7 e4 m6 _
327 S, I! W, g: Z3 i0 C
330 }8 c( _# G) K# x' q' ]% r- W
34$ Y$ f) k3 Q% i" S3 x% X& j
351 a( d, g% v; z* }- Y
36
7 G# ~  S( K. w3 K" f" d37
6 }; S+ Q6 g; s# b1 J$ G38+ `& R" @' ?* O5 P1 p7 e, U  P+ h
39$ T2 }  e7 ]7 D% Y6 G
40
- Z. K) v& f2 o0 }41
3 O7 O. r* T) m) c. y+ t& x42
. @2 j2 F# w- d" S$ K* `43
+ f& w$ t9 B! B# ^2 _44  [8 M5 D/ }; _' s4 ^0 i$ v" }
45( y& p, t  Y  O9 y
46
2 G) b' P. o- `$ j4 u47
- A1 i8 E4 t: U; T& m% r' ~48
' u( x% a, l/ }) u% c" [5 i49
/ L( w% Y: U# Z2 j503 Y7 B9 c; M, V( W8 h( T
51
& f* H9 c& J2 k! W( ~# I52* r# E( ~; c2 h$ F/ }2 `* d
53
4 b) a  Z; g" P' L$ U" a54  ~! \7 ^& {  m) _/ X" w
55
) k$ s2 F2 {+ k3 X) x569 p' A- ?; Y4 B) [* f! S
57
! h: z9 k* d, T0 d' B' C* b58
% v( r3 d! B/ t  I5 W* U# G0 C' U59
. H5 }+ o2 Y8 ^# @# m60
: n% T! k. ^2 L, |8 p' a: D5 r9 b61
% @/ ^5 R3 [/ d0 z62
4 ~- U1 o& s  ]634 M- ~; m6 M  p; U; P3 U# I
64
# S# w5 Y+ ~4 O; v  w2 ~/ r& L6 J65
/ f3 s6 \5 T" i; a6 q5 F66" X$ R4 R4 W  E. O+ {0 a
67
- I; Q/ d7 I# D* d: h+ G7 F" _68  G- I9 M! P) x) i( o5 A, \
69& n2 m' t/ [8 _" Q( t
70. D3 y' j4 p0 |/ R! M5 w4 u
71
* e) ]  N) r' ]0 `72, ]7 K$ |; a& r% A& N
73
8 ^+ J! W9 j, [2 J* W74
0 P4 @0 m* F4 @75
- K, l$ y3 U$ R  E4 v0 G76
, W2 A7 E4 W3 |( L7 C0 g/ P2 V77* z7 {, e9 d1 ?6 q  l( N. z
78
791 ]- l9 C) @+ e, P/ ?
80
+ `8 J+ W9 Q2 y& f! Z81; e5 p8 |2 Z" b- p' h8 J" U/ k5 x
82! @" D. f+ P* s7 `: U5 n7 E
83
0 Z  _- G5 Y/ i( G+ n: N84
0 d2 b* s$ z8 n+ F85
% [9 A( f. W9 E% Z86" }: z. e: h" B: t  I
87# `( A/ Z- t( j& k% G, q: z/ k
88
2 C0 ^8 P+ \0 Q9 [( B% c  g9 Q89
, Z' z7 V* T/ H2 k90
# Z7 K6 ?# n! o: L8 b91! q0 Q/ p$ i# q
92
3 A  C& b# x$ i% e93" `5 X& ?* f: y. r- x. T
94
9 [4 k0 n$ m" o1 u* z, b) M95
4 f/ M) w' x8 K; {96
' a! c2 a" [0 x4 o" m' ~! x% q8 b97( [* x, Q: l* N
981 g# O7 X8 U' b* `1 ~7 Y& D
99
# L/ t9 h/ |% P5 O5 @* k1000 m, I8 @. e1 l5 [
101" }% D) N% `! L9 p/ B( T( t! `. i( q
102) y) u; R! X/ m9 D4 t
103
* P7 F4 j6 {6 A# h( ]104' u% }$ V8 g3 q. Q  h  V
105& q/ j3 Z% r+ N- Y' V
1062 S% g( e1 X5 `! j0 a/ w  i9 [: r
107
+ z2 \2 I% B; D2 d1 ]* Z1089 o3 x0 C+ C) Q9 o# @. }; z% R) T
109/ g3 N. A5 u9 @- Z  i
110
7 A6 k7 Q& \# V, L) c/ j0 {111
; t+ v9 M! w: i# V1121 g3 \& W7 `8 [2 V0 p4 X- t
113- z+ \9 C  }( L0 y, }0 Q
1144 U* I$ `$ v0 s# f' ?- O; ^
115; ~* f% o& c" ?0 T0 q* o
116& y- s3 P: Z; P
117( J+ g! U. [+ V+ f: S+ b$ h) `
118" I) {% ]# f8 K6 @: l
119
( y- S0 ^& n3 Z+ f: _6 |1205 j5 G$ y4 v- P2 z2 Q
121
3 l+ G( o- d  f+ Y, z122
/ D  U: i- R$ D  c, A123
7 R( X2 m; L% ]124* X- {, N  O0 X' H" |) E1 y
125/ T! e8 O+ j$ W( c- v
126
' q6 d1 g9 m) j* }7 F; T127
/ i. i6 e0 F6 J( v7 d/ L& a( c1288 j$ |0 |8 X" D6 l+ i! X( A
129
" J5 x1 o8 ]5 p* n0 w1 x130, p6 ~, E8 }8 ]
131+ C2 I+ I2 |! H; g# S  H
132
, {8 g, ]% G9 S1 M2 d1337 Y7 [$ s1 r4 l
134' `* r2 {8 G: I/ C. _1 Q4 Q. ^: `
135
8 n3 j. t& e6 |4 t% B2 B136
3 Z3 w( n# x6 s, `8 k+ i137- e# X7 K' d: \) z% X
138, _0 u' e- h4 M  }% U
139. K: j( J, F3 M/ E5 w' `! B. K
140: x7 S; Y7 J, Y/ Z) [+ @
匹配问题：, _8 g8 Y3 @+ S

+ B# O1 R/ o7 W3 t6 i5 ~) z最大匹配：
5 ]4 U2 J2 T, _  }
9 E, r0 f2 A1 X' U5 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