14268934115

Committed 04 Apr 2025 03:36PM UTC coverage: 92.828% (-0.2%) from 93.071%

Build # 14268934115

Build Type

push

github

Committed by

swryan

Commit Message

Merge branch 'master' into docs_on_tag

Run Details

5 of 13 new or added lines in 3 files covered. (38.46%)

557 existing lines in 26 files now uncovered.

33343 of 35919 relevant lines covered (92.83%)

5.64 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

86.67

/dymos/grid_refinement/error_estimation.py

"""
Generic utilities for use by the grid refinement schemes.
"""
import numpy as np
from scipy.linalg import block_diag

import openmdao.api as om

from ..transcriptions import GaussLobatto, Radau
from ..utils.lagrange import lagrange_matrices
from ..utils.misc import get_rate_units
from ..utils.introspection import get_targets

from .grid_refinement_ode_system import GridRefinementODESystem


def interpolation_lagrange_matrix(old_grid, new_grid):
    """
    Evaluate lagrange matrix to interpolate state and control values from the solved grid onto the new grid.

    Parameters
    ----------
    old_grid : <GridData>
        GridData object representing the grid on which the problem has been solved.
    new_grid : <GridData>
        GridData object representing the new, higher-order grid.

    Returns
    -------
    ndarray
        The lagrange interpolation matrix.
    """
    L_blocks = []
    D_blocks = []

    for iseg in range(old_grid.num_segments):
        i1, i2 = old_grid.subset_segment_indices['all'][iseg, :]
        indices = old_grid.subset_node_indices['all'][i1:i2]
        nodes_given = old_grid.node_stau[indices]

        i1, i2 = new_grid.subset_segment_indices['all'][iseg, :]
        indices = new_grid.subset_node_indices['all'][i1:i2]
        nodes_eval = new_grid.node_stau[indices]

        L_block, D_block = lagrange_matrices(nodes_given, nodes_eval)

        L_blocks.append(L_block)
        D_blocks.append(D_block)

    L = block_diag(*L_blocks)
    D = block_diag(*D_blocks)

    return L, D


def integration_matrix(grid):
    """
    Evaluate the Integration matrix of the given grid.

    Parameters
    ----------
    grid : <GridData>
        GridData object containing the grid where the integration matrix will be evaluated.

    Returns
    -------
    ndarray
        The integration matrix used to propagate initial states over segments.
    """
    I_blocks = []

    for iseg in range(grid.num_segments):
        i1, i2 = grid.subset_segment_indices['all'][iseg, :]
        indices = grid.subset_node_indices['all'][i1:i2]
        nodes_given = grid.node_stau[indices]

        i1, i2 = grid.subset_segment_indices['all'][iseg, :]
        indices = grid.subset_node_indices['all'][i1:i2]
        nodes_eval = grid.node_stau[indices][1:]

        _, D_block = lagrange_matrices(nodes_given, nodes_eval)
        I_block = np.linalg.inv(D_block[:, 1:])
        I_blocks.append(I_block)

    return block_diag(*I_blocks)


def eval_ode_on_grid(phase, transcription):
    """
    Evaluate the ODE from the given phase at all nodes of the given transcription.

    Parameters
    ----------
    phase : Phase
        A Phase object which has been executed and whose ODE is to be run at all nodes
        of the given transcription.
    transcription : Radau or GaussLobatto transcription instance
        The transcription object at which to execute the ODE of the given phase at all nodes.

    Returns
    -------
    dict of (str: ndarray)
        A dictionary of the state values from the phase interpolated to the new transcription.
    dict of (str: ndarray)
        A dictionary of the control values from the phase interpolated to the new transcription.
    dict of (str: ndarray)
        A dictionary of the state rates computed in the phase's ODE at the new transcription points.
    """
    x = {}
    u = {}
    u_rate = {}
    u_rate2 = {}
    p = {}
    param = {}
    f = {}

    # Build the interpolation matrix which interpolates from all nodes on the old grid to
    # all nodes on the new grid.
    grid_data = transcription.grid_data
    L, _ = interpolation_lagrange_matrix(old_grid=phase.options['transcription'].grid_data,
                                         new_grid=grid_data)

    # Create a new problem for the grid_refinement
    # For this test, use the same grid as the original problem.
    p_refine = om.Problem(model=om.Group())
    grid_refinement_system = GridRefinementODESystem(grid_data=grid_data,
                                                     time=phase.time_options,
                                                     states=phase.state_options,
                                                     controls=phase.control_options,
                                                     parameters=phase.parameter_options,
                                                     ode_class=phase.options['ode_class'],
                                                     ode_init_kwargs=phase.options[
                                                         'ode_init_kwargs'],
                                                     calc_exprs=phase._calc_exprs)
    p_refine.model.add_subsystem('grid_refinement_system', grid_refinement_system, promotes=['*'])
    p_refine.setup()

    # Set the values in the refinement problem using the outputs from the first
    ode = p_refine.model.grid_refinement_system.ode

    t_name = phase.time_options['name']

    t_prev = phase.get_val(f'timeseries.{t_name}', units=phase.time_options['units'])
    t_phase_prev = t_prev - t_prev[0]
    t_initial = np.repeat(t_prev[0, 0], repeats=transcription.grid_data.num_nodes, axis=0)
    t_duration = np.repeat(t_prev[-1, 0], repeats=transcription.grid_data.num_nodes, axis=0)
    t = np.dot(L, t_prev)
    t_phase = np.dot(L, t_phase_prev)
    targets = get_targets(ode, 'time', phase.time_options['targets'])
    t_phase_targets = get_targets(ode, 't_phase', phase.time_options['time_phase_targets'])
    t_initial_targets = get_targets(ode, 't_initial', phase.time_options['t_initial_targets'])
    t_duration_targets = get_targets(ode, 't_duration', phase.time_options['t_duration_targets'])
    if targets:
        p_refine.set_val('time', t)
    if t_phase_targets:
        p_refine.set_val('t_phase', t_phase)
    if t_initial_targets:
        p_refine.set_val('t_initial', t_initial)
    if t_duration_targets:
        p_refine.set_val('t_duration', t_duration)

    state_prefix = 'states:' if phase.timeseries_options['use_prefix'] else ''
    control_prefix = 'controls:' if phase.timeseries_options['use_prefix'] else ''

    for name, options in phase.state_options.items():
        x_prev = phase.get_val(f'timeseries.{state_prefix}{name}', units=options['units'])
        x[name] = np.dot(L, x_prev)
        targets = get_targets(ode, name, options['targets'])
        if targets:
            p_refine.set_val(f'states:{name}', x[name])

    for name, options in phase.control_options.items():
        targets = get_targets(ode, name, options['targets'])
        rate_targets = get_targets(ode, f'{name}_rate', options['rate_targets'])
        rate2_targets = get_targets(ode, f'{name}_rate2', options['rate2_targets'])

        u_prev = phase.get_val(f'timeseries.{control_prefix}{name}', units=options['units'])
        u[name] = np.dot(L, u_prev)
        if targets:
            p_refine.set_val(f'controls:{name}', u[name])

        if phase.timeseries_options['include_control_rates']:
            if rate_targets:
                u_rate_prev = phase.get_val(f'timeseries.control_rates:{name}_rate')
                u_rate[name] = np.dot(L, u_rate_prev)
                p_refine.set_val(f'control_rates:{name}_rate', u_rate[name])

            if rate2_targets:
                u_rate2_prev = phase.get_val(f'timeseries.control_rates:{name}_rate2')
                u_rate2[name] = np.dot(L, u_rate2_prev)
                p_refine.set_val(f'control_rates:{name}_rate2', u_rate2[name])

    # Configure the parameters
    for name, options in phase.parameter_options.items():
        targets = get_targets(ode, name, options['targets'])
        # The value of the parameter at one node
        param[name] = phase.get_val(f'parameter_vals:{name}', units=options['units'])[0, ...]
        if targets:
            p_refine.set_val(f'parameters:{name}', param[name], units=options['units'])

    # Execute the model
    p_refine.run_model()

    # Assign the state rates on the new grid to f
    for name, options in phase.state_options.items():
        rate_units = get_rate_units(options['units'], phase.time_options['units'])
        rate_source = options['rate_source']
        rate_source_class = phase.classify_var(rate_source)
        if rate_source_class in {'time'}:
            src_units = phase.time_options['units']
            f[name] = om.convert_units(t, src_units, rate_units)
        elif rate_source_class in {'time_phase'}:
            src_units = phase.time_options['units']
            f[name] = om.convert_units(t_phase, src_units, rate_units)
        elif rate_source_class in {'state'}:
            src_units = phase.state_options[rate_source]['units']
            f[name] = om.convert_units(x[rate_source], src_units, rate_units)
        elif rate_source_class in {'control'}:
            src_units = phase.control_options[rate_source]['units']
            f[name] = om.convert_units(u[rate_source], src_units, rate_units)
        elif rate_source_class in {'control_rate'}:
            u_name = rate_source[:-5]
            u_units = phase.control_options[u_name]['units']
            src_units = get_rate_units(u_units, phase.time_options['units'], deriv=1)
            f[name] = om.convert_units(u_rate[rate_source], src_units, rate_units)
        elif rate_source_class in {'control_rate2'}:
            u_name = rate_source[:-6]
            u_units = phase.control_options[u_name]['units']
            src_units = get_rate_units(u_units, phase.time_options['units'], deriv=2)
            f[name] = om.convert_units(u_rate2[rate_source], src_units, rate_units)
        elif rate_source_class in {'parameter'}:
            src_units = phase.parameter_options[rate_source]['units']
            shape = phase.parameter_options[rate_source]['shape']
            f[name] = om.convert_units(param[rate_source], src_units, rate_units)
            f[name] = np.broadcast_to(f[name], (grid_data.num_nodes,) + shape)
        elif rate_source_class in {'ode'}:
            f[name] = p_refine.get_val(f'ode.{rate_source}', units=rate_units)

        if len(f[name].shape) == 1:
            f[name] = np.reshape(f[name], (f[name].shape[0], 1))

    return x, u, p, f


def compute_state_quadratures(x_hat, f_hat, t_duration, transcription):
    """
    Compute the integral of the given states at each node in the given transcription.

    This estimation of the states is computed with a quadrature.

    Parameters
    ----------
    x_hat : dict of (str: float)
        The interpolated state values at the nodes for each state.
    f_hat : dict of (str: float)
        State rates computed using the interpolated state values at each node.
    t_duration : float
        The time duration of the phase.
    transcription : Radau or GaussLobatto
        The transcription instance used to compute f_hat.

    Returns
    -------
    dict
        A dictionary keyed by state name containing the estimated state values at each node,
        computed using a quadrature.
    """
    x_prime = {}
    gd = transcription.grid_data

    # Build the integration matrix which integrates state values at all nodes on the new grid.
    I = integration_matrix(gd)  # noqa: E741, allow ambiguous variable name `I`

    left_end_idxs = gd.subset_node_indices['segment_ends'][0::2]
    all_idxs = gd.subset_node_indices['all']
    not_left_end_idxs = np.array(sorted(set(all_idxs).difference(left_end_idxs)))

    dt_dstau = np.atleast_2d(0.5 * t_duration * gd.node_dptau_dstau[not_left_end_idxs]).T

    if x_hat.keys() != f_hat.keys():
        raise ValueError('x_hat and f_hat don\'t contain the same states.\n'
                         f'x_hat states are: {list(x_hat.keys())}\n'
                         f'f_hat states are: {list(f_hat.keys())}')

    for state_name in x_hat:
        x_prime[state_name] = np.zeros_like(x_hat[state_name])
        x_prime[state_name][left_end_idxs, ...] = x_hat[state_name][left_end_idxs, ...]
        nnps = np.array(gd.subset_num_nodes_per_segment['all']) - 1
        left_end_idxs_repeated = np.repeat(left_end_idxs, nnps)
        x_prime[state_name][not_left_end_idxs, ...] = \
            x_hat[state_name][left_end_idxs_repeated, ...] \
            + dt_dstau * np.dot(I, f_hat[state_name][not_left_end_idxs, ...])

    return x_prime


def check_error(phases):
    """
    Compute the error in every solved segment of the given phases.

    Parameters
    ----------
    phases : dict
        Dict of phase paths and phases.

    Returns
    -------
    dict
        Indicator for which segments of each phase require grid refinement.
    """
    refine_results = {}

    for phase_path, phase in phases.items():
        refine_results[phase_path] = {}

        # Save the original grid to the refine results
        tx = phase.options['transcription']
        gd = tx.grid_data
        numseg = gd.num_segments

        refine_results[phase_path]['num_segments'] = numseg
        refine_results[phase_path]['order'] = gd.transcription_order
        refine_results[phase_path]['segment_ends'] = gd.segment_ends
        refine_results[phase_path]['need_refinement'] = np.zeros(numseg, dtype=bool)
        refine_results[phase_path]['max_rel_error'] = np.zeros(numseg, dtype=float)  # Eq. 21
        refine_results[phase_path]['error_state'] = ['' for _ in phase.state_options]

        # Instantiate a new phase as a copy of the old one, but first up the transcription order
        # by 1 for Radau and by 2 for Gauss-Lobatto
        new_num_segments = tx.options['num_segments']
        new_segment_ends = tx.options['segment_ends']
        new_compressed = tx.options['compressed']
        if isinstance(tx, GaussLobatto):
            new_order = tx.options['order'] + 2
            new_tx = GaussLobatto(num_segments=new_num_segments, order=new_order,
                                  segment_ends=new_segment_ends, compressed=new_compressed)
        elif isinstance(tx, Radau):
            new_order = tx.options['order'] + 1
            new_tx = Radau(num_segments=new_num_segments, order=new_order,
                           segment_ends=new_segment_ends, compressed=new_compressed)
        else:
            # Only refine GuassLobatto or Radau transcription phases
            continue

        # Let x be the interpolated states on the new transcription
        # Let f by the evaluated state rates given the interpolation of the states and controls
        # onto the new grid.
        x, _, _, f = eval_ode_on_grid(phase=phase, transcription=new_tx)

        # x_hat is the state value at each node computed using a quadrature
        # from the initial state value in each segment and the computed state rates
        # at each node in the new transcription.
        x_hat = compute_state_quadratures(x, f, phase.get_val('t_duration'), new_tx)

        E = {}  # The absolute error computed in each state at each node (Eq. 20 pt 1)
        e = {}  # The relative error computed in each state at each node (Eq. 20 pt 2)

        for state_name in phase.state_options:
            E[state_name] = np.abs(x_hat[state_name] - x[state_name])  # Equation 20.1
            e[state_name] = np.zeros_like(E[state_name])               # Equation 20.2
            for k in range(numseg):
                i1, i2 = new_tx.grid_data.subset_segment_indices['all'][k, :]
                k_idxs = new_tx.grid_data.subset_node_indices['all'][i1:i2]
                e[state_name][k_idxs, ...] = E[state_name][k_idxs] \
                    / (1.0 + np.max(np.abs(x[state_name][k_idxs])))
                if np.any(np.max(e[state_name][k_idxs]) > refine_results[phase_path]['max_rel_error'][k]):
                    refine_results[phase_path]['max_rel_error'][k] = np.max(e[state_name][k_idxs])
                    refine_results[phase_path]['error_state'] = state_name
                    if refine_results[phase_path]['max_rel_error'][k] > phase.refine_options['tolerance']:
                        refine_results[phase_path]['need_refinement'][k] = True

    return refine_results

1	"""
2	Generic utilities for use by the grid refinement schemes.
3	"""
4	import numpy as np	5✔
5	from scipy.linalg import block_diag	5✔
6
7	import openmdao.api as om	5✔
8
9	from ..transcriptions import GaussLobatto, Radau	5✔
10	from ..utils.lagrange import lagrange_matrices	5✔
11	from ..utils.misc import get_rate_units	5✔
12	from ..utils.introspection import get_targets	5✔
13
14	from .grid_refinement_ode_system import GridRefinementODESystem	5✔
15
16
17	def interpolation_lagrange_matrix(old_grid, new_grid):	5✔
18	"""
19	Evaluate lagrange matrix to interpolate state and control values from the solved grid onto the new grid.
20
21	Parameters
22	----------
23	old_grid : <GridData>
24	GridData object representing the grid on which the problem has been solved.
25	new_grid : <GridData>
26	GridData object representing the new, higher-order grid.
27
28	Returns
29	-------
30	ndarray
31	The lagrange interpolation matrix.
32	"""
33	L_blocks = []	7✔
34	D_blocks = []	7✔
35
36	for iseg in range(old_grid.num_segments):	7✔
37	i1, i2 = old_grid.subset_segment_indices['all'][iseg, :]	7✔
38	indices = old_grid.subset_node_indices['all'][i1:i2]	7✔
39	nodes_given = old_grid.node_stau[indices]	7✔
40
41	i1, i2 = new_grid.subset_segment_indices['all'][iseg, :]	7✔
42	indices = new_grid.subset_node_indices['all'][i1:i2]	7✔
43	nodes_eval = new_grid.node_stau[indices]	7✔
44
45	L_block, D_block = lagrange_matrices(nodes_given, nodes_eval)	7✔
46
47	L_blocks.append(L_block)	7✔
48	D_blocks.append(D_block)	7✔
49
50	L = block_diag(*L_blocks)	7✔
51	D = block_diag(*D_blocks)	7✔
52
53	return L, D	7✔
54
55
56	def integration_matrix(grid):	5✔
57	"""
58	Evaluate the Integration matrix of the given grid.
59
60	Parameters
61	----------
62	grid : <GridData>
63	GridData object containing the grid where the integration matrix will be evaluated.
64
65	Returns
66	-------
67	ndarray
68	The integration matrix used to propagate initial states over segments.
69	"""
70	I_blocks = []	7✔
71
72	for iseg in range(grid.num_segments):	7✔
73	i1, i2 = grid.subset_segment_indices['all'][iseg, :]	7✔
74	indices = grid.subset_node_indices['all'][i1:i2]	7✔
75	nodes_given = grid.node_stau[indices]	7✔
76
77	i1, i2 = grid.subset_segment_indices['all'][iseg, :]	7✔
78	indices = grid.subset_node_indices['all'][i1:i2]	7✔
79	nodes_eval = grid.node_stau[indices][1:]	7✔
80
81	_, D_block = lagrange_matrices(nodes_given, nodes_eval)	7✔
82	I_block = np.linalg.inv(D_block[:, 1:])	7✔
83	I_blocks.append(I_block)	7✔
84
85	return block_diag(*I_blocks)	7✔
86
87
88	def eval_ode_on_grid(phase, transcription):	5✔
89	"""
90	Evaluate the ODE from the given phase at all nodes of the given transcription.
91
92	Parameters
93	----------
94	phase : Phase
95	A Phase object which has been executed and whose ODE is to be run at all nodes
96	of the given transcription.
97	transcription : Radau or GaussLobatto transcription instance
98	The transcription object at which to execute the ODE of the given phase at all nodes.
99
100	Returns
101	-------
102	dict of (str: ndarray)
103	A dictionary of the state values from the phase interpolated to the new transcription.
104	dict of (str: ndarray)
105	A dictionary of the control values from the phase interpolated to the new transcription.
106	dict of (str: ndarray)
107	A dictionary of the state rates computed in the phase's ODE at the new transcription points.
108	"""
109	x = {}	7✔
110	u = {}	7✔
111	u_rate = {}	7✔
112	u_rate2 = {}	7✔
113	p = {}	7✔
114	param = {}	7✔
115	f = {}	7✔
116
117	# Build the interpolation matrix which interpolates from all nodes on the old grid to
118	# all nodes on the new grid.
119	grid_data = transcription.grid_data	7✔
120	L, _ = interpolation_lagrange_matrix(old_grid=phase.options['transcription'].grid_data,	7✔
121	new_grid=grid_data)
122
123	# Create a new problem for the grid_refinement
124	# For this test, use the same grid as the original problem.
125	p_refine = om.Problem(model=om.Group())	7✔
126	grid_refinement_system = GridRefinementODESystem(grid_data=grid_data,	7✔
127	time=phase.time_options,
128	states=phase.state_options,
129	controls=phase.control_options,
130	parameters=phase.parameter_options,
131	ode_class=phase.options['ode_class'],
132	ode_init_kwargs=phase.options[
133	'ode_init_kwargs'],
134	calc_exprs=phase._calc_exprs)
135	p_refine.model.add_subsystem('grid_refinement_system', grid_refinement_system, promotes=['*'])	7✔
136	p_refine.setup()	7✔
137
138	# Set the values in the refinement problem using the outputs from the first
139	ode = p_refine.model.grid_refinement_system.ode	7✔
140
141	t_name = phase.time_options['name']	7✔
142
143	t_prev = phase.get_val(f'timeseries.{t_name}', units=phase.time_options['units'])	7✔
144	t_phase_prev = t_prev - t_prev[0]	7✔
145	t_initial = np.repeat(t_prev[0, 0], repeats=transcription.grid_data.num_nodes, axis=0)	7✔
146	t_duration = np.repeat(t_prev[-1, 0], repeats=transcription.grid_data.num_nodes, axis=0)	7✔
147	t = np.dot(L, t_prev)	7✔
148	t_phase = np.dot(L, t_phase_prev)	7✔
149	targets = get_targets(ode, 'time', phase.time_options['targets'])	7✔
150	t_phase_targets = get_targets(ode, 't_phase', phase.time_options['time_phase_targets'])	7✔
151	t_initial_targets = get_targets(ode, 't_initial', phase.time_options['t_initial_targets'])	7✔
152	t_duration_targets = get_targets(ode, 't_duration', phase.time_options['t_duration_targets'])	7✔
153	if targets:	7✔
UNCOV 154	p_refine.set_val('time', t)	×
155	if t_phase_targets:	7✔
UNCOV 156	p_refine.set_val('t_phase', t_phase)	×
157	if t_initial_targets:	7✔
UNCOV 158	p_refine.set_val('t_initial', t_initial)	×
159	if t_duration_targets:	7✔
UNCOV 160	p_refine.set_val('t_duration', t_duration)	×
161
162	state_prefix = 'states:' if phase.timeseries_options['use_prefix'] else ''	7✔
163	control_prefix = 'controls:' if phase.timeseries_options['use_prefix'] else ''	7✔
164
165	for name, options in phase.state_options.items():	7✔
166	x_prev = phase.get_val(f'timeseries.{state_prefix}{name}', units=options['units'])	7✔
167	x[name] = np.dot(L, x_prev)	7✔
168	targets = get_targets(ode, name, options['targets'])	7✔
169	if targets:	7✔
170	p_refine.set_val(f'states:{name}', x[name])	7✔
171
172	for name, options in phase.control_options.items():	7✔
173	targets = get_targets(ode, name, options['targets'])	7✔
174	rate_targets = get_targets(ode, f'{name}_rate', options['rate_targets'])	7✔
175	rate2_targets = get_targets(ode, f'{name}_rate2', options['rate2_targets'])	7✔
176
177	u_prev = phase.get_val(f'timeseries.{control_prefix}{name}', units=options['units'])	7✔
178	u[name] = np.dot(L, u_prev)	7✔
179	if targets:	7✔
180	p_refine.set_val(f'controls:{name}', u[name])	7✔
181
182	if phase.timeseries_options['include_control_rates']:	7✔
183	if rate_targets:	×
184	u_rate_prev = phase.get_val(f'timeseries.control_rates:{name}_rate')	×
185	u_rate[name] = np.dot(L, u_rate_prev)	×
UNCOV 186	p_refine.set_val(f'control_rates:{name}_rate', u_rate[name])	×
187
188	if rate2_targets:	×
189	u_rate2_prev = phase.get_val(f'timeseries.control_rates:{name}_rate2')	×
190	u_rate2[name] = np.dot(L, u_rate2_prev)	×
UNCOV 191	p_refine.set_val(f'control_rates:{name}_rate2', u_rate2[name])	×
192
193	# Configure the parameters
194	for name, options in phase.parameter_options.items():	7✔
195	targets = get_targets(ode, name, options['targets'])	7✔
196	# The value of the parameter at one node
197	param[name] = phase.get_val(f'parameter_vals:{name}', units=options['units'])[0, ...]	7✔
198	if targets:	7✔
199	p_refine.set_val(f'parameters:{name}', param[name], units=options['units'])	7✔
200
201	# Execute the model
202	p_refine.run_model()	7✔
203
204	# Assign the state rates on the new grid to f
205	for name, options in phase.state_options.items():	7✔
206	rate_units = get_rate_units(options['units'], phase.time_options['units'])	7✔
207	rate_source = options['rate_source']	7✔
208	rate_source_class = phase.classify_var(rate_source)	7✔
209	if rate_source_class in {'time'}:	7✔
210	src_units = phase.time_options['units']	×
UNCOV 211	f[name] = om.convert_units(t, src_units, rate_units)	×
212	elif rate_source_class in {'time_phase'}:	7✔
213	src_units = phase.time_options['units']	×
UNCOV 214	f[name] = om.convert_units(t_phase, src_units, rate_units)	×
215	elif rate_source_class in {'state'}:	7✔
216	src_units = phase.state_options[rate_source]['units']	7✔
217	f[name] = om.convert_units(x[rate_source], src_units, rate_units)	7✔
218	elif rate_source_class in {'control'}:	7✔
219	src_units = phase.control_options[rate_source]['units']	7✔
220	f[name] = om.convert_units(u[rate_source], src_units, rate_units)	7✔
221	elif rate_source_class in {'control_rate'}:	7✔
222	u_name = rate_source[:-5]	×
223	u_units = phase.control_options[u_name]['units']	×
224	src_units = get_rate_units(u_units, phase.time_options['units'], deriv=1)	×
UNCOV 225	f[name] = om.convert_units(u_rate[rate_source], src_units, rate_units)	×
226	elif rate_source_class in {'control_rate2'}:	7✔
227	u_name = rate_source[:-6]	×
228	u_units = phase.control_options[u_name]['units']	×
229	src_units = get_rate_units(u_units, phase.time_options['units'], deriv=2)	×
UNCOV 230	f[name] = om.convert_units(u_rate2[rate_source], src_units, rate_units)	×
231	elif rate_source_class in {'parameter'}:	7✔
232	src_units = phase.parameter_options[rate_source]['units']	4✔
233	shape = phase.parameter_options[rate_source]['shape']	4✔
234	f[name] = om.convert_units(param[rate_source], src_units, rate_units)	4✔
235	f[name] = np.broadcast_to(f[name], (grid_data.num_nodes,) + shape)	4✔
236	elif rate_source_class in {'ode'}:	7✔
237	f[name] = p_refine.get_val(f'ode.{rate_source}', units=rate_units)	7✔
238
239	if len(f[name].shape) == 1:	7✔
240	f[name] = np.reshape(f[name], (f[name].shape[0], 1))	7✔
241
242	return x, u, p, f	7✔
243
244
245	def compute_state_quadratures(x_hat, f_hat, t_duration, transcription):	5✔
246	"""
247	Compute the integral of the given states at each node in the given transcription.
248
249	This estimation of the states is computed with a quadrature.
250
251	Parameters
252	----------
253	x_hat : dict of (str: float)
254	The interpolated state values at the nodes for each state.
255	f_hat : dict of (str: float)
256	State rates computed using the interpolated state values at each node.
257	t_duration : float
258	The time duration of the phase.
259	transcription : Radau or GaussLobatto
260	The transcription instance used to compute f_hat.
261
262	Returns
263	-------
264	dict
265	A dictionary keyed by state name containing the estimated state values at each node,
266	computed using a quadrature.
267	"""
268	x_prime = {}	7✔
269	gd = transcription.grid_data	7✔
270
271	# Build the integration matrix which integrates state values at all nodes on the new grid.
272	I = integration_matrix(gd) # noqa: E741, allow ambiguous variable name `I`	7✔
273
274	left_end_idxs = gd.subset_node_indices['segment_ends'][0::2]	7✔
275	all_idxs = gd.subset_node_indices['all']	7✔
276	not_left_end_idxs = np.array(sorted(set(all_idxs).difference(left_end_idxs)))	7✔
277
278	dt_dstau = np.atleast_2d(0.5 * t_duration * gd.node_dptau_dstau[not_left_end_idxs]).T	7✔
279
280	if x_hat.keys() != f_hat.keys():	7✔
UNCOV 281	raise ValueError('x_hat and f_hat don\'t contain the same states.\n'	×
282	f'x_hat states are: {list(x_hat.keys())}\n'
283	f'f_hat states are: {list(f_hat.keys())}')
284
285	for state_name in x_hat:	7✔
286	x_prime[state_name] = np.zeros_like(x_hat[state_name])	7✔
287	x_prime[state_name][left_end_idxs, ...] = x_hat[state_name][left_end_idxs, ...]	7✔
288	nnps = np.array(gd.subset_num_nodes_per_segment['all']) - 1	7✔
289	left_end_idxs_repeated = np.repeat(left_end_idxs, nnps)	7✔
290	x_prime[state_name][not_left_end_idxs, ...] = \	7✔
291	x_hat[state_name][left_end_idxs_repeated, ...] \
292	+ dt_dstau * np.dot(I, f_hat[state_name][not_left_end_idxs, ...])
293
294	return x_prime	7✔
295
296
297	def check_error(phases):	5✔
298	"""
299	Compute the error in every solved segment of the given phases.
300
301	Parameters
302	----------
303	phases : dict
304	Dict of phase paths and phases.
305
306	Returns
307	-------
308	dict
309	Indicator for which segments of each phase require grid refinement.
310	"""
311	refine_results = {}	7✔
312
313	for phase_path, phase in phases.items():	7✔
314	refine_results[phase_path] = {}	7✔
315
316	# Save the original grid to the refine results
317	tx = phase.options['transcription']	7✔
318	gd = tx.grid_data	7✔
319	numseg = gd.num_segments	7✔
320
321	refine_results[phase_path]['num_segments'] = numseg	7✔
322	refine_results[phase_path]['order'] = gd.transcription_order	7✔
323	refine_results[phase_path]['segment_ends'] = gd.segment_ends	7✔
324	refine_results[phase_path]['need_refinement'] = np.zeros(numseg, dtype=bool)	7✔
325	refine_results[phase_path]['max_rel_error'] = np.zeros(numseg, dtype=float) # Eq. 21	7✔
326	refine_results[phase_path]['error_state'] = ['' for _ in phase.state_options]	7✔
327
328	# Instantiate a new phase as a copy of the old one, but first up the transcription order
329	# by 1 for Radau and by 2 for Gauss-Lobatto
330	new_num_segments = tx.options['num_segments']	7✔
331	new_segment_ends = tx.options['segment_ends']	7✔
332	new_compressed = tx.options['compressed']	7✔
333	if isinstance(tx, GaussLobatto):	7✔
334	new_order = tx.options['order'] + 2	7✔
335	new_tx = GaussLobatto(num_segments=new_num_segments, order=new_order,	7✔
336	segment_ends=new_segment_ends, compressed=new_compressed)
337	elif isinstance(tx, Radau):	7✔
338	new_order = tx.options['order'] + 1	7✔
339	new_tx = Radau(num_segments=new_num_segments, order=new_order,	7✔
340	segment_ends=new_segment_ends, compressed=new_compressed)
341	else:
342	# Only refine GuassLobatto or Radau transcription phases
UNCOV 343	continue	×
344
345	# Let x be the interpolated states on the new transcription
346	# Let f by the evaluated state rates given the interpolation of the states and controls
347	# onto the new grid.
348	x, _, _, f = eval_ode_on_grid(phase=phase, transcription=new_tx)	7✔
349
350	# x_hat is the state value at each node computed using a quadrature
351	# from the initial state value in each segment and the computed state rates
352	# at each node in the new transcription.
353	x_hat = compute_state_quadratures(x, f, phase.get_val('t_duration'), new_tx)	7✔
354
355	E = {} # The absolute error computed in each state at each node (Eq. 20 pt 1)	7✔
356	e = {} # The relative error computed in each state at each node (Eq. 20 pt 2)	7✔
357
358	for state_name in phase.state_options:	7✔
359	E[state_name] = np.abs(x_hat[state_name] - x[state_name]) # Equation 20.1	7✔
360	e[state_name] = np.zeros_like(E[state_name]) # Equation 20.2	7✔
361	for k in range(numseg):	7✔
362	i1, i2 = new_tx.grid_data.subset_segment_indices['all'][k, :]	7✔
363	k_idxs = new_tx.grid_data.subset_node_indices['all'][i1:i2]	7✔
364	e[state_name][k_idxs, ...] = E[state_name][k_idxs] \	7✔
365	/ (1.0 + np.max(np.abs(x[state_name][k_idxs])))
366	if np.any(np.max(e[state_name][k_idxs]) > refine_results[phase_path]['max_rel_error'][k]):	7✔
367	refine_results[phase_path]['max_rel_error'][k] = np.max(e[state_name][k_idxs])	7✔
368	refine_results[phase_path]['error_state'] = state_name	7✔
369	if refine_results[phase_path]['max_rel_error'][k] > phase.refine_options['tolerance']:	7✔
370	refine_results[phase_path]['need_refinement'][k] = True	7✔
371
372	return refine_results	7✔

swryan / dymos / 14268934115

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous