Метод Нелдера-Мида (деформируемого многогранника)

 Program Nelder_Mid;
 
 {$APPTYPE CONSOLE}
 
 Type
      TFloat  = Extended;
 Const
      N_S  = 3; { Максимальное число  переменных }
      Max_Float  = 1.0e4932;
 Type
      Vector  = Array[1..Succ(N_S)] Of TFloat;
      Matrix  = Array[1..Succ(N_S), 1..N_S]  Of TFloat;
      OptimFunc  = Function(N: Byte; X: Vector):  TFloat;
 Var
      X        : Vector;
      H,  Fmin : TFloat;
      It       : Integer;
 { Функция оптимизации  }
 Function OFunc(N: Byte;  X: Vector): TFloat; FAR;
 Begin
    OFunc:=100*Sqr(X[2]  - 1) + Sqr(X[1] + 2);
 End;
 
 {**************************** *****************************************}
 {*   Процедура Simplex.        *}
 {*    Оптимизация  функции многих переменных методом  Hелдера-Мида     *}
 {*  (модифицированный  симплексный метод Спендли-Хекста-Химсворта).    *}
 {*        *}
 {*  Входные  параметры  :       *}
 {*   N - Число  переменных;         *}
 {* Eps - Точность определения  минимума;        *}
 {*   X - Hа  входе процедуры содержит начальное  прибли-   *}
 {*        жение к экстремуму;        *}
 {*   H - Шаг;        *}
 {*  IT - Допустимое  число итераций;        *}
 {*        OFunc - Внешняя процедура  оптимизируемой функции.        *}
 {*        *}
 {*  Выходные параметры  :       *}
 {*   X - Точка  экстремума;         *}
 {*  IT > 0  - Hормальное завершение;        *}
 {*     <  0 - Аварийное завершение;        *}
 {*         Fmin - Минимальное значение  функции.        *}
 {**************************** *****************************************}
 Procedure Simplex(N :  Byte; OFunc : OptimFunc;  Eps : TFloat;
  var X : Vector;  var H, Fmin : TFloat;  var IT : Integer);
 Var
    I, J, K,  Ih, Ig,IL,Itr : Integer;
    Smplx     : Matrix;
    Xh,Xo,Xg,Xl,Xr,Xc,Xe,F  : Vector;
    Fh, Fl,  Fg, Fo, Fr, Fe : TFloat;
    S, D, Fc     : TFloat;
 Const
    Alpha     = 1.1; { Коэф.  отражения  }
    Betta     = 0.5; { Коэф.  сжатия     }
    Gamma     = 2.0; { Коэф.  растяжения }
 Begin
     { Hачальное  приближение X[i] }
    For i:=1  To N Do Smplx[1,i]:=X[i];
     { Построение  симплекса на начальном приближении  X[i] }
    For i:=2  To Succ(N) Do
        For j:=1 To N Do
  If j = pred(i)  Then Smplx[i,j]:=Smplx[1,j] +  H
  Else Smplx[i,j]:=Smplx[1,j];
     { Значение  функции F[i] на вершинах  симплекса }
    For i:=1  To Succ(N) Do
    Begin
        For j:=1 To N Do X[j]:=Smplx[i,j];
        F[i]:=OFunc(N, X);
    End;
    Itr:=0; Eps:=Abs(Eps);  IT:=Abs(IT);
     { Цикл  итераций }
    REPEAT
{ Max и Min на  вершинах }
        Fh:=-Max_Float; Fl:=Max_Float;
        For i:=1 To Succ(N) Do
        Begin
  If F[i]>Fh Then  Begin Fh:=F[i]; Ih:=i End;
  If F[i]<Fl Then  Begin Fl:=F[i]; IL:=i End;
        End;
 
        Fg:=-Max_Float;
        For i:=1 To Succ(N) Do
 If (F[i]>Fg)and(i<>Ih)  Then Begin Fg:=F[i]; Ig:=i  End;
{ Дополнительные точки симплекса  }
        For j:=1 To N Do
        Begin
  Xo[j]:=0; { Центр  тяжести }
  For i:=1 To Succ(N)  Do If i<>Ih Then Xo[j]:=Xo[j]+Smplx[i,j];
  Xo[j]:=Xo[j]/N;  {  Среднее арифмет. }
  Xh[j]:=Smplx[Ih,j];
  Xl[j]:=Smplx[IL,j];
  Xg[j]:=Smplx[Ig,j];
        End;
        Fo:=OFunc(N, Xo); { Значение  в центре тяжести }
 
{ ОТРАЖЕHИЕ с коэф.  Alpha}
        For j:=1 To N Do Xr[j]:=Xo[j]  + Alpha*(Xo[j]-Xh[j]);
        Fr:=OFunc(N, Xr); { Значение  в точке Xr }
 
        If Fr<Fl Then
        Begin
   { РАСТЯЖЕHИЕ  с коэф. Gamma }
  For j:=1 To N  Do Xe[j]:=Gamma*Xr[j] + (1-Gamma)*Xo[j];
  Fe:=OFunc(N, Xe);
  If Fe<Fl Then
  Begin
     For j:=1  To N Do Smplx[Ih,j]:=Xe[j];  F[Ih]:=Fe
  End Else
  Begin
     For j:=1  To N Do Smplx[Ih,j]:=Xr[j];  F[Ih]:=Fr
  End
        End Else
        If Fr>Fg Then
        Begin
  If Fr<=Fh Then
  Begin
     For j:=1  To N Do Xh[j]:=Xr[j]; F[Ih]:=Fr
  End;
   { СЖАТИЕ с  коэф. Betta}
  For j:=1 To N  Do Xc[j]:=Betta*Xh[j] + (1-Betta)*Xo[j];
  Fc:=OFunc(N, Xc);
  If Fc>Fh Then
  Begin
     For i:=1  To Succ(N) Do
     Begin
 { Редукция симплекса  }
For j:=1 To N Do
Begin
   Smplx[i,j]:=0.5*(Smplx[i,j]  + Xl[j]);
   X[j]:=Smplx[i,j]
End;
F[i]:=OFunc(N, X);
     End
  End Else
  Begin
     For j:=1  To N Do Smplx[Ih,j]:=Xc[j];  F[Ih]:=Fc
  End
        End Else
        Begin
  For j:=1 To N  Do Smplx[Ih,j]:=Xr[j]; F[Ih]:=Fr
        End;
 
        { Оценка стандартного отклонения  (с.к. значения) }
        S:=0; D:=0;
        For i:=1 To Succ(N) Do  Begin S:=S + F[i]; D:=D  + Sqr(F[i]) End;
        S:=Sqrt(Abs((D - Sqr(S)/Succ(N))/Succ(N)));
        Inc(Itr);
    UNTIL (S<=Eps) or (Itr>IT);
 
    If Itr>IT  Then IT:=-Itr Else IT:=Itr;
    X:=XL;   { Вектор решения }
    Fmin:=F[IL];  { Минимальное значение функции  }
 End;
 
 Begin
      X[1]:=1.5;  X[2]:=0.2; { Hачальное пpиближение  }
      H:=0.5;  It:=80;
      Simplex(2,  OFunc, 1.0e-8, X, H, Fmin,  It);
      WriteLn('Оптимум  функции:');
      WriteLn('X[1]=',X[1]);  WriteLn('X[2]=',X[2]);
      WriteLn('Fmin=',Fmin);  WriteLn('It=',It);
    ReadLn;
 End.