Minimize F_total = Σ F_i(P_i)  (i=1,2,...,N)

F_i(P_i) = a_i × P_i² + b_i × P_i + c_i

% 實數編碼 - 每個基因代表一台發電機的出力
function population = initialize_population(pop_size, n_gen, P_min, P_max)
    population = zeros(pop_size, n_gen);
    for i = 1:pop_size
        % 在約束範圍內隨機生成個體
        population(i,:) = P_min + rand(1, n_gen) .* (P_max - P_min);
    end
end

function fitness = calculate_fitness(individual, a, b, c, P_load, penalty_factor)
    % 計算總成本
    total_cost = 0;
    for i = 1:length(individual)
        P_i = individual(i);
        total_cost = total_cost + a(i)*P_i^2 + b(i)*P_i + c(i);
    end
    
    % 功率平衡約束懲罰
    total_power = sum(individual);
    power_imbalance = abs(total_power - P_load);
    
    % 適應度值（最小化問題取倒數）
    fitness = 1 / (total_cost + penalty_factor * power_imbalance^2);
end

function parents = selection(population, fitness, method='roulette')
    pop_size = size(population, 1);
    
    switch method
        case 'roulette'
            % 輪盤賭選擇
            prob = fitness / sum(fitness);
            cum_prob = cumsum(prob);
            parents = zeros(size(population));
            
            for i = 1:pop_size
                r = rand();
                idx = find(cum_prob >= r, 1);
                parents(i,:) = population(idx,:);
            end
            
        case 'tournament'
            % 錦標賽選擇
            tournament_size = 2;
            parents = zeros(size(population));
            
            for i = 1:pop_size
                contestants = randperm(pop_size, tournament_size);
                [~, best_idx] = max(fitness(contestants));
                parents(i,:) = population(contestants(best_idx),:);
            end
    end
end

function offspring = crossover(parents, crossover_rate)
    pop_size = size(parents, 1);
    n_gen = size(parents, 2);
    offspring = parents;
    
    for i = 1:2:pop_size-1
        if rand() < crossover_rate
            % 算術交叉
            alpha = rand();
            offspring(i,:) = alpha * parents(i,:) + (1-alpha) * parents(i+1,:);
            offspring(i+1,:) = alpha * parents(i+1,:) + (1-alpha) * parents(i,:);
        end
    end
end

function mutated_pop = mutation(population, mutation_rate, P_min, P_max)
    pop_size = size(population, 1);
    n_gen = size(population, 2);
    mutated_pop = population;
    
    for i = 1:pop_size
        for j = 1:n_gen
            if rand() < mutation_rate
                % 高斯變異
                mutated_pop(i,j) = population(i,j) + randn() * 0.1 * (P_max(j)-P_min(j));
                % 邊界處理
                mutated_pop(i,j) = max(P_min(j), min(P_max(j), mutated_pop(i,j)));
            end
        end
    end
end

function [best_solution, best_cost, convergence] = GA_economic_dispatch()
    % 參數設置
    pop_size = 50;          % 種羣大小
    max_gen = 200;          % 最大迭代次數
    crossover_rate = 0.8;   % 交叉概率
    mutation_rate = 0.1;    % 變異概率
    penalty_factor = 1000;  % 懲罰因子
    
    % 發電機參數 (a, b, c, P_min, P_max)
    % 示例：3台發電機
    a = [0.0056, 0.0050, 0.0048];
    b = [4.5, 4.0, 3.8];
    c = [100, 120, 110];
    P_min = [50, 50, 50];
    P_max = [200, 200, 200];
    
    P_load = 400;  % 總負荷需求
    n_gen = length(a);
    
    % 初始化種羣
    population = initialize_population(pop_size, n_gen, P_min, P_max);
    
    % 記錄收斂過程
    convergence = zeros(max_gen, 1);
    
    for gen = 1:max_gen
        % 計算適應度
        fitness = zeros(pop_size, 1);
        for i = 1:pop_size
            fitness(i) = calculate_fitness(population(i,:), a, b, c, P_load, penalty_factor);
        end
        
        % 記錄最佳解
        [best_fit, best_idx] = max(fitness);
        convergence(gen) = 1/best_fit;  % 轉換為成本
        
        % 選擇
        parents = selection(population, fitness, 'tournament');
        
        % 交叉
        offspring = crossover(parents, crossover_rate);
        
        % 變異
        mutated_pop = mutation(offspring, mutation_rate, P_min, P_max);
        
        % 精英保留
        population = mutated_pop;
        population(1,:) = parents(best_idx,:);
    end
    
    % 提取最優解
    best_solution = population(best_idx,:);
    best_cost = convergence(end);
    
    % 顯示結果
    display_results(best_solution, best_cost, a, b, c);
end

function display_results(solution, cost, a, b, c)
    fprintf('最優負荷分配結果:\n');
    for i = 1:length(solution)
        P_i = solution(i);
        cost_i = a(i)*P_i^2 + b(i)*P_i + c(i);
        fprintf('發電機%d: %.2f MW, 成本: %.2f $\n', i, P_i, cost_i);
    end
    fprintf('總成本: %.2f $\n', cost);
    fprintf('總髮電量: %.2f MW\n', sum(solution));
end

% B係數法計算網損
function P_loss = calculate_losses(P, B_coeff)
    P_loss = P * B_coeff * P';
end

% 修改適應度函數考慮網損
function fitness = calculate_fitness_with_loss(individual, a, b, c, P_load, B_coeff, penalty_factor)
    total_cost = sum(a.*individual.^2 + b.*individual + c);
    total_power = sum(individual);
    P_loss = calculate_losses(individual, B_coeff);
    power_imbalance = abs(total_power - P_loss - P_load);
    
    fitness = 1 / (total_cost + penalty_factor * power_imbalance^2);
end

% 考慮閥點效應的成本函數
function cost = valve_point_cost(P, a, b, c, e, f)
    % e, f 為閥點效應係數
    cost = a*P^2 + b*P + c + abs(e * sin(f * (P_min - P)));
end