Skip to main content
`Ctrl`+`K`

# Traffic Assignment

Similar to other complex algorithms that handle a large amount of data through complex computations, traffic assignment procedures can always be subject to at least one very reasonable question: Are the results right?

For this reason, we have used all equilibrium traffic assignment algorithms available in AequilibraE to solve standard instances used in academia for comparing algorithm results.

Instances can be downloaded here.

All tests were performed with the AequilibraE version 1.1.0.

## Validation

As shown below, the results produced by AequilibraE are within expected, although some differences have been found, particularly for Winnipeg. We suspect that there are issues with the reference results and welcome further investigations.

- Links: 39,018

- Nodes: 12,982

- Zones: 1,790

- Links: 2,522

- Nodes: 1,020

- Zones: 110

- Links: 914

- Nodes: 416

- Zones: 38

- Links: 914

- Nodes: 416

- Zones: 38

- Links: 76

- Nodes: 24

- Zones: 24

## Convergence Study

Besides validating the final results from the algorithms, we have also compared how well they converge for the largest instance we have tested (Chicago Regional), as that instance has a comparable size to real-world models.

Not surprisingly, one can see that Frank-Wolfe far outperforms the Method of Successive Averages for a number of iterations larger than 25 in the case of Chicago, and is capable of reaching 1.0e-04 just after 800 iterations, while MSA is still at 3.5e-4 even after 1,000 iterations for that same case.

The actual show, however, is left for the biconjugate Frank-Wolfe implementation, which delivers a relative gap of under 1.0e-04 in under 200 iterations, and a relative gap of under 1.0e-05 in just over 700 iterations.

This convergence capability, allied to its computational performance described below suggest that AequilibraE is ready to be used in large real-world applications.

## Computational performance

All tests were run on a workstation equipped AMD Threadripper 3970X with 32 cores (64 threads) @ 3.7 GHz (memory use is trivial for these instances).

On this machine, AequilibraE performed 1,000 iterations of biconjugate Frank-Wolfe assignment on the Chicago Network in a little over 4 minutes, or a little less than 0.43s per iteration.

Compared with AequilibraE previous versions, we can notice a reasonable decrease in processing time.

Note

The biggest opportunity for performance in AequilibraE right now it to apply network contraction hierarchies to the building of the graph, but that is still a long-term goal

## Want to run your own convergence study?

If you want to run the convergence study in your machine, with Chicago Regional instance or any other instance presented here, check out the code block below! Please make sure you have already imported TNTP files into your machine.

In the first part of the code, we’ll parse TNTP instances to a format AequilibraE can understand, and then we’ll perform the assignment.

    # Imports
    from pathlib import Path
    from time import perf_counter

    import matplotlib.pyplot as plt
    import numpy as np
    import pandas as pd
    import seaborn as sns
    from sklearn.linear_model import LinearRegression
    from sklearn.metrics import r2_score

    from aequilibrae.matrix import AequilibraeMatrix
    from aequilibrae.paths import Graph
    from aequilibrae.paths import TrafficAssignment
    from aequilibrae.paths.traffic_class import TrafficClass

    # Helper functions
    def build_matrix(folder: Path, model_stub: str) -> AequilibraeMatrix:
        omx_name = folder / f"{model_stub}_trips.omx"
        if omx_name.exists():
            mat = AequilibraeMatrix()
            mat.load(omx_name)
            mat.computational_view()
            return mat

        matfile = str(folder / f"{model_stub}_trips.tntp")
        # Creating the matrix
        f = open(matfile, 'r')
        all_rows = f.read()
        blocks = all_rows.split('Origin')[1:]
        matrix = {}
        for k in range(len(blocks)):
            orig = blocks[k].split('\n')
            dests = orig[1:]
            orig = int(orig[0])

            d = [eval('{' + a.replace(';', ',').replace(' ', '') + '}') for a in dests]
            destinations = {}
            for i in d:
                destinations = {**destinations, **i}
            matrix[orig] = destinations
        zones = max(matrix.keys())
        index = np.arange(zones) + 1
        mat_data = np.zeros((zones, zones))
        for i in range(zones):
            for j in range(zones):
                mat_data[i, j] = matrix[i + 1].get(j + 1, 0)

        # Let's save our matrix in AequilibraE Matrix format
        mat = AequilibraeMatrix()
        mat.create_empty(zones=zones, matrix_names=['matrix'], memory_only=True)
        mat.matrix['matrix'][:, :] = mat_data[:, :]
        mat.index[:] = index[:]
        mat.computational_view(["matrix"])
        mat.export(str(omx_name))
        return mat

    # Now let's parse the network
    def build_graph(folder: Path, model_stub: str, centroids: np.array) -> Graph:
        net = pd.read_csv(folder / f"{model_stub}_net.tntp", skiprows=7, sep='\t')
        cols = ['init_node', 'term_node', 'free_flow_time', 'capacity', "b", "power"]
        if 'toll' in net.columns:
            cols.append('toll')
        network = net[cols]
        network.columns = ['a_node', 'b_node', 'free_flow_time', 'capacity', "b", "power", "toll"]
        network = network.assign(direction=1)
        network["link_id"] = network.index + 1
        network.free_flow_time = network.free_flow_time.astype(np.float64)

        # If you want to create an AequilibraE matrix for computation, then it follows
        g = Graph()
        g.cost = net['free_flow_time'].values
        g.capacity = net['capacity'].values
        g.free_flow_time = net['free_flow_time'].values

        g.network = network
        g.network.loc[(g.network.power < 1), "power"] = 1
        g.network.loc[(g.network.free_flow_time == 0), "free_flow_time"] = 0.01
        g.prepare_graph(centroids)
        g.set_graph("free_flow_time")
        g.set_skimming(["free_flow_time"])
        g.set_blocked_centroid_flows(False)
        return g

    def known_results(folder: Path, model_stub: str) -> pd.DataFrame:
        df = pd.read_csv(folder / f"{model_stub}_flow.tntp", sep='\t')
        df.columns = ["a_node", "b_node", "TNTP Solution", "cost"]
        return df

    # Let's run the assignment
    def assign(g: Graph, mat: AequilibraeMatrix, algorithm: str):
        assigclass = TrafficClass("car", g, mat)
        if "toll" in g.network.columns:
            assigclass.set_fixed_cost("toll")

        assig = TrafficAssignment()
        assig.set_classes([assigclass])
        assig.set_vdf("BPR")
        assig.set_vdf_parameters({"alpha": "b", "beta": "power"})
        assig.set_capacity_field("capacity")
        assig.set_time_field("free_flow_time")
        assig.max_iter = 1000
        assig.rgap_target = 1e-10 # Nearly guarantees that convergence won't be reached
        assig.set_algorithm(algorithm)
        assig.execute()
        return assig

    # We compare the results
    def validate(assig: TrafficAssignment, known_flows: pd.DataFrame, algorithm: str, folder: Path, model_name):
        modeled = g.network[["link_id", "a_node", "b_node"]].merge(assig.results().matrix_ab.reset_index(),
                                                                   on="link_id").rename(
            columns={"matrix_ab": "AequilibraE Solution"})
        merged = known_flows.merge(modeled, on=["a_node", "b_node"])

        # Scatter plot
        plt.figure(figsize=(10, 6))
        sns.scatterplot(data=merged, x="TNTP Solution", y="AequilibraE Solution", s=30)

        # Linear regression
        X = merged["TNTP Solution"].values.reshape(-1, 1)
        y = merged["AequilibraE Solution"].values
        reg = LinearRegression(fit_intercept=False).fit(X, y)
        y_pred = reg.predict(X)
        r_squared = r2_score(y, y_pred)

        # Plot regression line
        plt.plot(merged["TNTP Solution"], y_pred, color='red', label='Linear regression')

        # Customize the plot
        plt.title(f'Comparison of Known and AequilibraE Solutions - Algorithm: {algorithm}', fontsize=16)
        plt.xlabel('Known Solution', fontsize=14)
        plt.ylabel('AequilibraE Solution (1,000 iterations)', fontsize=14)

        # Display the equation and R-squared on the plot
        equation_text = f'y = {reg.coef_[0]:.2f}x\nR-squared = {r_squared:.5f}'
        plt.text(x=merged["TNTP Solution"].max() * 0.5, y=merged["AequilibraE Solution"].max() * 0.85, s=equation_text,
                 fontsize=12)

        plt.legend()
        plt.savefig(folder / f"{model_name}_{algorithm}-1000_iter.png", dpi=300)
        plt.close()

    def assign_and_validate(g: Graph, mat: AequilibraeMatrix, folder: Path, model_stub: str):
        known_flows = known_results(folder, model_stub)
        # We run the traffic assignment
        conv = None
        for algorithm in ["bfw", "cfw", "fw", "msa"]:
            t = -perf_counter()
            assig = assign(g, mat, algorithm)
            t += perf_counter()
            print(f"{model_stub},{algorithm},{t:0.4f}")

            res = assig.report()[["iteration", "rgap"]].rename(columns={"rgap": algorithm})
            validate(assig, known_flows, algorithm, folder, model_stub)

            conv = res if conv is None else conv.merge(res, on="iteration")
        df = conv.replace(np.inf, 1).set_index("iteration")
        convergence_chart(df, data_folder, model_stub)
        df.to_csv(folder / f"{model_stub}_convergence.csv")

    def convergence_chart(df: pd.DataFrame, folder: Path, model_name):
        import matplotlib.pyplot as plt

        plt.cla()
        df = df.loc[15:, :]
        for column in df.columns:
            plt.plot(df.index, df[column], label=column)
        # Customize the plot
        plt.title('Convergence Comparison')
        plt.xlabel('Iterations')
        plt.ylabel('RGap')
        plt.yscale("log")
        plt.legend(title='Columns')
        plt.savefig(folder / f"convergence_comparison_{model_name}.png", dpi=300)

    models = {"chicago": [Path(r'D:\src\TransportationNetworks\chicago-regional'), "ChicagoRegional"],
          "sioux_falls": [Path(r'D:\src\TransportationNetworks\SiouxFalls'), "SiouxFalls"],
            "anaheim": [Path(r'D:\src\TransportationNetworks\Anaheim'), "Anaheim"],
            "winnipeg": [Path(r'D:\src\TransportationNetworks\Winnipeg'), "Winnipeg"],
            "barcelona": [Path(r'D:\src\TransportationNetworks\Barcelona'), "Barcelona"],
          }

    convergence = {}
    for model_name, (data_folder, model_stub) in models.items():
        print(model_name)
        mat = build_matrix(data_folder, model_stub)
        g = build_graph(data_folder, model_stub, mat.index)
        assign_and_validate(g, mat, data_folder, model_stub)