%% Linear inequality constraints and quadratic cost
% Maximize 14x_1-x_1^2+6x_2-x_2^2+7
% subject to x_1+x_2\leq 2,\;x_1+2x_2\leq 3,\; x_1\geq 0,\; x_2\geq 0
%% Plotting the feasibility region given by a set of linear inequalities  and on top of this also the contours of the (quadratic) cost function
% This needs a convenience third-party m-file: Per Bergström (2020). Plot 2D/3D region (https://www.mathworks.com/matlabcentral/fileexchange/9261-plot-2d-3d-region), MATLAB Central File Exchange. Retrieved February 19, 2020.

A = [1 1; 1 2]; 
b = [2 3]';
lb = [0; 0];

plotregion(-A,-b,lb)
xlabel('x_1'), ylabel('x_2')

X1 = -1:0.1:3; X2 = -1:0.1:2; [X1,X2] = meshgrid(X1,X2);
L = 14*X1-X1.^2+6*X2-X2.^2+7;
hold on,
[C,h] = contour(X1,X2,L);
set(h,'ShowText','on','TextStep',get(h,'LevelStep')/10);

%
%% Computing the gradients of the cost function and the constraints symbolically

syms x1 x2 lambda1 lambda2 lambda3 lambda4

f = -(14*x1+x1^2+6*x2+x2^2+7);
h1 = x1+x2-2;
h2 = x1+1*x2-3;
h3 = -x1;
h4 = -x2;

gf = gradient(f,[x1,x2])
gh1 = gradient(h1,[x1,x2])
gh2 = gradient(h2,[x1,x2])
gh3 = gradient(h3,[x1,x2])
gh4 = gradient(h4,[x1,x2])
gh = [gh1 gh2 gh3 gh4]

%
%% lambda1 = 0, lambda2 = 0, lambda3 = 0, lambda4 = 0

sol = solve(gf)
xst = [sol.x1; sol.x2]

%
%% ...
%% lambda1 = 1, lambda2 = 0, lambda3 = 0, lambda4 = 1

sol = solve([gf+lambda1*gh1+lambda4*gh4;h1;h4])
xst = [sol.x1; sol.x2]

%
%% Computing the stationary points numerically

Q = [-2 0; 0 -2]
r = [14; 6]
s = 7
Q = -Q
r = -r
s = -s

%
%% Using quadprog from Optimization Toolbox

xqp = quadprog(Q,r,A,b,[],[],lb,[])