Catenary - Solving a QCQP

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

using namespace xpress;
using namespace xpress::objects;

/*
 * QCQP problem (linear objective, convex quadratic constraints) Based on AMPL
 * model catenary.mod (Source:
 * http://www.orfe.princeton.edu/~rvdb/ampl/nlmodels/) This model finds the
 * shape of a hanging chain by minimizing its potential energy.
 */

const int N = 100;              // Number of chainlinks
const int L = 1;                // Difference in x-coordinates of endlinks
const double H = 2.0 * L / N;   // Length of each link

int main() {
    try {
        XpressProblem prob;

        /* VARIABLES */

        // Continuous variables to denote x and y positions of each chainlink
        std::vector<Variable> x = prob.addVariables(N + 1)
                .withName("x(%d)")
                .withLB(XPRS_MINUSINFINITY) // Overwrite default lowerbound from 0 to minus infinity
                .toArray();

        std::vector<Variable> y = prob.addVariables(N + 1)
                .withName("y(%d)")
                .withLB(XPRS_MINUSINFINITY) // Overwrite default lowerbound from 0 to minus infinity
                .toArray();

        /* CONSTRAINTS */

        // Position of left endpoint / anchor
        x[0].fix(0);
        y[0].fix(0);
        // Position of right endpoint / anchor
        x[N].fix(L);
        y[N].fix(0);

        // Positions of chainlinks: forall(j in 1..N) (x(j)-x(j-1))^2 + (y(j)-y(j-1))^2 <= H^2
        for (int i=1; i<=N; i++) {
            // TODO: Change `SumExpression` to `QuadExpression` (using QuadExpression::create().addTerm())
            SumExpression link_i = utils::sum(
                LinExpression::create().addTerm(x[i]).addTerm(x[i - 1], -1).square(),
                LinExpression::create().addTerm(y[i]).addTerm(y[i - 1], -1).square()
            );
            prob.addConstraint(link_i <= H * H).setName(xpress::format("link_%d", i));
        }

        /* OBJECTIVE */

        // Minimise the potential energy: sum(j in 1..N) (y(j-1) + y(j)) / 2
        LinTermList obj;
        for (int i=1; i<=N; i++) {
            obj.addTerm(0.5, y[i-1]).addTerm(0.5, y[i]); // Add ( y[j-1] + y[j] ) / 2
        }
        obj.compress(); // Add coefficients of the same variable together: 0.5*y[1] + 0.5*y[1] ---> 1*y[1]
        prob.setObjective(obj, ObjSense::Minimize);

        /* SOLVE & PRINT */

        prob.writeProb("Catenary.lp");
        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());
        }

        std::vector<double> sol = prob.getSolution();

        std::cout << "*** Solution ***" << std::endl;
        std::cout << "Objective value: " << prob.attributes.getObjVal() << std::endl;
        for (int i=0; i<=N; i++) {
            std::cout << i << ": " << x[i].getValue(sol) << ", " << y[i].getValue(sol) << std::endl;
        }
        return 0;
    }
    catch (std::exception& e) {
        std::cout << "Exception: " << e.what() << std::endl;
        return -1;
    }
}