// example of using:
// nagcpp::opt::handle_solve_ipopt (e04st)
// nagcpp::opt::handle_set_simplebounds (e04rh)
// nagcpp::opt::handle_set_linobj (e04re)
// nagcpp::opt::handle_set_linconstr (e04rj)
// nagcpp::opt::handle_set_nlnconstr (e04rk)
// nagcpp::opt::handle_set_nlnhess (e04rl)
// nagcpp::opt::handle_​print (e04ry)
// in addition, the following routines are implicitly
// called via the communication class:
// nagcpp::opt::handle_init (e04ra)
// nagcpp::opt::handle_opt_set (e04zm)
// nagcpp::opt::handle_free (e04rz)
#include <fstream>
#include <iomanip>
#include <iostream>
#include <math.h>
#include <vector>
#include "e04/nagcpp_class_CommE04RA.hpp"
#include "e04/nagcpp_e04re.hpp"
#include "e04/nagcpp_e04rh.hpp"
#include "e04/nagcpp_e04rj.hpp"
#include "e04/nagcpp_e04rk.hpp"
#include "e04/nagcpp_e04rl.hpp"
#include "e04/nagcpp_e04ry.hpp"
#include "e04/nagcpp_e04st.hpp"
void confun(const std::vector<double> &x,
const nagcpp::types::f77_integer ncnln, std::vector<double> &gx,
nagcpp::types::f77_integer inform) {
gx[0] = 12.0 * x[0] + 11.9 * x[1] + 41.8 * x[2] + 52.1 * x[3] -
1.645 * std::sqrt(0.28 * pow(x[0], 2) + 0.19 * pow(x[1], 2) +
20.5 * pow(x[2], 2) + 0.62 * pow(x[3], 2));
}
void congrd(const std::vector<double> &x,
std::vector<double> &gdx, const nagcpp::types::f77_integer nnzgd) {
double tmp = std::sqrt(0.62 * pow(x[3], 2) + 20.5 * pow(x[2], 2) +
0.19 * pow(x[1], 2) + 0.28 * pow(x[0], 2));
gdx[0] = (12.0 * tmp - 0.4606 * x[0]) / tmp;
gdx[1] = (11.9 * tmp - 0.31255 * x[1]) / tmp;
gdx[2] = (41.8 * tmp - 33.7225 * x[2]) / tmp;
gdx[3] = (52.1 * tmp - 1.0199 * x[3]) / tmp;
}
void hess(const std::vector<double> &x, const nagcpp::types::f77_integer idf,
const double sigma, const std::vector<double> &lamda,
std::vector<double> &hx, nagcpp::types::f77_integer inform) {
bool terminate = true;
std::fill(hx.begin(), hx.end(), 0.0);
if (idf == 0) {
terminate = false;
} else if (idf == 1 || idf == -1) {
double tmp = std::sqrt(0.62 * pow(x[3], 2) + 20.5 * pow(x[2], 2) +
0.19 * pow(x[1], 2) + 0.28 * pow(x[0], 2));
tmp = tmp * (pow(x[3], 2) + 33.064516129032258064 * pow(x[2], 2) +
0.30645161290322580645 * pow(x[1], 2) +
0.45161290322580645161 * pow(x[0], 2));
// 1,1..4
hx[0] = (-0.4606 * pow(x[3], 2) - 15.229516129032258064 * pow(x[2], 2) -
0.14115161290322580645 * pow(x[1], 2)) /
tmp;
hx[1] = (0.14115161290322580645 * x[0] * x[1]) / tmp;
hx[2] = (15.229516129032258064 * x[0] * x[2]) / tmp;
hx[3] = (0.4606 * x[0] * x[3]) / tmp;
// 2,2..4
hx[4] = (-0.31255 * pow(x[3], 2) - 10.334314516129032258 * pow(x[2], 2) -
0.14115161290322580645 * pow(x[0], 2)) /
tmp;
hx[5] = (10.334314516129032258 * x[1] * x[2]) / tmp;
hx[6] = (0.31255 * x[1] * x[3]) / tmp;
// 3,3..4
hx[7] = (-33.7225 * pow(x[3], 2) - 10.334314516129032258 * pow(x[1], 2) -
15.229516129032258065 * pow(x[0], 2)) /
tmp;
hx[8] = (33.7225 * x[2] * x[3]) / tmp;
// 4,4
hx[9] = (-33.7225 * pow(x[2], 2) - 0.31255 * pow(x[1], 2) -
0.4606 * pow(x[0], 2)) /
tmp;
if (idf == -1) {
for (auto ihx = hx.begin(); ihx != hx.end(); ++ihx) {
(*ihx) *= lamda[0];
}
}
terminate = false;
}
if (terminate) {
throw nagcpp::error_handler::CallbackEarlyTermination("hess", "unhandled "
"value of "
"idf");
}
}
typedef std::vector<double>::size_type vsize_t;
int main(void) {
std::cout << "nagcpp::opt::handle_solve_lp_ipm Example" << std::endl
<< std::endl;
try {
vsize_t k = 0;
nagcpp::types::f77_integer nvar = 4;
// initialise the problem handle
nagcpp::opt::CommE04RA handle(nvar);
double bigbnd = 1.0e20;
// set simple bounds, x_i>=0 (e04rh)
std::vector<double> bl(nvar, 0.0);
std::vector<double> bu(nvar, bigbnd);
nagcpp::opt::handle_set_simplebounds(handle, bl, bu);
// set linear objective (e04re)
std::vector<double> coeffs = {24.55, 26.75, 39.00, 40.50};
nagcpp::opt::handle_set_linobj(handle, coeffs);
// set linear constraints (e04rj)
std::vector<double> linbl = {5.0, 1.0};
std::vector<double> linbu = {bigbnd, 1.0};
std::vector<nagcpp::types::f77_integer> irowb = {1, 1, 1, 1, 2, 2, 2, 2};
std::vector<nagcpp::types::f77_integer> icolb = {1, 2, 3, 4, 1, 2, 3, 4};
std::vector<double> b = {2.3, 5.6, 11.1, 1.3, 1.0, 1.0, 1.0, 1.0};
nagcpp::opt::handle_set_linconstr(handle, linbl, linbu, irowb, icolb, b);
// set nonlinear constraints (e04rk)
std::vector<double> nlnbl = {21.0};
std::vector<double> nlnbu = {bigbnd};
std::vector<nagcpp::types::f77_integer> irowgb = {1, 1, 1, 1};
std::vector<nagcpp::types::f77_integer> icolgb = {1, 2, 3, 4};
nagcpp::opt::handle_set_nlnconstr(handle, nlnbl, nlnbu, irowgb, icolgb);
// hessian of nonlinear constraint (e04rl)
std::vector<nagcpp::types::f77_integer> irowh(nvar * (nvar + 1) / 2);
std::vector<nagcpp::types::f77_integer> icolh(nvar * (nvar + 1) / 2);
// NB: i and j are 1-based
for (vsize_t i = 1; i <= static_cast<vsize_t>(nvar); ++i) {
for (vsize_t j = i; j <= static_cast<vsize_t>(nvar); ++j) {
irowh[k] = i;
icolh[k] = j;
++k;
}
}
nagcpp::opt::handle_set_nlnhess(handle, 1, irowh, icolh);
// open some output streams
std::ofstream printing_fs("ex_e04st.out");
std::ofstream monitoring_fs("ex_e04st.mon");
// register the streams to the global IO manager
// this example uses three streams, std::cout, and two files
// an alternative to using the global IO manager is to create
// a local instance and pass that to the routines via the
// optional parameters argument
nagcpp::types::f77_integer nout =
nagcpp::iomanager::GLOBAL_IOMANAGER->register_ostream(std::cout);
nagcpp::types::f77_integer printing_unit_number =
nagcpp::iomanager::GLOBAL_IOMANAGER->register_ostream(printing_fs);
nagcpp::types::f77_integer monitoring_unit_number =
nagcpp::iomanager::GLOBAL_IOMANAGER->register_ostream(monitoring_fs);
// turn on monitoring
handle.MonitoringLevel(5);
handle.MonitoringFile(monitoring_unit_number);
// turn on printing
handle.PrintLevel(2);
handle.PrintFile(printing_unit_number);
// set some control parameters
handle.StopTolerance1(2.5e-8);
handle.OuterIterationLimit(26);
handle.TimeLimit(60);
// initial values for variable
std::vector<double> x(nvar, 1.0);
// nnzu = 2 * (number of variables +
// number linear constraints + number non-linear constraints)
nagcpp::types::f77_integer nnzu =
2 * nvar + 2 * linbl.size() + 2 * nlnbl.size();
std::vector<double> u(nnzu);
std::vector<double> rinfo;
std::vector<double> stats;
// we are going to be using the default values for the optional
// arguments, however we want to print out any warning messages
nagcpp::opt::OptionalE04ST opt;
// summarise the problem being solved (e04ry)
nagcpp::opt::handle_print(handle, nout, "Overview");
// call the solver
handle_solve_ipopt(handle, nullptr, nullptr, confun, congrd, hess, nullptr,
x, u, rinfo, stats, opt);
if (opt.fail.warning_thrown) {
std::cout << std::endl << opt.fail.msg << std::endl << std::endl;
}
std::cout.precision(6);
std::cout.setf(std::ios::scientific);
std::cout << std::endl << "Variables" << std::endl;
for (vsize_t i = 0; i < x.size(); ++i) {
std::cout << " x(" << std::setw(10) << i << ")" << std::setw(17)
<< "=" << std::setw(20) << x[i] << std::endl;
}
std::cout << "Variable bound Lagrange multipliers" << std::endl;
k = 0;
for (vsize_t i = 0; i < static_cast<vsize_t>(nvar); i++) {
std::cout << " zL(" << std::setw(10) << i << ")" << std::setw(17)
<< "=" << std::setw(20) << u[k] << std::endl;
++k;
std::cout << " zU(" << std::setw(10) << i << ")" << std::setw(17)
<< "=" << std::setw(20) << u[k] << std::endl;
++k;
}
std::cout << "Linear constraints Lagrange multipliers" << std::endl;
for (vsize_t i = 0; i < linbl.size(); i++) {
std::cout << " l-(" << std::setw(10) << i << ")" << std::setw(17)
<< "=" << std::setw(20) << u[k] << std::endl;
++k;
std::cout << " l+(" << std::setw(10) << i << ")" << std::setw(17)
<< "=" << std::setw(20) << u[k] << std::endl;
++k;
}
std::cout << "Nonlinear constraints Lagrange multipliers" << std::endl;
for (vsize_t i = 0; i < nlnbl.size(); i++) {
std::cout << " l-(" << std::setw(10) << i << ")" << std::setw(17)
<< "=" << std::setw(20) << u[k] << std::endl;
++k;
std::cout << " l+(" << std::setw(10) << i << ")" << std::setw(17)
<< "=" << std::setw(20) << u[k] << std::endl;
++k;
}
std::cout.precision(7);
std::cout << "At solution, Objective minimum =" << std::setw(20)
<< rinfo[0] << std::endl;
std::cout.precision(2);
std::cout << " Constraint violation =" << std::setw(20)
<< rinfo[1] << std::endl;
std::cout << " Dual infeasibility =" << std::setw(20)
<< rinfo[2] << std::endl;
std::cout << " Complementarity =" << std::setw(20)
<< rinfo[3] << std::endl;
std::cout << " KKT error =" << std::setw(20)
<< rinfo[4] << std::endl;
// round time to the nearest 10 seconds, to avoid suprious differences in results file
std::cout << "Solution took less than "
<< static_cast<int>(std::ceil(stats[7] / 10) * 10) << " seconds"
<< std::endl;
std::cout << " after iterations : "
<< std::setw(10) << int(stats[0]) << std::endl;
std::cout << " after objective evaluations : "
<< std::setw(10) << int(stats[18]) << std::endl;
std::cout << " after objective gradient evaluations : "
<< std::setw(10) << int(stats[4]) << std::endl;
std::cout << " after constraint evaluations : "
<< std::setw(10) << int(stats[19]) << std::endl;
std::cout << " after constraint gradient evaluations : "
<< std::setw(10) << int(stats[20]) << std::endl;
std::cout << " after hessian evaluations : "
<< std::setw(10) << int(stats[3]) << std::endl;
// deregister the streams as we no longer need them
nagcpp::iomanager::GLOBAL_IOMANAGER->deregister_ostream(printing_fs);
nagcpp::iomanager::GLOBAL_IOMANAGER->deregister_ostream(monitoring_fs);
printing_fs.close();
monitoring_fs.close();
} catch (nagcpp::error_handler::Exception &e) {
std::cout << e.msg << std::endl;
return 1;
}
return 0;
}