园区综合能源系统容量优化配置全流程解析：从业务逻辑到 MATLAB 实现

figure;  

labels = {'光伏发电', '风力发电', '储能电池', '三联供系统'};  

sizes = [spareM, spareN, spareZ, triGenCap];  

colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728'];  

pie(sizes, colors);  

legend(labels, 'Location', 'southeast');  

text(0, -1.2, ['总容量：', num2str(sum(sizes)), 'kW'], 'HorizontalAlignment', 'center');  

for d = 1:3  

    dayData = daysData{d};  

    loadElec = dayData.load(:,1);  

    gasElec = gasTurbineElec{d};  

    pvOutput = calculatePV(spareM, dayData.solar);  

    windOutput = calculateWind(spareN, dayData.wind);  

    netPower = pvOutput + windOutput - (loadElec - gasElec);  

    [~, chargeDischarge, ~] = calculateStorage(spareZ, netPower);  

    figure;  

    bar(1:24, [-(loadElec - gasElec), chargeDischarge(chargeDischarge>0), pvOutput, windOutput, gasElec, -chargeDischarge(chargeDischarge<0)], 'stacked');  

    ylim([-max(loadElec) max(pvOutput + windOutput + gasElec)]);  

    title(['典型日 - ', dayData.name]);  

    ylabel('功率 (kW)');  

    xlabel('小时');  

    legend('电负荷缺额', '储能充电', '光伏发电', '风力发电', '燃气轮机电出力', '储能放电', 'Location', 'eastoutside');  

end  

for d = 1:3  

    dayData = daysData{d};  

    netPower = calculatePV(spareM, dayData.solar) + calculateWind(spareN, dayData.wind) - (dayData.load(:,1) - gasTurbineElec{d});  

    [soc, ~, ~] = calculateStorage(spareZ, netPower);  

    figure;  

    plot(1:24, soc(1:24), 'b-o', 'LineWidth', 1.5);  

    title(['储能SOC - ', dayData.name]);  

    xlabel('小时');  

    ylabel('SOC (kWh)');  

    ylim([0 spareZ]);  

    grid on;  

end  

figure;  

plot(1:iter, bestFit(1:iter));  

xlabel('迭代次数');  

ylabel('适应度值');  

title('适应度变化曲线');  

grid on;  

%% 数据读取与预处理  

clc; clear; close all;  

% 替换为实际文件路径  

dataPath = 'C:\Users\13301\OneDrive\桌面\EnergyData.xlsx';  

solarData = xlsread(dataPath, 'Solar', 'A:C');       % 光伏辐射数据（3列，过渡日、夏日、冬日）  

windData = xlsread(dataPath, 'Wind', 'A:C');         % 风速数据  

loadTrans = xlsread(dataPath, 'Load过渡日表', 'A:C'); % 过渡日负荷（电、热、冷）  

loadSummer = xlsread(dataPath, 'Load夏日表', 'A:C');  

loadWinter = xlsread(dataPath, 'Load冬日表', 'A:C');  

% 整理三类典型日数据  

daysData = {  

    struct('name', '过渡日', 'solar', solarData(:,1), 'wind', windData(:,1), 'load', loadTrans),  

    struct('name', '夏日',   'solar', solarData(:,2), 'wind', windData(:,2), 'load', loadSummer),  

    struct('name', '冬日',   'solar', solarData(:,3), 'wind', windData(:,3), 'load', loadWinter)  

};  

%% 三联供系统容量计算  

triGenCap = 0;            % 初始化三联供系统容量  

gasTurbineElec = cell(3,1); % 存储各典型日燃气轮机电出力  

for d = 1:3  

    dayData = daysData{d};  

    coldLoad = dayData.load(:,3); % 冷负荷  

    heatLoad = dayData.load(:,2); % 热负荷  

    % 计算余热锅炉出力  

    wasteHeatOutput = coldLoad/0.7 + heatLoad/0.9;  

    currentMaxWaste = max(wasteHeatOutput);  

    % 更新三联供系统容量（取最大值）  

    triGenCap = max(triGenCap, currentMaxWaste/0.8 * 4.5);  

    % 计算燃气轮机电出力（24小时数据）  

    gasTurbineElec{d} = 3.5 * (wasteHeatOutput / 0.8);  

end  

disp(['三联供系统容量：', num2str(triGenCap), ' kW']);  

%% 遗传算法优化单个典型日  

function [bestM, bestN, bestZ] = optimizeDay(solarData, windData, loadElec, gasElec)  

    %% 算法参数  

    maxIter = 50;         % 最大迭代次数  

    popSize = 50;         % 种群规模  

    crossProb = 0.6;      % 交叉概率  

    mutProb = 0.01;       % 变异概率  

    fitThresh = 0.0005;   % 适应度变化阈值  

    chromLen = 20;        % 染色体总长度（7+7+6）  

    Mbits = 7; Nbits = 7; Zbits = 6; % 各变量二进制位数  

    Mmin = 0; Mmax = 500; % 光伏容量范围（kW）  

    Nmin = 0; Nmax = 500; % 风电容量范围（kW）  

    Zmin = 0; Zmax = 200; % 储能容量范围（kWh）  

    %% 初始化种群（二进制编码）  

    population = randi([0,1], popSize, chromLen);  

    bestFit = zeros(maxIter, 1); % 记录每代最优适应度  

    for iter = 1:maxIter  

        %% 解码染色体  

        [M, N, Z] = decodeChromosome(population, Mbits, Nbits, Zbits, Mmin, Mmax, Nmin, Nmax, Zmin, Zmax);  

        %% 计算光伏/风电出力  

        pvOutput = calculatePV(M, solarData);  

        windOutput = calculateWind(N, windData);  

        %% 计算电负荷缺额与净功率  

        elecDeficit = loadElec - gasElec; % 电负荷缺额（需补充的电力）  

        netPower = pvOutput + windOutput - elecDeficit; % 净功率（正=盈余，负=缺额）  

        %% 储能充放电计算与约束检查  

        [soc, valid] = calculateStorage(Z, netPower);  

        %% 目标函数与适应度计算  

        objFunc = 0.4*M + 0.6*N + 0.1*Z;  

        penalty = 100 * (1 - valid); % 约束不满足时惩罚项为100  

        fitness = 1 ./ (objFunc + penalty + 1);  

        %% 记录最优适应度  

        [currBestFit, idx] = max(fitness);  

        bestFit(iter) = currBestFit;  

        %% 终止条件检查  

        if iter > 1 && abs(currBestFit - bestFit(iter-1)) < fitThresh  

            disp(['迭代', num2str(iter), '代，适应度变化小于阈值，提前终止']);  

            break;  

        end  

        %% 遗传操作：选择、交叉、变异  

        population = geneticOperators(population, fitness, crossProb, mutProb, chromLen, Mbits, Nbits);  

    end  

    %% 解码最优解  

    [bestM, bestN, bestZ] = decodeChromosome(population(idx,:), Mbits, Nbits, Zbits, Mmin, Mmax, Nmin, Nmax, Zmin, Zmax);  

    %% 绘制适应度变化图  

    figure;  

    plot(1:iter, bestFit(1:iter), 'b-', 'LineWidth', 1.5);  

    xlabel('迭代次数'); ylabel('适应度值');  

    title('适应度变化曲线'); grid on;  

end  

%% 染色体解码函数  

function [M, N, Z] = decodeChromosome(chrom, Mbits, Nbits, Zbits, Mmin, Mmax, Nmin, Nmax, Zmin, Zmax)  

    % 分割染色体片段  

    Mbin = chrom(:, 1:Mbits);  

    Nbin = chrom(:, Mbits+1:Mbits+Nbits);  

    Zbin = chrom(:, Mbits+Nbits+1:end);  

    % 二进制转十进制并映射到变量范围  

    M = bin2dec(num2str(Mbin')) * (Mmax-Mmin)/(2^Mbits-1) + Mmin;  

    N = bin2dec(num2str(Nbin')) * (Nmax-Nmin)/(2^Nbits-1) + Nmin;  

    Z = bin2dec(num2str(Zbin')) * (Zmax-Zmin)/(2^Zbits-1) + Zmin;  

end  

%% 光伏出力计算函数  

function pvOutput = calculatePV(capacity, solarData)  

    pvOutput = zeros(size(solarData));  

    for t = 1:length(solarData)  

        G = solarData(t);  

        if G < 60 % 低于10%额定辐射  

            pvOutput(t) = 0;  

        elseif G <= 600 % 正常范围  

            pvOutput(t) = capacity * G / 600;  

        else % 超过额定辐射  

            pvOutput(t) = capacity;  

        end  

    end  

end