Efficient Operations of Micro-Grids with Meshed Topology and Under Uncertainty through Exact Satisfaction of AC-PF, Droop Control and Tap-Changer Constraints

Abstract

Micro-grids’ operations offer local reliability; in the event of faults or low voltage/frequency events on the utility side, micro-grids can disconnect from the main grid and operate autonomously while providing a continued supply of power to local customers. With the ever-increasing penetration of renewable generation, however, operations of micro-grids become increasingly complicated because of the associated fluctuations of voltages. As a result, transformer taps are adjusted frequently, thereby leading to fast degradation of expensive tap-changer transformers. In the islanding mode, the difficulties also come from the drop in voltage and frequency upon disconnecting from the main grid. To appropriately model the above, non-linear AC power flow constraints are necessary. Computationally, the discrete nature of tap-changer operations and the stochasticity caused by renewables add two layers of difficulty on top of a complicated AC-OPF problem. To resolve the above computational difficulties, the main principles of the recently developed “$l_1$-proximal” Surrogate Lagrangian Relaxation are extended. Testing results based on the nine-bus system demonstrate the efficiency of the method to obtain the exact feasible solutions for micro-grid operations, thereby avoiding approximations inherent to existing methods; in particular, fast convergence of the method to feasible solutions is demonstrated. It is also demonstrated that through the optimization, the number of tap changes is drastically reduced, and the method is capable of efficiently handling networks with meshed topologies.

1. Introduction

As the name suggests, a micro-grid is a relatively small geographically localized electricity grid intended to provide uninterrupted service to local communities such as campuses, business centers, hospital complexes, and other critical infrastructure. To enable self-sufficiency to hedge against blackouts caused by natural disasters, micro-grids typically include generators, combined heat and power, and batteries. With the recent push for clean, green, and renewable energy, micro-grids may also include solar panels. Other grid devices may also include electric vehicle (EV) charging stations.

Under normal conditions, micro-grids are typically connected to the main grid (e.g., a power distribution system) and may exchange power. The normal operations, therefore, include proactive power generation at the least cost in anticipation of the increase/decrease in customer demand, as well as in anticipation of fluctuations of power generation from renewables. The intermittency of renewables can, however, lead to fluctuations of voltage, which may lead to frequent adjustments of taps within on-load tap-changers, thereby leading to fast degradation of expensive equipment and adversely impacting micro-grid economic viability. Stochasticity, as well as the discrete nature of the underlying problem, needs to be explicitly captured to design the optimal (or near-optimal) control to ensure a low-cost power supply to local communities.

Under faulty conditions or upon detection of low voltage/frequency on the main grid side, an interconnecting device such as a circuit breaker opens to switch the micro-grid to the islanded mode. In this mode, the micro-grid’s goal is to serve the local loads by using locally available distributed energy resources (DERs). The micro-grid’s operations are also complicated by the fact that disconnection may result in a drop in voltage and frequency. Therefore, in addition to the considerations described in the previous paragraph, droop controls need to be used to restore voltage and frequency to nominal ranges.

From the modeling standpoint, the main goal of the paper is to formulate a micro-grid optimization problem while including AC power flows to explicitly capture voltage fluctuations. To control voltage fluctuations, tap-changer, as well as droop-control constraints, will be included in the formulation. Uncertainties will be captured through the use of Markov processes. The resulting problem, while generally smaller than the optimization problems solved at the main grid level, is complicated due to non-convexities caused by non-linear AC power flows and the discrete nature of tap changes (tap positions). Rather than using AC power flow approximations (such as DistFlow or Second-Order Cone Relaxation), an exact AC power flow in rectangular coordinates will be used.

From a methodological standpoint, the main goal of the paper is to develop an algorithm to coordinate islanded micro-grid resources at a high level while overcoming the above-mentioned difficulties, rather than to provide low-level detailed modeling of micro-grid devices. Therefore, in Section 2, important micro-grid components/features such as tap-changers, droop controls, AC power flow, as well as a Markovian approach to modeling uncertainties will be reviewed. Modeling of batteries and EV charging stations is out of the scope of this paper, although the method developed in Section 3 is intended to support plug-and-play capabilities; therefore, the consideration of storage and other DERs will only affect the way subproblems are solved at the low level and will not affect the overall coordination at the high level. The gist of the “$l_1$-proximal” Surrogate Lagrangian Relaxation decomposition and coordination methodology is the relaxation of coupling constraints (e.g., nodal flow balance) to reduce the complexity of the problem while coordinating nodal subproblems. Non-linear constraints are linearized to both further reduce complexity and enable the use of mixed-integer linear programming (MILP) solvers. To ensure the overall feasibility, “$l_1$-proximal” terms, also amenable to the use of MILP solvers, are introduced.

In Section 4, by considering a nine-bus micro-grid, the ability of the method to efficiently handle non-linearity, non-convexity, and stochasticity of the underlying problem while guaranteeing feasibility is demonstrated. It is also demonstrated that the number of tap changes is drastically reduced, and that the method is capable of handling networks with meshed topologies while overcoming non-linearity difficulties caused by AC power flow, tap changer, and droop-control constraints.

2. Literature Review

2.1. AC Power Flow

Because of the high $r/x$-ratios within distribution systems, the popular DC power flow model [1] is no longer suitable for formulating the power flows within micro-grids as well. In the following, AC Power Flow modeling methodologies will be reviewed.

1.

AC Power Flow for Radial Topologies.

The AC Power Flows are known for their non-linearity and non-convexity and the associated AC optimal power flow (AC-OPF) problem is known to be extremely difficult. One of the popular ways to handle the non-linearities is the so-called DistFlow model [2,3,4,5,6]. The model has been popular within distribution and micro-grids due to the general ease of handling linear constraints and fairly high accuracy of the approximation. The exactness of solutions, however, cannot be guaranteed. To guarantee the exactness, the model has to be non-linear. The so-called Second-Order Cone Relaxation (SOCR) [6,7] addressed the exactness issue and the resulting model is convex, which is amenable for the use of commercial solvers unlike the original non-convex AC Power Flow.

2.

AC Power Flow for Meshed Topologies.

The distribution networks (and micro-grids as their part) are not necessarily radial, as acknowledged for the next generation of power distribution systems [8]. The SOCR is, however, not exact to meshed topologies. The AC power flow in rectangular coordinates is appropriate for both radial and meshed topologies. One way to handle the associated non-convexity is by defining a convex hull of the solutions; however, even if the tightest convex relaxation is attained [9], the solution may still be inside of the convex hull rather than at its boundary, which could only be guaranteed for linear programming problems. Another approach to handling non-linearity and non-convexity is through dynamic linearization [10]—The linearization around the current operation point, which may change from iteration to iteration. To guarantee feasibility, the so-called “$l_1$-proximal” terms have been used to gradually penalize the violations of current solutions from previously obtained ones until convergence to a steady-state solution.

2.2. Intermittent Renewables

In several studies on the operation of micro-grids, deterministic approaches were adopted, whereby mean values have been used for renewable generation, which is inherently intermittent and uncertain. For example in [11], solar irradiation of Photovoltaic (PV) panels was calculated offline using a deterministic approach. Without explicitly capturing uncertainties, the solutions obtained using deterministic approaches are not robust against realizations of renewable generation. To explore the intermittent and uncertain nature, stochastic programming has also been used based on representative scenarios (e.g., [12]). It is, however, difficult to select an appropriate scenario number while balancing modeling accuracy, computational efficiency, and solution feasibility. To overcome these difficulties, a Markovian approach was developed [13]. Without considering power-transmission constraints, Markov chains were used to model renewable generation (e.g., wind) through aggregation; within the method, renewable states at a particular time period probabilistically capture the past information, thereby reducing the complexity compared to a total enumeration of all possible states; in particular, the complexity is significantly reduced compared to scenario-based methods. In our previous work [14], a Markov-based model was established to integrate intermittent and uncertain PV generation into micro-grids.

2.3. Tap-Changers

Tap changers [15] are used to reduce voltage deviations due to intermittent renewables by keeping voltage amplitudes within pre-specified limits. However, high levels of renewable penetration, voltage fluctuations, and amplitudes may increase greatly, thereby leading to frequent adjustment of the transformer tap position to preserve power quality and reliability. As a result, the tap-changer transformers may reach the end of their useful lives or experience premature failures. Hawaii utilities, for example, reported that their tap changers, previously maintenance-free for over 40 years, require quarterly maintenance with a prospect of retirement within a mere two years of exploitation, since the taps frequently need to be adjusted over 300 times a day due to voltage fluctuations [16].

3. Micro-Grid Model

Consider a micro-grid with a partly connected mesh topology operated by a micro-grid system operator (MSO) planning ahead T units of time (e.g., a unit of time could be 5, 15 or 60 min). Let buses be denoted by b, generators by i, and power lines by l. Solar generation is modeled by using N discrete states (denoted by n) with associated probabilities ${\varphi }_{n,t}$ at each time period t and transition probabilities ${\pi }_{n,m}$ from state n to state m.

3.1. Objective

The goal of MSO is to minimize the expected generation cost together with the expected tap-changer cost:

$$\begin{array}{cc}& \underset{\begin{array}{c}\mathbf{F},\mathbf{G},\mathbf{V},\mathbf{X}\end{array}}{min}O\left(\mathbf{G}\right)=\underset{\begin{array}{c}\mathbf{F},\mathbf{G},\mathbf{V},\mathbf{X}\end{array}}{min}\left\{\mathbb{E}(\mathbf{C}·{\mathbf{G}}^{{}^{\prime }})+\mathbb{E}({\mathbf{C}}^{tap}·{\mathbf{D}}^{{}^{\prime }})\right\},\hfill \end{array}$$

(1)

where $\mathbf{G}=\left\{g_{i,n,t}^p,g_{i,n,t}^q\right\}$ is a vector consisting of generation levels (both active (denoted by p) and reactive (denoted by q)) with the corresponding costs $\mathbf{C}=\left\{C_{i,t}^p,C_{i,t}^q\right\}$, and $\mathbf{D}=\left\{d_{b,n,t}^{up},d_{b,n,t}^{down}\right\}$ is a vector of the change of tap position up and down at PV state n with the corresponding tap-changing costs ${\mathbf{C}}^{tap}$. The expectation operator is defined as $\mathbb{E}(\mathbf{C}·{\mathbf{G}}^{{}^{\prime }})\equiv {\sum }_{i,n,t}{\varphi }_{n,t}·C_{i,t}^p·g_{i,n,t}^p+{\sum }_{i,n,t}{\varphi }_{n,t}·C_{i,t}^q·g_{i,n,t}^q$. The second term in (1) is similarly defined.

Other decision variables involved in (1) include $\mathbf{X}=\left\{x_{i,t}\right\}$—a vector of binary commitment decision variables, $\mathbf{F}=\left\{f_{l,n,t}^p,f_{l,n,t}^q\right\}$—a vector of continuous power flow decision variables and $\mathbf{V}=\left\{v_{b,n,t}^{Re},v_{b,n,t}^{Im}\right\}$—a vector of continuous voltage decision variables. The optimization (1) is subject to the following constraints.

3.2. Intra-Nodal (Local) Constraints

Generation Capacity Constraints. Generation levels satisfy the following generation capacity constraints:

$$\begin{array}{cc}& \underline{\mathbf{G}}·\mathbf{X}\le \mathbf{G}\le \overline{\mathbf{G}}·\mathbf{X},\hfill \end{array}$$

(2)

where $\underline{\mathbf{G}}=\left\{{\underline{g}}_{i,n,t}^p,{\underline{g}}_{i,n,t}^q\right\}$ and $\overline{\mathbf{G}}=\left\{{\overline{g}}_{i,t}^p,{\overline{g}}_{i,t}^q\right\}$ are the minimum and maximum generation levels, respectively.

Ramp-Rate Constraints. For probable transitions, ramp-rate constraints require that the change of generation levels between two consecutive time periods does not exceed ramp rates $\mathbf{R}=\left\{r_i^p,r_i^q\right\}$:

$$\begin{array}{cc}& -\mathbf{R}\le {\mathbf{G}}_{n,t}-{\mathbf{G}}_{m,t-1}\le \mathbf{R},\forall {\pi }_{n,m}\ne 0.\hfill \end{array}$$

(3)

Voltage Restrictions. The voltages are restricted as follows:

$$\begin{array}{cc}& {\underline{\mathbf{V}}}^2\le \mathbf{V}·{\mathbf{V}}^{{}^{\prime }}\le {\overline{\mathbf{V}}}^2,\hfill \end{array}$$

(4)

where ${\underline{\mathbf{V}}}^2=\left\{{\underline{v_b}}^2\right\}$ and ${\overline{\mathbf{V}}}^2=\left\{{\overline{v_b}}^2\right\}$ with $\underline{v_b}$ and $\overline{v_b}$ being minimum and maximum voltage magnitudes.

Droop-Control Constraints. It is assumed that generator $\widehat{i}\in {\mathcal{I}}_{\widehat{b}}$ employs a droop control strategy and the corresponding droop-control constraints are:

$$f_{n,t}=f_{n,t}^{ref}-k_f·(g_{\widehat{i},n,t}^{p,ref}-g_{\widehat{i},n,t}^p),$$

(5)

$$\begin{array}{cc}& \sqrt{\left(v_{\widehat{b},n,t}\right)·{\left(v_{\widehat{b},n,t}\right)}^{{}^{\prime }}}=\sqrt{\left(v_{\widehat{b},n,t}^{ref}\right)·{\left(v_{\widehat{b},n,t}^{ref}\right)}^{{}^{\prime }}}-k_v·(g_{\widehat{i},n,t}^{q,ref}-g_{\widehat{i},n,t}^q).\hfill \end{array}$$

(6)

PV generation. Following our previous work [14], a Markov-based model is adopted for PV generation. Uncertainties due to whether changes follow a Markovian process with N discrete states, with state n denoted as $W_n$. Based on historical data, the transition probability to state n from state m is ${\pi }_{m,n}$ as shown below in Figure 1.

The uncertain PV Generation $g_{n,t}^{PV}$ is also a Markovian process, as follows:

$$\begin{array}{c}{g}_{n,t}^{PV}=g_t^{IPV}{W}_n,\end{array}$$

(7)

where $g_t^{IPV}$ is the ideal PV generation at time n. The probability ${\varphi }_{n,t}$ that the PV generation is $g_{n,t}^{PV}$ at time t is the weighted summation of probabilities at time $t-1$:

$$\begin{array}{c}{\varphi }_{n,t}=\sum _{m=1..N}{\pi }_{n,m}{\varphi }_{n,t-1}.\end{array}$$

(8)

Generally, probabilities of PV levels at subsequent time-periods are obtained based on PV levels at previous time-periods and the transition matrix consisting of $\left\{{\pi }_{m,n}\right\}$.

Tap changer. Assume that tap changers are associated with the corresponding solar farms. Following [15], and assuming that bus indices are omitted for brevity, tap-changer constraints are:

$$\begin{array}{c}{V}_{n,t}=\frac{1}{a_{n,t}}{V}^{in}-\frac{a_{n,t}{Z}_{t\left(a\right)}{S}_n}{{\left(V^{in}\right)}^{*}},\end{array}$$

(9)

where $V^{in}$ is voltage input, a is a transformer-turn ratio under no-load conditions as: $V^{in}$/$V_n$, $Z_{t\left(a\right)}$ is transformer leakage impedance, generally a function of a, but here $Z_{t\left(a\right)}$ is assumed to be a complex number, and S is the transformer load. For mathematical expediency, (9) is equivalently written as:

$$\begin{array}{c}{a}_{n,t}{\left(V^{in}\right)}^{*}{V}_n=V^{in}{\left(V^{in}\right)}^{*}+{\left(a_{n,t}\right)}^2{Z}_{t\left(a\right)}{S}_n.\end{array}$$

(10)

The transformer turn ratios $a_{n,t}$, which are decision variables, are controlled as:

$$\begin{array}{c}{a}_{n,t}=a_0+d_{n,t}·\mathsf{\Delta }a,\end{array}$$

(11)

where $a_0$ is the “rated turn ratio”, $\mathsf{\Delta }a$ is the single tap position change, and $\left\{d_{n,t}\right\}$ are integer decision variables denoting tap positions as:

$$\begin{array}{c}{d}_{n,t}=d_{n,t-1}+d_{n,t}^{up}-d_{n,t}^{down}.\end{array}$$

(12)

3.3. Inter-Nodal (Global) Constraints

AC Power Flow in Rectangular Coordinates. Following [17,18,19,20,21], AC power flows are modeled by using complex voltages $v_{b,n,t}=\left(v_{b,n,t}^{Re},v_{b,n,t}^{Im}\right)$ as:

$$\begin{array}{cc}& f_{l,n,t}^p=v_{s\left(l\right),n,t}·\left(\begin{array}{cc}{g}_{s\left(l\right),r\left(l\right)}& -b_{s\left(l\right),r\left(l\right)}\\ b_{s\left(l\right),r\left(l\right)}& g_{s\left(l\right),r\left(l\right)}\end{array}\right)·{\left(v_{r\left(l\right),n,t}\right)}^{{}^{\prime }},\hfill \end{array}$$

(13)

$$\begin{array}{cc}& f_{l,n,t}^q=v_{s\left(l\right),n,t}·\left(\begin{array}{cc}-b_{s\left(l\right),r\left(l\right)}& -g_{s\left(l\right),r\left(l\right)}\\ g_{s\left(l\right),r\left(l\right)}& -b_{s\left(l\right),r\left(l\right)}\end{array}\right)·{\left(v_{r\left(l\right),n,t}\right)}^{{}^{\prime }}.\hfill \end{array}$$

(14)

Here, $b_{s\left(l\right),r\left(l\right)}$ is susceptance and $g_{s\left(l\right),r\left(l\right)}$ is conductance of line $\left(s\left(l\right),r\left(l\right)\right)$. Node $s\left(l\right)$ stands for a “sending” node, and $r\left(l\right)$ stands for a “receiving” node of line l.

Nodal Power Flow Balance Constraints. For every PV state n and at every node b, the net “inflow” should be equal to the net “outflow”:

$$\begin{array}{cc}& \sum _{i\in {\mathcal{I}}_b}{g}_{i,n,t}^p+g_{b,n,t}^{PV,p}+\sum _{\begin{array}{c}l=1:\\ r\left(l\right)=b\end{array}}^L{f}_{l,n,t}^p=L_{b,t}^p+\sum _{\begin{array}{c}l=1:\\ s\left(l\right)=b\end{array}}^L{f}_{l,n,t}^p,\hfill \end{array}$$

(15)

$$\begin{array}{cc}& \sum _{i\in {\mathcal{I}}_b}{g}_{i,n,t}^q+g_{b,n,t}^{PV,q}+\sum _{\begin{array}{c}l=1:\\ r\left(l\right)=b\end{array}}^L{f}_{l,n,t}^q=L_{b,t}^q+\sum _{\begin{array}{c}l=1:\\ s\left(l\right)=b\end{array}}^L{f}_{l,n,t}^q.\hfill \end{array}$$

(16)

If bus b contains no generators, then ${\sum }_{i\in {\mathcal{I}}_b}{g}_{i,n,t}^{p/q}=0$ (here, ${\mathcal{I}}_b$ is a set of generators at bus b); if bus b contains no load, then $L_{b,t}^{p/q}$ = 0, and if there is no PV generation, then $g_{b,n,t}^{PV,p/q}=0$.

Line-Capacity Constraints. For every PV state n and every line l, power flows satisfy the following line-capacity constraints:

$$\begin{array}{cc}& \sqrt{{\left(f_{l,n,t}^p\right)}^2+{\left(f_{l,n,t}^q\right)}^2}\le {\overline{f}}_l.\hfill \end{array}$$

(17)

The above problem belongs to a class of Mixed-Integer Non-linear Programming problems notable for non-convexities brought by non-linear AC power flows as well as by discrete decision variables.

4. Solution Methodology

To solve the above problem, the methodology “$l_1$-proximal” Surrogate Lagrangian Relaxation method [10] is extended to handle tap-changer and droop-control constraints. After linearizing the problem and relaxing nodal flow balance constraints, Lagrangian multipliers are updated based on constraint violations, and the violations are penalized using “absolute-value” $l_1$ penalties; to ensure overall feasibility “$l_1$-proximal” terms are used.

4.1. Surrogate Absolute-Value Lagrangian Relaxation

Relaxed Problem. After relaxing nodal flow balance (15)–(16), and penalizing their violations, the relaxed problem becomes:

$$\begin{array}{cc}& \underset{\begin{array}{c}\mathbf{F},\mathbf{G},\mathbf{V},\mathbf{X}\end{array}}{min}{L}_c(\mathbf{G};\Lambda )=\underset{\begin{array}{c}\begin{array}{c}\mathbf{F},\mathbf{G},\mathbf{V},\mathbf{X}\end{array}\end{array}}{min}\left\{O\left(\mathbf{G}\right)+\Lambda ·\mathbf{R}+c·\parallel \mathbf{R}{\parallel }_1\right\},\mathrm{s}.\mathrm{t}.,\left(2\right)–\left(14\right),(17),\hfill \end{array}$$

where ${\Lambda }_j=\left({\Lambda }_j^P,{\Lambda }_j^Q\right)$ are multipliers relaxing power flow balance constraints (16). The vector, $\mathbf{R}={({\mathbf{B}}^P,{\mathbf{B}}^Q)}^{{}^{\prime }}$ denotes a vector of power flow balance constraint violations.

Following [10], multipliers are updated as:

$$\begin{array}{c}{\Lambda }^k={\Lambda }^{k-1}+s^k·{\mathbf{R}}^k,\end{array}$$

(18)

where

$$\begin{array}{cc}& s^k={\alpha }^k·s^{k-1}·\frac{\parallel {\mathbf{R}}^{k-1}{\parallel }_2}{\parallel {\mathbf{R}}^k{\parallel }_2},\hfill \end{array}$$

(19)

where ${\alpha }^k$ is defined as

$$\begin{array}{cc}& {\alpha }^k=1-\frac{1}{M·k^{1-\frac{1}{k^r}}},M>1,r>0.\hfill \end{array}$$

(20)

In the beginning, $c^k$ is increased by a constant $\beta >1$:

$$\begin{array}{cc}& c^k=c^{k-1}·\beta .\hfill \end{array}$$

(21)

After constraint violations become zero, a feasible solution is obtained, and $c^k$ needs to decrease as:

$$\begin{array}{cc}& c^k=c^{k-1}·{\beta }^{-1}\hfill \end{array}$$

(22)

Subsequently, the $c^k$ is not increased.

4.2. Practical Implementation

Voltage restrictions (4), droop control for voltage (6), tap-changer constraints (10), AC power flows (13) and (14), as well as line-capacity constraints (17) are non-linear. Moreover, the left-hand side constraint of (4) and (13) and (14) delineate non-convex regions. Following the work [10], these difficulties are briefly addressed next. The method of [10] is then extended to improve the linearization of (4) as well as to resolve non-linearity difficulties brought by the newly considered droop control (6) and tap-changer constraints (10).

In the following, the main linearization principles will be delineated first. Then, the feasibility will be established.

1.

Linearization of Cross-Product Terms within (13) and (14). To linearize AC Power Flow while updating all the voltages, the following formula is used:

$$\begin{array}{cc}& {\widehat{f}}_{l,t}^p=\frac{1}{2}{v}_{s\left(l\right),t}^{k-1}·\left(\begin{array}{cc}{g}_{s\left(l\right),r\left(l\right)}& -b_{s\left(l\right),r\left(l\right)}\\ b_{s\left(l\right),r\left(l\right)}& g_{s\left(l\right),r\left(l\right)}\end{array}\right)·\left(v_{r\left(l\right),t}\right)+\frac{1}{2}{v}_{s\left(l\right),t}·\left(\begin{array}{cc}{g}_{s\left(l\right),r\left(l\right)}& -b_{s\left(l\right),r\left(l\right)}\\ b_{s\left(l\right),r\left(l\right)}& g_{s\left(l\right),r\left(l\right)}\end{array}\right)·\left(v_{r\left(l\right),t}^{k-1}\right).\hfill \end{array}$$

(23)

Reactive power is linearized in the same fashion.

2.

Linearization of Voltage Restrictions (4). First, the squared terms within (4) are linearized as:

$$\begin{array}{cc}& {\underline{\mathbf{V}}}^2\le {\mathbf{V}}^{k-1}·{\mathbf{V}}^{{}^{\prime }}\le {\overline{\mathbf{V}}}^2.\hfill \end{array}$$

(24)

To avoid infeasibility along the iterative procedure, “soft" penalization is introduced through the non-negative penalty variables ${\underline{\mathbf{V}}}^{pen}=\left\{{\underline{v}}_{b,t}^{pen}\right\}$ and ${\overline{\mathbf{V}}}^{pen}=\left\{{\overline{v}}_{b,t}^{pen}\right\}$ as:

$$\begin{array}{cc}& {\underline{\mathbf{V}}}^2-{\underline{\mathbf{V}}}^{pen}\le {\mathbf{V}}^{k-1}·{\mathbf{V}}^{{}^{\prime }}\le {\overline{\mathbf{V}}}^2+{\overline{\mathbf{V}}}^{pen}.\hfill \end{array}$$

(25)

To enforce the feasibility of (25), ${\underline{\mathbf{V}}}^{pen}$ and ${\overline{\mathbf{V}}}^{pen}$ are penalized.

3.

Linearization of Tap-Changer Constraints (10). Tap-changer constraints contain both the cross-products and the squared terms. Following the ideas of above-presented linearization, the linearization is performed as follows:

$$\begin{array}{c}\frac{1}{2}{a}_{n,t}^{k-1}{\left(V^{in}\right)}^{\ast }{V}_n+\frac{1}{2}{a}_{n,t}{\left(V^{in}\right)}^{\ast }{V}_n^{k-1}=V^{in}{\left(V^{in}\right)}^{\ast }+a_{n,t}^{k-1}{a}_{n,t}{Z}_{t\left(a\right)}{S}_n.\end{array}$$

(26)

The linearization of line-capacity constraints (17) is operationalized in the same way as described in point 2 and 3 above.

4.

Linearization Voltage Droop-Control Constraints (6). The first step of linearization of (6) is squaring both sides of the equation as such:

$$\begin{array}{cc}& {\left|v_{\widehat{b},n,t}\right|}^2={\left|v_{\widehat{b},n,t}^{ref}\right|}^2-2·\left|v_{\widehat{b},n,t}^{ref}\right|·\left(k_v·(g_{\widehat{i},n,t}^{q,ref}-g_{\widehat{i},n,t}^q)\right)+k_v^2·{\left(g_{\widehat{i},n,t}^{q,ref}-g_{\widehat{i},n,t}^q\right)}^2.\hfill \end{array}$$

(27)

The left-hand side and the last term of the right-hand side of (27) are squared terms, which are linearized as described above. Other terms of (27) are either constant or linear and require no further linearization.

To ensure that solutions to the linearized problem satisfy the corresponding original (non-linear) problem terms $\parallel \mathbf{V}-{\mathbf{V}}^{k-1}{\parallel }_1$, $\parallel \widehat{\mathbf{F}}-{\widehat{\mathbf{F}}}^{k-1}{\parallel }_1$, and $\parallel \mathbf{A}-{\mathbf{A}}^{k-1}{\parallel }_1$, which capture the deviations of voltages $\mathbf{V}$, linearized power flows $\widehat{\mathbf{F}}$ and transformer turns ratios $\mathbf{A}\equiv \left\{a_{n,t}\right\}$ from previously obtained values are introduced and are penalized by $c_p^k$. The intention here is to discourage oscillations of solutions while encouraging their approach to common values. To avoid getting trapped at previously obtained values, $c_p^k$ is chosen to be lower than $c^k$.

Linearized Relaxed Problem.

The resulting MSO relaxed problem then becomes:

$$\begin{array}{cc}& \underset{\begin{array}{c}\mathbf{A},\widehat{\mathbf{F}},\mathbf{G},\mathbf{V},\mathbf{X}\end{array}}{min}\left\{\begin{array}{c}{L}_{c^{k-1}}(\mathbf{G};{\Lambda }^{k-1})+c_p^{k-1}\left(\parallel \mathbf{V}-{\mathbf{V}}^{k-1}{\parallel }_1+\parallel \widehat{\mathbf{F}}-{\widehat{\mathbf{F}}}^{k-1}{\parallel }_1+\parallel \mathbf{A}-{\mathbf{A}}^{k-1}{\parallel }_1\right)+\\ {\underline{\Lambda }}^{pen}·\left({\underline{\mathbf{V}}}^2-{\mathbf{V}}^{k-1}·{\mathbf{V}}^{{}^{\prime }}\right)+{\overline{\Lambda }}^{pen}·\left({\mathbf{V}}^{k-1}·{\mathbf{V}}^{{}^{\prime }}-{\overline{\mathbf{V}}}^2\right)+c^{k-1}·{\underline{\mathbf{V}}}^{pen}+c^{k-1}·{\overline{\mathbf{V}}}^{pen}\end{array}\right\},\hfill \\ & \mathrm{s}.\mathrm{t}.,\left(2\right)\mathrm{and}\left(3\right),\left(5\right),\left(11\right)\mathrm{and}\left(12\right),\left(15\right)\mathrm{and}\left(16\right),\left(23\right)\left(25\right)–\left(27\right).\hfill \end{array}$$

(28)

Multipliers are updated as in (18), with the difference that the vector of constraint violations $\mathbf{R}$ is appended by $\left({\underline{\mathbf{V}}}^2-{\mathbf{V}}^{k-1}·{\mathbf{V}}^{{}^{\prime }}\right)$ and $\left({\mathbf{V}}^{k-1}·{\mathbf{V}}^{{}^{\prime }}-{\overline{\mathbf{V}}}^2\right)$, and the vector of multipliers $\Lambda$ is appended by ${\underline{\Lambda }}^{pen}$ and ${\overline{\Lambda }}^{pen}$. Multipliers corresponding to inequality constraints cannot be negative, and are projected onto a positive orthant. Piece-wise linear $l_1$-norms within (28) are linearized by using standard special ordered sets [22,23].

Feasibility. As $c^k$ increases, violation levels of relaxed constraints decrease. However, the total lack of constraint violations does not imply feasibility because, for example, linearized power flows do not coincide with original power flows. To ensure that ${\widehat{f}}_{l,t}^p\to f_{l,t}^p$ and ${\widehat{f}}_{l,t}^q\to f_{l,t}^q$, $c_p^k$ increases as in (21):

$$\begin{array}{cc}& c_p^k=c_p^{k-1}·{\beta }_p.\hfill \end{array}$$

(29)

The intention here is to force “$l_1$-proximal” terms to zero to ensure feasibility, after which $c_p^k$ is decreased as

$$\begin{array}{cc}& c_p^k=c_p^{k-1}·{\beta }_p^{-1}.\hfill \end{array}$$

(30)

After constraint violations and proximal terms are both zero, the feasible solution is obtained. Forcing the violations exactly to zero may require a considerable number of iterations, so tolerances $ϵ$ and ${ϵ}_p$ are introduced for constraint violations and proximal terms, respectively. The algorithm is shown in Algorithm 1 and the flow chart in shown in Figure 2 below.

Algorithm 1 Surrogate Lagrangian Relaxation.

Initialize ${\Lambda }^0$, $c^0$, $c_p^0$ and $s^0$ while stopping criteria are not satisfied 1 solve MSO problem 2 If the total constraints violations $>ϵ$ and total violation of proximal terms is $>{ϵ}_p$ update multipliers and increase c 3 If the total constraints violations $<ϵ$ and total violation of proximal terms is $>{ϵ}_p$ update multipliers, stop increasing c and increase $c_p$ 4 If the total constraints violations $<ϵ$ and total violation of proximal terms is $<{ϵ}_p$ reduce c and stop increasing $c_p.$

Stopping Criteria. Within the above Algorithm, the “if” condition in step $\mathbf{4}$ can by itself be used as a stopping criterion. Alternatively, after c is reduced per $\mathbf{4}$, the algorithm can continue with the multipliers update until a predefined CPU limit is reached.

5. Results

The method is implemented using IBM ILOG CPLEX Optimization Studio V 12.8.0.0 [24] on a PC with 2.40 GHz Intel Xeon E-2286M CPU and 32G RAM. To demonstrate the coordination aspect and the near-optimal performance of the method, a simple nine-bus system is considered.

5.1. System Description

Consider a nine-bus system with meshed topology as shown in Figure 3 below [25]. Assume that at buses 2 and 3 there is a generator, respectively; at bus 4 there is a solar farm of three stochastic states, and a tap changer; while at buses 5, 7, and 9, there is a load. The micro-grid is connected to the main grid through bus 1. The time interval is 15 min and the planning horizon is 3 h.

5.2. Results and Discussions

5.2.1. Droop Control

To demonstrate the performance of droop control, the micro-grid is disconnected from the grid at interval 12; thus turning the local operations into the islanded mode. By using the approach developed above, the problem (1)–(17) is solved, and a feasible solution is obtained in 202 seconds. The results show that the frequency stays within the range of [59 Hz, 61 Hz] due to effective droop control.

5.2.2. Tap-Changer

To show the impacts of tap changing, the problem in Section 5.2.1 is solved again without tap-change-related costs in the objective function. The results are compared in

Table 1. Results on tap-changing.

	Expected Number of Tap Changes	Solving Time
Without tap-changing cost	7	202 s
With tap-changing cost	0	262 s

below. It can be seen that the expected number of tap changes is greatly reduced when the corresponding cost is considered in the objective function.

5.2.3. Costs

To show the benefits of being connected to the main grid all the time, the problem in Section 5.2.1 is solved again under the grid-connected mode. The results are compared in

Table 2. Results of two different grid-connection scenarios.

	Expected Cost	Expected Number of Tap Changes	Solving Time
Grid connected until interval 11	$37,126	0	202 s
Grid connected all the time	$36,724	0	154 s

below. It can be seen that the total cost is reduced by connecting to the main grid, taking advantage of the time-varying costs of the grid power.

To provide insights into the convergence, the calculated trajectories of voltages at selected buses 2 and 5 at interval 12 are plotted on a complex plane in Figure 4. The feasible region delineated by constraints (4) is in-between and on the boundaries of the two concentric circles. Arrows schematically indicate the directions voltages move along to escape from infeasibility. It is demonstrated that voltages approach near-feasible solutions within a few iterations.

6. Conclusions

Motivation. Recent blackouts in California, New York, and Britain in 2019, which left millions of customers without power, pushed for an increase in the power system’s resilience, stability, and reliability. DERs, such as solar panels and wind turbines, have provided a promising solution. In order to prompt the applications of renewable energy or DERs, micro-grid is a promising paradigm. Informative Methodology. This paper addresses the issues of sustainability of micro-grid operations. Through the consideration of Markov processes to capture stochasticity of renewable generation; through tap-changer constraints penalizing frequent changes of taps, with which the lifespan of the expensive equipment will be greatly extended; through consideration of AC power-flow constraints appropriate for distribution systems; and through droop-control constraints for restoration of frequency and voltage after disconnection from the main grid. The above-mentioned constraints lead to non-convexities and difficulties solving the resulting micro-grid operation optimization problems. The methodology based on dynamic linearization and $l_1$-norm penalization is exact and is amenable to the use of MILP solvers. After linearization, the only source of non-convexities is the presence of discrete variables, which are efficiently handled by MILP solvers. Key Results. It is demonstrated that voltages approach feasible values within a few iterations, leading to exact satisfaction of AC power flows ensuring the feasibility of operations. As a result of the above, the overall performance of the method leads to convergence within a matter of minutes. It is also demonstrated that micro-grid operations benefit from the new solution methodology; specifically, the number of tap changes is drastically reduced, thereby ensuring higher sustainability in the presence of voltage fluctuations caused by uncertainties. Main Advantages over Previous Studies. The main advantage over previous studies is the consideration of key elements that lead to the feasibility of sustainability of micro-grid operations; namely, exact AC power flow, droop control, and tap-changer constraints are considered without simplification, and the methodology ensures their exact satisfaction. Limitations and Challenges of the Experimental Study. When considering multiple micro-grids, PV (or wind) generation at different locations cannot be aggregated and has to be modeled as a Markov chain per PV bus. The resulting global states are a large number of combinations of nodal states, leading to much-increased complexity. Suggestions for Future Research. With an increase in the number of micro-grids, which are expected to be interconnected, it is suggested that the methodology is extended to coordinate several micro-grids; on the bright side, upon the decomposition inherent in the method developed, the micro-grid complexity associated with both the number of discrete variables and the number of stochastic states can be significantly reduced. The present work: 1. has enabled efficient resolution of micro-grid subproblems and 2. paves the way for efficient coordination of micro-grids. The coupling constraints that connect adjacent micro-grids can be dualized by the Surrogate Lagrangian Relaxation to efficiently coordinate micro-grids for high resilience and reliability spanning larger geographical areas.