How can OpenMDAO be used to solve a linear system of equations without inverting the A matrix?

Question

I have a system of equations that is in the form:

Ax = b

Where A and b are a mixture of known states and state rates derived from earlier components and x is a vector of four yet unknown state rates. I've used Matlab to linearise the problem, all I need to do now is to create some components to find x. However, the inverse of A is large in terms of the number of variables in each index, so I can't just turn these into a straightforward linear equation. Could someone suggest a route to go?

Answer 1

I don't fully understand what you mean by "the inverse of A is large in terms of the number of variables in each index", however I think mean that the inverse of A is to larger and dense to compute and store in memory.

OpenMDAO or not, When you run into this situation you are forced to use an iterative linear solver such as gmres. So that is broadly the approach that is needed here too.

OpenMDAO does have a LinearSystemComponent that you can use as a rough blueprint here. However, it does compute a factorization and store it which is not what you want. Regardless, it gives you the blueprint for how to represent a linear system as an implicit component in OpenMDAO.

Broadly, you have to think of defining a linear residual: R = Ax-b = 0

Your component will have two inputs A and b, and and one output x.

The two key methods here are apply_nonlinear and solve_nonlinear. I realize that the word nonlinear in the method names is confusing. OpenMDAO assumes that the analysis is nonlinear. In your case it happens to be linear, but you use the nonlinear methods all the same.

I will assume that, although you can't compute/store [A] inverse you can compute/store A (perhaps in a sparse format). In that case you might pass the sparse data array of [A] as the input and fill the sparse matrix as needed from that.

the apply_nonlinear method would look like this:

def apply_nonlinear(self, inputs, outputs, residuals):
        """
        R = Ax - b.
        Parameters
        ----------
        inputs : Vector
            unscaled, dimensional input variables read via inputs[key]
        outputs : Vector
            unscaled, dimensional output variables read via outputs[key]
        residuals : Vector
            unscaled, dimensional residuals written to via residuals[key]
        """
        
        residuals['x'] = inputs['A'].dot(outputs['x']) - inputs['b']

The key to your question is really the solve_nonlinear method. It would look something like this (using scipy gmres):

def solve_nonlinear(self, inputs, outputs):
        """
        Use numpy to solve Ax=b for x.
        Parameters
        ----------
        inputs : Vector
            unscaled, dimensional input variables read via inputs[key]
        outputs : Vector
            unscaled, dimensional output variables read via outputs[key]
        """
    
        
        x, exitCode = gmres(inputs['A'], inputs['b'])

        outputs['x'] = x