www.gusucode.com > optim 案例源码 matlab代码程序 > optim/FactoryExample.m
%% Factory, Warehouse, Sales Allocation Model % 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. % % 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$. % Copyright 2014-2016 The MathWorks, Inc. %% 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 Objective and Constraint Matrices and Vectors % The objective function vector |obj| in |intlincon| consists of the % coefficients of the variables $x(p,f,w)$ and $y(s,w)$. So there are % naturally |P*F*W + S*W| coefficients in |obj|. % % One way to generate the coefficients is to begin with a |P-by-F-by-W| % array |obj1| for the $x$ coefficients, and an |S-by-W| array |obj2| for % the $y(s,w)$ coefficients. Then convert these arrays to two vectors and % combine them into |obj| by calling % % |obj = [obj1(:);obj2(:)];| obj1 = zeros(P,F,W); % Allocate arrays obj2 = 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 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 %% % Generate the entries of |obj1| and |obj2|. for ii = 1:P for jj = 1:F for kk = 1:W obj1(ii,jj,kk) = pcost(jj,ii) + tcost(ii)*distfw(jj,kk); end end end for ii = 1:S for jj = 1:W obj2(ii,jj) = distsw(ii,jj)*sum(d(ii,:).*tcost); end end %% % 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 the |obj| % vector. matwid = length(obj); %% % There are two types of linear inequalities: the production capacity % constraints, and the warehouse capacity constraints. % % There are |P*F| production capacity constraints, and |W| warehouse % 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 Aeq bineq = 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 constraints counter = 1; for ii = 1:F for jj = 1:P xtemp = clearer1; xtemp(jj,ii,:) = 1; % Sum over warehouses for each product and factory xtemp = sparse([xtemp(:);clearer22]); % Convert to sparse Aineq(counter,:) = xtemp'; % Fill in the row bineq(counter) = pcap(ii,jj); counter = counter + 1; end end % Now the warehouse capacity constraints vj = zeros(S,1); % The multipliers for jj = 1:S vj(jj) = sum(d(jj,:)./turn); % A sum of P elements end for ii = 1:W xtemp = clearer2; xtemp(:,ii) = vj; xtemp = sparse([clearer12;xtemp(:)]); % Convert to sparse Aineq(counter,:) = xtemp'; % Fill in the row bineq(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 sparse beq = zeros(P*W + S,1); % Allocate vectors as full counter = 1; % Demand is satisfied: for ii = 1:P for jj = 1:W xtemp = clearer1; xtemp(ii,:,jj) = 1; xtemp2 = clearer2; xtemp2(:,jj) = -d(:,ii); xtemp = sparse([xtemp(:);xtemp2(:)]'); % Change to sparse row Aeq(counter,:) = xtemp; % Fill in row counter = counter + 1; end end % Only one warehouse for each sales outlet: for ii = 1:S xtemp = clearer2; xtemp(ii,:) = 1; xtemp = sparse([clearer12;xtemp(:)]'); % Change to sparse row Aeq(counter,:) = xtemp; % Fill in row beq(counter) = 1; counter = counter + 1; end %% Bound Constraints and Integer Variables % The integer variables are those from |length(obj1) + 1| to the end. intcon = P*F*W+1:length(obj); %% % The upper bounds are from |length(obj1) + 1| to 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); if isempty(solution) % If the problem is infeasible or you stopped early with no solution disp('intlinprog did not return a solution.') return % Stop the script because there is nothing to examine end %% Examine the Solution % The solution is feasible, to within the given tolerances. exitflag infeas1 = max(Aineq*solution - bineq) infeas2 = norm(Aeq*solution - beq,Inf) %% % 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, see % <matlab:helpview(fullfile(docroot,'toolbox','optim','helptargets.map'),'milp_not_integer'); % the documentation>. diffint = norm(solution(intcon) - round(solution(intcon)),Inf) %% % 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) infeas2 = norm(Aeq*solution - beq,Inf) diffrounding = norm(fval - obj(:)'*solution,Inf) %% % 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 variables solution2 = solution(intcon); % The integer variables solution1 = 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 %% % Plot the connection between each sales outlet and its warehouse. figure(h); hold on for ii = 1:S jj = find(solution2(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.