clc
%%====== Setup =========
S = [[3000 1000 2000];[1000 1000 4000];[2000 500 3000]];
sol = [];
c = S(:,4);
B = S; B(:,4) = [];
[sr sc] = size(B);
Axes = eye(sc);
%% ============ Find Interceptions ===============
A = [B(1,:);B(2,:)];
b = [c(1);c(2)];
sol = [sol linsolve(A,b)];
A = [B(1,:);B(3,:)];
b =[c(1);c(3)];
sol = [sol linsolve(A,b)];
A = [B(2,:);B(3,:)];
b = [c(2);c(3)];
sol = [sol linsolve(A,b)];
%% === Check for axes intersections ==
for i=1:sc
for j=1:sr
A = [Axes(i,:);B(j,:)];
b = [0;c(j)];
sol = [sol linsolve(A,b)];
end
end
%% === Find corner points of feasible
%% === region ======
[r c] = size(sol);
feasible = []; tolerance = 0;
for i=1:c
v = (sol(:,i));
x = v(1); y = v(2); z = v(3);
if (3000*x + 1000*y + 2000*z) <= 24000 + tolerance
&& (1000*x + 1000*y + 4000*z) <= 16000 + tolerance &&
(2000*x + 500*y + 3000*z) <= 48800 + tolerance
feasible = [feasible sol(:,i)];
end
end
disp(feasible)
%% ====== Find optimum =======