电力市场的输电阻塞模型带matlab程序-学位论文.doc

电力市场的输电阻塞模型

摘要:
本文针对第一问，建立了多元线性均值回归模型，用于描述各线路上的有功潮流与各发电机组出力的关系。模型中利用均值定理，运用matlab软件求得解，此模型的误差较小，较符合实际情况。
对于第二问，我们引入基本补偿价概念，利用限上和限下的方法设计阻塞费用计算规则，在阻塞费用计算规则里面我们认为序外容量在低于对应报价的清算价上出力的时候，网方已经按清算价给予了发电商一部分的费用，这部分费用不作为经济补偿，即其是不包含在阻塞费用中，对于序外容量的经济补偿我们就用对应段容量的报价和清算价之差作为补偿价格来计算阻塞费用。
对于第三问，在考虑爬坡速率和负荷需求的前提下，根据电力市场规则，考虑最小购电成本为目标函数，利用最优化的方法，我们得到了各机组的出力分配预案（如下表所示）。
机组
1
2
3
4
5
6
7
8
出力
150
79
180
99.5
125
140
95
113.9
此时各机组所对应的潮流值（如下表所示）
线路
1
2
3
4
5
6
潮流值
173.3157
96.8095
-150.5051
101.8151
136.8380
168.5319
对于第四问，利用第一问求得的结果，检验第三问各机组的出力分配预案是否会引起输电阻塞，如果引起输电阻塞，根据输电阻塞管理原则并考虑经济原则，用计算机编程搜索调整预案，得出阻塞费用为405.8004元，新的方案见下表
机组
1
2
3
4
5
6
7
8
出力
153
88
228
90.1775
152
95
60.1
116.1225
此时各线路所对应的潮流值（如下表所示）
线路
1
2
3
4
5
6
潮流值
165
100.6408
-154.8553
105.5524
131.5169
159.5725
对于第五问，重复三，四问的工作，我们得出方案（如下表所示），阻塞费用为94.0377元。
机组
1
2
3
4
5
6
7
8
出力
153
88
228
99.5
152
155
60.3
117
此时各线路对应的潮流值（如下表所示）
线路
1
2
3
4
5
6
潮流值
173.4619
94.1143
-154.9968
104.9631
135.3452
160.4779

问题的重述：

电网公司在组织交易，调度和配送时,必须遵循电网“安全第一”的原则，同时制定一个电力市场交易规则，按照购电费用最小的经济目标来运作。市场交易-调度中心根据负荷预报和交易规则制定满足电网安全运行的调度计划――各发电机组的处理(发电功率)分配方案；在执行调度计划的过程中，还需实时调度承担AGC(自动发电控制)辅助服务的机组出力，以跟踪电网中实时变化的负荷。
设某电网有若干台发电机组和若干条主要线路，每条线路上的有功潮流（输电功率和方向）取决于电网结构和各发电机组的出力。电网每条线路上的有功潮流的绝对值有一安全限值，限值还具有一定的相对安全裕度（即在应急情况下潮流绝对值可以超过限值的百分比的上限）。如果各机组出力分配方案使某条线路上的有功潮流的绝对值超出限值，称为输电阻塞。当发生输电阻塞时，需要研究如何制订既安全又经济的调度计划。
电力市场交易规则：
1. 以15分钟为一个时段组织交易，每台机组在当前时段开始时刻前给出下一个时段的报价。各机组将可用出力由低到高分成至多10段报价，每个段的长度称为段容量，每个段容量报一个价（称为段价），段价按段序数单调不减。在最低技术出力以下的报价一般为负值，表示愿意付费维持发电以避免停机带来更大的损失。
2. 在当前时段内，市场交易-调度中心根据下一个时段的负荷预报,每台机组的报价、当前出力和出力改变速率，按段价从低到高选取各机组的段容量或其部分（见下面注释），直到它们之和等于预报的负荷，这时每个机组被选入的段容量或其部分之和形成该时段该机组的出力分配预案（初始交易结果）。最后一个被选入的段价（最高段价）称为该时段的清算价，该时段全部机组的所有出力均按清算价结算。
注释：
（a） 每个时段的负荷预报和机组出力分配计划的参照时刻均为该时段结束时刻。
（b） 机组当前出力是对机组在当前时段结束时刻实际出力的预测值。
（c） 假设每台机组单位时间内能增加或减少的出力相同，该出力值称为该机组的爬坡速率。由于机组爬坡速率的约束，可能导致选取它的某个段容量的部分。
（d） 为了使得各机组计划出力之和等于预报的负荷需求，清算价对应的段容量可能只选取部分。
市场交易-调度中心在当前时段内要完成的具体操作过程如下：
1、 监控当前时段各机组出力分配方案的执行，调度AGC辅助服务，在此基础上给出各机组的当前出力值。
2、 作出下一个时段的负荷需求预报。
3、 根据电力市场交易规则得到下一个时段各机组出力分配预案。
4、 计算当执行各机组出力分配预案时电网各主要线路上的有功潮流，判断是否会出现输电阻塞。如果不出现，接受各机组出力分配预案；否则，按照如下原则实施阻塞管理：
输电阻塞管理原则：
（1） 调整各机组出力分配方案使得输电阻塞消除。
（2） 如果（1）做不到，还可以使用线路的安全裕度输电，以避免拉闸限电（强制减少负荷需求），但要使每条线路上潮流的绝对值超过限值的百分比尽量小。
（3） 如果无论怎样分配机组出力都无法使每条线路上的潮流绝对值超过限值的百分比小于相对安全裕度，则必须在用电侧拉闸限电。
（4） 当改变根据电力市场交易规则得到的各机组出力分配预案时，一些通过竞价取得发电权的发电容量（称序内容量）不能出力；而一些在竞价中未取得发电权的发电容量（称序外容量）要在低于对应报价的清算价上出力。因此，发电商和网方将产生经济利益冲突。网方应该为因输电阻塞而不能执行初始交易结果付出代价，网方在结算时应该适当地给发电商以经济补偿，由此引起的费用称之为阻塞费用。网方在电网安全运行的保证下应当同时考虑尽量减少阻塞费用。

现要解决的问题如下：（表格见附录1）
1. 某电网有8台发电机组，6条主要线路，表1和表2中的方案0给出了各机组的当前出力和各线路上对应的有功潮流值，方案1~32给出了围绕方案0的一些实验数据，试用这些数据确定各线路上有功潮流关于各发电机组出力的近似表达式。
2. 设计一种简明、合理的阻塞费用计算规则，除考虑上述电力市场规则外，还需注意：在输电阻塞发生时公平地对待序内容量不能出力的部分和报价高于清算价的序外容量出力的部分。
3. 假设下一个时段预报的负荷需求是982.4MW，表3、表4和表5分别给出了各机组的段容量、段价和爬坡速率的数据，试按照电力市场规则给出下一个时段各机组的出力分配预案。
4. 按照表6给出的潮流限值，检查得到的出力分配预案是否会引起输电阻塞，并在发生输电阻塞时，根据安全且经济的原则，调整各机组出力分配方案，并给出与该方案相应的阻塞费用。
5. 假设下一个时段预报的负荷需求是1052.8MW，重复3~4的工作。
问题的分析：
1 初步分析题意，根据表1和表2所给的数据，我们利用数理统计知识和线性回归分析方法，拟建立多元线性均值回归模型，推算出各线路上有功潮流关于各发电机组出力的近似表达式；
2 设计简明，合理的阻塞费用的计算规则，关键在于体现网方是怎样充分补偿发电厂的损失的，在这里，我们引入基本补偿价的概念，利用限上和限下的方法得出阻塞费用的计算规则；
3 仅仅考虑电力市场交易规则和各机组的爬坡速率的影响，要得出下一时段各机组的出力分配预案，就是在各机组的总段容量不超出爬坡速率所限制的极限值的前提条件下，选取价格最低的全部机组的段容量总和，使得这个总和等于负荷需求；
4 用第一问所得出来的有功潮流的表达式检查预案是否会引起输电阻塞，要是会发生输电阻塞，则根据输电阻塞管理原则并考虑电力市场规则的情况下调整预案得出一个新的方案，并根据设计的阻塞费用计算规则求出新方案对于预案的阻塞费用。
模型的假设：

1 各线路上的有功潮流只受到这个电网里面八个发电机组所发的电量的影响，与外界的电网无关，即电网结构不会改变；
2 发电机组出力是按照爬坡速率一直增加或减少到符合出力分配方案的要求后，出力就是一直恒定不变的；
3 发电机出力的增量是由低价格的容量增加到高价格的容量的，而出力的减量就是由高价格的容量减少到低价格的容量；
符号的说明：

： 第i条线路的潮流值（i＝1，2，……6）；
： 序外容量要出力的段容量或其部分；
： 序内容量不能出力的段容量或其部分；
： 第i个机组的出力（i＝1，2，……8）；
： 当前方案0中第i个机组的出力（i＝1，2，……8）；
： 第i个机组在单位时间内的功率增量（i＝1，2，……8）；
： 第i个机组的爬坡速率（i＝1，2，……8）；
： 一个时间段（15分钟）；
： 对应段容量/对应段容量所在机组的爬坡速率；
： 对应段容量的前面容量总和变化到方案0的出力所需的时间；
模型的建立：
问题1
多元线性均值回归模型
根据题中所给统计数据，我们用回归分析方法，考虑线性分析模型，讨论得到：运用全回归分析虽然计算简单，但是最终的计算结果误差比较大。所以我们建立多元线性均值回归模型。
则有：
————（1）

对（1）式的算法如下：
————（2）
————（3）
由表1和表2的1－32方案中知，为稀疏矩阵，元素只有一项i不为零（此项机组表现为相对0方案的变动），且同类的有4组，刚好8个机组的单独变动都有4组数，很明显我们可以均值就可求出i组的k的列向量，运用循环即可全部求得8组k的列向量，合并得k值。（程序附录2和3）

系数矩阵k
＝
观察以上求出的系数k可以发现：k值比较小，而且相邻k值的变化也比较小。
求c时（程序见附录4），利用以上求得k值，有：
————（4）
与对应于1－32方案中的某一方案，由此运用循环当求36组，求均值可得。
常数c＝
所以写出y和p的近似表达式为：
对上述方法进行误差分析，由求出的表达式返回32组方案与原方案对比，得出y的相对误差图（如图一所示）：

如图一可知，此模型的误差较小。

问题2

对于第二问，考虑的是当阻塞发生的时候网方所要补偿给发电商的费用，这个费用称为阻塞费用S，一共包括两个部分：

一部分是序内容量不能出力的部分，决定这部分阻塞费用的主要因素是不能出力的机组停发或少发电的度数以及少发一度电所需的补偿价格；
另一部分是报价高于清算价的序外容量出力的部分，决定这部分阻塞费用的主要因素是序外容量机组多发电的度数以及多发一度电所需的补偿价格。
模型一：限上和限下的方法
根据上面的思想，我们应用限上和限下的方法[1]来计算阻塞费用，规则如下：
1 当出现输电阻塞时，在竞价中未取得发电权的发电容量（序外容量）按报价由低到高的顺序安排机组发电(限上),在竞价中本已经取得发电权的发电容量（序内容量）按报价由高到低的顺序安排机组停发或少发电(限下)。
2 在竞价中未取得发电权的发电容量（序外容量）是在低于其对应报价的清算 价上出力的，这样就会使得发电商利益亏损而网方从中盈利，网方应该对这 些发电商给与经济补偿，这里补偿我们就用对应段容量的报价和清算价之差作为补偿价格；而那些在竞价中本已经取得发电权的发电容量（序内容量）却不能出力的发电商，由于不能出力或部分不能出力，使得这些发电商由于网方线路原因造成经济利益的损失，这时网方也应该给与这些发电商经济补偿，这里我们就以清算价和其不能出力的对应的段容量的报价之差作为补偿价格。
根据计算规则，我们可以列出阻塞费用S的式子：
————（5）
其中 ＝
＝清算价－对应段容量的报价
＝
＝对应段容量的报价－清算价
模型二：
模型二是对于模型一所进行的改进，在模型一中可以发现：对于序内容量，当某个的段容量的报价等于它的清算价的时候，此时＝0，这表明这个段容量是得不到经济补偿的。由于每一个发电厂发的每度电都存在一个利润空间，这个利润是包含在报价里面的，按照模型一的结论表明发电厂按清算价作为报价的段容量的利润得不到经济补偿，显然这是不合理的。
为了解决这个问题，我们在这里引入一个概念：基本补偿价B，根据文献资料[ 2] ，我们定 ————（6），
于是我们就以清算价和其不能出力的对应的段容量的报价之差再加上基本补偿价作为序内容量不能出力的发电商的经济补偿。
于是式子（5）变为： ————（7）

问题3
由题意可知: ————（8）
在下一个时间段的结束时刻各机组的出力大小为
————（9）
根据式子（8）（9），结合表1中方案0各机组的数据和表5中各机组的爬坡速率，相对于方案0，我们可以算出在下一个时间段的结束时刻各机组所允许能选取的出力范围（如下表3－1所示）
表3－1 各机组所允许能选取的出力范围
机组
1
2
3
4
5
6
7
8
最大出力
153
88
228
99.5
152
155
102.1
117
最少出力
87
58
132
60.5
98
95
60.1
63

根据电力市场规则，我们对附录1中的表3，表4进行处理。（按照段价由低到高排列如下表格群所示）表格群3－2
价格
-800
-610
-607
-590
-560
-505
-500
-500
机组
8
3
6
5
1
1
7
4
段容量
70
110
95
75
70
70
50
55

价格
0
0
0
0
0
116
120
124
机组
1
2
3
5
6
5
7
1
段容量
0
0
0
5
0
15
15
50

价格
146
150
152
153
159
168
170
173
机组
5
4
3
8
6
1
4
6
段容量
0
5
40
0
10
0
10
20

价格
180
182
183
188
189
200
203
205
机组
7
2
8
5
3
4
2
6
段容量
5
20
20
15
0
10
8
0

价格
210
215
233
233
245
250
251
252
机组
1
5
3
8
2
5
7
1
段容量
0
15
30
0
15
0
15
30

价格
252
253
255
258
260
283
300
302
机组
6
8
4
3
7
8
2
4
段容量
15
20
10
0
10
0
6
10

价格
303
305
306
308
310
312
315
318
机组
8
6
7
3
5
1
7
8
段容量
20
10
10
20
10
0
5
10

价格
320
325
330
335
348
356
360
363
机组
2
4
1
7
7
3
2
1
段容量
2
15
0
10
3
40
0
0

价格
380
380
396
400
405
410
415
435
机组
6
4
5
8
6
2
3
4
段容量
20
0
10
15
0
0
0
0

价格
489
495
500
510
520
548
800
800
机组
1
2
3
5
6
7
8
4
段容量
40
8
40
10
10
2
5
1

按照以上表格，仅仅考虑电力市场规则，而不考虑各机组所允许的出力范围的时候，选取到第7个表格中的303价格位的段容量就能够满足负荷需求（这个价格位之前的段容量全部选取），这时候我们得到各个机组的出力分别为：
表3－3
机组
1
2
3
4
5
6
7
8
出力
150
79
180
100
125
140
95
113.4

由于受到爬坡速率的限制，第四机组的出力最大值为99.5，不能达到100，将第四组不能达到的0.5段容量按照电力市场规则选取还未选取的且对应段容量价格最低的，
由此，我们得出下一个时段各机组的出力分配预案为（如下表3－4所示）：
表3－4 下一时段各机组的出力分配预案（负荷需求为982.4MW时）
机组
1
2
3
4
5
6
7
8
出力
150
79
180
99.5
125
140
95
113.9

问题4
将上题所得的各机组的出力分配预案代入第一问的模型二中求出的y和p的关系式中，得到各线路的潮流值如下表4－1所示：
表4－1 各线路的潮流值
线路
1
2
3
4
5
6
潮流值
173.3157
96.8095
-150.5051
101.8151
136.8380
168.5319
上表4－1与附录1中的表6相比较可以发现：线路1，5，6会引起输电阻塞。

图 二
根据算法（如图二所示）编写程序（见附录6），由计算机搜索得到：
改变后的出力分配方案如下表4－2所示：

表4－2 各机组的出力分配方案
机组
1
2
3
4
5
6
7
8
出力
153
88
228
90.1775
152
95
60.1
116.1225

此时各线路的潮流值如下表4－3所示：
表4－3 改变方案后各线路的潮流值
线路
1
2
3
4
5
6
潮流值
165
100.6408
-154.8553
105.5524
131.5169
159.5725
并得到对应的阻塞费用为405.8004元。

问题5
利用第三问和第四问同样的方法，得出在下一个时段预报的负荷需求是1052.8MW时，各机组的出力分配预案如下（如下表5－1所示）
表5－1 下一时段各机组的出力分配预案（负荷需求为1052.8MW时）
机组
1
2
3
4
5
6
7
8
出力
150
81
218.2
99.5
135
150
102.1
117

此时各线路的潮流值如下表5－2所示
表5－2 各线路的潮流值
线路
1
2
3
4
5
6
潮流值
177.2608
94.6307
-155.7527
110.0612
134.8521
167.0870

调整之后各机组的出力分配方案为：（如表5－3所示）
表5－3 调整后的各机组出力分配方案
机组
1
2
3
4
5
6
7
8
出力
153
88
228
99.5
152
155
60.3
117
各线路的潮流值为：（如表5－4所示）
表5－4 调整后的各线路的潮流值
线路
1
2
3
4
5
6
潮流值
173.4619
94.1143
-154.9968
104.9631
135.3452
160.4779
从表5－4可以知道，第1，5条线路都是超过限值的，但是都在安全裕度内，第一条线路超出限值的百分比为5.128％，第五条线路超出限值的百分比为2.5342％。
并得出相应的阻塞费用为：94.0377元。
模型的优缺点和改进方向：
优缺点：
本论文中的平均值分析法能在一定范围内很好地吻合数据，经误差分析亦可得相对误差0.0001左右，而且机组系数比较小，能很好地说明此电网有比较好的稳定性，不会因一个机组的小变化而有大的改动，但是，此模型只适应机组出力都较大的情况下，在所有机组都比较小时，此模型有较大误差，但此范围并不在此题的计算范围内，此模型已经能非常很好地解决一般的电网调动。
另外，本文中的计算较多应用matlab编程，对大量的数据计算提供很好的桥梁，针对每一个问题都有一个程序来解决，为以后的利用提供方便。桥梁
改进方向：
由于给出数据中有突变数据，所以在计算得出的相对误差中有6个是大于0.05的，当踢除这些值代以估算的均值，能更完美地完成数据的拟合。
对于赔偿算法中有发电商的隐藏利润，此利润应该赔偿过加入更加合理准确，对此中的专业估计，是改进的重点。

参考文献：
[1] 柯进,管霖 电力市场下的输电阻塞管理技术 电力系统自动化 第26卷 第14期 第四页 2002年
[2] 西部水利水电投资项目介绍
http://159.226.136.229/xbshj/xbbshd/xbbshd_indx.htm 访问时间：2004-9-19
韩忠诚 吕伯雄 金昌市招标项目
访问时间：2004-9-19

附录1：

表1 各机组出力方案 （单位：兆瓦，记作MW）
方案\机组
1
2
3
4
5
6
7
8
0
120
73
180
80
125
125
81.1
90
1
133.02
73
180
80
125
125
81.1
90
2
129.63
73
180
80
125
125
81.1
90
3
158.77
73
180
80
125
125
81.1
90
4
145.32
73
180
80
125
125
81.1
90
5
120
78.596
180
80
125
125
81.1
90
6
120
75.45
180
80
125
125
81.1
90
7
120
90.487
180
80
125
125
81.1
90
8
120
83.848
180
80
125
125
81.1
90
9
120
73
231.39
80
125
125
81.1
90
10
120
73
198.48
80
125
125
81.1
90
11
120
73
212.64
80
125
125
81.1
90
12
120
73
190.55
80
125
125
81.1
90
13
120
73
180
75.857
125
125
81.1
90
14
120
73
180
65.958
125
125
81.1
90
15
120
73
180
87.258
125
125
81.1
90
16
120
73
180
97.824
125
125
81.1
90
17
120
73
180
80
150.71
125
81.1
90
18
120
73
180
80
141.58
125
81.1
90
19
120
73
180
80
132.37
125
81.1
90
20
120
73
180
80
156.93
125
81.1
90
21
120
73
180
80
125
138.88
81.1
90
22
120
73
180
80
125
131.21
81.1
90
23
120
73
180
80
125
141.71
81.1
90
24
120
73
180
80
125
149.29
81.1
90
25
120
73
180
80
125
125
60.582
90
26
120
73
180
80
125
125
70.962
90
27
120
73
180
80
125
125
64.854
90
28
120
73
180
80
125
125
75.529
90
29
120
73
180
80
125
125
81.1
104.84
30
120
73
180
80
125
125
81.1
111.22
31
120
73
180
80
125
125
81.1
98.092
32
120
73
180
80
125
125
81.1
120.44

表2 各线路的潮流值（各方案与表1相对应，单位：MW）
方案\线路
1
2
3
4
5
6
0
164.78
140.87
-144.25
119.09
135.44
157.69
1
165.81
140.13
-145.14
118.63
135.37
160.76
2
165.51
140.25
-144.92
118.7
135.33
159.98
3
167.93
138.71
-146.91
117.72
135.41
166.81
4
166.79
139.45
-145.92
118.13
135.41
163.64
5
164.94
141.5
-143.84
118.43
136.72
157.22
6
164.8
141.13
-144.07
118.82
136.02
157.5
7
165.59
143.03
-143.16
117.24
139.66
156.59
8
165.21
142.28
-143.49
117.96
137.98
156.96
9
167.43
140.82
-152.26
129.58
132.04
153.6
10
165.71
140.82
-147.08
122.85
134.21
156.23
11
166.45
140.82
-149.33
125.75
133.28
155.09
12
165.23
140.85
-145.82
121.16
134.75
156.77
13
164.23
140.73
-144.18
119.12
135.57
157.2
14
163.04
140.34
-144.03
119.31
135.97
156.31
15
165.54
141.1
-144.32
118.84
135.06
158.26
16
166.88
141.4
-144.34
118.67
134.67
159.28
17
164.07
143.03
-140.97
118.75
133.75
158.83
18
164.27
142.29
-142.15
118.85
134.27
158.37
19
164.57
141.44
-143.3
119
134.88
158.01
20
163.89
143.61
-140.25
118.64
133.28
159.12
21
166.35
139.29
-144.2
119.1
136.33
157.59
22
165.54
140.14
-144.19
119.09
135.81
157.67
23
166.75
138.95
-144.17
119.15
136.55
157.59
24
167.69
138.07
-144.14
119.19
137.11
157.65
25
162.21
141.21
-144.13
116.03
135.5
154.26
26
163.54
141
-144.16
117.56
135.44
155.93
27
162.7
141.14
-144.21
116.74
135.4
154.88
28
164.06
140.94
-144.18
118.24
135.4
156.68
29
164.66
142.27
-147.2
120.21
135.28
157.65
30
164.7
142.94
-148.45
120.68
135.16
157.63
31
164.67
141.56
-145.88
119.68
135.29
157.61
32
164.69
143.84
-150.34
121.34
135.12
157.64

表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
速率
2.2
1
3.2
1.3
1.8
2
1.4
1.8

表6 各线路的潮流限值（单位：MW）和相对安全裕度
线路
1
2
3
4
5
6
限值
165
150
160
155
132
162
安全裕度
13%
18%
9%
11%
15%
14%
附录2
function k=zhen(fa8,r6) %求k初值
p0=[164.78 140.87 -144.25 119.09 135.44 157.69];
pj=[120 73 180 80 125 125 81.1 90];
j=0;
for i=1:8
j=fa8(i)-pj(i)+j;
end
pi=r6-p0;
k=pi./j;
附录3
function ok=zhen1(fa8,r6) %求k系数
clc
for i=1:32
a=fa8(i,:);b=r6(i,:);
k(i,:)=zhen(a,b);
end
options=[];
for i=1:8
for j=1:6
ok(i,j)=mean([k(4*i,j),k(4*i-1,j),k(4*i-2,j),k(4*i-3,j)]);
end

end

附录4
function c=zhen4(ok,fa8,r6) %求c值

for k=1:32
if k==11
c(k,:)=c(k-1,:);
end
for j=1:6
c(k,j)=r6(k,j)-fa8(k,:)*ok(:,j);
end
end

c=mean(c);

附录5
function pout=zhen5(pin) %求潮流值
ok=[0.08284 0.04828 0.05297 0.11993 -0.02544 0.12201 0.12158 -0.00123;...
-0.05456 0.12785 -0.00003 0.03328 -0.08685 -0.11244 -0.01893 -0.09873;...
-0.06954 0.06105 -0.15662 -0.00992 0.12449 0.00212 0.00251 -0.20139;...
-0.03446 -0.10241 0.20516 -0.02083 -0.01183 -0.00595 0.14492 -0.07655;...
0.00053 0.24329 -0.06455 -0.04113 -0.06522 0.07034 -0.00426 -0.00891;...
0.23781 -0.06017 -0.07787 0.09298 0.04690 0.00008 0.16593 0.00069];
c=[110.2965 131.2289 -108.873 77.4817 132.97447 120.6633];
pout=ok*pin'+c';

附录6
function f=zhenfun2(x)
x0=[120 73 180 80 125 125 81.1 90];
x1=[150 81 218.2 9