# Adding noise to nonlinear control system using NDSolve

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.

-
do I understand correctly, that you would like x'[t] == u[t] m[t] + rn[t]? –  user21 Jun 27 '13 at 10:00

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

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}];


-
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

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.

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].

-
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