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 |
|
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 |
|
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 |
|
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¶
-
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. ↩