Adding SDE based noise process as default.

Author: matt-graham

pompy/models.py

@@ -290,8 +290,9 @@ class WindModel(object):
     """
 
     def __init__(self, sim_region=None, n_x=21, n_y=21, u_av=1., v_av=0.,
-                 k_x=2., k_y=2., noise_gain=20., noise_damp=0.1,
-                 noise_bandwidth=0.2, rng=np.random):
+                 k_x=20., k_y=20., noise_gain=2., noise_damp=0.1,
+                 noise_bandwidth=0.2, use_original_noise_updates=False,
+                 rng=np.random):
         """
         Parameters
         ----------
@@ -323,6 +324,10 @@ def __init__(self, sim_region=None, n_x=21, n_y=21, u_av=1., v_av=0.,
         noise_bandwidth : float
             Bandwidth for boundary condition noise generation (dimension:
             angle / time).
+        use_original_noise_updates : boolean
+            Whether to use the original non-SDE based updates for the noise
+            process as defined in Farrell et al. (2002), see notes in
+            `ColouredNoiseGenerator` documentation.
         rng : RandomState
             Random number generator to use in generating input noise. Defaults
             to `numpy.random` global generator however a seeded `RandomState`
@@ -342,7 +347,8 @@ def __init__(self, sim_region=None, n_x=21, n_y=21, u_av=1., v_av=0.,
         # for both components of the wind velocity field so (2,8) state
         # vector (2 as state includes first derivative)
         self.noise_gen = ColouredNoiseGenerator(
-            np.zeros((2, 8)), noise_damp, noise_bandwidth, noise_gain, rng)
+            np.zeros((2, 8)), noise_damp, noise_bandwidth, noise_gain,
+            use_original_noise_updates, rng)
         # compute grid node spacing
         self.dx = sim_region.w / (n_x - 1)  # x grid point spacing
         self.dy = sim_region.h / (n_y - 1)  # y grid point spacing
@@ -486,12 +492,44 @@ def _centred_second_diffs(self, f):
 class ColouredNoiseGenerator(object):
     """Generator of coloured (correlated) Gaussian noise process.
 
-    Generates a coloured noise output via integration of a state space
-    system formulation.
+    Generates a coloured noise output via numerical integration of a stochastic
+    differential equation formulation. The system is assumed to be defined by
+    the system of SDEs::
+
+        dx_0  = x_1 * dt
+        dx_1  = -(a * x_0 + b * x_1) * dt + c * dn
+
+    where `a = bandwidth**2` and `b = 2 * damping * bandwidth`,
+    `c = gain * bandwidth**2` and `dn` is a standard Gaussian white noise
+    process. This is numerically integrated using an Euler-Maruyama scheme::
+
+        for t in range(n_timestep):
+            x[t+1,0] = x[t,0] + dt * x[t,1]
+            x[t+1,1] = x[t,1] - dt * (a*x[t,0] + b*x[t,1]) + dt**0.5 * c * n[t]
+
+    where `x` is an array of shape `(n_timestep, 2)` and `n` is an array of
+    shape `(n_timestep)` filled with random standard normal draws.
+
+    This differs from the code accompanying Farrell et al. (2002) which applies
+    an Euler integration scheme to a state space formulation of a second-order
+    linear system with Gaussian noise input at each time step, resulting in
+    updates of the form::
+
+        for t in range(n_timestep):
+            x[t+1,0] = x[t,0] + dt * x[t,1]
+            x[t+1,1] = x[t,1] + dt * (-a * x[t,0] - b * x[t,1] + c * n[t])
+
+    This differs from the scheme implemented here by scaling the noise input
+    by the timestep `dt` rather than its square root `dt**0.5`. This introduces
+    an implicit dependence of the amplitude of the process on `dt`, in
+    particular that the amplitude scales roughly as `dt**0.5`.
+
+    Updates consistent with the Farrell et al. (2002) implementation can be
+    achieved by setting the `use_original_updates` flat to `True`.
     """
 
-    def __init__(self, init_state, damping=0.1, bandwidth=0.2, gain=20.,
-                 rng=np.random):
+    def __init__(self, init_state, damping=0.1, bandwidth=0.2, gain=1.,
+                 use_original_updates=False,  rng=np.random):
         """
         Parameters
         ----------
@@ -507,14 +545,16 @@ def __init__(self, init_state, damping=0.1, bandwidth=0.2, gain=20.,
         bandwidth : float
             Bandwidth or equivalently undamped natural frequency of system,
             affects system reponsiveness to variations in (noise) input
-            (dimension = angular measure / time).
+            (dimension = 1 / time).
         gain : float
-            Input gain of system, affects scaling of (noise) input
-            (dimensionless).
+            Input gain of system, affects scaling of (noise) input.
         rng : RandomState
             Random number generator to use in generating input noise. Defaults
             to `numpy.random` global generator however a seeded `RandomState`
             object can be passed if it is desired to have reproducible output.
+        use_original_updates : boolean
+            Whether to use the original non-SDE based updates for the noise
+            process as defined in Farrell et al. (2002), see above notes.
         """
         # set up state space matrices
         self.a_mtx = np.array([
@@ -523,6 +563,7 @@ def __init__(self, init_state, damping=0.1, bandwidth=0.2, gain=20.,
         # initialise state
         self.state = init_state
         self.rng = rng
+        self.use_original_updates = use_original_updates
 
     @property
     def output(self):
@@ -539,6 +580,10 @@ def update(self, dt):
         """
         # get normal random input
         n = self.rng.normal(size=(1, self.state.shape[1]))
-        # apply update with Euler-Maruyama integration
-        self.state += (
-            self.a_mtx.dot(self.state) * dt + self.b_mtx * n * dt**0.5)
+        if self.use_original_updates:
+            # apply Farrell et al. (2002) update
+            self.state += dt * (self.a_mtx.dot(self.state) + self.b_mtx.dot(n))
+        else:
+            # apply update with Euler-Maruyama integration
+            self.state += (
+                dt * self.a_mtx.dot(self.state) + self.b_mtx.dot(n) * dt**0.5)