The Van der Pol Oscillator#
In dynamics, the Van Der Pol oscillator [Wikipediacontributors20] is a non-conservative oscillator with non-linear damping. It evolves in time according to the second-order differential equation:
(16)#\[\begin{align}
\frac{d^2x}{dt^2} - u (1 - x^2) \frac{dx}{dt} + x &= 0
\end{align}\]
where \(x\) is the position coordinate (a function of the time \(t\)), and \(u\) is a scalar parameter indicating the nonlinearity and the strength of the damping.
To make this an optimal control problem, we want to find the smallest control that will dampen the oscillation (drive the state variables to zero). We can express this as an objective function \(J\) to minimize:
(17)#\[\begin{align}
J &= \int x^2_0 + x^2_1 + u^2
\end{align}\]
In other words, we want to find the optimal (smallest) trajectory of the control \(u\) such that the oscillation and the oscillation’s rate of change are driven to zero.
State Variables#
There are three state variables are used to define the configuration of the system at any given instant in time.
\(x_1\): The primary output of the oscillator.
\(x_0\): The rate of change of the primary output.
\(J\): The objective function to be minimized.
The objective function is included as a state variable so that Dymos will do the integration.
The \(x_1\) and \(x_0\) state variables are also inputs to the system, along with the control \(u\).
System Dynamics#
The evolution of the state variables is given by the following ordinary differential equations (ODE):
(18)#\[\begin{align}
\frac{dx_0}{dt} &= (1 - x^2_1) x_0 - x_1 + u \\
\frac{dx_1}{dt} &= x_0 \\
\frac{dJ}{dt} &= x^2_0 + x^2_1 + u^2
\end{align}\]
Control Variables#
This system has a single control variable:
\(u\): The control input.
The control variable has a constraint: \(-0.75 \leq u \leq 1.0\)
The initial and final conditions#
The initial conditions are:
(19)#\[\begin{align}
t_i &= 0 \\
x_{0i} &= 1 \\
x_{1i} &= 1 \\
u_i &= -0.75
\end{align}\]
The final conditions are:
(20)#\[\begin{align}
t_f &= 10 \\
x_{0f} &= 0 \\
x_{1f} &= 0 \\
u_f &= 0
\end{align}\]
Defining the ODE as an OpenMDAO System#
In Dymos, the ODE is an OpenMDAO System (a Component, or a Group of components). The following ExplicitComponent computes the state rates for the Van der Pol problem.
More detail on the workings of an ExplicitComponent can be found in the OpenMDAO documentation. In summary:
initialize: Called at setup, and used to define options for the component. ALL Dymos ODE components should have the property
num_nodes, which defines the number of points at which the outputs are simultaneously computed.setup: Used to add inputs and outputs to the component, and declare which outputs (and indices of outputs) are dependent on each of the inputs.
compute: Used to compute the outputs, given the inputs.
compute_partials: Used to compute the derivatives of the outputs with respect to each of the inputs analytically. This method may be omitted if finite difference or complex-step approximations are used, though analytic is recommended.
Note
Things to note about the Van der Pol ODE system
Only the vanderpol_ode class below is important for defining the basic problem. The other classes are used to demonstrate Message Passing Interface (MPI) parallel calculation of the system. They can be ignored.
\(x_1\), \(x_0\), and \(u\) are inputs.
\(\dot{x_1}\), \(\dot{x_0}\), and \(\dot{J}\) are outputs.
declare_partials is called for every output with respect to every input.
For efficiency, partial derrivatives that are constant have values specified in the setup method rather than the compute_partials method. So although 7 partials are declared, only 5 are computed in compute_partials.
This ODE includes some additional complexity to allow it to run in parallel. See the Van der Pol example in the dymos examples directory for an example of how to run the problem under MPI.
import numpy as np
import openmdao.api as om
import time
from openmdao.utils.array_utils import evenly_distrib_idxs
class VanderpolODE(om.ExplicitComponent):
"""intentionally slow version of vanderpol_ode for effects of demonstrating distributed component calculations
MPI can run this component in multiple processes, distributing the calculation of derivatives.
This code has a delay in it to simulate a longer computation. It should run faster with more processes.
"""
def __init__(self, *args, **kwargs):
self.progress_prints = False
super().__init__(*args, **kwargs)
def initialize(self):
self.options.declare('num_nodes', types=int)
self.options.declare('distrib', types=bool, default=False)
def setup(self):
nn = self.options['num_nodes']
comm = self.comm
rank = comm.rank
sizes, offsets = evenly_distrib_idxs(comm.size, nn) # (#cpus, #inputs) -> (size array, offset array)
self.start_idx = offsets[rank]
self.io_size = sizes[rank] # number of inputs and outputs managed by this distributed process
self.end_idx = self.start_idx + self.io_size
# inputs: 2 states and a control
self.add_input('x0', val=np.ones(nn), desc='derivative of Output', units='V/s')
self.add_input('x1', val=np.ones(nn), desc='Output', units='V')
self.add_input('u', val=np.ones(nn), desc='control', units=None)
# outputs: derivative of states
# the objective function will be treated as a state for computation, so its derivative is an output
self.add_output('x0dot', val=np.ones(self.io_size), desc='second derivative of Output',
units='V/s**2', distributed=self.options['distrib'])
self.add_output('x1dot', val=np.ones(self.io_size), desc='derivative of Output',
units='V/s', distributed=self.options['distrib'])
self.add_output('Jdot', val=np.ones(self.io_size), desc='derivative of objective',
units='1.0/s', distributed=self.options['distrib'])
# self.declare_coloring(method='cs')
# # partials
r = np.arange(self.io_size, dtype=int)
c = r + self.start_idx
self.declare_partials(of='x0dot', wrt='x0', rows=r, cols=c)
self.declare_partials(of='x0dot', wrt='x1', rows=r, cols=c)
self.declare_partials(of='x0dot', wrt='u', rows=r, cols=c, val=1.0)
self.declare_partials(of='x1dot', wrt='x0', rows=r, cols=c, val=1.0)
self.declare_partials(of='Jdot', wrt='x0', rows=r, cols=c)
self.declare_partials(of='Jdot', wrt='x1', rows=r, cols=c)
self.declare_partials(of='Jdot', wrt='u', rows=r, cols=c)
def compute(self, inputs, outputs):
# The inputs contain the entire vector, be each rank will only operate on a portion of it.
x0 = inputs['x0'][self.start_idx:self.end_idx]
x1 = inputs['x1'][self.start_idx:self.end_idx]
u = inputs['u'][self.start_idx:self.end_idx]
outputs['x0dot'] = (1.0 - x1**2) * x0 - x1 + u
outputs['x1dot'] = x0
outputs['Jdot'] = x0**2 + x1**2 + u**2
def compute_partials(self, inputs, jacobian):
x0 = inputs['x0'][self.start_idx:self.end_idx]
x1 = inputs['x1'][self.start_idx:self.end_idx]
u = inputs['u'][self.start_idx:self.end_idx]
jacobian['x0dot', 'x0'] = 1.0 - x1 * x1
jacobian['x0dot', 'x1'] = -2.0 * x1 * x0 - 1.0
jacobian['Jdot', 'x0'] = 2.0 * x0
jacobian['Jdot', 'x1'] = 2.0 * x1
jacobian['Jdot', 'u'] = 2.0 * u
Defining the Dymos Problem#
Once the ODEs are defined, they are used to create a Dymos Problem object that allows solution.
Note
Things to note about the Van der Pol Dymos Problem definition
The vanderpol function creates and returns a Dymos Problem instance that can be used for simulation or optimization.
The vanderpol function has optional arguments for specifying options for the type of transcription, number of segments, optimizer, etc. These can be ignored when first trying to understand the code.
The Problem object has a Trajectory object, and the trajectory has a single Phase. Most of the problem setup is performed by calling methods on the phase (set_time_options, add_state, add_boundary_constraint, add_objective).
The add_state and add_control calls include the target parameter for \(x_0\), \(x_1\), and \(u\). This is required so that the inputs are correctly calculated.
Initial (linear) guesses are supplied for the states and control.
import openmdao.api as om
import dymos as dm
def vanderpol(transcription='gauss-lobatto', num_segments=15, transcription_order=3,
compressed=True, optimizer='SLSQP', use_pyoptsparse=False):
"""Dymos problem definition for optimal control of a Van der Pol oscillator"""
# define the OpenMDAO problem
p = om.Problem(model=om.Group())
if not use_pyoptsparse:
p.driver = om.ScipyOptimizeDriver()
else:
p.driver = om.pyOptSparseDriver(print_results=False)
p.driver.options['optimizer'] = optimizer
if use_pyoptsparse:
if optimizer == 'SNOPT':
p.driver.opt_settings['iSumm'] = 6 # show detailed SNOPT output
elif optimizer == 'IPOPT':
p.driver.opt_settings['print_level'] = 4
p.driver.declare_coloring()
# define a Trajectory object and add to model
traj = dm.Trajectory()
p.model.add_subsystem('traj', subsys=traj)
# define a Transcription
if transcription == 'gauss-lobatto':
t = dm.GaussLobatto(num_segments=num_segments,
order=transcription_order,
compressed=compressed)
elif transcription == 'radau-ps':
t = dm.Radau(num_segments=num_segments,
order=transcription_order,
compressed=compressed)
# define a Phase as specified above and add to Phase
phase = dm.Phase(ode_class=VanderpolODE, transcription=t)
traj.add_phase(name='phase0', phase=phase)
t_final = 10
phase.set_time_options(fix_initial=True, fix_duration=True, duration_val=t_final, units='s')
# set the State time options
phase.add_state('x0', fix_initial=True, fix_final=True,
rate_source='x0dot',
units='V/s', ref=0.1, defect_ref=0.1) # target required because x0 is an input
phase.add_state('x1', fix_initial=True, fix_final=True,
rate_source='x1dot',
units='V', ref=0.1, defect_ref=0.1)
phase.add_state('J', fix_initial=True, fix_final=False,
rate_source='Jdot',
units=None)
# define the control
phase.add_control(name='u', units=None, lower=-0.75, upper=1.0, continuity=True,
rate_continuity=True)
# define objective to minimize
phase.add_objective('J', loc='final')
# setup the problem
p.setup(check=True)
p['traj.phase0.t_initial'] = 0.0
p['traj.phase0.t_duration'] = t_final
# add a linearly interpolated initial guess for the state and control curves
p['traj.phase0.states:x0'] = phase.interp('x0', [0, 0])
p['traj.phase0.states:x1'] = phase.interp('x1', [1, 0])
p['traj.phase0.states:J'] = phase.interp('J', [0, 1])
p['traj.phase0.controls:u'] = phase.interp('u', [0, 0])
return p
Simulating the Problem (without control)#
The following script creates an instance of the Dymos vanderpol problem and simulates it.
Since the problem was only simulated and not solved, the solution lines in the plots show only the initial guesses for \(x_0\), \(x_1\), and \(u\). The simulation lines shown in the plots are the system response with the control variable \(u\) held constant.
# Create the Dymos problem instance
p = vanderpol(transcription='gauss-lobatto', use_pyoptsparse=True, optimizer='IPOPT', num_segments=15)
# Run the problem (simulate only)
dm.run_problem(p, run_driver=False, simulate=True)
--- Constraint Report [traj] ---
--- phase0 ---
None
INFO: checking out_of_order
INFO:check_config:checking out_of_order
INFO: checking system
INFO:check_config:checking system
INFO: checking solvers
INFO:check_config:checking solvers
INFO: checking dup_inputs
INFO:check_config:checking dup_inputs
INFO: checking missing_recorders
INFO:check_config:checking missing_recorders
WARNING: The Problem has no recorder of any kind attached
WARNING:check_config:The Problem has no recorder of any kind attached
INFO: checking unserializable_options
INFO:check_config:checking unserializable_options
INFO: checking comp_has_no_outputs
INFO:check_config:checking comp_has_no_outputs
INFO: checking auto_ivc_warnings
INFO:check_config:checking auto_ivc_warnings
Model viewer data has already been recorded for Driver.
INFO: checking out_of_order
INFO:check_config:checking out_of_order
INFO: checking system
INFO:check_config:checking system
INFO: checking solvers
INFO:check_config:checking solvers
INFO: checking dup_inputs
INFO:check_config:checking dup_inputs
INFO: checking missing_recorders
INFO:check_config:checking missing_recorders
WARNING: The Problem has no recorder of any kind attached
WARNING:check_config:The Problem has no recorder of any kind attached
INFO: checking unserializable_options
INFO:check_config:checking unserializable_options
INFO: checking comp_has_no_outputs
INFO:check_config:checking comp_has_no_outputs
INFO: checking auto_ivc_warnings
INFO:check_config:checking auto_ivc_warnings
Simulating trajectory traj
Model viewer data has already been recorded for Driver.
Done simulating trajectory traj
The first two plots shows the variables \(x_0\) and \(x_1\) vs time. The third plots shows \(x_0\) vs. \(x_1\) (which will be mostly circular in the case of undamped oscillation). The final plot is the (fixed) control variable \(u\) vs time.
Show code cell source
from dymos.examples.plotting import plot_results
# Display the results
def vanderpol_plots():
sol = om.CaseReader('dymos_solution.db').get_case('final')
sim = om.CaseReader('dymos_simulation.db').get_case('final')
plot_results([('traj.phase0.timeseries.time',
'traj.phase0.timeseries.x1',
'time (s)',
'$x_1$ (V)'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.x0',
'time (s)',
'$x_0$ (V/s)'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.J',
'time (s)',
'J'),
('traj.phase0.timeseries.x0',
'traj.phase0.timeseries.x1',
'$x_0$ (V/s)',
'$x_1$ (V)'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.u',
'time (s)',
'control u'),
],
title='Van Der Pol Simulation',
p_sol=sol, p_sim=sim);
vanderpol_plots()
Solving the Optimal Control Problem#
The next example shows optimization followed by simulation.
With a successful optimization, the resulting plots show a good match between the simulated (with varying control) and optimized results. The state variables \(x_0\) and \(x_1\) as well as the control variable \(u\) are all driven to zero.
However one can notice that the distribution of segments/nodes in the phase results in a bit of error. If we take the simulation results as the truth, there is a bit of divergence in the state values of the solution and simulation towards the end.
import dymos as dm
from dymos.examples.plotting import plot_results
# Create the Dymos problem instance
p = vanderpol(transcription='gauss-lobatto', num_segments=15,
transcription_order=3, compressed=False, use_pyoptsparse=True, optimizer='IPOPT')
# Find optimal control solution to stop oscillation
dm.run_problem(p, simulate=True)
--- Constraint Report [traj] ---
--- phase0 ---
None
INFO: checking out_of_order
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmdao/recorders/sqlite_recorder.py:227: UserWarning:The existing case recorder file, dymos_solution.db, is being overwritten.
INFO:check_config:checking out_of_order
INFO: checking system
INFO:check_config:checking system
INFO: checking solvers
INFO:check_config:checking solvers
INFO: checking dup_inputs
INFO:check_config:checking dup_inputs
INFO: checking missing_recorders
INFO:check_config:checking missing_recorders
WARNING: The Problem has no recorder of any kind attached
WARNING:check_config:The Problem has no recorder of any kind attached
INFO: checking unserializable_options
INFO:check_config:checking unserializable_options
INFO: checking comp_has_no_outputs
INFO:check_config:checking comp_has_no_outputs
INFO: checking auto_ivc_warnings
INFO:check_config:checking auto_ivc_warnings
Model viewer data has already been recorded for Driver.
INFO: checking out_of_order
INFO:check_config:checking out_of_order
INFO: checking system
INFO:check_config:checking system
INFO: checking solvers
INFO:check_config:checking solvers
INFO: checking dup_inputs
INFO:check_config:checking dup_inputs
INFO: checking missing_recorders
INFO:check_config:checking missing_recorders
WARNING: The Problem has no recorder of any kind attached
WARNING:check_config:The Problem has no recorder of any kind attached
INFO: checking unserializable_options
INFO:check_config:checking unserializable_options
INFO: checking comp_has_no_outputs
INFO:check_config:checking comp_has_no_outputs
INFO: checking auto_ivc_warnings
INFO:check_config:checking auto_ivc_warnings
Full total jacobian was computed 3 times, taking 0.200369 seconds.
Total jacobian shape: (60, 130)
Jacobian shape: (60, 130) ( 5.15% nonzero)
FWD solves: 0 REV solves: 5
Total colors vs. total size: 5 vs 60 (91.7% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.200369 sec.
Time to compute coloring: 0.073859 sec.
Memory to compute coloring: 0.250000 MB.
List of user-set options:
Name Value used
file_print_level = 5 yes
hessian_approximation = limited-memory yes
linear_solver = mumps yes
nlp_scaling_method = user-scaling yes
output_file = IPOPT.out yes
print_level = 4 yes
print_user_options = yes yes
sb = yes yes
Total number of variables............................: 130
variables with only lower bounds: 0
variables with lower and upper bounds: 45
variables with only upper bounds: 0
Total number of equality constraints.................: 115
Total number of inequality cons
traints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
Number of Iterations....: 34
(scaled) (unscaled)
Objective...............: 2.8794727905908770e+00 2.8794727905908770e+00
Dual infeasibility......: 3.5078177740982541e-09 3.5078177740982541e-09
Constraint violation....: 1.0269562977782698e-14 1.0269562977782698e-14
Variable bound violation: 9.8074164345263171e-09 9.8074164345263171e-09
Complementarity.........: 1.0000000000000001e-11 1.0000000000000001e-11
Overall NLP error.......: 3.5078177740982541e-09 3.5078177740982541e-09
Number of objective function evaluations = 35
Number of objective gradient evaluations = 35
Number of equality constraint evaluations = 35
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 35
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 0
Total seconds in IPOPT = 0.653
EXIT: Optimal Solution Found.
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmdao/visualization/opt_report/opt_report.py:634: UserWarning: Attempting to set identical low and high ylims makes transformation singular; automatically expanding.
ax.set_ylim([ymin_plot, ymax_plot])
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmdao/recorders/sqlite_recorder.py:227: UserWarning:The existing case recorder file, dymos_simulation.db, is being overwritten.
Simulating trajectory traj
Model viewer data has already been recorded for Driver.
Done simulating trajectory traj
False
Show code cell source
vanderpol_plots()
Solving the Optimal Control Problem with Grid Refinement#
Repeating the optimization with grid refinement enabled requires changing only two lines in the code. For the sake of grid refinement demonstration, the initial number of segments is also reduced by a factor of 5.
Optimization with grid refinement gets results similar to the example without grid refinement, but runs faster and does not require supplying a good guess for the number segments.
# Create the Dymos problem instance
p = vanderpol(transcription='gauss-lobatto', num_segments=15,
transcription_order=3, compressed=True, optimizer='SLSQP')
# Enable grid refinement and find optimal control solution to stop oscillation
p.model.traj.phases.phase0.set_refine_options(refine=True, tol=1.0E-6)
dm.run_problem(p, simulate=True, refine_iteration_limit=5, refine_method='hp')
--- Constraint Report [traj] ---
--- phase0 ---
None
INFO: checking out_of_order
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmdao/recorders/sqlite_recorder.py:227: UserWarning:The existing case recorder file, dymos_solution.db, is being overwritten.
INFO:check_config:checking out_of_order
INFO: checking system
INFO:check_config:checking system
INFO: checking solvers
INFO:check_config:checking solvers
INFO: checking dup_inputs
INFO:check_config:checking dup_inputs
INFO: checking missing_recorders
INFO:check_config:checking missing_recorders
WARNING: The Problem has no recorder of any kind attached
WARNING:check_config:The Problem has no recorder of any kind attached
INFO: checking unserializable_options
INFO:check_config:checking unserializable_options
INFO: checking comp_has_no_outputs
INFO:check_config:checking comp_has_no_outputs
INFO: checking auto_ivc_warnings
INFO:check_config:checking auto_ivc_warnings
Model viewer data has already been recorded for Driver.
INFO: checking out_of_order
INFO:check_config:checking out_of_order
INFO: checking system
INFO:check_config:checking system
INFO: checking solvers
INFO:check_config:checking solvers
INFO: checking dup_inputs
INFO:check_config:checking dup_inputs
INFO: checking missing_recorders
INFO:check_config:checking missing_recorders
WARNING: The Problem has no recorder of any kind attached
WARNING:check_config:The Problem has no recorder of any kind attached
INFO: checking unserializable_options
INFO:check_config:checking unserializable_options
INFO: checking comp_has_no_outputs
INFO:check_config:checking comp_has_no_outputs
INFO: checking auto_ivc_warnings
INFO:check_config:checking auto_ivc_warnings
Full total jacobian was computed 3 times, taking 0.201899 seconds.
Total jacobian shape: (60, 74)
Jacobian shape: (60, 74) ( 8.74% nonzero)
FWD solves: 9 REV solves: 0
Total colors vs. total size: 9 vs 74 (87.8% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.201899 sec.
Time to compute coloring: 0.052703 sec.
Memory to compute coloring: 0.000000 MB.
Optimization terminated successfully (Exit mode 0)
Current function value: [2.87947342]
Iterations: 33
Function evaluations: 34
Gradient evaluations: 33
Optimization Complete
-----------------------------------
==================================================
Grid Refinement - Iteration 1
--------------------------------------------------
Phase: traj.phases.phase0
Refinement Options:
Allow Refinement = True
Tolerance = 1e-06
Min Order = 3
Max Order = 14
Original Grid:
Number of Segments = 15
Segment Ends = [-1.0, -0.8667, -0.7333, -0.6, -0.4667, -0.3333, -0.2, -0.0667, 0.0667, 0.2, 0.3333, 0.4667, 0.6, 0.7333, 0.8667, 1.0]
Segment Order = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]
Error = [0.009519, 0.003042, 0.0003774, 0.0006496, 0.0004868, 0.0002292, 0.0001781, 8.934e-05, 3.454e-05, 1.197e-05, 1.256e-05, 9.183e-06, 4.124e-06, 4.358e-06, 2.004e-05]
Phase: traj.phases.phase0
New Grid:
Number of Segments = 15
Segment Ends = [-1.0, -0.8667, -0.7333, -0.6, -0.4667, -0.3333, -0.2, -0.0667, 0.0667, 0.2, 0.3333, 0.4667, 0.6, 0.7333, 0.8667, 1.0]
Segment Order = [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5]
Refined: True
--- Constraint Report [traj] ---
--- phase0 ---
None
Model viewer data has already been recorded for Driver.
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmdao/core/system.py:4077: OpenMDAOWarning:Calling `list_inputs` before `final_setup` will only display the default values of variables and will not show the result of any `set_val` calls.
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/phase/phase.py:2774: OpenMDAOWarning:phase0.states:x0 specifies 'fix_final=True'. If the given restart file has a different final value this will overwrite the user-specified value
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/phase/phase.py:2774: OpenMDAOWarning:phase0.states:x1 specifies 'fix_final=True'. If the given restart file has a different final value this will overwrite the user-specified value
Full total jacobian was computed 3 times, taking 0.388824 seconds.
Total jacobian shape: (105, 149)
Jacobian shape: (105, 149) ( 6.78% nonzero)
FWD solves: 14 REV solves: 0
Total colors vs. total size: 14 vs 149 (90.6% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.388824 sec.
Time to compute coloring: 0.109039 sec.
Memory to compute coloring: 0.000000 MB.
Optimization terminated successfully (Exit mode 0)
Current function value: [2.87331094]
Iterations: 20
Function evaluations: 21
Gradient evaluations: 20
Optimization Complete
-----------------------------------
==================================================
Grid Refinement - Iteration 2
--------------------------------------------------
Phase: traj.phases.phase0
Refinement Options:
Allow Refinement = True
Tolerance = 1e-06
Min Order = 3
Max Order = 14
Original Grid:
Number of Segments = 15
Segment Ends = [-1.0, -0.8667, -0.7333, -0.6, -0.4667, -0.3333, -0.2, -0.0667, 0.0667, 0.2, 0.3333, 0.4667, 0.6, 0.7333, 0.8667, 1.0]
Segment Order = [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5]
Error = [7.669e-05, 0.0001911, 4.013e-05, 2.411e-06, 1.774e-06, 3.464e-06, 3.225e-07, 9.714e-07, 1.221e-07, 2.319e-07, 1.015e-08, 2.475e-08, 5.63e-08, 4.637e-08, 4.277e-07]
Phase: traj.phases.phase0
New Grid:
Number of Segments = 21
Segment Ends = [-1.0, -0.9333, -0.8667, -0.8, -0.7333, -0.6667, -0.6, -0.5333, -0.4667, -0.4, -0.3333, -0.2667, -0.2, -0.0667, 0.0667, 0.2, 0.3333, 0.4667, 0.6, 0.7333, 0.8667, 1.0]
Segment Order = [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5]
Refined: True
--- Constraint Report [traj] ---
--- phase0 ---
None
Model viewer data has already been recorded for Driver.
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmdao/core/system.py:4077: OpenMDAOWarning:Calling `list_inputs` before `final_setup` will only display the default values of variables and will not show the result of any `set_val` calls.
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/phase/phase.py:2774: OpenMDAOWarning:phase0.states:x0 specifies 'fix_final=True'. If the given restart file has a different final value this will overwrite the user-specified value
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/phase/phase.py:2774: OpenMDAOWarning:phase0.states:x1 specifies 'fix_final=True'. If the given restart file has a different final value this will overwrite the user-specified value
Full total jacobian was computed 3 times, taking 0.535802 seconds.
Total jacobian shape: (147, 209)
Jacobian shape: (147, 209) ( 4.88% nonzero)
FWD solves: 14 REV solves: 0
Total colors vs. total size: 14 vs 209 (93.3% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.535802 sec.
Time to compute coloring: 0.154523 sec.
Memory to compute coloring: 0.000000 MB.
Optimization terminated successfully (Exit mode 0)
Current function value: [2.87331097]
Iterations: 11
Function evaluations: 12
Gradient evaluations: 11
Optimization Complete
-----------------------------------
==================================================
Grid Refinement - Iteration 3
--------------------------------------------------
Phase: traj.phases.phase0
Refinement Options:
Allow Refinement = True
Tolerance = 1e-06
Min Order = 3
Max Order = 14
Original Grid:
Number of Segments = 21
Segment Ends = [-1.0, -0.9333, -0.8667, -0.8, -0.7333, -0.6667, -0.6, -0.5333, -0.4667, -0.4, -0.3333, -0.2667, -0.2, -0.0667, 0.0667, 0.2, 0.3333, 0.4667, 0.6, 0.7333, 0.8667, 1.0]
Segment Order = [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5]
Error = [1.619e-05, 3.825e-06, 2.313e-06, 1.966e-05, 2.497e-05, 5.579e-06, 1.382e-07, 4.725e-07, 5.89e-07, 1.296e-07, 1.736e-07, 2.793e-07, 1.346e-05, 1.342e-06, 1.148e-06, 5.102e-07, 2.147e-07, 1.072e-07, 2.596e-08, 2.332e-07, 3.348e-07]
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/grid_refinement/hp_adaptive/hp_adaptive.py:379: RuntimeWarning: divide by zero encountered in true_divide
q_smooth = (np.log(self.error[phase_path][smooth_need_refine_idxs] /
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/grid_refinement/hp_adaptive/hp_adaptive.py:411: RuntimeWarning: divide by zero encountered in true_divide
q_split = np.log((self.error[phase_path][split_seg_idxs] /
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmdao/core/system.py:4077: OpenMDAOWarning:Calling `list_inputs` before `final_setup` will only display the default values of variables and will not show the result of any `set_val` calls.
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/phase/phase.py:2774: OpenMDAOWarning:phase0.states:x0 specifies 'fix_final=True'. If the given restart file has a different final value this will overwrite the user-specified value
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/phase/phase.py:2774: OpenMDAOWarning:phase0.states:x1 specifies 'fix_final=True'. If the given restart file has a different final value this will overwrite the user-specified value
Phase: traj.phases.phase0
New Grid:
Number of Segments = 26
Segment Ends = [-1.0, -0.9667, -0.9333, -0.8667, -0.8, -0.7667, -0.7333, -0.7, -0.6667, -0.6333, -0.6, -0.5333, -0.4667, -0.4, -0.3333, -0.2667, -0.2, -0.1333, -0.0667, 0.0667, 0.2, 0.3333, 0.4667, 0.6, 0.7333, 0.8667, 1.0]
Segment Order = [5, 5, 13, 9, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 11, 9, 5, 5, 5, 5, 5, 5]
Refined: True
--- Constraint Report [traj] ---
--- phase0 ---
None
Model viewer data has already been recorded for Driver.
Full total jacobian was computed 3 times, taking 0.854925 seconds.
Total jacobian shape: (215, 314)
Jacobian shape: (215, 314) ( 4.16% nonzero)
FWD solves: 0 REV solves: 33
Total colors vs. total size: 33 vs 215 (84.7% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.854925 sec.
Time to compute coloring: 0.253186 sec.
Memory to compute coloring: 0.125000 MB.
Optimization terminated successfully (Exit mode 0)
Current function value: [2.87331113]
Iterations: 5
Function evaluations: 6
Gradient evaluations: 5
Optimization Complete
-----------------------------------
==================================================
Grid Refinement - Iteration 4
--------------------------------------------------
Phase: traj.phases.phase0
Refinement Options:
Allow Refinement = True
Tolerance = 1e-06
Min Order = 3
Max Order = 14
Original Grid:
Number of Segments = 26
Segment Ends = [-1.0, -0.9667, -0.9333, -0.8667, -0.8, -0.7667, -0.7333, -0.7, -0.6667, -0.6333, -0.6, -0.5333, -0.4667, -0.4, -0.3333, -0.2667, -0.2, -0.1333, -0.0667, 0.0667, 0.2, 0.3333, 0.4667, 0.6, 0.7333, 0.8667, 1.0]
Segment Order = [5, 5, 13, 9, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 11, 9, 5, 5, 5, 5, 5, 5]
Error = [9.881e-07, 5.651e-07, 2.643e-07, 2.474e-06, 1.029e-06, 3.243e-08, 1.71e-06, 1.8e-06, 2.133e-07, 2.365e-07, 3.385e-07, 7.385e-07, 1.061e-07, 1.627e-07, 1.072e-06, 1.285e-06, 1.692e-07, 1.029e-06, 7.225e-07, 2.977e-07, 5.927e-07, 2.915e-07, 8.94e-08, 4.391e-08, 1.129e-07, 1.193e-07]
Phase: traj.phases.phase0
New Grid:
Number of Segments = 27
Segment Ends = [-1.0, -0.9667, -0.9333, -0.8667, -0.8333, -0.8, -0.7667, -0.7333, -0.7, -0.6667, -0.6333, -0.6, -0.5333, -0.4667, -0.4, -0.3333, -0.2667, -0.2, -0.1333, -0.0667, 0.0667, 0.2, 0.3333, 0.4667, 0.6, 0.7333, 0.8667, 1.0]
Segment Order = [5, 5, 13, 9, 9, 7, 5, 9, 9, 5, 5, 5, 5, 5, 5, 7, 11, 5, 7, 11, 9, 5, 5, 5, 5, 5, 5]
Refined: True
--- Constraint Report [traj] ---
--- phase0 ---
None
Model viewer data has already been recorded for Driver.
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmdao/core/system.py:4077: OpenMDAOWarning:Calling `list_inputs` before `final_setup` will only display the default values of variables and will not show the result of any `set_val` calls.
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/phase/phase.py:2774: OpenMDAOWarning:phase0.states:x0 specifies 'fix_final=True'. If the given restart file has a different final value this will overwrite the user-specified value
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/phase/phase.py:2774: OpenMDAOWarning:phase0.states:x1 specifies 'fix_final=True'. If the given restart file has a different final value this will overwrite the user-specified value
Full total jacobian was computed 3 times, taking 1.117113 seconds.
Total jacobian shape: (258, 384)
Jacobian shape: (258, 384) ( 3.83% nonzero)
FWD solves: 0 REV solves: 33
Total colors vs. total size: 33 vs 258 (87.2% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 1.117113 sec.
Time to compute coloring: 0.306057 sec.
Memory to compute coloring: 0.000000 MB.
Optimization terminated successfully (Exit mode 0)
Current function value: [2.87331162]
Iterations: 2
Function evaluations: 3
Gradient evaluations: 2
Optimization Complete
-----------------------------------
==================================================
Grid Refinement - Iteration 5
--------------------------------------------------
Phase: traj.phases.phase0
Refinement Options:
Allow Refinement = True
Tolerance = 1e-06
Min Order = 3
Max Order = 14
Original Grid:
Number of Segments = 27
Segment Ends = [-1.0, -0.9667, -0.9333, -0.8667, -0.8333, -0.8, -0.7667, -0.7333, -0.7, -0.6667, -0.6333, -0.6, -0.5333, -0.4667, -0.4, -0.3333, -0.2667, -0.2, -0.1333, -0.0667, 0.0667, 0.2, 0.3333, 0.4667, 0.6, 0.7333, 0.8667, 1.0]
Segment Order = [5, 5, 13, 9, 9, 7, 5, 9, 9, 5, 5, 5, 5, 5, 5, 7, 11, 5, 7, 11, 9, 5, 5, 5, 5, 5, 5]
Error = [1.02e-06, 4.302e-07, 8.674e-08, 6.22e-07, 7.625e-07, 2.013e-07, 5.558e-08, 3.517e-07, 1.383e-07, 1.936e-07, 2.45e-07, 2.435e-07, 4.368e-07, 2.104e-07, 5.198e-07, 3.064e-07, 5.049e-07, 1.66e-07, 2.793e-07, 4.327e-07, 1.777e-07, 1.886e-07, 1.793e-07, 4.019e-08, 2.256e-08, 4.069e-08, 3.824e-08]
Phase: traj.phases.phase0
New Grid:
Number of Segments = 27
Segment Ends = [-1.0, -0.9667, -0.9333, -0.8667, -0.8333, -0.8, -0.7667, -0.7333, -0.7, -0.6667, -0.6333, -0.6, -0.5333, -0.4667, -0.4, -0.3333, -0.2667, -0.2, -0.1333, -0.0667, 0.0667, 0.2, 0.3333, 0.4667, 0.6, 0.7333, 0.8667, 1.0]
Segment Order = [7, 5, 13, 9, 9, 7, 5, 9, 9, 5, 5, 5, 5, 5, 5, 7, 11, 5, 7, 11, 9, 5, 5, 5, 5, 5, 5]
Refined: True
--- Constraint Report [traj] ---
--- phase0 ---
None
Model viewer data has already been recorded for Driver.
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmdao/core/system.py:4077: OpenMDAOWarning:Calling `list_inputs` before `final_setup` will only display the default values of variables and will not show the result of any `set_val` calls.
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/phase/phase.py:2774: OpenMDAOWarning:phase0.states:x0 specifies 'fix_final=True'. If the given restart file has a different final value this will overwrite the user-specified value
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/phase/phase.py:2774: OpenMDAOWarning:phase0.states:x1 specifies 'fix_final=True'. If the given restart file has a different final value this will overwrite the user-specified value
Full total jacobian was computed 3 times, taking 1.085524 seconds.
Total jacobian shape: (261, 389)
Jacobian shape: (261, 389) ( 3.79% nonzero)
FWD solves: 0 REV solves: 33
Total colors vs. total size: 33 vs 261 (87.4% improvement)
Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 1.085524 sec.
Time to compute coloring: 0.311173 sec.
Memory to compute coloring: 0.000000 MB.
Optimization terminated successfully (Exit mode 0)
Current function value: [2.87331173]
Iterations: 2
Function evaluations: 3
Gradient evaluations: 2
Optimization Complete
-----------------------------------
Successfully completed grid refinement.
==================================================
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmdao/recorders/sqlite_recorder.py:227: UserWarning:The existing case recorder file, dymos_simulation.db, is being overwritten.
Simulating trajectory traj
Model viewer data has already been recorded for Driver.
Done simulating trajectory traj
False
Show code cell source
vanderpol_plots()
References#
[Wikipediacontributors20]
Wikipedia contributors. Van der pol oscillator — Wikipedia, the free encyclopedia. 2020. [Online; accessed 12-June-2020]. URL: https://en.wikipedia.org/w/index.php?title=Van_der_Pol_oscillator&oldid=953857769.