# Make a signal that fits another the best possible with a limitation in the 2nd derivative

Consider this step function:

The signal that "fits" this should look like the following (in green):

The corners are now smooth because the maximum second derivative allowed is not infinite anymore.

The idea is to find the most close to the original signal given the maximum second derivative allowed, minimizing the difference between the two signals (in orange):

If your system can predict things before it happens it could do something like this:

Minimizing the total difference even more. But also could be achieved if it works in non real time or has some kind of look ahead. But that will be not always possible or desirable. This is not in the very core of my question, however is also a interesting option.

So, my question is: How can this be done? Is it possible (eg: exist more than one curve that minimizes error to the same value)?

Also to me looks like a filter, is there any filter which definition is something like this (hopefully a digital one)?

[Note: This is a second answer, split from the original, by request of the OP. See the comments on my original answer.]

The idea behind the parabolic solutions being the optimal solution is the specification was a limitation on the second derivative. Thus the second derivative should be as maximized as possible, so with the limitation, it means it will be constant for the durations, leading to parabolas. First in a maximum acceleration manner, and then in a maximum decelaration manner.

The model is that of a teenage driver: Full throttle or full brake if possible.

The next point is considered the target. If the target is reachable in one point, the "velocity" is set appropriately. In the case of an overshoot, or the target moving and the teenager can't stop in time, full brake is applied. If it is possible to go faster without having to brake yet, full throttle is applied.

The code is fairly simple. The formula for distance, being discrete not continuous, is

$$Distance = Thrust \cdot Steps(Steps+1) / 2$$

From this, the number of steps (theLimitTime) can be calculated using the quadratic formula.

Here is the plot. The path follows the optimal parabolas.

import matplotlib.pyplot as plt
import numpy             as np

#=============================================================================
def main():

#---- Set Parameters

a = .01  # Maximum Allowable Acceleration
N = 150

#---- Set up the Step Function

x = np.zeros( N )

x[30:N] = 1.0

#---- Allocate and Start "Filtered" Signal

y = np.zeros( N )

y[0] = x[0]
y[1] = x[1]

#---- Calculate "Filtered" Signal

for n in range( 1, N ):
v = CalculateVelocity( a, x[n], y[n-1], y[n-1]-y[n-2] )
y[n] = y[n-1] + v

#---- Show Result at Step Region

for n in range( 25, 60 ):
zeroth = "%8.4f" % ( y[n] )
first  = "%8.4f" % ( y[n] - y[n-1] )
second = "%8.4f" % ( y[n+1] - 2.0 * y[n] + y[n-1] )
print n, zeroth, first, second

#---- Plot the Results

plt.plot(x)
plt.plot(y)
plt.legend(['x[n]', 'y[n]'])
plt.show()

#=============================================================================
def CalculateVelocity( argMaxThrust, argTarget, argSpot, argVelocity ):

#---- Orient to a Positive Reference

if( argVelocity < 0 ):
theDirection      = -1.0
theSpeed          = -argVelocity
theTargetDistance = argSpot - argTarget
else:
theDirection      = 1.0
theSpeed          = argVelocity
theTargetDistance = argTarget - argSpot

#---- Bail on Full Reverse Thrust Needed

theMinSpeed = theSpeed - argMaxThrust

if( theTargetDistance < theMinSpeed ):
return theMinSpeed * theDirection

#---- Bail if Arrival Possible

theMaxSpeed = theSpeed + argMaxThrust

if( theTargetDistance < argMaxThrust ):
if( theSpeed < argMaxThrust ):
return theTargetDistance * theDirection

#---- Calculate Stop Time for Distance

theTemp = 2.0 * theTargetDistance / argMaxThrust

theLimitTime = np.sqrt( 0.25 + theTemp ) - 0.5

#---- Calculate Right Speed

theLimitSpeed = theLimitTime * argMaxThrust

#---- Exit

if( theLimitSpeed > theMaxSpeed ):
return theMaxSpeed * theDirection

if( theLimitSpeed < theMinSpeed ):
return theMinSpeed * theDirection

return theLimitSpeed * theDirection

#=============================================================================
main()

Hmmmmmmmmm, interesting question.

Since you want to use the second derivative as your criteria, it would seem that you would want to have the maximum second derivative absolutie value for as short of a duration as possible. This would suggest piecing together parabolas, matching the first derivatives at the joints.

How to do this algorithmically will take a little more thought. Somehow you have to identify regions of "violations" and then figure out how big of an interval you have to use to patch it piecewise with parabolas.

Just some thoughts. If I think of a sure fire way, I'll add it on.

Ced

Followup:

It turns out to be really simple to make a FIR that reacts to a step function with a upward parabola for the first half and a downward parabola for the second half. It just becomes a matter of selecting the width and the rise rate to set the value of the second derivative. I haven't done the math yet for the actual value, but here is a program which gives proof of concept.

h[k] = [ 1  2  3  4  3  2  1 ]  / 16

It appears the second derivative proxy $y[n-1] - 2 * y[n] + y[n+1]$ has an absolute value of the normalizing factor of the integer step vector, e.g. $1/16$.

import numpy as np

#======================================================================
def main():

P = 2

#---- Build the Weights

K = 2 * P - 1

h = np.zeros( K )

for k in range( P ):
h[k] = k + 1

for k in range( P - 1 ):
h[P+k] = P - k - 1

theRecip = 1.0 / sum( h )

print h, sum( h )

h *= theRecip

print h

#---- Set up the Step Function

x = np.zeros( 20 )

x[10:20] = 1.0

print x

#---- Apply Filter

y = np.zeros( 20 )

for n in range( 8, 20 ):
Sum = 0.0
for k in range( K ):
Sum += x[n-k] * h[k]
y[n] = Sum

#---- Show Result

for n in range( 6, 18 ):
second = y[n-1] - 2.0 * y[n] + y[n+1]
print n, y[n], second

#======================================================================
main()

Yet a better answer:

[Moved to a separate answer]

• Wow! Interesting aproach – Juan Molina Apr 1 '18 at 2:52
• @Juan Molina, Thanks. I'm pretty sure that is the only way to get a minimum duration alteration. If it is just a step function, the solution is easy. Obviously the parabolas are going to be alternating upward and downward. To match the joints, they should connect halfway between the extremes of the two parabolas. But matching entry and exit conditions seems a bit tricky. It would seem that internal parabolas would alway contain the extreme, but the edge ones may be partial parabolas. I am still thinking about it. Probably best to sleep on it. – Cedron Dawg Apr 1 '18 at 3:02
• Yes, this also apply if instead of a simple step, it suddenly canges the target in a stepped way, and if you think a discrete time signal as all litle steps you solved for every discrete time signal :) but I want to know if there are also some simpler algorithms for this like some sort of iir filter, PID controller or so. If exist, for sure the parabolas can lead you there. Btw I'm still playing with the parabolas in a piece of paper and still amazed how obious the solution was – Juan Molina Apr 1 '18 at 3:19
• If the system that smothes with the parabolas is linear, it has a well defined impulse response thus can be applied as a filter. Now I'm not totally sure if it is, but it looks like – Juan Molina Apr 1 '18 at 3:27
• @Juan Molina, I'm happy I could help. About your comment about disrete singnals. I wouldn't consider it a step unless you had two horizontal points going in and two horizontal points going out. I think this is the optimum theoretical solution given the constraints you set, but I'd also point out that you will have a step function in the second derivative at the joints, making a "real life" solution likely unattainable. – Cedron Dawg Apr 1 '18 at 3:29

You can make a discretized regularized linear equation system.

• $\bf d$ is the original signal
• $\bf v$ is what we add to the signal, (the additive change)
• $\bf v+d$ is the result.

So there are two things we want:

1. We want the result to have as small second derivative as possible.
2. We want the change to be as small as possible. $$\min_{\bf v}\left\{\|{\bf D_2}({\bf v+d})\|_2^2 + \|{\bf v}\|_2^2\right\}$$

Now what you need to do is to solve this. You can use your favourite method!

It will be slightly more challenging to impose the causality requirement, we can leave it as an exercise for the curious student.

EDIT turns out we need an extra term in excess to the one above:

$$+\epsilon\|{\bf Dv}\|_2^2$$

the reason being that the second order regularizer above only limits the absolute value of second derivative and does not prohibit any ringing effects.

As we can see, if $\epsilon = 0$ we get ringing/overshoot, but at a rather small $\epsilon = 0.006$ the ringing disappears.

Edit: For the causality we can add a diagonal weight matrix $\bf C$ working on prohibiting the change selectively for $t\leq 0$:

$$+\epsilon_C\|{\bf Cv}\|_2^2$$

• As an alternative, one could minimize the error under the constraint, that the second derivative of the result should be less than the given maximum value. – Uroc327 Apr 1 '18 at 11:40
• @Uroc327: Yes you are right, but I think it is a more difficult constraint. – mathreadler Apr 1 '18 at 12:14
• @mathreadler sounds like you'd approach this with Lagrange multiplication, then; in your vector notation, that'd be relatively easy to formulate! – Marcus Müller Apr 1 '18 at 13:58
• @Juan Molina, mathreader. Check out my followup for a much simpler solution. – Cedron Dawg Apr 1 '18 at 21:48
• Beautifull plot!. Does it mean that the least error possible is a curve with overshoot? I considered it as a possibility but never thought that was the real solution. Could you describe the process that leaded to the solution? Thanks!. Also I worked something for discrete time and no time prediction, but it's still incomplete. I commented the jupyter notebook link in @Cedron Dawg's answer also here. – Juan Molina Apr 2 '18 at 3:27

I think Ced is correct here; parabolas are going to get you the most bang for your buck where minimizing the second derivative is concerned.

This answer gives a slightly more generalized solution for achieving that parabolic shape, and other shapes (in case your requirements change a bit, maybe). I'm going to make a few assumptions based on the image in your question:

• Your input is something like a step function; that is, the same value appears several times in a row before a different value appears.

• Your maximum allowed second derivative is large enough that your output signal can usually "catch up" to the input signal before that value changes.

• You don't want any overshoot (ringing/clipping).

## Parabolic shape

First, let's express the parabolic shape as a simple equation. We'll have it begin at $(-1,-1)$ and end at $(1,1)$; in other words we'll have it cover two distance units (Y) over two time units (X).

$$y = x\left(2-\left|x\right|\right)$$

It looks like this, along with its first and second derivatives. From now on I'll refer to the first derivative as "velocity" and the second as "acceleration."

The maximum acceleration here is $2$; call that $A_0$. The distance along Y is also $2$; call that $D_0$.

We can scale this thing along X and Y to get a maximum acceleration of your choosing (call it $A_1$), over a desired distance (call it $D_1$). The X-scale and Y-scale can be calculated like this:

$$x_s = \sqrt{\frac{A_1D_0}{A_0D_1}}\\y_s = \frac{D_1}{D_0}$$

Now we'll modify that first equation to multiply all $x$ by $x_s$, and multiply the entire thing by $y_s$.

$$y = (x_sx)\left(2-\left|x_sx\right|\right)y_s$$

So, let's say you want the maximum acceleration to be $\frac{1}{4}$, and need to travel a distance of $1$.

$$x_s = \sqrt{\frac{A_1D_0}{A_0D_1}} = \sqrt{\frac{\frac{1}{4}\cdot2}{2\cdot1}} = \sqrt{\frac{1}{4}} = 0.5$$

$$y_s = \frac{D_1}{D_0} = \frac{1}{2} = 0.5$$

At this rate of acceleration, it takes four units of time (X axis) to travel one unit of distance (Y axis). You can calculate the time it will take like this:

$$t = 2\sqrt{\frac{A_0D_1}{A_1D_0}}$$

## Other shapes

The nice thing about this is you can use pretty much any shape you want and those calculations will still work. Let's say you want to use a sinusoidal curve instead of parabolic (for aesthetic reasons, maybe).

$$y = \sin\left(\frac{\pi}{2}x\right)$$

Again, we have it begin at $(-1,-1)$ and end at $(1,1)$.

The maximum acceleration here (our $A_0$) is $\frac{\pi^2}{4}$, or $\approx 2.4674$ as we can see from the plot. Our $D_0$ is always $2$. You can calculate the X-scale, Y-scale and time just like we did for the parabola; it'll still work fine.

Or instead of the sine shape, try something like this:

$$y = \frac{15x-10x^3+3x^5}{8}$$

Acceleration is zero at each end; no jerky "takeoffs" or "landings."

(work in progress...)

• Actually, the idea is to maximize the second derivative within the limitation constraint in order to minimize the duration of being off target. This leads to a constant maximum value of the 2nd which implies a parabola. The other big consideration is discrete vs continuous. The derivatives belong in the latter category and the problem is stated in the former. My accepted answer still doesn't have it quite right because there is a slight overshoot for non-commensurate values. Also, my solution doesn't require your assumptions and will still have no jerky "takeoffs" or "landings" – Cedron Dawg Apr 17 '18 at 3:12
• I am surprised that there isn't a stock answer to this. Since force is mass times acceleration, it would seem logical that many real world systems would have a force limitation on a mechanism and a requirement for no overshoots. – Cedron Dawg Apr 17 '18 at 3:14
• @CedronDawg, agreed on that first point, just thought it was interesting to look for a more generalized solution. By "jerky" I meant that jerk (third derivative) is infinite during takeoff and landing, because acceleration is undefined at those exact moments. That's also true of the halfway point... so if there is a "maximum second derivative" per OP, not sure whether sometimes-undefined second derivative should be allowed or not. If there's a stock "physical" answer to this, that discontinuity wouldn't be there... if you plot 2nd deriv of last eq. above, I think it might look similar to that. – Guest Apr 17 '18 at 6:04
• Hope that made sense, it's late here. Gonna update it tomorrow with the halfway point correction thing, haven't figured out if it's possible to do at other points with this method, but at halfway point we can take advantage of all even derivatives being zero – Guest Apr 17 '18 at 6:08
• I understand what you are saying. I made the same observation in my Apr 1st comment with "step function in the second derivative at the joints". The discrete vs continuous aspects come into play though. In the discrete case, the cap on the "jerk" is twice the maximum acceleration. If there needs to be a lower limit on that, then it only will have impact at the junctions. The rest of the time, for optimal results, the 2nd derivative should still be maximized, leaving the third at zero. – Cedron Dawg Apr 18 '18 at 1:13