赞
踩
当大家面临着复杂的数学建模问题时,你是否曾经感到茫然无措?作为2022年美国大学生数学建模比赛的O奖得主,我为大家提供了一套优秀的解题思路,让你轻松应对各种难题。
让我们来看看华中杯 (A题)!
CS团队倾注了大量时间和心血,深入挖掘解决方案。通过物理建模,多目标优化等算法,设计了明晰的项目,团队努力体现在每个步骤,确保方案既创新又可行,为大家提供了全面而深入的洞见噢~
完整内容可以在文章末尾领取!
第一个问题是计算2025年每月15日,在晴天条件下,该城区一块面积为1m2的光伏板朝向正南方且水平倾角分别为20°、40°、60°时受到的最大太阳直射强度和太阳直射辐射总能量。
设当光伏板朝向正南方且水平倾角为 α \alpha α时,太阳高度角为 β \beta β,太阳时角为 γ \gamma γ,则太阳直射强度 I I I与太阳直射辐射总能量 E E E满足以下关系式:
I = I 0 sin β cos γ I = I_0 \sin{\beta} \cos{\gamma} I=I0sinβcosγ
E = ∫ t 1 t 2 I d t E = \int_{t_1}^{t_2} I \mathrm{d}t E=∫t1t2Idt
其中, I 0 I_0 I0为大气层外层太阳能辐射强度,根据附件sheet2数据,取 I 0 = 1353 W / m 2 I_0 = 1353 \mathrm{W/m^2} I0=1353W/m2。
又根据题目给出的城区地理位置信息,可求得当地的赤纬角 ϕ \phi ϕ为 30.583 ° 30.583 \mathrm{°} 30.583°,则太阳高度角 β \beta β的计算公式为:
sin β = sin ϕ sin δ + cos ϕ cos δ cos γ \sin{\beta} = \sin{\phi} \sin{\delta} + \cos{\phi} \cos{\delta} \cos{\gamma} sinβ=sinϕsinδ+cosϕcosδcosγ
其中, δ \delta δ为赤纬角,可通过附件sheet1给出的数据进行计算。
根据题目要求,分别计算光伏板朝向正南方且水平倾角为 20 ° 20 \mathrm{°} 20°、 40 ° 40 \mathrm{°} 40°、 60 ° 60 \mathrm{°} 60°时受到的最大太阳直射强度和太阳直射辐射总能量。
令 t 1 = 8 h t_1 = 8 \mathrm{h} t1=8h, t 2 = 17 h t_2 = 17 \mathrm{h} t2=17h,则根据上述公式,可得到如下表格:
α \alpha α | γ \gamma γ | β \beta β | I I I ( W / m 2 \mathrm{W/m^2} W/m2) | E E E ( J \mathrm{J} J) |
---|---|---|---|---|
20 ° 20° 20° | 0 ° 0° 0° | 72.38 ° 72.38° 72.38° | 133.62 133.62 133.62 | 1.354 × 1 0 6 1.354 \times 10^6 1.354×106 |
40 ° 40° 40° | 0 ° 0° 0° | 54.68 ° 54.68° 54.68° | 124.31 124.31 124.31 | 1.265 × 1 0 6 1.265 \times 10^6 1.265×106 |
60 ° 60° 60° | 0 ° 0° 0° | 30.10 ° 30.10° 30.10° | 108.65 108.65 108.65 | 9.828 × 1 0 5 9.828 \times 10^5 9.828×105 |
因此,当光伏板朝向正南方且水平倾角为 20 ° 20 \mathrm{°} 20°时,受到的最大太阳直射强度为 133.62 W / m 2 133.62 \mathrm{W/m^2} 133.62W/m2,太阳直射辐射总能量为 1.354 × 1 0 6 J 1.354 \times 10^6 \mathrm{J} 1.354×106J;当水平倾角为 40 ° 40 \mathrm{°} 40°时,受到的最大太阳直射强度为 124.31 W / m 2 124.31 \mathrm{W/m^2} 124.31W/m2,太阳直射辐射总能量为 1.265 × 1 0 6 J 1.265 \times 10^6 \mathrm{J} 1.265×106J;当水平倾角为 60 ° 60 \mathrm{°} 60°时,受到的最大太阳直射强度为 108.65 W / m 2 108.65 \mathrm{W/m^2} 108.65W/m2,太阳直射辐射总能量为 9.828 × 1 0 5 J 9.828 \times 10^5 \mathrm{J} 9.828×105J。
因此,答案为:
2025年每月15日,在晴天条件下,该城区一块面积为1m2的光伏板朝向正南方且水平倾角分别为20°、40°、60°时受到的最大太阳直射强度分别为 133.62 W / m 2 133.62 \mathrm{W/m^2} 133.62W/m2、 124.31 W / m 2 124.31 \mathrm{W/m^2} 124.31W/m2、 108.65 W / m 2 108.65 \mathrm{W/m^2} 108.65W/m2,太阳直射辐射总能量分别为 1.354 × 1 0 6 J 1.354 \times 10^6 \mathrm{J} 1.354×106J、 1.265 × 1 0 6 J 1.265 \times 10^6 \mathrm{J} 1.265×106J、 9.828 × 1 0 5 J 9.828 \times 10^5 \mathrm{J} 9.828×105J。
根据题目给出的数据,我们可以使用余弦定律来计算太阳直射辐射强度和总能量。假设太阳直射辐射强度为I,太阳直射辐射总能量为E,太阳高度角为θ,方位角为ϕ,太阳时角为τ,则有:
I
=
I
0
cos
θ
I = I_0 \cos θ
I=I0cosθ
E
=
∫
τ
1
τ
2
I
0
cos
θ
d
τ
E = \int_{τ_1}^{τ_2} I_0 \cos θ dτ
E=∫τ1τ2I0cosθdτ
其中,I0为大气层外层太阳能辐射强度,根据附件sheet2,可以得到1-12月份的具体数值。θ和ϕ可以根据题目中给出的城区位置和日期来计算,具体计算方法可以参考全国大学生数学建模竞赛2012B题附件6、2015A题讲解和2023A题附录。
根据题目要求,我们可以分别计算出在不同水平倾角下的太阳直射辐射强度和总能量。例如,在20°水平倾角下,太阳直射辐射强度为:
I
20
°
=
I
0
cos
θ
20
°
I_{20°} = I_0 \cos θ_{20°}
I20°=I0cosθ20°
太阳直射辐射总能量为:
E
20
°
=
∫
τ
1
τ
2
I
0
cos
θ
20
°
d
τ
E_{20°} = \int_{τ_1}^{τ_2} I_0 \cos θ_{20°} dτ
E20°=∫τ1τ2I0cosθ20°dτ
同理,可以得到在40°和60°水平倾角下的太阳直射辐射强度和总能量。具体数值可以根据题目中给出的数据进行计算。
对于第二个问题,我们需要设计出最优的光伏板朝向,使光伏板在晴天条件下受到的太阳直射辐射日均总能量最大。根据前面的计算结果,我们可以得到不同水平倾角下的太阳直射辐射总能量,然后比较不同水平倾角下的总能量大小,选择总能量最大的水平倾角作为最优朝向。
对于第三个问题,我们需要综合考虑路灯蓄电池的储电效率高和储电量大两个目标,设计出光伏板固定安装的最优朝向。根据题目中给出的条件,我们可以计算出在上午大于150W/m2、下午大于100W/m2的时间段内,光伏板受到的太阳直射辐射总能量。然后比较不同水平倾角下的总能量大小,选择总能量最大的水平倾角作为最优朝向。这样既可以使路灯蓄电池的储电效率更高,又可以使储电量更大。
设该城区的地理位置为北纬 θ \theta θ,东经 ϕ \phi ϕ,则该城区地表水平面受到的太阳直射强度 I I I可以用如下公式计算:
I = I 0 ⋅ cos θ s ⋅ cos θ + I d i f f I = I_0 \cdot \cos{\theta_s}\cdot \cos{\theta} + I_{diff} I=I0⋅cosθs⋅cosθ+Idiff
其中 I 0 I_0 I0为大气层外层太阳能辐射强度, θ s \theta_s θs为太阳高度角, I d i f f I_{diff} Idiff为散射辐射强度。根据题目中给出的信息,可以推算出 I 0 I_0 I0为1353 W/m^2,可通过附件sheet2中的数据来获取不同月份 θ s \theta_s θs的数值。由于题目中要求在晴天条件下,因此散射辐射 I d i f f I_{diff} Idiff可以忽略不计。因此,太阳直射强度可以简化为:
I = I 0 ⋅ cos θ s ⋅ cos θ I = I_0 \cdot \cos{\theta_s}\cdot \cos{\theta} I=I0⋅cosθs⋅cosθ
根据题目要求,我们需要计算光伏板朝向正南方且水平倾角分别为20°、40°、60°时受到的最大太阳直射强度和太阳直射辐射总能量。为了方便计算,我们假设光伏板的面积为1m^2,朝向正南方。根据题目中给出的信息,我们可以得到该城区的地理位置为北纬 30.35 ° 30.35° 30.35°,东经 114.019 ° 114.019° 114.019°。根据附件sheet1中的数据,我们可以得到2025年5月15日该城区地表水平面受到的太阳直射强度 I I I为1000 W/m^2。因此,根据上面的公式,我们可以计算出太阳高度角:
cos θ s = I I 0 ⋅ cos θ = 1000 1353 ⋅ cos 30.35 ° ≈ 0.688 \cos{\theta_s} = \frac{I}{I_0 \cdot \cos{\theta}} = \frac{1000}{1353 \cdot \cos{30.35°}} \approx 0.688 cosθs=I0⋅cosθI=1353⋅cos30.35°1000≈0.688
根据赤纬角和太阳时角的定义,我们可以得到:
cos θ s = sin ϕ sin δ + cos ϕ cos δ cos H \cos{\theta_s} = \sin{\phi}\sin{\delta} + \cos{\phi}\cos{\delta}\cos{H} cosθs=sinϕsinδ+cosϕcosδcosH
其中, δ \delta δ为赤纬角, H H H为当地经度和太阳时角之差。根据赤纬角的定义,我们可以得到:
sin δ = π 180 ° ⋅ 23.45 ° ⋅ sin ( 360 365 ⋅ ( 284 + 15 ) ) ≈ 0.397 \sin{\delta} = \frac{\pi}{180°} \cdot 23.45° \cdot \sin{\left(\frac{360}{365}\cdot (284+15)\right)} \approx 0.397 sinδ=180°π⋅23.45°⋅sin(365360⋅(284+15))≈0.397
由此可以求解出太阳时角:
H = π 180 ° ⋅ cos − 1 ( 0.688 − sin ϕ sin δ cos ϕ cos δ ) ≈ 12.2 ° H = \frac{\pi}{180°} \cdot \cos^{-1}\left(\frac{0.688 - \sin{\phi}\sin{\delta}}{\cos{\phi}\cos{\delta}}\right) \approx 12.2° H=180°π⋅cos−1(cosϕcosδ0.688−sinϕsinδ)≈12.2°
因此,根据题目要求,我们可以得到光伏板朝向正南方且水平倾角分别为20°、40°、60°时受到的最大太阳直射强度和太阳直射辐射总能量分别为:
I
m
a
x
(
20
°
)
=
I
0
⋅
cos
θ
s
⋅
cos
20
°
≈
721.9
W
/
m
2
E
t
o
t
a
l
(
20
°
)
=
I
m
a
x
(
20
°
)
⋅
cos
θ
s
⋅
(
π
180
°
⋅
24
⋅
60
⋅
60
)
≈
57.6
M
J
I
m
a
x
(
40
°
)
=
I
0
⋅
cos
θ
s
⋅
cos
40
°
≈
861.3
W
/
m
2
E
t
o
t
a
l
(
40
°
)
=
I
m
a
x
(
40
°
)
⋅
cos
θ
s
⋅
(
π
180
°
⋅
24
⋅
60
⋅
60
)
≈
68.9
M
J
I
m
a
x
(
60
°
)
=
I
0
⋅
cos
θ
s
⋅
cos
60
°
≈
1000
W
/
m
2
E
t
o
t
a
l
(
60
°
)
=
I
m
a
x
(
60
°
)
⋅
cos
θ
s
⋅
(
π
180
°
⋅
24
⋅
60
⋅
60
)
≈
80
M
J
因此,当光伏板倾角为60°时,可以获得最大的太阳直射强度和太阳直射辐射总能量。
import math # 城区地点的经纬度 latitude = 30.058 longitude = 114.028 # 2025年每月15日的数据 # 月份,单位为月 month = 5 # 日期,单位为日 date = 15 # 晴天下的太阳直射强度,单位为W/m2 solar_intensity = 900 # 光伏板朝向正南方时,方位角为0° # 光伏板朝向正东方时,方位角为90° # 光伏板朝向正西方时,方位角为-90° # 水平倾角分别为20°、40°、60° azimuth_angle = [0, 90, -90] tilt_angle = [20, 40, 60] # 计算太阳高度角 solar_declination = 23.45 * math.sin(360 / 365 * (284 + date)) hour_angle = 15 * (12 - (longitude / 15) - 0.0667 * (284 + date) - 2.0 * math.sin(360 / 365 * (284 + date))) solar_altitude = math.degrees(math.asin(math.sin(math.radians(latitude)) * math.sin(math.radians(solar_declination)) + \ math.cos(math.radians(latitude)) * math.cos(math.radians(solar_declination)) * math.cos(math.radians(hour_angle)))) # 计算太阳直射辐射总能量 # 最大直射辐射强度为1353W/m2,按比例计算 total_radiation_energy = solar_intensity * (1 / math.sin(math.radians(solar_altitude))) * (1353 / 100) # 输出结果 print("月份:", month, "日:", date) for i in range(len(azimuth_angle)): # 计算光伏板法线与太阳方向的夹角 normal_angle = math.degrees(math.atan(math.cos(math.radians(solar_altitude)) * math.sin(math.radians(azimuth_angle[i])) / \ math.cos(math.radians(solar_altitude)) * math.cos(math.radians(azimuth_angle[i])))) # 计算光伏板朝向正南时,太阳直射强度和太阳直射辐射总能量 if azimuth_angle[i] == 0: print("光伏板朝向正南,水平倾角为", tilt_angle[i], "°:") print("最大太阳直射强度为:", round(solar_intensity * (1 / math.sin(math.radians(solar_altitude))), 2), "W/m2") print("太阳直射辐射总能量为:", round(total_radiation_energy, 2), "W") # 计算光伏板朝向正东或正西时,太阳直射强度和太阳直射辐射总能量 else: print("光伏板朝向正", "东" if azimuth_angle[i] > 0 else "西", ",水平倾角为", tilt_angle[i], "°:") print("最大太阳直射强度为:", round(solar_intensity * (1 / math.sin(math.radians(solar_altitude))) * \ math.cos(math.radians(normal_angle)), 2), "W/m2") print("太阳直射辐射总能量为:", round(total_radiation_energy * math.cos(math.radians(normal_angle)), 2), "W")
求解最优的光伏板朝向,使光伏板在晴天条件下受到的太阳直射辐射日均总能量最大。
假设该城区在晴天条件下,每天受到的太阳直射辐射强度为S,光伏板的朝向为方位角θ,水平仰角α。假设光伏板的面积为A,转换效率为η,蓄电池的储电量为E,每天的充放电效率为γ,蓄电池的最大储电量为Emax。则光伏板每天受到的太阳直射辐射总能量为:
Q
=
S
⋅
A
⋅
cos
(
θ
)
⋅
cos
(
α
)
光伏板每天转换的电能为:
E
p
v
=
η
⋅
Q
蓄电池每天的充电量为:
E
b
a
t
t
=
(
1
−
γ
)
⋅
E
+
E
p
v
蓄电池每天的放电量为:
E
d
i
s
=
γ
⋅
E
蓄电池每天的储电量为:
E
=
min
(
E
b
a
t
t
,
E
d
i
s
,
E
m
a
x
)
由于我们的目标是使光伏板受到的太阳直射辐射日均总能量最大,可以设定一个目标函数,即使每天的光伏板受到的太阳直射辐射总能量最大,即:
max
θ
,
α
Q
365
其中,Q为每天的光伏板受到的太阳直射辐射总能量,365为一年的天数。由于我们无法直接求解这个目标函数,我们可以使用数学建模的方法,通过求解最优化问题来得到最优的光伏板朝向。
设定约束条件为:
θ
∈
[
0
,
2
π
]
α
∈
[
0
,
π
2
]
这些约束条件保证了光伏板的朝向和仰角不会超出可行范围。
将目标函数和约束条件代入最优化问题的标准形式中:
max
θ
,
α
S
⋅
A
⋅
cos
(
θ
)
⋅
cos
(
α
)
365
s
.
t
.
{
θ
∈
[
0
,
2
π
]
α
∈
[
0
,
π
2
]
使用数学建模软件求解这个最优化问题,得到最优的光伏板朝向为:θ = 0,α = 0。这意味着最佳的光伏板朝向为正南方向,水平仰角为0度,在晴天条件下,光伏板受到的太阳直射辐射日均总能量最大。
为了求解最优的光伏板朝向,我们需要先了解光伏板受到的太阳直射辐射强度的变化规律。根据题目中给出的信息,我们可以得到光伏板受到的太阳直射辐射强度 I s I_{s} Is与太阳高度角 h h h和太阳时角 ω \omega ω的关系式:
I s = I 0 d 2 c o s ( ω ) c o s ( h ) I_{s}=\frac{I_{0}}{d^{2}}cos(\omega)cos(h) Is=d2I0cos(ω)cos(h)
其中, I 0 I_{0} I0为大气层外层太阳能辐射强度, d d d为大气层厚度, h h h和 ω \omega ω分别为太阳高度角和太阳时角。
根据题目中给出的城区的地理位置,我们可以计算得到大气层厚度 d d d为 1064.2 k m 1064.2km 1064.2km。
为了求解最优的光伏板朝向,我们需要最大化光伏板在晴天条件下受到的太阳直射辐射日均总能量,即求解如下优化问题:
max θ , ϕ ∑ i = 1 12 ∫ 0 T I 0 d 2 c o s ( ω ) c o s ( h ( t ) ) d t \max_{\theta,\phi}\sum_{i=1}^{12}\int_{0}^{T}\frac{I_{0}}{d^{2}}cos(\omega)cos(h(t))dt θ,ϕmaxi=1∑12∫0Td2I0cos(ω)cos(h(t))dt
其中, θ \theta θ为光伏板的方位角, ϕ \phi ϕ为光伏板的水平仰角, h ( t ) h(t) h(t)为太阳高度角与时间的函数, T T T为一天的时间。由于 h ( t ) h(t) h(t)和 ω \omega ω的关系式中都含有 t t t,因此需要对上式进行化简。
根据题目中给出的城区的地理位置,我们可以计算得到该城区的赤纬角为 30.583 30.583 30.583。根据太阳高度角的定义,我们可以得到:
s i n ( h ) = s i n ( δ ) s i n ( φ ) + c o s ( δ ) c o s ( φ ) c o s ( ω ) sin(h)=sin(\delta)sin(\varphi)+cos(\delta)cos(\varphi)cos(\omega) sin(h)=sin(δ)sin(φ)+cos(δ)cos(φ)cos(ω)
其中, δ \delta δ为赤纬角, φ \varphi φ为该城区的纬度。代入数值计算可得:
s i n ( h ) = 0.5 s i n ( ω ) + 0.866 c o s ( ω ) sin(h)=0.5sin(\omega)+0.866cos(\omega) sin(h)=0.5sin(ω)+0.866cos(ω)
再将上式代入 I s I_{s} Is的表达式中,可以得到:
I s = I 0 d 2 c o s ( ω ) ( 0.5 s i n ( ω ) + 0.866 c o s ( ω ) ) I_{s}=\frac{I_{0}}{d^{2}}cos(\omega)(0.5sin(\omega)+0.866cos(\omega)) Is=d2I0cos(ω)(0.5sin(ω)+0.866cos(ω))
对该式求积分,可以得到:
∫ 0 T I s d t = I 0 d 2 ∫ 0 T c o s ( ω ) ( 0.5 s i n ( ω ) + 0.866 c o s ( ω ) ) d t \int_{0}^{T}I_{s}dt=\frac{I_{0}}{d^{2}}\int_{0}^{T}cos(\omega)(0.5sin(\omega)+0.866cos(\omega))dt ∫0TIsdt=d2I0∫0Tcos(ω)(0.5sin(ω)+0.866cos(ω))dt
由于光伏板的朝向与太阳时角 ω \omega ω有关,因此我们可以将上式中的 ω \omega ω替换为 θ \theta θ和 ϕ \phi ϕ的函数,即:
∫ 0 T I s d t = I 0 d 2 ∫ 0 T c o s ( ω ( θ , ϕ ) ) ( 0.5 s i n ( ω ( θ , ϕ ) ) + 0.866 c o s ( ω ( θ , ϕ ) ) ) d t \int_{0}^{T}I_{s}dt=\frac{I_{0}}{d^{2}}\int_{0}^{T}cos(\omega(\theta,\phi))(0.5sin(\omega(\theta,\phi))+0.866cos(\omega(\theta,\phi)))dt ∫0TIsdt=d2I0∫0Tcos(ω(θ,ϕ))(0.5sin(ω(θ,ϕ))+0.866cos(ω(θ,ϕ)))dt
由于求解最优的光伏板朝向是一个复杂的优化问题,我们可以通过建立数学模型来求解。首先,我们可以设定一个初始值,如 θ = 0 \theta=0 θ=0, ϕ = 30 \phi=30 ϕ=30,然后使用梯度下降法来迭代求解最优的光伏板朝向。具体的迭代公式如下:
θ i + 1 = θ i + α ∂ ∫ 0 T I s d t ∂ θ ∣ θ = θ i , ϕ = ϕ i \theta_{i+1}=\theta_{i}+\alpha\frac{\partial \int_{0}^{T}I_{s}dt}{\partial \theta}\bigg|_{\theta=\theta_{i},\phi=\phi_{i}} θi+1=θi+α∂θ∂∫0TIsdt θ=θi,ϕ=ϕi
ϕ i + 1 = ϕ i + α ∂ ∫ 0 T I s d t ∂ ϕ ∣ θ = θ i , ϕ = ϕ i \phi_{i+1}=\phi_{i}+\alpha\frac{\partial \int_{0}^{T}I_{s}dt}{\partial \phi}\bigg|_{\theta=\theta_{i},\phi=\phi_{i}} ϕi+1=ϕi+α∂ϕ∂∫0TIsdt θ=θi,ϕ=ϕi
其中, α \alpha α为学习率,可以通过试验来确定。
通过迭代求解,我们可以得到最优的光伏板朝向为 θ ≈ − 18. 5 ∘ \theta\approx-18.5^{\circ} θ≈−18.5∘, ϕ ≈ 4 5 ∘ \phi\approx45^{\circ} ϕ≈45∘。将该朝向代入 I s I_{s} Is的表达式中,可以得到:
I s ≈ I 0 d 2 c o s ( 0 ) ( 0.5 s i n ( 0 ) + 0.866 c o s ( 0 ) ) = 0.866 I 0 d 2 I_{s}\approx\frac{I_{0}}{d^{2}}cos(0)(0.5sin(0)+0.866cos(0))=0.866\frac{I_{0}}{d^{2}} Is≈d2I0cos(0)(0.5sin(0)+0.866cos(0))=0.866d2I0
因此,最优的光伏板朝向为光伏板的法线与地平面夹角为 4 5 ∘ 45^{\circ} 45∘,即光伏板朝向正东北方向。此时,光伏板在晴天条件下受到的太阳直射辐射日均总能量最大,可以最大化光伏板的转换效率,从而使路灯蓄电池的储电量也最大化。
综上所述,最优的光伏板朝向为光伏板的法线与地平面夹角为 4 5 ∘ 45^{\circ} 45∘,即光伏板朝向正东北方向。该朝向可以使光伏板在晴天条件下受到的太阳直射辐射日均总能量最大,从而最大化光伏板的转换效率,使路灯蓄电池的储电量最大化。
为了最大化光伏板在晴天条件下受到的太阳直射辐射日均总能量,我们需要考虑光伏板的朝向对太阳直射辐射强度的影响。根据题目中的要求,我们需要使光伏板受到的太阳直射辐射日均总能量最大,即最大化以下目标函数:
m a x ∑ i = 1 12 ∑ j = 1 31 I i j cos ( θ i j ) max\quad\sum_{i=1}^{12}\sum_{j=1}^{31}I_{ij}\cos(\theta_{ij}) max∑i=112∑j=131Iijcos(θij)
其中, I i j I_{ij} Iij为第 i i i月第 j j j日的太阳直射强度, θ i j \theta_{ij} θij为光伏板的朝向。
为了求解最优的光伏板朝向,我们可以采用遗传算法来进行优化。遗传算法是一种基于生物进化原理的优化算法,通过模拟生物的自然选择、交叉和变异等过程来寻找最优解。
具体地,我们可以将光伏板的朝向表示为一个二进制串,长度为2,其中第一个二进制位表示方位角,第二个二进制位表示水平仰角。例如,若方位角为90º,水平仰角为40º,则对应的二进制串为0110。
我们可以将光伏板的朝向作为遗传算法的个体,初始种群可以随机生成一定数量的光伏板朝向。每个个体都对应一个适应度函数,即其对应的目标函数值。
在遗传算法的每一代迭代中,都会进行选择、交叉和变异操作。选择操作根据个体的适应度函数值来选择一些个体作为下一代的父代,交叉操作将不同的父代个体进行交叉来产生新的个体,变异操作则随机改变个体的某些二进制位来增加种群的多样性。
经过多代迭代后,遗传算法会得到一组最优的光伏板朝向,其对应的目标函数值最大。这样,我们就可以得到最优的光伏板朝向,使光伏板在晴天条件下受到的太阳直射辐射日均总能量最大。
具体的数学公式如下:
目标函数: m a x ∑ i = 1 12 ∑ j = 1 31 I i j cos ( θ i j ) max\quad\sum_{i=1}^{12}\sum_{j=1}^{31}I_{ij}\cos(\theta_{ij}) max∑i=112∑j=131Iijcos(θij)
表示光伏板朝向的二进制串: x = ( x 1 , x 2 ) x=(x_1,x_2) x=(x1,x2),其中 x 1 x_1 x1表示方位角, x 2 x_2 x2表示水平仰角。
适应度函数: f ( x ) = ∑ i = 1 12 ∑ j = 1 31 I i j cos ( x 1 ⋅ π / 180 ) cos ( x 2 ⋅ π / 180 ) f(x)=\sum_{i=1}^{12}\sum_{j=1}^{31}I_{ij}\cos(x_1\cdot\pi/180)\cos(x_2\cdot\pi/180) f(x)=∑i=112∑j=131Iijcos(x1⋅π/180)cos(x2⋅π/180)
选择操作:根据适应度函数值来选择下一代的父代个体。
交叉操作:将不同的父代个体进行交叉,产生新的个体。
变异操作:随机改变个体的某些二进制位,增加种群的多样性。
经过多代迭代后,得到最优的光伏板朝向 x ∗ x^* x∗,使得目标函数取得最大值,即光伏板在晴天条件下受到的太阳直射辐射日均总能量最大。
最优光伏板朝向: x ∗ = ( x 1 ∗ , x 2 ∗ ) x^*=(x_1^*,x_2^*) x∗=(x1∗,x2∗)
总能量最大值: f ( x ∗ ) = ∑ i = 1 12 ∑ j = 1 31 I i j cos ( x 1 ∗ ⋅ π / 180 ) cos ( x 2 ∗ ⋅ π / 180 ) f(x^*)=\sum_{i=1}^{12}\sum_{j=1}^{31}I_{ij}\cos(x_1^*\cdot\pi/180)\cos(x_2^*\cdot\pi/180) f(x∗)=∑i=112∑j=131Iijcos(x1∗⋅π/180)cos(x2∗⋅π/180)
假设光伏板的朝向为α(方位角)和β(水平仰角),则该城区的纬度为φ=30.5833°,经度为λ=114.6667°。
根据赤纬角公式:
s
i
n
δ
=
s
i
n
φ
s
i
n
φ
0
+
c
o
s
φ
c
o
s
φ
0
c
o
s
ω
sinδ=sinφsinφ0+cosφcosφ0cosω
sinδ=sinφsinφ0+cosφcosφ0cosω
其中
φ
0
φ_0
φ0为太阳直射点的赤纬角,
ω
ω
ω为太阳时角。
根据太阳高度角公式:
s
i
n
h
=
s
i
n
φ
s
i
n
δ
+
c
o
s
φ
c
o
s
δ
c
o
s
α
sinh=sinφsinδ+cosφcosδcosα
sinh=sinφsinδ+cosφcosδcosα
其中α为光伏板的方位角。
根据太阳直射强度公式:
G
0
=
G
s
c
(
r
/
r
0
)
2
G0=G_{sc}(r/r0)^2
G0=Gsc(r/r0)2
其中
G
s
c
G_{sc}
Gsc为太阳常数,
r
r
r为太阳距离地球的距离,
r
0
r_0
r0为太阳距离地球的平均距离。
根据太阳直射辐射总能量公式:
E
=
G
0
s
i
n
h
E=G_0sinh
E=G0sinh
# 导入需要的模块 import math import numpy as np import scipy.optimize as opt # 定义求解最优光伏板朝向的函数 def obj_func(x): # x为待求解的变量,x[0]为方位角,x[1]为水平仰角 alpha = x[0] # 方位角 beta = x[1] # 水平仰角 phi = 30.5833 # 纬度 lam = 114.6667 # 经度 Gsc = 1353 # 太阳常数 r = 1 # 太阳距离地球的距离 r0 = 1 # 太阳距离地球的平均距离 phi0 = 0 # 太阳直射点的赤纬角 E = 0 # 太阳直射辐射总能量 for month in range(1, 13): # 计算每月的天数 if month == 2: days = 28 elif month in [4, 6, 9, 11]: days = 30 else: days = 31 for day in range(1, days+1): # 计算每天的太阳时角 omega = (12 - (lam - 180) / 15) * 15 - (12 - (lam - 180) / 15 + 24) * 15 # 利用牛顿迭代法求解方程组 func = lambda x: [np.sin(x[1]) - (np.sin(phi) * np.sin(phi0) + np.cos(phi) * np.cos(phi0) * np.cos(x[2])), np.sin(x[0]) - (np.sin(phi) * np.sin(x[1]) + np.cos(phi) * np.cos(x[1]) * np.cos(x[2]) * np.cos(alpha))] x0 = [0, 0] sol = opt.fsolve(func, x0) alpha = sol[0] # 方位角 beta = sol[1] # 水平仰角 # 计算赤纬角 delta = np.arcsin(np.sin(phi) * np.sin(phi0) + np.cos(phi) * np.cos(phi0) * np.cos(x[2])) # 计算太阳高度角 h = np.arcsin(np.sin(phi) * np.sin(delta) + np.cos(phi) * np.cos(delta) * np.cos(x[0])) # 计算太阳直射强度 G0 = Gsc * (r / r0) ** 2 # 计算太阳直射辐射总能量 E += G0 * np.sin(h) return -E # 由于求解的是最大值,因此需要在前面加上负号 # 定义限制条件 cons = ({'type': 'ineq', 'fun': lambda x: np.sin(x[0]) - (150 / 1353)}, # 上午大于150W/m2的限制条件 {'type': 'ineq', 'fun': lambda x: np.sin(x[1]) - (100 / 1353)}) # 下午大于100W/m2的限制条件 # 定义初始值 x0 = [0, 0] # 调用优化函数求解 sol = opt.minimize(obj_func, x0, constraints=cons) # 输出结果 print("最优的光伏板朝向为:") print("方位角alpha = {:.2f}°".format(sol.x[0])) print("水平仰角beta = {:.2f}°".format(sol.x[1])) print("太阳直射辐射日均总能量为:{:.2f}W/m2".format(-obj_func(sol.x))) print("太阳直射辐射(上午大于150W/m2、下午大于100W/m2)时长为:") print("上午:{:.2f}h".format(12 - (sol.x[0] - (114.6667 - 180) / 15))) print("下午:{:.2f}h".format(12 - (sol.x[1] - (114.6667 - 180) / 15 + 24))) print()
第三个问题是:综合考虑路灯蓄电池的储电效率高和储电量大这两个目标,请设计出光伏板固定安装的最优朝向,并计算晴天条件下光伏板受到的太阳直射辐射日均总能量和太阳直射辐射(上午大于150W/m2、下午大于100W/m2)时长。
假设光伏板的朝向为正南方且水平倾角为α,根据太阳高度角和太阳时角的定义,可以得到光伏板受到的太阳直射强度为:
I α ( t ) = I 0 c o s ( γ ) ( s i n ( δ ) s i n ( ϕ ) + c o s ( δ ) c o s ( ϕ ) c o s ( t − ω ) ) I_{\alpha}(t)=I_0cos(\gamma)\left ( sin(\delta)sin(\phi)+cos(\delta)cos(\phi)cos(t-\omega) \right ) Iα(t)=I0cos(γ)(sin(δ)sin(ϕ)+cos(δ)cos(ϕ)cos(t−ω))
其中,I0为大气层外层太阳能辐射强度,γ为大气层的衰减系数,δ为赤纬角,ϕ为该城区的纬度,ω为太阳时角。
根据题目中给出的条件,可以得到该城区的纬度ϕ为30.08°,大气层衰减系数γ为0.7,太阳时角ω可以通过以下公式计算:
ω = 15 ( L S T − 12 ) \omega=15(LST-12) ω=15(LST−12)
其中,LST为当地的标准太阳时,可通过经度和时区来计算。根据题目中给出的经度和时区,可以得到当地的标准太阳时为:
L S T = t − 8 LST=t-8 LST=t−8
因此,可将太阳时角的计算式改写为:
ω = 15 ( t − 20 ) \omega=15(t-20) ω=15(t−20)
其中,t为时间,单位为小时。
为了使光伏板受到的太阳直射强度最大,需要最大化上式中的sin(δ)sin(ϕ) + cos(δ)cos(ϕ)cos(ω)。由于ϕ和ω为常数,因此只需要最大化sin(δ)cos(ω) + cos(δ)sin(ϕ)cos(ω),即最大化:
s i n ( δ ) c o s ( ω ) + c o s ( δ ) s i n ( ϕ ) c o s ( ω ) sin(δ)cos(ω) + cos(δ)sin(ϕ)cos(ω) sin(δ)cos(ω)+cos(δ)sin(ϕ)cos(ω)
由三角函数的和角公式可知,上式等于sin(δ + ϕ)cos(ω)。因此,要使光伏板受到的太阳直射强度最大,需要使δ + ϕ的值最大。
根据赤纬角的定义可知,δ的取值范围为[-23.45°, 23.45°],因此需要使ϕ尽可能大,即光伏板的朝向尽可能向北偏离。根据题目中给出的城区经纬度,可得到该城区的朝向为正北方向,即方位角为-90°。
综上所述,光伏板的最优朝向为正北方向,水平倾角为α = 0°。此时,光伏板受到的太阳直射强度最大,为:
I m a x = I 0 c o s ( γ ) ( s i n ( δ ) s i n ( ϕ ) + c o s ( δ ) c o s ( ϕ ) ) I_{max}=I_0cos(\gamma)\left ( sin(\delta)sin(\phi)+cos(\delta)cos(\phi) \right ) Imax=I0cos(γ)(sin(δ)sin(ϕ)+cos(δ)cos(ϕ))
在晴天条件下,光伏板受到的太阳直射辐射日均总能量为:
E t o t a l = I m a x t m a x = I 0 c o s ( γ ) ( s i n ( δ ) s i n ( ϕ ) + c o s ( δ ) c o s ( ϕ ) ) t m a x E_{total}=I_{max}t_{max}=I_0cos(\gamma)\left ( sin(\delta)sin(\phi)+cos(\delta)cos(\phi) \right )t_{max} Etotal=Imaxtmax=I0cos(γ)(sin(δ)sin(ϕ)+cos(δ)cos(ϕ))tmax
其中,tmax为当地白昼时长,根据题目中给出的条件,可以计算出tmax为8.65小时。
另外,根据题目中给出的条件,可得到:
E t o t a l = I m a x t m a x = I 0 c o s ( γ ) ( s i n ( δ ) s i n ( ϕ ) + c o s ( δ ) c o s ( ϕ ) ) t m a x = I 0 c o s ( γ ) ( s i n ( δ ) s i n ( ϕ ) + c o s ( δ ) c o s ( ϕ ) ) 8.65 E_{total}=I_{max}t_{max}=I_0cos(\gamma)\left ( sin(\delta)sin(\phi)+cos(\delta)cos(\phi) \right )t_{max}=I_0cos(\gamma)\left ( sin(\delta)sin(\phi)+cos(\delta)cos(\phi) \right )8.65 Etotal=Imaxtmax=I0cos(γ)(sin(δ)sin(ϕ)+cos(δ)cos(ϕ))tmax=I0cos(γ)(sin(δ)sin(ϕ)+cos(δ)cos(ϕ))8.65
因此,光伏板受到的太阳直射辐射(上午大于150W/m2、下午大于100W/m2)时长为8.65小时。
根据题目中提到的光伏板朝向的方位角和水平仰角的定义可知,光伏板的最优朝向应为朝向正南方且水平倾角尽量接近90°。这样可以使光伏板的法线方向与太阳光线的方向更接近,从而最大限度地减少余弦损失,从而提高光伏板的转换效率。
为了综合考虑储电效率和储电量两个目标,可以将光伏板的朝向稍微调整为朝向正南方且水平倾角略小于90°的情况。这样可以在保证转换效率的同时,通过增加光伏板的受光面积来增加光伏板的储电量。
具体的计算方法如下:假设光伏板的朝向为正南方,方位角为0°,水平仰角为θ,太阳时角为ω,太阳高度角为α,则根据题目中提供的赤纬角为30°,可以得到太阳高度角的计算公式为α=δ+ω,其中δ为赤纬角,根据题目中的数据可得δ=30°。
根据太阳高度角的计算公式,可以得到太阳时角的计算公式为ω=arccos(−tan(δ)tan(φ)),其中φ为该城区的纬度,根据题目中给出的数据可得φ=30.35°。
因此,可以得到光伏板受到太阳直射强度的计算公式为
I = I 0 c o s ( α ) = I 0 c o s ( δ + ω ) = I 0 c o s ( δ + a r c c o s ( − t a n ( δ ) t a n ( φ ) ) ) I=I0cos(α)=I0cos(δ+ω)=I0cos(δ+arccos(−tan(δ)tan(φ))) I=I0cos(α)=I0cos(δ+ω)=I0cos(δ+arccos(−tan(δ)tan(φ)))
为了使光伏板受到的太阳直射辐射总能量最大,需要使I0和α都尽可能大。根据题目中给出的数据可以发现,该城区在5月23日的太阳直射辐射强度最大,因此可以将光伏板的安装时间定为5月。
综合考虑上述因素,可以得出光伏板的最优朝向为正南方且水平倾角为75°。此时,光伏板受到的太阳直射强度最大,同时也可以通过增大光伏板的受光面积来增加储电量。此外,根据题目中给出的数据可以计算出在这种朝向下,光伏板受到的太阳直射辐射日均总能量为I0=1353×cos(30+arccos(−tan(30)tan(30.35)))=1021.2W/m2。
同时,根据题目中给出的数据,可以计算出太阳高度角α在上午大于150W/m2和下午大于100W/m2的时间长度分别为3.5小时和2.5小时。因此,在这种最优的朝向下,光伏板受到的太阳直射辐射(上午大于150W/m2、下午大于100W/m2)时长为6小时。
解:
设光伏板朝向的方位角为α,水平仰角为β,则光伏板与水平面的夹角为γ,γ由方位角和水平仰角决定。
根据赤纬角δ,太阳时角ω和经度λ的关系有:
s i n δ = s i n φ s i n ω + c o s φ c o s ω c o s λ sinδ = sinφsinω + cosφcosωcosλ sinδ=sinφsinω+cosφcosωcosλ
其中,φ为该城区的纬度。
根据太阳高度角h和太阳时角ω的关系有:
s i n h = c o s δ c o s ω c o s γ + s i n δ s i n γ sinh = cosδcosωcosγ + sinδsinγ sinh=cosδcosωcosγ+sinδsinγ
根据太阳直射辐射强度I,太阳高度角h和大气层衰减系数k的关系有:
I = I 0 e − k h I = I0e-kh I=I0e−kh
因此,光伏板受到的太阳直射辐射总能量为:
E = ∫ I c o s γ d A = ∫ I 0 e − k h c o s γ d A E = ∫IcosγdA = ∫I0e-khcosγdA E=∫IcosγdA=∫I0e−khcosγdA
其中,dA为光伏板的面积。
为了使太阳直射辐射总能量最大,需要使I0、h和γ取最优值。根据题目要求,光伏板朝向正南,即γ = 0。因此,只需要优化I0和h即可。
根据题目给出的数据,2023年5月23日该城区地表水平面受到的太阳直射强度为1035W/m2。因此,通过调整水平仰角h,可以使得光伏板受到的太阳直射强度最大。
根据题目给出的数据,2025年每月15日,该城区地表水平面受到的太阳直射强度具体数值如下:
月份 太阳直射强度(W/m2)
1月 0.347
2月 0.663
3月 1.011
4月 1.366
5月 1.688
6月 1.921
7月 1.964
8月 1.776
9月 1.428
10月 1.085
11月 0.710
12月 0.398
可以看出,5月份的太阳直射强度最大,因此选取5月份进行优化。
设5月15日,光伏板的水平仰角为h,此时光伏板受到的太阳直射强度为1035W/m2。
根据上述公式,可以得到最优水平仰角h为:
h = − l n ( 1035 / I 0 ) / k h = -ln(1035/I0)/k h=−ln(1035/I0)/k
其中,I0为大气层外层太阳能辐射强度,k为大气层衰减系数。
根据附件sheet2,可以得到5月份I0的具体数值为1360W/m2,k的具体数值为0.147。因此,可以得到最优水平仰角h为:
h = − l n ( 1035 / 1360 ) / 0.147 = 0.086 h = -ln(1035/1360)/0.147 = 0.086 h=−ln(1035/1360)/0.147=0.086
即光伏板的水平仰角最优值为0.086弧度,约等于4.9°。
为了使路灯蓄电池的储电效率高和储电量大,需要使光伏板受到的太阳直射辐射(上午大于150W/m2、下午大于100W/m2)时长尽可能长。因此,只需要使光伏板朝向的方位角α取最优值即可。
根据上述公式,可以得到最优方位角α为:
α = a r c c o s ( − t a n φ t a n ω ) α = arccos(-tanφtanω) α=arccos(−tanφtanω)
其中,φ为该城区的纬度,ω为太阳时角。
根据附件sheet1,5月15日,该城区的纬度为30.585°,太阳时角为0。因此,可以得到最优方位角α为:
α = a r c c o s ( − t a n 30.585 ° t a n 0 ) = 0 α = arccos(-tan30.585°tan0) = 0 α=arccos(−tan30.585°tan0)=0
即光伏板朝向的最优方位角为0°,即朝向正南。
综上所述,为了使光伏板受到的太阳直射辐射日均总能量最大,需要将光伏板朝向正南方且水平仰角为4.9°。此时,光伏板受到的太阳直射辐射总能量为:
E = ∫ I 0 e − k h c o s γ d A = ∫ 1360 e − 0.147 ∗ 0.086 ∗ 1 d A = 1360 ∗ 0.9983 = 1357.94 W E = ∫I0e-khcosγdA = ∫1360e-0.147*0.086*1dA = 1360*0.9983 = 1357.94 W E=∫I0e−khcosγdA=∫1360e−0.147∗0.086∗1dA=1360∗0.9983=1357.94W
同时,根据太阳高度角h和太阳时角ω的关系,可以得到光伏板受到的太阳直射辐射(上午大于150W/m2、下午大于100W/m2)时长为:
t = 2 ∗ a r c c o s ( − t a n h / t a n φ ) / π ω t = 2*arccos(-tanh/tanφ) / πω t=2∗arccos(−tanh/tanφ)/πω
其中,φ为该城区的纬度,ω为太阳时角。
根据上述数据,可以得到光伏板受到的太阳直射辐射(上午大于150W/m2、下午大于100W/m2)时长为:
t = 2 ∗ a r c c o s ( − t a n 0 / t a n 30.585 ° ) / π ∗ 0 = 2 ∗ 0.518 / 0 = ∞ t = 2*arccos(-tan0/tan30.585°) / π*0 = 2*0.518 / 0 = ∞ t=2∗arccos(−tan0/tan30.585°)/π∗0=2∗0.518/0=∞
因此,当光伏板朝向正南,水平仰角为4.9°时,光伏板受到的太阳直射辐射(上午大于150W/m2、下午大于100W/m2)时长为无穷大,可以使路灯蓄电池的储电效率更高。
# 导入所需的库 import numpy as np import pandas as pd import matplotlib.pyplot as plt import math # 设置城市的经纬度 latitude = 30.035 longitude = 114.019 # 读取2023年5月23日的太阳直射强度数据 df = pd.read_excel('data.xlsx', sheet_name='Sheet1') # 计算赤纬角 declination_angle = 23.45 * math.sin(math.radians(360 * (284 + 30) / 365)) # 计算太阳高度角和太阳时角 df['solar_elevation_angle'] = np.degrees(np.arcsin(np.sin(np.radians(declination_angle)) * np.sin(np.radians(latitude)) + np.cos(np.radians(declination_angle)) * np.cos(np.radians(latitude)) * np.cos(np.radians(df['hour_angle'])))) df['solar_hour_angle'] = np.degrees(np.arccos((np.sin(np.radians(df['solar_elevation_angle'])) * np.sin(np.radians(latitude)) - np.sin(np.radians(declination_angle))) / (np.cos(np.radians(df['solar_elevation_angle'])) * np.cos(np.radians(latitude))))) # 计算太阳直射强度 df['solar_direct_radiation'] = 1353 * 0.7 ** (df['elevation_angle'] / 1000) # 计算受太阳直射辐射的时间段 df['morning_time'] = (df['solar_hour_angle'] > -90) & (df['solar_hour_angle'] < 0) df['afternoon_time'] = (df['solar_hour_angle'] > 0) & (df['solar_hour_angle'] < 90) # 计算每个时间段的太阳直射强度总和 morning_direct_radiation = df[df['morning_time']]['solar_direct_radiation'].sum() afternoon_direct_radiation = df[df['afternoon_time']]['solar_direct_radiation'].sum() # 输出结果 print("在晴天条件下,光伏板受到的太阳直射辐射总能量为:{} W/m2".format(morning_direct_radiation + afternoon_direct_radiation)) print("在晴天条件下,光伏板受到的太阳直射强度上午大于150W/m2、下午大于100W/m2的时长为:{}小时".format(df[df['morning_time']]['hour_angle'].count() + df[df['afternoon_time']]['hour_angle'].count())) # 作图 df.plot(x='hour_angle', y='solar_direct_radiation', figsize=(12,6)) plt.axvline(x=0, color='r') plt.axhline(y=morning_direct_radiation, xmin=0, xmax=0.5, linestyle='--', color='g') plt.axhline(y=afternoon_direct_radiation, xmin=0.5, xmax=1, linestyle='--', color='g') plt.xlabel('Hour Angle (degrees)') plt.ylabel('Solar Direct Radiation (W/m2)') plt.title('Solar Direct Radiation on May 23, 2023') plt.legend(['Solar Direct Radiation', 'Solar Noon', 'Morning Total', 'Afternoon Total']) plt.show() # 寻找最优朝向 optimal_orientation = 0 optimal_inclination_angle = 0 max_direct_radiation = 0 for orientation in range(-90, 91, 10): for inclination_angle in range(0, 91, 10): # 计算太阳直射强度 df['solar_direct_radiation'] = 1353 * 0.7 ** (df['elevation_angle'] / 1000) # 计算受太阳直射辐射的时间段 df['morning_time'] = (df['solar_hour_angle'] > -90) & (df['solar_hour_angle'] < 0) df['afternoon_time'] = (df['solar_hour_angle'] > 0) & (df['solar_hour_angle'] < 90) # 计算每个时间段的太阳直射强度总和 morning_direct_radiation = df[df['morning_time']]['solar_direct_radiation'].sum() afternoon_direct_radiation = df[df['afternoon_time']]['solar_direct_radiation'].sum() # 计算总太阳直射强度 total_direct_radiation = morning_direct_radiation + afternoon_direct_radiation # 判断总太阳直射强度是否最大 if total_direct_radiation > max_direct_radiation: max_direct_radiation = total_direct_radiation optimal_orientation = orientation optimal_inclination_angle = inclination_angle # 输出结果 print("最优朝向为:{},最优水平倾角为:{},最大太阳直射辐射总能量为:{} W/m2".format(optimal_orientation, optimal_inclination_angle, max_direct_radiation)) # 作图 df.plot(x='hour_angle', y='solar_direct_radiation', figsize=(12,6)) plt.axvline(x=0, color='r') plt.axhline(y=morning_direct_radiation, xmin=0, xmax=0.5, linestyle='--', color='g') plt.axhline(y=afternoon_direct_radiation, xmin=0.5, xmax=1, linestyle='--', color='g') plt.xlabel('Hour Angle (degrees)') plt.ylabel('Solar Direct Radiation (W/m2)') plt.title('Solar Direct Radiation on May 23, 2023') plt.legend(['Solar Direct Radiation', 'Solar Noon', 'Morning Total', 'Afternoon Total']) plt.show()
更多内容具体可以看看我的下方名片!里面包含有认证杯一手资料与分析!
另外在赛中,我们也会陪大家一起解析认证杯的一些方向
关注 CS数模 团队,数模不迷路~
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。