%% 太阳能-热泵联合供暖系统优化控制 - 第一部分
% 完整参数设置和多目标优化框架
clc; clear; close all;
warning off;
% 创建必要的目录结构
if ~exist('functions', 'dir')
mkdir('functions');
end
if ~exist('results', 'dir')
mkdir('results');
end
if ~exist('results/figures', 'dir')
mkdir('results/figures');
end
addpath('functions');
rng(42); % 设置随机种子保证可重复性
%% 1. 参数设置 - 详细参数定义
fprintf('初始化系统参数...\n');
tic;
% ========== 时间参数 ==========
total_days = 7; % 总模拟天数
interval_min = 10; % 时间间隔(分钟)
interval_h = interval_min / 60; % 时间间隔(小时)
total_steps = total_days * 24 * 60 / interval_min;
time_vec = (0:total_steps-1) * interval_min / 60;
% ========== 地理位置参数 ==========
geo.lat = 40.81; % 纬度(呼和浩特)
geo.lon = 111.65; % 经度
geo.altitude = 1063; % 海拔高度(m)
geo.timezone = 8; % 时区
% ========== 建筑热力学参数 ==========
building = struct(...
'A_floor', 120, ... % 建筑面积(m^2)
'H_room', 3, ... % 房间高度(m)
'U_wall', 0.45, ... % 墙体传热系数(W/(m^2·K))
'U_window', 2.8, ... % 窗户传热系数
'A_wall', 180, ... % 墙体面积(m^2)
'A_window', 25, ... % 窗户面积(m^2)
'A_roof', 120, ... % 屋顶面积(m^2)
'U_roof', 0.5, ... % 屋顶传热系数
'Q_infiltration', 80, ...% 渗透风量热损失(W/K)
'air_density', 1.2, ... % 空气密度(kg/m³)
'air_heat_capacity', 1005); % 空气比热容(J/(kg·K))
% 建筑热容计算
building.volume = building.A_floor * building.H_room;
building.C_thermal = building.air_density * building.air_heat_capacity * building.volume;
% ========== 太阳能集热系统参数 ==========
solar = struct(...
'A_collector', 20, ... % 集热器面积(m^2)
'eta_max', 0.72, ... % 集热器最大效率
'a1', 3.5, ... % 热损失系数1(W/(m^2·K))
'a2', 0.008, ... % 热损失系数2(W/(m^2·K^2))
'min_radiation', 20, ... % 最小有效辐射(W/m²)
'tilt_angle', 40, ... % 集热器倾角(度)
'azimuth', 180, ... % 集热器方位角(度,南向)
'flow_rate', 0.02); % 集热器流量(kg/s)
% ========== 空气源热泵参数 ==========
hp = struct(...
'COP_base', 2.85, ... % 基准COP
'COP_temp_coeff', 0.048, ... % 温度对COP的影响系数
'COP_min', 1.6, ... % 最低COP
'COP_max', 3.9, ... % 最高COP
'P_min', 800, ... % 热泵最小功率(W)
'P_max', 4800, ... % 热泵最大功率(W)
'ramp_rate', 800, ... % 热泵功率变化速率(W/步长)
'efficiency', 0.92, ... % 热泵电气效率
'response_time', 5); % 热泵响应时间(分钟)
% ========== 储热水箱参数 ==========
tank = struct(...
'volume', 2.0, ... % 水箱容积(m^3)
'T_min', 40, ... % 水箱最低温度(℃)
'T_max', 60, ... % 水箱最高温度(℃)
'U_loss', 0.48, ... % 水箱散热系数(W/(m^2·K))
'surface_area', 4.8, ... % 水箱表面积(m^2)
'initial_temp', 52, ... % 初始水箱温度(℃)
'max_flow_rate', 0.04); % 最大流量(kg/s)
% 水的热物理参数
water.density = 1000; % 密度(kg/m³)
water.cp = 4186; % 比热容(J/(kg·K))
water.mass = tank.volume * water.density;
% ========== 供暖系统参数 ==========
heating = struct(...
'T_set', 18.5, ... % 设定室温(℃)
'T_min_comfort', 16.5, ... % 最低舒适室温
'T_max_comfort', 21.5, ... % 最高舒适室温
'Q_rad_min', 1500, ... % 暖气片最小散热量(W)
'Q_rad_max', 7500, ... % 暖气片最大散热量(W)
'radiator_eff', 0.93, ... % 暖气片效率
'pipe_loss_coeff', 0.05); % 管道热损失系数
% 人员在场系数(0-1) - 24小时分布
heating.occupancy = [0.05, 0.05, 0.05, 0.05, 0.05, 0.1, 0.3, ... % 0-6点
0.6, 0.8, 0.7, 0.6, 0.6, 0.7, ... % 7-12点
0.6, 0.5, 0.6, 0.8, 0.9, 1.0, ... % 13-18点
0.9, 0.8, 0.7, 0.4, 0.2]; % 19-23点
% ========== 电价参数 ==========
price = struct(...
'peak', 0.78, ... % 峰时电价(元/kWh)
'valley', 0.32, ... % 谷时电价
'flat', 0.52, ... % 平时电价
'emission_factor', 0.785); % 碳排放因子(kgCO2/kWh)
% ========== 环境参数生成 ==========
% 呼和浩特冬季典型日温度(℃)
temp_data = [-14.5, -15.2, -15.8, -16.3, -16.5, -16.2, -15.5, ...
-14.0, -12.0, -10.5, -9.0, -8.5, -8.0, -8.5, ...
-9.5, -10.5, -11.5, -12.5, -13.8, -14.5, -15.0, ...
-15.5, -15.0, -14.5];
% 扩展到一周(考虑日温差变化)
T_out_weekly = [];
for d = 1:total_days
daily_temp = temp_data + (d-4)*0.7 + 2*randn(1,24);
T_out_weekly = [T_out_weekly, daily_temp];
end
% 插值到10分钟间隔
T_out = interp1(0:23, T_out_weekly, linspace(0,23,total_steps), 'pchip');
% 太阳辐射强度(W/m^2) - 使用简化模型
G_solar = calculate_solar_radiation(geo, total_days, interval_min, solar.tilt_angle, solar.azimuth);
% ========== 峰谷时段设置 ==========
price_vec = ones(1, total_steps) * price.flat;
for i = 1:total_steps
hour = mod(time_vec(i), 24);
% 峰时: 8:00-11:00, 18:00-21:00
if (hour >= 8 && hour < 11) || (hour >= 18 && hour < 21)
price_vec(i) = price.peak;
% 谷时: 23:00-7:00
elseif (hour >= 23) || (hour < 7)
price_vec(i) = price.valley;
end
end
% ========== 控制时段划分 ==========
control.num_periods = 8;
control.period_hours = 24 / control.num_periods;
control.periods = cell(1, control.num_periods);
for p = 1:control.num_periods
start_hour = (p-1)*control.period_hours;
end_hour = p*control.period_hours;
control.periods{p} = [start_hour, end_hour];
end
fprintf('参数初始化完成,耗时 %.2f 秒\n', toc);
%% 2. 目标函数定义
fprintf('定义目标函数...\n');
% 保存目标函数到functions目录
if ~exist(fullfile('functions', 'heatingSystemObjective.m'), 'file')
fid = fopen(fullfile('functions', 'heatingSystemObjective.m'), 'w');
fprintf(fid, 'function [F, T_room, energy_dist] = heatingSystemObjective(x, G_solar, T_out, price_vec, solar, hp, tank, water, heating, building, control, time_vec, interval_min)\n');
fprintf(fid, '%% 计算供暖系统多目标\n');
fprintf(fid, 'interval_h = interval_min / 60; %% 转换为小时\n');
fprintf(fid, 'total_steps = length(time_vec);\n');
fprintf(fid, 'steps_per_hour = 60 / interval_min;\n\n');
% ... 完整的目标函数代码(由于长度限制,实际代码需完整实现)
fclose(fid);
end
%% 3. 约束函数定义
fprintf('定义约束函数...\n');
% 保存约束函数到functions目录
if ~exist(fullfile('functions', 'powerRampConstraint.m'), 'file')
fid = fopen(fullfile('functions', 'powerRampConstraint.m'), 'w');
fprintf(fid, 'function [c, ceq] = powerRampConstraint(x, control_periods, time_vec, ramp_rate, interval_min)\n');
fprintf(fid, '%% 热泵功率变化速率约束\n');
fprintf(fid, 'num_periods = length(control_periods);\n');
fprintf(fid, 'hp_power = x(num_periods+1:end);\n');
fprintf(fid, 'ceq = [];\n');
fprintf(fid, 'c = zeros(num_periods-1, 1); %% 初始化不等式约束\n\n');
% ... 完整的约束函数代码
fclose(fid);
end
%% 4. 多目标优化设置
fprintf('设置多目标优化参数...\n');
tic;
% 优化变量:每个时段的[太阳能比例, 热泵功率]
lb = [zeros(1, control.num_periods), hp.P_min * ones(1, control.num_periods)];
ub = [ones(1, control.num_periods), hp.P_max * ones(1, control.num_periods)];
% 约束函数
nonlcon = @(x) powerRampConstraint(x, control.periods, time_vec, hp.ramp_rate, interval_min);
% 设置NSGA-II算法参数
options = optimoptions('gamultiobj', ...
'PopulationSize', 250, ...
'MaxGenerations', 150, ...
'FunctionTolerance', 1e-5, ...
'CrossoverFraction', 0.8, ...
'MutationFcn', {@mutationadaptfeasible, 0.2}, ...
'SelectionFcn', {@selectiontournament, 6}, ...
'ParetoFraction', 0.4, ...
'Display', 'iter', ...
'PlotFcn', {@gaplotpareto, @gaplotdistance}, ...
'UseParallel', true);
fprintf('优化参数设置完成,耗时 %.2f 秒\n', toc);
%% 5. 运行多目标优化
fprintf('开始多目标优化...\n');
optim_vars = 2 * control.num_periods; % 优化变量数量
% 定义目标函数句柄
obj_fun = @(x) heatingSystemObjective(x, G_solar, T_out, price_vec, ...
solar, hp, tank, water, heating, building, ...
control, time_vec, interval_min);
% 运行优化
tic;
[X, F, exitflag, output] = gamultiobj(obj_fun, optim_vars, ...
[], [], [], [], lb, ub, nonlcon, options);
optim_time = toc;
fprintf('优化完成!耗时 %.2f 分钟\n', optim_time/60);
fprintf('找到 %d 个帕累托最优解\n', size(X,1));
%% 6. 结果初步分析
fprintf('分析优化结果...\n');
% 保存优化结果
save(fullfile('results', 'optim_results.mat'), 'X', 'F', 'output', ...
'time_vec', 'control', 'geo', 'solar', 'hp', 'tank', 'heating', 'building', ...
'T_out', 'G_solar', 'price_vec', 'water');
% 找出帕累托前沿解
front = F;
[front_sorted, idx] = sortrows(front, 1); % 按能耗排序
% 计算综合评分(加权和)
weights = struct(...
'energy', 0.3, ... % 能耗
'cost', 0.3, ... % 成本
'emission', 0.2, ... % 碳排放
'comfort', 0.2); % 舒适度
% 归一化目标值
norm_front = front - min(front);
norm_front = norm_front ./ max(norm_front);
% 计算综合评分
scores = norm_front * [weights.energy; weights.cost; weights.emission; weights.comfort];
[best_score, best_idx] = min(scores);
best_solution = X(best_idx, :);
fprintf('\n===== 优化结果总结 =====\n');
fprintf('最佳折衷解得分: %.4f\n', best_score);
fprintf('总能耗: %.2f kWh\n', front(best_idx,1));
fprintf('总成本: %.2f 元\n', front(best_idx,2));
fprintf('总碳排放: %.2f kgCO₂\n', front(best_idx,3));
fprintf('平均舒适度偏差: %.3f ℃\n', front(best_idx,4));
%% 7. 生成详细能源分配方案
fprintf('生成详细能源分配方案...\n');
[~, ~, energy_dist_optim] = obj_fun(best_solution);
% 保存详细结果
save(fullfile('results', 'energy_distribution.mat'), 'energy_dist_optim');
%% 8. 保存关键变量用于后续部分
save('transition_vars.mat', 'best_solution', 'total_days', ...
'interval_min', 'time_vec', 'heating', 'price_vec', ...
'control', 'T_out', 'G_solar', 'solar', 'hp', 'tank', ...
'water', 'building', 'geo');
fprintf('第一部分已完成,请运行第二部分文件\n');
%% 太阳能-热泵联合供暖系统优化控制 - 第二部分
% 结果分析与可视化
clc; close all;
%% 加载第一部分保存的数据
fprintf('加载第一部分数据...\n');
load('temp_vars_part1.mat'); % 加载临时变量
load('results/energy_distribution.mat'); % 加载详细能源分配
load('results/optim_results.mat'); % 加载优化结果
%% 7. 结果可视化
fprintf('生成可视化结果...\n');
% 创建结果目录
if ~exist('results/figures', 'dir')
mkdir('results/figures');
end
% 1. 帕累托前沿图
figure('Position', [100, 100, 800, 600], 'Name', '帕累托前沿');
plot(F(:,1), F(:,2), 'bo', 'MarkerSize', 8, 'LineWidth', 1.5);
hold on;
plot(F(best_idx,1), F(best_idx,2), 'ro', 'MarkerSize', 10, 'LineWidth', 2);
xlabel('总能耗 (kWh)');
ylabel('总成本 (元)');
title('多目标优化帕累托前沿');
grid on;
legend('帕累托解', '最佳折衷解', 'Location', 'northwest');
saveas(gcf, 'results/figures/pareto_front.png');
% 2. 温度变化曲线
figure('Position', [100, 100, 1000, 400], 'Name', '室外温度变化');
plot(time_vec, T_out, 'b-', 'LineWidth', 1.5);
xlabel('时间 (小时)');
ylabel('温度 (°C)');
title('一周室外温度变化');
xlim([0, total_days*24]);
grid on;
saveas(gcf, 'results/figures/temperature_variation.png');
% 3. 供暖季最冷日热源出力分配图
coldest_day = find(T_out == min(T_out));
cold_day_time = time_vec(coldest_day:coldest_day+steps_per_day-1);
cold_day_solar = energy_dist_optim.solar_heat(coldest_day:coldest_day+steps_per_day-1);
cold_day_hp = energy_dist_optim.hp_heat(coldest_day:coldest_day+steps_per_day-1);
cold_day_tank = energy_dist_optim.tank_heat(coldest_day:coldest_day+steps_per_day-1);
cold_day_demand = energy_dist_optim.heat_demand(coldest_day:coldest_day+steps_per_day-1);
figure('Position', [100, 100, 1000, 600], 'Name', '最冷日热源分配');
subplot(2,1,1);
plot(cold_day_time, cold_day_demand, 'k-', 'LineWidth', 1.5, 'DisplayName', '热需求');
hold on;
plot(cold_day_time, cold_day_solar, 'r-', 'LineWidth', 1.5, 'DisplayName', '太阳能供热');
plot(cold_day_time, cold_day_hp, 'b-', 'LineWidth', 1.5, 'DisplayName', '热泵供热');
plot(cold_day_time, cold_day_tank, 'g-', 'LineWidth', 1.5, 'DisplayName', '水箱供热');
xlabel('时间 (小时)');
ylabel('功率 (W)');
title('最冷日供热功率分配');
legend('Location', 'northwest');
grid on;
% 4. 水箱温度变化
subplot(2,1,2);
plot(cold_day_time, energy_dist_optim.tank_temp(coldest_day:coldest_day+steps_per_day-1), 'm-', 'LineWidth', 1.5);
hold on;
plot(cold_day_time, tank.T_min*ones(1, steps_per_day), 'r--', 'LineWidth', 1.5);
plot(cold_day_time, tank.T_max*ones(1, steps_per_day), 'r--', 'LineWidth', 1.5);
xlabel('时间 (小时)');
ylabel('水箱温度 (°C)');
title('最冷日水箱温度变化');
legend('水箱温度', '温度限制', 'Location', 'southwest');
grid on;
saveas(gcf, 'results/figures/cold_day_energy_dist.png');
% 5. 太阳能辐射与热泵COP关系
figure('Position', [100, 100, 1000, 400], 'Name', '太阳能辐射与热泵性能');
yyaxis left;
plot(time_vec, G_solar, 'r-', 'LineWidth', 1.5);
ylabel('太阳辐射 (W/m²)');
yyaxis right;
plot(time_vec, energy_dist_optim.COP, 'b-', 'LineWidth', 1.5);
ylabel('热泵COP');
xlabel('时间 (小时)');
title('太阳能辐射与热泵COP关系');
grid on;
saveas(gcf, 'results/figures/solar_radiation_cop.png');
% 6. 电价与热泵运行成本
figure('Position', [100, 100, 1000, 400], 'Name', '电价与运行成本');
yyaxis left;
plot(time_vec, price_vec, 'm-', 'LineWidth', 1.5);
ylabel('电价 (元/kWh)');
yyaxis right;
plot(time_vec, energy_dist_optim.hp_power_cost, 'c-', 'LineWidth', 1.5);
ylabel('热泵运行成本 (元)');
xlabel('时间 (小时)');
title('电价与热泵运行成本');
grid on;
saveas(gcf, 'results/figures/price_operation_cost.png');
%% 8. 室内温度舒适度分析
fprintf('分析室内温度舒适度...\n');
% 计算舒适度指标
comfort_deviation = abs(energy_dist_optim.T_room - heating.T_set);
comfort_score = 1./(1 + comfort_deviation.^2); % 舒适度评分函数
% 按时间段分析舒适度
comfort_by_period = zeros(1, control.num_periods);
for p = 1:control.num_periods
period_indices = find(time_vec >= control.periods{p}(1) & time_vec < control.periods{p}(2));
comfort_by_period(p) = mean(comfort_score(period_indices));
end
figure('Position', [100, 100, 800, 400], 'Name', '各时段舒适度分析');
bar(comfort_by_period, 'FaceColor', [0.5 0.7 0.9]);
set(gca, 'XTick', 1:control.num_periods, 'XTickLabel', arrayfun(@(x) sprintf('%d',x), 1:control.num_periods, 'UniformOutput', false));
xlabel('控制时段');
ylabel('平均舒适度评分');
title('各时段室内舒适度评分');
ylim([0.8, 1.0]);
grid on;
saveas(gcf, 'results/figures/comfort_analysis.png');
fprintf('第二部分完成,请运行第三部分\n');
%% 太阳能-热泵联合供暖系统优化控制 - 第三部分
% 与传统系统比较及结论
clc; close all;
%% 加载数据
fprintf('加载数据...\n');
load('results/energy_distribution.mat');
load('results/optim_results.mat');
load('temp_vars_part1.mat');
%% 9. 与传统系统比较
fprintf('与传统系统比较...\n');
% 传统系统1:仅使用太阳能集热器
solar_only_control = [ones(1, control.num_periods), zeros(1, control.num_periods)];
[F_solar_only, ~, energy_dist_solar] = heatingSystemObjective(solar_only_control, G_solar, T_out, price_vec, ...
solar, hp, tank, water, heating, building, control, time_vec, interval_min);
% 传统系统2:仅使用热泵(恒功率运行)
hp_only_control = [zeros(1, control.num_periods), (hp.P_min + hp.P_max)/2 * ones(1, control.num_periods)];
[F_hp_only, ~, energy_dist_hp] = heatingSystemObjective(hp_only_control, G_solar, T_out, price_vec, ...
solar, hp, tank, water, heating, building, control, time_vec, interval_min);
% 传统系统3:仅使用热泵(基于温度调节)
% 创建更智能的仅热泵控制策略
hp_smart_control = zeros(1, 2*control.num_periods);
for p = 1:control.num_periods
period_hours = control.periods{p};
period_indices = find(time_vec >= period_hours(1) & time_vec < period_hours(2));
avg_temp = mean(T_out(period_indices));
% 根据室外温度调整热泵功率
if avg_temp > -5
hp_power = hp.P_min + 0.3*(hp.P_max - hp.P_min);
elseif avg_temp > -10
hp_power = hp.P_min + 0.5*(hp.P_max - hp.P_min);
elseif avg_temp > -15
hp_power = hp.P_min + 0.7*(hp.P_max - hp.P_min);
else
hp_power = hp.P_min + 0.9*(hp.P_max - hp.P_min);
end
hp_smart_control(control.num_periods + p) = hp_power;
end
[F_hp_smart, ~, energy_dist_hp_smart] = heatingSystemObjective(hp_smart_control, G_solar, T_out, price_vec, ...
solar, hp, tank, water, heating, building, control, time_vec, interval_min);
% 优化系统(两者结合)
[F_optim, ~, energy_dist_optim] = heatingSystemObjective(best_solution, G_solar, T_out, price_vec, ...
solar, hp, tank, water, heating, building, control, time_vec, interval_min);
% 计算改进比例
improvement.energy_vs_solar = (F_solar_only(1) - F_optim(1)) / F_solar_only(1) * 100;
improvement.energy_vs_hp = (F_hp_only(1) - F_optim(1)) / F_hp_only(1) * 100;
improvement.energy_vs_hp_smart = (F_hp_smart(1) - F_optim(1)) / F_hp_smart(1) * 100;
improvement.cost_vs_solar = (F_solar_only(2) - F_optim(2)) / F_solar_only(2) * 100;
improvement.cost_vs_hp = (F_hp_only(2) - F_optim(2)) / F_hp_only(2) * 100;
improvement.cost_vs_hp_smart = (F_hp_smart(2) - F_optim(2)) / F_hp_smart(2) * 100;
improvement.emission_vs_solar = (F_solar_only(3) - F_optim(3)) / F_solar_only(3) * 100;
improvement.emission_vs_hp = (F_hp_only(3) - F_optim(3)) / F_hp_only(3) * 100;
improvement.emission_vs_hp_smart = (F_hp_smart(3) - F_optim(3)) / F_hp_smart(3) * 100;
improvement.comfort_vs_solar = (F_optim(4) - F_solar_only(4)) / abs(F_solar_only(4)) * 100;
improvement.comfort_vs_hp = (F_optim(4) - F_hp_only(4)) / abs(F_hp_only(4)) * 100;
improvement.comfort_vs_hp_smart = (F_optim(4) - F_hp_smart(4)) / abs(F_hp_smart(4)) * 100;
fprintf('\n===== 与传统系统比较 =====\n');
fprintf('仅太阳能系统能耗: %.2f kWh, 成本: %.2f 元\n', F_solar_only(1), F_solar_only(2));
fprintf('仅热泵(恒功率)系统能耗: %.2f kWh, 成本: %.2f 元\n', F_hp_only(1), F_hp_only(2));
fprintf('仅热泵(智能)系统能耗: %.2f kWh, 成本: %.2f 元\n', F_hp_smart(1), F_hp_smart(2));
fprintf('优化系统能耗: %.2f kWh (比太阳能降低 %.1f%%, 比智能热泵降低 %.1f%%)\n', ...
F_optim(1), improvement.energy_vs_solar, improvement.energy_vs_hp_smart);
fprintf('优化系统成本: %.2f 元 (比太阳能降低 %.1f%%, 比智能热泵降低 %.1f%%)\n', ...
F_optim(2), improvement.cost_vs_solar, improvement.cost_vs_hp_smart);
fprintf('优化系统碳排放: %.2f kgCO₂ (比太阳能降低 %.1f%%, 比智能热泵降低 %.1f%%)\n', ...
F_optim(3), improvement.emission_vs_solar, improvement.emission_vs_hp_smart);
fprintf('优化系统舒适度偏差: %.4f ℃ (比太阳能改善 %.1f%%, 比智能热泵改善 %.1f%%)\n', ...
F_optim(4), improvement.comfort_vs_solar, improvement.comfort_vs_hp_smart);
%% 10. 比较结果可视化
fprintf('生成比较结果可视化...\n');
% 创建对比数据结构
systems = {'仅太阳能', '仅热泵(恒功率)', '仅热泵(智能)', '优化系统'};
energy = [F_solar_only(1), F_hp_only(1), F_hp_smart(1), F_optim(1)];
cost = [F_solar_only(2), F_hp_only(2), F_hp_smart(2), F_optim(2)];
emission = [F_solar_only(3), F_hp_only(3), F_hp_smart(3), F_optim(3)];
comfort = [F_solar_only(4), F_hp_only(4), F_hp_smart(4), F_optim(4)];
% 能耗比较图
figure('Position', [100, 100, 900, 700]);
subplot(2,2,1);
bar(energy, 'FaceColor', [0.2 0.6 0.8]);
set(gca, 'XTickLabel', systems);
ylabel('总能耗 (kWh)');
title('系统能耗比较');
grid on;
% 成本比较图
subplot(2,2,2);
bar(cost, 'FaceColor', [0.9 0.6 0.1]);
set(gca, 'XTickLabel', systems);
ylabel('总成本 (元)');
title('系统运行成本比较');
grid on;
% 碳排放比较图
subplot(2,2,3);
bar(emission, 'FaceColor', [0.5 0.8 0.3]);
set(gca, 'XTickLabel', systems);
ylabel('碳排放 (kgCO₂)');
title('系统碳排放比较');
grid on;
% 舒适度比较图
subplot(2,2,4);
bar(comfort, 'FaceColor', [0.8 0.4 0.6]);
set(gca, 'XTickLabel', systems);
ylabel('平均舒适度偏差 (°C)');
title('系统舒适度比较');
grid on;
saveas(gcf, 'results/figures/system_comparison.png');
%% 11. 能源结构分析
fprintf('分析能源结构...\n');
% 计算各类能源占比
solar_energy = sum(energy_dist_optim.solar_heat) * interval_min / 60 / 1000; % kWh
hp_energy = sum(energy_dist_optim.hp_heat) * interval_min / 60 / 1000; % kWh
tank_energy = sum(energy_dist_optim.tank_heat) * interval_min / 60 / 1000; % kWh
total_energy = solar_energy + hp_energy + tank_energy;
solar_percent = solar_energy / total_energy * 100;
hp_percent = hp_energy / total_energy * 100;
tank_percent = tank_energy / total_energy * 100;
% 能源结构图
figure('Position', [100, 100, 600, 500]);
pie([solar_percent, hp_percent, tank_percent], ...
{sprintf('太阳能 %.1f%%', solar_percent), ...
sprintf('热泵 %.1f%%', hp_percent), ...
sprintf('水箱储热 %.1f%%', tank_percent)});
title('优化系统能源结构分布');
saveas(gcf, 'results/figures/energy_composition.png');
%% 12. 经济性分析
fprintf('进行经济性分析...\n');
% 计算节省费用
weekly_savings_vs_solar = F_solar_only(2) - F_optim(2);
weekly_savings_vs_hp_smart = F_hp_smart(2) - F_optim(2);
% 假设供暖季为5个月(22周)
annual_savings_vs_solar = weekly_savings_vs_solar * 22;
annual_savings_vs_hp_smart = weekly_savings_vs_hp_smart * 22;
% 系统增量成本(优化系统相比传统系统的额外投资)
incremental_investment = struct(...
'solar_collector', 15000, ... % 太阳能集热器投资(元)
'heat_pump', 10000, ... % 热泵投资(元)
'control_system', 5000, ... % 智能控制系统投资(元)
'installation', 3000); % 安装费用(元)
total_incremental = incremental_investment.solar_collector + ...
incremental_investment.heat_pump + ...
incremental_investment.control_system + ...
incremental_investment.installation;
% 投资回收期(年)
payback_period_vs_solar = total_incremental / annual_savings_vs_solar;
payback_period_vs_hp_smart = total_incremental / annual_savings_vs_hp_smart;
fprintf('\n===== 经济性分析 =====\n');
fprintf('相比太阳能系统周节省: %.2f 元, 年节省: %.2f 元\n', weekly_savings_vs_solar, annual_savings_vs_solar);
fprintf('相比智能热泵系统周节省: %.2f 元, 年节省: %.2f 元\n', weekly_savings_vs_hp_smart, annual_savings_vs_hp_smart);
fprintf('系统增量投资: %.2f 元\n', total_incremental);
fprintf('相比太阳能系统的投资回收期: %.1f 年\n', payback_period_vs_solar);
fprintf('相比智能热泵系统的投资回收期: %.1f 年\n', payback_period_vs_hp_smart);
% 投资回收期可视化
figure('Position', [100, 100, 700, 400]);
bar([payback_period_vs_solar, payback_period_vs_hp_smart], 0.6, 'FaceColor', [0.3 0.7 0.9]);
set(gca, 'XTickLabel', {'相比太阳能系统', '相比智能热泵系统'});
ylabel('投资回收期 (年)');
title('系统投资回收期比较');
grid on;
saveas(gcf, 'results/figures/payback_period.png');
%% 13. 结论与建议
fprintf('\n===== 结论与建议 =====\n');
fprintf('1. 优化控制系统在保证室内舒适度的同时,显著降低了能耗和运行成本\n');
fprintf(' - 能耗比太阳能系统降低 %.1f%%, 比智能热泵系统降低 %.1f%%\n', improvement.energy_vs_solar, improvement.energy_vs_hp_smart);
fprintf(' - 运行成本比太阳能系统降低 %.1f%%, 比智能热泵系统降低 %.1f%%\n', improvement.cost_vs_solar, improvement.cost_vs_hp_smart);
fprintf('2. 通过合理利用峰谷电价,在谷电时段增加热泵运行和水箱储热\n');
fprintf('3. 系统集成太阳能和热泵,充分发挥两种能源的优势\n');
fprintf('4. 投资回收期为 %.1f-%.1f 年,具有良好经济性\n', min(payback_period_vs_solar, payback_period_vs_hp_smart), max(payback_period_vs_solar, payback_period_vs_hp_smart));
fprintf('5. 建议在太阳能资源丰富的地区广泛推广该系统\n');
fprintf('6. 未来可考虑加入天气预报数据进一步优化控制策略\n');
%% 14. 保存最终结果并导出数据
fprintf('保存最终结果...\n');
final_results = struct(...
'optim_performance', F_optim, ...
'solar_performance', F_solar_only, ...
'hp_performance', [F_hp_only; F_hp_smart], ...
'improvement', improvement, ...
'energy_composition', [solar_percent, hp_percent, tank_percent], ...
'economic_analysis', struct(...
'weekly_savings', [weekly_savings_vs_solar, weekly_savings_vs_hp_smart], ...
'annual_savings', [annual_savings_vs_solar, annual_savings_vs_hp_smart], ...
'incremental_investment', incremental_investment, ...
'payback_period', [payback_period_vs_solar, payback_period_vs_hp_smart]), ...
'system_comparison', struct(...
'systems', {systems}, ...
'energy', energy, ...
'cost', cost, ...
'emission', emission, ...
'comfort', comfort));
save(fullfile('results', 'final_results.mat'), 'final_results');
% 导出关键数据到CSV
csv_data = table(time_vec', T_out', G_solar', energy_dist_optim.T_room', ...
energy_dist_optim.solar_heat', energy_dist_optim.hp_heat', ...
energy_dist_optim.tank_heat', energy_dist_optim.heat_demand', ...
energy_dist_optim.tank_temp', energy_dist_optim.COP', ...
'VariableNames', {'Time_h', 'Temp_out', 'Solar_rad', 'Temp_room', ...
'Solar_heat', 'HP_heat', 'Tank_heat', 'Heat_demand', 'Tank_temp', 'HP_COP'});
writetable(csv_data, fullfile('results', 'detailed_results.csv'));
fprintf('第三部分完成,如需㶲分析请运行第四部分\n');
%% 太阳能-热泵联合供暖系统优化控制 - 第四部分
% 㶲分析与能量质量评估
clc; close all;
%% 加载数据
fprintf('加载数据...\n');
load('results/final_results.mat');
load('results/detailed_results.csv'); % 加载第三部分导出的详细数据
load('temp_vars_part1.mat');
%% 15. 㶲分析理论基础
fprintf('进行㶲分析...\n');
% 㶲(Exergy) = 能量中可转化为有用功的部分
% 物理㶲计算: Ex = m * cp * [(T - T0) - T0 * ln(T/T0)]
% 参考环境温度T0 = 273.15 + (-10) = 263.15 K (呼和浩特冬季平均低温)
% 设置参数
T0 = 263.15; % 参考环境温度(K)
T0_C = T0 - 273.15; % 转换为℃
T_sun = 5800; % 太阳表面温度(K)
cp_water = 4186; % 水的比热容 J/(kg·K)
rho_water = 1000; % 水的密度 kg/m³
%% 16. 各组件㶲效率计算
fprintf('计算各组件㶲效率...\n');
% 预分配结果数组
exergy_solar_in = zeros(total_steps, 1); % 太阳能输入㶲(W)
exergy_solar_out = zeros(total_steps, 1); % 太阳能输出㶲(W)
exergy_hp_in = zeros(total_steps, 1); % 热泵输入㶲(W)
exergy_hp_out = zeros(total_steps, 1); % 热泵输出㶲(W)
exergy_tank_in = zeros(total_steps, 1); % 水箱输入㶲(W)
exergy_tank_out = zeros(total_steps, 1); % 水箱输出㶲(W)
% 计算每个时间步的㶲
for i = 1:total_steps
% 1. 太阳能集热器㶲分析
% 输入㶲 = 太阳辐射㶲 = G_solar(i) * solar.A * (1 - T0/T_sun)
exergy_solar_in(i) = G_solar(i) * solar.A * (1 - T0/T_sun);
% 水流质量流量 (kg/s)
m_dot_solar = solar.flow_rate * rho_water / 3600;
% 集热器输出热水温度 (℃)
T_solar_out = min(solar.T_out_max, energy_dist_optim.solar_out_temp(i));
% 输出㶲 = m_dot * cp * [(T_out - T_in) - T0 * ln((T_out+273.15)/(T_in+273.15))]
if T_solar_out > water.T_in
exergy_solar_out(i) = m_dot_solar * cp_water * ...
((T_solar_out - water.T_in) - T0 * log((T_solar_out+273.15)/(water.T_in+273.15)));
end
% 2. 热泵㶲分析
% 输入㶲 = 电能输入 (电能可完全转换为功)
exergy_hp_in(i) = energy_dist_optim.hp_power(i) * 1000; % kW → W
% 输出㶲 = 热量输出 * (1 - T0/(T_supply+273.15))
T_supply = hp.T_supply; % 热泵供水温度
exergy_hp_out(i) = energy_dist_optim.hp_heat(i) * (1 - T0/(T_supply+273.15));
% 3. 储热水箱㶲分析
% 水箱输入㶲 = 进入水箱的热水㶲值
T_in_tank = max(water.T_in, min(energy_dist_optim.solar_out_temp(i), tank.T_max));
if i == 1
T_tank_prev = tank.T_initial;
else
T_tank_prev = energy_dist_optim.tank_temp(i-1);
end
% 水箱输出㶲 = 流出水箱的热水㶲值
T_out_tank = energy_dist_optim.tank_out_temp(i);
% 计算水箱输入输出㶲
exergy_tank_in(i) = energy_dist_optim.tank_in_heat(i) * ...
(1 - T0/(T_in_tank+273.15));
exergy_tank_out(i) = energy_dist_optim.tank_out_heat(i) * ...
(1 - T0/(T_out_tank+273.15));
end
% 计算各组件㶲效率
eta_ex_solar = exergy_solar_out ./ (exergy_solar_in + eps); % 太阳能集热器㶲效率
eta_ex_hp = exergy_hp_out ./ (exergy_hp_in + eps); % 热泵㶲效率
eta_ex_tank = exergy_tank_out ./ (exergy_tank_in + eps); % 储热水箱㶲效率
% 系统总㶲效率
exergy_total_in = exergy_solar_in + exergy_hp_in;
exergy_total_out = exergy_solar_out + exergy_hp_out + exergy_tank_out;
eta_ex_total = exergy_total_out ./ (exergy_total_in + eps);
%% 17. 㶲分析结果可视化
fprintf('生成㶲分析图表...\n');
% 创建结果目录
if ~exist('results/figures/exergy', 'dir')
mkdir('results/figures/exergy');
end
% 1. 各组件㶲效率随时间变化
figure('Position', [100, 100, 1000, 600], 'Name', '㶲效率随时间变化');
subplot(2,1,1);
plot(time_vec, eta_ex_solar*100, 'r-', 'LineWidth', 1.5, 'DisplayName', '太阳能集热器');
hold on;
plot(time_vec, eta_ex_hp*100, 'b-', 'LineWidth', 1.5, 'DisplayName', '热泵');
plot(time_vec, eta_ex_tank*100, 'g-', 'LineWidth', 1.5, 'DisplayName', '储热水箱');
xlabel('时间 (小时)');
ylabel('㶲效率 (%)');
title('各组件㶲效率随时间变化');
legend('Location', 'northeast');
grid on;
subplot(2,1,2);
plot(time_vec, eta_ex_total*100, 'm-', 'LineWidth', 1.5);
xlabel('时间 (小时)');
ylabel('㶲效率 (%)');
title('系统总㶲效率随时间变化');
grid on;
saveas(gcf, 'results/figures/exergy/exergy_efficiency_time.png');
% 2. 平均㶲效率比较
avg_eta_ex_solar = mean(eta_ex_solar(energy_dist_optim.solar_heat > 0)) * 100;
avg_eta_ex_hp = mean(eta_ex_hp(energy_dist_optim.hp_heat > 0)) * 100;
avg_eta_ex_tank = mean(eta_ex_tank(energy_dist_optim.tank_heat > 0)) * 100;
avg_eta_ex_total = mean(eta_ex_total) * 100;
figure('Position', [100, 100, 800, 500], 'Name', '平均㶲效率比较');
components = {'太阳能集热器', '热泵', '储热水箱', '系统总效率'};
efficiencies = [avg_eta_ex_solar, avg_eta_ex_hp, avg_eta_ex_tank, avg_eta_ex_total];
bar(efficiencies, 'FaceColor', [0.4 0.6 0.8]);
text(1:length(efficiencies), efficiencies, num2str(efficiencies', '%.1f%%'),...
'HorizontalAlignment','center', 'VerticalAlignment','bottom');
set(gca, 'XTickLabel', components);
ylabel('平均㶲效率 (%)');
title('各组件及系统平均㶲效率');
ylim([0, max(efficiencies)*1.2]);
grid on;
saveas(gcf, 'results/figures/exergy/avg_exergy_efficiency.png');
% 3. 㶲损失分析
exergy_loss_solar = exergy_solar_in - exergy_solar_out;
exergy_loss_hp = exergy_hp_in - exergy_hp_out;
exergy_loss_tank = exergy_tank_in - exergy_tank_out;
exergy_loss_total = exergy_total_in - exergy_total_out;
% 累计㶲损失
cum_loss_solar = cumsum(exergy_loss_solar);
cum_loss_hp = cumsum(exergy_loss_hp);
cum_loss_tank = cumsum(exergy_loss_tank);
cum_loss_total = cumsum(exergy_loss_total);
figure('Position', [100, 100, 1000, 600], 'Name', '累计㶲损失');
plot(time_vec, cum_loss_solar/1e6, 'r-', 'LineWidth', 1.5, 'DisplayName', '太阳能损失');
hold on;
plot(time_vec, cum_loss_hp/1e6, 'b-', 'LineWidth', 1.5, 'DisplayName', '热泵损失');
plot(time_vec, cum_loss_tank/1e6, 'g-', 'LineWidth', 1.5, 'DisplayName', '水箱损失');
plot(time_vec, cum_loss_total/1e6, 'k--', 'LineWidth', 2, 'DisplayName', '总损失');
xlabel('时间 (小时)');
ylabel('累计㶲损失 (MJ)');
title('各组件累计㶲损失');
legend('Location', 'northwest');
grid on;
saveas(gcf, 'results/figures/exergy/cumulative_exergy_loss.png');
%% 18. 能量质量评估
fprintf('进行能量质量评估...\n');
% 计算能量效率和㶲效率的比值(能量质量指标)
energy_efficiency_solar = energy_dist_optim.solar_heat ./ (G_solar' * solar.A + eps);
quality_factor_solar = eta_ex_solar ./ (energy_efficiency_solar + eps);
energy_efficiency_hp = energy_dist_optim.hp_heat ./ (energy_dist_optim.hp_power*1000 + eps);
quality_factor_hp = eta_ex_hp ./ (energy_efficiency_hp + eps);
% 能量质量指标可视化
figure('Position', [100, 100, 1000, 400], 'Name', '能量质量指标');
subplot(1,2,1);
plot(time_vec, quality_factor_solar, 'r-', 'LineWidth', 1.5);
xlabel('时间 (小时)');
ylabel('质量指标');
title('太阳能集热器能量质量指标');
ylim([0, 1]);
grid on;
subplot(1,2,2);
plot(time_vec, quality_factor_hp, 'b-', 'LineWidth', 1.5);
xlabel('时间 (小时)');
ylabel('质量指标');
title('热泵能量质量指标');
ylim([0, 1]);
grid on;
saveas(gcf, 'results/figures/exergy/energy_quality_index.png');
%% 19. 结论与应用
fprintf('\n===== 㶲分析结论与建议 =====\n');
fprintf('1. 系统平均㶲效率: %.1f%%\n', avg_eta_ex_total);
fprintf('2. 主要㶲损失环节: 太阳能集热器 (平均效率 %.1f%%) 和热泵系统 (平均效率 %.1f%%)\n', ...
avg_eta_ex_solar, avg_eta_ex_hp);
fprintf('3. 改进建议:\n');
fprintf(' a. 采用选择性吸收涂层提高太阳能集热器效率\n');
fprintf(' b. 使用变频压缩机优化热泵部分负荷性能\n');
fprintf(' c. 添加热回收装置减少㶲损失\n');
fprintf(' d. 优化水箱保温减少热损失\n');
fprintf('4. 能量质量评估显示热泵系统的能量质量较高 (平均质量指标 %.2f)\n', mean(quality_factor_hp));
%% 20. 保存㶲分析结果
exergy_results = struct(...
'exergy_solar_in', exergy_solar_in, ...
'exergy_solar_out', exergy_solar_out, ...
'exergy_hp_in', exergy_hp_in, ...
'exergy_hp_out', exergy_hp_out, ...
'exergy_tank_in', exergy_tank_in, ...
'exergy_tank_out', exergy_tank_out, ...
'eta_ex_solar', eta_ex_solar, ...
'eta_ex_hp', eta_ex_hp, ...
'eta_ex_tank', eta_ex_tank, ...
'eta_ex_total', eta_ex_total, ...
'quality_factor_solar', quality_factor_solar, ...
'quality_factor_hp', quality_factor_hp, ...
'avg_eta_ex_solar', avg_eta_ex_solar, ...
'avg_eta_ex_hp', avg_eta_ex_hp, ...
'avg_eta_ex_tank', avg_eta_ex_tank, ...
'avg_eta_ex_total', avg_eta_ex_total);
save(fullfile('results', 'exergy_analysis.mat'), 'exergy_results');
% 导出㶲分析数据
exergy_data = table(time_vec', eta_ex_solar, eta_ex_hp, eta_ex_tank, eta_ex_total, ...
quality_factor_solar, quality_factor_hp, ...
'VariableNames', {'Time_h', 'ExEff_Solar', 'ExEff_HP', 'ExEff_Tank', 'ExEff_Total', ...
'Quality_Solar', 'Quality_HP'});
writetable(exergy_data, fullfile('results', 'exergy_analysis.csv'));
fprintf('第四部分完成,㶲分析结果已保存\n'); 检查上述代码并改正错误