diff --git a/src/progpy/models/battery_electrochem.py b/src/progpy/models/battery_electrochem.py index 497eb50f..b2fe2363 100644 --- a/src/progpy/models/battery_electrochem.py +++ b/src/progpy/models/battery_electrochem.py @@ -1,10 +1,9 @@ # Copyright © 2021 United States Government as represented by the Administrator of the # National Aeronautics and Space Administration. All Rights Reserved. -from copy import deepcopy import numpy as np from scipy.optimize import fsolve -import warnings +from copy import deepcopy from progpy import PrognosticsModel @@ -137,6 +136,132 @@ def update_qSBmax(params: dict) -> dict: 'qBMax': params['qMax']*(1.0-params['VolSFraction']), } +def calculate_temp_voltage(x, params, qSMax): + ''' + Calculate and return the temperature and voltage of the battery. + + :param x: battery state + :param params: battery parameters + :param qSMax: maximum charge at the surface + :return t: battery temperature + :return v: battery voltage + ''' + An = params['An'] + # Negative Surface + xnS = x['qnS']/qSMax + xnS2 = xnS+xnS # Note: in python x+x is more efficient than 2*x + + one_minus_xnS = 1 - xnS + xnS2_minus_1 = xnS2 - 1 + VenParts = [ + An[0] *xnS2_minus_1/F, # Ven0 + An[1] *(xnS2_minus_1**2 - (xnS2*one_minus_xnS))/F, # Ven1 + An[2] *(xnS2_minus_1**3 - (4 *xnS*one_minus_xnS)*xnS2_minus_1)/F, #Ven2 + An[3] *(xnS2_minus_1**4 - (6 *xnS*one_minus_xnS)*xnS2_minus_1**2) /F, #Ven3 + An[4] *(xnS2_minus_1**5 - (8 *xnS*one_minus_xnS)*xnS2_minus_1**3) /F, #Ven4 + An[5] *(xnS2_minus_1**6 - (10*xnS*one_minus_xnS)*xnS2_minus_1**4) /F, #Ven5 + An[6] *(xnS2_minus_1**7 - (12*xnS*one_minus_xnS)*xnS2_minus_1**5) /F, #Ven6 + An[7] *(xnS2_minus_1**8 - (14*xnS*one_minus_xnS)*xnS2_minus_1**6) /F, #Ven7 + An[8] *(xnS2_minus_1**9 - (16*xnS*one_minus_xnS)*xnS2_minus_1**7) /F, #Ven8 + An[9] *(xnS2_minus_1**10 - (18*xnS*one_minus_xnS)*xnS2_minus_1**8) /F, #Ven9 + An[10]*(xnS2_minus_1**11 - (20*xnS*one_minus_xnS)*xnS2_minus_1**9) /F, #Ven10 + An[11]*(xnS2_minus_1**12 - (22*xnS*one_minus_xnS)*xnS2_minus_1**10)/F, #Ven11 + An[12]*(xnS2_minus_1**13 - (24*xnS*one_minus_xnS)*xnS2_minus_1**11)/F #Ven12 + ] + Ven = params['U0n'] + R*x['tb']/F*np.log(one_minus_xnS/xnS) + sum(VenParts) + + # Positive Surface + Ap = params['Ap'] + xpS = x['qpS']/qSMax + one_minus_xpS = 1 - xpS + xpS2 = xpS + xpS + xpS2_minus_1 = xpS2 - 1 + VepParts = [ + Ap[0] *(xpS2_minus_1)/F, #Vep0 + Ap[1] *((xpS2_minus_1)**2 - (xpS2*one_minus_xpS))/F, #Vep1 + Ap[2] *((xpS2_minus_1)**3 - (4 *xpS*one_minus_xpS)*(xpS2_minus_1)) /F, #Vep2 + Ap[3] *((xpS2_minus_1)**4 - (6 *xpS*one_minus_xpS)*(xpS2_minus_1)**(2)) /F, #Vep3 + Ap[4] *((xpS2_minus_1)**5 - (8 *xpS*one_minus_xpS)*(xpS2_minus_1)**(3)) /F, #Vep4 + Ap[5] *((xpS2_minus_1)**6 - (10*xpS*one_minus_xpS)*(xpS2_minus_1)**(4)) /F, #Vep5 + Ap[6] *((xpS2_minus_1)**7 - (12*xpS*one_minus_xpS)*(xpS2_minus_1)**(5)) /F, #Vep6 + Ap[7] *((xpS2_minus_1)**8 - (14*xpS*one_minus_xpS)*(xpS2_minus_1)**(6)) /F, #Vep7 + Ap[8] *((xpS2_minus_1)**9 - (16*xpS*one_minus_xpS)*(xpS2_minus_1)**(7)) /F, #Vep8 + Ap[9] *((xpS2_minus_1)**10 - (18*xpS*one_minus_xpS)*(xpS2_minus_1)**(8)) /F, #Vep9 + Ap[10]*((xpS2_minus_1)**11 - (20*xpS*one_minus_xpS)*(xpS2_minus_1)**(9)) /F, #Vep10 + Ap[11]*((xpS2_minus_1)**12 - (22*xpS*one_minus_xpS)*(xpS2_minus_1)**(10))/F, #Vep11 + Ap[12]*((xpS2_minus_1)**13 - (24*xpS*one_minus_xpS)*(xpS2_minus_1)**(11))/F #Vep12 + ] + Vep = params['U0p'] + R*x['tb']/F*np.log(one_minus_xpS/xpS) + sum(VepParts) + + t = np.atleast_1d(x['tb'] - 273.15) + v = np.atleast_1d(Vep - Ven - x['Vo'] - x['Vsn'] - x['Vsp']) + + return t, v + +def calculate_EOD(x, params, qnMax, qSMax): + ''' + Calculate and return the battery EOD. + + :param x: battery state + :param params: battery parameters + :param qnMax: maximum charge at the negative electrode + :param qSMax: maximum charge at the surface + :return: dictionary with battery EOD + ''' + An = params['An'] + # Negative Surface + xnS = x['qnS']/qSMax + xnS2 = xnS+xnS # Note: in python x+x is more efficient than 2*x + one_minus_xnS = 1 - xnS + xnS2_minus_1 = xnS2 - 1 + VenParts = [ + An[0] *xnS2_minus_1/F, # Ven0 + An[1] *(xnS2_minus_1**2 - (xnS2*one_minus_xnS))/F, # Ven1 + An[2] *(xnS2_minus_1**3 - (4 *xnS*one_minus_xnS)*xnS2_minus_1)/F, #Ven2 + An[3] *(xnS2_minus_1**4 - (6 *xnS*one_minus_xnS)*xnS2_minus_1**2) /F, #Ven3 + An[4] *(xnS2_minus_1**5 - (8 *xnS*one_minus_xnS)*xnS2_minus_1**3) /F, #Ven4 + An[5] *(xnS2_minus_1**6 - (10*xnS*one_minus_xnS)*xnS2_minus_1**4) /F, #Ven5 + An[6] *(xnS2_minus_1**7 - (12*xnS*one_minus_xnS)*xnS2_minus_1**5) /F, #Ven6 + An[7] *(xnS2_minus_1**8 - (14*xnS*one_minus_xnS)*xnS2_minus_1**6) /F, #Ven7 + An[8] *(xnS2_minus_1**9 - (16*xnS*one_minus_xnS)*xnS2_minus_1**7) /F, #Ven8 + An[9] *(xnS2_minus_1**10 - (18*xnS*one_minus_xnS)*xnS2_minus_1**8) /F, #Ven9 + An[10]*(xnS2_minus_1**11 - (20*xnS*one_minus_xnS)*xnS2_minus_1**9) /F, #Ven10 + An[11]*(xnS2_minus_1**12 - (22*xnS*one_minus_xnS)*xnS2_minus_1**10)/F, #Ven11 + An[12]*(xnS2_minus_1**13 - (24*xnS*one_minus_xnS)*xnS2_minus_1**11)/F #Ven12 + ] + Ven = params['U0n'] + R*x['tb']/F*np.log(one_minus_xnS/xnS) + sum(VenParts) + + # Positive Surface + Ap = params['Ap'] + xpS = x['qpS']/qSMax + one_minus_xpS = 1 - xpS + xpS2 = xpS + xpS + xpS2_minus_1 = xpS2 - 1 + VepParts = [ + Ap[0] *(xpS2_minus_1)/F, #Vep0 + Ap[1] *((xpS2_minus_1)**2 - (xpS2*one_minus_xpS))/F, #Vep1 + Ap[2] *((xpS2_minus_1)**3 - (4 *xpS*one_minus_xpS)*(xpS2_minus_1)) /F, #Vep2 + Ap[3] *((xpS2_minus_1)**4 - (6 *xpS*one_minus_xpS)*(xpS2_minus_1)**(2)) /F, #Vep3 + Ap[4] *((xpS2_minus_1)**5 - (8 *xpS*one_minus_xpS)*(xpS2_minus_1)**(3)) /F, #Vep4 + Ap[5] *((xpS2_minus_1)**6 - (10*xpS*one_minus_xpS)*(xpS2_minus_1)**(4)) /F, #Vep5 + Ap[6] *((xpS2_minus_1)**7 - (12*xpS*one_minus_xpS)*(xpS2_minus_1)**(5)) /F, #Vep6 + Ap[7] *((xpS2_minus_1)**8 - (14*xpS*one_minus_xpS)*(xpS2_minus_1)**(6)) /F, #Vep7 + Ap[8] *((xpS2_minus_1)**9 - (16*xpS*one_minus_xpS)*(xpS2_minus_1)**(7)) /F, #Vep8 + Ap[9] *((xpS2_minus_1)**10 - (18*xpS*one_minus_xpS)*(xpS2_minus_1)**(8)) /F, #Vep9 + Ap[10]*((xpS2_minus_1)**11 - (20*xpS*one_minus_xpS)*(xpS2_minus_1)**(9)) /F, #Vep10 + Ap[11]*((xpS2_minus_1)**12 - (22*xpS*one_minus_xpS)*(xpS2_minus_1)**(10))/F, #Vep11 + Ap[12]*((xpS2_minus_1)**13 - (24*xpS*one_minus_xpS)*(xpS2_minus_1)**(11))/F #Vep12 + ] + Vep = params['U0p'] + R*x['tb']/F*np.log(one_minus_xpS/xpS) + sum(VepParts) + v = Vep - Ven - x['Vo'] - x['Vsn'] - x['Vsp'] + + charge_EOD = (x['qnS'] + x['qnB'])/qnMax + voltage_EOD = (v - params['VEOD'])/params['VDropoff'] + + return { + 'EOD': np.clip(min(charge_EOD, voltage_EOD), 0, 1) + } + class BatteryElectroChemEOD(PrognosticsModel): """ @@ -220,6 +345,8 @@ class BatteryElectroChemEOD(PrognosticsModel): Redlich-Kister parameter (- electrode) VEOD : float End of Discharge Voltage Threshold + VDropoff : float + Voltage above EOD after which voltage will be considered in SOC calculation x0 : dict[str, float] Initial :term:`state` @@ -423,119 +550,22 @@ def event_state(self, x) -> dict: # longer accurately captures this behavior, so voltage_EOD takes over as # the driving factor. 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 - one_minus_xnS = 1 - xnS - xnS2_minus_1 = xnS2 - 1 - VenParts = [ - An[0] *xnS2_minus_1/F, # Ven0 - An[1] *(xnS2_minus_1**2 - (xnS2*one_minus_xnS))/F, # Ven1 - An[2] *(xnS2_minus_1**3 - (4 *xnS*one_minus_xnS)*xnS2_minus_1)/F, #Ven2 - An[3] *(xnS2_minus_1**4 - (6 *xnS*one_minus_xnS)*xnS2_minus_1**2) /F, #Ven3 - An[4] *(xnS2_minus_1**5 - (8 *xnS*one_minus_xnS)*xnS2_minus_1**3) /F, #Ven4 - An[5] *(xnS2_minus_1**6 - (10*xnS*one_minus_xnS)*xnS2_minus_1**4) /F, #Ven5 - An[6] *(xnS2_minus_1**7 - (12*xnS*one_minus_xnS)*xnS2_minus_1**5) /F, #Ven6 - An[7] *(xnS2_minus_1**8 - (14*xnS*one_minus_xnS)*xnS2_minus_1**6) /F, #Ven7 - An[8] *(xnS2_minus_1**9 - (16*xnS*one_minus_xnS)*xnS2_minus_1**7) /F, #Ven8 - An[9] *(xnS2_minus_1**10 - (18*xnS*one_minus_xnS)*xnS2_minus_1**8) /F, #Ven9 - An[10]*(xnS2_minus_1**11 - (20*xnS*one_minus_xnS)*xnS2_minus_1**9) /F, #Ven10 - An[11]*(xnS2_minus_1**12 - (22*xnS*one_minus_xnS)*xnS2_minus_1**10)/F, #Ven11 - An[12]*(xnS2_minus_1**13 - (24*xnS*one_minus_xnS)*xnS2_minus_1**11)/F #Ven12 - ] - Ven = params['U0n'] + R*x['tb']/F*np.log(one_minus_xnS/xnS) + sum(VenParts) - - # Positive Surface - Ap = params['Ap'] - xpS = x['qpS']/params['qSMax'] - one_minus_xpS = 1 - xpS - xpS2 = xpS + xpS - xpS2_minus_1 = xpS2 - 1 - VepParts = [ - Ap[0] *(xpS2_minus_1)/F, #Vep0 - Ap[1] *((xpS2_minus_1)**2 - (xpS2*one_minus_xpS))/F, #Vep1 - Ap[2] *((xpS2_minus_1)**3 - (4 *xpS*one_minus_xpS)*(xpS2_minus_1)) /F, #Vep2 - Ap[3] *((xpS2_minus_1)**4 - (6 *xpS*one_minus_xpS)*(xpS2_minus_1)**(2)) /F, #Vep3 - Ap[4] *((xpS2_minus_1)**5 - (8 *xpS*one_minus_xpS)*(xpS2_minus_1)**(3)) /F, #Vep4 - Ap[5] *((xpS2_minus_1)**6 - (10*xpS*one_minus_xpS)*(xpS2_minus_1)**(4)) /F, #Vep5 - Ap[6] *((xpS2_minus_1)**7 - (12*xpS*one_minus_xpS)*(xpS2_minus_1)**(5)) /F, #Vep6 - Ap[7] *((xpS2_minus_1)**8 - (14*xpS*one_minus_xpS)*(xpS2_minus_1)**(6)) /F, #Vep7 - Ap[8] *((xpS2_minus_1)**9 - (16*xpS*one_minus_xpS)*(xpS2_minus_1)**(7)) /F, #Vep8 - Ap[9] *((xpS2_minus_1)**10 - (18*xpS*one_minus_xpS)*(xpS2_minus_1)**(8)) /F, #Vep9 - Ap[10]*((xpS2_minus_1)**11 - (20*xpS*one_minus_xpS)*(xpS2_minus_1)**(9)) /F, #Vep10 - Ap[11]*((xpS2_minus_1)**12 - (22*xpS*one_minus_xpS)*(xpS2_minus_1)**(10))/F, #Vep11 - Ap[12]*((xpS2_minus_1)**13 - (24*xpS*one_minus_xpS)*(xpS2_minus_1)**(11))/F #Vep12 - ] - Vep = params['U0p'] + R*x['tb']/F*np.log(one_minus_xpS/xpS) + sum(VepParts) - v = Vep - Ven - x['Vo'] - x['Vsn'] - x['Vsp'] - charge_EOD = (x['qnS'] + x['qnB'])/self.parameters['qnMax'] - voltage_EOD = (v - self.parameters['VEOD'])/self.parameters['VDropoff'] - return { - 'EOD': np.clip(min(charge_EOD, voltage_EOD), 0, 1) - } + return calculate_EOD(x, params, params['qnMax'], params['qSMax']) def output(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 - - one_minus_xnS = 1 - xnS - xnS2_minus_1 = xnS2 - 1 - VenParts = [ - An[0] *xnS2_minus_1/F, # Ven0 - An[1] *(xnS2_minus_1**2 - (xnS2*one_minus_xnS))/F, # Ven1 - An[2] *(xnS2_minus_1**3 - (4 *xnS*one_minus_xnS)*xnS2_minus_1)/F, #Ven2 - An[3] *(xnS2_minus_1**4 - (6 *xnS*one_minus_xnS)*xnS2_minus_1**2) /F, #Ven3 - An[4] *(xnS2_minus_1**5 - (8 *xnS*one_minus_xnS)*xnS2_minus_1**3) /F, #Ven4 - An[5] *(xnS2_minus_1**6 - (10*xnS*one_minus_xnS)*xnS2_minus_1**4) /F, #Ven5 - An[6] *(xnS2_minus_1**7 - (12*xnS*one_minus_xnS)*xnS2_minus_1**5) /F, #Ven6 - An[7] *(xnS2_minus_1**8 - (14*xnS*one_minus_xnS)*xnS2_minus_1**6) /F, #Ven7 - An[8] *(xnS2_minus_1**9 - (16*xnS*one_minus_xnS)*xnS2_minus_1**7) /F, #Ven8 - An[9] *(xnS2_minus_1**10 - (18*xnS*one_minus_xnS)*xnS2_minus_1**8) /F, #Ven9 - An[10]*(xnS2_minus_1**11 - (20*xnS*one_minus_xnS)*xnS2_minus_1**9) /F, #Ven10 - An[11]*(xnS2_minus_1**12 - (22*xnS*one_minus_xnS)*xnS2_minus_1**10)/F, #Ven11 - An[12]*(xnS2_minus_1**13 - (24*xnS*one_minus_xnS)*xnS2_minus_1**11)/F #Ven12 - ] - Ven = params['U0n'] + R*x['tb']/F*np.log(one_minus_xnS/xnS) + sum(VenParts) - - # Positive Surface - Ap = params['Ap'] - xpS = x['qpS']/params['qSMax'] - one_minus_xpS = 1 - xpS - xpS2 = xpS + xpS - xpS2_minus_1 = xpS2 - 1 - VepParts = [ - Ap[0] *(xpS2_minus_1)/F, #Vep0 - Ap[1] *((xpS2_minus_1)**2 - (xpS2*one_minus_xpS))/F, #Vep1 - Ap[2] *((xpS2_minus_1)**3 - (4 *xpS*one_minus_xpS)*(xpS2_minus_1)) /F, #Vep2 - Ap[3] *((xpS2_minus_1)**4 - (6 *xpS*one_minus_xpS)*(xpS2_minus_1)**(2)) /F, #Vep3 - Ap[4] *((xpS2_minus_1)**5 - (8 *xpS*one_minus_xpS)*(xpS2_minus_1)**(3)) /F, #Vep4 - Ap[5] *((xpS2_minus_1)**6 - (10*xpS*one_minus_xpS)*(xpS2_minus_1)**(4)) /F, #Vep5 - Ap[6] *((xpS2_minus_1)**7 - (12*xpS*one_minus_xpS)*(xpS2_minus_1)**(5)) /F, #Vep6 - Ap[7] *((xpS2_minus_1)**8 - (14*xpS*one_minus_xpS)*(xpS2_minus_1)**(6)) /F, #Vep7 - Ap[8] *((xpS2_minus_1)**9 - (16*xpS*one_minus_xpS)*(xpS2_minus_1)**(7)) /F, #Vep8 - Ap[9] *((xpS2_minus_1)**10 - (18*xpS*one_minus_xpS)*(xpS2_minus_1)**(8)) /F, #Vep9 - Ap[10]*((xpS2_minus_1)**11 - (20*xpS*one_minus_xpS)*(xpS2_minus_1)**(9)) /F, #Vep10 - Ap[11]*((xpS2_minus_1)**12 - (22*xpS*one_minus_xpS)*(xpS2_minus_1)**(10))/F, #Vep11 - Ap[12]*((xpS2_minus_1)**13 - (24*xpS*one_minus_xpS)*(xpS2_minus_1)**(11))/F #Vep12 - ] - Vep = params['U0p'] + R*x['tb']/F*np.log(one_minus_xpS/xpS) + sum(VepParts) - - return self.OutputContainer(np.array([ - np.atleast_1d(x['tb'] - 273.15), - np.atleast_1d(Vep - Ven - x['Vo'] - x['Vsn'] - x['Vsp']) - ])) + t, v = calculate_temp_voltage(x, params, params['qSMax']) + + return self.OutputContainer(np.array([t, v])) + def threshold_met(self, x) -> dict: z = self.output(x) # Return true if voltage is less than the voltage threshold return { - 'EOD': z['v'] < self.parameters['VEOD'] + 'EOD': z['v'] < self.parameters['VEOD'] } @@ -639,90 +669,279 @@ def merge_dicts(a: dict, b: dict) -> None: else: a[key] = b[key] -def OverwrittenWarning(params): - """ - Function to warn if overwritten changes - """ - warnings.warn("Ro, qMobile, and tDiffusion will be overwritten within the model as part of battery aging modeling. Use BatteryElectroChemEOD to remove this behavior.") - return {} - -class BatteryElectroChemEODEOL(BatteryElectroChemEOL, BatteryElectroChemEOD): +class BatteryElectroChemEODEOL(PrognosticsModel): """ - Prognostics :term:`model` for a battery degredation and discharge, represented by an electrochemical model as described in [Daigle2013]_ and [Daigle2016]_ + Prognostics :term:`model` for a battery degredation and discharge, represented by an electrochemical model as described in [Daigle2013]_ and [Daigle2016]_. The default model parameters included are for Li-ion batteries, specifically 18650-type cells. Experimental discharge curves for these cells can be downloaded from the Prognostics Center of Excellence Data Repository [DataRepo]_. :term:`Events`: (2) - | EOD: End of Discharge + | EOD: End of discharge | InsufficientCapacity: Insufficient battery capacity :term:`Inputs/Loading`: (1) i: Current draw on the battery :term:`States`: (11) - See BatteryElectroChemEOD, BatteryElectroChemEOL + | tb: Battery temperature (K) + | Vo: Voltage drops due to solid-phase ohmic resistances + | Vsn: Negative surface voltage (V) + | Vsp: Positive surface voltage (V) + | qnB: Amount of negative ions at the battery bulk + | qnS: Amount of negative ions at the battery surface + | qpB: Amount of positive ions at the battery bulk + | qpS: Amount of positive ions at the battery surface + | qMobile: Maximum battery capacity + | tDiffusion : Diffusion time constant (increasing this causes decrease in diffusion rate) + | Ro : Ohmic drop (current collector resistances plus electrolyte resistance plus solid phase resistances at anode and cathode) :term:`Outputs` (2) | t: Temperature of battery (°C) | v: Voltage supplied by battery + Keyword Args + ------------ + process_noise : Optional, float or dict[str, float] + :term:`Process noise` (applied at dx/next_state). + Can be number (e.g., .2) applied to every state, a dictionary of values for each + state (e.g., {'x1': 0.2, 'x2': 0.3}), or a function (x) -> x + process_noise_dist : Optional, str + Distribution for :term:`Process noise` (e.g., normal, uniform, triangular) + measurement_noise : Optional, float or dict[str, float] + :term:`Measurement noise` (applied in output eqn). + Can be number (e.g., .2) applied to every output, a dictionary of values for each + output (e.g., {'z1': 0.2, 'z2': 0.3}), or a function (z) -> z + measurement_noise_dist : Optional, str + Distribution for :term:`measurement noise` (e.g., normal, uniform, triangular) + xnMax : float + Maximum mole fraction (neg electrode) + xpMax : float + Maximum mole fraction (pos electrode). Typically 1. + alpha : float + Anodic/cathodic electrochemical transfer coefficient + Sn : float + Surface area (- electrode) + Sp : float + Surface area (+ electrode) + kn : float + Lumped constant for BV (- electrode) + kp : float + Lumped constant for BV (+ electrode) + Vol : float + Total interior battery volume/2 (for computing concentrations) + VolSFraction : float + Fraction of total volume occupied by surface volume + to : float + For Ohmic voltage + tsn : float + For surface overpotential (neg) + tsp : float + For surface overpotential (pos) + U0p : float + Redlich-Kister parameter (+ electrode) + Ap : float + Redlich-Kister parameter (+ electrode) + U0n : float + Redlich-Kister parameter (- electrode) + An : float + Redlich-Kister parameter (- electrode) + VEOD : float + End of discharge voltage threshold + VDropoff : float + Voltage above EOD after which voltage will be considered in SOC calculation + qMaxThreshold : float + Threshold for qMax (for threshold_met and event_state), after which the InsufficientCapacity event has occurred. Note: Battery manufacturers specify a threshold of 70-80% of qMax + wq : float + Wear rate for qMax + wr : float + Wear rate for Ro + wd : float + Wear rate for D + x0 : dict[str, float] + Initial :term:`state` + See Also -------- BatteryElectroChemEOL, BatteryElectroChemEOD, BatteryElectroChem - - Note - ---- - For keyword arguments, see BatteryElectroChemEOD, BatteryElectroChemEOL """ - inputs = BatteryElectroChemEOD.inputs - outputs = BatteryElectroChemEOD.outputs - states = BatteryElectroChemEOD.states + BatteryElectroChemEOL.states - events = BatteryElectroChemEOD.events + BatteryElectroChemEOL.events + events = ['EOD', 'InsufficientCapacity'] + inputs = ['i'] + states = ['tb', 'Vo', 'Vsn', 'Vsp', 'qnB', 'qnS', 'qpB', 'qpS', 'qMobile', 'tDiffusion', 'Ro'] + outputs = ['t', 'v'] + performance_metric_keys = ['max_i'] is_vectorized = False - default_parameters = deepcopy(BatteryElectroChemEOD.default_parameters) - merge_dicts(default_parameters, - BatteryElectroChemEOL.default_parameters) + state_limits = { + 'tb': (0, np.inf), # Limited by Absolute Zero (0 K) + 'qnB': (0, np.inf), + 'qnS': (0, np.inf), + 'qpB': (0, np.inf), + 'qpS': (0, np.inf), + 'qMobile': (0, np.inf) + } + + param_callbacks = { # Callbacks for derived parameters + 'VolSFraction': [update_vols], + 'Vol': [update_vols], + 'xpMax': [update_xnMin], + 'xnMax': [update_xpMin] + } + + default_parameters = { # Set to defaults + 'xnMax': 0.6, + 'xnMin': 0.0, + 'xpMax': 1.0, + 'xpMin': 0.4, + 'wq': -1e-2, + 'wr': 1e-6, + 'wd': 1e-2, + 'qMaxThreshold': 5320, + + # Li-ion parameters + 'alpha': 0.5, + 'Sn': 0.000437545, + 'Sp': 0.00030962, + 'kn': 2120.96, + 'kp': 248898, + 'Vol': 2e-5, + 'VolSFraction': 0.1, + + # Time constants + 'to': 6.08671, + 'tsn': 1001.38, + 'tsp': 46.4311, + + # Redlich-Kister parameters (+ electrode) + 'U0p': 4.03, + 'Ap': [-31593.7, 0.106747, 24606.4, -78561.9, 13317.9, 307387, 84916.1, -1.07469e+06, 2285.04, 990894, 283920, -161513, -469218], + + # Redlich-Kister parameters (- electrode) + 'U0n': 0.01, + 'An': [86.19, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - state_limits = deepcopy(BatteryElectroChemEOD.state_limits) - state_limits.update(BatteryElectroChemEOL.state_limits) + 'x0': { + 'Vo': 0, + 'Vsn': 0, + 'Vsp': 0, + 'tb': 292.1, # in K, about 18.95 C + 'qMobile': 7600, + 'tDiffusion': 7e6, + 'Ro': 0.117215, + 'qMax': 7600 # Kept in so EOL model will work + }, + + # End of discharge voltage threshold + 'VEOD': 3.0, + 'VDropoff': 0.1 # Voltage above EOD after which voltage will be considered in SOC calculation + } - def __init__(self, **kwargs): - self.param_callbacks['qMobile'].append(OverwrittenWarning) - self.param_callbacks['tDiffusion'] = [OverwrittenWarning] - self.param_callbacks['Ro'] = [OverwrittenWarning] - super().__init__(**kwargs) + def initialize(self, u=None, z=None): + params = self.parameters + x0 = deepcopy(self.parameters['x0']) + qMax = x0['qMobile']/(params['xnMax']-params['xnMin']) + x0['qpS'] = qMax*params['xpMin']*params['VolSFraction'], + x0['qpB'] = qMax*params['xpMin']*(1.0-params['VolSFraction']) + x0['qnS'] = qMax*params['xnMax']*params['VolSFraction'] + x0['qnB'] = qMax*params['xnMax']*(1.0-params['VolSFraction']) + return self.StateContainer(x0) def dx(self, x, u): - # Set EOD Parameters (corresponding to health) - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - self.parameters['qMobile'] = x['qMax'] - self.parameters['Ro'] = x['Ro'] - self.parameters['tDiffusion'] = x['D'] + params = self.parameters + + # EOD model dx + # Negative Surface + CnBulk = x['qnB']/params['VolB'] + CnSurface = x['qnS']/params['VolS'] + qMax = x['qMobile']/(params['xnMax']-params['xnMin']) + qSMax = qMax*params['VolSFraction'] + + xnS = x['qnS']/qSMax + + qdotDiffusionBSn = (CnBulk-CnSurface)/x['tDiffusion'] + qnBdot = -qdotDiffusionBSn + qnSdot = qdotDiffusionBSn - u["i"] + + Jn = u['i']/params['Sn'] + Jn0 = params['kn']*((1-xnS)*xnS)**params['alpha'] + + v_part = R_F*x['tb']/params['alpha'] + + VsnNominal = v_part*np.arcsinh(Jn/(Jn0 + Jn0)) + Vsndot = (VsnNominal-x['Vsn'])/params['tsn'] + + # Positive Surface + CpBulk = x['qpB']/params['VolB'] + CpSurface = x['qpS']/params['VolS'] + xpS = x['qpS']/qSMax - # Calculate - x_dot = BatteryElectroChemEOD.dx(self, x, u) - x_dot2 = BatteryElectroChemEOL.dx(self, x, u) - x_dot.matrix = np.vstack((x_dot.matrix, x_dot2.matrix)) - return x_dot + qdotDiffusionBSp = (CpBulk-CpSurface)/x['tDiffusion'] + qpBdot = -qdotDiffusionBSp + qpSdot = u['i'] + qdotDiffusionBSp - def output(self, x): - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - self.parameters['qMobile'] = x['qMax'] - return BatteryElectroChemEOD.output(self, x) + Jp = u['i']/params['Sp'] + Jp0 = params['kp']*((1-xpS)*xpS)**params['alpha'] + VspNominal = v_part*np.arcsinh(Jp/(Jp0+Jp0)) + Vspdot = (VspNominal-x['Vsp'])/params['tsp'] + + # Combined + VoNominal = u['i']*x['Ro'] + Vodot = (VoNominal-x['Vo'])/params['to'] + + # Thermal Effects + voltage_eta = x['Vo'] + x['Vsn'] + x['Vsp'] # (Vep - Ven) - V; + + Tbdot = voltage_eta*u['i']/mC + (params['x0']['tb'] - x['tb'])/tau # Newman + + # Additional states + Rodot = params['wr'] * abs(u['i']) + + # EOL model dx + qMobiledot = params['wq'] * abs(u['i']) + tDiffusiondot = params['wd'] * abs(u['i']) + + return self.StateContainer(np.array([ + np.atleast_1d(Tbdot), + np.atleast_1d(Vodot), + np.atleast_1d(Vsndot), + np.atleast_1d(Vspdot), + np.atleast_1d(qnBdot), + np.atleast_1d(qnSdot), + np.atleast_1d(qpBdot), + np.atleast_1d(qpSdot), + np.atleast_1d(qMobiledot), + np.atleast_1d(tDiffusiondot), + np.atleast_1d(Rodot), + ])) + def event_state(self, x) -> dict: - e_state = BatteryElectroChemEOD.event_state(self, x) - e_state.update(BatteryElectroChemEOL.event_state(self, x)) + params = self.parameters + + qMax = x['qMobile']/(params['xnMax']-params['xnMin']) + qnMax = qMax*params['xnMax'] + qSMax = qMax*params['VolSFraction'] + + e_state = calculate_EOD(x, params, qnMax, qSMax) + e_state.update(BatteryElectroChemEOL.event_state(self, {'qMax': x['qMobile']})) + return e_state + def output(self, x): + params = self.parameters + + qMax = x['qMobile']/(params['xnMax']-params['xnMin']) + qSMax = qMax*params['VolSFraction'] + + t, v = calculate_temp_voltage(x, params, qSMax) + + return self.OutputContainer(np.array([t, v])) + def threshold_met(self, x) -> dict: t_met = BatteryElectroChemEOD.threshold_met(self, x) - t_met.update(BatteryElectroChemEOL.threshold_met(self, x)) - return t_met + t_met.update(BatteryElectroChemEOL.threshold_met(self, {'qMax': x['qMobile']})) -BatteryElectroChem = BatteryElectroChemEODEOL + return t_met + +BatteryElectroChem = BatteryElectroChemEODEOL \ No newline at end of file diff --git a/tests/test_battery.py b/tests/test_battery.py index d4ea6936..cf08e527 100644 --- a/tests/test_battery.py +++ b/tests/test_battery.py @@ -3,6 +3,7 @@ from io import StringIO import sys import unittest +import numpy as np from progpy.models import BatteryCircuit, BatteryElectroChem, BatteryElectroChemEOL, BatteryElectroChemEOD, BatteryElectroChemEODEOL, SimplifiedBattery from progpy.loading import Piecewise @@ -36,15 +37,107 @@ def test_battery_electrochem(self): result = batt.simulate_to(200, future_loading, {'t': 18.95, 'v': 4.183}) self.assertEqual(BatteryElectroChem, BatteryElectroChemEODEOL) - # check warning raised when changing overwritten parameter - with self.assertWarns(UserWarning): - batt.parameters['Ro'] = 10 + ''' + Test that the current combined model has the same results as the initial implementation. + + Note that this initial implementation is no longer in ProgPy. While we have confidence in + the quantitative results, a development challenge related to different results if print is True or False + forced the need for an updated EODEOL model (see https://github.com/nasa/progpy/issues/199). In the + below test, some values from this initial implementation are hard-coded for testing purposes. + + The combined model changes can be found here: https://github.com/nasa/progpy/pull/207 + Note that states are not compared since the updated implementation uses different states and parameters. + ''' + def test_battery_electrochem_results(self): + config = { + 'save_freq': 1000, + 'dt': 2, + 'events': 'InsufficientCapacity' + } + + def future_loading(t, x=None): + load = 1 + + if x is not None: + event_state = batt.event_state(x) + if event_state["EOD"] > 0.95: + load = 1 # Discharge + elif event_state["EOD"] < 0.05: + load = -1 # Charge + + return batt.InputContainer({'i': load}) + + batt = BatteryElectroChem() + result = batt.simulate_to_threshold(future_loading, **config) + + # Check a middle result + self.assertEqual(result.times[-10], 219000.0) + self.assertEqual(result.inputs[-10], {'i': np.float64(-1.0)}) + self.assertEqual(result.outputs[-10], {'t': np.float64(18.775450787889213), 'v': np.float64(3.0165983707625728)}) + self.assertEqual(result.event_states[-10], {'EOD': np.float64(0.16598370762572756), 'InsufficientCapacity': np.float64(0.03947368418956007)}) + + # Check the last result + self.assertEqual(result.times[-1], 228000.0) + self.assertEqual(result.inputs[-1], {'i': np.float64(-1.0)}) + self.assertEqual(result.outputs[-1], {'t': np.float64(18.76922976507518), 'v': np.float64(3.0087596988706986)}) + self.assertEqual(result.event_states[-1], {'EOD': np.float64(0.08759698870698607), 'InsufficientCapacity': 0.0}) + + ''' + Test that the combined model has the same result if print is True or False. + + This check is necessary since the event_state and output are calculated during the simulation + when print is True but only as needed when print is False due to LazySim optimization. + For more details, refer to https://github.com/nasa/progpy/issues/199 + ''' + def test_battery_electrochem_printed(self): + config = { + 'save_freq': 1000, + 'dt': 2, + 'events': 'InsufficientCapacity', + 'print': True + } + + config2 = { + 'save_freq': 1000, + 'dt': 2, + 'events': 'InsufficientCapacity', + 'print': False + } + + def future_loading(t, x=None): + load = 1 + + if x is not None: + event_state = batt.event_state(x) + if event_state["EOD"] > 0.95: + load = 1 # Discharge + elif event_state["EOD"] < 0.05: + load = -1 # Charge + + return batt.InputContainer({'i': load}) - with self.assertWarns(UserWarning): - batt.parameters['qMobile'] = 10 + def future_loading2(t, x=None): + load = 1 + + if x is not None: + event_state = batt2.event_state(x) + if event_state["EOD"] > 0.95: + load = 1 # Discharge + elif event_state["EOD"] < 0.05: + load = -1 # Charge + + return batt2.InputContainer({'i': load}) + + batt = BatteryElectroChem() + batt2 = BatteryElectroChem() + + result = batt.simulate_to_threshold(future_loading, **config) + result2 = batt2.simulate_to_threshold(future_loading2, **config2) - with self.assertWarns(UserWarning): - batt.parameters['tDiffusion'] = 10 + self.assertEqual(result.times, result2.times) + self.assertEqual(result.states, result2.states) + self.assertEqual(result.event_states, result2.event_states) + self.assertEqual(result.outputs, result2.outputs) def test_battery_electrochem_EOD(self): batt = BatteryElectroChemEOD()