333. FORM - Structural Reliability
44==============================================
55
6- The benchmark problem is a simple structural reliability problem
6+ The benchmark problem is a simple structural reliability problem (example 7.1 in :cite:`FORM_XDu`)
77defined in a two-dimensional parameter space consisting of a resistance :math:`R` and a stress :math:`S`. The failure
88happens when the stress is higher than the resistance, leading to the following limit-state function:
99
10- .. math:: \t extbf{X}=\{R, S\}
10+ .. math:: \\ textbf{X}=\{R, S\}
1111
12- .. math:: g(\t extbf{X}) = R - S
12+ .. math:: g(\\ textbf{X}) = R - S
1313
1414The two random variables are independent and distributed
1515according to:
1919.. math:: S \sim N(150, 10)
2020"""
2121
22- #%% md
22+ # %% md
2323#
2424# Initially we have to import the necessary modules.
2525
26- #%%
27- import shutil
26+ # %%
2827
2928import numpy as np
3029import matplotlib .pyplot as plt
31- from UQpy .run_model .RunModel import RunModel
32- from UQpy .run_model .model_execution .PythonModel import PythonModel
30+ plt .style .use ('ggplot' )
3331from UQpy .distributions import Normal
3432from UQpy .reliability import FORM
33+ from UQpy .run_model .RunModel import RunModel
34+ from UQpy .run_model .model_execution .PythonModel import PythonModel
35+
36+
37+ # %% md
38+ #
39+ # Next, we initialize the :code:`RunModel` object.
40+ # The `local_pfn.py <https://github.com/SURGroup/UQpy/tree/master/docs/code/reliability/sorm>`_ file can be found on
41+ # the UQpy GitHub. It contains a simple function :code:`example1` to compute the difference between the resistence and the
42+ # stress.
43+
44+ # %%
45+
46+ model = PythonModel (model_script = 'local_pfn.py' , model_object_name = "example1" )
47+ runmodel_object = RunModel (model = model )
48+
49+ # %% md
50+ #
51+ # Now we can define the resistence and stress distributions that will be passed into :code:`FORM`.
52+ # Along with the distributions, :code:`FORM` takes in the previously defined :code:`runmodel_object` and tolerances
53+ # for convergences. Since :code:`tolerance_gradient` is not specified in this example, it is not considered.
3554
55+ # %%
3656
37- model = PythonModel (model_script = 'pfn.py' , model_object_name = "example1" )
38- RunModelObject = RunModel (model = model )
57+ distribution_resistance = Normal (loc = 200. , scale = 20. )
58+ distribution_stress = Normal (loc = 150. , scale = 10. )
59+ form = FORM (distributions = [distribution_resistance , distribution_stress ], runmodel_object = runmodel_object ,
60+ tolerance_u = 1e-5 , tolerance_beta = 1e-5 )
61+ # %% md
62+ #
63+ # With everything defined we are ready to run the first-order reliability method and print the results.
64+ # The analytic solution to this problem is :math:`\textbf{u}^*=(-2, 1)` with a reliability index of
65+ # :math:`\beta_{HL}=2.2361` and a probability of failure :math:`P_{f, \text{form}} = \Phi(-\beta_{HL}) = 0.0127`
3966
40- dist1 = Normal (loc = 200. , scale = 20. )
41- dist2 = Normal (loc = 150 , scale = 10. )
42- Q = FORM (distributions = [dist1 , dist2 ], runmodel_object = RunModelObject , tol1 = 1e-5 , tol2 = 1e-5 )
43- Q .run ()
67+ # %%
4468
69+ form .run ()
70+ print ('Design point in standard normal space:' , form .design_point_u )
71+ print ('Design point in original space:' , form .design_point_x )
72+ print ('Hasofer-Lind reliability index:' , form .beta )
73+ print ('FORM probability of failure:' , form .failure_probability )
74+ print ('FORM record of the function gradient:' , form .state_function_gradient_record )
4575
46- # print results
47- print ('Design point in standard normal space: %s' % Q .DesignPoint_U )
48- print ('Design point in original space: %s' % Q .DesignPoint_X )
49- print ('Hasofer-Lind reliability index: %s' % Q .beta )
50- print ('FORM probability of failure: %s' % Q .failure_probability )
51- print (Q .dg_u_record )
76+ # %% md
77+ #
78+ # This problem can be visualized in the following plots that show the FORM results in both :math:`\textbf{X}` and
79+ # :math:`\textbf{U}` space.
5280
81+ # %%
5382
54- # Supporting function
55- def multivariate_gaussian ( pos , mu , Sigma ):
83+ def multivariate_gaussian ( pos , mu , sigma ):
84+ """Supporting function"""
5685 n = mu .shape [0 ]
57- Sigma_det = np .linalg .det (Sigma )
58- Sigma_inv = np .linalg .inv (Sigma )
59- N = np .sqrt ((2 * np .pi ) ** n * Sigma_det )
60- fac = np .einsum ('...k,kl,...l->...' , pos - mu , Sigma_inv , pos - mu )
86+ sigma_det = np .linalg .det (sigma )
87+ sigma_inv = np .linalg .inv (sigma )
88+ N = np .sqrt ((2 * np .pi ) ** n * sigma_det )
89+ fac = np .einsum ('...k,kl,...l->...' , pos - mu , sigma_inv , pos - mu )
6190 return np .exp (- fac / 2 ) / N
6291
92+
6393N = 60
6494XX = np .linspace (150 , 250 , N )
6595YX = np .linspace (120 , 180 , N )
@@ -69,85 +99,65 @@ def multivariate_gaussian(pos, mu, Sigma):
6999YU = np .linspace (- 3 , 3 , N )
70100XU , YU = np .meshgrid (XU , YU )
71101
72- # Mean vector and covariance matrix in the original space
73- parameters = [[200 , 20 ], [150 , 10 ]]
74- mu_X = np .array ([parameters [0 ][0 ], parameters [1 ][0 ]])
75- Sigma_X = np .array ([[parameters [0 ][1 ] ** 2 , 0.0 ], [0.0 , parameters [1 ][1 ] ** 2 ]])
76102
77- # Mean vector and covariance matrix in the standard normal space
78- mu_U = np .array ([0. , 0. ])
79- Sigma_U = np .array ([[1. , 0.0 ], [0.0 , 1 ]])
103+ # %% md
104+ #
105+ # Define the mean vector and covariance matrix in the original :math:`\textbf{X}` space and the standard normal
106+ # :math:`\textbf{U}` space.
107+
108+ # %%
109+ mu_X = np .array ([distribution_resistance .parameters ['loc' ], distribution_stress .parameters ['loc' ]])
110+ sigma_X = np .array ([[distribution_resistance .parameters ['scale' ]** 2 , 0 ],
111+ [0 , distribution_stress .parameters ['scale' ]** 2 ]])
112+
113+ mu_U = np .array ([0 , 0 ])
114+ sigma_U = np .array ([[1 , 0 ],
115+ [0 , 1 ]])
80116
81117# Pack X and Y into a single 3-dimensional array for the original space
82118posX = np .empty (XX .shape + (2 ,))
83119posX [:, :, 0 ] = XX
84120posX [:, :, 1 ] = YX
85- ZX = multivariate_gaussian (posX , mu_X , Sigma_X )
121+ ZX = multivariate_gaussian (posX , mu_X , sigma_X )
86122
87123# Pack X and Y into a single 3-dimensional array for the standard normal space
88124posU = np .empty (XU .shape + (2 ,))
89125posU [:, :, 0 ] = XU
90126posU [:, :, 1 ] = YU
91- ZU = multivariate_gaussian (posU , mu_U , Sigma_U )
92-
93- # Figure 4a
94- plt .figure ()
95- plt .rcParams ["figure.figsize" ] = (12 , 12 )
96- plt .rcParams .update ({'font.size' : 22 })
97- plt .plot (parameters [0 ][0 ], parameters [1 ][0 ], 'r.' )
98- plt .plot ([0 , 200 ], [0 , 200 ], 'k' , linewidth = 5 )
99- plt .plot (Q .DesignPoint_X [0 ][0 ], Q .DesignPoint_X [0 ][1 ], 'bo' , markersize = 12 )
100- plt .contour (XX , YX , ZX , levels = 20 )
101- plt .xlabel (r'$X_1$' )
102- plt .ylabel (r'$X_2$' )
103- plt .text (170 , 182 , '$X_1 - X_2=0$' ,
104- rotation = 45 ,
105- horizontalalignment = 'center' ,
106- verticalalignment = 'top' ,
107- multialignment = 'center' )
108- plt .ylim ([120 , 200 ])
109- plt .xlim ([130 , 240 ])
110- plt .grid ()
111- plt .title ('Original space' )
112- plt .axes ().set_aspect ('equal' , 'box' )
113- plt .show ()
127+ ZU = multivariate_gaussian (posU , mu_U , sigma_U )
128+
129+ # %% md
130+ #
131+ # Plot the :code:`FORM` solution in the original :math:`\textbf{X}` space and the standard normal :math:`\text{U}`
132+ # space.
133+
134+ # %%
135+ fig , ax = plt .subplots ()
136+ ax .contour (XX , YX , ZX ,
137+ levels = 20 )
138+ ax .plot ([0 , 200 ], [0 , 200 ],
139+ color = 'black' , linewidth = 2 , label = '$G(R,S)=R-S=0$' , zorder = 1 )
140+ ax .scatter (mu_X [0 ], mu_X [1 ],
141+ color = 'black' , s = 64 , label = 'Mean $(\mu_R, \mu_S)$' )
142+ ax .scatter (form .design_point_x [0 ][0 ], form .design_point_x [0 ][1 ],
143+ color = 'tab:orange' , marker = '*' , s = 100 , label = 'Design Point' , zorder = 2 )
144+ ax .set (xlabel = 'Resistence $R$' , ylabel = 'Stress $S$' , xlim = (145 , 255 ), ylim = (115 , 185 ))
145+ ax .set_title ('Original $X$ Space ' )
146+ ax .set_aspect ('equal' )
147+ ax .legend (loc = 'lower right' )
148+
149+ fig , ax = plt .subplots ()
150+ ax .contour (XU , YU , ZU ,
151+ levels = 20 , zorder = 1 )
152+ ax .plot ([0 , - 3 ], [5 , - 1 ],
153+ color = 'black' , linewidth = 2 , label = '$G(U_1, U_2)=0$' , zorder = 2 )
154+ ax .arrow (0 , 0 , form .design_point_u [0 ][0 ], form .design_point_u [0 ][1 ],
155+ color = 'tab:blue' , length_includes_head = True , width = 0.05 , label = '$\\ beta=||u^*||$' , zorder = 2 )
156+ ax .scatter (form .design_point_u [0 ][0 ], form .design_point_u [0 ][1 ],
157+ color = 'tab:orange' , marker = '*' , s = 100 , label = 'Design Point $u^*$' , zorder = 2 )
158+ ax .set (xlabel = '$U_1$' , ylabel = '$U_2$' , xlim = (- 3 , 3 ), ylim = (- 3 , 3 ))
159+ ax .set_aspect ('equal' )
160+ ax .set_title ('Standard Normal $U$ Space' )
161+ ax .legend (loc = 'lower right' )
114162
115- # Figure 4b
116- plt .figure ()
117- plt .rcParams ["figure.figsize" ] = (12 , 12 )
118- plt .rcParams .update ({'font.size' : 22 })
119- plt .plot ([0 , Q .DesignPoint_U [0 ][0 ]], [0 , Q .DesignPoint_U [0 ][1 ]], 'b' , linewidth = 2 )
120- plt .plot ([0 , - 3 ], [5 , - 1 ], 'k' , linewidth = 5 )
121- plt .plot (Q .DesignPoint_U [0 ][0 ], Q .DesignPoint_U [0 ][1 ], 'bo' , markersize = 12 )
122- plt .contour (XU , YU , ZU , levels = 20 )
123- plt .axhline (0 , color = 'black' )
124- plt .axvline (0 , color = 'black' )
125- plt .plot (0 , 0 , 'r.' )
126-
127- plt .xlabel (r'$U_1$' )
128- plt .ylabel (r'$U_2$' )
129- plt .text (- 1.0 , 1.1 , '$U^\star$=({:1.2f}, {:1.2f})' .format (- 2.0 , 1.0 ),
130- rotation = 0 ,
131- horizontalalignment = 'center' ,
132- verticalalignment = 'top' ,
133- multialignment = 'center' )
134-
135- plt .text (- 2.1 , 2.05 , '$20U_1 - 10U_2 + 50=0$' ,
136- rotation = 63 ,
137- horizontalalignment = 'center' ,
138- verticalalignment = 'top' ,
139- multialignment = 'center' )
140-
141- plt .text (- 1.5 , 0.7 , r'$\overrightarrow{\beta}$' ,
142- rotation = 0 ,
143- horizontalalignment = 'center' ,
144- verticalalignment = 'top' ,
145- multialignment = 'center' )
146-
147- plt .text (0.02 , - 0.2 , '({:1.1f}, {:1.1f})' .format (0.0 , 0.0 ))
148- plt .ylim ([- 1 , 3 ])
149- plt .xlim ([- 3.5 , 2 ])
150- plt .grid ()
151- plt .title ('Standard Normal space' )
152- plt .axes ().set_aspect ('equal' , 'box' )
153163plt .show ()
0 commit comments