# The Van der Pol Oscillator

## Contents

# 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:

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:

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):

## 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:

The final conditions are:

## 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
```

```
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.

```
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.states:x1',
'time (s)',
'$x_1$ (V)'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.states:x0',
'time (s)',
'$x_0$ (V/s)'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.states:J',
'time (s)',
'J'),
('traj.phase0.timeseries.states:x0',
'traj.phase0.timeseries.states:x1',
'$x_0$ (V/s)',
'$x_1$ (V)'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.controls:u',
'time (s)',
'control u'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.state_rates:x1',
'time (s)',
'$\dot{x}_1$ (V)'),
('traj.phase0.timeseries.time',
'traj.phase0.timeseries.state_rates:x0',
'time (s)',
'$\dot{x}_0$ (V/s)'),
],
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.222382 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.222382 sec.
Time to compute coloring: 0.077306 sec.
Memory to compute coloring: 0.000000 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 constraints...............: 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.8794727905908775e+00 2.8794727905908775e+00
Dual infeasibility......: 3.5078177584857428e-09 3.5078177584857428e-09
Constraint violation....: 3.6082248300317588e-15 3.6082248300317588e-15
Variable bound violation: 9.8074164345263171e-09 9.8074164345263171e-09
Complementarity.........: 1.0000000000000001e-11 1.0000000000000001e-11
Overall NLP error.......: 3.5078177584857428e-09 3.5078177584857428e-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.726
EXIT: Optimal Solution Found.
```

```
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmdao/visualization/opt_report/opt_report.py:666: UserWarning: Attempting to set identical low and high ylims makes transformation singular; automatically expanding.
ax.set_ylim([ymin_plot, ymax_plot])
```

```
Simulating trajectory traj
```

```
/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.
```

```
Done simulating trajectory traj
```

```
False
```

```
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
```

```
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.204216 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.204216 sec.
Time to compute coloring: 0.055257 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.
```

```
Model viewer data has already been recorded for Driver.
```

```
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/load_case.py:157: 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/load_case.py:157: 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.409311 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.409311 sec.
Time to compute coloring: 0.114417 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.
Model viewer data has already been recorded for Driver.
```

```
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/load_case.py:157: 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/load_case.py:157: 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.582097 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.582097 sec.
Time to compute coloring: 0.158434 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:380: 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:412: RuntimeWarning: divide by zero encountered in true_divide
q_split = np.log((self.error[phase_path][split_seg_idxs] /
```

```
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.
```

```
Model viewer data has already been recorded for Driver.
```

```
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/load_case.py:157: 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/load_case.py:157: 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.918048 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.918048 sec.
Time to compute coloring: 0.267339 sec.
Memory to compute coloring: 0.000000 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.
Model viewer data has already been recorded for Driver.
```

```
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/load_case.py:157: 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/load_case.py:157: 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.196457 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.196457 sec.
Time to compute coloring: 0.342971 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.
Model viewer data has already been recorded for Driver.
```

```
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/load_case.py:157: 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/load_case.py:157: 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.197115 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.197115 sec.
Time to compute coloring: 0.317540 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.
==================================================
```

```
Simulating trajectory traj
```

```
Done simulating trajectory traj
```

```
False
```

```
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.