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