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

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

Вам нужно определить, когда генератор запускается после того, как быть off. Let

Очевидно, тебе нужен способ задать 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 переменные, автоматически совпадающие с периодами активный/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 переменные в целевой функции, решатель пытается снизить значения 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