Catenary - Solving a QCQP

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

  file xbcatena.java
  ``````````````````
  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
********************************************************/

import java.io.*;
import com.dashoptimization.*;

public class xbcatena
{
 static final int N = 100;	    // Number of chainlinks
 static final int L = 1;	    // Difference in x-coordinates of endlinks

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


 public static void main(String[] args) throws IOException
 {
  XPRB bcl;
  int i;
  XPRBvar [] x,y;              // x-/y-coordinates of endpoints of chainlinks
  XPRBexpr qe;
  XPRBctr cobj, c;
  XPRBprob prob;

  bcl = new XPRB();                     // Initialize BCL
  prob = bcl.newProb("catenary");       // Create a new problem in BCL

/**** VARIABLES ****/
  x = new XPRBvar[N+1];
  for(i=0;i<=N;i++) 
   x[i] = prob.newVar("x(" + i + ")", XPRB.PL, -XPRB.INFINITY,
                     XPRB.INFINITY);
  y = new XPRBvar[N+1];
  for(i=0;i<=N;i++) 
   y[i] = prob.newVar("y(" + 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=new XPRBexpr();
  for(i=1;i<=N;i++) qe .add( y[i-1].add(y[i]) );
  cobj = prob.newCtr("Obj", qe.mul(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("Link_"+i, (x[i].add(x[i-1].mul(-1))).sqr() 
                    .add( (y[i].add(y[i-1].mul(-1))).sqr() ).lEql(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");
*/

/* Start values:
  for(i=0;i<=N;i++) x[j].setInitval(j*L/N);
  for(i=0;i<=N;i++) y[j].setInitval(0);
*/

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

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

 }
}