|
| 1 | +// Copyright (c) FIRST and other WPILib contributors. |
| 2 | +// Open Source Software; you can modify and/or share it under the terms of |
| 3 | +// the WPILib BSD license file in the root directory of this project. |
| 4 | + |
| 5 | +#pragma once |
| 6 | + |
| 7 | +#include <benchmark/benchmark.h> |
| 8 | + |
| 9 | +#include <chrono> |
| 10 | + |
| 11 | +#include <sleipnir/optimization/problem.hpp> |
| 12 | + |
| 13 | +#include "wpi/math/system/NumericalIntegration.hpp" |
| 14 | + |
| 15 | +inline slp::VariableMatrix<double> CartPoleDynamics( |
| 16 | + const slp::VariableMatrix<double>& x, |
| 17 | + const slp::VariableMatrix<double>& u) { |
| 18 | + constexpr double m_c = 5.0; // Cart mass (kg) |
| 19 | + constexpr double m_p = 0.5; // Pole mass (kg) |
| 20 | + constexpr double l = 0.5; // Pole length (m) |
| 21 | + constexpr double g = 9.806; // Acceleration due to gravity (m/s²) |
| 22 | + |
| 23 | + auto q = x.segment(0, 2); |
| 24 | + auto qdot = x.segment(2, 2); |
| 25 | + auto theta = q[1]; |
| 26 | + auto thetadot = qdot[1]; |
| 27 | + |
| 28 | + // [ m_c + m_p m_p l cosθ] |
| 29 | + // M(q) = [m_p l cosθ m_p l² ] |
| 30 | + slp::VariableMatrix<double> M{{m_c + m_p, m_p * l * cos(theta)}, |
| 31 | + {m_p * l * cos(theta), m_p * std::pow(l, 2)}}; |
| 32 | + |
| 33 | + // [0 −m_p lθ̇ sinθ] |
| 34 | + // C(q, q̇) = [0 0 ] |
| 35 | + slp::VariableMatrix<double> C{{0, -m_p * l * thetadot * sin(theta)}, {0, 0}}; |
| 36 | + |
| 37 | + // [ 0 ] |
| 38 | + // τ_g(q) = [-m_p gl sinθ] |
| 39 | + slp::VariableMatrix<double> tau_g{{0}, {-m_p * g * l * sin(theta)}}; |
| 40 | + |
| 41 | + // [1] |
| 42 | + // B = [0] |
| 43 | + constexpr Eigen::Matrix<double, 2, 1> B{{1}, {0}}; |
| 44 | + |
| 45 | + // q̈ = M⁻¹(q)(τ_g(q) − C(q, q̇)q̇ + Bu) |
| 46 | + slp::VariableMatrix<double> qddot{4}; |
| 47 | + qddot.segment(0, 2) = qdot; |
| 48 | + qddot.segment(2, 2) = solve(M, tau_g - C * qdot + B * u); |
| 49 | + return qddot; |
| 50 | +} |
| 51 | + |
| 52 | +inline void BM_CartPole(benchmark::State& state) { |
| 53 | + using namespace std::chrono_literals; |
| 54 | + |
| 55 | + // NOLINTNEXTLINE(clang-analyzer-deadcode.DeadStores) |
| 56 | + for (auto _ : state) { |
| 57 | + constexpr std::chrono::duration<double> T = 5s; |
| 58 | + constexpr std::chrono::duration<double> dt = 50ms; |
| 59 | + constexpr int N = T / dt; |
| 60 | + |
| 61 | + constexpr double u_max = 20.0; // N |
| 62 | + constexpr double d_max = 2.0; // m |
| 63 | + |
| 64 | + constexpr Eigen::Vector<double, 4> x_initial{{0.0, 0.0, 0.0, 0.0}}; |
| 65 | + constexpr Eigen::Vector<double, 4> x_final{ |
| 66 | + {1.0, std::numbers::pi, 0.0, 0.0}}; |
| 67 | + |
| 68 | + slp::Problem<double> problem; |
| 69 | + |
| 70 | + // x = [q, q̇]ᵀ = [x, θ, ẋ, θ̇]ᵀ |
| 71 | + auto X = problem.decision_variable(4, N + 1); |
| 72 | + |
| 73 | + // Initial guess |
| 74 | + for (int k = 0; k < N + 1; ++k) { |
| 75 | + X(0, k).set_value( |
| 76 | + std::lerp(x_initial[0], x_final[0], static_cast<double>(k) / N)); |
| 77 | + X(1, k).set_value( |
| 78 | + std::lerp(x_initial[1], x_final[1], static_cast<double>(k) / N)); |
| 79 | + } |
| 80 | + |
| 81 | + // u = f_x |
| 82 | + auto U = problem.decision_variable(1, N); |
| 83 | + |
| 84 | + // Initial conditions |
| 85 | + problem.subject_to(X.col(0) == x_initial); |
| 86 | + |
| 87 | + // Final conditions |
| 88 | + problem.subject_to(X.col(N) == x_final); |
| 89 | + |
| 90 | + // Cart position constraints |
| 91 | + problem.subject_to(X.row(0) >= 0.0); |
| 92 | + problem.subject_to(X.row(0) <= d_max); |
| 93 | + |
| 94 | + // Input constraints |
| 95 | + problem.subject_to(U >= -u_max); |
| 96 | + problem.subject_to(U <= u_max); |
| 97 | + |
| 98 | + // Dynamics constraints - RK4 integration |
| 99 | + for (int k = 0; k < N; ++k) { |
| 100 | + problem.subject_to(X.col(k + 1) == |
| 101 | + wpi::math::RK4<decltype(CartPoleDynamics), |
| 102 | + slp::VariableMatrix<double>, |
| 103 | + slp::VariableMatrix<double>>( |
| 104 | + CartPoleDynamics, X.col(k), U.col(k), dt)); |
| 105 | + } |
| 106 | + |
| 107 | + // Minimize sum squared inputs |
| 108 | + slp::Variable J = 0.0; |
| 109 | + for (int k = 0; k < N; ++k) { |
| 110 | + J += U.col(k).T() * U.col(k); |
| 111 | + } |
| 112 | + problem.minimize(J); |
| 113 | + |
| 114 | + problem.solve({.diagnostics = true}); |
| 115 | + } |
| 116 | +} |
0 commit comments