import numpy as np

Nofrealization = 3 # Number of Phi Nofobservedwell = 1 # Number of obervation (In that case, Q_0) Observation = np.array([10]) # observation data, z phi = np.array([ [0.7, 2.2, 2.0, 0.4, 13], [0.6, 1.8, 2.4, 0.6, 8], [1.1, 1.7, 1.9, 0.5, 12] ]).transpose() En_mean = np.mean(phi,1)[np.newaxis,:].reshape(-1,1)

def Ensemble_Kalman_Filter(Observation, State_vector,En_mean, Nofrealization, Nofobservedwell): # 측정행렬 H = np.zeros([Nofobservedwell ,State_vector.shape[0]]) H for i in range(Nofobservedwell): H[-1 - i,-1 - i] = 1 # 추정오차공분산 L = State_vector - ( En_mean * np.ones([1,Nofrealization])) Cov_e = 1 / (Nofrealization-1) * L @ L.transpose() # 노이즈 추가 Observation = Observation.reshape(Nofobservedwell,-1) * np.ones([Nofobservedwell, Nofrealization]) + np.array([[-0.4, 0.5, 0.1]]).reshape([Nofobservedwell, Nofrealization]) C_d = pow(1,2) ''' Observation = Observation.reshape(Nofobservedwell,-1) * np.ones([Nofobservedwell, Nofrealization]) + pow(StdErr_Dynamic,2) * np.random.rand(Nofobservedwell, Nofrealization) C_d = np.zeros([Nofobservedwell,Nofobservedwell]) for kk in range(Nofobservedwell): C_d[kk, kk] = pow(StdErr_Dynamic,2) ''' # 칼만 게인 K = Cov_e @ H.transpose() @ np.linalg.inv(H @ Cov_e @ H.transpose() + C_d) # 교정 메인수식 State_vector_Assimilation = State_vector + K @ (Observation - H @ State_vector) return State_vector_Assimilation

phi_Assimilation = Ensemble_Kalman_Filter(Observation, phi, En_mean, Nofrealization, Nofobservedwell) print(phi_Assimilation)

phi_Assimilation[-1] = np.array([11, 10, 12]) phi_next = phi_Assimilation phi_next_mean = np.mean(phi_next,1)[:,np.newaxis].reshape([phi_next.shape[0],-1]) print('Phi for Next Step\n',phi_next) print('*'*40) print('Mean of Phi for Next\n', phi_next_mean)

phi_Assimilation_next = Ensemble_Kalman_Filter(Observation, phi_next, phi_next_mean, Nofrealization, Nofobservedwell) print(phi_Assimilation_next)