CN112711862B - 基于并发l-bfgs算法的水文模型参数率定方法及系统 - Google Patents
基于并发l-bfgs算法的水文模型参数率定方法及系统 Download PDFInfo
- Publication number
- CN112711862B CN112711862B CN202110056440.9A CN202110056440A CN112711862B CN 112711862 B CN112711862 B CN 112711862B CN 202110056440 A CN202110056440 A CN 202110056440A CN 112711862 B CN112711862 B CN 112711862B
- Authority
- CN
- China
- Prior art keywords
- model parameters
- groups
- model
- multiple groups
- main processor
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 66
- 238000000034 method Methods 0.000 title claims abstract description 66
- 238000012545 processing Methods 0.000 claims abstract description 10
- 230000008569 process Effects 0.000 claims description 31
- 238000004088 simulation Methods 0.000 claims description 31
- 230000008859 change Effects 0.000 claims description 23
- 238000004364 calculation method Methods 0.000 abstract description 14
- 230000007246 mechanism Effects 0.000 abstract description 3
- 239000011159 matrix material Substances 0.000 description 26
- 230000006870 function Effects 0.000 description 23
- 238000010586 diagram Methods 0.000 description 9
- 238000009795 derivation Methods 0.000 description 8
- 238000005516 engineering process Methods 0.000 description 5
- 125000004122 cyclic group Chemical group 0.000 description 3
- 241000282414 Homo sapiens Species 0.000 description 2
- 238000004590 computer program Methods 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000007726 management method Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 238000013473 artificial intelligence Methods 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 230000003628 erosive effect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000008020 evaporation Effects 0.000 description 1
- 238000001704 evaporation Methods 0.000 description 1
- 239000000835 fiber Substances 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000008595 infiltration Effects 0.000 description 1
- 238000001764 infiltration Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 238000013439 planning Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000002054 transplantation Methods 0.000 description 1
- 238000003911 water pollution Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F9/00—Arrangements for program control, e.g. control units
- G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
- G06F9/46—Multiprogramming arrangements
- G06F9/50—Allocation of resources, e.g. of the central processing unit [CPU]
- G06F9/5005—Allocation of resources, e.g. of the central processing unit [CPU] to service a request
- G06F9/5027—Allocation of resources, e.g. of the central processing unit [CPU] to service a request the resource being a machine, e.g. CPUs, Servers, Terminals
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2209/00—Indexing scheme relating to G06F9/00
- G06F2209/50—Indexing scheme relating to G06F9/50
- G06F2209/5018—Thread allocation
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Software Systems (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本申请提供了一种基于并发L‑BFGS算法的水文模型参数率定方法及系统。该方案应用于异构平台，包括：主处理器生成水文模型的多组模型参数并发送到协处理器；协处理器的不同线程分别获得各自划分的模型参数，计算获得模型参数对应的纳什效率系数并返回主处理器；主处理器判断多组模型参数对应的纳什效率系数是否满足预设条件，当满足时，输出所述更新后的多组模型参数，当不满足时，获得多组模型参数对应的梯度值，并利用并发L‑BFGS算法获得更新后的多组模型参数再次率定，直至满足结束条件。本申请的方案利用协处理的多线程并发机制，提高了水文模型参数率定的效率。
Description
技术领域
本申请涉及水文模拟技术领域，尤其涉及一种基于并发L-BFGS算法的水文模型参数率定方法及系统。
背景技术
水文模型是对大自然中复杂的水文现象在计算机中的抽象化和概念化的展现，是人们认识复杂的水循环运动过程和机制下的有效方法和手段，同时也是处理各种复杂的水文问题的重要工具。水文模型由众多的物理模型构成，如降雨入渗、蒸发散发、产流回流等具有特定物理意义的物理模型，其中大多数物理模型都包含一系列流域相关的待定参数，而对水文模型进行参数率定就是通过确定当前所研究流域的待定参数值，来完成对当前流域的水文建模。因此对模型参数率定是水文模拟过程的重要组成部分，参数选取的优劣直接影响到最终的模拟精度。
如今在大数据和云计算的发展推动下，用于构建水文模型的研究正逐步进入一个新的时代，使得水文建模无论在时间上还是空间上都有不断延伸和扩展的趋势。此外，人类对水文过程的认识水平也在不断地提高，加上不断增长的研究和发展需求，使得水文过程涉及到的物理、化学领域发展过程的叙述也更加全面。在这几大因素的推动下，使得水文模型运行所需的计算水平和时间需求也不断增长。
在采用水文模型进行水文模拟和预报时，水文工作者主要的任务是去寻找水文模型的最优参数集，水文模型模拟的好坏程度主要就是依赖于水文模型的参数集的选择。不同的流域包含自己的模型参数，要得到一套适合当前流域的参数集合，传统的方式要么依靠经验，要么通过率定算法(如粒子群算法和遗传算法等优化算法)进行搜索，而参数空间维度通常高达十几、二十维，高维参数空间搜索计算量巨大、耗时长，传统搜索方式已经难以实施。
发明内容
本申请实施例提供了一种基于并发L-BFGS算法的水文模型参数率定方法及系统，结合并发的L-BFGS算法和异构平台对水文模型参数进行率定，提高了参数率定的效率，解决了现有技术中计算量大、耗时长的问题。
第一方面，本申请提供一种基于并发L-BFGS算法的水文模型参数率定方法，应用于异构平台，异构平台包括：主处理器和协处理器，该方法包括：
步骤S1.所述主处理器生成水文模型的多组模型参数，并发送至所述协处理器；
步骤S2.所述协处理器将所述多组模型参数划分至不同的线程，各线程分别根据划分的模型参数运行所述水文模型对当前流域进行水文模拟，以及根据所述划分的模型参数对应的当前流域的水文模拟值，获得所述划分的模型参数对应的纳什效率系数，并发送至所述主处理器；
步骤S3.所述主处理器判断所述多组模型参数对应的纳什效率系数是否满足预设条件，当满足时，所述主处理器输出所述更新后的多组模型参数，模型率定结束，当不满足时，进入步骤S4；
步骤S4.所述主处理器根据所述多组模型参数对应的纳什效率系数，获得所述多组模型参数对应的梯度值，利用并发L-BFGS算法获得更新后的多组模型参数，并发送至所述协处理器，返回步骤S2。
由上，本申请利用异构平台中协处理器的并发机制，对水文模型进行多组参数率定，提高了水文模型率定的效率，减少了计算时间。
在一个可能的实施例中，所述主处理器根据所述多组模型参数对应的纳什效率系数，获得所述多组模型参数对应的梯度值包括：
所述主处理根据所述多组模型参数和预设的差商步长，获得叠加所述差商步长后的多组模型参数，并通过所述协处理器获得所述叠加所述差商步长后的多组模型参数对应的纳什效率系数；
所述主处理根据所述多组模型参数对应的纳什效率系数、所述叠加所述差商步长后的多组模型参数对应的纳什效率系数以及所述差商步长，获得所述多组模型参数对应的梯度值。
在一个可能的实施例中，所述利用并发L-BFGS算法获得更新后的多组模型参数包括：
根据当前更新前的多组模型参数对应的梯度值、以及获取的当前更新前的Z次历史更新后的多组模型参数的变化量和Z次历史更新后的多组模型参数对应的梯度值变化量，获得当前更新前的多组模型参数对应的搜索方向；
根据当前更新前的多组模型参数及其对应的搜索方向，获得当前更新后的多组模型参数；
其中，所述Z次历史更新后的多组模型参数的变化量根据所述Z次历史更新前的多组模型参数与所述Z次历史更新后的多组模型参数获得，所述Z次历史更新后的多组模型参数对应的梯度值变化量根据所述Z次历史更新前的多组模型参数对应的梯度值与所述Z次历史更新后的多组模型参数对应的梯度值获得。
在一个可能的实施例中，所述根据当前更新前的多组模型参数对应的梯度值、以及获取的当前更新前的Z次历史更新后的多组模型参数的变化量和Z次历史更新后的多组模型参数对应的梯度值变化量，获得当前更新前的多组模型参数对应的搜索方向包括：
根据所述多组模型参数对应的梯度值、所述Z次历史更新后的多组模型参数的变化量和所述Z次历史更新后的多组模型参数对应的梯度值变化量，进行后向循环，获得第Z次历史更新前的多组模型参数对应的搜索方向；
根据所述Z次历史更新前的多组模型参数及其对应的梯度值、所述Z次历史更新后的多组模型参数的变化量和所述Z次历史更新后的多组模型参数对应的梯度值变化量，进行前向循环，获得当前更新前的多组模型参数对应的搜索方向。
由上，本申请利用Z次历史中模型参数和梯度值计算当前的搜索方向，减少了方案中存储的数据量，在保证精度的前提下，提高了计算效率。
在一个可能的实施例中，所述根据所述多组模型参数及其对应的搜索方向，获得更新后的多组模型参数包括：
在实数范围内选择所述多组模型参数对应的多个步长因子；
根据所述多组模型参数对应的搜索方向和多个步长因子，获得多组候选的多组模型参数，并通过所述协处理器获得各组候选的多组模型参数对应的纳什效率系数；
根据所述各组候选的多组模型参数对应纳什效率系数，确定所述更新后的多组模型参数。
第二方面，本申请还提供一种基于并发L-BFGS算法的水文模型参数率定系统，该系统包括：
主处理器，用于生成水文模型的多组模型参数，并发送至所述协处理器；
协处理器，用于将所述多组模型参数划分至不同的线程，各线程分别根据划分的模型参数运行所述水文模型对当前流域进行水文模拟，以及根据所述划分的模型参数对应的当前流域的水文模拟值，获得所述划分的模型参数对应的纳什效率系数，并发送至所述主处理器；
主处理器，还用于判断所述多组模型参数对应的纳什效率系数是否满足预设条件，当满足时，所述主处理器输出所述更新后的多组模型参数，算法结束，当不满足时，所述主处理器根据所述多组模型参数对应的纳什效率系数，获得所述多组模型参数对应的梯度值，利用并发L-BFGS算法获得更新后的多组模型参数，并发送至所述协处理器。
在一个可能的实施例中，所述主处理器具体用于：
所述主处理根据所述多组模型参数和预设的差商步长，获得叠加所述差商步长后的多组模型参数，并通过所述协处理器获得所述叠加所述差商步长后的多组模型参数对应的纳什效率系数；
所述主处理根据所述多组模型参数对应的纳什效率系数、所述叠加所述差商步长后的多组模型参数对应的纳什效率系数以及所述差商步长，获得所述多组模型参数对应的梯度值。
在一个可能的实施例中，所述主处理器还具体用于：
根据当前更新前的多组模型参数对应的梯度值、以及获取的当前更新前的Z次历史更新后的多组模型参数的变化量和Z次历史更新后的多组模型参数对应的梯度值变化量，获得当前更新前的多组模型参数对应的搜索方向；
根据当前更新前的多组模型参数及其对应的搜索方向，获得当前更新后的多组模型参数；
其中，所述Z次历史更新后的多组模型参数的变化量根据所述Z次历史更新前的多组模型参数与所述Z次历史更新后的多组模型参数获得，所述Z次历史更新后的多组模型参数对应的梯度值变化量根据所述Z次历史更新前的多组模型参数对应的梯度值与所述Z次历史更新后的多组模型参数对应的梯度值获得。
在一个可能的实施例中，所述主处理器还具体用于：
根据所述多组模型参数对应的梯度值、所述Z次历史更新后的多组模型参数的变化量和所述Z次历史更新后的多组模型参数对应的梯度值变化量，进行后向循环，获得第Z次历史更新前的多组模型参数对应的搜索方向；
根据所述Z次历史更新前的多组模型参数及其对应的梯度值、所述Z次历史更新后的多组模型参数的变化量和所述Z次历史更新后的多组模型参数对应的梯度值变化量，进行前向循环，获得当前更新前的多组模型参数对应的搜索方向。
在一个可能的实施例中，所述主处理器还具体用于：
在实数范围内选择所述多组模型参数对应的多个步长因子；
根据所述多组模型参数对应的搜索方向和多个步长因子，获得多组候选的多组模型参数，并通过所述协处理器获得各组候选的多组模型参数对应的纳什效率系数；
根据所述各组候选的多组模型参数对应纳什效率系数，确定所述更新后的多组模型参数。
附图说明
图1是本申请提供的L-BFGS算法的执行流程图；
图2是本申请实施例提供的将L-BFGS应用于水文模型进行参数率定的原理示意图；
图3是本申请实施例提供的并发L-BFGS算法的原理示意图；
图4是本申请提供的水文模型的工作流程图；
图5是本申请实施例提供的水文模型在异构平台的工作流程图；
图6是本申请实施例提供的并发L-BFGS应用于异构平台对水文模型进行参数率定的流程图；
图7是本申请实施例提供的一种水文模型参数率定系统结构示意图。
具体实施方式
为了使本申请实施例的目的、技术方案和优点更加清楚，下面将结合附图，对本申请实施例中的技术方案进行描述。
在本申请实施例的描述中，“示例性的”、“例如”或者“举例来说”等词用于表示作例子、例证或说明。本申请实施例中被描述为“示例性的”、“例如”或者“举例来说”的任何实施例或设计方案不应被解释为比其它实施例或设计方案更优选或更具优势。确切而言，使用“示例性的”、“例如”或者“举例来说”等词旨在以具体方式呈现相关概念。
在本申请实施例的描述中，术语“和/或”，仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，A和/或B，可以表示：单独存在A，单独存在B，同时存在A和B这三种情况。另外，除非另有说明，术语“多个”的含义是指两个或两个以上。例如，多个系统是指两个或两个以上的系统，多个屏幕终端是指两个或两个以上的屏幕终端。
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。术语“包括”、“包含”、“具有”及它们的变形都意味着“包括但不限于”，除非是以其他方式另外特别强调。
图1示出了L-BFGS算法的执行流程图。如图1所示，传统的L-BFGS算法是对初始化的参数进行循环求解，判断求解结果是否满足循环终止条件。满足时，结束算法循环；不满足时，拟合海森矩阵(或海森逆矩阵)，从而更新迭代的搜索方向、模型参数等变量，直至满足终止条件，获得最优的参数。
传统的L-BFGS算法由于采用了线搜索策略，使算法具有较强的全局收敛能力。因此，将其应用在水文模型的参数率定过程中，可以得到较好的率定结果。
图2示出了本申请实施例将L-BFGS应用于水文模型进行参数率定的原理示意图。如图2所示，该过程包括：
首先初始化模型参数，并代入水文模型；
将历史水文数据代入水文模型，获得模型输出的模拟值，将水文模型输出的模拟值及水文数据的实际值代入Nash公式，获得Nash系数；
计算梯度值，将Nash系数取反后代入L-BFGS算法、以及将模型参数对应的梯度值代入算法进行循环求解；
判断是否满足循环终止条件，当不满足时，拟合海森矩阵，以更新搜索方向、模型参数及其Nash系数；循环此过程，即可得到水文模型的最优参数及其Nash系数。
然而，由于水文模型参数的不确定性，一次参数率定可能无法得到最优参数组，这样就需要水文工作者进行多次参数率定，才能得到较好的率定结果。不足的是，运用现有技术对高维参数进行多次率定时，计算量巨大，时间长，无法满足实际需求。高性能的异构并行计算技术为水文模型需多次率定的问题提供了解决方案，同时也为大尺度水文模拟提供了可能性。
并行计算是指同时使用多种计算资源解决计算问题的过程，且主要目的是快速解决大型且复杂的计算问题。其基本思想是用多个处理器来协同求解同一问题，即大问题分解为小问题，然后通过多个处理器来解决小问题以达到同时解决的效果。许多问题都会运用到并行计算，比如计算流体动力学、天文物理、生物学、遗传学、气候模拟、水文模拟、医学等。随着各个领域需要大量数据的计算与处理，并行计算能够有效的提高处理效率。因此，将异构并行计算技术引入水文模型中进行参数率定，可以解决计算量大时间长的问题。
将传统的L-BFGS算法与异构并行计算技术结合使用时，L-BFGS在CPU上运行并无问题。但是由于该算法在搜索过程中仅能进行一条搜索路径的寻优，而与CPU搭配的协处理器拥有几百上千线程。此时，传统L-BFGS算法面对协处理器，依然使用单一搜索路径就显得非常浪费线程资源。因此，本申请对L-BFGS进行深入分析，挖掘其搜索方向及步长的选取原理，在保留算法总体搜索方式的基础上对其进行并行化改造，使算法具备并发处理多个搜索线路的能力，实现与协处理器的架构相适应。改进后的算法可以称为并发L-BFGS算法。
图3示出了本申请实施例中并发L-BFGS算法的原理示意图。如图3所示，该过程包括：
将多组参数分别代入算法的目标函数中，获得每组参数对应的目标函数值；
通过差商求导法获得每组参数对应的梯度值；
将每组参数及其对应的目标函数值和梯度值代入并发L-BFGS算法中进行循环求解；
通过判断是否满足终止条件，结束算法循环。
当不满足时，计算参数迭代前后的模型参数的差值和梯度值的差值，进而拟合海森矩阵；根据拟合的海森矩阵更新搜索方向和每组参数，进而更新Nash系数。
本申请设计的并发L-BFGS算法运行在CPU上，参数率定方案主要是综合GPU的并行计算，协同CPU工作解决水文模型参数率定中计算量大、精度不高、耗时长等问题。GPU是协处理器的一种。GPU刚开始只是在图形图像中得到应用，如今由于GPU的计算能力不断提高，加上计算结构的不断完善，使得GPU在并行计算中也得到了应用。GPU的显著特点是计算单元众多，这决定了其非常适合于大规模的计算，如矩阵计算。近些年，由于GPU的可编程性、计算能力不断提高的优势，其在人工智能、大数据方面的应用也十分成熟。
改进后的并发L-BFGS算法同样符合拟牛顿法的思想，不使用海森矩阵，而是构造一个近似海森矩阵(或其逆矩阵)的正定对称阵来代替，在“拟牛顿”的条件下优化目标函数。
下面，示例性的给出并发L-BFGS算法的海森矩阵(或逆矩阵)的推导过程。
并发L-BFGS算法的目标函数可表示为f(X)，设Xi,k为第i组迭代k次(k＝1,2，...，下同)后计算得到目标函数值fi(Xk)的极小值点的估计值，在Xi,k处进行泰勒展开，得到公式(1)。
令为第i组参数对应的梯度向量，/>为第i组参数对应的海森矩阵。记Bi,k≈Hi,k，/>yi,k＝gi,k+1-gi,k，si,k＝Xi,k+1-Xi,k。由拟牛顿条件，可得近似公式：/>
由上可知，获得并发L-BFGS算法的关键是拟合近似海森矩阵(或逆矩阵)的公式。并发L-BFGS算法的拟合海森矩阵的过程与DFP算法类似。
传统的DFP算法是采用海森逆矩阵D计算，但并不直接计算D，而是计算每一步D的增量来间接求D，因为一般上一步的中间结果对下一步的计算仍有价值，将i组同时执行，则有Di,k+1＝Di,k+ΔDi,k。通常将初始的海森逆矩阵取单位矩阵Ii，DFP算法计算的关键便是推导出ΔDi,k。通过一系列的推导假设，得出推导结果为：
因此，对并发执行的第i组计算，最终的递推关系为：其中：/>ρi,k和Vi,k分别为第一中间变量和第二中间变量。
下面介绍并发BFGS算法的执行流程。该过程包括如下步骤：
1)给出迭代终止条件，初始化Di,0＝Ii、i＝1、k＝1以及并发执行的组数m和参数Xi,k，i∈[1,m]；
2)k≠0时，计算并发执行的第i组迭代k次后的参数值Xi,k进行迭代搜索时的搜索方向di,k＝-Di,kgi,k；
3)根据并发BFGS算法的目标函数，按照公式(2)在实数范围内计算步长因子λi,k，然后按照公式(3)和(4)计算出si,k、Xi,k+1；
si,k＝λi,kdi,k (3)
Xi,k+1＝Xi,k+si,k (4)
4)判断是否满足终止条件，若满足则停止迭代，若不满足则执行下一步；
5)先计算Xi,k+1对应的梯度值，然后计算yi,k＝gi,k+1-gi,k；
6)根据si,k和yi,k，按照公式(5)拟合海森逆矩阵；
7)令k＝k+1，返回步骤2)。
上述公式(2)的具体做法是：
在实数范围内选择多组模型参数对应的多个步长因子λ；
根据预先获得的多组模型参数对应的搜索方向和多个步长因子λ，获得多组模型参数对应的多个候选参数Xi,k+λdi,k；
将多组模型参数对应的多个候选参数Xi,k+λdi,k代入更新后的目标函数，获得所述多组模型参数对应的多个候选参数的纳什效率系数f(Xi,k+λdi,k)；
将最小纳什效率系数对应的候选参数作为更新后的多组模型参数。
公知的，L-BFGS算法的执行过程与BFGS算法类似，但L-BFGS算法与BFGS算法存在一定区别。前者不再存储完整的Di,k，因为Di,k比较大，计算较为复杂。根据前面的推导可知，Di,k只与Di,0以及并发执行的第i组序列{si,k}和{yi,k}有关。即我们知道了后者，也可求得前者。对海森逆矩阵近似时，我们只需要序列{si,k}和{yi,k}中最近的z个值。这样，计算机内存中只需要存储这两个序列即可。
因此，L-BFGS算法在计算搜索方向di,k与BFGS算法存在不同，前者是根据当前迭代之前的Z次历史更新后的多组模型参数变化量及对应的梯度值变化量，获得当前迭代的搜索方向。Z次历史更新后的多组模型参数变化量根据Z次历史更新前的多组模型参数与Z次历史更新后的多组模型参数获得，Z次历史更新后的多组模型参数对应的梯度值变化量根据Z次历史更新前的多组模型参数对应的梯度值与Z次历史更新后的多组模型参数对应的梯度值获得。
也就是说，只保存当前迭代次数之前的Z次迭代获得的si,k和yi,k，k∈[k-z,k-1]，因此计算di,k时，对m组参数执行如下具体过程：
令qa,k＝ga,k，先执行后向循环：
for j＝k-1 to k-z do
gi,j＝gi,j+1-αi,jyi,j
end for
执行后向循环之后，可获得di,k-z-1＝Hi，0gi,k-z，然后执行前向循环获得搜索方向di,k：
forj＝k-z to k-1 do
di,j＝di,j-1+si,j(αi,j-βi,j)。
end for
上述后向循环，具体是在j取k-1到k-z范围内，循环执行和gi,j＝gi,j+1-αi,jyi,j，获得j＝k-1～k-z时，αi,j和gi,j的值，用于前向循环。
上述前向循环，具体是在j取k-z到k-1范围内，循环执行和di,j＝di,j-1+si,j(αi,j-βi,j)，获得第k次迭代是各组参数对应的搜索方向di,k。
上述循环中，α和β分别为第一变量和第二变量。
要想将上述并发L-BFGS算法结合异构并行计算技术，应用于水文模型参数率定中，还需要对水文模型进行异构移植。水文模型以HIMS(hydroinformatic modelingsystem)模型为例，HIMS模型最初是面向CPU端开发，并运行于PC机环境。HIMS模型主要应用在流域内评价、规划和管理水资源，以及处理与水相关的生态环境保护工作，如洪水预报、管理水污染与侵蚀以及气候变化等。HIMS具有多源信息融合、数字流域分析、分布式模拟以及模型定制等功能，可以针对国内各个流域不同的水文尺度、流域空间非均一性以及不同的自然和人文环境进行水循环模拟、参数率定等水文工作。目前，HIMS模型已应用并部署于国内外众多流域的水文模拟过程中，并取得了良好的应用效果。
图4示出了HIMS模型的工作流程图。如图4所示，HIMS模型对一组N维的数组进行马斯京根过程，输出模拟值；然后计算模拟值与存储的实测值的纳什效率系数。
当利用L-BFGS算法结合异构并行计算技术对HIMS模型进行多组参数率定时，需要对HIMS整体架构进行调整。根据目前的协处理器GPU的体系结构，原HIMS模型中用于读取流域数据而开辟的存储空间需要在计算前全部转移到协处理器，并在计算结束后将纳什系数作为结果传回主处理器。从CPU到GPU的移植过程包括：将涉及的函数定义成device，原始HIMS模型中所有的global类型变量和指针全部通过核函数传入协处理器，模型内部所有变量的内存大小需调整为原先的a倍。
图5示出了HIMS模型在异构平台的工作流程图。如图5所示，该过程包括：
主处理器(host端)生成由a组n维参数组成的数组Xi,k(x1,1,x1,2,x1,n,...,xm,1,xm,2,...,xm,n)，将生成的数组传入协处理器；
协处理器(device端)使用线程ID作为标识将多组参数划分至不同线程，各线程分别独取模型存储的流域数据(如降雨量、温度、海拔、流域面积等)及其对应的水文观测值，运行HIMS模型，并进行马斯京根过程(包括产流和汇流两大部分)，获得m组模拟值；然后，使用得到的m组模拟值及其对应的实测值代入纳什效率系数计算公式，得到m个Nash值，将Nash值取反后传回主处理器。
在将水文模型移植到协处理器后，下面结合附图示例性的说明并发L-BFGS应用于异构平台对水文模型进行参数率定的过程。
图6示出了本申请实施例提供的并发L-BFGS应用于异构平台对水文模型进行参数率定的流程图。异构平台包括主处理器和协处理器，下面以主处理器为CPU、协处理器为GPU为例，介绍本实施例的率定过程。
如图6所示，该过程包括：
CPU生成多组模型参数，并传入GPU中运行水文模型，获得每组模型参数对应的Nash系数，取反后传回CPU；
CPU通过差商求导法模拟获得每组参数中每个参数对应的梯度值，将每组参数及其对应的Nash系数取反得到的值和梯度值代入并发L-BFGS算法进行循环求解；
判断是否达到循环终止条件；不满足时，拟合海森矩阵，计算当前迭代的搜索方向，更新每组参数及其Nash系数。
在本申请的实施例中，传入GPU中运行水文模型，包括：
GPU将多组模型参数划分至不同的线程，各线程分别根据划分的模型参数运行水文模型对当前流域进行水文模拟，以及根据划分的模型参数对应的当前流域的水文模拟值，获得划分的模型参数对应的纳什效率系数。其中，运行水文模型对当前流域进行水文模拟包括：将当前流域的历史水文数据输入模型，获得划分的模型参数对应的当前流域的水文模拟值。
历史水文数据包括：历史时刻的降雨量和温度，以及流域的海拔和流域面积。
在本申请的实施例中，水文模型参数率定时，梯度无法直接通过求导算出，所以这里用到的梯度值都是通过用差商计算导数的方法模拟算出。
差商求导的原理是将初值代入目标函数中，算出初始的函数值，之后将每一维度的初值均加上差商步长，再代入目标函数中算出每个维度的初值加上步长后的函数值，与未加步长算出的初始函数值作差再除以步长，就可以得到每个维度的梯度值。具体步骤如下：
将第一组参数(x1,1,x1,2,...x1,n)代入目标函数中得到第一组的函数值f1(x1,1,x1,2,...x1,n)；
在第一组中，将x1,1加上差商步长q后代入目标函数，得到此时的目标函数值f1(x1,1+q,x1,2,...x1,n)；
则第一组中x1,1处的梯度值为
同理，那么第一组中x1,n处的梯度值为
则第m组xm,n处的梯度值为
在本申请的实施例中，上述目标函数为纳什效率系数的计算公式，即 为当前流域t时刻的观测值，/>为当前流域t时刻的模拟值，/>为当前流域的观测值的平均值，T为观测周期。/>通过协处理器中各线程运行水文模型获得。
在本申请的实施例中，采用差商法计算梯度值包括：
主处理获得叠加差商步长后的多组模型参数，并通过协处理器获得叠加差商步长后的多组模型参数对应的纳什效率系数；
接着，主处理根据多组模型参数对应的纳什效率系数、叠加所述差商步长后的多组模型参数对应的纳什效率系数以及差商步长，按照上述计算梯度值的公式，获得多组模型参数对应的梯度值。
在本申请的实施例中，计算搜索方向的过程包括：
根据当前更新前的多组模型参数对应的梯度值、以及获取的当前更新前的Z次历史更新后的多组模型参数的变化量和Z次历史更新后的多组模型参数对应的梯度值变化量，获得当前更新前的多组模型参数对应的搜索方向。
其中，Z次历史更新后的多组模型参数的变化量根据Z次历史更新前的多组模型参数与Z次历史更新后的多组模型参数获得，Z次历史更新后的多组模型参数对应的梯度值变化量根据Z次历史更新前的多组模型参数对应的梯度值与Z次历史更新后的多组模型参数对应的梯度值获得。
在本申请的实施例中，计算搜索方向的过程具体可以是如下方法：
根据当前更新前的多组模型参数对应的梯度值、Z次历史更新后的多组模型参数的变化量和Z次历史更新后的多组模型参数对应的梯度值变化量，拟合当前更新前的多组模型参数对应的海森逆矩阵；
根据多组模型参数对应的海森逆矩阵和梯度值，获得多组模型参数对应的搜索方向。
在本申请的实施例中，计算搜索方向的过程还可以是如下方法：
根据所述多组模型参数对应的梯度值、所述Z次历史更新后的多组模型参数的变化量和所述Z次历史更新后的多组模型参数对应的梯度值变化量，进行后向循环，获得第Z次历史更新前的多组模型参数对应的搜索方向；
根据所述Z次历史更新前的多组模型参数及其对应的梯度值、所述Z次历史更新后的多组模型参数的变化量和所述Z次历史更新后的多组模型参数对应的梯度值变化量，进行前向循环，获得当前更新前的多组模型参数对应的搜索方向。
本实施例的上述方法中的前向循环和后向循环，可通过本申请前述的前向循环和后向循环的算法实现。相比前述通过直接拟合海森矩阵获得搜索方向的方法，可进一步提高参数率定的效率。
在本申请的实施例中，在获得搜索方向后，更新多组模型参数的过程，如本申请公式(2)所描述的方法，此处不再描述。
与上述方法实施例对应的，本申请实施例还提供一种水文模型参数率定系统。如图7所示，该系统包括：
主处理器，用于生成水文模型的多组模型参数，并发送至所述协处理器；
协处理器，用于将所述多组模型参数划分至不同的线程，各线程分别根据划分的模型参数运行所述水文模型对当前流域进行水文模拟，以及根据所述划分的模型参数对应的当前流域的水文模拟值，获得所述划分的模型参数对应的纳什效率系数，并发送至所述主处理器；
主处理器，还用于判断所述多组模型参数对应的纳什效率系数是否满足预设条件，当满足时，所述主处理器输出所述更新后的多组模型参数，算法结束，当不满足时，所述主处理器根据所述多组模型参数对应的纳什效率系数，获得所述多组模型参数对应的梯度值，利用并发L-BFGS算法获得更新后的多组模型参数，并发送至所述协处理器。
本申请的系统实施例中，各处理器的具体执行步骤如本申请前述的方法实施例中的各步骤，此处不再叙述。
本申请的实施例中的方法步骤可以通过硬件的方式来实现，也可以由处理器执行软件指令的方式来实现。软件指令可以由相应的软件模块组成，软件模块可以被存放于随机存取存储器(random access memory，RAM)、闪存、只读存储器(read-only memory，ROM)、可编程只读存储器(programmable rom，PROM)、可擦除可编程只读存储器(erasable PROM，EPROM)、电可擦除可编程只读存储器(electrically EPROM，EEPROM)、寄存器、硬盘、移动硬盘、CD-ROM或者本领域熟知的任何其它形式的存储介质中。一种示例性的存储介质耦合至处理器，从而使处理器能够从该存储介质读取信息，且可向该存储介质写入信息。当然，存储介质也可以是处理器的组成部分。处理器和存储介质可以位于ASIC中。
在上述实施例中，可以全部或部分地通过软件、硬件、固件或者其任意组合来实现。当使用软件实现时，可以全部或部分地以计算机程序产品的形式实现。所述计算机程序产品包括一个或多个计算机指令。在计算机上加载和执行所述计算机程序指令时，全部或部分地产生按照本申请实施例所述的流程或功能。所述计算机可以是通用计算机、专用计算机、计算机网络、或者其他可编程装置。所述计算机指令可以存储在计算机可读存储介质中，或者通过所述计算机可读存储介质进行传输。所述计算机指令可以从一个网站站点、计算机、服务器或数据中心通过有线(例如同轴电缆、光纤、数字用户线(DSL))或无线(例如红外、无线、微波等)方式向另一个网站站点、计算机、服务器或数据中心进行传输。所述计算机可读存储介质可以是计算机能够存取的任何可用介质或者是包含一个或多个可用介质集成的服务器、数据中心等数据存储设备。所述可用介质可以是磁性介质，(例如，软盘、硬盘、磁带)、光介质(例如，DVD)、或者半导体介质(例如固态硬盘(solid state disk，SSD))等。
可以理解的是，在本申请的实施例中涉及的各种数字编号仅为描述方便进行的区分，并不用来限制本申请的实施例的范围。
Claims (8)
1.一种基于并发L-BFGS算法的水文模型参数率定方法，应用于异构平台，所述异构平台包括：主处理器和协处理器，其特征在于，所述方法包括：
步骤S1.所述主处理器生成水文模型的多组模型参数，并发送至所述协处理器；
步骤S2.所述协处理器将所述多组模型参数划分至不同的线程，各线程分别根据划分的模型参数运行所述水文模型对当前流域进行水文模拟，以及根据所述划分的模型参数对应的当前流域的水文模拟值，获得所述划分的模型参数对应的纳什效率系数，并发送至所述主处理器；
步骤S3.所述主处理器判断所述多组模型参数对应的纳什效率系数是否满足预设条件，当满足时，所述主处理器输出更新后的多组模型参数，模型率定结束，当不满足时，进入步骤S4；
步骤S4.所述主处理器根据所述多组模型参数对应的纳什效率系数，获得所述多组模型参数对应的梯度值，利用并发L-BFGS算法获得更新后的多组模型参数，并发送至所述协处理器，返回步骤S2；
其中，所述利用并发L-BFGS算法获得更新后的多组模型参数包括：
根据当前更新前的多组模型参数对应的梯度值、当前更新前的Z次历史更新后的多组模型参数的变化量和所述Z次历史更新后的多组模型参数对应的梯度值变化量，进行后向循环，获得第Z次历史更新前的多组模型参数对应的搜索方向；
根据所述第Z次历史更新前的多组模型参数对应的搜索方向、所述Z次历史更新后的多组模型参数的变化量和所述Z次历史更新后的多组模型参数对应的梯度值变化量，进行前向循环，获得当前更新前的多组模型参数对应的搜索方向；
根据所述当前更新前的多组模型参数对应的搜索方向和多个步长因子，获得多组候选的多组模型参数，并通过所述协处理器获得各组候选的多组模型参数对应的纳什效率系数；
根据所述各组候选的多组模型参数对应纳什效率系数，确定当前更新后的多组模型参数。
2.根据权利要求1所述的方法，其特征在于，所述主处理器根据所述多组模型参数对应的纳什效率系数，获得所述多组模型参数对应的梯度值包括：
所述主处理根据所述多组模型参数和预设的差商步长，获得叠加所述差商步长后的多组模型参数，并通过所述协处理器获得所述叠加所述差商步长后的多组模型参数对应的纳什效率系数；
所述主处理根据所述多组模型参数对应的纳什效率系数、所述叠加所述差商步长后的多组模型参数对应的纳什效率系数以及所述差商步长，获得所述多组模型参数对应的梯度值。
3.根据权利要求1所述的方法，其特征在于，所述Z次历史更新后的多组模型参数的变化量根据所述Z次历史更新前的多组模型参数与所述Z次历史更新后的多组模型参数获得，所述Z次历史更新后的多组模型参数对应的梯度值变化量根据所述Z次历史更新前的多组模型参数对应的梯度值与所述Z次历史更新后的多组模型参数对应的梯度值获得。
4.根据权利要求1所述的方法，其特征在于，所述方法还包括：
在实数范围内选择所述多组模型参数对应的多个步长因子。
5.一种实施如权利要求1-4任一项所述的基于并发L-BFGS算法的水文模型参数率定系统，其特征在于，所述系统包括：
主处理器，用于生成水文模型的多组模型参数，并发送至协处理器；
协处理器，用于将所述多组模型参数划分至不同的线程，各线程分别根据划分的模型参数运行所述水文模型对当前流域进行水文模拟，以及根据所述划分的模型参数对应的当前流域的水文模拟值，获得所述划分的模型参数对应的纳什效率系数，并发送至所述主处理器；
主处理器，还用于判断所述多组模型参数对应的纳什效率系数是否满足预设条件，当满足时，所述主处理器输出更新后的多组模型参数，算法结束，当不满足时，所述主处理器根据所述多组模型参数对应的纳什效率系数，获得所述多组模型参数对应的梯度值，利用并发L-BFGS算法获得更新后的多组模型参数，并发送至所述协处理器；
其中，所述主处理器利用并发L-BFGS算法获得每次更新后的多组模型参数包括：
根据当前更新前的多组模型参数对应的梯度值、当前更新前的Z次历史更新后的多组模型参数的变化量和所述Z次历史更新后的多组模型参数对应的梯度值变化量，进行后向循环，获得第Z次历史更新前的多组模型参数对应的搜索方向；
根据所述第Z次历史更新前的多组模型参数对应的搜索方向、所述Z次历史更新后的多组模型参数的变化量和所述Z次历史更新后的多组模型参数对应的梯度值变化量，进行前向循环，获得当前更新前的多组模型参数对应的搜索方向；
根据所述当前更新前的多组模型参数对应的搜索方向和多个步长因子，获得多组候选的多组模型参数，并通过所述协处理器获得各组候选的多组模型参数对应的纳什效率系数；
根据所述各组候选的多组模型参数对应纳什效率系数，确定当前更新后的多组模型参数。
6.根据权利要求5所述的系统，其特征在于，所述主处理器具体用于：
所述主处理根据所述多组模型参数和预设的差商步长，获得叠加所述差商步长后的多组模型参数，并通过所述协处理器获得所述叠加所述差商步长后的多组模型参数对应的纳什效率系数；
所述主处理根据所述多组模型参数对应的纳什效率系数、所述叠加所述差商步长后的多组模型参数对应的纳什效率系数以及所述差商步长，获得所述多组模型参数对应的梯度值。
7.根据权利要求5所述的系统，其特征在于，所述Z次历史更新后的多组模型参数的变化量根据所述Z次历史更新前的多组模型参数与所述Z次历史更新后的多组模型参数获得，所述Z次历史更新后的多组模型参数对应的梯度值变化量根据所述Z次历史更新前的多组模型参数对应的梯度值与所述Z次历史更新后的多组模型参数对应的梯度值获得。
8.根据权利要求5所述的系统，其特征在于，所述主处理器还用于：
在实数范围内选择所述多组模型参数对应的多个步长因子。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110056440.9A CN112711862B (zh) | 2021-01-15 | 2021-01-15 | 基于并发l-bfgs算法的水文模型参数率定方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110056440.9A CN112711862B (zh) | 2021-01-15 | 2021-01-15 | 基于并发l-bfgs算法的水文模型参数率定方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112711862A CN112711862A (zh) | 2021-04-27 |
| CN112711862B true CN112711862B (zh) | 2024-04-09 |
Family
ID=75549144
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110056440.9A Active CN112711862B (zh) | 2021-01-15 | 2021-01-15 | 基于并发l-bfgs算法的水文模型参数率定方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112711862B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113761807B (zh) * | 2021-09-28 | 2022-11-15 | 科瑞特空调集团有限公司 | 一种基于混合优化的快速空调风水系统优化方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5729451A (en) * | 1995-12-01 | 1998-03-17 | Coleman Research Corporation | Apparatus and method for fusing diverse data |
| CN106682355A (zh) * | 2017-01-12 | 2017-05-17 | 中国水利水电科学研究院 | 一种基于pso‑ga混合算法的水文模型参数率定方法 |
| CN111259522A (zh) * | 2020-01-09 | 2020-06-09 | 河海大学 | 一种水文模型在地理空间上多流域并行率定的方法 |
| CN112113545A (zh) * | 2020-09-17 | 2020-12-22 | 中国科学院海洋研究所 | 一种基于多维海面信息的内波振幅反演方法 |
-
2021
- 2021-01-15 CN CN202110056440.9A patent/CN112711862B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5729451A (en) * | 1995-12-01 | 1998-03-17 | Coleman Research Corporation | Apparatus and method for fusing diverse data |
| CN106682355A (zh) * | 2017-01-12 | 2017-05-17 | 中国水利水电科学研究院 | 一种基于pso‑ga混合算法的水文模型参数率定方法 |
| CN111259522A (zh) * | 2020-01-09 | 2020-06-09 | 河海大学 | 一种水文模型在地理空间上多流域并行率定的方法 |
| CN112113545A (zh) * | 2020-09-17 | 2020-12-22 | 中国科学院海洋研究所 | 一种基于多维海面信息的内波振幅反演方法 |
Non-Patent Citations (1)
| Title |
|---|
| 并发L-BFGS异构率定算法设计与实现;田在荣等;青岛大学学报(自然科学版);第34卷(第3期);43-50 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112711862A (zh) | 2021-04-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Liu et al. | Sequential optimization using multi-level cokriging and extended expected improvement criterion | |
| CN114254561A (zh) | 一种内涝预测方法、系统及存储介质 | |
| CN110264471A (zh) | 一种图像分割方法、装置、存储介质及终端设备 | |
| Agoshkov et al. | Problems of variational assimilation of observational data for ocean general circulation models and methods for their solution | |
| Prayogo | Metaheuristic-based machine learning system for prediction of compressive strength based on concrete mixture properties and early-age strength test results | |
| CN114970302B (zh) | 一种基于地下水监测系统的区域地下水情预测方法 | |
| Monteil et al. | Multi-objective calibration by combination of stochastic and gradient-like parameter generation rules–the caRamel algorithm | |
| CN112711862B (zh) | 基于并发l-bfgs算法的水文模型参数率定方法及系统 | |
| CN112733997A (zh) | 基于woa-lstm-mc的水文时间序列预测优化方法 | |
| Xia et al. | Efficient parallel surrogate optimization algorithm and framework with application to parameter calibration of computationally expensive three-dimensional hydrodynamic lake PDE models | |
| Yuan et al. | Modelling and pathway identification involving the transport mechanism of a complex metabolic system in batch culture | |
| Cuomo et al. | Toward a multi-level parallel framework on GPU cluster with PetSC-CUDA for PDE-based Optical Flow computation | |
| CN112131794A (zh) | 基于lstm网络的水工建筑物多效应量优化预测及可视化方法 | |
| Hu et al. | Dynamically Optimized Unstructured Grid (DOUG) for Analog Ensemble of numerical weather predictions using evolutionary algorithms | |
| CN116629352A (zh) | 一种亿级参数寻优平台 | |
| CN115965160A (zh) | 一种数据中心能耗预测方法、装置、存储介质及电子设备 | |
| He et al. | Research on geometric features and point cloud properties for tree skeleton extraction | |
| Hu et al. | ℓ-DARTS: Light-weight differentiable architecture search with robustness enhancement strategy | |
| Obiols-Sales et al. | NUNet: Deep learning for non-uniform super-resolution of turbulent flows | |
| Xing et al. | Accelerating reliability-based topology optimization via gradient online learning and prediction | |
| Leng et al. | Variable-fidelity surrogate model based on transfer learning and its application in multidisciplinary design optimization of aircraft | |
| Han et al. | A kriging-based active learning algorithm for contour estimation of integrated response with noise factors | |
| CN116110630A (zh) | Dcs系统中反应堆堆芯热功率测量方法与装置 | |
| Peng et al. | Enhanced Adaptive Graph Convolutional Network for Long-term Fine-grained SST Prediction | |
| CN117436160A (zh) | 面向大坝变形的多点监测互验拟合方法、装置及存储介质 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |