P-step CVRP and VRPTW extension

Documentation index

baldes::extensions::pstep implements the compact continuous relaxation in Proposition 1 of A p-step Formulation for the Capacitated Vehicle Routing Problem, Dollevoet, Munari and Spliet.

The extension is available without Cargo features. Solving requires an implementation of Baldes's LpMaster, such as Clp (--features clp).

cargo build --release --features clp --bin baldes
./target/release/baldes pstep instances/cvrplib/P-n16-k8.vrp --p 2 --threads 4
./target/release/baldes pstep instances/cvrplib/P-n16-k8.vrp --p 3 --threads 4 --enumerate --json pstep.json

Use baldes pstep --help for time, iteration and label limits. --enumerate checks the bound against the complete column pool; use it on small cases. The JSON separates a converged full-model lp_bound from a potentially nonconverged restricted_objective. A work-limit error cannot certify a bound. Neither output is an integer routing solution. The fleet is unlimited even if the CVRPLIB instance name or header specifies a conventional truck count.

Rust callers can use Problem::from_instance on a parsed inputs::Instance. It preserves the input's distance rounding. Then call Problem::solve_column_generation with an empty LpMaster and ColumnGenerationOptions; pricing.p and pricing.threads control the step size and CPU parallelism; pricing.method selects Auto (default), Enumeration, or Bucket. The older examples/pstep.rs illustrates full-pool construction, and examples/pstep_cvrp.rs produces CSV benchmark records.

Create a Problem with customer demands, vehicle capacity and a directed cost matrix. Matrix indices are start depot 0, customers 1..=n, end depot n+1. For a conventional single-depot instance, copy depot-to-customer costs into row 0 and customer-to-depot costs into column n+1. None means no arc. Customer demands must be strictly positive; capacity and supplied costs must be finite, with nonnegative costs. The model has an unlimited homogeneous fleet.

Problem::enumerate builds the complete column pool for a supplied positive p. Every partial path is elementary and capacity feasible. Paths have exactly p arcs, except those starting at the depot, which may be shorter. Empty depot routes are omitted. At large p, the pool contains full routes and depot prefixes; flow conservation forces incomplete prefixes to zero.

The compact prior-demand choices are zero for depot starts, maximum feasible prior demand for depot ends, and both extremes for internal paths. Coincident extremes are stored once. A junction customer's demand belongs to both adjacent paths; the load coefficients account for this overlap.

Formulation::entries exposes sparse columns in this row order:

Formulation::into_master populates an empty LP backend with nonnegative continuous columns. Solve through the backend's status and time-limit API. An optimal solve gives the full p-step LP bound to solver tolerance. At p=1 this has the single-commodity flow bound; at p>=n+1 it has the set-partitioning bound. Bounds need not increase at every successive p; the paper guarantees nondecrease when p is replaced by an integer multiple.

The reference enumerate implementation has combinatorial cost. Both column count and search work are bounded by EnumerationOptions. Hitting a limit returns an error instead of returning a restricted pool as a complete model. It is intended as an exact column-pool reference for small instances and small p.

price_steps solves independent jobs for each start, end and prior-demand extreme. The default Auto method uses streaming elementary enumeration when p<=6 and the dense prefix estimate 1+n+...+n^(p-1) is at most 100,000; otherwise it uses the core bucket-graph labeling engine. This is a bounded size heuristic, not a prediction that one backend always wins. Capacity and missing arcs prune enumeration further. Enumeration retains just the best path per job, so it does not build the full column pool. Both methods solve the same exact pricing problem and enforce work/time limits.

Use --pricing enumeration or --pricing bucket to compare a specific case; --enumerate remains the independent full-pool LP cross-check. --labels limits generated labels for bucket pricing and attempted extensions for streaming enumeration. JSON reports these separately as generated_labels and enumerated_extensions, and records the requested and selected methods. Exceeding either budget aborts pricing without certifying convergence.

For bucket pricing, hop count is an equality resource, preventing artificial waiting from creating extra arcs. Capacity is disposable: every load window has lower bound zero, so waiting never helps and a cheaper label using less capacity safely dominates. Reduced costs use actual customer demands, not any resource waiting. Customer starts require exactly p arcs; depot starts allow 1..p. The fixed endpoints and extreme make the reduced cost additive along arcs, including degree, endpoint and load duals. Each returned path is replayed against the master coefficients.

Workers share immutable instance costs and duals. Each owns its search storage and results. Results are merged in fixed job order and all threads join before the main thread updates the LP. threads=1 runs the same pricing algorithm, making serial/parallel results reproducible. Parallel overhead may outweigh the benefit on very small cases.

The solver starts with degree artificials, uses exact pricing to complete phase I, and then switches to travel costs. Convergence requires every pricing job to finish without a negative reduced-cost step at the requested tolerance. Steps are deduplicated by node order and prior demand. LP columns are added in batches. Restricted-master objectives are not full-model lower bounds until pricing converges. Reported bounds are numerical, not independently certified.

For p<=3 the path topology itself enforces elementarity. For p=4, a relaxed bucket solve is sufficient when its minimum is elementary. Otherwise, pricing partitions by the first arc and forbids later returns to its head. Each resulting four-arc path is elementary; every feasible elementary step belongs to one partition. The relaxed minimum also bounds these partitions. Longer steps use exact cycle exclusion around bucket solves: a cyclic relaxed minimum branches on forbidding each arc of a detected cycle. Every elementary path lies in at least one child; relaxed costs bound each child. This avoids the core's 15-element local NG mask limit. It can still require exponential work for larger p. PricingOptions bounds total labels and bucket subproblems per endpoint job, as well as worker count and time. Exceeding a pricing limit returns an error.

Time windows, cuts, integer solving and integration with Baldes's branch-and-price tree are not implemented here. For integer optimization of this compact model, impose binary aggregate arc flows, as in equations (7)–(9); making extreme p-step columns binary is not equivalent. The exposed paths and coefficients provide the data needed for a separate arc-flow integrality implementation.

Run node scripts/benchmarks/benchmark_pstep.mjs ./target/release/baldes from the repository root to repeat the known-instance validation and serial/parallel timing checks. It writes raw JSON and CSV records under results/pstep/ (override with PSTEP_OUT) and compares every bound with full column enumeration.

Add --grid to run p=2,3,4 on all three bundled CVRPLIB instances, with four workers, 120 seconds per solve and five million labels per endpoint job. Grid records go under results/pstep-20261020/p234/. Reference enumeration is omitted for the larger p=4 pools and A-n32-k5 at p=3, whose pool exceeds the CLI's column limit. Pricing stops are recorded explicitly, without a claimed LP bound.

The validation runs gave these numerical bounds (four workers, development profile). The normal CVRP column uses default root NG-route pricing; the integer values were obtained separately with the normal tree solver.

The A-n32-k5 p=4 run stopped at 120 seconds without a full LP bound. The other eight cases converged; six also matched complete-column enumeration within 1e-5. The A-n32-k5 integer solve used capacity cuts, rank-one cuts and a one-second HGS warm start before closing the gap at 784.

Compare backends with node scripts/benchmarks/benchmark_pstep_pricing.mjs. It alternates backend order over three repetitions, checks converged bounds agree, and writes raw JSON and CSV to results/pstep-pricing/ (or the second positional argument). Set PSTEP_RUNS, PSTEP_SECONDS, or PSTEP_THREADS to change repetitions, per-solve limits, or worker count. An optional binary path and output directory are accepted.

The 2026-10-21 release comparison used an AMD Ryzen 9 9950X, four workers, three repetitions with alternating backend order, no warmup runs, and five million labels/extensions per job. The binary was built with cargo build --release --features clp --bin baldes (opt-level 3, thin LTO). These are end-to-end solve medians, including the LP; they are local timings, not a guarantee for larger or different instances.

The initial run used a 15-second limit; the extended run allowed 20 seconds. E-n22 p=4 needed the extended limit. Converged backend bounds agreed within 1e-5. Enumeration also converged for A-n32 p=4 at 691.463021739; this is a numerical p-step LP bound. The larger p=4 full-pool model was not independently solved for that instance. Raw records, binary/input hashes, and compiler metadata are under results/pstep-pricing-20261021/final/ and verified/.

Experimental CUDA pricing

Build with an NVIDIA CUDA toolkit (nvcc on PATH) and a working CUDA driver:

cargo build --release --features clp,cuda --bin baldes
./target/release/baldes pstep instances/cvrplib/P-n16-k8.vrp --p 3 --gpu --rounding none
cargo test --release --features cuda --lib extensions::pstep
cargo run --release --features cuda --example pstep_gpu_bench
python3 scripts/benchmarks/benchmark_pstep_cuda.py

--gpu is the CLI shorthand for --pricing cuda.

PricingMethod::Cuda runs exact elementary enumeration for each endpoint/load extreme on an NVIDIA GPU. It supports effective p <= 6, uses double precision, and preserves CPU enumeration traversal and tie order. It is opt-in; Auto keeps selecting CPU methods. The LP and reduced-cost replay remain on CPU. threads applies to CPU methods; CUDA reports zero CPU pricing workers.

The prototype allocates and transfers data each round. Timing includes these costs, synchronization, and result replay. Runtime errors and extension limits abort the entire round. The deadline is checked after synchronization: an active kernel is not preempted, and an over-budget result is discarded. Large p or large dense instances can therefore overrun a requested wall-clock limit.

The synthetic benchmark warms every configuration once, then alternates five measured rounds on identical problems and duals, comparing every returned path. The CVRPLIB script uses three fresh-process runs per configuration, alternates CPU/CUDA order, fixes rounding to none, and checks converged LP bounds.

Time windows

The CLI automatically enforces finite service-start windows in Solomon inputs:

cargo build --release --features clp,cuda --bin baldes
./target/release/baldes pstep instances/Solomon/C201.txt --p 3 --gpu --rounding none --json c201-gpu.json
./target/release/baldes pstep instances/Solomon/C201.txt --p 3 --pricing enumeration --threads 16 --rounding none --json c201-cpu.json
cargo test --features cuda,clp --test pstep_vrptw
cargo run --release --features cuda --example pstep_vrptw_bench -- instances/Solomon/C201.txt instances/Solomon/C202.txt

Rust callers use TimeWindowProblem::from_instance. TimeWindowProblem wraps the existing CVRP data and adds windows and service durations for both depot copies and every customer. Travel duration equals travel cost, in the same units as the windows. Finite ordered windows and nonnegative service times are required. Auto uses CPU enumeration; Cuda uses the GPU; bucket pricing is currently unsupported for this model. CUDA still requires effective p <= 6.

The master has four blocks of n rows: degree, endpoint flow, cumulative load, and service-start flow. A timed step contributes its start time positively at its first customer and its end time negatively at its last customer. Service starts of outgoing steps must therefore cover incoming end times. This follows the endpoint-time extension in Vaes (2022), sections 3.2–3.3.

For a fixed path, propagation represents its earliest end service time as max(start + transit, barrier). transit includes service on each arc's tail; barrier accounts for waiting. Latest feasible starts are bounded by every visited due time. Pricing considers the nondominated earliest scheduling breakpoint and the latest feasible start, together with the existing load extremes. This retains compatible time assignments in the master instead of only filtering individually feasible segments.

The result remains an unlimited-fleet LP relaxation, with no integer route certificate. The JSON identifies problem: VRPTW and time_windows_enforced. A converged LP bound is distinct from a restricted-master objective. GPU/CPU comparisons must use identical rounding; the examples explicitly use none.

Measured Solomon comparison

On a Ryzen 9 9950X and RTX 5090, three fresh-process runs with all 100 customers and --rounding none gave these median LP solve times. P-step CPU pricing used 16 workers; the normal solver used its default settings.

All CPU/GPU p-step pairs converged to matching bounds. The default solver obtained 591.55655667 for both cases. C201 p-step obtained the same bound; C202 p-step obtained 590.78102398 at p=2 and 591.54926228 at p=3, so its faster times do not represent equally strong bounds. GPU startup outweighed pricing savings at p=2. At p=3, warmed, matched-dual pricing rounds were 6.96x (C201) and 7.25x (C202) faster on GPU; whole-solve improvements were 1.27x and 1.33x because LP solving remained on CPU.

Machine details, input/binary hashes, timing ranges and LP bounds record the measurements. These timings precede the --gpu alias and integration of the unrelated imported-pricing update; the measured CUDA and VRPTW algorithms are identical. The normal solver comparison uses its stronger route LP.