Merge pull request #25080 from lynnmunday/Tikhonov

add Tikhonov regularization to optimization reporter

modules/optimization/doc/content/source/optimizationreporters/OptimizationReporter.md

@@ -7,9 +7,14 @@
 This optimization reporter provides a basic interface for inverse optimization with measurement data where all data is provided in the input file or from csv file. The objective function is defined as:
 
 !equation
-f(\mathbf{u}, \mathbf{p}) = \frac{1}{2}\sum^N_i \left(u_i - \tilde u_i\right)^2,
+f(\mathbf{u}, \mathbf{p}) = \frac{1}{2}\sum^N_i \left(u_i - \tilde u_i\right)^2 +R(\mathbf{u}, \mathbf{p}),
 
-where \mathbf{u} represents the simulation solution and $u_i$ is the solution value at measurement location $i$. $\tilde u_i$ is the measurement value at location $i$. $\mathbf{p}$ is the vector of parameters being optimized that the simulation solution depends on.
+where $\mathbf{u}$ represents the simulation solution with $u_i$ being the solution value at measurement location $i$. $\tilde u_i$ is the measurement value at location $i$. $\mathbf{p}$ is the vector of parameters being optimized that the simulation solution depends on.  $R(\mathbf{u}, \mathbf{p})$ provides regularization of the objective function.  We currently only support the Tikhonov regularization given by
+
+!equation
+R(\mathbf{u}, \mathbf{p}) = \frac{\alpha}{2}\|\mathbf{p}\|^2,
+
+where $\alpha$ is the Tikhonov coefficient specified in the input file by [!param](/OptimizationReporter/OptimizationReporter/tikhonov_coeff).
 
 ## Measurement Data id=sec:measure_data
 
@@ -21,7 +26,7 @@ Additionally, locations and values can be specified at input using [!param](/Opt
 
 ## Optimization Parameters
 
-`OptimizationReporter` is also responsible for creating parameter vector(s) for optimization, setting the initial condition for the optimization, and setting parameter bounds. Although the [Optimize.md] executioner holds a single vector for parameter values, this vector can be split into groups of parameters. This is done by specifying a name for each group with [!param](/OptimizationReporter/OptimizationReporter/parameter_names) and the number of parameters in each group with [!param](/OptimizationReporter/OptimizationReporter/num_values). The total number of parameters is ultimately defined by the sum of [!param](/OptimizationReporter/OptimizationReporter/num_values). The initial condition for the optimization can then be defined with [!param](/OptimizationReporter/OptimizationReporter/initial_condition), where a vector of data must defined for each group.  This vector an be a single value in which case all parameters in that group are set to that value or a value can be set for every parameter in that group.  The lower and upper bounds for the parameters can then specified by [!param](/OptimizationReporter/OptimizationReporter/lower_bounds) and [!param](/OptimizationReporter/OptimizationReporter/upper_bounds), respectively. The bounds follow the same input format rules as the `initial_condtion`.  If no initial conditions are provided, the parameters are initialized with 0.  Default values for `upper_bounds` and `lower_bounds` are std::numeric<Real>::max() and std::numeric<Real>::lower(), respectively.  These bounds are only applied if a bounded optimization algorithm is used.
+`OptimizationReporter` is also responsible for creating parameter vector(s) for optimization, setting the initial condition for the optimization, and setting parameter bounds. Although the [Optimize.md] executioner holds a single vector for parameter values, this vector can be split into groups of parameters. This is done by specifying a name for each group with [!param](/OptimizationReporter/OptimizationReporter/parameter_names) and the number of parameters in each group with [!param](/OptimizationReporter/OptimizationReporter/num_values). The total number of parameters is ultimately defined by the sum of [!param](/OptimizationReporter/OptimizationReporter/num_values). The initial condition for the optimization can then be defined with [!param](/OptimizationReporter/OptimizationReporter/initial_condition), where a vector of data must defined for each group.  This vector an be a single value in which case all parameters in that group are set to that value or a value can be set for every parameter in that group.  The lower and upper bounds for the parameters can then specified by [!param](/OptimizationReporter/OptimizationReporter/lower_bounds) and [!param](/OptimizationReporter/OptimizationReporter/upper_bounds), respectively. The bounds follow the same input format rules as the `initial_condtion`.  If no initial conditions are provided, the parameters are initialized with 0.  Default values for `upper_bounds` and `lower_bounds` are `std::numeric<Real>::max()` and `std::numeric<Real>::lower()`, respectively.  These bounds are only applied if a bounded optimization algorithm is used.
 
 ## Declared Data
 
@@ -38,7 +43,7 @@ Additionally, locations and values can be specified at input using [!param](/Opt
 | Simulation values | `simulation_values` | $N$ |
 | $u_i - \tilde{u}_i$ | `misfit_values` | $N$ |
 | Values of parameter group $g$ | [!param](/OptimizationReporter/OptimizationReporter/parameter_names)$_g$ | [!param](/OptimizationReporter/OptimizationReporter/num_values)$_g$ |
-| Parameter Gradient | `adjoint` | $\sum_{g}$[!param](/OptimizationReporter/OptimizationReporter/num_values)$_g$ |
+| Parameter Gradient of parameter group $g$ | `grad_`[!param](/OptimizationReporter/OptimizationReporter/parameter_names)$_g$  | [!param](/OptimizationReporter/OptimizationReporter/num_values)$_g$ |
 
 
 ## Example Input Syntax

modules/optimization/include/executioners/OptimizeSolve.h

@@ -132,6 +132,7 @@ class OptimizeSolve : public SolveObject
     BOUNDED_CONJUGATE_GRADIENT,
     NEWTON_LINE_SEARCH,
     BOUNDED_NEWTON_LINE_SEARCH,
+    BOUNDED_QUASI_NEWTON_TRUST_REGION,
     NEWTON_TRUST_LINE,
     BOUNDED_NEWTON_TRUST_LINE,
     QUASI_NEWTON,

modules/optimization/include/optimizationreporters/OptimizationReporterBase.h

@@ -86,6 +86,9 @@ class OptimizationReporterBase : public OptimizationData
   /// Gradient values declared as reporter data
   std::vector<std::vector<Real> *> _gradients;
 
+  /// Tikhonov Coefficient for regularization
+  const Real _tikhonov_coeff;
+
   /// Bounds of the parameters
   std::vector<Real> _lower_bounds;
   std::vector<Real> _upper_bounds;

modules/optimization/src/executioners/OptimizeSolve.C

@@ -18,8 +18,9 @@ InputParameters
 OptimizeSolve::validParams()
 {
   InputParameters params = emptyInputParameters();
-  MooseEnum tao_solver_enum("taontr taobntr taobncg taonls taobnls taontl taobntl taolmvm "
-                            "taoblmvm taonm taobqnls taoowlqn taogpcg taobmrm");
+  MooseEnum tao_solver_enum(
+      "taontr taobntr taobncg taonls taobnls taobqnktr taontl taobntl taolmvm "
+      "taoblmvm taonm taobqnls taoowlqn taogpcg taobmrm");
   params.addRequiredParam<MooseEnum>(
       "tao_solver", tao_solver_enum, "Tao solver to use for optimization.");
   ExecFlagEnum exec_enum = ExecFlagEnum();
@@ -95,6 +96,9 @@ OptimizeSolve::taoSolve()
     case TaoSolverEnum::BOUNDED_NEWTON_LINE_SEARCH:
       ierr = TaoSetType(_tao, TAOBNLS);
       break;
+    case TaoSolverEnum::BOUNDED_QUASI_NEWTON_TRUST_REGION:
+      ierr = TaoSetType(_tao, TAOBQNKTR);
+      break;
     case TaoSolverEnum::NEWTON_TRUST_LINE:
       ierr = TaoSetType(_tao, TAONTL);
       break;
@@ -348,7 +352,6 @@ OptimizeSolve::objectiveFunction()
     _inner_solve->solve();
 
   _obj_iterate++;
-
   return _obj_function->computeObjective();
 }
 

modules/optimization/src/optimizationreporters/OptimizationReporterBase.C

@@ -18,6 +18,8 @@ OptimizationReporterBase::validParams()
   params.registerBase("OptimizationReporterBase");
   params.addRequiredParam<std::vector<ReporterValueName>>(
       "parameter_names", "List of parameter names, one for each group of parameters.");
+  params.addRangeCheckedParam<Real>(
+      "tikhonov_coeff", 0.0, "tikhonov_coeff >= 0", "Coefficient for Tikhonov Regularization.");
   params.suppressParameter<std::vector<VariableName>>("variable");
   params.suppressParameter<std::vector<std::string>>("variable_weight_names");
   params.registerBase("OptimizationReporterBase");
@@ -29,7 +31,8 @@ OptimizationReporterBase::OptimizationReporterBase(const InputParameters & param
     _parameter_names(getParam<std::vector<ReporterValueName>>("parameter_names")),
     _nparams(_parameter_names.size()),
     _parameters(_nparams),
-    _gradients(_nparams)
+    _gradients(_nparams),
+    _tikhonov_coeff(getParam<Real>("tikhonov_coeff"))
 {
   for (const auto & i : make_range(_nparams))
   {
@@ -51,6 +54,16 @@ OptimizationReporterBase::computeObjective()
   for (auto & misfit : _misfit_values)
     val += misfit * misfit;
 
+  if (_tikhonov_coeff > 0.0)
+  {
+    Real param_norm_sqr = 0;
+    for (const auto & data : _parameters)
+      for (const auto & val : *data)
+        param_norm_sqr += val * val;
+
+    val += _tikhonov_coeff * param_norm_sqr;
+  }
+
   return val * 0.5;
 }
 
@@ -71,15 +84,26 @@ OptimizationReporterBase::setSimulationValuesForTesting(std::vector<Real> & data
 void
 OptimizationReporterBase::computeGradient(libMesh::PetscVector<Number> & gradient) const
 {
-  for (const auto & p : make_range(_nparams))
-    if (_gradients[p]->size() != _nvalues[p])
+  for (const auto & param_group_id : make_range(_nparams))
+  {
+    if (_gradients[param_group_id]->size() != _nvalues[param_group_id])
       mooseError("The gradient for parameter ",
-                 _parameter_names[p],
+                 _parameter_names[param_group_id],
                  " has changed, expected ",
-                 _nvalues[p],
+                 _nvalues[param_group_id],
                  " versus ",
-                 _gradients[p]->size(),
+                 _gradients[param_group_id]->size(),
                  ".");
+
+    if (_tikhonov_coeff > 0.0)
+    {
+      auto params = _parameters[param_group_id];
+      auto grads = _gradients[param_group_id];
+      for (const auto & param_id : make_range(_nvalues[param_group_id]))
+        (*grads)[param_id] += (*params)[param_id] * _tikhonov_coeff;
+    }
+  }
+
   OptUtils::copyReporterIntoPetscVector(_gradients, gradient);
 }
 

modules/optimization/test/tests/optimizationreporter/bimaterial/gold/main_out_reg_OptimizationReporter_0001.csv

@@ -0,0 +1,12 @@
+diffusivity_values,grad_diffusivity_values,measurement_time,measurement_values,measurement_xcoord,measurement_ycoord,measurement_zcoord,misfit_values,simulation_values
+5.0392293470097,-1.8258437961549e-07,0,1.5499999950466,0,-3.99,0,-0.012066421017338,1.5379335740293
+9.1392923375053,9.2441609922389e-11,0,122.15999960961,0,-3.192,0,-0.95098967192064,121.20900993769
+0,0,0,230.74999926258,0,-2.394,0,-1.7963397740259,228.95365948855
+0,0,0,326.45999895672,0,-1.596,0,-2.54142180988,323.91857714684
+0,0,0,409.18999869233,0,-0.798,0,-3.1854573006883,406.00454139164
+0,0,0,479.99999846604,0,0,0,-3.736698121486,476.26330034455
+0,0,0,508.43499922095,0,0.798,0,-1.058786683721,507.37621253723
+0,0,0,530.90999981763,0,1.596,0,1.0578322564506,531.96783207408
+0,0,0,546.895000242,0,2.394,0,2.5632451044535,549.45824534645
+0,0,0,556.44000049541,0,3.192,0,3.4621606899426,559.90216118535
+0,0,0,559.97500058926,0,3.99,0,3.7950749481626,563.77007553742

modules/optimization/test/tests/optimizationreporter/bimaterial/main.i

@@ -25,9 +25,9 @@
   type = OptimizationReporter
   parameter_names = diffusivity_values
   num_values = 2 # diffusivity in the bottom material and in the top material of model.i
-  initial_condition = '3 4' # the expected result is about '1 10' so this initial condition is not too bad
+  initial_condition = '15 15' # the expected result is about '1 10' so this initial condition is not too bad
   lower_bounds = '1'
-  upper_bounds = '20'
+  upper_bounds = '50'
   measurement_file = 'synthetic_data.csv'
   file_value = 'temperature'
 []
@@ -38,7 +38,7 @@
   # petsc_options_iname = '-tao_fd_gradient -tao_gatol'
   # petsc_options_value = ' true            0.001'
   type = Optimize
-  tao_solver = taoblmvm
+  tao_solver = taobqnktr
   petsc_options_iname = '-tao_gatol'
   petsc_options_value = '1e-3'
   ## THESE OPTIONS ARE FOR TESTING THE ADJOINT GRADIENT

modules/optimization/test/tests/optimizationreporter/bimaterial/regularization_for_main.i

@@ -0,0 +1,16 @@
+[OptimizationReporter]
+  type = OptimizationReporter
+  parameter_names = diffusivity_values
+  num_values = 2
+  initial_condition = '35 35'
+  tikhonov_coeff = 10
+  lower_bounds = '1'
+  upper_bounds = '50'
+  measurement_file = 'synthetic_data.csv'
+  file_value = 'temperature'
+[]
+
+[Outputs]
+  file_base = main_out_reg
+  csv = true
+[]

modules/optimization/test/tests/optimizationreporter/bimaterial/tests

@@ -2,23 +2,34 @@
   issues = '#21885'
   design = 'Optimize.md NearestReporterCoordinatesFunction.md ElementOptimizationDiffusionCoefFunctionInnerProduct.md'
   [bimaterial]
-    requirement = 'The system shall be able to invert for the value of a material property in two separate regions using'
-    [taoblmvm]
+    requirement = 'The system shall be able to invert for the value of a material property in two separate regions using '
+    [adjoint]
       type = CSVDiff
       max_threads = 1 # Optimize executioner does not support multiple threads
       input = main.i
-      abs_zero = 1e-2
+      abs_zero = 1 # test parameter output; ignore gradient and objectve
       csvdiff = main_out_OptimizationReporter_0001.csv
       recover = false
-      detail = 'a gradient evaluated in a separate multiapp or'
+      detail = 'a gradient evaluated in a separate multiapp with parameter initial conditons 1.5x higher than the true solution; '
+    []
+    [adjoint_reg]
+      type = CSVDiff
+      input = "main.i regularization_for_main.i"
+      csvdiff = main_out_reg_OptimizationReporter_0001.csv
+      max_threads = 1 # Optimize executioner does not support multiple threads
+      rel_err = 0.1
+      abs_zero = 1 # test parameter output; ignore gradient and objectve
+      # steady solve
+      recover = false
+      detail = 'a gradient evaluated in a separate multiapp using Tikhonov regularization and parameter initial conditons 3x higher than the true solution; or '
     []
     [auto_adjoint]
       type = CSVDiff
       input = main_auto_adjoint.i
       cli_args = 'Outputs/file_base=main_out'
       csvdiff = main_out_OptimizationReporter_0001.csv
       max_threads = 1 # Optimize executioner does not support multiple threads
-      abs_zero = 1e-2
+      abs_zero = 1 # test parameter output; ignore gradient and objectve
       recover = false
       detail = 'a gradient evaluated with an automatically computed adjoint.'
     []

modules/optimization/test/tests/optimizationreporter/point_loads/gold/main_out_reg_OptimizationReporter_0001.csv

@@ -0,0 +1,4 @@
+grad_parameter_results,measurement_time,measurement_values,measurement_xcoord,measurement_ycoord,measurement_zcoord,misfit_values,parameter_results,simulation_values
+-1.8041124150159e-16,0,293,0.5,0.28,0,2.1564416769367,-235.62098191604,295.15644167694
+-1.0026701691146e-15,0,320,0.5,1.1,0,-1.5143353633924,-257.85338047184,318.48566463661
+-5.1764148523148e-15,0,0,0,0,0,0,512.22832227015,0

modules/optimization/test/tests/optimizationreporter/point_loads/regularization_for_main.i

@@ -0,0 +1,14 @@
+[OptimizationReporter]
+  type = OptimizationReporter
+  parameter_names = 'parameter_results'
+  num_values = '3'
+  tikhonov_coeff = 0.0001
+  measurement_points = '0.5 0.28 0
+                        0.5 1.1 0'
+  measurement_values = '293 320'
+[]
+
+[Outputs]
+  file_base = main_out_reg
+  csv = true
+[]

modules/optimization/test/tests/optimizationreporter/point_loads/tests

@@ -1,6 +1,6 @@
 [Tests]
-  issues = '#21885'
-  design = 'ReporterPointSource.md'
+  issues = '#21885 #24101'
+  design = 'ReporterPointSource.md OptimizationReporter/index.md'
   [point_loads]
     requirement = "The system shall be able to invert for point loads using"
     [hessian]
@@ -17,9 +17,20 @@
       recover = false
       detail = 'Hessian-based optimization,'
     []
+    [hessian_reg]
+      type = CSVDiff
+      input = "main.i regularization_for_main.i"
+      csvdiff = main_out_reg_OptimizationReporter_0001.csv
+      max_threads = 1 # Optimize executioner does not support multiple threads
+      rel_err = 0.1
+      abs_zero = 10
+      # steady solve
+      recover = false
+      detail = 'ill-posed Hessian-based optimization with Tikhonov regularization,'
+    []
     [adjoint]
       type = CSVDiff
-      input = main.i
+      input =main.i
       cli_args = "Executioner/tao_solver=taolmvm "
                  "Executioner/petsc_options_iname='-tao_gttol' "
                  "Executioner/petsc_options_value='1e-5'"
@@ -29,7 +40,21 @@
       abs_zero = 10
       # steady solve
       recover = false
-      detail = "gradient-based optimization, or"
+      detail = 'gradient-based optimization, '
+    []
+    [adjoint_reg]
+      type = CSVDiff
+      input = "main.i regularization_for_main.i"
+      cli_args = "Executioner/tao_solver=taolmvm "
+                 "Executioner/petsc_options_iname='-tao_gttol' "
+                 "Executioner/petsc_options_value='1e-5'"
+      csvdiff = main_out_reg_OptimizationReporter_0001.csv
+      max_threads = 1 # Optimize executioner does not support multiple threads
+      rel_err = 0.1
+      abs_zero = 10
+      # steady solve
+      recover = false
+      detail = 'ill-posed gradient-based optimization with Tikhonov regularization, '
     []
     [auto_adjoint]
       type = CSVDiff
@@ -41,7 +66,19 @@
       abs_zero = 10
       # steady solve
       recover = false
-      detail = 'gradient-based optimization with an automatically computed adjoint.'
+      detail = 'gradient-based optimization with an automatically computed adjoint, or '
+    []
+    [auto_adjoint_reg]
+      type = CSVDiff
+      cli_args = "OptimizationReporter/initial_condition='1000 1000 1000'"
+      input = 'main_auto_adjoint.i regularization_for_main.i'
+      csvdiff = main_out_reg_OptimizationReporter_0001.csv
+      max_threads = 1 # Optimize executioner does not support multiple threads
+      rel_err = 0.1
+      abs_zero = 10
+      # steady solve
+      recover = false
+      detail = 'gradient-based optimization with an automatically computed adjoint, with Tikhonov regularization, and initial conditions do not affect the converged regularized solution.'
     []
   []
 []