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
采纳的回答
Walter Roberson
2021-1-16
0 个投票
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.
10 个评论
l l
2021-1-16
l l
2021-1-16
Walter Roberson
2021-1-16
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)
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]);
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
2021-1-16
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
2021-1-18
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')
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]);
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
2021-1-18
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
2021-1-20
Walter Roberson
2021-1-20
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.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Data Type Identification 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
