www.gusucode.com > optim 案例源码 matlab代码程序 > optim/OptimalDispatchExample.m

    %% Optimal Dispatch of Power Generators
% 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.

%   Copyright 2014 The MathWorks, Inc. 

%% 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
% |http://www.nemweb.com.au/REPORTS/CURRENT/| in mid-2013, and is used
% under their terms
% |http://www.aemo.com.au/About-AEMO/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 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)