Edit 1
Hi Justin,
Sorry for the confusion. I have edited & updated the question now after making some progress. Also adding an N2 diagram to make it more readable - N2 Diagram
I am trying to build a dymos model which takes in electrical power-time series (power_out_gearbox) and sizes the fuel cell system and battery system to deliver the required power at different time steps using a phase.parameter called split_ratio.
With respect to fuel cell system, it is divided into n_fcsysmodules. n_fcsysmodules are further split into nminus1_active_fcsysmodules and nth_active_fcsysmodule. Additionally, the fuel, oxidant, and water rates (dXdt:) are computed alongwith mass_fcsys and cost_fcsys
With respect to battery system, it outputs SoC rate (dXdt:), mass_batsys and cost_fcsys.
The mass_batsys, mass_fcsys, and the mass_fuel, mass_oxidant, and mass_water are added together to obtain takeoff mass. The total takeoff mass is then minimized at loc=initial.
I have made a sample code (below) but I am running into some errors. Could you please help with the following?
1 - I am also not sure about defining the input elecpower-time series. I saw a post on stackoverflow answered by rob and tried to construct as shown in below code (Code 1). But I am not sure if this is the right approach. Could you please help?
2 - Dymos documentation says I can set opt=False for a state variable and tried to do that with state_of_charge, but i ran into constraint does not depend on design variable error. Could you please share your thoughts on having a state variable with no targets and is not a design variable but only has rate source?
3 - With respect to phase.parameters, when i set static_target=True (Error 1) i get array indices error and when i set it to false, i get a Nan error (Error 2). Could you please help on what is appropriate in my case?
If you could also share your thoughts or possible improvements on other aspects of the code, that would be amazing. Thank you in advance.
best regards
# Code 1
import openmdao.api as om
import dymos as dm
from dymos.utils.lgl import lgl
# from ExtMC import ExtMCGroup
# from FCBatSysODE3ExtMCGroup import FCBatSysODEGroup
from FCBatSysODE4 import FCBatSysODEGroup
import matplotlib.pyplot as plt
import numpy as np
# -------------------------------Model Creation-------------------------------
model = om.Group()
prob = om.Problem(model)
# -------------------------------Model Creation-------------------------------
# -------------------------------Driver Settings-------------------------------
optimizer = prob.driver = om.ScipyOptimizeDriver()
# optimizer = prob.driver = om.SimpleGADriver()
optimizer.declare_coloring()
# prob.driver.declare_coloring()
optimizer.options['optimizer'] = 'SLSQP'
# -------------------------------Driver Settings-------------------------------
# -------------------------------Input Power Series-------------------------------
# power_profile = [0, 50000, 100000, 100000, 100000, 50000, 0]
# power_profile = np.array([0, 50000, 100000, 100000, 100000, 50000, 0])
# time_duration = len(power_profile)
# num_control_points = len(power_profile)
# # Note: LGL nodes are in the range [-1, 1] (this is tau-space)
# nodelocations, nodeweights = lgl(num_control_points)
# -------------------------------Input Power Series-------------------------------
# -------------------------------------IVC-------------------------------------
ivc = prob.model.add_subsystem('states_ivc', om.IndepVarComp(), promotes_outputs=['*'])
ivc.add_output('current_maxfcstack', val=642, units=None) # A
# ivc.add_output('pwr_max_out_gearbox', val=max(power_profile), units=None) # J/s
ivc.add_output('pwr_dens_fcstack', val=2, units=None) # kW/kg
# ivc.add_design_var('t_duration', units='s', lower=0.1, upper=10.)
# -------------------------------------IVC-------------------------------------
# -------------------------------External Mass and Cost Group-------------------------------
# prob.model.add_subsystem('mass_and_cost', ExtMCGroup(), promotes_inputs = ['current_maxfcstack', 'pwr_max_out_gearbox', 'split_ratio', 'pwr_dens_fcstack', 'egy_batsys', 'mass_fuel', 'mass_oxidant', 'mass_water'], promotes_outputs=['pwr_el_max_in_motors', 'pwr_el_max_fcsys', 'mass_fcsys', 'cost_fcsys', 'mass_batsys', 'cost_batsys', 'total_mass_fcbatsys','tot_cost_fcbatsys', 'tot_mass_takeoff'])
# -------------------------------External Mass and Cost Group-------------------------------
# ----------------Trajectory Creation and Transcription Settings---------------
traj = prob.model.add_subsystem('traj', dm.Trajectory())
# //[ ]: Check and change num_segments1
num_segments1 = 7
transcription1 = dm.GaussLobatto(num_segments=num_segments1, compressed=False)
# phase = dm.Phase(ode_class=FCBatSysODEGroup, transcription=transcription1) # ode_class=MinTimeClimbODE,
phase = traj.add_phase('phase0', dm.Phase(ode_class=FCBatSysODEGroup, transcription=transcription1))
######################## Radau Settings ########################
# num_seg = 5
# seg_ends, _ = lgl(num_seg + 1)
# transcription2 = dm.Radau(num_segments=num_seg, segment_ends=seg_ends, order=3, compressed=False)
# phase = traj.add_phase('phase0', dm.Phase(ode_class=FCBatSysODEGroup, transcription=transcription2))
# First phase: normal operation. # Multibranch Example
# transcription3 = dm.Radau(num_segments=num_seg, segment_ends=seg_ends, order=5, compressed=False)
# phase0 = dm.Phase(ode_class=BatteryODE, transcription=transcription3)
# traj_p0 = traj.add_phase('phase0', phase0)
# ----------------Trajectory Creation and Transcription Settings---------------
# -------------------------------Polynomial Controls-------------------------------
power_profile = [0, 50000, 100000, 100000, 100000, 50000, 0]
time_duration = len(power_profile) # Number of points
order_polynomial = time_duration - 1
# Polynomial controls are defined by values at the Legendre-Gauss-Lobatto nodes of a single polynomial, defined on [-1, 1] in phase tau space.
# For a polynomial control of a given order, the number of nodes used to define the polynomial is (order + 1). Nodes are like time steps, and if the number of nodes is order + 1, that is equal to number of time steps. Seems reasonable. But have to che
phase.add_polynomial_control('pwr_out_gearbox', opt=False, units=None, val=power_profile, order=order_polynomial, targets='pwr_out_gearbox')
# -------------------------------Polynomial Controls-------------------------------
# -------------------------------States and Boundary Constraints-------------------------------
# //[ ]: Check and adjust limits for states: mass_fuel, mass_oxidant, mass_water, mass_emission.
# //[ ]: Check the meaning of ref, detect_ref, and scaler for states and boundary constraints.
# //[x]: Limits () & //[ ]: Units [-] to be added.
phase.add_state('state_of_charge', opt=False, units=None, fix_initial=True, fix_final=True, targets=None, rate_source='dXdt:SoC')
phase.add_boundary_constraint('state_of_charge', loc='initial', equals=1, scaler=1.0E-2)
phase.add_boundary_constraint('state_of_charge', loc='final', equals=0.3, scaler=1.0E-2)
# //[x]: Limits () & //[ ]: Units [kg] to be added.
phase.add_state('mass_fuel', units=None, rate_source='dXdt:mass_fuel', targets='mass_fuel', fix_initial=False, fix_final=True, upper=3.0E7, lower=0.0, ref=1e2, defect_ref=1e2)
phase.add_boundary_constraint('mass_fuel', loc='final', equals=0, scaler=1.0E-2)
# //[x]: Limits () & //[ ]: Units [kg] to be added.
phase.add_state('mass_oxidant', units=None, rate_source='dXdt:mass_oxidant', targets='mass_oxidant', fix_initial=False, fix_final=True, upper=3.0E6, lower=0.0, ref=1e2, defect_ref=1e2)
phase.add_boundary_constraint('mass_oxidant', loc='final', equals=0, scaler=1.0E-2)
# //[x]: Limits () & //[ ]: Units [kg] to be added.
phase.add_state('mass_water', units=None, rate_source='dXdt:mass_water', targets='mass_water', fix_initial=False, fix_final=False, upper=4.0E6, lower=0.0, ref=1e2, defect_ref=1e2)
# //BUG: Actually mass_water must be some number at initial time step before fc runs. The fix_inital is set to True and the boundary condition is removed.
# phase.add_boundary_constraint('mass_water', loc='initial', equals=0, scaler=1.0E-2)
# //[x]: Limits () & //[ ]: Units [kg] to be added.
phase.add_state('mass_emission', units=None, rate_source='dXdt:mass_emission', fix_initial=True, fix_final=False, targets=None, upper=4.0E7, lower=0.0, ref=1e2, defect_ref=1e2)
phase.add_boundary_constraint('mass_emission', loc='initial', equals=0, scaler=1.0E-2)
# -------------------------------States and Boundary Constraints-------------------------------
# -------------------------------Controls-------------------------------
# //[ ]: Understand continuity, rate_continuity, rate2_continuity, rate_continuity_scaler for controls.
# //[ ]: See how to use rate_targets and rate2_targets to set the ramp-up and ramp-down rates working with a constraint in the ODE.
# //[ ]: See if the control values need fix_initial and fix_final. For example at the end of mission C_rate_dischg must be 0 or something like that.
# //[ ]: Limits () & //[ ]: Units [1/h] to be added.
phase.add_control('C_rate_chg', units=None, lower=0.0/1.0, upper=1.0/1.0, scaler=1.0E-4, opt=True, continuity=True, rate_continuity=False, rate2_continuity=False) # rate_continuity_scaler=100.0
# //[ ]: Limits () & //[ ]: Units [1/h] to be added.
phase.add_control('C_rate_dischg', units=None, lower=0.0/1.0, upper=1.0/1.0, scaler=1.0E-4, opt=True, continuity=True, rate_continuity=False, rate2_continuity=False) # rate_continuity_scaler=100.0,
# //[ ]: Limits () & //[ ]: Units [-] to be added.
phase.add_control('nminus1_active_fcsysmodules', units=None, lower=0.0, upper=1000, scaler=1.0, opt=True, continuity=False, rate_continuity=False, rate2_continuity=False) # rate_continuity_scaler=100.0,
# //[ ]: Limits () & //[ ]: Units [A] to be added.
phase.add_control('current_nminus1fcstack', units=None, lower=0.0, upper=642.0, scaler=1.0E-4, opt=True, continuity=True, rate_continuity=False, rate2_continuity=False) # rate_continuity_scaler=100.0
# //[ ]: Limits () & //[ ]: Units [A] to be added.
phase.add_control('current_nthfcstack', units=None, lower=0.0, upper=642.0, scaler=1.0E-4, opt=True, continuity=True, rate_continuity=False, rate2_continuity=False) # rate_continuity_scaler=100.0,
# -------------------------------Controls-------------------------------
# -------------------------------Parameters-------------------------------
# //[ ]: See if egy_batsys_default can be an ivc comp output.
# //[ ]: Limits (None) & //[ ]: Units [J/s] to be added.
phase.add_parameter('pwr_max_out_gearbox', units=None, val=max(power_profile), opt=False, targets=['pwr_max_out_gearbox'], static_target=True)
# //[x]: Limits (0.01 to 0.99) & //[ ]: Units [-] to be added.
phase.add_parameter('split_ratio', units=None, val=0.5, lower=0.01, upper=0.99, opt=True, targets=['split_ratio'], static_target=True)
egy_batsys_default = 1000 #kWh or 1 MWh = 3.6 GJ or 3.6E10J
def kWh_to_J(kWh):
"""
Converts kilowatt-hours (kWh) to joules (J).
Parameters:
kWh (float): The energy in kilowatt-hours.
Returns:
float: The energy in joules.
"""
joules = kWh * 3600 * 1000 # 1000 * [J/s] * 3600 [s]
return joules
# //[ ]: Limits (0.01 to 1E10) & //[ ]: Units [J] to be added.
# //[x] Add a function that takes in kWh and converts to J.
phase.add_parameter('egy_batsys', units=None, val=kWh_to_J(egy_batsys_default), lower=0.01, upper=1E12, opt=True, targets=['egy_batsys'], static_target=True)
# -------------------------------Parameters-------------------------------
prob.model.options['assembled_jac_type'] = 'csc'
prob.model.linear_solver = om.DirectSolver(assemble_jac=True)
# -------------------------------Objective Function-------------------------------
phase.add_objective('tot_mass_takeoff', loc='initial', ref=1.0)
# -------------------------------Objective Function-------------------------------
prob.setup(check=True)
phase.set_time_options(fix_initial=True, fix_duration=True, duration_val=time_duration)
prob['traj.phase0.t_initial'] = 0
prob['traj.phase0.t_duration'] = time_duration
prob.set_val('phase.polynomial_controls:pwr_out_gearbox', power_profile, units=None)
prob['traj.phase0.states:state_of_charge'] = phase.interp('state_of_charge', [1, 0.3])
prob['traj.phase0.states:mass_fuel'] = phase.interp('mass_fuel', [15000, 0])
prob['traj.phase0.states:mass_oxidant'] = phase.interp('mass_oxidant', [5000, 0])
prob['traj.phase0.states:mass_water'] = phase.interp('mass_water', [0, 15000])
prob['traj.phase0.states:mass_emission'] = phase.interp('mass_emission', [0, 25000])
prob['traj.phase0.controls:C_rate_chg'] = phase.interp('C_rate_chg', [0.0, 1])
prob['traj.phase0.controls:C_rate_dischg'] = phase.interp('C_rate_dischg', [1.0, 0.0])
prob['traj.phase0.controls:nminus1_active_fcsysmodules'] = phase.interp('nminus1_active_fcsysmodules', [0.0, 80])
prob['traj.phase0.controls:current_nminus1fcstack'] = phase.interp('current_nminus1fcstack', [0.0, 642])
prob['traj.phase0.controls:current_nthfcstack'] = phase.interp('current_nthfcstack', [0.0, 642])
prob.set_solver_print(level=1)
dm.run_problem(prob, simulate=False)
print('Optimization finished')
Error 1
{
"name": "RuntimeError",
"message": "
Collected errors for problem 'problem':
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:pwr_max_out_gearbox' to 'traj.phases.phase0.rhs_col.max_elecpower_inmotors.pwr_max_out_gearbox'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:pwr_max_out_gearbox' to 'traj.phases.phase0.rhs_col.power_input.pwr_max_out_gearbox'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:split_ratio' to 'traj.phases.phase0.rhs_col.power_max_fcsys.split_ratio'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:split_ratio' to 'traj.phases.phase0.rhs_col.power_splitter.split_ratio'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:egy_batsys' to 'traj.phases.phase0.rhs_col.batsys_massandcost.egy_batsys'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:egy_batsys' to 'traj.phases.phase0.rhs_col.batterysys.egy_batsys'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:egy_batsys' to 'traj.phases.phase0.rhs_col.power_splitter.egy_batsys'. The target shape is (7,) but indices are shape (14,).",
"stack": "---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
Cell In[1], line 188
184 phase.add_objective('tot_mass_takeoff', loc='initial', ref=1.0)
186 # -------------------------------Objective Function-------------------------------
--> 188 prob.setup(check=True)
190 phase.set_time_options(fix_initial=True, fix_duration=True, duration_val=time_duration)
192 prob['traj.phase0.t_initial'] = 0
File ~/anaconda3/envs/OM-dymos/lib/python3.11/site-packages/openmdao/core/problem.py:1044, in Problem.setup(self, check, logger, mode, force_alloc_complex, distributed_vector_class, local_vector_class, derivatives)
1040 self._logger = logger
1042 self._metadata['setup_status'] = _SetupStatus.POST_SETUP
-> 1044 self._check_collected_errors()
1046 return self
File ~/anaconda3/envs/OM-dymos/lib/python3.11/site-packages/openmdao/utils/hooks.py:131, in _hook_decorator.<locals>.execute_hooks(*args, **kwargs)
129 def execute_hooks(*args, **kwargs):
130 _run_hooks(pre_hooks, inst)
--> 131 ret = f(*args, **kwargs)
132 _run_hooks(post_hooks, inst)
133 return ret
File ~/anaconda3/envs/OM-dymos/lib/python3.11/site-packages/openmdao/core/problem.py:619, in Problem._check_collected_errors(self)
616 else:
617 raise exc_type('\
'.join(final_msg)).with_traceback(tback)
--> 619 raise RuntimeError('\
'.join(final_msg))
RuntimeError:
Collected errors for problem 'problem':
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:pwr_max_out_gearbox' to 'traj.phases.phase0.rhs_col.max_elecpower_inmotors.pwr_max_out_gearbox'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:pwr_max_out_gearbox' to 'traj.phases.phase0.rhs_col.power_input.pwr_max_out_gearbox'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:split_ratio' to 'traj.phases.phase0.rhs_col.power_max_fcsys.split_ratio'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:split_ratio' to 'traj.phases.phase0.rhs_col.power_splitter.split_ratio'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:egy_batsys' to 'traj.phases.phase0.rhs_col.batsys_massandcost.egy_batsys'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:egy_batsys' to 'traj.phases.phase0.rhs_col.batterysys.egy_batsys'. The target shape is (7,) but indices are shape (14,).
'traj.phases.phase0' <class Phase>: The source indices (array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]),) do not specify a valid shape for the connection 'traj.phases.phase0.param_comp.parameter_vals:egy_batsys' to 'traj.phases.phase0.rhs_col.power_splitter.egy_batsys'. The target shape is (7,) but indices are shape (14,)."
}
I get the following NaN errors if I set the parameters' static_target=False.
Error 2
File ~/anaconda3/envs/OM-dymos/lib/python3.11/site-packages/openmdao/core/driver.py:1079, in Driver._compute_totals(self, of, wrt, return_format, use_abs_names, driver_scaling)
1076 self._recording_iter.push(('_compute_totals', 0))
1078 try:
-> 1079 totals = total_jac.compute_totals()
1080 finally:
1081 self._recording_iter.pop()
File ~/anaconda3/envs/OM-dymos/lib/python3.11/site-packages/openmdao/core/total_jac.py:1515, in _TotalJacInfo.compute_totals(self)
1512 if ln_solver._assembled_jac is not None and \
1513 ln_solver._assembled_jac._under_complex_step:
1514 model.linear_solver._assembled_jac._update(model)
-> 1515 ln_solver._linearize()
1516 finally:
1517 model._tot_jac = None
File ~/anaconda3/envs/OM-dymos/lib/python3.11/site-packages/openmdao/solvers/linear/direct.py:288, in DirectSolver._linearize(self)
286 self._lu = scipy.sparse.linalg.splu(matrix)
287 except RuntimeError as err:
--> 288 raise RuntimeError(format_singular_error(system, matrix))
290 elif isinstance(matrix, np.ndarray): # dense
291 # During LU decomposition, detect singularities and warn user.
292 with warnings.catch_warnings():
RuntimeError: NaN entries found in <model> <class Group> for rows associated with states/residuals ['traj.phase0.rhs_disc.fuelcellsys.nminus1_fcsysmodulegroup.eff_nminus1fcsysmodule', 'traj.phase0.rhs_disc.fuelcellsys.nth_fcsysmodulegroup.eff_nthfcsysmodule'].
Edit 2 & 3 - Solved shape and np.any/np.all() errors.