Initializing help system before first use

# Production planning under energy supply uncertainty

Topics covered in this section:

- Problem description
- Mathematical formulation
- Implementation
- Input Data
- Results
- References

## 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 Plin and Plox for each period and the inventory capacity Slin and Slox for each gas. The demand is given by the parameter demtg 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:

∑i ∈ T χi ≤ K.

Let us define the variables prodtg and invtg 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:

0≤ prodtg ≤ Pg ∀ g, ∀ t∈ {1,2,...,N}
0≤ invtg ≤ Sg ∀ g, ∀ t∈ {1,2,...,N}
0≤ inv0g = Ig ∀ g;

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

invt-1,g + prodtg = invtg + demtg ∀ t ∈ T, g∈ G.

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:

invt-1,g + (1-χt)prodtg = invtg + demtg ∀ t ∈ T, g∈ G.

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:

inv0,g + ∑τ∈ T: τ ≤ t ((1-χτ)prodτ g - demτ g) ≥ 0 ∀ t ∈ T, g∈ G.

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:

χ1 = 0.

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:

invtg ≥ inv0,g + ∑τ∈ T: τ ≤ t (prodτ g - demτ g) ∀ t ∈ T, g∈ G.

## 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 18: Production and inventory capacity
Gas | Prod. capacity | Inv. capacity | Initial inv.
--- | --- | --- | ---
LIN | 29 | 60 | 20
LOX | 5.5 | 27 | 20

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

Table 19: Customer demand
Gas/Time | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
LIN | 6 | 14 | 10 | 8 | 11 | 15 | 10 | 9 | 10 | 9 | 10 | 12 | 11 | 15 | 9
LOX | 2 | 5 | 3 | 4 | 8 | 4 | 8 | 7 | 5 | 4 | 3 | 3 | 5 | 9 | 7

## 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.

Table 20: Production of LIN
Time | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
Dem | 6 | 14 | 10 | 8 | 11 | 15 | 10 | 9 | 10 | 9 | 10 | 12 | 11 | 15 | 9
Prod | 29 | 15 | 15 | 15 | 15 | 15 | 10 | 9 | 10 | 9 | 10 | 12 | 11 | 15 | 9
Inv* | 43 | 44 | 49 | 56 | 60 | 60 | 60 | 60 | 60 | 60 | 60 | 60 | 60 | 60 | 60

Table 21: Production of LOX
Time | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
Dem | 2 | 5 | 3 | 4 | 8 | 4 | 8 | 7 | 5 | 4 | 3 | 3 | 5 | 9 | 7
Prod | 5.5 | 5.2 | 5.2 | 5.2 | 5.2 | 5.2 | 5.2 | 5.2 | 5 | 5.2 | 5.2 | 5.2 | 5.2 | 5.2 | 5.2
Inv* | 23.5 | 23.7 | 25.8 | 27 | 24.2 | 25.3 | 22.5 | 20.7 | 20.7 | 21.8 | 24 | 26.2 | 26.3 | 22.5 | 20.7

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].

Parent Topic Robust Optimization Guide

© 2001-2026 Fair Isaac Corporation. All rights reserved. This documentation is the property of Fair Isaac Corporation (“FICO”). Receipt or possession of this documentation does not convey rights to disclose, reproduce, make derivative works, use, or allow others to use it except solely for internal evaluation purposes to determine whether to purchase a license to the software described in this documentation, or as otherwise set forth in a written software license agreement between you and FICO (or a FICO affiliate). Use of this documentation and the software described in it must conform strictly to the foregoing permitted uses, and no other use is permitted.