The increasing deployment of direct current (DC) sources and loads (e.g., solar power generation and electric vehicles), has promoted the evolution of conventional alternating current (AC) microgrids into hybrid AC/DC microgrids (HMGs) [1]. This transition to a hybrid structure capitalizes on the strengths of both AC and DC microgrids. It enables AC and DC loads to be directly powered by their corresponding sources within segregated sub-microgrids, thereby reducing conversion losses that would otherwise be incurred. At the heart of this synergy are bi-directional converters (BdCs) [2], which serve as conduits for power transfer between AC and DC sub-microgrids [3]. These converters enhance the overall robustness and reliability of both subsystems [4], [5]. Consequently, strategic interconnections for HMG networks are required to achieve a balance between reliability and economic efficiency in their operation.

Among existing studies focused on interconnection planning for HMGs, a handful have applied oversimplified planning models that ignore BdC efficiency characteristics altogether [8], [9]. Other studies have applied static BdC efficiency models in the planning process [10], [11]. However, static models fail to capture the dynamic nature of BdC efficiency arising from their highly non-convex characteristics, resulting in compromised interconnection planning solutions. This has been addressed to some extent by applying complex dynamic BdC efficiency models. For example, Li et al. [12] dynamically modeled BdC efficiency as high-order polynomials operating as functions of the BdC utilization rate. Moreover, the results of real-world examples demonstrated that these models effectively captured the dynamic nature of BdC efficiency. However, highly non-convex optimization problems are observed when using either the original high-order polynomial model or a mixed integer linear model that includes binary variables to ensure that each BdC does not engage in rectification and inverter operations simultaneously [13], [14], [15].

Efforts to address the high degree of uncertainty associated with RES outputs have widely relied on stochastic programming and robust optimization [16]. Stochastic programming is designed to optimize decisions under uncertainty by explicitly incorporating probabilistic information about uncertain RES outputs [17]. It models uncertainty through scenarios [18], each representing a possible realization of the uncertain RES outputs, with associated probabilities. The objective is to minimize the expected cost across all scenarios [19]. While it provides solutions that are optimal on average, it requires accurate probability distributions of RES outputs that are often difficult to estimate, and can become computationally intensive with a large number of scenarios.

In contrast, robust optimization handles RES uncertainty by assuming that uncertain RES outputs lie within a predefined uncertainty set rather than relying on probabilities [20]. It seeks solutions that remain feasible across all possible realizations within this uncertainty set, focusing on obtaining optimized planning solutions under worst-case scenarios within the pre-determined uncertainty set [21]. The main types of linear uncertainty sets currently used include box-based uncertainty sets (BbUSs) [22], budget-constrained uncertainty sets (BcUSs) [23], and convex-hull-based uncertainty sets (CHUSs) [24]. BbUSs define uncertainty as fixed intervals for each RES unit, forming a hyper-rectangle in the RES output space. While simple and computationally efficient, they often assume independence between RES outputs, making them overly conservative. To improve realism, BcUSs introduce a budget constraint on the total deviation of RES outputs from their nominal values. Compared with BbUSs, BcUSs often result in less conservative solutions while maintaining computational efficiency [25]. However, determining the appropriate budget value can be challenging and is often subjective. The budget is typically chosen based on expert judgment and preference, but inaccuracies in its selection can result in either overly conservative or overly optimistic solutions. This lack of objectivity can limit the robustness and economic performance of the obtained results. Compared with BbUSs and BcUSs, CHUSs offer superior capabilities for representing spatial and temporal correlations in RES outputs. Moreover, the conservativeness of CHUS solutions is much less than the solutions of BbUSs and BcUSs because the convex hull is the smallest convex set that encompasses all historical datasets. However, CHUSs have limited applicability because the number of vertices in a CHUS increases rapidly as the number of uncertainty sources increases, which greatly increases the computational complexity of the optimization process.

Efforts to reduce the conservativeness of robust optimization solutions without introducing excessive complexity include the development of distributionally robust optimization (DRO) methods [26]. These methods obtain optimized planning solutions under the worst-case distributions within distribution ambiguity sets [27]. The distribution ambiguity sets employed include continuous distribution uncertainty sets based on the Wasserstein or Kullback–Leibler metric [28], and discrete distribution uncertainty sets [29]. However, continuous distribution uncertainty sets can lead to optimization problems that are computationally intensive, especially when the size of the problem increases. Moreover, while discrete methods generate a robust solution under the worst-case probability distribution of predicted RES outputs, they fail to yield a solution that is robust against the actual worst-case scenario of RES output.

Despite the progress in robust planning for HMGs, there are still three main research gaps that need to be addressed. ① Although non-convex high-order polynomial models have been proposed to represent the dynamic efficiency of BdCs, no study has yet introduced a convex model that accurately captures this efficiency. The absence of such a convex model can result in convergence challenges and significant computational demands during the planning model solution process. ② The current uncertainty sets, such as BbUSs, BcUSs, and CHUSs, either tend to be too cautious and subjective, resulting in unnecessarily conservative strategies, or tend to be computationally intensive. This results in a challenging balance between producing robust solutions and keeping the computational load manageable. ③ When planning interconnections for multiple HMGs, the scale of the optimization problem can grow considerably. This complexity necessitates tailored algorithms designed specifically to speed up the solution process. To encapsulate these points, a comparison is included in Table 1 [8], [9], [10], [11], [12], [13], [15], [20], [21], [22], [23], [24], [27], [28], [29], which contrasts the features of our proposed approach against those of existing methods.

To actively tackle the challenges previously discussed in HMG interconnection planning, this present work proposes a tri-level robust BdC planning model that accounts for dynamic BdC power conversion losses and RES uncertainty. This model yields BdC investment strategies that are less conservative than those based on extreme scenarios within the BbUS. The key contributions are as follows.

Convexification of BdC efficiency model. This work pioneers the use of a convex model to capture the dynamic efficiency of BdCs, by employing the least-squares approximation (LSA) method to approximate the complex efficiency expression of BdCs. Additionally, the study establishes an exact convex relaxation of the lower-level model, representing a significant advancement in the field. This innovation enables the first incorporation of a convex representation of dynamic BdC efficiency within a robust BdC-based interconnection planning model, specifically tailored for HMGs.

The rest of the paper proceeds as follows. Section 2 formulates the tri-level robust BdC based interconnection planning model for solar power penetrated HMGs with dynamic BdC efficiency, while Section 3 models the DcUS of RES outputs. In Section 4, the robust interconnection planning model is tackled using a series of advanced methods: the LSA method, exact convex relaxation techniques, and the FPC&CG algorithm. Numerical results, showcasing the effectiveness of these methods, are discussed in Section 5. Finally, Section 6 summarizes the key findings and conclusions of this research.

This section first introduces the overall structure of an HMG and discusses the potential assumptions underlying the model. Subsequently, a tri-level robust interconnection planning model is formulated to address the uncertainties associated with RES units.

Fig. 1 illustrates a typical HMG structure involving two AC sub-microgrids (MGs) and two DC MGs with corresponding AC and DC loads, and distributed solar power generator (SPG), diesel engine (DG), and fuel cell (FC) units. The system is designed to directly support both AC and DC loads, minimizing conversion losses. For instance, DC loads (e.g., data centers or electric vehicles) are powered directly from the DC MG, while AC loads (e.g., residential appliances) are supplied directly by the AC MG. BdCs enable power flow between the AC and DC MGs, facilitating flexibility in power sharing across subsystems [15]. This structure leverages the benefits of both AC and DC MGs, supporting diverse energy resources and load types.

It is important to note that BdCs and multi-terminal energy routers (MERs) are two key technologies for interconnecting different microgrids [30]. While both enable power flow and energy exchange, they differ in terms of complexity, scalability, and cost. MERs are a promising solution for highly scalable and advanced multi-microgrid systems, whereas BdCs offer a practical and cost-effective alternative for many applications. This study focuses on BdC-based interconnection planning due to its lower cost, ease of implementation, and technological maturity. Nonetheless, the proposed tri-level robust interconnection planning model is flexible enough to incorporate MERs without requiring structural modifications.

To develop the microgrid interconnection planning model, several assumptions are made to simplify the analysis and focus on the study’s core objectives. These assumptions are as follows:

(1) The microgrids are assumed to primarily operate in islanded mode, meaning they are disconnected from the main grid. However, the framework can be easily adapted to include grid-connected operation when necessary.

(2) A single-stage planning approach is adopted for simplicity, where all decisions regarding resource allocation and interconnection design are made at a single point in time. Nonetheless, the methodology can be extended to a multi-stage planning framework to account for phased development or dynamic adjustments over time.

(3) The model focuses exclusively on steady-state interconnection design and does not consider dynamic stability aspects, such as transient or small-signal stability. These aspects are addressed in subsequent stages of microgrid control system optimization.

It is noteworthy that, interconnection planning typically prioritizes static characteristics such as capacity constraints, generator outputs, and steady-state operational efficiency, as these factors serve as the foundation for system design and ensure reliable operation under normal conditions. This study focuses specifically on long-term interconnection planning, which spans months to years and aims to achieve an economically efficient balance between generation and load. Unlike real-time control, this process emphasizes minimizing investment and operational costs over the long term, rather than addressing transient responses or instantaneous electrical phenomena.

Dynamic stability, which involves sub-second responses to disturbances such as sudden load changes or renewable power fluctuations [31], is excluded from this study for several reasons. First, the time scale of stability dynamics (milliseconds to seconds) is far shorter than the long-term planning horizon (months to years), making their inclusion irrelevant to the overall objectives of this model. Second, incorporating stability dynamics would significantly increase the complexity and computational demands of the planning process, as they involve highly non-convex behaviors requiring detailed time-domain simulations. Third, long-term planning decisions, such as investment costs (IC), generation scheduling, and load management, are not directly influenced by transient behaviors, making stability dynamics unnecessary for economic planning. Finally, in practice, real-time stability control and long-term planning are decoupled, with the latter focused on investment and expansion rather than immediate operational safety or control. For these reasons, this study excludes stability dynamics while concentrating on the core objectives of long-term microgrid interconnection planning.

At the lower level, the model minimizes HMG operation costs under both forecasted and extreme scenarios while ensuring operational constraints, such as power balance and BdC limits, are met. The tri-level structure integrates long-term investment decisions with short-term operational planning, ensuring robustness to uncertainties and adaptability to varying conditions. While this framework is centered on BdCs, it is flexible enough to incorporate alternative technologies like MERs, as the model’s structure remains unchanged, and only the input parameters (e.g., cost, scalability) need adjustment. This ensures the framework provides a practical and scalable solution for interconnecting microgrids. The rest of this section will present in detail each of the three levels of the planning model based on the system illustrated in Fig. 2.

The upper-level model aims to minimize the TC of annual BdC investments and annual HMG operation under the worst-case realization of the SPG output distributions, subject to BdC investment constraints.

where $\mathit{\Omega }$ represents the set of corridors; ${\mathcal{H}}_{\mathit{ij}}$ and ${\mathcal{H}}_{\mathit{ij}}^{+}$ are sets of all BdC-based interconnection lines and candidate interconnection lines, respectively; ${\xi }_1\left(\mathcal{P},\mathcal{Y}\right)$ denotes the value of the objective function obtained after maximizing the middle-level model. The decision variables are defined as follows: $\mathcal{P}$ represents the decision variable sets in the middle-level model; $\mathcal{Y}$ corresponds to the decision variable sets in the lower-level model; $I_{\mathit{ij},h}$ is the decision variable in the upper-level model that indicates whether BdC-based interconnection line h is to be constructed through corridor ij between microgrid i and microgrid j. The parameters in the upper-level model are defined as follows: $c_{\mathit{ij},h}^{\text{BDC}}$ represents the annualized IC of line h in corridor ij; $n_{\mathit{ij}}^{\max }$ denotes the maximum number of lines in corridor ij; $\sigma$ is the weight parameter applied to annualize the HMG operation cost; $\tau$ indicates the investment budget allocated for the construction of BdC-based interconnection lines. Eq. (1b) ensures that the investment budget is not exceeded; Eq. (1c) defines binary variables for each candidate BdC-based interconnection line $h\in {\mathcal{H}}_{\mathit{ij}}$; and Eq. (1d) sets $I_{\mathit{ij},h}=1$ for each existing BdC-based line $h\in {\mathcal{H}}_{\mathit{ij}}\backslash {\mathcal{H}}_{\mathit{ij}}^{+}$. Finally, Eq. (1e) limits the maximum number of BdC-based lines that can be installed in each corridor.

The middle-level model seeks the worst-case realization within $\mathcal{P}$ as follows.

where ${\xi }_{2,s}$ and ${\xi }_{2,s^{\prime }}$ denote the values of the objective function obtained after minimizing the lower-level model for the forecasted scenario $s^{\prime }\in {\mathcal{S}}^{\text{F}}$ and the extreme scenario $s\in {\mathcal{S}}^{\text{Dc}}$ in the DcUS, respectively, and the actual probability $p_{s^{\prime }}$ of the forecasted scenario is the decision variable in the middle-level model. As will be discussed in detail later, each extreme scenario corresponds to a single vertex in the DcUS. The parameters in this model are defined as follows: $p_{s^{\prime }}^0$ denotes the predicted probability in forecasted scenario $s^{\prime }\in {\mathcal{S}}^{\text{F}}$; $\theta$ is the weight factor enabling decision-makers to flexibly adjust their preference between forecasted scenarios and extreme scenarios; $\delta$ represents the uncertainty budget; $Z$ denotes the number of data points for historical SPG outputs; $\alpha$ is the confidence level.

Here, the l₂-norm metric is employed in Eq. (2b) to capture the correlation of SPG output probability distributions more accurately in uncertainty sets than that obtained using the l₁-norm metric in previous studies [29]. This approach aims to minimize unnecessary conservativeness in interconnection planning solutions. Although solar power has been used as an example, the same method can be extended to formulate distribution uncertainty sets $p_{s^{\prime }}$ for other sources of uncertainty, such as loads. Therefore, this process will not be elaborated on here.

The lower-level model minimizes the operation costs of the HMGs under the forecasted scenarios and DcUS extreme scenarios as follows:

In detail, the operation costs under both these scenarios include the operation costs of the AC microgrid (ACM) $C_{s^{\prime }}^{\text{AC}}$ and the DC microgrid (DCM) $C_{s^{\prime }}^{\text{DC}}$, as well as the energy loss cost (ELC) of BdCs $c_{s^{\prime }}^{\text{L}}{\vartheta }_{s^{\prime }}$, which are outlined as follows.

where $\text{AC}\left(i\right)$ is the set of DCMs that are connected with ACM i, while $\text{DC}\left(j\right)$ is the set of ACMs that are connected with DCM j. The decision variables in the lower-level model are defined as follows: ${\vartheta }_{s^{\prime }/s}$ is ELC through BdCs in scenario $s^{\prime }/s$; $P_{i,s^{\prime }/s}^{\text{DE}}$ denotes the output of DE unit in ACM i in scenario $s^{\prime }/s$; $P_{j,s^{\prime }/s}^{\text{FC}}$ represents the output of FC unit in DCM j in scenario $s^{\prime }/s$; $P_{\mathit{ij},h,s^{\prime }/s}^{\text{acdc}}$ is the power flow of line h in corridor ij from an ACM to a DCM in scenario $s^{\prime }/s$; $P_{\mathit{ij},h,s^{\prime }/s}^{\text{dcac}}$ denotes the power flow of line h in corridor ij from a DCM to an ACM in scenario $s^{\prime }/s$; $Z_{\mathit{ij},h,s^{\prime }/s}^{\text{acdc}}$ represents the binary variables indicating whether the BdC associated with line h in corridor ij operates in the rectifier state in scenario $s^{\prime }/s$; $Z_{\mathit{ij},h,s^{\prime }/s}^{\text{dcac}}$ denotes the binary variables indicating whether the BdC associated with line h in corridor ij operates in the inverter state in scenario $s^{\prime }/s$; ${\rho }_{i,s^{\prime }/s}^{\text{ac}}$ and ${\rho }_{j,s^{\prime }/s}^{\text{dc}}$ represent the SPG curtailment of ACM i or DCM j in scenario $s^{\prime }/s$, respectively; ${\eta }_{\mathit{ij},h}^{\text{acdc}}$ and ${\eta }_{\mathit{ij},h}^{\text{dcac}}$ denote the efficiency of BdCs associated with line h in corridor ij, respectively; ${\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{acdc}}$ and ${\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{dcac}}$ represent the utilization rate of line h in corridor ij in scenario $s^{\prime }/s$, respectively. The parameters in the lower-level model are defined as follows: $a_q$ represents the coefficients of BdC efficiency, where q is the exponent parameter representing the power in polynomial expansion; $c_{i,s^{\prime }/s}^{\text{DE}}$, $c_{j,s^{\prime }/s}^{\text{FC}}$, and $c_{s^{\prime }/s}^{\text{L}}$ are the cost coefficients of DE units, FC units, and BdC energy loss in ACM i or DCM j in scenario $s^{\prime }/s$, respectively; $c^{\text{SPG}}$ represents the penalty coefficients of SPG curtailment; $P_i^{\text{DE min}}$ and $P_i^{\text{DE}\max }$ denote the lower and upper power outputs of DE unit in ACM i, respectively; $P_j^{\text{FC min}}$ and $P_j^{\text{FC}\max }$ represent the lower and upper power outputs of FC unit in DCM j, respectively; $P_{\mathit{ij},h}^{\text{R}}$ is the capacity of line h in corridor ij; ${\widehat{\lambda }}_{i,s^{\prime }}$ and ${\widehat{\lambda }}_{j,s^{\prime }}$ denote the SPG outputs of ACM i and DCM j in forecasted scenario $s^{\prime }$, respectively; $\mu$ represents the allowable ratio of RES curtailment; $P_{\text{L},i,s^{\prime }}^{\text{ac}}$ and $P_{\text{L},j,s^{\prime }}^{\text{dc}}$ indicate the load demand of ACM i and DCM j in forecasted scenario $s^{\prime }$, respectively;${\lambda }_{i,s}^{\text{AC}}$ and ${\lambda }_{j,s}^{\text{DC}}$ denote the SPG outputs of ACM i and DCM j in the sth vertex of the DcUS (which corresponds to the sth scenario), respectively; $P_{\text{L},i,s}^{\text{ac}}$ and $P_{\text{L},j,s}^{\text{dc}}$ denote the load demand of ACM i and DCM j in scenario $s$, respectively.

The operation costs for each ACM and DCM are, respectively, calculated in Eqs. (3d), (3e), both of which include the fuel costs of conventional units and penalty costs associated with RES curtailment (RC) and load-shedding demands. Eq. (3f) evaluates the energy loss of BdCs due to rectifier and inverter efficiencies, which are defined in Eqs. (3g), (3h), (3i), (3j) based on previously proposed high-order polynomial relationships with respect to ${\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{acdc}}$ and ${\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{dcac}}$ [11], [12]. As can be seen, Eqs. (3k), (3l), (3m), (3n), (3o), (3p), (3q), (3r), (3s), (3t) are applicable under forecasted scenario $s^{\prime }$. Specifically, Eqs. (3k), (3l), respectively, limit the power of BdC rectification and inverter operations, while Eq. (3m) ensures that no BdC resides in rectifier and inverter states simultaneously. Eq. (3n) ensures that $Z_{\mathit{ij},h,s^{\prime }}^{\text{acdc}}=Z_{\mathit{ij},h,s^{\prime }}^{\text{dcac}}=0$ if $I_{\mathit{ij},h}=0$. Eqs. (3o), (3p) limit the power produced by DE and FC units, respectively. Eqs. (3q), (3r) ensure that generation and load demands in forecasted scenarios are balanced at each ACM and DCM, respectively. Eqs. (3s), (3t) limit the allowable SPG curtailment in forecasted scenarios at each ACM and DCM, respectively. As can be seen, Eqs. (3u)–(3zd) applicable under DcUS extreme scenarios have similar forms to those of Eqs. (3k), (3l), (3m), (3n), (3o), (3p), (3q), (3r), (3s), (3t).

Remark 1 Different from previous data-driven stochastic approaches that focus solely on constructing uncertainty sets for the probability distributions of RES generation [29], [32], the proposed model incorporates uncertainty sets for both the probability distributions and the generation outputs of RES units. This allows the model to produce solutions that are robust against uncertainties in both aspects. Moreover, as highlighted in previous studies [33], a solution remains robust against any fluctuations in RES generation within the predefined uncertainty set if it is robust against a limited number of extreme scenarios within that set. To leverage this finding, the proposed model accounts for constraints under extreme scenarios, ensuring the robustness of the planning strategy while simultaneously reducing computational complexity.

A novel data-correlated SPG uncertainty set is developed to balance the need for security in the planning solution while minimizing unnecessary conservatism. This model accurately depicts the spatial and temporal correlations of SPG uncertainty by cutting off the unnecessary power output region in the SPG uncertainty set. Subsequently, a minimum high-dimensional uncertainty set volume based model is proposed to optimize the cutting hyperplane or the vertices of the DcUS.

The BbUS can be expressed according to the actual SPG output ${\lambda }_i$, the predicted SPG output ${\widehat{\lambda }}_i$, and the upward/downward fluctuations in SPG outputs ${\lambda }_i^{\text{u}}$/${\lambda }_i^{\text{d}}$ as follows [21], [22].

where ${\epsilon }_i^{\text{u}}$ and ${\epsilon }_i^{\text{d}}$ are binary variables, respectively, indicating the upward and downward direction of SPG output fluctuations. Accordingly, it can be inferred that the BbUS directly bounds the predicted SPG outputs. As a result, a BbUS typically provides the most conservative planning solutions. This is addressed to some extent in the BcUS by excluding unlikely scenarios from the BbUS via the application of an uncertainty budget $\mathit{\Gamma }$, detailed as below [23].

Fig. 3(a) presents the historical SPG output data, predicted data, and BbUS and BcUS regions corresponding to two SPG clustering units. Note that the number of vertices in both the BbUS and BcUS is 2^I, where $I$ denotes the number of SPG clustering units. Additionally, the number of SPG clustering units is equal to the number of microgrids, as each microgrid contains a group of SPG units. As can be seen, the application of $\mathit{\Gamma }$ in the BcUC truncates the BbUS by omitting the less likely SPG outputs located in the upper right and lower left corners of the historical dataset, which can therefore reduce the conservativeness of the obtained planning strategy to some extent. However, this method of controlling BcUS conservativeness neglects the spatiotemporal correlations of SPG output uncertainties, which can produce non-robust solutions for scenarios that fall outside of the BcUS.

By contrast, the CHUS can be expressed as follows:

where $m$ is the index representing the set ${\mathcal{M}}_{\text{CH}}$, which is defined as the set of vertices, and ${\epsilon }_{\mathit{im}}$ denotes a continuous variable within the interval [0, 1]. Fig. 3(b) illustrates the CHUS region corresponding to the historical SPG output data and predicted data presented for the two SPG clustering units in Fig. 3(a). In contrast to the 2^I vertices observed for the BbUS, the number of vertices in the CHUS depends heavily on the historical SPG output data applied in its construction.

The proposed DcUS formulation for applying hyperplane boundary constraints to the BbUS is given as follows:

where $u$ is the index representing the vertices of the BbUS, and ${\kappa }_{u1}$, ${\kappa }_{u2}$, … , ${\kappa }_{\mathit{uI}}$ are the corresponding hyperplane coefficients applied for the uth vertex. The vertices where all solar power units reach their minimum and maximum output values are designated as the first vertex and the last (i.e., the 2^Ith) vertex, respectively. Note that each vertex of the BbUS corresponds to a hyperplane constraint, with the exception of the 1st and 2^Ith vertices. Thus, the number of constraints applied in Eq. (7a) is 2^I − 2. For convenience of expression, the hyperplane constraint in Eq. (7a) is rewritten in matrix form by Eq. (7b).

where ${\mathit{\kappa }}_u^{\text{T}}$ represents the coefficient matrix of the SPG output vector $\mathit{\lambda }$. Fig. 3(b) illustrates the corresponding DcUS region. As can be seen, the DcUS region substantially truncates the BbUS region by omitting unlikely scenarios not represented in the historical data, as observed for the CHUS region, and therefore reduces the conservativeness of planning solutions. However, the DcUS includes many fewer vertices than the CHUS and therefore offers enhanced solution efficiency.

Fig. 4 illustrates the process applied for selecting the hyperplanes that minimize the volume of the DcUS while enclosing all historical data points. This process can be reformulated as maximizing the remaining area of the BbUS given by the red shaded area in Fig. 4. In addition, each hyperplane in the DcUS must correspond to a vertex in the BbUS because the proposed DcUS is a subset of the BbUS. Thus, the objective function applied for selecting the hyperplanes can be formulated as follows:

where $d_{\mathit{ui}}$ is the geometric distance between the uth vertex ${\lambda }_u^{\text{B}}$ of the BbUS and the vth vertex ${\lambda }_v^{\text{Dc}}$ of the DcUS. For each vertex u $=2$, 3, …, 2^I − 1 in the BbUS, the corresponding adjacent vertex index v in the DcUS can be expressed in the following manner.

The geometric distances between vertices are defined as follows. The geometric distance between ${\lambda }_u^{\text{B}}$ and ${\lambda }_v^{\text{Dc}}$ is set to zero across all SPG clustering units, except in the ith SPG clustering unit, where the geometric distance is computed in the manner described below:

Here, ${\gamma }_{\mathit{ui}}$ is a binary directional parameter with a possible value of either 1 or −1, where a value of ${\gamma }_{\mathit{ui}}=1$ indicates that ${\lambda }_u^{\text{B}}$ resides at the upper boundary of the BbUS in the ith SPG clustering unit, and a value of ${\gamma }_{\mathit{ui}}=-1$ indicates that ${\mathit{\lambda }}_u^{\text{B}}$ is located at the lower boundary of the BbUS in the ith SPG clustering unit.

The constraints of the proposed uncertainty set model primarily consider the geometrical position relationships between the hyperplanes, the vertices of the BbUS and DcUS, and the historical SPG output data. These constraints can be expressed as follows:

Here, $\mathcal{R}$ and r refer to the set of all historical SPG output data and the corresponding index, respectively, and ${\mathcal{R}}^{\text{DcUS}}$ and ${\mathcal{R}}^{\text{BbUS}}$ denote the vertex sets of the DcUS and BbUS, respectively. In the above formulation, Eq. (8e) ensures that the vertices of the proposed DcUS adjacent to the vth vertex of the BbUS are on the hyperplane, while Eqs. (8b), (8c) restrict the historical SPG output data and the vertices of the BbUS to be located on different sides of the hyperplane. Eq. (8h) ensures that the coefficients of the hyperplane are not constantly zero.

However, the aforementioned uncertainty set model can represent a substantial computational burden for large HMG systems. This issue is mitigated by applying the following two modifications.

First, the volume of each resulting truncated polyhedron, along with the constraints delineated in Eqs. (8a), (8b), (8c), (8d), are independent of each other. Thus, the optimization process in Eq. (8) can be decomposed into 2^I − 2 smaller nonlinear optimization processes, which can be implemented independently in parallel. To this end, objective function Eq. (8a) is reformulated for each vertex $u=2,3,\dots ,2^I-1$, as shown in Eq. (9a).

In addition, the computational burden can be further reduced by applying the logarithm within objective function Eq. (9a) and omitting the constant term ($-\ln I!$), which yields the following.

Second, each historical SPG output data point generates a constraint in Eq. (8f). Therefore, the computational burden is reduced by replacing the vast amount of historical SPG output data in Eq. (8f) with the following constraint applied for the vertices of the CHUS.

The proposed substitution is formally validated according to the following proposition.

Proposition 1 Eq. (8f) will hold true for all historical SPG output data as long as all vertices of the CHUS satisfy Eq. (9c).

The detailed proof is provided in Appendix A. Based on this proposition, it can be deduced that the hyperplane constraints in the DcUS can be determined for each vertex $u=2,3,\dots ,2^I-1$ by solving the following nonlinear minimum uncertainty set volume-based programming model.

The proposition relies on the relationship between Eq. (8f) and the vertices of the CHUS. By ensuring that all vertices of the CHUS satisfy Eq. (9c), it guarantees that the entire convex hull (and thus any point within it) satisfies the condition. In practice, incorporating Eq. (8f) for all historical data into the uncertainty set optimization model can be computationally expensive. To address this, the proposition simplifies the model without sacrificing accuracy by considering only Eq. (9c) for the vertices of the CHUS, which is more computationally efficient since the number of vertices is significantly smaller compared to the vast amount of SPG output data.

It must be noted that the lower-level model in the proposed tri-level framework defined in Section 2 does not satisfy the strong duality condition due to the presence of highly nonlinear equations in Eqs. (3g), (3h), in addition to binary variables. Thus, the proposed model cannot be effectively solved using existing C&CG algorithms. As discussed, the highly nonlinear terms in Eqs. (3g), (3h) are first linearized using the LSA method. Then, the binary variables are removed by applying an exact convex relaxation of the lower-level model. Finally, an FPC&CG algorithm is developed based on the unique structure of the proposed planning model.

Note that a direct linearization of the nonlinear terms ${\eta }_{\mathit{ij},h}^{\text{acdc}}\left({\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{acdc}}\right)$/${\eta }_{\mathit{ij},h}^{\text{dcac}}\left({\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{dcac}}\right)$ in Eqs. (3g), (3h) leaves Eq. (3f) nonlinear. Thus, Eq. (3f) is linearized as a whole. The two terms in Eq. (3f) have similar expressions. Here, the proposed linearization method is illustrated in Fig. 5 for the linearization of $P_{\mathit{ij},h,s^{\prime }/s}^{\text{acdc}}\left(1-{\eta }_{\mathit{ij},h}^{\text{acdc}}\left({\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{acdc}}\right)\right)$ as an example. The other term $P_{\mathit{ij},h,s^{\prime }/s}^{\text{dcac}}\left(1-{\eta }_{\mathit{ij},h}^{\text{dcac}}\left({\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{dcac}}\right)\right)$ in Eq. (3f) can be linearized in a similar manner.

By defining a function $g\left(x\right)=x\left(1-{\eta }_{\mathit{ij},h}^{\text{acdc}}\left(x\right)\right)$, where $x$ corresponds to the variable ${\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{dcac}}$, $g\left(x\right)$ belongs to the space of continuous functions over the interval $\left[0,1\right]$ and $g\left(x\right)$ possesses a continuous first-order derivative over the interval [0, 1], $P_{\mathit{ij},h,s^{\prime }/s}^{\text{acdc}}\left(1-{\eta }_{\mathit{ij},h}^{\text{acdc}}\left({\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{acdc}}\right)\right)$ can be rewritten as follows.

Then, we seek to identify a linear Legendre expression $O\left(x\right)={\sum }_{y=0}^1{o}_y{\varphi }_y\left(x\right)$ (${\varphi }_y\left(x\right)=x^y$) that can be applied as a substitute for $g\left({\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{acdc}}\right)$. This is achieved by minimizing the approximation error $\delta$ over the interval [0, 1], as presented below:

where y is an index used in the polynomial expansion $O\left(x\right)$.

Because $g\left(x\right)$ possesses a continuous first-order derivative over the interval [0, 1], the following necessary conditions are obtained for achieving the extreme value of the multivariate function $\delta$:

where c is an index that can take values 0 or 1. This can be rewritten in Eq. (11c).

For convenience of expression, the inner product of ${\varphi }_c$, ${\varphi }_y$, and $g$ over the interval $\left[0,1\right]$ is defined as follows.

This yields the following Hilbert matrix.

Finally, Eq. (11c) can be reformulated in the following manner.

Here, $\mathit{O}={\left[o_0{o}_1\right]}^{\text{T}}$ and $\mathit{g}={\left[\left(g,{\varphi }_0\right)\left(g,{\varphi }_1\right)\right]}^{\text{T}}$. Because Cramer’s determinant $\left|{\mathit{H}}_{\text{b}}\right|\ne 0$, ${\varphi }_0$ and ${\varphi }_1$ are linearly independent. Therefore, a unique solution ${\mathit{O}}^{\ast }$ to Eq. (11g) can be found, which yields the following linear expression of $g\left(x\right)$.

As indicated in Eqs. (3g), (3h), the coefficient $a_q$ for the BdC-based transfer of power from ACMs to DCMs is identical to that of the power transferred from DCMs to ACMs. This yields the following outcome.

However, it is worth noting that a similar linearization process can be applied to derive the linear form of $P_{\mathit{ij},h}^{\text{R}}g\left({\varpi }_{\mathit{ij},h,s^{\prime }/s}^{\text{dcac}}\right)$ when the same coefficient $a_q$ cannot be applied to both the rectifier and inverter operations of a BdC.

Finally, a mixed integer linear programming (MILP) version of the lower-level model can be developed by substituting Eqs. (3g), (3h) with Eqs. (11e), (11f), respectively.

While the lower-level model in its MILP form could be solved by a nested C&CG algorithm, the solution process is computationally demanding. To address this issue, the following convex formulation of the lower-level model is presented.

Then, the solution of this convex formulation is formally proven to be exactly equivalent to the solution of the lower-level model in its MILP form according to the following proposition.

Proposition 2 Given the value of $I_{\text{BDC},ij,h}$, a solution to the lower level optimization Eq. (12) is sub-optimal in the case when $P_{\mathit{ij},h,s^{\prime }}^{\text{acdc}}{P}_{\mathit{ij},h,s^{\prime }}^{\text{dcac}}>0$, as long as $o_{\mathit{ij},h}^{\ast 1}>0$.

The detailed proof is presented in Appendix A. Intuitively, the condition $o_{\mathit{ij},h}^{\ast 1}>0$ indicates that the energy loss of a BdC is positively correlated with its rectifier/inverter power $P_{\mathit{ij},h,s^{\prime }}^{\text{acdc}}/P_{\mathit{ij},h,s^{\prime }}^{\text{dcac}}$. Thus, a sub-optimal solution must be obtained when a BdC operates in both its rectifier and inverter modes simultaneously because unnecessary BdC energy loss occurs. In addition, assuming that the lower-level model is in its linear form yields the following proposition.

Proposition 3 Lower-level model Eq. (12) is robust against any SPG output fluctuations enclosed within the DcUS as long as it is robust against the extreme scenarios or vertices of the DcUS.

The detailed proof is presented in Appendix A. This proposition indicates that applying the DcUS ensures that the proposed planning model is fully robust against the worst-case realization of SPG output fluctuations in the standard BbUS [13], but with much less conservativeness.

The DcUS-based tri-level robust BdC planning model obtained after applying LSA and convex relaxation can be compactly reformulated in the following matrix form.

where $\mathit{a}$, ${\mathit{b}}_1$, ${\mathit{b}}_2$, $\mathit{C}$, $\mathit{c}$, ${\mathit{A}}_1$, ${\mathit{A}}_2$, and ${\mathit{A}}_3$ are coefficient matrices or vectors corresponding to the tri-level robust model.

The middle and lower level models in Eq. (13) can be merged into a single-level bi-layer model using duality theory. However, this single-level model can be excessively large, particularly when the number of scenarios is substantial. In addition, the middle-level model contains second-order cone constraints that can potentially complicate the solution process. This issue is addressed by recasting tri-level model Eq. (13) into a master–slave problem within the advanced FPC&CG algorithm to solve the model in parallel without the need for applying duality theory or Karush–Kuhn–Tucker (KKT) conditions.

(1) Slave problem. A bi-layer maximization–minimization process is applied within the slave model to identify the worst-case scenario in the uncertainty set under a specific set of investment decision variables $\mathit{I}$. This yields the following model.

Given that ${\mathcal{S}}^{\text{F}}\cap {\mathcal{S}}^{\text{Dc}}=\mathrm{\varnothing }$ and that the feasible region of ${\mathit{x}}_{s^{\prime }}$ is separate from that of ${\mathit{x}}_s$, objective function Eq. (14a) is reformulated in Eq. (15a).

Note that the optimal value of the term $\underset{\mathcal{P}}{\text{max}}\underset{\mathcal{Y}}{\text{min}}{\sum }_{s\in {\mathcal{S}}^{\text{Dc}}}{\mathit{b}}_2^{\text{T}}{\mathit{x}}_s$ in Eq. (15a) is independent of the maximization process, as the variables $p_{s^{\prime }}\in \mathcal{P}$ in the maximization process are not involved in any terms included in the minimization process. Hence, Eq. (15a) can be simplified as follows:

With regard to the first term ${\max }_{\mathcal{P}}{\min }_{\mathcal{Y}}{\sum }_{s^{\prime }\in {\mathcal{S}}^{\text{F}}}{p}_{s^{\prime }}{\mathit{b}}_{\text{1}}^{\text{T}}{\mathit{x}}_{s^{\prime }}$ in Eq. (15b), the fact that its feasible region is disjoint with $\mathcal{P}$ enables Eq. (15b) to be further restructured as follows:

In addition, the two large-scale minimization problems in Eq. (15c) can be decomposed into $(\left|{\mathcal{S}}^{\text{F}}\right|+\left|{\mathcal{S}}^{\text{Dc}}\right|)$ small independent optimization models because the operational constraints of the small models under each forecasted $s^{\prime }$ and the operational constraints of the small models under extreme scenario $s$ are disjoint from each other.

The parallel strategy applied for solving the $(\left|{\mathcal{S}}^{\text{F}}\right|+\left|{\mathcal{S}}^{\text{Dc}}\right|)$ small linear programming models for $s^{\prime }\in {\mathcal{S}}^{\text{F}}$ and $s\in {\mathcal{S}}^{\text{Dc}}$ is given as follows:

After obtaining the optimal solution $\left(y_{s^{\prime }}^{\ast },y_s^{\ast }\right)$ for Eq. (16), Eq. (15c) is reformulated into the following second-order cone programming model:

The optimal solution of the slave problem obtained in the kth iteration $\left(p_{s^{\prime }}^{k\ast },{\mathit{x}}_{s^{\prime }}^{k\ast },{\mathit{x}}_s^{k\ast }\right)$ provides the upper bound (UB) for model Eq. (13), where $k=1,2,\dots ,K$. Subsequently, a set of slave problem constraints associated with variables $({\mathit{x}}_{s^{\prime }}^{k+1\ast },{\mathit{x}}_s^{k+1\ast })$ are generated by fixing the optimal probability value $p_{s^{\prime }}^{k\ast }$ ($s^{\prime }\in {\mathcal{S}}^{\text{F}}$) in the original slave problem Eq. (14). These constraints are then added to the master problem with a dummy continuous variable $\eta$ as follows:

(2) Master problem. The master problem serves as a relaxation of the original problem and thereby provides the lower bound (LB) for model Eq. (13). Following the kth iteration of the C&CG algorithm, the master model is formulated as the following MILP.

Finally, the master problem and slave problem are solved iteratively until the difference between the upper and lower bounds satisfies the optimality threshold $\epsilon$. The proposed FPC&CG algorithm is summarized in Algorithm 1.

Remark 2 The FPC&CG algorithm ensures the optimality of the obtained solutions by iteratively solving the master problem and subproblems, progressively refining the feasible region and objective value. At each iteration, the algorithm incorporates a new, previously unconsidered critical uncertainty scenario into the restricted master problem and updates the UB and LB of the objective value. Since the vertices of the uncertainty set $\mathcal{P}$ is finite, the algorithm converges after a finite number of steps. Convergence is achieved when the gap between the bounds falls below a predefined tolerance, and no newly identified adversarial scenario significantly degrades the solution, ensuring that the final decision is robustly optimal across all possible uncertainty realizations. Compared to typical C&CG algorithms [34], [35], the developed FPC&CG algorithm achieves faster convergence and guarantees the optimality of the solutions by leveraging parallelization and convex relaxation techniques.

The developed solution algorithm has the potential to effectively handle large test cases with thousands of variables by combining the LSA method and the FPC&CG algorithm. The LSA simplifies the complex efficiency expressions of BdCs into computationally efficient approximations, reducing the computational burden and enabling scalability. Meanwhile, the FPC&CG algorithm dynamically generates only relevant variables and constraints during the optimization process, significantly reducing problem size and computational overhead. Its fully parallel design allows simultaneous processing, further accelerating computation. By tailoring both methods to exploit problem-specific structures, your approach ensures efficient, scalable, and accurate optimization for large-scale BdC systems.

The performance of the proposed BdC-based interconnection planning method was evaluated by its application to a numerical case study involving the practical HMG system illustrated in Fig. 1. The system parameters have been reported in Ref. [36]. All numerical computations were executed using General Algebraic Modeling System (GAMS) software operating on a laptop computer outfitted with an Intel® i5-3230M central processing unit (CPU) and 4 GB random access memory (RAM). Nonlinear programming model Eq. (10) applied for determining the hyperplane constraints in the DcUS was solved using the Interior Point OPTimizer (IPOPT) solver, and the MILP model was solved using the CPLEX solver. The convergence thresholds for both these solvers and the optimality threshold $\epsilon$ of the proposed FPC&CG algorithm were all set to 10⁻⁶, ensuring high precision, consistency with solver tolerances and computational efficiency. The penalty coefficients for load shedding (LS) and RES output curtailment were both set to 1.5 USD∙(kW∙h)⁻¹.

The performance of the proposed method was evaluated from three perspectives, namely: ① An analysis of the error associated with the proposed LSA method for approximating BdC efficiency, compared to the error produced by a previously proposed approximation method [10], [15] that regards BdC efficiency ${\eta }_{\mathit{ij},h}^{\text{acdc}}$/${\eta }_{\mathit{ij},h}^{\text{dcac}}$ as a constant. ② The convex relaxation performance of the lower-level model. Comparisons of the costs, interconnection configurations, and robustness obtained using the BbUS [22], CHUS [24], and the employed DcUS within the tri-level robust model. ③ Additionally, results from the previously proposed advanced distributionally robust planning (ADRP) method [29], [37], [38] are presented, which seeks to identify an optimal solution that is robust against the uncertainty of SPG output probability distributions but without considering the uncertainty of SPG output scenarios.

The error associated with the proposed LSA method and the previously proposed method was analyzed based on the absolute error (AE), the relative error (RE), and the average error ($\overline{A}$), which are defined as follows.

The AE and RE results obtained for the proposed method and the previously proposed method [10], [15] are presented in Fig. 6 as functions of the BdC utilization rate.

As can be seen, the RE and AE values obtained using the proposed LSA method are generally less than those of the previously proposed method. The superior estimation performance of the LSA method is particularly apparent when comparing the average error of 9.1% obtained for this method with the average error of 33.4% obtained for the existing method.

The convex relaxation performance of lower-level model Eq. (12) was evaluated by examining the hourly violations in simultaneous rectifier and inverter operations by BdCs observed across different SPG output scenario sizes, namely, one day, one week, one month, and one year, according to the convex relaxation gap (CRG) defined as follows.

The results show that the CRG values are uniformly zero for all scenario sizes. That is, no violations in simultaneous rectifier and inverter operations were observed for BdCs over all periods considered. This can be attributed to the fact that the condition $o_{\mathit{ij},h}^{\ast 1}=0.01586>0$ was obtained, which strictly adheres to the condition established in Proposition 2 that the relaxation in Eq. (12) is exact as long as $o_{\mathit{ij},h}^{\ast 1}>0$.

The area and number of vertices of the BbUS, CHUS, and proposed DcUS obtained for two SPG units were compared based on 8760 h of historical data. The results are presented in Fig. 7.

As expected, the BbUS generates the largest uncertainty set with the least number of vertices, while the CHUS produces the smallest uncertainty set size with the greatest number of vertices. Specifically, the area of the CHUS is reduced by 78.8% relative to that of the BbUS, but the number of vertices is substantially increased from 4 to 15. By contrast, the proposed DcUS presents a striking balance between the area of the uncertainty set and its number of vertices, where the area of the DcUS increases by only 13.2% relative to that of the CHUS, but the number of vertices decreases from 15 to 6.

The BdC interconnection planning results obtained for the system illustrated in Fig. 1 using the proposed method with the different uncertainty sets and the previously proposed ADRP method are presented in Fig. 8. As can be seen, the plan developed by the existing ADRP method includes the smallest number of 10 newly built BdCs, while the BbUS solution results in the largest number of 14 newly built BdCs. By contrast, the CHUS solution introduces 12 newly built BdCs, and the proposed DcUS solution introduces 13 newly built BdCs. Thus, both approaches have alleviated the overly conservative solution obtained under the BbUS to some extent. It is noteworthy that most newly installed BdCs are located within the same HMG rather than across different HMGs. This arises because installing a BdC across different HMGs incurs a higher cost than doing so within the same HMG. In addition, the computation times required to obtain the BdC planning schemes using the existing ADRP method and the proposed method with the BbUS, CHUS, and DcUS were 681, 733, 6098, and 942 s, respectively. Noticeably, the use of the DcUS substantially reduced the computation time relative to that of the CHUS owing to its greatly reduced number of vertices.