% ROME Example: Robust Warehouse Management
%
% See: robust_warehouse.pdf
%
% Modified by:
% 1. Joel & Melvyn   (Created 16 Sep 2009)

% Make sure you have the latest version of ROME

% Display welcome message
clear;
disp('Robust Warehouse Management');

% Setup parameters
M = 3;                      % Number of products
N = 5;                      % Number of classes
T = 3;                     % Number of periods 

%Warehouse parameters
s = [1 1.5 2 3 100];       % Storage cost
r = [1 1.5 2 3 100];       % Retreve cost
c = [150 300 450 600 inf];     % Capacity
x0 = zeros(M,N);            % Initial number of items at storage class

% Arrival infromation
a = [350 350 350
     350 350 350
     100 100 250]; 
% Demand information
mu = [290 190 450
      100 470 370
      60  70  200];  % Average info
sigma = 50*ones(M,T);   % Demand deviation

    
% begin rome
h = rome_begin;   

% Define uncertain factors
newvar z(M,T) uncertain;
z.set_mean(0);
rome_box(z, -1, 1);
d = mu + sigma.*z; % Create demand uncertainty

% demand uncertainty set
rome_constraint(d>=0);
for tt = 1:T
      rome_constraint(sum(d(:,1:tt),2)<=sum(a(:,1:tt),2));
end

% num pallets of product i assigned to storage class j in period t 
% (after observing z_(t-1))
v = [];
for tt = 1:T
    newvar tmp(M,N,z(:, 1:tt-1)) linearrule nonneg
    v = cat(3,v,tmp);
end

% num pallets of product i retrieved from storage class j in period t (
% after observing z_t)
w = []; 
for tt = 1:T
    newvar tmp(M,N,z(:, 1:tt)) linearrule nonneg
    w = cat(3,w,tmp);
end

% num pallets of product i in storage class j at start of period t 
% (after observing z_t-1)
x = [];  
for tt = 1:T+1
    newvar tmp(M,N,z(:, 1:tt-1)) linearrule nonneg
    x = cat(3,x,tmp);
end


% No Initial Inventory Constraint Set
rome_constraint(x(:,:,1) == x0);  

% Define Constraints
% Goods Supply Constraint Set
rome_constraint(squeeze(sum(v, 2)) == a);

% Goods Demand Constraint Set
rome_constraint(squeeze(sum(w, 2)) == d);

% Goods Carry Forward Constraint Set
for tt =1:T
    rome_constraint(x(:,:,tt+1) == x(:,:,tt) + v(:,:,tt)-w(:,:,tt));
end

% Storage Capacity Constraint Set
cfinite = find(c<inf);  % Indices of finite capacity
c_ijt = repmat(c(cfinite)',1,T);
rome_constraint(squeeze(sum(x(:,cfinite,1:T) + v(:,cfinite,1:T),1)) ...
   <= c_ijt);  


% Formulate Objective
rome_minimize(mean(sum(s.*sum(sum(v,1),3)) + sum(r.*sum(sum(w,1),3))));
        

% solve
tic
h.solve;
toc
obj_val = h.objective;

% Display demand information
disp('Display demand');
h.eval(d)

vsol = linearpart(h.eval(v));
wsol = linearpart(h.eval(w));

% Display storage policy at t= 1
disp('Storage Policy at t=1');
vsol(:,:,1,1)   


% Display retrieval policy at t= 1
disp('Retrieval Policy at t=1');
wsol(:,:,1,1:4)   

% Display objective
disp('Objective value');
disp(obj_val);

rome_end;





% ROME: Copyright (C) 2009 by Joel Goh and Melvyn Sim
% See the file COPYING.txt for full copyright information.