Hi,
I have developped a code for the problem shown in this video https://www.youtube.com/watch?v=1mCyC1FyMhk.
The problem is:
minimize f(x) NL
subject to A*X <= B
Here is the code:
% Example 1 :
% minimize: f(X) = x1^2 + x2^2 + x3^2 + x4^2 - 2x^1 - 3x4
% g(X) = AX <= B
% Where A = [Matrix] B = [7; 6; 0; 0; 0; 0]
% and all X >= 0, size(X) = (4, 1)
clear;
close;
clc;
%% Initial State : Choose an initial point that satisfy all the constraints
X0 = zeros(4, 1);
x10 = X0(1);
x20 = X0(2);
x30 = X0(3);
x40 = X0(4);
% Evaluate the Objective function, the Gradient of the Objective function and the constraint function at X = X0
ObjFun_X0 = x10^2 + x20^2 + x30^2 + x40^2 - 2*x10 - 3*x40;
Minus_Grad_ObjFun_X0 = [2 - 2*x10; -2*x20; -2*x30; 3-2*x40];
d_X0 = Minus_Grad_ObjFun_X0;
B = [7; 6; 0; 0; 0; 0];
A = [2 1 1 4; % A is the gradient of g(X)
1 1 2 1;
-1 0 0 0;
0 -1 0 0;
0 0 -1 0;
0 0 0 -1];
g_X0 = A*X0 - B; % the constraints equations are described by g(X) = AX <= B
% Find LC and TC
Test_C = A * X0;
index_AT = 0;
index_AL = 0;
LC = 0;
TC = 0;
for i = 1:6
% Tight Constraints
if ( Test_C(i) >= B(i) )
index_AT = index_AT + 1;
TC(index_AT) = i;
% Loose Constraints
else
index_AL = index_AL + 1;
LC(index_AL) = i;
end
end
TC = TC';
LC = LC';
% Get AL and AT
AL = zeros(size(LC, 1), 4);
AT = zeros(size(TC, 1), 4);
for j = 1:size(LC, 1)
AL (j, :) = A(LC(j), :);
end
for k = 1:size(TC, 1)
AT (k, :) = A(TC(k), :);
end
% find the Direction vector
Test_G = AT * d_X0;
dir_X0 = d_X0;
for e = 1:size(Test_G, 1)
if ( Test_G (e) > 0 )
Projection = eye(4) - AT' * (inv(AT * AT')) * AT;
dir_X0 = Projection * d_X0;
break;
end
end
% Compute the Langrange Multipliers for the tight Constraints : They must all postives > 0 to satify the optimality of the solution
LM = inv(AT * AT') * AT * d_X0;
count_LM = 0;
for i = 1:size(LM, 1)
if ( LM(i) <= 0 )
LM_Condition = 0;
count_LM = 0;
break;
end
if ( LM(i) > 0)
count_LM = count_LM + 1;
if ( count_LM >= size(LM, 1))
LM_Condition = 1;
end
end
end
count = 1;
Tolerance = 1e-8;
All_X(:, count) = X0;
%% ***** Begin the Main Loop of the Gradient Prjection method ***** %%
while ( norm(dir_X0) > Tolerance || LM_Condition~= 1 )
count = count + 1;
% Compute Lamda using the Table
All_Lamda = zeros(size(LC, 1), 1);
for i = 1 : size(LC, 1)
index = LC(i);
All_Lamda(i) = -g_X0(index) / ( AL(i, :) * dir_X0 );
end
All_Lamda_corrected = All_Lamda(All_Lamda >= 0);
Lamda = min(All_Lamda_corrected);
% Compute the New X
X = X0 + Lamda * dir_X0;
X0 = X;
% Store the solution at each step
All_X(:, count) = X0;
% Evaluate the ObjFun, Grad of the ObjFun and the Constraint Func
x10 = X0(1);
x20 = X0(2);
x30 = X0(3);
x40 = X0(4);
ObjFun_X0 = x10^2 + x20^2 + x30^2 + x40^2 - 2*x10 - 3*x40;
Minus_Grad_ObjFun_X0 = [2 - 2*x10; -2*x20; -2*x30; 3-2*x40];
d_X0 = Minus_Grad_ObjFun_X0;
g_X0 = A*X0 - B;
% Find LC and TC, we always compare to initial original constraints
Test_C = A * X0;
index_AT = 0;
index_AL = 0;
LC = 0;
TC = 0;
for i = 1:6
% Tight Constraints
if ( round(Test_C(i), 4) >= B(i) )
index_AT = index_AT + 1;
TC(index_AT) = i;
% Loose Constraints
else
index_AL = index_AL + 1;
LC(index_AL) = i;
end
end
TC = TC';
LC = LC';
% Get AL and AT
AL = zeros(size(LC, 1), 4);
AT = zeros(size(TC, 1), 4);
for j = 1:size(LC, 1)
AL (j, :) = A(LC(j), :);
end
for k = 1:size(TC, 1)
AT (k, :) = A(TC(k), :);
end
% Test if we need to project the gradient or not
Test_G = AT * d_X0;
dir_X0 = d_X0;
for e = 1:size(Test_G, 1)
if ( Test_G (e) > 0 )
Projection = eye(4) - AT' * (inv(AT * AT')) * AT;
dir_X0 = Projection * d_X0;
break;
end
end
% Compute the Langrange Multipliers for the tight Constraints
LM = inv(AT * AT') * AT * d_X0;
count_LM = 0;
for i = 1:size(LM, 1)
if ( LM(i) <= 0 )
LM_Condition = 0;
count_LM = 0;
break;
end
if ( LM(i) > 0)
count_LM = count_LM + 1;
if ( count_LM >= size(LM, 1))
LM_Condition = 1;
end
end
end
end % End While
