Random-Fuzzy Chance-Constrained Programming Optimal Power Flow of Wind Integrated Power Considering Voltage Stability

Considering the random and fuzzy nature of wind speed, this paper proposes a multi-objective random-fuzzy chance-constrained programming optimal power flow (OPF) for wind integrated power systems. The proposed method is based on random-fuzzy chance-constrained programming. The optimization model aims at minimizing generation cost, carbon dioxide (CO2) emission, and system functional power loss, and P-Q-V steady-state voltage stability is included in the constraints. Based on random-fuzzy chance-constrained programming, the corresponding solution process of the proposed multi-objective OPF is proposed, which is a hybrid of random-fuzzy simulation, non-dominated sorting genetic algorithm-II (NSGA-II), and fuzzy satisfaction-maximizing decision-making method. The proposed approach is simulated on the IEEE 30-bus system to provide an example of its application. The simulation results demonstrate that the proposed random-fuzzy chance-constrained programming OPF has higher security and more economy than dynamic stochastic optimal power flow (DSOPF) and dynamic fuzzy optimal power flow (DFOPF).


I. INTRODUCTION
Wind power generation (WG) is currently a popular form of renewable wind energy. Since its power output is fluctuating and uncertain, high WG penetration into a power system brings the power scheduling challenges [1].
Traditionally, research studies on WG output mostly treat its uncertainty either as randomness or fuzziness. A model of WG output incorporating its random nature can be built generally in two ways. One is to depict the statistically uncertain characteristics by a probability distribution, usually representing wind speed by Weibull distribution in simulation followed by producing WG output based on the relationship between wind speed and power output [2]. The other is to utilize the random time sequence methods such as Markov chain [3] and autoregressive integrated moving The associate editor coordinating the review of this manuscript and approving it for publication was Ton Duc Do . average model [4]. Approaches to building a fuzzy model of WG output can also be classified into two kinds: direct data mining of historical statistics, or transferring WG's fuzziness into the fuzziness of forecast error [5].
However, for the uncertain wind speed, randomness and fuzziness both exist. Although wind speed is generally accepted to obey Weibull distribution, the parameters of the Weibull distribution are fuzzy and typically limited by historical statistics. In considering such phenomenon of a combination of randomness and fuzziness, wind speed is more accurately and rationally defined as a random-fuzzy variable [6] based on the fundamental theory introduced in [7], and the corresponding modelling and random-fuzzy simulation methods are utilized. There also exist other ways to consider randomness and fuzziness together. In [8], the electric load is described as an uncertain variable which is divided into the random part and fuzzy part, and random-fuzzy neural network is applied to tackle the uncertainties of VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ load forecasting. In [9], fuzzy-probabilistic modelling techniques and a hybrid method of fuzzy set and Monte Carlo simulation for power system risk assessment is proposed to capture both randomness and fuzziness of loads and component outage parameters. A random fuzzy model proposed in [10] is to accommodate the impact of various operation conditions and other factors on the failure probability of system components. Faced with the uncertainty of power system, probabilistic optimal power flow (POPF) [11], [12] and stochastic optimal power flow (SOPF) [13] are proposed to analyze and optimize generation power schedule, as well as distributionally robust chance-constrained OPF in [14]- [16]. For situations where the uncertainty distribution is not fully known, reference [14]solve the problem with distributionally robust optimization, to ensure that any distribution in the established fuzzy set satisfies the chance-constrained. SOPF based on chance-constrained programming [13], [17] can obtain optimal dispatch schedule, which meets with a certain confidence level of chance-constraint. Chance-constraint in SOPF is handled by probabilistic power flow (PPF). The uncertainty quantification approaches proposed for PPF can be technically divided into Monte Carlo simulation method, analytic method and approximate method [18]. As for the Monte Carlo simulation method, probabilities when state variables satisfy the constraint are calculated after obtaining repeated simulation samplings and stochastic power flow distribution of the system. Currently, research on SOPF mostly concentrates on a single time section. In this direction, reference [19] from the aspect of day-ahead time horizon and depicting WG output forecast via auto regression moving average, has researched chance-constrained DSOPF which considers randomness of forecast error and temporal-spatial correlation. In [20], the scenario tree model is applied to approximate the stochastic nature of WG in the 24-h operation horizon optimization problem. In [21], proposes a improved black hole (IBH) algorithm, to achieve practically coordinate the mobile PEVs in multiobjective constrained constrained dynamic optimal power flow (DOPF)problem. However, most of the above studies only focus on randomness or fuzziness to study OPF of system. Once the uncertain power injection of WG is depicted by the random-fuzzy variable, traditional stochastic chance-constrained programming should be extended to random-fuzzy chance-constrained programming.
When the power system is under uncertain condition, voltage stability plays an essential role in constraining the operation. In [13], SOPF is extended to include constraints for voltage stability as well as small-signal stability, based on approximating the voltage stability and small-signal stability constraint surfaces with second-order approximations in parameter space. In [22], voltage stability constrained optimal power flow model is presented, the text pointed out that the stability index can be incorporated via being added as a new voltage stability constraint in the OPF constraints or being used as a voltage stability objective function. These studies are mainly for post contingency condition and on single time section. Besides, high renewable energy penetrations in power systems impact the investigation of steady-state voltage stability [23]. In [24], a new index is proposed called area of the voltage stability region (AVSR) assess steady-state voltage stability margins via P-Q-V curve, which is superior in reflecting the boundary of the maximum real/reactive power and the voltage collapse magnitude in three-dimensional space. OPF considering P-Q-V steady-state voltage stability is seldom reported.
Where there are multiple conflictive objectives to be taken into account in OPF formulation, researchers have carried out many studies on solving multi-objective OPF [25], such as the multi-objective differential evolutionary algorithm based decomposition (MOEA/D) in [26] and non-dominated sorting genetic algorithm-II (NSGA-II) in [27]. In [28], the proposed dynamic multi-objective unit commitment problem is solved by NSGA-II to obtain Pareto optimal solutions, and fuzzy satisfaction-maximizing method is adopted in decisionmaking.
This paper considers the random and fuzzy nature of wind speed. It proposes a random-fuzzy chance-constrained programming OPF of wind integrated power system based on random-fuzzy constraint programming. The significant contributions are: 1) Current OPF study does not consider random-fuzzy nature of wind speed. Therefore, in the OPF study of this paper, to cover more information of wind uncertainty, and solve the problem of errors in the digital features of random variables due to insufficient historical data, WG power output is depicted by a random-fuzzy model and simulated by random-fuzzy simulation method; 2) Moreover, existed multi-objective OPF study does not incorporate P-Q-V steady-state voltage stability. So a multi-objective OPF model which incorporates P-Q-V steady-state voltage stability constraint is established based on random-fuzzy chance constraint; 3) Since a new multi-objective OPF model is proposed, traditional solution process is not available to solve the model. Based on random-fuzzy simulation and random-fuzzy chance-constrained programming, the progress is a hybrid of random-fuzzy simulation, non-dominated sorting genetic algorithm-II (NSGA-II) and fuzzy satisfaction-maximizing decision-making method.

II. RANDOM-FUZZY CHANCE-CONSTRAINT A. Random-FUZZY MODEL OF WG POWER OUTPUT
Based on the previous studies, it is reasonable to assume that wind speed v follows Weibull distribution, with the probability density function [12], where c is scale parameter and k is shape parameter. However, limited by historical stochastics, it is hard to obtain the precise probability density function, i.e. the parameters of the density function are of fuzzy nature. Through statistics and data mining, c and k can be depicted by some specific fuzzy variables respectively as ξ c and ξ k [6]. Taking trapezoid fuzzy variable ξ x = (r 1 , r 2 , r 3 , r 4 ) as an example (where r 1 < r 2 < r 3 < r 4 ), the membership function is When r 2 = r 3 , the trapezoid fuzzy variable is simplified to a triangle fuzzy variable represented by ξ x = (r 1 , r 2 , r 4 ).
With parameters ξ c and ξ k , wind speed is depicted by random-fuzzy variable ξ v [7], and the chance measure distribution function is [6] The active power output of a large scale wind farm can be modeled approximately as the aggregation of single wind turbines.
where v ci , v co and v r are respectively cut-in, cut-out and rated wind speed; P WGr is rated power output of each WG unit; N WG is the number of wind turbines. Considering the wind farm is regulated under constant power factor mode, working at leading phase ϕ WG , its reactive power injection is

B. RANDOM-FUZZY SIMULATION SAMPLING AND THE CHANCE CONSTRAINT ON SYSTEM STATE VARIABLES
Generally, chance constraint in SOPF is handled by PPF. The uncertainty quantification approaches proposed for PPF can be technically divided into Monte Carlo simulation, analytic and approximate methods [18]. As for analytic and approximate methods, the chance constraint can be converted to a certain constraint and then the optimization problem can be solved. However, the random-fuzzy variable is with double uncertainty, and its chance measure is a transformation from possibility space to probability space which can be expressed as a mathematical function from the interval (0,1] to [0,1]. Further considering the fuzzy variable with its membership and the random variable with its probability distribution, the complex relationship between input variables and state variables of the power system, it is hard to apply the analytic method.
In Monte Carlo simulation method of traditional SOPF, after Monte Carlo simulation sampling of the input variable, state variable samples and its statistical frequency are obtained for chance-constraint. Thousands of simulation can ensure the reliability of the results to a certain extent.
Inspired by these considerations, in this paper, we assume that wind speed follows those above random-fuzzy Weibull distribution, and daily active/reactive power output sequence samples of WG are obtained through random-fuzzy simulations. Then the Newton-Raphson power flow method brings out the random-fuzzy power flow samples of the power system, finally by calculating the frequency of satisfying the node constraint conditions in the power flow result samples by periods. This way, traditional stochastic chance-constraint is extended to random-fuzzy chance-constraint.
The detailed process is described as follows: 1) Based on the membership functions, extract a combination of k-c which satisfies a particular membership level λ; 2) Examine if k < c. If so, this k-c combination is taken as the Weibull parameters during this simulation time. If not, the process returns to step 1); 3) Inverse transformation of (3) Daily power output sequences of the WG are obtained by (4) and (5); 5) For power injection of WG in every time, period Newton-Raphson method is applied to calculate power flow, and power flow results such as bus voltage, generation reactive power output are obtained; 6) Steps 1) ∼ 5) proceed for N S times, and power flow samples of different periods are obtained; 7) The frequency that the samples satisfy the constraints on bus voltage k U,t , generation reactive power output and AVSR k AVSR,t of each time period are counted. Then the chance measurements are calculated by

III. P-Q-V STEADY-STATE VOLTAGE STABILITY CONSTRAIN
Steady-state voltage stability margin can be more comprehensively and accurately assessed using P-Q-V curve, which shows the boundary of the maximum real/reactive power and the minimum voltage magnitude in three-dimensional space. P-Q-V curve is obtained via continuation power flow method, with varied load increase direction indicated by different power factors. Then, for a certain operating point P i , Q i of VOLUME 8, 2020 a specific PQ bus indexed by i (i = 1, 2, · · · n PQ , n PQ is the number of PQ buses), the voltage stability margin can be estimated by calculating AVSR i where (P cr i , Q cr i ) is the voltage collapse point of d n direction; A i , B i ,and C i are parameters of voltage stability boundary obtained by curve fitting [17].
A bus's voltage steady-state stability can be ensured by its AVSR value larger than zero. Allowing some margin, it can be where P cr i 0 (A i P 2 + B i P + C i )dP is the maximum margin, i.e. when the bus load is zero; η is the AVSR spare percentage.

A. OBJECTIVES
The three objectives are the minimization of generation cost, CO 2 emission and active power loss of system.
where a i , b i , c i are cost coefficients for coal-fired generation i; α i , β i , γ i are emission coefficients;Ũ i,t is the voltage of bus i (i = 1, 2, · · · , n bus , n bus is the number of system buses) in period t;δ ij,t is the phase difference between bus i an j; G ij,t and B ij,t are the conductance and susceptance of branch between bus i and j; is cluster of all buses connected with bus i. The superscript ''∼'' in following content denotes uncertain variable. Note that among all bus voltages of system, voltages of PV buses are certain variables.

B. CONSTRAINTS 1) STATIC CONSTRAINTS
For ∀t = 1, 2, · · · T , the following static constraints need to be satisfied: power balance, generation power output, bus voltage, system spinning reserve and AVSR. Limits on control variables are expressed as certain inequalities, while those for state variables are random-fuzzy chance constraints.
where P Li,t and Q Li,t are active and reactive load on bus i; andQ Gi,t denote reactive power output of generation on bus i.
where µ is the system spinning reserve percentage.
c: CONSTRAINTS OF CONTROL VARIABLES P min Gi ≤P Gi,t ≤ P max Gi , ∀i = 1, 2, · · · n G (4.7) U min Gi ≤Ū Gi,t ≤ U max Gi , ∀i = 1, 2, · · · n G (4.8) where P min Gi and P max Gi are active power output limits of generation i; U min Gi and U max Ch Q min Gi ≤Q Gi,t ≤ Q max Gi (λ) ≥ ρ Q , ∀i = 1, 2, · · · n G (4.11) 217960 VOLUME 8, 2020 Ch P cr where U min i and U max i (i = 1, 2, · · · n PQ ) are the voltage limits of PQ bus i; Q min Gi and Q max Gi are reactive power output limits of generation i; λ and ρ are confidence levels.

2) DYNAMIC CONSTRAINTS
The active power ramp rate constraint of coal-fired generations is − P down Gi ≤P Gi,t+1 −P Gi,t ≤ P up Gi , ∀i = 1, 2, · · · n G , t = 1, 2, · · · , T (4.13) where P down Gi and P up Gi are the maximum ramp rate and descending rate of generation i.

V. SOLUTION METHOD OF MULTI-OBJECTIVE RANDOM-FUZZY CHANCE-CONSTRAINED PROGRAMMING OPF A. SOLUTION PROCESS
The solution process for the proposed optimization is shown in Fig. 1.
First, uncertainty of WG power output is neglected, and the random-fuzzy expectation is taken as power injection. By solving certain multi-objective DOPF via NSGA-II and fuzzy satisfaction-maximizing decision-making method, preliminary power schedule is obtained.
Then, under this scheme, WG power output is depicted by random-fuzzy uncertain model. Times of random-fuzzy simulations combining with Newton-Raphson power flow calculations bring out random-fuzzy chance measure distributions of each state variable in the system. Subsequently, random-fuzzy chances are evaluated based on the confidence level constraints. If confidence level is satisfied, the process ends and the solution of multi-objective random-fuzzy chance-constrained programming OPF and its corresponding objective values and power flow are output. If not, the upper and lower limits of chance constraints should be adjusted, and process restarted from the very beginning.
Iteration process would come to an end when the scheme which satisfies all chance constraint is found.

B. PENALTY FUNCTION METHOD OF STATE VARIABLE CONSTRAINT IN MULTI-OBJECTIVE OPF
While solving the multi-objective DOPF, constraints on state variables are achieved though dynamic penalty function method. Note that the active power output of coal-fired generation on the reference bus is regarded as a state variable during the optimization.
It is examined whether PV bus reactive power, PQ bus voltage, reference bus power and its generation ramp-rate, and AVSR violate their limits. If so, the corresponding penalty function of this section is set 0. Otherwise, it is calculated according to the out-of-limit value. Total penalty function is In NSGA-II optimization, non-dominated sort and crowding distant assignment proceed based on the calculation of synthesized objective cf m (u t ), where f m (u t ) denotes the value of objective m (m = 1, 2, · · · M , M is the number of objectives); k is the current VOLUME 8, 2020 iteration times in NSGA-II; h(k) = k √ k is the dynamic adjustment factor.

C. FUZZY SATISFACTION-MAXIMIZING DECISION-MAKING METHOD
The proposed model with three objectives is solved by a hybrid method of NSGA-II and fuzzy satisfaction-maximizing method, where NSGA-II is adopted as the solver and fuzzy satisfaction-maximizing method as the chooser of all feasible solutions. Rather than incorporating the objectives together into one objective by ways of simple sum or weigh factor, NSGA-II can obtain Pareto optimal frontier by applying non-dominated sorting. Note that, Pareto frontier can offer a number of optimal solutions for decision makers to select according to their preference. Here the fuzzy satisfaction-maximizing method is adopted only as an example to make decision.
For a single non-dominated solution indexed by n (n = 1, 2, · · · , N , N represents number of the non-dominated solutions), satisfaction of each objective value is calculated by where f m min and f m max are respectively the maximum and minimum value of objective m in the Pareto solution set. Then the total fuzzy satisfaction µ n of the nth nondominated solution in Pareto frontier is calculated by (5.4) and the one with the maximum value is selected.
In this way, the optimization of the objective is represented by fuzzy satisfaction. When making decision, it is fuzzy satisfaction that makes sense, not the objective value itself. So the units as well as magnitude order of the objectives are not influenced with each other.

D. METHOD TO DEAL WITH RANDOM-FUZZY CHANCE CONSTRAINTS
The random-fuzzy constraints shown as (4.11)-(4.13) can be expressed in a standard form PQ bus voltageŨ i,t is taken as an example. During the uncertainty process, if its random-fuzzy chance measure does not satisfy the confidence level, i.e. the chance thatŨ i,t ∈ [U min i , U max i ] is not greater than the set confidence level, there are two situations.Ũ i,t samples mainly fall in the lower part or the higher scope. Hardly occurred is the situation that samples scatter in the space lower than U min i and higher than U max i . So, whenx i,t is unsatisfied with chance constraint, i.e. Ch x min i ≤x i,t ≤ x max i (λ) < ρ x , the average of the samples fall out of [x min i , x max i ] which is represented byx av i,t is calculated.
1) Ifx av i,t ≤ x min i , the upper limit of this state variable in this time period remains unchanged and the lower one is adjust as are the current limits. The updatedlimits [x min ] are applied in the next certain DOPF process; 2) Ifx av i,t ≥ x max i , the lower limit of this state variable in this time period remains unchanged and the upper one is adjusted as The updated limits are applied in the next DOPF process.

VI. CASE STUDY
An illustrative case is studied based on IEEE 6-units \30-nodes system of MATPOWER. And set relevant parameters according to IEEE 30-nodes system of MATPOWER. The optimization is for a day split into 24 time periods. Load curve is shown in FIGURE 2. In the simulation, to reduce the impact of wind power load access on the transmission system, this paper chooses to connect the WG units to bus eight, which has a more massive load. Parameters of coal-fired and WG units are listed in TABLE 1 and TABLE 2,  respectively. W0odel is according to data mining of historical wind statistics from National Renewable Energy Laboratory in America. It is estimated that the Weibull distribution parameters of the August wind speeds are adequately represented by triangle fuzzy variable ξ k = (1.14, 1.75, 3.64) and trapezoid fuzzy variable ξ c = (2.95, 4.40, 6.40, 8.22) [6].
The random-fuzzy simulation times are set 500, each with 24 data samples. As for random-fuzzy chance constraint, λ = 0.9 and ρ U = ρ Q = ρ AVSR = 0.9. P-Q-V curve for bus eight is obtained as FIGURE 3. shows. The blue lines denote the variation of voltage magnitude against real power P and reactive power Q in each dot,; the red line is the P-Q-V voltage stability boundary [23]. AVSR spare is set 10%. Within NSGA-II algorithm, population size = 100, Pareto fraction = 0.3 and generations = 500.
During optimization, Pareto frontiers were obtained as the example shown in FIGURE 4. Within the figure, each dot denotes a feasible optimal schedule, so the Pareto optimal frontier can offer decision-makers options to choose from. In this paper, the fuzzy satisfaction-maximizing method is adopted to make decision.
Optimization finished in 8 iterations, and optimal power scheduling is shown in FIGURE 5. FIGURE 5, comprises     three sub-figures, which represent the active power output of coal-fired generation, the generation voltage and the shunt susceptance of shunt capacitor respectively. Random-fuzzy WG power output simulation samples set is shown in FIGURE 6.
When the iteration process is finished, the updated limits of state variables are obtained. By contrast with the original constraints limit, it is easy to get adjustments. Take voltage of bus 12 for example, the lower limits of voltage constraints in almost each time period have been adjusted, as shown in Table 3. It shows that with the random-fuzzy power injection of WG into the bus 12 is more vulnerable and the bus voltage would easily rise. Under the DOPF schemes in the first and last iteration respectively and with random-fuzzy power injection samples of WG, the corresponding bus voltage sample distributions of bus 12 are as FIGURE 7. shows (take 6th time period as an example). Apparently, by narrowing the limits of constraints, the random-fuzzy samples  finally fell mainly in the secure region, and the random-fuzzy chance constraints were satisfied. Random-fuzzy chance is 0.274 in the first iteration while 1 in the last one. Constraints on the reactive power of coal-fired units and AVSR were seldom adjusted, as the AVSR samples of 6th period shown in FIGURE 8. It indicated that the reactive power regulation capability of the coal-fired units is abundant, so is the load margin.
To test the advantages of the method proposed in this article, compare it with the DSOPF and DFOPF algorithms. For DSOPF case, we use Monte Carlo simulation processing, here we take k = 2.07 and c = 5.49. For DFOPF case, we use fuzzy processing on wind speed, set the confidence level to 0.9. Similarly, set the simulation time to 500 × 24.
The maximum satisfaction compromise is listed in Table 4. From the table, compared to DSOPF, the algorithm proposed in this paper has less cost, power loss and iteration  times, which has better economy. And compared to DFOPF, the CO2 emissions, power loss and iteration times of the algorithm proposed in this article have been improved, but the cost is less. And the calculated standard deviations of the random-fuzzy samples and fuzzy simulation samples are 5.1912 and 3.8762 respectively. It can be seen that the proposed method can describe richer uncertainties.
The severe fluctuation of wind power injection causes undulation of the power system's state variables, and the chance constraint is hard to be ensured. So the optimization problem is hard to be solved, which results in more iteration time in random-fuzzy chance-constrained programming OPF than in DFOPF. However, when testing the optimal scheme obtained in DSOPF and DFOPF with random-fuzzy WG power output samples, a few state variables breach the chance constraints. This phenomenon indicates that, when applying the more rational random-fuzzy model to depict and simulate WG power injection, the optimal schemes of DSOPF and DFOPF are incapable of meeting security requirement.
So the proposed random-fuzzy chance-constrained programming OPF is of necessity and significance. Besides, according to the optimization results, it can be known that the proposed method can obtain a safer and more economical solution, and can describe richer uncertainties.

VII. CONCLUSION
This paper has studied multi-objective OPF of wind integrated power system based on random-fuzzy chanceconstrained programming. Major conclusions are as follows: (1) Modeling wind speed to follow a probability distribution with its parameters represented as fuzzy, this paper has proposed to use the random-fuzzy variable model to depict wind speed and adopt random-fuzzy simulation to produce WG power output samples. By comparing with traditional Weibull distribution stochastic model and Monte Carlo simulation method, it was shown to be able to cover more information about wind uncertainty and uncertain scenarios of WG power output. So it is a more comprehensive methodology.
(2) Faced with a random-fuzzy power output of WG, multi-objective OPF model was established based on the random-fuzzy chance constraint. The model included three objectives, and P-Q-V steady-state voltage stability was incorporated into regulations.
(3) This paper has provided the solution process to the proposed optimization model. By adopting NSGA-II, Pareto optimal solution sets were obtained. By the fuzzy satisfaction-maximizing method, individual multi-objective DOPF solution was selected as the preliminary dispatch schedule. Then random-fuzzy power flow is calculated via combining random-fuzzy simulation and Newton-Raphson strategy, and chance measure distributions of state variables are obtained. According to this, the upper and lower limits of chance constraints are adjusted, and multi-objective DOPF is solved again. Iterations bring out a multi-objective optimal schedule which satisfies all random-fuzzy chance constraints. Optimization result shows that the optimization model and solution method is available.
(4) An illustrative case is presented based on the IEEE 30-bus system to verify the feasibility of the proposed method. The simulation results demonstrate that the proposed method can obtain a safer and more economical solution, and can describe richer uncertainties.
Add system parameters and constraints, such as power market, operating costs, load changes, etc. To meet the actual needs of energy construction, consider extending the model in the article to a larger node system.