Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a nonlinear control system

system = {x'[t] == u[t] m[t], 
 WhenEvent[And[Mod[t, 9], x[t] < 100], m[t] -> x[t]], 
 WhenEvent[u[t] == 0, m[t] -> 10],  x[t /; t <= 0] == 0, 
 m[t /; t <= 0] == 1, u[t /; t <= 0] == 0};

control = {e[t] == xref - x[t], 
 WhenEvent[Mod[t, \[Tau]], 
 u[t] -> kp (k1 Sign[e[t]] - k2 Sign[e'[t]])], 
 e[t /; t <= 0] == xref};
 params = {kp -> 1, td -> 1, \[Tau] -> 0.5, xref -> 810, k1 -> 11/2, k2 -> 9/2};

which I simulate using NDSolve

 sol = NDSolve[{system, control} /. params, {x, m, u}, {t, 0, 60}, 
  DiscreteVariables -> {u, m}];

I would like to add random noise to the system, but so far nothing seems to work. Using RandomReal[NormalDistribution[]] is problematic, as it only samples once for NDSolve. Using WienerProcess with ItoProcess also does not seem to work, as the system contains discrete variables. I cannot convert it to a state space model, as StateSpaceModel only works for linear systems.

share|improve this question
1  
do I understand correctly, that you would like x'[t] == u[t] m[t] + rn[t]? –  user21 Jun 27 '13 at 10:00
add comment

2 Answers

up vote 3 down vote accepted

I can suggest you to add noise like the following!

noise[n_] := noise[n] = RandomReal[{-10, 10}];
MyNoise1[t_?NumericQ] :=1.5 Evaluate@With[{dt = .4},
    noise[Floor[t/dt]] Sin[2 Pi Mod[t, dt]/dt]^2];

or you can use something like this

sigma = 4.5;
{start, end} = {0, 60};
n = 240;
MyNoise2 = Interpolation@Join[{{start, 0.}}, 
Rest@Table[{t, RandomReal@NormalDistribution[0, sigma]}, {t, 
   start, end, (end- start)/(n - 1)}]];

The noise looks like the following

enter image description here

Now we add the noise like this

With[{v = MyNoise1},
     system = {x'[t] == u[t] m[t] + v[t], 
     WhenEvent[And[Mod[t, 9], x[t] < 100], m[t] -> x[t]], 
     WhenEvent[u[t] == 0, m[t] -> 10], x[t /; t <= 0] == 0, 
     m[t /; t <= 0] == 1, u[t /; t <= 0] == 0}];

enter image description here

share|improve this answer
    
nice answer +1! –  user21 Jun 27 '13 at 11:12
    
@ruebenko Thx! Finding the sensitivity of the response with respect to such random noise will be great. –  PlatoManiac Jun 27 '13 at 11:35
    
I have never tried that, you could give it a go with ParameticNDSolve I guess. I'd use a parameter for the noise and compute the derivative. –  user21 Jun 27 '13 at 12:08
    
Great, that seems to work well. However, I had to increase MaxSteps -> 30000 to be able to plot it? –  Gerrit Jun 27 '13 at 13:17
add comment

Another way to generate the noise would be to use RandomVariate.

noise = Interpolation[Thread[{Range[0, 60, 0.5], 
Join[{0}, RandomVariate[NormalDistribution[], 120]]}], t]

This is how the noise looks.

enter image description here

Then everything follows exactly as the answer by PlatoManiac. (assuming that the noise enters the system through the input)

With[{v = noise}, system = {x'[t] == u[t] m[t] + v, 
WhenEvent[And[Mod[t, 9], x[t] < 100], m[t] -> x[t]], 
WhenEvent[u[t] == 0, m[t] -> 10], x[t /; t <= 0] == 0, 
m[t /; t <= 0] == 1, u[t /; t <= 0] == 0}];

etc

And the plot of x[t].

enter image description here

share|improve this answer
    
It seems that the noise value is constant for different values of t, or at least when I try to plot it? –  Gerrit Jun 27 '13 at 13:16
    
Are you plotting noise =... and getting a constant value? –  Suba Thomas Jun 27 '13 at 13:30
    
No, when I plot x[t] and e[t] the segments between events/state changes are straight lines. This shows up clearer when you plot from 0-30 instead of 0-60. Is it possible that it should be v[t] instead of just v? –  Gerrit Jun 27 '13 at 13:38
    
v is correct not v[t]. My guess that it appears to be joined by straight lines is because Interpolation is used? Not sure. But the noise value does change from instant to instant. –  Suba Thomas Jun 27 '13 at 14:04
    
Sorry, you're right. I increased the $\sigma$ of the distribution and I can see the changes. –  Gerrit Jun 27 '13 at 14:31
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.