タンクモデル

タンクモデルのパラメータ最適化に、SCE-UAでなく、EXCEL ソルバーを使おうと考えました。

理由は「簡単だから」＆「文献であったから」。

わざわざプログラミングしなくても、大域的最適解（らしきもの）をソルバーで求められます。四則演算のみで負荷のかかる計算ではありません。VBA も必要ないので、ペナルティ関数やタンクの段数、横穴の数などカスタマイズが容易になります。パラメータ推定が終わればグラフも自動で更新されるので EXCEL の方が数段便利。

気象庁のホームページに3段の式とパラメータがあります。まずはこれを手本に組んでみました。が、タイムステップの取り方がイマイチわかりにくい表現になっており、また配信されている土壌雨量指数は補正がかかっているので、正しく組めているのか検証できません。
何か検算できる資料はないか、と古い資料「流出計算例題集2」を引っ張り出してきました。が、この計算例も時間の取り方が怪しい。

ま、基本に返りましょうと、菅原正巳「流出解析法」を見直し。というか、以前は眺めただけでしたね。数式部分を少し補足しながら書き下しておきましょう。


*********************************************
降雨を0とすると貯留高Sの時間変化（減少）は
-dS/dt=αS+βS=(α+β)S

λ=α+βとおくと、放射性元素の半減期と全く同じ式となり、変数分離型の解法が使えます。
-dS/dt=λS
-dS=λSdt
∫(1/S)dS＝-λ∫dt
lnS=-λt+c

t=0のときS=S0とおくと
S=S0e^（-λt）
S0=e^c

半減期t=TのときS=S0/2で
1/2=e^（-λT）
2=e^λT
ln2=λT
T=0.693/λ

つまり、（α+β）によって半減期が決まる。
というか、タンクモデルは指数関数型で多段の場合はその重ね合わせになっているのね。今まで理解していませんでした。


Δtにおける貯留高の変化は
S0-S0e^-λΔｔ
=S0(1-e^(-λΔｔ)）

Δｔ当たりの流出高は、マクローリン展開で
S0(1-e^(-λΔｔ)）（α/λ）
=αS0（1/λ-(e^(-λΔｔ)/λ）
=αS0（1/λ-(1-λΔｔ/1!+λ^2Δｔ^2/2!-λ^3Δｔ^3/3!・・・・)/λ）
=αS0（1/λ-(1/λ-Δｔ/1+λΔｔ^2/2-λ^2Δｔ^3/6・・・・)）
=αS0（Δｔ-λΔｔ^2/2+λ^2Δｔ^3/6・・・・)）

α+βやΔｔが小さいときは2項以下を無視して
≒αS0Δｔ

単位時間当たりの流出高の変化は
≒αS0

ってことでしょう。実測流量にあわせこむ際には、数項分の影響を考慮したαが得られるということ。
***********************************************

理屈では、タイムステップの取り方に間違いはありません。組み込みの際のミスチェックに何らかの検証は必要ですが、ひとまず、ペナルティ付き・1~4段タンク・最大18パラメータの最適化付き EXCEL シートができました。