MadNCL: Robust Augmented-Lagrangian Solver

• MadNCL is a robust augmented-Lagrangian framework that reformulates constrained problems with auxiliary residual variables, enhancing stability in degenerate optimization scenarios.
• It integrates interior-point methods and warm-start strategies to efficiently solve nonlinearly constrained subproblems, even when standard constraint qualifications like LICQ fail.
• GPU implementations of MadNCL employ specialized KKT formulations and sparse LDL factorization, yielding significant speedups on large-scale benchmarks and power systems analyses.

MadNCL denotes implementations of Algorithm NCL, a nonlinearly constrained augmented-Lagrangian method for smooth constrained optimization. In the 2021 Julia work, the implementation is called NCL.jl, and the paper’s “MadNCL” wording refers to that Julia codebase; in 2025, MadNCL designates a GPU implementation of Algorithm NCL and, in a closely related power-systems study, a GPU-accelerated solver for security-constrained AC optimal power flow. Across these usages, the common core is the solution of large-scale nonlinear programs—especially degenerate ones in which standard constraint qualifications such as LICQ or MFCQ may fail—by reformulating augmented-Lagrangian subproblems with auxiliary residual variables and solving the resulting nonlinearly constrained subproblems by interior-point methods (Ma et al., 2021, Montoison et al., 7 Oct 2025, Pacaud et al., 15 Oct 2025).

1. Naming, lineage, and scope

The name MadNCL is attached to a family of NCL-based implementations rather than to a single invariant software artifact. The 2021 paper presents a Julia implementation of Algorithm NCL under the package name NCL.jl, with the “MadNCL” wording referring to that implementation. The 2025 GPU paper uses MadNCL for a GPU implementation of Algorithm NCL for large-scale, degenerate nonlinear programs. The 2025 AC-SCOPF paper presents MadNCL.jl as a GPU-accelerated augmented-Lagrangian solver specialized to security-constrained AC optimal power flow in polar form (Ma et al., 2021, Montoison et al., 7 Oct 2025, Pacaud et al., 15 Oct 2025).

Paper	Implementation/name	Reported focus
(Ma et al., 2021)	NCL.jl	Julia implementation of Algorithm NCL for smooth constrained optimization, including LICQ-violating cases
(Montoison et al., 7 Oct 2025)	MadNCL	GPU implementation of Algorithm NCL for large-scale, degenerate nonlinear programs
(Pacaud et al., 15 Oct 2025)	MadNCL.jl	GPU-accelerated NCL solver for large-scale corrective AC-SCOPF with complementarity constraints

Algorithmically, the lineage runs through the classic LANCELOT augmented-Lagrangian framework. The 2021 formulation states that Algorithm NCL is equivalent to the LANCELOT augmented Lagrangian method, but recasts it as a short sequence of nonlinearly constrained subproblems that can be handled efficiently by IPOPT and KNITRO with warm starts (Ma et al., 2021). The later GPU work preserves this augmented-Lagrangian structure while redesigning the KKT linear algebra for sparse GPU execution (Montoison et al., 7 Oct 2025).

A common misconception is to treat MadNCL as merely an interior-point solver. The cited works instead present it as an augmented-Lagrangian method whose subproblems are solved by interior-point techniques. This distinction is central because the robustness claims are tied to the outer augmented-Lagrangian mechanism, not to replacing one IPM by another (Montoison et al., 7 Oct 2025).

2. Core optimization model and NCL reformulation

The baseline problem class is the smooth constrained program

$$\min_x \ \phi(x) \qquad \text{subject to} \qquad c(x)=0,\ \ell \le x \le u,$$

where $\phi$ and $c$ are twice continuously differentiable, and first and second derivatives are available (Ma et al., 2021, Montoison et al., 7 Oct 2025). The key pathological regime is constraint degeneracy. In the 2021 presentation, if $J(x)$ is the constraint Jacobian, LICQ requires $\operatorname{rank}(J(x))=m$ at the solution; many standard NLP solvers have difficulty when this fails (Ma et al., 2021).

The classical augmented-Lagrangian subproblem minimizes

$\phi(x) - y_k^T c(x) + \tfrac12 \rho_k \|c(x)\|^2$

over the bounds $\ell \le x \le u$. Algorithm NCL introduces an auxiliary variable $r\equiv -c(x)$ and instead solves

$$\min_{x,r}\ \phi(x) + y_k^T r + \tfrac12 \rho_k \|r\|^2 \quad \text{subject to} \quad c(x)+r=0,\ \ell \le x \le u.$$

This reformulation is mathematically equivalent to the classical ALM subproblem, but numerically advantageous because the subproblem constraint Jacobian becomes

$[J(x)\ \ I],$

which is structurally full row rank even if the original Jacobian $\phi$0 is rank deficient (Ma et al., 2021, Montoison et al., 7 Oct 2025).

That structural property explains the method’s role in degenerate optimization. Interior-point solvers such as IPOPT and KNITRO can exploit the regularized subproblem more stably than when applied directly to the original model. The 2025 GPU paper makes the same point in different language: the added variables $\phi$1 regularize the constraint system and make the constraints linearly independent, while the quadratic penalty acts as a dual regularizer (Montoison et al., 7 Oct 2025).

This suggests that the defining idea of MadNCL is not the use of GPUs, Julia, or a specific linear solver, but the NCL transformation itself: a residual-variable augmentation that restores a favorable constraint structure without changing the underlying optimization problem.

3. Major-iteration policy, stopping rules, and warm starts

In the 2021 implementation, each major iteration solves an NCL subproblem approximately and returns

$\phi$2

where $\phi$3 are bound multipliers (Ma et al., 2021). The major-iteration logic is driven by feasibility. With

$\phi$4

the implementation uses $\phi$5 as feasibility tolerance, $\phi$6 as optimality tolerance, $\phi$7 as penalty parameter, and $\phi$8 as the initial barrier parameter for IPOPT or KNITRO (Ma et al., 2021).

If

$\phi$9

the method updates the multipliers and tightens the subproblem tolerances: $c$0 Otherwise it keeps $c$1 fixed and increases the penalty: $c$2 with $c$3 in the implementation (Ma et al., 2021). A natural stopping test is

$c$4

since $c$5 makes the subproblem equivalent to the original constrained problem (Ma et al., 2021).

The literature uses two sign conventions for the multiplier update. The 2021 paper writes $c$6, while the 2025 GPU paper gives the equivalent convention $c$7 because its outer update is written in terms of $c$8 rather than $c$9 (Ma et al., 2021, Montoison et al., 7 Oct 2025). This is a notational difference, not an algorithmic one.

Warm-starting is treated as essential rather than optional. For subproblem $J(x)$0, the 2021 solver is initialized from the previous solution $J(x)$1. For IPOPT, the implementation uses warm_start_init_point=yes for $J(x)$2, with mu_init reduced through the sequence $J(x)$3. For KNITRO, it uses algorithm=1 for $J(x)$4, and for $J(x)$5 employs bar_directinterval=0, bar_initpt=2, bar_murule=1, together with analogous reductions of bar_initmu and bar_slackboundpush (Ma et al., 2021).

The practical consequences are explicit in the reported experiments: there are typically only about 10 major iterations; warm starts make later subproblems cheaper; and gradually tightening $J(x)$6 from $J(x)$7 to $J(x)$8 saves work (Ma et al., 2021).

4. GPU reformulation and KKT system design

The 2025 GPU implementation extends the NCL idea from a robust algorithmic framework into a GPU-oriented solver architecture. MadNCL fuses the outer ALM loop and the inner IPM loop, rather than treating them as separate phases, and adds an extrapolation step intended to improve efficiency and convergence stability (Montoison et al., 7 Oct 2025).

With primal-dual variables $J(x)$9, the inner iterations solve the smooth system

$\operatorname{rank}(J(x))=m$0

where the components consist of the stationarity equations, the residual equation $\operatorname{rank}(J(x))=m$1, and the barrier-perturbed complementarity conditions for the bound constraints (Montoison et al., 7 Oct 2025). Newton’s method is applied to

$\operatorname{rank}(J(x))=m$2

with a filter line search. Near convergence, MadNCL attempts an extrapolation step

$\operatorname{rank}(J(x))=m$3

accepted when

$\operatorname{rank}(J(x))=m$4

The paper attributes the fast asymptotic behavior and superlinear convergence in the final phase to this mechanism (Montoison et al., 7 Oct 2025).

A central contribution is the comparison of several KKT formulations. Starting from the full Newton system, eliminating the bound multipliers yields the classical symmetric augmented KKT system

$\operatorname{rank}(J(x))=m$5

with $\operatorname{rank}(J(x))=m$6 (Montoison et al., 7 Oct 2025). For GPU execution, the paper introduces two alternatives:

1. Stabilized KKT system $\operatorname{rank}(J(x))=m$7,

$\operatorname{rank}(J(x))=m$8

obtained by eliminating $\operatorname{rank}(J(x))=m$9.

2. Condensed KKT system $\phi(x) - y_k^T c(x) + \tfrac12 \rho_k \|c(x)\|^2$0,

$\phi(x) - y_k^T c(x) + \tfrac12 \rho_k \|c(x)\|^2$1

which is SPD after regularization (Montoison et al., 7 Oct 2025).

The inertia analysis links these systems: $\phi(x) - y_k^T c(x) + \tfrac12 \rho_k \|c(x)\|^2$2 Thus, correct inertia of the stabilized form or positive definiteness of the condensed form guarantees a valid descent direction for the filter line search (Montoison et al., 7 Oct 2025).

On the GPU, linear algebra is performed with sparse $\phi(x) - y_k^T c(x) + \tfrac12 \rho_k \|c(x)\|^2$3 factorization using static pivoting in NVIDIA cuDSS. For $\phi(x) - y_k^T c(x) + \tfrac12 \rho_k \|c(x)\|^2$4, the implementation uses cuDSS with a small pivot epsilon, e.g. $\phi(x) - y_k^T c(x) + \tfrac12 \rho_k \|c(x)\|^2$5, and then recovers the original-system solution by iterative refinement, typically Richardson iterations. For $\phi(x) - y_k^T c(x) + \tfrac12 \rho_k \|c(x)\|^2$6, cuDSS uses default settings because the matrix is positive definite after regularization (Montoison et al., 7 Oct 2025). MadNCL uses MadNLP as the inner IPM solver and implements $\phi(x) - y_k^T c(x) + \tfrac12 \rho_k \|c(x)\|^2$7 and $\phi(x) - y_k^T c(x) + \tfrac12 \rho_k \|c(x)\|^2$8 as specialized AbstractKKTSystem variants, keeping allocations and computations on the GPU (Montoison et al., 7 Oct 2025).

The numerical interpretation is deliberately nuanced. The stabilized $\phi(x) - y_k^T c(x) + \tfrac12 \rho_k \|c(x)\|^2$9 formulation is described as more robust, whereas the condensed $\ell \le x \le u$0 formulation can be faster when conditioning is favorable. The condensed system can also become too ill-conditioned for OPF-like problems as $\ell \le x \le u$1 grows, so MadNCL does not claim a universally dominant GPU linearization (Montoison et al., 7 Oct 2025).

5. Empirical behavior on general nonlinear programming benchmarks

The 2021 study evaluates the Julia implementation on tax policy models, nonlinear least-squares problems, and large constrained problems from CUTEst. The tax policy models are particularly significant because they may not satisfy LICQ at the solution and involve many active constraints. The test family ranges from 1D to 5D, with the largest instance having $\ell \le x \le u$2 variables and $\ell \le x \le u$3 nonlinear constraints (Ma et al., 2021).

On these tax models, the reported pattern is that standard solvers can struggle badly, while NCL with KNITRO or IPOPT solves the problems robustly. The logs show early subproblems solved at low accuracy, monotonic decrease of $\ell \le x \le u$4, decrease of $\ell \le x \le u$5, and later subproblems that become much cheaper because of warm starts. Direct IPOPT can time out or fail on large instances, and NCL/KNITRO often outperforms direct KNITRO on these difficult models (Ma et al., 2021).

For nonlinear least squares, the paper rewrites

$\ell \le x \le u$6

as the feasibility problem

$\ell \le x \le u$7

With $\ell \le x \le u$8, $\ell \le x \le u$9, and a single major iteration, the first NCL subproblem essentially reproduces the least-squares objective structure (Ma et al., 2021). On 127 CUTEst NLS problems, KNITRO’s native nonlinear least-squares method solved 101, whereas NCL solved 119. Performance profiles show fewer residual evaluations, fewer Jacobian evaluations, lower runtime, and better overall robustness, despite NCL evaluating Hessians (Ma et al., 2021).

On 166 constrained CUTEst problems with at least 100 variables and nonlinear constraints, KNITRO solved 147 and IPOPT solved 126. The paper emphasizes that NCL is not always the fastest on general problems, but succeeds on some cases where plain IPOPT or KNITRO fail (Ma et al., 2021). This is an important corrective to the view of MadNCL as a general-purpose speed leader: in the general NLP setting, its principal identity is as a robust method for hard constrained problems, especially when degeneracy or LICQ failure is present.

The 2025 GPU paper broadens the empirical domain to large CUTEst problems, optimal power flow, COPS benchmarks, and degenerate SCOPF/MPCC cases. On the large OPF instance 78484epigrids, MadNCL with $r\equiv -c(x)$0 and cuDSS solves the problem in about 54 seconds, roughly an 18x speedup relative to a CPU implementation, and about 17x faster in the linear solver alone (Montoison et al., 7 Oct 2025). On the COPS set, both MadNCL-$r\equiv -c(x)$1 and MadNCL-$r\equiv -c(x)$2 perform well; on OPF problems with very sparse Jacobians, MadNCL-$r\equiv -c(x)$3 is usually reliable and faster than CPU baselines, but still generally slower than the fastest GPU IPM variants such as MadNLP with LiftedKKT (Montoison et al., 7 Oct 2025). The paper therefore treats MadNCL as extending GPU optimization to difficult degenerate problems rather than supplanting every GPU IPM on every workload.

6. Specialization to AC-SCOPF and MPCCs

The AC-SCOPF work presents MadNCL as a solver for security-constrained AC optimal power flow in polar form, with complementarity constraints modeling corrective actions such as AGC / droop control and PV/PQ bus switching (Pacaud et al., 15 Oct 2025). The problem is written as

$r\equiv -c(x)$4

This is described as a very large nonlinear program with thousands of buses, hundreds of contingencies, and millions of variables and constraints in large cases (Pacaud et al., 15 Oct 2025).

After slack variables are introduced, the MPCC is written in vertical form as

$r\equiv -c(x)$5

It can also be written with the product constraint

$r\equiv -c(x)$6

The paper explicitly states that this problem is degenerate: the feasible set has empty relative interior, MFCQ fails, and multipliers may become unbounded (Pacaud et al., 15 Oct 2025). This is the setting in which MadNCL’s robustness claims are most direct.

At outer iteration $r\equiv -c(x)$7, NCL solves

$r\equiv -c(x)$8

with

$r\equiv -c(x)$9

Here $$\min_{x,r}\ \phi(x) + y_k^T r + \tfrac12 \rho_k \|r\|^2 \quad \text{subject to} \quad c(x)+r=0,\ \ell \le x \le u.$$0 relaxes the equality constraints and $$\min_{x,r}\ \phi(x) + y_k^T r + \tfrac12 \rho_k \|r\|^2 \quad \text{subject to} \quad c(x)+r=0,\ \ell \le x \le u.$$1 relaxes the complementarity product inequality (Pacaud et al., 15 Oct 2025). If the original problem is infeasible, the associated feasibility problem minimizes $$\min_{x,r}\ \phi(x) + y_k^T r + \tfrac12 \rho_k \|r\|^2 \quad \text{subject to} \quad c(x)+r=0,\ \ell \le x \le u.$$2, which the paper identifies as particularly useful for contingency screening (Pacaud et al., 15 Oct 2025).

The linear algebra reflects the same augmented-Lagrangian regularization principle seen in the general GPU paper. The barrier-perturbed Newton system is written as

$$\min_{x,r}\ \phi(x) + y_k^T r + \tfrac12 \rho_k \|r\|^2 \quad \text{subject to} \quad c(x)+r=0,\ \ell \le x \le u.$$3

and the paper stresses that the $$\min_{x,r}\ \phi(x) + y_k^T r + \tfrac12 \rho_k \|r\|^2 \quad \text{subject to} \quad c(x)+r=0,\ \ell \le x \le u.$$4 block $$\min_{x,r}\ \phi(x) + y_k^T r + \tfrac12 \rho_k \|r\|^2 \quad \text{subject to} \quad c(x)+r=0,\ \ell \le x \le u.$$5 is non-zero because of the AL regularization. This stabilizes the Newton matrix when rows of $$\min_{x,r}\ \phi(x) + y_k^T r + \tfrac12 \rho_k \|r\|^2 \quad \text{subject to} \quad c(x)+r=0,\ \ell \le x \le u.$$6 lose rank, enabling a pivoting-free factorization, specifically a signed Cholesky factorization, compatible with NVIDIA cuDSS (Pacaud et al., 15 Oct 2025).

The implementation uses Julia 1.11, MadNLP for the inner IPM solves, ExaModels for derivative evaluation, HSL MA57 on CPU, NVIDIA cuDSS on GPU, a filter line-search globalization, Fiacco-McCormick barrier updates, tol = 1e-6, and $$\min_{x,r}\ \phi(x) + y_k^T r + \tfrac12 \rho_k \|r\|^2 \quad \text{subject to} \quad c(x)+r=0,\ \ell \le x \le u.$$7 (Pacaud et al., 15 Oct 2025). Warm starts are also used across outer NCL iterations (Pacaud et al., 15 Oct 2025).

The reported performance is substantial. For one-contingency feasibility screening, average time per contingency is 0.01 s versus 0.5 s on 118ieee, 0.2 s versus 5.5 s on 300ieee, 0.3 s versus 5.4 s on ACTIVSg500, 5.0 s versus 75.4 s on 1354pegase, and 14.1 s versus 238.4 s on 2869pegase, comparing MadNCL to Knitro (Pacaud et al., 15 Oct 2025). On ACTIVSg500 with $$\min_{x,r}\ \phi(x) + y_k^T r + \tfrac12 \rho_k \|r\|^2 \quad \text{subject to} \quad c(x)+r=0,\ \ell \le x \le u.$$8, yielding nearly 1 million variables, the paper reports 11,136 s for Knitro, 1,120 s for MadNCL-CPU, and 171 s for MadNCL-GPU; it states that MadNCL solves the case fully on the GPU in less than 3 minutes, whereas Knitro takes more than 3 hours to find an equivalent solution (Pacaud et al., 15 Oct 2025).

The limitations are also explicit. As problem size increases, degeneracy in the Newton system becomes more severe; the paper notes ill-conditioning, especially from the PV/PQ reformulation when $$\min_{x,r}\ \phi(x) + y_k^T r + \tfrac12 \rho_k \|r\|^2 \quad \text{subject to} \quad c(x)+r=0,\ \ell \le x \le u.$$9 and $[J(x)\ \ I],$0 are close; MadNCL may require many iterations on larger instances and does not converge consistently on the largest problems (Pacaud et al., 15 Oct 2025). This reinforces the broader interpretation of MadNCL: it is a robustness-oriented augmented-Lagrangian framework whose strongest results appear on infeasible, complementarity-constrained, and otherwise degenerate nonlinear programs, but whose performance remains problem dependent.