diff --git a/FieldOpt/Optimization/constraints/constraint_handler.cpp b/FieldOpt/Optimization/constraints/constraint_handler.cpp index f1f4cde4b..581543de7 100644 --- a/FieldOpt/Optimization/constraints/constraint_handler.cpp +++ b/FieldOpt/Optimization/constraints/constraint_handler.cpp @@ -219,8 +219,16 @@ long double ConstraintHandler::GetWeightedNormalizedPenalties(Case *c) { long double wnp = 0.0L; for (auto con : constraints_) { long double pen = con->PenaltyNormalized(c); + if (VERB_OPT >= 3 && pen > 0) { + Printer::ext_info("Penalty from constraint " + con->name() + + ": " + Printer::num2str(pen), "Optimization", "ConstraintHandler"); + } wnp += pen * con->GetPenaltyWeight(); } + if (VERB_OPT >= 2) { + Printer::ext_info("Weighted, normalized penalty for case " + c->id().toString().toStdString() + + ": " + Printer::num2str(wnp), "Optimization", "ConstraintHandler"); + } return wnp; } diff --git a/FieldOpt/Optimization/constraints/interwell_distance.cpp b/FieldOpt/Optimization/constraints/interwell_distance.cpp index 5c1e22df0..d88cd519a 100644 --- a/FieldOpt/Optimization/constraints/interwell_distance.cpp +++ b/FieldOpt/Optimization/constraints/interwell_distance.cpp @@ -21,6 +21,8 @@ #include "ConstraintMath/well_constraint_projections/well_constraint_projections.h" #include #include +#include "Utilities/verbosity.h" +#include "Utilities/printer.hpp" namespace Optimization { namespace Constraints { @@ -111,9 +113,9 @@ void InterwellDistance::InitializeNormalizer(QList cases) { long double minimum_distance = 1e20; for (auto c : cases) { vector endp_dist = endpointDistances(c); - double min_dist = *min_element(endp_dist.begin(), endp_dist.end()); - if (min_dist < minimum_distance){ - minimum_distance = min_dist; + for (double dist : endp_dist) { + if (abs(dist) < minimum_distance) + minimum_distance = abs(dist); } } normalizer_.set_max(1.0L); @@ -125,9 +127,13 @@ double InterwellDistance::Penalty(Case *c) { double violation = 0.0; for (auto distance : endpoint_distances) { if (distance < distance_) { + if (VERB_OPT >= 2) Printer::ext_info("Interwell distance penalty for case " + c->id().toString().toStdString() + + ": " + Printer::num2str(distance - distance_) + " m too close"); violation += abs(distance - distance_); } } + if (VERB_OPT >= 2) Printer::ext_info("Interwell distance total violation for case " + + c->id().toString().toStdString() + ": " + Printer::num2str(violation)); return violation; } @@ -148,7 +154,11 @@ vector InterwellDistance::endpointDistances(Case *c) { return endpoint_distances; } long double InterwellDistance::PenaltyNormalized(Case *c) { - return normalizer_.normalize(Penalty(c)); + double penalty = Penalty(c); + if (penalty > 0.0) + return normalizer_.normalize(penalty); + else + return 0.0; } } diff --git a/FieldOpt/Optimization/constraints/reservoir_boundary.cpp b/FieldOpt/Optimization/constraints/reservoir_boundary.cpp index f59e3bfd1..85bc86ea3 100644 --- a/FieldOpt/Optimization/constraints/reservoir_boundary.cpp +++ b/FieldOpt/Optimization/constraints/reservoir_boundary.cpp @@ -20,6 +20,8 @@ #include "ConstraintMath/well_constraint_projections/well_constraint_projections.h" #include #include "Utilities/math.hpp" +#include "Utilities/verbosity.h" +#include "Utilities/printer.hpp" namespace Optimization { namespace Constraints { @@ -455,14 +457,20 @@ void ReservoirBoundary::InitializeNormalizer(QList cases) { double ReservoirBoundary::Penalty(Case *c) { if (CaseSatisfiesConstraint(c)) return 0.0; - else + else { + if (VERB_OPT >= 2) Printer::ext_info("Reservoir boundary penalty (infty) for case " + c->id().toString().toStdString(), + "Optimization", "ReservoirBoundary"); return INFINITY; + } } long double ReservoirBoundary::PenaltyNormalized(Case *c) { if (CaseSatisfiesConstraint(c)) return 0.0L; - else + else { + if (VERB_OPT >= 2) Printer::ext_info("Reservoir boundary penalty (infty) for case " + c->id().toString().toStdString(), + "Optimization", "ReservoirBoundary"); return 1.0L; + } } } diff --git a/FieldOpt/Optimization/constraints/well_spline_length.cpp b/FieldOpt/Optimization/constraints/well_spline_length.cpp index cf7f070d3..7d2612117 100644 --- a/FieldOpt/Optimization/constraints/well_spline_length.cpp +++ b/FieldOpt/Optimization/constraints/well_spline_length.cpp @@ -18,6 +18,8 @@ ******************************************************************************/ #include "Utilities/math.hpp" +#include "Utilities/verbosity.h" +#include "Utilities/printer.hpp" #include "well_spline_length.h" #include "ConstraintMath/well_constraint_projections/well_constraint_projections.h" @@ -103,10 +105,16 @@ void WellSplineLength::InitializeNormalizer(QList cases) { double WellSplineLength::Penalty(Case *c) { auto endpts = GetEndpointValueVectors(c, affected_well_); double well_length = (endpts.first - endpts.second).norm(); - if (well_length > max_length_) + if (well_length > max_length_) { + if (VERB_OPT >= 2) Printer::ext_info("Well length penalty for case " + c->id().toString().toStdString() + + ": " + Printer::num2str(well_length - max_length_) + " m too long", "Optimization", "WellSplineLength"); return well_length - max_length_; - else if (well_length < min_length_) + } + else if (well_length < min_length_) { + if (VERB_OPT >= 2) Printer::ext_info("Well length penalty for case " + c->id().toString().toStdString() + + ":" + Printer::num2str(min_length_ - well_length) + " m too short", "Optimization", "WellSplineLength"); return min_length_ - well_length; + } else return 0.0; } diff --git a/FieldOpt/Optimization/optimizer.cpp b/FieldOpt/Optimization/optimizer.cpp index 9d6223407..0d2b20a78 100644 --- a/FieldOpt/Optimization/optimizer.cpp +++ b/FieldOpt/Optimization/optimizer.cpp @@ -47,6 +47,7 @@ Optimizer::Optimizer(Settings::Optimizer *settings, Case *base_case, else { Printer::ext_info("Using shared CaseHandler.", "Optimizer", "Optimization"); case_handler_ = case_handler; + is_hybrid_component_ = true; } if (constraint_handler == 0) { constraint_handler_ = new Constraints::ConstraintHandler(settings->constraints(), variables, grid); @@ -64,6 +65,25 @@ Optimizer::Optimizer(Settings::Optimizer *settings, Case *base_case, enable_logging_ = true; verbosity_level_ = 0; penalize_ = settings->objective().use_penalty_function; + + if (penalize_) { + if (!normalizer_ofv_.is_ready()) { + if (VERB_OPT >=1) { + Printer::ext_info("Initializing normalizers", "Optimization", "Optimizer"); + } + initializeNormalizers(); + + // penalize the base case + double org_ofv = tentative_best_case_->objective_function_value(); + double pen_ofv = PenalizedOFV(tentative_best_case_); + tentative_best_case_->set_objective_function_value(pen_ofv); + if (VERB_OPT >=1) { + Printer::ext_info("Penalized base case. " + "Original value: " + Printer::num2str(org_ofv) + "; " + + "Penalized value: " + Printer::num2str(pen_ofv), "Optimization", "Optimizer"); + } + } + } } Case *Optimizer::GetCaseForEvaluation() @@ -81,7 +101,7 @@ Case *Optimizer::GetCaseForEvaluation() void Optimizer::SubmitEvaluatedCase(Case *c) { evaluated_cases_++; - if (penalize_ && iteration_ > 0) { + if (penalize_) { double penalized_ofv = PenalizedOFV(c); c->set_objective_function_value(penalized_ofv); } @@ -243,13 +263,23 @@ double Optimizer::PenalizedOFV(Case *c) { long double norm_ofv = normalizer_ofv_.normalize(c->objective_function_value()); long double penalty = constraint_handler_->GetWeightedNormalizedPenalties(c); long double norm_pen_ovf = norm_ofv - penalty; + double denormalized_ofv = normalizer_ofv_.denormalize(norm_pen_ovf); + + if (VERB_OPT >= 2) { + Printer::ext_info("Penalized case " + c->id().toString().toStdString() + ". " + "Initial OFV: " + Printer::num2str(c->objective_function_value()) + "; " + "Normalized OFV :" + Printer::num2str(norm_ofv) + "; " + "Normalized penalty: " + Printer::num2str(penalty) + "; " + "Penalized, normalized OFV: " + Printer::num2str(norm_pen_ovf) + "; " + "Denormalized, penalized OFV: " + Printer::num2str(denormalized_ofv), "Optimization", "Optimizer"); + } if (norm_pen_ovf <= 0.0L) { cout << "RETURNING ZERO OFV" << endl; return 0.0; } else { - return normalizer_ofv_.denormalize(norm_pen_ovf); + return denormalized_ofv; } } diff --git a/FieldOpt/Optimization/optimizer.h b/FieldOpt/Optimization/optimizer.h index b887b4208..f78647929 100644 --- a/FieldOpt/Optimization/optimizer.h +++ b/FieldOpt/Optimization/optimizer.h @@ -177,6 +177,7 @@ class Optimizer : public Loggable bool penalize_; //!< Switch for whether or not to use penalty function to account for constraints. Case *tentative_best_case_; //!< The best case encountered thus far. int tentative_best_case_iteration_; //!< The iteration in which the current tentative best case was found. + bool is_hybrid_component_; //!< Indicates that this object is a hybrid optimization component. Normalizer normalizer_ofv_; //!< Normalizer for objective function values. diff --git a/FieldOpt/Optimization/optimizers/APPS.cpp b/FieldOpt/Optimization/optimizers/APPS.cpp index 45ea0ea2c..790a58d34 100644 --- a/FieldOpt/Optimization/optimizers/APPS.cpp +++ b/FieldOpt/Optimization/optimizers/APPS.cpp @@ -44,13 +44,6 @@ APPS::APPS(Settings::Optimizer *settings, assert(settings->parameters().max_queue_size >= 1.0); max_queue_length_ = directions_.size() * settings->parameters().max_queue_size; is_async_ = true; - if (penalize_) { - if (!normalizer_ofv_.is_ready()) initializeNormalizers(); - - // penalize the base case - double pen_ofv = PenalizedOFV(tentative_best_case_); - tentative_best_case_->set_objective_function_value(pen_ofv); - } if (enable_logging_) { logger_->AddEntry(this); } @@ -60,11 +53,16 @@ void APPS::iterate() { if (enable_logging_) { logger_->AddEntry(this); } - if (inactive().size() > 0) { + if (inactive().size() > 0 && evaluated_cases_ < max_evaluations_) { case_handler_->AddNewCases(generate_trial_points(inactive())); set_active(inactive()); } iteration_++; + if (is_hybrid_component_) { + // Increment this here if this object is a hybrid optimization component, + // as it will not be incremented elsewhere. + evaluated_cases_ = iteration_; + } } void APPS::handleEvaluatedCase(Case *c) { @@ -89,6 +87,7 @@ void APPS::unsuccessful_iteration(Case *c) { set_inactive(unsuccessful_direction); contract(unsuccessful_direction); } + if (VERB_OPT >= 3) print_state("Unsuccessful iteration"); if (!is_converged()) iterate(); } @@ -117,10 +116,13 @@ vector APPS::inactive() { } void APPS::prune_queue() { - if (case_handler_->QueuedCases().size() <= max_queue_length_ - directions_.size()) + if (case_handler_->QueuedCases().size() <= max_queue_length_ - directions_.size()) { return; + } else { - while (case_handler_->QueuedCases().size() > max_queue_length_ - directions_.size()) { + int queue_size = max_queue_length_ - directions_.size(); + if (evaluated_cases_ >= max_evaluations_) queue_size = 1; + while (case_handler_->QueuedCases().size() > queue_size) { auto dequeued_case = dequeue_case_with_worst_origin(); if (dequeued_case->origin_case()->id() == GetTentativeBestCase()->id()) set_inactive(vector{dequeued_case->origin_direction_index()}); @@ -132,10 +134,12 @@ void APPS::prune_queue() { void APPS::print_state(string header) { std::stringstream ss; ss << header << "|"; - ss << "Step lengths : " << vec_to_str(vector(step_lengths_.data(), step_lengths_.data() + step_lengths_.size())) << "|"; - ss << "Iteration: " << iteration_ << "|"; + ss << "Iteration: " << iteration_ << "|"; + ss << "Evaluated cases: " << evaluated_cases_ << "|"; + ss << "Queued cases: " << case_handler_->QueuedCases().size() << "|"; ss << "Current best case: " << tentative_best_case_->id().toString().toStdString() << "|"; - ss << " OFV: " << tentative_best_case_->objective_function_value(); + ss << "OFV: " << tentative_best_case_->objective_function_value(); + ss << "Step lengths : " << vec_to_str(vector(step_lengths_.data(), step_lengths_.data() + step_lengths_.size())) << "|"; Printer::ext_info(ss.str(), "Optimization", "APPS"); } } diff --git a/FieldOpt/Optimization/optimizers/GeneticAlgorithm.cpp b/FieldOpt/Optimization/optimizers/GeneticAlgorithm.cpp index 233497990..b177479ff 100644 --- a/FieldOpt/Optimization/optimizers/GeneticAlgorithm.cpp +++ b/FieldOpt/Optimization/optimizers/GeneticAlgorithm.cpp @@ -131,12 +131,6 @@ Case *GeneticAlgorithm::generateRandomCase() { new_case->SetRealVarValues(erands); return new_case; } -void GeneticAlgorithm::penalizeInitialGeneration() { - initializeNormalizers(); - for (auto chrom : population_) { - double pen_ofv = PenalizedOFV(chrom.case_pointer); - chrom.case_pointer->set_objective_function_value(pen_ofv); - } -} + } } diff --git a/FieldOpt/Optimization/optimizers/GeneticAlgorithm.h b/FieldOpt/Optimization/optimizers/GeneticAlgorithm.h index 8a09707ef..ccd4bbb39 100644 --- a/FieldOpt/Optimization/optimizers/GeneticAlgorithm.h +++ b/FieldOpt/Optimization/optimizers/GeneticAlgorithm.h @@ -120,11 +120,6 @@ class GeneticAlgorithm : public Optimizer { */ Case *generateRandomCase(); - /*! - * @brief Initialize normalizers and penalize the initial generation after - * the fact (should be start of the first iteration). - */ - void penalizeInitialGeneration(); }; } diff --git a/FieldOpt/Optimization/optimizers/RGARDD.cpp b/FieldOpt/Optimization/optimizers/RGARDD.cpp index 2633bedc9..f8e67e606 100644 --- a/FieldOpt/Optimization/optimizers/RGARDD.cpp +++ b/FieldOpt/Optimization/optimizers/RGARDD.cpp @@ -39,6 +39,7 @@ RGARDD::RGARDD(Settings::Optimizer *settings, stagnation_limit_ = settings->parameters().stagnation_limit; mating_pool_ = population_; if (enable_logging_) { + logger_->AddEntry(this); logger_->AddEntry(new ConfigurationSummary(this)); } } @@ -51,9 +52,6 @@ void RGARDD::iterate() { if (enable_logging_) { logger_->AddEntry(this); } - if (iteration_ == 0 && penalize_) { // If we're done evaluating the initial population ... - penalizeInitialGeneration(); - } population_ = sortPopulation(population_); if (is_stagnant()) { @@ -110,6 +108,9 @@ void RGARDD::handleEvaluatedCase(Case *c) { population_[index] = mating_pool_[index]; if (isImprovement(c)) { updateTentativeBestCase(c); + if (enable_logging_) { + logger_->AddEntry(this); + } if (VERB_OPT >= 1) { Printer::ext_info("New best in generation " + Printer::num2str(iteration_) + ": " + Printer::num2str(GetTentativeBestCase()->objective_function_value()), "RGARDD", "Optimization"); diff --git a/FieldOpt/Optimization/optimizers/bayesian_optimization/EGO.cpp b/FieldOpt/Optimization/optimizers/bayesian_optimization/EGO.cpp index 973123d17..eec2eaeda 100644 --- a/FieldOpt/Optimization/optimizers/bayesian_optimization/EGO.cpp +++ b/FieldOpt/Optimization/optimizers/bayesian_optimization/EGO.cpp @@ -44,6 +44,20 @@ EGO::EGO(Settings::Optimizer *settings, time_af_opt_ = 0; settings_ = settings; + initializeNormalizers(); + + // penalize the base case + if (penalize_) { + double org_ofv = tentative_best_case_->objective_function_value(); + double pen_ofv = PenalizedOFV(tentative_best_case_); + tentative_best_case_->set_objective_function_value(pen_ofv); + if (VERB_OPT >=1) { + Printer::ext_info("Penalized base case. " + "Original value: " + Printer::num2str(org_ofv) + "; " + + "Penalized value: " + Printer::num2str(pen_ofv), "Optimization", "Optimizer"); + } + } + int n_cont_vars = variables->ContinousVariableSize(); if (constraint_handler_->HasBoundaryConstraints()) { @@ -149,9 +163,6 @@ Optimization::Optimizer::TerminationCondition EGO::IsFinished() { return tc; } void EGO::handleEvaluatedCase(Case *c) { - if (!normalizer_ofv_.is_ready()) - initializeNormalizers(); - gp_->add_pattern(c->GetRealVarVector().data(), normalizer_ofv_.normalize(c->objective_function_value())); if (isImprovement(c)) { updateTentativeBestCase(c); @@ -159,12 +170,6 @@ void EGO::handleEvaluatedCase(Case *c) { } } void EGO::iterate() { - if (!normalizer_ofv_.is_ready()) { - initializeNormalizers(); - // Add base case to GP - gp_->add_pattern(tentative_best_case_->GetRealVarVector().data(), normalizer_ofv_.normalize(tentative_best_case_->objective_function_value())); - } - if (enable_logging_) { logger_->AddEntry(this); }