Random generator seed for parallel simulation using fast restart

Question

Niklas Larsson 2021-9-10

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1450814-random-generator-seed-for-parallel-simulation-using-fast-restart

评论： Niklas Larsson 2021-9-27

I'm going to perform a monte carlo simulation by simulate a simulink model thousands of times. The model are built by blocks from SimEvent-toolbox and some standard blocks.

I use one simulink function block that are calling a MATLAB-function block. In here I try to simulate the break down of sensors during one year of service in the machine (which is the rest of the model).

MATLAB-function:

function [outputPort,log] = failureRate(time,lambda,nSensors)
    persistent failed failedCount
    
    if isempty(failed)
        failed = false(1,nSensors);
        failedCount = 0;
        rng('shuffle')
    end
    % Failure rate function and random variable
    failRate = @(t,lam) 1 - exp(-t*lam);
    r = rand(1,nSensors);
    nBroken = 0;
    
    % Check if sensors break or not
    if any(r < failRate(time,lambda))
        tmp = failed;
        failed = max(failed, r < failRate(time,lambda));
        if failedCount < sum(failed)
            failedCount = sum(failed);
            nBroken = sum(failed - tmp);
        end
        
    end
    
    if all(failed)
        % All sensors broken
        outputPort = 2;
        log = [nBroken,time];
    else
        % Some sensors still working
        outputPort = 1;
        log = [0, time];
    end
end

This function is called each time an entity enters a "Entity Server"-block, which will decide which path it will go (failed sensor path or not).

I can run this model using a normal for loop using fast restart but as soon as I try to use Fast restart for any parallel computing the randomization becomes repeated.

I've tried following:

for i = 1:n
    simIn(i) = Simulink.SimulationInput('SensorFailure');
    simIn(i) = simIn(i).setPreSimFcn(@(x) rng('shuffle'));
end

but it does not work when Fast restart is actiaved using the parsim command. Also I've tried with and without the "SetupFnc" option in parsim (and all different combinations of these three.

out = parsim(simIn,...
    'ShowSimulationManager','off',...
    'ShowProgress','on',...
    'TransferBaseWorkspaceVariables','on',...
    'UseFastRestart','on',...
    'SetupFcn',@() rng(randi(9999999)),...
    'StopOnError','on');

I also tried using parfor loop but as long as fast restart is active the randomization will be repeated.

parfor i = 1:n
    rng('shuffle') %tried with shuffle and randi(10000)
    out = sim('SensorFailure','FastRestart','on');
    data{i} = out.sensorfailed(out.sensorfailed(:,1) > 0,2);
    if mod(i,10) == 0
        disp(i)
    end
end

Since I'm looking to run thousands of simulation and problably things will need to be redone once the results are analyzed any gain in simulation time is worth some time to investigate!

I've tried using rng('shuffle') and rng(randi(100000)) for all possible situations since I've understood that "shuffle" is not optimal for parallel computing. I've also tried these commands in the simulink models callback (initializing and start callbacks)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Niklas Larsson 2021-9-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1450814-random-generator-seed-for-parallel-simulation-using-fast-restart#answer_786279

在 MATLAB Online 中打开

After a while I managed to find a solution that solved my problem!

By adding one more input parameter to the matlab function and feeding that input with a constant value-block:

function [outputPort,log] = failureRate(time,lambda,nSensors,seed)
persistent failed failedCount
if isempty(failed)
    failed = false(1,nSensors);
    failedCount = 0;
    rng(seed) % Initiate with seed from constant block outside of function.
end
% Failure rate function and random variable
failRate = @(t,lam) 1 - exp(-t*lam);
r = rand(1,nSensors);
nBroken = 0;
if any(r < failRate(time,lambda))
    tmp = failed;
    failed = max(failed, r < failRate(time,lambda));
    
    log = [0,0];
    if failedCount < sum(failed)
        failedCount = sum(failed);
        nBroken = sum(failed - tmp);
    end
    
end
if all(failed)
    % All sensors broken
    outputPort = 2;
    log = [nBroken,time];
else
    % Still working
    outputPort = 1;
    log = [nBroken > 0, time];
end
end

Then for each simulation in the SimulationInput-object randomize the seed value as below:

    seedBlockPath = 'SensorFailure/Subsystem/Simulink Function/Seed';
    for i = 1:n
        % Set random seed for each simulation
        simIn(i) = Simulink.SimulationInput('SensorFailure');
        r = randi(99999); 
        simIn(i) = simIn(i).setBlockParameter(seedBlockPath,'Value',num2str(r));
    end
    out = parsim(simIn,...
        'ShowSimulationManager','off',...
        'ShowProgress','on',...
        'TransferBaseWorkspaceVariables','on',...
        'UseFastRestart','on',...
        'StopOnError','on');

13 个评论
显示 11更早的评论隐藏 11更早的评论

Paul 2021-9-21

编辑：Paul 2021-9-21

在 MATLAB Online 中打开

I understand that you are modeling the time-of-failure (T) as a random variable with an exponential distribution. Its CDF is

P(T < t) = F(t) = 1 - exp(-lambda*t), for t>=0, and 0 otherwise. The parameter lamda is constant and is the failure rate. The expected value of T is E(T) = 1/lambda.

Its PDF, the derivative of the CDF, is

f(t) = lambda*exp(-lambda*t), for t>=0, and 0 otherwise.

So the question becomes, how does one generate the failure time, T, in a simulation such that across many simulations T has an exponential distribution. Assuming (and this is a very important assumption) that the function failureRate() is called with a constant sample time, the code below emulates the process in the simulation and workflow as I understand it.

lambda = 0.5; % example failure rate

nruns = 1000; % number of Monte Carlo runs

dt = 0.01; % fixed sample time of failurRate()

failuretime = nan(nruns,1); % pre-allocate results

rng(100); % for repeatability

for ii = 1:nruns % Monte Carlo loop

% initialize simulation

time = 0;

failed = false;

while ~failed % run simulation

failed = rand < 1 - exp(-lambda*time); % the test

if failed

failuretime(ii) = time;

else

time = time + dt;

end

% check the mean time of failure. Should be 1/0.5 = 2

mean(failuretime)

ans = 0.1753

% plot the histogram and compare to theoretical pdf

figure

histogram(failuretime,'Normalization','pdf')

hold on

tvals = 0:.01:8;

plot(tvals,exppdf(tvals,1/lambda))

xlim([0 2])

We see that the results of the simulation do not match expected value nor the desired pdf. The reason is that the test for failure is using time since t = 0. However, at each time step where the item is determined to have not failed, the probability "resets," which is a property of the exponential random variable. So the test needs to change as in the following code:

failuretime = nan(nruns,1); % pre-allocate results

rng(100); % for repeatability

for ii = 1:nruns % Monte Carlo loop

% initialize simulation

time = 0;

failed = false;

while ~failed % run simulation

failed = rand < 1 - exp(-lambda*dt); % the modified test, use dt

if failed

failuretime(ii) = time;

else

time = time + dt;

end

% check the mean time of failure. Should be 1/0.5 = 2

mean(failuretime)

ans = 2.0654

% plot the histogram and compare to theoretical pdf

figure

histogram(failuretime,'Normalization','pdf')

hold on

tvals = 0:.01:8;

plot(tvals,exppdf(tvals,1/lambda))

Now we have a good match between the simulation and the theoretical results (increase nruns to make it match even better). Note that failuretime(ii) is not really the time of failure, it's really the trailing edge of the window, defined by failuretime(ii) - dt to failuretime(ii), in which the failure occurred. So this approach is really an approximation because it's not telling you whent the failure occurred. But it might be a reasonable approximation for your application depending on the value of dt relative to failure rate.

Another option that may be applicable (depends on what you're actually simulating) is to avoid the dynamic computation of failuretime in the simuation altogether. Instead, compute the failure times a priori with exprnd() and use them as input parameters to the simulation.

Niklas Larsson 2021-9-24

在 MATLAB Online 中打开

The main problem for getting the mean value correct was just as you said, the assumption that the function is called at each time step and changing the time to dt instead of the current time.

Before was it only called once an entity arrived at a "Entity Server"-block, which was far from each time step by natural cause. I just moved out the MATLAB-function and let it run "freely" in the simulation as well as change the solver to a fixed step solver.

Se result below. I would say it is reasonable that using 10 sensors instead of 1 would increase the risk of failure since more sensors are at risk of breaking at each time step while the machine is running. The reason I do this is due to my application. The Fail Rate units is [1/hours] and MTBF in [hours]

FYI: When I changed the simulation to run the function at each timestep the simulations of course run slower but also running multiplie simulation runs (like 1000 simulation for 1 sensor and a certan fail rate, then another 1000 simultion with 10 sensors and a certain fail rate) were not possible for my computer. To fix it I had to turn of the command

'TransferBaseWorkspaceVariables'

and just input all variables using SimulationInput commands. Since the saved data increased drastically due to all the simulation steps I guess I ran out of memory.

Paul 2021-9-24

Upon further thought, I'd like to change my very important assumption to be that the function failureRate() is called once per time step. So the variable step discrete solver should work fine as long as

the function failureRate keeps track of (or gets an input that is) the dt on each call, and

the dt on every time step is small enough so that the model still achieves a good approximation to the exponential distribution.

Niklas Larsson 2021-9-27

I guess that would work too if, as you say, the dt inbetween the calls are small enought. For my case I dont think it is.

Either way, thank you very much for your time and expertise!

请先登录，再进行评论。

Answer 2

Abdolkarim Mohammadi 2021-9-10

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1450814-random-generator-seed-for-parallel-simulation-using-fast-restart#answer_784804

编辑：Abdolkarim Mohammadi 2021-9-10

I do not know the solution to your problem, but I just want to mention that I faced a problem with SimEvents and parsim() with fast restart a while ago in R2020a, and I wrote about it here.

https://www.mathworks.com/matlabcentral/answers/593833-parsim-with-fast-restart-returns-different-results

They reported it as a bug report here.

Bug Reports (mathworks.com)

But the problem persists in R2021a. Maybe your problem relates to this bug.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Edric Ellis 2021-9-13

That bug reports indicates the problem was fixed in R2020b. If you can reproduce your problem in R2020b (or later), then something else must be going wrong - please contact MathWorks support with your reproduction steps.

请先登录，再进行评论。

Random generator seed for parallel simulation using fast restart

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

13 个评论
显示 11更早的评论隐藏 11更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Random generator seed for parallel simulation using fast restart

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

13 个评论 显示 11更早的评论隐藏 11更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

13 个评论
显示 11更早的评论隐藏 11更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论