brewpi/components/pid/temp_controller_smith.py

from components.plant.pot import Pot
from matplotlib.pyplot import plot, figure, subplot, grid, show, legend
from components.pid import Kalman
from components.pid.temp_controller_base import TempControllerBase
from components.pid.tc_constants import *
import numpy as np


class TempController(TempControllerBase):
	def __init__(self, dt, params, model_params):
		TempControllerBase.__init__(self, dt, params)
		self.kalman_model = Kalman(dt, params['Kalman'])
		self.kalman_model_delay = Kalman(dt, params['Kalman'])
		self.kalman_plant = Kalman(dt, params['Kalman'])
		self.model = Pot(dt, model_params)

		self.theta_ist_plant = 0
		self.dtheta_ist_plant = 0
		self.theta_ist_model = 0
		self.dtheta_ist_model = 0
		self.theta_ist_model_delay = 0
		self.dtheta_ist_model_delay = 0

	def set_model_power(self, power):
		self.model.set_power(power)

	def init_kalman(self, value):
		self.kalman_model.initial((value, 0))
		self.kalman_model_delay.initial((value, 0))
		self.kalman_plant.initial((value, 0))

	def on_state_entered(self, state):
		if state == States.HEAT:
			self.model.initial(self.theta_ist)
			self.kalman_model.initial((self.theta_ist, 0))
			self.kalman_model_delay.initial((self.theta_ist, 0))

	def post_pid(self):
		self.model.process()

	def process(self):
		# Process Kalman of Plant
		Z_plant = self.kalman_plant.process_measurement((self.theta_ist_set, 0), 0.0)
		xp_plant = self.kalman_plant.process(Z_plant)
		theta_ist_plant = xp_plant[0, 0]
		heatrate_ist_plant = xp_plant[1, 0] * 60

		# Process Kalman of Model
		k_model = self.kalman_model.process_measurement((self.model.get_temperature_intermediate(), 0), 0.0)
		xp_model = self.kalman_model.process(k_model)
		theta_ist_model = xp_model[0, 0]
		heatrate_ist_model = xp_model[1, 0] * 60

		# 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

		self.theta_ist_plant = theta_ist_plant
		self.dtheta_ist_plant = heatrate_ist_plant
		self.theta_ist_model = theta_ist_model
		self.dtheta_ist_model = heatrate_ist_model
		self.theta_ist_model_delay = theta_ist_model_delay
		self.dtheta_ist_model_delay = dtheta_ist_model_delay
		self.theta_ist = theta_ist_plant
		self.heatrate_ist = heatrate_ist_plant

		# Compensate for max heat rate to reduce overshoot
		if self.heatrate_soll_set > 0:
			self.pid_hold.scale(1.0/self.heatrate_soll_set)

		self.heatrate_soll = self.heatrate_soll_set * self.pid_hold.get_y()

		if 1:
			theta_err = self.theta_soll_set - (theta_ist_plant - theta_ist_model_delay + theta_ist_model)
		else:
			theta_err = self.theta_soll_set - theta_ist_plant

		heatrate_err = self.heatrate_soll - (heatrate_ist_plant - dtheta_ist_model_delay + heatrate_ist_model)
		diff = self.theta_soll_set - self.theta_ist

		self.process_fsm(diff)
		self.process_pid(theta_err, heatrate_err)


if __name__ == '__main__':
	dt = 1.0

	temp_soll = 20
	ctrl = TempController(dt, Test.tc_ctrl_params, Test.tc_model_params)
	plant = Pot(dt, Test.tc_pot_params)
	_y = np.empty(0)
	_fb = np.empty(0)
	_t = np.empty(0)
	_temp_soll = np.empty(0)
	_temp_ist_kalman = np.empty(0)
	_heatrate_ist_kalman = np.empty(0)
	_temp_ist_kalman_plant = np.empty(0)
	_heatrate_ist_kalman_plant = np.empty(0)
	_temp_ist_kalman_model = np.empty(0)
	_heatrate_ist_kalman_model = np.empty(0)
	a = 0.5
	fb = 0
	rho = 0.02
	temps = [{'Temp': 20, 'Duration': 1000}, {'Temp': 40, 'Duration': 1000}, {'Temp': 50, 'Duration': 1000}, {'Temp': 60, 'Duration': 1000}, {'Temp': 70, 'Duration': 1000}, {'Temp': 80, 'Duration': 1000}, {'Temp': 78, 'Duration': 1000}]

	t = 0

	for temp in temps:
		temp_soll = temp['Temp']
		hold_counter = temp['Duration']
		hold = False
		ctrl.set_theta_soll(temp_soll)
		ctrl.set_heatrate_soll(1.0)

		while True:
			if hold:
				if hold_counter == 0:
					break
				hold_counter -= 1
			plant.process()
			temp_ist = plant.get_temperature()
			ctrl.set_theta_ist(temp_ist)
			ctrl.process()

			y = 3500*ctrl.get_power()
			power = max(0, y)
			plant.set_power(power)
			fb = plant.get_power()
			if abs(temp_ist - temp_soll) < 0.1:
				hold = True

			_y = np.append(_y, y)
			_fb = np.append(_fb, fb)
			_t = np.append(_t, t)
			_temp_soll = np.append(_temp_soll, temp_soll)
			_temp_ist_kalman_plant = np.append(_temp_ist_kalman_plant, ctrl.theta_ist_plant)
			_heatrate_ist_kalman_plant = np.append(_heatrate_ist_kalman_plant, max(-1, min(3, ctrl.dtheta_ist_plant)))
			_temp_ist_kalman_model = np.append(_temp_ist_kalman_model, ctrl.theta_ist_model_delay)
			_heatrate_ist_kalman_model = np.append(_heatrate_ist_kalman_model, max(-1, min(3, ctrl.dtheta_ist_model_delay)))
			t += 1

	figure(1)
	subplot(3, 1, 1)
	plot(_t, _temp_ist_kalman_plant, _t, _temp_soll, 'r-', linewidth=1)
	legend(["ist", "soll"])
	grid(True)
	subplot(3, 1, 2)
	plot(_t, _y, '-b', _t, _fb, '-r', linewidth=1)
	legend(["y", "pot"])
	grid(True)
	subplot(3, 1, 3)
	plot(_t, _heatrate_ist_kalman_plant, '-b', linewidth=1)
	legend(["heatrate"])
	grid(True)

	figure(2)
	subplot(2, 1, 1)
	plot(_t, _temp_ist_kalman_plant, '-b', _t, _temp_ist_kalman_model, 'r-', linewidth=1)
	legend(["plant", "model"])
	grid(True)
	subplot(2, 1, 2)
	plot(_t, _heatrate_ist_kalman_plant, '-b', _t, _heatrate_ist_kalman_model, '-r', linewidth=1)
	legend(["plant", "model"])
	grid(True)
	show()

	print("End of program")