name = 'generationcost.xlsx';
poolPrice = xlsread(name,'AD3:AD50'); % Revenue in dollars per MWh in interval
demand = xlsread(name,'AC3:AC52');
month = numel(demand);
e = xlsread(name,'W2:W4')'; %Cost for a unit of fuel
fuelPrice = repmat(e,1,1,2);
nPeriods = length(poolPrice);
nGens = 3;
genhigh = xlsread(name,'F2:F4');
genlow = xlsread(name,'Q2:Q4');
gen = [genlow(:),genhigh(:),]; %MW

fuel = [50,330;45,325;55,750]; % Fuel consumption for generator
startCost = [10000;10000;10000]; %Cost in dollars to start a generator after it has been off

y = optimvar('y',nPeriods,nGens,{'Low','High'},'Type','integer','LowerBound',0,'UpperBound',1);
z = optimvar('z',nPeriods,nGens,'Type','integer','LowerBound',0,'UpperBound',1);

powercons = y(:,:,'Low') + y(:,:,'High') <= 1;


yFuel = zeros(nPeriods,nGens,2);
yFuel(:,1,1) = fuel(1,1);
yFuel(:,1,2) = fuel(1,2);
yFuel(:,2,1) = fuel(2,1);
yFuel(:,2,2) = fuel(2,2);
yFuel(:,3,1) = fuel(3,1);
yFuel(:,3,2) = fuel(3,2);

fuelUsed = sum(reshape((sum(y.*yFuel)),[3 2])');
%fuelcons = fuelUsed <= totalFuel;

w = optimexpr(nPeriods,nGens);
idx = 1:(nPeriods-1);
w(idx,:) = y(idx+1,:,'Low') - y(idx,:,'Low') + y(idx+1,:,'High') -
y(idx,:,'High');
w(nPeriods,:) = y(1,:,'Low') - y(nPeriods,:,'Low') + y(1,:,'High') -
y(nPeriods,:,'High');
switchcons = w - z <= 0;

generatorlevel  = zeros(size(yFuel));
generatorlevel(:,1,1) = gen(1,1);
generatorlevel(:,1,2) = gen(1,2);
generatorlevel(:,2,1) = gen(2,1);
generatorlevel(:,2,2) = gen(2,2);
generatorlevel(:,3,1) = gen(3,1);
generatorlevel(:,3,2) = gen(3,2);

revenue = optimexpr(size(y));
for ii = 1:nPeriods
    revenue(ii,:,:) = poolPrice(ii)*y(ii,:,:).*generatorlevel(ii,:,:);
end

val = optimexpr(size(y));
demands = optimexpr(nPeriods,1);
id = 1:nPeriods;
val(id,:,:)= y(id,:,:).*generatorlevel(id,:,:);
demands(id,:)= sum(sum(val,3),2);
demandcons = demands-demand ;

fuelCost = sum(sum(fuelPrice.*sum(y.*yFuel)));
startingCost = z*startCost;

profit = sum(sum(sum(revenue))) - fuelCost - sum(startingCost);

dispatch = optimproblem('ObjectiveSense','maximize');
dispatch.Objective = profit;

dispatch.Constraints.powercons = powercons;
dispatch.Constraints.demandcons = demands - demand <= 0 ;

options = optimoptions('intlinprog','Display','final');
[dispatchsol,fval,exitflag,output] = solve(dispatch,'options',options);

subplot(5,1,1)
bar(dispatchsol.y(:,1,1)*gen(1,1)+dispatchsol.y(:,1,2)*gen(1,2),.5,'g')
xlim([.5,48.5])
ylabel('MW')
title('Generator 1 Optimal Schedule','FontWeight','bold')
subplot(5,1,2)
bar(dispatchsol.y(:,2,1)*gen(2,1)+dispatchsol.y(:,2,2)*gen(2,2),.5,'c')
title('Generator 2 Optimal Schedule','FontWeight','bold')
xlim([.5,48.5])
ylabel('MW')
subplot(5,1,3)
bar(dispatchsol.y(:,3,1)*gen(3,1)+dispatchsol.y(:,3,2)*gen(3,2),.5,'c')
title('Generator 3 Optimal Schedule','FontWeight','bold')
xlim([.5,48.5])
ylabel('MW')
subplot(5,1,4)
bar(demand,.5)
xlim([.5,48.5])
title('Daily Demand','FontWeight','bold')
xlabel('Period')
ylabel('MW')

starttimes = find(round(dispatchsol.z) == 1);
[theperiod,thegenerator] = ind2sub(size(dispatchsol.z),starttimes)