基于Frank Wolfe算法，求解交通分配UE模型（Python&NetWorkX）

#!/usr/bin/env python3 # -*- coding: utf-8 -*- # @Time : 2024/3/9 23:42 # @Author : zy_chai # @FileName: UE.py # @Software: PyCharm import pandas as pd import numpy as np import networkx as nx from scipy.optimize import minimize_scalar def BPR(FFT, flow, capacity, alpha=0.15, beta=4.0): return FFT * (1 + alpha * (flow / capacity) ** beta) def all_none_initialize(G, od_df): # 这个是实际分配，在零流图上按全有全无分配，用于得到flow_real # 这个函数仅使用一次，用于初始化 for _, od_data in od_df.iterrows(): source = od_data["o"] target = od_data["d"] demand = od_data["demand"] # 计算最短路径 shortest_path = nx.shortest_path(G, source=source, target=target, weight="weight") # 更新路径上的流量 for i in range(len(shortest_path) - 1): u = shortest_path[i] v = shortest_path[i + 1] G[u][v]['flow_real'] += demand # 初始化流量后，更新阻抗 for _, _, data in G.edges(data=True): data['weight'] = BPR(data['FFT'], data['flow_real'], data['Capacity']) def all_none_temp(G, od_df): # 这个是虚拟分配，用于得到f_temp nx.set_edge_attributes(G, 0, 'flow_temp') for _, od_data in od_df.iterrows(): # 每次更新都得读OD，后面尝试优化这个 source = od_data["o"] target = od_data["d"] demand = od_data["demand"] # 计算最短路径 shortest_path = nx.shortest_path(G, source=source, target=target, weight="weight") # 更新路径上的流量 for i in range(len(shortest_path) - 1): u = shortest_path[i] v = shortest_path[i + 1] # 更新流量 G[u][v]['flow_temp'] += demand def get_descent(G): for _, _, data in G.edges(data=True): data['descent'] = data['flow_temp'] - data['flow_real'] def objective_function(temp_step, G): s, alpha, beta = 0, 0.15, 4.0 for _, _, data in G.edges(data=True): x = data['flow_real'] + temp_step * data['descent'] s += data["FFT"] * (x + alpha * data["Capacity"] / (beta + 1) * (x / data["Capacity"]) ** (beta + 1)) return s def update_flow_real(G): # 这个函数实际调整流量，用于调整flow_real，并更新weight best_step = get_best_step(G) # 获取最优步长 for _, _, data in G.edges(data=True): # 调整流量，更新路阻 data['flow_real'] += best_step * data["descent"] data['weight'] = BPR(data['FFT'], data['flow_real'], data['Capacity']) def get_best_step(G, tolerance=1e-4): result = minimize_scalar(objective_function, args=(G,), bounds=(0, 1), method='bounded', tol=tolerance) return result.x def build_network(Link_path): # 读取边数据 links_df = pd.read_csv(Link_path) # 需要注意使用from_pandas_edge，其读取的边的顺序和csv中边的顺序有差异 G = nx.from_pandas_edgelist(links_df, source='O', target='D', edge_attr=['FFT', 'Capacity'], create_using=nx.DiGraph()) nx.set_edge_attributes(G, 0, 'flow_temp') nx.set_edge_attributes(G, 0, 'flow_real') nx.set_edge_attributes(G, 0, 'descent') nx.set_edge_attributes(G, nx.get_edge_attributes(G, "FFT"), 'weight') return G def main(): G = build_network("Link.csv") # 构建路网 od_df = pd.read_csv("OD.csv") # 获取OD需求情况 all_none_initialize(G, od_df) # 初始化路网流量 print("初始化流量", list(nx.get_edge_attributes(G, 'flow_real').values())) epoch = 0 # 记录迭代次数 err, max_err = 1, 1e-4 # 分别代表初始值、最大容许误差 f_list_old = np.array(list(nx.get_edge_attributes(G, 'flow_real').values())) while err > max_err: epoch += 1 all_none_temp(G, od_df) # 全有全无分配，得到flow_temp get_descent(G) # 计算梯度，即flow_temp-flow_real update_flow_real(G) # 先是一维搜索获取最优步长，再调整流量，更新路阻 # 计算并更新误差err f_list_new = np.array(list(nx.get_edge_attributes(G, 'flow_real').values())) # 这个变量是新的路网流量列表 d = np.sum((f_list_new - f_list_old) ** 2) err = np.sqrt(d) / np.sum(f_list_old) f_list_old = f_list_new print("均衡流量", list(nx.get_edge_attributes(G, 'flow_real').values())) print("迭代次数", epoch) # 导出网络均衡流量 df = nx.to_pandas_edgelist(G) df = df[["source", "target", "flow_real"]].sort_values(by=["source", "target"]) df.to_csv("网络均衡结果.csv", index=False) if __name__ == '__main__': main()