Two-Phase Schwartz Problem#

Mathematical Formulation#

Problem Statement#

Find the optimal control \(u(t)\) that minimizes the final state penalty:

\[J = 5(x_0^2(t_f) + x_1^2(t_f))\]

Subject to the nonlinear dynamics for both phases:

\[\frac{dx_0}{dt} = x_1\]
\[\frac{dx_1}{dt} = u - 0.1(1 + 2x_0^2)x_1\]

Phase-Specific Constraints#

Phase 1 (\(t \in [0,1]\)):

  • Elliptical exclusion constraint: \(1 - 9(x_0 - 1)^2 - \left(\frac{x_1 - 0.4}{0.3}\right)^2 \leq 0\)

  • Control bounds: \(-1 \leq u \leq 1\)

Phase 2 (\(t \in [1,2.9]\)):

  • No path constraints

  • Control bounds: \(-50 \leq u \leq 50\)

Boundary Conditions#

Phase 1:

  • Initial conditions: \(x_0(0) = 1\), \(x_1(0) = 1\)

  • Final conditions: Continuous with Phase 2

Phase 2:

  • Initial conditions: Continuous from Phase 1

  • Final conditions: \(x_0(2.9) = \text{free}\), \(x_1(2.9) = \text{free}\)

  • State bounds: \(-20 \leq x_0 \leq 10\), \(-10 \leq x_1 \leq 10\) (Phase 1), \(-10 \leq x_1 \leq 10\) (Phase 2)

Physical Parameters#

  • Phase 1 duration: \(t_1 = 1\) s

  • Phase 2 duration: \(t_2 = 1.9\) s

  • System damping parameter: \(\alpha = 0.1\)

  • Final state penalty weight: \(W = 5\)

State Variables#

  • \(x_0(t)\): System state coordinate 1

  • \(x_1(t)\): System state coordinate 2

Control Variable#

  • \(u(t)\): Control input

Notes#

This problem features different feasible regions and control authority in each phase, demonstrating the flexibility of multiphase optimal control formulations. The elliptical constraint in Phase 1 creates a forbidden region that the trajectory must avoid.

Running This Example#

cd examples/two_phase_schwartz
python two_phase_schwartz.py

Code Implementation#

examples/two_phase_schwartz/two_phase_schwartz.py#
  1import numpy as np
  2
  3import maptor as mtor
  4
  5
  6# Problem setup
  7problem = mtor.Problem("Two-Phase Schwartz Problem")
  8
  9# Phase 1
 10phase1 = problem.set_phase(1)
 11phase1.time(initial=0.0, final=1.0)
 12x0_1 = phase1.state("x0", initial=1.0, boundary=(-20.0, 10.0))
 13x1_1 = phase1.state("x1", initial=1.0, boundary=(-0.8, 10.0))
 14u1 = phase1.control("u", boundary=(-1.0, 1.0))
 15
 16phase1.dynamics(
 17    {
 18        x0_1: x1_1,
 19        x1_1: u1 - 0.1 * (1 + 2 * x0_1**2) * x1_1,
 20    }
 21)
 22
 23# Path constraint: feasible region outside ellipse
 24elliptical_constraint = 1 - 9 * (x0_1 - 1) ** 2 - ((x1_1 - 0.4) / 0.3) ** 2
 25phase1.path_constraints(elliptical_constraint <= 0)
 26phase1.mesh([20], [-1.0, 1.0])
 27
 28# Phase 2
 29phase2 = problem.set_phase(2)
 30phase2.time(initial=1.0, final=2.9)
 31x0_2 = phase2.state("x0", initial=x0_1.final, boundary=(-20.0, 10.0))
 32x1_2 = phase2.state("x1", initial=x1_1.final, boundary=(-10.0, 10.0))
 33u2 = phase2.control("u", boundary=(-50.0, 50.0))
 34
 35phase2.dynamics(
 36    {
 37        x0_2: x1_2,
 38        x1_2: u2 - 0.1 * (1 + 2 * x0_2**2) * x1_2,
 39    }
 40)
 41phase2.mesh([20], [-1.0, 1.0])
 42
 43# Objective
 44objective_expr = 5 * (x0_2.final**2 + x1_2.final**2)
 45problem.minimize(objective_expr)
 46
 47# Guess
 48states_p1 = []
 49controls_p1 = []
 50states_p2 = []
 51controls_p2 = []
 52
 53for N in [20]:
 54    tau_states = np.linspace(-1, 1, N + 1)
 55    t_norm_states = (tau_states + 1) / 2
 56    x0_vals = 1.0 + 0.2 * t_norm_states
 57    x1_vals = 1.0 - 0.3 * t_norm_states
 58    states_p1.append(np.array([x0_vals, x1_vals]))
 59
 60    t_norm_controls = np.linspace(0, 1, N)
 61    u_vals = 0.3 * np.sin(np.pi * t_norm_controls)
 62    controls_p1.append(np.array([u_vals]))
 63
 64for N in [20]:
 65    tau_states = np.linspace(-1, 1, N + 1)
 66    t_norm_states = (tau_states + 1) / 2
 67    x0_end_p1 = 1.2
 68    x1_end_p1 = 0.7
 69    x0_vals = x0_end_p1 * (1 - 0.8 * t_norm_states)
 70    x1_vals = x1_end_p1 * (1 - 0.9 * t_norm_states)
 71    states_p2.append(np.array([x0_vals, x1_vals]))
 72
 73    t_norm_controls = np.linspace(0, 1, N)
 74    u_vals = -1.0 + 0.5 * t_norm_controls
 75    controls_p2.append(np.array([u_vals]))
 76
 77phase1.guess(
 78    states=states_p1,
 79    controls=controls_p1,
 80    initial_time=0.0,
 81    terminal_time=1.0,
 82)
 83
 84phase2.guess(
 85    states=states_p2,
 86    controls=controls_p2,
 87    initial_time=1.0,
 88    terminal_time=2.9,
 89)
 90
 91# Solve
 92solution = mtor.solve_adaptive(
 93    problem,
 94    error_tolerance=1e-5,
 95    max_iterations=30,
 96    min_polynomial_degree=3,
 97    max_polynomial_degree=20,
 98    nlp_options={
 99        "ipopt.print_level": 0,
100        "ipopt.max_iter": 500,
101        "ipopt.tol": 1e-4,
102        "ipopt.constr_viol_tol": 1e-4,
103        "ipopt.acceptable_tol": 1e-3,
104        "ipopt.mu_strategy": "adaptive",
105        "ipopt.linear_solver": "mumps",
106    },
107)
108
109
110# Results
111if solution.status["success"]:
112    print(f"Objective: {solution.status['objective']:.8f}")
113    print(f"Total mission time: {solution.status['total_mission_time']:.6f}")
114
115    # Variable access
116    x0_final = solution[(2, "x0")][-1]
117    x1_final = solution[(2, "x1")][-1]
118    print(f"Final state: x0={x0_final:.6f}, x1={x1_final:.6f}")
119
120    # Plotting
121    solution.plot()
122
123else:
124    print(f"Failed: {solution.status['message']}")