1 - ⚡ Zap for Network Utility Maximization Basics ⚡

In this notebook, we introduce how to construct and solve NUM problems in Zap. We will manually construct a small network with links and routes, solve the NUM problem using both CVXPY and PMP (ADMM), and analyze the results.

1. Creating a NUM Problem

2. Solving NUM

3. Analyzing Results

import zap
import torch

import numpy as np
import cvxpy as cp
from scipy.sparse import csc_matrix


from experiments.resource_opt_solve.benchmarks.nu_opt_benchmark import NUOptBenchmarkSet
from zap.resource_opt.nu_opt_bridge import NUOptBridge
from zap.admm import ADMMSolver

Creating NUM Problems

NUM problems in Zap consist of a link-route matrix, which defines the underlying links and streams, link capacities, stream weights, and a list of which traffic streams have linear utilities (Zap assumes by default that all streams have logarithmic utilities).

We provide a synthetic data generator which creates NUM test instances. We can, for example, generate a problem with the following parameters:

• 2000 links

• 1000 streams

• 10 links in the average stream

• Capacities uniformly in the range (0.1, 1)

• 0% of links congested

• Uniform stream weights of 1

• Logarithmic utilities for all streams

base_seed = 42
m = 2000 # Number of links
n = 1000 # Number of streams
avg_stream_length = 3 # Average stream uses 10 links
capacity_range = (0.1, 1)
link_congest_num_frac = 0
lin_util_frac = 0

benchmark = NUOptBenchmarkSet(num_problems=1,
        m=m,
        n=n,
        avg_route_length=avg_stream_length,
        capacity_range=capacity_range,
        link_congest_num_frac=link_congest_num_frac,
        base_seed=base_seed)

We can extract the desired problem parameters from the synthetic data benchmark generator.

R, capacities, w, linear_flow_idxs = benchmark.get_data(0)
nu_opt_params = {
    "R": R,
    "capacities": capacities,
    "w": w,
    "lin_device_idxs": linear_flow_idxs,
}

While we can use the synthetic data generator to create test instances, for demonstration purposes, we instead solve an instance of ‘A Small Example’ from “Large-Scale Network Utility Maximization via GPU-Accelerated Proximal Message Passing” (i.e. we use the same link-route matrix from Fig. 1).

We generate a very small problem with the following parameters:

• 3 links

• 4 streams

• Capacities uniformly in the range (0.1, 1)

• 0% of links congested

• Uniform stream weights of 1

• Logarithmic utilities for all streams

m = 3
n = 4
rng = np.random.default_rng(42)
capacities = rng.uniform(0.1, 1, size=m) # Capacity range (0.1, 1)
w = np.ones(n)
lin_device_idxs = None # No devices with linear utility components


# Create link-route matrix R as CSC matrix
data = np.array([1, 1, 1, 1, 1])
row  = np.array([0, 1, 1, 2, 2])
col  = np.array([0, 1, 2, 1, 3])
R = csc_matrix((data, (row, col)), shape=(m, n))

nu_opt_params = {
    "R": R,
    "capacities": capacities,
    "w": w,
    "lin_device_idxs": lin_device_idxs,
}

Before we construct the NUM problem to solve, we need to finally specify how we group streams. We create a separate group for each set of streams with the same number of terminals (i.e. streams that utilize the same number of links).

grouping_params = {"variable_grouping_strategy": "discrete_terminal_groups"}

We now create the NUM problem. NUOptBridge provides the interface for taking problem data and converting it into a NUM problem that we can solve. NuOptBridge converts the problem data into its corresponding bipartite graph representation which is used by the solver.

nu_opt_bridge = NUOptBridge(nu_opt_params, grouping_params)

Solving a NUM Problem

We can interface with the nu_opt_bridge object to solve our NUM problem in one of two ways:

1. cvxpy - This method will build a model in cvxpy and send it to an off-the-shelf solver, such as Clarabel. You can also use commerical solvers like MOSEK. This approach is best for small to medium problems and finds highly accurate solutions.

2. admm - This method will use a PyTorch implementation of the alternating direction method of multipliers (ADMM). This can be run on either a CPU or GPU and scales well to larger problems. In general, this method is only capable of finding medium accuracy solutions within a reasonable amount of time.

Solving with CVXPY

To solve with cvxpy, we simply call nu_opt_bridge.solve and (optionally) specify a solver.

outcome = nu_opt_bridge.solve(solver=cp.CLARABEL)

We will shortly address how to interpret the solver output.

If you are using the benchmark set, you can also extract the benchmark problem via CVXPY and solve it directly as well.

for i, prob in enumerate(benchmark):
        problem = prob

problem.solve(solver=cp.CLARABEL, verbose=False)

-1909.3491417000203

Solving with ADMM (Proximal Message Passing)

Solving with ADMM is a little more complicated and requires two steps:

1. Transfering device data to PyTorch

2. Initializing an ADMM solver object.

machine = "cpu" # Pick "cuda" for a machine with an NVIDIA GPU
dtype = torch.float32
admm_devices = [d.torchify(machine=machine, dtype=dtype) for d in nu_opt_bridge.devices]
admm = ADMMSolver(
    machine=machine,
    dtype=dtype,
    atol=1e-6,
    rtol=1e-6,
    tau=2,
    alpha=1.6,
    num_iterations=1000
)
solution_admm, history_admm = admm.solve(nu_opt_bridge.net, admm_devices, nu_opt_bridge.time_horizon)

ADMM converged in 104 iterations.

# ADMM solutions need to be cast to a standard DispatchOutcome
outcome_admm = solution_admm.as_outcome()

outcome_admm

DispatchOutcome(phase_duals=[None, None, None], local_equality_duals=None, local_inequality_duals=None, local_variables=[[tensor([[0.7966],
        [0.2918],
        [0.6695]])], [tensor([[0.2032]])], None], power=[[tensor([[0.7966],
        [0.2918],
        [0.6695]])], [tensor([[0.2032]]), tensor([[0.2032]])], [tensor([[-0.7966],
        [-0.4950],
        [-0.8727]])]], angle=[None, None, None], prices=tensor([[-1.2554],
        [-3.4273],
        [-1.4936]]), global_angle=None, problem=None, ground=None)

Analyzing Results

Results are packaged into a hierachically structured DispatchOutcome object. This allows us to directly access the terminal flows in the bipartite representation of the problem.

1. At the top level, a DispatchOutcome has fields such as power and prices. You can access these fields like any other Python field, e.g., outcome.power.

2. Each field contains either stream-specific information or global information. Stream-specific fields will contain a list of length len(streams), where the ith entry in the list contains information specific to the ith stream. Global fields contain a 2d array of size (num_nodes, 1). You can access the information for stream i using normal indexing, e.g., outcome.power[i].

3. For stream-specific information, each block of information is further broken down by the terminal of the stream. Slack streams have just a single terminal and will always be indexed as outcome.power[i][0]. Other utility streams may have two or more terminals. The data for terminal j is stored in outcome.power[i][j].

4. Finally, the data for terminal j of stream i is just a 2d array of size (num_devices, 1), where num_devices = devices[i].num_devices is the number of devices in device group i.

outcome.prices  # Global, link duals

array([[-1.25535188],
       [-3.42724732],
       [-1.4936407 ]])

outcome.power # Stream-specific (utility and slack streams)

[[array([[0.79656044],
         [0.29177735],
         [0.66952488]])],
 [array([[0.20321324]]), array([[0.20321324]])],
 [array([[-0.79656044],
         [-0.4949906 ],
         [-0.87273812]])]]

outcome.power[0] # Terminal flows for the three single terminal utility streams (Traffic Streams 1, 3, 4) (x1, x3, x4)

[array([[0.79656044],
        [0.29177735],
        [0.66952488]])]

outcome.power[1] # Terminal flows for the one two terminal utility stream (Traffic Stream 2) (x2)

[array([[0.20321324]]), array([[0.20321324]])]

outcome.power[2] # Terminal flows for the three single terminal slack streams (Slack Streams 1, 2, 3) (s1-c1, s2-c2, s3-c3)

[array([[-0.79656044],
        [-0.4949906 ],
        [-0.87273812]])]