24 #ifndef REDUCED_SNOPT__ALGORITHM_H_
25 #define REDUCED_SNOPT__ALGORITHM_H_
53 template<
typename PROBLEM,
typename VECTOR>
91 int rsa_func_(DOpEWrapper::SNOPT_FUNC_DATA& data);
93 std::string postindex_;
96 double func_prec_, feas_tol_, opt_tol_;
97 int max_inner_iter_, max_outer_iter_;
104 using namespace dealii;
108 template<
typename PROBLEM,
typename VECTOR>
113 param_reader.
SetSubsection(
"reduced_snoptalgorithm parameters");
115 Patterns::Double(0, 1),
116 "Declares how many digits we assume to have computed correctly, this should correspond to the tolerance used for the PDE solve");
118 Patterns::Double(0, 1),
119 "Tolerance with respect to the feasibility of the constraints.");
121 Patterns::Double(0, 1),
122 "Tolerance with respect to the optimality condition.");
124 Patterns::Integer(0),
125 "Maximal allowed number of inner iterations over all outer iterations");
126 param_reader.
declare_entry(
"max iterations",
"1000", Patterns::Integer(0),
127 "Maximal allowed number of outer iterations over all outer iterations");
130 "Select if the snopt output should be stored in log file");
136 template<
typename PROBLEM,
typename VECTOR>
146 param_reader.
SetSubsection(
"reduced_snoptalgorithm parameters");
148 func_prec_ = param_reader.
get_double(
"function precision");
149 feas_tol_ = param_reader.
get_double(
"feasibility tol");
150 opt_tol_ = param_reader.
get_double(
"optimality tol");
151 max_inner_iter_ = param_reader.
get_integer(
"max inner iterations");
152 max_outer_iter_ = param_reader.
get_integer(
"max iterations");
153 capture_out_ = param_reader.
get_bool(
"capture snopt output");
155 vector_behavior_ = vector_behavior;
156 postindex_ =
"_" + this->
GetProblem()->GetName();
162 "Reduced_SnoptAlgorithm::Reduced_SnoptAlgorithm");
168 template<
typename PROBLEM,
typename VECTOR>
176 template<
typename PROBLEM,
typename VECTOR>
182 throw DOpEException(
"To use this algorithm you need to have SNOPT installed! To use this set the WITH_SNOPT CompilerFlag.",
"Reduced_SnoptAlgorithm::Solve");
188 this->GetReducedProblem()->GetControlBoxConstraints(q_min,q_max);
190 ConstraintVector<VECTOR> constraints(this->GetReducedProblem()->GetProblem()->GetSpaceTimeHandler(),vector_behavior_);
194 double cost_start=0.;
195 std::stringstream out;
196 this->GetOutputHandler()->InitNewtonOut(out);
197 global_tol = std::max(opt_tol_,global_tol);
199 out <<
"**************************************************\n";
200 out <<
"* Starting Solution using SNOPT *\n";
201 out <<
"* Solving : "<<this->GetProblem()->GetName()<<
"\t*\n";
205 this->GetReducedProblem()->StateSizeInfo(out);
206 out <<
"* Constraints : ";
208 out <<
"**************************************************";
209 this->GetOutputHandler()->Write(out,1+this->GetBasePriority(),1,1);
211 this->GetOutputHandler()->SetIterationNumber(iter,
"Opt_Snopt"+postindex_);
213 this->GetOutputHandler()->Write(q,
"Control"+postindex_,
"control");
217 cost_start = cost = this->GetReducedProblem()->ComputeReducedCostFunctional(q);
221 this->GetExceptionHandler()->HandleCriticalException(e,
"Reduced_SnoptAlgorithm::Solve");
224 this->GetOutputHandler()->InitOut(out);
225 out<<
"CostFunctional: " << cost;
226 this->GetOutputHandler()->Write(out,2+this->GetBasePriority());
227 this->GetOutputHandler()->InitNewtonOut(out);
230 out<<
"************************************************\n";
231 out<<
"* Calling SNOPT *\n";
233 out<<
"* output will be written to logfile only! *\n";
235 out<<
"* output will not be written to logfile! *\n";
236 out<<
"************************************************\n\n";
237 this->GetOutputHandler()->Write(out,1+this->GetBasePriority());
239 this->GetOutputHandler()->DisallowAllOutput();
241 this->GetOutputHandler()->StartSaveCTypeOutputToLog();
245 DOpEWrapper::SNOPT_Problem RSAProb;
251 integer *iAfun = NULL;
252 integer *jAvar = NULL;
253 doublereal *A = NULL;
255 integer lenG = n*(neF);
258 integer *iGfun =
new integer[lenG];
259 integer *jGvar =
new integer[lenG];
261 doublereal *x =
new doublereal[n];
262 doublereal *xlow =
new doublereal[n];
263 doublereal *xupp =
new doublereal[n];
264 doublereal *xmul =
new doublereal[n];
265 integer *xstate =
new integer[n];
267 doublereal *F =
new doublereal[neF];
268 doublereal *Flow =
new doublereal[neF];
269 doublereal *Fupp =
new doublereal[neF];
270 doublereal *Fmul =
new doublereal[neF];
271 integer *Fstate =
new integer[neF];
275 char *xnames =
new char[nxnames*8];
276 char *Fnames =
new char[nFnames*8];
279 doublereal ObjAdd = 0;
283 const VECTOR& gv_qmin = q_min.GetSpacialVector();
286 for(
unsigned int i=0; i < n; i++)
289 xlow[i] = gv_qmin(i);
290 xupp[i] = gv_qmax(i);
298 for(
unsigned int j = 1; j<neF; j++)
305 for(
unsigned int i=0; i < n; i++)
314 RSAProb.setProblemSize( n, neF );
315 RSAProb.setObjective ( ObjRow, ObjAdd );
316 RSAProb.setA ( lenA, iAfun, jAvar, A );
317 RSAProb.setG ( lenG, iGfun, jGvar );
318 RSAProb.setX ( x, xlow, xupp, xmul, xstate );
319 RSAProb.setF ( F, Flow, Fupp, Fmul, Fstate );
320 RSAProb.setXNames ( xnames, nxnames );
321 RSAProb.setFNames ( Fnames, nFnames );
322 RSAProb.setProbName (
"RSA" );
323 RSAProb.setNeA ( neA );
324 RSAProb.setNeG ( neG );
327 RSAProb.setUserFun ( DOpEWrapper::SNOPT_A_userfunc_ );
329 RSAProb.setIntParameter(
"Derivative option", 1 );
330 RSAProb.setRealParameter(
"Function precison", func_prec_);
331 RSAProb.setRealParameter(
"Major optimality tolerance", global_tol);
332 RSAProb.setRealParameter(
"Major feasibility tolerance", feas_tol_);
333 RSAProb.setRealParameter(
"Minor feasibility tolerance", feas_tol_);
334 RSAProb.setIntParameter(
"Minor iterations limit", max_inner_iter_);
335 RSAProb.setIntParameter(
"Major iterations limit", max_outer_iter_);
337 ret_val = RSAProb.GetReturnStatus();
340 for(
unsigned int i=0; i < gv_q.size(); i++)
349 delete []iGfun;
delete []jGvar;
351 delete []x;
delete []xlow;
delete []xupp;
352 delete []xmul;
delete []xstate;
354 delete []F;
delete []Flow;
delete []Fupp;
355 delete []Fmul;
delete []Fstate;
357 delete []xnames;
delete []Fnames;
362 this->GetOutputHandler()->StopSaveCTypeOutputToLog();
363 this->GetOutputHandler()->ResumeOutput();
364 out<<
"\n************************************************\n";
365 out<<
"* SNOPT Finished *\n";
366 out<<
"* with Exit Code: "<<std::setw(3)<<ret_val;
368 out<<
" (success) *\n";
370 out<<
" (unknown error) *\n";
371 out<<
"************************************************\n";
373 this->GetOutputHandler()->Write(out,1+this->GetBasePriority());
376 this->GetOutputHandler()->SetIterationNumber(iter,
"Opt_Snopt"+postindex_);
378 this->GetOutputHandler()->Write(q,
"Control"+postindex_,
"control");
381 cost = this->GetReducedProblem()->ComputeReducedCostFunctional(q);
385 this->GetExceptionHandler()->HandleCriticalException(e,
"Reduced_SnoptAlgorithm::Solve");
388 this->GetOutputHandler()->InitOut(out);
389 out<<
"CostFunctional: " << cost;
390 this->GetOutputHandler()->Write(out,2+this->GetBasePriority());
391 this->GetOutputHandler()->InitNewtonOut(out);
394 this->GetReducedProblem()->ComputeReducedFunctionals(q);
398 this->GetExceptionHandler()->HandleCriticalException(e,
"Reduced_SnoptAlgorithm::Solve");
401 out <<
"**************************************************\n";
402 out <<
"* Stopping Solution Using SNOPT *\n";
403 out <<
"* Relative reduction in cost functional:"<<std::scientific << std::setw(11) << this->GetOutputHandler()->ZeroTolerance((cost-cost_start)/fabs(0.5*(cost_start+cost)),1.0) <<
" *\n";
405 out <<
"* Final value: "<<cost<<
" *\n";
406 out <<
"**************************************************";
407 this->GetOutputHandler()->Write(out,1+this->GetBasePriority(),1,1);
409 #endif //Endof ifdef WITHSNOPT
415 template <
typename PROBLEM,
typename VECTOR>
421 ConstraintVector<VECTOR> constraints(this->GetReducedProblem()->GetProblem()->GetSpaceTimeHandler(),vector_behavior_);
422 ControlVector<VECTOR> tmp(this->GetReducedProblem()->GetProblem()->GetSpaceTimeHandler(),vector_behavior_);
423 VECTOR& ref_x = tmp.GetSpacialVector();
424 assert((
int) ref_x.size() == *(data.n));
425 for(
unsigned int i=0; i < ref_x.size(); i++)
426 ref_x(i) = (data.x)[i];
428 assert(*(data.needF) > 0);
431 (data.F)[0] = this->GetReducedProblem()->ComputeReducedCostFunctional(tmp);
433 catch(DOpEException& e)
436 this->GetExceptionHandler()->HandleException(e,
"Reduced_SnoptAlgorithm::rsa_func_");
444 this->GetReducedProblem()->ComputeReducedConstraints(tmp,constraints);
446 catch(DOpEException& e)
449 this->GetExceptionHandler()->HandleException(e,
"Reduced_SnoptAlgorithm::rsa_func_");
452 const dealii::Vector<double>& gc = constraints.GetGlobalConstraints();
453 assert(*(data.neF) == (
int) gc.size()+1);
454 for(
unsigned int i=0; i < gc.size(); i++)
456 (data.F)[i+1] = gc(i);
460 if(*(data.needG) > 0)
462 ControlVector<VECTOR> gradient(tmp);
463 ControlVector<VECTOR> gradient_transposed(tmp);
467 this->GetReducedProblem()->ComputeReducedGradient(tmp,gradient,gradient_transposed);
469 catch(DOpEException& e)
472 this->GetExceptionHandler()->HandleException(e,
"Reduced_SnoptAlgorithm::rsa_func_");
475 assert(*(data.neG) == *(data.n)**(data.neF));
476 const VECTOR& ref_g = gradient_transposed.GetSpacialVector();
477 for(
unsigned int i=0; i < *(data.n); i++)
479 (data.G)[i] = ref_g(i);
482 const dealii::Vector<double>& gc = constraints.GetGlobalConstraints();
483 for(
unsigned int j=0; j < gc.size(); j++)
487 this->GetReducedProblem()->ComputeReducedGradientOfGlobalConstraints(j,tmp,constraints,gradient,gradient_transposed);
489 catch(DOpEException& e)
492 this->GetExceptionHandler()->HandleException(e,
"Reduced_SnoptAlgorithm::rsa_func_");
495 const VECTOR& ref_g = gradient_transposed.GetSpacialVector();
496 for(
unsigned int i=0; i < *(data.n); i++)
498 (data.G)[*(data.n)*(j+1)+i] = ref_g(i);
506 #endif //Endof ifdef WITHSNOPT
void PrintInfos(std::stringstream &out)
Definition: controlvector.cc:896
Definition: constraintvector.h:48
double get_double(const std::string &entry_name)
Definition: parameterreader.h:115
Definition: dopetypes.h:106
ControlType
Definition: dopetypes.h:103
bool get_bool(const std::string &entry_name)
Definition: parameterreader.h:148
Definition: parameterreader.h:36
Definition: optproblemcontainer.h:70
~Reduced_SnoptAlgorithm()
Definition: reduced_snopt_algorithm.h:169
const dealii::Vector< double > & GetGlobalConstraints() const
Definition: constraintvector.cc:174
void PrintInfos(std::stringstream &out)
Definition: constraintvector.cc:450
static void declare_params(ParameterReader ¶m_reader)
Definition: reduced_snopt_algorithm.h:110
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation=std::string())
Definition: parameterreader.h:98
Definition: dopetypes.h:105
Definition: controlvector.h:49
PROBLEM * GetProblem()
Definition: reducedalgorithm.h:187
PROBLEM * GetProblem()
Definition: reducedprobleminterface.h:490
Definition: reducedalgorithm.h:54
void ReInit()
Definition: controlvector.cc:96
virtual int Solve(ControlVector< VECTOR > &q, double global_tol=-1.)
Definition: reduced_snopt_algorithm.h:178
VECTOR & GetSpacialVector()
******************************************************/
Definition: controlvector.cc:204
Definition: reduced_snopt_algorithm.h:54
static void declare_params(ParameterReader ¶m_reader)
Definition: reducedalgorithm.h:241
int get_integer(const std::string &entry_name)
Definition: parameterreader.h:126
VectorStorageType
Definition: dopetypes.h:120
void SetSubsection(const std::string subsection)
Definition: parameterreader.h:93
Definition: reducedprobleminterface.h:335
Definition: dopeexception.h:35
Reduced_SnoptAlgorithm(PROBLEM *OP, ReducedProblemInterface< PROBLEM, VECTOR > *S, DOpEtypes::VectorStorageType vector_behavior, ParameterReader ¶m_reader, DOpEExceptionHandler< VECTOR > *Except=NULL, DOpEOutputHandler< VECTOR > *Output=NULL, int base_priority=0)
Definition: reduced_snopt_algorithm.h:137
std::string DOpEtypesToString(const C &)
Definition: dopetypes.h:134
Definition: optproblemcontainer.h:72