temperature_control.py

import numpy as np
import matplotlib.pyplot as plt

# --- 1. THE PLANT (HEATING SYSTEM) MODEL ---
def plant_model(current_temp, heat_input, ambient_temp, time_step, heat_capacity=1.0, loss_factor=0.05):
    """
    Simulates a simple first-order thermal system (e.g., heating a beaker of water).
    
    dT/dt = (Heat_Input / C) - (Loss_Factor * (T - T_ambient))
    """
    
    # Calculate heat change from input and losses
    heat_gain = heat_input / heat_capacity
    heat_loss = loss_factor * (current_temp - ambient_temp)
    
    # Change in temperature over the time step
    delta_temp = (heat_gain - heat_loss) * time_step
    
    # Update temperature
    new_temp = current_temp + delta_temp
    return new_temp

# --- 2. THE PID CONTROLLER ---
class PIDController:
    def __init__(self, Kp, Ki, Kd, setpoint):
        self.Kp = Kp
        self.Ki = Ki
        self.Kd = Kd
        self.setpoint = setpoint
        self._previous_error = 0
        self._integral = 0
        self._output_limit = (0, 10)  # Max heater power 0 to 10

    def calculate_output(self, current_value, dt):
        # 1. Calculate Error
        error = self.setpoint - current_value
        
        # 2. Proportional Term (P)
        P_term = self.Kp * error
        
        # 3. Integral Term (I)
        self._integral += error * dt
        I_term = self.Ki * self._integral
        
        # 4. Derivative Term (D)
        derivative = (error - self._previous_error) / dt
        D_term = self.Kd * derivative
        
        # 5. Calculate Total Output
        output = P_term + I_term + D_term
        
        # 6. Apply Saturation (Clamping the output to realistic limits)
        output = np.clip(output, self._output_limit[0], self._output_limit[1])
        
        # 7. Update for next iteration
        self._previous_error = error
        
        return output, P_term, I_term, D_term

# --- 3. SIMULATION SETUP ---
def run_simulation(Kp, Ki, Kd, setpoint):
    TIME_DURATION = 100.0  # seconds
    TIME_STEP = 0.5       # seconds (dt)
    
    # System parameters
    ambient_temp = 20.0
    initial_temp = ambient_temp
    
    # Initialize PID and tracking arrays
    pid = PIDController(Kp, Ki, Kd, setpoint)
    time_points = np.arange(0, TIME_DURATION, TIME_STEP)
    current_temp = initial_temp
    
    # Data logging
    temp_history = [initial_temp]
    input_history = [0]
    p_history, i_history, d_history = [], [], []

    # Simulation loop
    for t in time_points[1:]:
        # Calculate control output (heater power)
        output, p, i, d = pid.calculate_output(current_temp, TIME_STEP)
        
        # Update the plant model (system temp)
        current_temp = plant_model(current_temp, output, ambient_temp, TIME_STEP)
        
        # Log data
        temp_history.append(current_temp)
        input_history.append(output)
        p_history.append(p)
        i_history.append(i)
        d_history.append(d)
        
    return time_points, temp_history, input_history, p_history, i_history, d_history

# --- 4. VISUALIZATION AND WORKSHOP DEMO ---
def plot_results(time_points, temp_history, input_history, setpoint, title, Kp, Ki, Kd):
    
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
    
    # --- Top Plot: Temperature Response ---
    ax1.plot(time_points, temp_history, label='System Temperature', color='C0')
    ax1.axhline(setpoint, color='r', linestyle='--', label='Setpoint')
    ax1.set_title(f'Temperature Response | Kp={Kp}, Ki={Ki}, Kd={Kd}', fontsize=14)
    ax1.set_ylabel('Temperature (°C)')
    ax1.grid(True)
    ax1.legend()
    
    # --- Bottom Plot: Control Effort (Heater Input) ---
    ax2.plot(time_points, input_history, label='Heater Power (Control Effort)', color='C1')
    ax2.set_title('Controller Output (Heater Power)', fontsize=14)
    ax2.set_xlabel('Time (seconds)')
    ax2.set_ylabel('Power (Units)')
    ax2.grid(True)
    
    plt.tight_layout()
    plt.show()

# Set the target temperature
SETPOINT = 50.0

print(f"--- Running PID Simulation to stabilize at {SETPOINT}°C ---")


Kp_tuned, Ki_tuned, Kd_tuned = 1, 0.02, 0.1
t, temp, input, p, i, d = run_simulation(Kp=Kp_tuned, Ki=Ki_tuned, Kd=Kd_tuned, setpoint=SETPOINT)
plot_results(t, temp, input, SETPOINT, "PID Control: Quick & Stable Response", Kp_tuned, Ki_tuned, Kd_tuned)