Implementation of Sequential Uncertainty Fitting (SUFI-2) method using PEST/PEST++ utilities?

[sjepsen395]

Hi All,

I'm interested in trying out the Sequential Uncertainty Fitting (SUFI-2) method (Abbaspour et al., 2007) for our highly non-linear (swiss cheese) model calibration problem. SUFI-2 is a semi-automated method combining Latin Hypercube Sampling (LHS) and a "sensitivity matrix" governing the contraction of parameter ranges (more influence = more contraction). Perhaps many of you are familiar with it. I've had success with this method in the past for a surface water / GW (bucket-type) model, though that model was much simpler than the one I'm currently using. My impression was that SUFI-2 did well when initial parameter values are not well known and the objective function surface has many local minima (swiss cheese effect). Last I heard (several years ago), SUFI-2 was proprietary.

Abbaspour, K. C., Yang, J., Maximov, I., Siber, R., Bogner, K., Mieleitner, J., Zobrist, J., & Srinivasan, R. (2007). Modelling hydrology and water quality in the pre-alpine/alpine Thur watershed using SWAT. Journal of Hydrology, 333(2–4), 413–430. SWAT. https://doi.org/10.1016/j.jhydrol.2006.09.014

I wonder if one could effectively emulate SUFI-2 using PEST++/PEST utilities? For example, something like the following??:

(1) Obtain parameter samples (sets) using a PEST LHS utility
(2) Run PEST once on each parameter set (NOPTMAX=0)
(3) Identify best parameters from Step 3 and call these "preferred parameters"
(4) PARREP preferred parameters to Pest Control File
(5) Run PEST with NOPTMAX=-1 on preferred parameters to generate a Jacobian matrix
(6) Run PEST utility GENLINPRED to compute "post-calibration" parameter uncertainties
(7) Assign ranges to preferred parameters using results from Step 6
(8) Repeat steps

Does anyone see any theoretical issues with this implementation? Perhaps there is a better way using PEST++/PEST utilities (e.g., Iterative Ensemble Smoother applied in steps?). There is a "sensitivity matrix" in Abbaspour et al. (2007), an observation-parameter covariance matrix. Maybe there is a way to use this sensitivity matrix instead of the Jacobian matrix, which would bring the method closer to SUFI-2. The sensitivity matrix reminds me of the covariance matrix of the Iterative Ensemble Smoother method in that both involve d(response)/d(parameter) in a randomized sense.

Thanks,
Steve

[jtwhite79]

That makes sense to me. I think you could take advantage of pestpp-swp to parallelize the evaluation of the parameter samples and also the baked-in FOSM stuff in pestpp-glm to get the par uncertainty piece (similar to with genlinpred would do). So the steps might be:

0. LHS sampling
1. evaluate the LHS samples in parallel using pestpp-swp (which also calculates phi for you)
2. parep best pars
3. use pestpp-glm with noptmax=-1 to get the jco matrix (in parallel) and the par unc estimates
   4)some crafty action to update what to do next time around

If your model is highly nonlinear, you could turn to pestpp-mou, which is a global evolutionary optimizer with the likes of differential evolution and particle swarm optimization to evolve the population...and you set it to use a single objective, with that objective being some measure of observed vs simulated that you want to minimize (can be any measure you like, just needs to be an "observation" in the control file). currently pestpp-mou only has a uniform random sampler internally to initialize the population, but you can also feed it an external population like one from LHS sampling...

Let us know how it goes If you go down the SUFI2 path...

[sjepsen395]

Hi jtwhite79,

I'm interested in trying out particle swarm optimization since success of our calibration seems so dependent on initial parameters.

Do you think it would be OK to use pestpp-pso for this instead of the pestpp-mou that you mentioned above? I'm having an easier time following the documentation for pestpp-pso, and pestpp-pso sounds simpler to me than pestpp-mou. My objective function is just the sum of weighted residuals squared.

Also, do you know of an example set of input files out there for pestpp-pso? In particular, I'm having trouble understanding how to specify my objective function to be the sum of weighted residuals squared of all observations in the PEST Control File. Page 174 of the PEST++ user's manual says "This name [parameter OBJNME] must coincide with an observation group name specified in the PEST control file." Does that mean that the objective function will only include those observations whose group name matches the OBJNME string? If that is the case, would one need to rename all of their observation groups to be the same name, set to OBJNME, if they wanted to include all observations in the objective function?

Thanks,
Steve

[jtwhite79]

Im biased but I dont think mou isnt so hard to use, you just to be have one of the observations be the sum-of-weighted-residual (the standard objective function we use in pest) and after that it should "just work" with like one or two optional ++ args (like ++mou_objectives and ++mou_population_size) - PSO is now the default generator in mou...and with mou, you arent limited to just the standard sum-of-weighted-residuals objective function - you can use any measure of misfit you want. In fact, you could use a couple (like percent bias and MAE) and look at the trade off between them - you just need to specify the observations you want to use as objectives.

If you do want to use pestpp-pso, @siadicus is the expert to talk to (and is also the expert with PSO in general)...

But beware! these global methods do not scale to large numbers of parameters/decision variables. Like a hundred is pushing it. And you can expect to need lots of iterations/generations (like 20-30 or more). So computational demand is high...

[sjepsen395]

jtwhite79,

OK, I'll give pestpp-mou a try.

Sorry, I'm missing something. Say I only had 4 observations in my PEST Control File:

• observation data
  wsb01_0001_101983 4.5583000000E+02 1.0430719108E+00 wg_01
  wsb01_0001_031984 4.5553000000E+02 1.0430719108E+00 wg_01
  wsb02_0001_052002 1.6687000000E+02 1.0430719108E+00 wg_02
  wsb02_0001_062002 1.5724500000E+02 1.0430719108E+00 wg_02

What would I enter into the PEST Control File to set an observation to the sum of weighted residuals squared? Would it be possible to see an example? I don't know the modeled counterpart to observations before I run PEST, so I don't have the information to enter that into my sum of weighted residuals squared. Each observation has a weight, which I'd set to 1.0 for the objective? Obviously I'm missing something.

That makes sense about computational demand. I have about 120 parameters, so maybe I can just squeak by.

Thanks,
Steve

[jtwhite79]

Ok. So what I was thinking is that you could have a post-processing script that runs right after your model as part of the "forward run" process that mou is calling for each model evaluation and that this post-processor would calculate one or more metrics that summarize the misfit between observed and what was just simulated (could be the standard pest "phi" or others like MAE, pbias, etc). And the those metrics would added to the control file as additional "observations". In this way, mou would see the resulting metric(s) for each model run.

Then if you set the obs group name of the metric to start with "less_than" and pass the observation name of the metric to the mou_objectives arg, this will tell mou you want to minimize that quantity (assuming making the metric a smaller number is equivalent to improving the fit between observed and simulated). clear as mud?

[sjepsen395]

jtwhite79,

I follow the concept now. I'll weigh this approach in with other options on the table. Doing this scripting to compute an objective may be a little too much for me given that I don't know how well PSO will work for our model calibration. I'll mull over this.

Thanks,
Steve