%=====================================================================
% Plik: fuzzycontroller.m
% Cel:  Rozmyty regulator prędkości wentylatora w zależności od temperatury
%       + wykresy funkcji przynależności (wejścia i wyjścia), 
%       przykładowe przycięcia i agregacja z dodatkowymi pionowymi liniami
%       ilustrującymi miejsca implikacji Mamdaniego dla przykładowej T
% Data: czerwiec 2025
%=====================================================================

%% 1. PRZYGOTOWANIE ŚRODOWISKA
clear all;
clc;

%% 2. PARAMETRY FUNKCJI PRZYNALEŻNOŚCI DLA ZBIORÓW LINGWISTYCZNYCH

% 2.1. Zbiory lingwistyczne temperatury (x ∈ [0,200])
A1 = [-20,   0, 150];  % "niska"
A2 = [  0, 100, 200];  % "średnia"
A3 = [ 50, 200, 220];  % "wysoka"

% 2.2. Zbiory lingwistyczne prędkości wentylatora (y ∈ [0,3])
B1 = [ -3,   0,  1.1];  % "mała"
B2 = [0.5, 1.5, 2.5];   % "średnia"
B3 = [1.9,   3,  3.3];  % "duża"

%% 3. PRZYKŁADOWA TEMPERATURA DO ILUSTRACJI IMPLIKACJI
T_przyklad = 125;  % wartość temperatury, dla której pokażemy, jak działa implikacja
% Obliczamy stopnie przynależności dla tej temperatury
deg1 = trimf(T_przyklad, A1);
deg2 = trimf(T_przyklad, A2);
deg3 = trimf(T_przyklad, A3);

%% 4. WYKRESY FUNKCJI PRZYNALEŻNOŚCI WEJŚCIA (TEMPERATURA)
x = linspace(0, 200, 1001)';  
mu_A1 = trimf(x, A1);
mu_A2 = trimf(x, A2);
mu_A3 = trimf(x, A3);

figure('Name','Funkcje przynależności temperatury','NumberTitle','off');
plot(x, mu_A1, 'r', 'LineWidth', 1.5); hold on;
plot(x, mu_A2, 'g', 'LineWidth', 1.5);
plot(x, mu_A3, 'b', 'LineWidth', 1.5);
% Dodajemy pionową linię w miejscu T_przyklad, żeby zilustrować, gdzie następuje implikacja
xline(T_przyklad, '--k', 'LineWidth', 1.2, 'DisplayName', sprintf('T = %g', T_przyklad));
xlabel('Temperatura T');
ylabel('\mu(T)');
title('Funkcje przynależności: A1 (niska), A2 (średnia), A3 (wysoka)');
legend('A1 – niska', 'A2 – średnia', 'A3 – wysoka', sprintf('T = %g', T_przyklad), 'Location','northwest');
grid on;
hold off;

%% 5. WYKRESY FUNKCJI PRZYNALEŻNOŚCI WYJŚCIA (PRĘDKOŚĆ)
y = linspace(0, 3, 301)';  
mu_B1 = trimf(y, B1);
mu_B2 = trimf(y, B2);
mu_B3 = trimf(y, B3);

figure('Name','Funkcje przynależności prędkości','NumberTitle','off');
plot(y, mu_B1, 'r', 'LineWidth', 1.5); hold on;
plot(y, mu_B2, 'g', 'LineWidth', 1.5);
plot(y, mu_B3, 'b', 'LineWidth', 1.5);
xlabel('Prędkość wentylatora y');
ylabel('\mu(y)');
title('Funkcje przynależności: B1 (mała), B2 (średnia), B3 (duża)');
legend('B1 – mała', 'B2 – średnia', 'B3 – duża', 'Location','northeast');
grid on;
hold off;

%% 6. PRZYKŁAD PRZYCIĘĆ I AGREGACJA DLA T = T_przyklad
% Obliczamy przycięte funkcje i agregację na podstawie deg1, deg2, deg3
con1 = min(mu_B1, deg1);  % przycinamy B1 na poziomie deg1
con2 = min(mu_B2, deg2);  % przycinamy B2 na poziomie deg2
con3 = min(mu_B3, deg3);  % przycinamy B3 na poziomie deg3
mu_out = max(max(con1, con2), con3);

% Obliczamy wartość defuzyfikowaną
y_star_przyklad = oblicz_predkosc_wentylatora(T_przyklad, A1, A2, A3, B1, B2, B3);

% Wyświetlamy stopnie przynależności i y* w konsoli
fprintf('\nPrzykład (T = %g):\n', T_przyklad);
fprintf('  Stopień przynależności do A1 (niska):   %5.3f\n', deg1);
fprintf('  Stopień przynależności do A2 (średnia): %5.3f\n', deg2);
fprintf('  Stopień przynależności do A3 (wysoka):  %5.3f\n', deg3);
fprintf('  Defuzyfikowana prędkość y* = %6.4f\n\n', y_star_przyklad);

% Rysunek przycięć i agregacji
figure('Name','Przycięcia i agregacja dla T=125','NumberTitle','off');
plot(y, con1, 'r--', 'LineWidth', 1.2); hold on;
plot(y, con2, 'g--', 'LineWidth', 1.2);
plot(y, con3, 'b--', 'LineWidth', 1.2);
plot(y, mu_out, 'k-', 'LineWidth', 2);

% Dodajemy pionową linię w miejscu y_star_przyklad, aby zilustrować wartość ostrego wyjścia
xline(y_star_przyklad, '--k', 'LineWidth', 1.2, 'DisplayName', sprintf('y* = %.3f', y_star_przyklad));

xlabel('Prędkość wentylatora y');
ylabel('\mu(y)');
title(sprintf('Przycinanie i agregacja (T = %g)', T_przyklad));
legend('con1 (mała obcięta)', 'con2 (średnia obcięta)', 'con3 (duża obcięta)', 'agregacja', sprintf('y* = %.3f', y_star_przyklad), 'Location','northwest');
grid on;
hold off;

%% 7. TESTOWY BLOK OBLICZEŃ (WYŚWIETLENIE TABELI)

test_temperatury = [0, 50, 100, 125, 160, 200];

fprintf('----------------------------------------------\n');
fprintf('  T [°C]   |   y* [prędkość w [0,3]]\n');
fprintf('----------------------------------------------\n');

for i = 1:length(test_temperatury)
    T = test_temperatury(i);
    y_star = oblicz_predkosc_wentylatora(T, A1, A2, A3, B1, B2, B3);
    fprintf('   %3.0f     |         %6.4f\n', T, y_star);
end

fprintf('----------------------------------------------\n');


%% 8. DEFINICJA FUNKCJI: trimf (funkcja przynależności trójkątnej)
function mu = trimf(x, params)
    % trimf – oblicza wartość funkcji przynależności typu trójkątnego
    % Wejście:
    %   x      – skalar lub wektor wartości (np. temperatura lub y)
    %   params – wektor [a, b, c], gdzie:
    %            a = lewy punkt (μ=0), b = wierzchołek (μ=1), c = prawy punkt (μ=0)
    % Wyjście:
    %   mu     – wektor tej samej wielkości co x, z wartościami μ(x)
    
    a = params(1);
    b = params(2);
    c = params(3);
    
    mu = zeros(size(x));  % inicjalizacja zerami
    
    % Krawędź rosnąca: a < x <= b
    idx1 = (x > a) & (x <= b);
    if b ~= a
        mu(idx1) = (x(idx1) - a) ./ (b - a);
    end
    
    % Wierzchołek: x == b
    mu(x == b) = 1;
    
    % Krawędź malejąca: b < x < c
    idx2 = (x > b) & (x < c);
    if c ~= b
        mu(idx2) = (c - x(idx2)) ./ (c - b);
    end
    
    % Dla x <= a lub x >= c, mu pozostaje 0
end


%% 9. DEFINICJA FUNKCJI: oblicz_predkosc_wentylatora
function y_star = oblicz_predkosc_wentylatora(temperatura, A1, A2, A3, B1, B2, B3)
    % oblicz_predkosc_wentylatora – realizuje algorytm Mamdaniego:
    %   1) fuzyfikacja temperatury (stopnie przynależności do A1,A2,A3)
    %   2) przycinanie funkcji przynależności prędkości (MIN)
    %   3) agregacja (MAX)
    %   4) defuzyfikacja (centroid)
    %
    % Wejście:
    %   temperatura – wartość T w skali [0,200]
    %   A1, A2, A3  – parametry [a,b,c] dla zbioru „temperatury”
    %   B1, B2, B3  – parametry [a,b,c] dla zbioru „prędkości”
    % Wyjście:
    %   y_star – wyjściowa prędkość wentylatora w skali [0,3]
    
    % 9.1. FUZYFIKACJA: stopnie przynależności temperatury
    deg1 = trimf(temperatura, A1);  
    deg2 = trimf(temperatura, A2);  
    deg3 = trimf(temperatura, A3);  
    
    % 9.2. TWORZYMY SIATKĘ Y ∈ [0,3], krok ≈0.01 (301 punktów)
    y = linspace(0, 3, 301)';  
    
    % 9.3. FUNKCJE PRZYNALEŻNOŚCI PRĘDKOŚCI B1, B2, B3 NA SIATCE Y
    mu_B1 = trimf(y, B1);  
    mu_B2 = trimf(y, B2);
    mu_B3 = trimf(y, B3);
    
    % 9.4. INFERENCJA (implikacja Mamdaniego = przycinanie przez MIN)
    con1 = min(mu_B1, deg1);
    con2 = min(mu_B2, deg2);
    con3 = min(mu_B3, deg3);
    
    % 9.5. AGREGACJA (MAX ze wszystkich przycięć)
    mu_out = max(max(con1, con2), con3);
    
    % 9.6. DEFUZYFIKACJA (środek ciężkości)
    dy = y(2) - y(1);  
    numerator   = sum(y .* mu_out) * dy;
    denominator = sum(mu_out) * dy;
    
    if denominator == 0
        % Jeśli żadna reguła się nie aktywowała (wszystkie deg = 0), zwracamy 0
        y_star = 0;
    else
        y_star = numerator / denominator;
    end
end