diff --git a/src/Trim.cpp b/src/Trim.cpp index 8022c8a83c3af81f412d3f77134b2c1f066bd197..73065a7c9583fa19ed48ea7c6aa2930c362801f2 100644 --- a/src/Trim.cpp +++ b/src/Trim.cpp @@ -25,6 +25,7 @@ #include "models/propulsion/FGTurbine.h" #include "models/propulsion/FGTurboProp.h" #include "math/FGNelderMead.h" +#include <stdexcept> template <class varType> void prompt(const std::string & str, varType & var) @@ -43,7 +44,9 @@ class Callback : public JSBSim::FGNelderMead::Callback public: void eval(const std::vector<double> &v) { - for (int i=0;i<v.size();i++) std::cout << "v[" << i << "] = " << v[i] << std::endl; + std::cout << "v: "; + for (int i=0;i<v.size();i++) std::cout << v[i] << " "; + std::cout << std::endl; } } callback; @@ -185,12 +188,31 @@ int main (int argc, char const* argv[]) initialGuess[5] = 0; // beta // solve - FGTrimmer trimmer(fdm, constraints); - FGNelderMead solver(trimmer,initialGuess, lowerBound, upperBound, initialStepSize, - iterMax,rtol,abstol,speed,showConvergeStatus,showSimplex,pause,&callback); + FGTrimmer trimmer(fdm, constraints); + FGNelderMead * solver; + try + { + solver = new FGNelderMead(trimmer,initialGuess, + lowerBound, upperBound, initialStepSize,iterMax,rtol, + abstol,speed,showConvergeStatus,showSimplex,pause,&callback); + } + catch (const std::runtime_error & e) + { + std::cout << e.what() << std::endl; + return 1; + } // output - trimmer.printSolution(solver.getSolution()); // this also loads the solution into the fdm + try + { + trimmer.printSolution(solver->getSolution()); // this also loads the solution into the fdm + } + catch(std::runtime_error & e) + { + std::cout << "caught std::runtime error" << std::endl; + std::cout << "exception: " << e.what() << std::endl; + return 1; + } //std::cout << "\nsimulating flight to determine trim stability" << std::endl; diff --git a/src/initialization/FGTrimmer.cpp b/src/initialization/FGTrimmer.cpp index e77932f963f8a5904b97b7da99431bdae44a7685..ed46681379739043b1294c5c1c13e8a28cfd746b 100644 --- a/src/initialization/FGTrimmer.cpp +++ b/src/initialization/FGTrimmer.cpp @@ -25,6 +25,7 @@ #include "models/FGAircraft.h" #include <iomanip> #include <cstdlib> +#include <stdexcept> namespace JSBSim { @@ -365,13 +366,13 @@ double FGTrimmer::eval(const std::vector<double> & v) } else if (iter>100) { - std::cout << "\ncost failed to converge to steady value" - << std::scientific - << "\ndelta dvt: " << std::abs(dvt0-dvt) - << "\nmost likely out of the flight envelope" - << "\ncheck constraints and initial conditions" - << std::endl; - exit(1); + //std::cout << "\ncost failed to converge to steady value" + //<< std::scientific + //<< "\ndelta dvt: " << std::abs(dvt0-dvt) + //<< "\nmost likely out of the flight envelope" + //<< "\ncheck constraints and initial conditions" + //<< std::endl; + throw(std::runtime_error("FGTrimmer: cost failed to converge to steady value, most likely out of flight envelope, check constraints and initial conditions")); } else dvt0=dvt; } diff --git a/src/math/FGNelderMead.cpp b/src/math/FGNelderMead.cpp index 04ad985bf9e23cb114d87ee77679f9153df974eb..9917ee89b285b9e0699ccc835c7d0420e47f4b83 100644 --- a/src/math/FGNelderMead.cpp +++ b/src/math/FGNelderMead.cpp @@ -22,6 +22,7 @@ #include <cstdlib> #include <iomanip> #include <iostream> +#include <stdexcept> namespace JSBSim { @@ -36,7 +37,7 @@ FGNelderMead::FGNelderMead(Function & f, const std::vector<double> & initialGues m_nDim(initialGuess.size()), m_nVert(m_nDim+1), m_iMax(1), m_iNextMax(1), m_iMin(1), m_simplex(m_nVert), m_cost(m_nVert), m_elemSum(m_nDim), - m_showSimplex(showSimplex), m_callback(callback) + m_showSimplex(showSimplex), m_callback(callback), m_stopRequested(false) { // setup std::cout.precision(3); @@ -46,7 +47,7 @@ FGNelderMead::FGNelderMead(Function & f, const std::vector<double> & initialGues int iter = 0; // solve simplex - while (1) + while (!m_stopRequested) { // reinitialize simplex whenever rtol condition is met if ( rtolI < rtol || iter == 0) @@ -74,7 +75,15 @@ FGNelderMead::FGNelderMead(Function & f, const std::vector<double> & initialGues // find vertex costs for (int vertex=0;vertex<m_nVert;vertex++) { - m_cost[vertex] = m_f.eval(m_simplex[vertex]); + try + { + m_cost[vertex] = m_f.eval(m_simplex[vertex]); + } + catch (const std::exception & e) + { + throw; + return; + } if (m_callback) m_callback->eval(m_simplex[vertex]); } @@ -179,29 +188,38 @@ FGNelderMead::FGNelderMead(Function & f, const std::vector<double> & initialGues // costs - - // try inversion - double costTry = tryStretch(-1.0); - - // if lower cost than best, then try further stretch by speed factor - if (costTry < minCost) - { - costTry = tryStretch(speed); - } - // otherwise try a contraction - else if (costTry > nextMaxCost) - { - // 1d contraction - costTry = tryStretch(1./speed); - - // if greater than max cost, contract about min - if (costTry > maxCost) - { - if (showSimplex) - std::cout << "multiD contraction about: " << m_iMin << std::endl; - contract(); - } - } + + try + { + // try inversion + double costTry = tryStretch(-1.0); + + // if lower cost than best, then try further stretch by speed factor + if (costTry < minCost) + { + costTry = tryStretch(speed); + } + // otherwise try a contraction + else if (costTry > nextMaxCost) + { + // 1d contraction + costTry = tryStretch(1./speed); + + // if greater than max cost, contract about min + if (costTry > maxCost) + { + if (showSimplex) + std::cout << "multiD contraction about: " << m_iMin << std::endl; + contract(); + } + } + } + + catch (const std::exception & e) + { + throw; + stop(); + } // iteration iter++; @@ -213,6 +231,12 @@ FGNelderMead::FGNelderMead(Function & f, const std::vector<double> & initialGues std::cout << std::fixed; } +void FGNelderMead::stop() +{ + std::cout << "FGNelderMead: stop requested." << std::endl; + m_stopRequested = true; +} + std::vector<double> FGNelderMead::getSolution() { return m_simplex[m_iMin]; @@ -231,8 +255,17 @@ double FGNelderMead::tryStretch(double factor) } // find trial cost - double costTry0 = m_f.eval(tryVertex); - double costTry = m_f.eval(tryVertex); + double costTry0 = 0, costTry = 0; + try + { + costTry0 = m_f.eval(tryVertex); + costTry = m_f.eval(tryVertex); + } + catch(const std::exception & e) + { + throw; + return 0; + } if (std::abs(costTry0-costTry) > std::numeric_limits<float>::epsilon()) { diff --git a/src/math/FGNelderMead.h b/src/math/FGNelderMead.h index 822a0573ae3fe6cef9226fddb5442084db3fa186..56395aea4925f7d5f72dfcdc2a057e5cd0404713 100644 --- a/src/math/FGNelderMead.h +++ b/src/math/FGNelderMead.h @@ -52,6 +52,9 @@ public: bool pause=false, Callback * callback=NULL); std::vector<double> getSolution(); + + void stop(); + private: // attributes Function & m_f; @@ -63,6 +66,7 @@ private: std::vector<double> m_cost; std::vector<double> m_elemSum; bool m_showSimplex; + bool m_stopRequested; // methods double tryStretch(double factor);