test/box.c

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