1 | // Copyright (C) 2004, International Business Machines and others. |
---|
2 | // All Rights Reserved. |
---|
3 | // This code is published under the Common Public License. |
---|
4 | // |
---|
5 | // $Id: IpFilterLineSearch.cpp 412 2005-07-26 20:42:51Z andreasw $ |
---|
6 | // |
---|
7 | // Authors: Carl Laird, Andreas Waechter IBM 2004-08-13 |
---|
8 | |
---|
9 | #include "IpFilterLineSearch.hpp" |
---|
10 | #include "IpJournalist.hpp" |
---|
11 | #include "IpRestoPhase.hpp" |
---|
12 | #include "IpAlgTypes.hpp" |
---|
13 | |
---|
14 | #ifdef OLD_C_HEADERS |
---|
15 | # include <math.h> |
---|
16 | # include <ctype.h> |
---|
17 | #else |
---|
18 | # include <cmath> |
---|
19 | # include <cctype> |
---|
20 | #endif |
---|
21 | |
---|
22 | namespace Ipopt |
---|
23 | { |
---|
24 | |
---|
25 | DBG_SET_VERBOSITY(0); |
---|
26 | |
---|
27 | DefineIpoptType(FilterLineSearch); |
---|
28 | |
---|
29 | FilterLineSearch::FilterLineSearch(const SmartPtr<RestorationPhase>& resto_phase, |
---|
30 | const SmartPtr<PDSystemSolver>& pd_solver) |
---|
31 | : |
---|
32 | LineSearch(), |
---|
33 | resto_phase_(resto_phase), |
---|
34 | pd_solver_(pd_solver), |
---|
35 | theta_min_(-1.0), |
---|
36 | theta_max_(-1.0), |
---|
37 | filter_(2) |
---|
38 | { |
---|
39 | DBG_START_FUN("FilterLineSearch::FilterLineSearch", |
---|
40 | dbg_verbosity); |
---|
41 | } |
---|
42 | |
---|
43 | FilterLineSearch::~FilterLineSearch() |
---|
44 | { |
---|
45 | DBG_START_FUN("FilterLineSearch::~FilterLineSearch()", |
---|
46 | dbg_verbosity); |
---|
47 | } |
---|
48 | |
---|
49 | void FilterLineSearch::RegisterOptions(SmartPtr<RegisteredOptions> roptions) |
---|
50 | { |
---|
51 | roptions->AddLowerBoundedNumberOption( |
---|
52 | "theta_max_fact", |
---|
53 | "Determines upper bound for constraint violation in the filter.", |
---|
54 | 0.0, true, 1e4, |
---|
55 | "The algorithmic parameter theta_max is determined as theta_max_fact " |
---|
56 | "times the maximum of 1 and the constraint violation at initial point. " |
---|
57 | "Any point with a constraint violation larger than theta_max is " |
---|
58 | "unacceptable to the filter (see Eqn. (21) in implementation paper)."); |
---|
59 | roptions->AddLowerBoundedNumberOption( |
---|
60 | "theta_min_fact", |
---|
61 | "Determines constraint violation threshold in switching rule.", |
---|
62 | 0.0, true, 1e-4, |
---|
63 | "The algorithmic parameter theta_min is determined as theta_max_fact " |
---|
64 | "times the maximum of 1 and the constraint violation at initial point. " |
---|
65 | "The switching rules treats an iteration as h-type iteration whenever " |
---|
66 | "the current constraint violation is larger than theta_min (see " |
---|
67 | "paragraph before Eqn. (19) in implementation paper)."); |
---|
68 | roptions->AddBoundedNumberOption( |
---|
69 | "eta_phi", |
---|
70 | "Relaxation factor in the Armijo condition.", |
---|
71 | 0.0, true, 0.5, true, 1e-8, |
---|
72 | "(See Eqn. (20) in implementation paper)"); |
---|
73 | roptions->AddLowerBoundedNumberOption( |
---|
74 | "delta", "Multiplier for constraint violation in switching rule.", |
---|
75 | 0.0, true, 1.0, |
---|
76 | "(See Eqn. (19) in implementation paper)"); |
---|
77 | roptions->AddLowerBoundedNumberOption( |
---|
78 | "s_phi", |
---|
79 | "Exponent for linear barrier function model in switching rule.", |
---|
80 | 1.0, true, 2.3, |
---|
81 | "(See Eqn. (19) in implementation paper)"); |
---|
82 | roptions->AddLowerBoundedNumberOption( |
---|
83 | "s_theta", |
---|
84 | "Exponent for current constraint violation in switching rule.", |
---|
85 | 1.0, true, 1.1, |
---|
86 | "(See Eqn. (19) in implementation paper)"); |
---|
87 | roptions->AddBoundedNumberOption( |
---|
88 | "gamma_phi", |
---|
89 | "Relaxation factor in filter margin for barrier function.", |
---|
90 | 0.0, true, 1.0, true, 1e-8, |
---|
91 | "(See Eqn. (18a) in implementation paper)"); |
---|
92 | roptions->AddBoundedNumberOption( |
---|
93 | "gamma_theta", |
---|
94 | "Relaxation factor in filter margin for constraint violation.", |
---|
95 | 0.0, true, 1.0, true, 1e-5, |
---|
96 | "(See Eqn. (18b) in implementation paper)"); |
---|
97 | roptions->AddBoundedNumberOption( |
---|
98 | "alpha_min_frac", |
---|
99 | "Safety factor for minimal step size (switch to restoration phase).", |
---|
100 | 0.0, true, 1.0, true, 0.05, |
---|
101 | "(This is gamma_alpha in Eqn. (20) in implementation paper)"); |
---|
102 | roptions->AddBoundedNumberOption( |
---|
103 | "alpha_red_factor", |
---|
104 | "Fractional reduction of trial step size in the backtracking line search.", |
---|
105 | 0.0, true, 1.0, true, 0.5, |
---|
106 | "Determines the fraction by how much the trial step size is reduced in " |
---|
107 | "every step of the backtracking line search."); |
---|
108 | roptions->AddLowerBoundedIntegerOption( |
---|
109 | "max_soc", |
---|
110 | "Maximal number of second order correction trial steps.", |
---|
111 | 0, 4, |
---|
112 | "Determines the maximal number of second order correction trial steps " |
---|
113 | "that should be performed. Choosing 0 disables the second order " |
---|
114 | "corrections. (This is p^{max} of Step A-5.9 of " |
---|
115 | "Algorithm A in implementation paper.)"); |
---|
116 | roptions->AddLowerBoundedNumberOption( |
---|
117 | "kappa_soc", |
---|
118 | "Factor in sufficient reduction rule for second order correction.", |
---|
119 | 0.0, true, 0.99, |
---|
120 | "Determines by how much a second order correction step must reduce the " |
---|
121 | "constraint violation so that further correction steps are attempted. " |
---|
122 | "(See Step A-5.9 of Algorithm A in implementation paper.)"); |
---|
123 | roptions->AddLowerBoundedNumberOption( |
---|
124 | "obj_max_inc", |
---|
125 | "Determines upper bound on acceptable increase of barrier objective function.", |
---|
126 | 1.0, true, 5.0, |
---|
127 | "A trial point leading to more orders of magnitude increase in the " |
---|
128 | "barrier objective function are rejected."); |
---|
129 | roptions->AddStringOption2( |
---|
130 | "magic_steps", |
---|
131 | "Enables magic steps.", |
---|
132 | "no", |
---|
133 | "no", "don't take magic steps", |
---|
134 | "yes", "take magic steps", |
---|
135 | "DOESN'T REALLY WORK YET!"); |
---|
136 | roptions->AddStringOption3( |
---|
137 | "corrector_type", |
---|
138 | "Type of corrector steps.", |
---|
139 | "none", |
---|
140 | "none", "no corrector", |
---|
141 | "affine", "corrector step towards mu=0", |
---|
142 | "primal-dual", "corrector step towards current mu", |
---|
143 | "Determines what kind of corrector steps should be tried."); |
---|
144 | |
---|
145 | roptions->AddStringOption2( |
---|
146 | "skip_corr_if_neg_curv", |
---|
147 | "Skip the corrector step in negative curvature iteration.", |
---|
148 | "yes", |
---|
149 | "no", "don't skip", |
---|
150 | "yes", "skip", |
---|
151 | "The corrector step is not tried if during the computation of " |
---|
152 | "the search direction in the current iteration negative curvature has " |
---|
153 | "been encountered."); |
---|
154 | |
---|
155 | roptions->AddStringOption2( |
---|
156 | "skip_corr_in_monotone_mode", |
---|
157 | "Skip the corrector step during monotone barrier parameter mode.", |
---|
158 | "yes", |
---|
159 | "no", "don't skip", |
---|
160 | "yes", "skip", |
---|
161 | "The corrector step is not tried if the algorithm is currently in the " |
---|
162 | "monotone mode (see also option \"barrier_strategy\")."); |
---|
163 | |
---|
164 | roptions->AddStringOption2( |
---|
165 | "accept_every_trial_step", |
---|
166 | "always accept the frist trial step", |
---|
167 | "no", |
---|
168 | "no", "don't arbitrarily accept the full step", |
---|
169 | "yes", "always accept the full step", |
---|
170 | "Setting this option to \"yes\" essentially disables the line search " |
---|
171 | "and makes the algorithm take aggressive steps."); |
---|
172 | |
---|
173 | roptions->AddStringOption7( |
---|
174 | "alpha_for_y", |
---|
175 | "Step size for constraint multipliers.", |
---|
176 | "primal", |
---|
177 | "primal", "use primal step size", |
---|
178 | "bound_mult", "use step size for the bound multipliers", |
---|
179 | "min", "use the min of primal and bound multipliers", |
---|
180 | "max", "use the max of primal and bound multipliers", |
---|
181 | "full", "take a full step of size one", |
---|
182 | "min_dual_infeas", "choose step size minimizing new dual infeasibility", |
---|
183 | "safe_min_dual_infeas", "like \"min_dual_infeas\", but safeguarded by \"min\" and \"max\"", |
---|
184 | "Determines which step size (alpha_y) should be used to update the " |
---|
185 | "constraint multipliers."); |
---|
186 | |
---|
187 | roptions->AddLowerBoundedNumberOption( |
---|
188 | "corrector_compl_avrg_red_fact", |
---|
189 | "Complementarity tolerance factor for accepting corrector step", |
---|
190 | 0.0, true, 1.0, |
---|
191 | "Determines the factor by which complementarity is allowed to increase " |
---|
192 | "for a corrector step to be accepted."); |
---|
193 | |
---|
194 | roptions->AddStringOption2( |
---|
195 | "expect_infeasible_problem", |
---|
196 | "Enable heuristics to quickly detect an infeasible problem.", |
---|
197 | "no", |
---|
198 | "no", "the problem probably be feasible", |
---|
199 | "yes", "the problem has a good chance to be infeasible", |
---|
200 | "This options is meant to activate heuristics that may speed up the " |
---|
201 | "infeasibility determination if you expect the problem to be " |
---|
202 | "infeasible. In the filter line search procedure, the restoration " |
---|
203 | "phase is called more qucikly than usually, and more reduction in " |
---|
204 | "the constraint violation is enforced. If the problem is square, this " |
---|
205 | "is enabled automatically."); |
---|
206 | roptions->AddLowerBoundedNumberOption( |
---|
207 | "expect_infeasible_problem_ctol", |
---|
208 | "Threshold for disabling \"expect_infeasible_problem\" option", |
---|
209 | 0.0, false, 1e-3, |
---|
210 | "If the constraint violation becomes small than this threshold, " |
---|
211 | "the \"expect_infeasible_problem\" heuristics in the filter line " |
---|
212 | "search will are disabled. If the problem is square, this is set to 0."); |
---|
213 | roptions->AddLowerBoundedNumberOption( |
---|
214 | "soft_resto_pderror_reduction_factor", |
---|
215 | "Required reduction in primal-dual error in soft restoration phase.", |
---|
216 | 0.0, false, (1.0 - 1e-4), |
---|
217 | "For the soft restoration phase (which attempts to reduce the " |
---|
218 | "primal-dual error with regular steps), this indicates by which " |
---|
219 | "factor the primal-dual error has to be reduced in order to continue " |
---|
220 | "with the soft restoration phase. If the regular primal-dual step, " |
---|
221 | "damped onl to satisfty the fraction-to-the-boundary rule, is not " |
---|
222 | "decreasing the error by this factor, then the regular restoration " |
---|
223 | "phase is called. Choosing \"0\" here disables the soft " |
---|
224 | "restoration phase."); |
---|
225 | roptions->AddStringOption2( |
---|
226 | "start_with_resto", |
---|
227 | "Tells algorithm to switch to restoration phase in first iteration.", |
---|
228 | "no", |
---|
229 | "no", "don't force start in restoration phase", |
---|
230 | "yes", "force start in restoration phase", |
---|
231 | "Setting this option to yes forces the algorithm to switch to the " |
---|
232 | "restoration phase in the first iteration. If the initial point " |
---|
233 | "is feasible, the algorithm will abort with a failure."); |
---|
234 | roptions->AddLowerBoundedNumberOption( |
---|
235 | "tiny_step_tol", |
---|
236 | "Tolerance for detecting numerically insignificant steps.", |
---|
237 | 0.0, false, 10.0*std::numeric_limits<double>::epsilon(), |
---|
238 | "If the search direction in the primal variables (x and s) is, in " |
---|
239 | "relative terms for each component, less than this values, the " |
---|
240 | "algorithm accepts the full step without line search. The default " |
---|
241 | "value is 10 times machine precision."); |
---|
242 | roptions->AddLowerBoundedIntegerOption( |
---|
243 | "watchdog_shortened_iter_trigger", |
---|
244 | "Number of shortened iterations that trigger the watchdog.", |
---|
245 | 0, 10, |
---|
246 | "If the number of iterations in which the backtracking line search " |
---|
247 | "did not accept the first trial point exceedes this number, the " |
---|
248 | "watchdog procedure is activated. Choosing \"0\" here disables the " |
---|
249 | "watchdog procedure."); |
---|
250 | roptions->AddLowerBoundedIntegerOption( |
---|
251 | "watchdog_trial_iter_max", |
---|
252 | "Maximal number of watchdog iterations.", |
---|
253 | 1, 3, |
---|
254 | "Determines the number of trial iterations before the watchdog " |
---|
255 | "procedure is aborted and the algorithm returns to the stored point."); |
---|
256 | roptions->AddLowerBoundedNumberOption( |
---|
257 | "acceptable_tol", |
---|
258 | "Threshold for NLP error to consider iterate as acceptable.", |
---|
259 | 0.0, false, 1e-6, |
---|
260 | "Determines tolerance for which an iterate is considered as accepable " |
---|
261 | "solution. If the algorithm would trigger the restoration phase at " |
---|
262 | "such a point, it instead terminates, returning this acceptable " |
---|
263 | "point. Further, if the algorithm encounters \"acceptable_iter_max\" " |
---|
264 | "successive points satisfying this NLP tolerance, the algorithm " |
---|
265 | "terminates."); |
---|
266 | roptions->AddLowerBoundedIntegerOption( |
---|
267 | "acceptable_iter_max", |
---|
268 | "Number of acceptable iterates to trigger termination.", |
---|
269 | 0, 15, |
---|
270 | "if the algorithm encounters so many successive acceptable iterates " |
---|
271 | "(see \"acceptable_tol\"), it terminates, assuming that the problem " |
---|
272 | "has been solved to best possible accuracy given round-off."); |
---|
273 | } |
---|
274 | |
---|
275 | bool FilterLineSearch::InitializeImpl(const OptionsList& options, |
---|
276 | const std::string& prefix) |
---|
277 | { |
---|
278 | options.GetNumericValue("theta_max_fact", theta_max_fact_, prefix); |
---|
279 | options.GetNumericValue("theta_min_fact", theta_min_fact_, prefix); |
---|
280 | ASSERT_EXCEPTION(theta_min_fact_ < theta_max_fact_, OptionsList::OPTION_OUT_OF_RANGE, |
---|
281 | "Option \"theta_min_fact\": This value must be larger than 0 and less than theta_max_fact."); |
---|
282 | options.GetNumericValue("eta_phi", eta_phi_, prefix); |
---|
283 | options.GetNumericValue("delta", delta_, prefix); |
---|
284 | options.GetNumericValue("s_phi", s_phi_, prefix); |
---|
285 | options.GetNumericValue("s_theta", s_theta_, prefix); |
---|
286 | options.GetNumericValue("gamma_phi", gamma_phi_, prefix); |
---|
287 | options.GetNumericValue("gamma_theta", gamma_theta_, prefix); |
---|
288 | options.GetNumericValue("alpha_min_frac", alpha_min_frac_, prefix); |
---|
289 | options.GetNumericValue("alpha_red_factor", alpha_red_factor_, prefix); |
---|
290 | options.GetIntegerValue("max_soc", max_soc_, prefix); |
---|
291 | if (max_soc_>0) { |
---|
292 | ASSERT_EXCEPTION(IsValid(pd_solver_), OptionsList::OPTION_OUT_OF_RANGE, |
---|
293 | "Option \"max_soc\": This option is non-negative, but no linear solver for computing the SOC given to FilterLineSearch object."); |
---|
294 | } |
---|
295 | options.GetNumericValue("kappa_soc", kappa_soc_, prefix); |
---|
296 | options.GetNumericValue("obj_max_inc", obj_max_inc_, prefix); |
---|
297 | options.GetBoolValue("magic_steps", magic_steps_, prefix); |
---|
298 | Index enum_int; |
---|
299 | options.GetEnumValue("corrector_type", enum_int, prefix); |
---|
300 | corrector_type_ = CorrectorTypeEnum(enum_int); |
---|
301 | options.GetBoolValue("skip_corr_if_neg_curv", skip_corr_if_neg_curv_, prefix); |
---|
302 | options.GetBoolValue("skip_corr_in_monotone_mode", skip_corr_in_monotone_mode_, prefix); |
---|
303 | options.GetBoolValue("accept_every_trial_step", accept_every_trial_step_, prefix); |
---|
304 | options.GetEnumValue("alpha_for_y", enum_int, prefix); |
---|
305 | alpha_for_y_ = AlphaForYEnum(enum_int); |
---|
306 | options.GetNumericValue("corrector_compl_avrg_red_fact", corrector_compl_avrg_red_fact_, prefix); |
---|
307 | options.GetNumericValue("expect_infeasible_problem_ctol", expect_infeasible_problem_ctol_, prefix); |
---|
308 | options.GetBoolValue("expect_infeasible_problem", expect_infeasible_problem_, prefix); |
---|
309 | |
---|
310 | options.GetBoolValue("start_with_resto", start_with_resto_, prefix); |
---|
311 | |
---|
312 | bool retvalue = true; |
---|
313 | if (IsValid(resto_phase_)) { |
---|
314 | retvalue = resto_phase_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(), |
---|
315 | options, prefix); |
---|
316 | } |
---|
317 | |
---|
318 | options.GetNumericValue("soft_resto_pderror_reduction_factor", |
---|
319 | soft_resto_pderror_reduction_factor_, prefix); |
---|
320 | options.GetNumericValue("tiny_step_tol", tiny_step_tol_, prefix); |
---|
321 | options.GetIntegerValue("watchdog_trial_iter_max", watchdog_trial_iter_max_, prefix); |
---|
322 | options.GetIntegerValue("watchdog_shortened_iter_trigger", watchdog_shortened_iter_trigger_, prefix); |
---|
323 | options.GetNumericValue("acceptable_tol", acceptable_tol_, prefix); |
---|
324 | options.GetIntegerValue("acceptable_iter_max", acceptable_iter_max_, prefix); |
---|
325 | |
---|
326 | // ToDo decide if also the PDSystemSolver should be initialized here... |
---|
327 | |
---|
328 | rigorous_ = true; |
---|
329 | skipped_line_search_ = false; |
---|
330 | tiny_step_last_iteration_ = false; |
---|
331 | |
---|
332 | Reset(); |
---|
333 | |
---|
334 | count_successive_shortened_steps_ = 0; |
---|
335 | count_acceptable_iter_ = 0; |
---|
336 | |
---|
337 | acceptable_iterate_ = NULL; |
---|
338 | |
---|
339 | return retvalue; |
---|
340 | } |
---|
341 | |
---|
342 | void FilterLineSearch::FindAcceptableTrialPoint() |
---|
343 | { |
---|
344 | DBG_START_METH("FilterLineSearch::FindAcceptableTrialPoint", |
---|
345 | dbg_verbosity); |
---|
346 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
347 | "--> Starting filter line search in iteration %d <--\n", |
---|
348 | IpData().iter_count()); |
---|
349 | |
---|
350 | // If the problem is square, we want to enable the |
---|
351 | // expect_infeasible_problem option automatically so that the |
---|
352 | // restoration phase is entered soon |
---|
353 | if (IpCq().IsSquareProblem()) { |
---|
354 | expect_infeasible_problem_ = true; |
---|
355 | expect_infeasible_problem_ctol_ = 0.; |
---|
356 | } |
---|
357 | |
---|
358 | // Store current iterate if the optimality error is on acceptable |
---|
359 | // level to restored if things fail later |
---|
360 | if (IpCq().curr_nlp_error()<=acceptable_tol_) { |
---|
361 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
362 | "Storing current iterate as backup acceptable point.\n"); |
---|
363 | IpData().Append_info_string("A"); |
---|
364 | StoreAcceptablePoint(); |
---|
365 | } |
---|
366 | |
---|
367 | // First assume that line search will find an acceptable trial point |
---|
368 | skipped_line_search_ = false; |
---|
369 | |
---|
370 | // Set the values for the reference point |
---|
371 | if (!in_watchdog_) { |
---|
372 | reference_theta_ = IpCq().curr_constraint_violation(); |
---|
373 | reference_barr_ = IpCq().curr_barrier_obj(); |
---|
374 | reference_gradBarrTDelta_ = IpCq().curr_gradBarrTDelta(); |
---|
375 | } |
---|
376 | else { |
---|
377 | reference_theta_ = watchdog_theta_; |
---|
378 | reference_barr_ = watchdog_barr_; |
---|
379 | reference_gradBarrTDelta_ = watchdog_gradBarrTDelta_; |
---|
380 | } |
---|
381 | |
---|
382 | // Get the search directions (this will store the actual search |
---|
383 | // direction, possibly including higher order corrections) |
---|
384 | SmartPtr<IteratesVector> actual_delta = IpData().delta()->MakeNewContainer(); |
---|
385 | |
---|
386 | bool goto_resto = false; |
---|
387 | if (actual_delta->Asum() == 0. ) { |
---|
388 | // In this case, a search direction could not be computed, and |
---|
389 | // we should immediately go to the restoration phase. ToDo: Cue |
---|
390 | // off of a something else than the norm of the search direction |
---|
391 | goto_resto = true; |
---|
392 | } |
---|
393 | |
---|
394 | if (start_with_resto_) { |
---|
395 | // If the use requested to start with the restoration phase, |
---|
396 | // skip the lin e search and do exactly that. Reset the flag so |
---|
397 | // that this happens only once. |
---|
398 | goto_resto = true; |
---|
399 | start_with_resto_= false; |
---|
400 | } |
---|
401 | |
---|
402 | bool accept = false; |
---|
403 | bool corr_taken = false; |
---|
404 | bool soc_taken = false; |
---|
405 | Index n_steps = 0; |
---|
406 | Number alpha_primal = 0.; |
---|
407 | |
---|
408 | // Check if search direction becomes too small |
---|
409 | // ToDo: move this into place independent of this particular line search? |
---|
410 | bool tiny_step = (!goto_resto && DetectTinyStep()); |
---|
411 | |
---|
412 | if (in_watchdog_ && (goto_resto || tiny_step)) { |
---|
413 | // If the step could not be computed or is too small and the |
---|
414 | // watchdog is active, stop the watch dog and resume everything |
---|
415 | // from reference point |
---|
416 | StopWatchDog(actual_delta); |
---|
417 | goto_resto = false; |
---|
418 | tiny_step = false; |
---|
419 | } |
---|
420 | |
---|
421 | // Check if we want to wake up the watchdog |
---|
422 | if (watchdog_shortened_iter_trigger_ > 0 && |
---|
423 | !in_watchdog_ && !goto_resto && !tiny_step && |
---|
424 | !in_soft_resto_phase_ && !expect_infeasible_problem_ && |
---|
425 | watchdog_shortened_iter_ >= watchdog_shortened_iter_trigger_) { |
---|
426 | StartWatchDog(); |
---|
427 | } |
---|
428 | |
---|
429 | if (tiny_step) { |
---|
430 | alpha_primal = |
---|
431 | IpCq().curr_primal_frac_to_the_bound(IpData().curr_tau()); |
---|
432 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
433 | "Tiny step detected. Use step size alpha = %e unchecked\n", |
---|
434 | alpha_primal); |
---|
435 | IpData().SetTrialPrimalVariablesFromStep(alpha_primal, *IpData().delta()->x(), *IpData().delta()->s()); |
---|
436 | IpData().Set_info_ls_count(0); |
---|
437 | |
---|
438 | if (tiny_step_last_iteration_) { |
---|
439 | IpData().Set_info_alpha_primal_char('T'); |
---|
440 | IpData().Set_tiny_step_flag(true); |
---|
441 | } |
---|
442 | else { |
---|
443 | IpData().Set_info_alpha_primal_char('t'); |
---|
444 | } |
---|
445 | |
---|
446 | tiny_step_last_iteration_ = true; |
---|
447 | accept = true; |
---|
448 | } |
---|
449 | else { |
---|
450 | tiny_step_last_iteration_ = false; |
---|
451 | } |
---|
452 | |
---|
453 | if (!goto_resto && !tiny_step) { |
---|
454 | |
---|
455 | if (in_soft_resto_phase_) { |
---|
456 | bool satisfies_original_filter = false; |
---|
457 | // ToDo use tiny_step in TrySoftRestoStep? |
---|
458 | accept = TrySoftRestoStep(actual_delta, |
---|
459 | satisfies_original_filter); |
---|
460 | if (accept) { |
---|
461 | IpData().Set_info_alpha_primal_char('s'); |
---|
462 | if (satisfies_original_filter) { |
---|
463 | in_soft_resto_phase_ = false; |
---|
464 | IpData().Set_info_alpha_primal_char('S'); |
---|
465 | } |
---|
466 | } |
---|
467 | } |
---|
468 | else { |
---|
469 | bool done = false; |
---|
470 | bool skip_first_trial_point = false; |
---|
471 | bool evaluation_error; |
---|
472 | while (!done) { |
---|
473 | accept = DoBacktrackingLineSearch(skip_first_trial_point, |
---|
474 | alpha_primal, |
---|
475 | corr_taken, |
---|
476 | soc_taken, |
---|
477 | n_steps, |
---|
478 | evaluation_error, |
---|
479 | actual_delta); |
---|
480 | DBG_PRINT((1, "evaluation_error = %d\n", evaluation_error)); |
---|
481 | if (in_watchdog_) { |
---|
482 | if (accept) { |
---|
483 | in_watchdog_ = false; |
---|
484 | IpData().Append_info_string("W"); |
---|
485 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
486 | "Watch dog procedure successful!\n"); |
---|
487 | done = true; |
---|
488 | } |
---|
489 | else { |
---|
490 | watchdog_trial_iter_++; |
---|
491 | if (evaluation_error || |
---|
492 | watchdog_trial_iter_ > watchdog_trial_iter_max_) { |
---|
493 | StopWatchDog(actual_delta); |
---|
494 | skip_first_trial_point = true; |
---|
495 | } |
---|
496 | else { |
---|
497 | done = true; |
---|
498 | accept = true; |
---|
499 | } |
---|
500 | } |
---|
501 | } |
---|
502 | else { |
---|
503 | done = true; |
---|
504 | } |
---|
505 | } |
---|
506 | } /* else: if (in_soft_resto_phase_) { */ |
---|
507 | } /* if (!goto_resto && !tiny_step) { */ |
---|
508 | |
---|
509 | // If line search has been aborted because the step size becomes too small, |
---|
510 | // go to the restoration phase |
---|
511 | if (!accept) { |
---|
512 | // If we are not asked to do a rigorous line search, do no call |
---|
513 | // the restoration phase. |
---|
514 | if (!rigorous_) { |
---|
515 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Skipping call of restoration phase...\n"); |
---|
516 | skipped_line_search_=true; |
---|
517 | } |
---|
518 | else { |
---|
519 | // Check if we should start the soft restoration phase |
---|
520 | if (!in_soft_resto_phase_ && soft_resto_pderror_reduction_factor_>0. |
---|
521 | && !goto_resto && !expect_infeasible_problem_) { |
---|
522 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
523 | "--> Starting soft restoration phase <--\n"); |
---|
524 | // Augment the filter with the current point |
---|
525 | AugmentFilter(); |
---|
526 | |
---|
527 | // Try the current search direction for the soft restoration phase |
---|
528 | bool satisfies_original_filter; |
---|
529 | accept = TrySoftRestoStep(actual_delta, |
---|
530 | satisfies_original_filter); |
---|
531 | // If it has been accepted: If the original filter is also |
---|
532 | // satisfied, we can just take that step and continue with |
---|
533 | // the regular algorithm, otherwise we stay in the soft |
---|
534 | // restoration phase |
---|
535 | if (accept) { |
---|
536 | if (satisfies_original_filter) { |
---|
537 | IpData().Set_info_alpha_primal_char('S'); |
---|
538 | } |
---|
539 | else { |
---|
540 | in_soft_resto_phase_ = true; |
---|
541 | IpData().Set_info_alpha_primal_char('s'); |
---|
542 | } |
---|
543 | } |
---|
544 | } |
---|
545 | |
---|
546 | if (!accept) { |
---|
547 | if (!in_soft_resto_phase_) { |
---|
548 | // Augment the filter with the current point if we are |
---|
549 | // already in the soft restoration phase, this has been |
---|
550 | // done earlier |
---|
551 | AugmentFilter(); |
---|
552 | } |
---|
553 | if (IpCq().curr_nlp_error()<=acceptable_tol_) { |
---|
554 | THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED, |
---|
555 | "Restoration phase called at acceptable point."); |
---|
556 | } |
---|
557 | |
---|
558 | if (!IsValid(resto_phase_)) { |
---|
559 | THROW_EXCEPTION(IpoptException, "No Restoration Phase given to this Filter Line Search Object!"); |
---|
560 | } |
---|
561 | if (IpCq().curr_constraint_violation()<= |
---|
562 | 1e-2*Min(IpData().tol(),IpData().primal_inf_tol())) { |
---|
563 | bool found_acceptable = RestoreAcceptablePoint(); |
---|
564 | if (found_acceptable) { |
---|
565 | THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED, |
---|
566 | "Restoration phase called at almost feasible point, but acceptable point could be restore.\n"); |
---|
567 | } |
---|
568 | else { |
---|
569 | THROW_EXCEPTION(RESTORATION_FAILED, |
---|
570 | "Restoration phase called, but point is almost feasible."); |
---|
571 | } |
---|
572 | } |
---|
573 | |
---|
574 | // Set the info fields for the first output line in the |
---|
575 | // restoration phase which reflects why the restoration phase |
---|
576 | // was called |
---|
577 | IpData().Set_info_alpha_primal(alpha_primal); |
---|
578 | IpData().Set_info_alpha_dual(0.); |
---|
579 | IpData().Set_info_alpha_primal_char('R'); |
---|
580 | IpData().Set_info_ls_count(n_steps+1); |
---|
581 | |
---|
582 | accept = resto_phase_->PerformRestoration(); |
---|
583 | if (!accept) { |
---|
584 | bool found_acceptable = RestoreAcceptablePoint(); |
---|
585 | if (found_acceptable) { |
---|
586 | THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED, |
---|
587 | "Restoration phase failed, but acceptable point could be restore.\n"); |
---|
588 | } |
---|
589 | else { |
---|
590 | THROW_EXCEPTION(RESTORATION_FAILED, |
---|
591 | "Failed restoration phase!!!"); |
---|
592 | } |
---|
593 | } |
---|
594 | count_successive_shortened_steps_ = 0; |
---|
595 | count_acceptable_iter_ = 0; |
---|
596 | if (expect_infeasible_problem_) { |
---|
597 | expect_infeasible_problem_ = false; |
---|
598 | } |
---|
599 | in_soft_resto_phase_ = false; |
---|
600 | watchdog_shortened_iter_ = 0; |
---|
601 | } |
---|
602 | } |
---|
603 | } |
---|
604 | else if (!in_soft_resto_phase_ || tiny_step) { |
---|
605 | // we didn't do the restoration phase and are now updating the |
---|
606 | // dual variables of the trial point |
---|
607 | Number alpha_dual_max = |
---|
608 | IpCq().dual_frac_to_the_bound(IpData().curr_tau(), |
---|
609 | *actual_delta->z_L(), *actual_delta->z_U(), |
---|
610 | *actual_delta->v_L(), *actual_delta->v_U()); |
---|
611 | |
---|
612 | PerformDualStep(alpha_primal, alpha_dual_max, actual_delta); |
---|
613 | |
---|
614 | if (acceptable_iter_max_>0) { |
---|
615 | if (IpCq().curr_nlp_error()<=acceptable_tol_) { |
---|
616 | count_acceptable_iter_++; |
---|
617 | if (count_acceptable_iter_>=acceptable_iter_max_) { |
---|
618 | IpData().AcceptTrialPoint(); |
---|
619 | THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED, |
---|
620 | "Algorithm seems stuck at acceptable level."); |
---|
621 | } |
---|
622 | } |
---|
623 | else { |
---|
624 | count_acceptable_iter_=0; |
---|
625 | } |
---|
626 | } |
---|
627 | |
---|
628 | if (n_steps==0) { |
---|
629 | // accepted this if a full step was |
---|
630 | // taken |
---|
631 | count_successive_shortened_steps_ = 0; |
---|
632 | watchdog_shortened_iter_ = 0; |
---|
633 | } |
---|
634 | else { |
---|
635 | count_successive_shortened_steps_++; |
---|
636 | watchdog_shortened_iter_++; |
---|
637 | } |
---|
638 | |
---|
639 | if (expect_infeasible_problem_ && |
---|
640 | IpCq().curr_constraint_violation() <= expect_infeasible_problem_ctol_) { |
---|
641 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
642 | "Constraint violation is with %e less than expect_infeasible_problem_ctol.\nDisable expect_infeasible_problem_heuristic.\n", IpCq().curr_constraint_violation()); |
---|
643 | expect_infeasible_problem_ = false; |
---|
644 | } |
---|
645 | } |
---|
646 | } |
---|
647 | |
---|
648 | bool FilterLineSearch::DoBacktrackingLineSearch(bool skip_first_trial_point, |
---|
649 | Number& alpha_primal, |
---|
650 | bool& corr_taken, |
---|
651 | bool& soc_taken, |
---|
652 | Index& n_steps, |
---|
653 | bool& evaluation_error, |
---|
654 | SmartPtr<IteratesVector>& actual_delta) |
---|
655 | { |
---|
656 | evaluation_error = false; |
---|
657 | bool accept = false; |
---|
658 | |
---|
659 | DBG_START_METH("FilterLineSearch::DoBacktrackingLineSearch", |
---|
660 | dbg_verbosity); |
---|
661 | |
---|
662 | // Compute primal fraction-to-the-boundary value |
---|
663 | Number alpha_primal_max = |
---|
664 | IpCq().primal_frac_to_the_bound(IpData().curr_tau(), |
---|
665 | *actual_delta->x(), |
---|
666 | *actual_delta->s()); |
---|
667 | |
---|
668 | // Compute smallest step size allowed |
---|
669 | Number alpha_min = alpha_primal_max; |
---|
670 | if (!in_watchdog_) { |
---|
671 | alpha_min = CalculateAlphaMin(); |
---|
672 | } |
---|
673 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
674 | "minimal step size ALPHA_MIN = %E\n", alpha_min); |
---|
675 | |
---|
676 | // Start line search from maximal step size |
---|
677 | alpha_primal = alpha_primal_max; |
---|
678 | |
---|
679 | // Step size used in ftype and armijo tests |
---|
680 | Number alpha_primal_test = alpha_primal; |
---|
681 | if (in_watchdog_) { |
---|
682 | alpha_primal_test = watchdog_alpha_primal_test_; |
---|
683 | } |
---|
684 | |
---|
685 | if (skip_first_trial_point) { |
---|
686 | alpha_primal *= alpha_red_factor_; |
---|
687 | } |
---|
688 | |
---|
689 | filter_.Print(Jnlst()); |
---|
690 | |
---|
691 | if (corrector_type_!=NO_CORRECTOR && !skip_first_trial_point && |
---|
692 | (!skip_corr_if_neg_curv_ || IpData().info_regu_x()==0.) && |
---|
693 | (!skip_corr_in_monotone_mode_ || IpData().FreeMuMode()) ) { |
---|
694 | // Before we do the actual backtracking line search for the |
---|
695 | // regular primal-dual search direction, let's see if a step |
---|
696 | // including a higher-order correctior is already acceptable |
---|
697 | accept = TryCorrector(alpha_primal_test, |
---|
698 | alpha_primal, |
---|
699 | actual_delta); |
---|
700 | } |
---|
701 | if (accept) { |
---|
702 | corr_taken = true; |
---|
703 | } |
---|
704 | |
---|
705 | if (!accept) { |
---|
706 | // Loop over decreaseing step sizes until acceptable point is |
---|
707 | // found or until step size becomes too small |
---|
708 | |
---|
709 | while (alpha_primal>alpha_min || |
---|
710 | n_steps == 0) { // always allow the "full" step if it is |
---|
711 | // acceptable (even if alpha_primal<=alpha_min) |
---|
712 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
713 | "Starting checks for alpha (primal) = %8.2e\n", |
---|
714 | alpha_primal); |
---|
715 | |
---|
716 | try { |
---|
717 | // Compute the primal trial point |
---|
718 | IpData().SetTrialPrimalVariablesFromStep(alpha_primal, *actual_delta->x(), *actual_delta->s()); |
---|
719 | |
---|
720 | if (magic_steps_) { |
---|
721 | PerformMagicStep(); |
---|
722 | } |
---|
723 | |
---|
724 | // If it is acceptable, stop the search |
---|
725 | alpha_primal_test = alpha_primal; |
---|
726 | accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test); |
---|
727 | } |
---|
728 | catch(IpoptNLP::Eval_Error& e) { |
---|
729 | e.ReportException(Jnlst()); |
---|
730 | Jnlst().Printf(J_WARNING, J_LINE_SEARCH, |
---|
731 | "Warning: Cutting back alpha due to evaluation error\n"); |
---|
732 | accept = false; |
---|
733 | evaluation_error = true; |
---|
734 | } |
---|
735 | |
---|
736 | if (accept) { |
---|
737 | break; |
---|
738 | } |
---|
739 | |
---|
740 | if (in_watchdog_) { |
---|
741 | break; |
---|
742 | } |
---|
743 | |
---|
744 | // Decide if we want to go to the restoration phase in a |
---|
745 | // short cut to check if the problem is infeasible |
---|
746 | if (expect_infeasible_problem_) { |
---|
747 | if (count_successive_shortened_steps_>=5) { |
---|
748 | break; |
---|
749 | } |
---|
750 | } |
---|
751 | |
---|
752 | // try second order correction step if the function could |
---|
753 | // be evaluated |
---|
754 | // DoTo: check if we want to do SOC when watchdog is active |
---|
755 | if (!evaluation_error) { |
---|
756 | Number theta_curr = IpCq().curr_constraint_violation(); |
---|
757 | Number theta_trial = IpCq().trial_constraint_violation(); |
---|
758 | if (alpha_primal==alpha_primal_max && // i.e. first trial point |
---|
759 | theta_curr<=theta_trial && max_soc_>0) { |
---|
760 | // Try second order correction |
---|
761 | accept = TrySecondOrderCorrection(alpha_primal_test, |
---|
762 | alpha_primal, |
---|
763 | actual_delta); |
---|
764 | } |
---|
765 | if (accept) { |
---|
766 | soc_taken = true; |
---|
767 | break; |
---|
768 | } |
---|
769 | } |
---|
770 | |
---|
771 | // Point is not yet acceptable, try a shorter one |
---|
772 | alpha_primal *= alpha_red_factor_; |
---|
773 | n_steps++; |
---|
774 | } |
---|
775 | } /* if (!accept) */ |
---|
776 | |
---|
777 | char info_alpha_primal_char; |
---|
778 | if (!accept && in_watchdog_) { |
---|
779 | info_alpha_primal_char = 'w'; |
---|
780 | } |
---|
781 | else { |
---|
782 | // Augment the filter if required |
---|
783 | if (!IsFtype(alpha_primal_test) || |
---|
784 | !ArmijoHolds(alpha_primal_test)) { |
---|
785 | AugmentFilter(); |
---|
786 | info_alpha_primal_char = 'h'; |
---|
787 | } |
---|
788 | else { |
---|
789 | info_alpha_primal_char = 'f'; |
---|
790 | } |
---|
791 | } |
---|
792 | if (soc_taken) { |
---|
793 | info_alpha_primal_char = toupper(info_alpha_primal_char); |
---|
794 | } |
---|
795 | IpData().Set_info_alpha_primal_char(info_alpha_primal_char); |
---|
796 | IpData().Set_info_ls_count(n_steps+1); |
---|
797 | if (corr_taken) { |
---|
798 | IpData().Append_info_string("C"); |
---|
799 | } |
---|
800 | |
---|
801 | return accept; |
---|
802 | } |
---|
803 | |
---|
804 | bool FilterLineSearch::IsFtype(Number alpha_primal_test) |
---|
805 | { |
---|
806 | DBG_START_METH("FilterLineSearch::IsFtype", |
---|
807 | dbg_verbosity); |
---|
808 | DBG_ASSERT(reference_theta_>0. || reference_gradBarrTDelta_ < 0.0); |
---|
809 | return (reference_gradBarrTDelta_ < 0.0 && |
---|
810 | alpha_primal_test*pow(-reference_gradBarrTDelta_,s_phi_) > |
---|
811 | delta_*pow(reference_theta_,s_theta_)); |
---|
812 | } |
---|
813 | |
---|
814 | void FilterLineSearch::AugmentFilter() |
---|
815 | { |
---|
816 | DBG_START_METH("FilterLineSearch::AugmentFilter", |
---|
817 | dbg_verbosity); |
---|
818 | |
---|
819 | Number phi_add = reference_barr_ - gamma_phi_*reference_theta_; |
---|
820 | Number theta_add = (1.-gamma_theta_)*reference_theta_; |
---|
821 | |
---|
822 | filter_.AddEntry(phi_add, theta_add, IpData().iter_count()); |
---|
823 | } |
---|
824 | |
---|
825 | bool |
---|
826 | FilterLineSearch::CheckAcceptabilityOfTrialPoint(Number alpha_primal_test) |
---|
827 | { |
---|
828 | DBG_START_METH("FilterLineSearch::CheckAcceptabilityOfTrialPoint", |
---|
829 | dbg_verbosity); |
---|
830 | |
---|
831 | if (accept_every_trial_step_) { |
---|
832 | |
---|
833 | // We call the evaluation at the trial point here, so that an |
---|
834 | // exception will the thrown if there are problem during the |
---|
835 | // evaluation of the functions (in that case, we want to further |
---|
836 | // reduce the step size |
---|
837 | Number trial_barr = IpCq().trial_barrier_obj(); |
---|
838 | Number trial_theta = IpCq().trial_constraint_violation(); |
---|
839 | return true; |
---|
840 | } |
---|
841 | |
---|
842 | bool accept; |
---|
843 | |
---|
844 | // First compute the barrier function and constraint violation at the |
---|
845 | // current iterate and the trial point |
---|
846 | |
---|
847 | Number trial_theta = IpCq().trial_constraint_violation(); |
---|
848 | // Check if constraint violation is becoming too large |
---|
849 | if (theta_max_ < 0.0) { |
---|
850 | theta_max_ = theta_max_fact_*Max(1.0, reference_theta_); |
---|
851 | } |
---|
852 | if (theta_min_ < 0.0) { |
---|
853 | theta_min_ = theta_min_fact_*Max(1.0, reference_theta_); |
---|
854 | } |
---|
855 | |
---|
856 | if (theta_max_>0 && trial_theta>theta_max_) { |
---|
857 | return false; |
---|
858 | } |
---|
859 | |
---|
860 | Number trial_barr = IpCq().trial_barrier_obj(); |
---|
861 | DBG_ASSERT(FiniteNumber(trial_barr)); |
---|
862 | |
---|
863 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
864 | "Checking acceptability for trial step size alpha_primal_test=%13.6e:\n", alpha_primal_test); |
---|
865 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
866 | " New values of barrier function = %23.16e (reference %23.16e):\n", trial_barr, reference_barr_); |
---|
867 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
868 | " New values of constraint violation = %23.16e (reference %23.16e):\n", trial_theta, reference_theta_); |
---|
869 | |
---|
870 | // Check if point is acceptable w.r.t current iterate |
---|
871 | if (IsFtype(alpha_primal_test) && reference_theta_ <= theta_min_) { |
---|
872 | // Armijo condition for the barrier function has to be satisfied |
---|
873 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking Armijo Condition...\n"); |
---|
874 | accept = ArmijoHolds(alpha_primal_test); |
---|
875 | } |
---|
876 | else { |
---|
877 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking sufficient reduction...\n"); |
---|
878 | accept = IsAcceptableToCurrentIterate(trial_barr, trial_theta); |
---|
879 | } |
---|
880 | |
---|
881 | if (!accept) { |
---|
882 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Failed...\n"); |
---|
883 | return accept; |
---|
884 | } |
---|
885 | else { |
---|
886 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Succeeded...\n"); |
---|
887 | } |
---|
888 | |
---|
889 | // Now check if that pair is acceptable to the filter |
---|
890 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Checking filter acceptability...\n"); |
---|
891 | accept = IsAcceptableToCurrentFilter(trial_barr, trial_theta); |
---|
892 | if (!accept) { |
---|
893 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Failed...\n"); |
---|
894 | return accept; |
---|
895 | } |
---|
896 | else { |
---|
897 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Succeeded...\n"); |
---|
898 | } |
---|
899 | |
---|
900 | return accept; |
---|
901 | } |
---|
902 | |
---|
903 | bool FilterLineSearch::ArmijoHolds(Number alpha_primal_test) |
---|
904 | { |
---|
905 | /* |
---|
906 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
907 | "ArmijoHolds test with trial_barr = %25.16e reference_barr = %25.16e\n alpha_primal_test = %25.16e reference_gradBarrTDelta = %25.16e\n", IpCq().trial_barrier_obj(), reference_barr_,alpha_primal_test,reference_gradBarrTDelta_); |
---|
908 | */ |
---|
909 | return Compare_le(IpCq().trial_barrier_obj()-reference_barr_, |
---|
910 | eta_phi_*alpha_primal_test*reference_gradBarrTDelta_, |
---|
911 | reference_barr_); |
---|
912 | } |
---|
913 | |
---|
914 | Number FilterLineSearch::CalculateAlphaMin() |
---|
915 | { |
---|
916 | Number gBD = IpCq().curr_gradBarrTDelta(); |
---|
917 | Number curr_theta = IpCq().curr_constraint_violation(); |
---|
918 | Number alpha_min = gamma_theta_; |
---|
919 | |
---|
920 | if (gBD < 0) { |
---|
921 | alpha_min = Min( gamma_theta_, |
---|
922 | gamma_phi_*curr_theta/(-gBD)); |
---|
923 | if (curr_theta <= theta_min_) { |
---|
924 | alpha_min = Min( alpha_min, |
---|
925 | delta_*pow(curr_theta,s_theta_)/pow(-gBD,s_phi_) |
---|
926 | ); |
---|
927 | } |
---|
928 | } |
---|
929 | |
---|
930 | return alpha_min_frac_*alpha_min; |
---|
931 | } |
---|
932 | |
---|
933 | bool FilterLineSearch::IsAcceptableToCurrentIterate(Number trial_barr, |
---|
934 | Number trial_theta, |
---|
935 | bool called_from_restoration /*=false*/) const |
---|
936 | { |
---|
937 | DBG_START_METH("FilterLineSearch::IsAcceptableToCurrentIterate", |
---|
938 | dbg_verbosity); |
---|
939 | |
---|
940 | // Check if the barrier objective function is increasing to |
---|
941 | // rapidly (according to option obj_max_inc) |
---|
942 | if (!called_from_restoration && trial_barr > reference_barr_) { |
---|
943 | Number basval = 1.; |
---|
944 | if (fabs(reference_barr_)>10.) { |
---|
945 | basval = log10(fabs(reference_barr_)); |
---|
946 | } |
---|
947 | if (log10(trial_barr-reference_barr_)>obj_max_inc_+basval) { |
---|
948 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Rejecting trial point because barrier objective function increasing too rapidly (from %27.15e to %27.15e)\n",reference_barr_,trial_barr); |
---|
949 | return false; |
---|
950 | } |
---|
951 | } |
---|
952 | |
---|
953 | DBG_PRINT((1,"trial_barr = %e reference_barr = %e\n", trial_barr, reference_barr_)); |
---|
954 | DBG_PRINT((1,"trial_theta = %e reference_theta = %e\n", trial_theta, reference_theta_)); |
---|
955 | return (Compare_le(trial_theta, (1.-gamma_theta_)*reference_theta_, reference_theta_) |
---|
956 | || Compare_le(trial_barr-reference_barr_, -gamma_phi_*reference_theta_, reference_barr_)); |
---|
957 | } |
---|
958 | |
---|
959 | bool FilterLineSearch::IsAcceptableToCurrentFilter(Number trial_barr, Number trial_theta) const |
---|
960 | { |
---|
961 | return filter_.Acceptable(trial_barr, trial_theta); |
---|
962 | } |
---|
963 | |
---|
964 | bool FilterLineSearch::Compare_le(Number lhs, Number rhs, Number BasVal) |
---|
965 | { |
---|
966 | DBG_START_FUN("FilterLineSearch::Compare_le", |
---|
967 | dbg_verbosity); |
---|
968 | DBG_PRINT((1,"lhs = %27.16e rhs = %27.16e BasVal = %27.16e\n",lhs,rhs,BasVal)); |
---|
969 | // ToDo: Comparison based on machine precision |
---|
970 | return (lhs - rhs <= 1e-15*fabs(BasVal)); |
---|
971 | } |
---|
972 | |
---|
973 | bool FilterLineSearch::TrySoftRestoStep(SmartPtr<IteratesVector>& actual_delta, |
---|
974 | bool &satisfies_original_filter) |
---|
975 | { |
---|
976 | DBG_START_FUN("FilterLineSearch::TrySoftRestoStep", dbg_verbosity); |
---|
977 | |
---|
978 | satisfies_original_filter = false; |
---|
979 | |
---|
980 | // ToDo: Need to decide if we want to try a corrector step first |
---|
981 | |
---|
982 | // Compute the maximal step sizes (we use identical step sizes for |
---|
983 | // primal and dual variables |
---|
984 | Number alpha_primal_max = |
---|
985 | IpCq().primal_frac_to_the_bound(IpData().curr_tau(), |
---|
986 | *actual_delta->x(), |
---|
987 | *actual_delta->s()); |
---|
988 | Number alpha_dual_max = |
---|
989 | IpCq().dual_frac_to_the_bound(IpData().curr_tau(), |
---|
990 | *actual_delta->z_L(), |
---|
991 | *actual_delta->z_U(), |
---|
992 | *actual_delta->v_L(), |
---|
993 | *actual_delta->v_U()); |
---|
994 | Number alpha_max = Min(alpha_primal_max, alpha_dual_max); |
---|
995 | |
---|
996 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
997 | "Trying soft restoration phase step with step length %13.6e\n", |
---|
998 | alpha_max); |
---|
999 | |
---|
1000 | // Set the trial point |
---|
1001 | IpData().SetTrialPrimalVariablesFromStep(alpha_max, *actual_delta->x(), *actual_delta->s()); |
---|
1002 | PerformDualStep(alpha_max, alpha_max, |
---|
1003 | actual_delta); |
---|
1004 | |
---|
1005 | // Check if that point is acceptable with respect to the current |
---|
1006 | // original filter |
---|
1007 | |
---|
1008 | Number trial_barr; |
---|
1009 | Number trial_theta; |
---|
1010 | try { |
---|
1011 | trial_barr = IpCq().trial_barrier_obj(); |
---|
1012 | trial_theta = IpCq().trial_constraint_violation(); |
---|
1013 | } |
---|
1014 | catch(IpoptNLP::Eval_Error& e) { |
---|
1015 | e.ReportException(Jnlst()); |
---|
1016 | Jnlst().Printf(J_WARNING, J_LINE_SEARCH, |
---|
1017 | "Warning: Evaluation error during soft restoration phase step.\n"); |
---|
1018 | return false; |
---|
1019 | } |
---|
1020 | if (theta_max_<=0 || trial_theta<=theta_max_) { |
---|
1021 | if (IsAcceptableToCurrentIterate(trial_barr, trial_theta)) { |
---|
1022 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1023 | " Trial step acceptable with respect to original filter.\n"); |
---|
1024 | satisfies_original_filter = true; |
---|
1025 | return true; |
---|
1026 | } |
---|
1027 | } |
---|
1028 | |
---|
1029 | // Evaluate the optimality error at the new point |
---|
1030 | Number mu = .0; |
---|
1031 | if (!IpData().FreeMuMode()) { |
---|
1032 | mu = IpData().curr_mu(); |
---|
1033 | } |
---|
1034 | Number trial_pderror; |
---|
1035 | Number curr_pderror; |
---|
1036 | try { |
---|
1037 | trial_pderror = IpCq().trial_primal_dual_system_error(mu); |
---|
1038 | curr_pderror = IpCq().curr_primal_dual_system_error(mu); |
---|
1039 | } |
---|
1040 | catch(IpoptNLP::Eval_Error& e) { |
---|
1041 | e.ReportException(Jnlst()); |
---|
1042 | Jnlst().Printf(J_WARNING, J_LINE_SEARCH, |
---|
1043 | "Warning: Evaluation error during soft restoration phase step.\n"); |
---|
1044 | return false; |
---|
1045 | } |
---|
1046 | |
---|
1047 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1048 | " Primal-dual error at current point: %23.16e\n", curr_pderror); |
---|
1049 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1050 | " Primal-dual error at trial point : %23.16e\n", trial_pderror); |
---|
1051 | // Check if there is sufficient reduction in the optimality error |
---|
1052 | if (trial_pderror <= soft_resto_pderror_reduction_factor_*curr_pderror) { |
---|
1053 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1054 | " Trial step accepted.\n"); |
---|
1055 | return true; |
---|
1056 | } |
---|
1057 | |
---|
1058 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1059 | " Trial step rejected.\n"); |
---|
1060 | return false; |
---|
1061 | } |
---|
1062 | |
---|
1063 | void FilterLineSearch::StartWatchDog() |
---|
1064 | { |
---|
1065 | DBG_START_FUN("FilterLineSearch::StartWatchDog", dbg_verbosity); |
---|
1066 | |
---|
1067 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1068 | "Starting Watch Dog\n"); |
---|
1069 | |
---|
1070 | in_watchdog_ = true; |
---|
1071 | watchdog_iterate_ = IpData().curr(); |
---|
1072 | watchdog_delta_ = IpData().delta(); |
---|
1073 | watchdog_trial_iter_ = 0; |
---|
1074 | watchdog_alpha_primal_test_ = |
---|
1075 | IpCq().curr_primal_frac_to_the_bound(IpData().curr_tau()); |
---|
1076 | watchdog_theta_ = IpCq().curr_constraint_violation(); |
---|
1077 | watchdog_barr_ = IpCq().curr_barrier_obj(); |
---|
1078 | watchdog_gradBarrTDelta_ = IpCq().curr_gradBarrTDelta(); |
---|
1079 | } |
---|
1080 | |
---|
1081 | void FilterLineSearch::StopWatchDog(SmartPtr<IteratesVector>& actual_delta) |
---|
1082 | { |
---|
1083 | DBG_START_FUN("FilterLineSearch::StopWatchDog", dbg_verbosity); |
---|
1084 | |
---|
1085 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1086 | "Stopping Watch Dog\n"); |
---|
1087 | |
---|
1088 | IpData().Append_info_string("w"); |
---|
1089 | |
---|
1090 | in_watchdog_ = false; |
---|
1091 | |
---|
1092 | // Reset all fields in IpData to reference point |
---|
1093 | SmartPtr<IteratesVector> old_trial = watchdog_iterate_->MakeNewContainer(); |
---|
1094 | IpData().set_trial(old_trial); |
---|
1095 | IpData().AcceptTrialPoint(); |
---|
1096 | actual_delta = watchdog_delta_->MakeNewContainer(); |
---|
1097 | IpData().SetHaveAffineDeltas(false); |
---|
1098 | |
---|
1099 | // reset the stored watchdog iterates |
---|
1100 | watchdog_iterate_ = NULL; |
---|
1101 | watchdog_delta_ = NULL; |
---|
1102 | |
---|
1103 | watchdog_shortened_iter_ = 0; |
---|
1104 | |
---|
1105 | reference_theta_ = watchdog_theta_; |
---|
1106 | reference_barr_ = watchdog_barr_; |
---|
1107 | reference_gradBarrTDelta_ = watchdog_gradBarrTDelta_; |
---|
1108 | } |
---|
1109 | |
---|
1110 | void FilterLineSearch::Reset() |
---|
1111 | { |
---|
1112 | DBG_START_FUN("FilterLineSearch::Reset", dbg_verbosity); |
---|
1113 | in_soft_resto_phase_ = false; |
---|
1114 | |
---|
1115 | // Inactivate the watchdog and release all stored data |
---|
1116 | in_watchdog_ = false; |
---|
1117 | watchdog_iterate_ = NULL; |
---|
1118 | watchdog_delta_ = NULL; |
---|
1119 | watchdog_shortened_iter_ = 0; |
---|
1120 | |
---|
1121 | filter_.Clear(); |
---|
1122 | } |
---|
1123 | |
---|
1124 | void FilterLineSearch::PerformDualStep(Number alpha_primal, |
---|
1125 | Number alpha_dual, |
---|
1126 | SmartPtr<IteratesVector>& delta) |
---|
1127 | { |
---|
1128 | DBG_START_FUN("FilterLineSearch::PerformDualStep", dbg_verbosity); |
---|
1129 | |
---|
1130 | // set the bound multipliers from the step |
---|
1131 | IpData().SetTrialBoundMultipliersFromStep(alpha_dual, *delta->z_L(), *delta->z_U(), *delta->v_L(), *delta->v_U()); |
---|
1132 | |
---|
1133 | Number alpha_y; |
---|
1134 | switch (alpha_for_y_) { |
---|
1135 | case PRIMAL_ALPHA_FOR_Y: |
---|
1136 | alpha_y = alpha_primal; |
---|
1137 | break; |
---|
1138 | case DUAL_ALPHA_FOR_Y: |
---|
1139 | alpha_y = alpha_dual; |
---|
1140 | break; |
---|
1141 | case MIN_ALPHA_FOR_Y: |
---|
1142 | alpha_y = Min(alpha_dual, alpha_primal); |
---|
1143 | break; |
---|
1144 | case MAX_ALPHA_FOR_Y: |
---|
1145 | alpha_y = Max(alpha_dual, alpha_primal); |
---|
1146 | break; |
---|
1147 | case FULL_STEP_FOR_Y: |
---|
1148 | alpha_y = 1; |
---|
1149 | break; |
---|
1150 | case MIN_DUAL_INFEAS_ALPHA_FOR_Y: |
---|
1151 | case SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y: |
---|
1152 | // Here we compute the step size for y so that the dual |
---|
1153 | // infeasibility is minimized along delta_y |
---|
1154 | |
---|
1155 | // compute the dual infeasibility at new point with old y |
---|
1156 | SmartPtr<IteratesVector> temp_trial |
---|
1157 | = IpData().trial()->MakeNewContainer(); |
---|
1158 | temp_trial->Set_y_c(*IpData().curr()->y_c()); |
---|
1159 | temp_trial->Set_y_d(*IpData().curr()->y_d()); |
---|
1160 | IpData().set_trial(temp_trial); |
---|
1161 | SmartPtr<const Vector> dual_inf_x = IpCq().trial_grad_lag_x(); |
---|
1162 | SmartPtr<const Vector> dual_inf_s = IpCq().trial_grad_lag_s(); |
---|
1163 | |
---|
1164 | SmartPtr<Vector> new_jac_times_delta_y = |
---|
1165 | IpData().curr()->x()->MakeNew(); |
---|
1166 | new_jac_times_delta_y->AddTwoVectors(1., *IpCq().trial_jac_cT_times_vec(*delta->y_c()), |
---|
1167 | 1., *IpCq().trial_jac_dT_times_vec(*delta->y_d()), |
---|
1168 | 0.); |
---|
1169 | |
---|
1170 | Number a = pow(new_jac_times_delta_y->Nrm2(), 2.) + |
---|
1171 | pow(delta->y_d()->Nrm2(), 2.); |
---|
1172 | Number b = dual_inf_x->Dot(*new_jac_times_delta_y) - |
---|
1173 | dual_inf_s->Dot(*delta->y_d()); |
---|
1174 | |
---|
1175 | Number alpha = - b/a; |
---|
1176 | |
---|
1177 | if (alpha_for_y_==SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y) { |
---|
1178 | alpha_y = Min(Max(alpha_primal, alpha_dual), Max(alpha, Min(alpha_primal, alpha_dual))); |
---|
1179 | } |
---|
1180 | else { |
---|
1181 | alpha_y = Min(1., Max(0., alpha)); |
---|
1182 | } |
---|
1183 | break; |
---|
1184 | } |
---|
1185 | |
---|
1186 | // Set the eq multipliers from the step now that alpha_y |
---|
1187 | // has been calculated. |
---|
1188 | DBG_PRINT((1, "alpha_y = %e\n", alpha_y)); |
---|
1189 | DBG_PRINT_VECTOR(2, "delta_y_c", *delta->y_c()); |
---|
1190 | DBG_PRINT_VECTOR(2, "delta_y_d", *delta->y_d()); |
---|
1191 | IpData().SetTrialEqMultipliersFromStep(alpha_y, *delta->y_c(), *delta->y_d()); |
---|
1192 | |
---|
1193 | // Set some information for iteration summary output |
---|
1194 | IpData().Set_info_alpha_primal(alpha_primal); |
---|
1195 | IpData().Set_info_alpha_dual(alpha_dual); |
---|
1196 | } |
---|
1197 | |
---|
1198 | bool |
---|
1199 | FilterLineSearch::TrySecondOrderCorrection( |
---|
1200 | Number alpha_primal_test, |
---|
1201 | Number& alpha_primal, |
---|
1202 | SmartPtr<IteratesVector>& actual_delta) |
---|
1203 | { |
---|
1204 | DBG_START_METH("FilterLineSearch::TrySecondOrderCorrection", |
---|
1205 | dbg_verbosity); |
---|
1206 | |
---|
1207 | bool accept = false; |
---|
1208 | Index count_soc = 0; |
---|
1209 | |
---|
1210 | Number theta_soc_old = 0.; |
---|
1211 | Number theta_curr = IpCq().curr_constraint_violation(); |
---|
1212 | Number theta_trial = IpCq().trial_constraint_violation(); |
---|
1213 | Number alpha_primal_soc = alpha_primal; |
---|
1214 | |
---|
1215 | SmartPtr<Vector> c_soc = IpCq().curr_c()->MakeNew(); |
---|
1216 | SmartPtr<Vector> dms_soc = IpCq().curr_d_minus_s()->MakeNew(); |
---|
1217 | c_soc->Copy(*IpCq().curr_c()); |
---|
1218 | dms_soc->Copy(*IpCq().curr_d_minus_s()); |
---|
1219 | while (count_soc<max_soc_ && !accept && |
---|
1220 | (count_soc==0 || theta_trial<=kappa_soc_*theta_soc_old) ) { |
---|
1221 | theta_soc_old = theta_trial; |
---|
1222 | |
---|
1223 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1224 | "Trying second order correction number %d\n", |
---|
1225 | count_soc+1); |
---|
1226 | |
---|
1227 | // Compute SOC constraint violation |
---|
1228 | c_soc->AddOneVector(1.0, *IpCq().trial_c(), alpha_primal_soc); |
---|
1229 | dms_soc->AddOneVector(1.0, *IpCq().trial_d_minus_s(), alpha_primal_soc); |
---|
1230 | |
---|
1231 | // Compute the SOC search direction |
---|
1232 | SmartPtr<IteratesVector> delta_soc = actual_delta->MakeNewIteratesVector(true); |
---|
1233 | SmartPtr<IteratesVector> rhs = actual_delta->MakeNewContainer(); |
---|
1234 | rhs->Set_x(*IpCq().curr_grad_lag_with_damping_x()); |
---|
1235 | rhs->Set_s(*IpCq().curr_grad_lag_with_damping_s()); |
---|
1236 | rhs->Set_y_c(*c_soc); |
---|
1237 | rhs->Set_y_d(*dms_soc); |
---|
1238 | rhs->Set_z_L(*IpCq().curr_relaxed_compl_x_L()); |
---|
1239 | rhs->Set_z_U(*IpCq().curr_relaxed_compl_x_U()); |
---|
1240 | rhs->Set_v_L(*IpCq().curr_relaxed_compl_s_L()); |
---|
1241 | rhs->Set_v_U(*IpCq().curr_relaxed_compl_s_U()); |
---|
1242 | pd_solver_->Solve(-1.0, 0.0, *rhs, *delta_soc, true); |
---|
1243 | |
---|
1244 | // Compute step size |
---|
1245 | alpha_primal_soc = |
---|
1246 | IpCq().primal_frac_to_the_bound(IpData().curr_tau(), |
---|
1247 | *delta_soc->x(), |
---|
1248 | *delta_soc->s()); |
---|
1249 | |
---|
1250 | // Check if trial point is acceptable |
---|
1251 | try { |
---|
1252 | // Compute the primal trial point |
---|
1253 | IpData().SetTrialPrimalVariablesFromStep(alpha_primal_soc, *delta_soc->x(), *delta_soc->s()); |
---|
1254 | |
---|
1255 | // in acceptance tests, use original step size! |
---|
1256 | accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test); |
---|
1257 | } |
---|
1258 | catch(IpoptNLP::Eval_Error& e) { |
---|
1259 | e.ReportException(Jnlst()); |
---|
1260 | Jnlst().Printf(J_WARNING, J_MAIN, "Warning: SOC step rejected due to evaluation error\n"); |
---|
1261 | accept = false; |
---|
1262 | } |
---|
1263 | |
---|
1264 | if (accept) { |
---|
1265 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Second order correction step accepted with %d corrections.\n", count_soc+1); |
---|
1266 | // Accept all SOC quantities |
---|
1267 | alpha_primal = alpha_primal_soc; |
---|
1268 | actual_delta = delta_soc; |
---|
1269 | } |
---|
1270 | else { |
---|
1271 | count_soc++; |
---|
1272 | theta_trial = IpCq().trial_constraint_violation(); |
---|
1273 | } |
---|
1274 | } |
---|
1275 | return accept; |
---|
1276 | } |
---|
1277 | |
---|
1278 | bool |
---|
1279 | FilterLineSearch::TryCorrector( |
---|
1280 | Number alpha_primal_test, |
---|
1281 | Number& alpha_primal, |
---|
1282 | SmartPtr<IteratesVector>& actual_delta) |
---|
1283 | { |
---|
1284 | DBG_START_METH("FilterLineSearch::TryCorrector", |
---|
1285 | dbg_verbosity); |
---|
1286 | |
---|
1287 | Index n_bounds = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim() |
---|
1288 | + IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim(); |
---|
1289 | if (n_bounds==0) { |
---|
1290 | // Nothing to be done |
---|
1291 | return false; |
---|
1292 | } |
---|
1293 | |
---|
1294 | bool accept = false; |
---|
1295 | |
---|
1296 | // Compute the corrector step based on corrector_type parameter |
---|
1297 | // create a new iterates vector and allocate space for all the entries |
---|
1298 | SmartPtr<IteratesVector> delta_corr = actual_delta->MakeNewIteratesVector(true); |
---|
1299 | |
---|
1300 | switch (corrector_type_) { |
---|
1301 | case AFFINE_CORRECTOR : { |
---|
1302 | // 1: Standard MPC corrector |
---|
1303 | |
---|
1304 | if (!IpData().HaveAffineDeltas()) { |
---|
1305 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1306 | "Solving the Primal Dual System for the affine step\n"); |
---|
1307 | // First get the right hand side |
---|
1308 | SmartPtr<IteratesVector> rhs_aff = delta_corr->MakeNewContainer(); |
---|
1309 | |
---|
1310 | rhs_aff->Set_x(*IpCq().curr_grad_lag_x()); |
---|
1311 | rhs_aff->Set_s(*IpCq().curr_grad_lag_s()); |
---|
1312 | rhs_aff->Set_y_c(*IpCq().curr_c()); |
---|
1313 | rhs_aff->Set_y_d(*IpCq().curr_d_minus_s()); |
---|
1314 | rhs_aff->Set_z_L(*IpCq().curr_compl_x_L()); |
---|
1315 | rhs_aff->Set_z_U(*IpCq().curr_compl_x_U()); |
---|
1316 | rhs_aff->Set_v_L(*IpCq().curr_compl_s_L()); |
---|
1317 | rhs_aff->Set_v_U(*IpCq().curr_compl_s_U()); |
---|
1318 | |
---|
1319 | // create a new iterates vector (with allocated space) |
---|
1320 | // for the affine scaling step |
---|
1321 | SmartPtr<IteratesVector> step_aff = delta_corr->MakeNewIteratesVector(true); |
---|
1322 | |
---|
1323 | // Now solve the primal-dual system to get the step |
---|
1324 | pd_solver_->Solve(-1.0, 0.0, *rhs_aff, *step_aff, false); |
---|
1325 | |
---|
1326 | DBG_PRINT_VECTOR(2, "step_aff", *step_aff); |
---|
1327 | |
---|
1328 | IpData().set_delta_aff(step_aff); |
---|
1329 | IpData().SetHaveAffineDeltas(true); |
---|
1330 | } |
---|
1331 | |
---|
1332 | DBG_ASSERT(IpData().HaveAffineDeltas()); |
---|
1333 | |
---|
1334 | const SmartPtr<const IteratesVector> delta_aff = IpData().delta_aff(); |
---|
1335 | |
---|
1336 | delta_corr->Copy(*actual_delta); |
---|
1337 | |
---|
1338 | // create a rhs vector and allocate entries |
---|
1339 | SmartPtr<IteratesVector> rhs = actual_delta->MakeNewIteratesVector(true); |
---|
1340 | |
---|
1341 | rhs->x_NonConst()->Set(0.); |
---|
1342 | rhs->s_NonConst()->Set(0.); |
---|
1343 | rhs->y_c_NonConst()->Set(0.); |
---|
1344 | rhs->y_d_NonConst()->Set(0.); |
---|
1345 | IpNLP().Px_L()->TransMultVector(-1., *delta_aff->x(), 0., *rhs->z_L_NonConst()); |
---|
1346 | rhs->z_L_NonConst()->ElementWiseMultiply(*delta_aff->z_L()); |
---|
1347 | IpNLP().Px_U()->TransMultVector(1., *delta_aff->x(), 0., *rhs->z_U_NonConst()); |
---|
1348 | rhs->z_U_NonConst()->ElementWiseMultiply(*delta_aff->z_U()); |
---|
1349 | IpNLP().Pd_L()->TransMultVector(-1., *delta_aff->s(), 0., *rhs->v_L_NonConst()); |
---|
1350 | rhs->v_L_NonConst()->ElementWiseMultiply(*delta_aff->v_L()); |
---|
1351 | IpNLP().Pd_U()->TransMultVector(1., *delta_aff->s(), 0., *rhs->v_U_NonConst()); |
---|
1352 | rhs->v_U_NonConst()->ElementWiseMultiply(*delta_aff->v_U()); |
---|
1353 | |
---|
1354 | pd_solver_->Solve(1.0, 1.0, *rhs, *delta_corr, true); |
---|
1355 | |
---|
1356 | DBG_PRINT_VECTOR(2, "delta_corr", *delta_corr); |
---|
1357 | } |
---|
1358 | break; |
---|
1359 | case PRIMAL_DUAL_CORRECTOR : { |
---|
1360 | // 2: Second order correction for primal-dual step to |
---|
1361 | // primal-dual mu |
---|
1362 | |
---|
1363 | delta_corr->Copy(*actual_delta); |
---|
1364 | |
---|
1365 | // allocate space for the rhs |
---|
1366 | SmartPtr<IteratesVector> rhs = actual_delta->MakeNewIteratesVector(true); |
---|
1367 | |
---|
1368 | rhs->x_NonConst()->Set(0.); |
---|
1369 | rhs->s_NonConst()->Set(0.); |
---|
1370 | rhs->y_c_NonConst()->Set(0.); |
---|
1371 | rhs->y_d_NonConst()->Set(0.); |
---|
1372 | |
---|
1373 | Number mu = IpData().curr_mu(); |
---|
1374 | SmartPtr<Vector> tmp; |
---|
1375 | |
---|
1376 | rhs->z_L_NonConst()->Copy(*IpCq().curr_slack_x_L()); |
---|
1377 | IpNLP().Px_L()->TransMultVector(-1., *actual_delta->x(), |
---|
1378 | -1., *rhs->z_L_NonConst()); |
---|
1379 | tmp = actual_delta->z_L()->MakeNew(); |
---|
1380 | tmp->AddTwoVectors(1., *IpData().curr()->z_L(), 1., *actual_delta->z_L(), 0.); |
---|
1381 | rhs->z_L_NonConst()->ElementWiseMultiply(*tmp); |
---|
1382 | rhs->z_L_NonConst()->AddScalar(mu); |
---|
1383 | |
---|
1384 | rhs->z_U_NonConst()->Copy(*IpCq().curr_slack_x_U()); |
---|
1385 | IpNLP().Px_U()->TransMultVector(1., *actual_delta->x(), |
---|
1386 | -1., *rhs->z_U_NonConst()); |
---|
1387 | tmp = actual_delta->z_U()->MakeNew(); |
---|
1388 | tmp->AddTwoVectors(1., *IpData().curr()->z_U(), 1., *actual_delta->z_U(), 0.); |
---|
1389 | rhs->z_U_NonConst()->ElementWiseMultiply(*tmp); |
---|
1390 | rhs->z_U_NonConst()->AddScalar(mu); |
---|
1391 | |
---|
1392 | rhs->v_L_NonConst()->Copy(*IpCq().curr_slack_s_L()); |
---|
1393 | IpNLP().Pd_L()->TransMultVector(-1., *actual_delta->s(), |
---|
1394 | -1., *rhs->v_L_NonConst()); |
---|
1395 | tmp = actual_delta->v_L()->MakeNew(); |
---|
1396 | tmp->AddTwoVectors(1., *IpData().curr()->v_L(), 1., *actual_delta->v_L(), 0.); |
---|
1397 | rhs->v_L_NonConst()->ElementWiseMultiply(*tmp); |
---|
1398 | rhs->v_L_NonConst()->AddScalar(mu); |
---|
1399 | |
---|
1400 | rhs->v_U_NonConst()->Copy(*IpCq().curr_slack_s_U()); |
---|
1401 | IpNLP().Pd_U()->TransMultVector(1., *actual_delta->s(), |
---|
1402 | -1., *rhs->v_U_NonConst()); |
---|
1403 | tmp = actual_delta->v_U()->MakeNew(); |
---|
1404 | tmp->AddTwoVectors(1., *IpData().curr()->v_U(), 1., *actual_delta->v_U(), 0.); |
---|
1405 | rhs->v_U_NonConst()->ElementWiseMultiply(*tmp); |
---|
1406 | rhs->v_U_NonConst()->AddScalar(mu); |
---|
1407 | |
---|
1408 | DBG_PRINT_VECTOR(2, "rhs", *rhs); |
---|
1409 | |
---|
1410 | pd_solver_->Solve(1.0, 1.0, *rhs, *delta_corr, true); |
---|
1411 | |
---|
1412 | DBG_PRINT_VECTOR(2, "delta_corr", *delta_corr); |
---|
1413 | } |
---|
1414 | break; |
---|
1415 | default: |
---|
1416 | DBG_ASSERT("Unknown corrector_type value."); |
---|
1417 | } |
---|
1418 | |
---|
1419 | // Compute step size |
---|
1420 | Number alpha_primal_corr = |
---|
1421 | IpCq().primal_frac_to_the_bound(IpData().curr_tau(), |
---|
1422 | *delta_corr->x(), |
---|
1423 | *delta_corr->s()); |
---|
1424 | // Set the primal trial point |
---|
1425 | IpData().SetTrialPrimalVariablesFromStep(alpha_primal_corr, *delta_corr->x(), *delta_corr->s()); |
---|
1426 | |
---|
1427 | // Check if we want to not even try the filter criterion |
---|
1428 | Number alpha_dual_max = |
---|
1429 | IpCq().dual_frac_to_the_bound(IpData().curr_tau(), |
---|
1430 | *delta_corr->z_L(), *delta_corr->z_U(), |
---|
1431 | *delta_corr->v_L(), *delta_corr->v_U()); |
---|
1432 | |
---|
1433 | IpData().SetTrialBoundMultipliersFromStep(alpha_dual_max, *delta_corr->z_L(), *delta_corr->z_U(), *delta_corr->v_L(), *delta_corr->v_U()); |
---|
1434 | |
---|
1435 | Number trial_avrg_compl = IpCq().trial_avrg_compl(); |
---|
1436 | Number curr_avrg_compl = IpCq().curr_avrg_compl(); |
---|
1437 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1438 | "avrg_compl(curr) = %e, avrg_compl(trial) = %e\n", |
---|
1439 | curr_avrg_compl, trial_avrg_compl); |
---|
1440 | if (corrector_type_==AFFINE_CORRECTOR && |
---|
1441 | trial_avrg_compl>=corrector_compl_avrg_red_fact_*curr_avrg_compl) { |
---|
1442 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1443 | "Rejecting corrector step, because trial complementarity is too large.\n" ); |
---|
1444 | return false; |
---|
1445 | } |
---|
1446 | |
---|
1447 | // Check if trial point is acceptable |
---|
1448 | try { |
---|
1449 | // in acceptance tests, use original step size! |
---|
1450 | accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test); |
---|
1451 | } |
---|
1452 | catch(IpoptNLP::Eval_Error& e) { |
---|
1453 | e.ReportException(Jnlst()); |
---|
1454 | Jnlst().Printf(J_WARNING, J_MAIN, |
---|
1455 | "Warning: Corrector step rejected due to evaluation error\n"); |
---|
1456 | accept = false; |
---|
1457 | } |
---|
1458 | |
---|
1459 | if (accept) { |
---|
1460 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1461 | "Corrector step accepted with alpha_primal = %e\n", |
---|
1462 | alpha_primal_corr); |
---|
1463 | // Accept all SOC quantities |
---|
1464 | alpha_primal = alpha_primal_corr; |
---|
1465 | actual_delta = delta_corr; |
---|
1466 | |
---|
1467 | if (Jnlst().ProduceOutput(J_MOREVECTOR, J_MAIN)) { |
---|
1468 | Jnlst().Printf(J_MOREVECTOR, J_MAIN, |
---|
1469 | "*** Accepted corrector for Iteration: %d\n", |
---|
1470 | IpData().iter_count()); |
---|
1471 | Jnlst().PrintVector(J_MOREVECTOR, J_MAIN, |
---|
1472 | "delta_corr", *delta_corr); |
---|
1473 | } |
---|
1474 | } |
---|
1475 | |
---|
1476 | return accept; |
---|
1477 | } |
---|
1478 | |
---|
1479 | void |
---|
1480 | FilterLineSearch::PerformMagicStep() |
---|
1481 | { |
---|
1482 | DBG_START_METH("FilterLineSearch::PerformMagicStep", |
---|
1483 | dbg_verbosity); |
---|
1484 | |
---|
1485 | DBG_PRINT((1,"Incoming barr = %e and constrviol %e\n", |
---|
1486 | IpCq().trial_barrier_obj(), |
---|
1487 | IpCq().trial_constraint_violation())); |
---|
1488 | DBG_PRINT_VECTOR(2, "s in", *IpData().trial()->s()); |
---|
1489 | DBG_PRINT_VECTOR(2, "d minus s in", *IpCq().trial_d_minus_s()); |
---|
1490 | DBG_PRINT_VECTOR(2, "slack_s_L in", *IpCq().trial_slack_s_L()); |
---|
1491 | DBG_PRINT_VECTOR(2, "slack_s_U in", *IpCq().trial_slack_s_U()); |
---|
1492 | |
---|
1493 | SmartPtr<const Vector> d_L = IpNLP().d_L(); |
---|
1494 | SmartPtr<const Matrix> Pd_L = IpNLP().Pd_L(); |
---|
1495 | SmartPtr<Vector> delta_s_magic_L = d_L->MakeNew(); |
---|
1496 | delta_s_magic_L->Set(0.); |
---|
1497 | SmartPtr<Vector> tmp = d_L->MakeNew(); |
---|
1498 | Pd_L->TransMultVector(1., *IpCq().trial_d_minus_s(), 0., *tmp); |
---|
1499 | delta_s_magic_L->ElementWiseMax(*tmp); |
---|
1500 | |
---|
1501 | SmartPtr<const Vector> d_U = IpNLP().d_U(); |
---|
1502 | SmartPtr<const Matrix> Pd_U = IpNLP().Pd_U(); |
---|
1503 | SmartPtr<Vector> delta_s_magic_U = d_U->MakeNew(); |
---|
1504 | delta_s_magic_U->Set(0.); |
---|
1505 | tmp = d_U->MakeNew(); |
---|
1506 | Pd_U->TransMultVector(1., *IpCq().trial_d_minus_s(), 0., *tmp); |
---|
1507 | delta_s_magic_U->ElementWiseMin(*tmp); |
---|
1508 | |
---|
1509 | SmartPtr<Vector> delta_s_magic = IpData().trial()->s()->MakeNew(); |
---|
1510 | Pd_L->MultVector(1., *delta_s_magic_L, 0., *delta_s_magic); |
---|
1511 | Pd_U->MultVector(1., *delta_s_magic_U, 1., *delta_s_magic); |
---|
1512 | delta_s_magic_L = NULL; // free memory |
---|
1513 | delta_s_magic_U = NULL; // free memory |
---|
1514 | |
---|
1515 | // Now find those entries with both lower and upper bounds, there |
---|
1516 | // the step is too large |
---|
1517 | // ToDo this should only be done if there are inequality |
---|
1518 | // constraints with two bounds |
---|
1519 | // also this can be done in a smaller space (d_L or d_U whichever |
---|
1520 | // is smaller) |
---|
1521 | tmp = delta_s_magic->MakeNew(); |
---|
1522 | tmp->Copy(*IpData().trial()->s()); |
---|
1523 | Pd_L->MultVector(1., *d_L, -2., *tmp); |
---|
1524 | Pd_U->MultVector(1., *d_U, 1., *tmp); |
---|
1525 | SmartPtr<Vector> tmp2 = tmp->MakeNew(); |
---|
1526 | tmp2->Copy(*tmp); |
---|
1527 | tmp2->ElementWiseAbs(); |
---|
1528 | tmp->Axpy(-2., *delta_s_magic); |
---|
1529 | tmp->ElementWiseAbs(); |
---|
1530 | // now, tmp2 = |d_L + d_u - 2*s| and tmp = |d_L + d_u - 2*(s+Delta s)| |
---|
1531 | // we want to throw out those for which tmp2 > tmp |
---|
1532 | tmp->Axpy(-1., *tmp2); |
---|
1533 | tmp->ElementWiseSgn(); |
---|
1534 | tmp2->Set(0.); |
---|
1535 | tmp2->ElementWiseMax(*tmp); |
---|
1536 | tmp = d_L->MakeNew(); |
---|
1537 | Pd_L->TransMultVector(1., *tmp2, 0., *tmp); |
---|
1538 | Pd_L->MultVector(1., *tmp, 0., *tmp2); |
---|
1539 | tmp = d_U->MakeNew(); |
---|
1540 | Pd_U->TransMultVector(1., *tmp2, 0., *tmp); |
---|
1541 | Pd_U->MultVector(1., *tmp, 0., *tmp2); |
---|
1542 | DBG_PRINT_VECTOR(2, "tmp indicator", *tmp2) |
---|
1543 | // tmp2 now is one for those entries with both bounds, for which |
---|
1544 | // no step should be taken |
---|
1545 | |
---|
1546 | tmp = delta_s_magic->MakeNew(); |
---|
1547 | tmp->Copy(*delta_s_magic); |
---|
1548 | tmp->ElementWiseMultiply(*tmp2); |
---|
1549 | delta_s_magic->Axpy(-1., *tmp); |
---|
1550 | |
---|
1551 | Number delta_s_magic_max = delta_s_magic->Amax(); |
---|
1552 | Number mach_eps = std::numeric_limits<Number>::epsilon(); |
---|
1553 | if (delta_s_magic_max>0.) { |
---|
1554 | if (delta_s_magic_max > 10*mach_eps*IpData().trial()->s()->Amax()) { |
---|
1555 | IpData().Append_info_string("M"); |
---|
1556 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Magic step with max-norm %.6e taken.\n", delta_s_magic->Amax()); |
---|
1557 | Jnlst().PrintVector(J_MOREVECTOR, J_LINE_SEARCH, |
---|
1558 | "delta_s_magic", *delta_s_magic); |
---|
1559 | } |
---|
1560 | |
---|
1561 | // now finally compute the new overall slacks |
---|
1562 | delta_s_magic->Axpy(1., *IpData().trial()->s()); |
---|
1563 | SmartPtr<IteratesVector> trial = IpData().trial()->MakeNewContainer(); |
---|
1564 | trial->Set_s(*delta_s_magic); |
---|
1565 | IpData().set_trial(trial); |
---|
1566 | } |
---|
1567 | |
---|
1568 | DBG_PRINT((1,"Outgoing barr = %e and constrviol %e\n", IpCq().trial_barrier_obj(), IpCq().trial_constraint_violation())); |
---|
1569 | DBG_PRINT_VECTOR(2, "s out", *IpData().trial()->s()); |
---|
1570 | DBG_PRINT_VECTOR(2, "d minus s out", *IpCq().trial_d_minus_s()); |
---|
1571 | DBG_PRINT_VECTOR(2, "slack_s_L out", *IpCq().trial_slack_s_L()); |
---|
1572 | DBG_PRINT_VECTOR(2, "slack_s_U out", *IpCq().trial_slack_s_U()); |
---|
1573 | } |
---|
1574 | |
---|
1575 | bool |
---|
1576 | FilterLineSearch::DetectTinyStep() |
---|
1577 | { |
---|
1578 | DBG_START_METH("FilterLineSearch::DetectTinyStep", |
---|
1579 | dbg_verbosity); |
---|
1580 | |
---|
1581 | Number max_step_x; |
---|
1582 | Number max_step_s; |
---|
1583 | |
---|
1584 | if (tiny_step_tol_==0.) |
---|
1585 | return false; |
---|
1586 | |
---|
1587 | // ToDo try to find more efficient implementation |
---|
1588 | |
---|
1589 | SmartPtr<Vector> tmp = IpData().curr()->x()->MakeNew(); |
---|
1590 | tmp->Copy(*IpData().curr()->x()); |
---|
1591 | tmp->ElementWiseAbs(); |
---|
1592 | tmp->AddScalar(1.); |
---|
1593 | |
---|
1594 | SmartPtr<Vector> tmp2 = IpData().curr()->x()->MakeNew(); |
---|
1595 | tmp2->Copy(*IpData().delta()->x()); |
---|
1596 | tmp2->ElementWiseDivide(*tmp); |
---|
1597 | max_step_x = tmp2->Amax(); |
---|
1598 | Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH, |
---|
1599 | "Relative step size for delta_x = %e\n", |
---|
1600 | max_step_x); |
---|
1601 | if (max_step_x > tiny_step_tol_) |
---|
1602 | return false; |
---|
1603 | |
---|
1604 | tmp = IpData().curr()->s()->MakeNew(); |
---|
1605 | tmp->Copy(*IpData().curr()->s()); |
---|
1606 | tmp->ElementWiseAbs(); |
---|
1607 | tmp->AddScalar(1.); |
---|
1608 | |
---|
1609 | tmp2 = IpData().curr()->s()->MakeNew(); |
---|
1610 | tmp2->Copy(*IpData().delta()->s()); |
---|
1611 | tmp2->ElementWiseDivide(*tmp); |
---|
1612 | max_step_s = tmp2->Amax(); |
---|
1613 | Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH, |
---|
1614 | "Relative step size for delta_s = %e\n", |
---|
1615 | max_step_s); |
---|
1616 | if (max_step_s > tiny_step_tol_) |
---|
1617 | return false; |
---|
1618 | |
---|
1619 | Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, |
---|
1620 | "Tiny step of relative size %e detected.\n", |
---|
1621 | Max(max_step_x, max_step_s)); |
---|
1622 | |
---|
1623 | return true; |
---|
1624 | } |
---|
1625 | |
---|
1626 | void FilterLineSearch::StoreAcceptablePoint() |
---|
1627 | { |
---|
1628 | DBG_START_METH("FilterLineSearch::StoreAcceptablePoint", |
---|
1629 | dbg_verbosity); |
---|
1630 | |
---|
1631 | acceptable_iterate_ = IpData().curr(); |
---|
1632 | } |
---|
1633 | |
---|
1634 | bool FilterLineSearch::RestoreAcceptablePoint() |
---|
1635 | { |
---|
1636 | DBG_START_METH("FilterLineSearch::RestoreAcceptablePoint", |
---|
1637 | dbg_verbosity); |
---|
1638 | |
---|
1639 | if (!IsValid(acceptable_iterate_)) { |
---|
1640 | return false; |
---|
1641 | } |
---|
1642 | |
---|
1643 | SmartPtr<IteratesVector> prev_iterate = acceptable_iterate_->MakeNewContainer(); |
---|
1644 | IpData().set_trial(prev_iterate); |
---|
1645 | IpData().AcceptTrialPoint(); |
---|
1646 | |
---|
1647 | return true; |
---|
1648 | } |
---|
1649 | |
---|
1650 | |
---|
1651 | } // namespace Ipopt |
---|