Решение задач оптимизации градиентным методом

clear all; clc;
 
% f(x,y)=x^2+3y^2+4sin(x+5y)
% Значения коэффициентов
b0 = 22;
b1 = 42;
b2 = 54;
b3 = 7;
b4=150;
 
g = 0.04; % постоянная шага
d = 3; % дельта
% Начальная точка
x1 = -9;
x2 =  8;
k = 0.001;  % Счетчик шагов
kmax = 1000; % Предельное число шагов,
% задается для предотвращения зацикливания
% Массивы для хранения промежуточных координат
x1trace = [x1];
x2trace = [x2];
i = 2;
while k < kmax
% Спуск по обеим координатам сразу
% gr1 = b1 + b12*x2 + 2*b11*x1;
% gr2 = b2 + b12*x1 + 2*b22*x2;
gr1 = 2*b0*x1+b1+b2*x2^2+b3*x2;
gr2 = b0*x1^2+b1*x1+2*b2*x2+b3;
x1 = x1 - g*gr1;
x2 = x2 - g*gr2;
% Сохранение координат
x1trace(i) = x1;
x2trace(i) = x2;
i = i + 1;
% Проверка условия останова
sqrt(gr1^2 + gr2^2)
if sqrt(gr1^2 + gr2^2) <= d;
break; % Выход из цикла в случае выполнения условия
end
k = k + 1;
end
% Построение графика
x = -20:0.5:20;
y = -20:0.5:20;
[X, Y] = meshgrid(x, y);
% Z = b0 + b1*X + b2*Y + b12*X.*Y + b11*X.^2 + b22*Y.^2;
Z = b0*X.^2+b1*X+b2*Y.^2+b3*Y+b4;
[C, h] = contour(X, Y, Z);
clabel(C, h);  % Отображение меток на линиях уровня
hold on;
plot(x1trace, x2trace, '-x');
% Вывод начальной точки на график
text(x1trace(1) + 0.2, x2trace(1) + 0.5, 'M0');
% Вывод решения на график
text(x1 + 2, x2, ...
strvcat(['x1 = ' (num2str(x1))], ...
        ['x2 = ' (num2str(x2))], ...
        ['k = '  (num2str(k))]));