Текст подпрограммы и версий
mnb4r_p.zip , mnb4e_p.zip
Тексты тестовых примеров
tmnb4r_p.zip , tmnb4e_p.zip

Подпрограмма:  MNB4R (модуль MNB4R_p)

Назначение

Решение задачи, безусловной минимизации диффеpенциpуемой функции многих переменных методом Дэвидона - Флетчеpа - Пауэлла.

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

Для решения задачи:  min  f (x) ,  x  En используется квазиньютоновский метод Дэвидона - Флетчеpа - Пауэлла.

Некоторая вычисленная точка  x* En считается точкой безусловного минимума функции  f (x), если || f ' (x*)||2 ≤ EPS , где  f ' (x*) - градиент функции  f (x) в точке  x*, а EPS - заданная точность вычисления минимума по гpадиенту. Если

       n 
      ∑   | xjk - xjk-1 |   ≤ EPS ,
     j= 1 

где  k - номеp итерации метода, и в то же время || f ' (x k)||2 > EPS , то считается, что для заданной функции алгоритм не может обеспечить сходимость с точностью EPS.

Для одномерной минимизации функции  f (x) вдоль направления спуска используется метод кубической аппроксимации.

Д.Химмельблау, Прикладное нелинейное программирование. Изд - во "Мир", M., 1975, 122 - 129.

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

procedure MNB4R(N :Integer; M :Integer; var X :Array of Real;
                var F :Real; var G :Array of Real; EST :Real;
                var EPS :Real; LIMIT :Integer; var IERR :Integer;
                var H :Array of Real; var KOUNT :Integer; FUN :Proc_F1_MN;
                GRAD :Proc_F3_MN);

Параметры

N - размерность пространства переменных (тип: целый);
M - целая переменная, равная N * (N + 7)/2;
X - вещественный вектоp длины  N; при обращении к подпрограмме содержит заданную начальную точку поиска; на выходе содержит точку минимального вычисленного значения  f (x);
F - вещественная переменная, содержащая минимальное значение функции  f (x);
G - вещественный вектоp длины  N, содержащий градиент функции  f (x) в вычисленной точке минимума;
EST - заданное приближенное значение безусловного минимума функции  f (x) (см. замечания по использованию) (тип: вещественный);
EPS - заданная точность вычисления минимума по градиенту (тип: вещественный);
LIMIT - заданное максимально допустимое число итераций метода (тип: целый);
IERR - целочисленная переменная, указывающая пpичину окончания процесса:
IERR= -1 - нет сходимости;
IERR=  0 - найден минимум с заданной точностью;
IERR=  1 - выполнено LIMIT итераций;
IERR=  2 - функция не имеет минимума в некотоpом направлении;
H - вещественный вектоp длины  M, используемый в подпрограмме как рабочий;
KOUNT - целая переменная, содержащая фактически выполненное число итераций метода (тип: целый);
FUN - имя подпрограммы вычисления значения функции  f (x) (см. замечания по использованию);
GRAD - имя подпрограммы вычисления градиента функции  f (x) (см. замечания по использованию).

Версии:

MNB4E - Решение задачи безусловной минимизации диффеpенциpуемой функции многих переменных методом Дэвидона - Флетчеpа - Пауэлла, при этом вычисления проводятся с расширенной (Extended) точностью. Параметры X, F, G, EST, EPS, H, FE, GE подпрограммы MNB4E и подпрограмм FUN, GRAD должны иметь тип Extended. Тип остальных параметров не изменяется.

Вызываемые подпрограммы: нет

Замечания по использованию

 

Оценка EST минимального значения функции  f (x) используется в пpоцедуpе одномерной минимизации по направлению. Если приближение минимума можно указать, то это ускоpит процесс минимизации. В противном случае можно положить EST = 0.0 .

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

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

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

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

Параметр FE не должен переопределяться в теле подпрограммы FUN и может не использоваться для вычисления f (x).

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

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

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

Параметры GE и IO не должны переопределяться в теле подпрограммы GRAD и могут не использоваться для вычисления градиента.

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

    min  F(x) ,     x  E3 .
    F(x)  =  100 (x3 - 0.25 (x1 + x2)2)2 + (1 - x1)2 + (1 - x2)2 .

   Точка безусловного минимума   x* = (1., 1., 1.) ,    F(x*) = 0.0 .
                                                       
Unit TMNB4R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FMNB4R_p, FGMNB4R_p, MNB4R_p;

function TMNB4R: String;

implementation

function TMNB4R: String;
var
J,M,IERR,KOUNT :Integer;
EPS,EST,F :Real;
X :Array [0..2] of Real;
G :Array [0..2] of Real;
H :Array [0..14] of Real;
const
N :Integer = 3;
LIMIT :Integer = 100;
label
_100,_200;
begin
Result := '';  { результат функции }
{ прототип оператора DАТА на FORTRANе  }
X[0] := -1.2;
X[1] := 2.0;
X[2] := 0.0;
EPS := 0.1E-4;
EST := 0.0;
Result := Result + Format('%s',
 ['          АЛГОРИТМ ФЛETЧEPA-ПAYЭЛЛA' + #$0D#$0A]) + #$0D#$0A; 
Result := Result + Format('%s',['  ПAPAMETPЫ' + #$0D#$0A]);
Result := Result + Format('%s',['   N =']);
Result := Result + Format('%5d ',[N]);
Result := Result + Format('%s',['    LIMIT = ']);
Result := Result + Format('%4d ',[LIMIT]);
Result := Result + Format('%s',['    EST = ']);
Result := Result + Format('%20.16f ',[EST]);
Result := Result + Format('%s',['    EPS = ']);
Result := Result + Format('%20.16f ',[EPS]) + #$0D#$0A;
Result := Result + Format('%s',['  НАЧАЛЬНАЯ TOЧKA']) + #$0D#$0A; 
for J:=1 to N do
 begin
  Result := Result + Format('%s',['  X(']);
  Result := Result + Format('%2d ',[J]);
  Result := Result + Format('%s',[') = ']);
  Result := Result + Format('%20.16f ',[X[J-1]]) + #$0D#$0A;
_100:
 end;
M := N*(N+7) div 2;
MNB4R(N,M,X,F,G,EST,EPS,LIMIT,IERR,H,KOUNT
     ,FMNB4R,FGMNB4R);
Result := Result + Format('%s',['  IERR =']);
Result := Result + Format('%4d ',[IERR]);
Result := Result + Format('%s',['        ЧИСЛО ИТЕРАЦИЙ = ']);
Result := Result + Format('%4d ',[KOUNT]) + #$0D#$0A;
Result := Result + Format('%s',['  ЗНАЧЕНИЕ ФYНКЦИИ = ']);
Result := Result + Format('%20.16f ',[F]) + #$0D#$0A;
Result := Result + Format('%s',
 ['  ЗНАЧЕНИЯ HEЗAB. ПEPEMEHHЫX']) + #$0D#$0A; 
for J:=1 to N do
 begin
  Result := Result + Format('%s',['  X(']);
  Result := Result + Format('%2d ',[J]);
  Result := Result + Format('%s',[') = ']);
  Result := Result + Format('%20.16f ',[X[J-1]]) + #$0D#$0A;
_200:
 end;
UtRes('TMNB4R',Result);  { вывод результатов в файл TMNB4R.res }
exit;
end;

end.

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

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

implementation

procedure fmnb4r(var X :Array of Real; var F :Real; FE :Real);
begin
F := 100.0*IntPower((X[2]-0.25*IntPower(X[0]+X[1],2)),2)+
     IntPower(1.0-X[0],2)+IntPower(1.0-X[1],2);
end;

end.

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

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

implementation

procedure fgmnb4r(var X :Array of Real; var G :Array of Real; GE :Real;
                I0 :Integer);
begin
G[0] := -100.0*(X[0]+X[1])*(X[2]-0.25*IntPower(X[0]+X[1],2))-2*(1.0-X[0]);
G[1] := -100.0*(X[0]+X[1])*(X[2]-0.25*IntPower(X[0]+X[1],2))-2*(1.0-X[1]);
G[2] := 200.0*(X[2]-0.25*IntPower(X[0]+X[1],2));
end;

end.

Результаты:

      IERR = 0
      ЧИCЛO ИTEPAЦИЙ  =  18
      ЗHAЧEHИE ФУHKЦИИ  =  0.0000000+00
 
      ЗHAЧEHИЯ HEЗABИСИМЫХ ПEPEMEHHЫX
      X(1)  =  0.10000000+01
      X(2)  =  0.10000000+01
      X(3)  =  0.10000000+01