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)