Chance-Constrained Optimization for Integrated Local Energy Systems Operation considering Correlated Wind Generation

Energy hubs, which integrate multiple energy vectors through converters, can enhance the value of Integrated Local Energy Systems (ILES) via increased flexibility and reduced costs. However, uncertain renewable energy and the non-convex, non-linear properties of energy flows complicate the modelling and operation of energy hub systems. This paper develops chance-constrained optimization methods for planning and operation of energy hub systems under uncertainty. The non-linear formulations of power and gas flows are relaxed by convexification methods, leading to a formulation of Second Order Cone Problem (SOCP), which can be efficiently solved to global optimality. The correlation between geographically close wind generators connected to the hub systems is modelled by establishing their relation using Gaussian copula. The proposed chance-constrained optimization is demonstrated on a six-hub system within a multi-vector energy distribution network with 7 electrical buses and 7 gas nodes. The value of different levels of system integration through the installation of energy hubs is investigated. The results show that by combining system integration via energy hubs with chance constrained operation, the proposed method can reduce operating costs and increase renewable energy yields, thereby benefitting hub system operators and customers with reduced energy infrastructure investment and energy costs. Keywords— Chance-constrained programming, copula, correlation, distribution network, energy hub, Integrated Local Energy Systems, uncertainty.


Introduction
Integrated Local Energy Systems couple the electricity and gas distribution networks to meet electrical and heat demand, thereby unlocking new business models and enhancing the value of existing infrastructure. However, their operation faces challenges due to the complex interaction between the two networks, high penetrations of distributed energy resources, and uncertainty arising from both supply and demand [1]. The concept of energy hubs -nodes at which energy can be converted from one form to another and passed between energy vectors -can help address challenges caused by uncertain renewable generation and significant changes in customer demand profiles. Energy hubs can coordinate various energy vectors to meet demand for gas, heat, and electricity, thereby significantly improving the efficiency of, and customer interaction with, energy systems [2,3]. ILES operators, customers, and generator owners can benefit from the flexibility of energy hubs, but the added complexity requires optimally scheduling the energy flows across all vectors.
Optimization of the energy hub system involves taking optimal decisions for energy flows and operation of energy hub components, including converters which enable multi-vector operation, which have been investigated in the existing literature: [4] presents optimization for renewable-based integrated electricity-gas-thermal networks incorporating demand response. [5] proposes a bi-level model to optimal coordinate the operations between the energy hubs and distribution network. [6] integrates the concept of smart energy hubs into a distribution network using integrated information gap decision theory. [7] proposes capacity planning of energy hubs with a robust, data-driven approach. [8] develops a generalized heuristic approach to solve the optimal energy flow problem in the energy hub system. [9] and [10] propose a decomposed approach to optimize the behavior of, and flows between, interconnected energy hub systems using Particle Swarm Optimization. The above literature all investigated the optimization of multiple energy hubs considering the energy flows with different objectives. However, most of the solutions rely on linearization of either the non-convex AC power flow or gas flow, which reduces the accuracy of the results; others simulate the full power flow and gas flow equations, yielding a non-convex problem formulation which fails to ensure a global optimality and requires large computational efforts for solution using a heuristic solver. A convex representation of a multi-carrier network and the energy hub components would therefore be of significant value to ensure fast and accurate solution of the optimization problem.
Energy hub operations and energy flows between hubs could be significantly affected by uncertainty in production and consumption of energy, most notably arising from renewable electricity generation.
Several options exist to solve the energy hub optimization problem with uncertainty: robust optimization has been applied in [11] and [12] for planning and dispatching demand response; chance constrained programming and distributionally robust optimization are respectively applied in [13] and [14] to optimize the energy hub system; Monte Carlo simulation is applied in [15,16] to model the uncertain inputs for energy management; and scenario-generation methods are applied in [17,18] to simulate the uncertainty in load, price, and renewable generation; information gap decision theory is applied in [19] to optimize the system under uncertainty from electric vehicles of energy hub system. However, the uncertain outputs from renewable generators at a given hub may be correlated with those at other hubs which are geographically nearby [20]; these correlations are generally ignored in the above literature, and could lead to unplanned curtailment of renewable energy or violation of network limits.
This feature challenges optimal system operation, particularly when renewable power generation is high. Therefore, the modelling of correlated uncertainties may be of significant importance in ensuring optimal energy hub operation. The dependent uncertain variables can be formulated by joint distributions, but they lead to intractable optimization problems. Some expansion methods, introduced in [21], can estimate the probability density functions of the joint distributions with joint moments and cumulants. However, the calculation of joint cumulants becomes very complex as the number of uncertain variables increases [21]. Hence, expansion methods are inappropriate to solve the optimization problem with multiple correlated uncertain variables. [22] employs the Cholesky decomposition to capture the correlations, but the method can only be applied to reflect the linear relationships.
The stochastic electricity price under electricity market could also affect the operations of energy hubs, and hence the energy flows between them. Reference [23] proposes a bidding strategy for wind and thermal generators participating in the day-ahead market by stochastic programming, considering the uncertainty in energy prices and wind energy; reference [24] employs a cooperative game approach to present a coordinated bidding strategy for Power-to-Gas (P2G) facilities and wind generators by taking into account the uncertainty of forecast error of wind generation; reference [25] proposes an optimal day-ahead scheduling for P2G storage and gas load management in electricity and gas markets. This paper focusses on the effect of uncertainty in production and consumption to the energy hub system and its relevant distribution systems; it is reasonable to assume that these systems are sufficiently small to be price-takers, and therefore the influence of energy markets is not examined. Future research, considering systems in which energy hubs make up a substantial proportion of the supply and demand base, could address this.
Based on the existing literature, there is a need to develop an optimization model for energy hubs within the ILES which can: 1. Pursue the global minimum efficiently and with high accuracy; 2. Address the uncertainty within the system; and 3. Accounting for correlation structures which exist between sources of uncertainty, such as renewable generator outputs. This paper achieves the first of these by applying the convexification method proposed in [26] to relax the non-convex power flow in a radial network; relaxing the non-convex gas flow (for a radial network) [27] and thereby formulating the ILES optimization problem as a convex Second Order Cone Programming (SOCP) problem. The second is addressed through Chance-Constrained Programming (CCP), which is applied to optimize the energy hub system under uncertainty. Finally, the correlation of wind speed forecast errors in a multi-energy hub system is modelled by Gaussian copula, which can fit the joint distributions with high dimensionality and incomplete parameter estimation of uncertain variables [28].
Some distribution system constraints can safely be violated so long as this is infrequent and for a short duration; most notably, thermal constraints are set probabilistically, and in overhead line networks they are likely to have a higher thermal capacity during periods of high wind generation due to wind cooling effects [29]. Energy hubs can help to manage these constraints by using converters to mitigate the influence of uncertain renewable generation, while allowing an acceptably small probability of overload to significantly reduce operating costs. This paper proposes the use of chance constraints to restrict the power flow limits and applies CCP to minimize system cost while addressing the wind generation uncertainty with an acceptable probability level.
The main contributions of the paper are therefore:

i)
A novel method is presented to enable optimal planning of integrated local energy systems under uncertainty; ii) The gas and electricity networks are transformed to a convex, second order conic model, enabling the globally optimal operating regime to be quickly and accurately obtained; iii) Operation under uncertainty is enabled by applying chance constraints to the power flow limits in the electrical network; iv) Systems integration is enabled through energy hubs; the impact of the number of hubs installed is quantified; v) Correlations between sources of uncertainty are modelled by Gaussian copula, and the implication of neglecting these correlations is demonstrated.
The rest of the paper is organised as follows: the modelling method for energy hub systems with convex-relaxed electricity and gas networks are described in Section 2. The correlated outputs from renewable generators are presented in Section 3. In Section 4, the optimization problem is formulated and the methodology is introduced. Case studies are discussed in Section 5 and conclusions are presented in Section 6.

Energy Hub Modelling
A general energy hub model, which includes all the converters applied in this paper, is described in in this section and illustrated in Fig 1. The energy hub system utilizes electricity, gas, and wind energy to meet electrical and thermal demand. Each hub contains a Combined Heat and Power (CHP) system which converts gas to heat and electricity and a Ground Source Heat Pump (GSHP) which converts electricity to heat. Wind Turbines (WT), which convert kinetic wind energy into electric power, are also included in some energy hubs. The application of GSHP benefits from high and stable ground temperature, it will also be essential to decarbonize the energy system [30]. The outputs of GSHP (PHP_out) is quantified in (1) where, COP denotes the Coefficient of Performance of the GSHP and PHP is the electrical power input to GSHP and PHP_out is the heat output of the GSHP.
The wind generation is expressed in terms of the wind speed 3 by piecewise function in (3) [31], where Prated is the wind turbine rated power. vci, vrd, vct respectively represent the cut-in, rated, and cutout wind speeds.
The energy transformation through the energy hub system in Fig. 1 is expressed in (4), where the output matrix (L) is equal to the converter coupling matrix (C) multiplied by the input matrix (P).
Equation (4) shows that the electrical demand (Lele) and thermal demand (Lth) at time step t are satisfied by utilising electricity (Pele), natural gas (Pgas), and wind energy (Pwi) with the converters. νe is the dispatch factor, which represents the fraction of electricity input into the GSHP over the total electricity consumption.

Convex Relaxation of Electricity Network
Each energy hub is connected to a distinct bus within the electricity distribution network. The complex power equilibrium at bus i is expressed by the power balance equations in (5) and (6). in the admittance matrix, Ne is the total number of the buses that connect with bus i. These equations are non-linear and non-convex. However, according to [26], (5) and (6) can be transferred to the linear formulations in (10) and (11) by introducing ancillary variables ui and vi.

Convex Relaxation of Gas Network
Each energy hub is also connected to a specific node in the gas distribution network. The nodal gas flow balance for the gas node i is expressed in (14) [27].
Where qgw,i is the gas injection at node i, qd,i denotes the gas load demand, :: Z and :% &' represent the gas inflow and outflow between node i and k, Ng is number of nodes connected to node i, and :` is a coefficient which represents the compressor's consumption ahead of the gas inflows.
The average gas flow qik between node i and k can be formulated as in (15)- (16).
However, in this paper optimal scheduling the of energy hub operation is computed for one hour, and it is reasonable to assume that the gas flow directions do not change within a day [27]. Hence (16) can be re-written as (17).
Similar to the relaxation of power flow, (17) can be relaxed to (18), which is a convex SOCP formulation [32].

Modelling of Correlated Uncertain Variables
Uncertain variables are often not independent of one another, and failure to account for this correlation can result in underestimation of the impact of the uncertainty on the system. In this paper, the uncertainty studied is wind generator output, but the method described would be valid for any soucre of uncertainty. Wind generation is the key technology to decaronize energy system, and facilitate microgrid to meet consumers demand [33,34]. Wind generators in energy hub systems can be geographically close to the one another, which yields strong correlations between wind speed at each site considering the forecast error -this correlation has a significant impact on the operation of the system and is therefore explicitly modelled in this paper. By adding up the forecasted wind speed and the data of forecast error, wind speed data from two geographically close sites are plotted in Fig. 2, which shows that the wind speeds at the two sites are highly correlated. Additionally, their probability distributions could be different. As suggested in previous literature [20], wind speeds can follow Gaussian, Weibull, Beta, or lognormal distributions, which complicates the formulation of joint distributions, particularly if the wind speeds at each site follows a different distribution. Nevertheless, copula are capable of tackling this problem to fit joint distributions with high dimensions and non-linear correlations among uncertain variables [28]. Therefore, Gaussian copula are applied in this paper to simulate the correlation of wind speeds at different wind generator locations and generate correlated copula samples for solving the CCP. Sobol sequence based Quasi-Monte Carlo Simulation (QMCS) is applied to improve the computational efficiency in generating correlated copula samples in this paper [28].

Definition and Application of Copula
By applying the probability integral transform to multiple random variables( b , p , … , Z ), the uniform marginal distributions can be derived as ( b , p , … , Z ) = ( b ( b ), p ( p ), … , Z ( Z )) .
Copula C is defined as the joint cumulative distribution function of ( b , p , … , Z ) in (19).
According to Sklar's theorem, the cumulative distribution function of random variables  (20).
Gaussian copula is applied in this paper to model the correlations between uncertain wind speeds with Sklar's theorem as the basis. The procedure of generating copula samples for uncertain variables with arbitrary distributions is introduced as follows: Step 1) Determine the marginal distributions of the uncertain variables and estimate the distribution parameters based on historical data.
Step 2) Employ kernel density estimation method to estimate the cumulative distribution function (CDF) of the uncertain variables, and the correlation coefficient between them. Then compute the copula structure to capture the joint distribution of the uncertain variables, which is achieved by Sobol sequence-based QMCS to generate correlated samples on the unit scale.
Step 3) Transfer the correlated samples from unit scale to the original scale of the data, which is achieved by applying the corresponding inverse cumulative distribution function on the unit-scale correlated samples based on the distribution parameters estimated in step 1.
In terms of generating copula samples for wind speeds in this paper, 100 correlated samples of wind speed forecast errors in unit scale are generated as examples shown in Fig. 3 (a), which are derived by implementing steps 1) and 2) based on historical data of wind speeds. Fig. 3 (b) can then be accordingly obtained by using the inverse transform method on the samples in Fig. 3

Chance Constraints on Power Flow Restrictions
In this paper, chance constraints are applied to the power flow limits to allow a small but acceptable probability of the current exceeding the conductor rating. This reduces the need to schedule the operation of the system to cater for the most extreme uncertainty values, delivering lower operating costs with a predetermined and acceptable increase in risk. Resolving the resulting CCP problem yields a minimum system cost while addressing the wind generation uncertainty with an acceptable probability.
The power flow between bus i and k can be formulated as in (21) Chance constraint is applied to restrict the power flow between bus i and k as: Where Pr is the minimum acceptable probability the chance constraints will not be violated; α is the probability the power flow is less than the maximum power flow :`, † ‡ˆ, which indicates that the power flow is allowed to exceed its line capacity with known probability α. Nevertheless, allowing a higher

CCP Optimization Formulation
This section presents a generic CCP problem formulation for multiple energy hubs connected by The decision variables of the CCP problem include: average gas flow, gas inflow and outflow of each branch, gas pressure at each node, ancillary variables ui, Iik, Rik, power and gas injected to each energy hub, gas input to each CHP, power input to each GSHP, electricity injection into the distribution network, gas output from each gas source, wind generation curtailment Pwc. The above decision variables for each scenario s are included in the vector xs as in (25).
In (25) In summary, the energy hub CCP problem is formulated as follows: Objective: Resolving the CCP problem means to the find the optimal decision variables (25) to minimize (23) for all scenarios, while the uncertain wind generator outputs are randomly distributed according to their joint distribution.

Transforming Chance Constraints into Deterministic Constraints
The optimization problem is a stochastic problem with chance constraints (22). In this paper, the chance constraints are transformed into deterministic constraints based on the generated copula samples. The transformation is carried out using the following procedure: First, the power flow :` between bus i and k can be re-arranged and written in terms of the polynomial function of decision variables D(x) and uncertain wind generation 3 as shown in (27), where b is the coefficient related to the uncertain variable 3 in expressing the power flow :`. For example, within this context, the power flow between bus 5 and 4 in Fig. 6 can be expressed as in (28), in this example b is derived as 1.
The chance constraints (22) can be re-arranged as in (29)

Overall Flowchart of the Methodology
The complete optimization procedure is illustrated in Fig. 5. To formulate the optimization problem, the CDF of the wind generators' forecasting errors must first be obtained using Gaussian copula to calculate the quantile for the chance constraints. Additionally, both the power and gas networks are convex relaxed, load and energy price data are input to establish the optimization problem. The original non-convex stochastic problem is transformed to a convex deterministic problem according to Fig. 5; commercially available solvers can be applied to solve the problem with global optimality.

Overview of the ILES
The proposed optimization scheme has been demonstrated using an ILES of a multi-vector distribution network with 7 electrical buses and 7 gas nodes, based on the model in [27]. Schematic diagrams of the electricity and gas distribution networks are depicted in Fig. 6. in blue and green, respectively. The electrical network comprises two radial feeders, with three geographically close wind generation sites at buses 4 and 5. The wind speed forecast error of the wind generations sites at buses 4 and 5 are assumed to respectively follow Normal and Beta distributions. Buses 2-7 in the electrical network correspond to nodes 2-7 in the gas network. Each coupling point between the electrical and gas systems is considered as an energy hub. Each energy hub is equipped with a CHP and GSHP. Two gas sources -which could represent conventional gas wells or green gas sources such as anaerobic digestors -support the gas network at nodes N5 and N7. The pipelines linking N1 to N6 and N1 to N3 are active The optimization was implemented on the ILES using a time step of one hour. The transformation between the gas load (m 3 ) and energy (MWh) was achieved by using calorific value, which is taken as 40 MJ/m 3 in this paper. Additional input parameters, taken from [13,27,28], are shown in Table 1. The scenarios representing the uncertain wind generator outputs were generated by Sobol sequence based QMCS which was been applied to generate the copula samples in section 3.1. 100 scenarios were generated, which are able to represent the uncertainty with high accuracy [28]; each scenario has a probability of occurrence of 1%. The CCP problem was solved using Gurobi [36].  The optimization was carried out with the chance constraints restricted to different probability levels to investigate how the CCP affects the optimization and energy hub system operations, and examine the sensitivity of the optimization to the chance constraints' probability. 11 probability levels between 0.1% and 10% were adopted to restrict the probabilities of constraint violation (PoCV) (α in (22)),

Impact of Chance Constraint Probability Levels
indicating that the probability of the power flow remaining below the line capacity varies from 99.9% to 90%. When the PoCV is 10%, the computation time of implementing the CCP with 100 scenarios is 125 s on average, which is derived based on a MacBook with i5 processer and 8 GB memory.
The change in ILES operations against different POCV are shown in Fig. 7, Fig. 7(a) shows the expected total network cost, Fig. 7(b) depicts the expected electricity and gas inputs to the network, and the total wind generation curtailment, respectively in blue, green, and red, Fig. 7 ( X ) = ∑ Pr X • XX ab (31) Fig. 7 indicates that the operations of the multi-carrier network are slightly affected when the PoCV is higher than 5%. When the PoCV reduces from 5% to 0.1%, reducing the acceptability of exceeding the power flow constraints, significant changes in system operations can be observed from Fig. 7. When the PoCV reaches 0.1%, the capacity of transmitting the available wind generation at bus 4 and 5 to other load buses is reduced, as a result, the expected total curtailment of the wind generators at bus 5 and 4 increase from 0.61 MWh to 2.17 MWh as seen in Fig. 7(b); the expected gas injection to the network increases from 0.39 MWh to 1.42 MWh to fill up the demand by CHP since the outputs from the wind generators are reduced; similarly, the expected power injection to the network increases from 0.05 MWh to 0.38 MWh. Conclusively, the reduction of PoCV leads to the increase of wind generation curtailment, gas and electricity inputs, hence significantly increases the system cost from £9.73 to £42.82 as seen in Fig. 7(a) when the PoCV decrease from 5% to 0.1%. Additionally, due to the high efficiency of the CHP, and lower gas price compared with the electricity price, the ILES inputs more gas than electricity to meet the demand when the PoCV is reduced, to minimize the system cost. Although the cost and energy inputs increase, the system power loss decreases from 1.68 MW to 0.99 MW with the reduction of PoCV as seen in Fig. 7(c), because a high probability of power flow exceeding the line capacity leads to higher utilisation of assets and therefore higher losses. It is notable that the cost is lower when the PoCV is higher, in spite, rather than because, of the higher losses.
The variation in expected power flows between buses, gas input to each CHP, and power input to each GSHP against the changing of PoCV are depicted in Fig. 8. All expected power flows decline with the decrease of PoCV, which is in accordance with the definition of using chance constraints to restrict power flows. The larger variations in power flow indicate the larger influence of the chance constraints, as seen in Fig. 8 (a). also decrease when the PoCV decreases. P16, and P67 are less affected because the majority of uncertain wind generation is consumed by the energy hubs at bus/node 2-5, the remaining uncertainty is less likely to cause these power flows to reach the line limit. Overall, when the PoCV reaches 0.1%, less uncertain wind generation can be transmitted from buses 5 and 4 to other buses, the power inputs to the GSHP at the other hubs decrease due to the reduced availability of wind energy to support them. Correspondly, in order to collaboratively meet the demand with the GSHP, the gas inputs to the CHP increase as indicated in Fig. 8(b). The expected gas input to hub 7 is most affected, increasing from 0.34 MWh to 1.29 MWh when the PoCV is reduced to 0.1%. This is because more generated wind power is curtailed with the decrease of PoCV, the farthest hub is most affected since less wind power could be transmitted there to support the load, and consequently The above analysis demonstrates the high dependency between the electricity and gas networks: thanks to the energy hubs, more gas can be input to the CHP to support the electricity and thermal demands when the power flow restrictions are stricter.

Effect of Energy hubs
Energy hubs increases the system flexibility because the energy demands no longer rely on a single energy vector. When the electrical power flow is restricted or electricity prices are high, the electrical loads can be met through the use of CHPs and increased flows in the gas network. However, the capital cost of the energy hubs including the CHP and GSHP is generally expensive, and the operation of some energy hubs provide less value to the system. This section investigates how the existence of energy hubs affects the operational costs and carbon emissions of the ILES. 7 comparative cases are considered by employing different number of energy hubs in the ILES in Fig. 6. Where there is no coupling between the and gas and electricity networks, the electrical load is met through electricity import and/or wind generation and the thermal load is met by gas furnaces with an efficiency of 95%. The carbon emissions of using per unit grid electricity and gas are taken as 0.10 kg/MWh and 0.18 kg/MWh based on the average emission data in Scotland [37], the utilization of wind generation present zero carbon emission.
The 7 cases consider the number of energy hubs from 0 to 6, the sequence of deploying a new energy hub in each case is from the coupling point that is the nearest to the wind generation sites, namely hub 4, 5, 3, 2, 6, and 7. The optimal expected system total costs and expected carbon emission are shown in Fig. 9 (a) and (b) respectively. Fig. 9(a) suggests that the system cost decreases when the number of energy hubs rises in the multicarrier network. A relatively large cost decrease can be observed when the number of energy hubs increases from 3 to 4, corresponding to the installation of the energy hub at bus 2. The electricity and thermal demands at bus 2 are higher than the other buses, so adding the energy hub at this bus leads to a higher cost reduction. Similarly, Fig. 9(b) indicates that the carbon emission of the multi-carrier network decreases with the installation of energy hubs. This is because a higher proportion of thermal demand can be met by using the carbon-free wind generation through GSHP when there are more energy hubs deployed in the network. Both the wind generation curtailment and the import of electricity Fig. 9. Expected optimal cost with different number of energy hubs and gas are reduced, leading to the reduction of carbon emission. The greatest reduction in carbon emissions is observed with the installation of the energy hub at bus 2 due to its high thermal energy demand.

Effect of Correlation
The value of using copula to model the correlations between wind speeds is analysed in this section by considering a comparative case is in which the wind speeds are assumed to be independent. Similar to Fig. 3 (b), the distribution of independent wind speeds can be derived as shown in Fig. 10 by fitting the same data. The wind speeds observed at the two sites are independent, and randomly distributed according to their own distributions.
The impact of ignoring the correlation was investigated using the following process: 1. Carry out optimal scheduling for cases with independent wind speeds and PoCV varying from 0.1% to 10% 2. Carry out a Monte Carlo simulation to evaluate the performance of the optimal schedule which neglected the correlation but with the correlation present (i.e. the correlation exists but the scheduler does not consider it) 3. Carry out the same Monte Carlo simulation to evaluate the performance of the optimal scheduler when the correlation is included 4. Compare the specified PoCV to the frequency with which the constraints were violated during the Monte Carlo simulation; the results of this are shown in Fig. 11. As shown in Fig. 11, by modelling the wind speeds correlations with Gaussian copula, the optimal operations satisfy the chance constraints on power flow limitation for majority of cases. However, if the wind speeds are considered to be independent, the optimal operations could further violate the predefined PoCV of the optimization for 1.5% to 8%, posing a risk for the system.

Exactness of Flow Relaxation
The exactness of power and gas flow relaxations is examined in this section. According to [38], the Second-Order Conic Relaxation (SOCR) is exact if every solution of SOCR also solves the optimal power flow problem. Consequently, the results derived from SOCR should meet the original unrelaxed constraints to ensure its exactness. In this paper, the optimal solution from the relaxed forms of 2 :`≥ the accuracy of computing gas inflow, outflow, gas sources outputs, and inputs to CHPs are unaffected, because they follow linear relationships, which are met by the optimization of SOCR; only the accuracy of calculating gas node pressures, namely pi and pk, is affected, because they follow non-linear relation with gas inflow and outflow, and the relaxation is not exact.
To ensure the feasibility of the optimal solution, we propose to take the gas inflows and outflows Fig. 11. Comparisons between optimal results derived from SOCR as inputs, the gas node pressures from SOCR as the initial searching point, the same pressure limitation as constraints, and employ non-linear solver to optimize gas node pressures in order to satisfy (17). This ensures the accuracy of solution and adds only 0.32s of extra computation time for optimizing a scenario. This procedure does not change the overall optimal solution but updates the nodal pressures to follow (17). An example result is shown in Fig. 12. The SOCR results are all at or close to the lower bound, due to the pressure minimization term included in the objective function; the true, non-linear gas flows results indicate higher pressures at some nodes, which result in all solutions being within the acceptable pressure range. Feasible solutions for the gas flow with gas pressures that satisfy (17) could always be reached by employing the non-linear solver to all 1100 optimal results from SOCR.

Conclusion
This paper proposes a chance-constrained optimization scheme for Integrated Local Energy Systems which: 1) convex-relaxes both the electricity and gas network models, enabling efficient optimization which reaches global optimality with commercially available solvers; 2) Accounts for uncertainty and enables a compromise between risk and cost by allowing line ratings to be exceeded with a specified probability, and 3) captures the correlations between geographically close wind generation sites via copula. The main findings are as follows: • The operating cost of the ILES is significantly reduced by allowing even a low probability of chance constraint violation. However, this delivers diminishing returns as the probability increases, with little improvement between a 5% and 10% PoCV.
• The number and location of energy hubs in the system has a significant impact on the system operating cost and carbon emission. Specific energy hubs -particularly those connected to larger demands -led to a much greater reduction in operating cost and carbon emission than other hubs, some of which had almost no effect. • Copula effectively and efficiently represents the correlations between wind speeds. Ignoring the correlations yields an underestimation of uncertainty, which may put the safe operation of the system at risk.
By explicitly modelling the uncertainty from renewable generation, the proposed optimization benefits the system operator with optimal determination on system operations to reach the minimal cost depending on the acceptable tolerance for chance constraints. It is particularly beneficial for the optimal operation of large-scale, multi-carrier systems due to the convex SOCP formulations. Furthermore, the proposed optimization could be extended to solve optimal planning problems from an investor's perspective. The rated power of the wind generations can be re-evaluated to maximize the profits while avoiding curtailment considering wind speed uncertainty, position and capacity of energy hubs. To further improve the system flexibility in coping with uncertain renewable generation, future work could model power-to-gas technology to establish the conversion between abundant renewable generation and gas network, thereby avoiding the curtailment of renewable generation.