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法作為一種經(jīng)典的解無(wú)約束優(yōu)化問(wèn)題的方法，在20世紀(jì)80~90年代發(fā)展起來(lái)的解線性規(guī)劃和凸規(guī)劃的內(nèi)點(diǎn)法中起到了重要作用。Newton法最初是Newton提出用于解非線性方程的，Newton曾用該法求解Kepler方程x-asinx=b，并得到精度很高的近似解。通過(guò)《OPEN CASCADE Multiple Variable Function》對(duì)OPEN CASCADE中多元函數(shù)的表達(dá)有了一個(gè)認(rèn)識(shí)。多元函數(shù)如何應(yīng)用的呢？下面提出一個(gè)問(wèn)題及如何用程序來(lái)解決這個(gè)問(wèn)題。對(duì)于任意給定的曲線和曲面，如何求出曲線和曲面上距離最近的點(diǎn)，假設(shè)曲線和曲面都是至少C2連續(xù)的。關(guān)于參數(shù)連續(xù)性可參考《OPEN CASCADE Curve Continuity》。如下圖所示：
Figure 1.1 A Curve and A Surface
本文給出OPEN CASCADE中對(duì)此類問(wèn)題的一種解法，即應(yīng)用Newton法求解非線性無(wú)約束多元函數(shù)的極值。學(xué)習(xí)如何將實(shí)際問(wèn)題抽象成數(shù)學(xué)模型，從而使用數(shù)學(xué)的方法進(jìn)行求解。
2.Construct Function
在實(shí)際應(yīng)用中，一個(gè)問(wèn)題是不是可以表述為一個(gè)最優(yōu)化模型和怎么表示為一個(gè)最優(yōu)化模型，這是優(yōu)化方法是否可以應(yīng)用的前提，因而十分重要。但優(yōu)化問(wèn)題的建模和其他數(shù)學(xué)問(wèn)題的建模一樣，不屬于精確科學(xué)或數(shù)學(xué)的范疇，而是一項(xiàng)技術(shù)或技藝，沒(méi)有統(tǒng)一標(biāo)準(zhǔn)和方法。當(dāng)然建立的模型是否正確和模型的優(yōu)劣是可以通過(guò)實(shí)際效果來(lái)檢驗(yàn)的。
OPEN CASCADE使用從math_MultipleVarFunctionWithHessian派生的一個(gè)具體類Extrema_GlobOptFuncCS來(lái)計(jì)算C2連續(xù)的曲線和曲面之間的距離的平方值。抽象出來(lái)的數(shù)學(xué)模型為：
因?yàn)槭菑木哂蠬essian Matrix的多元函數(shù)派生，所以要求曲線曲面具有至少C2連續(xù)，即有至少有二階導(dǎo)數(shù)。且在類中分別實(shí)現(xiàn)計(jì)算函數(shù)值，計(jì)算一階導(dǎo)數(shù)值（梯度），計(jì)算二階導(dǎo)數(shù)值（Hessian Matrix）。計(jì)算函數(shù)值的代碼如下所示：
//======================================================================
//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));
}
其中參數(shù)cu為參數(shù)曲線的參數(shù)，su,sv分別為參數(shù)曲面的參數(shù)。根據(jù)參數(shù)曲線曲面上的參數(shù)計(jì)算出相應(yīng)的點(diǎn)，然后計(jì)算出兩個(gè)點(diǎn)之間的距離的平方值即為函數(shù)值。與上述公式對(duì)應(yīng)。
根據(jù)多元函數(shù)一階導(dǎo)數(shù)（梯度）的定義，可得出梯度的計(jì)算公式如下：
計(jì)算梯度的代碼如下所示，從程序代碼可見(jiàn)程序就是上公式的具體實(shí)現(xiàn)。
//======================================================================
//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();
}
根據(jù)Hessian Matrix的定義，得到計(jì)算Hessian Matrix的公式如下：
將函數(shù)積的求導(dǎo)法則應(yīng)用于求偏導(dǎo)數(shù)得到上述公式。同理求出Hessian Matrix的其他各項(xiàng)，如下公式所示：
計(jì)算多元函數(shù)的二階導(dǎo)數(shù)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();
}
根據(jù)高階偏導(dǎo)數(shù)的定理可知，當(dāng)f(X)在點(diǎn)X0處所有二階偏導(dǎo)數(shù)連續(xù)時(shí)，那末在該區(qū)域內(nèi)這兩個(gè)二階混合偏導(dǎo)數(shù)必相等。所以Hessian Matrix為一個(gè)對(duì)稱矩陣，故
H(2,1)=H(1,2)
H(3,1)=H(1,3)
H(3,2)=H(2,3)
由此完成一個(gè)具有二階偏導(dǎo)數(shù)的多元函數(shù)的數(shù)學(xué)模型，用面向?qū)ο蟮姆绞角逦乇硎玖顺鰜?lái)。和我見(jiàn)過(guò)的國(guó)內(nèi)一些程序相比，這種抽象思路還是很清晰，便于程序的理解和擴(kuò)展。國(guó)內(nèi)有個(gè)圖形庫(kù)是C風(fēng)格的，一個(gè)函數(shù)幾千行，光函數(shù)參數(shù)就很多，參數(shù)名也很隨意，函數(shù)內(nèi)部變量名稱更是無(wú)法理解，什么i，j，k，ii，jj之類。這種程序別說(shuō)可擴(kuò)展，就是維護(hù)起來(lái)也是讓人頭疼的啊！
有了函數(shù)表達(dá)式，下面就是計(jì)算這個(gè)函數(shù)的極值了。
3.Newton’s Method
關(guān)于應(yīng)用Newton法計(jì)算一元非線性方程的根已經(jīng)在《OpenCASCADE Root-Finding Algorithm》中進(jìn)行了說(shuō)明，這里要學(xué)習(xí)下如何使用Newton法應(yīng)用于多元函數(shù)極值的計(jì)算。對(duì)于一元函數(shù)f(x)的求極值問(wèn)題，當(dāng)f(x)連續(xù)可微時(shí)，最優(yōu)點(diǎn)x滿足f’(x)=0。于是當(dāng)f(x)二次連續(xù)可微時(shí)，求解f’(x)=0的Newton法為：
該方法稱為解無(wú)約束優(yōu)化問(wèn)題的Newton方法。由《數(shù)學(xué)分析》可知，當(dāng)f(x)是凸函數(shù)時(shí)，f’(x)=0的解是minf(x)的整體最優(yōu)解。將Newton法擴(kuò)展到多元函數(shù)的情況，也是利用二階Taylor級(jí)數(shù)將函數(shù)展開(kāi)，得到多元函數(shù)的極值迭代公式：
關(guān)于多元函數(shù)Newton法公式的推導(dǎo)，可參考《最優(yōu)化方法》等書(shū)籍。Newton法的算法步驟如下：
A. 給定初始點(diǎn)，及精度；
B. 計(jì)算函數(shù)f(xk)的一階導(dǎo)數(shù)（梯度），二階導(dǎo)數(shù)（Hessian Matrix）：若|梯度|<精度，則停止迭代，輸出近似極小點(diǎn)；否則轉(zhuǎn)C；
C. 根據(jù)Newton迭代公式，計(jì)算x(k+1)；
OPEN CASCADE中Newton法計(jì)算極值的類是math_NewtonMinimum，可參考其代碼學(xué)習(xí)。下面給出前面提出的曲線曲面極值求解的實(shí)現(xiàn)代碼：
/*
* 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;
}
上述代碼通過(guò)擬合點(diǎn)得到的任意一曲線和曲面，計(jì)算其極值。其中在生成函數(shù)類時(shí)需要曲線曲面適配器的指針，這里分別使用了由Handle管理的類和無(wú)Handle管理的類來(lái)對(duì)比，說(shuō)明Handle的用法。即Handle是個(gè)智能指針，和C#中的內(nèi)存管理一樣，只需要new，不用手工delete。而沒(méi)用Handle管理的類，在new了之后，不再使用時(shí)，需要手工將內(nèi)存釋放。
生成BREP文件是為了便于在Draw Test Harness中查看顯示效果，計(jì)算結(jié)果顯示如下：
Figure 3.1 The minimum between a curve and a surface
如上圖所示的青色的線即是曲線和曲面之間距離最短的線。
4.Conclusion
綜上所述，在學(xué)習(xí)了最優(yōu)化理論之后，應(yīng)該結(jié)合實(shí)際進(jìn)行應(yīng)用。從OPEN CASCADE的計(jì)算曲線和曲面之間的極值的類中可以學(xué)習(xí)如何將實(shí)際問(wèn)題抽象成數(shù)學(xué)模型，進(jìn)而使用數(shù)學(xué)工具對(duì)問(wèn)題進(jìn)行求解。
理解了多元函數(shù)的Newton法求極值，再回過(guò)頭去看看一元函數(shù)的求根或求極值問(wèn)題，感覺(jué)應(yīng)該會(huì)簡(jiǎn)單很多。如參數(shù)曲線曲面上根據(jù)三維點(diǎn)反求參數(shù)問(wèn)題，就可以歸結(jié)為一元函數(shù)和二元函數(shù)的極值問(wèn)題，可以很快寫(xiě)出函數(shù)方程為如下：
同樣可以應(yīng)用Newton法來(lái)求極值。特別地，當(dāng)曲線和曲面有交點(diǎn)時(shí)，那么極值點(diǎn)就是曲線和曲面的交點(diǎn)了。
通過(guò)學(xué)習(xí)和應(yīng)用math_MultipleVarFunction及其子類，借鑒其將抽象數(shù)學(xué)概念用清晰的類來(lái)表示的方法，使程序便于理解，從而方便維護(hù)及擴(kuò)展，提高程序質(zhì)量。例如，若從類math_MultipleVarFunctionWithHessian類派生，所以要求函數(shù)至少具有C2連續(xù)，才能使用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. 同濟(jì)大學(xué)數(shù)學(xué)教研室. 高等數(shù)學(xué). 高等教育出版社. 1996
6. 同濟(jì)大學(xué)應(yīng)用數(shù)學(xué)系. 線性代數(shù). 高等教育出版社. 2003
7. 易大義, 陳道琦. 數(shù)值分析引論. 浙江大學(xué)出版社. 1998
8. 《運(yùn)籌學(xué)》教材編寫(xiě)組. 運(yùn)籌學(xué). 清華大學(xué)出版社. 2012
9. 何堅(jiān)勇. 最優(yōu)化方法. 清華大學(xué)出版社. 2007
10. 楊慶之. 最優(yōu)化方法. 科學(xué)出版社. 2015
11. 王宜舉, 修乃華. 非線性最優(yōu)化理論與方法. 科學(xué)出版社. 2012
12. 趙罡, 穆國(guó)旺, 王拉柱. 非均勻有理B樣條. 清華大學(xué)出版社. 2010