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

Рынок электроэнергии имеет разные цены в разное время суток. Если у вас есть генераторы, вы можете воспользоваться этой переменной ценой, планируя ваши генераторы работать, когда цены высоки. Предположим, что есть два генератора, которыми вы управляете. Каждый генератор имеет три уровня мощности (выключенный, низкий и высокий). Каждый генератор имеет заданную скорость потребления топлива и выработки энергии на каждом уровне мощности. Конечно, расход топлива равен 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 автоматически на основе настроек 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 соответствующего размера, создайте нулевые матрицы размеров y и z матрицы.

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