In recent years, energy consumption and global environmental problems have become increasingly serious, and the energy transition toward a low-carbon energy system is highly desirable (Zhao et al., 2017). The energy internet provides great potential for reliable power supply and an improved energy efficiency by integrating various energy production/conversion components, including renewable energy sources (RESs), energy storage systems (ESSs), hydrogen fuel vehicles (HVs), and electric vehicles (EVs) (Strasser et al., 2015; Meng and Wang, 2017; Yang and Fang, 2017). At the same time, interactions between urban electrical, transportation, and natural gas systems have become increasingly common due to the utilization of coupling components such as microturbines, combined heat and power (CHP), and power to gas units (Mancarella, 2014; Shabanpour-Haghighi and Seifi, 2015). The coupling of urban multi-energy networks positively contributes to the improved energy utilization efficiency, multi-energy complement, and the construction of a low-carbon sustainable energy system (Xu et al., 2020).

Traditionally, coupled energy systems have been operated and optimized independently (Chen et al., 2014). The urban distribution network that is inactive generally has limited ability toregulate power flow, which results in the problem of a large difference between the daily peak and valley of the net load (Uddin et al., 2018). On the other hand, the traditional energy storage has relatively small storage capacity with a high investment cost, and hence, it is difficult to store energy effectively on a large scale and for a long time (Zhao et al., 2020). Considering the fact that the traditional urban distribution network cannot fully meet the electricity demand of new-type loads such as EVs and HVs, the development of day-ahead co-optimization of multi-energy systems is of considerable interest.

Owing to the aforementioned considerations, the present work develops a day-ahead co-optimization strategy for urban energy networks considering traffic flows.

To fill the research gaps summarized above, our model makes the following contributions:

1) We propose a mixed-UE traffic flow model based on an urban transportation system and divide the overall vehicular flows into five categories: GV traffic flow, EV traffic flow, HV traffic flow, HV charging flow, and EV charging flow. Our model, hence, accurately simulates the distribution of EVs and HVs in the urban transportation network.

The remainder of this article is organized as follows: Section 2 elaborates on the mathematical model of day-ahead operation for an urban energy system considering traffic flows. Section 3 presents the performance of the proposed optimization model. Finally, Section 4 presents the conclusions of the work.

Since AC power flow constraints are non-convex, the day-ahead dispatch of the urban energy system is a non-convex optimization problem, which is difficult to solve. The second-order cone (SOC) relaxation based on the branch flow model is implemented to improve its computational tractability. Therefore, we use the branch flow model to replace the nonlinear power flow model of the urban distribution network (Baran and Wu, 1989): U ̃ i , t = U i , t 2 , ∀ i ∈ Ω n o d e , ∀ t ∈ T , (1) I ̃ i j , t = I i j , t 2 , ∀ i j ∈ Ω l i n e , ∀ t ∈ T , (2) ∑ j ∈ M i P i j , t − ∑ k ∈ N i P k i , t − I ̃ k i , t R k i , t = P i , t i n j ∑ j ∈ M i Q i j , t − ∑ k ∈ N i Q k i , t − I ̃ k i , t X k i , t = Q i , t i n j ∀ i ∈ Ω n o d e , ∀ t ∈ T , (3) U ̃ j , t = U ̃ i , t − 2 P i j , t R i j , t + Q i j , t X i j , t + I ̃ i j , t R i j 2 + X i j 2 ∀ i j ∈ Ω n o d e ∀ i j ∈ Ω l i n e , ∀ t ∈ T (4) I ̃ i j , t U ̃ i , t = P i j , t 2 + Q i j , t 2 ∀ i ∈ Ω n o d e ∀ i j ∈ Ω l i n e , ∀ t ∈ T , (5) P i , t i n j = P i , t s u b + P i , t D G − P i , t L + P i , t d i s − P i , t c h − P i , t E H − P i , t G F − P i , t C H P − P i , t P 2 H , ∀ i ∈ Ω n o d e , ∀ t ∈ T , (6) Q i , t i n j = Q i , t s u b + Q i , t D G − Q i , t L ∀ i ∈ Ω n o d e , ∀ t ∈ T . (7)

Eqs 1, 2 replace the voltage and current quadratic terms in the original nonlinear power flow equation with new variables. The DistFlow model for radial electricity networks is formulated by 3–7, in which variables of the squared voltage magnitude and squared branch current are employed. Eqs 3–5 represent the active and reactive power flow balance. Eqs 6, 7 represent the injected active and reactive power of the node.

The operation of the urban AC distribution network should satisfy the following security constraints: P i min s u b ≤ P i , t s u b ≤ P i max s u b Q i min s u b ≤ Q i , t s u b ≤ Q i max s u b ∀ i ∈ Ω s u b , ∀ t ∈ T , (8) P i j , t 2 + Q i j , t 2 ≤ S i j max 2 ∀ i j ∈ Ω l i n e , ∀ t ∈ T , (9) P i , t D G = P i max D G ∀ i ∈ Ω D G N , ∀ t ∈ T , (10) Q i , t D G = P i , t D G ⁡ tan ⁡ φ ∀ i ∈ Ω D G N , ∀ t ∈ T , (11) U ̃ i min ≤ U ̃ i , t ≤ U ̃ i max ∀ i ∈ Ω n o d e , ∀ t ∈ T . (12)

Constraint (8) limits the transmission power capacity of the root node. Constraint (9) represents the feeder capacity constraint. Then, curtailment constraints with respect to distributed generation are given in 10, 11. Constraint (12) limits the operating range of the squared voltage magnitude.

ESS operation constraints include charge and discharge state constraints, storage capacity constraints, charge and discharge power constraints, and daily allowable charge and discharge times constraints (Gabash and Li, 2012), which are provided as follows: y i , t c h + y i , t d i s ≤ 1 , ∀ i ∈ Ω E S S , ∀ t ∈ T , (13) 0 ≤ P i , t d i s ≤ y i , t d i s P max e 0 ≤ P i , t c h ≤ y i , t c h P max e , ∀ i ∈ Ω E S S , ∀ t ∈ T , (14) E i , t + 1 E S S = E i , t E S S + α P i , t c h − β P i , t d i s 0.2 E i max ≤ E i , t E S S ≤ 0.9 E i max , ∀ i ∈ Ω E S S , ∀ t ∈ T , (15) ∑ t y i , t + 1 c h − y i , t c h ≤ λ max E S S ∑ t y i , t + 1 d i s − y i , t d i s ≤ λ max E S S , ∀ i ∈ Ω E S S , ∀ t ∈ T . (16)

Constraint (13) denotes the charging and discharging states of ESS, and constraint (14) represents upper limit of charge and discharge. Eqs 15, 16 are energy constraints of ESS and daily allowable charge and discharge times constraints, respectively. Then, the quadratic equality (5) is further relaxed into the following SOC inequality (Taylor and Hover, 2012; Farivar and Low, 2013): 2 P i j , t 2 Q i j , t I ̃ i j , t − U ̃ i , t 2 ≤ I ̃ i j , t + U ̃ i , t , ∀ i ∈ Ω n o d e ∀ i j ∈ Ω l i n e , ∀ t ∈ T , (17)

where . 2 is the mathematical expression of the Euclid norm.

The EH (Chen et al., 2019) depicted in Figure 1A represents the coupling between electricity and gas networks. CHP represents the key component of this EH, which generates electricity and heat simultaneously, and the urban distribution network provides the charging service for EVs through the EH. The energy conversion efficiency between the input and output energy flows of the EH is formulated in 18 L i , t e L i , t h = η e e v η C H P , e 0 v η C H P , h + 1 − v η g h P i , t E H F i , t E H . (18)

Eq. 19 represents the dispatch factor of the EH input gas flow 0 ≤ v ≤ 1 . (19)

Then, the EH also needs to satisfy the following operating constraints: F E H , min ≤ F i , t E H ≤ F E H , max , (20) S c h p , min ≤ v F i , t E H ≤ S c h p , max , (21) S g f , min ≤ 1 − v F i , t E H ≤ S g f , max . (22)

Constraint (20) gives the range of total gas purchases of EH. Constraints (20) and (21) represent the production capacity of CHP units and gas furnaces, respectively.

For a specific case, the nodes of the distribution network are connected to some components of EH, for example, only CHP units or gas furnaces, as shown in Figure 2. On the other hand, the power distribution network converts electric power into hydrogen through P2H devices to provide hydrogen supply services for hydrogen-fueled vehicles on the traffic side (Korpas and Holen, 2006; Pan et al., 2021). The structure of P2H is shown in Figure 3.

The energy conversion efficiency between the input and output energy flows of P2H is given in 23 and 24 H e l , t = η e 2 h y P i , t P 2 H Δ t / H c , (23) 0 ≤ H e l , t ≤ H e l max . (24)

Constraint (23) calculates the energy conversion efficiency between the input and output energy flows of P2H. Constraint (24) limits the hydrogen production.

The user equilibrium (UE) model has been widely used to describe the urban transportation system in transportation research. In the case of the given traffic network and traffic demand, the UE model can output the current traffic state, namely, the traffic flow distribution on the road. The transportation system model in this article is based on the mixed-UE traffic flow model (Wei et al., 2018). We divide the overall vehicular flows into five categories: GV traffic flow, EV traffic flow, HV traffic flow, HV charging flow, and EV charging flow. The difference between traffic flow and charging flow lies in that the former represents the route choice of EVs and HVs, while the latter represents the charging choice of EVs and HVs. Traffic link has regular link, charging link, and bypass link. The charging links correspond to EV and HV queuing and charging events in EVCSs and HVCSs, and bypass links denote bypass events at EVCSs and HVCSs (Lv et al., 2021b; Teng et al., 2021).

The improved mixed-UE transportation model is a nonlinear complementary problem (NCP) consisting of four constraint sets U E − N C P C o n s − F l o w , C o n s − T i m e C o n s − C o s t , C o n s − C P (25)

where constraint sets are presented as follows:

Cons − Flow: ∑ ρ ∈ λ g r s f ρ , r s , t g = π g q r s , t , (26) ∑ ρ ∈ λ e r s f ρ , r s , t e = π e q r s , t , (27) ∑ ρ ∈ λ h r s f ρ , r s , t h = π h q r s , t , (28) x a , t g = ∑ r s ∑ ρ ∈ λ g r s f ρ , r s , t g δ a , ρ , r s g + ∑ r s ∑ ρ ∈ λ e r s f ρ , r s , t e δ a , ρ , r s e + ∑ r s ∑ ρ ∈ λ h r s f ρ , r s , t h δ a , ρ , r s h , (29) x a , t e = ∑ r s ∑ ρ ∈ λ e r s f ρ , r s , t e δ a , ρ , r s e , (30) x a , t h = ∑ r s ∑ ρ ∈ λ h r s f ρ , r s , t h δ a , ρ , r s h , (31) π g + π e + π h = 1 , π g , π e , π h > 0 . (32)

Cons − Time: t a , t r g = t a 0 1 + 0.15 x a , t g c a r g 4 , (33) t a , t c h , e = E e p s e r , e + t a , c h , e max x a , t e c a c h , e 3 , x a , t e ≤ c a c h , e , (34) t a , t c h , h = E h p s e r , h + t a , c h , h max x a , t h c a c h , h 3 , x a , t h ≤ c a c h , h . (35)

Cons − Cost: c ρ , r s , t g = ∑ a ω t a , t r g δ a , ρ , r s g , (36) c ρ , r s , t e = ∑ a ω t a , t r g δ a , ρ , r s g + ∑ a ω t a , t c h , e + λ a , t e E e δ a , ρ , r s e , (37) c ρ , r s , t h = ∑ a ω t a , t r g δ a , ρ , r s g + ∑ a ω t a , t c h , h + λ a , t h E h δ a , ρ , r s h . (38)

Cons − Cp: 0 ≤ f ρ , r s , t g ⊥ c ρ , r s , t g − u r s , t g ≥ 0 , (39) 0 ≤ f ρ , r s , t e ⊥ c ρ , r s , t e − u r s , t e ≥ 0 , (40) 0 ≤ f ρ , r s , t h ⊥ c ρ , r s , t h − u r s , t h ≥ 0 . (41)

Constraint set Cons − Flow describes flow conservation in the TN. Specifically, constraints (29)–(31) indicate that the traffic flow on each link is equal to the sum of the flows on all paths that pass through the link. Constraints (26)–(28) explain that the traffic flow on all paths connecting each O-D pair is equal to the corresponding traffic demand. Constraint set Cons − Time gives link travel time and estimation of different link types, respectively. For regular link a ∈ T A r g , constraint (33) uses the Bureau of Public Roads (BPR) function [Bureau of Public Roads (1, 1964)] to represent the link travel time as a function of x a g . In constraint (7), the first and second terms denote the charging time and queuing time, respectively. Constraint set Cons − Cost describes the travel costs of GVs, EVs, and HVs that choose path ρ. Note that the cost for EVs and HVs contains both travel cost and charging cost, respectively. Constraint set Cons − Cp represents the complementarity condition of the Wardrop UE principle (Sheffi, 1985), which indicates that no GV/EV/HV travelers could decrease their travel costs by unilaterally switching their route choices.

By changing the electricity price and hydrogen price for each period of time, we calculate the transportation network model 24 times. Then, we can get the load curve of each charging station and hydrogen charging station in a day according to Eqs 42, 43. In addition, the heat load L i , t h can be obtained from the literature (Soroudi, 2017) L i , t e = E e x a , t e , (42) L i , t h = E h x a , t h . (43)

With the goal of minimizing the sum of gas purchase cost and electricity purchase cost, we established a day-ahead optimization dispatch model of an urban energy system, and the optimization model is summarized below (Lu et al., 2020; Bakeer and Hossam, 2021; Cheng et al., 2021; Lu et al., 2020; Bakeer and Hossam, 2021; Cheng et al., 2021).

The objective function is provided as min μ e , t P i , t s u b + μ g , t F i , t E H + F i , t C H P + F i , t G F . (44)

The objective function (44) is subject to Day-ahead operation constraints: ( 1 ) − ( 4 ) ( 6 ) − ( 17 ) ( 18 ) − ( 24 ) ( 42 ) − ( 43 )

Figure 4 shows the framework of the day-ahead optimization dispatch model of an urban energy system, which clarifies the interactions between power and transportation networks from the modeling perspective.

Figure 6 shows the daily heat load, the charging load of EVs, and the hydrogen load of HVs on the traffic side. As can be observed from this figure, the charging load curve of EVs has the same trend with the hydrogen load curve of HVs, and both have the characteristic of double peaks. During the morning rush hour from 8:00 to 10:00 and the evening rush hour from 17:00 to 19:00, the electric load and hydrogen load on the traffic side are relatively high, while the heat load is higher during the day and lower at night (Ye et al., 2021). The trend of these three loads is consistent with that of the total load.

Figure 7 shows the comparison of the wind power output curve, total original load curve, and total net load curve without considering energy storage. In this figure, load 0 represents the total original load curve, and load 1 represents the total net load curve without considering the energy storage. As can be observed from the figure, wind turbines produce more output at night, which results in a lower net load at the valley time. However, the output of wind turbines is less in the daytime, and hence, it cannot significantly reduce the peak load. Consequently, the daily wind power production results in a larger peak-to-valley difference of net load.

It is noteworthy that when the ESS participates in the power regulation of the power grid, the disadvantages brought by the anti-peak-shaving of wind power can be well compensated. As shown in Figure 8, load 3 represents the total net load considering energy storage and multi-energy regulation, and load 2 represents the total net load without considering multi-energy regulation. ESS charges at night when the power load is low and the wind power output is large and discharges during the day when the load is relatively high and the wind power output is small. As a result, peak cutting and valley filling of the power load are achieved.

Natural gas is simultaneously converted into heat and electricity through CHP units and only into heat through the gas furnace. Figure 9 shows the total purchase of natural gas at each time period, from which we can observe that the amount of gas purchase of gas is larger at the time period when the power load is high than that at the time period when the power load is low. As shown in Figure 8, the difference between curve load 2 and load 3 is mainly due to the multi-energy complement. In addition, the EV load is relatively high at time periods from 7:00 to 17:00. To summarize, increasing the power supply of CHP units can effectively reduce the peak load, but this impact on valley filling is not obvious.

As discussed above, ESS can play an important role in power systems. Reasonable ESS scheduling decisions can alleviate the impact of wind power on the distribution network and promote peak shaving and valley filling. Figure 10 and Figure 11 show the charge–discharge curve of ESS and the SOC state of charge of ESS, respectively. We can note that for time periods from 0:00 to 7:00, the net load is relatively with low power consumption and high wind power production, and ESS charges to satisfy the load electricity demand. The SOC state of charge rises in preparation for the next stage of the discharge process. For time periods from 8:00 to 11:00, the electricity price is high and the net load is relatively high. ESS reduces the load peak by discharging. From 16:00 to 19:00, the load reaches the peak again, and ESS continues to discharge, which reduces the power purchase cost from the main grid. At the same time, in order to ensure the normal operation of the next charge and discharge cycle, the ESS charges and returns to the original state of charge at time periods from 20:00 to 24:00.

We introduce an index δP _a to measure the average peak-valley difference: δ P a = ∑ t p e a k P i , t p e a k / N p e a k − ∑ t v a l l e y P i , t v a l l e y / N v a l l e y . (45)

Table 1 lists mean peak-valley differences and the total operation costs under different ESS capacities. As shown in this table, the increase in ESS capacity positively reduces the average peak-valley difference of power demands and reduces the total operating costs. This is because a higher energy storage capacity can provide more sufficient charging power in load valley and discharging power in load peak, which also reduces the purchase of power at the peak time, thus reducing the operating cost.

To control the load peak-valley difference to an ideal level, we add an additional linear constraint (46) into the original model. Meanwhile, we set the energy storage capacity to be fixed at the original value δ P a ≤ λ . (46) λ denotes the load peak-valley difference which we want to get. By modifying the value of λ, we can obtain the corresponding operating costs under different requirements of average peak-valley difference, and the results are summarized in Table 2.

Generally, a more ideal peak-valley difference increases the total operating cost, which indicates that the operating cost and the requirement of peak shaving and valley filling are contradictory. The reason is that if we set the value of λ to be small, we need to increase the net load during the valley time. Therefore, part of the wind power will be curtailed during the load valley, and additional electricity needs to be purchased from the main grid, which reduces the peak-to-valley difference and inevitably increase the total operating cost.