Описание задачи

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

Каждому генератору можно присвоить уровень степени для каждого получасового временного интервала в течение дня (24 часа, так 48 интервалы). Основываясь на исторических записях, предположим, что вы знаете доход за мегаватт-час (MWh), который вы получаете в каждом временном интервале. Данные для этого примера получены от Австралийского оператора энергетического рынка 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 переменные, автоматически совпадающие с периодами активный/off 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