/* nag_mip_sqp (h02dac) Example Program.
*
* Copyright 2021 Numerical Algorithms Group.
*
* Mark 27.2, 2021.
*/
#include <nag.h>
#include <stdio.h>
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL confun(Integer *mode, Integer ncnln, Integer n,
const Integer varcon[], const double x[],
double c[], double cjac[], Integer nstate,
Nag_Comm *comm);
static void NAG_CALL objfun(Integer *mode, Integer n, const Integer varcon[],
const double x[], double *objmip, double objgrd[],
Integer nstate, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
#define CJAC(I, J) cjac[(J - 1) * ncnln + I - 1]
#define A(I, J) a[(J - 1) * pda + I - 1]
int main(void) {
/* Integer scalar and array declarations */
const Integer liopts = 200, lopts = 100, lcvalue = 40;
Integer i, j, pda, maxit, n, nclin, ncnln, exit_status = 0;
Integer iopts[200], p, *varcon = 0, ivalue;
/* NAG structures and types */
Nag_VariableType optype;
NagError fail;
Nag_Comm comm;
/* Double scalar and array declarations */
double acc, accqp, objmip;
double *a = 0, *ax = 0, *bl = 0, *bu = 0, *c = 0, *cjac = 0;
double *d = 0, *objgrd = 0, *x = 0, opts[200], rho;
static double ruser[2] = {-1.0, -1.0};
/* Character declarations */
char cvalue[40];
/* Initialize the error structure */
INIT_FAIL(fail);
printf("nag_mip_sqp (h02dac) Example Program Results\n\n");
n = 8;
nclin = 5;
ncnln = 2;
pda = nclin;
if (!(a = NAG_ALLOC(n * pda, double)) || !(d = NAG_ALLOC(nclin, double)) ||
!(ax = NAG_ALLOC(nclin, double)) || !(bl = NAG_ALLOC(n, double)) ||
!(bu = NAG_ALLOC(n, double)) ||
!(varcon = NAG_ALLOC(n + nclin + ncnln, Integer)) ||
!(x = NAG_ALLOC(n, double)) || !(c = NAG_ALLOC(ncnln, double)) ||
!(cjac = NAG_ALLOC(ncnln * n, double)) ||
!(objgrd = NAG_ALLOC(n, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 0; i < 4; i++) {
/* Set variable types: continuous then binary */
varcon[i] = 0;
varcon[4 + i] = 1;
/* Set continuous variable bounds */
bl[i] = 0.0;
bu[i] = 1.0e3;
}
/* Bounds of binary variables need not be provided */
for (i = 4; i < 8; i++) {
bl[i] = 0.0;
bu[i] = 1.0;
}
/* Set linear constraint, equality first */
varcon[n] = 3;
varcon[n + 1] = varcon[n + 2] = varcon[n + 3] = varcon[n + 4] = 4;
/* Set Ax=d then Ax>=d */
for (i = 1; i <= nclin; i++) {
for (j = 1; j <= n; j++) {
A(i, j) = 0.0;
}
}
A(1, 1) = A(1, 2) = A(1, 3) = A(1, 4) = 1.0;
A(2, 1) = -1.0;
A(2, 5) = 1.0;
A(3, 2) = -1.0;
A(3, 6) = 1.0;
A(4, 3) = -1.0;
A(4, 7) = 1.0;
A(5, 4) = -1.0;
A(5, 8) = 1.0;
d[0] = 1.0;
d[1] = d[2] = d[3] = d[4] = 0.0;
/* Set constraints supplied by CONFUN, equality first */
varcon[n + nclin] = 3;
varcon[n + nclin + 1] = 4;
/* Initialize communication arrays */
nag_mip_optset("Initialize = nag_mip_sqp", iopts, liopts, opts, lopts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_mip_optset (h02zkc).\n%s\n", fail.message);
exit_status = -1;
goto END;
}
/* Optimisiation parameters */
maxit = 500;
acc = 1.0e-6;
/* Initial estimate (binary variables need not be given) */
x[0] = x[1] = x[2] = x[3] = 1.0;
x[4] = x[5] = x[6] = x[7] = 0.0;
/* Portfolio parameters */
p = 3;
rho = 10.0;
comm.iuser = &p;
ruser[0] = rho;
comm.user = ruser;
/* Call MINLP solver h02dac (nag_mip_sqp) */
nag_mip_sqp(n, nclin, ncnln, a, pda, d, ax, bl, bu, varcon, x, confun, c,
cjac, objfun, objgrd, maxit, acc, &objmip, iopts, opts, &comm,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_mip_sqp (h02dac).\n%s\n", fail.message);
exit_status = -1;
goto END;
}
/* Query the accuracy of the mixed integer QP solver */
nag_mip_optget("QP Accuracy", &ivalue, &accqp, cvalue, lcvalue, &optype,
iopts, opts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_mip_optget (h02zlc).\n%s\n", fail.message);
exit_status = -1;
goto END;
}
/* Results */
printf("\nFinal estimate:");
for (i = 0; i < n; i++) {
printf("\nx[%4" NAG_IFMT "] = %12.4f", i + 1, x[i]);
}
printf("\n\nRequested accuracy of QP subproblems = %12.4g\n", accqp);
printf("\nOptimised value = %12.4g\n", objmip);
END:
NAG_FREE(a);
NAG_FREE(d);
NAG_FREE(ax);
NAG_FREE(bl);
NAG_FREE(bu);
NAG_FREE(varcon);
NAG_FREE(x);
NAG_FREE(c);
NAG_FREE(cjac);
NAG_FREE(objgrd);
return (exit_status);
}
static void NAG_CALL confun(Integer *mode, Integer ncnln, Integer n,
const Integer varcon[], const double x[],
double c[], double cjac[], Integer nstate,
Nag_Comm *comm) {
Integer p;
double rho;
/* This problem has two nonlinear constraints.
* As an example of using the mode mechanism,
* terminate if any other size is supplied.
*/
if (ncnln != 2) {
*mode = -1;
return;
}
if (nstate == 1)
printf("\n(confun was just called for the first time)\n");
if (*mode == 0) {
/* Constraints */
/* Mean return at least rho: */
rho = comm->user[0];
c[0] = 8.0 * x[0] + 9.0 * x[1] + 12.0 * x[2] + 7.0 * x[3] - rho;
/* Maximum of p assets in portfolio: */
p = *(comm->iuser);
c[1] = (double)p - x[4] - x[5] - x[6] - x[7];
} else {
/* Jacobian */
/* c[0] */
CJAC(1, 1) = 8.0;
CJAC(1, 2) = 9.0;
CJAC(1, 3) = 12.0;
CJAC(1, 4) = 7.0;
/* c[1] does not include continuous variables which requires
that their derivatives are zero */
CJAC(2, 1) = CJAC(2, 2) = CJAC(2, 3) = CJAC(2, 4) = 0.0;
}
}
static void NAG_CALL objfun(Integer *mode, Integer n, const Integer varcon[],
const double x[], double *objmip, double objgrd[],
Integer nstate, Nag_Comm *comm) {
/* This is an 8-dimensional problem.
* As an example of using the mode mechanism,
* terminate if any other size is supplied.
*/
if (n != 8) {
*mode = -1;
return;
}
if (nstate == 1 || comm->user[1] == -1.0) {
printf("\n(objfun was just called for the first time)\n");
comm->user[1] = 0.0;
}
if (*mode == 0) {
/* Objective value */
*objmip = x[0] * (4.0 * x[0] + 3.0 * x[1] - x[2]) +
x[1] * (3.0 * x[0] + 6.0 * x[1] + x[2]) +
x[2] * (x[1] - x[0] + 10.0 * x[2]);
} else {
/* Objective gradients for continuous variables */
objgrd[0] = 8.0 * x[0] + 6.0 * x[1] - 2.0 * x[2];
objgrd[1] = 6.0 * x[0] + 12.0 * x[1] + 2.0 * x[2];
objgrd[2] = 2.0 * (x[1] - x[0]) + 20.0 * x[2];
objgrd[3] = 0.0;
}
}