Resolution of large scale block-structured robust optimization problems

This archive provides a simple generic implementation of the (regularized) algorithm described in the paper

"Regularized decomposition of large scale block-structured robust optimization problems"

by Wim van Ackooij, Nicolas Lebbe and Jérôme Malick.

Robust block-structured optimization

We consider block-structured optimization problems of the following form

  (1)  minₓ ∑ᵢ fᵢ(xᵢ) + Ψ(x)
       s.t. xᵢ ∈ Xᵢ, for all i = 1,…,m
       
  x ∈ X ⊂ Rⁿ, n = ∑ᵢ nᵢ,
  Πᵢ Xᵢ = X with Xᵢ ⊂ Rⁿⁱ compact,
   Ψ : Rⁿ  → R U {+∞} convex,
  fᵢ : Rⁿⁱ → R convex

Where the coupling function Ψ is subject to uncertainties. We could for example assume that it depends on a parameter d lying in a given uncertainty set D ⊂ Rᵖ which leads to

  Ψ(x) = sup(d ∈ D) φ(d,x)

The algorithm use a bundle-based decomposition methods to tackle large scale instances and only need to solve multiple times each block with an additional linear and quadratic terms, that is to say

  minₓᵢ fᵢ(xᵢ) + bᵀx + 1/2 xᵀAx

Getting Started

Prerequisites

This program use the CPLEX and Eigen libraries and must be linked in the given Makefile by modifying respectively the CPLEXPATH and EIGENPATH variables.

Installing

The tests are then compiled using a simple

 $ make

command in the root directory of the archive.

Running the tests

Two examples are available in the tests directory.

Quadratic sum

A simple example is given in tests/test_simple.cpp where we solve problem (1) for m quadratic functions fᵢ = aᵢxᵢ² + bᵢxᵢ + cᵢ with Xᵢ = [-M,M] and a two-branch penalization function Ψ = sup(d ∈ [μ-k,μ+k]) φ(d,x) with φ(d,x) = max(c₁(∑ᵢxᵢ-d),c₂(∑ᵢxᵢ-d)).

To run this example :

 $ make
 $ ./bin/test_simple

Unit commitment

A more usefull example is given in tests/test_uc.cpp solving a basic unit commitment problem containing thermal units described below. Each unit Uᵢ (i = 0,…,8) use T=24 (the next 24 hours) variables xᵢₜ which correspond to the production of the unit at time t

Functions fᵢ :

The global cost for the schedule xᵢₜ of unit Uᵢ is given using a simple linear formula

fᵢ(xᵢₜ) = bᵢᵀxᵢ

Domain Xᵢ ⊂ R²⁴ :

For each unit the production is bounded from 0 to Mᵢ, so

0 ≤ xᵢₜ ≤ Mᵢ

and between two consecutive time steps t, t+1 the production might be modified to at most Kᵢ.

|xᵢₜ-xᵢₜ₊₁| ≤ Kᵢ,  t > 1

Lastly, each unit start with a given production xᵢ₀

|xᵢ₀-xᵢ₁| ≤ Kᵢ

Oracle Ψ :

We decide to modelize the uncertainty set D by a cubic set around the average load μ and a width of k.

For a given demand d ∈ D the cost φ penalize underproduction linearly by a factor c₁ ≤ 0 and surproduction by a factor c₂ ≥ 0 using the following formula

φ(d,y) = max(c₁(y-d),c₂(y-d)),

the Ψ function is then

Ψ(x) = ∑ₜ sup(dₜ ∈ Dₜ) φ(dₜ,(Ax)ₜ)

with (Ax)ₜ = ∑ᵢ xᵢₜ the production at time step t.

To run this example :

 $ make
 $ ./bin/test_uc

Documentation

Including the file rrbsopb.h you have access to the class rbsopb and rrbsopb whose unique constructors requires 3 parameters :

rrbsopb(VectorXi &ni, double(*fi)(int,VectorXd&,VectorXd&,VectorXd&), double(*psi)(VectorXd&,VectorXd&));
//  ni : vector of m integers containing the number of variables of each block
//  fi : pointer to a function of 4 variables : int i, VectorXd& x, VectorXd& b, VectorXd& A
//        i : the number of a block
//        x : the variable x which minimize fᵢ(xᵢ) + bᵀx + 1/2 xᵀdiag(A)x
//        b : vector corresponding to an additional linear term
//        A : vector corresponding to the diagonal of the matrix of an additional quadratic term
//       the function fi should return the value of minₓᵢ fᵢ(xᵢ) + bᵀx + 1/2 xᵀdiag(A)x
// psi : pointer to a function of 2 variables : VectorXd& x, VectorXd& g
//        x : point of evaluation of the function Ψ
//        g : a subgradient of Ψ(x)
//       the function psi should return the value of Ψ(x)

The files in the tests directory gives two simple examples using this constructor.

Once a rrbsopb (or rbsopb) is defined you can specify some parameters using the following methods :

// set type of primal recovery
// 0 : best-iterate, 1 : Dantzig-Wolfe-like
rbsopb* setPrimalRecovery(bool b);

// set if warm start is on or off
rbsopb* setWarmStart(bool b);

// set max number of cuts for the warm start
rbsopb* setMaxCutsWS(int nb);

// set maximum number of iterations for the bundle
rbsopb* setMaxInternIt(int nb);

// set maximum number of iterations for the algorithm
rbsopb* setMaxOuterIt(int nb);

// set precision for the bundle maximizing Θₖ(μ)
rbsopb* setInternPrec(double gap);

// set precision for the algorithm
rbsopb* setOuterPrec(double gap);

// set a time limit for the solver of each bundle iteration
rbsopb* setTimeLimit(double time);

// set a time limit for the dantzig-wolfe recovery
rbsopb* setPrimalTimeLimit(double time);

// set verbosity level
// 0 : no text displayed
// 1 : general information for each steps
// 2 : complete information for each inner and outer steps
rbsopb* setVerbosity(int n);

Then using the method double solve(VectorXd& x) you can solve the optimization problem, returning the final objective with the corresponding variable x.