The objective of this library is to model and solve optimal control problems (OCPs) of the form
\begin{align} \min_{x,u}~& \int_0^T \ell(x, u)\, dt + \ell_\mathrm{f}(x(T)) \\ \subjectto & \dot x (t) = f(x(t), u(t)) \\ & g(x(t), u(t)) = 0 \\ & h(x(t), u(t)) \leq 0 \end{align}
1
2
3
4
5
6
7"""Formulating and solving a linear quadratic regulator with Aligator."""
8
9import matplotlib.pyplot as plt
10import numpy as np
11import tap
12
13import aligator
14from aligator import constraints, dynamics, manifolds
16
17
18class Args(tap.Tap):
19 term_cstr: bool = False
20 bounds: bool = False
21
22
23args = Args().parse_args()
24np.random.seed(42)
25
26TAG = "LQR"
27if args.bounds:
28 TAG += "_bounded"
29if args.term_cstr:
30 TAG += "_cstr"
31
32nx = 3
33nu = 3
34space = manifolds.VectorSpace(nx)
35x0 = space.neutral() + (0.2, 0.3, -0.1)
36
37
38A = np.eye(nx)
39A[0, 1] = -0.2
40A[1, 0] = 0.2
41B = np.eye(nx)[:, :nu]
42B[2, :] = 0.4
43c = np.zeros(nx)
44c[:] = (0.0, 0.0, 0.1)
45
46
47Q = 1e-2 * np.eye(nx)
48R = 1e-2 * np.eye(nu)
49N = 1e-5 * np.eye(nx, nu)
50
51Qf = np.eye(nx)
52if args.term_cstr:
53 Qf[:, :] = 0.0
54
55
56
57
58rcost0 = aligator.QuadraticCost(Q, R, N)
59term_cost = aligator.QuadraticCost(Qf, R)
60dynmodel = dynamics.LinearDiscreteDynamics(A, B, c)
61stage = aligator.StageModel(rcost0, dynmodel)
62if args.bounds:
63 u_min = -0.18 * np.ones(nu)
64 u_max = +0.18 * np.ones(nu)
65 ctrl_fn = aligator.ControlErrorResidual(nx, np.zeros(nu))
66 stage.addConstraint(ctrl_fn, constraints.BoxConstraint(u_min, u_max))
67
68print(stage)
69
70
71nsteps = 20
72problem = aligator.TrajOptProblem(x0, nu, space, term_cost)
73
74for i in range(nsteps):
75 problem.addStage(stage)
76
77xtar2 = 0.1 * np.ones(nx)
78if args.term_cstr:
79 term_fun = aligator.StateErrorResidual(space, nu, xtar2)
80 problem.addTerminalConstraint(term_fun, constraints.EqualityConstraintSet())
81
82
83mu_init = 1e-1 if args.bounds else 1e-4
84verbose = aligator.VerboseLevel.VERBOSE
85tol = 1e-8
86solver = aligator.SolverProxDDP(tol, mu_init, verbose=verbose)
87
88his_cb = aligator.HistoryCallback(solver)
89solver.registerCallback("his", his_cb)
90print("Registered callbacks:", solver.getCallbackNames().tolist())
91solver.max_iters = 20
92solver.rollout_type = aligator.ROLLOUT_LINEAR
93
94u0 = np.zeros(nu)
95us_i = [u0] * nsteps
97prob_data = aligator.TrajOptData(problem)
98problem.evaluate(xs_i, us_i, prob_data)
99
100solver.setup(problem)
101solver.run(problem, xs_i, us_i)
102res = solver.results
103ws = solver.workspace
104print(res)
105
106plt.subplot(121)
107fig1: plt.Figure = plt.gcf()
108fig1.set_figwidth(6.4)
109fig1.set_figheight(3.6)
110
111lstyle = {"lw": 0.9, "marker": ".", "markersize": 5}
112trange = np.arange(nsteps + 1)
113plt.plot(res.xs, ls="-", **lstyle)
114
115if args.term_cstr:
116 plt.hlines(
117 xtar2,
118 *trange[[0, -1]],
119 ls="-",
120 lw=1.0,
121 colors="k",
122 alpha=0.4,
123 label=r"$x_\mathrm{tar}$",
124 )
125plt.hlines(
126 0.0,
127 *trange[[0, -1]],
128 ls=":",
129 lw=0.6,
130 colors="k",
131 alpha=0.4,
132 label=r"$x=0$",
133)
134plt.title("State trajectory $x(t)$")
135plt.xlabel("Time $i$")
136plt.legend(frameon=False)
137
138plt.subplot(122)
139plt.plot(res.us, **lstyle)
140if args.bounds:
141 plt.hlines(
142 np.concatenate([u_min, u_max]),
143 *trange[[0, -1]],
144 ls="-",
145 colors="tab:red",
146 lw=1.8,
147 alpha=0.4,
148 label=r"$\bar{u}$",
149 )
150plt.xlabel("Time $i$")
151plt.title("Controls $u(t)$")
152plt.legend(frameon=False, loc="lower right")
153plt.tight_layout()
154
155
156fig2: plt.Figure = plt.figure(figsize=(6.4, 3.6))
157ax: plt.Axes = fig2.add_subplot()
158niter = res.num_iters
159ax.hlines(
160 tol,
161 0,
162 niter,
163 colors="k",
164 linestyles="-",
165 linewidth=2.0,
166)
167plot_convergence(his_cb, ax, res, show_al_iters=True)
168ax.set_title("Convergence (constrained LQR)")
169fig2.tight_layout()
170fig_dicts = {"traj": fig1, "conv": fig2}
171
172for name, _fig in fig_dicts.items():
173 _fig: plt.Figure
174 _fig.savefig(f"assets/{TAG}_{name}.pdf")
175
176plt.show()
math_types< Scalar >::VectorOfVectors rollout(const std::vector< xyz::polymorphic< DynamicsModelTpl< Scalar > > > &dyn_models, const typename math_types< Scalar >::VectorXs &x0, const typename math_types< Scalar >::VectorOfVectors &us, typename math_types< Scalar >::VectorOfVectors &xout)
Perform a rollout of the supplied dynamical models.