%% Factory, Warehouse, Sales Allocation Model: Problem-Based
% This example shows how to set up and solve a mixed-integer linear programming problem. The problem is to find the optimal production and distribution levels among a set of factories, warehouses, and sales outlets. For the solver-based approach, see Factory, Warehouse, Sales Allocation Model: Solver-Based.
% The example first generates random locations for factories, warehouses, and sales outlets. Feel free to modify the scaling parameter N, which scales both the size of the grid in which the production and distribution facilities reside, but also scales the number of these facilities so that the density of facilities of each type per grid area is independent of N.
%% Facility Locations
% For a given value of the scaling parameter N, suppose that there are the following:
% \lfloor fN^2\rfloor factories
% \lfloor wN^2\rfloor warehouses
% \lfloor sN^2\rfloor sales outlets
% These facilities are on separate integer grid points between 1 and N in the x and y directions. In order that the facilities have separate locations, you require that f+w+s \le 1. In this example, take N = 20, f = 0.05, w = 0.05, and s = 0.1.
%% Production and Distribution
% There are P products made by the factories. Take P = 20.
% The demand for each product p in a sales outlet s is d(s,p). The demand is the quantity that can be sold in a time interval. One constraint on the model is that the demand is met, meaning the system produces and distributes exactly the quantities in the demand.
% There are capacity constraints on each factory and each warehouse.
% The production of product p at factory f is less than pcap(f,p).
% The capacity of warehouse w is wcap(w).
% The amount of product p that can be transported from warehouse w to a sales outlet in the time interval is less than turn(p)*wcap(w), where turn(p) is the turnover rate of product p.
% Suppose that each sales outlet receives its supplies from just one warehouse. Part of the problem is to determine the cheapest mapping of sales outlets to warehouses.
%% Costs
% The cost of transporting products from factory to warehouse, and from warehouse to sales outlet, depends on the distance between the facilities, and on the particular product. If dist(a,b) is the distance between facilities a and b, then the cost of shipping a product p between these facilities is the distance times the transportation cost tcost(p):
% dist(a,b)*tcost(p).
% The distance in this example is the grid distance, also known as the L_1 distance. It is the sum of the absolute difference in x coordinates and y coordinates.
% The cost of making a unit of product p in factory f is pcost(f,p).
%% Optimization Problem
% Given a set of facility locations, and the demands and capacity constraints, find:
% A production level of each product at each factory
% A distribution schedule for products from factories to warehouses
% A distribution schedule for products from warehouses to sales outlets
% These quantities must ensure that demand is satisfied and total cost is minimized. Also, each sales outlet is required to receive all its products from exactly one warehouse.
%% Variables and Equations for the Optimization Problem
% The control variables, meaning the ones you can change in the optimization, are
% x(p,f,w) = the amount of product p that is transported from factory f to warehouse w
% y(s,w) = a binary variable taking value 1 when sales outlet s is associated with warehouse w
% The objective function to minimize is
% \sum_f \sum_p \sum_w x(p,f,w)\ \cdot\ (pcost(f,p) + tcost(p)\cdot
% dist(f,w))
%  + \sum_s \sum_w \sum_p (d(s,p)\ \cdot\ tcost(p)\ \cdot\ dist(s,w)\ \cdot\ y(s,w)).
% The constraints are
% \sum_w x(p,f,w) \le pcap(f,p) (capacity of factory).
% \sum_f x(p,f,w) = \sum_s (d(s,p)\ \cdot\ y(s,w)) (demand is met).
% \sum_p \sum_s \frac{d(s,p)}{turn(p)}\ \cdot\ y(s,w) \le wcap(w) (capacity of warehouse).
% \sum_w y(s,w) = 1 (each sales outlet associates to one warehouse).
% x(p,f,w) \ge 0 (nonnegative production).
% y(s,w)\ \epsilon\ \{0,1\} (binary y).
% The variables x and y appear in the objective and constraint functions linearly. Because y is restricted to integer values, the problem is a mixed-integer linear program (MILP).
%% Generate a Random Problem: Facility Locations
% Set the values of the N, f, w, and s parameters, and generate the facility locations.

rng(1) % for reproducibility
N = 20; % N from 10 to 30 seems to work. Choose large values with caution.
N2 = N*N;
f = 0.05; % density of factories
w = 0.05; % density of warehouses
s = 0.1; % density of sales outlets

F = floor(f*N2); % number of factories
W = floor(w*N2); % number of warehouses
S = floor(s*N2); % number of sales outlets

xyloc = randperm(N2,F+W+S); % unique locations of facilities
[xloc,yloc] = ind2sub([N N],xyloc);

%
% Of course, it is not realistic to take random locations for facilities. This example is intended to show solution techniques, not how to generate good facility locations.
% Plot the facilities. Facilities 1 through F are factories, F+1 through F+W are warehouses, and F+W+1 through F+W+S are sales outlets.

h = figure;
plot(xloc(1:F),yloc(1:F),'rs',xloc(F+1:F+W),yloc(F+1:F+W),'k*',...
    xloc(F+W+1:F+W+S),yloc(F+W+1:F+W+S),'bo');
lgnd = legend('Factory','Warehouse','Sales outlet','Location','EastOutside');
lgnd.AutoUpdate = 'off';
xlim([0 N+1]);ylim([0 N+1])

%
%% Generate Random Capacities, Costs, and Demands
% Generate random production costs, capacities, turnover rates, and demands.

P = 20; % 20 products

% Production costs between 20 and 100
pcost = 80*rand(F,P) + 20;

% Production capacity between 500 and 1500 for each product/factory
pcap = 1000*rand(F,P) + 500;

% Warehouse capacity between P*400 and P*800 for each product/warehouse
wcap = P*400*rand(W,1) + P*400;

% Product turnover rate between 1 and 3 for each product
turn = 2*rand(1,P) + 1;

% Product transport cost per distance between 5 and 10 for each product
tcost = 5*rand(1,P) + 5;

% Product demand by sales outlet between 200 and 500 for each
% product/outlet
d = 300*rand(S,P) + 200;

% These random demands and capacities can lead to infeasible problems. In other words, sometimes the demand exceeds the production and warehouse capacity constraints. If you alter some parameters and get an infeasible problem, during solution you will get an exitflag of -2.
%
%% Generate Variables and Constraints
% To begin specifying the problem, generate the distance arrays distfw(i,j) and distsw(i,j).

distfw = zeros(F,W); % Allocate matrix for factory-warehouse distances
for ii = 1:F
    for jj = 1:W
        distfw(ii,jj) = abs(xloc(ii) - xloc(F + jj)) + abs(yloc(ii) ...
            - yloc(F + jj));
    end
end

distsw = zeros(S,W); % Allocate matrix for sales outlet-warehouse distances
for ii = 1:S
    for jj = 1:W
        distsw(ii,jj) = abs(xloc(F + W + ii) - xloc(F + jj)) ...
            + abs(yloc(F + W + ii) - yloc(F + jj));
    end
end

%
% Create variables for the optimization problem. x represents the production, a continuous variable, with dimension P-by-F-by-W. y represents the binary allocation of sales outlet to warehouse, an S-by-W variable.

x = optimvar('x',P,F,W,'LowerBound',0);
y = optimvar('y',S,W,'Type','integer','LowerBound',0,'UpperBound',1);

%
% Now create the constraints. The first constraint is a capacity constraint on production.

capconstr = sum(x,3) <= pcap';

%
% The next constraint is that the demand is met at each sales outlet.

demconstr = squeeze(sum(x,2)) == d'*y;

%
% There is a capacity constraint at each warehouse.

warecap = sum(diag(1./turn)*(d'*y),1) <= wcap';

%
% Finally, there is a requirement that each sales outlet connects to exactly one warehouse.

salesware = sum(y,2) == ones(S,1);

%
%% Create Problem and Objective
% Create an optimization problem.

factoryprob = optimproblem;

%
% The objective function has three parts. The first part is the sum of the production costs.

objfun1 = sum(sum(sum(x,3).*(pcost'),2),1);

%
% The second part is the sum of the transportation costs from factories to warehouses.

objfun2 = 0;
for p = 1:P
    objfun2 = objfun2 + tcost(p)*sum(sum(squeeze(x(p,:,:)).*distfw));
end

%
% The third part is the sum of the transportation costs from warehouses to sales outlets.

r = sum(distsw.*y,2); % r is a length s vector
v = d*(tcost(:));
objfun3 = sum(v.*r);

%
% The objective function to minimize is the sum of the three parts.

factoryprob.Objective = objfun1 + objfun2 + objfun3;

%
% Include the constraints in the problem.

factoryprob.Constraints.capconstr = capconstr;
factoryprob.Constraints.demconstr = demconstr;
factoryprob.Constraints.warecap = warecap;
factoryprob.Constraints.salesware = salesware;

%
%% Solve the Problem
% Turn off iterative display so that you don't get hundreds of lines of output. Include a plot function to monitor the solution progress.

opts = optimoptions('intlinprog','Display','off','PlotFcn',@optimplotmilp);

%
% Call the solver to find the solution.

[sol,fval,exitflag,output] = solve(factoryprob,'options',opts);

%

if isempty(sol) % If the problem is infeasible or you stopped early with no solution
    disp('The solver did not return a solution.')
    return % Stop the script because there is nothing to examine
end

%
%% Examine the Solution
% Examine the exit flag and the infeasibility of the solution.

exitflag
infeas1 = max(max(infeasibility(capconstr,sol)))
infeas2 = max(max(infeasibility(demconstr,sol)))
infeas3 = max(infeasibility(warecap,sol))
infeas4 = max(infeasibility(salesware,sol))

%
% Round the y portion of the solution to be exactly integer-valued. To understand why these variables might not be exactly integers, see Some "Integer" Solutions Are Not Integers.

sol.y = round(sol.y); % get integer solutions

% How many sales outlets are associated with each warehouse? Notice that, in this case, some warehouses have 0 associated outlets, meaning the warehouses are not in use in the optimal solution.

outlets = sum(sol.y,1)

%
% Plot the connection between each sales outlet and its warehouse.

figure(h);
hold on
for ii = 1:S
    jj = find(sol.y(ii,:)); % Index of warehouse associated with ii
    xsales = xloc(F+W+ii); ysales = yloc(F+W+ii);
    xwarehouse = xloc(F+jj); ywarehouse = yloc(F+jj);
    if rand(1) < .5 % Draw y direction first half the time
        plot([xsales,xsales,xwarehouse],[ysales,ywarehouse,ywarehouse],'g--')
    else % Draw x direction first the rest of the time
        plot([xsales,xwarehouse,xwarehouse],[ysales,ysales,ywarehouse],'g--')
    end
end
hold off

title('Mapping of sales outlets to warehouses')

% The black * with no green lines represent the unused warehouses.
% Copyright 2014-2017 The MathWorks, Inc.