外国人写的,matlab官方网站下载的。程序里面计算成本矩阵用到的矩阵变换很亮,虽然一下子很难看明白。不过知道方法的原理后,直接调用函数计算就可以了
function [z , x , J , arg_J , C] = wargner_whitin(d , K , c , h);
%
% Wagner-Whitin Algorithm
%
% Build production plans versus time cycles k =1,...,N responding production queries dk
% and miniming the sum of following production costs
%
% i) Fix production cost K_k
% ii) Unitary production cost c_k
% iii) Stockage unitary cost h_k
%
% Assumption : initial and final stocks are null.
%
%
% Inputs
% ------
%
% d Queries of production througth time cycles (1 x N)
% K Fix production cost (1 x N)
% c Unitary production cost (1 x N)
% h Stockage unitary cost (1 x N)
%
%
% Outputs
% ------
%
%
% z Optimal cost at the end of the N cycles
% x Quantities to product (1 x N)
% J Partial politics (1 x N)
% arg_J Arguments of the partial politics (1 x N)
% C Cost transition matrix (N x N)
%
%
%
% Example1
% -------
%
%clear,close all hidden
%d = [8 10 7 9 6];
%K = [50 40 50 70 60];
%c = [5 5 6 4 5];
%h = [4 2 3 4 NaN];
%[z , x , J , arg_J , C] = wargner_whitin(d , K , c , h);
%
%
% Example2
% --------
%
%
%clear,close all hidden
%N = 100;
%d = ceil(N*rand(1 , N));
%K = ceil(N*rand(1 , N));
%c = ceil(N*rand(1 , N));
%h = ceil(N*rand(1 , N));
%[z , x , J , arg_J , C] = wargner_whitin(d , K , c , h);
%ind_x = find(x);
%figure(1)
%subplot(211)
%stem((1:N),J)
%xlabel('Period N','fontname','times','fontsize',12)
%ylabel('Production costs','fontname','times','fontsize',12)
%title(['N = ' num2str(N) ],'fontname','times','fontsize',13)
%subplot(212)
%stem(ind_x , x(ind_x) , 'ro')
%xlabel('P-N','fontname','times','fontsize',12)
%ylabel('Quantities to product','fontname','times','fontsize',12)
%figure(2)
%imagesc(C)
%colorbar
%title('Transition cost production matrix ','fontname','times','fontsize',13)
%
%
% Author : Sebastien PARIS (sebastien.paris@lsis.org) 08/30/2001.
% ------
%
N = length(d);
%------------------- Initialization ----------------------%
C = zeros(N);
FF = C;
J = zeros(1 , N + 1);
arg_J = J ;
ZN = zeros(1 , N);
ON = ones(1 , N);
vect_N = (1:N);
%------- Cost matrix C = {c(k,l)}, k,l = 1,...,N -----%
mat_tp = triu(vect_N(ON , :));
K1 = K(:);
c1 = c(:);
h1 = h(:);
h1 = h1(1:end - 1);
KK = triu(K1(: , ON));
CC = triu(c1(: , ON));
indice = (mat_tp ~= 0);
FF(indice) = d(mat_tp(indice));
YY = cumsum(FF , 2);
SS = YY(2:end , 1:end);
HH1 = triu(h1(: , ON));
EE = HH1.*SS;
ZZ = cumsum(EE(end : -1 : 1 , :) , 1);
ZZ = ZZ(end : -1 : 1 , :);
C = KK + (CC.*YY) + ([ZZ ; ZN]);
%------------------------ Dynamic Programming --------------------%
J(end) = 0;
arg_J(end) = 0;
for k = N : -1 : 1
tp = ZN;
tp(:) = NaN;
tp(k :end) = C(k , k:end ) + J(k + 1 : end);
[J(k) arg_J(k) ] = min(tp);
end
z = J(1);
%------------------------ Back-Tracking --------------------%
k = 1;
x = ZN;
while(k <= N)
x(k) = sum(d(k : arg_J(k)));
k = arg_J(k) + 1;
end
J = J(1:end - 1);