A ship battery exchange model is developed to account for the spatiotemporal characteristics of port vessels. As operational vessels move back and forth between the sea and the port, their routes are variable. Consequently, the distance of each journey segment, denoted as d, is modeled using a log-normal distribution, with the probability density function described as follows: f D d = 1 d σ D 2 π exp − ln ⁡ d − μ D 2 2 σ D 2 (1) where σ _D represents the standard deviation, μ _D represents the expected value of the travel distance.

Since the Weibull distribution is highly adaptable and flexible for modeling random data (Abubakar et al., 2024), it is employed to fit the starting time t _a for each journey segment. The probability density function of the Weibull distribution for the starting time t _a is expressed as follows: f t a , k , c , γ = k c t a − γ c k − 1 e − t a − γ c k , t a − γ > 0 (2) where t _a is the random variable, k is the shape parameter, c is the scale parameter, and γ is the location parameter.

A Poisson distribution is used to model the arrival behavior of EB at the port. The specific distribution is represented as follows: P X = x t b = λ t x t b e − λ t x t b ! x t = 0 , 1 , 2 , . . . , n (3) where P (X = x _tb ) refers to the probability of x _tb boats arriving at time t _b , λ is the average frequency of vessel arrivals within a unit time interval.

The EB sailing time t _f and docking time t _p are defined as follows: t f = t b − t a t p = 24 − t f (4)

The simulation and prediction of ET load need to account for factors such as the truck’s weight, driving speed, and the slope of the driving area. The electrical power is derived from the truck’s mechanical power. The prediction formula for the energy consumption of the ET is as follows: H v = d m g K r ⁡ cos θ + 1 2 K a ρ K d v 2 + m g ⁡ sin θ 1000 K p K h (5) where d is the driving distance of ET, m is the total weight of ET, g is the gravitational constant, K _r is the rolling resistance coefficient, K _d is the aerodynamic resistance coefficient, K _a is the external surface area of ET, ρ is the air density, θ is the average road gradient, v is the instantaneous speed of ET, K _p is the battery efficiency of ET, K _h is the regression coefficient.

The speed-flow relationship of ET is defined as follows: V t = V 0 1 + q t Q ε ε = ε 1 + ε 2 q t Q 3 (6) where V ₀ is the zero-flow speed of the road, q is the traffic flow at time t, Q is the road’s ET capacity per unit of time, ε ₁ and ε ₂ are the adaptive coefficients based on the road grade or quality.

The charging model for electric trucks is: S c t = S c t − 1 + λ c P c t C c (7) where S _c (t) is the battery charge percentage of ET at time t, λ _c is the charging efficiency, C _p is the battery capacity, P _c (t) is the charging power at time t.

The annual total radiation in the port’s maritime area is relatively high, with unobstructed sunlight. Additionally, Based on stochastic optimization theory, the uncertainty of offshore wind power output follows an empirical distribution extracted from historical prediction error data. By integrating stochastic optimization theory with robust optimization theory, a fuzzy set is constructed to capture the uncertainty set of the potential true distribution. This set encompasses the worst-case scenario of the true distribution. The fuzzy set is characterized using the Wasserstein distance, which quantifies the distance between the empirical distribution P ^ w and the true distribution P ~ w . Specifically, the first-order Wasserstein distance is defined in integral form as follows: W P ^ w , P ~ w = inf Q ∈ Q P ^ w , P ~ w E Q ρ e ^ , e ~ (8) where inf (⋅) represents the infimum function, Q P ^ w , P ~ w denotes the set of all joint distributions between P ^ w and P ~ w . E (⋅) is the expectation function, and ρ(⋅) refers to the base distance used for calculating the Wasserstein distance. e ^ = e ^ 1 , e ^ 2 , … , e ^ t , … , e ^ T T is a vector composed of random variables at time t derived from the empirical distribution of prediction errors, while T represents the dimensionality of the random variables. Similarly, e ~ = e ~ 1 , e ~ 2 , … , e ~ t , … , e ~ T T is a vector composed of random variables at time t derived from the true distribution of prediction errors.

The Wasserstein distance is incorporated into the construction of the fuzzy set, resulting in a fuzzy set that can be interpreted as a sphere with a radius of θ. This sphere encompasses all potential true distributions of wind power forecast errors. Therefore, the Wasserstein fuzzy set F can be defined as: F θ = P ~ ∈ P 0 R T + 1 × S e ~ , v ~ , s ~ ∼ P ~ E P ~ v ~ ∣ s ~ ∈ S ≤ θ P ~ e ~ , v ~ ∈ ε s ∣ s ~ = s = 1 P ~ s ~ = s = p s , s ∈ S (9)

Wind power forecast errors are estimated using non-parametric kernel density estimation. Historical wind power forecast error data are clustered using the K-means algorithm to generate scenarios and their corresponding probabilities. The historical data sample E _t is defined as: E t = P t w , r − P t w (10) where P t w represent the wind power generation at time t in the port, P t w , r are the actual grid-connected power of wind generation at time t.

The output of PV panels is influenced by factors such as solar irradiance and temperature. Solar irradiance follows a Beta distribution. The offshore PV model is described as follows: f p v P p v = Γ α + β Γ α Γ β × P p v / P max p v α − 1 1 − P p v / P max p v β − 1 (11) where Г(·) represents the gamma function, while α and β are the shape parameters of the PV system. P max p v and P p v denote the predicted photovoltaic generation and the actual photovoltaic power absorbed by the grid at time t, respectively.

The gas turbine, as a traditional energy device, generates electricity by consuming natural gas, thereby providing stable power to the port. Compared to coal or oil-based power generation, it produces lower carbon emissions and offers a safer production process. To further enhance energy conservation and emission reduction, gas turbine equipment will be integrated with carbon capture and storage (CCS) technology. Additionally, waste heat from gas turbine power generation will be recovered through a waste heat boiler for secondary power generation or directly used for port heating via a heat converter. The power model is as follows: G G T = β G T 1 P G T 2 + β G T 2 P G T + β G T 3 Q G T = β G T ′ P G T η G T (12) where G _GT is the natural gas consumption of the gas turbine during operation, P _GT and Q _GT are the electric power output of the gas turbine and the recovered waste heat from the unit, β _GT1 , β _GT2 and β _GT3 are the corresponding coefficients that represent the relationship between the gas turbine’s electric power output and natural gas consumption, β ^’ _GT is the waste heat recovery index, η _GT is the electricity-to-heat energy conversion ratio of the gas turbine.

The gas boiler satisfies the port’s continuous thermal load demand by consuming natural gas. The heating power model is as follows: Q G H = η G H G G H (13) where η _GH is the heating efficiency of the gas boiler, and G _GH is the natural gas consumption of the gas boiler.

The port’s temperature control system primarily comprises the temperature regulation system for the living area and the cold chain system for the cargo area.

The time-series model for temperature control in the living area (Lin and Yi, 2000) is presented as follows: T r , t = θ 1 T r , t − 1 + φ 1 T g , t − 1 + ω 1 T w , t − 1 (14) T r , t = ∑ j = 1 J α j T r , t − j + ∑ j = 0 J β j T g , t − j + ∑ j = 0 J γ j T w , t − j (15) where T _g,t is the supply water temperature of the heating network at time t, T _h,t is the return water temperature, T _r,t is the indoor temperature of the heating building, T _w,t is the outdoor temperature coefficient, α _j , β _j , γ _j , θ ₁ , φ ₁ and ω ₁ are the physical parameters related to the thermal inertia of the heating system, J representing the order of thermal inertia of the system.

For temperature control in the living area, the predicted mean vote (PMV) index is utilized to evaluate the relationship between environmental temperature and comfort (Li H. et al., 2025). The PMV formula used by VPP-Ports is as follows: λ P M V = 0.303 e − 0.036 M + 0.028 { M − W − 3.05 × 10 − 3 * 5733 − 6.99 M − W − P a − 0.42 M − W − 58.15 − 1.7 × 10 − 5 M · 5867 − P a − 1.4 × 10 − 3 M 34 − t a − 3.96 × 10 − 8 f c l t c l + 273 4 − t r + 273 4 − f c l h c t c l − t a (16) where M and W represent the human metabolic rate and mechanical work, f _ct is the ratio of covered to exposed body surface area, h _c is the surface heat transfer coefficient, P _a is the water vapor partial pressure around the body, t _a , t _r and t _cl are the air temperature, mean radiant temperature, and clothing surface temperature.

For the cold chain system in the cargo area, which has stricter upper and lower temperature limits, the equivalent thermal parameter (ETP) model (Zhang and Lu, 2013) is utilized for simulating temperature inertia. The specific model is as follows: Q t L = η L P t L (17) T t n = T t w − Q t L R − T t w − Q t L R − T t − 1 n e − Δ t / R C (18) where Q t L is the total cooling power of the refrigeration unit during time t, R and C are the equivalent thermal resistance and capacity of the cold storage, T t n and T t w are the temperatures inside and outside the cold storage during time t, Δt is the time interval.

Energy storage devices in ports primarily include energy storage batteries and thermal storage tanks.

The energy storage battery model is described as follows: E t b a t = E t − 1 b a t + β c b a t P t b a t , c Δ t − β d c b a t P t b a t , d c Δ t (19) where E t b a t and E t − 1 b a t represent the actual energy levels of the storage battery at time t and t-1, respectively; β c b a t and β d c b a t denote the charging and discharging efficiecies of the storage battery, respectively; P t b a t , c and P t b a t , d c are the charging and discharging power of the energy storage battery from the grid at time t, respectively.

The thermal storage tank model is described as follows: E t h = E t − 1 h + β c h H t t t , c Δ t − β d c h H t t t , d c Δ t (20) where E t h and E t − 1 h represent the thermal storage levels of the tank at time t and t-1, respectively; β c h and β d c h represent the thermal storage and discharge efficiencies of the thermal storage tank; H t t t , c and H t t t , d c are the actual thermal charging and discharging power of the storage tank at time t, respectively.

To support energy conservation and emission reduction efforts and to achieve a “zero-carbon port area”, it is essential to integrate the port energy system into the carbon trading mechanism. This integration will help improve carbon emission management and optimize the energy structure. Currently, carbon trading is primarily divided into traditional and tiered mechanisms. The tiered carbon trading mechanism, introduced to more accurately simulate port carbon emissions, assigns prices to each tier based on the carbon emission weights within that range. This approach aims to incentivize VPP-Ports to engage in strategic carbon emission planning. The formula for the tiered carbon trading cost F _c (Lv et al., 2023) is as follows: F c = − β e E , 0 ≤ E < D β e E − D , D ≤ E < D + d β e d + β e 1 + α e E − D − d , D + d ≤ E < D + 2 d ⋮ β e d n + α e n − 1 ! + β e 1 + n − 1 α e E − D − n − 1 d , D + n d ≤ E < D + n d (21) where F _c,n is the carbon trading cost for each stage, F _c is the total carbon trading cost, E is the total carbon emissions of the port, β _e is the standard carbon trading price, D is the total carbon quota, d is the carbon emissions per unit in each interval, α _e is the tiered carbon trading price growth coefficient, When the port’s carbon emissions E < D, it indicates a surplus in the carbon quota, allowing the excess quotas to be sold to reduce total costs while providing users with corresponding rebates to incentivize them, promoting energy conservation and emission reduction.

For the general load consumed by the port, carbon emissions primarily originate from gas turbines. The historical baseline can be established by adjusting historical data, allowing for the calculation of the corresponding carbon quota D _s for this segment. The specific formula is as follows: D s = π Q + β D P (22) where π is the standard carbon quota coefficient, Q, P are the thermal power and electrical power consumed by the port, β _D is the electricity-heat power conversion coefficient.

The carbon emissions E of the port can be calculated based on its electrical power P and thermal power Q, using the following formula: E = η 1 + η 2 P + Q + η 3 P + Q 2 (23) where η ₁, η ₂ and η ₃ represent the corresponding carbon emission coefficients for the port.

(1) Electrical network

(1) Electrical network

By developing an optimized dispatch strategy for the aggregation of port resources, prescheduling results for the load can be obtained.

As shown in Figure 6, to meet the requirements of the port’s population activity area, the indoor temperature is maintained between 23°C and 26°C, ensuring a suitable environment for work and daily life. Meanwhile, the cold storage temperature is kept around −15°C to satisfy refrigeration needs.

Based on Figure 7 it can be concluded that the gas turbine maintains a stable output throughout the day. The electric heavy-duty trucks primarily charge between 20:00 and 06:00 the following day, while the ship battery swapping load is mainly concentrated between 11:00-14:00 and 20:00-24:00, and the battery charge and discharge to balance the power load response. Due to the absence of photovoltaic power generation at night, additional electricity is purchased between 01:00-02:00 and 19:00-24:00.

As shown in Figure 8, due to the high and stable output of the gas turbine, the absorption chiller’s output proportion is relatively high. During the 03:00-10:00 period, the cooling load is lower while the heating load is higher, so the waste heat boiler’s thermal output primarily serves the basic heating demand. As a result, the cooling network is supplied by electric chillers. From 11:00 to 19:00, when the cooling load demand increases, both the electric chiller and the absorption chiller work together to meet the cooling load balance.

As shown in Figure 9, during most periods, the waste heat boiler from the absorption gas turbine’s waste heat provides the majority of the thermal energy, with the surplus heat energy supporting the stable output of the absorption chiller. However, during the 05:00-10:00 and 19:00-23:00 periods, the heat storage tank exhibits significant heat storage behavior, and the waste heat boiler alone cannot meet the higher thermal load demand. In these timeframes, the gas boiler consumes gas energy to provide the additional thermal energy required.

As shown in Figure 10, the state of the gas network is consistent with the balance between the cooling network and the heating network. To ensure the balance of demand while improving economic efficiency, the gas turbine’s gas consumption in the gas network remains relatively high. The gas boiler only consumes gas energy when required by the system to meet the demand.

Based on the aforementioned scenario, the impact of changes in the benchmark price β _e and price growth rate α _e within the tiered carbon trading mechanism on the low-carbon economy of ports is analyzed. The results are illustrated in Figures 11, 12.

As shown in Figure 11, with the increase in the benchmark price, the system’s total cost and carbon emission cost rise correspondingly. Under the premise of maintaining economic viability, the port demonstrates an increased willingness to reduce emissions, leading to a gradual decrease in carbon emissions. However, when the benchmark price of carbon trading exceeds 250 CNY/t, the emission reduction potential for individual systems becomes limited, and port carbon emissions stabilize. At this stage, further guiding system emission reductions through price adjustments becomes less effective.

As shown in Figure 12, when the price growth rate is in the range of [0, 0.6], the carbon emission cost increases correspondingly. Under the constraints of system operating costs, equipment output becomes restricted, leading to a reduction in carbon emissions. When the price growth rate exceeds 0.6, the outputs of various system components gradually reach equilibrium, and carbon emissions begin to stabilize. With further increases in the price growth rate, both carbon emission costs and the system’s total costs continue to rise.

It is evident that to effectively guide the clean operation of port energy networks through the carbon trading mechanism, both the benchmark price and the price growth rate must be adjusted simultaneously. Achieving a coordinated configuration between tiered carbon trading and port scheduling optimization is necessary.

Taking into account the impact of port-specific loads and the carbon trading mechanism on the port’s electricity dispatch, this simulation focuses on the load response within VPP-Ports. The optimization objective is to minimize the total operational cost of the VPP. Three scenarios are defined for comparison:

Scenario 1: Only the general port load is considered, excluding the involvement of the carbon trading mechanism in the operation of VPP-Ports.

Scenario 2: The loads from electric trucks and ship battery swapping are considered, but the carbon trading mechanism is not incorporated into the operation of VPP-Ports.

Scenario 3: The operation of VPP-Ports considers both the loads from electric heavy-duty trucks and ship battery swapping, as well as the involvement of the carbon trading mechanism.

To analyze the impact of the carbon trading mechanism on the economic and environmental benefits of the port energy system operation, the dispatch costs for each scenario are shown in Table 3.

From Table 3, it can be seen that Scenario 1, which does not consider the operation of electric heavy-duty truck loads and ship battery swapping, relies more on traditional energy sources to supply the port’s work and living needs. As a result, it has the highest overall operational cost, the largest carbon emissions, and relatively low renewable energy output. Scenario 2, which considers the electric heavy-duty truck and ship battery swapping loads based on the port’s resource characteristics, reduces the consumption of petroleum and other fuels by heavy-duty trucks and ships. The total cost is reduced by 11.11% compared to Scenario 1, and renewable energy output increases by 4.08% compared to Scenario 1. Scenario 3, which considers both the electric heavy-duty truck and ship battery swapping loads along with the carbon trading mechanism, still reduces the total cost by 11.97% compared to Scenario 2, despite the inclusion of carbon trading costs. Renewable energy consumption increases by 8.23% compared to Scenario 2. This scenario reduces carbon emissions while meeting the port’s work and living needs and effectively increases the utilization of renewable energy.

In conclusion, accounting for the loads from electric heavy-duty trucks and ship battery swapping enables better adaptation to the port’s resource characteristics and more effective aggregation of these resources. Incorporating the carbon trading mechanism significantly reduces carbon emissions from VPP-Ports operations and encourages the use of clean energy. Therefore, considering both the port’s resource characteristics and the carbon trading mechanism ensures the economic and environmental performance of the port energy network, optimizing the efficiency and sustainability of VPP-Ports.