Reliability of Grid Networks and Recursive Decomposition Algorithms

Rice University, Department of Civil and Environmental Engineering.

Structural & Infrasturctural System Reliability & Risk Assessment (SISRRA) Research Group, 2016.

In this fiile we show how to use a Matlab implementation of the recursive decomposition algorithm for computing the s-t reliability in squared network grids. For running this example you'll need Matlab implementations files from SISRRA. Contact infromation can be found at the website of the group: https://duenas-osorio.rice.edu/.

(For best results use Internet Explorer)

NxN grids and reliability estimates between farthest diagonal pairs (Benchmark - Within 1,000 seconds - Analytical - Edge Failure Probability of 10%)

The table below shows upper and lower bound estimates of the two-terminal reliability for square grids with N-by-N nodes. The two pairs are selected at the ends of the square grids diagonals (see example of a 5x5 grid). The table also shows gap or distance between upper and lower bounds (e.g. gap = 0 means exact estimate), and running times (in seconds) are shown for reference (note truncation at 1,000 seconds).

format long

readtable('grid_table_1000.txt')

ans =

     N       LowerBound        UpperBound           Gap          CPU_Time_secs
    ___    ______________    ______________    ______________    ______________

      3    0.972502171407    0.972502171407                 0       0.008303088
      4    0.975046349577    0.975046349577                 0       1.203886717
      5    0.975556585069    0.975556601318        1.6249e-08    1000.000030948
     10    0.966242133758       0.976891412    0.010649278241    1000.000066173
     20    0.404469395567    0.987050144558    0.582580748991    1000.000027598
     30    0.040080180363    0.988199973472    0.948119793109    1000.000277193
     40    0.003468759545    0.988199996775    0.984731237229    1000.000117628
     50    0.000320686568    0.988199996079    0.987879309511     1000.00149665
     75       9.12253e-07    0.989999999794     0.98999908754    1000.008133252
    100         4.338e-09              0.99    0.989999995661    1000.000047891

NxN grids and reliability estimates between farthest diagonal pairs (Benchmark - Within 7.8 Hours - Analytical - Edge Failure Probability of 10%)

Table below is identical to the previous except for timeout threshold that has been set to 7.8 Hours.

readtable('grid_table_28080.txt')

ans =

     N       LowerBound        UpperBound           Gap           CPU_Time_secs
    ___    ______________    ______________    ______________    _______________

      3    0.972502171407    0.972502171407                 0        0.008301122
      4    0.975046349577    0.975046349577                 0        1.254786477
      5    0.975556589561    0.975556589475                 0     1405.522161999
     10    0.968901066035    0.976891411735      0.0079903457    28080.000063171
     20    0.433024435707    0.986947343626    0.553922907919    28080.000050634
     30    0.054874120262    0.987884076283    0.933009956021    28080.000052926
     40    0.004457073472    0.988199996796    0.983742923324    28080.000916846
     50    0.000408723505    0.988199996079    0.987791272574    28080.001586879
     75      1.240808e-06     0.98999999998    0.989998759172    28080.004209904
    100         4.338e-09              0.99    0.989999995661     28080.00233887

NxN grids and reliability estimates between farthest diagonal pairs (Benchmark - Within 8 Hours - Monte Carlo Simulation - Edge Failure Probability of 10%)

Using Monte Carlo Simulation we have simulated enough samples to guarantee that the reliability estimate has relative error below 1% with 99% confidence.

readtable('grid_table_mcs.txt')

ans =

     N     Confidence    Deviation       Estimate       CPU_Time_secs
    ___    __________    _________    ______________    ______________

      3    0.99          0.01         0.972874806801       0.444983707
      4    0.99          0.01          0.97117724125       1.568695888
      5    0.99          0.01         0.978684763907       2.475756175
     10    0.99          0.01         0.977558262204      25.572117525
     20    0.99          0.01         0.976907811654     143.077300822
     30    0.99          0.01         0.978163668504     265.384622704
     40    0.99          0.01         0.974636363636     398.185984023
     50    0.99          0.01         0.979324603722     509.095017926
     75    0.99          0.01         0.976257362936     729.075298347
    100    0.99          0.01         0.975250385208    1060.408605218

NxN grids and Failure Probability estimates between farthest diagonal pairs (Benchmark - Within 1,000 seconds - Analytical - Edge Failure Probability of 1%)

Table below shows s-t network failure probability. Ttimeout threshold that has been set to 1,000 seconds.

readtable('grid_table_fp_1000.txt')

ans =

     N       UpperBound        LowerBound           Gap          CPU_Time_secs
    ___    ______________    ______________    ______________    ______________

      3    0.000207985977    0.000207985977                 0       0.008193354
      4    0.000204091162    0.000204091163                 0       1.198138911
      5     0.00020403122     0.00020403122             1e-12    1000.000005063
     10    0.000204033254     0.00020402989         3.364e-09    1000.000241555
     20    0.000820989809     0.00019744058    0.000623549229     1000.00017032
     30    0.019713732197     0.00015780667    0.019555925527    1000.000453033
     40    0.046537641873    0.000147641021    0.046390000852    1000.000287839
     50    0.219849550634    0.000100003809    0.219749546825    1000.000703444
     75    0.435765881216    0.000100002305    0.435665878911    1000.000160812
    100    0.633630456136       1.38039e-07    0.633630318097    1000.005915499

NxN grids and Failure Probability estimates between farthest diagonal pairs (Benchmark - Within 7.8 Hours - Analytical - Edge Failure Probability of 1%)

Table below is identical to the previous except for timeout threshold that has been set to 7.8 Hours.

readtable('grid_table_fp_28080.txt')

ans =

     N       UpperBound        LowerBound           Gap           CPU_Time_secs
    ___    ______________    ______________    ______________    ________________

      3    0.000207985977    0.000207985977                 0         0.008182296
      4    0.000204091162    0.000204091163                 0          1.21664125
      5    0.000204031226    0.000204031227                 0      1392.676070669
     10    0.000204033254     0.00020402989         3.364e-09    28080.0001378279
     20    0.000820989809     0.00019744058    0.000623549229     28080.000221736
     30    0.019713173875    0.000157806671    0.019555367204    28080.0002575959
     40    0.046537641873    0.000147641021    0.046390000852    28080.0000827279
     50    0.219849550634    0.000100003809    0.219749546825    28080.0006459049
     75    0.435765881216    0.000100002305    0.435665878911     28080.000368343
    100    0.633630456136       1.38039e-07    0.633630318097    28080.0092558949

NxN grids and Failure Probability estimates between farthest diagonal pairs (Benchmark - Within 8 Hours - Monte Carlo Simulation - Edge Failure Probability of 1%)

Using Monte Carlo Simulation we have attempted to simulate enough samples to guarantee that the failure probability estimate has relative error below 1% with 99% confidence. As shown in the table below, the rare event condition has prevented Monte Carlo Simulation from giving an answer with approximation guarantees within the time threshold.

readtable('grid_table_fp_mcs.txt')

ans =

     N     Confidence    Deviation    Estimate    CPU_Time_secs
    ___    __________    _________    ________    _____________

      3    0.99          0.01         NaN         28800
      4    0.99          0.01         NaN         28800
      5    0.99          0.01         NaN         28800
     10    0.99          0.01         NaN         28800
     20    0.99          0.01         NaN         28800
     30    0.99          0.01         NaN         28800
     40    0.99          0.01         NaN         28800
     50    0.99          0.01         NaN         28800
     75    0.99          0.01         NaN         28800
    100    0.99          0.01         NaN         28800

NxN grids and Failure Probability estimates between farthest diagonal pairs (Benchmark - Within 1,000 seconds - Monte Carlo Simulation - Edge Failure Probability of 1%)

Below we show mean and variance for s-t network failure probability estimates using Monte Carlo Simulation. Note that this estimates have no guarantees of quality of approximation.

readtable('grid_table_fp_mcs_1000.txt')

ans =

     N     Variance          Mean         CPU_Time_secs    n_iterations
    ___    _________    ______________    _____________    ____________

      3        5e-12    0.000205192972    1000             41775310
      4        9e-12    0.000203899659    1000             21834269
      5      1.6e-11     0.00020488716    1000             13148701
     10      6.6e-11    0.000192990688    1000              2906876
     20     3.07e-10    0.000206661949    1000               672596
     30     7.84e-10    0.000215154145    1000               274222
     40    1.226e-09    0.000194931774    1000               159030
     50    1.504e-09    0.000150229852    1000                99847
     60    3.223e-09    0.000219883315    1000                68218
     70    4.629e-09     0.00022567806    1000                48742
     80    9.467e-09    0.000350877193    1000                37050
     90    4.857e-09    0.000139387392    1000                28697
    100            0                 0    1000                22835

NxN grids and Failure Probability estimates between farthest diagonal pairs (Benchmark - Within 7.8 hours - Monte Carlo Simulation - Edge Failure Probability of 1%)

Below we show mean and variance for s-t network failure probability estimates using Monte Carlo Simulation. Note that this estimates have no guarantees of quality of approximation.

readtable('grid_table_fp_mcs_28080.txt')

ans =

     N      Variance        Mean       CPU_Time_secs    n_iterations
    ___    __________    __________    _____________    ____________

      3    1.6709e-13    0.00020819    28080            1.2458e+09
      4    3.3459e-13    0.00020442    28080            6.1082e+08
      5    5.5347e-13    0.00020358    28080            3.6775e+08
     10    2.5068e-12    0.00020405    28080            8.1385e+07
     20    1.0696e-11    0.00020239    28080            1.8919e+07
     30    2.3987e-11    0.00019488    28080            8.1229e+06
     40    4.4119e-11    0.00019784    28080            4.4834e+06
     50    7.2645e-11    0.00020493    28080            2.8204e+06
     60    1.0597e-10    0.00020435    28080            1.9281e+06
     70      1.45e-10     0.0002006    28080            1.3859e+06
     80      2.07e-10    0.00020815    28080            1.0041e+06
     90      2.56e-10    0.00020251    28080            7.9008e+05
    100      2.93e-10    0.00019228    28080            6.5528e+05

TUTORIAL: grid graph of size 5x5

For the case of a 5x5 grid we can create its adjacency as follows

N = 5; % 25 nodes + 40 edges

A = kron(spdiags(ones(N,2),[1 -1],N,N),...

eye(N))+kron(eye(N),spdiags(ones(N,2),[1 -1],N,N)); % Adjacency

[X,Y] = meshgrid(linspace(0,1,N),linspace(0,1,N)); % Coordinates

Node Failure probabilities (100% Reliable - Uncomment for different cases)

format short

P = zeros( size(A) );

% P ( 1 : length(A)+1 : length(A)^2 ) = rand(length(A),1); %Node failure probability

% P ( 1 : length(A)+1 : length(A)^2 ) = .10; %Node failure probability

P % Failure probability matrix (diagonal entries for nodes)

P = 25x25 double

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Edge failure probabilities (90% Reliable - Uncomment for different cases)

id = find(A);

P(id) = .1; % Same failure probability

%P(id) = rand(size(id)); %Random failure probabilities

P % Failure probability matrix (off diagonal entries for edges)

P = 25x25 double

0 0.1000 0 0 0 0.1000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1000 0 0.1000 0 0 0 0.1000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1000 0 0.1000 0 0 0 0.1000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1000 0 0.1000 0 0 0 0.1000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1000 0 0 0 0 0 0.1000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1000 0 0 0 0 0 0.1000 0 0 0 0.1000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1000 0 0 0 0.1000 0 0.1000 0 0 0 0.1000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1000 0 0 0 0.1000 0 0.1000 0 0 0 0.1000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1000 0 0 0 0.1000 0 0.1000 0 0 0 0.1000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1000 0 0 0 0.1000 0 0 0 0 0 0.1000 0 0 0 0 0 0 0 0 0 0

Plotting (Needs Matlab R2016a)

% Convert edges to nodes for unreliable edges

[Aout,~,~,XYout] = SSP_edge_splitting(A,[X(:),Y(:)]);

id = (1:size(Aout,1));

N1 = size(A,1);

figure

gp = plot(graph(Aout),'XData',XYout(:,1),'YData',XYout(:,2),'NodeLabel',id,...

'MarkerSize',(id<=N1)*4+eps);

% gp = plot(graph(Aout),'Layout','force','NodeLabel',id,...

% 'MarkerSize',(id<=N1)*4+eps);

title('Labeling')

axis off

Source and terminal nodes

s = 1;

t = 25;

Decomposition algorithm (default is path-based Selective RDA)

Note that below we run a few iterations for clarity but we can run this algorithms indefinitely until convergence, truncate until treshold gap, or timeout based on a treshold time.

Selective = 1; %

Recursive = 0; % DFS (1) or BFS(0)

algo = 'minpath'; % 'minpath', 'mincut', 'flow'

plotit = [1 3]; %0: None, 1: Bounds, 2:Bounds Simple, 3: Decomposition Tree, 4: Animated Tree, 5: Animated Bounds

MRP = 1; % Minimum cost shortest path (1) Dijkstra (0)

Max_k = 1e5; %Truncation by iterations

Max_gap = 0.05; %Truncation of bounds (e.g. 0 for exact)

eq = 0; % Output equation

[ R, Max_k, R_bounds, Eventk, Likelihoodk, isConnQ, minpathk, SetUp,...

SetDown, Rf_links, Rf_cuts, parentk, Prob, P] = SSP( A , P , s , t, 2, algo, Selective , Recursive , plotit , Max_k , Max_gap , MRP, eq );

Elapsed time is 0.057093 seconds.

Visualize decomposition process (Table iteration-by-iteration)

T = decomp_table(R_bounds,Eventk,minpathk,isConnQ,parentk,Prob, SetUp, SetDown, P)

Note that nodes are always assumed to be not failed!!

T =

    IterationQueue                                   Event                                                                                       MinimumSuccessfulEvent                                                                                                                               UnreliableElements                                                                                 ChildEventsSentToQueue         LinkSetQ    number_of_config    Likelihood    R_ubound    R_lbound
    ______________    ____________________________________________________________________    ____________________________________________________________________________________________________________________________    ___________________________________________________________________________________________________________________________________________________________________    _______________________________    ________    ________________    __________    ________    ________

      0               'Not Failed: (); Failed: ()'                                            'Needed to be not failed: 1   2   3   4   9  10  15  20  25  26  27  28  35  38  46  55  64'                                    'Non-fixed: 29  30  31  32  33  34  36  37  39  40  41  42  43  44  45  47  48  49  50  51  52  53  54  56  57  58  59  60  61  62  63  65'                            '1  2  3  4  5  6  7  8'           true         4.295e+09             0.43047    0.43047           1
      1               'Not Failed: (); Failed: (26)'                                          'Needed to be not failed: 1   6   7  12  17  22  23  24  25  30  32  40  49  58  61  63  65'                                    'Non-fixed: 27  28  29  31  33  34  35  36  37  38  39  41  42  43  44  45  46  47  48  50  51  52  53  54  55  56  57  59  60  62  64'                                '9  10  11  12  13  14  15  16'    true        2.1475e+09            0.043047    0.47351           1
      2               'Not Failed: (26); Failed: (27)'                                        'Needed to be not failed: 1   2   7  12  13  18  19  20  25  26  31  40  43  51  54  56  64'                                    'Non-fixed: 28  29  30  32  33  34  35  36  37  38  39  41  42  44  45  46  47  48  49  50  52  53  55  57  58  59  60  61  62  63  65'                                '17  18  19  20  21  22  23'       true        2.1475e+09            0.043047    0.51656           1
      3               'Not Failed: (26  27); Failed: (28)'                                    'Needed to be not failed: 1   2   3   8   9  14  15  20  25  26  27  33  36  44  47  55  64'                                    'Non-fixed: 29  30  31  32  34  35  37  38  39  40  41  42  43  45  46  48  49  50  51  52  53  54  56  57  58  59  60  61  62  63  65'                                '24  25  26  27  28  29'           true        2.1475e+09            0.043047    0.55961           1
      4               'Not Failed: (26  27  28); Failed: (35)'                                'Needed to be not failed: 1   2   3   4   5  10  15  20  25  26  27  28  29  37  46  55  64'                                    'Non-fixed: 30  31  32  33  34  36  38  39  40  41  42  43  44  45  47  48  49  50  51  52  53  54  56  57  58  59  60  61  62  63  65'                                '30  31  32  33  34'               true        2.1475e+09            0.043047    0.60265           1
      5               'Not Failed: (26  27  28  35); Failed: (38)'                            'Needed to be not failed: 1   2   3   4   9  14  19  24  25  26  27  28  35  44  53  62  65'                                    'Non-fixed: 29  30  31  32  33  34  36  37  39  40  41  42  43  45  46  47  48  49  50  51  52  54  55  56  57  58  59  60  61  63  64'                                '35  36  37  38'                   true        2.1475e+09            0.043047     0.6457           1
      6               'Not Failed: (26  27  28  35  38); Failed: (46)'                        'Needed to be not failed: 1   2   3   4   9  14  19  24  25  26  27  28  35  44  53  62  65'                                    'Non-fixed: 29  30  31  32  33  34  36  37  39  40  41  42  43  45  47  48  49  50  51  52  54  55  56  57  58  59  60  61  63  64'                                    '39  40  41  42'                   true        1.0737e+09            0.038742    0.68444           1
      7               'Not Failed: (26  27  28  35  38  46); Failed: (55)'                    'Needed to be not failed: 1   2   3   4   9  10  14  15  19  20  25  26  27  28  35  38  46  47  53  56  64'                    'Non-fixed: 29  30  31  32  33  34  36  37  39  40  41  42  43  44  45  48  49  50  51  52  54  57  58  59  60  61  62  63  65'                                        '43  44  45  46'                   true        5.3687e+08            0.034868    0.71931           1
      8               'Not Failed: (26  27  28  35  38  46  55); Failed: (64)'                'Needed to be not failed: 1   2   3   4   9  10  15  19  20  24  25  26  27  28  35  38  46  55  56  62  65'                    'Non-fixed: 29  30  31  32  33  34  36  37  39  40  41  42  43  44  45  47  48  49  50  51  52  53  54  57  58  59  60  61  63'                                        '47  48  49'                       true        5.3687e+08            0.034868    0.75418           1
      9               'Not Failed: (); Failed: (26  30)'                                      'Needed to be not failed: '                                                                                                     'Non-fixed: 27  28  29  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65'    ''                                 false       2.7488e+11                0.01    0.75418        0.99
     10               'Not Failed: (30); Failed: (26  32)'                                    'Needed to be not failed: 1   6  11  12  13  14  15  20  25  30  39  41  43  45  47  55  64'                                    'Non-fixed: 27  28  29  31  33  34  35  36  37  38  40  42  44  46  48  49  50  51  52  53  54  56  57  58  59  60  61  62  63  65'                                    '50  51  52  53  54  55  56'       true        1.0737e+09           0.0043047    0.75848        0.99
     11               'Not Failed: (30  32); Failed: (26  40)'                                'Needed to be not failed: 1   6   7   8   9  14  15  20  25  30  32  34  36  44  47  55  64'                                    'Non-fixed: 27  28  29  31  33  35  37  38  39  41  42  43  45  46  48  49  50  51  52  53  54  56  57  58  59  60  61  62  63  65'                                    '64  65  66  67  68  69'           true        1.0737e+09           0.0043047    0.76279        0.99
     12               'Not Failed: (30  32  40); Failed: (26  49)'                            'Needed to be not failed: 1   6   7  12  13  18  19  24  25  30  32  40  43  51  54  62  65'                                    'Non-fixed: 27  28  29  31  33  34  35  36  37  38  39  41  42  44  45  46  47  48  50  52  53  55  56  57  58  59  60  61  63  64'                                    '82  83  84  85  86'               true        1.0737e+09           0.0043047    0.76709        0.99
     13               'Not Failed: (30  32  40  49); Failed: (26  58)'                        'Needed to be not failed: 1   6   7  12  17  18  19  20  25  30  32  40  49  52  54  56  64'                                    'Non-fixed: 27  28  29  31  33  34  35  36  37  38  39  41  42  43  44  45  46  47  48  50  51  53  55  57  59  60  61  62  63  65'                                    '101  102  103'                    true        1.0737e+09           0.0043047     0.7714        0.99
     14               'Not Failed: (30  32  40  49  58); Failed: (26  61)'                    'Needed to be not failed: 1   6   7  12  17  18  23  24  25  30  32  40  49  52  60  63  65'                                    'Non-fixed: 27  28  29  31  33  34  35  36  37  38  39  41  42  43  44  45  46  47  48  50  51  53  54  55  56  57  59  62  64'                                        '116  117'                         true        5.3687e+08           0.0038742    0.77527        0.99
     15               'Not Failed: (30  32  40  49  58  61); Failed: (26  63)'                'Needed to be not failed: 1   6   7  12  17  18  19  22  23  24  25  30  32  40  49  54  58  60  61  62  65'                    'Non-fixed: 27  28  29  31  33  34  35  36  37  38  39  41  42  43  44  45  46  47  48  50  51  52  53  55  56  57  59  64'                                            '128'                              true        2.6844e+08           0.0034868    0.77876        0.99
     16               'Not Failed: (30  32  40  49  58  61  63); Failed: (26  65)'            'Needed to be not failed: 1   6   7  12  17  19  20  22  23  24  25  30  32  40  49  56  58  61  62  63  64'                    'Non-fixed: 27  28  29  31  33  34  35  36  37  38  39  41  42  43  44  45  46  47  48  50  51  52  53  54  55  57  59  60'                                            ''                                 true        2.6844e+08           0.0034868    0.78225        0.99

% Write excel table

try

delete('T.xls') % Because we rewrite potentially leaving previous results

writetable(T,'T.xls')

display('T.xls file created')

catch

display('Not possible to create T.xls')

end

T.xls file created

Some remarks from table T:

The first iteration (0) is the root event and there are no conditional events in the column Event.
Child events (disjoint) permute edge failures form parent's MinimumSuccessfulEvent and are created lexicographically. Note that nodes do not fail in this example however.
true values in the LinkSetsQ column are internal nodes (with childrens) and false LinkSetQ values are leaf nodes in the decomposition process (no childrens sent to iteration queue).