A New Solution for Stochastic Optimal Power Flow: Combining Limit Relaxation with Iterative Learning Control

Abstract

A stochastic optimal power flow (S-OPF) model considering uncertainties of load and wind power is developed based on chance constrained programming (CCP). The difficulties in solving the model are the nonlinearity and probabilistic constraints. In this paper, a limit relaxation approach and an iterative learning control (ILC) method are implemented to solve the S-OPF model indirectly. The limit relaxation approach narrows the solution space by introducing regulatory factors, according to the relationship between the constraint equations and the optimization variables. The regulatory factors are designed by ILC method to ensure the optimality of final solution under a predefined confidence level. The optimization algorithm for S-OPF is completed based on the combination of limit relaxation and ILC and tested on the IEEE 14-bus system.

1. Introduction

Deterministic approaches to optimal power flow (OPF) problems have been well developed and widely used in power systems [1-2]. However, as renewable energy, such as wind power, has developed at a fast speed in recent years, it is of great significance to study the effects of uncertainty on the optimal power flow of the system containing renewable power supplies and other uncertain factors existing in power systems. These uncertainties mainly come from changes in network configuration, outages of system elements and forecast inaccuracies [3]. At present, the studies on OPF considering the stochastic factors can be classified as probabilistic optimal power flow (P-OPF) and stochastic optimal power flow (S-OPF) [4]. P-OPF can obtain the probability distribution functions (PDFs) of state variables from the probability distributions of loads and other uncertain factors by Monte-Carlo simulation [5], cumulant method [6], point estimate method [7], etc. The uncertain factors generally do not affect the optimal results of P-OPF. S-OPF differs from P-OPF in that it considers the uncertainty as part of the model’s constraints and objective [4]. Thus, uncertain factors affect the calculation process and the final results of S-OPF. Present studies on S-OPF also include two kinds: one formulates the model as deterministic OPF (or separates the model into sub-problems as deterministic OPF) and mainly considers the uncertain changes in network configuration and component failures [8], the other usually models the S-OPF based on chance constrained programming (CCP) [9].

The advantage of CCP is the comprehensive consideration between the profitability and reliability under uncertainty using chance constraint as inequality constraint satisfied with a predefined probability value to ensure a certain level of reliability. Recently, CCP-based methods have been applied to power system under uncertainty for transmission network expansion planning [10], unit commitment [11], evaluation of available transfer capability [12], S-OPF [13], etc. For the CCP-based model with constraints including power flow equations and nonlinear chance constraints, different solving strategies have been developed. Hui Zhang et al. [13] solved the formulated CCP-based OPF under load power uncertainties by a back-mapping approach and a linear approximation of the nonlinear model equations. Wang, Shuang et al. [14] transformed the chance constrains using sample average approximation (SAA) and solved the optimal model by genetic algorithm (GA). Qianfan Wang et al. [15] also used a combined SAA algorithm to solve the CCP-based problem. The disadvantage of SAA-based method is that it needs lots of samples to ensure the solution converges to the optimal one. Zhipeng Liu et al. [16] employed a Monte Carlo simulation-embedded genetic-algorithmbased approach to solve the developed CCP model, which is rather computationally expensive for larger systems. Schellenberg, A. et al [6] and Z. Hu et al [17] combined the statistical step with the pure Newton step in a linear fashion to produce a step that was applied to optimization variables.

In this paper, the S-OPF model is formulated based on CCP with stochastic inputs including load fluctuation and the forecast error of wind farm output. The solution of the formulated S-OPF model is based on the classic interior point method of deterministic OPF with high convergence rate and less calculation. Since the S-OPF is a high-dimensional nonlinear model with a large number of uncertain inputs as well as multiple probabilistically constrained outputs, it is difficult to deal with the nonlinear chance constraints of CCP directly in S-OPF model. Iterative learning control (ILC) is used to regulate the limits of chance constraints in this paper. This approach is very attractive in control area, because of its performance improving capabilities along repetitive trajectories, often met in practice [18]. A limit relaxation approach combined with ILC is proposed to solve the formulated CCP-based S-OPF model. The proposed method can deal with multiple chance constraints simultaneously in the optimization process.

The remainder of this paper is organized as follows. In Section 2, the S-OPF model based on CCP is formulated and the stochastic inputs are described. Section 3 proposes an optimization algorithm for the S-OPF model. Results of case study are given in Section 4. The conclusion is presented in Section 5.

 

2. Model of S-OPF

2.1 OPF model based on CCP

Ignoring the wind power operation costs, the objective of the power system economic dispatch is to minimize the expected power generation cost of the system’s conventional generator sets. The deterministic OPF can be established as following [1-2]:

s.t

The objective is a quadratic cost function. It can also be the transmission losses, economic benefits, the environmental benefits, etc. Power flow Eqs. (2-3) can be rewritten as

where, X = [PG ,QG ,θ ,V ].

The conventional unit output constraint (4), generator ramp rate constraint (5), spinning reserve constraints (6), reactive power constraint of reactive source (7), nodal voltage constraint (8), and branch flow constraint (9) can be represented in the manner of vector as:

Inequality constraints (11) can be further expressed as:

where, H(X) = [H' (X) − hmax , hmin − H' (X)]T .

The deterministic OPF (D-OPF) is established with (1) as the objective and (10) and (12) as the constraints. However, the solution obtained with the deterministic OPF which does not consider stochastic factors may lead to inferior or even infeasible decision under stochastic environment. When considering the uncertainties of load and wind power, PW , PD, QW and QD are stochastic variables. Therefore, H(X ) is also stochastic according to (2) and (3). Chance constrained method is considered as a useful tool to evaluate the probability or reliability of holding boundary conditions [10-13]. In CCP-based OPF model, the inequality constraint is replaced by the chance constraint which is defined in the form of probability as:

The chance constraint requires the inequality should be satisfied at a predetermined level of probability βn ( 0 ≤ βn ≤ 1 ). βn = 1 means that the inequality should always keep true in the supposed stochastic conditions. The chance constraint does not work if its confidence level is zero. The elements of β can be set as the same or defined separately according to different reliability requirements of inequality constraints. As a result, the OPF under uncertainty can be reformulated into a chance constrained OPF problem with (1) as the objective function subjected to constraints (2)-(3) and (13).

2.2 Stochastic Inputs

The stochastic inputs discussed in this paper include the uncertainty of wind farm outputs and loads. The load forecast is more precise than wind power forecast of a wind farm. Therefore, we focus on analyzing the uncertainty of wind forecast here. The analysis can also be applied to describe loads.

Based on a specific forecast method, the actual wind power output can be divided into two parts: the forecasted wind power (deterministic part) and the forecast error (uncertain part), as shown in Eq. (14). From (14), the PDF of the actual wind power can be obtained by the probability distribution of the forecast error.

The forecast methods of includes: the physical forecast based on weather information, the forecast method based on historical data, the combined forecast method, etc. [19, 20]. Due to the random fluctuation of wind speed as well as the precision and adaptability of the forecast methods, precise forecast of wind power is impossible. Gaussian distribution is usually adopted to describe the distribution of forecast error in analysis. However, it has already been found that the distribution of actual forecast error has the characteristic of skewing either to the left or to the right. Bludszuweit, H. et al [21] and Fabbri, A. et al [22] have proposed that Gaussian distribution is not valid to describe the wind power forecast error and the Beta distribution gives more reasonable and accurate results. In this paper, cumulant is used to describe the probabilistic characteristics of stochastic inputs no matter which distribution the stochastic inputs should fit. According to the history data of wind power forecast, the cumulants of the measured forecast error can be estimated by maximum likelihood estimation (MLE) [23].

Usually, correlations of these stochastic inputs may exist [14, 24]. Correlations can be included in the model under study, but it is not the main concern in this paper. The correlated stochastic inputs can be represented as (15) applying a Gram–Schmidt orthogonalization.

The elements of ΔW is . ρ1 ρ2 … ρn is a group of uncorrelated variables.

 

3. A Combined Solution

3.1 Limit relaxation

The difficulty in solving the CCP-based S-OPF is the processing of chance constraints as shown in (13). Due to complicated nonlinearity of the chance constraints, the scope of X is impossible to be expressed explicitly in mathematical formula. Therefore, it is hard to adjust the optimization variables directly. This paper proposes a limit relaxation method to regulate the optimization variables indirectly, taking advantage of the correlation between optimization variables and constraints. The limit relaxation introduces regulatory factor a ( a > 0 ) as the margin to the limits of the inequality constraints to narrow the solution space. Usually, the optimal solution achieves (or closes to) one (or more than one) of the limits of the inequality constraints in deterministic OPF. The regulatory factor defined in advance can prevent the optimal solution from reaching the limits of the inequality constraints. And it is regulated iteratively by ILC, described in section 3.3, until all of the chance constraints are satisfied. Therefore, the problem of directly calculating the optimization variables is transformed into a problem on how to narrow the solution space to ensure the optimal solution of the formulated S-OPF.

Regulatory factor is added to the inequality constraints as shown in (16).

Eq. (16) can also be represented as

When a is added to the deterministic OPF, the OPF model applying barrier function can be written as follows [2]:

The interior point method can be adopted to solve thisdeterministic model. A given regulatory factor a is added in Eq. (20) to relax the limits, which is different from the interior point method in the traditional optimal power flow.

The Lagrange function of the optimization model (18)-(21) is as following:

The Newton iteration step of the Karush–Kuhn–Tucker (KKT) equation of (22) is

where, , [μ] is used to denote a diagonal matrix with vector μ on the diagonal. Eq. (23) shows the relationship of H(X ) , ΔX and a . The regulatory factor can indirectly change theoptimization variables through (23). The output of the optimization model (18)-(21) includes: X , Z , γ , λ and μ .

3.2 Calculation of hβ

The Cornish-Fisher series is developed based on cumulant theory to approximate the probability density function. The principles of cumulant and Cornish-Fisher series can be referred in [25] and [26].

The linearized representation of Hm is as:

The uncertainty of the inequality constraints can be represented as:

The Jacobian matrix Jm at each iteration can be derived from the elements of GX when the deterministic optimization process ends. According to the cumulant theory, cumulants of ΔHm at each order can be calculated based on (25). Then the probability density functions of ΔHm can be obtained with Cornish-Fisher series. A cross-boundary vector , hβ,m can be obtained with elements calculated as

where, is the point of at which the probability of equals to βn . indicates that nth chance constraint is not satisfied at mth iteration. For nonchance constraints (or the chance constraints which are not considered in the model), the corresponding elements of hβ are considered as zeros.

3.3 Iterative learning control method

Iterative Learning Control (ILC) method can track the reference exactly in a finite time interval without the knowledge of the system models. (It applied to nonlinear systems.) It finds a suitable control action using the information from previous iterative results. A typical learning law of ILC is of the form:

where, L(·) is the learning function in ILC method. It is usually a function of the tracking errors in the learning process. Many control algorithms, ranging from simple scalar gains to sophisticated optimization computations, have been designed and applied in different areas [27-30]. Qing Wu, et al [30] used a PID-type iterative learning algorithm to detect, isolate, and estimate faults effectively in the fault diagnosis of a satellite. The PID-Type ILC method is also adopted in this paper to regulate the chance constraints. The discrete expression of a at (m+1) th iteration can be written as

According to (28), regulatory factor at (m+1) th iteration is updated using the information of (m-1) th iteration and m th iteration. The PID parameters are also updated adaptively by ILC to ensure good convergence behavior. The elements of regulatory factor corresponding to nonchance constraints and the satisfying chance constraints at 1st iteration are zeros. The reference of the model in the study keeps zero. It is a class of time-invariant nonlinear system with fixed target trajectory. Usually, the optimal solution can be obtained in several iterations by PID-based ILC.

The optimization framework for S-OPF with ILC is shown in Fig.1.

Fig. 1.The optimization framework for S-OPF

The memory block stores hβ,1,…, hβ,m and a1 ,…, am for later calculation. The probabilistic calculation is based on cumulant method and Cornish-Fisher series. The model (18)-(21) is solved by the classic interior point method. ILC method calculates the regulatory factors for the model (18)-(21). Usually, the chance constraints can be satisfied in finite iterative steps.

3.4 Optimization process

The optimization process is as follows:

1) Initialize: set a0 = 0 ;2) Solve the deterministic optimization model (18)-(21) with forecast inputs based on (22) and (23);3) Calculate the cumulants of ΔHm from (25) with outputs of step 2) and uncertain inputs;4) Generate PDFs of Hm from (24) based on cumulants and Cornish-Fisher series;5) If the iterations reach the preset maximum number, go to step7); else, go to next step;6) Calculate hβ,m from (26). If max(hβ,m) ≤ le − 4 , go to the end; else, update a according to (28) and go to step 2). The number of iterations (m) increases with one in value.7) Choose one of the selections: (i) go to the end, output the results of the maximum iteration and show “the optimal solution that satisfies the given probability constraints does not exist”; (ii) decrease β and go to step 2) to find a sub-optimal solution. The confidence levels ( β ) are determined by the reliability requirements of studied system or by compromised consideration between objective and reliability.

 

4. Case Study

The proposed optimization method is tested on a modified IEEE 14-bus system with wind farms. P-OPF is also considered and compared with the S-OPF in this paper. The P-OPF is solved by the steps 1) to 4) of the proposed optimization process. The solution of P-OPF equals to that of deterministic OPF.

The 14-bus test system is developed based on the standard IEEE 14-bus system and has two voltage levels: 110kV and 220kV. The connection and parameters of lines, voltage upper and lower bounds, and parameters of conventional generators are the same as those given in the standard IEEE 14-bus system [31]. Bus 1 is considered as the slack bus. The rated power is 100 MVA. The mean values of load powers equal to their nominal values and the parameters of the wind farms are from the actual wind farm data. In Table A.1 (in Appendix), the forecast values of loads are given in p.u.. The Std. (%) value is the ratio between the Std. of load forecast error and . The load of node 9 is discrete with distribution as shown in Table A.2. The forced outage rate of generator 1 is 0.09 and that of generator 2 is 0.08. Two wind farms are connected to bus 9 and bus 12 with forecasted power 0.45 and 0.30, respectively. The wind farm is modeled from a statistical point of view to analyze its impact on the dispatch of electrical network. The cumulants of forecast error of wind powers are shown in Table A.3. No. in Table A.3 is the number of the bus the wind farm connected to.

The chance constraints on (8) and (9) are considered. Therefore there are sixty-eight chance constraints. In order to explain the proposed optimization algorithm, the upper bound of apparent power in the branch between node 2 and node 3 is supposed to be 50MVA. When all elements of are set to 95%, there are six constraints which do not meet Eq. (13) using the proposed P-OPF. That is to say, the solution of deterministic OPF may be infeasible in stochastic environment. The six constraints can be represented as follows:

After iteration 10, the optimal solution can be achieved by the proposed optimization algorithm for S-OPF. Elements of a corresponding the six constraints are 0.0681, 0.00462, 0.00705, 0.00683, 0.00808, and 0.0712, respectively (other elements of are zeros). pn of the six constraints on each iteration regulated by the ILC method are shown in Fig. 2. It can be seen form Fig. 2 that C1 and C6 are the key ones which affect the ending condition.

Fig. 2.pn of C1-C6 regulated by ILC

Optimal results obtained with P-OPF and S-OPF are compared in Tables 1-3. According to the proposed optimization processes for P-OPF and S-OPF, the results of P-OPF are the solution of model (18)-(21) when a = 0 and that of S-OPF are the solution of model (18)-(21) with the final solution of a optimized by ILC. Table 1 shows the regulation of active and reactive power of the generators by the S-OPF compared with the P-OPF. The voltage magnitudes are shown in Table 2 and power flow of branch 3 in Table 3. It can be concluded from Tables 1-3: (i) The voltage magnitudes obtained with S-OPF have a certain margin to the limit. (ii) The system’s active loss is lowered using S-OPF. (iii) Adjusted by the proposed method for S-OPF, the apparent power of branch 3 are smaller than that of P-OPF. The final probabilities of C1-C6 solved by P-OPF and S-OPF are shown in Table 4. It can be seen from Table 4 that the probabilities of S-OPF are equal to or greater than the set confidence level. That is to say, the chance constraints are satisfied by the regulation of ILC in proposed algorithm for S-OPF. Table 4 shows the probability of satisfying C2 is more than that of C1. The reason is that the constraints C1 and C2 are correlated and they cannot reach the limits simultaneously due to the line loss on branch 3. With consideration of accuracy comparison, Monte Carlo Simulation (MCS) with 10000 samples is used as reference. Fig. 3 shows the PDFs of C1 solved by P-OPF and S-OPF using the proposed method and MCS, respectively. The PDFs of C2 -C6 change in the same way as C1.

Table 1.Output power of generators

Table 2.Bus voltage

Table 3.Power flow of branch 3

Table 4.Probabilities of the two methods

Fig. 3.Probability density curves of C1

Average Root Mean Square (ARMS) error is computed as

where, PMi is the ith point’s value on the cumulative distribution curve calculated using the proposed method. MCi is the ith point’s value on the cumulative distribution curve calculated using the Monte Carlo method; N represents the number of points and N = 30 here.

The relevant ARMS results of C1-C6 are shown in Table 5. It can be seen from Table 5 that the proposed methods for P-OPF and S-OPF can approximate Cumulative Distribution Function (CDF) of voltage magnitudes (C3-C6) more accurately than that of line flows (C1-C2).

Table 5.ARMS of C1-C6

Ignoring the stochastic factors, including load fluctuation, the failure rate of conventional generator and the forecast error of wind farms, the minimum cost of the deterministic calculation (or P-OPF) is $4063.250. When considering the stochastic factors mentioned above, the minimum cost of S-OPF with β = 0.95 solved by the proposed method is $4189.763 which is higher than that of P-OPF. Fig. 4 shows the objective function values of S-OPF when all the predefined probability levels increase from 0.8 to 0.98. The cost increases as β becomes greater. Although the generation cost of the S-OPF is greater than that of the P-OPF, the scheduling scheme of S-OPF can better withstand the error and fluctuation of the forecast loads and wind power.

Fig. 4.Objective values of S-OPF with different confidence levels.

Although the sensitivity term is difficult to be expressed in mathematical formula, it can be evaluated based on the experimental results. Take C1 and C6 as examples to analyze the responses of the objective value to different confidence levels of a single constraint. Set the confidence level of C1 or C6 as a specific value and other elements of β as 0.95, the results are shown in Fig. 5. The objective value increases as the confidence level of C1 or C6 varying from 0.80 to 0.95. The main increment of objective value occurs in [0.84 0.87] when the confidence level of C6 changes and in [0.90 0.95] when that of C1 changes. That is to say, sensitivities of objective value to constraints’ confidence levels have different sensitive areas and changing rules. The results of Fig. 4 and Fig. 5 are useful in determining suitable confidence levels by trade-off between objective and reliability.

Fig. 5.Objective values when a single confidence level varies.

All of the calculation is performed by a computer equipped with Intel(R) Core(TM) i5-2410M CPU@ 2.3GHz. The calculation time of P-OPF is 0.932 seconds and that of S-OPF is 9.128 seconds. The IEEE-118 case is also tested as the 14-bus, the calculation time of P-OPF and S-OPF are 2.125 seconds and 23.159 seconds, respectively.

 

5. Conclusion

Due to the complexity in solving the CCP-based S-OPF model, an indirect method is proposed based on limit regulation and ILC. Taking into account the relationship between the constraint equations and the optimization variables, the direct calculation of optimal solution is replaced by narrowing the solution space with minimum regulation through adding regulatory factors to the inequality constraints. The regulatory factors are controlled by ILC method to ensure the optimality of the solution. The final results verify the effectiveness of the proposed algorithm and the related solution validation process shows that the proposed model can help increase the reliability of power system. The results are also conducive to further consider the economic and security-constraint compromised dispatch and the risk management with the actual operational data.

 

Nomenclature

Parameters:

C2,i , C1,i , C0,i Coefficients of production cost functionCG The group of adjustable generators Minimum active output of generator i Maximum active output of generator iURi Ramp-up rate limit of generator iDRi Ramp-down rate limit of generator iUSRt Upper limit on system spinning reserveDSRt Lower limit on system spinning reserve Minimum reactive output of generator i Maximum reactive output of generator i Lower limit on voltage magnitude of bus i Upper limit on voltage magnitude of bus i Lower limit on apparent power of transmission line between bus i and bus j Upper limit on apparent power of transmission line between bus i and bus jkp , ki , kd PID parameters of ILC method

Variables:

PGi Active outputs of CGPWi Active outputs of ith wind farmQGi Reactive outputs of CGQWi Reactive outputs of ith wind farmSij,t Apparent power of transmission line between bus i and bus jPw(t) Actual wind power output Forecast wind powerNb Number of boundary-crossing constraints at 1st iteration Percentile of Hn corresponding to βnΔ Stochastic part of variablet Time periodm Counter of iterationsn Indexing number of the chance constraintsum mth input in the iterative process of ILC method

Vectors:

X Optimization variableθ = [θit ] Angle vector of bus voltagesV =[Vit ] Magnitude vector of bus voltagesa Regulatory factoramin Sub-vector of a corresponding to lower limitsamax Sub-vector of a corresponding to upper limitshmin Lower limit vector of inequality constraintshmax Upper limit vector of inequality constraintsβ Predefined confidence levelγ Perturbation parameterZ Positive slack variables when the inequality constraints are converted into equalityλ Lagrange multipliers of equality constraintsμ Lagrange multipliers of inequality constraintse A vector of all onesPd Forecast active inputsQd Forecast reactive inputsXm Solution of the deterministic OPF at m th iterationΔW Stochastic part of inputsρ A group of uncorrelated variablesr Tracking error vector during the iterative process in ILC method

Matrix:

Jm Jacobian matrix at m th iterationT Correlation matrix