i-808190579

subplot(1, 2, 2)

%Построение поверхности целевой функции (трехмерный вид) z = zeros(length(x1), length(x2));

for i = 1 : length(x1)

for j = 1 : length(x2)

z(i, j) = Q([x1(i, j) x2(i, j)]);

end

end

surf(x1, x2, z, 'EdgeColor', 'flat'); hold on;

%Построение траектории поиска (трехмерный вид) z = zeros(size(X, 1), 1);

for i = 1 : size(X, 1) z(i) = Q(X(i, :));

end

plot3(X(:, axes(1)), X(:, axes(2)), z, '.-g'); hold off;

function [ ret ] = Q ( x )

%Целевая функция. Задавается пользователем. a0 = 2;

a = [-2.8 -2.0 0.5 1 0.5];

ret = a0 + a(1) * x(1) + a(2) * x(2) + a(3) * x(1) * x(2) + ...

a(4) * x(1)^2 + a(5) * x(2)^2;

end

function [ ret ] = fogr ( x )

%Функциональные ограничения. Задаются пользователем. b = [-1 -1];

b1 = -4;

ret = sum(b .* x) - b1; end

function [ X, Y ] = SimplexSearch ( kp, n, Q, X0, xmin, xmax, fogr )

%Точность поиска (задать) nu = 0.01;

%Число итераций для "вписывания" исходного симплекса в ограничения задачи %(задать)

N1 = 5;

%Расчѐт исходного и конечного ребра симплекса

[L0, LN, rmax] = CalculateEdges(n, xmin, xmax, nu);

%Формирование функции проверки позиционных ограничений (возвращает %количество нарушенных позиционных ограничений)

poz = @(x) CheckPosConstraints(n, xmin, xmax, x);

%Формирование новой целевой функции Q1. Q1(x) = -Q(x), если требуется найти

%минимум Q; иначе Q1(x) = Q(x). if kp == 0

Q1 = @(x) -Q(x);

else

Q1 = Q;

end

141

%Вызов основной процедуры поиска

[X, Y] = SimplexSearchBody(n, Q1, X0, xmin, xmax, poz, fogr, N1, L0, ...

LN, rmax);

%Значения функции Q1 умножаются на -1, если требуется найти минимум Q if kp == 0

Y = -Y;

end end

function [ L0, LN, rmax ] = CalculateEdges ( n, xmin, xmax, nu )

%Процедура расчѐта исходного и конечного ребра симплекса, а также %расстояния между точками xmin и xmax.

d = xmax - xmin; rmin = min(d);

L0 = (rmin / (2 + n)) * sqrt(n * (n + 1) / 2); rmax = sqrt(sum(d .^ 2));

LN = nu * sqrt(2 * n / (n + 1)); end

function [ ret ] = CheckPosConstraints ( n, xmin, xmax, x )

%Функция проверки позиционных ограничений. Возвращает количество нарушенных %позиционных ограничений.

ret = 0;

for i = 1 : n

if x(i) > xmax(i) || x(i) < xmin(i) ret = ret + 1;

end

end end

function [ X, Y ] = SimplexSearchBody ( n, Q, X0, xmin, xmax, poz, ...

fogr, N1, L0, LN, rmax )

%Расчет числа шагов поиска

N = 8 * rmax / L0; N = ceil(N);

%Расчет параметра изменения ребра симплекса mu = 1 / N * (log(L0) - log(LN));

gamma = @(k) exp(-mu * k);

%Выделение памяти для массивов

X = zeros(N * (n + 1), n); Y = zeros(N * (n + 1), 1); x = zeros(n + 1, n);

zz = zeros(1, n);

f = zeros(1, n + 1); m = zeros(1, n + 1);

xnov = zeros(n + 1, n); znov = zeros(1, n);

%Количество проведенных экспериментов в исходном симплексе k = n + 1;

%Начальное значение номера вершины, запрещенной к "возврату" p = n + 2;

%Счетчик записи координат вершин симплекса в массив X для построения %графика

142

w = n + 1;
%Построение	исходного		симплекса
for t = 1 :	N1
for j =	1 : n +	1
for	i = 1 : n
	if i < j - 1
	eta	=	0;
	elseif i > j - 1
	eta	=	-1 / sqrt(2 * i * (i + 1));
	else
	eta	=	sqrt(i / (2 * (i + 1)));
	end
	if poz(X0) > 0 || fogr(X0) > 0
	X0(i)		= xmin(i) + rand(1) * (xmax(i) - xmin(i));
	end
	x(j, i)	=	X0(i) + eta * L0;
	X(j, i)	=	x(j, i);
	zz(i) =	x(j, i);

end

end

end

%Проведение эксперимента в вершинах исходного симплекса for j = 1 : n + 1

f(j) = Q(X(j, :));

if poz(X(j, :)) > 0 || fogr(X(j, :)) > 0 f(j) = -900000 - j;

end

Y(j) = f(j); m(j) = 0;

end

%Начало цикла поиска while k <= N

%Поиск худшей вершины ft = sort(f);

fmin = ft(1);

for j = 1 : n + 1 if f(j) == fmin

s = j; break;

end

end

%Запрет "возврата" вершины, отраженной на предыдущем шаге if s == p

fmin = ft(2);

%Поиск второй худшей вершины for j = 1 : n + 1

if f(j) == fmin s = j; break;

end

end

end

p = s;

143

%Расчет координат новой вершины for i = 1 : n

xnov(s, i) = 2 / n * sum(x(:, i)) - (2 + n) / n * x(s, i); znov(i) = xnov(s, i);

end

%Эксперимент в новой вершине f(s) = Q(znov);

%Проверка ограничений задачи

if poz(znov) > 0 || fogr(znov) > 0 f(s) = -900000 - k;

end

%Увеличение счетчика числа шагов поиска, в которых вершина с номером j %не отражалась

for j = 1 : n + 1 m(j) = m(j) + 1;

end

%Обнуление счетчика в новой вершине m(s) = 0;

k = k + 1; for i = 1 : n

x(s, i) = xnov(s, i);

end

%Уменьшение ребра симплекса относительно новой вершины for j = 1 : n + 1

for i = 1 : n if j ~= s

x(j, i) = xnov(s, i) + (x(j, i) - xnov(s, i)) * ...

gamma(k) / gamma(k - 1);

end

end

end

%Формирование массива для построения графиков for j = 1 : n + 1

for i = 1 : n

X(w + j, i) = x(j, i);

end

end

Y(k) = f(s);

w = w + n + 1; mp = max(m);

%Выявление зацикливания симплекса относительно вершины m(j) > n + 2 if mp < n + 2

continue

end

%Проведение повторного эксперимента в вершине с m(j) > n + 2 for j = 1 : n + 1

if m(j) == mp d = j;

end

144

end

for i = 1 : n

zz(i) = x(d, i);

end

f(d) = Q(zz); k = k + 1; Y(k) = f(d); m(d) = 0;

end

%"Отбрасываем" неиспользованную часть памяти для X и Y

X = X(1:w, :); Y = Y(1:k); end

>> lab8

Конечный симплекс:

Координаты: 1.518551 2.541198

Значение целевой функции: 0.129971

Координаты: 1.483946 2.541198

Значение целевой функции: 0.078994

Координаты: 1.501248 2.496516

Значение целевой функции: 0.047460

145

ЛАБОРАТОРНАЯ РАБОТА № 9

КОМПЛЕКС-МЕТОД

Цель работы – исследование метода прямого поиска экстремума целевой функции.

КРАТКОЕ ТЕОРЕТИЧЕСКОЕ ОПИСАНИЕ

Теоретические сведения приведены в разделе 3.9 [1], а также в [3, 4]. Общая постановка задачи

				,	(9.1)
	Q(x) min, max
	x X	kp 1			, j ,...,m .
где X x : x E n , x x x ,i ,...,p, q		j	(x) b
i	i i			j

Алгоритм комплекс-метода.

Исходные данные для решения задачи (9.1): n-размерность вектора х; kp – константа, определяющая поиск максимума (kp=1) или минимума

(kp=0); x , x позиционные и q(x) b функциональные ограничения (9.1); целевая функция Q(x) , N – число шагов поиска; r,e n – номера координат

для построения графика перемещения координат вершин комплекса при поиске; S x , SQ – среднее квадратическое отклонение координат вершин и целе-

вой функции текущего комплекса; δ – шаг перемещения координат вершин Vr , нарушившей позиционные ограничения.

Порядок выполнения операций следующий:

1) сформировать программно координаты x1 первой вершины V1 ,

удовлетворяющей ограничениям задачи. Положить координаты центра тяже-
сти вершин текущего комплекса x0 x1				и вычислить Q x1 .
Построить «комплекс», удовлетворяющий позиционным и функцио-
нальным ограничениям задачи
x	ji	x	p x x , i ,...,n , j 2,...,l ,
		i	i	i

где p rand (1) – равномерно-распределенное случайное число в интервале

[0,1].

Если на j-ой итерации точка (вершина) «комплекса» не удовлетворяет хотя бы одному функциональному ограничению задачи (позиционные при таком формировании выполняются автоматически), то она смещается на половину расстояния до центра тяжести множества уже принятых точек (вер-

шин) «комплекса»			xsi x i , s j, p j, x
x	si			ji	x	si	,i ,...,n ,

146

где xsi – координаты вершины, не удовлетворяющей функциональным ограничениям,

p

x i p j x ji ,i ,...,n ,

где x i – центр тяжести множества уже принятых р вершин «комплекса». Если точка xsi все еще не удовлетворяет функциональному ограниче-

нию, то изложенная процедура «усечения» повторяется до полного выполнения функциональных ограничений задачи;

2) вычисление целевой функции для всех вершин комплекса

Qj Q x j , j ,...,l ;

3) сортировка вершин текущего «комплекса» в порядке возрастания Q j . Выбрать вершину с наибольшим (поиск минимума) Q j , с наименьшим

(поиск максимума) s j , где s – индекс худшей вершины;

4) вычисление координат центра тяжести x текущего «комплекса», за исключением координаты xs

					l
			x i		x ji ,i ,...,n ;
				l j ,
					j s
5)	отразить худшую		вершину		«комплекса» относительно точки x
(рис.9.1)
			xr ( )x xs
где – коэффициент растяжения, рекомендуется брать . ;
		x2			xr
			x0
		xS
					x1
xs	– худшая вершина «комплекса»; x0 – центр тяжести лучших вершин «комплек-
		са»; xr	– новая вершина «комплекса»

Рис. 9.1. Геометрия поиска комплекс-метода

147

6) проверить выполнение позиционных ограничений задачи
	x x	ri	x ,	i 1,...,n .
	i		i
Если нарушена нижняя граница i-ой переменной (рис. 9.2), то
	h i, xrh xri ,				xri xrh ,
где , рекомендуется , ,		если нарушена верхняя граница i-ой пе-
ременной, то
	xrh xri ,				xri xrh ;
	xr
x2
	xrн				x
					2
			x0
				xS
					x1
x0 – центр тяжести лучших вершин «комплекса»; xr – вершина «комплекса», на-
рушившая ограничение; xrн	– скорректированная вершина «комплекса»; xs – отражаемая
вершина «комплекса»; x – верхнее ограничение по переменной x
	2				2
Рис. 9.2. Нарушение позиционного ограничения

7) проверка выполнения функциональных ограничений
q j xr bj , j ,...,l .
Если нарушено хоть одно ограничение (рис. 9.3), то точку xr смещают
на половину расстояния между xr	и центром x
xн	xri x0i ,i 1,...,n .
ri	2

Затем производится повторная проверка точки xнr на выполнение

функциональных ограничений до тех пор, пока не будут выполнены все ограничения. После этого

xri xriн ,i 1,...,n ;

148

	x2	xr
	xн
	r	q(x) b
	x0
		x1
x0	– центр тяжести лучших вершин «комплекса»;	xr – вершина «комплекса», на-
рушившая ограничение; xrн – скорректированная вершина «комплекса»;
	q(x) b – функция ограничения
	Рис. 9.3. Нарушение функционального ограничения

8)	вычисление целевой функции в новой точке
			Qr Q xr ;
9)	если Qr Qs , то шаг признать неудачным, и точку xr следует сме-
стить к центру x на половину расстояния между ними
	xн		xri x i	,i ,...,n
	ri

идалее перейти к п.6), пока не будет получена допустимая точка;

10)если Qr Qs , то шаг признается удачным

xsi xri , ,...,n, Qs Qr ;

11) проверить условие останова. Для этого вычислим оценку среднего квадратического отклонения целевой функции для текущего «комплекса»

l

SQ

Q j

Q ,

l

j

l

где Q

l

Q j , и среднее квадратичное отклонение координат вершин теку-

j

щего «комплекса»

Sx

n

Sxi ,

n i

l

Sxi

x ji xi

, i 1,...,n ,

l j

149

где xi 1 l x ji , i 1,...,n . l j 1

Если SQ ε1 и Sx ε2 , то останов процедуры, иначе переход к п.3). Здесь и – достаточно малые константы, определяющие точность решения задачи.

Координаты экстремума x* определяем по минимальному (максимальному) значению целевой функции в вершинах последнего комплекса Q x* .

Для невыпуклых функций и многоэкстремальных функций рекомендуется повторить процедуру поиска несколько раз. Так как точки «комплекса» формируются случайно, то имеется вероятность поиска глобального экстремума.

ЗАДАНИЕ И ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ

1.Изучите листинг программы и сопоставьте его с описанием алго-

ритма.

2.В соответствии с вариантом выберите условия задачи из лабораторной работы № 8 и внесите соответствующие изменения в программе и проце-

дурах комплекс-метода: n= , векторы xmin= , xmax= размерности n 1, признак поиска максимума (kp=1), минимума (kp=0), function q_ogr и function x_nov (здесь будьте внимательны со знаком неравенства или ).

3.Проведите оптимизацию без функционального ограничения. Для этого выберите такое b1 в function q_ogr, чтобы q(x) X .

4.Повторите п. 3 для α 1.0 и α 1.5, внеся необходимые изменения

впрограмме.

5.Повторите п. 3 с функциональными ограничениями для α 1.3.

6.Решите задачу для противоположного kp.

Выводы по работе должны включать:

	1)	сравнительную				характеристику	эффективности	поиска
Q*	Q*		Q*	2	при разных α (в форме таблицы) и объяснение влияния ко-
	1
эффициента растяжения на скорость поиска (k=							);
	2)	объяснение траектории движения комплекса в условиях ограниче-
ний;

3) сравнительный анализ с результатами лабораторных работ 7 и 8.

КОНТРОЛЬНЫЕ ВОПРОСЫ

1. Что представляет собой комплекс для n 1 и n 3 ?

150