From 631e49fb6b7cdaa834277c6eddaecbaec9e8e414 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Isak=20=C3=85slund?= Date: Thu, 21 Nov 2024 11:13:23 +0100 Subject: [PATCH 1/4] feat: Add example for fast system with PI control --- ..._fast_system_with_PI_room+floor_control.py | 145 ++++++++++++++++++ src/thermal_sim/thermal_sim.py | 1 + 2 files changed, 146 insertions(+) create mode 100644 example_fast_system_with_PI_room+floor_control.py diff --git a/example_fast_system_with_PI_room+floor_control.py b/example_fast_system_with_PI_room+floor_control.py new file mode 100644 index 0000000..8ec2875 --- /dev/null +++ b/example_fast_system_with_PI_room+floor_control.py @@ -0,0 +1,145 @@ +from math import sin, pi +from src.thermal_sim.thermal_sim import ThermalMass, ThermalResistance, ThermalSystem, PowerSource, TemperatureSource, Solver + +def clip(n, min, max): + sat = n + if n < min: + sat = min + elif n > max: + sat = max + + return (sat, n != sat) + + +def PI(e, integral, Kp, Ki, dt, feedforward): + prev_integral = integral + out_min = 15 + out_max = 30 + + int_min = -2 + int_max = 9 + + # Proportional term + P = Kp * e + + # Update integral term + integral += e * dt + + # Calculate the integral contribution and clamp is within bounds + I = Ki * integral + integral_saturated = False + if I > int_max: + I = int_max + integral = int_max / Ki + integral_saturated = True + elif I < int_min: + I = int_min + integral = int_min / Ki + integral_saturated = True + + # Calculate the output and check for saturation + setpoint_ref_raw = P + I + feedforward + setpoint_ref_sat, output_saturated = clip(setpoint_ref_raw, out_min, out_max) + + # Prevent integral windup + if output_saturated and not integral_saturated: + integral = prev_integral + + return setpoint_ref_sat, integral + + + + +def bangbang(T, T_ref, T_hyst, prev_u): + T_max = T_ref + T_hyst + T_min = T_ref - T_hyst + + if T > T_max: + u = 0 + elif T < T_min: + u = 1 + else: + u = prev_u + + return u + +def get_temp_by_name(system: ThermalSystem, name: str): + for mass in system.thermal_masses: + if mass.__name__ == name: + return mass.T + return None + + + +if __name__ == "__main__": + # Simulate floor and air connected model with air and floor feedback (PI + Bangbang) + + T0 = 17 + GND = ThermalMass.get_ground(T=T0) + + floor_name = "Floor" + air_name = "Air" + + def controller(t: float, prev_Q: float, system: ThermalSystem) -> float: + if not hasattr(controller, "integral"): + controller.integral = 0 + controller.t_prev = 0 + + P = 2000 + u_prev = 1 if prev_Q > 0 else 0 + + # TRM can only run at 1s update rate so make sure we do that + dt_last_run = t - controller.t_prev + if (dt_last_run) < 1: + return u_prev * P + + T_air = get_temp_by_name(system, air_name) + T_floor = get_temp_by_name(system, floor_name) + T_hyst = 2 + T_ref = 24 + if t > 1800: + T_ref = 20 + + + e = T_ref - T_air + Kp = 3 + Ti = 100 # Time constant in seconds + Ki = Kp / Ti + + + setpoint_floor, integral = PI(e, controller.integral, Kp, Ki, dt_last_run, T_ref) + controller.integral = integral + u = bangbang(T_floor, setpoint_floor, T_hyst, u_prev) + + controller.t_prev = t + return u * P + + + + # Setup physcal system elements + m_floor = ThermalMass(C=1000*4, T_initial=T0, name=floor_name) + m_air = ThermalMass(C=1000*8, T_initial=T0, name=air_name) + + R_floor_gnd = ThermalResistance(R=0.01, name="R_f_gnd") + R_air_gnd = ThermalResistance(R=0.1, name="R_a_gnd") + R_floor_air = ThermalResistance(R=0.05, name="R_f_a") + R_roof_gnd = ThermalResistance(R=0.01, name="R_r_gnd") + R_air_roof = ThermalResistance(R=1, name="R_a_r") + + heat_source = PowerSource(controller, name="Heating cable") + + # Connect elements + m_floor.connect(heat_source) + R_floor_gnd.connect(m_floor, GND) + + # m_air.connect(sun_source) + R_air_gnd.connect(m_air, GND) + R_floor_air.connect(m_floor, m_air) + + # Simulate + thermal_masses = [m_floor, m_air] + system = ThermalSystem(thermal_masses) + system.generate_diagram('img/multi_mass.png') + dt = 1 # [s] + total_time = 3600 # [s] + system.simulate(dt, total_time, solver=Solver.RK4) \ No newline at end of file diff --git a/src/thermal_sim/thermal_sim.py b/src/thermal_sim/thermal_sim.py index f900501..b34ed3b 100644 --- a/src/thermal_sim/thermal_sim.py +++ b/src/thermal_sim/thermal_sim.py @@ -175,6 +175,7 @@ def custom_ode_solver(self, y0, t:float, solver:Solver): return solution def simulate(self, dt:float, duration:float, plot:bool=True, solver:Solver=Solver.EULER): + self.dt = dt time = np.arange(0, duration, dt) y0 = [mass.T for mass in self.thermal_masses] From 5998d8cae665880a5e73d7c5bda6b8e9efc646d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Isak=20=C3=85slund?= Date: Thu, 21 Nov 2024 15:15:32 +0100 Subject: [PATCH 2/4] feat: PI+PWM (room control) fast simulaton --- room_only.py | 150 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 150 insertions(+) create mode 100644 room_only.py diff --git a/room_only.py b/room_only.py new file mode 100644 index 0000000..34231e4 --- /dev/null +++ b/room_only.py @@ -0,0 +1,150 @@ +from math import sin, pi +from src.thermal_sim.thermal_sim import ThermalMass, ThermalResistance, ThermalSystem, PowerSource, TemperatureSource, Solver + +def clip(n, min, max): + sat = n + if n < min: + sat = min + elif n > max: + sat = max + + return (sat, n != sat) + + +def PI(e, integral, Kp, Ki, dt, feedforward): + prev_integral = integral + out_min = 0 + out_max = 1 + + int_min = 0 + int_max = 1 + + # Proportional term + P = Kp * e + + # Update integral term + integral += e * dt + + # Calculate the integral contribution and clamp is within bounds + I = Ki * integral + integral_saturated = False + if I > int_max: + I = int_max + integral = int_max / Ki + integral_saturated = True + elif I < int_min: + I = int_min + integral = int_min / Ki + integral_saturated = True + + # Calculate the output and check for saturation + setpoint_ref_raw = P + I + feedforward + setpoint_ref_sat, output_saturated = clip(setpoint_ref_raw, out_min, out_max) + + # Prevent integral windup + if output_saturated and not integral_saturated: + integral = prev_integral + + return setpoint_ref_sat, integral + + + + +def bangbang(T, T_ref, T_hyst, prev_u): + T_max = T_ref + T_hyst + T_min = T_ref - T_hyst + + if T > T_max: + u = 0 + elif T < T_min: + u = 1 + else: + u = prev_u + + return u + +def pwm(t: float, duty:float, period: float) -> float: + return 1 if (t % period) < duty*period else 0 + +def get_temp_by_name(system: ThermalSystem, name: str): + for mass in system.thermal_masses: + if mass.__name__ == name: + return mass.T + return None + + + +if __name__ == "__main__": + # Simulate floor and air connected model with air and floor feedback (PI + Bangbang) + + T0 = 17 + GND = ThermalMass.get_ground(T=T0) + + floor_name = "Floor" + air_name = "Air" + + def controller(t: float, prev_Q: float, system: ThermalSystem) -> float: + if not hasattr(controller, "integral"): + controller.integral = 0 + controller.t_prev = 0 + + P = 2000 + u_prev = 1 if prev_Q > 0 else 0 + + # TRM can only run at 1s update rate so make sure we do that + dt_last_run = t - controller.t_prev + if (dt_last_run) < 1: + return u_prev * P + + T_air = get_temp_by_name(system, air_name) + T_ref = 24 + if t > 3600*2: + T_ref = 20 + + + e = T_ref - T_air + Kp = 1/3 + Ti = 300*4 # Time constant in seconds + Ki = Kp / Ti + + # Kp = 1/3 + # Ki = 0 + pwm_period = 5*60 + + + setpoint_duty, integral = PI(e, controller.integral, Kp, Ki, dt_last_run, 0) + controller.integral = integral + u = pwm(t, duty=setpoint_duty, period=pwm_period) + + controller.t_prev = t + return u * P + + + + # Setup physcal system elements + m_floor = ThermalMass(C=1000*4*10, T_initial=T0, name=floor_name) + m_air = ThermalMass(C=1000*8*10, T_initial=T0, name=air_name) + + R_floor_gnd = ThermalResistance(R=0.01, name="R_f_gnd") + R_air_gnd = ThermalResistance(R=0.1, name="R_a_gnd") + R_floor_air = ThermalResistance(R=0.05, name="R_f_a") + R_roof_gnd = ThermalResistance(R=0.01, name="R_r_gnd") + R_air_roof = ThermalResistance(R=1, name="R_a_r") + + heat_source = PowerSource(controller, name="Heating cable") + + # Connect elements + m_floor.connect(heat_source) + R_floor_gnd.connect(m_floor, GND) + + # m_air.connect(sun_source) + R_air_gnd.connect(m_air, GND) + R_floor_air.connect(m_floor, m_air) + + # Simulate + thermal_masses = [m_floor, m_air] + system = ThermalSystem(thermal_masses) + system.generate_diagram('img/multi_mass.png') + dt = 1 # [s] + total_time = 3600*4 # [s] + system.simulate(dt, total_time, solver=Solver.RK4) \ No newline at end of file From 32e482282002f17b3c632d65e81f3984d561277b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Isak=20=C3=85slund?= Date: Thu, 21 Nov 2024 15:21:27 +0100 Subject: [PATCH 3/4] feat: PI+PWM with faster system 5min period --- ...le_fast_system_with_PI_room_pwm_control.py | 150 ++++++++++++++++++ 1 file changed, 150 insertions(+) create mode 100644 example_fast_system_with_PI_room_pwm_control.py diff --git a/example_fast_system_with_PI_room_pwm_control.py b/example_fast_system_with_PI_room_pwm_control.py new file mode 100644 index 0000000..fd231d2 --- /dev/null +++ b/example_fast_system_with_PI_room_pwm_control.py @@ -0,0 +1,150 @@ +from math import sin, pi +from src.thermal_sim.thermal_sim import ThermalMass, ThermalResistance, ThermalSystem, PowerSource, TemperatureSource, Solver + +def clip(n, min, max): + sat = n + if n < min: + sat = min + elif n > max: + sat = max + + return (sat, n != sat) + + +def PI(e, integral, Kp, Ki, dt, feedforward): + prev_integral = integral + out_min = 0 + out_max = 1 + + int_min = 0 + int_max = 1 + + # Proportional term + P = Kp * e + + # Update integral term + integral += e * dt + + # Calculate the integral contribution and clamp is within bounds + I = Ki * integral + integral_saturated = False + if I > int_max: + I = int_max + integral = int_max / Ki + integral_saturated = True + elif I < int_min: + I = int_min + integral = int_min / Ki + integral_saturated = True + + # Calculate the output and check for saturation + setpoint_ref_raw = P + I + feedforward + setpoint_ref_sat, output_saturated = clip(setpoint_ref_raw, out_min, out_max) + + # Prevent integral windup + if output_saturated and not integral_saturated: + integral = prev_integral + + return setpoint_ref_sat, integral + + + + +def bangbang(T, T_ref, T_hyst, prev_u): + T_max = T_ref + T_hyst + T_min = T_ref - T_hyst + + if T > T_max: + u = 0 + elif T < T_min: + u = 1 + else: + u = prev_u + + return u + +def pwm(t: float, duty:float, period: float) -> float: + return 1 if (t % period) < duty*period else 0 + +def get_temp_by_name(system: ThermalSystem, name: str): + for mass in system.thermal_masses: + if mass.__name__ == name: + return mass.T + return None + + + +if __name__ == "__main__": + # Simulate floor and air connected model with air and floor feedback (PI + Bangbang) + + T0 = 19 + GND = ThermalMass.get_ground(T=T0) + + floor_name = "Floor" + air_name = "Air" + + def controller(t: float, prev_Q: float, system: ThermalSystem) -> float: + if not hasattr(controller, "integral"): + controller.integral = 0 + controller.t_prev = 0 + + P = 2000 + u_prev = 1 if prev_Q > 0 else 0 + + # TRM can only run at 1s update rate so make sure we do that + dt_last_run = t - controller.t_prev + if (dt_last_run) < 1: + return u_prev * P + + T_air = get_temp_by_name(system, air_name) + T_ref = 24 + if t > system.t_max/2: + T_ref = 20 + + + e = T_ref - T_air + Kp = 1/3 + Ti = 300*4 # Time constant in seconds + Ki = Kp / Ti + + # Kp = 1/3 + # Ki = 0 + pwm_period = 5*60 + + + setpoint_duty, integral = PI(e, controller.integral, Kp, Ki, dt_last_run, 0) + controller.integral = integral + u = pwm(t, duty=setpoint_duty, period=pwm_period) + + controller.t_prev = t + return u * P + + + + # Setup physcal system elements + m_floor = ThermalMass(C=1000*4*5, T_initial=T0, name=floor_name) + m_air = ThermalMass(C=1000*8*5, T_initial=T0, name=air_name) + + R_floor_gnd = ThermalResistance(R=0.01, name="R_f_gnd") + R_air_gnd = ThermalResistance(R=0.1, name="R_a_gnd") + R_floor_air = ThermalResistance(R=0.05, name="R_f_a") + R_roof_gnd = ThermalResistance(R=0.01, name="R_r_gnd") + R_air_roof = ThermalResistance(R=1, name="R_a_r") + + heat_source = PowerSource(controller, name="Heating cable") + + # Connect elements + m_floor.connect(heat_source) + R_floor_gnd.connect(m_floor, GND) + + # m_air.connect(sun_source) + R_air_gnd.connect(m_air, GND) + R_floor_air.connect(m_floor, m_air) + + # Simulate + thermal_masses = [m_floor, m_air] + system = ThermalSystem(thermal_masses) + system.generate_diagram('img/multi_mass.png') + dt = 1 # [s] + total_time = 3600*2 # [s] + system.simulate(dt, total_time, solver=Solver.RK4) \ No newline at end of file From 3b2ef4291b4e041eb95354742b4abbaee7efe656 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Isak=20=C3=85slund?= Date: Thu, 21 Nov 2024 17:50:16 +0100 Subject: [PATCH 4/4] feat: add duration for simulation to thermal system object --- src/thermal_sim/thermal_sim.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/thermal_sim/thermal_sim.py b/src/thermal_sim/thermal_sim.py index b34ed3b..cc2091a 100644 --- a/src/thermal_sim/thermal_sim.py +++ b/src/thermal_sim/thermal_sim.py @@ -176,6 +176,7 @@ def custom_ode_solver(self, y0, t:float, solver:Solver): def simulate(self, dt:float, duration:float, plot:bool=True, solver:Solver=Solver.EULER): self.dt = dt + self.t_max = duration time = np.arange(0, duration, dt) y0 = [mass.T for mass in self.thermal_masses]