Converging an Implicit Model: Non-linear circuit analysis

Consider a simple electrical circuit made up from two resistors, a diode, and a constant current source. Our goal is to solve for the steady-state voltages at node 1 and node 2.

diagram of a simple circuit with two resistors and one diode

In order to find the voltages, we’ll employ Kirchoff’s voltage law, and solve for the voltages needed at each node to drive the net-current to 0.

import numpy as np

from openmdao.api import ExplicitComponent, ImplicitComponent,  NewtonSolver, DirectSolver, ArmijoGoldsteinLS
from openmdao.api import IndepVarComp, Problem, Group

class Resistor(ExplicitComponent):
    def initialize(self):
        self.metadata.declare('R', default=1., desc='Resistance in Ohms')

    def setup(self):
        self.add_input('V_in', units='V')
        self.add_input('V_out', units='V')
        self.add_output('I', units='A')

        # partial derivs are constant, so we can assign their values in setup
        R = self.metadata['R']
        self.declare_partials('I', 'V_in', val=1 / R)
        self.declare_partials('I', 'V_out', val=-1 / R)

    def compute(self, inputs, outputs):
        deltaV = inputs['V_in'] - inputs['V_out']
        outputs['I'] = deltaV / self.metadata['R']


class Diode(ExplicitComponent):
    def initialize(self):
        self.metadata.declare('Is', default=1e-15, desc='Saturation current in Amps')
        self.metadata.declare('Vt', default=.025875, desc='Thermal voltage in Volts')

    def setup(self):
        self.add_input('V_in', units='V')
        self.add_input('V_out', units='V')
        self.add_output('I', units='A')

        # non-linear component, so we'll declare the partials here but compute them in compute_partials
        self.declare_partials('I', 'V_in')
        self.declare_partials('I', 'V_out')

    def compute(self, inputs, outputs):
        deltaV = inputs['V_in'] - inputs['V_out']
        Is = self.metadata['Is']
        Vt = self.metadata['Vt']
        outputs['I'] = Is * np.exp(deltaV / Vt - 1)

    def compute_partials(self, inputs, J):
        deltaV = inputs['V_in'] - inputs['V_out']
        Is = self.metadata['Is']
        Vt = self.metadata['Vt']
        I = Is*np.exp(deltaV/Vt-1)

        J['I', 'V_in'] = I/Vt
        J['I', 'V_out'] = -I/Vt

class Node(ImplicitComponent):
    def initialize(self):
        self.metadata.declare('n_in', default=1, types=int, desc='number of connections with + assumed in')
        self.metadata.declare('n_out', default=1, types=int, desc='number of current connections + assumed out')

    def setup(self):
        self.add_output('V', val=5., units='V')

        for i in range(self.metadata['n_in']):
            i_name = 'I_in:{}'.format(i)
            self.add_input(i_name, units='A')
            self.declare_partials('V', i_name, val=1)

        for i in range(self.metadata['n_out']):
            i_name = 'I_out:{}'.format(i)
            self.add_input(i_name, units='A')
            self.declare_partials('V', i_name, val=-1)

            # note: we don't declare any partials wrt `V` here,
            #      because the residual doesn't directly depend on it

    def apply_nonlinear(self, inputs, outputs, residuals):
        residuals['V'] = 0.
        for i_conn in range(self.metadata['n_in']):
            residuals['V'] += inputs['I_in:{}'.format(i_conn)]
        for i_conn in range(self.metadata['n_out']):
            residuals['V'] -= inputs['I_out:{}'.format(i_conn)]

class Circuit(Group):

    def setup(self):
        self.add_subsystem('n1', Node(n_in=1, n_out=2), promotes_inputs=[('I_in:0', 'I_in')])
        self.add_subsystem('n2', Node())  # leaving defaults

        self.add_subsystem('R1', Resistor(R=100.), promotes_inputs=[('V_out', 'Vg')])
        self.add_subsystem('R2', Resistor(R=10000.))
        self.add_subsystem('D1', Diode(), promotes_inputs=[('V_out', 'Vg')])

        self.connect('n1.V', ['R1.V_in', 'R2.V_in'])
        self.connect('R1.I', 'n1.I_out:0')
        self.connect('R2.I', 'n1.I_out:1')

        self.connect('n2.V', ['R2.V_out', 'D1.V_in'])
        self.connect('R2.I', 'n2.I_in:0')
        self.connect('D1.I', 'n2.I_out:0')

        self.nonlinear_solver = NewtonSolver()
        self.linear_solver = DirectSolver()

        self.nonlinear_solver.options['iprint'] = 2
        self.nonlinear_solver.options['maxiter'] = 10
        self.nonlinear_solver.options['solve_subsystems'] = True
        self.nonlinear_solver.linesearch = ArmijoGoldsteinLS()
        self.nonlinear_solver.linesearch.options['maxiter'] = 10
        self.nonlinear_solver.linesearch.options['iprint'] = 2


p = Problem()
model = p.model

model.add_subsystem('ground', IndepVarComp('V', 0., units='V'))
model.add_subsystem('source', IndepVarComp('I', 0.1, units='A'))
model.add_subsystem('circuit', Circuit())

model.connect('source.I', 'circuit.I_in')
model.connect('ground.V', 'circuit.Vg')

p.setup()

# set some initial guesses
p['circuit.n1.V'] = 10.
p['circuit.n2.V'] = 1e-3

p.run_model()
print(p['circuit.n1.V'])
print(p['circuit.n2.V'])
print(p['circuit.R1.I'])
print(p['circuit.R2.I'])
print(p['circuit.D1.I'])
# sanity check: should sum to .1 Amps
print(p['circuit.R1.I'] + p['circuit.D1.I'])
=======
circuit
=======
NL: Newton 0 ; 0.00141407214 1
|  LS: AG 0 ; 2.56438674e+152 1.8134766e+155
|  LS: AG 1 ; 3.13138657e+68 2.21444612e+71
|  LS: AG 3 ; 4.42844936e+68 3.13169974e+71
|  LS: AG 5 ; 4.8935977e+26 3.46064218e+29
|  LS: AG 7 ; 514418.469 363785167
|  LS: AG 9 ; 0.00136161968 0.962906799
NL: Newton 1 ; 0.00136988234 0.968749966
|  LS: AG 0 ; 2.53928001e+152 1.85364826e+155
|  LS: AG 1 ; 1.30583764e+71 9.53248032e+73
|  LS: AG 3 ; 1.84673331e+71 1.34809629e+74
|  LS: AG 5 ; 4.18787274e+30 3.05710398e+33
|  LS: AG 7 ; 1.99428914e+10 1.45581053e+13
|  LS: AG 9 ; 1.37554356 1004.13263
NL: Newton 2 ; 0.00132137032 0.934443361
|  LS: AG 0 ; 3.02186757e+32 2.28691951e+35
|  LS: AG 1 ; 4.94250674e+13 3.74044024e+16
|  LS: AG 3 ; 6.98976006e+13 5.28978132e+16
|  LS: AG 5 ; 28267.9172 21392880.3
|  LS: AG 7 ; 0.566492469 428.715901
|  LS: AG 9 ; 0.00207462986 1.57005937
NL: Newton 3 ; 0.00121787427 0.861253278
|  LS: AG 0 ; 0.0794819095 65.2628204
|  LS: AG 1 ; 0.00267232126 2.19425054
|  LS: AG 3 ; 0.00312376017 2.5649283
|  LS: AG 5 ; 0.000928550905 0.762435771
NL: Newton 4 ; 0.00103096282 0.72907371
|  LS: AG 0 ; 0.00196603165 1.90698598
|  LS: AG 1 ; 0.000592516181 0.574721191
NL: Newton 5 ; 0.000399319798 0.282389976
|  LS: AG 0 ; 5.76512173e-07 0.00144373551
NL: Newton 6 ; 5.76512173e-07 0.000407696436
|  LS: AG 0 ; 1.80072036e-10 0.000312347326
NL: Newton 7 ; 1.80072036e-10 1.27342892e-07
|  LS: AG 0 ; 1.89434663e-17 1.05199379e-07
NL: Newton 8 ; 1.89434663e-17 1.33963931e-14
NL: Newton Converged
[ 9.90830282]
[ 0.73858486]
[ 0.09908303]
[ 0.00091697]
[ 0.00091697]
[ 0.1]