% 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))