简介2023年全国大学生数学建模竞赛A题“定日镜场的优化设计”国家一等奖配套源码与论文面向参赛选手、科研学习者可用于复现定日镜场建模与优化全过程。资源包约30MB主要包含Python脚本、完整论文及推导材料便于按步骤还原代码与公式。已有115人学习浏览适合数模备赛的进阶参考。内容覆盖两问问题一通过三维坐标系和矢量理论推导余弦效率结合定日镜方位角仰角求解阴影遮挡效率再利用三余弦定理导出反射椭圆光斑方程并积分计算截断效率问题二围绕额定功率60MW优化吸收塔位置、定日镜尺寸、安装高度、数量及布点使单位镜面面积年平均输出热功率最大。论文中对截断效率增速、春分与夏至光学效率分布等关键结论有专门分析代码模块化程度高能帮助读者从公式推导到程序落地完整掌握国赛一等奖做法。1. 定日镜场优化设计这道国赛A题的完整Python方案该怎么用定日镜场优化设计是2023年全国大学生数学建模竞赛A题的核心表面是算光学效率、比优化结果实际考的是三维几何建模、效率分解和带约束寻优这三块恰好是所有参赛队最容易翻车的地方。这份资源打包了国家一等奖论文和配套Python源码论文把余弦效率、阴影遮挡效率、截断效率的推导过程写得很完整代码把逐时刻仿真和镜场布局优化跑通了适合正在备赛想复现真题的人也适合做光热电站定日镜场建模、需要一套可改写的几何光学计算底座的工程师。下面按我拆解项目的顺序来讲模型怎么建、问题一代码怎么读、问题二优化怎么调、坑在哪些位置。2. 先把物理模型立住从坐标变换到四个效率2.1 为什么以吸热塔为原点建坐标系定日镜场建模的第一步不是写公式而是选坐标系。这个问题里最大的坐标难点是每一面镜子的位置都不同太阳位置又随时间变化如果坐标系选在镜场一角或者直接套用经纬度下的平面近似后面计算方位角、仰角、遮挡判断时会不断返工。论文最终以吸热塔为原点建立地面三维坐标系东向为x轴、北向为y轴、z轴垂直向上这等于给所有矢量找到了同一个参考基。把我自己的习惯说给你拿到原始数据后不要急着算效率先做一步坐标标准化。把每面镜子的x、y转成相对吸热塔的平面距离这个量后面算大气透射率时必用而且做完这步镜场落点图能直接画出来初步检查数据有没有异常。原题里定日镜的安装高度和镜面尺寸是给定的但对坐标来说真正的关键点只有两个镜子中心的三维坐标和吸热塔窗口中心的三维坐标。前者决定反射光起点后者决定反射光终点两者的矢量差就是反射方向。2.2 余弦效率入射光与镜面法向的夹角投影余弦效率是所有效率里最基础的它等于cos(θ)其中θ是太阳入射方向与镜面法向的夹角。定日镜是双轴跟踪的镜面始终转动保证反射光打向吸热窗口所以这个夹角不可能直接读出来必须一步步算。论文用的是欧氏几何的矢量理论镜面法向是入射光方向和反射光方向的角平分线把这两个单位矢量相加、再归一化就得到法向单位矢量余弦效率就是入射单位矢量与法向单位矢量的点积。实际写代码时这里有个容易搞混的细节入射光方向要从镜面指向太阳反射光方向要从镜面指向吸热塔窗口两个都是“从镜心出发”的矢量。很多初学者会把反射光方向跟塔位置矢量直接绑定忘记减掉镜心坐标导致所有镜子共用同一个反射方向。我会在代码里把incident、reflect、normal三个变量分开命名每个占总功率的贡献也单独打印跑单面镜子时先手算验证一遍后面排查问题会轻松很多。2.3 阴影遮挡效率镜面面积占比的几何求解阴影遮挡效率关注的是两块内容前排名次与镜面反射关系的累加遮挡以及后排镜面对反射光线的拦截。论文的思路是先求出定日镜的方位角和仰角再把镜面投影到入射光方向和反射光方向统计被相邻镜面遮挡的面积占自身面积的比重得到阴影遮挡效率。这块的常见误区是把镜面当质点处理那样阴影遮挡效率永远只能是0或1。竞赛题里的镜面是正方形或近似正方形被遮挡时参与计算的是一个多边形面积占比得靠多边形求交。我在复现时采用的做法是把被测镜面的四个顶点沿入射光方向做平行投影再判断这些投影点有多少落入相邻镜面的投影多边形中用落入区域的重叠面积除以镜面总面积。这个方法的优点是编程量小缺点是边界处误差较大所以采样点要加密每边取10到20个点结果就足够稳定。2.4 截断效率椭圆光斑的能量密度积分与三余弦定理截断效率是四个效率里最考验几何功底的。定日镜反射到吸热塔窗口的光斑是椭圆窗口是矩形或者圆形开口光斑落在窗口范围内的能量比例就是截断效率。论文的推导路径是利用三余弦定理求出反射光与窗口平面的交线得到椭圆光斑区域的方程再在窗口区域内做能量密度积分积分得到的能量除以光斑总能量得到截断效率。我复现时的简化做法是先在以窗口中心为原点的局部二维坐标系里写出椭圆方程然后把问题转化为求矩形窗口与椭圆的相交面积。注意太阳本身有约4.65毫弧度的角半径镜面面形也不是理想平面能量密度并非均匀但竞赛题的附件公式通常默认均匀分布所以实际计算时保留均匀密度假设重点是把椭圆参数求准。这个椭圆从哪里来——把镜面轮廓沿反射方向投影到窗口平面得到的就是椭圆边界三余弦定理在这里的作用是精确计算投影变换。2.5 大气透射率与光学效率的工程近似大气透射率由附件里的经验公式给出核心变量是镜面到吸热塔的平面距离距离越远光束在大气中走过的路径越长透射率越低。光学效率则是余弦效率、阴影遮挡效率、截断效率、大气透射率四项连乘有的赛题版本还会再乘镜面反射率。把四项效率分开写而不是直接写一个总效率是为了后面做优化时定位瓶颈如果某时刻整体效率偏低能直接看出是截断效率垮了还是余弦效率不行。一个值得注意的细节是附件公式里的大气透射率用的是平面距离还是空间距离要看题目的原始定义。我拆这份代码时发现论文里用的是镜面中心到吸热塔窗口中心的水平投影距离这个参数在2.1里已经提前算好所以衔接很顺。你在改自己的数据时距离的定义一定要跟附件公式保持一致否则透射率会整体偏移。3. 问题一落地把四个效率翻译成逐时刻的Python代码3.1 输入数据与关键参数初始化问题一的输入是给定的吸热塔位置、定日镜尺寸和安装高度需要产出年平均光学效率、年平均输出热功率和单位镜面面积年平均输出热功率。代码的第一步是读入镜场数据并初始化参数我建议把塔高、窗口尺寸、纬度、经度都单独拎出来做成配置项方便后续换场景。import numpy as np import pandas as pd # 读入定日镜参数x坐标, y坐标, 镜面边长, 安装高度 df pd.read_csv(mirror_field.csv) mirrors df[[x, y, side, height]].values # 吸热塔与窗口参数 tower_height 80.0 # 吸热塔高度, m receiver_width 6.0 # 窗口宽度, m receiver_height 8.0 # 窗口高度, m lat 40.5 # 纬度, deg春分附近用典型值 lon 117.0 # 经度, deg # 单次计算的代表时刻夏至日 12:30 太阳时 day_of_year 172 solar_time 12.5 hour_angle 15.0 * (solar_time - 12.0)这里把纬度、经度、塔高作为常量放在最前面是因为后面优化问题二时要反复改这些值。读入的镜面数据不做过多的格式转换直接用数组切片x和y是绝对坐标相对位置在计算反射矢量时再减掉塔坐标避免数据预处理阶段引入错误。3.2 太阳位置从赤纬到时角到高度角方位角太阳位置是定日镜场仿真的时间基准。竞赛题一般允许使用简化天文公式太阳赤纬用Cooper方程近似太阳时角用12点减当地太阳时换算。注意这里默认用的是太阳时而不是北京时间两者之间的时差如果忽略方位角会偏移截断效率也会跟着偏。# 太阳赤纬近似Cooper方程输入年中第几天 declination 23.45 * np.sin(np.deg2rad(360.0 * (284 day_of_year) / 365.0)) def sun_position(lat, dec, hour_angle): 返回太阳高度角和方位角(度)方位角从北向东为正 lat_r np.deg2rad(lat) dec_r np.deg2rad(dec) ha_r np.deg2rad(hour_angle) sin_e np.sin(lat_r) * np.sin(dec_r) \ np.cos(lat_r) * np.cos(dec_r) * np.cos(ha_r) sin_e np.clip(sin_e, -1.0, 1.0) elevation np.arcsin(sin_e) cos_az (np.sin(dec_r) - np.sin(elevation) * np.sin(lat_r)) / \ (np.cos(elevation) * np.cos(lat_r) 1e-12) cos_az np.clip(cos_az, -1.0, 1.0) azimuth np.rad2deg(np.arccos(cos_az)) # 上午时角为负方位角取负下午取正 if hour_angle 0: azimuth -azimuth return np.rad2deg(elevation), azimuth这段函数里最关键的是方位角的符号修正。arccos的值域只有0到180度对应北半球正午前后太阳都在南侧所以arccos给的是“偏离正南或正北”的锐角量必须根据时角正负区分上午和下午。我在第一次复现时漏了符号判断结果上午和下午的镜面方位角完全对称镜像全年效率曲线出现周期性抖动排查了很久才发现是这里的问题。3.3 反射矢量与法向镜面追踪的核心逻辑拿到太阳位置后下一步要构造三个矢量入射光方向、反射光方向、镜面法向。入射光方向从镜心指向太阳反射光方向从镜心指向吸热窗口中心法向由两者相加归一化得到。这里塔向量必须是“窗口中心坐标减镜心坐标”窗口中心是塔顶下方若干米处的接收面中心。def mirror_normal(sun_elev, sun_azim, mirror_pos, tower_pos): 入射单位矢量: 从镜心指向太阳 反射单位矢量: 从镜心指向吸热窗口中心 法向: 入射反射 归一化 e np.deg2rad(sun_elev) a np.deg2rad(sun_azim) incident np.array([ np.cos(e) * np.sin(a), np.cos(e) * np.cos(a), np.sin(e) ]) reflect_vec tower_pos - mirror_pos reflect reflect_vec / (np.linalg.norm(reflect_vec) 1e-12) normal incident reflect normal normal / (np.linalg.norm(normal) 1e-12) return incident, reflect, normal入射矢量的x分量用sin(方位角)、y分量用cos(方位角)是因为这里的方位角定义为从正北起顺时针正东是90度正北是0度。如果不统一这个约定后面跟镜场坐标的方向对照会非常别扭。归一化时加1e-12是为了防止零向量除零这在镜子紧贴塔基时会发生虽然实际不会出现但写防御性代码能省掉很多NaN排查时间。3.4 四个效率的组装与8760小时循环单时刻的效率计算在拿到法向后就基本完成了。余弦效率直接取入射光与法向的点积阴影遮挡需要额外写投影函数截断效率要做椭圆窗口相交大气透射率用距离查公式。组装成整体后最外层循环处理全年8760个小时把每个小时的效率按小时权重累加得到年平均。def compute_hourly_efficiency(mirror_pos, tower_pos, sun_elev, sun_azim, side_length, receiver_w, receiver_h): inc, ref, n mirror_normal(sun_elev, sun_azim, mirror_pos, tower_pos) # 余弦效率 cos_eff max(0.0, np.dot(inc, n)) # 阴影遮挡效率这里简化为1完整版用多边形投影 shadow_eff 1.0 # 截断效率把镜面投影到窗口平面求椭圆与窗口交叠 trunc_eff compute_truncation_eff( mirror_pos, ref, n, side_length, receiver_w, receiver_h ) # 大气透射率经验公式distance为镜塔平面距离 distance np.linalg.norm(mirror_pos[:2] - tower_pos[:2]) atm_eff 0.9932 - 0.0001176 * distance 1.97e-8 * distance**2 optical_eff cos_eff * shadow_eff * trunc_eff * atm_eff * 0.92 return optical_eff # 全年8760小时循环主框架 total_energy 0.0 for doy in range(1, 366): for hour in range(24): sun_elev, sun_az sun_position(lat, declination, 15.0 * (hour - 12)) if sun_elev 0: continue # 夜间不做计算 # 累加光学效率 ...这里的0.92是镜面反射率的典型取值竞赛题中一般会给出具体值有的版本直接取0.9。注意夜间判断不能省略太阳高度角为负时虽然定日镜不会工作但这个时刻如果不剔除效率会被零值稀释年平均光学效率会明显偏低和论文结果对不上。3.5 输出项对齐三个结果一次算清问题一要交三个数年平均光学效率、年平均输出热功率、单位镜面面积年平均输出热功率。第二个数等于总镜面面积乘以年平均光学效率再乘以当地法向直射辐照度第三个数再除以总镜面面积。代码里建议把总镜面面积存成全局变量因为问题二优化时这个值会变单独存储能减少重复计算。我拆代码时发现论文里对DNI的处理很细它不是一个常数而是随太阳高度角变化的函数高度角越低大气路径越长DNI越小。附件里通常给出的是水平面直射辐照度与法向直射辐照度的换算关系这部分直接套公式就行但要注意单位统一别把W/m²和MW混在一起我在核对结果时差点被单位换算坑了一次。4. 问题二优化设计60MW约束下的镜场布局怎么调4.1 优化目标的数学翻译问题二的描述里给出额定功率60MW所有定日镜尺寸及安装高度相同要求设计吸热塔位置坐标、定日镜尺寸、安装高度、定日镜数目和镜面位置让单位镜面面积年平均输出热功率尽量大。写成数学语言就是目标函数是输出热功率除以总镜面面积约束条件是输出热功率不低于60MW。这里的60MW很容易理解错。它不是要求最终计算结果精确等于60MW而是要求你的设计方案能达到这个额定功率水平超过它是可以的但超出太多意味着镜面面积浪费单位面积效率反而不高。所以实际优化时会发现最优解往往卡在“刚好满足60MW”附近也就是输出功率勉强越过红线同时镜面面积尽量小。这个边界特性决定了适应度函数的设计方式功率不足要重罚功率有余要适度奖励但不能奖励太多导致镜面面积失控。4.2 用截断效率增速确定镜面尺寸与安装高度论文里有个很聪明的处理先不管镜子位置分布单独分析截断效率随镜面尺寸和安装高度的变化规律。截断效率的增速到底是什么意思——当镜面尺寸从4米增大到5米时截断效率下降的幅度是陡还是缓当安装高度从4米升到8米时截断效率提升是否明显。把这两个维度组合起来扫描找到增速放缓的拐点在那个位置选尺寸和高度称之为“增速拐点选型”。import numpy as np sizes np.linspace(4.0, 15.0, 60) # 镜面边长范围, m heights np.linspace(3.0, 12.0, 40) # 安装高度范围, m growth_matrix np.zeros((len(sizes), len(heights))) for i, s in enumerate(sizes): for j, h in enumerate(heights): # 以最外圈镜子为基准计算截断效率得到单位面积增益 gain evaluate_truncation_growth( mirror_sides, install_heighth, tower_height80.0, field_radius200.0 ) growth_matrix[i, j] gain best_idx np.unravel_index(np.argmax(growth_matrix), growth_matrix.shape) best_size sizes[best_idx[0]] best_height heights[best_idx[1]] print(f最优尺寸: {best_size:.2f}m, 最优安装高度: {best_height:.2f}m)这段扫描代码的本质是二维网格搜索。用网格搜索而不是解析求导是因为截断效率本身含有椭圆-矩形相交的数值计算求导代价大且容易出错。二维扫描的分辨率要配合后续优化精度先粗扫定下大致区间再局部加密扫描一轮能把拐点位置控制在0.1米以内。注意这里评估基准选的是最外圈镜子因为最外圈镜子离塔最远截断效率最差用它做基准能保证整个镜场都不掉链子。4.3 镜场布局优化的算法选型定了尺寸和高度之后剩下的设计变量是吸热塔位置、镜面数目、每面镜子的位置。这是一个典型的高维混合整数优化问题因为镜面数是整数、坐标是连续值。论文里提到基于春分和夏至的余弦效率与光学效率分布来确定后续参数这意味着优化时代步是从8760小时压缩到几个代表时刻。我复现时采用的优化思路是分三步走第一步在极坐标下按径向距离和环向角度离散化候选镜位半径方向按30到80米间隔环向按7.5度到15度间隔第二步对离散镜位集合做启发式搜索变量是塔坐标、镜面边长、安装高度和镜面数但镜面边长和高度已经在上一步锁定所以实际搜索的高维变量只剩塔坐标和镜面数第三步把粗搜索结果带回8760小时完整仿真精算误差超过1%就回退调整。# 伪代码离散镜位上的启发式搜索框架 def fitness_snapshot(ind): tower_x, tower_y, num_mirrors ind # 按环形布局生成镜面候选位置的稀疏子集 candidates generate_ring_candidates(tower_x, tower_y, num_mirrors) # 用春分夏至的代表时刻快速估算 p_avg, s_total fast_estimate(candidates, tower_x, tower_y) if p_avg 60000: return -1e6 # 功率红线硬惩罚 return p_avg / s_total # 目标函数单位面积功率 # 搜索算法参数 pop_size 200 max_iter 120 mutation_rate 0.15这里用了硬惩罚而不是软惩罚原因是功率低于60MW的方案直接淘汰不参与选择能避免遗传算法把大量代际花在无效解上。种群数量取200迭代120代是我在类似镜场规模下的经验值再往上增长度有限反而容易浪费计算时间。4.4 春分夏至代表时刻的选取细节用代表时刻代替全年8760小时核心是怎么选时刻。论文用的是春分和夏至因为春分时太阳赤纬接近0昼夜等长能代表春秋季的平均状态夏至时太阳高度角最高是夏季的极限工况冬至虽然光照最差但竞赛题的镜场规模下冬至往往无法满足功率约束所以不作为代表。冬半年的影响主要体现在总功率是否达标上如果春分夏至的代表时刻算下来能超过60MW则全年平均大概率也能满足。每个代表日取几个时刻也有讲究。早晨7点、9点、正午、下午15点、17点五个时刻足够覆盖一天的效率曲线再多取对精度提升很小但计算量线性增长。我一般会对这五个时刻做加权正午权重最高早晨和傍晚权重递减这样估算的年平均效率更接近8760小时完整计算的结果。5. 避坑与常见问题定日镜场代码最典型的五个翻车点5.1 余弦效率算出来大于1这是新手最容易遇到的现象。余弦效率是入射光与法向夹角的余弦值值域必须在0到1之间一旦打印结果超过1说明矢量没有归一化。根因通常是incident或normal这两个矢量里有一个没有做单位化或者normal incident reflect之后直接拿去点积忘了除模长。解决方法是把入射、反射、法向三个矢量都单独打印出来观察模长确认都为1后再做点积。我在代码里加了np.linalg.norm的归一化并且在正午正对太阳时加了一个断言如果余弦效率不等于1就报错能第一时间发现问题。5.2 太阳方位角跨时跳变现象是上午和下午的方位角出现180度的台阶跳变效率曲线在正午前后不连续。根因是arccos的值域限制它只能返回0到180度而标准的太阳方位角从正北起算范围是-180到180度。解决方式是判断时角的符号上午时角为负方位角取负值下午时角为正方位角取正值。这个符号修正必须在方位角计算的函数内部完成不能在全局循环里补否则那些正好在正午时刻的采样点会反复横跳。5.3 阴影遮挡效率全部等于0或1现象是阴影遮挡效率永远是极端的两个值没有任何中间状态。这说明投影计算没有真正发生原因多半是投影方向取错了比如把镜面直接投影到地面而不是沿入射光方向投影。还有一种是采样点太少镜面四个顶点恰好全部落入或全部没落入邻镜投影区域中间过渡被跳过。解决方法是加密采样点把镜面每条边细分到10段以上并且用一个简单的两个矩形重叠案例先手算验证确认投影方向确实沿入射光方向而不是沿z轴垂直向下。5.4 截断效率对镜面位置过于敏感现象是靠近塔的镜子截断效率接近99%最外圈的镜子截断效率只有20%平均值被拉得很低。根因是光斑半径随镜塔距离增长而窗口尺寸固定没有随之放大。解决办法是先检查镜塔距离与塔高的比值经验上这个比值超过1.5之后截断效率会快速衰减如果论文结果里效率曲线明显偏斜大概率是镜场半径设计过大。调整方向是缩小最外圈镜子的径向距离或者在满足面积要求的前提下增大安装高度。5.5 优化过程收敛缓慢甚至震荡现象是遗传算法跑到200代目标函数还在区间里来回波动。根因通常是三个变量取值范围太宽导致搜索空间稀疏种群数量不足50导致选择压力太小功率约束的惩罚权重配比不合理导致大量个体都卡在边界附近。解决方法是先用第4.2节的网格扫描锁定镜面尺寸和高度把连续变量范围缩窄种群数量提升到200以上把功率不足的个体直接淘汰而不是乘以惩罚系数这样能大幅减少无效代际。6. 从复现到验证三分钟确认结果可信再改自己的模型拿到源码后直接跑8760小时循环得出一个数字就敢写进论文这是最危险的做法。我的习惯是强制走三步验证每一步都能在几分钟内完成。第一步是单镜手算验证取正对塔的一排镜子中离塔最近的一面上用某个已知时刻做手算太阳高度角45度、方位角0度时入射光方向与法向夹角应该等于反射光方向与法向夹角余弦效率要跟手算结果一致这一步能同时检验矢量构造和坐标系方向。第二步是对称性检查镜场布局如果关于x轴对称那么对称位置的镜子在同一时刻的效率应该完全一致程序里打印效率并按镜号排序如果发现不对称的镜子说明坐标读入或索引顺序有错。第三步是参数扰动对光斑半径和大气透射率公式里的系数做正负10%的扰动观察最终年平均输出热功率的变化幅度如果变化超过5%说明结果对这个参数过于敏感论文结论的可信度就要打个问号。遇到三处以上波动不一致时我不会继续调参而是回到单面镜子的计算函数里逐行打断点。最常抓到的问题不是公式错而是单位或坐标约定在某个函数边界上没对齐。从那以后我每次拿到光学几何类的建模题无论镜场、反射面还是天线阵列都强制走一遍单点手算加对称性检查再加参数扰动这三步整个排查流程不会超过二十分钟但省下的反复试错时间通常是几小时甚至一天。希望这份源码的拆解过程能帮你少走一段弯路也希望你在复现时保住“先验证再信任”的习惯。本文还有配套的精品资源点击获取