The Basic Algorithm for the Constrained Zero-One Quadratic Programming Problem with k -diagonal Matrix and Its Application in the Power System

: Zero-one quadratic programming is a classical combinatorial optimization problem that has many real-world applications. However, it is well known that zero-one quadratic programming is non-deterministic polynomial-hard (NP-hard) in general. On one hand, the exact solution algorithms that can guarantee the global optimum are very time consuming. And on the other hand, the heuristic algorithms that generate the solution quickly can only provide local optimum. Due to this reason, identifying polynomially solvable subclasses of zero-one quadratic programming problems and their corresponding algorithms is a promising way to not only compromise these two sides but also offer theoretical insight into the complicated nature of the problem. By combining the basic algorithm and dynamic programming method, we propose an effective algorithm in this paper to solve the general linearly constrained zero-one quadratic programming problem with a k -diagonal matrix. In our algorithm, the value of k is changeable that covers different subclasses of the problem. The theoretical analysis and experimental results reveal that our proposed algorithm is reasonably effective and efﬁcient. In addition, the placement of the phasor measurement units problem in the power system is adopted as an example to illustrate the potential real-world applications of this algorithm.


Introduction
Optimization problems normally fall into two categories, one for continuous variables and the other for discrete variables. The latter one is called combinatorial optimization, which is an active field in applied mathematics [1]. Common problems are the maximum flow problem, traveling salesman problem, matching problem, knapsack problem, etc. Among those famous combinatorial optimizations, zero-one quadratic programming (01QP), whose variables can only be either 0 or 1 [2], is very important and attracts a lot of attention.
Zero-one quadratic programming, which can be divided into constrained (01CQP) and unconstrained (01UQP), is a combinatorial optimization and has practical significance. For example, it can be applied in circuit design [3], pattern recognition [4], capital budgeting [5], portfolio optimization [6], etc. Except for these well known applications, 01QP also has potential applications in the nonlinear control related fields [7][8][9][10][11]. Among these applications, the phasor measurement unit (PMU) placement has been widely studied. Phasor measurement units can be used in dynamic monitoring, system protection, and system analysis and prediction. Therefore, the placement of PMU has become an important issue. As the scale of the electric grid grows, the PMU placement problem becomes more difficult and must be addressed considering certain requirements. In [12],
Based on past works [30,31], we can have the algorithm, as follows, to solve 01UQP. The feasibility and effect of the algorithm can be seen in [31]. Figure 1 shows the Algorithm 1 process intuitively.
(1) Adjacent terms are the combination of x n−m+1 ,. . . , x n , whose value is the same except for the value of x n .
(2) Get the corresponding f (x). (1) Compare every two adjacent states and take the result with the smaller constant term as the new state.
(2) Update all the 2 m states.
Step 3: Get the optimal solution x.
(1) Update the states based on Step 2 until only the constant term is in f (x). (2) The optimal value is the minimal one.
(3) Trace back and get the optimal solution x.

01CQP Algorithm Description
Consider the constrained k-diagonal matrix zero-one quadratic programming problem: where Q has the same meaning as that of the 01UQP formula in Section 2, a ∈ Z n + , b ∈ Z + . In this section, we utilize the dynamic programming method to solve the 01CQP problem. To apply the dynamic programming method, we introduce a state variable s k ( s k ∈ Z) and a stage variable k(0 < k ≤ n), which should satisfy the following iteration. s k+1 can be expressed as: We only need to consider the integer point of the state space since a ∈ Z n + and b ∈ Z + . Since s k satisfies 0 ≤ s k ≤ b, we define a set s k = {s k |0 ≤ s k ≤ b, s k ∈ Z}. Algorithms 2 and 3 show the detailed calculation process.

Algorithm 2
Calculation Method of f (s k ) Case 1: When k = 1, there are two cases.
. , x n ) We will get a series of functions f (x) after executing Case 1.

Case 2:
When k ≥ 2, there are also two cases.
(2a) s k < a k x k must be 0 to satisfy s k−1 = s k − a k x k . At this time, s k = s k−1 , by which we can obtain the function f (s k ) = f (s k−1 )| x k =0 . (2b) s k ≥ a k In this case, x k can be both 0 or 1, which generates two more situations: We can see that there is only one function f (s k ) corresponding to each state s k when k = 1, and there are several f (s k ) to each state s k when k > 1. To save storage and computational time, f (s k ) should be selected satisfactorily in the next step. We need to find the optimal f (s , that the optimizing process refers to for Algorithm 3.
, which are almost the same except for the constant term and pick up the smallest one.
There are more than one f (s k ) = f (s k−1 )| x k =1 : 1) Find the f (x) using a similar approach to that in state 0 Case 2.
According to the above algorithm, the maximum number of functions per State is shown in

Analysis on the Effectiveness and Rationality of 01CQP Algorithm
In this section, we analyze the properties of the polynomial time algorithm for 01CQP. To analyze the algorithm, we need to demonstrate the rationality of Algorithms 2 and 3. Suppose f (x) has the form: f (x) = q 12 x 1 x 2 + q 13 x 1 x 3 + · · · + q 1,1+m x 1 x 1+m + · · · + q n−m,n−m+1 x n−m x n−m+1 + · · · + q n−m,n x n−m x n + · · · + q n−1,n x n−1 x n + c 1 x 1 + · · · + c n x n .
Clearly, when k = 1, each state s 1 has only one function f 1 (s 1 ). The coefficient of x 2 , . . . , x m+1 and the constants are the only different terms in the function.
(const represents the constant term of the f (s k−1 ).) If s k ≥ a k−1 , f k (s k ) has the same form as function (4): Then, f k (s k ) can be shown as: At this point, the maximum number of f k (s k ) corresponding to each s k is 2 k−2 . (2) If s k ≥ a k , execute state 0 and state 1 .
Here, the state 0 is the same as the one in Case 1, so we do not need to repeat it. Consider state 1 , has the form of function (6), then, f k (s k ) can be expressed as: x i x i+2 + · · · + q i,i+m x i x i+m + · · · + q n−m,n−m+1 x n−m x n−m+1 + · · · + q n−m,n x n−m x n + · · · + q n−1,n x n−1 x n + c k+1 x k+1 + · · · + c k+m x k+m + c k+m+1 x k+m+1 All f k (s k ) are only different in the constant terms and the coefficient of x k+1 , ..., x m+2 .
Consider the most complex case, s k ≥ a k and the state variable s k corresponds to 2 m−1 f k (s k ) and 2 m−1 f k (s k ) = f k−1 (s k − a k ) respectively (2 m−1 is the number of both). Therefore, we need to execute state 0 first, and then execute state 1 .
(1) state 0 : Since f k (s k ) has the form as function (8), the state 0 is executed for the f k (s k ) respectively, and the following expression is obtained: m,n−m+1 x n−m x n−m+1 + · · · + q n−m,n x n−m x n + · · · + q n−1,n x n−1 x n + c n−m x n−m + c n−m+1 x n−m+1 Clearly, these f k (s k ) are only different in the constants and the coefficients of x n−m .
(2) state 1 : It is possible to obtain f k (s k ), which has the following form: q n−m,n x n−m x n + · · · + q n−1,n x n−1 x n + c n−m x n−m + c n−m+1 x n−m+1 For both case (1) and (2), they are only different in the constants and the coefficients of x k+1 .
Step 4: k = m + 2, ..., n Based on Algorithm 2 Case 2, we calculate the function f k (s k ), and the core is the implementation of state 0 and state 1 . When executing state 0 and state 1 , the difference between the two f k (s k ) is the constant and the coefficients of x k . Therefore, when executing state 0 , we only need to compare the constant terms. When executing state 1 , we only need to compare the sum of the constants and the coefficients of x k . In this way, we can avoid solving the excess f k (s k ).
The form of the function f k (s k ) and the maximum number of f k (s k ) in each stage k are similar to the one in stage k = n − m, which also ensures the repeatability of the algorithm. The algorithm is supplemented by concrete examples and numerical simulations.

Calculation Example of 01CQP
For example, the parameters Q, c, a, and b are: In this example, we will omit some of the states.
(3) Set k = 4, 5, 6, 7 According to Algorithm 2 Case 2 and Algorithm 3, we can gain f k (s k ) with k = 3, ..., 7. Tables 4  and 5 show part of the calculation process. From the table, we can see that only the coefficient of s k+1 and constant terms are different in the adjacent items.

Application of Zero-One Quadratic Programming in Phasor Measurement Units Placement
How to find the installation location and the number of installations is the focus of the phasor measurement units placement problem. The matrix in the PMU placement problem differs from the previous k-diagonal matrix in the elements of the main diagonal. However, since x is either 0 or 1, we can convert it into a k-diagonal matrix. Then, we can take advantage of the algorithms designed above to work out PMU placement without considering too many constraints and realistic requirements.

Modelling
For the power system composed of n nodes, the placement of PMU is represented by n-dimensional vector X = (x 1 , x 2 , . . . , x n ). Here, i = 1, 2, . . . , n, The matrix H represents the network graph structure, and the objective function is where λ is the weight. N ∈ R n represents the upper bound of the maximum redundancy of each bus. R ∈ R n×n and Q ∈ R n×n are diagonal arrays, representing the importance of each bus and the cost of placing the PMU on the bus respectively. We expand Equation (11), Then, Equation (12) can be expressed as integer quadratic programming as, where . A i and Ω represent a set of nodes adjacent to the bus i and the bus set respectively. The above constraints require at least one PMU unit to be placed between the bus and its adjacent nodes to ensure that each adjacent node of the bus can be observed. The above problem can be expressed as unconstrained problems by the weighted least square method: where V = diag(v i ).

Example and Experimental Result
Suppose the coefficient matrix H of a power system is Set Q as the unit matrix, λ = 0.5, The main diagonal of the matrix for this problem is 1, while the definition of k-diagonal matrix Q requires the main diagonal elements to be zero. Since x ∈ {0, 1} n , the diagonal of G can be converted into a linear term and can be considered as a seven-diagonal matrix. This indicates that the diagonal of G becomes zero and f is updated as seen in Equation (15). Using the algorithm in Section 2, we can find the optimal solution (0, 1, 1, 1, 0, 1) without considering the constraints. This is to install four PMUs in bus 2, 3, 4, and 6. Through observation, it can be seen that the configuration results meet the system observability and the constraint is reached.
The above problem considers the placement of PMU under the condition that the system is completely observable. When considering the PMU placement problem with the N − 1 principle, are added based on the definition of a single node observable. N is a set that includes the nodes required to have objectivity when a fault occurs in the system. P 2 i and P 1 i represents a node set connected to node i with two and one line respectively.
We can get the optimal solution is the same as that in the former example. As the connection between nodes 5 and 6 fails, the system after placing PMU will not lose its observability. We only consider the algorithm to solve the constrained k-diagonal zero-one matrix quadratic programming and consider that it can be applied to some form of PMU placement problems. Some of the conditions to the placement of PMU have not been taken into account here, which has yet to be studied by specialized engineers.

Algorithm Experimental Simulation
Now, the experimental results are provided for this algorithm to illustrate its performance. We implement the algorithm with C++ and run it on an Intel(R) Core(TM) i7-8550U CPU. For the problems we tested, the dimension of matrix Q ranges from 10 to 100 and k takes 5, 7, . . . , 25. All simulation data (Q, c, a, b) are generated randomly. Q, c is set up as in Section 2 with numbers ranging from −100 to 100, and a, b range between 0 and 20. Then, we obtain Table 7, which shows the detailed computation time for different dimensional problems with different diagonal numbers.

Experimental Results Discussion
Here, based on the experimental results, we investigate the influence of n and k on the algorithm and discuss the importance and implications of our study results. Firstly, we fix the values of k but increase n. Figure 2 shows the experimental situations when the matrix dimension n changes from 10 to 100 with k is fixed at 9, 13, 19 and 23 respectively. As can be seen, the calculation time increases slightly with the increase of dimension. Secondly, we fix the values of n but increase k. Figure 3 shows the experimental situations when k changes from 5 to 25 with n is fixed at 20, 40, 60 and 80 respectively. It can be observed clearly that time changes significantly with the increase of diagonal number. Finally, we increase n and k simultaneously. Figure 4 illustrates this situation. It is obvious that k has an obvious influence on the calculation speed. All these observations coincide with the time complexity we derived in Section 3. It means that if the diagonal number is within an appropriate range, our algorithm can perform effectively and efficiently even for very large scale problems. And these kind of problems cover a large portion of the whole problem set. Therefore, the algorithm we proposed will have great potential in many real-world applications.

Conclusions
In this paper, a novel exact algorithm to the general linearly constrained zero-one quadratic programming problem with k-diagonal matrix is proposed. The algorithm is designed by analyzing the property of matrix Q and then combining the famous basic algorithm and dynamic programming method. The complexity of the algorithm is analyzed and shows that it is polynomially solvable when m is fixed. The experimental results also illustrate the feasibility and efficiency of the algorithm. Designing efficient algorithm to this special class of problem 01CQP not only provides useful information for designing efficient algorithms for other special classes but also can provide hints and facilitate the derivation of efficient relaxations for the general problems. And finally, the phasor measurement units placement problem is used to demonstrate that the algorithm has wide potential applications in decision-making real-life problems.
Author Contributions: S.G. put forward ideas and algorithms. S.G. and X.C simulated the results and wrote important parts of the article. X.C. formatted the whole paper, and summarized the results in tables. All authors have read and agreed to the published version of the manuscript.

Funding:
The work described in the paper was supported by the National Science Foundation of China under Grants 61876105 and 61503233.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: 01QP Zero-One Quadratic Programming 01UQP Unconstrained Zero-One Quadratic Programming 01CQP Constrained Zero-One Quadratic Programming PMU Phasor Measurement Units NP-hard Non-deterministic polynomial-hard

Variables
The following variables are used in this manuscript: Q Q = (q ij ) n×n and only q ij =q ji (i = 1, 2, . . . , n − 1, j = i + 1, i + 2, . . . , i + m) are not zero c c = (c ij ) n×1 a a = (a ij ) n×1 b b ∈ Z + n Dimensions of matrix Q m m ∈ [0, n − 1] and m ∈ Z k k = 2m + 1 (m = 0, 1, 2, . . . , n − 1) s k State variable s k ∈ Z, stage variable k(0 < k ≤ n) f (x) f (x) = 1 2 x T Qx + c T x H H = (h ij ) n×n is a network graph structure N N ∈ R n represents the upper bound of the maximum redundancy of each bus R R ∈ R n×n is the importance of each bus Q Q ∈ R n×n in the application is the cost of placing the PMU on the bus M(x) Column vector consisted of m i (x) = (1 − x i ) ∏ j∈A i (1 − x j )(i ∈ Ω) Ω Bus set λ Weight A i Nodes adjacent to the bus i