From aa4c95594b3376b753c0259d396ba881a1cfa8bc Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Thu, 14 Jan 2021 08:02:02 +0100 Subject: [PATCH] Improvement and extension of EnKF algorithm (MLEF-T) ------------- --- src/daComposant/daAlgorithms/EnsembleKalmanFilter.py | 4 ++++ src/daComposant/daCore/NumericObjects.py | 4 ++++ 2 files changed, 8 insertions(+) diff --git a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py index 313ea75..8a56f1f 100644 --- a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py +++ b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py @@ -44,6 +44,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): "ETKF-N-16", "MLEF", "MLEF-B", + "MLEF-T", ], ) self.defineRequiredParameter( @@ -189,6 +190,9 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): elif self._parameters["Minimizer"] in ["MLEF-B", "MLEF"]: NumericObjects.mlef(self, Xb, Y, U, HO, EM, CM, R, B, Q, BnotT=False) # + elif self._parameters["Minimizer"] == "MLEF-T": + NumericObjects.mlef(self, Xb, Y, U, HO, EM, CM, R, B, Q, BnotT=True) + # #-------------------------- else: raise ValueError("Error in Minimizer name: %s"%self._parameters["Minimizer"]) diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py index 9d6ed56..7103817 100644 --- a/src/daComposant/daCore/NumericObjects.py +++ b/src/daComposant/daCore/NumericObjects.py @@ -1262,6 +1262,8 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, BnotT=False, _epsilon=1.e-1, _e=1 # if BnotT: EE = vx.reshape((__n,-1)) + _epsilon * EaX # 7: + else: + EE = vx.reshape((__n,-1)) + numpy.sqrt(__m-1) * EaX @ Ta # 8: # EZ = H( [(EE[:,i], Un) for i in range(__m)], argsAsSerie = True, @@ -1271,6 +1273,8 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, BnotT=False, _epsilon=1.e-1, _e=1 # if BnotT: EY = (EZ - ybar) / _epsilon # 11: + else: + EY = ( (EZ - ybar) @ numpy.linalg.inv(Ta) ) / numpy.sqrt(__m-1) # 12: # GradJ = numpy.ravel(vw.reshape((__m,1)) - EY.transpose() @ (RI * (Ynpu - ybar))) # 13: mH = numpy.eye(__m) + EY.transpose() @ (RI * EY) # 14: -- 2.39.2