电力市场的交易及阻塞管理模型
电子科技大学
指导老师： 杜鸿飞
参赛队员： 张 帆
任建忠
谭利忠
2004年9月20日
电力市场的交易及阻塞管理模型
摘 要
本文遵循电力市场交易规则和阻塞管理原则，基于我们提出的阻塞费用计算规则，建立模型，通过求解模型，制定出对应一定负荷预报的出力分配方案。
为了确定各条线路的有功潮流值关于各机组出力的近似表达式，对实验数据进行多元线性回归拟合，并检验了回归方程的显著性。我们分析出在现行的电力市场规则下，发电商供电的清算价往往高于相应的报价，故引进投机因子来描述这种额外利润。从简明、合理、公平出发，提出以投机因子在调整后不减的原则，确定了阻塞费用计算方法。
针对电力交易原则下的相应机组出力分配预案的问题，我们设计相应分配预案的算法，通过编程得到相应的清算价和分配预案。
针对阻塞管理出力分配调整方案的问题。将调整方案转化为带约束的目标最优化问题。通过分配预案计算各线路的潮流值，并且与线路限值及其裕度范围进行比较，判断采取何种调整方式。特别是使用裕度输电的时候，问题的出现了两个优化目标。将多目标优化转化单目标优化，从而简化了求解难度。通过编程求解最优化问题，得到了机组出力最终方案：
在负荷需求预报下，出力分配预案为：150、79、180、、125、140、95、（这里以及以下单位均为MW），清算价为303元/MWh。经过出力分配方案的调整，达到无阻塞输电，各机组出力为：150、87. 875、228、、152、、、117，阻塞费用为20406元。在负荷需求预报下，出力分配预案为：150、81、、、135、150、、117，清算价为356元/MWh；经过阻塞管理，在安全裕度范围内输电，各机组出力为：、、228、、152、、、117，阻塞费用为10358元。
问题的重述
我国电力系统的市场化改革正在积极、稳步地进行。2003年3月国家电力监管委员会成立，2003年6月该委员会发文列出了组建东北区域电力市场和进行华东区域电力市场试点的时间表，标志着电力市场化改革已经进入实质性阶段。可以预计，随着我国用电紧张的缓解，电力市场化将进入新一轮的发展，这给有关产业和研究部门带来了可预期的机遇和挑战。
电力从生产到使用的四大环节——发电、输电、配电和用电是瞬间完成的。我国电力市场初期是发电侧电力市场，采取交易与调度一体化的模式。电网公司在组织交易、调度和配送时，必须遵循电网“安全第一”的原则，同时要制订一个电力市场交易规则，按照购电费用最小的经济目标来运作。市场交易-调度中心根据负荷预报和交易规则制订满足电网安全运行的调度计划——各发电机组的出力（发电功率）分配方案；在执行调度计划的过程中，还需实时调度承担AGC（自动发电控制）辅助服务的机组出力，以跟踪电网中实时变化的负荷。
设某电网有若干台发电机组和若干条主要线路，每条线路上的有功潮流（输电功率和方向）取决于电网结构和各发电机组的出力。电网每条线路上的有功潮流的绝对值有一安全限值，限值还具有一定的相对安全裕度（即在应急情况下潮流绝对值可以超过限值的百分比的上限）。如果各机组出力分配方案使某条线路上的有功潮流的绝对值超出限值，称为输电阻塞。当发生输电阻塞时，需要研究如何制订既安全又经济的调度计划。
本文所要解决的问题如下：
某电网有8台发电机组，6条主要线路，在附件一中，表1和表2中的方案0给出了各机组的当前出力和各线路上对应的有功潮流值，方案1~32给出了围绕方案0的一些实验数据，试用这些数据确定各线路上有功潮流关于各发电机组出力的近似表达式。
设计一种简明、合理的阻塞费用计算规则，除考虑上述电力市场规则外，还需注意：在输电阻塞发生时公平地对待序内容量不能出力的部分和报价高于清算价的序外容量出力的部分。
假设下一个时段预报的负荷需求是，表3、表4和表5分别给出了各机组的段容量、段价和爬坡速率的数据，试按照电力市场规则给出下一个时段各机组的出力分配预案。
按照表6给出的潮流限值，检查得到的出力分配预案是否会引起输电阻塞，并在发生输电阻塞时，根据安全且经济的原则，调整各机组出力分配方案，并给出与该方案相应的阻塞费用。
假设下一个时段预报的负荷需求是，重复3~4的工作。
问题分析
求解各机组的出力与各线路的有功潮之间的关系，是解决其他问题的基础。他们之间的关系可以分析已知的数据得出。
在电力的交易过程遵循电力市场交易规则（详见附件二），即按照购电费用最小的经济目标来运行，但是受到各机组的当前出力和爬坡速率的约束。
当出现输电阻塞时，实施阻塞管理会直接造成一些序内容量不能出力，而一些序外容量（序内容量指通过竞价取得发电权的发电容量；序外容量指竞价中未取得发电权的发电容量）要在低于其对应报价的清算价上出力。显然发电商的利益会受到损害，网方该给予一定的经济补偿——阻塞费用。阻塞费用的设置应该按市场经济的运作规则来进行，公平的对待序内容量不出力的部分，以及序外容量出力的部分。
在一个时段内，电网公司的纯收益可以表示如下：
纯收益＝电网公司的总收入－总支出费用（购电费用和阻塞费用）
一个时段内，电力的调配是以时段末的预计需求功率来计算的，用户的用电总量基本不随发电机组的出力的调配而改变，所以电网公司的总收入（向用户收取的用电费）基本上是一定的，要获取最高的纯收益，就得使时段内的总支出费用最小。
阻塞管理也是一个优化问题——在保证安全的情况下追求经济效益（购电费用和阻塞费用的总和最少）。由阻塞管理规则可以知道，阻塞管理主要分三步。第一步，调整各机组的出力方案使得输电阻塞消除，问题的约束条件是各电线的功率潮流不超过限值；第二步，如果第一做不到，则使用线路的安全裕度输电，这里不但要最求经济利益，同时也应使每条线路上潮流的绝对值超过限值的百分比尽量小——双目标规划，问题的约束条件变为各线路的有功潮流小于各线路的允许的最大功率潮流值。第三步，如果无论怎样分配各机组的出力都无法使每条线路的潮流绝对值超过限值的百分比小于相对安全裕度，则必须在用电侧拉闸限电。
本文所要解决的问题中，电力的交易和阻塞管理都可以归结为带约束条件的优化问题。
基本假设
1、 交易－调度中心对下一个时段的负荷需求预报是准确的；
2、 用户的用电量是随时段而阶跃的，一个时段内用电量不随机组出力的调整而改变；
3、 各机组增加出力是以最大的爬坡速率进行的，同时爬坡结束在时段末；各机组减小出力是也以最大的爬坡速率进行的，但从时段开始时刻就减小出力——简称“降先升后”。
对假设3的说明：发电机组爬坡的示意图如下：
p p p p
t t t t
（a） （b） （c） （d）
图1 发电机组爬坡示意图
这样处理考虑了三方面的原因。第一、输电安全：所有的计算都是以时段末的情况来计算的，而在时段内进行调整，当有机组按图1 (a)和(c) 的方式爬坡时，虽然在段末可以满足安全条件，但机组的出力在时段内最大，可能会造成线路的不安全。第二、时段内的总电能最少，线路的出力曲线与时间轴围成的面积实际上是这段时间出力的总电能，“降先升后”会使总电能量减少，电网商的购电费用减少；第三、简化问题。
符号变量说明
： 第个机组爬坡速率 单位：MW/分钟
： 第个机组的当前时段的出力 单位：MW
： 第个机组的预案出力 单位：MW
： 第个机组阻塞管理调整后的出力 单位：MW
： 第j条线路的有功潮流 单位：MW
： 第j条线路的有功潮流限值 单位：MW
： 各线路的有功潮流与各机组出力的关系矩阵
： 电网公司对第个机组的购电量 单位：MWh
： 清算价 单位：元/MWh
： 购电费用 单位：元
： 阻塞费用 单位：元
： 电网公司对第个发电机组补偿的费用 单位：元
： 时段的时间长度 值：15分钟
模型建立与求解
有功潮流与各发电机组出力的关系
参数估计
对附件一中的线路有功潮流值分析发现，第3条线路的有功潮流值在各种方案下均为负。这是由于负号仅仅代表了该线路上的潮流方向，而并不是指向发电机组输送能量。
分析各机组不同出力方案与各线路的对应潮流值，发现线路的潮流值与机组的出力有着线性关系。
第条线路上的功率潮流与各机组的出力关系
j＝1,2,…,6
将已知的数据带入方程
j＝1,2，…,6
k=0,1,…,32表示方案号，表示第k方案，第个机组的出力，表示第k方案，第j条线路上的功率潮流。
方便起见，以下采用矩阵叙述。
则 (1)
方程(1)是超定方程，利用最小二乘估计，对参数B求解
从而
.
2、回归显著性检验
残差平方和
k=0,1,…,32，j＝1，2，…，6
回归离差平方和
k=0,1,…,32，j＝1，2，…，6
则
计算得到的分别为：
＝，＝，＝21788，＝24424，＝，＝16029.
当显著性水平为时，，远远小于，可见回归显著。
实际上，各发电机组与各线路实际组成了一个复杂的电路网络系统。由电路的基本知识可以知道，电力系统是一个非线性系统，功率是不满足线性叠加原理的，即线路上的功率潮流与各发电机的关系不是线性的。对这种关系线性处理，是当发电机组的出力在一定的范围内波动才成立，否则会出现各发电机组均不出力，而线路上仍存在有功潮流的情况，这显然是违背自然规律的。
（二）阻塞费用的计算规则
各发电商的收入可分为两块：劳动报酬与投机回报。这里，
劳动报酬＝时段总发电量×相应的段价
投机回报＝时段总发电量×清算价－劳动报酬
投机回报量度了发电商获得的在自己估价以外的报酬。不妨以第i个发电商为例。设其序内容量（或序外容量）对应的出力量为，段价为,该时段清算价为MCP。在调整后，其出力变为，相应的段价为。如图2所示：
图2 出力方案调整前后供电量与对应报价的对比示意图
其中，左图示意序内容量不能出力，右图示意序外容量在低于对应报价的清算价上出力。
调整出力分配方案之前出力为，对应发电量为 收入中投机回报为。由于阻塞做出的调整使得出力、发电量、投机回报分别变为、、。定义投机因子为投机回报与发电量的比值，为调整后投机回报与发电量的比值，即
对序内容量而言，、均大于或等于0；对序外容量而言，大于或等于0，小于0。作为供电商，当然希望尽量地高，因为这意味着单位出力所带来的收益高。然而，在阻塞调整之后，这个比例一般会发生变动。我们认为，调整后这个比例的值至少应维持在调整前的水平上，即
，
此时，发电商与网方没有经济利益冲突。网方需支付的阻塞费用
否则，如果
，
则发电商与网方将产生经济利益冲突。为此，网方将不得不支付给第i家机组一定的阻塞费用，直到补偿后的修正投机因子
；
或 ，
同时，注意到网方应在安全运行的前提下尽量减少阻塞费用，由上式可知，应取
综上所述，
，
将的表达式代入，即
. (2)
这种方案通过支付发电商一定数额的阻塞费用，使发电商的投机因子恢复到调整以前，从而能很公平地对待序内容量不能出力的部分和报价高于清算价的序外容量出力的部分。
（三）电力市场交易与调度
不阻塞
阻塞
可消除
无法消除阻塞
可调整
无法调整
尽量满
足用户
图3 整个过程流程示意图
制定出力分配预案
各机组的出力分配预案是按照电力市场交易规则进行的，以满足购电费用最小这一目标。在交易时段，各机组给出自己的段容量，以及下一个时段的段价。交易中心根据每台机组的当前出力和爬坡速率，按照段价从低到高的选取机组的段容量或其部分，直到所选之和等于预报的负荷。所选的结果构成下一个时段的机组的出力分配预案。将最后一个被选入的段价（最高段价）作为该时段的清算价。根据此规则，可利用计算机编程求解对应输入当前出力、爬坡速率、下一时段负荷需求，相应输出各机组的出力分配。简略的算法流程框图如下（程序见附件）
图4 出力预案制定示意图
判断是否输电阻塞
根据机组出力预案，通过拟合方程计算各条线路上的有功潮流的估计值：
将与各线路的有功潮流的限值比较，如果存在，使得 ，则第j条线路发生输电阻塞，应实行阻塞管理，调整发电机组的出力方案；如果对于任意的j，都使得，则该方案不会导致阻塞，按出力预案输电。
3、调整机组出力，消除输电阻塞
当发生输电阻塞时，调整发电机组的出力方案，消除输电阻塞。电网公司希望该时段的购电费用和阻塞费用之和最小。
调整发电机组的出力，这个时段电网公司对第个机组的购电量情况
所以购电费用
电网公司给第发电机组的补偿（阻塞费用）
其中表示交易时对应的段价，表示调整后时对应的段价
则电网公司总共支出的阻塞费用
电网公司支出的总费用为
但这时候各机组的出力同样受到其当前出力与爬坡速率的限值
其中分钟
将问题表示成标准的带约束条件的最优化问题
(3)
.
其中表示机组出力预案的出力的总和。
使用安全裕度输电
如果不能通过出力方案的调整来消除输电阻塞，则考虑使用线路的安全裕度输电，以避免拉闸限电给用户造成极大的不便。
从“安全第一”的原则出发，首先应尽量让每条线路上的潮流的绝对值超过限制的百分比尽量小。这里我们定义第j条线路上的超限比例为，即
，
出于安全，上述比例应不超过相对安全裕度，即
，
对每个相对安全裕度，都有其饱和度，即 占的百分比。定义饱和度
，
于是，要使每条线路上潮流的绝对值超过限值的百分比尽量小，只需让尽量小。即
。
另一方面，从支出的总费用最小考虑，应有
于是，安全裕度输电的阻塞管理可归结为以下多目标规划问题：
(4)
.
为求解以上多目标规划问题，考虑将两个目标合理地统一为一个目标，从而得到问题的合理解答。
首先，分别求解两个单目标规划问题：
s．t．
s．t．
得到两个单目标规划问题的目标函数的最小值、。
我们希望能对两个目标很好的统一，得到一个新的目标函数，它能同时较好地反映这两个目标。于是通过构造新的目标函数Q将多目标规划转化为单目标规划：
(5)
是在当前限制下可能的最小值。如果越接近（表示对的优化程度越高），对应的Q越小；同理，对也是这样。因此，Q确实能较好地同时反映这两个目标。同时，这样做还有助于将不同量纲的量归一化，做加权求和。综上所述，做出以下规划：
(6)
.
其中分别表示、在总目标中所占的比重。它可以由政府的政策、电网公司的决策决定。在这里，我们考虑到“安全第一”，同时考虑经济因素，我们取。
拉闸限电
如果无论怎样分配机组出力都无法使每条线路上的潮流绝对值超过限值的百分比小于相对安全裕度，则必须在用电侧拉闸限电。这时，不仅要考虑电网公司的总支出尽量小、尽量安全，还要考虑尽量满足用户。类似前面的分析可以得到下面的关系
.
同安全裕度输电的处理方法，将多目标规划转化为单目标规划
(7)
.
（四）针对两种预报负荷，求解模型
按照电力市场交易与调度的流程，我们通过编制Matlab程序（详见附件三）。通过上文可以看出，对于每一个调整阶段，我们都提出了其目标规划式，并将规划式中的多目标规划转换为了单目标规划。转化后的目标函数为非线性多元函数，而且目标式中存在线性约束，这些规划问题都是典型的带约束的多元函数非线性优化问题。我们使用当今比较成熟可靠的Matlab软件优化工具箱的有约束的非线性多元函数寻优函数fmincon进行参数的优化求解。
在约束条件下，如果有解则通过参数寻优搜索最优的出力方案；如果没有可行解，即不存在一个方案能满足该调整阶段的全部约束条件，这时我们就进入下一个调整阶段用相同的求解方法对下一阶段的新的目标规划式进行优化求解（优化程序请见附件）
我们考虑到这是一个比较复杂的多元函数，我们通过该方法求出的最优的函数值有可能只是局部最优解，而不一定是全局最优解。为了最大可能的得到调整方案的全局最优解，我们通过使寻优函数从很多不同的初始值开始迭代，再对其每次迭代选出的最优解进行比较，选出比较下的最优值与最优参数。这样就保证了每次优化后选出全局相对优解的可能性相当的大。
计算出的结果如下：
当下一个时段的预报负荷是时
表1 各机组出力情况与费用
机组
1
2
3
4
5
6
7
8
出力预案（MW）
150
79
180
125
140
95
阻塞管理调整(MW)
150
228
152
117
购电费用（元）
10123
15453
10491
阻塞费用（元）
0
6273
0
10214
0
0
0
清算价303元/MWh，总费用90227元，其中购电费69821元，阻塞费用20406元。
表2 线路的潮流值
线路
1
2
3
4
5
6
预案分配线路的潮流（MW）
141
调整后的线路潮流(MW)
165
-155
132
超过的限值的百分比
可以看出通过调整可以进行无阻塞送电。
2、当下一个时段的预报需求是时
表3 各机组出力情况与费用
机组
1
2
3
4
5
6
7
8
出力预案(MW)
150
81
135
150
117
阻塞管理调整(MW)
228
152
117
购电费用(元)
11624
18156
12327
11596
阻塞费用（元）
0
0
0
6925
0
0
0
清算价为356元/MWh，总费用95040元，其中购电费84682元，阻塞费用10358元。
表4 线路的潮流值
线路
1
2
3
4
5
6
预案线路的潮流(MW)
调整后的线路潮流(MW)
绝对值超过限值的百分比
％
0
0
0
％
％
结果分析
如果线路的限值足够大，那么电网商就不会因为输电阻塞而给发电商额外的补偿。提高线路的限值，消除输电阻塞，当预报需求是时，电网商将减少20406元的阻塞费用；当预报需求是时，电网商将减少10358元的阻塞费用。线路的潮流限值使得电网公司的利润减少了，如果对电网改造增加电网潮流限值，电网公司的利润将大大增加，同时也会增加输电安全性。
七、模型评价
各线路的潮流值与各机组的出力关系的线性关系是有一定的范围的，要在更大范围的描述二者的关系，就需要分析电网的物理结构。建立二者更加准确关系模型，对电网商调整机组出力，和保障电网的安全，提高电网商的利润将有重要的意义。
模型从合理、公平的角度，设计出了一种阻塞费用规则，使发电商的损失得到补偿，并且使电网商可以接受，这种想法可以推广到解决多方的经济利益冲突，有一定普适性。
同时给电网公司提供了一种解决阻塞管理的量化方法，使得阻塞管理更加的科学，在遵守输电管理原则的条件下，提高电网公司的利润。但是模型中一些权重系数是不太客观的，有待改进。
参考文献：
[1] 周全仁，张清益 电网分析与发电计划 湖南科技出版社 1996
[2] 姜启源，谢金星, 叶俊 数学模型（第三版） 高等教育出版社 2003. 8
[3] 刘承平 数学建模方法 高等教育出版社
附件清单：
附件一：问题的有关的已知数据
附件二：电力市场交易规则
附件三：有关的求解程序
附 件
一、附件一：问题的有关的已知数据
表1 各机组出力方案 （单位：兆瓦，记作MW）
方案\机组
1
2
3
4
5
6
7
8
0
120
73
180
80
125
125
90
1
73
180
80
125
125
90
2
73
180
80
125
125
90
3
73
180
80
125
125
90
4
73
180
80
125
125
90
5
120
180
80
125
125
90
6
120
180
80
125
125
90
7
120
180
80
125
125
90
8
120
180
80
125
125
90
9
120
73
80
125
125
90
10
120
73
80
125
125
90
11
120
73
80
125
125
90
12
120
73
80
125
125
90
13
120
73
180
125
125
90
14
120
73
180
125
125
90
15
120
73
180
125
125
90
16
120
73
180
125
125
90
17
120
73
180
80
125
90
18
120
73
180
80
125
90
19
120
73
180
80
125
90
20
120
73
180
80
125
90
21
120
73
180
80
125
90
22
120
73
180
80
125
90
23
120
73
180
80
125
90
24
120
73
180
80
125
90
25
120
73
180
80
125
125
90
26
120
73
180
80
125
125
90
27
120
73
180
80
125
125
90
28
120
73
180
80
125
125
90
29
120
73
180
80
125
125
30
120
73
180
80
125
125
31
120
73
180
80
125
125
32
120
73
180
80
125
125
表2 各线路的潮流值（各方案与表1相对应，单位：MW）
方案\线路
1
2
3
4
5
6
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
119
20
21
22
23
24
25
26
141
27
28
29
30
31
32
表3 各机组的段容量 （单位：MW）
机组\段
1
2
3
4
5
6
7
8
9
10
1
70
0
50
0
0
30
0
0
0
40
2
30
0
20
8
15
6
2
0
0
8
3
110
0
40
0
30
0
20
40
0
40
4
55
5
10
10
10
10
15
0
0
1
5
75
5
15
0
15
15
0
10
10
10
6
95
0
10
20
0
15
10
20
0
10
7
50
15
5
15
10
10
5
10
3
2
8
70
0
20
0
20
0
20
10
15
5
表4 各机组的段价（单位：元/兆瓦小时，记作元/MWh）
机组\段
1
2
3
4
5
6
7
8
9
10
1
-505
0
124
168
210
252
312
330
363
489
2
-560
0
182
203
245
300
320
360
410
495
3
-610
0
152
189
233
258
308
356
415
500
4
-500
150
170
200
255
302
325
380
435
800
5
-590
0
116
146
188
215
250
310
396
510
6
-607
0
159
173
205
252
305
380
405
520
7
-500
120
180
251
260
306
315
335
348
548
8
-800
153
183
233
253
283
303
318
400
800
表5 各机组的爬坡速率 （单位：MW/分钟）
机组
1
2
3
4
5
6
7
8
速率
1
2
表6 各线路的潮流限值（单位：MW）和相对安全裕度
线路
1
2
3
4
5
6
限值
165
150
160
155
132
162
安全裕度
13%
18%
9%
11%
15%
14%
二、附件二：电力市场交易规则
1. 以15分钟为一个时段组织交易，每台机组在当前时段开始时刻前给出下一个时段的报价。各机组将可用出力由低到高分成至多10段报价，每个段的长度称为段容量，每个段容量报一个价（称为段价），段价按段序数单调不减。在最低技术出力以下的报价一般为负值，表示愿意付费维持发电以避免停机带来更大的损失。
2. 在当前时段内，市场交易-调度中心根据下一个时段的负荷预报,每台机组的报价、当前出力和出力改变速率，按段价从低到高选取各机组的段容量或其部分（见下面注释），直到它们之和等于预报的负荷，这时每个机组被选入的段容量或其部分之和形成该时段该机组的出力分配预案（初始交易结果）。最后一个被选入的段价（最高段价）称为该时段的清算价，该时段全部机组的所有出力均按清算价结算。
注释：
每个时段的负荷预报和机组出力分配计划的参照时刻均为该时段结束时刻。
机组当前出力是对机组在当前时段结束时刻实际出力的预测值。
假设每台机组单位时间内能增加或减少的出力相同，该出力值称为该机组的爬坡速率。由于机组爬坡速率的约束，可能导致选取它的某个段容量的部分。
为了使得各机组计划出力之和等于预报的负荷需求，清算价对应的段容量可能只选取部分。
市场交易-调度中心在当前时段内要完成的具体操作过程如下：
监控当前时段各机组出力分配方案的执行，调度AGC辅助服务，在此基础上给出各机组的当前出力值。
作出下一个时段的负荷需求预报。
根据电力市场交易规则得到下一个时段各机组出力分配预案。
计算当执行各机组出力分配预案时电网各主要线路上的有功潮流，判断是否会出现输电阻塞。
三、附件三：有关的求解程序
附件程序清单：
1. 第一问求解程序
第二问求阻塞费用的函数
第三问按电力市场规则给出相应需求的机组出力分配预案
第四问阻塞管理下按照优化指标和相应约束优化求解出分配方案
时段总费用函数
由当前机组出力和下时段机组出力预案求下时段机组发电总电量函数
衡量出力求得潮流值超过限值的函数
多目标划单目标后的优化目标函数
程序运行所必须的数据文件清单
题目表一数据 （各机组出力方案及其实验方案）
题目表二数据 （表一方案对应的各线路的潮流值）
题目表三数据 （各机组段容量）
题目表四数据 （各机组的段价）
%%%%%%%%%%%%%%%%
%% 读取数据，
%% (如果需要试运行,
%% 请将 拷入C盘根目录，
%% 或更改下面的目录)
clear
fid=fopen('c:\','rt');
data1=myread(fid);
%% data1为各机组出力方案及其实验方案
fid=fopen('c:\','rt');
data2=myread(fid);
%% data2为各方案相应的各线路的有功潮流值
%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
%% 多元线性回归
X=data1;
X(:,9)=ones(33,1);%常数项
for k=1:6
[A(:,k),A1,R,R1,stats(k,:)]=regress(data2(:,k),X,);
end
% A为有功潮流与各发电机出力的线性系数矩阵A. stats为线性相关性检验矩阵
sprintf('有功潮流与各发电机出力的线性关系矩阵A为:')
A
function fee=getfee1(origin_p,now_p,price,qingsuan,data3)
% fee:阻塞费用
% origin_p:市场决定的分配后各机组的下时段末发电功率
% now_p:调整后分配后各机组的下时段末发电功率
% price:发电方报价
% qingsuan:清算价
% data3: 段容量矩阵
for k=1:8
for kk=1:10
if now_p(k)<=data3(k,kk)
nowprice(k)=price(k,kk);
break;
end
end
for kk=1:10
if origin_p(k)<=data3(k,kk)
oriprice(k)=price(k,kk);
break;
end
end
end
%% 获得调整前后发电方报价
%% oriprice为调整前相应段报价 nowprice为调整后相应报价
origin=getelec(origin_p);
now=getelec(now_p);
%% origin市场决定的分配后各机组的下时段末总发电量
%% now:调整后分配后各机组的下时段末发电量
for k=1:8
ori_extra(k)=(qingsuan-oriprice(k)).*origin(k);% 调整前相应额外收入
now_extra(k)=(qingsuan-nowprice(k)).*now(k);% 调整后相应额外收入
if ori_extra(k)/origin(k)>now_extra(k)/now(k)
fee(k)=now(k)*ori_extra(k)./origin(k)-now_extra(k);
else
fee(k)=0;
end
end
%%%%%%%%%%%%%%%%
%% 读取数据，
%% (如果需要试运行,
%% 请将 拷入C盘根目录，
%% 或更改下面的目录)
fid=fopen('c:\','rt');
data1=myread(fid);
fid=fopen('c:\','rt');
data3=myread(fid);
fid=fopen('c:\','rt');
d_price=myread(fid);
%% data1为各机组出力方案及其实验方案
ac=[ 1 2 ];%%ac为爬坡速率
l=ac*15;
now=data1(1,:);
for k=1:8
p(k,1)=now(k)-l(k);
p(k,2)=now(k)+l(k);
end
%% p为下个时段各机组可能到达的出力区间
for k=1:8
for kk=2:10
data3(k,kk)=data3(k,kk)+data3(k,kk-1);
end
end
for k=1:8
for kk=10:-1:2
if data3(k,kk)==data3(k,kk-1);
data3(k,kk)=0;
end
end
end
%%读入并处理段容量数据
money=zeros(8,10);
for k=1:8
for kk=1:10
if data3(k,kk)>p(k,1) & data3(k,kk)<=p(k,2) & data3(k,kk)~=0
money(k,kk)=d_price(k,kk);
end
end
end
% 求可行解的段价矩阵
need=;% need为下一时段预报。题3数据为。 题5数据为
%%
max_d=[10 10 8 6 10 8 6 7];%各机组下时段最大可行段
money= ...
[
0 0 124 0 0 252 0 0 0 489;
0 0 0 0 245 300 320 0 0 495;
0 0 152 0 233 0 308 356 0 0;
0 0 170 200 255 302 0 0 0 0;
0 0 0 0 188 215 0 310 396 510;
0 0 159 173 0 252 305 380 0 0;
0 120 180 251 260 306 0 0 0 0;
-800 0 183 0 253 0 303 0 0 0
];
%% money为计算出来的可行段的段价矩阵（通过计算出的Money调整得到）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 按电力市场规则给出下一时段出力分配预案
start=[3 5 3 3 5 3 2 1];
total=zeros(1,34);
last_price=total;
index_x=total;
index_y=total;
temp=0;
for k=1:8
now_min(k)=money(k,start(k));
end
for m=1:34
[the_min,index]=min(now_min);
index_x(m)=index;
index_y(m)=start(index);
last_price(m)=money(index_x(m),index_y(m));
total(m)=data3(index_x(m),index_y(m));
start(index)=start(index)+1;
if start(index)>max_d(index)
now_min(index)=1000;
start(index)=max_d(index);
total(m)=p(index,2);
else
for kk=start(index):max_d(index)
if money(index,kk)~=0
now_min(index)=money(index,kk);
start(index)=kk;
break;
end
end
end
end
for m=1:34
power=zeros(1,8);
for k=1:m
power(index_x(k))=total(k);
end
temp=sum(power);
if temp>=need
power(index_x(m))=power(index_x(m))-(temp-need);
break;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
qingsuan=last_price(m);
sprintf('清算价为 %d',qingsuan)%% qingsuan为清算价
for k=1:8
sprintf('第%d组机组出力%d MW',k,power(k))%% power(k)为第k组发电机的出力预案
end
%%%%%%%%%%%%
%% 由于需要数据，故请事先运行与文件
%%%%%%%%%%%
power1=[150 79 180 125 140 95 ];
power2=[150 81 135 150 117];
qingsuan1=303;qingsuan2=356;
%% power1,power2为问题3，5按市场规则的分配预案
yudu=[ ];%各线路裕度矩阵
C0=[165 150 160 155 132 162];%各线路限值矩阵
AA=A(1:8,:)';
C=(C0-A(9,:))';
lb=p(:,1);
ub=p(:,2);
Aeq=ones(1,8);
need=;
beq=need;
global origin_p price qs d_data;
origin_p=power1;price=d_price;qs=qingsuan1;d_data=data3;
x0=now;
x1=now+800;
[x,fval,exitflag,output] = fmincon('get_total_fee1',x1,AA,C,Aeq,beq,lb,ub)
%% 问题3的优化 x为调整后最佳出力分配预案
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 问题4的优化
need=;
beq=need;
CC=(C0.*yudu-A(9,:))';
origin_p=power2;price=d_price;qs=qingsuan2;d_data=data3;
global A C0 fval1 fval2 yudu
opt=optimset('fmincon');
opt=optimset(opt,'MaxFunEvals','500*numberOfVariables','LargeScale','on','MaxIter',500);
[x1,fval1,exitflag1,output]=fmincon('get_total_fee1',x0,AA,CC,Aeq,beq,lb,ub,[],opt)
[x2,fval2,exitflag2]=fmincon('fun_xian',x0,AA,CC,Aeq,beq,lb,ub);
%% fval1,fval2 为单目标最优值
[x3,fval3,exitflag3]=fmincon('muti_goal1',x0,AA,CC,Aeq,beq,lb,ub);
%% x3为调整后最佳出力分配预案
function total_fee=get_total_fee1(now_p)
global origin_p price qs d_data;
% total_fee:阻塞费用与购电费用之和
% origin_p:市场决定的分配后各机组的下时段末发电功率
% now_p:调整后分配后各机组的下时段末发电功率
% price:发电方报价
% qs:清算价
% d_data: 段容量矩阵
fee1=getfee1(origin_p,now_p,price,qs,d_data);
fee2=qs.*getElec(now_p);
total_fee=sum(fee1)+sum(fee2);
function elec=getElec(power)
% 此为得到每时段各机组提供的电量总数(elec)
% power:市场决定的分配后各机组的下时段末发电功率
pre=[120 73 180 80 125 125 90];% 当前时段出力
ac=[ 1 2 ];%%ac为爬坡速率
for k=1:8
if pre(k)~=power(k);
t=abs(pre(k)-power(k))./ac(k);
if pre(k)>power(k)
elec(k)=(pre(k)+power(k)).*t/2+power(k)*(15-t);
else
elec(k)=(power(k)+pre(k)).*t/2+pre(k)*(15-t);
end
else
elec(k)=power(k)*15;
end
end
elec=elec/60;
function xianzhi=fun_xian(now_p)
global A C0 yudu
current=abs([now_p,1]*A);
xianzhi=0;
for k=1:6
if current(k)>C0(k)
xianzhi=xianzhi+((current(k)-C0(k))/((yudu(k)-1)*C0(k))).^2;
end
end
function mg=muti_goal1(now_p)
w=[ ];% 权系数
global fval1 fval2
goal2=(fun_xian(now_p)-fval2)/fval2;
goal1=(get_total_fee1(now_p)-fval1)/fval1;
mg=w(1)*goal1+w(2)*goal2
PAGE
PAGE 23
