Skip to content

Instantly share code, notes, and snippets.

@moorepants
Created November 2, 2025 07:08
Show Gist options
  • Select an option

  • Save moorepants/8fdb6507d16d3c4027c6185b4c006cf5 to your computer and use it in GitHub Desktop.

Select an option

Save moorepants/8fdb6507d16d3c4027c6185b4c006cf5 to your computer and use it in GitHub Desktop.
import numpy as np
import sympy as sm
import matplotlib.pyplot as plt
from opty import Problem, create_objective_function
t = sm.symbols('t', real=True)
x1, x2, v = sm.Function('x1')(t), sm.Function('x2')(t), sm.Function('v')(t)
eqs = sm.Matrix([x1.diff() - x1 + x1*x2 + 0.4*x1*v,
x2.diff() + x2 - x1*x2 + 0.2*x2*v])
bounds = {v: (0.0, 1.0)}
instance = (x1.func(0.0) - 0.5, x2.func(0.0) - 0.7)
duration = 12.0
num_nodes = 48
h = duration/(num_nodes - 1)
obj_func = sm.Integral((x1 - 1)**2 + (x2 - 1)**2, t)
obj, obj_grad = create_objective_function(
obj_func,
(x1, x2),
(v,),
tuple(),
num_nodes,
h,
time_symbol=t,
)
p = Problem(
obj,
obj_grad,
eqs,
(x1, x2),
num_nodes,
h,
instance_constraints=instance,
bounds=bounds,
time_symbol=t,
)
sol, info = p.solve(np.ones(p.num_free))
p.plot_trajectories(sol)
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment