diff --git a/examples-gallery/plot_betts_10_144_145.py b/examples-gallery/plot_betts_10_144_145.py new file mode 100644 index 00000000..6c626b63 --- /dev/null +++ b/examples-gallery/plot_betts_10_144_145.py @@ -0,0 +1,220 @@ +""" +Van der Pol Oscillator +====================== + +https://en.wikipedia.org/wiki/Van_der_Pol_oscillator + +These are example 10.144 / 145 from Betts' book "Practical Methods for Optimal Control +Using NonlinearProgramming", 3rd edition, Chapter 10: Test Problems. +It is described in more detail in section 4.14. example 4.11 of the book. + +""" +import numpy as np +import sympy as sm +import sympy.physics.mechanics as me +from opty.direct_collocation import Problem +from opty.utils import create_objective_function +import matplotlib.pyplot as plt + +# %% +# First version of the problem (10.144). +# -------------------------------------- +# **States** +# :math:`y_1, y_2` : state variables +# +# **Controls** +# :math:`u` : control variables +# +# Equations of Motion. +# -------------------- +t = me.dynamicsymbols._t + +y1, y2 = me.dynamicsymbols('y1, y2') +u = me.dynamicsymbols('u') + +eom = sm.Matrix([ + -y1.diff() + y2, + -y2.diff(t) + (1 - y1**2)*y2 - y1 + u, +]) +sm.pprint(eom) + +# %% +# Define and Solve the Optimization Problem. +# ------------------------------------------ +num_nodes = 2001 +iterations = 1000 + +t0, tf = 0.0, 5.0 +interval_value = (tf - t0)/(num_nodes - 1) + +state_symbols = (y1, y2) +unkonwn_input_trajectories = (u, ) + +objective = sm.integrate(u**2 + y1**2 + y2**2, t) +obj, obj_grad = create_objective_function(objective, + state_symbols, + unkonwn_input_trajectories, + tuple(), + num_nodes, + interval_value +) + +instance_constraints = ( + y1.func(t0) - 1.0, + y2.func(t0), +) + +bounds = { + y2: (-0.4, np.inf), +} + +prob = Problem( + obj, + obj_grad, + eom, + state_symbols, + num_nodes, + interval_value, + instance_constraints= instance_constraints, + bounds=bounds, +) + +# %% +# Solve the optimization problem. +# Give some rough estimates for the trajectories. + +initial_guess = np.ones(prob.num_free) + +# Find the optimal solution. +for i in range(1): + solution, info = prob.solve(initial_guess) + initial_guess = solution + print(info['status_msg']) + print(f'Objectve is: {info['obj_val']:.8f}, ' + + f'as per the book it is {2.95369916}, so the deviation is: ' + + f'{(info['obj_val'] -2.95369916) /2.95369916*100 :.5e} %') +solution1 = solution +# %% +# Plot the optimal state and input trajectories. +prob.plot_trajectories(solution) + +# %% +# Plot the constraint violations. +prob.plot_constraint_violations(solution) + +# %% +# Plot the objective function as a function of optimizer iteration. +prob.plot_objective_value() + +# %% +# Second version of the problem (10.145). +# --------------------------------------- +# This is example 10.145 from Betts' book "Practical Methods for Optimal Control +# Using NonlinearProgramming", 3rd edition, Chapter 10: Test Problems. +# It is same as problem 10.144 but a bit reformulated. +# It has two control variables and one additional algebraic equation of motion. +# As opty needs as many state variables as equations of motion, I simply call +# the additional control variable v a state variable. +# It is described in more detail in section 4.14, example 4.11 of the book. +# +# This formulation seems to be more accurate compared to the above, when +# considering the violations of the constraints. +# +# **States** +# :math:`y_1, y_2, v` : state variables +# +# **Controls** +# :math:`u` : control variables +# +# Equations of Motion. +# -------------------- +y1, y2, v = me.dynamicsymbols('y1, y2, v') +u = me.dynamicsymbols('u') + +eom = sm.Matrix([ + -y1.diff() + y2, + -y2.diff(t) + v - y1 + u, + v - (1-y1**2)*y2, +]) +sm.pprint(eom) + +# %% +# Define and Solve the Optimization Problem. +# ------------------------------------------ +state_symbols = (y1, y2, v) +unkonwn_input_trajectories = (u,) + +objective = sm.integrate(u**2 + y1**2 + y2**2, t) +obj, obj_grad = create_objective_function(objective, + state_symbols, + unkonwn_input_trajectories, + tuple(), + num_nodes, + interval_value +) + +instance_constraints = ( + y1.func(t0) - 1.0, + y2.func(t0), +) + +bounds = { + y2: (-0.4, np.inf), +} + +prob = Problem( + obj, + obj_grad, + eom, + state_symbols, + num_nodes, + interval_value, + instance_constraints= instance_constraints, + bounds=bounds, +) + +# %% +# Solve the optimization problem. +# Give some rough estimates for the trajectories. + +initial_guess = np.ones(prob.num_free) + +# Find the optimal solution. +for i in range(1): + solution, info = prob.solve(initial_guess) + initial_guess = solution + print(info['status_msg']) + print(f'Objectve is: {info['obj_val']:.8f}, ' + + f'as per the book it is {2.95369916}, so the deviation is: ' + + f'{(info['obj_val'] -2.95369916) /2.95369916*100 :.5e} %') +# %% +# Plot the optimal state and input trajectories. +prob.plot_trajectories(solution) + +# %% +# Plot the constraint violations. +prob.plot_constraint_violations(solution) + +# %% +# Plot the objective function as a function of optimizer iteration. +prob.plot_objective_value() + +# %% +# Plot the Difference between the two Solutions. +# ---------------------------------------------- +diffy1 = solution1[: num_nodes] - solution[: num_nodes] +diffy2 = solution1[num_nodes: 2*num_nodes] - solution[num_nodes: 2*num_nodes] +diffu = solution1[2*num_nodes:] - solution[3*num_nodes:] +times = np.linspace(t0, tf, num_nodes) + +fig, ax = plt.subplots(3, 1, figsize=(7, 4), sharex=True, constrained_layout=True) +ax[0].plot(times, diffy1, label='Delta y1') +ax[0].legend() +ax[1].plot(times, diffy2, label='Delta y2') +ax[1].legend() +ax[2].plot(times, diffu, label='Delta u') +ax[2].legend() +ax[2].set_xlabel('Time') +ax[0].set_title('Differences between the two solutions') +avoid_printing = True +# sphinx_gallery_thumbnail_number = 2