test_bds.cpp

/*--------------------------------------------------------------------------*/
/*---------------------------- File test_bds.cpp ---------------------------*/
/*--------------------------------------------------------------------------*/
/** @file
 * The forms of one two-stage stochastic investment problem, one against the
 * other on the same data.
 *
 * On the instances the generator writes there are three of them: the
 * extensive form, given to a :MILPSolver, which is the reference; the ad hoc
 * Benders one, an InvestmentBlock whose inner Block is the whole
 * TwoStageStochasticBlock, its InvestmentFunction being the value function
 * and FIXING the design in every scenario rather than mapping it into the
 * right-hand side, so that it produces one aggregated linearization per
 * iteration; and the generic one, i.e., the structure
 * BendersDecompositionSolver asks for, with the here-and-now Variable in a
 * single copy in the root and one LP sub-Block per scenario coupled to them
 * by its capacity Constraint alone.
 *
 * Being three formulations of one problem, the optima have to coincide,
 * which is what is checked; the iterations and the times are reported.
 *
 * The same structure, assembled around an instance that is on file rather
 * than built out of the data, is what benders_form() gives and what the
 * tester of the suite cross-checks the Solver on.
 *
 * \author Antonio Frangioni \n
 *         Dipartimento di Informatica \n
 *         Universita' di Pisa \n
 *
 * \author Donato Meoli \n
 *         Dipartimento di Informatica \n
 *         Universita' di Pisa \n
 *
 * \copyright © by Antonio Frangioni, Donato Meoli
 */
/*--------------------------------------------------------------------------*/
/*------------------------------ INCLUDES ----------------------------------*/
/*--------------------------------------------------------------------------*/

#include <chrono>
#include <cmath>
#include <iostream>
#include <netcdf>

#include "AbstractBlock.h"
#include "BendersDecompositionSolver.h"
#include "BlockSolverConfig.h"
#include "ColVariable.h"
#include "FRealObjective.h"
#include "FRowConstraint.h"
#include "InvestmentBlock.h"
#include "InvestmentFunction.h"
#include "LinearFunction.h"
#include "Objective.h"
#include "OneVarConstraint.h"
#include "TwoStageStochasticBlock.h"

using namespace SMSpp_di_unipi_it;
using namespace netCDF;

/*--------------------------------------------------------------------------*/
/*----------------------------- CONSTANTS ----------------------------------*/
/*--------------------------------------------------------------------------*/

// the relative gap at which two of these solves are the same number: the ad
// hoc form reads its objective back through a bundle whose required relative
// accuracy is 1e-4, hence nothing tighter can be asked of the comparison
static constexpr double tolerance = 1e-6;

/* The instance, the one in the repository unless another is named on the
 * command line: gen_investment.py writes them of whatever size. */

static std::string instance = "tssb_investment.nc4";

/*--------------------------------------------------------------------------*/
/*------------------------------- TYPES ------------------------------------*/
/*--------------------------------------------------------------------------*/

/// the instance, read once and used to build every form of it
/** Whatever the form, the data is the same, and it is read from the very file
 * the ad hoc form is deserialized from, so that the two cannot drift apart. */

struct Data {
 double cost;                       ///< cost of one unit of design
 double lb;                         ///< lower bound on the design
 double ub;                         ///< upper bound on the design
 std::vector< double > weight;      ///< one probability per scenario
 std::vector< std::vector< double > > demand;   ///< [ scenario ][ time ]
 std::vector< double > unit_max;    ///< max power, per unit of design
 std::vector< double > unit_cost;   ///< cost of the design-carrying unit
 std::vector< double > slack_max;   ///< max power of the slack
 std::vector< double > slack_cost;  ///< cost of the slack
 };

/*--------------------------------------------------------------------------*/
/*----------------------------- READING ------------------------------------*/
/*--------------------------------------------------------------------------*/

static std::vector< double > get_vec( const NcGroup & g ,
				      const std::string & name )
{
 auto v = g.getVar( name );
 std::size_t n = 1;
 for( int i = 0 ; i < v.getDimCount() ; ++i )
  n *= v.getDim( i ).getSize();
 std::vector< double > out( n );
 v.getVar( out.data() );
 return( out );
 }

/*--------------------------------------------------------------------------*/

static Data read( const std::string & file )
{
 Data d;
 NcFile f( file , NcFile::read );

 auto inv = f.getGroup( "Block_0" );
 d.cost = get_vec( inv , "Cost" ).front();
 d.lb = get_vec( inv , "LowerBound" ).front();
 d.ub = get_vec( inv , "UpperBound" ).front();

 auto tssb = inv.getGroup( "InnerBlock" );
 auto set = tssb.getGroup( "DiscreteScenarioSet" );
 d.weight = get_vec( set , "PoolWeights" );

 auto sc = set.getVar( "Scenarios" );
 const auto K = sc.getDim( 0 ).getSize() , T = sc.getDim( 1 ).getSize();
 auto flat = get_vec( set , "Scenarios" );
 d.demand.resize( K );
 for( std::size_t k = 0 ; k < K ; ++k )
  d.demand[ k ].assign( flat.begin() + k * T , flat.begin() + ( k + 1 ) * T );

 auto uc = tssb.getGroup( "StochasticBlock" ).getGroup( "Block" );
 auto u0 = uc.getGroup( "UnitBlock_0" ) , u1 = uc.getGroup( "UnitBlock_1" );
 d.unit_max = get_vec( u0 , "MaxPower" );
 d.unit_cost = get_vec( u0 , "ActivePowerCost" );
 d.slack_max = get_vec( u1 , "MaxPower" );
 d.slack_cost = get_vec( u1 , "ActivePowerCost" );

 return( d );
 }

/*--------------------------------------------------------------------------*/
/*--------------------------- BLOCK BUILDERS -------------------------------*/
/*--------------------------------------------------------------------------*/

/* BundleSolver does not take a general bound l <= x <= u on the Variable of
 * the Block it is attached to, so the design is shifted: the Variable is
 * x' = x - l in [ 0 , u - l ], and l is added back wherever x appears, which
 * leaves a constant cost * l in the Objective. The ad hoc form does the very
 * same thing, through the BlockConfig of the InvestmentBlock, so the two keep
 * describing one problem. */

static std::vector< ColVariable > * make_x( const Data & d )
{
 auto x = new std::vector< ColVariable >( 1 );
 ( * x )[ 0 ].is_positive( true );
 ( * x )[ 0 ].set_value( 0 );
 return( x );
 }

/*--------------------------------------------------------------------------*/

// the box on the shifted design, the only Constraint of the master

static void add_design_bound( AbstractBlock * block , const Data & d ,
			      std::vector< ColVariable > * x )
{
 auto bnd = new std::vector< BoxConstraint >( 1 );
 ( * bnd )[ 0 ].set_variable( & ( * x )[ 0 ] );
 ( * bnd )[ 0 ].set_lhs( 0 );
 ( * bnd )[ 0 ].set_rhs( d.ub - d.lb );
 block->add_static_constraint( * bnd , "design bound" );
 }

/*--------------------------------------------------------------------------*/

/* add to block the dispatch of one scenario: the power of the design-carrying
 * unit and of the slack, the demand Constraint and the capacity Constraint,
 * the latter being the only place where the design Variable appears, and the
 * costs of the scenario, weighted by its probability, to obj */

static void add_scenario( AbstractBlock * block , const Data & d ,
			  std::size_t k , std::vector< ColVariable > * x ,
			  LinearFunction * obj )
{
 const auto T = d.demand[ k ].size();

 auto g = new std::vector< ColVariable >( T );
 for( auto & v : * g ) v.is_positive( true );
 block->add_static_variable( * g , "g" );

 auto s = new std::vector< ColVariable >( T );
 for( auto & v : * s ) v.is_positive( true );
 block->add_static_variable( * s , "s" );

 auto dem = new std::vector< FRowConstraint >( T );
 for( std::size_t t = 0 ; t < T ; ++t ) {
  auto f = new LinearFunction();
  f->add_variable( & ( * g )[ t ] , 1 );
  f->add_variable( & ( * s )[ t ] , 1 );
  ( * dem )[ t ].set_function( f );
  ( * dem )[ t ].set_both( d.demand[ k ][ t ] );
  }
 block->add_static_constraint( * dem , "demand" );

 auto cap = new std::vector< FRowConstraint >( T );
 for( std::size_t t = 0 ; t < T ; ++t ) {
  auto f = new LinearFunction();
  f->add_variable( & ( * g )[ t ] , 1 );
  f->add_variable( & ( * x )[ 0 ] , - d.unit_max[ t ] );
  ( * cap )[ t ].set_function( f );
  ( * cap )[ t ].set_lhs( -Inf< double >() );
  ( * cap )[ t ].set_rhs( d.unit_max[ t ] * d.lb );
  }
 block->add_static_constraint( * cap , "capacity" );

 auto slk = new std::vector< FRowConstraint >( T );
 for( std::size_t t = 0 ; t < T ; ++t ) {
  auto f = new LinearFunction();
  f->add_variable( & ( * s )[ t ] , 1 );
  ( * slk )[ t ].set_function( f );
  ( * slk )[ t ].set_lhs( -Inf< double >() );
  ( * slk )[ t ].set_rhs( d.slack_max[ t ] );
  }
 block->add_static_constraint( * slk , "slack capacity" );

 for( std::size_t t = 0 ; t < T ; ++t ) {
  obj->add_variable( & ( * g )[ t ] , d.weight[ k ] * d.unit_cost[ t ] );
  obj->add_variable( & ( * s )[ t ] , d.weight[ k ] * d.slack_cost[ t ] );
  }
 }

/*--------------------------------------------------------------------------*/

// the monolithic model, every scenario in one Block: the reference optimum

static AbstractBlock * build_monolithic( const Data & d )
{
 auto block = new AbstractBlock();
 auto x = make_x( d );
 block->add_static_variable( * x , "x" );
 add_design_bound( block , d , x );

 auto f = new LinearFunction();
 f->add_variable( & ( * x )[ 0 ] , d.cost );
 f->set_constant_term( d.cost * d.lb );
 for( std::size_t k = 0 ; k < d.demand.size() ; ++k )
  add_scenario( block , d , k , x , f );

 auto obj = new FRealObjective( block , f );
 obj->set_sense( Objective::eMin );
 block->set_objective( obj );
 return( block );
 }

/*--------------------------------------------------------------------------*/

// the 2-level model: the design in the root, one LP sub-Block per scenario

static AbstractBlock * build_structured( const Data & d )
{
 auto root = new AbstractBlock();
 auto x = make_x( d );
 root->add_static_variable( * x , "x" );
 add_design_bound( root , d , x );

 auto df = new LinearFunction();
 df->add_variable( & ( * x )[ 0 ] , d.cost );
 df->set_constant_term( d.cost * d.lb );
 auto robj = new FRealObjective( root , df );
 robj->set_sense( Objective::eMin );
 root->set_objective( robj );

 for( std::size_t k = 0 ; k < d.demand.size() ; ++k ) {
  auto sub = new AbstractBlock( root );
  auto sf = new LinearFunction();
  add_scenario( sub , d , k , x , sf );
  auto sobj = new FRealObjective( sub , sf );
  sobj->set_sense( Objective::eMin );
  sub->set_objective( sobj );
  root->add_nested_Block( sub );
  }

 return( root );
 }

/*--------------------------------------------------------------------------*/
/*------------------------------ SOLVING -----------------------------------*/
/*--------------------------------------------------------------------------*/

static double solve( Block * block , const std::string & cfg , double & secs ,
		     long & iters , int & status )
{
 auto c = Configuration::deserialize( cfg );
 auto bsc = dynamic_cast< BlockSolverConfig * >( c );
 if( ! bsc ) {
  std::cerr << cfg << " is not a BlockSolverConfig" << std::endl;
  std::exit( 1 );
  }
 /* Whatever goes wrong, from the configuration being refused onwards, the
  * Block is given back what the Solver had taken from it before the error
  * travels on, so that the caller can go on using the Block. */

 auto give_back = [ & ]() { bsc->clear(); bsc->apply( block ); delete bsc; };

 double value;
 try {
  bsc->apply( block );

  auto solver = block->get_registered_solvers().front();

  auto t0 = std::chrono::steady_clock::now();
  status = solver->compute();
  secs = std::chrono::duration< double >(
                              std::chrono::steady_clock::now() - t0 ).count();
  iters = solver->get_elapsed_iterations();
  value = solver->get_var_value();
  }
 catch( ... ) { give_back(); throw; }

 give_back();
 return( value );
 }

/*--------------------------------------------------------------------------*/

/* Both Benders forms drive their master with a BundleSolver, and the master
 * keeps the subproblems hanging from it as sub-Block; a bundle that wants
 * each of its sub-Block to be a bare function Block refuses it in
 * set_Block(), and so does one that does not carry the required parameters.
 * In either case the form is skipped rather than failed. */

static double solve_by_bundle( Block * block , const std::string & cfg ,
			       double & secs , long & iters , int & status ,
			       bool & ran , double dflt )
{
 try {
  const auto value = solve( block , cfg , secs , iters , status );
  ran = true;
  return( value );
  }
 catch( const std::exception & e ) {
  std::cout << cfg << ": skipped, " << e.what() << std::endl;
  ran = false;
  return( dflt );
  }
 }

/*--------------------------------------------------------------------------*/

/* the ad hoc form: the InvestmentBlock is deserialized from the file and
 * given its OBlockConfig [see IBOCfg.txt], which reformulates its bound
 * Constraint into the shifted box that a BundleSolver takes, and gives the
 * InvestmentFunction the BlockSolverConfig of the inner Block */

static double solve_ad_hoc( double & secs , long & iters , int & status ,
			    bool & ran , double dflt )
{
 auto root = Block::deserialize( instance );
 auto inv = dynamic_cast< InvestmentBlock * >( root );
 if( ! inv ) {
  std::cerr << instance << " does not hold an InvestmentBlock" << std::endl;
  std::exit( 1 );
  }

 auto config = dynamic_cast< BlockConfig * >(
                       Configuration::deserialize( "IBOCfg.txt" ) );
 if( ! config ) {
  std::cerr << "IBOCfg.txt is not a BlockConfig" << std::endl;
  std::exit( 1 );
  }

 // the OBlockConfig gives its ComputeConfig to the Objective, which has to
 // be there already
 inv->generate_abstract_variables();
 inv->generate_objective();
 config->apply( inv );
 delete config;

 const auto value = solve_by_bundle( inv , "BSPar-Inv.txt" , secs , iters ,
				     status , ran , dflt );
 delete root;
 return( value );
 }

/*--------------------------------------------------------------------------*/
/*--------------------------------- MAIN -----------------------------------*/
/*--------------------------------------------------------------------------*/

int main( int argc , char ** argv )
{
 if( argc > 1 )
  instance = argv[ 1 ];

 auto d = read( instance );
 std::cout << d.demand.size() << " scenarios, " << d.demand[ 0 ].size()
	   << " time steps, design cost " << d.cost << " in [ " << d.lb
	   << " , " << d.ub << " ]" << std::endl;

 double t_ref , t_bds = 0 , t_inv = 0;
 long i_ref , i_bds = 0 , i_inv = 0;
 int s_ref , s_bds = 0 , s_inv = 0;
 bool has_bds , has_inv;

 auto mono = build_monolithic( d );
 const auto ref = solve( mono , "BDSSCfg.txt" , t_ref , i_ref , s_ref );

 auto bend = build_structured( d );
 const auto bds = solve_by_bundle( bend , "BSPar-BDS.txt" , t_bds ,
				   i_bds , s_bds , has_bds , ref );

 const auto inv = solve_ad_hoc( t_inv , i_inv , s_inv , has_inv , ref );

 std::cout.precision( 12 );
 std::cout << "\nmonolithic        " << ref << "  status " << s_ref
	   << "  in " << t_ref << " s";
 if( has_bds )
  std::cout << "\nBenders, generic  " << bds << "  status " << s_bds
	    << "  in " << t_bds << " s, " << i_bds << " iterations";
 if( has_inv )
  std::cout << "\nBenders, ad hoc   " << inv << "  status " << s_inv
	    << "  in " << t_inv << " s, " << i_inv << " iterations";
 std::cout << std::endl;

 const auto scale = std::max( 1.0 , std::abs( ref ) );
 bool ok = true;

 if( std::abs( bds - ref ) > tolerance * scale ) {
  std::cout << "ERROR: the generic Benders form is off by "
	    << std::abs( bds - ref ) / scale << std::endl;
  ok = false;
  }

 if( std::abs( inv - ref ) > tolerance * scale ) {
  std::cout << "ERROR: the ad hoc Benders form is off by "
	    << std::abs( inv - ref ) / scale << std::endl;
  ok = false;
  }

 if( ! ok )
  std::cout << "\nthe forms do not describe the same problem" << std::endl;
 else
  if( has_bds && has_inv )
   std::cout << "\nOK, the two forms agree with the extensive optimum"
	     << std::endl;
  else
   std::cout << "\nOK, what could be run agrees with the extensive optimum"
	     << std::endl;

 delete mono;
 delete bend;
 return( ok ? 0 : 1 );
 }

/*--------------------------------------------------------------------------*/
/*------------------------- End File test_bds.cpp -------------------------*/
/*--------------------------------------------------------------------------*/