1 /* Copyright (c) 2010 Massachusetts Institute of Technology
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:
11 * The above copyright notice and this permission notice shall be
12 * included in all copies or substantial portions of the Software.
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.
29 /* Improved Stochastic Ranking Evolution Strategy (ISRES) algorithm
30 for nonlinearly-constrained global optimization, based on
31 the method described in:
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).
37 It is a refinement of an earlier method described in:
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).
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
48 static int key_compare(void *keys_, const void *a_, const void *b_)
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);
56 static unsigned imax2(unsigned a, unsigned b) { return (a > b ? a : b); }
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 */
65 int population) /* pop. size (= 0 for default) */
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 */
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 */
81 double minf_penalty = HUGE_VAL, minf_gpenalty = HUGE_VAL;
83 double *results = 0; /* scratch space for mconstraint results */
88 if (!population) population = 20 * (n + 1);
90 nlopt_stop_msg(stop, "population %d is too small", population);
91 return NLOPT_INVALID_ARGS;
93 survivors = (int) ceil(population * SURVIVOR);
95 taup = PHI / sqrt(2*n);
96 tau = PHI / sqrt(2*sqrt(n));
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;
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;
109 sigmas = (double*) malloc(sizeof(double) * (population*n*2
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;
119 irank = (int*) malloc(sizeof(int) * population);
120 if (!irank) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
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]);
128 memcpy(xs, x, sizeof(double) * n); /* use input x for xs_0 */
130 while (1) { /* each loop body = one generation */
131 int all_feasible = 1;
133 /* evaluate f and constraint violations for whole population */
134 for (k = 0; k < population; ++k) {
138 fval[k] = f(n, xs + k*n, NULL, f_data);
139 if (nlopt_stop_forced(stop)) {
140 ret = NLOPT_FORCED_STOP; goto done; }
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;
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;
166 if (penalty[k] > 0) all_feasible = 0;
168 /* convergence criteria (FIXME: improve?) */
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? */
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],
184 ret = NLOPT_FTOL_REACHED;
185 else if (nlopt_stop_x(stop, xs+k*n, x))
186 ret = NLOPT_XTOL_REACHED;
188 memcpy(x, xs+k*n, sizeof(double)*n);
190 minf_penalty = feasible ? 0 : penalty[k];
191 minf_gpenalty = feasible ? 0 : gpenalty;
192 if (ret != NLOPT_SUCCESS) goto done;
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;
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);
206 /* Runarsson & Yao's stochastic ranking of the population */
207 for (i = 0; i < population; ++i) {
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];
220 else if (penalty[irank[j]] > penalty[irank[j+1]]) {
221 int irankj = irank[j];
222 irank[j] = irank[j+1];
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;
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;
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]
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);
257 for (j = 0; j < n; ++j) {
258 double xi = xs[rk*n+j];
260 xs[rk*n+j] += GAMMA * (x0[j] - xs[(k+1)*n+j]);
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;
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);
284 if (irank) free(irank);
285 if (sigmas) free(sigmas);
286 if (results) free(results);