Документация

Проблемное определение

Рынок электроэнергии имеет различные цены в разное время дня. Если у вас есть генераторы, можно использовать в своих интересах эту переменную оценку путем планирования генераторов, чтобы действовать, когда цены высоки. Предположим, что существует два генератора, которыми вы управляете. Каждый генератор имеет три уровня мощности (прочь, низко, и высоко). Каждый генератор имеет заданный уровень расхода топлива и выработки энергии на каждом уровне мощности. Конечно, расход топлива 0, когда генератор выключен.

Можно присвоить уровень мощности каждому генератору во время каждого получасового временного интервала в течение дня (24 часа, таким образом, 48 интервалов). На основе хронологических записей можно принять, что вы знаете доход на мегаватт-час (МВт·ч), что вы входите в каждый временной интервал. Данные для этого примера являются от австралийского Оператора Энергетического рынка https://www.nemweb.com.au/REPORTS/CURRENT/ в середине 2013 и используются в соответствии с их условиями https://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')

Существует стоимость, чтобы запустить генератор после того, как это было выключено. Другое ограничение является максимальным топливным использованием в течение дня. Максимальное топливное ограничение состоит в том, потому что вы покупаете свое топливо день загодя, так может использовать только, что вы только купили.

Проблемное обозначение и параметры

Можно сформулировать проблему планирования как бинарную проблему целочисленного программирования можно следующим образом. Задайте индексирует i, j, и k и двоичный файл, планируя векторный y как:

Необходимо определить, когда генератор запускается будучи прочь. Пусть

Очевидно, вам нужен способ установить z автоматически на основе настроек y. Линейное ограничение ниже обрабатывает эту установку.

Вам также нужны параметры проблемы для затрат, уровней генерации для каждого генератора, уровней потребления генераторов и доступного топлива.

Вы получили poolPrice, когда вы выполнили load dispatchPrice;. Установите другие параметры можно следующим образом.

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

Эффективность генератора

Исследуйте эффективность этих двух генераторов в их двух рабочих точках.

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

Заметьте, что генератор 2 немного более эффективен, чем генератор 1 в его соответствующих рабочих точках (низко или высоко), но генератор 1 в его высокой рабочей точке более эффективен, чем генератор 2 в его низкой рабочей точке.

Переменные для решения

Чтобы настроить проблему, необходимо закодировать все проблемные данные и ограничения в форме, которой требует решатель intlinprog. У вас есть переменные y(i,j,k), которые представляют решение проблемы и z(i,j) вспомогательные переменные для зарядки, чтобы включить генератор. y является массивом nPeriods-by-nGens-by-2, и z является массивом nPeriods-by-nGens.

Чтобы поместить эти переменные в один длинный вектор, задайте переменную неизвестных x:

x = [y(:);z(:)];

Для границ и линейных ограничений, является самым легким использовать естественную формулировку массивов y и z, затем преобразовать ограничения в общую переменную решения, векторный x.

Границы

Вектор решения x состоит из бинарных переменных. Настройте границы lb и 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

Линейные ограничения

Для линейных ограничений A*x <= b количество столбцов в матрице A должно совпасть с длиной x, который совпадает с длиной lb. Чтобы создать строки A соответствующего размера, создайте нулевые матрицы размеров матрицы z и y.

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

Чтобы гарантировать, что уровень мощности имеет не больше, чем один компонент, равный 1, устанавливает линейное ограничение неравенства:

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

Производственные затраты на период являются стоимостью для топлива в течение того периода. Для генератора j, действующий на уровне k, стоимостью является fuelPrice * fuel(j,k).

Чтобы гарантировать, что генераторы не используют слишком большого количества топлива, создайте ограничение неравенства на сумму топливного использования.

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

Установите переменные индикатора запуска генератора

Как можно заставить решатель устанавливать переменные z автоматически совпадать с активными/от периодами, которые представляют переменные y? Вспомните, что условием удовлетворить является z(i,j) = 1 точно когда

sum_k y(i,j,k) = 0 и sum_k y(i+1,j,k) = 1.

Заметьте это

sum_k ( - y(i,j,k) + y(i+1,j,k) ) > 0 точно когда это необходимо z(i,j) = 1.

Поэтому включайте линейные ограничения неравенства

sum_k ( - y(i,j,k) + y(i+1,j,k) ) - z(i,j) < = 0

в проблемной формулировке, и включают переменные z в стоимость целевой функции. Включением переменных z в целевой функции решатель пытается понизить значения переменных z, означая, что это пытается установить их всех, равняются 0. Но для тех интервалов, когда генератор включает, линейное неравенство обеспечивает z(i,j), чтобы равняться 1.

Добавьте дополнительные строки в линейную матрицу ограничения неравенства A, чтобы представлять эти новые неравенства. Перенеситесь примерно во время то, так, чтобы интервал 1 логически следовал за интервалом 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

Разреженность ограничений

Если у вас есть большая проблема, использование разреженных ограничительных матриц сохраняет память и может сэкономить вычислительное время также. Матрица ограничений A довольно разреженна:

filledfraction = nnz(A)/numel(A)

filledfraction = 0.0155

intlinprog принимает разреженные линейные ограничительные матрицы A и Aeq, но требует, чтобы их соответствующие векторные ограничения b и beq были полны.

Задайте цель

Целевая функция включает затраты на топливо для выполнения генераторов, дохода от выполнения генераторов, и стоит за запуск генераторов.

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

Входящий доход = x.*generatorlevel.*poolPrice

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

Общая стоимость топлива = y.*yFuel*fuelPrice

fuelCost = yFuel*fuelPrice;

Стоимость запуска = z.*ones(size(z))*startCost

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

Векторный x = [y(:);z(:)]. Запишите общую прибыль с точки зрения x:

получите прибыль = Входящий доход - Общая стоимость топлива - стоимость Запуска

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

Решите проблему

Чтобы оставить свободное место, подавите итеративное отображение.

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

Optimal solution found.

Intlinprog stopped because the objective value is within a gap tolerance of the
optimal value, options.AbsoluteGapTolerance = 0 (the default value). The intcon
variables are integer within tolerance, options.IntegerTolerance = 1e-05 (the
default value).

Исследуйте решение

Самый легкий способ исследовать решение делит вектор решения x на его два компонента, y и 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]);

Постройте решение как функцию времени.

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

Генератор 2 выполнения дольше, чем генератор 1, который вы ожидали бы, потому что это более эффективно. Генератор 2 выполнения на его мощном уровне каждый раз, когда это включено. Генератор 1 выполнение в основном на его мощном уровне, но падениях вниз к малой мощности для одной единицы измерения времени. Каждые выполнения генератора для одного непрерывного набора периодов ежедневно, поэтому несет только одни расходы запуска.

Проверяйте, что переменная z 1 в течение периодов, когда генераторы запускаются.

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

theperiod = 2×1

    23
    16

thegenerator = 2×1

     1
     2

Периоды, когда генераторы начинают матч графики.

Сравните с более низким штрафом за запуск

Если вы выбираете маленькое значение startCost, решение включает несколько периодов генерации.

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

Optimal solution found.

Intlinprog stopped because the objective value is within a gap tolerance of the
optimal value, options.AbsoluteGapTolerance = 0 (the default value). The intcon
variables are integer within tolerance, options.IntegerTolerance = 1e-05 (the
default value).

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)

theperiod = 3×1

    22
    16
    45

thegenerator = 3×1

     1
     2
     2