Towards Smart Energy Grids: A Box-Constrained Nonlinear Underdetermined Model for Power System Observability Using Recursive Quadratic Programming

: This paper introduces an underdetermined nonlinear programming model where the equality constraints are fewer than the design variables deﬁned on a compact set for the solution of the optimal Phasor Measurement Unit (PMU) placement. The minimization model is e ﬃ ciently solved by a recursive quadratic programming (RQP) method. The focus of this work is on applying an RQP to attempt to ﬁnd guaranteed global minima. The proposed minimization model is conducted on IEEE systems. For all simulation runs, the RQP converges superlinearly towards optimality in a ﬁnite number of iterations without to be rejected the full step-length. The simulation results indicate that the RQP ﬁnds out the minimal number and the optimal locations of PMUs to make the power system wholly observable.


Introduction
Phasor measurement unit (PMU) is a metering device that can provide real-time voltage and current synchrophasor measurements with high accuracy. PMUs play an imperative role in state estimation, monitoring, protection and wide area control in power systems [1]. Due to cost reasons, a full deployment of PMUs in the network is not realistic. One of the most important addressed issues is the strategic choice of the minimum number and locations of PMUs, ensuring complete network observability. Therefore, the optimal PMU placement (OPP) problem should be solved in order to make the power system completely observable by optimally placed PMUs at network buses [2,3].
Deterministic and stochastic algorithms are implemented for the solution of the OPP problem [2,3]. Several deterministic algorithms have been published in the literature for the solution of the OPP problem. The dominant optimization method regarding the OPP problem is the Binary Integer Linear Programming (BILP) technique. The BILP model requires the minimization of the objective function with inequality constraints and binary-valued decision variables [4][5][6][7][8][9][10].
A Weighed Least Square (WLS) algorithm is developed in [11] for the OPP problem. Authors of [12] proposed a nonlinear programming model for the solution of the OPP problem. In [13][14][15] The proposed nonlinear model is solved using a Recursive Quadratic Programming (RQP) method with super-linear convergence properties avoiding the Maratos effect.

•
The innovation of the local search procedure is that the RQP converges super-linearly towards optimality, satisfying the binary restriction.

•
The RQP presents a fast convergence rate towards optimality.

•
The RQP method delivers multiple optimal solutions in a reasonable time with those consumed by a BILP model being solved by the branch-and-bound method (BBM).
The remainder of this paper is organized as follows. Section 2 presents the basics of the Recursive Quadratic Programming (RQP) method. Section 3 formulates the OPP problem as an RQP optimization model. Section 4 gives an overview of computational results, whereas a performance evaluation of the proposed RQP algorithm for the solution of the OPP problem is carried out in Section 5. Finally, Section 6 concludes the paper.

A Recursive Quadratic Programming Background
The Recursive Quadratic Programming Algorithm solves nonlinear optimization problem of the form [35]: min x∈ n J(x) The RQP methods belong to the most frequently used algorithms for the solution of practical optimization problems due to their robustness and their good convergence properties. RQP methods can be proven to achieve global convergence and locally superlinear convergence rate; even locally quadratic convergence can be attained if the Hessian of the Lagrangian ∇ 2 xx L(x, λ) is available. Note that, the Lagrangian is given by: where λ ∈ m is the vector of Lagrange multipliers. The fundamental principle of RQP methods is to find Karush-Kuhn-Tucker (KKT) points, i.e., points that satisfy the necessary optimality conditions [35]. RQP methods mainly apply Newton's method to find zeros of that satisfy the constraints G, using a local linearization.
In this basic form, quadratic subproblems of the form are solved to obtain a search direction d.
Often, the Hessian of the Lagrangian ∇ 2 xx L(x, λ) is replaced by an update of Broydon-Fletcher -Goldfarbo-Shanno (BFGS) type, which has the additional benefit that only strictly convex quadratic programs have to be solved [35]. The computation of derivatives is a crucial element in nonlinear optimization. Mainly the first derivatives, i.e., the gradient of the objective function and the Jacobian of the constraints, is necessary to find a descent direction. The RQP is an iterative process that defines a series of points x [k] (called iterates) that (ideally) converge to an optimal point. It consists very roughly of the following steps [35]:

1.
Check termination criteria: For testing an iterate x [k] , the first-order optimality conditions KKT have to be evaluated.

3.
Use the solution from 2, to define a new iterate employing a merit function to find a suitable step-length.
Merit functions are needed to obtain a scalar measure of goodness of the trial point that is introduced to achieve global convergence (since RQP is derived from Newton's method that also depends on line search for globalization). The merit function is given by [35]: After determining the search direction d [k] from the QP, we have to find a suitable step-length a [k] . The line search procedure can be used to determine a step-length parameter a [k] based on the Armijo rule. The step-length is chosen to yield a sufficient decrease of a suitable merit function, which measures progress towards optimality [35]. Figure 1 shows RQP's flowchart. A difficulty in the context of RQP methods is that the line search procedure may prevent the full step from being taken [37]. The step-length is chosen to yield a sufficient decrease of a suitable penalty function, which measures progress towards optimality [37]. The convergence rate depends upon the reduction on the penalty function at each iteration, which relies on the choices of the descent directions and the step-length a [37]. The penalty function may not allow a unity step-length near the solution and thus may prevent q-superlinear convergence [37]. This phenomenon is called as the Maratos effect [37]. By using the RQP method through the fmincon of MATLAB optimization toolbox [38], the unity step-length is being taken decreasing the nonsmooth exact penalty function whenever the current iterate is sufficiently close to a local minimum point. Thus, the Maratos effect can be obviated, ensuring superlinear convergence to the optimum point. The RQP is converged when the difference between the current objective value and the previous objective value is less than the optimality tolerance [35].
Additionally, the constraints are satisfied within the feasibility tolerance, whereas step-tolerance specifies the termination tolerance for the design variable [36]. Therefore, the termination is succeeded based on the tolerance criteria in terms of objective evaluation, maximum constraint violation and first-order optimality [36]. A difficulty in the context of RQP methods is that the line search procedure may prevent the full step from being taken [37]. The step-length a ∈ (0, 1] is chosen to yield a sufficient decrease of a suitable penalty function, which measures progress towards optimality [37]. The convergence rate depends upon the reduction on the penalty function at each iteration, which relies on the choices of the descent directions and the step-length a [37]. The penalty function may not allow a unity step-length near the solution and thus may prevent q-superlinear convergence [37].

Optimum Design Model Formulation
This phenomenon is called as the Maratos effect [37]. By using the RQP method through the fmincon of MATLAB optimization toolbox [38], the unity step-length is being taken decreasing the nonsmooth exact penalty function whenever the current iterate is sufficiently close to a local minimum point. Thus, the Maratos effect can be obviated, ensuring superlinear convergence to the optimum point. The RQP is converged when the difference between the current objective value and the previous objective value is less than the optimality tolerance [35].
Additionally, the constraints are satisfied within the feasibility tolerance, whereas step-tolerance specifies the termination tolerance for the design variable [36]. Therefore, the termination is succeeded based on the tolerance criteria in terms of objective evaluation, maximum constraint violation and first-order optimality [36].

Optimum Design Model Formulation
Let the continuous decision variables x i denote the presence (x i = 1) or absence (x i = 0) of PMU at bus i. The OPP problem is formulated as a nonlinear minimization model [12]: where, x = (x 1 , · · · , x n ) T is a decision variable vector, whose entry x i is set to either 1 or 0 to denote the presence or absence of a PMU at bus i, W = diag(w 1 , · · · , w n ) is the PMU cost or weight matrix, n is the number of buses,0,1 is a vector whose entries are all zeros and ones, respectively, and f (x) is a vector function whose ith entry defines an observability constraint for the ith bus [11,12]: where N is the set of all buses and N i is the set of buses adjacent to bus i. Each Equality Constraint (3) implies that at least one PMU should be installed at any one of the buses i and j ∈ N i in order to make bus i observable. The optimal values of decision variables x i will be 1 or 0 as it has been proved in [12]. The above minimization model is reformulated and extended to employ with the underdetermined case.
Since we want to solve problems with a large number of constraints, many of them maybe are redundant [33,34]. The constraint redundancy gives rise to Jacobian rank deficiency, i.e., a singularity in iterative estimates of the constraint-set [34]. The preprocessing phase aims to simplify a given optimization problem by detecting and removing redundant constraints [33,34]. The preprocessing phase is applied to the primary constraint function to remove redundancies from the given constraints [33,34]. A redundant constraint is not required to define the boundaries of the feasible set, so the solution space is not affected. When the redundant constraints are dropped, the feasible region is expanded and more feasible solutions are located, whereas more constraints shrink the feasible region [35]. Since the enlarged feasible region involves more feasible points than the feasible region of the original model, more local minima are detected.
If a pre-solve procedure is applied after the model has been built, but before solving it, then a reduction in size can be achieved [33]. During the pre-solve process, an entirely new optimization problem is constructed. After the algorithm model has been optimized, the optimal solution is valid for the original problem [34].

Illustrative Example of Deletion Presolve Using the IEEE-14 Bus System
The size-reducing transformation, which leads to an underdetermined nonlinear model, is shown using the IEEE-14 bus system ( Figure 2). The observability constraint function f : R n → R n in case of the IEEE 14-bus system is as follows [12]: (1 -)(1 - As can be observed, some equality constraints are dependent on others and they are redundant constraints for the system of concern. A redundant constraint can be chosen as the one to be eliminated, for it can be reconstructed from those remaining [33]. Such a non-binding constraint could just as well be left out of the model, as this would not affect the optimum point [33].
For example, f 2 is obtained by multiplying f 1 = 0 by additional terms yielding the Esq.
The f 2 is unnecessary since the constraint is always satisfied i.e., f i = 0 and it is removed. When a constraint is satisfied f i = 0 , it can be removed from the system of interest. Let us consider the equality constraints f 12 , f 13 so that: The product f 13 is resulting by multiplying the equation f 12 with an additional term and thus, it is always satisfied. Then, the constraint function involves only mutually independent constraints expanding the feasible region without changing the solution structure. The constraint function is as: where the constraint function h : R n → R m (n > m) is considered to be an underdetermined system of equations defined on a Box Constraint set, that is, A solution that verifies the nonlinear system of equations and optimizes the quadratic objective function has been determined.

Simulation Results and Discussion
In this section, we solve the proposed nonlinear model by using the RQP method. To prove the effectiveness of the RQP to achieve optimality, the obtained objective value is compared to the one found by solving the Integer Linear Programming (ILP) model with binary-valued variables for each benchmark test system. The pure ILP model is solved by using the BBM. In this method, a relaxed continuous linear programming (LP) problem is formulated by ignoring the integer constraints. The relaxed problem is comprised of continuous variables. If the solution to the above relaxed continuous problem contains only involves integers, it is the optimal solution. If the solution has non-integer variables, one non-integer variable is selected and two sub-problems (two branches) are generated [39,40]. For one of them, a ceil constraint is added to the selected non-integer variable and, for the other branch, a floor constraint is added to it. If a branch has no feasible solutions, this branch is terminated [39,40]. If the solution only has integers, it becomes a candidate for the final optimal solution. If the solution has non-integer values, the process of the branching is again solved with these additional constraints. Once the branching process is completed for all the variables, all the integer solutions obtained by the difference branches are compared. The best solution is considered the final optimal solution as shown in the flowchart of Figure 3.  Several MILP routines implement the static-BBM, namely a BBM, including Cut Generation and Heuristics approaches via MATLAB intlinprog optimizer [38], a spatial BBM implemented by solving constraint integer programs (SCIP) [41], an LP-BBM using a primal-and-dual simplex implemented by MOSEK [42] and Gurobi optimizer [43]. Starting from the initial default point, MILP solvers can get only one optimal solution. This is because MILP employs the BB algorithm, and the solution is only updated when a better one is found. This means that solutions with the same results (same objective function value) are overlooked. Without considering other possible optimal solutions, the existing MILP solver stops searching and delivers the result when one optimal solution is found [39,40].
The optimal solutions derived by static-BBM are listed in Tables A1 and A2 of Appendix A for each IEEE bus system [44]. As observed, each static-BBM results in a different optimal solution with the same number of PMUs for each IEEE bus system [44]. On the other hand, the RQP solves the proposed nonlinear program. The RQP is a gradient-based algorithm implemented through a fmincon optimizer routine [38]. To accelerate the convergence towards optimality, the gradient's vectors are given analytically to the optimizer routine [38]. The RQP algorithm model is applied to standard IEEE test systems [44]. The simulation results indicate that the RQP detects the minimum number and the optimal locations of PMUs to make the power system completely observable. Now let us apply the RQP algorithm to the IEEE-300 system [44]. The nonlinear program is as follows: Several MILP routines implement the static-BBM, namely a BBM, including Cut Generation and Heuristics approaches via MATLAB intlinprog optimizer [38], a spatial BBM implemented by solving constraint integer programs (SCIP) [41], an LP-BBM using a primal-and-dual simplex implemented by MOSEK [42] and Gurobi optimizer [43]. Starting from the initial default point, MILP solvers can get only one optimal solution. This is because MILP employs the BB algorithm, and the solution is only updated when a better one is found. This means that solutions with the same results (same objective function value) are overlooked. Without considering other possible optimal solutions, the existing MILP solver stops searching and delivers the result when one optimal solution is found [39,40].
The optimal solutions derived by static-BBM are listed in Tables A1 and A2 of Appendix A for each IEEE bus system [44]. As observed, each static-BBM results in a different optimal solution with the same number of PMUs for each IEEE bus system [44]. On the other hand, the RQP solves the proposed nonlinear program. The RQP is a gradient-based algorithm implemented through a fmincon optimizer routine [38]. To accelerate the convergence towards optimality, the gradient's vectors are given analytically to the optimizer routine [38]. The RQP algorithm model is applied to standard IEEE test systems [44]. The simulation results indicate that the RQP detects the minimum number and the optimal locations of PMUs to make the power system completely observable. Now let us apply the RQP algorithm to the IEEE-300 system [44]. The nonlinear program is as follows: We set an arbitrary initial point through MATLAB command "rand" [38]. Then, the RQP generates the exact solution after 36 iterations, as shown in Table A3 of the Appendix A. The most notable is that no constraint violation exists. The numerical results displayed in Table A3 include the number of iterations, the total number of function evaluations, the convergent sequence of points {x k } ∞ k=0 → x * to the optimum point and the optimal objective function value. The termination is succeeded based on the design criteria in terms of objective function evaluation, the measured first-order optimality at the solution point and the calculated constraint violation. Table A3 illustrates that the step-length is fully being taken to enforce global convergence to the optimum point. As a consequence of that, the unit step-length is accepted so that the Maratos effect does not occur, ensuring a superlinear convergence rate to that optimum point.
Also, the norm of step and the first-order optimality are presented. During the RQP's iterations, the penalty function maintains the feasibility at the solution point. The objective function is minimized on the product x * in running time equal to 17.796910 sec whereas the solution point satisfies the constraint function. The nonlinear model is computationally efficient in finding the required PMUs for the IEEE-300 bus system, as shown in Table A3 of Appendix A and Figure 4.
We set an arbitrary initial point through MATLAB command ''rand'' [38]. Then, the RQP generates the exact solution after 36 iterations, as shown in Table A3 of the Αppendix A. The most notable is that no constraint violation exists. The numerical results displayed in Table A3 include Table A3 illustrates that the step-length is fully being taken to enforce global convergence to the optimum point. As a consequence of that, the unit step-length is accepted so that the Maratos effect does not occur, ensuring a superlinear convergence rate to that optimum point. Also, the norm of step and the first-order optimality are presented. During the RQP's iterations, the penalty function maintains the feasibility at the solution point. The objective function is minimized on the product * x in running time equal to 17.796910 sec whereas the solution point satisfies the constraint function. The nonlinear model is computationally efficient in finding the required PMUs for the IEEE-300 bus system, as shown in Table A3 of Appendix A and Figure 4.
The optimization is completed because the relative first-order optimality measure, 8.8015×10 -11 , is less than the selected tolerance TolFun. Additionally, step-tolerance TolX specifies the termination tolerance for the design variable, as shown in Figure 5. Note that TolFun is a tolerance for both: the size of the latest change in the objective function value and the value of the first-order optimality  x ∈ {0, 1} n is satisfied. The optimization is completed because the relative first-order optimality measure, 8.8015 × 10 −11 , is less than the selected tolerance TolFun. Additionally, step-tolerance TolX specifies the termination tolerance for the design variable, as shown in Figure 4. Note that TolFun is a tolerance for both: the size of the latest change in the objective function value and the value of the first-order optimality measure [38]. The RQP detects a feasible solution satisfying all of the constraints in which the objective function takes it's a minimal value, whereas the Euclidean norm of step is less than 5.50082 × 10 −11 (Figure 4). The norm of step is the current step-size, that is the last displayment in x. Therefore, the termination is succeeded based on the tolerance criteria in terms of objective evaluation, maximum constraint violation i.e., TolCon results in zero and first-order optimality [38].
The RQP identifies alternative solution points over the feasible region constituted by the design variable bounds and constraints, as shown in Table A4 of Appendix A. The number of function evaluations measures the convergence rate [35]. Provided that the function evaluation is the most time-consuming part of the algorithm, this process characterizes the convergence speed [35]. Using approximations differences, the RQP detects the solution point but takes longer to achieve optimality as shown in Table A5 of the Appendix A. The best performance would be obtained by combining analytically given gradient constraints with the RQP method [35].
The RQP requires only a small number of iterations and function evaluations to terminate successfully. The elapsed time overall outcomes by RQP are comparable to those spent by BBM, as shown in Table A5 of the Appendix A.

Performance Evaluation and Comparisons
In this paper, for unraveling the OPP problem, a minimization model is constructed in which the objective function is minimized under an underdetermined nonlinear system of equality constraints defined on a compact set. The proposed algorithm model is solved based on a hybrid-method coupling, a BBM and a gradient-based RQP method. The first phase of the optimization involves a static-BBM that achieves global optimality in solving the convex MILP model [39].
BBM is an efficient enumeration procedure for examining all possible integer feasible solutions. The pure ILP is relaxed to linear programming (LP), namely LP relaxations by removing the integrality conditions [39]. The concepts of branching, bounding, and fathoming are used to obtain the final solution in the process of building the enumeration tree [40]. If the LP optimum solution satisfies the integer requirement, the IP problem is solved. Otherwise, the LP objective value becomes the initial upper bound on the IP optimal value and the root node is partitioned into two successor nodes (subproblems) by two branches [39,40]. When a feasible binary integer solution is obtained via LP relaxations, an optimum point is displayed [39,40]. The main disadvantage is that each BBM has an alternative strategy to built the inherent-tree structure leading to a nonunique global minimum for the OPP [39,40]. The static strategies are best-first, depth-first, best-estimate, best-projection and breadth-first search strategy [39,40]. The number of LP relaxations being derived affects the PMU locations and the convergence speed justifying the differences among the derived running time by different MILP routines as shown in Table A5.
The second phase involves an RQP method, embedded in the fmincon optimizer, that generates a convergent sequence of estimated points to the global minimum point [38]. The RQP algorithm builds a sequence of non-strict global minima over the feasible region that constitutes by the objective function, the bound constraints and the constraint function.
Each optimum point satisfies the binary restriction as shown in Table A4. The RQP detects the solution point, whereas the first-order optimality measure is less than the selected tolerance without constraint violation. The objective is non-decreasing in feasible directions within the value of the function tolerance without constraint violation. Feasible directions are vectors from the current point that locally satisfy the constraints [35]. At a given point → x , not necessarily feasible, the algorithm generates a new point such that the first-order optimality is satisfied at the point and the step-length parameter is determined to minimize the penalty function [35]. The need for the penalty function stems from the fact that, that RQP algorithms do not maintain feasibility at each step performed by the underdetermined nonlinear programming model [35].
Starting from an arbitrary initial estimate, the step-length is required to enforce global convergence of the RQP method, i.e., the approximation of a point satisfying the necessary Karush-Kuhn-Tucker (KKT) optimality conditions [35]. The global convergence is included via a line-search requiring at each step the decrease of the penalty function whose reduction leads to optimality [35]. The fmincon optimizer computed the step-length to satisfy a certain Armijo condition to decrease the penalty function [38].
When the iterations are outsized, even if the iterative point is sufficiently close to the optimum point, the algorithm cannot guarantee that the step-length is 1 [35]. As shown from the simulation results, the unit step-length is achieved so that the RQP does not suffer from slow convergence near to the solution point, ensuring a superlinear convergence rate.
Hence, the Maratos Effect, where the method makes slow progress because second-order changes in the constraints are magnified and outweigh reductions in the objective function, is avoided. It is observed the fast ultimate convergence due to the full accepted step-length by the algorithm for estimates remote from the optimum point. The RQP presents both global and local convergence properties concerned with the asymptotic rate of the convergence, which indicates the rapidity with which the discrepancy between the iterates and the solution goes to zero. The number of function evaluations measures the convergence rate [35].
Provided that the function evaluations are the most time-consuming part of the algorithm, this process characterizes the convergence speed [35]. Using approximations differences, the RQP detects the solution point but takes longer to achieve optimality. This suggests that the best performance would be obtained by combining analytical given gradient constraints with RQP. The RQP requires only a small number of iterations and function evaluations to terminate successfully.
The function evaluations are significantly reduced due to the fact that the gradients are given analytically to the optimizer avoiding finite-differencing of approximate derivatives [35]. As a result, the algorithm's efficiency, measured by function, gradients evaluations and iterations, is improved concerning past studies [11,12,16]. Thus, the size-reducing transformation, which leads to an underdetermined nonlinear programming model, reduces the total run-time significantly as the power network size grows up regarding past studies [11,12,16].

Conclusions
This paper presents an underdetermined nonlinear programming model having fewer equality constraints than the design variables on a closed and bounded set for the solution of the optimal PMU placement problem. The proposed formulation leads to a simpler algorithm model ensuring complete power system observability compared to existing well-determined methods. The minimization model is tested on standard IEEE test systems. An RQP method is proposed to solve the proposed nonlinear model. The RQP produces a global as well as local convergence towards globally optimal solutions quickly with guaranteed accuracy. Each optimum point is a set of binary values assigned to the design variable. The strategic implications of our proposed method cover a wide range of methodological and applied contributions: (1) It delivers a fully functional solution for the OPP problem, with efficiency and increased effectiveness. (2) It serves as a first test best for a novel methodological approach for cost-efficient solutions to improve the monitoring of the power network across large geographic areas. (3) It has the potential to be integrated with sensor networks and 5G networks as well as advanced Data Miners in order to promote increased performance in energy management. (4) This approach can also be integrated with energy hardware solutions for advanced, low and middle scale smart home and smart city projects.
The items (3) and (4) above also define the future research directions for our research study. We intend to use this methodology for smart home power management applications as well as to analyze the impact of our approach for sophisticated data mining services for smart allocation and storage of energy resources over power grids. We understand that the computational sophistication of our approach and it's technical character can also be very useful for other researchers that are dealing with the same research problems. Last but not least the issue of Smart Grids management will require more efforts in the near future with the evolution of smart homes and smart cities applications.