Transcript Identification of the Manning’s Roughness in the CCHE2D

Optimal Estimation of Bottom
Friction Coefficient for Free-Surface
Flow Models
Yan Ding, Ph.D.
Research Assistant Professor, National Center for Computational Hydroscience and
Engineering, The University of Mississippi, University, MS 38677
Presentation for ENGR 691-73, Numerical Optimization, Summer Session, 2007- 2008
Outline
• Introduction and Brief Review of Optimal Parameter
Estimation
• Optimal Estimation of Manning’s Roughness
Coefficient in Shallow Water Flows
• Variational Parameter Estimation of Manning’s
Roughness Coefficient Using Adjoint Sensitivity
Method
• Discussions
I. Introduction and Brief Review of Optimal
Parameter Estimation
Parameters in Mathematical Models
• Physical Parameters: measurable
describe physical properties and simple physical processes
have state equations which describe the physical process, e.g.,
=(T, P)
• Empirical Parameters: unmeasurable
complicated physical processes, unmeasurable from the physical
point of view, e.g., Manning Roughness n,
without state equation
have to be determined from historical observations
Complexity of Manning’s Roughness
Physical-process-based parameters
Depth (m)
0
250
500m
Bridge Pier
1.0m/s
Manning Roughness n = f(surface friction,
bed form friction, wave, flow unsteadiness
,vegetation (?), etc) [1] , or
n=n(x,y) : distributed parameters,
spatially dependent
Moyie River at Eastport, Idaho
H
10
9
8
7
6
5
4
3
2
1
0
Objectives of Presentation
• Theoretical Framework for estimation of the Manning’s
roughness by using optimal theories (Optimization)
• Estimate the Manning’s roughness so as to minimize the
discrepancy between computation and observation
• Compare the efficiency and accuracy of various minimization
procedures (High Performance of Identification)
• Demonstrate practical parameter estimation to determine the
Manning’s roughness for the CCHE1D and CCHE2D
hydrodynamic models (Application)
General Analysis Tools Based on Optimal Theories
Performance Function
Optimal Theory
Minimization
Procedure
e.g.,Weight Least-Square Method
1.
2.
3.
1.
2.
3.
4.
5.
6.
Sensitivity Analysis
Maximum likelihood
Extended Kalman Filter
CG Method
Newton Method
LMQN (BFGS, LBFGS, etc)
Trust
Region Method
Sakawa-Shindo
Method
Linear Programming
Quadratic Programming
Ill-Posedness of the Problem of Parameter Estimation
The inverse problem for parameter estimation is often illposed. This ill-posedness is characterized by:
• Non-uniqueness (Identifiability)
The concept of identifiability addresses the question whether
it is at all possible to obtain unique solutions of the inverse
problem for unknown parameters of interest in a
mathematical model from data collected in the spatial and
time domains
• Instability of the identified parameters
Instability here means that small errors in data will cause
serious errors in the identified parameters
History of Parameter Identification in
Computational Hydroscience and Engineering
• 1980sSakawa and Shindo (1980): a general optimal control algorithm
Bennett and McIntosh (1982) and Prevost and Salmon (1986) :
early work of variational method for tidal flow
Heemink (1987) : EKF to flow identification (English channel)
Panchang and O’Brien (1989) : Adjoint parameter identification
(API) for bottom friction in a tidal channel
Das and Lardner (1990, 1992) : bottom friction in a tidal
channel
History of Parameter Identification in
Computational Hydroscience and Engineering
• 1990sYu and O’Brien (1991) : wind stress coefficient
Richardson and Panchang (1992) : API for eddy viscosity,
Zhou et al (1993): LMQN for meteorology
Kawahara and Goda (1993) : eddy viscosity by maximum likelihood
Lardner and Song (1995): eddy viscosity profile in a 3-d channel
Hayakawa and Kawahara (1996) : velocity and phase by Kalman filter
Atanov et al (1999): Roughness in open channel (de St. Venant equation)
Ding & Wang (2003): API for Manning’s roughness coefficient in channel
network in the CCHE1D model
Ding & Wang (2004): Sensitivity-equation method for estimation of spatiallydistributed Manning’s n in shallow water flow models (CCHE2D)
Reviews and Books about Parameter
Identification
• Review
Ghil and Malanotte-Rizzoli (1991), Advances in Geophysics,
vol.33 : Review about data assimilation in meteorology and
oceanography,
Navon (1998), Dynamics of Atmosphere and Oceans, vol 27:
Review about parameter identification, identifiability and
nonuniqueness related to ill-poseness of parameter
identification problem
• Books:
Malanotte-Rizzoli, P. (1996): Modern Approaches to Data Assimilation in
Ocean Modelling, Elsevier Oceano. Series
Wunsch, C. (1996): The Ocean Circulation Inverse Problem, Cambridge
University Press.
Nocedal & Wright (1999), Numerical Optimization, Springer, NY
Journals and Online Resources about Optimal Parameter
Estimation
• Journals
SIAM J. on Control and Optimization /Optimization
J. Acoust. Soc.
J. Geophysical Research
Water Resource Research (American Geophysical Union)
Monthly Weather Review
Tellus
ASCE J. Hydr. Engrg.
(OOOOO)
(OOOO)
(OOOO)
(OOO)
(OOO)
(OO)
(OO)
• Web Resources
SPIN Database Online ( http://scitation.aip.org/jhtml/scitation/search/)
American Institute of Physics, NY
ISI Web of Science (http://scientific.thomson.com/products/wos/ )
EI (Engineering Village) http://daypass.engineeringvillage.com/
Output Error Measure:
Performance Function J
• To evaluate the discrepancy between computation and
observation, one can define a weighted least squares form as
1 tf
J ( X )   ( X  X obs )T W ( X  X obs )dt
2 t0
or its discrete form:
1 M N
obs T
obs
J ( X )    ( xi , j  xi , j ) W ( xi , j  xi, j ) t j
2 j 1 i 1
where X = physical variables, the weight W is related to confidence in
the observation data. The superscript ‘obs’ means measured data. t0 :
the starting time, tf : the final time during the period of parameter
identification
Mathematical Framework for Optimal Parameter Estimation
• The parameter identification is to find the parameter n
satisfying a dynamic system such that
1 tf


f (n)  min  J ( X ) : J ( X )   ( X  X obs )T W ( X  X obs )dt 
2 t0


subject to
X
 F ( X , n, t )
t
• Local minimum theory[3] :
Necessary Condition: If n* is the true value, then J(n*)=0;
Sufficient Condition: If the Hessian matrix  2J(n*) is
positive definite, then n* is a strict local minimizer of f
Establish an Iterative Procedure to Search for Optimal
Parameter Value: Sakawa-Shindo Method (1)
• Assumed that the computational domain  =l
(l=1,2,..L) in which the Manning roughness is nk,
by adding a penalty term, the performance
function can be modified as the following form,
J
*( k )
1 t f ( k 1)
( X , nl )  J ( X )   (nl
 nl( k ) )T C ( k ) (nl( k 1)  nl( k ) )dt
2 t0
(l )
where (k) is iteration step in the
Sakawa-Shindo procedure, C(k)
is a diagonal weighting matrix at
the (k)th step, the weight matrix
is renewed at every step.
n1
nl
nL
Fig. Spatially-distributed parameters
Sakawa-Shindo Procedure (2)
• Correction Formulation of nk
Considering the first-order necessary condition
 J(n)=0,
the parameter can be estimated from the previous value as
follow,
T
 X 
n
 n  (C )    W ( X  X obs )dt
t0
 nl 
• Change the weighting matrix C(l) :
( k 1)
l
(k )
l
( k ) 1
tf
If J(k+1)J(k), then C(k+1)=0.9C(k),otherwise C(k)=2.0C(k),
continue to do next estimation of n.
Sensitivity Analysis
• Evaluate the gradient of performance function
J ( X , n)  
tf
t0
 X 
obj
W
(
X

X
)dt
 n 
T
the sensitivity matrix can be computed by
1. Influence Coefficient Method [2] : Parameter perturbation
trial-and-error
2. Sensitivity Equation Method: Forward computation
Solve the sensitivity equations obtained by taking the partial
derivatives with respect to each parameter in the governing equations
and all of boundary conditions
3. Adjoint Equation Method: Backward computation
The sensitivity is evaluated by the adjoint variables governed
by the adjoint equations derived from the variational analysis.
Minimization Procedures for
Unconstrained Optimization
Trust Region Methods: change searching directions and step size
• Sakawa-Shindo method[4]
considering the first order derivative of performance function only, stable in most of
practical problems
Linear Search Approaches: directions are pre-determined, change step size
• Conjugate Gradient Methods (e.g. Fletcher-Reeves method) [5]
The convergence direction of minimization is considered as the gradient of
performance function.
• Truncated Newton Methods or Hessian-free Newton methods
• Limited-Memory Quasi-Newton Method (LMQN) [3]
Newton-like method, applicable for large-scale computation, considering the second
order derivative of performance function (the approximate Hessian matrix)
Unconstrained optimization methods
nk 1  nk   k dk
k :
Step length
1) Line search
dk :
Search direction
n2
Contour of J
.n*
2) Trust-Region algorithms
d1
n1
Searching in Parameter Space
Methods for Unconstrained Optimization
1) Steepest descent (Cauchy, 1847)
dk  J (nk )
2) Newton
1
d k   J (nk ) J (nk )
2
3) Quasi-Newton (Broyden, 1965; and many others)
dk  H k J (nk )
H k   2 J (nk )1
Methods for Unconstrained Optimization (cont.)
4) Conjugate Gradient Methods (1952)
dk 1   gk 1  k sk
 k is known as the conjugate gradient parameter
sk  nk 1  nk
gk  J (nk )
5) Truncated Newton method (Dembo, et al, 1982)
1
d k   J (nk ) g k
2
6) Trust Region methods
rk   2 J (nk )d k  g k
Step Length Computation
1) Armijo rule:
f ( xk )  f ( xk   m k d k )    m k f ( xk )T d k
  (0,1)
  (0,1/ 2)
 k  (f ( xk ) d k ) / d k
T
2
2) Goldstein rule:
1 k gkT dk  f ( xk   k d k )  f ( xk )  2 k g kT d k
0  2  12  1  1
3) Wolfe conditions:
f ( xk   k d k )  f ( xk )   k g d k
T
k
f ( xk   k d k ) d k   g d k
T
f ( xk  k dk ) dk   g dk
T
T
k
0    1
Implementations:
Shanno (1978)
Moré - Thuente (1992-1994)
T
k
Remarks:
1) The Newton method, the quasi-Newton or the limited
memory quasi-Newton methods has the ability
to accept unity step lengths along the iterations.
2) In conjugate gradient methods the step lengths
may differ from 1 in a very unpredictable manner.
They can be larger or smaller than 1 depending
on how the problem is scaled*.
*N. Andrei, (2007) Acceleration of conjugate gradient
algorithms for unconstrained optimization.
(submitted JOTA)
Limited-Memory Quasi-Newton Method
(LMQN) (Basic Concept 1)
Given the iteration of a line search method for parameter n
nk+1 = nk + kdk
k = the step length of line search
sufficient decrease and curvature conditions
dk = the search direction (descent direction)
dk  Bk1J k (nk )
n2
Bk = nn symmetric positive definite matrix
d1
For the Steepest Descent Method: Bk = I
Newton’s Method: Bk= 2J(nk)
Quasi-Newton Method:
Bk= an approximation of the Hessian  2J(nk)
Contour of J
.n*
n1
L-BFGS (One of LMQN method) (1)
• Difficulties of Newton’s method in large-scale optimization
problem: obtain the inverse Hessian matrix, because
the Hessian is fully dense, or,
the Hessian cannot be computed.
• BFGS (Broyden, Fletcher, Goldfarb, and Shanno, 1970)
Constructs the inverse Hessian approximation H k  Bk,1
Hk 1  VkT HkVk  k sk skT
1
k  T ,
Vk  I  k yk skT
yk sk
sk  nk 1  nk ,
yk  J k 1  J k
However, all of the vector pairs (sk, yk) have to be stored.
L-BFGS (2)
• Updating process of the Hessian by using the information
from the last m Q-N iterations
H k 1  (vkT  vkTmˆ ) H 0 (vk mˆ  vk )
  k mˆ (vkT  vkTmˆ 1 ) pk mˆ pkTmˆ (vk mˆ 1  vk )
  k mˆ 1 (vkT  vkTmˆ 2 ) pk mˆ 1 pkTmˆ 1 (vk mˆ 2  vk )

  k pk pkT
Where m is the number of Q-N updates supplied by the user,
Generally [3] ,
ˆ  min{k , m  1}
3  m  7 and m
Only m vector pairs (si, yi),i=1, m, need to be stored
Set the initial n=( n1,n2,…nL )T
k=0
Solve the initial state vector X0
Calculation of performance
function J0, gradient g0, and search
direction d0
L-BFGS
Calculation of
Sensitivity Analysis
n k 1  n k   k d k
Solve the state vector Xk+1
Flow chart of
parameter identification
by using L-BFGS
procedure
Shallow
Water
Equations
Line Search
Shallow
Water
Equations
Calculation of Jk+1
g k 1   J ( n k 1 )
Calculation of
Sensitivity Analysis
d k 1   H k 1 g k 1
Calculation of
Update Hessian matrix
by the recursive
iteration
Yes
||gk+1||max{1,||nk+1||}
nl  nl
Max ( k 1 l k )   n
nk
Yes
Stop
J k 1   J
Yes
No
L-BFGS-B
• The purpose of the L-BFGS-B method is to minimize the
performance function J(n) , i.e.,
min J(n),
subject to the following simple bound constraint,
nmin  n  nmax,
where the vectors nmin and nmax mean lower and upper bounds
on the Manning roughness.
• L-BFGS-B is an extension of the limited memory algorithm
(L-BFGS) for unconstrained problem proposed in [6].
• For the details about the limited memory algorithm for
bound constrained optimization method, see [7].
Scaling Problems
• When the difference of the identified parameters is very large,
the identification suffers from the poor scaling problem. 
lead to instability in the identification process
• Scaling : parameter transformation  n/=Dn
Gill et al [8] proposed the diagonal matrix D
1  J (n)
Di 
2 J (n)
But we found this approach could not do very well, even
destabilize the identification process (L-BFGS). Therefore, we
proposed the two approaches for forming the matrix D
Scaling 1
1
Di 
2 J (n) rms
Scaling 2
1  J (n)
Di 
2 J ( n ) rms
II. Optimal Estimation of Spatially-Distributed
Manning’s Roughness Coefficient in Shallow
Water Flows
Application to Shallow Water Equations
• Governing Equations

 [( h   )u j ], j  0
t
gn2 urur
ui
 u j ui , j   g,i   e (ui , j  u j ,i ), j 
u in ,
4/3 i
t
(h   )
where n is the Manning roughness coefficient, is to
be identified.
Boundary and Initial Conditions
Given the boundary S=S1S2, the boundary conditions are
subject to
ui  uˆi
on S1,
t   (u  u )n  tˆ on S2,
i
e
i, j
j ,i
j
i
The initial conditions are
ui  ui0
 
0
n1
at t=t0,
nl
nL
Sensitivity Equations for Forward
Computation
•
The sensitivity equations of sensitivity vector
  u v 
 n , n , n ,
k
k
 k
•
k  1,2, L,
can be derived, by taking the partial derivative with respect to each parameter
in the governing equations.
   
u

u


  (h   )( i ),i 
ui ,i  ( i )(h   ),i  ui (
),i  0
t  nk 
nk
nk
nk
nk
  
  ui 
  ui  u j ui


  u j

 
   e
 g 
t  nk 
x j  nk  nk x j
x j
 nk ,i
gn2 ur ur ui
2 gnk ur ur



ui  0 ,
4/3
4/3
h
nk
(h   )
   ui    u j 

 



 x j  nk  xi  nk 
1
0
 
within k
out of k
Boundary and Initial Conditions
• Boundary Conditions
ui
 0 on S1 ,
nk
ti
 0 on S2
nk
• Initial Conditions
 ui

 nk
0

  0,

 

 nk
0

  0 at t  t0

This becomes a well-posed problem for the identification of
Manning Roughness in shallow water flows
Numerical Techniques and Program Properties
• Numerical method to solve the sensitivity
equations: Efficient Element Method (CCHE2D)
• Structural modules for parameter identification are
compatible with CCHE2D;
• The user can choose an appropriate minimization
procedure (Sakawa-Shindo method, L-BFGS, or
L-BFGS-B) (Multifunction)
• Easy to prepare the input and control data files
Data Flows for Parameter Identification
Input Data
Output Data
Input data for flow model
Output data for flow model
Observation data:
Filename: *.obs
Initial control data for
parameter : Filename: *.man
Control data of L-BFGS-B:
Filename: *.lbfgsb
Model of Parameter
Identification of
Manning Roughness in
the CCHE2D
Results of identified
parameter: Filename: *.para
iterate.dat
Results of performance
function: Filename: *.per
History output for each
observed point: *.opt_his
Final sensitivity results:
Filename: *_sen.plt
Numerical Results (Verification)
• Strategy for verification
Specified n
Flow Model
(CCHE2D)
Computed
Flow Flied
Parameter
Identification
Create
observation
data
n=n* ?
Identified n*
Initial Values of parameters : far from the true values
Verify : identifiability, uniqueness of solution,
stability of process
Observation Data and Identified Parameters
(Open channel flow: N=3)
Water Elevation (m)
10.05
(a) Longitudinal profile of water elevation
10.04
10.03
n1=0.01
n2=0.02
10.02
n3=0.03
10.01
10.00
0
250
500
X (m)
0
Observed Point
n1=0.01
750
n2=0.02
1000
100
n3=0.03
(b) Velocity vectors and observed points
200m
Cases for verification of parameter
identification (Open channel flow)
CASE
Identification Method
No.
Objective
Initial
Initial Value
Total
Bound
Values of
Values
of Weighting
Observed
Constraint
Manning’s n
of n
Matrix C
Points*
of n
1
Sakawa-Shindo
0.01, 0.02, 0.03
0.005
10.0
6
No
2
LBFGS
0.01, 0.02, 0.03
0.005
N/A
6
No
3
LBFGS + Scaling 1
0.01, 0.02, 0.03
0.005
N/A
6
No
4
LBFGS + Scaling 2
0.01, 0.02, 0.03
0.005
N/A
6
No
5
L-BFGS-B
0.01, 0.02, 0.03
0.005
N/A
6
[10-5, 1.0]
6
L-BFGS-B+ Transition
0.01, 0.02, 0.03
0.005
N/A
6
[10-5, 1.0]
*) The observed variables at each point have three, i.e., two velocity components and water elevation.
Number of Q-N updates in L-BFGS: m=5
Case 1 : Sakawa-Shindo Method
n1
n2
n3
Objective Value
n=0.03
0.03
10-3
Performance function
0.04
10-2
Manning's n
10-4
n=0.02
0.02
n=0.01
0.01
0
10
20
30
40
50
60
Iterations of Sakawa-Shindo loop
(a) Iterations of Manning’s n
70
10-5
10
-6
10
-7
10
-8
10
20
30
40
50
60
70
Iterations of Sakawa-Shindo loop
(b) Iterations of performance function
Case 2 : L-BFGS Method
102
100
10-1
10
-2
10-4
Objective Value
n=0.03
n=0.02
10
10
Performance function
Manning's n
n1
n2
n3
Sakawa-Shindo
LBFGS without Scaling
0
n=0.01
-2
10
-6
10-8
10-10
10
20
30
40
50
Iterations of L-BFGS
(a) Iterations of Manning’s n
60
10
-12
10
20
30
40
50
60
70
Iterations of identification process
(b) Iterations of performance function
Case 3 and Case 4: L-BFGS +Scaling (1)
10
100
0
n1
n2
n3
10-1
Objective Value
n=0.03
Manning's n
Manning's n
n1
n2
n3
10-1
Objective Value
n=0.03
n=0.02
n=0.02
n=0.01
10-2
10
20
30
40
50
Iterations of L-BFGS
(a) L-BFGS+Scaling 1
60
n=0.01
10-2
10
20
30
40
Iterations of L-BFGS
(b) L-BFGS+Scaling 2
50
60
Case 3 and Case 4: L-BFGS +Scaling (2)
Norm of error = ||n -nobj||2
100
102
L-BFGS without scaling
Scaling 1
Scaling 2
10
10-1
-2
10-2
Norm of error
Performance function
100
10-4
10-3
10
-6
10-8
10-10
10
LBFGS without scaling
LBFGS with Scaling 1
LBFGS with Scaling 2
10
-4
10
-5
10
-6
-12
10
20
30
40
Iterations of L-BFGS
(a) Comparison of J
50
60
10
20
30
40
50
60
Iterations of L-BFGS
(b) Comparison of norm of error
70
Case 5 and Case 6: L-BFGS-B
Bound constraint: 0.00001  n 1.0
100
n1
n2
n3
10-1
Objective Value
n=0.03
n=0.02
n=0.01
10-2
10
20
30
40
Iterations of L-BFGS
(a) L-BFGS-B
50
60
n1
n2
n3
Manning's n
Manning's n
100
10-1
Objective Value
n=0.03
n=0.02
n=0.01
10-2
10
20
30
40
Iterations of L-BFGS
(b) L-BFGS-B+Transition
50
60
Comparison of performance function
10
2
Sakawa-Shindo
LBFGS
LBFGS + Scaling 1
LBFGS + Scaling 2
LBFGS-B
LBFGS-B+Transition
100
Performance function
10-2
10
-4
10-6
10
-8
10-10
10
-12
10
20
30
40
50
Iterations of Identification process
60
70
Comparison of norm error ||n -nobj||2
100
-1
10
-2
Norm of error
10
Sakawa-Shindo (Case 1)
LBFGS (Case 2)
LBFGS + Scaling 1 (Case 3)
LBFGS + Scaling 2 (Case 4)
LBFGS-B (Case 5)
LBFGS-B+Transition (Case 6)
10-3
10-4
10
-5
10
-6
10
20
30
40
50
Iterations of identification process
60
70
Comparison of norm of gradient J (n) 2
100
L-BFGS + Scaling 1
L-BFGS + Scaling 2
L-BFGS-B
L-BFGS-B+TRANSITION
10-1
Norm of gradient
10-2
10-3
10
-4
10
-5
10
-6
10
-7
5
10
15
20
25
30
35
Iterations of performance function evaluation
40
Summary of Verification for Open Channel
Cases
Case
No.
Total
log10 J
Total
Averaged
number of
CPU
CPU time of
function
time (s)
one iteration
evaluation
n  n obj
Identification
2
Method
(s)
1
67
-7.517
67053.6
1000.8
1.82610-3
Sakawa-Shindo
2
52
-10.169
60239.5
1158.5
1.93910-5
L-BFGS
3
32
-10.241
37122.6
1160.1
2.10410-5
L-BFGS+Scaling 1
4
41
-9.966
47540.0
1159.5
7.55010-5
L-BFGS+Scaling 2
5
30
-10.155
34832.8
1161.1
2.40210-5
L-BFGS-B
6
31
-9.680
35983.7
1160.8
1.31110-5
L-BFGS-B
+ Transition
Remark: The L-BFGS-B with bound constraints has excellent features for
stabilizing the identification process and accelerating the convergence.
Application: East Fork River
n5
n4
Partition of
Manning’s n
L=5
n3
n2
n1
Length of study reach = 3.3km
0.15
Manning's n
History of Manning’s
n at each river reach
0.05
0
Manning's n
Procedure: L-BFGS
Problem: n2 < 0
(unphysical meaning)
n1
0.1
5
10
15
0
-0.2
-0.4
n2
-0.6
-0.8
-1
5
10
15
Manning's n
0.15
0.05
0
5
10
15
Manning's n
0.15
n4
0.1
0.05
0
5
10
15
0.15
Manning's n
Initial Values
n1=0.05
n2=0.01
n3=0.02
n4=0.03
n5=0.08
n3
0.1
n5
0.1
0.05
0
5
10
Iterations of L-BFGS
15
Manning's n
History of Manning’s
n at each river reach
1
0.8
0.4
0.2
0
15
20
n2
0.01
0.005
0
5
10
15
20
1
0.8
n3
0.6
0.4
0.2
0
5
10
15
20
Manning's n
1
0.8
n4
0.6
0.4
0.2
0
5
10
15
20
1
Manning's n
Initial Values
n1=0.05
n2=0.01
n3=0.02
n4=0.03
n5=0.08
10
0.015
Manning's n
Bound constraint of n:
0.00001  n  1.0
5
0.02
Manning's n
Procedure: L-BFGS-B
n1
0.6
0.8
n5
0.6
0.4
0.2
0
5
10
15
Iterations of L-BFGS-B
20
Comparison of Water Elevation
8
Ideal Result
optimized Result
Computed Water Elevation (m)
7.5
7
6.5
6
5.5
5
5
5.5
6
6.5
7
Observed Water Elevation (m)
Procedure: L-BFGS-B
7.5
8
Comparison of Identified Parameters : East Fork River
Identification
n1
n2
0.05790
0.00001*
n3
n4
n5
Method
SakawaShindo
L-BFGS
L-BFGS-B
0.02144 0.03128 0.08249
0.05615 - 0.00026** 0.02493 0.03283 0.08289
0.05839
0.00489
0.02525 0.03995 0.08915
*) The parameter was forced to reach the lower bound of Manning’s n.
**) The parameter is unphysical.
Remarks on Optimization Algorithms
• Sakawa-Shindo Method
Stable in most of practical cases, ‘rigid’ feature of convergence
in multi-parameters identification, the solution may not be
optimal or minimizer in some cases
• LMQN methods:
Fast, effective for large-scale parameter identification
• Scaling Technique for L-BFGS
Necessary to stabilize the process of the identification process
• L-BFGS-B: the excellent minimization procedure
not only guarantee the demand of bound constraints, keep the
parameters from being unphysical,but also stabilize and speed
up the identification process.
Questions? and Discussion?
III. Variational Parameter Estimation of
Manning’s Roughness Coefficient Using
Adjoint Sensitivity Method
Outline
• Optimal Estimation of Manning Roughness
Coefficient in channel network (Theory
Framework)
Definition of Objective Function
Adjoint Sensitivity Analysis of 1D Saint
Venant Equations
Parameter Identification Procedure
• Numerical Techniques
• Verification and Application
• Conclusions and Further Topics
Objectives
Theoretically,
• Establish a theory framework for identification of the
Manning’s n so as to minimize the discrepancy
between computation and observation using adjoint
approach
• Find an efficient procedure (Ding et al 2003)
For Engineering Applications,
• Enable optimal estimation of the parameter in natural
river which may be single or channel network
• Give a distribution of the parameter along river
• Provide an analysis tool to determine the Manning’s
roughness for the CCHE1D model (Watershed Model)
de Saint Venant Equations
- Dynamic Wave
A Q
L1 

q0
t x
  Q    Q
L2      2
t  A  x  2 A
2

Z
  g
 gS f  0
x

n 2Q | Q |
S f  2 4/3
AR
where Q = discharge; Z=water stage;
A=Cross-sectional Area;
q=Side discharge;
=correction factor;
R=hydraulic radius
n = Manning’s roughness identified in the study
Initial Conditions and Boundary Conditions
I.C. (Base Flow)
Q t 0  Q ( x,0), x  [0, L]
A t 0  A( x,0), x  [0, L]
B.C.s
Q x0  Q(0, t ), t [0, T ]
Upstream
Downstream
or
(Hydrograph)
Z xL  Z ( L, t ), t [0, T ]
Q xL   (Z )
(Stage-discharge rating curve)
or open downstream boundary
Output Error Criterion
- Objective Function J
To evaluate the discrepancy between computations and
observations, a weighted least squares form is defined as
1 T L
J (Q, Z , n )    WQ (Q  Q obs ) 2  WZ ( Z  Z obs ) 2 dxdt ,
2 0 0
where Q = discharge; Z=water stage; the weight W is
related to confidence in the observation data. The
superscript ‘obs’ means measured data; T=the period of
parameter identification; L=length of channel
Mathematical Framework for Parameter
Identification
• The parameter identification is to find the parameter n
satisfying a dynamic system such that
f (n)  min(J (Q, Z , n)),
where Q and Z are satisfied with the continuity equation
and momentum equation, respectively (i.e., de Saint
Venant Equations)
• Local minimum theory :
Necessary Condition: If n* is the true value, then J(n*)=0;
Sufficient Condition: If the Hessian matrix  2J(n*) is
positive definite, then n* is a strict local minimizer of f
Variational Analysis
- To Obtain Adjoint Equations
Extended Objective Function
T
J  J 
*
0
L
 ( L   L )dxdt
0
A 1
Q 2
where A and Q are the Lagrangian multipliers
Necessary Condition (Stationary Condition)
J *  J  0
on the conditions that
T
L1 (Q,Z )0
L2 (Q,Z )0
t
D
C
O
L
B
{
A
x
Fig. Solution domain
Variation of Extended Objective Function
J * 
T
L
0
0

(
f Z
f
f
δA 
Q  n )dxdt
Z A
Q
n
2/3
A g Q Q Q Q 2 Q 2 g (1   ) R Q | Q | Q
  (
 *

 3

)Adxdt
0 0
t
B x A2 t
A x
nK 3
T L 
1 Q Q Q 2 g | Q | Q
  ( A 
 2

)Qdxdt
2
0 0
x A t
A x
K
T L 2 gQ | Q | Q
 
ndxdt
2
0 0
nK
Q
Q 2Q
Q
Q
  [( A  2 Q )A  Q ]dx  [ 

A

(


Q )Q ]dt
A
3
2
A
A
A
A
T
L
0
where
AR2 / 3
K
;
n
2 B*  4 R

;
*
3B
B*  Top Width
Adjoint Equations for
the Full Nonlinear Saint Venant Equations
According to the extremum condition, all terms multiplied by
A and Q can be set to zero, respectively, so as to obtain the
equations of the two Lagrangian multipliers, i.e, adjoint equations
QWQ
A Q A g Q 2 g (1   )Q | Q | Q WZ
obs
obs

 *


(
Z

Z
)

(
Q

Q
)
2
*
t
A x B x
AK
B
A
Q
Q Q
A 2 gAQ
obs

A



AW
(
Q

Q
)
Q
Q
2
t
A x
x
K
Variation of J with Respect to n
J (n)  
T
0

L
0
 f 2 gQ | Q | Q 
 
ndxdt
2
nK
 n

If the identified Manning’s roughness is a set of distributed parameter,
i.e., n = (n1, n2, …, nN)T, we have the chain law
 J
 J 
 J 
n2    
J   n1  
 n1 
 n2 
 nN
J
ni

T
0

xi 1
xi

nN

 f 2 gQ | Q | Q 
 
dxdt
2
nK
 x

where ni represent the Manning’s roughness in the reach Li [xi, xx+1]
Boundary Conditions - for Single Channel
Considering the contour integral in J*, This term I needs to be zero.
Q
Q 2Q
Q
Q
I   [(A  2 Q )A  Q ]dx  [ 3 A  (A  2 Q )Q ]dt
A



AB

BC


CD
A


A
0
DA
Final Conditions
Q ( x, T )  0, x  [0, L]
Backward Computation
t
A ( x, T )  0, x  [0, L]
Upstream
A
T
D
C
O
L
B
Q (0, t )  0, t [0, T ]
 A

Downstream Q ( L, t )  
 Q A ( L, t ), t  [0, T ]


2
A
x
Fig. Solution domain
Internal Boundary Conditions
– for Channel Network
I.B.C.s of Flow Model
Z1  Z 2  Z 3
Q3  Q1  Q2
I.B.C.s of Adjoint Equations
QQ  
QQ  
QQ 

 A 
   A 
   A 

2 
2 
2 
A 1 
A 2 
A 3

 B* Q 2Q   B* Q 2Q   B* Q 2Q 

 
 

3
3
3

 
 

A
A
A

3 
1 
2
Fig. Confluence
Numerical Techniques
1-D Time-Space Discretization (Preissmann, 1961)
 ( x, t )   [in11  (1  )in1 ]  (1   )[in1  (1  )in ]
 ( x, t )
in11  in1
in 1  in

 (1  )
t
t
t
 ( x, t )
in11  in 1
in1  in

 (1   )
x
x
x
where  and  are two weighting parameters in time and space, respectively;
t=time increment; x=spatial length
Solver of the resulting linear algebraic equations
(Pentadiagonal Matrix)
Double Sweep Algorithm based on the Gauss Elimination
Minimization Procedures
• Limited-Memory Quasi-Newton Method (LMQN)
Newton-like method, applicable for large-scale computation,
considering the second order derivative of objective function
(the approximate Hessian matrix) (Ding et al, 2002)
Algorithms:
BFGS (named after its inventors, Broyden, Fletcher,
Goldfarb, and Shanno)
L-BFGS (unconstrained optimization)
L-BFGS-B (bound constrained optimization)
Data Flows for Parameter Identification
Output Data
Input Data
Input data for the CCHE1D,
e.g., *.bc, *.bf
Observation data:
Filename: *.obs
Initial control data for
parameter : Filename: *.man
Control data of L-BFGS-B:
Filename: *.lbf
Output data from the
CCHE1D
Model of Parameter
Identification of
Manning Roughness in
the CCHE1D
Results of identified
parameter: Filename: *.para
iterate.dat
Results of Objective
Function: Filename: *.per
History output for each
observed point: *hispt.plt
Numerical Results (Verification)
• Approach for Verification
Specified n*
Flow Model
(CCHE1D)
Retrieving
Numerical
Results
n=n* ?
Identified n
Parameter
Identification
Observation
Data
Initial
Estimation of n0
Initial Estimations of Parameters : n0 << n*
Verify : identifiability, uniqueness of solution,
stability of identification process
Verification of Model ( Case 1)
– A Hypothetic Single Channel
Observation
Station*
Qp
1:2
Discharge
Channel
70m
Qb
+0.0m
Tp
Time
Td
Hydrograph
Parameter L
x
t


Unit
(km) (km) (min)
Value
10.0 0.5 10.0 1.0(0.55**) 0.5
*: Time interval of observation = 10t
**: This value is used for solving adjoint equations
1:1. 5
+2.0m
20m
Cross-section
QP
Qb
Tp
Td
(m3/s) (m3/s) (hour) (hour)
100.0 10.0 16.0
48.0
Verification of Model (Case 1)
– Single Parameter in A Hypothetic Single Channel
Table Control Parameters
in the L-BFGS-B
102
Manning's n
Relative error of n
nmax
1.0
nobj=
objective n in channel
n0 = initial estimation
nmin = lower bound
nmax = upper bound
m = store size in the LBFGS-B
Relative Error = |n-nobj|/n
10-1
10
1
10
0
10-1
10-2
10
Objective value :n = 0.03
-2
10
-3
10
-4
10-5
10-6
10
-3
10
20
Iterations of L-BFGS
30
10
-7
Relative error of n
100
Manning's n
nobj n0
nmin
0.03 0.01 0.001
m WQ W Z
7
1.0 1.0
History of Identified Manning's n in A Hypothetical Single Channel
Verification of Model (Case 1)
– Performance of Lagrangian Multipliers
Impulse response functions
Q(t)
0.0003
ITERATION= 2
A

Q

0
0
-1E-06
-0.0001
-2E-06
-0.0002
-3E-06
-0.0003
-4E-06
0
10
20
Time (hr)
30
40
50
2E-05
4E-07
1.5E-05
3E-07
ITERATION= 4
1E-05
A

Q

10
20
Time (hr)
30
40
50
10t
ITERATION= 4
2E-07
0
0
-5E-06
-1E-07
-1E-05
-2E-07
-1.5E-05
-3E-07
0
10
20
Time (hr)
30
40
50
-4E-07
0
10
20
Time (hr)
30
40
50
4E-07
2E-05
1.5E-05
ITERATION= 6
1E-05
ITERATION= 6
3E-07
2E-07
1E-07
A
5E-06
0

Q
0
1E-07
5E-06

ITERATION= 2
2E-06
1E-06
0.0001
0
-5E-06
-1E-07
-1E-05
-2E-07
-1.5E-05
-3E-07
-2E-05
Pulse Response
3E-06
0.0002
-2E-05
A(t)
4E-06
0
10
20
Time (hr) 30
40
50
-4E-07
0
10
20
Time (hr) 30
40
50
Iteration Processes of Lagrangian Multipliers Q and A at the Observation Station
Verification of Model (Case 2)
– Distributed Parameters in A Hypothetic Single Channel
Observation
Stations
Objective n
Initial Estimations:
Bound Constraints:
n0=(0.005, 0.005, 0.005)
0.001 n  1.0
Verification of Model (Case 2)
– Distributed Parameters in A Hypothetic Single Channel
History of identified parameters (CASE: Multi-parameters) by L-BFGS-B
100
n1
n2
n3
Relative Error of n
Manning's n
10-1
Objective Value
n3obj=0.03
100
n2obj=0.02
10
n1obj=0.01
-2
10-1
Relative Error of n
101
n  n obj
10
10
-3
0
10
20
30
40
Iterations of L-BFGS
50
60
70
-2
n obj
2
2
Verification of Model (Case 2)
– Distributed Parameters in A Hypothetic Single Channel
Iterative process of objective function and its gradient
(CASE: Multi-parameters)
Objective function
10-1
1
10
10-1
-3
-3
10
10
10-5
10-5
-7
10
-7
10
20
30
40
Iterations of L-BFGS
50
60
10
Norm of gradient
Objective function
Norm of gradient
1
10
J (n) 2
Verification of Model (Case 2)
– Distributed Parameters in A Hypothetic Single Channel
Simulation
Observation
60
40
4
3
2
Simulation
Observation
1
20
0
5
Water Stage (m)
80
3
Discharge (m /s)
100
0
10
20
30
Time (hr)
40
50
0
0
10
20
30
Time (hr)
40
Comparisons of Discharge and Stage at An Observation Station
50
Verification of Model (Case 3)
1
– Distributed Parameters in A Hypothetic Channel Network
3
L1 = 4,000
m
Channel
QP
Qb
No.
(m3/s) (m3/s)
1
50.0 2.0
2
50.0 2.0
3
60.0 6.0
L
3
=
13
Tp
(hour)
16.0
16.0
16.0
Td
(hour)
48.0
48.0
48.0
nobj=(0.01, 0.02, 0.03)
,0
00
m
n0=(0.005, 0.005, 0.005)
Channel N o . 2
L2 = 4,5 0 0 m
Observation Station
Confluence
0.001 n  1.0
A
Verification of Model (Case 3)
– Distributed Parameters in A Hypothetic Channel Network
History of identified parameters (CASE: Multi-parameters) by L-BFGS-B
100
n1
n2
n3
Relative Error of n
10-1
Manning's n
100
Objective Value
n3obj=0.03
n2obj=0.02
10
10
10
-1
10
-2
10
70
-3
n1obj=0.01
-2
-3
0
10
20
30
40
50
Iterations of L-BFGS-B
60
Relative Error of n
101
Verification of Model (Case 3)
– Distributed Parameters in A Hypothetic Channel Network
Iterative process of performance function and its gradient
(CASE 5 : Channel network)
3
3
10
10
Objective Function
Norm of Gradient
2
10
10
1
10
0
10
10-1
10-1
10-2
10-2
10-3
10-3
10-4
10-4
10-5
10-5
0
-6
10
-6
0
10
20
30
40
50
Iterations of Function Evaluation
60
10
70
Norm of Gradient
1
10
Objective Function
2
10
Verification of Model (Case 3)
– Distributed Parameters in A Hypothetic Channel Network
Searching Processes
Iterative process of hydrographs at station No.45 (St. A)
Iterative process of hydrographs at station No.45 (St. A)
ITERATION= 1
ITERATION= 2
ITERATION= 3
ITERATION= 4
ITERATION= 5
ITERATION= 6
ITERATION= 7
ITERATION= 20
ITERATION= 40
Observation at I = 45
140
130
120
3
Discharge (m /s)
110
100
90
Simulation: dt=10min
Observation: dt = 100min
80
70
60
50
ITERATION= 1
ITERATION= 2
ITERATION= 3
ITERATION= 4
ITERATION= 5
ITERATION= 6
ITERATION= 7
ITERATION= 20
ITERATION= 40
Observation at I = 45
9
8
7
6
Simulation: dt=10min
Observation: dt = 100min
5
4
3
40
2
30
20
1
10
0
10
Water Stage (m)
150
0
10
20
30
Time (hr)
40
50
0
0
10
20
30
40
50
Time (hr)
Comparisons of Discharge and Stage at An Observation Station (Downstream)
Application (Natural River)
– Distributed Parameters in East Fork River, Wyoming
OUTLET
St. A
INLET
Map of the Study Reach (3km), East Fork River, Wyomin
Application (East Fork River)
– Error Criterion and Control Parameters
1 T L
J (Q , Z , n )    ( Z  Z obs ) 2 dxdt ,
2 0 0
0.01  n  0.5
Total Number of Reaches: 83
Total Number of the Manning’s Roughness: N= 83
 Identify the profile of n along the river
The Size of Memory in the L-BFGS-BL: m = 10
CPU time in PC with 500MHz : nearly one hour
Application (East Fork River)
– Boundary Conditions at Inlet and Outlet
10
50
Discharge at Inlet
Stage at Outlet
First Flood:
Parameter Identification
The Other:
Prediction
45
40
7
35
6
30
5
25
4
20
3
15
2
10
1
0
5
For Identification
0
10
20
Days
0
30
3
8
Discharge (m /s)
Stage (m)
9
Application (East Fork River)
– Identification Processes of Profiles of n along the river
0
0.1
ITERATION= 1
ITERATION= 3
ITERATION= 6
ITERATION= 10
ITERATION= 20
Optimal Profile
0.09
0.08
Manning s n
0.07
0.06
0.05
0.04
0.03
0.02
nmin=0.01
0.01
0
0
1000
2000
Distance From Outlet (m)
3000
200m
Application (East Fork River)
– Objective Function and Norm of Its Gradient
10-1
10-1
Objective Function
Norm of Gradient
10-3
10
-4
10
-5
10-3
10
20
30
Iterations
40
50
10
-4
10
60
-5
Norm of Gradient
10-2
Objective Function
10-2
Application (East Fork River)
– Predictions of Stages
9.5
Observation at Inlet
Simulation at Inlet
Simulation at St. A
Observation at St. A
Water Stage (m)
9
8.5
8
7.5
Prediction
For Identification
7
0
10
Days
20
30
Conclusions
• An effective parameter identification tool has been
established, based on the adjoint sensitivity analysis for the
full nonlinear de Saint Venant Equations.
• The derived internal boundary conditions enable the tool to
identify the parameter in network channel.
• This tool can included engineers’ experience about the
parameter n in terms of bound constraints.
• As a result, combining with the CCHE1D model, this tool
can be utilized for identifications of n in natural channels
and branched network rivers.
• The adjoint equations can be directly applied to optimal
control for watershed management in the future.
Topics in the Future
• Apply to automatic calibrations in river flows;
• Enable the tool to handle other hydraulic structures, e.g.
culvert, bridge cross, …, so as to make the applications
widely.
• Further study the data requirements in the parameter
identification for possible research direction in optimally
selecting observation sites;
• Extend to establish optimal control procedures for flood
control and/or water quality management;
• Extend to the 2-D problems in flood waves and tidal flows
by using the CCHE2D
Variational Parameter Estimation in 2-D
Shallow Water Equations
Governing Equations:
Continuity:
 hU hU
L 


0
t
x
y
Momentum
hu huu huv
  xb
Lu 


 Fx  fhv  gh
 0
t
x
y
x 
hv huv hvv
  y
LV 


 Fy  fhu  gh
 0
t
x
y
y 
b
f: the Coriolis parameter
Eddy Viscosity and Bottom Friction

u   
u v 
Fx   2h t
   h t (  ) 
x 
x  y 
y x 
 
v   
u v 
Fx   2h t    h t (  ) 
y 
y  x 
y x 
 xb  Cb | u | u
 yb  Cb | u | v
Cb= bottom friction coefficient
Performance Function
1
J ( X )   ( X  X obs )T W ( X  X obs )dxdydt
2 xyt
X = physical variables, (u, v, η)
Variational Analysis
- To Obtain Adjoint Equations
Extended Objective Function
J * ( X , Cb )  J ( X )    L  u Lu  v Lv dxdydt
xyt
where A and Q are the Lagrangian multipliers
Necessary Condition
J *  J  0
on the conditions that
Lη = 0
Lu = 0
Lv = 0
Adjoint Equations
λη:
λu:

 hu hv 
obs
 g


W
(



)

t
y 
 x
     
u

   
  
 2hu u  hv  u  v    2h t u    h t ( u  v ) 
t
x
x  x 
x  y 
y
x 
 y

Cb  2u 2  v 2
2uv 2 
obs
 fhv   u
 v

h

W
(
u

u
)

 
|u |
|u | 
x
h
λv:
   
v

    
 
 
 hu  u  v   2hv v   h t ( u  v )    2h t v  
t
x 
y x 
y
x  y 
x 
 y

Cb  2uv 2
2u 2  v 2 
obs
 fhu   u
 v

h

W
(
v

v
)

  |u |
|u | 
y
h
Transversality (Final) Conditions of
Adjoint Variables
 ( x, y, t f )  0
u ( x, y, t f )  0
v ( x, y, t f )  0
Variation of J with Respect to n
 f

 J (Cb )   
 | u | (uu  vv )  Cb dxdydt
Cb

xyt 
If the identified bottom friction coefficient is a set of distributed
parameter, i.e., Cb = (Cb1, Cb2, …, CbN)T, we have the chain law
 J 
 J 
J 
  Cb1  
  Cb 2 
 Cb1 
 Cb 2 
J
C bi

tf
t0
 J

 CbN

  CbN

 f

  Cbi  | u | (uu  vv )  dxdydt
i
where Cbi represent the bottom friction over a subdomain Ωi Ω
References
[1] Yen, B. C. (1991), Hydraulic resistance in open channel, In: Channel
Flow Resistance: Centennial of Manning’s Formula, B.C. Yen, Ed., Water
Resources Publications, Highland Ranch, Colorado, 1-135.
[2] Yeh, W. W. –G. (1986), Review of parameter identification procedures
in groundwater hydrology: The inverse problem, Water Resour. Res.,
22(2), 95-108.
[3] Nocedal J. and S. J. Wright (1999), Numerical optimization, Springer
Verlag Series in Operation Research, Ed.: P. Glynn and S.M. Robinson,
Springer-Verlag.
[4] Sakawa Y. and Y. Shindo (1980), On global convergence of an
algorithm for optimal control, IEEE Transactions on Automatic Control,
vol. AC-25, No.6, 1149-1153.
References (cont.)
[5] Fletcher, R. (1987), Practical Methods of Optimization, John Wiley and
Sons, New York.
[6] Liu D.C. and J. Nocedal (1989), On the Limited Memory Method for
Large Scale Optimization, Mathematical Programming B, 45, 3, pp. 503528.
[7] Byrd, R.H., P. Lu and J. Nocedal (1995), A Limited Memory Algorithm
for Bound Constrained Optimization, SIAM Journal on Scientific and
Statistical Computing, 16, 5, pp.1190-1208.
[8] Gill P.E., W. Murray, and M. H. Wright (1981), Practical optimization,
Academic Press, pp346-353.
References (Cont.)
Byrd, R.H., P. Lu and J. Nocedal (1995), A Limited Memory Algorithm for Bound
Constrained Optimization, SIAM Journal on Scientific and Statistical Computing, 16, 5,
pp.1190-1208.
Ding, Y., Jia, Y., and Wang, S. S. Y., (2002), Identification of the Manning's roughness in
the CCHE2D model model with Sakawa-Shindo method and Limited-memory methods,
submitted to J. Hydr. Engrg..
Liu D.C. and J. Nocedal (1989), On the Limited Memory Method for Large Scale
Optimization, Mathematical Programming B, 45, 3, pp. 503-528.
Nocedal J. and S. J. Wright (1999), Numerical optimization, Springer Verlag Series in
Operation Reasearchn, Ed.: P. Glynn and S.M. Robinson, Springer-Verlag.
Vieira, D.~A., and Wu, W.-M., (2002), CCHE1D version 3.0-model capabilities and
applications, Technical Report No. NCCHE-TR-2002-05, NCCHE.
Yen, B. C. (1991), Hydraulic resistance in open channel, In: Channel Flow Resistance:
Centennial of Manning's Formula, B.C. Yen, Ed., Water Resources Publications,
Highland Ranch, Colorado, 1-135.