source: branches/dev/Algorithm/IpRestoMinC_1Nrm.cpp @ 416

Last change on this file since 416 was 416, checked in by andreasw, 14 years ago
  • revised handling of "acceptable level of accuracy" (now in ConvCheck?)
  • fixed uncaught evaluation error exceptions
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.0 KB
Line 
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: IpRestoMinC_1Nrm.cpp 416 2005-07-29 19:11:41Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpRestoMinC_1Nrm.hpp"
10#include "IpCompoundVector.hpp"
11#include "IpRestoIpoptNLP.hpp"
12#include "IpDefaultIterateInitializer.hpp"
13
14namespace Ipopt
15{
16  DBG_SET_VERBOSITY(0);
17
18  DefineIpoptType(MinC_1NrmRestorationPhase);
19
20  MinC_1NrmRestorationPhase::MinC_1NrmRestorationPhase
21  (IpoptAlgorithm& resto_alg,
22   const SmartPtr<EqMultiplierCalculator>& eq_mult_calculator)
23      :
24      resto_alg_(&resto_alg),
25      eq_mult_calculator_(eq_mult_calculator),
26      resto_options_(NULL)
27  {
28    DBG_ASSERT(IsValid(resto_alg_));
29  }
30
31  MinC_1NrmRestorationPhase::~MinC_1NrmRestorationPhase()
32  {}
33
34  void MinC_1NrmRestorationPhase::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
35  {
36    roptions->AddLowerBoundedNumberOption(
37      "bound_mult_reset_threshold",
38      "Threshold for resetting bound multipliers after restoration phase.",
39      0.0, false,
40      1e3,
41      "After returning from the restoration phase, the bound multipliers are "
42      "updated with a Newton step for complementarity.  Here, the "
43      "change in the primal variables during the entire restoration "
44      "phase is taken to be the corresponding primal Newton step. "
45      "However, if after the update the largest bound multiplier "
46      "exceeds the threshold specified by this option, the multipliers "
47      "are all reset to 1.");
48    roptions->AddLowerBoundedNumberOption(
49      "constr_mult_reset_threshold",
50      "Threshold for resetting equality and inequality multipliers after restoration phase.",
51      0.0, false,
52      0e3,
53      "After returning from the restoration phase, the constraint multipliers "
54      "are recomputed by a least square estimate.  This option triggers when "
55      "those least-square esimates should be ignored.");
56  }
57
58  bool MinC_1NrmRestorationPhase::InitializeImpl(const OptionsList& options,
59      const std::string& prefix)
60  {
61    // keep a copy of these options to use when setting up the
62    // restoration phase
63    resto_options_ = new OptionsList(options);
64
65    options.GetNumericValue("constr_mult_reset_threshold",
66                            constr_mult_reset_threshold_,
67                            prefix);
68    options.GetNumericValue("bound_mult_reset_threshold",
69                            bound_mult_reset_threshold_,
70                            prefix);
71    options.GetBoolValue("expect_infeasible_problem",
72                         expect_infeasible_problem_,
73                         prefix);
74
75    // ToDo take care of this somewhere else?  avoid that the
76    // restoration phase is trigged by user option in first iteration
77    // of the restoration phase
78    resto_options_->SetValue("resto.start_with_resto", "no");
79
80    count_restorations_ = 0;
81
82    bool retvalue = true;
83    if (IsValid(eq_mult_calculator_)) {
84      retvalue = eq_mult_calculator_->Initialize(Jnlst(), IpNLP(), IpData(),
85                 IpCq(), options, prefix);
86    }
87    return retvalue;
88  }
89
90  bool
91  MinC_1NrmRestorationPhase::PerformRestoration()
92  {
93    DBG_START_METH("MinC_1NrmRestorationPhase::PerformRestoration",
94                   dbg_verbosity);
95
96    // Increase counter for restoration phase calls
97    count_restorations_++;
98    Jnlst().Printf(J_DETAILED, J_MAIN,
99                   "Starting Restoration Phase for the %d. time\n",
100                   count_restorations_);
101
102    DBG_ASSERT(IpCq().curr_constraint_violation()>0.);
103
104    // Create the restoration phase NLP etc objects
105    SmartPtr<IpoptData> resto_ip_data = new IpoptData();
106    SmartPtr<IpoptNLP> resto_ip_nlp =
107      new RestoIpoptNLP(IpNLP(), IpData(), IpCq());
108    SmartPtr<IpoptCalculatedQuantities> resto_ip_cq =
109      new IpoptCalculatedQuantities(resto_ip_nlp, resto_ip_data);
110
111    // Determine if this is a square problem
112    bool square_problem = IpCq().IsSquareProblem();
113
114    // Decide if we want to use the original option or want to make
115    // some changes
116    SmartPtr<OptionsList> actual_resto_options = resto_options_;
117    if (expect_infeasible_problem_ && count_restorations_==1) {
118      actual_resto_options = new OptionsList(*resto_options_);
119      // Ask for significant reduction of infeasibility, in the hope
120      // that we do not return from the restoration phase is the
121      // problem is infeasible
122      actual_resto_options->SetNumericValue("resto.required_infeasibility_reduction", 1e-3);
123    }
124    else if(square_problem) {
125      actual_resto_options = new OptionsList(*resto_options_);
126      // If this is a square problem, the want the restoration phase
127      // never to be left until the problem is converged
128      actual_resto_options->SetNumericValue("resto.required_infeasibility_reduction", 0.);
129    }
130
131    // Initialize the restoration phase algorithm
132    resto_alg_->Initialize(Jnlst(), *resto_ip_nlp, *resto_ip_data,
133                           *resto_ip_cq, *actual_resto_options, "resto.");
134
135    // Set iteration counter and info field for the restoration phase
136    resto_ip_data->Set_iter_count(IpData().iter_count()+1);
137    resto_ip_data->Set_info_regu_x(IpData().info_regu_x());
138    resto_ip_data->Set_info_alpha_primal(IpData().info_alpha_primal());
139    resto_ip_data->Set_info_alpha_primal_char(IpData().info_alpha_primal_char());
140    resto_ip_data->Set_info_alpha_dual(IpData().info_alpha_dual());
141    resto_ip_data->Set_info_ls_count(IpData().info_ls_count());
142
143    // Call the optimization algorithm to solve the restoration phase
144    // problem
145    SolverReturn resto_status   = resto_alg_->Optimize();
146
147    int retval=-1;
148
149    if (resto_status == SUCCESS) {
150      if (Jnlst().ProduceOutput(J_DETAILED, J_LINE_SEARCH)) {
151        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
152                       "\nRESTORATION PHASE RESULTS\n");
153        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
154                       "\n\nOptimal solution found! \n");
155        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
156                       "Optimal Objective Value = %.16E\n", resto_ip_cq->curr_f());
157        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
158                       "Number of Iterations = %d\n", resto_ip_data->iter_count());
159      }
160      if (Jnlst().ProduceOutput(J_VECTOR, J_LINE_SEARCH)) {
161        Jnlst().PrintVector(J_VECTOR, J_LINE_SEARCH,
162                            "x", *resto_ip_data->curr()->x());
163        Jnlst().PrintVector(J_VECTOR, J_LINE_SEARCH,
164                            "y_c", *resto_ip_data->curr()->y_c());
165        Jnlst().PrintVector(J_VECTOR, J_LINE_SEARCH,
166                            "y_d", *resto_ip_data->curr()->y_d());
167        Jnlst().PrintVector(J_VECTOR, J_LINE_SEARCH,
168                            "z_L", *resto_ip_data->curr()->z_L());
169        Jnlst().PrintVector(J_VECTOR, J_LINE_SEARCH,
170                            "z_U", *resto_ip_data->curr()->z_U());
171        Jnlst().PrintVector(J_VECTOR, J_LINE_SEARCH,
172                            "v_L", *resto_ip_data->curr()->v_L());
173        Jnlst().PrintVector(J_VECTOR, J_LINE_SEARCH,
174                            "v_U", *resto_ip_data->curr()->v_U());
175      }
176
177      retval = 0;
178    }
179    else if (resto_status == STOP_AT_TINY_STEP ||
180             resto_status == STOP_AT_ACCEPTABLE_POINT) {
181      Number orig_primal_inf =
182        IpCq().curr_primal_infeasibility(NORM_MAX);
183      // ToDo make the factor in following line an option
184      if (orig_primal_inf <= 1e2*IpData().tol()) {
185        THROW_EXCEPTION(RESTORATION_CONVERGED_TO_FEASIBLE_POINT,
186                        "Restoration phase converged to a point with small primal infeasibility");
187      }
188      else {
189        THROW_EXCEPTION(LOCALLY_INFEASIBLE,
190                        "Restoration phase converged to a point of local infeasibility");
191      }
192    }
193    else if (resto_status == MAXITER_EXCEEDED) {
194      THROW_EXCEPTION(IpoptException, "Maximal number of iterations exceeded in restoration phase.");
195      retval = 1;
196    }
197    else if (resto_status == LOCAL_INFEASIBILITY) {
198      // converged to locally infeasible point - pass this on to the outer algorithm...
199      THROW_EXCEPTION(LOCALLY_INFEASIBLE, "Restoration phase converged to a point of local infeasibility");
200    }
201    else if (resto_status == RESTORATION_FAILURE) {
202      THROW_EXCEPTION(RESTORATION_FAILED, "Restoration phase in the restoration phase failed.");
203    }
204    else {
205      Jnlst().Printf(J_ERROR, J_MAIN, "Sorry, things failed ?!?!\n");
206      retval = 1;
207    }
208
209    if (retval == 0) {
210      // Copy the results into the trial fields;. They will be
211      // accepted later in the full algorithm
212      SmartPtr<const CompoundVector> cx =
213        dynamic_cast<const CompoundVector*>(GetRawPtr(resto_ip_data->curr()->x()));
214      DBG_ASSERT(IsValid(cx));
215      SmartPtr<IteratesVector> trial = IpData().trial()->MakeNewContainer();
216      trial->Set_primal(*cx->GetComp(0), *resto_ip_data->curr()->s());
217      IpData().set_trial(trial);
218
219      // If this is a square problem, we are done because a
220      // sufficiently feasible point has been found
221      if (square_problem) {
222        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
223                       "Recursive restoration phase algorithm termined successfully for square problem.\n");
224        IpData().AcceptTrialPoint();
225        THROW_EXCEPTION(FEASIBILITY_PROBLEM_SOLVED, "Restoration phase converged to sufficiently feasible point of original square problem.");
226      }
227
228      // Update the bound multiplers, pretending that the entire
229      // progress in x and s in the restoration phase has been one
230      // [rimal-dual Newton step (and therefore the result of solving
231      // an augmented system)
232      SmartPtr<IteratesVector> delta = IpData().curr()->MakeNewIteratesVector(true);
233      delta->Set(0.0);
234      ComputeBoundMultiplierStep(*delta->z_L_NonConst(), *IpData().curr()->z_L(),
235                                 *IpCq().curr_slack_x_L(),
236                                 *IpCq().trial_slack_x_L());
237      ComputeBoundMultiplierStep(*delta->z_U_NonConst(), *IpData().curr()->z_U(),
238                                 *IpCq().curr_slack_x_U(),
239                                 *IpCq().trial_slack_x_U());
240      ComputeBoundMultiplierStep(*delta->v_L_NonConst(), *IpData().curr()->v_L(),
241                                 *IpCq().curr_slack_s_L(),
242                                 *IpCq().trial_slack_s_L());
243      ComputeBoundMultiplierStep(*delta->v_U_NonConst(), *IpData().curr()->v_U(),
244                                 *IpCq().curr_slack_s_U(),
245                                 *IpCq().trial_slack_s_U());
246
247      DBG_PRINT_VECTOR(1, "delta_z_L", *delta->z_L());
248      DBG_PRINT_VECTOR(1, "delta_z_U", *delta->z_U());
249      DBG_PRINT_VECTOR(1, "delta_v_L", *delta->v_L());
250      DBG_PRINT_VECTOR(1, "delta_v_U", *delta->v_U());
251
252      Number alpha_dual = IpCq().dual_frac_to_the_bound(IpData().curr_tau(),
253                          *delta->z_L_NonConst(),
254                          *delta->z_U_NonConst(),
255                          *delta->v_L_NonConst(),
256                          *delta->v_U_NonConst());
257      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
258                     "Step size for bound multipliers: %8.2e\n", alpha_dual);
259
260      IpData().SetTrialBoundMultipliersFromStep(alpha_dual, *delta->z_L(), *delta->z_U(), *delta->v_L(), *delta->v_U() );
261
262      // ToDo: Check what to do here:
263      Number bound_mult_max = Max(IpData().trial()->z_L()->Amax(),
264                                  IpData().trial()->z_U()->Amax(),
265                                  IpData().trial()->v_L()->Amax(),
266                                  IpData().trial()->v_U()->Amax());
267      if (bound_mult_max > bound_mult_reset_threshold_) {
268        trial = IpData().trial()->MakeNewContainer();
269        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
270                       "Bound multipliers after restoration phase too large (max=%8.2e). Set all to 1.\n",
271                       bound_mult_max);
272        trial->create_new_z_L();
273        trial->create_new_z_U();
274        trial->create_new_v_L();
275        trial->create_new_v_U();
276        trial->z_L_NonConst()->Set(1.0);
277        trial->z_U_NonConst()->Set(1.0);
278        trial->v_L_NonConst()->Set(1.0);
279        trial->v_U_NonConst()->Set(1.0);
280        IpData().set_trial(trial);
281
282      }
283
284      DefaultIterateInitializer::least_square_mults(
285        Jnlst(), IpNLP(), IpData(), IpCq(),
286        eq_mult_calculator_, constr_mult_reset_threshold_);
287
288      DBG_PRINT_VECTOR(2, "y_c", *IpData().curr()->y_c());
289      DBG_PRINT_VECTOR(2, "y_d", *IpData().curr()->y_d());
290
291      IpData().Set_iter_count(resto_ip_data->iter_count()-1);
292      // Skip the next line, because it would just replicate the first
293      // on during the restoration phase.
294      IpData().Set_info_skip_output(true);
295    }
296
297    return (retval == 0);
298  }
299
300  void MinC_1NrmRestorationPhase::ComputeBoundMultiplierStep(Vector& delta_z,
301      const Vector& curr_z,
302      const Vector& curr_slack,
303      const Vector& trial_slack)
304  {
305    Number mu = IpData().curr_mu();
306
307    delta_z.Copy(curr_slack);
308    delta_z.Axpy(-1., trial_slack);
309    delta_z.ElementWiseMultiply(curr_z);
310    delta_z.AddScalar(mu);
311    delta_z.ElementWiseDivide(curr_slack);
312    delta_z.Axpy(-1., curr_z);
313  }
314
315} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.