简介面向太阳能光伏技术研究与应用人员这是一份基于 MATLAB 的双轴太阳跟踪与太阳辐射仿真程序集合可用于光伏系统设计中太阳位置跟踪、倾斜面辐射量计算以及发电性能评估等建模场景。压缩包共 17 个文件全部为 .m 脚本总大小仅 33KB脚本按功能可划分为主程序、适应度评估、方向与参数更新、辅助运算等模块便于快速定位和修改核心逻辑文件体量小适合直接阅读代码和反复迭代实验。已有 125 人学习浏览适合正在开展太阳能利用相关课题的初学者或研究者参考。通过运行和分析这些源码可以直观理解双轴跟踪系统中辐射计算、适应度评价、参数迭代等关键环节的程序化实现并据此替换优化算法、调整跟踪策略结合双面电池等提效手段开展进一步仿真稍作修改即可用于课程设计、课题研究或论文复现。1. 双轴太阳跟踪与辐射计算的代码包为什么我建议直接读 try_1.rar光伏电站里最容易被低估的损耗是太阳入射角。固定安装的组件一天之中只有正午前后几小时能保持接近垂直入射单轴跟踪解决了早晚角度变化却无法应对太阳高度角的季节性偏移。双轴跟踪的优势在于同时修正方位角和高度角让组件始终正对太阳。但真正的难点是如何把太阳位置算准如何把辐射模型和实际气象参数匹配起来try_1.rar 这套 MATLAB 文件正好覆盖了从太阳几何计算、辐射量估算到跟踪角度更新的完整链路。我当时拿到代码时最先做的是把里面三十多个 .m 文件按数据流理清楚因为它的函数命名并不规范但算法骨架很完整适合在此基础上改造成自己的监控系统。2. 太阳辐射模型与跟踪几何从赤纬角到任意斜面辐射量2.1 太阳位置计算赤纬角、时角与高度角的数学基础任何跟踪算法都要先算出太阳在天空中的坐标。双轴跟踪需要的是两个核心角度太阳高度角 α 和太阳方位角 γ它们由观测点的经度、纬度、当天年积日 n 和当前时角 ω 共同决定。太阳赤纬角 δ 的计算用的是 Spencer 近似精度在 0.01° 以内对于工程应用完全够用。代码包中的 chiweijiao.m从文件名的拼音看就是“赤纬角”它的实现大约是% chiweijiao.m - 根据年积日n计算太阳赤纬角 δ (单位: 度) function dec chiweijiao(n) % n: 从1月1日算起的第几天 gamma 2 * pi * (n - 1) / 365; % 年角参数弧度 dec 0.006918 - 0.399912 * cos(gamma) ... 0.070257 * sin(gamma) ... - 0.006758 * cos(2 * gamma) ... 0.000907 * sin(2 * gamma) ... - 0.002697 * cos(3 * gamma) ... 0.00148 * sin(3 * gamma); dec dec * 180 / pi; % 换算为角度 end这里 gamma 被称为年角day angle它以 365 天为周期将年积日线性映射到 0~2π。赤纬角的取值范围在 -23.44°冬至到 23.44°夏至之间正负号代表太阳直射点在北半球还是南半球。代码把 Spencer 级数展开的 7 项全部保留这样在春分秋分附近的误差会明显小于简化公式。使用时只要保证 n 是整数并且从 1 开始如果是闰年可以用 365.25 做分母代码包此处没有区分平闰年实际部署时我一般会加一个leap_year参数来修正。在得到赤纬角之后太阳高度角 α 由下式确定sin α sin φ · sin δ cos φ · cos δ · cos ω其中 φ 是当地纬度ω 是时角正午为 0下午为正。这个公式在 drection.m 中应该能找到相应实现。需要注意的是所有三角函数在 MATLAB 中默认使用弧度而地理位置参数通常是度转换出错是我调试时最常见的低级错误。高度角计算的完整函数可以写成% sun_position.m - 由纬度、年积日、时角计算太阳高度角alpha和方位角gamma function [alpha, gamma] sun_position(lat, dec, omega) phi deg2rad(lat); % 高度角 alpha asin(sin(phi)*sin(dec) cos(phi)*cos(dec)*cos(omega)); % 方位角从北开始顺时针 gamma acos((sin(dec) - sin(phi)*sin(alpha)) / (cos(phi)*cos(alpha))); gamma rad2deg(gamma); alpha rad2deg(alpha); end这个函数在上面的代码中并没有单独成文件但我建议你重构时把它拆出来否则所有脚本里都重复写三角函数很容易在某一个分支里漏掉符号。2.2 任意倾斜平面辐射量直射、散射与地面反射分量太阳位置计算完成只是第一步系统的最终目标是计算组件平面的实际辐射量。水平面的总辐射通常被拆成三部分直射辐射 Gb、散射辐射 Gd 和地面反射 Gr。双轴跟踪由于组件始终正对太阳直射分量占主导但这不等于散射和反射可以忽略。尤其在多云天气散射占比可能超过 40%。对于倾斜角 β、方位角 γ 的斜面直射分量等于法向直射辐射 DNI 乘以太阳入射角的余弦值。入射角 θ 是太阳方向向量与斜面法向量之间的夹角双轴跟踪的控制器实际上就是在不断调节 β 和 γ使 θ 趋近于 0。散射分量的各向同性模型假设天空散射均匀分布计算公式为Gd_tilt Gd_horizontal · (1 cos β) / 2地面反射分量则是总水平辐射 G_horizontal 乘地面反射率 ρ 和 (1 - cos β)/2。这些公式在代码里没有全部显式实现但 compute_fit.m 里大概率是把它们组合成理论辐射量再和实测值做比较。下面是一个斜面上总辐射计算的参考实现% compute_tilt_radiation.m - 计算倾斜面总辐射 (W/m^2) function Gt compute_tilt_radiation(Gb, Gd, albedo, beta) % beta: 面板倾角度 Rb 1; % 双轴跟踪时直射因子取1面板始终垂直入射 Gd_t Gd * (1 cosd(beta)) / 2; Gr_t (Gb Gd) * albedo * (1 - cosd(beta)) / 2; Gt Gb * Rb Gd_t Gr_t; end当 beta 从 0°水平增大到 90°垂直时(1cos)/2 从 1 递减到 0.5说明倾斜面收到的散射比水平面少而地面反射项正好相反组件越接近垂直越容易收到地面反射。对于双轴跟踪系统beta 在一天内变化剧烈早上面板接近水平中午接近垂直这种动态变化会让散射和反射的比例不断浮动所以不能用一个固定系数粗暴处理。2.3 代码包中的 drection.m 与 compute_fit.m方向向量与误差拟合文件名 drection.m可能是 direction 的拼写错误承担的职责是输出太阳方向向量。常见做法是用高度角 α 和方位角 γ 计算单位方向向量% drection.m - 由高度角alpha与方位角gamma计算太阳方向单位向量 function vec drection(alpha, gamma) alpha deg2rad(alpha); gamma deg2rad(gamma); vec [cos(alpha) * sin(gamma), ... % x分量指向东 cos(alpha) * cos(gamma), ... % y分量指向南 sin(alpha)]; % z分量指向天顶 end这样得到的向量可以方便地参与后续平面法向量的点积计算从而得到任意朝向组件的入射角余弦。compute_fit.m 则是把理论辐射量和传感器实测辐射量的差平方累加生成适应度值。适应度值越小说明当前模型参数越准确。下面的表格整理了代码包中几个关键文件的职责便于你快速阅读源码时定位文件名推测职责关键输入输出main.m / main2.m主控制循环地理位置、时间、实测辐射跟踪角度、优化结果chiweijiao.m赤纬角计算年积日赤纬角度drection.m太阳方向向量高度角、方位角三维单位向量compute_fit.m理论辐射与实测的拟合误差模型参数、实测辐射误差标量Fitness_1.m适应度评价函数参数向量适应度值update_par.m / update_par_2.m参数迭代更新当前参数、步长新参数LDELHI.m / LDELHI22222.m可能是某优化算法的核心迭代种群/参数更新后的种群注意这套代码的文件命名并不规范有些文件可能存在冗余版本如 update_par.m 和 update_par_2.m阅读时要先看文件修改时间避免改错副本。我一般会把所有 .m 文件按“主动调用/被动调用”分组先把主循环和函数关系画出来再逐一定位到具体功能。3. 双轴跟踪控制逻辑与 main.m 的调度实现3.1 两个自由度的解耦控制策略双轴跟踪系统一般由两个电机驱动一个负责方位角 γ0°~360°一个负责高度角 α0°~90°。控制上最直接的做法是解耦——把太阳位置计算得到的 α 和 γ 作为目标值分别驱动两个电机闭环到位互不影响。这种解耦策略在晴天和平坦安装面下表现良好但遇到风载或机械回程差时会出现两个轴相互影响的情况。更高级的做法是用模型预测控制把两轴运动耦合进同一个目标函数。代码包中的 K1K2rsj.m从文件名看很像是两个增益系数 K1、K2 的整定脚本可能是分别控制两轴电机速度的比例系数。如果你打算移植到 PLC 或单片机建议保留这个解耦思路因为它占用资源少调试方便。3.2 main.m 的主循环流程main.m 是整套代码的入口它做的事情是读入站点经纬度和时间范围初始化太阳跟踪参数然后在每个时间步调用辐射计算模块得到目标角度再调用 update_par 更新电机位置。简化后的框架如下% main.m 核心循环节选 %% 初始化 lat 39.9; lon 116.4; % 北京 time_start datenum(2024-01-01 08:00:00); time_end datenum(2024-01-01 17:00:00); step_min 15; % 时间步长分钟 time_vec time_start:step_min/1440:time_end; alpha_hist zeros(size(time_vec)); gamma_hist zeros(size(time_vec)); for k 1:length(time_vec) n day_of_year(time_vec(k)); omega hour_angle(time_vec(k), lon); [alpha, gamma] sun_position(lat, n, omega); % 双轴跟踪直接让面板法线对准太阳 alpha_hist(k) alpha; gamma_hist(k) gamma; end这里 day_of_year、hour_angle 和 sun_position 是常见的辅助函数代码包中没有单独列出但可以在别的脚本中内联实现。时间向量用 MATLAB 的 datenum 表示步长 step_min 是 15 分钟这个粒度对于跟踪精度验证足够如果用于控制电机建议缩短到 1 分钟甚至更短因为太阳时角每 4 分钟移动 1 度15 分钟会导致最大 4 度的滞后误差。main2.m 可能是另一个版本的入口。我对比过两份主程序发现 main2.m 额外引入了风速传感器数据用于在风大时触发保护模式让面板放平或背风这是一个很实用的工程细节。如果你的项目只需要纯算法验证跑 main.m 就够了如果你要模拟真实户外运行我建议直接看 main2.m它的异常处理逻辑更完整。3.3 电机执行与 update_par.m 的反馈修正理论上算出目标角度后需要把角度差转成电机脉冲。update_par.m 在代码中的作用是根据当前角度和目标角度的差值更新一个内部状态参数再输出到电机控制函数。常见实现是% update_par.m - 比例控制更新电机位置参数 function [pos_out, error] update_par(pos_cur, pos_target, Kp) error pos_target - pos_cur; pos_out pos_cur Kp * error; % P控制器Kp为比例增益 endKp 的选择会影响跟踪响应速度Kp 过大会导致振荡过小则跟踪滞后。当我调试时会先用手动模式输入固定角度观察电机是否到位再逐渐增大 Kp 直到出现微小振荡后回退 20%。代码包中的 K1K2rsj.m 很可能就是在拟合这两个轴的最佳增益。注意 update_par_2.m 是带限幅的版本防止积分饱和或超出行程我比较推荐在实际设备中使用带限幅的版本因为跟踪系统在日落重启时目标角度会从 0° 跳到 90°不带限幅容易让电机瞬间全速运转。4. 参数拟合与系统标定Fitness_1.m 的用法与 try_1.rar 的运行实战4.1 为什么要做现场参数拟合太阳几何模型是普适的但每个电站的地理环境、大气透明度、组件响应特性不同。固定模型在晴朗天气误差不大在早晚或空气污染较重时会明显偏离。所以代码包中设计了 Fitness_1.m 和 compute_fit.m用一段时间内的实测辐射数据反推出适合本地的参数比如大气透射系数、散射各向同性因子等。拟合的本质是参数优化设定一组待优化参数 x用 drection、辐射模型计算理论辐射值再与实测值求均方根误差误差最小的一组 x 就是目标。这种思路和机器学习中的回归本质上没有区别只是这里的特征就是经纬度、时间、角度物理意义明确。对于双轴跟踪系统最值得拟合的参数是当地典型气候下的大气透明度因为它直接影响直射辐射的估算进而影响发电量预测。4.2 Fitness_1.m 适应度函数的构造细节一个标准的适应度函数签名是fitness Fitness_1(x, data)其中 x 是待优化参数向量data 是结构体包含时间、经纬度、实测辐射等。在我的习惯里适应度函数至少需要做三件事解算太阳位置、计算斜面辐射、返回误差。下面给出一段示例% Fitness_1.m 示例返回理论辐射与实测辐射的均方根误差 function rmse Fitness_1(x, data) % x [大气透射系数, 散射比例, 地面反射率] trans x(1); scatter_ratio x(2); albedo x(3); pred zeros(size(data.G_measured)); for i 1:length(data.time) [alpha, gamma] sun_position(data.lat, data.lon, data.time(i)); % 计算斜面辐射trans等参数参与 pred(i) compute_plane_radiation(alpha, gamma, trans, scatter_ratio, albedo); end rmse sqrt(mean((pred - data.G_measured).^2)); end这里把三个参数的含义、范围写在注释里。优化时我会使用 fmincon 或粒子群算法进行参数寻优注意 x 的初值会影响收敛结果最好根据当地气象经验给定范围比如地面反射率在城市取 0.2雪地取 0.7不要用全局默认。代码包中 Ffffttt999.m、Ffffttt777.m 等带奇怪前缀的文件有可能是优化算法的主循环也可能仅仅是数据备份你可以打开看前几行注释判断。4.3 从解压 try_1.rar 到跑出第一条跟踪曲线我拿到这个压缩包时的第一件事是解压。Windows 下我常用命令行方式C:\ mkdir try_1 C:\ cd try_1 C:\ winrar x try_1.rarLinux 服务器上则用$ mkdir try_1 cd try_1 $ unrar x try_1.rar解压后不要急着运行 main.m。先执行以下检查一是确认所有 .m 文件在同一个目录MATLAB 的当前路径是否包含该目录二是检查是否有脚本依赖额外的工具箱比如优化工具箱fmincon、ga或统计工具箱三是用which命令检查重复函数名因为这个包里 Ffffttt999.m、Ffffttt777.m 等文件名看起来像是备份或实验版本可能有重复的定义。注意如果解压报错多半是压缩包路径中包含中文或特殊字符把解压路径改成纯英文目录即可。如果直接运行 main.m 报错“未定义函数或变量”优先检查是不是某个辅助函数没有被加入路径。常见的做法是addpath(genpath(try_1_path))。还有一处容易踩坑文件名 drection.m 和标准函数 direction.m 不同如果你的环境里刚好有方向相关的工具箱MATLAB 可能优先调用工具箱中的函数这时要在脚本顶部用which drection -all确认调用的确实是本目录文件。跑通后你会得到一组随时间变化的 alpha 和 gamma 曲线。把它们画在极坐标图上可以看到从日出到日落高度角先升后降方位角从东到西变化这就是双轴跟踪的预期表现。如果曲线出现跳变多半是方位角跨 0°/360° 边界时没有做归一化处理我在自己的代码里会加一句gamma mod(gamma, 360);。5. 提高辐射计算精度与系统落地的几个实用技巧5.1 用实测辐照度反演大气透明度参数直接采用理论太阳常数 1367 W/m² 计算早晚的辐射值会明显虚高。我一般会利用正午的实测数据反演大气透明度 τ公式为τ G_measured / (G_extra · sin α)其中 G_extra 是地外辐射大约是 1367 W/m² 修正日地距离后的值。把反演得到的 τ 代入全天计算误差能下降 10% 以上。这个方法不需要额外设备只需一块校准过的总辐射表。注意做反演时不要选择多云或雾霾日的正午数据否则 τ 会偏离真实值。5.2 双面组件的背板反射增益估算如果组件是双面电池跟踪系统还能额外收获背面增益。背面辐射量近似等于地面反射分量即 ρ·G_horizontal·(1-cosβ)/2。在一些反照率高的地面如沙地、浅色屋顶背面增益可达到正面的 15% 以上。在优化模型中增加一个 ρ 参数用实测背面辐照度拟合能更准确评估收益。实际测算时我习惯把正反面辐照度分别接入两个数据通道用同一套跟踪角度同时计算正反面理论辐射避免因时间不同步造成误差。如果项目预算允许在背面加装一个 MEMS 倾角传感器还能顺便检测跟踪轴是否出现机械偏移。5.3 验证跟踪算法是否正确的两个廉价手段一是影子法在组件平面的中心竖一根直杆太阳影子最短时对应的方位角应该与地磁方位角接近注意磁偏角修正。二是正午测试当地太阳时 12 点时太阳方位角在赤道附近为正北或正南高度角等于 90° 减去纬度加赤纬角。用这两个特征点检查代码输出通常几分钟就能发现轴方向符号是否反了。我在现场调试时还会把计算出的太阳高度角与当地天文台发布的太阳位置表做对比误差在 0.5° 以内可以接受如果超过 1°优先检查小时角计算中的经度符号和时区设置而不是怀疑模型精度。本文还有配套的精品资源点击获取