global c;
c = 0;
uslov_grad
c
function value = F(x, y)
    value = (conj(x)*(21 - 10i) - conj(y)*(10 + 21i) + (conj(x) - conj(y)*1i)^2*(4 - 18i) + (conj(x) - conj(y)*1i)^3*(10 - 10i) + (17 - 12i))*(x*(21 + 10i) - y*(10 - 21i) + (x + y*1i)^2*(4 + 18i) + (x + y*1i)^3*(10 + 10i) + 17 + 12i);
end

function value = Res_1(x, y)
    value = (x-5.3)^2+(y+0.9)^2-2.7^2; % ограничение1, внутренность окружности R = 2.7, центр в (5.3, -0.9) = (2*R, -R/3)
end
function value = Res_2(x, y) % ограничение2, верхняя часть круга
    value = -y-0.9;
end

function value = M2_F_and_Con(r, x, y)
    value = (conj(x)*(21 - 10i) - conj(y)*(10 + 21i) + (conj(x) - conj(y)*1i)^2*(4 - 18i) + (conj(x) - conj(y)*1i)^3*(10 - 10i) + (17 - 12i))*(x*(21 + 10i) - y*(10 - 21i) + (x + y*1i)^2*(4 + 18i) + (x + y*1i)^3*(10 + 10i) + 17 + 12i) + r*((max(0, Res_1(x,y)))^2 + max(0, Res_2(x,y))^2);
end

function value = M1_F_and_Con(r, x, y)
    value = (conj(x)*(21 - 10i) - conj(y)*(10 + 21i) + (conj(x) - conj(y)*1i)^2*(4 - 18i) + (conj(x) - conj(y)*1i)^3*(10 - 10i) + (17 - 12i))*(x*(21 + 10i) - y*(10 - 21i) + (x + y*1i)^2*(4 + 18i) + (x + y*1i)^3*(10 + 10i) + 17 + 12i) + r*((max(0, Res_1(x,y)))^2 + max(0, Res_2(x,y))^2);
end

function d = M1_Der_x(r, x, y) % производная вспомогательной функции F = f*f_сопр - r*(1/Res_1+1/Res_2) по х
    d = (x*(8 + 36i) - y*(36 - 8i) + (x + y*1i)^2*(30 + 30i) + 21 + 10i)*(conj(x)*(21 - 10i) - conj(y)*(10 + 21i) + (conj(x) - conj(y)*1i)^2*(4 - 18i) + (conj(x) - conj(y)*1i)^3*(10 - 10i) + (17 - 12i)) + (conj(x)*(8 - 36i) - conj(y)*(36 + 8i) + (conj(x) - conj(y)*1i)^2*(30 - 30i) + (21 - 10i))*(x*(21 + 10i) - y*(10 - 21i) + (x + y*1i)^2*(4 + 18i) + (x + y*1i)^3*(10 + 10i) + 17 + 12i) + (r*(2*x - 53/5))/((x - 53/10)^2 + (y + 9/10)^2 - 729/100)^2;
end

function d = M1_Der_y(r, x, y) % производная вспомогательной функции F = f*f_сопр - r*(1/Res_1+1/Res_2) по y
    d =  -(x*(36 - 8i) + y*(8 + 36i) + (x + y*1i)^2*(30 - 30i) + 10 - 21i)*(conj(x)*(21 - 10i) - conj(y)*(10 + 21i) + (conj(x) - conj(y)*1i)^2*(4 - 18i) + (conj(x) - conj(y)*1i)^3*(10 - 10i) + (17 - 12i)) - r*(1/(y + 9/10)^2 - (2*y + 9/5)/((x - 53/10)^2 + (y + 9/10)^2 - 729/100)^2) - (conj(x)*(36 + 8i) + conj(y)*(8 - 36i) + (conj(x) - conj(y)*1i)^2*(30 + 30i) + (10 + 21i))*(x*(21 + 10i) - y*(10 - 21i) + (x + y*1i)^2*(4 + 18i) + (x + y*1i)^3*(10 + 10i) + 17 + 12i);
end

function d = M2_Der_x(r, x, y) % производная вспомогательной функции F = f*f_сопр + r/2*(max(0,Res_1)^2+max(0,Res_2)^2) по х
    d = (x*(8 + 36i) - y*(36 - 8i) + (x + y*1i)^2*(30 + 30i) + 21 + 10i)*(conj(x)*(21 - 10i) - conj(y)*(10 + 21i) + (conj(x) - conj(y)*1i)^2*(4 - 18i) + (conj(x) - conj(y)*1i)^3*(10 - 10i) + (17 - 12i)) + (conj(x)*(8 - 36i) - conj(y)*(36 + 8i) + (conj(x) - conj(y)*1i)^2*(30 - 30i) + (21 - 10i))*(x*(21 + 10i) - y*(10 - 21i) + (x + y*1i)^2*(4 + 18i) + (x + y*1i)^3*(10 + 10i) + 17 + 12i) + r*((2*x - 53/5)*max(0, Res_1(x,y)));
end

function d = M2_Der_y(r, x, y) % производная вспомогательной функции F = f*f_сопр + r/2*(max(0,Res_1)^2+max(0,Res_2)^2) по y
    d = -(x*(36 - 8i) + y*(8 + 36i) + (x + y*1i)^2*(30 - 30i) + 10 - 21i)*(conj(x)*(21 - 10i) - conj(y)*(10 + 21i) + (conj(x) - conj(y)*1i)^2*(4 - 18i) + (conj(x) - conj(y)*1i)^3*(10 - 10i) + (17 - 12i)) - (conj(x)*(36 + 8i) + conj(y)*(8 - 36i) + (conj(x) - conj(y)*1i)^2*(30 + 30i) + (10 + 21i))*(x*(21 + 10i) - y*(10 - 21i) + (x + y*1i)^2*(4 + 18i) + (x + y*1i)^3*(10 + 10i) + 17 + 12i) + r*((2*y + 9/5)  - max(0, Res_2(x,y)));
end

function d = M3_DerF_x(x, y) % производная целевой функции по х
    global c;
    d = (x*(8 + 36i) - y*(36 - 8i) + (x + y*1i)^2*(30 + 30i) + 21 + 10i)*(conj(x)*(21 - 10i) - conj(y)*(10 + 21i) + (conj(x) - conj(y)*1i)^2*(4 - 18i) + (conj(x) - conj(y)*1i)^3*(10 - 10i) + (17 - 12i)) + (conj(x)*(8 - 36i) - conj(y)*(36 + 8i) + (conj(x) - conj(y)*1i)^2*(30 - 30i) + (21 - 10i))*(x*(21 + 10i) - y*(10 - 21i) + (x + y*1i)^2*(4 + 18i) + (x + y*1i)^3*(10 + 10i) + 17 + 12i);
    c = c + 1;
end

function d = M3_DerF_y(x, y) % производная целевой функции по х
    global c;
    d = -(x*(36 - 8i) + y*(8 + 36i) + (x + y*1i)^2*(30 - 30i) + 10 - 21i)*(conj(x)*(21 - 10i) - conj(y)*(10 + 21i) + (conj(x) - conj(y)*1i)^2*(4 - 18i) + (conj(x) - conj(y)*1i)^3*(10 - 10i) + (17 - 12i)) - (conj(x)*(36 + 8i) + conj(y)*(8 - 36i) + (conj(x) - conj(y)*1i)^2*(30 + 30i) + (10 + 21i))*(x*(21 + 10i) - y*(10 - 21i) + (x + y*1i)^2*(4 + 18i) + (x + y*1i)^3*(10 + 10i) + 17 + 12i);
    c = c + 1;
end

function d = M3_DerG1_x(x,y)
    global c;
    c = c + 1;
    d = 2*x - 53/5;
end

function d = M3_DerG1_y(x,y)
    global c;
    c = c + 1;
    d = 2*y + 9/5;
end

function d = M3_DerG2_x(x,y)
    d = 0;
end

function d = M3_DerG2_y(x,y)
    d = -1;
end


function bool = is_in_trust_region(x, y)  % проверка на принадлежность допустимому множеству
    if Res_1(x, y) <= 0 && Res_2(x,y) <= 0
        bool = 1;
    else
        bool = 0;
    end
end

function value = my_F1
    syms x y
    f(x, y) = (10+10*1i)*(x+y*1i).^3+(4+18*1i)*(x+y*1i).^2+(21+10*1i)*(x+y*1i)+(17+12*1i);
    Fun(x, y) = f*conj(f) +r*(1/Res_1(x, y) + 1/Res_2(x, y));
    dx = diff(Fun, x)
    dy = diff(Fun, y)

end
function my_F2
    syms x y
    f(x, y) = (10+10*1i)*(x+y*1i).^3+(4+18*1i)*(x+y*1i).^2+(21+10*1i)*(x+y*1i)+(17+12*1i);
    Fun(x,y) = f*conj(f);
    dRes_1(x, y) = (x-5.3)^2+(y+0.9)^2-2.7^2;
    dRes_2(x, y) = -y-0.9;
    dx = diff(dRes_1, x)
    dy = diff(dRes_1, y)
end
function arg_min = ConstStep(poly, begin_dot)
c = 0;
j = begin_dot; % x_k+1
alfa = 0.0001;
while true
    c = c + 1;
    k = j;     % x_k
    j = j - alfa*Grad(Fun(j, poly), DFun(j, poly)); %x_k+1
    value_j =  Fun(j, poly);
    value_k =  Fun(k, poly);
    if abs(Grad(Fun(j, poly), DFun(j, poly))) < 10^-4
        break;
    end
end
arg_min = j;
end
function min = project1(dot1, dot2, alfa)
    x = dot1;
    y = dot2;
    while true
        grad_fx = M3_DerF_x(x,y);
        grad_fy = M3_DerF_y(x,y);
        grad_gx = M3_DerG1_x(x,y);
        grad_gy = M3_DerG1_y(x,y);
        x = x + alfa*(-grad_fx-(-grad_fx*grad_gx - grad_fy*grad_gy)/norm([grad_gx grad_gy])^2*grad_gx);
        y = y + alfa*(-grad_fy-(-grad_fx*grad_gx - grad_fy*grad_gy)/norm([grad_gx grad_gy])^2*grad_gy);
        if abs(-grad_fx*grad_gx - grad_fy*grad_gy - norm([grad_fx grad_fy])*norm([grad_gx grad_gy])) <= 10^-4 %условие коллинеарности
            break;
        end
    end
    min = [x y];
end

function min = project2(dot1, dot2, alfa)
    x = dot1;
    y = dot2;
    while true
        grad_fx = M3_DerF_x(x,y);
        grad_fy = M3_DerF_y(x,y);
        grad_gx = M3_DerG2_x(x,y);
        grad_gy = M3_DerG2_y(x,y);
        x = x + alfa*(-grad_fx-(-grad_fx*grad_gx - grad_fy*grad_gy)/norm([grad_gx grad_gy])^2*grad_gx);
        y = y + alfa*(-grad_fy-(-grad_fx*grad_gx - grad_fy*grad_gy)/norm([grad_gx grad_gy])^2*grad_gy);
        if abs(-grad_fx*grad_gx - grad_fy*grad_gy - norm([grad_fx grad_fy])*norm([grad_gx grad_gy])) <= 10^-4 %условие коллинеарности
            break;
        end
    end
    min = [x y];
end

function min = shtraf %МЕТОД ВНУТРЕННИХ В АББАСОВЕ (Метод штрафных функций machinelearning)
x = 6; y = 1.5;
r = 10;
alfa = 10^-4;
c = 1;
xtrace = [x];
ytrace = [y];
while true
    while true
        c = c + 1;
        dx = M1_Der_x(r, x, y);
        dy = M1_Der_y(r, x, y);
        x_t = x; y_t = y;
        x = x - alfa*dx;
        y = y - alfa*dy;
        xtrace(c) = x;
        ytrace(c) = y;
        if is_in_trust_region(x, y) == 0
            x = x_t; y = y_t;
            alfa = alfa / 5;
        end
        if abs(x-x_t) <= 10^-2 && abs(y-y_t) <= 10^-2
            break;
        end
    end
if abs(-r*(1/Res_1(x,y) + 1/Res_2(x,y))) <= 10^-3
    break
end
    r = r / 5;
end
grad_fx = M3_DerF_x(x,y);
grad_fy = M3_DerF_y(x,y);
grad_gx = M3_DerG1_x(x,y);
grad_gy = M3_DerG1_y(x,y);
xt = x; yt = y;
x = x - grad_fx*Res_1(xt, yt)/(grad_fx*grad_gx + grad_fy*grad_gy);
y = y - grad_fy*Res_1(xt, yt)/(grad_fx*grad_gx + grad_fy*grad_gy);

% x = -5:0.1:6;
% y = -5:0.1:6;
% [X, Y] = meshgrid(x, y);
% Z = (conj(X)*(21 - 10i) - conj(Y)*(10 + 21i) + (conj(X) - conj(Y)*1i)^2*(4 - 18i) + (conj(X) - conj(Y)*1i)^3*(10 - 10i) + (17 - 12i))*(X*(21 + 10i) - Y*(10 - 21i) + (X + Y*1i)^2*(4 + 18i) + (X + Y*1i)^3*(10 + 10i) + 17 + 12i);
% [C, h] = contour(X, Y, real(Z));
% clabel(C, h);  % Отображение меток на линиях уровня
% hold on;
% plot(xtrace, ytrace, '-x');
% %Вывод начальной точки на график
% text(xtrace(1) + 0.2, ytrace(1) + 0.5, 'M0');
% %Вывод решения на график
% text(x + 2, y, ...
% strvcat(['x = ' (num2str(x))], ...
%         ['y = ' (num2str(y))], ...
%         ['c = '  (num2str(c))]));

min = [x y];
fval = F(x, y)
end

function min = vnesh_shtraf
x = 0; y = 0;
r = 10^-1;
alfa = 10^-4;
xtrace(1) = [x]; ytrace(1) = [y];
c = 1;
while true
    while true
        c = c + 1;
        dx = alfa*M2_DerFun_x(r, x, y);
        dy = alfa*M2_DerFun_y(r, x, y);
        norm = sqrt(dx^2+dy^2);
        step_x = dx / norm;
        step_y = dy / norm;
%         if abs(dy) > 10
%             dy = dy / norm;
%         end
%         if abs(dx) > 10
%             dx = dx / norm;
%         end
        x_t = x; y_t = y;
        x = x - step_x;
        y = y - step_y;
        xtrace(c) = x;
        ytrace(c) = y;
%         if is_in_trust_region(x, y) == 0
%             x = x_t; y = y_t;
%             alfa = alfa / 2;
%         end
        if sqrt(dx^2+dy^2) <= 10^-3
            break;
        else
            alfa = alfa / 2;
        end
    end
if r/2*(max(0,Res_1(x,y))^2+max(0,Res_2(x,y))^2) <= 10^-4
    break
end
r = r * 10;
end
bool = is_in_trust_region(x, y)
% x = -5:0.1:6;
% y = -5:0.1:6;
% [X, Y] = meshgrid(x, y);
% Z = (conj(X)*(21 - 10i) - conj(Y)*(10 + 21i) + (conj(X) - conj(Y)*1i)^2*(4 - 18i) + (conj(X) - conj(Y)*1i)^3*(10 - 10i) + (17 - 12i))*(X*(21 + 10i) - Y*(10 - 21i) + (X + Y*1i)^2*(4 + 18i) + (X + Y*1i)^3*(10 + 10i) + 17 + 12i);
% [C, h] = contour(X, Y, real(Z));
% clabel(C, h);  % Отображение меток на линиях уровня
% hold on;
% plot(xtrace, ytrace, '-x');
% %Вывод начальной точки на график
% text(xtrace(1), ytrace(1), 'M0');
% %Вывод решения на график
% text(x + 2, y, ...
% strvcat(['x = ' (num2str(x))], ...
%         ['y = ' (num2str(y))], ...
%         ['c = '  (num2str(c))]));
min = [x y];
fval = F(x, y)
gradx = M2_DerFun_x(r, x, y)
grady = M2_DerFun_y(r, x, y)
end

function min = uslov_grad
    x = 6; y = 1.5;
    alfa = 10^-4;
    xtrace = [x]; ytrace = [y];
    c = 1;
    while true
        c = c + 1;
        dx = M3_DerF_x(x, y);
        dy = M3_DerF_y(x, y);
        x_t = x; y_t = y; xt = x; yt = y;
        x = x - alfa*dx;
        xtrace(c) = x;
        y = y - alfa*dy;
        ytrace(c) = y;
        if is_in_trust_region(x, y) == 0
            x = x_t; y = y_t;
            xt = xt + 1; yt = yt + 1;
            alfa = alfa / 2;
        end
        if abs(x-xt) <= 10^-2 && abs(y-yt) <= 10^-2
            break;
        end
    end
    if abs(Res_1(x,y)) <= 10^-1 && abs(Res_2(x,y)) <= 10^-1
        min = [2.6 -0.9];
    elseif abs(Res_1(x,y)) <= 10^-1
        min = project1(x, y, alfa);
    else
        min = project2(x, y, alfa);
    end
%     x = -5:0.1:6;
%     y = -5:0.1:6;
%     [X, Y] = meshgrid(x, y);
%     Z = (conj(X)*(21 - 10i) - conj(Y)*(10 + 21i) + (conj(X) - conj(Y)*1i)^2*(4 - 18i) + (conj(X) - conj(Y)*1i)^3*(10 - 10i) + (17 - 12i))*(X*(21 + 10i) - Y*(10 - 21i) + (X + Y*1i)^2*(4 + 18i) + (X + Y*1i)^3*(10 + 10i) + 17 + 12i);
%     [C, h] = contour(X, Y, real(Z));
%     clabel(C, h);  % Отображение меток на линиях уровня
%     hold on;
%     plot(xtrace, ytrace, '-x');
%     %Вывод начальной точки на график
%     text(xtrace(1) + 0.2, ytrace(1) + 0.5, 'M0');
%     %Вывод решения на график
%     text(x + 2, y, ...
%     strvcat(['x = ' (num2str(x))], ...
%             ['y = ' (num2str(y))], ...
%             ['c = '  (num2str(c))]));
    fval = F(min(1), min(2))
    bool = is_in_trust_region(x, y);
end