Program SimplexMethod;

Uses Crt;

Type
  TFloat = Extended;

Const
  N_S = 3; { Максимальное число переменных }
  Max_Float = 1.0e+4932;

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:=4*sqr(X[1]-5)+sqr(X[2]-6);
  { OFunc:=2*sqr(X[1])+X[1]*X[2]+sqr(X[2]); }
End;

(*****
  Процедура Simplex.
  Оптимизация функции многих переменных методом Hелдера-Мида

  ***** Входные параметры: *****

  N - Число переменных;
  Eps - Точность определения минимума;
  X - Hа входе процедуры содержит начальное приближение к экстремуму;
  H - Шаг;
  IT - Допустимое число итераций;
  OFunc - Внешняя процедура оптимизируемой функции.

  ***** Выходные параметры: *****
  X - Точка экстремума;
  IT > 0 - Hормальное завершение;
  IT < 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.0; { Коэф. отражения }
  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
  ClrScr;
  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.