Commit 2a280d2d authored by webmanue's avatar webmanue Committed by Manuel Weberndorfer
Browse files

add option to set nonlinear equality constraints

- some solvers provided by TAO support these constraints
parent 3203fc82
// © 2020 ETH Zurich, Mechanics and Materials Lab
// © 2020, 2022 ETH Zurich, Mechanics and Materials Lab
// © 2020 California Institute of Technology
//
// Licensed under the Apache License, Version 2.0 (the "License");
......@@ -25,6 +25,7 @@
#include "ae108/cpppetsc/UniqueEntity.h"
#include "ae108/cpppetsc/Vector.h"
#include <functional>
#include <optional>
#include <petscmat.h>
#include <petscmath.h>
#include <petscsys.h>
......@@ -76,7 +77,8 @@ public:
ssfls,
asils,
asfls,
ipm
ipm,
pdipm,
};
/**
......@@ -132,6 +134,29 @@ public:
void setBounds(const distributed<vector_type> &lowerBounds,
const distributed<vector_type> &upperBounds);
struct EqualityConstraints {
/**
* @brief The solver attempts to reach a residual of zero.
*/
GradientFunctor residual;
/**
* @brief The jacobian of the residual.
*/
HessianFunctor jacobian;
/**
* @brief Used to store residuals and jacobians internally.
*/
std::pair<distributed<vector_type>, matrix_type> buffer;
};
/**
* @brief Set the equality constraints of the problem.
* @note Not all solvers consider these constraints.
*/
void setConstraints(EqualityConstraints constraints);
/**
* @brief Minimizes the objective function using optimization.
*
......@@ -164,6 +189,8 @@ private:
UniqueEntity<Tao> _tao;
matrix_type _buffer_matrix;
std::optional<EqualityConstraints> _constraints;
};
template <class Policy>
......@@ -288,6 +315,10 @@ template <class Policy> void TAOSolver<Policy>::setType(const Type type) {
taoType = TAOPOUNDERS;
break;
}
case Type::pdipm: {
taoType = TAOPDIPM;
break;
}
case Type::lcl: {
taoType = TAOLCL;
break;
......@@ -323,6 +354,19 @@ void TAOSolver<Policy>::setBounds(const distributed<vector_type> &lowerBounds,
_tao.get(), lowerBounds.unwrap().data(), upperBounds.unwrap().data()));
}
template <class Policy>
void TAOSolver<Policy>::setConstraints(EqualityConstraints constraints) {
_constraints.emplace(std::move(constraints));
Policy::handleError(TaoSetEqualityConstraintsRoutine(
_tao.get(), _constraints->buffer.first.unwrap().data(),
&TAOSolver::gradientAdapter, &_constraints->residual));
Policy::handleError(TaoSetJacobianEqualityRoutine(
_tao.get(), _constraints->buffer.second.data(),
_constraints->buffer.second.data(), &TAOSolver::hessianAdapter,
&_constraints->jacobian));
}
template <class Policy> void TAOSolver<Policy>::checkConvergence() const {
const auto errorCode = [this]() {
auto reason = TaoConvergedReason{};
......
// © 2020 ETH Zurich, Mechanics and Materials Lab
// © 2020, 2022 ETH Zurich, Mechanics and Materials Lab
// © 2020 California Institute of Technology
//
// Licensed under the Apache License, Version 2.0 (the "License");
......@@ -140,6 +140,54 @@ TYPED_TEST(TAOSolver_Test, raises_exception_on_nonconvergence) {
EXPECT_THROW(callSolve(), TAOSolverDivergedException);
}
TYPED_TEST(TAOSolver_Test, minimizes_with_equality_constraints) {
using real_type = typename TestFixture::real_type;
using vector_type = typename TestFixture::vector_type;
using matrix_type = typename TestFixture::matrix_type;
typename TestFixture::solver_type solver(
matrix_type(2, 2), TestFixture::solver_type::Type::pdipm);
const auto value = 5.;
solver.setConstraints({
[&](const distributed<vector_type> &input,
distributed<vector_type> *const output) {
const auto full = vector_type::fromDistributed(input);
output->unwrap().setZero();
output->unwrap().replace().element(0, full(0) - value);
},
[](const distributed<vector_type> &, matrix_type *const output) {
output->setZero();
const auto replacer =
output->assemblyView().replace().element(0, 0, 1.);
},
{tag<DistributedTag>(vector_type(1)), matrix_type(1, 2)},
});
const auto solution = solver.solve(
[](const distributed<vector_type> &input, real_type *const output) {
const auto full = vector_type::fromDistributed(input);
*output = std::pow(full(0) - 1., 2.) + std::pow(full(1) - 2., 2.);
},
[](const distributed<vector_type> &input,
distributed<vector_type> *const output) {
const auto full = vector_type::fromDistributed(input);
const auto replacer = output->unwrap().replace();
replacer(0) = 2. * (full(0) - 1.);
replacer(1) = 2. * (full(1) - 2.);
},
[](const distributed<vector_type> &, matrix_type *const output) {
const auto replacer = output->assemblyView().replace();
replacer(0, 0) = 2.;
replacer(1, 1) = 2.;
},
tag<DistributedTag>(vector_type::fromList({3., 4.})));
EXPECT_THAT(solution.unwrap(), ScalarEqIfLocal(0, value));
EXPECT_THAT(solution.unwrap(), ScalarEqIfLocal(1, 2.));
}
template <class Policy> struct TAOSolver_BoundsTest : Test {
using solver_type = TAOSolver<Policy>;
using value_type = typename solver_type::value_type;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment