# Drawing samples from an LP program

Say I have an LP program in standard form:

\begin{equation*} \begin{array}{rl} \mathbf{x}^* = \underset{\mathbf{x}}{\text{arg}\;\text{min}} & \mathbf{c}^T\mathbf{x} \\ \mbox{s.t.} & \mathbf{A}\mathbf{x} = \mathbf{b} \\ & \mathbf{x} \ge 0. \end{array} \end{equation*}

Say that the $\mathbf{c}$ follows a normal distribution $N \left(\mathbf{\mu}, \mathbf{\sigma} \right)$, and that I draw multiple samples from it. What distribution would the corresponding $\mathbf{x}^*$ realizations follow? How could I proceed to analytically describe $P\left(\mathbf{x}^*\right | \mathbf{c})$ as a function of $\mu$ and $\sigma$?

Finally, and more generally, is there a topic in the literature that studies this sampling problem? (maybe in the context of optimization in general) Perhaps there are known relationships between this sampling problem and Sensitivity Analysis of LP solutions?

-
c is a vector. Do you want each element of the vector to have an independent $N(\mu,\sigma)$ distribution? Or only one element? – Mike Spivey Sep 26 '11 at 3:10
It seems like very interesting problem. There are few topics "Linear programming under uncertainty". Hopefully it can be useful from Management Science and from Operations Reseach – com Sep 26 '11 at 4:19

If you want to allow all elements in c to vary, I think you're going to have to do something like the following.

1. Find all basic feasible solutions. (This can be done by enumeration, although it might be time-consuming; if $A$ is $m \times n$, simply consider all possible ways to select $m$ of the $n$ variables to be in the basis. Then strike out the basic solutions that don't satisfy the nonnegativity constraints.)
2. For each basic feasible solution x', find its reduced cost vector $\bar{c}_N$. This will be a function of c. (The worst part of this step will be inverting the basis matrix.)
3. In order for x' to be the optimal minimum, its reduced cost vector must satisfy $\bar{c}_N \geq 0$. This is a system of inequalities in $c_1, c_2, \ldots, c_n$. Integrating the pdf of c over the solution set to the system of inequalities will yield the probability that x' is the optimal solution.

The integration in the last step will almost have to be done numerically, and I'm pretty sure that the system of inequalities will be nonlinear, so this will probably be more computationally intensive than you're wanting. Again, though, I'm not sure how to do this any other way.

In the case of letting a single element $c_i$ of the vector c have a $N(\mu,\sigma^2)$ distribution, though, there's a fairly standard procedure for finding the optimal solutions ${\bf x}^*$ that correspond to the different values of $c_i$. It's called parametric programming. See, for example, Section 7.2 in Hillier and Lieberman's Introduction to Operations Research (9th edition) or Section 10.4 of David Rader's Deterministic Operations Research. Once you've got the set of ${\bf x}^*$s that correspond to the different values of $c_i$, finding the probability distribution for ${\bf x}^*$ is straightforward.

Briefly, the idea is this: Solve the original problem for a fixed value of $c_i$, say, $0$. Then replace $0$ with $\theta$ in the objective function and recalculate the reduced cost vector $\bar{c}_N$. Since (for a minimization problem) we must have $\bar{c}_N \geq 0$, you can solve the resulting inequalities to find the range of values of $\theta$ that keep $\bar{c}_N \geq 0$; i.e., the values of $\theta$ that keep the current solution optimal. This will be something like $\theta_L \leq \theta \leq \theta_U$. (Some LP solvers will give you this - the range of values of a particular objective coefficient for which the current solution remains optimal - as part of the basic sensitivity analysis output.) Then, find the alternative optimal solutions at $\theta_L$ and $\theta_U$. For each solution found, find the $\theta$ interval for which that solution is optimal. Repeat the procedure until all values of $\theta$ have been included in some interval.

You'll end up with something like the following (which is from an example I do in my LP class): $$\begin{matrix} \theta \text{ Range } & \text{ Optimal Solution} \\ \theta \leq -1 & (0,10) \\ -1 \leq \theta \leq -\frac{2}{7} & (6,10) \\ - \frac{2}{7} \leq \theta \leq \frac{1}{4} & (24,4) \\ \theta \geq \frac{1}{4} & (30,0) \end{matrix}$$

Then, since $\theta$ is $N(\mu,\sigma^2)$, the probability distribution for ${\bf x}^*$ can be found easily by calculating the probabilities that $\theta$ falls in each of the ranges you just found. For example, $P({\bf x}^* = (0,10)) = P(\theta \leq -1) = P(Z \leq \frac{-1-\mu}{\sigma})$, where $Z$ is a standard normal random variable.

-
Thanks @Mike! This is excellent. I'll go through your answer more carefully tomorrow, but just to quickly clarify the original question, I'm interested in the general case with $\mathbf{c}$, the full vector, following a normal distribution $N(\mu,\sigma)$ (i.e. with $\mu$ and $\sigma$ vectors) – Amelio Vazquez-Reina Sep 26 '11 at 4:24
@intrpc: I've added some thoughts on the general case. – Mike Spivey Sep 26 '11 at 5:05