1. Introduction

{Tower solar thermal power generation, as a new type of clean and storable energy technology, is of significant importance for countries to achieve "carbon peaking" and "carbon neutrality." Heliostats are the basic components for collecting solar energy. Through control, they can reflect sunlight onto a receiver, thereby converting solar energy into storable thermal energy. The various parameters of the heliostat field, which is composed of heliostats, have a crucial impact on optical efficiency. With the development of physical technology and data science, we can analyze the parameters of the heliostat field and their impact on thermal efficiency. By combining this analysis with actual conditions, it is possible to design a heliostat field that provides maximum working efficiency. This has profound significance for the efficient and high-quality utilization of solar energy resources. }

2. Model Assumptions

{ (1) The impact of extreme weather conditions, such as rainfall, on the operation of the heliostat field is ignored.
(2) Only mirror reflection on the heliostats is considered, ignoring refraction and diffuse reflection of light.
(3) The shadow caused by the absorber tower is ignored.
(4) The influence of terrain undulations is ignored, assuming the terrain is flat.
(5) The installation site of the heliostat field is circular. }

3. Model Establishment and Solution

3.1. Solving for Average Optical Efficiency and Output Thermal Power with Fixed Heliostat Parameters

{ To balance efficiency improvement and avoid land resource waste, we fix the heliostat dimensions at 6m * 6m and the installation height at 4m, with the initial positions of each heliostat field specified. For subsequent optimization design, it is essential first to calculate the annual average optical efficiency, annual average output power, and annual average output power per unit mirror area of the heliostat field. This paper establishes an angular solution model based on the geometric relationship between the sun, heliostats, and absorber tower. The optical efficiency and other indicators are then solved using the flat-plate projection method and the Monte Carlo method. }

3.1.1. Establishment of the Optical Indicators Model Based on the Flat-Plate Projection Method and the Monte Carlo Method

{ (1) Establishment of the Spatial Rectangular Coordinate System
A spatial rectangular coordinate system is established with the base of the absorber tower as the center (unit: m). The initial positions x and y of each heliostat field are given, with z = 4.

\( {{\alpha}}{}_{{s}} \)  represents the solar altitude angle, and  \( {{\gamma}}{}_{{s}} \)  represents the solar azimuth angle. The calculation formulas are: \[ \sin{{\alpha}}{}_{{s}}=\cos\delta\cos\varphi\cos\omega+\sin\delta\sin\varphi \] \[ \cos{{\gamma}}{}_{{s}}=\frac{\sin\delta-{{}}{}\sin{{\alpha}}{}_{{s}}\sin\varphi}{\cos{{\alpha}}{}_{{s}}\cos\varphi} \] Where:  \(\phi\)  is the local latitude, positive for north latitude;  \(\omega\)  is the hour angle of the sun, calculated as:  \( \omega=\frac{\pi}{12}(ST-12) \) . ST is the local time, and  \(\delta\)  is the solar declination angle:  \( \sin\delta=\sin\frac{2\pi D}{365}\sin(\frac{2\pi}{360}23.45) \) . Here, D is the number of days starting from the vernal equinox, counted as day 0. The vector expression of the incident sunlight is denoted as  \(\overrightarrow{i}\) , so: \[ \overrightarrow{i} =(\cos{{\alpha}}{}_{{s}}\cos{{\gamma}}{}_{{s_{{}}}},{{}}{}_{}\cos{{\alpha}}{}_{{s}}\sin\gamma{}_{{s}},\sin{{\alpha}}{}_{{s}}) \] Let the coordinates of the center of the i-th heliostat be  \(({{x}}{}_{{i}},{{y}}{}_{{i}},4) \) , and the coordinates of the center of the receiver be (0,0,80). The direction vector of the reflected sunlight from the center of the heliostat is denoted as  \(\overrightarrow{j} \) , so:  \(\overrightarrow{j}=(-{{x}}{}_{{i}},-{{y}}{}_{{i}},76) \) . (2) Solving for Cosine Efficiency Using the Angle Between Vectors
Let the angle between the incident and reflected light be 2 \(\theta\) , then:  \(\cos2\theta=\frac{\overrightarrow{i}\cdot \overrightarrow{j}}{\mid \overrightarrow{i}\mid\cdot\mid \overrightarrow{j}\mid} \) . Thus, we get:  \(\cos\theta=\sqrt{\frac{\cos2\theta+1}{2}} \) . Since  \(\theta\)  is always an acute angle,  \(\cos\theta\) >0. The cosine efficiency, which is the cosine of the angle between the incident light and the normal at the reflection point, is then:  \({{\eta}}{}_{{cos}} =\cos\theta\) . (3) Shadowing Efficiency Model Based on the Flat-Plate Projection Method and the Monte Carlo Method
To use the Monte Carlo method, an accurate pre-decision model must be established to determine whether any given heliostat will cast a shadow on the target heliostat. The main basis for establishing the pre-decision model is the principle of flat-plate projection, with appropriate improvements and innovations {[1]}. Suppose heliostat A1 potentially casts a shadow, and heliostat A2 blocks it; collectively, they are referred to as A. The target heliostat is B. Using the ray-tracing method, the projection positions of the centers of A1 and A2 on the surface coordinate system of B are determined as:  \(A1=({{x}}_{{A1}}^{{B}},{{y}}_{{A1}}^{{B}},{{z}}_{{A1}}^{{B}}) \) ,  \(A2=({{x}}_{{A2}}^{{B}},{{y}}_{{A2}}^{{B}},{{z}}_{{A2}}^{{B}}) \) . The shadowing between heliostats can be visualized as shown in Figure 1:

According to the classical flat-plate projection model, when  \(|{{{x}}_{{A1}}^{{B}}}|\)  is less than the width of B's mirror and  \( |{{{x}}_{{A2}}^{{B}}}| \)  is less than the height of B's mirror, A may cast a shadow on or block B. Since A and B are not parallel, the error is considerable. Based on this, we correct the model. Assume that the projection of heliostat A is enlarged from A'B'C'D' to a larger horizontal rectangle A''B''C''D''. Introduce a magnification factor  \(\psi\) , and when heliostat A is projected onto the surface of heliostat B, the projected figure is assumed to be  \( {{\psi}}{}_{{y}} *h\)  in length and  \( {{\psi}}{}_{{x}} *w\)  in width, where h is the height of the heliostat mirror and w is the width, both being 6m. Rewrite and simplify the judgment conditions:

where  \(\varPhi\)  is the correction coefficient, and the relationship between  \(\psi\)  and  \(\varPhi\)  is:  \({{\psi}}{}_{{x\:or\:y}}=\frac{{{\varPhi}}{}_{{x\:or\:y}}+1}{2} \)  .
Since the projection of A is irregular and its four edges are not necessarily all parallel, the calculation formula for  \(\psi\)  is: \[ {{\psi}}{}_{{x}}=\max({{x}}_{{A^{\prime}}}^{{B}}-{{x}}_{{A}}^{{B}},{{x}}_{{B^{\prime}}}^{{B}}-{{x}}_{{A}}^{{B}},{{x}}_{{C^{\prime}}}^{{B}}-{{x}}_{{A}}^{{B}},{{x}}_{{D^{\prime}}}^{{B}}-{{x}}_{{A}}^{{B}}){\frac{}{}} /(w/2)\] \[{{\psi}}{}_{{y}}=\max({{y}}_{{A^{\prime}}}^{{B}}-{{x}}_{{A}}^{{B}},{{y}}_{{B^{\prime}}}^{{B}}-{{y}}_{{A}}^{{B}},{{y}}_{{C^{\prime}}}^{{B}}-{{y}}_{{A}}^{{B}},{{y}}_{{D^{\prime}}}^{{B}}-{{y}}_{{A}}^{{B}}){\frac{}{}} /(h/2) \] Combining this formula with (1) and (2), we can pre-decide the relationship between heliostat B and heliostat A. When it is determined that shadowing or blocking occurs between two heliostats, it is necessary to record  \(({{x}}_{{A}}^{{B}},{{y}}_{{A}}^{{B}},{{z}}_{{A}}^{{B}})\)  and the values of  \(\psi\)  or  \(\varPhi\)  for subsequent calculations. Scholar Noone {[2]} proposed the boundary grid method based on computational fluid dynamics, which discretizes the four boundaries of the heliostat mirror. By deciding the positions on the heliostat where shadowing might occur, the boundary discretization is applied directly to these areas. The shadowing efficiency  \( {{\eta}}{}_{{sb}} \)  is the ratio of the unshaded or unblocked area to the total area. (4) Since the surface of the cylinder can be divided into several rectangles, the conical light reflected from the heliostat forms circular spots on the rectangular surface of the receiver. Using the previously established coordinate system, calculate the distance  \(d_{i}\)  from the center of the i-th heliostat to the center of the receiver:  \( {{d}}{}_{{i}}\:=\:|{.{\overset{arrow}{j}}|}. \) .
The light emitted by the sun is reflected by the heliostat to the center of the receiver, and the half-angle broadening of the conical light beam  \(\beta_{s} \)  is 4.65 \(mrad\) . We express the spot range (ring)  \(dS_{j}\)  formed by the heliostat microelement as a function of its half-angle  \(\beta\) : \[ dS_{j}\:=\:\pi{{d}}_{{i}}^{{2}}[{\tan^2}(\beta+d\beta)-\tan^2\beta] \] Since  \(d\beta\)  is very small, we can consider  \(d\beta = 0\) , thus:  \( dS_{j}\:=\:2\pi{{d}}_{{i}}^{{2}}\tan^{}({\beta})d\beta \) . Thus, the energy of the ring  \(dE_{ban}\)  is:  \(dE_{ban} = f(\beta)*dS_{j}\) , where  \(f(\beta)\)  represents the energy flux density function at a certain point caused by the sun's conical light beam's uneven energy distribution {[3]}. The energy of the spot is:  \({{E}}{}_{{ban}}\:=\:\int_0^{{{\beta}}{}_{{s}}}f(\beta)\cdot2\pi L^2\cdot\tan(\beta)d\beta \) . By discretely approximating  \(E_{ban}\) , we can obtain:  \( {{E}}{}_{{ban\:}}=\:\sum_{i=1}^N{\pi\tan^2}({{\alpha}}{}_{{s}})\frac{{{d}}{}_{{i}}^2}{n}\cdot{{S}}{}_{{0}}[{1-\lambda(\frac{{{\beta}}{}_{{i}}}{{{\beta}}{}_{{s}}})}^4] \) . By simulating different distances  \(d_{i}\) , we can obtain for the heliostat microelement  \(dS_{j}\) , which can be replaced by a conical beam, resulting in the energy flux density  \(S_{0}\)  of a single spot:  \( {{S}}{}_{{0}}\:=\:\:1.206\times\frac{DNI\times d{{S}}{}_{{j}}}{\pi\tan^2({{\alpha}}{}_{{s}}){{d}}{}_{{i}}} \) . By including all the spots in a rectangular grid with K rows and L columns, the energy flux density of each grid after spot overlapping can be calculated as:  \( S{{}}{}_{{m,ij}}\:=\:\sum_{k=1}^M{{f}}_{{i,j}}^{{k}}\times\frac{{{S}}{}_{{j}u}}{KL} \) . Where M is the total number of spots, and  \(S_{ju}\)  is the total area of the rectangle formed by the spots. Since the plane where the spots are located may not be parallel to the plane of the receiver, projecting  \(S{{}}{}_{{m,ij}}\)  and dividing the projected rectangle into grids can yield the energy flux density of each grid as:  \( {{S}}{}_{{r,pq\:}}\:=\:\sum{{{S}}{}_{{m,ij}}} \) . The energy received by the receiver is:  \( {{E}}{}_{{}jieshou}\:=\:\sum_{p,q}{{{S}}{}_{{r,pq}}} \) . The truncation efficiency is obtained as:  \( {{\eta}}{}_{{trunc}}\:=\:\frac{{{E}}{}_{{jieshou}}}{DNI\times A\times N\times{{\eta}}{}_{{nef}}\times{{\eta}}{}_{{sb}}} \) , where DNI is the direct normal irradiance (unit: kW/m \(^2\) ), calculated by the formula:  \( G_0 [ a+b \space exp(-\frac{c}{sin\space \alpha _s} ) ] \) , where  \(a=0.4237-0.00821(6-H)^2,b=0.5055+0.00595(6.5-H)^2,c=0.2711+0.01858(2.5-H)^2\)  ]. Here,  \({{G}}{}_{{0}} \)  is the solar constant, with a value of 1.366kW/m \(^2\) , and H is the altitude (km). (5) Calculation of Other Relevant Efficiency Coefficients The calculation formula for atmospheric transmittance  \( {{\eta}}{}_{{at}} \)  is: \[ {{\eta}}{}_{{at}}=0.99321-0.0001176{{d}}{}_{{HR}}+1.97\times10^{-8}\times{{d}}_{{HR}}^{{2}}({{d}}{}_{{HR\leq1000}}) \] where,  \( {{d}}{}_{{HR}} \)  is the distance from the center of the heliostat to the center of the receiver. The center coordinate of the i-th heliostat is ( \({{x}}{}_{{i}} \) ,  \({{x}}{}_{{i}} \) , 4), thus for the i-th heliostat:  \({{d}}{}_{{HR}}=\sqrt{{{x}}_{{i}}^{{2}}+{{y}}_{{i}}^2+76^2} \) . Using Python code, we can obtain the atmospheric transmittance  \({{\eta}}{}_{{at}} \)  for each heliostat. The reflectivity of the heliostat  \({{\eta}}{}_{{ref}} \)  is taken as 0.92. (6) Based on the above steps, we can obtain the optical efficiency of the heliostat at a fixed time point: \[ {{\eta={{\eta}}{}_{{sb}}{{\eta}}{}_{{\cos}}\eta}}{}_{{at}}{{\eta}}{}_{{trunc}}{{\eta}}{}_{{ref}} \] (7) Using the formula for calculating output thermal power:  \( {{E}}{}_{{field}}=DNI\cdot\sum_i^N{{{A}}{}_{{i}}}{{\eta}}{}_{{i}} \) . DNI is only related to the altitude, so it is the same for all heliostats.  \({{A}}{}_{{i}}\)  represents the aperture area of the i-th heliostat, which is  \(6m*6m=36m^2\) . (8) Using the average output thermal power obtained in (7), dividing by the total area of all heliostats in the heliostat field, yields the average output power per unit area of the heliostat. (9) Calculate the annual average of the above indicators:  \( G\:=\:\frac{\sum_{i=1}^{12}{{{\:{\sum_t^{}{{\omega}}{}_{{k}}}}}}}{60} \)  , where  \(G\)  is the annual average of different indicators,  \(\omega_{k}\)  represents different average efficiency indices, and  \(t\)  indicates different time points. }

3.2. Exploring the Optimal Design of Heliostat Field Using Simulation of Gravitational Search Algorithm

3.2.1. Establishment of Optimization Design Model of Heliostat Field Based on Simulation of Gravitational Search Algorithm

{ 1. Let the width of the heliostat be denoted as w, the height as h, and the installation height as H, with the coordinates of the absorber tower at (0, 0). 2. Consider the Inherent Restrictions of the Heliostat Field
(1) No heliostats are installed within the range of Am around the absorber tower, hence for any heliostat  \( ({{x}}{}_{{i}},{{y}}{}_{{i}},H) \) , it holds that  \( \sqrt{({{{x}}{}_{{i}}})^2+({{{y}}{}_{{i}}})^2}\geqslant A \)  .
(2) The mirror side length and installation height are between m and n meters, therefore:  \( m\leq w\leq n,m\leqslant h\leqslant n,m\leqslant H\leqslant n \)  .
(3) The heliostat does not touch the ground when rotating around the horizontal axis, hence the general height of the mirror cannot be greater than its installation height, i.e.:  \( \frac{h}{2}\leqslant H \)  .
(4) The distance between the bases of adjacent mirrors is more than the width of the mirror by lm, therefore, for the center coordinates of any two different mirrors  \( ({{x}}{}_{{i}},{{y}}{}_{{i}},H) \)  and  \( ({{x}}{}_{{j}},{{y}}{}_{{j}},H) \) , it holds that:  \( \sqrt{({{{x}}{}_{{i}}}-{{x}}{}_{{j}})^2+({{{y}}{}_{{i}}}-{{y}}{}_{{j}})^2} \geqslant w+l \)  .
(5) If the rated annual average output thermal power of the heliostat field reaches 60MW, then:  \( \frac{\sum_{l=1}^{60}{{E}}_{{field}}^{{l}}}{60}=60MW \) .
3. Use the Campo Layout Method to Determine the Initial Layout of the Heliostat Field {[3]} We adopt the Campo method {[4]} to establish the initial layout. The densest arrangement of heliostats in the Campo method is shown in Figure 3. % Please add the following required packages to your document preamble: %

%

Table 1: Initial Heliostat Field Numbering Results

Based on the above analysis, we take the dense heliostat field obtained from the Campo method as the initial layout field before optimization. Then, using the formulas analyzed in step 2 as constraints, we take the annual average output thermal power as the objective function. 4. Numbering the Heliostats Under these conditions, the optimization design is a high-dimensional optimization problem. To reduce the difficulty of inputting heliostat data, we use a combined numbering system for the heliostats, using the Campo method to initially number the heliostats. As previously discussed, Nrow  \( {{s}}{}_{{n}}\) , N,  \( {{R}}{}_{{n}}\)  form the initial heliostat field, where Nrow  \( {{s}}{}_{{n}}\)  can be calculated from N, and  \( {{R}}{}_{{n}}\)  will be the optimized input variable. In subsequent steps, the row radius will be optimized, so there is no need to worry about result accuracy at this stage. Using this numbering method on the initial layout, we obtain the initial design as shown in the figure. The figure shows three sets of data: the first set is the exact number of heliostat rows in each region, which is fixed; the second set is the number of heliostats in each row within the region; the third set is the radius of each row in the heliostat field. In the subsequent optimization process, we seek the optimal solution by changing the size of  \( {{R}}{}_{{n}}\) . 5. Discussion on the Principle of the Simulation-Based Gravitational Search Algorithm {[5]} (1) Principle of the Universal Gravitational Search Algorithm: We have improved the traditional gravitational search method by mimicking the characteristics of the human brain's adaptive environment and optimal selection strategy {[6]}, naming it the "Simulation-Based Gravitational Search Algorithm." First, we treat heliostats as particles in the population. In the process of finding the maximum annual average output power per unit mirror area, we consider particles with higher thermal power to have greater mass, occupying positions that receive more energy and attracting particles with lower thermal power.  \(M_{i}(t),M_{j}(t)\)  represent the masses of "particle i" and "particle j" respectively, R is the Euclidean distance between "particle i" and "particle j," and  \(\varepsilon\)  is a very small number to ensure the denominator is not zero. The gravitational constant G(t) is usually calculated using the following formula:  \( {{G}}{}_{{t}}\:=\:{{G}}{}_{{0}}\:e^{-\alpha\frac{t}{{{I}}{}_{{\max}}}} \) , where  \(\alpha\)  is a constant,  \(G_{0}\)  is the initial value of the gravitational constant, and  \(I_{max}\)  is the maximum number of iterations. (2) Analysis and Simulation Improvement of the Influence of Parameters  \(G_{0}\)  and  \(\alpha\) : Parameters  \(G_{0}\)  and  \(\alpha\)  both affect the iteration step size of the algorithm, where  \(G_{0}\)  shows a positive correlation, and  \(\alpha\)  shows a negative correlation. We improved the decay coefficient  \(\alpha\) , obtaining the following calculation formula: \[ \begin{cases}{{{\alpha}}{}(t)}\:=\:A\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:,0.6{{R}}{}_{{0}}\leq{{R}}{}_{{md}}\leq{{R}}{}_{{0}}
{\alpha(t)=\:A+B\lbrace\cos\lbrack(\frac{t}{{{I}}{}_{{\max}}})\cdot\pi+\pi\rbrack+1\rbrace}\:\:\:,0.{{2R}}{}_{{0}}\leq{{R}}{}_{{md}}\leq{{0.6R}}{}_{{0}}
{\alpha(t)\:}=\:A+\frac{B}{1+e^{-u(t)}}\:\&\:u(t)\:=\:\frac{t\cdot S}{{{I}}{}_{{\max}}}\:-\:\frac{S}{2}\:\:\:\:,0\leq{{R}}{}_{{md}}\leq{{0.2R}}{}_{{0}}\end{cases} \] A is the initial value to ensure that  \(\alpha\)  is not zero at the beginning of the algorithm iteration, B is the growth coefficient of  \(\alpha\) .  \(R_{0}\)  is the dispersion degree of the initial population, and  \(R_{md}\)  is the current population dispersion degree. Through this improvement,  \(\alpha\)  can continuously change based on the location and state during the optimization process, allowing the search to converge more quickly to the optimal direction. The improved calculation formula for the initial value of the gravitational constant  \(G_{0}\)  is: \[ {{G}}{}_{{0}}(t)\:=\:M+N\lbrack\cos(\frac{t}{{{I}}{}_{{\max}}}\cdot\pi)-1\rbrack \] where t is the current iteration number,  \(I_{max}\)  is the maximum number of iterations, M is the initial value ensuring that  \(G_{0}\)  is not zero at the start of the algorithm iteration, and N is the decay coefficient of  \(G_{0}\) .  \( \varDelta R\:=\:\varDelta{{R}}{}_{{\min}}\:+r\cdot\varDelta{{R}}{}_{{\min\:\:}} \) . Here, r is a random number between 0 and 1. (3) Set the maximum number of iterations of the algorithm to 100. The gravitational constant  \(G(t)\)  is calculated using the corrected formulas for  \(G_{0}\)  and  \(\alpha\)  obtained within 5 iterations. (4) Using the same parameters, repeat the algorithm 30 times to eliminate the influence of randomness. (5) After the iterations are completed, the algorithm can calculate the optimal value and the average of the mean values. (6) Obtain the optimal layout of the heliostat field. 6. Establish an Optimization Model for the Heliostat Field Based on the Simulation-Based Gravitational Search Algorithm (1) Set the row radius of the heliostat field as the input variable and use the annual average thermal efficiency of the heliostat field as the fitness function. Optimize the initial layout of the heliostat field using the adaptive gravitational search algorithm. (2) Set the population size Z to 40, generate the initial population, and stipulate that the minimum row spacing between consecutive rows of the corresponding heliostat field is:  \( \varDelta R\:=\:\varDelta{{R}}{}_{{\min}}\:+r\cdot\varDelta{{R}}{}_{{\min\:\:}} \) , where r is a random number between 0 and 1. (3) Set the maximum number of iterations of the algorithm to 100. Adjust the gravitational constant  \(G(t)\)  using the simulation-based strategy. (4) Repeat the above algorithm 30 times using the same parameters to eliminate randomness. (5) After each iteration, record the optimal value and average value obtained by the two algorithms, plot the iteration curves, and compare the two curves. (6) Through the simulation-based gravitational search algorithm, obtain the optimal layout of the heliostat field, and further calculate parameters and optical efficiency indices. }

4. Conclusion

The model established in this paper first fixes the parameters of the heliostats and solves for the optical efficiencies, providing conditions for the subsequent maximization of solar energy utilization. Then, by establishing an optimization design model for the heliostat field based on a simulation-based gravitational search algorithm, this paper explores the objective conditions for maximizing the utilization of the heliostat field, facilitating the green and environmentally friendly use of solar energy. We believe that this model can also be extended to the calculation of the production capacity of the heliostat field and can similarly include the construction and maintenance costs of the heliostat field, which is of great significance.