在“双碳”目标下，新能源快速增长、分布式资源规模化并网已然成为未来配电网发展的趋势．然而，以光伏(photovoltaic，PV)、风电(wind turbine，WT)和电动汽车(electric vehicle，EV)为代表的大规模波动性分布式资源的接入势必影响配电网的安全经济运行．园区综合能源系统(community integrated energy system，CIES)作为一种在用户侧建立的多能耦合系统[1-2]，通过系统优化调度对电/气/热等多种能源形式实现变换、传输和存储，能够提高配电网分布式资源接入的消纳能力[3]．因此，在考虑接入的分布式资源的不确定性下，实现对CIES的优化调度进而提升配电网的能源利用效率，得到了广泛的研究．目前在不确定性场景下CIES的优化调度已有较多的研究．文献[4]提出了一种在风光出力不确定性下的CIES主从博弈模型，有效降低了能源消费者的用能成本．文献[5]考虑了多CIES之间的能量分享和碳交易，基于分布式鲁棒优化方法，采用核密度估计和Wasserstein距离来描述风光不确定集，构建了两阶段分布鲁棒优化模型．文献[6]研究了基于源荷多重不确定性下的CIES鲁棒优化调度方法，采用一种智能优化算法获取最恶劣优化场景，从而求解两阶段鲁棒优化模型．文献[7]采用区间线性优化方法求解考虑分布式电源及电价不确定性的CIES调度模型．文献[8]通过深度强化学习方法构建源荷不确定性变化特征，实现在源荷不确定性环境下的CIES动态调度．文献[9]在风光出力不确定性环境下，研究了基于合作博弈的多CIES分布式优化模型．然而，上述文献所考虑的源荷不确定性环境为单一场景下的不确定性环境，没有考虑不同场景(如不同天气情况)对于风光出力不确定性的影响，从而影响模型中不确定性参数建模的准确性．文献[10-11]研究了在风光出力场景不确定下的CIES多阶段鲁棒优化调度模型．通过加入寻找最恶劣风光出力分布的中间层，完善了CIES鲁棒优化调度模型．文献[12-14]研究了EV接入情况下的CIES鲁棒优化模型．文献[12]在考虑EV接入的情况下，构造了CIES的双层优化模型，通过迭代方法最小化CIES的运行成本和EV充电站的运行成本．文献[13]研究了EV和燃气车接入CIES的两阶段鲁棒优化模型．文献[15]研究了CIES调度层与EV调度层之间的双层调度模型，其中CIES调度层为多目标优化模型，EV调度层以最大化用户满意度为目标．文献[16]基于信息间隙决策理论构建了CIES与EV集群之间的双层鲁棒优化模型，并采用二分迭代法求解优化模型．在前述文献的基础之上，本研究计及EV的接入，建立了多不确定性下CIES两阶段多层鲁棒优化调度模型．采用Stackelberg博弈将CIES与EV优化调度模型表示为双层优化模型，并在模型中考虑EV的数量不确定性，CIES与EV集群通过交互充电电价与充电策略实现优化运行．最后，采用列和约束生成方法(column and constraint generation method，C&CG)迭代求解所提出的两阶段四层鲁棒优化模型．1 计及EV接入的CIES优化运行模型如图1所示，CIES系统电网络部分主要包含PV、WT、储能(energy storage，ES)、电负荷(electrical load，EL)及EV，可以实现与外部配电网的电能交互．CIES可以从外部气网购气，通过天然气网络连接各能量转换设备，包括热电联产机组(combined heat and power，CHP)和燃气锅炉(gas boiler，GB)．碳捕集-电转气系统(carbon capture power plant-P2G，CCPP-P2G)可以捕集CO2并生成天然气．DOI：10.13245/j.hust.250666.F001图1CIES系统拓扑结构图CIES作为Stackelberg博弈的上层领导者，向下层追随者EV集群制定充电电价，EV集群则向上层CIES反馈充电策略．EV集群优化模型的目标是自身效用最大化，即充放电成本最小化，决策变量是自身的充放电行为．CIES优化模型的目标是调度日运行成本最小化，决策变量除了包括系统内能量转换设备的运行约束及不确定性集合外，还包括向下层EV集群制定的充电电价．1.1　CIES优化运行模型考虑多重不确定性，CIES的优化模型可以表示为两阶段四层鲁棒优化问题       minx∈X ∑t=1T(λbuytPbuyt-λselltPsellt)+maxσs∈Ω  maxu∈Uminy∈Y(x,u)∑s=1K∑t=1Tσs(Cgast+Cest+Cctt-Cevt)，(1)式中：σs为风光源荷出力场景概率分布；u为各场景内的不确定变量；y∈(x，u)为第二阶段决策变量；K为场景总数；T为优化时段；λbuyt和λsellt为t时刻CIES的购/售电价；Pbuyt和Psellt为t时刻CIES的购/售电量；Cgast为t时刻的购气成本；Cest为t时刻的储能退化成本；Cctt为t时刻的阶梯碳交易成本[17]；Cevt为t时刻CIES向EV集群的售电收益．具体地Cgast=λgas(QCHPt+QGBt-QP2Gt)Δt；(2)Cest=εes(Pchart+Pdist)Δt；(3)Cctt=εct1Ect1+εct2Ect2+εct3Ect3；(4)Cevt=Nevσevλevt(Pev,chart-Pev,dist)，(5)式中：λgas为CIES的购气价格；QCHPt，QGBt和QP2Gt分别为t时刻CHP和GB的天然气耗气量及P2G的天然气生成量；εes为ES的退化成本系数；Pchart和Pdist为t时刻ES的充放电功率；εct1~εct3为碳排放购额区间1~3的购买交易成本系数；Ect1~Ect3为区间1~3的碳排放购额；Nev为CIES充电桩平均每日充放电的EV数量；由于EV数量庞大，难以针对每一辆EV构建不确定集合，因此采用分布鲁棒优化方法，通过聚类方法将离散EV数据聚合为Kev类；σev为Kev类EV的数量概率分布[18]；λevt为CIES向EV集群制定的充电电价．式(1)约束条件包含以下几种约束．a．CIES购售电约束，0≤Pbuyt, Psellt≤Pmaxt，(6)式中Pmaxt为CIES的购售电量限值．b．CHP和GB机组运行约束，0≤PCHP,et,PCHP,ht≤PCHP,e,maxt,PCHP,h,maxt;0≤PGB,ht≤PGB,h,maxt,(7)式中：PCHP,et和PCHP,ht为CHP机组的发电和发热功率；PCHP,e,maxt和PCHP,h,maxt为发电和发热功率限值；PGB,ht为GB机组的发热功率；PGB,h,maxt为发热功率限值．c．CCPP-P2G机组运行约束，0≤PCCPPt≤PCCPP,maxt;0≤QCCPPt≤QCCPP,maxt;0≤PP2Gt≤PP2G,maxt;QCCPPt=ξCCPPPCCPPt;QP2Gt=ξP2GPP2Gt,(8)式中：PCCPPt和QCCPPt为CCPP机组耗电功率和捕碳量；PCCPP,maxt和QCCPP,maxt为功率限值和捕碳量限值；PP2Gt和PP2G,maxt为P2G机组耗电功率和限值；ξCCPP和ξP2G为CCPP-P2G机组的转化系数．d．储能系统运行约束，0≤Pchart,Pdist≤Pchar,max,Pdis,max;PchartPdist=0;Smin≤St≤Smax;St=St-1+(ηcharPchart-ηdisPdist)Δt,(9)式中：Pchart和Pdist为储能在t时刻的充放电功率；St为储能t时刻的荷电状态(SoC)．e．系统功率平衡约束，   Pbuyt+PCHP,et+PPVt+PWTt+Pdist+NevσevPev,dist= Ploadt+Pchart+Psellt+PCCPPt+PP2Gt+NevσevPev,chartPGB,ht+PCHP,ht=Hloadt,(10)式(10)中第一个等式为CIES的电功率平衡约束，第二个等式为CIES的热功率平衡约束．f．阶梯碳交易约束，Ece-Ecq=Ect1+Ect2+Ect3;0≤Ect1≤Ect,max1;0≤Ect2≤Ect,max2;Ect3≥0,(11)式中：Ece和Ecq为实际碳排放量和碳排放配额；Ect,max1和Ect,max2为阶梯碳交易区间限值．g．EV充电电价约束，∑tλevt/T≤∑tλbuyt/T;λev,mint≤λevt≤λev,maxt,(12)式中λev,maxt和λev,mint为充电电价的上下限值．h．EV种类概率分布σev约束，Ωev=σev≥0;∑sσev=1;σev-σev01≤θ1;σev-σev0∞≤θ∞,(13)式中：σev满足分布鲁棒优化问题的不确定集合Ωev约束[19]；σev0为初始种类概率分布；θ1和θ∞分别为1-范数和∞-范数的概率分布置信度，计算方式参考文献[20]．上述决策变量属于决策变量集合X和Y．集合U包含PV，WT出力及电热负荷不确定性，本研究采用盒式不确定集在第二阶段优化问题中来表示单一场景下源荷不确定性，UCIES=ξt=[PPVt,PWTt,Ploadt,Hloadt];ξt=ξ˜t+u+tΔξ-u-tΔξ;u+t+u-t≤1,∑tu+t+u-t≤τ,(14)式中：ξ˜t为不确定预测值；u+t和u-t为二元基准变量；Δξ为不确定波动范围；τ用来限制不确定性波动从而调节优化模型保守度．考虑不同场景下不确定性数据的不同特性，在传统的两阶段鲁棒优化min-max-min模型中加入场景分布不确定性，Ω=σs≥0,N=1,2,...,Ks;∑sσs=1;σs-σs01≤θ1;σs-σs0∞≤θ∞,(15)式中：σS为场景概率分布；σS0为初始场景概率分布；θ1为1-范数置信度约束；θ∞为∞-范数置信度约束Ω为概率分布集合；Ω表示方式与式(13)类似，同样满足1-范数和∞-范数的概率分布置信度约束．1.2　EV集群优化模型CIES作为Stackelberg博弈上层的领导者，向下层EV集群制定合理的充电电价从而引导EV的充电策略．考虑EV种类数量不确定性，EV集群的优化模型可以表示为     {Pevt}=argmax∑t=1TNevσev[λutt(ηev,charPev,chart-ηev,disPev,dist)-λevt(Pev,chart-Pev,dist),(16)式中：Pevt为EV集群的优化充电策略；λutt为单位效用系数[18]．式(16)表示EV集群作为博弈下层，优化目标为调度日内使得自身效用最大化的充放电计划．式(16)的约束条件为0≤Pchar,nt,Pdis,nt≤Pchar,n,max,Pdis,n,max;Pchar,ntPdis,nt=0;Sn,min≤Snt≤Sn,max;Snt=Snt-1+(ηchar,nPchar,nt-ηdis,nPdis,nt)Δt;Snt,end=Sn,obj,(17)式中：Pchar,nt和Pdis,nt为第n类EV在t时刻的充放电功率；Pchar,n,maxt和Pdis,n,maxt为第n类EV在t时刻的充放电功率限值；Sn,max和Sn,min为第n类EV的SoC限值；Snt,end和Sn,obj为第n类车离开充电桩时间的SoC和期望的SoC．式(1)和(16)包含了决策线性项乘积所构成的非线性项，可以通过KKT条件(卡罗需-库恩-塔克条件)和对偶理论处理双层优化问题及所包含的非线性项[21-22]．最终的优化结果即为博弈的均衡点．1.3　CIES与EV的博弈模型转化式(1)的紧凑形式可以表示为：minx∈XcTx+maxσs∈Ω  maxu∈Uminy∈Y(x,u)∑σs(dTy-bTz) ；(18)s.t.    X=Ax≤b,Y=Cx+Dy+Fz+Eu≤g,U=Hu≤k,Ω=Mσs≤n,(19)式中：c为第一阶段参数矩阵；x为第一阶段决策变量矩阵；d为第二阶段参数矩阵；y为第二阶段决策变量矩阵；b为电动汽车参数矩阵；z为电动汽车决策变量矩阵；A为第一阶段决策变量约束参数矩阵；C，D，F和E为第二阶段决策变量约束参数矩阵，g为其约束上限；H为不确定决策变量约束参数矩阵；u为不确定参数，k为其约束上限；M为场景概率分布约束参数矩阵；n为其约束上限．式(16)的紧凑形式可以表示为z=argmax(aTz-bTz);(20)s.t.   Gz≥h:μ,(21)式中：μ为式(21)的对偶变量；bTz为非线性项．式(20)和(21)的KKT条件可以表示为：aT-bT+μG=0;0≤(Gz-h)⊥μ≥0.(22)根据对偶理论，最优解处原问题和对偶问题的目标函数值相等，因此原目标函数可以通过对偶变量表示为bTz-aTz=μTh．(23)因此，在模型主问题中的非线性项可以由式(23)来替代．式(18)的约束条件除了包含原约束条件(19)外，还包括了KKT约束条件(22)．式(18)的最优解为CIES和EV集群博弈问题的均衡点．2 CIES多层鲁棒优化求解策略2.1　优化模型主问题主问题的物理意义：在子问题中找到最恶劣场景分布及源荷最恶劣出力和EV集群种类分布后，来优化CIES的调度成本．紧凑形式可以表示为minx,yn,η,λev      cTx+η；(24)s.t.    η≥∑sσs*dTyn+aTz+μTh,Cx+Dyn+Fz+Eu≤gn,∀n∈Niter,Gλev≤e,aT-bT+μG=0,0≤(Gz-h)⊥μ≥0.(25)式中：η为子问题的最优割；n和Niter为当前迭代次数和总迭代次数．主问题优化结果作为迭代上界UB．2.2　优化模型子问题子问题的物理意义为：找到最恶劣场景分布以及源荷最恶劣出力和EV集群种类分布后，返回到主问题．紧凑形式可以表示为：maxσs∈Ω  maxu∈Uminy∈Y(x,u)∑σs(dTy-bTz)；(26)s.t.    Cx*+Dy+Fz+Eu≤g,Hu≤k,Mσs≤n.(27)子问题优化模型(26)可以进一步分解为两个子问题去求解[10]：σs固定子问题和u固定子问题．σs固定子问题可以表示为：maxu∈Uminy∈Y(x,u)∑σs*(dTy-bTz);(28)s.t.    Cx*+Dy+Fz+Eu≤g:λ,Hu≤k,(29)式中σs*为定值．式(28)的双层优化问题同样可以通过KKT条件转化为单层优化问题：maxu∈U,y∈Y(x,u)∑σs*(dTy-bTz)；(30)s.t.    ∑σs*dT+λTD=0,λ≥0,0≥λ⊥Cx*+Dy+Fz+Eu-g≤0,Cx*+Dy+Fz+Eu≤g:λ,Hu≤k.(31)u固定子问题可以表示为：maxσs∈Ωminy∈Y(x,u)∑σs(dTy-bTz);(32)s.t.    Cx*+Dy+Fz+Eu*≤g,Mσs≤n,(33)式中u*为定值．式(32)的外层优化变量场景概率分布与内层优化变量没有耦合关系，因此无须KKT条件将(32)转化为单层优化模型求解．采用文献[10]的改进C&CG算法，当σs固定子问题的优化结果UB1,in和u固定子问题的优化结果UB2,in相等时，子问题停止迭代求解，输出子问题优化结果作为上界UB．3 仿真分析CIES仿真系统结构如图1所示．在本算例中，分布鲁棒优化的场景数聚合为5类，电动汽车类型同样也聚合为5类，设置CIES每天服务EV数量为50．仿真基于MATLAB R2021a平台并调用求解器GUROBI 10.0.1实现．CIES源荷数据的鲁棒调度优化结果如图2所示，图中虚线表示各个场景下的源荷数据的鲁棒优化结果．从图2可以看出：风光出力的最恶劣情况发生在电价峰时段和负荷峰时段，此时电热负荷需求量大且配电网电价较高．DOI：10.13245/j.hust.250666.F002图2CIES源荷数据的鲁棒调度优化结果3.1　结果分析C&CG算法求解CIES两阶段鲁棒优化问题的迭代次数与迭代结果如表1所示，收敛精度指标ε设置为1×104．C&CG算法良好的收敛性能，算法迭代流程在第三次完成收敛，从而证明了算法的有效性．DOI：10.13245/j.hust.250666.T001表1C&CG算法迭代次数与结果迭代次数迭代上界/万元迭代下界/万元优化间隙/%12.848 0722.627 6110.077 422.848 6022.847 2330.000 532.849 3542.849 4550.000 0图3所示是以场景1为例的CIES电系统和热系统的鲁棒调度优化结果．可以看出：在配电网电价谷时段24:00—07:00间，CIES为满足系统电负荷需求主要从电网购电，此时段为了降低系统运行DOI：10.13245/j.hust.250666.F003图3场景1下电系统和热系统的鲁棒调度优化结果成本，CHP机组的电出力为0，同时储能和EV集群主要以充电行为为主．在配电网电价峰时段08:00—23:00间则相反．因此，CIES通过制定EV充电电价从而引导EV集群进行充放电行为，相对提高了CIES的运行经济性与灵活性，但同时又由于EV集群的停靠时间与自身效用的因素，在个别的CIES负荷峰时段仍须充电从而满足自身用能需求．在CIES热系统的优化结果中可以看出：热负荷需求在场景1下全部由CHP机组的热出力满足，而GB机组的热出力为0，这是由于在CHP机组的热出力范围内，其发热成本低于GB机组．在热负荷高于CHP机组热出力范围的时段，GB机组出力．EV集群的充放电曲线图如图4所示，其中每条曲线为某一类型的一辆EV的充放电曲线．从图4可以看出：EV集群在停靠CIES充电站时间段内，负荷的峰谷时段参与CIES的充电策略引导，其中在CIES制定的较高电价下表现出放电行为从而满足CIES内部的电负荷用能需求，在CIES制定的较低电价下表现出充电行为来满足自身的用能需求．DOI：10.13245/j.hust.250666.F004图4场景1下EV集群充放电曲线图5所示为CIES运行场景概率的分布鲁棒优化结果．场景2和场景5为相对恶劣的运行场景，CIES的调度成本更高，因此其场景发生概率在最恶劣分布下增长．图5也证明了所提CIES两阶段鲁棒优化模型具有更强的鲁棒性，避免了单一场景过度乐观或过度保守的情况．DOI：10.13245/j.hust.250666.F005图5场景概率分布鲁棒优化结果3.2　算例对比分析设立4个系统运行算例对比分析场景概率分布、EV集群种类概率分布及源荷不确定性对CIES调度运行成本的影响．算例1：不考虑单独场景下的风光出力不确定性以及电热负荷不确定性，仅考虑场景概率分布不确定性的CIES运行情况．算例2：考虑单独场景下的风光出力不确定性及电热负荷不确定性，不考虑场景概率分布不确定性．算例3：考虑CIES的场景不确定性与风光出力不确定性，但不考虑EV集群的种类不确定性．算例4：最恶劣运行场景发生时不采用提出的最恶劣场景下的调度方案，仍然考虑风光出力不确定性等．算例5：所提的优化调度场景．5个场景的成本分别为2.633 870，2.841 475，2.846 873，2.859 711和2.849 455万元．算例5的CIES运行成本高于前3个运行算例，证明了所提出的优化模型在牺牲运行经济效益同时提高了模型的鲁棒性．而在算例4中，当最恶劣运行场景发生时不采用该场景下的优化调度方案，则其调度成本较高．在算例2中，单独场景下的最低运行成本和最高运行成本分别为2.723 024万元和2.871 035万元．所提的两阶段鲁棒优化模型避免了运行成本过度乐观或过度保守的情况．3 结语考虑了运行场景概率分布的不确定性，进一步提升了系统运行的鲁棒性．采用的分布鲁棒优化方法避免了单一场景下系统优化结果的相对乐观性和保守性．通过分解CIES的四层鲁棒优化问题为主问题和子问题，保证算法流程的快速收敛．