综合能源系统多目标优化与NSGA-II算法实践

1. 为什么综合能源系统需要多目标优化？

在能源系统规模不断扩大、用能需求日益复杂的今天，传统的单目标优化方法已经难以满足实际需求。我去年参与的一个工业园区能源改造项目就深刻印证了这一点——当我们仅考虑经济性最优时，系统碳排放量比行业标准高出37%；而单纯追求低碳排放，又会导致运行成本增加近50%。这种顾此失彼的困境，正是多目标优化算法大显身手的场景。

综合能源系统（Integrated Energy System, IES）本质上是一个包含电、热、冷、气等多种能源形式的复杂耦合系统。以典型的区域能源站为例，其核心设备包括：

• 燃气轮机（同时产生电能和热能）
• 电制冷机组
• 吸收式制冷机（利用余热制冷）
• 储电/储热装置
• 光伏发电系统

这些设备之间存在复杂的能量转换关系，比如燃气轮机的余热可以驱动吸收式制冷机，光伏发电可以优先供给电制冷机组等。系统运行需要同时考虑多个相互冲突的目标：

1. 经济性目标：最小化运行成本（燃料成本、设备维护成本等）
2. 环保性目标：最小化碳排放量
3. 能效目标：最大化能源利用率
4. 可靠性目标：最大化供电/供热可靠性

这些目标之间往往存在此消彼长的关系。例如提高燃气轮机出力可以降低运行成本，但会增加碳排放；增加储能系统充放电次数可能提高能效，但会降低设备寿命。传统的加权求和法难以准确表达这种复杂关系，而NSGA-II这类多目标优化算法则能给出完整的Pareto最优解集，为决策者提供全面的方案选择。

关键提示：在实际项目中，我们常发现决策者最初认为重要的目标，在看到Pareto前沿后往往会改变优先级。这就是为什么可视化呈现解集比强行确定权重更有价值。

2. NSGA-II算法核心原理拆解

2.1 非支配排序：解决方案的层次划分

NSGA-II（Non-dominated Sorting Genetic Algorithm II）的核心创新在于其分层筛选机制。我曾用这样一个类比向客户解释：假设你正在为一支足球队选拔队员，需要同时考虑技术、速度和体能三个指标。非支配排序就像先选出三项全优的球员（第一前沿），然后排除这些人，再从剩余球员中选出次优组合（第二前沿），以此类推。

数学上，对于最小化问题，解x支配解y的定义为：

• ∀i∈[1,M]: f_i(x) ≤ f_i(y)
• ∃j∈[1,M]: f_j(x) < f_j(y)

其中M是目标函数数量。算法实现时，每个解需要计算两个关键指标：

1. 支配计数np：被多少个其他解支配
2. 支配集合Sp：支配哪些其他解

下面是快速非支配排序的伪代码实现：

function [fronts] = fastNonDominatedSort(population) fronts = {}; for p = 1:length(population) Sp = []; np = 0; for q = 1:length(population) if dominates(population(p), population(q)) Sp = [Sp q]; elseif dominates(population(q), population(p)) np = np + 1; end end if np == 0 population(p).rank = 1; fronts{1} = [fronts{1} p]; end end i = 1; while ~isempty(fronts{i}) Q = []; for p = fronts{i} for q = Sp population(q).np = population(q).np - 1; if population(q).np == 0 population(q).rank = i+1; Q = [Q q]; end end end i = i+1; fronts{i} = Q; end end

2.2 拥挤度计算：保持解集多样性

在能源优化项目中，我们经常遇到解集过度集中在某些区域的问题。NSGA-II通过拥挤度距离（crowding distance）来解决这个问题。这个概念可以理解为：在目标空间中，某个解与其相邻解之间的"私人空间"大小。

计算步骤包括：

1. 对每个前沿层内的解按各目标函数值排序
2. 边界解（最大值和最小值）赋予无限拥挤度
3. 中间解的拥挤度为相邻解在各目标维度上的距离之和

Matlab实现示例：

function population = calculateCrowdingDistance(population, front) numObjectives = size(population(1).objectives, 2); for i = 1:numObjectives [~, order] = sort([population(front).objectives](i)); population(front(order(1))).distance = Inf; population(front(order(end))).distance = Inf; for j = 2:length(front)-1 population(front(order(j))).distance = ... population(front(order(j))).distance + ... (population(front(order(j+1))).objectives(i) - ... population(front(order(j-1))).objectives(i)) / ... (max([population(front).objectives](i)) - ... min([population(front).objectives](i))); end end end

2.3 精英保留策略：代际间的智慧传承

NSGA-II相比初代NSGA的最大改进就是引入了精英保留策略。在实际编码中，我通常采用以下步骤：

1. 合并父代和子代种群（大小2N）
2. 进行非支配排序
3. 按前沿层级从高到低填充新种群
4. 同一前沿层内按拥挤度从大到小选择
5. 直到填满N个个体为止

这种策略既保留了优秀基因，又避免了早熟收敛。在某个微网优化项目中，采用精英策略后，算法收敛所需的代数减少了约40%。

3. 综合能源系统建模关键点

3.1 设备数学模型构建

3.1.1 燃气轮机模型

燃气轮机是典型的电热联产设备，其数学模型需要同时考虑电效率和热效率：

P_gt = η_elec × Q_fuel H_gt = η_heat × Q_fuel

其中η_elec通常为25%-40%，η_heat可达40%-50%。在实际项目中，我发现采用二次曲线拟合厂家提供的性能数据比固定效率更准确：

% 某型号燃气轮机的拟合模型 function [P_gt, H_gt] = gasTurbineModel(fuelInput) % 电功率输出 (MW) P_gt = 0.32*fuelInput - 0.00018*fuelInput.^2; % 热功率输出 (MW) H_gt = 0.41*fuelInput - 0.00022*fuelInput.^2; % 约束条件 P_gt = min(max(P_gt, 2), 10); % 出力范围2-10MW H_gt = min(max(H_gt, 1.5), 8); % 热输出范围1.5-8MW end

3.1.2 储能系统模型

储能设备的建模需要特别注意SOC（State of Charge）约束：

SOC(t) = SOC(t-1) + (η_ch × P_ch - P_dis/η_dis) × Δt / Capacity

在Matlab中实现时，我通常会添加防止过充/过放的保护逻辑：

function [newSOC, actualPower] = batteryModel(SOC, power, capacity, eta_ch, eta_dis, dt) if power > 0 % 充电 possible = min(power, (capacity*0.95 - SOC)*eta_ch/dt); newSOC = SOC + possible*dt/eta_ch; actualPower = possible; else % 放电 possible = max(power, (SOC - capacity*0.05)*eta_dis/dt); newSOC = SOC + possible*dt*eta_dis; actualPower = possible; end end

3.2 多目标函数设计

3.2.1 经济性目标

运行成本通常包括：

• 燃料成本：∑(燃气轮机燃料消耗×燃气价格)
• 购电成本：∑(从电网购电×电价)
• 维护成本：∑(设备出力×单位维护系数)

function cost = economicObjective(schedule, gasPrice, elecPrice) fuelCost = sum(schedule.gtFuel * gasPrice); purchaseCost = sum(max(0, schedule.load - schedule.pv) * elecPrice); maintenance = 0.02*sum(schedule.gtPower) + 0.01*sum(abs(schedule.batteryPower)); cost = fuelCost + purchaseCost + maintenance; end

3.2.2 环保性目标

碳排放主要来自：

• 燃气轮机：燃料消耗×碳排放系数
• 电网购电：购电量×电网排放因子

function emission = environmentalObjective(schedule, gasEmission, gridEmission) gtEmission = sum(schedule.gtFuel * gasEmission); gridEmission = sum(max(0, schedule.load - schedule.pv) * gridEmission); emission = gtEmission + gridEmission; end

3.3 系统约束处理

3.3.1 能量平衡约束

电功率平衡：

P_gt + P_pv + P_grid + P_battery_dis = P_load + P_battery_ch + P_electricChiller

热功率平衡：

H_gt + H_heatPump = H_heating + H_absorptionChiller

在Matlab中通常转化为不等式约束：

function [c, ceq] = powerBalanceConstraints(schedule, load, pv) % 电功率不平衡量 imbalance = schedule.gtPower + pv + max(0, schedule.gridImport) ... - load - max(0, schedule.batteryCharge) ... - schedule.electricChiller; ceq = [imbalance]; c = []; end

3.3.2 设备运行约束

以燃气轮机为例：

P_gt_min ≤ P_gt ≤ P_gt_max Ramp_down ≤ P_gt(t) - P_gt(t-1) ≤ Ramp_up

在NSGA-II中，我通常采用罚函数法处理约束：

function penalizedFitness = applyPenalties(originalFitness, violations) penaltyFactor = 1e6; % 根据问题规模调整 penalizedFitness = originalFitness + penaltyFactor * sum(violations.^2); end

4. Matlab实现全流程解析

4.1 算法参数配置

经过多个项目实践，我总结出以下参数设置经验：

params.popSize = 100; % 种群大小：复杂问题需要更大种群 params.maxGen = 200; % 最大代数：通常100-500代 params.pCrossover = 0.9; % 交叉概率：0.8-0.95 params.pMutation = 0.1; % 变异概率：1/染色体长度 params.etaC = 20; % 交叉分布指数：10-30 params.etaM = 20; % 变异分布指数：15-30 params.eliteRatio = 0.1; % 精英保留比例：0.05-0.2

调试技巧：可以先运行少量代数（如50代）快速查看解集分布，再调整参数。我曾发现当pMutation超过0.15时，解集多样性会显著提高但收敛速度下降。

4.2 染色体编码设计

对于综合能源调度问题，我推荐采用实数编码。以24小时调度为例：

染色体结构： [GT_1, GT_2, ..., GT_24, % 燃气轮机出力 Bat_1, Bat_2, ..., Bat_24, % 电池充放电功率 Grid_1, ..., Grid_24] % 电网交互功率

编码示例：

function pop = initializePopulation(popSize, nVars, lb, ub) pop = zeros(popSize, nVars); for i = 1:popSize pop(i,:) = lb + (ub-lb).*rand(1,nVars); end end

4.3 遗传算子实现

4.3.1 模拟二进制交叉（SBX）

function [child1, child2] = sbxCrossover(parent1, parent2, etaC, lb, ub) u = rand(size(parent1)); beta = zeros(size(parent1)); beta(u<=0.5) = (2*u(u<=0.5)).^(1/(etaC+1)); beta(u>0.5) = (1./(2*(1-u(u>0.5)))).^(1/(etaC+1)); child1 = 0.5*((1+beta).*parent1 + (1-beta).*parent2); child2 = 0.5*((1-beta).*parent1 + (1+beta).*parent2); % 边界处理 child1 = min(max(child1, lb), ub); child2 = min(max(child2, lb), ub); end

4.3.2 多项式变异

function mutated = polynomialMutation(individual, etaM, lb, ub) r = rand(size(individual)); delta = zeros(size(individual)); ind = r <= 0.5; delta(ind) = (2*r(ind)).^(1/(etaM+1)) - 1; ind = r > 0.5; delta(ind) = 1 - (2*(1-r(ind))).^(1/(etaM+1)); mutated = individual + delta.*(ub-lb); mutated = min(max(mutated, lb), ub); end

4.4 结果可视化技巧

4.4.1 Pareto前沿展示

function plotParetoFront(population, front) objectives = [population(front).objectives]; scatter(objectives(1,:), objectives(2,:), 'filled'); xlabel('运行成本（万元）'); ylabel('碳排放量（吨）'); title('Pareto最优前沿'); grid on; % 标注典型解 [~, minCostIdx] = min(objectives(1,:)); [~, minEmissionIdx] = min(objectives(2,:)); text(objectives(1,minCostIdx), objectives(2,minCostIdx), ' 最低成本', 'Color','red'); text(objectives(1,minEmissionIdx), objectives(2,minEmissionIdx), ' 最低排放 ', 'Color','blue'); end

4.4.2 调度方案对比

function compareSchedules(schedule1, schedule2, time) figure; subplot(3,1,1); plot(time, schedule1.gtPower, 'r', time, schedule2.gtPower, 'b--'); ylabel('燃气轮机出力(MW)'); legend('低成本方案', '低碳方案'); subplot(3,1,2); plot(time, schedule1.batteryPower, 'r', time, schedule2.batteryPower, 'b--'); ylabel('电池功率(MW)'); subplot(3,1,3); plot(time, schedule1.gridImport, 'r', time, schedule2.gridImport, 'b--'); ylabel('电网购电(MW)'); xlabel('时间(h)'); end

5. 工程实践中的挑战与解决方案

5.1 计算效率优化

在参与某园区能源管理系统开发时，原始NSGA-II对24小时调度问题的计算时间长达6小时。通过以下优化措施，最终将时间缩短到45分钟：

1. 并行化评估：利用Matlab的parfor并行计算目标函数

parfor i = 1:popSize [f1(i), f2(i)] = evaluateIndividual(pop(i,:)); end

2. 向量化计算：避免循环操作，改用矩阵运算

% 优化前 for t = 1:24 cost = cost + gasPrice * fuel(t); end % 优化后 cost = sum(gasPrice * fuel);

3. 自适应参数调整：在进化过程中动态调整变异率

if gen > params.maxGen/2 params.pMutation = params.pMutation * 0.8; % 后期减少变异 end

5.2 解集决策支持

获得Pareto前沿后，如何选择最终实施方案是实际工程中的关键问题。我总结了几种常用方法：

1. 模糊隶属度法：

% 归一化目标值 normCost = (cost - min(cost)) / (max(cost) - min(cost)); normEmission = (emission - min(emission)) / (max(emission) - min(emission)); % 计算综合满意度 satisfaction = 0.5*(1-normCost) + 0.5*(1-normEmission); % 假设权重各50% [~, bestIdx] = max(satisfaction);

2. TOPSIS法：

% 构建决策矩阵 matrix = [cost; emission]'; % 归一化 normMatrix = matrix ./ sqrt(sum(matrix.^2)); % 定义理想解和负理想解 ideal = min(normMatrix); negativeIdeal = max(normMatrix); % 计算距离 dPlus = sqrt(sum((normMatrix - ideal).^2, 2)); dMinus = sqrt(sum((normMatrix - negativeIdeal).^2, 2)); % 计算接近度 closeness = dMinus ./ (dPlus + dMinus); [~, bestIdx] = max(closeness);

5.3 实际项目经验分享

在某医院综合能源系统优化项目中，我们遇到了几个教科书上没提过的问题：

1. 设备启停约束：燃气轮机每天最多启停2次，这需要在编码中加入特殊处理：

function valid = checkStartStopConstraints(schedule) changes = diff(schedule.gtPower > 0); startStopCount = sum(changes ~= 0); valid = startStopCount <= 2; end

2. 分时电价影响：电价峰谷差达3倍，导致最优解在电价谷时段集中充电：

electricityPrice = [0.25*ones(1,7), 0.8*ones(1,8), 1.2*ones(1,5), 0.8*ones(1,4)]; % 24小时电价

3. 天气不确定性：光伏预测误差可能达30%，解决方案是采用鲁棒优化：

% 考虑光伏出力下限 effectivePV = 0.7 * forecastPV; % 按70%的预测值计算

经过这些调整后，最终方案比原系统运行成本降低18%，碳排放减少27%。决策者最终选择了成本比最优解高5%，但碳排放低15%的折中方案。