Production planning under energy supply uncertainty

Topics covered in this section:

Problem description

A company produces liquid nitrogen and liquid oxygen (we will call them LIN and LOX in the remainder) and needs to plan production for the next N periods, identified as shifts of 8 hours each. Although the resource procurement is rather trivial (these two elements compose most of the troposphere), the process of obtaining these two liquid gases requires a vast amount of electricity, both for refrigerating the stored liquid gases and for powering the plant.

Given that the energy cost is the largest component of the production cost, power suppliers offer several types of energy supply contracts; one of these is known as Interruptible Load Contracts (or ILC for short) and allows the power supplier to interrupt, with a short notice (a few minutes), the provision of electricity to a large customer such as a plant for a limited period of time; most likely the power supplier will take advantage of this in times of high electricity demand, such as hot summer days. However, the ILC requires that there be no more than K interruptions throughout the planning horizon.

We are dealing therefore with a production planning problem under uncertainty in the power supply. We are given in input the production and inventory cost, the maximum production, the inventory capacity, the initial inventory, the demand in LIN and LOX for each period, and the maximum number K of interruptions. The problem consists of finding the amount of LIN and LOX to be produced every day so that the demands are satisfied irrespective of the interruptions that may occur within the clause of the contract.

This optimization problem requires the application of Robust Optimization for one reason among several: among the customers of this company there are hospitals, for which the fulfilment of the demand of LOX is imperative.

Mathematical formulation

We are given the set of gases G={lin,lox} and the set of time periods T={1,2,...,N}, the maximum production level P_lin and P_lox for each period and the inventory capacity S_lin and S_lox for each gas. The demand is given by the parameter dem_tg for t∈ T,g∈ G.

The uncertainty set is defined by every vector of interruptions (χ)_{i∈ T} such that χ_i ∈ [0,1] and not more than K of its elements are positive:

Let us define the variables prod_tg and inv_tg representing the production and inventory level, respectively, at period t for gas g. Without uncertainty in the input data, the model would have simple of constraints: bounds on the variables and fixing of the initial inventory level:

and the production planning constraint would consist of a mass conservation constraint for every gas g and time period t:

The modeling is more involved for this case. A first attempt would be to write the robust constraint by multiplying the production variable by (1-χ_t), so that the actual production at time period t is zero if χ_t=1:

However, this would yield no feasible solution: for each robust constraint related to production at time period t the uncertain χ_t would be set to one, and no production would occur. We need to aggregate the production constraints in such a way that the balance between production and demand yields an inventory that is nonnegative at every time period. The above constraint is replaced by the following one:

This robust constraint considers uncertainty in the interruptions but allows for a production plan that satisfies all demands even with K interruptions. In order to deal with contour conditions and take into account the initial inventory, we add a slight modification to the uncertainty set and require that no interruption happens at time period 1:

Finally, the robust value of the inventory must take into account the situation in which, after planning the production, no interruption actually occurs: this is simply enforced by bounding the inventory at every time period as the sum of all accumulated production discounted by the accumulated demand:

Implementation

Below is the implementation in Mosel: note that the option ROBUST_OVERLAP_UNCERTAIN is again used since the nonnegative-inventory constraints use the same set of uncertains multiple times.

model robust_prodplan
  uses "mmrobust"

  parameters
    plan_data = "prodplan_robust.dat"
  end-parameters

  declarations
    NDAYS:    integer                      ! Planning horizon
    PERDAY:   integer                      ! Number of periods per day
    NPERIODS: integer                      ! Total number of periods

    PERIODS,PERIODS0: range                ! Time periods
    GASES: set of string                   ! Set of products

    PROD_CAP: array(GASES) of real         ! Production capacity every day
    INV_CAP:  array(GASES) of real         ! Inventory capacity
    INV_0:    array(GASES) of real         ! Initial inventory

    PROD_COST: real                        ! Production cost
    INV_COST:  real                        ! Inventory cost

    DEMAND: array (PERIODS, GASES) of real ! Demand of each gas every day

    MAX_NINTERR: integer                   ! Maximum number of interruptions
                                           ! (as per contract)
  end-declarations

!**** Initialize data ****
  initializations from plan_data
    NDAYS PERDAY
    PROD_CAP  INV_CAP
    PROD_COST INV_COST INV_0
    DEMAND
    MAX_NINTERR
  end-initializations

  NPERIODS := NDAYS * PERDAY
  PERIODS0 := 0..NPERIODS

!**** Problem formulation ****
  declarations
    produce:   array(PERIODS, GASES) of mpvar   ! Production every day
    inventory: array(PERIODS0, GASES) of mpvar  ! Inventory level every day,
                                                ! including initial level

    interruption: array(PERIODS) of uncertain   ! Is power cut at this time?
    RobProd: array(PERIODS, GASES) of robustctr ! Robust constraint dealing with uncertainty
  end-declarations

!**** Constraints ****

  ! Start inventory
  forall(g in GASES)
    inventory (0,g) = INV_0 (g)

  ! Inventory balance
  forall(t in PERIODS, g in GASES) do

    RobProd(t,g) :=  inventory(0,g) + sum(tp in PERIODS | tp <= t)
      ((1 - interruption(tp)) * produce(tp,g) - DEMAND(tp,g)) >= 0

    inventory (0,g) + sum (tp in PERIODS | tp <= t)
      (produce(tp,g) - DEMAND(tp,g)) <= inventory(t,g)

    inventory(t,g) <= INV_CAP(g)
    produce(t,g)   <= PROD_CAP(g)

  end-do

  ! Interruptions of production
  forall (t in PERIODS) do
    interruption (t) <= 1
    interruption (t) >= 0
  end-do

  sum(t in PERIODS) interruption (t) <= MAX_NINTERR
  interruption(1) = 0

!**** Solving ****
  setparam("robust_uncertain_overlap", true)

  ! Set verbosity level
  setparam("xprs_verbose", true)

  ! Objective function: total cost of production and storage
  minimize(sum (t in PERIODS, g in GASES)
             (PROD_COST * produce (t,g) + INV_COST * inventory(t,g)))

!**** Solution reporting ****
  writeln("\nNumber of interruptions: ", MAX_NINTERR)
  writeln("\nOptimal solution has cost ", getobjval)

  COLWIDTH := 6

  forall(g in GASES) do

    writeln("\nProduction of ", g)

    write(strfmt ("Time",-COLWIDTH))
    forall(t in PERIODS0) write (strfmt(t,COLWIDTH))

    write("\n", strfmt("Dem",-2*COLWIDTH))
    forall(t in PERIODS) write(strfmt(DEMAND(t,g),COLWIDTH,1))

    write("\n", strfmt("Prod",-2*COLWIDTH))
    forall(t in PERIODS) write(strfmt (produce(t,g).sol,COLWIDTH,1))

    write("\n", strfmt("Inv*",-COLWIDTH))
    forall(t in PERIODS0)
      write (strfmt(inventory(0,g).sol +
            sum (tp in PERIODS | tp <= t)
              (produce(tp,g).sol - DEMAND(tp,g)),COLWIDTH,1))
    writeln("\n\nWorst-case interruptions for ", g, ": ")

    forall(t in PERIODS | getsol(interruption(t),RobProd(t,g)) > 0)
      write (t, " ")

    writeln
  end-do

end-model

Input Data

For this example, we consider a simple instance of the problem where the planning has a horizon of five days and three periods per day. The cost of production and inventory is 4 and 3, respectively, and the ILC contract allows for four interruptions. Table Production and inventory capacity below shows the production and inventory capacity and the initial inventory level for both gases.

Table Customer demand below lists the demand for each gas over all 15 periods.

Results

The optimal solution of this problem has cost 4727.17. Tables Production of LIN and Production of LOX give the optimal levels of production required in order to satisfy the demands regardless of when the interruptions occur, if ever. Note that the production levels are higher at the beginning and then, for liquid oxygen, they equal the demand. This depends on the fact that production at the initial time periods must assume for worst-case scenarios such as K=4 consecutive interruptions early in the planning horizon. The inventory levels equal the balance between demand and production as they must account for the event that no interruption occurs.

Note that this is the result of the optimization problem solved at the beginning of the time period. Using a rolling horizon approach, i.e., re-solving the model at every time period while taking into account the interruptions that already occurred, would allow us to obtain a less expensive production plan although the planned production levels might have to change.

Finally, since the interruption of a plant at a given period affects the production of both gases, the retrieval of uncertains needs to be done for each gas and getsol takes both interruption and the robust constraint RobProd as arguments. As would be expected, most worst-case interruptions are at the initial time periods, since those are critical moments when the LIN and LOX tanks need to be built up.

Worst-case interruptions for LIN:
2 3 4 5 6 14

Worst-case interruptions for LOX:
2 3 4 5 6 8 10 15 

References

This example is a simplified version of a real-world application that has been documented in [LBS2013].