Monte Carlos Method: Particular Equation Style Required

**johnact121** · July 29th 2011, 01:48 AM

I require a numerical program to generate a random equation for me with the requirements listed below. This curve is for a velocity profile of an object, must have positive velocity at all times and has limits on the maximum acceleration and deceleration. I desire to create a slightly different curve each time I generate one, still meeting all the requirements.

The requirements for the function, v(t) for 0 < t < t_max, is:

global max < b
global min = 0
v(0) = v(t_max) = 0
v ' (t)_max <= a
v ' (t)_min >= c
integral ( v(t) ) = d

The following two requirements, as seen on the image, add no new information to the problem:

v ' (0) <= a
v ' (t_max) >= c

An image with this information is shown below:

Monte Carlos Method: Particular Equation Style Required-mote_carlo_req.png

Any help on this problem would be appreciated. If it is very lengthy to explain how to do kindly point me in the right direction as I have some maths textbooks.

For the numerical program I plan to use either MATLAB or Mathematica.

Thanks again!

**Ackbeet** · July 29th 2011, 05:57 AM

An interesting problem. I have a few comments:

1. $\int_{0}^{t_{\text{max}}}v(t)\,dt = d$ is superfluous if $\min_{0\le t\le t_{\text{max}}}v(t)=0.$ Unless, as I suspect is the case, $d$ is a number you know ahead-of-time. Certainly, if $v(t)\ge 0$ on the interval in question, then $\int_{0}^{t_{\text{max}}}v(t)\,dt$ is the distance traveled, whatever that is.

2. Depending on the maximum slopes of the velocity, a and c, and the required distance traveled, d, you may not be able to produce a function satisfying all those requirements. In particular, if the required distance traveled, d, is large compared with a and c, you will not be able to exhibit the required function.

3. My solution: From the origin, start rising at the maximum slope a. Backwards from the endpoint, descend at the maximum slope c. Then fiddle with what's in-between those two lines to make the area equal to d. You could try to make the area under the extreme lines equal to half of d, and account for the other half of d by various shenanigans that allow you to change things and still keep the same area.

**SpringFan25** · July 29th 2011, 08:44 AM

Notice that if you have n functions which satisfy the constraints, then any normalised linear combinations of them will also match the constraints

ie, if you have a collection of n functions $f_i(t)$ which satisfy the contraints, then the following does as well:
$g(t) = \sum_{i=1}^{n} w_i f_i(t)$

with the resitrictions $0 \leq w_i \leq 1$ and $\sum w_i = 1$

So, you can generate your fluctuations by randomising the weights $w_i$ . My suggested algorithm for you would be:
(1) Manually find a few (say, 10) functions which satisfy your contraints using any method (eg, trial and error). This should be feasible enough, and youll only have to do it once.
(2) use a random number generator to assign a weight to each function you found in (1)
(3) then plot the function g(t) corresponding to those randomly chosen weights

to create more realisations, you can generate a new set of random numbers for step (2).

.

Notes

To Generate random weights you can generate $n$ random numbers $X_i \geq 0$ by any method (all maths programs shold be able to do this), then use $w_i=\frac{X_i}{\sum_{i=1}^{n} X_i}$
a particularly simple implementation is with n=2 since you could then use a single random number generator to get $w_1$ betgween 0 and 1, then use $w_2 = 1- w_1$
Obviously the n functions you find will need to be dsignificantly different if you want to see significant fluctions
the types of paths you observe will be heavily dependent on the shapes of the functions you choose, which may be unrealistic for your purpose. i would try to choose functions that look different to each other.

PS: Although i started this post by saying "notice that", you should check the result in case i made a mistake.

**johnact121** · July 29th 2011, 06:28 PM

Thank you Ackbeet and SpringFan25 for your prompt replies.

Ackbeet:

Unless, as I suspect is the case, $d$ is a number you know ahead-of-time.

Yes, d is known and the values a, b, c are also known. Yes, since this is a velocity profile the derivative is acceleration and the integration is the displacement of the object.

For example, let us choose the values of:

distance (integral) = d = 74m
max accn = a = 2m/s^2
max dccn = c = -4m/s^2
t_max = 20 seconds.

One function which achieves this is a max acceleration for 2 seconds, holding 4m/s for 17 seconds, and max deceleration for 1 second. This generates a function in a trapezium shape with no negative values, meeting all the requirements.

I wish, however, to find a function (continuous) which meets these conditions (notice I have removed the constraint for the global maximum (b) as this is taken care of by the maximum rates of change).

SpringFan25:

Notice that if you have n functions which satisfy the constraints, then any normalised linear combinations of them will also match the constraints

Yes, your idea of normalised linear combinations is a good idea. My first question would be: how to I generate the continuous functions in the first place?

I have the following ideas so far:

A numerical solver to solve y = a.x^(n) + b.x^(n-1) + ... &c. with the various constrains put into the numerical program.
Analytical method to solve for the equation.
Simple constant acceleration equations for trapezium shaped functions (not continuous but meet the conditions)
Others?

**SpringFan25** · July 30th 2011, 03:18 AM

Simple constant acceleration equations for trapezium shaped functions (not continuous but meet the conditions)

Minor point: I think that trapezium functions are continuous but not differentiable. I agree that they meet the conditions as written in post #1.

Major point: if you are going to allow non-differentiable functions then the problem is dramatically simplified. You can just start/finish your graph with two short lines with slope a and c respectively. For simplicity, choose the length of those lines so that the area underneath them is 1.

Now all you need to do is choose some random points in the middle and connect them with straight lines. Scale up/down these middle points by a constant factor so that the area under them is d-1. This will bring the total area under the graph to $d$ as required.

If re-scaling causes the global maximum constraint to be breached, discard your middle points points and randomly choose some more. Alternatively reduce the highest point and then repeat the scaling factor calculation.

Edit: I think ackbeet already suggested the above in post #2.

Many combinations will be discarded but a computer could make thousands/millions of attempts in a very short time.

Poosible Smooth functions
Notice that there is no constraint on v' except at t_max and t_min. Trigonometric functions (specifically, combinations of sine waves) look promising if you run out of simpler options (polynomials etc):

they can be combined to make many different shapes.
Because they are periodic they can be easily parameterized to meet the slope restrictions at t_min and t_max.
An upper bound on the global maximum is also easy to compute. (the global maximum cant be more than the sum of the amplitudes of your waves)
The integrals should also be tractable as long as you confine yourself to simple functions, eg combinations of sine waves with different periods/phases/amplitudes. The integral calculations will be slightly easier if you chose half-periods that fit into the period (t_max - t_min) a whole number of times.
One minor challenge will be ensuring that the global minimum is not breached but this can be overcome by trial and error.

PS: I expect the above is tractable but i haven't actually tried it myself! I wont look at it further unless you actually need smooth functions.

**SpringFan25** · July 30th 2011, 08:56 AM

it was too much fun to ignore so I've attached an excel implementation of the smooth trig functions approach i suggested above.

The main changes from my above post are

t_min is assumed 0
i was wrong to say that there was no restriction on v'() except at t_min and t_max. i had misread your first post.
All sine waves complete an odd number of half periods in the time interval. This makes all graphs symettric so v(0) = v(tmax) =0. We also have v'(0) = -v(t_max).
To simplify calculations all sine waves have a phase shift of 0 or $\pi$
Ive used a stronger gradient constraint than you need. I used |v'(t)| < min(|a|,|c|). This is useful because an upper bound on |v'(t)| can be easily calculated from the amplitude and period of the sine waves (derivation below)
There is a big "calculate me" button to generate a new random graph. it requires macros enabled to run. I dont accept responsibility if this breaks your computer/life etc.

Upper bound on |v'(t)|
I'll derive this. Suppose v(t) is the sum of $n$ sine waves with ammplitude $A_i$ , period $B_i$ , and phase-shift $C_i$ as follows:

$f_i(t) = A_i \sin(\frac{2 \pi}{B_i} t + C_i)$
$v(t) = \sum f_i(t)$

Differentiate the first equation. This gives:
$f'_i(t) = \frac{2 \pi A_i}{B_i} \cos (2 \pi B_i t + C_i)$

Since $|cos(x)| \leq 1$ it follows that
$|f'_i(t)| \leq \frac{2 \pi A_i}{B_i}$

And we can sum as many sine waves as we want, adding up the upper bounds for each. This is handy since v(t) is the sum of the sine waves:

$|v'(t)| \leq \sum |f'_i(t)| \leq 2 \pi \sum \frac{A_i}{B_i}$
(the left hand inequality is valid. Think "triangle inequality" with lots of variables.)

How does this help?
Note that $|v'(t)| \leq |a| \implies v'(t) \leq a ~~ \text{ since } a > 0$
and, $|v'(t)| \leq |c| \implies v'(t) \geq c ~~ \text{ since } c < 0$

Combine the above pair of results:

$|v'(t)| \leq |a| \text{ and } |v'(t)| \leq |c| \implies c \leq v'(t) \leq a$

it follows that
$|v'(t)| \leq min(|a|,|b|) \implies c \leq v'(t) \leq a$
So testing that $|v'(t)| \leq min(|a|,|c|)$ will ensure that the acceleration constraints are met.

Instead of evaluating |v'(t)| everywhere, It's sufficient to test the upper bound of |v'(t)| that was derived in the previous part of this post. ie the test is:

$2 \pi \sum \frac{A_i}{B_i} \leq \min{(|a|,|c|)}$

Meeting all the constraints
for the rest of this post assume that v(t) is constructed as follows:

$v(t) = \sum f_i(t)$
$f_i(t) = A_i \sin(\frac{2 \pi}{B_i} t + C_i) ~~~ \text{ with }A_i,B_i \geq 0$

v(0) = v(tmax) = 0
This is the easiest one. Consider time 0:
$f_i(0) = A_i \sin(\frac{2 \pi}{B_i} \cdot 0 + C_i)$
$f_i(0) = A_i \sin(C_i)$

so a sufficient condition for $f_i(0)$ is $\sin(C_i)=0$

so the phase shift should be $0, \pi, 2\pi, 3\pi$ . Since the period is $2 \pi$ anyway, so only unique phase shifts are 0 and $\pi$ .

What about $f_i (t_{max})$ ?
The sin function returns to its starting point every half-period. So its suffiient to ensure that the interval $[0,t_{max}]$ has a whole number of half periods in it. Then we can use $f_i(t_{max}) = f_i(0) = 0$ .

(in fact, we will want an odd number of half-periods later on.).

c < v'(t) < a
I already showed that a sufficient condition is
$2 \pi \sum \frac{A_i}{B_i} \leq \min{(|a|,|c|)}$

Total displacement = d
This is the tricky one. It's very convenient to restrict the period of each sine wave so that an odd number of half periods fit in the interval $[0,t_{max}]$ . This means that the number of completed oscilations will always be 1.5, 2.5, 3.5, 4.5 etc. This is important because the integral of a full oscillation is always zero.

So we only need to compute the integral of the remaining half oscillation of each wave. We already established the the phase shift $C_i$ is 0 or $\pi$ . For the case where the phase shift is 0, This area is equal to:

$A_i \int_0^{\pi} \sin(B_i t + 0) dt = A_i \left[ \frac{-\cos(B_i t)}{B_i}\right]_0^{\pi}$

For the case where the phase shift is $\pi$ the half-period area is equal to
$A_i \int_{0}^{\pi} \sin(B_i t + \pi) dt$

But remember that the sin function has this property: $sin(x + \pi) = -sin(x)$

So for the phase shifted case the area is just -1 times the area we just calculated:
$A_i \int_{0}^{\pi} \sin(B_i t + \pi) = -A_i \int_{0}^{pi} \sin(B_i t) dt = A_i \left[ \frac{\cos(B_i t)}{B_i}\right]_0^{\pi}$

Finally, a single formula covering both cases:

$\int_{0}^{\pi} f_i dt = A_i \left[ \frac{-\cos(B_i t)}{B_i}\right]_0^{\pi} \cdot (-1)^{\frac{C_i}{\pi}}$
(if you dont udnerstand what the -1 factor is doing at the end, it evaluates to 0 if there is no phase shift, and -1 if there is a phase shift)

Hence the area under the graph $\Omega$ is:
$\Omega = \int_0^{t_{max}} v(t) dt = \sum_i \left( \int_0^{t_{max}} f_i(t) dt \right)= \sum_i \left( A_i \left[ \frac{-\cos(B_i t)}{B_i}\right]_0^{\pi} \cdot (-1)^{\frac{C_i}{\pi}} \right)$

We require the area under the graph to be exactly d. This is easy (compared to all that integration!). Simply set new coefficients:

$\tilde{A}_i = A_i \frac{d}{\Omega}$

This scales up the entire function so that the total area will be d.

Obviously, this scale up/down may brake the acceleration constraints, since the slopes will also be scaled. The attached workbook checks this, and tries a new set of $A_i, B_i$ if there is a breach.

Algorithm
1) Randomly choose $A_i, B_i.C_i$ subject to the rules derived above, ie:

$C_i$ = 0 or $\pi$
$B_i = \frac{t_{max}}{k+0.5} ~~~ : ~ k \in \mathbb{N}$

2) Calcualte adjusted $\tilde{A}_i = A_i \frac{d}{\Omega}$

3) Test the adjusted coefficients meet the accelleration constraints. If they dont, discard and return to step 1. The test is $2 \pi \sum \frac{A_i}{B_i} \leq \min{(|a|,|c|)}$

The attached workbook does the above. I've chosen ranges on $A_i,B_i$ that produce reasonably varied output.

Sample output
here is some sample output

**johnact121** · August 2nd 2011, 11:41 PM

SpringFan25,

Thank you very much for your detailed reply. I will need some time to go through it and understand it.

One thing to point out, it seems that all the functions generated are even and thus symmetrical. That constraint is not required - is there a way to get around that?

**SpringFan25** · August 3rd 2011, 03:45 AM

The symmetry is convenient because it makes sure that the velocity graph returns to its starting value, ie v(0) = v(tmax) =0.

You can get rid of the symmetry by adding a further asymmetric function to v(t) which also satisfies v(0) = t(tmax) = 0.

One way is to add another pair of sine waves which are 180 degrees out of phase, and with one having double the frequency of the other:

examples of 2 sine waves and their sum are below to illustrate:
plot y=sin(x),y=0.5sin(2x-pi) on 0,pi - Wolfram|Alpha
plot y=sin(x) + 0.5sin(2x-pi) on 0,pi - Wolfram|Alpha

However these functions are not ideal as the low gradient at the start means that small fluctuations from any other waves you put on top can easily make the overall value of v(t) negative at low values of t.

Thread: Monte Carlos Method: Particular Equation Style Required

LinkBack

Thread Tools

Display

Monte Carlos Method: Particular Equation Style Required

Re: Monte Carlos Method: Particular Equation Style Required

Re: Monte Carlos Method: Particular Equation Style Required

Re: Monte Carlos Method: Particular Equation Style Required

Re: Monte Carlos Method: Particular Equation Style Required

Re: Monte Carlos Method: Particular Equation Style Required

Re: Monte Carlos Method: Particular Equation Style Required

Re: Monte Carlos Method: Particular Equation Style Required

Similar Math Help Forum Discussions

Help required? Solve the Differential equation by the Laplace transform method"

Simpson's method and Monte Carlo Method

Monte Carlo method

monte carlo method

Monte Carlo Method

Search Tags