Zadanie

Poniżej znajdziecie kod prostego symulatora - umożliwia on Wam obliczenie odpowiedzi układu regulacji na wymuszenie skokowe, i obliczenie podstawowych wskaźników jakości regulacji. Wykorzystajcie ten symulator do optymalizacji nastaw regulatora wybraną techniką. 

import matplotlib.pyplot as plt
import numpy as np

class PIDSecondOrderSystem:
    """
    Symulator obiektu inercyjnego drugiego rzędu z regulatorem PID.
    """

    def __init__(self, K_p, T_1, T_2, T_n, T_d, K_r, T_i, T_d_reg, dt=0.01):
        """
        Inicjalizacja symulatora.

        Args:
            K_p (float): Wzmocnienie obiektu.
            T_1 (float): Stała czasowa pierwszego członu inercyjnego.
            T_2 (float): Stała czasowa drugiego członu inercyjnego.
            T_n (float): Stała czasowa opóźnienia transportowego (przybliżona dyskretyzacja).
            T_d (float): Stała czasowa inercji pomiaru (przybliżona dyskretyzacja).
            K_r (float): Wzmocnienie regulatora PID.
            T_i (float): Czas całkowania regulatora PID.
            T_d_reg (float): Czas różniczkowania regulatora PID.
            dt (float): Krok czasowy symulacji.
        """
        self.K_p = K_p
        self.T_1 = T_1
        self.T_2 = T_2
        self.T_n = T_n
        self.T_d = T_d
        self.K_r = K_r
        self.T_i = T_i if T_i != 0 else 1e-9 # Zapobieganie dzieleniu przez zero
        self.T_d_reg = T_d_reg
        self.dt = dt

        # Stan obiektu
        self.x1 = 0.0
        self.x2 = 0.0
        self.u_prev = 0.0 # Poprzednia wartość sterowania

        # Stan regulatora PID
        self.integral_error = 0.0
        self.previous_error = 0.0
        self.previous_measurement = 0.0

        # Stan opóźnienia transportowego i inercji pomiaru
        self.transport_delay_buffer = []
        self.measurement_inertia_buffer = []

    def simulate_step(self, setpoint):
        """
        Wykonuje jeden krok symulacji.

        Args:
            setpoint (float): Wartość zadana.

        Returns:
            float: Aktualna wartość wyjściowa obiektu.
        """
        # Obliczenie uchybu regulacji
        measurement = self.x2 # Przyjmujemy x2 jako pomiar
        error = setpoint - measurement

        # Obliczenie składowych regulatora PID
        proportional = self.K_r * error
        self.integral_error += error * self.dt
        integral = self.K_r / self.T_i * self.integral_error
        derivative = self.K_r * self.T_d_reg * (error - self.previous_error) / self.dt

        # Ograniczenie składowej całkującej (anti-windup - prosty)
        if abs(self.u_prev) > 10: # Przykład prostego ograniczenia
             if np.sign(self.u_prev) == np.sign(integral):
                 integral = self.u_prev - proportional - derivative

        # Obliczenie sygnału sterującego
        u = proportional + integral + derivative

        # Ograniczenie sygnału sterującego (np. do zakresu +/- 10)
        u = np.clip(u, -10, 10)

        self.u_prev = u
        self.previous_error = error

        # Symulacja obiektu inercyjnego drugiego rzędu (metoda Eulera w przód)
        dx1dt = (self.K_p * u - self.x1) / self.T_1
        dx2dt = (self.x1 - self.x2) / self.T_2

        self.x1 += dx1dt * self.dt
        self.x2 += dx2dt * self.dt

        # Symulacja opóźnienia transportowego (bufor)
        self.transport_delay_buffer.append(self.x2)
        if len(self.transport_delay_buffer) > int(self.T_n / self.dt):
            output_with_delay = self.transport_delay_buffer.pop(0)
        else:
            output_with_delay = self.transport_delay_buffer[0] if self.transport_delay_buffer else 0

        # Symulacja inercji pomiaru (przybliżenie 1. rzędu)
        if not self.measurement_inertia_buffer:
            self.measurement_inertia_buffer.append(output_with_delay)
        else:
            dy_dt = (output_with_delay - self.measurement_inertia_buffer[-1]) / self.T_d
            self.measurement_inertia_buffer.append(self.measurement_inertia_buffer[-1] + dy_dt * self.dt)

        return self.measurement_inertia_buffer[-1]


    def simulate_step_response(self, setpoint, simulation_time):
        """
        Symuluje odpowiedź na skokową zmianę wartości zadanej.

        Args:
            setpoint (float): Wartość zadana.
            simulation_time (float): Czas trwania symulacji.

        Returns:
            tuple: Krotka zawierająca listy czasu, wartości zadanej i wartości wyjściowej.
        """
        time_points = np.arange(0, simulation_time, self.dt)
        setpoint_values = [setpoint] * len(time_points)
        output_values = []

        # Reset stanu symulatora przed nową symulacją
        self.__init__(self.K_p, self.T_1, self.T_2, self.T_n, self.T_d,
                      self.K_r, self.T_i, self.T_d_reg, self.dt)

        for t in time_points:
            output = self.simulate_step(setpoint)
            output_values.append(output)

        return time_points, setpoint_values, output_values

    def calculate_quality_indicators(self, time_points, setpoint_values, output_values):
        """
        Oblicza podstawowe wskaźniki jakości regulacji.

        Args:
            time_points (list): Lista punktów czasowych.
            setpoint_values (list): Lista wartości zadanych.
            output_values (list): Lista wartości wyjściowych.

        Returns:
            dict: Słownik zawierający wskaźniki jakości regulacji.
        """
        errors = [setpoint_values[i] - output_values[i] for i in range(len(time_points))]

        # Czas regulacji (przyjmujemy tolerancję 2% wartości zadanej)
        settling_time = None
        tolerance = 0.02 * setpoint_values[-1] # 2% wartości zadanej na końcu symulacji
        stable_region_start = None

        # Szukamy momentu, od którego sygnał pozostaje w tolerancji
        for i in range(len(time_points) - 1, -1, -1):
            if abs(errors[i]) > tolerance:
                stable_region_start = i + 1 # Czas regulacji jest po ostatnim przekroczeniu tolerancji
                break

        if stable_region_start is not None and stable_region_start < len(time_points):
             settling_time = time_points[stable_region_start]

        # Odchyłka statyczna (różnica między wartością zadaną a wyjściową w stanie ustalonym)
        # Przyjmujemy ostatnią wartość jako stan ustalony
        static_error = setpoint_values[-1] - output_values[-1]

        # Całka z kwadratu odchyłki (ISE - Integral Squared Error)
        ise = sum([e**2 * self.dt for e in errors])

        # Całka z wartości bezwzględnej odchyłki (IAE - Integral Absolute Error)
        iae = sum([abs(e) * self.dt for e in errors])

        return {
            "Czas regulacji (2% tolerancji)": f"{settling_time:.2f} s" if settling_time is not None else "Nie osiągnięto stanu ustalonego",
            "Odchyłka statyczna": f"{static_error:.4f}",
            "ISE (Integral Squared Error)": f"{ise:.4f}",
            "IAE (Integral Absolute Error)": f"{iae:.4f}"
        }

# Przykładowe użycie klasy

# Parametry obiektu inercyjnego drugiego rzędu
Kp_obj = 2.0  # Wzmocnienie obiektu
T1_obj = 5.0  # Stała czasowa 1
T2_obj = 2.0  # Stała czasowa 2
Tn_delay = 0.5 # Opóźnienie transportowe
Td_inertia = 0.1 # Inercja pomiaru

# Parametry regulatora PID (przykładowe wartości)
Kr_reg = 1.0  # Wzmocnienie regulatora
Ti_reg = 15  # Czas całkowania
Td_reg_val = 0.5 # Czas różniczkowania

# Krok i czas symulacji
simulation_step = 0.01
simulation_duration = 30.0 # Czas trwania symulacji w sekundach
setpoint_value = 5.0 # Wartość zadana (skok)

# Utworzenie instancji symulatora
simulator = PIDSecondOrderSystem(
    K_p=Kp_obj,
    T_1=T1_obj,
    T_2=T2_obj,
    T_n=Tn_delay,
    T_d=Td_inertia,
    K_r=Kr_reg,
    T_i=Ti_reg,
    T_d_reg=Td_reg_val,
    dt=simulation_step
)

# Wykonanie symulacji odpowiedzi na skok
time_points, setpoint_values, output_values = simulator.simulate_step_response(
    setpoint=setpoint_value,
    simulation_time=simulation_duration
)

# Obliczenie wskaźników jakości regulacji
quality_indicators = simulator.calculate_quality_indicators(
    time_points,
    setpoint_values,
    output_values
)

# Wyświetlenie wskaźników jakości regulacji
print("Wskaźniki jakości regulacji:")
for key, value in quality_indicators.items():
    print(f"- {key}: {value}")

# Wyświetlenie wykresu odpowiedzi układu
plt.figure(figsize=(12, 6))
plt.plot(time_points, setpoint_values, label='Wartość zadana', linestyle='--')
plt.plot(time_points, output_values, label='Wartość wyjściowa')
plt.title('Odpowiedź układu regulacji na skok wartości zadanej')
plt.xlabel('Czas [s]')
plt.ylabel('Wartość')
plt.grid(True)
plt.legend()
plt.show()