数学建模--图与网络（2）
最大流：https://img-blog.csdn.net/20170902231642951?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXFfMzQ4NjExMDI=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast注意Matlab 中的最大流问题是必须是单源和单汇问题，因此这里需要构建虚拟的源点S和汇点G。clc,clear
a = zeros(9,9);
a(1,) = ;
a(2,) = ;
a(3,) = ;
a(4,) = ;
a(,9) = ;
a = sparse(a);
= graphmaxflow(a,1,9)
最大流拓展：最小费用最大流仅仅是在求得最大流之后进行对最小费用的求解：https://img-blog.csdn.net/20170903005142549?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXFfMzQ4NjExMDI=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast
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); = graphmaxflow(a,1,5)
最大流量为11自定义Matlab代码：
最小费用求解
Lingo:
model:
sets:
nodes/s,1,2,3,t/:d;
arcs(nodes,nodes)/s 1,s 2,1 3,1 t,2 1,2 3,3 t/:b,c,f;
endsets
data:
b = 4 1 6 1 2 3 2;
c = 10 8 2 7 5 10 4;
d = 11 0 0 0 -11;
enddata
n = @size(nodes);
min = @sum(arcs:b*f);
@for(nodes(i):@sum(arcs(i,j):f(i,j))-@sum(arcs(j,i):f(j,i)) = d(i));
@for(arcs:@bnd(0,f,c));
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Matlab实现：
n = 5;
%弧容量
a = zeros(5);
a(1,) = ;
a(2,) = ;
a(3,) = ;
a(4,5) = 4;
C = a;
%C = ;
%弧上单元的费用
a(1,) = ;
a(2,) = ;
a(3,) = ;
a(4,5) = 2;
b = a;
%b = ;
%wf表示最大流量。wf0表示预定的流量值
wf = 0;
wf0 = Inf;
%取初始可行流f为零流
for i = 1:n
for j = 1:n
f(i,j) = 0;
end
end
while (1)
%构造有向赋权图
for i = 1:n
for j = 1:n
if j~=i
a(i,j) = Inf;
end
end
end
for i = 1:n
for j = 1:n
if C(i,j) > 0 && f(i,j) == 0
a(i,j) = b(i,j);
elseif C(i,j) > 0 && f(i,j) == C(i,j)
a(j,i) = -b(i,j);
elseif C(i,j) > 0
a(i,j) = b(i,j);
a(j,i) = -b(i,j);
end
end
end
%使用ford算法求最短路，赋初值
for i = 2:n
p(i) = Inf;
s(i) = i;
end
%求有向赋权图vs到vt的最短路,赋初值
for k = 1:n
pd = 1;
for i = 2:n
for j = 1:n
if p(i) > p(j) + a(j,i)
p(i) = p(j) + a(j,i);
s(i) = j;
pd = 0;
end
end
end
%求最短路的Ford算法结束
if pd
break;
end
end
%不存在vs到vt的最短路，算法终止，注意在求最小费用最大流时构造有向赋权图中不会含有负权回路，不会出现k = n
if p(n) == Inf
break;
end
%进入调整过程,dvt表示调整量
dvt = Inf;
dvtt = Inf;
t = n;
while(1)
if a(s(t),t) > 0
%前向弧调整量
dvtt = C(s(t),t)-f(s(t),t);
%后向弧调整量
elseif a(s(t),t) < 0
dvtt = f(t,s(t));
end
if dvt > dvtt
dvt = dvtt;
end
%当t的标号为vs时，终止计算调整量
if s(t) == 1
break;
end
%继续调整前一段弧上的流f
t = s(t);
end
pd = 0;
%如果最大流量大于或等于预定的流量值
if wf + dvt >= wf0
dvt = wf0 - wf;
pd = 1;
end
t = n;
%调整过程
while(1)
if a(s(t),t) > 0
%前向弧调整
f(s(t),t) = f(s(t),t) + dvt;
elseif a(s(t),t) < 0
%后向弧调整
f(t,s(t)) = f(t,s(t)) - dvt;
end
%当t的标号为vs时，终止调整过程
if s(t) == 1
break;
end
t = s(t);
end
%如果最大流量达到预定的流量值
if pd
break;
end
%计算最大流量
wf = 0;
for j = 1:n
wf = wf + f(1,j);
end
end
%计算最小费用
zwf = 0;
for i = 1:n
for j = 1:n
zwf = zwf + b(i,j)*f(i,j);
end
end
%最小费用最大流
f
%最小费用最大流量
wf
%显示最小费用
zwf
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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
匹配问题：
最大匹配：
https://img-blog.csdn.net/20170904114649233?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXFfMzQ4NjExMDI=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast
代码实现：
m = 5;
n = 5;
A = ;
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 = ;
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
页:
[1]