fmda/moisture_models.py

   1 import numpy as np
   2 import math
   3 import matplotlib.pyplot as plt
   4 from utils import vprint
   5
   6
   7 def model_decay(m0,E,partials=0,T1=0.1,tlen=1):
   8     # Arguments:
   9     #   m0          fuel moisture content at start dimensionless, unit (1)
  10     #   E           fuel moisture eqilibrium (1)
  11     #   partials=0: return m1 = fuel moisture contents after time tlen (1)
  12     #           =1: return m1, dm0/dm0
  13     #           =2: return m1, dm1/dm0, dm1/dE
  14     #           =3: return m1, dm1/dm0, dm1/dE dm1/dT1
  15     #   T1          1/T, where T is the time constant approaching the equilibrium
  16     #               default 0.1/hour
  17     #   tlen        the time interval length, default 1 hour
  18
  19     exp_t = np.exp(-tlen*T1)                  # compute this subexpression only once
  20     m1 = E + (m0 - E)*exp_t                   # the solution at end
  21     if partials==0:
  22         return m1
  23     dm1_dm0 = exp_t
  24     if partials==1:
  25         return m1, dm1_dm0          # return value and Jacobian
  26     dm1_dE = 1 - exp_t
  27     if partials==2:
  28         return m1, dm1_dm0, dm1_dE
  29     dm1_dT1 = -(m0 - E)*tlen*exp_t            # partial derivative dm1 / dT1
  30     if partials==3:
  31         return m1, dm1_dm0, dm1_dE, dm1_dT1       # return value and all partial derivatives wrt m1 and parameters
  32     raise('Bad arg partials')
  33
  34
  35 def ext_kf(u,P,F,Q=0,d=None,H=None,R=None):
  36     """
  37     One step of the extended Kalman filter.
  38     If there is no data, only advance in time.
  39     :param u:   the state vector, shape n
  40     :param P:   the state covariance, shape (n,n)
  41     :param F:   the model function, args vector u, returns F(u) and Jacobian J(u)
  42     :param Q:   the process model noise covariance, shape (n,n)
  43     :param d:   data vector, shape (m). If none, only advance in time
  44     :param H:   observation matrix, shape (m,n)
  45     :param R:   data error covariance, shape (n,n)
  46     :return ua: the analysis state vector, shape (n)
  47     :return Pa: the analysis covariance matrix, shape (n,n)
  48     """
  49     def d2(a):
  50         return np.atleast_2d(a) # convert to at least 2d array
  51
  52     def d1(a):
  53         return np.atleast_1d(a) # convert to at least 1d array
  54
  55     # forecast
  56     uf, J  = F(u)          # advance the model state in time and get the Jacobian
  57     uf = d1(uf)            # if scalar, make state a 1D array
  58     J = d2(J)              # if scalar, make jacobian a 2D array
  59     P = d2(P)              # if scalar, make Jacobian as 2D array
  60     Pf  = d2(J.T @ P) @ J + Q  # advance the state covariance Pf = J' * P * J + Q
  61     # analysis
  62     if d is None or not d.size :  # no data, no analysis
  63         return uf, Pf
  64     # K = P H' * inverse(H * P * H' + R) = (inverse(H * P * H' + R)*(H P))'
  65     H = d2(H)
  66     HP  = d2(H @ P)            # precompute a part used twice
  67     K   = d2(np.linalg.solve( d2(HP @ H.T) + R, HP)).T  # Kalman gain
  68     # print('H',H)
  69     # print('K',K)
  70     res = d1(H @ d1(uf) - d)          # res = H*uf - d
  71     ua = uf - K @ res # analysis mean uf - K*res
  72     Pa = Pf - K @ d2(H @ P)        # analysis covariance
  73     return ua, d2(Pa)
  74
  75 ### Define model function with drying, wetting, and rain equilibria
  76
  77 # Parameters
  78 r0 = 0.05                                   # threshold rainfall [mm/h]
  79 rs = 8.0                                    # saturation rain intensity [mm/h]
  80 Tr = 14.0                                   # time constant for rain wetting model [h]
  81 S = 250                                     # saturation intensity [dimensionless]
  82 T = 10.0                                    # time constant for wetting/drying
  83
  84 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  85
  86 def model_moisture(m0,Eqd,Eqw,r,t=None,partials=0,T=10.0,tlen=1.0):
  87     # arguments:
  88     # m0         starting fuel moistureb (%s
  89     # Eqd        drying equilibrium      (%)
  90     # Eqw        wetting equilibrium     (%)
  91     # r          rain intensity          (mm/h)
  92     # t          time
  93     # partials = 0, 1, 2
  94     # returns: same as model_decay
  95     #   if partials==0: m1 = fuel moisture contents after time 1 hour
  96     #              ==1: m1, dm1/dm0
  97     #              ==2: m1, dm1/dm0, dm1/dE
  98
  99     if r > r0:
 100         # print('raining')
 101         E = S
 102         T1 =  (1.0 - np.exp(- (r - r0) / rs)) / Tr
 103     elif m0 <= Eqw:
 104         # print('wetting')
 105         E=Eqw
 106         T1 = 1.0/T
 107     elif m0 >= Eqd:
 108         # print('drying')
 109         E=Eqd
 110         T1 = 1.0/T
 111     else: # no change'
 112         E = m0
 113         T1=0.0
 114     exp_t = np.exp(-tlen*T1)
 115     m1 = E + (m0 - E)*exp_t
 116     dm1_dm0 = exp_t
 117     dm1_dE = 1 - exp_t
 118     #if t>=933 and t < 940:
 119     #  print('t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE',
 120     #        t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE)
 121     if partials==0:
 122         return m1
 123     if partials==1:
 124         return m1, dm1_dm0
 125     if partials==2:
 126         return m1, dm1_dm0, dm1_dE
 127     raise('bad partials')
 128
 129 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 130
 131 ##  NOT TESTED
 132 def model_moisture_run(Eqd,Eqw,r,hours=None,T=10.0,tlen=1.0):
 133 # for arrays of FMC model input run the fuel moisture model
 134     if hours is None:
 135         hours = min(len(Eqd),len(Eqw),len(r))
 136     m = np.zeros(hours)
 137     m[0]=(Eqd[0]+Eqw[0])/2
 138     for k in range(hours-1):
 139         m[k+1]=model_moisture(m[k],Eqd[k],Eqw[k],r[k],T=T,tlen=tlen)
 140     return m
 141
 142 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 143
 144 def model_augmented(u0,Ed,Ew,r,t):
 145     # state u is the vector [m,dE] with dE correction to equilibria Ed and Ew at t
 146     #
 147     m0, Ec = u0  # decompose state u0
 148     # reuse model_moisture(m0,Eqd,Eqw,r,partials=0):
 149     # arguments:
 150     # m0         starting fuel moistureb (1)
 151     # Ed         drying equilibrium      (1)
 152     # Ew         wetting equilibrium     (1)
 153     # r          rain intensity          (mm/h)
 154     # partials = 0, 1, 2
 155     # returns: same as model_decay
 156     #   if partials==0: m1 = fuel moisture contents after time 1 hour
 157     #              ==1: m1, dm0/dm0
 158     #              ==2: m1, dm1/dm0, dm1/dE
 159     m1, dm1_dm0, dm1_dE  = model_moisture(m0,Ed + Ec, Ew + Ec, r, t, partials=2)
 160     u1 = np.array([m1,Ec])   # dE is just copied
 161     J =  np.array([[dm1_dm0, dm1_dE],
 162                    [0.     ,     1.]])
 163     return u1, J
 164
 165 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 166 ### Default Uncertainty Matrices
 167 Q = np.array([[1e-3, 0.],
 168             [0,  1e-3]]) # process noise covariance
 169 H = np.array([[1., 0.]])  # first component observed
 170 R = np.array([1e-3]) # data variance
 171
 172 def run_augmented_kf(dat,h2=None,hours=None, H=H, Q=Q, R=R):
 173     if h2 is None:
 174         h2 = int(dat['h2'])
 175     if hours is None:
 176         hours = int(dat['hours'])
 177
 178     d = dat['fm']
 179     Ed = dat['Ed']
 180     Ew = dat['Ew']
 181     rain = dat['rain']
 182
 183     u = np.zeros((2,hours))
 184     u[:,0]=[0.1,0.0]       # initialize,background state
 185     P = np.zeros((2,2,hours))
 186     P[:,:,0] = np.array([[1e-3, 0.],
 187                       [0.,  1e-3]]) # background state covariance
 188     # Q = np.array([[1e-3, 0.],
 189     #             [0,  1e-3]]) # process noise covariance
 190     # H = np.array([[1., 0.]])  # first component observed
 191     # R = np.array([1e-3]) # data variance
 192
 193     for t in range(1,h2):
 194       # use lambda construction to pass additional arguments to the model
 195         u[:,t],P[:,:,t] = ext_kf(u[:,t-1],P[:,:,t-1],
 196                                   lambda uu: model_augmented(uu,Ed[t],Ew[t],rain[t],t),
 197                                   Q,d[t],H=H,R=R)
 198       # print('time',t,'data',d[t],'filtered',u[0,t],'Ec',u[1,t])
 199     for t in range(h2,hours):
 200         u[:,t],P[:,:,t] = ext_kf(u[:,t-1],P[:,:,t-1],
 201                                   lambda uu: model_augmented(uu,Ed[t],Ew[t],rain[t],t),
 202                                   Q*0.0)
 203       # print('time',t,'data',d[t],'forecast',u[0,t],'Ec',u[1,t])
 204     return u
 205