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

最小费用求解! o* A6 ^% c2 j

' U8 u. q$ U, j- Y) ~Lingo:
. J5 z! r% I+ {3 c  u* Y
! L( u& w. R8 f5 n! W9 M& Nmodel:
# v+ \. H7 _. A, `: b( Csets:; O+ S4 u5 r" C% _6 q
nodes/s,1,2,3,t/:d;
) V0 T$ _* e) l% D! zarcs(nodes,nodes)/s 1,s 2,1 3,1 t,2 1,2 3,3 t/:b,c,f;
" J' V& f. x  \0 [7 o/ r& Kendsets# ^. f2 [# C" b6 o( l/ A& n
data:5 c9 X  O$ X4 S" W" @3 P6 B
b = 4 1 6 1 2 3 2;
9 o, i% s; F  ]7 C# W; x: yc = 10 8 2 7 5 10 4;
! R+ l; _' x0 z0 y8 t0 _! Ed = 11 0 0 0 -11;
2 W+ E* L7 V& Genddata0 ^, t2 Y" u: L9 \; g6 T3 \
n = @size(nodes);
. V5 e& W( O! d- u3 l8 E6 wmin = @sum(arcs:b*f);9 l9 R: X8 n) l+ E! k5 r1 q9 C( e
@for(nodes(i)sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i)) = d(i));% S3 [1 F" w3 e2 A# \: Q& K4 e5 |
@for(arcsbnd(0,f,c));
6 {2 }- t; q+ s# pend% f& V" [& d- Q5 Y/ W/ q6 K" R
1- ^2 L1 Q. Y0 C  a2 `- M/ l
2) \$ z" d- c$ d8 l
3. F! D( J" h( y7 [7 q, Q' }
4, D3 u' ~, I& T
5
% q# R: r" u' k' D# o' l- ~9 T6
5 G* N% ?# P5 ^7/ s+ O# y9 H0 R% D: R
88 }7 w3 ^; K/ g% a
9
2 G7 K- z( U9 ?; Z10- J% v: U+ R  A' \
11
+ F1 r0 R. X+ n$ `12/ C/ ^8 I* G) z9 T# Y
13
- m( }5 E$ J: b/ F& F8 E9 T; U14: l5 N- k* g9 O, C. D8 B9 R
15
- ~6 q) n/ M  ?* v6 B8 JMatlab实现：) N* Y% p& i! L; w2 ~3 h
: c2 `# y) X: s" E3 {& P' C$ K
: Z7 Z3 v# |9 G( Z. W: @2 R6 N
n = 5;
* m/ F. C" t" m3 ?* `$ n+ g( `%弧容量
  Y8 A% k4 O  {, n. F& ^/ O- Ca = zeros(5);
0 L! [6 N- j+ c- Pa(1,[2 3]) = [10 8];5 F7 Q# S8 q4 n2 k
a(2,[4,5]) = [2 7];5 [$ }! e# Q' ^. k
a(3,[2 4]) = [5 10];
& i5 Y4 K4 v0 Q! ~* g3 R2 Ba(4,5) = 4;
* W7 K, u/ G4 ^4 D) KC = a;
' @( @9 O7 h% 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];
, j3 U: C& Q, @* _9 Z) @: V%弧上单元的费用6 B" v5 `% }, g5 E- s' N* ]$ a! O
a(1,[2 3]) = [4 1];/ M7 r& K. \' G9 H
a(2,[4,5]) = [6 1];
. X1 A/ z% Q/ d; l; Z& a- ]# Ea(3,[2 4]) = [2 3];# S2 G: e$ W9 _) \. w; ]' m9 o
a(4,5) = 2;
+ Q  k& }' V( ^b = a;5 ]: [& [9 A7 G$ n0 }9 {
%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];" f. e" O# |0 w; w) r
%wf表示最大流量。wf0表示预定的流量值2 r0 Y' n, X) Q: Z+ C8 a, ^; n" D
wf = 0;
9 G5 P- u, w  D9 w* d# dwf0 = Inf;
0 M9 V6 B9 K) Z! Z- ?4 Q4 D%取初始可行流f为零流) W/ v# x3 E5 O, M: ^
for i = 1:n/ x: v) |; ?, p
for j = 1:n
8 i: k1 K' m, {" m% [: \5 D       f(i,j) = 0;
6 j9 r- r4 b1 m' r9 n end
7 a, c5 ~9 s1 z# Aend9 v6 X% }. M/ E# b. H
while (1)! @+ p6 p+ ^  U/ d4 ]
%构造有向赋权图
% D. |+ I- E! p- M- }1 Z for i = 1:n: [- }' Z% B5 R, U4 A+ ]
      for j = 1:n
! i% L& S( @, p2 m" J" a! G          if j~=i
, V: \/ O- J% _+ H6 c( \/ M             a(i,j) = Inf;
5 {9 K* o$ {+ D4 W9 |6 d5 ^          end
5 {3 Z, V8 S% i+ h       end' i( H" x! f* `- f* w; G# V
end6 t: f8 s' G% b: D& q5 T; d/ |6 b
for i = 1:n3 l0 c: I) q- O. y1 y; ?3 Y4 g) p+ R
      for j = 1:n/ c" u; a! f6 g2 M0 h9 L
         if C(i,j) > 0 && f(i,j) == 0( C& y6 a- l" g1 u; `
            a(i,j) = b(i,j);$ J. Z" N3 V; |, a% o
         elseif C(i,j) > 0 && f(i,j) == C(i,j)& J; \, `" O2 g- }, c
            a(j,i) = -b(i,j);
% y  b) b8 \9 m          elseif C(i,j) > 0
4 @; Q6 {& ~* v) v; g% f             a(i,j) = b(i,j);' W/ `4 Y/ N2 L/ l% \) _
            a(j,i) = -b(i,j);8 `$ v4 z7 L+ Q) W, h3 n* p
         end
, N7 o$ B+ D8 [& X# e) u       end
8 n& C1 {, m" z5 u4 J end" ?. M5 ?. n8 l& W3 e/ k
%使用ford算法求最短路，赋初值4 g; `* j7 b7 }2 z
for i = 2:n' P8 n! C, ~: q* z& T; N
      p(i) = Inf;
9 A/ D; h0 t+ @' V1 P( G       s(i) = i;
) R# ?% C" y, X9 K* y end3 N- v' A. u" X: O9 A: Y! D
%求有向赋权图vs到vt的最短路,赋初值5 ?+ w$ i  x+ |
for k = 1:n
0 L& q9 n! O* f, ~: Q3 D4 B       pd = 1;( \3 g( e& w$ E2 Q9 q4 K# j
      for i = 2:n
& x0 f! Y: Y( u, u* ^          for j = 1:n
/ l; B( r5 T2 T: F. P             if p(i) > p(j) + a(j,i)
6 t  S- |! G9 B  ]4 J, x3 X) c2 t                   p(i) = p(j) + a(j,i);" U: }/ [+ X6 e4 o% G0 X* a% J2 W
                  s(i) = j;
; b# E. T) x7 Q0 c" m+ O; i$ ?4 [                   pd = 0;1 _4 U/ e  l: [3 Z# L  d8 O7 Q
            end
9 n4 Y& k1 d- C3 \          end
! h) b6 M$ @+ w       end! n5 g* n4 k: u2 _9 N, _# [% z
      %求最短路的Ford算法结束; k1 F% D! ?$ k& D' L) \
      if pd+ N4 @. `' l- u( b0 Y  n+ {
         break;: y# f! b5 l$ t4 j
      end
* d0 x# Q9 u) U/ a# Q end3 S, m$ H$ E, F* @5 U! j
%不存在vs到vt的最短路，算法终止，注意在求最小费用最大流时构造有向赋权图中不会含有负权回路，不会出现k = n
3 K/ @2 p( `- e( B if p(n) == Inf# K# Z7 A+ c* U8 c! a
      break;
# J6 X$ x% @% C. P3 U! X3 C end5 G* |" d, M: G3 n
%进入调整过程,dvt表示调整量- h. g5 m$ y+ w5 g# w/ ~
dvt = Inf;
3 g. w2 L7 Z1 Q, w dvtt = Inf;( f% E- Y. h3 T4 ^( |8 s
t = n;
" \8 a1 {" @  F while(1)- S- g( i: d1 s. E! @; H' v
      if a(s(t),t) > 05 H* G* H- c9 {( q( q/ T- M
         %前向弧调整量& h* e$ O9 `! c- K4 H
         dvtt = C(s(t),t)-f(s(t),t);+ y0 n1 W3 p1 `6 c# i" J
         %后向弧调整量
5 I; r' J4 w4 q* U# v, A/ L       elseif a(s(t),t) < 0
2 t& z: y2 l/ K! X# J) n          dvtt = f(t,s(t));6 }' r1 j& l1 a+ b
      end
: y& p$ W6 z) M* O. m       if dvt > dvtt" B. i: y* O5 P( j$ `
         dvt = dvtt;
, G0 M6 x1 c7 L       end
& {0 h" O. x/ O       %当t的标号为vs时，终止计算调整量. J5 J9 Z* [0 O! {  n) e7 D1 U: M- G
      if s(t) == 1
+ [+ F3 E- O+ w7 O! \          break;
$ c' A8 v' C: ^$ r3 O       end$ A7 T3 \! F& B( @# l! U6 A
      %继续调整前一段弧上的流f
3 s: Z; V( [9 r% D* t' h       t = s(t);# C! D2 X* [( k
end
' L5 e3 G9 @0 ]1 G$ T. ^( ] pd = 0;
9 o1 S1 A$ K7 o4 X6 s' O7 c$ d %如果最大流量大于或等于预定的流量值
4 a, h0 O8 N6 D5 }3 N if wf + dvt >= wf0
/ j5 p/ K- o- X9 d+ t/ f       dvt = wf0 - wf;
; a7 z, p! q( W5 q6 R2 N( u       pd = 1;9 n; d8 D8 g# u: I( l/ x5 B8 z: y, ]
end
7 k2 s8 e% F& i; t8 H t = n;
. R, P0 H( Q) s( @0 S %调整过程
/ g8 X% H, {8 X) `2 N, u' S/ y while(1)
0 c; M- S4 T# U3 ^ if a(s(t),t) > 0
  Z/ u. J1 C7 A; v! S+ J       %前向弧调整
/ H# H3 G; P8 f) N4 O5 a+ r5 ~; I       f(s(t),t) = f(s(t),t) + dvt; . M% a, N  O8 K3 y
elseif a(s(t),t) < 0: e8 p1 f5 L  {$ |/ q
      %后向弧调整: {/ ^' i2 ~% v6 _2 T* y
      f(t,s(t)) = f(t,s(t)) - dvt;
+ v/ o5 ]9 w) q end
2 ^3 [0 J1 q3 s7 o %当t的标号为vs时，终止调整过程. p' [8 E( l" |+ H1 _
if s(t) == 16 S1 t9 _' H" N) X8 ^
      break;! Z: w# Z& E8 t( S& s! y
end0 ]8 z2 g( ^) i: v5 O
t = s(t);$ h! q' i. t* ?9 I
end; F; J) u7 W) R3 m8 |
%如果最大流量达到预定的流量值1 ]2 @; f7 j7 p" D' \' v
if pd
. T# l9 Z, f, J( g3 a8 o! I+ ~) G       break;% I$ d3 E9 [0 N- f5 Y0 ^. |% h
end
3 i0 e+ V# a1 d6 p %计算最大流量% e/ F5 K% j, p6 N6 d
wf = 0;2 G+ N; U- q( t. ^0 i
for j = 1:n0 V3 ]7 N& T9 ~% f
      wf = wf + f(1,j);* n8 J+ @( G+ T' q: @" |7 ~- ?, F
end
0 [3 j: z4 W4 Y" eend
3 {8 p6 h" {3 z$ Y: _. K%计算最小费用8 K2 X9 [7 I/ N( n
zwf = 0;1 X0 J8 @8 T) m
for i = 1:n( C0 M+ q5 B5 {6 L' O  i0 R' x
for j = 1:n$ S- q: ?6 G: \2 c$ Z% D, ]$ a
      zwf = zwf + b(i,j)*f(i,j);
0 _4 w  ]3 f2 F1 y' ^: E end( D4 i' |, M% Q% W& ?2 h4 q5 t8 Q
end: o1 u' h. @' q  @5 q& n
%最小费用最大流
# m* b- z1 d7 Ff
0 G0 v! o8 }" g& ?7 c* F- l2 W$ N$ Q% q%最小费用最大流量+ t+ [+ f" k1 K1 p1 H- ~% ^0 s
wf
. O3 J1 ], o9 Y; u%显示最小费用
" Z# W1 F) c# g9 H4 U4 kzwf( t1 [$ n3 z" d4 U& m
1
  q6 L4 g( q- @2
; D$ e5 s0 s( W/ N) X' M1 E, m32 u0 s' D! u/ l4 L
4+ t% _  a! f1 s0 l, u
5" z. d5 l3 V  ]* r. m
6
1 I' e4 p' ?0 }( T. c5 o72 d0 m. g: s  b3 v
8
* M7 J; A/ ]9 H6 q9
9 K* q! R, h5 A' Z5 M9 f10
$ h: |0 Q' r/ f; @6 Z9 T111 G+ j5 ^! m/ D
12
8 f+ I7 ]* q1 h4 A; G132 i+ _- p7 l% a. u  ~
14
7 @5 t$ h: Z8 U. U3 w15
* f+ N1 w# w, u: B& s16
" g: l0 v" P' u: @+ v& f& a& [3 `17
9 Q, T: b0 L- @, q% }* Q% l' d187 H, P5 o: S' b( d' L$ I
19$ r. x# U: J. D: |
20, @8 T/ L" f; L6 q8 G" g/ G4 r
21
2 v7 B' C$ h2 U+ X& f3 C8 n224 {  y4 k* |  r. b: Z6 d# i
235 P8 N+ o4 K6 U% q$ Q' u0 x# O
24
: n1 F6 f7 y; U2 _4 e; d5 g% K25# g2 K) q7 E% \  F4 Y3 a, K9 @
26# L" @* B- D; V  @7 [4 S
27
. ~/ o/ q; {* T287 Y7 G4 x% }, W$ b6 i3 q/ R3 R
29: d, x1 u  q0 o  A" |' ~* U) `# U
304 W1 [% R' V$ Y. i, p
31
  M. N2 v8 t  ]4 o7 J0 ~( P$ Q3 J32
9 c# w* D$ `* m336 z2 `% ]/ ^% T( T) S# U3 u
342 w4 h3 A: j. _7 D. Z
35. k: ?$ b  r& s7 U' E! g6 W
36( L6 g5 y6 B2 G9 U
37; a( k2 C/ D$ {  \3 J
38: d! w+ @# K8 S! m
39% z4 {) }8 A6 Z0 v& S( y; A, K
40
' n( l# q9 F$ b4 f- {+ w41
. P7 y/ Q; H0 w42. s1 m9 f$ m7 e# x
43% L. p7 ?6 A- K! x& u
44
7 I" o* F2 f/ c  M/ M" `45
* O1 C# B( ~! H7 x46. q9 T0 R! u. J3 j+ u
47
* i" n# w! h! |2 @7 a! m. Y/ e% D$ ]6 e48
8 W# Z4 \! {- X4 q. A49
! b+ T4 ]* b, u! Q, R50$ X% a2 k* O& t3 M6 l4 \
51
$ `' r4 V1 R1 x6 P/ R; e52
2 D  K2 D3 K5 b9 h' \7 n5 u1 H! f* Z53; Q# a5 O  N! H: D5 S# Q9 x6 Y2 @
54% @& \4 n( c& h# Z) b- f
55
9 v4 w1 b% R0 K# D2 m56
1 \) `% C. G6 r$ E" r# |! U: u( u57
! X0 K. h2 D( h& |, q58! h* w0 G* W, u6 G* g( o# t
596 _7 a5 D/ ~' l) @
60, S4 q  m; E# {, I
61
5 F# y: i- `2 `# b, z62; |& Z4 g8 ^) w
63
3 U0 R( ]' J0 @) D9 C/ ?8 s, U% l64
; C1 H: h+ ^5 i& C6 X65. P2 L' \8 H" `* a
66
1 b) p7 F6 {5 x: I( y$ |, D67' E- w  @* L, h- S" w% g
68. l+ Z9 Q, x! K8 {/ v" ]/ h$ x
695 D+ C( K. s4 ?- i: |+ i: j
70( G- {# O) T& R3 D* Q1 X
719 X' g: y! R+ \8 H  m7 C
72
4 c6 O5 x) ^8 s! h* X73
0 ?- S: R: K* J) ]. s74
# ~& u( ^- a: d) Q9 }/ {75( @# B5 w2 [" k1 Z
76; V& E# N4 L  U% X" Q/ O
77
6 K& e0 \/ w) V* v# s1 a) p* y5 E78
3 y+ G* Z2 n: d4 [5 Q79
# q$ V  t* N  m0 Y% _7 c80/ m- I9 e: y/ w2 @! U
81
+ X5 H. \9 Y. O82
% J# V# m3 Q3 Q& n83- u6 U& b8 u: ?; J: \
84! R3 f# M9 U/ U
85
! }2 Q: P3 y$ R! X7 x86
, o% u$ e. |4 V4 z. R- V879 r# n' o! W% p. G9 ^& c6 H
88
2 E0 _" r- C. a) r- Z$ X5 r89# V( E2 G9 B. D3 _
90( X( y2 ~8 N' m3 f# ^% i3 z8 \  v
91
9 g* T0 h+ W, l. s" h/ e$ S; J" h92
& n9 X  t0 t3 Y9 m! r1 m: U8 d) [93
3 t+ s; i7 y2 j  N4 g) ?94
' L  O& {4 D: o2 e7 O" @95
4 X  X, w% {; S( b' p96: [* r4 a4 S( Q( t5 i
97
% }' E8 a4 b7 I) t98
' M$ e1 }! \; P99
2 a0 Z; i3 ^" F3 S100
3 s' m) _5 h5 a+ t0 a) b1 [5 }101
8 E* X6 q7 [; v- a; u; M102
( U& ]. x% L. q$ q2 B( p1035 R6 d9 a% c9 U
104* i0 r% \  A! G5 w0 Y
105$ W4 G' I" p" B0 B+ R1 ]
106' z7 z: h3 a: o* e- f+ `
107
9 {3 _. e7 g% T$ A; i2 I. G108
* {. M2 h+ n( e# X- v7 r# o1 r& U/ c( T109
5 Q; M* s) F) G/ }4 q( X9 _110
1 k- T, k6 s# k1 X+ ~  h3 U111! ~% z: H% S" k* P' K2 P' D
112. D- W% C( k3 S
113/ R- C( w+ e5 C1 J1 v. t4 U
1149 g. v2 G7 z# r! J6 t8 T9 u
115. a8 A3 L( ]) B1 G3 S, k, k
1160 [) P2 ^3 Q; S& l% G$ U
117
) Q2 A- @( i: M118. m& ^, G- I% ]; i) ^, P) H& V- A. j
119
+ _5 s! o+ w6 t6 U# X  L% ~+ _120
8 P) O* T3 K5 N0 [8 U" r. M1211 C) S0 ]3 @' D8 _, _; {( ~* }
122
: T/ K; `8 h0 E: F( I123' Q4 w! a; M% x2 V: z8 s* Y
124
+ T% a8 W: x" H  f) X; x8 K$ |" j125% W0 R: n# x9 h; P; R  D- d, B
126
1 {6 [! l7 o# T127
9 x2 p8 l, a1 C128
+ t5 w! ~9 O) ^4 q( B  s129; R$ T9 E! P$ l6 U2 B  D- u" X  W
130% [; ^( D" o: c
131
% |& ~( W, {* q, H) M& B1322 |9 |/ a& i6 J) u: o
133* R4 d9 U! {# E) K, E
134
# f# b7 S! w0 s: @. e$ S3 ~6 S135  M1 w* B2 V5 c1 i# e
1363 N& \8 h, r( E# k; g
137
6 Q# p6 r! @" p, v138, N" t- Q, Q: ~3 c9 e
139
: ?& |, P  i, g9 g. O8 x. K140
2 d. X) O# c8 m1 A! F1 U匹配问题：1 h! n7 s, ?1 M- ^$ T  a; \
" k% A4 c* f( C& E, h+ F
最大匹配：* G- T) n% u& \5 g5 Y# g. N, n. Y6 ~" F

2 A- ~) _2 {5 Y! b& o" y
代码实现：

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