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

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):

# partial derivs are constant, so we can assign their values in setup
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):

# 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']
outputs['I'] = Is * np.exp(deltaV / Vt - 1)

def compute_partials(self, inputs, J):
deltaV = inputs['V_in'] - inputs['V_out']
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):

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

for i in range(self.metadata['n_out']):
i_name = 'I_out:{}'.format(i)
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('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.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]