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

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.