function  TRM4 
% All Rights Reserved lilongduzhi@yahoo.com.cn hehe,:) 
% Optimal Function : min = -15;X = [0 ,3 ,0 ,4];
% myfun = x1-x2-x3-x1*x3+x1*x4+x2*x3-x2*x4;
% st:  x1 + 2x2 <= 8 ;
%     4x1 + x2  <= 12 ;
%     3x1 + 4x2 <= 12 ;
%     2x3 - x4  <= 8 ;
%     x3 + 2x4  <= 8 ;
%     x3 + x4   <= 5 ;
%     x1,x2,x3  >= 0 ;
% diff_g(X) = [ 1-x3+x4 ,-1+x3-x4 ,-1-x1+x2 ,x1-x2];
nvars = 4; % 變量個數 ；
A = [ 1  ,2 ,0 ,0;
      4  ,1 ,0 ,0;
      3  ,4 ,0 ,0;
      0  ,0 ,2 ,-1;
      0  ,0 ,1 ,2;
      0  ,0 ,1 ,1;];
A = [-A;eye(3),zeros(3,1)] ; 
b = [ -8 ,-12 ,-12 ,-8 ,-8 ,-5 ,0 ,0 ,0]';
% step 1:
delta = 0.9; % [0 ,1]
epsilon = 1e-6;
X = [1 ,1 ,1 ,1]'; % X=[x1 ,x2];
B = eye(nvars);
k = 0;
%  step 2:
grad = diff_g(X);
sqp_A = -A;
sqp_b = A*X-b;
upper = delta*X;
downer = -delta*X;
d = quadprog(B ,grad ,sqp_A ,sqp_b ,[],[] ,-delta*X ,delta*X);
while max(abs(d)) > epsilon
    k = k+1 ;
    % step 2:
    norm_d = max(abs(d))
    sqp_b =  A*X-b;
    lower = -abs(delta*X);
    upper =  abs(delta*X);
    % d = quadprog(B ,grad ,sqp_A,sqp_b,[],[] ,-delta*X ,delta*X);
      d = quadprog(B ,grad ,sqp_A,sqp_b,[],[] ,lower ,upper);
    % step 2.2:
    grad = diff_g(X);
    ared = myfun(X)-myfun(X+d);
    pred = q(zeros(nvars,1) ,grad ,B)-q(d ,grad ,B);
    r = ared/pred;
    % step 2.1:
    if r > 0
        X = X+d;
    else
        X = X;
    end
    % step 3:
    if r>=0.25 ,
        % goto step 4:
    else delta = delta/2;
    end
    if r<0.75 || max(abs(d)) < delta ,
        % goto step 5:
    else
        delta = min(2*delta ,1-1e-10);
    end
    % step 5:
    delta = delta ;
    % 5.1 : mod B
    B = BFGS(B ,X ,d);
    k = k+1
end
disp('求解結果：');
X
min_value = myfun(X)
% =====================> sub function <=====================%
function y = myfun(X)
% myfun = x1-x2-x3-x1*x3+x1*x4+x2*x3-x2*x4;
y = X(1)-X(2)-X(3)-X(1)*X(3)+X(1)*X(4)+ X(2)*X(3)-X(2)*X(4);

function grad = diff_g(X)
% diff_g(X) = [ 1-x3+x4 ,-1+x3-x4 ,-1-x1+x2 ,x1-x2];
grad = [ 1-X(3)+X(4) ,-1+X(3)-X(4) ,-1-X(1)+X(2) ,X(1)-X(2)]';

function Q = q(d ,g ,B)
Q = g'*d+1/2*(d'*B*d);

function B = BFGS(B ,X,d)
y = diff_g(X+d)-diff_g(X);
s = d;
if y'*s > 0
    B = B - B*s*s'*B/(s'*B*s)+y*y'/(y'*s);
end
%=============================================================%