% Mixed Integer Linear Programming - example - version based on yalmip
% script designed by Dr. Marcello Torchio and Prof. Davide M. Raimondo
% (davide.raimondo@unipv.it)

% inputs:
%  iso_plot -> 1 plots isocost lines - any other value does not plot

% outputs:
%  x_out -> provides the optimal value of the continuous variables
%  delta_out -> provides the optimal value of the binary variables
%  fval -> optimal cost
%  exitflag -> 0 Successfully solved
%              1 Infeasible problem
%              2 Unbounded objective function
%              ...

function [x_out,delta_out,fval,exitflag]=MILP_example_yalmip(iso_plot)

% Show decimals according to the short format
format short

%% Definition of the problem variables

cost_to_sell=90;

% Washing machines production costs
manual_cost = 45;
semi_cost   = 40;
auto_cost   = 35;

%fixed costs
fixed_cost_semi=4000;
fixed_cost_auto=6000;

% 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      = 4000;
nonqualified_pool   = 8000;
assemblage_pool     = 12000;

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

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

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

%% Define the problem variables

% Declare the yalmip continous variables x1,x2, and x3 
x = sdpvar(3,1);
% Declare the yalmip binary variables delta2,delta3
delta=binvar(2,1);

%% Definition of the problem contraints and cost function

% "implies" is used to associate delta(1)=1 to x(2)>0. 
% the same is done for delta(2) and x(3)
constraints = [implies(x(2) >=eps,delta(1));implies(x(3) >= eps,delta(2))];

% 2*x1+4*x2+8*x3          <= 4000
constraints = [constraints;minutes_matrix(1,1)*x(1)+minutes_matrix(2,1)*x(2)+minutes_matrix(3,1)*x(3)<=qualified_pool];
% 40*x1+20*x2+10*x3     <= 8000
constraints = [constraints;minutes_matrix(1,2)*x(1)+minutes_matrix(2,2)*x(2)+minutes_matrix(3,2)*x(3)<=nonqualified_pool];
% 4*x1+8*x2+12*x3       <= 12000
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;x(2)<=qualified_pool/minutes_matrix(2,1);x(2)<=nonqualified_pool/minutes_matrix(2,2);x(2)<=assemblage_pool/minutes_matrix(2,3);x(3)<=qualified_pool/minutes_matrix(3,1);x(3)<=nonqualified_pool/minutes_matrix(3,2);x(3)<=assemblage_pool/minutes_matrix(3,3)];


% Define the cost array
% note that since the optimization will minimize the cost we negated the
% cost which we want to maximize
% note the fixed costs associated to delta(1) and delta(2)
cost_function =- manual_profit*x(1) - semi_profit*x(2) - auto_profit *x(3)+fixed_cost_semi*delta(1)+fixed_cost_auto*delta(2);

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

% Choose the solver to be used: glpk solves MILPs
yalmipopt.solver = 'glpk';
yalmipopt.verbose=0; %this is to avoid yalmip to show the iterations of glpk

%% 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); %here we separate the optimal values of the continous and binary variables
delta_out=double(delta);

% 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);
    for i=1:length(feasible_set) %we loop over the 4 possible values of delta1delta2
        subplot(2,2,i)
        fs1(i)=plot(feasible_set(i));
        title(sprintf('delta_1delta_2=%s',char(dec2bin(i-1,2))));
        xlabel('Manual units')
        ylabel('Semi-automatic units')
        zlabel('Automatic units')
        hold on
        
    end
    if exitflag>=-2&&exitflag<=1
        try
            for i=1:length(feasible_set)
                set(fs1(i),'FaceAlpha',0.3); %this command allows to make the figure a bit transparent
            end
        catch
        end
    end
    
    
    if iso_plot&&exitflag>=-2&&exitflag<=1%if we decided to plot the isocost lines we get into this loop
        for iii=1:length(feasible_set)
            subplot(2,2,iii)
            i=0;
            for ii=0:2:10
                i=i+1;
                cost_function =- manual_profit*x(1) - semi_profit*x(2) - auto_profit *x(3);
                isocost_constraints{i}=[-cost_function==ii*5000;x>=0];
                isocost{i}=Polyhedron(isocost_constraints{i});
                id_p{i}=plot(isocost{i},'Color','y');
                set(id_p{i},'FaceAlpha',0.2);
                box on
                grid on
                
            end
        end
    end
    % Plot a marker for the optimal vertex. this gets plotted only for the
    % optimal value of the binary variables
    if exitflag==0
        subplot(2,2,bin2dec([num2str(delta_out(1)) num2str(delta_out(2))])+1)
        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
end
