Fix Pneumatic Valve

Author: teubert

src/prog_models/models/pneumatic_valve.py

@@ -240,16 +240,32 @@ def next_state(self, x, u, dt):
         vdot = pistonForces/params['m']
 
         new_x = x['x']+x['v']*dt
-        if (x['x']==0 and pistonForces<0) or (new_x<0):
-            vel = 0
-            pos = 0
-        elif (x['x']==params['Ls'] and pistonForces>0) or (new_x>params['Ls']):
-            vel = 0
-            pos = params['Ls']
+
+        pos = minimum(maximum(new_x, 0.0), params['Ls'])
+
+        def calc_x(x, forces, Ls, new_x):
+            lower_wall = (x==0 and forces<0) or (new_x<0)
+            upper_wall = (x==Ls and forces>0) or (new_x>Ls)
+            if lower_wall:
+                return 0
+            if upper_wall:
+                return Ls
+            return new_x
+
+        def calc_v(x, v, dv, forces, Ls, new_x):
+            lower_wall = (x==0 and forces<0) or (new_x<0)
+            upper_wall = (x==Ls and forces>0) or (new_x>Ls)
+            if lower_wall or upper_wall:
+                return 0
+            return v + dv
+
+        if isscalar(pistonForces):
+            vel = calc_v(x['x'], x['v'], vdot*dt, pistonForces, params['Ls'], new_x)
+            pos = calc_x(x['x'], pistonForces, params['Ls'], new_x)
         else:
-            # moving
-            vel = x['v'] + vdot*dt
-            pos = new_x
+            # If array- run for each element
+            vel = [calc_v(xi, vi, vdot_i*dt, force, params['Ls'], new_x_i) for xi, vi, vdot_i, force, new_x_i in zip(x['x'], x['v'], vdot, pistonForces, new_x)]
+            vel = [calc_x(xi, force, params['Ls'], new_x_i) for xi, force, new_x_i in zip(x['x'], pistonForces, new_x)]
 
         return {
             'x': pos,

tests/test_pneumatic_valve.py

@@ -2,8 +2,52 @@
 
 import unittest
 from prog_models.models.pneumatic_valve import PneumaticValve, PneumaticValveWithWear, PneumaticValveBase
+from numpy import array
 
 class TestPneumaticValve(unittest.TestCase):
+    def test_pneumatic_valve_vectorized(self):
+        m = PneumaticValveWithWear(process_noise= 0)    
+
+        cycle_time = 20
+        def future_loading(t, x=None):
+            t = t % cycle_time
+            if t < cycle_time/2:
+                return {
+                    'pL': 3.5e5,
+                    'pR': 2.0e5,
+                    # Open Valve
+                    'uTop': False,
+                    'uBot': True
+                }
+            return {
+                'pL': 3.5e5,
+                'pR': 2.0e5,
+                # Close Valve
+                'uTop': True,
+                'uBot': False
+            }
+
+        x0 = {
+            'x': array([0]*4),
+            'v': array([0]*4),
+            'Ai': array([0]*4),
+            'r': array([6000]*4),
+            'k': array([48000]*4),
+            'Aeb': array([1e-5]*4),
+            'Aet': array([1e-5]*4),
+            'condition': array([1]*4),
+            'mTop': array([0.069, 0.048, 0.027, 0.06]),
+            'mBot': array([9.45e-4, 0.008, 0.016, 0.024]),
+            'pDiff': array([1.5e5]*4),
+            'wb': array([0]*4),
+            'wi': array([0]*4),
+            'wk': array([0]*4),
+            'wr': array([0]*4),
+            'wt': array([0]*4)
+        }
+
+        x = m.next_state(x0, future_loading(0), 0.1)
+
     def test_pneumatic_valve_with_wear(self):
         # Test using PneumaticValveWithWear
         m = PneumaticValveWithWear(process_noise= 0)