diff --git a/src/Trim.cpp b/src/Trim.cpp
index 8326f119f6da9a0a99a3509b19c7558522c40466..a62c0027cc1ba889a4073e6ecaa31b317ce717bf 100644
--- a/src/Trim.cpp
+++ b/src/Trim.cpp
@@ -69,6 +69,7 @@ int main (int argc, char const* argv[])
double rtol = std::numeric_limits<float>::epsilon();
double abstol = std::numeric_limits<double>::epsilon();
double speed = 2.0;
+ double random = 0.1; // 10% random factor added to all simplex calcs
int iterMax = 2000;
bool showConvergeStatus = false;
bool pause = false;
@@ -194,7 +195,7 @@ int main (int argc, char const* argv[])
{
solver = new FGNelderMead(trimmer,initialGuess,
lowerBound, upperBound, initialStepSize,iterMax,rtol,
- abstol,speed,showConvergeStatus,showSimplex,pause,&callback);
+ abstol,speed,random,showConvergeStatus,showSimplex,pause,&callback);
while(solver->status()==1) solver->update();
}
catch (const std::runtime_error & e)
diff --git a/src/math/FGNelderMead.cpp b/src/math/FGNelderMead.cpp
index f62a3a8831d11d039effa29f80126ad4e8f93303..2d93944883aa9e4c63eed7dfbfca49c515b804fd 100644
--- a/src/math/FGNelderMead.cpp
+++ b/src/math/FGNelderMead.cpp
@@ -23,6 +23,7 @@
#include <iomanip>
#include <iostream>
#include <stdexcept>
+#include <ctime>
namespace JSBSim
{
@@ -31,19 +32,22 @@ FGNelderMead::FGNelderMead(Function & f, const std::vector<double> & initialGues
const std::vector<double> & lowerBound,
const std::vector<double> & upperBound,
const std::vector<double> & initialStepSize, int iterMax,
- double rtol, double abstol, double speed, bool showConvergeStatus,
+ double rtol, double abstol, double speed, double randomization,
+ bool showConvergeStatus,
bool showSimplex, bool pause, Callback * callback) :
m_f(f), m_lowerBound(lowerBound), m_upperBound(upperBound),
+ m_randomization(randomization),
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_status(1),
+ m_callback(callback), m_status(1),
initialGuess(initialGuess), initialStepSize(initialStepSize),
iterMax(iterMax), rtol(rtol), abstol(abstol),
speed(speed), showConvergeStatus(showConvergeStatus), showSimplex(showSimplex),
pause(pause), rtolI(), minCostPrevResize(1), minCost(), minCostPrev(), maxCost(),
nextMaxCost(), iter()
{
+ srand ( time(NULL) ); // seed random number generator
}
void FGNelderMead::update()
@@ -197,17 +201,27 @@ void FGNelderMead::update()
{
// try inversion
double costTry = tryStretch(-1.0);
+ //std::cout << "cost Try 0: " << costTry << std::endl;
- // if lower cost than best, then try further stretch by speed factor
+ // if lower cost than best, then try further stretch by double speed factor
if (costTry < minCost)
{
+ double costTry0 = costTry;
costTry = tryStretch(speed);
+ //std::cout << "cost Try 1: " << costTry << std::endl;
+
+ if (showSimplex)
+ {
+ if (costTry < costTry0) std::cout << "inversion about: " << m_iMax << std::endl;
+ else std::cout << "inversion and stretch about: " << m_iMax << std::endl;
+ }
}
// otherwise try a contraction
else if (costTry > nextMaxCost)
{
// 1d contraction
costTry = tryStretch(1./speed);
+ //std::cout << "cost Try 2: " << costTry << std::endl;
// if greater than max cost, contract about min
if (costTry > maxCost)
@@ -216,6 +230,11 @@ void FGNelderMead::update()
std::cout << "multiD contraction about: " << m_iMin << std::endl;
contract();
}
+ else
+ {
+ if (showSimplex)
+ std::cout << "contraction about: " << m_iMin << std::endl;
+ }
}
}
@@ -244,6 +263,13 @@ int FGNelderMead::status()
return m_status;
}
+double FGNelderMead::getRandomFactor()
+{
+ double randFact = 1+(float(rand() % 1000)/500-1)*m_randomization;
+ //std::cout << "random factor: " << randFact << std::endl;;
+ return randFact;
+}
+
std::vector<double> FGNelderMead::getSolution()
{
return m_simplex[m_iMin];
@@ -251,6 +277,9 @@ std::vector<double> FGNelderMead::getSolution()
double FGNelderMead::tryStretch(double factor)
{
+ // randomize factor so we can avoid locking situations
+ factor = factor*getRandomFactor();
+
// create trial vertex
double a= (1.0-factor)/m_nDim;
double b = a - factor;
@@ -276,8 +305,8 @@ double FGNelderMead::tryStretch(double factor)
if (std::abs(costTry0-costTry) > std::numeric_limits<float>::epsilon())
{
- std::cout << "\twarning: dynamics not stable!" << std::endl;
- return 2*m_cost[m_iMax];
+ //std::cout << "\twarning: dynamics not stable!" << std::endl;
+ //return 1000000*m_cost[m_iMax];
}
// if trial cost lower than max
@@ -290,7 +319,7 @@ double FGNelderMead::tryStretch(double factor)
for (int dim=0;dim<m_nDim;dim++) m_simplex[m_iMax][dim] = tryVertex[dim];
// update the cost
m_cost[m_iMax] = costTry;
- if (m_showSimplex) std::cout << "stretched\t" << m_iMax << "\tby : " << factor << std::endl;
+ if (showSimplex) std::cout << "stretched\t" << m_iMax << "\tby : " << factor << std::endl;
}
return costTry;
}
@@ -302,7 +331,7 @@ void FGNelderMead::contract()
for (int vertex=0;vertex<m_nVert;vertex++)
{
m_simplex[vertex][dim] =
- 0.5*(m_simplex[vertex][dim] +
+ getRandomFactor()*0.5*(m_simplex[vertex][dim] +
m_simplex[m_iMin][dim]);
}
}
@@ -314,16 +343,15 @@ void FGNelderMead::constructSimplex(const std::vector<double> & guess,
for (int vertex=0;vertex<m_nVert;vertex++)
{
m_simplex[vertex] = guess;
- std::vector<double> upperBound(guess.size());
- for (int dim=0;dim<m_nDim;dim++) upperBound[dim] = m_upperBound[dim]-stepSize[dim];
- boundVertex(m_simplex[vertex],m_lowerBound,upperBound);
- }
+ }
+
for (int dim=0;dim<m_nDim;dim++)
{
int vertex = dim + 1;
- m_simplex[vertex][dim] += stepSize[dim];
+ m_simplex[vertex][dim] += stepSize[dim]*getRandomFactor();
+ boundVertex(m_simplex[vertex],m_lowerBound,m_upperBound);
}
- if (m_showSimplex)
+ if (showSimplex)
{
std::cout << "simplex: " << std::endl;;
for (int j=0;j<m_nVert;j++) std::cout << "\t\t" << j;
diff --git a/src/math/FGNelderMead.h b/src/math/FGNelderMead.h
index 5ded8614da109e34df6507d9f6537981514303dc..e21053c108019607c9b6f63cb7294a07762d4acb 100644
--- a/src/math/FGNelderMead.h
+++ b/src/math/FGNelderMead.h
@@ -48,6 +48,7 @@ public:
double rtol=std::numeric_limits<float>::epsilon(),
double abstol=std::numeric_limits<float>::epsilon(),
double speed = 2.0,
+ double randomization=0.1,
bool showConvergeStatus=true,bool showSimplex=false,
bool pause=false,
Callback * callback=NULL);
@@ -60,13 +61,13 @@ private:
// attributes
Function & m_f;
Callback * m_callback;
+ double m_randomization;
const std::vector<double> & m_lowerBound;
const std::vector<double> & m_upperBound;
int m_nDim, m_nVert, m_iMax, m_iNextMax, m_iMin;
std::vector< std::vector<double> > m_simplex;
std::vector<double> m_cost;
std::vector<double> m_elemSum;
- bool m_showSimplex;
int m_status;
const std::vector<double> & initialGuess;
const std::vector<double> & initialStepSize;
@@ -77,6 +78,7 @@ private:
maxCost, nextMaxCost;
// methods
+ double getRandomFactor();
double tryStretch(double factor);
void contract();
void constructSimplex(const std::vector<double> & guess, const std::vector<double> & stepSize);