MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_TangentialStepWithInequStd_Step.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <math.h>
43 
44 #include <ostream>
45 #include <sstream>
46 
47 #include "MoochoPack_TangentialStepWithInequStd_Step.hpp"
48 #include "MoochoPack_moocho_algo_conversion.hpp"
49 #include "MoochoPack_Exceptions.hpp"
50 #include "IterationPack_print_algorithm_step.hpp"
51 #include "ConstrainedOptPack_MatrixIdentConcat.hpp"
52 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
53 #include "AbstractLinAlgPack_VectorOut.hpp"
54 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
55 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
56 #include "Teuchos_dyn_cast.hpp"
57 
58 namespace {
59 template< class T >
60 inline
61 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
62 template< class T >
63 inline
64 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
65 } // end namespace
66 
67 namespace MoochoPack {
68 
69 TangentialStepWithInequStd_Step::TangentialStepWithInequStd_Step(
70  const qp_solver_ptr_t &qp_solver
71  ,const qp_tester_ptr_t &qp_tester
72  ,value_type warm_start_frac
73  ,EQPTesting qp_testing
74  ,bool primal_feasible_point_error
75  ,bool dual_feasible_point_error
76  )
77  :qp_solver_(qp_solver)
78  ,qp_tester_(qp_tester)
79  ,warm_start_frac_(warm_start_frac)
80  ,qp_testing_(qp_testing)
81  ,primal_feasible_point_error_(primal_feasible_point_error)
82  ,dual_feasible_point_error_(dual_feasible_point_error)
83  ,dl_iq_(dl_name)
84  ,du_iq_(du_name)
85 {}
86 
88  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
89  ,poss_type assoc_step_poss
90  )
91 {
92 
93  using Teuchos::RCP;
94  using Teuchos::dyn_cast;
95  using ::fabs;
96  using LinAlgOpPack::Vt_S;
97  using LinAlgOpPack::V_VpV;
98  using LinAlgOpPack::V_VmV;
100  using LinAlgOpPack::Vp_V;
101  using LinAlgOpPack::V_StV;
102  using LinAlgOpPack::V_MtV;
103 // using ConstrainedOptPack::min_abs;
105  typedef VectorMutable::vec_mut_ptr_t vec_mut_ptr_t;
106 
107  NLPAlgo &algo = rsqp_algo(_algo);
108  NLPAlgoState &s = algo.rsqp_state();
109  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
110  EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
111  std::ostream &out = algo.track().journal_out();
112  //const bool check_results = algo.algo_cntr().check_results();
113 
114  // print step header.
115  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
116  using IterationPack::print_algorithm_step;
117  print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
118  }
119 
120  // problem dimensions
121  const size_type
122  //n = algo.nlp().n(),
123  m = algo.nlp().m(),
124  r = s.equ_decomp().size();
125 
126  // Get the iteration quantity container objects
127  IterQuantityAccess<value_type>
128  &alpha_iq = s.alpha(),
129  &zeta_iq = s.zeta(),
130  &eta_iq = s.eta();
131  IterQuantityAccess<VectorMutable>
132  &dl_iq = dl_iq_(s),
133  &du_iq = du_iq_(s),
134  &nu_iq = s.nu(),
135  *c_iq = m > 0 ? &s.c() : NULL,
136  *lambda_iq = m > 0 ? &s.lambda() : NULL,
137  &rGf_iq = s.rGf(),
138  &w_iq = s.w(),
139  &qp_grad_iq = s.qp_grad(),
140  &py_iq = s.py(),
141  &pz_iq = s.pz(),
142  &Ypy_iq = s.Ypy(),
143  &Zpz_iq = s.Zpz();
144  IterQuantityAccess<MatrixOp>
145  &Z_iq = s.Z(),
146  //*Uz_iq = (m > r) ? &s.Uz() : NULL,
147  *Uy_iq = (m > r) ? &s.Uy() : NULL;
148  IterQuantityAccess<MatrixSymOp>
149  &rHL_iq = s.rHL();
150  IterQuantityAccess<ActSetStats>
151  &act_set_stats_iq = act_set_stats_(s);
152 
153  // Accessed/modified/updated (just some)
154  VectorMutable *Ypy_k = (m ? &Ypy_iq.get_k(0) : NULL);
155  const MatrixOp &Z_k = Z_iq.get_k(0);
156  VectorMutable &pz_k = pz_iq.set_k(0);
157  VectorMutable &Zpz_k = Zpz_iq.set_k(0);
158 
159  // Comupte qp_grad which is an approximation to rGf + Z'*HL*Y*py
160 
161  // qp_grad = rGf
162  VectorMutable
163  &qp_grad_k = ( qp_grad_iq.set_k(0) = rGf_iq.get_k(0) );
164 
165  // qp_grad += zeta * w
166  if( w_iq.updated_k(0) ) {
167  if(zeta_iq.updated_k(0))
168  Vp_StV( &qp_grad_k, zeta_iq.get_k(0), w_iq.get_k(0) );
169  else
170  Vp_V( &qp_grad_k, w_iq.get_k(0) );
171  }
172 
173  //
174  // Set the bounds for:
175  //
176  // dl <= Z*pz + Y*py <= du -> dl - Ypy <= Z*pz <= du - Ypz
177 
178  vec_mut_ptr_t
179  bl = s.space_x().create_member(),
180  bu = s.space_x().create_member();
181 
182  if(m) {
183  // bl = dl_k - Ypy_k
184  V_VmV( bl.get(), dl_iq.get_k(0), *Ypy_k );
185  // bu = du_k - Ypy_k
186  V_VmV( bu.get(), du_iq.get_k(0), *Ypy_k );
187  }
188  else {
189  *bl = dl_iq.get_k(0);
190  *bu = du_iq.get_k(0);
191  }
192 
193  // Print out the QP bounds for the constraints
194  if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
195  out << "\nqp_grad_k = \n" << qp_grad_k;
196  }
197  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
198  out << "\nbl = \n" << *bl;
199  out << "\nbu = \n" << *bu;
200  }
201 
202  //
203  // Determine if we should perform a warm start or not.
204  //
205  bool do_warm_start = false;
206  if( act_set_stats_iq.updated_k(-1) ) {
207  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
208  out << "\nDetermining if the QP should use a warm start ...\n";
209  }
210  // We need to see if we should preform a warm start for the next iteration
211  ActSetStats &stats = act_set_stats_iq.get_k(-1);
212  const size_type
213  num_active = stats.num_active(),
214  num_adds = stats.num_adds(),
215  num_drops = stats.num_drops();
216  const value_type
217  frac_same
218  = ( num_adds == ActSetStats::NOT_KNOWN || num_active == 0
219  ? 0.0
220  : my_max(((double)(num_active)-num_adds-num_drops) / num_active, 0.0 ) );
221  do_warm_start = ( num_active > 0 && frac_same >= warm_start_frac() );
222  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
223  out << "\nnum_active = " << num_active;
224  if( num_active ) {
225  out << "\nmax(num_active-num_adds-num_drops,0)/(num_active) = "
226  << "max("<<num_active<<"-"<<num_adds<<"-"<<num_drops<<",0)/("<<num_active<<") = "
227  << frac_same;
228  if( do_warm_start )
229  out << " >= ";
230  else
231  out << " < ";
232  out << "warm_start_frac = " << warm_start_frac();
233  }
234  if( do_warm_start )
235  out << "\nUse a warm start this time!\n";
236  else
237  out << "\nDon't use a warm start this time!\n";
238  }
239  }
240 
241  // Use active set from last iteration as an estimate for current active set
242  // if we are to use a warm start.
243  //
244  // ToDo: If the selection of dependent and independent variables changes
245  // then you will have to adjust this or not perform a warm start at all!
246  if( do_warm_start ) {
247  nu_iq.set_k(0,-1);
248  }
249  else {
250  nu_iq.set_k(0) = 0.0; // No guess of the active set
251  }
252  VectorMutable &nu_k = nu_iq.get_k(0);
253 
254  //
255  // Setup the reduced QP subproblem
256  //
257  // The call to the QP is setup for the more flexible call to the QPSolverRelaxed
258  // interface to deal with the three independent variabilities: using simple
259  // bounds for pz or not, general inequalities included or not, and extra equality
260  // constraints included or not.
261  // If this method of calling the QP solver were not used then 4 separate
262  // calls to solve_qp(...) would have to be included to handle the four possible
263  // QP formulations.
264  //
265 
266  // The numeric arguments for the QP solver (in the nomenclatrue of QPSolverRelaxed)
267 
268  const value_type qp_bnd_inf = NLP::infinite_bound();
269 
270  const Vector &qp_g = qp_grad_k;
271  const MatrixSymOp &qp_G = rHL_iq.get_k(0);
272  const value_type qp_etaL = 0.0;
273  vec_mut_ptr_t qp_dL = Teuchos::null;
274  vec_mut_ptr_t qp_dU = Teuchos::null;
276  qp_E = Teuchos::null;
278  vec_mut_ptr_t qp_b = Teuchos::null;
279  vec_mut_ptr_t qp_eL = Teuchos::null;
280  vec_mut_ptr_t qp_eU = Teuchos::null;
282  qp_F = Teuchos::null;
284  vec_mut_ptr_t qp_f = Teuchos::null;
285  value_type qp_eta = 0.0;
286  VectorMutable &qp_d = pz_k; // pz_k will be updated directly!
287  vec_mut_ptr_t qp_nu = Teuchos::null;
288  vec_mut_ptr_t qp_mu = Teuchos::null;
289  vec_mut_ptr_t qp_Ed = Teuchos::null;
290  vec_mut_ptr_t qp_lambda = Teuchos::null;
291 
292  //
293  // Determine if we can use simple bounds on pz.
294  //
295  // If we have a variable-reduction null-space matrix
296  // (with any choice for Y) then:
297  //
298  // d = Z*pz + (1-eta) * Y*py
299  //
300  // [ d(var_dep) ] = [ D ] * pz + (1-eta) * [ Ypy(var_dep) ]
301  // [ d(var_indep) ] [ I ] [ Ypy(var_indep) ]
302  //
303  // For a cooridinate decomposition (Y = [ I ; 0 ]) then Ypy(var_indep) ==
304  // 0.0 and in this case the bounds on d(var_indep) become simple bounds on
305  // pz even with the relaxation. Also, if dl(var_dep) and du(var_dep) are
306  // unbounded, then we can also use simple bounds since we don't need the
307  // relaxation and we can set eta=0. In this case we just have to subtract
308  // from the upper and lower bounds on pz!
309  //
310  // Otherwise, we can not use simple variable bounds and implement the
311  // relaxation properly.
312  //
313 
314  const MatrixIdentConcat
315  *Zvr = dynamic_cast<const MatrixIdentConcat*>( &Z_k );
316  const Range1D
317  var_dep = Zvr ? Zvr->D_rng() : Range1D::Invalid,
318  var_indep = Zvr ? Zvr->I_rng() : Range1D();
319 
320  RCP<Vector> Ypy_indep;
321  const value_type
322  Ypy_indep_norm_inf
323  = ( m ? (Ypy_indep=Ypy_k->sub_view(var_indep))->norm_inf() : 0.0);
324 
325  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
326  out
327  << "\nDetermine if we can use simple bounds on pz ...\n"
328  << " m = " << m << std::endl
329  << " dynamic_cast<const MatrixIdentConcat*>(&Z_k) = " << Zvr << std::endl
330  << " ||Ypy_k(var_indep)||inf = " << Ypy_indep_norm_inf << std::endl;
331 
332  const bool
333  bounded_var_dep
334  = (
335  m > 0
336  &&
337  num_bounded( *bl->sub_view(var_dep), *bu->sub_view(var_dep), qp_bnd_inf )
338  );
339 
340  const bool
341  use_simple_pz_bounds
342  = (
343  m == 0
344  ||
345  ( Zvr != NULL && ( Ypy_indep_norm_inf == 0.0 || bounded_var_dep == 0 ) )
346  );
347 
348  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
349  out
350  << (use_simple_pz_bounds
351  ? "\nUsing simple bounds on pz ...\n"
352  : "\nUsing bounds on full Z*pz ...\n")
353  << (bounded_var_dep
354  ? "\nThere are finite bounds on dependent variables. Adding extra inequality constrints for D*pz ...\n"
355  : "\nThere are no finite bounds on dependent variables. There will be no extra inequality constraints added on D*pz ...\n" ) ;
356 
357  if( use_simple_pz_bounds ) {
358  // Set simple bound constraints on pz
359  qp_dL = bl->sub_view(var_indep);
360  qp_dU = bu->sub_view(var_indep);
361  qp_nu = nu_k.sub_view(var_indep); // nu_k(var_indep) will be updated directly!
362  if( m && bounded_var_dep ) {
363  // Set general inequality constraints for D*pz
364  qp_E = Teuchos::rcp(&Zvr->D(),false);
365  qp_b = Ypy_k->sub_view(var_dep);
366  qp_eL = bl->sub_view(var_dep);
367  qp_eU = bu->sub_view(var_dep);
368  qp_mu = nu_k.sub_view(var_dep); // nu_k(var_dep) will be updated directly!
369  qp_Ed = Zpz_k.sub_view(var_dep); // Zpz_k(var_dep) will be updated directly!
370  }
371  else {
372  // Leave these as NULL since there is no extra general inequality constraints
373  }
374  }
375  else if( !use_simple_pz_bounds ) { // ToDo: Leave out parts for unbounded dependent variables!
376  // There are no simple bounds! (leave qp_dL, qp_dU and qp_nu as null)
377  // Set general inequality constraints for Z*pz
378  qp_E = Teuchos::rcp(&Z_k,false);
379  qp_b = Teuchos::rcp(Ypy_k,false);
380  qp_eL = bl;
381  qp_eU = bu;
382  qp_mu = Teuchos::rcp(&nu_k,false);
383  qp_Ed = Teuchos::rcp(&Zpz_k,false); // Zpz_k will be updated directly!
384  }
385  else {
387  }
388 
389  // Set the general equality constriants (if they exist)
390  Range1D equ_undecomp = s.equ_undecomp();
391  if( m > r && m > 0 ) {
392  // qp_f = Uy_k * py_k + c_k(equ_undecomp)
393  qp_f = s.space_c().sub_space(equ_undecomp)->create_member();
394  V_MtV( qp_f.get(), Uy_iq->get_k(0), BLAS_Cpp::no_trans, py_iq.get_k(0) );
395  Vp_V( qp_f.get(), *c_iq->get_k(0).sub_view(equ_undecomp) );
396  // Must resize for the undecomposed constriants if it has not already been
397  qp_F = Teuchos::rcp(&Uy_iq->get_k(0),false);
398  qp_lambda = lambda_iq->set_k(0).sub_view(equ_undecomp); // lambda_k(equ_undecomp), will be updated directly!
399  }
400 
401  // Setup the rest of the arguments
402 
403  QPSolverRelaxed::EOutputLevel
404  qp_olevel;
405  switch( olevel ) {
406  case PRINT_NOTHING:
407  qp_olevel = QPSolverRelaxed::PRINT_NONE;
408  break;
409  case PRINT_BASIC_ALGORITHM_INFO:
410  qp_olevel = QPSolverRelaxed::PRINT_NONE;
411  break;
412  case PRINT_ALGORITHM_STEPS:
413  qp_olevel = QPSolverRelaxed::PRINT_BASIC_INFO;
414  break;
415  case PRINT_ACTIVE_SET:
416  qp_olevel = QPSolverRelaxed::PRINT_ITER_SUMMARY;
417  break;
418  case PRINT_VECTORS:
419  qp_olevel = QPSolverRelaxed::PRINT_ITER_VECTORS;
420  break;
421  case PRINT_ITERATION_QUANTITIES:
422  qp_olevel = QPSolverRelaxed::PRINT_EVERY_THING;
423  break;
424  default:
426  }
427  // ToDo: Set print options so that only vectors matrices etc
428  // are only printed in the null space
429 
430  //
431  // Solve the QP
432  //
433  qp_solver().infinite_bound(qp_bnd_inf);
435  solution_type =
436  qp_solver().solve_qp(
437  int(olevel) == int(PRINT_NOTHING) ? NULL : &out
438  ,qp_olevel
439  ,( algo.algo_cntr().check_results()
440  ? QPSolverRelaxed::RUN_TESTS : QPSolverRelaxed::NO_TESTS )
441  ,qp_g, qp_G, qp_etaL, qp_dL.get(), qp_dU.get()
442  ,qp_E.get(), qp_trans_E, qp_E.get() ? qp_b.get() : NULL
443  ,qp_E.get() ? qp_eL.get() : NULL, qp_E.get() ? qp_eU.get() : NULL
444  ,qp_F.get(), qp_trans_F, qp_F.get() ? qp_f.get() : NULL
445  ,NULL // obj_d
446  ,&qp_eta, &qp_d
447  ,qp_nu.get()
448  ,qp_mu.get(), qp_E.get() ? qp_Ed.get() : NULL
449  ,qp_F.get() ? qp_lambda.get() : NULL
450  ,NULL // qp_Fd
451  );
452 
453  //
454  // Check the optimality conditions for the QP
455  //
456  std::ostringstream omsg;
457  bool throw_qp_failure = false;
458  if( qp_testing() == QP_TEST
459  || ( qp_testing() == QP_TEST_DEFAULT && algo.algo_cntr().check_results() ) )
460  {
461  if( int(olevel) >= int(PRINT_ALGORITHM_STEPS) ) {
462  out << "\nChecking the optimality conditions of the reduced QP subproblem ...\n";
463  }
464  if(!qp_tester().check_optimality_conditions(
465  solution_type,qp_solver().infinite_bound()
466  ,int(olevel) == int(PRINT_NOTHING) ? NULL : &out
467  ,int(olevel) >= int(PRINT_VECTORS) ? true : false
468  ,int(olevel) >= int(PRINT_ITERATION_QUANTITIES) ? true : false
469  ,qp_g, qp_G, qp_etaL, qp_dL.get(), qp_dU.get()
470  ,qp_E.get(), qp_trans_E, qp_E.get() ? qp_b.get() : NULL
471  ,qp_E.get() ? qp_eL.get() : NULL, qp_E.get() ? qp_eU.get() : NULL
472  ,qp_F.get(), qp_trans_F, qp_F.get() ? qp_f.get() : NULL
473  ,NULL // obj_d
474  ,&qp_eta, &qp_d
475  ,qp_nu.get()
476  ,qp_mu.get(), qp_E.get() ? qp_Ed.get() : NULL
477  ,qp_F.get() ? qp_lambda.get() : NULL
478  ,NULL // qp_Fd
479  ))
480  {
481  omsg << "\n*** Alert! at least one of the QP optimality conditions did not check out.\n";
482  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
483  out << omsg.str();
484  }
485  throw_qp_failure = true;
486  }
487  }
488 
489  //
490  // Set the solution
491  //
492  if( !use_simple_pz_bounds ) {
493  // Everything is already updated!
494  }
495  else if( use_simple_pz_bounds ) {
496  // Just have to set Zpz_k(var_indep) = pz_k
497  *Zpz_k.sub_view(var_indep) = pz_k;
498  if( m && !bounded_var_dep ) {
499  // Must compute Zpz_k(var_dep) = D*pz
500  LinAlgOpPack::V_MtV( &*Zpz_k.sub_view(var_dep), Zvr->D(), BLAS_Cpp::no_trans, pz_k );
501  // ToDo: Remove the compuation of Zpz here unless you must
502  }
503  }
504  else {
506  }
507 
508  // Set the solution statistics
509  qp_solver_stats_(s).set_k(0) = qp_solver().get_qp_stats();
510 
511  // Cut back Ypy_k = (1-eta) * Ypy_k
512  const value_type eps = std::numeric_limits<value_type>::epsilon();
513  if( fabs(qp_eta - 0.0) > eps ) {
514  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
515  out
516  << "\n*** Alert! the QP was infeasible (eta = "<<qp_eta<<"). Cutting back Ypy_k = (1.0 - eta)*Ypy_k ...\n";
517  }
518  Vt_S( Ypy_k, 1.0 - qp_eta );
519  }
520 
521  // eta_k
522  eta_iq.set_k(0) = qp_eta;
523 
524  //
525  // Modify the solution if we have to!
526  //
527  switch(solution_type) {
528  case QPSolverStats::OPTIMAL_SOLUTION:
529  break; // we are good!
530  case QPSolverStats::PRIMAL_FEASIBLE_POINT:
531  {
532  omsg
533  << "\n*** Alert! the returned QP solution is PRIMAL_FEASIBLE_POINT but not optimal!\n";
534  if( primal_feasible_point_error() )
535  omsg
536  << "\n*** primal_feasible_point_error == true, this is an error!\n";
537  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
538  out << omsg.str();
539  }
540  throw_qp_failure = primal_feasible_point_error();
541  break;
542  }
543  case QPSolverStats::DUAL_FEASIBLE_POINT:
544  {
545  omsg
546  << "\n*** Alert! the returned QP solution is DUAL_FEASIBLE_POINT"
547  << "\n*** but not optimal so we cut back the step ...\n";
548  if( dual_feasible_point_error() )
549  omsg
550  << "\n*** dual_feasible_point_error == true, this is an error!\n";
551  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
552  out << omsg.str();
553  }
554  // Cut back the step to fit in the bounds
555  //
556  // dl <= u*(Ypy_k+Zpz_k) <= du
557  //
558  vec_mut_ptr_t
559  zero = s.space_x().create_member(0.0),
560  d_tmp = s.space_x().create_member();
561  V_VpV( d_tmp.get(), *Ypy_k, Zpz_k );
562  const std::pair<value_type,value_type>
563  u_steps = max_near_feas_step( *zero, *d_tmp, dl_iq.get_k(0), du_iq.get_k(0), 0.0 );
564  const value_type
565  u = my_min( u_steps.first, 1.0 ); // largest positive step size
566  alpha_iq.set_k(0) = u;
567  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
568  out << "\nFinding u s.t. dl <= u*(Ypy_k+Zpz_k) <= du\n"
569  << "max step length u = " << u << std::endl
570  << "alpha_k = u = " << alpha_iq.get_k(0) << std::endl;
571  }
572  throw_qp_failure = dual_feasible_point_error();
573  break;
574  }
575  case QPSolverStats::SUBOPTIMAL_POINT:
576  {
577  omsg
578  << "\n*** Alert!, the returned QP solution is SUBOPTIMAL_POINT!\n";
579  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
580  out << omsg.str();
581  }
582  throw_qp_failure = true;
583  break;
584  }
585  default:
586  TEUCHOS_TEST_FOR_EXCEPT(true); // should not happen!
587  }
588 
589  //
590  // Output the final solution!
591  //
592  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
593  out << "\n||pz_k||inf = " << s.pz().get_k(0).norm_inf()
594  << "\nnu_k.nz() = " << s.nu().get_k(0).nz()
595  << "\nmax(|nu_k(i)|) = " << s.nu().get_k(0).norm_inf()
596 // << "\nmin(|nu_k(i)|) = " << min_abs( s.nu().get_k(0)() )
597  ;
598  if( m > r ) out << "\n||lambda_k(undecomp)||inf = " << s.lambda().get_k(0).norm_inf();
599  out << "\n||Zpz_k||2 = " << s.Zpz().get_k(0).norm_2()
600  ;
601  if(qp_eta > 0.0) out << "\n||Ypy||2 = " << s.Ypy().get_k(0).norm_2();
602  out << std::endl;
603  }
604 
605  if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
606  out << "\npz_k = \n" << s.pz().get_k(0);
607  if(var_indep.size())
608  out << "\nnu_k(var_indep) = \n" << *s.nu().get_k(0).sub_view(var_indep);
609  }
610 
611  if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
612  if(var_indep.size())
613  out << "\nZpz(var_indep)_k = \n" << *s.Zpz().get_k(0).sub_view(var_indep);
614  out << std::endl;
615  }
616 
617  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
618  if(var_dep.size())
619  out << "\nZpz(var_dep)_k = \n" << *s.Zpz().get_k(0).sub_view(var_dep);
620  out << "\nZpz_k = \n" << s.Zpz().get_k(0);
621  out << std::endl;
622  }
623 
624  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
625  out << "\nnu_k = \n" << s.nu().get_k(0);
626  if(var_dep.size())
627  out << "\nnu_k(var_dep) = \n" << *s.nu().get_k(0).sub_view(var_dep);
628  if( m > r )
629  out << "\nlambda_k(equ_undecomp) = \n" << *s.lambda().get_k(0).sub_view(equ_undecomp);
630  if(qp_eta > 0.0) out << "\nYpy = \n" << s.Ypy().get_k(0);
631  }
632 
633  if( qp_eta == 1.0 ) {
634  omsg
635  << "TangentialStepWithInequStd_Step::do_step(...) : Error, a QP relaxation parameter\n"
636  << "of eta = " << qp_eta << " was calculated and therefore it must be assumed\n"
637  << "that the NLP's constraints are infeasible\n"
638  << "Throwing an InfeasibleConstraints exception!\n";
639  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
640  out << omsg.str();
641  }
642  throw InfeasibleConstraints(omsg.str());
643  }
644 
645  if( throw_qp_failure )
646  throw QPFailure( omsg.str(), qp_solver().get_qp_stats() );
647 
648  return true;
649 }
650 
652  const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type
653  ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
654  ) const
655 {
656  out
657  << L << "*** Calculate the null-space step by solving a constrained QP\n"
658  << L << "qp_grad_k = rGf_k\n"
659  << L << "if w_k is updated then set qp_grad_k = qp_grad_k + zeta_k * w_k\n"
660  << L << "bl = dl_k - Ypy_k\n"
661  << L << "bu = du_k - Ypy_k\n"
662  << L << "etaL = 0.0\n"
663  << L << "*** Determine if we can use simple bounds on pz or not\n"
664  << L << "if num_bounded(bl_k(var_dep),bu_k(var_dep)) > 0 then\n"
665  << L << " bounded_var_dep = true\n"
666  << L << "else\n"
667  << L << " bounded_var_dep = false\n"
668  << L << "end\n"
669  << L << "if( m==0\n"
670  << L << " or\n"
671  << L << " ( Z_k is a variable reduction null space matrix\n"
672  << L << " and\n"
673  << L << " ( ||Ypy_k(var_indep)||inf == 0 or bounded_var_dep==false ) )\n"
674  << L << " ) then\n"
675  << L << " use_simple_pz_bounds = true\n"
676  << L << "else\n"
677  << L << " use_simple_pz_bounds = false\n"
678  << L << "end\n"
679  << L << "*** Setup QP arguments\n"
680  << L << "qp_g = qp_grad_k\n"
681  << L << "qp_G = rHL_k\n"
682  << L << "if (use_simple_pz_bounds == true) then\n"
683  << L << " qp_dL = bl(var_indep), qp_dU = bu(var_indep))\n"
684  << L << " if (m > 0) then\n"
685  << L << " qp_E = Z_k.D, qp_b = Ypy_k(var_dep)\n"
686  << L << " qp_eL = bl(var_dep), qp_eU = bu(var_dep)\n"
687  << L << " end\n"
688  << L << "elseif (use_simple_pz_bounds == false) then\n"
689  << L << " qp_dL = -inf, qp_dU = +inf\n"
690  << L << " qp_E = Z_k, qp_b = Ypy_k\n"
691  << L << " qp_eL = bl, qp_eU = bu\n"
692  << L << "end\n"
693  << L << "if (m > r) then\n"
694  << L << " qp_F = V_k, qp_f = Uy_k*py_k + c_k(equ_undecomp)\n"
695  << L << "else\n"
696  << L << " qp_F = empty, qp_f = empty\n"
697  << L << "end\n"
698  << L << "Given active-set statistics (act_set_stats_km1)\n"
699  << L << " frac_same = max(num_active-num_adds-num_drops,0)/(num_active)\n"
700  << L << "Use a warm start when frac_same >= warm_start_frac\n"
701  << L << "Solve the following QP to compute qp_d, qp_eta, qp_Ed = qp_E * qp_d\n"
702  << L << ",qp_nu, qp_mu and qp_lambda (" << typeName(qp_solver()) << "):\n"
703  << L << " min qp_g' * qp_d + 1/2 * qp_d' * qp_G * qp_d + M(eta)\n"
704  << L << " qp_d <: R^(n-r)\n"
705  << L << " s.t.\n"
706  << L << " etaL <= eta\n"
707  << L << " qp_dL <= qp_d <= qp_dU [qp_nu]\n"
708  << L << " qp_eL <= qp_E * qp_d + (1-eta)*qp_b <= qp_eU [qp_mu]\n"
709  << L << " qp_F * d_qp + (1-eta) * qp_f = 0 [qp_lambda]\n"
710  << L << "if (qp_testing==QP_TEST) or (fd_deriv_testing==QP_TEST_DEFAULT\n"
711  << L << "and check_results==true) then\n"
712  << L << " Check the optimality conditions of the above QP\n"
713  << L << " if the optimality conditions do not check out then\n"
714  << L << " set throw_qp_failure = true\n"
715  << L << " end\n"
716  << L << "end\n"
717  << L << "*** Set the solution to the QP subproblem\n"
718  << L << "pz_k = qp_d\n"
719  << L << "eta_k = qp_eta\n"
720  << L << "if (use_simple_pz_bounds == true) then\n"
721  << L << " nu_k(var_dep) = qp_mu, nu_k(var_indep) = qp_nu\n"
722  << L << " Zpz_k(var_dep) = qp_Ed, Zpz_k(var_indep) = pz_k\n"
723  << L << "elseif (use_simple_pz_bounds == false) then\n"
724  << L << " nu_k = qp_mu\n"
725  << L << " Zpz_k = qp_Ed\n"
726  << L << "end\n"
727  << L << "if m > r then\n"
728  << L << " lambda_k(equ_undecomp) = qp_lambda\n"
729  << L << "end\n"
730  << L << "if (eta_k > 0) then set Ypy_k = (1-eta_k) * Ypy_k\n"
731  << L << "if QP solution is suboptimal then\n"
732  << L << " throw_qp_failure = true\n"
733  << L << "elseif QP solution is primal feasible (not optimal) then\n"
734  << L << " throw_qp_failure = primal_feasible_point_error\n"
735  << L << "elseif QP solution is dual feasible (not optimal) then\n"
736  << L << " find max u s.t.\n"
737  << L << " dl_k <= u*(Ypy_k+Zpz_k) <= du_k\n"
738  << L << " alpha_k = u\n"
739  << L << " throw_qp_failure = dual_feasible_point_error\n"
740  << L << "end\n"
741  << L << "if (eta_k == 1.0) then\n"
742  << L << " The constraints are infeasible!\n"
743  << L << " throw InfeasibleConstraints(...)\n"
744  << L << "end\n"
745  << L << "if (throw_qp_failure == true) then\n"
746  << L << " throw QPFailure(...)\n"
747  << L << "end\n"
748  ;
749 }
750 
751 } // end namespace MoochoPack
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
T_To & dyn_cast(T_From &from)
Thrown if a the QP failed and was not corrected.
rSQP Algorithm control class.
T * get() const
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
void print_step(const Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss, std::ostream &out, const std::string &leading_str) const
Class for storing statistics about the changes in the active set of an SQP algorithm.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual std::ostream & journal_out() const
void V_VpV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
Reduced space SQP state encapsulation interface.
size_t size_type
AlgorithmTracker & track()
size_type num_bounded(const Vector &xl, const Vector &xu, value_type inf_bound)
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
std::pair< value_type, value_type > max_near_feas_step(const Vector &x, const Vector &d, const Vector &xl, const Vector &xu, value_type max_bnd_viol)
NLPAlgoState & rsqp_state()
<<std aggr>="">> members for algo_cntr
Transp
void V_VmV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
std::string typeName(const T &t)
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)