Optimization¶

anvil includes a Sequential Quadratic Programming (SQP) solver that, like NumericalFunction, can be both called from Python (JIT-compiled) and exported to standalone C++. The full SQP loop -- function evaluation, Hessian computation, QP solve via PIQP, primal/dual updates -- runs in native code with zero Python overhead per iteration.

Problem formulation¶

The SQP solver handles parameterized nonlinear programs of the form:

\[ \begin{align*} \min_x & \quad f(x, p) \\ \text{s.t.} & \quad g_\mathrm{eq}(x, p) = b_\mathrm{eq}, \\ &\quad g_\mathrm{lb} \leq g_\mathrm{ineq}(x, p) \leq g_\mathrm{ub}, \\ &\quad x_\mathrm{lb} \leq x \leq x_\mathrm{ub} \end{align*} \]

where:

\(x \in \mathbb{R}^n\) is the decision variable
\(p \in \mathbb{R}^{n_p}\) is a runtime parameter vector (fixed data, never differentiated)
\(f: \mathbb{R}^n \times \mathbb{R}^{n_p} \to \mathbb{R}\) is the cost function
\(g_\mathrm{eq}: \mathbb{R}^n \times \mathbb{R}^{n_p} \to \mathbb{R}^{m_\mathrm{eq}}\) are equality constraints
\(g_\mathrm{ineq}: \mathbb{R}^n \times \mathbb{R}^{n_p} \to \mathbb{R}^{m_\mathrm{ineq}}\) are inequality constraints
\(x_\mathrm{lb}, x_\mathrm{ub}\) are variable bounds

All callbacks take (x, p) as inputs. All derivatives are taken with respect to x only.

Defining a problem¶

Dense¶

import anvil as av
from anvil.opti import DenseProblem, DenseSQPFunction, SQPSettings, HessianApproximation

# Cost: (x, p) -> scalar
cost_fn = av.NumericalFunction(
    "cost",
    lambda x, p: (x[1] - 2) ** 2 - x[2] + x[3],
    (av.Arg(4), av.Arg(0)),  # Arg(0) = no runtime parameters
)

# Equality constraint: (x, p) -> R^m_eq
eq_fn = av.NumericalFunction(
    "eq_constraints",
    lambda x, p: (x[0] ** 2 + x[1] - 1).reshape(1),
    (av.Arg(4), av.Arg(0)),
)

problem = DenseProblem(
    cost_fn=cost_fn,
    eq_constraints_fn=eq_fn,
    ineq_constraints_fn=None,
    x_lb=av.Tensor([-1.0, -1e6, -1.0, -1.0]),
    x_ub=av.Tensor([1.0, 1e6, 1.0, 1.0]),
    eq_rhs=av.Tensor([0.0]),
)

When nparam == 0 (i.e. the last input is Arg(0)), the p argument is a zero-length tensor -- callbacks still accept it but can ignore it.

Sparse¶

For large problems, use SparseProblem which computes sparse constraint Jacobians via spjacobian:

from anvil.opti import SparseProblem, SparseSQPFunction

sparse_problem = SparseProblem(
    cost_fn=cost_fn,
    eq_constraints_fn=eq_fn,
    x_lb=x_lb,
    x_ub=x_ub,
    eq_rhs=av.Tensor([0.0]),
)

The interface is identical to DenseProblem. Internally, constraint Jacobians and the Lagrangian Hessian are computed as sparse matrices.

Solver settings¶

from anvil.opti import GlobalizationStrategy

settings = SQPSettings(
    max_iter=1000,
    verbose=False,
    compute_timings=False,
    eps_prim=1e-6,                    # primal feasibility tolerance
    eps_dual=1e-4,                    # dual feasibility tolerance
    hessian_approximation=HessianApproximation.EXACT,
    globalization_strategy=GlobalizationStrategy.LINE_SEARCH_L1,
    line_search_max_iter=20,
    tau=0.5,                          # line search backtracking factor
    eta=1e-4,                         # Armijo condition parameter
    rho=1.0,                          # L1 penalty parameter
)

Hessian approximation options:

EXACT: full Lagrangian Hessian (cost + constraints)
EXACT_NO_CONSTRAINTS: only cost Hessian (ignores constraint curvature)

Creating and calling a solver¶

import numpy as np

solver = DenseSQPFunction(
    name="my_solver",
    problem=problem,
    settings=settings,
)

result = solver(
    initial_x=np.zeros(4),
    initial_lam_bounds=np.zeros(4),  # optional
    initial_lam_eq=np.zeros(1),      # optional
)

print(result.status)     # SQPStatus.SOLVED
print(result.x)          # optimal solution
print(result.iter)       # number of SQP iterations
print(result.qp_iter)    # total QP iterations

With runtime parameters:

# nparam=2: inferred from Arg(2) on the last input
cost_fn = av.NumericalFunction(
    "cost",
    lambda x, p: (x[0] - p[0]) ** 2 + x[1] ** 2,
    (av.Arg(2), av.Arg(2)),
)

# When calling, p is required:
result = solver(initial_x=np.zeros(2), p=np.array([2.0, 1.0]))

Generated C++ interface¶

av.generate_module("my_module", [solver])

The generated header contains:

namespace my_solver {

enum class SQPStatus {
    SOLVED = 1, MAX_ITER_REACHED = -1, INFEASIBLE = -2,
    NON_CONVEX_QP = -3, QP_SOLVER_ERROR = -4,
    UNSOLVED = -9, INVALID_SETTINGS = -10
};

struct SQPSettings {
    int max_iter = 1000;
    bool verbose = false;
    bool compute_timings = false;
    double eps_prim = 1e-6;
    double eps_dual = 1e-4;
    // ... line search parameters
};

struct SQPTimings { /* per-component timing breakdown */ };
struct SQPResult {
    SQPStatus status;
    int iter, qp_iter;
    SQPTimings timings;
};

// Problem dimension constants
constexpr int my_solver_n = 4;
constexpr int my_solver_m_eq = 1;

// Buffer types
typedef Buffer<double, 4> X_t;
typedef Buffer<double, 4> LAM_BOUNDS_t;
typedef Buffer<double, 1> LAM_EQ_t;

// Workspace (contains PIQP data, all function workspaces, QP intermediates)
struct Ws { /* ... */ };

Ws init_ws();
void deinit_ws(Ws& ws);

SQPResult call(
    const SQPSettings& settings, Ws& ws,
    const X_t& initial_x,
    const LAM_BOUNDS_t& initial_lam_bounds,
    const LAM_EQ_t& initial_lam_eq
);

} // namespace my_solver

Usage:

#include "my_module.hpp"

int main() {
    using namespace my_solver;

    auto ws = init_ws();
    auto x0 = X_t::alloc();
    auto lam_bounds = LAM_BOUNDS_t::alloc();
    auto lam_eq = LAM_EQ_t::alloc();

    std::fill(x0.data, x0.data + X_t::size, 0.0);
    std::fill(lam_bounds.data, lam_bounds.data + LAM_BOUNDS_t::size, 0.0);
    std::fill(lam_eq.data, lam_eq.data + LAM_EQ_t::size, 0.0);

    SQPSettings settings;
    settings.max_iter = 100;
    settings.verbose = true;

    SQPResult result = call(settings, ws, x0, lam_bounds, lam_eq);

    if (result.status == SQPStatus::SOLVED)
        printf("Solved in %d iterations\n", result.iter);

    deinit_ws(ws);
    X_t::free(x0);
    LAM_BOUNDS_t::free(lam_bounds);
    LAM_EQ_t::free(lam_eq);
}

PIQP integration¶

anvil ships with the PIQP C interface. To integrate generated SQP solvers in your project:

from anvil.piqp_paths import LIB_DIR, INCLUDE_DIR, get_library_path
print(f"Library: {LIB_DIR}")
print(f"Headers: {INCLUDE_DIR}")

CMake¶

# Get paths from: python -c "from anvil.piqp_paths import *; print(LIB_DIR, INCLUDE_DIR)"
set(PIQP_LIB_DIR "/path/from/python")
set(PIQP_INCLUDE_DIR "/path/from/python")

add_executable(my_app main.cpp my_module.cpp)
target_include_directories(my_app PRIVATE ${PIQP_INCLUDE_DIR} .)
target_link_directories(my_app PRIVATE ${PIQP_LIB_DIR})
target_link_libraries(my_app PRIVATE piqpc)