A Novel Interval Programming Method and Its Application in Power System Optimization Considering Uncertainties in Load Demands and Renewable Power Generation

Abstract

This paper expresses the output power of renewable generators and load demand as interval data and develops the interval economic dispatch (IED), as well as interval reactive power optimization (IRPO) models. The two models are generalized into a specific type of linear interval programming (LIP) and nonlinear interval programming (NLIP), respectively. A security limits method (SLM) is proposed to solve LIP and NLIP problems. As for the LIP, the maximum radii of the interval variables are first calculated by the optimizing-scenarios method (OSM) for defining security limits, and the LIP is transformed into deterministic linear programming (LP), for which its constraints are the security limits, which can be solved by the simplex method. As for the NLIP, Monte Carlo simulations were used to obtain the maximum radii of the interval variables, and the average interval ratio of the interval variables is defined to compute the security limits for transforming the NLIP to deterministic nonlinear programming (NLP), which can be solved by using the interior point method. Finally, the IED and IRPO are used to verify the effectiveness and engineering of the proposed SLM.

1. Introduction

The development and utilization of renewable energy are important measures for reducing carbon emissions and confronting the energy crisis [1]. However, there is high uncertainty in the output of wind power and photovoltaic power, as well as the load demand that will affect the stable operation of the electric power system (EPS) [2,3]. Therefore, it is necessary to study the economic and stable operation of the EPS by incorporating renewable power generation. Two important optimization problems in the EPS are active power economic dispatch and reactive power optimization [4,5]. The former dispatches the active power output of the generator so as to minimize the operation cost of the power grid [6]. The latter aims to minimize the active power loss so as to keep the voltage of the entire power grid within security constraints [7].

The aforementioned two EPS optimization models can be generalized into the LP and NLP models, which both consist of the objective function, decision variables, equality, and inequality constraints. The equality constraint is the power flow equation, which is based on a certain value of power injection in order to obtain the voltage and the line power flow [8,9,10]. The value of power injection becomes uncertain due to the fluctuation in power generation and load demand with the renewable energy connected to the EPS [11]. The original optimization problem is transformed into solving a special type of uncertain LP and NLP models; i.e., the value of power injection is uncertain on the right side of the power flow equation [12,13,14].

There are two main approaches for solving uncertain LP and NLP problems, i.e., stochastic programming (SP) and robust programming (RP) [15,16]. As for the SP, it is required that the probability distribution or the probability of scenario occurrence of the uncertain parameters should be given in advance. The aforementioned two types of data are counted and simulated based on massive historical data, which are quite difficult to achieve in practical engineering due to various reasons such as data confidentiality and the great cost of data collection [17]. In contrast, the RP uses uncertain sets to represent uncertainties in LP and NLP models without obtaining the probability distribution. However, the RP is available only for solving convex models, and the equality constraints posed by linear and nonlinear equations in the LP and NLP models should be convexified prior. This leads to the limited applicability of the RP [18].

To address the abovementioned problems associated with SP and RP approaches, interval programming (IP) is proposed and applied to solve the LP and NLP models with uncertainties. The IP represents uncertainties as intervals, and there is no need to assume the probability distribution function or uncertainty sets of the uncertainties [19]. A large number of research studies have been conducted on LIP and NLIP problems [20,21]. In the research of the LIP, Mráz [22] first introduced the concept of the LIP and emphasized that it is an NP-hard problem of a two-layer programming model because the linear equation is solved for each set of coefficients in the right coefficient vector interval and coefficient matrix interval. Mráz [23] solved LP problems with interval parameters by obtaining the maximum and minimum values of the objective function, and two extreme constraints were used for substitution when dealing with constraints. In detail, linear programming was constructed to find the upper and lower bounds of each interval variable, and it combined them into an interval vector; the shell of the solution set of linear interval equations was then obtained. In the research of the NLIP, Chen [24] studied the NLP model with all variables as intervals, and its objective and constraint functions are expanded into functions about the middle point and the radius of the state variable by using Taylor’s formula. Jiang [25] studied methods for solving NLP problems with interval parameters in both the objective function and the constraints. Wu [26] derived Wolfe’s duality theory in interval optimization, and the optimality conditions for the models of four types of NLIP problems were given. Wu [27] provided the Karush–Kuhn–Tucker (KKT) optimality condition with interval parameters in the objective function. Hosseinzade [28] proposed and proved the corresponding KKT optimality conditions theoretically for multi-objective optimization problems with interval parameters in both objective and constraint functions. However, existing IP methods still have some problems in solving LIP and NLIP models. As for the LIP, the existing research mainly focuses on the LIP model with a few interval parameters, and the sensitivity analysis method is generally chosen for solving the LIP model. Although sensitivity analyses perform well in solving the LIP model with few interval parameters, it is not effective in dealing with LIP problems and includes many interval parameters because of great approximation errors. With respect to the NLIP, KKT optimality conditions are not suitable for handling NLIP with wider ranges of interval parameters, which restricts the KKT optimality condition-based methods to solving a few NLIP models with narrow interval parameters. Meanwhile, the intelligent algorithms (e.g., the particle swarm optimization methods and the genetic algorithm) applied for solving the NILP problems are not applicable to large-scale NLIP problems due to their poor solution efficiency.

To overcome the aforementioned problems of the existent IP methods in solving LP and NLP models, the SLM is proposed for solving a special type of LIP and NLIP problems in this paper. Specifically, the present study proposes an SLM and a modified SLM (MSLM), which can transform the LIP and NLIP into deterministic models that can be easily solved using linear and nonlinear programming methods. As to the LIP problem, first, the radii of the interval variables are determined using the OSM to define security limits. Second, the LIP model is transformed into a deterministic LP model with security limits as constraints, which can be solved using the simplex method. As to the NLIP problem, first, Monte Carlo simulations are used in NLIP to obtain the radii of interval variables. Next, the position ratio of the interval variables is defined to compute the security limits for converting the NLIP to a deterministic model that can be solved by using interior point method. At the end of this paper, the proposed SLM and MSLM are separately tested using the interval economic dispatch and the interval reactive power optimization in the EPS. Our contributions to this paper are as follows:

• The SLM is proposed to solve a specific type of LIP and NLIP model with a determined system of equations. The SLM first defines the security limits to switch the LIP and the NLIP to deterministic LP and NLP models, respectively, and then the switched models are solved by traditional optimization approaches to obtain the optimal control variables (i.e., real-valued variable), which are able to ensure the feasibility of interval variables residing in the LIP and the NLIP models.
• The MLSM is established to reduce the conservatism of the SLM in solving the NLIP. To expand the feasible domain of the switched deterministic model in the SLM, the interval position rate is defined to redefine the less conservative security limits of the NLIP model. Meanwhile, a predictor-corrector procedure is applied here to ensure the effectiveness of the redefined security limits so as to guarantee the feasibility of the SLM and to reduce conservatism.
• Two cases in the EPS are tested to demonstrate the effectiveness of the SLM in solving the LIP and the MSLM in solving the NLIP. In detail, the SLM is applied to solve the IED model of the IEEE 118-bus test system, and the MSLM is utilized to solve the IRPO model of the IEEE 14-bus test system.

The rest part of the paper has the following structure. The mathematical models of the LIP and the NLIP are introduced in Section 2. Section 3 presents the optimization procedure of the SLM for the LIP and the NLIP. Section 4 describes the applications and the numerical simulations. Finally, conclusions are presented in Section 5.

2. State of the Problem

The LIP and NLIP models with a determined system of equations are formulated in this section. The interval is expressed in two kinds of forms as $[\underset{¯}{\mathit{A}},\overline{\mathit{A}}]$ and $[{\mathit{A}}^L,{\mathit{A}}^U]$, where $\underset{¯}{\mathit{A}}$ and ${\mathit{A}}^L$ denote the upper boundaries of the variables, and $\overline{\mathit{A}}$ and ${\mathit{A}}^U$ denote the lower boundaries. In addition, subtraction, multiplication, and division of interval operations are provided in the literature [29].

The LIP model is given as follows [30]:

$$\begin{array}{l}\min f(\mathit{X})\\ \mathrm{s}.\mathrm{t}.\{\begin{array}{l}\mathit{B}\mathit{X}+\mathit{C}\mathit{u}=[\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]\\ {\mathit{X}}^{\min }\le \mathit{X}\le {\mathit{X}}^{\max }\\ {\mathit{u}}^{\min }\le \mathit{u}\le {\mathit{u}}^{\max }\end{array},\end{array}$$

(1)

where

• $\mathit{X}$ is an n-dimensional state variable with an interval value;
• $\mathit{u}$ is the control variable including only the discontinuous real number value;
• $\mathit{B}$ and $\mathit{C}$ are coefficient matrices of $\mathit{X}$ and $\mathit{u}$, respectively;
• $f$ is the objective function, which is a linear function of $\mathit{X}$;
• $\mathit{B}\mathit{X}+\mathit{C}\mathit{u}=[\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]$ is the linear equality constraint;
  the values of the objective function are interval values instead of real numbers for each specific $\mathit{X}$.

The NLIP model is given as follows [25]:

$$\begin{array}{l}\min f(\mathit{X},\mathit{u})\\ \mathrm{s}.\mathrm{t}.\{\begin{array}{l}\mathit{h}(\mathit{X},\mathit{u})=[\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]\hfill \\ \begin{array}{l}{\mathit{X}}^{\min }\le \mathit{X}\le {\mathit{X}}^{\max }\\ {\mathit{u}}^{\min }\le \mathit{u}\le {\mathit{u}}^{\max }\end{array}\hfill \end{array},\end{array}$$

(2)

where

• $\mathit{X}$ is an n-dimensional state variable with an interval value;
• $\mathit{u}$ is the control variable with the discrete real number value;
• $f$ and $\mathit{h}$ are the objective function and constraint, respectively, and they are both nonlinear functions of $\mathit{X}$ and $\mathit{u}$;
  the values of the objective function are interval values instead of real numbers for each specific $\mathit{X}$.

3. Solutions of the LIP and NLIP

3.1. Solution of the LIP Using the SLM

To solve the LIP model, we proposed an SLM in this section. In detail, we firstly compute the maximum radii and midpoints of the interval variables. Then, the security limits of the interval variables are defined based on the results of the interval variable’s maximum radii. By performing these steps, the LIP model is transformed into a deterministic LP model, which is solved by the simplex method. Furthermore, the solution of the transformed deterministic LP is verified as a feasible solution of the LIP by using a simple example. Finally, the computational procedures of the proposed SLM for solving the LIP are provided at the end of this section.

3.1.1. Definition of the Security Limits for the LIP

For defining the security limits of the interval variables, it is required that the maximum radii should be known in advance. Notice that maximum radii and midpoints of the LIP’s interval variables are fixed because the equations of the LIP model are linear. Therefore, we have the following propositions.

Proposition 1.

The maximum radius of an interval variable $X_i$ is a constant in the LIP model denoted by Equation (1):

$$\Delta X_{i,\max }=\underset{{\mathit{u}}_{\min }\le \mathit{u}\le {\mathit{u}}_{\max }}{\max }\{\Delta X_i|\mathit{B}\mathit{X}+\mathit{C}\mathit{u}=[\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]\}=\mathrm{constant},$$

(3)

where $X_i$ is an interval variable of $\mathit{X}$ under a given ${\mathit{u}}_i$; $\Delta X_i$ denotes the radius of an interval variable $X_i$; $\Delta X_{i,\max }$ is the maximum radius of an interval variable ${\mathit{X}}_i$ under all $\mathit{u}\in [{\mathit{u}}_{\min },{\mathit{u}}_{\max }]$.

Proof of Proposition 1.

We first convert the constraint equation $\mathit{B}\mathit{X}+\mathit{C}\mathit{u}=[\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]$ of Equation (1) to the following:

$$\begin{array}{l}\mathit{X}=\left|{\mathit{B}}^{-1}\right|\cdot ([\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]-\mathit{C}\mathit{u})\\ =\left|{\mathit{B}}^{-1}\right|\cdot [\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]-\left|{\mathit{B}}^{-1}\right|\cdot \mathit{C}\mathit{u},\end{array}$$

(4)

where |·| represents an arithmetic for switching all elements of a matrix to their absolute value. Obviously, $\left|{\mathit{B}}^{-1}\right|\cdot \mathit{C}\mathit{u}$ is a real-valued vector, and $\left|{\mathit{B}}^{-1}\right|\cdot [\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]$ is an interval-valued vector. Thereafter, it is concluded that $\mathit{u}$ only affects the upper and lower boundary values of the interval variable $\mathit{X}$ and not the interval’s widths. That is, the radii of $\mathit{X}$ are only related to the radii of the interval $[\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]$, for which its values are constant. Therefore, the maximum radii of $\mathit{X}$ in the LIP are constant, and Proposition 1 is proven. □

Proposition 2.

The real-valued vector ${\mathit{X}}_{\mathrm{C}}$ calculated by the constraint equation $\mathit{B}{\mathit{X}}_{\mathrm{C}}+\mathit{C}{\mathit{u}}_0=(\underset{-}{\mathit{h}}+\stackrel{-}{\mathit{h}})/2$ is the midpoint of the interval-valued vector $\mathit{X}$.

Proof of Proposition 2.

According to Proposition 1, $\mathit{u}$ does not affect the widths of the interval variable $\mathit{X}$, and here, the value of $\mathit{u}$ is set as ${\mathit{u}}_0=\frac{{\mathit{u}}_{\min }+{\mathit{u}}_{\max }}{2}$. The solution of the deterministic constraint Equation (6) is denoted as ${\mathit{X}}_{\mathrm{C}}$, which is a real-valued vector. To show the relationship between $\mathit{X}$ and ${\mathit{X}}_{\mathrm{C}}$, we list the interval equation constraints and the deterministic equation constraints as follows, respectively.

$$\mathit{B}\mathit{X}+\mathit{C}\mathit{u}=[\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}],$$

(5)

$$\mathit{B}{\mathit{X}}_{\mathrm{C}}+\mathit{C}{\mathit{u}}_0=(\underset{-}{\mathit{h}}+\stackrel{-}{\mathit{h}})/2.$$

(6)

The difference between Equations (5) and (6) will produce the following.

$$\mathit{B}(\mathit{X}-{\mathit{X}}_{\mathrm{C}})+\mathit{C}(\mathit{u}-{\mathit{u}}_0)=[\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]-(\underset{-}{\mathit{h}}+\stackrel{-}{\mathit{h}})/2,$$

(7)

By setting $\mathit{u}={\mathit{u}}_0$, Equation (7) is rewritten as follows:

$$\mathit{B}(\mathit{X}-{\mathit{X}}_{\mathrm{C}})=[-\Delta \mathit{h},\Delta \mathit{h}],$$

(8)

$$\mathit{X}-{\mathit{X}}_{\mathrm{C}}=\left|{\mathit{B}}^{-1}\right|\cdot \left[-\Delta \mathit{h},\Delta \mathit{h}\right],$$

(9)

where $\Delta \mathit{h}=(\stackrel{-}{\mathit{h}}-\underset{-}{\mathit{h}})/2$, which represents the maximum radii of the input interval $[\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]$. The midpoint of $(\mathit{X}-{\mathit{X}}_{\mathrm{C}})$ is equal to zero, because the midpoint of $\left|{\mathit{B}}^{-1}\right|\cdot [-\Delta \mathit{h},\Delta \mathit{h}]$ is zero. Therefore, ${\mathit{X}}_{\mathrm{C}}$ is the midpoint of the interval variable $\mathit{X}$, and Proposition 2 is proven. □

According to Propositions 1 and 2, we define the security limits of the LIP as follows:

$$A{\mathit{X}}_{}^{\min }={\mathit{X}}_{}^{\min }+\Delta \mathit{X},$$

(10)

$$A{\mathit{X}}_{}^{\max }={\mathit{X}}_{}^{\max }-\Delta \mathit{X},$$

(11)

where $\Delta \mathit{X}$ is the vector composed of all interval variables’ maximum radii. $A{\mathit{X}}_{}^{\min }$ and $A{\mathit{X}}_{}^{\max }$ represent lower and upper security limits, respectively, and they will be adopted to define the state variables’ limits of the switched deterministic LP.

3.1.2. Transformation of the LIP into a Deterministic LP

Based on the results of the security limits, the next process is that the LIP model is transformed into a deterministic LP model, and the detailed transformation mainly includes two aspects. Firstly, the security limits obtained by Equations (10) and (11) are set as the new inequality constraint of the deterministic variable $\mathit{x}$ in the LP. Secondly, $(\underset{-}{\mathit{h}}+\stackrel{-}{\mathit{h}})/2$ is taken as the input value of the deterministic equality constraint. By performing these steps, the LIP model is transformed into a deterministic LP model:

$$\begin{array}{l}\min f(\mathit{x})\\ \mathrm{s}.\mathrm{t}.\{\begin{array}{l}\mathit{A}\mathit{x}+\mathit{B}\mathit{u}=(\underset{-}{\mathit{h}}+\stackrel{-}{\mathit{h}})/2\\ \mathit{A}{\mathit{X}}_{}^{\min }\le \mathit{x}\le \mathit{A}{\mathit{X}}_{}^{\max }\\ {\mathit{u}}^{\min }\le \mathit{u}\le {\mathit{u}}^{\max }\end{array},\end{array}$$

(12)

where $\mathit{x}$ is a solution of the deterministic constraint equation $\mathit{A}\mathit{x}+\mathit{B}\mathit{u}=(\underset{-}{\mathit{h}}+\stackrel{-}{\mathit{h}})/2$, for which its value is a real-valued vector. $(\underset{-}{\mathit{h}}+\stackrel{-}{\mathit{h}})/2$ is the midpoint value of the input interval $[\underset{-}{\mathit{h}},\stackrel{-}{\mathit{h}}]$, and $(\underset{-}{\mathit{h}}+\stackrel{-}{\mathit{h}})/2$ can be taken as the input value of the deterministic constraint equation because of Proposition 2. By solving Equation (12) by using the simplex algorithm, an optimal solution is obtained, and it has been proven as an optimal solution for the LIP by Propositions 1 and 2.

3.1.3. The Procedures of the SLM in Solving the LIP

After introducing the two processes in Section 3.1.1 and Section 3.1.2, the SLM is constructed and proposed for solving the LIP model. In detail, the solving process of LIP is divided into four steps:

• Step 1: Input data of the LIP and set the parameters of the SLM;
• Step 2: Calculate the radii of the interval state variable by the OSM [31] and compute the security limits by Equations (10) and (11);
• Step 3: Set the security limits as the constraints of the deterministic state variable and transform the LIP model (i.e., (1)) into a deterministic LP model (i.e., Equation (12));
• Step 4: Apply the simplex method to solve Equation (12) and obtain the control variables and state variables, which have been proven as an optimal solution for the LIP.

3.1.4. A Simple Example for Testing the Effectiveness of the SLM

To test the effectiveness and facilitate the understanding of the SLM, a simple example of a three-variable LIP model is introduced here. The expressions of this example are given as follows:

$$\begin{array}{l}\min X_2+X_3\\ \{\begin{array}{l}{x}_1+X_2+2X_3=[4,6]\\ x_1+2X_2+X_3=[8,10]\\ 0\le x_1\le 3\\ 2.5\le X_2\le 5\\ -1.5\le X_3\le 1\end{array},\end{array}$$

(13)

where $X_2$ and $X_3$ are the state variables with interval values, and $x_1$ is the control variable with a real-valued number.

We first obtain the maximum interval radii of $X_2$ and $X_3$ by the OSM, and the results are $\Delta X_2=1$ and $\Delta X_3=1$, which support Proposition 1. As shown in Figure 1a, the values of $X_2$ and $X_3$ enclose a central symmetric graph, and the four central symmetric figures’ center points can be obtained when the input value takes the midpoint (5, 9), which supports Proposition 2. According to the results of the maximum interval radii, the upper security limits of $x_2$ and $x_3$ are obtained by subtracting the maximum radii from the upper bounds of $X_2$ and $X_3$, while the lower limits are obtained by adding the maximum radii to the lower bounds. Based on the results of the security limits shown in Figure 1b, the original LIP model is transformed into a deterministic LP model:

$$\begin{array}{l}\min x_2+x_3\\ \{\begin{array}{l}{x}_1+x_2+2x_3=5\\ x_1+2x_2+x_3=9\\ 0\le x_1\le 3\\ 3.5\le x_2\le 4\\ -0.5\le x_3\le 0\end{array}.\end{array}$$

(14)

where $x_2$ and $x_3$ are deterministic state variables with real-valued numbers. Equation (14) can be solved by the simplex method on the platform MATLAB R2018b.

3.2. Solution of the NLIP Using the MSLM

As for the NLIP, Propositions 1 and 2 suitable for the LIP are not applicable anymore, because its nonlinearity will affect the radii and midpoint of the interval variables, which leads to the ineffectiveness of the security limits defined by Equations (10) and (11).

Hence, an MSLM is proposed to solve the aforementioned problems. The MSLM first defines a special type of security limit, denoted as the absolute security limits, for keeping the feasibility of the interval variables to deal with the nonlinearity of the NLIP. To reduce the conservativeness of absolute security limits, the proposed MSLM redefines the modified security limits by applying the average interval ratio in the absolute security limits, and a predictor-corrector procedure is used to correct the modified security limits and to ensure optimized interval variables within their original limits. By performing these steps, the NLIP model is transformed into a deterministic NLP model, which can be solved by the interior point method. The detailed processes of the proposed MSLM in solving the NLIP are introduced in the following section.

3.2.1. Definition of the Absolute Security Limits for the NLIP

To define the absolute security limits, we first determine the maximum radii of the interval variables; here, the maximum radius of $X_i$ is given as follows:

$$\Delta X_i^{\max }=\underset{{\mathit{u}}^{\min }\le \mathit{u}\le {\mathit{u}}^{\max }}{\max }\{\Delta X_i|\mathit{h}(\mathit{X},\mathit{u})=[{\mathit{h}}^L,{\mathit{h}}^U]\},$$

(15)

where $X_i$ is the $it\mathit{h}$ dimension of $\mathit{X}$, and $\Delta X_i^{\max }$ is the maximum radius of $X_i$, which can be obtained by using a Monte Carlo simulation under all control variables $\mathit{u}$ satisfying ${\mathit{u}}^{\min }\le \mathit{u}\le {\mathit{u}}^{\max }$. Obviously, there is a control variable ${\mathit{u}}^{(i)}$ satisfying equation $\mathit{h}(X_i,{\mathit{u}}^{(i)})=[{\mathit{h}}^L,{\mathit{h}}^U]$, which ensures that the radius of an interval variable $X_i$ is at the maximum. For the sake of simplicity, we denote $\stackrel{-}{\Delta \mathit{X}}$ as a vector consisting of the maximum radii of all the intervals variables in $\mathit{X}$:

$$\stackrel{-}{\Delta \mathit{X}}=(\Delta X_1^{\max },\Delta X_2^{\max },\dots ,\Delta X_p^{\max }),$$

(16)

where $p$ is the length of the vector $\mathit{X}$.

Based on the results of the interval variables’ maximum radii, the absolute security limits are defined as follows:

$$\mathit{A}{\mathit{X}}_{}^{\min }={\mathit{X}}_{}^{\min }+2\stackrel{-}{\Delta \mathit{X}},$$

(17)

$$\mathit{A}{\mathit{X}}_{}^{\max }={\mathit{X}}_{}^{\max }-2\stackrel{-}{\Delta \mathit{X}},$$

(18)

where $\mathit{A}{\mathit{X}}_{}^{\max }$ and $\mathit{A}{\mathit{X}}_{}^{\min }$ are the upper and lower absolute security limits, respectively, and they will be adopted to define the state variable’s limits of the switched deterministic NLP. To demonstrate the effectiveness of Equations (17) and (18), a simple derivation for defining the absolute security limits is given. Given a fixed ${\mathit{u}}^{(i)}$ and $\xi \in [{\mathit{h}}^L,{\mathit{h}}^U]$, the following equations are obtained:

$$\mathit{h}(\mathit{X},{\mathit{u}}^{(i)})=[{\mathit{h}}^L,{\mathit{h}}^U],$$

(19)

$$\mathit{h}(\mathit{x},{\mathit{u}}^{(i)})=\xi ,(\xi \in [{\mathit{h}}^L,{\mathit{h}}^U],$$

(20)

where $\xi$ is an arbitrary vector of $[{\mathit{h}}^L,{\mathit{h}}^U]$, and $\mathit{x}$ is a deterministic real-valued variable satisfying $\mathit{x}\in \mathit{X}$. Obviously, $[\mathit{x}-2\stackrel{-}{\Delta \mathit{X}},\mathit{x}+2\stackrel{-}{\Delta \mathit{X}}]\supseteq \mathit{X}$ due to $\stackrel{-}{\Delta \mathit{X}}\ge \Delta \mathit{X}$. If $[\mathit{x}-2\stackrel{-}{\Delta \mathit{X}},\mathit{x}+2\stackrel{-}{\Delta \mathit{X}}]$ satisfies ${\mathit{X}}_{}^{\min }\le [\mathit{x}-2\stackrel{-}{\Delta \mathit{X}},\mathit{x}+2\stackrel{-}{\Delta \mathit{X}}]\le {\mathit{X}}_{}^{\max }$, it is guaranteed that $\mathit{X}$ satisfies constraint ${\mathit{X}}_{}^{\min }\le \mathit{X}\le {\mathit{X}}_{}^{\max }$.

Intuitively, Figure 2 exhibits the upper and lower absolute security limits expressed by Equations (17) and (18), and the interval variable results after being optimized. As shown in Figure 2, the upper bound of the optimized interval variables $\mathit{X}$ exceeded the absolute security limit $A{\mathit{X}}_{}^{\max }$ but did not exceed limit ${\mathit{X}}_{}^{\max }$, and this is because $\mathit{x}$ obtained by the deterministic NLP model is within the absolute security limits. Observe also that the width between the upper (or lower) absolute security limit and the upper (or lower) bound ${\mathit{X}}^{\max }$ is $2\stackrel{-}{\Delta \mathit{X}}$, which ensures that $\mathit{x}$ ($\mathit{x}\in \mathit{X}$) in addition to the width of the optimized interval variable $\mathit{X}$ is within $[{\mathit{X}}_{}^{\min },{\mathit{X}}_{}^{\max }]$.

3.2.2. Reducing the Conservativeness of the Absolute Security Limits

Obviously, a conservative problem resides in the absolute security limits obtained by Equations (17) and (18), because $[\mathit{x}-2\stackrel{-}{\Delta \mathit{X}},\mathit{x}+2\stackrel{-}{\Delta \mathit{X}}]$ represents the worst-case fluctuation ranges of interval variables. To solve this problem, we first adopt the smaller radii $\Delta {\mathit{X}}_0$ instead of $\stackrel{-}{\Delta \mathit{X}}$ and then define a special parameter ${\mathit{k}}^I$ named as the “average interval ratio” to reduce the conservativeness of security limits in Equations (17) and (18). The modified security limits are given as follows:

$$S{\mathit{X}}_{}^{\min }={\mathit{X}}_{}^{\min }+2{\mathit{k}}_{}^I\cdot \Delta {\mathit{X}}_0,$$

(21)

$$S{\mathit{X}}_{}^{\max }={\mathit{X}}_{}^{\max }-2(1-{\mathit{k}}_{}^I)\cdot \Delta {\mathit{X}}_0,$$

(22)

where $S{\mathit{X}}_{}^{\max }$ and $S{\mathit{X}}_{}^{\min }$ are the upper and lower modified security limits, respectively, and the division operation is a point division operation of vectors. $\Delta {\mathit{X}}_0$ denotes the radii of ${\mathit{X}}_0=[\underset{-}{{\mathit{X}}_0},\stackrel{-}{{\mathit{X}}_0}]$, and their values are obtained when control variables are set at the midpoint of the variation ranges, i.e., ${\mathit{u}}_0=({\mathit{u}}^{\min }+{\mathit{u}}^{\max })/2$ and $\mathit{h}({\mathit{X}}_0,{\mathit{u}}_0)=[{\mathit{h}}^L,{\mathit{h}}^U]$. ${\mathit{k}}^I$ denotes the average interval ratio, and its expressions are as follows:

$${\mathit{k}}^I=\frac{{\mathit{x}}_0-\underset{-}{{\mathit{X}}_0}}{\stackrel{-}{{\mathit{X}}_0}-\underset{-}{{\mathit{X}}_0}},$$

(23)

where ${\mathit{x}}_0$ denotes the solutions of the equations $\mathit{h}({\mathit{x}}_0,{\mathit{u}}_0)=({\mathit{h}}^L+{\mathit{h}}^U)/2$. Notice that ${\mathit{x}}_0\ne (\underset{-}{{\mathit{X}}_0}+\stackrel{-}{{\mathit{X}}_0})/2$ under most conditions because $\mathit{h}({\mathit{x}}_0,{\mathit{u}}_0)$ is a nonlinear equation. ${\mathit{k}}^I$ complies with the constraint $0\le {\mathit{k}}^I\le 1$, and ${\mathit{k}}^I=0.5$ when ${\mathit{x}}_0=(\underset{-}{{\mathit{X}}_0}+\stackrel{-}{{\mathit{X}}_0})/2$.

Although Equations (21) and (22) are defined to limit the control variables, there still exists a situation in which the optimized interval variable $\mathit{X}$ exceeds constraints $[{\mathit{X}}_{}^{\min },{\mathit{X}}_{}^{\max }]$. This problem occurs mainly because the average interval ratio is fixed at ${\mathit{u}}_0=({\mathit{u}}^{\min }+{\mathit{u}}^{\max })/2$, while the value of ${\mathit{k}}^I$ is not the same under different values of $\mathit{u}$. Therefore, a predictor-corrector procedure is applied to solve the aforementioned problem by squeezing the original modified security limits, and the corrected security limits are given as follows:

$$S^{\prime }{\mathit{X}}_{}^{\min }=S{\mathit{X}}_{}^{\min }+{\delta }_l\mathit{X},$$

(24)

$$S^{\prime }{\mathit{X}}_{}^{\max }=S{\mathit{X}}_{}^{\max }-{\delta }_u\mathit{X},$$

(25)

where $S^{\prime }{\mathit{X}}_{}^{\min }$ and $S^{\prime }{\mathit{X}}_{}^{\max }$ are the corrected lower and upper security limits. ${\delta }_u\mathit{X}$ and ${\delta }_l\mathit{X}$ are the magnitude of the interval variable $\mathit{X}$ exceeding the upper security limit and the lower security limit, respectively, and the values of ${\delta }_u{X}_i$ or ${\delta }_l{X}_i$ are set as 0 when an interval variable $X_i$ in $\mathit{X}$ does not exceed upper or lower constraints. Notice that the predictor-corrector operation will repeat until the optimized interval variable $\mathit{X}$ is within $[{\mathit{X}}_{}^{\min },{\mathit{X}}_{}^{\max }]$.

3.2.3. Transformation of the NLIP into a Deterministic NLP

To solve the NLIP problem, the NLIP model should be transformed into an NLP model, and the transformation is completed by the following three steps. Firstly, the modified security limits obtained by Equations (21) and (22) are set as new constraints of the deterministic state variable $\mathit{x}$ in the NLIP. Secondly, $({\mathit{h}}^L+{\mathit{h}}^U)/2$ is set as the input value of the equality constraint in the deterministic NLP. Finally, the NLIP model is transformed into a deterministic NLP model, which is expressed as follows:

$$\begin{array}{l}\min f(\mathit{x},\mathit{u})\\ \mathrm{s}.\mathrm{t}.\{\begin{array}{l}\mathit{h}(\mathit{x},{\mathit{u}}_0)=({\mathit{h}}^L+{\mathit{h}}^U)/2\\ S{\mathit{X}}_{}^{\min }\le \mathit{x}\le S{\mathit{X}}_{}^{\max }\\ {\mathit{u}}^{\min }\le \mathit{u}\le {\mathit{u}}^{\max }\end{array},\end{array}$$

(26)

where ${\mathit{u}}_0=({\mathit{u}}^{\min }+{\mathit{u}}^{\max })/2$, and $f(\mathit{x},\mathit{u})$ is assumed to be the midpoint value of $f(\mathit{X},\mathit{u})$ [32]. It is verified by a large number of arithmetic cases that the difference between $f(\mathit{x},\mathit{u})$ and the actual midpoint value of $f(\mathit{X},\mathit{u})$ is small [2]. Equation (26) is an NLP problem that can be solved by the interior point method [33]. Notice that new constraints of $\mathit{x}$ in Equation (26) are constructed from corrected security limits $S^{\prime }{\mathit{X}}_{}^{\min }$ and $S^{\prime }{\mathit{X}}_{}^{\max }$, when $\mathit{X}$ obtained by the deterministic NLP model exceeds $[{\mathit{X}}_{}^{\min },{\mathit{X}}_{}^{\max }]$.

3.2.4. The Procedures of the MSLM in Solving the NLIP

Based on the aforementioned introduction of MSLM in solving the NLIP model, its solution procedures are introduced as follows:

• Step 1: Input data of the NLIP and set parameters of the interior point method. The parameters of the interior point method include the convergence precision and the central parameters.
• Step 2: Obtain the maximum radius $\stackrel{-}{\Delta \mathit{X}}$ of the interval variable $\mathit{X}$ under all control variables ${\mathit{u}}^{\min }\le \mathit{u}\le {\mathit{u}}^{\max }$ by the Monte Carlo simulation. Compute the average interval ratios by (23) and the modified security limits by Equations (21) and (22). Here, values of the deterministic state variables are calculated by solving the equation $\mathit{h}({\mathit{x}}_0,{\mathit{u}}_0)=({\mathit{h}}^L+{\mathit{h}}^U)/2$ when computing the average interval ratio by Equation (23).
• Step 3: Transform the NLIP model into a deterministic NLP model, which is solved by the interior point method, and obtain the optimal control variables. Use the OSM to calculate the intervals of the state variables $\mathit{X}$ and the value of the objective function under the optimal control variables.
• Step 4: Judge whether state variables satisfy the constraints or not. If the interval variable $\mathit{X}$ does not exceed the limits, the optimization result is printed out and the procedures are completed. Otherwise, compute the amount by which the interval variables exceed the modified security limits and correct the security limits by Equations (24) and (25), and then switch to step 3.

According to the above procedures, the flow chart of the MSLM for solving the NLIP model is described in Figure 3.

4. Simulation Results

In this section, the IED model and the IRPO model in the EPS are established and tested to demonstrate the effectiveness and applicability of the proposed SLM and MSLM. In the first case, the SLM is applied to solve the IED model of the IEEE 118-bus test system. In the second case, the MSLM is utilized to solve the IRPO model of the IEEE 14-bus test system.

4.1. Simulation Results of the SLM in Solving the IED Model

4.1.1. The IED Model

The reactive power demand of loads in the IED model comprises intervals denoted as $[\underset{-}{P_{Li,t}},\stackrel{-}{P_{Li,t}}]$, and then the IED model of the LIP problem can be given as follows [34]:

$$\min \mathit{C}(P_G)={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in S_G}(a_i{P}_{Gi,t}^{}+c_i}}),$$

(27)

s.t.

$$\{\begin{array}{l}-{\displaystyle \sum _{j\in S}{B}_{ij}{\theta }_{j,t}=P_{Gi,t}-}[\underset{-}{P_{Li,t}},\stackrel{-}{P_{Li,t}}]\\ P_{Gi}^{\min }\le P_{Gi,t}\le P_{Gi}^{\max }\\ -P_{d,Gi}\le P_{Gi,t}-P_{Gi,t-1}\le P_{u,Gi}\\ -P_{ij}^{\min }\le B_{ij}({\theta }_{i,t}-{\theta }_{j,t})\le P_{ij}^{\max }\end{array},$$

(28)

where $S_G$ is the set of thermal power electric generators. $a$ and $c$ are the cost coefficients of the generators. $T=24 h$, and $P_{Gi,t}$ denotes the active power output of the generators. $B_{ij}$ and ${\theta }_{j,t}$ are the bus admittance matrix’s elements and the phase angle, respectively. $P_{d,Gi}$, $P_{u,Gi}$, $P_{ij}^{\max }$, and $P_{ij}^{\min }$ are the downward climb rate, upward climb rate, the maximum transmitted power, and the minimum transmitted power, respectively. All variables of the IED are classified as state variables (denoted by $\mathit{X}$, which include ${\theta }_{j,t}$ and $P_{Gi,t}$) and control variables (denoted by $\mathit{u}$, which include the active power output of the generators). Then, the IED model can be simplified as Equation (1), and it can be solved by the procedures for solving the LIP model proposed in Section 3.1.

4.1.2. Simulation Result

This section discusses the results of solving the IEEE 118-bus test system [35]. As shown in Figure 4, the IEEE 118-bus system contains 54 generators, 64 load buses, 178 branches including 169 transmission lines and 9 transformers, and 9 capacitors. Parameters of nine capacitors are revealed in

Table 1. The capacitor parameters of the IEEE 118-bus test system.

Control Variable Type	Location	Lower Range	Upper Range	Step Size
Reactive compensation capacitor	Bus 5	−0.5	0	0.1
	Bus 17	0	0.1	0.02
	Bus 37	0	0.2	0.04
	Bus 44	−0.3	0	0.06
	Bus 45	0	0.2	0.04
	Bus 48	0	0.2	0.04
	Bus 79	0	0.3	0.06
	Bus 82	0	0.3	0.06
	Bus 83	0	0.2	0.04

, and transformer ratios are all limited within the interval [0.9, 1.1], with a variation step of 0.05. A tolerance of ±10% is assumed in the active load power demand. Input data mainly include generator parameters, load parameters, transformer parameters, and line parameters. Moreover, all variables are valued in the per-unit system.

The results of the SLM when solving the IED model for the IEE-118 bus system are presented in Figure 5, Figure 6 and Figure 7. Figure 5 shows the power flow of 178 branches at 24 time periods before and after optimization. The transmission power constraint for each branch at different periods is calculated by the OSM [36], and the upper and lower limits of line transmission power are ±1.8 MW [2]. In Figure 5a, it is obvious that not all branches satisfy the upper and lower limits of transmission power before the IED optimization. In contrast, the optimized transmission power of each branch at 24 time periods is obtained in Figure 5b, which shows that all optimized line transmission power satisfies the upper and lower constraints. Figure 6 represents the total load curve at 24 time periods, and the load is divided into three types, i.e., peak load, valley load, and flat load, according to the value of total load in different periods of the day. The 24 periods of a day were divided into valley periods, peak periods, and flat periods according to the load value shown in Figure 6. Figure 7 illustrates the generator’s active output power at peak and valley periods and the generator’s cost factors. Valley periods and peak periods are selected to analyze scheduling effects because the economic dispatch is mainly for cutting the peak load and valley filling. In the valley load period, the output of four generators is relatively large, while in the peak load period, two other generators have relatively large output power in addition to the above four units. This is because the cost coefficients of the aforementioned six generators are relatively small.

In summary, the power dispatch schedule obtained by the proposed SLM is able to guarantee the transmission power of each branch within the upper and lower limits under the interval uncertainty, and the cost of power generation is minimized. Therefore, it is concluded that the proposed SLM is effective in solving the LIP model.

4.2. Simulation Results of the MSLM in Solving the IRPO Model

4.2.1. The IRPO Model

The IRPO model only considers the uncertainties of the input data in the active power generation and load demand, which are described as intervals, i.e., $[{P_{Gi}}^L,{P_{Gi}}^U]$, $[{P_{Li}}^L,{P_{Li}}^U]$, and $[{Q_{Li}}^L,{Q_{Li}}^U]$. Note that the active and reactive power loads of slack and generator buses are assumed as deterministic values, because they are mainly derived from power plants. Thus, the IRPO model can be expressed as follows:

$$\min P_{loss}={\displaystyle \sum _{i\in S}{\displaystyle \sum _{j\in S}{V}_i{V}_j{G}_{ij}\cos {\theta }_{ij}}},$$

(29)

s.t.

$$\{\begin{array}{l}[{P_{Gi}}^L,{P_{Gi}}^U]-[{P_{Li}}^L,{P_{Li}}^U]-V_i{\displaystyle \sum _{j\in S}{V}_j(G_{ij}\cos {\theta }_{ij}+B_{ij}\sin {\theta }_{ij})}=0,i\in {S^{\prime }}_G\\ Q_{Gi}-[{Q_{Li}}^L,{Q_{Li}}^U]-V_i{\displaystyle \sum _{j\in S}{V}_j(G_{ij}\sin {\theta }_{ij}-B_{ij}\cos {\theta }_{ij})}=0,i\in {S^{\prime }}_G\end{array},$$

(30)

$$\{\begin{array}{c}-[{P_{Li}}^L,{P_{Li}}^U]-V_i{\displaystyle \sum _{j\in S}{V}_j(G_{ij}\cos {\theta }_{ij}+B_{ij}\sin {\theta }_{ij})}=0,i\in S_L\\ Q_{Ci}-[{Q_{Li}}^L,{Q_{Li}}^U]-V_i{\displaystyle \sum _{j\in S}{V}_j(G_{ij}\sin {\theta }_{ij}-B_{ij}\cos {\theta }_{ij})}=0,i\in S_L\end{array},$$

(31)

$$\{\begin{array}{l}{V}_i^{\min }\le V_i\le V_i^{\max },i\in S_G\cup S_L\\ Q_{Gi}^{\min }\le Q_{Gi}\le Q_{Gi}^{\max },i\in S_G\\ T_l^{\min }\le T_l\le T_l^{\max },l\in S_T\\ Q_{Ci}^{\min }\le Q_{Ci}\le Q_{Ci}^{\max },i\in S_C\end{array},$$

(32)

where $S$, $S_G$, ${S^{\prime }}_G$, $S_L$, $S_C$, and $S_T$ are the indexed sets of system buses, generator buses, non-balance generator buses, load buses, load buses with static VAR compensators, and transformer branches, respectively. Equation (29) represents real power losses; Equations (30) and (31) are power flow equations; Equation (32) represents the security and system operational constraints, which include the limits of the reactive power generation of generators, static VAR compensators outputs, voltage magnitudes (except for the slack bus), and transformer ratios. ${\theta }_i$, $V_i$, $Q_{Ci}$, and $Q_{Gi}$ are the bus angle, voltage magnitude, reactive power compensator output, and reactive power generation, respectively. ${\theta }_{ij}={\theta }_i-{\theta }_j$, and $T_l$ is the ratio of the $l^{th}$ transformer. $G=\left\{G_{ij}\right\}$ and $B=\left\{B_{ij}\right\}$ are real and imaginary components of the admittance matrix [26].

All variables of the IRPO are classified as state variables and control variables. In detail, the state variable $\mathit{X}$ includes load bus voltage magnitudes, the reactive power generation of generator buses, the bus angles of non-balance buses, and the active power generation of the balance generator (assuming that the active output constraint of the balance generator is not considered). The control variable $\mathit{u}$ includes transformer ratios, reactive power compensations, and generator terminal voltages. Then, the IRPO model can be simplified as Equation (2), and it can be solved by the procedures for solving the NLIP proposed in Section 3.2.

4.2.2. Simulation Results

This section tests the convergence effect of the algorithm under different injected power data fluctuations to verify the engineering application capability of the proposed MSLM, and the IEEE 14-bus system shown in Figure 8 is adopted. The IEEE 14-bus system includes five power generators, one capacitor, and three transformers. The voltage magnitude of the generator buses and load buses are limited in [0.9, 1.1] and [0.95, 1.05], respectively. The transformer ratio is considered to be limited to [0.9, 1.1] with a variation step of 0.05, as shown in

Table 2. The capacitor and transformer parameters of the IEEE 14-bus test system.

Control Variable Type	Location	Lower Range	Upper Range	Step Size
Reactive compensation capacitor	Bus 9	0	0.5	0.1
Transformer ratio	Branch 4–7	0.9	1.1	0.05
	Branch 4–9	0.9	1.1	0.05
	Branch 5–6	0.9	1.1	0.05

. The reactive power generation variables are all constrained to interval [−2.0, 3.0], and initial penalty factors of the penalty function for the interior point method to handle discrete variables are set to $v_{bC}^{(0)}=5$, $v_{bT}^{(0)}=2$, and $V_{b\max }=500$. Tolerances of ±10%, ±15%, and ±20% are assumed in the load power demand and active power generation to verify the ability of the algorithm in handling variations in the input data. The iteration’s precisions of the interior point method are set as $\epsilon ={10}^{-6}$. All variables are valued in the per-unit system.

The results of the MSLM when solving the IRPO model for the IEE-14 bus system are presented in Figure 9, Figure 10 and Figure 11. Figure 9 shows midpoints of real power loss with different tolerances. Observe that the midpoints of real power losses are the smallest when the input data tolerance is ±10%, and this is because the feasible regions of the state variable, i.e., the widths of the security limits, is the largest among the three input data tolerances. In detail, the widths of the state variable’s security limits depend on the interval variables’ radii, and the radius is the smallest when the fluctuation interval of the input power data is ±10%. Figure 10 represents the security limits and optimized interval bounds of the load voltage magnitudes at different input data tolerances, while Figure 11 represents the reactive power generation of the generators. It is shown that the optimized intervals of the state variables are all within the limits. Observe also that the larger the data tolerance, the narrower the range of the security limits. In other words, the smaller the fluctuation range of the injected power data, the wider the security limits of the state variables.

These results indicate that all variable intervals are subject to their security limits, which demonstrates that the proposed MSLM is effective. Therefore, it is concluded from the above analysis that the present MSLM has good applicability to large-scale grids, and it provides a solution with good robustness relative to various tolerances of input data.

Furthermore, the EPS, which has different sizes, is tested to evaluate the efficiency of the MSLM, including the IEEE 14-bus system, IEEE 30-bus system, IEEE 57-bus system, IEEE 118-bus system, IEEE 300-bus system, and C703-bus system (from an actual system of Southern Power Grid in China). Then, the execution time of the algorithm and the time to obtain the maximum radius are counted shown in

Table 3. The IRPO algorithm based on the MSLM to test the calculation time of different EPS.

Test Systems	Algorithm Execution Time (s)	The Computation Time of the Maximum Interval Radii of the State Variables (s)
IEEE14	1.55	1.08
IEEE30	8.4	6.4
IEEE57	57.9	55
IEEE118	502.7	500
IEEE300	3002	2990
C703	23,250	23,200

, while the number of Monte Carlo simulations in the algorithm during the test is the same, which is 10. The computation time of the algorithm increases sharply when the number of system buses increases to 300, as observed in Table 1. The fact causing drastic changes in computation time is that the algorithm uses the interval power flow algorithm to obtain the maximum radius of the state variables, and the computation time of the interval power flow algorithm is closely related to the number of buses in EPS [31]. Thus, it is necessary to investigate and reduce the computation time needed when obtaining the interval maximum radii of the state variables and to meet the requirement of engineering online applications.

As far as the computational information is concerned, all tests are realized using MATLAB R2018b on a CPU i5-2.4GHz with 16 GB of RAM and eight cores.

5. Conclusions

In this paper, the SLM is proposed to solve a specific type of interval linear and nonlinear programming problems with equality constraints posed by linear and nonlinear equations, respectively. As for the LIP model, the radii of the interval variables are first calculated by the OSM for defining security limits, and the LIP model is transformed into a deterministic LP model for which its constraints include security limits, which can be solved using the simplex method. As to the NLIP, the Monte Carlo simulation is used to obtain the radii of the interval variables, and the position ratio of the interval variables is defined to compute the security limits for transforming the NLIP model to a deterministic NLP model that can be solved by the interior point method. Furthermore, the IED and the IRPO in the EPS are used to verify the effectiveness of the proposed SLM and MSLM. Firstly, the numerical simulation results of the IED solved by the SLM are obtained for the IEEE 118-bus test system, and the proposed SLM ensures that the line power flow is within constraints. Similarly, the analysis demonstrates the advantages and effectiveness of the proposed SLM. Secondly, the numerical simulation consequences of the IRPO with different input data tolerances are gained by the proposed MSLM for the IEEE 14-bus test system, which shows that midpoints of real power losses are the smallest when the input data tolerance is ±10%. It is concluded that the MSLM algorithm is able to solve the large-scale EPS with different input power data tolerances and has good robustness and engineering application capability.