이 번역 페이지는 최신 내용을 담고 있지 않습니다. 최신 내용을 영문으로 보려면 여기를 클릭하십시오.

발전기 최적 운용법: 솔버 기반

이 예제에서는 수입에서 비용을 뺀 값을 최대화할 수 있도록 가스 연소식 발전기 2대의 가동 일정을 최적화하는 방법을 다룹니다. 이 예제가 완전히 현실적인 상황은 아니지만, 결정 시점에 따라 변하는 비용을 염두에 두는 점은 분명히 다루고 있습니다.

이 문제에 대한 문제 기반 접근법은 발전기 최적 운용법: 문제 기반 항목을 참조하십시오.

문제 정의

전력 시장에서는 하루에도 시간대에 따라 가격이 계속 달라집니다. 발전기를 보유하고 있을 경우 가격이 높게 형성되는 시점에 발전기를 가동하도록 일정을 설정함으로써 이러한 가변적 가격 책정을 유리하게 이용할 수 있습니다. 두 대의 발전기를 제어한다고 가정하겠습니다. 각 발전기에는 끄기, 낮음, 높음의 세 가지 발전 수준이 있습니다. 각 발전기는 각 발전 수준마다 연료 소비율과 전력 생산율이 정해져 있습니다. 물론, 발전기가 꺼져 있을 때는 연료 소비율이 0입니다.

하루 중 30분의 시간 간격마다(하루는 24시간이므로 총 48개의 시간 간격이 있음) 각 발전기의 발전 수준을 지정할 수 있습니다. 과거 기록을 기준으로 각 시간 간격에서 얻을 수 있는 시간당 발전 메가와트(MWh)당 수입을 알고 있다고 가정할 수 있습니다. 이 예제에 사용할 데이터는 2013년도 중반 AEMO(Australian Energy Market Operator, https://www.nemweb.com.au/REPORTS/CURRENT/)에서 얻은 것으로, 사용약관(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입니다.

• i번째 주기에 발전기 j가 발전 수준 k로 가동할 때 y(i,j,k) = 1입니다. 낮은 전력이 k = 1이 되도록 하고 높은 전력은 k = 2가 되도록 합니다. sum_k y(i,j,k) = 0일 때는 발전기가 꺼진 것입니다.

발전기가 꺼진 후 시동되는 시점을 결정해야 합니다. 다음과 같이 가정하겠습니다.

• i번째 주기에 발전기 j가 꺼지고 i + 1번째 주기에 켜지면 z(i,j) = 1입니다. 그렇지 않으면 z(i,j) = 0입니다. 다시 말해, sum_k y(i,j,k) = 0이고 sum_k y(i+1,j,k) = 1일 때 z(i,j) = 1입니다.

분명히, y의 설정을 기준으로 z를 자동으로 설정할 방법이 필요합니다. 아래의 선형 제약 조건에서 이 설정을 다룹니다.

비용, 각 발전기의 발전 수준, 발전기의 소비 수준, 사용 가능한 연료에 대한 문제의 파라미터도 필요합니다.

• poolPrice(i) -- 간격 i에서 MWh당 수입(단위: 달러)

• gen(j,k) -- 발전기 j가 발전 수준 k에서 발전하는 전력(단위: MW)

• 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

주기당 가동 비용은 해당 기간 동안의 연료비입니다. k 수준에서 가동하는 발전기 j의 경우, 소요되는 비용은 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일 때.

다음 사실을 알 수 있습니다.

z(i,j) = 1을 원하는 경우 sum_k ( - y(i,j,k) + y(i+1,j,k) ) > 0입니다.

따라서 문제 정식화에 다음 선형 부등식 제약 조건

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.RelativeGapTolerance = 0.0001. The intcon variables are integer within tolerance, options.IntegerTolerance = 1e-05.

해 검토하기

가장 쉽게 해를 검토하는 방법은 해 벡터 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(2,1)+ysolution(:,2,2)*gen(2,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가 더 효율적이므로 예상할 수 있는 바입니다. 발전기 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.RelativeGapTolerance = 0.0001. The intcon variables are integer within tolerance, options.IntegerTolerance = 1e-05.

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(2,1)+ysolutionnew(:,2,2)*gen(2,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

관련 항목

• 발전기 최적 운용법: 문제 기반
• Integer and Logical Modeling