More efficient representation of stochastic programming problems

[vitaut]

Currently second stage subproblems are replicated for each scenario which is fine if the matrix is completely random, but not very efficient if only a small fraction of problem elements are random. For example, the size of the text .nl file for the airlift problem with the first data set is 12.4 KiB.

Getting rid of scenario indexing in second-stage variables and constraints and changing the satisfy_demand constraint to the following reduces the .nl file size to 1.7 KiB which is 2x smaller than the equivalent SMPS:

s.t. satisfy_demand{j in Routes}:
  load[j] - capacity_out[j] + capacity_in[j] +
    contracted[j] - unused[j] = random({s in Scen} Demand[j, s]);

The resulting problem should even be easier to convert to the solver form or SMPS because the core problem is already explicitly given.

[vitaut]

As an additional advantage, smpswriter will no longer have to process names.

[vitaut]

Since second-stage variables are no longer indexed over scenario set, probabilities should be passed differently, for example using explicit expectation:

maximize profit:
  expectation({s in Scen} P[s],
    ExcessSellingPrice * sell_excess +
    sum{c in Crops} (SellingPrice[c] * sell[c] -
                     PurchasePrice[c] * buy[c])) -
  sum{c in Crops} PlantingCost[c] * area[c]);

[vitaut]

Passing probabilities via expectation doesn't play well with independent random variables. Another option is to specify them together with realizations for each random variable/parameter:

function random;
demand{j in Routes}: random(
  {s in 1..NumScen} P[s], RandomDemand[j], {s in 1..NumScen} Demand[j, s]);

where P[s] is a parameter giving probability of scenario s, RandomDemand[j] is a random variable/parameter and {s in 1..NumScen} Demand[j, s] are its realizations.

Then expectation will only contain the expression:

minimize expected_cost:
  sum{i in AircraftTypes, j in Routes} AssignCost[i, j] * flights[i, j] +
  expectation(
    sum{(i, j, k) in Switches}
      (SwitchCost[i, j, k] -
        AssignCost[i, j] * (SwitchedHours[i, j, k] / Hours[i, j])) *
      switched_flights[i, j, k] +
    sum{j in Routes} ContractedCost[j] * contracted[j] +
    sum{j in Routes} UnusedCost[j] * unused[j]);

[vitaut]

The new approach is mostly implemented as of https://github.com/ampl/mp/commit/4d1de0535f73815d81b4163b8ea70ddab5a9568b and works for random RHS.

[vitaut]

Now need to handle random constraint matrix: #99