Оптимальная отправка производителей электроэнергии: основанный на решателе

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

Для подхода, основанного на проблеме к этой проблеме смотрите Оптимальную Отправку Производителей электроэнергии: основанный на проблеме.

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

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

Figure contains an axes object. The axes object contains an object of type bar.

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

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

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

  • nPeriods = количество периодов времени, 48 в этом случае.

  • i = период времени, 1 <= i <= 48.

  • j = индекс генератора, 1 <= j <= 2 для этого примера.

  • y(i,j,k) = 1 когда период i, генератор j работает с уровнем мощности k. Позвольте малой мощности быть k = 1, и большая мощность быть k = 2. Генератор выключен когда sum_k y(i,j,k) = 0.

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

  • z(i,j) = 1 когда генератор j выключено в период i, но включен в период i + 1. z(i,j) = 0 в противном случае. Другими словами, z(i,j) = 1 когда sum_k y(i,j,k) = 0 и sum_k y(i+1,j,k) = 1.

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

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

  • poolPrice(i) - Доход в долларах на МВт·ч в интервале i.

  • gen(j,k) - MW, сгенерированный генератором j на уровне мощности k.

  • fuel(j,k) - Топливо, используемое генератором j на уровне мощности k.

  • totalfuel - Топливо, доступное за один день.

  • startCost - Стойте в долларах, чтобы запустить генератор после того, как это будет выключено.

  • fuelPrice - Стойте за модуль топлива.

Вы получили 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')

Figure contains an axes object. The axes object contains 2 objects of type bar. These objects represent Generator 1, Generator 2.

Заметьте, что генератор 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')

Figure contains 3 axes objects. Axes object 1 with title Generator 1 optimal schedule contains an object of type bar. Axes object 2 with title Generator 2 optimal schedule contains an object of type bar. Axes object 3 with title Energy price contains an object of type bar.

Генератор 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')

Figure contains 3 axes objects. Axes object 1 with title Generator 1 optimal schedule contains an object of type bar. Axes object 2 with title Generator 2 optimal schedule contains an object of type bar. Axes object 3 with title Energy price contains an object of type bar.

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

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте