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 nlopt_result isres_minimize(int n, nlopt_func f, void *f_data,
  57                             int m, nlopt_constraint *fc, /* fc <= 0 */
  58                             int p, nlopt_constraint *h, /* h == 0 */
  59                             const double *lb, const double *ub, /* bounds */
  60                             double *x, /* in: initial guess, out: minimizer */
  61                             double *minf,
  62                             nlopt_stopping *stop,
  63                             int population) /* pop. size (= 0 for default) */
  64 {
  65      const double ALPHA = 0.2; /* smoothing factor from paper */
  66      const double GAMMA = 0.85; /* step-reduction factor from paper */
  67      const double PHI = 1.0; /* expected rate of convergence, from paper */
  68      const double PF = 0.45; /* fitness probability, from paper */
  69      const double SURVIVOR = 1.0/7.0; /* survivor fraction, from paper */
  70      int survivors;
  71      nlopt_result ret = NLOPT_SUCCESS;
  72      double *sigmas = 0, *xs; /* population-by-n arrays (row-major) */
  73      double *fval; /* population array of function vals */
  74      double *penalty; /* population array of penalty vals */
  75      double *x0;
  76      int *irank = 0;
  77      int k, i, j, c;
  78      int mp = m + p;
  79      double minf_penalty = HUGE_VAL, minf_gpenalty = HUGE_VAL;
  80      double taup, tau;
  81
  82      *minf = HUGE_VAL;
  83
  84      if (!population) population = 20 * (n + 1);
  85      if (population < 1) return NLOPT_INVALID_ARGS;
  86      survivors = ceil(population * SURVIVOR);
  87
  88      taup = PHI / sqrt(2*n);
  89      tau = PHI / sqrt(2*sqrt(n));
  90
  91      /* we don't handle unbounded search regions */
  92      for (j = 0; j < n; ++j) if (nlopt_isinf(lb[j]) || nlopt_isinf(ub[j]))
  93                                   return NLOPT_INVALID_ARGS;
  94
  95      sigmas = (double*) malloc(sizeof(double) * (population*n*2
  96                                                  + population
  97                                                  + population
  98                                                  + n));
  99      if (!sigmas) return NLOPT_OUT_OF_MEMORY;
 100      xs = sigmas + population*n;
 101      fval = xs + population*n;
 102      penalty = fval + population;
 103      x0 = penalty + population;
 104
 105      irank = (int*) malloc(sizeof(int) * population);
 106      if (!irank) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 107
 108      for (k = 0; k < population; ++k) {
 109           for (j = 0; j < n; ++j) {
 110                sigmas[k*n+j] = (ub[j] - lb[j]) / sqrt(n);
 111                xs[k*n+j] = nlopt_urand(lb[j], ub[j]);
 112           }
 113      }
 114      memcpy(xs, x, sizeof(double) * n); /* use input x for xs_0 */
 115
 116      while (1) { /* each loop body = one generation */
 117           int all_feasible = 1;
 118
 119           /* evaluate f and constraint violations for whole population */
 120           for (k = 0; k < population; ++k) {
 121                int feasible = 1;
 122                double gpenalty;
 123                stop->nevals++;
 124                fval[k] = f(n, xs + k*n, NULL, f_data);
 125                penalty[k] = 0;
 126                for (c = 0; c < m; ++c) { /* inequality constraints */
 127                     double gval = fc[c].f(n, xs + k*n, NULL, fc[c].f_data);
 128                     if (gval > fc[c].tol) feasible = 0;
 129                     if (gval < 0) gval = 0;
 130                     penalty[k] += gval*gval;
 131                }
 132                gpenalty = penalty[k];
 133                for (c = m; c < mp; ++c) { /* equality constraints */
 134                     double hval = h[c-m].f(n, xs + k*n, NULL, h[c-m].f_data);
 135                     if (fabs(hval) > h[c-m].tol) feasible = 0;
 136                     penalty[k] += hval*hval;
 137                }
 138                if (penalty[k] > 0) all_feasible = 0;
 139
 140                /* convergence criteria (FIXME: improve?) */
 141
 142                /* FIXME: with equality constraints, how do
 143                   we decide which solution is the "best" so far?
 144                   ... need some total order on the solutions? */
 145
 146                if ((penalty[k] <= minf_penalty || feasible)
 147                    && (fval[k] <= *minf || minf_gpenalty > 0)
 148                    && ((feasible ? 0 : penalty[k]) != minf_penalty
 149                        || fval[k] != *minf)) {
 150                     if (fval[k] < stop->minf_max && feasible)
 151                          ret = NLOPT_MINF_MAX_REACHED;
 152                     else if (!nlopt_isinf(*minf)) {
 153                          if (nlopt_stop_f(stop, fval[k], *minf)
 154                              && nlopt_stop_f(stop, feasible ? 0 : penalty[k],
 155                                              minf_penalty))
 156                               ret = NLOPT_FTOL_REACHED;
 157                          else if (nlopt_stop_x(stop, xs+k*n, x))
 158                               ret = NLOPT_XTOL_REACHED;
 159                     }
 160                     memcpy(x, xs+k*n, sizeof(double)*n);
 161                     *minf = fval[k];
 162                     minf_penalty = feasible ? 0 : penalty[k];
 163                     minf_gpenalty = feasible ? 0 : gpenalty;
 164                     if (ret != NLOPT_SUCCESS) goto done;
 165                }
 166
 167                if (nlopt_stop_forced(stop)) ret = NLOPT_FORCE_STOP;
 168                else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 169                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 170                if (ret != NLOPT_SUCCESS) goto done;
 171           }
 172
 173           /* "selection" step: rank the population */
 174           for (k = 0; k < population; ++k) irank[k] = k;
 175           if (all_feasible) /* special case: rank by objective function */
 176                nlopt_qsort_r(irank, population, sizeof(int), fval,key_compare);
 177           else {
 178                /* Runarsson & Yao's stochastic ranking of the population */
 179                for (i = 0; i < population; ++i) {
 180                     int swapped = 0;
 181                     for (j = 0; j < population-1; ++j) {
 182                          double u = nlopt_urand(0,1);
 183                          if (u < PF || (penalty[irank[j]] == 0
 184                                         && penalty[irank[j+1]] == 0)) {
 185                               if (fval[irank[j]] > fval[irank[j+1]]) {
 186                                    int irankj = irank[j];
 187                                    irank[j] = irank[j+1];
 188                                    irank[j+1] = irankj;
 189                                    swapped = 1;
 190                               }
 191                          }
 192                          else if (penalty[irank[j]] > penalty[irank[j+1]]) {
 193                               int irankj = irank[j];
 194                               irank[j] = irank[j+1];
 195                               irank[j+1] = irankj;
 196                               swapped = 1;
 197                          }
 198                     }
 199                     if (!swapped) break;
 200                }
 201           }
 202
 203           /* evolve the population:
 204              differential evolution for the best survivors,
 205              and standard mutation of the best survivors for the rest: */
 206           for (k = survivors; k < population; ++k) { /* standard mutation */
 207                double taup_rand = taup * nlopt_nrand(0,1);
 208                int rk = irank[k], ri;
 209                i = k % survivors;
 210                ri = irank[i];
 211                for (j = 0; j < n; ++j) {
 212                     double sigmamax = (ub[j] - lb[j]) / sqrt(n);
 213                     sigmas[rk*n+j] = sigmas[ri*n+j]
 214                          * exp(taup_rand + tau*nlopt_nrand(0,1));
 215                     if (sigmas[rk*n+j] > sigmamax)
 216                          sigmas[rk*n+j] = sigmamax;
 217                     do {
 218                          xs[rk*n+j] = xs[ri*n+j]
 219                               + sigmas[rk*n+j] * nlopt_nrand(0,1);
 220                     } while (xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]);
 221                     sigmas[rk*n+j] = sigmas[ri*n+j] + ALPHA*(sigmas[rk*n+j]
 222                                                            - sigmas[ri*n+j]);
 223                }
 224           }
 225           memcpy(x0, xs, n * sizeof(double));
 226           for (k = 0; k < survivors; ++k) { /* differential variation */
 227                double taup_rand = taup * nlopt_nrand(0,1);
 228                int rk = irank[k];
 229                for (j = 0; j < n; ++j) {
 230                     double xi = xs[rk*n+j];
 231                     if (k+1 < survivors)
 232                          xs[rk*n+j] += GAMMA * (x0[j] - xs[(k+1)*n+j]);
 233                     if (k+1 == survivors
 234                         || xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]) {
 235                          /* standard mutation for last survivor and
 236                             for any survivor components that are now
 237                             outside the bounds */
 238                          double sigmamax = (ub[j] - lb[j]) / sqrt(n);
 239                          double sigi = sigmas[rk*n+j];
 240                          sigmas[rk*n+j] *= exp(taup_rand
 241                                                + tau*nlopt_nrand(0,1));
 242                          if (sigmas[rk*n+j] > sigmamax)
 243                               sigmas[rk*n+j] = sigmamax;
 244                          do {
 245                               xs[rk*n+j] = xi
 246                                    + sigmas[rk*n+j] * nlopt_nrand(0,1);
 247                          } while (xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]);
 248                          sigmas[rk*n+j] = sigi
 249                               + ALPHA * (sigmas[rk*n+j] - sigi);
 250                     }
 251                }
 252           }
 253      }
 254
 255 done:
 256      if (irank) free(irank);
 257      if (sigmas) free(sigmas);
 258      return ret;
 259 }