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 {
        int i;
        XPRBvar [] x,y;              // x-/y-coordinates of endpoints of chainlinks
        XPRBexpr qe;
        XPRBctr cobj, c;

        try (XPRBprob prob = new XPRBprob("catenary")) { // Initialize BCL and create a new problem

            /**** 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() );

        }
    }
}