考虑设备变工况特性的园区综合能源系统两阶段规划优化方法研究

2022-12-09 08:58:18华北电力大学控制与计算机工程学院北京102206
综合智慧能源 2022年10期
关键词:负荷能源阶段

(华北电力大学控制与计算机工程学院,北京 102206)

温港成,石鑫*,张怡,房方

0 引言

随着“双碳”能源发展目标的提出,能源转型加速推进,综合能源系统因其在促进清洁能源消纳、提升能源利用率和降低用能成本方面的显著作用,已成为满足多能供需的新的能源利用实现方式。从地理范围上讲,综合能源系统可以划分为局域级[1-2](如楼宇、社区、工厂等)、区域级[3-4](如城区)和广域级[5-6](如城市、跨区域、跨国等)。园区综合能 源 系 统(Community Integrated Energy System,CIES)[7-8]作为典型的局域级综合能源系统,通过建立电力、天然气、可再生等多种能源之间的耦合利用关系,可以显著提高能源综合利用效率、降低用能成本和减少碳排放量,为实现低碳园区建设提供了一种有效途径。因此,开展CIES规划研究具有重要意义。现有CIES 规划研究主要集中于能源站内设备选型和定容[9],大致可以划分为2个阶段。

第1阶段为单纯的规划模型优化研究阶段。文献[10-12]中以电-气-储综合能源系统作为研究对象,通过建立电网络拓扑、气网络拓扑和储能数学模型,以系统扩建成本最小优化得到能源站容量配置。文献[13]中同时对多个能源站进行规划,充分考虑区域供能的经济性和环保性,优化得到各能源站的容量配置。文献[14]中建立了综合能源系统主要设备的能量流模型,并在此基础上提出基于信息间隙决策理论(IGBT)对系统进行规划优化,计算得到系统最佳规划方案等。以上文献虽对综合能源系统进行规划,但都是只考虑系统的规划优化,未涉及系统运行优化。第2阶段为协同规划优化研究阶段,在进行规划的同时考虑系统运行经济性、环保性等,使得规划的结果更具实用性。文献[15]基于能量梯级利用原理,构建了CIES双层规划优化模型,上层规划求解能源网络布局、能源站的位置和容量,下层通过建立多能源站规划运行联合优化模型对系统运行进行优化,优化结果返回上层进行更新。文献[16-17]考虑区域冷热负荷互补特性,建立了CIES 站-网协同规划优化模型,为供能范围划分和管网布局交替优化提供优化思路。文献[18]对多区域综合能源系统规划进行研究,提出了一种两阶段优化方法,通过优化系统的经济和环境成本求解各区域能源容量配置。上述文献在进行系统规划的同时考虑了系统运行优化,但也存在一定的局限性,如规划模型建模普遍简单,为单目标优化问题;系统供能结构单一,候选设备类型较少;普遍只考虑了系统运行时的经济性和环保性,而对于系统运行时的其他重要因素(如风光出力不确定性)考虑较少等。

针对现有CIES规划研究方面存在的不足,本文充分考虑系统运行时多项因素(经济性、环保性、不确定性),建立了CIES 两阶段优化模型。第1 阶段研究了能源站设备类型及定容问题,第2 阶段研究了系统运行优化问题。第1阶段利用带精英保留策略的非支配排序遗传算法(NSGA-II),分别以系统规划年经济成本和年环境成本最小为优化目标,规划求解系统设备类型及容量。第2阶段主要针对第1 阶段得到的多组Pareto 前沿解,考虑风、光出力不确定性,以系统运行经济、环境总成本、用户舒适度为优化目标,利用Benders 分解算法优化求解各规划方案对应系统运行时设备出力及运行成本。在此基础上针对不同规划方案计算其综合成本,最终实现考虑设备变工况特性的综合能源系统协同规划优化。

1 园区综合能源系统

本文以河北雄安新区某园区冬季供暖期为例进行源-荷协同优化研究。园区所在位置附近建有1 座220 kV 变电站可为其供电,霸保管线是目前新区的主要天然气管线设施,可为园区供气,同时雄安地区具有较好的中深层地热资源可用于园区供热系统,太阳能和风能资源条件中等,可与园区其他能源耦合互补。园区供暖期系统供能结构如图1所示。在供电系统中,除直接从外部购电外,园区可以利用热电联产装置(Combined Heat and Power,CHP)、光伏(Photovoltaic,PV)、风机(Wind Turbin,WT)联合互补运行供电,保障供电稳定性;在供热系统中,电可用于地源热泵(Ground Source Heat Pump,GSHP)实现地热供热,天然气可用于CHP 和燃气锅炉(Gas Boiler,GB)实现燃气供热,通过电、气联供供热保障园区供热的可靠性;同时,通过安装电储(Electric Storage,ES)和热储(Heat Storage,HS)实现用能低谷时存储、用能高峰时释放的功能,提高系统的经济性和灵活性。

图1 园区供暖期系统供能结构Fig.1 Structure of the heating system in the CIES

2 CIES变工况模型

2.1 光伏模型

光伏机组的输出功率和辐照度、温度有关,其最大输出功率为

式中:PPV,t为t时刻机组输出功率;fPV为功率降额因数;PPV,max为装机容量;E0为标准条件下的辐照度;Et为t时刻的辐照度;αPV为功率温度系数;T0为标准条件下光伏机组表面温度;Tt为t时刻光伏机组表面温度。

由式(1)可以看出,光伏出力和光伏表面温度有关,光伏机组的表面温度与环境温度的关系为

式中:Ta为环境温度;Tc,NOCT为标准测试条件下的光伏机组表面温度;Tα,NOCT为标准环境温度;ET,NOCT为温度T下的辐照度;ηs为标准测试条件下的光伏效率;ταs的默认值为0.9。

2.2 风机模型

风机的输出功率和风速有关,当风速小于切入风速时,风机输出功率为0;当风速大于切入风速小于额定风速时,系统采用最大功率追踪控制,追求风能利用效率最大化;当风速大于额定风速时,风机以恒定功率运行;当风速等于额定风速时,存在过渡区间以满足平稳切换,该模型为

式中:PWT,t为t时刻机组输出功率;PWT,max为机组装机容量;vin,vout为切入、切出风速;vR+,vR-为风速过渡区间的上、下界;aWT,i,bWT,i为多项式拟合参数;vw,t为t时刻风速。

2.3 热电联产机组模型

热电联产机组包括燃气发电机和余热回收装置,可以供电和供热,其电转换效率可表示为

热电联产机组的热电比会随着机组工况变化而改变,通过多项式拟合可得到热电比和负载率的函数表达式为

2.4 燃气锅炉模型

燃气锅炉是通过消耗燃气来进行供热,相比于热电联产机组,燃气锅炉具有热转化效率高、响应速度快等特点,燃气锅炉数学模型为

2.5 地源热泵模型

地源热泵从土壤获取低品位热能,通过电力做功输出可用高品位热能。其输出功率与热泵转换效率和负荷率、运行环境温度有关,即

2.6 蓄电设备模型

蓄电设备采用通用储能模型,约束条件主要包括能量平衡约束、存储能量上下限约束以及充放功率约束,即

式中:PES,t为t时刻蓄电池储电量;ηesc和ηesd分别为蓄电池充、放电能效系数;PES,min和PES,max分别为蓄电池储电量上、下限;Pesc,max和Pesd,max分别为蓄电池最大充、放电功率。

2.7 蓄热设备模型

热储能同样采用通用储能模型,约束条件主要包括能量平衡约束、存储能量上下限约束以及充放功率约束,具体为

式中:PHS,t为t时刻热储能装置储热量;ηhsc和ηhsd分别为热储能充、放热能效系数;NHS为热储能装置数量;PHS,min和PHS,max分 别 为 单 组 热 储 能 装 置 储 热 量上、下限;Phsc,max和Phsd,max分别为热储能装置最大充、放电功率。

2.8 柔性负荷模型

在经济学理论中,当商品价格下降时消费者需求会增加。能源市场也同样满足经济学理论,具体体现在能源价格与用能需求的关系上[19]。本文用价格弹性系数来表达用户用能需求与能源单价变化率的关系,可表示为

式中:Δπele,i,πNele,i分别为i时段价格变化量和改变前的价格;Δpj,pNj分别为j时段能源需求变化量和改变前的能源需求;eij为弹性系数,当i=j时为自弹性系数,否则为交叉弹性系数;弹性矩阵E由弹性系数和交叉弹性系数组成,即

2.9 风、光不确定性模型

为应对风电机组、光伏机组出力的不确定性,本文将园区综合能源系统的不确定运行优化问题转化为鲁棒问题描述。采用盒形不确定性集描述风、光出力的不确定性[19],用区间表示风、光发电功率的预测误差,采用不确定性调节参数Γ 来调节模型计算结果的保守度,为

式中:u为0~1 变量,0 表示该时段存在偏差,1 表示该时段不存在偏差;δ为该时段的预测误差;t0为预测初始值,+表示正偏差,-表示负偏差;整数Γ 作为不确定性调节参数,其取值集合为{0,1,…,N},若Γ取值为0 则为确定性优化,若Γ 取值为N,则在调度周期内,风、光出力潜在值为可行区间内任意值。

3 CIES两阶段协同规划优化模型

本文提出的规划优化模型主要分为2 个阶段。第1 阶段为规划阶段,基于园区多能用户用能需求提出设备规划方案,在优化过程中兼顾经济性指标和环境指标,为解决环境成本和经济成本量纲不匹配问题以及目标函数合并不合理导致的计算误差问题,在第1阶段采用多目标优化,将经济优化和环境优化看作2 个独立子问题进行优化求解。第2 阶段为运行模拟阶段,基于第1 阶段得到的规划方案获取设备运行边界条件,考虑风、光出力的不确定性以及用户的需求响应,在典型日场景下进行精细化运行模拟。最后综合方案的经济性指标、环保性指标评价出最佳方案。

3.1 第1阶段规划优化模型

3.1.1 目标函数

第1 阶段规划优化模型为多目标函数优化模型,目标函数包括年经济成本和年环境成本。年经济成本主要由能源站年投资成本和年能源消费成本组成,即

式中:f1为能源站年经济成本;Cinv为能源站年投资成本;Cene为年能源消费成本;下标i,n,d,t分别为设备、设备寿命、典型日和时刻;CRFn为资本回收系数;r为资本年利率;τd为典型日的天数;πele,t和πgas,t分别为典型日t时刻的电价和气价;Pgrid,t为典型日t时刻购买的电功率。

年环境成本主要是年消费能源的污染物排放成本,即

式中:f2为系统年环境成本;下标j,d,t分别为污染物类别、典型日和时刻;θj为第j类污染物排放成本;wgrid,j和wgas,j分别为消耗电和天然气产生的第j类污染物单位排放系数。

3.1.2 约束条件

针对图1 中供能结构,第1 阶段规划优化模型中的约束条件主要为系统能量平衡。园区供暖季典型日电、热能量平衡约束为

式中:Pele,t为电负荷;Ph,t为热负荷。

3.2 第2阶段运行优化模型

3.2.1 目标函数

3.2.1.1 经济性

3.2.1.3 舒适性CIES 在为用户供能时也要注重用户用能体验,不能过渡调控柔性负荷,本文采用负荷曲线变动率来评价用户舒适度,即

3.2.2.2 电价约束

在调度柔性负荷时,电价变化要在一定范围内,否则将会出现过度调控柔性负荷或响应不足等现象,即

此外,设备运行还需满足各类设备数学模型,具体可见第2节。

4 模型求解算法

第1 阶段规划优化模型是多目标优化问题,目标函数分别为最小化年经济成本f1和最小化年环境成本f2,可以直接通过Pareto最优解进行求解而无需转化为单目标优化问题,避免了因目标函数合并不合理(如权重设置不合理等)而引起的计算误差。

NSGA-II 是一种基于Pareto 最优解的多目标优化算法,在迭代过程中采用了快速非支配排序算法,降低了计算复杂度,引入了拥挤度和拥挤度比较算子,可在排序后的同级比较中作为胜负标准,使得Pareto 域中的个体能扩展到整个Pareto 域,保持了种群的多样性,引入了精英策略,扩大了采样空间,对于难以求解的优化问题通常仍能找到较优解且算法具有较好的收敛性,所以本文采用带有精英保留策略的NSGA-II算法进行模型求解。

第1阶段规划模型求解的目的是确定能源站中各设备的配置容量P=Px,maxNx。x表示设备,根据图1中的候选能源设备种类,该解可以表示为长度为9的向量x=[x1,x2,…,x9]。设种群为R,种群大小为K,种群进化代数为g,目标函数为fm(m=1,2),个体拥挤度为nr(r=1,…,K),其计算公式为

式中:fm(i+ 1)为该个体排序后后一位目标函数值;fmaxm和fminm分别为该Pareto 等级下个体目标函数fm的最大值和最小值。

具体算法流程如下。步骤1:系统参数初始化。输入电热负荷、设备参数、辐照度、风速等系统参数。步骤2:初始化种群R,生成种群Rg(g=0)。步骤3:对种群Rg快速非支配排序,划分为不同Pareto等级。步骤4:计算年经济成本f1(g)、年环境成本f2(g),并按照式(24)计算Rg中各个体拥挤度n(g)r。步骤5:对种群Rg进行遗传操作。选择、交叉、变异,产生一代子群Qg。步骤6:将Rg与Qg合成种群Wg,对Wg进行快速非支配排序,划分为不同Pareto 等级,并计算Wg中各个体拥挤度n(g)w。步骤7:按以下规则从种群Wg生成新一代种群Rg+1。(1)根据Pareto等级将整层Wg放入Rg+1,直到某一层个体不能全部放入Rg+1;(2)将该层个体根据拥挤度从大到小排序,依次放入Rg+1中直到填满为止。步骤8:g=g+1,判断是否满足迭代终止条件:如不满足则返回步骤3;满足则输出规划结果x及年经济成本f1、环境成本f2。

第1 阶段规划模型求解通常会得到多组Pareto最优解,即多个规划方案,规划方案是CIES 设备的容量配置结果,为确定最佳规划方案。第2 阶段基于各组解(规划方案)确定CIES 设备的运行边界条件,如设备运行上下限、负载率等,开展系统运行优化分析,在兼顾经济性、环保性和舒适性的基础上考虑风、光出力不确定性,为了便于表达,将优化模型简化为

式中:H为第1 阶段各设备出力计划;U为第1 阶段优化的设备启停状态;Φ为第2阶段最恶劣风、光出力场景集合;O为第2 阶段各设备出力计划,其他符号为目标函数中和约束条件中所用到的常数。

式(25)为鲁棒优化,即在最恶劣的风、光出力场景下找到使系统综合运行成本最小的CIES 调度方案。由于这类优化涉及max-min 问题,导致很难直接求解,通常采用列和约束生成算法或Benders分解算法求解,列和约束生成算法在迭代过程中会不断的增加辅助约束和变量,增加计算压力。

本文采用Benders 分解算法进行求解,将原问题分解为主问题和子问题进行交替求解的形式得到原问题的最优解。对式(25)进行分解,得到主问题的形式为

子问题为max-min 问题,其中内层最小化是线性问题,在此根据强对偶理论[19]及式(25)的对应关系,将内层转化为max 形式,并与外层max 问题合并,得到如下对偶问题

通过对各规划方案进行综合量化评估以确定最佳方案,模型求解算法流程如图2 所示。基于上述分解得到的主子问题,下面给出Benders 分解算法求解流程[20]。

图2 模型求解算法流程Fig.2 Flow of the model solving algorithm

步骤1:求解主问题,满足主问题式(25)约束的初始解,代入子问题优化模型,获得子问题式(28)当前最优O',将初始下限设为无穷小,上限设为无穷大,设置当前迭代次数K=1,设置迭代过程中的允许误差参数ζ。

步骤2:求解第K代主问题式(26),得到二元变量U'的最优决策,并将下界(LB)更新为η。

步骤3:求解子问题式(28),得到最坏情景并更新上界(UB)为v。

步骤4:检查收敛准则UB-LB≤ζ是否满足。如果是,则终止解决方案,输出最糟方案以及设备启停状态。否则,继续迭代。

5 算例分析

本文以河北雄安新区某园区综合能源系统为研究对象,该园区主要包括会议中心、行政办公、生活服务中心、专家公寓和商务酒店共5种不同业态。园区的日负荷可以聚类为3类典型日,即工作日、高峰日和休息日,如图3、图4所示(图中横坐标分别对应00:00,01:00,…,23:00)。

图3 冬季典型日电负荷曲线Fig.3 Typical daily power load curve in winter

图4 冬季典型日热负荷曲线Fig.4 Typical daily heat load curve in winter

模型求解涉及的其他参数,如风速、温度、辐照度如图5—7 所示(图中横坐标分别对应00:00—01:00,01:00—02:00,…,23:00—24:00)。

图5 典型日风速曲线Fig.5 Typical daily wind speed curve in winter

图6 典型日温度曲线Fig.6 Typical daily temperature curve

图7 典型日辐照度Fig.7 Typical daily illumination intensity

各个设备变工况特性如图8 所示,设定调度周期为24 h,单位调度时间为1 h。考虑到该园区可利用资源形式,其冬季供暖期供能结构如图1 所示。系统设备的经济、技术参数见表1。

表1 CIES设备参数Table 1 Parameters of the CIES devices

图8 CIES 设备效率的变工况特性Fig.8 Performances of CIES devices under off-design conditions

5.1 规划结果分析

针对场景1 进行规划,利用带精英保留策略的NSGA-II 算法对第1 阶段规划模型求解,实现年经济成本f1和年环境成本f2的协同优化,算法参数设置如下:种群数量为100,迭代最大次数为100,交叉概率为0.85,变异概率为0.50。得到的Pareto 最优解集(前沿解)如图9 所示。由图9 可知,Pareto 前沿解的分布可以看出碳排放量会随着经济成本的增加逐渐降低,二者之间呈反比关系。图9中根据Pareto最优解集碳排放最小的解和经济成本最低的解可以得到理想情况的最优解和最劣解,基于优劣解距离法(TOPSIS)决策理论[21-22]可知,Pareto 解集中最折中的解方案往往在最优解和最劣解连线的交点,由于第1 阶段为混合整数规划模型,其可行解空间不是连续的,导致连线的交点可能不满足约束条件。同时由于第1 阶段模型是规划模型,主要作用是得出规划方案,并未考虑实际情况中的风、光等不确定因素,所以需要第2 阶段模型对规划方案进行更深层次的分析。在此选取交点附近的3组最优解(即高经济成本低环境成本、中经济成本中环境成本和低经济成本高环境成本)进行分析,对应形成的3种规划方案容量配置情况见表2。

图9 第1阶段规划模型Pareto前沿解Fig.9 Pareto front of the first-stage planning model

表2 第1阶段不同规划方案容量配置情况Table 2 Capacity allocation schemes of different plans in the first stage

由表2可以看出,随着年投资成本的增加,系统会倾向于配置更多的光伏和风电机组从而降低运行成本和碳排放。由于候选设备中同时存在电锅炉(EB)以及GSHP,且EB 的能量转换效率远不如GSHP,故3 个方案都选择以GSHP 供热为主。对比3 个方案容量配置结果可知,3 个方案中CHP,GSHP,ES 和HS 配置容量依次减少而GB 配置容量依次增加,这是由于燃气锅炉造价、运行成本低导致的。虽然CHP 机组也消耗天然气产生碳排放,但其每度电产生的碳排放低于电网的碳排放,同时还会产生余热以供应热负荷的需求,所以系统在倾向于低碳规划时会增加CHP 机组的投资成本。雄安地区光伏、风力资源中等,结合电负荷需求特性曲线,3 个方案均配置了风电机组和光伏机组,3 个方案中电储设备的容量随风电、光伏机组的容量减少而减少,这是由于风、光出力具有间歇性和波动性。系统为了实现安全可靠供能,需要搭配一定容量的储能系统。总体来看,方案1 倾向于低碳排放而高经济成本,方案3 倾向于低经济成本而碳排放增加较大,方案2是二者折中的一种方案。

5.2 运行优化分析

针对第1 阶段得到的3 个规划方案,根据文中所设计的方法,第2 阶段对各方案进行运行优化。本文在运行优化阶段采用盒型不确定集合来描述风光出力的不确定,采用Benders 算法对其进行求解。通过聚类将负荷分为高峰日、工作日和休息日3类,各方案运行成本见表3。结合各个方案各项成本,方案2年综合最低,故方案2为最优方案。

表3 基于冬季典型日负荷的源-荷互动优化调度结果Table 3 Source-load interactive optimal dispatching result based on typical daily power in winter

由表3 可以看出,3 个方案都采用CHP 作为主要的供电、供热方式,其最低供电占比为44.4%。3个方案风、光出力占比逐步降低,方案3倾向于经济性,风电、光伏供电占最低,这导致方案3 的碳排放成本高于其他方案。为降低系统整体运行成本,3个方案都吸收了燃气轮机的余热来供应热负荷需求,余热在总供热量中的占比最低为39.1%,其余热需求大部分来自于热泵,对于尖峰负荷则由燃气锅炉供应。3 个方案的燃气锅炉供热占比逐步增加,这是由于方案3倾向于经济性,在配置阶段燃气锅炉的容量相对较高。相比于方案2,方案1通过增加CHP,ES,HS 的安装容量进一步降低环境成本,提升了系统利用清洁能源满足负荷供应的能力,在表3 中的3 种场景供电比例中也可以看到系统减少了电网购电。以高峰日负荷为基础,场景2 下的优化调度结果如图10 所示(图中横坐标分别对应00:00—01:00,01:00—02:00,…,23:00—24:00)。

图10 工作日供需平衡曲线Fig.10 Supply and demand balance curve

由分析可知,外界环境对清洁能源机组最大输出功率的影响会辐射到整个CIES,间接影响不同时段内各类设备的运行工况。在清洁能源出力的低谷时段(如04:00—05:00,18:00—21:00),电负荷主要由热电联产机组提供,由此产生大量的余热,因此在供热过程中提高热电联产机组供热所占的比重,以减少能量浪费,降低供热成本;而在清洁能源较为丰富的时段(如06:00—09:00,10:00—15:00),电力供应较为充足,热电联产矩阵的负载率处于中等水平,热负荷缺额由燃气锅炉和地源热泵机组共同承担。地源热泵出力受到电价、清洁能源出力综合影响,在电价较高、清洁能源出力较低的时段(16:00—19:00),地源热泵停止工作,提高CHP机组负载率,以满足电负荷、热负荷需求。

5.3 设备变工况特性分析

该部分以工作日为例,对固定工况设备模型和变工况设备模型进行对比分析。采用定工况设备模型调度得到的地源热泵出力结果如图11所示(图中横坐标分别对应00:00—01:00,01:00—02:00,…,23:00—24:00)。由图11可知,固定工况下锅炉效率与实际效率差别较大,当地源热泵处在低负载运行时(如13:00—19:00),优化调度模型高估了地源热泵的出力功率,地源热泵误差功率为负,表示模型计划出力低于实际出力,导致实际调度计划不能满足负荷需求;当地源热泵处在高负载运行时(00:00—11:00),调度模型低估了地源热泵出力功率,导致实际调度计划高于负荷需求,造成了能源浪费。由此可知CIES 设备的负载率会随工况变化出现大幅度波动,固定工况模型对调度结果产生的影响不容忽视。

图11 地源热泵功率曲线Fig.11 Power curve of the GSHP

6 结论

本文对河北雄安新区CIES 规划优化问题进行了研究,充分考虑地区可利用资源形式,设计了园区供能结构,在此基础上提出了一种考虑设备变工况特性及风光出力不确定性的园区综合能源系统两阶段规划优化方法,并提出利用带精英保留策略的NSGA-II 算法进行求解,为运行优化阶段提供可行方案,规划结果更加符合实际情况。

算例分析结果表明,相比于单阶段定参数模型易引起优化后的运行方案与实际运行情况产生偏移,本文所采用的两阶段规划方法可以在第2 阶段考虑系统运行时的不确定性,做到精细化运行模拟,使得规划结果更加符合实际情况。

同时,第2 阶段将柔性电力负荷看作可调度资源,利用分时电价的手段对需求侧进行有效管理,增加可再生能源消纳的同时有效提高系统运行的经济性。

猜你喜欢
负荷能源阶段
关于基础教育阶段实验教学的几点看法
科学与社会(2022年1期)2022-04-19 11:38:42
第六章意外的收获
小学科学(2020年5期)2020-05-25 02:58:24
在学前教育阶段,提前抢跑,只能跑得快一时,却跑不快一生。
莫愁(2019年36期)2019-11-13 20:26:16
用完就没有的能源
————不可再生能源
家教世界(2019年4期)2019-02-26 13:44:20
福能源 缓慢直销路
防止过负荷时距离保护误动新判据
主动降负荷才是正经事
负荷跟踪运行下反应堆一回路控制系统仿真与验证
大热的O2O三个阶段,你在哪?
营销界(2015年22期)2015-02-28 22:05:18
两岸婚恋迈入全新阶段
海峡姐妹(2015年6期)2015-02-27 15:11:19