test/box.c

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include "nlopt.h"
   4
   5 /****************************************************************************/
   6 /* test function from M. J. Box, "A new method of constrained optimization
   7    and a comparison with other methods," Computer J. 8 (1), 42-52 (1965) */
   8
   9 int testfuncs_verbose = 1;
  10 int testfuncs_counter = 0;
  11
  12 static double testfuncs_status(int n, const double *x, double f)
  13 {
  14     ++testfuncs_counter;
  15     if (testfuncs_verbose) {
  16         int i;
  17         printf("f_%d (%g", testfuncs_counter, x[0]);
  18         for (i = 1; i < n; ++i)
  19             printf(", %g", x[i]);
  20         printf(") = %g\n", f);
  21     }
  22     return f;
  23 }
  24
  25 #define RETURN(f) return testfuncs_status(n, x, f);
  26
  27 static const double k1 = -145421.402, k2 = 2931.1506, k3 = -40.427932, k4 = 5106.192,
  28     k5 = 15711.36, k6 = -161622.577, k7 = 4176.15328, k8 = 2.8260078,
  29     k9 = 9200.476, k10 = 13160.295, k11 = -21686.9194, k12 = 123.56928,
  30     k13 = -21.1188894, k14 = 706.834, k15 = 2898.573, k16 = 28298.388,
  31     k17 = 60.81096, k18 = 31.242116, k19 = 329.574, k20 = -2882.082,
  32     k21 = 74095.3845, k22 = -306.262544, k23 = 16.243649, k24 = -3094.252,
  33     k25 = -5566.2628, k26 = -26237, k27 = 99, k28 = -0.42, k29 = 1300, k30 = 2100, k31 = 925548.252, k32 = -61968.8432, k33 = 23.3088196, k34 = -27096.648, k35 = -50843.766;
  34
  35 static double box(int n, const double *x, double *grad, void *data)
  36 {
  37     double x1 = x[0], x2 = x[1], x3 = x[2], x4 = x[3], x5 = x[4];
  38     double b, x6, y1, y2, y3, y4, u;
  39     const double a0 = 9, a1 = 15, a2 = 50, a3 = 9.583, a4 = 20, a5 = 15, a6 = 6, a7 = 0.75;
  40     b = x2 + 0.01 * x3;
  41     x6 = (k1 + k2 * x2 + k3 * x3 + k4 * x4 + k5 * x5) * x1;
  42     y1 = k6 + k7 * x2 + k8 * x3 + k9 * x4 + k10 * x5;
  43     y2 = k11 + k12 * x2 + k13 * x3 + k14 * x4 + k15 * x5;
  44     y3 = k16 + k17 * x2 + k18 * x3 + k19 * x4 + k20 * x5;
  45     y4 = k21 + k22 * x2 + k23 * x3 + k24 * x4 + k25 * x5;
  46     u = a2 * y1 + a3 * y2 + a4 * y3 + a5 * y4 + 7840 * a6 - 100000 * a0 - 50800 * b * a7 + k31 + k32 * x2 + k33 * x3 + k34 * x4 + k35 * x5;
  47     if (grad) {
  48         int i;
  49         grad[0] = u + a1 * (k1 + k2 * x2 + k3 * x3 + k4 * x4 + k5 * x5);
  50         grad[1] = x1 * (a2 * k7 + a3 * k12 + a4 * k17 + a5 * k22 - 50800 * a7 + k32) + a1 * (k2 * x1);
  51         grad[2] = x1 * (a2 * k8 + a3 * k13 + a4 * k18 + a5 * k23 - 50800 * a7 * 0.01) + a1 * x1 * k3;
  52         grad[3] = x1 * (a2 * k9 + a3 * k14 + a4 * k19 + a5 * k24) + a1 * x1 * k4;
  53         grad[4] = x1 * (a2 * k10 + a3 * k15 + a4 * k20 + a5 * k25) + a1 * x1 * k5;
  54         for (i = 0; i < 5; ++i)
  55             grad[i] = -grad[i];
  56     }
  57     RETURN(-(u * x1 - 24345 + a1 * x6));
  58 }
  59
  60 static double box_constraint(int n, const double *x, double *grad, void *data)
  61 {
  62     int which_constraint = *((int *) data);
  63     double x1 = x[0], x2 = x[1], x3 = x[2], x4 = x[3], x5 = x[4];
  64     double x6, y1, y2, y3, x7, x8;
  65     int i;
  66
  67     x6 = (k1 + k2 * x2 + k3 * x3 + k4 * x4 + k5 * x5) * x1;
  68     y1 = k6 + k7 * x2 + k8 * x3 + k9 * x4 + k10 * x5;
  69     y2 = k11 + k12 * x2 + k13 * x3 + k14 * x4 + k15 * x5;
  70     y3 = k16 + k17 * x2 + k18 * x3 + k19 * x4 + k20 * x5;
  71     x7 = (y1 + y2 + y3) * x1;
  72     x8 = (k26 + k27 * x2 + k28 * x3 + k29 * x4 + k30 * x5) * x1 + x6 + x7;
  73
  74     if (grad) {
  75         grad[0] = grad[1] = grad[2] = grad[3] = grad[4] = 0;
  76
  77         if (which_constraint != 2 && which_constraint != -2) {
  78             /* grad x6 */
  79             grad[0] += (k1 + k2 * x2 + k3 * x3 + k4 * x4 + k5 * x5);
  80             grad[1] += x1 * k2;
  81             grad[2] += x1 * k3;
  82             grad[3] += x1 * k4;
  83             grad[4] += x1 * k5;
  84         }
  85         if (which_constraint != 1 && which_constraint != -1) {
  86             /* grad x7 */
  87             grad[0] += (y1 + y2 + y3);
  88             grad[1] += x1 * (k7 + k12 + k17);
  89             grad[2] += x1 * (k8 + k13 + k18);
  90             grad[3] += x1 * (k9 + k14 + k19);
  91             grad[4] += x1 * (k10 + k15 + k20);
  92         }
  93         if (which_constraint == 3 || which_constraint == -3) {
  94             /* grad (x8 - x6 - x7) */
  95             grad[0] += k26 + k27 * x2 + k28 * x3 + k29 * x4 + k30 * x5;
  96             grad[1] += x1 * k27;
  97             grad[2] += x1 * k28;
  98             grad[3] += x1 * k29;
  99             grad[4] += x1 * k30;
 100         }
 101     }
 102
 103     if (which_constraint == -1) {
 104         if (grad)
 105             for (i = 0; i < 5; ++i)
 106                 grad[i] = -grad[i];
 107         return -x6;
 108     } else if (which_constraint == 1) {
 109         return x6 - 294000;
 110     } else if (which_constraint == -2) {
 111         if (grad)
 112             for (i = 0; i < 5; ++i)
 113                 grad[i] = -grad[i];
 114         return -x7;
 115     } else if (which_constraint == 2) {
 116         return x7 - 294000;
 117     } else if (which_constraint == -3) {
 118         if (grad)
 119             for (i = 0; i < 5; ++i)
 120                 grad[i] = -grad[i];
 121         return -x8;
 122     } else if (which_constraint == 3) {
 123         return x8 - 277200;
 124     }
 125     return 0;
 126 }
 127
 128 int main(void)
 129 {
 130     const double box_lb[5] = { 0, 1.2, 20, 9, 6.5 };
 131     const double box_ub[5] = { HUGE_VAL, 2.4, 60, 9.3, 7.0 };
 132     const double box_xmin[5] = { 4.53743, 2.4, 60, 9.3, 7.0 };
 133
 134     int cdata[6] = { -1, 1, -2, 2, -3, 3 };
 135     double x[5] = { 2.52, 2, 37.5, 9.25, 6.8 };
 136     double minf;
 137
 138     nlopt_minimize_constrained(NLOPT_LN_COBYLA, 5, box, NULL, 6, box_constraint, cdata, sizeof(int), box_lb, box_ub, x, &minf, -HUGE_VAL, 1e-10, 0, 0, NULL, 0, 0);
 139     printf("found f(%g,%g,%g,%g,%g) = %0.8f after %d iters\n", x[0], x[1], x[2], x[3], x[4], minf, testfuncs_counter);
 140     return 0;
 141 }