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

最小费用求解- G; ?1 [0 |' C5 n6 V$ y

  g, }- x$ G% Q, C2 F2 H7 uLingo:
! }6 B$ C/ n2 A1 y/ q$ Y" n3 R0 `6 f- P6 n
model:
% v$ q* v4 X% m: h% ^+ csets:% \# |. h8 |5 L6 M6 u0 l9 `, Y
nodes/s,1,2,3,t/:d;: L  t( L* g2 [, u. Q. f
arcs(nodes,nodes)/s 1,s 2,1 3,1 t,2 1,2 3,3 t/:b,c,f;- ^" t3 |5 r! t# `0 I2 s6 Q% z
endsets
+ ~8 J! H- U$ tdata:
7 o3 {  p& h0 V: S; Kb = 4 1 6 1 2 3 2;0 _# a0 c, f: m/ g
c = 10 8 2 7 5 10 4;. L$ }! j( v. Z1 ^' q$ Q  M) c4 \
d = 11 0 0 0 -11;
6 u8 Y( M+ K+ J: f6 _8 n7 t* }enddata
' A1 G# n) O) ]# B  N, q8 `n = @size(nodes);8 l: d9 X5 @( D/ [) `& e2 s
min = @sum(arcs:b*f);! R9 {4 ^0 r- F" J  D7 Q( q" V
@for(nodes(i)sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i)) = d(i));7 Q  d  N& P; q: a( m
@for(arcsbnd(0,f,c));  ?, b! Q/ K/ A4 s
end' w! V3 }; z: `$ B5 [
1, Q* \' C  f# o7 |
2$ s  b+ Q! r" f" w% x. m& M( B+ h
3
, ^. I" Y% z$ ~/ E! K& y" X4) W; C. a7 P# U1 n; g- N, C7 Z6 \
5  ~8 `/ X. c5 h3 ^; d' U
6* E: j( c. ]5 k! g; b3 Z# _
7/ c4 O' k! W0 j# e4 Q  c* |
8- I+ H% E4 c# }/ r' }
9
# t# s% \) J; x" s' I6 M10& {5 G; s1 u2 X7 M, I3 ?2 |
118 i* n; i- F* r- [. M
12
$ n6 |0 f  B  {3 e1 F13. X4 R7 ?- J: q' T0 x
14: {8 x; h3 s, \* w1 b
15" o+ b6 n% K  S  I3 M
Matlab实现：  x; [. C' C5 w
0 v$ A! v. P3 E
: k& _! {! q4 W5 [
n = 5;* v  }& I; {1 h7 C$ Z$ D) a0 v
%弧容量
7 Z- ]  g* L# J/ F" ha = zeros(5);
7 _6 b$ s: c" Ua(1,[2 3]) = [10 8];
) i' B' I7 C/ Z" la(2,[4,5]) = [2 7];1 U9 E" a7 \" U! Z6 b+ [
a(3,[2 4]) = [5 10];8 n+ N8 o1 r9 {4 N6 K- _
a(4,5) = 4;
5 ?9 a& r2 R! T7 G. KC = a;+ [5 a  S# z6 [; b% n' f( w! V
%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];' j  l* l+ f. q+ I
%弧上单元的费用7 k( {( f0 [( r+ B  p7 u- H/ z
a(1,[2 3]) = [4 1];& |0 }: R4 X- j. f
a(2,[4,5]) = [6 1];( }5 L+ d; s& i& r1 Y# h9 U9 L
a(3,[2 4]) = [2 3];
6 U  r/ N2 V. u+ v6 na(4,5) = 2;
9 L0 \% Q6 l- Tb = a;6 m( a8 R: K" S; j4 }" \
%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];
+ m; G" Y% G! ?) `/ |%wf表示最大流量。wf0表示预定的流量值
wf = 0;- P- {7 }5 Z& u6 x0 S- |
wf0 = Inf;
5 E: C7 ]# O6 p4 b1 c3 }& I%取初始可行流f为零流
# ]/ Z2 f6 e' u6 ]+ |for i = 1:n
; O$ M7 X# H+ g& K! G- \+ ?3 j for j = 1:n
/ v0 z9 H2 _+ n, T       f(i,j) = 0;
. {0 y% ^, `# |9 I end
/ Z+ S) t  [: |9 xend# a9 ^2 A! K. ^& {
while (1)$ T) I* Y# E  W$ M& m+ z
%构造有向赋权图3 d* W' w  h8 X8 ]2 R/ M
for i = 1:n# k# X  S$ L0 ^. A( o4 L
      for j = 1:n* C  R" c  R5 j0 m  u. O
         if j~=i
& i5 s1 e- m7 D: j  Y             a(i,j) = Inf;# ^2 R% x0 t# G' \; g$ A3 o  l
         end
! [9 S% u. w% L2 x8 i: {       end
6 L0 e% ]# x. H5 [, S0 d end
. ~7 j; R9 Y6 D for i = 1:n- z; C- D1 z0 U' }2 F
      for j = 1:n" }' [6 U& N6 F) U9 J2 [5 ?
         if C(i,j) > 0 && f(i,j) == 0* \: d* M# B, y* p
            a(i,j) = b(i,j);
6 ?# u8 O0 Y% \& s* x' e          elseif C(i,j) > 0 && f(i,j) == C(i,j)  ^& {1 r( M' _5 J
            a(j,i) = -b(i,j);
- c: R7 z/ b6 Z# w/ S" y, {          elseif C(i,j) > 0# X, d) P3 j: m  j, g# Q
            a(i,j) = b(i,j);
. D1 [0 [# I* r0 ?& f$ M$ S9 R             a(j,i) = -b(i,j);
# ]- S7 z7 D- x          end6 }% N8 y6 G8 a0 n* v3 b% I0 E$ Q
      end5 X  Z  j7 I: S; g2 I
end
. ~6 l, C6 ]. f- k. X, W' Q- P. z %使用ford算法求最短路，赋初值' K6 Q* g- p8 q! j
for i = 2:n' x, U6 ~( P& i% I- U2 L$ f
      p(i) = Inf;1 d. V7 I& |$ M2 h. m& T. U
      s(i) = i;
. W5 k  f" y  a# Y( F6 E end
, W" W6 \4 ^7 f  E4 A: e %求有向赋权图vs到vt的最短路,赋初值" u0 L$ _- E  n" a
for k = 1:n$ M8 ?9 S' U# d5 ^0 t0 H; X" k7 j
      pd = 1;
6 N; L* H0 p8 B( f, n       for i = 2:n( O$ }& \( m5 a" c. v
         for j = 1:n
) S' `  i6 R+ ^6 ?' \             if p(i) > p(j) + a(j,i)% {1 I7 z' {8 H8 b; j
                  p(i) = p(j) + a(j,i);
" T1 a, W: R' j7 j" u" Y                   s(i) = j;
! A0 R* F$ v, p' \& G2 j# d! }6 Q                   pd = 0;
7 A$ r( z1 o: T2 k. r             end
0 f6 j! P0 S6 W2 T" I          end0 F6 [  Y5 F* {/ c% Q
      end
- B" Q3 U' O* m2 s% B- J       %求最短路的Ford算法结束
& x) R$ e2 G- c  o; j' ]+ t       if pd
( R  R3 G7 b" R  E4 s4 Z          break;
: V1 I' ]  n4 d2 ]: j       end
( `/ u* d+ U! [- K) `6 o5 |- u end
3 w; c) o+ q, y( F %不存在vs到vt的最短路，算法终止，注意在求最小费用最大流时构造有向赋权图中不会含有负权回路，不会出现k = n
+ D: s  z2 x( f- T4 h" f if p(n) == Inf
, G. |3 O  Q$ X: I' B       break;
- T7 Z; H1 E2 J! W& c4 J3 j end4 b6 G; P$ A, _0 t  }: ?1 D
%进入调整过程,dvt表示调整量
5 q( E) b7 }: m% r0 l dvt = Inf;
  e- P# ~4 j" L dvtt = Inf;6 [3 X- m, }5 g' G
t = n;6 ]* n: i/ P) V6 u+ O! F; q( o5 D
while(1)
, y) t; ~0 f3 h9 k( Q7 g5 b       if a(s(t),t) > 0% G* _4 e+ M& ^# p' x4 z
         %前向弧调整量
3 h1 \3 K5 U8 s6 P/ \          dvtt = C(s(t),t)-f(s(t),t);
% W5 W" p3 @2 Y1 Q# ^4 j3 d- W          %后向弧调整量5 s. o& `# W' v+ }& r  d9 V
      elseif a(s(t),t) < 0; i9 S3 A. g) ~, V, }
         dvtt = f(t,s(t));
/ t. L7 q( d0 p" ^0 O) m       end0 q1 O% m- \5 w% z( ^
      if dvt > dvtt1 T6 v: g9 o. q6 e; P& \. [
         dvt = dvtt;
( l6 J9 ~0 h1 F* I4 t2 T1 Q       end3 o+ {) N  G* h. q% g; F: R) k! O
      %当t的标号为vs时，终止计算调整量
+ {! f8 n6 d9 c0 d' z" w* ?       if s(t) == 1
6 E+ N, u3 O) ?7 h: L3 N2 v: A          break;
0 X* D) @. j! H3 K3 N$ e       end
6 Q% ]" V# W! E% w       %继续调整前一段弧上的流f8 e& M; o* @( Z2 A& x0 d2 t
      t = s(t);
9 q. w! C2 g! p7 F end% d+ k' l! l3 b+ I0 w% _$ t
pd = 0;
  ~. o' H* j5 L %如果最大流量大于或等于预定的流量值2 ~; t' t, s) U
if wf + dvt >= wf03 G$ T4 l, ~* \& j
      dvt = wf0 - wf;
. o) f/ X- k" e& i) x       pd = 1;$ g% ]7 n) W+ Y  D
end* f6 S$ e# S' e5 o! i% P1 H* f: e
t = n;7 }% }3 m/ G, w7 G# Q5 P+ Z
%调整过程5 Q! R. ]7 d# A" W5 k2 y4 r
while(1)' ]- Y' ~8 b5 p" X& l% K
if a(s(t),t) > 06 k+ k: g3 M7 r! ^0 \* Y
      %前向弧调整
1 ?8 P# r8 y+ l) p3 H       f(s(t),t) = f(s(t),t) + dvt; % \8 @8 i- y- o  ]3 [
elseif a(s(t),t) < 0
( ]! s9 i0 A9 c2 S$ b       %后向弧调整
- {; R' u( k" i0 D4 ]! ^* p       f(t,s(t)) = f(t,s(t)) - dvt;
; o$ \6 @& z8 p( J: o, w end
/ ]6 ?2 k; }/ ^ %当t的标号为vs时，终止调整过程: e3 f- w: b* L7 A9 ^
if s(t) == 13 c8 u  R: ~9 g7 [) Z: q
      break;
) @5 I+ g& G  c0 U9 Y end
9 `8 p8 Q+ Z2 l2 p0 W& L* D t = s(t);" \8 o1 f7 Y: g; e( s
end
6 e1 ?' q  p# o+ `; `( I2 S %如果最大流量达到预定的流量值9 F. ?3 V5 S& u. Z/ F, O
if pd
* o* t8 N; k8 i. l       break;
( O8 c' I! m$ ~6 l! @, R* o end
7 f) ^) I6 l, C3 F$ T %计算最大流量0 f6 u6 t  S  N  k1 m) A4 U. B
wf = 0;
' D/ @& C. u* t+ Z0 V: t  Z for j = 1:n
+ b! E& W3 M4 A5 B+ c       wf = wf + f(1,j);
+ E% ]' b5 E2 Z' n- Z/ l end5 J# y, ]: ?2 K9 z; T! h/ K- B% X
end
, V1 P# d9 N) G, y5 c3 e%计算最小费用
, W8 ~' C( ^, q7 a: d5 ozwf = 0;
6 O- E8 t! j$ h/ [( B0 rfor i = 1:n1 g2 p' A: P7 l
for j = 1:n, B' p* a7 Y. E% A' ]- x& G* `" b
      zwf = zwf + b(i,j)*f(i,j);% I  l) k& a4 Y( {4 R% O
end3 Z( Q) c& j) }/ m# |2 y
end7 m) ?, b: F' s9 i2 w+ i) E/ G* X
%最小费用最大流
6 ^6 Y$ G/ A' ~/ P7 s1 I, If' d# `- Y/ y0 j( z* h- B2 \7 P* P) C) r
%最小费用最大流量
3 r$ I6 H/ ]+ R" A7 J5 ?wf- [& C8 R2 {( C7 N) d
%显示最小费用0 {$ Y! q! g% ]0 L
zwf
6 r6 t) c' M( S* g+ ~# H6 m1
9 U, D8 \* |  G8 @22 f9 d9 ]- e* m+ [) j9 o  N- |
3. Q/ `+ G; n3 \9 p1 s* P5 a
44 {& u% h0 f: p3 C( p0 S
5
( v% C1 v8 e: I# i+ d; ?63 A/ k6 v2 t1 n5 m3 U' v& r5 y, F
7$ b' Q+ S- l' c; I# e  {
8$ K! J. b0 p! f  ]' {/ v- E+ X. W0 }
94 s! W7 q6 H6 P1 Y! q' _1 }; z5 E
100 r, e+ I2 k: n  {* T$ C
11
4 J) `2 Z  {9 p% G' L1 p12" N8 d2 T' V0 d/ b
13$ S# R% ^, n2 x! A8 q) n
14
/ F* \  h7 h* o8 a! E+ `6 k15
/ S) D8 W7 F# R4 W1 M8 k16- k& ^, R4 i: I/ T
179 H6 [# d! n# s& z& m5 a$ \
18
( E+ o/ h4 A3 D/ u19, |! @0 u3 ~/ I  n
20: t6 ~) J0 |+ U% s1 ^, ]1 w8 l
21
% o( t8 G1 `6 m1 b: {22
' }$ `/ j1 w/ ~6 N230 [8 s5 I; |) T8 v, @
24
7 ]+ t  b: Z- s$ A4 L9 f& S25; V3 X3 R) o+ a+ d/ s$ v9 N& W
26+ ]( F, _& Q2 T: [
27
& R- a; V9 Y" h6 g& U) ^28: c# T2 {$ S/ v* b9 }/ R
29- }7 P$ ^! Y. [1 c. `
30
- v. l; \7 H7 J) `315 k' T/ U9 v3 h( a2 h& u
32
4 o5 O# R; `; d& `9 r( }; Z" P3 U0 F33+ J" i7 Z( V% p3 q, T1 S1 @
34% o& Q5 \% {$ n
35* Y7 M, j: Z! G. ]$ A
36
. L& E/ P' N  f37
0 Z+ P0 L9 O2 f, l38
: P5 j$ s' p9 S4 Y393 k2 D& y! Q- ~( `9 b4 o, }5 z2 E
40
0 O; J- @& ~1 Z+ g' O; m) ?41
' k% c6 d8 c  w1 n2 L9 ?42
% R# d  u7 ^' C( K6 |432 r8 K" ]8 C$ C' {4 G. |$ _
446 Q  l. G' ?: X& a
455 i7 |5 j5 e4 ?5 i- a7 K* {6 ?
46
/ L; C; [; R5 F+ L0 g, R9 O47
, Q5 a0 _, F9 Y; j7 f48' `) ?; X& ~4 H3 z) \/ p
49: g; e  ~, B. T4 W. P+ u$ R
50
3 I' N+ u9 b( e1 @& C9 \51: A9 D7 k$ M; v5 q" m. D0 b5 x9 H
52
  }' v# g' U  k1 a5 B53
4 Z1 ], y; L$ [: G+ C/ M) @4 t54/ h5 ^6 [  U4 [. z. M: K
55
7 u$ W9 [# L* [2 `# s& ~* F56
+ P0 }; g3 j8 ]2 b57& e; r# F& \. N, \- y
58
8 f( e' \8 h, U. J% z$ C/ T59; ]: l5 N( H- c; r8 o3 o
60
. `8 {7 y7 V# W& ~+ ~61) G  J: L7 W! Y% x: ^
62
) E; s0 n0 D, i& u63
# e8 ^. z7 V  s64
! r  h$ o" K3 a. g' B65
( `0 P. a% L4 k; C! s66
* W# k1 ]$ x8 H% X( J4 ~672 H2 Z$ f5 R* p
68
# Y; ~& \- B% c* K* n/ P( e1 F: r697 k0 J# l* T7 t: b8 o
701 R' G6 R  ]) y" L9 Z
714 B! A, Y8 ^3 j
72" F* ?4 _/ D9 h* D8 z
732 q- m' C' L$ v; ^7 [* a. |4 a
74
% ^, H1 q& [$ `4 m9 b75
. k  w1 b0 a1 W% p76
, Y) q; L; X% Y, U7 g$ h  p' q775 X% b0 L$ U( o: i  F: |, q- o+ U6 l, x- a
78
$ m; ~* w% B& y796 A! n1 G0 r5 D# a. E
802 k. x( |! |7 r% _& T* _. \0 d
81
2 N1 L2 O, ?9 l82
# ^+ I; y/ m# d7 D" @' ^( _83  i( V1 }0 \' Z
84" q9 U! [6 E& S- e; `2 S: N
85
7 J6 e* r! @  r; ?+ `" I86) w% q1 b: D+ Z" u- d' G: c
87- G' y; h0 y  {7 c
88' _% W: h- v: X- [/ [$ U4 U
89, @3 F; E3 e# W$ K
90
/ S7 |' a7 U+ i91
% n/ u5 \8 j! S$ b' o6 U92) ]" e& K8 G! q/ b3 b: z" [
932 w  V! u( B" M& T* j+ N2 Z% g. h
94
0 t& U6 n& Y7 c( {5 ]9 V95
% B9 F6 \2 ]' [4 b) {* i96
7 T: w2 h: Z, y& \6 R5 Z97
9 C- y4 z; |: V0 K98
% t' U2 F0 y, b; ?8 e# r. d% B5 Z# K( |99
& C0 h# t: W9 M8 S+ ^: G1004 o, Y5 B4 r- Z6 q+ q; B
101: E0 |2 B& B* l! ]% q1 K; T1 Y
102
: t+ f" s8 Z4 }$ B! R8 k" i! v103
' b. E2 r( R( Z3 _& |. W0 J- r1044 I. ]' P1 z$ ?
105& h9 D. ]3 u. T' X& e& _) O
106% v$ t3 Y$ t6 E- Y+ s/ P( A
1070 G; s% c/ g) h* A8 I% ]
108/ {* e4 T4 E" n8 r& b
109
* f0 n/ K+ |2 i. M6 L* g110
* f2 ?) L0 o3 \( H111+ y2 L8 T. s3 x  G
112- s0 ~; {+ m8 K4 e3 S/ g! F$ s4 v  A% G
113
  U) O) N( C5 i+ r" g: M7 V  W  z1146 m) d, X% u* W/ p4 a
115
  i/ `8 U# J3 G: K116/ z% S$ z: }7 g/ c2 L2 k& d
117
! a8 f: M" K4 k! S118
* @  w& u0 W. u$ Y( v119: }" Z# F# n- |6 s0 c
1207 K  K* O! T! k1 `- _" x& A* ^
121
7 r8 y+ q" d' |/ x+ k0 U3 p% Q% r122
9 i* H+ L+ w! \' u- O3 B% M123
9 a' C* ]& ?2 L# h124; _3 g. y, P' C6 a- k! I; P# I
125* m" T) [" x! m' K6 m7 t  h, G
1264 O: k. _$ a& }0 Y" n: b% v
127! v: K: A! B/ @  Y( b
128# M1 J( ~; W) r  p
129( P3 A4 A( k: U2 w. q. }7 d
130' O5 ?9 @0 _% T2 x+ R
131: G8 q; o: r& W# g
132
) e. d; c$ l& j+ s$ ~! T( y133
2 \- F+ `2 s( V4 X5 P7 i( @134
+ P3 h6 m' a/ T0 I8 t: ], c; s135
& ?4 [8 ~+ c; j2 B9 u136
& }+ |( i$ x0 U2 ]* a137
; v3 q  y" m' ?$ I( Q+ F1 M% C( h138
: X- a  \. k# O" k, y7 s4 O139
, x+ O  C: r$ d$ n# `) F$ l8 P140. @4 ]! v* `9 e% p' o" C& ]% I
匹配问题：8 A  ]) b9 o% k4 s1 c

" W7 m4 h, `% \最大匹配：
7 [( ~- _: S2 L% w& {: l( k/ q$ M* L. D* j

代码实现：

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