From: Jean-Philippe ARGAUD Date: Thu, 4 Feb 2021 18:29:13 +0000 (+0100) Subject: Improvement and extension of EnKF algorithm (EnKF) X-Git-Tag: V9_7_0b1~36 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=3757a00bd02f2e32f86201f5b2d53e55cda60c7e;p=modules%2Fadao.git Improvement and extension of EnKF algorithm (EnKF) --- diff --git a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py index a7e178e..a1d3598 100644 --- a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py +++ b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py @@ -34,7 +34,6 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): typecast = str, message = "Minimiseur utilisé", listval = [ - "StochasticEnKF", "EnKF", "ETKF", "ETKF-N", @@ -42,6 +41,9 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): "IEnKF", ], listadv = [ + "StochasticEnKF", + "EnKF-05", + "EnKF-16", "ETKF-KFF", "ETKF-VAR", "ETKF-N-11", @@ -170,9 +172,12 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): self._pre_run(Parameters, Xb, Y, U, HO, EM, CM, R, B, Q) # #-------------------------- - # Default EnKF = StochasticEnKF - if self._parameters["Minimizer"] in ["StochasticEnKF", "EnKF"]: - NumericObjects.senkf(self, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula") + # Default EnKF = EnKF-16 = StochasticEnKF + if self._parameters["Minimizer"] in ["EnKF-05"]: + NumericObjects.senkf(self, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula05") + # + elif self._parameters["Minimizer"] in ["EnKF-16", "StochasticEnKF", "EnKF"]: + NumericObjects.senkf(self, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16") # #-------------------------- # Default ETKF = ETKF-KFF diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py index 6ae0705..2c6b000 100644 --- a/src/daComposant/daCore/NumericObjects.py +++ b/src/daComposant/daCore/NumericObjects.py @@ -695,7 +695,7 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): else: Rn = R if hasattr(Q,"asfullmatrix"): Qn = Q.asfullmatrix(__n) else: Qn = Q - Xn = numpy.asmatrix(numpy.dot( Xb.reshape((__n,1)), numpy.ones((1,__m)) )) + Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m ) # if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]: selfA.StoredVariables["Analysis"].store( numpy.ravel(Xb) ) @@ -705,15 +705,11 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): # previousJMinimum = numpy.finfo(float).max # - # Predimensionnement - Xn_predicted = numpy.asmatrix(numpy.zeros((__n,__m))) - HX_predicted = numpy.asmatrix(numpy.zeros((__p,__m))) - # for step in range(duration-1): if hasattr(Y,"store"): - Ynpu = numpy.asmatrix(numpy.ravel( Y[step+1] )).T + Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,-1)) else: - Ynpu = numpy.asmatrix(numpy.ravel( Y )).T + Ynpu = numpy.ravel( Y ).reshape((__p,-1)) # if U is not None: if hasattr(U,"store") and len(U)>1: @@ -732,10 +728,11 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): ) # if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast - EMX = M( [(Xn[:,i], Un) for i in range(__m)], argsAsSerie = True ) - for i in range(__m): - qi = numpy.random.multivariate_normal(numpy.zeros(__n), Qn) - Xn_predicted[:,i] = (numpy.ravel( EMX[i] ) + qi).reshape((__n,-1)) + EMX = M( [(Xn[:,i], Un) for i in range(__m)], + argsAsSerie = True, + returnSerieAsArrayMatrix = True ) + qi = numpy.random.multivariate_normal(numpy.zeros(__n), Qn, size=__m).T + Xn_predicted = EMX + qi HX_predicted = H( [(Xn_predicted[:,i], Un) for i in range(__m)], argsAsSerie = True, returnSerieAsArrayMatrix = True ) @@ -750,25 +747,37 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): returnSerieAsArrayMatrix = True ) # # Mean of forecast and observation of forecast - Xfm = Xn_predicted.mean(axis=1, dtype=mfp).astype('float') - Hfm = HX_predicted.mean(axis=1, dtype=mfp).astype('float') + Xfm = Xn_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__n,-1)) + Hfm = HX_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1)) # #-------------------------- - if VariantM == "KalmanFilterFormula": + if VariantM == "KalmanFilterFormula05": PfHT, HPfHT = 0., 0. for i in range(__m): - Exfi = Xn_predicted[:,i] - Xfm.reshape((__n,-1)) - Eyfi = (HX_predicted[:,i] - Hfm).reshape((__p,1)) + Exfi = Xn_predicted[:,i] - Xfm + Eyfi = HX_predicted[:,i] - Hfm PfHT += Exfi * Eyfi.T HPfHT += Eyfi * Eyfi.T PfHT = (1./(__m-1)) * PfHT HPfHT = (1./(__m-1)) * HPfHT - K = PfHT * ( R + HPfHT ).I + Kn = PfHT * ( R + HPfHT ).I del PfHT, HPfHT # for i in range(__m): ri = numpy.random.multivariate_normal(numpy.zeros(__p), Rn) - Xn[:,i] = Xn_predicted[:,i] + K @ (numpy.ravel(Ynpu) + ri - HX_predicted[:,i]).reshape((__p,1)) + Xn[:,i] = Xn_predicted[:,i] + numpy.ravel(Kn @ (numpy.ravel(Ynpu) + ri - HX_predicted[:,i])).T + #-------------------------- + elif VariantM == "KalmanFilterFormula16": + EpY = EnsembleOfCenteredPerturbations(Ynpu, Rn, __m) + EpYm = EpY.mean(axis=1, dtype=mfp).astype('float').reshape((__p,-1)) + # + EaX = EnsembleOfAnomalies( Xn_predicted ) / numpy.sqrt(__m-1) + EaY = (HX_predicted - Hfm - EpY + EpYm) / numpy.sqrt(__m-1) + # + Kn = EaX @ EaX.T @ numpy.linalg.inv( EaY @ EaY.T) + # + for i in range(__m): + Xn[:,i] = Xn_predicted[:,i] + Kn @ (numpy.ravel(EpY[:,i]) - HX_predicted[:,i]) #-------------------------- else: raise ValueError("VariantM has to be chosen in the authorized methods list.") @@ -779,7 +788,7 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): selfA._parameters["InflationFactor"], ) # - Xa = Xn.mean(axis=1, dtype=mfp).astype('float') + Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,-1)) #-------------------------- # if selfA._parameters["StoreInternalVariables"] \