2021年亚太杯APMCM数学建模大赛
B题 热光电发电技术中热发射器的优化设计
为方便各位阅览及了解掌握亚太杯的写作技巧，这里非技术使用中文，公式部分由于翻译过程繁琐使用英文来撰写此文章，由不知名同学传稿
原题再现：
近年来，世界各大国一个接一个地转向“星海”，制定了各种太空探索计划。2020年，中国的“天1号”发射并通过太空进入火星；2021年，“朱融”号完成了计划任务，留在火星上，仍在寻找更多广阔宇宙的发现。为了确保月球车携带的各种仪器和设备能够在不晒太阳的情况下运行良好，并为其长期工作提供必要的技术支持，科学家们探索并开发了热光电技术。下图显示了一个热光电装置的测试原型。
热光电技术是一种利用各种热源加热热发射器（吸收器），然后通过光伏电池将热发射器的红外辐射转化为电能的技术。热源有很多种类型，包括化学能、太阳能、核能等。系统中的热发射器主要利用不同的材料结构来调节所吸收热的发射，使发射的光子大多在光伏电池的带隙波长以下。光伏电池主要转换在特定带隙波长以下的高能光子。它具有一定的带隙能量，因此也具有相应的带隙波长。例如，一个带隙波长为1100纳米的硅太阳能电池只能吸收上述波长以下的高能光子，并将其转换为电能，而电池吸收波长以上的低能光子不能通过光电效应转化为电能。相反，它们只能转化为热能，从而降低了电池的光电转换效率。因此，为了提高热光电系统的热电转换效率，必须调节热发射器的发射光谱。发射光谱的计算方法主要包括传输矩阵法(TMM)[1-2]、有限差分时域法(FDTD)和严格耦合波分析方法(RCWA)。而影响热发射器发射光谱的主要因素是材料的光学性质（折射率或介电常数）和结构性质（厚度）。王等人开发了一种亚微量厚的多层选择性太阳能吸收器，它由钨、二氧化硅和氮化硅组成，在太阳波段的吸收率可达0.95。2014年，麻省理工学院的王伊芙琳团队设计了一种光子控制的太阳能热光电装置，在实验中运行良好。在他们的工作中，热发射器是由硅和二氧化硅组成的多层薄膜结构。优化了每一层的厚度，使其发射光谱对应于铟砷化镓(InGaAssb)电池的带隙。不同材料的折射率或介电常数可以通过检索文档或参考材料[5]的光学性质数据库来找到，该数据库提供了普通材料的折射率。请按上述介绍的方式来解决以下问题。
(1) 请说明单层结构的发射光谱与材料性能（折射率、厚度）的关系，并计算出0.3-5微米范围内50纳米厚钨（如图所示）的发射光谱。
(2) 请说明多层结构的发射光谱与材料性能（折射率、厚度）的关系，并计算钨（50nm）和二氧化硅（50nm）（如下图所示）在0.3-5微米范围内形成的复合结构的发射光谱。
(3) 为了提高散热器的光谱控制能力，有时热发射器设计为窄带形式发射，即发射集中在很小的带内，从而提高热光电器件的热电转换效率，如樱井等设计的多层窄带发射器。[6]与硅、二氧化硅和锗结合。请选择合理的材料，设计一个多人热发射器，使其发射尽可能窄、高，并给出多层结构的设计参数（包括层数、每层材料和厚度）及其发射光谱。需要注意的是，这个问题中的热发射器具有1.5微米的锐而高的热发射，计算出的波长范围为0.3-5微米。
(4) 锑化镓(GaSb)电池目前还比较先进。假设其带隙波长为1.71微米。其理想化的热发射器的发射光谱大致用下图中的红色虚线表示。蓝虚线表示外部量子效率(EQE)，可以适当地考虑其影响。请选择合理的材料，为GaSb电池设计一个多层热发射器，以达到尽可能高的热电转换效率，并给出多层结构的设计参数（包括层数、每层材料和厚度）及其发射光谱。
注：这里的发射光谱都为垂直发射。
专业术语说明：
热电光电：包括热源、散热器、光伏电池和散热系统等的装置。用于将热能转化为电能。
热发射器：在加热时发出热辐射的部件。
整体求解过程概述(摘要)
在问题1中，我们可以知道波长的变化总是在增加，但折射率的变化却在增加和减少。我们首先分别进行分析波长、折射率和消光系数之间的关系，然后对波长、折射率、波长和消光系数进行多元回归建模。根据这种情况，我们建立了多元回归模型，得到了波长、折射率和消光系数之间的隐函数关系。通过代入数据，我们得到50纳米厚的钨在0.3-5 μ米中的折射率为144.354。
在问题2中，我们模拟了多层结构的发射光谱与材料特性之间的关系。基于基于转移矩阵法（TMM）的复合结构发射光谱模型，我们替换了二氧化硅和钨的数据，我们可以清楚地看到钨（50 nm）和二氧化硅（50 nm）形成的复合结构在0.3-5微米范围内的发射光谱。
在问题3中，针对该问题，我们训练了一个正确的BP神经网络，以层数、每层材料和厚度、发射光谱污染指数为输入层，作为训练样本，综合热效率作为输出层，得到多层结构综合热效率与设计参数的网络映射关系。然后，通过粒子群优化算法和遗传算法，得到6层，采用二氧化硅硅锗硅锗结构，折射率为21.7832，厚度为130.819。目前，相应的热电转换效率最高。
在问题4中，我们建立了一个基于傅里叶变换的曲线方程模型。通过建立方程并结合神经网络模型，计算出在最高热电转换效率下GaSb电池设计多层热发射器的层数、材料和每层厚度以及发射光谱为4层，4层为GaSb、GA、GA和Sb。 厚度为94.732。
问题分析：
问题分析 1
我们可以知道，波长的变化总是在增加，但折射率的变化却在增加和减少。首先分别分析波长、折射率和消光系数之间的关系，然后对波长、折射率、波长和消光系数进行多元回归分析和建模。根据这种情况，我们建立了多元回归模型，得到波长、折射率和消光系数之间的隐式函数关系。通过代入数据，我们得到50纳米厚的钨在0.3-5 μ米中的折射率为144.354。
问题分析 2
标题要求我们表达多层结构的发射光谱与材料性能（折射率和厚度）之间的关系，并计算钨（50nm）和二氧化硅（50nm）（如图所示）在0.3-5微米内形成的复合结构的发射光谱。我们模拟了多层结构的发射光谱与材料特性之间的关系。基于基于转移矩阵法（TMM）的复合结构发射光谱模型，我们替换了二氧化硅和钨的数据，我们可以清楚地看到钨（50 nm）和二氧化硅（50 nm）形成的复合结构在0.3-5微米内的发射光谱。
问题分析 3
该主题要求我们设计一个多层热发射器，使其发射尽可能窄和高，并给出多层结构的设计参数（包括层数，每层的材料和厚度）以及其发射光谱。以层数、各层材料厚度和发射光谱污染指数为输入层，作为训练样本，综合热效率作为输出层，训练出正确的BP神经网络，得到多层结构综合热效率与设计参数的网络映射关系。然后，通过粒子群优化算法和遗传算法，得到6层，采用二氧化硅硅锗硅锗结构，折射率为21.7832，厚度为130.819。目前，相应的热电转换效率最高.
问题分析 4
该主题要求我们选择合理的材料，为GaSb电池设计多层热发射器，以实现尽可能高的热电转换效率，并给出多层结构的设计参数（包括层数，材料和每层厚度）及其发射光谱。我们建立了基于傅里叶变换的曲线方程模型。通过建立方程并结合神经网络模型，计算出在最高热电转换效率下，GaSb电池设计多层热发射器的层数、材料和厚度以及发射光谱为4层，4层分别为GaSb、GA、GA和Sb，厚度为94.732。
模型假设：
（1）假定它严格遵守光的传播和反射规律。
（2）复合材料之间没有间隙。
（3）不考虑相对论效应。
（4）不考虑其他不可抗力因素的影响。
模型的建立与求解
基于多元回归分析模型的光谱分析
波长与折射率、波长和消光系数的关系
从题目中我们可以知道，波长的变化总是在增加，但折射率的变化却在增加和减少。分析后，我们使用一个函数来表达它。首先，分别分析波长、折射率和消光系数之间的关系，并使用软件进行它们之间的关系。
这里可以得到两个拟合函数，即波长和折射率、波长和消光系数。根据这种情况，我们建立了多元回归模型，得到波长、折射率和消光系数之间的二元函数关系。
基于傅里叶变换的曲线方程模型
我们知道锑化镓（GaSb）电池目前比较先进。假定带隙波长为1.71微米。理想化热发射器的发射光谱大致用下图中的红色虚线表示。蓝色虚线代表外部量子效率（EQE），我们希望使GaSb电池设计多层热发射器，以实现尽可能高的热电转换效率，并给出多层结构的设计参数（包括层数，材料和每层厚度），以及其发射光谱。
综上所述，（35）是推导的多层热辐射体的热电转换效率方程。不难发现，我们可以结合第三个问题直接解决，用计算机软件来解决。研究发现，层数、材料、厚度和每层以及发射光谱分别为4层，4层分别为GaSb、GA、GA和Sb，厚度为94.732。在这种情况下，GaSb电池设计的多层热发射器具有相对较高的热电转换效率。至此，问题4已解决。
论文缩略图：
程序代码：
close all;
clear all;
clc;
load asd.txt;
aa=asd;
[m,n]=size(aa);
for i=1:m
s(i,:)=aa(i,:);
F(i)=sqrt((mean(s(i,:))^2+max(s(i,:))^2)/2);
end
F;
x1=aa(:,1)';
x2=aa(:,2)';
x3=aa(:,3)';
x4=aa(:,4)';
x5=aa(:,5)';
x6=aa(:,6)';
x7=aa(:,7)';
x8=aa(:,8)';
x9=aa(:,9)';
x10=aa(:,10)';
x11=aa(:,11)';
x12=aa(:,12)';
x13=aa(:,13)';
x14=aa(:,14)';
p=[x1;x2;x3;x4;x5;x6;x7;x8;x9;x10;x11;x12;x13;x14];
t=F;
% 变量筛选 MIV算法的初步实现（增加或者减少自变量）
p=p';
[m,n]=size(p);
yy_temp=p;
% p_increase为增加10%的矩阵 p_decrease为减少10%的矩阵
for i=1:n
p=yy_temp;
pX=p(:,i);
pa=pX*1.8;
p(:,i)=pa;
aa=['p_increase' int2str(i) '=p;'];
eval(aa);
end
for i=1:n
p=yy_temp;
pX=p(:,i);
pa=pX*0.2;
p(:,i)=pa;
aa=['p_decrease' int2str(i) '=p;'];
eval(aa);
end
% 利用原始数据训练一个正确的神经网络
nntwarn off;
p=p';
% bp网络建立
net=newff(minmax(p),[5,1],{'tansig','purelin'},'traingdm');
% 初始化bp网络
net=init(net);
% 网络训练参数设置
net.trainParam.show=50;
net.trainParam.lr=0.075;
net.trainParam.mc=0.9;
net.trainParam.epochs=1000;
% bp网络训练
net=train(net,p,t);
% 变量筛选 MIV算法的后续实现（差值计算）
% 转置后sim
for i=1:n
eval(['p_increase',num2str(i),'=transpose(p_increase',num2str(i),');'])
end
for i=1:n
eval(['p_decrease',num2str(i),'=transpose(p_decrease',num2str(i),');'])
end
% result_in为增加10%后的输出 result_de为减少10%后的输出
for i=1:n
eval(['result_in',num2str(i),'=sim(net,','p_increase',num2str(i),');'])
end
for i=1:n
eval(['result_de',num2str(i),'=sim(net,','p_decrease',num2str(i),');'])
end
for i=1:n
eval(['result_in',num2str(i),'=transpose(result_in',num2str(i),');'])
end
for i=1:n
eval(['result_de',num2str(i),'=transpose(result_de',num2str(i),');'])
end
% MIV的值为各个项网络输出的MIV值,绝对值大小代表影响的相对重要性。
for i=1:n
IV= ['result_in',num2str(i), '-result_de',num2str(i)];
eval(['MIV_',num2str(i) ,'=mean(',IV,');'])
end