点击蓝字 关注我们
珍惜当下，把握现在。
——王婷
引言：从本科的课程《交通规划》开始，我们就学习了UE平衡和交通分配(Traffic Assignment)的概念。而Frank-wolfe算法作为求解用户平衡交通分配问题的基本算法，是学习交通分配的重中之重，也是学习交通类优化算法的重点内容。本文介绍了用户平衡和Frank-wolfe算法的基本原理，并给出了非常详细的编程实现过程。文中的程序既可当作一次编程练习，也可以作为编写这一类算法的程序框架，为今后用程序求解线性约束下的优化问题提供了思路。
(注：文中重要的公式和算法用红色五角星★标明，以方便查询)
1 背景知识补充
在讲解Frank-Wolfe算法的基本原理之前，先补充一些背景知识。
交通分配(Traffic Assignment)：是将预测出来的OD(Origin-Destination)需求按照一定的规则符合实际地分配到网络中的各条道路上，并求出各条道路的交通流量。交通分配是四阶段交通规划方法中的最后一个阶段。
BPR函数：由美国联邦公路局(BPR—Bureau of public road)开发的阻抗函数，用来计算路段阻抗。
Beckmann模型：用户平衡交通分配问题的数学规划问题。
间隙函数(Gap Function)：在UE交通分配问题的求解过程中，用于评价解的精度。
2 Frank-Wolfe算法基本原理
3 在UE交通分配问题中的应用
4 Frank-Wolfe算法详细过程
全有全无算法是UE交通分配算法中的重要步骤，可以构造成函数供Frank-Wolfe算法调用。具体算法步骤如下：
采用Frank-wolf算法求解UE交通分配问题的详细算法步骤如下：
5 Frank-Wolfe算法代码详解
考虑算法的计算量和网络规模，本文以效率更高的C#语言为例进行代码讲解。编程语言是相通的，擅长使用Python等其他语言的读者可以自行切换。
5.1 代码文件Network的编写
(1)准备工作
首先在项目中新建名为“Network.cs”的代码文件(主函数main在默认的“Program.cs”文件中)，Network文件开头引用如下命名空间。
using System;using System.IO;using System.Collections.Generic;(2)创建描述网络基本属性的类
构建三个类，分别对应网络的节点(CNode)、路段(CLink)和OD对(COrigin)，每读取一个网络的节点，就相当于在类CNode中创建一个新的对象。这些类是网络数据结构的基础，在类中分别声明节点、路段和OD对的属性，具体属性见代码的注释部分。
// 节点public class CNode{ public int ID; // 节点的编号，从零开始编号 public double PositionX; // 节点的X坐标 public double PositionY; // 节点的Y坐标 public int Origin_ID = -1; // 节点对应的起点编号，-1表示不是起点 public List<int> IncomingLink = new List<int>(); // 进入节点的路段编号集合 public List<int> OutgoingLink = new List<int>(); // 离开节点的路段编号集合 }// 路段public class CLink{ public int ID; // 路段的编号，从零开始编号 public CNode pInNode; // 路段的起节点 public CNode pOutNode; // 路段的终节点 public double FreeFlowTravelTime; // 自由流走行时间 public double TravelTime; // 走行时间 public double Capacity; // 路段通行能力 public double Alpha = 0.15; // BPR函数参数，一般取0.15 public double Power = 4.0; // BPR函数参数，一般取4.0}// OD对public class COrigin{ public int ID; // 起点的编号，从零开始编号 public CNode pOriginNode; // 起点对应的节点 public List<int> DestinationNode = new List<int>(); // OD对，只记录有需求的OD public List<double> ODDemand = new List<double>(); // OD需求，只记录有需求的OD}(3)创建算法中变量和函数的类
在这个类中，定义了Frank-Wolfe算法中的一些重要的变量作为类CNetwork的属性，如路段流量、路段走行时间、UE误差等。此外，还在类中构造了实现算法的函数，这些构造函数是算法的主体。
创建类，并定义类属性的代码如下。
需要注意的是，后续所有的构造函数都在类CNetwork中编写。因此，这些函数将代替以下代码中“……”的部分。
public class CNetwork{ public List m_Node = new List(); //网络节点集合 public List m_Link = new List(); //网络路段集合 public List m_Origin = new List(); //网络的起点集合 public int m_nNode; //节点数 public int m_nLink; //路段数量 public int m_nOrigin; //起点数量 public double[] LinkFlow; //路段流量 public double[] LinkTravelTime; //路段走行时间 public double MaxUEGap = 1.0e-5; //UE的最大误差 public double UEGap; //UE误差 double[] ShortestPathCost;//临时变量，所有节点到起点的最短路int[] ShortestPathParent;//最短路上，所有节点到起点的在最短路上的前继路段……}(4)构造读取文件的函数
采用如表所示的字段描述交通网络。节点、路段和OD对信息分别储存在三个TXT文档中，每行为一个节点(路段或OD对)，不同字段之间用制表符分隔。三个文档包含的字段信息如下表：
构造三个无返回值的函数读取文件，并将文件中的信息分别读入CNode，CLink和COrigin的代码如下：
// 函数：读取节点文件 public void ReadNode(String DataPath) { if (!File.Exists(DataPath)) // 文件是否存在 { Console.WriteLine("{0} does not exist!", DataPath); return; } StreamReader f1 = File.OpenText(DataPath); string row; string[] Data; m_Node.Clear(); m_nNode = 0; while ((row = f1.ReadLine()) != null) { if (row == "") continue; Data = row.Split('\t'); CNode pNode = new CNode(); pNode.ID = m_nNode; pNode.PositionX = Convert.ToDouble(Data[1]); pNode.PositionY = Convert.ToDouble(Data[2]); m_nNode++; m_Node.Add(pNode); } f1.Close(); } // 函数：读取路段文件 public void ReadLink(string DataPath) { if (!File.Exists(DataPath)) // 文件是否存在 { Console.WriteLine("{0} does not exist!", DataPath); return; } StreamReader f2 = File.OpenText(DataPath); string row; string[] Data; m_Link.Clear(); m_nLink = 0; CNode pNode = new CNode(); while ((row = f2.ReadLine()) != null) { if (row == "") continue; Data = row.Split('\t'); CLink pLink = new CLink(); pLink.ID = m_nLink; pLink.pInNode = (CNode)m_Node[Convert.ToInt32(Data[0]) - 1]; pLink.pOutNode = (CNode)m_Node[Convert.ToInt32(Data[1]) - 1]; pLink.FreeFlowTravelTime = Convert.ToDouble(Data[2]); pLink.Capacity = Convert.ToDouble(Data[3]); pLink.pInNode.OutgoingLink.Add(pLink.ID); pLink.pOutNode.IncomingLink.Add(pLink.ID); m_nLink++; m_Link.Add(pLink); } f2.Close(); } // 函数：读取OD对文件 public void ReadODpairs(string DataPath) { if (!File.Exists(DataPath)) // 文件是否存在 { Console.WriteLine("{0} does not exist!", DataPath); return; } StreamReader f3 = File.OpenText(DataPath); string row; string[] Data; m_Origin.Clear(); CNode pNode; COrigin pOrigin; m_nOrigin = 0; while ((row = f3.ReadLine()) != null) { if (row == "") continue; Data = row.Split('\t'); pNode = (CNode)m_Node[Convert.ToInt32(Data[0]) - 1]; if (pNode.Origin_ID == -1) { pOrigin = new COrigin(); m_nOrigin++; pOrigin.ID = m_Origin.Count; pOrigin.pOriginNode = pNode; pNode.Origin_ID = pOrigin.ID; m_Origin.Add(pOrigin); } else pOrigin = m_Origin[pNode.Origin_ID]; pNode = (CNode)m_Node[Convert.ToInt32(Data[1]) - 1]; pOrigin.DestinationNode.Add(pNode.ID); pOrigin.ODDemand.Add(Convert.ToDouble(Data[2])); } f3.Close(); }(5)构造计算最短路的函数
在Frank-Wolfe算法的编写过程中遇到两种计算最短路的情况。首先是在算法1计算全有全无交通分配的line4中，需要返回最短路上的所有路段集合。其次是在间隙函数式(3)中，计算当前路段流量下的最小出行费用需要用到最短路。
分别构造两个名为GetShortestPath()和GetShortestCost()的函数用于返回最短路上的路段集合和最小出行费用。这里采用的是单队列Label-correcting的最短路算法。
// 计算最短路：返回Start到End最短路径上的 *路段集合* public List<int> GetShortestPath(int Start, int End) { CNode pNode; CLink pLink; // 初始化 int startposition = 0, endposition = 1; ShortestPathCost = new double[m_nNode]; //当前最短路 ShortestPathParent = new int[m_nNode]; //前驱节点 int[] checkList = new int[m_nNode];//队列，循环使用 bool[] binCheckList = new bool[m_nNode]; //是否在队列中 bool[] bscanStatus = new bool[m_nNode]; for (int node = 0; node < m_nNode; node++) { ShortestPathParent[node] = -1; ShortestPathCost[node] = double.MaxValue; binCheckList[node] = false; } ShortestPathCost[Start] = 0; checkList[0] = Start; //开始while循环 while (startposition != endposition) { if (startposition >= m_nNode) startposition = 0; int i = checkList[startposition]; startposition++; pNode = m_Node[i]; for (int index = 0; index < pNode.OutgoingLink.Count; index++) { pLink = m_Link[pNode.OutgoingLink[index]]; int j = pLink.pOutNode.ID; double value = pLink.TravelTime; if (ShortestPathCost[j] > ShortestPathCost[i] + value) { ShortestPathCost[j] = ShortestPathCost[i] + value; ShortestPathParent[j] = i; // 添加到队列尾部 if (endposition >= m_nNode) endposition = 0; checkList[endposition] = j; endposition++; bscanStatus[j] = true; } } } // 返回最短路径集合 List<int> ShortestPath_Link = new List<int>(); int point_out = End; for (; ; ) { int i = 0; int point_in = ShortestPathParent[point_out]; for (int link = 0; link < m_nLink; link++) { pLink = m_Link[link]; if ((point_in == pLink.pInNode.ID) && (point_out == pLink.pOutNode.ID)) { ShortestPath_Link.Insert(0, pLink.ID); point_out = point_in; } } i++; if (point_in == Start) break; } return ShortestPath_Link; } // 计算最短距离：返回Start到End最短路径上的 *最短距离* public double GetShortestCost(int Start, int End) { CNode pNode; CLink pLink; // 初始化 int startposition = 0, endposition = 1; ShortestPathCost = new double[m_nNode]; //当前最短路 ShortestPathParent = new int[m_nNode]; //前驱节点 int[] checkList = new int[m_nNode];//队列，循环使用 bool[] binCheckList = new bool[m_nNode]; //是否在队列中 bool[] bscanStatus = new bool[m_nNode]; for (int node = 0; node < m_nNode; node++) { ShortestPathParent[node] = -1; ShortestPathCost[node] = double.MaxValue; binCheckList[node] = false; } ShortestPathCost[Start] = 0; checkList[0] = Start; //开始while循环 while (startposition != endposition) { if (startposition >= m_nNode) startposition = 0; int i = checkList[startposition]; startposition++; pNode = m_Node[i]; for (int index = 0; index < pNode.OutgoingLink.Count; index++) { pLink = m_Link[pNode.OutgoingLink[index]]; int j = pLink.pOutNode.ID; double value = pLink.TravelTime; if (ShortestPathCost[j] > ShortestPathCost[i] + value) { ShortestPathCost[j] = ShortestPathCost[i] + value; ShortestPathParent[j] = i; // 添加到队列尾部 if (endposition >= m_nNode) endposition = 0; checkList[endposition] = j; endposition++; bscanStatus[j] = true; } } } // 返回最短路距离 return ShortestPathCost[End]; }(6)构造全有全无交通分配的函数
全有全无交通分配直接按照算法1编写。
// 全有全无交通分配：输入为初始路段流量(索引对应CLink中的ID) public double[] AllorNothingAssignment(double[] LinkFlow) // int[] LinkFlow是现有的路段流量 { double[] ANLinkFlow = new double[m_Link.Count]; // 全有全无路段流量 CLink pLink; COrigin pOrinin; for (int link = 0; link < m_nLink; link++) { pLink = m_Link[link]; // 计算路段走行时间 pLink.TravelTime = pLink.FreeFlowTravelTime * (1 + pLink.Alpha * Math.Pow(LinkFlow[link] / pLink.Capacity, pLink.Power)); // 初始化全有全无路段流量 ANLinkFlow[link] = 0; } for (int origin = 0; origin < m_nOrigin; origin++) // 遍历O的循环 { pOrinin = m_Origin[origin]; int OriginNode = pOrinin.pOriginNode.ID; // 任意一个O的编号 for (int i = 0; i < pOrinin.DestinationNode.Count;i++) // 遍历D的循环 { int DestinationNode = pOrinin.DestinationNode[i]; // 对应D的编号 double Demand = pOrinin.ODDemand[i]; // 计算最短路上的路径集合 List<int> ShortestPath_Link = new List<int>(GetShortestPath(OriginNode, DestinationNode)); foreach (int link_index in ShortestPath_Link) { ANLinkFlow[link_index] += Demand; } } } // 输出全有全无结果 // foreach (int index in ANLinkFlow) Console.WriteLine(index); return ANLinkFlow; }(7)构造Frank-Wolfe算法的函数
Frank-Wolfe算法直接按照算法2编写。需要注意的是，在这个函数中引用了间隙函数和计算最优步长时式(14)等号左边的两个辅助函数，将在下一部分给出详细代码。计算最优步长时采用二分法求解式(14)在[0,1]区间内的近似根。
// Frank-Wolfe算法 public void FrankWolfeAlgorithm() { // 初始化： int k = 0; // 迭代次数 LinkFlow = new double[m_nLink]; for (int link = 0; link < m_nLink; link++) LinkFlow[link] = 0; // 所有路段流量为0 LinkFlow = AllorNothingAssignment(LinkFlow); // 做全有全无分配，得到初始路段流量 UEGap = GetUEGap(LinkFlow); Console.WriteLine("第{0}次迭代", k); Console.WriteLine("Gap函数值：" + UEGap); Console.Write("LinkFlow："); foreach (int index in LinkFlow) Console.Write(index + ","); Console.WriteLine(); Console.WriteLine(); // 开始while循环 using (TextWriter tw = File.CreateText("C:/Users/王婷/Desktop/Gap.txt")) { while (UEGap > MaxUEGap) { double[] newLinkFlow = AllorNothingAssignment(LinkFlow); // 计算可行下降方向 double[] FeasibleDescentDirection = new double[m_nLink]; for (int link = 0; link < m_nLink; link++) { FeasibleDescentDirection[link] = newLinkFlow[link] - LinkFlow[link]; } // 计算最优步长(二分法) double paraLamuda; // 最优步长 double left = 0; double right = 1; double mid = 0; double f_left = OptimalStepFunction(FeasibleDescentDirection, LinkFlow, left); // left的求根函数值 double f_right = OptimalStepFunction(FeasibleDescentDirection, LinkFlow, right); // right的求根函数值 double f_mid; if ((f_left * f_right) > 0) // f_left和f_right同号，说明根不在[0,1]区间 { if (Math.Abs(f_left) > Math.Abs(f_right)) paraLamuda = right; else paraLamuda = left; } else // f_left和f_right异号，说明[0,1]区间存在方程的根 { // 用二分法求根 while (right - left > MaxUEGap) { mid = (left + right) / 2; f_left = OptimalStepFunction(FeasibleDescentDirection, LinkFlow, left); f_right = OptimalStepFunction(FeasibleDescentDirection, LinkFlow, right); f_mid = OptimalStepFunction(FeasibleDescentDirection, LinkFlow, mid); if ((f_left * f_mid) > 0) left = mid; else right = mid; } paraLamuda = (left + right) / 2; } // 更新路段流量 for (int link = 0; link < m_nLink; link++) { LinkFlow[link] = LinkFlow[link] + paraLamuda * FeasibleDescentDirection[link]; } k += 1; UEGap = GetUEGap(LinkFlow); // 计算Gap函数 // 将Gap函数值写入文档 tw.WriteLine(UEGap); Console.WriteLine(UEGap); // 输出结果 /* Console.WriteLine("第{0}次迭代", k); Console.WriteLine("Gap函数值：" + UEGap); Console.Write("LinkFlow："); foreach (int index in LinkFlow) Console.Write(index + ","); Console.WriteLine(); Console.WriteLine(); */ } } Console.WriteLine("UE平衡的路段流量："); foreach (int index in LinkFlow) Console.Write(index + ","); }(8)构造其他辅助函数
辅助函数包括间隙函数和计算最优步长等式的左边，具体公式分别见式(3)和(14)。
// 间隙函数 public double GetUEGap(double[] LinkFlow) //获取UE的误差 { COrigin pOrigin; CLink pLink; // 计算被减项分母 + 更新路段走行时间(出行费用：TravelTime) double num1 = 0; for (int link = 0; link < m_nLink; link++) { pLink = m_Link[link]; // 更新TravelTime pLink.TravelTime = pLink.FreeFlowTravelTime * (1 + pLink.Alpha * Math.Pow(LinkFlow[link] / pLink.Capacity, pLink.Power)); // 单一求和项的计算 num1 += (pLink.TravelTime * LinkFlow[link]); } // 计算被减项分子 double num2 = 0; for (int origin = 0; origin < m_nOrigin; origin++) // 遍历O的循环 { pOrigin = m_Origin[origin]; int OriginNode = pOrigin.pOriginNode.ID; // 任意一个O的编号 for (int i = 0; i < pOrigin.DestinationNode.Count; i++) // 遍历D的循环 { int DestinationNode = pOrigin.DestinationNode[i]; // 对应D的编号 // 单一求和项的计算 double Demand = pOrigin.ODDemand[i]; double ODCost = GetShortestCost(OriginNode, DestinationNode); num2 += (Demand * ODCost); } } // 计算间隙函数 UEGap = 1 - num2 / num1; return UEGap; } // 计算最优步长的求根函数的左边 public double OptimalStepFunction(double[] FeasibleDescent, double[] LinkFlow, double paraLamuda) { CLink pLink; double Value = 0; for (int link = 0; link < m_nLink; link++) { pLink = m_Link[link]; double para = LinkFlow[link] + paraLamuda * FeasibleDescent[link]; Value += FeasibleDescent[link] * (pLink.FreeFlowTravelTime * (1 + pLink.Alpha * Math.Pow(para / pLink.Capacity, pLink.Power))); } return Value; }5.2 代码文件Program的编写
在前述Network文件中编写的是Frank-Wolfe算法的计算过程，要使Network文件中的代码块能被调用且顺利运行，还需要设置程序的入口，即main()函数。Program文件中main()函数的操作包括用ReadNode()，ReadLink()，ReadODPairs()读取具体文件，和运行Frank-Wolfe算法。
using System;using System.Collections.Generic;using System.Linq;using System.Text;using System.Threading.Tasks;namespace WT_TrafficAssignment{ class Program { static void Main(string[] args) { CNetwork Network = new CNetwork(); // 读取文件 Network.ReadNode("C:\\……\\Node.txt"); Network.ReadLink("C:\\……\\Link.txt"); Network.ReadODpairs("C:\\……\\ODPairs.txt"); // 运行FW算法 Network.FrankWolfeAlgorithm(); Console.ReadKey(); } }}6 案例与结果
点击阅读原文获取案例中的数据！！！
6.1 样例网络
先在一个简单的样例网络上检验Frank-Wolf算法程序的计算结果，样例网络如下图所示。
(1)计算结果
代码输出的迭代情况、Gap函数值、及UE平衡的路段流量如下。可见模型经过一次迭代，Gap函数的精度满足要求而停止计算。
UE平衡时的路段流量见下图。
(2)正确性检验
OD对(2，3)：只有唯一的路径，不做检验；
OD对(1，4)：有3条可行路径，分别为：1-2-4，1-3-4，1-2-3-4，计算平衡路段流量条件下的走行时间：
被采用的路径1和路径2平衡路段流量条件下的走行时间近似相等，而未被采用的路径3走行时间T3明显大于T1和T2，满足UE平衡条件。
6.2 Sioux Falls网络
算法在比较复杂的Sioux Falls网络中经过10044次迭代达到UE平衡，平衡时的路段流量如下图所示(流量大小的数字标签更接近O点)：
用Frank-wolf算法求解UE交通分配问题的收敛过程如下。
参考：
《城市动态交通流分配模型与算法》高自友，任华玲著。
《交通网络优化计算机编程讲义》
《最优化与最优控制》
阅读原文链接：
提取码：qy7s
编辑：庄桢
“交通科研Lab”：分享学习点滴，期待科研交流！
如果觉得还不错
点这里！???
716
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
