smooth step function with rise time

Hi everyone,

I know this question is more related to math than Matlab, but I know there are many gurus out there who can easily solve this problem.

I would like to have a smooth step function with specified rise time, with some control of the output. I try to explain better. It is pretty simple to create a piece-wise linear step function and the correspondent smooth function with controlled slope using a parametrized logistic function:

% time axis
t = linspace(0,5,10000);

% step function params
Tr = 1;
a = 2;
t50 = 2;

% piece-wise linear step function
x = @(t)(t > t50-Tr/2 & t < t50+Tr/2).*(a/Tr.*(t-t50)+a/2)+(t >= t50+Tr/2).*a;

% logistic function
y = @(t)a./(1+exp(-4*(t-t50)/Tr));

% plot
figure, hold on
plot(t,x(t))
plot(t,y(t))
grid on

However I would like to control the shape of the rounded corners, making them more or less close to the original function according to an additional parameter (but keeping the slope of the original function) as suggested by the arrows in the following picture

I'm not in love with the logistic function, so any other function is welcome.

Thank you in advance for your help

Fabio

Answer 1

[Fabio Freschi]

My colleague Sofia have found a paper [1] where authors report the analytic approximation of the piece-wise linear step function with the so-called squashing function firstly introduced in [2].

This function has two parameters related to the position of the center γand the amplitude and one parameter that controls the curvature β:

The code below show this function at work

clear variables, close all

% time axis
t = linspace(-2,4,1000).';

% step function params
Tr = 2;    % rise time
A = 2;     % max value
T50 = 2;   % half rise point

% piece-wise linear step function
x = @(t)(t > T50-Tr/2 & t < T50+Tr/2).*(A/Tr.*(t-T50)+A/2)+(t >= T50+Tr/2).*A;

% squashing function
y = @(t,beta)A/Tr*log(((1+exp(beta*(t-T50+Tr/2)))./(1+exp(beta*(t-T50-Tr/2)))).^(1/beta));

% plot
figure, hold on, grid on
plot(t,x(t))
plot(t,y(t,5))
plot(t,y(t,20))
legend('linear','squashing \beta = 5','squashing \beta = 20','Location','best')

[1] A.I. Iliev, N. Kyurkchiev, S. Markov, On the Approximation of the Cut and Step Functions by Logistic and Gompertz Functions, Biomath Vol 4, No 2 (2015) 1-12, dx.doi.org/10.11145/j.biomath.2015.10.101

[2] J. Dombi and Z. Gera, The Approximation of Piecewise Linear Membership Functions and Lukasiewicz Operators, Fuzzy Sets and Systems 154(2) (2005) 275–286, http://dx.doi.org/10.1016/j.fss.2005.02.016

EDIT

The exponential term of the squashing function can become very large for large (positive) times (i.e. if the argument of the exp function is larger than about 700). I suggest to edit the function using simple algebra

Answer 2

[John D'Errico]

You cannot solve the problem you claim to want to solve using a simple logistic. It does not have the flexiibility to control both the amplitude of the rise, as well as the location and the slope at the midpoint. There are not sufficient free parameters in that functional form for the task at hand.

A simple variation on the question would allow you to write the problem in the form of a piecewise polynomial though.

Break it down into the subproblem: Can you find a low order polynomial that passes through the point with (x,y) coordinates (T50,1), with a given slope? As well, require that the curve must pass through a point (T100,2), with a slope of 0 at that point. See that since you could control how close T100 is to T50, then it would allow you to control the shape of the curve over that rise.

Once you have done the above task, then you would do exactly the same thing, effectively reflecting the form you just found to complete the bottom half of this curve.

A straight line is not sufficiently complex. But perhaps a 2nd or 3rd order polynomial would suffice. In fact, you can show that the polynomial must be at least a cubic to satisfy all of those conditions.

So, write down the requirements. A cubic polynomial has 4 coefficients. Call them a,b,c,d. I'll use syms here to do the work, since I'm too lazy myself. We would have:

syms a b c d T50 T100 t dT
P(t) = a*t^3 + b*t^2 + c*t + d;
DPdt = diff(P,t);
eq(1) = P(T50) == 1; % force the curve to pass through (T50,1)
eq(2) = P(T100) == 2; % force it through (T100,2)
eq(3) = DPdt(T100) == 0; % force the slope at the upper endpoint to be 0

Finally, what is the slope at the midpoint? It must match the slope of a straight line segment. So pick some dT, such that dT is greater than 0, but less than T100-T50. The line segment will be chosen so that it passes through the points (T50 - dT,0) and (T50 + dT,2). (This is the linear segment you have drawn. The slope of that line segment is 2/(2*dT) = 1/dT.

eq(4) = DPdt(T50) == 1/dT;

Now we can solve for a,b,c,d.

abcdsol = solve(eq,[a,b,c,d])

abcdsol = struct with fields:
    a: (T50 - T100 + 2*dT)/(dT*(T50 - T100)^3)
    b: -(T50^2 + T50*T100 + 3*dT*T50 - 2*T100^2 + 3*dT*T100)/(dT*(T50 - T100)^3)
    c: -(- 2*T50^2*T100 + T50*T100^2 - 6*dT*T50*T100 + T100^3)/(dT*(T50 - T100)^3)
    d: (2*dT*T50^3 - T50^2*T100^2 - 6*dT*T50^2*T100 + T50*T100^3 + 3*dT*T50*T100^2 - dT*T100^3)/(dT*(T50 - T100)*(T50^2 - 2*T50*T100 + T100^2))

And yes, the solution looks a bit messy, but that is why we have computers. It is also why I did not do this with pencil and paper.

Pupper = subs(P,abcdsol)

Pupper(t) = 

Again, a mess, but still easy enough to use.

First, I'll pick some values for T50, T100, and dT. In this example, I'll pick T50=2, and then we will leave T100 and dT to be controlled at a later time.

Pupperfun = matlabFunction(subs(Pupper,T50,2))

Pupperfun = function_handle with value:
    @(t,T100,dT)(t.^2.*1.0./(T100-2.0).^3.*(T100.*2.0+dT.*6.0+T100.*dT.*3.0-T100.^2.*2.0+4.0))./dT-(dT.*1.6e+1-T100.*dT.*2.4e+1+T100.^2.*dT.*6.0-T100.^3.*dT-T100.^2.*4.0+T100.^3.*2.0)./(dT.*(T100-2.0).*(T100.*-4.0+T100.^2+4.0))-(t.^3.*1.0./(T100-2.0).^3.*(-T100+dT.*2.0+2.0))./dT-(t.*1.0./(T100-2.0).^3.*(T100.*8.0+T100.*dT.*1.2e+1-T100.^2.*2.0-T100.^3))./dT

Now, plot the curve.

% line segment for comparison. Here, I'll choose 
T50 = 2;
dT = 1;
T100 = 4;
plot([T50-dT, T50, T50+dT],[0 1 2],'-r')
hold on

fplot(@(T) Pupperfun(T,3.5,dT),[T50,3.5],'g')
fplot(@(T) Pupperfun(T,4,dT),[T50,4],'b')
xlabel T
grid on

legend('Line segment','T100 = 3.5','T100 = 4','Location','SouthEast')

See how you can control the shape of the curve by adjusting T100 and dT. You can now reflect that upper functional form to the bottom half of the curve. I won't do that part as this is surely a student project of some form.

[Fabio Freschi]

@John D'Errico thank you for your answer. It really has enlightened me, so I spent some time in working with this idea.

I re-implemented your method, using a slightly different notation and observing that the coefficients of the polynomial function can be calculated solving a linear system

clear variables, close all

% time axis
t = linspace(1,4,1000).';

% step function params
Tr = 2;  % rise time
A = 2;     % max value
T50 = 2; % half rise point

% piece-wise linear step function
x = @(t)(t > T50-Tr/2 & t < T50+Tr/2).*(A/Tr.*(t-T50)+A/2)+(t >= T50+Tr/2).*A;

% control parameter
dT = 0.5; 

% polynomial of order 3: p0+p1*t+p2*t^2+p3*t^3
P = @(t)[ones(size(t)) t t.^2 t.^3];
% derivative of polynomial
dP = @(t)[zeros(size(t)) ones(size(t)) 2*t 3*t.^2];

% build linear system to impose
% 1. passage through (T50+Tr/2+dT,A). Note: T50+Tr/2 = T100
% 2. passage through (T50,A/2)
% 3. null derivative in T50+Tr/2+dT
% 4. prescribes slope in T50

% coefficient matrix
K = [P(T50+Tr/2+dT)  % 1
    P(T50)           % 2
    dP(T50+Tr/2+dT)  % 3
    dP(T50)];        % 4

b = [A               % 1
    A/2;             % 2
    0;               % 3
    A/Tr];           % 4

% solve to find the p_i coefficients of the polynomial
p = K\b;

% polynomial function
y = @(t)(t >= T50 & t < T50+Tr/2+dT).*P(t)*p;

figure, hold on, grid on
plot(t,x(t))
plot(t,y(t))
legend('linear','cubic')

And the results is identical to your green curve. However there is a problem. For values of the parameter dT, the maximum derivative is not located at T50. For example, by using dT = 0.3, che cubic polynomial intersects the original line segment.

I also realized that a similar methodology can be applied to write a higher order polynomial that covers the entire line segment, imposing a suitable number of constraints. In this case a quintic polynomial seems to be appropriate. The problem, again, is that for some values for the parameter dT the maximum slope can be not located at T50. Below I show a case where the derivative is not maximum at T50, generating multiple intersections with the straight line.

clear variables, close all

% time axis
t = linspace(1,4,1000).';

% step function params
Tr = 2;    % rise time
A = 2;     % max value
T50 = 2;   % half rise point

% piece-wise linear step function
x = @(t)(t > T50-Tr/2 & t < T50+Tr/2).*(A/Tr.*(t-T50)+A/2)+(t >= T50+Tr/2).*A;

% control parameter
dT = 0.15;

% polynomial of order 5: p0+p1*t+p2*t^2+p3*t^3+p4*t^4+p5*t^5
P = @(t)[ones(size(t)) t t.^2 t.^3 t.^4 t.^5];
% derivative of polynomial
dP = @(t)[zeros(size(t)) ones(size(t)) 2*t 3*t.^2 4*t.^3 5*t.^4];

% build linear system to impose
% 1. passage through (T50-Tr/2-dT,0)
% 2. passage through (T50+Tr/2+dT,A)
% 3. null derivative in T50-Tr/2-dT
% 4. null derivative in T50+Tr/2+dT
% 5. passage through (T50,A/2)
% 6. prescribed slope in T50: A/Tr

% coefficient matrix
K = [P(T50-Tr/2-dT)  % 1
    P(T50+Tr/2+dT)   % 2
    dP(T50-Tr/2-dT)  % 3
    dP(T50+Tr/2+dT)  % 4
    P(T50)           % 5
    dP(T50)];        % 6

% rhs
b = [0; A; 0; 0; A/2; A/Tr];

% solve to find the p_i coefficients of the polynomial
p = K\b;

% smooth-step function
y = @(t)(t > T50-Tr/2-dT & t < T50+Tr/2+dT).*P(t)*p+(t >= T50+Tr/2+dT).*A;
% derivative of smooth-step function
dy = @(t)(t > T50-Tr/2-dT & t < T50+Tr/2+dT).*dP(t)*p;

figure, hold on, grid on
plot(t,x(t))
plot(t,y(t))
plot(t,dy(t))
legend('linear','quintic','derivative')

Any suggestions to overcome this inconvenient?

Thank you once again for your help

Fabio