Types of robust constraints

A robust constraint can only contain uncertains either multiplied by a variable or at the right-hand side. The type of a robust constraint is defined by the type of uncertainty set (or sets) which define the uncertain parameters used in the constraint. We shall work with a small knapsack problem example and explore alternative formulations of the uncertainty set.

model Knapsack
 uses "mmrobust"                                 ! Load the robust library

 parameters
  NUM=5                                          ! Number of items
  MAXVAL=100                                     ! Maximum value
  MAXWEIGHT=80                                   ! Maximum weight
  WTMAX=102                                      ! Max weight allowed for haul
 end-parameters

 declarations
  Items=1..NUM                                   ! Index range for items
  VALUE: array(Items) of real                    ! Value of items
  WEIGHT: array(Items) of real                   ! Weight of items
  x: array(Items) of mpvar                       ! Decision variables
 end-declarations

 setrandseed(5);
 forall(i in Items) do
  VALUE(i):=50+random*MAXVAL
  WEIGHT(i):=1+random*MAXWEIGHT
 end-do
 forall(i in Items) x(i) is_binary               ! All x are 0/1

 MaxVal:= sum(i in Items) VALUE(i)*x(i)          ! Objective: maximize total value
 WtMax:= sum(i in Items) WEIGHT(i)*x(i) <= WTMAX ! Weight restriction

 maximize(MaxVal)

 writeln("Solution:\n Objective: ", getobjval)
 writeln("Item  Weight  Value")
 forall(i in Items)
   writeln(i, ": ", getsol(x(i)), strfmt(WEIGHT(i),8,2), strfmt(VALUE(i),8,2))

end-model

This basic model without uncertainty solves to 256.601. Its execution results in the following output display.

 Objective: 256.601
Item  Weight  Value
1: 0   74.37  118.04
2: 1   62.34  141.75
3: 1   27.74  114.85
4: 0   53.17  100.60
5: 0   74.02   65.82

Simple bounds on the uncertain coefficients

It is possible to impose simple box constraints on the uncertain values. This is often an immediate approach as it simply defines the range of values a coefficient can take, on a per-coefficient basis. Uncertains are introduced to the model as follows.

 parameters
  NUM=5                                        ! Number of items
  MAXVAL=100                                   ! Maximum value
  MAXWEIGHT=80                                 ! Maximum weight
  WTMAX=102                                    ! Max weight allowed for haul
  WTPERCENT=0.3                                ! Uncertainty as a percentage
 end-parameters

 declarations
  Items=1..NUM                                 ! Index range for items
  VALUE: array(Items) of real                  ! Value of items
  WEIGHT: array(Items) of real                 ! Weight of items
  x: array(Items) of mpvar
  WeightUncertainty: array(Items) of uncertain ! Uncertains representing
                                               ! deviation from weight
 end-declarations

 forall(i in Items) x(i) is_binary             ! All x are 0/1

 setrandseed(5);
 forall(i in Items) do
  VALUE(i):=50+random*MAXVAL
  WEIGHT(i):=1+random*MAXWEIGHT
 end-do

 MaxVal:= sum(i in Items) VALUE(i)*x(i)        ! Objective: maximize total value

 forall(i in Items) do
  WeightUncertainty(i) <=  WTPERCENT*WEIGHT(i) ! Uncertainty is a percentage of
                                               ! the expected weight
  WeightUncertainty(i) >=  0                   ! and only expected to go up this time
 end-do
 WtMax:= sum(i in Items) (WEIGHT(i) + WeightUncertainty(i))*x(i) <= WTMAX
                                               ! Weight restriction

 maximize(MaxVal)

Note the uncertain type for the array WeightUncertainty. This type is used in Mosel to define a parameter that is not known a priori by the user and is not under the user's control. The constraint WtMax that previously was of the type linctr now becomes a robustctr (the declaration of constraints is optional). The new model allows a 30% maximum deviation from the expected weight. Solving it we get 141.751.

 Objective: 141.751
Item  Weight  Value
1: 0   74.37  118.04
2: 1   62.34  141.75
3: 0   27.74  114.85
4: 0   53.17  100.60
5: 0   74.02   65.82

The solution has shifted to pick the most valuable item only. Increasing the weight of the selected item by the allowed 30% it becomes 81.042. The optimal solution is robust, i.e., for any binary solution with a larger objective function there exists a vector of uncertains that would violate the robust constraint. The objective has decreased significantly because of uncertainty; the difference in objective function values between the original problem and the robust version is often referred to as the price of robustness [BS2004]. If a subset of the uncertainty set is used, for example with a smaller value of WTPERCENT, the effect of uncertainty decreases and we can obtain a (still robust) solution with a larger objective function. On the contrary, increasing WTPERCENT is equivalent to increasing uncertainty, which reduces the objective function value of the robust optimal solution.

Note the explicit non-negativity restrictions on the uncertain coefficients. There is no default lower bound imposed on uncertains, in contrast to traditional decision variables; this is to reflect that uncertainty is not typically biased toward positive values.

Even if a robust modelling feature makes solving such a model simple, it is important to realize that when each individual uncertain is independently bounded by simple bounds only, the same model can be derived by modifying all coefficients in the way that restricts the constraint the most. In the knapsack example, as all variables must be nonnegative and all weights are also non-negative, this means adding the upper bound of the uncertains to the original coefficients.

 declarations
  Items=1..NUM                                 ! Index range for items
  VALUE: array(Items) of real                  ! Value of items
  WEIGHT: array(Items) of real                 ! Weight of items
  x: array(Items) of mpvar                     ! Decision variables
 end-declarations

 forall(i in Items) x(i) is_binary             ! All x are 0/1

 setrandseed(5);
 forall(i in Items) do
  VALUE(i):=50+random*MAXVAL
  WEIGHT(i):=1+random*MAXWEIGHT
 end-do

 MaxVal:= sum(i in Items) VALUE(i)*x(i)        ! Objective: maximize total value
 WtMax:= sum(i in Items) WEIGHT(i)*(1+WTPERCENT)*x(i) <= WTMAX
                                               ! Weight restriction

 maximize(MaxVal)

Indeed, the solution remains 141.751.

 Objective: 141.751
Item  Weight  Value
1: 0   74.37  118.04
2: 1   62.34  141.75
3: 0   27.74  114.85
4: 0   53.17  100.60
5: 0   74.02   65.82

Realization of uncertain parameters Recall that an uncertain is an unknown quantity that is not under our control, but is rather chosen by a third party ``opponent''. Solving a robust optimization problem amounts to solving a leader-follower game in which the user is the leader and the opponent obtains the variables from the user and chooses the uncertain parameters.

In some cases, after solving a robust optimization problem it is possible to obtain the value of the uncertain parameters that correspond to the robust solution. For this, it is necessary to have an optimal solution ^-x to the robust problem. Consider the robust constraint

∑_i=1 a_i u_i x_i ≤ b

and assume that the uncertains u_i are within the ellipsoidal uncertainty set defined by

∑_i=1 u_i² ≤ 4

The opponent implicitly chooses the value of the uncertain so that the left-hand side of the robust constraint is maximized (against the ≤ constraint). Hence a realization of the uncertains for this robust constraint can be found by maximizing ∑_i=1 a_i ^- x_i u_i subject to the quadratic uncertainty constraint ∑_i=1 u_i² ≤ 4. It is possible to obtain the values of the uncertain parameters through the same function used to obtain the variable, getsol. The modalities for use of getsol are explained in the Mosel manual. Note that getsol can be used on uncertains only after the robust counterpart has been solved.

Linear constraints on the uncertainty sets

Bound constraints are only one example of the possible uncertainty sets we can model. Let us suppose that some or all of the uncertains can be aggregated and that we can write a single uncertain constraint with all of them. Let us refine the example from the simple bounds case, and require instead that the sum of the uncertains is at most 30% of the total weight.

 parameters
  NUM=5                                        ! Number of items
  MAXVAL=100                                   ! Maximum value
  MAXWEIGHT=80                                 ! Maximum weight
  WTMAX=102                                    ! Max weight allowed for haul
  WTPERCENT=0.3                                ! Uncertainty as a percentage
 end-parameters

 declarations
  Items=1..NUM                                 ! Index range for items
  VALUE: array(Items) of real                  ! Value of items
  WEIGHT: array(Items) of real                 ! Weight of items
  x: array(Items) of mpvar                     ! Decision variables
  WeightUncertainty: array(Items) of uncertain ! Uncertains representing
                                               ! deviation from weight
 end-declarations

 forall(i in Items) x(i) is_binary             ! All x are 0/1

 setrandseed(5);
 forall(i in Items) do
  VALUE(i):=50+random*MAXVAL
  WEIGHT(i):=1+random*MAXWEIGHT
 end-do

 MaxVal:= sum(i in Items) VALUE(i)*x(i)        ! Objective: maximize total value

 forall(i in Items) do
  WeightUncertainty(i) >=  0
 end-do

 sum(i in Items) WeightUncertainty(i) <= WTPERCENT * sum(i in Items) WEIGHT(i)
 WtMax:= sum(i in Items) (WEIGHT(i) + WeightUncertainty(i))*x(i) <= WTMAX
                                               ! Weight restriction

 maximize(MaxVal)

While this sounds reasonable, because of the large right-hand side (30% of the total weight) the modified model is subject to a very conservative uncertainty set and will solve to a solution of all zeros.

 Objective: 0
Item  Weight  Value
1: 0   74.37  118.04
2: 0   62.34  141.75
3: 0   27.74  114.85
4: 0   53.17  100.60
5: 0   74.02   65.82

This behavior highlights how robust optimization works. When moving from the previous uncertainty set, with a small upper bound on the uncertains, to the new uncertainty set, we migh in fact relax the uncertainty set: if WTPERCENT * sum (i in Items) WEIGHT (i) is large enough, we have made the feasible space of the opponent significantly larger. This makes it possible that all uncertainty can be placed on a single item. Since each item can be too heavy even individually, the all zero solution becomes the only feasible one. This is a very important note to the behavior of robust optimization. When the uncertainty set is relaxed, the robust counterpart becomes more restrictive. Vice versa, to relax a robust constraint (and make it less conservative), its corresponding uncertainty set needs to be tightened.

A solution to the model in the example is then to restrict the feasible space of the uncertain values, that is, to tighten the uncertainty set by adding new uncertain constraints. In fact, these constraints would need to limit the amount of uncertainty per coefficient and will resemble the simple bound case in effect.

An uncertainty set defined by a set of linear constraints is often referred to as a polyhedral uncertainty set. It is important to note that polyhedral uncertainty sets have very important applications, but care must be taken when applying them as the example above has shown.

Ellipsoidal uncertainty sets

As we have seen with the previous example, the robust counterpart or opponent can take full advantage of the vertex solutions of a polyhedral uncertainty set making the model more conservative than possibly intended. This could be countered by either restricting such constraints or by refining the set, for example adding new linear constraints. It is easy to see that this is an approximation only to the underlying problem, and may be impractical in larger dimensions.

If the uncertainty set is known to be described by a quadratic constraint, for instance we know that the norm of the vector of the uncertains is bounded from above, then we can use ellipsoidal uncertainty sets. Another use for ellipsoidal uncertainty comes when the uncertainty set can be described by a confidence ellipsoid that is defined by a mean vector a, a covariance matrix Q, and a level of confidence α. In this case, the uncertainty set is of the form (u-a)^T Q (u-a)≤α. In the following model, we set α to WTPERCENT * sum(i in Items) WEIGHT(i).

 parameters
  NUM=5                                        ! Number of items
  MAXVAL=100                                   ! Maximum value
  MAXWEIGHT=80                                 ! Maximum weight
  WTMAX=102                                    ! Max weight allowed for haul
  WTPERCENT=0.3                                ! Uncertainty as a percentage
 end-parameters

 declarations
  Items=1..NUM                                 ! Index range for items
  VALUE: array(Items) of real                  ! Value of items
  WEIGHT: array(Items) of real                 ! Weight of items
  x: array(Items) of mpvar                     ! Decision variables
  WeightUncertainty: array(Items) of uncertain ! Uncertains representing
                                               ! deviation from weight
 end-declarations

 forall(i in Items) x(i) is_binary             ! All x are 0/1

 setrandseed(5);
 forall(i in Items) do
  VALUE(i):=50+random*MAXVAL
  WEIGHT(i):=1+random*MAXWEIGHT
 end-do

 MaxVal:= sum(i in Items) VALUE(i)*x(i)        ! Objective: maximize total value

 sum(i in Items) WeightUncertainty(i)^2 <= WTPERCENT * sum(i in Items) WEIGHT(i)
 WtMax:= sum(i in Items) (WEIGHT(i) + WeightUncertainty(i))*x(i) <= WTMAX
                                               ! Weight restriction

 maximize(MaxVal)

This model gives the desired solution as shown below. Note the value of the uncertains as obtained after maximization.

 Objective: 215.445
Item  Weight  Value  Uncertain
1: 0   74.37  118.04    0.00
2: 0   62.34  141.75    0.00
3: 1   27.74  114.85    6.61
4: 1   53.17  100.60    6.61
5: 0   74.02   65.82    0.00

Note that polyhedral uncertainty creates a robust counterpart that is of the same class as the original problem, as it amounts to adding a set of linear constraints. This is not the case with ellipsoidal uncertainty, which creates a robust counterpart that contains one second-order conic constraint for each robust constraint subject to ellipsoidal uncertainty, and hence may change the nature of the problem: if the original problem is a Linear Programming (LP) problem, the robust counterpart is a Second Order Conic Programming (SOCP) problem; if the original problem is a Mixed Integer Linear Programming (MILP) problem, the robust counterpart is a MISOCP. While solvers for MISOCP are now more powerful, a MISOCP problem is generally more difficult to solve than a MILP.

We conclude this section with a note about the use of getsol on uncertains: the last column of the output table shows the value of the uncertain WeightUncertainty for this particular solution. It is easy to notice that the nonzero uncertain is in accordance with the nonzero of the solution. This squarely aligns with viewing robust optimization as a leader-follower game: although the robust counterpart is solved just once, it emulates the dynamics of the opponent, who chose to increase the uncertain to make the robust constraint as tight as possible, and this can only be done by increasing WeightUncertainty(3) and WeightUncertainty(4).

Equality constraints and uncertain values with no uncertain constraints

At this point it is worth mentioning some pitfalls of robust optimization that one needs to be careful about. First, equality robust constraints are very restrictive in robust optimization.

 declarations
  x, y, z : mpvar
  e : uncertain
 end-declarations

 e <= 1
 e >= 0

 x + y + e*z = 10

 maximize(z)

In the solution to this problem, z will always be zero. The reason is simple: whenever z is equal to a non-zero value, if the value of variables x and y are fixed, any change in the value of the uncertain e will make the equality constraint infeasible. In general, the effect of uncertainty in equality robust constraints is to project the feasible space of the original problem to a smaller dimensional one.

Second, it is always beneficial to make sure that the uncertainty set is bounded from both below and above. The small example below shows the same effect as the case of the equality constraint.

 declarations
  x, y, z : mpvar
  e : uncertain
 end-declarations

 x + y + e*z <= 10

 maximize(z)

Here z will again solve to zero, as for any nonzero value of z a large enough e would make the constraint infeasible. Often, such missing bounds lead to an infeasible robust optimization problem.

Mixing uncertainty sets and types

Due to the theory of robust optimization, there are limitations on how uncertains can be used in uncertainty sets and robust constraints.

The most prominent restriction prohibits, by default, the use of the same uncertain values in different robust constraints. This restriction applies also to a model where two uncertains are used independently in two robust constraints if these two uncertains are contained in a constraint defining an uncertainty set. This restriction is enforced by virtue of how the corresponding robust counterparts are created, which allows the opponent to select its strategy, i.e. the value of an uncertain, on a per-constraint basis. Therefore, if the same uncertain appears in two robust constraints its value is allowed to take different values in each robust constraint. However, it is often convenient to make use of this behavior if the same uncertain set description is used for multiple robust constraints.

Consider the following example.

 declarations
  x, y: mpvar
  e, f : uncertain
 end-declarations

 x*e <= 1
 y*f <= 1

 e >= 0
 f >= 0
 e+f <= 1

 maximize(x+y)

The example will not solve using default settings, as the uncertain constraint e + f ≤ 1 connects the two robust constraints to the same uncertainty set, thereby they overlap. Setting ROBUST_UNCERTAIN_OVERLAP to true will allow for the problem to be solved.

 declarations
  x, y: mpvar
  e, f : uncertain
 end-declarations

 RobC1:= x*e <= 1
 RobC2:= y*f <= 1

 e >= 0
 f >= 0
 e+f <= 1

 setparam("XPRS_PROBNAME","h")

 setparam("ROBUST_UNCERTAIN_OVERLAP", true)

 maximize(x+y)

However, the solution returned is x=1=y, which can easily be seen as non-optimal if e + f ≤ 1 is satisfied: the solution is assuming that e=1 to satisfy the first robust constraint and that f=1 to satisfy the second one. Allowing the overlap will result in the opponent being able to optimize its decisions independently for the two robust constraints. This is the reason why this mode is not allowed by default; although often this is exactly what the modeller wants, setting ROBUST_UNCERTAIN_OVERLAP to true allows for reusing the same uncertainty set definition multiple times without the need to repeat the uncertain constraints. However, care must be taken as the uncertain values do not take a defined value upon solving the problem and it might be different for different robust constraints.

The retrieval of uncertains that overlap across different constraints adds a further twist: the two constraints x e≤ 1 and y f≤1 admit different combinations of the uncertains (e,f) as can be seen from the output.

x = 1  y = 1
constraint RobC1: x*e <= 1, uncertains (e,f) = (1,0)
constraint RobC2: y*f <= 1, uncertains (e,f) = (0,1)

The first constraint can only be made tighter by increasing uncertain e, while the second, whose only term is multiplied by f, can be made tight by setting the latter to its maximum value. For this reason, overlapping uncertains can only be retrieved by calling getsol passing one of the robust constraints where they appear. The value of that uncertain will then be equal to what the opponent would choose if that were the only robust constraint.

Cardinality restrictions for uncertains

Cardinality restrictions allow for limiting the number of uncertain values that are different from zero. In other words, a cardinality constraint limits the number of uncertains that are non-zero or, more generally, not at their nominal value. Due to the way in which the robust counterpart is created, uncertains that appear in a cardinality restriction should always have reasonable lower and upper bounds.

In the next example we show how to use uncertains in multiple robust constraints (i.e. in uncertains that overlap). The example is a small knapsack type problem, where only 1 of the 3 items is allowed to differ in weight from zero.

  declarations
    x,y,z: mpvar
    ex, ey, ez: uncertain
  end-declarations

  x is_binary ; y is_binary ; z is_binary

  20*x + 10*y + 5*z >= 20

  ex>=0 ; ex<=1
  ey>=0 ; ey<=1
  ez>=0 ; ez<=1

  cardinality({ex,ey,ez},1)   ! Allow only one of the uncertains to be nonzero

  setparam("ROBUST_UNCERTAIN_OVERLAP",true)   ! Allow overlaps

  10*(x-x*ex) + 10*(y-y*ey) + 10*(z-z*ez) >= 10
  10*(x-x*ex) + 10*(y-y*ey) + 10*(z-z*ez) <= 20

  maximize(x+y+z)

The solution returned is x=y=1 and z=0. It is important to emphasize that, as for the other constraints, the cardinality restriction is taken into consideration on a per robust constraint basis, and the model does not solve unless the ROBUST_UNCERTAIN_OVERLAP parameter is set to true.

Using historical data: scenarios

Scenarios are an efficient way of building robust optimization problems using historical realizations of the uncertains. Suppose we have a vector u of uncertains. We do not know an uncertainty set that would fit our model, but for the uncertain vector we have a database of records containing the value of each uncertain for each month of the past 20 years. Therefore we would like a model where one or more constraints are satisfied for each record of the uncertain vector in our database. Although we may not have a good approximation of the uncertainty set, robustness against historical data might be sufficient.

Declaring a scenario is equivalent to manually adding all the corresponding constraints for all realizations to the problem; this capability is made available to make sure that the scenario data are efficiently handled by the solver and to avoid the creation of an overly large problem.

Scenarios can be driven by the historical data, as shown in the next extension of the knapsack problem analyzed before where a number of historical measurements have been taken.

 parameters
  NUM=5                                  ! Number of items
  MAXVAL=100                             ! Maximum value
  MAXWEIGHT=80                           ! Maximum weight
  WTMAX=102                              ! Max weight allowed for haul
  UNCERTAINTY_LEVEL=0.3                  ! How much we are uncertain about the weight
  HISTORIC_PERIODS=100                   ! Number of scenarios
 end-parameters

 declarations
  Items=1..NUM                           ! Index range for items
  VALUE: array(Items) of real            ! Value of items
  WEIGHT: array(Items) of real           ! Weight of items
  UncertainWeight:array(Items) of uncertain
  x: array(Items) of mpvar               ! 1 if we take item i; 0 otherwise
  historical_weights: array(range, set of uncertain) of real
 end-declarations

 forall(i in Items) x(i) is_binary       ! All x are 0/1

 setrandseed(5);
 forall(i in Items) do
  VALUE(i):=50+random*MAXVAL
  WEIGHT(i):=1+random*MAXWEIGHT
 end-do

 MaxVal:= sum(i in Items) VALUE(i)*x(i)  ! Objective: maximize total value
 WtMax:= sum(i in Items) (WEIGHT(i)+UncertainWeight(i))*x(i) <= WTMAX

 ! Generate historical data, this would be data collected from actual realizations
 forall(period in 1..HISTORIC_PERIODS, i in Items)
   historical_weights(period, UncertainWeight(i)) :=
     WEIGHT(i)*UNCERTAINTY_LEVEL*random
 ! Generate a solution that would be feasible for ALL historic realizations
 scenario(historical_weights)

 maximize(MaxVal)

It is beneficial to examine the deterministic equivalent to a scenario based robust model. Consider the following simple scenario based model.

 declarations
  x, y : mpvar
  e, f : uncertain
  HISTDATA: array(range, set of uncertain) of real
 end-declarations

 ! Load historical data for e and f
 HISTDATA(1, e) := 1
 HISTDATA(1, f) := 2
 HISTDATA(2, e) := 3
 HISTDATA(2, f) := 1
 HISTDATA(3, e) := 3
 HISTDATA(3, f) := 2

 e*x + f * y <= 1

 ! Generate a solution that would be feasible for ALL historic realizations
 scenario(HISTDATA)

 maximize(x+y)

Its deterministic equivalent is:

 declarations
  x, y : mpvar
 end-declarations

 Scenario1 := 1 * x + 2 * y <= 1
 Scenario2 := 3 * x + 1 * y <= 1
 Scenario3 := 3 * x + 2 * y <= 1

 maximize(x+y)

Uncertainty in the objective function

When uncertain coefficients are introduced in the objective, the solution to the robust counterpart reflects the most conservative objective. More specifically, in a minimization problem where the objective function is f(x;u), where x is a vector of variables and u is a vector of uncertains, the optimal solution is one that minimizes the function g(x) = max_uf(x;u). This is equivalent to converting the objective to the robust constraint z≥f(x;u), where a new auxiliary variable z is used to represent the objective.