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

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

#include <stdio.h>
#include "xprb.h"
#include "xprs.h"

#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;
 XPRBprob prob;
 XPRBvar x[N+1], y[N+1];       /* x-/y-coordinates of endpoints of chainlinks */
 XPRBctr cobj, c;

 prob=XPRBnewprob("catenary");           /* Initialize a new problem in BCL */

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

/****OBJECTIVE****/
/* Minimise the potential energy:  sum(j in 1..N) (y(j-1)+y(j))/2  */
 cobj = XPRBnewctr(prob, "Obj", XPRB_N);
 for(i=1;i<=N;i++)
 {
  XPRBaddterm(cobj, y[i-1], 0.5);
  XPRBaddterm(cobj, y[i], 0.5);
 }
 XPRBsetobj(prob, cobj);                 /* Set the 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++)
 {
  c = XPRBnewctr(prob, XPRBnewname("Link_%d",i), XPRB_L);
  XPRBaddqterm(c, x[i], x[i], 1);
  XPRBaddqterm(c, x[i], x[i-1], -2);
  XPRBaddqterm(c, x[i-1], x[i-1], 1);
  XPRBaddqterm(c, y[i], y[i], 1);
  XPRBaddqterm(c, y[i], y[i-1], -2);
  XPRBaddqterm(c, y[i-1], y[i-1], 1);
  XPRBaddterm(c, NULL, H*H);
 }

/****SOLVING + OUTPUT****/
 XPRBsetsense(prob, XPRB_MINIM);         /* Choose the sense of optimization */

/* Problem printing and matrix output: */
/*
 XPRBprintprob(prob);
 XPRBexportprob(prob, XPRB_MPS, "catenary");
 XPRBexportprob(prob, XPRB_LP, "catenary");
*/

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

 XPRBlpoptimize(prob,"");               /* Solve the problem */

 printf("Solution: %g\n", XPRBgetobjval(prob));
 for(i=0;i<=N;i++)
  printf(" %d: %g, %g\n", i, XPRBgetsol(x[i]), XPRBgetsol(y[i]));

 return 0;
}