交通运输系统模型与方法——FW法配流

• Wardrop第一平衡——用户均衡：当网络达到平衡状态时， 每个 OD 对的各条被使用的径路具有相等而且最小的行驶时间； 没有被使用的径路的行驶时间大于或等于最小行驶时间。
• 以Beckman模型表达平衡的概念与定义
• Frank-Wolfe算法：在每次迭代中， 将目标函数 f(x)线性化， 通过解线性规划来求得可行下降方向， 进而沿此方向在可行域内作一维搜索以得到新的迭代点。

苏福尔斯网络不被视为现实网络。 然而， 该网络在许多出版物中都使用。 它有利于代码调试。

1.SiouxFalls_net:主要包含计算阻抗信息的相关量

2.SiouxFalls_node:主要为整个网络各个节点位置信息

3.SiouxFalls_trips:主要包含整个网络各个节点 OD 信息

4.Sioux-Falls-Network： 主要为网络图

5.阻抗函数:link travel time=freeflowtime

（1+B*(flow/capacity)^power）*

> https://github.com/bstabler/TransportationNetworks

# coding=utf-8

import math
import re
import torch
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import os
from datetime import datetime
import copy
import pylab as pl

matplotlib.rc('font', family='MicroSoft YaHei', weight='bold')

os.environ['KMP_DUPLICATE_LIB_OK'] = 'TRUE'
np.set_printoptions(suppress=True)

# 阻抗函数:link travel time=freeflowtime*（1+B*(flow/capacity)^power）
# xn1 = (1-a)*xn + a*yn
infinity = float("inf")
file_path_net = 'SiouxFalls_net_1.txt'
file_path_trips = 'SiouxFalls_trips_1.txt'


class Section():
    # 定义路段类
    def __init__(self, list, flow):
        self.init_node = list[0]
        self.term_node = list[1]
        self.capacity = float(list[2])
        self.length = int(list[3])
        self.freeflowtime = int(list[4])
        self.B = float(list[5])
        self.power = int(list[6])
        self.flow = flow
        f_c_p = math.pow(self.flow / self.capacity, self.power)
        self.travel_time = round(self.freeflowtime * (1 + self.B * f_c_p), 3)

    def update_flow_time(self, flow_new):  # 流量更新方法
        self.flow = flow_new
        f_c_p = math.pow(self.flow / self.capacity, self.power)
        self.travel_time = round(self.freeflowtime * (1 + self.B * f_c_p), 3)


def build_section(filePath):
    # 建立有路长、有时间花费的原始图，利用路网获取、更新、保存每个路段。
    section_list = []
    with open(filePath, 'r', encoding='utf-8') as f:
        for line in f:
            line_list = re.split(r'\t', line.strip(' ; \t \n'))
            section_i = Section(line_list, 0)
            section_list.append(section_i)
    # print(len(section_list))
    for sec in section_list:
        # print(sec.__dict__)
        pass
    return section_list


def build_pointpair(filePath):
    # 建立有OD点对之间流量需求的原始图，利用点对获取每个路段、路径的流量。
    time_list = []
    sim_list = []
    dou_list = []
    num = 0
    with open(filePath, 'r', encoding='utf-8') as f:
        for line in f:
            line_list = re.findall(r':( .*?);', line)
            for i in line_list:
                t = i.strip(' ')
                time_list.append(int(float(t)))
    for time in time_list:
        sim_list.append(time)
        num += 1
        while num >= 24:
            dou_list.append(sim_list)
            num = 0
            sim_list = []
    # print(dou_list)
    return dou_list


def update_section(section_list, flow_list):
    # 原地更新流量！！！已修改
    # 根据新的输入流量，更新成该状态下的路网（即路段列表）
    new_section_list = []
    for i in range(len(section_list)):
        flow_tmp = section_list[i].flow
        section_list[i].update_flow_time(flow_list[i])
        sec_tmp = copy.deepcopy(section_list)  # 临时路段
        new_section_list.append(sec_tmp[i])
        section_list[i].update_flow_time(flow_tmp)
    for sec in new_section_list:
        # print(sec.__dict__)
        pass
    return new_section_list


def get_flow(section_list):
    # 路网各路段流量
    flow_list = []
    for sec in section_list:
        flow_list.append(sec.flow)
    # print(flow_list)
    return flow_list


def get_time(section_list):
    # 路网各路段阻抗时间
    time_list = []
    for sec in section_list:
        time_list.append(sec.travel_time)
    # print(time_list)
    return time_list


def cal_cost(section_list):
    # 计算花费
    cost = 0
    for sec in section_list:
        cost += sec.flow * sec.travel_time
    # print(cost)
    return cost


def search_best_way(section_list, o, d):
    # 在上一迭代流量的阻抗下，Floyd算法寻找最短路
    way_list = [o]
    dis = []
    for hang in range(24):
        hang_jl_list = []
        for gezi in range(24):
            if hang == gezi:
                hang_jl_list.append(0)
                continue
            for sec in section_list:
                if hang == int(sec.init_node) - 1 and gezi == int(sec.term_node) - 1:
                    hang_jl_list.append(sec.travel_time)
            if len(hang_jl_list) < gezi + 1:
                hang_jl_list.append(infinity)
        dis.append(hang_jl_list)
    # print(dis)
    prior = [[0 for i in range(1, 25)] for j in range(1, 25)]  # 初始化一个空矩阵用来存放所经过的节点
    for i in range(24):
        for j in range(24):
            prior[i][j] = j + 1
    for m in range(24):  # m相当于是中间点
        for i in range(24):  # i相当于是起始点
            for j in range(24):  # j相当于是终点
                if (dis[i][j] > dis[i][m] + dis[m][j]):

                    dis[i][j] = dis[i][m] + dis[m][j]
                    prior[i][j] = prior[i][m]
                else:
                    continue
    # print(prior)
    if prior[o - 1][d - 1] == d:
        way_list.append(d)
    else:
        x = prior[o - 1][d - 1]
        mid = [x]  # 创建一个空列表装载路径中间点
        while (prior[x - 1][d - 1] != d):
            x = prior[x - 1][d - 1]
            mid.append(x)
        way_list.extend(mid)
        way_list.append(d)
    # print(dis)
    # print(way_list, '==>总路程：', dis[o - 1][d - 1])
    return way_list


def product_flow_list(way_list, section_list, flow):
    # 根据最短路全有全无，分配路径流量至各路段
    sec_tmp = copy.deepcopy(section_list)
    for sec in sec_tmp:
        change = 0
        for i in range(len(way_list) - 1):
            if sec.init_node == str(way_list[i]) and sec.term_node == str(way_list[i + 1]):
                sec.update_flow_time(flow)
                change = 1
        if change == 0:
            sec.update_flow_time(0)
        # print(sec.__dict__)
    new_flow = get_flow(sec_tmp)
    # print(new_flow)
    return new_flow


def graph_new_flow(section_list):
    #将全部OD点对的流量分配至某一路网
    flow_new = flow_init
    for o in range(1, 25):
        for d in range(1, 25):
            way_list = search_best_way(section_list, o, d)  # 根据上次迭代的阻抗，获得最短路径
            flow_new = np.array(flow_new) + np.array(
                product_flow_list(way_list, section_list, flow_given[o - 1][d - 1]))  # 获得新的流量分配
    flow_new.tolist()
    # print(flow_new)
    return flow_new


def find_best_alpha(flow_1, flow_2):
    c = infinity

    for alpha in np.arange(0.01, 1, 0.01):
        flow_y_x = (np.array(flow_2) - np.array(flow_1)).tolist()

        flow_a_y_x = [alpha * i for i in flow_y_x]
        flow_n1 = (np.array(flow_1) + np.array(flow_a_y_x)).tolist()
        section_list_1 = update_section(section_list_0, flow_n1)  # 获得新的路段对象
        time_list = get_time(section_list_1)  # 获得新的路段阻抗

        all = zip(flow_y_x, time_list)
        cost_n1 = 0
        for i, j in all:
            cost_n1 += i * j
        if abs(cost_n1) < c:
            c = abs(cost_n1)
            a = alpha
    print(a)
    # print(c)
    return a


# ---------------------主程序----------------------------
print(datetime.now().strftime('%Y-%m-%d %H:%M:%S'))

flow_given = build_pointpair(file_path_trips)

flow_init = [0 for i in range(76)]  # 获得初始流量分配
# print(flow_init)
section_list_0 = build_section(file_path_net)  # 获得路段对象
get_flow(section_list_0)
get_time(section_list_0)
cal_cost(section_list_0)

flow_1 = graph_new_flow(section_list_0)
print(flow_1)
n = 1

cost_list = []
for i in range(10):

    print('第%d次迭代' % (i + 1))
    section_list_1 = update_section(section_list_0, flow_1)  # 获得新的路段对象
    time_1 = get_time(section_list_1)  # 获得新的路段阻抗
    cost_1 = cal_cost(section_list_1)
    # print(flow_1)  # flow_1为xn
    # print(time_1)  # time_1为t(xn)
    # print(cost_1)

    # 搜索可行方向
    flow_2 = graph_new_flow(section_list_1)
    way_list = search_best_way(section_list_1, 2, 8)
    # print(flow_2)  # flow_2为yn

    alpha = find_best_alpha(flow_1, flow_2)  # 寻找迭代步长
    # alpha = 0.618     #0.618法，可进行尝试

    x = [alpha * i for i in flow_1]
    y = [alpha * i for i in flow_2]
    flow_n1 = (np.array(flow_1) + np.array(y) - np.array(x)).tolist()
    print(flow_n1)  # flow_n1为xn+1

    section_list_2 = update_section(section_list_1, flow_n1)  # 获得新的路段对象
    cost_n1 = cal_cost(section_list_2)
    for sec in section_list_1:
        # print(sec.__dict__)
        pass
    # print(get_time(section_list_2))
    print(cost_n1)
    cost_list.append(cost_n1)

    # 迭代，令n=n+1
    section_list_0 = update_section(section_list_2, flow_n1)
    flow_1 = flow_n1
    # print(get_flow(section_list_0))

    print('################################################################')

print(cost_list)
aaa = [i for i in range(len(cost_list))]
plt.plot(aaa, cost_list)

ax = pl.gca()
ax.get_yaxis().get_major_formatter().set_scientific(False)

plt.xlabel('迭代次数')
plt.ylabel('路网总花费')
plt.title('平衡Frank-Wolfe算法30步')
plt.show()

str = '\n'
f = open("baojy1.txt", "w")
f.write(float.join(flow_n1))
f.write(float.join(cost_n1))
f.close()