The stochastic response surface method (SRSM) is a technique for the reliability analysis of complex systems with low failure probabilities, for which Monte Carlo simulation (MCS) is too computationally intensive and for which other approximate methods may be inaccurate. Typically, the SRSM approximates a limit state function with a multidimensional quadratic polynomial by fitting the polynomial to a number of sampling points from the limit state function. This method can give biased approximations of the failure probability for cases in which the quadratic response surface can not conform to the true limit state function. In contrast to recently proposed algorithms that focus on the positions of sample points to improve the accuracy of the quadratic SRSM, this paper describes the use of higher order polynomials in order to approximate the true limit state more accurately. The use of higher order polynomials has received relatively little attention to date because of problems associated with ill-conditioned systems of equations and an approximated limit state which is very inaccurate outside the domain of the sample points. To address these problems, an algorithm using orthogonal polynomials is proposed to determine the necessary polynomial orders. Four numerical examples compare the proposed algorithm with the conventional quadratic polynomial SRSM and a detailed MCS.

  [Pf,g] = HOSRSM(limit_state,n,g_idx,tol,h_ord,h_reg,N,Xmcs)
 
  High Order Stochastic Response Surface Method (HOSRSM)
 
  This program approximates the probability of failure of a system using
  the high order stochastic response surface method (HOSRSM). The limit
  state function is returned by the matlab function called limit_state.m
  as follows:
      function G = limit_state(X)
  where X is the vector X of standardized normal random variables.
  Transformations from the standard normal form to the required
  distribution for each random variable is carried out in the
  limit_state() matlab function.

  G need not be a scalar.   If G is not a scalar, then the desired element
  of G is specified by the scalar integer input variable "g_idx".  
 
  The HOSRSM approximates the true limit state function with a
  polynomial reponse surface of arbitrary order,
     g(X) = a + \sum_{i=1}^n \sum_{j=1}^{k_i) b_ij X_i^j +
            \sum_{q=1}^m c_q \prod_{i=1}^n X_i^{p_iq}.
  The first stage of the algorithm determines the correct polynomial order,
  k_i in the response surace. Then the formulation of the mixed terms
  \sum_{q=1}^m c_q \prod_{i=1}^n X_i^{p_iq} are derived by the second stage
  based on previous results. In the third stage, the response surface
  is approximated. In the last stage, a Monte Carlo simulation is carried
  out to determine the probability of failure.
 
  INPUT
  ------------
  limit_state  the limit state function in the form
                 function g = limit_state(z)
               where z are sample values of the standardized normal
               random variables involved in the limit state function
  n            number of random variables in the limit state function
  g_idx        index of the desired element of the limit state function
  tol          accuracy tolerance for response surface mode  default = 0.01
  h_ord        parameter controlling the sampling domain in the polynomial
               order determination (1st stage)        default = 3
  h_reg        parameter controlling the sampling domain in the response
               surface approximation (3rd stage)      default = 3
  N            sample size for full scale Monte Carlo simulation (4th stage)
                                                      default = 1e6
  Zmcs         optional matrix of Monte-Carlo samples of standardized r.v's
               i.e., Zmcs = randn(N,n);
 
  OUTPUT
  ------------
  Pf           probability of failure
   g           values of the limit state function evaluated using
               Monte Carlo simulatiaon of the response surface

  Reference:
     Gavin, HP and Yau SC, ``High order limit state functions in the
     response surface method for structural reliability analysis,''
     submitted to Structural Safety, December 2005.
 
  Department of Civil and Environmental Engineering
  Duke University
  Siu Chung Yau, Henri P. Gavin, January 2006, Updated October 2009