%% Optimal Dispatch of Power Generators: Solver-Based
% This example shows how to schedule two gas-fired electric generators optimally, meaning to get the most revenue minus cost. While the example is not entirely realistic, it does show how to take into account costs that depend on decision timing.
% For the problem-based approach to this problem, see Optimal Dispatch of Power Generators: Problem-Based.
%% Problem Definition
% The electricity market has different prices at different times of day. If you have generators, you can take advantage of this variable pricing by scheduling your generators to operate when prices are high. Suppose that there are two generators that you control. Each generator has three power levels (off, low, and high). Each generator has a specified rate of fuel consumption and power production at each power level. Of course, fuel consumption is 0 when the generator is off.
% You can assign a power level to each generator during each half-hour time interval during a day (24 hours, so 48 intervals). Based on historical records, you can assume that you know the revenue per megawatt-hour (MWh) that you get in each time interval. The data for this example is from the Australian Energy Market Operator https://www.nemweb.com.au/REPORTS/CURRENT/ in mid-2013, and is used under their terms https://www.aemo.com.au/privacy-and-legal-notices/copyright-permissions.

load dispatchPrice; % Get poolPrice, which is the revenue per MWh
bar(poolPrice,.5)
xlim([.5,48.5])
xlabel('Price per MWh at each period')

% There is a cost to start a generator after it has been off. The other constraint is a maximum fuel usage for the day. The maximum fuel constraint is because you buy your fuel a day ahead of time, so can use only what you just bought.
%
%% Problem Notation and Parameters
% You can formulate the scheduling problem as a binary integer programming problem as follows. Define indexes i, j, and k, and a binary scheduling vector y as:
% nPeriods = the number of time periods, 48 in this case.
% i = a time period, 1 <= i <= 48.
% j = a generator index, 1 <= j <= 2 for this example.
% y(i,j,k) = 1 when period i, generator j is operating at power level k. Let low power be k = 1, and high power be k = 2. The generator is off when sum_k y(i,j,k) = 0.
% You need to determine when a generator starts after being off. Let
% z(i,j) = 1 when generator j is off at period i, but is on at period i + 1. z(i,j) = 0 otherwise. In other words, z(i,j) = 1 when sum_k y(i,j,k) = 0 and sum_k y(i+1,j,k) = 1.
% Obviously, you need a way to set z automatically based on the settings of y. A linear constraint below handles this setting.
% You also need the parameters of the problem for costs, generation levels for each generator, consumption levels of the generators, and fuel available.
% poolPrice(i) -- Revenue in dollars per MWh in interval i.
% gen(j,k) -- MW generated by generator j at power level k.
% fuel(j,k) -- Fuel used by generator j at power level k.
% totalfuel -- Fuel available in one day.
% startCost -- Cost in dollars to start a generator after it has been off.
% fuelPrice -- Cost for a unit of fuel.
% You got poolPrice when you executed load dispatchPrice;. Set the other parameters as follows.

fuelPrice = 3;
totalfuel = 3.95e4;
nPeriods = length(poolPrice); % 48 periods
nGens = 2; % Two generators
gen = [61,152;50,150]; % Generator 1 low = 61 MW, high = 152 MW
fuel = [427,806;325,765]; % Fuel consumption for generator 2 is low = 325, high = 765
startCost = 1e4; % Cost to start a generator after it has been off

%
%% Generator Efficiency
% Examine the efficiency of the two generators at their two operating points.

efficiency = gen./fuel; % Calculate electricity per unit fuel use
rr = efficiency'; % for plotting
h = bar(rr);
h(1).FaceColor = 'g';
h(2).FaceColor = 'c';
legend(h,'Generator 1','Generator 2','Location','NorthEastOutside')
ax = gca;
ax.XTick = [1,2];
ax.XTickLabel = {'Low','High'};
ylim([.1,.2])
ylabel('Efficiency')

% Notice that generator 2 is a bit more efficient than generator 1 at its corresponding operating points (low or high), but generator 1 at its high operating point is more efficient than generator 2 at its low operating point.
%% Variables for Solution
% To set up the problem, you need to encode all the problem data and constraints in the form that the intlinprog solver requires. You have variables y(i,j,k) that represent the solution of the problem, and z(i,j) auxiliary variables for charging to turn on a generator. y is an nPeriods-by-nGens-by-2 array, and z is an nPeriods-by-nGens array.
% To put these variables in one long vector, define the variable of unknowns x:
% x = [y(:);z(:)];
% For bounds and linear constraints, it is easiest to use the natural array formulation of y and z, then convert the constraints to the total decision variable, the vector x.
%% Bounds
% The solution vector x consists of binary variables. Set up the bounds lb and ub.

lby = zeros(nPeriods,nGens,2); % 0 for the y variables
lbz = zeros(nPeriods,nGens); % 0 for the z variables
lb = [lby(:);lbz(:)]; % Column vector lower bound
ub = ones(size(lb)); % Binary variables have lower bound 0, upper bound 1

%
%% Linear Constraints
% For linear constraints A*x <= b, the number of columns in the A matrix must be the same as the length of x, which is the same as the length of lb. To create rows of A of the appropriate size, create zero matrices of the sizes of the y and z matrices.

cleary = zeros(nPeriods,nGens,2);
clearz = zeros(nPeriods,nGens);

%
% To ensure that the power level has no more than one component equal to 1, set a linear inequality constraint:
% x(i,j,1) + x(i,j,2) <= 1

A = spalloc(nPeriods*nGens,length(lb),2*nPeriods*nGens); % nPeriods*nGens inequalities
counter = 1;
for ii = 1:nPeriods
    for jj = 1:nGens
        temp = cleary;
        temp(ii,jj,:) = 1;
        addrow = [temp(:);clearz(:)]';
        A(counter,:) = sparse(addrow);
        counter = counter + 1;
    end
end
b = ones(nPeriods*nGens,1); % A*x <= b means no more than one of x(i,j,1) and x(i,j,2) are equal to 1

%
% The running cost per period is the cost for fuel for that period. For generator j operating at level k, the cost is fuelPrice * fuel(j,k).
% To ensure that the generators do not use too much fuel, create an inequality constraint on the sum of fuel usage.

yFuel = lby; % Initialize fuel usage array
yFuel(:,1,1) = fuel(1,1); % Fuel use of generator 1 in low setting
yFuel(:,1,2) = fuel(1,2); % Fuel use of generator 1 in high setting
yFuel(:,2,1) = fuel(2,1); % Fuel use of generator 2 in low setting
yFuel(:,2,2) = fuel(2,2); % Fuel use of generator 2 in high setting

addrow = [yFuel(:);clearz(:)]';
A = [A;sparse(addrow)];
b = [b;totalfuel]; % A*x <= b means the total fuel usage is <= totalfuel

%
%% Set the Generator Startup Indicator Variables
% How can you get the solver to set the z variables automatically to match the active/off periods that the y variables represent? Recall that the condition to satisfy is z(i,j) = 1 exactly when
% sum_k y(i,j,k) = 0 and sum_k y(i+1,j,k) = 1.
% Notice that
% sum_k ( - y(i,j,k) + y(i+1,j,k) ) > 0 exactly when you want z(i,j) = 1.
% Therefore, include the linear inequality constraints
% sum_k ( - y(i,j,k) + y(i+1,j,k) ) - z(i,j) < = 0
% in the problem formulation, and include the z variables in the objective function cost. By including the z variables in the objective function, the solver attempts to lower the values of the z variables, meaning it tries to set them all equal to 0. But for those intervals when a generator turns on, the linear inequality forces the z(i,j) to equal 1.
% Add extra rows to the linear inequality constraint matrix A to represent these new inequalities. Wrap around the time so that interval 1 logically follows interval 48.

tempA = spalloc(nPeriods*nGens,length(lb),2*nPeriods*nGens);
counter = 1;
for ii = 1:nPeriods
    for jj = 1:nGens
        temp = cleary;
        tempy = clearz;
        temp(ii,jj,1) = -1;
        temp(ii,jj,2) = -1;
        if ii < nPeriods % Intervals 1 to 47
            temp(ii+1,jj,1) = 1;
            temp(ii+1,jj,2) = 1;
        else % Interval 1 follows interval 48
            temp(1,jj,1) = 1;
            temp(1,jj,2) = 1;
        end
        tempy(ii,jj) = -1;
        temp = [temp(:);tempy(:)]'; % Row vector for inclusion in tempA matrix
        tempA(counter,:) = sparse(temp);
        counter = counter + 1;
    end
end
A = [A;tempA];
b = [b;zeros(nPeriods*nGens,1)]; % A*x <= b sets z(i,j) = 1 at generator startup

%
%% Sparsity of Constraints
% If you have a large problem, using sparse constraint matrices saves memory, and can save computational time as well. The constraint matrix A is quite sparse:

filledfraction = nnz(A)/numel(A)

% intlinprog accepts sparse linear constraint matrices A and Aeq, but requires their corresponding vector constraints b and beq to be full.
%
%% Define Objective
% The objective function includes fuel costs for running the generators, revenue from running the generators, and costs for starting the generators.

generatorlevel = lby; % Generation in MW, start with 0s
generatorlevel(:,1,1) = gen(1,1); % Fill in the levels
generatorlevel(:,1,2) = gen(1,2);
generatorlevel(:,2,1) = gen(2,1);
generatorlevel(:,2,2) = gen(2,2);

%
% Incoming revenue = x.*generatorlevel.*poolPrice

revenue = generatorlevel; % Allocate revenue array
for ii = 1:nPeriods
    revenue(ii,:,:) = poolPrice(ii)*generatorlevel(ii,:,:);
end

%
% Total fuel cost = y.*yFuel*fuelPrice

fuelCost = yFuel*fuelPrice;

%
% Startup cost = z.*ones(size(z))*startCost

starts = (clearz + 1)*startCost;
starts = starts(:); % Generator startup cost vector

%
% The vector x = [y(:);z(:)]. Write the total profit in terms of x:
% profit = Incoming revenue - Total fuel cost - Startup cost

f = [revenue(:) - fuelCost(:);-starts]; % f is the objective function vector

%
%% Solve the Problem
% To save space, suppress iterative display.

options = optimoptions('intlinprog','Display','final');
[x,fval,eflag,output] = intlinprog(-f,1:length(f),A,b,[],[],lb,ub,options);

%
%% Examine the Solution
% The easiest way to examine the solution is dividing the solution vector x into its two components, y and z.

ysolution = x(1:nPeriods*nGens*2);
zsolution = x(nPeriods*nGens*2+1:end);
ysolution = reshape(ysolution,[nPeriods,nGens,2]);
zsolution = reshape(zsolution,[nPeriods,nGens]);

%
% Plot the solution as a function of time.

subplot(3,1,1)
bar(ysolution(:,1,1)*gen(1,1)+ysolution(:,1,2)*gen(1,2),.5,'g')
xlim([.5,48.5])
ylabel('MWh')
title('Generator 1 optimal schedule','FontWeight','bold')
subplot(3,1,2)
bar(ysolution(:,2,1)*gen(1,1)+ysolution(:,2,2)*gen(1,2),.5,'c')
title('Generator 2 optimal schedule','FontWeight','bold')
xlim([.5,48.5])
ylabel('MWh')
subplot(3,1,3)
bar(poolPrice,.5)
xlim([.5,48.5])
title('Energy price','FontWeight','bold')
xlabel('Period')
ylabel('$ / MWh')

%
% Generator 2 runs longer than generator 1, which you would expect because it is more efficient. Generator 2 runs at its high power level whenever it is on. Generator 1 runs mainly at its high power level, but dips down to low power for one time unit. Each generator runs for one contiguous set of periods daily, so incurs only one startup cost.
% Check that the z variable is 1 for the periods when the generators start.

starttimes = find(round(zsolution) == 1); % Use round for noninteger results
[theperiod,thegenerator] = ind2sub(size(zsolution),starttimes)

% The periods when the generators start match the plots.
%
%% Compare to Lower Penalty for Startup
% If you choose a small value of startCost, the solution involves multiple generation periods.

startCost = 500; % Choose a lower penalty for starting the generators
starts = (clearz + 1)*startCost;
starts = starts(:); % Start cost vector
fnew = [revenue(:) - fuelCost(:);-starts]; % New objective function
[xnew,fvalnew,eflagnew,outputnew] = ...
    intlinprog(-fnew,1:length(fnew),A,b,[],[],lb,ub,options);

ysolutionnew = xnew(1:nPeriods*nGens*2);
zsolutionnew = xnew(nPeriods*nGens*2+1:end);
ysolutionnew = reshape(ysolutionnew,[nPeriods,nGens,2]);
zsolutionnew = reshape(zsolutionnew,[nPeriods,nGens]);

subplot(3,1,1)
bar(ysolutionnew(:,1,1)*gen(1,1)+ysolutionnew(:,1,2)*gen(1,2),.5,'g')
xlim([.5,48.5])
ylabel('MWh')
title('Generator 1 optimal schedule','FontWeight','bold')
subplot(3,1,2)
bar(ysolutionnew(:,2,1)*gen(1,1)+ysolutionnew(:,2,2)*gen(1,2),.5,'c')
title('Generator 2 optimal schedule','FontWeight','bold')
xlim([.5,48.5])
ylabel('MWh')
subplot(3,1,3)
bar(poolPrice,.5)
xlim([.5,48.5])
title('Energy price','FontWeight','bold')
xlabel('Period')
ylabel('$ / MWh')

starttimes = find(round(zsolutionnew) == 1); % Use round for noninteger results
[theperiod,thegenerator] = ind2sub(size(zsolution),starttimes)

% Copyright 2014-2017 The MathWorks, Inc.