2. MOTIVATION FOR IMPORTANCE SAMPLING

Suppose we wish to estimate ρ = E{φ(Z)}, where Z is a random variable describing some observation on a random system. Usually φ is the indicator function of some set implying that ρ is the probability of the set. Suppose that the observation random variable Z has probability density function p(·). The direct (Monte Carlo) simulation method would be to generate a sequence of independent and identically distributed (i.i.d.) random numbers Z₁,Z₂,…,Z_N from the density p(·) and form the estimate

By the law of large numbers,

. Therefore, as the number of observations approaches infinity, we converge to the true value. Suppose instead that we generate a sequence of i.i.d. random numbers Z₁,Z₂,…,Z_N using a possibly different density q(·). We then form the estimate

The ratio p(·)/q(·) will be called the weight function of the importance sampling estimator. It is simple to verify that the expected value of

under the density q(·) is precisely ρ.

Clearly, the estimate

is unbiased, and as N → ∞, we also expect it to be converging (by the law of large numbers) to its mean value ρ. The obvious question is: Are there better choices for q(·) than just p(·)? The answer is that by making a good choice for q(·), orders of magnitude decrease in the estimator variance can be achieved over a direct Monte Carlo simulation. It is this fact that has spurred most, if not all, of the recent interest in importance sampling techniques.

The important point is that the rate of convergence of the variance to zero as a function of N is precisely 1/N regardless of the importance sampling strategy chosen. The next section demonstrates that there are other estimation strategies that can give better rates of convergence in certain types of problems.

3. ANALYSIS OF RANDOM SEARCH STRATEGIES

Suppose we have some countable state space Ω. In the power system security problem [1], an element of Ω would be an ordered sequence of edges deleted from the network due to an initial failure event. Since a power system is a finite graph, Ω would also be finite (although usually very large). Even though the cardinality of Ω may be a very large finite number, on a short-time simulation scale due to the size of these state spaces they may behave as if they were infinite. Therefore, it is instructive to understand the more general countably infinite setting.

Suppose we have some probability measure P over the subsets of Ω and we are interested in estimating P(A) for some particular subset of Ω. For example, A could be the set of all ordered sequences that cause the graph to be disconnected. To estimate P(A), we choose a “search” probability measure F over the (Borel) subsets of Ω. In many situations, this measure F is not known explicitly. In general, it would be very hard to write down explicitly what this measure is for some arbitrary random path search on a graph.

We choose N independent samples {x_i ∈ Ω : i = 1,…,N} from the F distribution. We suppose that we can evaluate P(x_i) for each of these chosen values. We keep track of the total number of unique points that we find in set A. We suppose that in N trials, we find K such points.

Define

We can then write our estimator as

Of course, this estimator is biased. In fact (as pointed out in [1]), the estimator will always be a lower bound (a useful fact in itself) to the true value of P(A). To get a better notion of its performance, let us compute the mean squared error of the estimate:

where 1_{B} = 1 if B is true, and zero otherwise. Without loss of generality, we can assume P(x_i) > 0 for all x_i ∈ A (if not, just delete those points from A). From the above equation, we can immediately draw some important conclusions. If the state space Ω is finite, then we obtain that the mean squared error decreases exponentially (or geometrically) fast in N. This is an important point and we will return to it later.

When Ω is finite, let us define f₁ = min_xi∈A F(x_i) with minimum occurring at the value x_m1 [i.e., F(x_m1) = f₁]. Of course, this minimum could occur at several points and so x_m1 is not necessarily unique. For any x_i ∈ A and x_j ∈ A, we have

Suppose |Ω| < ∞; it is obvious that

where lim_N→∞ ε(N) = 0. Then we have the following theorem.

Theorem 3.1: Suppose |Ω| < ∞. Then

The limit log[1 − f₁] indicates how fast the mean squared error of the estimate converges to zero. The larger the value of f₁, the faster the mean squared error converges to zero. When Ω is finite, the distribution that has maximal f₁ is the uniform distribution. An immediate consequence of the theorem is the following.

Lemma 3.1: Suppose |Ω| < ∞. The optimal search distribution is the uniform distribution.

When the cardinality of the Ω space is countably infinite, determining the behavior of the random search estimator is quite a bit more complicated. Let us first consider a thought experiment. For simplicity (and without loss of generality), let us suppose that Ω = Z⁺, the nonnegative integers. Define P(j) = p_j and F(j) = f_j for all j ∈ Z⁺. Suppose p_i = C/i^1+γ (where C is a constant). Suppose we are interested in A = {i:i ≥ T}. Further, suppose that we decide to sample deterministically from Ω (instead of choosing random samples from some F). The best possible choice (which will not usually be possible) would be to choose the values T, T + 1, T + 2,…,T + N − 1 and form the estimate

. Thus, for the mean squared error,

This, of course, is the best possible rate for the mean squared error. Therefore, when γ < ½, the possible rate is worse than a relative frequency-type estimator, which always has a rate of 1/N.

One may ask: How much does one have to give up to use a random search strategy? Let us consider the γ = 1 case in a little more detail. We get the following useful upper bound:

Suppose that we choose f_i = C′/i^1+ε for some ε > 0. We can approximate the above upper bounding sum with the integral

Thus, by choosing δ (and ε) small enough, we can get arbitrarily close to the best possible rate 1/N². (Interestingly enough, when ε = 1 we can show that the integral in question has exactly rate 1/N. Therefore, a finer analysis than the above may be possible.)

The lesson to be gained from this analysis is that the sampling distribution should be in some sense as “heavy tailed” or as maximally spread out as we can get it. Of course, in the finite state space case, this is the uniform distribution.

In the power systems problem [1], the random search strategy works because the state space Ω is finite and fairly small. As a result, the guaranteed asymptotic convergence rate is observed early. The exponential rate promise (in the finite state space setting) is really conditioned upon visiting all the points in Ω with a high probability. Many graph problems have their complexity; hence the number of points in Ω increases exponentially fast as the size increases. Consequently, random search strategies are doomed to perform more poorly as the size of Ω increases. Obviously, the random search strategy would fail completely for an uncountably infinite state space where each point in the space has probability zero of occurrence. Another major drawback is the necessity to store all previously chosen samples to make sure that a given point in Ω is not used twice.

4. A UNIFORM SEARCH STRATEGY

As we saw in the previous section, the optimal strategy for a random search is to search uniformly over the state space. In a path search over a large graph, it is not at all clear how one goes about searching uniformly over the set of all paths leading to some failure set. Indeed, in most applications we do not even know the size of the failure set or the total number of paths. In this section, we propose a search technique based on the Metropolis Algorithm of Markov chain Monte Carlo, which allows us to easily specify a uniform search distribution.

Suppose we have an undirected graph G = (V,E). We are interested in the phenomenological effects of cascading link (link = edge) failures. We suppose that an initial event triggers a failure on one of the edges, say j ∈ E. Define the neighborhood edge set N(j) = {e ∈ E : e and j share a common vertex}. (Each element of N(j) represents a possible failure following j in practical power system cascading failure problems.) With the failure of j, we suppose that each of the elements of N(j) is at “risk” to fail also with given a priori probabilities p_e,e ∈ N(j). In the model we may allow these probabilities to change with the “state” of the graph after a chain of previous failures have occurred. We are interested in sequences of failures that disconnect the graph, cause large load losses, violate security criteria, have large numbers of lines lost—any or all of these. We denote the set of subgraphs of G that have these properties of large load lost, massive disconnection, and so forth that result from these unfortunate chains of events as B. Quite simply then, we wish to estimate the probability that the original graph G moves into one of the states B through a sequence of link failures. We suppose (although it is not completely necessary) that the disconnected graphs are always at least a subset of B.

Typical random searches suffer from the problem that they tend to find the short-path-length elements of B much more readily than the longer paths associated perhaps with greater calamities. We propose an algorithm based on techniques from Markov chain Monte Carlo (MCMC) to allow us to sample uniformly (or with any other distribution!) from the set B regardless of the individual element's path length.

We can think of the sequence of line outages as forming a tree. The original lost link j is the first node of the tree. The neighborhood set of that link then gives rise to the second level of the tree, which will have |N(j)| nodes. Each of these nodes will have a neighborhood set that will give rise to the third level of the tree and so on. Eventually, a sequence of link outages will terminate with a graph in B.

We first propose a Markov chain with state space Z, the nodes of the tree. We wish this Markov chain to have stationary distribution π. Usually we choose π to be uniform over Z. Let Z denote the nodes of the tree. We specify the transition structure of the Markov chain as

Q = (q_ij) is called the “candidate generating” transition matrix and is arbitrary as long as it is chosen to be irreducible. The idea is simple: When the present state is i, the next “tentative” state j is chosen with probability q_ij. When j ≠ i, this new state is “accepted” with probability α_ij. Hence, the probability of moving from i to j when j ≠ i is given as above.

We need to specify Q for our tree search. At node k ∈ Z, there are m_k nodes directly below it. We denote these as k_d,1,k_d,2,…,k_d,mk. There is, of course, only one node directly above node k, which we denote as k_u. We will only allow transitions between k and these neighboring nodes. Let 0 < p < 1 be chosen. We can then choose

When k is a leaf, there are no “down” nodes, and so we may choose

These choices, of course, give rise to an irreducible Q. Furthermore, we use the Metropolis algorithm choice for α_ij as

When π is uniform, this simplifies to

Thus, we have constructed a Markov chain with a uniform stationary distribution on the nodes of the tree. The leaf nodes are visited uniformly, and because every subgraph in B is associated with a unique leaf, we are also visiting these subgraphs uniformly.

5. A MODIFIED UNIFORM SEARCH STRATEGY

The strategy described in Section 4 constructs a Markov chain with a uniform stationary distribution on the nodes of the tree. Since we only need to produce a uniform sampling distribution on the leaf nodes, we can modify the algorithm to reduce the number of visits of nonleaf nodes while keeping a uniform stationary distribution on the leaf nodes. We denote L as the state space of all leaf nodes in Z. We start by modifying the transition matrix Q = (q_ij). In the previous section, q_ij is time invariant (i.e., it is not time dependent). Now we let q_ij also depend on the last state s. When the last state is s and the present state is i, the next “tentative” state j is chosen with probability q_ij(s). As in Section 4, we only allow transitions between neighboring nodes. At node k ∈ Z, let 0 < p_k < 1 be chosen. p_k can be either constant or k dependent, but it must be time invariant. If k is the root, let p_k ≡ 0. Using the same denotation as in Section 4, we then choose

If m_k = 1, we choose p_k = 1. If k is a leaf, then we choose q_{k ku} = 1.

Suppose we generate a random sequence {x_m|x_m ∈ Z} by applying transition matrix Q = (q_ij(s)) and assuming α_ij ≡ 1. Because Q is not time invariant, the sequence {x_m} is not a Markov chain. Let {y_n} be the subsequence of {x_m} that contains all of the leaf nodes in {x_m} (i.e., y_n ∈ L). Let

It is easy to prove that the matrix T = (t_ij) is time invariant and irreducible [7].

Now, we can construct a Markov chain with state space L, the leaf nodes of the tree. We specify the transition structure of the Markov chain as follows:

where t_ij is as defined earlier and β_ij is chosen as

The implementation is simple. Start from a leaf y_n, walk in the tree by the rule of Q until another leaf

is reached. Let

by the rule of β. Then start from y_n+1 to generate y_n+2 by the same rules. By the theory of the Metropolis algorithm, the generated sequence is the Markov chain with a uniform stationary distribution on the leaf nodes of the tree.

6. IMPORTANCE SAMPLING IS NOT DEAD

Denote the sequence of subgraphs from B generated by the above Metropolis algorithm as S₁,S₂,…,S_N with associated probabilities of occurrence p(S₁),p(S₂),…, p(S_N). Of course, this sequence of subgraphs is not independent, but the variance reduction technique of importance sampling does not require that in order to give an unbiased estimate.

We propose the importance sampling estimator

An exhaustive search algorithm would find the K << N unique graphs from the set {S₁,S₂,…,S_N} denoted as S⁽¹⁾,S⁽²⁾,…,S^(K) and form the estimate

The exhaustive search requires the saving of all the search paths and it becomes impractical when the size of the system is very large. The importance sampling estimator is unbiased but, unfortunately, requires the computation of |B|, which can be difficult. In applications applied to power systems, the computation of |B| is often unnecessary, especially when the study is to compare different dispatch strategies [7]. Reference 7 also gives a method to estimate |B| through simulation.

We simulated on an artificial tree to compare different important sampling strategies. The artificial tree is constructed to mimic the propagation of disturbances in large power systems. The details of the construction of the tree can be found in [7]. Figure 1 shows simulation results on a tree, which has about 2 × 10¹¹ leaf nodes. One of the reasons we choose to use an artificial tree is that we can construct the tree in a way that all of the statistical properties of the tree are computable; therefore we can evaluate the performance of our simulations conveniently. Simulation results from three sampling strategies are shown in Figure 1. IS1 (the solid line) represents the Metropolis algorithm described in Section 5. IS2 (the dashed line) represents a simple importance sampling strategy. In IS2, each leaf node sample is reached by starting a search from the root. At each node in the search path, a topology study is performed and a child node is chosen randomly. In this strategy, each child node has equal probability to be chosen. IS3 (the dot dashed line) differs from IS2 in that some child nodes have a higher priority to be chosen in the simulation so that some search paths have significant higher probabilities than other paths. Figure 1 shows that IS3 does not perform as well as the other two. This illustrates the importance of the proper choice of sample strategy in importance sampling techniques. The simulation results are consistent with our theory that the sampling distribution should be in some sense as “heavy tailed” or as maximally spread out as we can get it.

Simulation results.

In our simulation, IS1 and IS2 showed similar performance. This is because the structure of the artificial tree makes IS2 behave similarly to uniform sampling. We expect IS1 and IS2 could perform differently when the structure of the tree becomes more complex. We also want to point out that IS1 has less average computation cost per sample than IS2 does. For power system simulations, this means that, on average, IS1 requires fewer power flow calculations to get the same number of samples. However, the tuning of a Metropolis algorithm can be tricky [5].