diff --git a/include/spii/constrained_function.h b/include/spii/constrained_function.h new file mode 100644 index 0000000..82d81e1 --- /dev/null +++ b/include/spii/constrained_function.h @@ -0,0 +1,102 @@ +// Petter Strandmark 2013. +#ifndef SPII_CONSTRAINED_FUNCTION_H +#define SPII_CONSTRAINED_FUNCTION_H + +#include + +namespace spii { + +class Solver; +struct SolverResults; + +// +// Builds an augmented lagrangian of the terms added +// to the objective function and constraints. +// +// See Nocedal & Wright, chapter 17. +// +class SPII_API ConstrainedFunction +{ +public: + ConstrainedFunction(); + ~ConstrainedFunction(); + + // Function improvement tolerance. The constrained solver stops + // if |df| / (|f| + tol) < tol + // between outer iterations. + double function_improvement_tolerance = 0; + + // When updating the dual variables, the constrained solver + // stop if the maximum norm of the change is less than this + // tolerance. + double dual_change_tolerance = 1e-8; + + // The maximum number of (outer) iterations the solver will + // perform. + int max_number_of_iterations = 100; + + // Access to the objective function for e.g. evaluation. + const Function& objective() const; + + // Whether the provided point is feasible. + bool is_feasible() const; + + // + // Adds a term to the objective function. + // + void add_term(std::shared_ptr term, const std::vector& arguments); + + template + void add_term(std::shared_ptr term, PointerToDouble... args) + { + add_term(term, {args...}); + } + + template + void add_term(PointerToDouble... args) + { + add_term(std::make_shared(), {args...}); + } + + // + // Adds a term as a constraint. + // + // This will add the constraint c(x) ≤ 0. + // + void add_constraint_term(const std::string& constraint_name, + std::shared_ptr term, + const std::vector& arguments); + + template + void add_constraint_term(const std::string& constraint_name, + std::shared_ptr term, + PointerToDouble... args) + { + add_constraint_term(constraint_name, term, {args...}); + } + + template + void add_constraint_term(const std::string& constraint_name, + PointerToDouble... args) + { + add_constraint_term(constraint_name, std::make_shared(), {args...}); + } + + + // Minimized the constrained function. In the future, this method + // may be extraced to e.g. a ConstrainedSolver class. + void solve(const Solver& solver, SolverResults* results); + +private: + // Disable copying for now. + ConstrainedFunction(const ConstrainedFunction&); + + class Implementation; + // unique_pointer would have been nice, but there are issues + // with sharing these objects across DLL boundaries in VC++. + Implementation* impl; +}; + +} // namespace spii. + +#endif diff --git a/source/constrained_function.cpp b/source/constrained_function.cpp new file mode 100644 index 0000000..3367ad3 --- /dev/null +++ b/source/constrained_function.cpp @@ -0,0 +1,282 @@ +// Petter Strandmark 2013. +// +// Constrained minimization. +// +// See Nocedal & Wright, chapter 17. + +#include + +#include +#include + +using namespace std; + +namespace spii +{ + +class Constraint + : public Function +{ +public: + double lambda = 0; + double cached_value = 0; +}; + +// Nocedal & Wright p. 524 +// Equation 17.65. +// +// But this implementation uses c(x) ≤ 0 instead +// of c(x) ≥ 0 as in Nocedal & Wright. +class Phi + : public Term +{ +public: + Phi(std::shared_ptr term_, double* sigma_, double* mu_) + : term(term_), sigma(sigma_), mu(mu_) + { + } + + int number_of_variables() const override + { + return term->number_of_variables(); + } + + int variable_dimension(int var) const override + { + return term->variable_dimension(var); + } + + double evaluate(double * const * const variables) const override + { + double t = term->evaluate(variables); + if (- t - *sigma / *mu <= 0) { + return *sigma * t + (*mu / 2) * t*t; + } + else { + return - 1.0 / (2.0 * *mu) * (*sigma)*(*sigma); + } + } + + double evaluate(double * const * const variables, + std::vector* gradient) const override + { + double t = term->evaluate(variables, gradient); + + if (- t - *sigma / *mu <= 0) { + for (size_t i = 0; i < number_of_variables(); ++i) { + for (size_t j = 0; j < variable_dimension(i); ++j) { + (*gradient)[i][j] = *sigma * (*gradient)[i][j] + + (*mu) * t * (*gradient)[i][j]; + } + } + return *sigma * t + (*mu / 2) * t*t; + } + else { + for (size_t i = 0; i < number_of_variables(); ++i) { + for (size_t j = 0; j < variable_dimension(i); ++j) { + (*gradient)[i][j] = 0; + } + } + return - 1.0 / (2.0 * *mu) * (*sigma)*(*sigma); + } + } + + double evaluate(double * const * const variables, + std::vector* gradient, + std::vector< std::vector >* hessian) const override + { + check(false, "Phi: hessian not supported."); + return 0; + } + +private: + const shared_ptr term; + double* const sigma; + double* const mu; +}; + +class ConstrainedFunction::Implementation +{ +public: + double mu; + + Function augmented_lagrangian; + Function objective; + map constraints; +}; + +ConstrainedFunction::ConstrainedFunction() + : impl{new Implementation} +{ +} + +ConstrainedFunction::~ConstrainedFunction() +{ + delete impl; +} + +void ConstrainedFunction::add_term(shared_ptr term, + const vector& arguments) +{ + impl->objective.add_term(term, arguments); + impl->augmented_lagrangian.add_term(term, arguments); +} + +void ConstrainedFunction::add_constraint_term(const string& constraint_name, + shared_ptr term, + const vector& arguments) +{ + check(impl->constraints.find(constraint_name) == impl->constraints.end(), + "add_constraint_term: Term already added."); + auto& constraint = impl->constraints[constraint_name]; + + constraint.add_term(term, arguments); + + auto phi = make_shared(term, &constraint.lambda, &impl->mu); + impl->augmented_lagrangian.add_term(move(phi), arguments); +} + +const Function& ConstrainedFunction::objective() const +{ + return impl->objective; +} + +bool ConstrainedFunction::is_feasible() const +{ + for (auto& itr: impl->constraints) { + auto c_x = itr.second.evaluate(); + if (c_x > 1e-12) { + return false; + } + } + return true; +} + +void ConstrainedFunction::solve(const Solver& solver, SolverResults* results) +{ + results->exit_condition = SolverResults::NA; + + auto& mu = impl->mu; + mu = 10; + double nu = 1.0 / std::pow(mu, 0.1); + + auto log_function = solver.log_function; + + double f_prev = numeric_limits::quiet_NaN(); + + int iterations = 0; + while (true) { + + // Minimize the smooth approximation of the Lagrangian. + solver.solve(impl->augmented_lagrangian, results); + double f = impl->objective.evaluate(); + + // Update all lambdas. + double infeasibility = - numeric_limits::infinity(); + double max_violation = 0; + for (auto& itr: impl->constraints) { + auto c_x = itr.second.evaluate(); + itr.second.cached_value = c_x; + auto& lambda = itr.second.lambda; + + infeasibility = max(infeasibility, c_x * lambda); + max_violation = max(max_violation, c_x); + } + + if (log_function) { + char str[1024]; + std::sprintf(str, + " ___________________________________________________________\n" + "| mu | nu | objective | infeas. | max viol. |\n" + "+--------+--------+---------------+------------+------------+\n" + "|%7.1e |%7.1e | %+10.6e | %10.3e | %10.3e |\n" + "|________|________|_______________|____________|____________|", + mu, nu, f, infeasibility, max_violation); + log_function(str); + } + + if (std::abs(f - f_prev) / (std::abs(f) + this->function_improvement_tolerance) < + this->function_improvement_tolerance) { + results->exit_condition = SolverResults::FUNCTION_TOLERANCE; + break; + } + + if (max_violation <= nu) { + // The maximum constraint violation is small enough. + // Update the dual variables (explicit formula). + + double max_change = 0; + for (auto& itr: impl->constraints) { + auto c_x = itr.second.cached_value; + auto& lambda = itr.second.lambda; + + double prev = lambda; + if (c_x + lambda / mu <= 0) { + lambda = 0; + } + else { + lambda = lambda + mu * c_x; + } + max_change = std::max(max_change, std::abs(prev - lambda)); + } + + if (log_function) { + stringstream sout; + sout << "Updating dual variables. " + << "Maximum change: " << max_change << "."; + log_function(sout.str()); + } + + if (max_change < this->dual_change_tolerance) { + results->exit_condition = SolverResults::GRADIENT_TOLERANCE; + break; + } + + // Nocedal & Wright. + nu = nu / std::pow(mu, 0.9); + } + else { + // The maximum constraint violation too big. + // Update the penalty variable to decrease it. + + if (log_function) { + log_function("Updating penalty parameter."); + } + + // Increase penalty parameter. + mu *= 100; + + // Nocedal & Wright. + nu = 1.0 / std::pow(mu, 0.1); + } + + if (log_function) { + log_function(""); + for (auto& itr: impl->constraints) { + auto c_x = itr.second.cached_value; + auto& lambda = itr.second.lambda; + + stringstream sout_name, sout; + sout_name << "lambda[" << itr.first << "]"; + sout << left << setfill('.') << setw(25) << sout_name.str() + << ": " << setfill(' ') << setw(10) << lambda; + if (c_x > 0) { + sout << " Violation : " << c_x; + } + log_function(sout.str()); + } + log_function(""); + } + + iterations++; + if (iterations >= this->max_number_of_iterations) { + results->exit_condition = SolverResults::NO_CONVERGENCE; + break; + } + + // Go to next iteration. + f_prev = f; + } +} + +} // namespace spii. diff --git a/tests/test_constrained_function.cpp b/tests/test_constrained_function.cpp new file mode 100644 index 0000000..56f4c99 --- /dev/null +++ b/tests/test_constrained_function.cpp @@ -0,0 +1,114 @@ +// Petter Strandmark 2013. + +#define CATCH_CONFIG_MAIN +#include + +#include +#include +#include + +using namespace spii; +using namespace std; + +class ObjectiveTerm +{ +public: + template + R operator()(const R* const x) const + { + auto dx = x[0] - 2.0; + auto dy = x[1] - 2.0; + return dx*dx + dy*dy; + } +}; + +// x·x + y·y ≤ 1 +class ConstraintTerm +{ +public: + template + R operator()(const R* const x) const + { + auto dx = x[0]; + auto dy = x[1]; + return dx*dx + dy*dy - 1; + } +}; + +// x ≤ 100 +class LessThan100 +{ +public: + template + R operator()(const R* const x) const + { + auto dx = x[0]; + return dx - 100; + } +}; + +void test_simple_constrained_function(const std::vector x_start, bool feasible) +{ + ConstrainedFunction function; + stringstream log_stream; + LBFGSSolver solver; + solver.log_function = + [&log_stream](const string& str) + { + log_stream << str << endl; + }; + + auto x = x_start; + + function.add_term(make_shared>(), &x[0]); + function.add_constraint_term("circle", make_shared>(), &x[0]); + function.add_constraint_term("less than 100", make_shared>(), &x[0]); + CHECK(function.is_feasible() == feasible); + + SolverResults results; + function.solve(solver, &results); + log_stream << results << endl; + + INFO(log_stream.str()); + CAPTURE(x[0]); + CAPTURE(x[1]); + CHECK((x[0]*x[0] + x[1]*x[1]) <= (1.0 + 1e-8)); + + auto dx = x[0] - 2; + auto dy = x[1] - 2; + CHECK(abs(function.objective().evaluate() - (dx*dx + dy*dy)) <= 1e-14); +} + +TEST_CASE("feasible_start") +{ + test_simple_constrained_function({0.0, 0.0}, true); +} + + +TEST_CASE("infeasible_start") +{ + test_simple_constrained_function({10.0, 4.0}, false); +} + +TEST_CASE("max_iterations") +{ + ConstrainedFunction function; + function.max_number_of_iterations = 1; + vector x = {0, 0}; + function.add_term(make_shared>(), &x[0]); + function.add_constraint_term("circle", make_shared>(), &x[0]); + function.add_constraint_term("less than 100", make_shared>(), &x[0]); + LBFGSSolver solver; + solver.log_function = nullptr; + SolverResults results; + function.solve(solver, &results); + CHECK(results.exit_condition == SolverResults::NO_CONVERGENCE); +} + +TEST_CASE("no_readding_constraints") +{ + ConstrainedFunction function; + double x[2] = {0, 0}; + function.add_constraint_term("circle", make_shared>(), &x[0]); + CHECK_THROWS(function.add_constraint_term("circle", make_shared>(), &x[0])); +}