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