📄 optim.pas
字号:
{ **********************************************************************
* Unit OPTIM.PAS *
* Version 2.4d *
* (c) J. Debord, February 2003 *
**********************************************************************
This unit implements the following methods for function minimization:
* Golden search for a function of one variable
* Simplex, Marquardt, BFGS for a function of several variables
**********************************************************************
References:
1) 'Numerical Recipes' by Press et al.
2) D. W. MARQUARDT, J. Soc. Indust. Appl. Math., 1963, 11, 431-441
3) J. A. NELDER & R. MEAD, Comput. J., 1964, 7, 308-313
4) R. O'NEILL, Appl. Statist., 1971, 20, 338-345
********************************************************************** }
unit optim;
interface
uses
fmath, matrices;
{ **********************************************************************
Error codes
********************************************************************** }
const
OPT_OK = 0; { No error }
OPT_SING = - 1; { Singular hessian matrix }
OPT_BIG_LAMBDA = - 2; { Too high Marquardt's parameter }
OPT_NON_CONV = - 3; { Non-convergence }
{ **********************************************************************
Functional types
********************************************************************** }
type
{ Procedure to compute gradient vector }
TGradient = procedure(Func : TFuncNVar;
X : TVector;
Lbound, Ubound : Integer;
G : TVector);
{ Procedure to compute gradient vector and hessian matrix }
THessGrad = procedure(Func : TFuncNVar;
X : TVector;
Lbound, Ubound : Integer;
G : TVector;
H : TMatrix);
{ **********************************************************************
Log file
********************************************************************** }
const
WriteLogFile : Boolean = False; { Write iteration info to log file }
LogFileName : String = 'optim.log'; { Name of log file }
{ **********************************************************************
Minimization routines
********************************************************************** }
function GoldSearch(Func : TFunc;
A, B : Float;
MaxIter : Integer;
Tol : Float;
var Xmin, Ymin : Float) : Integer;
{ ----------------------------------------------------------------------
Performs a golden search for the minimum of function Func
----------------------------------------------------------------------
Input parameters : Func = objective function
A, B = two points near the minimum
MaxIter = maximum number of iterations
Tol = required precision (should not be less than
the square root of the machine precision)
----------------------------------------------------------------------
Output parameters : Xmin, Ymin = coordinates of minimum
----------------------------------------------------------------------
Possible results : OPT_OK
OPT_NON_CONV
---------------------------------------------------------------------- }
function LinMin(Func : TFuncNVar;
X, DeltaX : TVector;
Lbound, Ubound : Integer;
MaxIter : Integer;
Tol : Float;
var F_min : Float) : Integer;
{ ----------------------------------------------------------------------
Minimizes function Func from point X in the direction specified by
DeltaX
----------------------------------------------------------------------
Input parameters : Func = objective function
X = initial minimum coordinates
DeltaX = direction in which minimum is searched
Lbound,
Ubound = indices of first and last variables
MaxIter = maximum number of iterations
Tol = required precision
----------------------------------------------------------------------
Output parameters : X = refined minimum coordinates
F_min = function value at minimum
----------------------------------------------------------------------
Possible results : OPT_OK
OPT_NON_CONV
---------------------------------------------------------------------- }
function Simplex(Func : TFuncNVar;
X : TVector;
Lbound, Ubound : Integer;
MaxIter : Integer;
Tol : Float;
var F_min : Float) : Integer;
{ ----------------------------------------------------------------------
Minimization of a function of several variables by the simplex method
of Nelder and Mead
----------------------------------------------------------------------
Input parameters : Func = objective function
X = initial minimum coordinates
Lbound,
Ubound = indices of first and last variables
MaxIter = maximum number of iterations
Tol = required precision
----------------------------------------------------------------------
Output parameters : X = refined minimum coordinates
F_min = function value at minimum
----------------------------------------------------------------------
Possible results : OPT_OK
OPT_NON_CONV
---------------------------------------------------------------------- }
procedure NumGradient(Func : TFuncNVar;
X : TVector;
Lbound, Ubound : Integer;
G : TVector);
{ ----------------------------------------------------------------------
Computes the gradient vector of a function of several variables by
numerical differentiation
----------------------------------------------------------------------
Input parameters : Func = function of several variables
X = vector of variables
Lbound,
Ubound = indices of first and last variables
----------------------------------------------------------------------
Output parameter : G = gradient vector
---------------------------------------------------------------------- }
procedure NumHessGrad(Func : TFuncNVar;
X : TVector;
Lbound, Ubound : Integer;
G : TVector;
H : TMatrix);
{ ----------------------------------------------------------------------
Computes gradient vector & hessian matrix by numerical differentiation
----------------------------------------------------------------------
Input parameters : as in NumGradient
----------------------------------------------------------------------
Output parameters : G = gradient vector
H = hessian matrix
---------------------------------------------------------------------- }
function Marquardt(Func : TFuncNVar;
HessGrad : THessGrad;
X : TVector;
Lbound, Ubound : Integer;
MaxIter : Integer;
Tol : Float;
var F_min : Float;
H_inv : TMatrix) : Integer;
{ ----------------------------------------------------------------------
Minimization of a function of several variables by Marquardt's method
----------------------------------------------------------------------
Input parameters : Func = objective function
HessGrad = procedure to compute gradient & hessian
X = initial minimum coordinates
Lbound,
Ubound = indices of first and last variables
MaxIter = maximum number of iterations
Tol = required precision
----------------------------------------------------------------------
Output parameters : X = refined minimum coordinates
F_min = function value at minimum
H_inv = inverse hessian matrix
----------------------------------------------------------------------
Possible results : OPT_OK
OPT_SING
OPT_BIG_LAMBDA
OPT_NON_CONV
---------------------------------------------------------------------- }
function BFGS(Func : TFuncNVar;
Gradient : TGradient;
X : TVector;
Lbound, Ubound : Integer;
MaxIter : Integer;
Tol : Float;
var F_min : Float;
H_inv : TMatrix) : Integer;
{ ----------------------------------------------------------------------
Minimization of a function of several variables by the
Broyden-Fletcher-Goldfarb-Shanno method
----------------------------------------------------------------------
Parameters : Gradient = procedure to compute gradient vector
Other parameters as in Marquardt
----------------------------------------------------------------------
Possible results : OPT_OK
OPT_NON_CONV
---------------------------------------------------------------------- }
implementation
var
Eps : Float; { Fractional increment for numer. derivation }
X1 : TVector; { Initial point for line minimization }
DeltaX1 : TVector; { Direction for line minimization }
Xt : TVector; { Minimum found by line minimization }
Lbound1, Ubound1 : Integer; { Bounds of X1 and DeltaX1 }
LinObjFunc : TFuncNVar; { Objective function for line minimization }
LogFile : Text; { Stores the result of each minimization step }
procedure MinBrack(Func : TFunc; var A, B, C, Fa, Fb, Fc : Float);
{ ----------------------------------------------------------------------
Given two points (A, B) this procedure finds a triplet (A, B, C)
such that:
1) A < B < C
2) A, B, C are within the golden ratio
3) Func(B) < Func(A) and Func(B) < Func(C).
The corresponding function values are returned in Fa, Fb, Fc
---------------------------------------------------------------------- }
begin
if A > B then
Swap(A, B);
Fa := Func(A);
Fb := Func(B);
if Fb > Fa then
begin
Swap(A, B);
Swap(Fa, Fb);
end;
C := B + GOLD * (B - A);
Fc := Func(C);
while Fc < Fb do
begin
A := B;
B := C;
Fa := Fb;
Fb := Fc;
C := B + GOLD * (B - A);
Fc := Func(C);
end;
if A > C then
begin
Swap(A, C);
Swap(Fa, Fc);
end;
end;
function GoldSearch(Func : TFunc;
A, B : Float;
MaxIter : Integer;
Tol : Float;
var Xmin, Ymin : Float) : Integer;
var
C, Fa, Fb, Fc, F1, F2, MinTol, X0, X1, X2, X3 : Float;
Iter : Integer;
begin
MinTol := Sqrt(MACHEP);
if Tol < MinTol then Tol := MinTol;
MinBrack(Func, A, B, C, Fa, Fb, Fc);
X0 := A;
X3 := C;
if (C - B) > (B - A) then
begin
X1 := B;
X2 := B + CGOLD * (C - B);
F1 := Fb;
F2 := Func(X2);
end
else
begin
X1 := B - CGOLD * (B - A);
X2 := B;
F1 := Func(X1);
F2 := Fb;
end;
Iter := 0;
while (Iter <= MaxIter) and (Abs(X3 - X0) > Tol * (Abs(X1) + Abs(X2))) do
if F2 < F1 then
begin
X0 := X1;
X1 := X2;
F1 := F2;
X2 := X1 + CGOLD * (X3 - X1);
F2 := Func(X2);
Inc(Iter);
end
else
begin
X3 := X2;
X2 := X1;
F2 := F1;
X1 := X2 - CGOLD * (X2 - X0);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -