Apply Newton Method to Find Extrema in OPEN CASCADE

Apply Newton Method to Find Extrema in OPEN CASCADE

eryar@163.com

Abstract. In calculus, Newton’s method is used for finding the roots of a function. In optimization, Newton’s method is applied to find the roots of the derivative. OPEN CASCADE implement Newton method to find the extrema for a multiple variables function, such as find the extrema point for a curve and a surface.

Key Words. Nonlinear Programming, Newton Method, Extrema, OPEN CASCADE

1. Introduction

Newton法作为一种经典的解无约束优化问题的方法,在20世纪80~90年代发展起来的解线性规划和凸规划的内点法中起到了重要作用。Newton法最初是Newton提出用于解非线性方程的,Newton曾用该法求解Kepler方程x-asinx=b,并得到精度很高的近似解。通过《OPEN CASCADE Multiple Variable Function》对OPEN CASCADE中多元函数的表达有了一个认识。多元函数如何应用的呢?下面提出一个问题及如何用程序来解决这个问题。对于任意给定的曲线和曲面,如何求出曲线和曲面上距离最近的点,假设曲线和曲面都是至少C2连续的。关于参数连续性可参考《OPEN CASCADE Curve Continuity》。如下图所示:

Figure 1.1 A Curve and A Surface

本文给出OPEN CASCADE中对此类问题的一种解法,即应用Newton法求解非线性无约束多元函数的极值。学习如何将实际问题抽象成数学模型,从而使用数学的方法进行求解。

2.Construct Function

在实际应用中,一个问题是不是可以表述为一个最优化模型和怎么表示为一个最优化模型,这是优化方法是否可以应用的前提,因而十分重要。但优化问题的建模和其他数学问题的建模一样,不属于精确科学或数学的范畴,而是一项技术或技艺,没有统一标准和方法。当然建立的模型是否正确和模型的优劣是可以通过实际效果来检验的。

OPEN CASCADE使用从math_MultipleVarFunctionWithHessian派生的一个具体类Extrema_GlobOptFuncCS来计算C2连续的曲线和曲面之间的距离的平方值。抽象出来的数学模型为:

因为是从具有Hessian Matrix的多元函数派生,所以要求曲线曲面具有至少C2连续,即有至少有二阶导数。且在类中分别实现计算函数值,计算一阶导数值(梯度),计算二阶导数值(Hessian Matrix)。计算函数值的代码如下所示:

//======================================================================
//function : value
//purpose  : 
//======================================================================
void Extrema_GlobOptFuncCS::value(Standard_Real cu,
                                  Standard_Real su,
                                  Standard_Real sv,
                                  Standard_Real &F)
{
  F = myC->Value(cu).SquareDistance(myS->Value(su, sv));
}

其中参数cu为参数曲线的参数,su,sv分别为参数曲面的参数。根据参数曲线曲面上的参数计算出相应的点,然后计算出两个点之间的距离的平方值即为函数值。与上述公式对应。

根据多元函数一阶导数(梯度)的定义,可得出梯度的计算公式如下:

计算梯度的代码如下所示,从程序代码可见程序就是上公式的具体实现。

 

//======================================================================
//function : gradient
//purpose  : 
//======================================================================
void Extrema_GlobOptFuncCS::gradient(Standard_Real cu,
                                     Standard_Real su,
                                     Standard_Real sv,
                                     math_Vector &G)
{
  gp_Pnt CD0, SD0;
  gp_Vec CD1, SD1U, SD1V;

  myC->D1(cu, CD0, CD1);
  myS->D1(su, sv, SD0, SD1U, SD1V);

  G(1) = + (CD0.X() - SD0.X()) * CD1.X()
         + (CD0.Y() - SD0.Y()) * CD1.Y()
         + (CD0.Z() - SD0.Z()) * CD1.Z();
  G(2) = - (CD0.X() - SD0.X()) * SD1U.X()
         - (CD0.Y() - SD0.Y()) * SD1U.Y()
         - (CD0.Z() - SD0.Z()) * SD1U.Z();
  G(3) = - (CD0.X() - SD0.X()) * SD1V.X()
         - (CD0.Y() - SD0.Y()) * SD1V.Y()
         - (CD0.Z() - SD0.Z()) * SD1V.Z();
}

根据Hessian Matrix的定义,得到计算Hessian Matrix的公式如下:

将函数积的求导法则应用于求偏导数得到上述公式。同理求出Hessian Matrix的其他各项,如下公式所示:

计算多元函数的二阶导数Hessian Matrix的程序代码如下所示:

 

//======================================================================
//function : hessian
//purpose  : 
//======================================================================
void Extrema_GlobOptFuncCS::hessian(Standard_Real cu,
                                    Standard_Real su,
                                    Standard_Real sv,
                                    math_Matrix &H)
{
  gp_Pnt CD0, SD0;
  gp_Vec CD1, SD1U, SD1V, CD2, SD2UU, SD2UV, SD2VV;

  myC->D2(cu, CD0, CD1, CD2);
  myS->D2(su, sv, SD0, SD1U, SD1V, SD2UU, SD2VV, SD2UV);

  H(1,1) = + CD1.X() * CD1.X()
           + CD1.Y() * CD1.Y()
           + CD1.Z() * CD1.Z()
           + (CD0.X() - SD0.X()) * CD2.X()
           + (CD0.Y() - SD0.Y()) * CD2.Y()
           + (CD0.Z() - SD0.Z()) * CD2.Z();

  H(1,2) = - CD1.X() * SD1U.X()
           - CD1.Y() * SD1U.Y()
           - CD1.Z() * SD1U.Z();

  H(1,3) = - CD1.X() * SD1V.X()
           - CD1.Y() * SD1V.Y()
           - CD1.Z() * SD1V.Z();

  H(2,1) = H(1,2);

  H(2,2) = + SD1U.X() * SD1U.X()
           + SD1U.Y() * SD1U.Y()
           + SD1U.Z() * SD1U.Z()
           - (CD0.X() - SD0.X()) * SD2UU.X()
           - (CD0.Y() - SD0.Y()) * SD2UU.Y()
           - (CD0.Z() - SD0.Z()) * SD2UU.Z();

  H(2,3) = + SD1U.X() * SD1V.X()
           + SD1U.Y() * SD1V.Y()
           + SD1U.Z() * SD1V.Z()
           - (CD0.X() - SD0.X()) * SD2UV.X()
           - (CD0.Y() - SD0.Y()) * SD2UV.Y()
           - (CD0.Z() - SD0.Z()) * SD2UV.Z();

  H(3,1) = H(1,3);

  H(3,2) = H(2,3);

  H(3,3) = + SD1V.X() * SD1V.X()
           + SD1V.Y() * SD1V.Y()
           + SD1V.Z() * SD1V.Z()
           - (CD0.X() - SD0.X()) * SD2VV.X()
           - (CD0.Y() - SD0.Y()) * SD2VV.Y()
           - (CD0.Z() - SD0.Z()) * SD2VV.Z();
}

根据高阶偏导数的定理可知,当f(X)在点X0处所有二阶偏导数连续时,那末在该区域内这两个二阶混合偏导数必相等。所以Hessian Matrix为一个对称矩阵,故

H(2,1)=H(1,2)

H(3,1)=H(1,3)

H(3,2)=H(2,3)

由此完成一个具有二阶偏导数的多元函数的数学模型,用面向对象的方式清晰地表示了出来。和我见过的国内一些程序相比,这种抽象思路还是很清晰,便于程序的理解和扩展。国内有个图形库是C风格的,一个函数几千行,光函数参数就很多,参数名也很随意,函数内部变量名称更是无法理解,什么i,j,k,ii,jj之类。这种程序别说可扩展,就是维护起来也是让人头疼的啊!

有了函数表达式,下面就是计算这个函数的极值了。

3.Newton’s Method

关于应用Newton法计算一元非线性方程的根已经在《OpenCASCADE Root-Finding Algorithm》中进行了说明,这里要学习下如何使用Newton法应用于多元函数极值的计算。对于一元函数f(x)的求极值问题,当f(x)连续可微时,最优点x满足f’(x)=0。于是当f(x)二次连续可微时,求解f’(x)=0的Newton法为:

该方法称为解无约束优化问题的Newton方法。由《数学分析》可知,当f(x)是凸函数时,f’(x)=0的解是minf(x)的整体最优解。将Newton法扩展到多元函数的情况,也是利用二阶Taylor级数将函数展开,得到多元函数的极值迭代公式:

关于多元函数Newton法公式的推导,可参考《最优化方法》等书籍。Newton法的算法步骤如下:

A. 给定初始点,及精度;

B. 计算函数f(xk)的一阶导数(梯度),二阶导数(Hessian Matrix):若|梯度|<精度,则停止迭代,输出近似极小点;否则转C;

C. 根据Newton迭代公式,计算x(k+1);

OPEN CASCADE中Newton法计算极值的类是math_NewtonMinimum,可参考其代码学习。下面给出前面提出的曲线曲面极值求解的实现代码:

 

/*
*    Copyright (c) 2015 Shing Liu All Rights Reserved.
*
*           File : main.cpp
*         Author : Shing Liu(eryar@163.com)
*           Date : 2015-12-05 21:00
*        Version : OpenCASCADE6.9.0
*
*    Description : Learn Newton's Method for multiple variables
*                  function.
*/

#define WNT
#include <TColgp_Array1OfPnt.hxx>
#include <TColgp_Array2OfPnt.hxx>

#include <math_NewtonMinimum.hxx>

#include <GeomTools.hxx>
#include <BRepTools.hxx>

#include <GC_MakeSegment.hxx>

#include <GeomAdaptor_HCurve.hxx>
#include <GeomAdaptor_Surface.hxx>

#include <Extrema_GlobOptFuncCS.hxx>

#include <GeomAPI_PointsToBSpline.hxx>
#include <GeomAPI_PointsToBSplineSurface.hxx>

#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>

#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")
#pragma comment(lib, "TKG2d.lib")
#pragma comment(lib, "TKG3d.lib")
#pragma comment(lib, "TKBRep.lib")
#pragma comment(lib, "TKGeomBase.lib")
#pragma comment(lib, "TKGeomAlgo.lib")
#pragma comment(lib, "TKTopAlgo.lib")

void testNewtonMethod(void)
{
    // approximate curve from the points
    TColgp_Array1OfPnt aCurvePoints(1, 5);
    aCurvePoints.SetValue(1, gp_Pnt(0.0, 0.0, -2.0));
    aCurvePoints.SetValue(2, gp_Pnt(1.0, 2.0, 2.0));
    aCurvePoints.SetValue(3, gp_Pnt(2.0, 3.0, 3.0));
    aCurvePoints.SetValue(4, gp_Pnt(4.0, 3.0, 4.0));
    aCurvePoints.SetValue(5, gp_Pnt(5.0, 5.0, 5.0));

    GeomAPI_PointsToBSpline aCurveApprox(aCurvePoints);

    // approximate surface from the points.
    TColgp_Array2OfPnt aSurfacePoints(1, 5, 1, 5);
    aSurfacePoints(1, 1) = gp_Pnt(-4,-4,5);
    aSurfacePoints(1, 2) = gp_Pnt(-4,-2,5);
    aSurfacePoints(1, 3) = gp_Pnt(-4,0,4);
    aSurfacePoints(1, 4) = gp_Pnt(-4,2,5);
    aSurfacePoints(1, 5) = gp_Pnt(-4,4,5);

    aSurfacePoints(2, 1) = gp_Pnt(-2,-4,4);
    aSurfacePoints(2, 2) = gp_Pnt(-2,-2,4);
    aSurfacePoints(2, 3) = gp_Pnt(-2,0,4);
    aSurfacePoints(2, 4) = gp_Pnt(-2,2,4);
    aSurfacePoints(2, 5) = gp_Pnt(-2,5,4);

    aSurfacePoints(3, 1) = gp_Pnt(0,-4,3.5);
    aSurfacePoints(3, 2) = gp_Pnt(0,-2,3.5);
    aSurfacePoints(3, 3) = gp_Pnt(0,0,3.5);
    aSurfacePoints(3, 4) = gp_Pnt(0,2,3.5);
    aSurfacePoints(3, 5) = gp_Pnt(0,5,3.5);

    aSurfacePoints(4, 1) = gp_Pnt(2,-4,4);
    aSurfacePoints(4, 2) = gp_Pnt(2,-2,4);
    aSurfacePoints(4, 3) = gp_Pnt(2,0,3.5);
    aSurfacePoints(4, 4) = gp_Pnt(2,2,5);
    aSurfacePoints(4, 5) = gp_Pnt(2,5,4);

    aSurfacePoints(5, 1) = gp_Pnt(4,-4,5);
    aSurfacePoints(5, 2) = gp_Pnt(4,-2,5);
    aSurfacePoints(5, 3) = gp_Pnt(4,0,5);
    aSurfacePoints(5, 4) = gp_Pnt(4,2,6);
    aSurfacePoints(5, 5) = gp_Pnt(4,5,5);

    GeomAPI_PointsToBSplineSurface aSurfaceApprox(aSurfacePoints);

    // construct the function.
    Handle_Adaptor3d_HCurve aAdaptorCurve = new GeomAdaptor_HCurve(aCurveApprox.Curve());
    Adaptor3d_Surface* aAdaptorSurface = new GeomAdaptor_Surface(aSurfaceApprox.Surface());

    Extrema_GlobOptFuncCS aFunction(&(aAdaptorCurve->Curve()), aAdaptorSurface);

    math_Vector aStartPoint(1, 3, 0.2);
    math_NewtonMinimum aSolver(aFunction, aStartPoint);
    aSolver.Perform(aFunction, aStartPoint);

    if (aSolver.IsDone())
    {
        aSolver.Dump(std::cout);
        math_Vector aLocation = aSolver.Location();

        gp_Pnt aPoint1 = aAdaptorCurve->Value(aLocation(1));
        gp_Pnt aPoint2 = aAdaptorSurface->Value(aLocation(2), aLocation(3));
        GC_MakeSegment aSegmentMaker(aPoint1, aPoint2);

        BRepBuilderAPI_MakeEdge anEdgeMaker(aSegmentMaker.Value());
        BRepTools::Write(anEdgeMaker.Shape(), "d:/tcl/min.brep");
    }

    GeomTools::Dump(aCurveApprox.Curve(), std::cout);
    GeomTools::Dump(aSurfaceApprox.Surface(), std::cout);

    BRepBuilderAPI_MakeEdge anEdgeMaker(aCurveApprox.Curve());
    BRepBuilderAPI_MakeFace aFaceMaker(aSurfaceApprox.Surface(), Precision::Approximation());

    BRepTools::Write(anEdgeMaker.Shape(), "d:/tcl/edge.brep");
    BRepTools::Write(aFaceMaker.Shape(), "d:/tcl/face.brep");

    // need release memory for the adaptor surface pointer manually.
    // whereas do not need release memory for the adaptor curve.
    // because it is mamanged by handle.
    delete aAdaptorSurface;
}

int main(int argc, char* argv[])
{
    testNewtonMethod();

    return 0;
}

上述代码通过拟合点得到的任意一曲线和曲面,计算其极值。其中在生成函数类时需要曲线曲面适配器的指针,这里分别使用了由Handle管理的类和无Handle管理的类来对比,说明Handle的用法。即Handle是个智能指针,和C#中的内存管理一样,只需要new,不用手工delete。而没用Handle管理的类,在new了之后,不再使用时,需要手工将内存释放。

生成BREP文件是为了便于在Draw Test Harness中查看显示效果,计算结果显示如下:

Figure 3.1 The minimum between a curve and a surface

如上图所示的青色的线即是曲线和曲面之间距离最短的线。

4.Conclusion

综上所述,在学习了最优化理论之后,应该结合实际进行应用。从OPEN CASCADE的计算曲线和曲面之间的极值的类中可以学习如何将实际问题抽象成数学模型,进而使用数学工具对问题进行求解。

理解了多元函数的Newton法求极值,再回过头去看看一元函数的求根或求极值问题,感觉应该会简单很多。如参数曲线曲面上根据三维点反求参数问题,就可以归结为一元函数和二元函数的极值问题,可以很快写出函数方程为如下:

同样可以应用Newton法来求极值。特别地,当曲线和曲面有交点时,那么极值点就是曲线和曲面的交点了。

通过学习和应用math_MultipleVarFunction及其子类,借鉴其将抽象数学概念用清晰的类来表示的方法,使程序便于理解,从而方便维护及扩展,提高程序质量。例如,若从类math_MultipleVarFunctionWithHessian类派生,所以要求函数至少具有C2连续,才能使用Newton方法。

5. References

1. http://mathworld.wolfram.com/NewtonsMethod.html

2. https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization

3. Shing Liu. OpenCASCADE Root-Finding Algorithm

4. http://tutorial.math.lamar.edu/Classes/CalcI/NewtonsMethod.aspx

5. 同济大学数学教研室. 高等数学. 高等教育出版社. 1996

6. 同济大学应用数学系. 线性代数. 高等教育出版社. 2003

7. 易大义, 陈道琦. 数值分析引论. 浙江大学出版社. 1998

8. 《运筹学》教材编写组. 运筹学. 清华大学出版社. 2012

9. 何坚勇. 最优化方法. 清华大学出版社. 2007

10. 杨庆之. 最优化方法. 科学出版社. 2015

11. 王宜举, 修乃华. 非线性最优化理论与方法. 科学出版社. 2012

12. 赵罡, 穆国旺, 王拉柱. 非均匀有理B样条. 清华大学出版社. 2010

PDF Version and Source Code Apply Newton Method to Find Extrema in OPEN CASCADE

时间: 2024-10-26 10:04:15

Apply Newton Method to Find Extrema in OPEN CASCADE的相关文章

js中的apply、call

问题描述 this.method().method.apply().method.call()有什么本质上的区别?在js中有只能用其中一两种形式的情况吗? 解决方案 解决方案二:不太清楚...解决方案三:眼长在手上面!解决方案四:搜一下不就什么都有了解决方案五:随便g.cn一下就很多了啊..http://blog.csdn.net/liuxiaoyi666/article/details/1839875

OpenCASCADE Root-Finding Algorithm

OpenCASCADE Root-Finding Algorithm eryar@163.com Abstract. A root-finding algorithm is a numerical method, or algorithm, for finding a value x such that f(x)=0, for a given function f. Such an x is called a root of the function f. In OpenCASCADE math

Refactoring Notes-Refactoring Methods(3)

5.Introduce Explaining Variable    If you have a complicated expression,put the result of the expression, or parts of the expression , in a temporary variable with a name that explains the purpose.   Introduce Explaining Variable is particulaly valua

AI人工智能专业词汇集

作为最早关注人工智能技术的媒体,机器之心在编译国外技术博客.论文.专家观点等内容上已经积累了超过两年多的经验.期间,从无到有,机器之心的编译团队一直在积累专业词汇.虽然有很多的文章因为专业性我们没能尽善尽美的编译为中文呈现给大家,但我们一直在进步.一直在积累.一直在提高自己的专业性. 两年来,机器之心编译团队整理过翻译词汇对照表「红宝书」,编辑个人也整理过类似的词典.而我们也从机器之心读者留言中发现,有些人工智能专业词汇没有统一的翻译标准,这可能是因地区.跨专业等等原因造成的.举个例子,Deep

Python求解平方根的方法_python

本文实例讲述了Python求解平方根的方法.分享给大家供大家参考.具体如下: 主要通过SICP的内容改写而来.基于newton method求解平方根.代码如下: #!/usr/bin/python def sqrt_iter(guess,x): if(good_enough(guess, x)): print guess else: sqrt_iter(improve(guess, x),x) def improve(guess, x): return average(guess, x/gue

Learning Machine Learning, Part 2: Algorithms and Techniques

The previous blog post, Introduction to Machine Learning, presented the Machine Learning concept. Now, let's discuss representative methods used in the technology. Regression Algorithms In most Machine Learning courses, regression algorithms are the

Codeforces Round #327 (Div. 2) B. Rebranding C. Median Smoothing

B. Rebranding The name of one small but proud corporation consists of n lowercase English letters. The Corporation has decided to try rebranding - an active marketing strategy, that includes a set of measures to change either the brand (both for the

[置顶].NET平台开源项目速览(13)机器学习组件Accord.NET框架功能介绍

    Accord.NET Framework是在AForge.NET项目的基础上封装和进一步开发而来.因为AForge.NET更注重与一些底层和广度,而Accord.NET Framework更注重与机器学习算法以及提供计算机视频.音频.信号处理以及统计应用相关的解决方案.该项目使用C#语言编写,项目主页:http://accord-framework.net/     说明:该文章只是一个基本介绍,主要内容是翻译的官方文档和介绍,部分英文表述个人能力有限,不太熟悉,所以直接照搬原文,有比较

深入JUnit源码之Runner

初次用文字的方式记录读源码的过程,不知道怎么写,感觉有点贴代码的嫌疑.不过中间还是加入了一些自己的理解和心得,希望以后能够慢慢的改进,感兴趣的童鞋凑合着看吧,感觉JUnit这个框架还是值得看的,里面有许多不错的设计思想在,更何况它是Kent Beck和Erich Gamma这样的大师写的..... 写在前面的话 不知道是因为第一份工作的影响还是受在博客园上看到的那句"源代码里没有秘密"的影响,总之,近来对很多框架的源码都很感兴趣,拿到一个都想看看.其实自从学习Java以来也看过不少了,