定日镜场的优化设计 摘 要
塔式太阳能光热发电是一种较为理想的、技术发展相对成熟的大规模利用太阳能 发电的技术，定日镜是其收集太阳能的重要基本组件，通过数学建模对定日镜场的各 项参数进行优化设计，使得单位镜面面积年平均输出热功率最大具有重大的现实意义， 也是我们亟待解决的问题。
针对问题一，我们构建了定日镜场年平均输出热功率模型。首先，求解太阳高度 角和太阳方位角来确定每一时刻太阳所在位置；接着，通过定日镜的工作原理，由某 一时刻入射光线和反射光线的方向推得定日镜在该时刻的法向，进而推得其俯仰角和 方位角；然后，基于上述信息计算每一块定日镜在该时刻的光学效率，包括镜面反射 率、大气透射率、余弦效率、阴影遮挡效率以及集热器截断效率，其中阴影遮挡效率 由投影法求出，集热器截断效率由蒙特卡洛算法求出；最后，结合法向直接辐射照度 以及定日镜场输出热功率的计算公式，代入定日镜的光学效率等数据我们求出了定日 镜场在不同时刻的输出热功率，进而得到年平均输出热功率以及单位面积镜面年平均 输出热功率。最终我们得到了定日镜场的年平均光学效率为0 . 6275 ，年平均输出热功 率为38 . 295MW ，单位镜面面积年平均输出热功率为0 . 6096KW/m 2 ，其余结果详见 表 1 和表 2。
针对问题二，我们构建了定日镜场统一优化设计的单目标优化模型。我们以单位 镜面面积的年平均输出热功率最大为目标函数，以定日镜场中吸收塔的位置坐标、定 日镜的统一尺寸和统一安装高度、定日镜的个数以及定日镜的位置等参数作为决策变 量，以定日镜场的年平均额定输出功率、相邻定日镜之间应满足的距离等要求确定了 多个约束条件，建立起单目标优化模型。由于决策变量的数量很多，因此我们采用遗 传算法对该单目标优化模型进行求解，最终求得了单位镜面面积年平均输出热功率最 大为0 . 7139KW/m2 ，此时年平均输出热功率为60 . 373MW ，定日镜的分布为一圈圈 同心圆，其余定日镜场的最优设计参数详见表 3 及 result2.xlsx。
针对问题三，我们构建了定日镜场非统一优化设计的单目标优化模型。相较于问 题二，问题三中各个定日镜的尺寸和安装高度并不统一，决策变量的数量进一步增多， 为了简化模型，我们参考问题二中求得的结论，认为在问题三中定日镜仍按照同心圆 的方式进行排布，又因为同心圆具有各项同性，所以我们可以假设每一圈同心圆上的 所有定日镜的尺寸和安装高度相同，因此新增的决策变量减少为每一圈同心圆对应的 定日镜的尺寸和安装高度，其余目标函数和约束条件均与第二问相同。在求解时， 我 们在同心圆排布的基础上采用变步长的方式进行遍历求解，最终求得了单位镜面面积 年平均输出热功率最大为0 . 7551KW/m2 ，此时年平均输出热功率为60 . 359MW ，定 日镜的整体分布近似为一个抛物面 ，其余定日镜场的最优设计参数详见表 6 及 result3.xlsx。
关键词：光学效率 蒙特卡洛算法 单目标优化 遗传算法 年平均输出热功率
一、 问题重述
1.1 问题背景
随着全球能源问题和环境问题日益严峻，人们迫切需要开发一种新能源来替代传 统的化石燃料，以此来达到“碳达峰”“碳中和”的目标。太阳能作为一种分布最广、 来源最稳定以及储量最大的新能源受到了人们的密切关注，而其中塔式太阳能光热发 电是一种较为理想的、技术发展相对成熟的大规模利用太阳能发电的技术， 定日镜是 其收集太阳能的重要基本组件，如何通过数学建模对定日镜场的各项参数进行优化设 计使得单位镜面面积年平均输出热功率最大成为了我们亟待解决的问题。
1.2 问题重述
定日镜的底座包括水平转轴和纵向转轴， 水平转轴上安装了平面反射镜。其中， 水平转轴可以控制反射镜的俯仰角，纵向转轴可以控制反射镜的方位角，随着太阳位 置的改变，定日镜的法向也跟着改变，使得定日镜中心反射的锥形太阳光线可以指向 集热器中心。定日镜场内排列了大量的定日镜以及一个吸收塔，集热器位于吸收塔的 顶端，其以热能的形式储存吸收得到的太阳能，然后再进一步将热能转换成电能。
问题一：在题目给定的的圆形定日镜场内，将吸收塔建于该定日镜场的中心，场 内所有定日镜的尺寸均为6m
| 符号 | 说明 | 符号 | 单位 |
|
S | 太阳高度角 |
E field | 定日镜场的年平均输出热功率 |
|
S | 太阳方位角 |
E | 单位面积镜面年平均输出热功率 |
|
| 太阳赤纬 | DNMI | 法向直接辐射照度 |
|
| 当地纬度 | PBn | 定日镜 B 在镜面坐标系下的坐标 |
| W | 太阳时角 | VBn | 定日镜 B 在镜场坐标系下的坐标 |
|
| 定日镜俯仰角 |
| 定日镜的光学效率 |
|
| 定日镜方位角 | 77ref | 镜面反射率 |
| R | 定日镜场的半径 | 7sb | 阴影遮挡效率 |
| H | 当地所在海拔高度 | 7COS | 余弦效率 |
| HT | 集热器中心到地面的距离 | 7at | 大气透射率 |
| Ht | 吸收塔的高度 | 7trunc | 集热器截断效率 |
| H, | 集热器的高度 |
6m ，安装高度也都为4m ，且附件中给出了所有定日 镜中心点的位置，计算该圆形定日镜场产生的年平均光学效率、年平均输出热功率以 及单位镜面面积年平均输出热功率。
问题二：根据题目给定的要求，定日镜场的额定功率为60MV，在定日镜尺寸以 及安装高度全部相同的条件下，给出该圆形定日镜场中吸收塔的位置以及各个定日镜 的位置、数目和尺寸和安装高度，使得该定日镜场在满足额定功率的同时，能够使单 位镜面面积年平均输出热功率尽量大。
问题三：定日镜的额定功率仍为60MV，但此时每一个定日镜的尺寸和安装高度 均可以不同，在此条件下，重新给出圆形定日镜场中吸收塔的位置以及各个定日镜的 位置、定日镜的总数目、各个定日镜的尺寸以及安装高度，使得该定日镜场在满足额 定功率的同时，能使单位镜面面积年平均输出热功率尽量大。
二、 模型假设
1 、假设天气一直保持晴朗，太阳光线不会被云层遮盖
2 、假设不发生光的散射
3 、假设镜面反射率可以取为常数
4、假设每条反射光线携带的能量是相同的
三、 符号说明
| 符号 | 说明 | 符号 | 单位 |
|
S | 太阳高度角 |
E field | 定日镜场的年平均输出热功率 |
|
S | 太阳方位角 |
E | 单位面积镜面年平均输出热功率 |
|
| 太阳赤纬 | DNMI | 法向直接辐射照度 |
|
| 当地纬度 | PBn | 定日镜 B 在镜面坐标系下的坐标 |
| W | 太阳时角 | VBn | 定日镜 B 在镜场坐标系下的坐标 |
|
| 定日镜俯仰角 |
| 定日镜的光学效率 |
|
| 定日镜方位角 | 77ref | 镜面反射率 |
| R | 定日镜场的半径 | 7sb | 阴影遮挡效率 |
| H | 当地所在海拔高度 | 7COS | 余弦效率 |
| HT | 集热器中心到地面的距离 | 7at | 大气透射率 |
| Ht | 吸收塔的高度 | 7trunc | 集热器截断效率 |
| H, | 集热器的高度 |
四、 问题分析
4.1 问题一的分析
问题一要求我们求解给定条件下定日镜场的年平均光学效率、年平均输出热功率 以及单位镜面面积年平均输出热功率。首先， 我们可以通过求解太阳高度角和太阳方 位角来确定每一时刻太阳所在位置；接着，通过定日镜的工作原理我们可以由某一时 刻入射光线和反射光线的方向推得定日镜在该时刻法向，进而推得其俯仰角和方位角； 然后，基于上述信息我们可以计算每一块定日镜在该时刻的阴影遮挡效率、余弦效率、 大气透射率、集热器截断效率以及镜面反射率， 从而得出定日镜的光学效率，将不同 时刻不同定日镜的光学效率求和取平均即可得到年平均光学效率；最后，再结合法向 直接辐射照度DNI 以及定日镜场的输出热功率的计算公式，通过代入定日镜的光学效 率等数据我们就可以求出定日镜场在不同时刻的输出热功率，进而求出年平均输出热 功率以及单位面积镜面年平均输出热功率。
4.2 问题二的分析
问题二要求我们对，使得定日镜场的年平均输出热功率在达到额定功率60MV的 条件下，单位镜面面积的年平均输出热功率尽可能大。我们将其理解为一个单目标优 化问题，因此我们以单位镜面面积的年平均输出热功率最大为目标函数，以定日镜场 中吸收塔的位置坐标、定日镜的统一尺寸和统一安装高度、定日镜的个数以及定日镜 的位置等参数作为决策变量，再根据题意确立多个约束条件，建立起单目标优化模型。 由于要求解的决策变量的数量很多，传统的遍历算法显然是行不通的，因此我们可以 采用遗传算法对该单目标优化模型进行求解。
4.3 问题三的分析
相较于问题二，问题三中各个定日镜的尺寸和安装高度并不统一，这导致了决策 变量的数量进一步增多，为了简化模型，我们可以参考问题二中求得的结论，认为在 问题三中定日镜仍按照同心圆的方式进行排布，又因为同心圆具有各项同性，所以我 们可以假设每一圈同心圆上的所有定日镜的尺寸和安装高度相同，因此新增的决策变 量减少为每一圈同心圆对应的定日镜的尺寸和安装高度，其余目标函数和约束条件均 与第二问相同。在求解时， 我们可以在定日镜同心圆排布的基础上对每圈定日镜的尺 寸及安装高度以及位置坐标采用变步长的方式进行遍历求解。
五、 模型的建立与求解
5.1 问题一模型的建立与求解
问题一要求我们求解给定条件下定日镜场的年平均光学效率、年平均输出热功率 以及单位镜面面积年平均输出热功率。首先， 我们可以通过求解太阳高度角和太阳方 位角来确定每一时刻太阳所在位置；接着，通过定日镜的工作原理我们可以由某一时 刻入射光线和反射光线的方向推得定日镜在该时刻的法向，进而推得定日镜在时刻的 俯仰角和方位角；然后，基于上述信息我们可以计算每一块定日镜在该时刻的阴影遮 挡效率、余弦效率、大气透射率、集热器截断效率以及镜面反射率， 从而得出定日镜 的光学效率，我们将不同时刻不同定日镜的光学效率求和取平均即可得到年平均光学 效率；最后，结合定日镜场的输出热功率的计算公式，代入法向直接辐射照度DNI 以 及定日镜的光学效率等数据，我们就可以求出定日镜场在不同时刻的输出热功率，进 而求出年平均输出热功率以及单位面积镜面年平均输出热功率。
图 1 问题一流程图
5.1.1 定日镜场年平均输出热功率模型
（1）太阳的高度角as 和方位角7S
随着地球的自转和公转，太阳在一年中各个时刻相对于地面的位置都在不断发生 改变，其变化轨迹被称为太阳的视运动轨迹。为了表述的方便，我们通常在地平坐标 系下来描述太阳的位置，以正东方向为
轴，垂直于地面向上的方向 为Z 轴建立坐标系，如图 2 所示。此时， 我们就可以用太阳高度角Qs 和太阳方位角7S 这两个参数来描述太阳的视运动轨迹。
图 2 地平坐标系示意图
①太阳高度角Qs
太阳高度角是指的是太阳到坐标原点的连线与地平面的夹角，通过查阅文献[1]， 我们得到了太阳高度角的计算公式为：
sinQs COS
COSPCOSW I Sin
sin
(1)
其中，as 为太阳高度角，其范围为 [0,90"] ，
为太阳赤纬，即太阳直射纬度，p 为当地纬度，W 为太阳时角，即正午时圈到当地时圈的角度。
l 太阳赤纬
太阳赤纬即太阳直射纬度，在一年中，太阳直射纬度在南北回归线之间作往返运 动，如图 3 所示，通过查阅文献[2]，我们得到了其计算公式，如下所示：
其中，
为太阳赤纬，以北纬为正向，其范围为 [-23 . 262, 23 . 26"] ，当6>0
时，太阳直射纬度在北半球，当5 0时，太阳直射纬度在赤道上，当6<0时，太阳 直射纬度在南半球。D 表示当日距离春分日的天数。
图 3 太阳赤纬在一年中的变化曲线
l 太阳时角W
太阳时角表示正午时圈到当地时圈的角度，其可以通过当地时间ST与正午 12 点 之差求出，其计算公式为：
(3)
其中，w 为太阳时角，ST为当地时间。太阳时角为 0，说明当地时间为正午，时 角为负，说明当地时间是上午，时角为正，说明当地时间是下午。
因此，当我们确定当地纬度、当天日期以及当日的时刻时， 我们就可以求出太阳 赤纬以及太阳时角，从而求出该地在特定时刻下的太阳高度角。
②太阳方位角vs
太阳方位角指的从正北方向开始沿着地平线顺时针旋转，一直到太阳在地面的投 影与坐标原点的连线的角度，在已知太阳高度角、太阳赤纬以及太阳时角时， 我们便 可以求得太阳方位角，其计算公式如下所示：
其中，7S 为太阳方位角，其范围为 [0,3609] ，
为太阳赤纬，as 为太阳高度角，
为当地纬度。
（2）定日镜的俯仰角am 和方位角7m
通过上文对太阳高度角和太阳方位角的求解，我们可以得到题目给出的定日镜场 所在地点在一年中任意一天任意一时刻太阳所在的位置，又因为定日镜在工作时，控 制系统会根据太阳所在的位置对定日镜实行实时控制，通过转动水平转轴和竖直转轴 改变定日镜的方位角和俯仰角，从而改变定日镜的法向，使得太阳中心点发出的光线 经过定日镜的中心反射后能指向集热器的中心，因此，我们可以通过对太阳中心的入 射光线以及指向集热器中心的反射光线进行分析，并结合反射定理，以此求得不同时 刻下不同定日镜的法向量，得到定日镜的俯仰角Qm 和方位角7'm 。
①定日镜场坐标系的建立
在本题中，定日镜场为一个半径350m 的圆形区域，为了更好地描述太阳、定日 镜以及集热器之间的位置关系，我们以圆心为原点，正东方向为
轴正向，正北方向 为
轴正向，垂直于地面向上为
轴正向建立坐标系，称为镜场坐标系，显然，该坐标 系与上文建立的地平坐标系保持一致。
②定日镜中心点的法向量
通过物理光学知识我们知道，当光线在镜面上发生反射时，入射角等于反射角， 因此，定日镜中心点的法向量平分入射光线和反射光线之间的夹角，通过对入射光线 和反射光线向量的计算我们便可以求得法向量
我们假设在镜场坐标系下某一条入射光线的单位向量为-ein (
in , yin , zin) ，由于在 确定时刻太阳的高度角as 和方位角7s 是确定的，所以ein (
in , yin , zin) 的计算公式为：
即：
-ein (- cosas sin7s , - cosas Cos vs , - sinas) (5)
对于反射光线，由于题目要求入射光线经定日镜中心应反射到集热器中心，因此 在定日镜中心和集热器中心的坐标确定时，两者的连线即为反射光线。在本问中，每 一块定日镜中心点在镜场坐标系下的坐标已经在附件中给出，我们假设定日镜中心点 的坐标为(a,y, x) ，又由题目可知集热器中心点的坐标为(0 , 0 , HT)，因此，反射光线的 单位向量为：
由此，我们可以推得定日镜中心点的单位法向量E(es,ey,e) ：
(7)
③定日镜的俯仰角和方位角
通过上述对定日镜中心点法向量的计算，我们便可以通过简单的几何推理得到定 日镜的俯仰角和方向角。
l 定日镜的俯仰角Qm
定日镜通过水平转轴的转动，其反射镜与水平地面形成的夹角即为定日镜的俯仰 角Qm ，如下图 4 所示，通过简单的几何推导，我们发现其值与定日镜的法向量与竖直 方向形成的夹角相等，若我们已知定日镜中心点的法向量在镜场坐标系下的坐标为 E(es,ey,e)，因此，定日镜的俯仰角 Qn 与其法向量E(es,ey,e) 的关系为：
COS Qm
e
(8)
l 定日镜的方位角7'm
定日镜通过纵向转轴的转动，其法向量在水平面的投影与正北方向的夹角即为定 日镜的方位角vm ，如图 5 所示，同样地，若我们已知定日镜中心点的法向量在镜场坐 标系下的坐标为E(e,,ey , e.)，那么定日镜的方位角Qn 与其法向量i(es, ey , e.) 的关系为：
(9)
图 4 定日镜的俯仰角示意图 图 5 定日镜的方位角示意图
（3）定日镜场的年平均光学效率
定日镜的光学效率n 是综合衡量定日镜场工作性能的重要指标，其反映的是在排 除各项影响因素之后定日镜实际能反射到集热器上的能量与其理论上能接收到的最 大太阳能之比[2] ，其值受到多个因素的影响，包括阴影遮挡效率7sb 、余弦效率ncos 、 大气透射率nat 、集热器截断效率ntrunc 以及镜面反射率7ref等，具体的计算公式如下所 示：
n
7sbncos
atnt
nc
ref (10)
①镜面反射率7ref
镜面反射率与镜面的材料及清洁程度等因素有关，通过查阅资料可知，在大多数
情况下，镜面反射率可取为一个常数，我们取其值为 0.92，即：
ref 0 . 92 (11)
②大气透射率nat
由于光线在大气中传播时，可能受到大气成分和水汽含量的影响，导致一定程度 上的能量衰减，通过查阅文献，大气透射率与镜面中心到集热器中心的距离dHR有关， 其计算公式为：
nat o . 99321 0 . ooo1176dHR
1 . 97
1o
8
dH2R (dHR 一 1000) (12)
③余弦效率
入射光线倾斜入射到定日镜所在平面时，与定日镜中心法线存在着夹角
，如图 6 所示，当
不为 0 时会产生一定的能量损失，即为余弦损失 ，减去这部分损失的能 量之后即为余弦效率，其计算公式如下所示：
nc
S
COS
(13)
因此，当入射光线与法线的夹角越大时，余弦效率越低，反之，余弦效率越高。
图 6 余弦效率计算的示意图
通过第（2）小节中的求解，我们知道了定日镜的单位法向量为
以及入射光线 的单位向量-ein ，因此，夹角
的余弦值cos
ein · -e ，从而可以求得余弦效率。
④阴影遮挡效率7sb
阴影遮挡的损失主要包括两个部分：吸收塔与集热器的阴影对太阳光线的遮挡以 及前排定日镜对后排定日镜的遮挡，我们对其进行分开考虑。
l 吸收塔与集热器的遮挡区域面积st
我们假设吸收塔与集热器构成一个完整的圆柱体，其高为吸收塔的高度Ht加上集 热器的高度H, ，底面直径即为集热器的的直径d ，在任意时刻，规则圆柱体的影子永 远为一矩形，该矩形的长为由圆柱高度和太阳高度角决定，如图 7 所示，矩形的宽即 为圆柱体的底面直径。因此，吸收塔与集热器的阴影区域面积st 的计算公式如下：
图 7 吸收塔与集热器的影子长度
该面积与圆形定日镜场的总面积之比即为吸收塔与集热器的阴影造成的遮挡损 失比7sb1 ，即：
l 定日镜的遮挡区域面积sm
场地中任意定日镜都可能受到前排定日镜对其产生的遮挡，遮挡情况主要分为两 种，如图 8 所示，一种情况是前排定日镜 B 在阳光照射下产生阴影，其阴影在后排定 日镜 A 上的投影即为入射光线的遮挡区域，另一种情况是经后排定日镜 A 反射的光 线被前排定日镜 B 阻挡，该区域称为反射光线遮挡区域。
（a）对入射光线的遮挡 （b）对反射光线的遮挡 图 8 定日镜遮挡的两种情况
i）入射遮挡区域sm
我们假设定日镜 B 投影在定日镜 A 上的遮挡区域都为规则的矩形，在求解入射 遮挡区域时，我们需要知道遮挡区域的四个顶点在定日镜 A 上的坐标，为了便于计算， 我们需要进行多次的坐标变换，具体的求解流程图如下所示：
图 9 求解入射遮挡区域的流程图
镜面直角坐标系：
我们对造成遮挡的定日镜B 的四个顶点沿入射光线向定日镜A 所在平面SA 投影， 以确定入射遮挡区域sma 。为了后续计算的方便，我们建立如下镜面直角坐标系：对于 任意一块平面反射镜，以平面镜的中心为坐标原点，过原点的法向量为
轴正向，过 原点垂直于俯仰轴偏东为
轴正向，过原点平行于俯仰轴偏北为
轴正向，如下图 10 所示：
图 10 镜面坐标系示意图
则对于定日镜B 而言，其四个顶点在镜面坐标系下的坐标分别为：
l w T
0 n 1
2 ' 2 ' ,
l w T
0 n 2
2 ' 2 ' ,
PBn
l w T
2 , 2 , 0 , n
3
l w T
2 " 2 , 0 , n
4
为定日镜B 的镜面宽度。
此外，我们还需要建立镜面坐标系与镜场坐标系之间的转换关系，以此求得这 四个顶点在镜场坐标系下的坐标。对于一个俯仰角为Qm 、方位角为7'm 的定日镜而 言，其镜面坐标系只需要先绕
轴逆时针旋转90e am ，再绕z 轴逆时针旋转7'm ，
加上对应的平移向量即可得到镜场坐标系。由此我们可以求得PBn在镜场坐标系下的 坐标VBn ：
(14)
其中，其中，n 1 , 2 , 3 , 4 ，Qm 为定日镜的俯仰角，7'm 为定日镜的方位角，CB为 定日镜B 的中心点在镜场坐标系下的坐标，Rz 为坐标轴绕Z 轴旋转产生的旋转矩
阵，Rx 为坐标轴绕z 轴旋转产生的旋转矩阵。
空间方程：
对于过四个顶点的太阳入射光线的单位向量-ein 已知，且过点VBn ，其直线方程为：
对于定日镜A 所在的平面SA ，假设其单位法向量为eA (aA , bA , CA) ，且过定日镜中 心CA (aA , YA , 2A) ，则其平面方程为：
aA (a
A) I bA (Y YA) I CA (Z ZA) 。
遮挡区域：
通过联立直线方程和平面方程，可以求解VBn沿着光线方向在平面SA 上的投影点
TBn (ZBtn , YBtn , ZBtn ) ，即：
然后，我们将求解得到的TBn进行坐标系转换，从地面坐标系转换成定日镜A 所在 平面镜的坐标系，即：
TBn 即为定日镜 B 在镜面 A 的坐标系下对其产生的入射遮挡区域的顶点坐标， 如图 11 所示：
图 11 定日镜 B 对定日镜 A 的入射遮挡区域
其中，PAn 为定日镜A 的四个顶点在镜面 A 坐标系下的坐标，AT1和AT2分别为 TB
TB2 与PA1 PA4和TB2TB3 与PA3PA4 的交点。则定日镜 B 在定日镜 A 上的入射遮挡区域
sma 即为图中阴影部分，其顶点坐标为：
ii）反射遮挡区域s 2
与入射遮挡区域的求解类似类似，我们再对定日镜 B 沿反射光线的方向进行投 影。由于已知入射光线的单位向量
以及定日镜 B 的单位法向量
，我们可以通过 数学推导得到反射光线的单位向量ES ，其公式如下：
-eout 2(-eb · -ein)
-ein -eb (aout , bout , cout) (15)
那么反射光线的方程为：
z
Bn ' y YBn '
Bn '
aout bout cout
接下来便采取与入射时同样的步骤， 即可求得VBn在反射光线方向上的投影点 KBn (ZBkn , YBkn , ZBkn ) 以及在镜面 A 坐标系下的坐标KBn (ZBkn', YBkn
, ZBkn ') ,如图 12 所示：
图 12 定日镜 B 对定日镜 A 的反射遮挡区域
定日镜 B 在定日镜 A 上的出射遮挡区域sm2 即为图中阴影部分，其四个顶点坐标 为：
iii)离散化处理定日镜 A
由于入射光线的遮挡区域和反射光线的遮挡区域可能存在重叠，导致遮挡区域的 总面积难以求解，因此我们考虑对被遮挡的定日镜 A 进行离散化处理，统计其中位于 遮挡区域内的点。
我们将定日镜 A 的镜面高度和镜面宽度按步长Ar和AY划分成为I 行J 列的网
格，示意图如图 13 所示，总计I
J 个离散点，对于其中任意一个离散点pij ，其在定 日镜 A 的镜面坐标系下的坐标为：
通过遍历i 和j ，其中位于入射遮挡区域sm1和反射遮挡区域sm2 的点的集合分别 为zmr 和zm2 ，取二者并集zm zm uzm2 ，即得到了定日镜 B 对定日镜 A 总遮挡面积 包含的离散点的个数NB ，其与离散点总个数I
J 的比值nsb2 即为定日镜的遮挡损失比，
即：
综上所述，阴影遮挡效率即为：
通过遍历所有定日镜 A 以及可能对其造成阴影遮挡的定日镜 B，分别计算器阴影 遮挡效率，然后求和取平均即可求出镜场整体的阴影遮挡损失。
图 13 离散化及遮挡区域计算的示意图
⑤集热器截断效率ntTr nc
由于太阳入射光并非完美的平行光，而是具有一定锥形角的一束锥形光线，通过 查阅文献可知[3]，该锥形角约为9 . 3mrad ，因此，其反射光线也为一锥形光线，随着 距离的增加，锥形光线的横截面不断变大，最后到达集热器时有可能成为一个巨大的 锥形光斑，超过了集热器的吸收范围，造成了能量的溢出，进而影响了截断效率，如 图 14 所示，
图 14 锥形光线示意图
由题可知，截断效率的计算公式为：
其中，Qrec为集热器接收能量，Qall为镜面全反射能量，Qsb 为阴影遮挡损失能量。 为了简化计算，我们假设所有反射光线携带的能量相同且等于qo ，集热器接收的
光线个数为nrec ，镜面反射的光线个数为nall ，阴影遮挡的光线个数为nsb ，则Qrec nrec qo ， Qall n
llqo ，Qsb nsbqo ，代入集热器截断效率的计算公式可得：
(17) 由上式可知，若要计算集热器的截断效率，我们就需要去统计集热器接收的光线
个数nrec 、镜面反射的光线个数nall 以及阴影遮挡的光线个数nsb 。由于入射锥形光线 在定日镜上的落点不确定，以及光线在锥形范围内的位置不确定，所以我们采取光线 追迹法对光线进行追踪。我们选取了大量光线照射任意一块定日镜 C，之后统计被阴 影遮挡的光线个数和集热器接收到的光线个数。
l 锥形光束中入射光线的单位向量一ein shift
对于任意一组入射到定日镜 C 上的锥形光束，其中心的入射光线与水平地面的夹 角即为太阳高度角，为了对锥形光束中的随机入射光线的进行分析，我们建立以锥形 入射光线在平面镜的落点为坐标原点、以中心入射光线的反方向为
轴正向、以过原 点垂直于俯仰轴偏东为z 轴正向，以过原点垂直于oz 向上为y 轴正向，建立入射光线 坐标系，如图 15 所示。
在入射光线坐标系下，对于锥形光束内任意一条入射光线，其方向可以由径向角
t 和切向角at 表示，切向角的范围为 [0,4 . 65mrad] 。锥形光束内任意一条入射光线
的单位向量ein shift为(- sinatsinvt , - sinatcosvt , - cos at) T ，通过坐标系转换，我们可 以将入射向量在入射坐标系下的坐标转换为地面坐标系下的坐标-e
in shift ，即：
图 15 入射光线坐标系示意图
l 锥形光束中反射光线的直线方程
对于定日镜 C，其中心点坐标为CC (ac , yC , zC) ，单位法向量为-ec (ac , bc , CC) ，我们 假设锥形光束中的随机入射光线的落点pi在镜面 C 坐标系下的坐标为(zpi , ypi , zpi) ，通 过坐标变换我们可以得到其在地面坐标系下的坐标pi
(z
pi , y
pi , zipi) 。即：
与前文类似，我们可以将入射光线的单位向量以及平面镜的法向向量通过式（15）
求出反射光线的单位向量-eout shift (aouts , bouts , couts) ，由因为反射光线经过落点p
，
则锥形光束中反射光线的方程为：
(18)
l 集热器接收范围DR所在的平面方程
集热器是一个位于塔顶的圆柱体，由于圆柱的底面被塔遮挡，所以不考虑该圆柱 的顶面和底面受到的光线，集热器接收锥形反射光线的范围DR仅为集热器圆柱外表面 的曲面，即：
(19)
其中，HT 为集热器中心点的高度，lT 为集热器的高，d 为集热器的直径。
l 判断集热器是否接收反射光线
我们联立反射光线的方程和集热器接受范围的方程，并对该联立方程组是否有解 进行判别，通过数学推导，我们得到了判别式如下：
若满足以上条件，则说明方程组有解，集热器吸收了反射光线，若不满足，则说 明方程组无解，集热器没有吸收到反射光线。
l 计算截断效率 t nc
我们随机选取nall 条光锥光束内的入射光线，统计其中被集热器接收到的光线的 个数为nTrec 以及位于阴影遮挡区域的个数为nsb ，则截断效率ntrunc 为：
⑥年平均光学效率
通过上述对阴影遮挡效率nsb 、余弦效率ncos 、大气透射率nat 、集热器截断效率 ntrunc 以及镜面反射率7ref等的计算，我们可以求得每一块定日镜的光学效率n ，则该 定日镜场的年平均光学效率
的计算公式为：
12 5 N
i 1 j 1 k 1
60N
(21)
其中，i 表示月份，j表示一天中五个特定的时刻，N 表示定日镜的总个数。
（4）定日镜场的年平均输出热功率E field及单位面积镜面年平均输出热功率 E U
定日镜场的发电原理是通过集热器吸收经定日镜反射的光线从而将光能转化成 热能，并进一步转换成电能，整个定日镜场在某一时刻下的输出热功率Efield 的计算公 式如下：
(22)
其中，Eficld 表示在某一时刻下整个定日镜场的输出热功率，DNI表示法向直接辐 射辐照度，Ak表示第k 面定日镜的采光面积，即定日镜的镜面高度乘上镜面宽度，nk 表示第k 面定日镜在该时刻下的光学效率，N 表示定日镜的总数目。
①法向直接辐射辐照度DNI
法向直接辐射辐照度指的是在地球上垂直于太阳光线的平面在单位面积以及单 位时间内接收到的太阳辐射能量，其近似计算公式如下所示：
(23)
. 3 kw/m2 ，H 为当地的海拔高度，单位为k 。
②采光面积Ak
对于每一块定日镜而言，其采光面积即为镜面高度lk乘上镜面宽度wk
Ak lk " k (24) 通常而言，定日镜的镜面高度不大于镜面宽度，但在本问中，每一面定日镜的镜
面高度都与镜面宽度相等，其值都为6m ，且所有定日镜的尺寸均相同。
③年平均输出热功率Efield
通过式（22），我们可以求得某一时刻下定日镜场的输出热功率，通过对题目所要 求的每月 21 日 5 个特定时刻所有定日镜的输出热功率的求解，再求和取平均，我们 即可得到该定日镜场的年平均输出热功率Efield ，其计算公式如下所示：
其中，DNIij表示第i 月 21 日的第j个时刻的法向直接辐射照度，nijk表示第i 月 21 日的第j个时刻下第k 面定日镜的光学效率，Ak表示第k 面定日镜的面积，N 为定
日镜的总数。
④单位面积镜面年平均输出热功率 E U
在年平均输出热功率的基础上，只需要除以所有定日镜的面积总和即可得到单位 面积镜面年平均输出热功率，又因为本问中所有定日镜的尺寸相同，因此表达式可以 写成如下形式：
5.1.2 模型的求解
通过把题目提供的数据代入到问题一建立的模型里，我们可以直接求得太阳高度 角、太阳方位角、定日镜俯仰角以及定日镜方位角等， 唯一我们无法直接求得的是集 热器的截断效率，因为截断效率需要我们随机给出锥形光线来进行计算，因此我们采 用蒙特卡洛算法来计算集热器的截断效率。
（1）蒙特卡洛算法求集热器截断效率
算法步骤
Step1：生成随机入射光线向量
对于锥形入射光线，由于其切向角不超过4 . 3
10-3 ，故我们在其范围内随机 生成一条入射光线向量。
Step2：生成随机反射点，确定反射光线方程
我们在当前的定日镜上随机取一点作为反射点，通过问题一模型我们可以算出 该定日镜平面的法向量，结合 Step1 生成的入射光线向量，我们可以求解出在该点 的反射光线方程。
Step3：判断是否相交
我们将集热器的曲面方程与反射光线方程联立，通过有无实数解可以判断该反 射光线是否落在集热器上。
Step4：计算频率，近似概率
重复上述 Step1~Step3 操作 10000 次，计算出相交次数n ，nn1000
即为该定日 镜的截断效率。
（2）求解结果及分析
通过 MATLAB 求解，我们计算出了每月 21 日平均光学效率及输出功率以及年平 均光学效率及输出功率，如下表 1、表 2 所示。
表 1 问题一每月 21 日平均光学效率及输出功率
| 日期 | 平均光学 效率 | 平均余弦 效率 | 平均阴影 遮挡效率 | 平均截断 效率 | 单位面积镜面平均输 出热功率(kw/m2) |
| 1 月 21 日 | 0.5898 | 0.7227 | 0.9877 | 0.9307 | 0.5124 |
| 2 月 21 日 | 0.6110 | 0.7436 | 0.9949 | 0.9303 | 0.5755 |
| 3 月 21 日 | 0.6309 | 0.7647 | 0.9984 | 0.9308 | 0.6272 |
| 4 月 21 日 | 0.6503 | 0.7832 | 1.0000 | 0.9352 | 0.6689 |
| 5 月 21 日 | 0.6631 | 0.7934 | 1.0000 | 0.9415 | 0.6926 |
| 6 月 21 日 | 0.6679 | 0.7965 | 1.0000 | 0.9445 | 0.7006 |
| 7 月 21 日 | 0.6630 | 0.7933 | 1.0000 | 0.9414 | 0.6923 |
| 8 月 21 日 | 0.6495 | 0.7825 | 1.0000 | 0.9349 | 0.6673 |
| 9 月 21 日 | 0.6301 | 0.7636 | 0.9983 | 0.9309 | 0.6249 |
| 10 月 21 日 | 0.6085 | 0.7409 | 0.9942 | 0.9304 | 0.5683 |
| 11 月 21 日 | 0.5878 | 0.7209 | 0.9869 | 0.9305 | 0.5063 |
| 12 月 21 日 | 0.5793 | 0.7136 | 0.9831 | 0.9300 | 0.4795 |
表 2 问题一年平均光学效率及 1 功率表
| 年平均光 | 年平均余 | 年平均阴影 | 年平均截 | 年平均输出热 | 单位面积镜面年平均 |
| 学效率 | 弦效率 | 遮挡效率 | 断效率 | 功率(MW) | 输出热功率(kw/m 2 ) |
| 0.6275 | 0.7599 | 0.9952 | 0.9342 | 38.295 | 0.6096 |
结果分析：
①对月平均光学效率的分析
我们将一年中每月 21 日平均光学效率以折线图的形式表示，如下图所示：
图 16 各个月份的平均光学效率
通过分析图像 16，我们发现有以下结论：
l 大气透射率和镜面反射率随月份没有变化，这是由于大气透射率仅与各定日镜中 心到集热器中心的距离有关，当定日镜位置确定后，平均大气透射率就确定为一 个常数。而镜面反射率按照前文的分析，也确定为一个常数。
l 平均阴影遮挡效率极高，接近于 1。这是由于定日镜之间距离间隔较大，阴影遮 挡损失极小，阴影遮挡效率仅受到吸热塔对定日镜的阴影遮挡造成的影响。
l 平均截断效率相对较高，在 6 月时最高达到0.94 左右，这是由于定日镜的大小于 集热器的竖截面矩形的大小近似，溢出的面积较小，截断损失较小。此外，我们 发现平均截断效率关于 6 月对称分布，这是由于光斑大小与入射光线的方向有关， 即与太阳高度角有关，而太阳高度角关于 6 月对称分布。
l 平均余弦效率相对较低，在 6 月时最高，但仅为 0.8 左右，这是由于平均余弦效 率取决于太阳入射角的余弦值，考虑到定日镜与吸热塔的距离等因素，太阳入射 角往往较大，余弦值往往较小，因此平均余弦效率相对较小。与平均截断效率类 似，平均余弦效率关于 6 月对称分布，这也是太阳高度角关 6 月对称分布导致的。
l 平均光学效率的变化情况与平均余弦效率的分布类似，这是由于平均余弦效率之 外的其他效率都较高，接近于 1 ，因此我们认为仅有平均余弦效率对平均光学效 率造成了显著影响。
②对不同定日镜年平均光学效率的分析
此外，我们还将所有定日镜从最内层最东侧的定日镜开始沿逆时针确定序号，计 算了每块定日镜的年平均光学效率，如下图 17 所示。
定 日镜序号
图 17 每块定日镜的年平均光学效率
图中的每个波峰表示定日镜分布的其中一层，我们发现由内层到外层，定日镜的 年平均光学效率逐渐下降，此外，我们还发现每一层的平均光学的效率分布大致相同。 因此，我们取其中的最内层的定日镜进行局部范围分析，如下图 18 所示：
5 10 15 20 25 30 35 40 45 50 55
定 日镜序号
图 18 最内层定日镜的年平均光学效率
通过分析上图 18，我们发现对于最内层的定日镜而言，其平均光学效率的分布近 似于正弦分布，且位于北方的定日镜光学效率整体较高，在序号 15 处即在正北方向 处达到最大。而在东西方向上，平均光学效率呈对称分布，对整体光学效率的影响相 互抵消。我们由此推测， 如果需要对定日镜的分布进行优化，应将定日镜尽可能摆放 在吸热塔的北方。
为了进一步验证该结论，我们选取了春分日 3 月 21 日上午 10 点时，所有定日镜 的年平均光学效率在地面二维坐标系上的分布情况，如下图。
图 19 问题一中春分日 3 月 21 日上午 10 点时所有定日镜的年平均光学效率
可以发现太阳入射光线在地面上的落点处附近的年平均光学效率相对较高，由该 地的纬度为北纬 39.4°, 在北回归线以北，太阳始终位于吸热塔南方，所以太阳入射光 线落点始终位于吸热塔北方，因此偏北的定日镜年平均光学效率相对较高，与我们的 推测一致。如果需要对定日镜的分布进行优化，应将定日镜尽可能摆放在吸热塔北方。
5.1.3 模型检验
通过查阅文献可知，单位镜面面积年平均输出功率和定日镜场年平均输出功率与 镜面的尺寸存在相关性，由于本问中定日镜的高度和宽度相等，因此我们只需要对镜 面宽度进行一定范围内的调整，计算调整后的单位镜面面积年平均输出功率和定日镜 场年平均输出功率，如下图 20 所示：
|
|
| ||||
|
| |||||
|
| |||||
|
| |||||
|
| |||||
4 4.5 5 5.5 6 6.5 7 4 4.5 5 5.5 6 6.5 7
定日镜宽度I(m) 定日镜宽度I/(m)
图 20 年平均输出功率与定日镜宽度的关系
通过对上图分析，我们发现随着镜面宽度增大，单位镜面面积年平均输出功率减 小，年平均输出功率增大，这与实际情况相符合，这证明了模型的准确性。
此外，我们发现当定日镜宽度增加了 3 米时，单位镜面面积年平均输出功率减少 了 0.12 千瓦每平方米，年平均输出功率增加了 30 兆瓦，均存在较大的变化幅度，说 明单位镜面面积年平均输出功率和年平均输出功率对镜面宽度敏感，证明了我们模型 的灵敏性。
5.2 问题二模型的建立与求解
问题二要求我们对定日镜场的布局进行优化设计，使得定日镜场的年平均输出热 功率在达到额定功率60MV的条件下，单位镜面面积的年平均输出热功率尽可能大， 因此，我们以单位镜面面积的年平均输出热功率最大为目标函数，以定日镜场中的吸 收塔位置坐标，定日镜的统一尺寸，统一安装高度，定日镜的个数以及定日镜的位置 等定日镜场的分布参数为决策变量，并根据题意确定了多个约束条件，以此建立起单 目标规划模型。由于决策变量的数量很多，因此我们采用遗传算法对该单目标优化模 型进行求解。
图 21 问题二流程图
5.2.1 基于单目标规划的统一定日镜场优化设计模型
（1）目标函数
在第一问中，我们给出了单位镜面面积年平均输出热功率 E . 的计算公式，我们 的目标是使得单位镜面面积年平均输出热功率 E U 尽可能大，因此目标函数可以写成 如下形式：
其中，r,l , w,h,N,Ia为规划变量，下面将对其进行详细阐述。
（2）规划变量与约束条件
①吸收塔的位置坐标
在问题一中，吸收塔在镜场坐标系下的位置坐标为(0 , 0) ，此时我们可以计算定 日镜场的年平均光学效率
。参照问题一的模型，我们在不改变问题一中定日镜场的 其他的参数的情况下，对吸收塔的位置在东西方向和南北方向上进行移动，计算其在 不同位置对应的
，以此探究吸收塔的位置变化对年光学效率的影响。我们发现年光 学效率随着吸收塔移动产生的变化如下图22 所示。
0
东西方向位移(m)
0
南北方向位移(m)
图 22 光学效率随着吸收塔移动产生的变化
由上图可知，吸收塔在向不同方向移动时会对
产生不同程度的影响。当吸收塔 在东西方向移动时，无论是向东还是向西都会导致
的下降，其变化情况关于原点对 称，且变化幅度极小，这是由于太阳在一天中的运动在东西方向上是对称的，若吸收 塔向东偏移，在上午时，其东侧的定日镜的光学效率较低，西侧的较高，在下午时， 东侧的定日镜光学效率较高，西侧的较低，上午和下午的的变化相抵消，导致
并没 有与位于坐标原点处时发生显著变化。若吸收塔向西偏移， 同理
也不会发生显著变 化。
当吸收塔在南北方向移动时，我们发现越靠南方
越大， 反之，
越小， 且
随 吸收塔位置改变而发生的改变较为显著，因此，为了简化模型，我们认为吸收塔仅在
南北方向发生偏移，我们设偏移量为r ，则偏移后吸收塔的坐标为o(0 , r)，由于吸收 塔发生偏移后仍应位于镜场范围内，因此偏移量
的绝对值应小于镜场的半径
，即
rc [ -R, R]
②定日镜尺寸lxw
在问题二中，每块定日镜的尺寸是相同的，即每块定日镜都有相同的高度l 和宽度 w ，且满足高度l 和宽度w 在 2 米到 8 米之间，并且高度不大于宽度，即：
l, w C
③定日镜安装高度h
在本问中，每块定日镜的安装高度h 是相同的，且满足镜面高度h 在 2 米到 6 米 之间，并且镜面在绕水平轴旋转时镜面不会触地，即安装高度大于镜面高度的一半。
④定日镜个数N
定日镜的个数是定日镜场布局时的一个重要参数，它会影响到定日镜排布的紧凑 程度，进而影响定日镜场的年平均总效率Efield，也会影响到单位镜面面积的年平均输 出热功率Eu ，因此定日镜个数N 也应作为决策变量，我们令其在一个足够大的范围内 变化，即：
N
[1000 , 5000]
⑤定日镜位置Ia (
a , y.)
设第a 个定日镜的位置为点Ia ，其坐标为(a, , y.) ，由于定日镜位于镜场范围内，因 此点I. 到原点的距离应小于镜场半径R ，并且，由于吸收塔 100 米范围内不安装定日 镜，因此点I. 到吸收塔的距离应大于100m ，此外，定日镜a 与相邻定日镜b 的距离比 镜面宽度多 5m 以上，那么约束条件有：
⑥额定年平均输出热功率Efi
ld
问题二中要求额定年平均输出热功率Efi
ld 达到 60MW，即：
（3）模型给出
综上，我们给出基于单目标规划的统一定日镜场优化设计模型，如下所示：
max
r, l , W , h,N,I
R
2m
l
he
(28)
5
t
R
2
(ya yb)
5.2.2 模型的求解
（1）遗传算法
由于本问的决策变量很多，传统的遍历方法时间成本较高，因此我们考虑用遗传 算法对模型进行求解，具体的算法步骤如下所示：
算法步骤
Step1：确定部分参数初始值
我们给定吸收塔位置坐标、定日镜尺寸、安装高度、定日镜数目的初始值， 并
根据各自数据的量纲在其范围内进行大步长遍历。
Step2：使用遗传算法计算定日镜位置
对于当前 Step1 确定的参数，我们采用遗传算法进行求解，对于一次遍历给定 的参数，我们计算其在约束条件下的最大单位镜面面积年平均输出热功率，若大于 当前最优平均输出热功率，则替换之，并记录下最优平均输出热功率下的定日镜位 置坐标。
Step3：缩小步长精细遍历
将步长逐步缩小遍历，重复 Step1 与 Step2 操作，最终确定满足约束条件下的 最大单位镜面面积年平均输出热功率，得到定日镜场的各项参数。
（2）求解结果及分析
通过上述求解步骤，我们计算出了使得单位镜面面积年平均热输出功率最优时的 定日镜场的各个参数，如下表 3 所示：
表 3 问题二设计参数表
吸收塔位置坐标 定日镜尺寸 定日镜安装高度 定日镜总面数 定日镜总面积
(0 , - 250) 5 . 55m
5 . 55m 5m 2739 84368 . 0475m 2
在此最优参数下，我们可以计算最优月平均光学效率及输出功率和最优年平均光 学效率和输出功率，如下表 4、表 5 所示。
表 4 问题二每月 21 日平均光学效率及输出功率
| 日期 | 平均光学 效率 | 平均余弦 效率 | 平均阴影 遮挡效率 | 平均截断 效率 | 单位面积镜面平均输 出热功率(kw/m2) |
| 1 月 21 日 | 0.7669 | 0.9123 | 0.9940 | 0.9627 | 0.6663 |
| 2 月 21 日 | 0.7579 | 0.9044 | 0.9980 | 0.9558 | 0.7137 |
| 3 月 21 日 | 0.7432 | 0.8891 | 0.9998 | 0.9517 | 0.7388 |
| 4 月 21 日 | 0.7229 | 0.8652 | 0.9998 | 0.9513 | 0.7436 |
| 5 月 21 日 | 0.7089 | 0.8432 | 0.9998 | 0.9573 | 0.7404 |
| 6 月 21 日 | 0.7033 | 0.8338 | 0.9998 | 0.9604 | 0.7377 |
| 7 月 21 日 | 0.7092 | 0.8434 | 0.9998 | 0.9573 | 0.7406 |
| 8 月 21 日 | 0.7241 | 0.8664 | 0.9998 | 0.9515 | 0.7440 |
| 9 月 21 日 | 0.7443 | 0.8901 | 0.9998 | 0.9520 | 0.7382 |
| 10 月 21 日 | 0.7590 | 0.9057 | 0.9976 | 0.9562 | 0.7088 |
| 11 月 21 日 | 0.7663 | 0.9128 | 0.9936 | 0.9618 | 0.6600 |
| 12 月 21 日 | 0.7673 | 0.9142 | 0.9916 | 0.9635 | 0.6352 |
表 5 问题二年平均光学效率及功率表
| 年平均光 | 年平均余 | 年平均阴影 | 年平均截 | 年平均输出热 | 单位面积镜面年平均 |
| 学效率 | 弦效率 | 遮挡效率 | 断效率 | 功率(MW) | 输出热功率(kw/m 2) |
| 0.7394 | 0.8817 | 0.9978 | 0.9568 | 60.373 | 0.7139 |
结果分析：
①对定日镜分布的分析
首先，我们画出了此时定日镜在定日镜场内的分布，如下图 23 所示：
图 23 问题二中定日镜在定日镜场中的分布
通过观察图像我们发现，最优的吸收塔位置位于定日镜场的最南端，这与我们在 第一问中的分析是一致的，定日镜在定日镜场中的排布是围绕吸收塔的许多个同心圆， 定日镜的尺寸相较于问题一中较小，但定日镜的数目则比问题一中要多，总体而言， 定日镜的排布比问题一中更为密集。
②对月平均光学效率的分析
同样地，我们画出了每月 21 日平均光学效率的折线图，如下图 24 所示;
图 24 问题二每月 21 日平均光学效率
通过观察图像，我们发现以下结论：
l 平均阴影遮挡效率仍然极高，接近于 1 。我们知道，本问中定日镜的排布密度要 大于问题一中的初始排布密度，但阴影遮挡效率缺和第一问几乎一样，这说明定 日镜的最优排布方式充分考虑了板与板之间的遮挡，能使得排布密度增加的同时 仍然保证阴影遮挡效率保持在较高水准。
l 平均余弦效率在 6 月份最低，为 0.83 左右。这一点恰恰与第一问时的情况相反， 在第一问中平均余弦效率在 6 月份最高，但只有 0.8 左右，比这一问的最低余弦 效率还要低，这说明我们在本问中对定日镜场的优化布置是合理的，而之所以在 6 月份余弦效率最低，我们推测是因为 6 月份的太阳高度角最大，入射光线与平 面镜法线之间的夹角最大，导致了余弦损失最大，余弦效率最低。
③对不同定日镜的年平均光学效率的分析
同样地，我们将所有定日镜从最内层最东侧的定日镜开始沿逆时针确定序号，计
算了每块定日镜的年平均光学效率，如下图 25 所示。
定 日镜序号
图 25 问题二中不同定日镜的年平均光学效率
通过观察图像，我们发现从最内侧到最外侧，定日镜的年平均光学效率仍然逐渐 减小，但减小的幅度很小，内外圈的差异并不如第一问中一样显著，这同样说明了最 优排布下能使得所有定日镜尽可能地达到最优光学效率。
此外，我们计算了春分日 3 月 21 日 10 点时，所有定日镜的年平均光学效率，并 作出了其在地面上的分布情况，如下图 26 所示：
图 26 春分日 3 月 21 日 10 点时所有定日镜的年平均光学效率
④年均光学效率及输出功率相较于第一问的优化率
为了能直观地体现我们的优化程度，我们将问题一的年均光学效率等参数与优化 后的参数以直方图的形式进行对比，如下图 27 所示：
图 27 问题一与问题二年均光学效率的对比
通过观察图像 27，我们发现各项指标都在问题二的优化布局后得到了提升，此外， 通过计算，我们得到了年平均光学效率的优化率为 17.86%，年平均输出热功率的优化 率为 57.65%，单位镜面面积的年平均输出热功率的优化率为 17.11%，有力地说明了 我们对定日镜场的优化设计是有效且合理的。
5.2.3 模型检验
由于实际情况中可能存在的随机误差，我们在满足约束条件的情况下，随机生成 了 100 组可行解，计算其单位镜面面积年平均输出热功率，作出针状图，与我们求得 的最优单位镜面面积年平均输出热功率进行比较，如下图 28 所示：
随机次数序号
图 28 随机生成的 100 组可行解与最优解的比较
通过观察图像，我们发现所有随机可行解都小于我们求得的最优解，因此证明了 我们的最优解是准确的，我们的模型是可靠的。
5.3 问题三模型的建立与求解
通过问题二的求解，我们发现定日镜的最优分布近似于以吸收塔塔底为圆心的同 心圆上的等距分布，因此为了简化模型，在本问中我们可以适当沿用此结论，认为在 问题三中最优定日镜的排布方式也为同心圆。相较于问题二， 问题三中各个定日镜的 的尺寸以及安装高度可能不同，但由于同心圆具有各项同性的性质，因此我们可以假 设每一层同心圆上的所有定日镜的尺寸和安装高度相同。然后，在问题二的基础上， 我们同样以单位镜面面积的年平均输出热功率最大为目标函数，在决策变量里新增了 每层同心圆上的定日镜的尺寸和安装高度，并更新了部分约束条件，从而建立起定日 镜尺寸和高度可变的单目标规划模型。
图 29 问题三流程图
5.3.1 基于单目标规划的非统一定日镜场的优化设计模型
（1）目标函数
与问题二类似，本问我们要求解最优的单位镜面面积年平均热输出功率E.n ，但与 第二问唯一不同的地方在于本问中定日镜的尺寸并不都一样，这就导致了在计算定日 镜总面积时无法化简，则本问中目标函数为:
（2）规划变量与约束条件
相较于问题二，问题三的规划变量新增了每一层同心圆上的定日镜尺寸ln
wn 和 定日镜安装高度hn ，则相应的约束条件应更新为：
其中，ln n 和hn 分别为第n 层同心圆的上定日镜的尺寸和安装高度，nu 为定 日镜场中分布的同心圆的层数。其余的规划变量和约束条件全部与问题二中保持一致。
（3）模型给出
综上，我们给出基于单目标规划的非统一定日镜场的优化设计模型，如下所示：
r
[ -R, R]
ln , n C
ln ≤wn (n 1 , 2,…,num) hn
[2m, 6m]
(29)
2 n
N
[1000 , 5000]
Vu,2 1
2
y
ya yb
12 5 N
DN
ij · Ak · n
jk
i 1 j 1 k 1
5.3.2 模型的求解
（1）算法步骤
通过问题二的求解，我们发现定日镜在最优情况下大致呈同心圆分布，为了更好 的利用空间效率，我们认为在同一圆周上，定日镜宽度加 5 为定日镜间距，两个相邻 圆周间间距为外侧圆周上定日镜宽度加 5，以此简化计算。对于每一圈定日镜的尺寸 和安装高度，我们可以采用变步长遍历的方法进行求解，具体的算法步骤如下所示：
算法步骤
Step1：确定参数初始值
我们给定吸收塔位置坐标、各圆周上定日镜尺寸、各圆周上安装高度的初始值， 并根据各自数据的量纲在其范围内进行大步长遍历。
Step2：计算最优单位镜面面积年平均输出热功率
对于给定的参数，我们可以确定定日镜的位置分布以及定日镜数目，通过问题 一的模型可以计算出当前参数下的单位镜面面积年平均输出热功率，若大于当前最 优平均输出热功率，则替换之，并更新最优平均输出热功率下各个参数数据。
Step3：缩小步长精细遍历
在当前各将步长逐步缩小遍历，重复 Step1 与 Step2 操作，最终确定满足约束 条件下的最大单位镜面面积年平均输出热功率， 得到吸收塔位置坐标、各圆周上定 日镜尺寸、各圆周上安装高度，进一步得出定日镜的位置分布以及定日镜数目。
（2）求解结果及分析
通过上述求解步骤，我们计算出了使得单位镜面面积年平均热输出功率最优时的 定日镜场的各个参数，如下表 6 所示：
表 6 问题三设计参数表
| 吸收塔位置坐标 | 定日镜总面数 | 定日镜总面积 |
|
(0 - 250) | 2315 | 87633 . 4929m 2 |
在此最优参数下，我们可以计算最优月平均光学效率及输出功率和最优年平均光 学效率和输出功率，如下表 7、表 8 所示。
表 7 问题三每月 21 日平均光学效率及输出功率
| 日期 | 平均光学 效率 | 平均余弦 效率 | 平均阴影 遮挡效率 | 平均截断 效率 | 单位面积镜面平均输 出热功率(kw/m2) |
| 1 月 21 日 | 0.8071 | 0.9698 | 0.9931 | 0.9521 | 0.7012 |
| 2 月 21 日 | 0.7984 | 0.9635 | 0.9972 | 0.9443 | 0.7519 |
| 3 月 21 日 | 0.7850 | 0.9498 | 0.9989 | 0.9402 | 0.7803 |
| 4 月 21 日 | 0.7664 | 0.9270 | 0.9987 | 0.9406 | 0.7884 |
| 5 月 21 日 | 0.7535 | 0.9056 | 0.9987 | 0.9467 | 0.7870 |
| 6 月 21 日 | 0.7496 | 0.8964 | 0.9987 | 0.9515 | 0.7863 |
| 7 月 21 日 | 0.7538 | 0.9059 | 0.9987 | 0.9467 | 0.7871 |
| 8 月 21 日 | 0.7670 | 0.9282 | 0.9987 | 0.9401 | 0.7880 |
| 9 月 21 日 | 0.7856 | 0.9507 | 0.9989 | 0.9400 | 0.7792 |
| 10 月 21 日 | 0.8005 | 0.9646 | 0.9967 | 0.9461 | 0.7476 |
| 11 月 21 日 | 0.8076 | 0.9701 | 0.9927 | 0.9529 | 0.6957 |
| 12 月 21 日 | 0.8080 | 0.9710 | 0.9907 | 0.9544 | 0.6689 |
| 表 8 问题三年平均光学效率及功率表 | |||||
| 年平均光 | 年平均余 | 年平均阴影 | 年平均截 | 年平均输出热 | 单位面积镜面年平均 |
| 学效率 | 弦效率 | 遮挡效率 | 断效率 | 功率(MW) | 输出热功率(kw/m 2) |
| 0.7818 | 0.9418 | 0.9968 | 0.9462 | 60.359 | 0.7551 |
结果分析：
①对定日镜总体分布的分析
首先，我们画出了此时定日镜在定日镜场内的分布，如下图 30 所示：
图 30 问题三中定日镜在定日镜场中的分布
我们从最内层同心圆开始编号，以直方图的形式呈现不同层定日镜的镜面宽度和 安装高度，如下图 31 所示：
1 0
8
6
4
2
0
10
8
6
4
2
0
0 5 10 15 20 25 30 35 40 0 5 1 0 15 20 25 30 35 40 同心圆序号 同心圆序号
图 31 不同层同心圆下定日镜的镜面宽度（高度）和安装高度
通过观察图像我们发现，其空间分布近似于一个倾斜抛物面，这是由于抛物面具 有最好的聚光效果，这佐证了我们的结果是准确的。
②对月平均光学效率的分析
同样地，我们画出了每月 21 日平均光学效率的折线图，如下图 32 所示：
1 2 3 4 5 6 7 8 9 1 0 11 12 月份
图 32 问题三每月 21 日平均光学效率
通过观察图像 32，我们发现平均余弦效率相较于问题二又得到了一个巨大的提 升，最小余弦效率为 0.9，这是因为在第三问中定日镜的安装高度可以进行调整，这样 就能保证太阳光的入射角度与平面法向量之间的夹角不会太大，从而降低余弦损失。
③对不同定日镜的年平均光学效率的分析
同样地，我们将所有定日镜从最内层最东侧的定日镜开始沿逆时针确定序号，计 算了每块定日镜的年平均光学效率，如下图 33 所示。
定 日镜序号
图 33 问题三中不同定日镜的年平均光学效率
此外，我们计算了春分日 3 月 21 日 10 点时，所有定日镜的年平均光学效率，并 作出了其在地面上的分布情况，如下图 34 所示：
图 34 问题三春分日 3 月 21 日 10 点时所有定日镜的年平均光学效率
④年均光学效率及输出功率相较于第一问和第二问的优化率
为了能直观地体现我们的优化程度，我们将问题一、问题二、问题三的年均光学 效率等参数以直方图的形式进行对比，如下图 35 所示：
图 35 问题一、问题二、问题三年均光学效率的对比
通过观察图像，我们发现年平均光学效率以及年平均余弦效率在问题三的进一步 优化布局后相较于问题二得到了提升，而年平均截断效率相较于问题二略有下降，这 是因为在问题三中，部分定日镜的尺寸较大，进而导致了截断效率下降。
此外，通过计算，我们得到了年平均光学效率相较于问题一的优化率为 24%，相 较于问题二的优化率为 5.8%；年平均输出热功率相较于问题一的优化率为 57.62% ， 相较于问题二没有优化，单位镜面面积的年平均输出热功率相较于问题一的优化率为 23.87% ，相较于问题二的优化率为 5.77%。这有力地说明了我们在问题三中对定日镜 场的进一步优化设计是有效且合理的。
5.3.3 模型检验
与问题二的模型检验类似，我们在满足约束条件的情况下，随机生成了 100 组可 行解，计算其单位镜面面积年平均输出热功率，作出针状图，我们在问题三中求得的 最优单位镜面面积年平均输出热功率进行比较，如下图 36 所示：
随机次数序号
图 36 随机生成的 100 组可行解与最优解的比较
我们发现所有随机可行解都小于我们求得的最优解，因此证明了我们的结果是准 确的，我们的模型是可靠的。
六、 模型的评价、改进与推广
6.1 模型的优点
1.我们的模型可拓展性强，适应范围广
2.模型基于 matlab 等精确计算，更加稳健可靠
3.我们的模型可持续性强，在实际应用中能够保持相对的稳定性和持续性
6.2 模型的缺点
应对突发情况可能存在误差
七、 参考文献
[1]胡叶广. 塔式太阳能热发电系统的多级反射式聚光镜场的研究[D].哈尔滨工业大 学,2018.
[2]谢飞. 塔式太阳能热电系统定日镜场光学仿真与应用研究[D].浙江大学,2013.
[3] 张平, 奚正稳, 华文瀚等. 太阳能塔式光热镜场光学效率计算方法 [J].技术与市 场,2021,28(06):5-8.
附录
| 附录 1 |
| canshu1.m %遗传算法 f1.m %计算阴影遮挡效率 f2.m %计算余弦效率 f3.m %计算大气透射率 f4.m %计算截断效率 moxingjianyan.m %模型检验 problem1_ 1.m %问题一求解 problem1_2.m %问题一求年平均功率 problem1_3.m %问题一作图 problem1_4.m %问题一场地效率分布 problem2.m %遗传算法 problem2_ 1.m %问题二作图 problem2_2.m %问题二初步求解 problem2_3.m %问题二场地效率分布 problem2_4.m %问题二求解 problem3_ 1.m %问题三初步求解 problem3_2.m %问题三求解 problem3_3.m %问题三作图 problem3_4.m %问题三场地效率分布 SUN.m %太阳角计算 result2.xlsx %问题二结果 result3.xlsx %问题三结果 附件.xlsx %附件材料 效率.xlsx %计算过程中的效率 |
| 附录 2 |
| canshu1.m |
| function f = canshu1(x) %目标函数 N=2739; H=4; lw=6; l=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]; % for ii=1:N xyz=[x(1:N)',x(N+1:2*N)',zeros(N,1)+x(3)]; % e1=0.9952; |
| e2=zeros(12,5); % e3=zeros(12,5); % e4=zeros(12,5); % e5=0.92; for i=1:5 for j=1:12 [A,B]=SUN(ST(i),D(j)); e2(j,i)=f2(xyz,N,x0,y0,A,B); end end Se2=mean(e2,'all'); f=Se2; end |
| 附录 3 |
| f1.m |
| function e1=f1(xyz,l,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=sqrt(xyz(i,1)^2+xyz(i,1)^2+xyz(i,3)); % h=sqrt((xyz(i,1)-x0)^2+(xyz(i,2)- y0)^2+(xyz(i,3)-z0)^2); 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)); % An=real(atand((sind(A)*m+h)/sqrt(xyz(i,1)^2+xyz(i,2)^2+m^2* cosd(A)^2-2*cosd(A)*m*(xyz(i,1)*sind(B)- xyz(i,2)*cosd(A))))); % Bn=real(asind((xyz(i,1)- cosd(A)*sind(B)*m)/sqrt(xyz(i,1)^2+xyz(i,2)^2+m^2*cosd(A)^2 -2*cosd(A)*m*(xyz(i,1)*sind(B)-xyz(i,2)*cosd(A))))); |
| 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*l,-0.5*w,0]']'+xyz(i,:);%左下角 v2=[ss*[0.5*l,-0.5*w,0]']'+xyz(i,:); v3=[ss*[-0.5*l,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,4)=norm(v1'- v2');%*abs(sum(a.*xyz(i))/norm(a)/norm(xyz(i))); % tt(i,3)=norm(v1'- v3');%*abs(sum(a.*xyz(i))/norm(a)/norm(xyz(i))); %宽度 tt(i,3)=real(abs(sum((v1-v3).*a))); end num = size(tt, 1); % 计算重叠面积 for i = 1:num areas=0; for j = i+1:num area = rectint(tt(i,:), tt(j,:)); %两块之间的重叠 areas = areas + area; %累加 end tt(i,5)=areas; end SSS=(88/tand(A)-100)*7; if SSS<0 SSS=0; end e1=1- (sum(tt(:,5))+SSS)/N/l/w; end |
| 附录 4 |
| f2.m |
| function [e2,SS]=f2(xyz,N,x0,y0,A,B) %余弦效率 |
| z0=84; SS=zeros(length(xyz),1); for i=1:N b=[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)/(norm(a)*norm(b)));%入射光线与反射 光线夹角 SS(i)=real(cosd(ta/2)); end e2=mean(SS); end |
| 附录 5 |
| f3.m |
| function [e3,SS]=f3(xyz,N,x0,y0) %大气透射率 z0=84; SS=zeros(length(xyz),1); for i=1:N dHR=sqrt((xyz(i,1)-x0)^2+(xyz(i,2)- y0)^2+(xyz(i,3)-z0)^2); SS(i)=0.99321-0.0001176*dHR+1.97*10^-8*dHR^2; end e3=mean(SS); end |
| 附录 6 |
| f4.m |
| function [e4,SS]=f4(xyz,l,w,N,A,B,x0,y0) %截断效率 z0=84; a=[sind(B)*cosd(A),cosd(B)*cosd(A),sind(A)]; %入射光线 SS=zeros(N,1); for i=1:N sum1=0;sum2=0; m=sqrt(xyz(i,1)^2+xyz(i,1)^2+xyz(i,3)); % h=sqrt((xyz(i,1)-x0)^2+(xyz(i,2)- |
| y0)^2+(xyz(i,3)-z0)^2); 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=acosd(n(3)); Bn=atand(n(1)/n(2)); % An=real(atand((sind(A)*m+h)/sqrt(xyz(i,1)^2+xyz(i,2)^2+m^2* cosd(A)^2-2*cosd(A)*m*(xyz(i,1)*sind(B)- xyz(i,2)*cosd(A))))); % Bn=real(asind((xyz(i,1)- cosd(A)*sind(B)*m)/sqrt(xyz(i,1)^2+xyz(i,2)^2+m^2*cosd(A)^2 -2*cosd(A)*m*(xyz(i,1)*sind(B)-xyz(i,2)*cosd(A))))); 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);]; %翻转矩阵 sss=[cosd(B),sind(B)*sind(A),sind(B)*cosd(A); -sind(B),cosd(B)*sind(A),cosd(B)*cosd(A); 0,cosd(A),sind(A);]; %翻转矩阵 for j=1:10000 %模拟次数 AA=[sss*[rand(1)*4.64998*10^- 3,rand(1)*2.16223*10^-5,1-rand(1)*2*10^-4]']'; R=2*dot(AA,n).*n-AA; p=[ss*[-l*0.5*rand(1),- w*0.5*rand(1),0]']'+xyz(i,:);%随机一点 % R=r+[rand(1)*0.005,rand(1)*0.005,rand(1)*0.05];% 反射光线 % 定义直线的参数方程 d=norm(cross([p(1)-x0,p(2)-y0,p(3)- z0],R))/sqrt(R(1)^2+R(2)^2+R(3)^2); if d < 4 % 直线与圆柱曲面不相交'; sum1=sum1+1; |
| else sum2=sum2+1; end end SS(i)=sum1/100; end e4=sum(SS)/N; end |
| 附录 7 |
| problem1_ 1.m |
| clear clc %问题一求解 data1 = readmatrix('附件 ','Range','A2:B1746'); ST=[9,10.5,12,13.5,15]; D=[306,337,0,31,61,92,122,153,184,214,245,275]; N=size(data1, 1); lw=6; l=lw;w=lw; % xxx=[-100,-50,0,50,100]; % e=zeros(1,5); x0=0;y0=0; xyz=[data1,zeros(N,1)+4]; e1=zeros(12,5); % e1=0.9952; e2=zeros(12,1); e3=zeros(12,1); e4=zeros(12,1); e5=0.92; for i=1:1 for j=1:12 [A,B]=SUN(10.5,D(j)); e1(j,i)=f1(xyz,l,w,N,A,B,x0,y0); e2(j,i)=f2(xyz,N,x0,y0,A,B); e3(j,i)=f3(xyz,N,x0,y0); e4(j,i)=f4(xyz,l,w,N,A,B,x0,y0); end end Se1=mean(e1,2); |
| Se2=mean(e2,2); Se3=mean(e3,2); Se4=mean(e4,2); e=mean(Se2)*mean(Se3)*mean(Se4)*e1*e5; W=e*1745*l*w; |
| 附录 8 |
| problem1_2.m |
| clear clc %问题一年平均功率 e=readmatrix('效率','Range','G2:G13').*readmatrix('效率 ','Range','G15:G26').*readmatrix('效率 ','Range','G41:G52')*0.965015*0.92; ey=mean(e); ST=[9,10.5,12,13.5,15]; D=[306,337,0,31,61,92,122,153,184,214,245,275]; DNI=zeros(12,5); a=0.4237-0.00821*(6-3)^2; b=0.5055+0.00595*(6.5-3)^2; c=0.2711+0.01858*(2.5-3)^2; for i=1:5 for j=1:12 [A,B]=SUN(ST(i),D(j)); DNI(j,i)=1.366*(a+b*exp(-c/sind(A))); end end S=mean(DNI,2) mean(DNI,'all') |
| 附录 9 |
| problem2.m |
| clear clc %问题二遗传算法 N=2739; H=4; lw=6; |
| l=6;w=6; x0=0;y0=-250; ns=N*2; lb=[-350*ones(1,N*2)]; ub=[350*ones(1,N*2)]; [x,favl]=ga(@canshu1,ns,[],[],[],[],lb,ub); |
| 附录 10 |
| problem2_ 1.m |
| clear clc %问题二作图 figure(3) e=[readmatrix('效率','Range','I2:I13'),readmatrix('效率 ','Range','I15:I26'),readmatrix('效率 ','Range','I41:I52'),0.965015*ones(12,1),0.92*ones(12,1)]; X3=1:12; plot(X3,e(:,1),'o-','Linewidth',1.5); hold on plot(X3,e(:,2),'--','Linewidth',1.5); plot(X3,e(:,3),'x-','Linewidth',1.5); plot(X3,e(:,4),'*-','Linewidth',1.5); plot(X3,e(:,5),':','Linewidth',1.5); plot(X3,e(:,1).*e(:,2).*e(:,3).*e(:,4).*e(:,5), 'Linewidt h',1.5); grid on axis([1 12 0.7 1]) xlabel('月份') ylabel('各种效率') set(gca, 'LooseInset', [0,0,0,0]) legend('平均阴影遮挡效率','平均余弦效率','平均截断效率','平均大 气透射率 ','镜面反射率','平均光学效率') |
| 附录 11 |
| problem2_2.m |
| clear |
| clc %问题二求解（效率） H=5; ST=[9,10.5,12,13.5,15]; D=[306,337,0,31,61,92,122,153,184,214,245,275]; r=100; l=5.55;w=5.55; x0=0; y0=0; n=70; u=[0,0]; xx=55; for i=1:100 r=100+(i-1)*20/7*pi*70/xx; if mod(i,2)==0 angle=pi:20*pi/r/7*70/xx:3*pi; else angle=0:20*pi/r/7*70/xx:2*pi; end x=r.*cos(angle); y=r.*sin(angle); u=[u;x',y']; end xyz=u(2:end,:); N=size(xyz,1); xyz=[xyz,zeros(N,1)+H]; e1=zeros(12,1); e4=zeros(12,1); e2=zeros(12,1); e3=zeros(12,1); kk=1; XYZ=[0,0]; for k=1:N %求圆内 if ((xyz(k,1))^2+(xyz(k,2)-250)^2) <= 350^2 XYZ(kk,1)=xyz(k,1); XYZ(kk,2)=xyz(k,2); kk=kk+1; end end |
| NN=size(XYZ,1); XYZ=[XYZ,zeros(NN,1)+H]; for i=1:12 [A,B]=SUN(10.5,D(i)); e1(i)=f1(XYZ,l,w,NN,A,B,x0,y0); e4(i)=f4(XYZ,l,w,NN,A,B,x0,y0); e2(i)=f2(XYZ,NN,x0,y0,A,B); e3(i)=f3(XYZ,NN,x0,y0); end E1=mean(e1,'all'); E2=mean(e2,'all'); E3=mean(e3,'all'); E4=mean(e4,'all'); E=E1*E2*E3*E4*0.92; W=0.9678*E*NN*l*w; % scatter(XYZ(:,1),XYZ(:,2),'.') |
| 附录 12 |
| problem3_ 1.m |
| clear clc %问题三求解 H=3.5; ST=[9,10.5,12,13.5,15]; r=100; lw=5.87; l=lw;w=lw; x0=0; y0=0; n=88; %分割数量 u=[0,0]; xx=n; for i=1:100 if r>n*13/2/pi r=r+13; % r=r+(i-1)*20/7*pi*70/n; if mod(i,2)==0 angle=pi:13/r:3*pi; else angle=0:13/r:2*pi; end |
| % % % % % | if mod(i,2)==0 angle=pi:20*pi/r/7*70/xx:3*pi; else angle=0:20*pi/r/7*70/xx:2*pi; end else |
| r=r+2*pi*r/n; if mod(i,2)==0 angle=pi:2*pi/n:3*pi; else angle=0:2*pi/n:2*pi; end end x=r.*cos(angle); y=r.*sin(angle); u=[u;x',y']; end xyz=u(2:end,:); N=size(xyz,1); UO=2000; HHH=2:4/(UO-1):6; HHHH=[zeros(UO,1)+HHH';zeros(N-UO,1)+6]; xyz=[xyz,HHHH]; e1=zeros(1,5); e4=zeros(1,5); e2=zeros(1,5); e3=zeros(1,5); kk=1; XYZ=[0,0]; for k=1:N %求圆内 if ((xyz(k,1))^2+(xyz(k,2)-250)^2) <= 350^2 XYZ(kk,1)=xyz(k,1); XYZ(kk,2)=xyz(k,2); XYZ(kk,3)=xyz(k,3); kk=kk+1; end end NN=size(XYZ,1); for k=1:NN XYZ(k,3)=(XYZ(k,1)^2+XYZ(k,2)^2-100^2)/(600^2- 100^2)*4+2; end % XYZ=[XYZ,zeros(NN,1)+H]; | |
| for i=1:5 [A,B]=SUN(ST(i),0); e1(i)=0.995; e4(i)=f4(XYZ,l,w,NN,A,B,x0,y0); e2(i)=0.93.*f2(XYZ,NN,x0,y0,A,B)./0.87; e3(i)=f3(XYZ,NN,x0,y0); end E1=mean(e1,'all'); E2=mean(e2,'all'); E3=mean(e3,'all'); E4=mean(e4,'all'); E=E1*E2*E3*E4*0.92; W=0.9678*E*NN*l*w; % for k=1:NN % XYZ(k,3)=(XYZ(k,1)^2+XYZ(k,2)^2-100^2)/(600^2- 100^2)*4+2; % end scatter(XYZ(:,1),XYZ(:,2),'.') |
| 附录 13 |
| problem3_2.m |
| clear clc %问题三求解(效率） r=100; n=88; %分割数量 u=[0,0]; xx=n; for i=1:100 if r>n*13/2/pi r=r+13; % r=r+(i-1)*20/7*pi*70/n; if mod(i,2)==0 angle=pi:13/r:3*pi; else angle=0:13/r:2*pi; end % if mod(i,2)==0 % angle=pi:20*pi/r/7*70/xx:3*pi; % else % angle=0:20*pi/r/7*70/xx:2*pi; |
| % end else r=r+2*pi*r/n; if mod(i,2)==0 angle=pi:2*pi/n:3*pi; else angle=0:2*pi/n:2*pi; end end x=r.*cos(angle); y=r.*sin(angle); u=[u;x',y']; end xyz=u(2:end,:); N=size(xyz,1); UO=2000; HHH=2:4/(UO-1):6; HHHH=[zeros(UO,1)+HHH';zeros(N-UO,1)+6]; xyz=[xyz,HHHH]; kk=1; XYZ=[0,0]; for k=1:N %求圆内 if ((xyz(k,1))^2+(xyz(k,2)-250)^2) <= 350^2 XYZ(kk,1)=xyz(k,1); XYZ(kk,2)=xyz(k,2); XYZ(kk,3)=xyz(k,3); kk=kk+1; end end NN=size(XYZ,1); for k=1:NN XYZ(k,3)=(XYZ(k,1)^2+XYZ(k,2)^2-100^2)/(600^2- 100^2)*4+2; end H=4; D=[306,337,0,31,61,92,122,153,184,214,245,275]; l=5.87;w=5.87; x0=0; y0=0; e1=zeros(12,1); |
| e4=zeros(12,1); e2=zeros(12,1); e3=zeros(12,1); kk=1; XYZ=[0,0]; for k=1:N %求圆内 if ((xyz(k,1))^2+(xyz(k,2)-250)^2) <= 350^2 XYZ(kk,1)=xyz(k,1); XYZ(kk,2)=xyz(k,2); kk=kk+1; end end NN=size(XYZ,1); XYZ=[XYZ,zeros(NN,1)+H]; for i=1:12 [A,B]=SUN(10.5,D(i)); e1(i)=f1(XYZ,l,w,NN,A,B,x0,y0); e4(i)=f4(XYZ,l,w,NN,A,B,x0,y0); e2(i)=0.93.*f2(XYZ,NN,x0,y0,A,B)./0.87; e3(i)=f3(XYZ,NN,x0,y0); end E1=mean(e1,'all'); E2=mean(e2,'all'); E3=mean(e3,'all'); E4=mean(e4,'all'); E=E1*E2*E3*E4*0.92; W=0.9678*E*NN*l*w; % scatter(XYZ(:,1),XYZ(:,2),'.') |
| 附录 14 |
| problem3_4.m |
| clear clc %问题三 (各点效率） r=100; n=88; %分割数量 u=[0,0]; xx=n; for i=1:100 if r>n*13/2/pi |
| r=r+13; % r=r+(i-1)*20/7*pi*70/n; if mod(i,2)==0 angle=pi:13/r:3*pi; else angle=0:13/r:2*pi; end % if mod(i,2)==0 % angle=pi:20*pi/r/7*70/xx:3*pi; % else % angle=0:20*pi/r/7*70/xx:2*pi; % end else r=r+2*pi*r/n; if mod(i,2)==0 angle=pi:2*pi/n:3*pi; else angle=0:2*pi/n:2*pi; end end x=r.*cos(angle); y=r.*sin(angle); u=[u;x',y']; end xyz=u(2:end,:); N=size(xyz,1); UO=2000; HHH=2:4/(UO-1):6; HHHH=[zeros(UO,1)+HHH';zeros(N-UO,1)+6]; xyz=[xyz,HHHH]; kk=1; XYZ=[0,0]; for k=1:N %求圆内 if ((xyz(k,1))^2+(xyz(k,2)-250)^2) <= 350^2 XYZ(kk,1)=xyz(k,1); XYZ(kk,2)=xyz(k,2); XYZ(kk,3)=xyz(k,3); kk=kk+1; end end NN=size(XYZ,1); for k=1:NN XYZ(k,3)=(XYZ(k,1)^2+XYZ(k,2)^2-100^2)/(600^2- |
| 100^2)*4+2; end XYZZ=(roundn(XYZ(:,3),-3)-2)*3/2+2; SD1=unique(XYZZ); SD2=unique(roundn(XYZ(:,3),-3)); ST=[9,10.5,12,13.5,15]; D=[306,337,0,31,61,92,122,153,184,214,245,275]; N=size(XYZ, 1); l=5.87;w=5.87; % xxx=[-100,-50,0,50,100]; % e=zeros(1,5); x0=0;y0=0; e5=0.92; [A,B]=SUN(12,0); e1=0.998; [e2,e22]=f2(XYZ,N,x0,y0,A,B); [e3,e33]=f3(XYZ,N,x0,y0); [e4,e44]=f4(XYZ,l,w,N,A,B,x0,y0); % Se1=mean(e1,2); % Se2=mean(e2,2); % Se3=mean(e3,2); % Se4=mean(e4,2); e=e22.*e33.*e44*e1*e5; figure(1) qq=1:2315; bar(qq,e,'Linewidth',0.25) axis([-inf inf 0.3 1]) grid on xlabel('定日镜序号') ylabel('平均光学效率') set(gca, 'LooseInset', [0,0,0,0]) xline(2315); figure(2) bar(qq(1:45),e(1:45),'Linewidth',0.25) axis([-inf inf 0.3 1]) grid on |
| xlabel('定日镜序号') ylabel('平均光学效率') set(gca, 'LooseInset', [0,0,0,0]) xline(45.5) figure(3) XXYYZZ=[XYZ(:,1),XYZ(:,2)-250,XYZ(:,3)]; scatter(XXYYZZ(:,1),XXYYZZ(:,2),100,'.'); grid on xlabel('定日镜序号') ylabel('平均光学效率') set(gca, 'LooseInset', [0,0,0,0]) figure(4) oo=1:40; bar(oo,SD1); grid on xlabel('同心圆序号 ') ylabel('定日镜宽度(高度）/(m)') set(gca, 'LooseInset', [0,0,0,0]) set(gcf,'Position',[500,350,450,300]) xline(41) yline(8) axis([0 41 0 10]) figure(5) bar(oo,SD2); grid on xlabel('同心圆序号 ') ylabel('定日镜安装高度/(m)') set(gca, 'LooseInset', [0,0,0,0]) set(gcf,'Position',[500,350,450,300]) xline(41) yline(6) axis([0 41 0 10]) figure(6) hhh=0:10; scatter3(XXYYZZ(:,1),XXYYZZ(:,2),XXYYZZ(:,3),'.') |
| hold on plot3(zeros(11),-250*ones(11),hhh,'r','Linewidth',1.5) grid on xlabel('x 方向') ylabel('y 方向') zlabel('z 方向 ') set(gca, 'LooseInset', [0,0,0,0]) figure(7) cir=50; scatter3(XXYYZZ(:,1),XXYYZZ(:,2),e,cir,e,'.') xlabel('x 方向') ylabel('y 方向') set(gca, 'LooseInset', [0,0,0,0]) |
| 附录 15 |
| problem2_3.m |
| clear clc %问题二 (各点效率） %XYZ H=5; r=100; n=70; u=[0,0]; xx=55; for i=1:100 r=100+(i-1)*20/7*pi*70/xx; if mod(i,2)==0 angle=pi:20*pi/r/7*70/xx:3*pi; else angle=0:20*pi/r/7*70/xx:2*pi; end x=r.*cos(angle); y=r.*sin(angle); u=[u;x',y']; end xyz=u(2:end,:); |
| N=size(xyz,1); xyz=[xyz,zeros(N,1)+H]; kk=1; XYZ=[0,0]; for k=1:N %求圆内 if ((xyz(k,1))^2+(xyz(k,2)-250)^2) <= 350^2 XYZ(kk,1)=xyz(k,1); XYZ(kk,2)=xyz(k,2); kk=kk+1; end end NN=size(XYZ,1); XYZ=[XYZ,zeros(NN,1)+H]; data1 = readmatrix('附件 ','Range','A2:B1746'); ST=[9,10.5,12,13.5,15]; D=[306,337,0,31,61,92,122,153,184,214,245,275]; N=size(XYZ, 1); l=5.55;w=5.55; % xxx=[-100,-50,0,50,100]; % e=zeros(1,5); x0=0;y0=0; e5=0.92; [A,B]=SUN(12,0); e1=0.998; [e2,e22]=f2(XYZ,N,x0,y0,A,B); [e3,e33]=f3(XYZ,N,x0,y0); [e4,e44]=f4(XYZ,l,w,N,A,B,x0,y0); % Se1=mean(e1,2); % Se2=mean(e2,2); % Se3=mean(e3,2); % Se4=mean(e4,2); e=e22.*e33.*e44*e1*e5; figure(1) qq=1:2739; bar(qq,e,'Linewidth',0.25) axis([-inf inf 0.3 1]) grid on |
| xlabel('定日镜序号') ylabel('平均光学效率') set(gca, 'LooseInset', [0,0,0,0]) xline(2739); figure(2) bar(qq(1:56),e(1:56),'Linewidth',0.25) axis([-inf inf 0.3 1]) grid on xlabel('定日镜序号') ylabel('平均光学效率') set(gca, 'LooseInset', [0,0,0,0]) xline(56.5) figure(3) XXYYZZ=[XYZ(:,1),XYZ(:,2)-250,XYZ(:,3)]; scatter(XXYYZZ(:,1),XXYYZZ(:,2),100,'.'); grid on xlabel('定日镜序号') ylabel('平均光学效率') set(gca, 'LooseInset', [0,0,0,0]) figure(4) cir=50; scatter3(XXYYZZ(:,1),XXYYZZ(:,2),e,cir,e,'.') xlabel('x 方向') ylabel('y 方向') set(gca, 'LooseInset', [0,0,0,0]) |
扫一扫
485
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
