Add files from Micah Black on SoC refactor

displays/SOC_velocity_graph.py

@@ -4,7 +4,7 @@
 sys.path.append(os.path.join(os.path.dirname(__file__), '..'))
 
 from optimization.car_model import Car
-from soc.SoCEstimation import CoulombCounter
+from soc.soc_deprecated.SoCEstimation import CoulombCounter
 
 def calculate_SOC_values(v_profile, e_profile, distance, initial_soc, min_speed=None, max_speed=None):
     if min_speed is not None and max_speed is not None:

displays/tests/test_SOC_velocity_graph.py

@@ -4,8 +4,6 @@
 import pytest
 import mock
 from displays.SOC_velocity_graph import calculate_SOC_values
-from matplotlib.testing.decorators import check_figures_equal
-import matplotlib.pyplot as plt
 
 def test_calculate_SOC_values_zero_values():
     velocities = [0, 0, 0]
@@ -41,7 +39,7 @@ def test_calculate_SOC_values_zero_length_inputs():
         calculate_SOC_values([], [], [], 1)
 
 def test_calculate_SOC_values_single_element_inputs():
-    with mock.patch("soc.SoCEstimation.CoulombCounter.get_soc", return_value=0.9):
+    with mock.patch("soc.soc_deprecated.SoCEstimation.CoulombCounter.get_soc", return_value=0.9):
         expected_result = [1, 0.9]
         assert(expected_result == (calculate_SOC_values([1], [(1, 1)], [1], 1)))
 

soc/BatTestPlot.py

@@ -0,0 +1,74 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import itertools
+from sklearn import datasets, linear_model, metrics
+from sklearn.preprocessing import PolynomialFeatures
+from sklearn.pipeline import make_pipeline
+
+#fName = "D:\Documents\MSXIV\MSXIV\Strategy\CellDataFileTestMJ1.txt"
+fName = r"C:\Users\micke\Documents\MSXIV\Strategy\CellDataFileTestMJ1.txt"
+
+
+class BatTestPlot:
+	def __init__(self):
+		#includes 0.050 cell and 0.035 cable, from tests done on single cell
+		self.IR = 0.085
+
+	def plot_SoCOCV(self):
+		print ("plotting SoC OCV")
+
+		x = np.arange(0,10,0.2)
+		y = np.sin(x)
+
+		socfig, ax = plt.subplots()
+		ax.plot(x,y)
+		plt.title('MJ1 SoC OCV @3.5A')
+		plt.show()
+
+	def read_data(self, filename):
+		#voltage, capacity = np.loadtxt(filename, delimiter = '\t', skiprows = 1, usecols = (1,7), unpack = True)
+		with open(filename) as f_in:
+			voltage, current, soc = np.loadtxt(itertools.islice(f_in, 1, None, 2), delimiter = '\t', usecols = (1,3,9), unpack = True)
+		v_ocv = voltage + self.IR*current
+
+		fig, ax = plt.subplots()
+
+		ax.plot(soc, voltage, 'bo',label = 'v')
+		ax.plot(soc, v_ocv, 'm-', label = 'ocv')
+
+		soc = soc.reshape(-1,1)
+		voltage = voltage.reshape(-1,1)
+
+		#model = linear_model.LinearRegression()
+		#model.fit(soc, voltage)
+
+		#voltage_test = model.predict(soc)
+
+		#ax.plot(soc, voltage_test, 'r-')
+
+		#print('Coefficients: \n', model.coef_)
+		#print('Variance Score: %.2f' % metrics.r2_score(voltage, voltage_test))
+
+		poly_model = make_pipeline(PolynomialFeatures(5), linear_model.Ridge())
+		poly_model.fit(soc, v_ocv)
+
+		print ("coefficients:")
+		print (poly_model[1].coef_)
+		print ("Intercept:")
+		print (poly_model[1].intercept_)
+
+		poly_voltage_plot = poly_model.predict(soc)
+		ax.plot(soc, poly_voltage_plot, 'g-', label = 'v_prediction', linewidth = 6)
+
+		plt.title('Soc vs Voltage LG MJ1 @3A DSC')
+		plt.legend(loc = 'lower right')
+		plt.grid(True)
+		plt.show()
+
+
+
+print ("Started")
+
+test = BatTestPlot()
+test.read_data(fName)
+#test.plot_SoCOCV()

soc/PackEfficiency.py

@@ -1,6 +1,6 @@
-#this is programmed in notepad++ and not tested.
 
 import cmath
+import SoC_OCV
 
 class PackEfficiency:
 	'''
@@ -18,27 +18,112 @@ class PackEfficiency:
 	power loss in a pack is due primarily to ohmic heating
 	'''
 
-
 	def __init__(self):
-		self._pack_resistance = 0.055
+		self.pack_resistance = 0.055
+
+		#initialize SoC-OCV model curve
+		self.soc_ocv_curve = SoC_OCV.SoC_OCV()
 
 	#if we want say 1000W at the motors, we will need to draw more power (1050W) from the cells
-	def draw_power(self, power_outside_pack_W):
-		pack_ocv = get_pack_ocv() #from SoC-OCV curve or telemetry, assume cells are balanced
+	def draw_power(self, power_outside_pack_W, soc, max = 100): #max is the maximum value for 100% SoC
+		if(power_outside_pack_W == 0):
+			return (0, 1)
+
+		#pack_ocv = 36*4 #from SoC-OCV curve or telemetry, assume cells are balanced
+		pack_ocv = self.soc_ocv_curve.get_cell_ocv(soc, max) * 36
+		print("Pack OCV: {}".format(pack_ocv))
 
 		#solve quadratic equation to find current
 		#I^2Ri - IVo + Po = 0
-		current_draw = (pack_ocv - cmath.sqrt((pack_ocv ** 2) - (4 * self._pack_resistance * power_outside_pack_W))) / (2 * (self._pack_resistance))
+		current_draw = (pack_ocv - cmath.sqrt((pack_ocv**2) - (4 * self.pack_resistance * power_outside_pack_W))) / (2*(self.pack_resistance))
+		print("Current Draw: {}".format(current_draw))
+
 		#this power loss is wasted as heat through the resistive elements
-		power_loss = (current_draw ** 2) * self._pack_resistance
+		power_loss = (current_draw**2)*self.pack_resistance
 
 		#the power drawn from the cell itself (before resistive losses from internal resistance), the actual power draw for SoC estimation
 		power_inside_pack = power_outside_pack_W + power_loss
 
 		#1 = 100% efficiency
-		efficiency = (power_inside_pack - power_loss) / power_inside_pack
+		#discharging efficiency - useable energy / given energy
+		if(power_outside_pack_W > 0):
+			efficiency = (power_inside_pack - power_loss) / power_inside_pack
+		#charging efficiency - power loss will be negative in this case - useable energy / given energy
+		else:
+			efficiency = (power_inside_pack) / power_outside_pack_W
+
+		return (power_inside_pack, efficiency)
+
 
-		return power_inside_pack
+#Make sure we get values that make sense
+class Test_Pack_Efficiency:
+	def __init__(self):
+		self.PE = PackEfficiency()
+	def test(self):
+		power, efficiency = self.PE.draw_power(100, 100)
+		print("100W Test: {} Eff: {}".format(power, efficiency))
+		power, efficiency = self.PE.draw_power(500, 90)
+		print("500W Test: {} Eff: {}".format(power, efficiency))
+		power, efficiency = self.PE.draw_power(5000, 80)
+		print("5000W Test: {} Eff: {}".format(power, efficiency))
+		power, efficiency = self.PE.draw_power(20000, 10)
+		print("20000W Test: {} Eff: {}".format(power, efficiency))
+		power, efficiency = self.PE.draw_power(0, 10)
+		print("0W Test: {} Eff: {}".format(power, efficiency))
+		power, efficiency = self.PE.draw_power(-500, 50)
+		print("-500W Test: {} Eff: {}".format(power, efficiency))
+		power, efficiency = self.PE.draw_power(-5000, 80)
+		print("-5000W Test: {} Eff: {}".format(power, efficiency))
+
+		power, efficiency = self.PE.draw_power(20000, 100)
+		print("\n\n20000W  100% Test: {} Eff: {}".format(power, efficiency))
+		power, efficiency = self.PE.draw_power(20000, 50)
+		print("20000W 50% Test: {} Eff: {}".format(power, efficiency))
+		power, efficiency = self.PE.draw_power(20000, 0)
+		print("20000W 0% Test: {} Eff: {}".format(power, efficiency))
+
+
+
 
-def get_pack_ocv():
-		return 0.0
+'''
+#Testing Code	
+print ("Testing Started")
+pe = Test_Pack_Efficiency()
+pe.test()
+'''
+
+if __name__ == "__main__":
+	#plot efficiency and power loss vs power draw
+	#make sure we get reasonable values
+	import matplotlib.pyplot as plt
+	import numpy as np
+
+	pet = Test_Pack_Efficiency()
+	pet.test()
+
+	pe = PackEfficiency()
+
+	p_array = []
+	e_array = []
+	w_array = []
+
+	for i in range(0, 20000, 1000):
+		power, efficiency = pe.draw_power(i, 60)
+		p_array.append(power)
+		e_array.append(efficiency)
+		w_array.append(power - i)
+
+	print(p_array)
+	print(e_array)
+
+	fig, ax = plt.subplots()
+	ax.plot(p_array, e_array, 'b-',label = '%e')
+	plt.ylabel('Efficiency')
+	plt.xlabel('Power Draw')
+
+	ax2 = ax.twinx()
+	ax2.set_ylabel('Power Loss')
+	ax2.plot(p_array, w_array)
+
+	plt.show()
+

soc/SoCEstimation.py

@@ -1,11 +1,11 @@
-#Cell SOC Estimation class - Coulomb counting by energy (Amp draw over telemetry)
+#Cell SOC Estimation class - Coulomb counting by pack data available
 #Based off Karl's coulomb counting method
 #https://uwmidsun.atlassian.net/wiki/spaces/~karlding/pages/620036100/BMS+SOC+Algorithm+v1.0.0+aka.+We+live+in+a+SOCiety
-#not really a coulomb counter, more of an energy meter.
 
-#how is this being accessed, etc??
-#Will it be different when used on telemetry data / modelling
-#this is programmed in notepad++ and not tested.
+
+#This class operates on the assumptions that our pack is perfectly balanced, and there is no evident capacity fade
+#SoC returned is between 0 (0%) and 1 (100%)
+import PackEfficiency as PE
 
 class CoulombCounter:
 	'''
@@ -15,37 +15,81 @@ class CoulombCounter:
 	3.635V nominal
 	'''
 	def __init__(self):
-		self._energy_available = 0
-		self._pack_energy = 36 * 36 * 3.45 * 3.635
-		self._SoC = self._energy_available / self._pack_energy
+		self.energy_available = 0
+		self.pack_energy = 36*36*3.45*3.635 #36S * 36P * 3.45Ah * 3.635V = 16 252.812 Wh
+		self.SoC = self.energy_available / self.pack_energy
+		self.maxSoC = 1 #the value of self.SoC when at 100% full
+		self.PackEff = PE.PackEfficiency()
+		#cell ocv is available with PackEff.soc_ocv_curve.get_cell_ocv(self.SoC, max = 1)
 
 	def set_SoC(self, soc):
-		self._SoC = soc
-		self._energy_available = soc * self._pack_energy
+		self.SoC = soc
+		self.energy_available = soc * self.pack_energy
 
 	def set_energy(self, energy):
-		self._energy_available = energy
-		self._SoC = self._energy_available / self._pack_energy
+		self.energy_available = energy
+		self.SoC = self.energy_available / self.pack_energy
 
 	def get_soc(self):
-		return self._SoC
+		return self.SoC
+
+	def _print(self):
+		print("SoC: {}  Energy Remaining: {}  Voltage: {}".format(self.SoC, self.energy_available, self.PackEff.soc_ocv_curve.get_cell_ocv(self.SoC, max = self.maxSoC)*36))
 
 	#given telemetry updates, discharge the cells
 	def telemetry_discharge(self, telemetry_pack_current, telemetry_pack_voltage, telemetry_interval_ms):
 		self.discharge(self, telemetry_pack_current * telemetry_pack_voltage, telemetry_interval_ms*1000)
 
 	#use a certain amount of power for a certain amount of time
 	#dirOUT true for discharge, false for charge
-	def discharge(self, power_W, time_S, dirOUT):
-		#total energy in kWh
-		energy = power_W * time_S / 3600 / 1000
-
-		if dirOUT:
-			#discharge
-			self._energy_available -= energy
-		else:
-			#charge
-			self._energy_available += energy
+	def discharge(self, power_W, time_S, dirOUT = True):
+		if not dirOUT:
+			power_W *= -1
+
+		#power calculation taking the power loss into account
+		power_total, efficiency = self.PackEff.draw_power(power_W, self.SoC, max = self.maxSoC)
+		print(power_total)
+
+		#total energy in Wh
+		energy = power_total * time_S/3600
+
+		#discharge
+		self.energy_available -= energy
 
 		#update SoC
-		self._SoC = self._energy_available / self._pack_energy
+		self.SoC = self.energy_available / self.pack_energy
+
+
+
+class TestCoulombCounter:
+	def __init__(self):
+		self.CC = CoulombCounter()
+		self.CC._print()
+
+	def test(self):
+		self.CC.set_SoC(1)
+		self.CC._print()
+		self.CC.set_SoC(0.75)
+		self.CC._print()
+
+	def test2(self):
+		self.CC.set_SoC(0.9)
+		self.CC._print()
+		for x in range(4):
+			self.CC._print()
+			self.CC.discharge(20000, 500)
+		self.CC._print()
+
+		#expected without other power loss
+		print("Expected without Power loss")
+		print(0.9 * self.CC.pack_energy - 20000 * 500 / 3600 * 4)
+
+#next step: plot this against another discharge curve
+if __name__ == "__main__":
+	#print out to console and make sure we get reasonable values
+	print ("Starting")
+	TCC = TestCoulombCounter()
+	print("\n\n Test1: \n")
+	TCC.test()
+	print("\n\n Test 2: \n")
+	TCC.test2()