Optimal sensor placement for linear systems

Bachelorarbeit im Studiengang Mathematik Lehrstuhl für angewandte Mathematik Optimal Sensor Placement for linear Systems eingereicht von: Maximilian Pfister eingereicht am: 17.08.2012 Betreuer: Prof. Dr. Tobias Damm Dr.-Ing. Ulrich Münz (Siemens AG) Abstract The aim of sensor placement is to observe the state of a dynamical system while using only a small part of the available output information. Thus, the observer does not need sensors at every possible node of the system. We use sensor placement because it is not practical for large-scale networks, such as power grids, to place sensors at each node. With an optimal sensor placement we obtain a subset of sensors which minimizes the observer error in comparison to any other subset of the same size. This means we generate an optimal observation with the given number of sensors. We compute the observer error, for the linear dynamical systems we consider, with the H2 -norm of the observer error system. In this approach, we optimize both the subset of selected sensors and the observer gain matrix in parallel. The optimization problem is non-convex both in a constraint, which bounds the H2 -norm, as well as in the objective function which uses a `0 -norm to count the used sensors. To obtain a semidefinite program, we first relax the `0 -norm by an iterative reweighted `1 -norm. Second, we use a reformulation of the H2 -norm with linear matrix inequalities to replace an occuring bilinear and therefore non-convex term. We use this computationally efficient formulation of the sensor placement problem to derive three algorithms. Furthermore, existing algorithms, which do not use the convex reformulation of the optimization problem, were implemented. The algorithms are compared extensively relating to execution time, performance of the chosen sensors, and the applicability on a practical problem. The practical problem is a model of a high-voltage power grid with the aim to measure the phase angles and the frequencies at every node. The result of the comparison is that a algorithm with a greedy approach solves the optimization problem fast and usually with a good solution. However, this algorithm is problematic because the shortsighted greedy approach cannot exclude that a worst case solution is generated. The best results in general were produced by a novel approach made in this thesis. This novel algorithm iteratively solves the relaxed optimization problem and finds near-optimal sensor subsets. Contents 1 Introduction 7 2 Problem Formulation 9 2.1 Problem Formulation for Discrete-Time Systems . . . . . . . . . . . . . . . . . . 9 2.2 Problem Formulation for Continuous-Time Systems . . . . . . . . . . . . . . . . . 10 3 Convex Relaxation of the Optimization Problem 11 3.1 Convex Relaxation of the `0-Objective........................ 11 3.2 Reformulation of the H2-Performance Constraint . . . . . . . . . . . . . . . . . . 12 3.3 Algorithms for Optimal Sensor Placement for Discrete-Time Systems . . . . . . . 14 3.3.1 Algorithm with Substitution . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.3.2 Algorithm with Cone Complementarity Linearization . . . . . . . . . . . . 17 3.4 Algorithm for Optimal Sensor Placement for Continuous-Time Systems . . . . . 19 4 Comparison of Different Algorithms 23 4.1 State of the art Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.1.1 GreedyAlgorithm ............................... 23 4.1.2 Branch and Bound Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.2 RuntimeComparison.................................. 24 4.3 PerformanceComparison ............................... 26 4.4 ApplicationExample.................................. 28 4.5 Concluding Comparison of the Algorithms . . . . . . . . . . . . . . . . . . . . . . 31 5 Conclusion and Future Directions 33 6 Appendix - Mathematical Definitions and Notation 35 6.1 `0-Vector-Norm..................................... 35 6.2 H2-System-Norm.................................... 35 6.3 MatricesandLMIs................................... 36 Bibliography 37 5 1 Introduction In this thesis, we consider the problem of observing the state of dynamical systems, such as energy networks, as accurate as possible with a small number of sensors. This problem is interesting since there exist large networks, where it is not cost-efficient to place a sensor at every node of the system. Consequently, we can not measure all the states of the system which results in an observer error. The aim of this thesis is to find an optimal subset of sensors, i.e. a subset, such that every other subset of the same size has a bigger observer error. The problem of sensor selection appears in several areas of application like robotics, sensor placement for structures, target tracking, chemical plant control and wireless networks as listed in Joshi and Boyd [2009]. The field of application of this thesis are energy networks, especially high-voltage networks, as seen in the example of Section 4.4. We consider the problem of choosing an optimal subset from among n potential sensors of a time-invariant linear dynamic system for state estimation subject to white input noise. Each sensor can measure one component of the output vector y . Thus the process of sensor selection reduces the available information for the observer. The sensor selection consequently has an influence on the observer error. The aim of this thesis is to choose an H2-optimal subset, i.e. a subset that minimizes the H2 -norm of the observer error system. The H2 -norm describes the total output energy of the impulse response of the error system. A simple approach to evaluate the best k sensor subset would be to calculate the H2 -norm with an optimal observer for all possible n k combinations of subsets. But since n k grows rapidly with increasing n and k , this method is not practical. For example, with n = 50 potential sensors and k = 25 sensors to choose there are over 1014 possible tuples, so a sequential evaluation is obviously not possible. In this thesis, we provide several algorithms, partly based on convex optimization, for a near optimal solution of the sensor selection problem. And we use one combinatorial algorithm with a branch and bound technique for an optimal solution in order to compare the performance of the relaxed algorithms to this optimal solution. When comparing the execution time of the algorithms, it can be seen that the branch and bound algorithm is a NP-hard problem, while the other algorithms have a polynomial computational complexity. The problem we consider consists of two components. First we have to choose a subset of sensors which have to be used. Second, we have to calculate an optimal observer gain matrix for these settings in order to evaluate the objective function. In this approach, we provide a setting that solves these two problems in parallel. Various papers treat related problems to the one discussed here. In Mo et al. [2011], we find a general approach for sensor selection strategies for wireless sensor networks. However, the framework does not hold for our case since the strategy is chosen for a finite time horizon. This finiteness is used in the included manipulation of the optimization problem for a recursively defined equation of the state estimation. Given an infinite time horizon, we could not perform the same reformulations. 7 In Schuler et al. [2012], we find an almost dual problem to the observer design in this thesis. This paper with the title Decentralized State Feedback Control for Interconnected Process Systems seeks to minimize the number of measurement links between sensors and controllers and creates a convex optimization problem, that is similar to the problem derived in this thesis. Nevertheless, the derivation of the convex optimization problem is not applicable to the problem discussed here, because, on the one hand, the controller problem has slightly different bilinearities which have to be substituted for convex optimization, and on the other hand Schuler et al. [2012] use the H∞-norm in contrast to the here applied H2-norm. The approach in this thesis starts with describing the system of the estimation error. By introducing a design matrix, we create the possibility to vary the setting of the problem, i.e. to punish estimation errors in certain states more than in others. In the next step, we use a LMI-characterization of the H2 -norm to reformulate this constraint to a semidefinite one. The mentioned constraint is relaxed in two different ways that afterwards result in different algorithms. The `0 -norm for observer-sparsity in the objective is a non-convexity that is relaxed by an iterative weighted `1-norm. The rest of the thesis is organized as follows. We present the problem formulation in Chapter 2. In Chapter 3, we describe convex relaxations of the optimization problem and the derived algorithms both for discrete-time and continuous-time systems. In Chapter 4, extensive comparison of the different algorithms is presented, including a runtime and a performance comparison as well as an example of an energy-network model to test the algorithms in praxis. The thesis concludes with a summary and an outlook in Chapter 5. 8 2 Problem Formulation 2.1 Problem Formulation for Discrete-Time Systems Consider the following discrete-time dynamical system xk+1 =Axk+Bwk(2.1a) yk=Cxk, (2.1b) where k∈N+ 0 , xk∈Rn is the state of a system with n states and yk∈Rm , m≤n , is the set of possible sensors. The unknown input wk∈Rp is zero-mean white noise with unit variance. Assuming that the matrices A∈Rn×n , B∈Rn×p , and C∈Rm×n are known, we define a Luenberger observer. ˆxk+1 =Aˆxk+L(yk−ˆyk)(2.2a) ˆyk=Cˆxk, (2.2b) where L∈Rn×m describes the observer gain matrix that has to be designed. An optimal L could be found by using the principle of the Kalman filter, which uses the same setup as the one described here. We can express the observer error e=x−ˆx: ek+1 = (A−LC)ek+Bwk(2.3a) zk=Wek. (2.3b) The matrix W∈Rn×n can be used as a design matrix when the desired accuracy of the observations of certain states differs or if a scaling has to be done. Elsewise W could simply be defined as an identity matrix. Let us assume there exists a stabilizing L . Thus we can guarantee that the H2 -norm of the system exists and characterize the observer error with that norm, again as in the setup of Kalman filter. With the aim to minimize the number of used sensors with a bounded observer error we can define the following simplified optimization problem: min L”no. of sensors”(2.4a) s.t. ||Σe(L)||2 2< γ, (2.4b) where Σ e ( L )is a short form for the error system depending on L as described in (2.3). || ·||2 denotes the H2 -norm of the system and γ is a predefined upper bound to the squared H2 - norm. The number of used sensors in (2.4) is the number of non-zero columns in L , i.e. if every element of a column of L is zero, the corresponding entry in yk has no influence on the observer system, 9 Algorithm 1. Algorithm with substitution : Relaxed sensor placement algorithm with linearization via substitution for discrete-time systems 1. Set the iteration count µ to zero and choose an appropriate value for the parameter α . Initialize the weights vector ω(0) to ω(0) j = 1 for j = 1 , . . . , n and choose a sufficiently small ε > 0. 2. Solve the optimization problem (3.6). 3. Update the weights: ω(µ+1) j=1 ||˜ L∗j||1+ε,j= 1, . . . , n. 4. Terminate the iteration on convergence with L = X−1˜ L . Otherwise increase µ by 1 and return to Step 2. 5. Improve the choice of Lby solving the optimization problem: min γ, ˜ L,X γ(3.7a) s.t. X=XT0(3.7b) trace(BTXB)< γ (3.7c) −X+WTW ATX−˜ CT˜ LT XA −˜ L˜ C−X!≺0, (3.7d) where ˜ C := C with predefined zero lines. The j th line is set to zero if the j th sensor was removed in the previous solution of Step 2 (i.e. if ||L∗j||1ε). The weights vector ω(0) could also be initialized as a zero vector in Step 2 of the algorithm. Thus, the first iteration would yield to a optimal observer gain matrix with no zero columns in L and ω(1) would counteract the actual size of the column in an optimal non-sparse observer. With the current implementation it is possible that a sensor is set to zero in the first iteration because the column of L is smaller than the other columns. This is in general not equivalent to a bad performance of the sensor. In Step 3 of Algorithm 1, ˜ L or L could be used to perform the update on the weights. The reason why ˜ L is used here is that we want to counteract the size of ˜ L in the objective. As mentioned in Remark 1, it is equivalent to reduce either the non-zero columns of Lor ˜ L. In Step 3 of Algorithm 1, we introduced the variable ε to improve the reweighting with ω(µ) and avoid numerical problems, such as dividing by zero, as proposed in Candes et al. [2008]. Step 5 of Algorithm 1 is necessary, since the observer, that is calculated in the optimization problem (3.6), has a suboptimal observer gain matrix, because even the reweighted `1 -norm in the objective minimizes the absolute value of the entries of ˜ L . After Step 5, the exact H2 -norm of the observer error system can be calculated. The inequation ||L∗j||1ε in Step 5 is equivalent to the statement that the j th column is counteracted with a maximum weight. This means that the corresponding sensor is removed. Therefore, the jth sensor is not considered when computing the improvement of L. Remark 2. As written in equation (2.5), it is also possible to choose a constant γ and remove γ from the objective function. We choose the other way because it is on the one hand more 16 difficult to estimate an appropriate γ before starting the optimization, and on the other hand the algorithm has no incentive to choose the best subset of k sensors, if another subset of k sensors also satisfies the bound on the H2-norm. 3.3.2 Algorithm with Cone Complementarity Linearization As mentioned before, we will consider the matrix inequality (3.3) to perform a cone complementarity linearization. Here, you can see the matrix inequality again: −X+WTW AT−CTLT A−LC −X−1!≺0. (3.8) The following Lemma is necessary for formulating the algorithm with cone complementarity linearization. Lemma 3. With X0the following statements hold: (i) The following two expressions are equivalent: a) KX =I b) trace(KX) = nand K I I X!0 (ii) K I I X!0⇒trace(KX)≥n. Proof: Ad (i): Under condition of X0the following steps are valid: K I I X!0 ⇔Reformulation of the LMI with the Schur complement: K−IX−1I0 ⇔K−X−10. Postmultiplying with X leads to the next matrix inequality. ⇔KX −I0 With this reformulation the proof of "a) ⇒ b)" is obvious since KX = I implies trace( KX ) = trace(I) = nand KX −I0holds. For the proof of "b) ⇒ a)" we will look at a characteristic of the matrix inequality KX −I 0 and use the premise trace(KX) = trace(I) = n. trace(KX −I)=0is valid, because of trace(KX) = trace(I) = n. ⇒KX −I has all eigenvalues equal to zero. Because assuming that KX −I has a positive eigenvalue, which is possible since KX −I 0, would result in KX −I having a negative eigenvalue so that the sum of all eigenvalues is zero again. That leads to a contradiction to the statement KX −I0. Ad (ii): Under condition of X0the following step is valid: K I I X!0⇔KX −I0follows from the first steps of the proof of (i). Since the trace of a matrix is the sum of its eigenvalues, it results that trace( KX ) ≥trace ( I ) = n .  17 The lemma is vital for our approach to the second algorithm because we now have the possibility to introduce a new variable K that behaves as the inverse of X . The LMI mentioned in Lemma 3 can be used as a constraint, while trace( KX )can be part of the objective function. With minimizing trace( KX ), we will receive K as the inverse of X because trace( KX ) ≥n and K=X−1at the minimum, where trace(KX) = n. Applying Lemma 3 to the matrix inequality (3.8), the optimization problem (2.6) results in: min γ,L,X,K γ+α n X j=1 ω(µ) j||L∗j||1+βtrace(XK)(3.9a) s.t. X=XT0(3.9b) trace(BTXB)< γ (3.9c) −X+WTW AT−LTC A−LC −K!≺0(3.9d) K I I X!0, (3.9e) where β is another parameter for adjusting the objective function. As we can see, linearizing the matrix inequality resulted in a bilinear term in the objective function. Our aim is to have semidefinite program so we need to replace the bilinear part. We could use a linearization method as seen in El Ghaoui et al. [1997] to solve such a problem. At a given point ( Xold, Kold ), a linear approximation of trace(XK)is φlin(X, K) = constant +trace(XKold +XoldK). The optimization problem for the algorithm with cone complementarity linearization (CCL) is then: min γ,L,X,K γ+α n X j=1 ω(µ) j||L∗j||1+βtrace(XKold +XoldK)(3.10a) s.t. X=XT0(3.10b) trace(BTXB)< γ (3.10c) −X+WTW AT−LTC A−LC −K!≺0(3.10d) K I I X!0. (3.10e) Algorithm 2. Algorithm with CCL : Relaxed sensor placement algorithm with cone complementarity linearization for discrete-time systems 1. Find a feasible solution for X in (3.4) and set Xold = X . If there are none, exit. Set the outer iteration count µand the inner iteration count kto zero. 2. Set Kold = X−1 old and choose appropriate values for the parameters α , β , γ and δ . Initialize the weights vector ω(0) to ω(0) j= 1 for j= 1, . . . , n and choose a sufficiently small ε > 0. 3. Solve the optimization problem (3.10). 18 4. If |trace ( XK ) −n|< δ , set k = 0 and go to Step 5, else set k = k + 1, Xold = X and Kold =Kand go to Step 3. 5. Update the weights: ω(µ+1) j=1 ||L∗j||1+ε,j= 1, . . . , n. 6. Terminate the outer iteration on convergence. Otherwise increase µ by 1 and return to Step 3. 7. Improve the choice of L by solving the optimization problem (3.7). Define ˜ C := C with predefined constant zero lines. The j th line is set to zero if the j th sensor was removed in the previous solution of Step 3 (i.e. if ||L∗j||1ε). Remark 3. In Step 1 of Algorithm 2, there are several possibilities to find a feasible solution in (3.4). You could perform one iteration of (3.6) and use the generated X . In this case you could make the first update for the weights vector afterwards. Whereas a faster method is to find X by calculating a non-sparse Kalman filter and the corresponding Lyapunov matrix of the observer error system. But in this case the matrix X has to change much during the next iterations to remove sensors. This means that a lot of iterations are needed because the linearization method punishes any changes in X and K . Therefore, the algorithm removes sensors only very slowly. The objective function of this algorithm is divided up in three different parts, which are the bound on the H2 -norm γ , the `1 -relaxation of the number of sensors, and the punishing term for the inverse matrix of X . These three terms have totally different functions in the optimization problem and have a diverse scale for each term. Therefore, it is hard to find the right parameterization for the terms. Depending on the application of the algorithm, it is possible to reduce the complexity of the objective function. If only a reduced accuracy is needed for the result, one could take the term βtrace ( XKold + XoldK )from the objective function and use it as a constraint. With an appropriate constant ζthe constraint could look like: trace(XKold +XoldK)<2n+ζ. As mentioned in Remark 3, a very small trace ( XKold + XoldK )slows down the removal of sensors. Therefore, it could even be an advantage to put this term into the constraints and allow greater inaccuracy in order to have a result with less sensors. 3.4 Algorithm for Optimal Sensor Placement for Continuous-Time Systems In this section, we derive the approach for the algorithm with substitution for continuoustime systems. The algorithm is very similar to the one for discrete-time systems but the LMI characterization of the H2-norm is different since the Lyapunov equation is of another form. The following Lemma is like Lemma 2, but it is not necessary to show the second equivalence since the substitution even works when applied on the variables in the Lyapunov inequality. 19 Furthermore, in this thesis the cone complementarity linearization is used only for discrete-time systems as in the literature. Consequently, the augmented LMI for continuous-time systems is of no use for our approach. Lemma 4. (As in Rieber [2006]) Consider a continuous-time system with transfer function G(s) = "A B C0#:= C(sI −A)−1B. The following two statements are equivalent: (i) ||G(s)||2 2< γ and Ais asymptotically stable. (ii) There exists a positive definite symmetric solution P of the Lyapunov equation ATP + PA = −CTCand trace(BTPB)< γ. Proof: (partly as in Smith [2010]) ||G(s)||2 2< γ ⇔1 2πZ∞ −∞ trace(G(jω)G(jω)∗)dω < γ, which is the definition of the H2-norm on the frequency domain. In the next step we use an equivalent definition on the time domain. ⇔Z∞ 0 trace(g(t)Tg(t))dt < γ, where g(t)is the impulse response of the system G(s). ⇔trace[BT(Z∞ 0 eATtCTCAtdt)B]< γ ⇔trace[BTPB]< γ, where Pis the observability Gramian (possible, since Ais stable) which is defined as P=Z∞ t=0 eATtCTCeAtdt. ⇔trace[BTPB]< γ, P 0and ATP+PA =−CTC  We now consider the optimization problem (previously described in (2.6)) min L,γ γ+α n X j=1 || n X i=1 |Lij| ||`0(3.11a) s.t. ||Σe(L)||2 2< γ (3.11b) and apply the optimization problem to a continuous-time system (as described in (2.8)) ˙e= (A−LC)e+Bw (3.12a) z=We. (3.12b) 20 Using Lemma 4, the resulting optimization problem is the following: min γ,L,P γ+α n X j=1 ω(µ) j||L∗j||`1(3.13a) s.t. P=PT0(3.13b) trace(BTPB)< γ (3.13c) ATP+PA −CTLTP−PLC +I= 0, (3.13d) where P is the solution of the Lyapunov equation as in Lemma 4 (ii). The constraint (3.13d) could also be written as a LMI, if the matrix P is replaced by a matrix Xcont which is defined as the continuous-time equivalent of the matrix Xmentioned in the proof of Lemma 2. We could now perform the substitution ˜ L = PL as in Section 3.3.1. Consequently, Remark 1 holds also for this case and we replace L in the objective function with ˜ L . This results in the following optimization problem: min γ, ˜ L,P γ+α n X j=1 ω(µ) j||˜ L∗j||`1(3.14a) s.t. P=PT0(3.14b) trace(BTPB)< γ (3.14c) ATP+PA −CT˜ LT−˜ LC +I= 0. (3.14d) The algorithm with substitution is then formulated as: Algorithm 3. Algorithm with substitution : Relaxed sensor placement algorithm with linearization via substitution for continuous-time systems 1. Set the iteration count µ to zero and choose an appropriate value for the parameter α . Initialize the weights vector ω(0) to ω(0) j = 1 for j = 1 , . . . , n and choose a sufficiently small ε > 0. 2. Solve the optimization problem (3.14). 3. Update the weights: ω(µ+1) j=1 ||˜ L∗j||1+ε,j= 1, . . . , n. 4. Terminate the iteration on convergence with L = P−1˜ L . Otherwise increase µ by 1 and return to Step 2. 5. Improve the choice of Lby solving the optimization problem: min γ, ˜ L,P γ(3.15a) s.t. P=PT0(3.15b) trace(BTPB)< γ (3.15c) ATP+PA −˜ CT˜ LT−˜ L˜ C+I= 0, (3.15d) where ˜ C := C with predefined zero lines. The j th line is set to zero if the j th sensor was removed in the previous solution of Step 2 (i.e. if ||L∗j||1ε). 21 The comments concerning ˜ L , L , ε , and ω made about the algorithm with substitution for discrete-time systems are also applicable to this algorithm. 22 4 Comparison of Different Algorithms In this chapter, the algorithms of Chapter 3 are compared with state of the art algorithms. One is a greedy approach in two different versions as seen in Bach et al. [2010]. The other algorithm uses a branch and bound technique to provide an optimal solution which is also presented in the first part of this chapter. In the second part of the chapter, the algorithms are compared in terms of computation time and H2 -performance of the chosen subset of sensors. Additionally, the algorithms have to solve an application example consisting of a model of a power grid. 4.1 State of the art Algorithms In the following, we discuss three algorithms that use the original optimization problem without any relaxation. Hence, these algorithms cannot use interior point methods as the two algorithms mentioned before. These algorithms do not optimize the chosen subset of sensors and the observer gain matrix in parallel. In fact, the algorithms compute the H2 -norm of different subsets and compare this H2 -performance. Hence, there is no change in difficulty of the computation when the algorithms are implemented for discrete-time or continuous-time models. 4.1.1 Greedy Algorithm The first state of the art algorithm, which is described, is a greedy approach that is very intuitive. The algorithm starts with all sensors and compares all subsets with n− 1sensors. After turning-off the worst sensor the algorithm consequently searches for the worst sensor of the remaining subset. Consequently, the greedy algorithm finds subsets of every size. This greedy algorithm generates a fast solution and is able to give an overview of how many sensors are necessary for which H2 performance. But as the algorithm removes all the sensors separately, it could not account for correlation between certain sensor subsets. This could lead to a very poor performance. The algorithm searches for a subset of k sensors, where k is a predefined number with k < n . This means that the algorithm has to remove n−ksensors. Define: J ( k1, k2, . . . , km )is the H2 -norm of the error system with an optimal reduced observer which uses the sensors k1, k2, . . . , km with k1, k2, . . . , km∈ { 1 , 2 , . . . , n} . The evaluation of this function could be done with the predefined system matrix C , in which one replaces every column except the columns k1, . . . , km with zero columns in order to remove the corresponding sensors. One way to evaluate the function is to solve the optimization problem (4.1) that is presented 23 below. min γ, ˜ L,X γ(4.1a) s.t. X=XT0(4.1b) trace(BTXB)< γ (4.1c) −X+WTW ATX−CT˜ LT XA −˜ LC −X!≺0(4.1d) Algorithm 4. Greedy algorithm : Greedy approach to the sensor selection problem as in Bach et al. [2010]. 1. Initialize M:= {1,2, . . . , n}and iteration variable i:= 1. 2. Compute m=arg min lJ(M\ {l})with l∈M. 3. Define M:= M\m. 4. If i = n−k , end the algorithm with the solution M . If not, increase i by one and go back to Step 2. Mis the determined subset of sensors and J(M)the corresponding value of the H2norm. The greedy algorithm could also be implemented backwards. This means, the algorithm starts with zero sensors and then iteratively adds the best sensor. This implementation of the algorithm is called reverse greedy algorithm in the following. 4.1.2 Branch and Bound Algorithm The next algorithm is a combinatorial algorithm, which finds the optimal subset of k sensors, where k is a predefined number with k < n . The algorithm is called branch and bound algorithm, because it arranges the sensors in branches and looks for an upper bound at the beginning of the algorithm. If the algorithm finds a good upper bound, it is possible that a lot of branches, i.e. tuples of sensors, will not have to be evaluated. The exact algorithm as it was implemented during the writing of this thesis and more details about the algorithm can be found in Narendra and Fukunaga [1977]. 4.2 Runtime Comparison For comparison discrete-time dynamical systems are generated randomly with B, C, and, W as identity matrices. The dynamic matrix A is chosen randomly with each entry uniformly, independently, and identically distributed on [0 , 1] and 80% zero elements to have a sparse structure as it is common in praxis for example when considering power grids. The continuous-time systems were generated by converting the randomly generated discrete-time systems to a continuous-time model. The loss of the sparse structure during the converting process was not considered in this comparison. Algorithms 1, 2 and 3 have objective functions in which the number of chosen sensors is controlled only indirectly by chosing the parameters α and β . For the test, the parameters are set so that 24 about one fourth of the sensors were removed. However for these algorithms, the paramterization has little influence on the execution time if we assume convergence of the algorithms. Although a bad set of parameters could easily result in divergence, especially when handling the algorithm with CCL, and thus to no useful solution. The branch and bound algorithm and the greedy algorithms had to choose subsets of n 2 sensors. In contrast to the other algorithms, the execution time is dependent on the number of chosen sensors, since both start with all sensors and consequently remove sensors. When choosing very few sensors, the branch and bound algorithm has a very large execution time. The algorithm could even be outperformed by an algorithm, which tries all possible tuples, because the upper bound used in the branch and bound algorithm would be too big and therefore nearly useless. The formulation of the optimization problem for the algorithms with substitution and the algorithm with CCL is a semidefinite programming problem, since the constraints are linear matrix inequalities (see also in Vandenberghe and Boyd [1999]) and the objective is a convex function and can even be recast as a linear function. This allows us to use interior point methods which have polynomial complexity as can be seen in Sturm [1999]. The reweighted `1 relaxation usually needs only few steps to converge when the parameter ε is chosen appropriate (as in Candes et al. [2008]). In the (reverse) greedy algorithm, the number of computations of an optimal observer gain matrix and the corresponding H2 -norm increases with n2 , where n is the number of states, assuming that it has to reduce a certain percentage of the senors. The optimal observer gain matrix could be evaluated as a Kalman filter, which has the complexity of n3 . So the greedy algorithm also has polynomial complexity. The reverse greedy algorithm computes the optimal observer gain matrix as often as the greedy algorithm when choosing a subset of n 2 sensors. Hence, the execution time of these two greedy algorithms is the same. In the following, we do not distinguish between these two algorithms when considering the runtime. The branch and bound algorithm is NP-hard because at a worst-case behavior it has to compute more than all possible subsets of k sensors, which are n k . For this reason the branch and bound algorithm was mainly used here to have an optimal subset in order to compare the performance of the other algorithms. This comparison can be seen in the next section. For each size, we generated 10 random examplary systems as described above. Every algorithm had to find the optimal sensors for all these systems. In Figure 4.1 we see the execution time in seconds. The displayed execution time is the average of the different computations. The dashed lines show the greedy algorithm and the branch and bound algorithm with a slightly changed implementation. While the solid lines were produced with algorithms that used an implementation of the kalman function in Matlab, the dashed lines solved the optimization problem (4.1), which was implemented with SeDuMi (Sturm [1999]) and YALMIP (Lofberg [2004]). Both of the implementations achieve the same results, but the implementation with SeDuMi and YALMIP needs more execution time and allows a better comparison with the other three algorithms that also use SeDuMi and YALMIP. As we can see, the greedy algorithm, which uses the kalman function method of Matlab, is the fastest algorithm. The algorithm with substitution and the algorithm with CCL are very similar, but the algorithm with CCL needs more execution time, which could be the effect of the second loop in the algorithm. Moreover, the algorithm with CCL quite often diverged which was not considered for the execution time here. 25 The algorithm with cone complementarity linearization has the same benefits as the algorithm with substitution, except for the good execution-time and the easy paramterization. The algorithm with CCL has a second loop inside of the iterative reweighted `1 -minimization which slows the algorithm down. The great disadvantage of this algorithm is that it is very hard to parameterize the algorithm correctly. In the application example no appropriate parameterization was found although over 50 different sets of parameters were tried. Nevertheless, the algorithm only lead to useless results. The branch and bound algorithm generates the optimal solution and is easy to parameterize. In contrast to the other algorithms, the branch and bound algorithm has a non-polynomial complexity and therefore an execution-time that is not usable in praxis. The greedy algorithm is an intuitive and fast approach which generates usually a good solution. However, the shortsighted approach cannot consider the importance of certain sensors in small subsets. Another application example was generated where the greedy algorithm first removes the sensor at the highest branched node, although this sensor would belong to the optimal subset at the end. It seems that highly branched nodes are easy to observe with a lot active sensors so the sensor at this node is removed at the beginning of the algorithm. However, with few active sensors, the highly branched sensors are important to observe the states at other nodes. It is possible that this property is more important when considering larger networks and did not affect the results in Section 4.3. However, in our examples the greedy algorithm shows a good performance and generates an overview of the H2-performance for different sizes of subsets. The reverse greedy algorithm is also a fast algorithm with a good solution to our examples. In contrast to our examples, more complex grids with instable parts can not be solved by the reverse greedy algorithm because it needs an observer with only one sensor that has a stable observer error system. If such a sensor could not be found, the algorithm has no starting point and ends with no result. As the other greedy approach, the reverse greedy algorithm can only perform shortsighted choices and cannot consider correlation of the sensors. Therefore, it is possible that a very poor solution is provided. When solving a sensor placement problem in praxis, the algorithm with substitution and the (reverse) greedy algorithm should be chosen. The (reverse) greedy algorithm can provide a fast overview of the performance of the desired number of sensors and a first subset. After that, the algorithm with substitution can be used to improve the chosen subset or confirm its performance. 32 5 Conclusion and Future Directions In this thesis, we considered a H2 -optimal sensor placement for linear dynamical systems for an infinite time horizon. The model of the system can either be a discrete-time or a continuous-time model. Furthermore, the influence of errors on certain states of the system can be parameterized as well as the subset of all possible sensors. The problem of sensor placement arises in various fields of application, for example at observing power grids as described in chapter 4. The aim is here to reduce costs for expensive sensors and still have the best possible observer system for a given number of sensors. The original problem of sensor placement is a non-convex problem, which can either be approached by combinatorial algorithms like the branch and bound algorithm or with greedy algorithms in different variations, as described in Chapter 4. In contrast to these methods, the main algorithms in this thesis use another approach. The presented algorithms optimize the structure, i.e. the sensor placement, and the observer gain matrix in parallel. Due to non-convexities, which arise when counting the used sensors and when calculating the H2 -norm with an unknown observer, we had to reformulate the optimization problem. The discrete `0 -norm was relaxed with an iterative reweighted `1 -norm. The H2 -norm was written with a LMI characterization, which simplified replacing the non-convex bilinear part either with a substitution or with a cone complementarity linearization. The proposed approaches were then developed into algorithms and implemented to efficiently and accurately solve the sensor placement problem. The mentioned algorithms were tested in Chapter 4. The first test was a comparison of the execution time depending on the system size. The second test checked on the H2 -performance of the chosen sensor subsets of the different algorithms, while comparing the ability to parameterize the algorithms so that different numbers of sensors were selected. Additionally, a practical example was solved. Here, the algorithms had to determine an optimal sensor subset. The performance of the chosen subsets were compared. The example had slightly other requirements for the algorithms since the problem was more ill-conditioned than the randomly generated examples. This revealed further difficulties when parameterizing the algorithm with cone complementarity linearization. The conclusion of all the tests is that the algorithm with substitution is one of the best of the presented algorithms, since it has no severe lacks, e.g. execution time or parameterization, and a satisfying performance. Similar to the greedy algorithm, the algorithm with substitution for continuous-time systems has a very short execution time despite of its not advanced implementation. The greedy approach and the reverse greedy approach both lead to good solutions at the considered examples and are the fastest algorithms. The shortsighted procedure does not influence the performance at the described examples. Future research concerning this topic could expand this approach into the following directions. The assumption that the system matrix Ais completely known may not be fulfilled. It is more realistic to assume that one knows A with certain inaccuracies. Another possibility is that a non-linear system can be estimated in a defined sector and that the algorithms use bounds of this sector to derive an optimal subset of sensors. 33 Another direction of future research could be the expansion of the model to actuator placement which is the dual problem to the one described here. Moreover, for actuator placement usually the H∞ -norm is used. Consequently, the approach would have to be adapted to the LMI characterization of this norm. 34 6 Appendix - Mathematical Definitions and Notation 6.1 `0-Vector-Norm The `0 -norm for vectors has more than one commonly used definition like ||x||0 := limp→0||xi||p p for p∈R+ and x∈Rn . Precisely said, the `0 -norm refers to the number of non-zero elements in a vector. Therefore it is straightforward to define also: ||x||0:= n X i=1 |sign(xi)| with sign(xi) :=      −1for xi<0 0for xi= 0 1for xi>0 Since the positive scalability (i.e. |α| ||x||0 = ||αx||0 for α∈R+ ) is not satisfied in general, the `0-norm is in fact not a real norm, however it is often referred to as the 0-norm. The `0-norm can be relaxed using a `1-norm, which is defined as follows: ||x||1:= n X i=1 |xi|. 6.2 H2-System-Norm Define G ( s ) = "A B C0# := C ( sI −A ) −1B as the state-space realization of the transfer matrix, QC := R∞ 0eAtBBTeATtdt as the controllability Gramian and QO := R∞ 0eATtCCTeAtdt as the observability Gramian for continuous-time systems. The H2-norm of a system G(s)is defined as (according to Scherer and Weiland [2000]): ||G(s)||2:= s1 2πZ∞ −∞ trace(G(jω)G(jω)∗)dω. (6.1) This system norm is used to specify the performance of the systems in this thesis, because it is the typical norm for observers of dynamical systems, e.g. the well-established Kalman filter (see Herrnberger [2010]). A typical interpretation of the norm is to see ||G ( s ) ||2 as the total output energy of an impulse response of the system(as in Scherer and Weiland [2000]). 35 The relation of the H2-norm to the Gramians for a system G(s) = "A B C0#is ||G||2 2=trace(CQCCT) = trace(BTQOB). The definition of the H2 -norm for discrete-time systems is very similar. Define G ( z ) = "A B C0# := C ( zI −A ) −1B as the state-space realization of the transfer matrix for a discrete-time system. Consequently, the H2-norm for discrete-time systems is defined as: ||G(z)||2:= s1 2πZπ −π trace(G(ejω)G(ejω)∗)dω. (6.2) 6.3 Matrices and LMIs Next, we will give some notational specifications of this thesis. Given a matrix L∈Rn×n we will denote L−1 and LT as the inverse and the transpose of the matrix. L∗j with j∈1,2, . . . , m defines the jth column of the matrix L. A positive or negative definite (semidefinite) matrix L is written as L 0( L 0) or L≺ 0 (L0) respectively, while >, <, ≤,≥denote an elementwise comparison. 36 Bibliography F. Bach, S.D. Ahipasaoglu, and A. d’Aspremont. Convex relaxations for subset selection. ArXiv preprint ArXiv:1006.3601, 2010. E.J. Candes, M.B. Wakin, and S.P. Boyd. Enhancing sparsity by reweighted `1 minimization. Journal of Fourier Analysis and Applications, 14(5):877–905, 2008. L. El Ghaoui, F. Oustry, and M. AitRami. A cone complementarity linearization algorithm for static output-feedback and related problems. Automatic Control, IEEE Transactions on, 42 (8):1171–1176, 1997. M. Fazel. Matrix rank minimization with applications. PhD thesis, Stanford University, 2002. M. Herrnberger. Moderne Methoden der Regelungstechnik III. Lecture at the Institute of Automatic Control in SS 2010, Technichal University Munich, Germany, October 2010. S. Joshi and S. Boyd. Sensor selection via convex optimization. Signal Processing, IEEE Transactions on, 57(2):451–462, 2009. P. Kundur, N.J. Balu, and M.G. Lauby. Power system stability and control, volume 4. McGraw-hill New York, 1994. J. Lofberg. Yalmip: A toolbox for modeling and optimization in matlab. In Computer Aided Control Systems Design, 2004 IEEE International Symposium on, pages 284–289. IEEE, 2004. Y. Mo, R. Ambrosino, and B. Sinopoli. Sensor selection strategies for state estimation in energy constrained wireless sensor networks. Automatica, pages 1330–1338, 2011. P.M. Narendra and K. Fukunaga. A branch and bound algorithm for feature subset selection. Computers, IEEE Transactions on, 100(9):917–922, 1977. J. Rieber. Lecture robust contol. Lecture at the Institute for Systems Theory and Automatic Control in WS 2006/2007, University of Stuttgart, Germany, October 2006. C. Scherer and S. Weiland. Linear matrix inequalities in control. Lecture Notes, Dutch Institute for Systems and Control, Delft, The Netherlands, 2000. S. Schuler, U. Münz, and F Allgöwer. Decentralized state feedback control for interconnected process systems. In AdChem, 2012. R. Smith. Robust control & convex optimization. Lecture at the department of electrical & computer engineering in SS 2010, University of California, Santa Barbara, USA, June 2010. J.F. Sturm. Using sedumi 1.02, a matlab toolbox for optimization over symmetric cones. Optimization methods and software, 11(1-4):625–653, 1999. L. Vandenberghe and S. Boyd. Applications of semidefinite programming. Applied Numerical Mathematics, 29(3):283–299, 1999. 37 Danksagung An dieser Stelle möchte ich mich recht herzlich bei Herrn Dr . Ulrich Münz für die intensive Betreuung bedanken, die teilweise viel Zeit in Anspruch nahm. Ohne die zahlreichen fachlichen Diskussionen und interessanten Anregungen wäre diese Arbeit so nicht zustande gekommen. Außerdem bedanke ich mich auch bei der Siemens AG für die Ermöglichung der Arbeit und die finanzielle Unterstützung währenddessen. Herrn Prof . Dr . Tobias Damm danke ich für die sehr angenehme Kooperation und die sehr hilfsbereite Beantwortung bei allen fachlichen oder organisatorischen Fragestellungen. Schließlich bedanke ich mich noch bei Herrn Prof . Dr . Lars Grüne für die Bereitschaft, die Arbeit des Zweitgutachters zu übernehmen. Erklärung Hiermit erkläre ich, dass ich diese Bachelorarbeit selbstständig verfasst und keine anderen als die von mir angegebenen Quellen und Hilfsmittel benutzt habe und dass ich diese Arbeit nicht bereits zur Erlangung eines akademischen Grades eingereicht habe. Bayreuth, 17. August 2012