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