Текст подпрограммы и версий
mnk3r_p.zip 		Тексты тестовых примеров
tmnk3r_p.zip

Подпрограмма:  MNK3R (модуль MNK3R_p)

Назначение

Решение задачи минимизации диффеpенциpуемой функции многих переменных по заданному направлению (задача одномерной минимизации) методами линейной и квадратичной аппроксимации производной по направлению.

Математическое описание

Для решения задачи 

    min  f(x0 - h p0) ,   0 ≤ h ≤ hL < +∞ ,   x0 ∈ EN ,   p0 ∈ EN ,
      h

x0 и  p0 заданы, используются последовательно методы линейной и квадратичной аппроксимации производной по направлению  ( f ' (x0 - h p0), p0 ) .

Минимум  f (x) по направлению  p0 считается найденным, если выполнено одно из следующих условий:

1.  | hK - hK + 1| < EPS1 , где  hK, hK + 1 - два последовательных шага по направлению  p0,  EPS1 - заданная абсолютная точность одномерной минимизации;
2.  | hK - hK + 1| / hK < EPS2, где EPS2 - заданная относительная точность одномерной минимизации;
3.  | fi' (x0 - hK p0)| < EPSGi для всех  i = 1, 2, ..., N, где EPSG - вектоp покомпонентной точности вычисления градиента в точке  x0 - hK p0;
4.  
| ( fi' (x0 - hK p0), p0) | < max   { max  EPSGi, 10 -1 -EM/2} ,
                                                  i = 1, 2, ..., N
   где  EM - максимальный порядок абсолютной величины
   вещественной константы;
5.  K ≥ Kmax ,  где  Kmax - максимально допустимое количество вычислений градиента.

Алгоритм является модификацией алгоритма одномерной минимизации, изложенного в работе:

B.A.Cкоков, Стандартная программа минимизации диффеpенциpуемой функции многих переменных методом сопряженных градиентов, серия "Стандартные программы решения задач математического программирования", вып.10, Изд - во МГУ, 1968, ротапринт.

Использование 

procedure MNK3R(var N :Integer; var I1 :Array of Integer;
                var P :Array of Real; HB :Real; var HT :Real;
                var XT :Array of Real; var FT :Real; var FTE :Real;
                var G1 :Array of Real; var G1E :Array of Real;
                FUNC :Proc_F1_MN; GRAD :Proc_F2_MN; var UP :Array of Real;
                var WP :Array of Integer; var RM :Array of Real);


   Параметры
N - размерность пространства переменных (тип: целый);
I0 - целый вектоp длины  N, задающий фиксированные на время минимизации компоненты вектоpа  X: если  I0 (I) = 0, то  I - ая компонента вектоpа  X остается равной своему начальному значению;
P - вещественный вектоp длины  N, задающий направление одномерной минимизации  p0;
HL - правый конец отрезка одномерной минимизации  hL (тип: вещественный);
H - вещественная переменная, содержащая на входе заданный начальный шаг одномерного поиска, а на выходе - расстояние от начальной точки поиска до вычисленной точки минимума;
X - вещественный вектоp длины  N, содержащий на входе начальную точку  x0 поиска минимума  f (x) по направлению  p0, а на выходе - точку минимального вычисленного значения  f (x);
F - вещественная переменная, содержащая на входе значение функции в начальной точке  x0, а на выходе - вычисленное минимальное значение функции вдоль направления  p0;
FE - вещественная переменная, содержащая на входе заданную точность вычисления значений функции, а на выходе - точность вычисления функции в точке минимума;
G - вещественный вектоp длины  N, содержащий на входе градиент функции в начальной точке  x0, а на выходе - градиент функции в вычисленной точке минимума;
GE - вещественный вектоp длины  N, содержащий на входе заданную покомпонентную точность вычисления градиента, а на выходе - точность вычисления градиента в точке минимума;
FUNC - имя подпрограммы вычисления значения функции  f (x) (см. замечания по использованию);
GRAD - имя подпрограммы вычисления градиента функции (см. замечания по использованию);
UP - вещественный вектоp длины 6, задающий упpавляющие параметры алгоритма:
UP(1) - заданный признак; если  UP (1) = 2, то для построения линейной аппроксимации используется алгоритм удвоения пробных шагов; если  UP (1) = 1, то используется постоянный шаг (см. замечания по использованию);
UP(2) - заданная точность одномерной минимизации (относительная при  UP (3) = 1, абсолютная при  UP (3) ≠ 1);
UP(3) - признак:  UP (3) = 1, если задана относительная точность одномерной минимизации;  UP (3) ≠ 1, если задана абсолютная точность;
UP(4) - заданное максимально допустимое количество вычислений градиента;
UP(5) - максимальный порядок абсолютной величины вещественной константы;
UP(6) - максимальное количество десятичных значащих цифр в вещественной константе;
IP - целый вектоp длины 4, содержащий на выходе:
IP(1) - количество вычислений функции при одномерной минимизации;
IP(2) - количество вычислений градиента при одномерной минимизации;
IP(3) - признак: если  IP (3) = 0, то вычисление градиента  G в точке минимума не производилось; если  IP (3) = 1, то в  G содержится вычисленное значение градиента в точке минимума;
IP(4) - признак: если  IP (4) = 0, то одномерный минимум функции найден и лежит внутpи отрезка [0, HL]; если  IP (4) = 1, то минимум достигается в точке  HL; если  IP (4) = 2, то функция  f (x) не убывает по направлению  p0 и алгоритм не может обеспечить сходимость (см. замечания по использованию);
RM - вещественный вектоp длины 2 * N, используемый в подпрограмме как рабочий. 
   Версии:  нет 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
 

Подпрограммы FUNC и GRAD составляются пользователем.

Первый оператор подпрограммы вычисления функции должен иметь вид: 

      procedure FUNC (var X :Array of Real; var F :Real; FE :Real);

        Параметры
        X   - вещественный вектор длины  N, задающий точку
                пространства, в которой вычисляется значение функции;
        F   - вещественная переменная, содержащая
                вычисленное значение функции в точке  X;
        FE - вещественная переменная, содержащая на входе
                заданную точность вычисления значения функции
                в точке  X; на выходе - достигнутую точность.

Если значение достигнутой точности вычисления функции не известно, то в теле подпрограммы FUNC параметр FE не должен переопределяться.

Первый оператор подпрограммы вычисления градиента функции  f (x) должен иметь вид: 

      procedure GRAD (var X :Array of Real; var G :Array of Real;
                var GE :Array of Real; var I0 :Array of Integer);

        Параметры
        X   - вещественный вектор длины  N, задающий точку
                пространства, в которой вычисляется градиент;
        G   - вещественный вектоp длины  N, содержащий
                градиент функции в точке  X;
        GE - вещественный вектоp длины  N, содержащий  на
                входе заданную покомпонентную точность
                вычисления градиента функции; на выходе -
                достигнутую точность вычисления градиента;
        I0  - целый вектоp фиксированных компонент,
                упpавляющий вычислением компонент градиента:
                если  I0(I) = 0, то полагается  G(I) = 0 .

Если значение достигнутой точности GE (I) для некоторого  I не известно, то в теле подпрограммы GRAD параметр GE не должен переопределяться.

Если решается задача безусловной минимизации функции по направлению  p0, то следует задать  HL = C, где  C - достаточно большое представимое в машине вещественное число.

B общем случае наибольшее влияние на эффективность программы оказывает выбор начального шага одномерной минимизации  H и значения параметра  UP (1).

При одномерной минимизации функции  f (x) на направлении  p0 вначале ищется отрезок для построения линейной аппроксимации производной по направлению с помощью постоянного шага (при  UP (1) = 1) или алгоритма удвоения пробных шагов (при  UP (1) = 2). Затем последовательно осуществляется линейная и квадратичные аппроксимации производной по направлению.

Если функция многоэкстремальна, то при некоторых значениях  H алгоритм может не обеспечить сходимость к ближайшей к  x0 точке локального одномерного минимума. В этом случае, а также в случае неубывания функции по направлению  p0 (IP (4) = 2) следует уменьшить начальный шаг  H и задать  UP (1) = 1 .

Используются служебные подпрограммы: MNKU1, MNKU3, MNKU5, MNKU6, MNKO2, MNKO3.

Пример использования 

    min  F(x) ,  x ∈ E1 ,  x0 = 0.0 ,  p0 = -1.0 ,
    F(x)  =  100 e -x + x .
    Точка минимума  F(x) вдоль  p0 pавна   x*  =  - ln(0.01) ,
    F(x*)  =  1 - ln(0.01)


Unit TMNK3R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FMNK3R_p, FGMNK3R_p, MNK3R_p;

function TMNK3R: String;

implementation

function TMNK3R: String;
var
F :Real;
P :Array [0..0] of Real;
X :Array [0..0] of Real;
G :Array [0..0] of Real;
GE :Array [0..0] of Real;
RM :Array [0..1] of Real;
I0 :Array [0..0] of Integer;
IP :Array [0..3] of Integer;
const
N :Integer = 1;
H :Real = 0.5;
HL :Real = 1.0E+10;
FE :Real = 1.0E-07;
UP :Array [0..5] of Real = ( 2.0,1.0E-06,2.0,50.0,18.0,12.0 );
begin
Result := '';  { результат функции }
{ прототип оператора DАТА на FORTRANе  }
I0[0] := 1;
P[0] := -1.0;
X[0] := 0.0;
GE[0] := 1.0E-07;
FMNK3R(X,F,FE);
FGMNK3R(X,G,GE,I0);
MNK3R(N,I0,P,HL,H,X,F,FE,G,GE,FMNK3R,FGMNK3R,UP,IP,RM);
Result := Result + Format('%s',['  BЫЧИCЛ. ФYHKЦИИ']);
Result := Result + Format('%2d ',[IP[0]]);
Result := Result + Format('%s',['    BЫЧИCЛ. ГPAДИEHTA']);
Result := Result + Format('%3d ',[IP[1]]);
Result := Result + Format('%s',['   ' + #$0D#$0A]);
Result := Result + Format('%s',[' X=']);
Result := Result + Format('%20.16f ',[X[0]]);
Result := Result + Format('%s',['   F=']);
Result := Result + Format('%20.16f ',[F]) + #$0D#$0A;
UtRes('TMNK3R',Result);  { вывод результатов в файл TMNK3R.res }
exit;
end;

end.

Unit fmnk3r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

procedure fmnk3r(var X :Array of Real; var F :Real; FE :Real);

implementation

procedure fmnk3r(var X :Array of Real; var F :Real; FE :Real);
begin
F := 100.0*Exp(-X[0])+X[0];
end;

end.

Unit fgmnk3r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

procedure fgmnk3r(var X :Array of Real; var G :Array of Real;
                var GE :Array of Real; var I0 :Array of Integer);

implementation

procedure fgmnk3r(var X :Array of Real; var G :Array of Real;
                var GE :Array of Real; var I0 :Array of Integer);
begin
G[0] := -100.0*Exp(-X[0])+1.0;
end;

end.

Результаты счета:

     BЫЧИCЛЕНИЙ ФУHKЦИИ  1        BЫЧИCЛЕНИЙ ГPAДИEHTA  10

     X  =  4.6051700 + 00              F  =  5.6051700 + 00