1 #ifndef SimTK_SIMMATH_INTEGRATOR_H_ 2 #define SimTK_SIMMATH_INTEGRATOR_H_ 3 4 /* -------------------------------------------------------------------------- * 5 * Simbody(tm): SimTKmath * 6 * -------------------------------------------------------------------------- * 7 * This is part of the SimTK biosimulation toolkit originating from * 8 * Simbios, the NIH National Center for Physics-Based Simulation of * 9 * Biological Structures at Stanford, funded under the NIH Roadmap for * 10 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. * 11 * * 12 * Portions copyright (c) 2006-12 Stanford University and the Authors. * 13 * Authors: Michael Sherman * 14 * Contributors: * 15 * * 16 * Licensed under the Apache License, Version 2.0 (the "License"); you may * 17 * not use this file except in compliance with the License. You may obtain a * 18 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * 19 * * 20 * Unless required by applicable law or agreed to in writing, software * 21 * distributed under the License is distributed on an "AS IS" BASIS, * 22 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * 23 * See the License for the specific language governing permissions and * 24 * limitations under the License. * 25 * -------------------------------------------------------------------------- */ 26 27 /** @file 28 * This is the header file that user code should include to pick up the 29 * SimTK Simmath numerical integration tools. 30 */ 31 32 #include "SimTKcommon.h" 33 #include "simmath/internal/common.h" 34 35 namespace SimTK { 36 37 38 class IntegratorRep; 39 40 /** An Integrator is an object that can advance the State of a System 41 through time. This is an abstract class. Subclasses implement a variety of 42 different integration methods, which vary in their speed, accuracy, stability, 43 and various other properties. 44 45 <h3>Usage</h3> 46 47 An Integrator is most often used in combination with a TimeStepper. The 48 TimeStepper automates much of the work of using an Integrator: invoking it 49 repeatedly to advance time, calling event handlers, reinitializing the 50 Integrator as necessary, and so on. A typical use of an Integrator generally 51 resembles the following: 52 53 @code 54 // Instantiate an Integrator subclass which is appropriate for your problem. 55 VerletIntegrator integ(system); 56 // Set configuration options on the Integrator. 57 integ.setAccuracy(1e-3); 58 // Create a TimeStepper and use it to run the simulation. 59 TimeStepper stepper(system, integ); 60 stepper.initialize(initialState); 61 stepper.stepTo(finalTime); 62 @endcode 63 64 <h3>Mathematical Overview</h3> 65 66 Given a continuous system of differential equations for state variables y, and 67 optionally a manifold (set of algebraic equations) on which the solution must 68 lie, an Integrator object will advance that system through time. If the full 69 system is continuous, this is sufficient. However, most interesting systems 70 consist of both continuous and discrete equations, in which case the overall 71 time advancement is handled by a TimeStepper object which will use an 72 Integrator as an "inner loop" for advancing the system across continuous 73 intervals between discrete updates. In that case, in addition to continuous 74 state variables y there will be a set of discrete variables d which are held 75 constant during an integration interval. 76 77 The continuous part of the system is an ODE-on-a-manifold system suitable for 78 solution via coordinate projection[1], structured like this: 79 (1) y' = f(d;t,y) differential equations 80 (2) 0 = c(d;t,y) algebraic equations (manifold is c=0) 81 (3) e = e(d;t,y) event triggers (watch for zero crossings) 82 with initial conditions t0,y0 such that c=0. The "d;" is a reminder that the 83 overall system is dependent on discrete variables d as well as y, but d cannot 84 change during a continuous interval. 85 86 By "ODE on a manifold" we mean that the ODE (1) automatically satisfies the 87 condition that IF c==0, THEN c'=0, where 88 c'=partial(c)/partial(t) + [partial(c)/partial(y)]*y'. 89 This is a less stringent condition than an ODE with "first integral" invariant, 90 in which c'=0 regardless of whether c=0. 91 92 [1] Hairer, Lubich, Wanner, "Geometric Numerical Integration: 93 Structure-Preserving Algorithms for Ordinary Differential Equations", 2nd ed., 94 section IV.4, pg 109ff, Springer, 2006. 95 96 The discrete variables d are updated by the time stepper upon occurrence of 97 specific events, which terminate a continuous interval. The Integrator detects 98 these using a set of scalar-valued event trigger functions shown as equation 99 (3) above. An event trigger function for a particular event must be designed so 100 that it has a zero crossing when the event occurs. The integrator can thus 101 watch for sign changes in event triggers and terminate the current step when a 102 zero crossing occurs, notifying the system and giving it a chance to handle the 103 event; that is, update its state variables discontinuously. 104 105 The zero crossings of continuous event trigger functions will be isolated 106 quickly; discontinuous ones have to be "binary chopped" which is more expensive 107 but they will still work. 108 109 We are given a set of weights W for the y's, and a set of tolerances T for the 110 constraint errors. Given an accuracy specification (like 0.1%), the integrators 111 here are expected to solve for y(t) such that the local error 112 |W*yerr|_RMS <= accuracy, and |T*c(t,y)|_RMS <= accuracy at all times. 113 114 TODO: isolation tolerances for witnesses; dealing with simultaneity. 115 **/ 116 class SimTK_SIMMATH_EXPORT Integrator { 117 public: Integrator()118 Integrator() : rep(0) { } 119 ~Integrator(); 120 121 // These are the exceptions that can be thrown by this class. 122 123 class InitializationFailed; 124 class StepSizeTooSmall; 125 class StepFailed; 126 class TriedToAdvancePastFinalTime; 127 class CantAskForEventInfoWhenNoEventTriggered; 128 129 /// Get the name of this integration method 130 const char* getMethodName() const; 131 /// Get the minimum order this Integrator may use 132 int getMethodMinOrder() const; 133 /// Get the maximum order this Integrator may use 134 int getMethodMaxOrder() const; 135 /// Get whether this Integrator provides error control. An error controlled Integrator will dynamically adjust its 136 /// step size to maintain the level of accuracy specified with setAccuracy(). An Integrator which does not provide 137 /// error control cannot do this, and will usually ignore the value specified with setAccuracy(). 138 bool methodHasErrorControl() const; 139 140 /// Supply the integrator with a starting state. This must be called before the first call to stepBy() or stepTo(). 141 /// The specified state is copied into the Integrator's internally maintained current state; subsequent changes to 142 /// the State object passed in here will not affect the integration. 143 void initialize(const State& state); 144 145 /// After an event handler has made a discontinuous change to the Integrator's 146 /// "advanced state", this method must be called to reinitialize the Integrator. 147 /// Event handlers can do varying amounts of damage and some events will require no 148 /// reinitialization, or minimal reinitialization, depending on details 149 /// of the particular integration method. So after the handler has 150 /// mangled our State, we tell the Integrator the lowest Stage which was 151 /// changed and allow the Integrator to figure out how much reinitialization 152 /// to do. 153 /// 154 /// If "shouldTerminate" is passed in true, the Integrator will wrap 155 /// things up and report that the end of the simulation has been reached. 156 void reinitialize(Stage stage, bool shouldTerminate); 157 158 /// Return a State corresponding to the "current" time at the end of the last call to 159 /// stepTo() or stepBy(). This may be an interpolated value earlier than getAdvancedState(). 160 const State& getState() const; 161 /// Get the time of the current State. This is equivalent to calling getState().getTime(). getTime()162 Real getTime() const {return getState().getTime();} 163 164 /// Get whether getState() will return an interpolated state or just the same thing as getAdvancedState() does. 165 /// In most cases, you should not have reason to care whether the state is interpolated or not. 166 bool isStateInterpolated() const; 167 168 /// Return the state representing the trajectory point to which the 169 /// integrator has irreversibly advanced. This may be later than the 170 /// state return by getState(). 171 const State& getAdvancedState() const; 172 /// Get the time of the advanced State. This is equivalent to calling getAdvancedState().getTime(). getAdvancedTime()173 Real getAdvancedTime() const {return getAdvancedState().getTime();} 174 175 /// Get a non-const reference to the advanced state. 176 State& updAdvancedState(); 177 178 /// Get the accuracy which is being used for error control. Usually this is the same value that was 179 /// specified to setAccuracy(). 180 Real getAccuracyInUse() const; 181 /// Get the constraint tolerance which is being used for error control. Usually this is the same value that was 182 /// specified to setConstraintTolerance(). 183 Real getConstraintToleranceInUse() const; 184 185 /// When a step is successful, it will return an indication of what 186 /// caused it to stop where it did. When unsuccessful it will throw 187 /// an exception so you won't see any return value. 188 /// 189 /// When return of control is due ONLY to reaching a report time, 190 /// (status is ReachedReportTime) the integrator's getState() method may return an 191 /// interpolated value at an earlier time than its getAdvancedState() 192 /// method would return. For the other returns, and whenever the report 193 /// time coincides with the end of an internal step, getState() and 194 /// getAdvancedState() will be identical. 195 /// 196 /// Note: we ensure algorithmically that no report time, scheduled time, 197 /// or final time t can occur *within* an event window, that is, we will 198 /// never have t_low < t < t_high for any interesting t. Further, t_report, 199 /// t_scheduled and t_final can coincide with t_high but only t_report can 200 /// be at t_low. The interior of t_low:t_high is a "no man's land" where we 201 /// don't understand the solution, so must be avoided. 202 enum SuccessfulStepStatus { 203 /// stopped only to report; might be interpolated 204 ReachedReportTime =1, 205 /// localized an event; this is the *before* state (interpolated) 206 ReachedEventTrigger =2, 207 /// reached the limit provided in stepTo() (scheduled event) 208 ReachedScheduledEvent=3, 209 /// user requested control whenever an internal step is successful 210 TimeHasAdvanced =4, 211 /// took a lot of internal steps but didn't return control yet 212 ReachedStepLimit =5, 213 /// termination; don't call again 214 EndOfSimulation =6, 215 /// the beginning of a continuous interval: either the start of the 216 /// simulation, or t_high after an event handler has modified the state. 217 StartOfContinuousInterval=7, 218 InvalidSuccessfulStepStatus = -1 219 }; 220 /// Get a human readable description of the reason a step returned. 221 static String getSuccessfulStepStatusString(SuccessfulStepStatus); 222 223 /// Integrate the System until something happens which requires outside processing, and return a status code describing what happened. 224 /// @param reportTime the time of the next scheduled report 225 /// @param scheduledEventTime the time of the next scheduled event 226 SuccessfulStepStatus stepTo(Real reportTime, Real scheduledEventTime=Infinity); 227 /// Integrate the System until something happens which requires outside processing, and return a status code describing what happened. 228 /// @param interval the interval from the current time (as returned by getTime()) until the next scheduled report 229 /// @param scheduledEventTime the time of the next scheduled event 230 SuccessfulStepStatus stepBy(Real interval, Real scheduledEventTime=Infinity); 231 232 233 /// Get the window (tLow, tHigh] within which one or more events have been localized. This may only be called 234 /// when stepTo() or stepBy() has returned ReachedEventTrigger. 235 Vec2 getEventWindow() const; 236 /// Get the IDs of all events which have been localized within the event window. This may only be called 237 /// when stepTo() or stepBy() has returned ReachedEventTrigger. 238 const Array_<EventId>& getTriggeredEvents() const; 239 /// Get the estimated times of all events which have been localized within the event window. This may only be called 240 /// when stepTo() or stepBy() has returned ReachedEventTrigger. 241 const Array_<Real>& getEstimatedEventTimes() const; 242 /// Get EventTriggers describing the events which have been localized within the event window. This may only be called 243 /// when stepTo() or stepBy() has returned ReachedEventTrigger. 244 const Array_<Event::Trigger>& getEventTransitionsSeen() const; 245 246 247 // TERMINATION // 248 249 /// Once the simulation has ended, getTerminationReason() may be called to find out what caused it to end. 250 enum TerminationReason { 251 /// The simulation reached the time specified by setFinalTime(). 252 ReachedFinalTime = 1, 253 /// An error occurred which the Integrator was unable to handle. 254 AnUnrecoverableErrorOccurred = 2, 255 /// An event handler function requested that the simulation terminate immediately. 256 EventHandlerRequestedTermination = 3, 257 /// This will be returned if getTerminationReason() is called before the simulation has ended. 258 InvalidTerminationReason = -1 259 }; 260 261 /// Get whether the simulation has terminated. If this returns true, you should not 262 /// call stepTo() or stepBy() again. 263 bool isSimulationOver() const; 264 265 /// Get the reason the simulation terminated. This should only be called if 266 /// isSimulationOver() returns true. 267 TerminationReason getTerminationReason() const; 268 269 /// Get a human readable description of the termination reason. 270 static String getTerminationReasonString(TerminationReason); 271 272 /// Reset all statistics to zero. 273 void resetAllStatistics(); 274 275 /// Get the size of the first successful step after the last initialize() call. 276 Real getActualInitialStepSizeTaken() const; 277 278 /// Get the size of the most recent successful step. 279 Real getPreviousStepSizeTaken() const; 280 281 /// Get the step size that will be attempted first on the next call to stepTo() or stepBy(). 282 Real getPredictedNextStepSize() const; 283 284 /// Get the total number of steps that have been attempted (successfully or unsuccessfully) since 285 /// the last call to resetAllStatistics(). 286 int getNumStepsAttempted() const; 287 /// Get the total number of steps that have been successfully taken since the last call to resetAllStatistics(). 288 int getNumStepsTaken() const; 289 /// Get the total number of state realizations that have been performed since the last call to resetAllStatistics(). 290 int getNumRealizations() const; 291 /// Get the total number of times a state positions Q have been projected 292 /// since the last call to resetAllStatistics(). 293 int getNumQProjections() const; 294 /// Get the total number of times a state velocities U have been projected 295 /// since the last call to resetAllStatistics(). 296 int getNumUProjections() const; 297 /// Get the total number of times a state has been projected (counting 298 /// both Q and U projections) since the last call to resetAllStatistics(). 299 int getNumProjections() const; 300 /// Get the number of attempted steps that have failed due to the error being unacceptably high since 301 /// the last call to resetAllStatistics(). 302 int getNumErrorTestFailures() const; 303 /// Get the number of attempted steps that failed due to non-convergence of 304 /// internal step iterations. This is most common with iterative methods 305 /// but can occur if for some reason a step can't be completed. It is reset 306 /// to zero by resetAllStatistics. 307 int getNumConvergenceTestFailures() const; 308 /// Get the number of attempted steps that have failed due to an error when realizing the state since 309 /// the last call to resetAllStatistics(). 310 int getNumRealizationFailures() const; 311 /// Get the number of attempted steps that have failed due to an error 312 /// when projecting the state positions (Q) since the last call to 313 /// resetAllStatistics(). 314 int getNumQProjectionFailures() const; 315 /// Get the number of attempted steps that have failed due to an error 316 /// when projecting the state velocities (U) since the last call to 317 /// resetAllStatistics(). 318 int getNumUProjectionFailures() const; 319 /// Get the number of attempted steps that have failed due to an error 320 /// when projecting the state (either a Q- or U-projection) since 321 /// the last call to resetAllStatistics(). 322 int getNumProjectionFailures() const; 323 /// For iterative methods, get the number of internal step iterations in steps that led to 324 /// convergence (not necessarily successful steps). Reset to zero by resetAllStatistics(). 325 int getNumConvergentIterations() const; 326 /// For iterative methods, get the number of internal step iterations in steps that did not lead 327 /// to convergence. Reset to zero by resetAllStatistics(). 328 int getNumDivergentIterations() const; 329 /// For iterative methods, this is the total number of internal step iterations taken regardless 330 /// of whether those iterations led to convergence or to successful steps. This is the sum of 331 /// the number of convergent and divergent iterations which are available separately. 332 int getNumIterations() const; 333 334 /// Set the time at which the simulation should end. The default is infinity. Some integrators may 335 /// not support this option. 336 void setFinalTime(Real tFinal); 337 /// Set the initial step size that should be attempted. The default depends on the integration method. 338 /// Some integrators may not support this option. 339 void setInitialStepSize(Real hinit); 340 /// Set the minimum step size that should ever be used. The default depends on the integration method. 341 /// Some integrators may not support this option. 342 void setMinimumStepSize(Real hmin); 343 /// Set the maximum step size that should ever be used. The default depends on the integration method. 344 /// Some integrators may not support this option. 345 void setMaximumStepSize(Real hmax); 346 347 /// Set the integrator to use a single fixed step size for all steps. This 348 /// is exactly equivalent to calling setInitialStepSize(), 349 /// setMinimumStepSize(), and setMaximumStepSize(), passing the same value 350 /// to each one. This will therefore not work correctly if the integrator 351 /// does not support minimum and/or maximum step sizes. 352 void setFixedStepSize(Real stepSize); 353 354 /// Set the overall accuracy that should be used for integration. If the 355 /// Integrator does not support error control (methodHasErrorControl() 356 /// returns false), this may have no effect. 357 void setAccuracy(Real accuracy); 358 /// Set the tolerance within which constraints must be satisfied. 359 void setConstraintTolerance(Real consTol); 360 /// (Advanced) Use infinity norm (maximum absolute value) instead of 361 /// default RMS norm to evaluate whether accuracy has been achieved for 362 /// states and for constraint tolerance. The infinity norm is more strict 363 /// but may permit use of a looser accuracy request. 364 void setUseInfinityNorm(bool useInfinityNorm); 365 /// (Advanced) Are we currently using the infinity norm? 366 bool isInfinityNormInUse() const; 367 368 369 /// Set the maximum number of steps that may be taken within a single call 370 /// to stepTo() or stepBy(). If this many internal steps occur before 371 /// reaching the report time, it will return control with a ReachedStepLimit 372 /// status. If nSteps <= 0, the number of steps will be unlimited. 373 void setInternalStepLimit(int nSteps); 374 375 /// Set whether the Integrator should return from stepTo() or stepBy() after 376 /// every internal step, even if no event has occurred and the report time 377 /// has not been reached. The default is false. 378 void setReturnEveryInternalStep(bool shouldReturn); 379 380 /// Set whether the system should be projected back to the constraint 381 /// manifold after every step. If this is false, projection will only be 382 /// performed when the constraint error exceeds the allowed tolerance. 383 void setProjectEveryStep(bool forceProject); 384 /// Set whether the Integrator is permitted to return interpolated states 385 /// for reporting purposes which may be less accurate than the "real" 386 /// states that form the trajectory. Setting this to false may 387 /// significantly affect performance, since the Integrator will be forced 388 /// to decrease its step size at every scheduled reporting time. 389 /// 390 /// This option is generally only meaningful if interpolated states are 391 /// less accurate than other states on the trajectory. If an Integrator can 392 /// produce interpolated states that have the same accuracy as the rest of 393 /// the trajectory, it may ignore this option. 394 void setAllowInterpolation(bool shouldInterpolate); 395 /// Set whether interpolated states should be projected back to the 396 /// constraint manifold after interpolation is performed. The default is 397 /// "true". 398 void setProjectInterpolatedStates(bool shouldProject); 399 /// (Advanced) Constraint projection may use an out-of-date iteration 400 /// matrix for efficiency. You can force strict use of a current iteration 401 /// matrix recomputed at each iteration if you want. 402 void setForceFullNewton(bool forceFullNewton); 403 404 /// OBSOLETE: use getSuccessfulStepStatusString(). successfulStepStatusString(SuccessfulStepStatus stat)405 static String successfulStepStatusString(SuccessfulStepStatus stat) 406 { return getSuccessfulStepStatusString(stat); } 407 408 protected: getRep()409 const IntegratorRep& getRep() const {assert(rep); return *rep;} updRep()410 IntegratorRep& updRep() {assert(rep); return *rep;} 411 412 // opaque implementation for binary compatibility 413 IntegratorRep* rep; 414 friend class IntegratorRep; 415 }; 416 417 } // namespace SimTK 418 419 #endif // SimTK_SIMMATH_INTEGRATOR_H_ 420