0
\$\begingroup\$

I'm simulating a simple cart-pendulum system for purposes of testing different control strategies and data-driven approaches. I need to simulate the dynamical system one time-step at a time.

I'm trying to make the simulation as realistic as possible (non-linear, with disturbances and noise etc.)

At the moment I'm just using the Euler method:

y_dot = cartpend_dydt(t, y)

# Simple state update (Euler method)
y += tau*y_dot

where y_dot is the time-derivative of y, the state-vector (length 4) and tau is the time-step.

Since it's a very fast calculation I can run at 50 steps/second and so I'm sure Euler method is accurate enough.

But what if I want to reduce the sample efficiency to say 10 steps/second (more typical real-world scenario)?

Then I assume I should do something more accurate (and a bit more computationally expensive) each time step. But what's 'the next step up' from Euler?

I tried using a numerical integrator (in Python):

from scipy.integrate import solve_ivp

sol = solve_ivp(cartpend_dydt, [t, t + tau], y)
y = sol.y[:, 1]

where cartpend_dydt(t, y) is a function that returns the time derivative of y.

I think scipy.integrate.solve_ivp is an "explicit Runge-Kutta method of order 5(4)" by default.

But is this over-kill for a relatively short time-step? Is there a simpler method that would be better than Euler but not too computationally intensive?

\$\endgroup\$
3
  • \$\begingroup\$ There's a lot of off-topic topics for this question. First of all, Python scripting and their libraries... The only on-topic questions for Python is either circuit simulation or embedded scripting. Second, a cart-pendulum is a control system, yes, but not an electrical control system. Third, you're asking for a more accurate mathematical algorithm than Euler method. \$\endgroup\$ – user103380 May 16 '19 at 20:14
  • \$\begingroup\$ Yes, where to post control engineering questions is an open question. Thanks for clarifying that I'm asking for a more accurate mathematical algorithm than Euler method. Are you implying I should be posting on Mathematics? I'm not clear what you are recommending. \$\endgroup\$ – Bill May 16 '19 at 20:19
  • \$\begingroup\$ I'm not entirely sure where this question would belong since you're talking about mathematical algorithms and accuracy optimization (implying both Math SE and Stackoverflow respectively). However, if you did this by hand, you could ask this on Math SE. Once you get a more accurate algorithm, then you can apply it to Python. \$\endgroup\$ – user103380 May 16 '19 at 20:26
3
\$\begingroup\$

The next step is the trapezoidal rule; i.e., approximate the integrand with line segments between the sample points. If you draw it, you see that the correction is graphically the area of a rectangular triangle. Numerically one integration step gives an increment (delta t) x (the average of consequent samples).

The next step is to use Simpson's rule.

Numerical integration has been a must in science and engineering for hundreds of years. Before computers using effective methods was more important than today. Then the mathematicians also tried to estimate the error limit for the result. Their theories are still valid although lousy programming today only decreases the time increment until the result seems to stabilize.

\$\endgroup\$
2
  • \$\begingroup\$ Thanks, this is very helpful. Based on your description of the trapezoidal rule, is it the same as multiplying the time step-size by the average of the derivatives of y at t and t-1? \$\endgroup\$ – Bill May 17 '19 at 0:02
  • 1
    \$\begingroup\$ No, I mean you should calculate at first next point Y(k+1) with Euler's method from point Y(k) ie. Y(k+1)=Y(k)+(Ydot(k))dT. Then calculate the next derivative estimate Ydot(k+1). Finally recalculate Y(k+1) with the average 0.5(Ydot(k)+Ydot(k+1)). Strictly this is not the trapezoidal rule because Ydot (k+1) is not known explicitly, \$\endgroup\$ – user287001 May 17 '19 at 8:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.