11#include <Eigen/Sparse>
12#include <smooth/diff.hpp>
17namespace smooth::feedback {
19template<u
int8_t Deriv>
43 bool allocated{
false};
55 Eigen::SparseMatrix<double>
dF;
70 Eigen::SparseMatrix<double>
d2F;
79template<u
int8_t Deriv>
83 if constexpr (Deriv >= 1) {
set_zero(mv.dF); }
84 if constexpr (Deriv >= 2) {
set_zero(mv.d2F); }
114template<u
int8_t Deriv, diff::Type DT = diff::Type::Default>
122 std::ranges::range
auto && xs,
123 std::ranges::range
auto && us,
126 using utils::zip, std::views::iota;
127 using X = PlainObject<std::decay_t<std::ranges::range_value_t<
decltype(xs)>>>;
128 using U = PlainObject<std::decay_t<std::ranges::range_value_t<
decltype(us)>>>;
130 static constexpr auto nx = Dof<X>;
131 static constexpr auto nu = Dof<U>;
132 static constexpr auto nf = Dof<std::invoke_result_t<
decltype(f),
double, X, U>>;
134 static_assert(nx > -1,
"State dimension must be static");
135 static_assert(nu > -1,
"Input dimension must be static");
136 static_assert(nf > -1,
"Output size must be static");
138 const auto N = m.N_colloc();
140 const Eigen::Index numOuts = nf * N;
141 const Eigen::Index numVars = 2 + nx * (N + 1) + nu * N;
145 if (!out.allocated) {
146 out.F.resize(numOuts);
148 if constexpr (Deriv >= 1) {
149 Eigen::VectorXi pattern = Eigen::VectorXi::Zero(numVars);
150 pattern.segment(0, 1).setConstant(numOuts);
151 pattern.segment(1, 1).setConstant(numOuts);
152 pattern.segment(2, nx * N).setConstant(nf);
153 pattern.segment(2 + nx * (N + 1), nu * N).setConstant(nf);
155 out.dF.resize(numOuts, numVars);
156 out.dF.reserve(pattern);
159 if constexpr (Deriv >= 2u) {
160 Eigen::VectorXi pattern = Eigen::VectorXi::Zero(numVars);
161 pattern.segment(0, 1).setConstant(1);
162 pattern.segment(1, 1).setConstant(2);
163 for (
auto i = 0u; i < N; ++i) {
164 for (
auto j = 0u; j < nx; ++j) {
165 pattern(2 + nx * i + j) = 2 + (j + 1);
168 for (
auto i = 0u; i < N; ++i) {
169 for (
auto j = 0u; j < nu; ++j) {
170 pattern(2 + nx * (N + 1) + nu * i + j) = 2 + nx + (j + 1);
173 out.d2F.resize(numVars, numVars);
174 out.d2F.reserve(pattern);
177 out.allocated =
true;
182 for (
const auto & [i, tau, w_quad, x, u] : zip(iota(0u, N), m.all_nodes(), m.all_weights(), xs, us)) {
183 const double ti = t0 + (tf - t0) * tau;
187 const double mtau = 1. - tau;
188 const double w = scale ? w_quad : 1.;
190 const auto fvals = diff::dr<Deriv, DT>(f, wrt(ti, xi, ui));
192 out.F.segment(i * nf, nf) = w * std::get<0>(fvals);
194 if constexpr (Deriv >= 1u) {
195 const auto & df = std::get<1>(fvals);
196 const auto row0 = nf * i;
198 block_add(out.dF, row0, 0, df.middleCols(0, 1), w * mtau);
199 block_add(out.dF, row0, 1, df.middleCols(0, 1), w * tau);
200 block_add(out.dF, row0, 2 + i * nx, df.middleCols(1, nx), w);
201 block_add(out.dF, row0, 2 + nx * (N + 1) + nu * i, df.middleCols(1 + nx, nu), w);
203 if constexpr (Deriv >= 2u) {
204 assert(out.lambda.size() == numOuts);
206 const auto & d2f = std::get<2>(fvals);
208 for (
auto j = 0u; j < nf; ++j) {
209 const double wl = w * out.lambda(row0 + j);
214 const auto x_d = 2 + i * nx;
215 const auto u_d = 2 + (N + 1) * nx + i * nu;
218 const auto b_s = (1 + nx + nu) * j;
221 const auto u_s = 1 + nx;
225 block_add(out.d2F, t0_d, t0_d, d2f.block(t_s, b_s + t_s, 1, 1 ), wl * mtau * mtau,
true);
226 block_add(out.d2F, t0_d, tf_d, d2f.block(t_s, b_s + t_s, 1, 1 ), wl * mtau * tau,
true);
227 block_add(out.d2F, t0_d, x_d, d2f.block(t_s, b_s + x_s, 1, nx), wl * mtau,
true);
228 block_add(out.d2F, t0_d, u_d, d2f.block(t_s, b_s + u_s, 1, nu), wl * mtau,
true);
231 block_add(out.d2F, tf_d, tf_d, d2f.block(t_s, b_s + t_s, 1, 1 ), wl * tau * tau,
true);
232 block_add(out.d2F, tf_d, x_d, d2f.block(t_s, b_s + x_s, 1, nx), wl * tau,
true);
233 block_add(out.d2F, tf_d, u_d, d2f.block(t_s, b_s + u_s, 1, nu), wl * tau,
true);
236 block_add(out.d2F, x_d, x_d, d2f.block(x_s, b_s + x_s, nx, nx), wl,
true);
237 block_add(out.d2F, x_d, u_d, d2f.block(x_s, b_s + u_s, nx, nu), wl,
true);
240 block_add(out.d2F, u_d, u_d, d2f.block(u_s, b_s + u_s, nu, nu), wl,
true);
273template<u
int8_t Deriv, diff::Type DT = diff::Type::Default>
281 std::ranges::range
auto && xs,
282 std::ranges::range
auto && us)
284 using utils::zip, std::views::iota;
285 using X = PlainObject<std::decay_t<std::ranges::range_value_t<
decltype(xs)>>>;
286 using U = PlainObject<std::decay_t<std::ranges::range_value_t<
decltype(us)>>>;
288 static constexpr auto nx = Dof<X>;
289 static constexpr auto nu = Dof<U>;
290 static constexpr auto nf = Dof<std::invoke_result_t<
decltype(f),
double, X, U>>;
292 static_assert(nx > -1,
"State dimension must be static");
293 static_assert(nu > -1,
"Input dimension must be static");
294 static_assert(nf > -1,
"Output size must be static");
296 const auto N = m.N_colloc();
298 const Eigen::Index numOuts = nf;
299 const Eigen::Index numVars = 2 + nx * (N + 1) + nu * N;
303 if (!out.allocated) {
304 out.F.resize(numOuts);
306 if constexpr (Deriv >= 1) {
307 Eigen::VectorXi pattern = Eigen::VectorXi::Constant(numVars, numOuts);
308 pattern.segment(2 + nx * N, nx).setZero();
310 out.dF.resize(numOuts, numVars);
311 out.dF.reserve(pattern);
314 if constexpr (Deriv >= 2) {
315 Eigen::VectorXi pattern = Eigen::VectorXi::Zero(numVars);
318 for (
auto i = 0u; i < N; ++i) {
319 for (
auto j = 0u; j < nx; ++j) {
320 pattern(2 + nx * i + j) = 2 + (j + 1);
323 for (
auto i = 0u; i < N; ++i) {
324 for (
auto j = 0u; j < nu; ++j) {
325 pattern(2 + nx * (N + 1) + nu * i + j) = 2 + nx + (j + 1);
329 out.d2F.resize(numVars, numVars);
330 out.d2F.reserve(pattern.replicate(numOuts, 1).reshaped());
333 out.allocated =
true;
338 for (
const auto & [i, tau, w, x, u] : zip(iota(0u, N), m.all_nodes(), m.all_weights(), xs, us)) {
339 const double ti = t0 + (tf - t0) * tau;
343 const double mtau = 1. - tau;
345 const auto fvals = diff::dr<Deriv, DT>(f, wrt(ti, xi, ui));
347 const auto & fval = std::get<0>(fvals);
348 out.F.noalias() += w * (tf - t0) * fval;
350 if constexpr (Deriv >= 1u) {
351 const auto & df = std::get<1>(fvals);
353 block_add(out.dF, 0, 0, df.middleCols(0, 1), w * (tf - t0) * mtau);
356 block_add(out.dF, 0, 1, df.middleCols(0, 1), w * (tf - t0) * tau);
359 block_add(out.dF, 0, 2 + i * nx, df.middleCols(1, nx), w * (tf - t0));
361 block_add(out.dF, 0, 2 + nx * (N + 1) + nu * i, df.middleCols(1 + nx, nu), w * (tf - t0));
363 if constexpr (Deriv >= 2u) {
364 assert(out.lambda.size() == numOuts);
366 const auto & d2f = std::get<2>(fvals);
368 for (
auto j = 0u; j < nf; ++j) {
369 const double wl = w * out.lambda(j);
372 const auto b_s = (1 + nx + nu) * j;
375 const auto u_s = 1 + nx;
380 const auto x_d = 2 + nx * i;
381 const auto u_d = 2 + nx * (N + 1) + nu * i;
385 block_add(out.d2F, t0_d, t0_d, d2f.block(t_s, b_s + t_s, 1, 1), wl * (tf - t0) * mtau * mtau,
true);
386 block_add(out.d2F, t0_d, t0_d, df.block(j, t_s, 1, 1), -wl * 2 * mtau,
true);
388 block_add(out.d2F, t0_d, tf_d, d2f.block(t_s, b_s + t_s, 1, 1), wl * (tf - t0) * mtau * tau,
true);
389 block_add(out.d2F, t0_d, tf_d, df.block(j, t_s, 1, 1), wl * (1 - 2 * tau),
true);
391 block_add(out.d2F, t0_d, x_d, d2f.block(t_s, b_s + x_s, 1, nx), wl * (tf - t0) * mtau,
true);
392 block_add(out.d2F, t0_d, x_d, df.block(j, x_s, 1, nx), -wl,
true);
394 block_add(out.d2F, t0_d, u_d, d2f.block(t_s, b_s + u_s, 1, nu), wl * (tf - t0) * mtau,
true);
395 block_add(out.d2F, t0_d, u_d, df.block(j, u_s, 1, nu), -wl,
true);
398 block_add(out.d2F, tf_d, tf_d, d2f.block(t_s, b_s + t_s, 1, 1), wl * (tf - t0) * tau * tau,
true);
399 block_add(out.d2F, tf_d, tf_d, df.block(j, t_s, 1, 1), wl * 2 * tau,
true);
401 block_add(out.d2F, tf_d, x_d, d2f.block(t_s, b_s + x_s, 1, nx), wl * (tf - t0) * tau,
true);
402 block_add(out.d2F, tf_d, x_d, df.block(j, x_s, 1, nx), wl,
true);
404 block_add(out.d2F, tf_d, u_d, d2f.block(t_s, b_s + u_s, 1, nu), wl * (tf - t0) * tau,
true);
405 block_add(out.d2F, tf_d, u_d, df.block(j, u_s, 1, nu), wl,
true);
408 block_add(out.d2F, x_d, x_d, d2f.block(x_s, b_s + x_s, nx, nx), wl * (tf - t0),
true);
410 block_add(out.d2F, x_d, u_d, d2f.block(x_s, b_s + u_s, nx, nu), wl * (tf - t0),
true);
413 block_add(out.d2F, u_d, u_d, d2f.block(u_s, b_s + u_s, nu, nu), wl * (tf - t0),
true);
450template<u
int8_t Deriv, diff::Type DT = diff::Type::Default>
458 std::ranges::range
auto && xs,
459 std::ranges::range
auto && us)
461 using utils::zip, std::views::iota;
462 using X = PlainObject<std::decay_t<std::ranges::range_value_t<
decltype(xs)>>>;
463 using U = PlainObject<std::decay_t<std::ranges::range_value_t<
decltype(us)>>>;
465 static constexpr auto nx = Dof<X>;
466 static constexpr auto nu = Dof<U>;
467 static constexpr auto nf = Dof<std::invoke_result_t<
decltype(f),
double, X, U>>;
469 static_assert(nx > -1,
"State dimension must be static");
470 static_assert(nu > -1,
"Input dimension must be static");
471 static_assert(nf > -1,
"Output size must be static");
472 static_assert(nx == nf,
"Output dimension must be same as state dimension");
474 const auto N = m.N_colloc();
475 const Eigen::Index numOuts = nx * N;
476 const Eigen::Index numVars = 2 + nx * (N + 1) + nu * N;
480 if (!out.allocated) {
481 out.F.resize(numOuts);
483 if constexpr (Deriv >= 1) {
484 Eigen::VectorXi pattern = Eigen::VectorXi::Zero(numVars);
485 pattern(0) = numOuts;
486 pattern(1) = numOuts;
490 for (
auto ival = 0u; ival < m.N_ivals(); ++ival) {
491 const std::size_t K = m.N_colloc_ival(ival);
492 pattern.segment(idx0, K * nx) += Eigen::VectorXi::Constant(K * nx, nx + K - 1);
493 pattern.segment(idx0 + K * nx, nx) += Eigen::VectorXi::Constant(nx, K);
498 pattern.segment(2 + nx * (N + 1), nu * N).setConstant(nx);
500 out.dF.resize(numOuts, numVars);
501 out.dF.reserve(pattern);
504 if constexpr (Deriv >= 2) {
505 Eigen::VectorXi pattern = Eigen::VectorXi::Zero(numVars);
508 for (
auto i = 0u; i < N; ++i) {
509 for (
auto j = 0u; j < nx; ++j) {
510 pattern(2 + nx * i + j) = 2 + (j + 1);
513 for (
auto i = 0u; i < N; ++i) {
514 for (
auto j = 0u; j < nu; ++j) {
515 pattern(2 + nx * (N + 1) + nu * i + j) = 2 + nx + (j + 1);
519 out.d2F.resize(numVars, numVars);
520 out.d2F.reserve(pattern);
523 out.allocated =
true;
534 for (
const auto & [i, tau, w, x, u] : zip(iota(0u, N), m.all_nodes(), m.all_weights(), xs, us)) {
535 const double ti = t0 + (tf - t0) * tau;
539 const auto row0 = nx * i;
540 const double mtau = 1. - tau;
542 const auto fvals = diff::dr<Deriv, DT>(f, wrt(ti, xi, ui));
543 const auto fval = std::get<0>(fvals);
545 out.F.segment(row0, nx) += w * (tf - t0) * fval;
547 if constexpr (Deriv >= 1) {
548 const auto & df = std::get<1>(fvals);
553 block_add(out.dF, row0, 0, df.col(0), w * (tf - t0) * mtau);
556 block_add(out.dF, row0, 1, df.col(0), w * (tf - t0) * tau);
558 block_add(out.dF, row0, 2 + nx * i, df.middleCols(1, nx), w * (tf - t0));
560 block_add(out.dF, row0, 2 + nx * (N + 1) + nu * i, df.middleCols(1 + nx, nu), w * (tf - t0));
563 if constexpr (Deriv >= 2) {
564 assert(out.lambda.size() == numOuts);
566 const auto & d2f = std::get<2>(fvals);
568 for (
auto j = 0u; j < nx; ++j) {
569 const double wl = w * out.lambda(row0 + j);
574 const auto x_d = 2 + i * nx;
575 const auto u_d = 2 + (N + 1) * nx + i * nu;
578 const auto b_s = (1 + nx + nu) * j;
581 const auto u_s = 1 + nx;
585 block_add(out.d2F, t0_d, t0_d, d2f.block(t_s, b_s + t_s, 1, 1 ), wl * (tf - t0) * mtau * mtau,
true);
586 block_add(out.d2F, t0_d, t0_d, df.block(j, t_s, 1, 1 ), -wl * 2 * mtau,
true);
588 block_add(out.d2F, t0_d, tf_d, d2f.block(t_s, b_s + t_s, 1, 1 ), wl * (tf - t0) * mtau * tau,
true);
589 block_add(out.d2F, t0_d, tf_d, df.block(j, t_s, 1, 1 ), wl * (1. - 2 * tau),
true);
591 block_add(out.d2F, t0_d, x_d, d2f.block(t_s, b_s + x_s, 1, nx), wl * (tf - t0) * mtau,
true);
592 block_add(out.d2F, t0_d, x_d, df.block(j, x_s, 1, nx), -wl,
true);
594 block_add(out.d2F, t0_d, u_d, d2f.block(t_s, b_s + u_s, 1, nu), wl * (tf - t0) * mtau,
true);
595 block_add(out.d2F, t0_d, u_d, df.block(j, u_s, 1, nu), -wl,
true);
598 block_add(out.d2F, tf_d, tf_d, d2f.block(t_s, b_s + t_s, 1, 1), wl * (tf - t0) * tau * tau,
true);
599 block_add(out.d2F, tf_d, tf_d, df.block(j, t_s, 1, 1), wl * 2 * tau,
true);
601 block_add(out.d2F, tf_d, x_d, d2f.block(t_s, b_s + x_s, 1, nx), wl * (tf - t0) * tau,
true);
602 block_add(out.d2F, tf_d, x_d, df.block(j, x_s, 1, nx), wl,
true);
604 block_add(out.d2F, tf_d, u_d, d2f.block(t_s, b_s + u_s, 1, nu), wl * (tf - t0) * tau,
true);
605 block_add(out.d2F, tf_d, u_d, df.block(j, u_s, 1, nu), wl,
true);
608 block_add(out.d2F, x_d, x_d, d2f.block(x_s, b_s + x_s, nx, nx), wl * (tf - t0),
true);
610 block_add(out.d2F, x_d, u_d, d2f.block(x_s, b_s + u_s, nx, nu), wl * (tf - t0),
true);
613 block_add(out.d2F, u_d, u_d, d2f.block(u_s, b_s + u_s, nu, nu), wl * (tf - t0),
true);
624 auto Nival = m.N_colloc_ival(ival);
626 for (
const auto & [i, x] : zip(iota(0u, N + 1), xs)) {
627 if (i == ival_idx0 + Nival) {
629 const auto [alpha, Dus] = m.interval_diffmat_unscaled(ival);
630 for (
const auto & [j, w] : zip(iota(0u, Nival), m.interval_weights(ival))) {
631 const auto row0 = (ival_idx0 + j) * nx;
632 const double coef = -w * alpha * Dus(i - ival_idx0, j);
634 out.F.segment(row0, nx) += coef * x;
636 if constexpr (Deriv >= 1) {
638 for (
auto k = 0u; k < nx; ++k) { out.dF.coeffRef(row0 + k, 2 + nx * i + k) += coef; }
644 if (ival < m.N_ivals()) {
646 Nival = m.N_colloc_ival(ival);
651 const auto [alpha, Dus] = m.interval_diffmat_unscaled(ival);
652 for (
const auto & [j, w] : zip(iota(0u, Nival), m.interval_weights(ival))) {
653 const auto row0 = (ival_idx0 + j) * nx;
654 const double coef = -w * alpha * Dus(i - ival_idx0, j);
656 out.F.segment(row0, nx) += coef * x;
658 if constexpr (Deriv >= 1) {
660 for (
auto k = 0u; k < nx; ++k) { out.dF.coeffRef(row0 + k, 2 + nx * i + k) += coef; }
MeshType is a specialization of Mesh.
Refinable Legendre-Gauss-Radau mesh of time interval [0, 1].
void mesh_eval(MeshValue< Deriv > &out, const MeshType auto &m, auto &f, const double t0, const double tf, std::ranges::range auto &&xs, std::ranges::range auto &&us, bool scale=false)
Evaluate function over a mesh.
void mesh_dyn(MeshValue< Deriv > &out, const MeshType auto &m, auto &f, const double t0, const double tf, std::ranges::range auto &&xs, std::ranges::range auto &&us)
Evaluate dynamic constraints over a mesh (differentiation version).
void mesh_integrate(MeshValue< Deriv > &out, const MeshType auto &m, auto &f, const double t0, const double tf, std::ranges::range auto &&xs, std::ranges::range auto &&us)
Evaluate integral over a mesh.
void set_zero(MeshValue< Deriv > &mv)
Reset MeshValue to zeros.
void block_add(Eigen::SparseMatrix< double, Options > &dest, Eigen::Index row0, Eigen::Index col0, const Source &source, double scale=1, bool upper_only=false)
Add block into a sparse matrix.
Eigen::VectorXd F
Function value (size M)
Eigen::SparseMatrix< double > dF
Size M x numVar.
Eigen::SparseMatrix< double > d2F
Size numVar x numVar.
Eigen::VectorXd lambda
Multipliers (must be set before)