From 37a18797d235d7a973720bbdcc4dcb327e32c58c Mon Sep 17 00:00:00 2001 From: Anirudh Subramanyam Date: Mon, 14 Oct 2019 23:28:17 -0500 Subject: [PATCH] initial commit --- .gitignore | 47 + CMakeLists.txt | 66 + LICENSE | 11 + README.md | 55 + inc/Constants.h | 30 + inc/constraintExpr.hpp | 988 +++++++++++ inc/indexingTools.h | 182 ++ inc/instance_knp.hpp | 166 ++ inc/instance_psp.hpp | 109 ++ inc/instance_spp.hpp | 134 ++ inc/problemInfo.hpp | 317 ++++ inc/problemInfo_knp.hpp | 72 + inc/problemInfo_psp.hpp | 72 + inc/problemInfo_spp.hpp | 72 + inc/robustSolver.hpp | 502 ++++++ inc/uncertainty.hpp | 260 +++ src/indexingTools.cpp | 529 ++++++ src/problemInfo.cpp | 108 ++ src/problemInfo_knp.cpp | 278 +++ src/problemInfo_psp.cpp | 400 +++++ src/problemInfo_spp.cpp | 210 +++ src/robustSolver.cpp | 3644 +++++++++++++++++++++++++++++++++++++++ src/test.cpp | 74 + src/uncertainty.cpp | 452 +++++ 24 files changed, 8778 insertions(+) create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 LICENSE create mode 100644 README.md create mode 100644 inc/Constants.h create mode 100644 inc/constraintExpr.hpp create mode 100644 inc/indexingTools.h create mode 100644 inc/instance_knp.hpp create mode 100644 inc/instance_psp.hpp create mode 100644 inc/instance_spp.hpp create mode 100644 inc/problemInfo.hpp create mode 100644 inc/problemInfo_knp.hpp create mode 100644 inc/problemInfo_psp.hpp create mode 100644 inc/problemInfo_spp.hpp create mode 100644 inc/robustSolver.hpp create mode 100644 inc/uncertainty.hpp create mode 100644 src/indexingTools.cpp create mode 100644 src/problemInfo.cpp create mode 100644 src/problemInfo_knp.cpp create mode 100644 src/problemInfo_psp.cpp create mode 100644 src/problemInfo_spp.cpp create mode 100644 src/robustSolver.cpp create mode 100644 src/test.cpp create mode 100644 src/uncertainty.cpp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f68a5cb --- /dev/null +++ b/.gitignore @@ -0,0 +1,47 @@ +# Prerequisites +*.d + +# Compiled Object files +*.slo +*.lo +*.o +*.obj + +# Precompiled Headers +*.gch +*.pch + +# Compiled Dynamic libraries +*.so +*.dylib +*.dll + +# Fortran module files +*.mod +*.smod + +# Compiled Static libraries +*.lai +*.la +*.a +*.lib + +# Executables +*.exe +*.out +*.app + +# do track lib.a, even though you're ignoring .a files above +# !lib.a + +# only ignore the TODO file in the current directory, not subdir/TODO +# /TODO + +# ignore all files in the build/ directory +build/ + +# ignore doc/notes.txt, but not doc/server/arch.txt +# doc/*.txt + +# ignore all .pdf files in the doc/ directory +# doc/**/*.pdf \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..ce3d1e8 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,66 @@ +# CHECK CMAKE VERSION +cmake_minimum_required(VERSION 2.8) + +# SPECIFY PROJECT NAME +project(K_ADAPTABILITY) + +# SPECIFY SOURCE FILES (.CPP, .C) +set(SOURCE_FILES + ./src/uncertainty.cpp + ./src/problemInfo.cpp + ./src/problemInfo_knp.cpp + ./src/problemInfo_spp.cpp + ./src/problemInfo_psp.cpp + ./src/indexingTools.cpp + ./src/robustSolver.cpp + ./src/test.cpp +) + +# SPECIFY CUSTOM HEADER FILE DIRECTORIES (OR SET "DEFAULT" IF NONE) +set(HDRDIR ./inc) + +# SPECIFY NAME OF EXECUTABLE +set(EXEC kadaptability) + +# SET COMPILER DIRECTORY (OR SET "DEFAULT") +set(COMPILER_DIR /usr/local/bin/g++) + +# SET CPLEX DIRECTORY +set(CPLEX_DIR /opt/ibm/ILOG/CPLEX_Studio127) + +# CHECK OS +if(APPLE) + set(mySystem "osx") +elseif(UNIX) + set(mySystem "linux") +else() + message(FATAL_ERROR "ONLY UNIX IS SUPPORTED AT THIS TIME") +endif() + +# SET COMPILER DIRECTORY +if(NOT "${COMPILER_DIR}" STREQUAL "DEFAULT") + set(CMAKE_CXX_COMPILER ${COMPILER_DIR}) +endif() + +# SET COMPILER FLAGS +set(COMP_FLAGS_DEBUG " -g -std=c++14 -m64 -O0 -DIL_STD -pthread -D_GLIBCXX_DEBUG -Wall -Wextra -fPIC -fno-strict-aliasing -fexceptions") +set(COMP_FLAGS_RELEASE " -std=c++14 -m64 -O2 -DIL_STD -pthread -DNDEBUG -Wall -Wextra -fPIC -fno-strict-aliasing -fexceptions") +set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} ${COMP_FLAGS_DEBUG}") +set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} ${COMP_FLAGS_RELEASE}") + +# SET INCLUDE DIRECTORIES (.h .hpp) +include_directories(${EXEC} ${CPLEX_DIR}/cplex/include ${CPLEX_DIR}/concert/include) + +# INCLUDE THE CUSTOM INC DIRECTORIES IF DEFINED +if(NOT "${HDRDIR}" STREQUAL "DEFAULT") + include_directories(${EXEX} ${HDRDIR}) +endif() + +# SET LINK DIRECTORIES +link_directories(${CPLEX_DIR}/cplex/lib/x86-64_${mySystem}/static_pic ${CPLEX_DIR}/concert/lib/x86-64_${mySystem}/static_pic) + +# CREATE THE TARGET (executable or library) +add_executable(${EXEC} ${SOURCE_FILES}) + +# LINK THE DESIRED LIBRARIES +target_link_libraries(${EXEC} concert ilocplex cplex m pthread) diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..b23ff04 --- /dev/null +++ b/LICENSE @@ -0,0 +1,11 @@ +Copyright 2019 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann + +Anirudh Subramanyam and Chrysanthos Gounaris are affiliated with Carnegie Mellon University. Wolfram Wiesemann is affiliated with Imperial College London. + +Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..8ebf5ee --- /dev/null +++ b/README.md @@ -0,0 +1,55 @@ + +k-adaptability-solver +==== + +k-adaptability-solver is a numerical optimization package, written in C++, for solving K-adaptability counterparts of two-stage mixed-integer robust optimization problems, based on our paper [K-Adaptability in Two-Stage Mixed-Integer Robust Optimization](https://arxiv.org/abs/1706.07097). + +If you use this code in your research, please cite: +``` +@article{sgw:19, + title = {K-adaptability in two-stage mixed-integer robust optimization}, + author = {Subramanyam, Anirudh and Gounaris, Chrysanthos E and Wiesemann, Wolfram}, + journal = {Mathematical Programming Computation}, + year = {2019} +} +``` +## Requirements +* Cmake 2.8 or higher +* CPLEX 12.6 or higher +* OS X or Linux OS + +## Installation on OS X and Linux +* Obtain a copy of the CPLEX software and license [here](https://www.ibm.com/analytics/cplex-optimizer) +* Specify the CPLEX installation directory in CMakeLists.txt +* Enter the following commands on the terminal + ``` + $ mkdir build + $ cd build/ + $ cmake Unix Makefiles -DCMAKE_BUILD_TYPE=Release .. + $ make + ``` +* You should find an executable `kadaptability` in the build/ directory + +## Note +* You may also compile in debug mode by changing `Release` to `Debug` in the above command. +* You may change the `TIME_LIMIT`, `MEMORY_LIMIT` and `NUM_THREADS` parameters in inc/Constants.h. If you do modify the defaults, please re-build for the changes to take effect. + +## Usage +The problem (variables, objective function, constraints, uncertainty set etc) is modeled in the abstract class `KAdaptableInfo` defined in problemInfo.hpp/cpp. + +To solve a new problem, you must create a new class which will inherit from this abstract class. +You must then override the following functions: + +* `makeVars` +* `makeUncSet` +* `makeConsX` +* `makeConsY` +* `create` +* `clone` + +As an example, have a look at problemInfo_knp.hpp/cpp. This file implements examples 4.2 and 4.3 from the above paper. + +After you have implemented the above class, you can solve it by creating a `KAdaptableSolver` object and calling the function `solve_KAdaptability(K, heuristic, x)`. The solution vector is returned in `x`, and a summary of the run will be written to standard output. + +As an example, please follow the `main` function defined in src/test.cpp. + diff --git a/inc/Constants.h b/inc/Constants.h new file mode 100644 index 0000000..60bd4bd --- /dev/null +++ b/inc/Constants.h @@ -0,0 +1,30 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + +#ifndef CONSTANTS_H_ +#define CONSTANTS_H_ + +// Exception numbers +#define EXCEPTION_CPXINIT 1 // CPLEX could not be initialized +#define EXCEPTION_CPXEXIT 2 // CPLEX could not be closed +#define EXCEPTION_CPXNEWCOLS 3 // Error while trying to create new columns +#define EXCEPTION_CPXNEWROWS 4 // Error while trying to create new rows +#define EXCEPTION_X 5 // Invalid size of solution vector +#define EXCEPTION_K 6 // Invalid policy index +#define EXCEPTION_Q 7 // Invalid size of uncertain parameters + +// Defaults +#define TIME_LIMIT 7200 +#define MEMORY_LIMIT 3072 +#define NUM_THREADS 1 + +#endif diff --git a/inc/constraintExpr.hpp b/inc/constraintExpr.hpp new file mode 100644 index 0000000..3a11198 --- /dev/null +++ b/inc/constraintExpr.hpp @@ -0,0 +1,988 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#ifndef CONSTRAINT_EXPRESSION_HPP +#define CONSTRAINT_EXPRESSION_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "uncertainty.hpp" + + +#define CONSTRAINT_EXPRESSION_OUTPUTLEVEL 0 + + +typedef std::vector indices_t; +typedef std::vector coefficients_t; +typedef std::pair pair_t; +typedef std::vector pairIndices_t; + + + + +/*--------------------------------------------------------- + This class supports constraints of the following type: + + cx + bq + xSq (sense) d + + x are decision variables (varIndices) + c are const doubles with same size as x (varCoeffs) + q are uncertain parameters (paramIndices) + b are const doubles with same size as q (paramCoeffs) + S is a (nx times nq) matrix (sMatrix) + d is the right hand side double value (rhs) + + *---------------------------------------------------------*/ +class ConstraintExpression { + friend class KAdaptableExpression; +private: + char sense;/* Default: 'L', i.e., <= */ + double rhs;/* Default: 0 */ + std::string name; + indices_t varIndices; + coefficients_t varCoeffs; + + + /* Following are required only for Robust counterparts */ + pairIndices_t bilinearIndices; /* pairs of (indices of) bilinear terms (xi,qj), i.e., products of decision variables with uncertain parameters */ + coefficients_t bilinearCoeffs; /* S(xi,qj) := coefficient of term xi*qj */ + indices_t paramIndices; + coefficients_t paramCoeffs; + + inline bool isConsistentWithDesign() const { + if (std::set(varIndices.begin(), varIndices.end()).size() != varIndices.size()) return false; + if (std::set(paramIndices.begin(), paramIndices.end()).size() != paramIndices.size()) return false; + for (unsigned i = 0; i < bilinearIndices.size(); ++i) + for (unsigned j = i + 1; j < bilinearIndices.size(); ++j) + if (bilinearIndices[i] == bilinearIndices[j])return false; + if (varIndices.size() != varCoeffs.size()) return false; + if (paramIndices.size() != paramCoeffs.size()) return false; + if (bilinearIndices.size() != bilinearCoeffs.size()) return false; + if (sense != 'L' && sense != 'E' && sense != 'G') return false; + /* For every bilinear term xi*qj, ci*xi and bj*qj must exist (possibly with ci=0,bj=0) */ + /* Inverse need not be true */ + for (const auto& xq : bilinearIndices) { + if (std::find(varIndices.begin(), varIndices.end(), xq.first) == varIndices.end()) return false; + if (std::find(paramIndices.begin(), paramIndices.end(), xq.second) == paramIndices.end()) return false; + } + return true; + } + /* Add the term cx (vec1 = varIndices, coeffs = paramCoeffs) + * or bq (vec1 = paramIndices, coeffs = paramCoeffs) + */ + inline void addTerm(const int ind, const double val, indices_t& vec1, coefficients_t& coeffs) { + + if (val == 0.0) { if(CONSTRAINT_EXPRESSION_OUTPUTLEVEL) std::cerr << "Warning: Attempting to add term with zero coefficient. Ignoring. \n"; return; } + + const int Pos = std::find(vec1.begin(), vec1.end(), ind) - vec1.begin(); + + if (Pos >= (int)vec1.size()) { + vec1.push_back(ind); + coeffs.push_back(val); + } + else { + assert(vec1.at(Pos) == ind); + coeffs.at(Pos) += val; + } + + assert(isConsistentWithDesign()); + return; + } + inline void reserve(size_t n) { + varIndices.reserve(n); + varCoeffs.reserve(n); + paramIndices.reserve(n); + paramCoeffs.reserve(n); + bilinearIndices.reserve(n); + bilinearCoeffs.reserve(n); + + return; + } + +public: + ConstraintExpression() : sense('L'), rhs(0) { reserve(1000); } + ConstraintExpression(const std::string& argName) : sense('L'), rhs(0), name(argName) { reserve(1000); } + ConstraintExpression(const std::string& argName, const char argSense, const double argRhs) + : sense(argSense), rhs(argRhs), name(argName) { reserve(1000); assert(isConsistentWithDesign()); } + ~ConstraintExpression() {/**/} + ConstraintExpression (const ConstraintExpression&) = default; + ConstraintExpression (ConstraintExpression&&) = default; + ConstraintExpression& operator=(const ConstraintExpression&) = default; + ConstraintExpression& operator=(ConstraintExpression&&) = default; + + inline void sign(const char argSense) { sense = argSense; assert(isConsistentWithDesign()); return; } + inline void RHS(const double argRhs) { rhs = argRhs; return; } + inline void rowname(const std::string& argName) { name = argName; return; } + inline void addTermX(const int ind, const double val) { addTerm(ind, val, varIndices, varCoeffs); return; } + inline void addTermQ(const int ind, const double val) { addTerm(ind, val, paramIndices, paramCoeffs); return; } + inline void clear() { + sense = 'L'; + rhs = 0; + name.clear(); + varIndices.clear(); + varCoeffs.clear(); + bilinearIndices.clear(); + bilinearCoeffs.clear(); + paramIndices.clear(); + paramCoeffs.clear(); + return; + } + /* Add the term xSq */ + inline void addTermProduct(const int xind, const int qind, const double coeff = 1.0) { + if (coeff == 0.0) { if(CONSTRAINT_EXPRESSION_OUTPUTLEVEL) std::cerr << "Warning: Attempting to add term with zero coefficient. Ignoring. \n"; return; } + + /* Add x to cx with c = 0 if it does not already exist */ + if (std::find(varIndices.begin(), varIndices.end(), xind) == varIndices.end()) { + varIndices.push_back(xind); + varCoeffs.push_back(0); + } + /* Add q to bq with b = 0 if it does not already exist */ + if (std::find(paramIndices.begin(), paramIndices.end(), qind) == paramIndices.end()) { + paramIndices.push_back(qind); + paramCoeffs.push_back(0); + } + + const pair_t xq(xind, qind); + /* Add xi*qj with coefficient "coeff" if this bilinear term does not already exist */ + const int curPos = std::find(bilinearIndices.begin(), bilinearIndices.end(), xq) - bilinearIndices.begin(); + if (curPos >= (int)bilinearIndices.size()) { + bilinearIndices.emplace_back(xq); + bilinearCoeffs.emplace_back(coeff); + } + /* Else, update S(xi,qj) with "coeff" */ + else { + assert(bilinearIndices.at(curPos) == xq); + bilinearCoeffs.at(curPos) += coeff; + } + + assert(isConsistentWithDesign()); + return; + } + + + + inline bool isEmpty() const { + for (const auto& i : varCoeffs) if (i != 0.0) return false; + for (const auto& i : bilinearCoeffs) if (i != 0.0) return false; + return true; + } + inline bool existConstQTerms() const { + for (const auto& i : paramCoeffs) if (i != 0.0) return true; + return false; + } + inline bool existBilinearTerms() const { + for (const auto& i : bilinearCoeffs) if (i != 0.0) return true; + return false; + } + inline std::vector getVarNames(CPXCENVptr env = nullptr, CPXCLPptr lp = nullptr) const { + + std::vector varNames; + for (unsigned i = 0; i < varIndices.size(); i++) { + if (env == nullptr || lp == nullptr) + varNames.emplace_back("x(" + std::to_string(varIndices[i]) + ")"); + else { + CPXSIZE surplus = 0; CPXXgetcolname(env, lp, nullptr, nullptr, 0, &surplus, varIndices[i], varIndices[i]); + CPXSIZE colnamespace = -surplus; + if (colnamespace > 0) { + char **colname = (char **)malloc(sizeof(char *)); + char *colnamestore = (char *) malloc(colnamespace); + CPXXgetcolname(env, lp, colname, colnamestore, colnamespace, &surplus, varIndices[i], varIndices[i]); + + varNames.emplace_back(colname[0], std::find(colname[0], colname[0] + strlen(colname[0]), '\0')); + if (colname != NULL) { free(colname); colname = nullptr; } + if (colnamestore != NULL) { free(colnamestore); colnamestore = nullptr; } + } + else varNames.emplace_back("xFake(" + std::to_string(varIndices[i]) + ")"); + } + } + + return varNames; + } + inline void print(CPXCENVptr env = nullptr, CPXCLPptr lp = nullptr, std::ostream& out = std::cout) const { + if (isEmpty()) { out << "(empty expression)\n"; return; } + + /* Determine variable names */ + std::vector varNames = getVarNames(env, lp); + + + if (!name.empty()) out << "Printing " << name << "-> "; + for (unsigned i = 0; i < varCoeffs.size(); i++) + if(varCoeffs[i] != 0.0) + out << ( (varCoeffs[i] > 0) ? "+" : "-") << std::abs(varCoeffs[i]) << " " << varNames[i] << " " ; + for (unsigned i = 0; i < paramCoeffs.size(); i++) + if (paramCoeffs[i] != 0.0) + out << ((paramCoeffs[i] > 0) ? "+" : "-") << std::abs(paramCoeffs[i]) << " q(" << paramIndices[i] << ") "; + for (unsigned i = 0; i < bilinearCoeffs.size(); i++) + if (bilinearCoeffs[i] != 0.0) { + out << ((bilinearCoeffs[i] > 0) ? "+" : "-") << std::abs(bilinearCoeffs[i]) << " "; + out << varNames[std::find(varIndices.begin(), varIndices.end(), bilinearIndices[i].first) - varIndices.begin()] << " q(" << bilinearIndices[i].second << ") "; + } + out << (sense=='E' ? "=" : (sense=='G' ? ">=" : "<=")) << " " << rhs << "\n"; + + return; + } + inline void getDeterministicConstraint( const std::vector& paramValues, + CPXNNZ& nzcnt, + double& trueRhs, + char& trueSense, + std::vector& rmatind, + std::vector& rmatval) const + { + trueSense = sense; + + /* Compute right-hand-size as supplied rhs - bq */ + trueRhs = rhs; + for (unsigned i = 0; i < paramCoeffs.size(); i++) + if (paramCoeffs.at(i) != 0.0) + trueRhs -= (paramCoeffs.at(i)*paramValues.at(paramIndices.at(i))); + + /* Compute non-zeros, i.e., coefficients of decision variables */ + nzcnt = 0; + rmatind.clear(); + rmatval.clear(); + for (unsigned i = 0; i < varIndices.size(); ++i) { + bool inExpr = 0; /* Does this decision variable appear in this expression? */ + double coef = 0.0;/* It's coefficient should be non-zero */ + + /* Check terms involving x only: cx */ + if (varCoeffs[i] != 0.0) { + coef += varCoeffs[i]; + inExpr = 1; + } + /* Check product terms: xSq */ + for (unsigned k = 0; k < bilinearIndices.size(); ++k) if (bilinearIndices[k].first == varIndices[i]) { + if (bilinearCoeffs[k] != 0.0) { + coef += paramValues.at(bilinearIndices[k].second)*bilinearCoeffs[k]; + inExpr = 1; + } + } + + /* This variable must appear somewhere */ + if (inExpr) { + nzcnt++; + rmatind.emplace_back(varIndices[i]); + rmatval.emplace_back(coef); + } + else if(CONSTRAINT_EXPRESSION_OUTPUTLEVEL) + std::cerr << "Warning: Expression contains empty variable. Ignoring. \n"; + } + assert((CPXNNZ)rmatind.size() == nzcnt); + assert(rmatval.size() == rmatind.size()); + + return; + } + inline void getStochasticConstraint(const std::vector& varValues, + CPXNNZ& nzcnt, + double& trueRhs, + char& trueSense, + std::vector& rmatind, + std::vector& rmatval) const + { + trueSense = sense; + + /* Compute right-hand-size as supplied rhs - cx */ + trueRhs = rhs; + for (unsigned i = 0; i < varCoeffs.size(); i++) + if (varCoeffs.at(i) != 0.0) + trueRhs -= (varCoeffs.at(i)*varValues.at(varIndices.at(i))); + + /* Compute non-zeros, i.e., coefficients of decision variables */ + nzcnt = 0; + rmatind.clear(); + rmatval.clear(); + for (unsigned i = 0; i < paramIndices.size(); ++i) { + bool inExpr = 0; /* Does this parameter appear in this expression? */ + double coef = 0.0;/* It's coefficient should be non-zero */ + + /* Check terms involving q only: bq */ + if (paramCoeffs[i] != 0.0) { + coef += paramCoeffs[i]; + inExpr = 1; + } + /* Check product terms: xSq */ + for (unsigned k = 0; k < bilinearIndices.size(); ++k) if (bilinearIndices[k].second == paramIndices[i]) { + if (bilinearCoeffs[k] != 0.0) { + coef += varValues.at(bilinearIndices[k].first)*bilinearCoeffs[k]; + inExpr = 1; + } + } + + /* This variable must appear somewhere */ + if (inExpr) { + nzcnt++; + rmatind.emplace_back(paramIndices[i]); + rmatval.emplace_back(coef); + } + else if(CONSTRAINT_EXPRESSION_OUTPUTLEVEL) + std::cerr << "Warning: Expression contains empty uncertain parameter. Ignoring. \n"; + } + assert((CPXNNZ)rmatind.size() == nzcnt); + assert(rmatval.size() == rmatind.size()); + + return; + } + inline int addToCplex(CPXCENVptr env, + CPXLPptr lp, + UNCSetCPtr UncSet = nullptr, + const bool isDualized = false, + const std::vector& paramValues = {}) const + { + if (isEmpty()) { + if(CONSTRAINT_EXPRESSION_OUTPUTLEVEL) std::cerr << " Warning: Attempting to add empty expression to CPLEX. Ignoring. \n"; + return 0; + } + + /* Do the uncertain parameters appear as coefficients of decision variables? */ + const bool qTermsInProduct = existBilinearTerms(); + /* Do the uncertain parameters appear independently by themselves? */ + const bool qTermsConst = existConstQTerms(); + /* Is this a "deterministic" constraint that does not require any dualization? */ + bool deterministic = !isDualized; + + if (isDualized) if (!qTermsInProduct) if (!qTermsConst) deterministic = true; + if (deterministic) if (qTermsInProduct || qTermsConst) if (paramValues.empty()) { std::cerr << " Error: Parameter values not supplied for deterministic model. \n"; return 1; } + + + + + //===================================================================== + /* THIS CONSTRAINT DOES NOT REQUIRE REFORMULATION IN THE MAIN MODEL */ + //===================================================================== + if (deterministic) { + char trueSense = sense; + const CPXNNZ rmatbeg = 0; + const char *rowname = name.c_str(); + + /* Compute right-hand-size as supplied rhs - bq */ + double trueRhs = rhs; + for (unsigned i = 0; i < paramCoeffs.size(); i++) + if (paramCoeffs.at(i) != 0.0) + trueRhs -= (paramCoeffs.at(i)*paramValues.at(paramIndices.at(i))); + + /* Compute non-zeros, i.e., coefficients of decision variables */ + CPXNNZ nzcnt = 0; + std::vector rmatind; + std::vector rmatval; + for (unsigned i = 0; i < varIndices.size(); ++i) { + bool inExpr = 0; /* Does this decision variable appear in this expression? */ + double coef = 0.0;/* It's coefficient should be non-zero */ + + /* Check terms involving x only: cx */ + if (varCoeffs[i] != 0.0) { + coef += varCoeffs[i]; + inExpr = 1; + } + /* Check product terms: xSq */ + for (unsigned k = 0; k < bilinearIndices.size(); ++k) if (bilinearIndices[k].first == varIndices[i]) { + if (bilinearCoeffs[k] != 0.0) { + coef += paramValues.at(bilinearIndices[k].second)*bilinearCoeffs[k]; + inExpr = 1; + } + } + + /* This variable must appear somewhere */ + if (inExpr) { + nzcnt++; + rmatind.emplace_back(varIndices[i]); + rmatval.emplace_back(coef); + } + else if(CONSTRAINT_EXPRESSION_OUTPUTLEVEL) + std::cerr << "Warning: Expression contains empty variable. Ignoring. \n"; + } + assert((CPXNNZ)rmatind.size() == nzcnt); + assert(rmatval.size() == rmatind.size()); + + return CPXXaddrows(env, lp, 0, 1, nzcnt, &trueRhs, &trueSense, &rmatbeg, &rmatind[0], &rmatval[0], nullptr, &rowname); + } + //===================================================================== + /* DUALIZE CONSTRAINT BEFORE ADDING BACK TO THE (REFORMULATED) MODEL */ + //===================================================================== + else { + if (UncSet == nullptr) { + std::cerr << " Error: Null pointer supplied for uncertainty set in dualized constraint. \n"; + return 1; + } + if (!UncSet->isUncertain()) { + std::cerr << " Error: Attempting to dualize constraint for deterministic problem. \n"; + return 1; + } + if (!qTermsInProduct) if (qTermsConst) { /* "Effectively" deterministic constraint */ + if (sense == 'E') { + std::cerr << "Error: Attempting to dualize equality constraint in which uncertain parameters do not appear as coefficients of decision variables. \n"; + return 1; + } + //if(CONSTRAINT_EXPRESSION_OUTPUTLEVEL) std::cerr << " Warning: Not dualizing constraint because uncertain parameters do not appear as coefficients of decision variables. \n"; + + /* If this is a >= constraint, move the q_terms to the right hand side and then maximize */ + std::vector vecCoeffs(paramCoeffs); if (sense == 'G') for (auto& qj : vecCoeffs) qj *= (-1.0); + + /* Obtain maximum value of the q_terms over the uncertainty set */ + const double max_qTermsConst = UncSet->getMaximumValue(paramIndices, vecCoeffs); + const double NewRhs = rhs + ((sense == 'G') ? +max_qTermsConst : -max_qTermsConst); + + /* Create a new constraint which is a copy of *this but with modified RHS + * and add that constraint deterministically + */ + ConstraintExpression NewExpr(name, sense, NewRhs); + for (unsigned i = 0; i < varIndices.size(); ++i) + NewExpr.addTermX(varIndices[i], varCoeffs[i]); + return NewExpr.addToCplex(env, lp); + } + + /*************************************************************/ + /* Before you add the "single" dualized constraint, + * create as many new dual variables as there are rows in your + * "W matrix" of the uncertainty set + */ + /*************************************************************/ + assert(qTermsInProduct); + + + + int status = 0; /* return value */ + + if (sense == 'E') { + std::cerr << "Error!!!! To-do!! Coefficient matching in equality constraints!! \n"; + exit(100); + } + + /* r = # of new (dual) variables created on top of existing model */ + const int r = UncSet->addVariables_DualVars(env, lp, name); + if (r == 0) { std::cerr << "Error: No dual variables created while reformulating/dualizing constraints.\n"; return 1; } + + /* UncSet: W*q + V*Ksee [sense] h */ + /* Depending on sense [>=, <=, =], the dual variables have appropriate signs and bounds */ + const auto W = UncSet->getMatrixW(); + const auto V = UncSet->getMatrixV(); + const auto h = UncSet->getMatrixH(); + assert(W.size() == h.size()); + assert(V.size() == h.size()); + assert(r + 1 == (int)h.size()); + + /* Last variable index of the original model + * before new variables were created + */ + const int lastIndBeforeDual = CPXXgetnumcols(env, lp) - r - 1; + + + + + /****/ + /****/ + + + + /* Add the constraint corresponding to the objective of the inner maximization problem */ + ConstraintExpression ReformObjExpr("DualObj(" + name + ")", sense, rhs); + + /* Terms involving only x (e.g., cx) carry over as-is */ + for (unsigned i = 0; i < varIndices.size(); ++i) if(varCoeffs[i] != 0.0) ReformObjExpr.addTermX(varIndices[i], varCoeffs[i]); + + /* Add objective function expression of the dual of the inner problem + * Only add those coefficients which are true non-zeros + * SUM[s = 1 to r, h(s)Ksee(s)] + * Note: if the original constraint is a >= constraint, then multiply above term by (-1) + */ + for (int s = 1; s <= r; s++) if (h.at(s) != 0.0) { + const int dualVarIndex = lastIndBeforeDual + s; + ReformObjExpr.addTermX(dualVarIndex, ((sense == 'G') ? -1 : +1)*h.at(s)); + } + + /* Add constraint */ + status += ReformObjExpr.addToCplex(env, lp); + + + + /****/ + /****/ + + + + /* Add the constraints of the dual problem */ + /* Assuming that rows and columns of W, V and h are 1-indexed */ + for (int ll = 1; ll <= UncSet->getNoOfUncertainParameters(); ll++) { + const int l = UncSet->getParamIndex('q', ll); + + /* To get the rhs of the constraint of the dual problem, compute the coefficient of q(l) in the original constraint */ + /* This coefficient has an x component (coming from xSq) and a constant component (coming from bq) */ + /* The "rhs" (which is a constant) comes from the constant component bq */ + double dualConstraintRhs = 0; + + /* If b(l)*q(l) is in this expression, then rhs = b(l) if <=, -b(l) if >= */ + const int lPos = std::find(paramIndices.begin(), paramIndices.end(), l) - paramIndices.begin(); + if (lPos < (int)paramIndices.size()) dualConstraintRhs += paramCoeffs[lPos]; + if (sense == 'G') dualConstraintRhs *= (-1.0); + + /* Constraint of dual problem corresponding to this column of the uncertainty set */ + ConstraintExpression ReformConstrExpr("DualCon(" + name + "," + std::to_string(l) + ")"); + ReformConstrExpr.sign('E'); + ReformConstrExpr.RHS(dualConstraintRhs); + + /* Add terms involving products of q(l) with x(i), i.e., xSq */ + for (unsigned i = 0; i < bilinearIndices.size(); i++) if (bilinearIndices[i].second == l) { + assert(lPos < (int)paramIndices.size()); + if (bilinearCoeffs[i] != 0.0) + ReformConstrExpr.addTermX(bilinearIndices[i].first, ((sense == 'G') ? +1 : -1)*bilinearCoeffs[i]); + } + + /* Add nonzeros corresponding to constraints of the uncertainty set */ + /* Only add those coefficients which are true non-zeros */ + /* SUM[s = 1 to r, W(s,l)Ksee(s)] */ + for (int s = 1; s <= r; s++) if (W.at(s).at(l) != 0.0) { + const int dualVarIndex = lastIndBeforeDual + s; + ReformConstrExpr.addTermX(dualVarIndex, W.at(s).at(l)); + } + + /* Add constraint */ + status += ReformConstrExpr.addToCplex(env, lp); + } + + + /****/ + /****/ + + + + /* If using factor models, then we need an additional constraint + * Note that for any other uncertainty set, V = 0 and we'll end up adding no new constraints + */ + for (int f = 1; f <= UncSet->getNumBudgetsFactors(); f++) { + bool nonzero = 0; + for (int s = 1; s <= r; s++) if (V.at(s).at(f) != 0.0) { nonzero = 1; break; } + + + /* Attempt to add the constraint only if at least one non-zero exists */ + if (nonzero) { + ConstraintExpression ReformFactorExpr("DualConFactor(" + name + ",F" + std::to_string(f) + ")", 'E', 0); + for (int s = 1; s <= r; s++) if (V.at(s).at(f) != 0.0) { + const int dualVarIndex = lastIndBeforeDual + s; + ReformFactorExpr.addTermX(dualVarIndex, V.at(s).at(f)); + } + + /* Add constraint */ + status += ReformFactorExpr.addToCplex(env, lp); + } + + /****/ + /****/ + } + + return status; + } + } + inline char getSense() const { return sense; } + inline double getRHS() const { return rhs; } + inline std::string getName() const { return name; } + inline double getViolation_fixedXQ(const std::vector& varValues, + const std::vector& paramValues = {}, + const bool normalized = 0) const + { + double viol = 0; + + const bool xTerms = !isEmpty(); + const bool qTermsInProduct = existBilinearTerms(); + const bool qTermsConst = existConstQTerms(); + + + if (xTerms) { + if(varValues.empty()) { + std::cerr << "Error: getViolation() requires values for primal variables\n"; + return 1.0; + } + } + if (qTermsInProduct || qTermsConst) { + if (paramValues.empty()) { + std::cerr << "Error: getViolation() requires values for uncertain parameters\n"; + return 1.0; + } + } + else if (!paramValues.empty()) { + if (CONSTRAINT_EXPRESSION_OUTPUTLEVEL) + std::cerr << "Warning: Ignoring supplied parameter values in getViolation()\n"; + } + + // normalizing factor + double f = (normalized ? 0.0 : 1.0); + + // -d + viol += -rhs; + + // +bq* + if (qTermsConst) { + for (unsigned i = 0; i < paramCoeffs.size(); i++) { + viol += paramCoeffs[i] * paramValues.at(paramIndices.at(i)); + if (normalized) f += paramCoeffs[i] * paramCoeffs[i]; + } + } + + // +cx* + if (xTerms) { + for (unsigned i = 0; i < varCoeffs.size(); i++) { + viol += varCoeffs[i] * varValues.at(varIndices.at(i)); + } + } + + // +x*Sq* + if (qTermsInProduct) { + for (unsigned i = 0; i < bilinearCoeffs.size(); i++) { + viol += bilinearCoeffs[i] * varValues.at(bilinearIndices.at(i).first) * paramValues.at(bilinearIndices.at(i).second); + if (normalized) f += bilinearCoeffs[i] * varValues.at(bilinearIndices.at(i).first) * bilinearCoeffs[i] * varValues.at(bilinearIndices.at(i).first); + } + } + + // normalizing factor + if (normalized) { + assert(f != 0.0); + viol /= std::sqrt(f); + } + + // Make changes for >= , = + switch (sense) { + case 'G': viol *= -1.0; break; + case 'E': if (viol < 0) viol *= -1.0; break; + } + + return viol; + } + inline ConstraintExpression flip0() const { + ConstraintExpression N; + switch (sense) { + case 'G': N.sense = 'L';break; + case 'L': N.sense = 'G';break; + case 'E': N.sense = 'E';break; + } + N.rhs = rhs; + N.name = name; + N.varIndices = paramIndices; + N.varCoeffs = paramCoeffs; + N.paramIndices = varIndices; + N.paramCoeffs = varCoeffs; + N.bilinearCoeffs = bilinearCoeffs; + N.bilinearIndices = bilinearIndices; + for (unsigned int i = 0; i < bilinearIndices.size(); i++) { + N.bilinearIndices[i].first = bilinearIndices[i].second; + N.bilinearIndices[i].second = bilinearIndices[i].first; + } + assert(N.isConsistentWithDesign()); + + const int curPos = std::find(N.paramIndices.begin(), N.paramIndices.end(), 0) - N.paramIndices.begin(); + assert(curPos < (int)N.paramIndices.size()); + if (curPos < (int)N.paramIndices.size()) { + assert(N.paramCoeffs[curPos] == 1); + N.addTermX(0, N.paramCoeffs[curPos]); + N.paramIndices.erase(N.paramIndices.begin() + curPos); + N.paramCoeffs.erase(N.paramCoeffs.begin() + curPos); + } + assert(N.isConsistentWithDesign()); + return N; + } +}; + + + + + + + + +/*---------------------------------------------------------------------- + This class supports constraints of the following type: + + max_q min_p {cx(p) + bq + x(p)Sq - d(p)} <= rhs + <=> max_{q,t} t <= rhs + s.t. t <= cx(p) + bq + x(p)Sq - d(p) \forall p + q \in UncSet + +[ cx(p) + bq(p) + x(p)Sq <= d(p) ] are objects of class ConstraintExpression +There must be p such objects (hence the parametrization on p) + +*-----------------------------------------------------------------------*/ +class KAdaptableExpression { +private: + std::vector Expr; /* Collection of constraints */ + std::string name; /* Name of this "meta-constraint" */ + double rhs; /* right-hand side of K-adaptable expression */ + + inline bool isConsistentWithDesign() const { + if(Expr.empty()) return false; + for (const auto& exp : Expr) + if (!exp.isConsistentWithDesign()) return false; + for (const auto& exp : Expr) + if (exp.sense != 'L') return false; + + return true; + } + +public: + + KAdaptableExpression(const std::vector& arg, const std::string arg_name, const double arg_rhs = 0) + : Expr(arg), name(arg_name), rhs(arg_rhs) + { + /* Convert to <= form */ + for(auto& exp:Expr) { + if (exp.sense == 'G') { + exp.sign('L'); + exp.rhs *= (-1.0); + for (auto& c : exp.varCoeffs) c *= (-1.0); + for (auto& b : exp.paramCoeffs) b *= (-1.0); + for (auto& S : exp.bilinearCoeffs) S *= (-1.0); + } + else if (exp.sense == 'E') { + std::cerr << "Cannot use equality constraints in K-adaptable expression.\n"; + exit(-2); + } + } + assert(isConsistentWithDesign()); + } + + inline bool existConstXTerms() const { + for (const auto& exp : Expr) if (exp.isEmpty()) return false; + return true; + } + inline bool existConstQTerms() const { + for (const auto& exp : Expr) if(exp.existConstQTerms()) return true; + return false; + } + inline bool existBilinearTerms() const { + for (const auto& exp : Expr) if(exp.existBilinearTerms()) return true; + return false; + } + + + + + + /*------------------------------------------ + | This function evaluates the value of: + | + | max_q min_p {cx*(p) + bq + x*(p)Sq - d(p)} + | + *------------------------------------------*/ + inline std::vector evaluate(UNCSetCPtr UncSet, + const std::vector& paramValues, + const std::vector& varValues = {}, + const bool normalized = 0) const + { + const bool qTermsInProduct = existBilinearTerms(); + const bool qTermsConst = existConstQTerms(); + const bool xTermsConst = existConstXTerms(); + + const bool donotSolveLP = (UncSet == nullptr) ? 1 : (!UncSet->isUncertain() || (UncSet->isUncertain() && !qTermsConst && !qTermsInProduct)); + + /* Deterministic case: set q equal to supplied value */ + if(donotSolveLP) { + + + if(qTermsInProduct || qTermsConst) if(paramValues.empty()) { std::cerr << " Error: Parameter values not supplied for evaluation of K-adaptable expression. \n"; exit(-1); } + if(xTermsConst) if(varValues.empty()) { std::cerr << " Error: Primal variable values not supplied for evaluation of K-adaptable expression. \n"; exit(-1); } + + + /* Value of the {...} term in max_q min_p {...} */ + std::vector ExprEval(Expr.size(), 0.0); + + /* Normalized value --> quantity dividing the {...} term in max_q min_p {...} */ + std::vector Normalization(Expr.size(), (normalized ? 0.0 : 1.0)); + + /* -d(p) */ + for(size_t p=0; p eval{*std::min_element(ExprEval.begin(), ExprEval.end())}; eval.insert(eval.end(), paramValues.begin(), paramValues.end()); + + return eval; + } + else { + + assert(UncSet->isUncertain()); + assert(qTermsInProduct || qTermsConst); + if (qTermsInProduct || xTermsConst) if (varValues.empty()) { std::cerr << " Error: Primal variable values not supplied for evaluation of K-adaptable expression. \n"; exit(-1); } + + + + int stat; + + /* Get ENV object from Uncertainty Set */ + CPXENVptr env = UncSet->getENVObject(); + + /* Clone LP object from Uncertainty Set */ + CPXLPptr lp = UncSet->getLPObject(&stat); + if (stat) {std::cerr << "Failed to clone CPLEX LP object inside evaluation of K-adaptable expression. \n"; exit(-1);} + + + + + int ind = 0; double coef = 1; char lb = 'L'; char ub = 'U'; double lbVal = -CPX_INFBOUND; double ubVal = +CPX_INFBOUND; + + /* Make it a maximization problem */ + CPXXchgobjsen(env, lp, CPX_MAX); + /* Make it an LP */ + CPXXchgprobtype(env, lp, CPXPROB_LP); + /* Replace current objective function */ + CPXXchgobj(env, lp, 1, &ind, &coef); coef = 0; for (int l = 1; l <= UncSet->getNoOfUncertainParameters(); ++l) CPXXchgobj(env, lp, 1, &l, &coef); + /* Change bounds on objective function variable to make it (-Inf, +Inf) */ + CPXXchgbds(env, lp, 1, &ind, &lb, &lbVal); + CPXXchgbds(env, lp, 1, &ind, &ub, &ubVal); + + /* Add constraints corresponding to each of the P constraints in Expr */ + for(const auto& exp: Expr) { + ConstraintExpression TauConstr("TauCon(" + exp.name + ")"); + + /* -d(p) + cx*(p) */ + double this_rhs = -exp.rhs; if (xTermsConst) for (size_t i = 0; i < exp.varCoeffs.size(); i++) this_rhs += exp.varCoeffs[i] * varValues.at(exp.varIndices.at(i)); + + TauConstr.sign('L'); + TauConstr.RHS(this_rhs); + + /* Normalized value --> quantity dividing the {...} term in max_q min_p {...} */ + double Normalization = (normalized ? 0.0 : 1.0); + + /* -bq */ + for (size_t i = 0; i < exp.paramCoeffs.size(); i++) if (exp.paramCoeffs[i] != 0.0) { + TauConstr.addTermX( exp.paramIndices[i], -exp.paramCoeffs[i]); + if (normalized) Normalization += exp.paramCoeffs[i] * exp.paramCoeffs[i]; + } + + /* -x*Sq */ + if (qTermsInProduct) for (size_t i = 0; i < exp.bilinearCoeffs.size(); i++) if (exp.bilinearCoeffs[i] != 0.0) if (varValues.at(exp.bilinearIndices[i].first) != 0.0) { + const double t = (exp.bilinearCoeffs[i] * varValues.at(exp.bilinearIndices[i].first)); + TauConstr.addTermX( exp.bilinearIndices[i].second, -t ); + if (normalized) Normalization += t * t; + } + + /* \tau */ + TauConstr.addTermX(0, std::sqrt(Normalization)); + + TauConstr.addToCplex(env, lp); + } + + + /* Solve LP */ + stat = CPXXlpopt(env, lp); + if (stat) std::cerr << "Could not solve LP inside evaluation of K-adaptable expression.\n"; + + + /* Assert that LP was solved to optimality */ + if (CPXXgetstat(env, lp) != CPX_STAT_OPTIMAL) { + std::cerr << "\n\n Could not solve LP to optimality inside the evaluation problem of K-adaptable expression.\n\n"; + exit(-2); + } + + /* Get optimal RHS value */ + double objval; CPXXgetobjval(env, lp, &objval); + + /* Get primal solution vector */ + std::vector primalX(1 + UncSet->getNoOfUncertainParameters(), objval); CPXXgetx(env, lp, &primalX[1], UncSet->getParamIndex('q', 1), UncSet->getParamIndex('q', UncSet->getNoOfUncertainParameters())); + + /* Free memory */ + if (lp != nullptr) { + stat = CPXXfreeprob(env, &lp); + if (stat) std::cerr << "CPXfreeprob failed inside K-adaptable expression, error code " << stat << "\n"; + } + + + return primalX; + } + } +}; + + + + + + + + +static inline double getViolation(const ConstraintExpression& con, + const std::vector& varValues, + const std::vector& paramValues = {}, + const bool normalized = 0) +{ + if (con.getSense() == 'E') { + return con.getViolation_fixedXQ(varValues, paramValues, normalized); + } + + // Construct K-adaptable expression with single policy + std::vector exp{con}; + KAdaptableExpression CExpr(exp, con.getName()); + + // evaluate worst-case value + const auto q = CExpr.evaluate(nullptr, paramValues, varValues, normalized); + assert(!q.empty()); + + return q[0]; +} + +static inline std::vector getViolation(const ConstraintExpression& con, + UNCSetCPtr UncSet, + const std::vector& varValues = {}, + const bool normalized = 0) +{ + // temporary var + std::vector paramValues; + + if (UncSet == nullptr && con.getSense() == 'E') { + const double viol = con.getViolation_fixedXQ(varValues, paramValues, normalized); + std::vector q{viol}; + q.insert(q.end(), paramValues.begin(), paramValues.end()); + return q; + } + + // Construct K-adaptable expression with single policy + std::vector exp{con}; + KAdaptableExpression CExpr(exp, con.getName()); + + // evaluate worst-case value + const auto q = CExpr.evaluate(UncSet, paramValues, varValues, normalized); + assert((int)q.size() == 1 + UncSet->getNoOfUncertainParameters()); + + return q; +} + +#endif \ No newline at end of file diff --git a/inc/indexingTools.h b/inc/indexingTools.h new file mode 100644 index 0000000..082eafe --- /dev/null +++ b/inc/indexingTools.h @@ -0,0 +1,182 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Nikolaos Lappas and Chrysanthos Gounaris */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ +#ifndef IND_TOOLS +#define IND_TOOLS + +// Variable Indexing and Querying Class +// supporting up to 4 indices + +#include +#include +#include +#include +#include + +// v2: re-arranged order of indices in VarType +// v3: includes support for a fifth index +class VarInfo { + +public: + + // basic internal data structure. Map with key the variable name, and a tuple + // with fundamental sizes + // "VarName", < First, Last, numDimensions, size1, size2, size3, size4, size5> + // get 0 1 2 3 4 5 6 7 + std::map > VarType; + int totalVars; // stores the number of elements of these vars + std::vector UpperBound; // stores the upper bounds + std::vector LowerBound; // stores the lower bounds + std::vector ObjCoefficient; // stores the objective coefficients + std::vector ColumnType; // stores the type of variable: continuous 'C' or binary 'B' + std::vector UndefinedVar; // stores if a specific variable is undefined or not + std::vector UndefinedVarCount; // stores the number of variables undefined up to a specific index + + bool isVarTypeConsistent(const int & ind1, const int & ind2, const int & ind3, const int & ind4, const int & ind5); + + void linTo1dIndex(const int & linIndex, const std::tuple & VarType, int & ind1) const; + + void linTo2dIndex(const int & linIndex, const std::tuple & VarType, int & ind1, int & ind2) const; + + void linTo3dIndex(const int & linIndex, const std::tuple & VarType, int & ind1, int & ind2, int & ind3) const; + + void linTo4dIndex(const int & linIndex, const std::tuple & VarType, int & ind1, int & ind2, int & ind3, int & ind4) const; + + void linTo5dIndex(const int & linIndex, const std::tuple & VarType, int & ind1, int & ind2, int & ind3, int & ind4, int & ind5) const; + + void linToIndex(const int & linIndex, const std::tuple & VarType, int & ind1, int & ind2, int & ind3, int & ind4, int & ind5) const; + + int xdIndicesToLinIndex(const std::tuple & VarType, const int & ind1, const int & ind2, const int & ind3, const int & ind4, const int & ind5) const; + +public: + + // Clear the object + void clear(); + + // Add a new variable type + void addVarType(const std::string & name, const char & type, const double & lb, const double & ub, const int & ind1, const int & ind2 = -1, const int & ind3 = -1, const int & ind4 = -1, const int & ind5 = -1); + + // set the given var to be undefined + void setUndefinedVar(const int & index); + + // set the UB of the given var + void setVarUB(const double & Ub, const int & index); + + // set the LB of the given var + void setVarLB(const double & Lb, const int & index); + + // set the column type of the given var + void setVarColType(const char & Ctype, const int & index); + + // set the objective coefficient of the given var + void setVarObjCoeff(const double & Ocoeff, const int & index); + + // set the given var to be undefined + void setUndefinedVar(const std::string & type, const int ind1 = -1, const int ind2 = -1, const int ind3 = -1, const int ind4 = -1, const int ind5 = -1); + + // set the UB of the given var + void setVarUB(const double & Ub, const std::string & type, const int ind1 = -1, const int ind2 = -1, const int ind3 = -1, const int ind4 = -1, const int ind5 = -1); + + // set the LB of the given var + void setVarLB(const double & Lb, const std::string & type, const int ind1 = -1, const int ind2 = -1, const int ind3 = -1, const int ind4 = -1, const int ind5 = -1); + + // set the column type of the given var + void setVarColType(const char & Ctype, const std::string & type, const int ind1 = -1, const int ind2 = -1, const int ind3 = -1, const int ind4 = -1, const int ind5 = -1); + + // set the objective coefficient of the given var + void setVarObjCoeff(const double & Ocoeff, const std::string & type, const int ind1 = -1, const int ind2 = -1, const int ind3 = -1, const int ind4 = -1, const int ind5 = -1); + + // set the vector of undefined vars - must be of size getTotalVarSize() + // void setUndefinedVars(const std::vector & Undefined); + + // // set the vector with the UB of the vars - must be of size getTotalVarSize() + // void setVarsUB(const std::vector & Ubounds); + + // // set the vector with the LB of the vars - must be of size getTotalVarSize() + // void setVarsLB(const std::vector & Lbounds); + + // // set the vector with the column types - must be of size getTotalVarSize() + // void setVarsColType(const std::vector & Ctype); + + // // set the vector with the objective coefficients - must be of size getTotalVarSize() + // void setVarsObjCoeff(const std::vector & Ocoeff); + + + + + + + // get the number of vars that this object holds + int getTotalVarSize() const; + + // get the name and the indices, given the linear index + void getVarInfo(const int & linIndex, std::string &type, int & ind1, int & ind2, int & ind3, int & ind4, int & ind5) const; + + // get the linear index, given the name and the indices + int getVarLinIndex(const std::string & type, const int ind1 = -1, const int ind2 = -1, const int ind3 = -1, const int ind4 = -1, const int ind5 = -1) const; + + // get the number of variables of the requested type + int getVarTypeSize(const std::string & name) const; + + // get the linear index of the first var of a given var type (iterating tool) + int getFirstOfVarType(const std::string & name) const; + + // get the linear index of the last var of a given var type (iterating tool) + int getLastOfVarType(const std::string & name) const; + + // get the variable name, given the linear index as returned by getVarLinIndex() + std::string getVarName(const int & index) const; + + // query if a var is undefined, given the linear index as returned by getVarLinIndex() + bool isUndefVar(const int & index) const; + + // get the UB of a specific var, given the linear index as returned by getVarLinIndex() + double getVarUB(const int & index) const; + + // get the LB of a specific var, given the linear index as returned by getVarLinIndex() + double getVarLB(const int & index) const; + + // get the type of a specific var, given the linear index as returned by getVarLinIndex() + char getVarColType(const int & index) const; + + // get the objective coefficient of a specific var, given the linear index as returned by getVarLinIndex() + double getVarObjCoeff(const int & index) const; + + + + + + + // get the true number of defined vars that this object holds + int getTotalDefVarSize() const; + + // get the true linear index of defined var, given the original index as returned by getVarLinIndex() + int getDefVarLinIndex(const int & index) const; + + // get the true linear index of defined var, given the name and the indices + int getDefVarLinIndex(const std::string & type, const int ind1 = -1, const int ind2 = -1, const int ind3 = -1, const int ind4 = -1, const int ind5 = -1) const; + + // get the number of defined vars of the requested type + int getDefVarTypeSize(const std::string & name) const; + + // get the true linear index of the first defined var of a given var type (iterating tool) + int getFirstDefOfVarType(const std::string & name) const; + + // get the true linear index of the last defined var of a given var type (iterating tool) + int getLastDefOfVarType(const std::string & name) const; + +}; + + +#endif + + + diff --git a/inc/instance_knp.hpp b/inc/instance_knp.hpp new file mode 100644 index 0000000..9b7ff0f --- /dev/null +++ b/inc/instance_knp.hpp @@ -0,0 +1,166 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#ifndef KNP_INSTANCE_HPP +#define KNP_INSTANCE_HPP + + +#include +#include +#include +#include +#include //std::sort, std::random_shuffle +#include +#include + +#define KNP_NFACTORS 4 + +/** + * Class represesnting an instance of the Knapsack Problem + * (as defined in the K-Adaptability paper) + * + * Note: All arrays and matrices are 1-indexed + */ +struct KNP { + /** # of items */ + int N; + + /** Cost of items */ + std::vector cost; + + /** Budget (size of knapsack) */ + double B; + + /** Profit of items */ + std::vector profit; + + /** Late-stage profit penalty */ + double theta; + + /** Option of taking loans */ + bool loan; + + /** Cost of borrowing in the first-stage */ + double ell; + + /** Cost of borrowing in the second-stage */ + double lambda; + + /** Risk factor coefficients for costs */ + std::vector > phi; + + /** Risk factor coefficients for profits */ + std::vector > ksi; + + /** Solution file name */ + std::string solfilename; + +}; + + + + +/** + * Construct an instance of the Knapsack Problem + * @param data (reference to) instance of Knapsack Problem + * @param n number of items + * @param seed seed that is unique to this instance + * @param mixed_integer true if mixed-integer knapsack problem, otherwise pure integer + */ +static inline void gen_KNP(KNP& data, unsigned int n, int seed = 1, bool mixed_integer = false) { + + if (n < 5) { + fprintf(stderr, "warning: N < 5 in Knapsack Problem. Setting N = 5.\n"); + n = 5; + } + + // seed for re-producibility + std::default_random_engine gen (1111 + seed); + std::uniform_real_distribution interval_1 (0.0, 1.0); + std::uniform_real_distribution interval_10 (0.0, 10.0); + + // # of items + data.N = n; + data.B = 0.0; + data.theta = 0.8; + data.loan = mixed_integer; + data.ell = 0.12; + data.lambda = 1.2; + + // Cost of items + data.cost.assign(1 + data.N, 0.0); + for (n = 1; (int)n <= data.N; n++) { + data.cost[n] = interval_10(gen); + data.B += data.cost[n]; + } + + // Budget (size of knapsack) + data.B /= 2.0; + + + // Profit of items + data.profit.assign(1 + data.N, 0.0); + for (n = 1; (int)n <= data.N; n++) { + data.profit[n] = data.cost[n] / 5.0; + } + + + // Risk factor coefficients + data.phi.resize(1 + data.N); + data.ksi.resize(1 + data.N); + for (n = 1; (int)n <= data.N; n++) { + // generate F-1 numbers between 0 and 1 + std::array xn, yn; + for (size_t f = 0; f < xn.size(); f++) { + xn[f] = interval_1(gen); + yn[f] = interval_1(gen); + } + + // sort these numbers + std::sort(xn.begin(), xn.end()); + std::sort(yn.begin(), yn.end()); + + // generate F numbers such that their sum = 1 + data.phi[n][0] = 0.0; + data.phi[n][1] = xn[0]; + data.ksi[n][0] = 0.0; + data.ksi[n][1] = yn[0]; + for (size_t f = 2; f < KNP_NFACTORS; f++) { + data.phi[n][f] = xn[f-1] - xn[f-2]; + data.ksi[n][f] = yn[f-1] - yn[f-2]; + } + data.phi[n][KNP_NFACTORS] = 1 - xn[KNP_NFACTORS-2]; + data.ksi[n][KNP_NFACTORS] = 1 - yn[KNP_NFACTORS-2]; + + // randomize in order to eliminate bias + std::random_shuffle(data.phi[n].begin() + 1, data.phi[n].end()); + std::random_shuffle(data.ksi[n].begin() + 1, data.ksi[n].end()); + + // check (just to be safe) + #ifndef NDEBUG + double check_phi = 0.0; + double check_ksi = 0.0; + for (size_t f = 1; f <= KNP_NFACTORS; f++) { + check_phi += data.phi[n][f]; + check_ksi += data.ksi[n][f]; + } + assert(std::abs(1 - check_phi) <= 1.E-12); + assert(std::abs(1 - check_ksi) <= 1.E-12); + #endif + } + + // set solution file name + data.solfilename = "KNP-n" + std::to_string(data.N) + "-s" + std::to_string(seed) + "-t" + std::to_string(data.loan) + ".opt"; +} +#undef KNP_NFACTORS +#endif \ No newline at end of file diff --git a/inc/instance_psp.hpp b/inc/instance_psp.hpp new file mode 100644 index 0000000..b730e88 --- /dev/null +++ b/inc/instance_psp.hpp @@ -0,0 +1,109 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#ifndef PSP_INSTANCE_HPP +#define PSP_INSTANCE_HPP + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include //std::sort, std::random_shuffle +#include + +#define PSP_NFACTORS 4 + +/** + * Class represesnting an instance of the Project Scheduling Problem + * (as defined in the K-Adaptability paper) + * + * Note: All arrays and matrices are 1-indexed + */ +struct PSP { + /** # of nodes in project graph = # of activities including dummy start and end nodes */ + int N; + + /** Budget of resources */ + double B; + + /** Use affine decision rules */ + bool affine; + + /** Special case of project graphs */ + bool special; + + /** Activity durations */ + std::vector duration; + + /** All nodes in the graph which succeed this node */ + std::vector > successor; + + /** Risk factor coefficients */ + std::vector > phi; + + /** Solution file name */ + std::string solfilename; + +}; + + + +/** + * Construct an instance of the Project Scheduling Problem + * See Example 2.2 in "Robust Resource Allocations in Temporal Networks", Math. Prog. (2012) + * @param data (reference to) instance of Project Scheduling Problem + * @param k parameter to generate instances of Example 2.2 + * @param seed seed that is unique to this instance + */ +static inline void gen_PSP(PSP& data, unsigned int k, int seed = 1) { + + if (k == 0) { + fprintf(stderr, "warning: k = 0 in Resource Allocation Problem (special case). Setting k = 1.\n"); + k = 1; + } + + // initialize + data.N = (3 * k) + 1; + data.B = 0; + data.affine = (seed != 0); + data.special = true; + data.duration.resize(1 + data.N, 0.0); + data.successor.resize(1 + data.N); + data.phi.resize(data.N); + data.solfilename = "PSP-n" + std::to_string(data.N) + "-s" + std::to_string(seed) + "-d0.opt"; + + // generate project graph + for (unsigned int l = 0; l < k; l++) for (unsigned int p = 2; p <= 3; p++) { + data.successor[(3 * l) + 1].emplace_back((3 * l) + p); + data.successor[(3 * l) + p].emplace_back((3 * l) + 4); + } + + // set durations + for (unsigned int l = 0; l < k; l++) { + data.duration[(3 * l) + 1] = 0; + data.duration[(3 * l) + 2] = 1; + data.duration[(3 * l) + 3] = -1; + } + data.duration[data.N] = 0; + + // dummy end node has no successors + assert(data.successor[data.N].empty()); +} +#undef PSP_NFACTORS +#endif \ No newline at end of file diff --git a/inc/instance_spp.hpp b/inc/instance_spp.hpp new file mode 100644 index 0000000..174e7db --- /dev/null +++ b/inc/instance_spp.hpp @@ -0,0 +1,134 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#ifndef SPP_INSTANCE_HPP +#define SPP_INSTANCE_HPP + + +#include +#include +#include +#include +#include + +/** + * Class represesnting an instance of the Shortest Path Problem + * (as defined in the K-Adaptability paper) + * + * Note: All arrays and matrices are 1-indexed + */ +struct SPP { + /** # of nodes in graph */ + int N; + + /** # of arcs in graph */ + int A; + + /** Arc-weights (an entry of -1 indicates the arc does not exist) */ + std::vector > costs; + + /** start node */ + int src; + + /** end node */ + int tgt; + + /** Solution file name */ + std::string solfilename; + +}; + + + + +/** + * Construct an instance of the Shortest Path Problem + * @param data (reference to) instance of Shortest Path Problem + * @param n number of nodes + * @param seed seed that is unique to this instance + */ +static inline void gen_SPP(SPP& data, unsigned int n, int seed = 1) { + + if (n < 20) { + fprintf(stderr, "warning: N < 20 in Shortest Path Problem. Setting N = 20.\n"); + n = 20; + } + + // seed for re-producibility + std::default_random_engine gen (1111 + seed); + + // # of nodes, # of arcs + data.N = n; + data.A = n * (n - 1); + + // costs(i,j) = weight of arc from i to j for i, j = 1, ..., N + data.costs.assign(1 + data.N, std::vector(1 + data.N, -1.0)); + + // coordinates(i) = co-ordinates of the i'th node + std::vector > coordinates (1); + std::uniform_real_distribution interval (0.0, 10.0); + for (n = 1; (int)n <= data.N; ++n) { + double xn = interval(gen); + double yn = interval(gen); + coordinates.emplace_back(std::vector{xn, yn}); + } + + // compute arc-weights + for (int i = 1; i <= data.N; ++i) for (int j = i + 1; j <= data.N; ++j) { + double xd = std::abs(coordinates[i][0] - coordinates[j][0]); + double yd = std::abs(coordinates[i][1] - coordinates[j][1]); + data.costs[i][j] = std::sqrt((xd * xd) + (yd * yd)); + data.costs[j][i] = data.costs[i][j]; + } + + // delete 70\% of the arcs + int D = (7 * data.A) / 10; + for (int count = 1; count <= D; ++count) { + + // Compute max element in cost matrix + double max = -1.0; + int s = 0; + int t = 0; + for (int i = 1; i <= data.N; ++i) { + data.costs[i][i] = -1.0; + for (int j = 1; j <= data.N; ++j) if (i != j) { + if (data.costs[i][j] > max) { + max = data.costs[i][j]; + s = i; + t = j; + } + } + } + assert(max >= 0.0); + assert(s > 0 && t > 0); + assert(s != t); + + // Kill this arc by setting cost = -1 + data.costs[s][t] = -1.0; + data.A--; + + // Source and terminal nodes + if (count == 1) { + data.src = s; + data.tgt = t; + } + } + assert(data.src >= 1 && data.src <= data.N); + assert(data.tgt >= 1 && data.tgt <= data.N); + assert(data.src != data.tgt); + + // set solution file name + data.solfilename = "SPP-n" + std::to_string(data.N) + "-s" + std::to_string(seed) + ".opt"; +} + +#endif \ No newline at end of file diff --git a/inc/problemInfo.hpp b/inc/problemInfo.hpp new file mode 100644 index 0000000..0c3473d --- /dev/null +++ b/inc/problemInfo.hpp @@ -0,0 +1,317 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#ifndef KADAPTABLEINFO_HPP +#define KADAPTABLEINFO_HPP + +#include "indexingTools.h" +#include "constraintExpr.hpp" +#include "uncertainty.hpp" +#include +#include + +/** + * Abstract class representing information that any + * two-/multi-stage RO instance must provide. + * All problems must be minimization problems. + */ +class KAdaptableInfo { + +protected: + + /** Continuous or (mixed) integer problem? */ + bool hasInteger; + + /** Uncertainty in objective only */ + bool objectiveUnc; + + /** Any 1st-stage decisions to be made? */ + bool existsFirstStage; + + /** # of 1st-stage variables including worst-case objective */ + unsigned int numFirstStage; + + /** # of 2nd-stage variables */ + unsigned int numSecondStage; + + /** # of 2nd-stage policies */ + unsigned int numPolicies; + + /** Solution file name */ + std::string solfilename; + + /** Resident uncertainty set */ + UncertaintySet U; + + /** 1st-stage variables only */ + VarInfo X; + + /** 2nd-stage variables only */ + VarInfo Y; + + /** Constraints with 1st-stage variables only (excluding bounds) */ + std::vector C_X; + + /** Constraints with 2nd-(and, if applicable, 1st-) stage variables BUT NO uncertain parameters (excluding bounds) */ + std::vector > C_XY; + + /** Constraints with 1st-stage variables AND uncertain parameters */ + std::vector C_XQ; + + /** Constraints with 2nd-(and, if applicable, 1st-) stage variables AND uncertain parameters */ + std::vector > C_XYQ; + + /** Bound constraints on 1st-stage variables only */ + std::vector B_X; + + /** Bound constraints on 2nd-stage variables only */ + std::vector > B_Y; + + /** + * Define 1st-stage and 2nd-stage variables + */ + virtual void makeVars() = 0; + + /** + * Define uncertainty set + */ + virtual void makeUncSet() = 0; + + /** + * Define constraints involving 1st-stage variables only + */ + virtual void makeConsX() = 0; + + /** + * Define constraints involving 2nd-stage variables + * @param k policy number to make constraints for + */ + virtual void makeConsY(unsigned int k = 0) = 0; + +public: + /** + * Virtual destructor does nothing + */ + virtual ~KAdaptableInfo() {/**/} + + /** + * Virtual (default) constructor + */ + virtual KAdaptableInfo* create() const = 0; + + /** + * Virtual (copy) constructor + */ + virtual KAdaptableInfo* clone() const = 0; + + /** + * Resizes data structures to contain at least K 2nd-stage policies + * @param K Number of 2nd-stage policies to be supported + */ + void resize(unsigned int K); + + /** + * Check consistency with class design + * @return true if object is conistent with class design + */ + bool isConsistentWithDesign() const; + + /** + * Linear variable index of 1st-stage variable + * [All indices start from 0] + * + * @param type name of variable + * @param ind1 first index + * @param ind2 second index + * @param ind3 third index + * @param ind4 fourth index + * @return the linear index of the variable + */ + inline int getVarIndex_1(const std::string type, const int ind1, const int ind2 = -1, const int ind3 = -1, const int ind4 = -1) const { + return X.getDefVarLinIndex(type, ind1, ind2, ind3, ind4); + } + + /** + * Linear variable index of 2nd-stage variable + * [All indices start from 0] + * [Policy numbering starts from 0] + * + * @param k policy number in which you want linear index + * @param type name of variable + * @param ind1 first index + * @param ind2 second index + * @param ind3 third index + * @param ind4 fourth index + * @return the linear index of the variable + */ + inline int getVarIndex_2(const unsigned int k, const std::string type, const int ind1, const int ind2 = -1, const int ind3 = -1, const int ind4 = -1) const { + return numFirstStage + (k * numSecondStage) + Y.getDefVarLinIndex(type, ind1, ind2, ind3, ind4); + } + + /** + * Get total # of 1st- and 2nd-stage decisions + * @param K Total number of 2nd-stage policies to consider + * @return Total # of 1st- and 2nd-stage decisions for a K-policy problem + */ + inline int getNumVars(const unsigned int K = 1) const { + return numFirstStage + (K * numSecondStage); + } + + /** + * Get # of 1st-stage decisions + * @return # of 1st-stage decisions + */ + inline int getNumFirstStage() const { + return numFirstStage; + } + + /** + * Get # of 2nd-stage decisions + * @return # of 2nd-stage decisions + */ + inline int getNumSecondStage() const { + return numSecondStage; + } + + /** + * Get current upper bound on # of 2nd-stage policies + * @return current upper bound on # of 2nd-stage policies + */ + inline int getNumPolicies() const { + return numPolicies; + } + + /** + * Get solution file name + * @return solution file name + */ + inline std::string getSolFileName() const { + return solfilename; + } + + /** + * Get nominal parameter values + * @return nominal parameter values + */ + inline std::vector getNominal() const { + return U.getNominal(); + } + + /** + * Get # of uncertain parameters + * @return # of uncertain parameters + */ + inline int getNoOfUncertainParameters() const { + return U.getNoOfUncertainParameters(); + } + + /** + * Indicates if problem has continuous variables only + * @return true if problem has continuous variables only + */ + inline bool isContinuous() const { + return !hasInteger; + } + + /** + * Indicates if problem has objective uncertainty only + * @return true if problem has objective uncertainty only + */ + inline bool hasObjectiveUncOnly() const { + return objectiveUnc; + } + + /** + * Indicates if problem has 2nd-stage decisions only + * @return true if problem has 2nd-stage decisions only + */ + inline bool isSecondStageOnly() const { + return !existsFirstStage; + } + + /** + * Get resident uncertainty set + * @return return the uncertainty set + */ + inline const decltype(U)& getUncSet() { + return U; + } + + /** + * Get 1st-stage variables only + * @return the 1st-stage variables + */ + inline const decltype(X)& getVarsX() { + return X; + } + + /** + * Get 2nd-stage variables only + * @return the 2nd-stage variables + */ + inline const decltype(Y)& getVarsY() { + return Y; + } + + /** + * Get constraints with 1st-stage variables only (excluding bounds) + * @return constraints with 1st-stage variables only (excluding bounds) + */ + inline const decltype(C_X)& getConstraintsX() { + return C_X; + } + + /** + * Get constraints with 2nd-(and, if applicable, 1st-) stage variables BUT NO uncertain parameters (excluding bounds) + * @return constraints with 2nd-(and, if applicable, 1st-) stage variables BUT NO uncertain parameters (excluding bounds) + */ + inline const decltype(C_XY)& getConstraintsXY() { + return C_XY; + } + + /** + * Get constraints with 1st-stage variables AND uncertain parameters + * @return constraints with 1st-stage variables AND uncertain parameters + */ + inline const decltype(C_XQ)& getConstraintsXQ() { + return C_XQ; + } + + /** + * Get constraints with 2nd-(and, if applicable, 1st-) stage variables AND uncertain parameters + * @return constraints with 2nd-(and, if applicable, 1st-) stage variables AND uncertain parameters + */ + inline const decltype(C_XYQ)& getConstraintsXYQ() { + return C_XYQ; + } + + /** + * Get bound constraints on 1st-stage variables only + * @return bound constraints on 1st-stage variables only + */ + inline const decltype(B_X)& getBoundsX() { + return B_X; + } + + /** + * Get bound constraints on 2nd-stage variables only + * @return bound constraints on 2nd-stage variables only + */ + inline const decltype(B_Y)& getBoundsY() { + return B_Y; + } + +}; + + +#endif \ No newline at end of file diff --git a/inc/problemInfo_knp.hpp b/inc/problemInfo_knp.hpp new file mode 100644 index 0000000..e6c5ca3 --- /dev/null +++ b/inc/problemInfo_knp.hpp @@ -0,0 +1,72 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#ifndef KADAPTABLEINFO_KNP_HPP +#define KADAPTABLEINFO_KNP_HPP + +#include "instance_knp.hpp" +#include "problemInfo.hpp" + +class KAdaptableInfo_KNP : public KAdaptableInfo { +protected: + /** Specific instance data */ + KNP data; + + /** + * Define 1st-stage and 2nd-stage variables + */ + void makeVars() override; + + /** + * Define uncertainty set + */ + void makeUncSet() override; + + /** + * Define constraints involving 1st-stage variables only + */ + void makeConsX() override; + + /** + * Define constraints involving 2nd-stage variables + * @param k policy number to make constraints for + */ + void makeConsY(unsigned int k = 0) override; + +public: + + /** + * Set data of the instance that this object refers to + * @param data problem instance + */ + void setInstance(const KNP& data); + + /** + * Virtual (default) constructor + */ + inline KAdaptableInfo_KNP* create() const override { + return new KAdaptableInfo_KNP(); + } + + /** + * Virtual (copy) constructor + */ + inline KAdaptableInfo_KNP* clone() const override { + return new KAdaptableInfo_KNP (*this); + } + + + +}; + +#endif \ No newline at end of file diff --git a/inc/problemInfo_psp.hpp b/inc/problemInfo_psp.hpp new file mode 100644 index 0000000..afb4504 --- /dev/null +++ b/inc/problemInfo_psp.hpp @@ -0,0 +1,72 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#ifndef KADAPTABLEINFO_PSP_HPP +#define KADAPTABLEINFO_PSP_HPP + +#include "instance_psp.hpp" +#include "problemInfo.hpp" + +class KAdaptableInfo_PSP : public KAdaptableInfo { +protected: + /** Specific instance data */ + PSP data; + + /** + * Define 1st-stage and 2nd-stage variables + */ + void makeVars() override; + + /** + * Define uncertainty set + */ + void makeUncSet() override; + + /** + * Define constraints involving 1st-stage variables only + */ + void makeConsX() override; + + /** + * Define constraints involving 2nd-stage variables + * @param k policy number to make constraints for + */ + void makeConsY(unsigned int k = 0) override; + +public: + + /** + * Set data of the instance that this object refers to + * @param data problem instance + */ + void setInstance(const PSP& data); + + /** + * Virtual (default) constructor + */ + inline KAdaptableInfo_PSP* create() const override { + return new KAdaptableInfo_PSP(); + } + + /** + * Virtual (copy) constructor + */ + inline KAdaptableInfo_PSP* clone() const override { + return new KAdaptableInfo_PSP (*this); + } + + + +}; + +#endif \ No newline at end of file diff --git a/inc/problemInfo_spp.hpp b/inc/problemInfo_spp.hpp new file mode 100644 index 0000000..853fc6f --- /dev/null +++ b/inc/problemInfo_spp.hpp @@ -0,0 +1,72 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#ifndef KADAPTABLEINFO_SPP_HPP +#define KADAPTABLEINFO_SPP_HPP + +#include "instance_spp.hpp" +#include "problemInfo.hpp" + +class KAdaptableInfo_SPP : public KAdaptableInfo { +protected: + /** Specific instance data */ + SPP data; + + /** + * Define 1st-stage and 2nd-stage variables + */ + void makeVars() override; + + /** + * Define uncertainty set + */ + void makeUncSet() override; + + /** + * Define constraints involving 1st-stage variables only + */ + void makeConsX() override; + + /** + * Define constraints involving 2nd-stage variables + * @param k policy number to make constraints for + */ + void makeConsY(unsigned int k = 0) override; + +public: + + /** + * Set data of the instance that this object refers to + * @param data problem instance + */ + void setInstance(const SPP& data); + + /** + * Virtual (default) constructor + */ + inline KAdaptableInfo_SPP* create() const override { + return new KAdaptableInfo_SPP(); + } + + /** + * Virtual (copy) constructor + */ + inline KAdaptableInfo_SPP* clone() const override { + return new KAdaptableInfo_SPP (*this); + } + + + +}; + +#endif \ No newline at end of file diff --git a/inc/robustSolver.hpp b/inc/robustSolver.hpp new file mode 100644 index 0000000..68694cc --- /dev/null +++ b/inc/robustSolver.hpp @@ -0,0 +1,502 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#ifndef KADAPTABLESOLVER_HPP +#define KADAPTABLESOLVER_HPP + +#include "problemInfo.hpp" +#include +#include + + + +/** + * Class to solve all CPLEX-based computational examples + * in K-adaptability paper. + * + * Note that the problem is stored in pInfo (temporary variable), + * which must be in the same scope as the solver object. + */ +class KAdaptableSolver { +protected: + + /** Pointer to const information */ + KAdaptableInfo *pInfo = NULL; + + /** Best known solution for K-adaptable problem (in MILP representation) */ + std::vector xsol; + + + + + + +public: + /** + * Delete default constructor + */ + KAdaptableSolver() = delete; + + + /** + * Constructor takes a specific instance as input + */ + KAdaptableSolver(const KAdaptableInfo& pInfoData); + + + /** + * Define pre-C++0x methods + */ + ~KAdaptableSolver(); + KAdaptableSolver (const KAdaptableSolver&); + KAdaptableSolver& operator=(const KAdaptableSolver&); + + + /** + * Delete C++11 methods + */ + KAdaptableSolver (KAdaptableSolver&&) = delete; + KAdaptableSolver& operator=(KAdaptableSolver&&) = delete; + + + /** + * Reset solver with information about a new problem + * @param pInfoData reference to a new problem instance + */ + void setInfo(const KAdaptableInfo& pInfoData); + + + /** + * Reset the solver with a new problem + * @param pInfoData reference to a new problem instance + * @param K upper bound on # of 2nd-stage policies to be supported + */ + void reset(const KAdaptableInfo& pInfoData, const unsigned int K = 1); + + + + + + +public: + + /** + * Set the best known solution for K-adaptable problem + * @param x candidate best known solution for K-adaptable problem + * @param K # of second-stage policies + */ + void setX(const std::vector& x, const unsigned int K); + + /** + * Return # of 2nd-stage policies in candidate solution + * @param x candidate solution + * @return # of 2nd-stage policies in x + */ + unsigned int getNumPolicies(const std::vector& x) const; + + /** + * Resize the K'-adaptable solution x to the size of a K-adaptable solution + * @param x candidate solution + * @param K # of 2nd-stage policies + */ + void resizeX(std::vector& x, const unsigned int K) const; + + /** + * Remove policy number k of K-adaptable solution x (remove only 2nd-stage variables) + * @param x candidate solution + * @param K # of 2nd-stage policies + * @param k policy number of interest + */ + void removeXPolicy(std::vector& x, const unsigned int K, const unsigned int k) const; + + /** + * Return policy number k of K-adaptable solution x (1st-stage and 2nd-stage variables) + * @param x candidate solution + * @param K # of 2nd-stage policies + * @param k policy number of interest + * @return policy k of candidate solution + */ + const std::vector getXPolicy(const std::vector& x, const unsigned int K, const unsigned int k) const; + + + + + + +public: + /** + * Add 1st-stage variables and constraints (not involving uncertain parameters) + * [Does not check if model is of the correct size] + * + * @param env solver environment object + * @param lp solver model object + */ + void updateX(CPXCENVptr env, CPXLPptr lp) const; + + /** + * Add 1st-stage constraints involving uncertain parameters + * [Does not check if model is of the correct size] + * + * @param env solver environment object + * @param lp solver model object + * @param q candidate scenario (if empty, then attempt to reformulate via dualization) + * @param lazy add constraints in a lazy manner + */ + void updateXQ(CPXCENVptr env, CPXLPptr lp, const std::vector& q = {}, const bool lazy = false) const; + + /** + * Add 2nd-stage variables and constraints (not involving uncertain parameters) for policy k + * [Does not check if model is of the correct size] + * + * @param env solver environment object + * @param lp solver model object + * @param k policy to which to add 2nd-stage variables and constraints + */ + void updateY(CPXCENVptr env, CPXLPptr lp, const unsigned int k) const; + + /** + * Add 2nd-stage constraints involving uncertain parameters in policy k + * [Does not check if model is of the correct size] + * + * @param env solver environment object + * @param lp solver model object + * @param k policy to which to add constraints + * @param q candidate scenario (if empty, then attempt to reformulate via dualization) + * @param lazy add constraints in a lazy manner + */ + void updateYQ(CPXCENVptr env, CPXLPptr lp, const unsigned int k, const std::vector& q = {}, const bool lazy = false) const; + + /** + * Get all constraints involving uncertain parameters and first-stage variables using realization q + * @param q candidate scenario + * @param rcnt # of constraints will be returned here (size of rhs, sense, rmatbeg) + * @param nzcnt # of non-zeros across all constraints (size of rmatind, rmatval) + * @param rhs rhs of each constraint + * @param sense sense of each constraint + * @param rmatbeg rmatbeg[i] = starting position of constraint i in rmatind/rmatval, rmatbeg[0] = 0 + * @param rmatind linear indices of variables appearing in constraints + * @param rmatval coefficients of variables appearing in constraints + */ + void getXQ_fixedQ(const std::vector& q, CPXDIM& rcnt, CPXNNZ& nzcnt, std::vector& rhs, std::vector& sense, std::vector& rmatbeg, std::vector& rmatind, std::vector& rmatval) const; + + /** + * Get all constraints involving uncertain parameters and first-stage variables using decisions x + * @param x candidate solution + * @param rcnt # of constraints will be returned here (size of rhs, sense, rmatbeg) + * @param nzcnt # of non-zeros across all constraints (size of rmatind, rmatval) + * @param rhs rhs of each constraint + * @param sense sense of each constraint + * @param rmatbeg rmatbeg[i] = starting position of constraint i in rmatind/rmatval, rmatbeg[0] = 0 + * @param rmatind linear indices of parameters appearing in constraints + * @param rmatval coefficients of parameters appearing in constraints + */ + void getXQ_fixedX(const std::vector& x, CPXDIM& rcnt, CPXNNZ& nzcnt, std::vector& rhs, std::vector& sense, std::vector& rmatbeg, std::vector& rmatind, std::vector& rmatval) const; + + /** + * Get all constraints involving uncertain parameters in policy k using realization q + * @param k policy for which to get constraints + * @param q candidate scenario + * @param rcnt # of constraints will be returned here (size of rhs, sense, rmatbeg) + * @param nzcnt # of non-zeros across all constraints (size of rmatind, rmatval) + * @param rhs rhs of each constraint + * @param sense sense of each constraint + * @param rmatbeg rmatbeg[i] = starting position of constraint i in rmatind/rmatval, rmatbeg[0] = 0 + * @param rmatind linear indices of variables appearing in constraints + * @param rmatval coefficients of variables appearing in constraints + */ + void getYQ_fixedQ(const unsigned int k, const std::vector& q, CPXDIM& rcnt, CPXNNZ& nzcnt, std::vector& rhs, std::vector& sense, std::vector& rmatbeg, std::vector& rmatind, std::vector& rmatval) const; + + /** + * Get all constraints involving uncertain parameters in policy k using decisions x + * @param k policy for which to get constraints + * @param x candidate solution + * @param rcnt # of constraints will be returned here (size of rhs, sense, rmatbeg) + * @param nzcnt # of non-zeros across all constraints (size of rmatind, rmatval) + * @param rhs rhs of each constraint + * @param sense sense of each constraint + * @param rmatbeg rmatbeg[i] = starting position of constraint i in rmatind/rmatval, rmatbeg[0] = 0 + * @param rmatind linear indices of parameters appearing in constraints + * @param rmatval coefficients of parameters appearing in constraints + */ + void getYQ_fixedX(const unsigned int k, const std::vector& x, CPXDIM& rcnt, CPXNNZ& nzcnt, std::vector& rhs, std::vector& sense, std::vector& rmatbeg, std::vector& rmatind, std::vector& rmatval) const; + + + + + + +public: + + /** + * Get a lower bound on the fully adaptive solution value + * @return lower bound on fully adaptive solution value + */ + double getLowerBound() const; + + /** + * Get the worst-case deterministic objective value + * @param q scenario which generates the worst-case deterministic value is returned here + * @return worst-case deterministic objective value + */ + double getLowerBound_2(std::vector& q) const; + + /** + * Get the worst-case objective value of the given solution + * @param x candidate solution + * @param K # of 2nd-stage policies + * @param q scenario which generates the worst-case objective value is returned here + * @return worst-case objective value + */ + double getWorstCase(const std::vector& x, const unsigned int K, std::vector& q) const; + + + + + + +public: + + /** + * Check if the given solution is 'structurally' feasible for the K-adaptable problem + * Structural feasibility = satisfaction of all uncertainty-independent constraints. + * This includes bounds and constraints linking first- and second-stage variables. + * + * @param x candidate k-adaptable solution with k <= K + * @param K # of 2nd-stage policies + * @param xx structurally feasible solution with (possibly) reduced number of policies will be returned here (if feasible) + * @return true if feasible + */ + bool feasible_DET_K(const std::vector& x, const unsigned int K, std::vector& xx) const; + + /** + * Check if the given solution satisfies 1st-stage constraints involving uncertain parameters + * [to be used by internal routines only] + * + * @param x candidate solution + * @param q scenario from the uncertainty set will be returned here + * @return true if solution is feasible + */ + bool feasible_XQ(const std::vector& x, std::vector& q) const; + + /** + * Check if the given solution satisfies 1st-stage constraints involving uncertain parameter vectors contained in 'samples' + * [to be used by internal routines only] + * + * @param x candidate solution + * @param samples collection of scenarios to check feasibility against + * @param label sequence number of violating scenario in samples will be returned here (if not feasible) + * @return true if solution is feasible + */ + bool feasible_XQ(const std::vector& x, const std::vector >& samples, int & label) const; + + /** + * Check if the given K-adaptable solution satisfies 2nd-stage constraints involving uncertain parameters + * [to be used by internal routines only] + * + * @param x candidate solution + * @param K # of 2nd-stage policies + * @param q scenario from the uncertainty set will be returned here + * @return true if solution is feasible + */ + bool feasible_YQ(const std::vector& x, const unsigned int K, std::vector& q) const; + + /** + * Check if the given K-adaptable solution satisfies 2nd-stage constraints involving uncertain parameter vectors contained in 'samples' + * [to be used by internal routines only] + * + * @param x candidate solution + * @param K # of 2nd-stage policies + * @param samples collection of scenarios to check feasibility against + * @param label sequence number of violating scenario in samples will be returned here (if not feasible) + * @return true if feasible + */ + bool feasible_YQ(const std::vector& x, const unsigned int K, const std::vector >& samples, int & label, bool heur = 0) const; + + + + + + +public: + + /** + * Check if the given solution is deterministically feasible for the given scenario + * @param x candidate solution + * @param q candidate scenario + * @return true if feasible + */ + bool feasible_DET(const std::vector& x, const std::vector& q) const; + + /** + * Check if the given solution is static robust feasible + * @param x candidate solution + * @param q scenario from the uncertainty set will be returned here (if not feasible) + * @return true if feasible + */ + bool feasible_SRO(const std::vector& x, std::vector& q) const; + + /** + * Check if the given solution is feasible for the collection of scenarios contained in 'samples' + * @param x candidate solution + * @param samples collection of scenarios to check feasibility against + * @param label sequence number of violating scenario in samples will be returned here (if not feasible) + * @return true if feasible + */ + bool feasible_SRO(const std::vector& x, const std::vector >& samples, int & label) const; + + /** + * Check if the given solution is feasible for the K-adaptable problem + * @param x candidate k-adaptable solution with k <= K + * @param K # of 2nd-stage policies + * @param q scenario from the uncertainty set will be returned here (if not feasible) + * @return true if feasible + */ + bool feasible_KAdaptability(const std::vector& x, const unsigned int K, std::vector& q) const; + + /** + * Check if the given solution is K-adaptable feasible for the collection of scenarios contained in 'samples' + * @param x candidate k-adaptable solution with k <= K + * @param K # of 2nd-stage policies + * @param samples collection of scenarios to check feasibility against + * @param label sequence number of violating scenario in samples will be returned here (if not feasible) + * @return true if feasible + */ + bool feasible_KAdaptability(const std::vector& x, const unsigned int K, const std::vector >& samples, int & label) const; + + + + + + +public: + /** + * Solve the deterministic problem for the given scenario + * @param q scenario from the uncertainty set + * @param x optimal solution (if any) will be returned here + * @return solve status (non-zero value indicates unsuccessful termination) + */ + int solve_DET(const std::vector& q, std::vector& x) const; + + /** + * Solve the static robust problem via duality-based reformulation + * @param x optimal solution (if any) will be returned here + * @return solve status (non-zero value indicates unsuccessful termination) + */ + int solve_SRO_duality(std::vector& x) const; + + /** + * Solve the static robust problem via cutting plane method implemented using solver callbacks + * @param x optimal solution (if any) will be returned here + * @return solve status (non-zero value indicates unsuccessful termination) + */ + int solve_SRO_cuttingPlane(std::vector& x); + + /** + * Solve the scenario-based static robust problem (uncertainty set is a finite set of points) + * @param samples collection of scenarios to be robust against + * @param x optimal solution (if any) will be returned here + * @return solve status (non-zero value indicates unsuccessful termination) + */ + int solve_ScSRO(const std::vector >& samples, std::vector& x) const; + + + + + + +public: + /** Library of samples (temporary var -- to be used by solve_KAdaptability() only) */ + std::vector > bb_samples; + + /** Feasible (ideally optimal) static robust solution -- to be used by solve_KAdaptability() only) */ + std::vector xstatic; + + /** (K - 1)-adaptable solution -- to be used by solve_KAdaptability() in heuristic_mode only */ + std::vector xfix; + + /** Run in heuristic mode -- to be used by solve_KAdaptability() only */ + bool heuristic_mode; + + /** Current value of K in K-Adaptability -- to be used by solve_KAdaptability() only) */ + unsigned int NK; + + /** + * Get best solution found + * @return best solution found + */ + inline const std::vector& getX() const { + return xsol; + } + + /** + * Get nominal scenario + * @return nominal scenario from uncertainty set + */ + inline std::vector getNominal() const { + return pInfo->getNominal(); + } + + /** + * Is the underlying problem a pure LP? + * @return true if the problem is continuous + */ + inline bool isContinuous() const { + return pInfo->isContinuous(); + } + + /** + * Indicates if problem has 2nd-stage decisions only + * @return true if problem has 2nd-stage decisions only + */ + inline bool isSecondStageOnly() const { + return pInfo->isSecondStageOnly(); + } + + /** + * Indicates if problem has objective uncertainty only + * @return true if problem has objective uncertainty only + */ + inline bool hasObjectiveUncOnly() const { + return pInfo->hasObjectiveUncOnly(); + } + + /** + * Solve the K-adaptability problem using solver callbacks + * @param K # of 2nd-stage policies + * @param h run in heuristic mode + * @param x optimal solution (if any) will be returned here + * @return solve status (non-zero value indicates unsuccessful termination) + */ + int solve_KAdaptability(const unsigned int K, const bool h, std::vector& x); + + /** + * Solve the separation problem arising in solve_KAdaptability() + * @param x candidate solution + * @param K # of 2nd-stage policies + * @param label sequence number of violating scenario in bb_samples will be returned here + * @param heur separation is not guaranteed to be exact if this parameter is true + * @return true if separating scenario was found, false otherwise + */ + bool solve_separationProblem(const std::vector& x, const unsigned int K, int& label, bool heur = 0); + + /** + * Solve the min-max-min robust combinatorial optimization problem (as in Buchheim & Kurtz) + * @param K # of 2nd-stage policies + * @param x optimal solution (if any) will be returned here + * @return solve status (non-zero value indicates unsuccessful termination) + */ + int solve_KAdaptability_minMaxMin(const unsigned int K, std::vector& x); + +}; + +#endif \ No newline at end of file diff --git a/inc/uncertainty.hpp b/inc/uncertainty.hpp new file mode 100644 index 0000000..09190c1 --- /dev/null +++ b/inc/uncertainty.hpp @@ -0,0 +1,260 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + + +//////////////////////////////////////////// +// File: Uncertainty.hpp // +//////////////////////////////////////////// +// Class that represents a generic // +// uncertainty set // +//////////////////////////////////////////// + + + +#ifndef UNCERTAINTY_HPP +#define UNCERTAINTY_HPP + + +#include +#include +#include +#include + + + + +/** + * Generic uncertainty set + * All arrays are 1-indexed (indexing starts from 1, not 0) + */ +class UncertaintySet { +protected: + + /** # of uncertain parameters */ + int N; + + /** Nominal realization */ + std::vector nominal; + + /** Lower bounds on uncertain parameters */ + std::vector low; + + /** Upper bounds on uncertain parameters */ + std::vector high; + + /** Constraint matrix defining uncertainty set (including bounds) */ + std::vector > polytope_W; + + /** Matrix of zeros (backward compatibility) */ + std::vector > polytope_V; + + /** Corresponding right hand side vector */ + std::vector polytope_h; + + /* Corresponding sense of each constraint */ + std::vector polytope_sense; + + /** Solver environment to carry out optimizations on uncertainty set */ + CPXENVptr env; + + /** Solver problem object to carry out above optimizations */ + CPXLPptr lp; + + + +public: + + /** + * Constructs an UncertaintySet object + */ + UncertaintySet(); + + /** + * Copy constructor + */ + UncertaintySet(const UncertaintySet& U); + + /** + * Copy assignment operator + */ + UncertaintySet& operator=(const UncertaintySet& U); + + /** + * Free memory allocated + */ + ~UncertaintySet() noexcept(false); + + /** + * Clear all data structures that this uncertainty set holds + */ + void clear(); + + /** + * Add an uncertain parameter to the set + * @param nom the nominal value of the parameter + * @param lo the lower bound of the parameter + * @param hi the upper bound of the parameter + */ + void addParam(const double nom, const double lo, const double hi); + + /** + * Add a single facet to the uncertainty set + * @param input the list of uncertain parameters and their coefficients in the facet + * @param sense the sense of the constraint ('L', 'G', 'E') + * @param rhs the rhs of the facet + */ + void addFacet(const std::vector >& input, const char sense, const double rhs); + + /** + * Compute the maximum of a linear function of uncertain parameters over the uncertainty set + * @param input the list of uncertain parameters and their coefficients in the linear function + * @param result argmax is returned here (1-indexed) + * @return the result of the optimization (objective) + */ + double max(const std::vector >& input, std::vector& result) const; + + /** + * Compute the minimum of a linear function of uncertain parameters over the uncertainty set + * @param input the list of uncertain parameters and their coefficients in the linear function + * @param result argmin is returned here (1-indexed) + * @return the result of the optimization (objective) + */ + double min(const std::vector >& input, std::vector& result) const; + + /** + * Overload of previous function (to maintain backward compatibility) + * @param ind the list of indices of the uncertain parameters + * @param coef the list of coefficients of the above uncertain parameters + * @return the result of the optimization (objective) + */ + double getMaximumValue(const std::vector& ind, const std::vector& coef) const; + + /** + * This function adds an entire set of (empty) columns to the given lp object. + * These columns correspond to dual variables associated with + * each of the constraints of the inner optimization problem (to be dualized), + * i.e., constraints of the uncertainty set. + * + * THIS FUNCTION MUST BE CALLED WHENEVER A CONSTRAINT (WITH UNCERTAIN PARAMETERS) + * IS TO BE DUALIZED AND ADDED TO THE MAIN MILP. + * + * @param env solver environment in which master MILP resides + * @param lp solver model object in which master MILP resides + * @param dualName name of constraint to be attached to new columns + * @return number of new dual variables added + */ + int addVariables_DualVars(CPXCENVptr env, CPXLPptr lp, const std::string& dualName = "") const; + + /** + * Get nominal realization + * @return the nominal realization + */ + inline std::vector getNominal() const { return nominal; } + + /** + * Get lower bounds of uncertain parameters + * @return the lower bounds + */ + inline std::vector getLowerBounds() const { return low; } + + /** + * Get upper bounds of uncertain parameters + * @return the upper bounds + */ + inline std::vector getUpperBounds() const { return high; } + + /** + * Get clone of solver model object representing uncertainty set + * @param env_ pointer to environment to be used in conjunction with cloned solver model object + * @param stat pointer to solve status of clone operation + * @return pointer to clone of solver model object + */ + inline CPXLPptr getLPObject(CPXCENVptr env_, int *stat) const { return CPXXcloneprob(env_, lp, stat); } + + /** + * Get clone of solver model object representing uncertainty set + * @param stat pointer to solve status of clone operation + * @return pointer to clone of solver model object + */ + inline CPXLPptr getLPObject(int *stat) const { return getLPObject(env, stat); } + + /** + * Get pointer to solver environment to carry out optimizations on uncertainty set + * @return pointer to solver environment + */ + inline CPXENVptr getENVObject() const { return env; } + + /** + * Is the uncertainty set empty? + * @return bool indicating if uncertainty set is empty + */ + inline bool isUncertain() const { return (N != 0); } + + /** + * Get number of constraints defining uncertainty set (other than bound constraints) + * @return number of facets of uncertainty set + */ + inline int getNoOfFacets() const { return polytope_h.size(); } + + /** + * Get number of uncertain parameters + * @return number of uncertain parameters + */ + inline int getNoOfUncertainParameters() const { return N; } + + /** + * Get index of uncertain parameter inside uncertainty set + * @param char character specifying uncertain parameter + * @param ind index of uncertain parameter + * @return index of uncertain parameter + */ + inline int getParamIndex(const char, const int ind) const { return ind; } + + /** + * Get number of factors (backward compatibility) + * @return number of factors + */ + inline int getNumBudgetsFactors() const { return 0; } + + /** + * Get constraint matrix W in Wq [sense] h + * @return constraint matrix W + */ + inline std::vector > getMatrixW() const { return polytope_W; } + + /** + * Get constraint matrix V + * @return constraint matrix V + */ + inline std::vector > getMatrixV() const { return polytope_V; } + + /** + * Get vector h in Wq [sense] h + * @return right hand side vector h + */ + inline std::vector getMatrixH() const { return polytope_h; } + + /** + * Get vector h in Wq [sense] h + * @return right hand side vector h + */ + inline std::vector getMatrixSense() const { return polytope_sense; } +}; + +typedef const UncertaintySet* UNCSetCPtr; +typedef UncertaintySet* UNCSetPtr; + + +#endif + + diff --git a/src/indexingTools.cpp b/src/indexingTools.cpp new file mode 100644 index 0000000..1dfcf90 --- /dev/null +++ b/src/indexingTools.cpp @@ -0,0 +1,529 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Nikolaos Lappas and Chrysanthos Gounaris */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ +// Variable Indexing and Querying Class +// supporting up to 4 indices + +#include "indexingTools.h" +#include +#include +//################################ +static inline void abortProgram(const std::string &customMessage = "", const std::string &doNotUseThis = ""){ +#define abortProgram(...) abortProgram(std::string(__FILE__) + " -- Line " + std::to_string(__LINE__), std::string(__VA_ARGS__)) + std::cerr << " == ABORTING PROGRAM >>> " << (doNotUseThis.empty() ? customMessage : doNotUseThis) << " == " << std::endl; + abort(); + return; +} +//################################ + + + +bool VarInfo::isVarTypeConsistent(const int & ind1, const int & ind2, const int & ind3, const int & ind4, const int & ind5){ + if(ind1 < 0) return false; + if(ind2 > -1 && ind1 <0 ) return false; + if((ind3 > -1 && ind2 <0) || (ind3 > -1 && ind1 <0)) return false; + if((ind4 > -1 && ind3 < 0) || (ind4 > -1 && ind2 <0) || (ind4 > -1 && ind1 <0)) return false; + if((ind5 > -1 && ind4 < 0) || (ind5 > -1 && ind3 <0) || (ind5 > -1 && ind2 <0) || (ind5 > -1 && ind1 <0)) return false; + return true; +} + + +void VarInfo::linTo1dIndex(const int & linIndex, const std::tuple & VarType, int & ind1) const { + int diff = linIndex - std::get<0>(VarType); assert(diff >= 0); assert(diff <= std::get<1> (VarType)); + ind1 = diff; +} + +void VarInfo::linTo2dIndex(const int & linIndex, const std::tuple & VarType, int & ind1, int & ind2) const { + int diff = linIndex - std::get<0>(VarType); assert(diff >= 0); assert(diff <= std::get<1> (VarType)); + ind1 = diff / (std::get<4>(VarType) ); diff -= ind1 * std::get<4>(VarType); + ind2 = diff; +} + +void VarInfo::linTo3dIndex(const int & linIndex, const std::tuple & VarType, int & ind1, int & ind2, int & ind3) const { + int diff = linIndex - std::get<0>(VarType); assert(diff >= 0); assert(diff <= std::get<1> (VarType)); + ind1 = diff / (std::get<4>(VarType) * std::get<5>(VarType) ); diff -= ind1 * std::get<4>(VarType) * std::get<5>(VarType); + ind2 = diff / std::get<5>(VarType); diff -= ind2 * std::get<5>(VarType); + ind3 = diff; +} + +void VarInfo::linTo4dIndex(const int & linIndex, const std::tuple & VarType, int & ind1, int & ind2, int & ind3, int & ind4) const { + int diff = linIndex - std::get<0>(VarType); assert(diff >= 0); assert(diff <= std::get<1> (VarType)); + ind1 = diff / (std::get<4>(VarType) * std::get<5>(VarType) * std::get<6>(VarType) ); diff -= ind1 * std::get<4>(VarType) * std::get<5>(VarType) * std::get<6>(VarType); + ind2 = diff / (std::get<5>(VarType) * std::get<6>(VarType) ); diff -= ind2 * std::get<5>(VarType) * std::get<6>(VarType); + ind3 = diff / std::get<6>(VarType); diff -= ind3 * std::get<6>(VarType); + ind4 = diff; +} + +void VarInfo::linTo5dIndex(const int & linIndex, const std::tuple & VarType, int & ind1, int & ind2, int & ind3, int & ind4, int & ind5) const { + int diff = linIndex - std::get<0>(VarType); assert(diff >= 0); assert(diff <= std::get<1> (VarType)); + ind1 = diff / (std::get<4>(VarType) * std::get<5>(VarType) * std::get<6>(VarType) * std::get<7>(VarType) ); diff -= ind1 * std::get<4>(VarType) * std::get<5>(VarType) * std::get<6>(VarType) * std::get<7>(VarType); + ind2 = diff / (std::get<5>(VarType) * std::get<6>(VarType) * std::get<7>(VarType) ); diff -= ind2 * std::get<5>(VarType) * std::get<6>(VarType) * std::get<7>(VarType); + ind3 = diff / (std::get<6>(VarType) * std::get<7>(VarType) ); diff -= ind3 * std::get<6>(VarType) * std::get<7>(VarType); + ind4 = diff / std::get<7>(VarType); diff -= ind4 * std::get<7>(VarType); + ind5 = diff; +} + +void VarInfo::linToIndex(const int & linIndex, const std::tuple & VarType, int & ind1, int & ind2, int & ind3, int & ind4, int & ind5) const { + int dimensions = std::get<2> (VarType); + if(dimensions == 1) { + linTo1dIndex(linIndex, VarType, ind1); + ind2 = -1; + ind3 = -1; + ind4 = -1; + ind5 = -1; + } + else if(dimensions == 2) { + linTo2dIndex(linIndex, VarType, ind1, ind2); + ind3 = -1; + ind4 = -1; + ind5 = -1; + } + else if(dimensions == 3) { + linTo3dIndex(linIndex, VarType, ind1, ind2, ind3); + ind4 = -1; + ind5 = -1; + } + else if(dimensions == 4) { + linTo4dIndex(linIndex, VarType, ind1, ind2, ind3, ind4); + ind5 = -1; + } + else if(dimensions == 5) { + linTo5dIndex(linIndex, VarType, ind1, ind2, ind3, ind4, ind5); + } + else { + abortProgram("invalid number of dimensions for this Variable"); + } +} + +int VarInfo::xdIndicesToLinIndex(const std::tuple & VarType, const int & ind1, const int & ind2, const int & ind3, const int & ind4, const int & ind5) const { + int index; + int dimensions = std::get<2> (VarType); + if(dimensions == 1) { + assert(ind1 > -1); + assert(ind2 < 0 && ind3 < 0 && ind4 < 0 && ind5 < 0); + index = ind1; + } + else if(dimensions == 2) { + assert(ind1 > -1 && ind2 > -1); + assert(ind3 < 0 && ind4 < 0 && ind5 < 0); + index = ind1 * std::get<4>(VarType) + ind2; + } + else if(dimensions == 3) { + assert(ind1 > -1 && ind2 > -1 && ind3 > -1); + assert(ind4 < 0 && ind5 < 0); + index = ind1 * std::get<4>(VarType) * std::get<5>(VarType) + ind2 * std::get<5>(VarType) + ind3; + } + else if(dimensions == 4) { + assert(ind1 > -1 && ind2 > -1 && ind3 > -1 && ind4 > -1); + assert(ind5 < 0); + index = ind1 * std::get<4>(VarType) * std::get<5>(VarType) * std::get<6>(VarType) + ind2 * std::get<5>(VarType) * std::get<6>(VarType) + ind3 * std::get<6>(VarType) + ind4; + } + else if(dimensions == 5) { + assert(ind1 > -1 && ind2 > -1 && ind3 > -1 && ind4 > -1 && ind5 > -1); + index = ind1 * std::get<4>(VarType) * std::get<5>(VarType) * std::get<6>(VarType) * std::get<7>(VarType) + ind2 * std::get<5>(VarType) * std::get<6>(VarType) * std::get<7>(VarType) + ind3 * std::get<6>(VarType) * std::get<7>(VarType) + ind4 * std::get<7>(VarType) + ind5; + } + else { + abortProgram("invalid number of dimensions for this Variable (xdIndicesToLinIndex)"); + } + + index += std::get<0>(VarType); + assert(index <= std::get<1>(VarType)); + return index; +} + +// Clear the object +void VarInfo::clear() { + VarType.clear(); + totalVars = 0; + UpperBound.clear(); + LowerBound.clear(); + ColumnType.clear(); + ObjCoefficient.clear(); + UndefinedVar.clear(); + UndefinedVarCount.clear(); +} + +// Add a new variable type +void VarInfo::addVarType(const std::string & name, const char & type, const double & lb, const double & ub, const int & ind1, const int & ind2, const int & ind3, const int & ind4, const int & ind5) { + int first = 0; + if(VarType.size() > 0) { + first = totalVars; + } + if (type != 'B' && type != 'C') abortProgram("Column type must be 'B' or 'C'"); + // Calculate the size + int size = 0; + int numDimensions = 0; + if(ind1 > -1) { + size = ind1; + numDimensions++; + if(ind2 > -1 ) { + size *= ind2; + numDimensions++; + if(ind3 > -1) { + size *= ind3; + numDimensions++; + if(ind4 > -1 ) { + size *= ind4; + numDimensions++; + if(ind5 > -1 ) { + size *= ind5; + numDimensions++; + } + } + } + } + } + totalVars = first + size; + int last = first + size -1; + auto temp = std::make_tuple(first, last, numDimensions, ind1, ind2, ind3, ind4, ind5); + VarType.emplace(name, temp); + UpperBound.resize(totalVars, ub); + LowerBound.resize(totalVars, lb); + ColumnType.resize(totalVars, type); + ObjCoefficient.resize(totalVars, 0.0); + UndefinedVar.resize(totalVars, 0); + UndefinedVarCount.resize(totalVars, (UndefinedVarCount.empty() ? 0 : UndefinedVarCount.back()) ); + + if(isVarTypeConsistent(ind1, ind2, ind3, ind4, ind5) != true) abortProgram("Var Type not valid"); +} + +// // set the vector of undefined vars +// void VarInfo::setUndefinedVars(const std::vector & Undefined) { +// UndefinedVar.clear(); UndefinedVar.resize(totalVars, 0); +// assert((int)Undefined.size() == totalVars); +// UndefinedVar = Undefined; +// UndefinedVarCount.resize(totalVars, 0); +// if (totalVars > 0) { +// UndefinedVarCount[0] = static_cast(UndefinedVar[0]); +// } +// for (int i = 1; i < totalVars; i++) { +// UndefinedVarCount[i] = UndefinedVarCount[i-1] + static_cast(UndefinedVar[i]); +// assert(std::accumulate(UndefinedVar.begin(), UndefinedVar.begin() + i + 1, int(0)) == UndefinedVarCount[i]); +// } +// } + +// // set the vector with the LB of the vars +// void VarInfo::setVarsLB(const std::vector & Lbounds) { +// LowerBound.clear(); LowerBound.resize(totalVars, 0.0); +// assert((int)Lbounds.size() == totalVars); +// LowerBound = Lbounds; +// } + +// // set the vector with the UB of the vars +// void VarInfo::setVarsUB(const std::vector & Ubounds) { +// UpperBound.clear(); UpperBound.resize(totalVars, 0.0); +// assert((int)Ubounds.size() == totalVars); +// UpperBound = Ubounds; +// } + +// // set the vector with the column types +// void VarInfo::setVarsColType(const std::vector & Ctype) { +// ColumnType.clear(); ColumnType.resize(totalVars, 0.0); +// assert((int)Ctype.size() == totalVars); +// ColumnType = Ctype; +// } + +// // set the vector with the objective coefficients +// void VarInfo::setVarsObjCoeff(const std::vector & Ocoeff) { +// ObjCoefficient.clear(); ObjCoefficient.resize(totalVars, 0.0); +// assert((int)Ocoeff.size() == totalVars); +// ObjCoefficient = Ocoeff; +// } + +// set the given var to be undefined +void VarInfo::setUndefinedVar(const int & index) { + assert(index > -1 && index < totalVars); + if (!UndefinedVar[index]) { + UndefinedVar[index] = 1; + for (int i = index; i < totalVars; i++) UndefinedVarCount[i]++; + } +} + +// set the UB of the given var +void VarInfo::setVarUB(const double & Ub, const int & index) { + assert(index > -1 && index < totalVars); + UpperBound[index] = Ub; +} + +// set the LB of the given var +void VarInfo::setVarLB(const double & Lb, const int & index) { + assert(index > -1 && index < totalVars); + LowerBound[index] = Lb; +} + +// set the column type of the given var +void VarInfo::setVarColType(const char & Ctype, const int & index) { + assert(index > -1 && index < totalVars); + ColumnType[index] = Ctype; +} + +// set the objective coefficient of the given var +void VarInfo::setVarObjCoeff(const double & Ocoeff, const int & index) { + assert(index > -1 && index < totalVars); + ObjCoefficient[index] = Ocoeff; +} + +// set the given var to be undefined +void VarInfo::setUndefinedVar(const std::string & type, const int ind1, const int ind2, const int ind3, const int ind4, const int ind5) { + setUndefinedVar(getVarLinIndex(type, ind1, ind2, ind3, ind4, ind5)); +} + +// set the UB of the given var +void VarInfo::setVarUB(const double & Ub, const std::string & type, const int ind1, const int ind2, const int ind3, const int ind4, const int ind5) { + setVarUB(Ub, getVarLinIndex(type, ind1, ind2, ind3, ind4, ind5)); +} + +// set the LB of the given var +void VarInfo::setVarLB(const double & Lb, const std::string & type, const int ind1, const int ind2, const int ind3, const int ind4, const int ind5) { + setVarLB(Lb, getVarLinIndex(type, ind1, ind2, ind3, ind4, ind5)); +} + +// set the column type of the given var +void VarInfo::setVarColType(const char & Ctype, const std::string & type, const int ind1, const int ind2, const int ind3, const int ind4, const int ind5) { + setVarColType(Ctype, getVarLinIndex(type, ind1, ind2, ind3, ind4, ind5)); +} + +// set the objective coefficient of the given var +void VarInfo::setVarObjCoeff(const double & Ocoeff, const std::string & type, const int ind1, const int ind2, const int ind3, const int ind4, const int ind5) { + setVarObjCoeff(Ocoeff, getVarLinIndex(type, ind1, ind2, ind3, ind4, ind5)); +} + + +// get the number of vars that this object holds +int VarInfo::getTotalVarSize() const { + return totalVars; +} + +// get the number of defined vars that this object holds +int VarInfo::getTotalDefVarSize() const { + return totalVars - (totalVars > 0 ? UndefinedVarCount.back() : 0); +} + +// get the name and the indices, given the linear index +void VarInfo::getVarInfo(const int & linIndex, std::string &type, int & ind1, int & ind2, int & ind3, int & ind4, int & ind5) const { + bool found = false; + for(auto it = VarType.cbegin(); it != VarType.cend(); ++it) { + int firstOf = std::get<0>((*it).second); + int lastOf = std::get<1>((*it).second); + if(linIndex >= firstOf && linIndex <= lastOf) { + found = true; + type = (*it).first; + linToIndex(linIndex, (*it).second, ind1, ind2, ind3, ind4, ind5); + } + } + if(!found) abortProgram("Variable not found (getVarInfo)"); +} + +// get the linear index, given the name and the indices +int VarInfo::getVarLinIndex(const std::string & type, const int ind1, const int ind2, const int ind3, const int ind4, const int ind5) const { + if(VarType.find(type) == VarType.end()) { // key not found + abortProgram("Unknown Variable key (getVarLinIndex)"); + return 0; + } + else { + return xdIndicesToLinIndex(VarType.at(type), ind1, ind2, ind3, ind4, ind5); + } +} + +// get the linear index of defined variable, given the name and the indices +int VarInfo::getDefVarLinIndex(const std::string & type, const int ind1, const int ind2, const int ind3, const int ind4, const int ind5) const { + if(VarType.find(type) == VarType.end()) { // key not found + abortProgram("Unknown Variable key (getDefVarLinIndex)"); + } + const int linIndex = xdIndicesToLinIndex(VarType.at(type), ind1, ind2, ind3, ind4, ind5); + assert(!UndefinedVar[linIndex]); + return linIndex - UndefinedVarCount[linIndex]; +} + +// get the linear index of defined variable, given the original index +int VarInfo::getDefVarLinIndex(const int & index) const { + assert(index > -1 && index < totalVars); + return index - UndefinedVarCount[index]; +} + +// get the number of variables of the requested type +int VarInfo::getVarTypeSize(const std::string & name) const { + if(VarType.find(name) == VarType.end()) { + abortProgram("Invalid key (getVarTypeSize)"); + } + + const int elements = std::get<1>(VarType.at(name)) - std::get<0>(VarType.at(name)) + 1; + return elements; +} + +// get the number of defined variables of the requested type +int VarInfo::getDefVarTypeSize(const std::string & name) const { + if(VarType.find(name) == VarType.end()) { + abortProgram("Invalid key (getVarTypeSize)"); + } + const int last = std::get<1>(VarType.at(name)); + const int first = std::get<0>(VarType.at(name)); + int elements = last - first + 1 - UndefinedVarCount[last]; + if (first >= 1) { + elements += UndefinedVarCount[first - 1]; + } + return elements; +} + +// get the linear index of the first var of a given var type (iterating tool) +int VarInfo::getFirstOfVarType(const std::string & name) const { + if(VarType.find(name) == VarType.end()) abortProgram("invalid var key (getFirstOfVarType)"); + return std::get<0> (VarType.at(name)); +} + +// get the linear index of the last var of a given var type (iterating tool) +int VarInfo::getLastOfVarType(const std::string & name) const { + if(VarType.find(name) == VarType.end()) abortProgram("invalid var key (getLastOfVarType)"); + return std::get<1> (VarType.at(name)); +} + +// get the linear index of the first defined var of a given var type (iterating tool) +int VarInfo::getFirstDefOfVarType(const std::string & name) const { + if(VarType.find(name) == VarType.end()) abortProgram("invalid var key (getFirstOfDefVarType)"); + int first = std::get<0> (VarType.at(name)); + const int last = std::get<1> (VarType.at(name)); + while (UndefinedVar[first] && first <= last) first++; + if (first > last) abortProgram("All variables of type are undefined (getFirstOfDefVarType)"); + return getDefVarLinIndex(first); +} + +// get the linear index of the last defined var of a given var type (iterating tool) +int VarInfo::getLastDefOfVarType(const std::string & name) const { + if(VarType.find(name) == VarType.end()) abortProgram("invalid var key (getLastOfLastVarType)"); + int last = std::get<1> (VarType.at(name)); + const int first = std::get<0> (VarType.at(name)); + while(UndefinedVar[last] && last >= first) last--; + if (last < first) abortProgram("All variables of type are undefined (getLastOfDefVarType)"); + return getDefVarLinIndex(last); +} + +// get the variable name, given the linear index +std::string VarInfo::getVarName(const int & index) const { + bool found = false; + std::string type; + for(auto it = VarType.cbegin(); it != VarType.cend(); ++it) { + int firstOf = std::get<0>((*it).second); + int lastOf = std::get<1>((*it).second); + if(index >= firstOf && index <= lastOf) { + found = true; + type = (*it).first; + return type; + } + } + if(!found) abortProgram("Variable not found (getVarName)"); + return ""; +} + +// query if a var is undefined +bool VarInfo::isUndefVar(const int & index) const { + assert(index > -1 && index < totalVars); + if(UndefinedVar.size() == 0) { + return false; // default is 0 + } + else { + return UndefinedVar[index]; + } +} + +// get the UB of a specific var +double VarInfo::getVarUB(const int & index) const { + assert(index > -1 && index < totalVars); + if(UpperBound.size() == 0) { + return 0.0; // default bound is 0 + } + else { + return UpperBound[index]; + } +} + + +// get the LB of a specific var +double VarInfo::getVarLB(const int & index) const { + assert(index > -1 && index < totalVars); + if(LowerBound.size() == 0) { + return 0.0; // default bound is 0 + } + else { + return LowerBound[index]; + } +} + + +// get the type of a specific var, given the linear index as returned by getVarLinIndex() +char VarInfo::getVarColType(const int & index) const { + assert(index > -1 && index < totalVars); + if(ColumnType.size() == 0) { + return 'C'; // default type is 'C' + } + else { + return ColumnType[index]; + } +} + +// get the LB of a specific var +double VarInfo::getVarObjCoeff(const int & index) const { + assert(index > -1 && index < totalVars); + if(ObjCoefficient.size() == 0) { + return 0.0; // default coefficient is 0 + } + else { + return ObjCoefficient[index]; + } +} + + + + + +// void addToCPX(const VarInfo& X) { + +// } + + + + +// int main(int argc, char** argv) { + +// VarInfo X; +// std::vector UV; +// X.addVarType("z", 1); UV.emplace_back(0); +// X.addVarType("x", 6, 6); +// for (int i = 0; i <= 5; i ++) for (int j = 0; j <= 5; j++) { +// UV.emplace_back(0); +// X.setVarUB(X.getVarLinIndex("x",i,j), 1.0); +// if (i == j || (i == 5) || (j==0) || (i==0 && j==5)) UV.back() = 1; +// } +// X.setUndefinedVars(UV); + + +// for (int i = 0; i <= 5; i ++) for (int j = 0; j <= 5; j++) if (!X.isUndefVar(X.getVarLinIndex("x", i, j))) { +// std::cout << "(" << i << "," << j << ") = " << X.getDefVarLinIndex("x", i, j) << "\n"; +// } + +// std::cout << X.getDefVarLinIndex("z", 0) << std::endl; + +// for (int i = 0; i < X.getTotalVarSize(); i++) if (!X.isUndefVar(i)) { +// std::cout << X.getVarName(i) << "(" << X.getDefVarLinIndex(i) << ") in [" << X.getVarLB(i) << "," << X.getVarUB(i) <<"]\n"; +// } + +// X.addVarType("y", 0); + +// std::cout << X.getFirstOfVarType("x") << " " << X.getFirstOfVarType("z") << "\n"; +// std::cout << X.getFirstDefOfVarType("x") << " " << X.getFirstDefOfVarType("z") << "\n"; +// std::cout << X.getLastOfVarType("x") << " " << X.getLastOfVarType("z") << "\n"; +// std::cout << X.getLastDefOfVarType("x") << " " << X.getLastDefOfVarType("z") << "\n"; +// std::cout << X.getDefVarTypeSize("x") << " " << X.getDefVarTypeSize("z") << " " << X.getDefVarTypeSize("y") << "\n"; +// std::cout << X.getTotalVarSize() << " " << X.getTotalDefVarSize() << "\n"; + +// return 0; +// } + + diff --git a/src/problemInfo.cpp b/src/problemInfo.cpp new file mode 100644 index 0000000..899aaac --- /dev/null +++ b/src/problemInfo.cpp @@ -0,0 +1,108 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#include "problemInfo.hpp" +#include + + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo::resize(unsigned int K) { + assert(isConsistentWithDesign()); + unsigned int l = numPolicies; + if (numPolicies < K) { + numPolicies = K; + } + for (unsigned int k = l; k < numPolicies; k++) { + makeConsY(k); + } + assert(isConsistentWithDesign()); +} + +//----------------------------------------------------------------------------------- + + +bool KAdaptableInfo::isConsistentWithDesign() const { + if (numPolicies < 1) return 0; + if (C_XY.size() != numPolicies) return 0; + if (C_XYQ.size() != numPolicies) return 0; + if (B_Y.size() != numPolicies) return 0; + if (X.getTotalDefVarSize() <= 0) return 0; + if (Y.getTotalDefVarSize() < 0) return 0; + if (X.getTotalDefVarSize() != getNumFirstStage()) return 0; + if (Y.getTotalDefVarSize() != getNumSecondStage()) return 0; + if (existsFirstStage) { + if (getNumFirstStage() <= 1) return 0; + } + if (U.getNoOfUncertainParameters() == 0) { + if (!C_XQ.empty()) return 0; + for (unsigned int p = 0; p < numPolicies; p++) { + if (!C_XYQ[p].empty()) return 0; + } + } + if (objectiveUnc) { + if (!C_XQ.empty()) return 0; + for (unsigned int p = 0; p < numPolicies; p++) { + if (C_XYQ[p].size() != 1) return 0; + } + } + bool check = 0; + for (int x = 0; x < X.getTotalVarSize(); x++) { + if (X.getVarColType(x) == 'B') { + check = 1; + if (!hasInteger) return 0; + } + } + for (int y = 0; y < Y.getTotalVarSize(); y++) { + if (Y.getVarColType(y) == 'B') { + check = 1; + if (!hasInteger) return 0; + } + } + if (hasInteger && !check) return 0; + for (const auto& con: C_X) { + if (con.isEmpty()) return 0; + if (con.existBilinearTerms()) return 0; + if (con.existConstQTerms()) return 0; + } + for (const auto& con: C_XQ) { + if (con.isEmpty()) return 0; + if (!con.existBilinearTerms() && !con.existConstQTerms()) return 0; + } + for (const auto& con_K: C_XY) { + for (const auto& con: con_K) { + if (con.isEmpty()) return 0; + if (con.existBilinearTerms()) return 0; + if (con.existConstQTerms()) return 0; + } + } + for (const auto& con_K: C_XYQ) { + for (const auto& con: con_K) { + if (con.isEmpty()) return 0; + if (!con.existBilinearTerms() && !con.existConstQTerms()) return 0; + } + } + for (const auto& con: B_X) { + if (con.isEmpty()) return 0; + if (con.existBilinearTerms()) return 0; + if (con.existConstQTerms()) return 0; + } + for (const auto& con_K: B_Y) { + for (const auto& con: con_K) { + if (con.isEmpty()) return 0; + if (con.existBilinearTerms()) return 0; + if (con.existConstQTerms()) return 0; + } + } + return 1; +} \ No newline at end of file diff --git a/src/problemInfo_knp.cpp b/src/problemInfo_knp.cpp new file mode 100644 index 0000000..cdac595 --- /dev/null +++ b/src/problemInfo_knp.cpp @@ -0,0 +1,278 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#include "problemInfo_knp.hpp" +#include + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_KNP::makeUncSet() { + U.clear(); + + // define uncertain risk factors + for (unsigned int f = 1; f < data.phi[0].size(); ++f) { + U.addParam(0, -1, 1); + } +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_KNP::makeVars() { + X.clear(); + Y.clear(); + + // objective variable (considered to be 1st-stage) + // should always be defined and should always be declared first + X.addVarType("O", 'C', -CPX_INFBOUND, +CPX_INFBOUND, 1); + + // x(i) : invest in project i before observing risk factors + X.addVarType("x", 'B', 0, 1, 1 + data.N); + + // y(i) : invest in project i after observing risk factors + Y.addVarType("y", 'B', 0, 1, 1 + data.N); + + // All arrays/matrices are 1-indexed except if option to borrow + if (!data.loan) { + X.setUndefinedVar("x", 0); + Y.setUndefinedVar("y", 0); + } + else { + X.setVarColType('C', "x", 0); + Y.setVarColType('C', "y", 0); + X.setVarUB(+CPX_INFBOUND, "x", 0); + Y.setVarUB(+CPX_INFBOUND, "y", 0); + } + + // define objective variable + X.setVarObjCoeff(1.0, "O", 0); +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_KNP::makeConsX() { + ConstraintExpression temp; + + ///////// + // B_X // + ///////// + B_X.clear(); + + // bounds on x(0) + if (data.loan) { + temp.clear(); + temp.addTermX(getVarIndex_1("x", 0), 1); + temp.rowname("LB_x(0)"); + temp.sign('G'); + temp.RHS(0); + B_X.emplace_back(temp); + } + + // bounds on x(i) + for (int i = 1; i <= data.N; ++i) { + temp.clear(); + temp.addTermX(getVarIndex_1("x", i), 1); + + temp.rowname("LB_x(" + std::to_string(i) + ")"); + temp.sign('G'); + temp.RHS(0); + B_X.emplace_back(temp); + + temp.rowname("UB_x(" + std::to_string(i) + ")"); + temp.sign('L'); + temp.RHS(1); + B_X.emplace_back(temp); + } + + + + ///////// + // C_X // + ///////// + C_X.clear(); + + if (data.loan) { + temp.clear(); + temp.rowname("BUDGET_1"); + temp.sign('L'); + temp.RHS(data.B); + temp.addTermX(getVarIndex_1("x", 0), -1.0); + for (int i = 1; i <= data.N; ++i) if (data.cost[i] != 0.0) { + temp.addTermX(getVarIndex_1("x", i), data.cost[i]); + + for (unsigned int f = 1; f < data.phi[i].size(); f++) { + double coef = data.phi[i][f] * data.cost[i] * 0.5; + + if (coef != 0.0) { + temp.addTermX(getVarIndex_1("x", i), coef); + } + } + } + C_X.emplace_back(temp); + } + + + ////////// + // C_XQ // + ////////// + C_XQ.clear(); + + // nothing to do +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_KNP::makeConsY(unsigned int l) { + assert(C_XY.size() == B_Y.size()); + assert(C_XY.size() == C_XYQ.size()); + assert(numPolicies >= l); + if (l == 0) { + B_Y.clear(); + C_XY.clear(); + C_XYQ.clear(); + } + ConstraintExpression temp; + + ///////// + // B_Y // + ///////// + for (unsigned int k = B_Y.size(); k <= l; k++) { + if (B_Y.size() < k + 1) B_Y.resize(k + 1); + + B_Y[k].clear(); + + // bounds on y(0) + if (data.loan) { + temp.clear(); + temp.addTermX(getVarIndex_2(k, "y", 0), 1); + temp.rowname("LB_y(0," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(0); + B_Y[k].emplace_back(temp); + } + + // bounds on y(i) + for (int i = 1; i <= data.N; ++i) { + temp.clear(); + temp.addTermX(getVarIndex_2(k, "y", i), 1); + + temp.rowname("LB_y(" + std::to_string(i) + "," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(0); + B_Y[k].emplace_back(temp); + + temp.rowname("UB_y(" + std::to_string(i) + "," + std::to_string(k) + ")"); + temp.sign('L'); + temp.RHS(1); + B_Y[k].emplace_back(temp); + } + } + + + ////////// + // C_XY // + ////////// + for (unsigned int k = C_XY.size(); k <= l; k++) { + if (C_XY.size() < k + 1) C_XY.resize(k + 1); + + C_XY[k].clear(); + + // invest early or invest late + for (int i = 1; i <= data.N; ++i) { + temp.clear(); + temp.rowname("EITHER(" + std::to_string(i) + "," + std::to_string(k) + ")"); + temp.sign('L'); + temp.RHS(1.0); + temp.addTermX(getVarIndex_1("x", i), 1.0); + temp.addTermX(getVarIndex_2(k, "y", i), 1.0); + + C_XY[k].emplace_back(temp); + } + } + + + /////////// + // C_XYQ // + /////////// + for (unsigned int k = C_XYQ.size(); k <= l; k++) { + if (C_XYQ.size() < k + 1) C_XYQ.resize(k + 1); + + C_XYQ[k].clear(); + + // objective function + temp.clear(); + temp.rowname("OBJ_CONSTRAINT(" + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(0); + temp.addTermX(getVarIndex_1("O", 0), 1); + if (data.loan) temp.addTermX(getVarIndex_1("x", 0), -data.ell); + if (data.loan) temp.addTermX(getVarIndex_2(k, "y", 0), -(data.ell * data.lambda)); + for (int i = 1; i <= data.N; ++i) if (data.profit[i] != 0.0) { + temp.addTermX(getVarIndex_1("x", i), data.profit[i]); + temp.addTermX(getVarIndex_2(k, "y", i), data.theta * data.profit[i]); + + for (unsigned int f = 1; f < data.ksi[i].size(); f++) { + double coef = data.ksi[i][f] * data.profit[i] * 0.5; + + if (coef != 0.0) { + temp.addTermProduct(getVarIndex_1("x", i), U.getParamIndex('q', f), coef); + temp.addTermProduct(getVarIndex_2(k, "y", i), U.getParamIndex('q', f), data.theta * coef); + } + } + } + C_XYQ[k].emplace_back(temp); + + + // budget + temp.clear(); + temp.rowname("BUDGET(" + std::to_string(k) + ")"); + temp.sign('L'); + temp.RHS(data.B); + if (data.loan) temp.addTermX(getVarIndex_1("x", 0), -1.0); + if (data.loan) temp.addTermX(getVarIndex_2(k, "y", 0), -1.0); + for (int i = 1; i <= data.N; ++i) if (data.cost[i] != 0.0) { + temp.addTermX(getVarIndex_1("x", i), data.cost[i]); + temp.addTermX(getVarIndex_2(k, "y", i), data.cost[i]); + + for (unsigned int f = 1; f < data.phi[i].size(); f++) { + double coef = data.phi[i][f] * data.cost[i] * 0.5; + + if (coef != 0.0) { + temp.addTermProduct(getVarIndex_1("x", i), U.getParamIndex('q', f), coef); + temp.addTermProduct(getVarIndex_2(k, "y", i), U.getParamIndex('q', f), coef); + } + } + } + C_XYQ[k].emplace_back(temp); + } +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_KNP::setInstance(const KNP& d) { + data = d; + hasInteger = 1; + objectiveUnc = 0; + existsFirstStage = 1; + numFirstStage = 1 + data.N + (data.loan); + numSecondStage = data.N + (data.loan); + numPolicies = 1; + solfilename = data.solfilename; + + makeUncSet(); + makeVars(); + makeConsX(); + makeConsY(0); + assert(isConsistentWithDesign()); +} + +//----------------------------------------------------------------------------------- \ No newline at end of file diff --git a/src/problemInfo_psp.cpp b/src/problemInfo_psp.cpp new file mode 100644 index 0000000..e0ba1cf --- /dev/null +++ b/src/problemInfo_psp.cpp @@ -0,0 +1,400 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#include "problemInfo_psp.hpp" +#include + +//----------------------------------------------------------------------------------- + +#define FACTOR_ALPHA 0.5 +#define FACTOR_BETA 1.0 +#define DO_MAX 1 +#define SPECIAL_RESTRICTIVE_LDR 0 + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_PSP::makeUncSet() { + U.clear(); + + if (data.special) { + // compute k + int k = (data.N - 1)/3; + + // define primary uncertain parameters + for (int i = 1; i <= k; ++i) { + U.addParam(0.5, 0, 1); + } + + // define secondary parameters representing absolute deviation from 0.5 + for (int i = 1; i <= k; ++i) { + U.addParam(0, 0, 0.5); + } + + // define budget constraint + std::vector > constraint; + for (int i = k + 1; i <= (2 * k); ++i) { + constraint.emplace_back(i, 1); + } + U.addFacet(constraint, 'L', 0.5); + + // enforce absolute value definition + for (int i = 1; i <= k; i++) { + constraint.clear(); + constraint.emplace_back(i + k, 1); + constraint.emplace_back(i, -1); + U.addFacet(constraint, 'G', -0.5); + + constraint.clear(); + constraint.emplace_back(i + k, 1); + constraint.emplace_back(i, +1); + U.addFacet(constraint, 'G', +0.5); + } + } + else { + // define uncertain risk factors + for (unsigned int f = 1; f < data.phi[0].size(); ++f) { + U.addParam(0, -1, 1); + } + + // define beta-net-alpha + std::vector > constraint; + for (unsigned int f = 1; f < data.phi[0].size(); ++f) { + constraint.emplace_back(f, 1); + } + + U.addFacet(constraint, 'L', FACTOR_BETA * ((int)data.phi[0].size() - 1) ); + U.addFacet(constraint, 'G', -FACTOR_BETA * ((int)data.phi[0].size() - 1)); + } +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_PSP::makeVars() { + X.clear(); + Y.clear(); + + // objective variable (considered to be 1st-stage) + X.addVarType("O", 'C', -CPX_INFBOUND, +CPX_INFBOUND, 1); + + // x(i) : resource allocated to task i + if (!data.special) X.addVarType("x", 'C', 0, 0.5, 1 + data.N); + + // y(i) : start time of task i + Y.addVarType("y", 'C', (data.affine ? (data.special ? (SPECIAL_RESTRICTIVE_LDR ? 0 : -CPX_INFBOUND) : -CPX_INFBOUND) : 0), +CPX_INFBOUND, 1 + data.N); + + // gamma(i,j) : slopes of affine decision rules + if (data.affine) { + if (!data.special) { + Y.addVarType("g", 'C', -CPX_INFBOUND, +CPX_INFBOUND, data.N, data.phi[0].size()); + } else { + Y.addVarType("g", 'C', (SPECIAL_RESTRICTIVE_LDR ? 0 : -CPX_INFBOUND), +CPX_INFBOUND, data.N, 1 + ((data.N - 1) / 3)); + } + } + + // All arrays/atrices are 1-indexed + if (!data.special) X.setUndefinedVar("x", 0); + Y.setUndefinedVar("y", 0); + if (data.affine) { + for (int i = 0; i < data.N; i++) Y.setUndefinedVar("g", i, 0); + for (int f = 0, flim = (data.special ? (1 + ((data.N - 1)/3)) : data.phi[0].size()); f < flim; f++) Y.setUndefinedVar("g", 0, f); + } + + // cannot adjust uncertain durations of dummy nodes + if (data.affine) for (int f = 0, flim = (data.special ? (1 + ((data.N - 1)/3)) : data.phi[0].size()); f < flim; f++) { + Y.setUndefinedVar("g", 1, f); + } + + // define objective variable + X.setVarObjCoeff(1.0, "O", 0); + + // fix allocations of dummy nodes + if (!data.special) X.setVarUB(0, "x", 1); + if (!data.special) X.setVarUB(0, "x", data.N); + + // set early-start schedule + Y.setVarLB(0, "y", 1); + Y.setVarUB(0, "y", 1); + if (data.affine) { + Y.setVarLB(0, "y", data.N); + } + + // restrict slopes of affine decision rules + if (data.special && data.affine) { + int k = (data.N - 1)/3; + for (int i = 1; i <= k; i++) for (int f = i + 1; f <= k; f++) { + Y.setVarUB(0, "g", (3*i) - 1, f); + Y.setVarUB(0, "g", (3*i), f); + if (i != k) Y.setVarUB(0, "g", (3*i) + 1, f); + } + } +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_PSP::makeConsX() { + ConstraintExpression temp; + + ///////// + // B_X // + ///////// + B_X.clear(); + + // bounds on x(i) + if (!data.special) for (int i = 1; i <= data.N; ++i) { + temp.clear(); + temp.addTermX(getVarIndex_1("x", i), 1); + + temp.rowname("LB_x(" + std::to_string(i) + ")"); + temp.sign('G'); + temp.RHS(0); + B_X.emplace_back(temp); + + temp.rowname("UB_x(" + std::to_string(i) + ")"); + temp.sign('L'); + temp.RHS(0.5 * ((i > 1) && (i < data.N))); + B_X.emplace_back(temp); + } + + + + ///////// + // C_X // + ///////// + C_X.clear(); + + // Sum(x) <= B + if (!data.special) { + temp.clear(); + temp.rowname("BUDGET"); + temp.sign('L'); + temp.RHS(data.B); + for (int i = 1; i <= data.N; ++i) { + temp.addTermX(getVarIndex_1("x", i), 1); + } + C_X.emplace_back(temp); + } + + + + ////////// + // C_XQ // + ////////// + C_XQ.clear(); + + // nothing to do +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_PSP::makeConsY(unsigned int l) { + assert(C_XY.size() == B_Y.size()); + assert(C_XY.size() == C_XYQ.size()); + assert(numPolicies >= l); + if (l == 0) { + B_Y.clear(); + C_XY.clear(); + C_XYQ.clear(); + } + ConstraintExpression temp; + + ///////// + // B_Y // + ///////// + for (unsigned int k = B_Y.size(); k <= l; k++) { + if (B_Y.size() < k + 1) B_Y.resize(k + 1); + + B_Y[k].clear(); + + // early-start schedule + temp.clear(); + temp.addTermX(getVarIndex_2(k, "y", 1), 1); + + temp.rowname("LB_y(1," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(0); + B_Y[k].emplace_back(temp); + + temp.rowname("UB_y(1," + std::to_string(k) + ")"); + temp.sign('L'); + temp.RHS(0); + B_Y[k].emplace_back(temp); + + // non-negative makespan + temp.clear(); + temp.addTermX(getVarIndex_2(k, "y", data.N), 1); + temp.rowname("LB_y(" + std::to_string(data.N) + "," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(0); + B_Y[k].emplace_back(temp); + + // bounds on y(i) start times + if (!data.affine || (data.special && SPECIAL_RESTRICTIVE_LDR)) for (int i = 2; i < data.N; i++) { + temp.clear(); + temp.addTermX(getVarIndex_2(k, "y", i), 1); + temp.rowname("LB_y(" + std::to_string(i) + "," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(0); + B_Y[k].emplace_back(temp); + + // bound slopes of decision rules g(i,f) to be non-negative + if (data.affine && data.special && SPECIAL_RESTRICTIVE_LDR) for (int f = 1; f <= ((data.N - 1)/3); f++) { + temp.clear(); + temp.addTermX(getVarIndex_2(k, "g", i, f), 1); + temp.rowname("LB_g(" + std::to_string(i) + "," + std::to_string(f) + "," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(0); + B_Y[k].emplace_back(temp); + } + } + } + + + ////////// + // C_XY // + ////////// + for (unsigned int k = C_XY.size(); k <= l; k++) { + if (C_XY.size() < k + 1) C_XY.resize(k + 1); + + C_XY[k].clear(); + + // objective function + temp.clear(); + temp.rowname("OBJ_CONSTRAINT(" + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(0); + temp.addTermX(getVarIndex_1("O", 0), 1); + temp.addTermX(getVarIndex_2(k, "y", data.N), -1); + C_XY[k].emplace_back(temp); + + + // precedence relationships + if (!data.special && !data.affine && DO_MAX) for (int i = 2; i < data.N; i++) for (int j : data.successor[i]) { + temp.clear(); + temp.rowname("PRE(" + std::to_string(i) + "," + std::to_string(j) + "," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(data.duration[i]); + temp.addTermX(getVarIndex_2(k, "y", j), 1); + temp.addTermX(getVarIndex_2(k, "y", i), -1); + temp.addTermX(getVarIndex_1("x", i), data.duration[i]); + C_XY[k].emplace_back(temp); + } + } + + + /////////// + // C_XYQ // + /////////// + for (unsigned int k = C_XYQ.size(); k <= l; k++) { + if (C_XYQ.size() < k + 1) C_XYQ.resize(k + 1); + + C_XYQ[k].clear(); + + + // precedence relationships + for (int i = 2; i < data.N; i++) for (int j : data.successor[i]) { + temp.clear(); + temp.rowname("PRE(" + std::to_string(i) + "," + std::to_string(j) + "," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(data.special ? 0 : data.duration[i]); + temp.addTermX(getVarIndex_2(k, "y", j), 1); + if (data.affine) if (j < data.N) for (int f = 1, flim = (data.special ? (1 + ((data.N - 1)/3)) : data.phi[i].size()); f < flim; f++) { + temp.addTermProduct(getVarIndex_2(k, "g", j, f), U.getParamIndex('q', f), 1); + } + temp.addTermX(getVarIndex_2(k, "y", i), -1); + if (data.affine) for (int f = 1, flim = (data.special ? (1 + ((data.N - 1)/3)) : data.phi[i].size()); f < flim; f++) { + temp.addTermProduct(getVarIndex_2(k, "g", i, f), U.getParamIndex('q', f), -1); + } + if (!data.special) { + temp.addTermX(getVarIndex_1("x", i), data.duration[i]); + for (unsigned int f = 1; f < data.phi[i].size(); f++) { + double coef = data.phi[i][f] * data.duration[i] * FACTOR_ALPHA; + + if (coef != 0.0) { + temp.addTermProduct(getVarIndex_1("x", i), U.getParamIndex('q', f), coef); + temp.addTermQ(U.getParamIndex('q', f), -coef); + } + } + } + else { + if (data.duration[i] > 0.5) { + int l = (i + 1)/3; + temp.addTermQ(U.getParamIndex('q', l), -1); + } + else if (data.duration[i] < -0.5) { + int l = i / 3; + temp.RHS(1); + temp.addTermQ(U.getParamIndex('q', l), 1); + } + else if (!data.affine) { + C_XY[k].emplace_back(temp); + continue; + } + } + C_XYQ[k].emplace_back(temp); + } + + + if (!data.special && data.affine && DO_MAX) for (int i = 2; i < data.N; i++) for (int j : data.successor[i]) { + temp.clear(); + temp.rowname("PRE(" + std::to_string(i) + "," + std::to_string(j) + "," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(data.duration[i]); + temp.addTermX(getVarIndex_2(k, "y", j), 1); + if (j < data.N) for (unsigned int f = 1; f < data.phi[i].size(); f++) { + temp.addTermProduct(getVarIndex_2(k, "g", j, f), U.getParamIndex('q', f), 1); + } + temp.addTermX(getVarIndex_2(k, "y", i), -1); + for (unsigned int f = 1; f < data.phi[i].size(); f++) { + temp.addTermProduct(getVarIndex_2(k, "g", i, f), U.getParamIndex('q', f), -1); + } + temp.addTermX(getVarIndex_1("x", i), data.duration[i]); + C_XYQ[k].emplace_back(temp); + } + + + + // non-negative start times + if (data.affine) if (data.special ? (SPECIAL_RESTRICTIVE_LDR ? 0 : 1) : 1) for (int i = 2; i < data.N; i++) { + temp.clear(); + temp.rowname("GREATER(" + std::to_string(i) + "," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(0); + temp.addTermX(getVarIndex_2(k, "y", i), 1); + for (int f = 1, flim = (data.special ? (1 + ((data.N - 1)/3)) : data.phi[i].size()); f < flim; f++) { + temp.addTermProduct(getVarIndex_2(k, "g", i, f), U.getParamIndex('q', f), 1); + } + C_XYQ[k].emplace_back(temp); + } + } +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_PSP::setInstance(const PSP& d) { + data = d; + hasInteger = 0; + objectiveUnc = 0; + existsFirstStage = !data.special; + numFirstStage = 1 + (data.special ? 0 : data.N); + numSecondStage = data.N + (data.affine * (data.N - 2) * (data.special ? ((data.N - 1)/3) : ((int)data.phi[0].size() - 1))); + numPolicies = 1; + solfilename = data.solfilename; + + makeUncSet(); + makeVars(); + makeConsX(); + makeConsY(0); + assert(isConsistentWithDesign()); +} + +//----------------------------------------------------------------------------------- \ No newline at end of file diff --git a/src/problemInfo_spp.cpp b/src/problemInfo_spp.cpp new file mode 100644 index 0000000..ead72ab --- /dev/null +++ b/src/problemInfo_spp.cpp @@ -0,0 +1,210 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#include "problemInfo_spp.hpp" +#include + +//----------------------------------------------------------------------------------- + +#define SPP_UPPER_BOUND_BUDGET 3 + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_SPP::makeUncSet() { + U.clear(); + + // define uncertain costs + for (int i = 1; i <= data.A; ++i) { + U.addParam(0, 0, 1); + } + + // define budget constraint + std::vector > constraint; + for (int i = 1; i <= data.A; ++i) { + constraint.emplace_back(i, 1); + } + U.addFacet(constraint, 'L', SPP_UPPER_BOUND_BUDGET); +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_SPP::makeVars() { + X.clear(); + Y.clear(); + + // objective variable (considered to be 1st-stage) + X.addVarType("O", 'C', -CPX_INFBOUND, +CPX_INFBOUND, 1); + + // y(i,j) : travel along arc (i,j) + Y.addVarType("y", 'B', 0, 1, 1 + data.N, 1 + data.N); + + // All arrays/atrices are 1-indexed + Y.setUndefinedVar("y", 0, 0); + for (int i = 1; i <= data.N; ++i) { + Y.setUndefinedVar("y", 0, i); + Y.setUndefinedVar("y", i, 0); + Y.setUndefinedVar("y", i, i); + } + + // delete non-existent arcs + for (int i = 1; i <= data.N; ++i) for (int j = 1; j <= data.N; ++j) if (i != j) { + assert(data.costs[i][j] == -1.0 || data.costs[i][j] >= 0); + if (data.costs[i][j] < 0) { + Y.setUndefinedVar("y", i, j); + } + } + assert(Y.getTotalDefVarSize() == data.A); + + // define objective variable + X.setVarObjCoeff(1.0, "O", 0); +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_SPP::makeConsX() { + + + ///////// + // B_X // + ///////// + B_X.clear(); + + + + ///////// + // C_X // + ///////// + C_X.clear(); + + + + ////////// + // C_XQ // + ////////// + C_XQ.clear(); + + // nothing to do +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_SPP::makeConsY(unsigned int l) { + assert(C_XY.size() == B_Y.size()); + assert(C_XY.size() == C_XYQ.size()); + assert(numPolicies >= l); + if (l == 0) { + B_Y.clear(); + C_XY.clear(); + C_XYQ.clear(); + } + ConstraintExpression temp; + + ///////// + // B_Y // + ///////// + for (unsigned int k = B_Y.size(); k <= l; k++) { + if (B_Y.size() < k + 1) B_Y.resize(k + 1); + + B_Y[k].clear(); + + // bounds on y(i,j) + for (int i = 1; i <= data.N; ++i) for (int j = 1; j <= data.N; ++j) { + if (data.costs[i][j] > 0.0) { + temp.clear(); + temp.addTermX(getVarIndex_2(k, "y", i, j), 1); + + temp.rowname("LB_y(" + std::to_string(i) + "," + std::to_string(j) + "," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(0); + B_Y[k].emplace_back(temp); + + temp.rowname("UB_y(" + std::to_string(i) + "," + std::to_string(j) + "," + std::to_string(k) + ")"); + temp.sign('L'); + temp.RHS(1); + B_Y[k].emplace_back(temp); + } + } + } + + + ////////// + // C_XY // + ////////// + for (unsigned int k = C_XY.size(); k <= l; k++) { + if (C_XY.size() < k + 1) C_XY.resize(k + 1); + + C_XY[k].clear(); + + // Flow balances + for (int j = 1; j <= data.N; ++j) { + temp.clear(); + temp.rowname("FLOW(" + std::to_string(j) + "," + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS((j == data.src) ? 1.0 : (j == data.tgt ? -1.0 : 0)); + for (int m = 1; m <= data.N; ++m) if (data.costs[j][m] >= -0.1) { + temp.addTermX(getVarIndex_2(k, "y", j, m), 1.0); + } + for (int i = 1; i <= data.N; ++i) if (data.costs[i][j] >= -0.1) { + temp.addTermX(getVarIndex_2(k, "y", i, j), -1.0); + } + C_XY[k].emplace_back(temp); + } + } + + + /////////// + // C_XYQ // + /////////// + for (unsigned int k = C_XYQ.size(); k <= l; k++) { + if (C_XYQ.size() < k + 1) C_XYQ.resize(k + 1); + + C_XYQ[k].clear(); + + // objective function + temp.clear(); + temp.rowname("OBJ_CONSTRAINT(" + std::to_string(k) + ")"); + temp.sign('G'); + temp.RHS(0); + temp.addTermX(getVarIndex_1("O", 0), 1); + for (int i = 1; i <= data.N; ++i) for (int j = 1; j <= data.N; ++j) { + if (data.costs[i][j] > 0.0) { + const int qind = getVarIndex_2(0, "y", i, j); + const int yind = getVarIndex_2(k, "y", i, j); + temp.addTermX(yind, -data.costs[i][j]); + temp.addTermProduct(yind, qind, -data.costs[i][j] / 2.0); + } + } + C_XYQ[k].emplace_back(temp); + } +} + +//----------------------------------------------------------------------------------- + +void KAdaptableInfo_SPP::setInstance(const SPP& d) { + data = d; + hasInteger = 1; + objectiveUnc = 1; + existsFirstStage = 0; + numFirstStage = 1; + numSecondStage = data.A; + numPolicies = 1; + solfilename = data.solfilename; + + makeUncSet(); + makeVars(); + makeConsX(); + makeConsY(0); + assert(isConsistentWithDesign()); +} + +//----------------------------------------------------------------------------------- diff --git a/src/robustSolver.cpp b/src/robustSolver.cpp new file mode 100644 index 0000000..09ef778 --- /dev/null +++ b/src/robustSolver.cpp @@ -0,0 +1,3644 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#include "robustSolver.hpp" +#include "Constants.h" +#include +#include +#include +#include +#include +#include +#include +///* +#include +#include +#include +static double get_wall_time(){ + struct timeval time; + if (gettimeofday(&time,NULL)){ + // Handle error + return 0; + } + return (double)time.tv_sec + (double)time.tv_usec * .000001; +} +//*/ + +#define CPX_TOL_INT (1.E-5) +#define CPX_TOL_FEA (1.E-5) +#define CPX_TOL_OPT (1.E-5) +#define OUTPUTLEVEL 2 +#define enterCallback(callback_type) {if(OUTPUTLEVEL >= 3) std::cout << " >>>--->>> ENTER " << #callback_type << " Callback <<<---<<< " << std::endl;} +#define exitCallback(callback_type) {if(OUTPUTLEVEL >= 3) std::cout << " >>>--->>> EXIT " << #callback_type << " Callback <<<---<<< " << std::endl; return(0);} +const double EPS_INFEASIBILITY_Q = 1.E-4; +const double EPS_INFEASIBILITY_X = 1.E-4; +static std::vector Q_TEMP; +static std::vector X_TEMP; +static int LABEL_TEMP; +static std::vector > ZT_VALUES, BOUND_VALUES; + +// static global members +static int NUM_DUMMY_NODES = 0, TX = 0; +static double START_TS = 0; +static int CPXPUBLIC cutCB_solve_SRO_cuttingPlane(CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, int *useraction_p); +static int CPXPUBLIC cutCB_solve_KAdaptability_cuttingPlane(CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, int *useraction_p); +static int CPXPUBLIC incCB_solve_KAdaptability_cuttingPlane(CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, double objval, double *x, int *isfeas_p, int *useraction_p); +static int CPXPUBLIC heurCB_solve_KAdaptability_cuttingPlane (CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, double *objval_p, double *x, int *checkfeas_p, int *useraction_p); +static int CPXPUBLIC branchCB_solve_KAdaptability_cuttingPlane(CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, int type, CPXDIM sos, int nodecnt, CPXDIM bdcnt, const CPXDIM *nodebeg, const CPXDIM *indices, const char *lu, const double *bd, const double *nodeest, int *useraction_p); +static int CPXPUBLIC nodeCB_solve_KAdaptability_cuttingPlane(CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, CPXCNT *nodeindex_p, int *useraction_p); +static void CPXPUBLIC deletenodeCB_solve_KAdaptability_cuttingPlane(CPXCENVptr env, int wherefrom, void *cbhandle, CPXCNT seqnum, void *handle); + +// Algorithmic control options +const bool GET_MAX_VIOL = 1; +const int SEPARATION_STRATEGY = 2; +const bool COLLECT_RESULTS = 1; +const bool BB_IMPLEMENT_LAZY_CON = 0; +const bool BNC_BRANCH_ALL_CONSTR = 1; +const bool BNC_DO_STRONG_BRANCH = 0; +const bool SEPARATE_FROM_SAMPLES = 1; +const bool SEPARATE_ALTERNATE = 0; +const bool SEPARATE_ALTERNATE_AVG= 0; +const int BRANCHING_STRATEGY = 2; +// 1 = always branch as per CPLEX, unless necessary to resort to K-adaptability +// 2 = alternate between CPLEX and K-adaptability branching +// 3 = branch as per K-adaptability until gap <= BNC_GAP_VALUE, then switch to strategy 1 +// 4 = always branch as per K-adaptability whenever possible +// 5 = branch as per K-adaptability branching whenever depth modulo K == 0 + +const double BNC_GAP_VALUE = 10; + + +//----------------------------------------------------------------------------------- + +#ifndef NDEBUG +#define MYERROR(_n) do {\ + std::cerr << __FILE__ << " -- Line " << __LINE__ << std::endl;\ + throw(_n);\ +} while (0) +#else +#define MYERROR(_n) do {throw(_n);} while(0) +#endif + +//----------------------------------------------------------------------------------- + +#define MY_SIZE_X(_K) (pInfo->getNumFirstStage() + ((_K) * pInfo->getNumSecondStage())) + +//----------------------------------------------------------------------------------- + +#define MY_SIZE_Q (1 + pInfo->getNoOfUncertainParameters()) + +#define SUCCESS_STATUS ((pInfo->isContinuous() ? ())) + +//----------------------------------------------------------------------------------- + +static inline void setCPXoptions(CPXENVptr& env) { + CPXXsetintparam(env, CPXPARAM_ScreenOutput, CPX_OFF); + CPXXsetintparam(env, CPXPARAM_Threads, NUM_THREADS); + CPXXsetintparam(env, CPXPARAM_Parallel, CPX_PARALLEL_DETERMINISTIC); + CPXXsetintparam(env, CPXPARAM_ClockType, 2); + CPXXsetdblparam(env, CPXPARAM_TimeLimit, TIME_LIMIT); + CPXXsetdblparam(env, CPXPARAM_MIP_Limits_TreeMemory, MEMORY_LIMIT); + // CPXXsetdblparam(env, CPX_PARAM_EPGAP, 0.0); + // CPXXsetdblparam(env, CPX_PARAM_EPAGAP, 1.E-5); + // CPXXsetdblparam(env, CPX_PARAM_EPINT, CPX_TOL_INT); + // CPXXsetdblparam(env, CPX_PARAM_EPRHS, CPX_TOL_FEA); + // CPXXsetdblparam(env, CPX_PARAM_EPOPT, CPX_TOL_OPT); +} + +//----------------------------------------------------------------------------------- + +static inline void addVariable(CPXCENVptr env, CPXLPptr lp, const char xctype = 'C', const double lb = 0, const double ub = +CPX_INFBOUND, const double obj = 0, const char *cname = "") { + if (CPXXnewcols(env, lp, 1, &obj, &lb, &ub, &xctype, &cname)) { + MYERROR(EXCEPTION_CPXNEWCOLS); + } +} + +//----------------------------------------------------------------------------------- + +static inline void addVariable(CPXCENVptr env, CPXLPptr lp, const char xctype = 'C', const double lb = 0, const double ub = +CPX_INFBOUND, const double obj = 0, const std::string cname = "") { + addVariable(env, lp, xctype, lb, ub, obj, cname.c_str()); +} + +static inline void fixVariable(CPXCENVptr env, CPXLPptr lp, const int varIndex, const double val) { + char lu = 'B'; + CPXXtightenbds(env, lp, 1, &varIndex, &lu, &val); +} + +//----------------------------------------------------------------------------------- + +static inline std::pair checkViol(CPXCENVptr env, void *cbdata, int wherefrom, const double rhs, const char sense, const std::vector& cutind, const std::vector& cutval) { + + // get node solution + CPXCLPptr lp = NULL; CPXXgetcallbacklp(env, cbdata, wherefrom, &lp); + CPXDIM numcols = CPXXgetnumcols(env, lp); + std::vector x(numcols); CPXXgetcallbacknodex(env, cbdata, wherefrom, &x[0], 0, numcols - 1); + + // compute left-hand side + double lhs = 0; + for (unsigned int i = 0; i < cutind.size(); i++) { + lhs += x[cutind[i]] * cutval[i]; + } + + // compute violation + double viol = (lhs - rhs) * ((sense == 'G') ? -1.0 : +1.0); + + return std::make_pair(viol, ((std::trunc(1.E+6 * viol) / 1.E+6) >= -1.E-8)); +} + +//----------------------------------------------------------------------------------- +// GENERATE ALL P-TUPLES OF THE INDEX SET {startInd, startInd + 1,...,lastInd} WITH REPETITION +static inline std::vector > generatePTuples(int startInd, int lastInd, const int P) { + assert(lastInd >= startInd); + assert(P >= 1); + + // return value + std::vector > PTuples; + + // Temporary + std::vector Pt(P, startInd); + + + int k = 0; do { + PTuples.emplace_back(Pt); + k = P - 1; + while (k >= 0) { + Pt.at(k)++; + if (Pt.at(k) == 1+lastInd) { + Pt.at(k) = startInd; + k--; + } + else { + break; + } + } + } while (k >= 0); + + // Must have K^P different tuples + assert(static_cast(PTuples.size()) == std::pow(lastInd - startInd + 1, P)); + + return PTuples; +} + +//----------------------------------------------------------------------------------- + +static inline void write(std::ostream& out, std::string N, unsigned int K, std::string seed, std::string status, double final_objval, double total_solution_time, double final_gap) { + out << " N " << " K " << " Seed " << " Status " << " Obj Value " << " Time (sec) " << " Gap (%) \n"; + out << " --- " << " --- " << " ---- " << " -------- " << " --------- " << " ---------- " << " ------- \n"; + out << std::setw(4) << N; + out << std::setw(4) << K; + out << " "; + out << std::setw(4) << std::fixed << seed; + out << " "; + out << std::setw(8) << status; + out << " "; + out << std::setw(8) << std::setprecision(4) << std::fixed << final_objval; + out << " "; + out << std::setw(7) << std::setprecision(2) << std::fixed << total_solution_time; + out << " "; + out << std::setw(5) << std::setprecision(2) << std::fixed << final_gap; + out << " \n\n"; +} + +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- + +KAdaptableSolver::KAdaptableSolver(const KAdaptableInfo& pInfoData) { + pInfo = pInfoData.clone(); +} + +//----------------------------------------------------------------------------------- + +KAdaptableSolver::~KAdaptableSolver() { + if (pInfo) { + delete pInfo; + pInfo = NULL; + } +} + +//----------------------------------------------------------------------------------- + +KAdaptableSolver::KAdaptableSolver (const KAdaptableSolver& S) { + if (S.pInfo) { + pInfo = S.pInfo->clone(); + } + xsol = S.xsol; +} + +//----------------------------------------------------------------------------------- + +KAdaptableSolver& KAdaptableSolver::operator=(const KAdaptableSolver& S) { + if (this == &S) { + return *this; + } + if (pInfo) { + delete pInfo; + pInfo = NULL; + } + if (S.pInfo) { + pInfo = S.pInfo->clone(); + } + xsol = S.xsol; + return *this; +} + +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::setInfo(const KAdaptableInfo& pInfoData) { + if (pInfo) { + delete pInfo; + pInfo = NULL; + } + pInfo = pInfoData.clone(); +} + +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::reset(const KAdaptableInfo& pInfoData, const unsigned int K) { + setInfo(pInfoData); + if (K >= 1) pInfo->resize(K); + xsol.clear(); +} + +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::setX(const std::vector& x, const unsigned int K) { + std::vector q; + if (feasible_KAdaptability(x, K, q)) { + if (xsol.empty()) { + xsol = x; + } + else if (x[0] < xsol[0]) { + xsol = x; + } + } + + return; +} + +//----------------------------------------------------------------------------------- + +unsigned int KAdaptableSolver::getNumPolicies(const std::vector& x) const { + const int num_second_stage_cols = (int)x.size() - pInfo->getNumFirstStage(); + if (num_second_stage_cols < 0) return 0; + if (num_second_stage_cols % pInfo->getNumSecondStage()) return 0; + return static_cast(num_second_stage_cols / pInfo->getNumSecondStage()); +} + +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::resizeX(std::vector& x, const unsigned int K) const { + assert(pInfo); + + const unsigned int Kx = getNumPolicies(x); + + // vectors must be of appropriate size + if (Kx < 1) MYERROR(EXCEPTION_X); + if (x.size() != MY_SIZE_X(Kx)) MYERROR(EXCEPTION_X); + + // remove extra policies + if (Kx > K) { + x.resize(MY_SIZE_X(K)); + } + + // duplicate the last policy + if (Kx < K) { + const auto x2 = x.begin() + pInfo->getNumFirstStage() + ((Kx - 1) * pInfo->getNumSecondStage()); + const std::vector xk(x2, x2 + pInfo->getNumSecondStage()); + for (unsigned int k = Kx; k < K; k++) { + x.insert(x.end(), xk.begin(), xk.end()); + } + } + + assert(x.size() == MY_SIZE_X(K)); + assert(getNumPolicies(x) == K); +} + +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::removeXPolicy(std::vector& x, const unsigned int K, const unsigned int k) const { + assert(pInfo); + + // vectors must be of appropriate size + if (K < 1 || k >= K) MYERROR(EXCEPTION_X); + if (x.size() != MY_SIZE_X(K)) MYERROR(EXCEPTION_X); + + const auto x2 = x.begin() + pInfo->getNumFirstStage() + (k * pInfo->getNumSecondStage()); + x.erase(x2, x2 + pInfo->getNumSecondStage()); + + assert(x.size() == MY_SIZE_X(K - 1)); +} + +//----------------------------------------------------------------------------------- + +const std::vector KAdaptableSolver::getXPolicy(const std::vector& x, const unsigned int K, const unsigned int k) const { + assert(pInfo); + + // vectors must be of appropriate size + if (K < 1 || k >= K) MYERROR(EXCEPTION_K); + if (x.size() != MY_SIZE_X(K)) MYERROR(EXCEPTION_X); + + const auto x2 = x.begin() + pInfo->getNumFirstStage(); + + // 1st-stage solution of policy # k + std::vector xk(x.begin(), x2); + + // 2nd-stage solution of policy # k + xk.insert(xk.end(), x2 + (k * pInfo->getNumSecondStage()), x2 + ((k + 1) * pInfo->getNumSecondStage())); + + assert((int)xk.size() == MY_SIZE_X(1)); + + return xk; +} + +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::updateX(CPXCENVptr env, CPXLPptr lp) const { + assert(pInfo); + + // 1st-stage variables and constraints + auto& X = pInfo->getVarsX(); + auto& C_X = pInfo->getConstraintsX(); + + // define variables + for (int i = 0; i < X.getTotalVarSize(); i++) if (!X.isUndefVar(i)) { + std::string cname; + int ind1, ind2, ind3, ind4, ind5; + X.getVarInfo(i, cname, ind1, ind2, ind3, ind4, ind5); + if (ind1 > -1) { + cname += "(" + std::to_string(ind1); + if (ind2 > -1) cname += "," + std::to_string(ind2); + if (ind3 > -1) cname += "," + std::to_string(ind3); + if (ind4 > -1) cname += "," + std::to_string(ind4); + if (ind5 > -1) cname += "," + std::to_string(ind5); + cname += ")"; + } + addVariable(env, lp, X.getVarColType(i), X.getVarLB(i), X.getVarUB(i), X.getVarObjCoeff(i), cname); + } + + // define constraints + for (const auto& con: C_X) { + con.addToCplex(env, lp); + } +} + +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::updateXQ(CPXCENVptr env, CPXLPptr lp, const std::vector& q, const bool lazy) const { + assert(pInfo); + + // reformulate constraints? + const bool reformulate = q.empty(); + if (!reformulate) if ((int)q.size() != MY_SIZE_Q) MYERROR(EXCEPTION_Q); + + // 1st-stage constraints involving uncertain parameters + auto& C_XQ = pInfo->getConstraintsXQ(); + + // define constraints + if (lazy) { + CPXDIM rcnt = 0; + CPXNNZ nzcnt = 0; + std::vector rhs; + std::vector sense; + std::vector rmatbeg; + std::vector rmatind; + std::vector rmatval; + getXQ_fixedQ(q, rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval); + + CPXXaddlazyconstraints(env, lp, rcnt, nzcnt, &rhs[0], &sense[0], &rmatbeg[0], &rmatind[0], &rmatval[0], NULL); + } + else for (const auto& con: C_XQ) { + con.addToCplex(env, lp, &pInfo->getUncSet(), reformulate, q); + } +} + +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::updateY(CPXCENVptr env, CPXLPptr lp, const unsigned int k) const { + assert(pInfo); + + // vectors must be of appropriate size + if (k >= pInfo->getConstraintsXY().size()) MYERROR(EXCEPTION_K); + + // 1st-stage variables and constraints + auto& Y = pInfo->getVarsY(); + auto& C_XY = pInfo->getConstraintsXY()[k]; + + // define variables + for (int i = 0; i < Y.getTotalVarSize(); i++) if (!Y.isUndefVar(i)) { + std::string cname; + int ind1, ind2, ind3, ind4, ind5; + Y.getVarInfo(i, cname, ind1, ind2, ind3, ind4, ind5); + cname += "_" + std::to_string(k); + if (ind1 > -1) { + cname += "(" + std::to_string(ind1); + if (ind2 > -1) cname += "," + std::to_string(ind2); + if (ind3 > -1) cname += "," + std::to_string(ind3); + if (ind4 > -1) cname += "," + std::to_string(ind4); + if (ind5 > -1) cname += "," + std::to_string(ind5); + cname += ")"; + } + addVariable(env, lp, Y.getVarColType(i), Y.getVarLB(i), Y.getVarUB(i), Y.getVarObjCoeff(i), cname); + } + + // define constraints + for (const auto& con: C_XY) { + con.addToCplex(env, lp); + } +} + +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::updateYQ(CPXCENVptr env, CPXLPptr lp, const unsigned int k, const std::vector& q, const bool lazy) const { + assert(pInfo); + + // vectors must be of appropriate size + if (k >= pInfo->getConstraintsXYQ().size()) MYERROR(EXCEPTION_K); + + // reformulate constraints? + const bool reformulate = q.empty(); + if (!reformulate) if ((int)q.size() != MY_SIZE_Q) MYERROR(EXCEPTION_Q); + + // 1st-stage constraints involving uncertain parameters + auto& C_XYQ = pInfo->getConstraintsXYQ()[k]; + + // define constraints + if (lazy) { + CPXDIM rcnt = 0; + CPXNNZ nzcnt = 0; + std::vector rhs; + std::vector sense; + std::vector rmatbeg; + std::vector rmatind; + std::vector rmatval; + getYQ_fixedQ(k, q, rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval); + + CPXXaddlazyconstraints(env, lp, rcnt, nzcnt, &rhs[0], &sense[0], &rmatbeg[0], &rmatind[0], &rmatval[0], NULL); + } + else for (const auto& con: C_XYQ) { + con.addToCplex(env, lp, &pInfo->getUncSet(), reformulate, q); + } +} + +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::getXQ_fixedQ(const std::vector& q, CPXDIM& rcnt, CPXNNZ& nzcnt, std::vector& rhs, std::vector& sense, std::vector& rmatbeg, std::vector& rmatind, std::vector& rmatval) const { + assert(pInfo); + + // vectors must be of appropriate size + if ((int)q.size() != MY_SIZE_Q) MYERROR(EXCEPTION_Q); + + // initialize members + auto& C_XQ = pInfo->getConstraintsXQ(); + rcnt = 0; + nzcnt = 0; + rhs.clear(); + sense.clear(); + rmatbeg.clear(); + rmatind.clear(); + rmatval.clear(); + + // Statically add constraints involving uncertain parameters + for (const auto& con : C_XQ) { + + rcnt ++; + CPXNNZ nzcnt_p = 0; + rhs.emplace_back(0); + sense.emplace_back('L'); + rmatbeg.emplace_back(rmatind.size()); + std::vector indices; + std::vector values; + + con.getDeterministicConstraint(q, nzcnt_p, rhs.back(), sense.back(), indices, values); + + nzcnt += nzcnt_p; + rmatind.insert(rmatind.end(), indices.begin(), indices.end()); + rmatval.insert(rmatval.end(), values.begin(), values.end()); + } + + assert(rhs.size() == rmatbeg.size()); + assert(rhs.size() == sense.size()); + assert((CPXDIM)rhs.size() == rcnt); + assert(rmatind.size() == rmatval.size()); + assert((CPXNNZ)rmatind.size() == nzcnt); +} + +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::getXQ_fixedX(const std::vector& x, CPXDIM& rcnt, CPXNNZ& nzcnt, std::vector& rhs, std::vector& sense, std::vector& rmatbeg, std::vector& rmatind, std::vector& rmatval) const { + assert(pInfo); + + // vectors must be of appropriate size + if ((int)x.size() < MY_SIZE_X(1)) MYERROR(EXCEPTION_X); + + // initialize members + auto& C_XQ = pInfo->getConstraintsXQ(); + rcnt = 0; + nzcnt = 0; + rhs.clear(); + sense.clear(); + rmatbeg.clear(); + rmatind.clear(); + rmatval.clear(); + + // Statically add constraints involving uncertain parameters + for (const auto& con : C_XQ) { + + rcnt ++; + CPXNNZ nzcnt_p = 0; + rhs.emplace_back(0); + sense.emplace_back('L'); + rmatbeg.emplace_back(rmatind.size()); + std::vector indices; + std::vector values; + + con.getStochasticConstraint(x, nzcnt_p, rhs.back(), sense.back(), indices, values); + + nzcnt += nzcnt_p; + rmatind.insert(rmatind.end(), indices.begin(), indices.end()); + rmatval.insert(rmatval.end(), values.begin(), values.end()); + } + + assert(rhs.size() == rmatbeg.size()); + assert(rhs.size() == sense.size()); + assert((CPXDIM)rhs.size() == rcnt); + assert(rmatind.size() == rmatval.size()); + assert((CPXNNZ)rmatind.size() == nzcnt); +} + +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::getYQ_fixedQ(const unsigned int k, const std::vector& q, CPXDIM& rcnt, CPXNNZ& nzcnt, std::vector& rhs, std::vector& sense, std::vector& rmatbeg, std::vector& rmatind, std::vector& rmatval) const { + assert(pInfo); + + // vectors must be of appropriate size + if (k >= pInfo->getConstraintsXYQ().size()) MYERROR(EXCEPTION_K); + if ((int)q.size() != MY_SIZE_Q) MYERROR(EXCEPTION_Q); + + // initialize members + auto& C_XYQ = pInfo->getConstraintsXYQ()[k]; + rcnt = 0; + nzcnt = 0; + rhs.clear(); + sense.clear(); + rmatbeg.clear(); + rmatind.clear(); + rmatval.clear(); + + // Statically add constraints involving uncertain parameters + for (const auto& con : C_XYQ) { + + rcnt ++; + CPXNNZ nzcnt_p = 0; + rhs.emplace_back(0); + sense.emplace_back('L'); + rmatbeg.emplace_back(rmatind.size()); + std::vector indices; + std::vector values; + + con.getDeterministicConstraint(q, nzcnt_p, rhs.back(), sense.back(), indices, values); + + nzcnt += nzcnt_p; + rmatind.insert(rmatind.end(), indices.begin(), indices.end()); + rmatval.insert(rmatval.end(), values.begin(), values.end()); + } + + assert(rhs.size() == rmatbeg.size()); + assert(rhs.size() == sense.size()); + assert((CPXDIM)rhs.size() == rcnt); + assert(rmatind.size() == rmatval.size()); + assert((CPXNNZ)rmatind.size() == nzcnt); +} + +//----------------------------------------------------------------------------------- + +void KAdaptableSolver::getYQ_fixedX(const unsigned int k, const std::vector& x, CPXDIM& rcnt, CPXNNZ& nzcnt, std::vector& rhs, std::vector& sense, std::vector& rmatbeg, std::vector& rmatind, std::vector& rmatval) const { + assert(pInfo); + + // vectors must be of appropriate size + if (k >= pInfo->getConstraintsXYQ().size()) MYERROR(EXCEPTION_K); + if (x.size() < MY_SIZE_X(k+1)) MYERROR(EXCEPTION_X); + + // initialize members + auto& C_XYQ = pInfo->getConstraintsXYQ()[k]; + rcnt = 0; + nzcnt = 0; + rhs.clear(); + sense.clear(); + rmatbeg.clear(); + rmatind.clear(); + rmatval.clear(); + + // Statically add constraints involving uncertain parameters + for (const auto& con : C_XYQ) { + + rcnt ++; + CPXNNZ nzcnt_p = 0; + rhs.emplace_back(0); + sense.emplace_back('L'); + rmatbeg.emplace_back(rmatind.size()); + std::vector indices; + std::vector values; + + con.getStochasticConstraint(x, nzcnt_p, rhs.back(), sense.back(), indices, values); + + nzcnt += nzcnt_p; + rmatind.insert(rmatind.end(), indices.begin(), indices.end()); + rmatval.insert(rmatval.end(), values.begin(), values.end()); + } + + assert(rhs.size() == rmatbeg.size()); + assert(rhs.size() == sense.size()); + assert((CPXDIM)rhs.size() == rcnt); + assert(rmatind.size() == rmatval.size()); + assert((CPXNNZ)rmatind.size() == nzcnt); +} + +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- + +double KAdaptableSolver::getLowerBound() const { + assert(pInfo); + + + double LB, LB1; + std::vector xtemp; + + int status; + CPXENVptr env; + CPXLPptr lp; + + std::vector indices; + std::vector xctype; + + + + // 1> You can always use any deterministic solution as a lower bound + // 2> In case of objective uncertainty, you can do better + + + // 1> Get deterministic objective value for nominal q + if (solve_DET(pInfo->getNominal(), xtemp)) { + // did not terminate successfully + LB = -std::numeric_limits::max(); + } + else { + LB = xtemp[0]; + } + + + // 2> Objective uncertainty only + // min_x max_q min_y + if (pInfo->hasObjectiveUncOnly()) { + // (1) Make y continuous (lower bound) + // (2) Flip min_y with max_q + // This is valid because of equality in max-min inequality + // See, for e.g., the minimax theorem + // (3) You get classical (static) RO with x as before and y continuous + + + // initialize CPLEX objects + env = NULL; + lp = NULL; + env = CPXXopenCPLEX (&status); + lp = CPXXcreateprob(env, &status, "LowerBoundingProblem"); + + // set options + setCPXoptions(env); + CPXXchgprobtype(env, lp, (pInfo->isContinuous() ? CPXPROB_LP : CPXPROB_MILP)); + CPXXchgobjsen(env, lp, CPX_MIN); + + // define variables and constraints + updateX(env, lp); + updateXQ(env, lp); + updateY(env, lp, 0); + updateYQ(env, lp, 0); + + + // change y variables to continuous + indices.resize(pInfo->getNumSecondStage()); + xctype.resize(indices.size(), CPX_CONTINUOUS); + std::iota(indices.begin(), indices.end(), static_cast(pInfo->getNumFirstStage())); + CPXXchgctype(env, lp, indices.size(), &indices[0], &xctype[0]); + + + // solve problem + status = CPXXmipopt(env, lp); + if (!status) { + // Update lower bound + CPXXgetbestobjval(env, lp, &LB1); + if (LB < LB1) { + LB = LB1; + } + } + else { + MYERROR(status); + } + + // Free memory + CPXXfreeprob(env, &lp); + CPXXcloseCPLEX (&env); + } + + + return LB; +} + +//----------------------------------------------------------------------------------- + +double KAdaptableSolver::getLowerBound_2(std::vector& q) const { + assert(pInfo->hasObjectiveUncOnly()); + assert(pInfo->isSecondStageOnly()); + assert(pInfo->getConstraintsXYQ()[0].size() == 1); + + int status; + CPXENVptr env = NULL; + CPXLPptr lp = NULL; + + env = pInfo->getUncSet().getENVObject(); + lp = pInfo->getUncSet().getLPObject(&status); + + int ind = 0; double coef = 1; char lb = 'L'; char ub = 'U'; double lbVal = -CPX_INFBOUND; double ubVal = +CPX_INFBOUND; + + + CPXXchgobjsen(env, lp, CPX_MAX); + CPXXchgprobtype(env, lp, CPXPROB_LP); + CPXXchgobj(env, lp, 1, &ind, &coef); + for (int l = 1; l <= pInfo->getUncSet().getNoOfUncertainParameters(); ++l) { + coef = 0; + CPXXchgobj(env, lp, 1, &l, &coef); + } + CPXXchgbds(env, lp, 1, &ind, &lb, &lbVal); + CPXXchgbds(env, lp, 1, &ind, &ub, &ubVal); + + UncertaintySet MyU; + + // 1st-stage variables and constraints + auto& Y = pInfo->getVarsY(); + auto& C_XY = pInfo->getConstraintsXY()[0]; + + // define variables + for (int i = 0; i < Y.getTotalVarSize(); i++) if (!Y.isUndefVar(i)) { + MyU.addParam(Y.getVarLB(i), Y.getVarLB(i), Y.getVarUB(i)); + } + + + // define constraints + for (const auto& con: C_XY) { + std::vector qtemp; + CPXNNZ nzcnt; + double rhs; + char sense; + std::vector rmatind; + std::vector rmatval; + con.getDeterministicConstraint(qtemp, nzcnt, rhs, sense, rmatind, rmatval); + std::vector > facet; + for (CPXNNZ i = 0; i < nzcnt; i++) { + facet.emplace_back(rmatind[i], rmatval[i]); + } + MyU.addFacet(facet, sense, rhs); + } + + // define constraint to be dualized + const ConstraintExpression myexp = pInfo->getConstraintsXYQ()[0][0].flip0(); + myexp.addToCplex(env, lp, &MyU, true); + + CPXXchgobjsen(env, lp, CPX_MAX); + CPXXchgprobtype(env, lp, CPXPROB_LP); + + // solve problem + status = CPXXlpopt(env, lp); + if (status) MYERROR(status); + + // get status + assert(CPXXgetstat(env, lp) == CPX_STAT_OPTIMAL); + + // Get primal solution vector + q.resize(1 + pInfo->getUncSet().getNoOfUncertainParameters(), 100); + CPXXgetx(env, lp, &q[0], 0, q.size() - 1); + + // Free memory + if (lp) CPXXfreeprob(env, &lp); + + + + return q[0]; +} + +//----------------------------------------------------------------------------------- + +double KAdaptableSolver::getWorstCase(const std::vector& x, const unsigned int K, std::vector& q) const { + assert(pInfo); + + // must be of appropriate size + if (K < 1 || (int)K > pInfo->getNumPolicies()) MYERROR(EXCEPTION_K); + if (!feasible_DET_K(x, K, X_TEMP)) MYERROR(EXCEPTION_X); + + // if infeasible, then any violating scenario is okay + if (!feasible_KAdaptability(x, K, q)) return +std::numeric_limits::max(); + + + CPXDIM rcnt = 0; + CPXNNZ nzcnt = 0; + int objective_idx = -1; + std::vector rhs; + std::vector sense; + std::vector rmatbeg; + std::vector rmatind; + std::vector rmatval; + + + ///////////////////////////////////////////////////////////////////////// + // NOTE: I'M ASSUMING THAT THE ORIGINAL PROBLEM IS ALWAYS REFORMULATED // + // USING AN EPIGRAPHICAL VARIABLE FOR THE OBJECTIVE FUNCTION // + ///////////////////////////////////////////////////////////////////////// + getYQ_fixedQ(0, pInfo->getNominal(), rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval); + for (CPXDIM j = 0; j < rcnt; j++) { + auto begin = rmatind.begin() + rmatbeg[j]; + auto end = rmatind.begin() + ((j+1 == rcnt) ? nzcnt : rmatbeg[j+1]); + if (std::find(begin, end, 0) != end) { + objective_idx = j; + break; + } + } + + + + // 1> In case of pure objective uncertainty, solve separation problem + // 2> In case of pure constraint uncertainty, get any scenario and return deterministic objective + // 3> In case of both objective and constraint uncertainty, need more work + + + // 1> Separation problem already solved + if (pInfo->hasObjectiveUncOnly()) return x[0] + q[0]; + + + // 2> Any scenario is okay + if (objective_idx < 0) return x[0]; + + + // 3> Solve (approximately) sup_q inf_k [objective(q, k) : such that policy k is feasible for q] + int status; + CPXENVptr env = NULL; + CPXLPptr lp = NULL; + char b; + double bd; + int index; + auto& U = pInfo->getUncSet(); + const auto PTuples = generatePTuples(0, 1, K); + std::vector max_q(1, -std::numeric_limits::max()); + bool single_constraint = (rcnt == 2); + + + // tuple[k] denotes whether policy k is feasible in the current partition + for (const auto& tuple : PTuples) { + // ignore the case where all policies are infeasible + int sum = 0; for (auto t : tuple) sum += t; + if (sum < 1) continue; + + // get solver objects from uncertainty set + env = U.getENVObject(); + lp = U.getLPObject(&status); + + // convert problem to maximization MILP with objective = [tau] + CPXXchgobjsen(env, lp, CPX_MAX); + CPXXchgprobtype(env, lp, single_constraint ? CPXPROB_LP : CPXPROB_MILP); + for (index = 0; index <= U.getNoOfUncertainParameters(); ++index) { + double coef = static_cast((index == 0)); + CPXXchgobj(env, lp, 1, &index, &coef); + } + + // Change bounds on objective function variable to make it (-Inf, +Inf) + index = 0; b = 'L'; bd = -CPX_INFBOUND; CPXXchgbds(env, lp, 1, &index, &b, &bd); + index = 0; b = 'U'; bd = +CPX_INFBOUND; CPXXchgbds(env, lp, 1, &index, &b, &bd); + index = CPXXgetnumcols(env, lp) - 1; + + // Add constraints for all policies + for (unsigned int k = 0; k < K; ++k) { + + // get all uncertain constraints in policy k for fixed x + getYQ_fixedX(k, x, rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval); + + // only for infeasible policies: sum(j, z_jk) = 1 + ConstraintExpression zConstraint("z(" + std::to_string(k) + ")"); + zConstraint.sign('E'); + zConstraint.RHS(1); + + for (CPXDIM j = 0; j < rcnt; j++) { + assert(sense[j] == 'G' || sense[j] == 'L'); + CPXNNZ ilim = (j+1 == rcnt) ? nzcnt : rmatbeg[j+1]; + double linrhs = rhs[j]; + char linsense = sense[j]; + std::vector linind(rmatind.begin() + rmatbeg[j], rmatind.begin() + ilim); + std::vector linval(rmatval.begin() + rmatbeg[j], rmatval.begin() + ilim); + if (linind.empty()) continue; + + // Policy k must be feasible + if (tuple[k]) { + // IF objective function, maximize violation: + // -tau + [violation expression of constraint j in terms of q for fixed x in policy k] >= 0 + // ELSE + // add as a regular constraint + if (j == objective_idx) { + linsense = (sense[j] == 'G') ? 'L' : 'G'; + linind.emplace_back(0); + linval.emplace_back((sense[j] == 'G') ? 1 : -1); + } + CPXXaddrows(env, lp, 0, 1, linind.size(), &linrhs, &linsense, &rmatbeg[0], &linind[0], &linval[0], NULL, NULL); + } + + // Policy k must be infeasible + else { + if (j == objective_idx) continue; + + // z_jk indicates whether constraint j is violated + // z_jk => [reversed expression of constraint j in terms of q for fixed x in policy k] + linrhs += ((sense[j] == 'G') ? -1.0 : 1.0)*EPS_INFEASIBILITY_X; + linsense = (sense[j] == 'G') ? 'L' : 'G'; + + if (single_constraint) { + CPXXaddrows(env, lp, 0, 1, linind.size(), &linrhs, &linsense, &rmatbeg[0], &linind[0], &linval[0], NULL, NULL); + } + else { + // define z_jk + const int z_index = ++index; + addVariable(env, lp, 'B', 0, 1, 0, "z(" + std::to_string(j) + "," + std::to_string(k) + ")"); + + // add indicator constraint + CPXXaddindconstr(env, lp, z_index, 0, linind.size(), linrhs, linsense, &linind[0], &linval[0], NULL); + + // UPDATE original constraints + zConstraint.addTermX(z_index, 1.0); + } + } + } + + // add constraints only if at least one non-zero was added + if (!tuple[k]) if (!zConstraint.isEmpty()) if (!single_constraint) { + zConstraint.addToCplex(env, lp); + } + } + + + // solve MILP + status = single_constraint ? CPXXlpopt(env, lp) : CPXXmipopt(env, lp); + if (!status) { + status = CPXXgetstat(env, lp); + if (status == CPX_STAT_OPTIMAL || status == CPXMIP_OPTIMAL || status == CPXMIP_OPTIMAL_TOL) { + status = 0; + q.resize(1 + U.getNoOfUncertainParameters()); + CPXXgetx(env, lp, &q[0], 0, U.getNoOfUncertainParameters()); + if (q[0] > max_q[0]) { + max_q = q; + } + } + else if (status == CPX_STAT_INFEASIBLE || status == CPX_STAT_INForUNBD || status == CPXMIP_INFEASIBLE || status == CPXMIP_INForUNBD) { + CPXXfreeprob(env, &lp); + continue; + } + else { + MYERROR(status); + } + } + else { + MYERROR(status); + } + + // Free memory + CPXXfreeprob(env, &lp); + } + + // something went wrong + if (max_q.size() == 1) MYERROR(999); + + q = max_q; + return q[0] + x[0]; +} + +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- + +bool KAdaptableSolver::feasible_DET_K(const std::vector& x, const unsigned int K1, std::vector& xx) const { + assert(pInfo); + + unsigned int K = K1; + unsigned int Kx = getNumPolicies(x); + + // vectors must be of appropriate size + if (K < 1 || (int)K > pInfo->getNumPolicies()) MYERROR(EXCEPTION_K); + if (Kx < 1) MYERROR(EXCEPTION_X); + + // Reduce # of policies to check (if possible) + if (Kx < K) K = Kx; assert(x.size() >= MY_SIZE_X(K)); + + + // Check deterministic constraints + for (const auto& con: pInfo->getConstraintsX()) { + if (getViolation(con, x) > EPS_INFEASIBILITY_X) return 0; + } + + // Check for feasible 2nd-stage policies + std::vector infeas; + for (unsigned int i = 0; i < K; ++i) { + for (const auto& con: pInfo->getConstraintsXY()[i]) { + if (getViolation(con, x) > EPS_INFEASIBILITY_X) { + infeas.emplace_back(i); + break; + } + } + } + + // No 2nd-stage policies are feasible + if (infeas.size() == K) { + return 0; + } + // Discard infeasible 2nd-stage policies; + else { + + // modified x-vector (includes only 'feasible' 2nd-stage policies) + xx = x; + + for (unsigned int i = 0; i < infeas.size(); i++) { + const unsigned int k = infeas[i]; // to be discarded + const auto it = xx.begin() + pInfo->getNumFirstStage() + (k * pInfo->getNumSecondStage()); + + xx.erase(it, it + pInfo->getNumSecondStage()); + for (unsigned int j = i + 1; j < infeas.size(); j++) { + infeas[j]--; + } + K--; + } + } + assert(K); + assert(xx.size() >= MY_SIZE_X(K)); + + return true; +} + +//----------------------------------------------------------------------------------- + +bool KAdaptableSolver::feasible_XQ(const std::vector& x, std::vector& q) const { + assert ((int)x.size() >= pInfo->getNumFirstStage()); + + // Obtain worst violation? + std::vector max_q(1, -std::numeric_limits::max()); + + // Loop through 1st-stage uncertain constraints + for (const auto& con: pInfo->getConstraintsXQ()) { + q = getViolation(con, &pInfo->getUncSet(), x); + if (GET_MAX_VIOL) { + if (q[0] > max_q[0]) { + max_q = q; + } + } + else if (q[0] > EPS_INFEASIBILITY_Q) { + return false; + } + } + + // return maximum violation + if (GET_MAX_VIOL) { + q = max_q; + if (q[0] > EPS_INFEASIBILITY_Q) { + return false; + } + } + + return true; +} + +//----------------------------------------------------------------------------------- + +bool KAdaptableSolver::feasible_XQ(const std::vector& x, const std::vector >& samples, int & label) const { + assert ((int)x.size() >= pInfo->getNumFirstStage()); + + // vectors must be of appropriate size + if ((int)x.size() != MY_SIZE_X(1)) MYERROR(EXCEPTION_X); + for (const auto& q : samples) { + if ((int)q.size() != MY_SIZE_Q) MYERROR(EXCEPTION_Q); + } + + // Obtain worst violation? + double maxViol = -std::numeric_limits::max(); + int label_max = 0; + + // Check each label + for (label = 0; label < (int)samples.size(); label++) { + for (const auto& con: pInfo->getConstraintsXQ()) { + const double viol = getViolation(con, x, samples[label]); + if (GET_MAX_VIOL) { + if (viol > maxViol) { + maxViol = viol; + label_max = label; + } + } + else if (viol > EPS_INFEASIBILITY_Q) { + return false; + } + } + } + + // return maximum violation + if (GET_MAX_VIOL) { + label = label_max; + if (maxViol > EPS_INFEASIBILITY_Q) { + return false; + } + } + + return true; +} + +//----------------------------------------------------------------------------------- + +bool KAdaptableSolver::feasible_YQ(const std::vector& x, const unsigned int K, std::vector& q) const { + assert(K >= 1); + assert(x.size() >= MY_SIZE_X(K)); + + // Obtain worst violation? + std::vector max_q(1, -std::numeric_limits::max()); + + + ////////////////////////////////////////////////////////////////////////////// + // Objective uncertainty only: construct appropriate K-Adaptable expression // + ////////////////////////////////////////////////////////////////////////////// + if (pInfo->hasObjectiveUncOnly()) { + // temporary + std::vector CExpr(K); + + for (unsigned int i = 0; i < K; ++i) { + assert(pInfo->getConstraintsXYQ()[i].size() == 1); + CExpr[i] = pInfo->getConstraintsXYQ()[i][0]; + } + const KAdaptableExpression KExpr (CExpr, "feasible_KAdaptability"); + + // compute violation + q = KExpr.evaluate(&pInfo->getUncSet(), pInfo->getUncSet().getNominal(), x); + if (GET_MAX_VIOL) { + if (q[0] > max_q[0]) { + max_q = q; + } + } + else if (q[0] > EPS_INFEASIBILITY_Q) { + return false; + } + } + + + /////////////////////////////////////// + // Explicitly solve all possible LPs // + /////////////////////////////////////// + else if (SEPARATION_STRATEGY == 0 || K == 1) { + + // temporaries + std::vector CExpr(K); + unsigned int i; + + // get # of 2nd-stage constraints + const auto CI = pInfo->getConstraintsXYQ()[0].size(); + const auto PTuples = generatePTuples(0, CI - 1, K); + + for (const auto& tuple : PTuples) { + assert(tuple.size() == K); + for (i = 0; i < K; ++i) { + if (tuple[i] < (int)pInfo->getConstraintsXYQ()[i].size()) { + CExpr[i] = pInfo->getConstraintsXYQ()[i][tuple[i]]; + } + else { + break; + } + } + assert(i == K); + if (i < K) continue; + + const KAdaptableExpression KExpr(CExpr, "Check"); + + // compute violation + q = KExpr.evaluate(&pInfo->getUncSet(), pInfo->getUncSet().getNominal(), x); + if (GET_MAX_VIOL) { + if (q[0] > max_q[0]) { + max_q = q; + } + } + else if (q[0] > EPS_INFEASIBILITY_Q) { + return false; + } + } + } + + + ///////////////////////////// + // Generic: formulate MILP // + ///////////////////////////// + else if (SEPARATION_STRATEGY == 1) { + + /////////////////////////////////////////// + // TODO: normalize constraint violations // + /////////////////////////////////////////// + + int status; + CPXENVptr env = NULL; + CPXLPptr lp = NULL; + char b; + double bd; + int index; + auto& U = pInfo->getUncSet(); + const auto qLB = U.getLowerBounds(); + const auto qUB = U.getUpperBounds(); + + // get solver objects from uncertainty set + env = U.getENVObject(); + lp = U.getLPObject(&status); + + + + // convert problem to maximization MILP with objective = [tau] + CPXXchgobjsen(env, lp, CPX_MAX); + CPXXchgprobtype(env, lp, CPXPROB_MILP); + for (index = 0; index <= U.getNoOfUncertainParameters(); ++index) { + double coef = static_cast((index == 0)); + CPXXchgobj(env, lp, 1, &index, &coef); + } + + + // Change bounds on objective function variable to make it (-Inf, +Inf) + index = 0; b = 'L'; bd = -CPX_INFBOUND; CPXXchgbds(env, lp, 1, &index, &b, &bd); + index = 0; b = 'U'; bd = +CPX_INFBOUND; CPXXchgbds(env, lp, 1, &index, &b, &bd); + index = CPXXgetnumcols(env, lp) - 1; + + // Add constraints maximizing violation + for (unsigned int k = 0; k < K; ++k) { + CPXDIM rcnt = 0; + CPXNNZ nzcnt = 0; + std::vector rhs; + std::vector sense; + std::vector rmatbeg; + std::vector rmatind; + std::vector rmatval; + + // get all uncertain constraints in policy k for fixed x + getYQ_fixedX(k, x, rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval); + + // -tau + sum(j, z_jk * [violation expression in terms of q for fixed x in policy k]) >= 0 + ConstraintExpression tauConstraint("tau(" + std::to_string(k) + ")"); + tauConstraint.sign('G'); + tauConstraint.RHS(0); + + // sum(j, z_jk) = 1 + ConstraintExpression zConstraint("z(" + std::to_string(k) + ")"); + zConstraint.sign('E'); + zConstraint.RHS(1); + + for (CPXDIM j = 0; j < rcnt; j++) { + assert(sense[j] == 'G' || sense[j] == 'L'); + const double f = (sense[j] == 'G') ? (-1.0) : (+1.0); + + // define z_jk + const std::string z_name = "z(" + std::to_string(j) + "," + std::to_string(k) + ")"; + const int z_index = ++index; + addVariable(env, lp, 'B', 0, 1, 0, z_name); + + for (CPXNNZ i = rmatbeg[j], ilim = ((j+1 == rcnt) ? nzcnt : rmatbeg[j+1]); i < ilim; i++) { + + // lower and upper bounds of q_i + const double qlb = qLB[rmatind[i]]; + const double qub = qUB[rmatind[i]]; + + // define bilinear = z_jk * q_i + const std::string bl_name = "BL(" + z_name + ",q(" + std::to_string(rmatind[i]) + "))"; + const int bl_index = ++index; + addVariable(env, lp, 'C', qlb, qub, 0, bl_name); + + // define bilinear constraints + ConstraintExpression bl; + + // under-estimator 1: bl >= z_jk * qlb + bl.clear(); + bl.sign('G'); + bl.RHS(0); + bl.rowname(bl_name + "_l1"); + bl.addTermX(bl_index, 1.0); + if (qlb != 0.0) bl.addTermX(z_index, -qlb); + bl.addToCplex(env, lp); + + // over-estimator 1: bl <= z_jk * qub + bl.clear(); + bl.sign('L'); + bl.RHS(0); + bl.rowname(bl_name + "_u1"); + bl.addTermX(bl_index, 1.0); + if (qub != 0.0) bl.addTermX(z_index, -qub); + bl.addToCplex(env, lp); + + // under-estimator 2: bl >= (z_jk * qub) + q_i - qub + bl.clear(); + bl.sign('G'); + bl.RHS(-qub); + bl.rowname(bl_name + "_l2"); + bl.addTermX(bl_index, 1.0); + if (qub != 0.0) bl.addTermX(z_index, -qub); + bl.addTermX(rmatind[i], -1.0); + bl.addToCplex(env, lp); + + // over-estimator 1: bl <= (z_jk * qlb) + q_i - qlb + bl.clear(); + bl.sign('L'); + bl.RHS(-qlb); + bl.rowname(bl_name + "_u2"); + bl.addTermX(bl_index, 1.0); + if (qlb != 0.0) bl.addTermX(z_index, -qlb); + bl.addTermX(rmatind[i], -1.0); + bl.addToCplex(env, lp); + + // UPDATE original constraints + if (rmatval[i] != 0.0) tauConstraint.addTermX(bl_index, (f * rmatval[i])); + } + + // UPDATE original constraints + zConstraint.addTermX(z_index, 1.0); + if (rhs[j] != 0.0) tauConstraint.addTermX(z_index, -(f * rhs[j])); + } + + // add constraints only if at least one non-zero was added + if (!tauConstraint.isEmpty()) { + tauConstraint.addTermX(0, -1.0); + tauConstraint.addToCplex(env, lp); + } + if (!zConstraint.isEmpty()) { + zConstraint.addToCplex(env, lp); + } + } + + + // solve MILP + status = CPXXmipopt(env, lp); + if (!status) { + status = CPXXgetstat(env, lp); + if (status == CPXMIP_OPTIMAL || status == CPXMIP_OPTIMAL_TOL) { + status = 0; + q.resize(1 + U.getNoOfUncertainParameters()); + CPXXgetx(env, lp, &q[0], 0, U.getNoOfUncertainParameters()); + if (GET_MAX_VIOL) { + if (q[0] > max_q[0]) { + max_q = q; + } + } + else if (q[0] > EPS_INFEASIBILITY_Q) { + CPXXfreeprob(env, &lp); + return false; + } + } + else { + MYERROR(status); + } + } + else { + MYERROR(status); + } + + // Free memory + CPXXfreeprob(env, &lp); + } + + + ///////////////////////////////////////////////////////// + // Generic: formulate MILP using indicator constraints // + ///////////////////////////////////////////////////////// + else { + + /////////////////////////////////////////// + // TODO: normalize constraint violations // + /////////////////////////////////////////// + + int status; + CPXENVptr env = NULL; + CPXLPptr lp = NULL; + char b; + double bd; + int index; + auto& U = pInfo->getUncSet(); + + // get solver objects from uncertainty set + env = U.getENVObject(); + lp = U.getLPObject(&status); + + + + // convert problem to maximization MILP with objective = [tau] + CPXXchgobjsen(env, lp, CPX_MAX); + CPXXchgprobtype(env, lp, CPXPROB_MILP); + for (index = 0; index <= U.getNoOfUncertainParameters(); ++index) { + double coef = static_cast((index == 0)); + CPXXchgobj(env, lp, 1, &index, &coef); + } + + + // Change bounds on objective function variable to make it (-Inf, +Inf) + index = 0; b = 'L'; bd = -CPX_INFBOUND; CPXXchgbds(env, lp, 1, &index, &b, &bd); + index = 0; b = 'U'; bd = +CPX_INFBOUND; CPXXchgbds(env, lp, 1, &index, &b, &bd); + index = CPXXgetnumcols(env, lp) - 1; + + // Add constraints maximizing violation + for (unsigned int k = 0; k < K; ++k) { + CPXDIM rcnt = 0; + CPXNNZ nzcnt = 0; + std::vector rhs; + std::vector sense; + std::vector rmatbeg; + std::vector rmatind; + std::vector rmatval; + + // get all uncertain constraints in policy k for fixed x + getYQ_fixedX(k, x, rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval); + + // sum(j, z_jk) = 1 + ConstraintExpression zConstraint("z(" + std::to_string(k) + ")"); + zConstraint.sign('E'); + zConstraint.RHS(1); + + for (CPXDIM j = 0; j < rcnt; j++) { + assert(sense[j] == 'G' || sense[j] == 'L'); + + // define linear portion of indicator constraint + // z_jk => -tau + [violation expression of constraint j in terms of q for fixed x in policy k] >= 0 + CPXNNZ ilim = (j+1 == rcnt) ? nzcnt : rmatbeg[j+1]; + double linrhs = rhs[j]; + int linsense = (sense[j] == 'G') ? 'L' : 'G'; + std::vector linind(rmatind.begin() + rmatbeg[j], rmatind.begin() + ilim); + std::vector linval(rmatval.begin() + rmatbeg[j], rmatval.begin() + ilim); + std::string indname = "tau(" + std::to_string(j) + "," + std::to_string(k) + ")"; + + if (!linind.empty()) { + // add tau term + linind.emplace_back(0); + linval.emplace_back((sense[j] == 'G') ? 1 : -1); + + // define z_jk + const int z_index = ++index; + addVariable(env, lp, 'B', 0, 1, 0, "z(" + std::to_string(j) + "," + std::to_string(k) + ")"); + + // add indicator constraint + CPXXaddindconstr(env, lp, z_index, 0, linind.size(), linrhs, linsense, &linind[0], &linval[0], indname.c_str()); + + // UPDATE original constraints + zConstraint.addTermX(z_index, 1.0); + } + } + + // add constraints only if at least one non-zero was added + if (!zConstraint.isEmpty()) { + zConstraint.addToCplex(env, lp); + } + } + + + // solve MILP + status = CPXXmipopt(env, lp); + if (!status) { + status = CPXXgetstat(env, lp); + if (status == CPXMIP_OPTIMAL || status == CPXMIP_OPTIMAL_TOL) { + status = 0; + q.resize(1 + U.getNoOfUncertainParameters()); + CPXXgetx(env, lp, &q[0], 0, U.getNoOfUncertainParameters()); + if (GET_MAX_VIOL) { + if (q[0] > max_q[0]) { + max_q = q; + } + } + else if (q[0] > EPS_INFEASIBILITY_Q) { + CPXXfreeprob(env, &lp); + return false; + } + } + else { + MYERROR(status); + } + } + else { + MYERROR(status); + } + + // Free memory + CPXXfreeprob(env, &lp); + } + + + + + // return maximum violation + if (GET_MAX_VIOL) { + q = max_q; + if (q[0] > EPS_INFEASIBILITY_Q) { + return false; + } + } + + return true; +} + +//----------------------------------------------------------------------------------- + +bool KAdaptableSolver::feasible_YQ(const std::vector& x, const unsigned int K, const std::vector >& samples, int & label, bool heur) const { + assert(K >= 1); + assert(x.size() >= MY_SIZE_X(K)); + const double eps = (heur ? 1.E-2 : EPS_INFEASIBILITY_Q); + + // Obtain worst violation? + double maxViol = -std::numeric_limits::max(); + int label_max = 0; + + // Check each label + for (label = 0; label < (int)samples.size(); label++) { + double labelViol = +std::numeric_limits::max(); + unsigned int k; + for (k = 0; k < K; k++) { + bool isPolicyFeasible = true; + double policyViol = -std::numeric_limits::max(); + for (const auto& con: pInfo->getConstraintsXYQ()[k]) { + const double v = getViolation(con, x, samples[label]); + + if (GET_MAX_VIOL) { + if (v > policyViol) { + policyViol = v; + } + } + else if (v > eps) { + isPolicyFeasible = false; + if (v < labelViol) { + labelViol = v; + } + break; + } + } + + if (GET_MAX_VIOL) { + if (policyViol > eps) { + isPolicyFeasible = false; + if (policyViol < labelViol) { + labelViol = policyViol; + } + } + } + + if (isPolicyFeasible) { + break; + } + } + assert(k <= K); + + // This means that this sample is not insured by any policy + if (k >= K) { + assert(labelViol > eps); + assert(labelViol < +std::numeric_limits::max()); + if (GET_MAX_VIOL) { + if (labelViol > maxViol) { + maxViol = labelViol; + label_max = label; + } + } + else if (labelViol > eps) { + return false; + } + } + + } + + // return maximum violation + if (GET_MAX_VIOL) { + label = label_max; + if (maxViol > eps) { + return false; + } + } + + return true; +} + +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- + +bool KAdaptableSolver::feasible_DET(const std::vector& x, const std::vector& q) const { + assert(pInfo); + + // vectors must be of appropriate size + if ((int)x.size() != MY_SIZE_X(1)) MYERROR(EXCEPTION_X); + if ((int)q.size() != MY_SIZE_Q) MYERROR(EXCEPTION_Q); + + // Loop through all constraints and get violation + std::vector xx; + if (!feasible_DET_K(x, 1, xx)) { + return 0; + } + assert(xx == x); + for (const auto& con: pInfo->getConstraintsXQ()) { + if (getViolation(con, x, q) > EPS_INFEASIBILITY_Q) return 0; + } + for (const auto& con: pInfo->getConstraintsXYQ()[0]) { + if (getViolation(con, x, q) > EPS_INFEASIBILITY_Q) return 0; + } + + return true; +} + +//----------------------------------------------------------------------------------- + +bool KAdaptableSolver::feasible_SRO(const std::vector& x, std::vector& q) const { + assert(pInfo); + + // vectors must be of appropriate size + if ((int)x.size() != MY_SIZE_X(1)) return 0; + + // SRO = 1-adaptable ARO + return feasible_KAdaptability(x, 1, q); +} + +//----------------------------------------------------------------------------------- + +bool KAdaptableSolver::feasible_SRO(const std::vector& x, const std::vector >& samples, int & label) const { + assert(pInfo); + + // vectors must be of appropriate size + if ((int)x.size() != MY_SIZE_X(1)) MYERROR(EXCEPTION_X); + for (const auto& q : samples) { + if ((int)q.size() != MY_SIZE_Q) MYERROR(EXCEPTION_Q); + } + + // SRO = 1-adaptable ARO + return feasible_KAdaptability(x, 1, samples, label); +} + +//----------------------------------------------------------------------------------- + +bool KAdaptableSolver::feasible_KAdaptability(const std::vector& xx, const unsigned int K1, std::vector& q) const { + + // Get true solution + auto x(xx); + + // Check deterministic feasibility + if (!feasible_DET_K(xx, K1, x)) { + return false; + } + + // Get true # of policies + unsigned int K = K1; + unsigned int Kx = getNumPolicies(x); + if (Kx < K) K = Kx; + assert(K); + + // Obtain worst violation? + std::vector max_q(1, -std::numeric_limits::max()); + + + // Check 1st-stage uncertain constraints + if (!feasible_XQ(x, q)) { + if (GET_MAX_VIOL) { + if (q[0] > max_q[0]) { + max_q = q; + } + } + else { + assert(q[0] > EPS_INFEASIBILITY_Q); + return false; + } + } + else { + max_q = q; + } + + // Check 2nd-stage uncertain constraints + if (!feasible_YQ(x, K, q)) { + if (GET_MAX_VIOL) { + if (q[0] > max_q[0]) { + max_q = q; + } + } + else { + assert(q[0] > EPS_INFEASIBILITY_Q); + return false; + } + } + else { + max_q = q; + } + + // return maximum violation + if (GET_MAX_VIOL) { + q = max_q; + if (q[0] > EPS_INFEASIBILITY_Q) { + return false; + } + } + + return true; +} + +//----------------------------------------------------------------------------------- + +bool KAdaptableSolver::feasible_KAdaptability(const std::vector& xx, const unsigned int K1, const std::vector >& samples, int & label) const { + + // Get true solution + auto x(xx); + + // Check deterministic feasibility + if (!feasible_DET_K(xx, K1, x)) { + return false; + } + + // Get true # of policies + unsigned int K = K1; + unsigned int Kx = getNumPolicies(x); + if (Kx < K) K = Kx; + assert(K); + + + // Check 1st-stage uncertain constraints + if (!feasible_XQ(x, samples, label)) { + return false; + } + + // Check 2nd-stage uncertain constraints + if (!feasible_YQ(x, K, samples, label)) { + return false; + } + + return true; +} + +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- + +int KAdaptableSolver::solve_DET(const std::vector& q, std::vector& x) const { + assert(pInfo); + + // vectors must be of appropriate size + if ((int)q.size() != MY_SIZE_Q) MYERROR(EXCEPTION_Q); + + int status; + CPXENVptr env; + CPXLPptr lp; + + // initialize CPLEX objects + env = NULL; + lp = NULL; + env = CPXXopenCPLEX (&status); + lp = CPXXcreateprob(env, &status, "DETERMINISTIC"); + + // define variables and constraints + updateX(env, lp); + updateXQ(env, lp, q); + updateY(env, lp, 0); + updateYQ(env, lp, 0, q); + + // set options + setCPXoptions(env); + CPXXchgprobtype(env, lp, (pInfo->isContinuous() ? CPXPROB_LP : CPXPROB_MILP)); + CPXXchgobjsen(env, lp, CPX_MIN); + + // solve problem + status = (pInfo->isContinuous() ? CPXXlpopt(env, lp) : CPXXmipopt(env, lp)); + if (!status) { + status = CPXXgetstat(env, lp); + if (status == CPX_STAT_OPTIMAL || status == CPXMIP_OPTIMAL || status == CPXMIP_OPTIMAL_TOL) { + status = 0; + x.resize(CPXXgetnumcols(env, lp)); + CPXXgetx(env, lp, &x[0], 0, x.size() - 1); + assert(feasible_DET(x, q)); + } + } + else { + MYERROR(status); + } + + // Free memory + CPXXfreeprob(env, &lp); + CPXXcloseCPLEX (&env); + + return status; +} + +//----------------------------------------------------------------------------------- + +int KAdaptableSolver::solve_SRO_duality(std::vector& x) const { + assert(pInfo); + + + int status; + CPXENVptr env; + CPXLPptr lp; + + // initialize CPLEX objects + env = NULL; + lp = NULL; + env = CPXXopenCPLEX (&status); + lp = CPXXcreateprob(env, &status, "STATIC_ROBUST"); + + // define variables and constraints + updateX(env, lp); + updateXQ(env, lp); + updateY(env, lp, 0); + updateYQ(env, lp, 0); + + // set options + setCPXoptions(env); + CPXXchgprobtype(env, lp, (pInfo->isContinuous() ? CPXPROB_LP : CPXPROB_MILP)); + CPXXchgobjsen(env, lp, CPX_MIN); + + // solve problem + status = (pInfo->isContinuous() ? CPXXlpopt(env, lp) : CPXXmipopt(env, lp)); + if (!status) { + status = CPXXgetstat(env, lp); + if (status == CPX_STAT_OPTIMAL || status == CPXMIP_OPTIMAL || status == CPXMIP_OPTIMAL_TOL) { + status = 0; + x.resize(MY_SIZE_X(1)); + CPXXgetx(env, lp, &x[0], 0, x.size() - 1); + #ifndef NDEBUG + std::vector q; + assert(feasible_SRO(x, q)); + #endif + } + } + else { + MYERROR(status); + } + + // Free memory + CPXXfreeprob(env, &lp); + CPXXcloseCPLEX (&env); + + return status; +} + +//----------------------------------------------------------------------------------- + +int KAdaptableSolver::solve_SRO_cuttingPlane(std::vector& x) { + assert(pInfo); + + + int status; + CPXENVptr env; + CPXLPptr lp; + + // initialize CPLEX objects + env = NULL; + lp = NULL; + env = CPXXopenCPLEX (&status); + lp = CPXXcreateprob(env, &status, "STATIC_ROBUST_CP"); + + // define variables and constraints using nominal parameters to prevent unboundedness + updateX(env, lp); + updateXQ(env, lp, pInfo->getUncSet().getNominal()); + updateY(env, lp, 0); + updateYQ(env, lp, 0, pInfo->getUncSet().getNominal()); + + // set options + setCPXoptions(env); + CPXXchgprobtype(env, lp, CPXPROB_MILP); // to use callbacks + CPXXchgobjsen(env, lp, CPX_MIN); + + // callbacks + CPXXsetintparam(env, CPX_PARAM_MIPSEARCH, CPX_MIPSEARCH_TRADITIONAL); + CPXXsetintparam(env, CPX_PARAM_MIPCBREDLP, CPX_OFF); + CPXXsetintparam(env, CPX_PARAM_REDUCE, CPX_PREREDUCE_PRIMALONLY); + CPXXsetintparam(env, CPX_PARAM_PRELINEAR, CPX_OFF); + CPXXsetusercutcallbackfunc(env, cutCB_solve_SRO_cuttingPlane, this); + CPXXsetlazyconstraintcallbackfunc(env, cutCB_solve_SRO_cuttingPlane, this); + + // solve problem + status = CPXXmipopt(env, lp); + if (!status) { + status = CPXXgetstat(env, lp); + if (status == CPXMIP_OPTIMAL || status == CPXMIP_OPTIMAL_TOL) { + status = 0; + x.resize(MY_SIZE_X(1)); + CPXXgetx(env, lp, &x[0], 0, x.size() - 1); + #ifndef NDEBUG + std::vector q; + assert(feasible_SRO(x, q)); + #endif + } + } + else { + MYERROR(status); + } + + // Free memory + CPXXfreeprob(env, &lp); + CPXXcloseCPLEX (&env); + + return status; +} + +//----------------------------------------------------------------------------------- + +int KAdaptableSolver::solve_ScSRO(const std::vector >& samples, std::vector& x) const { + assert(pInfo); + + // vectors must be of appropriate size + for (const auto& q : samples) { + if ((int)q.size() != MY_SIZE_Q) MYERROR(EXCEPTION_Q); + } + + int status; + CPXENVptr env; + CPXLPptr lp; + + // initialize CPLEX objects + env = NULL; + lp = NULL; + env = CPXXopenCPLEX (&status); + lp = CPXXcreateprob(env, &status, "SAMPLE_BASED_STATIC_ROBUST"); + + // define variables and constraints + updateX(env, lp); + updateY(env, lp, 0); + for (const auto& q: samples) { + updateXQ(env, lp, q); + updateYQ(env, lp, 0, q); + } + + // set options + setCPXoptions(env); + CPXXchgprobtype(env, lp, (pInfo->isContinuous() ? CPXPROB_LP : CPXPROB_MILP)); + CPXXchgobjsen(env, lp, CPX_MIN); + + // solve problem + status = (pInfo->isContinuous() ? CPXXlpopt(env, lp) : CPXXmipopt(env, lp)); + if (!status) { + status = CPXXgetstat(env, lp); + if (status == CPX_STAT_OPTIMAL || status == CPXMIP_OPTIMAL || status == CPXMIP_OPTIMAL_TOL) { + status = 0; + x.resize(CPXXgetnumcols(env, lp)); + CPXXgetx(env, lp, &x[0], 0, x.size() - 1); + #ifndef NDEBUG + int label; + assert(feasible_SRO(x, samples, label)); + #endif + } + } + else { + MYERROR(status); + } + + // Free memory + CPXXfreeprob(env, &lp); + CPXXcloseCPLEX (&env); + + return status; +} + +//----------------------------------------------------------------------------------- + +int KAdaptableSolver::solve_KAdaptability(const unsigned int K, const bool h, std::vector& x) { + assert(pInfo); + if (K < 1) MYERROR(EXCEPTION_K); + + // Create K set of constraints + pInfo->resize(K); + NK = K; + + // Preconditions for heuristic + heuristic_mode = h; + if (heuristic_mode && K > 1) { + if (getNumPolicies(xfix) + 1 != K) MYERROR(EXCEPTION_X); + } + + int solstat; + int status; + CPXENVptr env; + CPXLPptr lp; + + // initialize CPLEX objects + env = NULL; + lp = NULL; + env = CPXXopenCPLEX (&status); + lp = CPXXcreateprob(env, &status, "KADAPTABILITY_CUTTING_PLANE"); + + + + + + + + // temporary + std::vector qtemp = pInfo->getNominal(), xnom; + + + // solve SRO to attempt to get a tighter primal bound + if (!solve_SRO_duality(xstatic)) { + assert(feasible_KAdaptability(xstatic, K, qtemp)); + if (xstatic[0] < (xsol.empty() ? +std::numeric_limits::max() : xsol[0])) { + xsol = xstatic; + } + } + + + // Get initial sample + double initialObj = -std::numeric_limits::max(); + + + // Get qtemp from SRO + if (!feasible_SRO(xstatic, qtemp)) assert(0); + if (!solve_DET(qtemp, xnom)) initialObj = xnom[0]; + if (!solve_DET(pInfo->getNominal(), xnom)) if (xnom[0] > initialObj) { + initialObj = xnom[0]; + qtemp = pInfo->getNominal(); + } + if (0) if (hasObjectiveUncOnly() && isSecondStageOnly()) { + getLowerBound_2(qtemp); + assert(getLowerBound_2(qtemp) >= initialObj); + } + + // Re-compute qtemp for heuristic + if (heuristic_mode && K > 1) getWorstCase(xfix, K - 1, qtemp); + + // Generate initial scenario + bb_samples.assign(1, qtemp); + + + // assign sample to root node + updateX(env, lp); + for (unsigned int k = 0; k < K; k++) { + updateY(env, lp, k); + } + updateXQ(env, lp, qtemp); + updateYQ(env, lp, 0, qtemp); + + + // set options + setCPXoptions(env); + CPXXsetintparam(env, CPXPARAM_ScreenOutput, CPX_ON); + CPXXchgprobtype(env, lp, CPXPROB_MILP); // to use callbacks + CPXXchgobjsen(env, lp, CPX_MIN); + + + // callbacks + CPXXsetintparam(env, CPX_PARAM_MIPSEARCH, CPX_MIPSEARCH_TRADITIONAL); + CPXXsetintparam(env, CPX_PARAM_MIPCBREDLP, CPX_OFF); + CPXXsetintparam(env, CPX_PARAM_REDUCE, CPX_PREREDUCE_PRIMALONLY); + CPXXsetintparam(env, CPX_PARAM_PRELINEAR, CPX_OFF); + if (!hasObjectiveUncOnly() && !BNC_BRANCH_ALL_CONSTR){ + CPXXsetusercutcallbackfunc(env, cutCB_solve_KAdaptability_cuttingPlane, this); + CPXXsetlazyconstraintcallbackfunc(env, cutCB_solve_KAdaptability_cuttingPlane, this); + } + CPXXsetbranchcallbackfunc(env, branchCB_solve_KAdaptability_cuttingPlane, this); + CPXXsetincumbentcallbackfunc(env, incCB_solve_KAdaptability_cuttingPlane, this); + if (heuristic_mode && K > 1) CPXXsetheuristiccallbackfunc(env, heurCB_solve_KAdaptability_cuttingPlane, this); + CPXXsetdeletenodecallbackfunc(env, deletenodeCB_solve_KAdaptability_cuttingPlane, this); + if (COLLECT_RESULTS && K > 2) { + CPXXsetnodecallbackfunc(env, nodeCB_solve_KAdaptability_cuttingPlane, this); + } + + + + + + + + // add MIP start + if (!xsol.empty()) { + const unsigned int Kx = getNumPolicies(xsol); assert(Kx); + const auto x2 = xsol.begin() + pInfo->getNumFirstStage(); + + // mip start to be added + std::vector x_mipstart(xsol); + + // All second-stage policies are identical + for (unsigned int k = Kx; k < K; k++) { + x_mipstart.insert(x_mipstart.end(), x2, x2 + pInfo->getNumSecondStage()); + } + assert(x_mipstart.size() == MY_SIZE_X(K)); + assert(feasible_KAdaptability(x_mipstart, K, Q_TEMP)); + + // CPLEX structures + std::vector effortlevel{CPX_MIPSTART_CHECKFEAS}; + std::vector beg{0}; + std::vector varindices(x_mipstart.size(), 0); + std::vector xctype(x_mipstart.size(), 'C'); + std::iota(varindices.begin(), varindices.end(), 0); + CPXXgetctype(env, lp, &xctype[0], 0, x_mipstart.size() - 1); + + // pass to CPLEX + CPXXaddmipstarts(env, lp, 1, x_mipstart.size(), &beg[0], &varindices[0], &x_mipstart[0], &effortlevel[0], NULL); + } + + // Fix policies for heuristic + if (heuristic_mode && K > 1) { + assert(getNumPolicies(xfix) + 1 == K); + std::vector xctype(CPXXgetnumcols(env, lp), 'C'); + CPXXgetctype(env, lp, &xctype[0], 0, xctype.size() - 1); + for (int init = pInfo->getNumFirstStage(), i = init; i < (int)xfix.size(); i++) { + double value = xfix[i]; + if (xctype[i] == CPX_BINARY) { + if (std::abs(value - 0.0) <= 1.E-5) value = 0; + if (std::abs(value - 1.0) <= 1.E-5) value = 1.0; + } + fixVariable(env, lp, i + pInfo->getNumSecondStage(), value); + } + } + + // (time, incumbent) data + ZT_VALUES.clear(); + NUM_DUMMY_NODES = TX = 0; + START_TS = get_wall_time(); + BOUND_VALUES.clear(); + + + + + + + + + double start_time = get_wall_time(); + + // solve problem + status = CPXXmipopt(env, lp); + if (!status) { + solstat = CPXXgetstat(env, lp); + if (solstat == CPXMIP_OPTIMAL || solstat == CPXMIP_OPTIMAL_TOL) { + solstat = 0; + } + x.resize(MY_SIZE_X(K)); + if (CPXXgetx(env, lp, &x[0], 0, x.size() - 1) == 0) { + assert(feasible_KAdaptability(x, K, Q_TEMP)); + } + } + else { + MYERROR(status); + } + + double end_time = get_wall_time(); + + + // Resize xfix for heuristic + if (heuristic_mode) { + xfix = this->xsol; + resizeX(xfix, K); + } + + + + // get results + if (COLLECT_RESULTS) { + double total_solution_time = end_time - start_time; + double global_lb; CPXXgetbestobjval(env, lp, &global_lb); + double final_gap = INFINITY, final_objval = INFINITY; + if (CPXXgetobjval(env, lp, &final_objval) == 0) { + final_gap = 100*(final_objval - global_lb)/(1E-10 + std::abs(final_objval)); + } + + + + size_t pos_keyword = 0, pos_dash = 0; + const std::string s = pInfo->getSolFileName(); + const std::string problemType(s, 0, 3); + pos_keyword = s.find_first_of("n", 3); + pos_dash = s.find_first_of("-", pos_keyword); + const std::string n(s, pos_keyword + 1, pos_dash - pos_keyword - 1); + pos_keyword = s.find_first_of("s", pos_dash); + pos_dash = s.find_first_of("-.opt", pos_keyword); + const std::string seed(s, pos_keyword + 1, pos_dash - pos_keyword - 1); + std::string stat = "Unknown"; + if (solstat == 0) stat = "Optimal"; + if (solstat == CPXMIP_INFEASIBLE || solstat == CPXMIP_INForUNBD) stat = "Infeas"; + if (solstat == CPXMIP_TIME_LIM_FEAS || solstat == CPXMIP_TIME_LIM_INFEAS) stat = "Time Lim"; + if (solstat == CPXMIP_MEM_LIM_FEAS || solstat == CPXMIP_MEM_LIM_INFEAS) stat = "Mem Lim"; + if (heuristic_mode) stat = "Heur"; + write(std::cout, n, K, seed, stat, final_objval, total_solution_time, final_gap); + + + + } + + + // clear samples + bb_samples.clear(); + + + // clear (time, incumbent) data + ZT_VALUES.clear(); + + + // Free memory + CPXXfreeprob(env, &lp); + CPXXcloseCPLEX (&env); + + return solstat; +} + +//----------------------------------------------------------------------------------- + +bool KAdaptableSolver::solve_separationProblem(const std::vector& x, const unsigned int K, int& label, bool heur) { + + if (SEPARATE_FROM_SAMPLES ? (!feasible_YQ(x, K, bb_samples, label, heur)) : false) { + assert(label < static_cast(bb_samples.size())); + return true; + } + + if (heur) return false; + + if (!feasible_YQ(x, K, Q_TEMP)) { + bb_samples.emplace_back(Q_TEMP); + label = static_cast(bb_samples.size()) - 1; + return true; + } + + return false; +} + +//----------------------------------------------------------------------------------- + +int KAdaptableSolver::solve_KAdaptability_minMaxMin(const unsigned int K, std::vector& x) { + assert(K > 0); + assert(pInfo); + assert(pInfo->hasObjectiveUncOnly()); + assert(pInfo->isSecondStageOnly()); + assert(pInfo->getConstraintsXYQ()[0].size() == 1); + if (!pInfo->hasObjectiveUncOnly()) return 1; + if (!pInfo->isSecondStageOnly()) return 1; + if (pInfo->getConstraintsXYQ()[0].size() != 1) return 1; + if (solve_SRO_duality(xstatic)) MYERROR(1); + + + // Algorithm 1 from Buccheim & Kurtz, Math Prog. (2017) + int status; + unsigned int i = 0; + double zstar = +std::numeric_limits::max(); + std::vector qstar = pInfo->getNominal(), xtemp, xm3{xstatic[0] + 1.0}, lambda; + do { + // line 3 + if (i) { + pInfo->resize(i); + zstar = getWorstCase(xm3, i, qstar); + } + + // line 4 + status = solve_DET(qstar, xtemp); + if (status) MYERROR(status); + + // line 5 + i++; + xm3.insert(xm3.end(), xtemp.begin() + 1, xtemp.end()); + assert(xm3.size() == MY_SIZE_X(i)); + } while (xtemp[0] < zstar - EPS_INFEASIBILITY_Q); + + + //////////////////////////////////////////////////////////////////////// + // line 7 // + // Note that the paper is a bit vague // + // The true problem to be solved is as follows: // + // find a continuous solution x* = conv(x[j], j = 1,...,i) // + // such that the worst-case objective evaluated at x* = z* // + //////////////////////////////////////////////////////////////////////// + CPXENVptr env = NULL; + CPXLPptr lp = NULL; + env = CPXXopenCPLEX(&status); + lp = CPXXcreateprob(env, &status, "MIN_MAX_MIN"); + addVariable(env, lp, 'C', zstar, zstar, 0.0, "zstar"); + for (int n = 1; n <= pInfo->getNumSecondStage(); n++) { + addVariable(env, lp, 'C', -CPX_INFBOUND, +CPX_INFBOUND, 0.0, "xstar(" + std::to_string(n) + ")"); + } + assert(CPXXgetnumcols(env, lp) == (int)xtemp.size()); + for (unsigned int j = 0; j < i; j++) { + addVariable(env, lp, 'C', 0.0, +CPX_INFBOUND, 0.0, "lambda(" + std::to_string(j) + ")"); + } + ConstraintExpression constraint; + constraint.rowname("convex"); + constraint.sign('E'); + constraint.RHS(1.0); + for (unsigned int j = 0; j < i; j++) { + constraint.addTermX(xtemp.size() + j, 1.0); + } + constraint.addToCplex(env, lp); + for (int n = 1; n <= pInfo->getNumSecondStage(); n++) { + constraint.clear(); + constraint.RHS(0); + constraint.sign('E'); + constraint.rowname("xdef(" + std::to_string(n) + ")"); + for (unsigned int j = 0; j < i; j++) { + double coef = getXPolicy(xm3, i, j)[n]; + if (coef != 0.0) constraint.addTermX(xtemp.size() + j, coef); + } + if (constraint.isEmpty()) { + fixVariable(env, lp, n, 0); + } + else { + constraint.addTermX(n, -1.0); + constraint.addToCplex(env, lp); + } + } + // This is where we add the worst-case constraint + pInfo->getConstraintsXYQ()[0][0].addToCplex(env, lp, &pInfo->getUncSet(), true); + + + // set options + setCPXoptions(env); + CPXXchgprobtype(env, lp, CPXPROB_LP); + CPXXchgobjsen(env, lp, CPX_MIN); + + // solve problem + status = CPXXdualopt(env, lp); + if (status == 0 && CPXXgetstat(env, lp) == CPX_STAT_OPTIMAL) { + lambda.resize(i); + CPXXgetx(env, lp, &lambda[0], xtemp.size(), xtemp.size() + i - 1); + assert(std::abs(std::accumulate(lambda.begin(), lambda.end(), 0.0) - 1.0) <= EPS_INFEASIBILITY_X); + } + else { + MYERROR((int)(status ? status : CPXXgetstat(env, lp))); + } + + // Free memory + CPXXfreeprob(env, &lp); + CPXXcloseCPLEX (&env); + + // Get K largest lambda + x.assign(1, xstatic[0] + 1.0); + std::vector order(i, 0); + std::iota(order.begin(), order.end(), 0); + std::sort(order.begin(), order.end(), [&](const unsigned int& a, const unsigned int& b) { + return lambda[b] < lambda[a]; + }); + for (unsigned int j = 0; j < std::min(K, i); j++) { + // std::cout << j << " ---> " << lambda[order[j]] << "\n"; + const auto xk = getXPolicy(xm3, i, order[j]); + x.insert(x.end(), xk.begin() + 1, xk.end()); + } + if (i <= K) { + pInfo->resize(K); + resizeX(x, K); + } + x[0] = getWorstCase(x, K, qstar); + + return status; +} + +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- +//----------------------------------------------------------------------------------- + +/** + * Data structure to be attached to every node of the branch-and-bound tree + * to implement K-ary branching + * + * The structure keeps track of whether the node is a 'dummy node' + * as well scenario information to complete the branching step at descendant nodes + * + * There will always be max(K-2, 0) dummy nodes created at every branching step + * + * For example, 3-way branching will be implemented as follows: + * N + * +-+-+ + * c3 d1 + * +-+-+ + * c2 c1 + * where N denotes the current node to be branched on, + * c1, c2, c3 denote the 3 children nodes we wish to create + * d1 denotes a dummy node which will be further branched to create nodes c1 and c2 + * + * Dummy nodes are created without defining any additional branching constraints + */ +struct CPLEX_CB_node { + /** + * Measures true depth of nodes -- for K <= 2, this must be equal to depth returned by CPLEX + */ + unsigned int trueDepth; + + /** + * Indicates that the last branching decision was a CPLEX 0-1 branching decision + */ + bool br_flag; + + /** + * Indicates node be immediately branched on -- recall that CPLEX does not support multi-way branching + */ + bool isDummy; + + /** + * Node is a direct descendant of incumbent CB - accept or reject by branching accordingly + */ + bool fromIncumbentCB; + + /** + * For dummy nodes -- scenario label to be attached to children nodes + * From incumbent CB -- if feasible then -1, otherwise equal to violating scenario label + * For other cases it does not have any meaning + */ + int br_label; + + /** + * Makes sense only for dummy nodes -- objective value predicted for children nodes + */ + double nodeObjective; + + /** + * Makes sense only for dummy nodes -- # of "true" children nodes of this node + */ + unsigned int numNodes; + + /** + * # of active policies in this node + */ + unsigned int numActivePolicies; + + /** + * [k] = Scenario labels for policy k in this node -- to be used in cut and lazy constraint callbacks + * if objective uncertainty then always empty (for efficiency reasons), otherwise has size K + */ + std::vector > labels; +}; + +//----------------------------------------------------------------------------------- + +static int CPXPUBLIC cutCB_solve_SRO_cuttingPlane(CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, int *useraction_p) { + + enterCallback(cutCB_solve_SRO_cuttingPlane); + + // Get K-Adaptability solver + auto S = static_cast(cbhandle); + + // Get # of columns + CPXCLPptr lp = NULL; CPXXgetcallbacklp(env, cbdata, wherefrom, &lp); + CPXDIM numcols = CPXXgetnumcols(env, lp); + + // get node solution + std::vector x(numcols); CPXXgetcallbacknodex(env, cbdata, wherefrom, &x[0], 0, numcols - 1); + std::vector q; + *useraction_p = CPX_CALLBACK_DEFAULT; + + + // violated cuts + CPXDIM rcnt = 0; + CPXNNZ nzcnt; + std::vector rhs; + std::vector sense; + std::vector rmatbeg; + std::vector rmatind; + std::vector rmatval; + + // check constraints (x, q) + if (!S->feasible_XQ(x, q)) { + + // get all constraints + S->getXQ_fixedQ(q, rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval); + } + // check constraints (x, y, q) + else if (!S->feasible_YQ(x, 1, q)) { + + // get all constraints + S->getYQ_fixedQ(0, q, rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval); + } + else { + assert(S->feasible_SRO(x, q)); + } + + + + // Add only the most violated constraint + double maxViol = 0; + double rhs_cut = 0; + char sense_cut = 'L'; + std::vector cutind; + std::vector cutval; + for (CPXDIM i = 0; i < rcnt; i++) { + + // coefficients and indices of this constraint + CPXDIM next = (i + 1 == rcnt) ? nzcnt : rmatbeg[i+1]; + std::vector cutind_i(rmatind.begin() + rmatbeg[i], rmatind.begin() + next); + std::vector cutval_i(rmatval.begin() + rmatbeg[i], rmatval.begin() + next); + + // check violation + const auto W = checkViol(env, cbdata, wherefrom, rhs[i], sense[i], cutind_i, cutval_i); + if (W.first > maxViol) { + maxViol = W.first; + rhs_cut = rhs[i]; + sense_cut = sense[i]; + cutind = cutind_i; + cutval = cutval_i; + } + } + + // Add local cut to be consistent with K-Adaptability implementation + if (maxViol > EPS_INFEASIBILITY_Q) { + CPXXcutcallbackaddlocal(env, cbdata, wherefrom, cutind.size(), rhs_cut, sense_cut, &cutind[0], &cutval[0]);//, CPX_USECUT_FORCE); + *useraction_p = CPX_CALLBACK_SET; + } + + + exitCallback(cutCB_solve_SRO_cuttingPlane); +} + +//----------------------------------------------------------------------------------- + +static void CPXPUBLIC deletenodeCB_solve_KAdaptability_cuttingPlane(CPXCENVptr, int, void *, CPXCNT, void *handle) { + if (handle) { + auto oldInfo = static_cast(handle); + delete oldInfo; + } +} + +//----------------------------------------------------------------------------------- + +static int CPXPUBLIC incCB_solve_KAdaptability_cuttingPlane(CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, double objval, double *x_, int *isfeas_p, int *useraction_p) { + enterCallback(inc); + *useraction_p = CPX_CALLBACK_SET; + + // Get K-Adaptability solver + auto S = static_cast(cbhandle); + const unsigned int K = S->NK; + + + // Get node data + void *nodeData = NULL; + + + // Get node solution + CPXCLPptr lp = NULL; CPXXgetcallbacklp(env, cbdata, wherefrom, &lp); + CPXDIM numcols = CPXXgetnumcols(env, lp); + std::vector x(x_, x_ + numcols); + int label = 0; + + + // Solution must be structurally feasible + // assert(S->feasible_DET_K(x, K, X_TEMP)); + // assert(x == X_TEMP); + + + /////////////////////////////////////////////////////////////////////////// + // NOTE -- I AM ASSUMING THAT ALL CONSTRAINTS (X,Q) HAVE BEEN DUALIZED!! // + /////////////////////////////////////////////////////////////////////////// + // assert(S->feasible_XQ(x, Q_TEMP)); + + + // Exit if node is a dummy node + if (wherefrom == CPX_CALLBACK_MIP_INCUMBENT_NODESOLN) { + CPXXgetcallbacknodeinfo(env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_USERHANDLE, &nodeData); + if (nodeData) if (static_cast(nodeData)->isDummy) { + assert(K > 2); + assert(static_cast(nodeData)->numNodes > 1); + assert(static_cast(nodeData)->numNodes < K); + assert(!static_cast(nodeData)->fromIncumbentCB); + assert(!S->feasible_YQ(x, K, Q_TEMP)); + *isfeas_p = 0; + exitCallback(inc); + } + } + + + + // Check feasibility + *isfeas_p = !S->solve_separationProblem(x, K, label); + if (*isfeas_p) { + S->setX(x, K); + if (COLLECT_RESULTS) ZT_VALUES.emplace_back(objval, get_wall_time()); + } + + + // Tag this node -- to be identified in branch CB + if (wherefrom == CPX_CALLBACK_MIP_INCUMBENT_NODESOLN) { + #ifndef NDEBUG + CPXLONG depth = 0; CPXXgetcallbacknodeinfo(env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_DEPTH_LONG, &depth); + #endif + CPXXgetcallbacknodeinfo(env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_USERHANDLE, &nodeData); + + // new node data to be attached + CPLEX_CB_node* newInfo = new CPLEX_CB_node; + newInfo->trueDepth = 0; + newInfo->br_flag = 0; + newInfo->isDummy = 0; + newInfo->fromIncumbentCB = 1; + newInfo->br_label = (*isfeas_p ? -1 : label); + newInfo->nodeObjective = objval; + newInfo->numNodes = 0; + newInfo->numActivePolicies = (S->heuristic_mode ? K : 1); + if (S->hasObjectiveUncOnly()) { + newInfo->labels.clear(); + } else { + newInfo->labels.assign(K, std::vector{}); + newInfo->labels[0].emplace_back(0); + } + + CPXXcallbacksetuserhandle(env, cbdata, wherefrom, newInfo, &nodeData); + + // delete old data + if (nodeData) { + auto oldInfo = static_cast(nodeData); + newInfo->trueDepth = oldInfo->trueDepth; + assert((K > 2) || (oldInfo->trueDepth == depth)); + assert(!oldInfo->isDummy); + assert(oldInfo->numActivePolicies <= K); + assert(oldInfo->labels.size() == (S->hasObjectiveUncOnly() ? 0 : K)); + newInfo->numActivePolicies = oldInfo->numActivePolicies; + newInfo->labels = oldInfo->labels; + delete oldInfo; + } + // no data here -- can only happen for root node + // This is the first 'tag' of the B&B tree + else { + assert(depth == 0); + } + } + + exitCallback(inc); +} + +//----------------------------------------------------------------------------------- + +static int CPXPUBLIC heurCB_solve_KAdaptability_cuttingPlane (CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, double *objval_p, double *x_, int *checkfeas_p, int *useraction_p) { + enterCallback(heur); + + // Get K-Adaptability solver + auto S = static_cast(cbhandle); + const unsigned int K = S->NK; + + + // Get node solution + CPXCLPptr lp = NULL; CPXXgetcallbacklp(env, cbdata, wherefrom, &lp); + CPXDIM numcols = CPXXgetnumcols(env, lp); + std::vector x(x_, x_ + numcols); + + + // Check integer feasibility + std::vector feas(numcols, 0); CPXXgetcallbacknodeintfeas(env, cbdata, wherefrom, &feas[0], 0, numcols - 1); + if (std::accumulate(feas.begin(), feas.end(), 0) == 0) { + double ub = 0; CPXXgetcallbackinfo(env, cbdata, wherefrom, CPX_CALLBACK_INFO_BEST_INTEGER, &ub); + x[0] = ub; + x[0] = S->getWorstCase(x, K, Q_TEMP); + + // worst-case objective value is less than current incumbent + if (x[0] < ub - EPS_INFEASIBILITY_X) { + x_[0] = x[0] + EPS_INFEASIBILITY_X; + *objval_p = x_[0]; + *checkfeas_p = 1; + *useraction_p = CPX_CALLBACK_SET; + } + else { + *useraction_p = CPX_CALLBACK_DEFAULT; + } + } + + exitCallback(heur); +} + +//----------------------------------------------------------------------------------- + +static int CPXPUBLIC branchCB_solve_KAdaptability_cuttingPlane(CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, int type, CPXDIM, int nodecnt, CPXDIM bdcnt, const CPXDIM *nodebeg, const CPXDIM *indices, const char *lu, const double *bd, const double *nodeest, int *useraction_p) { + enterCallback(branch); + *useraction_p = CPX_CALLBACK_DEFAULT; + + // Get K-Adaptability solver + auto S = static_cast(cbhandle); + const unsigned int K = S->NK; + + + // Get node data + void *nodeData = NULL; CPXXgetcallbacknodeinfo(env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_USERHANDLE, &nodeData); + CPXLONG depth = 0; CPXXgetcallbacknodeinfo(env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_DEPTH_LONG, &depth); + double nodeobjval = 0; CPXXgetcallbacknodeinfo(env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_OBJVAL, &nodeobjval); + int existsfeas = 0; CPXXgetcallbackinfo(env, cbdata, wherefrom, CPX_CALLBACK_INFO_MIP_FEAS, &existsfeas); + double bestinteger = 0;CPXXgetcallbackinfo(env, cbdata, wherefrom, CPX_CALLBACK_INFO_BEST_INTEGER, &bestinteger); + + double gap = +std::numeric_limits::max(); + if (existsfeas) { + gap = 100*(bestinteger - nodeobjval)/(1E-10 + std::abs(bestinteger)); + } + if (COLLECT_RESULTS) { + double time_elapsed = get_wall_time() - START_TS; + if (time_elapsed > (double)TX*10) { + double lb = 0; CPXXgetcallbackinfo(env, cbdata, wherefrom, CPX_CALLBACK_INFO_BEST_REMAINING, &lb); + BOUND_VALUES.emplace_back(lb, (existsfeas ? bestinteger : 0)); + TX++; + } + } + + + // Get node solution + CPXCLPptr lp = NULL; CPXXgetcallbacklp(env, cbdata, wherefrom, &lp); + CPXDIM numcols = CPXXgetnumcols(env, lp); + std::vector x(numcols); CPXXgetcallbacknodex(env, cbdata, wherefrom, &x[0], 0, numcols - 1); + + + + /////////////////////////////////////////////////////////////////////////// + // NOTE -- I AM ASSUMING THAT ALL CONSTRAINTS (X,Q) HAVE BEEN DUALIZED!! // + /////////////////////////////////////////////////////////////////////////// + // assert(S->feasible_XQ(x, Q_TEMP)); + + + + + + + + + + + + + + /////////////////////////////////////////////////////////////////// + // DECIDE WHETHER TO RUN THE SEPARATION PROBLEM IN THE BRANCH CB // + /////////////////////////////////////////////////////////////////// + // non-zero value for label will indicate + // the sequence number of violating scenario to branch upon + int label = 0; + bool compute_separation = false; + bool possible_branching = false; + + + ///////////////////////////////// + // CASE 1: NO NODE DATA EXISTS // + ///////////////////////////////// + if (!nodeData) { + // NOTE: this can only happen for root node + // NOTE: this cannot be from an incumbent CB -- because the latter always tags each node + assert(depth == 0); + compute_separation = (BRANCHING_STRATEGY >= 2); + } + ////////////////////////////////////// + // CASE 2: NODE DATA ALREADY EXISTS // + ////////////////////////////////////// + else { + auto oldInfo = static_cast(nodeData); + assert((K > 2) || (oldInfo->trueDepth == depth)); + assert(oldInfo->numActivePolicies >= 1); + assert(oldInfo->numActivePolicies <= K); + assert(oldInfo->labels.size() == (S->hasObjectiveUncOnly() ? 0 : K)); + + + + // if the incumbent CB rejected the solution and + // if the label suggested does not actually violate x + // then re-compute a new label + if (oldInfo->fromIncumbentCB && oldInfo->br_label >= 0) { + + // check for infeasibility with this label or re-compute + std::vector > qcheck{S->bb_samples.at(oldInfo->br_label)}; + compute_separation = S->feasible_YQ(x, K, qcheck, LABEL_TEMP); + if (!compute_separation) label = oldInfo->br_label; + } + + + // if the node is not a dummy node and + // if the node is not coming directly from the incumbent CB + // then attempt to branch as per k-adaptability + if (!oldInfo->isDummy && !oldInfo->fromIncumbentCB) { + switch (BRANCHING_STRATEGY) { + case 1: compute_separation = false; break; + case 2: compute_separation = oldInfo->br_flag; break; + case 3: compute_separation = (gap > BNC_GAP_VALUE); break; + case 4: compute_separation = true; break; + case 5: compute_separation = ((K == 1) ? 1 : !((int)oldInfo->trueDepth % (1 + (int)std::ceil(K/2)))); break; + } + if (compute_separation) possible_branching = true; + } + } + + + ///////////////////////////////////////////// + // SOLVE THE SEPARATION PROBLEM FOR THIS X // + ///////////////////////////////////////////// + if (compute_separation) { + if (S->solve_separationProblem(x, K, label, possible_branching)) { + assert(label); + } + else { + label = 0; + + + // if the incumbent CB rejected the solution, + // then we must be able to branch as per k-adaptability + #ifndef NDEBUG + if (nodeData) { + CPXINT numiinf; CPXXgetcallbacknodeinfo(env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_NIINF, &numiinf) ; + auto oldInfo = static_cast(nodeData); + assert(!(oldInfo->fromIncumbentCB && oldInfo->br_label >= 0) || numiinf); + } + #endif + } + } + + + + + + + + + + + + + + + + + + /////////////////////////////////////////////////////////////// + // CHECK IF IT IS WORTH TRYING OUT THE K-ADAPTABILITY BRANCH // + /////////////////////////////////////////////////////////////// + if (possible_branching && label) { + // don't use K-Adaptability branch if CPLEX is proposing something 'unique' + if (type != CPX_TYPE_VAR) label = 0; + } + if (BNC_DO_STRONG_BRANCH && gap <= BNC_GAP_VALUE) if (possible_branching && label) { + assert(nodecnt); + int status = 0; + CPXLPptr nodelp = NULL; CPXXgetcallbacknodelp(env, cbdata, wherefrom, &nodelp); + + // don't mess with original nodelp -- remember to free memory at the end + CPXLPptr lpcopy = CPXXcloneprob(env, nodelp, &status); + + // get current basis + std::vector cstat(CPXXgetnumcols(env, nodelp)); + std::vector rstat(CPXXgetnumrows(env, nodelp)); + status = CPXXgetbase(env, nodelp, &cstat[0], &rstat[0]); + + // get K-Adaptability branches + std::vector rmatbeg_cut(K, 0); + std::vector rhs_cut(K, 0); + std::vector sense_cut(K, 'L'); + std::vector > cutind(K); + std::vector > cutval(K); + for (unsigned int k = 0; k < K; k++) { + auto xk = S->getXPolicy(x, K, k); + + // violated constraints + CPXDIM rcnt = 0; + CPXNNZ nzcnt; + std::vector rhs; + std::vector sense; + std::vector rmatbeg; + std::vector rmatind; + std::vector rmatval; + + S->getYQ_fixedQ(k, S->bb_samples[label], rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval); + + + + // Add only the most violated constraint + double maxViol = 0; + for (CPXDIM i = 0; i < rcnt; i++) { + + // coefficients and indices of this constraint + CPXDIM next = (i + 1 == rcnt) ? nzcnt : rmatbeg[i+1]; + std::vector cutind_i(rmatind.begin() + rmatbeg[i], rmatind.begin() + next); + std::vector cutval_i(rmatval.begin() + rmatbeg[i], rmatval.begin() + next); + + // check violation + const auto W = checkViol(env, cbdata, wherefrom, rhs[i], sense[i], cutind_i, cutval_i); + if (W.first > maxViol) { + maxViol = W.first; + rhs_cut[k] = rhs[i]; + sense_cut[k] = sense[i]; + cutind[k] = cutind_i; + cutval[k] = cutval_i; + } + } + + // Violation must exist + assert (maxViol > EPS_INFEASIBILITY_Q); + } + + + // EVALUATE CPLEX branches first + double minChildrenBound_CPLEX = +std::numeric_limits::max(); + double maxChildrenBound_CPLEX = -std::numeric_limits::max(); + for (int i = 0; i < nodecnt; i++) { + + // copy basis + if (!status) CPXXcopybase(env, lpcopy, &cstat[0], &rstat[0]); + + // add CPLEX sepecified rows + CPXDIM numrowsadded = 0; + for (CPXDIM j = nodebeg[i], next = ((i + 1 == nodecnt) ? bdcnt : nodebeg[i+1]); j < next; j++) { + ConstraintExpression bound; + bound.addTermX(indices[j], 1); + bound.RHS(bd[j]); + switch (lu[j]) { + case 'L': bound.sign('G'); break; + case 'U': bound.sign('L'); break; + case 'B': bound.sign('E'); break; + } + bound.addToCplex(env, lpcopy); + numrowsadded++; + } + + // solve the LP + CPXXdualopt(env, lpcopy); + if (CPXXgetstat(env, lpcopy) == CPX_STAT_OPTIMAL) { + double objval; CPXXgetobjval(env, lpcopy, &objval); + objval = objval - nodeobjval; + if (objval < minChildrenBound_CPLEX) { + minChildrenBound_CPLEX = objval; + } + if (objval > maxChildrenBound_CPLEX) { + maxChildrenBound_CPLEX = objval; + } + } + + // Delete last few rows + if (numrowsadded) CPXXdelrows(env, lpcopy, rstat.size(), (CPXDIM)rstat.size() + numrowsadded - 1); + } + + + // NEXT, evaluate K-Adaptability branches + double minChildrenBound_KAd = +std::numeric_limits::max(); + double maxChildrenBound_KAd = -std::numeric_limits::max(); + for (unsigned int k = 0; k < K; k++) { + + // copy basis + if (!status) CPXXcopybase(env, lpcopy, &cstat[0], &rstat[0]); + + // add K-Adaptability branches + CPXXaddrows(env, lpcopy, 0, 1, cutind[k].size(), &rhs_cut[k], &sense_cut[k], &rmatbeg_cut[k], &cutind[k][0], &cutval[k][0], NULL, NULL); + + // solve the LP + CPXXdualopt(env, lpcopy); + if (CPXXgetstat(env, lpcopy) == CPX_STAT_OPTIMAL) { + double objval; CPXXgetobjval(env, lpcopy, &objval); + objval = objval - nodeobjval; + if (objval < minChildrenBound_KAd) { + minChildrenBound_KAd = objval; + } + if (objval > maxChildrenBound_KAd) { + maxChildrenBound_KAd = objval; + } + } + + // Delete last row + CPXXdelrows(env, lpcopy, rstat.size(), rstat.size()); + } + + + // free allocated memory + CPXXfreeprob(env, &lpcopy); + + + + // compute CPLEX score + double mu = 0.5; + double score_CPLEX = (mu * minChildrenBound_CPLEX) + ((1 - mu) * maxChildrenBound_CPLEX); + double score_KAd = (mu * minChildrenBound_KAd) + ((1 - mu) * maxChildrenBound_KAd); + + static int count1 = 0; + if (score_CPLEX >= score_KAd) label = 0; + else { + if (++count1%1000 == 0) std::cout << count1 << "\n"; + } + } + + + + + + + + + + + + + + + + + + //////////////////////////////////////////////////////////////////////// + // COMPUTE NEW VALUES FOR NODE DATA TO BE ATTACHED TO CHILDREN NODES // + // THIS DATA WILL BE ATTACHED IRRESPECTIVE OF THE BRANCHING RULE USED // + //////////////////////////////////////////////////////////////////////// + CPXCNT seqnum = 0; + unsigned int k_max = K; + CPLEX_CB_node *newInfo0 = NULL; + CPLEX_CB_node *newInfo1 = NULL; + + ///////////////////////////////// + // CASE 1: NO NODE DATA EXISTS // + ///////////////////////////////// + if (!nodeData) { + // no option but to branch as CPLEX would + if (!label) { + for (int i = 0; i < nodecnt; i++) { + CPLEX_CB_node *newInfo = new CPLEX_CB_node; + newInfo->trueDepth = depth + 1; + newInfo->br_flag = 1; + newInfo->isDummy = 0; + newInfo->fromIncumbentCB = 0; + newInfo->br_label = 0; + newInfo->nodeObjective = nodeest[i]; + newInfo->numNodes = 0; + newInfo->numActivePolicies = 1; + if (S->hasObjectiveUncOnly()) { + newInfo->labels.clear(); + } else { + newInfo->labels.assign(K, std::vector{}); + newInfo->labels[0].emplace_back(0); + } + CPXXbranchcallbackbranchasCPLEX(env, cbdata, wherefrom, i, newInfo, &seqnum); + } + *useraction_p = CPX_CALLBACK_SET; + exitCallback(branch); + } + + + // create your own branches as per k-adaptability (heuristic_mode only) + else if (S->heuristic_mode && K > 1) { + k_max = K - 1; + + // Node data for policy 0 branch + newInfo1 = new CPLEX_CB_node; + newInfo1->trueDepth = depth + 1; + newInfo1->br_flag = 0; + newInfo1->isDummy = 0; + newInfo1->fromIncumbentCB = 0; + newInfo1->br_label = 0; + newInfo1->nodeObjective = x[0]; + newInfo1->numNodes = 0; + newInfo1->numActivePolicies = K; + if (S->hasObjectiveUncOnly()) { + newInfo1->labels.clear(); + } else { + newInfo1->labels.assign(K, std::vector{}); + newInfo1->labels[0].emplace_back(0); + newInfo1->labels[k_max].emplace_back(label); + } + + + // + // Next, create the second node data + // This is either a dummy node too + // or + // it is the last true child node corresponding to policy number 0 + // + if (K > 2) { + // dummy node + newInfo0 = new CPLEX_CB_node(*newInfo1); + newInfo0->isDummy = 1; + newInfo0->br_label = label; + newInfo0->numNodes = k_max; + if (!S->hasObjectiveUncOnly()) { + newInfo0->labels.assign(K, std::vector{}); + newInfo0->labels[0].emplace_back(0); + } + } + else { + // last child node corresponding to policy number 0 + newInfo0 = new CPLEX_CB_node(*newInfo1); + if (!S->hasObjectiveUncOnly()) { + newInfo0->labels.assign(K, std::vector{}); + newInfo0->labels[0].emplace_back(0); + newInfo0->labels[0].emplace_back(label); + } + } + } + + + // create your own branches as per k-adaptability + else { + k_max = ((K < 2) ? 0 : 1); + + // Node data for policy 0 branch + newInfo0 = new CPLEX_CB_node; + newInfo0->trueDepth = depth + 1; + newInfo0->br_flag = 0; + newInfo0->isDummy = 0; + newInfo0->fromIncumbentCB = 0; + newInfo0->br_label = 0; + newInfo0->nodeObjective = x[0]; + newInfo0->numNodes = 0; + newInfo0->numActivePolicies = 1; + if (S->hasObjectiveUncOnly()) { + newInfo0->labels.clear(); + } else { + newInfo0->labels.assign(K, std::vector{}); + newInfo0->labels[0].emplace_back(0); + newInfo0->labels[0].emplace_back(label); + } + + // Node data for policy 1 branch + if (k_max) { + newInfo1 = new CPLEX_CB_node(*newInfo0); + newInfo1->numActivePolicies = 2; + if (!S->hasObjectiveUncOnly()) { + newInfo1->labels.assign(K, std::vector{}); + newInfo1->labels[0].emplace_back(0); + newInfo1->labels[1].emplace_back(label); + } + } + } + } + ////////////////////////////////////// + // CASE 2: NODE DATA ALREADY EXISTS // + ////////////////////////////////////// + else { + auto oldInfo = static_cast(nodeData); + + + + // ====================== + // SUBCASE 1: DUMMY NODE + // ====================== + if (oldInfo->isDummy) { + assert(K > 2); + assert(oldInfo->numNodes > 1); + assert(oldInfo->numNodes < K); + assert(oldInfo->numActivePolicies >= oldInfo->numNodes); + assert(!oldInfo->fromIncumbentCB); + + + // First create the "true node" + // This corresponds to attaching label to policy number [numNodes - 1] + k_max = oldInfo->numNodes - 1; + label = oldInfo->br_label; + + + // "True node" data + newInfo1 = new CPLEX_CB_node(*oldInfo); + newInfo1->isDummy = 0; + newInfo1->fromIncumbentCB = 0; + newInfo1->br_label = 0; + newInfo1->numNodes = 0; + if (S->hasObjectiveUncOnly()) { + assert(newInfo1->labels.empty()); + } else { + newInfo1->labels[k_max].emplace_back(label); + } + + + + + // + // Next, create the second node data + // This is either a dummy node too + // or + // it is the last true child node corresponding to policy number 0 + // + if (oldInfo->numNodes > 2) { + // dummy node + newInfo0 = new CPLEX_CB_node(*oldInfo); + newInfo0->numNodes--; + } + else { + assert(k_max == 1); + + // last child node corresponding to policy number 0 + newInfo0 = new CPLEX_CB_node(*oldInfo); + newInfo0->isDummy = 0; + newInfo0->fromIncumbentCB = 0; + newInfo0->br_label = 0; + newInfo0->numNodes = 0; + if (S->hasObjectiveUncOnly()) { + assert(newInfo0->labels.empty()); + } else { + newInfo0->labels[0].emplace_back(label); + } + } + } + + + + // ============================================= + // SUBCASE 2: COMING DIRECTLY FROM INCUMBENT CB + // AND FLAGGED AS FEASIBLE + // ============================================= + else if (oldInfo->fromIncumbentCB && oldInfo->br_label == -1) { + // feasible node + assert((K > 2) || (oldInfo->trueDepth == depth)); + assert(!oldInfo->isDummy); + assert(nodecnt == 0); + assert(S->feasible_KAdaptability(x, K, Q_TEMP)); + } + + // ========================================================================== + // SUBCASE 3: INFEASIBLE NODE EITHER BECAUSE + // (A) COMING DIRECTLY FROM INCUMBENT CB AND FLAGGED AS INFEASIBLE OR + // (B) IT DOES NOT SATISFY INTEGRALITY + // ========================================================================== + else { + // infeasible node + if (oldInfo->fromIncumbentCB) assert(!oldInfo->isDummy); + + // no option but to branch as CPLEX would + if (!label) { + for (int i = 0; i < nodecnt; i++) { + CPLEX_CB_node *newInfo = new CPLEX_CB_node(*oldInfo); + newInfo->trueDepth++; + newInfo->br_flag = 1; + CPXXbranchcallbackbranchasCPLEX(env, cbdata, wherefrom, i, newInfo, &seqnum); + } + } + + + // create your own branches as per k-adaptability + else { + k_max = oldInfo->numActivePolicies; + if (k_max >= K) k_max--; + + // SRO case + if (K < 2) { + assert(!k_max); + assert(oldInfo->numActivePolicies == 1); + + // node data corresponding to policy number 0 branch + newInfo0 = new CPLEX_CB_node(*oldInfo); + newInfo0->trueDepth++; + newInfo0->br_flag = 0; + newInfo0->isDummy = 0; + newInfo0->fromIncumbentCB = 0; + newInfo0->br_label = 0; + newInfo0->numNodes = 0; + if (S->hasObjectiveUncOnly()) { + assert(newInfo0->labels.empty()); + } else { + newInfo0->labels[0].emplace_back(label); + } + } + else { + assert(k_max); + + // "True node" data + newInfo1 = new CPLEX_CB_node(*oldInfo); + newInfo1->trueDepth++; + newInfo1->br_flag = 0; + newInfo1->isDummy = 0; + newInfo1->fromIncumbentCB = 0; + newInfo1->br_label = 0; + newInfo1->numNodes = 0; + newInfo1->numActivePolicies = k_max + 1; + if (S->hasObjectiveUncOnly()) { + assert(newInfo1->labels.empty()); + } else { + newInfo1->labels[k_max].emplace_back(label); + } + + + // + // Next, create the second node data + // This is either a dummy node too + // or + // it is the last true child node corresponding to policy number 0 + // + if (k_max > 1) { + // dummy node + newInfo0 = new CPLEX_CB_node(*oldInfo); + newInfo0->trueDepth++; + newInfo0->br_flag = 0; + newInfo0->isDummy = 1; + newInfo0->fromIncumbentCB = 0; + newInfo0->br_label = label; + newInfo0->numNodes = k_max; + } + else { + assert(k_max == 1); + + // last child node corresponding to policy number 0 + newInfo0 = new CPLEX_CB_node(*oldInfo); + newInfo0->trueDepth++; + newInfo0->br_flag = 0; + newInfo0->isDummy = 0; + newInfo0->fromIncumbentCB = 0; + newInfo0->br_label = 0; + newInfo0->numNodes = 0; + if (S->hasObjectiveUncOnly()) { + assert(newInfo0->labels.empty()); + } else { + newInfo0->labels[0].emplace_back(label); + } + } + } + } + } + } + + + + + + + + + + + + + + + + + + + + ///////////////////// + // CREATE BRANCHES // + ///////////////////// + if (k_max == K) { + assert(!label); + assert(!newInfo0); + assert(!newInfo1); + } + else { + assert(label); + assert(label < (int)S->bb_samples.size()); + for (int k = k_max, k_min = ((k_max <= 1) ? 0 : k_max); k >= k_min; k--) { + auto xk = S->getXPolicy(x, K, k); + auto newInfo = (k == 0) ? newInfo0 : newInfo1; assert(newInfo); + + // violated constraints + CPXDIM rcnt = 0; + CPXNNZ nzcnt; + std::vector rhs; + std::vector sense; + std::vector rmatbeg; + std::vector rmatind; + std::vector rmatval; + + S->getYQ_fixedQ(k, S->bb_samples[label], rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval); + + + + // Add only the most violated constraint + if (!BNC_BRANCH_ALL_CONSTR && !S->hasObjectiveUncOnly()) { + double maxViol = -std::numeric_limits::max(); + std::vector rhs_cut = {0}; + std::vector sense_cut = {'L'}; + std::vector cutind; + std::vector cutval; + for (CPXDIM i = 0; i < rcnt; i++) { + + // coefficients and indices of this constraint + CPXDIM next = (i + 1 == rcnt) ? nzcnt : rmatbeg[i+1]; + std::vector cutind_i(rmatind.begin() + rmatbeg[i], rmatind.begin() + next); + std::vector cutval_i(rmatval.begin() + rmatbeg[i], rmatval.begin() + next); + + // check violation + const auto W = checkViol(env, cbdata, wherefrom, rhs[i], sense[i], cutind_i, cutval_i); + if (W.first > maxViol) { + maxViol = W.first; + rhs_cut[0] = rhs[i]; + sense_cut[0] = sense[i]; + cutind = cutind_i; + cutval = cutval_i; + } + } + + // Violation must exist + // UNLESS node is a dummy node and solution changed because of local cuts (for e.g.) + #ifndef NDEBUG + if (maxViol <= EPS_INFEASIBILITY_Q) { + std::cout << maxViol << " !!!!!!!!!!!!!!!!!\n"; + assert(static_cast(nodeData)->isDummy); + } + #endif + + // Create branch + CPXXbranchcallbackbranchconstraints(env, cbdata, wherefrom, 1, cutind.size(), &rhs_cut[0], &sense_cut[0], &rmatbeg[0], &cutind[0], &cutval[0], newInfo->nodeObjective, newInfo, &seqnum); + } + // Add all constraints + else { + CPXXbranchcallbackbranchconstraints(env, cbdata, wherefrom, rcnt, nzcnt, &rhs[0], &sense[0], &rmatbeg[0], &rmatind[0], &rmatval[0], newInfo->nodeObjective, newInfo, &seqnum); + } + } + + // Create branch for dummy node + if (k_max > 1) { + assert(newInfo0); + CPXXbranchcallbackbranchconstraints(env, cbdata, wherefrom, 0, 0, NULL, NULL, NULL, NULL, NULL, newInfo0->nodeObjective, newInfo0, &seqnum); + } + } + + + + + *useraction_p = CPX_CALLBACK_SET; + exitCallback(branch); +} + +//----------------------------------------------------------------------------------- + +static int CPXPUBLIC cutCB_solve_KAdaptability_cuttingPlane(CPXCENVptr env, void *cbdata, int wherefrom, void *cbhandle, int *useraction_p) { + enterCallback(cut); + *useraction_p = CPX_CALLBACK_DEFAULT; + + // Get K-Adaptability solver + auto S = static_cast(cbhandle); + const unsigned int K = S->NK; + + + // Attempt to add cuts for insured scenarios + // Note that in case of objective only uncertainty, + // all constraints are defined at the time of branching + if (S->hasObjectiveUncOnly() || BNC_BRANCH_ALL_CONSTR) exitCallback(cut); + + + // Get node data + void *nodeData = NULL; CPXXgetcallbacknodeinfo(env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_USERHANDLE, &nodeData); + + // no node data -- can happen only for root node + if (!nodeData) exitCallback(cut); + + // cannot be from incumbent CB + auto nodeInfo = static_cast(nodeData); + + + // exit if node is a dummy node + if (nodeInfo->isDummy) { + assert(K > 2); + assert(nodeInfo->numNodes > 1); + assert(nodeInfo->numNodes < K); + assert(nodeInfo->numActivePolicies <= K); + assert(nodeInfo->numActivePolicies >= nodeInfo->numNodes); + assert(!nodeInfo->fromIncumbentCB); + assert(nodeInfo->labels.size() == K); + exitCallback(cut); + } + + + // Get node solution + CPXCLPptr lp = NULL; CPXXgetcallbacklp(env, cbdata, wherefrom, &lp); + CPXDIM numcols = CPXXgetnumcols(env, lp); + std::vector x(numcols); CPXXgetcallbacknodex(env, cbdata, wherefrom, &x[0], 0, numcols - 1); + int label = 0; + + + // Solution must be structurally feasible + // assert(S->feasible_DET_K(x, K, X_TEMP)); + // assert(x == X_TEMP); + + + /////////////////////////////////////////////////////////////////////////// + // NOTE -- I AM ASSUMING THAT ALL CONSTRAINTS (X,Q) HAVE BEEN DUALIZED!! // + /////////////////////////////////////////////////////////////////////////// + // assert(S->feasible_XQ(x, Q_TEMP)); + + + // construct policy k solution + for (unsigned int k = 0; k < K; k++) { + auto xk = S->getXPolicy(x, K, k); + + // scenarios that xk must insure against + std::vector > samples_k; + for (const auto& l : nodeInfo->labels[k]) { + samples_k.emplace_back(S->bb_samples[l]); + } + if (samples_k.empty()) continue; + + + // violated cuts + CPXDIM rcnt = 0; + CPXNNZ nzcnt; + std::vector rhs; + std::vector sense; + std::vector rmatbeg; + std::vector rmatind; + std::vector rmatval; + + + // check constraints (x, y, q) + if (!S->feasible_YQ(xk, 1, samples_k, label)) { + + // get all constraints + S->getYQ_fixedQ(k, samples_k[label], rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval); + } + else { + assert(S->feasible_SRO(xk, samples_k, label)); + } + + + + // Add only the most violated constraint + double maxViol = 0; + double rhs_cut = 0; + char sense_cut = 'L'; + std::vector cutind; + std::vector cutval; + for (CPXDIM i = 0; i < rcnt; i++) { + + // coefficients and indices of this constraint + CPXDIM next = (i + 1 == rcnt) ? nzcnt : rmatbeg[i+1]; + std::vector cutind_i(rmatind.begin() + rmatbeg[i], rmatind.begin() + next); + std::vector cutval_i(rmatval.begin() + rmatbeg[i], rmatval.begin() + next); + + // check violation + const auto W = checkViol(env, cbdata, wherefrom, rhs[i], sense[i], cutind_i, cutval_i); + if (W.first > maxViol) { + maxViol = W.first; + rhs_cut = rhs[i]; + sense_cut = sense[i]; + cutind = cutind_i; + cutval = cutval_i; + } + } + + // Add local cut + if (maxViol > EPS_INFEASIBILITY_Q) { + CPXXcutcallbackaddlocal(env, cbdata, wherefrom, cutind.size(), rhs_cut, sense_cut, &cutind[0], &cutval[0]); + *useraction_p = CPX_CALLBACK_SET; + } + } + + exitCallback(cut); +} + +//----------------------------------------------------------------------------------- + +static int CPXPUBLIC nodeCB_solve_KAdaptability_cuttingPlane(CPXCENVptr env, void *cbdata, int wherefrom, void *, CPXCNT *, int *useraction_p) { + *useraction_p = CPX_CALLBACK_DEFAULT; + + // Get node data + void *nodeData = NULL; CPXXgetcallbacknodeinfo(env, cbdata, wherefrom, 0, CPX_CALLBACK_INFO_NODE_USERHANDLE, &nodeData); + if (nodeData) { + NUM_DUMMY_NODES += (static_cast(nodeData)->isDummy); + } + + return 0; +} \ No newline at end of file diff --git a/src/test.cpp b/src/test.cpp new file mode 100644 index 0000000..ab1acb8 --- /dev/null +++ b/src/test.cpp @@ -0,0 +1,74 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#include "problemInfo_knp.hpp" +#include "problemInfo_spp.hpp" +#include "problemInfo_psp.hpp" +#include "robustSolver.hpp" + +int main (int, char*[]) { + + const bool heuristic_mode = false; + const unsigned int Kmax = 4; + KAdaptableInfo *pInfo; + + try { + // Generate the instance data + KNP data; + KAdaptableInfo_KNP knpInfo; + gen_KNP(data, 10, 0, false); // set to 'true' to allow option of loans + knpInfo.setInstance(data); + + // Clone the KAdaptableInfo object + pInfo = knpInfo.clone(); + + + // Other examples from the paper + /* + SPP data; + KAdaptableInfo_SPP sppInfo; + gen_SPP(data, 20, 0); + sppInfo.setInstance(data); + pInfo = sppInfo.clone(); + */ + /* + PSP data; + KAdaptableInfo_PSP pspInfo; + gen_PSP(data, 3, false); // set to 'true' for piecewise affine decision rules + pspInfo.setInstance(data); + pInfo = pspInfo.clone(); + */ + + // CALL THE SOLVER + KAdaptableSolver S(*pInfo); + std::vector x; + // solving in a loop the 1-,2-,...,K-adaptability problems + // will automatically warm-start from the previous solve + for (unsigned int K = 1; K <= Kmax; K++) { + S.solve_KAdaptability(K, heuristic_mode, x); + } + + // Uncomment to also solve without warm-start -- only in exact mode + // S.solve_KAdaptability(K, false, x); + + delete pInfo; + pInfo = NULL; + } + catch (const int& e) { + std::cerr << "Program ABORTED: Error number " << e << "\n"; + } + + + return 0; +} + diff --git a/src/uncertainty.cpp b/src/uncertainty.cpp new file mode 100644 index 0000000..572084c --- /dev/null +++ b/src/uncertainty.cpp @@ -0,0 +1,452 @@ +/***************************************************************************************/ +/* */ +/* Copyright 2018 by Anirudh Subramanyam, Chrysanthos Gounaris and Wolfram Wiesemann */ +/* */ +/* Licensed under the FreeBSD License (the "License"). */ +/* You may not use this file except in compliance with the License. */ +/* You may obtain a copy of the License at */ +/* */ +/* https://www.freebsd.org/copyright/freebsd-license.html */ +/* */ +/***************************************************************************************/ + + +#include "uncertainty.hpp" +#include "Constants.h" +#include +#include +#include // std::transform +#include + +//---------------------------------------------------------------------------// + +static inline void setCPXoptions(CPXENVptr& env) { + CPXXsetdefaults(env); + CPXXsetintparam(env, CPXPARAM_ScreenOutput, CPX_OFF); + CPXXsetintparam(env, CPXPARAM_Threads, NUM_THREADS); + CPXXsetintparam(env, CPXPARAM_Parallel, CPX_PARALLEL_DETERMINISTIC); + CPXXsetintparam(env, CPXPARAM_ClockType, 2); + CPXXsetdblparam(env, CPXPARAM_TimeLimit, TIME_LIMIT); + CPXXsetdblparam(env, CPXPARAM_MIP_Limits_TreeMemory, MEMORY_LIMIT); +} + +//---------------------------------------------------------------------------// + +static inline void initCPX(CPXENVptr& env) { + int status; + + assert(!env); + + // Try to initialize solver environment + env = CPXXopenCPLEX(&status); + if (!env) throw(EXCEPTION_CPXINIT); + + + // Set solver options + setCPXoptions(env); +} + +//---------------------------------------------------------------------------// + +UncertaintySet::UncertaintySet() { + + int status; + + // Initialize members + N = 0; + nominal.reserve(1000); + nominal.emplace_back(0); + low = nominal; + high = nominal; + + polytope_W.emplace_back(std::vector(1, 0)); + polytope_V.emplace_back(std::vector(1, 0)); + polytope_h.emplace_back(0); + polytope_sense.emplace_back('L'); + + + // initialize CPLEX objects + env = NULL; + lp = NULL; + initCPX(env); + + // Create problem object/model + lp = CPXXcreateprob(env, &status, "UncertaintySet"); + if (!lp) throw(EXCEPTION_CPXINIT); + + // Create objective function variable + const double obj = 1; + const double lo = -CPX_INFBOUND; + const double hi = +CPX_INFBOUND; + const char xctype = 'C'; + const char *colname = "O"; + if (CPXXnewcols(env, lp, 1, &obj, &lo, &hi, &xctype, &colname)) { + throw(EXCEPTION_CPXNEWCOLS); + } + + // Make it linear + CPXXchgprobtype(env, lp, CPXPROB_LP); + + // Make it a maximization problem + CPXXchgobjsen(env, lp, CPX_MAX); +} + + +//---------------------------------------------------------------------------// + + +UncertaintySet::UncertaintySet(const UncertaintySet& U) : + N(U.N), + nominal(U.nominal), + low(U.low), + high(U.high), + polytope_W(U.polytope_W), + polytope_V(U.polytope_V), + polytope_h(U.polytope_h), + polytope_sense(U.polytope_sense) +{ + int status; + + // Try to initialize solver environment + env = NULL; + lp = NULL; + initCPX(env); + + + // Copy solver model object + lp = U.getLPObject(env, &status); + if (status) { + throw(EXCEPTION_CPXINIT); + } +} + +//---------------------------------------------------------------------------// + + +UncertaintySet& UncertaintySet::operator=(const UncertaintySet& U) { + if (this == &U) { + return *this; + } + + N = U.N; + nominal = U.nominal; + low = U.low; + high = U.high; + polytope_W = U.polytope_W; + polytope_V = U.polytope_V; + polytope_h = U.polytope_h; + polytope_sense = U.polytope_sense; + + int status; + + // env, lp are already allocated + assert(env); + assert(lp); + + // delete lp + if (lp) if (CPXXfreeprob (env, &lp)) { + throw(EXCEPTION_CPXEXIT); + } + + // Copy solver model object + lp = U.getLPObject(env, &status); + if (status) { + throw(EXCEPTION_CPXINIT); + } + + return *this; +} + + + + +//---------------------------------------------------------------------------// + + +UncertaintySet::~UncertaintySet() noexcept(false) { + this->clear(); + + assert(env); + assert(lp); + + // Free problem object memory + if (lp) if (CPXXfreeprob(env, &lp)) { + throw(EXCEPTION_CPXEXIT); + } + + + // Free solver memory + if (env) if (CPXXcloseCPLEX(&env)) { + char errmsg[CPXMESSAGEBUFSIZE]; + int status = 0; + CPXXgeterrorstring (env, status, errmsg); + std::cerr << errmsg << std::endl; + throw(EXCEPTION_CPXEXIT); + } +} + + +//---------------------------------------------------------------------------// + +void UncertaintySet::clear() { + assert(env); + assert(lp); + + // get # of cols + const int cur_numcols = CPXXgetnumcols(env, lp); + + // delete all columns except the objective + assert(cur_numcols >= 1); + if (cur_numcols > 1) { + CPXXdelcols(env, lp, 1, CPXXgetnumcols(env, lp) - 1); + } + assert(CPXXgetnumcols(env, lp) == 1); + + // set default options + setCPXoptions(env); + + // Zero all members + N = 0; + nominal.assign(1, 0.0); + low = nominal; + high = nominal; + polytope_W.assign(1, std::vector(1, 0)); + polytope_V.assign(1, std::vector(1, 0)); + polytope_h.assign(1, 0.0); + polytope_sense.assign(1, 'L'); +} + +//---------------------------------------------------------------------------// + + + + +void UncertaintySet::addParam(const double nom, const double lo, const double hi) { + if (nom < lo) + std::cerr << "Warning. Nominal value of uncertain parameter is lower than the lower bound.\n"; + + if (nom > hi) + std::cerr << "Warning. Nominal value of uncertain parameter is greater than the upper bound.\n"; + + N++; + nominal.emplace_back(nom); + low.emplace_back(lo); + high.emplace_back(hi); + + assert(nominal.size() == low.size()); + assert(high.size() == low.size()); + assert(1 + N == (int)low.size()); + + + // update matrices + for (unsigned i = 0; i < polytope_W.size(); ++i) { + polytope_W[i].emplace_back(0); + assert((int)polytope_W[i].size() == 1 + N); + } + + // Upper bound + polytope_W.emplace_back(std::vector(1 + N, 0)); + polytope_V.emplace_back(std::vector(1, 0)); + polytope_W.back().back() = 1; + polytope_sense.emplace_back('L'); + polytope_h.emplace_back(hi); + + // Lower bound + polytope_W.emplace_back(std::vector(1 + N, 0)); + polytope_V.emplace_back(std::vector(1, 0)); + polytope_W.back().back() = 1; + polytope_sense.emplace_back('G'); + polytope_h.emplace_back(lo); + + // Matrix sizes must match + assert(polytope_h.size() == polytope_sense.size()); + assert(polytope_W.size() == polytope_sense.size()); + assert(polytope_V.size() == polytope_sense.size()); + + // Add variable to solver model object + double obj = 0; + const char xctype = 'C'; + const char *colname = std::string("q(" + std::to_string(N) + ")").c_str(); + if (CPXXnewcols(env, lp, 1, &obj, &lo, &hi, &xctype, &colname)) { + throw(EXCEPTION_CPXNEWCOLS); + } + assert(CPXXgetnumcols(env, lp) == 1 + N); + + return; + +} + +//---------------------------------------------------------------------------// + + +void UncertaintySet::addFacet(const std::vector >& data, const char sense, const double rhs) { + + for (unsigned i = 0; i < data.size(); ++i) { + if(data[i].first < 1 || data[i].first > N) { + std::cerr << "Warning: Facet description contains invalid indices of uncertain parameters. "; + std::cerr << "Will not add facet.\n"; + return; + } + } + + // update matrices + polytope_W.emplace_back(std::vector(1 + N, 0)); + polytope_V.emplace_back(std::vector(1, 0)); + polytope_h.emplace_back(rhs); + polytope_sense.emplace_back(sense); + + for (const auto& d : data) + polytope_W.back().at(d.first) = d.second; + + + // Matrix sizes must match + assert(polytope_h.size() == polytope_sense.size()); + assert(polytope_W.size() == polytope_sense.size()); + assert(polytope_V.size() == polytope_sense.size()); + + // Update solver model object + std::vector indices; + std::vector coeffs; + + for (const auto& d : data) { + indices.emplace_back(d.first); + coeffs.emplace_back(d.second); + } + + const std::string rname = "row(" + std::to_string(CPXXgetnumrows(env, lp) + 1) + ")"; + const std::vector rowname({rname.c_str()}); + CPXNNZ rmatbeg = 0; + + if (CPXXaddrows(env, lp, 0, 1, indices.size(), &rhs, &sense, &rmatbeg, &indices[0], &coeffs[0], nullptr, &rowname[0])) + throw(EXCEPTION_CPXNEWROWS); + + return; +} + + +//---------------------------------------------------------------------------// + + +int UncertaintySet::addVariables_DualVars(CPXCENVptr env_, CPXLPptr lp_, const std::string& dualName) const { + /* Assumes that polytope_W, polytope_V, polytope_h and polytope_sense + * have all been allocated and have the same number of rows. + * + * Assumes that these matrices/vectors are all 1-indexed. + */ + const CPXDIM ccnt = polytope_sense.size() - 1; + const std::vector obj(ccnt, 0); + const std::vector xctype(ccnt, 'C'); + std::vector lb, ub; + std::vector cname; + for (int s = 1; s <= ccnt; s++) { + cname.emplace_back("UncSetDual(" + (std::string)dualName + "," + std::to_string(s) + ")"); + switch (polytope_sense.at(s)){ + case 'L': ub.push_back(+CPX_INFBOUND); lb.push_back(0); break; + case 'G': ub.push_back(0); lb.push_back(-CPX_INFBOUND); break; + case 'E': ub.push_back(+CPX_INFBOUND); lb.push_back(-CPX_INFBOUND); break; + } + } + + /* ADD COLUMNS TO CPLEX */ + std::vector colname; for (size_t i = 0; i >& data, std::vector& result) const { + + if(N == 0) return 0; + + int probtype = CPXXgetprobtype(env, lp); + std::vector indices; + std::vector coeffs; + + for (const auto& d : data) { + indices.emplace_back(d.first); + coeffs.emplace_back(d.second); + } + + + // return value + double val = 0; result.assign(1 + N, 0); + + // Set objective function of LP over the uncertainty set + std::vector CplexIndices(1 + N, 0); + std::vector values(1 + N, 0); + std::iota(CplexIndices.begin(), CplexIndices.end(), 0); + for (unsigned i = 0; i < indices.size(); ++i) { + values.at(indices.at(i)) = coeffs.at(i); + } + + // Replace current objective function + CPXXchgobj(env, lp, (CPXDIM)CplexIndices.size(), &CplexIndices[0], &values[0]); + + // Set objective sense + CPXXchgobjsen(env, lp, CPX_MAX); + + + // Solve MILP/LP + switch(probtype) { + case CPXPROB_LP: CPXXlpopt(env, lp); assert(CPXXgetstat(env, lp) == CPX_STAT_OPTIMAL); break; + case CPXPROB_MILP: CPXXmipopt(env, lp); assert(CPXXgetstat(env, lp) == CPXMIP_OPTIMAL || CPXXgetstat(env, lp) == CPXMIP_OPTIMAL_TOL); break; + } + + + // Get optimal value + CPXXgetobjval(env, lp, &val); + + // Get solution vector + CPXXgetx(env, lp, &result[0], 0, N); + + + return val; +} + +//---------------------------------------------------------------------------// + + + +double UncertaintySet::min(const std::vector >& data, std::vector& result) const { + + auto d(data); + for (unsigned i = 0; i < data.size(); ++i) d[i].second *= -1.0; + + return (-1.0*UncertaintySet::max(d, result)); +} + +//---------------------------------------------------------------------------// + +double UncertaintySet::getMaximumValue(const std::vector& ind, const std::vector& coef) const { + std::vector > input; + input.reserve(ind.size()); + + std::transform(ind.begin(), ind.end(), coef.begin(), std::back_inserter(input), + [](int i, double c) { + return std::make_pair(i, c); + } + ); + + std::vector result; + + return UncertaintySet::max(input, result); +} + + +//---------------------------------------------------------------------------// + + + + + + + + + + +