SED_code_LS_2.txt

% Initialization
%model = 'katup_zn.slx';
model = 'fullmodel_IPC_A_bintang_ps_regulated.slx';

% Simulation configurations
dt = 0.001;

% !!! Define LOG bounds for Kp, Ki, Kd 
%
Kp_min = log10(0.01); Kp_max = log10(3);
Ki_min = -6; Ki_max = log10(0.1);
Kd_min = -6; Kd_max = log10(0.1);
%}
% !!! Define step size Kg and exploration probability
Kg = 0.05; 
E = 0.8;

% Initialize LOG params w random values within bounds
Kp_log = log10(0.03);
Ki_log = log10(0.00001);
Kd_log = log10(0.0044);

% Define best params and best obj func value
Kp_best = 10^Kp_log;
Ki_best = 10^Ki_log;
Kd_best = 10^Kd_log;

% Arrays to store data
iterations = [];
current_f = [];
best_f = [];
r1_k = [];
% SED optimization
Kp = 10^Kp_log;
Ki = 10^Ki_log;
Kd = 10^Kd_log;
f_best = objectiveFunction(model, Kp, Ki, Kd);

% Optimization loop
k_max = 100; % !!!
for k = 1:k_max
    % Evaluate obj func w current params
    f_current = objectiveFunction(model, Kp, Ki, Kd)
    current_f = [current_f; f_current];
	
    % Check if the current params provide better solution
    if f_current < f_best
        Kp_best = Kp;
        Ki_best = Ki;
        Kd_best = Kd;
        f_best = f_current;
    end

    % Generate random number r1 to decide exploration or exploitation
    r1 = rand() % Between 0 and 1
    r1_k = [r1_k; r1];
    if r1 < E
        % Exploration: generate new random params
        r2_Kp = Kp_min + (Kp_max - Kp_min)*rand();
        r2_Ki = Ki_min + (Ki_max - Ki_min)*rand();
        r2_Kd = Kd_min + (Kd_max - Kd_min)*rand();
        
        % Update params w random exploration
        Kp_log = h(log10(Kp) - Kg*r2_Kp, Kp_min, Kp_max);
        Ki_log = h(log10(Ki) - Kg*r2_Ki, Ki_min, Ki_max);
        Kd_log = h(log10(Kd) - Kg*r2_Kd, Kd_min, Kd_max);
    else 
        % Exploitation using best known params
        Kp_log = log10(Kp_best);
        Ki_log = log10(Ki_best);
        Kd_log = log10(Kd_best);
    end 

    % Evaluate obj func w new params
    Kp = 10^Kp_log;
    Ki = 10^Ki_log;
    Kd = 10^Kd_log;
    f_new = objectiveFunction(model, Kp, Ki, Kd);

    % Update best solution if the new one is better
    if f_new < f_best
        Kp_best = Kp;
        Ki_best = Ki;
        Kd_best = Kd;
        f_best = f_new;
    end 
    best_f = [best_f; f_best];
    iterations = [iterations; k]; 
end 

% Output the best found params
fprintf('Optimized PID parameters:\n');
fprintf('Kp: %f\n', Kp_best);
fprintf('Ki: %f\n', Ki_best);
fprintf('Kd: %f\n', Kd_best);

% Tabled data
dataSED = table(iterations, current_f, best_f, r1_k)
%writetable(dataSED,'HORE_katup_step_SED_100.xlsx','Sheet',1)
writetable(dataSED,'HORE_IPC_step_SED_100.xlsx','Sheet',6)