OpenMDAO implicit relationship between components

Question

I have the following system, that I am trying to solve in OpenMDAO, as a part of a larger system:

C*x = F1(x) + F2(x)

where F1 and F2 are some responses (6x1 vectors) computed for the given x (6x1 vector) in separate explicit components (F1, F2 being the outputs of the components, and x being the input in both components). The relationship between F and x cannot be described analytically. C does not depend on x, and is a 6x6 matrix the inverse of which happens to be singular.

What would be the best way to solve this system for x in OpenMDAO?

Justin Gray · Answer

You have the right instincts here, that this is am implicit problem. You only need to rearrange things to form a residual --- i.e. you want the left hand side to be equal to 0. Once in that form, it becomes more clear how to handle this in OpenMDAO.

Option 1) restructure the components a little and use a BalanceComp

0 = F1(x)/Cx - (F2(x)/Cx-1)

You may not want to express it in that form, but lets assume for now that you were willing to modify you two components that output F1 and F2 into the above form.

Then you could add the BalanceComp to create the necessary residual equation. Like you pointed out, the residual owning component (the balance comp) by itself is singular. You need a group with all three components in it, and you need to use the Newton solver at that group level to create a non singular Jacobian for the solver to work with.

Option 2) Keep F1 and F2 the way they are and define a custom balance comp

What if you didn't want to rearrange those components to us the balance comp? Then you will need to write you own simple implicit component instead. You still need to move the C*x term over to the right hand side in order to form a proper residual though.

import openmdao.api as om

class CustomBalance(om.ImplicitComponent):
    def setup(self):
        self.add_input('F1')
        self.add_input('F2')
        self.add_input('C')
        self.add_output('x')

        self.declare_partials('*', '*', method='cs')

    def apply_nonlinear(self, inputs, outputs, residuals):

        residuals['x'] = inputs['F1'] + inputs['F2'] - inputs['C']*outputs['x']


prob = om.Problem()

ec1 = om.ExecComp('F1=2*x**2')
ec2 = om.ExecComp('F2=-x/2')

prob.model.add_subsystem(name='func1', subsys=ec1, promotes=['x', 'F1'])
prob.model.add_subsystem(name='func2', subsys=ec2, promotes=['x', 'F2'])
prob.model.add_subsystem(name='balance', subsys=CustomBalance(), promotes=['x', 'F1', 'F2', 'C'])

prob.model.linear_solver = om.DirectSolver(assemble_jac=True)
prob.model.nonlinear_solver = om.NewtonSolver(solve_subsystems=False, maxiter=100, iprint=2)

prob.setup()

prob.set_val('C', 2.55)

# initial guess
prob['x'] = 1.0


prob.run_model()

print(prob.get_val('balance.x'))

That gives you a model structure that looks like this:

And a convergence history like this:

NL: Newton 0 ; 1.88480768 1
NL: Newton 1 ; 2.4432133 1.29626663
NL: Newton 2 ; 0.413841411 0.219566916
NL: Newton 3 ; 0.0271563582 0.014408026
NL: Newton 4 ; 0.000154934262 8.22016293e-05
NL: Newton 5 ; 5.16021093e-09 2.73779175e-09
NL: Newton 6 ; 4.4408921e-16 2.35615131e-16
NL: Newton Converged
[1.525]

OpenMDAO implicit relationship between components

1 Answers