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

Question

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)?

Answer 1

[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()

Answer 2

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]

Answer 3

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$$

Answer 4

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...)