图与网络模型及方法(二)

        在我图与网络模型及方法(一)的文章中介绍了图论相关模型与方法的原理,这篇文章将详细介绍一些经典的应用问题。

        python中有个专门处理图、网络相关问题的库:networkx,以下部分应用代码实现部分给出直接调库的代码实现方法。

一.匹配问题

 

代码示例:

import numpy as np
import networkx as nx
# 匹配问题 二分图
worker_n = 5  # 工人数目
work_n = 5  # 工作数目
a = np.array([[3, 5, 5, 4, 1], [2, 2, 0, 2, 2], [2, 4, 4, 1, 0],
              [0, 2, 2, 1, 0], [1, 2, 1, 3, 3]])
b = np.zeros((10, 10))
b[0:worker_n, work_n:] = a
G = nx.Graph(b)
s0 = nx.max_weight_matching(G)  # 返回值为(人员、工作的集合)
s = [sorted(w) for w in s0]
L1 = np.array([x[0] for x in s])  # 人员编号
L2 = np.array([x[1] for x in s])  # 工作编号
c = []
for i in range(len(L1)):
    c.append(a[L1[i]][L2[i] - worker_n])
d = sum(c)  # 计算总的效益
print("工作分配对应关系为:")
print("人员编号:", L1 + 1)
print("工作编号:", L2 - worker_n + 1)
print("总效益:", d)
print('-------------------------------------------------')

二.Euler图和Hamilton图

2.1 基本概念

定义

        经过G的每条边的迹叫做G 的 Euler 迹;闭的 Euler 迹叫做 Euler回路或E回路;含 Euler回路的图叫做Euler图。直观地讲,Euler图就是从一顶点出发每边恰通过一次能回到出发点的那种图,即不重复地行遍所有的边再回到出发点。

        Euler图最初是为了解决戈尼斯堡七桥问题

 

 定义

        包含G 的每个顶点的轨叫做 Hamilton(哈密顿)轨;闭的 Hamilton 轨叫做Hamilton 圈或 H 圈;含 Hamilton 圈的图叫做 Hamilton 图。

        直观地讲,Hamilton 图就是从一顶点出发每顶点恰通过一次能回到出发点的那种图,即不重复地行遍所有的顶点再回到出发点

2.2  Euler回路的Fleury算法

 python实现代码:Python实现Fleury算法_Paul_barnabas的博客-CSDN博客_fleury算法python

2.3 应用

2.3.1 邮递员问题

中国邮递员问题

        一位邮递员从邮局选好邮件去投递,然后返回邮局,当然他必须经过他负责投递的每条街道至少一次,为他设计一条投递路线,使得他行程最短。

         上述中国邮递员问题的数学模型是:在一个赋权连通图上求一个含所有边的回路,且使此回路的权最小。

        显然,若此连通赋权图是 Euler 图,则可用 Fleury 算法求 Euler 回路,此回路即为所求。

        对于非 Euler 图,1973 年,Edmonds 和 Johnson 给出下面的解法

多邮递员问题(kpp)

        邮局有 k(k ≥ 2) 位投递员,同时投递信件,全城街道都要投递,完成任务返回邮 局,如何分配投递路线,使得完成投递任务的时间最早?

 

 

参考:中国邮递员问题的深入剖析与算法实现(附例题及MATLAB、LINGO代码)_第二号的博客-CSDN博客_中国邮递员问题例题

2.3.2 旅行商问题TSP

 旅行商问题TSP

         一名推销员准备前往若干城市推销产品,然后回到他的出发地。如何为他设计一条最短的旅行路线(从驻地出发,经过每个城市恰好一次,最后返回驻地)?这个问题称为旅行商问题。用图论的术语说,就是在一个赋权完全图中,找出一个有最小权的 Hamilton 圈。称这种圈为最优圈。

        与最短路问题及连线问题相反,目前还没有求解旅行 商问题的有效算法。所以希望有一个方法以获得相当好(但不一定最优)的解。 

function main 
clc,clear 
global a 
a=zeros(6); 
a(1,2)=56;a(1,3)=35;a(1,4)=21;a(1,5)=51;a(1,6)=60; 
a(2,3)=21;a(2,4)=57;a(2,5)=78;a(2,6)=70; 
a(3,4)=36;a(3,5)=68;a(3,6)=68; a(4,5)=51;a(4,6)=61; 
a(5,6)=13; a=a+a'; L=size(a,1); 
c1=[5 1:4 6]; 
[circle,long]=modifycircle(c1,L); 
c2=[5 6 1:4];%改变初始圈,该算法的最后一个顶点不动
[circle2,long2]=modifycircle(c2,L); 
if long2<long 
    long=long2; 
    circle=circle2; 
end 
circle,long 
%******************************************* 
%修改圈的子函数
%******************************************* 
function [circle,long]=modifycircle(c1,L); 
global a 
flag=1; 
while flag>0 
    flag=0; 
    for m=1:L-3 
        for n=m+2:L-1 
            if a(c1(m),c1(n))+a(c1(m+1),c1(n+1))<... a(c1(m),c1(m+1))+a(c1(n),c1(n+1)) 
                flag=1; 
                c1(m+1:n)=c1(n:-1:m+1); 
            end 
        end 
    end 
end 
long=a(c1(1),c1(L)); 
for i=1:L-1 
    long=long+a(c1(i),c1(i+1)); 
end 
circle=c1;

 其他求解方式有模拟退火、遗传算法等现代智能算法,如

import random
import math
import numpy as np
import matplotlib.pyplot as plt
# 随机初始化城市坐标
number_of_citys = 100
citys = []
for i in range(number_of_citys):
    citys.append([random.randint(1, 100), random.randint(1, 100)])
citys = np.array(citys)

# 由城市的坐标计算距离矩阵
distance = np.zeros((number_of_citys, number_of_citys))
for i in range(number_of_citys):
    for j in range(number_of_citys):
        distance[i][j] = math.sqrt((citys[i][0] - citys[j][0])**2 +
                                   (citys[i][1] - citys[j][1])**2)
# 初始化参数
iter1 = 2000  # 外循环迭代次数
T0 = 100000  # 初始温度
Tf = 1  # 截止温度
T_now = T0  # 当前温度
alpha = 0.95  # 温度更新因子
iter2 = 10  # 内循环迭代次数
fbest = 0  # 最佳距离

# 初始化初解
x = []
for i in range(number_of_citys):
    x.append(i)
np.random.shuffle(x)
x = np.array(x)
fbest = float('inf')
xbest = x.copy()
f_now = fbest
x_now = xbest.copy()

for m in range(iter1):
    for k in range(iter2):
        # 生成新解
        x1 = [0 for q in range(number_of_citys)]
        n1 = random.randint(0, number_of_citys - 1)
        n2 = random.randint(0, number_of_citys - 1)  # 对路径进行“微扰”
        n = [n1, n2]
        n.sort()
        n1, n2 = n
        # n1 == 0 单独写
        if n1 > 0:
            x1[0:n1] = x_now[0:n1]
            x1[n1:n2 + 1] = x_now[n2:n1 - 1:-1]
            x1[n2 + 1:number_of_citys] = x_now[n2 + 1:number_of_citys]
        else:
            x1[0:n1] = x_now[0:n1]
            x1[n1:n2 + 1] = x_now[n2::-1]
            x1[n2 + 1:number_of_citys] = x_now[n2 + 1:number_of_citys]
        s = 0
        for j in range(len(x1) - 1):
            s = s + distance[x1[j]][x1[j + 1]]
        s += distance[x1[-1]][x1[0]]
        # 判断是否有更新解
        if s <= f_now:
            f_now = s
            x_now = x1.copy()
        else:
            if random.random() < math.exp((f_now - s) / T_now):
                f_now = s
                x_now = x1.copy()
        if s < fbest:  # 历史最优解
            fbest = s
            xbest = x1.copy()
    T_now = T_now * alpha

print("最佳路线:", xbest)
print("最佳距离:", fbest)

plt.title('SA_TSP')
plt.xlabel('x')
plt.ylabel('y')
plt.plot(citys[..., 0], citys[..., 1], 'ob', ms=3)
citys_x = []
citys_y = []
for i in xbest:
    citys_x.append(citys[i][0])
    citys_y.append(citys[i][1])
plt.plot(citys_x, citys_y)
plt.plot([citys[xbest[-1]][0], citys[xbest[0]][0]],
         [citys[xbest[-1]][1], citys[xbest[0]][1]],
         ms=2)
plt.show()

结果:

 三.最大流问题

3.1 最大流问题的数学描述

3.1.1 网络中的流

定义

 在以V 为节点集, A 为弧集的有向图G = (V, A) 上定义如下的权函数:

(i) L : A → R 为孤上的权函数,弧 (i, j)∈ A 对应的权 L(i, j) 记为lij ,称为孤 (i, j) 的容量下界(lower bound);

(ii)U : A → R 为弧上的权函数,弧(i, j)∈ A对应的权U(i, j) 记为uij ,称为孤 (i, j) 的容量上界,或直接称为容量(capacity);

(iii) D :V → R 为顶点上的权函数,节点i ∈V 对应的权 D(i) 记为 di ,称为顶点i 的供需量(supply/demand);

此时所构成的网络称为流网络,可以记为

N = (V, A,L,U,D)

 

 定义

 

 

 

 

3.1.2 最大流问题

 

定义

        如果一个矩阵 A 的任何子方阵的行列式的值都等于0,1或 −1,则称 A 是全幺模的(totally unimodular TU,又译为全单位模的),或称A是全幺模矩阵。

定理

        (整流定理) 最大流问题所对应的约束矩阵是全幺模矩阵。若所有弧容量均为正整数,则问题的最优解为整数解

        最大流问题是一个特殊的线性规划问题。我们将会看到利用图的特点,解决这个问 题的方法较之线性规划的一般方法要方便、直观得多。 

3.1.3 单源和单汇运输网络

 3.2 最大流和最小割关系

 3.3 最大流的一种算法——标号法

 

clc,clear 
u(1,2)=1;u(1,3)=1;u(1,4)=2;u(2,3)=1;u(2,5)=2; 
u(3,5)=1;u(4,3)=3;u(4,5)=3; 
f(1,2)=1;f(1,3)=0;f(1,4)=1;f(2,3)=0;f(2,5)=1; 
f(3,5)=1;f(4,3)=1;f(4,5)=0; 
n=length(u);list=[];maxf(n)=1; 
while maxf(n)>0 
    maxf=zeros(1,n);pred=zeros(1,n); 
    list=1;record=list;maxf(1)=inf; 
 %list是未检查邻接点的标号点,record是已标号点
    while (~isempty(list))&(maxf(n)==0) 
        flag=list(1);list(1)=[]; 
        label1= find(u(flag,:)-f(flag,:)); 
        label1=setdiff(label1,record); 
        list=union(list,label1); 
        pred(label1)=flag; 
        maxf(label1)=min(maxf(flag),u(flag,label1)... -f(flag,label1)); 
        record=union(record,label1); 
        label2=find(f(:,flag)); 
        label2=label2'; 
        label2=setdiff(label2,record); 
        list=union(list,label2); 
        pred(label2)=-flag; 
        record=union(record,label2); 
    end 
    if maxf(n)>0 
        v2=n; v1=pred(v2); 
        while v2~=1 
            if v1>0 
                f(v1,v2)=f(v1,v2)+maxf(n); 
            else 
                v1=abs(v1); 
            f(v2,v1)=f(v2,v1)-maxf(n); 
            end 
            v2=v1; v1=pred(v2); 
        end 
    end 
end 
f

或者通过使用networkx中maximum_flow函数直接求解最大流问题

import networkx as nx
# 最大流问题
L = [(1, 2, 5), (1, 3, 3), (2, 4, 2), (3, 2, 1), (3, 5, 4), (4, 3, 1),
     (4, 5, 3), (5, 6, 5)]
G = nx.DiGraph()
for k in range(len(L)):
    G.add_edge(L[k][0] - 1, L[k][1] - 1, capacity=L[k][2])
value, flow_dict = nx.maximum_flow(G, 0, 5)
print("最大流的流量为:", value)
print("最大流为:", flow_dict)  # 最大流情形下, 每条边实际流量

四.最小费用流及其求法

4.1 最小费用流

        上面我们介绍了一个网络上最短路以及最大流的算法,但是还没有考虑到网络上流的费用问题,在许多实际问题中,费用的因素很重要。例如,在运输问题中,人们总是希望在完成运输任务的同时,寻求一个使总的运输费用最小的运输方案。这就是下面要 介绍的最小费用流问题。

 4.2 求最小费用流的一种方法——迭代法

 

 matlab代码如下:

function mainexample19 
clear;clc; 
global M num 
c=zeros(6);u=zeros(6); 
c(1,2)=2;c(1,4)=8;c(2,3)=2;c(2,4)=5; 
c(3,4)=1;c(3,6)=6;c(4,5)=3;c(5,3)=4;c(5,6)=7; 
u(1,2)=8;u(1,4)=7;u(2,3)=9;u(2,4)=5; 
u(3,4)=2;u(3,6)=5;u(4,5)=9;u(5,3)=6;u(5,6)=10; 
num=size(u,1);M=sum(sum(u))*num^2; 
[f,val]=mincostmaxflow(u,c) 
%求最短路径函数
function path=floydpath(w); 
global M num 
w=w+((w==0)-eye(num))*M; 
p=zeros(num); 
for k=1:num 
 for i=1:num 
 for j=1:num 
 if w(i,j)>w(i,k)+w(k,j) 
 w(i,j)=w(i,k)+w(k,j); 
 p(i,j)=k; 
 end 
 end 
 end 
end 
if w(1,num) ==M 
 path=[]; 
else 
 path=zeros(num); 
 s=1;t=num;m=p(s,t); 
 while ~isempty(m) 
 if m(1) 
 s=[s,m(1)];t=[t,t(1)];t(1)=m(1); 
 m(1)=[];m=[p(s(1),t(1)),m,p(s(end),t(end))]; 
 else 
 path(s(1),t(1))=1;s(1)=[];m(1)=[];t(1)=[]; 
 end 
 end 
end 
%最小费用最大流函数
function [flow,val]=mincostmaxflow(rongliang,cost,flowvalue); 
%第一个参数:容量矩阵;第二个参数:费用矩阵;
%前两个参数必须在不通路处置零
%第三个参数:指定容量值(可以不写,表示求最小费用最大流)
%返回值 flow 为可行流矩阵,val 为最小费用值
global M 
flow=zeros(size(rongliang));allflow=sum(flow(1,:)); 
if nargin<3 
 flowvalue=M; 
end 
while allflow<flowvalue 
 w=(flow<rongliang).*cost-((flow>0).*cost)'; 
 path=floydpath(w);%调用 floydpath 函数
 if isempty(path) 
 val=sum(sum(flow.*cost)); 
 return; 
 end 
 theta=min(min(path.*(rongliang-flow)+(path.*(rongliang-flow)==0).*M)); 
 theta=min([min(path'.*flow+(path'.*flow==0).*M),theta]); 
 flow=flow+(rongliang>0).*(path-path').*theta; 
 allflow=sum(flow(1,:)); 
end 
val=sum(sum(flow.*cost));

调用python库函数方法:

# 最小费用流
L = [(1, 2, 5, 3), (1, 3, 3, 6), (2, 4, 2, 8), (3, 2, 1, 2), (3, 5, 4, 2),
     (4, 3, 1, 1), (4, 5, 3, 4), (4, 6, 2, 10), (5, 6, 5, 2)]
G = nx.DiGraph()
for k in range(len(L)):
    G.add_edge(L[k][0] - 1, L[k][1] - 1, capacity=L[k][2], weight=L[k][3])
mincostFlow = nx.max_flow_min_cost(G, 0, 5)
print("所求流为:", mincostFlow)
mincost = nx.cost_of_flow(G, mincostFlow)
print("最小费用为:", mincost)

  • 1
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值