简介这是一份面向数学建模参赛者与太阳能光热利用方向学习者的定日镜场优化设计完整方案基于2023年高教社杯全国大学生数学建模竞赛A题整理兼顾光学、热力学与工程约束帮助读者理解从太阳轨迹建模、镜面朝向求解到成本效益评估的整体建模路径。资源包共23个文件以docx正文、xml数据记录、png示意图和jpg图形为主便于对照文档查看模型结构、结果图表与原始工程文件压缩包大小约844KB便于快速下载部署。目前已有44人学习使用适合备赛演练与课外拓展开源模型。包内包含完整的定日镜场布局优化思路、Matlab仿真计算代码、目标函数与约束处理细节并对镜面形状、排列方式、能量收集效率等关键因素做了可视化呈现读者可借助这些材料快速复用算法框架按需调整参数后应用于同类聚光系统设计问题提升建模实操与工程分析能力。1. 定日镜场优化设计模型这道 A 题到底在算什么拿到“定日镜场优化设计模型”这个题目时很多人的第一反应是“这不就是算镜子反射吗”。真正动手才发现它是一道把光学效率、站址几何、离散布局优化和全年辐照积分全串起来的综合优化题。2023 年数模国赛 A 题就是典型代表给定塔高、集热器尺寸和镜面尺寸让你在有限地面区域里布置定日镜目标函数是单位镜面面积上的年平均输出热功率最大。适合正在备赛的队伍也适合刚接触塔式光热电站优化的工程师——这篇笔记会把从建模到 Matlab 实现的关键环节一次讲透。先说结论这道题的难点不在某个公式而在“优化变量怎么编码”和“效率函数怎么写才不回环震荡”。下面按我自己的解题顺序展开从几何建模、效率计算到优化器选型、代码实现最后把最容易翻车的几个坑单独列出来。2. 定日镜场建模从太阳位置到光学效率的四个乘子2.1 太阳位置与镜面法线先把方向向量统一到一套坐标系定日镜场里所有几何计算都围绕两个向量太阳光入射方向以及镜面中心指向集热器的反射方向。很多队伍一开始就在坐标正负号上栽跟头原因是太阳方位角的定义在不同教材里不一样。我的习惯是先把所有角度转换成“镜面中心指向太阳”的三维单位向量后续全部用向量点积避免反复翻公式。假定镜场平面是 x 轴向东、y 轴向北、z 轴向上太阳高度角为 alpha方位角 gamma 自正南向西为正。那么镜面中心指向太阳的单位向量可以写成% sunVector.m % 输入alpha 太阳高度角(度)gamma 太阳方位角(度正南为0向西为正) % 输出sunVec 三维单位向量方向为“镜面中心 - 太阳” function sunVec sunVector(alpha, gamma) alphaRad deg2rad(alpha); gammaRad deg2rad(gamma); % x 方向向东正南方向是 y 轴负方向 % gamma 向西为正所以投影到 x-y 平面时x 分量为负 sunVec [ -cos(alphaRad) * sin(gammaRad), ... -cos(alphaRad) * cos(gammaRad), ... sin(alphaRad) ]; % 归一化防止角度查表后微小误差 sunVec sunVec / norm(sunVec); end这里最容易被忽略的是符号。如果 gamma 按正南为 0、向东为正那分母上的负号位置要换。代码本身不重要重要的是整套模型里只能有一种约定而且后续所有阴影遮挡、法线计算都沿用同一个 sunVec。有了入射方向还需要镜面法线。定日镜是双轴追踪的镜面法线方向要保证反射光线正好指向集热器中心。反射方向 refl 是从镜面中心指向塔顶也就是towerTop - mirrorPos。镜面法线就是“太阳入射方向”和“反射方向”的角平分线方向% 镜面法线与余弦效率 refl towerTop - mirrorPos; refl refl / norm(refl); % 镜面 - 塔顶的单位向量 normal sunVec refl; % 角平分线方向 normal normal / norm(normal); eta_cos dot(sunVec, normal); % 余弦效率余弦效率是定日镜场里最基础的损失项。入射角为零时法线与太阳方向完全一致eta_cos 等于 1入射角越大镜面有效开口面积越小。实际镜场中远离塔的镜子余弦效率通常只有 0.85 到 0.95这部分损失靠镜面面积补不回来所以布局优化必须在“远离塔但遮挡少”和“靠近塔但余弦效率高”之间找平衡。2.2 四个效率乘子余弦效率只是第一层后面还有三个单面镜子反射到集热器的热功率工程上通常写成P DNI * A_mirror * eta_cos * eta_atm * eta_shadow * eta_truncDNI 是法向直射辐照度A_mirror 是镜面面积。eta_cos 刚才算过剩下三个效率分别对应不同物理过程缺一不可。大气透射效率 eta_atm 是镜面到集热器距离 d 的函数。空气会吸收和散射一部分反射能量距离越远损失越大。一个常见经验公式是eta_atm 0.99321 - 0.0001176 * d 1.97e-8 * d^2这里 d 的单位是米。注意这个公式在 d 超过 1000 米时可能算出不合理值实际镜场半径一般不会超过三四百米所以可以直接用。如果题目给了自己的公式优先用题目公式。截断效率 eta_trunc 描述的是反射光斑是否完全落在集热器窗口内。太阳光不是绝对平行光太阳半张角大约 4.65 mrad加上镜面面形误差和追踪误差光斑会随距离扩散。简化计算可以用r_spot d * tan(4.65e-3 epsilon)其中 epsilon 是镜面斜率误差通常取 1~3 mrad。集热器半径是 receiverR工程上常用高斯分布近似eta_trunc 1 - exp(-0.5 * (receiverR / r_spot)^2)这个公式的好处是连续可导适合优化器迭代。缺点是光斑远小于集热器时eta_trunc 会接近 1但不会超过 1不会出现逻辑错误。真正要小心的是 r_spot 接近零时除零所以代码里要对 r_spot 做下限保护。阴影遮挡效率 eta_shadow 是最难算的一项。它包含两层含义一是相邻定日镜在入射光方向上挡住当前镜面叫阴影二是当前镜面反射出去的光线被前排镜子挡住叫遮挡。精确做法是把所有镜面投影到与光线方向垂直的平面上用多边形求交面积。这个放在后面代码部分细说。3. 优化模型设计决策变量、目标函数与求解器选择3.1 决策变量编码与约束设计定日镜场优化不能真的把每一面镜子的 x、y 坐标都作为自由变量那样变量维度太高优化器根本搜不完。我见过的大多数获奖思路以及工程上常用的初始设计都采用同心圆环布局镜场由若干圈半径不同的环组成同一圈上的镜子按方位角均匀分布。决策变量可以编码为三个向量变量含义范围示例类型r_i第 i 圈镜场的半径60 ~ 300 m连续n_i第 i 圈的镜子数量6 ~ 40整数h_i第 i 圈镜面中心离地高度4 ~ 12 m连续如果某个圈不用就直接令 n_i 0。这样变量总数为 3 乘以最大圈数通常取 5 到 7 圈变量数量在 15 到 21 之间对遗传算法来说非常友好。约束条件至少要包含三条一是半径必须递增避免两个圈重叠二是相邻半径之差要大于镜面的横向尺寸否则遮挡率会高到无法接受三是所有镜面必须落在给定地面区域内。用 Matlab 的 ga 工具箱时线性约束可以写成 Aineq 矩阵半径递增这种约束也能写进去但更省事的办法是在目标函数里加惩罚项半径不满足要求时直接返回一个很大的正数。目标函数是单位镜面面积上的年平均热功率maximize P_avg / A_totalP_avg 的计算需要遍历全年代表性时刻。比赛题目通常会给每月 21 日 9:00 到 15:00 的太阳位置数据那么就对这 12 个月、若干个小时分别计算镜场瞬时功率再按每月天数加权平均。DNI 如果题目给了就用题目数据没给就按常数处理比如 800 W/m^2。优化问题本质上是求比值最大化所以 DNI 取常数不影响最优布局。3.2 算法选型ga 做全域搜索fmincon 做局部精调这个问题有整数变量有非线性目标函数而且目标函数里嵌着阴影遮挡计算噪声很大。fmincon 这类基于梯度的算法很容易陷进局部最优尤其是初始布局给得不好的时候。我的方案是先用遗传算法跑出大致区域再用模式搜索或 fmincon 做一轮精调前提是变量里没有整数。Matlab 的 ga 参数设置参考如下% gaOptions.m options optimoptions(ga, ... PopulationSize, 80, ... % 种群越大搜索越充分但耗时线性增长 MaxGenerations, 40, ... % 迭代代数先跑少一点看收敛曲线 Display, iter, ... % 打开迭代过程方便判断是否卡住 UseParallel, true, ... % 并行计算目标函数多核机器能快很多 PlotFcn, gaplotbestf, ... % 画最优值收敛曲线 FunctionTolerance, 1e-6, ... ConstraintTolerance, 1e-3);这里有两个参数最值得调PopulationSize 和 MaxGenerations。对于这种带阴影遮挡的目标函数每次评估可能几十毫秒种群太大会让单次迭代非常慢但种群太小又容易早熟。我自己一般先跑 20 代看曲线走势如果最优值还在下降就加大代数重启一轮。变量里如果有整数镜数 n_i需要设置 IntCon 参数ga 会使用特殊变异策略保证整数约束。等 ga 跑完可以把最优解里半径和高度提取出来固定整数镜数再用 fmincon 对连续变量做一次梯度精调往往还能提升 1% 到 3% 的目标值。4. Matlab 实现一份可以直接改参数跑通的代码骨架4.1 主脚本与初始参数所有参数集中放在一个 struct 里避免函数之间传递大量散装变量。下面是最小可运行的脚本骨架核心逻辑注释都写在代码里% main_heliofield.m % 定日镜场优化设计模型——可运行骨架 clear; clc; param.lat 40; % 纬度题目给定 param.towerH 80; % 塔高(m) param.receiverR 3.5; % 集热器半径(m) param.mirrorW 6; % 镜面宽(m) param.mirrorH 6; % 镜面高(m) param.heightRange [4, 12]; % 镜面中心离地高度范围(m) param.ringR1 60; % 最内圈半径(m) param.ringR2 300; % 最外圈半径(m) param.nRingMax 6; % 最多圈数 param.nMirrorMin 6; % 每圈最少镜数 param.nMirrorMax 40; % 每圈最多镜数 param.DNI 800; % 法向直射辐照度(W/m^2)题目有数据则替换 nvars param.nRingMax * 3; % 依次是 r(1..6), n(1..6), h(1..6) lb [repmat(param.ringR1, 1, param.nRingMax), ... repmat(param.nMirrorMin, 1, param.nRingMax), ... repmat(param.heightRange(1), 1, param.nRingMax)]; ub [repmat(param.ringR2, 1, param.nRingMax), ... repmat(param.nMirrorMax, 1, param.nRingMax), ... repmat(param.heightRange(2), 1, param.nRingMax)]; intcon param.nRingMax 1 : 2 * param.nRingMax;逻辑说明变量顺序是先把所有圈的半径排完再排所有圈的镜数最后排安装高度。这种编码方式在优化器眼里只是一个 1×18 的向量但我们的目标函数里会把它拆成三个向量分别构造镜场。intcon 只约束中间 6 个变量为整数半径和高度保持连续。参数说明nRingMax不要设太大6 圈已经能覆盖大多数镜场ringR1和ringR2要参考塔高塔高 80 米时最内圈取 50~60 米是合理的否则镜子会离塔太近反射光线角度过陡截断效率反而变差。4.2 布局生成与目标函数目标函数是整个优化的核心。它接收决策变量 x先拆出半径、镜数、高度然后生成所有镜面的位置再遍历代表时刻累加功率最后返回单位面积热功率的相反数。ga 是最小化所以目标函数前要加负号。% objFun.m function f objFun(x, param) nR param.nRingMax; r x(1:nR); nMir round(x(nR1 : 2*nR)); h x(2*nR1 : 3*nR); % 生成镜面结构体数组 mirrors buildLayout(r, nMir, h, param); % 简单惩罚没有有效镜面时直接返回大数 if isempty(mirrors) f 1e6; return; end % 计算年平均输出热功率 P_avg averagePower(mirrors, param); A_total 0; for i 1:numel(mirrors) A_total A_total mirrors(i).a; end % ga 最小化所以取负 f -P_avg / A_total; endbuildLayout 函数负责把半径、镜数、高度转换成每面镜子的位置和面积% buildLayout.m function mirrors buildLayout(r, nMir, h, param) mirrors struct(x, {}, y, {}, z, {}, w, {}, h, {}, a, {}); k 0; for i 1:length(r) if nMir(i) 0 || r(i) 0 continue; end for j 0 : nMir(i)-1 theta j * 2 * pi / nMir(i); k k 1; mirrors(k).x r(i) * cos(theta); mirrors(k).y r(i) * sin(theta); mirrors(k).z h(i); mirrors(k).w param.mirrorW; mirrors(k).h param.mirrorH; mirrors(k).a param.mirrorW * param.mirrorH; end end end这里把镜面简化为一个带坐标和面积的矩形x、y 是镜面中心在地面的投影z 是镜面中心离地高度。镜面本身的朝向不在结构体里保存因为朝向是随时变化的应当由太阳位置实时计算。averagePower 函数里要调太阳位置表下面给一个简化版% averagePower.m function P_avg averagePower(mirrors, param) towerTop [0, 0, param.towerH]; % 简化示例只取每年 5 个代表小时实际按题目数据扩展 P_total 0; for month 1:12 for hour 9:3:15 [alpha, gamma] solarPosition(param.lat, month, hour); sunVec sunVector(alpha, gamma); P evalPower(mirrors, sunVec, towerTop, param); P_total P_total param.timeWeight(month) * P; end end P_avg P_total; end注意这里 9:3:15 表示 9、12、15 三个时刻代表日取每月 21 日按月份天数加权后才能代表全年。实际比赛中题目给的是每个时刻的太阳高度角、方位角可能还有 DNI 数据直接把对应数值读进来替换solarPosition就行。4.3 阴影遮挡效率的近似实现阴影遮挡是计算量最大的环节也最容易在简化时出错。下面这段代码用“采样点 射线拦截”的方式估算遮挡比例逻辑直观适合比赛阶段的快速迭代% shadowBySampling.m % 对第 idx 面镜子在镜面上均匀采 N*N 个点 % 每个点沿入射光反向延伸到其他镜面统计被拦截比例 function etaSh shadowBySampling(idx, mirrors, sunVec, towerTop) N 5; % 采样密度5x5 已经是速度和精度折中 P [mirrors(idx).x, mirrors(idx).y, mirrors(idx).z]; normal sunVec (towerTop - P) / norm(towerTop - P); normal normal / norm(normal); t1 cross(normal, [0, 0, 1]); if norm(t1) 1e-6 t1 cross(normal, [0, 1, 0]); end t1 t1 / norm(t1); t2 cross(normal, t1); t2 t2 / norm(t2); blocked 0; total N * N; for i 1:N for j 1:N u (i - 0.5) / N - 0.5; v (j - 0.5) / N - 0.5; q P u * mirrors(idx).w * t1 v * mirrors(idx).h * t2; if isRayBlocked(q, sunVec, mirrors, idx) blocked blocked 1; end end end etaSh 1 - blocked / total; end function flag isRayBlocked(q, sunVec, mirrors, idx) flag false; rayDir -sunVec; % 从镜面出发逆着阳光方向找遮挡物 for j 1:numel(mirrors) if j idx continue; end Pj [mirrors(j).x, mirrors(j).y, mirrors(j).z]; % 简化用镜面中心和法线判断射线是否穿过矩形平面 d Pj - q; denom dot(rayDir, mirrors(j).normal); if abs(denom) 1e-8 continue; end t dot(d, mirrors(j).normal) / denom; if t 0 || t 500 continue; end hit q t * rayDir - Pj; if abs(dot(hit, mirrors(j).t1)) mirrors(j).w/2 ... abs(dot(hit, mirrors(j).t2)) mirrors(j).h/2 flag true; break; end end end这段代码的关键是t1、t2和normal构成了镜面局部坐标系。每个镜面的normal在真实运行时要预先计算好并且随太阳位置每步更新否则isRayBlocked里的平面方程不对。采样点取 5×5 已经能反映遮挡趋势如果要提交论文级结果建议改成多边形裁剪精确计算而不是采样。这里最大的坑是镜面数量超过 300 时每次目标评估都要执行 300×300×25 次射线判断速度会非常慢必须先用近邻搜索把候选遮挡镜面缩小到相邻几面。5. 定日镜场优化避坑效率计算翻车与迭代不收敛的排查5.1 阴影遮挡效率算出来恒等于 1现象无论布局怎么变目标函数里的阴影遮挡项始终是 1效率曲线没有任何波动。原因多半是isRayBlocked用的镜面法线没有随太阳位置更新。很多代码在初始化时算了一次法线之后就不再改变。但定日镜法线是太阳向量和反射向量的角平分线太阳每小时都在动法线必须每步重算。另一个常见原因是采样点q的局部坐标系方向设错了导致所有点都落不到其他镜面平面内。解决先单步调试固定某个时刻打印任意两面镜子的normal、t1、t2看投影点是否落在合理范围。最直接的办法是画三维镜场图把太阳方向画成箭头肉眼确认入射、反射、法线三者关系是否像一面真实的镜子。5.2 优化结果像随机数不同次运行布局完全不一样现象ga 每次运行得到的目标值差不多但最优布局细节差异很大半径忽大忽小镜数分布没有规律。原因这通常不是算法问题而是目标函数里有大量平坦区域。每圈镜数变化 1 面总功率变化很小目标函数对整数变量不敏感ga 在后期缺少压力。另一个原因是约束惩罚太轻无效布局和有效布局的目标值差距不够大。解决给目标函数加一个“镜面间距惩罚项”相邻圈半径差小于某个阈值时返回一个随着违规程度放大的惩罚值。另外把 PopulationSize 加大到 100 以上并设置随机种子至少在比赛期间能复现结果。也可以用rng(2023)固定随机种子。5.3 ga 迭代巨慢一次要跑几个小时现象代码逻辑没问题但每次迭代时间从几十秒涨到几分钟根本没法调参。原因第一阴影遮挡采样对每一面镜子都遍历所有镜子复杂度是 O(N^2)第二ga 每代评估种群 60 个个体也就是 60 次完整镜场功率计算第三全年小时数如果全遍历比如 8760 小时每次评估就是 8760 次效率和遮挡。三者相乘计算量瞬间爆炸。解决先建 KDTree 或按极角排序只检查相邻角度的镜子再把全年采样从每小时降为代表时刻。我的经验是每天取 9、12、15 点三个时刻每月取 21 日为代表日12 个月一共 36 次瞬时计算优化结果已经足够稳定。如果要求更细腻再加密到每 15 分钟一个点。5.4 截断效率远距离时反而增大不符合物理直觉现象镜场边缘镜子离塔 300 米光斑半径大于集热器但算出来的 eta_trunc 比中间区域还高显然不对。原因光斑半径公式用d * tan(4.65e-3 epsilon)epsilon 如果取 2 mrad总角在 300 米处光斑半径是 2 米左右而集热器半径 3.5 米-0.5*(R/r_spot)^2指数值接近 1截断效率接近 100%看似正常。但如果你把d * tan错写成了d / tan距离越远 r_spot 越小效率就会异常变大。这类低级错误很难用眼睛发现建议画一条“截断效率随距离变化”的曲线看是否单调递减。解决单独写一个测试函数固定 eta_cos、eta_atm 为 1扫描距离从 60 到 300 米输出 eta_trunc 随距离的曲线。正常应该从 1 缓慢下降不会出现上升段。5.5 年平均功率算出来比单时刻功率还高权重和采样没配对现象P_avg 算出来大于代表时刻的最大瞬时功率目标函数分母越大布局越倾向于堆面积。原因时间加权系数和月份天数没对齐。比如 2 月用了 31 天权重或者代表日采样点没覆盖整个白天导致加权和超过实际日照时数。题目给的数据集如果只有 9:00 到 15:00那 P_avg 应该理解为“这些代表时刻的平均值”而不是全天能量。解决把权重向量改为monthDays / sum(monthDays) * numMonths确保所有权重加起来等于采样时刻数。更稳妥的办法是只求“平均功率”而不是“累计能量”这样时间权重只影响相对占比不改变量级。6. 验证与进阶让结果经得起评委推敲优化做完了不能只给一个目标值就收工。我习惯在提交前做三件事效率分布可视化、参数敏感性扫描、局部精调对比。第一件事是把最优布局画成彩色散点图x、y 坐标对应镜子位置颜色对应单面镜的年平均热功率。这样能一眼看出哪些镜子贡献大、哪些镜子是“拖后腿”的。效率分布如果呈现明显的径向带状特征说明布局还有调整空间如果出现随机离散的高亮点大概率是阴影遮挡计算有误。第二件事是敏感性扫描。固定镜场布局不变把塔高连续性从 60 米扫到 100 米记录目标值变化曲线。你会发现目标值对塔高非常敏感塔高增加能明显降低遮挡但也会增加成本。比赛题通常不会要求你做经济性分析但评委很愿意看到这种“目标值随关键参数的变化趋势”这说明你不是只跑出一个数而是真的理解模型行为。第三件事是拿最优布局和一个均匀环形布局做对比。均匀布局的半径等差、每圈镜数按周长比例设置目标值差距通常在 8% 到 15% 之间。把这组对比写进论文比任何文字解释都有说服力。最后说一个我自己的教训。第一次跑这个题时我为了追求快把阴影遮挡效率直接设成常数 0.95结果优化器疯狂把镜子堆在塔附近目标值看起来很高实际上完全不可用。后来把遮挡项加回来最优布局才变成“中间密、边缘疏”的合理形态。效率和布局是强耦合的任何用常数代替物理项的地方最后都会在结果上露出马脚。希望这篇文章能让你少走这一步弯路也帮到你把这个模型从“能跑”做到“能信”。本文还有配套的精品资源点击获取