趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

河道内貯留を考える:Dynamic wave VS. Kinematic wave

河道内貯留のついて勉強のためにも一度計算してみました。

河道内貯留とは

河道内貯留とは、下図1のように下流への伝播に伴いピーク流量(水深)が減衰する現象を示す。

f:id:SedimentHydraulics:20200919224847p:plain

このような伝播を計算するためには洪水波の移流拡散を考慮する必要があるため、次式のDynamic wave modelを用いる必要がある。


\begin{align}
   & \frac{\partial Q}{\partial t} + \frac{\partial }{\partial x}\left(\frac{Q^2}{A}\right)+gA\frac{\partial H}{\partial x}+gAi_e = 0 \\
   & \dfrac{\partial A}{\partial t}+\dfrac{\partial Q}{\partial x} = 0 
\end{align}

また、次の移流拡散項を省略したKinematic wave modelではこの伝播特性を計算することはできない。


\begin{align}
   & Q = \dfrac{1}{n} R^{2/3} i_b^{1/2}A \\
   & \dfrac{\partial A}{\partial t}+\dfrac{\partial Q}{\partial x} = 0 
\end{align}

数値計算:Dynamic waveとKinematic waveの比較

計算条件

  • 河道断面:矩形断面
  • 河床勾配1/1000
  • 川幅100m
  • 縦断距離20km
  • マニングの粗度係数0.03
  • 上流端流量は実洪水波形が良いと思い、2019年度台風19号相模川の城山ダムの放流量を参照しました。

f:id:SedimentHydraulics:20200919180409p:plain

出典 https://www.mlit.go.jp/river/shinngikai_blog/damchousetsu_kentoukai/dai01kai/2-3_R1T19_dam_taiou.pdf

元データがなかったのでトレースしました。

計算結果

計算結果より上流から0km、10km、20km地点の流量の時系列変化を示す。

見易いようにピーク付近の拡大図を示す。Dynamic waveではピーク流量が200m3/s程度減衰していることが確認できる。洪水全体でみるとわずかにも見えるが、200m3/sはダム1個分くらいの容量である。

f:id:SedimentHydraulics:20200920002330p:plain

水理学の教科書等では河床勾配1/3000より急な河川では河道内貯留はそれほど重要でないという記載がしばしば見られますが、私の経験上、実河川の流量波形は大洪水時にはかなり急になるので急流河川でも十分に河道内貯留は生じると考えています。

余談ですが、

日本の治水計画では河道内貯留は考慮していない。これは計画論なので仕方無いこともあるが、実洪水の再現計算にKinematic wave(貯留関数法の河道モデルも同様)を用いることは現象把握の観点からは適切でないと考える。

Github

Jupyter Notebook Viewer


コードもぜひ。一般断面(自然河道断面)のkinematic wave法の勾配はどう設定すればいいのかな。。。


  1. 椿:水理学演習 下

    水理学演習 下 POD版

    水理学演習 下 POD版

Manningの粗度係数を再考する

頭の整理に書いてみました。ご指摘があれば


Manningの公式の概要

Manningの公式は、Robert Manningによって提案された横断面内で平均化された流速Vを示す実験公式であり、次式のように示される。


\begin{align}
 V = \dfrac{1}{n}i_e^{1/2}R^{2/3}
\end{align}

ここで、n:Manningの粗度係数(m^{-1/3} \cdot s)、i_e:エネルギー勾配、R:径深(m)である。

本式は現在世界で最も使われている平均流速公式である。 上式の実験定数はManningの粗度係数のみであるが次元を有しており定数として適切でない。 それにもかかわらず、式形が非常に単純なため、実験水路から実河川まで様々な条件下での豊富な観測資料を元にManningの粗度係数は整理されている。そのため、容易に概略値を把握することが可能である。

例えば、Chowの「OPEN-CHANNEL HYDRAULICS」では水路条件を100種類以上に分類し、Manningの粗度係数を一般値を示している。1

日本では「水理公式集(平成11年版)」2や「河川砂防技術基準(案)同解説調査編」3に次表のようにまとめられている。

  • 「水理公式集(平成11年版)」
  • 「河川砂防技術基準(案)同解説調査編」

Manningの公式への水理学的な意味付け

流速係数\phi=V/u_*をManningの公式を用いて整理すると次式となる。


\begin{align}
\phi &= \frac{V}{u_*} = \frac{k_s^{1/6}}{n \sqrt{g}} \left(\frac{R}{k_s} \right)^{1/6} \nonumber 
\end{align}

ここにk_s:相当粗度(m)である。 一方開水路完全粗面の対数則による流速係数は次式となる。


\begin{align}
\phi &= \frac{1}{\kappa}\log_e{\frac{h}{k}} - \frac{1}{\kappa} + C  
\end{align}

両式より、\phik_s^{1/6}/(n \sqrt{g})の関係を整理したものが次図の点線である。

図 Manning・Stricklerの式と粗面対数則との比較4

管路や開水路では流速係数\phiは8~25程度であるため、k_s^{1/6}/(n \sqrt{g})は図中直線の一定値に近似できる。 これより、流速係数は以下のように近似できる。


\begin{align}
\frac{k_s^{1/6}}{n \sqrt{g}} & \sim 7.66 \nonumber  \\
\phi &= \frac{V}{u_*} = 7.66 \left(\frac{R}{k_s} \right)^{1/6} 
\end{align}

これはManning・Stricklerの公式と呼ばれる。 上記のように近似できることは実用上は径深によってManningの粗度係数が変わらないことを示している。 今日では本式を使用することは少ないが、Manningの公式に水理学的な意味を与えた点で重要である。

沖積河川の河床抵抗

多くの沖積河川では河床上に河床波(小規模河床形態)が形成されている。河床波とは砂漣、砂堆、反砂堆等の水深スケールの河床形状を示して、日本土木学会水理委員会によると次表のように分類される。5 また、Simonsはより詳細な分類を示している。6

河床波(ここでは砂堆)上の流れを模式的に示すと下図のとおりである。7

ここで、(水深平均)流れに対する河床が与える抵抗を考えると「摩擦抵抗」と「剥離に伴う形状抵抗」に分類できる。 Manningの公式のような平均流速公式ではこの両者の影響が混ざって含まれることになる。

「摩擦抵抗」は河床材料によって決定するため対数則によって評価する。一方、「剥離に伴う形状抵抗」は河床波の形状を考慮して急拡損失の考え方により評価する。 ここで河床波の形状が問題となる。模型実験では二次元水路を使用することにより河床波の形状を計測することは比較的容易であるが、実河川では難しく、近年の研究で高解像度のセンサーを用いた実河川の河床波の計測事例が報告されている8が、まだ一般的とは言えず、実河川の河床波形状を計測すること難しい。

さらに、実河川の流れは非定常にであり、流れの規模によって河床波の形状が変化する。このような実河川における出水中の河床波の変化を計測事例はほとんど報告されておらず9,10、実態は不明である。

これらより実河川の河床抵抗を厳密に推定することは難しいことが理解できる。そのため、Manningの公式が積極的に代用されている。前項にはManningの粗度係数は水深によってほとんど変化しないと示したが、これは河床波形状が変化しない条件で成立し、変化する場合には形状抵抗による損失が変化するため、Manningの粗度係数も変化する。

なお、Manningの粗度係数は一定値として取り扱うことが一般的である。観測資料が十分でなく、大規模出水時に観測データから同定された値を使うためである。

河床波の変化を考慮した河床抵抗の評価方法

前項のとおり、実河川での検証が難しいため、模型実験結果を元に河床抵抗の評価を行う。

Engelund11は河床抵抗を以下のようにモデル化した。


\begin{align}
 \tau_* = \tau_*^{\prime} + \tau_*^{\prime\prime}
\end{align}

ここで、 \tau_*:全無次元掃流力、\tau_*^{\prime}:摩擦抵抗による無次元掃流力(=無次元有効掃流力)、 \tau_*^{\prime\prime}:形状抵抗による無次元掃流力である。

相似則より\tau_*\tau_*^{\prime}のみの関数であることを示した。これより実験結果を使用して次図12の関係式を示した。

さらに、岸・黒木[^5]はEngelundの相似則を修正し、\tau_*\tau_*^{\prime},R/dの関数となることを示し、次図[^12]の関係を示した。

現時点ではこれが最も有用な研究成果と考えられる。

参考までにEngelundと岸・黒木の研究成果を比較する。両者による無次元掃流力-有効無次元掃流力の図を重ね合わせると次図のとおりである。なお、岸・黒木の図はR/dが1000、500、100の3ケースを示す。

両者を比較するとdune領域、flat-bed・antidune領域ではそれほど大きな差異はないが、transition領域の取り扱いが異なっている。 Engelundはtransition領域の変化が不明なため定式化していないが、岸・黒木は定式化しており、この領域の取り扱いがR/dによって異なっている。 よって、transition領域の影響を考慮するためには岸・黒木の方法を使う他はない。

河床波の変化による流れ場への影響を理解しやすくするため、岸・黒木の図を\phi\tau_*の関係に変換した図[^12]を以下に示す。

縦軸の\phi=V/u_*は流れやすさを示すため、transition領域で急激に流れやすくなる、つまり、河床抵抗が減少することが理解できる。

岸・黒木の方法による河床抵抗の評価

岸・黒木の方法による河床抵抗の評価の計算方法について示す。 次図(再記)の各領域区分の無次元掃流力\tau_*と有効無次元掃流力\tau_*^{\prime}の関係式は以下のとおりとなっている。


\begin{align}
\rm{dune1} : \tau_*^{\prime} &= 0.21\tau_*^{1/2} \\
\rm{dune2} : \tau_*^{\prime} &= 1.49(R/d)^{-1/4}\tau_* \\
\rm{flat\text{-}bed} : \tau_*^{\prime} &=  \tau_* \\
\rm{antidune} : \tau_*^{\prime} &= 0.264(R/d)^{1/5}\tau_*^{1/2} \\
\rm{transition1} : \tau_*^{\prime} &= 6.5 \times 10^7(R/d)^{-5/2} \tau_*^{11/2}
\end{align}

また、dune領域とtransition領域、flat-bed領域とantidune領域の区分は次式となる。


\begin{align}
\rm{dune \text{ -- } transition} : \tau_* &= 0.02(R/d)^{1/2} \\
\rm{flat\text{-}bed \text{ -- } antidune} : \tau_* &= 0.07(R/d)^{2/5} 
\end{align}

上図中のtransition2領域は関係式が設定されていない。transition2はdune2とantiduneを直線で繋いでおり、dune領域の上限値を示している。 transition領域では複雑な変化を示すため定式化は難しいが、岸・黒木の方法より、dune1 -- transition1 -- flat-bed -- antiduneの過程で遷移する(図中赤線)ものと考える。 なお、antiduneは射流場のみで発生するため、ほとんど沖積河川の場合、flat-bedまでの変化を考慮すれば良い。

岸・黒木の方法の課題

河床波の遷移の履歴

ここまでに示した考え方は基本的に定常、等流、平衡を前提に構築されている。 また、模型実験のデータを基に関数をフィッティングしている。

実河川の適用性は、Engelund[^11]はRio Grande川での観測データと一致することを示しており、岸・黒木[^5]も同じデータを用いて妥当性を示している。 しかし、岸・黒木は石狩川の観測データとの比較により水理量と河床形状の関係が常に平衡にはならずに履歴の影響を受けること示しており、手法の適用上の課題を指摘している。

今後はより多くの河川で観測を実施して手法の適用性を検証する必要がある。

河床波の形状抵抗の直接評価

Engelundおよび岸・黒木の方法は、相似を仮定して無次元掃流力と有効無次元掃流力の関係を用いて河床波の影響を評価しているが、河床波の形状損失を直接評価する方が理想的である。

この方法は、これまでの研究では一次元流れを仮定した急拡損失による評価が行われてきたが、乱流モデルを用いた鉛直二次元または三次元解析によって評価する必要がある。 また、現時点では実河川による検証データの取得が難しいため、模型実験を基本に検証を進めることが望ましい。

参考:流砂量評価における有効掃流力

掃流力は流れにとっては抵抗になるが河床土砂にとっては駆動力となる。その力は全抵抗から形状抵抗を差し引いた有効掃流力である。 そのため、河床波の遷移は掃流砂量にも大きな影響を与える。

日本で最も使用される芦田・道上の掃流砂量式には有効掃流力が含まれている。これはdune上の流れを想定していることを意味しており、本来であれば、河床波形状に応じてこの式形を変更する必要がある。

また、使用頻度は低いが、佐藤・吉川・芦田の式(土研公式)ではManningの粗度係数によって式形を変更している。これは、河床波形状の影響を示しているが、河床波形状はManningの粗度係数のみでは決定しないため、使用の際には注意が必要である。

Github

Jupyter Notebook Viewer


  1. Chow,V.,T.: Open-channel hydraulics, pp.110-113, McGraw-Hill, 1959.

  2. 土木学会水理委員会水理公式集改訂小委員会 : 水理公式集, p.89, 土木学会, 1999.

  3. 建設省河川局,日本河川協会 : 建設省河川砂防技術基準(案)同解説調査編, p.133, 技報堂出版, 1997.

  4. 椿東一郎 : 水理学I, p.110, 森北出版, 1973.

  5. 水理委員会移動床流れの抵抗と河床形状研究小委員会 : 移動床流れにおける河床形態と粗度, 土木学会論文報告集第210号, pp.65-91, 1973.

  6. Simons,D.B.,Şentürk, F. : Sediment Transport Technology, Water Resources Publications, 1977.

  7. 芦田和男,江頭進治,中川一 : 21世紀の河川学:安全で自然豊かな河川を目指して, p.121, 京都大学学術出版会, 2008.

  8. 例えば、秋田麗子,西口亮太,野間口芳希 : 水中の河床地形の面的計測とその活用方策について, 河川技術論文集, 第23巻, pp.173-178, 2017.

  9. 高部一彦,人見寿,坂野章,山本浩一 : 涸沼川における河床波観測及び解析, 河川技術論文集, 第11巻, pp.387-392, 2005.

  10. 末次忠司, 日下部隆昭, 坊野聡子 : 土砂管理施策のためのキーノート, 国土技術政策総合研究所資料, No.231, 2005.

  11. Engelund F. : Hydraulic Resistance of Alluvial Streams, Journal of the Hydraulics Division, Vol. 92, Issue 2, pp.315-326, 1966.

  12. 河村三郎 : 土砂水理学, pp.226-230, 森北出版, 1982.

Engelundの有効掃流力の図をトレースしてみる

有名なEngelundによる有効掃流力の上図(元図はEngelund、上図は河村のp.226、実験値はUSGSのレポートらしい)は反砂堆領域の曲線の式が設定されてないです。

何かと不便なのでトレースしてみました。

こんな感じです。

トレースしたデータはこちら

この図は実験値と異常なくらいフィットしているなと思っていたのですが、Engelund and Hansenの関連論文を読んでいると、

To avoid confusion the points corresponding to the transition between upper and lower flow regime are omitted.

と書いてました。 この手の実験だとtransition領域の結果はぐちゃぐちゃでフィットしないのでそれを載せたほうが良い気がします。

Github

Jupyter Notebook Viewer


トレースのツール

図のトレースにはこちらを使いました。

WebPlotDigitizer 4.3

https://apps.automeris.io/wpd/

日本語の解説記事はこちら

ダウンロード不要!論文のグラフデータを数値化する技 – てつたま~哲学のたまご~


  • グラフの作成はHoloviewsを使用

computational-sediment-hyd.hatenablog.jp

  • はてブロへのhtmlグラフの埋め込みは以下を参照

computational-sediment-hyd.hatenablog.jp

備忘録:MinGWのmakeはmingw32-make

MinGWのmakeはmingw32-make

MinGWのmakeコマンドについて - (仮)

最近fortranを書く機会が減りましたが、年に数回書く際にmakeコマンドを忘れてしまいます。 windowsはやっぱり特殊ですね。個人的にはWSLを標準にして欲しいのですがなかなか完全には浸透しないですね。

備忘録:開水路において水深平均化された対数則

なんとなく気分転換に。この手の話って最近の本だとさらっとしか書いてないですよね。


開水路において水深平均化された対数則の導出

まず、完全粗面の円管流れの対数則による流速分布を示すと次のとおりとなる。


\begin{align}
\frac{u}{u_*}  = \frac{1}{\kappa}\log_e{\frac{y}{k}} + C = \frac{1}{\kappa \log_{10}{e}}\log_{10}{\frac{y}{k}}+ C 
\end{align}

ここに、\kappa:カルマン定数、k:粗度高さ、y:壁面からの距離である。 また、CはNikuradseの実験より8.5である。

次に水深流速Vに対する対数則を導出する。上式を水深積分により平均化すると次のとおりである。


\begin{align}
\frac{V}{u_*} &= \frac{1}{h}\int_0^h \frac{u}{u_*} dy \nonumber \\
 &= \frac{1}{\kappa}\log_e{\frac{h}{k}} - \frac{1}{\kappa} + C  \\
 &= \frac{1}{\kappa \log_{10}{e}}\log_{10}{\frac{h}{k}}- \frac{1}{\kappa}+ C 
\end{align}

これが開水路において水深平均化された対数則である。

取り扱いやすいように次式のように変形することが多い。


\begin{align}
\frac{V}{u_*} &= \frac{1}{\kappa}\log_e{\frac{h}{k}} - \frac{1}{\kappa} + C  = \frac{1}{\kappa}\log_e \left({e^{\kappa C-1}\frac{h}{k}} \right)
\end{align}

また、自然河道を対象とする場合水深hを径深Rとする。

manning-stricklerによる近似形

manning則を用いて式変形を行うと流速係数\phi = V/u_*は以下のように示される。 本式は自然河道の場合\phiが8~25であることを条件とした近似式である。


\begin{align}
\phi &= \frac{V}{u_*} = \frac{k_s^{1/6}}{n \sqrt{g}} \left(\frac{R}{k_s} \right)^{1/6} \nonumber \\
\frac{k_s^{1/6}}{n \sqrt{g}} & \sim 7.66
% \frac{1}{\kappa}\log_e{\frac{h}{k}} - \frac{1}{\kappa} + C  = \frac{1}{\kappa}\log_e \left({e^{\kappa C-1}\frac{h}{k}} \right)
\end{align}

まとめというか愚痴を

教科書等では開水路の対数則は以下のような式形が多く見られる。


\begin{align}
\frac{V}{u_*}  =  5.75\log_{10}{\frac{h}{k}}+ 6.0
\end{align}

この書き方をすると定数がたくさんあるように見えるのでやめて欲しいです。前項の書き方のほうが実験定数はNikuradseの8.5であることがわかりやすいと思う。

参考文献

水理学 1 (基礎土木工学全書6)

水理学 1 (基礎土木工学全書6)

  • 作者:椿東一郎
  • 発売日: 1973/02/01
  • メディア: 単行本(ソフトカバー)

Github

Jupyter Notebook Viewer

L-BFGS法(記憶制限BFGS法)の実装

ちょこちょこ使用していますが半年使わないと忘れてしまうので備忘録的にきっちりまとめておきます。


はじめに

非線形問題の最適化に最もよく使用されるL-BFGS法(記憶制限BFGS法、BFGS法は準ニュートン法の解法の一つ)について書きます。準ニュートン法なので勾配法ですね。

ディープラーニングでは確率的勾配降下法(力技!!)が使用されますが、最適化問題ディープラーニング以外の機械学習でもL-BFGS法はよく使われています。

しかし、L-BFGS法のオープンライブラリはあまり多くないです。 その理由としては、私が思うのものは、

  • 汎用的なものが書きづらい。用途に合わせたチューニングが必要。
  • コーディングが泥臭い。不安定なモデルのためエラー処理が山程必要。

です。

そんな感じで、本記事では汎用的なものでなく、L-BFGS法の主たる部分であるヘッセ行列の更新の部分を主に書いておきます。

L-BFGS法の導出

ググると山程引っかるのでここでは省略します。いくつか有用な記事をリンクしておきます。 私は上の2つを参照しました。

L-BFGS法はだからメモリが節約できるのか! - あらびき日記

Limited-memory BFGS - Wikipedia

L-BFGS法の更新式を導出 - ysk24ok.github.io

3分でわかるL-BFGS - Kotaro's blog

ちょっと今回の話とは異なりますが、L-BFGS法の高速化の論文も面白いので貼っておきます。

https://papers.nips.cc/paper/5333-large-scale-l-bfgs-using-mapreduce.pdf

L-BFGS法の実装

上記の数式を元にそのままコーディングします。いつもどおりpython使います。

import numpy as np

def lbfgs(x, f, g, stepsize, maxiterate, memorysize, epsilon):

    def searchdirection(s, y, g):
        q = -np.copy(g)
        num = len(s)
        a = np.zeros(num)
        
        if num==0 : return q
        
        for i in np.arange(num)[::-1]:
            a[i] = np.dot(s[i], q)/np.dot(y[i], s[i])
            q -= a[i]*y[i]
        
        q = np.dot(s[-1], y[-1]) / np.dot(y[-1], y[-1]) * q
            
        for i in range(num) : 
            b = np.dot(y[i], q) / np.dot(y[i], s[i])
            q += ( a[i] - b ) * s[i]
            
        return q
    
    outx = []
    s = np.empty((0, len(x)))
    y = np.empty((0, len(x)))
    xold = x.copy() 
    gold = g(xold) 
    J1 = f(xold)
    mmax = memorysize
    sp = stepsize
    print( 'f= ', J1 )
    
    outx.append(xold)
    for num in range(maxiterate):
        if np.linalg.norm(gold) < epsilon : 
            print('g=',np.linalg.norm(gold))
            break
            
        d = searchdirection(s, y, gold)
        
        sp = stepsize*0.01 if num == 0 else stepsize           
        
        xnew = xold + sp * d
        gnew = g(xnew) 
        J2 = f(xnew)
        
        J1 = J2
        si, yi = xnew - xold, gnew - gold
        if len(s) == mmax :
            s = np.roll(s, -1, axis=0)
            y = np.roll(y, -1, axis=0)
            s[-1] = si
            y[-1] = yi
        else:
            s = np.append(s, [si], axis=0)
            y = np.append(y, [yi], axis=0)
        
        xold, gold = xnew, gnew
        
        print( 'iterate= ', num, 'f= ', J1, ' stepsize= ', sp )
        outx.append(xold)
        
    return xold, outx

引数は、最適化する変数xと評価関数fとその勾配g、それに加えて、ステップ幅やイテレーションの最大値や収束条件、ヘッセ行列更新時に過去何回分の履歴を用いるかのパラメータmemorysizeを与えます。

コーディングの難しいところは特に無いですが、面白いところは配列の大きさがmemorysizeを超えた時点でnumpyのroll関数で配列の順番を更新しているところです。

私はイテレーションの1回目はステップ幅を十分に小さく取ります。今回の計算例の場合は大きくても問題ないですが、評価関数が微分方程式の場合は発散する可能性があるためです。

計算例

次式を例にL-BFGS法を適用してみます。

f:id:SedimentHydraulics:20200815185723p:plain

Xの初期値はx=-2.5,y=0で、最適値はx=1,y=0です。パラメータは適当です。ステップ幅は1です。

%%time
f = lambda x : 5.0*x[0]**2.0 - 6.0*x[0]*x[1] + 5.0*x[1]**2 - 10.0*x[0] + 6.0*x[1]
g = lambda x : np.array( [ 10.0 * x[0] - 6.0 *x[1] -10.0, 10.0 * x[1] - 6.0 *x[0] + 6.0 ] )

# initial value
x = np.array( [-2.5, 0.0] )

xx, out = lbfgs(x, f, g, stepsize=1.0, maxiterate=20, memorysize=5, epsilon=10.0**(-8))
print(xx)
> f=  56.25
> iterate=  0 f=  40.864  stepsize=  0.01
> iterate=  1 f=  -1.3307052922247742  stepsize=  1.0
> iterate=  2 f=  -3.1273465839434973  stepsize=  1.0
> iterate=  3 f=  -4.999999466916309  stepsize=  1.0
> iterate=  4 f=  -4.999999999941274  stepsize=  1.0
> iterate=  5 f=  -4.999999999999999  stepsize=  1.0
> iterate=  6 f=  -5.0  stepsize=  1.0
> g= 8.580967415593408e-10
> [1.00000000e+00 1.53469307e-10]
> Wall time: 4 ms

イテレーションは6回で収束します。その過程を図化するとこんな感じです。

scipy.optimize.fmin_l_bfgs_bとの比較

scipyには、L-BFGS法scipy.optimize.fmin_l_bfgs_bが、用意されています。一応それと比較しておきます。

%%time
f = lambda x : 5.0*x[0]**2.0 - 6.0*x[0]*x[1] + 5.0*x[1]**2 - 10.0*x[0] + 6.0*x[1]
g = lambda x : np.array( [ 10.0 * x[0] - 6.0 *x[1] -10.0, 10.0 * x[1] - 6.0 *x[0] + 6.0 ] )

out2 = []
def callback(xk):
    out2.append(xk)
    print(xk)
    
# initial value
x = np.array( [-2.5, 0.0] )
out2.append(x)
xx, y, d = opt.fmin_l_bfgs_b(f, x, fprime=g, m=5, iprint=0, epsilon=10.0**(-8), maxiter=20, maxls=1, callback=callback)
> [-1.64250707 -0.51449576]
> [ 0.10902424 -1.00977252]
> [ 0.4446687  -0.73580975]
> [ 1.00115755 -0.00273162]
> [ 1.00005759e+00 -1.68279906e-05]
> [1.00000095e+00 7.64484423e-07]
> Wall time: 3 ms

計算時間、イテレーションの回数、収束の過程はほとんど一緒ですね。

メモ

scipy.optimize.fmin_l_bfgs_b の中身は、次のライブラリで数千行に及ぶごりごりのFortran90です。(かつてこれを使ってFortran03のコードを書きました。難解すぎて死にそうでした。)

L-BFGS-B Nonlinear Optimization Code

もちろん高機能で幅広いケースに対応してますが、単純な問題なら自作の数十行のコードで十分です。

CG法との比較

単純な最急降下法との比較も示しておきます。ステップ幅を0.05、0.1と2ケースを書きました。 2ケースともにイテレーションは20回まで計算しました。

%%time
f = lambda x : 5.0*x[0]**2.0 - 6.0*x[0]*x[1] + 5.0*x[1]**2 - 10.0*x[0] + 6.0*x[1]
g = lambda x : np.array( [ 10.0 * x[0] - 6.0 *x[1] -10.0, 10.0 * x[1] - 6.0 *x[0] + 6.0 ] )

# initial value
x = np.array( [-2.5, 0.0] )

out3 = []
out3.append(x)
for _ in range(20):
    xx = x - 0.05*g(x)
    out3.append(xx)
    x = xx.copy()
    
# initial value
x = np.array( [-2.5, 0.0] )

out4 = []
out4.append(x)
for _ in range(20):
    xx = x - 0.1*g(x)
    out4.append(xx)
    x = xx.copy()

2ケースともにいいところまで収束しますが、その経路は大きく異なります。 このように最急降下法の場合はステップ幅が重要になります。

考察とまとめ

  • 単純なL-BFGS法の実装について書きました。実問題に通用するコードを書こうとするとはもうちょっと複雑になりますが、基本はわかるかと思います。
  • L-BFGS法の利点の一つは勾配降下法の重要なステップ幅の設定が適当で良い点である。これは非常に重要な点で、汎用コードを見てるとコードの大半はステップ幅の設定に関する部分(ラインサーチ)になるので、これを省略できることは有利になります。
  • また、ラインサーチではArmijo条件とWolf条件をよく使いますが、その際にfやgの計算を複数回実施します。単純な例では問題ないですが、流体問題の最適化ではfやgの計算1回に数十分~数時間かかることもよくあります。そのため、L-BFGS法は計算時間短縮の点からもメリットがあります。
  • 今後は、流体問題のL-BFGS法の適用例(ラグランジュの未定乗数法かな)も書いてみようと思います。

Github

github.com


  • グラフの作成はHoloviewsを使用

computational-sediment-hyd.hatenablog.jp

  • はてブロへのhtmlグラフの埋め込みは以下を参照

computational-sediment-hyd.hatenablog.jp

河川流一次元計算における河川横断面形状の簡略化について

モチベーション

  • 横断面平均一次元河川流計算で実河川を対象とする場合、地形条件は河川横断測量データを使用することが一般的である。
  • このデータより、水位に対する河積、川幅、潤辺等を計算するが計算負荷が大きく、特に長期計算を行う場合には計算速度に大きく影響する。
  • この改善方法として、河川横断測量データの簡略化とその影響について検証する。

河川横断データの簡素化

  • 基本の横断図は以下の過去記事のものを使用する。

computational-sediment-hyd.hatenablog.jp

  • 横断面の簡略化には、pythonのモジュールshapelysimplifyメソッドを使います。(手法の詳細は勉強中のため後日まとめます。)
  • 簡略化のパラメータはtoleranceのみ。
  • 簡略化した横断図と測点数は次のとおりとなる。
  • スケールによるが、見た目だけで判断するとtolerance = 1.0および2.0は粗く感じる。tolerance = 0.2や0.5はoriginalとほとんど変わらないように思う。
  • 測点数は、originalの80に対して、tolerance = 0.2でも24と大幅に減少している。
number of points
original 80
tolerance = 0.2 24
tolerance = 0.5 19
tolerance = 1.0 12
tolerance = 2.0 7

簡略化による物理量への影響

  • 河川流の場合は、水位Hと流量Qの関係が重要となるのでこれへの影響を確認する。
  • 等流を仮定して以下のマニングの式により評価する。


Q = \frac{1}{n}i_e^{1/2}R^{2/3}A

  • 右辺の径深R、河積Aは水位Hの関数になるので、H \sim R^{2/3}Aを比較する。
  • Hが17~24の範囲で図化すると次のとおりとなる。
  • tolerance = 2.0でoriginalとの誤差は10%と相対的に大きいが、それ以外では局所を除くと誤差は非常に小さい。

簡略化による計算時間への影響

computational time
original 686 µs ± 11.5 µs
tolerance = 0.2 288 µs ± 5.87 µs
tolerance = 0.5 265 µs ± 12.1 µs
tolerance = 1.0 207 µs ± 7.21 µs
tolerance = 2.0 174 µs ± 2.27 µs

まとめ

  • tolerance = 0.2の場合でも計算時間は十分に高速し、計算精度は保証されるため、河川流計算の前処理として簡略化を行ったほうが合理的である。

  • グラフの作成はHoloviewsを使用

computational-sediment-hyd.hatenablog.jp

  • はてブロへのhtmlグラフの埋め込みは以下を参照

computational-sediment-hyd.hatenablog.jp