1 /* _______________________________________________________________________
2
3 DAKOTA: Design Analysis Kit for Optimization and Terascale Applications
4 Copyright 2014-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Dakota directory.
7 _______________________________________________________________________ */
8
9 //- Class: NOWPACOptimizer
10 //- Description: Implementation code for the NOWPACOptimizer class
11 //- Owner: Mike Eldred
12
13 #include "DakotaModel.hpp"
14 #include "DakotaResponse.hpp"
15 #include "NOWPACOptimizer.hpp"
16 #include "PostProcessModels.hpp"
17 #include "ProblemDescDB.hpp"
18 #include "ParallelLibrary.hpp"
19
20 static const char rcsId[]="@(#) $Id: NOWPACOptimizer.cpp 7029 2010-10-22 00:17:02Z mseldre $";
21
22 namespace Dakota {
23
24
NOWPACOptimizer(ProblemDescDB & problem_db,Model & model)25 NOWPACOptimizer::NOWPACOptimizer(ProblemDescDB& problem_db, Model& model):
26 Optimizer(problem_db, model, std::shared_ptr<TraitsBase>(new NOWPACTraits())),
27 nowpacSolver(numContinuousVars, "nowpac_diagnostics.dat"),
28 nowpacEvaluator(iteratedModel)
29 {
30 nowpacEvaluator.allocate_constraints();
31 nowpacSolver.set_blackbox(nowpacEvaluator,
32 nowpacEvaluator.num_ineq_constraints());
33 initialize_options();
34 }
35
36
NOWPACOptimizer(Model & model)37 NOWPACOptimizer::NOWPACOptimizer(Model& model):
38 Optimizer(MIT_NOWPAC, model, std::shared_ptr<TraitsBase>(new NOWPACTraits())),
39 nowpacSolver(numContinuousVars, "nowpac_diagnostics.dat"),
40 nowpacEvaluator(iteratedModel)
41 {
42 nowpacEvaluator.allocate_constraints();
43 nowpacSolver.set_blackbox(nowpacEvaluator,
44 nowpacEvaluator.num_ineq_constraints());
45 initialize_options();
46 }
47
48
~NOWPACOptimizer()49 NOWPACOptimizer::~NOWPACOptimizer()
50 { }
51
52
initialize_options()53 void NOWPACOptimizer::initialize_options()
54 {
55 // Refer to bit bucket docs: https://bitbucket.org/fmaugust/nowpac
56
57 // Optional: note that we are overridding NOWPAC defaults with Dakota defaults
58 // May want to leave NOWPAC defaults in place if there is no user spec.
59 nowpacSolver.set_option("eta_1",
60 probDescDB.get_real("method.trust_region.contract_threshold") );
61 nowpacSolver.set_option("eta_2",
62 probDescDB.get_real("method.trust_region.expand_threshold") );
63 // Criticality measures:
64 //nowpacSolver.set_option("eps_c" , 1e-6 );
65 //nowpacSolver.set_option("mu" , 1e1 );
66 // Upper bound on poisedness constant augmented with distance penalty:
67 //nowpacSolver.set_option("geometry_threshold" , 5e2 );
68 nowpacSolver.set_option("gamma_inc",
69 probDescDB.get_real("method.trust_region.expansion_factor") );
70 nowpacSolver.set_option("gamma",
71 probDescDB.get_real("method.trust_region.contraction_factor") );
72 // Reduction factors:
73 //nowpacSolver.set_option("omega" , 0.8 );
74 //nowpacSolver.set_option("theta" , 0.8 );
75 // Inner boundary path constant: (default should be good and will be adapted)
76 //if (nowpacEvaluator.num_ineq_constraints() > 0)
77 // nowpacSolver.set_option("eps_b" , 1e1 );
78
79 // NOWPAC output verbosity is 0 (least) to 3 (most)
80 nowpacSolver.set_option("verbose", std::min(outputLevel, (short)3) );
81
82 //nowpacSolver.set_max_trustregion( 1e0 ); // scale dependent
83
84 // This option activates SNOWPAC which turns on 3 features:
85 // (a) links TR size to noise returned from evaluator
86 // (b) feasibility restoration (on but not active in deterministic mode)
87 // (c) outer Gaussian process approximation (smooths noisy evaluations)
88 bool stochastic = (methodName == MIT_SNOWPAC);
89 nowpacSolver.set_option("stochastic_optimization", stochastic);
90 // This is tied to the other BlackBoxBaseClass::evaluate() fn redefinition.
91 if (stochastic) {
92 // SNOWPAC picks random points in the trust region to improve the
93 // distribution of the Gaussian Process regression
94 int random_seed = probDescDB.get_int("method.random_seed");
95 if (random_seed) // default for no user spec is zero
96 nowpacSolver.set_option("seed", random_seed);
97 //else SNOWPAC uses a machine generated seed and is non-repeatable
98 const int adaptionSteps = 5*numContinuousVars;
99 nowpacSolver.set_option("GP_adaption_factor", adaptionSteps );
100 nowpacSolver.set_option("use_analytic_smoothing", false);
101 }
102
103 // Maximum number of total accepted steps
104 nowpacSolver.set_option("max_nb_accepted_steps", maxIterations);// inf default
105 // Within special context of meta-iteration like MG/Opt, ensure that we
106 // have at least 2 successful steps
107 //if (subIteratorFlag && ...)
108 // nowpacSolver.set_option("max_nb_accepted_steps", 2 );
109
110 // This is also a termination criterion, but not one of the required ones.
111 // NOWPAC stops if the Frobenius norm of the Hessian blows up, indicating
112 // that the optimizer is now in the noise.
113 //nowpacSolver.set_option("noise_termination" , false);
114 // This is not tied to the other BlackBoxBaseClass::evaluate() function
115 // redefinition, but rather is tied to the Hessian of the computed local
116 // surrogate. Generally, this is only advisable in the deterministic
117 // optimization case. There are additional parameters associated with
118 // termination from noise-detection (e.g., consecutive iterations...)
119
120 // Required:
121 // Must specify one or more stopping criteria: min TR size or maxFnEvals
122
123 nowpacSolver.set_max_number_evaluations(maxFunctionEvals); // default is +inf
124
125 // NOWPAC trust region controls are not relative to global bounds; rather,
126 // they are absolute values for a hyper-sphere of constant dimensional radius.
127 // Therefore, we present a scaled problem to NOWPAC as consistent with these
128 // trust region controls (see use of {un,}scale() within this file).
129 const RealVector& tr_init
130 = probDescDB.get_rv("method.trust_region.initial_size");
131 size_t num_factors = tr_init.length();
132 Real min_factor = probDescDB.get_real("method.trust_region.minimum_size"),
133 tr_factor = (num_factors) ? tr_init[0] : 0.5;
134 if (num_factors > 1)
135 Cerr << "\nWarning: ignoring trailing trust_region initial_size content "
136 << "for NOWPACOptimizer.\n" << std::endl;
137 // domain is [0,1]; default max radius is 1; {init,min}_radius are required
138 // NOWPAC inputs (no defaults) --> use Dakota defaults for {init,min}_radius
139 if (min_factor < 0.) min_factor = 0.; // Dakota default is 1.e-6
140 if (tr_factor < min_factor) tr_factor = min_factor;
141 else if (tr_factor > 1.) tr_factor = 1.;
142 nowpacSolver.set_trustregion(tr_factor, min_factor);
143 //nowpacSolver.set_trustregion(tr_factor); // if only initial
144
145 // Scale the design variables since the TR size controls are absolute, not
146 // relative. Based on the default max TR size of 1., scale to [0,1].
147 RealArray l_bnds, u_bnds; size_t num_v = iteratedModel.cv();
148 l_bnds.assign(num_v, -1.); u_bnds.assign(num_v, 1.);
149 nowpacSolver.set_lower_bounds(l_bnds);
150 nowpacSolver.set_upper_bounds(u_bnds);
151
152 // NOTES from 7/29/15 discussion:
153 // For Lagrangian minimization within MG/Opt:
154 // Option 1. Minimize L s.t. dummy constraints on objective + all constraints
155 // --> NOWPAC will still carry them along, even though inactive, and allow
156 // retrieval of optimal values and derivatives.
157 // --> This will require additional logic in the BBEvaluator in terms of
158 // alternate subproblem formulation
159 // --> allow infeasibility in meta-algorithm
160 // Option 2. If meta-algorithm should enforce feasibility. can change
161 // --> sub-problem to min L s.t. g <= 0. No need for dummy constraints
162 // assuming we back out grad f from grad L, grad g, lambda.
163 }
164
165
166 //void NOWPACOptimizer::initialize_run()
167 //{
168 // Optimizer::initialize_run();
169 //}
170
171
core_run()172 void NOWPACOptimizer::core_run()
173 {
174 //////////////////////////////////////////////////////////////////////////
175 // Set bound constraints at run time to catch late updates
176 nowpacEvaluator.set_unscaled_bounds(iteratedModel.continuous_lower_bounds(),
177 iteratedModel.continuous_upper_bounds());
178
179 // allocate arrays passed to optimization solver
180 RealArray x_star; RealArray obj_star;
181 nowpacEvaluator.scale(iteratedModel.continuous_variables(), x_star);
182 // create data object for nowpac output ( required for warm start )
183 BlackBoxData bb_data(numFunctions, numContinuousVars);
184
185 //////////////////////////////////////////////////////////////////////////
186 // start optimization (on output: bbdata contains data that allows warmstart
187 // and enables post-processing to get model values, gradients and hessians)
188 nowpacSolver.optimize(x_star, obj_star, bb_data);
189 if (outputLevel >= DEBUG_OUTPUT) {
190 // create post-processing object to compute surrogate models
191 PostProcessModels<> PPD( bb_data );
192 Cout << "\n----------------------------------------"
193 << "\nSolution returned from nowpacSolver:\n optimal value = "
194 << obj_star[0] << "\n optimal point =\n" << x_star
195 << "\nData from PostProcessModels:\n"
196 << " tr size = " << PPD.get_trustregion() << '\n';
197 // model value = c + g'(x-x_c) + (x-x_c)'H(x-x_c) / 2
198 // model gradient = g + H (x-x_c)
199 // model Hessian = H
200 for ( int i = 0; i < numFunctions; ++i)
201 Cout << "\n model number " << i+1
202 << "\n value = " << PPD.get_c(i, x_star)
203 << "\n gradient = [\n" << PPD.get_g(i, x_star)
204 << " ]\n hessian = [\n" << PPD.get_H(i) << " ]\n";
205 Cout << "----------------------------------------" << std::endl;
206 }
207
208 //////////////////////////////////////////////////////////////////////////
209 // Publish optimal variables
210 RealVector c_vars = bestVariablesArray.front().continuous_variables_view();
211 nowpacEvaluator.unscale(x_star, c_vars);
212 // Publish optimal response
213 if (!localObjectiveRecast) {
214 RealVector best_fns(numFunctions);
215 const BoolDeque& max_sense = iteratedModel.primary_response_fn_sense();
216 best_fns[0] = (!max_sense.empty() && max_sense[0]) ? -obj_star[0] : obj_star[0];
217
218 // objective and mapped nonlinear inequalities returned from optimize()
219 const SizetList& nln_ineq_map_indices
220 = nowpacEvaluator.nonlinear_inequality_mapping_indices();
221 const RealList& nln_ineq_map_mult
222 = nowpacEvaluator.nonlinear_inequality_mapping_multipliers();
223 const RealList& nln_ineq_map_offsets
224 = nowpacEvaluator.nonlinear_inequality_mapping_offsets();
225 StLCIter i_iter; RLCIter m_iter, o_iter;
226 size_t cntr = 0;//numEqConstraints;
227 for (i_iter = nln_ineq_map_indices.begin(),
228 m_iter = nln_ineq_map_mult.begin(),
229 o_iter = nln_ineq_map_offsets.begin();
230 i_iter != nln_ineq_map_indices.end(); ++i_iter, ++m_iter, ++o_iter)
231 best_fns[(*i_iter)+1] = (obj_star[++cntr] - (*o_iter))/(*m_iter);
232 /*
233 size_t i, offset = iteratedModel.num_nonlinear_ineq_constraints() + 1,
234 num_nln_eq = iteratedModel.num_nonlinear_eq_constraints();
235 const RealVector& nln_eq_targets
236 = iteratedModel.nonlinear_eq_constraint_targets();
237 for (i=0; i<num_nln_eq; i++)
238 best_fns[i+offset] = fn_star[++cntr] + nln_eq_targets[i];
239 */
240
241 bestResponseArray.front().function_values(best_fns);
242 }
243 // else local_objective_recast_retrieve() used in Optimizer::post_run()
244 }
245
246
allocate_constraints()247 void NOWPACBlackBoxEvaluator::allocate_constraints()
248 {
249 // NOWPAC handles 1-sided inequalities <= 0. Equalities cannot be mapped to
250 // two oppositely-signed inequalities due to the interior path requirement.
251 // Hard error for now...
252 bool constraint_err = false;
253 if (iteratedModel.num_nonlinear_eq_constraints()) {
254 Cerr << "Error: NOWPAC does not support nonlinear equality constraints."
255 << std::endl;
256 constraint_err = true;
257 }
258 if (iteratedModel.num_linear_eq_constraints()) {
259 Cerr << "Error: NOWPAC does not support linear equality constraints."
260 << std::endl;
261 constraint_err = true;
262 }
263 if (constraint_err)
264 abort_handler(METHOD_ERROR);
265
266 nonlinIneqConMappingIndices.clear();
267 nonlinIneqConMappingMultipliers.clear();
268 nonlinIneqConMappingOffsets.clear();
269
270 linIneqConMappingIndices.clear();
271 linIneqConMappingMultipliers.clear();
272 linIneqConMappingOffsets.clear();
273
274 // Compute number of 1-sided inequalities to pass to NOWPAC and the mappings
275 // (indices, multipliers, offsets) between DAKOTA and NOWPAC constraints.
276 numNowpacIneqConstr = 0;
277 size_t i, num_nln_ineq = iteratedModel.num_nonlinear_ineq_constraints();
278 const RealVector& nln_ineq_lwr_bnds
279 = iteratedModel.nonlinear_ineq_constraint_lower_bounds();
280 const RealVector& nln_ineq_upr_bnds
281 = iteratedModel.nonlinear_ineq_constraint_upper_bounds();
282 for (i=0; i<num_nln_ineq; i++) {
283 if (nln_ineq_lwr_bnds[i] > -BIG_REAL_BOUND) {
284 ++numNowpacIneqConstr;
285 // nln_ineq_lower_bnd - dakota_constraint <= 0
286 nonlinIneqConMappingIndices.push_back(i);
287 nonlinIneqConMappingMultipliers.push_back(-1.);
288 nonlinIneqConMappingOffsets.push_back(nln_ineq_lwr_bnds[i]);
289 }
290 if (nln_ineq_upr_bnds[i] < BIG_REAL_BOUND) {
291 ++numNowpacIneqConstr;
292 // dakota_constraint - nln_ineq_upper_bnd <= 0
293 nonlinIneqConMappingIndices.push_back(i);
294 nonlinIneqConMappingMultipliers.push_back(1.);
295 nonlinIneqConMappingOffsets.push_back(-nln_ineq_upr_bnds[i]);
296 }
297 }
298 size_t num_lin_ineq = iteratedModel.num_linear_ineq_constraints();
299 const RealVector& lin_ineq_lwr_bnds
300 = iteratedModel.linear_ineq_constraint_lower_bounds();
301 const RealVector& lin_ineq_upr_bnds
302 = iteratedModel.linear_ineq_constraint_upper_bounds();
303 for (i=0; i<num_lin_ineq; i++) {
304 if (lin_ineq_lwr_bnds[i] > -BIG_REAL_BOUND) {
305 ++numNowpacIneqConstr;
306 // lin_ineq_lower_bnd - Ax <= 0
307 linIneqConMappingIndices.push_back(i);
308 linIneqConMappingMultipliers.push_back(-1.);
309 linIneqConMappingOffsets.push_back(lin_ineq_lwr_bnds[i]);
310 }
311 if (lin_ineq_upr_bnds[i] < BIG_REAL_BOUND) {
312 ++numNowpacIneqConstr;
313 // Ax - lin_ineq_upper_bnd <= 0
314 linIneqConMappingIndices.push_back(i);
315 linIneqConMappingMultipliers.push_back(1.);
316 linIneqConMappingOffsets.push_back(-lin_ineq_upr_bnds[i]);
317 }
318 }
319 }
320
321
322 void NOWPACBlackBoxEvaluator::
evaluate(RealArray const & x,RealArray & vals,void * param)323 evaluate(RealArray const &x, RealArray &vals, void *param)
324 {
325 // NOWPACOptimizer enforces an embedded scaling: incoming x is scaled on [0,1]
326 // --> unscale for posting to iteratedModel
327 RealVector& c_vars
328 = iteratedModel.current_variables().continuous_variables_view();
329 unscale(x, c_vars);
330
331 iteratedModel.evaluate(); // no ASV control, use default
332
333 const RealVector& dakota_fns
334 = iteratedModel.current_response().function_values();
335 // If no mappings...
336 //copy_data(dakota_fns, vals);
337
338 // apply optimization sense mapping. Note: Any MOO/NLS recasting is
339 // responsible for setting the scalar min/max sense within the recast.
340 const BoolDeque& max_sense = iteratedModel.primary_response_fn_sense();
341 Real obj_fn = dakota_fns[0];
342 vals[0] = (!max_sense.empty() && max_sense[0]) ? -obj_fn : obj_fn;
343
344 // apply nonlinear inequality constraint mappings
345 StLIter i_iter; RLIter m_iter, o_iter; size_t cntr = 0;
346 for (i_iter = nonlinIneqConMappingIndices.begin(),
347 m_iter = nonlinIneqConMappingMultipliers.begin(),
348 o_iter = nonlinIneqConMappingOffsets.begin();
349 i_iter != nonlinIneqConMappingIndices.end();
350 ++i_iter, ++m_iter, ++o_iter) // nonlinear ineq
351 vals[++cntr] = (*o_iter) + (*m_iter) * dakota_fns[(*i_iter)+1];
352
353 // apply linear inequality constraint mappings
354 const RealMatrix& lin_ineq_coeffs
355 = iteratedModel.linear_ineq_constraint_coeffs();
356 size_t j, num_cv = x.size();
357 for (i_iter = linIneqConMappingIndices.begin(),
358 m_iter = linIneqConMappingMultipliers.begin(),
359 o_iter = linIneqConMappingOffsets.begin();
360 i_iter != linIneqConMappingIndices.end();
361 ++i_iter, ++m_iter, ++o_iter) { // linear ineq
362 size_t index = *i_iter;
363 Real Ax = 0.;
364 for (j=0; j<num_cv; ++j)
365 Ax += lin_ineq_coeffs(index,j) * x[j];
366 vals[++cntr] = (*o_iter) + (*m_iter) * Ax;
367 }
368 }
369 // TO DO: asynchronous evaluate_nowait()/synchronize()
370
371
372 void NOWPACBlackBoxEvaluator::
evaluate(RealArray const & x,RealArray & vals,RealArray & noise,void * param)373 evaluate(RealArray const &x, RealArray &vals, RealArray &noise, void *param)
374 {
375 // NOWPACOptimizer enforces an embedded scaling: incoming x is scaled on [0,1]
376 // --> unscale for posting to iteratedModel
377 RealVector& c_vars
378 = iteratedModel.current_variables().continuous_variables_view();
379 unscale(x, c_vars);
380
381 iteratedModel.evaluate(); // no ASV control, use default
382
383 const RealVector& dakota_fns
384 = iteratedModel.current_response().function_values();
385 // NonD implements std error estimates for selected QoI statistics (see, e.g.,
386 // Harting et al. for MC error estimates). NestedModel implements
387 // sub-iterator mappings for both fn vals & std errors.
388 const RealVector& errors = iteratedModel.error_estimates();
389
390 // apply optimization sense mapping. Note: Any MOO/NLS recasting is
391 // responsible for setting the scalar min/max sense within the recast.
392 const BoolDeque& max_sense = iteratedModel.primary_response_fn_sense();
393 Real obj_fn = dakota_fns[0], mult;
394 size_t index, cntr = 0;
395 vals[cntr] = (!max_sense.empty() && max_sense[0]) ? -obj_fn : obj_fn;
396 noise[cntr] = 2. * errors[0]; // for now; TO DO: mapping of noise for MOO/NLS...
397 ++cntr;
398
399 // apply nonlinear inequality constraint mappings
400 StLIter i_iter; RLIter m_iter, o_iter;
401 for (i_iter = nonlinIneqConMappingIndices.begin(),
402 m_iter = nonlinIneqConMappingMultipliers.begin(),
403 o_iter = nonlinIneqConMappingOffsets.begin();
404 i_iter != nonlinIneqConMappingIndices.end();
405 ++i_iter, ++m_iter, ++o_iter, ++cntr) { // nonlinear ineq
406 index = (*i_iter)+1; // offset single objective
407 mult = (*m_iter);
408 vals[cntr] = (*o_iter) + mult * dakota_fns[index];
409 noise[cntr] = 2. * std::abs(mult) * errors[index];
410 }
411
412 // apply linear inequality constraint mappings
413 const RealMatrix& lin_ineq_coeffs
414 = iteratedModel.linear_ineq_constraint_coeffs();
415 size_t j, num_cv = x.size();
416 for (i_iter = linIneqConMappingIndices.begin(),
417 m_iter = linIneqConMappingMultipliers.begin(),
418 o_iter = linIneqConMappingOffsets.begin();
419 i_iter != linIneqConMappingIndices.end();
420 ++i_iter, ++m_iter, ++o_iter, ++cntr) { // linear ineq
421 size_t index = *i_iter;
422 Real Ax = 0.;
423 for (j=0; j<num_cv; ++j)
424 Ax += lin_ineq_coeffs(index,j) * x[j];
425 vals[cntr] = (*o_iter) + (*m_iter) * Ax;
426 noise[cntr] = 0.; // no error in linear case
427 }
428 }
429
evaluate_samples(std::vector<double> const & samples,const unsigned int index,std::vector<double> const & x)430 double NOWPACBlackBoxEvaluator::evaluate_samples ( std::vector<double> const &samples, const unsigned int index, std::vector<double> const &x){
431 //TO DO implement if smoothing works at some point.
432 return -1;
433 }
434 // TO DO: asynchronous evaluate_nowait()/synchronize()
435
436
437 #ifdef HAVE_DYNLIB_FACTORIES
new_NOWPACOptimizer(ProblemDescDB & problem_db,Model & model)438 NOWPACOptimizer* new_NOWPACOptimizer(ProblemDescDB& problem_db, Model& model)
439 {
440 #ifdef DAKOTA_DYNLIB
441 not_available("NOWPAC");
442 return 0;
443 #else
444 return new NOWPACOptimizer(problem_db, model);
445 #endif // DAKOTA_DYNLIB
446 }
447
new_NOWPACOptimizer(Model & model)448 NOWPACOptimizer* new_NOWPACOptimizer(Model& model)
449 {
450 #ifdef DAKOTA_DYNLIB
451 not_available("NOWPAC");
452 return 0;
453 #else
454 return new NOWPACOptimizer(model);
455 #endif // DAKOTA_DYNLIB
456 }
457 #endif // HAVE_DYNLIB_FACTORIES
458
459 } // namespace Dakota
460