[X,fval,i]=frank_wolfe(f,X0,e,A,b,Aeq,beq,lb,ub)

function [X,fval,i]=frank_wolfe(f,X0,e,A,b,Aeq,beq,lb,ub)
%
% Yinxiao Li
%
% function function [X,fval,i]=frank_wolfe(f,X0,e,A,b,Aeq,beq,lb,ub)
%
% order = 3;
%
% Input:    f       cost function
%           X0      starting feasible point
%           e       stopping criteria
%           A       defined in linprog
%           b      
%           Aeq  
%           beq  
%           lb   
%           ub      defined in linprog
%
% Output:   X       optimal point
%           fval    cost at the optimal point
%           i       iterations

X=X0;
syms x1 x2 x3 lamida;
df_dx1=diff(f,x1);
df_dx2=diff(f,x2);
df_dx3=diff(f,x3);
x=[100;100;100];
g_f=[100;100;100];
i=0;
while(abs(g_f'*(x-X))>e)%stopping criteria
g_f1=subs(df_dx1,{x1,x2,x3},X);
g_f2=subs(df_dx2,{x1,x2,x3},X);
g_f3=subs(df_dx3,{x1,x2,x3},X);
g_f=[double(g_f1);double(g_f2);double(g_f3)];
[x,feval]=linprog(g_f,A,b,Aeq,beq,lb,ub);
f0=subs(f,{x1,x2,x3},X+lamida*(x-X));
f1=inline(char(f0));
[lamida,fval]=fminbnd(f1,0,1);
X1=X;
X=X+lamida*(x-X);
clear lamida;
syms lamida;
i=i+1;
end
fval=subs(f,{x1,x2,x3},X);