2020-08-09

【最適化】SCE-UA法

SCE-UA法について
数値実験
実験結果
引用文献

SCE-UA法について

SCE-UA法（Shuffled Complex Evolution Method developed at the University of Arizona）は，シンプレックス法，ランダム探索，遺伝的アルゴリズムに類似した競争進化，集団混合の概念を組み合わせたパラメータの大域的探索手法である（Duan et al., 1992）．以下にSCE-UA法のアルゴリズムを示す（杉原ら, 2011）．

f:id:okome29taberu:20200809221811p:plain — SCE-UA法のアルゴリズム（杉原ら, 2011）

f:id:okome29taberu:20200809221814p:plain — SCE-UA法内のCCE法のアルゴリズム（杉原ら, 2011）

数値実験

SCE-UA法でのパラメータ探索において設定の必要な値は，探索するパラメータの数n，集団の個数p，各集団における個体（点）の数m，親個体の数q，反復回数α，βである．本稿では，nとpについては藤井（2011）と同様にn=10，p=10とした．ほかの値についてDuan et al. (1992) は，m=2n+1，q=n+1，α=1，β=2n+1を推奨している．この推奨値に準じて，m=21，q=11，α=1，β=21に設定した．

数値実験に用いるテスト関数として，藤井（2011）と同様に次の8関数を選択した．下表にこれらの最適解とその時のパラメータ値および定義域を示す．

f:id:okome29taberu:20200809221702p:plain — 数値実験に用いるテスト関数とその最適解，パラメータ値，定義域

収束判定基準は最良個体の関数値が10^-8を下回った場合に，パラメータは最適解に近づいたとして，探索成功とみなした．一方で，目的関数により評価した回数（評価回数）が840,000回を超えた場合は探索失敗とみなした（藤井, 2011）．各テスト関数での試行回数は100回とした．

実験結果

以下に探索成功回数と平均評価回数をテスト関数ごとに示す．参考に同様の条件で数値実験を行った藤井（2011）の結果も併記した．どのテスト関数についても100試行すべてで探索に成功した．また本稿で作成したプログラムと藤井（2011）の間で結果に大差はなく，再現できたといえる．

No	テスト関数	探索成功回数		平均評価回数
		本稿	藤井(2011)	本稿	藤井(2011)
1	Sphere	100	100	8077	7745
2	Ridge	100	100	10870	9966
3	Rosenbrock	100	100	18232	14662
4	Bohachevsky	100	100	10773	9325
5	Rastrigin	100	100	40891	37099
6	Schwefel	100	100	393456	423574
7	Griewank	100	100	14561	13071
8	Griewank-d	100	100	14813	13344

実装したプログラムはこちら（GitHub - Okome-chan/SCEUA）

引用文献

Duan, Q., Soroshin, S. and Gupta, V. K. (1992) Effective and Efficient Global Optimization for Conceptual Rainfall-Runoff Models, Water Resources Research, 28(4), pp.1015-1031.
杉原成満, 福田慎哉, 倉本和正, 荒木義則, 朝位孝二, 古川浩平 (2011) SCE-UA法を用いたタンクモデルの構築とそれを用いた土砂災害発生危険基準線の設定, 土木学会論文集F6 (安全問題), 67(1), pp.1-13.
藤井厚紀 (2011) 関数最適化アルゴリズム SCE-UA法の性能評価と改良, 福岡工業大学研究論集, 44(1), pp.45-52.

2020-05-02

【流出解析】タンクモデルのプログラム作成

はじめに
タンクモデルについて
基礎式
パラメータ
計算方法
プログラムと計算結果
引用文献

はじめに

先週は貯留関数法を用いて河川流域からの流出量を計算するプログラムを作成した（【流出解析】貯留関数法のプログラムを作ってみる - 雑記）．流出量を計算するモデルは貯留関数法以外にもあり，タイトル中のタンクモデルもその一つである．貯留関数法に続いて，本稿ではタンクモデルのプログラムを作成した．

タンクモデルについて

タンクモデルは，流域を，側面にいくつかの流出孔を持つタンクに置き換えて考える流出計算手法である．降雨はタンクモデルの最上段のタンクに入り，2段目以下のタンクは1段上のタンクの孔から水を受ける．各段のタンク内の水の一部は側面の孔から外部へ流出し，これらの和が河川流域からの流出量となる．

f:id:okome29taberu:20200502201131p:plain:w300

タンクモデルでは一般に直列3～4段のものが利用される．最上段が洪水の表面流出，2段目が中間流出，3，4段目が地下水流出に対応すると考えられている．この理由から，洪水解析用タンクモデルは2～3段，低水解析用タンクモデルは3～4段で構成されることが多い．

基礎式

上図で考えると，まず，1段目のタンクに降った雨 $r$ mm/hはタンクに一時的に貯留され， $e_1$ mm/hで蒸発散し，貯留量 $S_1$ mmと流出孔の高さ $h_1$ および $h_2$ mmに応じた流量 $q_{11}$ mm/h， $q_{12}$ mm/h， $g_1$ mm/hが，タンクからそれぞれ外部もしくは下部タンクへ流出する．このことから， $Δt$ 時間における1段目のタンクの貯留量変化 $ΔS_1$ は，

$ΔS_{1}=rΔt-e_{1} Δt-q_{11} Δt-q_{12} Δt-g_{1} Δt$

と表せられるので，以下の式が得られる．

$\frac{dS_1}{dt}=r-e_1-Q_1-g_1$

ただし， $Q_1=q_{11}+q_{12}$ であり，これが1段目のタンクの基礎式である．一方，2段目以降のタンクでは，水の供給源は降雨の代わりに上部タンクからの流入（浸透） $g_i$ であるので，上式と同様に考えると以下の基礎式が得られる．

$\frac{dS_i}{dt}=g_{i-1}-e_{i}-Q_{i}-g_{i}$

$Q_{i}=q_{i1}+q_{i2}+⋯+q_{ij}$

ここで， $S_i$ は $i$ 段目のタンクの貯留量， $e_{i}$ は $i$ 段目のタンクからの蒸発散量， $Q_{i}$ は $i$ 段目のタンクからの流出量， $q_{ij}$ は $i$ 段目のタンクにおける $j$ 番目の孔からの流出量である．一方で，流出量 $q_{ij}$ と貯留量 $S_i$ および孔の高さ $h_{ij}$ の間には以下の式の関係があると仮定する．

${\displaystyle q_{ij}= \begin{eqnarray} \left\{ \begin{array}{l} 0\hspace{38mm}(S_{i}≤h_{ij}) \\ a_{ij}(S_{i}-h_{ij})\hspace{12mm}(S_i>h_{ij}) \end{array} \right. \end{eqnarray} }$

ここで， $a_{ij}$ は $i$ 段目のタンクにおける $j$ 番目の孔のパラメータである．この式から， $i$ 段目のタンクにおいて，貯留量 $S_{i}$ が孔の高さ $h_{ij}$ よりも低い場合， $j$ 個目の孔からは流出しないということが分かる．同様に下部タンクへの流出量（浸透量） $g_i$ は，

$g_{i}=b_{i} S_{i}$

であらわされる．ここで， $b_{i}$ は $i$ 段目のタンクにおけるパラメータである．

パラメータ

本稿では，下図の菅原ら（1986）に記載されている河川用パラメータを用いる．

f:id:okome29taberu:20200502191032p:plain:w150

流出孔の高さ $h$ については，以下の通りである．

流出孔の位置	$h$ [mm]
1段目タンクの上側孔	25~60
1段目タンクの下側孔	5~15
2段目タンクの孔	0~30
3段目タンクの孔	0~60
4段目タンクの孔	0

計算方法

流出量 $Q_{i}$ および $q_{ij}$ ，浸透量 $g_{i}$ の計算式は前述したので，貯留量 $S_{i}$ の計算式だけ書き記す． $S_{i}$ についての基礎式において時間方向に前進差分をとると，

$\frac{S_{i}^{t+1}-S_{i}^{t}}{Δt}=g_{i-1}^{t}-e_{i}^{t}-Q_{i}^{t}-g_{i}^{t}$

となる．ここで， $t$ はタイムステップである． $S_{i}^{t+1}$ ついてまとめると，

$S_{i}^{t+1}=S_{i}^{t}+(g_{i-1}^{t}-e_{i}^{t}-Q_{i}^{t}-g_{i}^{t})Δt$

となる．ただし，1断面目のタンクにおいては $g^{t}_{i-1}$ の代わりに降水量 $r^{t}$ を用いる．実際の計算の際は，1段目のタンクから順に各孔からの流出量を計算した後に，各タンクの貯留量を計算するという手順となる．

プログラムと計算結果

本稿では，【流出解析】貯留関数法のプログラムを作ってみる - 雑記の際と同様に，降雨量データとして，財）北海道河川防災研究センター・研究所（2006）に記載されている湧別川・丸瀬布地点における2001年9月10日～15日に発生した洪水のデータを利用した．ただしモデルでの蒸発散量は0とした（用いた降雨量データは有効雨量であるから，蒸発散がすでに実際の降雨量から引かれていると仮定した）．また各流出孔のパラメータは菅原ら（1986）に従い，流出孔の高さはその中央値とした．

プログラムは以下の通りである．一応C++で書いた．

#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

using namespace std;

// Parameter setting
// 各タンクの流出孔の係数
double a11=0.20;    // 1段目タンクの上側孔
double a12=0.20;    // 1段目タンクの下側孔
double b1=0.20;     // 1段目タンクの下向き孔
double a21=0.050;   // 2段目タンクの孔
double b2=0.050;    // 2段目タンクの下向き孔
double a31=0.010;   // 3段目タンクの孔
double b3=0.010;    // 3段目タンクの下向き孔
double a41=0.001;   // 4段目タンクの孔

// 各タンクにおける流出孔の高さ(mm)
double h11=42.5;    // 1段目タンクの上側孔
double h12=7.5;     // 1段目タンクの下側孔
double h21=15.;     // 2段目タンクの孔
double h31=2.5;     // 3段目タンクの孔
double h41=0.;      // 4段目タンクの孔

// その他のパラメータ
double dt=1.;       // time step [h]

// 4段タンクモデル
void TankModel4steps(double& r, double& q, double& s1, double& s2, double& s3, double& s4){

  double q11=0.;                     // 1段目タンクの上側孔
  double q12=0.;                     // 2段目タンクの下側孔
  double g1=0.;                      // 1段目タンクの下向き孔
  double q21=0.;                     // 2段目タンクの孔
  double g2=0.;                      // 2段目タンクの下向き孔
  double q31=0.;                     // 3段目タンクの孔
  double g3=0.;                      // 3段目タンクの下向き孔
  double q41=0.;                     // 4段目タンクの孔

  // 各流出孔からの流出量の計算
  if(s1>h11) q11=a11*(s1-h11);    // 1段目タンクの上側孔
  if(s1>h12) q12=a12*(s1-h12);    // 1段目タンクの下側孔
  if(s2>h21) q21=a21*(s2-h21);    // 2段目タンクの孔
  if(s3>h31) q31=a31*(s3-h31);    // 3段目タンクの孔
  if(s4>h41) q41=a41*(s4-h41);    // 4段目タンクの孔
  g1=b1*s1;                       // 1段目タンクの下向き孔
  g2=b2*s2;                       // 2段目タンクの下向き孔
  g3=b3*s3;                       // 3段目タンクの下向き孔

  // 合計流出量の計算
  q=q11+q12+q21+q31+q41;

  // 各タンクの貯留量の計算
  s1+=dt*(r-q11-q12-g1);          // 1段目タンク
  s2+=dt*(g1-q21-g2);             // 2段目タンク
  s3+=dt*(g2-q31-g3);             // 3段目タンク
  s4+=dt*(g3-q41);                // 4段目タンク

}

int main (){

  char fileName[256];     // File name
  string ss;                  // Dammy line
  double s1,s2,s3,s4;   // 貯留量[mm]
  double q;             // 流出量[mm/h]
  double qobs;          // 実測流出量[mm/h]
  double r;             // 時間tからt+1における有降雨量[mm/h]
  double tmax=94;       // 計算時間[h]

  // 入力データファイルの設定
  ifstream ifs("./data.txt");
  if(!ifs){
    cout<<"cannot open file"<<endl;
    exit(1);
  }
  getline(ifs,ss);

  // 出力データファイルの設定
  ofstream ofs("./out.txt");
  ofs<<"time[h]\train[mm/h]\tQobs[mm/h]\tQsim[mm/h]"<<endl;

  //貯留量の初期値
  s1=0.;
  s2=0.;
  s3=0.;
  s4=0.;

  // 計算開始
  for(int t=1;t<tmax; t++){
    ifs>>ss>>r>>qobs;                   // 時間tからt+1における有降雨量の読み込み
    TankModel4steps(r,q,s1,s2,s3,s4);   // 4段タンクモデルでの計算

    // 計算結果の出力
    ofs<<std::fixed<<std::setprecision(4)<<t<<"\t"<<r<<"\t"<<qobs<<"\t"<<q<<endl;
    cout<<t<<endl;
  }
  return 0;
}

計算結果は下図の通りである．明らかに計算値は実測値に一致していないが，計算プログラムは合っている（はず）．これは用いたパラメータが今回の対象流域に対して適切でないためである．パラメータの最適化は今後の課題とする．

f:id:okome29taberu:20200502194041p:plain:w600

引用文献

菅原正巳, 渡辺一郎, 尾崎睿子, 勝山ヨシ子 (1986) パーソナル・コンピュータのためのタンク・モデル・プログラムとその使い方, 国立防災科学技術センター研究報告, 37.
（財）北海道河川防災研究センター・研究所（2006）実践流出解析ゼミQ&A編

2020-04-26

【流出解析】貯留関数法のプログラムを作ってみる

はじめに
基礎方程式の導出
数値解法
ニュートン法
入力データ
プログラムと計算結果
引用文献

はじめに

もともと流出解析はやっていたので，貯留関数法は知っていた（聞いたことがある程度だ）が，実際にその計算をしたことはなかった．そこで，少し時間と気持ちに余裕があるので，貯留関数法のプログラムを作ることにした．

貯留関数法とは，流域および河道を一つの貯水池と考えて，そこでの貯留量と流出量の間に一義的な関係性を仮定し，貯留量を媒介変数として，降雨量（流入量）から流出量を求める手法のことである．

基礎方程式の導出

貯留関数法では，河川流域に降った雨が流域に貯留された後に流出する過程を，流域を下図のような貯留タンクに見立てて考える．

f:id:okome29taberu:20200426202628p:plain:w400

流域（タンク）に降った降雨 $r$ mm/sはタンクに一時的に貯留され，貯留量 $S$ mmに応じた流量 $q$ mm/sがタンクから流出する．この時の降雨というのは，流域に降った雨のうち，樹木による遮断や地中による貯留の効果によって流域からの流出とならない損失分を除いたもの（有効雨量）である．ここで，貯留量 $S$ と流出量 $q$ の間には以下の式の関係があると仮定する．

$S=kq^p$ 　　（運動方程式）

ここで，k，pはパラメータである．次に，Δt時間におけるタンクの貯留量変化ΔSは，

$ΔS=rΔt-qΔt$

と表せられるので，以下の式が得られる．

$\frac{dS}{dt}=r-q$ 　　（連続式）

これら貯留関数法の運動方程式と連続式を連立して解くと，流域からの流出量を得られる．

数値解法

はじめに，貯留関数法の連続式を前進差分法で近似すると，

$\frac{S_{i+1}-S_{i}}{Δt}=r_{i+1}-\frac{q_{i+1}+q_{i}}{2}$

となり， $S=kq^{p}$ であるから，

$\frac{kq_{i+1}^{p}-kq_{i}^{p}}{Δt}=r_{i+1}-\frac{q_{i+1}+q_{i}}{2}$

となる．さらに， $q_{i+1}$ について整理すると，

$\frac{2k}{Δt }q_{i+1}^{p}+q_{i+1}-2(r_{i+1}-\frac{q_{i}}{2}+\frac{k}{Δt}q_{i}^{p})= 0$

となる． $q_{i+1}$ を求めるためには，この式を数値的に解けばよいので，本稿ではニュートン法で求める．

ニュートン法

ニュートン法は，主として多項式の実根を求めるために使用される．例えば， $x$ に関する方程式 $f(x)=0$ の根の近似値を $x_{1}$ とし，これを最初の試行値とすると，次の試行値 $x_{2}$ は $x_{1}$ における関数 $f(x_{1})=0$ の接線を求めればよい．そして，接線とx軸の交点が次の試行値 $x_{2}$ となる．接線の勾配を $m$ とすると，

$m=\frac{f(x_{1})}{x_{1}-x_{2}}=f'(x_{1})$

となるので，

$x_{2}=x_{1}-\frac{f(x_{1} )}{m}= x_{1}-\frac{f(x_{1})}{f'(x_{1})}$

以下同様に，

$x_{n+1}= x_{n}-\frac{f(x_{n})}{f'(x_{n})}$

として（ $x_{1}，x_{2}，…，x_{n+1}$ ）を求めていけば，最終的に真値に収束する．この式を用いて $q_{i+1}$ を求めるとすると， $q_{i+1}=x$ ， $r_{i+1}-\frac{q_{i}}{2}+\frac{k}{Δt}q_{i}^{p}=C$ とおいて，

$f(x)=\frac{2k}{Δt}x^{p}+x-2C$

$f'(x)=\frac{2kp}{Δt}x^{p-1}+1$

をそれぞれ使えばよい．

入力データ

貯留関数法による流出量の計算に必要なデータは，降雨量 $r$ およびパラメータ $k$ ， $p$ である．これらの値は解析対象とする洪水によって降雨量が異なるのはもちろんのこと，パラメータも流域に固有であるため各流域で異なる．本稿では，降雨量データの入手やパラメータの設定の手間を省くため，財）北海道河川防災研究センター・研究所（2006）に記載されている湧別川・丸瀬布地点における2001年9月10日～15日に発生した洪水のデータを利用する．そこでは，パラメータは $k=22.7602$ ， $p=0.6$ と設定されている．

プログラムと計算結果

本稿では，以下のようなプログラムをC++で書いた．ただC++は勉強し始めたばかりであるため，適切な書き方でないかもしれないが，ゆるしてほしい．

#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

using namespace std;

// Parameter setting
double k=22.7602;
double p=0.6;
double dt=1.;       // time step [h]
double lim=1.0e-6;  // 収束判断値

// f(x)
double f (double x, double r, double xold){
  return 2.*k/dt*pow(x,p) + x -2.*(r-xold/2.+k/dt*pow(xold,p));
}

// f'(x)
double df (double x){
  return 2.*k*p/dt*pow(x,p-1.)+1.;
}

// Newton Raphson method
double NewtonRaphson (double x, double r){

  int i=0;
  int imax=10000;   // 最大反復回数
  double xn;        // 近似値＝次の試行値
  double xold=x;    // 時間tにおける流出量[mm/h]
  
  // 反復開始
  while(i<=imax){
    i++;
    xn=x-f(x,r,xold)/df(x); 
    if(fabs(xn-x)<lim) break;   // 値が収束していた場合は反復から抜ける
    x=xn;                       // 近似値を次の試行値に置き換える
  }
  x=xn;
  return x;
}

int main (){

    char fileName[256];       // File name
    string ss;                // Dammy line
  double q;             // 時間tにおける流出量[mm/h]
  double qn;            // 時間t+1における流出量[mm/h]
  double qobs;          // 実測流出量[mm/h]
  double r;             // 時間tからt+1における有降雨量[mm/h]
  double tmax=94;       // 計算時間[h]

  // 入力データファイルの設定
  ifstream ifs("./data.txt");
  if(!ifs){
    cout<<"cannot open file"<<endl;
    exit(1);
  }
  getline(ifs,ss);

  // 出力データファイルの設定
  ofstream ofs("./out.txt");
  ofs<<"time[h]\train[mm/h]\tQsim[mm/h]\tQobs[mm/h]"<<endl;


  //流出量の初期値
  q=0.0001;
  
  // 計算開始
  for(int t=1;t<tmax; t++){
    ifs>>ss>>r>>qobs;       // 時間tからt+1における有降雨量の読み込み
    qn=NewtonRaphson(q,r);
    q=qn;                   // 流出量の値の更新

    // 計算結果の出力
    ofs<<std::fixed<<std::setprecision(4)<<t<<"\t"<<r<<"\t"<<q<<"\t"<<qobs<<endl;
    cout<<t<<endl;
  }

  return 0;
}

上記のプログラムを用いた計算結果は以下の通りである．定量的な評価を本稿では実施していないので，強くは言えないが，計算値は，大まかではあるものの実測値を再現できているように見える．この計算結果以上に再現性を上昇させるためには，パラメータを最適化するなどの手続きが必要であろう．

f:id:okome29taberu:20200426211816p:plain:w500

引用文献

（財）北海道河川防災研究センター・研究所（2006）実践流出解析ゼミQ&A編