本文是参加工作之初，研究水电站优化运行时所写，已成为遥远的过去了。记得当时好像是要投下稿，所以煞有介事的写了中英文摘要。学生时代曾经研究过用演化类算法求解负荷分配问题，那是刻意往“新”上靠，因为动态规划方法是写进教科书的传统方法，不宜做毕业论文。实际上，只要算力可及，传统动态规划方法比引入随机因素的演化算法要“靠谱”，演化算法的好处在于“简单粗暴”，按它的套路总能够处理，维数灾的“灾点”比较高。

一个问题，当自己不再关注时，热门也成了冷门。这个问题本就比较小众，现在或许没什么人研究它了吧。压箱底好多年，贴出来，供“万一需要的人”参考。

负荷分配问题的动态规划算法递归实现

摘要

Abstract

1 引言

2 算法设计

2.1 基本思想

2.2 离散化与可行域

2.3 递归过程

3 数值试验

4 结语

参考文献

源程序：（基于 Delphi 实现）

负荷分配问题的动态规划算法递归实现

摘要

本文研究了水电站优化运行中的负荷分配问题，给出了可行解的充分必要条件暨构造方法，在此基础上设计了动态规划算法的一个递归方案，并通过数值试验，验证了算法设计的有效性。该递归设计思路清晰、简单易行，可以在工程实际中参考实施。

关键词：水电站；优化运行；负荷分配；动态规划；递归

Abstract

This paper presents a study on the power dispatch problem of hydropower station economic operation. A rule to construct feasible solutions is proposed. Then a recursion approach of dynamic programming (DP) is designed based on the rule. After all, the efficiency of this new method is showed through a numerical example. The recursion design of DP on the power dispatch problem is clear, easy to implement and deserves referring to in real project.

Key words: hydropower station; economic operation; dynamical programming; recursion

1 引言

水电站厂内优化运行作为整个水电站优化运行的一部分，主要研究机组特性曲线，负荷的优化分配等，以在有限资源条件下获得尽可能多的经济效益。本文考察如下基本的负荷分配问题：

设有n台机组，总的发电负荷任务为N，已知各机组耗流量特性曲线 $Q_{i}(N_{i})$ 以及出力范围，要求将总的发电任务在n台机组间进行分配，使得总的耗水量最小。表达为数学模型即

其中，N满足可行性条件。

实际中求解上述问题的常用方法有等微增率法、动态规划法、演化算法等，动态规划法是其中较突出的一种。[1]中介绍了动态规划求解问题 (1) 的思想方法，但其流程复杂，可读性较差，在研究动态规划求解负荷分配问题的过程中，我们找到了一种简单易行的递归实现方法。下文分别对该方案的基本思想、技术要点、递归过程进行详细介绍。

2 算法设计

2.1 基本思想

动态规划适用于求解多阶段决策过程，使用动态规划方法求解负荷分配问题 (1)，需把原问题转化为多阶段决策过程。为此，定义k阶段状态变量

其中，满足可行性条件。则由动态规划基本原理，对，有如下递推关系式

成立。从而，问题 (1) 转化为通过 (2)-(3) 式求。为求某k阶段状态变量，若对每一个可行的，相应的k-1阶段变量已求出，则不难得到。

2.2 离散化与可行域

使用动态规划求解负荷分配问题 (1)，需要作变量的离散化处理，如果直接在变量的给定区间上离散化，将包含进对许多不可行解的计算尝试，浪费大量计算资源。关于问题 (1) 的可行域，我们给出如下定理。

定理1（可行解的充分必要条件1）是问题(1)的一个可行解的充分必要条件是满足下式

证明：必要性。若是可行解，即满足 (1) 的约束条件，亦即满足

从而

进而易见 (4) 式成立。

充分性。若满足 (4) 式，则只需证明，即可说明是可行解。由 (4) 式知

令，得

即有成立。

证毕。

定理2（可行解的充分必要条件2）是问题(1)的一个可行解的充分必要条件是满足下式

注：将定理1中作“倒序置换”，即得本定理形式。定理2可视作定理1的“对偶定理”。

根据上述定理，可以逐次确定变量的可行区间，比如负荷S在1~k号机组间分配，首先第k号机组的可行区间可以确定，一旦k机组分得一个负荷值，第k-1号机组的可行区间也可以确定，如此下去即可保证始终在可行域内搜索。后面设计的递推过程，将以定理2为理论基础，在整个可行域上对变量动态地进行逐次离散化。

2.3 递归过程

由动态规划的递推公式及上文关于可行域的讨论，设计递归过程StateRecursion：StateRecursion(

) 确定在1~k号机组间分配 $S_{k}$ 的负荷时的最小耗流量，并记录k号机组此时分得的负荷。具体过程如下：

算法 StateRecursion()

1) 若k >1，即在多台机组间分配

1.1) 根据定理2对在可行域内以一定步长进行离散化；

1.2) 对每一个可行的，递归计算k-1阶段状态变量：

1.3) 计算K阶段状态变量（最小耗流量）

1.4) 保存K号机组的最优负荷，即保存所有取得最小耗流量的值；

2) 若 k = 1，即分配给一台机组

2.1) 计算K阶段状态变量（最小耗流量）

2.2) 保存K号机组的最优负荷

3) 返回；

3 数值试验

根据文献[1，P133]表4-1，即下表1机组流量特性曲线数据，用前面介绍的递归方案进行计算，得到“优化分配总表”如表2所示。

表 1 机组流量特性曲线

N Q(k, N) k	1	2	3	4
0	25	24	25	23
1	30	31	29	28
2	34	34	35	35
3	38	37	38	37
4		40

表1注：k为机组编号，N为机组出力，Q(k, N)为k号机组出力为N时的耗流量。例如1号机组出力为3时的耗流量为38，2号机组出力为2时耗流量为34。

表 2 优化运行总表

S N(k, S) k	1	2	3	4	Q(S)
0	0	0	0	0	97
1	0	0	1	0	101
2	2	0	0	0	106
2	1	0	1	0	106
2	0	0	1	1	106
3	3	0	0	0	110
3	0	3	0	0	110
3	2	0	1	0	110
3	0	0	3	0	110
4	0	4	0	0	113
5	0	4	1	0	117
6	2	4	0	0	122
6	1	4	1	0	122
6	0	4	1	1	122
7	3	4	0	0	126
7	2	4	1	0	126
7	0	4	3	0	126
8	3	4	1	0	130
9	2	4	3	0	135
9	3	4	1	1	135
10	3	4	3	0	139
11	3	4	3	1	144
11	3	4	1	3	144
12	2	4	3	3	149
13	3	4	3	3	153

表2注：k为机组编号（k=1~n，这里n=4），S为总待分配负荷，N(k,S)为总负荷S在n台机组间分配时k号机组分得的负荷。例如，总负荷为12时，计算得最优负荷分配方案为1~4号机组的负荷依次为2、4、3、3，总耗水量为149。

将表2结果与[1, P134]中相应的表4-3对比，并通过验算发现，[1]中总负荷为2、7、8、9、10、11时结果或有遗漏或有错误，而其余结果是一致的。关于[1]中的错误，可能是编辑与印刷带来，也有可能是根据[1]中表4-2计算整理不当导致。此例初步表明，本方案的设计是有效的。

4 结语

本文所关注的问题 (1) 考虑的是负荷在固定机组间的优化分配，这是水电站优化运行的一个基本问题。求解此问题，计算得优化运行总表后，任意机组间的最优分配可以通过查表计算获得，若要考虑机组的气蚀振动约束，在程序中增加是否在气蚀振动区的判断或者通过惩罚躲避气蚀振动区即可实现；进一步，在优化运行总表基础上，还可以根据日负荷特性制定出日开停机计划[1]。基于问题的重要性，本文介绍了动态规划求解问题(1)的一个递归方案，首先给出对可行域进行离散化的技术，然后设计了具体的递归过程，最后选用[1]中的例子对我们的方案进行了测试。本动态规划递归实现方案，逻辑清晰、形式简洁，易于实现和维护，可在工程实际中参考实施。

参考文献

[1] 张勇传. 水电站经济运行原理（第二版）[M]. 北京：中国水利水电出版社，1998.

源程序：（基于 Delphi 实现）

(*动态规划算法实现================2011-12-9,hll
适用问题形式：min sum_{i=1}^{n}f_i(x_i)s.t. lb_i <= x_i <= ub_ix_1 + x_2 + ... + x_n = Sum
应用例子：水电站厂内优化运行负荷分配问题；
使用方法：
*)

unit DP;

interface


typeTFunc =  function(k: Integer; x: Double):Double of object;  // 特性曲线函数；TProbPara = record              // 问题参数结构;n: Integer;                   // 问题中的n值，变量x_i, i=1..n;ub: array of Double;          // 变量上界（Upper Bounds）, Ub[0..n];lb: array of Double;          // 变量下界（Lower Bounds），Lb[0..n];dx: Double;                   // 变量的分辨率（Delta X），即离散化步长；f: TFunc;                     // 特性曲线；end;

TAuxRec = recordS: Double;X: array of array of Double;U: array of Double;end;TAuxTable = array of TAuxRec;  // DP 计算辅助表;
TResRec = recordS: Double;X: array of array of Double;U: Double;end;TResTable = array of TResRec;  // DP结果表（优化运行总表）；
TDP = class(TObject)privateFProbPara: TProbPara;          // 问题信息；FUbSum: Double;                // 变量上界和；FLbSum: Double;                // 变量下界和；FM: Integer;                   // 离散化段数（有效待分配负荷的个数）；FAuxTable: TAuxTable;  // 动态规划计算辅助表；FResTable: TResTable;  // 动态规划结果表（在优化运行中称优化运行总表）；protectedfunction GetSum(AnArray: array of Double; LowPos, HighPos: Integer): Double;function StateRecursion(K: Integer; Sum: Double): Double; // 递归计算；publicconstructor Create(AProbPara: TProbPara);function CompAuxTable: Integer;  // 计算得辅助表；function CompResTable: Integer;  // 生成结果表；function CompResRec(K,M: Integer;var ResRec: TResRec): Integer; overload;// 生成一个结果记录，将第M(0..FM)个离散值,分配给1~K个自变量的最优解；function CompResRec(K: Integer; Sum: Double;var ResRec: TResRec): Integer; overload;// 生成一个结果记录，将Sum分配给1~K个子变量的最优解；publishedproperty ResTable: TResTable read FResTable;end;

implementation

usesMath, Dialogs, SysUtils;

{ 计算辅助表 FAuxTable；}
function TDP.CompAuxTable: Integer;
varM:Integer;
beginfor M:= 0 to FM doStateRecursion(FProbPara.n, FLbSum + M * FProbPara.dx);   // 2012-3-27, 修正；Result :=0;
end;

{ 计算结果表 FResTable；返回值：0正常，1不正常；}
function TDP.CompResTable: Integer;
varM: Integer;
beginResult := 0;SetLength(FResTable, FM + 1);for M:=0 to FM doif CompResRec(FProbPara.n, M, FResTable[M]) <> 0 thenbeginResult := 1;Break;end;
end;

{ 计算一组结果记录；参数：K: 变量个数，1~K；Sum: 对应总和FAuxTable[?].S；ResRec：将对应于M的和,在1~K个变量之间分配的最优结果记录；返回值：0正常，1不正常；
}
function TDP.CompResRec(K: Integer; Sum: Double; var ResRec: TResRec): Integer;
varM: Integer;
beginM :=  Round((Sum-FLbSum) / FProbPara.dx);Result := CompResRec(K, M, ResRec);
end;

{ 计算一组结果记录；参数：K: 变量个数，1~K；M: 对应 FAuxTable[M].S 的总和；ResRec：将对应于M的和,在1~K个变量之间分配的最优结果记录；返回值：0正常，1不正常；
}
function TDP.CompResRec(K, M: Integer; var ResRec: TResRec): Integer;
vari, j, MM, KK: Integer;Trace: Array of array of Integer;      // 解析AuxTable的辅助用表；
{ 辅助过程：横向补全Trace表；}
procedure HTraceComplete(aK: Integer);
vari: Integer;
beginif (aK >= 2) and (aK <= K) thenfor i := aK downto 2 dobeginMM := Trace[0,i];j := Trace[1,i];Trace[0, i-1] := MM - Floor(FAuxTable[MM].X[i,j] / FProbPara.dx);Trace[1, i-1] := 1; // 初始位置为1，范围为1~FAuxTable[MM].X[aK,0]；end;
end;
{ 辅助过程：纵向计算Trace表下一个格局；}
function VTraceNext: Integer;
vari: Integer;
beginResult := 0;  // Trace表格局更新；kk := 0;for i:= 1 to K dobeginMM := Trace[0,i];j := Trace[1,i];if j >= FAuxTable[MM].X[i,0] thenContinueelsebeginInc(Trace[1,i]);kk := i; // 更新的节点位置；Break;end;end;if kk > 1 thenHTraceComplete(kk)else if kk = 0 thenResult := 1; // 已经是最后格局，不能再更新；
end;
{辅助过程：Trace表格局转换为结果记录；}
procedure TraceToResRec(aRow: Integer);
varKK: Integer;
beginfor KK := 1 to K dobeginMM := Trace[0, KK];j := Trace[1, KK];ResRec.X[aRow, KK] := FAuxTable[MM].X[KK, j];end;
end;
beginResRec.S := FAuxTable[M].S;ResRec.U := FAuxTable[M].U[K];{ Trace[0,k] = m: 记录ResRec.X[*,k]在FAuxTable中的行位置m；Trace[1,k] = j: 记录ResRec.X[*,k]取FAuxTable[m].X[k,j]；}SetLength(Trace, 2, K+1);// Trace 表初始格局处理；Trace[0,K] := M;Trace[1,K] := 1;HTraceComplete(K);SetLength(ResRec.X, 1, K+1); // 列确定为K+1，1~K对应K个变量，行可能增加；TraceToResRec(0);            // Trace表格局转化为相应的第1条结果；// Trace 表其他格局处理；while VTraceNext = 0 do      // Trace 表有新格局；begini := High(ResRec.X);SetLength(ResRec.X, i+2, K+1); // 行数增1，列数维持 K+1 不变；TraceToResRec(i+1);            // Trace表格局转化为相应的第i+1条结果；end;Result := 0;
end;

constructor TDP.Create(AProbPara: TProbPara);
vari: Integer;
begin{ 问题参数信息初始化；}FProbPara.n := AProbPara.n;SetLength(FProbPara.Ub, AProbPara.n + 1);// 分配存储，Ub[1..n] 存自变量上界，Ub[0] 空缺；SetLength(FProbPara.Lb, AProbPara.n + 1);// 分配存储，Lb[1..n] 存自变量下界，Lb[0] 空缺；FProbPara.Ub := Copy(AProbPara.Ub);FProbPara.Lb := Copy(AProbPara.Lb);FProbPara.dx := AProbPara.dx;FProbPara.f := AProbPara.f;
FUbSum := GetSum(FProbPara.Ub, 1, FProbPara.n);FLbSum := GetSum(FProbPara.Lb, 1, FProbPara.n);FM := Round((FUbSum-FLbSum) / FProbPara.Dx); // 2012-03-27,Floor改为Round；// 离散化段数（有效待分配负荷的个数）；
{ 最优解信息表FAuxTable预处理；}SetLength(FAuxTable, FM+1);for i := 0 to FM dobeginSetLength(FAuxTable[i].U, FProbPara.n + 1);// 状态变量，0空留，1~n 对应各阶段；SetLength(FAuxTable[i].X, FProbPara.n + 1);// 自变量最优值，0空留，1~n 对应各阶段；FAuxTable[i].S := i * FProbPara.Dx;// 自变量和的等式约束值，end;
end;

{ 功能：递归计算K阶段最优（最小）状态变量值；参数:K：阶段数；Sum：变量等式约束值，即阶段之间的；返回：最优K阶段状态变量值；
}
function TDP.StateRecursion(K: Integer; Sum: Double): Double;
varUk: Double;Uk1: array of Double;Xk: Double;iMin, iMax: Integer;m: Integer;Count: array of Integer;i: Integer;
beginif (FLbSum <= Sum) and (Sum <= FUbSum) then   // 验证可行性条件；beginif K > 1 then     // 变量个数大于1；beginiMin := Ceil(Max(FProbPara.Lb[k], Sum - GetSum(FProbPara.Ub,1,K-1)) /FProbPara.Dx);iMax := Floor(Min(FProbPara.Ub[k], Sum - GetSum(FProbPara.Lb,1,K-1)) /FProbPara.Dx);SetLength(Uk1, iMax - iMin + 1);// 临时存储 iMax-iMin + 1 个k-1阶段状态Uk1；for i:=iMin to iMax dobeginXk := i * FProbPara.Dx;Uk1[i-iMin] := StateRecursion(K-1, Sum-Xk);// 由Uk1作为辅助，来计算最优Uk；Uk1[i-iMin] := FProbPara.f(K, Xk) + Uk1[i-iMin];// 第k阶段状态变量，其最小的作为最优Uk；if i > iMin thenbeginif Uk1[i-iMin] < Uk thenUk := Uk1[i-iMin]; // 找到Uk[i-Min]最小值,作为Uk;endelseUk := Uk1[i-iMin];   // 以Uk1[0] 作为Uk初值；end;{ 找出Uk1[i-iMin] := FProbPara.f(K, Xk) + Uk1[i-iMin]中最小的1个或几个（等于Uk的）；}SetLength(Count, iMax-iMin+1 + 1); // Count[*] 存储Uk1[i-iMin]最小时的i值；Count[0] := 0;                     // count[0] 记录最小i值个数；for i:=iMin to iMax dobeginif Uk1[i-iMin] = Uk thenbeginCount[0] := Count[0] + 1;Count[Count[0]] := i;end;end;{ K阶段最优状态变量取Sum，对应的第k号自变量的最优取值信息，填入FAuxTable表；}m := Round((Sum-FLbSum) / FProbPara.dx);  //2012-03-27,Floor改为Round；SetLength(FAuxTable[m].X[K], Count[0]+1);FAuxTable[m].X[K, 0] := Count[0]; // 第k号自变量的取值个数；for i := 1 to Count[0] do         // 第k号自变量的取值；beginFAuxTable[m].X[K, i] := Count[i] * FProbPara.dx;end;FAuxTable[m].U[K] := Uk;endelse if K = 1 then                 // 一个变量情形；beginUk := FProbPara.f(K, Sum);       // 可行性已有算法机制保证，直接分配;m := Round((Sum-FLbSum) / FProbPara.dx);  //2012-03-27,Floor改为Round；SetLength(FAuxTable[m].X[K], 2);FAuxTable[m].X[K, 0] := 1;       // 最优解个数，此时只有一个；FAuxTable[m].X[K, 1] := Sum;     // 最优解；FAuxTable[m].U[K] := Uk;endelse        // K < 1;Exit;endelseExit;       // 不满足可行性条件；
Result := Uk;
end;

function TDP.GetSum(AnArray: array of Double; LowPos, HighPos: Integer): Double;
vari: Integer;L: Integer;H: Integer;
beginL := Low(AnArray);H := High(AnArray);if (LowPos < L) or (HighPos > H) or (LowPos > HighPos) thenExit;Result := 0;for i:= LowPos to HighPos doResult := Result + AnArray[i];
end;

end.

代码不多，但算法的实现细节，尤其是根据动态规划递归过程得到的“不太友好”的最优负荷分配表，生成优化运行总表的“格局”推演方法，还是颇精细或者说晦涩的，记得仔细捋了一两天才捋顺。源码包中，另有一个文档有详细描述，一并放在百度网盘，供参考。

下载地址：https://pan.baidu.com/s/1rzipSa0xrLBzuhY3KmxTow

提取码：dv2y