OptimizationAnalysis.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.1 $ 00026 // $Date: 2003/10/27 23:45:41 $ 00027 // $Source: /usr/local/cvs/OpenSees/SRC/reliability/analysis/analysis/OptimizationAnalysis.cpp,v $ 00028 00029 00030 // 00031 // Written by Terje Haukaas (haukaas@ce.berkeley.edu) 00032 // 00033 00034 #include <OptimizationAnalysis.h> 00035 #include <ReliabilityDomain.h> 00036 #include <ReliabilityAnalysis.h> 00037 #include <Vector.h> 00038 #include <Matrix.h> 00039 #include <MatrixOperations.h> 00040 #include <NormalRV.h> 00041 #include <RandomVariable.h> 00042 #include <math.h> 00043 00044 #include <fstream> 00045 #include <iomanip> 00046 #include <iostream> 00047 using std::ifstream; 00048 using std::ios; 00049 using std::setw; 00050 using std::setprecision; 00051 using std::setiosflags; 00052 00053 00054 OptimizationAnalysis::OptimizationAnalysis(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 OptimizationAnalysis::~OptimizationAnalysis() 00065 { 00066 if (fileName != 0) 00067 delete [] fileName; 00068 } 00069 00070 00071 00072 int 00073 OptimizationAnalysis::analyze(void) 00074 { 00075 00076 // Alert the user that the FORM analysis has started 00077 opserr << "OptimizationAnalysis is running ... " << endln; 00078 00079 00080 // Assumptions 00081 // - Have start values for the design parameters x 00082 00083 00084 // Algorithm parameters 00085 double alphapar = 0.5; 00086 double betapar = 0.8; 00087 double gammapar = 2.0; 00088 double deltapar = 1.0; 00089 00090 // Bound on probability of failure 00091 double p0 = 0.00134989803163; 00092 double beta0 = 3.0; 00093 00094 // Related probability quantities... 00095 double pmax = p0; 00096 double pmin = 0.00001*p0; 00097 double betavarmax = -5.56; 00098 double betavarmin = beta0; 00099 00100 double covtarget = 0.11; 00101 00102 // Number of limit-state functions 00103 int K = 4; 00104 00105 // Declar 'tvector' 00106 Vector tvector(K,1); 00107 for (int k=0; k<K; k++) { 00108 for (int j=0; j<K; j++) { 00109 tvector(k,j) = 1.0; 00110 } 00111 } 00112 00113 // Even more parameters... 00114 double targetcost = 0.0; 00115 int maxouteriter = 25; 00116 int maxinneriter = 150; 00117 int N = 50000; 00118 00119 00120 // Number of random variables 00121 int size_u = 8; 00122 00123 // Discretization of u-space (one matrix per limit-state function!) 00124 Matrix Omega1(1,size_u); 00125 Matrix Omega2(1,size_u); 00126 Matrix Omega3(1,size_u); 00127 Matrix Omega4(1,size_u); 00128 00129 double betavar = beta0; 00130 00131 // Evaluate cost (objective function) 00132 cost = f0(x,betavar); 00133 00134 // Evaluate deterministic constraints 00135 convio = fjs(x,betavar,betavarmax,betavarmin); 00136 00137 termflag = 1; 00138 for (int i=1; i<1; i++) { 00139 00140 int ok = mooa(x_out, betavar, maxconstr, maxoverballs,u, 00141 x, 00142 alphapar, 00143 betapar, 00144 gammapar, 00145 deltapar, 00146 tvector, 00147 Omega, 00148 maxinneriter, 00149 size_u, 00150 assump, 00151 betavar, 00152 betavarmax, 00153 betavarmin); 00154 00155 00156 00157 // compute prob of failure 00158 [pf,beta,cov_pf,numsim,cputime,counter] = impsamp(x,u,i*N,0,covtarget - (0.004*i),assump); 00159 00160 00161 // Heuristic update of parameters t 00162 for (int kk=0; kk<K; kk++) { 00163 tvector(kk) = tvector(kk)*betavar/beta; 00164 } 00165 00166 } 00167 00168 00169 // Open output file 00170 ofstream outputFile( fileName, ios::out ); 00171 00172 outputFile << "#######################################################################" << endln; 00173 outputFile << "# FORM ANALYSIS RESULTS, LIMIT-STATE FUNCTION NUMBER " 00174 <<setiosflags(ios::left)<<setprecision(1)<<setw(4)<<lsf <<" #" << endln; 00175 outputFile << "# #" << endln; 00176 outputFile << "# No convergence! #" << endln; 00177 outputFile << "# #" << endln; 00178 outputFile << "#######################################################################" << endln << endln << endln; 00179 00180 // Clean up 00181 outputFile.close(); 00182 00183 // Print summary of results to screen (more here!!!) 00184 opserr << "OptimizationAnalysis completed." << endln; 00185 00186 return 0; 00187 } 00188 00189 00190 00191 00192 00193 00194 00195 00196 double 00197 OptimizationAnalysis::f0(Vector x, double betavar) 00198 { 00199 double Cs = 50.0; 00200 double Cc = 1.0; 00201 double alpha = 0.1; 00202 double Y3 = 18.30; 00203 00204 double X19 = 2*Y3*(1/x(7-1) + 1/x(8-1) + 1/x(9-1))/6; 00205 00206 double initialcost = 3*Cs*Y3*x(1-1)/4 + Cs*X19*x(6-1)*(x(3-1) + x(5-1) - alpha + 0.5*x(4-1)) + Cc*Y3*(x(2-1)*x(3-1) + x(4-1)*x(5-1)); 00207 00208 double costoffailure = 500*initialcost; 00209 00210 double cost = initialcost + costoffailure*normcdf(-betavar); 00211 00212 return cost; 00213 } 00214 00215 00216 Vector 00217 OptimizationAnalysis::fjs(Vector x, double betavar, double betavarmax, double betavarmin) 00218 { 00219 00220 double alpha = 0.1; 00221 00222 double X20b = (x(3-1) + x(5-1) - alpha)/2.0; 00223 double X20c = 24.0*0.0254; 00224 00225 Vector fjvector(25); 00226 00227 fjvector(1-1) = x(7-1) - X20b; 00228 fjvector(2-1) = x(7-1) - X20c; 00229 fjvector(3-1) = x(8-1) - X20b; 00230 fjvector(4-1) = x(8-1) - X20c; 00231 fjvector(5-1) = x(9-1) - X20b; 00232 fjvector(6-1) = x(9-1) - X20c; 00233 00234 fjvector(7-1) = 0.5*x(4-1) - x(3-1); 00235 fjvector(8-1) = x(2-1) - 4*x(4-1); 00236 fjvector(9-1) = x(4-1) - x(2-1); 00237 fjvector(10-1) = x(2) - 1.22; 00238 fjvector(11-1) = 0.15 - x(3-1); 00239 fjvector(12-1) = 0.15 - x(4-1); 00240 fjvector(13-1) = x(5-1)/x(4-1) - 4; 00241 fjvector(14-1) = 1 - x(6-1)/0.0001; 00242 fjvector(15-1) = -x(7-1); 00243 fjvector(16-1) = -x(8-1); 00244 fjvector(17-1) = -x(9-1); 00245 fjvector(18-1) = x(3-1) + x(5-1) - 1.2; 00246 00247 fjvector(19-1) = -g(x,zeros(8,1),9,assump); 00248 00249 // Assuming tension force in the steel is balanced by a concrete compression zone IN THE FLANGE 00250 // fjvector(20-1) = -g(x,zeros(8,1),5,assump); 00251 // fjvector(21-1) = -g(x,zeros(8,1),7,assump); 00252 // fjvector(22-1) = g(x,zeros(8,1),10,assump); 00253 00254 // Assuming tension force in steel is balanced by a concrete compression zone that propagates into the THE WEB 00255 fjvector(20-1) = -g(x,zeros(8,1),6,assump); 00256 fjvector(21-1) = -g(x,zeros(8,1),8,assump); 00257 fjvector(22-1) = g(x,zeros(8,1),11,assump); 00258 00259 00260 fjvector(23-1) = -x(5-1); 00261 fjvector(24-1) = betavar - betavarmax; 00262 fjvector(25-1) = -betavar + betavarmin; 00263 00264 00265 } 00266 00267 00268 double 00269 OptimizationAnalysis::g(Vector x, Vector u, int k) 00270 { 00271 00272 // Probability transformations 00273 Vector v(8); 00274 v(1-1) = u(1-1)*0.15*413.4*pow(10.0,6.0) + 413.4*pow(10.0,6.0); 00275 v(2-1) = u(2-1)*0.15*27.56*pow(10.0,6.0) + 27.56*pow(10.0,6.0); 00276 v(3-1) = u(3-1)*0.2*13.57*pow(10.0,3.0) + 13.57*pow(10.0,3.0); 00277 v(4-1) = u(4-1)*0.243*929*pow(10.0,3.0) + 929*pow(10.0,3.0); 00278 v(5-1) = u(5-1)*0.243*138.31*pow(10.0,3.0) + 138.31*pow(10.0,3.0); 00279 v(6-1) = u(6-1)*0.243*183.39*pow(10.0,3.0) + 183.39*pow(10.0,3.0); 00280 v(7-1) = u(7-1)*0.243*228.51*pow(10.0,3.0) + 228.51*pow(10.0,3.0); 00281 v(8-1) = u(8-1)*0.1*22.74*pow(10.0,3.0) + 22.74*pow(10.0,3.0); 00282 00283 double Y3 = 18.30; 00284 double alpha = 0.1; 00285 00286 double X15 = (2*x(4-1)/0.0254*(x(3-1)/0.0254 + x(5-1)/0.0254 - alpha/0.0254)*sqrt(v(2-1)/(6.89*pow(10.0,3.0))))*4.45; 00287 double X16 = x(6-1)*v(1-1)*(x(3-1) + x(5-1) - alpha)/x(7-1); 00288 double X17 = x(6-1)*v(1-1)*(x(3-1) + x(5-1) - alpha)/x(8-1); 00289 double X18 = x(6-1)*v(1-1)*(x(3-1) + x(5-1) - alpha)/x(9-1); 00290 00291 00292 double gval; 00293 00294 switch (k) { 00295 case 1: 00296 double X14 = 0.85*v(2-1)*(x(2-1) - x(4-1))*x(3-1)/v(1-1); 00297 double X10 = (x(1-1) - X14)*v(1-1)/(0.85*v(2-1)*x(4-1)); 00298 double X11 = (x(1-1) - X14)*v(1-1)*((x(3-1) + x(5-1) - alpha) - X10/2) + X14*v(1-1)*((x(3-1) + x(5-1) - alpha) - x(3-1)/2.0); 00299 gval = 1 - v(4-1)/X11 - (v(3-1)*Y3^2/8)/X11 - ((x(2-1)*x(3-1) + x(4-1)*x(5-1))*v(8-1)*Y3^2/8)/X11; 00300 break; 00301 case 2: 00302 double cap = X15 + X16; 00303 gval = 1 - v(5-1)/cap - (v(3-1)*Y3/6)/cap - ((x(2-1)*x(3-1) + x(4-1)*x(5-1))*v(8-1)*Y3/6)/cap; 00304 break; 00305 case 3: 00306 double cap = X15 + X17; 00307 gval = 1 - v(6-1)/cap - (v(3-1)*Y3/3)/cap - ((x(2-1)*x(3-1) + x(4-1)*x(5-1))*v(8-1)*Y3/3)/cap; 00308 break; 00309 case 4: 00310 double cap = X15 + X18; 00311 gval = 1 - v(7-1)/cap - (v(3-1)*Y3/2)/cap - ((x(2-1)*x(3-1) + x(4-1)*x(5-1))*v(8-1)*Y3/2)/cap; 00312 break; 00313 case 5: 00314 double X12 = x(1-1)/(x(2-1)*(x(3-1) + x(5-1) - alpha)); 00315 double X13 = (0.85*0.85*v(2-1)/v(1-1))*(87/(87 + v(1-1)/(6.89*pow(10.0,6.0)))); 00316 gval = -X12 + 0.75*X13; 00317 break; 00318 case 6: 00319 double X12 = x(1-1)/(x(4-1)*(x(3-1) + x(5-1) - alpha)); 00320 double X14 = 0.85*v(2-1)*(x(2-1) - x(4-1))*x(3-1)/v(1-1); 00321 double X13 = ((0.85*0.85*v(2-1)/v(1-1))*87/(87 + v(1-1)/(6.89*pow(10.0,6.0))) + X14/(x(4-1)*(x(3-1) + x(5-1) - alpha)))*x(4-1)/x(2-1); 00322 gval = -X12 + 0.75*X13; 00323 break; 00324 case 7: 00325 double X21 = 200/(1000*v(1-1)/(6.89*pow(10.0,6.0))); 00326 double X12 = x(1-1)/(x(2-1)*(x(3-1) + x(5-1) - alpha)); 00327 gval = -X21 + X12; 00328 break; 00329 case 8: 00330 double X21 = 200/(1000*v(1-1)/(6.89*pow(10.0,6.0))); 00331 double X12 = x(1-1)/(x(4-1)*(x(3-1) + x(5-1) - alpha)); 00332 gval = -X21 + X12; 00333 break; 00334 case 9: 00335 gval = -x(6-1)*(v(1-1)/pow(10.0,3.0)/6.89)/x(9-1)/2/x(4-1)*(1/sqrt(v(2-1)/(6.89*pow(10.0,3.0)))) + 4; 00336 break; 00337 case 10: 00338 gval = 1 - (0.85*v(2-1)*x(2-1)*x(3-1))/(v(1-1)*x(1-1)); 00339 break; 00340 case 11: 00341 gval = -1 + (0.85*v(2-1)*x(2-1)*x(3-1))/(v(1-1)*x(1-1)); 00342 break; 00343 default: 00344 break; 00345 } 00346 00347 return gval; 00348 00349 } 00350 00351 00352 00353 int 00354 OptimizationAnalysis::mooa(Vector x_out, double betavar, double maxconstr, double maxoverballs, Vector u, 00355 Vector x, 00356 double alpha, 00357 double beta, 00358 double gammapar, 00359 double delta, 00360 Vector tvector, 00361 Matrix Omega, 00362 int maxiter, 00363 int size_u, 00364 int assump, 00365 double betavar, 00366 double betavarmax, 00367 double betavarmin) 00368 { 00369 00370 00371 00372 00373 00374 Omega_basic = Omega; 00375 00376 epsilon = 0.1/(N*N); 00377 00378 K = length(tvector); // Number of random variables... 00379 00380 innersol = []; 00381 psi_0_bar = []; 00382 psi_omega = []; 00383 00384 Omeganames = fieldnames(Omega); 00385 00386 Vector uinitial(nrv); 00387 00388 00389 for (int N=1; N<maxiter+1; N++) { 00390 00391 // Compute approximate solution to inner problem(s) 00392 eps_inner = epsilon/K; 00393 00394 // For each limit-state function... 00395 // for (int kk=1; kk<=K; kk++) { 00396 innersol_temp = phinner(x, 00397 uinitial, 00398 alpha, 00399 beta, 00400 gammapar, 00401 delta, 00402 tvector(kk), 00403 eps_inner, 00404 kk, 00405 assump, 00406 betavar); 00407 // } 00408 00409 uinitial = innersol_temp(1:m,:); 00410 innersol = [innersol innersol_temp]; 00411 00412 // Construct new set omega_N 00413 fjvalues = fjs(x,assump,betavar,betavarmax,betavarmin); 00414 psi_0_bar = [psi_0_bar max(fjvalues)]; 00415 00416 Omega = Omega_basic; 00417 00418 eps_vector = eps_constdrop(N); 00419 00420 psi_omega = [psi_omega max([0; max(innersol(m + 1,(N - 1)*K + 1:N*K)); psi_0_bar(N)])]; 00421 00422 00423 for n = 1:N 00424 00425 psival = psi_omega(n); 00426 00427 for kk = 1:K 00428 00429 ctr = innersol(size_u + 1,(n - 1)*K + kk); 00430 00431 if abs(ctr - psival) < 10^(-15) & psival > eps_vector(n) 00432 00433 kfield = char(Omeganames(kk)); 00434 00435 oldvalues = getfield(Omega,kfield); 00436 00437 Omega = setfield(Omega,kfield,[oldvalues innersol(1:m,(n - 1)*K + kk)]); 00438 00439 end 00440 end 00441 end 00442 00443 [x,betavar,constr,theta,maxoverballs] = 00444 phouter(x, 00445 alpha, 00446 beta, 00447 gammapar, 00448 Omega, 00449 K, 00450 epsilon, 00451 assump, 00452 betavar, 00453 betavarmax, 00454 betavarmin, 00455 tvector); 00456 00457 00458 00459 newres = [x;betavar;constr;maxoverballs;theta]; 00460 00461 H = [H newres]; 00462 00463 N = N + 1; 00464 epsilon = eps_outer(N); 00465 00466 } 00467 00468 last_u = uinitial; 00469 00470 } 00471 00472
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