太阳能小屋.doc

太阳能小屋的设计

摘要

在能源紧缺的当下，太阳能的开发利用已经是能源领域的热门话题。本文主要研究太阳能小屋安装光伏组件的过程中所涉及到的组件的选配、角度的选择、逆变器的挑选等问题。
针对问题一，本文选用系统抽样方法对数据进行筛选，在贴附式安装的前提下选用光伏组件使小屋的年发电量尽可能大，同时使费用尽量小，为此本文为小屋的顶面、东、西、南四个外表面分别选取A3、B2、C1、A3电池安装，在该种选配下，逆变器的选取分别为SN4、SN4、SN7、SN3四种型号。最后得出总发电量为22361.54千瓦时，经济效益为353540.00元，投资回收年限为31年。
针对问题二，本文首先探究了太阳板的朝向与倾角对光伏电池的工作效率影响，采用Hay的异质分布模型，对一些参数进行量化，着重对屋顶电池构架进行探讨，得出山西省大同地区光伏电池方阵倾角≈46.17度的结论。
针对问题三，本文尽可能多的采用问题一和问题二的结论，综合考虑房屋的地面投影面积的限制，设计出了太阳能小屋（如图1），在小屋表面铺设A3电池效率最高，同种条件下，比问题一中的小屋效率提高了近8%。
由于太阳能发展时间短，电池效率低等缘故未被大众认可，但从本文中仍然可以看出太阳能依然是有着巨大潜力的清洁能源。

关键词：系统抽样 多目标规划 Hay的异质分布模型 太阳能

一、问题重述
在设计太阳能小屋时，需在建筑物外表面（屋顶及外墙）铺设光伏电池，光伏电池组件所产生的直流电需要经过逆变器转换成220V交流电才能供家庭使用，并将剩余电量输入电网。不同种类的光伏电池每峰瓦的价格差别很大，且每峰瓦的实际发电效率或发电量还受诸多因素的影响，如太阳辐射强度、光线入射角、环境、建筑物所处的地理纬度、地区的气候与气象条件、安装部位及方式（贴附或架空）等。因此，在太阳能小屋的设计中，研究光伏电池在小屋外表面的优化铺设是很重要的问题。
对于问题一，根据山西省大同市的气象数据，仅考虑贴附安装方式，选定光伏电池组件，对小屋的部分外表面进行铺设，并根据电池组件分组数量和容量，选配相应的逆变器的容量和数量。
对于问题二，考虑到电池板的朝向与倾角对光伏电池工作效率的影响，需要选择恰当架空方式安装光伏电池，并重新考虑问题1。
对于问题三，根据题目给出的小屋建筑要求，为大同市重新设计了一个小屋，画出小屋的外形图，并对所设计小屋的外表面优化铺设光伏电池，算出铺设及分组连接方式，选配逆变器，计算得出结果。

二、模型假设

1、假设光伏组件在寿命期内不发生故障；
2、假设逆变器将电输入电网过程中不产生额外费用；
3、假设逆变器与电池相连不产生电能损耗。

三、符号说明

：表示倾斜面上的太阳辐射;
：表示水平面上总辐射量；
：表示水平面上直接辐射和总辐射量的比值；
：表示透明系数；
：表示陆地表面的平均反射率；
：表示各面所花费的总费用；
：表示各面电池板全年总发电量；

四、模型的建立及求解
问题一
1、问题分析
在太阳能的开发利用过程中，如何在现有设施条件下对太阳能的充分利用成为了一个无法回避的问题，本问题中，在小屋外表面附贴式铺设光伏组件，由于外表面积有限及阳光因地理位置不同而分布不均，必须对小屋不同的地方做不同的搭配，使得发电量尽可能大，而单位发电量的费用尽可能小，考虑到在多目标规划中，两个目标不可以同时达到，本文考虑将发电量尽可能大作为优先实现的目标，这样就可以先选定光伏组件，再在逆变器的选择中实现成本尽可能低的目标，这种算法的可行解有很多，本文只选取其中一种来探讨。
论文题目中已经给出了小木屋所承受的东西南北四面墙上的总辐射强度，由于顶面与水平面形成夹角且由两块倾斜面组成，总辐射强度未知，所以下面进行顶面总辐射强度的求解。
据《某职工宿舍屋顶光伏系统的设计方案研究》文章中的其倾斜面上所接受到的太阳总辐射量的近似算法得公式：

式中，为倾斜面上的太阳辐射；为水平面上总辐射量；为水平面上直接辐射和总辐射量的比值；、、分别表示倾斜面上的直接辐射、天空散射辐射、地面散射辐射与水平面上的辐射量之间的比值；为地面反射率。由于建筑物是座北朝南，方位角为,可表示为：

其中，
，
式中：、、分别表示赤纬角、纬度角、斜面倾角；、分别表示水平面和倾斜面上各月代表日的日出、日落的时角。
采用Hay的异质分布模型，把水平面上的散射辐射看成是由直接来自太阳方位的绕日分量和来自天空其余部分同质分布散射分量之和，此模型更接近实际的分布。因此，可表示为：

式中，为透明系数。
在小木屋的三维图中，人为规定面积较大的倾斜面的倾斜角为，另一个倾斜面的倾斜角为，根据附图，结合三角函数及反三角函数公式：
，;,
其中，，，，计算得，，查资料发现陆地表面的平均反射率为10%—35%，大同市地处中部无特殊地形地貌，取20%为较准确的值；、。
下面求解小木屋的东西南北四面墙上的平均辐射强度，由于一年当中一个特定地区的光照量存在不确定性，例如，某年大同市空气污染严重，太阳辐射强度减弱，经过改造后，大同市天空清澈，如果采用一年全年的辐射强度记录不仅数据量大，但人为误差也相当大，为此采取系统抽样法抽取25天的纪录作为样本，以此估算全年数据。
系统抽样法以15为步长从总体365天中抽取25天（不够步长多取一天）作为样本得样本为（单位：天）
1 16 31 46 61 76 91 106 121 136 151 166 181 196 211 226 241 256 271 286 301 316 331 346 361.
在样本的25天中，超过的日照时数为213小时，将日总辐射强度低于的值平均，得，可估算得全年超过的日照时数为3110小时；将日总辐射强度低于的值平均，得，日照时数为60小时，可估算得全年日总辐射强度低于的日照时数为876小时。
据上述方法求得样本倾面总辐射强度（见附件三）及各面的平均辐射强度（见表格1）。

东（）
南（）
西（）
北（）
顶面（）
153.46
269.31
227.19
82.65
457.17
表格 1（各面的平均辐射强度）

经比较可知，最大总辐射强度为顶面，依次是南面、西面、东面和北面。
结合光伏组件的特性，本文认为总辐射强度大的应该安装效率较高的光伏组件，比较光伏组件的转换效由高到低排列应是：A3>A1>A2>A4>B2>B1>B3>B5>B6>A6>B7>A5>B4>C1>C5>C3>C2>C4>C11>C10>C8 >C9>C6>C7。
按照接入同一个逆变器的光伏组件电压相差不低于10%的标准将三种电池归为十类（见表格2）。

类别
型号
平均效率
1
C1
6.99%
2
C4
5.84%
3
C3,C5
6.42%
4
C2
6.17%
5
C11
4.27%
6
A1,A2,A3,A6,B2,B5,B6
16.4%
7
A4,A5,B1,B4,B7
15.49%
8
B3
15.98%
9
C6,C8,C9,C10
3.77%
10
C7
3.63%
表格 2（按相差不超过10%的标准分类）

使用多目标整数规划方法将要使全年发电总量最大作为第一目标，将单位发电量的费用尽可能小作为第二目标，则应选用平均效率最大的三组（即6、7、8三组）进行安装。
本文把房屋的五个外露面分开考虑以解决问题一。
1.1 顶面设计
由于顶面中的北倾面根本没有光线直射，而光伏组件最大能量来源为太阳，如果在北倾面铺设电池，效率是相当低的，所以本文不考虑北倾面的电池铺设问题。
据附图数据，南倾面顶面面积，归纳样本，统计分析可知25天中，总辐射能超过的日照时数为213小时，总辐射能低于且高于的日照时数为47小时，根据估算得全年超过的日照时数为3110小时（平均辐射强度为）及全年日总辐射强度低于且高于日照时数为686小时（平均辐射强度为），进一步得到顶面单位面积上的全年总辐射能（见表格3）。
全年日总辐射强度超过的总辐射能（）

全年日总辐射强度低于的总辐射能（）

表格 3（顶面单位面积上的全年总辐射能）

下面逐步优化，在光伏发电量尽可能大的情况下，选择陈本较低且符合实际要求的逆变器，。
先建立方程如下：

式中，为最大发电量，为对应下标数量，为对应下标转换效率，为对应下标的面积。编程求解的如下解：
0.0000，0.0000，47.6422，0.0000，0.0000，0.0000，0.0000 ，0.0000，0.0000，0.0000，0.0000，0.0000，0.0000；
说明满足第一目标时，仅需选取A3电池，此时逆变器SN3、SN4、SN5、SN6均为价格较低的类型，结合逆变器的功率额定输入电压电流的限定，本文认为在上述类型中选择较实惠。对上述类型逆变器与电池进行配置，得到如下三种组合（见表格4）。

型号
费用（元）
组合一
SN5*3+SN3
35100
组合二
SN4*6
41400
组合三
SN8*2+SN5
40800
表格 4（三种组合）

比较可得，组合一最低廉但组合一电池功率超过了逆变器额定功率，同样组合三也超过了，仅有组合二可用，此时总费用为：
；
全年总发电量为：
；
代入数值求得：
；
。
跟据以上两个数据，结合题目给出的电池老化的理想计算方法，得到如下表格（见表格5）。
电池组件图见附件一。
全年总发电量（千瓦时）
成本（元）
35年发电量（千瓦时）
经济效益（元）
投资回收年限（年）
19668.24
169540
619549.56
309774.78
18

表格 5（顶面的数据结果）
1.2 东立面的设计
据附图所示数据，东立面的面积为，由样本数据观察得出在样本25天的273个日照时数中，超过的辐射强度时数为155小时，辐射强度在到的时数为75小时，低于的时数为43小时， 根据估算得全年超过的日照时数为2263小时（平均辐射强度为）及全年日总辐射强度低于且高于的日照时数为1095小时（平均辐射强度为），进一步得到顶面单位面积上的全年总辐射能（见表格6）。

全年日总辐射强度超过的总辐射能（）

全年日总辐射强度低于高于的总辐射能（）

表格 6（顶面单位面积全年总辐射能）
由于单晶硅电池的效率对辐射强度有要求，即低于时，效率将低于原来的5%，而东立面的超过的辐射强度时数所占比例较少，说明用单晶硅电池并不理想。所以选用多晶硅电池及薄膜电池来作优化。为此，建立方程如下：

符号前面已经介绍，此处不再赘述。
编程求解得选用B2电池，计算可知需要14块，但东立面的空间并不规则，实际分析可以算出东立面仅需要9块。
逆变器的选择：满足发电量最大后，对逆变器进行合理的选用以期成本尽量小，为此，可找出条件的逆电器组合（见表格7）。

组别
型号
费用（元）
组合一
SN4*2
13800
组合二
SN6
15000
组合三
SN9
35000
表格 7（逆变器的选用方案）

选SN4，此时总费用为：
；
全年总发电量为：

代入数据，算得：
；
。
全年总发电量（千瓦时）
成本（元）
35年发电量（千瓦时）
经济效益（元）
投资回收年限（年）
1406.23
49800
45702.47
22851.28
——
表格 8（东立面的各项数据）
电池组建见附件一。
1.3 南立面的设计
南立面的可利用面积仅为，由样本数据观察得出在样本25天的273个日照时数中，辐射强度在到的时数为71小时，超过的辐射强度时数为223小时，超过的辐射强度时数所占比例为，由此比例，南立面不用考虑薄膜电池的安装，而南立面正是居家的正门所在，考虑成本问题，结合实际情况，本文认为可参考电池效率和尺寸来选取电池，不给出优化模型，只对所用电池方案进行比较，选出较优方案，根据估算得全年超过的日照时数为3255小时（平均辐射强度为），全年日总辐射强度超过的日照时数为1036小时（平均辐射强度为），低于的时数不考虑，

全年日总辐射强度在到的总辐射能（）

全年日总辐射强度超过的总辐射能（）

表格 9（不同条件下全年日总辐射强度）

按照效率的高低，结合尺寸的大小，本文给出三种方案（见表格10）。
方案
型号
电池费用（元）
方案一
A1*8
25624
方案二
A3*8
23840
方案三
A1*4+A3*4
24732
表格 10（电池方案）

同样，逆变器要考虑费用因素，以达到成本尽可能低的要求，下表给出了三种组合。

组别
型号
费用（元）
组合一
SN3*3
13500
组合二
SN6
15000
组合三
SN8
15300
表格 11（逆电器方案）

从表格十一的费用一栏可清楚得到：组合一为最佳方案，选用SN3，此时最少总费用为：
；此时电池选A3.
全年最大总发电量为：

算得 ，。
全年总发电量（千瓦时）
成本（元）
35年发电量（千瓦时）
经济效益（元）
投资回收年限（年）
71.32
3740
2317.9
1158.95
——
表格 12（南立面的数据）
据上表，南立面的收益报表为亏本。
1.4 北立面的设计
北立面为完全被光，没有太阳直射辐射，仅有天空散射等间接辐射，如果安装电池，将难以使电池工作，为节约成本，本文中不考虑北立面的设计。
1.5西立面的设计
按照上述步骤逐一算出各参数量，，由样本数据观察得出在样本25天的273个日照时数中，超过的辐射强度时数为157小时，低于且高于的日照时数为104小时，根据估算得全年超过的日照时数为2292小时（平均辐射强度为）及全年日总辐射强度低于且高于的日照时数为1518小时（平均辐射强度为），进一步得到顶面单位面积上的全年总辐射能（见表格13）
单位面积全年日总辐射强度超过的总辐射能（）

单位面积全年日总辐射强度超过的总辐射能（）

表格 13（不同条件下单位面积全年日总辐射强度）

根据上面的统计数据，可得出单晶硅电池的工作条件难以满足，如果选用单晶硅电池，将使效率极大降低，且费用成本极大提高，故选用多晶硅电池和薄膜电池较为理想。
建立方程如下：

编程求解得选用C1电池，计算可知需要18块，但西立面的空间同样不规则，实际分析可以算出东立面仅需要15块。当电池成本无法降低时，只能是降低逆变器等其他部件的费用。
逆变器同样按照上述方法选取，得到下表：

组别
型号
费用（元）
组合一
SN7
10200
组合二
SN8
15300
组合三
SN9
35000
表格 14（逆变器的备选方案）
比较费用一栏可清楚得到SN7型号为最佳选择，此时总费用为：
；
全年总发电量为：

代值算得：
；
。
进一步可得出西立面的各项数值（见表15）：

全年总发电量（千瓦时）
成本（元）
35年发电量（千瓦时）
经济效益（元）
投资回收年限（年）
1215.75
17400
39511.88
19755.93
31
表格 15（西立面的各项数据）

综上五个面的成本收益表，按照以上方法计算出的最大发电量为 ：
问题二
2.1太阳板的朝向问题
定义：赤纬角
赤纬角也称为太阳赤纬，即太阳直射纬度，其计算公式近似为

其中为日期序号。令，做关于太阳直射纬度，日期序号和的表格，将数据1月1号至12月31号的日期序号代入上述公式，当时，求解得-23.01，当，求解得-22.93，等等（保留两位有效数字）。
由于数据过于密集且它们之间波动微小，我们取日期序列号步长为5，并不影响结果，做关于上述表格16。

n
T

n
T

1
4.91
-23.01
186
8.09
22.80
6
4.99
-22.54
191
8.18
22.24
11
5.08
-21.90
196
8.26
21.52
16
5.16
-21.10
201
8.35
20.64
21
5.25
-20.14
206
8.43
19.60
26
5.34
-19.03
211
8.52
18.42
31
5.42
-17.78
216
8.01
17.11
36
5.51
-16.40
221
8.69
15.67
41
5.59
-14.90
226
8.78
14.11
46
5.68
-13.29
231
8.87
12,45
51
5.77
-11.58
236
8.95
10.69
56
5.85
-9.78
241
9.04
8.86
61
5.94
-7.91
246
9.12
6.90
66
6.02
-5.99
251
9.21
5.01
71
6.11
-4.02
256
9.30
3.02
76
6.20
-2.02
261
9.38
1.01
81
6.28
-0.00
266
9.47
-1.01
86
6.37
2.02
271
9.55
-3.02
91
6.46
4.02
276
9.64
-5.01
96
6.53
5.99
281
9.73
-6.96
101
6.63
7.91
286
9.80
-8.86
106
6.71
9.78
291
9.90
-10.69
111
6.80
11.58
296
9.98
-12.45
116
6.89
13.29
301
10.07
-14.11
121
6.97
14.90
306
10.16
-15.67
126
7.06
16.40
311
10.24
-17.11
131
7.14
17.78
316
10.33
-18.42
136
7.23
19.03
321
10.41
-19.50
141
7.32
20.04
326
10.50
-20.64
146
7.40
21.10
331
10.59
-21.52
151
7.49
21.90
336
10.67
-22.24
156
7.57
22.54
341
10.76
-22.80
161
7.66
23.01
346
10.84
-23.18
166
7.75
23.31
351
10.93
-23.40
171
7.83
23.43
356
11.02
-23.44
176
7.92
23.40
361
11.10
-23.31
181
8.00
23.18
365
11.17
-23.09
表格 16

随着序号的增大，的值是随周期变化，其值在23.45至-23.45之间变化，所以太阳赤纬角，即太阳直射纬度，一直小于大同的纬度为，也就是说太阳永远不会直射大同。用matlab作图得
横坐标表示以5步长日期序号，纵坐标表示太阳直射纬度。

可得结论：太阳板的朝向一直是向南面的。
2.2太阳板的倾角大小问题
斜面上太阳总辐射量的近似算法得公式：

式中，为倾斜面上的太阳辐射；为水平面上总辐射量；为水平面上直接辐射和总辐射量的比值；、、分别表示倾斜面上的直接辐射、天空散射辐射、地面散射辐射与水平面上的辐射量之间的比值；为地面反射率。 经查资料陆地地表平均反射率为10%-35%，大同地区透明系数，根据大同实际环境和方便计算，所以取20%，
根据附件4，已知水平面总辐射强度=水平面直射辐射强度+水平面散射辐射强度，可得水平面直射辐射强度，由统计学知识，用水平面年直射幅射强度均值表示水平面上直接辐射强度，水平面年总辐射强度均值表示水平面总辐射强度，水平面年散射辐射强度表示水平面散射辐射强度，通过计算，制作表格17。

月份
1
2
3
4
5
6
太阳直射
183.04
186.42
264.35
277.37
299.61
309.58
太阳散射
100.37
150.85
158.29
187.99
166.64
166.09
H总辐射
283.41
337.27
422.64
465.36
466.25
475.67
月份
7
8
9
10
11
12
太阳直射
266.80
264.25
234.25
257.21
196.50
149.33
太阳散射
206.52
173.45
170.52
128.45
101.10
103.42
H总辐射
473.32
437.70
404.77
385.66
297.60
252.75
表格 17

由上表及相关数据得表18.

太阳直射
太阳散射
总辐射
全年
240.73
151.14
391.87
表格 18
其中，，的单位是。因此水平面直接辐射强度为240.73，水平面总辐射强度为391.87，水平面散射辐射强度为151.14。
在方位角不为0的情况下，斜面日射量普遍降低，避免考虑太多，所以我们假设建筑物在北半球理想条件下，座北朝南（正南），方位角为,可表
示为：

其中，
， ，
式中：、、分别表示赤纬角、纬度角、斜面倾角；、分别表示水平面和倾斜面上各月代表日的日出、日落的时角。
采用Hay的异质分布模型，把水平面上的散射辐射看成是由直接来自太阳方位的绕日分量和来自天空其余部分同质分布散射分量之和，此模型更接近实际的分布。因此，可表示为：

式中：为透明系数。
引用定义：时角
时角是以正午12点为0度开始算，每一小时为15度，上午为负下午为正，即10点和14点分别为-30度和30度。因此，时角的计算公式为：

其中为太阳时（单位：小时）。
真太阳时，即太阳时，和平太阳时相差不多，且为减少繁琐计算，用平太阳时代替太阳时（平太阳时是24小时计时式）。根据附件4，分别统计出一年中每个月日出和日落代表时间：

月份
时间
一
二
三
四
五
六
日出
8
8
7
7
6
6
日落
5
6
6
7
7
7
月份
时间
七
八
九
十
十一
十二
日出
6
6
7
7
8
8
日落
7
7
6
5
5
5
表格 19
故的下列结论：

因此：

根据上述公式，对倾角求导得：

令，即可算出山西省大同地区光伏电池方阵倾角≈46.17度，取整数得=46度。
假设可能遮挡物或前排组阵列的高度为H，任意时刻可能遮挡物或前排阵列在地面上的阴影长度为L，阵列与阵列或与遮挡物之间的距离为D。一给定的地方，任意时刻的太阳高度角表示为：

其中为太阳高度角，为时角，为当时的太阳赤纬，为当地的纬度。在北半球遮挡最严重的时间是12月21日至12月23日，一般光伏方阵发电时间在8:00至17:00，在这段时间内系统不能有遮挡。任意时刻遮挡物或前排阵列在地面上的影长为：

则影长在阵列底边上的垂直距离D为：

在大同，度，度，8:00时ω=-60度，17:00时ω=75度。把数据代入上式中求解得，，，由于距离为正数，所以，阵列之间的间距应为，该工程取。
经过分析，无法对小屋四周进行倾角设置，因此对小屋四周不给考虑。对小屋顶层进行分析，由第一问题可知选取A3光伏电池，为最优化选择，根据A3电池的长度1580mm,可以计算出h=1139mm，再根据顶层面积是固定不变的，电池之间可能产生遮挡阴影，影响发电量，电池板之间应保持在h=1139mm的距离，综上考虑，只能在顶层面积架设A3光伏电池最多25块。
对于逆变器的选择，下面给出五种组合

型号
费用（元）
组合一
SN3*7
13500
组合二
SN4*4
27600
组合三
SN4*3+SN3
25200
组合四
SN5*2+SN3
24700
组合五
SN6+SN3
19500
表格 20
选择组合一是最合理的。此时总费用为：

全年总发电量为：

求得： 1387.25
169540元
跟据以上两个数据，考虑电池老化得到下面表格21：
全年总发电量（千瓦时）
成本（元）
35年发电量（千瓦时）
经济效益（元）
投资回收年限（年）
19668.24
169540
619549.56
309774.78
18
表格 21

问题三：根据建筑小屋的要求，附件7，建立目标规划模型，把最低净空高度距地面作为恒定标量，建立函数，

，，，
x代表窗户的面积，代表房屋地面的长，代表房屋地面的宽，代表最短的高，代表最长的高。
根据上式优化方案的，由山西大同地理条件，假设小屋的倾斜角顶为40度，即接近大同的纬度，通过求解的=3.98，
=15，=2.8，=5.4，=31.68。
由第一问，浅显考虑，由于四周铺设光伏电池不经济，仅考虑倾斜面铺设，倾斜面的面积约等于60平方米，再用A3光伏电池，进行铺设，得到最50块。

型号
费用（元）
组合一
SN4*4
27600
组合二
SN4*3+SN3
25200
组合三
SN5*2+SN3
24700
选择组合三是最合理的。
此时总费用为
全年总发电量为
，
41248.26
跟据以上两个数据，考虑电池老化的到下面的表格：
全年总发电量（千瓦时）
成本（元）
35年发电量（千瓦时）
经济效益（元）
投资回收年限（年）
41248.26
25445
1443689.1
721844.5
1.2

五、模型评价
模型的优点
1、 模型直观清晰，可操作性强，可读性强；
2、 采用系统抽样法对数据处理，基本保持了数据的可信度。
模型的缺点
1、 模型在细节处仍不够完善，参数选择有失忖度；

六、参考文献
[1] 方荣生. 太阳能应用技术. 北京：中国农业机械出版社，1985。
[2] 姜启源 谢金星 叶俊，数学模型，高等教育出版社，2011年。
[3] Stine , W. B. and Geyer , M. Power From The Sun.
网址:
[4] http://zh.wikipedia.org/zh-cn/ （维基百科）。
[5] 中华人民共和国国家标准（GB12936.1-91），太阳能热利用术语 第一部分。
[6] 韩中庚 宋明武 邵广纪，数学建模竞赛，科学出版社，2007年。

附件一：

SN4(48V,1.28KW)
B2(45.98V,320W)

B2(45.98V,320W)

B2(45.98V,320W)

B2(45.98V,320W)

B2(45.98V,320W)
B2(45.98V,320W)

SN4(48V,1.6KW)
B2(45.98V,320W)

B2(45.98V,320W)

B2(45.98V,320W)

附件二：matlab程序
b=12.31;
w0=75;
p=0.2;
f=40.1;
for i=1:8760
j=ceil(i/24);
Q=23.45*sin(2*3.14*(284+j)/365);
ws=sin(f)*sin(Q)+cos(f)*cos(Q)*cos(75);
htb=(3.14/180)*ws*sin(27.79)*sin(Q)+cos(27.79)*cos(Q)*sin(ws);
htbo=(75*3.14/180)*sin(40.1)*cos(f)+cos(f)*cos(Q)*sin(75);
rb=htb/htbo;
H=load('q.txt');
H0=load('q0.txt');
H00=load('q00.txt');
rg=H0(i)/H(i);
kb=H00(i)/H(i);
rd=kb*rb+(1-kb)*(1+cos(12.31))/2;
HT=H(i)*[kb*rb+(1+kb)*rd+0.2*rg]
end

附件三：