/********************************************************
 * 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-2024 Fair Isaac Corporation
 * author: S.Heipcke, June 2008, rev. Mar. 2011
 ********************************************************/

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

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 {
    try (XPRBprob prob = new XPRBprob("catenary"); // Initialize BCL and create a new problem
        XPRBexprContext context =
            new XPRBexprContext() /* Release XPRBexpr instances at end of block. */) {
      int i;
      XPRBvar[] x, y; // x-/y-coordinates of endpoints of chainlinks
      XPRBexpr qe;
      XPRBctr cobj, c;

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