Implements BsplineTrajectory::AsLinearInControlPoints()

Resolves RobotLocomotion#22533.

common/trajectories/bspline_trajectory.h

@@ -47,6 +47,42 @@ class BsplineTrajectory final : public trajectories::Trajectory<T> {
 
   ~BsplineTrajectory() final;
 
+  /** Supports writing optimizations using the control points as decision
+   variables.  This method returns the matrix, `M`, defining the control points
+   of the `order` derivative in the form:
+   <pre>
+   derivative.control_points() = this.control_points() * M
+   </pre>
+   For instance, since we have
+   <pre>
+   derivative.control_points().col(k) = this.control_points() * M.col(k),
+   </pre>
+   constraining the kth control point of the `n`th derivative to be in [ub,
+   lb], could be done with:
+   @code
+   auto M = curve.AsLinearInControlPoints(n);
+   for (int i=0; i<curve.rows(); ++i) {
+     auto c = std::make_shared<solvers::LinearConstraint>(
+       M.col(k).transpose(), Vector1d(lb(i)), Vector1d(ub(i)));
+     prog.AddConstraint(c, curve.row(i).transpose());
+   }
+   @endcode
+   Iterating over the rows of the control points is the natural sparsity pattern
+   here (since `M` is the same for all rows).  For instance, we also have
+   <pre>
+   derivative.control_points().row(k).T = M.T * this.control_points().row(k).T,
+   </pre> or
+   <pre>
+   vec(derivative.control_points().T) = blockMT * vec(this.control_points().T),
+   blockMT = [ M.T,   0, .... 0 ]
+             [   0, M.T, 0, ... ]
+             [      ...         ]
+             [  ...    , 0, M.T ].
+   </pre>
+   @pre derivative_order >= 0. */
+  Eigen::SparseMatrix<double> AsLinearInControlPoints(
+      int derivative_order = 1) const;
+
   /** Evaluates the BsplineTrajectory at the given time t.
   @param t The time at which to evaluate the %BsplineTrajectory.
   @return The matrix of evaluated values.