C++/FortranProblem/main.cpp

This example shows how to define optimization problems using Fortran routines.

This example shows how to define optimization problems using Fortran routines.

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 Fortran

module problem
 
   use, intrinsic::iso_c_binding, only: c_double, c_ptrdiff_t
   implicit none (type,external)
 
   integer(c_ptrdiff_t), parameter :: n = 2
   integer(c_ptrdiff_t), parameter :: m = 1
   real(c_double) :: Q(n, n)
   data q /3, -1, -1, 3/
   real(c_double) :: A(m, n)
   data a /2, 1/
 
contains
 
   ! Evaluate the cost, f(x)
   pure real(c_double) function problem_eval_objective(x) bind(C)
      real(c_double), intent(in) :: x(n)
 
      problem_eval_objective = 0.5_c_double * dot_product(x, matmul(q, x))
   end function
 
   ! Evaluate the gradient of the cost, ∇f(x)
   pure subroutine problem_eval_objective_gradient(x, grad_fx) bind(C)
      real(c_double), intent(in) :: x(n)
      real(c_double), intent(out) :: grad_fx(n)
 
      grad_fx = matmul(q, x)
   end subroutine
 
   ! Evaluate the constraints, g(x)
   pure subroutine problem_eval_constraints(x, gx) bind(C)
      real(c_double), intent(in) :: x(n)
      real(c_double), intent(out) :: gx(m)
 
      gx = matmul(a, x)
   end subroutine
 
   ! Evaluate the matrix-vector product of the gradient of the constraints and
   ! the vector y, ∇g(x) y
   pure subroutine problem_eval_constraints_gradient_product(x, y, grad_gxy) bind(C)
      real(c_double), intent(in) :: x(n)
      real(c_double), intent(in) :: y(m)
      real(c_double), intent(out) :: grad_gxy(n)
 
      grad_gxy = matmul(transpose(a), y)
   end subroutine
 
   ! Get the number of variables of the problem
   pure integer(c_ptrdiff_t) function problem_get_num_vars() bind(C)
      problem_get_num_vars = n
   end function
 
   ! Get the number of constraints of the problem
   pure integer(c_ptrdiff_t) function problem_get_num_constr() bind(C)
      problem_get_num_constr = m
   end function
 
end module

Problem solution using alpaqa

#include <alpaqa/example-util.hpp>
#include <alpaqa/problem/box-constr-problem.hpp>
#include <alpaqa/problem/problem-with-counters.hpp>
#include <alpaqa/structured-panoc-alm.hpp>
 
#include <iostream>
 
// Double precision, same as in Fortran
USING_ALPAQA_CONFIG(alpaqa::EigenConfigd);
 
// External Fortran routines
extern "C" {
ptrdiff_t problem_get_num_vars(void);
ptrdiff_t problem_get_num_constr(void);
double problem_eval_objective(const double *);
void problem_eval_objective_gradient(const double *, double *);
void problem_eval_constraints(const double *, double *);
void problem_eval_constraints_gradient_product(const double *, const double *,
                                               double *);
}
 
// Problem specification by inheriting from alpaqa::Problem
struct FortranProblem : alpaqa::BoxConstrProblem<config_t> {
    using alpaqa::BoxConstrProblem<config_t>::BoxConstrProblem;
 
    real_t eval_objective(crvec x) const {
        return problem_eval_objective(x.data());
    }
    void eval_objective_gradient(crvec x, rvec fx) const {
        problem_eval_objective_gradient(x.data(), fx.data());
    }
    void eval_constraints(crvec x, rvec gx) const {
        problem_eval_constraints(x.data(), gx.data());
    }
    void eval_constraints_gradient_product(crvec x, crvec y,
                                           rvec grad_gxy) const {
        problem_eval_constraints_gradient_product(x.data(), y.data(),
                                                  grad_gxy.data());
    }
};
 
int main() {
    alpaqa::init_stdout();
 
    // Instantiate a problem
    FortranProblem problem{problem_get_num_vars(), problem_get_num_constr()};
    // Wrap the problem to count the function evaluations
    auto counted_problem = alpaqa::problem_with_counters_ref(problem);
 
    // Specify the bounds
    const auto inf = alpaqa::inf<config_t>;
    const auto n = problem.num_variables, m = problem.num_constraints;
    vec b                         = vec::Constant(m, -1);
    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;
}