MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_FeasibilityStepReducedStd_Strategy.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 "MoochoPack_FeasibilityStepReducedStd_Strategy.hpp"
43 #include "MoochoPack_NLPAlgo.hpp"
44 #include "MoochoPack_NLPAlgoState.hpp"
45 #include "MoochoPack_Exceptions.hpp"
46 #include "Teuchos_dyn_cast.hpp"
47 
48 namespace MoochoPack {
49 
51  const quasi_range_space_step_ptr_t &quasi_range_space_step
52  ,const qp_solver_ptr_t &qp_solver
53  ,const qp_tester_ptr_t &qp_tester
54  ,EQPObjective qp_objective
55  ,EQPTesting qp_testing
56  )
57  :quasi_range_space_step_(quasi_range_space_step)
58  ,qp_solver_(qp_solver)
59  ,qp_tester_(qp_tester)
60  ,qp_objective_(qp_objective)
61  ,qp_testing_(qp_testing)
62  ,dl_iq_(dl_name)
63  ,du_iq_(du_name)
64  ,current_k_(-1)
65 {}
66 
68  std::ostream& out, EJournalOutputLevel olevel, NLPAlgo *algo, NLPAlgoState *s
69  ,const Vector& xo, const Vector& c_xo, VectorMutable* w
70  )
71 {
72  using Teuchos::dyn_cast;
73 
74 /* Todo: UPdate below code!
75 
76  // problem dimensions
77  const size_type
78  n = algo->nlp().n(),
79  m = algo->nlp().m(),
80  r = s->equ_decomp().size();
81 
82  // Compute the quasi-range space step Ywy
83  Workspace<value_type> Ywy_ws(wss,xo.size());
84  DVectorSlice Ywy(&Ywy_ws[0],Ywy_ws.size());
85  if(!quasi_range_space_step().solve_quasi_range_space_step(
86  out,olevel,algo,s,xo,c_xo,&Ywy ))
87  return false;
88 
89  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
90  out << "\n||Ywy||2 = " << DenseLinAlgPack::norm_2(Ywy);
91  out << std::endl;
92  }
93  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
94  out << "\nYwy = \n" << Ywy;
95  }
96 
97  //
98  // Set the bounds on the null space QP subproblem:
99  //
100  // d_bounds_k.l <= (xo - x_k) + (1-eta) * Ywy + Z*wz <= d_bounds_k.u
101  // =>
102  // bl <= Z*wz - eta * Ywy <= bu
103  //
104  // where: bl = d_bounds_k.l - (xo - x_k) - Ywy
105  // bu = d_bounds_k.u - (xo - x_k) - Ywy
106  //
107  // Above, fix the variables that are at an active bound as equalities
108  // to maintain the same active-set.
109  //
110  const SparseBounds
111  &d_bounds = d_bounds_(*s).get_k(0);
112  const SpVectorSlice
113  &dl = d_bounds.l,
114  &du = d_bounds.u;
115  const DVector
116  &x_k = s->x().get_k(0).v();
117  const SpVector
118  &nu_k = s->nu().get_k(0);
119  TEUCHOS_TEST_FOR_EXCEPT( !( nu_k.is_sorted() ) );
120  SpVector bl(n,n), bu(n,n);
121  sparse_bounds_itr
122  d_bnds_itr(dl.begin(),dl.end(),dl.offset(),du.begin(),du.end(),du.offset());
123  SpVector::const_iterator
124  nu_itr = nu_k.begin(),
125  nu_end = nu_k.end();
126  for( ; !d_bnds_itr.at_end(); ++d_bnds_itr ) {
127  typedef SpVectorSlice::element_type ele_t;
128  const size_type i = d_bnds_itr.indice();
129  while( nu_itr != nu_end && nu_itr->indice() + nu_k.offset() < i )
130  ++nu_itr;
131  if( nu_itr != nu_end && nu_itr->indice() + nu_k.offset() == i ) {
132  const value_type
133  act_bnd = nu_itr->value() > 0.0 ? d_bnds_itr.ubound() : d_bnds_itr.lbound();
134  bl.add_element(ele_t( i, act_bnd - xo(i) + x_k(i) - Ywy(i) ));
135  bu.add_element(ele_t( i, act_bnd - xo(i) + x_k(i) - Ywy(i) ));
136  }
137  else {
138  if( d_bnds_itr.lbound() != -d_bnds_itr.big_bnd() )
139  bl.add_element(ele_t(i,d_bnds_itr.lbound() - xo(i) + x_k(i) - Ywy(i) ));
140  if( d_bnds_itr.ubound() != +d_bnds_itr.big_bnd() )
141  bu.add_element(ele_t(i, d_bnds_itr.ubound() - xo(i) + x_k(i) - Ywy(i) ));
142  }
143  }
144  bl.assume_sorted(true);
145  bu.assume_sorted(true);
146  //
147  // Setup the objective function for the null space QP subproblem
148  //
149  //
150  // OBJ_MIN_FULL_STEP
151  // min 1/2 * (Y*wy + Z*wz)'*(Y*wy + Z*wz) = 1/2*wy'*Y'*Y*wy + (Z'*Y*wy)'*wz + 1/2*wz'*(Z'*Z)*wz
152  // => grad = (Z'*Y*wy), Hess = Z'*Z
153  //
154  // OBJ_MIN_WZ
155  // min 1/2 * wz'*wz => grad = 0, Hess = I
156  //
157  // OBJ_RSQP
158  // min qp_grad_k'*wz + 1/2 * wz'*rHL_k*wz
159  // => grad = qp_grad, Hess = rHL_k
160  //
161  const MatrixOp
162  &Z_k = s->Z().get_k(0);
163  if( current_k_ != s->k() ) {
164  if( qp_objective() != OBJ_RSQP )
165  grad_store_.resize(n-r);
166  if( qp_objective() == OBJ_MIN_FULL_STEP )
167  Hess_store_.resize(n-r+1,n-r+1);
168  }
169  DVectorSlice grad;
170  switch(qp_objective())
171  {
172  case OBJ_MIN_FULL_STEP: // grad = (Z'*Ywy), Hess = Z'*Z
173  {
174  grad.bind( grad_store_() );
175  if( current_k_ != s->k() ) {
176  // grad = (Z'*Ywy)
177  LinAlgOpPack::V_MtV( &grad, Z_k, BLAS_Cpp::trans, Ywy );
178  // Hess = Z'*Z
179  DMatrixSliceSym S(Hess_store_(2,n-r+1,1,n-r),BLAS_Cpp::lower); // Must be strictly lower triangular here!
180  Z_k.syrk( BLAS_Cpp::trans, 1.0, 0.0, &S ); // S = 1.0*Z'*Z + 0.0*S
181  MatrixSymPosDefCholFactor
182  *H_ptr = NULL;
183  if( Hess_ptr_.get() == NULL || dynamic_cast<const MatrixSymPosDefCholFactor*>(Hess_ptr_.get()) == NULL )
184  Hess_ptr_ = new MatrixSymPosDefCholFactor;
185  H_ptr = const_cast<MatrixSymPosDefCholFactor*>(dynamic_cast<const MatrixSymPosDefCholFactor*>(Hess_ptr_.get()));
186  TEUCHOS_TEST_FOR_EXCEPT( !( H_ptr ) ); // Should not be null!
187  H_ptr->init_setup(
188  &Hess_store_() // The original matrix is stored in the lower triangular part (below diagonal)!
189  ,NULL // Nothing to deallocate
190  ,n-r
191  ,true // maintain the original factor
192  ,false // don't maintain the factor
193  ,true // allow the factor to be computed if needed
194  ,true // Set the view
195  ,1.0 // Scale the matrix by one
196  );
197  }
198  break;
199  }
200  case OBJ_MIN_NULL_SPACE_STEP: // grad = 0, Hess = I
201  {
202  grad.bind( grad_store_() );
203  MatrixSymIdent
204  *H_ptr = NULL;
205  if( Hess_ptr_.get() == NULL || dynamic_cast<const MatrixSymIdent*>(Hess_ptr_.get()) == NULL )
206  Hess_ptr_ = new MatrixSymIdent;
207  if( current_k_ != s->k() ) {
208  H_ptr = const_cast<MatrixSymIdent*>(dynamic_cast<const MatrixSymIdent*>(Hess_ptr_.get()));
209  TEUCHOS_TEST_FOR_EXCEPT( !( H_ptr ) ); // Should not be null!
210  H_ptr->init_setup(n-r,1.0);
211  grad = 0.0;
212  }
213  break;
214  }
215  case OBJ_RSQP: // grad = qp_grad, Hess = rHL_k
216  {
217  grad.bind( s->qp_grad().get_k(0)() );
218  Hess_ptr_ = Hess_ptr_t( &s->rHL().get_k(0), false ); // don't delete memory!
219  break;
220  }
221  defaut:
222  TEUCHOS_TEST_FOR_EXCEPT(true); // Not a valid option
223  }
224 
225  //
226  // Solve the null space subproblem
227  //
228 
229  Workspace<value_type> wz_ws(wss,n-r),Zwz_ws(wss,n);
230  DVectorSlice wz(&wz_ws[0],wz_ws.size());
231  DVectorSlice Zwz(&Zwz_ws[0],Zwz_ws.size());
232  value_type qp_eta = 0;
233 
234  bool throw_qp_failure = false;
235 
236  if( bl.nz() == 0 && bu.nz() == 0 && m-r == 0 ) {
237  //
238  // Just solve the unconstrainted QP
239  //
240  // wz = - inv(Hess)*grad
241 #ifdef _WINDOWS
242  const MatrixFactorized &Hess = dynamic_cast<const MatrixFactorized&>(*Hess_ptr_);
243 #else
244  const MatrixFactorized &Hess = dyn_cast<const MatrixFactorized>(*Hess_ptr_);
245 #endif
246  AbstractLinAlgPack::V_InvMtV( &wz, Hess, BLAS_Cpp::no_trans, grad );
247  DenseLinAlgPack::Vt_S(&wz,-1.0);
248  // Zwz = Z*wz
249  LinAlgOpPack::V_MtV( &Zwz, Z_k, BLAS_Cpp::no_trans, wz );
250  }
251  else {
252 
253  //
254  // Set the arguments to the QP subproblem
255  //
256 
257  DVectorSlice qp_g = grad;
258  const MatrixOp& qp_G = *Hess_ptr_;
259  const value_type qp_etaL = 0.0;
260  SpVectorSlice qp_dL(NULL,0,0,n-r); // If nz() == 0 then no simple bounds
261  SpVectorSlice qp_dU(NULL,0,0,n-r);
262  const MatrixOp *qp_E = NULL;
263  BLAS_Cpp::Transp qp_trans_E = BLAS_Cpp::no_trans;
264  DVectorSlice qp_b;
265  SpVectorSlice qp_eL(NULL,0,0,n);
266  SpVectorSlice qp_eU(NULL,0,0,n);
267  const MatrixOp *qp_F = NULL;
268  BLAS_Cpp::Transp qp_trans_F = BLAS_Cpp::no_trans;
269  DVectorSlice qp_f;
270  DVectorSlice qp_d = wz;
271  SpVector *qp_nu = NULL;
272  SpVector *qp_mu = NULL;
273  DVectorSlice qp_Ed;
274  DVectorSlice qp_lambda;
275 
276  SpVector _nu_wz, _nu_Dwz, // Possible storage for multiplers for separate inequality
277  _nu; // constriants for wz.
278  DVector _Dwz; // Possible storage for D*wz computed by QP solver?
279 
280  //
281  // Determine if we can use simple bounds on wz.
282  //
283  // If we have a variable reduction null space matrix
284  // (with any choice for Y) then:
285  //
286  // w = Z*wz + (1-eta) * Y*wy
287  //
288  // [ w(dep) ] = [ D ] * wz + (1-eta) * [ Ywy(dep) ]
289  // [ w(indep) ] [ I ] [ Ywy(indep) ]
290  //
291  // For a cooridinate decomposition (Y = [ I ; 0 ]) then Ywy(indep) = 0 and
292  // in this case the bounds on d(indep) become simple bounds on pz even
293  // with the relaxation.
294  //
295  // Otherwise, we can not use simple variable bounds and implement the
296  // relaxation properly.
297  //
298 
299  const ZVarReductMatrix
300  *Zvr = dynamic_cast<const ZVarReductMatrix*>( &Z_k );
301  Range1D
302  indep = Zvr ? Zvr->indep() : Range1D(),
303  dep = Zvr ? Zvr->dep() : Range1D();
304 
305  const bool
306  use_simple_wz_bounds = ( Zvr!=NULL && DenseLinAlgPack::norm_inf(Ywy(indep))==0.0 );
307 
308  if( use_simple_wz_bounds ) {
309 
310  // Set simple bound constraints on pz
311  qp_dL.bind( bl(indep) );
312  qp_dU.bind( bu(indep) );
313  qp_nu = &( _nu_wz = s->nu().get_k(0)(indep) ); // warm start?
314 
315  // Set general inequality constraints for D*pz
316  qp_E = &Zvr->D();
317  qp_b.bind( Ywy(dep) );
318  qp_eL.bind( bl(dep) );
319  qp_eU.bind( bu(dep) );
320  qp_mu = &( _nu_Dwz = s->nu().get_k(0)(dep) ); // warm start?
321  _Dwz.resize(r);
322  qp_Ed.bind(_Dwz()); // Part of Zwz will be updated directly!
323 
324  }
325  else {
326 
327  // Set general inequality constraints for Z*pz
328  qp_E = &Z_k;
329  qp_b.bind( Ywy() );
330  qp_eL.bind( bl() );
331  qp_eU.bind( bu() );
332  qp_mu = &(_nu = s->nu().get_k(0)); // warm start??
333  qp_Ed.bind(Zwz); // Zwz
334  }
335 
336  // Set the general equality constriants (if they exist)
337  DVector q(m-r);
338  Range1D undecomp = s->con_undecomp();
339  if( m > r ) {
340  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement when needed!
341  }
342 
343  // Setup the rest of the arguments
344 
345  QPSolverRelaxed::EOutputLevel
346  qp_olevel;
347  switch( olevel ) {
348  case PRINT_NOTHING:
349  qp_olevel = QPSolverRelaxed::PRINT_NONE;
350  break;
351  case PRINT_BASIC_ALGORITHM_INFO:
352  qp_olevel = QPSolverRelaxed::PRINT_BASIC_INFO;
353  break;
354  case PRINT_ALGORITHM_STEPS:
355  qp_olevel = QPSolverRelaxed::PRINT_BASIC_INFO;
356  break;
357  case PRINT_ACTIVE_SET:
358  qp_olevel = QPSolverRelaxed::PRINT_ITER_SUMMARY;
359  break;
360  case PRINT_VECTORS:
361  qp_olevel = QPSolverRelaxed::PRINT_ITER_VECTORS;
362  break;
363  case PRINT_ITERATION_QUANTITIES:
364  qp_olevel = QPSolverRelaxed::PRINT_EVERY_THING;
365  break;
366  default:
367  TEUCHOS_TEST_FOR_EXCEPT(true);
368  }
369 
370  //
371  // Solve the QP
372  //
373  const QPSolverStats::ESolutionType
374  solution_type =
375  qp_solver().solve_qp(
376  int(olevel) == int(PRINT_NOTHING) ? NULL : &out
377  , qp_olevel
378  , algo->algo_cntr().check_results()
379  ? QPSolverRelaxed::RUN_TESTS : QPSolverRelaxed::NO_TESTS
380  , qp_g, qp_G, qp_etaL, qp_dL, qp_dU
381  , qp_E, qp_trans_E, qp_E ? &qp_b : NULL
382  , qp_E ? &qp_eL : NULL, qp_E ? &qp_eU : NULL
383  , qp_F, qp_trans_F, qp_F ? &qp_f : NULL
384  , NULL
385  , &qp_eta, &qp_d
386  , qp_nu
387  , qp_mu, qp_E ? &qp_Ed : NULL
388  , qp_F ? &qp_lambda : NULL, NULL
389  );
390 
391  //
392  // Check the optimality conditions for the QP
393  //
394  std::ostringstream omsg;
395  if( qp_testing() == QP_TEST
396  || ( qp_testing() == QP_TEST_DEFAULT && algo->algo_cntr().check_results() ) )
397  {
398  if( int(olevel) >= int(PRINT_ALGORITHM_STEPS) ) {
399  out << "\nChecking the optimality conditions of the reduced QP subproblem ...\n";
400  }
401  if(!qp_tester().check_optimality_conditions(
402  solution_type
403  , int(olevel) == int(PRINT_NOTHING) ? NULL : &out
404  , int(olevel) >= int(PRINT_VECTORS) ? true : false
405  , int(olevel) >= int(PRINT_ITERATION_QUANTITIES) ? true : false
406  , qp_g, qp_G, qp_etaL, qp_dL, qp_dU
407  , qp_E, qp_trans_E, qp_E ? &qp_b : NULL
408  , qp_E ? &qp_eL : NULL, qp_E ? &qp_eU : NULL
409  , qp_F, qp_trans_F, qp_F ? &qp_f : NULL
410  , NULL
411  , &qp_eta, &qp_d
412  , qp_nu
413  , qp_mu, qp_E ? &qp_Ed : NULL
414  , qp_F ? &qp_lambda : NULL, NULL
415  ))
416  {
417  omsg << "\n*** Alert! at least one of the QP optimality conditions did not check out.\n";
418  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
419  out << omsg.str();
420  }
421  throw_qp_failure = true;
422  }
423  }
424 
425  if( solution_type != QPSolverStats::OPTIMAL_SOLUTION ) {
426  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
427  out << "\nCould not solve the QP!\n";
428  }
429  return false;
430  }
431 
432  //
433  // Set the solution
434  //
435  if( use_simple_wz_bounds ) {
436  // Set Zwz
437  Zwz(dep) = _Dwz;
438  Zwz(indep) = wz;
439  }
440  else {
441  // Everything should already be set!
442  }
443 
444  // Cut back Ywy = (1-eta) * Ywy
445  const value_type eps = std::numeric_limits<value_type>::epsilon();
446  if( fabs(qp_eta - 0.0) > eps ) {
447  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
448  out
449  << "\n*** Alert! the QP was infeasible (eta = "<<qp_eta<<"). Cutting back Ywy_k = (1.0 - eta)*Ywy ...\n";
450  }
451  DenseLinAlgPack::Vt_S( &Ywy , 1.0 - qp_eta );
452  }
453  }
454 
455  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
456  out << "\n||wz||inf = " << DenseLinAlgPack::norm_inf(wz);
457  out << "\n||Zwz||2 = " << DenseLinAlgPack::norm_2(Zwz);
458  if(qp_eta > 0.0) out << "\n||Ypy||2 = " << DenseLinAlgPack::norm_2(Ywy);
459  out << std::endl;
460  }
461  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
462  out << "\nwz = \n" << wz;
463  out << "\nZwz = \n" << Zwz;
464  if(qp_eta > 0.0) out << "\nYwy = \n" << Ywy;
465  }
466  if( qp_eta == 1.0 ) {
467  std::ostringstream omsg;
468  omsg
469  << "FeasibilityStepReducedStd_Strategy::compute_feasibility_step(...) : "
470  << "Error, a QP relaxation parameter\n"
471  << "of eta = " << qp_eta << " was calculated and therefore it must be assumed\n"
472  << "that the NLP's constraints are infeasible\n"
473  << "Throwing an InfeasibleConstraints exception!\n";
474  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
475  out << omsg.str();
476  }
477  throw_qp_failure = true;
478 // throw InfeasibleConstraints(omsg.str());
479  }
480 
481  //
482  // Set the full step
483  //
484  // w = Ywy + Zwz
485  //
486  DenseLinAlgPack::V_VpV( w, Ywy, Zwz );
487 
488  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
489  out << "\n||w||inf = " << DenseLinAlgPack::norm_inf(*w);
490  out << std::endl;
491  }
492  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
493  out << "\nw = \n" << *w;
494  }
495 
496  current_k_ = s->k();
497 
498  if( throw_qp_failure )
499  return false;
500 
501 */
503 
504  return true;
505 }
506 
507 void FeasibilityStepReducedStd_Strategy::print_step( std::ostream& out, const std::string& L ) const
508 {
509  out << L << "*** Computes feasibility steps by solving a constrained QP using the range and null\n"
510  << L << "*** space decomposition\n"
511  << L << "begin quais-range space step: \"" << typeName(quasi_range_space_step()) << "\"\n";
512 
513  quasi_range_space_step().print_step( out, L + " " );
514 
515  out << L << "end quasi-range space step\n"
516  << L << "if quasi-range space step failed then\n"
517  << L << " this feasibility step fails also!\n"
518  << L << "end\n"
519  << L << "Ywy = v\n"
520  << L << "*** Set the bounds for bl <= Z*wz <= bu\n"
521  << L << "bl = d_bounds_k.l - (xo - x_k) - Ywy\n"
522  << L << "bu = d_bounds_k.u - (xo - x_k) - Ywy\n"
523  << L << "Set bl(i) = bu(i) for those nu_k(i) != 0.0\n"
524  << L << "if (qp_objective == OBJ_MIN_FULL_STEP) and (current_k != k) then\n"
525  << L << " grad = Z_k'*Ywy\n"
526  << L << " Hess = Z_k'*Z_k\n"
527  << L << "elseif (qp_objective == OBJ_MIN_NULL_SPACE_STEP) and (current_k != k) then\n"
528  << L << " grad = 0\n"
529  << L << " Hess = I\n"
530  << L << "elseif (qp_objective == OBJ_RSQP) and (current_k != k) then\n"
531  << L << " grad = qp_grad_k\n"
532  << L << " Hess = rHL_k\n"
533  << L << "end\n"
534  << L << "if check_results == true then\n"
535  << L << " assert that bl and bu are valid and sorted\n"
536  << L << "end\n"
537  << L << "etaL = 0.0\n"
538  << L << "*** Determine if we can use simple bounds on pz or not\n"
539  << L << "if Z_k is a variable reduction null space matrix and norm(Ypy_k(indep),0) == 0 then\n"
540  << L << " use_simple_wz_bounds = true\n"
541  << L << "else\n"
542  << L << " use_simple_wz_bounds = false\n"
543  << L << "end\n"
544  << L << "*** Setup QP arguments\n"
545  << L << "qp_g = qp_grad_k\n"
546  << L << "qp_G = rHL_k\n"
547  << L << "if use_simple_wz_bounds == true then\n"
548  << L << " qp_dL = bl(indep), qp_dU = bu(indep)\n"
549  << L << " qp_E = Z_k.D, qp_b = Ywy(dep)\n"
550  << L << " qp_eL = bl(dep), qp_eU = bu(dep)\n"
551  << L << "else\n"
552  << L << " qp_dL = -inf, qp_dU = +inf\n"
553  << L << " qp_E = Z_k, qp_b = Ywy\n"
554  << L << " qp_eL = bl, qp_eU = bu\n"
555  << L << "end\n"
556  << L << "if m > r then\n"
557  << L << " qp_F = V_k, qp_f = c_k(undecomp) + Gc_k(undecomp)'*Ywy\n"
558  << L << "else\n"
559  << L << " qp_F = empty, qp_f = empty\n"
560  << L << "end\n"
561  << L << "Use a warm start given the active-set in nu_k\n"
562  << L << "Solve the following QP to compute qp_d, qp_eta, qp_Ed = qp_E * qp_d\n"
563  << L << ",qp_nu, qp_mu and qp_lambda (" << typeName(qp_solver()) << "):\n"
564  << L << " min qp_g' * qp_d + 1/2 * qp_d' * qp_G * qp_d + M(eta)\n"
565  << L << " qp_d <: R^(n-r)\n"
566  << L << " s.t.\n"
567  << L << " etaL <= qp_eta\n"
568  << L << " qp_dL <= qp_d <= qp_dU [qp_nu]\n"
569  << L << " qp_eL <= qp_E * qp_d + (1-eta)*qp_b <= qp_eU [qp_mu]\n"
570  << L << " qp_F * d_qp + (1-eta) * qp_f = 0 [qp_lambda]\n"
571  << L << "if (qp_teing==QP_TEST) or (fd_deriv_testing==QP_TEST_DEFAULT\n"
572  << L << "and check_results==true) then\n"
573  << L << " Check the optimality conditions of the above QP\n"
574  << L << " if the optimality conditions do not check out then\n"
575  << L << " set throw_qp_failure = true\n"
576  << L << " end\n"
577  << L << "end\n"
578  << L << "*** Set the solution to the QP subproblem\n"
579  << L << "wz = qp_d\n"
580  << L << "eta = qp_eta\n"
581  << L << "if use_simple_wz_bounds == true then\n"
582  << L << " Zwz(dep) = qp_Ed, Zwz(indep) = wz\n"
583  << L << "else\n"
584  << L << " Zwz = qp_Ed\n"
585  << L << "end\n"
586  << L << "if eta > 0 then set Ywy = (1-eta) * Ywy\n"
587  << L << "if QP solution is suboptimal then\n"
588  << L << " throw_qp_failure = true\n"
589  << L << "elseif QP solution is primal feasible (not optimal) then\n"
590  << L << " throw_qp_failure = primal_feasible_point_error\n"
591  << L << "elseif QP solution is dual feasible (not optimal) then\n"
592  << L << " find max u s.t.\n"
593  << L << " d_bounds_k.l <= (xo - x) + u*(Ywy+Zwz) <= d_bounds_k.u\n"
594  << L << " alpha_k = u\n"
595  << L << " throw_qp_failure = true\n"
596  << L << "end\n"
597  << L << "if eta == 1.0 then\n"
598  << L << " The constraints are infeasible!\n"
599  << L << " throw_qp_failure = true\n"
600  << L << "end\n"
601  << L << "current_k = k\n"
602  << L << "w = Zwz + Ywy\n"
603  << L << "if (throw_qp_failure == true) then\n"
604  << L << " The feasibility step computation has failed!\n"
605  << L << "end\n"
606  ;
607 }
608 
609 } // end namespace MoochoPack
void print_step(std::ostream &out, const std::string &leading_str) const
T_To & dyn_cast(T_From &from)
rSQP Algorithm control class.
Reduced space SQP state encapsulation interface.
bool compute_feasibility_step(std::ostream &out, EJournalOutputLevel olevel, NLPAlgo *algo, NLPAlgoState *s, const Vector &xo, const Vector &c_xo, VectorMutable *w)
Computes a feasibility step by computing simple quasi-range and null space components.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
std::string typeName(const T &t)
FeasibilityStepReducedStd_Strategy(const quasi_range_space_step_ptr_t &quasi_range_space_step, const qp_solver_ptr_t &qp_solver, const qp_tester_ptr_t &qp_tester, EQPObjective qp_objective=OBJ_MIN_NULL_SPACE_STEP, EQPTesting qp_testing=QP_TEST_DEFAULT)
Construct and initialize.