rtetc.m

物流分析工具包。Facility location: Continuous minisum facility location, alternate location-allocation (ALA)

   
   if ~isempty(ld) && all(length(ld) ~= [1 n n+1])
      error('Length of "ld" must equal 1, n, or n + 1.')
   elseif ~isempty(TW)
      if iscell(TW), TW = cell2padmat(TW); end
      if all(size(TW,1) ~= [n n+1]) || mod(size(TW,2),2) ~= 0 || ...
            any(all(isnan(TW'))) || ...
            any(any(TW(:,1:2:end-1) > TW(:,2:2:end))) || ...
            any(any(xor(isnan(TW(:,1:2:end-1)),isnan(TW(:,2:2:end)))))
         error('TW not valid time windows.')
      end
   elseif ~isempty(st) && ~isempty(rte) && length(st(:)) ~= m
      error('Starting time "st" must be an m-element vector.')
   end
end

% Route feasibilty function
if ~isempty(rtefeas) &&  strcmp(rtefeas{1},'maxTCfeas') && ...
      all(isinf(rtefeas{2}))
   rtefeas = [];  % maxTC = Inf => always feasible
end
if isfirstcall && ~isempty(rtefeas)
   if length(rtefeas(:)) < 1
      error('"rtefeas" must be at least a one element cell array.')
   end
   if ~ischar(rtefeas{1})
      error('First element of "rtefeas" must be a string.')
   elseif ~strcmp(rtefeas{1},'maxTCfeas') && ~exist(rtefeas{1},'file')
      error(['Function "' rtefeas{1} '" not found.'])
   end
end

% Empty "rte" used for error checking and to store input arguments
if isempty(rte)
   if nargout > 0, TC = []; XFlg = []; out = []; end
   return
end
% End (Input Error Checking) **********************************************

% Initial timing output structure
if nargout > 2 || ~isempty(rtefeas)
   out = struct('Rte',[],'Cost',[],'Demand',[],'Arrive',[],'Wait',[],...
      'Start',[],'LD',[],'Depart',[],'Total',[],'EarlySF',[],'LateSF',[]);
   out(1:m) = out;
end

% Evaluate each route
for i = 1:m
    
   if iscell(rte), r = rte{i}; else r = rte(:)'; end
   
   % 1. Capacity feasibility
   if ~isempty(q) && ~isinf(Q)
      if nargout > 2 || ~isempty(rtefeas), out(i).Demand = q(r); end
      if sum(q(r)) > Q; XFlg(i) = -1; continue, end
   end
   
   % Calculate route cost
   c = diag(C(r(1:end-1),r(2:end)))';
   if isempty(ld) || all(ld == 0)
      ldi = zeros(1,length(r));
   else
      if length(ld) == 1
         ld = [0 ld*ones(1,n-1)];
         if r(1) == r(end), ld = [ld 0]; end
      end
      ldi = ld(r);
      if r(1) == r(end)
         if length(ld) ~= n+1
            error('Length of "ld" must equal n + 1.')
         end
         ldi(end) = ld(n+1);
      end
   end
   TC(i) = sum(c) + sum(ldi);
   if nargout > 2 || ~isempty(rtefeas)
      out(i).Rte = r;
      out(i).Cost = [0 c];
      if ~isempty(ld), out(i).LD = ldi; end
      out(i).Total = [0 c] + ldi;
   end
   
   % 2. Time window feasibility
   if ~isempty(TW)
      B = TW(:,1:2:end-1);
      E = TW(:,2:2:end);
      Br = B(r,:); Er = E(r,:);
      if r(1) == r(end)
         if size(B,1) ~= n+1, error('Length of TW must equal n + 1.'), end
         Br(end,:) = B(n+1,:); Er(end,:) = E(n+1,:);
      end
      if ~isempty(st), sti = st(i); else sti = []; end
      if nargout < 3 && isempty(rtefeas)
         TC(i) = rteTW(c,ldi,Br,Er,sti);
      else
         [TC(i),s,w] = rteTW(c,ldi,Br,Er,sti);
         
         if ~isinf(TC(i))
            s_late = latestart(c,ldi,Br,Er);
         else
            s_late = NaN;
         end        
         out(i).Arrive = [0 s(2:end)-w(2:end)];
         out(i).Wait = w;
         out(i).Start = s;
         out(i).Depart = s + ldi;
         out(i).Total = [0 c] + w + ldi;
         out(i).EarlySF = [s(1) s(end) + ldi(end)];
         out(i).LateSF = [s_late(1) s_late(1) + TC(i)];
      end
      if isinf(TC(i)), XFlg(i) = -2; continue, end
   end
   
   % 3. User defined feasibility function
   if ~isempty(rtefeas)
      isfeas = feval(rtefeas{1},out(i),rtefeas{2:end});
      if length(isfeas(:)) ~= 1 % | ~islogical(isfeas)
         error('Output argument "isfeas" must be a scalar logical value.')
      end
      if isfeas == 0, 
         TC(i) = Inf; XFlg(i) = -3; continue, end
   end
   
end % FOR loop


% *************************************************************************
% *************************************************************************
% *************************************************************************
function [TC,s,w] = rteTW(t,ld,B,E,st)
%RTETW Single route time window.

tol = 1e-8;
n = size(B,1);

if isempty(st), s = min(B(1,:)); else s = st; end

s = s + ld(1);
for i = 2:n  % Forward scan to determine earliest finish time
   s = s + t(i-1) + ld(i); 
   Bi = B(i,:) + ld(i);
   if ~any(s + tol >= Bi && s - tol <= E(i,:))
      s = min(Bi(Bi >= s));
      if isempty(s)
         TC = Inf; s = NaN; w = NaN; return
      end
   end
end
f = s;

s = f - ld(n);
for i = n-1:-1:1  % Reverse scan to determine latest start time for the
   % earliest finish
   s = s - t(i) - ld(i);
   Ei = E(i,:) - ld(i);
   if ~any(s + tol >= B(i,:) & s - tol <= Ei)
      s = max(Ei(Ei <= s));
   end
end
TC = f - s;
if isnan(TC), TC = sum(t) + sum(ld); end % If all Br == -Inf and all Er = Inf

if nargout > 1
   s = [s zeros(1,n-1)];
   w = zeros(1,n);
   for i = 2:n  % Second forward scan to delay waits as much as possible
      % to the end of the route in case unexpected events occur
      s(i) = s(i-1) + ld(i-1) + t(i-1); 
      Bi = B(i,:);
      if ~any(s(i) + tol >= Bi & s(i) + ld(i) - tol <= E(i,:))
         w(i) = s(i);
         s(i) = min(Bi(Bi >= s(i)));
         w(i) = s(i) - w(i);
      end
   end
end


% *************************************************************************
% *************************************************************************
% *************************************************************************
function s = latestart(t,ld,B,E)
%LATESTART Determine latest start time.

tol = 1e-8;
n = size(B,1);

s = max(E(end,:)) - ld(n);
for i = n-1:-1:1  % Reverse scan to determine latest start time
   s = s - t(i) - ld(i);
   Ei = E(i,:) - ld(i);
   if ~any(s + tol >= B(i,:) && s - tol <= Ei)
      s = max(Ei(Ei <= s));
   end
end


% *************************************************************************
% *************************************************************************
% *************************************************************************
function isfeas = maxTCfeas(outi,maxTC)
%MAXTCFEAS Maximum total cost route feasibility function.
% isfeas = maxwaitfeas(outi,maxTC)
%   outi = struct array of outputs from RTETC for single route i
%          (automatically passed to function)
%  maxTC = scalar maximum total cost (including un/loading times) of route
%
% Route is feasible if sum(outi.Total) <= maxTC
%
% This function can be used as a template for developing other
% route feasibility functions.

% Input error check
if ~isnumeric(maxTC) || length(maxTC(:)) ~= 1 || maxTC < 0
   error('"maxTC" must be a nonnegative scalar.')
end

% Feasibility test
if sum(outi.Total) <= maxTC
   isfeas = true;
else
   isfeas = false;
end