任务：开展三峡水库中长期调度研究
任务要求(4.14)
进行水库长系列常规调度模拟计算（常规，非优化）
- 一般采用常规调度图作为调度依据。
进行水库长系列优化调度模拟计算
- 约束：末水位为常规调度末水位
- 目标函数：
- 算法：
需要通过惩罚系数调整发电保证率，得出惩罚系数、发电保证率和发电量三者的关系。
在指定年份运行常规调度和优化调度的模型，获得结果并对比各特征值
个性化数据
- 计算初期水位：173.45m
- 优化调度最小下泄流量：5690m³/s
- 计划年：1975
明日任务：完成任务一思路构建
- idea：根据95%的设计保证率进行年平均流量排频处理，输入P-Ⅲ型曲线代码得到相近年份，对照调度图进行常规调度模拟。
使用方法（参照）(4.15)
- 经验适线法
- 点绘样本经验点据
将实测资料由大到小排列，计算各年流量经验频率，然后将经验点据绘制在机率格纸上。
- 估计参数的初值并绘制频率曲线
采用矩法或其他方法，估计分布的3个参数，记为（\(\overline{x}\), \(C_{v}\) , \(C_{s}\) ),作为适线法的初值。根据该参数查P-Ⅲ型分布\(\Phi\)值表，可以求得一组不同频率p对应的设计值\(x_{p}\)，即
\[ x_{p}=\overline{x}*\left[1+C_{p}*\Phi(p,C_{s})\right] \]
根据(p,\(x_{p}\))绘制频率曲线（也称为理论频率曲线),并将此线画在绘有经验点据的机率格纸上。
- 调整适线
检查频率曲线与经验点据的拟合情况，若不理想，则调整参数（主要是调整\(C_{s}\),\(C_{v}\))，再重新计算频率曲线。
point
在计算时需要根据已有的\(\Phi\)值表读出不同\(C_{v}\),\(C_{s}\)下的\(K_{p}\)值进行择优选取曲线，当比值大于3时，则需要慎重考虑取值，可以反算前面的数据是否出现问题。
选年
经过排频分析，得出95%设计保证率选年可以选为1959年。
任务一存在问题：
已知发电量，如何确定水头
由出力公式，\(N=kQH\)，由水力学关系，有k(H),h(Q),\(H_{下}\)(Q),故所有参数均为下泄流量Q的函数。
如何使已知的一切数据表通过代码联动起来，实现导入运行即可导出
- 为试算过程。
需要解决程序如何在各阶段都能够线性插值，表格的使用顺序是什么样的，如何放置各阶段的约束条件。
明日任务
完成程序的基本调试。
单时段常规调度程序(4.16)
目前已经完成常规调度基本程序，需要进一步导入整个表格实现计算。以下为单个时段计算程序。
思路存在问题：试算的结果无法进行及时的更新，无法确定具体的值。故不应将N设为试算结果评判标准，或应设置适当的Q的自增值。
- 以下为上述思路代码
代码一
1 | {%%程序一 |
最终结果误差小于0.0001，故接受。
尚未解决的问题
- 如果初始流量值大于实际流量值，如何收敛到正确结果；
- 初始填入的数据仍然需要通过读图读表得到，如何将其加入代码实现自动读图；
- 在得出实际下泄流量后，还需要计算时段末水位和库容。库容还需设置边界条件(这里是很重要的数据边界条件)。
明日任务
完善程序细节，尽量解决上述问题。
任务一终稿（有问题）(4.17)
最终联动程序设计思路： - 根据已得出的时段初水位读出当前时段发电量； - 根据入库流量假设出库流量，列入程序一进行试算； - 利用试算结果得出时段末水位，得出下一时段计算初始条件，重复计算。
需要逐步解决的问题：
1.如何设置数组条件；
2.如何在自定义函数中使用工作区变量；
3.三峡汛限水位和死水位为同一高度，即在6月-9月之间无规定发电量，来多少泄多少发多少，如何设置程序具体功能；
4.如何限制最大发电量，最大下泄流量，以及计算弃水
目前已解决1、2、3问题
遗留问题：如何确定下泄流量初值；如何限制最大发电量，最大下泄流量，以及计算弃水。
代码二
以下为主程序代码：
1 | {%changguidiaodu.m |
读取调度图实际出力功能代码：
1 | {%NFun.m |
以下为已知出力试算流量代码：
1 | {%%计算发电量函数 |
存在问题
- 全功能代码过于复杂，并不容易调试，现尝试分函数块逐情况处理。
明日任务
优化处理任务一程序，完善功能。
课后新思路(4.19)
- 155m为非汛期最低水位，即真正的死水位，故在枯水期调度时水位不能低于155；
- 可以考虑对每一时段进行以电定水试算（查阅水文水利计算课本）
存在问题：
- 两个主要功能函数存在循环嵌套的问题，需要细改
以下为新代码：
代码三
1 | {%changguidiaodu |
现在出现的问题关键在于想实现全自动读取数据，也许有贪心不足的情况，接下来的任务，还是继续跑第一次的程序，手动算出一组数据再进行改进。目前思路出现了一些循环论证的问题，需要搞清楚怎么把所需要试算的数据拆开再循环，应该是逻辑本身有问题，下次推进的时候先解决一下整个计算的试算步骤。
明日任务
解决函数逻辑问题。
今日无具体推进。(4.20)
上课间隙简单跑了一下昨天最后的思路，发现一开始就出现了问题，因为水头本身很低，所以最后循环出来的时段末水位甚至成为了负值，这就导致以电定水的循环无法结束，故回到最初的思路，即在计算完出力后在循环内立即计算时段末水位，若不低于145，则继续迭代流量值，若低于145，则转为来多少发多少，计算结束后直接退出循环。
这里也是上文所述的出现了循环嵌套的问题，今天所产生的新思路为，将正常插值计算时段末水位和来多少发多少计算时段末水位分开成两个函数进行使用，在出现水位过低的情况时，调用后者，并在计算完成后直接退出循环。
仍需实践探索。
明日任务
争取改善时段循环的关键代码。
正式完成水电站常规调度。最终代码如下：(4.22)
代码四
主程序代码：
1 | {changguidiaodu.m |
计算当前水位设计出力代码：
1 | {NFun.m |
计算实际发电出力
1 | {Nqfun.m |
计算正常情况下时段末水位： 1
2
3
4
5
6
7
8
9
10{ZVendnormal.m
function w = ZVendnormal(x)%x为时段末库容
for i = 1:55
if x >= ZV(i,2) && x <= ZV(i+1,2)
w = ((ZV(i+1,1)-ZV(i,1))/(ZV(i+1,2)-ZV(i,2)))*(x-ZV(i,2))+ZV(i,1);
elseif x > ZV(56,2)
w = ((ZV(56,1)-ZV(55,1))/(ZV(56,2)-ZV(55,2)))*(x-ZV(55,2))+ZV(55,1);
end
end
}
计算发电水头： 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46{head.m
%%水头计算函数
function [h51,kcl] = head(x,y)%下泄流量,时段初水位
%插值下游水位流量曲线
for i = 1:29
if x < ZO(1,1)
hxy = ((ZO(1,2) - ZO(2,2))/(ZO(1,1) - ZO(2,1)))*(ZO(1,1)-x)+ZO(2,2);
break;
elseif x >= ZO(i,1) && x <= ZO(i+1,1)
hxy = ((ZO(i+1,2) - ZO(i,2))/(ZO(i+1,1) - ZO(i,1)))*(x-ZO(i,1))+ZO(i,2);
break;
end
end
%%插值发电水头损失曲线
for i = 1:13
if x < Odh(1,1)
hss = ((Odh(1,2) - Odh(2,2))/(Odh(1,1) - Odh(2,1)))*(Odh(1,1)-x)+Odh(2,2);
break;
elseif x >= Odh(i,1) && x <= Odh(i+1,1)
hss = ((Odh(i+1,2) - Odh(i,2))/(Odh(i+1,1) - Odh(i,1)))*(x-Odh(i,1))+Odh(i,2);
break;
end
if x >= Odh(15,1)
hss = ((Odh(15,2) - Odh(14,2))/(Odh(15,1) - Odh(14,1)))*(x-Odh(14,1))+Odh(14,2);
break;
end
end
%%插值出力系数水头曲线
h51 = y - hss - hxy;
for i = 1:5
if h51 < K(1,1)
kcl = ((K(1,2) - K(2,2))/(K(1,1) - K(2,1)))*(K(1,1)-h51)+K(2,2);
break;
elseif h51 >= K(i,1) && h51 <= K(i+1,1)
kcl = ((K(i+1,2) - K(i,2))/(K(i+1,1) - K(i,1)))*(h51-K(i,1))+K(i,2);
break;
end
if h51 >= K(6,1)
kcl = ((K(6,2) - K(5,2))/(K(6,1) - K(5,1)))*(h51-K(5,1))+K(5,2);
break;
end
end
}
最终得到枯水年枯水期发电水位变化曲线为：
如图，在全年除汛期应正常维持在死水位外，其他时期的水头对发电来说过低，整体发电量很少。所选择的枯水年来水极少是一方面的原因，还有一方面，即直接使用调度图可能会导致竭泽而渔，彻底不给后面的来水情况做预备。
接下来的任务是在已知初水位末水位的情况下，在汛限水位的约束下对发电和水位曲线进行优化调度。即任务二：
进行水库长系列优化调度模拟计算
- 约束：末水位为常规调度末水位
- 目标函数：
- 算法：
需要通过惩罚系数调整发电保证率，得出惩罚系数、发电保证率和发电量三者的关系。
明日任务
初步寻找所使用的优化调度方法和方程建立思路。
(4.23) 上午上课发现之前计算的初始水位选的有点太低，所以回来又改了一下常规调度的代码，矫正了一些参数的输出，变动不大，就不放出来了。
晚上在查找优化调度使用的方法，在多约束又多阶段的情况下，决定使用微分动态规划进行优化计算。需要具体学习一下微分动态规划的理论和建立方程的方法，作为明天的任务。
查到有关微分动态规划理论的网址如下：
(4.24) 昨晚观察表格发现后三个月来水并不少，所以决定今天首先把时段末水位重新蓄回正常蓄水位，水位蓄高之后发电量与之前相比也大了很多。修改的主程序代码如下，其余功能函数没有发生改变：
代码五
1 | { |
明日任务
首先尝试使用最基本的动态规划解决问题。先建立好动态规划方程并设置计算的数据结构。
（5.6）
首先建立逆序递推方程
1.阶段与阶段变量
根据中长期规划，按旬划分为36个阶段，以t表示阶段变量，则t为面临时段，t+1~T为余留时期；
2.状态变量
选择水库蓄水量为状态变量，记\(V(t)\)为t时段初蓄水量；
3.决策变量
取水电站发电引用流量\(Qfd(t)\)为决策变量；
4.状态转移方程
即水量平衡方程，\(V(t+1) = V(t) + [Q_{rk}(t) - Q_{fd}(t) - Q_{qs}(t)]*\Delta T\);
5.指标函数
1 | { |
6.递推方程
\[\left\{ \begin{array} Pbest(t+1) = max(Pcl(t)+Pbest(t)); \\ V(t+1) = V(t) + (Qrk(t)-Qfd(t)-Qqs(t))*T; \end{array} \right.\]
7.约束条件
- 上下游水头差计算（函数约束）
- 水电站预想出力限制（范围约束）
- 水电站最大过流能力限制（范围约束）
- 综合利用约束（范围约束）
- 水库水量平衡（函数约束）
- 库容曲线约束（函数约束）
- 库水位或库蓄水量约束（范围约束）
- 下游水位流量关系约束（值约束）
即得如下代码：
代码六
1 | { |
预想出力限制计算函数：
1 |
|
库容曲线约束计算函数：
1 | { |
写出计算方程后，出现如下问题
- 问题1：只能计算下一个时段，如何实现全时段计算？如何将一个时段内的计算初始条件划分为多个段进行计算？如何设置数组？如何设置循环次数？
- 问题2：约束条件很多，如何写入程序？写在哪里？如何实现变量约束条件内的改变？
代码七
1 | { |
以上代码基本解决上述问题，出现新问题：
- 为什么总是跳出循环？如何解决？
- 为什么计算结果里依然出现了不符合约束条件的值？
- 如何确保最后一个水位为145？[增加if语句，若j=15则进入，单独计算流量值]
明日任务
解决上述问题。
（5.7） ## 问题：在哪里以及如何改变Q？
解决：Qfd为事先定义的数组。现在将约束条件变为if语句，加入Qfd初值设置，即将Q设为决策变量，相比于课堂所讲方法，他将水位作为决策变量。
问题：出现维数灾，不知道需要计算多久。
代码八
1 | { |
明日任务：等待结果。
（5.8）
发现问题：Qfd的改变过于死板，导致程序会陷入一个时段的死循环
完成优化调度，代码如下：
代码九
1 | { |
得到水位变化如下图：
使用1975年数据进行演算，得到优化调度保证率为77%，全年发电量13562.95wkW,常规调度保证率为45%，全年发电处理7916.28wkW,有显著的优化效果，可见调度图调度的保守性。
（5.9） 需要在程序中加入弃水计算部分，经过查找资料，得到弃水的出现条件是在水头一定的条件下电站出力不得大于水电站的预想出力，若大于，则多余水量为弃水.代码如下：
代码十
1 | { |
代码十一
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
%%水库优化调度代码
Qst = 5690;%生态流量限制
T = 12*24*3600/100000000;%时段间隔
Pb = 499;%保证出力
V(1,1) = 393;
flg = 0;%标志位
m = 0;%计算保证率指标
Pbest = zeros(15000,15);
for c = 1:46
b = punish(c,1);%设置惩罚系数
%优化计算过程
for j = 1:15
Z1(1,1) = 175;
%下泄流量范围
l = length(Qst:1:Qst+25000);
Qfd(1:l,j) = (Qst:1:Qst+25000)';
flg = flg+1;%循环标志位
for i = 1:15000
if V(1,j) < 171.5
if i-1 == 0
f = 1;
Qfd(i,j) = Qfd(i-1,j);
break;
end
if f == 1
faulce(1,j) = 1;
break;
end
%%约束条件编制
%发电水头约束
[H(i,1) K] = head(Qfd(i,j),Z1(1,j));%输入下泄流量及上游水头
%水电站预想出力限制
yxcl(i,1) = Yuxcl(H(i,1));%输入当前水头
%水电站最大过流能力限制
xxnl(i,1) = Xiaxnl(Z1(1,j));%输入当前上游水位
%指标函数（含惩罚函数）
output(i,j) = K * Qfd(i,j) * H(i,1)/10000;
if output(i,j) > yxcl(i,1)
continue;
end
if output(i,j) >= Pb
Pcl1(i,j) = output(i,j);
else
Pcl1(i,j) = output(i,j) + b*(output(i,j)-Pb);
end
%递推方程
Pbest(1,j+1) = max(Pcl1(i,:)+Pbest(1,j));
end
%建立状态转移方程
V(1,j+1) = V(1,j) + (Qrk1(j,1)-Qfd(1,j))*T;
%库容曲线约束
Z1(1,j+1) = Kurqx(V(1,j+1));%当前时段末水位
Z1(1,16) = 145;
end
%统计发电出力
n = size(Pcl1,1);
for i = 1:15
for j = 1:n
if Pcl1(j,i) == 0
P1(1,i) = Pcl1(j-1,i);
break;
end
if Pcl1(j,i) <= 0
break;
end
if Pcl1(j,i) >= 0 && j == n
P1(1,i) = Pcl1(j,i);
end
end
end
%
%优化计算过程2
V(1,1) = 171.5;
for j = 1:9
Z2(1,1) = 164.27;
%下泄流量范围
l = length(Qst:1:Qst+15000);
Qfd(1:l,j) = (Qst:1:Qst+15000)';
flg = flg+1;%循环标志位
for i = 1:15000
if V(1,j) < 171.5
Qfd(i,j) = Qfd(i-1,j);
break;
end
%%约束条件编制
%发电水头约束
[H(i,1) K] = head(Qfd(i,j),Z2(1,j));%输入下泄流量及上游水头
%水电站预想出力限制
yxcl(i,1) = Yuxcl(H(i,1));%输入当前水头
%水电站最大过流能力限制
xxnl(i,1) = Xiaxnl(Z2(1,j));%输入当前上游水位
%指标函数（含惩罚函数）
output(i,j) = K * Qfd(i,j) * H(i,1)/10000;
if output(i,j) > yxcl(i,1)
continue;
end
if output(i,j) >= Pb
Pcl2(i,j) = output(i,j);
else
Pcl2(i,j) = output(i,j) + b*(output(i,j)-Pb);
end
%递推方程
Pbest(i+1,j) = max(Pcl2(i,j)+Pbest(i,j));
end
%建立状态转移方程
V(1,j+1) = V(1,j) + (Qrk2(j,1)-Qfd(1,j))*T;
%库容曲线约束
Z2(1,j+1) = Kurqx(V(1,j+1));%当前时段末水位
Z2(1,16) = 175;
end
%统计发电出力
n = size(Pcl2,1);
for i = 1:9
for j = 1:n
if Pcl2(j,i) == 0
P2(1,i) = Pcl2(j-1,i);
break;
end
if Pcl1(j,i) <= 0
break;
end
if Pcl2(j,i) >= 0 && j == n
P2(1,i) = Pcl2(j,i);
end
end
end
for i=1:15
endZ(i,1) = Z1(1,i);
end
for i = 16:27
endZ(i,1) = 145;
end
for i = 28:36
endZ(i,1) = Z2(1,i-27);
end
plot(1:36,endZ(1:36,1));
xlim([1 36])
ylim([145 185])
ylabel("水位/m");
xlabel("旬");
title("1975年来流优化调度水位变化曲线");
%}
%%计算发电保证率
for k = 1:15
P(1,k) = P1(1,k);
end
for k = 16:24
P(1,k) = P2(1,k-15);
end
for d = 1:24
if P(1,d) < 499
m = m+1;
end
end
percent(c,1) = 1-m/24;
end
end
1 | %%水库优化调度代码 |
基本完成课程设计任务。从最终得到的数据曲线来看，动态规划的问题比较大，另外两种方法均无法达到保证出力，是还没有解决的问题。常规调度第二段若完全按照调度图调度，整体发电甚至会更少。如果仍有时间，可以考虑改变出力目标，尝试全部改为保证出力计算。动态规划需要先做一些简单的小题再理解一下决策变量如何改变。
先这样。over