% Mixed Integer Linear Programming - example - version based on glpk
% 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 ->  1 solution is undefined
%               2 solution is feasible
%               3 solution is infeasible
%               4 no feasible solution exists
%               5 solution is optimal
%               6 solution is unbounded
%              ...

function [x_out,delta_out,fval,exitflag]=MILP_example_glpk(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

%% Definition of the problem contraints and cost function

% Define the constraints of the problem as Ax<=b
% note that here we have extra 2 columns associated with the binary
% variables
% moreover, the last 4 lines of the constraints are used to define
% delta(1)=1 <-> x(2)>0 and delta(2)=1 <-> x(3)>0 (see the slides of the
% related lectures for further details)

A = [minutes_matrix(1,1) minutes_matrix(2,1) minutes_matrix(3,1) 0 0;     % 2x1+4x2+8x3
    minutes_matrix(1,2) minutes_matrix(2,2) minutes_matrix(3,2) 0 0;     % 40x1+20x2+10x3
    minutes_matrix(1,3) minutes_matrix(2,3) minutes_matrix(3,3) 0 0;     % 4x1+8x2+12x3
    -1 0 0 0 0;                                                          % -x1
    0 -1 0 0 0;                                                          % -x2
    0 0 -1 0 0;                                                          % -x3
    0 -1  0 eps 0;                                                       % -x2+eps*delta_2 
    0 1  0 -120000 0;                                                    % x2-12000*delta_2 (chosen as upper bound of x2)
    0 0 -1 0  eps;                                                       % -x3+eps*delta_3 
    0 0 1 0  -120000];                                                   % x3-12000*delta_3 (chosen as upper bound of x3)                      

b = [qualified_pool;                                                  % 4000
    nonqualified_pool;                                               % 8000
    assemblage_pool;                                                 % 12000
    0;                                                               % 0
    0;                                                               % 0
    0;                                                               % 0
    0;                                                               % 0
    0;                                                               % 0
    0;                                                               % 0
    0];                                                              % 0

% Define the cost array - note: the last 2 coefficients are related to the
% binary variables
c = [-manual_profit -semi_profit -auto_profit fixed_cost_semi fixed_cost_auto];


%% Solve the MILP with GLPK

%msglev is one of the parameters of GLPK 
%(=0 for no visualization of intermediate results. help glpk for further details and further parameters)
param.msglev=0;


tic
% glpk is very similar to linprog. the main difference is in the string
% which allows to specify which optimization variables are continuous and
% which are binary
[xopt, fmin, status, extra]=glpk(c,A,b,[],[],[],'CCCBB',[],param);
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

% here we separate the optimal values of the continuous and discrete
% variables
x_out=xopt(1:3);
delta_out=xopt(4:5);

% since yalmip and glpk have different names, here we rename the flags and
% the optimal cost 
exitflag=status;
fval=fmin;



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

figure %here we plot
if exitflag>=2&&exitflag<100
    for i=1:4
        values=dec2bin(i-1,2);
        % the command below contains extra 2 equality constraints to fix
        % the value of the binary variables. This allows to show the
        % feasible domain of x1x2x3 for the different values of the
        % binaries
        poly(i)=Polyhedron('A',A,'b',b,'Ae',[0 0 0 1 0;0 0 0 0 1],'be',[bin2dec(values(1));bin2dec(values(2))]);
        % this command allows to project the obtained set in the 3
        % continuous variables x1x2x3 (note that in any case delta1 and
        % delta2 have been fixed to a single value)
        polyP(i)=poly(i).projection(1:3);
        subplot(2,2,i)
        fs1(i)=plot(polyP(i));
        if exitflag>=2&&exitflag<100&&not(exitflag==3)
            set(fs1(i),'FaceAlpha',0.3); %this command allows to make the figure a bit transparent
        end
        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 iso_plot&&exitflag>=2&&exitflag<100&&not(exitflag==3)%if we decided to plot the isocost lines we get into this loop
        for iii=1:4
            subplot(2,2,iii)
            i=0;
            for ii=0:2:10
                i=i+1;
                isocost{i}=Polyhedron('A',[-1 0 0;0 -1 0;0 0 -1],'b',[0;0;0],'Ae',-c(1:3),'be',ii*5000);
                id_p{i}=plot(isocost{i},'Color','y');
                set(id_p{i},'FaceAlpha',0.2);
                box on
                grid on
                
            end
        end
    end
    
end
% Plot a marker for the optimal vertex. this gets plotted only for the
% optimal value of the binary variables
if exitflag==5
    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


