%% Q4.m —— 多目标遗传算法 + 传输矩阵 + 斯特藩–玻尔兹曼
clear; clc;
%% 1. 读取 nk 数据并清洗（按你的本地路径）
dataStruct.PDMS = cleanData(readmatrix('C:\Users\Administrator\Desktop\亚太杯代码\Zhang-10-1-0.4-20μm.csv')); % [wl, n, k]
dataStruct.PMMA = cleanData(readmatrix('C:\Users\Administrator\Desktop\亚太杯代码\Zhang-MitsubishiPMMA.csv'));
dataStruct.Ag = cleanData(readmatrix('C:\Users\Administrator\Desktop\亚太杯代码\CiesielskiAg.csv'));
dataStruct.SiO2 = cleanData(readmatrix('C:\Users\Administrator\Desktop\亚太杯代码\Herguedas-SiO1.92.csv'));
dataStruct.TiO2 = cleanData(readmatrix('C:\Users\Administrator\Desktop\亚太杯代码\KischkatTIO2.csv'));
%% 2. 读取成本、密度、加工费
costTbl = readtable('C:\Users\Administrator\Desktop\成本-密度.xlsx', ...
'Sheet',1, ...
'VariableNamingRule','preserve');
% Excel 里的真实列名（注意带下标 ₂）
rawNames = {'Ag','PMMA','SiO₂','TiO₂','PDMS'};
% 在代码中统一使用的名字（不带下标）
stdNames = {'Ag','PMMA','SiO2','TiO2','PDMS'};
for i = 1:numel(rawNames)
colName = rawNames{i}; % 表格中的列名，例如 'SiO₂'
fieldName = stdNames{i}; % matProp 中用的字段名，例如 'SiO2'
% 三行分别是：成本 / 密度 / 加工费
costStr = costTbl{1,colName}; % 第一行：成本
densStr = costTbl{2,colName}; % 第二行：密度
procStr = costTbl{3,colName}; % 第三行：加工费
% 提取数字
matProp.(fieldName).cost = extractNumber(costStr); % $/kg
densVal = extractNumber(densStr); % g/mL
matProp.(fieldName).rho = densVal * 1000; % kg/m^3
matProp.(fieldName).proc = extractNumber(procStr); % $/m^2
end
%% 3. 基本物理参数与波长网格
sigma = 5.670374e-8; % 斯特藩–玻尔兹曼常数
Ts = 320; % K，器件表面温度
Tenv = 300; % K，环境等效温度
lambda = linspace(8,13,101); % μm，大气窗口 8–13 μm
%% 4. 三层膜材料组合：PDMS 固定为中间层，Ag 只做基底
materialCombos = {
{'PMMA','PDMS','SiO2'};
{'PMMA','PDMS','TiO2'};
{'SiO2','PDMS','PMMA'};
{'SiO2','PDMS','TiO2'};
{'TiO2','PDMS','PMMA'};
{'TiO2','PDMS','SiO2'};
};
%% 5. 遗传算法参数（厚度范围在循环里按材料设置）
options = optimoptions('gamultiobj',...
'PopulationSize', 80, ...
'MaxGenerations', 80, ...
'UseParallel', false, ...
'Display','iter');
%% 6. 对每个材料组合进行多目标优化
for c = 1:numel(materialCombos)
mats = materialCombos{c}; % 比如 {'PMMA','PDMS','SiO2'}
% === 6.0 根据材料类型设置每一层的厚度上下限 ===
lb = zeros(1,3);
ub = zeros(1,3);
for j = 1:3
matName = mats{j};
switch matName
case 'PDMS'
% PDMS 厚度 50–80 μm
lb(j) = 50;
ub(j) = 80;
case 'PMMA'
% PMMA 厚度 50–100 μm
lb(j) = 50;
ub(j) = 100;
case {'SiO2','TiO2'}
% SiO2 / TiO2 厚度 0.1–0.5 μm
lb(j) = 0.1;
ub(j) = 0.5;
otherwise
% 理论上不会进入，因为 Ag 不在材料组合里
lb(j) = 0.1;
ub(j) = 1.0;
end
end
% ====================================================
% 多目标目标函数 [成本, -冷却功率]
fitnessFcn = @(d) multiObjFun(d, mats, lambda, ...
matProp, dataStruct, sigma, Ts, Tenv);
fprintf('=== 正在优化组合 %d: %s / %s / %s ===\n', ...
c, mats{1}, mats{2}, mats{3});
[xPareto, fPareto] = gamultiobj(fitnessFcn, 3, ...
[],[],[],[], lb, ub, [], options);
% 存储原始 Pareto 结果（可选）
result(c).materials = mats;
result(c).xPareto = xPareto;
result(c).fPareto = fPareto;
%% 6.1 不再限制冷却功率，直接输出所有 Pareto 解到 Excel
Ctotal = fPareto(:,1); % 成本
Pcool = -fPareto(:,2); % 冷却功率（目标里是 -Pcool）
rows = size(xPareto,1);
mat1 = repmat(string(mats{1}), rows, 1);
mat2 = repmat(string(mats{2}), rows, 1);
mat3 = repmat(string(mats{3}), rows, 1);
T = table(mat1,mat2,mat3, ...
xPareto(:,1),xPareto(:,2),xPareto(:,3), ...
Ctotal, Pcool, ...
'VariableNames',{'Layer1','Layer2','Layer3', ...
'd1_um','d2_um','d3_um', ...
'Cost_USD_m2','Pcool_W_m2'});
filename = sprintf('Combo_%d_results.xlsx', c);
writetable(T, filename);
fprintf('组合 %d 已输出到 %s（所有 Pareto 解）。\n', c, filename);
end
disp('所有组合的 Pareto 结果已输出到当前文件夹的 Excel 文件。');
%% ========= 辅助函数区域 =========
function data = cleanData(raw)
% 去掉 NaN 行并按波长升序排序
data = raw;
data = data(~any(isnan(data),2),:);
data = sortrows(data,1);
end
function num = extractNumber(str)
% 从类似 "2.1美元" 这种字符串里提取数字 2.1
if iscell(str)
str = str{1};
end
str = char(str); % 转为 char 方便正则处理
pat = '[0-9.]+';
tok = regexp(str, pat, 'match');
if isempty(tok)
num = NaN;
else
num = str2double(tok{1});
end
end
function f = multiObjFun(d, mats, lambda, matProp, dataStruct, sigma, Ts, Tenv)
% 多目标： f = [成本, -冷却功率]
% 1. 用 TMM 计算 8–13 μm 光谱发射率
eps_lambda = calcEmissivityTMM(mats, d, lambda, dataStruct);
% 2. 有效发射率（简单平均）
eps_eff = mean(eps_lambda);
% 3. 冷却功率（Stefan–Boltzmann 简化）
deltaP_bb = sigma*(Ts^4 - Tenv^4);
Pcool = eps_eff * deltaP_bb;
% 4. 成本计算（按 1 m^2 面积）
C_total = 0;
for j = 1:3
matName = mats{j}; % 如 'SiO2'
p = matProp.(matName);
rho = p.rho; % kg/m^3
c = p.cost; % $/kg
cpr = p.proc; % $/m^2
dj_m = d(j)*1e-6; % μm -> m
mass = rho * dj_m * 1; % 面积 = 1 m^2
C_total = C_total + c*mass + cpr;
end
% 5. 返回多目标
f = [C_total, -Pcool];
end
function eps_lambda = calcEmissivityTMM(mats, d, lambda, dataStruct)
% mats: {'PDMS','SiO2','TiO2'} 等
% d : [d1,d2,d3]，单位 μm
% lambda: 波长数组，单位 μm
Nlambda = numel(lambda);
eps_lambda = zeros(1,Nlambda);
for k = 1:Nlambda
lam = lambda(k);
% 1) 各层复折射率
n_layers = zeros(1,numel(mats));
for j = 1:numel(mats)
matName = mats{j};
nkData = dataStruct.(matName); % [wl, n, k]
wl = nkData(:,1);
ncol = nkData(:,2);
kcol = nkData(:,3);
n_real = myInterp1(wl, ncol, lam);
k_imag = myInterp1(wl, kcol, lam);
n_layers(j) = n_real + 1i*k_imag;
end
% 2) 基底 Ag 的复折射率（只做底层，不出现在 mats 里）
nkAg = dataStruct.Ag;
wlA = nkAg(:,1);
nAcol = nkAg(:,2);
kAcol = nkAg(:,3);
nA = myInterp1(wlA, nAcol, lam);
kA = myInterp1(wlA, kAcol, lam);
n_sub = nA + 1i*kA;
% 3) 传输矩阵
M = eye(2);
for j = 1:numel(mats)
n_j = n_layers(j);
eta_j = 1/n_j;
delta_j = 2*pi*n_j*d(j)/lam; % 注意 d 和 lam 都是 μm
Mj = [cos(delta_j), 1i*sin(delta_j)/eta_j;
1i*eta_j*sin(delta_j), cos(delta_j)];
M = M * Mj;
end
% 4) 计算等效反射系数
n0 = 1; % 空气
eta0 = 1/n0;
etas = 1/n_sub; % 基底 Ag
num = eta0*M(1,1) + eta0*etas*M(1,2) - M(2,1) - etas*M(2,2);
den = eta0*M(1,1) + eta0*etas*M(1,2) + M(2,1) + etas*M(2,2);
r = num / den;
R = abs(r)^2;
% 5) 发射率 = 1 - 反射率（无透射）
eps_lambda(k) = max(0, min(1, 1-R)); % 截断数值误差到 [0,1]
end
end
function vq = myInterp1(x, v, xq)
% 简单的一维线性插值 + 线性外推，避免使用内置 interp1
% x: N×1 单调（升序），v: N×1，xq: 标量
x = x(:);
v = v(:);
% 如果只一个点，直接返回
if numel(x) == 1
vq = v(1);
return;
end
% 保证 x 升序
[x, idx] = sort(x);
v = v(idx);
if xq <= x(1)
% 左侧外推
vq = v(1) + (v(2)-v(1)) * (xq - x(1)) / (x(2) - x(1));
elseif xq >= x(end)
% 右侧外推
vq = v(end-1) + (v(end)-v(end-1)) * (xq - x(end-1)) / (x(end) - x(end-1));
else
% 区间内线性插值
i = find(x <= xq, 1, 'last');
if i == numel(x)
i = numel(x) - 1;
end
x1 = x(i); x2 = x(i+1);
v1 = v(i); v2 = v(i+1);
t = (xq - x1)/(x2 - x1);
vq = v1 + t*(v2 - v1);
end
end
将这个代码的图表绘制全部改掉，根据计算数据或者上传数据，画出归一化光谱热图，对比多种材料（如不同配方、厚度）的光谱发射率 / 反射率分布，用颜色梯度体现数值差异，替代多折线图的视觉混乱。冷却功率类数据，用响应面 3D 图，展示冷却功率随两个关键变量（如材料厚度 d、填充率 f）的三维分布，直观呈现 “最优参数区域”。还画出小提琴图，展示多组冷却功率数据的 “分布特征”（如不同批次实验、不同环境条件下的冷却功率分布），替代传统误差棒图，更全面体现数据离散度和概率密度。雷达图，综合评价材料的 “多维度性能”（如太阳反射率、大气窗口发射率、冷却功率、机械强度、成本），替代传统的 “多列数据表格”。对于光谱发射率反射率采用堆叠面积 - 折线组合图，展示光谱参数（ε/ρ）在不同波长区间的 “贡献占比”（如大气窗口 / 太阳光谱区的积分贡献），同时保留原始光谱趋势。 将以上集中图全部根据上传的 代码和数据进行改写代码，编纂出matlab完整代码，并且绘图要求能够复合sci和nature的绘图质量
最低0.47元/天 解锁文章
扫一扫
614
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
