2023年国赛高教杯数学建模
A题 定日镜场的优化设计
原题再现
构建以新能源为主体的新型电力系统，是我国实现“碳达峰”“碳中和”目标的一项重要措施。塔式太阳能光热发电是一种低碳环保的新型清洁能源技术[1]。
定日镜是塔式太阳能光热发电站（以下简称塔式电站）收集太阳能的基本组件，其底座由纵向转轴和水平转轴组成，平面反射镜安装在水平转轴上。纵向转轴的轴线与地面垂直，可以控制反射镜的方位角。水平转轴的轴线与地面平行，可以控制反射镜的俯仰角，定日镜及底座示意图见图 1。两转轴的交点（也是定日镜中心）离地面的高度称为定日镜的安装高度。塔式电站利用大量的定日镜组成阵列，称为定日镜场。定日镜将太阳光反射汇聚到安装在镜场中吸收塔顶端上的集热器，加热其中的导热介质，并将太阳能以热能形式储存起来，再经过热交换实现由热能向电能的转化。太阳光并非平行光线， 而是具有一定锥形角的一束锥形光线，因此太阳入射光线经定日镜任意一点的反射光线也是一束锥形光线[2]。定日镜在工作时，控制系统根据太阳的位置实时控制定日镜的法向，使得太阳中心点发出的光线经定日镜中心反射后指向集热器中心。集热器中心的离地高度称为吸收塔高度。
现计划在中心位于东经 98.5∘，北纬 39.4∘，海拔 3000 m，半径 350 m 的圆形区域内建设一个圆形定日镜场（图 2）。以圆形区域中心为原点，正东方向为 𝑥 轴正向，正北方向为 𝑦 轴正向，垂直于地面向上方向为 z 轴正向建立坐标系，称为镜场坐标系。
规划的吸收塔高度为 80 m，集热器采用高 8 m、直径 7 m 的圆柱形外表受光式集热器。吸收塔周围 100 m 范围内不安装定日镜，留出空地建造厂房，用于安装发电、储能、控制等设备。定日镜的形状为平面矩形，其上下两条边始终平行于地面，这两条边之间的距离称为镜面高度，镜面左右两条边之间的距离称为镜面宽度，通常镜面宽度不小于镜面高度。镜面边长在 2 m 至8 m 之间，安装高度在 2 m 至 6 m 之间，安装高度必须保证镜面在绕水平转轴旋转时不会触及地面。由于维护及清洗车辆行驶的需要，要求相邻定日镜底座中心之间的距离比镜面宽度多 5 m以上。
为简化计算，本问题中所有“年均”指标的计算时点均为当地时间每月 21 日 9:00、10:30、12:00、13:30、15:00。
请建立模型解决以下问题：
问题 1 若将吸收塔建于该圆形定日镜场中心，定日镜尺寸均为 6 m×6 m，安装高度均为4 m，且给定所有定日镜中心的位置（以下简称为定日镜位置，相关数据见附件），请计算该定日镜场的年平均光学效率、年平均输出热功率，以及单位镜面面积年平均输出热功率（光学效率及输出热功率的定义见附录）。请将结果分别按表 1 和表 2 的格式填入表格。
问题 2 按设计要求，定日镜场的额定年平均输出热功率（以下简称额定功率）为 60 MW。若所有定日镜尺寸及安装高度相同，请设计定日镜场的以下参数：吸收塔的位置坐标、定日镜尺寸、安装高度、定日镜数目、定日镜位置，使得定日镜场在达到额定功率的条件下，单位镜面面积年平均输出热功率尽量大。请将结果分别按表 1、2、3 的格式填入表格，并将吸收塔的位置坐标、定日镜尺寸、安装高度、位置坐标按模板规定的格式保存到 result2.xlsx 文件中。
问题 3 如果定日镜尺寸可以不同，安装高度也可以不同，额定功率设置同问题 2，请重新设计定日镜场的各个参数，使得定日镜场在达到额定功率的条件下，单位镜面面积年平均输出热功率尽量大。请将结果分别按表 1、表 2 和表 3 的格式填入表格，并将吸收塔的位置坐标、各定日镜尺寸、安装高度、位置坐标按模板规定的格式保存到 result3.xlsx 文件中。
整体求解过程概述(摘要)
塔式太阳能光热发电是一种较为理想的、技术发展相对成熟的大规模利用太阳能发电的技术，定日镜是其收集太阳能的重要基本组件，通过数学建模对定日镜场的各项参数进行优化设计，使得单位镜面面积年平均输出热功率最大具有重大的现实意义也是我们亟待解决的问题。
针对问题一，我们构建了定日镜场年平均输出热功率模型。首先，求解太阳高度角和太阳方位角来确定每一时刻太阳所在位置；接着，通过定日镜的工作原理，由某一时刻入射光线和反射光线的方向推得定日镜在该时刻的法向，进而推得其俯仰角和方位角；然后，基于上述信息计算每一块定日镜在该时刻的光学效率，包括镜面反射率、大气透射率、余弦效率、阴影遮挡效率以及集热器截断效率，其中阴影遮挡效率由投影法求出，集热器截断效率由蒙特卡洛算法求出;最后，结合法向直接辐射照度以及定日镜场输出热功率的计算公式，代入定日镜的光学效率等数据我们求出了定日镜场在不同时刻的输出热功率，进而得到年平均输出热功率以及单位面积镜面年平均输出热功率。最终我们得到了定日镜场的年平均光学效率为0.6275，年平均输出热功率为38.295MW，单位镜面面积年平均输出热功率为06096kW/㎡，其余结果详见表1表2。
针对问题二，我们构建了定日镜场统一优化设计的单目标优化模型。我们以单位镜面面积的年平均输出热功率最大为目标函数，以定日镜场中吸收塔的位置坐标、定日镜的统一尺寸和统一安装高度、定日镜的个数以及定日镜的位置等参数作为决策变量，以定日镜场的年平均额定输出功率、相邻定日镜之间应满足的距离等要求确定了多个约束条件，建立起单目标优化模型。由于决策变量的数量很多，因此我们采用遗传算法对该单目标优化模型进行求解，最终求得了单位镜面面积年平均输出热功率最大为0.7139kW/㎡，此时年平均输出热功率为60.373MW，定日镜的分布为一圈圈同心圆，其余定日镜场的最优设计参数详见表3及result2xlsx。
针对问题三，我们构建了定日镜场非统一优化设计的单目标优化模型。相较于问题二，问题三中各个定日镜的尺寸和安装高度并不统一，决策变量的数量进一步增多，为了简化模型，我们参考问题二中求得的结论，认为在问题三中定日镜仍按照同心圆的方式进行排布，又因为同心圆具有各项同性，所以我们可以假设每一圈同心圆上的所有定日镜的尺寸和安装高度相同，因此新增的决策变量减少为每一圈同心圆对应的定日镜的尺寸和安装高度，其余目标函数和约束条件均与第二问相同。在求解时，我们在同心圆排布的基础上采用变步长的方式进行遍历求解，最终求得了单位镜面面积年平均输出热功率最大为0.7551KW/㎡，此时年平均输出热功率为60.359MW，定日镜的整体分布近似为一个抛物面，其余定日镜场的最优设计参数详见表 6及result3.xlsx。
模型假设：
1、假设天气一直保持晴朗，太阳光线不会被云层遮盖
2、假设不发生光的散射
3、假设镜面反射率可以取为常数
4、假设每条反射光线携带的能量是相同的
问题分析：
问题一的分析
问题一要求我们求解给定条件下定日镜场的年平均光学效率、年平均输出热功率以及单位镜面面积年平均输出热功率。首先，我们可以通过求解太阳高度角和太阳方位角来确定每一时刻太阳所在位置;接着，通过定日镜的工作原理我们可以由某一时刻入射光线和反射光线的方向推得定日镜在该时刻法向，进而推得其俯仰角和方位角然后，基于上述信息我们可以计算每一块定日镜在该时刻的阴影遮挡效率、余弦效率、大气透射率、集热器截断效率以及镜面反射率，从而得出定日镜的光学效率，将不同时刻不同定日镜的光学效率求和取平均即可得到年平均光学效率:最后，再结合法向直接辐射照度DNI以及定日镜场的输出热功率的计算公式，通过代入定日镜的光学效率等数据我们就可以求出定日镜场在不同时刻的输出热功率，进而求出年平均输出热功率以及单位面积镜面年平均输出热功率。
问题二的分析
问题二要求我们对，使得定日镜场的年平均输出热功率在达到额定功率60MV的条件下，单位镜面面积的年平均输出热功率尽可能大。我们将其理解为一个单目标优化问题，因此我们以单位镜面面积的年平均输出热功率最大为目标函数，以定日镜场中吸收塔的位置坐标、定日镜的统一尺寸和统一安装高度、定日镜的个数以及定日镜的位置等参数作为决策变量,再根据题意确立多个约束条件，建立起单目标优化模型。由于要求解的决策变量的数量很多，传统的遍历算法显然是行不通的，因此我们可以采用遗传算法对该单目标优化模型进行求解。
问题三的分析
相较于问题二，问题三中各个定日镜的尺寸和安装高度并不统一，这导致了决策变量的数量进一步增多，为了简化模型，我们可以参考问题二中求得的结论，认为在问题三中定日镜仍按照同心圆的方式进行排布，又因为同心圆具有各项同性，所以我们可以假设每一圈同心圆上的所有定日镜的尺寸和安装高度相同，因此新增的决策变量减少为每一圈同心圆对应的定日镜的尺寸和安装高度，其余目标函数和约束条件均与第二问相同。在求解时，我们可以在定日镜同心圆排布的基础上对每圈定日镜的尺寸及安装高度以及位置坐标采用变步长的方式进行遍历求解。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：(代码和文档not free)
%遗传算法
function f=canshul(x)
%目标函数N=2739;H-4;lw=6;1=6;w=6;x0=0;y0=-250;ST=[9,10.5,12,13.5,15];D=[306,337,0,31,61,92,122,153,184,214,245,275];xyz=[x(1:N)',x(N+1:2*N)',zeros(N,1)+x(3)];e2=zeros（12，5）；for i-1:5for j=1:12[A,B]=SUN(ST(i),D(j))；e2(j,i)=f2(xyz,N,x0,y0,A,B);endendSe2=mean(e2.'al1')；f=Se2;
end
%计算阴影遮挡效率function el=fl(xyz,1,w,N,A,B,x0,y0)Z0=84;a=[sind(B)*cosd(A),cosd(B)*cosd(A),sind(A)];tt=zeros(length(xyz),5);for i=1:size(xyz，1)m=sgrt(xyz(i,1)^2+xyz(1,1)^2+xyz(1,3));r=[-xyz(i,1)+x0,-xyz(i,2)+y0,-xyz(i,3)+z0]/sqrt((xyz(i,1)-x0)^2+(xyz(i,2)-y0)^2+(xyz(i,3)-z0)^2)n=real((r-a)/norm(r-a));An-acos(n(3));Bn=atan(n(1)/n(2));ss=[cosd(Bn],sind(Bn)*sind(An],-sind(Bn]*cosd(An];-sind(Bn],cosd(Bn)*sind (An],-cosd(Bn)*cosd(An];0,cosd(An)，sind[An];];v1=[ss*[-0.5*1,-0.5*w,o]']+xyz(i,:};v2=[ss*[0.5*1,-0.5w,0]"]'+xyz(i,:);v3=[ss*[-0.51,0.5*w,0]"]'+xyz(i,:);tt[i,1)=v1(1);tt[i,2)=v1(2);tt[i,4)=real(abs[sum[(v1-v2).*a)));tt[i,3)-real(abs[sum(v1-v3).*a)));endnum . size(tt, 1);for i = i:numareas=0;for j = i+1:numarea = rectint(tt(i,:)，tt(j,:));areas = areas + area;endtt(1,5)-areas;
endSSS-(88/tand(A)-100)*7;if sss<0sss=0;
ende1=1-(sum(tt(:,5))+SSS)/N/1/w;
end
function [e2,SS]=f2(xyz,N,x0,y0,A,B)
%余弦效率z0=84;SS=zeros(length(xyz),1);for i=1:Nb=(xyz(i,1)=x0,xyz(i,2)=y0,xyz(i,3)=z0);%反射光线b=-1.*b;a=[sind(B)*cosd(A],cosd(B]*cosd(A],sind(A]];&入射光线ta = acosd(dot{a,b)/(normla]*norm(b)));入射光线与反射光线夹角SS(i)=real(cosd(ta/2)];ende2=mean(ss);
end
