Examples

The following examples illustrate the functionality of Coek.

Model Syntax

Rosenbrock

#include <coek/coek.hpp>

coek::Model rosenbr()
{
    coek::Model m;

    auto x0 = m.add(coek::variable().value(1));
    auto x1 = m.add(coek::variable().value(2));

    m.add_objective(100 * pow(x1 - pow(x0, 2), 2) + pow(x0 - 1, 2));

    return m;
}

//   Source:  problem 1 in
//   J.J. More', B.S. Garbow and K.E. Hillstrom,
//   "Testing Unconstrained Optimization Software",
//   ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981.

The function rosenbr constructs a coek::Model object and returns it. In Lines 7 and 8, two decision variables are defined, initialized and added to the model. On Line 10, an expression is constructed for the optimization objective, and this expression is added to the model.

Multidimentional Rosenbrock

#include <coek/coek.hpp>

coek::Model srosenbr_vector()
{
    coek::Model m;

    size_t N = 10000;
    std::vector<coek::Variable> x(N);
    for (size_t i = 0; i < N; i++) {
        if (i % 2 == 0)
            m.add(x[i].value(-1.2));
        else
            m.add(x[i].value(1));
    }

    auto obj = coek::expression();
    for (size_t i = 0; i < N / 2; i++)
        obj += 100 * pow(x[2 * i] - pow(x[2 * i + 1], 2), 2) + pow(x[2 * i] - 1, 2);
    m.add_objective(obj);

    return m;
}

// Source:  problem 21 in
// J.J. More', B.S. Garbow and K.E. Hillstrom,
// "Testing Unconstrained Optimization Software",
// ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981.

The function rosenbr_vector constructs a coek::Model object and returns it. In Lines 7-14, a std::vector of Coek variables is declared, initialized, and added to the model. In Lines 16-19, an expression is constructed for the optimization objective, and this expression is added to the model.

#include <coek/coek.hpp>

coek::Model srosenbr_array()
{
    coek::Model m;

    size_t N = 10000;
    auto x = m.add(coek::variable("x", N));
    for (size_t i : coek::range(N)) {
        if (i % 2 == 0)
            x(i).value(-1.2);
        else
            x(i).value(1);
    }

    auto obj = coek::expression();
    for (size_t i : coek::range(N / 2))
        obj += 100 * pow(x(2 * i) - pow(x(2 * i + 1), 2), 2) + pow(x(2 * i) - 1, 2);
    m.add_objective(obj);

    return m;
}

// Source:  problem 21 in
// J.J. More', B.S. Garbow and K.E. Hillstrom,
// "Testing Unconstrained Optimization Software",
// ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981.

The function rosenbr_array constructs a coek::Model object using the more compact array syntax supported by Coek. Line 8 declares an array of variables and adds it to the model. Lines 9-14 initialize the variables, using the coek::range function to simplify the loop. Similarly, Lines 16-19, construct an expression for the objective using the coek::range function, which is then added to the model.

Applying Optimizers

Solving a Linear Program

#include <coek/coek.hpp>

void simplelp1_solve()
{
    coek::Model m;

    auto x = m.add(coek::variable("x").bounds(0, m.inf));
    auto y = m.add(coek::variable("y").bounds(0, m.inf));

    m.add_objective(50 * x + 40 * y).sense(m.maximize);
    m.add_constraint(2 * x + 3 * y <= 1500);
    m.add_constraint(2 * x + y <= 1000);

    //
    // Optimize the model
    //
    coek::Solver solver("gurobi");
    auto status = solver.solve(m);

    std::cout << "Value of " + x.name() + ": " << x.value() << std::endl;
    std::cout << "Value of " + y.name() + ": " << y.value() << std::endl;
}

The function simplelp1_solve constructs a simple linear program that is solved using the ipopt nonlinear optimization solver.

Solving a Quadratic Problem

#include <coek/coek.hpp>
#include <vector>

void invquad_array_solve()
{
    // Create model
    coek::Model m;

    size_t N = 10;
    auto p = coek::parameter("p", N);
    for (auto& param : p) param.value(0.5);

    // Initialize variables and add them to the model
    auto x = coek::variable("x", N).bounds(-10, 10).value(0.0);
    m.add(x);

    // Create objective and add it to the model
    auto e = coek::expression();
    for (size_t i : coek::range(N)) e -= (x(i) - p(i)) * (x(i) - p(i));
    m.add_objective(e);

    // Optimize the model
    coek::NLPModel nlp(m, "cppad");
    coek::NLPSolver solver("ipopt");
    solver.set_option("print_level", 0);
    solver.solve(nlp);

    // x^*_i = -10
    for (size_t i : coek::range(N))
        std::cout << "Value of " << x(i).name() << ": " << x(i).value() << std::endl;
}

The function invquad_array_solve constructs an inverse quadratic problem that is optimized using the ipopt nonlinear optimization solver. The optimal solution is in the opposite corner of the feasible space from the constant values passed into this function.

#include <coek/coek.hpp>
#include <vector>

void invquad_array_resolve()
{
    // Create model
    coek::Model m;

    size_t N = 5;
    auto p = coek::parameter("p", N);
    for (auto& param : p) param.value(0.5);

    // Initialize variables and add them to the model
    auto x = coek::variable("x", N).bounds(-10, 10).value(0.0);
    m.add(x);

    // Create objective and add it to the model
    auto e = coek::expression();
    for (size_t i : coek::range(N)) e -= (x(i) - p(i)) * (x(i) - p(i));
    m.add_objective(e);

    // Optimize the model
    coek::NLPModel nlp(m, "cppad");
    coek::NLPSolver solver("ipopt");
    solver.set_option("print_level", 0);
    solver.load(nlp);

    solver.resolve();
    // x^*_i = -10
    for (size_t i : coek::range(N))
        std::cout << "Value of " << x(i).name() << ": " << x(i).value() << std::endl;

    for (auto& param : p) param.value(-0.5);
    solver.resolve();
    // x^*_i = -10
    for (size_t i : coek::range(N))
        std::cout << "Value of " << x(i).name() << ": " << x(i).value() << std::endl;
}

The function invquad_array_resolve constructs an inverse quadratic problem that is successively re-optimized using the ipopt nonlinear optimization solver. The initial optimal solution is [-10, -10, -10, -10, -10], which is maximally far away from the initial parameter values. Resolving after changing the parameter values does not move the solution, since the resolve uses an initial value matching the last solution.