mirror of
https://git.dev.opencascade.org/repos/occt.git
synced 2025-08-04 13:13:25 +03:00
0024608: Development of methods of global optimization of multivariable function
math_GlobOptMin - new global optimization minimization algorithm Extrema_GlobOptFuncCC, Extrema_ExtCC, Extrema_ExtCC2d - implementation of GlobOptMin algorithm to extrema curve / curve Extrema_CurveCache - deleted as obsolete code ChFi3d_Builder.cxx - fixed processing of extrema math_NewtonMinimum.cxx - fixed step to avoid incorrect behavior Test cases modification to meet new behavior.
This commit is contained in:
aml
apn
@@ -8,4 +8,6 @@ math_Vector.hxx
|
||||
math_Vector.cxx
|
||||
math_IntegerVector.hxx
|
||||
math_IntegerVector.cxx
|
||||
math_SingleTab.hxx
|
||||
math_SingleTab.hxx
|
||||
math_GlobOptMin.hxx
|
||||
math_GlobOptMin.cxx
|
||||
@@ -51,6 +51,7 @@ is
|
||||
deferred class MultipleVarFunctionWithHessian;
|
||||
deferred class FunctionSet;
|
||||
deferred class FunctionSetWithDerivatives;
|
||||
imported GlobOptMin;
|
||||
class IntegerRandom;
|
||||
class Gauss;
|
||||
class GaussLeastSquare;
|
||||
|
||||
386
src/math/math_GlobOptMin.cxx
Normal file
386
src/math/math_GlobOptMin.cxx
Normal file
@@ -0,0 +1,386 @@
|
||||
// Created on: 2014-01-20
|
||||
// Created by: Alexaner Malyshev
|
||||
// Copyright (c) 2014-2014 OPEN CASCADE SAS
|
||||
//
|
||||
// This file is part of Open CASCADE Technology software library.
|
||||
//
|
||||
// This library is free software; you can redistribute it and/or modify it under
|
||||
// the terms of the GNU Lesser General Public License version 2.1 as published
|
||||
// by the Free Software Foundation, with special exception defined in the file
|
||||
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
|
||||
// distribution for complete text of the license and disclaimer of any warranty.
|
||||
//
|
||||
// Alternatively, this file may be used under the terms of Open CASCADE
|
||||
// commercial license or contractual agreement
|
||||
|
||||
#include <math_GlobOptMin.hxx>
|
||||
|
||||
#include <math_BFGS.hxx>
|
||||
#include <math_Matrix.hxx>
|
||||
#include <math_MultipleVarFunctionWithGradient.hxx>
|
||||
#include <math_MultipleVarFunctionWithHessian.hxx>
|
||||
#include <math_NewtonMinimum.hxx>
|
||||
#include <math_Powell.hxx>
|
||||
#include <Precision.hxx>
|
||||
#include <Standard_Integer.hxx>
|
||||
#include <Standard_Real.hxx>
|
||||
|
||||
const Handle(Standard_Type)& STANDARD_TYPE(math_GlobOptMin)
|
||||
{
|
||||
static Handle(Standard_Type) _atype = new Standard_Type ("math_GlobOptMin", sizeof (math_GlobOptMin));
|
||||
return _atype;
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : math_GlobOptMin
|
||||
//purpose : Constructor
|
||||
//=======================================================================
|
||||
math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
|
||||
const math_Vector& theA,
|
||||
const math_Vector& theB,
|
||||
Standard_Real theC)
|
||||
: myN(theFunc->NbVariables()),
|
||||
myA(1, myN),
|
||||
myB(1, myN),
|
||||
myGlobA(1, myN),
|
||||
myGlobB(1, myN),
|
||||
myX(1, myN),
|
||||
myTmp(1, myN),
|
||||
myV(1, myN)
|
||||
{
|
||||
Standard_Integer i;
|
||||
|
||||
myFunc = theFunc;
|
||||
myC = theC;
|
||||
myZ = -1;
|
||||
mySolCount = 0;
|
||||
|
||||
for(i = 1; i <= myN; i++)
|
||||
{
|
||||
myGlobA(i) = theA(i);
|
||||
myGlobB(i) = theB(i);
|
||||
|
||||
myA(i) = theA(i);
|
||||
myB(i) = theB(i);
|
||||
}
|
||||
|
||||
myDone = Standard_False;
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : SetGlobalParams
|
||||
//purpose : Set params without memory allocation.
|
||||
//=======================================================================
|
||||
void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
|
||||
const math_Vector& theA,
|
||||
const math_Vector& theB,
|
||||
Standard_Real theC)
|
||||
{
|
||||
Standard_Integer i;
|
||||
|
||||
myFunc = theFunc;
|
||||
myC = theC;
|
||||
myZ = -1;
|
||||
mySolCount = 0;
|
||||
|
||||
for(i = 1; i <= myN; i++)
|
||||
{
|
||||
myGlobA(i) = theA(i);
|
||||
myGlobB(i) = theB(i);
|
||||
|
||||
myA(i) = theA(i);
|
||||
myB(i) = theB(i);
|
||||
}
|
||||
|
||||
myDone = Standard_False;
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : SetLocalParams
|
||||
//purpose : Set params without memory allocation.
|
||||
//=======================================================================
|
||||
void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
|
||||
const math_Vector& theLocalB)
|
||||
{
|
||||
Standard_Integer i;
|
||||
|
||||
myZ = -1;
|
||||
mySolCount = 0;
|
||||
|
||||
for(i = 1; i <= myN; i++)
|
||||
{
|
||||
myA(i) = theLocalA(i);
|
||||
myB(i) = theLocalB(i);
|
||||
}
|
||||
|
||||
myDone = Standard_False;
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : ~math_GlobOptMin
|
||||
//purpose :
|
||||
//=======================================================================
|
||||
math_GlobOptMin::~math_GlobOptMin()
|
||||
{
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : Perform
|
||||
//purpose : Compute Global extremum point
|
||||
//=======================================================================
|
||||
// In this algo indexes started from 1, not from 0.
|
||||
void math_GlobOptMin::Perform()
|
||||
{
|
||||
Standard_Integer i;
|
||||
|
||||
// Compute parameters range
|
||||
Standard_Real minLength = RealLast();
|
||||
Standard_Real maxLength = RealFirst();
|
||||
for(i = 1; i <= myN; i++)
|
||||
{
|
||||
Standard_Real currentLength = myB(i) - myA(i);
|
||||
if (currentLength < minLength)
|
||||
minLength = currentLength;
|
||||
if (currentLength > maxLength)
|
||||
maxLength = currentLength;
|
||||
}
|
||||
|
||||
Standard_Real myTol = 1e-2;
|
||||
|
||||
myE1 = minLength * myTol / myC;
|
||||
myE2 = 2.0 * myTol / myC * maxLength;
|
||||
myE3 = - maxLength * myTol / 4;
|
||||
|
||||
// Compure start point.
|
||||
math_Vector aPnt(1,2);
|
||||
for(i = 1; i <= myN; i++)
|
||||
{
|
||||
Standard_Real currCentral = (myA(i) + myB(i)) / 2.0;
|
||||
aPnt(i) = currCentral;
|
||||
myY.Append(currCentral);
|
||||
}
|
||||
|
||||
myFunc->Value(aPnt, myF);
|
||||
mySolCount++;
|
||||
|
||||
computeGlobalExtremum(myN);
|
||||
|
||||
myDone = Standard_True;
|
||||
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : computeLocalExtremum
|
||||
//purpose :
|
||||
//=======================================================================
|
||||
Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt,
|
||||
Standard_Real& theVal,
|
||||
math_Vector& theOutPnt)
|
||||
{
|
||||
Standard_Integer i;
|
||||
|
||||
//Newton method
|
||||
if (dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
|
||||
{
|
||||
math_MultipleVarFunctionWithHessian* myTmp =
|
||||
dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
|
||||
|
||||
math_NewtonMinimum newtonMinimum(*myTmp, thePnt);
|
||||
if (newtonMinimum.IsDone())
|
||||
{
|
||||
newtonMinimum.Location(theOutPnt);
|
||||
theVal = newtonMinimum.Minimum();
|
||||
}
|
||||
else return Standard_False;
|
||||
} else
|
||||
|
||||
// BFGS method used.
|
||||
if (dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
|
||||
{
|
||||
math_MultipleVarFunctionWithGradient* myTmp =
|
||||
dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
|
||||
math_BFGS bfgs(*myTmp, thePnt);
|
||||
if (bfgs.IsDone())
|
||||
{
|
||||
bfgs.Location(theOutPnt);
|
||||
theVal = bfgs.Minimum();
|
||||
}
|
||||
else return Standard_False;
|
||||
} else
|
||||
|
||||
// Powell method used.
|
||||
if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
|
||||
{
|
||||
math_Matrix m(1, myN, 1, myN, 0.0);
|
||||
for(i = 1; i <= myN; i++)
|
||||
m(1, 1) = 1.0;
|
||||
|
||||
math_Powell powell(*myFunc, thePnt, m, 1e-10);
|
||||
|
||||
if (powell.IsDone())
|
||||
{
|
||||
powell.Location(theOutPnt);
|
||||
theVal = powell.Minimum();
|
||||
}
|
||||
else return Standard_False;
|
||||
}
|
||||
|
||||
if (isInside(theOutPnt))
|
||||
return Standard_True;
|
||||
else
|
||||
return Standard_False;
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : ComputeGlobalExtremum
|
||||
//purpose :
|
||||
//=======================================================================
|
||||
void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
|
||||
{
|
||||
Standard_Integer i;
|
||||
Standard_Real d; // Functional in moved point.
|
||||
Standard_Real val = RealLast(); // Local extrema computed in moved point.
|
||||
Standard_Real stepBestValue = RealLast();
|
||||
math_Vector stepBestPoint(1,2);
|
||||
Standard_Boolean isInside = Standard_False;
|
||||
Standard_Real r;
|
||||
|
||||
for(myX(j) = myA(j) + myE1; myX(j) < myB(j) + myE1; myX(j) += myV(j))
|
||||
{
|
||||
if (myX(j) > myB(j))
|
||||
myX(j) = myB(j);
|
||||
|
||||
if (j == 1)
|
||||
{
|
||||
isInside = Standard_False;
|
||||
myFunc->Value(myX, d);
|
||||
r = (d - myF) * myZ;
|
||||
if(r > myE3)
|
||||
{
|
||||
isInside = computeLocalExtremum(myX, val, myTmp);
|
||||
}
|
||||
stepBestValue = (isInside && (val < d))? val : d;
|
||||
stepBestPoint = (isInside && (val < d))? myTmp : myX;
|
||||
|
||||
// Solutions are close to each other.
|
||||
if (Abs(stepBestValue - myF) < Precision::Confusion() * 0.01)
|
||||
{
|
||||
if (!isStored(stepBestPoint))
|
||||
{
|
||||
if ((stepBestValue - myF) * myZ > 0.0)
|
||||
myF = stepBestValue;
|
||||
for(i = 1; i <= myN; i++)
|
||||
myY.Append(stepBestPoint(i));
|
||||
mySolCount++;
|
||||
}
|
||||
}
|
||||
|
||||
// New best solution.
|
||||
if ((stepBestValue - myF) * myZ > Precision::Confusion() * 0.01)
|
||||
{
|
||||
mySolCount = 0;
|
||||
myF = stepBestValue;
|
||||
myY.Clear();
|
||||
for(i = 1; i <= myN; i++)
|
||||
myY.Append(stepBestPoint(i));
|
||||
mySolCount++;
|
||||
}
|
||||
|
||||
myV(1) = myE2 + Abs(myF - d) / myC;
|
||||
}
|
||||
else
|
||||
{
|
||||
myV(j) = RealLast() / 2.0;
|
||||
computeGlobalExtremum(j - 1);
|
||||
}
|
||||
if ((j < myN) && (myV(j + 1) > myV(j)))
|
||||
{
|
||||
if (myV(j) > (myB(j + 1) - myA(j + 1)) / 3.0) // Case of too big step.
|
||||
myV(j + 1) = (myB(j + 1) - myA(j + 1)) / 3.0;
|
||||
else
|
||||
myV(j + 1) = myV(j);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : IsInside
|
||||
//purpose :
|
||||
//=======================================================================
|
||||
Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt)
|
||||
{
|
||||
Standard_Integer i;
|
||||
|
||||
for(i = 1; i <= myN; i++)
|
||||
{
|
||||
if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i))
|
||||
return Standard_False;
|
||||
}
|
||||
|
||||
return Standard_True;
|
||||
}
|
||||
//=======================================================================
|
||||
//function : IsStored
|
||||
//purpose :
|
||||
//=======================================================================
|
||||
Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
|
||||
{
|
||||
Standard_Integer i,j;
|
||||
Standard_Boolean isSame = Standard_True;
|
||||
|
||||
for(i = 0; i < mySolCount; i++)
|
||||
{
|
||||
isSame = Standard_True;
|
||||
for(j = 1; j <= myN; j++)
|
||||
{
|
||||
if ((Abs(thePnt(j) - myY(i * myN + j))) > (myB(j) - myA(j)) * Precision::Confusion())
|
||||
{
|
||||
isSame = Standard_False;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (isSame == Standard_True)
|
||||
return Standard_True;
|
||||
|
||||
}
|
||||
return Standard_False;
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : NbExtrema
|
||||
//purpose :
|
||||
//=======================================================================
|
||||
Standard_Integer math_GlobOptMin::NbExtrema()
|
||||
{
|
||||
return mySolCount;
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : GetF
|
||||
//purpose :
|
||||
//=======================================================================
|
||||
Standard_Real math_GlobOptMin::GetF()
|
||||
{
|
||||
return myF;
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : IsDone
|
||||
//purpose :
|
||||
//=======================================================================
|
||||
Standard_Boolean math_GlobOptMin::isDone()
|
||||
{
|
||||
return myDone;
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
//function : Points
|
||||
//purpose :
|
||||
//=======================================================================
|
||||
void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
|
||||
{
|
||||
Standard_Integer j;
|
||||
|
||||
for(j = 1; j <= myN; j++)
|
||||
theSol(j) = myY((theIndex - 1) * myN + j);
|
||||
}
|
||||
101
src/math/math_GlobOptMin.hxx
Normal file
101
src/math/math_GlobOptMin.hxx
Normal file
@@ -0,0 +1,101 @@
|
||||
// Created on: 2014-01-20
|
||||
// Created by: Alexaner Malyshev
|
||||
// Copyright (c) 2014-2014 OPEN CASCADE SAS
|
||||
//
|
||||
// This file is part of Open CASCADE Technology software library.
|
||||
//
|
||||
// This library is free software; you can redistribute it and/or modify it under
|
||||
// the terms of the GNU Lesser General Public License version 2.1 as published
|
||||
// by the Free Software Foundation, with special exception defined in the file
|
||||
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
|
||||
// distribution for complete text of the license and disclaimer of any warranty.
|
||||
//
|
||||
// Alternatively, this file may be used under the terms of Open CASCADE
|
||||
// commercial license or contractual agreement.
|
||||
|
||||
#ifndef _math_GlobOptMin_HeaderFile
|
||||
#define _math_GlobOptMin_HeaderFile
|
||||
|
||||
#include <math_MultipleVarFunction.hxx>
|
||||
#include <NCollection_Sequence.hxx>
|
||||
#include <Standard_Type.hxx>
|
||||
|
||||
//! This class represents Evtushenko's algorithm of global optimization based on nonuniform mesh.<br>
|
||||
//! Article: Yu. Evtushenko. Numerical methods for finding global extreme (case of a non-uniform mesh). <br>
|
||||
//! U.S.S.R. Comput. Maths. Math. Phys., Vol. 11, N 6, pp. 38-54.
|
||||
|
||||
class math_GlobOptMin
|
||||
{
|
||||
public:
|
||||
|
||||
Standard_EXPORT math_GlobOptMin(math_MultipleVarFunction* theFunc,
|
||||
const math_Vector& theA,
|
||||
const math_Vector& theB,
|
||||
Standard_Real theC = 9);
|
||||
|
||||
Standard_EXPORT void SetGlobalParams(math_MultipleVarFunction* theFunc,
|
||||
const math_Vector& theA,
|
||||
const math_Vector& theB,
|
||||
Standard_Real theC = 9);
|
||||
|
||||
Standard_EXPORT void SetLocalParams(const math_Vector& theLocalA,
|
||||
const math_Vector& theLocalB);
|
||||
|
||||
Standard_EXPORT ~math_GlobOptMin();
|
||||
|
||||
Standard_EXPORT void Perform();
|
||||
|
||||
//! Get best functional value.
|
||||
Standard_EXPORT Standard_Real GetF();
|
||||
|
||||
//! Return count of global extremas. NbExtrema <= MAX_SOLUTIONS.
|
||||
Standard_EXPORT Standard_Integer NbExtrema();
|
||||
|
||||
//! Return solution i, 1 <= i <= NbExtrema.
|
||||
Standard_EXPORT void Points(const Standard_Integer theIndex, math_Vector& theSol);
|
||||
|
||||
private:
|
||||
|
||||
math_GlobOptMin & operator = (const math_GlobOptMin & theOther);
|
||||
|
||||
Standard_Boolean computeLocalExtremum(const math_Vector& thePnt, Standard_Real& theVal, math_Vector& theOutPnt);
|
||||
|
||||
void computeGlobalExtremum(Standard_Integer theIndex);
|
||||
|
||||
//! Check that myA <= pnt <= myB
|
||||
Standard_Boolean isInside(const math_Vector& thePnt);
|
||||
|
||||
Standard_Boolean isStored(const math_Vector &thePnt);
|
||||
|
||||
Standard_Boolean isDone();
|
||||
|
||||
// Input.
|
||||
math_MultipleVarFunction* myFunc;
|
||||
Standard_Integer myN;
|
||||
math_Vector myA; // Left border on current C2 interval.
|
||||
math_Vector myB; // Right border on current C2 interval.
|
||||
math_Vector myGlobA; // Global left border.
|
||||
math_Vector myGlobB; // Global right border.
|
||||
|
||||
// Output.
|
||||
Standard_Boolean myDone;
|
||||
NCollection_Sequence<Standard_Real> myY;// Current solutions.
|
||||
Standard_Integer mySolCount; // Count of solutions.
|
||||
|
||||
// Algorithm data.
|
||||
Standard_Real myZ;
|
||||
Standard_Real myC; //Lipschitz constant
|
||||
Standard_Real myE1; // Border coeff.
|
||||
Standard_Real myE2; // Minimum step size.
|
||||
Standard_Real myE3; // Local extrema starting parameter.
|
||||
|
||||
math_Vector myX; // Current modified solution
|
||||
math_Vector myTmp; // Current modified solution
|
||||
math_Vector myV; // Steps array.
|
||||
|
||||
Standard_Real myF; // Current value of Global optimum.
|
||||
};
|
||||
|
||||
const Handle(Standard_Type)& TYPE(math_GlobOptMin);
|
||||
|
||||
#endif
|
||||
@@ -143,12 +143,19 @@ void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F,
|
||||
}
|
||||
|
||||
LU.Solve(TheGradient, TheStep);
|
||||
*suivant = *precedent - TheStep;
|
||||
Standard_Boolean hasProblem = Standard_False;
|
||||
do
|
||||
{
|
||||
*suivant = *precedent - TheStep;
|
||||
|
||||
// Gestion de la convergence
|
||||
hasProblem = !(F.Value(*suivant, TheMinimum));
|
||||
|
||||
// Gestion de la convergence
|
||||
|
||||
F.Value(*suivant, TheMinimum);
|
||||
if (hasProblem)
|
||||
{
|
||||
TheStep /= 2.0;
|
||||
}
|
||||
} while (hasProblem);
|
||||
|
||||
if (IsConverged()) { NbConv++; }
|
||||
else { NbConv=0; }
|
||||
|
||||
Reference in New Issue
Block a user