1

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?

Kasia
  • 105
  • 6

1 Answers1

0

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:

N2 diagram of sample model

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]
Justin Gray
  • 5,605
  • 1
  • 11
  • 16
  • Thanks a lot! I have implemented option 2 and it seems to work fine. I tried both NewtonSolver and BroydenSolver, and although both converge, they converge at slightly different values. I can see that for both solvers the default rtol and atol is 1e-10, but at the end I'm getting residuals as large as 1e-3. What is the tolerance defined for, is it residuals, or some other quantity maybe? If it is not on the residuals, is there any way to force the solvers to continue until the residuals are closer to zero? – Kasia Sep 07 '22 at 14:17
  • without seeing your code, its hard to answer exactly why this is happening. I would guess that your initial residual is massive and the rtol is getting triggered first. Try setting rtol to 1e-100 to effectively turn it off and rely on the atol instead. – Justin Gray Sep 08 '22 at 16:00
  • Thank you, it works fine now. One question: in the example you provided, we have 1 implicit and 2 explicit components, grouped into 1 group. We define the NL solver for the group. Does the implicit component not need its own NL solver? The way I understood it, we drive its residual to zero by changing x, keeping F1 and F2 fixed, and then we need another group-level solver to converge between x, F1, and F2. So the group solver and the implicit component's solvers would be called repeatedly until all three states converge. But this is not the case, could you help me understand why, please? – Kasia Sep 14 '22 at 19:14
  • Nonlinear Block Gauss Seidel fundamentally (because of the math and algorithm) requires implicit components to solve themselves. It fundamentally won't work without that happening. Explicit components are also required to "solve" themselves, but that happens by default. This is why I used a newton solver around the CustomBalance. The balance is singular, because the residual does not depend on its own state variables at all. It fundamentally can NOT solve itself, and hence NLBS can't be used for balances. – Justin Gray Sep 15 '22 at 21:01