Added current limit estimation algorithm to electrochem model

examples/observables_check.py

@@ -0,0 +1,51 @@
+# Copyright © 2021 United States Government as represented by the Administrator of the
+# National Aeronautics and Space Administration.  All Rights Reserved.
+
+"""
+Example of implementation of CLE algorithm to determine maximum current possible from battery on time interval 
+"""
+
+from timeit import timeit
+from prog_models.models import BatteryElectroChemEOD
+
+def run_example():
+    # Step 1: Create a model object
+    batt = BatteryElectroChemEOD()
+
+    # Step 2: Define future loading function 
+    def future_loading(t, x=None):
+        # Variable (piece-wise) future loading scheme 
+        if (t < 600):
+            i = 2
+        elif (t < 900):
+            i = 1
+        elif (t < 1800):
+            i = 4
+        elif (t < 3000):
+            i = 2
+        else:
+            i = 3
+        return {'i': i}
+
+    # Define parameters for CLE test calculations: 
+    time_to_sim_to = 300 # time cutoff value/time interval of simulation 
+    V_cutoff = 3.2 # voltage cutoff value 
+
+    # Initialize to default parameters
+    initial_state = batt.parameters['x0']
+
+    # Run CLE algorithm from default inital state
+    i_start = future_loading(0)
+    CLE_estimate = batt.current_limit_est(initial_state, i_start, delta_t=time_to_sim_to, VCutoff=V_cutoff)
+
+    # Simulate to 600 sec
+    (times, inputs, states, outputs, event_states) = batt.simulate_to(600, future_loading)
+
+    # Using state at 600 sec, run CLE algorithm 
+    CLE_estimate_2 = batt.current_limit_est(states[-1], inputs[-1], delta_t=time_to_sim_to, VCutoff=V_cutoff)
+
+    debug_check = 1
+
+# This allows the module to be executed directly 
+if __name__=='__main__':
+    run_example()

src/prog_models/models/battery_circuit.py

@@ -78,6 +78,7 @@ class BatteryCircuit(PrognosticsModel):
     states = ['tb', 'qb', 'qcp', 'qcs']
     outputs = ['t',  'v']
     is_vectorized = True
+    observables_keys = ['currentMin', 'currentMax']
 
     default_parameters = {  # Set to defaults
         'V0': 4.183,
@@ -187,3 +188,13 @@ def threshold_met(self, x):
         return {
              'EOD': V < parameters['VEOD']
         }
+
+    def observables(self, x) -> dict:
+        params = self.parameters
+        nomCapacity = params['nomCapacity']
+        CRateMin = params['CRateMin']
+        CRateMax = params['CRateMax']
+        return {
+            'currentMin': nomCapacity * CRateMin,
+            'currentMax': nomCapacity * CRateMax,
+        }

src/prog_models/models/battery_electrochem.py

@@ -190,6 +190,7 @@ class BatteryElectroChemEOD(PrognosticsModel):
     events = ['EOD']
     inputs = ['i']
     states = ['tb', 'Vo', 'Vsn', 'Vsp', 'qnB', 'qnS', 'qpB', 'qpS']
+    observables_keys = ['currentMin', 'currentMax']
     outputs = ['t', 'v']
     is_vectorized = True
 
@@ -235,7 +236,12 @@ class BatteryElectroChemEOD(PrognosticsModel):
 
         # End of discharge voltage threshold
         'VEOD': 3.0, 
-        'VDropoff': 0.1 # Voltage above EOD after which voltage will be considered in SOC calculation
+        'VDropoff': 0.1, # Voltage above EOD after which voltage will be considered in SOC calculation
+        # current ratings
+        'nomCapacity': 2.2,  # nominal capacity, Ah
+        'CRateMin': 0.7,  # current necessary for cruise,
+        'CRateMax': 2.5   # current necessary for hover
+        # CRateMin, CRateMax based on values determined in `C. Silva and W. Johnson, "VTOL Urban Air Mobility Concept Vehicles for Technology Development" Aviation and Aeronautics Forum (Aviation 2018),June 2018. https://arc.aiaa.org/doi/abs/10.2514/6.2018-3847`
     }
 
     state_limits = {
@@ -327,6 +333,118 @@ def event_state(self, x):
             'EOD': min(charge_EOD, voltage_EOD)
         }
 
+    # for use with current_limit_est below
+    def get_OCV(self, x):
+        params = self.parameters
+        An = params['An']
+        # Negative Surface
+        xnS = x['qnS']/params['qSMax']
+        xnS2 = xnS+xnS  # Note: in python x+x is more efficient than 2*x
+        VenParts = [
+            An[0] * (xnS2-1)/F,  # Ven0
+            An[1] * ((xnS2-1)**2 - ((xnS + xnS)*(1-xnS)))/F,  # Ven1
+            An[2] * ((xnS2-1)**3 - (4 * xnS*(1-xnS))*(xnS2-1))/F,  # Ven2
+            An[3] * ((xnS2-1)**4 - (6 * xnS*(1-xnS))*(xnS2-1)**2) / F,  # Ven3
+            An[4] * ((xnS2-1)**5 - (8 * xnS*(1-xnS))*(xnS2-1)**3) / F,  # Ven4
+            An[5] * ((xnS2-1)**6 - (10*xnS*(1-xnS))*(xnS2-1)**4) / F,  # Ven5
+            An[6] * ((xnS2-1)**7 - (12*xnS*(1-xnS))*(xnS2-1)**5) / F,  # Ven6
+            An[7] * ((xnS2-1)**8 - (14*xnS*(1-xnS))*(xnS2-1)**6) / F,  # Ven7
+            An[8] * ((xnS2-1)**9 - (16*xnS*(1-xnS))*(xnS2-1)**7) / F,  # Ven8
+            An[9] * ((xnS2-1)**10 - (18*xnS*(1-xnS))*(xnS2-1)**8) / F,  # Ven9
+            An[10]*((xnS2-1)**11 - (20*xnS*(1-xnS))*(xnS2-1)**9) / F,  # Ven10
+            An[11]*((xnS2-1)**12 - (22*xnS*(1-xnS))*(xnS2-1)**10)/F,  # Ven11
+            An[12]*((xnS2-1)**13 - (24*xnS*(1-xnS))*(xnS2-1)**11)/F  # Ven12
+        ]
+
+        Ven = params['U0n'] + R*x['tb']/F*log((1-xnS)/xnS) + sum(VenParts)
+
+        # Positive Surface
+        Ap = params['Ap']
+        xpS = x['qpS']/params['qSMax']
+        xpS2 = xpS + xpS
+        VepParts = [
+            Ap[0] * (xpS2-1)/F,  # Vep0
+            Ap[1] * ((xpS2-1)**2 - (xpS2*(1-xpS)))/F,  # Vep1
+            Ap[2] * ((xpS2-1)**3 - (4 * xpS*(1-xpS)) / \
+                     (xpS2-1)**(-1)) / F,  # Vep2
+            Ap[3] * ((xpS2-1)**4 - (6 * xpS*(1-xpS)) / \
+                     (xpS2-1)**(-2)) / F,  # Vep3
+            Ap[4] * ((xpS2-1)**5 - (8 * xpS*(1-xpS)) / \
+                     (xpS2-1)**(-3)) / F,  # Vep4
+            Ap[5] * ((xpS2-1)**6 - (10*xpS*(1-xpS)) / \
+                     (xpS2-1)**(-4)) / F,  # Vep5
+            Ap[6] * ((xpS2-1)**7 - (12*xpS*(1-xpS)) / \
+                     (xpS2-1)**(-5)) / F,  # Vep6
+            Ap[7] * ((xpS2-1)**8 - (14*xpS*(1-xpS)) / \
+                     (xpS2-1)**(-6)) / F,  # Vep7
+            Ap[8] * ((xpS2-1)**9 - (16*xpS*(1-xpS)) / \
+                     (xpS2-1)**(-7)) / F,  # Vep8
+            Ap[9] * ((xpS2-1)**10 - (18*xpS*(1-xpS)) / \
+                     (xpS2-1)**(-8)) / F,  # Vep9
+            Ap[10]*((xpS2-1)**11 - (20*xpS*(1-xpS)) / \
+                    (xpS2-1)**(-9)) / F,  # Vep10
+            Ap[11]*((xpS2-1)**12 - (22*xpS*(1-xpS)) / \
+                    (xpS2-1)**(-10))/F,  # Vep11
+            Ap[12]*((xpS2-1)**13 - (24*xpS*(1-xpS))/(xpS2-1)**(-11))/F  # Vep12
+        ]
+        Vep = params['U0p'] + R*x['tb']/F*log((1-xpS)/xpS) + sum(VepParts)
+
+        return {
+            'OCV': Vep - Ven
+        }
+
+    # This method implements a current limit estimate algorithm as described in the following paper:
+    # S. Bharathraj et al., "Accessing Current Limits in Lithium Ion Batteries: Analysis of Propensity for Unexpected Power Loss
+    # as a Function of Depth of Discharge, Temperature, and Pulse Duration," Journal of Power Sources, 494, 2021
+    def current_limit_est(self, x, u, delta_t, VCutoff):
+        params = self.parameters
+        z = self.output(x)  # provides present voltage and temp
+        w = self.get_OCV(x)  # provides present open current voltage
+
+        # Define approximate for R_lump
+            # Currently estimated with R_ohm, defined as R0
+            # Potential future work (from Elizabeth Hale):
+                # VCell = z['v']
+                # VOpen = w['OCV']
+                # RLump = (VOpen - VCell)/u['i']
+        R_lump = params['Ro']
+
+        # Initialize values 
+        CLE_val = u['i']  # initialize CLE_val
+        relative_change = 1   # initialize relative change
+        last_state = x  # initialize last state
+        original_state = params['x0'] # save original initialization parameters to restore later
+        params['x0'] = last_state  # start model simulation at present state
+        first_output = z  # first output based on present state
+        time_to_simulate_to = delta_t  # time cutoff
+        sim_config = {'save_freq': delta_t, 'dt': 0.1}
+        while relative_change > 0.001:  # convergence criteria in Bharathraj et al.
+            def cle_loading(t, x=None):  # load for simulation is CLE_val
+                return {'i': CLE_val}
+            (times, inputs, states, outputs, event_states) = self.simulate_to(time_to_simulate_to,
+                                                                              cle_loading, first_output, **sim_config)
+            last_state = states[1]
+            V_OCV_temp = self.get_OCV(last_state)
+            V_OCV = V_OCV_temp['OCV']
+            CLE_old = CLE_val
+            CLE_val = (V_OCV - VCutoff)/R_lump
+            relative_change = abs((CLE_val - CLE_old)/CLE_old)
+        params['x0'] = original_state
+        a=1
+        return {
+            'CLE': CLE_val,
+            'final_voltage': outputs[1]['v']
+        }
+
+    def observables(self, x) -> dict:
+        params = self.parameters
+        nomCapacity = params['nomCapacity']
+        CRateMin = params['CRateMin']
+        CRateMax = params['CRateMax']
+        return {
+            'currentMin': nomCapacity * CRateMin,
+            'currentMax': nomCapacity * CRateMax,
+        }
     def output(self, x):
         params = self.parameters
         An = params['An']