Linear-Quadratic Regulator

Problem

The Linear-Quadratic Regulator (LQR) is one of the most classic optimal control problems. This section demonstrates the usage of pockit by solving the most basic LQR problem.

The system contains only one phase. There is \(n_x = 1\) state variable \(x\) and \(n_u = 1\) control variable \(u\). The dynamic equation is

$$ \dot x = ax + bu $$

where \(a = -1, b = 1\) are given parameters. There is a single integral

$$ \mathbb{I} = \int_{t_0}^{t_f} \left( q x^2 + r u^2 \right) \mathrm dt $$

with parameters \(q = 1, r = 0.1\). The time interval of the phase is fixed as \(t_0 = 0, t_f = 1\), with boundary conditions

$$ x_0 = 1, x_f = \text{free} $$

At the system level, there is an objective function

$$ F = \mathbb I + \frac{s x_f^2}{2} $$

with parameter \(s = 1\).

Code

from pockit.lobatto import System, constant_guess
from pockit.optimizer import ipopt
import matplotlib.pyplot as plt

# LQR problem:
# min ∫_0^1 (q * x^2 + r * u^2) dt + s * x_f^2 / 2
# s.t. x' = a * x + b * u, x(0) = 1

# Set parameters
a, b, s, q, r = -1, 1, 1, 1, 0.1

# Set up the system
system = System(["x_f"])  # the system has one free parameter x_f
x_f, = system.s  # extract the Sympy symbol x_f
phase = system.new_phase(["x"], ["u"])  # the phase has one state x and one control u
x, = phase.x  # extract the Sympy symbol x
u, = phase.u  # extract the Sympy symbol u
phase.set_dynamics([a * x + b * u])  # x' = a * x + b * u
phase.set_integral([q * x ** 2 + r * u ** 2])  # I = ∫ q * x^2 + r * u^2 dt
phase.set_boundary_condition([1], [x_f], 0, 1)  # x(0) = 1, x(t_f) = x_f, t_0 = 0, t_f = 1
phase.set_discretization(10, 10)  # 10 subintervals with 10 collocation points in each subinterval
system.set_phase([phase])  # bind the phase to the system
system.set_objective(phase.I[0] + s * x_f ** 2 / 2)  # define the objective function

# Solve the problem
guess_p = constant_guess(phase, 0)  # initial guess for the phase
guess_s = [0.]  # initial guess for the free parameter x_f
# var_p, var_s: the optimal solution for the phase and the free parameter
[var_p, var_s], info = ipopt.solve(system, [guess_p, guess_s])

# Print the results
print("status:", info["status_msg"].decode())
print("objective:", info["obj_val"])  # 0.2319139744522318

# Plot the results
plt.plot(var_p.t_x, var_p.x[0], label="x")  # only one state variable and one control variable, 
plt.plot(var_p.t_u, var_p.u[0], label="u")  # so the indices are 0 for both
plt.legend()
plt.minorticks_on()
plt.grid(linestyle='--')
plt.show()

Results