MindOpt——优化虚拟电厂智能调度问题(二)

上一案例光储荷经济性调度只考虑到经济性调度的单一目标，这里我们考虑综合多个评测指标，同时以智慧楼宇全天各时段调度为目标，构建一个完整的日调度场景。

本案例来源：https://opt.aliyun.com/example/aqIv2uBnQZfe/src/02_%E6%99%BA%E6%85%A7%E6%A5%BC%E5%AE%87%E5%A4%9A%E7%9B%AE%E6%A0%87%E8%B0%83%E5%BA%A6.ipynb
作者：蒋蔚(蔚绛)

智慧楼宇多目标调度

1. 问题描述

智慧楼宇调度，是在保证社区负荷需求的情况下，通过储能设备的指令控制，以用电经济性、环保性和对电网稳定性为综合目标的一种调度场景。

社区用电来源可以包括：（1）光伏设备发电；（2）储能设备放电；（3）向公网购买。我们考虑多时段多目标优化问题，通过资源协调利用，在满足设备功率约束和能量状态约束的条件下，提供调度决策指令参考值。

本案例部分内容参考自竞赛 NeurIPS 2022 CityLearn Challenge。在这个比赛中，达摩院决策智能实验室团队提出的电力调度AI解决方案，创新性地融合预测技术、优化技术和强化学习，使得实验项目在用电量不变的前提下，碳排放量减少13.6%，电费减低28.2%。，可参阅文章。

2. 数学规划模型

2.1 数学建模

考虑$I$个楼宇在$T$时间段内的调度问题，输入以下参数：

• 楼宇的净负荷需求（负荷和光伏功率的差值）$L=(l_{it})\in {\mathbb{R}}^{I\times T}$ （单位：kWh）
• 电价 $P=(p_t)\in {[0,\infty )}^T$ （单位: $/kWh）
• 碳排价格 $E=(e_t)\in {[0,\infty )}^T$ (单位: kg CO2/kWh)
• 储能充电效率 $c_1\in [0,1]$
• 储能放电效率 $c_2\in [0,1]$
• 储能单次充放电限制 $\alpha \in [0,1]$

为了实现电费、碳排和电网稳定性的多目标，引入决策变量归一化的储能动作：

• 充电变量 $X=(x_{it})\in [0,\alpha {]}^{I\times T}$
• 放电变量 $Y=(y_{it})\in [-\alpha ,0{]}^{I\times T}$
• 储能设备能量状态 $S=(s_{it})\in [0,1{]}^{I\times T}$

最小化以下目标：

$\mathrm{minimize}\text{Direct Cost + Emission Charge + Ramping Charge}$

$\iff$

${\mathrm{minimize}}_{X,Y,S}{\sum }_{t=1}^T\max \left({\sum }_{i=1}^I\left(l_{it}+x_{it}+y_{it}\right),0\right)\cdot p_t+{\sum }_{t=1}^T\left({\sum }_{i=1}^I\max \left(l_{it}+x_{it}+y_{it},0\right)\right)\cdot e_t+{\sum }_{t=2}^T\left|\max \left({\sum }_{i=1}^I(l_{it}+x_{it}),0\right)-\max \left({\sum }_{i=1}^I(l_{i,t-1}+x_{i,t-1}),0\right)\right|$

满足以下约束条件：

• 储能设备功率约束：$0\le x_{it}\le \alpha ,-\alpha \le y_{it}\le 0,i=1,\dots ,I,t=1,\dots ,T$
• 储能能量状态更新：$s_{i0}=0,s_{it}=s_{i,t-1}+c_1\cdot x_{it}+\frac{1}{c_2}\cdot y_{it},i=1,\dots ,I,t=1,\dots ,T$
• 储能能量状态约束：$0\le s_{it}\le 1,i=1,\dots ,I,t=1,\dots ,T$

2.2 线性化模型

在本案例中，我们引入了变量替换来实现问题线性化。

• 楼宇净消费电量 $\mathrm{N}=({\nu }_{it})\in [0,\infty {)}^{I\times T}$
• 区域净消费电量 $\mathrm{M}=({\mu }_t)\in [0,\infty {)}^T$
• 区域爬坡电量 $\Omega =({\omega }_t)\in [0,\infty {)}^{T-1}$

则最小化目标等价为：

${\mathrm{minimize}}_{X,Y,S,\text{N},\text{M},\Omega }{\sum }_{t=1}^T{\mu }_t\cdot p_t+{\sum }_{t=1}^T{\sum }_{i=1}^I{\nu }_{it}\cdot e_t+{\sum }_{t=1}^{T-1}{\omega }_t$

满足以下约束条件：

• 储能设备功率约束：$0\le x_{it}\le \alpha ,-\alpha \le y_{it}\le 0,i=1,\dots ,I,t=1,\dots ,T$
• 储能能量状态更新：$s_{i0}=0,s_{it}=s_{i,t-1}+c_1\cdot x_{it}+\frac{1}{c_2}\cdot y_{it},i=1,\dots ,I,t=1,\dots ,T$
• 储能能量状态约束：$0\le s_{it}\le 1,i=1,\dots ,I,t=1,\dots ,T$
• 区域净功率约束：${\mu }_t\ge 0,{\mu }_t\ge {\sum }_i^I(l_{it}+x_{it}+y_{it}),t=1,\dots ,T$
• 楼宇净功率约束：${\nu }_{it}\ge 0,{\nu }_{it}\ge l_{it}+x_{it}+y_{it},i=1,\dots ,I,t=1,\dots ,T$
• 区域爬坡状态：$-{\omega }_t\le {\mu }_{t+1}-{\mu }_t\le {\omega }_t,t=1,\dots ,T$

下面使用MAPL建模语言求解该问题。

3. 使用MAPL建模和求解

根据MindOpt APL的建模语法，建立此经济调度的完整模型，可如下两种方式运行。

3.1. 方法1：cell中直接输入代码

在下面的cell中运行如下的代码来完建模和求解:

clear model;   #清除环境内存，多次run的时候使用
option modelname ../model/smart_building; #定义存储文件名

set I := 1..2;
set T := 1..24;
set T0 := 1..23; # T时段内部共T-1个状态点
set T1 := 1..25; # T时段前后共T+1个状态点


#### 负荷矩阵作为参数输入
param load_matrix[I*T] := 
    | 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24|
|1| 0.10844531,  0.10140885,  0.10111719,  0.10541146,  0.29224219,
         0.31547396,  0.12083073,  0.03450521, -0.10114063, -0.2266276 ,
        -0.31739323, -0.38168229, -0.40365365, -0.39797917, -0.30016406,
        -0.25292187, -0.13258333, -0.00121875,  0.5913776 ,  0.53110677,
         0.27235938,  0.21871354,  0.19486198,  0.55161979 |
|2| 0.06417708,  0.14725521,  0.05614323,  0.05621615,  0.05623177,
         0.05604688,  0.06448958,  0.06498698,  0.136375  , -0.0548151 ,
        -0.15175781, -0.22945833, -0.28111719, -0.27846615, -0.05311458,
        -0.06278125,  0.17883594,  0.4074974 ,  0.18320052,  0.22567187,
         0.27401302,  0.28129688,  0.172125  ,  0.46530469|;

#### 单位电价和碳排价格
param price[T] := <17> 0.54, <18> 0.54, <19> 0.54, <20> 0.54, <21> 0.54  default 0.22; 
param emission[T] := <1> 00.10920568, <2> 0.10281591, <3> 0.10226245, <4> 0.09849792, <5> 0.09972523,
       <6> 0.10239134, <7> 0.10813721, <8> 0.11276986, <9> 0.12070297, <10> 0.13228371,
       <11> 0.13734146, <12> 0.13988744, <13> 0.13140874, <14> 0.12927128, <15> 0.13609823,
       <16> 0.13802977, <17> 0.14356089, <18> 0.14756169, <19> 0.15092846, <20> 0.16013497,
       <21> 0.15563251, <22> 0.14269217, <23> 0.13361955, <24> 0.12968476;

#### 充放电损耗系数
param charge_efficiency := 0.92;
param discharge_efficiency := 0.95;

#### 储能单次充放电限制
param nominal_power := 0.8;

#### 在不使用储能情况下的各指标，作为评测的对比基准值
param direct_cost_base := sum{t in T} max(sum{i in I} load_matrix[i,t], 0) * price[t];
param emission_charge_base := sum{t in T, i in I}  max(load_matrix[i,t], 0) * emission[t];
param ramping_charge_base := sum{t in T0} abs(max(sum{i in I} load_matrix[i,t],0) - max(sum{i in I} load_matrix[i,t+1],0));

#### 定义优化问题的变量
var x[I*T] >= 0 <= nominal_power; # 充电变量
var y[I*T] <= 0 >= - nominal_power; # 放电变量
var soc[I*T1] >= 0 <= 1; # State-of-Charge 储能设备状态
var vu[I*T] >= 0; # 楼宇净消费电量
var mu[T] >= 0; # 区域净消费电量
var w[T0] >= 0; # 区域爬坡电量

#### 定义优化目标函数
minimize Average_Cost:
    1/3 * sum{t in T} mu[t] * price[t] / direct_cost_base +
    1/3 * sum{t in T, i in I} vu[i,t] * emission[t] / emission_charge_base +  
    1/3 * sum{t in T0} w[t] / ramping_charge_base;

#### 定义约束条件
subto District_Net_Consumption: # 区域净功率约束
    forall {t in T} do
        mu[t] - sum {i in I} x[i,t] - sum{i in I} y[i,t] >= sum{i in I} load_matrix[i, t];

subto Building_Net_Consumption: # 楼宇净功率约束
    forall {t in T ,i in I} do
        vu[i,t] - (x[i,t] + y[i,t]) >= load_matrix[i, t];

subto Rampling: # 区域爬坡状态
    forall {t in T0} do
        w[t] - (mu[t+1] - mu[t]) >= 0 and 
        w[t] - (mu[t] - mu[t+1]) >= 0;

subto State_of_Charge_Init: # 储能设备初始状态
    forall {i in I} do soc[i,1] == 0;

subto State_of_Charge_Update: # 储能设备状态更新
    forall {i in I, t in T} do
        soc[i,t+1] == soc[i,t] + x[i,t] * charge_efficiency + y[i,t] / discharge_efficiency;    
      
option solver mindopt;     # 指定求解用的求解器
solve;         # 求解

#### 展示结果
print "-----------------Display---------------"; 
print "Charging(+) / Discharging (-) Solution (only print non-zero results):";
forall {i in I, t in T with x[i,t] > 0 or y[i,t] < 0} print 'Time {}, Building {}, Solution: {}'  % t, i, x[i,t] + y[i,t];

#### 评测指标
print "-----------------Evaluation Metrics---------------";
print "Average Score = "; 
print 1/3 * sum{t in T, i in I} vu[i,t] * emission[t] / emission_charge_base +
        1/3 * sum{t in T} mu[t] * price[t] / direct_cost_base +
        1/3 * sum{t in T0} w[t] / ramping_charge_base;
print "Direct Cost = ";
print sum{t in T} mu[t] * price[t] / direct_cost_base;
print "Emission Charge = ";
print sum{t in T, i in I} vu[i,t] * emission[t] / emission_charge_base;
print "Ramping Charge =";
print sum{t in T0} w[t] / ramping_charge_base;

运行代码得到结果如下：

Running mindoptampl
wantsol=1
mip_integer_tolerance=1e-9
MindOpt Version 0.23.0 (Build date: 20221129)
Copyright (c) 2020-2022 Alibaba Cloud.

Start license validation (current time : 02-MAR-2023 13:49:38).
License validation terminated. Time : 0.003s

Model summary.
 - Num. variables     : 241
 - Num. constraints   : 166
 - Num. nonzeros      : 594
 - Bound range        : [1.2e-03,1.0e+20]
 - Objective range    : [3.4e-02,1.3e-01]
 - Matrix range       : [9.2e-01,1.1e+00]

Presolver started.
Presolver terminated. Time : 0.001s

Simplex method started.
Model fingerprint: ==gZ3dmbidHYkV2dhFmZ

    Iteration       Objective       Dual Inf.     Primal Inf.     Time
            0     0.00000e+00      0.0000e+00      1.0646e+01     0.00s  
          226     5.18050e-01      0.0000e+00      0.0000e+00     0.00s  
Postsolver started.
Simplex method terminated. Time : 0.003s


OPTIMAL; objective 0.52
226 simplex iterations

Completed.
-----------------Display---------------
Charging(+) / Discharging (-) Solution (only print non-zero results):
Time 1, Building 1, Solution: 0.07608708880640472
Time 2, Building 1, Solution: 4.54188064047123e-05
Time 3, Building 1, Solution: 0.09144905880640467
Time 4, Building 1, Solution: 0.022955003817846265
Time 5, Building 1, Solution: -0.0997644811935953
Time 6, Building 1, Solution: -0.0667644811935953
Time 13, Building 1, Solution: 0.40365365
Time 14, Building 1, Solution: 0.13021694173913048
Time 15, Building 1, Solution: 0.30016406
Time 16, Building 1, Solution: 0.25292187
Time 19, Building 1, Solution: -0.4430963499999996
Time 20, Building 1, Solution: -0.3828255199999995
Time 21, Building 1, Solution: -0.12407812999999951
Time 4, Building 2, Solution: 0.06412686498855844
Time 6, Building 2, Solution: -0.05604688
Time 10, Building 2, Solution: 0.0548151
Time 11, Building 2, Solution: 0.12720392173913048
Time 12, Building 2, Solution: 0.22945833
Time 13, Building 2, Solution: 0.28111719
Time 14, Building 2, Solution: 0.27846615
Time 15, Building 2, Solution: 0.05311458
Time 16, Building 2, Solution: 0.06278125
Time 18, Building 2, Solution: -0.25799739999999954
Time 19, Building 2, Solution: -0.18320052
Time 20, Building 2, Solution: -0.22567187
Time 21, Building 2, Solution: -0.27401302
Time 22, Building 2, Solution: 0.14247976107714755
Time 23, Building 2, Solution: 0.2755032010771475
Time 24, Building 2, Solution: -0.37443429892285235
-----------------Evaluation Metrics---------------
Average Score = 
0.5180502138863339
Direct Cost = 
0.5626249732899371
Emission Charge = 
0.6740691491806657
Ramping Charge =
0.3174565191883989

---

3.2. 方法2：导入.mapl文件来运行

可以把中间的建模代码存在一个 smart_building.mapl的文件里，然后cell中用 model ../model/smart_building.mapl替代建模环节。如下cell中运行：

clear model;  #清除model，多次run的时候使用

#--------------------------
model ../model/smart_building.mapl;
#--------------------------

option solver mindopt;     # （可选）指定求解用的求解器，默认是MindOpt
solve;         # 求解

运行结果如同方法一

4 求解结果分析

在执行以上命令后，MAPL会打印如下的结果。这里我们选择只打印出非零的结果：

---

-----------------Display---------------
Charging(+) / Discharging (-) Solution (only print non-zero results):
Time 1, Building 1, Solution: 0.01196022381784631
Time 2, Building 1, Solution: 4.541880640474005e-05
Time 3, Building 1, Solution: 0.09144905880640472
Time 4, Building 1, Solution: 0.08708186880640473
Time 5, Building 1, Solution: -0.09976448119359539
Time 6, Building 1, Solution: -0.06676448119359527
Time 11, Building 1, Solution: 0.31739323
Time 12, Building 1, Solution: 0.38168229
Time 13, Building 1, Solution: 0.38788100173913037
Time 19, Building 1, Solution: -0.4430963499999996
Time 20, Building 1, Solution: -0.3828255199999995
Time 21, Building 1, Solution: -0.12407812999999951
Time 1, Building 2, Solution: 0.06412686498855844
Time 6, Building 2, Solution: -0.05604688
Time 10, Building 2, Solution: 0.0548151
Time 11, Building 2, Solution: 0.12720392173913048
Time 12, Building 2, Solution: 0.22945833
Time 13, Building 2, Solution: 0.28111719
Time 14, Building 2, Solution: 0.27846615
Time 15, Building 2, Solution: 0.05311458
Time 16, Building 2, Solution: 0.06278125
Time 18, Building 2, Solution: -0.25799739999999954
Time 19, Building 2, Solution: -0.18320052
Time 20, Building 2, Solution: -0.22567187
Time 21, Building 2, Solution: -0.27401302
Time 22, Building 2, Solution: 0.14247976107714755
Time 23, Building 2, Solution: 0.2755032010771475
Time 24, Building 2, Solution: -0.37443429892285235
-----------------Evaluation Metrics---------------
Average Score = 
0.5180502138863339
Direct Cost = 
0.5626249732899371
Emission Charge = 
0.6740691491806657
Ramping Charge =
0.3174565191883989

---

此结果展示了各时刻各楼宇的最优调度策略，正向代表充电，负向代表放电。同时，我们展示了多维度的评测指标，我们观察到在使用储能之后，电费降低44%，碳排可以优化33%，电网稳定性提升了68%。

结语

该智慧楼宇多目标调度案例考虑的是给定负荷和功率真实值情况下的调度。在实际中，我们会遇到更多的挑战，例如，在日前时间点，我们无法提前得知实际的负荷和功率，因此在调度决策时需要使用预测值。而预测值和实际值之间的偏差向问题中引入了不确定性。面对这些挑战，我们可以结合预测、优化、强化学习等技术来应对不确定性。