SystemAnalysis.cppGo to the documentation of this file.00001 /* ****************************************************************** ** 00002 ** OpenSees - Open System for Earthquake Engineering Simulation ** 00003 ** Pacific Earthquake Engineering Research Center ** 00004 ** ** 00005 ** ** 00006 ** (C) Copyright 2001, The Regents of the University of California ** 00007 ** All Rights Reserved. ** 00008 ** ** 00009 ** Commercial use of this program without express permission of the ** 00010 ** University of California, Berkeley, is strictly prohibited. See ** 00011 ** file 'COPYRIGHT' in main directory for information on usage and ** 00012 ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** 00013 ** ** 00014 ** Developed by: ** 00015 ** Frank McKenna (fmckenna@ce.berkeley.edu) ** 00016 ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** 00017 ** Filip C. Filippou (filippou@ce.berkeley.edu) ** 00018 ** ** 00019 ** Reliability module developed by: ** 00020 ** Terje Haukaas (haukaas@ce.berkeley.edu) ** 00021 ** Armen Der Kiureghian (adk@ce.berkeley.edu) ** 00022 ** ** 00023 ** ****************************************************************** */ 00024 00025 // $Revision: 1.6 $ 00026 // $Date: 2003/10/27 23:45:41 $ 00027 // $Source: /usr/local/cvs/OpenSees/SRC/reliability/analysis/analysis/SystemAnalysis.cpp,v $ 00028 00029 00030 // 00031 // Written by Terje Haukaas (haukaas@ce.berkeley.edu) 00032 // 00033 00034 #include <SystemAnalysis.h> 00035 #include <ReliabilityDomain.h> 00036 #include <ReliabilityAnalysis.h> 00037 #include <LimitStateFunction.h> 00038 #include <MatrixOperations.h> 00039 #include <NormalRV.h> 00040 #include <Vector.h> 00041 #include <Matrix.h> 00042 #include <math.h> 00043 #include <string.h> 00044 00045 #include <fstream> 00046 #include <iomanip> 00047 #include <iostream> 00048 using std::ifstream; 00049 using std::ios; 00050 using std::setw; 00051 using std::setprecision; 00052 using std::setiosflags; 00053 00054 SystemAnalysis::SystemAnalysis(ReliabilityDomain *passedReliabilityDomain, 00055 TCL_Char *passedFileName) 00056 :ReliabilityAnalysis() 00057 { 00058 theReliabilityDomain = passedReliabilityDomain; 00059 fileName = new char[256]; 00060 strcpy(fileName,passedFileName); 00061 } 00062 00063 00064 SystemAnalysis::~SystemAnalysis() 00065 { 00066 if (fileName != 0) 00067 delete [] fileName; 00068 } 00069 00070 00071 double 00072 SystemAnalysis::getLowerBound(void) 00073 { 00074 return minLowerBound; 00075 } 00076 00077 00078 00079 double 00080 SystemAnalysis::getUpperBound(void) 00081 { 00082 return maxUpperBound; 00083 } 00084 00085 00086 00087 int 00088 SystemAnalysis::analyze(void) 00089 { 00090 00091 // Alert the user that the system analysis has started 00092 opserr << "System Reliability Analysis is running ... " << endln; 00093 00094 00095 // Initial declarations 00096 double beta; 00097 double pf1; 00098 Vector alpha; 00099 LimitStateFunction *theLimitStateFunction; 00100 NormalRV aStdNormalRV(1, 0.0, 1.0, 0.0); 00101 int i, j, k, m, n; 00102 00103 // Number of limit-state functions 00104 int numLsf = theReliabilityDomain->getNumberOfLimitStateFunctions(); 00105 00106 // Number of random variables 00107 int nrv = theReliabilityDomain->getNumberOfRandomVariables(); 00108 00109 // Allocate vectors to store ALL the betas and alphas 00110 Vector allBetas(numLsf); 00111 Vector allPf1s(numLsf); 00112 Matrix allAlphas(nrv,numLsf); 00113 00114 // Loop over number of limit-state functions and collect results 00115 for (i=0; i<numLsf; i++ ) { 00116 00117 // Get FORM results from the limit-state function 00118 theLimitStateFunction = theReliabilityDomain->getLimitStateFunctionPtr(i+1); 00119 beta = theLimitStateFunction->FORMReliabilityIndexBeta; 00120 pf1 = theLimitStateFunction->FORMProbabilityOfFailure_pf1; 00121 alpha = theLimitStateFunction->normalizedNegativeGradientVectorAlpha; 00122 00123 // Put FORM results into vector of all betas and alphas 00124 // Note that the first index of 'allAlphas' here denote 'nrv' 00125 allBetas(i) = beta; 00126 allPf1s(i) = pf1; 00127 for (j=0; j<nrv; j++ ) { 00128 allAlphas(j,i) = alpha(j); 00129 } 00130 } 00131 00132 00133 // Compute vector of 'rhos', that is, dot products between the alphas 00134 // Note that only the upper diagonal is covered 00135 Matrix rhos(numLsf,numLsf); 00136 double dotProduct; 00137 for (i=0; i<numLsf; i++ ) { 00138 for (j=i+1; j<numLsf; j++ ) { 00139 dotProduct = 1.0; 00140 for (k=0; k<nrv; k++ ) { 00141 dotProduct=dotProduct*allAlphas(k,i)*allAlphas(k,j); 00142 } 00143 rhos(i,j) = dotProduct; 00144 } 00145 } 00146 00147 00148 // Compute the bi-variate normal distribution for all the pairs 00149 double beta1, beta2, integral; 00150 double a = 0.0; 00151 double h, fa, fb, sum_fx2j, sum_fx2j_1; 00152 int n_2; 00153 Matrix Pmn(numLsf,numLsf); 00154 for (m=0; m<numLsf; m++ ) { 00155 for (n=m+1; n<numLsf; n++ ) { 00156 beta1 = allBetas(m); 00157 beta2 = allBetas(n); 00158 // Composite Simpson's numerical integration 00159 n_2 = 400; // (half the number of intervals) 00160 h = rhos(m,n)/(2.0*n_2); 00161 fa = functionToIntegrate(a,beta1,beta2); 00162 fb = functionToIntegrate(rhos(m,n),beta1,beta2); 00163 sum_fx2j = 0.0; 00164 sum_fx2j_1 = 0.0; 00165 for (int j=1; j<=n_2; j++) { 00166 sum_fx2j = sum_fx2j + functionToIntegrate((a+(j*2.0)*h) ,beta1,beta2); 00167 sum_fx2j_1 = sum_fx2j_1 + functionToIntegrate((a+(j*2.0-1.0)*h),beta1,beta2); 00168 } 00169 sum_fx2j = sum_fx2j - fb; 00170 integral = h/3.0*(fa + 2.0*sum_fx2j + 4.0*sum_fx2j_1 + fb); 00171 Pmn(m,n) = aStdNormalRV.getCDFvalue(-beta1) 00172 *aStdNormalRV.getCDFvalue(-beta2) + integral; 00173 } 00174 } 00175 00176 // Arrange-envelope 00177 minLowerBound = 1.0; 00178 maxUpperBound = 0.0; 00179 Vector arrangement(numLsf); 00180 Matrix Pmn_arranged = Pmn; 00181 Vector allPf1s_arranged = allPf1s; 00182 00183 for (i=0; i<(factorial(numLsf)); i++) { 00184 00185 arrangement = arrange(numLsf); 00186 00187 for (j=0; j<numLsf; j++) { 00188 allPf1s_arranged(j) = allPf1s(((int)arrangement(j)-1)); 00189 } 00190 for (j=0; j<numLsf; j++) { 00191 for (k=0; k<numLsf; k++) { 00192 Pmn_arranged(j,k) = Pmn(((int)arrangement(j)-1),((int)arrangement(k)-1)); 00193 } 00194 } 00195 00196 // Compute lower probability bound 00197 double lowerBound = allPf1s_arranged(0); 00198 double temp1, temp2; 00199 for (m=2; m<=numLsf; m++) { 00200 temp1 = 0.0; 00201 for (n=1; n<=n-1; n++) { 00202 temp1 += Pmn_arranged(m-1,n-1); 00203 } 00204 temp2 = allPf1s_arranged(m-1) - temp1; 00205 if (temp2 < 0.0) { 00206 temp2 = 0.0; 00207 } 00208 lowerBound += temp2; 00209 } 00210 00211 00212 // Compute upper probability bound 00213 double upperBound = allPf1s_arranged(0); 00214 double maximus; 00215 for (m=2; m<=numLsf; m++) { 00216 maximus = 0.0; 00217 for (n=1; n<=m-1; n++) { 00218 if (Pmn_arranged(m-1,n-1) > maximus) { 00219 maximus = Pmn(m-1,n-1); 00220 } 00221 } 00222 upperBound += allPf1s_arranged(m-1) - maximus; 00223 } 00224 00225 // Update bounds 00226 if (lowerBound < minLowerBound) { 00227 minLowerBound = lowerBound; 00228 } 00229 if (upperBound > maxUpperBound) { 00230 maxUpperBound = upperBound; 00231 } 00232 } 00233 00234 00235 // Print results (should do this over all defined systems) 00236 ofstream outputFile( fileName, ios::out ); 00237 00238 outputFile << "#######################################################################" << endln; 00239 outputFile << "# #" << endln; 00240 outputFile << "# SYSTEM ANALYSIS RESULTS #" << endln; 00241 outputFile << "# #" << endln; 00242 outputFile.setf(ios::scientific, ios::floatfield); 00243 outputFile << "# Lower probability bound: ........................... " 00244 <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<minLowerBound 00245 << " #" << endln; 00246 outputFile << "# Upper probability bound: ........................... " 00247 <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<maxUpperBound 00248 << " #" << endln; 00249 outputFile << "# #" << endln; 00250 outputFile << "# Warning: All permutations of the limit-state functions should #" << endln; 00251 outputFile << "# be checked. Currently, this is not done automatically. #" << endln; 00252 outputFile << "# The user is therefore advised to do this manually. #" << endln; 00253 outputFile << "# Contact the developer for further details/plans. #" << endln; 00254 outputFile << "# #" << endln; 00255 outputFile << "#######################################################################" << endln << endln << endln; 00256 00257 outputFile.close(); 00258 00259 00260 // INform the user on screen 00261 opserr << "System reliability analysis completed." <<endln; 00262 00263 return 0; 00264 } 00265 00266 00267 00268 double 00269 SystemAnalysis::functionToIntegrate(double rho, double beta1, double beta2) 00270 { 00271 double pi = 3.14159265358979; 00272 return 1.0/(2.0*pi*sqrt(1.0-rho*rho)) 00273 * exp(-(beta1*beta1+beta2*beta2-2*rho*beta1*beta2) 00274 /(2.0*(1.0-rho*rho))); 00275 } 00276 00277 int 00278 SystemAnalysis::factorial(int num) 00279 { 00280 int i = num-1; 00281 int result = num; 00282 00283 while (i>0) { 00284 result = result * i; 00285 i--; 00286 } 00287 00288 return result; 00289 } 00290 00291 Vector 00292 SystemAnalysis::arrange(int num) 00293 { 00294 Vector arrang(num); 00295 00296 // This is not yet implemented 00297 00298 for (int i=0; i<num; i++) { 00299 arrang(i) = i+1; 00300 } 00301 00302 return arrang; 00303 }
Generated on Mon Oct 23 15:05:25 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents