function [Best_pos, Best_score, NFE] = MUS_SGO(ObjectionFunc, type, Max_NFE)
    addpath('./CauchyFunction');
    
    [LB, UB, D, CLB, CUB] = RParameter_s(type);
    
    c = 0.2;               % Self-reflection
    NP = 50;               % Population size
    NFE = 0;               % Evaluation counter
    
    replacement_interval = 10;  % Elimination intervals
    initial_replacement_ratio = 0.3;  
    memory_size = 5;            % Initial elimination ratio
    forget_rate = 0.9;          % Forgetting Factor
    
    %% Population initialization
    X = zeros(NP, D);  
    C_X = zeros(NP, length(LB) + length(CLB));  

    for i = 1:NP
        
        for j = 1:D
            X(i, j) = LB(j) + rand * (UB(j) - LB(j));
        end
        
        x2 = CalculateDE_s(type, X(i,:));
        x2(3) = 1 / x2(3);
        for jj = 1:length(CLB)
            if x2(jj) < CLB(jj) || x2(jj) > CUB(jj)
                x2(jj) = CLB(jj) + rand * (CUB(jj) - CLB(jj));
            end
        end
        Iph = x2(1);
        Io = x2(2);
        Rsh = x2(3);
        C_X(i,:) = [Iph, Io, X(i,1), Rsh, X(i,2)];
    end

    
    fitness = zeros(1, NP);
    for i = 1:NP
        fitness(i) = ObjectionFunc(type, C_X(i,:));
        NFE = NFE + 1;
    end

    [Best_score, idx] = min(fitness);
    Best_pos = C_X(idx,:);
    GBestX = X(idx,:);  
    
    
    memory = struct('positions', [], 'scores', []);
    generation_counter = 0;  

    while NFE < Max_NFE
        generation_counter = generation_counter + 1;
        
        % -------- Improving phase --------
        for i = 1:NP
            Temp = c * X(i,:) + rand(1,D) .* (GBestX - X(i,:));
            
         
            for j = 1:D
                if Temp(j) > UB(j) || Temp(j) < LB(j)
                    Temp(j) = LB(j) + rand * (UB(j) - LB(j));
                end
            end

            
            x2 = CalculateDE_s(type, Temp);
            x2(3) = 1 / x2(3);
            for jj = 1:length(CLB)
                if x2(jj) < CLB(jj) || x2(jj) > CUB(jj)
                    x2(jj) = CLB(jj) + rand * (CUB(jj) - CLB(jj));
                end
            end
            new_CX = [x2(1), x2(2), Temp(1), x2(3), Temp(2)];
            fTemp = ObjectionFunc(type, new_CX);
            NFE = NFE + 1;

            
            if fTemp < fitness(i)
                X(i,:) = Temp;
                C_X(i,:) = new_CX;
                fitness(i) = fTemp;
                if fTemp < Best_score
                    Best_score = fTemp;
                    Best_pos = new_CX;
                    GBestX = Temp;
                end
            end
            
            
            if size(memory.positions, 1) >= memory_size
                memory.positions(1,:) = [];
                memory.scores(1) = [];
            end
            memory.positions = [memory.positions; Temp];
            memory.scores = [memory.scores; fTemp];
        end

        % -------- Acquiring phase --------
        for i = 1:NP
            
            k = i;
            while k == i
                k = randi(NP);
            end
            
            
            if ~isempty(memory.positions)
                [~, best_mem_idx] = min(memory.scores);
                mem_best = memory.positions(best_mem_idx,:);
                mem_weight = forget_rate ^ (memory_size - best_mem_idx);
                
                if fitness(i) < fitness(k)
                    Temp = X(i,:) + rand(1,D) .*  (X(i,:) - X(k,:)) + ...
                           mem_weight .* (mem_best - X(i,:)) + ...
                           (1 - mem_weight) .* (GBestX - X(i,:));
                else
                    Temp = X(i,:) + rand(1,D) .* (X(k,:) - X(i,:)) + ...
                           mem_weight .* (mem_best - X(i,:)) + ...
                           (1 - mem_weight) .* (GBestX - X(i,:));
                end
            else
                if fitness(i) < fitness(k)
                    Temp = X(i,:) + rand(1,D) .* (X(i,:) - X(k,:)) + rand(1,D) .* (GBestX - X(i,:));
                else
                    Temp = X(i,:) + rand(1,D) .* (X(k,:) - X(i,:)) + rand(1,D) .* (GBestX - X(i,:));
                end
            end

           
            for j = 1:D
                if Temp(j) > UB(j) || Temp(j) < LB(j)
                    Temp(j) = LB(j) + rand * (UB(j) - LB(j));
                end
            end

        
            x2 = CalculateDE_s(type, Temp);
            x2(3) = 1 / x2(3);
            for jj = 1:length(CLB)
                if x2(jj) < CLB(jj) || x2(jj) > CUB(jj)
                    x2(jj) = CLB(jj) + rand * (CUB(jj) - CLB(jj));
                end
            end
            new_CX = [x2(1), x2(2), Temp(1), x2(3), Temp(2)];
            fTemp = ObjectionFunc(type, new_CX);
            NFE = NFE + 1;

          
            if fTemp < fitness(i)
                X(i,:) = Temp;
                C_X(i,:) = new_CX;
                fitness(i) = fTemp;
                if fTemp < Best_score
                    Best_score = fTemp;
                    Best_pos = new_CX;
                    GBestX = Temp;
                end
            end
            
            
            if size(memory.positions, 1) >= memory_size
                memory.positions(1,:) = [];
                memory.scores(1) = [];
            end
            memory.positions = [memory.positions; Temp];
            memory.scores = [memory.scores; fTemp];
        end

        if generation_counter >= replacement_interval
            generation_counter = 0;
            
            progress = NFE / Max_NFE; 
            
           
            dynamic_ratio = initial_replacement_ratio * exp(-5 * progress);
            
           
            replace_count = max(1, round(dynamic_ratio * NP));
            
            
            [~, sorted_idx] = sort(fitness, 'descend');  
            worst_idx = sorted_idx(1:replace_count);
            worst_idx = worst_idx(:); 
            
            
            for i = 1:length(worst_idx)
                idx = worst_idx(i);  
                
               
                new_X = LB + rand(1, D) .* (UB - LB);
                
               
                x2 = CalculateDE_s(type, new_X);
                x2(3) = 1 / x2(3);
                for jj = 1:length(CLB)
                    if x2(jj) < CLB(jj) || x2(jj) > CUB(jj)
                        x2(jj) = CLB(jj) + rand * (CUB(jj) - CLB(jj));
                    end
                end
                new_CX = [x2(1), x2(2), new_X(1), x2(3), new_X(2)];
                
                
                new_fitness = ObjectionFunc(type, new_CX);
                NFE = NFE + 1;
                
               
                X(idx,:) = new_X;
                C_X(idx,:) = new_CX;
                fitness(idx) = new_fitness;
                
               
                if size(memory.positions, 1) >= memory_size
                    memory.positions(1,:) = [];
                    memory.scores(1) = [];
                end
                memory.positions = [memory.positions; new_X];
                memory.scores = [memory.scores; new_fitness];
            end
        end

        if NFE >= Max_NFE
            break;
        end
    end
    
    rmpath('./CauchyFunction');
end