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. 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.

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 \(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 \(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

\[\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 \(P_{i}\) is measurement \(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 setValues 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<MyNetwork>(boost::shared_ptr<MyNetwork> network)

and the methods

boost::shared_ptr<gridpack::math::Matrix> mapToMatrix(void)
boost::shared_ptr<gridpack::math::RealMatrix> mapToRealMatrix(void)

void mapToMatrix(boost::shared_ptr<gridpack::math::Matrix> matrix)
void mapToRealMatrix(boost::shared_ptr<gridpack::math::RealMatrix> matrix)

void mapToMatrix(gridpack::math::Matrix &matrix)
void mapToRealMatrix(gridpack::math::RealMatrix &matrix)

void overwriteMatrix(boost::shared_ptr<gridpack::math::Matrix> matrix)
void overwriteRealMatrix(boost::shared_ptr<gridpack::math::RealMatrix> matrix)

void overwriteMatrix(gridpack::math::Matrix &matrix)
void overwriteRealMatrix(gridpack::math::RealMatrix &matrix)

void incrementMatrix(boost::shared_ptr<gridpack::math::Matrix> matrix)
void incrementRealMatrix(boost::shared_ptr<gridpack::math::RealMatrix> 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<MyNetwork>(boost::shared_ptr<MyNetwork> network)

and supports the methods

boost::shared_ptr<gridpack::math::Vector> mapToVector(void)
boost::shared_ptr<gridpack::math::RealVector> mapToRealVector(void)

void mapToVector(boost::shared_ptr<gridpack::math::Vector> &vector)
void mapToRealVector(boost::shared_ptr<gridpack::math::RealVector> &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<gridpack::math::Vector> &vector)
mapToNetwork(boost::shared_ptr<gridpack::math::RealVector> &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.