Main Content

发电机的最佳调度：基于求解器

打开实时脚本

此示例说明如何以最佳方式调度两台燃气发电机，即收益减去成本的值最大。虽然此示例并不完全切合实际，但它能够说明如何考虑取决于决策时机的成本。

要了解如何通过基于问题的方法处理此问题，请参阅发电机的最佳调度：基于问题。

问题定义

电力市场在一天的不同时段有不同用电价格。如果您有发电机，您可以利用这种可变定价，通过安排发电机在价格高的时候运行。假设您控制两台发电机。每台发电机有三个功率水平（关、低和高）。每台发电机在每个功率水平下都有指定的燃料消耗率和发电量。当然，当发电机关闭时，燃料消耗为 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 定义为：

  • 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 时，发电机关闭。

您需要确定发电机关闭后何时启动。假设

  • 发电机 j 在时间段 i 关闭，但在时间段 i + 1 打开，表示为 z(i,j) = 1。在其他时间段 z(i,j) = 0。换句话说，当 sum_k y(i,j,k) = 0 时 z(i,j) = 1，且 sum_k y(i+1,j,k) = 1。

显然，您需要一种方法来根据 y 的设置自动设置 z。以下线性约束可处理此设置。

您还需要成本问题的参数、每台发电机的发电水平、发电机的燃料消耗水平和可用燃料。

  • poolPrice(i) - i 间隔内每兆瓦时收入（美元）。

  • gen(j,k) - 由发电机 j 在功率等级 k 下产生的兆瓦。

  • fuel(j,k) - 发电机 j 在功率等级 k 下使用的燃料。

  • totalfuel - 一天内可用的燃料。

  • startCost - 发电机关闭后再启动的成本（美元）。

  • fuelPrice - 每单位燃料的成本。

当您执行 load dispatchPrice; 时，您得到 poolPrice。按如下所示设置其他参数。

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.ConstraintTolerance = 1e-06.

检查解

检查解的最简单方法是将解向量 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 只要处于运行状态就以其高功率水平运行。发电机 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.ConstraintTolerance = 1e-06.

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

相关主题