doc/en/ref_algorithm_3DVAR.rst

   1 ..
   2    Copyright (C) 2008-2015 EDF R&D
   3
   4    This file is part of SALOME ADAO module.
   5
   6    This library is free software; you can redistribute it and/or
   7    modify it under the terms of the GNU Lesser General Public
   8    License as published by the Free Software Foundation; either
   9    version 2.1 of the License, or (at your option) any later version.
  10
  11    This library is distributed in the hope that it will be useful,
  12    but WITHOUT ANY WARRANTY; without even the implied warranty of
  13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14    Lesser General Public License for more details.
  15
  16    You should have received a copy of the GNU Lesser General Public
  17    License along with this library; if not, write to the Free Software
  18    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19
  20    See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21
  22    Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  23
  24 .. index:: single: 3DVAR
  25 .. _section_ref_algorithm_3DVAR:
  26
  27 Calculation algorithm "*3DVAR*"
  28 -------------------------------
  29
  30 Description
  31 +++++++++++
  32
  33 This algorithm performs a state estimation by variational minimization of the
  34 classical :math:`J` function in static data assimilation:
  35
  36 .. math:: J(\mathbf{x})=(\mathbf{x}-\mathbf{x}^b)^T.\mathbf{B}^{-1}.(\mathbf{x}-\mathbf{x}^b)+(\mathbf{y}^o-\mathbf{H}.\mathbf{x})^T.\mathbf{R}^{-1}.(\mathbf{y}^o-\mathbf{H}.\mathbf{x})
  37
  38 which is usually designed as the "*3D-VAR*" function (see for example
  39 [Talagrand97]_).
  40
  41 Optional and required commands
  42 ++++++++++++++++++++++++++++++
  43
  44 .. index:: single: Background
  45 .. index:: single: BackgroundError
  46 .. index:: single: Observation
  47 .. index:: single: ObservationError
  48 .. index:: single: ObservationOperator
  49 .. index:: single: Minimizer
  50 .. index:: single: Bounds
  51 .. index:: single: MaximumNumberOfSteps
  52 .. index:: single: CostDecrementTolerance
  53 .. index:: single: ProjectedGradientTolerance
  54 .. index:: single: GradientNormTolerance
  55 .. index:: single: StoreSupplementaryCalculations
  56 .. index:: single: Quantiles
  57 .. index:: single: SetSeed
  58 .. index:: single: NumberOfSamplesForQuantiles
  59 .. index:: single: SimulationForQuantiles
  60
  61 The general required commands, available in the editing user interface, are the
  62 following:
  63
  64   Background
  65     *Required command*. This indicates the background or initial vector used,
  66     previously noted as :math:`\mathbf{x}^b`. Its value is defined as a
  67     "*Vector*" or a *VectorSerie*" type object.
  68
  69   BackgroundError
  70     *Required command*. This indicates the background error covariance matrix,
  71     previously noted as :math:`\mathbf{B}`. Its value is defined as a "*Matrix*"
  72     type object, a "*ScalarSparseMatrix*" type object, or a
  73     "*DiagonalSparseMatrix*" type object.
  74
  75   Observation
  76     *Required command*. This indicates the observation vector used for data
  77     assimilation or optimization, previously noted as :math:`\mathbf{y}^o`. It
  78     is defined as a "*Vector*" or a *VectorSerie* type object.
  79
  80   ObservationError
  81     *Required command*. This indicates the observation error covariance matrix,
  82     previously noted as :math:`\mathbf{R}`. It is defined as a "*Matrix*" type
  83     object, a "*ScalarSparseMatrix*" type object, or a "*DiagonalSparseMatrix*"
  84     type object.
  85
  86   ObservationOperator
  87     *Required command*. This indicates the observation operator, previously
  88     noted :math:`H`, which transforms the input parameters :math:`\mathbf{x}` to
  89     results :math:`\mathbf{y}` to be compared to observations
  90     :math:`\mathbf{y}^o`. Its value is defined as a "*Function*" type object or
  91     a "*Matrix*" type one. In the case of "*Function*" type, different
  92     functional forms can be used, as described in the section
  93     :ref:`section_ref_operator_requirements`. If there is some control :math:`U`
  94     included in the observation, the operator has to be applied to a pair
  95     :math:`(X,U)`.
  96
  97 The general optional commands, available in the editing user interface, are
  98 indicated in :ref:`section_ref_assimilation_keywords`. In particular, the
  99 optional command "*AlgorithmParameters*" allows to choose the specific options,
 100 described hereafter, of the algorithm. See
 101 :ref:`section_ref_options_AlgorithmParameters` for the good use of this command.
 102
 103 The options of the algorithm are the following:
 104
 105   Minimizer
 106     This key allows to choose the optimization minimizer. The default choice is
 107     "LBFGSB", and the possible ones are "LBFGSB" (nonlinear constrained
 108     minimizer, see [Byrd95]_, [Morales11]_ and [Zhu97]_), "TNC" (nonlinear
 109     constrained minimizer), "CG" (nonlinear unconstrained minimizer), "BFGS"
 110     (nonlinear unconstrained minimizer), "NCG" (Newton CG minimizer). It is
 111     strongly recommended to stay with the default.
 112
 113     Example : ``{"Minimizer":"LBFGSB"}``
 114
 115   Bounds
 116     This key allows to define upper and lower bounds for every state variable
 117     being optimized. Bounds have to be given by a list of list of pairs of
 118     lower/upper bounds for each variable, with possibly ``None`` every time
 119     there is no bound. The bounds can always be specified, but they are taken
 120     into account only by the constrained optimizers.
 121
 122     Example : ``{"Bounds":[[2.,5.],[1.e-2,10.],[-30.,None],[None,None]]}``
 123
 124   MaximumNumberOfSteps
 125     This key indicates the maximum number of iterations allowed for iterative
 126     optimization. The default is 15000, which is very similar to no limit on
 127     iterations. It is then recommended to adapt this parameter to the needs on
 128     real problems. For some optimizers, the effective stopping step can be
 129     slightly different of the limit due to algorithm internal control
 130     requirements.
 131
 132     Example : ``{"MaximumNumberOfSteps":100}``
 133
 134   CostDecrementTolerance
 135     This key indicates a limit value, leading to stop successfully the
 136     iterative optimization process when the cost function decreases less than
 137     this tolerance at the last step. The default is 1.e-7, and it is
 138     recommended to adapt it to the needs on real problems.
 139
 140     Example : ``{"CostDecrementTolerance":1.e-7}``
 141
 142   ProjectedGradientTolerance
 143     This key indicates a limit value, leading to stop successfully the iterative
 144     optimization process when all the components of the projected gradient are
 145     under this limit. It is only used for constrained optimizers. The default is
 146     -1, that is the internal default of each minimizer (generally 1.e-5), and it
 147     is not recommended to change it.
 148
 149     Example : ``{"ProjectedGradientTolerance":-1}``
 150
 151   GradientNormTolerance
 152     This key indicates a limit value, leading to stop successfully the
 153     iterative optimization process when the norm of the gradient is under this
 154     limit. It is only used for non-constrained optimizers.  The default is
 155     1.e-5 and it is not recommended to change it.
 156
 157     Example : ``{"GradientNormTolerance":1.e-5}``
 158
 159   StoreSupplementaryCalculations
 160     This list indicates the names of the supplementary variables that can be
 161     available at the end of the algorithm. It involves potentially costly
 162     calculations or memory consumptions. The default is a void list, none of
 163     these variables being calculated and stored by default. The possible names
 164     are in the following list: ["APosterioriCovariance", "BMA", "CostFunctionJ",
 165     "CurrentState", "OMA", "OMB", "Innovation", "SigmaObs2",
 166     "MahalanobisConsistency", "SimulatedObservationAtBackground",
 167     "SimulatedObservationAtCurrentState", "SimulatedObservationAtOptimum",
 168     "SimulationQuantiles"].
 169
 170     Example : ``{"StoreSupplementaryCalculations":["BMA","Innovation"]}``
 171
 172   Quantiles
 173     This list indicates the values of quantile, between 0 and 1, to be estimated
 174     by simulation around the optimal state. The sampling uses a multivariate
 175     gaussian random sampling, directed by the *a posteriori* covariance matrix.
 176     This option is useful only if the supplementary calculation
 177     "SimulationQuantiles" has been chosen. The default is a void list.
 178
 179     Example : ``{"Quantiles":[0.1,0.9]}``
 180
 181   SetSeed
 182     This key allow to give an integer in order to fix the seed of the random
 183     generator used to generate the ensemble. A convenient value is for example
 184     1000. By default, the seed is left uninitialized, and so use the default
 185     initialization from the computer.
 186
 187     Example : ``{"SetSeed":1000}``
 188
 189   NumberOfSamplesForQuantiles
 190     This key indicates the number of simulation to be done in order to estimate
 191     the quantiles. This option is useful only if the supplementary calculation
 192     "SimulationQuantiles" has been chosen. The default is 100, which is often
 193     sufficient for correct estimation of common quantiles at 5%, 10%, 90% or
 194     95%.
 195
 196     Example : ``{"NumberOfSamplesForQuantiles":100}``
 197
 198   SimulationForQuantiles
 199     This key indicates the type of simulation, linear (with the tangent
 200     observation operator applied to perturbation increments around the optimal
 201     state) or non-linear (with standard observation operator applied to
 202     perturbated states), one want to do for each perturbation. It changes mainly
 203     the time of each elementary calculation, usually longer in non-linear than
 204     in linear. This option is useful only if the supplementary calculation
 205     "SimulationQuantiles" has been chosen. The default value is "Linear", and
 206     the possible choices are "Linear" and "NonLinear".
 207
 208     Example : ``{"SimulationForQuantiles":"Linear"}``
 209
 210 Information and variables available at the end of the algorithm
 211 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 212
 213 At the output, after executing the algorithm, there are variables and
 214 information originating from the calculation. The description of
 215 :ref:`section_ref_output_variables` show the way to obtain them by the method
 216 named ``get`` of the variable "*ADD*" of the post-processing. The input
 217 variables, available to the user at the output in order to facilitate the
 218 writing of post-processing procedures, are described in the
 219 :ref:`subsection_r_o_v_Inventaire`.
 220
 221 The unconditional outputs of the algorithm are the following:
 222
 223   Analysis
 224     *List of vectors*. Each element is an optimal state :math:`\mathbf{x}*` in
 225     optimization or an analysis :math:`\mathbf{x}^a` in data assimilation.
 226
 227     Example : ``Xa = ADD.get("Analysis")[-1]``
 228
 229   CostFunctionJ
 230     *List of values*. Each element is a value of the error function :math:`J`.
 231
 232     Example : ``J = ADD.get("CostFunctionJ")[:]``
 233
 234   CostFunctionJb
 235     *List of values*. Each element is a value of the error function :math:`J^b`,
 236     that is of the background difference part.
 237
 238     Example : ``Jb = ADD.get("CostFunctionJb")[:]``
 239
 240   CostFunctionJo
 241     *List of values*. Each element is a value of the error function :math:`J^o`,
 242     that is of the observation difference part.
 243
 244     Example : ``Jo = ADD.get("CostFunctionJo")[:]``
 245
 246 The conditional outputs of the algorithm are the following:
 247
 248   APosterioriCovariance
 249     *List of matrices*. Each element is an *a posteriori* error covariance
 250     matrix :math:`\mathbf{A}*` of the optimal state.
 251
 252     Example : ``A = ADD.get("APosterioriCovariance")[-1]``
 253
 254   BMA
 255     *List of vectors*. Each element is a vector of difference between the
 256     background and the optimal state.
 257
 258     Example : ``bma = ADD.get("BMA")[-1]``
 259
 260   CurrentState
 261     *List of vectors*. Each element is a usual state vector used during the
 262     optimization algorithm procedure.
 263
 264     Example : ``Xs = ADD.get("CurrentState")[:]``
 265
 266   Innovation
 267     *List of vectors*. Each element is an innovation vector, which is in static
 268     the difference between the optimal and the background, and in dynamic the
 269     evolution increment.
 270
 271     Example : ``d = ADD.get("Innovation")[-1]``
 272
 273   MahalanobisConsistency
 274     *List of values*. Each element is a value of the Mahalanobis quality
 275     indicator.
 276
 277     Example : ``m = ADD.get("MahalanobisConsistency")[-1]``
 278
 279   OMA
 280     *List of vectors*. Each element is a vector of difference between the
 281     observation and the optimal state in the observation space.
 282
 283     Example : ``oma = ADD.get("OMA")[-1]``
 284
 285   OMB
 286     *List of vectors*. Each element is a vector of difference between the
 287     observation and the background state in the observation space.
 288
 289     Example : ``omb = ADD.get("OMB")[-1]``
 290
 291   SigmaObs2
 292     *List of values*. Each element is a value of the quality indicator
 293     :math:`(\sigma^o)^2` of the observation part.
 294
 295     Example : ``so2 = ADD.get("SigmaObs")[-1]``
 296
 297   SimulatedObservationAtBackground
 298     *List of vectors*. Each element is a vector of observation simulated from
 299     the background :math:`\mathbf{x}^b`.
 300
 301     Example : ``hxb = ADD.get("SimulatedObservationAtBackground")[-1]``
 302
 303   SimulatedObservationAtCurrentState
 304     *List of vectors*. Each element is an observed vector at the current state,
 305     that is, in the observation space.
 306
 307     Example : ``Ys = ADD.get("SimulatedObservationAtCurrentState")[-1]``
 308
 309   SimulatedObservationAtOptimum
 310     *List of vectors*. Each element is a vector of observation simulated from
 311     the analysis or optimal state :math:`\mathbf{x}^a`.
 312
 313     Example : ``hxa = ADD.get("SimulatedObservationAtOptimum")[-1]``
 314
 315   SimulationQuantiles
 316     *List of vectors*. Each element is a vector corresponding to the observed
 317     state which realize the required quantile, in the same order than the
 318     quantiles required by the user.
 319
 320     Example : ``sQuantiles = ADD.get("SimulationQuantiles")[:]``
 321
 322 See also
 323 ++++++++
 324
 325 References to other sections:
 326   - :ref:`section_ref_algorithm_Blue`
 327   - :ref:`section_ref_algorithm_ExtendedBlue`
 328   - :ref:`section_ref_algorithm_LinearityTest`
 329
 330 Bibliographical references:
 331   - [Byrd95]_
 332   - [Morales11]_
 333   - [Talagrand97]_