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

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

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

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

Можно сформулировать проблему планирования как бинарную проблему целочисленного программирования можно следующим образом. Задайте индексирует iJ, и 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