Skip to content

Instantly share code, notes, and snippets.

@moorepants

moorepants/so_68266571.py

Created November 2, 2025 07:08
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save moorepants/8fdb6507d16d3c4027c6185b4c006cf5 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save moorepants/8fdb6507d16d3c4027c6185b4c006cf5 to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. moorepants created this gist Nov 2, 2025.
    45 changes: 45 additions & 0 deletions so_68266571.py
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,45 @@
    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()