MadNLP.jl: Julia Solver for Large-Scale NLP

• MadNLP.jl is a Julia-based nonlinear programming solver designed for large-scale, sparse optimization problems that exploit graph structures and physics constraints.
• It implements advanced interior-point methods—including Schur decomposition and GPU-optimized condensed-space formulations—to significantly accelerate solution times.
• The solver integrates seamlessly with Julia’s ecosystem (e.g., Plasmo.jl and ExaModels.jl) to support diverse workflows in energy systems, optimal control, and generative sampling.

MadNLP.jl is a Julia-based nonlinear programming solver for large-scale sparse optimization problems. In the published literature it is presented as a general-purpose NLP solver, a graph-aware solver for networked and space–time models, a condensed-space and GPU-enabled interior-point solver for optimal power flow and optimal control, an interior-point subproblem engine inside augmented-Lagrangian methods for degenerate problems, and the solver used to enforce hard physics constraints in projection-based generative sampling (Shin et al., 2020, Cole et al., 2022, Shin et al., 2023, Montoison et al., 7 Oct 2025, Kolli et al., 30 Jun 2026).

1. Origins, scope, and design goals

MadNLP.jl was introduced as a general-purpose nonlinear programming solver in work on graph-based modeling and decomposition of energy infrastructures, where it was interfaced to Plasmo.jl and used to accelerate transient gas network optimization and multi-period AC optimal power flow (Shin et al., 2020). From the outset, the solver was designed for large, sparse, physics-based models rather than for small dense problems. The stated target class includes graph-structured models with tens or hundreds of thousands of variables and constraints, and later work explicitly frames the solver as applicable to “millions of variables/constraints” when structure can be exposed effectively (Shin et al., 2020).

A recurrent design theme is that MadNLP.jl is not treated merely as a monolithic sparse solver. Instead, the literature presents it as a solver framework whose behavior depends strongly on how problem structure is exposed. In the graph-structured setting, the relevant structure is the decomposition of a model into local node problems plus linking constraints; in GPU-oriented work, it is the repeated sparse pattern induced by SIMD-style formulations, batched samples, multi-period dynamics, or direct transcription over time grids (Cole et al., 2022, Shin et al., 2023, Montoison et al., 4 Oct 2025).

The reported empirical motivation is substantial. In the early energy-infrastructure study, solution times were reduced from 72.36 sec to 23.84 sec for the gas problem and from 515.81 sec to 149.45 sec for the power flow problem, with the paper characterizing the gain as over 300% improvement in wall-clock time (Shin et al., 2020). Later work extends that scope from multicore CPU decomposition to GPU-resident sparse NLP, high-memory multi-period AC OPF, and physics-constrained generative modeling (Shin et al., 2024, Kolli et al., 30 Jun 2026).

2. Nonlinear-programming formulation and interior-point machinery

At its core, MadNLP.jl solves nonlinear programs of the form

$$\begin{aligned} \min_x \quad & f(x) \ \text{s.t.}\quad & c(x)=0,\ & x \ge 0, \end{aligned}$$

or, in bounded form,

$$\min_z\quad f(z)\qquad \text{s.t.}\quad g(z)=0,\qquad h_L \le h(z) \le h_U,$$

with inequalities handled through slacks and barrier terms (Shin et al., 2020, Montoison et al., 4 Oct 2025).

The solver implements a primal–dual interior-point method with a filter line-search, in the same spirit as IPOPT and Wächter–Biegler-type algorithms (Shin et al., 2020, Montoison et al., 4 Oct 2025). In barrier form, one representative formulation is

$$\min_x \; \varphi(x):= f(x)-\mu\,\mathbf{e}^\top \log(x) \quad \text{s.t.}\quad c(x)=0,$$

with KKT conditions linearized into a Newton system of the form

$$\begin{bmatrix} W + \Sigma + \delta_w I & A^\top \ A & -\delta_c I \end{bmatrix} \begin{bmatrix} d^x\ d^\lambda \end{bmatrix} = - \begin{bmatrix} \nabla_x \varphi(x)+A^\top \lambda\ c(x) \end{bmatrix},$$

where $W=\nabla^2_{xx}\mathcal{L}(x,\lambda,z)$, $\Sigma=X^{-1}Z$, and $\delta_w,\delta_c>0$ are regularization parameters (Shin et al., 2020).

Later GPU-oriented and optimal-control papers retain the same basic interior-point logic but emphasize several implementation details: Newton steps are computed from symmetric indefinite or condensed symmetric positive definite KKT systems; second-order information is used explicitly through the Hessian of the Lagrangian; primal–dual regularization with inertia correction is employed to preserve numerical stability and the correct inertia; and optional iterative refinement can be applied when condensed systems become ill-conditioned near convergence (Shin et al., 2023, Pacaud et al., 2024, Montoison et al., 4 Oct 2025).

The algorithmic role of MadNLP.jl therefore spans both classical large-scale NLP and more specialized formulations. In some applications it solves the original constrained problem directly; in others it solves equality-relaxed or augmented-Lagrangian subproblems whose KKT systems have been altered to better match GPU sparse factorization or degeneracy-robust outer algorithms (Shin et al., 2023, Montoison et al., 7 Oct 2025).

3. Graph structure, decomposition, and KKT-system organization

A defining contribution of the early MadNLP literature is the treatment of graph-structured optimization. In this setting, the model is written over a graph $G=(V,E)$, with local variables $x_i$, local objective terms $f_i(x_i)$, inner constraints $$\min_z\quad f(z)\qquad \text{s.t.}\quad g(z)=0,\qquad h_L \le h(z) \le h_U,$$0, and linking constraints $$\min_z\quad f(z)\qquad \text{s.t.}\quad g(z)=0,\qquad h_L \le h(z) \le h_U,$$1 coupling neighboring nodes (Shin et al., 2020). Plasmo.jl exposes this problem graph explicitly, while MadNLP.jl uses the resulting structure to organize function evaluation and linear algebra (Shin et al., 2020, Cole et al., 2022).

For graph-structured models, the Newton systems inherit a block-bordered structure: node-local Hessian and Jacobian blocks appear on the diagonal, while linking constraints create a border block. MadNLPGraph.jl exploits this through Schur complement decomposition. The local KKT blocks are factorized independently and in parallel, contributions to a smaller Schur complement are accumulated, the Schur system is solved, and node-local directions are recovered by back-substitution (Cole et al., 2022). This is especially effective in two-stage stochastic or multi-block structures where the number of linking constraints remains modest even as the number of scenarios or local subproblems grows (Cole et al., 2022).

A second structure-exploiting approach is the restricted additive Schwarz scheme for KKT systems. In that formulation, the graph is partitioned into non-overlapping subdomains, enlarged by an overlap parameter $$\min_z\quad f(z)\qquad \text{s.t.}\quad g(z)=0,\qquad h_L \le h(z) \le h_U,$$2, and the global KKT system is approximated through parallel local solves

$$\min_z\quad f(z)\qquad \text{s.t.}\quad g(z)=0,\qquad h_L \le h(z) \le h_U,$$3

When $$\min_z\quad f(z)\qquad \text{s.t.}\quad g(z)=0,\qquad h_L \le h(z) \le h_U,$$4, the method reduces to block-Jacobi; when the overlap is maximal, it degenerates to one large direct solve (Shin et al., 2020). The reported purpose is to bridge centralized and decentralized solution strategies while reusing local sparse factorizations and Julia multithreading (Shin et al., 2020).

The same structural viewpoint reappears in later non-graph settings. Batched physics-constrained generative sampling yields Jacobians and KKT systems that are block-diagonal across samples; direct transcription of optimal control produces banded or block-banded KKT structure along time; and multi-period AC OPF yields repeated per-period constraint blocks plus simple inter-temporal couplings (Montoison et al., 4 Oct 2025, Kolli et al., 30 Jun 2026). Across these domains, MadNLP.jl is consistently presented as a solver that benefits when sparsity and decomposition are preserved at the modeling level rather than obscured by dense tensor algebra.

4. Condensed-space methods and GPU-oriented formulations

A major later development is the introduction of condensed-space interior-point methods designed specifically for GPUs. The motivating observation is that classical sparse symmetric indefinite KKT factorization with numerical pivoting is difficult to parallelize on GPUs, whereas condensed symmetric positive definite systems can be factorized by Cholesky or LDL$$\min_z\quad f(z)\qquad \text{s.t.}\quad g(z)=0,\qquad h_L \le h(z) \le h_U,$$5 with static pivoting (Shin et al., 2023, Pacaud et al., 2024).

One formulation relaxes equality constraints through additional slacks:

$$\min_z\quad f(z)\qquad \text{s.t.}\quad g(z)=0,\qquad h_L \le h(z) \le h_U,$$6

with $$\min_z\quad f(z)\qquad \text{s.t.}\quad g(z)=0,\qquad h_L \le h(z) \le h_U,$$7 and $$\min_z\quad f(z)\qquad \text{s.t.}\quad g(z)=0,\qquad h_L \le h(z) \le h_U,$$8 in the lifted-KKT approach (Shin et al., 2023, Shin et al., 2024). Eliminating slack and dual blocks yields a condensed system in the primal variables,

$$\min_z\quad f(z)\qquad \text{s.t.}\quad g(z)=0,\qquad h_L \le h(z) \le h_U,$$9

which is enforced to be symmetric positive definite through inertia-correction regularization (Shin et al., 2023, Shin et al., 2024).

A second GPU-oriented study compares two condensed-space methods implemented inside MadNLP.jl: HyKKT and LiftedKKT. HyKKT augments the condensed matrix with a term $$\min_x \; \varphi(x):= f(x)-\mu\,\mathbf{e}^\top \log(x) \quad \text{s.t.}\quad c(x)=0,$$0 and solves the equality block through a hybrid direct/iterative strategy, while LiftedKKT replaces equalities by tight inequalities and solves a fully direct condensed SPD system (Pacaud et al., 2024). Both formulations are implemented through MadNLP’s AbstractKKTSystem abstraction and interfaced with NVIDIA cuDSS for sparse LDL$$\min_x \; \varphi(x):= f(x)-\mu\,\mathbf{e}^\top \log(x) \quad \text{s.t.}\quad c(x)=0,$$1/Cholesky factorization on GPUs (Pacaud et al., 2024).

GPU optimal-control work extends the same pattern to direct transcription. ExaModels.jl compiles objectives, constraints, gradients, Jacobians, and Hessians into GPU kernels; MadNLP.jl assembles the KKT matrix directly on the GPU; symbolic analysis and reordering are done once on the host; and all interior-point iterations then run on the GPU via cuDSS (Montoison et al., 4 Oct 2025). In the multi-period AC OPF setting, the same strategy is combined with high-memory GH200 hardware and unified memory to solve instances with more than 10 million variables (Shin et al., 2024).

These formulations do not replace the classical interior-point implementation; they coexist with it. The literature therefore describes MadNLP.jl as a solver family supporting unreduced, augmented, condensed, lifted, hybrid, and application-specific KKT systems, selected according to hardware constraints, conditioning, and problem structure rather than according to a single fixed algorithm (Pacaud et al., 2024, Montoison et al., 7 Oct 2025).

5. Software ecosystem and compositional workflows

MadNLP.jl is embedded in a broader Julia ecosystem in which modeling, automatic differentiation, structure exposure, and sparse linear algebra are split across specialized packages. In graph-based work, Plasmo.jl provides OptiGraph models on top of JuMP.jl, and MadNLP consumes local node oracles together with graph partitions (Shin et al., 2020, Cole et al., 2022). In GPU work, ExaModels.jl serves as the modeling and AD layer, generating GPU-resident sparse derivatives that MadNLP then uses for KKT assembly and sparse factorization (Shin et al., 2023, Pacaud et al., 2024).

The same compositional design extends to domain-specific layers. ExaModelsPower.jl targets static and multi-period AC OPF models on GPUs; OptimalControl.jl transcribes continuous-time control problems into ExaModels NLPs; DynamicNLPModels.jl builds sparse or dense reduced linear-quadratic MPC models; and MadNCL.jl wraps MadNLP inside Algorithm NCL, using custom KKT-system types K2rAuglagKKTSystem and K1sAuglagKKTSystem (Shin et al., 2024, Montoison et al., 4 Oct 2025, Cole et al., 2022, Montoison et al., 7 Oct 2025).

Setting	Modeling layer	MadNLP.jl role
Graph-structured energy systems	Plasmo.jl + JuMP.jl	Interior-point solver with Schur or RAS decomposition
AC OPF and multi-period AC OPF	ExaModels.jl + ExaModelsPower.jl	Condensed/lifted GPU NLP solver with cuDSS
Optimal control and MPC	OptimalControl.jl or DynamicNLPModels.jl	Sparse NLP solver or dense condensed MPC solver
Degenerate NLPs	MadNCL.jl wrapper	Interior-point subproblem solver inside augmented Lagrangian
Physics-constrained generative sampling	ExaModels.jl in PCFM/SNAP-FM	Projection NLP solver enforcing hard constraints

This ecosystem-level organization is technically significant because it lets the same solver participate in markedly different workflows without changing the basic user-facing optimization model. A Plasmo OptiGraph, an ExaModels GPU model, a control problem discretized by OptimalControl.jl, or a PCFM projection model compiled by ExaModels can all become MadNLP problems, but with different KKT formulations, linear solvers, and hardware backends selected by multiple dispatch and solver configuration (Shin et al., 2020, Montoison et al., 4 Oct 2025).

6. Application domains and empirical record

The energy-infrastructure papers establish the original empirical profile of MadNLP.jl. In graph-based transient gas network optimization and multi-period AC optimal power flow, the combined Plasmo–MadNLP framework reduced solution times from 72.36 sec to 23.84 sec and from 515.81 sec to 149.45 sec, respectively (Shin et al., 2020). A later graph-structured framework study scaled to a stochastic gas network problem containing over 1.7 million variables and reported that the Schur-based solver achieved a lower time per interior-point iteration than monolithic MA57 or PardisoMKL at the largest scales (Cole et al., 2022).

In AC optimal power flow, the GPU-oriented results are more pronounced. The ACOPF study based on ExaModels.jl and condensed-space IPM reports that MadNLP.jl and ExaModels.jl on NVIDIA GPUs achieve an order of magnitude speedup compared with state-of-the-art tools running on contemporary CPUs (Shin et al., 2023). The condensed-space comparison on PGLIB and COPS reports that GPUs can attain up to a tenfold speed increase compared to CPUs when solving large-scale instances (Pacaud et al., 2024). For multi-period AC OPF on an NVIDIA GH200, the reported headline result is a multi-period OPF instance with more than 10 million variables solved up to $$\min_x \; \varphi(x):= f(x)-\mu\,\mathbf{e}^\top \log(x) \quad \text{s.t.}\quad c(x)=0,$$2 precision in less than 10 minutes, and the same study reports successful runs with over 20 million constraints (Shin et al., 2024).

In optimal control, MadNLP.jl functions as the numerical engine for sparse NLPs produced by direct transcription. The GPU workflow paper reports that, for the Goddard problem, the GPU becomes faster after approximately $$\min_x \; \varphi(x):= f(x)-\mu\,\mathbf{e}^\top \log(x) \quad \text{s.t.}\quad c(x)=0,$$3 on A100 or $$\min_x \; \varphi(x):= f(x)-\mu\,\mathbf{e}^\top \log(x) \quad \text{s.t.}\quad c(x)=0,$$4 on H100, and that the largest H100 instance reaches approximately 2 million unknowns in approximately 15 s (Montoison et al., 4 Oct 2025). For the Quadrotor problem, the crossover occurs near $$\min_x \; \varphi(x):= f(x)-\mu\,\mathbf{e}^\top \log(x) \quad \text{s.t.}\quad c(x)=0,$$5 on A100 or $$\min_x \; \varphi(x):= f(x)-\mu\,\mathbf{e}^\top \log(x) \quad \text{s.t.}\quad c(x)=0,$$6 on H100, with the largest H100 problem at approximately $$\min_x \; \varphi(x):= f(x)-\mu\,\mathbf{e}^\top \log(x) \quad \text{s.t.}\quad c(x)=0,$$7 variables+constraints solved in approximately 13 s (Montoison et al., 4 Oct 2025). In linear-quadratic MPC, MadNLP is used as a condensed-space primal–dual interior-point solver after state elimination; for a 3D temperature-control test at $$\min_x \; \varphi(x):= f(x)-\mu\,\mathbf{e}^\top \log(x) \quad \text{s.t.}\quad c(x)=0,$$8 and $$\min_x \; \varphi(x):= f(x)-\mu\,\mathbf{e}^\top \log(x) \quad \text{s.t.}\quad c(x)=0,$$9, the paper reports CPU sparse $$\begin{bmatrix} W + \Sigma + \delta_w I & A^\top \ A & -\delta_c I \end{bmatrix} \begin{bmatrix} d^x\ d^\lambda \end{bmatrix} = - \begin{bmatrix} \nabla_x \varphi(x)+A^\top \lambda\ c(x) \end{bmatrix},$$0 s total, CPU dense $$\begin{bmatrix} W + \Sigma + \delta_w I & A^\top \ A & -\delta_c I \end{bmatrix} \begin{bmatrix} d^x\ d^\lambda \end{bmatrix} = - \begin{bmatrix} \nabla_x \varphi(x)+A^\top \lambda\ c(x) \end{bmatrix},$$1 s, and GPU dense $$\begin{bmatrix} W + \Sigma + \delta_w I & A^\top \ A & -\delta_c I \end{bmatrix} \begin{bmatrix} d^x\ d^\lambda \end{bmatrix} = - \begin{bmatrix} \nabla_x \varphi(x)+A^\top \lambda\ c(x) \end{bmatrix},$$2 s (Cole et al., 2022).

In degenerate nonlinear programming, MadNCL uses MadNLP as its interior-point subproblem solver and reports strong robustness on MPCC-like security-constrained OPF instances where Ipopt and plain MadNLP fail, while also achieving large GPU speedups for stabilized KKT solves via cuDSS (Montoison et al., 7 Oct 2025). This usage positions MadNLP.jl not only as a standalone solver but also as a reusable interior-point engine inside more elaborate augmented-Lagrangian algorithms.

In physics-constrained generative modeling, MadNLP.jl solves the projection NLP

$$\begin{bmatrix} W + \Sigma + \delta_w I & A^\top \ A & -\delta_c I \end{bmatrix} \begin{bmatrix} d^x\ d^\lambda \end{bmatrix} = - \begin{bmatrix} \nabla_x \varphi(x)+A^\top \lambda\ c(x) \end{bmatrix},$$3

at every PCFM sampling step, enforcing PDE-derived hard constraints without retraining the generative model (Kolli et al., 30 Jun 2026). The paper’s proposed stack, ExaModels + MadNLP (GPU), reports 11.73 ± 0.34 s with infeasibility $$\begin{bmatrix} W + \Sigma + \delta_w I & A^\top \ A & -\delta_c I \end{bmatrix} \begin{bmatrix} d^x\ d^\lambda \end{bmatrix} = - \begin{bmatrix} \nabla_x \varphi(x)+A^\top \lambda\ c(x) \end{bmatrix},$$4 for the heat equation with simple constraints, 110.26 ± 1.92 s with infeasibility 0.222 for a nonlinear heat-equation constraint set where CPU variants either take 3501.11 ± 89.36 s or do not finish, and 38.83 ± 0.55 s with infeasibility 0.920 for Burgers’ equation with flux constraints versus 2048.83 ± 187.58 s on CPU (Kolli et al., 30 Jun 2026). The same paper reports approximately linear runtime growth with sample count in the scaling study, consistent with block-diagonal batching structure (Kolli et al., 30 Jun 2026).

7. Limitations, numerical trade-offs, and open directions

The reported advantages of MadNLP.jl are strongly structure-dependent. The graph-based decomposition papers explicitly state that gains are largest when the problem is graph-structured and sparse, that dense fully coupled problems are less suitable, and that RAS overhead only pays off for sufficiently large problems (Shin et al., 2020). They also note sensitivity to partition quality, communication overhead, and overlap-size tuning, with spatial decomposition identified as future work (Shin et al., 2020).

GPU-oriented formulations introduce a separate trade-off between parallel factorization and attainable accuracy. The multi-period AC OPF study adopts $$\begin{bmatrix} W + \Sigma + \delta_w I & A^\top \ A & -\delta_c I \end{bmatrix} \begin{bmatrix} d^x\ d^\lambda \end{bmatrix} = - \begin{bmatrix} \nabla_x \varphi(x)+A^\top \lambda\ c(x) \end{bmatrix},$$5 and states that convergence is reliable only up to $$\begin{bmatrix} W + \Sigma + \delta_w I & A^\top \ A & -\delta_c I \end{bmatrix} \begin{bmatrix} d^x\ d^\lambda \end{bmatrix} = - \begin{bmatrix} \nabla_x \varphi(x)+A^\top \lambda\ c(x) \end{bmatrix},$$6, lower than $$\begin{bmatrix} W + \Sigma + \delta_w I & A^\top \ A & -\delta_c I \end{bmatrix} \begin{bmatrix} d^x\ d^\lambda \end{bmatrix} = - \begin{bmatrix} \nabla_x \varphi(x)+A^\top \lambda\ c(x) \end{bmatrix},$$7 on CPUs (Shin et al., 2024). By contrast, the condensed-space comparison reports current implementations reliable down to about $$\begin{bmatrix} W + \Sigma + \delta_w I & A^\top \ A & -\delta_c I \end{bmatrix} \begin{bmatrix} d^x\ d^\lambda \end{bmatrix} = - \begin{bmatrix} \nabla_x \varphi(x)+A^\top \lambda\ c(x) \end{bmatrix},$$8, while still emphasizing the need for better modified Cholesky/LDL$$\begin{bmatrix} W + \Sigma + \delta_w I & A^\top \ A & -\delta_c I \end{bmatrix} \begin{bmatrix} d^x\ d^\lambda \end{bmatrix} = - \begin{bmatrix} \nabla_x \varphi(x)+A^\top \lambda\ c(x) \end{bmatrix},$$9 strategies and customized GPU sparse factorizations for tighter tolerances (Pacaud et al., 2024). A plausible implication is that the practical accuracy ceiling is formulation-dependent and remains an active numerical issue rather than a settled property of the solver as a whole.

Memory is another recurrent constraint. The GH200 study shows that very large multi-period AC OPF instances require high-memory GPUs and that A100-class hardware runs out of memory on the biggest benchmarks (Shin et al., 2024). The SNAP-FM study similarly treats GPU memory as an implicit limitation for larger PDE grids and higher-dimensional problems, noting sparse-factorization overhead beyond raw matrix storage (Kolli et al., 30 Jun 2026). On the algorithmic side, the MadNCL work observes that condensed systems can become extremely ill-conditioned as penalty parameters increase, while stabilized KKT formulations require careful tuning of static-pivot LDL$W=\nabla^2_{xx}\mathcal{L}(x,\lambda,z)$0, pivot epsilon, and iterative refinement (Montoison et al., 7 Oct 2025).

These limitations define the main research directions visible across the literature: improved robustness for highly nonlinear and degenerate constraints, more memory-efficient GPU sparse factorizations, better user-facing interfaces for constrained scientific machine learning workflows, extensions of graph-based decomposition beyond temporal settings, and continued development of KKT formulations that are either symmetric quasi-definite or symmetric positive definite under modest regularization (Shin et al., 2020, Montoison et al., 7 Oct 2025, Kolli et al., 30 Jun 2026). Taken together, the published record suggests that MadNLP.jl is best understood not as a single immutable solver configuration but as a rapidly evolving solver platform for large-scale structured nonlinear optimization.