LQR Problem#

Mathematical Formulation#

Problem Statement#

Find the optimal control \(u(t)\) that minimizes the quadratic cost functional:

\[J = \int_0^1 \left( 0.625 x^2 + 0.5 x u + 0.5 u^2 \right) dt\]

Subject to the linear dynamic constraint:

\[\frac{dx}{dt} = 0.5 x + u\]

Boundary Conditions#

  • Initial condition: \(x(0) = 1.0\)

  • Final condition: \(x(1) = \text{free}\)

  • Control bounds: \(u \in \mathbb{R}\) (unconstrained)

Physical Parameters#

  • Time horizon: \(t_f = 1.0\) s

  • System parameter: \(a = 0.5\)

State Variables#

  • \(x(t)\): System state

Control Variable#

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

Running This Example#

cd examples/lqr
python lqr.py

Code Implementation#

examples/lqr/lqr.py#
  1import numpy as np
  2
  3import maptor as mtor
  4
  5
  6# Problem setup
  7problem = mtor.Problem("LQR Problem")
  8phase = problem.set_phase(1)
  9
 10# Variables
 11t = phase.time(initial=0, final=1)
 12x = phase.state("x", initial=1.0)
 13u = phase.control("u")
 14
 15# Dynamics
 16phase.dynamics({x: 0.5 * x + u})
 17
 18# Objective
 19integrand = 0.625 * x**2 + 0.5 * x * u + 0.5 * u**2
 20integral_var = phase.add_integral(integrand)
 21problem.minimize(integral_var)
 22
 23# Mesh and guess
 24phase.mesh([4, 4], [-1.0, 0.0, 1.0])
 25
 26# Initial guess - explicit arrays for each mesh interval
 27states_guess = []
 28controls_guess = []
 29
 30# First interval: 5 state points, 4 control points
 31states_guess.append([[1.0, 0.8, 0.6, 0.4, 0.2]])
 32controls_guess.append([[-0.5, -0.4, -0.3, -0.2]])
 33
 34# Second interval: 5 state points, 4 control points
 35states_guess.append([[0.2, 0.15, 0.1, 0.05, 0.0]])
 36controls_guess.append([[-0.1, -0.05, 0.0, 0.05]])
 37
 38phase.guess(
 39    states=states_guess,
 40    controls=controls_guess,
 41    initial_time=0.0,
 42    terminal_time=1.0,
 43    integrals=0.5,
 44)
 45
 46# Solve with adaptive mesh
 47solution = mtor.solve_adaptive(
 48    problem,
 49    error_tolerance=1e-6,
 50    max_iterations=20,
 51    min_polynomial_degree=3,
 52    max_polynomial_degree=10,
 53    nlp_options={"ipopt.print_level": 0, "ipopt.max_iter": 200},
 54)
 55
 56# Results
 57if solution.status["success"]:
 58    print(f"Adaptive objective: {solution.status['objective']:.12f}")
 59    print("Literature reference: 0.380797077977481140")
 60    print(f"Error: {abs(solution.status['objective'] - 0.380797077977481140):.2e}")
 61
 62    # Final state value
 63    x_final = solution["x"][-1]
 64    print(f"Final state: x(1) = {x_final:.6f}")
 65
 66    solution.plot()
 67
 68    # Compare with fixed mesh
 69    phase.mesh([8, 8], [-1.0, 0.0, 1.0])
 70
 71    states_guess_fixed = []
 72    controls_guess_fixed = []
 73
 74    # First interval: 9 state points, 8 control points
 75    x_vals_1 = np.linspace(1.0, 0.2, 9)
 76    states_guess_fixed.append(x_vals_1.reshape(1, -1))
 77    u_vals_1 = np.linspace(-0.5, -0.1, 8)
 78    controls_guess_fixed.append(u_vals_1.reshape(1, -1))
 79
 80    # Second interval: 9 state points, 8 control points
 81    x_vals_2 = np.linspace(0.2, 0.0, 9)
 82    states_guess_fixed.append(x_vals_2.reshape(1, -1))
 83    u_vals_2 = np.linspace(-0.1, 0.05, 8)
 84    controls_guess_fixed.append(u_vals_2.reshape(1, -1))
 85
 86    phase.guess(
 87        states=states_guess_fixed,
 88        controls=controls_guess_fixed,
 89        initial_time=0.0,
 90        terminal_time=1.0,
 91        integrals=0.5,
 92    )
 93
 94    fixed_solution = mtor.solve_fixed_mesh(
 95        problem, nlp_options={"ipopt.print_level": 0, "ipopt.max_iter": 200}
 96    )
 97
 98    if fixed_solution.status["success"]:
 99        print(f"Fixed mesh objective: {fixed_solution.status['objective']:.12f}")
100        print(
101            f"Difference: {abs(solution.status['objective'] - fixed_solution.status['objective']):.2e}"
102        )
103    else:
104        print(f"Fixed mesh failed: {fixed_solution.status['message']}")
105
106else:
107    print(f"Failed: {solution.status['message']}")