Boxes - Nonlinear constraints

#include <xpress.hpp>
#include <stdexcept>  // For throwing exceptions

using namespace xpress;
using namespace xpress::objects;
using xpress::objects::utils::sum;
using xpress::objects::utils::pow;
using xpress::objects::utils::mul;
using xpress::objects::utils::div;
using xpress::objects::utils::scalarProduct;

/**
 * A craftsman makes small wooden boxes for sale. He has four different types of box,
 * and can make each type in any size (keeping all the dimensions in proportion), but
 * all boxes of the same type must have the same size. The profit he makes on a box
 * depends on the size. He has only a limited amount of the necessary wood available
 * and a limited amount of time in the week to do the work. How many boxes should he
 * make, and what size should they be, in order to maximize his profit?
 */

// A box.
class Box {
public:
    const std::string name;    // The name of this box.
    const double lengthCoeff;  // Relative length of this box.
    const double widthCoeff;   // Relative width of this box.
    const double heightCoeff;  // Relative height of this box.
    const double profitCoeff;  // Coefficient for the profit of this box, relative to its size.
    const double timeCoeff;    // Coefficient for the production time of this box, relative to its ply.

    // Constructor
    Box(std::string name,
        double lengthCoeff, double widthCoeff, double heightCoeff,
        double profitCoeff, double timeCoeff)
        : name(name),
          lengthCoeff(lengthCoeff), widthCoeff(widthCoeff), heightCoeff(heightCoeff),
          profitCoeff(profitCoeff), timeCoeff(timeCoeff) {}

    // Override toString() method
    std::string toString() const {
        return name;
    }

    // resources required: nrBattens = size * 4 * (lengthCoeff + widhtCoeff + heightCoeff)
    double getBattensCoeff() const {
        return 4 * (lengthCoeff + widthCoeff + heightCoeff);
    }

    // resources required: ply = size^2 * 2 * (lengthCoeff * widthCoeff + widthC * heightC + heightC * lengthC)
    double getPlyCoeff() const {
        return 2 * (lengthCoeff * widthCoeff + widthCoeff * heightCoeff + heightCoeff * lengthCoeff);
    }
};

// The boxes used in this example.
const std::vector<Box> boxTypeArray = {
    Box("Cube", 1, 1, 1, 20, 1),
    Box("Oblong", 1, 2, 1, 27.3, 1),
    Box("Flat", 4, 4, 1, 90, 1),
    Box("Economy", 1, 2, 1, 10, 0.2)
};

// The resource constraints used in this example.
const double maxSize = 2.0;
const int maxNumProducedPerBoxType = 6;
const double maxNrBattens = 200.0;
const double maxPly = 210.0;
const double maxTime = 35.0;

int main() {
    try {
        // Create a problem instance with verbose messages printed to Console
        XpressProblem prob;
        prob.callbacks.addMessageCallback(XpressProblem::console);

        // The number of each of the type of boxes that should be produced
        std::vector<Variable> numProduced = prob.addVariables(static_cast<int>(boxTypeArray.size()))
                .withType(ColumnType::Integer)
                .withName("numProduced_%d")
                .withUB(maxNumProducedPerBoxType)
                .toArray();

        // The relative size (a multiplier of length/width/height) of each of the box types to produce
        std::vector<Variable> size = prob.addVariables(static_cast<int>(boxTypeArray.size()))
                .withType(ColumnType::Continuous)  // This is the default value so not required to specify
                .withName("size_%d")
                .withUB(maxSize)
                .toArray();

        /* Resource Availabililty Constraints */

        // Construct Expression for ply required for each type of box
        std::vector<Expression> plyPerBoxType(static_cast<int>(boxTypeArray.size()));
        for (std::size_t i=0; i<boxTypeArray.size(); i++) {
            // plyPerBoxType = numProduced * size^2 * plyCoeff
            plyPerBoxType[i] = mul(numProduced[i], pow(size[i], 2.0)) * boxTypeArray[i].getPlyCoeff();
        }
        // Inspect the constructed Expression
        std::cout << "Total ply: " << sum(plyPerBoxType).toString() << std::endl;
        prob.addConstraint(sum(plyPerBoxType) <= maxPly);

        /* Next, we construct a same type of resource constraint as above, but then for the total battens
         * instead of ply. Except, instead of using sum and mul, we use scalarProduct.
         */
        // First collect the BattensCoefficients into one vector. We do this using the std::transform function
        // along with a lambda function from Box to double (similar to `map` function in Python if you are familiar)
        std::vector<double> battensCoeffs(boxTypeArray.size());
        std::transform(boxTypeArray.begin(), boxTypeArray.end(), battensCoeffs.begin(), [](Box box) {
                       return box.getBattensCoeff();
        });
        // battensPerBox = numProduced * size * battensCoeff
        prob.addConstraint(scalarProduct(numProduced, size, battensCoeffs) <= maxNrBattens);


        // Again a similar resource constraint as above, now for total time.
        // We use a different way to use sum (now with a lambda function)
        Expression totalTime = sum(static_cast<int>(boxTypeArray.size()), [&](int i) {
            // timeNeededForBox = 1 + timeCoeff * 1.5^(ply/10)
            Expression plyNeeded = pow(size[i], 2.0) * boxTypeArray[i].getPlyCoeff();
            Expression timeNeeded = 1 + mul(boxTypeArray[i].timeCoeff, pow(1.5, div(plyNeeded, 10)));
            return numProduced[i] * timeNeeded;
        });
        // Inspect the constructed Expression
        std::cout << "Total time: " << totalTime.toString() << std::endl;
        prob.addConstraint(totalTime <= maxTime);

        /* Construct the objective */
        Expression totalProfit = sum(static_cast<int>(boxTypeArray.size()), [&](int i) {
            // profit = numProduced * size^1.5 * profitCoeff
            return numProduced[i] * pow(size[i], 1.5) * boxTypeArray[i].profitCoeff;
        });
        // Inspect the constructed Expression
        std::cout << "Total profit: " << totalProfit.toString() << std::endl;


        // To make the objective linear, just maximize the objtransfercol and add non-linear constraint
        Variable objtransfercol = prob.addVariable(ColumnType::Continuous, "objTransferCol");
        prob.setObjective(objtransfercol,  xpress::ObjSense::Maximize);
        // Add the objective transfer row: objtransfercol = totalProfit
        prob.addConstraint(objtransfercol == totalProfit);

        // Dump the problem to disk so that we can inspect it.
        prob.writeProb("Boxes02.lp");


        // By default we will solve to global optimality, uncomment for an MISLP solve to local optimality
        // prob.setNLPSolver(XPRS_NLPSOLVER_LOCAL);

        // Solve
        prob.optimize();

        // Check the solution status
        if (prob.attributes.getSolStatus() != SolStatus::Optimal && prob.attributes.getSolStatus() != SolStatus::Feasible) {
            std::ostringstream oss; oss << prob.attributes.getSolStatus(); // Convert xpress::SolStatus to String
            throw std::runtime_error("Optimization failed with status " + oss.str());
        }

        // Get the solution and print it
        std::vector<double> sol = prob.getSolution();
        std::cout << std::endl << "*** Solution ***" << std::endl;
        // Printing the objective is currently not supported for nonlinear
        // std::cout << "Objective: " + prob.getNLPObjVal();
        // Print out the variables
        for (std::size_t i = 0; i < boxTypeArray.size(); ++i) {
            if (numProduced[i].getValue(sol) > 0.0) {
                std::cout << "Producing " << numProduced[i].getValue(sol) << " "
                          << boxTypeArray[i].toString() << " boxes of size "
                          << size[i].getValue(sol) << "." << std::endl;
            }
        }
        return 0;
    }
    catch (std::exception& e) {
        std::cout << "Exception: " << e.what() << std::endl;
        return -1;
    }
}