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
