Глава 3. Двухмерная оптимизация с применением пакета MATLAB

Метод сканирования

В области исследуемой функции двух переменных выбирается значение одной из переменных, соответствующее краю области поиска. По другой переменной интервал поиска разбивается на равные участки, длина каждого из которых равна шагу поиска. Далее значение функции последовательно определяется во всех точках разбиения и на краях отрезка, наименьшее из них запоминается. Затем увеличивают переменную, значение которой было постоянным, на шаг поиска и повторют процедуру. И так до тех пор, пока не будет исследована вся область поиска с заранее выбранным шагом по каждой из переменных. Минимальное из всех сохраненных значений минимумов и есть глобальный минимум.

X₁

40 30 20 40 30

X₂

Рис. 3.1. Графическое представление поиска экстремума для случая двух переменных методом сканирования

Далее приведен текст программы, позволяющий реализовать описанный выше алгоритм.

function[]=Scanirovanie2D_050809();

function z=fz_xy2(x, y)

nx=length(x); %функция вычисляется от двух аргументов и возвращает одно значение

ny=length(y); %или функция вычисляется от векторов аргументов

% и возвращает вектор, каждый элемент которого

% вычислен от соответствующих элементов векторов аргументов

for i=1: nx

for j=1: ny

z(i, j)=-sqrt(256-x(i)*x(i)-y(j)*y(j));

end

%метод сканирования 2-мерный

function[x, y, z, masx, masy, masz, xleft, xright, yleft, yright, hx]=scan2();

disp('МЕТОД ПОИСКА МИНИМУМА ФУНКЦИЙ ДВУХ ПЕРЕМЕННЫХ МЕТОДОМ СКАНИРОВАНИЯ');

xleft=input('введите левую границу диапазона X поиска! ');

xright=input('введите правую границу диапазона X поиска! ');

hx=input('введите допустимую погрешность X');

yleft=input('введите левую границу диапазона Y поиска! ');

yright=input('введите правую границу диапазона Y поиска! ');

hy=input('введите допустимую погрешность Y');

%создание массива точек

masx=xleft: hx: xright;

masy=yleft: hy: yright;

masz=fz_xy2(masx, masy);

%поиск минимума

min=masz(1, 1); % стандартный алгоритм поиска минимума

numx=1;

numy=1;

nx=length(masx);

ny=length(masy);

for i=1: nx

for j=1: ny

if masz(i, j)< min

min=masz(i, j);

numx=i;

numy=j;

end

%возвращаемые значения (минимум)

x=masx(numx);

y=masy(numy);

z=min;

end

%обращение к функции сканирования

[x, y, z, masx, masy, masz, xleft, xright, yleft, yright, hx]=scan2();

%выбор оформления вывода на экран

choiceTab=input('ВОПРОС выводить пошагово как таблицу? Да=1 нет=0 ');

choiceGraf=input('ВОПРОС выводить трехмерный график? Да=1 нет=0');

%вывод таблицы

if choiceTab==1

disp('вывод результатов');

for i=1: length(masx)

for j=1: length(masy)

a=(length(masy))*(i-1)+j;

string=print(‘номер= %4i\tаргумент x= %7.3f\tаргумент y= %7.3f\tфункция= %7.3f’, a, masx(i), masy(j), masz(i, j));

disp(string);

end

% подготовка к выводу текста

% num2str преобразует число в строку символов

disp('НАЙДЕН МИНИМУМ');

sxy=strcat('оптимальные значения аргументов: х=', num2str(x), ' у=', num2str(y));

sz=strcat('минимум функции=', num2str(z));

disp(sxy) % вывод на экран

disp(sz)

%подготовка данных для графика

if choiceGraf==1

n=floor(abs((xright-xleft)/hx));

hy=(yright-yleft)/n;

% подготовка к построению графика

masx=xleft: hx: xright;

masy=yleft: hy: yright;

masz=fz_xy2(masx, masy);

mesh(masx, masy, masz); %построена сеточная поверхность

grid on;

title(‘z=-sqrt(256-x.^2-y.^2)’);

xlabel(‘X’);

ylabel(‘Y’);

zlabel(‘Z’);

text(x, y, z, ’\leftarrow Minimum’);

zeroMas=masx*0;

hold on;

%surf(zeroMas, zeroMas, masz); % построены оси

%surf(zeroMas, masy, zeroMas); % построены оси

%surf(masx, zeroMas, zeroMas); % построены оси

legend(‘z(i, j)=-sqrt(256-x(i)*x(i)-y(j)*y(j))’, 0); % подписи к линиям на графике

end

После выполнения этой программы пользователь видит на экране соответствующий текст, вводит нужные значения и получает график:

МЕТОД ПОИСКА МИНИМУМА ФУНКЦИЙ ДВУХ ПЕРЕМЕННЫХ МЕТОДОМ СКАНИРОВАНИЯ

введите левую границу диапазона X поиска! -10

введите правую границу диапазона X поиска! 10

введите допустимую погрешность X.5

введите левую границу диапазона Y поиска! -10

введите правую границу диапазона Y поиска! 10

введите допустимую погрешность Y.5

ВОПРОС выводить пошагово как таблицу? Да=1 нет=0 1

ВОПРОС выводить трехмерный график? Да=1 нет=0 1

вывод результатов

номер= 1 аргумент x=-10.000 аргумент y=-10.000 функция= -7.483

номер= 2 аргумент x=-10.000 аргумент y= -9.500 функция= -8.109

номер= 3 аргумент x=-10.000 аргумент y= -9.000 функция= -8.660

номер= 4 аргумент x=-10.000 аргумент y= -8.500 функция= -9.152

номер= 5 аргумент x=-10.000 аргумент y= -8.000 функция= -9.592

номер= 6 аргумент x=-10.000 аргумент y= -7.500 функция= -9.987

Часть вывода вырезана

номер=1677 аргумент x= 10.000 аргумент y= 8.000 функция= -9.592

номер=1678 аргумент x= 10.000 аргумент y= 8.500 функция= -9.152

номер=1679 аргумент x= 10.000 аргумент y= 9.000 функция= -8.660

номер=1680 аргумент x= 10.000 аргумент y= 9.500 функция= -8.109

номер=1681 аргумент x= 10.000 аргумент y= 10.000 функция= -7.483

НАЙДЕН МИНИМУМ

оптимальные значения аргументов: х=0 у=0

минимум функции=-16

Рис. 3.2. Пример применения метода сканирования

Метод Гаусса–Зейделя

Движение к оптимуму происходит поочередно по каждой из осей до нахождения частного экстремума (рис. 3.3). После нахождения частного экстремума по одной из осей начинается поиск экстремума по другой оси до частного минимума при условии, что значение первой переменной равно найденному минимуму (максимуму) на первой оси. И так поочередно. Процесс последовательно продолжается до тех пор, пока не будет достигнута заданная точность локализации экстремума, т. е. если шаг по каждой из осей приводит к возрастанию (уменьшению) функции, а величина шага меньше или равна заданной точности поиска, то расчет закончен.

Рис. 3.3. Путь поиска экстремума методом поочерёдного изменения переменных

Далее приведен текст программы, позволяющий реализовать описанный выше алгоритм.

Function[]=GaussZeidel2D_170809();

function z=fz_xy2(x, y);

nx=length(x);

ny=length(y);

for i=1: nx

for j=1: ny

z(i, j)=-sqrt(256-(x(i))^2-(y(j))^2);

end

disp('ДВУМЕРНЫЙ ПОИСК МИНИМУМА ФУНКЦИИ ');

disp('С ИСПОЛЬЗОВАНИЕМ МЕТОДА ГАУССА–ЗЕЙДЕЛЯ');

disp('применен для указанной функции: z(i, j)=-sqrt(256-(x(i))^2-(y(j))^2);

disp('ВВОД ИСХОДНЫХ ДАННЫХ точность по аргументу, края диапазона поиска');

function[x, z, masx, masz]=scanx1(xleft, xright, h, y);

%создание массива точек

n=ceil(abs((xright-xleft)/h)); % определяем число шагов

% функция ceil округляет число в сторону плюс бесконечности

masx(1)=xleft;

masz(1)=fz_xy2(xleft, y);

for i=2: n

masx(i)=masx(i-1)+h;

masz(i)=fz_xy2(masx(i), y);

end

% поиск минимума

min=masz(1); % стандартный алгоритм поиска минимума

num=1;

for i=2: n

if masz(i)< min

min=masz(i);

num=i;

end

%возвращаемые значения (минимум)

x=masx(num);

z=masz(num);

end

function[y, z, masy, masz]=scany1(yleft, yright, hy, x);

%создание массива точек

n=ceil(abs((yright-yleft)/hy));

masy(1)=yleft;

masz(1)=fz_xy2(x, yleft);

for i=2: n

masy(i)=masy(i-1)+hy;

masz(i)=fz_xy2(x, masy(i));

end

% поиск минимума

min=masz(1); % стандартный алгоритм поиска минимума

num=1;

for i=2: n

if masz(i)< min

min=masz(i);

num=i;

end

%возвращаемые значения (минимум)

y=masy(num);

z=masz(num);

end

%метод сканирования 2-мерный

function[x, y, z, xleft, xright, yleft, yright, xn]=scan2();

xleft=input('введите левую границу диапазона X поиска=');

xright=input('введите правую границу диапазона X поиска=');

hx=input('введите допустимую погрешность X=');

yleft=input('введите левую границу диапазона Y поиска=');

yright=input('введите правую границу диапазона Y поиска=');

hy=input('введите допустимую погрешность Y=');

eps=input('введите минимальное изменение функции (можно ноль): ');

choiceTab=input('ВОПРОС выводить пошагово как таблицу? Да=1 нет=0 ');

%создание массива точек

xn=ceil(abs((xright-xleft)/hx)); % определяем число шагов

%объявляем переменные для возвращаемых значений

x=xleft;

y=yleft;

z=-1;

disp('ВЫВОД ПРОМЕЖУТОЧНЫХ РЕЗУЛЬТАТОВ ПОШАГОВО ');

for i=1: xn

%обращение к функции сканирования

[x, z]=scanx1(xleft, xright, hx, y); % по х

%вывод таблицы

if choiceTab==1

string=print(‘номер= %4i\tаргумент x= %7.3f\tаргумент y= %7.3f\tфункция= %7.3f\tпоиск по х’, i, x, y, z);

disp(string);

end

zx=z;

[y, z]=scany1(yleft, yright, hy, x); % по у

%вывод таблицы

if choiceTab==1

string=print(‘номер= %4i\tаргумент x= %7.3f\tаргумент y= %7.3f\tфункция= %7.3f\tпоиск по y’, i, x, y, z);

disp(string);

end

zy=z;

dz(i)=abs(zx-zy);

if dz(i)< eps

disp(‘dz< eps’);

break;

end

%обращение к функции сканирования

[x, y, z, xleft, xright, yleft, yright, xn]=scan2();

%выбор оформления вывода на экран

choiceGraf=input('ВОПРОС выводить трехмерный график? Да=1 нет=0 ');

% подготовка к выводу текста

% num2str преобразует число в строку символов

disp('НАЙДЕН МИНИМУМ');

sxy=strcat('оптимальные значения аргументов: х=', num2str(x), ' у=', num2str(y));

sz=strcat('минимум функции=', num2str(z));

disp(sxy) % вывод на экран

disp(sz)

%подготовка данных для графика

if choiceGraf==1

hx=(xright-xleft)/xn;

hy=(yright-yleft)/xn;

% подготовка к построению графика

masx=xleft: hx: xright;

masy=yleft: hy: yright;

masz=fz_xy2(masx, masy);

surf(masx, masy, masz);

grid on;

title(‘z=-sqrt(256-x.^2-y.^2)’);

xlabel(‘X’);

ylabel(‘Y’);

zlabel(‘Z’);

smin=strcat(‘\leftarrow Minimum (‘, num2str(x), ’, ’, num2str(y), ’, ’, num2str(z), ’)’);

text(x, y, z, smin);

zeroMas=masx*0;

hold on;

%surf(zeroMas, zeroMas, masz); % построены оси

%surf(zeroMas, masy, zeroMas); % построены оси

%surf(masx, zeroMas, zeroMas); % построены оси

legend(‘z(i, j)=-sqrt(256-x(i)*x(i)-y(j)*y(j))’, 0);

end

ДВУМЕРНЫЙ ПОИСК МИНИМУМА ФУНКЦИИ

С ИСПОЛЬЗОВАНИЕМ МЕТОДА ГАУССА–ЗЕЙДЕЛЯ

Применен для указанной функции: z(i, j)=-sqrt(256-(x(i))^2-(y(j))^2)

ВВОД ИСХОДНЫХ ДАННЫХ точность по аргументу, края диапазона поиска

введите левую границу диапазона X поиска=-10

введите правую границу диапазона X поиска=10

введите допустимую погрешность X=.5

введите левую границу диапазона Y поиска=-10

введите правую границу диапазона Y поиска=10

введите допустимую погрешность Y=.5

введите минимальное изменение функции (можно ноль):.1

ВОПРОС выводить пошагово как таблицу? Да=1 нет=0 1

ВЫВОД ПРОМЕЖУТОЧНЫХ РЕЗУЛЬТАТОВ ПОШАГОВО

номер= 1 аргумент x= 0.000 аргумент y=-10.000 функция=-12.490 поиск по х

номер= 1 аргумент x= 0.000 аргумент y= 0.000 функция=-16.000 поиск по y

номер= 2 аргумент x= 0.000 аргумент y= 0.000 функция=-16.000 поиск по х

номер= 2 аргумент x= 0.000 аргумент y= 0.000 функция=-16.000 поиск по y

dz< eps

ВОПРОС выводить трехмерный график? Да=1 нет=0 1

НАЙДЕН МИНИМУМ

оптимальные значения аргументов: х=0 у=0

минимум функции=-16

Рис. 3.4. Пример применения метода Гаусса–Зейделя

Метод пробных движений

По каждой из переменных:

· из исходной точки делается маленький пробный шаг;

· находится значение функции;

· шаг делается обратно, чтобы вернуться в исходную точку.

Из всех опробованных направлений выбирается то, в котором уменьшение (увеличение) целевой функции наибольшее. Делается шаг (возможно, больше пробного) в выбранном таким образом направлении. Процедура повторяется до достижения заданной точности поиска.

Далее приведен текст программы, позволяющий реализовать описанный выше алгоритм.

function[]=ProbShag2D_170809();

function [f]=Function(xVector);

nVar=length(xVector); % nVar означает число варьируемых переменных

f=1;

for i=1: nVar

f=f+(xVector(i)-2)^2;

end

%метод пробных движений

disp('ДВУМЕРНЫЙ ПОИСК МИНИМУМА ФУНКЦИИ ');

disp('С ИСПОЛЬЗОВАНИЕМ МЕТОДА ПРОБНОГО ШАГА');

disp('применен для указанной функции: z(i, j)=-sqrt(256-(x(i))^2-(y(j))^2)');

disp('ВВОД ИСХОДНЫХ ДАННЫХ точность по аргументу, края диапазона поиска');

function[xOptVector, fMinValue, leftBorderVector, rightBorderVector, stepNumberVector]=FunctionPS();

numberOfVariables=input('введите число переменных: ');

%ввод вектора аргументов, их левых и правых границ, их шагов

for i=1: numberOfVariables

lBVstring=strcat('введите левую границу диапазона поиска по х(', num2str(i), '): ');

leftBorderVector(i)=input(lBVstring);

rBVstring=strcat('введите правую границу диапазона поиска по х(', num2str(i), '): ');

rightBorderVector(i)=input(rBVstring);

sNVstring=strcat(‘введите величину шага по х(‘, num2str(i), ’): ’);

stepNumberVector(i)=input(sNVstring);

OK=input(‘перейти к следующей переменной=1 переделать ввод текущей переменной=0’);

if OK==0

i=i-1;

continue;

end

k=1;

for i=1: numberOfVariables

xVector(i)=(leftBorderVector(i)+rightBorderVector(i))/2; %все переменные по середине диапазонов поиска

end

funVector(k)=Function(xVector); % вычислили функцию

df=-1;

%основной цикл поиска минимума функции

control=1;

while (df< 0)& (control> 0)

if k> 1000

disp('k> 1000');

break;

end

%делаем пробные шаги в сторону возрастания переменных по всем

%переменным и создаем первую половину вектора изменений функции

l=0;

for i=1: (numberOfVariables*2)

funChangeVector(i)=0; % массив, хранящих изменения значений функций при следанных шагах

end

disp('делаем пробные шаги в сторону возрастания переменных');

for i=1: numberOfVariables

if (leftBorderVector< =xVector(i))& (xVector(i)< =rightBorderVector(i))

xVector(i)=xVector(i)+stepNumberVector(i);

l=l+1;

funChangeVector(l)=Function(xVector)-funVector(k);

disp(strcat(‘funChangeVector(‘, num2str(l), ’)=’, num2str(funChangeVector(l))));

xVector(l)=xVector(l)-stepNumberVector(l);

else

control=-1;

disp(‘control=-1; ’);

end

%делаем пробные шаги в сторону убывания переменных по всем переменным и создаем вектор %изменений функции

l=0;

disp('делаем пробные шаги в сторону убывания переменных');

for i=1: numberOfVariables

if (leftBorderVector< =xVector(i))& (xVector(i)< =rightBorderVector(i))

xVector(i)=xVector(i)-stepNumberVector(i);

l=l+1;

funChangeVector(l+numberOfVariables)=Function(xVector)-funVector(k);

disp(strcat(‘funChangeVector(‘, num2str(l), ’)=’, num2str(funChangeVector(l+numberOfVariables))));

xVector(l)=xVector(l)+stepNumberVector(l);

else

control=-1;

end

disp(strcat(‘xVector после шага ‘, num2str(k), ’ равен ‘, num2str(xVector), ’ funChangeVector после шага ‘, num2str(k), ’ равен ‘, num2str(funChangeVector)));

minCF=0;

for i=1: (numberOfVariables*2)

if minCF> funChangeVector(i)

minCF=funChangeVector(i);

if i< =numberOfVariables

indexOfArg=i;

znak=1;

else

indexOfArg=i-numberOfVariables;

znak=-1;

end

disp(strcat(‘Наибольшее убывание функциии по переменной номер ‘, num2str(indexOfArg), ’ эта переменная =’, num2str(xVector(indexOfArg))));

%нашли переменную, по которой наибольшее убывание функции

%увеличиваем (или уменьшаем) ее на полную величину шага

if znak==1

xVector(indexOfArg)=xVector(indexOfArg)+stepNumberVector(indexOfArg);

else

xVector(indexOfArg)=xVector(indexOfArg)-stepNumberVector(indexOfArg);

end

k=k+1;

funVector(k)=Function(xVector); % вычислили функцию

df=funVector(k)-funVector(k-1); % вычислили ее приращение

if df> 0 %если пройден минимум и функция начала возрастать

%отступаем на один шаг назад

if znak==1

xVector(indexOfArg)=xVector(indexOfArg)-stepNumberVector(indexOfArg);

else

xVector(indexOfArg)=xVector(indexOfArg)+stepNumberVector(indexOfArg);

end

funVector(k)=Function(xVector); % вычислили функцию

disp(‘df> 0’);

end

disp(strcat(‘xVector: ‘, num2str(xVector), ’ function= ‘, num2str(funVector(k))));

end

%возвращаемые значения

xOptVector=xVector;

fMinValue=funVector(k);

end %конец функции поиска методом пробных движений

%обращение к функции

[xOptVector, fMinValue, leftBorderVector, rightBorderVector, stepNumberVector] =FunctionPS();

%вывод ответа в виде текста

disp('НАЙДЕН МИНИМУМ ФУНКЦИИ');

for i=1: length(xOptVector)

string=sprint(‘номер переменной= %4i\t переменная= %7.3f\t’, i, xOptVector(i));

disp(string);

end

string=sprint(‘минимум функции= %7.3f’, fMinValue);

disp(string);

% вывод на экран

%подготовка данных для графика

if length(xOptVector)==5

hx=(rightBorderVector(1)-leftBorderVector(1))/100;

hy=(rightBorderVector(2)-leftBorderVector(2))/100;

x=leftBorderVector(1): hx: rightBorderVector(1);

y=leftBorderVector(2): hy: rightBorderVector(2);

z(2, 2)=0;

for i=2: 100

for j=2: 100

x(i)=x(i-1)+hx;

y(j)=y(j-1)+hy;

xV(1)=x(i);

xV(2)=y(j);

z(i, j)=Function(xV);

end

surf(x, y, z);

title(‘f=f+(xVector(i)-2)^2; ’);

xlabel(‘X’);

ylabel(‘Y’);

zlabel(‘Z’);

text(xOptVector(1), xOptVector(2), fMinValue, ’\leftarrow Minimum ‘);

legend(‘f=f+(xVector(i)-2)^2’, 0);

end

В данном случае рисункок к данному примеру поиска экстремума функции не приведен из-за его сходства с приведенными выше (например, с рис. 3.4).

После выполнения этой программы пользователь видит на экране соответствующий текст, вводит нужные значения и получает результат:

ДВУМЕРНЫЙ ПОИСК МИНИМУМА ФУНКЦИИ

С ИСПОЛЬЗОВАНИЕМ МЕТОДА ПРОБНОГО ШАГА

Применен для указанной функции: z(i, j)=-sqrt(256-(x(i))^2-(y(j))^2)

ВВОД ИСХОДНЫХ ДАННЫХ точность по аргументу, края диапазона поиска

введите число переменных: 2

введите левую границу диапазона поиска по х(1): -10

введите правую границу диапазона поиска по х(1): 10

введите величину шага по х(1): 0.5

перейти к следующей переменной=1 переделать ввод текущей переменной=01

введите левую границу диапазона поиска по х(2): -10

введите правую границу диапазона поиска по х(2): 10

введите величину шага по х(2): 0.5

перейти к следующей переменной=1 переделать ввод текущей переменной=01

делаем пробные шаги в сторону возрастания переменных

изменение функции после шага по: 1 переменной =-1.75

изменение функции после шага по: 2 переменной =-1.75

делаем пробные шаги в сторону убывания переменных

изменение функции после шага по: 1 переменной =2.25

изменение функции после шага по: 2 переменной =2.25

вектор переменных после шага: 1 равен: 0 0

Наибольшее убывание функциии по переменной номер: 1 эта переменная =0

вектор переменных: 0.5 0 function=7.25

делаем пробные шаги в сторону возрастания переменных

изменение функции после шага по: 1 переменной =-1.25

изменение функции после шага по: 2 переменной =-1.75

делаем пробные шаги в сторону убывания переменных

изменение функции после шага по: 1 переменной =1.75

изменение функции после шага по: 2 переменной =2.25

вектор переменных после шага: 2 равен: 0.5 0

Наибольшее убывание функциии по переменной номер: 2 эта переменная =0

вектор переменных: 0.5 0.5 function=5.5

Часть вывода вырезана

делаем пробные шаги в сторону возрастания переменных

изменение функции после шага по: 1 переменной =-0.75

изменение функции после шага по: 2 переменной =-0.75

делаем пробные шаги в сторону убывания переменных

изменение функции после шага по: 1 переменной =1.25

изменение функции после шага по: 2 переменной =1.25

вектор переменных после шага: 5 равен: 1 1

Наибольшее убывание функциии по переменной номер: 1 эта переменная =1

вектор переменных: 1.5 1 function=2.25

делаем пробные шаги в сторону возрастания переменных

изменение функции после шага по: 1 переменной =-0.25

изменение функции после шага по: 2 переменной =-0.75

делаем пробные шаги в сторону убывания переменных

изменение функции после шага по: 1 переменной =0.75

изменение функции после шага по: 2 переменной =1.25

вектор переменных после шага: 6 равен: 1.5 1

Наибольшее убывание функциии по переменной номер: 2 эта переменная =1

вектор переменных: 1.5 1.5 function=1.5

Часть вывода вырезана

делаем пробные шаги в сторону возрастания переменных

изменение функции после шага по: 1 переменной =0.25

изменение функции после шага по: 2 переменной =-0.25

делаем пробные шаги в сторону убывания переменных

изменение функции после шага по: 1 переменной =0.25

изменение функции после шага по: 2 переменной =0.75

вектор переменных после шага: 8 равен: 2 1.5

Наибольшее убывание функциии по переменной номер: 2 эта переменная =1.5

вектор переменных: 2 2 function=1

делаем пробные шаги в сторону возрастания переменных

изменение функции после шага по: 1 переменной =0.25

изменение функции после шага по: 2 переменной =0.25

делаем пробные шаги в сторону убывания переменных

изменение функции после шага по: 1 переменной =0.25

изменение функции после шага по: 2 переменной =0.25

⇐ Предыдущая 1 2 3 4 5 67