/* nag_opt_handle_solve_ipopt (e04stc) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.3, 2022.
*/
/* NLP example: linear objective + linear constraint + nonlinear constraint.
* For illustrative purposes, the linear objective will be treated as a
* nonlinear function to show the usage of objfun and objgrd user call-backs.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL objfun(Integer nvar, const double x[], double *fx,
Integer *flag, Nag_Comm *comm);
static void NAG_CALL objgrd(Integer nvar, const double x[], Integer nnzfd,
double fdx[], Integer *flag, Nag_Comm *comm);
static void NAG_CALL confun(Integer nvar, const double x[], Integer ncnln,
double gx[], Integer *flag, Nag_Comm *comm);
static void NAG_CALL congrd(Integer nvar, const double x[], Integer nnzgd,
double gdx[], Integer *flag, Nag_Comm *comm);
static void NAG_CALL hess(Integer nvar, const double x[], Integer ncnln,
Integer idf, double sigma, const double lambda[],
Integer nnzh, double hx[], Integer *flag,
Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
int main(void) {
#define BIGBND 1.0E20
/* Scalars */
Integer exit_status = 0;
Integer i, idlc, idx, j, nnzu, nvar, nclin, ncnln, nnzfd, nnzgd, flag;
double lmtest;
/* Arrays */
Integer iuser[2];
Integer idxfd[4] = {1, 2, 3, 4};
double ruser[4] = {24.55, 26.75, 39.00, 40.50};
double rinfo[100], stats[100], *x = 0, *u = 0;
double *lambda = 0, *fdx = 0, *gdx = 0;
double bl[4] = {0, 0, 0, 0};
double bu[4] = {BIGBND, BIGBND, BIGBND, BIGBND};
double linbl[2] = {5.0, 1.0};
double linbu[2] = {BIGBND, 1.0};
double nlnbl[1] = {21.0};
double nlnbu[1] = {BIGBND};
Integer irowb[8] = {1, 1, 1, 1, 2, 2, 2, 2};
Integer icolb[8] = {1, 2, 3, 4, 1, 2, 3, 4};
Integer irowgd[4] = {1, 1, 1, 1};
Integer icolgd[4] = {1, 2, 3, 4};
Integer irowh[10], icolh[10];
double b[8] = {2.3, 5.6, 11.1, 1.3, 1.0, 1.0, 1.0, 1.0};
void *handle = 0;
/* Nag Types */
NagError fail;
Nag_FileID nout;
Nag_Comm comm;
/* nag_file_open (x04acc).
* Open unit number for reading, writing or appending, and
* associate unit with named file
*/
nag_file_open("", 1, &nout, NAGERR_DEFAULT);
/* nag_file_line_write (x04bac).
* Write formatted record to external file
*/
nag_file_line_write(nout, "nag_opt_handle_solve_ipopt (e04stc) "
"Example Program Results");
/* Problem size */
nvar = 4;
/* Counter for Lagrange multipliers */
nnzu = 0;
/* Objective gradient nonzero elements quantity */
nnzfd = 4;
/* Constraint jacobian nonzero elements quantity */
nnzgd = 4;
/* nag_opt_handle_init (e04rac).
* Initialize an empty problem handle with NVAR variables.
*/
nag_opt_handle_init(&handle, nvar, NAGERR_DEFAULT);
/* nag_opt_handle_set_simplebounds (e04rhc).
* Define bounds on the variables
*/
nag_opt_handle_set_simplebounds(handle, nvar, bl, bu, NAGERR_DEFAULT);
/* Specify the amount of Lagrange mult. required */
nnzu += 2 * nvar;
/* nag_opt_handle_set_nlnobj (e04rgc).
* Define nonlinear objective
*/
nag_opt_handle_set_nlnobj(handle, nnzfd, idxfd, NAGERR_DEFAULT);
/* Add two linear constraints */
nclin = 2;
idlc = 0;
/* nag_opt_handle_set_linconstr (e04rjc).
* Define a block of linear constraints
*/
nag_opt_handle_set_linconstr(handle, nclin, linbl, linbu, nclin * nvar, irowb,
icolb, b, &idlc, NAGERR_DEFAULT);
/* Update the Lagrange mult. count */
nnzu += 2 * nclin;
/* Add one nonlinear constraint */
ncnln = 1;
/* nag_opt_handle_set_nlnconstr (e04rkc).
* Define a block of nonlinear constraints
*/
nag_opt_handle_set_nlnconstr(handle, ncnln, nlnbl, nlnbu, nnzgd, irowgd,
icolgd, NAGERR_DEFAULT);
/* Update the Lagrange mult. count */
nnzu += 2 * ncnln;
/* Define structure of the Hessian, dense upper triangle */
for (idx = 0, i = 1; i <= nvar; i++)
for (j = i; j <= nvar; j++, idx++) {
icolh[idx] = j;
irowh[idx] = i;
}
/* nag_opt_handle_set_nlnhess (e04rlc).
* Define structure of Hessian of objective, constraints or the Lagrangian
*/
nag_opt_handle_set_nlnhess(handle, 1, idx, irowh, icolh, NAGERR_DEFAULT);
nag_opt_handle_set_nlnhess(handle, 0, idx, irowh, icolh, NAGERR_DEFAULT);
if (!(x = NAG_ALLOC(nvar, double)) || !(u = NAG_ALLOC(nnzu, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 0; i < nvar; i++)
x[i] = 1.0;
iuser[0] = 0;
nag_opt_handle_opt_set(handle, "Print Level = 0", NAGERR_DEFAULT);
nag_opt_handle_opt_set(handle, "Outer Iteration Limit = 26", NAGERR_DEFAULT);
nag_opt_handle_opt_set(handle, "Stop Tolerance 1 = 2.5e-8", NAGERR_DEFAULT);
comm.iuser = iuser;
comm.user = ruser;
INIT_FAIL(fail);
/* nag_opt_handle_solve_ipopt (e04stc).
* Call the solver
*/
nag_opt_handle_solve_ipopt(handle, objfun, objgrd, confun, congrd, hess,
NULLFN, nvar, x, nnzu, u, rinfo, stats, &comm,
&fail);
if (fail.code == NE_NOERROR) {
char msg[80];
nag_file_line_write(nout, "\nVariables");
for (i = 0; i < nvar; i++) {
sprintf(msg, " x(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "", x[i]);
nag_file_line_write(nout, msg);
}
nag_file_line_write(nout, "Variable bound Lagrange multipliers");
for (i = 0; i < nvar; i++) {
sprintf(msg, " zL(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "", u[2 * i]);
nag_file_line_write(nout, msg);
sprintf(msg, " zU(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "",
u[2 * i + 1]);
nag_file_line_write(nout, msg);
}
nag_file_line_write(nout, "Linear constraints Lagrange multipliers");
for (i = 0; i < nclin; i++) {
sprintf(msg, " zL(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "",
u[2 * nvar + 2 * i]);
nag_file_line_write(nout, msg);
sprintf(msg, " zU(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "",
u[2 * nvar + 2 * i + 1]);
nag_file_line_write(nout, msg);
}
nag_file_line_write(nout, "Nonlinear constraints Lagrange multipliers");
for (i = 0; i < ncnln; i++) {
sprintf(msg, " zL(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "",
u[2 * (nvar + nclin) + 2 * i]);
nag_file_line_write(nout, msg);
sprintf(msg, " zU(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "",
u[2 * (nvar + nclin) + 2 * i + 1]);
nag_file_line_write(nout, msg);
}
if (!(lambda = NAG_ALLOC(7, double)) || !(fdx = NAG_ALLOC(nnzfd, double)) ||
!(gdx = NAG_ALLOC(nnzgd, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Check of Lagrange multipliers (complementariety)
* Gradient of the objective (also available in data%coeffs(1:nnzfd))
*/
objgrd(nvar, x, nnzfd, fdx, &flag, &comm);
if (flag != 0) {
exit_status = -1;
goto END;
}
/* Gradient of the nonlinear constraint */
congrd(nvar, x, nnzgd, gdx, &flag, &comm);
if (flag != 0) {
exit_status = -1;
goto END;
}
/* Obtain the multipliers with correct sign
* 4 bound constraints, 2 linear constraints, and 1 nonlinear constraint
*/
for (i = 0; i < 7; i++) {
lambda[i] = u[2 * i + 1] - u[2 * i];
}
/* lambda (0..3) mult. related to bounds
* lambda (4..5) mult. related to linear constraints
* lambda (6) mult. related to the nonlinear constraint
*/
nag_file_line_write(nout, "Stationarity test for Lagrange multipliers");
for (i = 0; i < nvar; i++) {
lmtest = fdx[i] + lambda[i] + lambda[4] * b[i] + lambda[5] * b[4 + i] +
lambda[6] * gdx[i];
if (lmtest < 2.5e-8)
sprintf(msg, " lx(%10" NAG_IFMT ")%17s=%16.2E%5s%6s", i + 1, "",
lmtest, "", "Ok ");
else
sprintf(msg, " lx(%10" NAG_IFMT ")%17s=%16.2E%5s%6s", i + 1, "",
lmtest, "", "Failed");
nag_file_line_write(nout, msg);
}
NAG_FREE(lambda);
NAG_FREE(fdx);
NAG_FREE(gdx);
sprintf(msg, "\nAt solution, Objective minimum =%16.4E", rinfo[0]);
nag_file_line_write(nout, msg);
sprintf(msg, " Constraint violation =%6s%10.2E", "", rinfo[1]);
nag_file_line_write(nout, msg);
sprintf(msg, " Dual infeasibility =%6s%10.2E", "", rinfo[2]);
nag_file_line_write(nout, msg);
sprintf(msg, " Complementarity =%6s%10.2E", "", rinfo[3]);
nag_file_line_write(nout, msg);
sprintf(msg, " KKT error =%6s%10.2E", "", rinfo[4]);
nag_file_line_write(nout, msg);
sprintf(msg, " after iterations :%11.0f",
stats[0]);
nag_file_line_write(nout, msg);
sprintf(msg, " after objective evaluations :%11.0f",
stats[18]);
nag_file_line_write(nout, msg);
sprintf(msg, " after objective gradient evaluations :%11.0f",
stats[4]);
nag_file_line_write(nout, msg);
sprintf(msg, " after constraint evaluations :%11.0f",
stats[19]);
nag_file_line_write(nout, msg);
sprintf(msg, " after constraint gradient evaluations :%11.0f",
stats[20]);
nag_file_line_write(nout, msg);
sprintf(msg, " after hessian evaluations :%11.0f",
stats[3]);
nag_file_line_write(nout, msg);
} else {
printf("Error from nag_opt_handle_solve_ipopt (e04stc).\n%s\n",
fail.message);
exit_status = 1;
}
if (handle)
/* nag_opt_handle_free (e04rzc).
* Destroy the problem handle and deallocate all the memory used
*/
nag_opt_handle_free(&handle, NAGERR_DEFAULT);
NAG_FREE(x);
NAG_FREE(u);
END:
return exit_status;
}
/* Subroutine */
static void NAG_CALL objfun(Integer nvar, const double x[], double *fx,
Integer *flag, Nag_Comm *comm) {
*flag = 0;
/* nag_blast_ddot (f16eac).
* Dot product of two vectors, allows scaling and accumulation
*/
nag_blast_ddot(Nag_NoConj, 4, 1.0, x, 1, 0.0, comm->user, 1, fx,
NAGERR_DEFAULT);
}
static void NAG_CALL objgrd(Integer nvar, const double x[], Integer nnzfd,
double fdx[], Integer *flag, Nag_Comm *comm) {
*flag = 0;
/* nag_blast_dge_copy (f16qfc).
* Matrix copy, real rectangular matrix
*/
nag_blast_dge_copy(Nag_ColMajor, Nag_NoTrans, nnzfd, 1, comm->user, nnzfd,
fdx, nnzfd, NAGERR_DEFAULT);
}
static void NAG_CALL confun(Integer nvar, const double x[], Integer ncnln,
double gx[], Integer *flag, Nag_Comm *comm) {
*flag = 0;
gx[0] = 12.0 * x[0] + 11.9 * x[1] + 41.8 * x[2] + 52.1 * x[3] -
1.645 * sqrt(.28 * x[0] * x[0] + .19 * x[1] * x[1] +
20.5 * x[2] * x[2] + .62 * x[3] * x[3]);
}
static void NAG_CALL congrd(Integer nvar, const double x[], Integer nnzgd,
double gdx[], Integer *flag, Nag_Comm *comm) {
double tmp;
*flag = 0;
tmp = sqrt(0.62 * x[3] * x[3] + 20.5 * x[2] * x[2] + 0.19 * x[1] * x[1] +
0.28 * x[0] * x[0]);
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;
}
static void NAG_CALL hess(Integer nvar, const double x[], Integer ncnln,
Integer idf, double sigma, const double lambda[],
Integer nnzh, double hx[], Integer *flag,
Nag_Comm *comm) {
double tmp;
*flag = 0;
/* nag_blast_dload (f16fbc).
* Broadcast scalar into real vector
*/
nag_blast_dload(nnzh, 0.0, hx, 1, NAGERR_DEFAULT);
if (idf == 0)
return; /* objective is linear */
tmp = sqrt(0.62 * x[3] * x[3] + 20.5 * x[2] * x[2] + 0.19 * x[1] * x[1] +
0.28 * x[0] * x[0]);
tmp = tmp * (x[3] * x[3] + 33.064516129032258064 * x[2] * x[2] +
0.30645161290322580645 * x[1] * x[1] +
0.45161290322580645161 * x[0] * x[0]);
/* Col 1 */
hx[0] = (-0.4606 * x[3] * x[3] - 15.229516129032258064 * x[2] * x[2] -
0.14115161290322580645 * x[1] * x[1]) /
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;
/* Col 2 */
hx[4] = (-0.31255 * x[3] * x[3] - 10.334314516129032258 * x[2] * x[2] -
0.14115161290322580645 * x[0] * x[0]) /
tmp;
hx[5] = (10.334314516129032258 * x[1] * x[2]) / tmp;
hx[6] = (0.31255 * x[1] * x[3]) / tmp;
/* Col 3 */
hx[7] = (-33.7225 * x[3] * x[3] - 10.334314516129032258 * x[1] * x[1] -
15.229516129032258065 * x[0] * x[0]) /
tmp;
hx[8] = (33.7225 * x[2] * x[3]) / tmp;
/* Col 4 */
hx[9] =
(-33.7225 * x[2] * x[2] - 0.31255 * x[1] * x[1] - 0.4606 * x[0] * x[0]) /
tmp;
/* nag_blast_daxpby (f16ecc).
* Real weighted vector addition
*/
if (idf == -1)
nag_blast_daxpby(nnzh, 0.0, hx, 1, lambda[0], hx, 1, NAGERR_DEFAULT);
}