On the minimal set of controllers and sensors for linear power flow

We consider a linear power flow model with interval-bounded nodal power injections and limited line power flows. We determine the minimal number of power injections to control based on a minimal set of measurements, such that the overall system is feasible for all assignments of the non-controlled power injections. For the important case where the possible measurements are the nodal power injections, we show that the problem can be solved efficiently as a mixed-integer linear program (MILP). When also line power flows are considered as potential measurements, we derive an iterative, greedy algorithm that provides a feasible, but potentially conservative solution. We apply the developed algorithms to both a small microgrid and a modified version of the IEEE 118 bus test power system. We show that in both cases a sparse solution in terms of the number of required controllers and measurements can be obtained. Moreover, the number of required measurements can be reduced significantly if line flow measurements are considered additionally to nodal power injections.


I. INTRODUCTION
Volatile renewable energies are transforming classical power grids with few large generators into complex cyber-physical networks. These contain a large number of distributed generators and controllable loads, and power lines are often operated close to their limits. In this context, we ask: What is the smallest set of generators and/or loads that must be controlled based on the values of a minimal number of measurements, such that (s.t.) the entire system state is feasible for all possible values of the remaining elements?
Being able to identify the (optimally small) set of critical elements in complex power grids reduces the cost and effort for their control. Moreover, it is an important ingredient to reduce such systems' potentially high vulnerability with respect to (w.r.t.) natural disasters or cyber-attacks [1], enhancing their operational resilience. An increased protection status could be mandated for the identified critical elements, to keep the number of outages and failures in this group at a minimum, see [2] where the hardening of power systems to minimize system damage in case of disasters is examined. This work was sponsored by the German Federal Ministry of Education and Research in project AlgoRes, grant no. 01|S18066A. It has been performed in the context of the LOEWE center emergenCITY.
Our research question is an instance of the well-known optimal input/output selection problem, also known as the optimal actuator/sensor placement problem. Starting with classical work on controllability [3] this problem has attracted long-term research attention, in particular, for linear timeinvariant systems. The problem has recently become very active again in the study of complex networks, see, e.g., [4]. While most formulations of the problem are NP-Hard due to its combinatorial nature, finding only the minimum set of actuators is possible in polynomial-time [5]. This finding is based on structural controllability theory [6] and can be used to develop distributed algorithms for finding the minimum number of controlled and measured nodes [4]. Structural controllability theory can also be used to analyze cyber-security aspects in distributed power grids [7], e.g., for evaluating the detectability and identifiability of hacked nodes. Another line of research aims at designing control structures that minimize the control effort, using controllability metrics derived from the controllability Gramian of the system [8]. Many of the related input/ouput selection problems are submodular which implies that greedy algorithms using these metrics, e.g., for the optimal placement of High-Voltage direct current lines in a simplified model of the European power transmission network, have provable suboptimality bounds [9]. Time-varying minimal configurations of sensors and actuators can be computed with the help of semi-definite programming [10]. All these works are valid for linear (dynamical, algebraic) systems without state or input/output restrictions.
In this contribution, we propose an alternative, novel approach based on the steady-state representation of the system only, but considering constrained variables. This is an important step towards real applications where power injections and line flows are always subject to physical limits.
Our approach extends current work on the distributed control of power systems [11]. For instance, the robust optimal power flow algorithm by [12] allows computing set points and droop constants for some generators while guaranteeing feasible grid operation for all power injections of other uncertain producers and consumers. While we use similar modeling, we focus on identifying the minimal sets of controllers and measurements that are required for computing such set points.
The rest of the paper is organized as follows. Section II introduces the employed linear power flow model. The feasibility of a given set of controllers and sensors is defined in Section III. We also give a formal problem statement there as well as further computationally advantageous conditions for testing feasibility. In Section IV, we exploit those conditions for developing two efficient algorithms that minimize the number of controllers and sensors. In Section V, we apply the proposed algorithms to find the smallest number of controllers and sensors for 1) a simple microgrid consisting of 4 buses and 2) a modified version of the IEEE 118 bus test case. Finally, concluding remarks and an outlook for future research are provided in Section VI.

II. LINEAR POWER FLOW
We analyze an electrical network with N electrical buses connected by T transmission lines under the common DC power flow assumptions [13]. The voltage phase angles θ ∈ R N determine the nodal active power injections p I ∈ R N and the active power line flows p F ∈ R T as where the entries of B I ∈ R N ×N and B F ∈ R T ×N are defined element-wise as B I,jk = −b jk if j = k, B I,jj = k b jk and B F,jk = b jk , with b jk the susceptance of the line connecting buses j and k.
Without loss of generality, we assume that exactly one generator or load is connected to each bus, with an externally defined active power set point x i . If the sum of the set points in the grid is not balanced, a droop-based primary control scheme [13] adjusts power injections p I under adaptation of the frequency to achieve this balance, such that in steady state we obtain Here, k ∈ R N represents the vector of droop constants, k i ≥ 0 and i k i > 0, and ∆ω ∈ R the frequency deviation with respect to its nominal value. This common setup implies that the measurable quantities p I , p F , and ∆ω are linearly determined by the controllable quantities x. The kernel of the Laplacian matrix B I contains only the constant vectors for connected graphs, that is, a constant shift of the phase angles has no impact on p I . We thus fix θ 1 = 0 and delete the first column of B I to obtainB I . The remaining dimensions of θ are denoted byθ. We similarly reduce B F toB F . The image ofB I moreover contains all vectors with balanced nodal injections. To handle unbalanced set points x, we add k as the last column. This lets us compute for all x with · denoting zero entries In real systems the nodal injections p I will be limited above and below by the technical capabilities of the connected generator or load. Valid set points x might be restricted to smaller intervals than the p I , to leave some space for power generation scheduled by the primary controller. Similarly, line power flows p F and the frequency deviation ∆ω are typically subject to upper and lower bounds.

III. FEASIBLE SETS OF CONTROLLERS AND MEASUREMENTS A. Feasibility Conditions & Problem Statement
The power flow model of the previous section can be abstracted as follows: let x ∈ X ⊆ R N be the variables that can be set externally. X is assumed to be a product of intervals, i.e., X = [x 1 , x 1 ] × · · · × [x N , x N ]. Variables x can be partitioned into the controlled variables x c , for which we will design a controller in the following, and the free variables x f , that are left free to be determined either by other users, cooperative or malicious, by fixed external conditions, such as e.g. the weather, or at random. The index set of the controlled variables is denoted by C and the corresponding partitions of X as X c and X f . We assume that the variables x determine the system state uniquely and that the set of feasible system states X * can be characterized via a set of linear inequalities, where A ∈ R K×N and b ∈ R K . Similarly, we assume a set of possible measurements y ⊆ R L to be linearly related to the system state, i.e., y = Mx with M ∈ R L×N . We partition these possible measurements into the monitored measurements y m , that are used as inputs to the control law, and the unmonitored variables y u , that are not required for the controller and may or may not be recorded in practice. The index set of the monitored variables is denoted by M.
The defined partitions of x and y allow to partition the matrices A and M along their columns or rows as well, yielding The aim of the paper is to determine the minimal set of controllers C and measurements M that allows for the design of a control law x c (y m ) that can guarantee a feasible system state, independently of the state of the free variables x f . This can be formalized as follows.
Definition 1 (Condition C 1 ). Sets C and M are feasible if The idea behind this definition is that the control x c (y m f ) chosen for y m f should be valid for the x f from which y m f originated. Note that we consider only control values for the steady state of the system in this paper. We do not examine whether and how it is possible to get there from arbitrary initial positions. Moreover, in order to simplify the notation of the involved sets, we have used only a part of y m as input to the control law x c (y m f ). However, since y m = M m c x c + y m f one could easily rewrite the controller into the form x c (y m ), i.e., directly using the measurements that are actually available to the controller.
Since x f uniquely determines y m f , we can also express the control law as x c (x f ). The formulation x c (y m f ) implies that x c will attain the same value for all values of x f that lead to the 21st Power Systems Computation Conference same measurements. We thus obtain the following equivalent condition.
These definitions allow us to state the optimization task we aim to solve in this work.
Problem statement. Find the set of controllers C and measurements M that solves |C| and |M| denote the cardinality of C and M. The cost of placing a sensor is weighted by 0 ≤ γ ≤ 1 since it will typically be smaller than implementing a full actuator.
One could additionally incorporate into the objective the varying efforts and costs for controlling certain elements or acquiring certain measurements. Instead of just weighting the total number of controllers and measurements we would then determine an individual weight for each element separately. While we do not follow this idea below, all algorithms could straightforwardly be adapted.

B. Related Conditions
Verifying conditions C 1 and C 1 based on their definition requires checking infinitely many values of y m f or x f , respectively. We therefore derive two related conditions that are testable with finite computational resources. The relation of all derived conditions is presented in Fig. 1. In the next section we then show how to exploit them to efficiently solve our problem.
Condition C 1 requests the existence of a mapping x c : M m f (X f ) → X c yielding valid control values. One possibility is that this mapping is affine-linear. where It is, however, not necessary as can be shown by counterexample, where piecewise linear control laws sometimes allow for fewer sensors and controllers. The condition is testable with finite efforts, as we show in Section IV.
The conditions presented so far are continuous in the sense that testing their validity requires checking an infinite set of possible realizations of x f or y m f . However, since the possible values of x f and y m f are restricted to bounded polytopes, i.e., X f and M m f (X f ), we can derive a necessary condition for C 1 based only on the corners of such polytopes. In contrast Fig. 1: Relation of the desired conditions C 1 and C 1 to the conditions C * 1 , C 2 , which are testable with finite resources.
to C 2 , such necessary condition will not assume the control law u c (y m f ) to be affine-linear. Definition 4 (Corner). z ∈ Z is an extreme point or corner of the convex set Z if there are no two distinct points z 1 , z 2 ∈ Z and λ ∈ (0, 1) such that z = λz 1 + (1 − λ)z 2 .
Denote C(X f ) as the set containing the corners of X f . The number of corners of X f , denoted as |C(X f )|, is finite, but grows exponentially with the number of free variables. A condition based on all corners of X f would therefore be computationally prohibitive for larger dimensions of X f . Instead, we focus on a subset of corners only, namely those ones which have the maximum impact on the constraints Remark 1. Note that the optimization problem (9) defining the elements of C A (X f ) can be solved analytically for row A i as The optimal values thus depend only on the sign of the corresponding elements of A. In many cases the optimal vectors for different rows of A will therefore coincide and the cardinality of C A (X f ) is even smaller than its maximum possible value K.
The reverse is not always true, as can be shown by counterexample. However, we will show below that this condition can be exploited for a very efficient computation of approximate sets C and M, at least for the case when the set of possible measurements consists of the power set points at each node, i.e., when M is an identity matrix of appropriate dimensions, here denoted by I. For nodes with zero droop constant, e.g., typical loads or small-scale generators, the measurement of the power set points is equivalent to measuring nodal power injections.
Minimizing the objective C+γ|M| with respect to condition C * 1 or C 2 will provide a lower or upper bound for the optimal solution of problem (7), respectively. In our experiments we found that for minimal sets C and M fulfilling C * 1 it was often possible to determine valid affine-linear control realizations by testing C 2 for such sets, i.e., the upper and lower bound coincided. In this case, C and M are optimal solutions of (7).

IV. ALGORITHMS
The feasibility conditions formulated above enable us to develop two methods for addressing optimization task (7): The first, derived from condition C * 1 , leads to a mixedinteger linear program (MILP) that finds the smallest feasible sets C and M, provided that M = I. Since C * 1 is necessary for C 1 but not sufficient, the obtained sets C and M may be too small to be feasible. While we often obtained feasible results anyway, the algorithm can also be used to generate a good initial solution for the second approach.
The second method for solving problem (7) is designed for all possible measurement matrices M. It is a greedy procedure based on hill climbing (HC) and condition C 2 . Recalling that condition C 2 is sufficient for C 1 but not necessary, the obtained sets C and M may possibly be too large, but are guaranteed to be feasible.

A. MILP-Based Approach
In this section, we develop a MILP for finding the smallest feasible sets C and M based on condition C * 1 , provided that M = I. The key is to formulate condition C * 1 as a set of linear inequalities that holds for all choices of sets C and M.
To this end, consider the binary decision variables u c ∈ {0, 1} N and u m ∈ {0, 1} N , defined element-wise as The decision variables u c and u m encode the elements of C and M, respectively. Finding the smallest number of elements of C and M is thus equivalent to minimizing the cost u c 1 + γ u m 1 . Letx i ∈ X , i = 1, . . . , K, be defined element-wise as x i is the A i -optimal analytical solution of (9) for the case when all variables are assumed to be free. Moreover, for any given set of controlled variables C, the elements of where 1 is a vector of ones of appropriate dimension and • represents the Hadamard product. Since we assume here that M = I, we can further partition the free variables into monitored and unmonitored variables, i.e., we can write 1 − u c = u m + u u , with u u ∈ {0, 1} N being the binary vector that encodes the elements of the unmonitored variables. The A j -optimal corners of X f can then be identified withx j • u u , j = 1, . . . , K. Similarly tõ x i • u m , a vector of length N whose non-measured entries are zero, we now consider an associated control vectorx i c ∈ X for which i.e.,x i c is a vector of length N whose non-controlled entries are zero.
Condition C * 1 states that for all i = 1, . . . , K, the control vectorx i c for the A i -optimal cornerx i f of X f should be valid and that it should be the identical to the control vector for all other corners in C A (X f ) that cannot be distinguished given the measurements. We thus can consider only the worst case of the unknown elements and write compactly withÃ i m and A u defined element-wise as, k = 1, . . . , K, Thus, the mixed-integer linear program that solves (7) when Remark 2. Note that, particularly in large scale applications, there may be several constraints, i.e., rows of A and corresponding entries of b, that are not violated for any realization of x. Hence, when optimizing (12), we only take into account the rows of A, for which a violation of (4) is possible, i.e., where A ixi − b i > 0. This preprocessing is also utilized by the greedy search proposed below.

B. Greedy Approach
In this section, we first show how to check condition C 2 efficiently via a linear program (LP) for fixed sets C and M. Thereafter we describe an iterative algorithm to choose and adapt these sets in order to find minimal feasible sets.
For given C and M, condition C 2 mandates to check if there exists a valid affine-linear control law that makes the system feasible for every possible value x f ∈ X f . More precisely, there should exist an affine-linear control law defined via S and w such that for all x f ∈ X f we have where we introduce η ∈ R as an indicator of how far the system is from being infeasible. A control law is valid if η ≤ 0.
To tackle condition (13) for all x f ∈ X f we only need to consider the maximum of the left hand side expression. Let K = K + 2|C| be the number of rows ofÂ(S) and N f = N − |C| the number of free variables. We can introduce an upper bound onÂ(S)x f via a matrix H ∈ RK ×N f , whose entries fulfill H ij ≥Â ij (S)x fj , 21st Power Systems Computation Conference PSCC 2020 for all i = 1, ...,K and j = 1, ..., N f . The upper bound of A(S)x f is then given by H1 and condition (13) is equivalent to H1 + Fw − l ≤ ηv.
Putting these results together allows us to compute the minimum possible value of η for given M and C via the following linear program The above described algorithm for testing the validity of C 2 for fixed C and M can now be used as a subroutine to minimize over the sets C and M as well. To do this, we proceed iteratively from initial sets C and M adapting them one element at a time. Since we want to measure the optimization progress also for non-feasible combinations C and M, we extend the minimization objective to where η is the feasibility indicator obtained from solving problem (16). µ > 0 is a weighting factor that penalizes the infeasibility of C and M. We choose µ 1 to steer the iteration quickly towards feasible solutions.
The cost function (17) is minimized via a greedy hill climbing procedure. In each iteration we compute the objective value for all sets M or C that can be generated by adding one element to either M or C. We then choose the step which yields the largest improvement of the objective value (17). As soon as the sets of controllers and measurements are feasible, we stop the iteration.
It is well known that the solution of this greedy approach depends on the selection of the starting point. A natural option is to start with empty sets, selecting the most important controllers and measurements during the first iterations. Alternatively, we propose to use the MILP (12) formulation as an initial guess. More specifically, we solve the MILP (12) for M = I first. We then use the found controller set C as a starting point for the greedy approach, while disregarding the found measurements. Instead, we start with an empty M. This way the measurements resulting from general M, which potentially allow for more compact control systems than the identity measurements, can be integrated well, but the critical controllers are already identified. Remark 3. Since general MILP has exponential worst-case time complexity, this is an upper bound on the complexity of our first approach (12). In contrast, LP as used for our second approach (16) is known to have polynomial worstcase time complexity, and the hill climbing procedure only adds polynomial factors. However, for the realistic examples discussed in the next section we found the MILP approach to be more efficient than the hill climbing procedure. The latter's computation time depends strongly on the starting point. For the examined medium to large problem instances, it allowed finding small, guaranteed to be feasible solutions for M and C with very reasonable efforts. For cases when M = I the MILP solution could often be verified to be feasible (and thus also optimal) by solving the small LP (16) only once without further adaptation of M or C. We thus see both algorithms as an important contribution for solving real control design problems with state constraints.

V. NUMERICAL EXAMPLES
The algorithms developed in IV-A and IV-B are now applied to find the minimal feasible configuration of controllers and measurements for two exemplary power systems. We first demonstrate our setup and typical effects on a simple microgrid of 4 buses connected in a line. Subsequently, a modified version of the IEEE 118 bus test case is addressed. The experiments were performed using an i5 notebook with 8 GB of RAM. The algorithms were implemented in Matlab R2018b, using YALMIP [14] as modeling language and CPLEX 12.9 as LP and MILP solver. Fig. 2 shows the considered microgrid consisting of three generators supplying a demand of 5 MW. It gives the topology of the grid together with the capacity limits of each transmission line and each generator/load. The generator located at bus 4 provides primary reserve, initially with a droop of 12 MW/Hz and later with 4 MW/Hz. The maximum allowed frequency deviation is ±0.1 Hz. We first assume that all transmission lines have a power transfer capacity of ±10 MW, which is adequate to avoid grid limitations. In scenario (d) we add an active line constraint in the middle.

A. Simple Microgrid
In scenario (a) where only the power set point at each bus may be measured, it is sufficient to control the large generator located at bus 4 for achieving feasible grid operation. The set points of the remaining smaller generators can be chosen freely and no additional measurement devices are required.
In scenario (b) we reduce the droop of the generator at bus 4 to 4 MW/Hz. This makes the measurement of the power injections at buses 1 and 2 necessary. Although the power injections at buses 1 and 2 can be chosen arbitrarily, they must be monitored so that the power produced by the generator located at bus 4 can be set appropriately to balance the system within the given frequency tolerance.
In scenarios (a) and (b), where M = I, the solutions of the MILP were feasible (and optimal) without further adaptation of M and C and the greedy approach, starting from empty sets M and C, produced the same results.
In scenario (c) the measurement of the line flows and the grid frequency is added to the set of potential measurements, when performing the greedy optimization. This allows to reduce the number of measurements to only one. For this scenario the solution is not unique: one possibility is to take the measurement of the frequency deviation as controller input, yielding an adapted primary control scheme. An alternative solution that is shown in the figure is to monitor the sum 21st Power Systems Computation Conference PSCC 2020 Porto, Portugal -June 29 -July 3, 2020   x 1 ;  For each scenario, the resulting (non-unique) affine-linear control realization is provided below together with the behavior of the potentially active constraints for all x f ∈ X f on the right. For scenario (c) with multiple, equivalent optimal solutions, the colored frames denote alternative optimal solutions. The rationale behind the scenarios is as follows: In (b) the primary control droop is reduced compared to (a). In (c) we allow for additional measurements. In (d) we reduce the transfer capacity of the middle link to form an additional active constraint.
of the outputs of generators 1 and 2 by measuring the line flow between bus 2 and 3 for controlling the set point of the generator at bus 4. This situation will be very common in future active distribution grids, where individual small scale loads or generators are not able to violate local grid constraints, but their aggregated effect is important to the system. Since the load is fixed, measuring the line between buses 3 and 4 would be equally informative. The feasibility of all these solution candidates was verified via LP (16), obtaining valid affine-linear control realizations in all cases.
In scenario (d), we constrain the capacity of the transmission line connecting buses 2 and 3 to the interval [−1, 1] MW. This represents an active grid constraint if the generators at buses 1 and 2 produce at maximum power. The solution obtained via hill climbing optimization consists of additionally controlling the power injection at bus 1. Again, several alternative solutions are possible.
In scenarios (a), (b), and (c), the frequency deviation represents an active constraint to the operation of the system. Observe in Fig. 2   linear control law keeps the frequency deviation inside the feasible region for all values of the non-controlled injections.
In scenario (d), the designed controller also ensures feasible system operation despite the limited power capacity of the middle line. While for the demonstrated example all solutions can readily be verified manually, it shows that the situation may become much more complex in larger grids. The topological location of generators and loads in the grid is important as well as their capacity and their neighborhood. An automated algorithm for selecting critical elements to control and/or measure is thus very beneficial for complex networks with distributed generation and transmission lines that are operated close to their technical limits.
In scenarios (c) and (d) the use of the greedy approach is required to deal with M = I. Using the MILP solution as an initial guess for C or starting with empty sets led to the same optimal objective function value. The solutions for M and C did not always agree exactly, but could be shown to be equally optimal.
The total solver time for all scenarios is shown in Table  I. As expected, the MILP optimization performs faster than the hill climbing optimization for the same instances. When computing the optimal sets for scenarios (a) and (b), the MILP algorithm was more than 2 times faster than the hill climbing with empty sets. It was also 1.2 times faster than the hill climbing that uses the MILP solution for C as initial guess, which corroborates the benefits of such concatenated optimization procedure.

B. IEEE 118 Bus Test Case
We now analyze the modified version of the IEEE 118 bus test case, see Fig. 3. This power system is composed of 54 generators, 99 loads, and 186 transmission lines. The topology of the power system, the load values and the line and generator capacities were taken from [15]. We assume that each generator can be scheduled in the range of 10%-90% of its available capacity. In addition, we admit 10% of uncertainty of each load in both directions. The maximum allowed frequency deviation is taken as ±0.2 Hz.
We first consider the case when only the power set points may be measured, i.e., M = I, see Fig. 3a. We obtain an optimal set of 12 controller and 20 measurement devices to guarantee feasible grid operation. The remaining 96 injections can be left operating free and/or be manipulated deliberately and do not require any monitoring equipment.
To obtain this result, we first use the MILP algorithm and then validate its solution via LP (16). The obtained η is smaller than zero, thereby proving the feasibility and optimality of the MILP solution. When we initialize the greedy search with empty sets, we obtain a feasible solution consisting of 23 controllers and 9 measurements. As expected, the obtained solution in this case is larger than the one provided via MILP optimization. This confirms that taking the MILP solution as initial guess is beneficial for the greedy search.
We now add the measurements of the line flows and the frequency deviation into the set of potential measurements, see Fig. 3b. This yields in total 305 possible sensor devices.  We first apply the MILP algorithm and then the greedy one, starting with the controllers identified via the MILP. As expected, the solution is much sparser than before. The total number of required sensors is reduced from 20 to 3. The selected line flows confer a large amount of information that help avoiding grid capacity violations. It is insightful to observe the progress of the hill climbing procedure: buses with major generators connected are selected as controlled nodes first. The procedure is thus initially reducing the impact of the free variables on the system by controlling the highest uncertain injections first. When enough controlled nodes were selected, the selection of measurements starts to be significant for the minimization of the cost. Selected measurements are often related to nodes connected either to large non-controlled generators or to high uncertain loads. Remaining buses with smaller injections are mostly left unobserved. Table II shows the obtained solver time for all studied cases. The solution for case (a) using MILP optimization was found in about 2.57 minutes. Observe that the solution for case (b) was computed in about 28 min for the concatenated execution of both algorithms, compared to the ca. 154 minutes needed by the solver when starting hill climbing with empty sets. A single verification step using LP (16) took less than a second.
The computation time could further be improved, e.g., by testing not all possible set extensions in each step of the greedy search but using only a representative subset, selected by proximity in the graph. Another idea would be to add more than one element in each iteration. For the control design task described in this paper, however, the achieved computation time seemed acceptable even without these extensions.

VI. OUTLOOK
The theoretic framework and the algorithms developed in this work allow for the efficient identification of critical controllers and measurements in complex power systems with uncertain producers and consumers. Unlike previous work, we take specific power limitations of lines, generators, and loads into account. This step strongly improves the applicability in practice, where our approach will help reducing control costs and efforts and increasing power systems' resilience.
While we have only considered active power in this work, the approach can straightforwardly be applied to linearized power flow models taking into account also reactive power and voltages. Developing a MILP formulation for condition C 2 is also possible, but our experiments so far have not yielded satisfying run times.