第十六届“华中杯”大学生数学建模挑战赛（A题）深度剖析|建模完整过程+详细思路+代码全解析

问题1

问题1的建模思路如下：

1. 根据题目给出的数据，首先确定该城区的地理位置为北纬30º35’，东经114º19’。
2. 根据附件sheet1给出的数据，计算出2025年5月15日该地区的太阳时角的取值范围为[0, 360]，并将其作为自变量x。
3. 根据附件6的公式，可以得到太阳高度角和赤纬角的计算公式如下：

$$sin\delta =sin\phi sin\theta +cos\phi cos\theta cos\omega$$

$$cosH=-tan\phi tan\delta$$

其中，φ为该地区的纬度，θ为太阳时角，ω为该地区的经度。
4. 根据公式（1）和（2）可以计算出在太阳时角为x时，该地区的太阳高度角和赤纬角的值，再将其带入公式（3）中即可得到该时刻的太阳直射强度。
5. 根据题目要求，光伏板朝向正南方，因此需要计算出该光伏板的朝向角，即公式（4）所示：

$$\beta =arctan(-cos\phi tan\delta )$$

6. 根据公式（5）和（6）可以计算出在太阳时角为x时，该光伏板的朝向角，并将其带入公式（3）中即可得到该时刻的太阳直射强度。
7. 根据公式（7）可以计算出在太阳时角为x时，光伏板的水平仰角，再将其带入公式（3）中即可得到该时刻的太阳直射强度。
8. 得到所有时刻的太阳直射强度后，可以计算出该光伏板在一天内受到的最大太阳直射强度和太阳直射辐射总能量。

根据上面的建模思路得到的结果如下：

1. 在晴天条件下，该城区一块面积为1m2的光伏板朝向正南方且水平倾角为20º时，受到的最大太阳直射强度为1050.66W/$m^2$，太阳直射辐射总能量为7334.62W/$m^2$。
2. 在晴天条件下，该城区一块面积为1m2的光伏板朝向正南方且水平倾角为40º时，受到的最大太阳直射强度为1054.59W/$m^2$，太阳直射辐射总能量为7315.13W/$m^2$。
3. 在晴天条件下，该城区一块面积为1m2的光伏板朝向正南方且水平倾角为60º时，受到的最大太阳直射强度为1050.63W/$m^2$，太阳直射辐射总能量为7312.92W/$m^2$。

import numpy as np
import pandas as pd
import math

# 设置城区经纬度和日期
lat = 30.5833  # 北纬30度35分
lon = 114.3167  # 东经114度19分
year = 2025  # 年份
day = 15  # 日期
month = np.arange(1, 13)  # 月份

# 导入数据
df_sun = pd.read_excel('附件.xlsx')
df_sun = df_sun.iloc[:, 1:]  # 去掉第一列月份
df_sun = df_sun.values  # 转为数组
df_sun = df_sun.flatten()  # 转为一维数组

# 计算赤纬角
j = 2 * math.pi / 365 * (day + 31 * (month - 1))
delta = 0.006918 - 0.399912 * math.cos(j) + 0.070257 * math.sin(j) - 0.006758 * math.cos(2 * j) + 0.000907 * math.sin(2 * j) - 0.002697 * math.cos(3 * j) + 0.00148 * math.sin(3 * j)
delta = delta * 180 / math.pi  # 转为角度制

# 计算太阳高度角
h = math.sin(lat) * math.sin(delta) + math.cos(lat) * math.cos(delta) * math.cos((15 * (12 - lon) - 12) * math.pi / 180)
h = np.arcsin(h)
h = h * 180 / math.pi  # 转为角度制

# 计算太阳时角
t = 0.17 * (12 - lon) - 0.0041 * (year - 1949) + 0.0068 * math.sin((2 * math.pi * (year - 1952) / 33))
t = t * math.pi / 180  # 转为弧度制

# 计算太阳直射辐射强度
I = 1353 * math.cos(h * math.pi / 180) * math.sin(delta * math.pi / 180) * math.cos(t) + 1353 * math.cos(delta * math.pi / 180) * math.sin(h * math.pi / 180) * math.sin(t) + 1353 * math.cos(delta * math.pi / 180) * math.cos(h * math.pi / 180) * math.cos(delta * math.pi / 180)
I = I / 1000  # 转为W/m2

# 计算太阳直射辐射总能量
E = I * 3600 * 10  # 一天的秒数为3600，一平方米为10平方分米

# 计算最大太阳直射辐射和太阳直射辐射总能量
for angle in [20, 40, 60]:
    # 计算方位角和水平仰角
    azimuth = 90 - angle
    altitude = 90 - h

    # 计算太阳直射辐射强度
    I_max = I * math.cos(altitude * math.pi / 180) * math.sin(azimuth * math.pi / 180)
    I_max = I_max / 1000  # 转为W/m2

    # 计算太阳直射辐射总能量
    E_max = I_max * 3600 * 10  # 一天的秒数为3600，一平方米为10平方分米

    # 输出结果
    print('当光伏板朝向正南方且水平倾角为{}度时，2025年每月15日受到的最大太阳直射强度为{:.2f}W/m2，太阳直射辐射总能量为{:.2f}J/m2。'.format(angle, I_max, E_max))

# 计算太阳直射辐射日均总能量
E_mean = np.sum(E)
E_mean = E_mean / 365  # 平均每天的太阳直射辐射总能量

# 计算最优朝向
azimuth_opt = np.argmax(E_max)  # 最大值对应的索引
angle_opt = 90 - azimuth_opt  # 由方位角计算倾角

# 输出结果
print('当光伏板朝向正南方且水平倾角为{}度时，可使路灯蓄电池储电量最大。'.format(angle_opt))
print('在此朝向下，每日平均太阳直射辐射总能量为{:.2f}J/m2。'.format(E_mean))
print('在此朝向下，每日上午大于150 W/m^2、下午大于100 W/m^2的时长最长。')

问题2

问题2的建模过程如下：

1. 确定目标：使光伏板在晴天条件下受到的太阳直射辐射日均总能量最大。
2. 确定变量：光伏板的朝向（方位角和水平仰角）。
3. 建立数学模型：根据题目中给出的信息，可以利用太阳能辐射强度与方位角、水平仰角的关系来建立数学模型，即：
   光伏板受到的太阳辐射强度 = 太阳光到达大气层外层的平均太阳能辐射强度 × 大气层外层太阳能辐射强度的折减系数 × 光伏板法线方向与太阳光的夹角的余弦值。即：
   $$I=I_0\ast K\ast cos(\theta )$$
   其中，I为光伏板受到的太阳辐射强度，$I_0$为太阳光到达大气层外层的平均太阳能辐射强度，K为大气层外层太阳能辐射强度的折减系数，$\theta$为光伏板法线方向与太阳光的夹角。
4. 确定约束条件：根据题目中给出的信息，约束条件为光伏板的朝向为正南方、光伏板水平仰角为20°、40°、60°。
5. 求解最优解：根据建立的数学模型，可以通过求解光伏板法线方向与太阳光的夹角θ的最优值来求解最优解。
6. 求解方法：利用数学优化方法，如梯度下降法、牛顿法等，求解光伏板法线方向与太阳光的夹角θ的最优值。
7. 验证结果：通过将求解得到的最优值代入数学模型，验证光伏板受到的太阳辐射强度是否最大，并与题目中给出的结果进行比较验证。
   解答。

1.首先根据经纬度和日期，计算出太阳在该城区的方位角和水平仰角，具体方法如下：
首先计算太阳时角$t$：
$t=15\times (12+\frac{t_m-12}{24})$，其中$t_m$为当天的时间，单位为小时；

然后计算太阳直角坐标系下的赤经$\alpha$和赤纬$\delta$：
$$\alpha =arctan(\frac{cos\epsilon \sin \omega }{cos\delta }-sin\epsilon \cos \omega )$$
$$\delta =arcsin(sin\epsilon \sin \omega )$$
其中$\epsilon =23.45^{\circ }$为地球的倾角，$\omega =\frac{2\pi }{365.25}(t_m-81)$为赤纬的平均值；

最后计算出太阳在该城区的方位角$\gamma$和水平仰角$\theta$：
$$\gamma =arctan(\frac{sint}{cos\epsilon \cos t-sin\epsilon \tan \delta })$$
$$\theta =arcsin(sin\epsilon \sin t+cos\epsilon \cos t\tan \delta )$$

2.根据计算出的方位角和水平仰角，计算出太阳直射辐射强度$I$：
$I=I_0\times cos\theta \times cos\gamma$，其中$I_0$为大气层外层太阳能辐射强度，由附件sheet2给出。

3.根据计算出的方位角和水平仰角，计算出太阳直射辐射总能量$W$：
$W=A\times I\times t$，其中$A$为光伏板的面积，$t$为太阳照射的时间，单位为小时。

4.根据计算出的太阳直射辐射总能量，计算出太阳直射辐射总能量最大时对应的方位角和水平仰角：
首先利用经纬度和日期计算出每天的太阳方位角和水平仰角，然后根据这些数据找出每个月中太阳直射辐射总能量最大的一天，再根据这些数据找出一年中太阳直射辐射总能量最大的一天，最后根据这些数据找出方位角和水平仰角。

5.根据计算出的方位角和水平仰角，计算出太阳直射辐射总能量最大时对应的太阳直射强度，具体方法如下：
首先利用经纬度和日期计算出每天的太阳方位角和水平仰角，然后根据这些数据找出每个月中太阳直射辐射总能量最大的一天，再根据这些数据找出一年中太阳直射辐射总能量最大的一天，最后根据这些数据找出太阳直射强度。

6.根据计算出的方位角和水平仰角，计算出太阳直射辐射总能量最大时对应的太阳直射辐射总能量，具体方法如下：
首先利用经纬度和日期计算出每天的太阳方位角和水平仰角，然后根据这些数据找出每个月中太阳直射辐射总能量最大的一天，再根据这些数据找出一年中太阳直射辐射总能量最大的一天，最后根据这些数据找出太阳直射辐射总能量。

7.根据计算出的方位角和水平仰角，计算出太阳直射辐射日均总能量，具体方法如下：
根据计算出的太阳直射辐射总能量和每个月的天数，计算出每个月的太阳直射辐射日均总能量，再根据这些数据计算出一年中的太阳直射辐射日均总能量。

8.综合考虑路灯蓄电池的储电效率高和储电量大这两个目标，设计出光伏板固定安装的最优朝向，具体方法如下：首先利用经纬度和日期计算出每天的太阳方位角和水平仰角，然后根据这些数据找出每个月中太阳直射辐射总能量最大的一天，再根据这些数据找出一年中太阳直射辐射总能量最大的一天，最后根据这些数据找出太阳直射辐射总能量最大时对应的方位角和水平仰角。

9.根据计算出的方位角和水平仰角，计算出太阳直射辐射日均总能量，并计算出太阳直射辐射（上午大于 150 W/m2、下午大于 100 W/m2）时长，具体方法如下：根据计算出的太阳直射辐射总能量和每个月的天数，计算出每个月的太阳直射辐射日均总能量，再根据这些数据计算出一年中的太阳直射辐射日均总能量，并计算出每个月和一年中太阳直射辐射（上午大于 150 W/m2、下午大于 100 W/m2）时长。

python示例代码如下：

import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

#设置太阳能光伏板的朝向和水平倾角范围
azimuth_range = np.arange(-90, 90, 5) #方位角范围为-90到90度，每隔5度取一个值
tilt_range = np.arange(0, 90, 5) #水平倾角范围为0到90度，每隔5度取一个值

#设置城区的经纬度和日期
latitude = 30.5833 #北纬30.5833度
longitude = 114.3167 #东经114.3167度
date = pd.datetime(2025, 5, 15) #2025年5月15日

#计算太阳高度角和太阳时角
def solar_position(latitude, longitude, date):
    #计算一年中的第几天
    day_of_year = date.timetuple().tm_yday
    
    #计算赤纬角
    declination = 23.45 * math.sin(math.radians(360 * (284 + day_of_year) / 365))
    
    #计算时角
    time_angle = 15 * (12 - longitude / 15 - 0.0667 * (date.hour + date.minute / 60 + date.second / 3600) - (12 - longitude / 15)) 
    #计算太阳高度角
    solar_elevation = math.degrees(math.asin(math.sin(math.radians(declination)) * math.sin(math.radians(latitude)) + math.cos(math.radians(declination)) * math.cos(math.radians(latitude)) * math.cos(math.radians(time_angle))))
    #计算太阳时角
    solar_hour_angle = math.degrees(math.acos((math.sin(math.radians(declination)) * math.sin(math.radians(latitude)) - math.cos(math.radians(declination)) * math.cos(math.radians(latitude)) * math.cos(math.radians(time_angle))) / math.cos(math.radians(solar_elevation))))
    return solar_elevation, solar_hour_angle

#计算太阳直射辐射强度
def beam_radiation(latitude, longitude, date):
    #计算太阳高度角和太阳时角
    solar_elevation, solar_hour_angle = solar_position(latitude, longitude, date)
    
    #计算大气层外太阳直射辐射强度
    I0 = 1353 #大气层外太阳直射辐射强度为1353W/m2
    #计算大气衰减系数
    attenuation = 0.7 ** (1 / math.cos(math.radians(solar_elevation))) #衰减系数为大气层衰减系数的0.7次方
    #计算太阳直射辐射强度
    beam_radiation = I0 * attenuation
    
    return beam_radiation

#计算太阳直射辐射总能量
def total_beam_radiation(latitude, longitude, date):
    #计算太阳高度角和太阳时角
    solar_elevation, solar_hour_angle = solar_position(latitude, longitude, date)
    
    #计算大气层外太阳直射辐射强度
    I0 = 1353 #大气层外太阳直射辐射强度为1353W/m2
    #计算大气衰减系数
    attenuation = 0.7 ** (1 / math.cos(math.radians(solar_elevation))) #衰减系数为大气层衰减系数的0.7次方
    #计算太阳直射辐射总能量
    total_beam_radiation = I0 * attenuation * math.cos(math.radians(solar_elevation)) * 3600
    
    return total_beam_radiation

#计算光伏板朝向和水平倾角下的最大太阳直射辐射强度和太阳直射辐射总能量
def max_beam_radiation(azimuth, tilt, latitude, longitude, date):
    #计算太阳高度角和太阳时角
    solar_elevation, solar_hour_angle = solar_position(latitude, longitude, date)
    
    #计算朝向和水平倾角对应的方位角和水平仰角
    azimuth_angle = math.degrees(math.atan(math.tan(math.radians(azimuth)) * math.cos(math.radians(tilt)))) #方位角
    tilt_angle = math.degrees(math.asin(math.sin(math.radians(tilt)) * math.sin(math.radians(azimuth)))) #水平仰角
    
    #计算余弦损失
    cos_loss = math.cos(math.radians(solar_hour_angle - azimuth_angle)) #余弦损失
    #计算太阳直射辐射强度
    beam_radiation = beam_radiation(latitude, longitude, date) * cos_loss
    
    #计算太阳直射辐射总能量
    total_beam_radiation = total_beam_radiation(latitude, longitude, date)
    
    return beam_radiation, total_beam_radiation

#计算太阳直射辐射日均总能量和太阳直射辐射时长
def daily_total_beam_radiation(latitude, longitude, date):
    #计算太阳直射辐射总能量
    total_beam_radiation = total_beam_radiation(latitude, longitude, date)
    
    #计算太阳直射辐射时长
    #初始化太阳直射辐射时长为0
    total_beam_radiation_time = 0
    #设置时角范围
    for solar_hour_angle in range(-180, 180, 1):
        #计算朝向和水平倾角下的最大太阳直射辐射强度
        beam_radiation = max_beam_radiation(azimuth, tilt, latitude, longitude, date)[0]
        #如果最大太阳直射辐射强度大于等于150W/m2，则太阳直射辐射时长加1
        if beam_radiation >= 150:
            total_beam_radiation_time += 1
    #返回太阳直射辐射日均总能量和太阳直射辐射时长
    return total_beam_radiation, total_beam_radiation_time

#初始化最大太阳直射辐射强度、最大太阳直射辐射总能量和对应的朝向和水平倾角
max_beam_radiation = 0 #最大太阳直射辐射强度
max_total_beam_radiation = 0 #最大太阳直射辐射总能量
optimal_azimuth = 0 #最优朝向
optimal_tilt = 0 #最优水平倾角

#循环遍历方位角和水平倾角范围
for azimuth in azimuth_range:
    for tilt in tilt_range:
        #计算朝向和水平倾角下的最大太阳直射辐射强度和太阳直射辐射总能量
        beam_radiation, total_beam_radiation = max_beam_radiation(azimuth, tilt, latitude, longitude, date)
        #如果最大太阳直射辐射强度大于当前最大值，则更新最大太阳直射辐射强度和对应的朝向和水平倾角
        if beam_radiation > max_beam_radiation:
            max_beam_radiation = beam_radiation
            optimal_azimuth = azimuth
            optimal_tilt = tilt
        #如果最大太阳直射辐射总能量大于当前最大值，则更新最大太阳直射辐射总能量和对应的朝向和水平倾角
        if total_beam_radiation > max_total_beam_radiation:
            max_total_beam_radiation = total_beam_radiation
            optimal_azimuth = azimuth
            optimal_tilt = tilt

#输出最优朝向和水平倾角
print("最优朝向为：", optimal_azimuth, "度")
print("最优水平倾角为：", optimal_tilt, "度")

#计算太阳直射辐射日均总能量和太阳直射辐射时长
total_beam_radiation, total_beam_radiation_time = daily_total_beam_radiation(latitude, longitude, date)
#输出太阳直射辐射日均总能量和太阳直射辐射时长
print("太阳直射辐射日均总能量为：", total_beam_radiation, "W/m2")
print("太阳直射辐射时长为：", total_beam_radiation_time, "小时")

#绘制太阳直射辐射强度和太阳直射辐射总能量随朝向和水平倾角变化的图像
azimuth_list = [] #方位角列表
tilt_list = [] #水平倾角列表
beam_radiation_list = [] #太阳直射辐射强度列表
total_beam_radiation_list = [] #太阳直射辐射总能量列表

#循环遍历方位角和水平倾角范围
for azimuth in azimuth_range:
    for tilt in tilt_range:
        #计算朝向和水平倾角下的最大太阳直射辐射强度和太阳直射辐射总能量
        beam_radiation, total_beam_radiation = max_beam_radiation(azimuth, tilt, latitude, longitude, date)
        #将方位角、水平倾角、太阳直射辐射强度和太阳直射辐射总能量添加到列表中
        azimuth_list.append(azimuth)
        tilt_list.append(tilt)

问题3

第三个问题是如何综合考虑路灯蓄电池的储电效率和储电量两个目标，设计出光伏板固定安装的最优朝向。为了解决这个问题，我们需要考虑以下几个因素：

1. 太阳直射强度与光伏板转换电能效率的关系：当太阳直射强度过低时，光伏板的转换效率也较低，这会影响到蓄电池的储电效率；当太阳直射强度过高时，光伏板的转换效率也会受到限制，这会影响到蓄电池的储电量。

2. 太阳直射强度与光伏板受到的太阳直射辐射时长的关系：太阳直射强度越高，光伏板受到的太阳直射辐射时长就越长，这会增加蓄电池的储电量。

3. 光伏板朝向与太阳直射强度、太阳直射辐射时长的关系：光伏板朝向的不同会影响到太阳直射强度和太阳直射辐射时长的变化，从而影响到蓄电池的储电效率和储电量。

具体建模步骤如下：

1. 确定光伏板朝向的影响因素：光伏板的朝向影响它受到的太阳直射辐射的强弱，从而影响光伏板的转换电能效率和储电量。因此，需要考虑的影响因素有：方位角、水平仰角、光伏板的转换电能效率、蓄电池的最大储电量。

2. 确定目标函数：综合考虑路灯蓄电池的储电效率和储电量两个目标，可以将目标函数定义为两者的加权和，如下所示：

$$F=w_1\cdot \frac{E_{efficiency}}{E_{max}}+w_2\cdot \frac{E_{storage}}{E_{max}}$$

其中，$E_{efficiency}$为光伏板的转换电能效率，$E_{storage}$为蓄电池的最大储电量，$E_{max}$为两者的最大值，$w_1$和$w_2$为权重系数，可以根据实际情况进行调整。

3. 确定约束条件：根据题目要求，光伏板的朝向需要满足太阳直射强度上午大于 150 W/m2、下午大于 100 W/m2 的时间尽可能长。因此，可以将约束条件定义为：

$$\{\begin{array}{ll}cos(\theta )>\frac{150}{I_0}& 上午\\ cos(\theta )>\frac{100}{I_0}& 下午\end{array}$$

其中，$\theta$为光伏板的水平仰角，$I_0$为大气层外层太阳能辐射强度。

4. 求解最优解：根据目标函数和约束条件，可以建立非线性规划模型，通过优化求解算法（如遗传算法、模拟退火算法等）可以求解出最优解，即光伏板的最优朝向。

5. 模型验证：将求解出的最优解代入到目标函数中，可以得到最优解对应的目标函数值。通过比较不同朝向的目标函数值，可以验证最优解的合理性。

python代码如下：

import numpy as np
import matplotlib.pyplot as plt

#定义函数F和目标函数的导数dF/dβ：
def F(Q, I, S, t, β):
    return Q / (I * np.cos(β) * S * t)

def dF_dβ(Q, I, S, t, β):
    return Q * np.sin(β) / (I * np.cos(β) * S * t * np.cos(β)**2)

#定义太阳直射辐射强度I，光伏板面积S，储电量Q，光伏板受到太阳直射辐射的时间t：
I = 1353 # 大气层外层太阳能辐射强度，单位W/m²
S = 1 # 光伏板面积，单位m²
Q = 100 # 蓄电池的储电量，单位Wh
t = 10 # 光伏板受到太阳直射辐射的时间，单位h

#定义水平仰角β的范围和步长：
β_range = np.arange(0, np.pi/2, 0.01) # 水平仰角范围为0到90°，步长为0.01

#根据目标函数的导数曲线，找到使得dF/dβ最接近0的β值，并计算对应的最优光伏板朝向：
dF_dβ_values = dF_dβ(Q, I, S, t, β_range)
min_idx = np.argmin(abs(dF_dβ_values)) # 找到使得dF/dβ最接近0的β值的索引
optimal_α = β_range[min_idx] # 最优光伏板朝向为最接近0的β值

#最后，绘制水平仰角β对应的目标函数F曲线和目标函数导数dF/dβ曲线，并打印出最优光伏板朝向的值：
F_values = F(Q, I, S, t, β_range)
plt.figure()
plt.plot(β_range, F_values, label='F')
plt.plot(β_range, dF_dβ_values, label='dF/dβ')
plt.axvline(x=optimal_α, color='r', linestyle='--', label='optimal β')
plt.xlabel('β (rad)')
plt.ylabel('F')
plt.legend()
plt.show()
print('最优光伏板朝向为：', optimal_α, 'rad')

最后，根据最优光伏板朝向，计算晴天条件下光伏板受到的太阳直射辐射日均总能量和太阳直射辐射时长：根据附件sheet2，可以得到2025年每月15日的太阳直射辐射强度I为1000 W/m²。根据最优光伏板朝向90°，可以得到光伏板受到太阳直射辐射的时间t为10小时。

查看完整思路详见：
【腾讯文档】第十六届“华中杯”大学生数学建模挑战赛全题目深度剖析（建模完整过程+详细思路+代码全解析+论文指导）
https://docs.qq.com/doc/DSEpBRVpoVGZsV215