How can I code pseudo random binary sequence input current signal made by idinput?

I tried to define a prbs current signal input as function of time, something like:

u=idinput(100,'prbs',[0 1],[1 3]);
u=iddata([],u,1);
I(t)=u;

so when I give t get I(t) at that time to use in my ode equation:

function rate=rate_equ(t,y)

global   q Vg Nsd Ssd Beta Gamma A0 N0 Epsil Tp Tc Vact Alpha 

%---- Computing the Rate Equations ------------------------------
rate = zeros(2,1);
rate(1) = (Gamma*A0*Vg*((y(2)-N0)/(1+Epsil*y(1)))-1/Tp)*y(1)+Beta*Gamma*y(2)/Tc + ...
           randn*sqrt(2*Beta*Vact*Nsd*(Vact*Ssd+1)^3/Tc) ;
rate(2) = I(t)/(q*Vact)-y(2)/Tc-Vg*A0*y(1)*(y(2)-N0)/(1+Epsil*y(1)) + ... 
          randn*sqrt(2*Vact*Nsd/Tc*(Beta*Vact*Ssd+1)) ;
end

Answer 1

[Walter Roberson]

When your functions use random numbers, such as your call to randn(), then they are never ODE (Ordinary Differential Equations).

The tools for Stocastic Differential Equations are in the Finance toolbox; https://www.mathworks.com/help/finance/sde-objects.html

If you were to remove the randn() then you would still have problems. The ode*() routines are only valid when the functions you provide are twice differentiable, which is never the case when you move abruptly between two states as in your PRBS data.

The ode*() routines are not fixed step-size routines: they are variable step size. So the function will not just be executed with respect to time 0, dt, 2*dt, 3*dt, and so on. You therefore need some way of mapping from the input time (first value passed to the ode function) and the sample number, such as by using interp1(). But! Interp1() use 'linear' by default, and the values produced by 'linear' or 'previous' or 'next' interpolation methods are not twice differentiable. And if you use interp1() with spline, then that is twice differentiable, but it will also have intermediate values between 1 and 3 (your signal bounds), which is incompatible with the idea of using a pseudo random binary signal.

[l l]

Thank you for your reply,

But yet I don't uderstand how to solve this problem and give such a current to the equation!

[l l]

I want first make prbs current signal and then use that signal in ode.

[Walter Roberson]

In your mathematics, what is the solution to this equation:

randn*x == 5

? And what is

f = @(x) randn(size(x)).*x
integral(f, 0, 1)

[Walter Roberson]

rng('shuffle');

%poster did not supply values for the global variables, and
%gave no documentation about their range. Toss in random values,
%just to have SOMETHING to test with.
P.q = rand;
P.Vg = rand;
P.Nsd = rand;
P.Ssd = rand;
P.Beta = rand;
P.Gamma = rand;
P.A0 = rand;
P.N0 = rand;
P.Epsil = rand;
P.Tp = rand;
P.Tc = rand;
P.Vact = rand;
P.Alpha = rand;

%prbs bandwidth of [0 1] specifies that the signal is to change every
%clock period. How long is a clock period? The iddata() call
%tells us that the clock sample time is 1.

u = idinput(100,'prbs',[0 1],[1 3]);

Warning: The PRBS signal delivered is the 100 first values of a full sequence of length 127.

u = iddata([],u,1);
P.I = u.u;
P.T = 0:length(P.I)-1;  %sample time 1 second

stairs(P.T, P.I); ylim([0.9 3.1]); title('PRBS');

y0 = rand(1,2);

[t, y] = ode45(@(t,y) rate_equ(t, y, P), linspace(0,10,50), y0);
plot(t, y); legend('y1', 'y2'); title('try 1')

[t, y] = ode45(@(t,y) rate_equ(t, y, P), linspace(0,10,50), y0);
plot(t, y); legend('y1', 'y2'); title('try 2')

function rate=rate_equ(t, y, P)

%---- Computing the Rate Equations ------------------------------
    rate = zeros(2,1);
    It = interp1(P.T, P.I, t);
    rate(1) = (P.Gamma * P.A0 * P.Vg * ((y(2)- P.N0) / (1 + P.Epsil * y(1))) - 1/P.Tp) * y(1) + P.Beta * P.Gamma * y(2)/P.Tc + ...
               randn * sqrt(2 * P.Beta * P.Vact * P.Nsd * (P.Vact * P.Ssd +1)^3 / P.Tc) ;
    rate(2) = It / (P.q * P.Vact) - y(2) / P.Tc - P.Vg * P.A0 * y(1) * (y(2) - P.N0) / (1 + P.Epsil * y(1)) + ... 
              randn * sqrt(2 * P.Vact * P.Nsd / P.Tc * (P.Beta * P.Vact * P.Ssd + 1)) ;
end

[Walter Roberson]

I am a bit surprised that ode45() does not complain about singularities. As some of my runs timed out, I suspect it can happen.

I am also surprised that despite the randn() the two attempts come out looking fairly close to each other.

[l l]

I tried the code but it doesn't work. when I use interp1 on (P.T,P.I) it only gives me one for each value of t, but I want as stairs plot for example from 0 to 6.5 one and from 6.5 to 8 three and so on.

plot(P.T,P.I);hold on
stairs(P.T,P.I,'k');
t=linspace(0,100,500);
It = interp1(P.T,P.I,t);
plot(It,'bo');legend('P.T-P.I','stairs','interp1')

[Walter Roberson]

rng('shuffle');

%poster did not supply values for the global variables, and
%gave no documentation about their range. Toss in random values,
%just to have SOMETHING to test with.
P.q = rand;
P.Vg = rand;
P.Nsd = rand;
P.Ssd = rand;
P.Beta = rand;
P.Gamma = rand;
P.A0 = rand;
P.N0 = rand;
P.Epsil = rand;
P.Tp = rand;
P.Tc = rand;
P.Vact = rand;
P.Alpha = rand;

%prbs bandwidth of [0 1] specifies that the signal is to change every
%clock period. How long is a clock period? The iddata() call
%tells us that the clock sample time is 1.

u = idinput(100,'prbs',[0 1],[1 3]);

Warning: The PRBS signal delivered is the 100 first values of a full sequence of length 127.

u = iddata([],u,1);
P.I = u.u;
P.T = 0:length(P.I)-1;  %sample time 1 second

plot(P.T, P.I); hold on
stairs(P.T,P.I,'k');
t=linspace(0,100,500);
It = interp1(P.T,P.I,t,'prev');
plot(t, It,'bo');legend('P.T-P.I','stairs','interp1')
xlim([0 16]); ylim([0.9 3.1]); title('PRBS');
hold off

y0 = rand(1,2);

[t, y] = ode45(@(t,y) rate_equ(t, y, P), linspace(0,10,50), y0);
plot(t, y); legend('y1', 'y2'); title('try 1')

[t, y] = ode45(@(t,y) rate_equ(t, y, P), linspace(0,10,50), y0);
plot(t, y); legend('y1', 'y2'); title('try 2')

function rate=rate_equ(t, y, P)

%---- Computing the Rate Equations ------------------------------
    rate = zeros(2,1);
    It = interp1(P.T, P.I, t, 'previous');
    rate(1) = (P.Gamma * P.A0 * P.Vg * ((y(2)- P.N0) / (1 + P.Epsil * y(1))) - 1/P.Tp) * y(1) + P.Beta * P.Gamma * y(2)/P.Tc + ...
               randn * sqrt(2 * P.Beta * P.Vact * P.Nsd * (P.Vact * P.Ssd +1)^3 / P.Tc) ;
    rate(2) = It / (P.q * P.Vact) - y(2) / P.Tc - P.Vg * P.A0 * y(1) * (y(2) - P.N0) / (1 + P.Epsil * y(1)) + ... 
              randn * sqrt(2 * P.Vact * P.Nsd / P.Tc * (P.Beta * P.Vact * P.Ssd + 1)) ;
end

[Walter Roberson]

Reminder: ode45() used with a function that has discontinuous jumps will always produce the wrong answer. Sometimes you get unlucky and the ode does not produce an error message, and you are left with the mistaken impression that you got a valid result.

[l l]

My thanks and appreciation for your guidance.It works well for me.

Another question that I have is that I actually want to use this procedure in dde23 function

but I don't know how to set 'Jumps' to avoid wrong answer of discontinuous.

[Walter Roberson]

You are using random numbers. The system is discontinuous. If you are fortunate, you would get explicit error messages telling you that your functions are not suitable for the ode functions that you are using.