CN110070213B - 一种电-气综合能源系统的日前调度方法 - Google Patents

Abstract

本发明涉及属于对能源互联网的电力系统运行和控制技术领域，更具体地，涉及一种电‑气综合能源系统的日前调度方法。包括以下步骤：S1.建立电‑气综合能源系统日前调度模型：包括建立电力网络模型、建立天然气网络模型、建立燃气机组模型、建立非燃气机组模型、建立系统运行成本最低模型；S2.对S1步骤建立的电‑气综合能源系统日前调度模型进行约束处理；S3.对模型进行求解：首先利用罚函数求解模型获取粗糙解，然后利用泰勒级数法求解模型，并采用序贯二阶锥规划形成迭代优化。本发明将罚函数法与泰勒级数法结合，利用序贯二阶锥规划方法，设计了一套电‑气综合能源系统的日前调度算法，该算法具有较好的收敛性与精确性，提高了计算时间。

Description

一种电-气综合能源系统的日前调度方法

技术领域

本发明涉及属于对能源互联网的电力系统运行和控制技术领域，更具体地，涉及一种电-气综合能源系统的日前调度方法。

背景技术

作为第三次工业革命的核心技术，能源互联网力图结合可再生能源技术与互联网技术，推动分布式可再生能源的大规模利用与分享，促进电力、交通、天然气等多种复杂网络系统的相互融合，最终实现改变能源利用模式，推动经济与社会可持续发展的目的。其中天然气相比于其他一次能源，经济环保、储量丰富，联合循环燃气轮机发电比重近年来逐步提升，电力与天然气系统的融合日益深厚。能源互联网中，以电力系统与天然气系统为核心的“电-气”综合能源系统应用十分广泛，而在日前调度中同时考虑电力网络与天然气网络有助于提高系统稳定性与可靠性，降低整体运营成本，同时也推动了综合能源系统的持续健康发展。

“电-气”综合能源系统的日前调度分为两部分，电力系统日前调度与天然气系统日前调度。

在电力系统日前调度中，往往考虑网络安全约束的机组组合优化或经济调度，而前者在保证电力系统安全性的前提下，以运行成本最低为目标。机组组合模型是一个混合整数规划(MIP)模型，其中整数变量来源于机组的开停状态，而机组组合模型中的发电机成本曲线可以是线性形式，也可以是非线性的二次函数形式。机组组合常用的求解方法是通过成熟的商用求解器，求解该混合整数规划问题。

在天然气系统日前调度中，往往以天然气网络稳态潮流优化问题为主，在保证管道与节点气压不越限的前提下，以运行成本最低为目标。天然气网络稳态潮流模型是一个非线性非凸的模型，其中管道-气压规则由Weymouth方程表示，该方程是一个非线性非凸的方程，这也是求解天然气网络稳态潮流的难点。常用的天然气网络稳态潮流的求解方法包括内点法、分段线性化等。

传统的电力系统日前调度方法，大多把天然气网络约束造成的影响简单化，并反映到燃气轮机上，进而体现到电力系统日前调度中。然而这并不准确，燃气轮机在“电-气”综合能源系统中属于耦合部件，不能完全体现天然气网络对电力网络调度造成的影响。此外，现有“电-气”综合能源系统的模型不够完善，求解方法的效率与精度较低，会造成额外的调度成本。

发明内容

本发明为克服上述现有技术所述的至少一种缺陷，提供一种电-气综合能源系统的日前调度方法，将罚函数法与泰勒级数法结合，利用序贯二阶锥规划方法，设计了一套电-气综合能源系统的日前调度算法，该算法具有较好的收敛性与精确性，提高了计算时间。

为解决上述技术问题，本发明采用的技术方案是：一种电-气综合能源系统的日前调度方法，包括以下步骤：

S1.建立电-气综合能源系统日前调度模型：包括建立电力网络模型、建立天然气网络模型、建立燃气机组模型、建立非燃气机组模型、建立系统运行成本最低模型；

S2.对S1步骤建立的电-气综合能源系统日前调度模型进行约束处理；

S3.对模型进行求解：首先利用罚函数求解模型获取粗糙解，然后利用泰勒级数法求解模型，并采用序贯二阶锥规划形成迭代优化。

进一步的，所述的S1步骤中，所述的建立电力网络模型包括：

1)建立电力平衡方程：

2)建立发电机组出力约束方程：

3)建立机组爬坡约束方程：

4)建立机组最小开停时间约束方程：

5)建立机组启停状态逻辑约束方程：

y_i,t-z_i,t＝u_i,t-u_i,t-1 公式(7)

y_i,t+z_i,t≤1 公式(8)

6)建立启动与停止成本约束方程：

SU_i,t＝su_iy_i,t 公式(9)

SD_i,t＝sd_iz_i,t 公式(10)

7)建立直流潮流约束方程：

E·pf＝C·P_i,t-L·D_b,t 公式(11)

8)建立线路传输容量约束方程：

θ_ref＝0 公式(14)

式中：P_i，t指机组i在t时刻的出力；D_b，t指节点b在t时刻的负荷；

指机组i的最小与最大出力；u_i，t指机组i在t时刻的启停状态变量(1为开，0为关)；y_i，t、z_i，t指机组i在t时刻的启动、停止动作状态变量；UR_i、DR_i指机组i每小时上、下爬坡的最大值；

指机组i在t-1时刻启动、停止的时间长；

指机组i的最小启动、停止时间；su_i、sd_i指机组i的启动、停止费用常数；SU_i，t、SD_i，t指机组i的启动、停止费用；E、C、L分别指节点-线路关联矩阵、节点-发电机组关联矩阵、节点-负荷关联矩阵；θ_j指节点j的电压相角；x_j，k指节点j到节点k的线路电抗；pf_br为线路br的潮流；

为线路br的潮流最大值；θ_ref为参考节点的电压相角；

所述的建立天然气网络模型具体包括：

1)建立Weymouth约束方程：

2)建立节点气压上下限约束方程：

3)建立管道容量约束方程：

4)建立节点流量平衡约束方程：

T^sW^s-T^dW^d-T^gW^g＝T^pF^p 公式(18)

5)建立气源出气量上下限约束方程：

式中，t在此指时间，

指节点m到节点n的天然气管道流量；sgn(·)是一个信号函数；C_mn是管道的固定参数，与环境温度等外界因素有关；NN表示天然气网络节点集合；ω_m指节点气压；

分别指节点气压的上下限；

分别指天然气管道流量的上下限；T^s、T^d、T^g、T^p分别指节点-气源关联矩阵、节点-非燃气轮机气负荷关联矩阵、节点-燃气轮机气负荷关联矩阵、节点-管道关联矩阵；W^s、W^d分别指气源的出力和非燃气轮机气负荷；W^s，max、W^s，min分别为气源出力的上、下限；GW指气源集合；

所述的建立燃气机组模型如以下约束方程：

式中，

指燃气机组i在时刻t消耗的天然气量，

指燃气机组i的耗量系数，GU指燃气机组集合；

所述的建立非燃气机组模型如以下约束方程：

式中，

指非燃气机组的成本函数，

指非燃气机组的成本系数；

所述的建立系统运行成本最低目标模型如以下方程：

式中，λ_i指天然气网络中气源的成本系数。

进一步的，所述的S2步骤具体包括以下步骤：

S21.燃气机组模型约束处理：根据S1步骤中建立的综合能源日前调度模型，燃气机组模型作为天然气网络与电力网络的耦合部分，其耗气量方程即公式(20)可以转换为二阶锥规划松弛约束方程，如公式(23)所示，并且由于成本最小的目标函数限制，约束方程公式(23)保持紧致；约束方程如下所示：

S22.天然气网络模型约束处理：对天然气网络模型约束的处理首先引入一组整数变量表示天然气管道潮流的方向，将约束公式(15)转换为约束公式(24)表示，约束公式(17)转换为公式(25)表示，并增加公式(26)和公式(27)两个约束条件；公式(24)至公式(27)如下所示：

式中，π_m指节点气压的平方值，

为一组整数变量用于表示天然气管道潮流的方向；

S23.约束公式(24)符合二阶锥的基本形式，将公式(24)转换为如下约束公式(28)至公式(32)：

式中，Γ_mn为一组辅助变量，在物理上的含义为管道两端的节点气压平方差；在这一步中，约束公式(28)是一个不等式约束，当其保持紧致的时候才会与约束公式(24)相等。

进一步的，所述的S3步骤具体包括：

S31.罚函数法求解获取初值：在目标函数增加一个罚函数项，具体如公式(33)所示：

其中

为极小的正常数，在目标函数公式(33)与约束公式(1)至公式(14)、公式(16)、公式(18)、公式(19)、公式(21)、公式(23)、公式(25)至公式(32)下求解模型；

该步骤可以收紧约束公式(28)，但是只能得到一个较为粗糙的天然气网络潮流解即

作为接下来算法的初始点，并增加一个辅助松弛变量S_mn，令其初值为0，该步骤为迭代的初始值获取，因此k＝0；

S32.收敛参数设置：设置惩罚因子φ⁰，惩罚因子增长率v>1，最大惩罚因子φ^max，收敛残差值ξ^Z与ξ^S，同时，通过增加一个约束条件即公式(34)将约束条件公式(28)收紧，使其保持松弛的严格性；公式(34)如下所示：

S33.利用泰勒级数展开的数学定理，将约束公式(34)转换为约束公式(35)：

S34.在目标函数公式(36)与约束公式(1)至公式(14)、公式(16)、公式(18)、公式(19)、公式(21)、公式(23)、公式(25)至公式(32)、公式(35)下求解模型：

S35.收敛判断，具体判断方法如公式37和公式38所示：

gas^z＝z^k-z^k-1≤ξ^Z 公式(37)

S36.如果满足收敛条件，那么算法结束，得到电-气综合能源系统的日前调度结果；如果不满足收敛条件，那么更新惩罚因子和迭代步数如公式(39)和公式(40)所示，并重复步骤S34至S35；

φ^k＝min(vφ^k-1,φ^max) 公式(39)

k＝k+1 公式(40)。

与现有技术相比，有益效果是：

1.本发明将罚函数法与泰勒级数法结合，利用序贯二阶锥规划方法，设计了一套“电-气”综合能源系统的日前调度算法，该算法具有较好的收敛性与精确性，并提高了计算时间，提高“电-气”综合能源系统管理者管理效率；

2.在“电-气”综合能源系统领域内，针对本发明的优化目标进行了模型的处理与转换，将燃气机组模型与天然气网络模型转换成便于求解与表达的形式，为新算法提供更好的效果和针对性；

3.与传统的“电-气”综合能源系统的日前调度方法相比，本发明在不简化模型的前提下，通过适当的处理将模型转换为面向本发明所提方法的“电-气”综合能源系统模型。其中结合罚函数法与泰勒级数法的算法，吸收了各自的优点，具有快速收敛、精度较高等特点；

4.本发明面向“电-气”综合能源系统模型所提的日前调度方法，具有普适性，在保证系统安全性的前提下，使得系统运营者效益更优。

附图说明

图1是本发明方法流程图。

具体实施方式

附图仅用于示例性说明，不能理解为对本发明的限制；为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。附图中描述位置关系仅用于示例性说明，不能理解为对本发明的限制。

如图1所示，一种电-气综合能源系统的日前调度方法，包括以下步骤：

步骤1.建立电-气综合能源系统日前调度模型：包括建立电力网络模型、建立天然气网络模型、建立燃气机组模型、建立非燃气机组模型、建立系统运行成本最低模型；

(1)建立电力网络模型包括：

1)建立电力平衡方程：

2)建立发电机组出力约束方程：

3)建立机组爬坡约束方程：

4)建立机组最小开停时间约束方程：

5)建立机组启停状态逻辑约束方程：

y_i,t-z_i,t＝u_i,t-u_i,t-1 公式(7)

y_i,t+z_i,t≤1 公式(8)

6)建立启动与停止成本约束方程：

SU_i,t＝su_iy_i,t 公式(9)

SD_i,t＝sd_iz_i,t 公式(10)

7)建立直流潮流约束方程：

E·pf＝C·P_i,t-L·D_b,t 公式(11)

8)建立线路传输容量约束方程：

θ_ref＝0 公式(14)

式中：P_i，t指机组i在t时刻的出力；D_b，t指节点b在t时刻的负荷；

指机组i的最小与最大出力；u_i，t指机组i在t时刻的启停状态变量(1为开，0为关)；y_i，t、z_i，t指机组i在t时刻的启动、停止动作状态变量；UR_i、DR_i指机组i每小时上、下爬坡的最大值；

指机组i在t-1时刻启动、停止的时间长；

指机组i的最小启动、停止时间；su_i、sd_i指机组i的启动、停止费用常数；SU_i，t、SD_i，t指机组i的启动、停止费用；E、C、L分别指节点-线路关联矩阵、节点-发电机组关联矩阵、节点-负荷关联矩阵；θ_j指节点j的电压相角；x_j，k指节点j到节点k的线路电抗；pf_br为线路br的潮流；

为线路br的潮流最大值；θ_ref为参考节点的电压相角。

(2)建立天然气网络模型具体包括：

1)建立Weymouth约束方程：

2)建立节点气压上下限约束方程：

3)建立管道容量约束方程：

4)建立节点流量平衡约束方程：

T^sW^s-T^dW^d-T^gW^g＝T^pF^p 公式(18)

5)建立气源出气量上下限约束方程：

式中，t在此指时间，

指节点m到节点n的天然气管道流量；sgn(·)是一个信号函数；C_mn是管道的固定参数，与环境温度等外界因素有关；NN表示天然气网络节点集合；ω_m指节点气压；

分别指节点气压的上下限；

分别指天然气管道流量的上下限；T^s、T^d、T^g、T^p分别指节点-气源关联矩阵、节点-非燃气轮机气负荷关联矩阵、节点-燃气轮机气负荷关联矩阵、节点-管道关联矩阵；W^s、W^d分别指气源的出力和非燃气轮机气负荷；W^s，max、W^s，min分别为气源出力的上、下限；GW指气源集合。

(3)建立燃气机组模型如以下约束方程：

式中，

指燃气机组i在时刻t消耗的天然气量，

指燃气机组i的耗量系数，GU指燃气机组集合。

(4)建立非燃气机组模型如以下约束方程：

式中，

指非燃气机组的成本函数，

指非燃气机组的成本系数。

(5)建立系统运行成本最低目标模型如以下方程：

式中，λ_i指天然气网络中气源的成本系数。

步骤2.对S1步骤建立的电-气综合能源系统日前调度模型进行约束处理；

具体包括以下步骤：

S21.燃气机组模型约束处理：根据S1步骤中建立的综合能源日前调度模型，燃气机组模型作为天然气网络与电力网络的耦合部分，其耗气量方程即公式(20)可以转换为二阶锥规划松弛约束方程，如公式(23)所示，并且由于成本最小的目标函数限制，约束方程公式(23)保持紧致；约束方程如下所示：

S22.天然气网络模型约束处理：对天然气网络模型约束的处理首先引入一组整数变量表示天然气管道潮流的方向，将约束公式(15)转换为约束公式(24)表示，约束公式(17)转换为公式(25)表示，并增加公式(26)和公式(27)两个约束条件；公式(24)至公式(27)如下所示：

式中，π_m指节点气压的平方值，

为一组整数变量用于表示天然气管道潮流的方向；

S23.约束公式(24)符合二阶锥的基本形式，将公式(24)转换为如下约束公式(28)至公式(32)：

式中，Γ_mn为一组辅助变量，在物理上的含义为管道两端的节点气压平方差；在这一步中，约束公式(28)是一个不等式约束，当其保持紧致的时候才会与约束公式(24)相等。

步骤3.对模型进行求解：通过上述对模型的处理可以将原本的非线性非凸的“电-气”综合能源系统模型转换为混合整数二阶锥模型，但是需要收紧约束公式(28)。通过罚函数、泰勒分段线性化的手段将约束公式(28)收紧，并且采用序贯二阶锥规划形成迭代优化。

具体包括以下步骤：

S31.罚函数法求解获取初值：在目标函数增加一个罚函数项，具体如公式(33)所示：

其中

为极小的正常数，具体范围小于0.1^4，在目标函数公式(33)与约束公式(1)至公式(14)、公式(16)、公式(18)、公式(19)、公式(21)、公式(23)、公式(25)至公式(32)下求解模型；

该步骤可以收紧约束公式(28)，但是只能得到一个较为粗糙的天然气网络潮流解即

作为接下来算法的初始点，并增加一个辅助松弛变量S_mn，令其初值为0，该步骤为迭代的初始值获取，因此k＝0；

S32.收敛参数设置：设置惩罚因子φ⁰，惩罚因子增长率v>1，最大惩罚因子φ^max，收敛残差值ξ^Z与ξ^S，同时，通过增加一个约束条件即公式(34)将约束条件公式(28)收紧，使其保持松弛的严格性；公式(34)如下所示：

S33.利用泰勒级数展开的数学定理，将约束公式(34)转换为约束公式(35)：

S34.在目标函数公式(36)与约束公式(1)至公式(14)、公式(16)、公式(18)、公式(19)、公式(21)、公式(23)、公式(25)至公式(32)、公式(35)下求解模型：

S35.收敛判断，具体判断方法如公式37和公式38所示：

gas^z＝z^k-z^k-1≤ξ^Z 公式(37)

S36.如果满足收敛条件，那么算法结束，得到电-气综合能源系统的日前调度结果；如果不满足收敛条件，那么更新惩罚因子和迭代步数如公式(39)和公式(40)所示，并重复步骤S34至S35；

φ^k＝min(vφ^k-1,φ^max) 公式(39)

k＝k+1 公式(40)。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (3)

1.一种电-气综合能源系统的日前调度方法，其特征在于，包括以下步骤：

S1.建立电-气综合能源系统日前调度模型：包括建立电力网络模型、建立天然气网络模型、建立燃气机组模型、建立非燃气机组模型、建立系统运行成本最低模型；

S2.对S1步骤建立的电-气综合能源系统日前调度模型进行约束处理；

S3.对模型进行求解：首先利用罚函数求解模型获取粗糙解，然后利用泰勒级数法求解模型，并采用序贯二阶锥规划形成迭代优化；

其中，S1步骤中，所述的建立电力网络模型包括：

1)建立电力平衡方程：

2)建立发电机组出力约束方程：

P_i ^minu_i,t≤P_i,t≤P_i ^maxu_i,t 公式(2)

3)建立机组爬坡约束方程：

P_i,t-P_i,t-1≤UR_i(1-y_i,t)+P_i ^miny_i,t 公式(3)

P_i,t-1-P_i,t≤DR_i(1-z_i,t)+P_i ^minz_i,t 公式(4)

4)建立机组最小开停时间约束方程：

5)建立机组启停状态逻辑约束方程：

y_i,t-z_i,t＝u_i,t-u_i,t-1 公式(7)

y_i,t+z_i,t≤1 公式(8)

6)建立启动与停止成本约束方程：

SU_i,t＝su_iy_i,t 公式(9)

SD_i,t＝sd_iz_i,t 公式(10)

7)建立直流潮流约束方程：

E·pf＝C·P_i,t-L·D_b,t 公式(11)

8)建立线路传输容量约束方程：

θ_ref＝0 公式(14)

式中：P_i，t指机组i在t时刻的出力；D_b，t指节点b在t时刻的负荷；

指机组i的最小与最大出力；u_i，t指机组i在t时刻的启停状态变量(1为开，0为关)；y_i，t、z_i，t指机组i在t时刻的启动、停止动作状态变量；UR_i、DR_i指机组i每小时上、下爬坡的最大值；

指机组i在t-1时刻启动、停止的时间长；

指机组i的最小启动、停止时间；su_i、sd_i指机组i的启动、停止费用常数；SU_i，t、SD_i，t指机组i的启动、停止费用；E、C、L分别指节点-线路关联矩阵、节点-发电机组关联矩阵、节点-负荷关联矩阵；θ_j指节点j的电压相角；x_j，k指节点j到节点k的线路电抗；pf_br为线路br的潮流；

为线路br的潮流最大值；θ_ref为参考节点的电压相角；

所述的建立天然气网络模型具体包括：

1)建立Weymouth约束方程：

2)建立节点气压上下限约束方程：

3)建立管道容量约束方程：

4)建立节点流量平衡约束方程：

T^sW^s-T^dW^d-T^gW^g＝T^pF^p 公式(18)

5)建立气源出气量上下限约束方程：

式中，t在此指时间，

指节点m到节点n的天然气管道流量；sgn(·)是一个信号函数；C_mn是管道的固定参数，与环境温度等外界因素有关；NN表示天然气网络节点集合；ω_m指节点气压；

分别指节点气压的上下限；

分别指天然气管道流量的上下限；T^s、T^d、T^g、T^p分别指节点-气源关联矩阵、节点-非燃气轮机气负荷关联矩阵、节点-燃气轮机气负荷关联矩阵、节点-管道关联矩阵；W^s、W^d分别指气源的出力和非燃气轮机气负荷；W^s.max、W^s.min分别为气源出力的上、下限；GW指气源集合；

所述的建立燃气机组模型如以下约束方程：

式中，

指燃气机组i在时刻t消耗的天然气量，

指燃气机组i的耗量系数，GU指燃气机组集合；

所述的建立非燃气机组模型如以下约束方程：

式中，

指非燃气机组的成本函数，

指非燃气机组的成本系数；

所述的建立系统运行成本最低目标模型如以下方程：

式中，λ_i指天然气网络中气源的成本系数。

2.根据权利要求1所述的一种电-气综合能源系统的日前调度方法，其特征在于，所述的S2步骤具体包括以下步骤：

S21.燃气机组模型约束处理：根据S1步骤中建立的综合能源日前调度模型，燃气机组模型作为天然气网络与电力网络的耦合部分，其耗气量方程即公式(20)可以转换为二阶锥规划松弛约束方程，如公式(23)所示，并且由于成本最小的目标函数限制，约束方程公式(23)保持紧致；约束方程如下所示：

S22.天然气网络模型约束处理：对天然气网络模型约束的处理首先引入一组整数变量表示天然气管道潮流的方向，将约束公式(15)转换为约束公式(24)表示，约束公式(17)转换为公式(25)表示，并增加公式(26)和公式(27)两个约束条件；公式(24)至公式(27)如下所示：

式中，π_m指节点气压的平方值，

为一组整数变量用于表示天然气管道潮流的方向；

S23.约束公式(24)符合二阶锥的基本形式，将公式(24)转换为如下约束公式(28)至公式(32)：

式中，Γ_mn为一组辅助变量，在物理上的含义为管道两端的节点气压平方差；在这一步中，约束公式(28)是一个不等式约束，当其保持紧致的时候才会与约束公式(24)相等。

3.根据权利要求2所述的一种电-气综合能源系统的日前调度方法，其特征在于，所述的S3步骤具体包括：

S31.罚函数法求解获取初值：在目标函数增加一个罚函数项，具体如公式(33)所示：

其中

为极小的正常数，在目标函数公式(33)与约束公式(1)至公式(14)、公式(16)、公式(18)、公式(19)、公式(21)、公式(23)、公式(25)至公式(32)下求解模型；

该步骤可以收紧约束公式(28)，但是只能得到一个较为粗糙的天然气网络潮流解即

作为接下来算法的初始点，并增加一个辅助松弛变量S_mn，令其初值为0，该步骤为迭代的初始值获取，因此k＝0；

S32.收敛参数设置：设置惩罚因子φ⁰，惩罚因子增长率v＞1，最大惩罚因子φ^max，收敛残差值ξ^Z与ξ^S，同时，通过增加一个约束条件即公式(34)将约束条件公式(28)收紧，使其保持松弛的严格性；公式(34)如下所示：

S33.利用泰勒级数展开的数学定理，将约束公式(34)转换为约束公式(35)：

S34.在目标函数公式(36)与约束公式(1)至公式(14)、公式(16)、公式(18)、公式(19)、公式(21)、公式(23)、公式(25)至公式(32)、公式(35)下求解模型：

S35.收敛判断，具体判断方法如公式37和公式38所示：

gas^z＝z^k-z^k-1≤ξ^Z 公式(37)

S36.如果满足收敛条件，那么算法结束，得到电-气综合能源系统的日前调度结果；如果不满足收敛条件，那么更新惩罚因子和迭代步数如公式(39)和公式(40)所示，并重复步骤S34至S35；

φ^k＝min(vφ^k-1,φ^max) 公式(39)

k＝k+1 公式(40)。

Images