yalmip example-unit commitment代码注释（1）

来自于本人自学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');

  生成图片如下：

  笔者从自己的理解出发对代码加以注解，如有错误，恳请指正。