定日镜场建模过程

第一问思路

求解目标

若将吸收塔建于该圆形定日镜场中心，定日镜尺寸均为 6 m×6 m，安装高度均为4 m，且给定所有定日镜中心的位置（以下简称为定日镜位置，相关数据见附件），请计算该定日镜场的年平均光学效率、年平均输出热功率，以及单位镜面面积年平均输出热功率（光学效率及输出热功率的定义见附录）。请将结果分别按表 1 和表 2 的格式填入表格。

复述为：

在定日镜尺寸，安装高度与所有定日镜中心位置均固定的情况下，计算该定日镜场的年平均光学效率、年平均输出热功率，以及单位镜面面积年平均输出热功率。

需要设计一个算法能够通过定日镜的具体的参数去计算输出热功率这一指标。

计算过程

符号说明

符号	说明	单位
${\alpha }_s$	太阳高度角
${\gamma }_s$	太阳方位角
$\delta$	太阳赤纬
$\phi$	当地纬度
$\omega$	太阳时角
${\alpha }_p$	定日镜俯仰角
${\gamma }_p$	定日镜方位角
$R$	定日镜场的半径
$H$	当地所在海拔高度
$H_T$	集热器中心到地面的距离
$H_t$	吸收塔的高度
$H_r$	集热器的高度
$E_{\text{year}}$	定日镜场的年平均输出热功率
$\overline{E_u}$	单位面积镜面年平均输出热功率
$\text{DNI}$	法向直接辐射照度
$P_i$	第$i$面定日镜在镜面坐标系下的坐标
$V_i$	第$i$面定日镜在镜场坐标系下的坐标
$\eta$	定日镜的光学效率
${\eta }_{\text{ref}}$	镜面反射率
${\eta }_{\text{sb}}$	阴影遮挡效率
${\eta }_{\text{cos}}$	余弦效率
${\eta }_{\text{at}}$	大气透射率
${\eta }_{\text{trunc}}$	集热器截断效率
${\vec{e}}_i$	入射到定日镜的光线向量
${\vec{e}}_o$	定日镜反射的光线向量
${\vec{e}}_n$	定日镜法向量
$W_a$	定日镜镜面长度
$W_b$	定日镜镜面宽度

具体计算

太阳方位的确定

太阳高度角:
$$\sin {\alpha }_s=\cos \delta \cos \phi \cos \omega +\sin \delta \sin \phi$$

其中，${\alpha }_s$为太阳高度角，其范围为$[0,90^{\circ }],\delta$为太阳赤纬，即太阳直射纬度，$\phi$为当地纬度，$\omega$为太阳时角，即正午时圈到当地时圈的角度。

太阳赤纬:
$$\sin \delta =\sin \frac{2\pi D}{365}\sin \left(\frac{2\pi }{360}23.45\right)(2)$$
其中，$\delta$为太阳赤绯，以北绯为正向，其范围为[-23.26°,23.26°],当$\delta >0$时，太阳直射纬度在北半球，当$\delta =0$时，太阳直射绯度在赤道上，当$\delta <0$时，太阳直射纬度在南半球。$D$表示当日距离春分日的天数。

太阳方位角${\gamma }_s$
$$\cos {\gamma }_s=\frac{\sin \delta -\sin {\alpha }_s\sin \phi }{\cos {\alpha }_s\cos \phi }$$
其中$\phi$为当地纬度，北纬为正；$\omega$为太阳时角 $\omega =\frac{\pi }{12}(ST-12)$，$\delta$为太阳赤纬$\sin \delta =\sin \frac{2\pi D}{365}\sin \left(\frac{2\pi }{360}23.45\right)$

根据${\alpha }_s$，${\gamma }_s$计算得到太阳的具体方位，示意图如下图所示（A是太阳方位角，h是太阳高度角，两者可确定太阳位置）：

定日镜光学效率计算

建立定日镜坐标系

首先需要建立定日镜坐标系，包括镜场坐标系和镜面坐标系

1. 镜场坐标系，以集热器吸收塔为原点，正东方向为$X$轴正向，正北方向 为$Y$轴正向，垂直于地面向上为$Z$轴正向建立镜场坐标系$XYZ$
   

2. 对任意定日镜，令定日镜中心为原点$o$,平行于镜面左右两边向上方向为$x$轴的正方向，平行于镜面上下两边向右方向为 $y$ 轴的正方向，以垂直于镜面向外方向为$z$轴的正方向，建立镜面坐标系 $xyz$。

确定入射光线与出射光线方程

由题定义入射光线为${\vec{e}}_i\left(x_i,y_i,z_i\right)$由于在镜场坐标系下太阳方位角与高度角已经确定，故根据几何关系入射光线计算公式为写为：
$$\{\begin{array}{l}{x}_i=-\cos {\alpha }_s\cos (90°-{\gamma }_s)\\ y_i=-\cos {\alpha }_s\sin (90°-{\gamma }_s)\\ z_i=-\sin {\alpha }_s\end{array}$$
即此时入射光线方程为：${\vec{e}}_i=(-\cos {\alpha }_s\sin {\gamma }_s,-\cos {\alpha }_s\cos {\gamma }_s,-\sin {\alpha }_s)$

对于反射光线，由于题目要求入射光线经定日镜中心应反射到集热器中心，因此 在定日镜中心和集热器中心的坐标确定时，两者的连线即为反射光线。同时由于入射光线与反射光关于镜面法向量对称，基于上述分析，设任意定日镜在镜场坐标系下中心点的坐标为：$P_i(X_i,Y_i,h_i)$,其中$h_i$为定日镜中心点距离地面的高度，集热器在镜场坐标系下的坐标为：$I(0,0,H_T)$因此可确定反射光线的单位向量为：
$${\vec{e}}_o=\frac{(-X_i,-Y_i,H_T-h_i)}{\sqrt{x^2+y^2+(H_T-h_i{)}^2}}$$
联立${\vec{e}}_i$与${\vec{e}}_o$两式，根据入射光线，反射光线与平面法向量的关系可推导出镜面单位法向量的公式为：
$$\vec{e}=\frac{{\vec{e}}_o-{\vec{e}}_i}{|{\vec{e}}_o-{\vec{e}}_i|}$$

定日镜姿态确定

根据已经推导得到的定日镜法向量$\overrightarrow{e_n}\left(e_x,e_y,e_z\right)$，可对定日镜姿态进行进一步推导，以确定定日镜的俯仰角和方向角：

• 定日镜的俯仰角

  定日镜通过水半转轴的转动，其反射镜与水半地面形成的夹角即为定日镜的俯仰角${\alpha }_p$，其中根据简单几何推导可知，其值与定日镜的法向量与竖直 向形成的夹角相等，公式表达为：
  $$\cos {\alpha }_p=e_z$$

• 定日镜的方向角

  定日镜通过纵向转轴的转动，其法向量在水平面的投影与正北方向的夹角即为定日镜的方位角${\gamma }_p$，同样根据简单的几何推导，可知其方位角的表示公式为：

$$\tan {\gamma }_p=\frac{e_x}{e_y}$$

定日镜镜面反射率计算

镜面反射率与镜面的材料及清洁程度等因素有关，通过查阅资料可知，在大多数情况下，镜面反射率可取为一个常数，即：
$${\eta }_{ref}=0.92$$

定日镜大气透射率计算

当太阳光在大气中传播时，部分太阳辐射能量会被大气吸收，造成实际到达集热器 的太阳辐射能量小于定日镜反射光线的太阳辐射能量。大气透射率能够刻画反射光线的 传播过程中大气吸收对太阳辐射能传递的影响。根据附录，大气透射率是关于镜面中心至集热器中心距离的二次函数，对于任意定日镜的大气透射率计算公式为：
$${\eta }_{at}=0.99321-0.0001176d_{HR}+1.97\times 10^{-8}\times d_{HR}^2(d_{HR}⩽1000)$$
其中，$d_{HR}$为该定日镜镜面中心与集热器中心之间的距离，其公式表示为：
$$d_{\mathrm{HR},i}=|P_i-I|=\sqrt{X_i^2+Y_i+(h_i-H_T{)}^2}$$
此处$I$代指集热器中心坐标。

定日镜余弦效率的计算

入射光线倾斜入射到定日镜所在平面时，与定日镜中心法线存在着夹角 ，当 不为 0 时会产生一定的能量损失，即为余弦损失 ，减去这部分损失的能量之后即为余弦效率。查阅资料可得，$cos\theta$在数值上恰好与余弦效率相等,可得到余弦效率的计算公式为：
$${\eta }_{\cos }=\cos \theta =-\frac{{\vec{e}}_i\cdot {\vec{e}}_n}{\left|{\vec{e}}_i\right|\cdot \left|{\vec{e}}_n\right|}$$

定日镜阴影遮挡效率计算

阴影遮挡的损失主要包括两个部分：吸收塔与集热器的阴影对太阳光线的遮挡以 及前排定日镜对后排定日镜的遮挡，对其进行分开考虑。

• 吸收塔遮挡面积$S_t$

  此处假设吸收塔与集热器构成一个完整的圆柱体，则总高度$H_s$可表示为：$H_s=H_t+H_r$，底面半径为集热器半径$d$，在任意时刻，规则圆柱体的影子永 远为一矩形，该矩形的长为由圆柱高度和太阳高度角决定，如图 7 所示，矩形的宽即 为圆柱体的底面直径。因此，吸收塔与集热器的阴影区域面积$S_t$的计算公式如下：
  $$S_t=\frac{(H_t+H_r)}{\tan {\alpha }_s}\cdot d$$
  同时定义0-1变量$b_0$对定日镜进行标记，判断其是否收到吸收塔遮挡的影响，其中$b_0$的计算公式为：
  $$b_0^i=\{\begin{array}{ll}1,& x_i^2+y_i^2\le {\left(\frac{H_t+H_r}{\tan {\alpha }_s}\right)}^2\text{ 且 }{z}_0>z_1^{\prime },\\ 0,& \text{其他}.\end{array}(8)$$

• 定日镜遮挡面积

  场地中的任何定日镜都有可能受到前排定日镜的遮挡。遮挡情况主要有两种：第一种是前排定日镜 B 在阳光照射下产生的阴影，其在后排定日镜 A 上的投影形成入射光线的遮挡区域；第二种是后排定日镜 A 反射的光线被前排定日镜 B 阻挡，该区域称为反射光线的遮挡区域。如下图所示：

  

首先定义在镜面坐标系下定日镜四个顶点的坐标表示为：$P_i=\left\{\left(\frac{W_a}{2},\frac{W_b}{2},0\right),\left(\frac{W_a}{2},-\frac{W_b}{2},0\right),\left(-\frac{W_a}{2},\frac{W_b}{2},0\right),\left(-\frac{W_a}{2},-\frac{W_b}{2},0\right)\right\}$，通过坐标变换将其投影到镜场坐标系下，其转换公式为：
$$V_i=R_Z\cdot R_X\cdot \left(\begin{array}{c}{x}_n\\ y_n\\ z_n\end{array}\right)+C_B$$
其中： $R_Z$是绕$Z$轴的旋转矩阵， $R_X$是绕$X$轴的旋转矩阵，$C_B$是平移向量，表示镜面中心在镜场坐标系中的位置。依据前面推导得到的入射光线向量方程与镜面法向量方程，可求解出三者分别为：

$R_Z=\left(\begin{array}{ccc}\cos {\gamma }_p& -\sin {\gamma }_p& 0\\ \sin {\gamma }_p& \cos {\gamma }_p& 0\\ 0& 0& 1\end{array}\right)$ $R_X=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos {\alpha }_p& -\sin {\alpha }_p\\ 0& \sin {\alpha }_p& \cos {\alpha }_p\end{array}\right)$ $C_B=\left(\begin{array}{c}{X}_i\\ Y_i\\ h_i\end{array}\right)$

通过求解上式可得$V_i=(x_{i1},y_{i1},z_{i1})$,构造在入射光线方向且过点$V_i$的直线方程:
$$\frac{x-x_{i1}}{\cos {\alpha }_s\sin {\gamma }_s}=\frac{y-y_{i1}}{\cos {\alpha }_s\cos {\gamma }_s}=\frac{z-z_{i1}}{\sin {\alpha }_s}$$
联立该直线方程与定日镜平面方程，即可求解$V_i$沿着光线方向在镜场坐标系上的坐标$T_i$，再将其进行反向坐标变换，从镜场坐标系投影到被遮挡定日镜A平面,得到最终投影点$T_i^{\prime }$：
$$T_i^{\prime }=\left(\begin{array}{c}{x}_i^{\prime }\\ y_i^{\prime }\\ z_i^{\prime }\end{array}\right)=\left(\begin{array}{ccc}\cos {\gamma }_A& -\sin {\gamma }_A& 0\\ \sin {\alpha }_A\sin {\gamma }_A& \sin {\alpha }_A\cos {\gamma }_A& \cos {\alpha }_A\\ -\cos {\alpha }_A\sin {\gamma }_A& -\cos {\alpha }_A\cos {\gamma }_A& \sin {\alpha }_A\end{array}\right)(T_i-C_A)$$
最终确定入射遮挡区域为：
$$S_{m1}=\left\{\left(-\frac{W_a}{2},y_{i1}^{\prime },0\right),\left(x_{i1}^{\prime },y_{i1}^{\prime },0\right),\left(x_{i1}^{\prime },-\frac{W_b}{2},0\right),\left(-\frac{W_a}{2},-\frac{W_b}{2},0\right)\right\}$$
同时定义决策变量$b_1$,用来判断入射光线是否被遮挡，其公式为：

$b_1^i=\{\begin{array}{ll}1,& x\in \left[-\frac{W_a}{2},\frac{W_a}{2}\right],y\in \left[-\frac{W_b}{2},\frac{W_b}{2}\right],z_0>z_{i1},\\ 0,& \text{其他}.\end{array}$

即判断第i根入射光线是否被遮挡。

反射光线计算过程与入射光线计算过程类似，但是在光线方向与变换矩阵上存在不同。由前文对反射光线的计算中可得：${\vec{e}}_o=2({\vec{e}}_n\cdot {\vec{e}}_i)\times {\vec{e}}_i-{\vec{e}}_n=(a_o,b_o,c_o)$,依此推导出对任意定日镜，其反射光线方程为：
$$\frac{x-x_o^{\prime }}{a_{\mathrm{o}}}=\frac{y-y_o^{\prime }}{b_{\mathrm{o}}}=\frac{z-z_o^{\prime }}{c_{\mathrm{o}}}$$
剩下计算步骤与入射遮挡一致，最终计算得到结果为：
$$S_{m1}=\left\{\left(-\frac{W_a}{2},y_{o1}^{\prime },0\right),\left(x_{o1}^{\prime },y_{o1}^{\prime },0\right),\left(x_{o1}^{\prime },-\frac{W_b}{2},0\right),\left(-\frac{W_a}{2},-\frac{W_b}{2},0\right)\right\}$$
同样定义决策变量$b_2$,用来判断反射光线是否被遮挡，其公式为：

$b_1^i=\{\begin{array}{ll}1,& x\in \left[-\frac{W_a}{2},\frac{W_a}{2}\right],y\in \left[-\frac{W_b}{2},\frac{W_b}{2}\right],z_0>z_{o1},\\ 0,& \text{其他}.\end{array}$

根据光线追踪法，若光线传递的某一过程中被遮挡，则 不需要再计算后续过程．由此，根据上述三部分的阴影遮挡损失，得到阴影遮挡效率为：
$${\eta }_{\mathrm{sb}}=1-\sum _{i=1}^N\left[\begin{array}{c}{b}_0^i+(1-b_0^i)b_1^i+(1-b_0^i)(1-b_1^i)b_2^i\end{array}\right]/N.$$

定日镜集热器截断效率计算

由于太阳入射光并非完美的平行光，而是具有一定锥形角的一束锥形光线，因此，其反射光线也为一锥形光线，随着距离的增加，锥形光线的横截面不断变大，最后到达集热器时有可能成为一个巨大的锥形光斑，超过了集热器的吸收范围，造成了能量的溢出，进而影响了截断效率。

KaTeX parse error: Can't use function '$' in math mode at position 2: $̲n_r$为集热器接收到的光线，…
其中，$Q_r$为集热器接收能量，$Q_a$为镜面全反射能量，$Q_{sb}$为阴影遮挡损失能量。

为了简化计算，我们假设所有的光线携带的能量相同，则上述公式可化简为：
$${\eta }_{trunc}=\frac{n_r}{n_a-n_{sb}}$$
其中$n_r$为集热器接收到的光线，$n_a$为镜面反射的光线个数，$n_{sb}$为阴影遮挡的光线个数。为了能够便于计算，此处结合阴影遮挡部分计算结果，仍然使用光线追踪法进行计算，入射光线与反射光线计算公式如上$(15)(19)$所示，再分析集热器的吸收范围。

集热器接收锥形反射光线的范围$D_{{}_R}$仅为集热器圆柱外表面的曲面，即：
$$\{\begin{array}{l}{x}^2+y^2={\left(\frac{d}{2}\right)}^2\\ H_T-\frac{1}{2}{l}_T⩽z⩽H_T+\frac{1}{2}{l}_T\end{array}$$

其中，$H_T$为集热器中心点的高度，$l_T$为集热器的高，$d$ 为集热器的直径。

联立反射光线的方程和集热器接受范围的方程，并对该联立方程组是否有解进行判别，通过数学推导，得到了判别式如下：
$$\{\begin{array}{c}\Delta =4{\left(a_o{x}_o^{\prime }+b_o{y}_o^{\prime }\right)}^2-4\left(a_o+b_o\right)\left(x_o^{\prime }{}^2+y_o^{\prime }{}^2-R^2\right)\\ \frac{H_T-\frac{1}{2}{l}_T-z_o^{\prime }}{c_o}⩽\frac{-\left(a_o{x}_o^{\prime }+b_{}{y}_o^{\prime }\right)\pm \sqrt{\Delta }}{2\left(a_o+b_o\right)}⩽\frac{H_T+\frac{1}{2}{l}_T-z_o^{\prime }}{c_o}\end{array}$$
其中${\vec{e}}_o=(a_o,b_o,c_o)$为反射光线，$(x_o^{\prime },y_o^{\prime },z_o^{\prime })$为锥形入射光束在镜面上的落点在镜场坐标系中的表示坐标。

若满足以上条件，则说明方程组有解，集热器吸收了反射光线，若不满足，则说明方程组无解，集热器没有吸收到反射光线。

最后统计$n_r$，$n_a$与$n_{sb}$，带入公式$(22)$即可确定截断效率${\eta }_{trunc}$

年平均光学效率与输出热功率

通过上述对阴影遮挡效率${\eta }_{s}b$、余弦效率${\eta }_{\cos }$、大气透射率${\eta }_{a}t$、集热器截断效率${\eta }_{trunc}$以及镜面反射率${\eta }_{r}ef$等的计算，最终可求得每一块定日镜的光学效率$\eta$，根据附件其公式为：
$$\eta ={\eta }_{\mathrm{sb}}{\eta }_{\cos }{\eta }_{\mathrm{at}}{\eta }_{\mathrm{trunc}}{\eta }_{\mathrm{ref}}$$
则年平均光学功率为：
$$\overline{\eta }=\frac{{\sum }_{i=1}^{12}{\sum }_{j=1}^5{\sum }_{k=1}^N{\eta }_{ijk}}{60N}$$
其中$i$表示月份，$j$表示一天的五个时刻，$k$表示第$k$个定日镜，$N$表示定日镜的总个数

根据附件与查阅资料，某时刻的输出热功率的计算公式为：
$$E_{field}=DNI\cdot \sum _{k=1}^N{A}_k{\eta }_k$$
其中$A_k$表示第$k$面的采光面积，$DNI$表示法向直接辐射辐照度，题目中给出的计算公式为：
$$\{\begin{array}{l}DNI=G_0[a+b\exp (-\frac{c}{\sin {\alpha }_s}\left)\right]\\ a=0.4237-0.00821{\left(6-H\right)}^2\\ b=0.5055+0.00595{\left(6.5-H\right)}^2\\ c=0.2711+0.01858{\left(2.5-H\right)}^2\end{array}$$
基于上述公式，可写出年平均输出热功率$\overline{E_{\mathrm{year}}}$为：
$$\overline{E_{\mathrm{year}}}=\frac{1}{12}\sum _{i=1}^{12}\frac{1}{5}\sum _{j=1}^5\sum _{k=1}^{N}DN{I}_{ij}\cdot A_k\cdot {\eta }_{ijk}$$
故单位面积镜面年平均输出热功率$\overline{E_{\mathrm{u}}}$为：
$$\overline{E_u}=\frac{\overline{E_{\mathrm{year}}}}{{\sum }_{k=1}^N{A}_k}=\frac{{\sum }_{i=1}^{12}{\sum }_{j=1}^5{\sum }_{k=1}^{N}DN{I}_{ij}\cdot A_k\cdot {\eta }_{ijk}}{60\cdot {\sum }_{k=1}^N{A}_k}=\frac{{\sum }_{i=1}^{12}{\sum }_{j=1}^5{\sum }_{k=1}^{N}DN{I}_{ij}\cdot {\eta }_{ijk}}{60\cdot N}$$

第一问求解

第一问中如镜面反射率，大气透射率与余弦效率等可单独依靠公式进行计算，但是如阴影遮挡效率与集热器截断效率则必须通过具体的光线入射与出射进行计算，于是最终选择采取3D全仿真建模与蒙特卡洛法进行求解，通过使用Python中的Vedo来实现3D建模，并且模拟了1744*255束锥形光线入射到定日镜阵列的场景计算阴影效率和其集热器截断效率。其具体的求解步骤如下：

1. 场景构建

   吸收塔和集热器：吸收塔被建模为一个立柱形的物体，集热器位于吸收塔顶部。吸收塔的高度固定，集热器

   也被简化为一个圆柱体，放置在吸收塔的顶端。

   镜子： 每个镜子在3D空间中被表示为一个带有法向量的平面。镜子的大小、位置和朝向由输入文件中的数据提供。为了更精确地计算反射，选取了将进行划分为16*16的网格，从而模拟每个小部分入射，反射和遮挡的情况。

2. 光线的路径计算

   太阳光线： 由于锥形光线切向角不超过$4.3\times 10^{-3}$，故在其范围内随机生成一条入射光线向量，入射方向由太阳高度角（${\alpha }_s$）和方位角（${\gamma }_s$）决定。

   光线反射： 每个镜子根据其法向量将入射的太阳光线反射到吸收塔的方向。将集热器的曲面方程与反射光线方程联立，通过有无实数解可以判断该反 射光线是否落在集热器上

3. 遮挡情况计算：

   在光线传播过程中，需要判断它们是否会被其他镜子或吸收塔本身遮挡，依据光线追踪法，只要光线与其他物体相交则判断其无法继续下一步光路，将光线方程，平面镜方程与吸收塔阴影方程联立，即可判断该光线是否被遮挡。

4. 重复步骤2，模拟1744*255条锥形光线随机入射的情况，统计被遮挡光线与抵达光线的数量，带入公式即可计算出阴影效率和集热器截断效率。

下面是具体建模场景与光线入射情况可视化图，通过3D建模的方式可详细把握每一个光线入射的过程：

第一问结论

根据可视化结果，可得出如下结论：

• 平均余弦效率、平均光学效率具有季节性，平均阴影遮挡效率、平均阶段效率与 日期关系较小。
• 定日镜在夏季的输出功率大于冬季。从 1 月至 12 月，定日镜的输出功率先增大， 后减小，并在 6 月取得最大值。上述情况产生的原因可能是夏季太阳直射点位于北半球， 辐射至定日镜场的太阳辐射能比冬季多。

第二问思路

问题2要求在吸收塔位置坐标，定日镜尺寸，安装高度，定日镜数目与位置均可变的条件下，确定最佳的参数使得定日镜场的年平均输出热功率达到60MW额定功率的条件下，其单位面面积年平均输出热功率尽量大。将其视为对于单位镜面面积的年平均输出热功率的优化问题，以题目中给出可变的参数为决策变量，同时基于题意和实际情况给出多个约束条件，转换为一个多变量非线性规划问题，由于题目具有较多的决策变量与非线性的特点，经过分析，采取某某算法结合某某算法对该问题进行求解。

基于非线性单目标规划的定日镜场优化模型

目标函数

在问题1中，已经推导得到了单位镜面面积的年平均输出热功率的表达式，问题2目标是使得单位镜面面积的年平均输出热功率达到最大，故目标函数可定义为如下形式：
$$\underset{L,W_a,W_b,h_i,N,M_i}{\max }{\overline{E}}_u\left(L,W_a,W_b,h_i,N,M_i\right)=\frac{{\sum }_{i=1}^{12}{\sum }_{j=1}^5{\sum }_{k=1}^{N}DN{I}_{ij}\cdot {\eta }_{ijk}}{60N}$$
其中$L(x,y)$为吸收塔所在位置坐标，$W_a$与$W_b$为定日镜的长度和宽度，$h_i$为定日镜统一按照的高度，$N$为定日镜的数目，$M_i(x_i,y_i)$为第$i$个定日镜所在的位置坐标。

约束条件

1. 根据题意，吸收塔周围 100m 以内不设置定日镜，故吸收塔位置$L(x,y)$与任意定日镜位置$M_i(x_i,y_i)$之间的关系应该满足：

∀ i ∈ { 1 , 2 , . . . , N } , ( x i − x ) 2 + ( y i − y ) 2 ≥ 100 \forall i\in\{1,2,...,N\},\quad\sqrt{(x_i-x)^2+(y_i-y)^2}\geq100 ∀i∈{1,2,...,N},(xi​−x)2+(yi​−y)2 ​≥100

​ 其向量定义为：
∀ i ∈ { 1 , 2 , … , N } , ∥ M i − L ∥ ≥ 100 \forall i\in\{1,2,\ldots,N\},\quad\|\mathbf{M}_i-\mathbf{L}\|\geq100 ∀i∈{1,2,…,N},∥Mi​−L∥≥100

2. 镜面宽度不得小于镜面高度，故约束方程为：

$$\{\begin{array}{rl}& W_a,W_b\in [2,8]\\ & W_a⩽W_b\end{array}$$

3. 镜面边长在2m到8m之间且定日镜的安装高度在2m至6m之间，故有

$$\{\begin{array}{rl}& h_i\in [2,6]\\ & h_i>\frac{1}{2}{W}_a\end{array}$$

4. 定日镜个数，虽然题目未明确指出个数的范围，但是其数量应该控制在一个较大且合理的范围内，故定义定日镜个数约束方程为：

$$N\in [1000\text{,}5000]$$

5. 根据题意，为保持清洁，相邻的定日镜底座中心之间的距离比镜面宽度多 5m 以上，故：

∀ i , j ∈ { 1 , 2 , … , N } , i ≠ j , ∥ M i − M j ∥ ≥ 100 \forall i, j \in \{1, 2, \ldots, N\}, \quad i \neq j, \quad \|\mathbf{M}_i - \mathbf{M}_j\| \geq 100 ∀i,j∈{1,2,…,N},i=j,∥Mi​−Mj​∥≥100

6. 定日镜均位于定日镜场的范围内，即在半径为350m的圆内，故：

∀ i ∈ { 1 , 2 , . . . , N } , x i 2 + y i 2 ≤ 350 \forall i\in\{1,2,...,N\},\quad\sqrt{x_i^2+y_i^2}\leq350 ∀i∈{1,2,...,N},xi2​+yi2​ ​≤350

7. 定日镜场的额定年平均输出热功率还需要达到 60MW，故：

$$\overline{E_{\mathrm{year}}}=\frac{1}{12}\sum _{i=1}^{12}\frac{1}{5}\sum _{j=1}^5\sum _{k=1}^{N}DN{I}_{ij}\cdot A_k\cdot {\eta }_{ijk}\ge 60MW$$

综上所述，非线性单目标规划的定日镜场优化模型建立如下：
$$\begin{array}{rl}& \\ & \underset{L,W_a,W_b,h_i,N,M_i}{\max }{\overline{E}}_u\left(L,W_a,W_b,h_i,N,M_i\right)=\frac{{\sum }_{i=1}^{12}{\sum }_{j=1}^5{\sum }_{k=1}^{N}DN{I}_{ij}\cdot {\eta }_{ijk}}{60N}\\ & \{\begin{array}{llc}{W}_a,W_b\in [2,8]& & \\ W_a⩽W_b& & \\ h_i\in [2,6]& & \\ h_i>\frac{1}{2}{W}_a& & \\ N\in [1000,5000]& & \\ \sqrt{{x_i}^2+{y_i}^2}⩽R& & \\ \sqrt{{(x_i-x)}^2+(y_i-y{)}^2}⩾100& & \\ {\sqrt{{(x_i-x_j)}^2+(y_i-y_j)}}^2>W_b+5& & \\ {\overline{E}}_{year}=\frac{{\sum }_{i=1}^{12}{\sum }_{j=1}^5{\sum }_{k=1}^{N}DN{I}_{ij}\cdot A_k\cdot {\eta }_{ijk}}{60}⩾60MW& & \end{array}\end{array}$$

模型求解

由于本题决策变量较多，方程较为复杂，直接求解数值最优解难度较大，综合分析后，采取GA-BFGS混合优化算法求解上述模型最优解。具体的算法步骤如下所示：

• Step1

  给定吸收塔位置坐标、定日镜尺寸、安装高度、定日镜数目的初始值作为初始种群，根据各自数据的量纲在其范围内进行大步长遍历。

• Step2

  定义适应度函数为目标函数${\overline{E}}_u(L,W_a,W_b,h_i,N,M_i)$，根据遍历指标计算适应度高低，选择较优个体进入下一代，并且对部分个体进行变异操作，随机调整个体的某些变量值，进行适应度比较，选择较高适应度个体生存。反复迭代直至选出最终种群，作为GA部分的中间最优解。

• Step3

  从遗传算法得到的最终种群中，选择适应度最高的个体作为BFGS算法的初始解，根据${\overline{E}}_u$计算初始解的梯度与其Hessian矩阵，带入BFGS更新公式确定迭代方向和迭代步长，进行迭代得到最新解。

• Step4

  重复计算梯度和更新解的步骤，直到梯度小于设定的最终值。

第三问思路

过问题二的求解和查阅相关资料，我们发现定日镜的最优分布接近于以吸收塔塔底为圆心的同心圆上的等距分布。因此，为了简化模型，在问题三中提出了EB分区域同心圆规划策略，假设定日镜的最优排布方式为同心圆分布。相比于问题二，问题三中各定日镜的尺寸和安装高度可能会有所不同，但由于同心圆具有各向同性的特性，我们可以进一步假设每层同心圆上的所有定日镜具有相同的尺寸和安装高度。

在此基础上，问题三中的目标仍然是最大化单位镜面面积的年平均输出热功率。为此，我们在问题二的模型基础上，新增了关于每层同心圆上定日镜的尺寸和安装高度的决策变量，并更新了相关约束条件，从而构建了一个包含可变尺寸和高度的定日镜单目标优化模型。

EB分区域同心圆规划策略

EB分区域同心圆规划策略为：以吸收塔中心为圆心绘制同心圆，将定日镜场划分为多 个不同的区域，每个区域内划分为多层，定日镜底座安装在同心圆上。在每个区域内采 取不同的定日镜部署策略，使得每一层的定日镜均匀排布。

其具体布局如下图所示：

根据题目要求，相邻定日镜之间的间距至少比镜面宽度长 5m 以上，因此径向间距$\Delta D$满足条件：
$$\Delta D=W_a+5$$
在第一个部署区域，即镜场半径较小，不会出现遮挡损失的区域，定日镜的部署策略是同一层的定日镜应尽可能地紧密排列，每一层的定日镜应该以径向间距$\Delta D$为直径的圆始终保持相切。

故此时每一层之间的间距满足的关系是：
$$\Delta R_{i,j}=\Delta D$$
在第二个到第$N_c$个部署区域中，定日镜的部署策略为在同一个区域内，每一层布置的定 日镜数量相同，不同层之间定日镜交错排列；不同区域间第一层的周向间距相等。并且定义方位间距因子$A_{sf}$与周距因子$A_r$来确定每个区域的层数。

根据上述关系，镜场中各区域首环定日镜方位间距计算公式为：
$$\Delta A_{2,1}=\Delta A_{3,1}=\cdots =\Delta A_{z,1}=A_{sf}\cdot W_a$$
区域内后层定日镜的方位间距满足关系：
$$\Delta A_{i,j}=2R_{i,j}sin(\frac{\Delta {\alpha }_{i,j}}{2})$$
此时根据几何关系，除第1个区域外的区域内定日镜的排布规律如下图：

显然两层之间满足关系:
$$\Delta R_{i,j}\ge \sqrt{(\Delta D{)}^2-(\frac{\Delta A_{i,j}}{2}{)}^2}$$
基于上述同心圆规划策略，在给定具体的$W_a,W_{,}L,h_i$等决策变量的条件下可唯一确定定日镜坐标$M$这一决策变量，从而实现求解的简化。

目标函数

$$\underset{L,W_a,W_b,h_i,N,M_i}{\max }{\overline{E}}_u\left(L,W_{an},W_{bn},h_n,N,M_i\right)=\frac{{\sum }_{i=1}^{12}{\sum }_{j=1}^5{\sum }_{k=1}^{N}DN{I}_{ij}\cdot A_k\cdot {\eta }_{ijk}}{60{\sum }_{k=1}^N{A}_k}$$

其中$L(x,y)$为吸收塔所在位置坐标，$W_{an}$与$W_{bn}$为第n圈中定日镜的长度和宽度，$h_n$为第n圈定日镜的高度，$N$为定日镜的数目，$M_i(x_i,y_i)$为第$i$个定日镜所在的位置坐标。

规划变量与约束条件

在问题二的约束条件的基础上，需要修改的内容如下：

1. 第n层同心圆上的定日镜的高度满足：

$$\{\begin{array}{rl}& h_n\in [2,6]\\ & h_n>\frac{1}{2}{W}_a\end{array}$$

2. 第n层同心圆上的定日镜的尺寸满足：
   $$\{\begin{array}{rl}& W_{an},W_{bn}\in [2,8]\\ & W_{an}⩽W_{bn}\end{array}$$

3. 每一层定日镜的数量满足：

$$N_n=\frac{2\pi r_n}{d_n}$$

​ 其中，$r_n$为第$n$层同心圆的半径，$d_n$为定日镜之间的最小间距，由此可推导同心圆总层数$N_c$满足：
$$N_c=\max \left\{n\mid \sum _{i=1}^n\lfloor \frac{2\pi r_n}{d_n}\rfloor \le N\right\}$$
综上所述，问题三非线性单目标规划的定日镜场修改优化模型建立如下：
$$\begin{array}{rl}& \\ & \underset{L,W_{an},W_{bn},h_n,N,M_i}{\max }{\overline{E}}_u\left(L,W_{an},W_{bn},h_n,N,M_i\right)=\frac{{\sum }_{i=1}^{12}{\sum }_{j=1}^5{\sum }_{k=1}^{N}DN{I}_{ij}\cdot A_k\cdot {\eta }_{ijk}}{60{\sum }_{k=1}^N{A}_k}\\ & \{\begin{array}{llc}{W}_{an},W_{bn}\in [2,8]& & \\ W_{an}⩽W_{bn}& & \\ h_n\in [2,6]& & \\ h_n>\frac{1}{2}{W}_a& & \\ N\in [1000,5000]& & \\ \sqrt{{x_i}^2+{y_i}^2}⩽R& & \\ \sqrt{{(x_i-x)}^2+(y_i-y{)}^2}⩾100& & \\ {\sqrt{{(x_i-x_j)}^2+(y_i-y_j)}}^2>W_b+5& & \\ {\overline{E}}_{year}=\frac{{\sum }_{i=1}^{12}{\sum }_{j=1}^5{\sum }_{k=1}^{N}DN{I}_{ij}\cdot A_k\cdot {\eta }_{ijk}}{60{\sum }_{k=1}^N{A}_k}⩾60MW& & \end{array}\end{array}$$

模型求解

由于在考虑EB布局的情况下，模型仍然过于复杂，直接套用第二问的GA-BFGS算法，在大区域的搜索能力仍然需要进行很长时间的求解计算，为了能够简化计算且较好的搜索到最优解，本文提出了PSO-GABF算法，即在第二问的基础上加入粒子群搜索算法，引入了粒子群算法的大规模搜索能力，与GA-BFGS的小规模精细搜索相结合，从而更好的搜索到最优解，下面是模型的具体求解步骤：

• Step1 初始化粒子群

  依据EB策略和前文获得的布局结论，给出一个较好解作为每个粒子的初始位置，包括$L,W_{an},W_{b,n}{h}_n,N$以及各层同心圆的相关决策变量，并且定义每个粒子的初始速度代表其在全空间中的搜索方向和步长。

• Step2 粒子群搜索

  根据粒子当前的位置和速度，更新粒子的下一位置：$X_i^{t+1}=X_i^t+V_i^t$，并且结合结合当前速度、个体最佳位置、全局最佳位置以及随机权重，对粒子速度进行更新：$V_i^{t+1}=\omega \cdot V_i^t+c_1\cdot r_1\cdot (P_i^{best}-X_i^t)+c_2\cdot r_2\cdot (G^{best}-X_i^t)$其中，$\omega$是惯性权重，$c_1,c_2$是学习因子，$r_1,r_2$是随机数，$P_i^{best}$是个体最佳位置，$G^{best}$是全局最佳位置。再将最新坐标带入适应度函数${\overline{E}}_u$中，比较更新个体最佳和全局最佳适应度，直到迭代到规定次数或全局最佳适应度不再显著提升则停止粒子群搜索。

• Step3 遗传算法迭代

  将粒子群算法得到的最终解作为遗传算法的初始解进行迭代，给定较大步长进行步进，反复迭代直至选出最终种群，作为GA部分的中间最优解。

• Step4 BFGS算法搜索

  将GA得到的最优解作为BFGS算法的初始解，根据目标函数的梯度和Hessian矩阵，确定迭代方向，反复迭代直至收敛。

结合PSO的全局搜索和GA-BFGS的局部精细搜索，输出最终的最优解，即满足约束条件下单位镜面面积年平均输出热功率最大化的参数组合。

更新：$V_i^{t+1}=\omega \cdot V_i^t+c_1\cdot r_1\cdot (P_i^{best}-X_i^t)+c_2\cdot r_2\cdot (G^{best}-X_i^t)$其中，$\omega$是惯性权重，$c_1,c_2$是学习因子，$r_1,r_2$是随机数，$P_i^{best}$是个体最佳位置，$G^{best}$是全局最佳位置。再将最新坐标带入适应度函数${\overline{E}}_u$中，比较更新个体最佳和全局最佳适应度，直到迭代到规定次数或全局最佳适应度不再显著提升则停止粒子群搜索。

• Step3 遗传算法迭代

  将粒子群算法得到的最终解作为遗传算法的初始解进行迭代，给定较大步长进行步进，反复迭代直至选出最终种群，作为GA部分的中间最优解。

• Step4 BFGS算法搜索

  将GA得到的最优解作为BFGS算法的初始解，根据目标函数的梯度和Hessian矩阵，确定迭代方向，反复迭代直至收敛。

结合PSO的全局搜索和GA-BFGS的局部精细搜索，输出最终的最优解，即满足约束条件下单位镜面面积年平均输出热功率最大化的参数组合。