- use heat diffusion for water model instead of delay line

components/pid/temp_controller_smith.py

@@ -1,4 +1,4 @@
-from components.plant.pot import Pot, Delay
+from components.plant.pot import Pot, HeatDiffusion
 from matplotlib.pyplot import plot, figure, subplot, grid, show, legend
 from components.pid import Pid, Kalman
 from components import APid
@@ -28,6 +28,7 @@ class TempController(APid):
 		self.params = params
 		self.model_params = model_params
 		self.kalman_model = Kalman(dt, params['Kalman'])
+		self.kalman_model_delay = Kalman(dt, params['Kalman'])
 		self.kalman_plant = Kalman(dt, params['Kalman'])
 		self.y = -1
 		self.state = States.IDLE
@@ -35,8 +36,6 @@ class TempController(APid):
 		self.pid_hold.set_params(params['Hold'])
 		self.pid_rate.set_params(params['Heat'])
 		self.model = Pot(dt, model_params)
-		self.theta_ist_delay_model = Delay(dt, model_params['Td'], 0)
-		self.dtheta_ist_delay_model = Delay(dt, model_params['Td'], 0)
 		self.is_startup = True
 
 	def set_theta_ist(self, value):
@@ -44,6 +43,7 @@ class TempController(APid):
 		if self.is_startup:
 			self.is_startup = False
 			self.kalman_model.initial((value, 0))
+			self.kalman_model_delay.initial((value, 0))
 			self.kalman_plant.initial((value, 0))
 
 	def get_theta_ist(self):
@@ -75,12 +75,16 @@ class TempController(APid):
 
 	def process(self):
 		# Process Kalman of Model
-		Z_model = self.kalman_model.process_measurement((self.model.get_temperature_intermediate(), 0), 0.0)
+		k_model = self.kalman_model.process_measurement((self.model.get_temperature_intermediate(), 0), 0.0)
-		xp_model = self.kalman_model.process(Z_model)
+		xp_model = self.kalman_model.process(k_model)
 		theta_ist_model = xp_model[0, 0]
 		heatrate_ist_model = xp_model[1, 0] * 60
-		self.theta_ist_delay_model.put(theta_ist_model)
+
-		self.dtheta_ist_delay_model.put(heatrate_ist_model)
+		# Process Kalman of delayed Model
+		k_model_delay = self.kalman_model_delay.process_measurement((self.model.get_temperature(), 0), 0.0)
+		xp_model_delay = self.kalman_model_delay.process(k_model_delay)
+		theta_ist_model_delay = xp_model_delay[0, 0]
+		dtheta_ist_model_delay = xp_model_delay[1, 0] * 60
 
 		# Process Kalman of Plant
 		Z_plant = self.kalman_plant.process_measurement((self.theta_ist_set, 0), 0.0)
@@ -97,8 +101,6 @@ class TempController(APid):
 		self.theta_ist = theta_ist_plant
 		self.heatrate_ist = heatrate_ist_plant
 
-		theta_ist_model_delay = self.theta_ist_delay_model.get()
-		dtheta_ist_model_delay = self.dtheta_ist_delay_model.get()
 		print ("Model    : T_ist={:2.2f}, dT_ist={:2.2f}".format(theta_ist_model, heatrate_ist_model))
 		print ("Model*z-1: T_ist={:2.2f}, dT_ist={:2.2f}".format(theta_ist_model_delay, dtheta_ist_model_delay))
 		print ("Plant    : T_ist={:2.2f}, dT_ist={:2.2f}".format(theta_ist_plant, heatrate_ist_plant))

components/plant/heat_diffusion.py

@@ -0,0 +1,23 @@
+import numpy as np
+
+
+class HeatDiffusion:
+	def __init__(self, dt, delay, ic=0):
+		self.alpha = dt/delay
+		self.beta = 1 - self.alpha
+		self.state = ic
+
+	def initial(self, ic):
+		self.state = ic
+
+	def process(self, x):
+		self.state = self.beta * self.state + self.alpha*x
+		return self.state
+
+
+if __name__ == '__main__':
+	d = HeatDiffusion(0.1, 1)
+
+	for i in range(1, 200):
+		y = d.process(1)
+		print("In = {}, out = {}".format(i, y))

components/plant/pot.py

@@ -1,5 +1,5 @@
 from components.aplant import APlant
-from components.plant.delay import Delay
+from components.plant.heat_diffusion import HeatDiffusion
 
 
 class Pot(APlant):
@@ -22,7 +22,7 @@ class Pot(APlant):
 		self.theta_amb = theta_amb
 
 		self.power_set = 0
-		self.delay = Delay(self.dt, self.Td, params['theta'])
+		self.delay = HeatDiffusion(self.dt, self.Td, params['theta'])
 
 	def initial(self, temp):
 		self.temp_intermediate = temp
@@ -33,12 +33,12 @@ class Pot(APlant):
 		pass
 
 	def process(self):
+		power = self.gain * self.power_set
 		leak = 1/self.M * self.L * (self.theta_amb - self.temp)/self.theta_amb
-		self.temp_intermediate += (self.power_set / (self.M * self.C) + leak) * self.dt
+		self.temp_intermediate += (power / (self.M * self.C) + leak) * self.dt
 
 		# Delay
-		self.delay.put(self.temp_intermediate)
+		self.temp = self.delay.process(self.temp_intermediate)
-		self.temp = self.delay.get()
 
 	def is_activated(self):
 		return True