Skip to content

Commit

Permalink
Tested ResourcesFunction with MasterSolver and Solver.
Browse files Browse the repository at this point in the history
  • Loading branch information
allanleal committed Jul 30, 2021
1 parent 29c79fd commit 9af1209
Show file tree
Hide file tree
Showing 7 changed files with 46 additions and 11 deletions.
2 changes: 2 additions & 0 deletions Optima/Problem.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#include <Optima/Index.hpp>
#include <Optima/Matrix.hpp>
#include <Optima/ObjectiveFunction.hpp>
#include <Optima/ResourcesFunction.hpp>

namespace Optima {

Expand All @@ -34,6 +35,7 @@ class Problem
{
public:
Dims const dims; ///< The dimensions of the variables and constraints in the optimization problem.
ResourcesFunction r; ///< The optional function that precomputes shared resources for objective and constraint functions.
ObjectiveFunction f; ///< The objective function \eq{f(x, p)} of the optimization problem.
ConstraintFunction he; ///< The nonlinear equality constraint function \eq{h_{\mathrm{e}}(x, p)=0}.
ConstraintFunction hg; ///< The nonlinear inequality constraint function \eq{h_{\mathrm{g}}(x, p)\geq0}.
Expand Down
3 changes: 3 additions & 0 deletions Optima/Solver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,9 @@ struct Solver::Impl
// Initialize the dimensions of the master optimization problem
mproblem.dims = MasterDims(nxbar, np, ny, nz);

// Initialize the resources function in the master optimization problem
mproblem.r = problem.r;

// Create the objective function for the master optimization problem
mproblem.f = [&](ObjectiveResultRef resbar, VectorView xbar, VectorView p, VectorView c, ObjectiveOptions opts)
{
Expand Down
1 change: 1 addition & 0 deletions python/bindings/MasterProblem.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ void exportMasterProblem(py::module& m)
py::class_<MasterProblem>(m, "MasterProblem")
.def(py::init<>())
.def_readwrite("dims" , &MasterProblem::dims)
.def_readwrite("r" , &MasterProblem::r)
.def_readwrite("f" , &MasterProblem::f)
.def_readwrite("h" , &MasterProblem::h)
.def_readwrite("v" , &MasterProblem::v)
Expand Down
2 changes: 2 additions & 0 deletions python/bindings/Optima.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ void exportOptions(py::module& m);
void exportProblem(py::module& m);
void exportResidualFunction(py::module& m);
void exportResidualVector(py::module& m);
void exportResourcesFunction(py::module& m);
void exportResult(py::module& m);
void exportSensitivity(py::module& m);
void exportSensitivitySolver(py::module& m);
Expand Down Expand Up @@ -105,6 +106,7 @@ PYBIND11_MODULE(optima4py, m)
exportProblem(m);
exportResidualFunction(m);
exportResidualVector(m);
exportResourcesFunction(m);
exportResult(m);
exportSensitivity(m);
exportSensitivitySolver(m);
Expand Down
3 changes: 3 additions & 0 deletions python/bindings/Problem.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ using namespace Optima;

void exportProblem(py::module& m)
{
auto set_r = [](Problem& self, ResourcesFunction r) { self.r = r; };
auto set_f = [](Problem& self, ObjectiveFunction f) { self.f = f; };
auto set_he = [](Problem& self, ConstraintFunction he) { self.he = he; };
auto set_hg = [](Problem& self, ConstraintFunction hg) { self.hg = hg; };
Expand All @@ -38,6 +39,7 @@ void exportProblem(py::module& m)
auto set_bec = [](Problem& self, MatrixView4py bec) { self.bec = bec; };
auto set_bgc = [](Problem& self, MatrixView4py bgc) { self.bgc = bgc; };

auto get_r = [](Problem& self) { return self.r; };
auto get_f = [](Problem& self) { return self.f; };
auto get_he = [](Problem& self) { return self.he; };
auto get_hg = [](Problem& self) { return self.hg; };
Expand All @@ -52,6 +54,7 @@ void exportProblem(py::module& m)
py::class_<Problem>(m, "Problem")
.def(py::init<const Dims&>())
.def_readonly("dims" , &Problem::dims)
.def_property("r" , get_r, set_r)
.def_property("f" , get_f, set_f)
.def_property("he" , get_he, set_he)
.def_property("hg" , get_hg, set_hg)
Expand Down
19 changes: 15 additions & 4 deletions tests/MasterSolver.py
Original file line number Diff line number Diff line change
Expand Up @@ -74,10 +74,18 @@ def testMasterSolver(nx, np, ny, nz, nl, nul, nuu, diagHxx):
cp = ones(np)
cz = ones(nz)

class Resources:
dx, dp = None, None
def __iter__(self): return iter((self.dx, self.dp))

resources = Resources()

def resourcesfn_r(x, p, c, fopts, hopts, vopts):
resources.dx = x - cx
resources.dp = p - cp

def objectivefn_f(res, x, p, c, opts):
dx = x - cx
dp = p - cp
dx, dp = resources
res.f = 0.5 * dx.T @ Hxx @ dx + dx.T @ Hxp @ dp
res.fx = Hxx @ dx + Hxp @ dp
res.fxx = Hxx
Expand All @@ -87,14 +95,16 @@ def objectivefn_f(res, x, p, c, opts):
res.succeeded = True

def constraintfn_h(res, x, p, c, opts):
res.val = Jx @ (x - cx) + Jp @ (p - cp)
dx, dp = resources
res.val = Jx @ dx + Jp @ dp
res.ddx = Jx
res.ddp = Jp
res.ddx4basicvars = False
res.succeeded = True

def constraintfn_v(res, x, p, c, opts):
res.val = Vpx @ (x - cx) + Vpp @ (p - cp)
dx, dp = resources
res.val = Vpx @ dx + Vpp @ dp
res.ddx = Vpx
res.ddp = Vpp
res.ddx4basicvars = False
Expand All @@ -113,6 +123,7 @@ def constraintfn_v(res, x, p, c, opts):

problem = MasterProblem()
problem.dims = dims
problem.r = resourcesfn_r
problem.f = objectivefn_f
problem.h = constraintfn_h
problem.v = constraintfn_v
Expand Down
27 changes: 20 additions & 7 deletions tests/Solver.py
Original file line number Diff line number Diff line change
Expand Up @@ -121,26 +121,38 @@ def testMasterSolver(nx, np, ny, nz, nl, nul, nuu, diagHxx):

cx[ju] = +1.0e4 # large positive number to ensure x variables with ju indices are indeed unstable!

def objectivefn_f(res, x, p, c, opts):
class Resources:
f, fx, v, h = None, None, None, None

resources = Resources()

def resourcesfn_r(x, p, c, fopts, hopts, vopts):
# Precompute here f, fx, v and h, to demonstrate shared resources among functions are correctly precalculated
cx = c[:nx]
res.f = 0.5 * (x.T @ Hxx @ x) + x.T @ Hxp @ p + cx.T @ x
res.fx = Hxx @ x + Hxp @ p + cx
cp = c[nx:][:np]
cz = c[nx:][np:][ny:]
resources.f = 0.5 * (x.T @ Hxx @ x) + x.T @ Hxp @ p + cx.T @ x
resources.fx = Hxx @ x + Hxp @ p + cx
resources.v = Vpx @ x + Vpp @ p + cp
resources.h = Jx @ x + Jp @ p + cz

def objectivefn_f(res, x, p, c, opts):
res.f = resources.f # already computed in the resources function
res.fx = resources.fx # already computed in the resources function
res.fxx = Hxx
res.fxp = Hxp
if opts.eval.fxc: # compute d(fx)/dc = [d(fx)/d(cx), d(fx)/d(cp), d(fx)/d(cy), d(fx)/d(cz)]
res.fxc = npy.hstack([Ixx, Oxp, Oxy, Oxz])

def constraintfn_v(res, x, p, c, opts):
cp = c[nx:][:np]
res.val = Vpx @ x + Vpp @ p + cp
res.val = resources.v # already computed in the resources function
res.ddx = Vpx
res.ddp = Vpp
if opts.eval.ddc: # compute dv/dc = [dv/d(cx), dv/d(cp), dv/d(cy), dv/d(cz)]
res.ddc = npy.hstack([Opx, Ipp, Opy, Opz])

def constraintfn_h(res, x, p, c, opts):
cz = c[nx:][np:][ny:]
res.val = Jx @ x + Jp @ p + cz
res.val = resources.h # already computed in the resources function
res.ddx = Jx
res.ddp = Jp
if opts.eval.ddc: # compute dv/dc = [dv/d(cx), dv/d(cp), dv/d(cy), dv/d(cz)]
Expand All @@ -163,6 +175,7 @@ def constraintfn_h(res, x, p, c, opts):
dims.c = nc

problem = Problem(dims)
problem.r = resourcesfn_r
problem.f = objectivefn_f
problem.he = constraintfn_h
problem.v = constraintfn_v
Expand Down

0 comments on commit 9af1209

Please sign in to comment.