Zoutendijk's Method

clc;
clear;
 
% Intitial x and y values
% Test with (3;4), (1;-.5), (-.6,-.5)
x_val = [-.6;-.5];
x = x_val(1,1);
y = x_val(2,1);
i = 1;
alpha = .5;
 
% Evaluate constraint values at x and y 
g1 = con_g1(x,y)
g2 = con_g2(x,y)
g3 = con_g3(x,y)
g4 = con_g4(x,y)
 
% set initial search direction as gradient
s(:,1) = [diff_x(x,y),diff_y(x,y)];
 
% Start iterations and stop after 50
for k = 1:50
% Check which constraints are violated and set direction
    % Reset x and y values
    if k > 1
        x = x_val(1,i);
        y = x_val(2,i);
    end
    
    % All constraints satisfied, search direction is gradient
    if (g1 < 0) && (g2 < 0) && (g3 < 0) && (g4 < 0)
        disp("Feasible region")
        f_x = diff_x(x,y);
        f_y = diff_y(x,y);
        s(:,i) = [f_x,f_y];
    % Point violates the second constraint 
    % First constraint not considered since the second is more restrictive
    elseif g2 >= 0
        if g2 == 0
            disp("On constraint 2")
            s(:,i) = diff_g2(x,y);
        elseif g2 > 0
            disp("Violates constraint 2")
            s(:,i) = [1,5/3];
        end 
    % Point violates the third constraint
    elseif g3 >= 0 
        if g3 == 0
            disp("On constraint 3")
            s(:,i) = diff_g3(x,y);
        elseif g3 > 0
            disp("Violates constraint 3")
            s(:,i) = [-1,0];
        end
    % Point violates the fourth constraint   
    elseif g4 >= 0 & g3 > 0 
        if g4 == 0
            disp("On constraint 4")
            s(:,i) = diff_g4(x,y);
        elseif g4 > 0
            disp("Violates constraint 4")
            s(:,i) = [0,-1];
        end
    else 
        f_x = diff_x(x,y);
        f_y = diff_y(x,y);
        s(:,i) = [f_x,f_y];
    end
 
    i = i+1;
    
  %  Update x and y values and evaluate function and constraints
        x_val(:,i)=x_val(:,i-1)+alpha*-s(:,i-1);
        f_new(i) = myfunc(x_val(1,i),x_val(2,i));
        g1 = con_g1(x_val(1,i),x_val(2,i));
        g2 = con_g2(x_val(1,i),x_val(2,i));
        g3 = con_g3(x_val(1,i),x_val(2,i));
        g4 = con_g4(x_val(1,i),x_val(2,i));
        % Choose smaller alpha is within feasible region and improving
        if f_new(i) <= myfunc(x,y) && (g1 <= 0 && g2<0 && g3 <= 0 && g4<0)
            disp("Improve within feasible region")
            alpha = .001;
        end
        alpha = alpha*.9;
end
x_val
 
% Original function to minimize
function func=myfunc(x,y)
    func = x^2+y^2-6*x-8*y+10;
end
 
% First constraint g1
function con1 = con_g1(x,y)
    con1 = 4*x^2+y^2-16;
end
 
% Second constraint g2
function con2 = con_g2(x,y)
    con2 = 3*x+5*y-4; 
end
 
% X is postive constraint
function con3 = con_g3(x,y)
    con3 = -x;
end
 
% Y is positive constraint
function con4 = con_g4(x,y)
    con4 = -y;
end
 
% Gradient evaluations:
function grad_x=diff_x(x,y)
    grad_x = 2*x - 6 
end
 
function grad_y=diff_y(x,y)
    grad_y = 2*y-8
end
 
function grad_g1 = diff_g1(x,y)
    grad_g1 = [8*x,2*y];
end
 
function grad_g2 = diff_g12(x,y)
    grad_g1 = [3,5];
end
 
function grad_g3 = diff_g3(x,y)
    grad_g1 = [-1,0];
end
 
function grad_g4 = diff_g4(x,y)
    grad_g1 = [0,-1];
end