Skip to content

Multibranch Trajectory

This example demonstrates the use of a Trajectory to encapsulate a series of branching phases.

Overview

For this example, we build a system that contains two components: the first component represents a battery pack that contains multiple cells in parallel, and the second component represents a bank of DC electric motors (also in parallel) driving a gearbox to achieve a desired power output. The battery cells have a state of charge that decays as current is drawn from the battery. The open circuit voltage of the battery is a function of the state of charge. At any point in time, the coupling between the battery and the motor component is solved with a Newton solver in the containing group for a line current that satisfies the equations.

Both the battery and the motor models allow the number of cells and the number of motors to be modified by setting the n_parallel option in their respective options dictionaries. For this model, we start with 3 cells and 3 motors. We will simulate failure of a cell or battery by setting n_parallel to 2.

Branching phases are a set of linked phases in a trajectory where the input ends of multiple phases are connected to the output of a single phase. This way you can simulate alternative trajectory paths in the same model. For this example, we will start with a single phase (phase0) that simulates the model for one hour. Three follow-on phases will be linked to the output of the first phase: phase1 will run as normal, phase1_bfail will fail one of the battery cells, and phase1_mfail will fail a motor. All three of these phases start where phase0 leaves off, so they share the same initial time and state of charge.

Battery and Motor models

The models are loosely based on the work done in Chin1.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
"""
Simple dynamic model of a LI battery.
"""
import numpy as np
from scipy.interpolate import Akima1DInterpolator

import openmdao.api as om
# Data for open circuit voltage model.
train_SOC = np.array([0., 0.1, 0.25, 0.5, 0.75, 0.9, 1.0])
train_V_oc = np.array([3.5, 3.55, 3.65, 3.75, 3.9, 4.1, 4.2])


class Battery(om.ExplicitComponent):
    """
    Model of a Lithium Ion battery.
    """

    def initialize(self):
        self.options.declare('num_nodes', default=1)
        self.options.declare('n_series', default=1, desc='number of cells in series')
        self.options.declare('n_parallel', default=3, desc='number of cells in parallel')
        self.options.declare('Q_max', default=1.05,
                             desc='Max Energy Capacity of a battery cell in A*h')
        self.options.declare('R_0', default=.025,
                             desc='Internal resistance of the battery (ohms)')

    def setup(self):
        num_nodes = self.options['num_nodes']

        # Inputs
        self.add_input('I_Li', val=np.ones(num_nodes), units='A',
                       desc='Current demanded per cell')

        # State Variables
        self.add_input('SOC', val=np.ones(num_nodes), units=None, desc='State of charge')

        # Outputs
        self.add_output('V_L',
                        val=np.ones(num_nodes),
                        units='V',
                        desc='Terminal voltage of the battery')

        self.add_output('dXdt:SOC',
                        val=np.ones(num_nodes),
                        units='1/s',
                        desc='Time derivative of state of charge')

        self.add_output('V_oc', val=np.ones(num_nodes), units='V',
                        desc='Open Circuit Voltage')
        self.add_output('I_pack', val=0.1*np.ones(num_nodes), units='A',
                        desc='Total Pack Current')
        self.add_output('V_pack', val=9.0*np.ones(num_nodes), units='V',
                        desc='Total Pack Voltage')
        self.add_output('P_pack', val=1.0*np.ones(num_nodes), units='W',
                        desc='Total Pack Power')

        # Derivatives
        row_col = np.arange(num_nodes)

        self.declare_partials(of='V_oc', wrt=['SOC'], rows=row_col, cols=row_col)
        self.declare_partials(of='V_L', wrt=['SOC'], rows=row_col, cols=row_col)
        self.declare_partials(of='V_L', wrt=['I_Li'], rows=row_col, cols=row_col)
        self.declare_partials(of='dXdt:SOC', wrt=['I_Li'], rows=row_col, cols=row_col)
        self.declare_partials(of='I_pack', wrt=['I_Li'], rows=row_col, cols=row_col)
        self.declare_partials(of='V_pack', wrt=['SOC', 'I_Li'], rows=row_col, cols=row_col)
        self.declare_partials(of='P_pack', wrt=['SOC', 'I_Li'], rows=row_col, cols=row_col)

        self.voltage_model = Akima1DInterpolator(train_SOC, train_V_oc)
        self.voltage_model_derivative = self.voltage_model.derivative()

    def compute(self, inputs, outputs):
        opt = self.options
        I_Li = inputs['I_Li']
        SOC = inputs['SOC']

        V_oc = self.voltage_model(SOC, extrapolate=True)

        outputs['V_oc'] = V_oc
        outputs['V_L'] = V_oc - (I_Li * opt['R_0'])
        outputs['dXdt:SOC'] = -I_Li / (3600.0 * opt['Q_max'])

        outputs['I_pack'] = I_Li * opt['n_parallel']
        outputs['V_pack'] = outputs['V_L'] * opt['n_series']
        outputs['P_pack'] = outputs['I_pack'] * outputs['V_pack']

    def compute_partials(self, inputs, partials):
        opt = self.options
        I_Li = inputs['I_Li']
        SOC = inputs['SOC']

        dV_dSOC = self.voltage_model_derivative(SOC, extrapolate=True)
        partials['V_oc', 'SOC'] = dV_dSOC
        partials['V_L', 'SOC'] = dV_dSOC

        partials['V_L', 'I_Li'] = -opt['R_0']

        partials['dXdt:SOC', 'I_Li'] = -1./(3600.0*opt['Q_max'])

        n_parallel = opt['n_parallel']
        n_series = opt['n_series']
        V_oc = self.voltage_model(SOC, extrapolate=True)
        V_L = V_oc - (I_Li * opt['R_0'])

        partials['I_pack', 'I_Li'] = n_parallel
        partials['V_pack', 'I_Li'] = -opt['R_0']
        partials['V_pack', 'SOC'] = n_series * dV_dSOC
        partials['P_pack', 'I_Li'] = n_parallel * n_series * (V_L - I_Li * opt['R_0'])
        partials['P_pack', 'SOC'] = n_parallel * I_Li * n_series * dV_dSOC


if __name__ == '__main__':

    import openmdao.api as om
    num_nodes = 1

    prob = om.Problem(model=Battery(num_nodes=num_nodes))
    model = prob.model

    prob.setup()
    prob.set_solver_print(level=2)

    prob.run_model()

    derivs = prob.check_partials(compact_print=True)

    print('done')

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
"""
Simple model for a set of motors in parallel where efficiency is a function of current.
"""
import numpy as np

import openmdao.api as om


class Motors(om.ExplicitComponent):
    """
    Model for motors in parallel.
    """

    def initialize(self):
        self.options.declare('num_nodes', default=1)
        self.options.declare('n_parallel', default=3, desc='number of motors in parallel')

    def setup(self):
        num_nodes = self.options['num_nodes']

        # Inputs
        self.add_input('power_out_gearbox', val=3.6*np.ones(num_nodes), units='W',
                       desc='Power at gearbox output')
        self.add_input('current_in_motor', val=np.ones(num_nodes), units='A',
                       desc='Total current demanded')

        # Outputs
        self.add_output('power_in_motor', val=np.ones(num_nodes), units='W',
                        desc='Power required at motor input')

        # Derivatives
        row_col = np.arange(num_nodes)

        self.declare_partials(of='power_in_motor', wrt=['*'], rows=row_col, cols=row_col)

    def compute(self, inputs, outputs):
        current = inputs['current_in_motor']
        power_out = inputs['power_out_gearbox']
        n_parallel = self.options['n_parallel']

        # Simple linear curve fit for efficiency.
        eff = 0.9 - 0.3 * current / n_parallel

        outputs['power_in_motor'] = power_out / eff

    def compute_partials(self, inputs, partials):
        current = inputs['current_in_motor']
        power_out = inputs['power_out_gearbox']
        n_parallel = self.options['n_parallel']

        eff = 0.9 - 0.3 * current / n_parallel

        partials['power_in_motor', 'power_out_gearbox'] = 1.0 / eff
        partials['power_in_motor', 'current_in_motor'] = 0.3 * power_out / (n_parallel * eff**2)


if __name__ == '__main__':

    import openmdao.api as om
    num_nodes = 1

    prob = om.Problem(model=Motors(num_nodes=num_nodes))
    model = prob.model

    prob.setup()

    prob.run_model()

    derivs = prob.check_partials(compact_print=True)

    print('done')

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
"""
ODE for example that shows how to use multiple phases in Dymos to model failure of a battery cell
in a simple electrical system.
"""
import numpy as np
import openmdao.api as om

from dymos.examples.battery_multibranch.batteries import Battery
from dymos.examples.battery_multibranch.motors import Motors


class BatteryODE(om.Group):

    def initialize(self):
        self.options.declare('num_nodes', default=1)
        self.options.declare('num_battery', default=3)
        self.options.declare('num_motor', default=3)

    def setup(self):
        num_nodes = self.options['num_nodes']
        num_battery = self.options['num_battery']
        num_motor = self.options['num_motor']

        self.add_subsystem(name='pwr_balance',
                           subsys=om.BalanceComp(name='I_Li', val=1.0*np.ones(num_nodes),
                                                 rhs_name='pwr_out_batt',
                                                 lhs_name='P_pack',
                                                 units='A', eq_units='W', lower=0.0, upper=50.))

        self.add_subsystem('battery', Battery(num_nodes=num_nodes, n_parallel=num_battery),
                           promotes_inputs=['SOC'],
                           promotes_outputs=['dXdt:SOC'])

        self.add_subsystem('motors', Motors(num_nodes=num_nodes, n_parallel=num_motor))

        self.connect('battery.P_pack', 'pwr_balance.P_pack')
        self.connect('motors.power_in_motor', 'pwr_balance.pwr_out_batt')
        self.connect('pwr_balance.I_Li', 'battery.I_Li')
        self.connect('battery.I_pack', 'motors.current_in_motor')

        self.nonlinear_solver = om.NewtonSolver(solve_subsystems=False, maxiter=20)
        self.linear_solver = om.DirectSolver()

Building and running the problem

import matplotlib.pyplot as plt

import openmdao.api as om
from openmdao.utils.assert_utils import assert_near_equal

import dymos as dm
from dymos.examples.battery_multibranch.battery_multibranch_ode import BatteryODE
from dymos.utils.lgl import lgl

prob = om.Problem()

opt = prob.driver = om.ScipyOptimizeDriver()
opt.declare_coloring()
opt.options['optimizer'] = 'SLSQP'

num_seg = 5
seg_ends, _ = lgl(num_seg + 1)

traj = prob.model.add_subsystem('traj', dm.Trajectory())

# First phase: normal operation.
transcription = dm.Radau(num_segments=num_seg, order=5, segment_ends=seg_ends, compressed=False)
phase0 = dm.Phase(ode_class=BatteryODE, transcription=transcription)
traj_p0 = traj.add_phase('phase0', phase0)

traj_p0.set_time_options(fix_initial=True, fix_duration=True)
traj_p0.add_state('state_of_charge', fix_initial=True, fix_final=False,
                  targets=['SOC'], rate_source='dXdt:SOC')

# Second phase: normal operation.

phase1 = dm.Phase(ode_class=BatteryODE, transcription=transcription)
traj_p1 = traj.add_phase('phase1', phase1)

traj_p1.set_time_options(fix_initial=False, fix_duration=True)
traj_p1.add_state('state_of_charge', fix_initial=False, fix_final=False,
                  targets=['SOC'], rate_source='dXdt:SOC')
traj_p1.add_objective('time', loc='final')

# Second phase, but with battery failure.

phase1_bfail = dm.Phase(ode_class=BatteryODE, ode_init_kwargs={'num_battery': 2},
                        transcription=transcription)
traj_p1_bfail = traj.add_phase('phase1_bfail', phase1_bfail)

traj_p1_bfail.set_time_options(fix_initial=False, fix_duration=True)
traj_p1_bfail.add_state('state_of_charge', fix_initial=False, fix_final=False,
                        targets=['SOC'], rate_source='dXdt:SOC')

# Second phase, but with motor failure.

phase1_mfail = dm.Phase(ode_class=BatteryODE, ode_init_kwargs={'num_motor': 2},
                        transcription=transcription)
traj_p1_mfail = traj.add_phase('phase1_mfail', phase1_mfail)

traj_p1_mfail.set_time_options(fix_initial=False, fix_duration=True)
traj_p1_mfail.add_state('state_of_charge', fix_initial=False, fix_final=False,
                        targets=['SOC'], rate_source='dXdt:SOC')

traj.link_phases(phases=['phase0', 'phase1'], vars=['state_of_charge', 'time'])
traj.link_phases(phases=['phase0', 'phase1_bfail'], vars=['state_of_charge', 'time'])
traj.link_phases(phases=['phase0', 'phase1_mfail'], vars=['state_of_charge', 'time'])

prob.model.options['assembled_jac_type'] = 'csc'
prob.model.linear_solver = om.DirectSolver(assemble_jac=True)

prob.setup()

prob['traj.phase0.t_initial'] = 0
prob['traj.phase0.t_duration'] = 1.0*3600

prob['traj.phase1.t_initial'] = 1.0*3600
prob['traj.phase1.t_duration'] = 1.0*3600

prob['traj.phase1_bfail.t_initial'] = 1.0*3600
prob['traj.phase1_bfail.t_duration'] = 1.0*3600

prob['traj.phase1_mfail.t_initial'] = 1.0*3600
prob['traj.phase1_mfail.t_duration'] = 1.0*3600

prob.set_solver_print(level=0)
dm.run_problem(prob)

soc0 = prob['traj.phase0.states:state_of_charge']
soc1 = prob['traj.phase1.states:state_of_charge']
soc1b = prob['traj.phase1_bfail.states:state_of_charge']
soc1m = prob['traj.phase1_mfail.states:state_of_charge']

# Final value for State of Chrage in each segment should be a good test.
print('State of Charge after 1 hour')

print('State of Charge after 2 hours')

print('State of Charge after 2 hours, battery fails at 1 hour')

print('State of Charge after 2 hours, motor fails at 1 hour')


# Plot Results
t0 = prob['traj.phases.phase0.time.time']/3600
t1 = prob['traj.phases.phase1.time.time']/3600
t1b = prob['traj.phases.phase1_bfail.time.time']/3600
t1m = prob['traj.phases.phase1_mfail.time.time']/3600

plt.subplot(2, 1, 1)
plt.plot(t0, soc0, 'b')
plt.plot(t1, soc1, 'b')
plt.plot(t1b, soc1b, 'r')
plt.plot(t1m, soc1m, 'c')
plt.xlabel('Time (hour)')
plt.ylabel('State of Charge (percent)')

I_Li0 = prob['traj.phases.phase0.rhs_all.pwr_balance.I_Li']
I_Li1 = prob['traj.phases.phase1.rhs_all.pwr_balance.I_Li']
I_Li1b = prob['traj.phases.phase1_bfail.rhs_all.pwr_balance.I_Li']
I_Li1m = prob['traj.phases.phase1_mfail.rhs_all.pwr_balance.I_Li']

plt.subplot(2, 1, 2)
plt.plot(t0, I_Li0, 'b')
plt.plot(t1, I_Li1, 'b')
plt.plot(t1b, I_Li1b, 'r')
plt.plot(t1m, I_Li1m, 'c')
plt.xlabel('Time (hour)')
plt.ylabel('Line Current (A)')

plt.legend(['Phase 1', 'Phase 2', 'Phase 2 Battery Fail', 'Phase 2 Motor Fail'], loc=2)

plt.show()

--- Linkage Report [traj] ---
    --- phase0 - phase1 ---
        state_of_charge [final]  ==  state_of_charge [initial]
        time            [final]  ==  time            [initial]
----------------------------
    --- phase0 - phase1_bfail ---
        state_of_charge [final]  ==  state_of_charge [initial]
        time            [final]  ==  time            [initial]
----------------------------
    --- phase0 - phase1_mfail ---
        state_of_charge [final]  ==  state_of_charge [initial]
        time            [final]  ==  time            [initial]
----------------------------
Full total jacobian was computed 3 times, taking 0.127775 seconds.
Total jacobian shape: (107, 122) 


Jacobian shape: (107, 122)  ( 4.63% nonzero)
FWD solves: 6   REV solves: 0
Total colors vs. total size: 6 vs 122  (95.1% improvement)

Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.127775 sec.
Time to compute coloring: 0.002988 sec.
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 7200.0
            Iterations: 3
            Function evaluations: 4
            Gradient evaluations: 3
Optimization Complete
-----------------------------------
State of Charge after 1 hour
State of Charge after 2 hours
State of Charge after 2 hours, battery fails at 1 hour
State of Charge after 2 hours, motor fails at 1 hour

References

  1. Jeff Chin, Sydney L Schnulo, Thomas Miller, Kevin Prokopius, and Justin S Gray. Battery performance modeling on sceptor x-57 subject to thermal and transient considerations. In AIAA Scitech 2019 Forum, 0784. 2019. ↩