Catenary - Solving a QCQP

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

  file xbcatena.cs
  ````````````````
  QCQP test problem
  Based on AMPL model catenary.mod
  (Source: http://www.orfe.princeton.edu/~rvdb/ampl/nlmodels/ )
 
  This model finds the shape of a hanging chain.
   
  (c) 2008 Fair Isaac Corporation
      authors: S.Heipcke, D.Brett, June 2008
****************************************************************/

using System;
using System.Text;
using System.IO;
using Optimizer;
using BCL;


namespace Examples
{

    public class TestBCatena
	{
        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

        public static void Main()
        {
            XPRB.init();
            int i;
            XPRBvar[] x = new XPRBvar[N + 1];
            XPRBvar[] y = new XPRBvar[N + 1];       // x-/y-coordinates of endpoints of chainlinks
            XPRBexpr qe;
            XPRBctr cobj;
            XPRBprob prob = new XPRBprob("catenary");  // Initialize a new problem in BCL 
            XPRSprob xprob;

            prob.setDictionarySize(BCLconstant.XPRB_DICT_NAMES,0);

            /**** VARIABLES ****/
            for(i=0;i<=N;i++) 
                x[i] = prob.newVar("x("+i+")", BCLconstant.XPRB_PL, -BCLconstant.XPRB_INFINITY, 
                         BCLconstant.XPRB_INFINITY);
            for(i=0;i<=N;i++) 
                y[i] = prob.newVar("y("+i+")", BCLconstant.XPRB_PL, -BCLconstant.XPRB_INFINITY, 
                         BCLconstant.XPRB_INFINITY);
            // Left anchor
            x[0].fix(0);  y[0].fix(0);
            // Right anchor
            x[N].fix(L);  y[N].fix(0);


            /****OBJECTIVE****/
            /*  sum(j in 1..N) (y(j-1)+y(j))/2  */
            qe = new XPRBexpr();
            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 ****/
            /*  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("Link_"+i, (x[i]-x[i-1]).sqr() + (y[i]-y[i-1]).sqr() <= H*H);

            /****SOLVING + OUTPUT****/
            prob.setSense(BCLconstant.XPRB_MINIM);   // Choose the sense of optimization

            /* Problem printing and matrix output: */
            /*
            prob.print(); 
            prob.exportProb(BCLconstant.XPRB_MPS, "caternary");
            prob.exportProb(BCLconstant.XPRB_LP, "caternary");
            prob.loadMat();
            */

            //Disable the convexivity check
            xprob = prob.getXPRSprob();
            xprob.IfCheckConvexity = 0;

            //Solve the problem
            prob.lpOptimize();                       

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

            return;
        }
    }
}