2024年华中杯数学建模A题太阳能路灯光伏板的朝向设计问题解题全过程文档及程序

2024年华中杯数学建模

A题 太阳能路灯光伏板的朝向设计问题

原题再现

  太阳能路灯由太阳能电池板组件部分（包括支架）、LED灯头、控制箱（包含控制器、蓄电池）、市电辅助器和灯杆几部分构成。太阳能电池板通过支架固定在灯杆上端。太阳能电池板也叫光伏板， 它利用光伏效应接收太阳辐射能并转化为电能输出，经过充放电控制器储存在蓄电池中。太阳能辐射由直射辐射和散射辐射组成，其中直射辐射对聚集太阳能系统起到了至关重要的影响。大气层对太阳能直射辐射的衰减变化量与其辐射强度、所穿过的大气层厚度成正比，其中衰减系数（W/(m2/km)）反映了一个地区大气层的透光性能。通常地球表面大气层厚度按1000公里计算，大气层可视为包裹地球的球壳。太阳光到达大气层外层上的平均太阳能辐射强度I0为1353W/m2。受地球运行轨道及太阳光传播的距离影响，大气层外层太阳能辐射强度随时间发生改变。附件sheet2给出了1-12月份大气层外层太阳能辐射强度具体数值。
  安装光伏板的朝向直接影响到光伏板获得太阳辐射能量的多少。光伏板的朝向包括方位角和水平仰角，方位角为光伏板的法线在水平面上的投影与正南方向的夹角。并按如下方法规定：如果一个光伏板朝向正南，那么它的方位角为零；如果一个光伏板朝向正东，那么它的方位角为90º ；如果一个光伏板朝向正西，那么它的方位角为−90º；水平仰角为光伏电池板平面与水平面的夹角。当太阳光线和光伏板的法线方向一致时，光伏板瞬时受到的太阳照射能量最大，否则会有余弦损失。
  某城区地处北纬30‘，35’，东经114‘19’，附件sheet1给出了该城区2023年5月23日晴天状况下测得地表水平面受到的太阳直射强度值。关于赤纬角、太阳高度角、太阳时角等相关概念，可参见全国大学生数学建模竞赛2012B 题附件6、2015A题讲解和2023A题附录。请在仅考虑太阳直射辐射的情况下建立数模，回答如下问题：
  1. 请计算2025年每月15日，在晴天条件下，该城区一块面积为1m2的光伏板朝向正南方且水平倾角分别为20‘、40’、60‘时受到的最大太阳直射强度和太阳直射辐射总能量；
  2. 如果光伏板受到的太阳直射辐射总能量最大时，可使路灯蓄电池储电量最大。请设计该城区固定安装太阳能光伏板的朝向，使光伏板在晴天条件下受到的太阳直射辐射日均总能量最大；
  3. 当光板受到太阳直射强度过低时，它转换电能的效率也很低；而当光伏板受到太阳直射强度过高时，它转换电能实现储电的效率也会受到限制。理想的情况是，光伏板受到太阳直射强度上午大于150W/m2、下午大于100W/m2的时间尽可能长，这样可以使路灯蓄电池的储电效率更高。综合考虑路灯蓄电池的储电效率高和储电量大这两个目标，请设计出光伏板固定安装的最优朝向，并计算晴天条件下光伏板受到的太阳直射辐射日均总能量和太阳直射辐射（上午大于150 W/m2、下午大于100W/m2）时长。

整体求解过程概述(摘要)

  对于问题一，大气层对太阳能直射辐射的衰减变化量与其辐射强度、所穿过的大气层厚度成正比。对此，我们通过分析太阳光线变化过程中的几何关系，建立了太阳辐射衰减模型。我们利用所给的5月23日地表太阳直射强度数据，确定了衰减系数的值。然后，我们通过分析太阳光线与倾斜面的关系，建立了倾斜面上的辐射模型。通过分析斜面上的瞬时太阳辐射强度变化规律，可以得到相应的最大辐射强度。每日太阳辐射总能量即为所建立的函数对时间的积分。但由于函数关系较为复杂，我们使用离散型数值积分方法进行近似计算，得到了比较符合实际的结果。
  对于问题二，需要设计出固定安装太阳能光伏板的朝向，使光伏板在晴天条件下受到的太阳直射辐射日均总能量最大。我们分别基于遗传算法和模拟退火算法建立了优化模型。遗传算法和模拟退火算法都是启发式搜索算法，在处理一些复杂的的优化问题时具有良好的鲁棒性和通用性。基于两种模型的计算结果，通过多维比较，我们确定了在该题背景下基于模拟退火算法的优化模型得出的解更优，确定了太阳能光伏板的水平方位角和倾角分别为−0.23◦,27.86◦。
  对于问题三，为了解决该多目标优化问题，我们在遗传算法的基础上，设计出基于精英策略的多目标非支配排序遗传算法NSGA-II。通过快速非支配排序算法识别出非支配个体，同时通过精英策略，将父代种群和子代种群合并，从合并后的种群中根据非支配排序和拥挤度选择个体组成新的父代种群，确保了优秀个体能够保留到下一代中，从而提高了算法的收敛速度和解的精度。根据NSGA-II 得出的最优解分布以及 Pareto 前沿面，我们采用 CRITIC 权重法来确定储电效率和储电量的权重，得到方位角为0.84◦，倾角 28.26◦, 日均总能量和相应时长为 4396.54W·h，9.27h。

模型假设：

  假设1 我们只考虑太阳直射辐射对太阳能光伏板工作的影响。
  解释. 直射辐射是太阳能辐射中最主要的部分，为了简化问题，我们不再考虑其他类型的辐射。
  假设2 我们忽略温度对太阳能光伏板发电效率的影响。
  解释. 温度过高会显著影响太阳能光伏板的发电效率，但由于缺乏相关信息且为使问题尽可能简化，我们在此问题中暂不考虑温度的影响。
  假设3 我们仅在晴天的情况下进行分析，计算并由此确定最佳朝向。
  解释. 由于阴雨天气及极端天气出现较为随机，难以准确预测未来的天气状况，故我们仅在晴天的前提下进行计算。
  假设4 我们所获得的5月23日地表太阳直射强度和1-12月份大气层外层太阳能辐射强度的数据是可靠的，且每年相同日期的太阳辐射强度规律是相同的。
  解释. 参考资料中建立的公式仅与一年中的具体日期与时间点有关，未涉及年份的计算，我们暂且忽略年度差异。

问题分析：

  对于问题一，首先我们假定我们所获得的5月23日地表太阳直射强度和1-12月份大气层外层太阳能辐射强度的数据是可靠的，且每年相同日期的太阳辐射强度规律是相同的。由此我们可以利用题目中所给的两个数据进行分析与计算。由题目背景信息可知，大气层对太阳能直射辐射的衰减变化量与其辐射强度、所穿过的大气层厚度成正比。随着时间的变化，太阳辐射强度和太阳光线所穿过的大气层厚度也会随之改变。对此，我们通过分析太阳光线变化过程中的几何关系，建立了太阳辐射衰减模型。我们利用所给的5月23日地表太阳直射强度数据，确定了衰减系数的值，并由此计算出特定时间穿过大气层后的太阳辐射强度。然后，我们通过分析太阳光线与倾斜面的关系，建立了倾斜面上的辐射模型。对于问题中给定斜面的倾角，我们得到了任意时刻斜面上的瞬时太阳辐射强度。通过分析斜面上的瞬时太阳辐射强度变化规律，我们可以得到相应的最大辐射强度。而我们发现每日太阳辐射总能量即为所建立的函数对时间的积分。但由于函数关系较为复杂，我们使用离散型数值积分方法进行近似计算，得到了比较符合实际的结果。
  对于问题二，需要在问题一建立的模型基础上，设计出固定安装太阳能光伏板的朝向，使光伏板在晴天条件下受到的太阳直射辐射日均总能量最大。依据问题背景，这是一个双变量单目标优化问题，又考虑到此题非线性等实际情况，我们分别基于遗传算法和模拟退火算法建立了优化模型。遗传算法和模拟退火算法都是启发式搜索算法，在处理一些复杂的的优化问题时具有良好的鲁棒性和通用性。基于两种模型的计算结果，通过多维比较，我们确定了在该题背景下基于模拟退火算法的优化模型得出的解更优，从而确定了符合实际的太阳能光伏板的水平方位角和倾角。
  对于问题三，需要在前两个问题建立的模型基础上，设计出同时能满足路灯蓄电池的储电效率高和储电量大两个目标的光伏板固定安装最优朝向。根据问题背景，储电效率与光伏板受到太阳直射强度上午大于150W/m2、下午大于100W/m2的时间呈正相关关系，储电量同光伏板受到的太阳直射辐射日均总能量呈正相关关系。因此为了解决该多目标优化问题，我们在遗传算法的基础上，设计出基于精英策略的多目标非支配排序遗传算法NSGA-II。算法的主要策略是通过快速非支配排序算法识别出非支配个体，即Pareto最优解。同时通过精英策略，即通过将父代种群和子代种群合并，然后从合并后的种群中根据非支配排序和拥挤度选择个体组成新的父代种群，确保了优秀个体能够保留到下一代中，从而提高了算法的收敛速度和解的精度。根据NSGAII 得出的最优解分布以及Pareto 前沿面，我们采用 CRITIC 权重法（一种比熵权法和标准离差法处理效果更优的客观赋权法）来确定储电效率和储电量的权重，从而得出比较符合实际问题答案。

模型的建立与求解整体论文缩略图

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

部分程序代码：

import numpy as np
import math
K_legal = 0.0005624701263657791
I0=1334
b=6371
B=7371
phi=30.35
pi=3.1415926
standard_meridian=120
#计算光线直射穿过大气层长度
def calculate_length(phi, delta, omega, b, B):
phi_rad = np.radians(phi)
delta_rad = np.radians(delta)
omega_rad = np.radians(omega)
r_double_prime = b * np.cos(phi_rad - delta_rad)
L_1 = np.sqrt(B**2 - r_prime**2)
r_1 = r_double_prime * np.sin(omega_rad)
r_2 = r_double_prime * np.cos(omega_rad)
L_3 = np.sqrt(L_1**2 - r_1**2)
L_4 = L_3 - r_2
return L_4
#计算赤纬角
def calculate_declination(local_time):
day_of_year = local_time.timetuple().tm_yday
earth_tilt = 23.45
day_angle = 2 * np.pi * day_of_year / 365
sin_declination = np.sin(day_angle) * np.sin(np.radians(earth_tilt)) declination = np.degrees(np.arcsin(sin_declination))
return declination
#计算时角
def calculate_hour_angle(local_time):
time_difference = local_time - local_time.replace(hour=12, minute=0, second=0, microsecond=0) time_difference_in_hours = time_difference.total_seconds() / 3600
hour_angle = time_difference_in_hours * 15 # degrees
return hour_angle
#函数计算出日均太阳直射辐射强度比值
def Ra(x, y):
angle_x = np.radians(x) # 将角度转换为弧度
angle_y = np.radians(y) # 将角度转换为弧度
omega = calculate_hour_angle(local_time)
delta=calculate_declination(local_time)
phi1=math.radians(phi)
delta1=math.radians(delta)
omega1=math.radians(omega)
alpha=math.asin(math.sin(phi1)*math.sin(delta1)+math.cos(phi1)*math.cos(delta1)*math.cos( 
omega1))
if local_time.hour<12:
A = math.acos((math.sin(delta1) - math.sin(phi1) * math.sin(alpha)) / (math.cos(alpha) * math.cos(phi1)))
else:
A = math.acos((math.sin(delta1) - math.sin(phi1) * math.sin(alpha) ) / (math.cos(alpha)* math.cos(phi1)))
A=2*pi-A
Ra=-math.cos(alpha)*math.cos(A)*np.sin(angle_y)*np.cos(angle_x)+math.cos(alpha)*math.sin(A)*np
.sin(angle_x)*np.sin(angle_y)+np.cos(angle_y)*math.sin(alpha) 
if Ra<0:
Ra=0
return Ra

import mathfrom random import randomimport matplotlib.pyplot as pltfrom datetime import datetime, timedeltaimport numpy as np
K_legal = 0.0005624701263657791b=6371B=7371phi=30.35pi=3.1415926standard_meridian=120#计算光线直射穿过大气层长度def calculate_length(phi, delta, omega, b, B):phi_rad = np.radians(phi)delta_rad = np.radians(delta)omega_rad = np.radians(omega)r_prime = b * np.sin(phi_rad - delta_rad)r_double_prime = b * np.cos(phi_rad - delta_rad)L_1 = np.sqrt(B**2 - r_prime**2)r_1 = r_double_prime * np.sin(omega_rad)r_2 = r_double_prime * np.cos(omega_rad)L_3 = np.sqrt(L_1**2 - r_1**2)L_4 = L_3 - r_2return L_4#计算赤纬角def calculate_declination(local_time):day_of_year = local_time.timetuple().tm_ydayearth_tilt = 23.45 # Earth's tilt in degrees# Day angle in radiansday_angle = 2 * np.pi * day_of_year / 365# Calculating sin of declinationsin_declination = np.sin(day_angle) * np.sin(np.radians(earth_tilt))# Convert sin of declination to declination angle in degrees declination = np.degrees(np.arcsin(sin_declination))return declination
#计算时角def calculate_hour_angle(local_time):time_difference = local_time - local_time.replace(hour=12, minute=0, second=0, microsecond=0) time_difference_in_hours = time_difference.total_seconds() / 3600hour_angle = time_difference_in_hours * 15 # degreesreturn hour_angledef func(x,y):# 设置起始日期和结束日期sum2 = 0angle_x = np.radians(x) # 将角度转换为弧度angle_y = np.radians(y) # 将角度转换为弧度start_date = datetime(2024, 1, 1)end_date = datetime(2024, 12, 31)# 设置每天的起始时间和结束时间start_time_of_day = datetime(2024, 1, 1, 6, 00)end_time_of_day = datetime(2024, 1, 1, 19, 00)# 初始化一个列表来保存每一天的日期dates = []# 遍历每一天local_time = start_datewhile local_time <= end_date:dates.append(local_time.date()) # 只保存日期部分local_time += timedelta(days=1)# 现在，我们遍历每一天，并为每一天生成从早上8点到下午17点的每30分钟的时间点for date in dates:sum1 = 0local_time = datetime.combine(date, start_time_of_day.time())# print(local_time)while local_time <= datetime.combine(date, end_time_of_day.time()):# print(local_time)month_to_I0 = {1: 1405,2: 1394,3: 1378,4: 1353,5: 1334,6: 1316,7: 1308,
8: 1315,9: 1330,10: 1350,11: 1372,12: 1392}omega = calculate_hour_angle(local_time)delta = calculate_declination(local_time)phi1 = math.radians(phi)delta1 = math.radians(delta)omega1 = math.radians(omega)alpha = math.asin(math.sin(phi1) * math.sin(delta1) + math.cos(phi1) * math.cos(delta1) * math.cos( omega1))if local_time.hour < 12:A = math.acos((math.sin(delta1) - math.sin(phi1) * math.sin(alpha)) / (math.cos(alpha) * math. cos(phi1)))else:domain = (math.sin(delta1) - math.sin(phi1) * math.sin(alpha)) / (math.cos(alpha) * math.cos(phi1))if domain > 1:domain = 1if domain < -1:domain = -1A = math.acos(domain)A = 2 * pi - ARa1 = -math.cos(alpha) * math.cos(A) * np.sin(angle_y) * np.cos(angle_x) + math.cos( alpha) * math.sin(A) * np.sin(angle_x) * np.sin(angle_y) + np.cos(angle_y) * math.sin(alpha)if Ra1 < 0:Ra1 = 0if Ra1 > 1:Ra1 = 1I0 = month_to_I0.get(local_time.month)L = calculate_length(phi, delta, omega, b, B)I_values = I0 * np.exp(-K_legal * L) * Ra1# print(I_values)sum1 += I_values /2# 将当前时间增加30分钟local_time += timedelta(minutes=30)sum2 += sum1
return -sum2/365class SA:def __init__(self, func, iter=10, T0=100, Tf=2, alpha=0.9):self.func = funcself. iter = iter # 内循环迭代次数,即为L =100self.alpha = alpha # 降温系数，alpha=0.99self.T0 = T0 # 初始温度T0为100self.Tf = Tf # 温度终值Tf为0.01self.T = T0 # 当前温度self.x = [random() * 180 - 90 for i in range( iter)] # 随机生成100个x的值self.y = [random() * 90 for i in range( iter)] # 随机生成100个y的值self.most_best = []self.history = {'f': [], 'T': []}def generate_new(self, x, y): # 扰动产生新解的过程while True:x_new = x + self.T * (random() - random())y_new = y + self.T * (random() - random())if (-90 <= x_new <= 90) & (0 <= y_new <= 90):break # 重复得到新解，直到产生的新解满足约束条件return x_new, y_newdef Metrospolis(self, f, f_new): # Metropolis准则if f_new <= f:return 1else:p = math.exp((f - f_new) / self.T)if random() < p:return 1else:return 0def best(self): # 获取最优目标函数值f_list = [] # f_list数组保存每次迭代之后的值for i in range(self. iter):f = self.func(self.x[i], self.y[i])f_list.append(f)f_best = min(f_list)idx = f_list.index(f_best)return f_best, idx # f_best,idx分别为在该温度下，迭代L次之后目标函数的最优解和最优解的下 标
def run(self):count = 0# 外循环迭代，当前温度小于终止温度的阈值while self.T > self.Tf:# 内循环迭代100次for i in range(self. iter):f = self.func(self.x[i], self.y[i]) # f为迭代一次后的值x_new, y_new = self.generate_new(self.x[i], self.y[i]) # 产生新解 f_new = self.func(x_new, y_new) # 产生新值if self.Metrospolis(f, f_new): # 判断是否接受新值self.x[i] = x_new # 如果接受新值，则把新值的x,y存入x数组和y数组 self.y[i] = y_new# 迭代L次记录在该温度下最优解ft, _ = self.best()self.history['f'].append(ft)self.history['T'].append(self.T)# 温度按照一定的比例下降（冷却）self.T = self.T * self.alphacount += 1# 得到最优解print(self.T)f_best, idx = self.best()print(f"F={f_best}, x={self.x[idx]}, y={self.y[idx]}")sa = SA(func)sa.run()plt.plot(sa.history['T'], sa.history['f'])plt.title('SA')plt.xlabel('T')plt.ylabel('f')plt.gca().invert_xaxis()plt.show()

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可