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) {
  90          nlopt_stop_msg(stop, "population %d is too small", population);
  91          return NLOPT_INVALID_ARGS;
  92      }
  93      survivors = (int) ceil(population * SURVIVOR);
  94
  95      taup = PHI / sqrt(2*n);
  96      tau = PHI / sqrt(2*sqrt(n));
  97
  98      /* we don't handle unbounded search regions */
  99      for (j = 0; j < n; ++j) if (nlopt_isinf(lb[j]) || nlopt_isinf(ub[j])) {
 100              nlopt_stop_msg(stop, "isres requires a finite search region");
 101              return NLOPT_INVALID_ARGS;
 102          }
 103
 104      ires = imax2(nlopt_max_constraint_dim(m, fc),
 105                   nlopt_max_constraint_dim(p, h));
 106      results = (double *) malloc(ires * sizeof(double));
 107      if (ires > 0 && !results) return NLOPT_OUT_OF_MEMORY;
 108
 109      sigmas = (double*) malloc(sizeof(double) * (population*n*2
 110                                                  + population
 111                                                  + population
 112                                                  + n));
 113      if (!sigmas) { free(results); return NLOPT_OUT_OF_MEMORY; }
 114      xs = sigmas + population*n;
 115      fval = xs + population*n;
 116      penalty = fval + population;
 117      x0 = penalty + population;
 118
 119      irank = (int*) malloc(sizeof(int) * population);
 120      if (!irank) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 121
 122      for (k = 0; k < population; ++k) {
 123           for (j = 0; j < n; ++j) {
 124                sigmas[k*n+j] = (ub[j] - lb[j]) / sqrt(n);
 125                xs[k*n+j] = nlopt_urand(lb[j], ub[j]);
 126           }
 127      }
 128      memcpy(xs, x, sizeof(double) * n); /* use input x for xs_0 */
 129
 130      while (1) { /* each loop body = one generation */
 131           int all_feasible = 1;
 132
 133           /* evaluate f and constraint violations for whole population */
 134           for (k = 0; k < population; ++k) {
 135                int feasible = 1;
 136                double gpenalty;
 137                ++ *(stop->nevals_p);
 138                fval[k] = f(n, xs + k*n, NULL, f_data);
 139                if (nlopt_stop_forced(stop)) {
 140                     ret = NLOPT_FORCED_STOP; goto done; }
 141                penalty[k] = 0;
 142                for (c = 0; c < m; ++c) { /* inequality constraints */
 143                     nlopt_eval_constraint(results, NULL,
 144                                           fc + c, n, xs + k*n);
 145                     if (nlopt_stop_forced(stop)) {
 146                          ret = NLOPT_FORCED_STOP; goto done; }
 147                     for (ires = 0; ires < fc[c].m; ++ires) {
 148                          double gval = results[ires];
 149                          if (gval > fc[c].tol[ires]) feasible = 0;
 150                          if (gval < 0) gval = 0;
 151                          penalty[k] += gval*gval;
 152                     }
 153                }
 154                gpenalty = penalty[k];
 155                for (c = m; c < mp; ++c) { /* equality constraints */
 156                     nlopt_eval_constraint(results, NULL,
 157                                           h + (c-m), n, xs + k*n);
 158                     if (nlopt_stop_forced(stop)) {
 159                          ret = NLOPT_FORCED_STOP; goto done; }
 160                     for (ires = 0; ires < h[c-m].m; ++ires) {
 161                          double hval = results[ires];
 162                          if (fabs(hval) > h[c-m].tol[ires]) feasible = 0;
 163                          penalty[k] += hval*hval;
 164                     }
 165                }
 166                if (penalty[k] > 0) all_feasible = 0;
 167
 168                /* convergence criteria (FIXME: improve?) */
 169
 170                /* FIXME: with equality constraints, how do
 171                   we decide which solution is the "best" so far?
 172                   ... need some total order on the solutions? */
 173
 174                if ((penalty[k] <= minf_penalty || feasible)
 175                    && (fval[k] <= *minf || minf_gpenalty > 0)
 176                    && ((feasible ? 0 : penalty[k]) != minf_penalty
 177                        || fval[k] != *minf)) {
 178                     if (fval[k] < stop->minf_max && feasible)
 179                          ret = NLOPT_MINF_MAX_REACHED;
 180                     else if (!nlopt_isinf(*minf)) {
 181                          if (nlopt_stop_f(stop, fval[k], *minf)
 182                              && nlopt_stop_f(stop, feasible ? 0 : penalty[k],
 183                                              minf_penalty))
 184                               ret = NLOPT_FTOL_REACHED;
 185                          else if (nlopt_stop_x(stop, xs+k*n, x))
 186                               ret = NLOPT_XTOL_REACHED;
 187                     }
 188                     memcpy(x, xs+k*n, sizeof(double)*n);
 189                     *minf = fval[k];
 190                     minf_penalty = feasible ? 0 : penalty[k];
 191                     minf_gpenalty = feasible ? 0 : gpenalty;
 192                     if (ret != NLOPT_SUCCESS) goto done;
 193                }
 194
 195                if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 196                else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 197                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 198                if (ret != NLOPT_SUCCESS) goto done;
 199           }
 200
 201           /* "selection" step: rank the population */
 202           for (k = 0; k < population; ++k) irank[k] = k;
 203           if (all_feasible) /* special case: rank by objective function */
 204                nlopt_qsort_r(irank, population, sizeof(int), fval,key_compare);
 205           else {
 206                /* Runarsson & Yao's stochastic ranking of the population */
 207                for (i = 0; i < population; ++i) {
 208                     int swapped = 0;
 209                     for (j = 0; j < population-1; ++j) {
 210                          double u = nlopt_urand(0,1);
 211                          if (u < PF || (penalty[irank[j]] == 0
 212                                         && penalty[irank[j+1]] == 0)) {
 213                               if (fval[irank[j]] > fval[irank[j+1]]) {
 214                                    int irankj = irank[j];
 215                                    irank[j] = irank[j+1];
 216                                    irank[j+1] = irankj;
 217                                    swapped = 1;
 218                               }
 219                          }
 220                          else if (penalty[irank[j]] > penalty[irank[j+1]]) {
 221                               int irankj = irank[j];
 222                               irank[j] = irank[j+1];
 223                               irank[j+1] = irankj;
 224                               swapped = 1;
 225                          }
 226                     }
 227                     if (!swapped) break;
 228                }
 229           }
 230
 231           /* evolve the population:
 232              differential evolution for the best survivors,
 233              and standard mutation of the best survivors for the rest: */
 234           for (k = survivors; k < population; ++k) { /* standard mutation */
 235                double taup_rand = taup * nlopt_nrand(0,1);
 236                int rk = irank[k], ri;
 237                i = k % survivors;
 238                ri = irank[i];
 239                for (j = 0; j < n; ++j) {
 240                     double sigmamax = (ub[j] - lb[j]) / sqrt(n);
 241                     sigmas[rk*n+j] = sigmas[ri*n+j]
 242                          * exp(taup_rand + tau*nlopt_nrand(0,1));
 243                     if (sigmas[rk*n+j] > sigmamax)
 244                          sigmas[rk*n+j] = sigmamax;
 245                     do {
 246                          xs[rk*n+j] = xs[ri*n+j]
 247                               + sigmas[rk*n+j] * nlopt_nrand(0,1);
 248                     } while (xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]);
 249                     sigmas[rk*n+j] = sigmas[ri*n+j] + ALPHA*(sigmas[rk*n+j]
 250                                                            - sigmas[ri*n+j]);
 251                }
 252           }
 253           memcpy(x0, xs, n * sizeof(double));
 254           for (k = 0; k < survivors; ++k) { /* differential variation */
 255                double taup_rand = taup * nlopt_nrand(0,1);
 256                int rk = irank[k];
 257                for (j = 0; j < n; ++j) {
 258                     double xi = xs[rk*n+j];
 259                     if (k+1 < survivors)
 260                          xs[rk*n+j] += GAMMA * (x0[j] - xs[(k+1)*n+j]);
 261                     if (k+1 == survivors
 262                         || xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]) {
 263                          /* standard mutation for last survivor and
 264                             for any survivor components that are now
 265                             outside the bounds */
 266                          double sigmamax = (ub[j] - lb[j]) / sqrt(n);
 267                          double sigi = sigmas[rk*n+j];
 268                          sigmas[rk*n+j] *= exp(taup_rand
 269                                                + tau*nlopt_nrand(0,1));
 270                          if (sigmas[rk*n+j] > sigmamax)
 271                               sigmas[rk*n+j] = sigmamax;
 272                          do {
 273                               xs[rk*n+j] = xi
 274                                    + sigmas[rk*n+j] * nlopt_nrand(0,1);
 275                          } while (xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]);
 276                          sigmas[rk*n+j] = sigi
 277                               + ALPHA * (sigmas[rk*n+j] - sigi);
 278                     }
 279                }
 280           }
 281      }
 282
 283 done:
 284      if (irank) free(irank);
 285      if (sigmas) free(sigmas);
 286      if (results) free(results);
 287      return ret;
 288 }