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_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 168                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 169                if (ret != NLOPT_SUCCESS) goto done;
 170           }
 171
 172           /* "selection" step: rank the population */
 173           for (k = 0; k < population; ++k) irank[k] = k;
 174           if (all_feasible) /* special case: rank by objective function */
 175                nlopt_qsort_r(irank, population, sizeof(int), fval,key_compare);
 176           else {
 177                /* Runarsson & Yao's stochastic ranking of the population */
 178                for (i = 0; i < population; ++i) {
 179                     int swapped = 0;
 180                     for (j = 0; j < population-1; ++j) {
 181                          double u = nlopt_urand(0,1);
 182                          if (u < PF || (penalty[irank[j]] == 0
 183                                         && penalty[irank[j+1]] == 0)) {
 184                               if (fval[irank[j]] > fval[irank[j+1]]) {
 185                                    int irankj = irank[j];
 186                                    irank[j] = irank[j+1];
 187                                    irank[j+1] = irankj;
 188                                    swapped = 1;
 189                               }
 190                          }
 191                          else if (penalty[irank[j]] > penalty[irank[j+1]]) {
 192                               int irankj = irank[j];
 193                               irank[j] = irank[j+1];
 194                               irank[j+1] = irankj;
 195                               swapped = 1;
 196                          }
 197                     }
 198                     if (!swapped) break;
 199                }
 200           }
 201
 202           /* evolve the population:
 203              differential evolution for the best survivors,
 204              and standard mutation of the best survivors for the rest: */
 205           for (k = survivors; k < population; ++k) { /* standard mutation */
 206                double taup_rand = taup * nlopt_nrand(0,1);
 207                int rk = irank[k], ri;
 208                i = k % survivors;
 209                ri = irank[i];
 210                for (j = 0; j < n; ++j) {
 211                     double sigmamax = (ub[j] - lb[j]) / sqrt(n);
 212                     sigmas[rk*n+j] = sigmas[ri*n+j]
 213                          * exp(taup_rand + tau*nlopt_nrand(0,1));
 214                     if (sigmas[rk*n+j] > sigmamax)
 215                          sigmas[rk*n+j] = sigmamax;
 216                     do {
 217                          xs[rk*n+j] = xs[ri*n+j]
 218                               + sigmas[rk*n+j] * nlopt_nrand(0,1);
 219                     } while (xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]);
 220                     sigmas[rk*n+j] = sigmas[ri*n+j] + ALPHA*(sigmas[rk*n+j]
 221                                                            - sigmas[ri*n+j]);
 222                }
 223           }
 224           memcpy(x0, xs, n * sizeof(double));
 225           for (k = 0; k < survivors; ++k) { /* differential variation */
 226                double taup_rand = taup * nlopt_nrand(0,1);
 227                int rk = irank[k];
 228                for (j = 0; j < n; ++j) {
 229                     double xi = xs[rk*n+j];
 230                     if (k+1 < survivors)
 231                          xs[rk*n+j] += GAMMA * (x0[j] - xs[(k+1)*n+j]);
 232                     if (k+1 == survivors
 233                         || xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]) {
 234                          /* standard mutation for last survivor and
 235                             for any survivor components that are now
 236                             outside the bounds */
 237                          double sigmamax = (ub[j] - lb[j]) / sqrt(n);
 238                          double sigi = sigmas[rk*n+j];
 239                          sigmas[rk*n+j] *= exp(taup_rand
 240                                                + tau*nlopt_nrand(0,1));
 241                          if (sigmas[rk*n+j] > sigmamax)
 242                               sigmas[rk*n+j] = sigmamax;
 243                          do {
 244                               xs[rk*n+j] = xi
 245                                    + sigmas[rk*n+j] * nlopt_nrand(0,1);
 246                          } while (xs[rk*n+j] < lb[j] || xs[rk*n+j] > ub[j]);
 247                          sigmas[rk*n+j] = sigi
 248                               + ALPHA * (sigmas[rk*n+j] - sigi);
 249                     }
 250                }
 251           }
 252      }
 253
 254 done:
 255      if (irank) free(irank);
 256      if (sigmas) free(sigmas);
 257      return ret;
 258 }