% Linear Programming example - version based on yalmip
% 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_out,fval,exitflag]=LP_example_yalmip(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

%% Define the problem variables

% Declare the yalmip variables x1,x2, and x3 (the optimization variables)
x = sdpvar(3,1);
% x = intvar(3,1);

%% Definition of the problem contraints and cost function

% Initialize the constraints array as an empty variable.
constraints = [];

% x1+x2+x3 >= minimum_number_of_machines
constraints = [constraints;x(1)+x(2)+x(3)>=minimum_number_of_machines];
% x1+4*x2+8*x3          <= 4500
constraints = [constraints;minutes_matrix(1,1)*x(1)+minutes_matrix(2,1)*x(2)+minutes_matrix(3,1)*x(3)<=qualified_pool];
% 40*x1+30*x2+20*x3     <= 36000
constraints = [constraints;minutes_matrix(1,2)*x(1)+minutes_matrix(2,2)*x(2)+minutes_matrix(3,2)*x(3)<=nonqualified_pool];
% 3*x1+2*x2+40*x3       <= 2700
constraints = [constraints;minutes_matrix(1,3)*x(1)+minutes_matrix(2,3)*x(2)+minutes_matrix(3,3)*x(3)<=assemblage_pool];
% x1,x2,x3>=0
constraints = [constraints;x>=0];


% Define the cost array
if strcmp(opt,'min')==1 % minimize costs
    cost_function = manual_cost*x(1) + semi_cost*x(2) + auto_cost *x(3);
else % maximize profits (which is equal to minimize -profit)
    cost_function =- manual_profit*x(1) - semi_profit*x(2) - auto_profit *x(3);
end

%% Yalmip options
% Define the yalmip options structure
yalmipopt = sdpsettings;

% Choose the solver to be used
yalmipopt.solver = 'linprog';

% For the chosen solver, define the numerical algorithm
yalmipopt.linprog.Algorithm = 'dual-simplex';

%% Solve the problem

tic
result      = solvesdp(constraints,cost_function,yalmipopt);
timeout     = toc; % tic and toc allow to compute the time required by linprog to solve the optimization

% Obtain the yalmip exitflag. See help yalmiperror for further details
exitflag    = result.problem;

% Display the text message corresponding to the yalmip exit flag
disp(yalmiperror(exitflag))

% Yalmip has a different set of error flags
%0 Successfully solved
%1 Infeasible problem
%2 Unbounded objective function
%...
% help yalmiperror for the entire list

% to obtain the optimal cost and value of the optimization variables one
% has to use the double()
fval= double(cost_function);
x_out=double(x);

% Print the solution time
disp(sprintf('Solution time %f',timeout));



figure
if exitflag>=-2&&exitflag<=2
    % Define the polyhedron representing the feasible set
    feasible_set =Polyhedron(constraints);
    fs1=plot(feasible_set);
    if exitflag>=-2&&exitflag<=1
        try
            set(fs1,'FaceAlpha',0.3); %this command allows to make the figure a bit transparent
        catch
        end
    end
    
    xlabel('Manual units')
    ylabel('Semi-automatic units')
    zlabel('Automatic units')
    hold on
    
    if iso_plot&&exitflag>=-2&&exitflag<=1
        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_constraints{i}=[cost_function==ii*10000;x>=0];
                isocost{i}=Polyhedron(isocost_constraints{i});
            else % the constraints below provide -c'x=ii*10000 and [x1;x2;x3]>=0
                isocost_constraints{i}=[-cost_function==ii*10000;x>=0];
                isocost{i}=Polyhedron(isocost_constraints{i});
            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);
keyboard
    % the command feasible_set.V provides the list containing all the
    % vertices of feasible_set
    
    if exitflag== ( 0 )
        plot3(x_out(1),x_out(2),x_out(3),'*')
        plot3(x_out(1),x_out(2),x_out(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