Interior Point Algorithm Applied to the Optimization of the Power Supplied by a Wind Farm with a BESS

. Wind power constitutes a variable energy source that introduces unbalance in electrical network management because it cannot be programmed. Then, the possibility of storing wind energy becomes very important. The lack of control is a drawback that disappears when the combination of a wind farm (WF) and a battery energy storage system (BESS) is considered. In that case, the goal is to adjust the power plant output and the load requirements of electrical network, i.e., to contribute to system adequacy as much as possible. Considering the features of the problem, it can be defined as an optimization problem. Two algorithms are proposed to solve it: the primal dual algorithm and the Mehrotra predictor-corrector one. In both cases, the best solution of the proposed problem is reached in an efficient manner. The primal dual algorithm performs better in terms of time and the Mehrotra predictor-corrector one needs fewer iterations.


Introduction
The lack of control in electrical network management due to wind power grows with the installed wind power. Therefore, in order to improve the balance between generation and consumption in the electrical network, storage systems need to be installed altogether with wind farms. The main objective of having the capability to store wind power is to inject the power into the electrical network when needed, combining periods of decreasing wind farm production and increasing ones. The power surplus has to be stored somewhere and is given back when wind farm production is not enough to supply the required power.
Among others, a BESS is a profitable mean to store wind power due to its efficiency, adaptability and no need of specific installation conditions. When dealing with a system formed by a WF and a BESS, the first issue to consider is the battery size [1,2], which is critical. However, it is not dealt with here and the value used in the case study is based on a real power plant. The equations and constraints that define the behavior of the system were obtained from [3][4][5][6][7][8][9][10][11].
Another issue is the application time interval of the BESS, that can be in an hourly, daily, weekly or monthly term [12], but in this paper the objective is particularized to a daily term. Besides, the prediction is not considered here, it is assumed that the wind speed is known with enough anticipation to apply the algorithm, which is realistic due to the daily term analyzed.
It is also important to consider the model used for the BESS. Several models appear in the related literature [7,13,14] and they have been simplified here, considering certain efficiency value during the charging and discharging processes.
Finally, it is also important to decide the optimal allocation of the BESS in the network and its configuration [15][16][17]. In this paper, the BESS is located adjacently to the WF and the chosen configuration is dual inverter [17].
Then, in order to carry out a strategy based on the adaptation of the wind power produced by a power plant made up by a WF and a BESS, to the power requested along a period, an optimization problem arises. The objective should be that the power plant generates more power when it is more valuable. However, the power plant will have some technical constraints that limit its operating capacity, most of them related with the BESS but also with the whole installation, turning the optimization problem into a non-linear one.
Linear programming [18,19] is a part of mathematical optimization, which deals with the kind of problems where a linear objective function must be optimized, being subject to a set of linear constraints. Due to the specific circumstances of this type of problem, its optimal solution has to be located on the border of the set of feasible solutions. In those cases, the simplex algorithm is a good option to obtain the optimal solution. However, this algorithm goes over the border, trying to find the optimal one and, when dealing with a lot of variables and restrictions computation, the time consumption results very high. As an alternative to this algorithm there is the possibility of designing an algorithm that goes through the interior of the set of feasible solutions until the optimal one is reached. This algorithm is known as interior point algorithm [20].
In 1979, Leonid Kachiyan [21] designed an interior point algorithm, the ellipsoid algorithm, and proved that a linear programming problem can be solved efficiently in proper time by means of it. Later on, in 1984, Narendra Karmarkar introduced the interior point algorithm of projective transformation with great success [22]. In Karmarkar's work, primal-dual algorithm [23] and Mehrotra predictor-corrector one [24] are used to solve this type of optimization problem.
The organization of this paper is the following: section 2 includes the concepts, expressions and constraints to be considered in a hybrid WF-BESS system, section 3 explains the primal-dual algorithm, section 4 is an outline of the Mehrotra predictorcorrector algorithm, in section 5 the application of both algorithms to the problem proposed is given and the conclusions are stated in section 6.

Power Supplied by a Wind Farm with a BESS
A power plant consisting of a WF and a BESS is considered in this section. The equations and constraints that describe the behavior of this type of power plants are presented.
Generally speaking, the power plant output in a time interval is the sum of the power supplied by both the WF and the BESS in the same interval, as in (1).
where n is the number of intervals assessed, P is the power provided by the power plant in the interval i, P is the power provided by the BESS to the electrical network in the interval i and the power supplied by the WF in the interval i is named P . Notice that, in case that P is negative, it means that the total power provided to the electrical network is lower than the power supplied by the WF, but during this time interval, the BESS receives the rest of the power.
On the other hand, in eq. (2) the charging process of the BESS is expressed. The energy stored in the BESS in an interval is obtained as the sum of the energy previously stored in the BESS and the energy received during the interval. Notice that, the energy received is the product of the power provided by the BESS and the length of the interval with negative sign as expressed in (2). The reason for the negative sign is that P is defined as positive when supplied.
where E is the energy stored in the BESS in the interval i and T is the interval length. Therefore, the possibilities are: E > E during the charging process, E < E if the BESS is discharging or = in case there is no power from or to the BESS during the interval. There are losses during the charging and discharging processes and some kind of efficiency has to be considered. However, those losses are not taken into account here due to the complexity of the problem. Once the problem is solved, those losses can be assessed.
The constraints that must be considered to apply equations (1) and (2) are the following: i) BESS minimum and maximum charge levels (3): The reason for establishing a minimum charge level is to extend its useful life. The maximum charge level of the BESS is given by its maximum energy storage capacity. So, on each interval, constraints shown on (3) have to be taken into account.
where E is the minimum charge level and E is the maximum charge level of the BESS. ii) Maximum power transferred to or from the BESS (4): These limits are due to the heating processes that appear when charging or discharging the BESS. Therefore, there are constraints in both directions and the minimum one is negative.
where P is the maximum power that can be transferred to the BESS in absolute value. P is the maximum value of power that can be provided by the BESS. iii) Maximum and minimum level of power that can be supplied by the power plant (5): The reason for those levels is the installation configuration, because the limits depend on the whole power plant: it cannot consume power and the output power has to be between zero and its rated power.
where the rated power of the power plant is P .
Therefore, the optimization problem to be solved can be summarized as in (6), considering the relationships in (1) and (2).
where are the weights of each time interval i and define comparatively the need to supply power.
They are squared in order to adapt the solution to the characteristics of the problem.

Primal Dual Algorithm
In this section, a full explanation of the primal dual algorithm is given in order to let the reader know the process to obtain the solution of the proposed optimization problem.

Main idea
The primal-dual algorithm is based on the use of a logarithmic barrier. The idea is to replace the border of the feasible set by an element in the objective function that penalizes its value as it is approached. For example, it is supposed a standard linear programming problem, such as the one defined in (7): where f is the objective function, x is the vector of variables of the primal problem and Ax = b is the constraint of the minimization problem.
The equivalent problem given by (8) can be obtained.
where the parameter µ is called the barrier parameter.
Considering the standard format of a linear programming problem and its dual, (9), [25,26] min = = 0 min = where y is the vector of variables of the dual problem and s is the slack vector.
The following assumptions are made for the primaldual algorithm: i) The set of solutions of the primal problem, named S, is not empty.
ii) The set of solutions of the dual problem, named T, is not empty.
= { , , : iii) The range of A is maximum. Considering these constraints, it can be stated that both problems, primal and dual, have optimal coincident solutions. Besides, the feasible regions of both problems are bounded.
The logarithmic barrier function is applied to the sign constraints in both problems, following the process explained above as in (10).
A y + s = c s 0 y not constrained (10) As 0, it can be expected that the optimum solution of the problem P (Eq. 10-top) converges to the optimum solution of the original primal problem. In order to prove it, it is considered that the objective function is strictly convex and additionally using the duality theorems (Karush-Kuhn-Tucker conditions, KKT [18]), the system (11) is obtained.
where A and S are diagonal matrices with the values of the variables x and s, and e is a vector of ones.
From the assumptions i, ii and supposing that the feasible region of the primal problem is bounded, it can be derived that the problem P is feasible and it has just one minimum in x( ) for each > 0. Therefore, the system in (11) has just one solution.
The system in (11) also provides the necessary and sufficient conditions (KKT) for y( ) and s( ) to be the maximal of the problem D (Eq. 10-bottom). When the condition iii) is true, variable x determines uniquely the variable y from (11). Besides, the duality gap can be expressed as in (12). c x( )b y( ) = c -y( ) A x( ) = n (12) Therefore, as 0 the duality gap tends to zero and it implies that they are the optimal solutions of the primal and dual problem.
So the set of optimum solutions of the types of problem P and D form a path to the optimum solution of the original problem. Therefore, this algorithm starts from an initial point (x , y , s ) , after which a sequence of points is generated, (x , y , s ) S, T, choosing a movement direction and right step size, .

Movement direction and step size
This is one of the fundamental stages of the algorithm and Newton method is used on it. Newton method consists of a procedure to treat the solution of a system of nonlinear equations by means of the consecutive approximation of the system by linear equations. For example, dealing with a nonlinear function, F(z), defined from R to R , treating to find a vector z that nullifies the function. Using the multivariable series of Taylor, a linear approximation can be obtained, as in (13) ( where J is the Jacobian matrix of function F( ), with the rows formed by the partial derivatives of the respective functions of F( ) . Each element can be expressed as in (14).
As the intention is to assess the root F(z) = 0, the system in (15) is obtained.
The solution to this system permits to iterate from z to z with a unit step. This method converges quadratically when the starting point z is located close to the solution, but this good convergence is just a local behaviour. In a nonlinear generic case, if the starting point is located far from the solution, the Newton method can fluctuate permanently without reaching the solution When applying the Newton method to the system in (11), in order to move from a given point (x , y , s ) to another (x , y , s ), so that x and y are no negatives, the system in (16) is obtained, providing the Newton directions corresponding to each variable.
Developing the system (16), the relationships in (17) are obtained.

Ad = r A d + d = r
S d + X d = r (17) where: In case that points x, y, s are feasible, the rests r and r will be null.

Practical implementation
The practical version of the algorithm does not need that the initial point satisfies all the constraints. As long as the iterations grow, the sequence of points enters the feasible region.
Solving the system, the Newton direction for x shown in (19) is obtained.

d = D P D X e-D P D c + CA (AD A )r
where D = X S and P = I-D A (AD A ) -AD , that is the projection matrix in the kernel of matrix AD .
The elements in (19) have an interesting geometric interpretation: The first one is the centering direction, because it is the projection of the push vector (1/x ), and it helps keep the algorithm far from the border of the primal feasible region.
The second element is the direction of reduction of the objective function, because it is the projection of the negative gradient of the objective function of the primal problem, and it helps reduce quickly its value.
The third one is the feasibility direction, because r is a measurement of the feasibility of the problem.
It can be pointed out that the centering direction and the reduction direction of the objective function belong to the kernel of A. Therefore, the third direction is the only one that joins in the feasibility of the primal problem.
At the beginning, one of the most important directions is the feasibility one, because the initial point chosen does no need to belong to the feasibility region. As long as the iterations move along, these elements are losing its importance with respect to the other two.
Once the Newton directions are obtained when solving the system, the variables have to be updated in order to move to the next iteration. The values are obtained in eq. (20).
x = x + d y = y + d z = z + d (20) In order to obtain the coefficients that determine the step size used on each iteration, equations (21)    The direction of movement in an iteration (k) is determined by the value of the barrier parameter. In fact, this translation should be made several times for a fixed value of , so that the Newton steps converge to the central path corresponding to that . However, this could demand a large number of calculations to the algorithm, therefore, on each iteration the Newton steps are calculated just once and the value of the barrier parameter is reduced.
In order to reduce the barrier parameter (23) is used.
= (23) As stopping criteria for the algorithm, elements in (24) are used. norm(r ) norm(r ) (24) where the function norm( ) provides the absolute value of the vector. The first criteria in (24) ensures that the duality interval is small enough to verify that the optimum is close, the other two criteria ensure that the problem is feasible.

Mehrotra Predictor-Corrector Algorithm
As an alternative to the primal dual algorithm, the Mehrotra predictor-corrector algorithm is explained pursuing the same objective as the presented in previous section. This algorithm is a refined version of the primaldual algorithm. It consists of two stages, the predictor stage and the corrector one.
Applying the step of the Newton directions in the former algorithm to the system, the system in (25) is obtained.
e-X S e (25) reduction of the barrier parameter.
Here, the term d d has not been neglected. The general idea is to make a first estimation of d and d , and then to replace the product d d in (25) and to solve the system.
In this algorithm it is very important the Cholesky decomposition, because it can be used in both stages and it speeds up the resoluteness of the system.

Predictor stage
In this stage, the system in (26) is going to be solved.
The solution of the system provides the affined movement directions, used to estimate the product d d of the system in (25). Additionally, these values provide a specific , using the expression in (27).

= ( )
where the values of and are obtained using (28) and (29).
The parameters are calculated using the same criteria as in the primal-dual algorithm.

Corrector stage
Once the values of d and d are calculated, the system in (30) is solved in order to obtain the corrected movement directions.
where D means a diagonal matrix formed by the elements of d . The step sizes are obtained again after the final movement directions, and the points are updated to move to the next iteration. The stopping criteria are the same as in the primal-dual algorithm, i.e., those given in (24).

Algorithms assessment.
As explained above, two algorithms have been used to solve the proposed problem. Primal-Dual and Mehrotra Predictor-Corrector algorithms have been tested with the aim of comparing results for the same cases.
Initially, ten cases are considered consisting of randomly generating ten different wind power vectors, corresponding to a period of 24 hours. Besides, an interval of 10 minutes has been taken, resulting in a total amount of 144 variables.
The results of the time performance tests for both algorithms can be seen in table 1, where the number of iterations and the objective function have been included for the ten cases. As can be seen, no differences have been detected in terms of the optimal value of the objective function. In order to compare the elapsed time and number of iterations, a graphical comparison is given in figure 1, where the mean values of both algorithms are compared. In figure 1, it can be seen that the primal-dual algorithm performs better in terms of time and the predictor-corrector one does the same when referring to the number of iterations. Using this result as a reference, it can be inferred that the primal-dual algorithm is the best choice in order to solve the problem described in this paper. However, the predictor-corrector algorithm should not be discarded completely due to the low number of iterations, i.e., it may be of application when the operations involved in an iteration are higher than in the case of considering a period of one week.

Estimation of the optimal capacity of a wind farm with BESS
In this paragraph a way is shown to determine a range of optimal capacities of a BESS system in terms of network stability (power generation when needed). Considering the results of the previous analysis, the primal-dual algorithm is chosen to solve this specific problem.
In order to obtain the results, a specific case has been taken. The related information can be seen in the table 2.  The BESS system used in this paper will be the NAS battery system. The information of this system is provided in table 4. Grouping different modules, the desired power and capacity can be reached. The power (MW) has been used as an index of the capacity of the system, as it is easier to identify and operate within the calculations.
To carry on the research, a representative week per season has been considered and different values of capacities have been evaluated in order to analyze their behavior.
Ten different experiments per week have been done and for the following different power values: 4, 6, 10, 14,16,20,30,50 MW This results in a total of 320 experiments. After all, these average values of the indexes have been obtained.
With the aim of rating the good behavior of the batteries the relationship in (31) has been defined: This simply shows the relation between the cases with and without batteries. So the higher it is, the better the behavior of the power plant.
Vector c will be composed with the different weight parameters obtained from demand data. The representative weeks that have been taken can be seen below: Spring -13/04/2015 Summer -3/08/2015 Autumn -12/10/2015 Winter -19/01/2015 The information needed to model the wind resource was taken from the Wind Atlas of the IDAE. Table 5 shows the Weibull parameters used. All these input data gave the results provided in table 6. Here the positive parameter is represented by the ratio, while the negative one is the power, as it is considered as an expense. If the relation between this ratio and the installed power is plotted against the installed power, the graph presented in figure 2 is obtained.
The behavior is very similar to a decreasing exponential curve. So, as installed power continues being added the amount of improvement achieved is lowered.
Finally, an optimal range of power/capacity can be set between 20 and 30 MW, representing the 40 and 60% of the rated output power.

Conclusions
The primal-dual algorithm and the Mehrotra predictorcorrector algorithm have been successfully implemented to optimize the power supplied by a hybrid WF-BESS system. For a period of one day and intervals of 10 min, the primal-dual algorithm performs better in terms of time and, for the same conditions, the predictor-corrector performs better in terms of number of iterations. Therefore, it can be established that the first method is the proper one for the case described here while the second one can be of application for a higher number of intervals.
As a collateral result, using the primal-dual algorithm, it can be concluded that the rated power of a BESS should be around 40-60% of the rated power of the wind farm, because a higher rated power would be underused.