Generalized Matrix-Vector Interface =================================== The matrix-vector interface described earlier is suitable for problems where the independent and dependent variables are both associated with buses. This is usually true for systems where Kirchhoff’s law is strictly enforced. However, this interface does not work for systems where some variables are associated with branches. This can occur in optimization problems such as state estimation, where measurements are made on both buses and branches. Every measurement contributes an equation to the state-estimation optimization, which results in dependent variables associated with branches. To handle these types of problems, a more general approach to creating matrices and vectors is required. This is implemented via the ``GenMatVecInterface`` class described below. In general, the ``GenMatVecInterface`` is much more complicated to use than the matrix-vector interface described in Section `4.7 <#mapper>`__. However, we believe that it is still simpler than managing the index calculations needed to construct the matrix equations required by methods such as state estimation and dynamic state estimation (Kalman filter). User that are interested in implementing calculations using the ``GenMatVecInterface`` are encouraged to follow the documentation in the user manual and to also study the implementation in existing GridPACK applications. .. _gen_matvec: Generalized Matrix-Vector Mappers --------------------------------- Unlike the ``MatVecInterface`` class, there is no definitive way to map which elements are contributed by a branch or bus, and the number of elements contributed by a branch or bus does not reduce to simple blocks. Thus, the idea that buses and branches contribute simple blocks of data must be abandoned. The ``GenMatVecInterface`` just assumes that buses and branches contribute some number of equations (dependent variables) to the matrix and that they also contribute some number of independent variables to the matrix. As illustrated in Figure 5, the ``BaseComponent`` class directly inherits from this interface, along with the ``MatVecInterface``. The information on number of equations and number of variables is embedded in the function calls :: virtual int matrixNumRows(void) virtual int matrixNumCols(void) These two functions specify how many dependent variables (rows) and how many independent variables (columns) are associated with a bus or branch. For the state estimation module that is currently available in the GridPACK release, the dependent variables are the number of measurements that are associated with the bus or branch and the independent variables are the voltage magnitude and phase angle, which are only associated with buses. Thus, if the state estimation Jacobian is being built, the ``matrixNumRows`` function returns the number of measurements on each bus and branch. The ``matrixNumCols`` only returns a non-zero value for buses since the branches have no independent variables. This value is generally 2, if the bus has any measurements associated with it or is attached to a bus or branch that has measurements, otherwise the value is 0. If the bus has measurements and is the reference bus, then the function returns 1. These functions allow the generalized mappers to determine the dimensions of the matrix (for state estimation, the Jacobian is not necessarily square). Unlike the original matrix-vector interface, the user has to assign the row and column indices to each matrix element. The actual values of these indices are evaluated by the mapper but it is up to the user to take the row index for a particular dependent variable (measurement) and the column index for a particular independent variable (voltage magnitude or phase angle) and pair them with a matrix element (contribution to the Jacobian). The functions that are used for this purpose are :: virtual void matrixSetRowIndex(int irow, int idx) virtual void matrixSetColIndex(int icol, int idx) virtual int matrixGetRowIndex(int irow) virtual int matrixGetColIndex(int icol) The first two functions are used by the mapper to assign indices for each of the rows and columns contributed by a component. The values of the indices need to be stored in the component so that they can be accessed by other components when evaluating matrix elements. Although these functions are only called by the mapper, they need to be implemented by the user, since multiple matrices may be generated by the application. The variables ``irow`` and ``icol`` refer to the list of rows and columns contributed by the component, while the index ``idx`` is the global index for that row or column in the full matrix. The point of the first two functions is to create a map between the local index of the row or column and the global index of the corresponding row or column in the full matrix. This map is needed because matrix elements constructed on one component may refer to rows or columns on other components. The second pair of functions allow users to recover the global index from the local index. For example, the state estimation calculation needs to be able to build the Jacobian matrix plus a diagonal matrix that represents the inverse of the uncertainties in all the measurements. The state estimation components have two modes, ``Jacobian_H`` and ``R_inv`` for each of these calculations. The ``matrixSetRowIndex`` method for the buses has the form :: void SEBus::matrixSetRowIndex(int irow, int idx) { if (p_mode == Jacobian_H) { if (irow < p_rowJidx.size()) { p_rowJidx[irow] = idx; } else { p_rowJidx.push_back(idx); } } else if (p_mode == R_inv) { if (irow < p_rowRidx.size()) { p_rowRidx[irow] = idx; } else { p_rowRidx.push_back(idx); } } } The row indices for the Jacobian and :math:`R^{-1}` are stored in two separate STL arrays ``p_rowJidx`` and ``p_rowRidx``. For the state estimation example, the number of rows (for both the Jacobian and :math:`R^{-1}`) is equal to the number of measurements associated with the component. These measurements are held in an internal list in some order. If the number of measurements on the bus is M then the ``irow`` index will run from 0,..,M-1, with the ``irow`` index corresponding to the ``irow`` element in the list of measurements. The independent variables are also assumed to be ordered in some fashion. Again, for the state estimation example, the phase angle is indexed by 0 and the voltage magnitude is indexed by 1. The function for accessing the row indices is implemented as :: int gridpack::state_estimation::SEBus::matrixGetRowIndex(int idx) { if (p_mode == Jacobian_H) { return p_rowJidx[idx]; } else if (p_mode == R_inv) { return p_rowRidx[idx]; } } Again, depending on the mode, this function will return different values and for this reason, these functions need to be implemented by the user. They cannot be implemented as part of the framework because the number of modes is application-specific and controlled by the developer. However, the form of these functions is fairly standard and can be generalized in a straightforward way from the examples presented here. The functions that are used to actually evaluate matrix elements are :: virtual int matrixNumValues(void) const virtual void matrixGetValues(ComplexType *values, int *rows, int *cols) virtual void matrixGetValues(double *values, int *rows, int *cols) The first function returns the total number of matrix elements that will be evaluated by the component. This is used inside the mapper to allocate arrays that hold matrix elements coming from the components. The second function is used to evaluate actual matrix elements, along with their row and column indices. The real-valued version of ``matrixGetValues`` replaces ``ComplexType`` with ``double``. The ``matrixGetRowIndex`` and ``matrixGetColIndex`` functions are used in implementing the ``matrixGetValues`` function. The evaluation of the ``matrixNumValues`` function can be quite complicated. For the state estimation Jacobian matrix, the number of matrix elements contributed by a component depends on the number of measurements associated with that component and the number of variables that couple to that measurement. A measurement on a bus will usually contribute two values for the independent variables on the bus, plus an additional two values for each bus that is attached to the center bus via a branch. This number will be modified slightly if one of the buses in this group is a reference bus. For branches, the number of matrix elements contributed by each measurement is approximately four, two elements for each bus at either end of the branch. This number may drop if one of the buses is a reference bus. The ``matrixGetValues`` function is used to evaluate each of the matrix elements. It also gets the matrix indices for this element from the appropriate network component. The number of matrix elements returned by this function must correspond to the number returned by the ``matrixNumValues`` function. To see how the assignment of the indices works, we can look at the matrix element of the Jacobian corresponding to the gradient of a real power injection measurement *P\ :math:`{}_{i}`* on bus *i* with respect to the phase angle on another bus *j* that is connected to *i* via a single branch. The contribution to the Jacobian from this measurement is given by the formula .. math:: \frac{\partial P_i}{\partial {\theta }_j}=V_iV_j(G_{ij}{\mathrm{sin} \left({\theta }_i-{\theta }_j\right)-B_{ij}{\mathrm{cos} ({\theta }_i-{\theta }_j)\ }\ }) Suppose :math:`P_{i}` is measurement :math:`k` on the bus. Then the row index ``im`` for this matrix element can be evaluated by calling the function :: im = matrixGetRowIndex(k); The column index is associated with the phase angle variable on the remote bus *j*. Assuming that a pointer (``bus_j``) to the remote bus is already available, then the column index ``jm`` for this matrix element could be obtained by calling :: jm = bus_j->matrixGetColIndex(0); This function is called with the argument 0 since the dependent variables are always ordered as phase angle (0) followed by voltage magnitude (1). The full list of Jacobian matrix elements can be obtained by looping over all measurements. For each bus measurement, there are contributions from the dependent variables on each connected bus plus two contributions from the calling bus. Similarly, for each branch measurement there are approximately four contributions coming from the independent variables associated with the buses at each end of the branch. A simple counter variable can be used to make sure that the matrix element value and the corresponding row and column indices stored in the same location of the ``values``, ``rows``\ ````\ and ``cols`` arrays that are returned by the ``getMatrixValues`` function. The ``GenMatVecInterface`` also includes functions for setting up vectors. These work in a very similar way to the generalized matrix functions, so they will only be described briefly. The two functions :: virtual void vectorSetElementIndex(int ielem, int idx) virtual void vectorGetElementIndices(int *idx) can be used to set and retrieve vector indices. The index ``ielem`` is the local index within the element while ``idx`` is the global index within the distributed vector. In this case it is usually more convenient to get all indices associate with a component at once, so the ``vectorGetElementIndices`` returns an array instead of a single value. The function :: virtual int vectorNumElements() const returns the number of vector elements contributed by a component and the function :: virtual void vectorGetElementValues(ComplexType *values, int *idx) virtual void vectorGetElementValues(double *values, int *idx) returns a list of the values along with their global indices. For real vectors, the ``ComplexType`` array is replaced an array of type ``double``. Again, the index value can be obtained by first calling the ``vectorGetElementIndices`` function and using this to obtain the correct index for each element. The vector interface includes one additional function that does not have a counterpart in the matrix interface. This is the function :: virtual void vectorSetElementValues(ComplexType *values) virtual void vectorSetElementValues(double *values) This function can be used to push values from a solution vector back into the network components. The values are ordered in the same way as the values in the corresponding ``vectorGetElementValues`` call, so it is possible to unpack them and assign them to the correct internal variables for each component. This function is analogous to the ``setValue``\ s call in the regular ``MatVecInterface``. The functions in the ``GenMatVecInterface`` are invoked in the generalized mappers. These reside in the ``GenMatrixMap`` and ``GenVectorMap`` classes. Like the standard mappers, these classes are relatively simple and contain only a few methods. The ``GenMatrixMap`` class consists of the constructor :: GenMatrixMap(boost::shared_ptr network) and the methods :: boost::shared_ptr mapToMatrix(void) boost::shared_ptr mapToRealMatrix(void) void mapToMatrix(boost::shared_ptr matrix) void mapToRealMatrix(boost::shared_ptr matrix) void mapToMatrix(gridpack::math::Matrix &matrix) void mapToRealMatrix(gridpack::math::RealMatrix &matrix) void overwriteMatrix(boost::shared_ptr matrix) void overwriteRealMatrix(boost::shared_ptr matrix) void overwriteMatrix(gridpack::math::Matrix &matrix) void overwriteRealMatrix(gridpack::math::RealMatrix &matrix) void incrementMatrix(boost::shared_ptr matrix) void incrementRealMatrix(boost::shared_ptr matrix) void incrementMatrix(gridpack::math::Matrix &matrix) void incrementRealMatrix(gridpack::math::RealMatrix &matrix) These functions all have the same behaviors as the analogous functions in the standard ``FullMatrixMap``. The ``GenVectorMap`` class has the constructor :: GenVectorMap(boost::shared_ptr network) and supports the methods :: boost::shared_ptr mapToVector(void) boost::shared_ptr mapToRealVector(void) void mapToVector(boost::shared_ptr &vector) void mapToRealVector(boost::shared_ptr &vector) void mapToVector(gridpack::math::Vector &vector) void mapToRealVector(gridpack::math::RealVector &vector) These functions have the same interpretations as the analogous functions in the ``BusVectorMap``\ ````\ class. This interface also hosts a new function :: mapToNetwork(boost::shared_ptr &vector) mapToNetwork(boost::shared_ptr &vector) mapToNetwork(const gridpack::math::Vector &vector) mapToNetwork(const gridpack::math::RealVector> &vector) which can be used to push data from a vector back into the network components (both buses and branches). Generalized Slab Mapper ----------------------- The generalized slab mapper also uses functions in the generalized matrix-vector interface to build dense matrices. These matrices are dense since they are generated by taking a typical vector, corresponding to a set of variables on the buses and branches, and inserting copies of the vector into the matrix for different values of the variables. An example would be a matrix formed from a time series of values for a set of variables on the buses and branches. One set of indices for the matrix corresponds to the set of variables and the other set of indices corresponds to the time series. In a certain sense, these matrices are “fat” vectors since instead of each variable having only one value, they have multiple values. In general, slab matrices are not square. The slab matrices are used in the Kalman filter application, but they may have applicability elsewhere. The slab mappers use additional functions from the ``GenMatVecInterface`` in order to construct matrices. These functions are analogous to the functions for setting up vectors using the ``GenVectorMap``. The main difference is that instead of describing a list of values, the functions describe a matrix block. The row dimension corresponds to a list of variables and the column dimension describes the number of values taken by each variable. The column dimension must be the same across all variables. The contribution to the matrix from each network component is given by the function :: void slabSize(int *rows, int *cols) const The index for each row can be stored using the function :: void slabSetRowIndex(int irow, int idx) This function is called by the mapper and is analogous to the ``vectorSetElementIndex`` function. For the slab matrices, there is no corresponding call for columns since the matrices are dense and all rows have the same number of non-zero columns. The indices can be retrieved by the function :: void slabGetRowIndices(int *idx) which is similar to the ``vectorGetElementIndices`` function.