We are working with the NVIDIA cuOpt team to evaluate GPU-accelerated linear programming for stochastic capacity expansion. The largest instance in this experiment has 100 scenarios, 185 million variables, 140 million constraints, and 471 million nonzeros. We tested whether cuOpt’s first-order PDLP method could solve the full model as a single linear program, and how solve time changed from one to eight GPUs. Faster solves at this scale translate directly to more scenarios evaluated per planning cycle, earlier visibility into capacity risks, better capital allocation decisions under uncertainty, and the ability to respond to demand shifts or supply disruptions before they become unmanageable.
Why linear programming was late to GPUs
GPUs transformed scientific computing and machine learning long before they became practical for large linear programs. Simplex and interior-point methods rely on sparse matrix factorization, an operation that is hard to distribute across thousands of parallel cores. For years, the honest answer to “can I use my GPU for solving LPs?” was not really.
That answer has changed for two reasons. First, GPU sparse direct solvers such as NVIDIA cuDSS allow interior-point methods to perform factorization on the GPU: MadIPM, CuClarabel, and cuOpt’s barrier method follow this general path. Second, first-order methods avoid factorization altogether. Our experiment takes the second path, using a method called PDLP, a primal-dual hybrid gradient (PDHG) method for LP.
PDLP solves the saddle-point formulation of a linear program without matrix factorization. Instead, it repeatedly performs sparse matrix-vector products, an operation well suited to GPU parallelism. We tested NVIDIA cuOpt on configurations of one, two, four, and eight NVIDIA B200 GPUs.
Julia is a common thread in this work. Both the original PDLP prototype and cuPDLP.jl were written in Julia, as were the GPU interior-point solvers MadIPM and CuClarabel. At PSR, our experience with first-order methods includes ProxSDP.jl, an open-source semidefinite programming solver based on the same Chambolle–Pock primal–dual hybrid gradient scheme. It addresses a different problem class, but shares the idea of using inexpensive iterations without matrix factorization.
Simplex has been refined since 1947, and interior-point methods have benefited from four decades of development. GPU implementations of PDLP are much newer. That difference in maturity leaves room for further improvements, and makes the next few years particularly interesting.
Problem overview
Capacity expansion planning asks a simple question with a difficult answer: given demand growth, technology costs, and the existing power system, what should we build, where, and when?
The harder modeling question is uncertainty. Hydrology, wind and solar availability, demand, and fuel prices all shape the value of an investment. A plan that looks optimal in an average year can perform poorly when exposed to a dry year or a low-wind week. A plan is only as good as the range of conditions it was chosen against.
We therefore solve a two-stage stochastic program:
- First stage: Decide how much of each project to build and where.
- Second stage: Operate that fixed plan across many realistic scenarios and measure its cost.
The optimizer minimizes investment cost plus the expected operating cost across the complete scenario set. It decides investment and operation jointly in one model.
We model candidate assets as continuous rather than discrete expansion decisions. This reflects the modular nature of renewable power plants and battery systems, where developers add panels, turbines, and racks incrementally. It also keeps the problem linear instead of turning it into a mixed-integer program.
Writing every scenario out explicitly turns that model into a single large linear program. This formulation is called the deterministic equivalent, and it is what we hand to the solver.
Three modeling requirements drive the size of this formulation:
- Hourly resolution. Renewables and storage operate on intraday timescales. A battery’s value comes from shifting energy between the solar peak and the evening ramp. Aggregate a day into load blocks, and you erase the very behavior you are trying to size. So, the operational subproblem preserves chronological hourly detail and links storage states between consecutive hours.
- Full-year chronology. Long-duration storage investments such as pumped-storage hydropower, as well as conventional reservoirs, can shift energy or water across days and seasons. To represent that behavior, the operational subproblem keeps all 8,760 hours in sequence, with storage states linked throughout the year.
- Many scenarios. A handful of scenarios cannot represent the dry years or low-wind weeks that drive reliability decisions. Our target formulation uses 100 scenarios to represent this uncertainty.
Multiply that through: 8,760 hours 100 scenarios = 876,000 operational time steps, each carrying its own dispatch decisions.
The deterministic equivalent grows linearly with the number of scenarios, reaching 185 million variables and 140 million constraints at 100 scenarios. Factorization-based methods can require substantially more storage than the original constraint matrix because of fill-in. This is one reason large stochastic expansion models are often solved with Benders or scenario decomposition, which avoid forming the complete problem at once.
PDLP takes a different approach. It operates on the constraint matrix and a small set of vectors without constructing matrix factors, making the monolithic deterministic equivalent worth testing directly.
Experiments
We compared a GPU first-order method against a CPU interior-point method on identical inputs.
| NVIDIA cuOpt | HiGHS | |
|---|---|---|
| Version | 26.8.0 | 1.15.1 |
| Method | PDLP, first-order | HiPO, interior point |
| Hardware | NVIDIA B200 (1/2/4/8 GPUs) | 20-core Arm processor (Cortex-X925 + Cortex-A725) |
Both solvers ran to a relative tolerance of 1e-6 under a four-hour limit, and both received the same deterministic equivalent problem, which is available here. Decomposition methods such as Benders and scenario decomposition are deliberately out of scope. What we want to know here is how solvers behave on identical inputs.
Each solver ran on its target architecture, so this is not a like-for-like hardware comparison. The algorithms are also very different. We use HiGHS as a reproducible CPU reference because it is open source, reliable across operating systems, and available through several language interfaces. Its other interior-point option, IPX, reached the time limit even on the smallest instance.
Results
cuOpt solved every instance in the set, including the 100-scenario deterministic equivalent with 471 million nonzeros. At 1e-6, the single-GPU configuration did not finish the two largest cases within four hours. Eight GPUs solved the 50-scenario instance in 3,337 seconds and the 100-scenario instance in 2,492 seconds.
HiGHS solved the five instances from one to ten scenarios, while the 20-scenario run reached the four-hour limit on the available hardware. Across the completed runs, the speedup varies with scenario count: against a single GPU it ranges from 1.7x to 6x. It stays consistently in the GPU’s favor, and it widens against eight GPUs, reaching 14.9x at ten scenarios. At 50 and 100 scenarios, only multi-GPU cuOpt completed within the four-hour limit.
TThe completed solves also provide an important accuracy check. Across all five instances with a CPU solution, HiGHS and cuOpt at 1e-6 reach objectives within 0.03% of each other, and on the one-scenario instance they agree to the displayed digit. A first-order method invites the question of whether the fast answer is also the right one. The close objective values suggest that both solvers reached solutions of similar quality.
The five-scenario instance also benefited from multiple GPUs. Four GPUs reduced solve time from 1,724 to 556 seconds, a 3.1x speedup; eight GPUs were slightly slower at 589 seconds. Moving from one to eight GPUs reduced solve time by 2.5x at 10 scenarios and 4.2x at 20 scenarios.
At 50 and 100 scenarios, the single-GPU runs did not converge within the four-hour limit, so we do not report speedup ratios.
GPU optimization adds another practical option alongside barrier solvers and decomposition. Barrier methods remain highly competitive, while decomposition is still essential for many planning models. For large deterministic equivalents, the results show that a multi-GPU first-order solve is a credible third option.
Although outside the scope of these experiments, Benders decomposition offers another application: a mixed-integer master problem chooses investments, while scenario LPs evaluate operation and provide cuts to refine those decisions at each iteration. cuOpt’s batch PDLP can solve related LPs together on one GPU when they share a constraint matrix, while multi-GPU PDLP can distribute a large LP across GPUs. These capabilities could accelerate the repeated LP solves that drive a Benders workflow.
Solving the full deterministic equivalent directly avoids the additional implementation and tuning work that decomposition approaches, like Benders or scenario decomposition, can require. For energy planners, shorter solve times could make room for more scenarios, finer temporal resolution, and more sensitivity runs within the same study window. That would let teams explore a wider range of conditions before committing to an investment plan.
