% min S(x) = -6*x1 - 2*x2 + 4 * x1^2 + 8 * x2^2
%  s.a. g1(x) = 2*x1 + x2 - 8 <= 0
%       g2(x) = 4*x1 + x2 - 10 <= 0
%        x >= 0
%  x' = [1  6]

Q=[4 0;0 8]*2
c=[-6; -2]
A=[2 1;4 1]
b=[8; 10]
lb=[0; 0]
ub=[]

op=optimset('LargeScale','off');
[x,s,ex,out,lambda]=quadprog(Q,c,A,b,[],[],lb,ub,[],op)

[X1,X2]=meshgrid([0:0.1:4],[0:0.1:10]);
S=c(1)*X1+c(2)*X2 + (Q(1,1)*X1.^2 + Q(2,2)* X2.^2)/2;
l=contour(X1,X2,S,[-2.3 -1.5 0 5 10 20 40 80 160 320 500 640]); hold on
clabel(l);
plot(X1(1,:),b(1)-A(1,1)*X1(1,:),'g',X1(1,:),b(2)-A(2,1)*X1(1,:),'r')
xlabel('x_1'); ylabel('x_2');

Al=[A eye(2,2) zeros(2,2) zeros(2,2);Q zeros(2,2) -A' eye(2,2)]
bl=[b;-c]
cl=0.5*[c;0;0;b;0;0]
LB=[0;0;0;0;0;0;0;0];
UB=[];
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=linprog(cl,[],[],Al,bl,LB,UB,[],op)
