% Backward DDP to solve Capacity Expansion Problem

% define inputs
nstages=3
costfrom=zeros(22,22)
% discount rate
rate=0.08
% cost array
exptab=[0:0.5:2.5]
costtab=[0 20 30 36 40 42]
expansion=[0.0:0.1:2.1]
capacity=[0.7:0.1:2.8]
cost=interp1(exptab,costtab,expansion)
% set up time array
demandtab=[0.7 0.9 1.1 1.7 2.2 2.5 2.7 2.8]
timetab=[0:5:35]
tau=interp1(demandtab,timetab,capacity)
demand=interp1(timetab,demandtab,tau)
close all
figure
plot(tau, capacity)

% backward sweep
vopt=zeros(22,nstages+1);
for stage=nstages:-1:1
   level2=22;
   if(stage==1) level2=1; end
   for level=1:level2
      discount(level)=(1/(1+rate))^tau(level);
      nextlevel1=level;
      if(stage==nstages) nextlevel1=22; end
      for nextlevel=nextlevel1:22
         increment=nextlevel-level+1;
         costfrom(level,nextlevel)=discount(level)*cost(increment);
      end
      [vopt(level,stage),kopt]= ...
         min(costfrom(level,(nextlevel1:22))+vopt((nextlevel1:22),stage+1)');
      uopt(level,stage)=nextlevel1-1+kopt-level;
   end
end
vopt=vopt
uopt=uopt
costopt=vopt(1,1)
%forward sweep
levelf(1)=1;
for stage=1:nstages;
   construct(stage)=0.1*(uopt(levelf(stage),stage)');
   levelf(stage+1)=levelf(stage)+uopt(levelf(stage),stage);
end
construct=construct
optcapacity=(levelf+6)*0.1
optcapacity=[optcapacity  optcapacity(length(optcapacity))]
exptime=tau(levelf(1:nstages))
exptime=[exptime 35]
hold on
stairs(exptime,optcapacity(2:5))