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;
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;