Added Dynamic step size feature

examples/__init__.py

@@ -1,4 +1,4 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the
 # National Aeronautics and Space Administration.  All Rights Reserved.
 
-__all__ = ['sim', 'sim_valve', 'sim_pump', 'sim_battery_eol', 'model_gen', 'benchmarking', 'new_model', 'sensitivity', 'noise', 'future_loading', 'param_est', 'derived_params', 'state_limits']
+__all__ = ['sim', 'sim_valve', 'sim_pump', 'sim_battery_eol', 'model_gen', 'benchmarking', 'new_model', 'sensitivity', 'noise', 'future_loading', 'param_est', 'derived_params', 'state_limits', 'dynamic_step_size']

examples/dynamic_step_size.py

@@ -0,0 +1,56 @@
+# Copyright © 2021 United States Government as represented by the Administrator of the
+# National Aeronautics and Space Administration.  All Rights Reserved.
+
+"""
+Example demonstrating ways to use the dynamic step size feature. This feature allows users to define a time-step that changes with time or state. Run using the command `python -m examples.dynamic_step_size`
+"""
+
+import prog_models
+from .new_model import ThrownObject
+
+def run_example():
+    print("EXAMPLE 1: dt of 1 until 8 sec, then 0.5\n\nSetting up...\n")
+    # Step 1: Create instance of model
+    m = ThrownObject()
+
+    # Step 2: Setup for simulation 
+    def future_load(t, x=None):
+        return {}
+
+    # Step 3: Define dynamic step size function
+    # This `next_time` function will specify what the next step of the simulation should be at any state and time. 
+    # f(x, t) -> (t, dt)
+    def next_time(t, x):
+        # In this example dt is a function of time. We will use a dt of 1 for the first 8 seconds, then 0.5 
+        if t < 8:
+            return 1
+        return 0.5
+
+    # Step 4: Simulate to impact
+    # Here we're printing every time step so we can see the step size change
+    print('\n\n------------------------------------------------')
+    print('Simulating to threshold\n\n')
+    (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(future_load, {'x':m.parameters['thrower_height']}, save_freq=1e-99, print=True, dt=next_time, threshold_keys=['impact'])
+
+    # Example 2
+    print("EXAMPLE 2: dt of 1 until impact event state 0.5, then 0.25 \n\nSetting up...\n")
+
+    # Step 3: Define dynamic step size function
+    # This `next_time` function will specify what the next step of the simulation should be at any state and time. 
+    # f(x, t) -> (t, dt)
+    def next_time(t, x):
+        # In this example dt is a function of state. Uses a dt of 1 until impact event state 0.5, then 0.25
+        event_state = m.event_state(x)
+        if event_state['impact'] < 0.5:
+            return 0.25
+        return 1
+
+    # Step 4: Simulate to impact
+    # Here we're printing every time step so we can see the step size change
+    print('\n\n------------------------------------------------')
+    print('Simulating to threshold\n\n')
+    (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(future_load, {'x':m.parameters['thrower_height']}, save_freq=1e-99, print=True, dt=next_time, threshold_keys=['impact'])
+
+# This allows the module to be executed directly 
+if __name__ == '__main__':
+    run_example()

sphinx_config/getting_started.rst

@@ -66,6 +66,9 @@ See the below examples for examples of use. Run these examples using the command
 * :download:`examples.state_limits <../examples/state_limits.py>`
     .. automodule:: examples.state_limits
     |
+* :download:`examples.dynamic_step_size <../examples/dynamic_step_size.py>`
+    .. automodule:: examples.dynamic_step_size
+    |
 
 There is also an included tutorial (:download:`tutorial <../tutorial.ipynb>`).
 

src/prog_models/prognostics_model.py

@@ -706,13 +706,13 @@ def simulate_to_threshold(self, future_loading_eqn, first_output, threshold_keys
             Configuration options for the simulation \n
             Note: configuration of the model is set through model.parameters.\n
             Supported parameters:\n
-             * dt (Number0: time step (s), e.g. {'dt': 0.1} \n
+             * dt (Number or function): time step (s), e.g. {'dt': 0.1} or function (t, x) -> dt\n
              * save_freq (Number): Frequency at which output is saved (s), e.g., save_freq = 10 \n
              * save_pts ([Number]): Additional ordered list of custom times where output is saved (s), e.g., save_pts= [50, 75] \n
              * horizon (Number): maximum time that the model will be simulated forward (s), e.g., horizon = 1000 \n
              * x (dict): optional, initial state dict, e.g., x= {'x1': 10, 'x2': -5.3}\n
              * thresholds_met_eqn (function/lambda): optional, custom equation to indicate logic for when to stop sim f(thresholds_met) -> bool\n
-             * print_inter (bool): optional, toggle intermediate printing, e.g., print_inter = True\n
+             * print (bool): optional, toggle intermediate printing, e.g., print_inter = True\n
             e.g., m.simulate_to_threshold(eqn, z, dt=0.1, save_pts=[1, 2])
 
         Returns
@@ -768,9 +768,9 @@ def simulate_to_threshold(self, future_loading_eqn, first_output, threshold_keys
         config.update(kwargs)
 
         # Configuration validation
-        if not isinstance(config['dt'], Number):
-            raise ProgModelInputException("'dt' must be a number, was a {}".format(type(config['dt'])))
-        if config['dt'] <= 0:
+        if not isinstance(config['dt'], Number) and not callable(config['dt']):
+            raise ProgModelInputException("'dt' must be a number or function, was a {}".format(type(config['dt'])))
+        if isinstance(config['dt'], Number) and config['dt'] < 0:
             raise ProgModelInputException("'dt' must be positive, was {}".format(config['dt']))
         if not isinstance(config['save_freq'], Number):
             raise ProgModelInputException("'save_freq' must be a number, was a {}".format(type(config['save_freq'])))
@@ -818,7 +818,6 @@ def check_thresholds(thresholds_met):
         states = []  
         saved_outputs = []
         saved_event_states = []
-        dt = config['dt']  # saving to optimize access in while loop
         save_freq = config['save_freq']
         horizon = config['horizon']
         next_save = save_freq
@@ -846,11 +845,20 @@ def update_all():
                 times.append(t)
                 inputs.append(u)
                 states.append(deepcopy(x))  # Avoid optimization where x is not copied
+
+        # configuring next_time function to define prediction time step, default is constant dt
+        if callable(config['dt']):
+            next_time = config['dt']
+        else:
+            dt = config['dt']  # saving to optimize access in while loop
+            def next_time(t, x):
+                return dt
 
         # Simulate
         update_all()
         while t < horizon:
-            t += dt
+            dt = next_time(t, x)
+            t = t + dt
             u = future_loading_eqn(t, x)
             x = next_state(x, u, dt)
             if (t >= next_save):

tests/test_base_models.py

@@ -570,6 +570,23 @@ def load(t, x=None):
         (times, inputs, states, outputs, event_states) = m.simulate_to(6, load, {'o1': 0.8}, **{'dt': 0.5, 'save_freq': 1.0})
         self.assertAlmostEqual(times[-1], 6.0, 5)
 
+    def test_next_time_fcn(self):
+        m = MockProgModel(process_noise = 0.0)
+        def load(t, x=None):
+            return {'i1': 1, 'i2': 2.1}
+
+
+        # Any event, default
+        (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(load, {'o1': 0.8}, **{'dt': 1, 'save_freq': 1e-99})
+        self.assertEqual(len(times), 6)
+
+        def next_time(t, x):
+            return 0.5
+
+        # With next_time
+        (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(load, {'o1': 0.8}, **{'save_freq': 1e-99, 'dt': next_time})
+        self.assertEqual(len(times), 11)
+
     def test_sim_prog(self):
         m = MockProgModel(process_noise = 0.0)
         def load(t, x=None):

tests/test_examples.py

@@ -110,6 +110,17 @@ def test_noise_example(self):
         # Reset stdout 
         sys.stdout = _stdout
 
+    def test_dynamic_step(self):
+        # set stdout (so it wont print)
+        _stdout = sys.stdout
+        sys.stdout = StringIO()
+
+        # Run example
+        dynamic_step_size.run_example()
+
+        # Reset stdout 
+        sys.stdout = _stdout
+
     def test_future_loading(self):
         # set stdout (so it wont print)
         _stdout = sys.stdout