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

  file xbairport.c
  ````````````````
  QCQP problem by
      Rodrigo de Barros Nabholz & Maria Aparecida Diniz Ehrhardt
      November 1994, DMA - IMECC- UNICAMP.
  Based on AMPL model airport.mod by Hande Y. Benson
  (Source: http://www.orfe.princeton.edu/~rvdb/ampl/nlmodels/ )

  (c) 2008-2024 Fair Isaac Corporation
      author: S.Heipcke, June 2008, rev. Mar. 2011
********************************************************/

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

#define N 42

double CX[] = {-6.3, -7.8, -9, -7.2, -5.7, -1.9, -3.5, -0.5, 1.4, 4,
               2.1, 5.5, 5.7, 5.7, 3.8, 5.3, 4.7, 3.3, 0, -1, -0.4, 4.2,
               3.2, 1.7, 3.3, 2, 0.7, 0.1, -0.1, -3.5, -4, -2.7, -0.5, -2.9,
               -1.2, -0.4, -0.1, -1, -1.7, -2.1, -1.8, 0};
double CY[] = {8, 5.1, 2, 2.6, 5.5, 7.1, 5.9, 6.6, 6.1, 5.6, 4.9, 4.7,
               4.3, 3.6, 4.1, 3, 2.4, 3, 4.7, 3.4, 2.3, 1.5, 0.5, -1.7, -2,
               -3.1, -3.5, -2.4, -1.3, 0, -1.7, -2.1, -0.4, -2.9, -3.4, -4.3,
               -5.2, -6.5, -7.5, -6.4, -5.1, 0};
double R[] = {0.09, 0.3, 0.09, 0.45, 0.5, 0.04, 0.1, 0.02, 0.02, 0.07, 0.4,
              0.045, 0.05, 0.056, 0.36, 0.08, 0.07, 0.36, 0.67, 0.38, 0.37,
              0.05, 0.4, 0.66, 0.05, 0.07, 0.08, 0.3, 0.31, 0.49, 0.09,
              0.46, 0.12, 0.07, 0.07, 0.09, 0.05, 0.13, 0.16, 0.46, 0.25, 0.1};

int main(int argc, char **argv)
{
 int i,j;
 XPRBprob prob;
 XPRBvar x[N],y[N];                     /* x-/y-coordinates to determine */
 XPRBctr cobj, c;

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

/**** VARIABLES ****/
 for(i=0;i<N;i++)
  x[i] = XPRBnewvar(prob, XPRB_PL, XPRBnewname("x(%d)",i), -10, 10);
 for(i=0;i<N;i++)
  y[i] = XPRBnewvar(prob, XPRB_PL, XPRBnewname("y(%d)",i), -10, 10);

/****OBJECTIVE****/
/* Minimize the total distance between all points */
/*  sum(i in 1..N-1,j in i+1..N) ((x(i)-x(j))^2+(y(i)-y(j))^2)  */
 cobj = XPRBnewctr(prob, "TotDist", XPRB_N);
  for(i=0;i<N-1;i++)
   for(j=i+1;j<N;j++)
   {
    XPRBaddqterm(cobj, x[i], x[i], 1);
    XPRBaddqterm(cobj, x[i], x[j], -2);
    XPRBaddqterm(cobj, x[j], x[j], 1);
    XPRBaddqterm(cobj, y[i], y[i], 1);
    XPRBaddqterm(cobj, y[i], y[j], -2);
    XPRBaddqterm(cobj, y[j], y[j], 1);
   }
 XPRBsetobj(prob, cobj);                /* Set the objective function */

/**** CONSTRAINTS ****/
/* All points within given distance of their target location */
/*  (x(i)-CX(i))^2+(y(i)-CY(i))^2 <= R(i)  */
 for(i=0;i<N;i++)
 {
  c = XPRBnewctr(prob, XPRBnewname("LimDist_%d",i), XPRB_L);
  XPRBaddqterm(c, x[i], x[i], 1);
  XPRBaddterm(c, x[i], -2*CX[i]);
  XPRBaddterm(c, NULL, -CX[i]*CX[i]);
  XPRBaddqterm(c, y[i], y[i], 1);
  XPRBaddterm(c, y[i], -2*CY[i]);
  XPRBaddterm(c, NULL, -CY[i]*CY[i]);
  XPRBaddterm(c, NULL, R[i]);
 }

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

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

 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;
}