来自于本人自学yalmip优化的一些理解,主要参考:
1.yalmipcommand手册Commands (reference manual) - YALMIP
2.Yalmip工具箱使用教程_配电网和matlab的博客-CSDN博客
引入官网正文(主要参考deepl机翻)。
1、机组组合模型的基本参数
机组组合问题是调度和整数编程中的一个经典问题。在这个问题中,我们的任务是开启和关闭发电厂,以满足预测的未来电力需求,同时最大限度地降低成本。我们有几个不同的发电厂,它们具有不同的特性和运行成本,在如何使用这些发电厂方面也有各种限制。我们将从一个非常简单的模型开始,然后用更高级的功能来扩展这个模型。
首先,我们假设有 N = 3 个不同类型的发电厂。每个发电厂都有一个最大发电量和一个最小发电量。我们的调度问题在 48 个时间单位(小时、天......)内求解
电力需求预测由一个周期函数给出。
发电厂在一个时间单位内的运行成本由二次函数 给出,其中
为
号发电厂输出功率。
Nunits=3;
Horizon=48;
Pmax=[100;50;25];
Pmin=[20;40;1];
%% The cost of running plant i for one time unit is given by a quadratic function
% QiiPi^2+CiPi
% where Pi is the delivered powerfrom plant Pi
Q=diag([.04 .01 .02]);
C=[10 20 20];
Pforecast=100+50*sin((1:Horizon)*2*pi/24);
%% the forecasted power demand is given by a periodic function.
2、线性开关模型
模型由两个变量组成。第一个变量是二进制(逻辑)变量,用于控制 i 号机组在 k 时刻是否开启;第二个变量是连续变量,代表 i 号机组在 k 时刻的输出功率。
onoff = binvar(Nunits,Horizon,'full');
P = sdpvar(Nunits,Horizon,'full');
目前,在对这些问题进行建模时,一个典型的错误是开始将开断变量与功率变量 P 相乘,其逻辑是发电厂的功率由onoff与 P 的乘积给出。如果 onoff 为 0,则 P 应为零;如果 onoff 为 1,则应允许 P 介于下限和上限之间。
(本文以48个时间单元组成一个运行周期,使用for循环以K=1开始逐个对每个时间单元进行操作,内嵌套for unit =1:Nunit在每一个时间单元下对不同机组逐个操作。)
for k = 1:Horizon
Constraints = [Constraints, onoff(:,k).*Pmin <= P(:,k) <= onoff(:,k).*Pmax];
end
%%onoff(:,k) 取onoff矩阵的第k列
%%onoff是逻辑变量矩阵,与P点乘可得P的原值或0,以此模拟控制电厂开关。
%%加入for 1:horizon模拟一个周期。
%%P每一列代表三个电厂在一时间单元下的发电量。Pmin<=P<=Pmax。
每次交付的输出功率必须满足Pforecast函数。
for k = 1:Horizon
Constraints = [Constraints, sum(P(:,k)) >= Pforecast(k)];
end
参考上述提到的运行成本二次函数,可计算出一个时间周期下的运行成本。
for k = 1:Horizon
Objective = Objective + P(:,k)'*Q*P(:,k) + C*P(:,k);
end
%%计算一个周期下的运行成本
%% 线性代数知识P(:,k)'*Q*P(:,k)=QiiPi^2,P(:,k)'为P的转置。
至此,基本模型已经建立完毕,我们可以解题并显示结果了(这需要安装一个混合整数二次规划(QP)求解器,笔者安装的是Gurobi)。
ops = sdpsettings('verbose',1,'debug',1);
optimize(Constraints,Objective,ops)
stairs(value(P)');%%绘制阶梯图
legend('Unit 1','Unit 2','Unit 3');%%标注图例
可生成如下图像。
3、考虑最小运行/停止时间
上述模型与实际情况相去甚远,其中一个主要缺陷是它假定我们可以任意开启和关闭工厂。在大多数情况下,如果开启设备,它必须运行若干个个时间单位,而当关闭时,它又会消失若干个时间单位(重新启动很复杂等)。我们定义了两个新变量来描述这些设备特性。
%%设置最小运行/停止时间
minup = [6;30;1];
mindown = [3;6;3];%%unit1、2、3的最小启动/停止时间
如code所示,unit1的最小运行时间为6个时间单位,最小停止时间为3个时间单位。
现在,让我们通过涉及二进制变量的线性约束来描述这些特性。我们从最短开机时间约束开始。
如果设备在i时刻处于开启状态,由于受到最小运行时间的限制,它必须在 minup(i) 个长度的时间单位保持运行状态。这可以表述为:如果设备在 k-1 时关闭,在 k 时开启,那么它必须在k+1、k+2...、k+minup(i)-1,保持运行状态。
现在的诀窍是得到一个开启的指标(indicator)。在机组在K时刻刚启动情况下,onoff(k)-onoff(k-1) 值为 1,否则就是 0 (持续运行或停止)或-1(刚停止)。因此下面的模型就可以用了。
%%设置最小启动/停止时间
minup = [6;30;1];
mindown = [3;6;3];%%unit1、2、3的最小启动/停止时间
for k = 2:Horizon
for unit = 1:Nunits
% indicator will be 1 only when switched on
indicator = onoff(unit,k)-onoff(unit,k-1);
%%matrix(row,col)访问row行,col列的元素,修改可直接加上赋值号。
%%根据不同情况indicator可能取值为-1 0 1
range = k:min(Horizon,k+minup(unit)-1);
%%每次循环至新的时间单元重新赋值range的起始范围,使其范围不变
% Constraints will be redundant unless indicator = 1
Constraints = [Constraints, onoff(unit,range) >= indicator];
%%此约束当且仅当indicator=1时生效,使当前机组range范围内的onoff置1,处于运行状态。
%%当indicator=0或-1时,constraints无作用。
end
end
%%最小停止时间写法思路同上
for k = 2:Horizon
for unit = 1:Nunits
% indicator will be 1 only when switched off
indicator = onoff(unit,k-1)-onoff(unit,k);
range = k:min(Horizon,k+mindown(unit)-1);
% Constraints will be redundant unless indicator = 1
Constraints = [Constraints, onoff(unit,range) <= 1-indicator];
end
end
ops = sdpsettings;
optimize(Constraints,Objective,ops);
stairs(value(P)');
legend('Unit 1','Unit 2','Unit 3');
生成图像如下:
4、量化功率水平
有些发电厂只能在有限的几种配置下运行,从而使输送的电能功率是一个量化变量。在此,我们假定 3 号电站受到这种限制,并使用 ismember 指令对量化进行建模。
在此,我们假设 3 号电站受到这种限制,并使用等式算子对其进行量化。
Unit3Levels = [0 1 6 10 12 20];
for k = 1:Horizon
Constraints = [Constraints, ismember(P(3,k),Unit3Levels)];
%%ismember函数用于约束变量的取值位于某个离散集中,P(3,k)只取unit3level中的值
end
optimize(Constraints,Objective);
stairs(value(P)');
legend('Unit 1','Unit 2','Unit 3');
生成图片如下:
笔者从自己的理解出发对代码加以注解,如有错误,恳请指正。