1 // -*- C++ -*-
2 /**
3 * @brief The test file for Waarts discontinuous limit state function
4 *
5 * Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca
6 *
7 * This library is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with this library. If not, see <http://www.gnu.org/licenses/>.
19 *
20 */
21 #include "openturns/OT.hxx"
22 #include "openturns/OTtestcode.hxx"
23
24 using namespace OT;
25 using namespace OT::Test;
26
clean(Scalar in)27 Scalar clean(Scalar in)
28 {
29 if (std::abs(in) < 1.e-10) in = 0.0;
30 return in;
31 }
32
main(int,char * [])33 int main(int, char *[])
34 {
35 TESTPREAMBLE;
36 OStream fullprint(std::cout);
37 setRandomGenerator();
38 try
39 {
40
41 Description inputFunction(2);
42 inputFunction[0] = "x1";
43 inputFunction[1] = "x2";
44
45 Description outputFunction(1);
46 outputFunction[0] = "G";
47
48 Description formulas(outputFunction.getSize());
49 formulas[0] = "if( x2 <= x1,-0.5+sqrt(x1-x2),-0.5 )";
50
51 Function EtatLimite(inputFunction, outputFunction, formulas);
52
53 UnsignedInteger dim = EtatLimite.getInputDimension();
54 fullprint << dim << std::endl;
55
56 // #########################################################################################################;
57 // # Probabilistic model;
58 // #########################################################################################################;
59
60 Point mean(dim, 0.0);
61 mean[0] = 15.0;
62 mean[1] = 5.0;
63
64 Point sigma(dim, 0.0);
65 sigma[0] = 2.5;
66 sigma[1] = 0.5;
67
68 CorrelationMatrix CORR(dim);
69 Normal myDistribution(mean, sigma, CORR);
70
71 Point start(myDistribution.getMean());
72 CovarianceMatrix Covariance = myDistribution.getCovariance();
73
74 // #########################################################################################################;
75 // # Limit state;
76 // #########################################################################################################;
77
78 RandomVector vect(myDistribution);
79
80 CompositeRandomVector output(EtatLimite, vect);
81
82 ThresholdEvent myEvent(output, Less(), 0.0);
83
84 // #########################################################################################################;
85 // # Calculs;
86 // #########################################################################################################;
87 // ;
88 // #########################################################################################################;
89 // # FORM/SORM Cobyla;
90 Cobyla myCobyla;
91 myCobyla.setMaximumEvaluationNumber(1000 * dim);
92 myCobyla.setMaximumAbsoluteError(1.0e-10);
93 myCobyla.setMaximumRelativeError(1.0e-10);
94 myCobyla.setMaximumResidualError(1.0e-10);
95 myCobyla.setMaximumConstraintError(1.0e-10);
96
97 FORM myAlgoC(myCobyla, myEvent, start);
98 SORM myAlgoC2(myCobyla, myEvent, start);
99
100 myAlgoC.run();
101 myAlgoC2.run();
102
103 FORMResult resultC(myAlgoC.getResult());
104 SORMResult resultC2(myAlgoC2.getResult());
105
106 // #########################################################################################################;
107 // # FORM/SORM Abdo Rackwitz;
108 AbdoRackwitz myAbdoRackwitz;
109 myAbdoRackwitz.setMaximumIterationNumber(100 * dim);
110 myAbdoRackwitz.setMaximumAbsoluteError(1.0e-10);
111 myAbdoRackwitz.setMaximumRelativeError(1.0e-10);
112 myAbdoRackwitz.setMaximumResidualError(1.0e-10);
113 myAbdoRackwitz.setMaximumConstraintError(1.0e-10);
114
115 FORM myAlgoAR(myAbdoRackwitz, myEvent, start);
116 SORM myAlgoAR2(myAbdoRackwitz, myEvent, start);
117
118 myAlgoAR.run();
119 myAlgoAR2.run();
120
121 FORMResult resultAR(myAlgoAR.getResult());
122 SORMResult resultAR2(myAlgoAR2.getResult());
123
124 // #########################################################################################################;
125 // # Monte Carlo;
126 Scalar CoV_MC = 0.5;
127 MonteCarlo myMC(myEvent);
128 myMC.setMaximumOuterSampling(10000000);
129 myMC.setBlockSize(1000);
130 myMC.setMaximumCoefficientOfVariation(CoV_MC);
131 myMC.run();
132
133 // #########################################################################################################;
134 // # LHS;
135 Scalar CoV_LHS = 0.1;
136 LHS myLHS(myEvent);
137 myLHS.setMaximumOuterSampling(100000);
138 myLHS.setBlockSize(100);
139 myLHS.setMaximumCoefficientOfVariation(CoV_LHS);
140 myLHS.run();
141 /* ;
142
143 #########################################################################################################;
144
145 #########################################################################################################;
146 # Results;
147 #########################################################################################################;
148
149 #########################################################################################################;
150 # FORM/SORM Cobyla*/;
151 Scalar PfC = resultC.getEventProbability();
152 Scalar Beta_generalizedC = resultC.getGeneralisedReliabilityIndex();
153 Point u_starC = resultC.getStandardSpaceDesignPoint();
154 Point x_starC = resultC.getPhysicalSpaceDesignPoint();
155 Bool PtC = resultC.getIsStandardPointOriginInFailureSpace();
156 Point gammaC = resultC.getImportanceFactors();
157 Point gammaCC = resultC.getImportanceFactors(true);
158 Scalar beta_hasoferC = resultC.getHasoferReliabilityIndex();
159 Analytical::Sensitivity SensitivityC = resultC.getEventProbabilitySensitivity();
160
161 Scalar PFBreitC2 = resultC2.getEventProbabilityBreitung();
162 Scalar BetaBreitC2 = resultC2.getGeneralisedReliabilityIndexBreitung();
163 Scalar PFHBC2 = resultC2.getEventProbabilityHohenbichler();
164 Scalar BetaHBC2 = resultC2.getGeneralisedReliabilityIndexHohenbichler();
165 Scalar PFTvedtC2 = resultC2.getEventProbabilityTvedt();
166 Scalar BetaTvedtC2 = resultC2.getGeneralisedReliabilityIndexTvedt();
167 Point CurvC2 = resultC2.getSortedCurvatures();
168 Point u_starC2 = resultC2.getStandardSpaceDesignPoint();
169 Point x_starC2 = resultC2.getPhysicalSpaceDesignPoint();
170 Bool PtC2 = resultC2.getIsStandardPointOriginInFailureSpace();
171 Point gammaC2 = resultC2.getImportanceFactors();
172 Point gammaCC2 = resultC2.getImportanceFactors(true);
173 Scalar beta_hasoferC2 = resultC2.getHasoferReliabilityIndex();
174
175 // #########################################################################################################;
176 // # FORM/SORM Abdo Rackwitz;
177 Scalar PfAR = resultAR.getEventProbability();
178 Scalar Beta_generalizedAR = resultAR.getGeneralisedReliabilityIndex();
179 Point u_starAR = resultAR.getStandardSpaceDesignPoint();
180 Point x_starAR = resultAR.getPhysicalSpaceDesignPoint();
181 Bool PtAR = resultAR.getIsStandardPointOriginInFailureSpace();
182 Point gammaAR = resultAR.getImportanceFactors();
183 Point gammaCAR = resultAR.getImportanceFactors(true);
184 Scalar beta_hasoferAR = resultAR.getHasoferReliabilityIndex();
185 Analytical::Sensitivity SensitivityAR = resultAR.getEventProbabilitySensitivity();
186
187 Scalar PFBreitAR2 = resultAR2.getEventProbabilityBreitung();
188 Scalar BetaBreitAR2 = resultAR2.getGeneralisedReliabilityIndexBreitung();
189 Scalar PFHBAR2 = resultAR2.getEventProbabilityHohenbichler();
190 Scalar BetaHBAR2 = resultAR2.getGeneralisedReliabilityIndexHohenbichler();
191 Scalar PFTvedtAR2 = resultAR2.getEventProbabilityTvedt();
192 Scalar BetaTvedtAR2 = resultAR2.getGeneralisedReliabilityIndexTvedt();
193 Point CurvAR2 = resultAR2.getSortedCurvatures();
194 Point u_starAR2 = resultAR2.getStandardSpaceDesignPoint();
195 Point x_starAR2 = resultAR2.getPhysicalSpaceDesignPoint();
196 Bool PtAR2 = resultAR2.getIsStandardPointOriginInFailureSpace();
197 Point gammaAR2 = resultAR2.getImportanceFactors();
198 Point gammaCAR2 = resultAR2.getImportanceFactors(true);
199 Scalar beta_hasoferAR2 = resultAR2.getHasoferReliabilityIndex();
200
201 // ######################################/*###################################################################;
202 // # Monte Carlo*/;
203 SimulationResult ResultMC = myMC.getResult();
204 Scalar PFMC = ResultMC.getProbabilityEstimate();
205 Scalar CVMC = ResultMC.getCoefficientOfVariation();
206 Scalar Variance_PF_MC = ResultMC.getVarianceEstimate();
207 Scalar length90MC = ResultMC.getConfidenceLength(0.90);
208
209 // #########################################################################################################;
210 // # LHS;
211 SimulationResult ResultLHS = myLHS.getResult();
212 Scalar PFLHS = ResultLHS.getProbabilityEstimate();
213 Scalar CVLHS = ResultLHS.getCoefficientOfVariation();
214 Scalar Variance_PF_LHS = ResultLHS.getVarianceEstimate();
215 Scalar length90LHS = ResultLHS.getConfidenceLength(0.90);
216
217 // #########################################################################################################;
218 // # Printting;
219 // #########################################################################################################;
220 fullprint << "" << std::endl;
221 fullprint << "" << std::endl;
222 fullprint << "************************************************************************************************" << std::endl;
223 fullprint << "***************************************** FORM COBYLA *****************************************" << std::endl;
224 fullprint << "************************************************************************************************" << std::endl;
225 fullprint << "event probability =" << PfC << std::endl;
226 fullprint << "generalized reliability index =" << Beta_generalizedC << std::endl;
227 fullprint << "************************************************************************************************" << std::endl;
228 for (UnsignedInteger i = 0; i < u_starC.getDimension(); i++)
229 fullprint << "standard space design point =" << u_starC[i] << std::endl;
230 fullprint << "************************************************************************************************" << std::endl;
231 for (UnsignedInteger i = 0; i < x_starC.getDimension(); i++)
232 fullprint << "physical space design point =" << x_starC[i] << std::endl;
233 fullprint << "************************************************************************************************" << std::endl;
234 fullprint << "is standard point origin in failure space? " << PtC << std::endl;
235 fullprint << "************************************************************************************************" << std::endl;
236 for (UnsignedInteger i = 0; i < gammaC.getDimension(); i++)
237 fullprint << "importance factors =" << gammaC[i] << std::endl;
238 for (UnsignedInteger i = 0; i < gammaCC.getDimension(); i++)
239 fullprint << "importance factors (classical)=" << gammaCC[i] << std::endl;
240 fullprint << "************************************************************************************************" << std::endl;
241 fullprint << "Hasofer reliability index =" << beta_hasoferC << std::endl;
242 fullprint << "************************************************************************************************" << std::endl;
243 for (UnsignedInteger i = 0; i < SensitivityAR.getSize(); ++i)
244 for (UnsignedInteger j = 0; j < SensitivityAR[i].getDimension(); ++j)
245 fullprint << "Pf sensitivity =" << i << j << SensitivityC[i][j] << std::endl;
246 fullprint << "************************************************************************************************" << std::endl;
247 fullprint << "" << std::endl;
248 fullprint << "************************************************************************************************" << std::endl;
249 fullprint << "************************************** FORM ABDO RACKWITZ **************************************" << std::endl;
250 fullprint << "************************************************************************************************" << std::endl;
251 fullprint << "event probability =" << PfAR << std::endl;
252 fullprint << "generalized reliability index =" << Beta_generalizedAR << std::endl;
253 fullprint << "************************************************************************************************" << std::endl;
254 for (UnsignedInteger i = 0; i < u_starAR.getDimension(); i++)
255 fullprint << "standard space design point =" << u_starAR[i] << std::endl;
256 fullprint << "************************************************************************************************" << std::endl;
257 for (UnsignedInteger i = 0; i < x_starAR.getDimension(); i++)
258 fullprint << "physical space design point =" << x_starAR[i] << std::endl;
259 fullprint << "************************************************************************************************" << std::endl;
260 fullprint << "is standard point origin in failure space? " << PtAR << std::endl;
261 fullprint << "************************************************************************************************" << std::endl;
262 for (UnsignedInteger i = 0; i < gammaAR.getDimension(); i++)
263 fullprint << "importance factors =" << gammaAR[i] << std::endl;
264 for (UnsignedInteger i = 0; i < gammaCAR.getDimension(); i++)
265 fullprint << "importance factors (classical)=" << gammaCAR[i] << std::endl;
266 fullprint << "************************************************************************************************" << std::endl;
267 fullprint << "Hasofer reliability index =" << beta_hasoferAR << std::endl;
268 fullprint << "************************************************************************************************" << std::endl;
269 for (UnsignedInteger i = 0; i < SensitivityAR.getSize(); ++i)
270 for (UnsignedInteger j = 0; j < SensitivityAR[i].getDimension(); ++j)
271 fullprint << "Pf sensitivity =" << i << j << SensitivityAR[i][j] << std::endl;
272 fullprint << "************************************************************************************************" << std::endl;
273 fullprint << "" << std::endl;
274 fullprint << "************************************************************************************************" << std::endl;
275 fullprint << "***************************************** SORM COBYLA *****************************************" << std::endl;
276 fullprint << "************************************************************************************************" << std::endl;
277 fullprint << "Breitung event probability =" << PFBreitC2 << std::endl;
278 fullprint << "Breitung generalized reliability index =" << BetaBreitC2 << std::endl;
279 fullprint << "Hohenbichler event probability =" << PFHBC2 << std::endl;
280 fullprint << "Hohenbichler generalized reliability index =" << BetaHBC2 << std::endl;
281 fullprint << "Tvedt event probability =" << PFTvedtC2 << std::endl;
282 fullprint << "Tvedt generalized reliability index =" << BetaTvedtC2 << std::endl;
283 fullprint << "************************************************************************************************" << std::endl;
284 for (UnsignedInteger i = 0; i < CurvC2.getDimension(); i++)
285 fullprint << "sorted curvatures =" << clean(CurvC2[i]) << std::endl;
286 fullprint << "************************************************************************************************" << std::endl;
287 for (UnsignedInteger i = 0; i < u_starC2.getDimension(); i++)
288 fullprint << "standard space design point =" << u_starC2[i] << std::endl;
289 fullprint << "************************************************************************************************" << std::endl;
290 for (UnsignedInteger i = 0; i < x_starC2.getDimension(); i++)
291 fullprint << "physical space design point =" << x_starC2[i] << std::endl;
292 fullprint << "************************************************************************************************" << std::endl;
293 fullprint << "************************************************************************************************" << std::endl;
294 fullprint << "is standard point origin in failure space? " << PtC2 << std::endl;
295 fullprint << "************************************************************************************************" << std::endl;
296 for (UnsignedInteger i = 0; i < gammaC2.getDimension(); i++)
297 fullprint << "importance factors =" << gammaC2[i] << std::endl;
298 for (UnsignedInteger i = 0; i < gammaCC2.getDimension(); i++)
299 fullprint << "importance factors (classical)=" << gammaCC2[i] << std::endl;
300 fullprint << "************************************************************************************************" << std::endl;
301 fullprint << "Hasofer reliability index =" << beta_hasoferC2 << std::endl;
302 fullprint << "************************************************************************************************" << std::endl;
303 fullprint << "" << std::endl;
304 fullprint << "************************************************************************************************" << std::endl;
305 fullprint << "************************************** SORM ABDO RACKWITZ **************************************" << std::endl;
306 fullprint << "************************************************************************************************" << std::endl;
307 fullprint << "Breitung event probability =" << PFBreitAR2 << std::endl;
308 fullprint << "Breitung generalized reliability index =" << BetaBreitAR2 << std::endl;
309 fullprint << "Hohenbichler event probability =" << PFHBAR2 << std::endl;
310 fullprint << "Hohenbichler generalized reliability index =" << BetaHBAR2 << std::endl;
311 fullprint << "Tvedt event probability =" << PFTvedtAR2 << std::endl;
312 fullprint << "Tvedt generalized reliability index =" << BetaTvedtAR2 << std::endl;
313 fullprint << "************************************************************************************************" << std::endl;
314 for (UnsignedInteger i = 0; i < CurvAR2.getDimension(); i++)
315 fullprint << "sorted curvatures =" << clean(CurvAR2[i]) << std::endl;
316 fullprint << "************************************************************************************************" << std::endl;
317 for (UnsignedInteger i = 0; i < u_starAR2.getDimension(); i++)
318 fullprint << "standard space design point =" << u_starAR2[i] << std::endl;
319 fullprint << "************************************************************************************************" << std::endl;
320 for (UnsignedInteger i = 0; i < x_starAR2.getDimension(); i++)
321 fullprint << "physical space design point =" << x_starAR2[i] << std::endl;
322 fullprint << "************************************************************************************************" << std::endl;
323 fullprint << "************************************************************************************************" << std::endl;
324 fullprint << "is standard point origin in failure space? " << PtAR2 << std::endl;
325 fullprint << "************************************************************************************************" << std::endl;
326 for (UnsignedInteger i = 0; i < gammaAR2.getDimension(); i++)
327 fullprint << "importance factors =" << gammaAR2[i] << std::endl;
328 for (UnsignedInteger i = 0; i < gammaCAR2.getDimension(); i++)
329 fullprint << "importance factors (classical)=" << gammaCAR2[i] << std::endl;
330 fullprint << "************************************************************************************************" << std::endl;
331 fullprint << "Hasofer reliability index =" << beta_hasoferAR2 << std::endl;
332 fullprint << "************************************************************************************************" << std::endl;
333 fullprint << "" << std::endl;
334 fullprint << "************************************************************************************************" << std::endl;
335 fullprint << "**************************************** MONTE CARLO *******************************************" << std::endl;
336 fullprint << "************************************************************************************************" << std::endl;
337 fullprint << "Pf estimation =" << PFMC << std::endl;
338 fullprint << "Pf Variance estimation =" << Variance_PF_MC << std::endl;
339 fullprint << "CoV =" << CVMC << std::endl;
340 fullprint << "90% Confidence Interval =" << length90MC << std::endl;
341 fullprint << "CI at 90% =[" << PFMC - 0.5 * length90MC << ";" << PFMC + 0.5 * length90MC << "]" << std::endl;
342 fullprint << "************************************************************************************************" << std::endl;
343 fullprint << "" << std::endl;
344 fullprint << "************************************************************************************************" << std::endl;
345 fullprint << "******************************************* L H S **********************************************" << std::endl;
346 fullprint << "************************************************************************************************" << std::endl;
347 fullprint << "Pf estimation =" << PFLHS << std::endl;
348 fullprint << "Pf Variance estimation =" << Variance_PF_LHS << std::endl;
349 fullprint << "CoV =" << CVLHS << std::endl;
350 fullprint << "90% Confidence Interval =" << length90LHS << std::endl;
351 fullprint << "CI at 90% =[" << PFLHS - 0.5 * length90LHS << ";" << PFLHS + 0.5 * length90LHS << "]" << std::endl;
352 fullprint << "************************************************************************************************" << std::endl;
353 }
354 catch (TestFailed & ex)
355 {
356 std::cerr << ex << std::endl;
357 return ExitCode::Error;
358 }
359
360 return ExitCode::Success;
361
362 }
363