isres/isres.c

   1 /* Copyright (c) 2010 Massachusetts Institute of Technology
   2  *
   3  * Permission is hereby granted, free of charge, to any person obtaining
   4  * a copy of this software and associated documentation files (the
   5  * "Software"), to deal in the Software without restriction, including
   6  * without limitation the rights to use, copy, modify, merge, publish,
   7  * distribute, sublicense, and/or sell copies of the Software, and to
   8  * permit persons to whom the Software is furnished to do so, subject to
   9  * the following conditions:
  10  *
  11  * The above copyright notice and this permission notice shall be
  12  * included in all copies or substantial portions of the Software.
  13  *
  14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  16  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  17  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  18  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  19  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  20  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  21  */
  22
  23 #include <stdlib.h>
  24 #include <math.h>
  25 #include <string.h>
  26
  27 #include "isres.h"
  28
  29 /* Improved Stochastic Ranking Evolution Strategy (ISRES) algorithm
  30    for nonlinearly-constrained global optimization, based on
  31    the method described in:
  32
  33    Thomas Philip Runarsson and Xin Yao, "Search biases in constrained
  34    evolutionary optimization," IEEE Trans. on Systems, Man, and Cybernetics
  35    Part C: Applications and Reviews, vol. 35 (no. 2), pp. 233-243 (2005).
  36
  37    It is a refinement of an earlier method described in:
  38
  39    Thomas P. Runarsson and Xin Yao, "Stochastic ranking for constrained
  40    evolutionary optimization," IEEE Trans. Evolutionary Computation,
  41    vol. 4 (no. 3), pp. 284-294 (2000).
  42
  43    This is an independent implementation by S. G. Johnson (2009) based
  44    on the papers above.  Runarsson also has his own Matlab implemention
  45    available from his web page: http://www3.hi.is/~tpr
  46 */
  47
  48 static int key_compare(void *keys_, const void *a_, const void *b_)
  49 {
  50      const double *keys = (const double *) keys_;
  51      const int *a = (const int *) a_;
  52      const int *b = (const int *) b_;
  53      return keys[*a] < keys[*b] ? -1 : (keys[*a] > keys[*b] ? +1 : 0);
  54 }
  55
  56 static unsigned imax2(unsigned a, unsigned b) { return (a > b ? a : b); }
  57
  58 nlopt_result isres_minimize(int n, nlopt_func f, void *f_data,
  59                             int m, nlopt_constraint *fc, /* fc <= 0 */
  60                             int p, nlopt_constraint *h, /* h == 0 */
  61                             const double *lb, const double *ub, /* bounds */
  62                             double *x, /* in: initial guess, out: minimizer */
  63                             double *minf,
  64                             nlopt_stopping *stop,
  65                             int population) /* pop. size (= 0 for default) */
  66 {
  67      const double ALPHA = 0.2; /* smoothing factor from paper */
  68      const double GAMMA = 0.85; /* step-reduction factor from paper */
  69      const double PHI = 1.0; /* expected rate of convergence, from paper */
  70      const double PF = 0.45; /* fitness probability, from paper */
  71      const double SURVIVOR = 1.0/7.0; /* survivor fraction, from paper */
  72      int survivors;
  73      nlopt_result ret = NLOPT_SUCCESS;
  74      double *sigmas = 0, *xs; /* population-by-n arrays (row-major) */
  75      double *fval; /* population array of function vals */
  76      double *penalty; /* population array of penalty vals */
  77      double *x0;
  78      int *irank = 0;
  79      int k, i, j, c;
  80      int mp = m + p;
  81      double minf_penalty = HUGE_VAL, minf_gpenalty = HUGE_VAL;
  82      double taup, tau;
  83      double *results = 0; /* scratch space for mconstraint results */
  84      unsigned ires;
  85
  86      *minf = HUGE_VAL;
  87
  88      if (!population) population = 20 * (n + 1);
  89      if (population < 1) return NLOPT_INVALID_ARGS;
  90      survivors = ceil(population * SURVIVOR);
  91
  92      taup = PHI / sqrt(2*n);
  93      tau = PHI / sqrt(2*sqrt(n));
  94
  95      /* we don't handle unbounded search regions */
  96      for (j = 0; j < n; ++j) if (nlopt_isinf(lb[j]) || nlopt_isinf(ub[j]))
  97                                   return NLOPT_INVALID_ARGS;
  98
  99      ires = imax2(nlopt_max_constraint_dim(m, fc),
 100                   nlopt_max_constraint_dim(p, h));
 101      results = (double *) malloc(ires * sizeof(double));
 102      if (ires > 0 && !results) return NLOPT_OUT_OF_MEMORY;
 103
 104      sigmas = (double*) malloc(sizeof(double) * (population*n*2
 105                                                  + population
 106                                                  + population
 107                                                  + n));
 108      if (!sigmas) { free(results); return NLOPT_OUT_OF_MEMORY; }
 109      xs = sigmas + population*n;
 110      fval = xs + population*n;
 111      penalty = fval + population;
 112      x0 = penalty + population;
 113
 114      irank = (int*) malloc(sizeof(int) * population);
 115      if (!irank) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 116
 117      for (k = 0; k < population; ++k) {
 118           for (j = 0; j < n; ++j) {
 119                sigmas[k*n+j] = (ub[j] - lb[j]) / sqrt(n);
 120                xs[k*n+j] = nlopt_urand(lb[j], ub[j]);
 121           }
 122      }
 123      memcpy(xs, x, sizeof(double) * n); /* use input x for xs_0 */
 124
 125      while (1) { /* each loop body = one generation */
 126           int all_feasible = 1;
 127
 128           /* evaluate f and constraint violations for whole population */
 129           for (k = 0; k < population; ++k) {
 130                int feasible = 1;
 131                double gpenalty;
 132                stop->nevals++;
 133                fval[k] = f(n, xs + k*n, NULL, f_data);
 134                if (nlopt_stop_forced(stop)) {
 135                     ret = NLOPT_FORCED_STOP; goto done; }
 136                penalty[k] = 0;
 137                for (c = 0; c < m; ++c) { /* inequality constraints */
 138                     nlopt_eval_constraint(results, NULL,
 139                                           fc + c, n, xs + k*n);
 140                     if (nlopt_stop_forced(stop)) {
 141                          ret = NLOPT_FORCED_STOP; goto done; }
 142                     for (ires = 0; ires < fc[c].m; ++ires) {
 143                          double gval = results[ires];
 144                          if (gval > fc[c].tol[ires]) feasible = 0;
 145                          if (gval < 0) gval = 0;
 146                          penalty[k] += gval*gval;
 147                     }
 148                }
 149                gpenalty = penalty[k];
 150                for (c = m; c < mp; ++c) { /* equality constraints */
 151                     nlopt_eval_constraint(results, NULL,
 152                                           h + (c-m), n, xs + k*n);
 153                     if (nlopt_stop_forced(stop)) {
 154                          ret = NLOPT_FORCED_STOP; goto done; }
 155                     for (ires = 0; ires < h[c-m].m; ++ires) {
 156                          double hval = results[ires];
 157                          if (fabs(hval) > h[c-m].tol[ires]) feasible = 0;
 158                          penalty[k] += hval*hval;
 159                     }
 160                }
 161                if (penalty[k] > 0) all_feasible = 0;
 162
 163                /* convergence criteria (FIXME: improve?) */
 164
 165                /* FIXME: with equality constraints, how do
 166                   we decide which solution is the "best" so far?
 167                   ... need some total order on the solutions? */
 168
 169                if ((penalty[k] <= minf_penalty || feasible)
 170                    && (fval[k] <= *minf || minf_gpenalty > 0)
 171                    && ((feasible ? 0 : penalty[k]) != minf_penalty
 172                        || fval[k] != *minf)) {
 173                     if (fval[k] < stop->minf_max && feasible)
 174                          ret = NLOPT_MINF_MAX_REACHED;
 175                     else if (!nlopt_isinf(*minf)) {
 176                          if (nlopt_stop_f(stop, fval[k], *minf)
 177                              && nlopt_stop_f(stop, feasible ? 0 : penalty[k],
 178                                              minf_penalty))
 179                               ret = NLOPT_FTOL_REACHED;
 180                          else if (nlopt_stop_x(stop, xs+k*n, x))
 181                               ret = NLOPT_XTOL_REACHED;
 182                     }
 183                     memcpy(x, xs+k*n, sizeof(double)*n);
 184                     *minf = fval[k];
 185                     minf_penalty = feasible ? 0 : penalty[k];
 186                     minf_gpenalty = feasible ? 0 : gpenalty;
 187                     if (ret != NLOPT_SUCCESS) goto done;
 188                }
 189
 190                if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 191                else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 192                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 193                if (ret != NLOPT_SUCCESS) goto done;
 194           }
 195
 196           /* "selection" step: rank the population */
 197           for (k = 0; k < population; ++k) irank[k] = k;
 198           if (all_feasible) /* special case: rank by objective function */
 199                nlopt_qsort_r(irank, population, sizeof(int), fval,key_compare);
 200           else {
 201                /* Runarsson & Yao's stochastic ranking of the population */
 202                for (i = 0; i < population; ++i) {
 203                     int swapped = 0;
 204                     for (j = 0; j < population-1; ++j) {
 205                          double u = nlopt_urand(0,1);
 206                          if (u < PF || (penalty[irank[j]] == 0
 207                                         && penalty[irank[j+1]] == 0)) {
 208                               if (fval[irank[j]] > fval[irank[j+1]]) {
 209                                    int irankj = irank[j];
 210                                    irank[j] = irank[j+1];
 211                                    irank[j+1] = irankj;
 212                                    swapped = 1;
 213                               }
 214                          }
 215                          else if (penalty[irank[j]] > penalty[irank[j+1]]) {
 216                               int irankj = irank[j];
 217                               irank[j] = irank[j+1];
 218                               irank[j+1] = irankj;
 219                               swapped = 1;
 220                          }
 221                     }
 222                     if (!swapped) break;
 223                }
 224           }
 225
 226           /* evolve the population:
 227              differential evolution for the best survivors,
 228              and standard mutation of the best survivors for the rest: */
 229           for (k = survivors; k < population; ++k) { /* standard mutation */
 230                double taup_rand = taup * nlopt_nrand(0,1);
 231                int rk = irank[k], ri;
 232                i = k % survivors;
 233                ri = irank[i];
 234                for (j = 0; j < n; ++j) {
 235                     double sigmamax = (ub[j] - lb[j]) / sqrt(n);
 236                     sigmas[rk*n+j] = sigmas[ri*n+j]
 237                          * exp(taup_rand + tau*nlopt_nrand(0,1));
 238                     if (sigmas[rk*n+j] > sigmamax)
 239                          sigmas[rk*n+j] = sigmamax;
 240                     do {
 241                          xs[rk*n+j] = xs[ri*n+j]
 242                               + sigmas[rk*n+j] * nlopt_nrand(0,1);
 243                     } while (xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]);
 244                     sigmas[rk*n+j] = sigmas[ri*n+j] + ALPHA*(sigmas[rk*n+j]
 245                                                            - sigmas[ri*n+j]);
 246                }
 247           }
 248           memcpy(x0, xs, n * sizeof(double));
 249           for (k = 0; k < survivors; ++k) { /* differential variation */
 250                double taup_rand = taup * nlopt_nrand(0,1);
 251                int rk = irank[k];
 252                for (j = 0; j < n; ++j) {
 253                     double xi = xs[rk*n+j];
 254                     if (k+1 < survivors)
 255                          xs[rk*n+j] += GAMMA * (x0[j] - xs[(k+1)*n+j]);
 256                     if (k+1 == survivors
 257                         || xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]) {
 258                          /* standard mutation for last survivor and
 259                             for any survivor components that are now
 260                             outside the bounds */
 261                          double sigmamax = (ub[j] - lb[j]) / sqrt(n);
 262                          double sigi = sigmas[rk*n+j];
 263                          sigmas[rk*n+j] *= exp(taup_rand
 264                                                + tau*nlopt_nrand(0,1));
 265                          if (sigmas[rk*n+j] > sigmamax)
 266                               sigmas[rk*n+j] = sigmamax;
 267                          do {
 268                               xs[rk*n+j] = xi
 269                                    + sigmas[rk*n+j] * nlopt_nrand(0,1);
 270                          } while (xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]);
 271                          sigmas[rk*n+j] = sigi
 272                               + ALPHA * (sigmas[rk*n+j] - sigi);
 273                     }
 274                }
 275           }
 276      }
 277
 278 done:
 279      if (irank) free(irank);
 280      if (sigmas) free(sigmas);
 281      if (results) free(results);
 282      return ret;
 283 }