Catenary - Solving a QCQP

/********************************************************
  BCL Example Problems
  ====================

  file xbcatena.cxx
  `````````````````
  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.
   
  (c) 2008 Fair Isaac Corporation
      author: S.Heipcke, June 2008, rev. Mar. 2011
********************************************************/

#include <iostream>
#include "xprb_cpp.h"
#include "xprs.h"

using namespace std;
using namespace ::dashoptimization;

#define N 100	               // Number of chainlinks
#define L 1	               // Difference in x-coordinates of endlinks

double H = 2.0*L/N;	               // Length of each link


int main(int argc, char **argv)
{
 int i;
 XPRBvar x[N+1],y[N+1];       // x-/y-coordinates of endpoints of chainlinks
 XPRBexpr qe;
 XPRBctr cobj, c;
 XPRBprob prob("catenary");              // Initialize a new problem in BCL

 prob.setDictionarySize(XPRB_DICT_NAMES,0);

/**** VARIABLES ****/
 for(i=0;i<=N;i++) 
  x[i] = prob.newVar(XPRBnewname("x(%d)",i), XPRB_PL, -XPRB_INFINITY, 
                     XPRB_INFINITY);
 for(i=0;i<=N;i++) 
  y[i] = prob.newVar(XPRBnewname("y(%d)",i), XPRB_PL, -XPRB_INFINITY, 
                     XPRB_INFINITY);
// Bounds: positions of endpoints
// Left anchor
 x[0].fix(0);  y[0].fix(0);
// Right anchor
 x[N].fix(L);  y[N].fix(0);


/****OBJECTIVE****/
/* Minimise the potential energy:  sum(j in 1..N) (y(j-1)+y(j))/2  */
 qe=0;
 for(i=1;i<=N;i++) qe+= y[i-1]+y[i];
 cobj = prob.newCtr("Obj", qe*0.5 );
 prob.setObj(cobj);                     // Set objective function 

/**** CONSTRAINTS ****/
/* Positions of chainlinks:
    forall(j in 1..N) (x(j)-x(j-1))^2+(y(j)-y(j-1))^2 <= H^2  */
 for(i=1;i<=N;i++)
  prob.newCtr(XPRBnewname("Link_%d",i), sqr(x[i]-x[i-1]) + sqr(y[i]-y[i-1]) <= H*H);
   
/****SOLVING + OUTPUT****/
 prob.setSense(XPRB_MINIM);             // Choose the sense of optimization
 
/* Problem printing and matrix output: */
/*
 prob.print(); 
 prob.exportProb(XPRB_MPS, "catenary");
 prob.exportProb(XPRB_LP, "catenary");
 prob.loadMat();
*/

/* Start values:
 for(i=0;i<=N;i++) x[j].setInitval(j*L/N)
 for(i=0;i<=N;i++) y[j].setInitval(0)
*/
                                       // Disable convexity check
 XPRSsetintcontrol(prob.getXPRSprob(), XPRS_IFCHECKCONVEXITY, 0);

 prob.lpOptimize("");                  // Solve the problem

 cout << "Solution: " << prob.getObjVal() << endl;
 for(i=0;i<=N;i++)
  cout << i << ": " << x[i].getSol() << ", " << y[i].getSol() << endl;

 return 0;
}