风电随机性动态经济调度模型（Matlab代码实现）

💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
📋📋📋本文内容如下：🎁🎁🎁

⛳️赠与读者

👨‍💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能解答你胸中升起的一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。
或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎

💥内容介绍

0 写在前面

大规模风电并网给经济调度模型带来的主要挑战是如何处理风电的随机性，如何降低其对电力系统优化调度的影响。建模是大规模风电并网电力系统优化调度的核心步骤，因此，我们回顾了风电并网经济调度的几种典型建模方法。此外，还用了随机经济调度程序。

1 引言

风电并网对经济调度的影响主要在于如何处理随机风电出力，以减少其不确定性对调度结果的影响。因此，在相关模型中，现有研究的研究重点主要集中在处理风电随机性方法的差异上。根据风电随机性处理方式的不同，目前涉及风电经济调度问题的常用建模方法主要分为机会约束规划方法、基于场景的方法和模糊模式.

1.1 机会约束规划

偶然约束编程在经济调度问题中，由于风力发电的不确定性，对约束条件进行确定性处理后，计算结果肯定是保守的。随机机会的约束编程是随机编程的一个重要分支，它主要是针对含有随机变量的约束条件的优化问题，其决策必须在观察到随机变量的实现之前做出。考虑到在不利条件下做出的决策可能不满足约束条件，允许决策在一定程度上不满足约束条件。但是，决策应保证约束条件出现的概率不低于一定的置信度。机会约束编程的一般形式是：

在公式(1)中，f(x)是目标函数；X是决策向量；ξ是参数向量；g(x,ξ) ≤0是约束条件；α是置信度；Pr{~}是概率。

1.2 基于场景的方法

基于场景的方法是解决随机过程问题的一个重要工具。德国柏林洪堡大学的Heitsch H、Romisch W等人在证明情景树及其使用方法方面进行了许多研究，他们讨论了情景树方法在多阶段随机过程中的应用、详细介绍了情景树的生成、缩减和重建方法。

1）目前，传统的思路是在保证一定精度的前提下缩减情景，分为情景生成、情景缩减和情景树重建三个步骤。

(1)相关文章提出的安全约束单元组合方法考虑了风电机组的无功特性，通过引入发电机功率分配系数以及将基于蒙特卡洛仿真的情景法替换为基于极值情景集的情景法，提高了计算时间和计算精度。

(2)基于情景的方法基础上提出了动态减少多切割方法也被提出。首先，风电场一体化的随机机组承诺是分解为主问题和子问题。主问题用于确定发电机组的开/关状态以及与预测的风电情况相对应的输出调度，而子问题则用于确定与采样情况相对应的发电量。然后，主问题和子问题通过最优切割联系起来，交替解决。

2）第二种思路是利用分解加速算法和高性能计算设备直接求解包括大规模风电场景的优化模型。美国阿贡国家实验室基于第二种思路的一些研究成果走在了前列。他们将抽样平均近似法应用于伊利诺伊州电网的随机编程问题。

首先，他们用SAA方法将随机问题转化为确定性的最优问题；然后，用内点法来解决这个确定性的模型。在迭代过程中，根据预测情景和偏离情景的顺序，修正方程的系数矩阵被重新排列为块边对角线形式（BBDF）。采用含有Schur补数的高斯消除法来加速修正方程的求解。特别是在一些研究中，计算涉及19.5亿个变量和16384个场景，算法在超级计算机Titan Cray XK7中编译和运行。含有Schur补数的高斯消除法可以使内点法具有并行计算的效果，在商业软件包OOPS和中的IPOPT都采用了这种方法。

1.3 模糊模型

模糊建模的基本思想是将不确定的风电输出作为一个模糊数，通过成员函数对风电场输出进行模糊处理，建立风电场输出成员函数。该值越大，意味着决策者对风电场输出的满意度越高。同时，建立目标函数和约束条件的成员函数，包括总成本的成员函数，负荷平衡的成员函数，以及旋转储备的成员函数。然后，根据模糊集理论的最大和最小定律，将全局成员定义为整个模型的模糊满足指数：

这样，寻求模型的解就转化为寻求满足约束条件的满足指数ω0的最大值。 Ω0的值反映了调度员对整个优化结果的总体满意度，包括总购电成本、负荷平衡和旋转储备。当Ω0的值接近或等于1时，说明优化后的调度结果满足了调度员的要求。

2 六种处理风电随机性动态经济调度模型描述

2.1 考虑正负旋转储备约束的确定性模型

模型 I 是考虑风电场输出的确定性模型。它将预测的风电场出力作为已知因素，通过预留足够的系统旋转备用容量来应对风电出力不确定性和负荷预测误差的影响。正旋转备用容量用于补偿高估风电出力或低估系统负荷造成的影响，负旋转备用容量用于补偿低估风电出力或高估系统负荷造成的影响。这样，通过优化得到的发电调度方案就有足够的系统旋转备用容量来应对次日风电场实际出力或负荷功率与预测值的偏差。数学模型如下：

其中，公式（1a）是目标函数，定义为购电成本。约束条件（1b）-（1g）与约束条件（4b）-（4g）相同，对于常规动力装置的输出变量和风机输出变量，去掉与方案有关的上标，相当于对单一预测方案的优化。约束条件（1h）是系统的正负旋转储备约束。在目前的优化模型中，利用纺纱储备应对风电并网带来的系统调度的不确定性是一种广泛使用的方法。为了保证几个模型比较的标准，我们只考虑风电预测的误差，不考虑负荷预测的误差。其中，第一和第二式为系统的正、负旋转储备约束，Sut 和,Sdt 为系统在时间间隔t的正、负旋转储备容量；为风电场w的集成功率，Pwt为风电场w的实际调度输出功率；wu%和wd%分别为正、负旋转储备的风电集成需求系数，本文选择30%。每个单元在响应时间内提供的旋转备用约束，和 ,分别是单元 g 在时间间隔 t 的正和负旋转备用容量，T10 是任意两个相邻时间间隔之间的周期，即 10 分钟。

2.2 机会约束的随机模型【机会】

模型二采用机会约束模型来反映风力发电的随机性，认为风力发电量的预测误差被视为服从正态分布的随机变量，采用机会约束规划处理风电出力约束。具体如下：

其中，c为置信度。采用参数估计的方法，将随机产生的S个误差样本代入Matlab的正态分布参数估计函数中，估计出样本的均值，即期望值。取这个平均值作为风电出力的上限，代入模型，于是机会约束模型转化为确定性模型：
其中，样本平均值。模型二不包括旋转储备约束，目标函数和其他约束与模型一相同。

2.3 基于概率最优功率流的随机模型（POPF）【机会】

第三个模型仍然使用机会约束模型。受[16]中概率最优功率流约束的启发，采用考虑风电不确定性的概率最优功率流来建立机会约束的编程模型。基于风力输出预测误差有上下限的假设，机会约束编程模型被转换为确定性模型。具体如下：

其中，Pwt,maxlimit和Pwt,minlimit值等于模型2.5中的两个极限误差情况，用Matlab正态分布函数的反函数来解决上述不等式，并将上述约束条件转化为。
其中，σ为标准差； αc 是该置信水平下对应的变量值。模型2.3不包括旋转储备约束，目标函数和其他约束与模型2.1相同。

2.4 基于场景法的确定性模型【场景】

第四个模型采用情景法，由于情景转移约束的存在，预测的情景输出可以转化为不同的误差情景输出，以及对应于不同风电场输出情景的、满足电力系统有功约束的常规动力装置的最优输出。

其中，ag 是常规单位的上网电价； awind 是风电上网电价。 NG、NL、ND、NW和NPS分别是常规机组、输电线路、负荷、风电场和抽水蓄能电站的集合。 N 是场景集，NS={0,1,2,..., S }，s=0 是预测情景； T为调度间隔集合，T={1, 2,... , 96}； T15为机组正常用电调度的时间间隔，本文选取15min。 ΔT 为场景转移重调度的时间间隔；为时间间隔 t 对应第 s 场景的单位 g 的有功输出，其他与有功功率相关的变量按此定律命名； rd,g 和 ru,g 分别为单位 g 的滑坡和斜坡速率（MW/min）。
式(4b)为功率平衡约束；公式（4c）和（4e）是常规机组和风力发电机的输出下限和下限约束。式(4d)为常规单元滑坡/坡道约束；式(4f)为抽水蓄能电站的上下限约束和抽水平衡约束[15]，ξ为抽水平衡的相关系数，取75%。公式（4g）为线路输电约束，本文采用直流潮流模型[17]对其进行描述。式(4h)为误差情景与预测情景之间的转移约束，表示常规机组应能在规定时间内将预测情景下的产出调整为误差情景下的产出。

2.5 考虑极限情景约束的确定性模型【场景】

模型采用极限误差情景模型,并考虑情景转移约束。该模型与模型2.4相似，不同的是情景集NS' ，NS'={0，1，2} 的情景数，s=0 是预测情景，s=1 是最大正误差情景，s=2 是最大负误差情景。下图是极限误差情景生成图。

最大正/负误差情景

如图所示，两个极限误差情景是S误差样本的上下限边界，模型V的目标函数和约束条件与模型IV相同。

2.6 基于样本平均法的确定性模型

模型六采用六种样本平均法（SAA），即分别求解每种误差情况，最后得到不同情况下购电成本的平均值，即预期购电成本。具体如下：

其中，式(6b)-6g)同式(4b)-(4g)，NS' ={1,2,...,S}；公式（6h是不可预测性抑制约束，其物理意义是假设t=1是已经过去的时间，所以所有场景的输出是相同的。根据伯努利大数定律，将(6h)代入(6a-1)，模型VI的目标函数转化为：

必然规律往往表现为大量的随机事件，而这种现象就是大数定律。通俗地讲，定理就是在不变的条件下进行大量的测试，随机事件的频率接近于它的概率。模型2.2的变量个数与模型2.3相同，为V=T·（NG+NW+2·NPS）。模型2.1的变量个数与模型2.3相同，均为3V；模型IV的变量个数与模型VI相同，为(S+1)·V。计算规模由变量个数和约束个数决定，计算规模会影响计算时间。

3 总结与展望

建模是包含风电并网的电力系统优化调度问题的最重要部分。本文采用机会约束规划法和随机情景法进行建模，将确定性模型的计算结果与SAA模型的计算结果进行了对比，并对各模型的计算规模、速度和优化效果进行了评述。

各级调度机构的职能不同，对优化调度精度的要求不同，各个调度中心的编译环境也不同。如果能快速有效地找到一个好的次优解（例如对风电随机性的适应性较强的解），在实际工程应用中可能更有价值。此类问题的建模方法主要集中在机会约束规划场景法和模糊建模上。目前尚无兼顾目标函数精度、计算速度和常规机组出力对风电随机性的适应性的建模方法。

因此，未来需要对建模方法进行进一步的研究，使其能够很好地适应风电的随机性，同时平衡计算模型的准确性和计算速度。

4 随机经济调度（Matlab实现）

4.1 代码前的准备

安装：CVX包，方法如下：

首先去CVX工具包官方网站：CVX: Matlab Software for Disciplined Convex Programming | CVX Research, Inc.

工具包安装说明在Installation — CVX Users' Guide

Unpack the file anywhere you like; a directory called cvx will be created. There are two important exceptions:
(1)Do notplace CVX in Matlab’s own toolbox directory, Octave’s built-in scripts directory.(2)Do notunpack a new version of CVX on top of an old one. We recommend moving the old version out of the way, but do not delete it until you are sure the new version is working as you expect.
StartMatlabor Octave.Do not add CVX to your path by hand.
Change directories to the top of the CVX distribution, and run the cvx_setup command. For example, if you installed CVX into C\personal\cvx on Windows, type these commands：

cd C:\personal\cvx cvx_setup

我翻译成中国话就是：

1.将文件解压到你喜欢的任何地方；一个名为cvx的目录将被创建。有两个重要的例外。
(1)不要把CVX放在Matlab自己的工具箱目录下，也就是Octave的内置脚本目录下。(2)不要把新版本的CVX解压在旧版本的上面。我们建议把旧版本移开，但在确定新版本能按你的期望工作之前，不要删除它。
2.启动Matlab或Octave。不要手工添加CVX到你的路径中。
3.改变目录到CVX发行版的顶部，然后运行cvx_setup命令。例如，如果你在Windows上把CVX安装到C\personal\cvx中，键入这些命令：

cd C:\personal\cvx cvx_setup

5 Matlab代码实现

5.1 Matlab代码

clc close all clear all format short g tic % 开始时间 load=[% 一天的负荷系数 0.6700;0.6300;0.6000;0.5900;0.5900;0.6000; 0.7400;0.8600;0.9499;0.9600;0.9600;0.9500; 0.9500;0.9500;0.9300;0.9300;0.9900;1.0000; 1.0000;0.9600;0.9100;0.8300;0.7300;0.6300;]; demand=load*500; % 平均功率需求 (kW) Y=cdf('Weibull',1:10,0.3,0.15); % 生成场景( 10个场景) demand=demand*Y; % 场景10和时间24的功率需求 cvx_begin % cvx启动 time=24; % 一天24小时 scenario=10; % 场景数 number=20; % 发电机台数 variables g(number,time,scenario) % 发电机i在场景s和时间t下的发电功率。 variables s(number,time,scenario) % 场景s和时间t中启动发电机i的二进制变量。 variables u(number,time,scenario) % 在场景 s 和时间 t 中停机发电机 i 的二进制变量。 %% 发电机的可变成本 a=[0;0;0.0141;0.0141;0;0;0.0141; 0.0141;0.0527;0.0527;0.0527;0.0072; 0.0072;0.0072;0;0.3284;0.3284;0.3284; 0.3284;0.3284;];a=a'; %% 发电机的固定成本 b=[130.0000;130.0000;16.0811;16.0811;130.0000;130.0000; 16.0811;16.0811;43.6615;43.6615;43.6615;48.5804;48.5804; 48.5804;0;56.5640;56.5640;56.5640;56.5640;56.5640;];b=b'; %% 发电机空载成本 c=[400.6849;400.6849;212.3076;212.3076;400.6849;400.6849; 212.3076;212.3076;781.5210;781.5210;781.5210;832.7575; 832.7575;832.7575;0;86.3852;86.3852;86.3852;86.3852;86.3852;];c=c'; %% 发电机的最小发电量（kW） pmin=[16;16;15.2;15.2;16;16;15.2;15.2;25;25;25;69; 69;69;0;2.4;2.4;2.4;2.4;2.4;];pmin=pmin'; %% 发电机最大发电量(千瓦) pmax=[20;20;76;76;20;20;76;76;100;100;100;197;197; 197;0;12;12;12;12;12;];pmax=pmax'; %% 发电机的启动成本 SU=[1430;1430;1430;1430;1430;1430;1430;1430;1725; 1725;1725;3056;3056;3056;437;437;437;437;437;437;];SU=SU'; %% 发电机停运成本 SD=[2860;2860;2860;2860;2860;2860;2860;2860;3450; 3450;3450;6112;6112;6112;874;874;874;874;874;874;];SD=SD'; %% 目标函数( f = a * g + b * g + c * su * s + sd * u )为线性模型 for number=1:20 for time=1:24; for scenario=1:10 f(number,time,scenario)=g(number,time,scenario)*a(number)+g(number,time,scenario)*b(number)+c(number)*SU(number)*s(number,time,scenario)+SD(number)*u(number,time,scenario); end end end %% 最小化目标函数 minimize sum(sum(sum (f)/scenario)) %% 约束条件 subject to for number=1:20; for time=1:24; for scenario=1:10; sum(g(1:20,time,scenario))==demand(time,scenario); % 功率平衡 g(number,time,scenario)>=pmin(number)*s(number,time,scenario); %发电机最小出力限制 g(number,time,scenario)<=pmax(number)*s(number,time,scenario); % 发电机最大出力限制 s(number,time,scenario)<=1;s(number,time,scenario)>=0; % 启动二进制变量的约束 u(number,time,scenario)<=1;u(number,time,scenario)>=0; % 停机二进制变量的约束 s(number,time,scenario)+u(number,time,scenario)==1; % 二元变量之间关系的约束 end end end cvx_end % 优化结束 %结果 disp('发电机产量(千瓦)') toc

5.2 结果

Calling SDPT3 4.0: 29040 variables, 9600 equality constraints For improved efficiency, SDPT3 is solving the dual problem. ------------------------------------------------------------ num. of constraints = 9600 dim. of linear var = 28800 dim. of free var = 240 480 linear variables from unrestricted variable. *** convert ublk to lblk ******************************************************************* SDPT3: Infeasible path-following algorithms ******************************************************************* version predcorr gam expon scale_data NT 1 0.000 1 0 it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime ------------------------------------------------------------------- 0|0.000|0.000|9.8e+02|1.6e+02|2.5e+14| 4.072437e+11 0.000000e+00| 0:0:00| spchol 1 1 1|1.000|0.729|2.5e-05|4.2e+01|7.0e+13| 9.630226e+10 1.075389e+10| 0:0:00| spchol 1 1 2|1.000|0.986|8.8e-10|5.8e-01|1.0e+12| 9.553163e+10 3.687323e+08| 0:0:00| spchol 1 1 3|1.000|0.975|1.8e-11|1.9e-02|8.6e+10| 6.012158e+10 2.257725e+08| 0:0:00| spchol 1 1 4|0.940|1.000|6.8e-11|2.2e-03|9.1e+09| 8.519096e+09 2.206476e+08| 0:0:01| spchol 1 1 5|0.944|1.000|1.2e-11|1.1e-03|5.5e+08| 7.463306e+08 2.266807e+08| 0:0:01| spchol 1 1 6|1.000|0.297|3.1e-11|9.5e-04|2.4e+08| 4.719523e+08 2.443795e+08| 0:0:01| spchol 1 1 7|1.000|0.453|2.8e-12|6.4e-04|1.4e+08| 4.129027e+08 2.787648e+08| 0:0:01| spchol 1 1 8|1.000|0.510|1.4e-12|3.9e-04|6.2e+07| 3.629994e+08 3.035048e+08| 0:0:01| spchol 1 1 9|1.000|0.743|1.3e-12|1.5e-04|1.4e+07| 3.414131e+08 3.276410e+08| 0:0:01| spchol 1 1 10|0.991|0.541|5.8e-13|8.8e-05|6.1e+06| 3.383346e+08 3.325221e+08| 0:0:01| spchol 1 1 11|1.000|0.427|3.8e-13|5.8e-05|3.5e+06| 3.376513e+08 3.344051e+08| 0:0:01| spchol 1 1 12|1.000|0.523|1.4e-13|3.2e-05|1.7e+06| 3.373867e+08 3.358193e+08| 0:0:01| spchol 1 1 13|0.984|0.549|5.1e-14|1.7e-05|7.5e+05| 3.372686e+08 3.365766e+08| 0:0:01| spchol 1 1 14|1.000|0.459|1.7e-14|1.0e-05|4.1e+05| 3.372470e+08 3.368660e+08| 0:0:01| spchol 1 1 15|1.000|0.328|5.8e-15|7.2e-06|2.8e+05| 3.372368e+08 3.369800e+08| 0:0:01| spchol 1 1 16|1.000|0.584|1.8e-15|3.3e-06|1.2e+05| 3.372293e+08 3.371198e+08| 0:0:01| spchol 1 1 17|1.000|0.510|9.0e-16|1.8e-06|5.8e+04| 3.372251e+08 3.371724e+08| 0:0:01| spchol 1 1 18|0.928|0.778|4.7e-16|3.9e-07|1.3e+04| 3.372240e+08 3.372121e+08| 0:0:01| spchol 1 1 19|1.000|0.090|8.2e-16|3.5e-07|1.2e+04| 3.372242e+08 3.372131e+08| 0:0:01| spchol 1 1 20|1.000|0.795|2.5e-15|7.3e-08|2.6e+03| 3.372238e+08 3.372214e+08| 0:0:01| spchol 1 1 21|0.987|0.933|1.3e-14|4.9e-09|1.7e+02| 3.372237e+08 3.372235e+08| 0:0:01| spchol 2 2 22|0.988|0.988|1.4e-15|6.2e-11|2.1e+00| 3.372237e+08 3.372237e+08| 0:0:01| stop: max(relative gap, infeasibilities) < 1.49e-08 ------------------------------------------------------------------- number of iterations = 22 primal objective value = 3.37223657e+08 dual objective value = 3.37223655e+08 gap := trace(XZ) = 2.14e+00 relative gap = 3.17e-09 actual relative gap = 2.88e-09 rel. primal infeas (scaled problem) = 1.41e-15 rel. dual " " " = 6.16e-11 rel. primal infeas (unscaled problem) = 0.00e+00 rel. dual " " " = 0.00e+00 norm(X), norm(y), norm(Z) = 1.5e+06, 1.7e+03, 1.7e+03 norm(A), norm(b), norm(C) = 4.9e+03, 8.0e+06, 4.8e+03 Total CPU time (secs) = 1.29 CPU time per iteration = 0.06 termination code = 0 DIMACS: 4.4e-14 0.0e+00 1.4e-09 0.0e+00 2.9e-09 3.2e-09 ------------------------------------------------------------------- ------------------------------------------------------------ Status: Solved Optimal value (cvx_optval): +3.23158e+07 发电机产量(千瓦) 时间已过 63.485817 秒。 >>