15. Uncertainty modeling: Stochastic optimization models

Dealing with unpredictable events is very hard, but worthwhile: many decisions can be significantly improved if uncertainty and, in particular, variability is taken into correct account. The field is very large and indeed complex, and it involves modeling aspects as well as algorithmic and deep theoretical analysis. In these pages we will just scrap the surface of this important theme, mainly to show that it is indeed possible to take into consideration stochastic events and also that, in some cases, the resulting models are still linear ones, obtained thanks to modification of standard, deterministic models.

15.1. Probabilistic constraints

In many applications some data are not known with precision but, at best, an assumption can be made on a probabilistic model. In other words, some data can be viewed as the realization of a random variable. In these cases the traditional constraints cannot be used anymore. Consider the following case:

\begin{align*} \sum_j \param{a}_j \var{x}_j \leq \param{\textbf{b}} \end{align*}

where \(\param{\textbf{b}}\) is a random variable whose probability distribution

\[ F_{\param{\textbf{b}}} (y): = P (\param{\textbf{b}} \leq y) \]

is assumed to be known.

When, as in this case, \(\param{\textbf{b}}\) is a random variable, the constraint is no more well defined, since its right hand side is no more a constant, but a function (recall that a random variable is a real-valued function from the event space).

Often in these cases we can subsitute the constraint with a probabilistic one, called a chance constraint:

\begin{align*} P \left (\sum_j \param{a}_j \var{x}_j \leq \param{\textbf{b}} \right) \geq \param{\beta} \end{align*}

where \(\param{\beta} \in [0,1]\). This constraint corresponds to the request that the constraint is satisfied with probability at least equal to \(\param{\beta}\), a given threshold. In general chance constraints are nonlinear and often it is impossible to write them in an explicit form. In some cases, however, it is possible to convert such a constraint into an equivalent regular (deterministic) one. This happens in particular in the present case, that is, when the only random variable is the right hand side. In fact, assume that the probability distribution function of the random variable is known and that, for the sake of simplicity, \(\param{\textbf{b}}\) is a continuous random variable and, thus, that its distribution function \(F_\param{\textbf{b}}(\cdot)\) is invertible. We obtain:

\begin{align*} P \left (\sum_j \param{a}_j \var{x}_j \leq \param{\textbf{b}} \right) & = 1-P \left (\param{\textbf{b}} <\sum_j \param{a}_j \var{x}_j \right) \\ & = 1-P \left (\param{\textbf{b}} \leq \sum_j \param{a}_j \var{x}_j \right) \\ & = 1-F_{\param{\textbf{b}}} \left (\sum_j \param{a}_j \var{x}_j \right) \geq \param{\beta} & \text{iff} \\ \sum_j \param{a}_j \var{x}_j & \leq F_{\param{\textbf{b}}}^{-1} (1- \param{\beta}) \end{align*}

The right-hand side of this equation is deterministic, so the probabilistic constraint has been transformed into a traditional (and linear) constraint. As an example, if the distribution of \(\param{\textbf{b}}\) can be assumed to be uniformly distributed between two values​ \(\param{b}_{\ell}\) and \(\param{b}_{u}\), after easy calculations we can deduce the deterministic equivalent constraint

\begin{align*} \sum_j \param{a}_j \var{x}_j \leq \param{\beta} \param{b}_{\ell} + (1- \param{\beta}) \param{b}_u \end{align*}

If the term \(\param{\textbf{b}}\) was, instead, exponentially distributed with parameter \(\lambda\):

\begin{align*} P (\param{\textbf{b}} \leq y) & = 1 - \exp \left (- \param{\lambda} y \right) \end{align*}

for which it is easy to see that

\begin{align*} F_{\param{\textbf{b}}}^{-1} (\param{\beta}) & = - \frac{1} {\param{\lambda}} \log (1 - \param{\beta}) \end{align*}

then the deterministic equivalent form of the constraint would become

\begin{align*} \sum_j \param{a}_j \var{x}_j \leq - \frac{1} {\param{\lambda}} \log (1 - \param{\beta}). \end{align*}

In the frequent case in which the random variable is assumed to be normal (gaussian) the procedure is exactly the same, although we need to resort to a numerical approximation of the inverse gaussian distribution, whose analytical expression is not available.

Of course the situation becomes significantly more complex when more than one coefficient of the same constraint (or of different constraints) is random; even more complex is the case in which it is not possible to assume stochastic independence among these variables. In these cases, it is often impossible to obtain a closed form of a deterministic equivalent of the original constraint.

15.2. Stochastic linear optimization models

It is not possible, in this volume, to present even a very superficial analysis of probabilistic optimization models, even in the linear case. In this section we will just introduce elementary models from which it will be possible to understand how, in some cases, it is possible to deal with stochastic decision problems using traditional optimization models of the kind presented in this volume and how to solve them with standard optimization tools. Of course, this is not always possible: in fact, in general, the optimization models which arise when uncertainty is present can only be transformed into complex nonlinear optimization models.

Consider any of the problems described up to here in the volume: as a typical case, consider a product mix problem, as described in section Product Mix Models. Assume, differently from what has been done there, that the demand for a product is not known a priori but can at best be estimated assuming that a finite set of possible scenarios might happen in the future. As an example, assume that a pessimistic scenario, a neutral, and an optimistic one are the only possibilities we can foresee. In cases like this one, in which uncertainty is associated to a finite set of possible scenarios, it is possible to formulate a linear optimization model in order to optimize the expected value of the objective function (in this case, the expected gain). In order to be able to formulate a deterministic equivalent model, the usual technique consists in introducing additional variables whose index is associated with each of the possible scenarios.

We start with the simplest possible case: consider a situation in which there is a finite set of different scenarios, indexed by \(s \in \set{Scenarios}\) and, to each scenario, we associate a known probability \(\param{Prob}_s\). In the product mix example, let us also assume that the demand for a certain product is a function of the scenario: \(\param{Dem}_s\). The decision on how much to produce needs to be made before knowing which of the possible scenarios will indeed occur, otherwise the problem would not be a stochastic one anymore. Thus we have to face a decision problem under uncertainty. As in the deterministic case, we let \(\var{x}\) denote the quantity to be produced: this variable is not indexed by the scenario, as its value has to be decided a priori and in an independent way with respect to the future actual scenario. After the decision has been taken and the real scenario reveals itself, it might happen that the quantity produced is insufficient, or it might be too much. In the first case, a penalty \(\param{PenIns}\) will be paid for each unit of product which is not available to satisfy the actual demand; in the other case, a penalty \(\param{PenEcc}\) will be paid for the disposal of the quantity of product which is in excess with respect to the demand. As it is not known which of the scenarios will occur, we need to take into account the randomness of the demand by choosing a production level which minimizes a measure of the future uncertain cost. It is usual to employ, in these cases, the expected value of the cost of production, which is the sum of the cost associated to each possible scenario weighted by the probability that such scenario will actually occur. In this case, the product mix model will remain exactly the same, except that in the objective function an additional cost term appears:

\begin{align*} \param{PenIns} \sum_{s \in \set{Scenarios}} \max (\param{Dem}_s- \var{x}, 0) \param{Prob}_s & + \\ \param{PenEcc} \sum_{s \in \set{Scenarios}} \max (\var{x} - \param{Dem}_s, 0) \param{Prob}_s \end{align*}

If, as in the product mix general model, the objective is to maximize gain, the above cost should be subtracted from the original objective.

Although the unknown variable appears in a non-linear (even non-differentiable) function, it is easy to transform the expression into a fully linear one. It is indeed possible to “linearize” the model by introducing additional variables associated with the different scenarios. We can also easily extend the model to the usual case of more than one product. Let \(\var{Surplus}_{ps}\) and \(\var{Slack}_{ps}\) denote two new variable sets which, in the various scenarios, represent excess production or missing production with respect to actual demand of a specific product \(p\). We can add the following constraints to the model:

\begin{align*} \var{Surplus}_{ps} & \geq \var{x}_p - \param{Dem}_{ps} & \forall \, s \in \set{Scenarios}, p \in \set{Products} \\ \var{Surplus}_{ps} & \geq 0 \\ \var{Slack}_{ps} & \geq \param{Dem}_{ps} - \var{x}_p & \forall \, s \in \set{Scenarios}, p \in \set{Products} \\ \var{Slack}_{ps} & \geq 0 \end{align*}

and the cost term in the objective can now be reformulated in a linear way:

\begin{align*} \sum_{s \in \set{Scenarios}, p \in \set{Products}} \param{Prob}_s \left( \param{PenEcc} \var{Surplus}_{ps} + \param{PenIns} \var{Slack}_{ps}\right) \end{align*}

As an example consider a stochastic product mix problem with the following data:


param: RESOURCES: Avail:=
carpentry   1600
painting    1200
assembly    1500
verification 300
packing      500

param: PRODUCTS: Gain:=
prod-1 10
prod-2 18
prod-3 20
prod-4 25
prod-5 27
prod-6 28
prod-7 35

param Consumption:
                prod-1 prod-2 prod-3 prod-4 prod-5 prod-6 prod-7 :=
carpentry         7      5      9     10     10     12     15
painting          2      2      2      3      3      3      3
assembly          2      2      4      7      9     15     18
verification      1      1      1      2      1      2      2 
packing           1      1      1      1      2      1      0

param PenIns := 1;
param PenEcc := 0.1;

param Demand:
        OPT  AVG   NEG :=
prod-1   45   30    20
prod-2   60   30    25
prod-3   30   30    30
prod-4   60   45    40
prod-5   40   20    10
prod-6   40   20     0
prod-7   20   5      2

param Prob:=
OPT   0.3
AVG   0.4
NEG   0.3

In this example there are three possible future scenarios which influence the demand for most products (indeed, product 3 has a constant forecasted demand, whichever the scenario). We assume that, once a quantity is chosen for each product, the scenario becomes known and there might be surplus or shortages. The a priori optimization problem can be implemented as in the following model:


param Gain {PRODUCTS} >= 0;
param Avail{RESOURCES} >= 0;
param PenIns;
param PenEcc;
param Prob{SCENARIOS} >=0, <=1;
param MinQty{PRODUCTS} default 0;
param MaxQty{p in PRODUCTS} >= MinQty[p], default Infinity;

param Consumption {RESOURCES,PRODUCTS} >= 0;
#     resource consumption per product unit

var Qty{p in PRODUCTS} >= MinQty[p], <= MaxQty[p]; 
var Surplus{PRODUCTS, SCENARIOS} >= 0;
var Slack{PRODUCTS, SCENARIOS} >= 0;

maximize expected_gain:  
    sum {p in PRODUCTS, s in SCENARIOS} Gain[p] * Sold[p,s] * Prob[s]
    - sum{p in PRODUCTS, s in SCENARIOS} (Prob[s]*(Surplus[p,s]*PenEcc 
        + Slack[p,s]*PenIns));

subject to max_consumption {r in RESOURCES}:
   sum {p in PRODUCTS} Consumption[r,p] * Qty[p] <= Avail[r];

subject to def_sold_1{p in PRODUCTS, s in SCENARIOS}:
	Sold[p,s] <= Qty[p]; 

subject to def_sold_2{p in PRODUCTS, s in SCENARIOS}:
	Sold[p,s] <= Demand[p,s]; 

subject to def_surplus{p in PRODUCTS, s in SCENARIOS}:
    Surplus[p,s] >= Qty[p]-Demand[p,s];

subject to def_slack{p in PRODUCTS, s in SCENARIOS}:
    Slack[p,s] >= Demand[p,s]-Qty[p];

and, when run and solved through a linear optimization solver, we obtain the solution which maximizes the expected gain:

Qty [*] :=
prod-1  20
prod-2  45
prod-3  30
prod-4  45
prod-5  20
prod-6  20
prod-7   5

:          Sold Demand Slack Surplus    :=
prod-1 AVG   20    30    10      0
prod-1 NEG   20    20     0      0
prod-1 OPT   20    45    25      0
prod-2 AVG   30    30     0     15
prod-2 NEG   25    25     0     20
prod-2 OPT   45    60    15      0
prod-3 AVG   30    30     0      0
prod-3 NEG   30    30     0      0
prod-3 OPT   30    30     0      0
prod-4 AVG   45    45     0      0
prod-4 NEG   40    40     0      5
prod-4 OPT   45    60    15      0
prod-5 AVG   20    20     0      0
prod-5 NEG   10    10     0     10
prod-5 OPT   20    40    20      0
prod-6 AVG   20    20     0      0
prod-6 NEG    0     0     0     20
prod-6 OPT   20    40    20      0
prod-7 AVG    5     5     0      0
prod-7 NEG    2     2     0      3
prod-7 OPT    5    20    15      0

expected_gain = 3436.66

If, as it is quite common a practice, the demand in the different scenarios is averaged, we might try to solve a product mix problem using the expected demand in place of the real one. However, differently from the generic product mix problem introduced earlier in this volume, we should take into account the fact that some demand might be unsatisfied, as the limited availability of resources does not allow to produce all what is (expected) to be required. In this case we can modify the basic model and introduce a penalty for unsatisfied demand; on the other hand, in a deterministic product mix model there will be no surplus, as there would be no reason to produce in excess of the (expected) demand. Running this model we apparently obtain a gain which is superior to the expected one in the model seen before. But if we fix the order quantities and compute the expected gain under the three scenarios, we obtain the following:

Qty [*] :=
prod-1   9.07143
prod-2  37.5
prod-3  30
prod-4  48
prod-5  23
prod-6  20
prod-7   8.6

:              Sold   Demand    Slack  Surplus    :=
prod-1 AVG    9.07143    30    20.9286    0
prod-1 NEG    9.07143    20    10.9286    0
prod-1 OPT    9.07143    45    35.9286    0
prod-2 AVG   30          30     0         7.5
prod-2 NEG   25          25     0        12.5
prod-2 OPT   37.5        60    22.5       0
prod-3 AVG   30          30     0         0
prod-3 NEG   30          30     0         0
prod-3 OPT   30          30     0         0
prod-4 AVG   45          45     0         3
prod-4 NEG   40          40     0         8
prod-4 OPT   48          60    12         0
prod-5 AVG   20          20     0         3
prod-5 NEG   10          10     0        13
prod-5 OPT   23          40    17         0
prod-6 AVG   20          20     0         0
prod-6 NEG    0           0     0        20
prod-6 OPT   20          40    20         0
prod-7 AVG    5           5     0         3.6
prod-7 NEG    2           2     0         6.6
prod-7 OPT    8.6        20    11.4       0

expected_gain = 3361.03

and we can see that indeed this produces a smaller expected gain: in general it is not wise to substitute the expected demand with the stochastic one. Among many other motivations, substituting the stochastic nature of a random variable with its expectation completely neglects the variability, which is one of the most relevant characteristics of stochastic variables. In general, taking into proper account probability distributions when dealing with stochastic phenomena generates much more reliable and efficient decisions. Of course, this improved quality comes at the price of an increased complexity, both in data collection (as we need to retrieve information on all possible scenarios) and in model development.

The example just presented is very simple and concerns a deterministic decision which must be taken now, given possible future scenarios; in the example, no other decision is taken after the scenario reveals itself. There are indeed variables which depend on the scenario, but these are just direct consequences of the original choices and of the characteristics of the scenario.

In many practical applications decisions are taken sequentially: a first decision is taken, as we did above, then a scenario occurs and a second-stage decision, possibly containing a correction of the previous one, is taken; after that another scenario will occur, and so on. This situation frequently occurs, e.g., in prime material procurement: in textile production, as an example, we might need to decide how much of some special tissue to buy right now, taking into account long lead and transportation times. Then, after the buying agreement has been signed and possibly paid, some events, like fairs or exhibitions, take place and the situation of the future market might become less uncertain. At this point if, as an example, it is likely that the market is ready to buy more finished product than originally planned, new and usually more expensive prime material can be bought from nearby producers. This second stage decision depends on the original one and on the scenario which has been revealed. Eventually, the real demand is observed.

Quite a similar situation happens in agriculture. Consider the following atficial example:


Stochastic farm planning

This application is taken from the introductory chapters of [6] to which we address the reader interested in a very good introduction to the theory and methods of stochastic optimization. Assume a farmer has a quantity of acres of land available and that a decision has to be taken about how much of the available land to devote to Wheat, Corn, Beans. This decision has to be taken now. In order to run the farm and feed the cattle it would be necessary to harvest a minimum quantity of wheat and corn (200t and 240t, respectively). One of the difficulties is that the exact yield is uncertain, due to possibly different future weather conditions; however, in any case, the exact future yield is related to the quantity of land devoted to that cultivation. If needed, it is always possible to buy additional wheat and corn, or it is possible to sell it, if it is in excess. For what concerns beans, these are possibly grown only for sale (they are not needed for the cattle), so there is no minimum requirement.

We assume that future scenarios will be bad, average, good, and that the scenario will directly influence the yield of each product. In this simplified example, we assume prices are not sensible to the scenarios, although it might be interesting to add this possibility. The three scenarios are assumed to occur with equal probability. The relevant data is contained in the following data file:

set PRODUCTS := Wheat Corn Beans;
set SCENARIOS := Bad Avg Good;

param TotalLand := 500;

param: MinReq  SellPrice BuyingCost  PlantingCost :=
Wheat    200     170       238           150
Corn     240     150       210           230
Beans      0      36         .           260

param MaxBuy:=
    Beans 0
param MaxSell:=
    Beans 6000

param Yield:
    Bad   Avg   Good :=
Wheat    2     2.5   3
Corn     2.4   3     3.6
Beans   16    20    24

param Prob:=
Bad  0.333333
Avg  0.333334
Good 0.333333

The model is again quite similar to a product mix. We choose to add an extra variable to the model whose value is the effective product quantity, as a function of the yield. This is not necessary, indeed, but just increases the readability of the model. Moreover we introduce another variable, \(\var{ExtraSell}\), in the balance equation to model possible extra production which cannot be sold, as the market is not capable of absorbing it. In the example, this happens for Beans, for which a maximum demand is imposed equal to 6000; if, based on the acres devoted to Beans, in a scenario the yield turns out to be more than this quantity, the excess production will not contribute to any profit in the objective function: it will thus be a waste, a pure cost.

As it can be seen in the implementation, the main variable is not indexed by the scenario, as the decision about land subdivision is a first stage decision, to be taken now:


param Prob{SCENARIOS};

param TotalLand;

param MinReq{PRODUCTS};  # tons
param MaxBuy{PRODUCTS} default Infinity;
param MaxSell{PRODUCTS} default Infinity;

param SellPrice{PRODUCTS} default 0; # euro pr ton
param BuyingCost{PRODUCTS} default Infinity; # euro per ton
param Yield{PRODUCTS,SCENARIOS}; #tone per acre
param PlantingCost{PRODUCTS}; # euro per acre

var x{PRODUCTS} >=0, <= TotalLand;  # acres per product
var Sell{p in PRODUCTS, s in SCENARIOS} >= 0, <= MaxSell[p];
# excess tons of products to be sold
var Buy{p in PRODUCTS, s in SCENARIOS} >= 0, <= MaxBuy[p];
# excess tons of products to be bought
var Production{PRODUCTS,SCENARIOS} >=0 ;
var ExtraSell{PRODUCTS,SCENARIOS} >= 0;

maximize TotalGain: sum{p in PRODUCTS, s in SCENARIOS}
    Prob[s]*(SellPrice[p] * Sell[p,s] - x[p] * PlantingCost[p])
    - sum{p in PRODUCTS, s in SCENARIOS: BuyingCost[p] != Infinity}
    Prob[s]* BuyingCost[p] * Buy[p,s]; 

s.t. TotalAcres:
	sum{p in PRODUCTS} x[p] <= TotalLand;

s.t. Balance{p in PRODUCTS, s in SCENARIOS}:
        Production[p,s] - Sell[p,s] -ExtraSell[p,s] + Buy[p,s]
        = MinReq[p] ;

s.t. defProd{p in PRODUCTS, s in SCENARIOS}:
	Production[p,s] = x[p] * Yield[p,s];

After the main decision is taken and the scenario becomes known, second stage or recourse actions need to be taken: depending on the actual yield, some product quantity will be purchased or sold. The constraints of the model are simply the request not to exceed the total available land and a set of balance equations for the definition of possibly excess or defect production, with respect to the demand. The objective function is the weighted average of the total gain (revenues minus costs) in the scenarios, weighted with the probability of each scenario.

With these data, the optimal solution, taking into account the randomness of scenarios, is the following:

CBC 2.10.3: CBC 2.10.3 optimal, objective 108390.001
9 iterations

Scenario: Bad

	Prod 	Sell	Buy	Demand	Acres 	Yield 	
Wheat:	 340	 140	 0	 200	 170	 2	
Corn:	 192	 0	 48	 240	 80	 2	
Beans:	 4000	 4000	 0	 0	 250	 16	

Scenario: Avg

	Prod 	Sell	Buy	Demand	Acres 	Yield 	
Wheat:	 425	 225	 0	 200	 170	 2	
Corn:	 240	 0	 0	 240	 80	 3	
Beans:	 5000	 5000	 0	 0	 250	 20	

Scenario: Good

	Prod 	Sell	Buy	Demand	Acres 	Yield 	
Wheat:	 510	 310	 0	 200	 170	 3	
Corn:	 288	 48	 0	 240	 80	 4	
Beans:	 6000	 6000	 0	 0	 250	 24	

It is possible to see that the second stage decision on how much to buy or sell is quite easily understandable: we sell excess production, and buy any product which is needed. So there is no real difficulty, in this case, in planning the second stage decision. What is really strategic is to be able to take these recourse actions into due considerations so that the initial choice is made in a profitable way. If we changed the conditions, e.g., to that of an average scenario (this is simple to obtain by letting the probability of this scenario to be 1), we obtain a land destination of 120 acres to Wheat, 80 to Corn, the remaining 300 to Beans, with an overproduction of 100t of Wheat to be sold. The profit obtainable in this case is \(118\,600\), which apparently is more convenient that the optimal one, obtained in the stochastic case, \(108\,390\). However this is not true. If, in fact, the unpredictable nature of scenarios is taken into correct account, there is no more an advantage in using the average scenario. In fact, fixing the quantity of land as in the average scenario, taking into account the three different possibilities, it turns out that the expected return will be only \(103\,240.0154\) with an expected loss of roughly 5% with respect to the stochastic decision.


© Fabio Schoen 2023