C++/DLProblem/main.cpp

This example shows how to load optimization problems from an external library using dlopen.

This example shows how to load optimization problems from an external library using dlopen.

It solves a simple quadratic program of the form:

\[ \begin{aligned} & \underset{x}{\text{minimize}} && \tfrac12 \tp x Q x \\ & \text{subject to} && Ax \le b \\ \end{aligned} \]

Problem definition in C

#include <alpaqa/dl/dl-problem.h>
#include <problem-c/export.h>
 
#include <stdlib.h>
#include <string.h>
 
#include "matmul.h"
 
typedef alpaqa_real_t real_t;
typedef alpaqa_length_t length_t;
typedef alpaqa_index_t index_t;
 
struct ProblemData {
    real_t *Q;
    real_t *A;
    real_t *work;
 
    alpaqa_problem_functions_t functions;
};
 
static void eval_objective_gradient(void *instance, const real_t *x,
                                    real_t *grad_f) {
    struct ProblemData *problem = instance;
    length_t n                  = problem->functions.num_variables;
    matmul(n, n, 1, problem->Q, x, grad_f); // grad_f = Q x
}
 
static real_t eval_objective_and_gradient(void *instance, const real_t *x,
                                          real_t *grad_f) {
    struct ProblemData *problem = instance;
    length_t n                  = problem->functions.num_variables;
    real_t result;
    eval_objective_gradient(instance, x, grad_f);
    matmul(1, n, 1, x, grad_f, &result); // result = xᵀ grad_f
    return (real_t)0.5 * result;
}
 
static real_t eval_objective(void *instance, const real_t *x) {
    struct ProblemData *problem = instance;
    return eval_objective_and_gradient(instance, x, problem->work);
}
 
static void eval_constraints(void *instance, const real_t *x, real_t *gx) {
    struct ProblemData *problem = instance;
    length_t n                  = problem->functions.num_variables;
    length_t m                  = problem->functions.num_constraints;
    matmul(m, n, 1, problem->A, x, gx); // gx = A x
}
 
static void eval_constraints_gradient_product(void *instance, const real_t *x,
                                              const real_t *y,
                                              real_t *grad_gxy) {
    (void)x;
    struct ProblemData *problem = instance;
    length_t n                  = problem->functions.num_variables;
    length_t m                  = problem->functions.num_constraints;
    matvec_transp(m, n, problem->A, y, grad_gxy); // grad_gxy = Aᵀ y
}
 
static void eval_constraints_jacobian(void *instance, const real_t *x,
                                      real_t *J_values) {
    (void)x;
    struct ProblemData *problem = instance;
    size_t n                    = (size_t)problem->functions.num_variables;
    size_t m                    = (size_t)problem->functions.num_constraints;
    memcpy(J_values, problem->A, sizeof(real_t) * m * n);
}
 
static void initialize_general_bounds(void *instance, real_t *lb, real_t *ub) {
    (void)instance;
    (void)lb; // -inf by default
    ub[0] = -1;
}
 
static struct ProblemData *create_problem(alpaqa_register_arg_t user_data) {
    (void)user_data;
    struct ProblemData *problem = malloc(sizeof(*problem));
    size_t n = 2, m = 1;
    problem->functions = (alpaqa_problem_functions_t){
        .num_variables                     = (length_t)n,
        .num_constraints                   = (length_t)m,
        .name                              = "example problem",
        .eval_objective                    = &eval_objective,
        .eval_objective_gradient           = &eval_objective_gradient,
        .eval_constraints                  = &eval_constraints,
        .eval_constraints_gradient_product = &eval_constraints_gradient_product,
        .eval_constraints_jacobian         = &eval_constraints_jacobian,
        .eval_objective_and_gradient       = &eval_objective_and_gradient,
        .initialize_general_bounds         = &initialize_general_bounds,
    };
    problem->Q    = malloc(sizeof(real_t) * n * n);
    problem->A    = malloc(sizeof(real_t) * m * n);
    problem->work = malloc(sizeof(real_t) * n);
    // Note: column major order
    problem->Q[0] = 3;
    problem->Q[1] = -1;
    problem->Q[2] = -1;
    problem->Q[3] = 3;
    problem->A[0] = 2;
    problem->A[1] = 1;
    return problem;
}
 
static void cleanup_problem(void *instance) {
    struct ProblemData *problem = instance;
    free(problem->Q);
    free(problem->A);
    free(problem->work);
    free(problem);
}
 
PROBLEM_C_EXPORT alpaqa_problem_register_t
register_alpaqa_problem(alpaqa_register_arg_t user_data) {
    struct ProblemData *problem = create_problem(user_data);
    alpaqa_problem_register_t result;
    ALPAQA_PROBLEM_REGISTER_INIT(&result);
    result.instance  = problem;
    result.cleanup   = &cleanup_problem;
    result.functions = &problem->functions;
    return result;
}
 
PROBLEM_C_EXPORT alpaqa_dl_abi_version_t register_alpaqa_problem_version(void) {
    return ALPAQA_DL_ABI_VERSION;
}

Problem solution using alpaqa

#include <alpaqa/dl/dl-problem.hpp>
#include <alpaqa/example-util.hpp>
#include <alpaqa/problem/problem-with-counters.hpp>
#include <alpaqa/structured-panoc-alm.hpp>
 
#include <filesystem>
#include <iostream>
 
namespace fs = std::filesystem;
 
// Double precision, same as in C
USING_ALPAQA_CONFIG(alpaqa::DefaultConfig);
 
int main(int argc, char *argv[]) {
    alpaqa::init_stdout();
 
    // Find the problem to load
    fs::path so_name = DLPROBLEM_DLL;
    if (argc > 1)
        so_name = fs::path(argv[1]);
    else if (argc > 0)
        so_name = fs::canonical(fs::path(argv[0])).parent_path() / so_name;
    std::cout << "Loading " << so_name << std::endl;
 
    // Instantiate a problem
    alpaqa::dl::DLProblem problem{so_name.string()};
    // Wrap the problem to count the function evaluations
    auto counted_problem = alpaqa::problem_with_counters_ref(problem);
 
    // Specify the bounds (optional, overrides any bounds set in problem.c)
    length_t n                    = problem.get_num_variables();
    length_t m                    = problem.get_num_constraints();
    vec b                         = vec::Constant(m, -1);
    const auto inf                = alpaqa::inf<config_t>;
    problem.variable_bounds.lower = vec::Constant(n, -inf);
    problem.variable_bounds.upper = vec::Constant(n, +inf);
    problem.general_bounds.lower  = vec::Constant(m, -inf);
    problem.general_bounds.upper  = b;
 
    // Define the solvers to use
    using Direction   = alpaqa::StructuredLBFGSDirection<config_t>;
    using InnerSolver = alpaqa::PANOCSolver<Direction>;
    using OuterSolver = alpaqa::ALMSolver<InnerSolver>;
 
    // Settings for the outer augmented Lagrangian method
    OuterSolver::Params almparam;
    almparam.tolerance             = 1e-8; // tolerance
    almparam.dual_tolerance        = 1e-8;
    almparam.penalty_update_factor = 10; // penalty update factor
    almparam.max_iter              = 20;
    almparam.print_interval        = 1;
 
    // Settings for the inner PANOC solver
    InnerSolver::Params panocparam;
    panocparam.max_iter       = 500;
    panocparam.print_interval = 10;
    // Settings for the L-BFGS algorithm used by PANOC
    Direction::AcceleratorParams lbfgsparam;
    lbfgsparam.memory = 2;
 
    // Create an ALM solver using PANOC as inner solver
    OuterSolver solver{
        almparam,                 // params for outer solver
        {panocparam, lbfgsparam}, // inner solver
    };
 
    // Initial guess
    vec x(2);
    x << 2, 2; // decision variables
    vec y(1);
    y << 1; // Lagrange multipliers
 
    // Solve the problem
    auto stats = solver(counted_problem, x, y);
    // y and x have been overwritten by the solution
 
    // Print the results
    std::cout << '\n' << *counted_problem.evaluations << '\n';
    std::cout //
        << "status: " << stats.status << '\n'
        << "solver: " << solver.get_name() << '\n'
        << "f = " << problem.eval_objective(x) << '\n'
        << "inner iterations: " << stats.inner.iterations << '\n'
        << "outer iterations: " << stats.outer_iterations << '\n'
        << "ε = " << stats.ε << '\n'
        << "δ = " << stats.δ << '\n'
        << "elapsed time:     "
        << std::chrono::duration<double>{stats.elapsed_time}.count() << " s"
        << '\n'
        << "x = " << x.transpose() << '\n'
        << "y = " << y.transpose() << '\n'
        << "avg τ = " << (stats.inner.sum_τ / stats.inner.count_τ) << '\n'
        << "L-BFGS rejected = " << stats.inner.direction_update_rejected << '\n'
        << "L-BFGS failures = " << stats.inner.direction_failures << '\n'
        << "Line search failures = " << stats.inner.linesearch_failures << '\n'
        << std::endl;
}