001// Copyright (c) FIRST and other WPILib contributors.
002// Open Source Software; you can modify and/or share it under the terms of
003// the WPILib BSD license file in the root directory of this project.
004
005package edu.wpi.first.math.trajectory;
006
007import java.util.Objects;
008
009/**
010 * A exponential curve-shaped velocity profile.
011 *
012 * <p>While this class can be used for a profiled movement from start to finish, the intended usage
013 * is to filter a reference's dynamics based on state-space model dynamics. To compute the reference
014 * obeying this constraint, do the following.
015 *
016 * <p>Initialization:
017 *
018 * <pre><code>
019 * ExponentialProfile.Constraints constraints =
020 *   ExponentialProfile.Constraints.fromCharacteristics(kMaxV, kV, kA);
021 * ExponentialProfile.State previousProfiledReference =
022 *   new ExponentialProfile.State(initialReference, 0.0);
023 * ExponentialProfile profile = new ExponentialProfile(constraints);
024 * </code></pre>
025 *
026 * <p>Run on update:
027 *
028 * <pre><code>
029 * previousProfiledReference =
030 * profile.calculate(timeSincePreviousUpdate, previousProfiledReference, unprofiledReference);
031 * </code></pre>
032 *
033 * <p>where `unprofiledReference` is free to change between calls. Note that when the unprofiled
034 * reference is within the constraints, `calculate()` returns the unprofiled reference unchanged.
035 *
036 * <p>Otherwise, a timer can be started to provide monotonic values for `calculate()` and to
037 * determine when the profile has completed via `isFinished()`.
038 */
039public class ExponentialProfile {
040  private final Constraints m_constraints;
041
042  /** Profile timing. */
043  public static class ProfileTiming {
044    /** Profile inflection time. */
045    public final double inflectionTime;
046
047    /** Total profile time. */
048    public final double totalTime;
049
050    /**
051     * Constructs a ProfileTiming.
052     *
053     * @param inflectionTime Profile inflection time.
054     * @param totalTime Total profile time.
055     */
056    protected ProfileTiming(double inflectionTime, double totalTime) {
057      this.inflectionTime = inflectionTime;
058      this.totalTime = totalTime;
059    }
060
061    /**
062     * Decides if the profile is finished by time t.
063     *
064     * @param t The time since the beginning of the profile.
065     * @return if the profile is finished at time t.
066     */
067    public boolean isFinished(double t) {
068      return t >= inflectionTime;
069    }
070  }
071
072  /** Profile constraints. */
073  public static class Constraints {
074    /** Maximum unsigned input voltage. */
075    public final double maxInput;
076
077    /** The State-Space 1x1 system matrix. */
078    public final double A;
079
080    /** The State-Space 1x1 input matrix. */
081    public final double B;
082
083    /**
084     * Constructs constraints for an ExponentialProfile.
085     *
086     * @param maxInput maximum unsigned input voltage
087     * @param A The State-Space 1x1 system matrix.
088     * @param B The State-Space 1x1 input matrix.
089     */
090    private Constraints(double maxInput, double A, double B) {
091      this.maxInput = maxInput;
092      this.A = A;
093      this.B = B;
094    }
095
096    /**
097     * Computes the max achievable velocity for an Exponential Profile.
098     *
099     * @return The seady-state velocity achieved by this profile.
100     */
101    public double maxVelocity() {
102      return -maxInput * B / A;
103    }
104
105    /**
106     * Constructs constraints for an ExponentialProfile from characteristics.
107     *
108     * @param maxInput maximum unsigned input voltage
109     * @param kV The velocity gain.
110     * @param kA The acceleration gain.
111     * @return The Constraints object.
112     */
113    public static Constraints fromCharacteristics(double maxInput, double kV, double kA) {
114      return new Constraints(maxInput, -kV / kA, 1.0 / kA);
115    }
116
117    /**
118     * Constructs constraints for an ExponentialProfile from State-Space parameters.
119     *
120     * @param maxInput maximum unsigned input voltage
121     * @param A The State-Space 1x1 system matrix.
122     * @param B The State-Space 1x1 input matrix.
123     * @return The Constraints object.
124     */
125    public static Constraints fromStateSpace(double maxInput, double A, double B) {
126      return new Constraints(maxInput, A, B);
127    }
128  }
129
130  /** Profile state. */
131  public static class State {
132    /** The position at this state. */
133    public double position;
134
135    /** The velocity at this state. */
136    public double velocity;
137
138    /** Default constructor. */
139    public State() {}
140
141    /**
142     * Constructs a state within an exponential profile.
143     *
144     * @param position The position at this state.
145     * @param velocity The velocity at this state.
146     */
147    public State(double position, double velocity) {
148      this.position = position;
149      this.velocity = velocity;
150    }
151
152    @Override
153    public boolean equals(Object other) {
154      if (other instanceof State) {
155        State rhs = (State) other;
156        return this.position == rhs.position && this.velocity == rhs.velocity;
157      } else {
158        return false;
159      }
160    }
161
162    @Override
163    public int hashCode() {
164      return Objects.hash(position, velocity);
165    }
166  }
167
168  /**
169   * Constructs an ExponentialProfile.
170   *
171   * @param constraints The constraints on the profile, like maximum input.
172   */
173  public ExponentialProfile(Constraints constraints) {
174    m_constraints = constraints;
175  }
176
177  /**
178   * Calculates the position and velocity for the profile at a time t where the current state is at
179   * time t = 0.
180   *
181   * @param t How long to advance from the current state toward the desired state.
182   * @param current The current state.
183   * @param goal The desired state when the profile is complete.
184   * @return The position and velocity of the profile at time t.
185   */
186  public State calculate(double t, State current, State goal) {
187    var direction = shouldFlipInput(current, goal) ? -1 : 1;
188    var u = direction * m_constraints.maxInput;
189
190    var inflectionPoint = calculateInflectionPoint(current, goal, u);
191    var timing = calculateProfileTiming(current, inflectionPoint, goal, u);
192
193    if (t < 0) {
194      return new State(current.position, current.velocity);
195    } else if (t < timing.inflectionTime) {
196      return new State(
197          computeDistanceFromTime(t, u, current), computeVelocityFromTime(t, u, current));
198    } else if (t < timing.totalTime) {
199      return new State(
200          computeDistanceFromTime(t - timing.totalTime, -u, goal),
201          computeVelocityFromTime(t - timing.totalTime, -u, goal));
202    } else {
203      return new State(goal.position, goal.velocity);
204    }
205  }
206
207  /**
208   * Calculates the point after which the fastest way to reach the goal state is to apply input in
209   * the opposite direction.
210   *
211   * @param current The current state.
212   * @param goal The desired state when the profile is complete.
213   * @return The position and velocity of the profile at the inflection point.
214   */
215  public State calculateInflectionPoint(State current, State goal) {
216    var direction = shouldFlipInput(current, goal) ? -1 : 1;
217    var u = direction * m_constraints.maxInput;
218
219    return calculateInflectionPoint(current, goal, u);
220  }
221
222  /**
223   * Calculates the point after which the fastest way to reach the goal state is to apply input in
224   * the opposite direction.
225   *
226   * @param current The current state.
227   * @param goal The desired state when the profile is complete.
228   * @param input The signed input applied to this profile from the current state.
229   * @return The position and velocity of the profile at the inflection point.
230   */
231  private State calculateInflectionPoint(State current, State goal, double input) {
232    var u = input;
233
234    if (current.equals(goal)) {
235      return current;
236    }
237
238    var inflectionVelocity = solveForInflectionVelocity(u, current, goal);
239    var inflectionPosition = computeDistanceFromVelocity(inflectionVelocity, -u, goal);
240
241    return new State(inflectionPosition, inflectionVelocity);
242  }
243
244  /**
245   * Calculates the time it will take for this profile to reach the goal state.
246   *
247   * @param current The current state.
248   * @param goal The desired state when the profile is complete.
249   * @return The total duration of this profile.
250   */
251  public double timeLeftUntil(State current, State goal) {
252    var timing = calculateProfileTiming(current, goal);
253
254    return timing.totalTime;
255  }
256
257  /**
258   * Calculates the time it will take for this profile to reach the inflection point, and the time
259   * it will take for this profile to reach the goal state.
260   *
261   * @param current The current state.
262   * @param goal The desired state when the profile is complete.
263   * @return The timing information for this profile.
264   */
265  public ProfileTiming calculateProfileTiming(State current, State goal) {
266    var direction = shouldFlipInput(current, goal) ? -1 : 1;
267    var u = direction * m_constraints.maxInput;
268
269    var inflectionPoint = calculateInflectionPoint(current, goal, u);
270    return calculateProfileTiming(current, inflectionPoint, goal, u);
271  }
272
273  /**
274   * Calculates the time it will take for this profile to reach the inflection point, and the time
275   * it will take for this profile to reach the goal state.
276   *
277   * @param current The current state.
278   * @param inflectionPoint The inflection point of this profile.
279   * @param goal The desired state when the profile is complete.
280   * @param input The signed input applied to this profile from the current state.
281   * @return The timing information for this profile.
282   */
283  private ProfileTiming calculateProfileTiming(
284      State current, State inflectionPoint, State goal, double input) {
285    var u = input;
286
287    double inflectionT_forward;
288
289    // We need to handle 5 cases here:
290    //
291    // - Approaching -maxVelocity from below
292    // - Approaching -maxVelocity from above
293    // - Approaching maxVelocity from below
294    // - Approaching maxVelocity from above
295    // - At +-maxVelocity
296    //
297    // For cases 1 and 3, we want to subtract epsilon from the inflection point velocity.
298    // For cases 2 and 4, we want to add epsilon to the inflection point velocity.
299    // For case 5, we have reached inflection point velocity.
300    double epsilon = 1e-9;
301    if (Math.abs(Math.signum(input) * m_constraints.maxVelocity() - inflectionPoint.velocity)
302        < epsilon) {
303      double solvableV = inflectionPoint.velocity;
304      double t_to_solvable_v;
305      double x_at_solvable_v;
306      if (Math.abs(current.velocity - inflectionPoint.velocity) < epsilon) {
307        t_to_solvable_v = 0;
308        x_at_solvable_v = current.position;
309      } else {
310        if (Math.abs(current.velocity) > m_constraints.maxVelocity()) {
311          solvableV += Math.signum(u) * epsilon;
312        } else {
313          solvableV -= Math.signum(u) * epsilon;
314        }
315
316        t_to_solvable_v = computeTimeFromVelocity(solvableV, u, current.velocity);
317        x_at_solvable_v = computeDistanceFromVelocity(solvableV, u, current);
318      }
319
320      inflectionT_forward =
321          t_to_solvable_v
322              + Math.signum(input)
323                  * (inflectionPoint.position - x_at_solvable_v)
324                  / m_constraints.maxVelocity();
325    } else {
326      inflectionT_forward = computeTimeFromVelocity(inflectionPoint.velocity, u, current.velocity);
327    }
328
329    var inflectionT_backward = computeTimeFromVelocity(inflectionPoint.velocity, -u, goal.velocity);
330
331    return new ProfileTiming(inflectionT_forward, inflectionT_forward - inflectionT_backward);
332  }
333
334  /**
335   * Calculates the position reached after t seconds when applying an input from the initial state.
336   *
337   * @param t The time since the initial state.
338   * @param input The signed input applied to this profile from the initial state.
339   * @param initial The initial state.
340   * @return The distance travelled by this profile.
341   */
342  private double computeDistanceFromTime(double t, double input, State initial) {
343    var A = m_constraints.A;
344    var B = m_constraints.B;
345    var u = input;
346
347    return initial.position
348        + (-B * u * t + (initial.velocity + B * u / A) * (Math.exp(A * t) - 1)) / A;
349  }
350
351  /**
352   * Calculates the velocity reached after t seconds when applying an input from the initial state.
353   *
354   * @param t The time since the initial state.
355   * @param input The signed input applied to this profile from the initial state.
356   * @param initial The initial state.
357   * @return The distance travelled by this profile.
358   */
359  private double computeVelocityFromTime(double t, double input, State initial) {
360    var A = m_constraints.A;
361    var B = m_constraints.B;
362    var u = input;
363
364    return (initial.velocity + B * u / A) * Math.exp(A * t) - B * u / A;
365  }
366
367  /**
368   * Calculates the time required to reach a specified velocity given the initial velocity.
369   *
370   * @param velocity The goal velocity.
371   * @param input The signed input applied to this profile from the initial state.
372   * @param initial The initial velocity.
373   * @return The time required to reach the goal velocity.
374   */
375  private double computeTimeFromVelocity(double velocity, double input, double initial) {
376    var A = m_constraints.A;
377    var B = m_constraints.B;
378    var u = input;
379
380    return Math.log((A * velocity + B * u) / (A * initial + B * u)) / A;
381  }
382
383  /**
384   * Calculates the distance reached at the same time as the given velocity when applying the given
385   * input from the initial state.
386   *
387   * @param velocity The velocity reached by this profile
388   * @param input The signed input applied to this profile from the initial state.
389   * @param initial The initial state.
390   * @return The distance reached when the given velocity is reached.
391   */
392  private double computeDistanceFromVelocity(double velocity, double input, State initial) {
393    var A = m_constraints.A;
394    var B = m_constraints.B;
395    var u = input;
396
397    return initial.position
398        + (velocity - initial.velocity) / A
399        - B * u / (A * A) * Math.log((A * velocity + B * u) / (A * initial.velocity + B * u));
400  }
401
402  /**
403   * Calculates the velocity at which input should be reversed in order to reach the goal state from
404   * the current state.
405   *
406   * @param input The signed input applied to this profile from the current state.
407   * @param current The current state.
408   * @param goal The goal state.
409   * @return The inflection velocity.
410   */
411  private double solveForInflectionVelocity(double input, State current, State goal) {
412    var A = m_constraints.A;
413    var B = m_constraints.B;
414    var u = input;
415
416    var U_dir = Math.signum(u);
417
418    var position_delta = goal.position - current.position;
419    var velocity_delta = goal.velocity - current.velocity;
420
421    var scalar = (A * current.velocity + B * u) * (A * goal.velocity - B * u);
422    var power = -A / B / u * (A * position_delta - velocity_delta);
423
424    var a = -A * A;
425    var c = (B * B) * (u * u) + scalar * Math.exp(power);
426
427    if (-1e-9 < c && c < 0) {
428      // Numerical stability issue - the heuristic gets it right but c is around -1e-13
429      return 0;
430    }
431
432    return U_dir * Math.sqrt(-c / a);
433  }
434
435  /**
436   * Returns true if the profile should be inverted.
437   *
438   * <p>The profile is inverted if we should first apply negative input in order to reach the goal
439   * state.
440   *
441   * @param current The initial state (usually the current state).
442   * @param goal The desired state when the profile is complete.
443   */
444  @SuppressWarnings("UnnecessaryParentheses")
445  private boolean shouldFlipInput(State current, State goal) {
446    var u = m_constraints.maxInput;
447
448    var xf = goal.position;
449    var v0 = current.velocity;
450    var vf = goal.velocity;
451
452    var x_forward = computeDistanceFromVelocity(vf, u, current);
453    var x_reverse = computeDistanceFromVelocity(vf, -u, current);
454
455    if (v0 >= m_constraints.maxVelocity()) {
456      return xf < x_reverse;
457    }
458
459    if (v0 <= -m_constraints.maxVelocity()) {
460      return xf < x_forward;
461    }
462
463    var a = v0 >= 0;
464    var b = vf >= 0;
465    var c = xf >= x_forward;
466    var d = xf >= x_reverse;
467
468    return (a && !d) || (b && !c) || (!c && !d);
469  }
470}