Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Stochastic differential equations with Simulink?

Subject: Stochastic differential equations with Simulink?

From: Aino

Date: 6 Aug, 2013 14:03:09

Message: 1 of 5

Hi,

is it possible to simulate stochastic differential equations with Simulink? The one I'm working with is an inverted pendulum with two PD-controllers. The other controller has time delay in input and works only when certain criteria is met. This is the equation:

I*A''(t)=m*g*h*A(t)-[P1*A(t)+D1*A'(t)+f1(A(t-d))+f2(A'(t-d))]+s*n(t),

f1(A(t-d))=P2*A(t-d) & f2(A'(t-d))=D2*A'(t-d), when A(t-d)*[A'(t-d)-a*A(t-d)]>0 and
f1(A(t-d))=0 & f2(A'(t-d))=0 otherwise.

A(t) is the pendulum angle, I,m,g,h,s,P1,D1,P2,D2,a are constants, n(t) is Gaussian noise (mean=0, variance=1), t is time and d is time delay.

After some reformulation I believe this is what should be simulated:

x(n+1)=x(n)+f(x(n),x(n-k))dt+s*n(t)*sqrt(dt)

k=d/dt. The initial state is x0=[A(0),A'(0)]=[0.01,0].

Apparently the sqrt(dt) instead of just dt is the key problem here. If this is not doable with Simulink, any hint about how to make this happen is appreciated.

Thank you,
Aino

Subject: Stochastic differential equations with Simulink?

From: Phil Goddard

Date: 6 Aug, 2013 14:36:11

Message: 2 of 5


> I*A''(t)=m*g*h*A(t)-[P1*A(t)+D1*A'(t)+f1(A(t-d))+f2(A'(t-d))]+s*n(t),
>
> f1(A(t-d))=P2*A(t-d) & f2(A'(t-d))=D2*A'(t-d), when A(t-d)*[A'(t-d)-a*A(t-d)]>0 and
> f1(A(t-d))=0 & f2(A'(t-d))=0 otherwise.

Isn't the above just a non-linear ODE with Gaussian noise?
In which case it's relatively simple to implement.

>
> x(n+1)=x(n)+f(x(n),x(n-k))dt+s*n(t)*sqrt(dt)
> k=d/dt. The initial state is x0=[A(0),A'(0)]=[0.01,0].

I don't understand the term k=d/dt (a time delay specified as a partial derivative?), but in general if you use a fixed step solver then you know dt in advance and the sqrt(dt) term is just a (constant) modifier to the variance of the noise.

Phil.

Subject: Stochastic differential equations with Simulink?

From: Phil Goddard

Date: 6 Aug, 2013 16:09:19

Message: 3 of 5

Correction: I now see that the d in k = d/dt is the original time delay.

Phil.

Subject: Stochastic differential equations with Simulink?

From: Aino

Date: 7 Aug, 2013 10:16:07

Message: 4 of 5

The equation comes from "A model of Postural Control in Quiet Standing: Robust Compensation of Delay-Induced Instability Using Intermittent Activation of Feedback Control" (Asai et al. 2009). The authors call this SDE, with additive (not multiplicative) noise and say that the integration should be performed in terms of Ito integral.

I hope my equation is not too confusing, for example in the article the time delay d is actually '\Delta' and the dt is '\Deltat'.

I'm still studying these creatures, but does anyone have any idea how to implement stochastic (delay) differential equations with Matlab? I can do this in Simulink without the noise part, but not with it. If(!) I have understood correctly I cannot just include noise to my simulation and be happy.

-Aino

Subject: Stochastic differential equations with Simulink?

From: Aino

Date: 29 Aug, 2013 12:15:08

Message: 5 of 5

I'll try to answer myself. :)

I don't think Simulink can be used for SDDEs, and if(!) I got this right, the answer is so simple that Simulink wouldn't be the best choice anyway. In fact, I think this question is now more "how to solve SDDEs with Euler's method" than how to use Matlab to do this, but I'll post my code anyway to get some closure.. Also, I haven't verified the code below (I created it by combining bits and pieces of information on how to solve ODEs, SDEs, DDEs.. etc.), so if some SDDE wiz would take a look at it, it would be greatly appreciated.

So here is the code:

%-------------------------------------------------
clear all;close all;clc;

%%%%%%%%%%%%%%%%%
%Parameters
%%%%%%%%%%%%%%%%%
TrialLength=100;%sec
SampleFrequency=1000;%Hz
dt=1/SampleFrequency;%sec
TimeVec=(dt:dt:TrialLength);
m=60;%kg
g=9.81;%m/s2
h=1;%m
J=60;%kgm2
Slope_a=-0.4;
P=0.25*m*g*h;%Nm/rad
D=10;%Nms/rad
K=0.8*m*g*h;%Nm/rad
B=4;%Nms/rad
TimeDelay=0.2;%s
NoiseAmp=0.2;%Nm
r=0.004;%rad-rad/s
k=TimeDelay/dt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Initial conditions and allocation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X(1:TrialLength*SampleFrequency)=0;%Angle
X(k+1)=0.01;
V(1:TrialLength*SampleFrequency)=0;%Angle velocity

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Sway Angle with Euler's method:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n=k+1:TrialLength*SampleFrequency+k
    
    a=((m*g*h)-K)/J;
    b=-B/J;
    c=-P/J;
    d=-D/J;
    e=NoiseAmp/J;
    
    N=randn(1,1);%Gaussian random number, mean=0, var=1
    
    %Euler's method (hopefully correct..)
    %Angle
    X(n+1)=X(n)+V(n)*dt;
    %Angular velocity
    V(n+1)=V(n)+[a*X(n)+b*V(n)+c*X(n-k)+d*V(n-k)]*dt+e*N*sqrt(dt);

    %PD active control on/off
    if X(n-k)*(V(n-k)-Slope_a*X(n-k))>0 && X(n-k)^2+V(n-k)^2>r^2
        P=0.25*m*g*h;
        D=10;
    else
        P=0;
        D=0;
    end
    
end

%Remove initial transient:
X=X(k+2:end);
V=V(k+2:end);

%Plotting
figure;plot(X,V);ylabel('Angular velocity, (rad/sec)');xlabel('Angle (rad)');
figure;plot(TimeVec,X);ylabel('Angle (rad)');xlabel('Time (s)')

%-------------------------------------------------

-Aino

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us