Non-Iterative Solution for Coordinated Optimal Dispatch via Equivalent Projection-Part II: Method and Applications

TL;DR

This paper introduces a non-iterative coordinated optimal dispatch framework using equivalent projection theory, featuring a novel projection algorithm with explicit error guarantees, applicable to multi-area and transmission-distribution systems.

Contribution

It develops a new projection algorithm with error control for equivalent projection models, enabling efficient non-iterative optimization in power system dispatch.

Findings

Algorithm achieves conservative approximation within error tolerance.

Method improves computational efficiency and scalability.

Validated effectiveness through diverse case studies.

Abstract

This two-part paper develops a non-iterative coordinated optimal dispatch framework, i.e., free of iterative information exchange, via the innovation of the equivalent projection (EP) theory. The EP eliminates internal variables from technical and economic operation constraints of the subsystem and obtains an equivalent model with reduced scale, which is the key to the non-iterative coordinated optimization. In Part II of this paper, a novel projection algorithm with the explicit error guarantee measured by the Hausdorff distance is proposed, which characterizes the EP model by the convex hull of its vertices. This algorithm is proven to yield a conservative approximation within the prespecified error tolerance and can obtain the exact EP model if the error tolerance is set to zero, which provides flexibility to balance the computation accuracy and effort. Applications of the EP-based coordinated dispatch are demonstrated based on the multi-area coordination and transmission-distribution coordination. Case studies with a wide range of system scales verify the superiority of the proposed projection algorithm in terms of computational efficiency and scalability, and validate the effectiveness of the EP-based coordinated dispatch in comparison with the joint optimization.

Tables

Table 1. TABLE II: Computational Efficiency of EP-based MACOD

System	Time of joint optimization (s)	Time of EP-based coordination (s)
		System reduction	Coordinated optimization	Total
$20\times \text{SG-200}$	0.42	0.38	0.11	0.49
$40\times \text{SG-200}$	0.61	0.32	0.10	0.42
$20\times \text{SG-500}$	0.51	0.49	0.09	0.58
$40\times \text{SG-500}$	0.78	0.51	0.10	0.61

Table 2. TABLE III: Computational Efficiency of EP Calculation

System	Number of DERs	Time of FME (s)	Time of PVE (s)	Rate of model reduction
DN-13	6	3.75	0.34	98.4%
DN-25	12	41.29	0.32	99.5%
DN-49	24	>1200	0.37	>99.9%
DN-241	120	>1200	0.62	>99.9%
DN-2401	1200	>1200	2.37	>99.9%

Table 3. TABLE IV: Computational Efficiency of EP-based TDCOD

Num of DNs	Time of joint optimization (s)	Time of EP-based coordination (s)
		System reduction	Coordinated optimization	Total
$50$	0.26	0.31	0.08	0.39
$100$	0.47	0.30	0.09	0.39
$200$	0.68	0.31	0.12	0.43
$400$	0.97	0.31	0.20	0.51

Equations

$$\Omega :=\left\{(x,y)\in {\mathbb{R}}^{N_x}\times {\mathbb{R}}^{N_y}:Ax+By\le c\right\}.$$

$$\Phi :=\left\{x\in {\mathbb{R}}^{N_x}:\exists y,\text{s.t.}(x,y)\in \Omega \right\}.$$

$$\Phi =\left\{x\in {\mathbb{R}}^{N_x}:\tilde{A}x\le \tilde{c}\right\}.$$

$$\begin{array}{rl}\Phi =& \text{conv}(V)\\ :=& \left\{\sum _{i=1}^{N_v}{\lambda }_i{\hat{x}}_i:{\hat{x}}_i\in V,{\lambda }_i\ge 0,\sum _{i=1}^{N_v}{\lambda }_i=1\right\}.\end{array}$$

$$\begin{array}{rl}\underset{x,y}{\max }& h={\alpha }^{\top }x\\ \text{s.t.}& Ax+By\le c\end{array}$$

$$\Delta h^{(k,j)}=\frac{({\tilde{A}}_j^{(k)}{)}^{\top }{\hat{x}}^{(k,j)}-{\tilde{d}}^{(k,j)}}{\Vert {\tilde{A}}_j^{(k)}\Vert }.$$

$$D^{(k)}:=\underset{x_1\in \Phi }{\max }\underset{x_2\in {\hat{\Phi }}^{(k)}}{\min }\Vert x_2-x_1{\Vert }_2.$$

$$\begin{array}{rl}{D}^{(k)}& =\underset{x_1\in V}{\max }\underset{x_2\in {\hat{\Phi }}^{(k)}}{\min }\Vert x_2-x_1{\Vert }_2\\ & =\underset{j}{\max }\Delta h^{(k,j)}\\ & =\max H^{(k)}.\end{array}$$

$$\min {\pi }^{MA}=\sum _{r\in {\mathcal{R}}^{RN}}{\pi }_r^{RN}$$

$${\overline{\pi }}_r^{RN}\ge {\pi }_r^{RN}\ge \sum _{n\in {\mathcal{N}}_r^{RN}}{\pi }_{r,n}^{RN,G},$$

$${\pi }_{r,n}^{RN,G}\ge a_{r,n,i}^{RN,G}{P}_{r,n}^{RN,G}+b_{r,n,i}^{RN,G},\forall n,i,$$

$$\sum _{n\in {\mathcal{N}}_r^{RN}}{P}_{r,n}^{RN,G}=\sum _{n\in {\mathcal{N}}_r^{RN}}{P}_{r,n}^{RN,D}+\sum _{n\in {\mathcal{N}}_r^{RN,B}}{P}_{r,n}^{RN,B},$$

$$\begin{array}{rl}-{\overline{F}}_{r,l}^{RN}\le & \sum _{n\in {\mathcal{N}}_r^{RN}}{T}_{r,l,n}{P}_{r,n}^{RN,G}-\sum _{n\in {\mathcal{N}}_r^{RN}}{T}_{r,l,n}{P}_{r,n}^{RN,D}\\ & -\sum _{n\in {\mathcal{N}}_r^{RN,B}}{T}_{r,l,n}{P}_{r,n}^{RN,B}\le {\overline{F}}_{r,l}^{RN},\forall l,\end{array}$$

$${\underline{P}}_{r,n}^{RN,G}\le P_{r,n}^{RN,G}\le {\overline{P}}_{r,n}^{RN,G},\forall n.$$

$$P_s^{TL}=\frac{1}{x_s}\left({\theta }_{r_s^F}-{\theta }_{r_s^T}\right),\forall s,$$

$$P_{r,n}^{RN,B}=\sum _{\begin{array}{c}{r}_s^F=r,\\ n_s^F=n\end{array}}{P}_s^{TL}-\sum _{\begin{array}{c}{r}_s^T=r,\\ n_s^T=n\end{array}}{P}_s^{TL},\forall n,r,$$

$${\underline{F}}_s^{TL}\le P_s^{TL}\le {\overline{F}}_s^{TL},\forall s.$$

$${\Omega }_r^{RN}=\left\{(P_r^{RN,B},{\pi }_r^{RN},P_r^{RN,G}):\text{Eq.}\text{eqref}macd:cs{t}_{a}rea\right\}.$$

$$\begin{array}{rl}{\Phi }_r^{RN}=\{& (P_r^{RN,B},{\pi }_r^{RN}):\exists P_r^{RN,G},\\ & \text{s.t. }(P_r^{RN,B},{\pi }_r^{RN},P_r^{RN,G})\in {\Omega }_r^{RN}\}.\end{array}$$

$$\min {\pi }^{TD}=\sum _{n\in {\mathcal{N}}^{TN}}{C}_n(P_n^{TN,G})+\sum _{r\in {\mathcal{R}}^{DN}}{\pi }_r^{DN}$$

$$\sum _{n\in {\mathcal{N}}^{TN}}{P}_n^{TN,G}+\sum _{r\in {\mathcal{R}}^{DN}}{P}_{r,0}^{DN}=\sum _{n\in {\mathcal{N}}^{TN}}{P}_n^{TN,D},$$

$$\begin{array}{rl}-{\overline{F}}_l^{TN}\le & \sum _{n\in {\mathcal{N}}^{TN}}{T}_{l,n}^{TN}\left(P_n^{TN,G}-P_n^{TN,D}\right)\\ & +\sum _{r\in {\mathcal{R}}^{DN}}{T}_{l,r}^{TN}{P}_{r,0}^{DN}\le {\overline{F}}_l^{TN},\forall l,\end{array}$$

$${\underline{P}}_n^{TN,G}\le P_n^{TN,G}\le {\overline{P}}_n^{TN,G},\forall n.$$

$${\overline{\pi }}_r^{DN}\ge {\pi }_r^{DN}\ge \sum _{n\in {\mathcal{N}}_r^{DN}}{\pi }_{r,n}^{DN,G},$$

$${\pi }_{r,n}^{DN,G}\ge a_{r,n,i}^{DN,G}{P}_{r,n}^{DN,G}+b_{r,n,i}^{DN,G},$$

$$\begin{array}{rl}& P_{r,n}^{DN,G}+\sum _{m\in {\mathcal{A}}_n}{P}_{r,mn}^{DN,F}=P_{r,n}^{DN,D},P_{r,0}^{DN}=\sum _{m\in {\mathcal{A}}_0}{P}_{r,m0}^{DN,F},\\ & Q_{r,n}^{DN,G}+\sum _{m\in {\mathcal{A}}_n}{Q}_{r,mn}^{DN,F}=Q_{r,n}^{DN,D},Q_{r,0}^{DN}=\sum _{m\in {\mathcal{A}}_0}{Q}_{r,m0}^{DN,F},\end{array}$$

$$V_{r,m}^2-V_{r,n}^2=2(r_{r,mn}{P}_{r,mn}^{DN,F}+x_{r,mn}{Q}_{r,mn}^{DN,F}),$$

$$(\cos \frac{2k\pi }{N_s})P_{r,mn}^{DN,F}+(\sin \frac{2k\pi }{N_s})Q_{r,mn}^{DN,F}\le (\cos \frac{\pi }{N_s}){\overline{F}}_{r,mn}^{DN},$$

$${\underline{V}}_{r,n}^2\le V_{r,n}^2\le {\overline{V}}_{r,n}^2,$$

$${\underline{P}}_{r,n}^{DN,G}\le P_{r,n}^{DN,G}\le {\overline{P}}_{r,n}^{DN,G},{\underline{Q}}_{r,n}^{DN,G}\le Q_{r,n}^{DN,G}\le {\overline{Q}}_{r,n}^{DN,G}.$$

Full text

I Introduction

I-A Background and Literature Review

The coordinated optimal dispatch (COD) of power systems in different regions, voltage levels, and communities is vital for improving the operational economy and security. Conventional coordinated optimization methods rely on iterative information exchange among subsystems, which results in drawbacks including the convergence issue, communication burden, scalability issue, and incompatibility with the serial coordination scheme in practice. In this regard, realizing the COD in a non-iterative fashion, i.e., without iterative information exchange among subsystems, is of great interest to both academia and industry. This two-part paper aims at this issue and proposes the equivalent projection (EP)-based solution. The basic idea is to make external equivalence of technical and economic features of the subsystem for the upper-level coordinated optimization, which requires only a single-round exchange of some boundary information. Following the theory and framework introduced in Part I, the calculation method and typical applications of the EP will be discussed in this present paper.

The EP is mathematically a geometric projection problem, which eliminates internal decision variables from the secure and economic operation constraints of the subsystem and yields a low-dimensional region regarding coordination variables to represent the system. Projection calculation is a challenging task even for linear systems due to the exponential complexity in worst-case scenarios. The most classic projection algorithm is the Fourier-Motzkin Elimination (FME). The FME eliminates variables one by one through linear combinations of inequalities with opposite signs before the coefficient of each variable to be eliminated [1]. However, numerous redundant constraints will be generated after eliminating each variable, which not only complicates the projection result, but also increases the number of constraints needed to be processed for eliminating the next variable and thus aggravates the computation burden dramatically [2]. To over these drawbacks, existing studies propose the redundancy identification technique to accelerate the FME calculation and obtain the minimal representation of the projection result [3, 4]. Another basic algorithm for polyhedral projection is the block elimination [5], which identifies facets of the projected polytope based on the projection lemma. Though one-by-one elimination of variables is avoided, the block elimination may also generate redundant constraints. The polyhedral projection problem is theoretically proven to be equivalent to the multi-parametric programming (MPP) problem [6], which enables the solution of the projection problem via MPP algorithms. This method also leverages the projection lemma to identify irredundant inequalities of the projection.

The application of above projection methods to the EP calculation in the power system is limited due to features of the coordinated dispatch problem. First, these methods have exponential complexity to the number of variables to be eliminated and thus, are inefficient for power system dispatch problems where the number of internal variables is large. Second, a conservative approximation of the exact EP model is usually desired in engineering applications to reduce the computation effort. However, the aforementioned methods cannot yield an inner approximation since they characterize the projection by generating inequalities and will overestimate the projection result when terminated prematurely without finding all the inequalities. In this regard, existing studies in the power system literature propose different models to approximate the projection region, e.g., the box model [7], the Zonotope model [8], the ellipse model [9], and the robust optimization-based model [10]. These methods approximate the projection region with simple shapes at the cost of accuracy. As the approximation shape is pre-fixed, these methods are only applicable to projection regions with specific geometric structures and the approximation accuracy cannot be adjusted flexibly.

Another critical issue regarding the EP calculation is the error metric for controlling the calculation accuracy of the algorithm. The error metric for the EP model measures the difference between the approximated region and the exact one. Reference [8] uses the distances between parallel facets of the region to measure the approximation quality. However, this metric can only be applied to a certain kind of polytope namely the Zonotope. The volume is a more general metric for comparing two regions. Reference [11] measures the calculation accuracy of the active and reactive flexibility region of the distribution network by comparing areas of regions. Reference [11] uses the same metric to measure the flexibility region enlargement brought by power flow routers. Based on the volume metric, the Jaccard similarity can be used to measure the relative error of the approximation [12]. However, the evaluation of volume metrics is time-consuming, especially for high-dimensional problems. Furthermore, The evaluation of the volume-based error metric also relies on the knowledge of the exact projected region. Hence, this metric can only be used to test the approximation quality when the ground truth of the projection is known, but can hardly be used to control the approximation accuracy in the projection algorithm.

I-B Contributions and Paper Organization

To meet practical requirements of the EP calculation for power system applications, the Part II of this paper proposes a novel projection algorithm with an explicit accuracy guarantee. With this algorithm, the EP theory and the non-iterative COD method are instantiated for multi-area system coordination and transmission-distribution coordination. The contribution of the Part II is threefold,

1. The progressive vertex enumeration (PVE) algorithm is proposed for the EP calculation in power system applications. Compared to existing methods that solve the projection problem by eliminating internal variables, the proposed method directly identifies the projection region in the lower-dimensional coordination space and thus, has lower computational complexity. Compared to existing approximation methods, the proposed algorithm is proven to be accurate for the polyhedral projection problem. The proposed method can also be used to get an inner approximation with an explicit and adjustable error tolerance according to practical requirements.

2. The Hausdorff distance is employed to measure the approximation error of the PVE algorithm. Compared to the volume metric and the Jaccard similarity, the proposed error metric is a byproduct of the vertex identification in each step of the algorithm and thus, is computationally efficient and can be used to balance the accuracy and computation effort of the algorithm.

3. The non-iterative COD for the multi-area system and transmission-distribution system is realized based on the EP calculation, which thoroughly overcomes the disadvantages of conventional coordinated optimization methods brought by iterations. The EP-based COD is verified to yield identical solutions as the joint optimization and consumes less computation time in scenarios with numerous subsystems.

The remainder of this paper is organized as follows. Section II introduces the PVE algorithm and its properties. Sections III and IV apply the EP-based COD method to the multi-area coordinated dispatch problem and transmission-distribution coordination problem, respectively. Section V concludes this paper and discusses future research.

II Calculation Method for EP

II-A Problem Setup

In power system applications, the linear model is mostly utilized for optimal dispatch for the sake of computational efficiency, reliability, and clear economic interpretation. A broad spectrum of literature has investigated the linearized modeling of the transmission network [13, 14], distribution network [15, 16], and multi-energy system [17]. Nonlinear models, in contrast, have intractable computation and may lead to incentive issues when non-convexity exists [18]. The application of nonlinear optimization is difficult even for centralized dispatch, let alone coordinated dispatch. Part II of this paper focuses on the practice-oriented coordinated dispatch method and thus, the widely used linear dispatch model is considered.

As introduced in Part I, the EP eliminates internal variables from the operation feasible region of the subsystem and makes external equivalence of the subsystem model. The operation feasible region of the subsystem is enforced by both technical constraints and the epigraph of the objective function. Without loss of generality, represent the operation feasible region of the subsystem in the following compact form,

[TABLE]

In the model, $x$ and $y$ denote the coordination variable and internal variable, respectively. Constant $N_{x}$ and $N_{y}$ are dimensions of $x$ and $y$. The partition of variables is introduced in Part I of this paper. Matrix $A$, $B$, and $c$ are coefficients defining operation constraints of the subsystem. Note that variables are all bounded in power system dispatch problems due to operation limits of devices. Hence, $\Omega$ formulated in (1) is a polytope in the space of $\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{y}}$.

As per Definition 2 in Part I of this paper, the EP model of the subsystem is expressed as follows,

[TABLE]

The EP model depicts the technical and economic operation characteristics in the subspace of the coordination variable. The EP model can replace the original model of the subsystem in the coordinated optimization and can ensure the equivalent optimality in a non-iterative and privacy-protected manner, as proven in Part I.

II-B Representations of the EP model

As defined in (2), EP model $\Phi$ is the geometric projection of $\Omega$ onto the space of $\mathbb{R}^{N_{x}}$. Since $\Omega$ is a polytope, its projection $\Phi$ is also a polytope [19]. According to the Minkowski-Weyl Theorem [20], a polytope has the following two equivalent representations, known as the double descriptions.

• •

Hyperplane-representation (H-rep)

[TABLE]

• •

Vertex-representation (V-rep)

[TABLE]

In equation (3), $\tilde{A}$ and $\tilde{c}$ respectively denote the coefficient and the right-hand-side of the constraints that represent hyperplanes of the polytope. In equation (4), function $\text{conv}(\cdot)$ denote the convex hull of a set of vectors. Set $V$ contains all the vertices of $\Phi$ and $N_{v}$ is the number of vertices. $\hat{x}_{i}$ is the $i$th vertex of $\Phi$. The H-rep and V-rep of a polytope are convertible [20].

The EP calculation is to determine either the H-rep or the V-rep of $\Phi$. Though the two representations are equivalent, they will lead to different philosophies for designing the projection algorithm. In power system optimization problems, the H-rep is mostly used to model the operation feasible region. When calculating the projection, however, we find that the V-rep is more suitable due to some features of power systems. First, the dimension of the coordination variable is lower compared with internal variables of subsystems, since networks among subsystems are relatively weak-connected compared with networks inside subsystems. Hyperplane-oriented projection methods, e.g., the FME and block elimination, calculate the projection by eliminating internal variables and thus have exponential complexity regarding the dimension of the internal variable. By identifying vertices, in contrast, the projected polytope is directly characterized in the lower-dimensional coordination space, which can lower the computational difficulty. Second, from the perspective of practical implementation, the projection algorithm is desired to terminate within a given time and yield a conservative approximation of the EP model. If the hyperplane-oriented method is terminated prematurely, some critical hyperplanes enclosing the EP model will be missed and the calculation result may contain infeasible operating points. On the contrary, if the EP is calculated by identifying vertices, the intermediate output of the algorithm will always be a subset of the exact result due to the property of the convex set. In this regard, this paper develops a vertex-oriented projection method namely the PVE to calculate the EP efficiently and practically.

II-C PVE Algorithm

The PVE algorithm calculates the projected polytope by identifying its vertices. The key of the algorithm is twofold. First, the error metric has to be carefully designed to provide a clear interpretation of the calculation results and allow the pre-specification of the error tolerance. Second, the vertices have to be identified in a proper sequence to improve computational efficiency. To this end, the Hausdorff distance is employed to measure the approximation error and a double-loop framework is developed to identify vertices of the projection along a path with the steepest descent of the approximation error. The basic idea of the algorithm is first finding vertices that are critical to the overall shape of the projection, and then expanding the convex hull of existing vertices to find new vertices outside the current approximation. The flowchart is summarized in Algorithm 1. The PVE algorithm contains two layers of iterative loops; however, this algorithm is run locally by each subsystem to calculate the EP model and its loops do not conflict with that the proposed coordinated optimization method does not require iterative information exchange among subsystems. Details of the PVE algorithm will be introduced as follows along with a numerical example illustrated in Fig. 1.

II-C1 Vertex identification problem

Each vertex of the projected polytope $\Phi$ is an extreme point, which can be identified by solving the following problem,

[TABLE]

The value vector $\alpha$ in the objective function represents the direction for identifying the vertex.

II-C2 Initialization

At least $N_{x}+1$ initial vertices are required to ensure the convex hull of vertices is degenerate. These vertices can be searched along each axis, i.e., solving problem (5) with $\alpha=\pm e_{i}$, where $e_{i}$ is the $i$th standard basis whose $i$th component equals 1 while others equal 0. The set of initial vertices is denoted as $V^{(0)}$. In the illustrative example in Fig. 1 (a), the blue polygon is the exact projection of a randomly-generated polytope. Four initial vertices of the projection, i.e., $\hat{x}^{(0,1)},\hat{x}^{(0,2)},\hat{x}^{(0,3)}$, and $\hat{x}^{(0,4)}$, are identified for initialization.

II-C3 Inner loop

The inner loop of the PVE algorithm constructs convex hull of existing vertices and identifies new vertices outside the convex hull. The convex hull of existing vertices in the $k$th outer loop is denoted as $\hat{\Phi}^{(k)}$, which is an intermediate approximation of the projection. Assume $\hat{\Phi}^{(k)}$ is enclosed by $J^{(k)}$ facets and the affine hull of the $j$th facet is $\tilde{A}_{j}^{(k)}x\leq\tilde{d}_{j}^{(k)}$. New vertices are searched along outer normal directions of facets of $\hat{\Phi}^{(k)}$. For the $j$th facet, its outer normal vector is $\tilde{A}_{j}^{(k)}$. Thereby, the $j$th inner loop identifies the new vertex by solving problem (5) with $\alpha^{\top}=\tilde{A}^{(k)}_{j}$. The newly identified vertex is the optimal solution of problem (5), denoted as $\hat{x}^{(k,j)}$. The improvement ratio (IR) is calculated as (6) to measure the contribution of vertex $\hat{x}^{(k,j)}$ to improve the current approximation.

[TABLE]

The IR measures the distance from vertex $\hat{x}^{(k,j)}$ to the $j$th facet of $\hat{\Phi}^{(k)}$ and is non-negative. If $\Delta h^{(k,j)}>0$, indicating that $\hat{x}^{(k,j)}$ is outside $\hat{\Phi}^{(k)}$ and will contribute to improve the current approximation, then $\hat{x}^{(k,j)}$ will be appended to the vertex set $V$ and $\Delta h^{(k,j)}$ will be recorded in set $H^{(k)}$. If $\Delta h^{(k,j)}=0$, indicating that $\hat{x}^{(k,j)}$ is on the facet of $\hat{\Phi}^{(k)}$ and will not contribute to improve the current approximation, then this vertex will be omitted.

The example in Fig. 1 (a) exhibits the vertex identification in the first outer loop ($k=1$). The convex hull of existing vertices is $\hat{\Phi}^{(1)}$, as the red polygon shows. Vertices $\hat{x}^{(1,1)},\hat{x}^{(1,2)},\hat{x}^{(1,3)}$, and $\hat{x}^{(1,4)}$ are identified along the outer normal vector of the four facets of $\hat{\Phi}^{(1)}$, respectively. The Distance of each vertex to the corresponding facet of $\hat{\Phi}^{(1)}$ is the IR, as the dotted line mark in Fig. 1 (a).

II-C4 Outer loop and termination criterion

The outer loop evaluates the error of the current approximation and compares it with the pre-specified error tolerance to decide whether to terminate the algorithm. The Hausdorff distance between the real projection $\Phi$ and the current polytope $\hat{\Phi}^{(k)}$ is employed as the error metric, which is defined as follows [21],

[TABLE]

The interpretation of $D^{(k)}$ is the maximum distance from points in $\Phi$ to polytope $\hat{\Phi}^{(k)}$. Note that $\hat{\Phi}^{(k)}$ is the convex hull of existing vertices of $\Phi$, therefore $\hat{\Phi}^{(k)}\subset\Phi$. Hence, $D^{(k)}$ will be the maximum distance from vertices of $\Phi$ to polytope $\hat{\Phi}^{(k)}$. According to the definition in (6), the distance from each identified vertex of $\Phi$ to polytope $\hat{\Phi}^{(k)}$ is actually the corresponding IR. Thereby, the error metric can be evaluated as follows,

[TABLE]

From the above derivation, the Hausdorff distance can be evaluated by selecting the largest value from a set of scalars, which will consume much less computation effort compared to error metrics based on the region volume. If $D^{(k)}$ is no larger than the pre-specified error tolerance $\varepsilon$, the algorithm terminates and outputs the convex hull of existing vertices as the approximation of the real projection. The error metric based on the Hausdorff distance ensures that the maximum deviation of the approximation is within $\varepsilon$. If the termination criterion is not met, the algorithm moves into the next outer loop to identify more vertices outside the current approximation.

In Fig. 1 (a), the error metric $D^{(1)}$ for outer loop $k=1$ equals the IR of vertex $\hat{x}^{(1,3)}$, which has the farthest distance to $\hat{\Phi}^{(1)}$. Fig. 1 (b) exhibits the number of new vertices and the error metric in each outer loop of the PVE algorithm. The algorithm terminates at the $4$th outer loop when $D^{(4)}=0$. A total of 22 vertices of the projection are identified. As shown in the figure, the error metric keeps decreasing along the steepest path during the algorithm, indicating that vertices that have significant influences on the projection are identified with higher priority.

With the Hausdorff distance-based error metric, the process of the PVE algorithm can be interpreted as searching vertices that make the approximation error decline fast. This error metric also has another two advantages. First, the Hausdorff distance can be evaluated conveniently according to (8), which does not bring additional computational complexity. Second, the evaluation of the Hausdorff distance does not rely on the ground truth of the projection result and thus, the proposed error metric can be used to control the accuracy during the projection algorithm. In contrast, other error metrics, e.g., the volume difference metric and the Jaccard similarity, rely on the knowledge of the exact projection result and consequently, are hard to calculate and cannot be used to control the accuracy of the projection algorithm.

II-D Discussions

II-D1 Properties of the PVE algorithm

The convergence of the PVE algorithm and the conservatism of the approximation are guaranteed by the following theorem.

Theorem 1**.**

If error tolerance $\varepsilon>0$, the output of the PVE algorithm is a conservative approximation of the exact projection. If $\varepsilon=0$, the exact projection can be obtained via the PVE algorithm within finite calculations.

Proof.

In each round of the outer loop of the PVE algorithm, the convex hull of existing vertices is an intermediate approximation of the projection. Since the projected region is a polytope, the convex hull of its vertices is a subset of the exact projection, i.e., $\hat{\Phi}^{(k)}\subset\Phi$ for any $k$. Hence, if $\varepsilon>0$, a conservative approximation of the exact projection will be obtained with the maximum Hausdorff distance no larger than $\varepsilon$. If $\varepsilon$ is set 0, indicating that no point of $\Phi$ is outside $\hat{\Phi}^{(k)}$ when the algorithm terminates, then $\Phi\subset\hat{\Phi}^{(k)}$. Note that $\hat{\Phi}^{(k)}\subset\Phi$ and thus, $\hat{\Phi}^{(k)}=\Phi$. Since the number of vertices of a polytope is finite, the algorithm will terminate within finite steps. ∎

II-D2 Accelerating strategies

Three strategies are proposed to further accelerate the PVE algorithm. The first is to accelerate the solving of the vertex identification problem in (5). Note that only the value vector in the objective function is varying for identifying different vertices, the structure and parameters of constraints are invariant. Thereby, the redundant constraint elimination technique and the warm simplex basis technique can be used to accelerate the solution of (5). The redundant constraint elimination technique removes constraints that do not impact the feasible region to reduce the problem scale. Typical methods can be found in literature [22] and [23]. The warm simplex basis can be specified for problem (5) using the solution of the adjacent vertex that is on the same facet of which the outer normal vector is used as the searching direction. The specification of the start basis can reduce the iteration number of the simplex algorithm and thus accelerating the vertex identification. The second strategy is to accelerate the inner loop of the PVE algorithm by parallel computing. In line 7 of the algorithm, the identification of each vertex is independent, which can be calculated in parallel to save the computation time. The third is to dynamically update the convex hull when adding new vertices in each outer loop (line 4 of Algorithm 1), instead of building a new convex hull for all vertices. The updating of the convex hull can be realized by the quick hull method [24], which is not only computationally efficient but also capable of constructing high-dimensional convex hulls.

II-D3 Analysis of adaptability

First, for the problem where the projection region is a general convex set, the proposed PVE algorithm can also be applied. Perimeter points of the projection region will be identified by the algorithm, and the convex hull of these perimeter points is an inner approximation of the projection result. The approximation error measured by the Hausdorff distance is less than $\varepsilon$. Second, if the projection region is non-convex, the PVE algorithm cannot be directly applied. A potential solution is decomposing the non-convex projection region into the union of several convex sub-regions and using the PVE algorithm to identify each sub-region. Third, for the problem involving multiple time intervals, the coordination variable may have a high dimension and the direct application of the PVE algorithm will be inefficient. There are strategies to decompose the time-coupled projection problem into a series of lower-dimensional subproblems regarding each single time interval and adjacent intervals [25]. After decomposition, each low-dimensional projection can be calculated efficiently by the PVE algorithm.

III Application to Multi-Area Coordinated Optimal Dispatch

III-A Problem Formulation

Power systems in different areas are interconnected physically through inter-area transmission lines, but different areas may be operated by different system operators. The multi-area coordinated optimal dispatch (MACOD) is thus becoming indispensable for the economic and secure operation of large-scale power systems. The MACOD minimizes the total operation cost of the multi-area system in a decomposed manner.

Use $P$, $\pi$, and $\theta$ to represent variables of active power, operation cost, and voltage phase angle, respectively. Use superscript $G$, $D$, $B$, $RN$, $MA$, and $TL$ to label generation, load, power exchange at boundary node, regional network, multi-area system, and tie-line, respectively. Let $r\in\mathcal{R}^{RN}$ and $s\in\mathcal{S}$ index areas and tie-lines, respectively. Let $n\in\mathcal{N}^{RN}_{r}$ and $l\in\mathcal{L}^{RN}_{r}$ index network nodes and internal branches of the $r$th area, respectively. Then the objective function of a basic MACOD is as follows,

[TABLE]

where $\pi^{RN}$ is the operation cost of area $r$ and $\pi^{MA}$ is the total cost of the multi-area system.

The following constraints need to be satisfied for area $r$.

[TABLE]

In the above model, equation (10a) models the operation cost of area $r$ in the epigraph form, where $\pi^{RN,G}_{r,n}$ is the cost of generator $n$ and $\overline{\pi}^{RN}_{r}$ is a large enough constant to bound $\pi^{RN}_{r}$. Equation (10b) is the epigraph of the piece-wise linear cost function of generator $n$, where $a^{RN,G}_{r,n,i}$ and $b^{RN,G}_{r,n,i}$ are coefficients of each bidding segment, respectively. Equation (10c), (10d), and (10e) enforce the power balance, transmission limits, and generation limits of area $r$. In the equations, $T_{r,l,n}$ is the power transfer distribution factor of node $n$ to branch $l$, $\overline{F}^{RN}_{r,l}$ is the capacity of branch $l$, $\underline{P}^{RN,G}_{r,n}$ and $\overline{P}^{RN,G}_{r,n}$ are lower and upper bounds of generation output, respectively.

The operating constraints of tie-lines are as follows,

[TABLE]

Equation (11a) is the direct-current power flow model for tie-lines, where $x_{s}$ is the reactance of tie-line $s$, $r^{F}_{s}$ and $r^{T}_{s}$ denote the from and end region of the tie-line, respectively. In this work, each area is aggregated as a node when evaluating power flows on tie-lines, which is widely accepted in the literature [26] and is also implemented in the flow-based market integration in Europe [27]. Equation (11b) maps tie-line power flows to boundary power injections of each area, where $n^{F}_{s}$ and $n^{T}_{s}$ denote the from and end node of tie-line $s$, respectively. Equation (11c) enforces transmission limits of tie-lines.

III-B EP-Based Solution for MACOD

In the MACOD problem, variable $P^{RN,B}_{r,n}$ appears in constraints of both the regional system and tie-line system, which restricts the independent optimization of each area. Based on the EP theory proposed in Part I of this paper, the MACOD problem can be solved in a decomposed manner without iterative information exchange.

From the perspective of primal decomposition, $P^{RN,B}_{r}=(P^{RN,B}_{r,1},\cdots,P^{RN,B}_{r,|\mathcal{N}^{RN,B}_{r}|})$ and $\pi^{RN}_{r}$ are seen as the coordination variable corresponding to area $r$, while $P^{RN,G}_{r}=(P^{RN,G}_{r,1},\cdots,P^{RN,G}_{r,|\mathcal{N}^{RN,G}_{r}|})$ is the internal variable. According to (1), the operation feasible region of each area is as follows,

[TABLE]

According to (2), the EP model of the area is the projection of $\Omega^{RN}_{r}$ onto the subspace of coordination variable, i.e.,

[TABLE]

The EP model of the regional system contains all values of $(P^{RN,B}_{r},\pi^{RN}_{r})$ that can be met by at least one generation schedule subject to regional operation constraints with generation cost no larger than $\pi^{RN}_{r}$. Note that the objective of the MACOD is to minimize the sum of $\pi^{RN}_{r}$. Using the EP model as a substitute for the regional system model (10) in the MACOD, the cost variable $\pi^{RN}_{r}$ will be minimized and the decision result of $P^{RN,B}_{r}$ will be ensured to be feasible for the regional system, which leads to the coordinated optimality of multiple areas. Theoretical proof of the optimality is in Theorem 3 in Part I of this paper.

Then the EP-based MACOD is realized following three steps.

Equivalent projection. Each area calculates the EP model $\Phi^{RN}_{r}$ of the local system according to the definition in (13), which can be realized by the PVE algorithm introduced in Section II-C. Whereafter, $\Phi^{RN}_{r}$ is submitted to the multi-area coordinator. 2. 2)

Coordinated optimization. The coordinator solves the EP-based MACOD problem formulated as $\min\{\sum_{r}\pi^{RN}_{r}:\text{Eq. }\eqref{macd:cst_tie},(P^{RN,B}_{r},\pi^{RN}_{r})\in\Phi^{RN}_{r},\forall r\in\mathcal{R}^{RN}\}$. The optimal value $(\hat{P}^{RN,B}_{r},\hat{\pi}^{RN}_{r})$ is then published to each area. 3. 3)

Regional system operation. Each area fixes $\hat{P}^{RN,B}_{r}$ as the boundary condition and solves the local optimal dispatch problem formulated as $\min\{\pi^{RN}_{r}:\text{Eq.}\eqref{macd:cst_area},P^{RN,B}_{r}=\hat{P}^{RN,B}_{r}\}$.

The above process is executed sequentially, which naturally overcomes drawbacks caused by repetitive iterations of conventional coordinated optimization methods.

III-C Case Study

The EP of the regional transmission system is visualized, and the computational performance of the EP calculation and the EP-based MACOD is tested. Case studies in this paper are simulated on a Lenovo ThinkPad X13 laptop with an Intel Core i7-10510U 1.80-GHz CPU. Algorithms are programmed in MATLAB R2020a with GUROBI V9.0.0.

III-C1 EP of regional system

We employ the IEEE-24 system for visualizing the EP model of the regional transmission system. Two tie-lines are connected to node 1 and node 3 of the system, and each of them has the capacity of $510\text{\,}\mathrm{M}\mathrm{W}$ (15% of the total generation capacity). The results are shown in Fig 2, where subplots (a) and (b) exhibit EP models at the valley-load hour (77% of the peak load) and peak-load hour, respectively. Red regions are EP models of the system, which are 3-dimensional polytopes. As can be seen, the volume of the EP model at the peak hour is smaller than that at the valley hour. This is because the increase of load occupies the export capacity and upraising the power supply cost. The projection of the EP model on the subspace of $(P^{RN,B}_{r,1},P^{RN,B}_{r,2})$, denoted as $\Xi^{RN}_{r}$, is also plotted, which characterizes the admissible set of tie-lie flows that can be executed by the regional system.

III-C2 EP calculation

We test the computational performance of the proposed PVE algorithm for EP calculation based on the IEEE-24 test system, the 200 node-synthetic grid (SG-200), and the 500 node-synthetic grid (SG-500) [28]. Cases with 3 and 6 tie-lines are considered for each test system. The FME is employed as the benchmark method. As summarized in TABLE I, the FME fails to obtain the EP model within 1200s even for the small-scale IEEE-24 system. The proposed PVE algorithm, in contrast, can yield the EP model within 4s for all the 6 cases, which is acceptable for practical application. Measure the scale of the regional system model (10) by the production of numbers of variables and constraints. With the EP, the model scale of the regional system is reduced by 92.7%$\sim$99.9%, which will greatly alleviate the communication and computation burden of the coordinated dispatch of multiple areas as well as protecting the private information of regional systems.

III-C3 EP-based MACOD

Four test systems composed of 20 and 40 synthetic grids of different scales are constructed according to the tie-line topology in [29]. The MACOD is solved by the joint optimization and the EP-based decomposed solution, respectively. The optimization results from the two methods are verified to be identical, which validates the accuracy of the proposed coordination method. According to the process introduced in Section III-B, the time for system reduction depends on the region that takes the longest computation time, since EP models of different regions are calculated in parallel. The total time to obtain the optimal tie-line scheduling is the sum of the system reduction time and the coordinated optimization time. Results are summarized in TABLE II. As is shown, the total computation time of the EP-based MACOD is comparable with that of the joint optimization. With the EP, the time consumed by the coordinated optimization is reduced by more than 73.8% compared with the joint optimization. This is because the scale of the coordination problem is significantly reduced with the EP of each region. In addition to the computation efficiency, the primary advantages of the proposed coordination method are avoiding private data disclosure compared with the joint optimization and avoiding iterative information exchange compared with conventional coordinated optimization methods.

IV Application to Transmission-Distribution Coordinated Optimal Dispatch

IV-A Problem Formulation

Distribution networks installed with DERs can provide flexibility to the transmission system by adjusting the power exchange at the substation, which calls for the transmission-distribution coordinated optimal dispatch (TDCOD). In this study, optimal active power dispatch is considered for the transmission system, while the co-optimization of active and reactive power is considered for the distribution system with voltage limits. Let $P$, $Q$, $V$, and $\pi$ denote variables of active power, reactive power, voltage magnitude, and operation cost, respectively. Similar to the notation in Section III, use superscript $G$, $D$, $F$ to label generation, load, and power flow, respectively. Use superscript $TN$ and $DN$ to label variables of transmission and distribution systems, respectively. Let $r\in\mathcal{R}^{DN}$ index distribution networks participate in the coordinated dispatch. Let $n\in\mathcal{N}^{DN}_{r}$ and $l\in\mathcal{L}^{DN}_{r}$ index nodes and branches of distribution network $r$. Let $n\in\mathcal{N}^{TN}$ and $l\in\mathcal{L}^{TN}$ index nodes and branches of the transmission network.

The objective of the TDCOD is minimizing the total operation cost of the transmission and distribution systems,

[TABLE]

The power balance constraint, transmission limits, and generation limits are considered for the transmission network,

[TABLE]

The operation constraints of the $r$th distribution network is as follows,

[TABLE]

Equation (16a) and (16b) model operation costs of the distribution system and DER $n$ in the epigraph form, where $\overline{\pi}^{DN}_{r}$ is a large enough constant, $a^{DN,G}_{r,n,i}$ and $b^{DN,G}_{r,n,i}$ are coefficients of the $i$th segment of the piece-wise linear cost function of DER $n$. Equation (16c) enforces nodal power balance, where $P^{DN,F}_{r,mn}$ and $Q^{DN,F}_{r,mn}$ respectively denote active and reactive power flow on the feeder from node $m$ to node $n$, $\mathcal{A}_{n}$ is the set of nodes adjacent to node $n$. $P^{DN}_{r,0}$ and $Q^{DN}_{r,0}$ are net active and reactive power the distribution network injects into the transmission system. Equation (16d) is the simplified DistFlow model of the distribution network [30], which incorporates the reactive power and can be used to the distribution network with large R/X ratio. Equation (16e) enforces the transmission limit of branch $mn$ with piece-wise inner approximation and $\overline{F}^{DN}_{r,mn}$ is the capacity of the branch. Equation (16f) and (16g) enforce limits on nodal voltage magnitudes and power output of DERs, respectively. Taking $V^{2}_{r,n}$ as the independent variable, the above constraints are linear.

IV-B EP-based Solution for TDCOD

In problem (14)-(16), the optimization of transmission and distribution systems are coupled by $P^{DN}_{r,0}$ and $\pi^{DN}_{r}$. For the $r$th distribution network, take $x^{DN}_{r}=(P^{DN}_{r,0},\pi^{DN}_{r})$ as the coordination variable and take $y^{DN}_{r}=(Q^{DN}_{r,0},P^{DN,G}_{r},Q^{DN,G}_{r},P^{DN,F}_{r},Q^{DN,F}_{r},V^{2}_{r})$ as the internal variable. The operation feasible region of the distribution network, denoted as $\Omega^{DN}_{r}$, is the set of $(x^{DN}_{r},y^{DN}_{r})$ subject to constraints in (16). Then the EP model of the distribution network can be formulated according to the definition in (2), denoted as $\Phi^{DN}_{r}$. The EP model is the projection of $\Omega^{DN}_{r}$ onto the subspace of $x^{DN}_{r}$. Note that constraints (16) are linear and thus, both $\Omega^{DN}_{r}$ and its projection $\Phi^{DN}_{r}$ are polytopes, and $\Phi^{DN}_{r}$ can be calculated by the proposed PVE algorithm. The EP model $\Phi^{DN}_{r}$ contains all possible combinations of $P^{DN}_{r,0}$ and $\pi^{DN}_{r}$ that can be executed by the distribution network with at least one internal generation schedule satisfying constraints in (16). Using EP models of distribution networks to replace their original model in the TDCOD problem, the operation cost $\pi^{DN}_{r}$ of the distribution network will be minimized as in (14), and the active power exchange $P^{DN}_{r,0}$ determined by the transmission system operator will be ensured to be feasible for the distribution network. Hence, the EP-based coordinated optimization of transmission and distribution systems will be identical to the joint optimization.

The process of the EP-based TDCOD contains three successive steps.

Equivalent projection. Each distribution system calculates the EP model $\Phi^{DN}_{r}$ and reports it to the transmission system operator. 2. 2)

Coordinated optimization. The transmission system operator solves the EP-based TDCOD problem, i.e., minimizing the objective function in (14) subject to constraint (15) and $(P^{DN}_{r,0},\pi^{DN}_{r})\in\mathcal{R}^{DN}$. The optimal value $(\hat{P}^{DN}_{r,0},\hat{\pi}^{DN}_{r})$ is then published to each distribution system. 3. 3)

Distribution system operation. Each DSO fixes $\hat{P}^{DN}_{r,0}$ as the boundary condition and optimally dispatches the local system by solving the problem $\min\{\pi^{DN}_{r}:\text{Eq.}\eqref{tdcd:dn},P^{DN}_{r,0}=\hat{P}^{DN}_{r,0}\}$.

The above coordination process does not require iterations between the transmission and distribution systems. This coordination scheme is also compatible with the existing transmission system dispatch that schedules the optimal output of resources with their submitted information.

IV-C Case Study

IV-C1 EP of the distribution network

The modified IEEE-13 test feeder with 6 DERs in reference [31] is used to demonstrate the EP of the distribution network. Red polygons in Fig. 3 exhibit EP models of the test system at the valley and peak hours, which are 2-dimensional regions regarding the net power exchange $P^{DN}_{r,0}$ and the operation cost $\pi^{DN}_{r}$. The projection of the EP model on the subspace of $P^{DN}_{r,0}$ is denoted as $\Xi^{DN}_{r}$, which characterizes the range of flexibility that the distribution network can provide. The green curves with triangular markers are cumulative cost functions of DERs. As can be seen, the cumulative cost curve is beneath the EP model and has a wider range of net power exchange. This is because network limits are omitted when generating the cumulative cost curve. Hence, only aggregating cost curves of DERs for transmission-level dispatch is not enough, the EP of the distribution system incorporating both cost functions and network constraints is needed for the optimal coordination between transmission and distribution systems.

IV-C2 EP calculation

Multiple IEEE-13 feeders are connected at the root node to create larger-scale test systems with 25, 49, 241, and 2401 nodes to examine the computation performance of the EP calculation. As shown in TABLE III, the proposed PVE algorithm successfully obtains EP results of all five test systems within 3 seconds. In contrast, the conventional FME method takes more than 10 and 129 times of computational effort to calculate the EP for the 13-node and 25-node test systems, respectively. For test systems with more than 49 nodes, the FME fails to yield the EP results within 1200s. This result shows the superiority of the PVE-based EP calculation in terms of computational efficiency and scalability. As for the reduction of the model scale, the EP model of the distribution network is described by a group of 2-dimensional linear constraints, which reduces the model scale of the distribution network by more than 98.4%.

IV-C3 EP-based TDCOD

The IEEE-24 system is connected with different numbers of IEEE-13 distribution feeders to simulate the coordinated dispatch between transmission and distribution systems. The coordinated dispatch is solved by the joint optimization and the EP-based coordination method respectively. Optimization results from the two methods are verified to be identical for all the test cases. The computation time of the joint optimization and the EP-based coordination is summarized in TABLE IV. The computation time of the EP-based coordination is dominated by the system reduction process. With the EP, the time of the coordinated optimization is reduced by more than 69.2%, which may help to relieve the computation burden of the transmission system. The total computation time of the EP-based coordination is also reduced compared with the joint optimization in cases with more than 100 distribution networks. This is because the EP distributes the computation burden of solving a large-scale centralized optimization to each subsystem in parallel and thus, is more efficient for the coordinated management of numerous distribution networks.

V Conclusion

This two-part paper proposes the EP theory to make external equivalence of the subsystem and realizes the COD of power systems in a non-iterative fashion. To calculate the EP efficiently and accurately, Part II of this paper proposes a novel polyhedral projection algorithm termed the PVE, which characterizes the EP model by identifying its vertices and building the convex hull of vertices. The vertex identification process in the PVE algorithm is finely designed to give higher priority to vertices that are critical to the overall shape of the projection, which makes the approximation error decreases along the steepest path. The Hausdorff distance is employed to measure and control the calculation accuracy of the PVE algorithm, which can provide flexibility to balance the computation accuracy and computation effort of the EP calculation in practice. The EP theory and the PVE algorithm are applied to the non-iterative COD of multi-area systems and transmission-distribution systems. Case studies of different scales testify the superiority of the proposed PVE algorithm in terms of computational efficiency and scalability compared with the conventional FME algorithm. The EP is verified to reduce the subsystem model scale by more than 92.7% for regional transmission systems and more than 98.4% for distribution networks, which can alleviate the computation and communication burden in coordinated dispatch. The effectiveness of the EP-based non-iterative coordinated dispatch is also validated in comparison with the joint optimization.

The proposed EP theory and the non-iterative coordinated optimization framework are general and may find a broad spectrum of applications in addition to instances in this paper, e.g., the coordinated dispatch of the multi-energy network, the coupling of power and traffic networks, and the integration of user-side flexible resources. Applications of the PVE algorithm to other projection problems, e.g., flexibility region aggregation, loadability set [2], and the projection-based robust optimization [32], are also worthy of future investigation. Incorporating the robust feasibility and chance constraints to deal with the uncertainty and incorporating inter-temporal constraints in the projection calculation should be future extensions of the proposed method.