Factory, Warehouse, Sales Allocation Model: Solver-Based

Facility Locations

For a given value of the scaling parameter $$N$$, suppose that there are the following:

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.

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)$$:

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:

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

The objective function to minimize is

The constraints are

$$\underset{w}{\sum }x(p,f,w)\le pcap(f,p)$$(capacity of factory).

$$\underset{f}{\sum }x(p,f,w)=\underset{s}{\sum }(d(s,p)\cdot y(s,w))$$(demand is met).

$$\underset{p}{\sum }\underset{s}{\sum }\frac{d(s,p)}{turn(p)}\cdot y(s,w)\le wcap(w)$$(capacity of warehouse).

$$\underset{w}{\sum }y(s,w)=1$$(each sales outlet associates to one warehouse).

$$x(p,f,w)\ge 0$$(nonnegative production).

$$y(s,w)ϵ\{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 reproducibilityN = 20;% N从10到30似乎工作。选择大的价值es with caution.N2 = N*N; f = 0.05;% density of factoriesw = 0.05;% density of warehousess = 0.1;% density of sales outletsF = floor(f*N2);% number of factoriesW = floor(w*N2);% number of warehousesS = floor(s*N2);% number of sales outletsxyloc = 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 100pcost = 80*rand(F,P) + 20;% Production capacity between 500 and 1500 for each product/factorypcap = 1000*rand(F,P) + 500;% Warehouse capacity between P*400 and P*800 for each product/warehousewcap = P*400*rand(W,1) + P*400;% Product turnover rate between 1 and 3 for each productturn = 2*rand(1,P) + 1;% Product transport cost per distance between 5 and 10 for each producttcost = 5*rand(1,P) + 5;% Product demand by sales outlet between 200 and 500 for each% product/outletd = 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 Objective and Constraint Matrices and Vectors

The objective function vectorobjinintlinconconsists of the coefficients of the variables $$x(p,f,w)$$and $$y(s,w)$$．因此,re are naturallyP*F*W + S*Wcoefficients inobj．

One way to generate the coefficients is to begin with aP-by-F-by-Warrayobj1for the $$x$$coefficients, and anS-by-Warrayobj2for the $$y(s,w)$$coefficients. Then convert these arrays to two vectors and combine them intoobjby calling

obj = [obj1(:);obj2(:)];

obj1 = zeros(P,F,W);% Allocate arraysobj2 = zeros(S,W);

Throughout the generation of objective and constraint vectors and matrices, we generate the $$(p,f,w)$$array or the $$(s,w)$$array, and then convert the result to a vector.

To begin generating the inputs, generate the distance arraysdistfw(i,j)anddistsw(i,j)．

distfw = zeros(F,W);% Allocate matrix for factory-warehouse distancesforii = 1:Fforjj = 1:W distfw(ii,jj) = abs(xloc(ii) - xloc(F + jj)) + abs(yloc(ii)．..- yloc(F + jj));endenddistsw = zeros(S,W);% Allocate matrix for sales outlet-warehouse distancesforii = 1:Sforjj = 1:W distsw(ii,jj) = abs(xloc(F + W + ii) - xloc(F + jj))．..+ abs(yloc(F + W + ii) - yloc(F + jj));endend

Generate the entries ofobj1andobj2．

forii = 1:Pforjj = 1:Fforkk = 1:W obj1(ii,jj,kk) = pcost(jj,ii) + tcost(ii)*distfw(jj,kk);endendendforii = 1:Sforjj = 1:W obj2(ii,jj) = distsw(ii,jj)*sum(d(ii,:).*tcost);endend

Combine the entries into one vector.

obj = [obj1(:);obj2(:)];% obj is the objective function vector

Now create the constraint matrices.

The width of each linear constraint matrix is the length of theobjvector.

matwid = length(obj);

There are two types of linear inequalities: the production capacity constraints, and the warehouse capacity constraints.

There areP*Fproduction capacity constraints, andWwarehouse capacity constraints. The constraint matrices are quite sparse, on the order of 1% nonzero, so save memory by using sparse matrices.

Aineq = spalloc(P*F + W,matwid,P*F*W + S*W);% Allocate sparse Aeqbineq = zeros(P*F + W,1);% Allocate bineq as full% Zero matrices of convenient sizes:clearer1 = zeros(size(obj1)); clearer12 = clearer1(:); clearer2 = zeros(size(obj2)); clearer22 = clearer2(:);% First the production capacity constraintscounter = 1;forii = 1:Fforjj = 1:P xtemp = clearer1; xtemp(jj,ii,:) = 1;% Sum over warehouses for each product and factoryxtemp = sparse([xtemp(:);clearer22]);% Convert to sparseAineq(counter,:) = xtemp';% Fill in the rowbineq(counter) = pcap(ii,jj); counter = counter + 1;endend%现在仓库容量约束vj = zeros(S,1);% The multipliersforjj = 1:S vj(jj) = sum(d(jj,:)./turn);% A sum of P elementsendforii = 1:W xtemp = clearer2; xtemp(:,ii) = vj; xtemp = sparse([clearer12;xtemp(:)]);% Convert to sparseAineq(counter,:) = xtemp';% Fill in the rowbineq(counter) = wcap(ii); counter = counter + 1;end

There are two types of linear equality constraints: the constraint that demand is met, and the constraint that each sales outlet corresponds to one warehouse.

Aeq = spalloc(P*W + S,matwid,P*W*(F+S) + S*W);% Allocate as sparsebeq = zeros(P*W + S,1);% Allocate vectors as fullcounter = 1;% Demand is satisfied:forii = 1:Pforjj = 1:W xtemp = clearer1; xtemp(ii,:,jj) = 1; xtemp2 = clearer2; xtemp2(:,jj) = -d(:,ii); xtemp = sparse([xtemp(:);xtemp2(:)]');% Change to sparse rowAeq(counter,:) = xtemp;% Fill in rowcounter = counter + 1;endend% Only one warehouse for each sales outlet:forii = 1:S xtemp = clearer2; xtemp(ii,:) = 1; xtemp = sparse([clearer12;xtemp(:)]');% Change to sparse rowAeq(counter,:) = xtemp;% Fill in rowbeq(counter) = 1; counter = counter + 1;end

Bound Constraints and Integer Variables

The integer variables are those fromlength(obj1) + 1to the end.

intcon = P*F*W+1:length(obj);

The upper bounds are fromlength(obj1) + 1to the end also.

lb = zeros(length(obj),1); ub = Inf(length(obj),1); ub(P*F*W+1:end) = 1;

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);

Solve the Problem

You generated all the solver inputs. Call the solver to find the solution.

[solution,fval,exitflag,output] = intlinprog(obj,intcon,．..Aineq,bineq,Aeq,beq,lb,ub,opts);

ifisempty(solution)% If the problem is infeasible or you stopped early with no solutiondisp('intlinprog did not return a solution.')return% Stop the script because there is nothing to examineend

Examine the Solution

The solution is feasible, to within the given tolerances.

exitflag

exitflag= 1

infeas1 = max(Aineq*solution - bineq)

infeas1 = 2.5011e-12

infeas2 = norm(Aeq*solution - beq,Inf)

infeas2 = 1.7621e-12

Check that the integer components are really integers, or are close enough that it is reasonable to round them. To understand why these variables might not be exactly integers, seeSome “Integer” Solutions Are Not Integers．

diffint = norm(solution(intcon) - round(solution(intcon)),Inf)

diffint = 1.4433e-15

Some integer variables are not exactly integers, but all are very close. So round the integer variables.

solution(intcon) = round(solution(intcon));

Check the feasibility of the rounded solution, and the change in objective function value.

infeas1 = max(Aineq*solution - bineq)

infeas1 = 2.5011e-12

infeas2 = norm(Aeq*solution - beq,Inf)

infeas2 = 1.7621e-12

diffrounding = norm(fval - obj(:)'*solution,Inf)

diffrounding = 2.2352e-08

Rounding the solution did not appreciably change its feasibility.

You can examine the solution most easily by reshaping it back to its original dimensions.

solution1 = solution(1:P*F*W);% The continuous variablessolution2 = solution(intcon);% The integer variablessolution1 = reshape(solution1,P,F,W); solution2 = reshape(solution2,S,W);

For example, 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(solution2,1)% Sum over the sales outlets

outlets =1×203 0 3 2 2 2 3 2 3 1 1 0 0 3 4 3 2 3 2 1

Plot the connection between each sales outlet and its warehouse.

figure(h); holdonforii = 1:S jj = find(solution2(ii,:));% Index of warehouse associated with iixsales = xloc(F+W+ii); ysales = yloc(F+W+ii); xwarehouse = xloc(F+jj); ywarehouse = yloc(F+jj);ifrand(1) < .5% Draw y direction first half the timeplot([xsales,xsales,xwarehouse],[ysales,ywarehouse,ywarehouse],'g--')else% Draw x direction first the rest of the timeplot([xsales,xwarehouse,xwarehouse],[ysales,ysales,ywarehouse],'g--')endendholdofftitle('Mapping of sales outlets to warehouses')

The black * with no green lines represent the unused warehouses.