1. Kiến trúc dự án

microgrid_pso/
 ├─ execute.m                 % Kịch bản khởi chạy chính
 ├─ adaptive_pso.m            % Lõi thuật toán PSO cải tiến
 ├─ grid_system.m             % Mô hình tham số vi lưới 24h
 ├─ objective_eval.m          % Hàm tính toán chi phí & ràng buộc
 ├─ render_plots.m            % Trực quan hóa kết quả điều phối
 └─ interface/
     └─ control_panel.mlapp   % Giao diện tham số hóa

2. Mô hình hệ thống vi lưới (Chu kỳ 24 giờ)

File grid_system.m định nghĩa các tham số kỹ thuật và dữ liệu dự báo tải, năng lượng mặt trời, gió cùng cấu hình thiết bị lưu trữ và máy phát dự phòng.

function sys = define_grid_params()
    sys.horizon = 24;
    sys.load_kw = [120,110,105,100,95,100,110,130,150,160,170,180,190,185,180,170,160,150,140,130,125,120,115,118];
    sys.pv_kw = [0,0,0,0,5,30,70,120,180,200,190,170,140,100,60,30,10,0,0,0,0,0,0,0];
    sys.wt_kw = [40,35,30,25,20,25,35,50,60,70,65,55,45,40,35,30,25,20,15,20,25,30,35,45];
    
    sys.gen_max = 200; sys.ramp_rate = 50;
    sys.batt_cap = 100; sys.batt_pwr = 50; sys.batt_eff = 0.95;
    sys.elec_tariff = [0.4*ones(1,7), 0.8*ones(1,10), 0.4*ones(1,7)];
end

3. Thuật toán PSO với trọng số thích nghi & nhiễu loạn hỗn loạn

File adaptive_pso.m triển khai cơ chế cập nhật vận tốc kết hợp hệ số quán tính suy giảm tuyến tính và ánh xạ Logistic để tránh mắc kẹt tại cực trị địa phương. Không gian tìm kiếm bao gồm 48 chiều (công suất máy phát + công suất sạc/xả pin).

function [best_pos, best_val] = run_adaptive_pso(sys, cfg)
    dim = 48; pop_size = cfg.pop_size; max_gen = cfg.max_iter;
    lb = [zeros(1,24), -sys.batt_pwr*ones(1,24)];
    ub = [sys.gen_max*ones(1,24), sys.batt_pwr*ones(1,24)];
    
    pos = lb + (ub - lb) .* rand(pop_size, dim);
    vel = randn(pop_size, dim) * 0.5;
    p_best_pos = pos; p_best_val = inf(1, pop_size);
    
    [g_best_val, idx] = min(arrayfun(@(i) objective_eval(pos(i,:), sys), 1:pop_size));
    g_best_pos = pos(idx, :);
    
    for gen = 1:max_gen
        w = cfg.w_max - (cfg.w_max - cfg.w_min) * (gen / max_gen);
        for i = 1:pop_size
            r1 = rand(1, dim); r2 = rand(1, dim);
            chaos_map = 4 * r1 .* (1 - r1);
            
            vel(i,:) = w * vel(i,:) + cfg.c1 * chaos_map .* (p_best_pos(i,:) - pos(i,:)) ...
                                  + cfg.c2 * r2 .* (g_best_pos - pos(i,:));
            pos(i,:) = pos(i,:) + vel(i,:);
            pos(i,:) = max(min(pos(i,:), ub), lb);
            
            f_curr = objective_eval(pos(i,:), sys);
            if f_curr < p_best_val(i)
                p_best_val(i) = f_curr; p_best_pos(i,:) = pos(i,:);
                if f_curr < g_best_val
                    g_best_val = f_curr; g_best_pos = pos(i,:);
                end
            end
        end
    end
    best_pos = g_best_pos; best_val = g_best_val;
end

4. Hàm mục tiêu đa thành phần

File objective_eval.m tổng hợp chi phí nhiên liệu, tiền mua điện lưới, cùng các hệ số phạt vi phạm cân bằng công suất, giới hạn SOC và tốc độ thay đổi công suất máy phát.

function obj = objective_eval(vars, sys)
    p_gen = vars(1:24); p_stor = vars(25:48);
    
    deficit = sys.load_kw - sys.pv_kw - sys.wt_kw - p_gen - p_stor;
    cost_balance = sum(abs(deficit));
    
    cost_fuel = sum(0.3 * p_gen.^2 + 0.5 * p_gen);
    grid_import = max(deficit, 0);
    cost_grid = sum(sys.elec_tariff .* grid_import);
    
    soc_trace = cumsum([0, p_stor * sys.batt_eff]); soc_trace = soc_trace(2:end);
    pen_soc = 1e3 * sum(soc_trace < 0 | soc_trace > sys.batt_cap);
    
    pen_ramp = 1e3 * sum(abs(diff(p_gen)) > sys.ramp_rate);
    
    obj = cost_balance + cost_fuel + cost_grid + pen_soc + pen_ramp;
end

5. Kịch bản thực thi & Trực quan hóa

Tệp execute.m thiết lập tham số tối ưu, gọi hàm PSO, giải mã kết quả và xuất biểu đồ vận hành.

clear; clc; close all;
sys = define_grid_params();
cfg.pop_size = 40; cfg.max_iter = 200;
cfg.w_max = 0.9; cfg.w_min = 0.4; cfg.c1 = 1.5; cfg.c2 = 1.5;

tic
[opt_vars, min_cost] = run_adaptive_pso(sys, cfg);
toc

p_gen_opt = opt_vars(1:24); p_bat_opt = opt_vars(25:48);
render_plots(sys, p_gen_opt, p_bat_opt, min_cost);

Hàm render_plots.m vẽ đồ thị phân bổ công suất theo giờ:

function render_plots(sys, p_gen, p_bat, cost)
    t = 1:sys.horizon;
    figure; 
    plot(t, sys.load_kw, 'k-', t, p_gen, 'b--', t, p_bat, 'r:', t, sys.pv_kw+sys.wt_kw, 'g-.', 'LineWidth', 1.5);
    legend('Tải', 'Máy phát Diesel', 'Hệ thống Lưu trữ', 'NL Tái tạo');
    grid on; title(sprintf('Kế hoạch điều phối - Tổng chi phí: %.2f đơn vị', cost));
    xlabel('Thời gian (giờ)'); ylabel('Công suất (kW)');
end

6. Tích hợp giao diện & Mở rộng MOPSO

Giao diện control_panel.mlapp cho phép điều chỉnh trực tiếp kích thước quần thể, số vòng lặp và hệ số học, sau đó kích hoạt callback để cập nhật đồ thị thời gian thực:

function StartCallback(app, event)
    cfg.pop_size = app.PopEditField.Value;
    cfg.max_iter = app.IterEditField.Value;
    [opt_sol, val] = run_adaptive_pso(app.sys, cfg);
    refresh_canvas(app, opt_sol, val);
end

Để chuyển sang tối ưu đa mục tiêu (MOPSO), cấu trúc hàm thích nghi được tách thành vector [cost_fuel, cost_grid, pen_soc], bổ sung kho lưu trữ ngoài (external archive) và cơ chế xếp hạng Pareto. Các module giao diện đã được chuẩn bị sẵn cho việc tích hợp phương pháp NSGA-II hoặc MOPSO tiêu chuẩn.