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);
89 if (population < 1) return NLOPT_INVALID_ARGS;
90 survivors = ceil(population * SURVIVOR);
92 taup = PHI / sqrt(2*n);
93 tau = PHI / sqrt(2*sqrt(n));
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;
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;
104 sigmas = (double*) malloc(sizeof(double) * (population*n*2
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;
114 irank = (int*) malloc(sizeof(int) * population);
115 if (!irank) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
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]);
123 memcpy(xs, x, sizeof(double) * n); /* use input x for xs_0 */
125 while (1) { /* each loop body = one generation */
126 int all_feasible = 1;
128 /* evaluate f and constraint violations for whole population */
129 for (k = 0; k < population; ++k) {
133 fval[k] = f(n, xs + k*n, NULL, f_data);
134 if (nlopt_stop_forced(stop)) {
135 ret = NLOPT_FORCED_STOP; goto done; }
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;
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;
161 if (penalty[k] > 0) all_feasible = 0;
163 /* convergence criteria (FIXME: improve?) */
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? */
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],
179 ret = NLOPT_FTOL_REACHED;
180 else if (nlopt_stop_x(stop, xs+k*n, x))
181 ret = NLOPT_XTOL_REACHED;
183 memcpy(x, xs+k*n, sizeof(double)*n);
185 minf_penalty = feasible ? 0 : penalty[k];
186 minf_gpenalty = feasible ? 0 : gpenalty;
187 if (ret != NLOPT_SUCCESS) goto done;
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;
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);
201 /* Runarsson & Yao's stochastic ranking of the population */
202 for (i = 0; i < population; ++i) {
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];
215 else if (penalty[irank[j]] > penalty[irank[j+1]]) {
216 int irankj = irank[j];
217 irank[j] = irank[j+1];
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;
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;
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]
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);
252 for (j = 0; j < n; ++j) {
253 double xi = xs[rk*n+j];
255 xs[rk*n+j] += GAMMA * (x0[j] - xs[(k+1)*n+j]);
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;
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);
279 if (irank) free(irank);
280 if (sigmas) free(sigmas);
281 if (results) free(results);