% Linear Programming example - version based on linprog
% script designed by Dr. Marcello Torchio and Prof. Davide M. Raimondo
% (davide.raimondo@unipv.it)
% inputs: 
%  minimum_number_of_machines -> lower bound on the number of machines to be produced
%  iso_plot -> 1 plots isocost lines - any other value does not plot 
%  opt -> 'min' solves the problem which minimizes costs min 70x1+80x2+85x3
%         'max' solves the problem which minimizes profits max (130-70)x1+(130-80)x2+(130-85)x3
%

% outputs: 
%  x -> provides the optimal value of the optimization variables
%  fval -> optimal cost
%  exitflag -> 1 optimization converged to a solution
%              0  Maximum number of iterations reached.
%             -2  No feasible point found.
%             -3  Problem is unbounded.
%              ...
%

function [x,fval,exitflag]=LP_example_linprog(minimum_number_of_machines,iso_plot,opt)

% Show decimals according to the short format
format short

%% Definition of the problem variables

cost_to_sell=130;

% Washing machines production costs
manual_cost = 70;
semi_cost   = 80;
auto_cost   = 85;

% Washing machines production profits
manual_profit = cost_to_sell-manual_cost;
semi_profit   = cost_to_sell-semi_cost;
auto_profit   = cost_to_sell-auto_cost;

% Pools of minutes for the different types of work.
qualified_pool      = 4500;
nonqualified_pool   = 36000;
assemblage_pool     = 2700;

% Minutes for the manual production
minutes_matrix(1,1) = 1; % qualified work
minutes_matrix(1,2) = 40;% non-qualified work
minutes_matrix(1,3) = 3; % assemblage

% Minutes for the semi-automatic production
minutes_matrix(2,1) = 4; % qualified work
minutes_matrix(2,2) = 30;% non-qualified work
minutes_matrix(2,3) = 2; % assemblage

% Minutes for the automatic production
minutes_matrix(3,1) = 8; % qualified work
minutes_matrix(3,2) = 20;% non-qualified work
minutes_matrix(3,3) = 4;% assemblage

%% Definition of the problem contraints and cost function

% Define the constraints of the problem as Ax<=b

A = [-1 -1 -1                                                        % x1+x2+x3
    minutes_matrix(1,1) minutes_matrix(2,1) minutes_matrix(3,1);     % x1+4x2+8x3
    minutes_matrix(1,2) minutes_matrix(2,2) minutes_matrix(3,2);     % 40x1+30x2+20x3
    minutes_matrix(1,3) minutes_matrix(2,3) minutes_matrix(3,3);     % 3x1+2x2+40x3
    -1 0 0;                                                          % -x1
    0 -1 0;                                                          % -x2
    0 0 -1;];                                                        % -x3

b = [-minimum_number_of_machines;
    qualified_pool;                                                  % 4500
    nonqualified_pool;                                               % 36000
    assemblage_pool;                                                 % 2700
    0;                                                               % 0
    0;                                                               % 0
    0];                                                              % 0

% Define the cost array
if strcmp(opt,'min')==1 % minimize costs
    c = [manual_cost semi_cost auto_cost];
else % maximize profits (which is equal to minimize -profit)
    c = -[manual_profit semi_profit auto_profit];
end


%% Solve the linprog

% Define the options structure for the linprog solver. Choose the simplex
% method to solve the optimization problem.
opts = optimoptions(@linprog,'Algorithm','dual-simplex');
% opts=optimoptions(@linprog,'Algorithm','interior-point');

% Solve the linear problem
% x is the optimal value of the optimization variables
% fval contains the optimal value of the cost
% exitflag tells whether the solution was optimal, or the problem feasible,
% unbounded, etc.

tic
[x,fval,exitflag] = linprog(c,A,b,[],[],[],[],[],opts);
timeout=toc; % tic and toc allow to compute the time required by linprog to solve the optimization

disp(sprintf('Solution time %f',timeout)); %disp diplays on the screen the value of timeout


% Display a text warning according to the exitflag value. See help linprog
% for furhter details
if exitflag==1
    disp('Optimal solution found');
elseif exitflag==0
    disp('Maximum number of iterations reached');
elseif  exitflag==-2
    disp('No feasible point found')
elseif exitflag==-3
    disp('Problem is unbounded');
else
    disp('Some problems occured, please check the value of the variable exitflag')
end


figure
if exitflag>=-3
    % Define the polyhedron representing the feasible set
    feasible_set = Polyhedron(A,b);
    fs1          = plot(feasible_set);
    
    if exitflag>=-2
        try
        set(fs1,'FaceAlpha',0.3); %this command allows to make the figure a bit transparent
        catch
        end
    end
    hold on
    xlabel('Manual units')
    ylabel('Semi-automatic units')
    zlabel('Automatic units')
    % Plot the isocost lines when required
    if iso_plot&&exitflag>=-2
        i=0;
        for ii=0:2:10
            i=i+1;
            if strcmp(opt,'min') % the constraints below provide c'x=ii*10000 and [x1;x2;x3]>=0
            isocost{i}=Polyhedron('A',[-1 0 0;0 -1 0;0 0 -1],'b',[0;0;0],'Ae',c,'be',ii*10000);
            else % the constraints below provide -c'x=ii*10000 and [x1;x2;x3]>=0
                isocost{i}=Polyhedron('A',[-1 0 0;0 -1 0;0 0 -1],'b',[0;0;0],'Ae',-c,'be',ii*10000);
            end
            id_p{i}=plot(isocost{i},'Color','y');
            set(id_p{i},'FaceAlpha',0.2);
        end
    end
    
    % Plot a marker for each vertex of the feasible set
    plot3(feasible_set.V(:,1),feasible_set.V(:,2),feasible_set.V(:,3),'og','Markersize',10,'LineWidth',4);
    % the command feasible_set.V provides the list containing all the
    % vertices of feasible_set
    
    
    if exitflag<=-2
        ;
    else
        % Plot with a different color and marker the vertex for which we have the optimal solution
        plot3(x(1),x(2),x(3),'*')
        plot3(x(1),x(2),x(3),'or','Markersize',20,'LineWidth',6)
        
    end
        
    axis tight
    box on
    grid on
    
    % Put the title to the plot
    if strcmp(opt,'min')
        title('min')
    else
        title('max')
    end
end