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