Определение проблемы

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

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

Кроме того, существует ограничение на максимальное потребление топлива за день. Это ограничение существует, потому что вы покупаете топливо на день раньше времени, так что вы можете использовать только то, что вы только что купили.

Нотация и параметры проблемы

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

Определение момента запуска генератора после отключения. Для этого определите вспомогательную двоичную переменную z(i,j) указывает, следует ли взимать плату за включение генератора j в период i.

Вам нужен способ установить 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 в своей нижней рабочей точке.

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

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

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

Линейные зависимости

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

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

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

Создание выражения fuelUsed это учитывает все использованное топливо.

yFuel = zeros(nPeriods,nGens,2);
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

fuelUsed = sum(sum(sum(y.*yFuel)));

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

fuelcons = fuelUsed <= 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 переменные в целевой функции, решатель пытается понизить их значения, то есть пытается установить их все равными 0. Но для тех интервалов, когда включается генератор, линейное неравенство заставляет z(i,j) равным 1.

Создание вспомогательной переменной w который представляет собой y(i+1,j,k) - y(i,j,k). Представление неравенства при запуске генератора с точки зрения w.

w = optimexpr(nPeriods,nGens); % Allocate w
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); % Fill in the levels
generatorlevel(:,1,2) = gen(1,2);
generatorlevel(:,2,1) = gen(2,1);
generatorlevel(:,2,2) = gen(2,2); 

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

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

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

fuelCost = fuelUsed*fuelPrice;

Стоимость запуска генератора = z*startCost.

startingCost = z*startCost;

Прибыль = входящий доход - общая стоимость топлива - стоимость запуска.

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

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

Создайте задачу оптимизации и включите цель и ограничения.

dispatch = optimproblem('ObjectiveSense','maximize');
dispatch.Objective = profit;
dispatch.Constraints.switchcons = switchcons;
dispatch.Constraints.fuelcons = fuelcons;
dispatch.Constraints.powercons = powercons;

Для экономии места подавьте итеративное отображение.

options = optimoptions('intlinprog','Display','final');

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

[dispatchsol,fval,exitflag,output] = solve(dispatch,'options',options);

Solving problem using intlinprog.

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

Анализ решения

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

subplot(3,1,1)
bar(dispatchsol.y(:,1,1)*gen(1,1)+dispatchsol.y(:,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(dispatchsol.y(:,2,1)*gen(1,1)+dispatchsol.y(:,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(dispatchsol.z) == 1); % Use round for noninteger results
[theperiod,thegenerator] = ind2sub(size(dispatchsol.z),starttimes)

theperiod = 2×1

    23
    16

thegenerator = 2×1

     1
     2

Периоды начала работы генераторов соответствуют графикам.

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

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

startCost = 500; % Choose a lower penalty for starting the generators
startingCost = z*startCost;
profit = sum(sum(sum(revenue))) - fuelCost - sum(sum(startingCost));
dispatch.Objective = profit;
[dispatchsolnew,fvalnew,exitflagnew,outputnew] = solve(dispatch,'options',options);

Solving problem using intlinprog.

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

subplot(3,1,1)
bar(dispatchsolnew.y(:,1,1)*gen(1,1)+dispatchsolnew.y(:,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(dispatchsolnew.y(:,2,1)*gen(1,1)+dispatchsolnew.y(:,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(dispatchsolnew.z) == 1); % Use round for noninteger results
[theperiod,thegenerator] = ind2sub(size(dispatchsolnew.z),starttimes)

theperiod = 3×1

    22
    16
    45

thegenerator = 3×1

     1
     2
     2